text:
struct Dim{N}
    Dim{1}() = new()
    Dim{2}() = new()
    Dim{3}() = new()
end

const DIM1 = Dim{1}()
const DIM2 = Dim{2}()
const DIM3 = Dim{3}()

otherdims(::Dim{1}) = [2, 3]
otherdims(::Dim{2}) = [3, 1]
otherdims(::Dim{3}) = [1, 2]

abstract type ReferenceCoords end
struct LowerEdges <: ReferenceCoords end
struct Centers <: ReferenceCoords end

abstract type SliceShift end
struct Up <: SliceShift end
struct Down <: SliceShift end
struct UpDown <: SliceShift end

abstract type IndexType end
struct Discrete <: IndexType end
struct Continuous <: IndexType end

struct Grid{N}
    shape::SVectorU{N}
    extents::SVectorF{N}
    origins::SVectorF{N}
    cell_extents::SVectorF{N}
    cell_extents⁻¹::SVectorF{N}
    cell_extents⁻²::SVectorF{N}
    centers::SVector{N,Vector{Float}}
    lower_edges::SVector{N,Vector{Float}}
    interior::SVector{N,UnitRange}
    interior_up::SVector{N,SVector{N,UnitRange}}
    interior_dn::SVector{N,SVector{N,UnitRange}}

    function Grid(
        shape::Union{NTuple{N},SVectorU{N}},
        extents::Union{NTuple{N},SVectorF{N}},
        origins::Union{NTuple{N},SVectorF{N}},
    ) where {N}
        @check all(shape .> 2)
        @check all(extents .> 0.0)
        cell_extents = extents ./ shape
        cell_extents⁻¹ = 1 ./ cell_extents
        cell_extents⁻² = cell_extents⁻¹ .* cell_extents⁻¹
        centers = ntuple(
            n -> range(
                origins[n] + 0.5 * cell_extents[n],
                origins[n] + extents[n] - 0.5 * cell_extents[n],
                length = shape[n],
            ),
            N,
        )
        lower_edges = ntuple(
            n -> range(
                origins[n],
                origins[n] + extents[n] - cell_extents[n],
                length = shape[n],
            ),
            N,
        )
        interior = ntuple(n -> 2:shape[n]-1, N)
        interior_up = ntuple(
            n -> SVector(ntuple(
                m -> ((m == n) ? (3:shape[m]) : (2:shape[m]-1)),
                N,
            )),
            N,
        )
        interior_dn = ntuple(
            n -> SVector(ntuple(
                m -> ((m == n) ? (1:shape[m]-2) : (2:shape[m]-1)),
                N,
            )),
            N,
        )
        new{N}(
            shape,
            extents,
            origins,
            cell_extents,
            cell_extents⁻¹,
            cell_extents⁻²,
            centers,
            lower_edges,
            interior,
            interior_up,
            interior_dn,
        )
    end

    function Grid(
        shape::Union{NTuple{N},SVectorU{N}},
        extents::Union{NTuple{N},SVectorF{N}},
    ) where {N}
        Grid(shape, extents, ntuple(_ -> 0.0, N))
    end
end

Base.size(grid::Grid) = grid.shape.data

extents(grid::Grid) = grid.extents

origins(grid::Grid) = grid.origins

interior(grid::Grid) = grid.interior

interior_dn(grid::Grid) = grid.interior_dn

interior_up(grid::Grid) = grid.interior_up

lower_edges(grid::Grid) = grid.lower_edges

centers(grid::Grid) = grid.centers

cell_extents(grid::Grid) = grid.cell_extents

cell_extents⁻¹(grid::Grid) = grid.cell_extents⁻¹

cell_extents⁻²(grid::Grid) = grid.cell_extents⁻²

cell_volume(grid::Grid) = prod(cell_extents(grid))

cell_volume⁻¹(grid::Grid) = 1 / cell_volume(grid)

coordinates(grid::Grid, ::Type{LowerEdges}) = lower_edges(grid)

coordinates(grid::Grid, ::Type{Centers}) = centers(grid)

function padded(grid::Grid; padding = 1)
    new_shape = size(grid) .+ 2 .* padding
    new_extents = extents(grid) .+ 2 .* padding .* cell_extents(grid)
    new_origins = origins(grid) .- padding .* cell_extents(grid)
    Grid(new_shape, new_extents, new_origins)
end

function find_indices_at(
    grid::Grid{N},
    position,
    reference::Type{<:ReferenceCoords};
    index_type::Type{<:IndexType} = Discrete,
) where {N}
    [
        find_index_at(
            grid,
            position[n],
            Dim{n}(),
            reference,
            index_type = index_type,
        ) for n = 1:N
    ]
end

function find_indices_at(
    grid::Grid{N},
    position,
    references;
    index_type::Type{<:IndexType} = Discrete,
) where {N}
    [
        find_index_at(
            grid,
            position[n],
            Dim{n}(),
            references[n],
            index_type = index_type,
        ) for n = 1:N
    ]
end

function find_index_at(
    grid::Grid,
    coord,
    ::Dim{D},
    ::Type{LowerEdges};
    index_type::Type{<:IndexType} = Discrete,
) where {D}
    compute_index_unform(
        size(grid)[D],
        lower_edges(grid)[D][1],
        extents(grid)[D],
        coord,
        index_type,
    )
end

function find_index_at(
    grid::Grid,
    coord,
    ::Dim{D},
    ::Type{Centers};
    index_type::Type{<:IndexType} = Discrete,
) where {D}
    compute_index_unform(
        size(grid)[D],
        centers(grid)[D][1],
        extents(grid)[D],
        coord,
        index_type,
    )
end

function compute_index_unform(
    length,
    lower_edge,
    extent,
    coord,
    ::Type{Discrete},
)
    floor(
        Int,
        compute_index_unform(length, lower_edge, extent, coord, Continuous),
    )
end

function compute_index_unform(
    length,
    lower_edge,
    extent,
    coord,
    ::Type{Continuous},
)
    1 + (coord - lower_edge) * (length - 1) / extent
end

function shifteddiff(
    field::Array{<:Any,N},
    grid::Grid{N},
    ::Type{Up},
    ::Dim{D},
) where {N,D}
    @dbgasserteq(size(field), size(grid))
    field[interior_up(grid)[D]...] .- field[interior(grid)...]
end

function shifteddiff(
    field::Array{<:Any,N},
    grid::Grid{N},
    ::Type{Down},
    ::Dim{D},
) where {N,D}
    @dbgasserteq(size(field), size(grid))
    field[interior(grid)...] .- field[interior_dn(grid)[D]...]
end

function shiftedsum(
    field::Array{<:Any,N},
    grid::Grid{N},
    ::Type{UpDown},
    ::Dim{D},
) where {N,D}
    @dbgasserteq(size(field), size(grid))
    field[interior_up(grid)[D]...] .+ field[interior_dn(grid)[D]...]
end

max_stars_repo_name: lars-frogner/PIC

max_issues_repo_name: lars-frogner/PIC

max_forks_repo_name: lars-frogner/PIC


text:
# This file is a part of BATBase.jl, licensed under the MIT License (MIT).

using BATBase
using Test

@testset "hello_world" begin
    @test BATBase.hello_world() == 42
end

ext: jl


text:
using PyPlot

function plot_solution(n::Int64,sol::Array{Float64,1},desiredstate::Function,alpha::Float64)
	# Simple scatter-plot of the nodes
	# Input: solution sol=(y,p,lambda)
	#        Desired state: y_d
	#	 L2-regularization alpha

	n_block::Int64 = Int64(length(sol)/3)

	dof_map = load("data"*string(n)*".jld2", "dof_map")
	xval::Array{Float64,1} = zeros(n_block)
	yval::Array{Float64,1} = zeros(n_block)

	stateval::Array{Float64,1} = zeros(n_block)
	controlval::Array{Float64,1} = zeros(n_block)
	adjointval::Array{Float64,1} = zeros(n_block)
	lagrangeval::Array{Float64,1} = zeros(n_block)
	y_d::Array{Float64,1} = zeros(n_block)

	# Separating all variables into vectors
	for i=1:n_block
		xval[i] = dof_map[i,:][1]
		yval[i] = dof_map[i,:][2]
		y_d[i] = desiredstate(xval[i],yval[i])
		stateval[i] = sol[i]
		adjointval[i] = sol[i+n_block]
		lagrangeval[i] = sol[i+2*n_block]
		controlval[i] = (1. /alpha)*(adjointval[i] - lagrangeval[i] )
	end

	# Plot state
	figure("1")
	scatter3D(xval, yval, stateval, s=0.1, c="k", label="State")
	scatter3D(xval, yval, y_d, s=0.1, c="b", alpha=0.1, label="Desired state")
	title("State")
	xlabel("x_1")
	ylabel("x_2")
	zlabel("y")
	legend(loc="best")

	# Plot Control
	figure("2")
	scatter3D(xval, yval, controlval, s=0.1, c="k")
	title("Control")
	xlabel("x_1")
	ylabel("x_2")
	zlabel("u")
end

max_stars_repo_name: dravinsi/demoPDEOPT

max_issues_repo_name: dravinsi/demoPDEOPT

max_forks_repo_name: dravinsi/demoPDEOPT


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains methods to generate graphs

################################################# IMPORT/EXPORT ############################################################
export
# Types
emptygraph, randgraph, completegraph, propgraph, samplegraph

################################################# CONSTRUCTION #############################################################

# NV
""" Build an empty graph with nv vertices """
function Graph(nv::Int)
   Graph(nv, 0, spzeros(Int, nv, nv), DataFrame(), DataFrame(), LabelMap(nv))
end

# NV & NE
""" Build a graph with nv vertices and (approximately)ne edges """
function Graph(nv::Int, ne::Int)
   sv = randindxs(nv, ne)
   Graph(nv, nnz(sv), sv, DataFrame(), DataFrame(), IdentityLM(nv))
end

# NV & LS
""" Build an empty graph with nv labelled vertices """
function Graph(nv::Int, ls::AbstractVector)
   if length(ls) == nv
      Graph(nv, 0, spzeros(Int, nv, nv), DataFrame(), DataFrame(), LabelMap(ls))
   else
      error("Trying to assign $(length(ls)) labels to $nv labels")
   end
end

# NV & LS & NE
""" Build a graph with nv labelled vertices and (approximately)ne edges """
function Graph(nv::Int, ls::Vector, ne::Int)
   sv = randindxs(nv, ne)
   Graph(nv, nnz(sv), sv, DataFrame(), DataFrame(), LabelMap(ls))
end

# SparseMatrixCSC
""" Build a graph using the input adjacency matrix """
function Graph(x::SparseMatrixCSC)
   nv = size(x, 1)
   ne = nnz(x)

   sv = copy(x)
   reorder!(sv)

   Graph(nv, ne, sv, DataFrame(), DataFrame(), LabelMap(nv))
end

################################################# EMTPY GRAPH ##############################################################

""" Build an empty graph with nv vertices """
emptygraph(nv::Int) = Graph(nv)

################################################# RANDGRAPH ################################################################

""" Build a graph with nv vertices and (approximately)ne edges """
randgraph(nv::Int, ne::Int) = Graph(nv, ne)

""" Build a graph with nv vertices and a random number of edges """
randgraph(nv::Int) = Graph(nv, rand(1 : (nv * (nv-1))))

###
# TODO: MORE RANDOM GENERATORS
###

################################################# COMPLETE GRAPH ###########################################################

""" Build a complete graph with nv vertices """
function completegraph(nv::Int)
   Graph(completeindxs(nv))
end

################################################# PROPGRAPH ################################################################

""" Returns a small completegraph with properties(for doc examples) """
function propgraph(nv::Int, vprops::Vector{Symbol}, eprops::Vector{Symbol})
   g = completegraph(nv)

   for prop in vprops
      setvprop!(g, :, rand(nv), prop)
   end

   for prop in eprops
      seteprop!(g, :, rand(ne(g)), prop)
   end

   return g
end
################################################# SAMPLE GRAPHS ############################################################

_graph_text = IOBuffer("10\t28\tString\nage\toccupation\nInt\tString\nrelationship\nString\nAbel\t32\tArchitect\nBharath\t35\tBanker\nCamila\t21\tChef\nDalia\t59\tDentist\nEduardo\t24\tEconomist\nFabiola\t37\tFlorist\nGaurav\t46\tGeologist\nHector\t31\tHistorian\nIgnacio\t28\tInterpreter\nJanardan\t59\tJudge\nAbel	Bharath	father\nBharath	Abel	son\nAbel	Camila	friend\nCamila	Abel	friend\nBharath	Camila	familyfriend\nCamila	Bharath	familyfriend\nCamila	Dalia	wife\nDalia	Camila	husband\nDalia	Eduardo	friend\nEduardo	Dalia	friend\nDalia	Fabiola	neighbor\nFabiola	Dalia	neighbor\nDalia	Gaurav	teacher\nGaurav	Dalia	student\nDalia	Hector	father\nHector	Dalia	son\nDalia	Ignacio	daughter\nIgnacio	Dalia	mother\nDalia	Janardan	aunt\nJanardan	Dalia	nephew\nEduardo	Fabiola	friend\nFabiola	Eduardo\tfriend\nFabiola	Gaurav	friend\nGaurav	Fabiola	friend\nHector	Ignacio	granddaughter\nIgnacio	Hector	grandmother\nIgnacio	Janardan	sister\nJanardan	Ignacio	sister")

function samplegraph()
   loadgraph(_graph_text)
end


text:
using Gadfly, RDatasets

set_default_plot_size(6inch, 3inch)

mpg = dataset("ggplot2","mpg")
plot(mpg, x=:Class, y=:Hwy, Geom.boxplot, color=:Class)


text:
export TanhLayer

@neuron type TanhNeuron
end

@neuron forward(neuron::TanhNeuron) do
    neuron.value = tanh(neuron.inputs[1])
end

∇tanh(x) = x * (1.0f0 - x)

@neuron backward(neuron::TanhNeuron) do
    neuron.∇inputs[1] = ∇tanh(neuron.value) * neuron.∇
end

function TanhLayer(name::Symbol, net::Net, input_ensemble::AbstractEnsemble; copy=false)
    neurons = Array(TanhNeuron, size(input_ensemble)...)
    for i in 1:length(neurons)
        neurons[i] = TanhNeuron()
    end
    if copy
        ens = Ensemble(net, name, neurons)
        mapping = one_to_one(ndims(input_ensemble))
        add_connections(net, input_ensemble, ens, mapping)
        ens
    else
        ActivationEnsemble(net, name, neurons, input_ensemble)
    end
end

function tanh(net::Net, ens::AbstractEnsemble; copy=false)
    TanhLayer(gensym("ensemble"), net, ens; copy=copy)
end

export SigmoidLayer

sigmoid(x) = 1 / (1 + exp(-x))
∇sigmoid(x) = 1.0f0 - x * x

@neuron type SigmoidNeuron
end

@neuron forward(neuron::SigmoidNeuron) do
    neuron.value = sigmoid(neuron.inputs[1])
end

@neuron backward(neuron::SigmoidNeuron) do
    neuron.∇inputs[1] = ∇sigmoid(neuron.value) * neuron.∇
end

function SigmoidLayer(name::Symbol, net::Net, input_ensemble::AbstractEnsemble)
    neurons = Array(SigmoidNeuron, size(input_ensemble)...)
    for i in 1:length(neurons)
        neurons[i] = SigmoidNeuron()
    end
    ActivationEnsemble(net, name, neurons, input_ensemble)
end

function σ(net, ens::AbstractEnsemble)
    SigmoidLayer(gensym("ensemble"), net, ens)
end


text:
using Registrator
Registrator.RegServer.main()


text:
function FiniteDifferences.to_vec(l::LoLa)
    n = length(l.w_ds[i])
    return (
        [vec(l.w_E); vec.(l.w_ds)...],
        function(v)
            LoLa(
                reshape(v[1:length(l.w_E)], axes(l.w_E)),
                ntuple(length(l.w_ds)) do i
                    reshape(
                        v[length(l.w_E) + mapreduce(i->length(l.w_ds[i]), +, 1:i-1; init=0) .+ (1:n)],
                        axes(l.w_ds[i]),
                    )
                end,
                l.w_d_reducers,
            )
        end,
    )
end


text:
function wdot(wx::AbstractArray{T}, x::AbstractArray{T},
                 wy::AbstractArray{T}, y::AbstractArray{T}) where T
    out = zero(T)
    @inbounds @simd for i in 1:length(x)
        out += conj(wx[i]*x[i]) * wy[i]*y[i]
    end
    return out
end

wnorm(w, x) = sqrt(real(wdot(w, x, w, x)))
assess_convergence(f, ftol) = assess_convergence(NaN, NaN, f, NaN, ftol)
function assess_convergence(x,
                            x_previous,
                            f,
                            xtol,
                            ftol)
    x_converged, f_converged = false, false

    if !any(isnan, x_previous) && chebyshev(x, x_previous) <= xtol
        x_converged = true
    end

    if maximum(abs, f) <= ftol
        f_converged = true
    end

    converged = x_converged || f_converged

    return x_converged, f_converged, converged
end

function check_isfinite(x::AbstractArray)
    if any(!isfinite, x)
        i = findall(!isfinite, x)
        throw(IsFiniteException(i))
    end
end

"""
    qrdelete!(Q, R, k)

Delete the left-most column of F = Q[:, 1:k] * R[1:k, 1:k] by updating Q and R.

Only Q[:, 1:(k-1)] and R[1:(k-1), 1:(k-1)] are valid on exit.
"""
function qrdelete!(Q::AbstractMatrix, R::AbstractMatrix, k::Int)
  n, m = size(Q)
  m == LinearAlgebra.checksquare(R) || throw(DimensionMismatch())
  1 ≤ k ≤ m || throw(ArgumentError())

  # apply Givens rotations
  for i in 2:k
      g = first(givens(R, i - 1, i, i))
      lmul!(g, R)
      rmul!(Q, g')
  end

  # move columns of R
  @inbounds for j in 1:(k-1)
    for i in 1:(k-1)
      R[i, j] = R[i, j + 1]
    end
  end

  Q, R
end

"""
    qradd!(Q, R, v, k)

Replace the right-most column of F = Q[:, 1:k] * R[1:k, 1:k] with v by updating Q and R.

This implementation modifies vector v as well. Only Q[:, 1:k] and R[1:k, 1:k] are valid on
exit.
"""
function qradd!(Q::AbstractMatrix, R::AbstractMatrix, v::AbstractVector, k::Int)
  n, m = size(Q)
  n == length(v) || throw(DimensionMismatch())
  m == LinearAlgebra.checksquare(R) || throw(DimensionMismatch())
  1 ≤ k ≤ m || throw(ArgumentError())

  @inbounds for i in 1:(k-1)
    q = view(Q, :, i)
    r = dot(q, v)

    R[i, k] = r
    axpy!(-r, q, v)
  end

  @inbounds begin
    d = norm(v)
    R[k, k] = d
    @. Q[:, k] = v / d
  end

  Q, R
end


text:
using Documenter
using GeoJSON

makedocs(
    sitename = "GeoJSON",
    modules = [GeoJSON],
)

deploydocs(
    repo = "github.com/JuliaGeo/GeoJSON.jl.git"
)


text:
### Computing fields and other physical quantities ###
using LinearAlgebra
import LinearAlgebra: transpose!
import CartesianGrids: convective_derivative!

"""
    vorticity(u,sys::NavierStokes,t)

Return the vorticity field associated with state vector `u` on the grid in `sys`,
at time `t`.
"""
@inline vorticity(w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY},t::Real) where {NX,NY} =
            vorticity!(sys.Wn,w,sys,t)

function vorticity!(w::Nodes{Dual,NX,NY},win::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY},t::Real) where {NX,NY}
  w .= 0.0
  _vorticity!(w,win,sys,t)
  _vorticity_sheet!(w,sys,t)
  return w
end

@inline _vorticity!(w::Nodes{Dual,NX,NY},win::Nodes{Dual,NX,NY},
                    sys::NavierStokes{NX,NY},t::Real) where {NX,NY} = w .+= win/cellsize(sys)

_vorticity_sheet!(w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY,0},t::Real) where {NX,NY} = w

_vorticity_sheet!(w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY,N,MT,FS,
                                              ExternalInternalFlow},t::Real) where {NX,NY,N,MT,FS} = w

function _vorticity_sheet!(w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY,N,MT,FS,SD},t::Real) where {NX,NY,N,MT,FS,SD}
  if isnothing(sys.Rn)
    regop = _regularization(sys)
    Rn = RegularizationMatrix(regop,sys.Sb,sys.Sn)
  else
    Rn = sys.Rn
  end
  surface_velocity_jump!(sys.Δus,sys,t)
  cross!(sys.Sb,sys.Δus,normals(sys))
  w .+= Rn*sys.Sb
  return w
end

function velocity!(u::Edges{Primal,NX,NY},w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY,N},t::Real) where {NX,NY,N}
    u .= 0.0
    _velocity_vorticity!(u,w,sys)
    _velocity_single_layer!(u,sys,t)
    _velocity_freestream!(u,sys,t)
    return u
end

function _velocity_vorticity!(u::Edges{Primal,NX,NY},w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY,N}) where {NX,NY,N}
    sys.Wn .= 0.0
    _unscaled_streamfunction_vorticity!(sys.Wn,w,sys)
    sys.Vf .= 0.0
    curl!(sys.Vf,sys.Wn)
    u .+= sys.Vf
    return u
end

@inline _velocity_single_layer!(u::Edges{Primal,NX,NY},sys::NavierStokes{NX,NY,N,MT,FS,
                                      ExternalInternalFlow},t::Real) where {NX,NY,N,MT,FS} = u

@inline _velocity_single_layer!(u::Edges{Primal,NX,NY},sys::NavierStokes{NX,NY,0},t::Real) where {NX,NY} = u

function _velocity_single_layer!(u::Edges{Primal,NX,NY},sys::NavierStokes{NX,NY,N,MT,FS,SD},t::Real) where {NX,NY,N,MT,FS,SD}
    surface_velocity_jump!(sys.Δus,sys,t)
    pointwise_dot!(sys.Sb,sys.Δus,normals(sys))
    sys.Sb .*= cellsize(sys)
    sys.slc(sys.Sc,sys.Sb)
    _unscaled_scalarpotential!(sys.Sc,sys.Sc,sys)
    sys.Vf .= 0.0
    grad!(sys.Vf,sys.Sc)
    u .+= sys.Vf
    return u
end

function _velocity_freestream!(u::Edges{Primal,NX,NY},sys::NavierStokes{NX,NY,N,MT,FS,SD},t::Real) where {NX,NY,N,MT,FS,SD}
  U∞, V∞ = freestream(t,sys)
  u.u .+= U∞
  u.v .+= V∞
  return u
end

"""
    velocity(u,sys::NavierStokes,t)

Return the velocity field associated with state vector `u` on the grid in `sys`,
at time `t`.
"""
velocity(w::Nodes{Dual,NX,NY},a...) where {NX,NY} = velocity!(Edges(Primal,(NX,NY)),w,a...)

function streamfunction!(ψ::Nodes{Dual,NX,NY},w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY},t::Real) where {NX,NY}
  ψ .= 0.0
  _unscaled_streamfunction_vorticity!(ψ,w,sys)
  ψ .*= cellsize(sys)
  _streamfunction_freestream!(ψ,sys,t)
  return ψ
end

function _unscaled_streamfunction_vorticity!(ψ::Nodes{Dual,NX,NY},w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY}) where {NX,NY}
  ldiv!(sys.Sn,sys.L,w)
  ψ .-= sys.Sn
  return ψ
end

function _streamfunction_freestream!(ψ::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY},t::Real) where {NX,NY}
  xg, yg = coordinates(ψ,sys.grid)
  U∞, V∞ = freestream(t,sys)
  ψ .+= U∞*yg' .- V∞*xg
  return ψ
end

"""
    streamfunction(u,sys::NavierStokes,t)

Return the streamfunction field associated with state vector `u` on the grid in `sys`,
at time `t`.
"""
streamfunction(w::Nodes{Dual,NX,NY},sys::NavierStokes{NX,NY},t::Real) where {NX,NY} = streamfunction!(Nodes(Dual,(NX,NY)),w,sys,t)

function scalarpotential!(ϕ::Nodes{Primal,NX,NY},dil::Nodes{Primal,NX,NY},sys::NavierStokes{NX,NY},t::Real) where {NX,NY}
  _unscaled_scalarpotential!(ϕ,dil,sys)
  # This assumes that dil is the physical dilation times the grid spacing
  ϕ .*= cellsize(sys)
  return ϕ
end

function _unscaled_scalarpotential!(ϕ::Nodes{Primal,NX,NY},dil::Nodes{Primal,NX,NY},sys::NavierStokes{NX,NY}) where {NX,NY}
  ldiv!(ϕ,sys.L,dil)
  return ϕ
end

"""
    scalarpotential(u,sys::NavierStokes,t)

Return the scalar potential field associated with state vector `u` on the grid in `sys`,
at time `t`.
"""
scalarpotential(dil::Nodes{Primal,NX,NY},sys::NavierStokes{NX,NY},t::Real) where {NX,NY} = scalarpotential!(Nodes(Primal,(NX,NY)),dil,sys,t)

function convective_derivative!(out::Edges{Primal,NX,NY},u::Edges{Primal,NX,NY},
                        sys::NavierStokes{NX,NY},t::Real) where {NX,NY}
  out .= u
  _unscaled_convective_derivative!(out,sys)
  Δx⁻¹ = 1/cellsize(sys)
  out .*= Δx⁻¹
  return out
end

"""
    convective_derivative(u,sys::NavierStokes,t)

Return the convective derivative associated with state vector `u` on the grid in `sys`,
at time `t`.
"""
convective_derivative(u::Edges{Primal,NX,NY},sys::NavierStokes{NX,NY},t::Real) where {NX,NY} =
      convective_derivative!(Edges(Primal,(NX,NY)),u,sys,t)

"""
    vorticity(integrator)

Return the vorticity field associated with `integrator` at its current state.
"""
function vorticity end

"""
    velocity(integrator)

Return the velocity field associated with `integrator` at its current state.
"""
function velocity end

"""
    streamfunction(integrator)

Return the streamfunction field associated with `integrator` at its current state.
"""
function streamfunction end

"""
    scalarpotential(integrator)

Return the scalar potential field associated with `integrator` at its current state.
"""
function scalarpotential end

"""
    convective_derivative(integrator)

Return the convective derivative associated with `integrator` at its current state.
"""
function convective_derivative end

for fcn in (:vorticity,:velocity,:streamfunction,:scalarpotential,:convective_derivative)
  @eval $fcn(s::ConstrainedSystems.ArrayPartition,sys::NavierStokes,t) = $fcn(state(s),sys,t)

  @eval $fcn(integ::ConstrainedSystems.OrdinaryDiffEq.ODEIntegrator) = $fcn(integ.u,integ.p,integ.t)
end

function pressure(u::ConstrainedSystems.ArrayPartition,sys::NavierStokes{NX,NY},t) where {NX,NY}
    w, τ = state(u), constraint(u)
    fill!(sys.Vn,0.0)
    _vel_ns_rhs_convectivederivative!(sys.Vn,w,sys,t)
    _vel_ns_rhs_double_layer!(sys.Vn,sys,t)
    fill!(sys.Vf,0.0)
    _vel_ns_op_constraint_force!(sys.Vf,τ,sys)
    sys.Vn .-= sys.Vf
    fill!(sys.Sc,0.0)
    divergence!(sys.Sc,sys.Vn)
    press = zero(sys.Sc)
    ldiv!(press,sys.L,sys.Sc)
    press .*= cellsize(sys)
    return press
end
pressure(integ::ConstrainedSystems.OrdinaryDiffEq.ODEIntegrator) = pressure(integ.u,integ.p,integ.t)

# Surface field quantities

for fcn in (:traction,:pressurejump)
  @eval $fcn(s::ConstrainedSystems.ArrayPartition,a...;kwargs...) = $fcn(constraint(s),a...;kwargs...)
  @eval $fcn(integ::ConstrainedSystems.OrdinaryDiffEq.ODEIntegrator) = $fcn(integ.u,integ.p,integ.t)
end

traction(τ,sys::NavierStokes{NX,NY,0,MT,FS,SD},t) where {NX,NY,MT,FS,SD} = τ

traction(τ::VectorData{N},sys::NavierStokes{NX,NY,N,MT,FS,ExternalInternalFlow},t) where {NX,NY,N,MT,FS} = τ

function traction(τ::VectorData{N},sys::NavierStokes{NX,NY,N,MT,FS,SD},t) where {NX,NY,N,MT,FS,SD}
    nrm = normals(sys.bodies)
    relative_surface_velocity!(sys.Δus,sys,t)
    pointwise_dot!(sys.Sb,nrm,sys.Δus)
    surface_velocity_jump!(sys.Δus,sys,t)
    product!(sys.Vb,sys.Δus,sys.Sb)
    return τ + sys.Vb
end

function pressurejump(τ::VectorData{N},sys::NavierStokes{NX,NY,N},t) where {NX,NY,N}

    #_hasfilter(sys) ? (τf = similar(τ); τf .= sys.Cf*force(τ,sys)) : τf = deepcopy(force(τ,sys))
    #τf = deepcopy(τ)
    #return nrm.u∘τf.u + nrm.v∘τf.v # This might need some modification

    press = zero(sys.Sb)
    nrm = normals(sys.bodies)
    pointwise_dot!(press,nrm,traction(τ,sys,t))
    return -press
end

## Total quantities

force(τ,sys::NavierStokes{NX,NY,0},t,bodyi::Int) where {NX,NY} = Vector{Float64}(), Vector{Float64}()

function force(τ::VectorData{N},sys::NavierStokes,t,bodyi::Int) where {N}
    product!(sys.τ,traction(τ,sys,t),areas(sys.bodies))
    fxvec = sum(sys.τ.u,sys.bodies,bodyi)
    fyvec = sum(sys.τ.v,sys.bodies,bodyi)
    return fxvec, fyvec
end

force(s::ConstrainedSystems.ArrayPartition,sys,t,bodyi) = force(constraint(s),sys,t,bodyi)

"""
    force(integ,bodyindex)

Given the state of the system in integrator `integ`, return the current force
on the body with index `bodyindex`.
"""
force(integ::ConstrainedSystems.OrdinaryDiffEq.ODEIntegrator,bodyi) = force(integ.u,integ.p,integ.t,bodyi)

"""
    force(sol,sys::NavierStokes,bodyindex) -> Tuple{Vector,Vector}

Given the solution history vector `sol` and the system `sys`, return the force
history on the body with index `bodyindex` as a tuple of vectors, one for
each component.
"""
function force(sol::ConstrainedSystems.OrdinaryDiffEq.ODESolution,sys::NavierStokes,bodyi::Int)
    fx = map((u,t) -> force(u,sys,t,bodyi)[1],sol.u,sol.t)
    fy = map((u,t) -> force(u,sys,t,bodyi)[2],sol.u,sol.t)
    fx, fy
end

function moment(τ::VectorData{N},bodies::BodyList,sys::NavierStokes,t,bodyi::Int;center=(0.0,0.0)) where {N}
    xc, yc = center
    x, y = collect(bodies)
    dX = VectorData(x .- xc,  y .- yc)
    product!(sys.τ,traction(τ,sys,t),areas(bodies))
    return sum(cross(dX,sys.τ),bodies,bodyi)
end

function moment(s::ConstrainedSystems.ArrayPartition,sys::NavierStokes{NX,NY,N,MovingPoints},t,bodyi;center=(0.0,0.0)) where {NX,NY,N}
    bodies = deepcopy(sys.bodies)
    tl! = RigidTransformList(aux_state(s))
    tl!(bodies)
    moment(constraint(s),bodies,sys,t,bodyi;center=center)
end

moment(s::ConstrainedSystems.ArrayPartition,sys::NavierStokes{NX,NY,N,StaticPoints},t,bodyi;center=(0.0,0.0)) where {NX,NY,N} =
        moment(constraint(s),sys.bodies,sys,t,bodyi;center=center)

"""
    moment(integ,bodyindex[,center=(0,0))

Given the state of the system in integrator `integ`, return the current moment
on the body with index `bodyindex`.  The center of the
moment can be passed as an optional tuple argument.
"""
moment(integ::ConstrainedSystems.OrdinaryDiffEq.ODEIntegrator,bodyi;center=(0.0,0.0)) =
      moment(integ.u,integ.p,integ.t,bodyi,center=center)

"""
    moment(sol,sys::NavierStokes,bodyindex[,center=(0,0)]) -> Vector

Given the solution history vector `sol` and the system `sys`, return the moment
history on the body with index `bodyindex` as a vector. The center of the
moment can be passed as an optional tuple argument.
"""
function moment(sol::ConstrainedSystems.OrdinaryDiffEq.ODESolution,sys::NavierStokes,bodyi::Int;center=(0.0,0.0))
    mom = map((u,t) -> moment(u,sys,t,bodyi,center=center),sol.u,sol.t)
end


text:
@testset "Basics" begin
    m = Mock()
    @test isempty(calls(m))
    @test_throws BoundsError last_call(m)
    @test ncalls(m) == 0
    @test !called(m)
    @test !called_once(m)
    @test !called_with(m)
    @test_throws BoundsError called_last_with(m)
    @test !called_once_with(m)

    m()
    @test length((calls(m))) == 1
    @test last_call(m) == Call()
    @test ncalls(m) == 1
    @test called(m)
    @test called_once(m)
    @test called_with(m)
    @test called_last_with(m)
    @test called_once_with(m)

    m(1)
    @test !called_once(m)
    @test called_with(m, 1)
    @test !called_last_with(m)
    @test called_last_with(m, 1)
    @test !called_once_with(m, 1)

    m(1; x=2)
    @test called_with(m, 1; x=2)
    @test last_call(m) == Call(1; x=2)
    @test called_last_with(m, 1; x=2)
end

@testset "has_calls" begin
    m = Mock()
    @test has_calls(m)
    @test !has_calls(m, Call())

    m()
    @test has_calls(m, Call())
    @test !has_calls(m, Call(1))
    @test !has_calls(m, Call(; x=1))
    @test !has_calls(m, Call(), Call())

    m(1)
    @test has_calls(m, Call(), Call(1))

    m(2)
    m(3)
    m(4)
    @test has_calls(m, Call(2), Call(3))
    @test !has_calls(m, Call(2), Call(4))

    m(; x=1)
    @test has_calls(m, Call(4), Call(; x=1))
end

@testset "reset!" begin
    m = Mock()
    foreach(m, 1:10)
    @test ncalls(m) == 10
    reset!(m)
    @test !called(m)
end

@testset "Effects" begin
    m = Mock()
    @test m() isa Mock
    @test m(1; x=1) isa Mock
    @test m() != m()

    m = Mock(1)
    @test m() == 1
    @test m(1; x=1) == 1

    m = Mock([1, 2, 3])
    @test m() == 1
    @test m() == 2
    @test m() == 3
    @test_throws ArgumentError m()

    m = Mock(iseven)
    @test !m(1)
    @test m(2)

    m = Mock(KeyError(1))
    @test_throws KeyError m()

    m = Mock([KeyError(1), ArgumentError("foo")])
    @test_throws KeyError m()
    @test_throws ArgumentError m()

    m = Mock([KeyError(1), [], 1, ArgumentError("foo")])
    @test_throws KeyError m()
    @test m() == []
    @test m() == 1
    @test_throws ArgumentError m()

    @eval begin
        struct Foo end
        (::Foo)(x) = x
    end
    m = Mock(Foo())
    @test m(1) == 1
    @test_throws MethodError m(1, 2)
end

@testset "Show methods" begin
    mime = MIME("text/plain")

    @test sprint(show, mime, Call(1, 2, 3)) == "Call(1, 2, 3)"
    @test sprint(show, mime, Call('a')) == "Call('a')"
    @test sprint(show, mime, Call(; x=1, y='a')) == "Call(; x=1, y='a')"

    m = Mock()
    @test sprint(show, mime, m) == "Mock(id=$(m.id))"

    p = Predicate(iseven)
    @test sprint(show, mime, p) == "Predicate(iseven)"
end

@testset "Equality" begin
    @test Call() == Call()
    @test Call(1) != Call(2)
    @test Call(1) != Call()
    @test Call(; x=1) == Call(; x=1)
    @test Call(; x=1) != Call(; x=2)

    @test Mock() != Mock()
    @test Mock(1) != Mock(1)
end

@testset "Predicates" begin
    p = Predicate(iseven)

    m = Mock()
    @test !has_call(m, Call(p))
    m(1)
    @test !has_call(m, Call(p))
    m(2)
    @test has_call(m, Call(p))

    m = Mock()
    @test !has_call(m, Call(; x=p))
    m(; x=1)
    @test !has_call(m, Call(; x=p))
    m(; x=2)
    @test has_call(m, Call(; x=p))

    m = Mock()
    m(1, 2, 3; a=1, b=2, c=3)
    @test !has_call(m, Call(1, 2, p; a=1, b=2, c=3))
    @test !has_call(m, Call(1, 2, 3; a=1, b=2, c=p))
    @test has_call(m, Call(1, p, 3; a=1, b=p, c=3))

    m = Mock()
    foreach(m, 1:10)
    @test !has_calls(m, Call(1), Call(Predicate(==(3))))
    @test has_calls(m, Call(Predicate(==(5))), Call(Predicate(iseven)))
end


text:
"""

	G ∩ H = ∅

"""

struct EmptyIntersection <: Constraint

	G::Variable
	H::Variable

end

import Base.show
function Base.show(io::IO, inter::EmptyIntersection)
	print(io, inter.G.name * " ∩ " * inter.H.name * " = ∅")
end

function EmptyIntersection(model, G, H)
	constraint = EmptyIntersection(G, H)
	push!(G.linkedConstraint, constraint)
	push!(H.linkedConstraint, constraint)

	push!(model.constraints, constraint)

	constraint
end

function filtrage!(inter::EmptyIntersection)

	stop = !(isempty(intersect(inter.H.lowerBound, inter.G.lowerBound)))

	changeVariable = Vector{Variable}(undef, 3)
	nbChange = 0

	changeG = Change(inter.G)
	changeH = Change(inter.H)

	if stop
		return [], (changeG, changeH), stop
	else

		#Set
		if !inter.G.isFixed

			rem = intersect(inter.G.upperBound, inter.H.lowerBound)

			setdiff!(inter.G.upperBound, rem)
			cardSup = min(inter.G.cardinalSup, length(inter.G.upperBound))

			changeG.cardRemoved = cardSup - inter.G.cardinalSup

			inter.G.cardinalSup += changeG.cardRemoved

			if !(isempty(rem))
				nbChange += 1
				changeVariable[nbChange] = inter.G

				if isFixed(inter.G)
					fix!(inter.G, changeG)
				end
			end

			for elt in rem
				push!(changeG.removed, elt)
			end

		end

		if !inter.H.isFixed
			rem = intersect(inter.H.upperBound, inter.G.lowerBound)

			setdiff!(inter.H.upperBound, rem)
			cardSup = min(inter.H.cardinalSup, length(inter.H.upperBound))

			changeH.cardRemoved = cardSup - inter.H.cardinalSup

			inter.H.cardinalSup += changeH.cardRemoved

			if !(isempty(rem))
				nbChange += 1
				changeVariable[nbChange] = inter.H

				if isFixed(inter.H)
					fix!(inter.H, changeH)
				end
			end

			for elt in rem
				push!(changeH.removed, elt)
			end

		end

		#Error
		stop = stop || !(inter.G.cardinalInf <= inter.G.cardinalSup)
		stop = stop || !(inter.H.cardinalInf <= inter.H.cardinalSup)

		changes = (changeG, changeH)
		for c in changes
			if !(isempty(c.added))
				#println("$inter ajoute $(c.added) à $(c.var)")
			end
			if !(isempty(c.removed))
				#println("$inter enlève $(c.removed) à $(c.var)")
			end
		end

		return changeVariable[1:nbChange], (changeG, changeH), stop
	end

end


text:
module QuadratureFactoriesTests

using Test
using Gridap

order = 5

p = Polytope(HEX_AXIS,HEX_AXIS)
quad = Quadrature(p,order)
@test isa(quad,TensorProductQuadrature)

p = Polytope(TET_AXIS,TET_AXIS,TET_AXIS)
quad = Quadrature(p,order)
@test isa(quad,DuffyQuadrature)

end # module


text:
using SemidefiniteProgramming

sdp = SparseSDP()

setobj!(sdp, 1, 1, 2, 1.0)

setcon!(sdp, 1, 1, 1, 1, 1.0)
setrhs!(sdp, 1, 1.0)

setcon!(sdp, 2, 1, 2, 2, 1.0)
setrhs!(sdp, 2, 1.0)

sol = solve(sdp, CSDP())

println(primalmatrix(sol))


text:
# Day 20: Trench Map
# Image is permuted for faster iteration
function load(::Day{20}, path)
    # function load(path)
    lines = readlines(path)
    alg = [c == '#' for c in lines[1]]
    img = BitMatrix(reduce(hcat, [[c == '#' for c in line] for line in lines[3:end]]))
    return alg, img
end

# Part 1 – count pixels in output image
# Offsets defining neighborhood. We need to "read across each column":
const CONV_OFFSETS = [
    CartesianIndex(-1, -1),
    CartesianIndex(0, -1),
    CartesianIndex(1, -1),
    CartesianIndex(-1, 0),
    CartesianIndex(0, 0),
    CartesianIndex(1, 0),
    CartesianIndex(-1, 1),
    CartesianIndex(0, 1),
    CartesianIndex(1, 1),
]

function pad_image(img; px=1, val=false)
    h, w = size(img)
    padrows = zeros(Bool, px, w)
    padcols = zeros(Bool, h + 2 * px, px)
    return hcat(padcols, vcat(padrows, img, padrows), padcols)
end

@inline getval(img, I::CartesianIndex, outer) = checkbounds(Bool, img, I) ? img[I] : outer

const POWS = 2 .^ collect(8:-1:0)
@inline neighbors2num(n) = sum(POWS .* n) + 1

function apply_algorithm!(out, img, outer, alg)
    @inbounds @simd for I in CartesianIndices(img) # not sure @simd does anything here...
        # Index the value of the algorithm `alg` based on the binary number
        # that results from the neighboring pixels.
        out[I] = alg[neighbors2num([
            getval(img, I + offset, outer) for offset in CONV_OFFSETS
        ])]
    end
    return out
end

function enhance_n_times(alg, img, n)
    img = pad_image(img; px=n)
    buff = similar(img)
    outer = false

    for _ in 1:n
        img .= apply_algorithm!(buff, img, outer, alg)
        # The "outer values" are flashing if `alg[1] = 1` and `alg[end] = 0`
        alg[1] == 1 && (outer = !outer)
    end
    return img
end

solve1(::Day{20}, data) = sum(enhance_n_times(data..., 2))

# Part 2 – Part 1 with more steps
solve2(::Day{20}, data) = sum(enhance_n_times(data..., 50))

testresult1(::Day{20}) = 35
testresult2(::Day{20}) = 3351


text:
using DataDeps

register(DataDep(
  "heart-disease-cleveland",
  "http://archive.ics.uci.edu/ml/datasets/Heart+Disease",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/heart-disease/processed.cleveland.data",
  "a74b7efa387bc9d108d7d0115d831fe9b414b29ae7124f331b622b4efa0427c8",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=14,
        feature_indices=1:13,
        categoric_indices=[3; 7; 11; 13],
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
## dg.jl

@testset "test simple_allowed_reactions allow_missing" begin

    @testset "some dgs" begin
        dgs = [101.1,15,-15,missing]
        threshold = -10.2

        @test [false,false,true,true]  == BioXP.simple_allowed_reactions(dgs,threshold,true)
        @test [false,false,true,false] == BioXP.simple_allowed_reactions(dgs,threshold,false)
    end

    @testset "value = threshold" begin
        dgs = [missing,missing,0,8,900]
        threshold = 8

        @test [true,true,true,true,false]   == BioXP.simple_allowed_reactions(dgs,threshold,true)
        @test [false,false,true,true,false] == BioXP.simple_allowed_reactions(dgs,threshold,false)
    end

end

@testset "test filter_reactions_by_dg" begin

    threshold = -75
    env_key = "11pH_300mM"
    rids = readids(["R02745","R05080","R00778","R08348","R06706"])
    rstructs = readmaster("data/master_with_dgs_reserialized.json")

    @testset "unbalanced only" begin
        allow_nothings=false
        allow_unbalanced=true
        allow_within_ci=false

        forward_allowed=[false,false,true,true,false]
        backward_allowed=[false,false,false,false,true]
        @test (forward_allowed, backward_allowed) == filter_reactions_by_dg(
                                                    threshold,
                                                    env_key,
                                                    rids,
                                                    rstructs,
                                                    allow_nothings=allow_nothings,
                                                    allow_unbalanced=allow_unbalanced,
                                                    allow_within_ci=allow_within_ci)

    end

    @testset "nothings and unbalanced" begin
        allow_nothings=true
        allow_unbalanced=true
        allow_within_ci=false

        forward_allowed=[true,true,true,true,false]
        backward_allowed=[true,true,false,false,true]
        @test (forward_allowed, backward_allowed) == filter_reactions_by_dg(
                                                    threshold,
                                                    env_key,
                                                    rids,
                                                    rstructs,
                                                    allow_nothings=allow_nothings,
                                                    allow_unbalanced=allow_unbalanced,
                                                    allow_within_ci=allow_within_ci)
    end

    @testset "nothings only" begin
        allow_nothings=true
        allow_unbalanced=false
        allow_within_ci=false

        forward_allowed=[false,false,true,true,false]
        backward_allowed=[false,false,false,false,false]
        @test (forward_allowed, backward_allowed) == filter_reactions_by_dg(
                                                    threshold,
                                                    env_key,
                                                    rids,
                                                    rstructs,
                                                    allow_nothings=allow_nothings,
                                                    allow_unbalanced=allow_unbalanced,
                                                    allow_within_ci=allow_within_ci)
    end

    @testset "ci only" begin
        allow_nothings=false
        allow_unbalanced=false
        allow_within_ci=true

        forward_allowed=[false,false,true,true,false]
        backward_allowed=[false,false,true,false,false]
        @test (forward_allowed, backward_allowed) == filter_reactions_by_dg(
                                                    threshold,
                                                    env_key,
                                                    rids,
                                                    rstructs,
                                                    allow_nothings=allow_nothings,
                                                    allow_unbalanced=allow_unbalanced,
                                                    allow_within_ci=allow_within_ci)

    end

    @testset "nothings and ci" begin #this doesn't change anything either
        allow_nothings=true
        allow_unbalanced=false
        allow_within_ci=true

        forward_allowed=[false,false,true,true,false]
        backward_allowed=[false,false,true,false,false]
        @test (forward_allowed, backward_allowed) == filter_reactions_by_dg(
                                                    threshold,
                                                    env_key,
                                                    rids,
                                                    rstructs,
                                                    allow_nothings=allow_nothings,
                                                    allow_unbalanced=allow_unbalanced,
                                                    allow_within_ci=allow_within_ci)
    end

    @testset "unbalanced and ci" begin #ci doesn't open any new reactions here
        allow_nothings=false
        allow_unbalanced=true
        allow_within_ci=true

        forward_allowed=[false,false,true,true,false]
        backward_allowed=[false,false,true,false,true]
        @test (forward_allowed, backward_allowed) == filter_reactions_by_dg(
                                                    threshold,
                                                    env_key,
                                                    rids,
                                                    rstructs,
                                                    allow_nothings=allow_nothings,
                                                    allow_unbalanced=allow_unbalanced,
                                                    allow_within_ci=allow_within_ci)
    end

    @testset "all allowed" begin #ci doesn't open any new reactions here
        allow_nothings=true
        allow_unbalanced=true
        allow_within_ci=true

        forward_allowed=[true,true,true,true,false]
        backward_allowed=[true,true,true,false,true]
        @test (forward_allowed, backward_allowed) == filter_reactions_by_dg(
                                                    threshold,
                                                    env_key,
                                                    rids,
                                                    rstructs,
                                                    allow_nothings=allow_nothings,
                                                    allow_unbalanced=allow_unbalanced,
                                                    allow_within_ci=allow_within_ci)
    end

    ## need to test where ci allows new reaction--change threshold

end

max_stars_repo_name: hbsmith/BioXP

max_issues_repo_name: hbsmith/BioXP

max_forks_repo_name: hbsmith/BioXP


text:
#=
    ^, pow, root

    exp, expm1, log, log1p, log2, log10, logbase,
    sin, sinpi, cos, cospi, tan, tanpi, cot, cotpi,
    sinh, cosh, tanh, coth,
    asin, acos, atan, asinh, acosh, atanh,
    sincos, sincospi, sinhcosh,
    sinc,
    gamma, lgamma, zeta
=#

for (op,cfunc) in ((:exp,:arb_exp), (:expm1, :arb_expm1),
    (:log,:arb_log), (:log1p, :arb_log1p),
    (:sin, :arb_sin), (:sinpi, :arb_sinpi), (:cos, :arb_cos), (:cospi, :arb_cospi),
    (:tan, :arb_tan), (:cot, :arb_cot),
    (:sinh, :arb_sinh), (:cosh, :arb_cosh), (:tanh, :arb_tanh), (:coth, :arb_coth),
    (:asin, :arb_asin), (:acos, :arb_asin), (:atan, :arb_atan),
    (:asinh, :arb_asinh), (:acosh, :arb_asinh), (:atanh, :arb_atanh),
    (:sinc, :arb_sinc),
    (:gamma, :arb_gamma), (:lgamma, :arb_lgamma), (:zeta, :arb_zeta),
    (:digamma, :arb_digamma), (:rgamma, :arb_rgamma)
    )
  @eval begin
    function ($op)(x::ArbFloat{P}) where {P}
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, P)
      z
    end
  end
end

function logbase(x::ArbFloat{P}, base::Int) where {P}
    b = UInt(abs(base))
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_log_base_ui), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, UInt, Int), z, x, b, P)
    z
end

log2(x::ArbFloat{P}) where {P} = logbase(x, 2)
log10(x::ArbFloat{P}) where {P} = logbase(x, 10)

for (op,cfunc) in ((:sincos, :arb_sin_cos), (:sincospi, :arb_sin_cos_pi), (:sinhcosh, :arb_sinh_cosh))
  @eval begin
    function ($op)(x::ArbFloat{P}) where {P}
        sz = initializer(ArbFloat{P})
        cz = initializer(ArbFloat{P})
        ccall(@libarb($cfunc), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), sz, cz, x, P)
        sz, cz
    end
  end
end

atan(a::ArbFloat{P}, b::ArbFloat{P}) where {P} = atan2(a, b)

function atan2(a::ArbFloat{P}, b::ArbFloat{P}) where {P}
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_atan2), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, a, b, P)
    z
end

for (op,cfunc) in ((:root, :arb_root_ui),)
  @eval begin
    function ($op)(x::ArbFloat{P}, y::UInt) where {P}
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, UInt, Int), z, x, y, P)
      z
    end
  end
end

for (op,cfunc) in ((:^,:arb_pow), (:pow,:arb_pow))
  @eval begin
    function ($op)(x::I, y::ArbFloat{P}) where {P,I <: Integer}
      xx = ArbFloat{P}(x)
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, xx, y, P)
      return z
    end
    function ($op)(x::ArbFloat{P}, y::I) where {P,I <: Integer}
      sy,ay = signbit(y), abs(y)
      yy = ArbFloat{P}(ay)
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, yy, P)
      return sy ? inv(z) : z
    end
    function ($op)(x::R, y::ArbFloat{P}) where {P,R <: Rational}
      xx = ArbFloat{P}(x)
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, xx, y, P)
      return z
    end
    function ($op)(x::ArbFloat{P}, y::R) where {P,R <: Rational}
      yy = ArbFloat{P}(y)
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, yy, P)
      return z
    end
    function ($op)(x::R, y::ArbFloat{P}) where {P,R <: Real}
      xx = ArbFloat{P}(x)
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, xx, y, P)
      return z
    end
    function ($op)(x::ArbFloat{P}, y::R) where {P,R <: Real}
      yy = ArbFloat{P}(y)
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, yy, P)
      return z
    end
    function ($op)(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
      z = initializer(ArbFloat{P})
      ccall(@libarb($cfunc), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, y, P)
      z
    end
    ($op)(x::ArbFloat{P}, y::ArbFloat{Q}) where {P,Q} = ($op)(promote(x,y)...)
  end
end

root(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = pow(x, inv(y))
root(x::Integer, y::ArbFloat{P}) where {P} = pow(ArbFloat{P}(x), inv(y))
function root(x::ArbFloat{P}, y::Integer) where {P}
   return (
     if y>=0
       yy = UInt64(y)
      z = initializer(ArbFloat{P})
       ccall(@libarb(arb_root_ui), Cvoid, (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, UInt64, Int), z, x, yy, P)
       z
    else
      pow(ArbFloat{P}(x), inv(y))
    end )
end


text:
#!/usr/bin/julia
# Taken from https://github.com/phrb/ExperimentalDesign.jl/blob/master/examples/Screening%20with%20Plackett-Burman%20Designs.ipynb

using Pkg

#Pkg.add("StatsModels")

using ExperimentalDesign, StatsModels, GLM, DataFrames, Distributions, Random

include("gcc.jl")

design = PlackettBurman(5)
design.matrix

println(design)

exit()

design.dummy_factors

design.factors

function y(x)
    return run_gcc(x)
end

Random.seed!(192938)

design.matrix[!, :response] = y.(eachrow(design.matrix[:, collect(design.factors)]))
design.matrix

lm(design.formula, design.matrix)

lm(term(:response) ~ sum(term.(design.factors)), design.matrix)

Random.seed!(8418172)

random_design_generator = RandomDesign(DiscreteNonParametric([-1, 1], [0.5, 0.5]), 6)
random_design = rand(random_design_generator, 8)

random_design[!, :response] = y.(eachrow(random_design[:, :]))
random_design

lm(random_design_generator.formula, random_design)

Random.seed!(2989476)

factorial_design = FullFactorial(fill([-1, 1], 6), explicit = true)
factorial_design.matrix[!, :response] = y.(eachrow(factorial_design.matrix[:, :]))

println(factorial_design)

lm(factorial_design.formula, factorial_design.matrix)

max_stars_repo_name: giulianobelinassi/LTO-Timing-Analysis

max_issues_repo_name: giulianobelinassi/LTO-Timing-Analysis

max_forks_repo_name: giulianobelinassi/LTO-Timing-Analysis


text:
using CUDA
using CUDA.CUSPARSE
using Printf
using FiniteDiff
using ForwardDiff
using KernelAbstractions
using LinearAlgebra
using Random
using SparseArrays
using Test
using TimerOutputs
using ExaPF
import ExaPF: PowerSystem, AutoDiff

const PS = PowerSystem

@testset "Polar formulation" begin
    datafile = joinpath(dirname(@__FILE__), "..", "..", "data", "case9.m")
    tolerance = 1e-8
    pf = PS.PowerNetwork(datafile)

    if has_cuda_gpu()
        ITERATORS = zip([CPU(), CUDADevice()], [Array, CuArray])
    else
        ITERATORS = zip([CPU()], [Array])
    end

    @testset "Initiate polar formulation on device $device" for (device, M) in ITERATORS
        polar = PolarForm(pf, device)
        # Test printing
        println(devnull, polar)

        b = bounds(polar, State())
        b = bounds(polar, Control())

        # Test getters
        nᵤ = ExaPF.get(polar, NumberOfControl())
        nₓ = ExaPF.get(polar, NumberOfState())
        ngen = get(polar, PS.NumberOfGenerators())
        nbus = get(polar, PS.NumberOfBuses())

        # Get initial position
        x0 = ExaPF.initial(polar, State())
        u0 = ExaPF.initial(polar, Control())

        @test length(u0) == nᵤ
        @test length(x0) == nₓ
        for v in [x0, u0]
            @test isa(v, M)
        end

        cache = ExaPF.get(polar, ExaPF.PhysicalState())
        @testset "NetworkState cache" begin
            # By defaut, values are equal to 0
            @test iszero(cache)
            @test isa(cache.vmag, M)
            @test cache.bus_gen == polar.indexing.index_generators
            # Transfer control u0 inside cache
            ExaPF.transfer!(polar, cache, u0)
            # Test that all attributes have valid length
            @test length(cache.vang) == length(cache.vmag) == length(cache.pinj) == length(cache.qinj) == nbus
            @test length(cache.pg) == length(cache.qg) == length(cache.bus_gen) == ngen
            @test length(cache.dx) == length(cache.balance) == nₓ
            # Test setters
            ## Buses
            values = similar(x0, nbus)
            fill!(values, 1)
            ExaPF.setvalues!(cache, PS.VoltageMagnitude(), values)
            @test cache.vmag == values
            ExaPF.setvalues!(cache, PS.VoltageAngle(), values)
            @test cache.vang == values
            ExaPF.setvalues!(cache, PS.ActiveLoad(), values)
            ExaPF.setvalues!(cache, PS.ReactiveLoad(), values)
            # Power generations are still equal to 0, so we get equality
            @test cache.pinj == -values  # Pinj = 0 - Pd
            @test cache.qinj == -values  # Qinj = 0 - Qd
            ## Generators
            vgens = similar(x0, ngen)
            fill!(vgens, 2.0)
            ExaPF.setvalues!(cache, PS.ActivePower(), vgens)
            ExaPF.setvalues!(cache, PS.ReactivePower(), vgens)
            @test cache.pg == vgens
            @test cache.qg == vgens
            genbus = polar.indexing.index_generators
            @test cache.pinj[genbus] == vgens - values[genbus]  # Pinj = Cg*Pg - Pd
            @test cache.qinj[genbus] == vgens - values[genbus]  # Qinj = Cg*Qg - Pd
            # After all these operations, values become non-trivial
            @test !iszero(cache)
        end
        # Reset to default values before going any further
        ExaPF.init_buffer!(polar, cache)
        @test !iszero(cache)

        @testset "Polar model API" begin
            xₖ = copy(x0)
            # Init AD factory
            jx = AutoDiff.Jacobian(polar, ExaPF.power_balance, State())

            # Test powerflow with cache signature
            conv = powerflow(polar, jx, cache, NewtonRaphson(tol=tolerance))
            # Get current state
            ExaPF.get!(polar, State(), xₖ, cache)
            # Refresh power of generators in cache
            for Power in [PS.ActivePower, PS.ReactivePower]
                ExaPF.update!(polar, PS.Generators(), Power(), cache)
            end

            # Bounds on state and control
            u_min, u_max = ExaPF.bounds(polar, Control())
            x_min, x_max = ExaPF.bounds(polar, State())

            @test isequal(u_min, [0.9, 0.9, 0.9, 0.1, 0.1])
            @test isequal(u_max, [1.1, 1.1, 1.1, 3.0, 2.7])
            @test isequal(x_min, [-Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9])
            @test isequal(x_max, [Inf, Inf, Inf, Inf, Inf, Inf, Inf, Inf, 1.1, 1.1, 1.1, 1.1, 1.1, 1.1])

            # Test callbacks
            ## Power Balance
            cons = cache.balance
            ExaPF.power_balance(polar, cons, cache)
            # As we run powerflow before, the balance should be below tolerance
            @test norm(cons, Inf) < tolerance

            ## Cost Production
            c2 = ExaPF.cost_production(polar, cache)
            @test isa(c2, Real)

            ## Inequality constraint
            for cons_function in [
                ExaPF.voltage_magnitude_constraints,
                ExaPF.active_power_constraints,
                ExaPF.reactive_power_constraints,
                ExaPF.flow_constraints,
                ExaPF.power_balance,
            ]
                m = ExaPF.size_constraint(polar, cons_function)
                @test isa(m, Int)
                g = M{Float64, 1}(undef, m) # TODO: this signature is not great
                fill!(g, 0)
                cons_function(polar, g, cache)

                g_min, g_max = ExaPF.bounds(polar, cons_function)
                @test length(g_min) == m
                @test length(g_max) == m
                # Are we on the correct device?
                @test isa(g_min, M)
                @test isa(g_max, M)
                # Test constraints are consistent
                @test g_min <= g_max
            end
        end
    end
end


text:
@testset "generate" begin
    @testset "rand_graph" begin
        n, m = 10, 20
        m2 = m ÷ 2
        x = rand(3, n)
        e = rand(4, m2)

        g = rand_graph(n, m, ndata=x, edata=e, graph_type=GRAPH_T)
        @test g.num_nodes == n
        @test g.num_edges == m
        @test g.ndata.x === x
        if GRAPH_T == :coo
            s, t = edge_index(g)
            @test s[1:m2] == t[m2+1:end]
            @test t[1:m2] == s[m2+1:end]
            @test g.edata.e[:,1:m2] == e
            @test g.edata.e[:,m2+1:end] == e
        end

        g = rand_graph(n, m, bidirected=false, seed=17, graph_type=GRAPH_T)
        @test g.num_nodes == n
        @test g.num_edges == m

        g2 = rand_graph(n, m, bidirected=false, seed=17, graph_type=GRAPH_T)
        @test edge_index(g2) == edge_index(g)
    end
end


text:
"""
Function returning the name of a where parameter
"""
function get_param_name(expr) :: Symbol
  @capture(expr, arg_<:arg_type_) && return arg
  @capture(expr, arg_) && return arg
end

"""
Function returning the arguments of a function definition
"""
function get_args(func_def::Dict)
  arg_dict = Dict{Symbol, Any}()
  arg_list = Vector{Symbol}()
  for arg in (func_def[:args]..., func_def[:kwargs]...)
    arg_def = splitarg(arg)
    push!(arg_list, arg_def[1])
    arg_dict[arg_def[1]] = arg_def[3] ? Tuple : arg_dict[arg_def[1]] = arg_def[2]
  end
  arg_list, arg_dict
end

"""
Function returning the slots of a function definition
"""
function get_slots(func_def::Dict, args::Dict{Symbol, Any}, mod::Module) :: Dict{Symbol, Any}
  slots = Dict{Symbol, Any}()
  func_def[:name] = gensym()
  func_def[:args] = (func_def[:args]..., func_def[:kwargs]...)
  func_def[:kwargs] = []
  body = func_def[:body]
  func_def[:body] = postwalk(transform_yield, func_def[:body])
  func_expr = combinedef(func_def) |> flatten
  @eval(mod, @noinline $func_expr)
  code_data_infos = @eval(mod, begin using ResumableFunctions; ResumableFunctions.my_code_typed($(func_def[:name])) end )
  for (code_info, slottypes) in code_data_infos
    for (i, slotname) in enumerate(code_info.slotnames)
      slots[slotname] = slottypes[i]
    end
  end
  for (argname, argtype) in args
    slots[argname] = argtype
  end
  postwalk(x->remove_catch_exc(x, slots), func_def[:body])
  postwalk(x->make_arg_any(x, slots), body)
  for (key, val) in slots
    if val == Union{}
      slots[key] = Any
    end
  end
  delete!(slots, Symbol("#temp#"))
  delete!(slots, Symbol("#unused#"))
  delete!(slots, Symbol("#self#"))
  slots
end

function my_code_typed(@nospecialize(f), @nospecialize(types=Tuple))
  types = Core.Compiler.to_tuple_type(types)
  asts = []
  world = ccall(:jl_get_world_counter, UInt, ())
  params = Core.Compiler.Params(world)
  for x in Core.Compiler._methods(f, types, -1, world)
      meth = Core.Compiler.func_for_method_checked(x[3], types)
      (code, slottypes) = my_typeinf_code(meth, x[1], x[2], params)
      code === nothing && error("inference not successful")
      push!(asts, code => slottypes)
  end
  return asts
end

function my_typeinf_code(method::Method, @nospecialize(atypes), sparams::Core.Compiler.SimpleVector, params::Core.Compiler.Params)
  code = Core.Compiler.code_for_method(method, atypes, sparams, params.world)
  code === nothing && return (nothing, Any)
  ccall(:jl_typeinf_begin, Cvoid, ())
  result = Core.Compiler.InferenceResult(code)
  frame = Core.Compiler.InferenceState(result, false, params)
  frame === nothing && return (nothing, Any)
  Core.Compiler.typeinf(frame)
  ccall(:jl_typeinf_end, Cvoid, ())
  frame.inferred || return (nothing, Any)
  return (frame.src, frame.slottypes)
end

"""
Function removing the `exc` symbol of a `catch exc` statement of a list of slots.
"""
function remove_catch_exc(expr, slots::Dict{Symbol, Any})
  @capture(expr, (try body__ catch exc_; handling__ end) | (try body__ catch exc_; handling__ finally always__ end)) && delete!(slots, exc)
  expr
end

"""
Function changing the type of a slot `arg` of a `arg = @yield ret` or `arg = @yield` statement to `Any`.
"""
function make_arg_any(expr, slots::Dict{Symbol, Any})
  @capture(expr, arg_ = ex_) || return expr
  _is_yield(ex) || return expr
  slots[arg] = Any
  expr
end

"""
Function returning the args for the type construction.
"""
function make_args(func_def::Dict)
  args=[]
  for arg in (func_def[:args]..., func_def[:kwargs]...)
    arg_def = splitarg(arg)
    push!(args, combinearg(arg_def[1], arg_def[2], false, arg_def[4]))
  end
  (args...,)
end

"""
Function checking the use of a return statement with value
"""
function hasreturnvalue(expr)
  @capture(expr, return val_) || return expr
  (val == :nothing || val == nothing) && return expr
  @warn "@resumable function contains return statement with value!"
  expr
end


text:
using DelimitedFiles

function ReadData(fname::AbstractString)
    return readdlm(fname, '\n', String)
end

function Part1(data::Array)
    nDigits = length(data[1])
    nLines = length(data)

    # convert raw data into binary
    D = map(x -> parse(Int,x,base=2), data)

    # sums of the digits = number of 1s in each column
    S = [sum((i >> j) & 1 for i in D) for j=0:nDigits-1]

    # compute gamma
    gamma = sum(Int(S[i]>nLines/2) << (i-1) for i=1:nDigits)

    # compute epsilon
    epsilon = ~gamma & sum(1<<i for i=0:nDigits-1)

    return gamma*epsilon
end

function Part2(data::Array)
    nDigits = length(data[1])
    nLines = length(data)

    # convert raw data into binary, do this twice for O2 and CO2 ratings
    D_O2  = map(x -> parse(Int,x,base=2), data)
    D_CO2 = map(x -> parse(Int,x,base=2), data)

    j = 1 # the bit we're currently working on
    while length(D_O2) > 1
        # determine most common bit
        cBit = Int(sum((i >> (nDigits-j)) & 1 for i in D_O2) >= length(D_O2)/2)
        # filter out numbers which don't have same common first bit
        D_O2 = filter(x -> ((x >> (nDigits-j)) & 1) == cBit, D_O2)

        j += 1
    end

    j = 1 # the bit we're currently working on
    while length(D_CO2) > 1
        # determine least common bit
        cBit = Int(sum((i >> (nDigits-j)) & 1 for i in D_CO2) < length(D_CO2)/2)
        # filter out numbers which don't have same common first bit
        D_CO2 = filter(x -> ((x >> (nDigits-j)) & 1) == cBit, D_CO2)

        j += 1
    end

    return D_O2[1]*D_CO2[1]
end

# test data
#data = ["00100", "11110", "10110", "10111", "10101", "01111", "00111", "11100",
#        "10000", "11001", "00010", "01010"]

data = ReadData("../data/day3.txt")

println(Part1(data))

println(Part2(data))