text:
@testset "fma" begin

  @test HILO(fma(Float64(pi), Float64(phi), Double64(gamma))) == (5.660419357216792, 4.1344240638440936e-16)
  @test HILO(fma(Float64(pi), Double64(phi), Float64(gamma))) == (5.660419357216792, 2.477303893634162e-16)
  @test HILO(fma(Double64(pi), Float64(phi), Float64(gamma))) == (5.660419357216793, -2.7164108363986695e-16)
  @test HILO(fma(Double64(pi), Double64(phi), Float64(gamma))) == (5.660419357216793, -4.422960158132908e-16)
  @test HILO(fma(Double64(pi), Float64(phi), Double64(gamma))) == (5.660419357216793, -2.765839987922976e-16)
  @test HILO(fma(Double64(pi), Double64(phi), Double64(gamma))) == (5.660419357216792, 4.4093948873440384e-16)

  @test HILO(muladd(Float64(pi), Float64(phi), Double64(gamma))) == (5.660419357216792, 4.1344240638440936e-16)
  @test HILO(muladd(Float64(pi), Double64(phi), Float64(gamma))) == (5.660419357216792, 2.477303893634162e-16)
  @test HILO(muladd(Double64(pi), Float64(phi), Float64(gamma))) == (5.660419357216793, -2.7164108363986695e-16)
  @test HILO(muladd(Double64(pi), Double64(phi), Float64(gamma))) == (5.660419357216793, -4.422960158132908e-16)
  @test HILO(muladd(Double64(pi), Float64(phi), Double64(gamma))) == (5.660419357216793, -2.765839987922976e-16)
  @test HILO(muladd(Double64(pi), Double64(phi), Double64(gamma))) == (5.660419357216792, 4.4093948873440384e-16)

end


text:
f(theta)               = simloglik!(grad, info_matrix, tmpvars, model, theta, data, false)
g!(grad, theta)        = simloglik!(grad, info_matrix, tmpvars, model, theta, data, true)
h!(info_matrix, theta) = simloglik!(grad, info_matrix, tmpvars, model, theta, data, true)
fg! = g!
newton!(f, g, h, theta, x0)

ext: jl


text:
export @checked, @debug_ccall

"""
    @checked function foo(...)
        rv = ...
        return rv
    end

Macro for wrapping a function definition returning a status code. Two versions of the
function will be generated: `foo`, with the function body wrapped by an invocation of the
`@check` macro (to be implemented by the caller of this macro), and `unsafe_foo` where no
such invocation is present and the status code is returned to the caller.
"""
macro checked(ex)
    # parse the function definition
    @assert Meta.isexpr(ex, :function)
    sig = ex.args[1]
    @assert Meta.isexpr(sig, :call)
    body = ex.args[2]
    @assert Meta.isexpr(body, :block)

    # generate a "safe" version that performs a check
    safe_body = quote
        @check $body
    end
    safe_sig = Expr(:call, sig.args[1], sig.args[2:end]...)
    safe_def = Expr(:function, safe_sig, safe_body)

    # generate a "unsafe" version that returns the error code instead
    unsafe_sig = Expr(:call, Symbol("unsafe_", sig.args[1]), sig.args[2:end]...)
    unsafe_def = Expr(:function, unsafe_sig, body)

    return esc(:($safe_def, $unsafe_def))
end

macro debug_ccall(target, rettyp, argtyps, args...)
    @assert Meta.isexpr(target, :tuple)
    f, lib = target.args

    quote
        # get the call target, as e.g. libcuda() triggers initialization, even though we
        # can't use the result in the ccall expression below as it's supposed to be constant
        $(esc(target))

        print($f, '(')
        for (i, arg) in enumerate(($(map(esc, args)...),))
            i > 1 && print(", ")
            render_arg(stdout, arg)
        end
        print(')')
        rv = ccall($(esc(target)), $(esc(rettyp)), $(esc(argtyps)), $(map(esc, args)...))
        println(" = ", rv)
        for (i, arg) in enumerate(($(map(esc, args)...),))
            if arg isa Base.RefValue
                println(" $i: ", arg[])
            end
        end
        rv
    end
end

render_arg(io, arg) = print(io, arg)
render_arg(io, arg::Union{<:Base.RefValue, AbstractArray}) = summary(io, arg)


text:
"""
    uncamel(s)

Convert a camel case string to a underscored spaced / lowercased string.

# Examples
```julia
julia> uncamel("CamelCase")
camel_case
```
"""
function uncamel(s)
    final = ""
    for item in s
        if item == uppercase(item)
            final *= ("_" * string(lowercase(item)))
        else
            final *= string(item)
        end
    end
    if final[1] == '_'
        final = final[2:end]
    end
    final
end

"""
    cleanup(s)

Cleanup string from prefix

# Examples
```julia
julia> cleanup("optInCamelCase")
CamelCase
```

"""
function cleanup_prefix(s)
    for prefix in ["in", "optIn", "outReal", "outInteger"]
        if startswith(s, prefix)
            if s != prefix
                N = length(prefix)
                s = s[N+1:end]
            else
                d = Dict(
                    "outInteger"=>"Integer",
                    "outReal"=>"Real",
                )
                s = d[s]
            end
        end
    end
    s
end

"""
    fix_varname(s)

Lower case and replace spaces by underscores

# Examples
```julia
julia> fix_varname("Time Period")
time_period
```

"""
function fix_varname(s)
    s = lowercase(s)
    s = replace(s, " " => "_")
    s = replace(s, "-" => "_")
    s
end

"""
    replace_var(s)

Replace a variable with a more appropriate name (using a dict)

# Examples
```julia
julia> replace_var("Open")
:price_open
```
"""
function replace_var(s::AbstractString)
    d = Dict{String, Symbol}(
        "Open" => :price_open,
        "High" => :price_high,
        "Low" => :price_low,
        "Close" => :price_close,
        "Volume" => :volume,
    )
    get(d, s, Symbol(s))
end

"""
    replace_output(s)

Replace a column name with a more appropriate name (using a dict)

# Examples
```julia
julia> replace_output("outReal")
:Value
```
"""
function replace_output(s::AbstractString)
    d = Dict{String, Symbol}(
        "outInteger" => :Value,
        "outReal" => :Value,
    )
    get(d, s, Symbol(s))
end


text:
module Paulis

import ..AbstractOps
import ..AbstractOps: op_index, unsafe_op, symplectic, _AbstractOp
using ..AbstractPaulis: AbstractPauli

export Pauli, unsafe_op
import IsApprox

"""
    struct Pauli <: AbstractPauli

An implementation of `AbstractPauli`.
"""
struct Pauli <: AbstractPauli{Complex{Int}}
    hi::Bool
    lo::Bool
end

####
#### Constructors
####

const I = Pauli(0, 0)
const X = Pauli(0, 1)
const Y = Pauli(1, 0)
const Z = Pauli(1, 1)

# Indexing into static Tuple seems faster
"""
    Pauli(ind::Union{Integer, Symbol, AbstractString, AbstractChar})

Return a `Pauli` indexed by `[0, 3]` or a representation of `[I, X, Y, Z]`.
"""
Pauli(ind::Integer)::Pauli = (I, X, Y, Z)[ind + 1]

## The unsafe_op, unsafe_pauli thing does not seem to be more efficient.

@inline AbstractOps.unsafe_op(::Type{Pauli}, ind::Integer) = unsafe_pauli(ind)

@inline function unsafe_pauli(ind::Integer)::Pauli
    return  @inbounds (I, X, Y, Z)[ind + 1]
end

Pauli(s::Union{Symbol, AbstractString, AbstractChar}) = _AbstractOp(Pauli, s)

Base.copy(p::Pauli) = p

op_index(p::Pauli) = 2 * p.hi + p.lo

function symplectic(p::Pauli)
    if p.hi == 1
        if p.lo == 1
            return (false, true)
        else
            return (true, true)
        end
    else
        if p.lo == 1
            return (true, false)
        else
            return (false, false)
        end
    end
end

####
#### Compare / predicates
####

Base.:(==)(p1::Pauli, p2::Pauli) = p1.hi == p2.hi && p1.lo == p2.lo

## This is maybe not really necessary, but this is stricter than the fallback method.
## Also, this would have prevented a perf regression bug.
Base.isone(p::Pauli) = p === one(Pauli)

####
#### Algebra / mathematical operations
####

"""
    *(p1::Pauli, p2::Pauli)

Multiplication that returns a `Pauli`, but ignores the phase.
"""
@inline Base.:*(p1::Pauli, p2::Pauli) = Pauli(p1.hi ⊻ p2.hi, p1.lo ⊻ p2.lo)

end # module Paulis


text:
# Statistical Models

abstract StatisticalModel

coef(obj::StatisticalModel) = error("coef is not defined for $(typeof(obj)).")
coeftable(obj::StatisticalModel) = error("coeftable is not defined for $(typeof(obj)).")
confint(obj::StatisticalModel) = error("coefint is not defined for $(typeof(obj)).")
deviance(obj::StatisticalModel) = error("deviance is not defined for $(typeof(obj)).")
loglikelihood(obj::StatisticalModel) = error("loglikelihood is not defined for $(typeof(obj)).")
nobs(obj::StatisticalModel) = size(model_response(obj), 1)
stderr(obj::StatisticalModel) = sqrt(diag(vcov(obj)))
vcov(obj::StatisticalModel) = error("vcov is not defined for $(typeof(obj)).")

abstract RegressionModel <: StatisticalModel

residuals(obj::RegressionModel) = error("residuals is not defined for $(typeof(obj)).")
model_response(obj::RegressionModel) = error("model_response is not defined for $(typeof(obj)).")
predict(obj::RegressionModel) = error("predict is not defined for $(typeof(obj)).")
predict!(obj::RegressionModel) = error("predict! is not defined for $(typeof(obj)).")


text:
using LightGraphs
using SimpleWeightedGraphs
using StaticArrays
using Distances
using NearestNeighbors

function getindexmatrix(dists)
    M = fill(0, size(dists))
    for (i, (u, v)) in enumerate(getsortededges(dists))
        M[u, v] = i
        M[v, u] = i
    end
    M
end

struct PersistenceState{T}
    # Precomputed properties.
    dists       ::Matrix{T}
    parents     ::Matrix{Int}
    indexmatrix ::Matrix{Int}
    triangles   ::EquilateralIterator{T}
    # Internal state
    reduced     ::Vector{BitSet}
    generators  ::Vector{Vector{Int}}
    diameters   ::Vector{T}
    # Buffers
    σ           ::BitSet
    cycle       ::Vector{Int}
end

coverdiameter(g::AbstractGraph) = maximum(weights(g))
coverdiameter(g::GeodesicComplex) = 2g.radius

function PersistenceState(gc::AbstractGraph)
    d = coverdiameter(gc)
    # todo: split connected components
    floyd = floyd_warshall_shortest_paths(gc)
    dists = floyd.dists
    parents = floyd.parents
    indexmatrix = getindexmatrix(dists)
    triangles = equilaterals(dists, d)

    reduced = fill(BitSet(), binomial(nv(gc), 2))

    PersistenceState(dists, parents, indexmatrix, triangles,
                     reduced, Vector{Int}[], eltype(dists)[], BitSet(), Int[])
end

Base.show(io::IO, st::PersistenceState) = print(io, "PersistenceState")

"""
    low(σ::BitSet)

Get the largest index in `σ`.
"""
low(σ::BitSet) = length(σ) > 0 ? last(σ) : 0

"""
    reduce!(st::PersistenceState)

Use the standard reduction algorithm to reduce `st.σ`.
"""
function reduce!(st::PersistenceState)
    l = low(st.σ)
    while l ≠ 0 && !isempty(st.reduced[l])
        symdiff!(st.σ, st.reduced[l])
        l = low(st.σ)
    end
    l
end

"""
    movecycle!(st::PersistenceState, Δ)

Move cycle spanned by `Δ` into `st.cycle` and `st.σ`.
"""
function movecycle!(st::PersistenceState, (i, j, k))
    empty!(st.σ)
    empty!(st.cycle)
    movepath!(st, i, j)
    movepath!(st, j, k)
    movepath!(st, k, i)
end

"""
    movepath!(st::PersistenceState, u, v)

Move shortest path from `u` to `v` into `st.cycle` and `st.σ`.
"""
function movepath!(st::PersistenceState, u, v)
    while u ≠ v
        u′ = st.parents[v, u]
        push!(st.σ, st.indexmatrix[u, u′])
        push!(st.cycle, u)
        u = u′
    end
end

function isgeodesic(st, Δ, cycle)
    true
end

function processtriangle!(st::PersistenceState, triangle, diameter, showprogress)
    movecycle!(st, triangle)
    # Skip cycles that visit a node more than once.
    length(st.σ) == length(st.cycle) || return
    # Geodesic check
    isgeodesic(st, Δ, cycle) || return

    l = reduce!(st)
    if l ≠ 0
        st.reduced[l] = copy(st.σ)
        if true || length(st.cycle) > 3
            push!(st.generators, copy(st.cycle))
            push!(st.diameters, diameter)
        end
    end
end

function persistence(gc::AbstractGraph, showprogress = false)
    showprogress && println("Calculating intrinsic persistence...")
    st = PersistenceState(gc)

    for (i, (Δ, diam)) in enumerate(st.triangles)
        processtriangle!(st, Δ, diam, showprogress)
    end
    #st
    #filter(x->x≢nothing, map(x->formatresults(gc,x), st.generators))
    showprogress && println("Postprocessing...")
    res = []
    for g in st.generators
        α, m = contract(gc, g)
        if α ≡ nothing
            #println("splat at $(m)!")
        else
            push!(res, (g, α, m))
        end
    end
    res
end

# TODO: optimize
function singlecontract!(α, gc::GeodesicComplex{T}, cycle) where {T}
    changed = false
    ε = zero(T)
    for (i, v) in enumerate(α)
        N = inrange(gc.tree, points(gc, v), gc.radius)
        D = diststocycle(gc, N, cycle)
        ε, j = findmin(vec(maximum(D, dims = 2)))
        changed |= α[i] ≠ N[j]
        α[i] = N[j]
    end
    unique!(α)
    changed, ε
end

function singlecontract_noalloc!(α, gc::GeodesicComplex{T}, cycle) where {T}
    changed = false
    ε = zero(T)
    for (i, u) in enumerate(α)
        N = inrange(gc.tree, points(gc, u), gc.radius)
        # Spaghetti mapreduce avoids allocating distance matrices.
        ε, j = mapreduce(min, enumerate(N)) do (i, v)
            (maximum(evaluate(gc.metric, points(gc, v), landmarks(gc, w))
                     for w in cycle), i)
        end
        changed |= α[i] ≠ N[j]
        α[i] = N[j]
    end
    unique!(α)
    changed, ε
end

function contract(gc::GeodesicComplex{T}, cycle) where {T}
    α = gc.landmarks[cycle]
    ε = zero(T)
    changed = true

    while changed
        changed, ε = singlecontract!(α, gc, cycle)
    end
    length(α) ≥ 3 ? α : nothing, ε
end

diststocycle(gc::GeodesicComplex{T, D}, is, cycle) where {T, D} =
    pairwise(gc.metric,
             reshape(reinterpret(T, points(gc, is)), (D, length(is))),
             reshape(reinterpret(T, landmarks(gc, cycle)), (D, length(cycle))))

#=
struct Cycle
end

struct IntrinsicPersistenceResults{G<:GeodesicComplex}
    complex ::G
    generators ::Vector{...}
    contracted ::Vector{...}
end
=#


text:
nvilli = 30000
latencyvals = [2.7137,0.4568]
refracvals = [log(150),log(7)]
bumpvals = [4.0,3.0,2.5]
paramvals = vcat(latencyvals,refracvals,bumpvals)

variance1 = [10.0] #variance for the noise model

#specify the values that determine the priors on the parameters
latencylocation = (2.0,3.5) #uniform distribution with (low,high) values
latencyscale = (0.2,0.7) #uniform distribution with (low,high) values
refractorylocation = (4.0,6.0) #uniform distribution with (low,high) values
refractoryscale = (1.5,2.5) #uniform distribution with (low,high) values
bumpamplitude = (3.0,5.0) #uniform distribution with (low,high) values
bumpshape = (log(3.0),0.3) #lognormal distribution with (location,scale) values; variable can vary roughly between 2.0 and 4.0, but becomes increasingly penalized outside
bumpscale = (log(2.5),0.3) #lognormal distribution with (location,scale) values; variable can vary roughly between 1.5 and 3.5, but becomes increasingly penalized outside

photodatafolder = joinpath(splitdir(@__FILE__())[1],"../../data")
photodatatestfile = joinpath(photodatafolder,"naturallight.jld")
photons1 = photonsequence(photodatatestfile)
zerocurrent1 = data(:array,dataindex(photons1),zeros(Float32,numvalues(photons1)))

policy1 = policy(:photoreceptor,bump=:sample) #stochastic lognormal latency and refractory parameters, sampled bumpshape parameters (with lognormal priors, see above)
params1 = parameters(:photoreceptor,policy1,latencylocation,latencyscale,refractorylocation,refractoryscale,bumpamplitude,bumpshape,bumpscale)

numevaluations = 100
nummodels = 10

meanresults1 = zeros(numvalues(photons1),nummodels)
varresults1 = zeros(numvalues(photons1),nummodels)
for i = 1:nummodels
    println("Creating a model with a new distribution of photons over microvilli")
    @time model1 = model(:photoreceptor,params1,photons1,zerocurrent1,variance1,nvilli,policy1)
    println("Evaluating the model $numevaluations times")
    @time results1 = evaluate!(model1,paramvals,numevaluations)
    meanresults1[:,i] = mean(results1,2)
    varresults1[:,i] = var(results1,2)
end
mmresults1 = mean(meanresults1,2)
mvresults1 = mean(varresults1,2)

PyPlot.figure("functionality/photoreceptor")
PyPlot.clf()
PyPlot.subplot(221)
tindex = dataindex(photons1)
PyPlot.plot(tindex,datavalues(photons1))
PyPlot.xlabel("Time")
PyPlot.ylabel("Number of Photons")
PyPlot.xlim(tindex[1],tindex[end])
PyPlot.title("Photon Sequence")
PyPlot.grid("on")

PyPlot.subplot(223)
PyPlot.plot(tindex,meanresults1)
PyPlot.plot(tindex,mmresults1,linewidth=4)
PyPlot.xlabel("Time")
PyPlot.ylabel("Current (nA)")
PyPlot.xlim(tindex[1],tindex[end])
PyPlot.title("Average light-induced current")
PyPlot.grid("on")

PyPlot.subplot(224)
PyPlot.plot(tindex,sqrt(varresults1))
PyPlot.plot(tindex,mean(sqrt(varresults1),2),linewidth=4)
PyPlot.xlabel("Time")
PyPlot.ylabel("
of Current (nA)")
PyPlot.xlim(tindex[1],tindex[end])
PyPlot.title("Standard deviation of light-induced current")
PyPlot.grid("on")

modelcurrent1 = data(:array,dataindex(photons1),mmresults1)
variance1 = [mean(mvresults1)]

println("Variance obtained from model: $variance1")

println("Creating a model with experimentally-obtained macrocurrent and variance for the noise")
@time model2 = model(:photoreceptor,params1,photons1,modelcurrent1,variance1,nvilli,policy1)
println("Evaluating the model $numevaluations times with loglikelihood calculation")
loglikelihood1 = zeros(numevaluations)
@time for i=1:numevaluations
    l1 = evaluate!(model2,paramvals)
    loglikelihood1[i] = loglikelihood(model2,l1)
end
loglikelihood2 = zeros(numevaluations)
paramvals2 = vcat([2.6,0.4,log(145),log(6.5)],bumpvals)
@time for i=1:numevaluations
    l2 = evaluate!(model2,paramvals2)
    loglikelihood2[i] = loglikelihood(model2,l2)
end

PyPlot.subplot(222)
PyPlot.plot(1:numevaluations,loglikelihood1,label="ll1")
PyPlot.plot(1:numevaluations,loglikelihood2;label="ll2")
PyPlot.xlabel("Evaluation")
PyPlot.ylabel("Log-Likelihood")
PyPlot.xlim(1,numevaluations)
PyPlot.title("Log-Likelihood across evaluations")
PyPlot.grid("on")
PyPlot.legend(loc="upper right",fancybox="true")


text:
using MathOptInterface

const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

using Test

using JuMP
using SeaPearl

@testset "MOI_wrapper.jl" begin

    include("optimizer_accessors.jl")
    include("homemade_bridging.jl")
    include("variables.jl")
    include("constraints.jl")
    include("jump.jl")

    # @testset "Creating an optimizer" begin
    #     model = SeaPearl.Optimizer()
    #     @test MOI.get(model, MOI.SolverName()) == "SeaPearl Solver"
    # end

    # @testset "Giving parameters to the optimizer" begin
    #     model = SeaPearl.Optimizer()
    #     MOI.set(model, MOI.RawParameter("Test"), "test")
    #     @test model.options["Test"] == "test"
    #     MOI.empty!(model)
    #     @test MOI.is_empty(model)
    #     @test model.options["Test"] == "test"
    # end

    # @testset "Adding constrained variables" begin
    #     model = SeaPearl.Optimizer()
    #     MOI.add_constrained_variable(model, MOI.Interval(1, 4))
    #     MOI.add_constrained_variable(model, MOI.Interval(1, 4))
    #     MOI.add_constrained_variable(model, MOI.Interval(1, 4))
    #     MOI.add_constrained_variable(model, MOI.Interval(1, 4))
    #     @test Set(keys(model.cpmodel.variables)) == Set(["1", "2", "3", "4"])
    # end

    # @testset "Adding constraints" begin
    #     model = SeaPearl.Optimizer()
    #     MOI.add_constrained_variable(model, MOI.Interval(1, 4))
    #     MOI.add_constrained_variable(model, MOI.Interval(1, 4))
    #     MOI.add_constrained_variable(model, MOI.Interval(1, 4))
    #     MOI.add_constrained_variable(model, MOI.Interval(1, 4))

    #     # add new constraints
    #     MOI.add_constraint(model, MOI.SingleVariable(MOI.VariableIndex(1)), MOI.LessThan(2))
    #     MOI.add_constraint(model, MOI.SingleVariable(MOI.VariableIndex(2)), MOI.GreaterThan(2))
    #     MOI.add_constraint(model, MOI.SingleVariable(MOI.VariableIndex(4)), MOI.Interval(2, 3))

    #     # use fixPoint
    #     SeaPearl.fixPoint!(model.cpmodel, model.cpmodel.constraints)

    #     # test if it had effect on the variables' domains
    #     @test SeaPearl.maximum(model.cpmodel.variables[string(MOI.VariableIndex(1).value)].domain) == 2
    #     @test SeaPearl.minimum(model.cpmodel.variables[string(MOI.VariableIndex(2).value)].domain) == 2
    #     @test SeaPearl.minimum(model.cpmodel.variables[string(MOI.VariableIndex(4).value)].domain) == 2
    #     @test SeaPearl.maximum(model.cpmodel.variables[string(MOI.VariableIndex(4).value)].domain) == 3

    #     # add some new constraints again
    #     MOI.add_constraint(model, MOI.SingleVariable(MOI.VariableIndex(3)), MOI.EqualTo(2))
    #     MOI.add_constraint(model, MOI.VectorOfVariables([MOI.VariableIndex(1), MOI.VariableIndex(3)]), SeaPearl.VariablesEquality(false))

    #     # perform another fixPoint
    #     SeaPearl.fixPoint!(model.cpmodel, model.cpmodel.constraints)

    #     # new bunch of test
    #     @test SeaPearl.isbound(model.cpmodel.variables[string(MOI.VariableIndex(1).value)])
    #     @test SeaPearl.minimum(model.cpmodel.variables[string(MOI.VariableIndex(3).value)].domain) == 2
    #     @test SeaPearl.maximum(model.cpmodel.variables[string(MOI.VariableIndex(3).value)].domain) == 2
    # end

    ### Not working yet ###
    # @testset "JuMP interface" begin
    #     model = Model()
    #     set_optimizer(model, SeaPearl.Optimizer, bridge_constraints = true)

    #     @variable(model, 1 <= x[1:3] <= 4)
    #     @constraint(model, x[1] in SeaPearl.NotEqualTo(2))
    #     @constraint(model, x[1] in SeaPearl.NotEqualTo(3))
    #     @constraint(model, x[2] in MOI.EqualTo(4))
    #     @constraint(model, x[1:2] in SeaPearl.VariablesEquality(false))
    #     @constraint(model, [x[3], x[2]] in SeaPearl.VariablesEquality(false))
    #     @constraint(model, 2x[1] + 3x[2] in MOI.GreaterThan(2))
    # end

    # @testset "JuMP interface full" begin
    #     model = Model()
    #     set_optimizer(model, SeaPearl.Optimizer, bridge_constraints = true)

    #     @variable(model, 1 <= x[1:3] <= 4)
    #     @constraint(model, x[1] != 2)
    #     @constraint(model, x[1] != 3)
    #     @constraint(model, x[2] == 4)
    #     @constraint(model, x[1:2] in SeaPearl.VariablesEquality(false))
    #     @constraint(model, [x[3], x[2]] in SeaPearl.VariablesEquality(false))
    #     @constraint(model, 2x[1] + 3x[2] >= 2)

    #     println(model)
    # end

end


text:
module StructuralInheritance

export @protostruct, totuple

const FieldType = Union{Symbol,Expr}
const SymbolTuple = Tuple{Vararg{Symbol,N}} where {N}

#Stores prototype field definitions
const fieldBacking = IdDict{Type, Vector{FieldType}}()

#prototype -> self
#concrete -> prototype
const shadowMap = IdDict{Type,Type}()

#store parametric type information
const parameterMap = IdDict{Type,Vector{FieldType}}()

const mutabilityMap = IdDict{Type,Bool}()

"""
Creates an abstract type with the given name
"""
abstracttype(name) = :(abstract type $name end)

"""
attaches parameters to a name
"""
function addparams(name,params)
    if length(params) == 0
        name
    else
        :($name{$(params...)})
    end
end

"""
returns an array with only the field definitions
"""
function filtertofields(unfiltered)
    filter(x->x isa Symbol || (x isa Expr && x.head == :(::)),unfiltered)
end

"""
flattens the scope of the fields
"""
flattenfields(x) = FieldType[]
flattenfields(x::Symbol) = FieldType[x]
function flattenfields(x::Expr)::Vector{FieldType}
    if x.head == :block
        vcat(FieldType[],flattenfields.(x.args)...)
    else
        FieldType[x]
    end
end

isfunctiondefinition(x) = false
isfunctiondefinition(x::Expr) = (x.head == :(=) || x.head == :function)

"""
returns array with only the field constructors
"""
function extractconstructors(_quote)
    filter(isfunctiondefinition,_quote.args[3].args)
end

"""
extracts the fields from a struct definition
"""
extractfields(leaf) = filtertofields(flattenfields(leaf.args[3]))

function newnames(nameNode::Symbol,module_,prefix)
    protoName = Symbol(prefix,nameNode)
    return (:($nameNode <: $protoName),protoName,nameNode,protoName)
end
function newnames(nameNode,module_,prefix)
    """
    handle inheritence conversions
    """
    function rectify(@nospecialize(x))::FieldType
        val = module_.eval(deparametrize(x))::Type
        if isabstracttype(val)
            x
        elseif haskey(shadowMap,val)
            addparams(:($(qualifyname(shadowMap[val]))),getparameters(x))
        else
            throw("inheritence from concrete types is limited to those defined by @protostruct, $val not found")
        end
    end

    if nameNode.head == :<:
        inheritFrom = rectify(nameNode.args[2])::Union{Type,FieldType}
        structHead = deepcopy(nameNode.args[1])::FieldType

        if isparametric(structHead)
            protoName = nameNode
            protoChild = protoName.args[1]::Expr
            protoChild.args[1] = Symbol(prefix,structHead.args[1])
            protoName.args[2] = inheritFrom
            return (:( $(structHead) <: $(detypevar(protoChild))),
                    protoName,
                    structHead,
                    protoChild)

        elseif structHead isa Symbol
            protoName = deepcopy(nameNode)
            lightProto = Symbol(prefix,structHead)
            protoName.args[1]  = lightProto
            protoName.args[2] = inheritFrom
            return (:( $(structHead) <: $(lightProto)),
                    protoName,
                    structHead,
                    lightProto)
        end
    end

    if isparametric(nameNode)
        protoName = deepcopy(nameNode)
        protoName.args[1] = Symbol(prefix,nameNode.args[1]::Symbol)
        return (:( $(nameNode) <: $(detypevar(protoName))),
                protoName,
                nameNode,
                protoName)
    end

    throw("structure of strucure name not identified")
end

"""
from Foo{A<:B}
returns Foo{A}
"""
detypevar(x) = x
function detypevar(x::Expr)
    if x.head == :<:
        detypevar(x.args[1]::FieldType)
    else
        Expr(x.head,detypevar.(x.args)...)::Expr
    end
end

isparametric(x) = false
isparametric(x::Expr) = x.head == :curly || any(isparametric,x.args)

ispath(x) = false
ispath(x::Expr) = x.head == :.

iscontainerlike(x) = false
iscontainerlike(x::Expr) = (x.head in (:vect,:hcat,:row,
                                       :vcat, :call,
                                       :tuple,:curly,:macrocall))

function getpath(x)
    oldpath = Symbol[]
    while ispath(x)
        push!(oldpath,(x.args[2]::QuoteNode).value)
        x = x.args[1]::FieldType
    end
    oldpath = push!(oldpath,x)
    oldpath[2:end], oldpath[1]
end

function get2parameters(x)
    if x isa Expr && x.head == :<:
        (stablegetparameters(x.args[1]),getparameters(x.args[2]))
    else
        (stablegetparameters(x),Symbol[])
    end
end

getparameters(x)::Vector = isparametric(x) ? detypevar.(x.args[2:end]) : Symbol[]
stablegetparameters(x)::Vector{Symbol} = isparametric(x) ? detypevar.(x.args[2:end]) : Symbol[]

function getfieldnames(x)
    f(x::Symbol) = x
    f(x::Expr) = x.args[1]::Symbol
    f.(x)
end

"""
creates AST for expanding a struct into a tuple with the fields given
"""
function tupleexpander(x,fields)
    Expr(:tuple,((y)->:($x.$y)).(getfieldnames(fields))...)
end

"""
throws an error is the fields contain overlapping symbols
"""
function assertcollisionfree(x,y)
    xf,yf = getfieldnames(x), getfieldnames(y)
    if any(x->(x in xf),yf) #faster then sets for expected field counts
        throw("Field defined in multiple locations")
    end
end

"""
returns a copy with replacement fields
"""
function replacefields(struct_,fields)
    out = deepcopy(struct_)
    out.args[3].args = fields
    out
end

"""
adds source module information to the type name
"""
fulltypename(x,modulePath,inhibit) = x #is a literal

function fulltypename(x::Union{Expr,Symbol},modulePath,inhibit)
    if x in inhibit
        return x
    end
    if iscontainerlike(x)
        fullargs = (fulltypename(y,modulePath,inhibit) for y in x.args)
        return Expr(x.head,fullargs...)
    end
    qualifyname(x,modulePath)
end

function qualifyname(@nospecialize(x::Type))
    qualifyname(nameof(x),Symbol[fullname(parentmodule(x))...])
end

function qualifyname(x::Expr,modulePath)
    oldpath,x = getpath(x)

    isWrapped = length(modulePath) == 1
    out = Expr(:.,modulePath[1],QuoteNode(isWrapped ? oldpath[1] : modulePath[2]))
    for i = 3:length(modulePath)
        out = Expr(:.,out,QuoteNode(modulePath[i]))
    end
    for i = length(oldpath):-1:(isWrapped + 1)
        out = Expr(:.,out,QuoteNode(oldpath[i]))
    end
    Expr(:.,out,QuoteNode(x))
end
function qualifyname(x::Symbol,modulePath)
    if length(modulePath) == 1
        return Expr(:.,modulePath[1],QuoteNode(x))
    end
    out = Expr(:.,modulePath[1],QuoteNode(modulePath[2]))
    for i = 3:length(modulePath)
        out = Expr(:.,out,QuoteNode(modulePath[i]))
    end
    Expr(:.,out,QuoteNode(x))
end

"""
annotates module information to unanotated typed fields
"""
function sanitize!(modulePath,fields,inhibit)
    for i = eachindex(fields)
        temp = fields[i]
        if temp isa Expr
            if temp.args[2] isa FieldType
                temp.args[2] = fulltypename(temp.args[2]::FieldType,modulePath,inhibit)
            end
        end
    end
end

"""
update parameters from old fields
"""
function updateParameters(oldFields,oldParams,parameters,parentType,modulePath)
    update(x) = x
    function update(x::Symbol)
        if x in oldParams
            loc = findfirst(y->(y==x),oldParams)
            newParam = parameters[2][loc]
            if newParam in parameters[1]
                return newParam
            else
                return fulltypename(newParam,modulePath,parameters[1])
            end
        end
        x
    end
    function update(x::Expr)
        Expr(x.head,(update(z) for z in x.args)...)
    end
    update.(oldFields)
end

"""
turns an object into a tuple of its fields
"""
function totuple(x) #low efficiency version
    if fieldcount(typeof(x)) > 0
        tuple((getfield(x,y) for y in fieldnames(typeof(x)))...)
    else
        (x,)
    end
end

"""
returns a renamed struct
"""
function rename(struct_,name)
    newStruct = deepcopy(struct_)
    newStruct.args[2] = name
    newStruct
end

"""
strips parameterization off of a name that does
not include inheritence information
"""
deparametrize(name) = name
deparametrize(name::Expr) = name.head == :curly ? name.args[1] : name

"""
registers a new struct and abstract type pair
"""
function register(@nospecialize(concrete),@nospecialize(proto),fields,parameters,mutability)

    fieldBacking[proto] = fields
    shadowMap[concrete] = proto
    parameterMap[proto] = parameters
    mutabilityMap[proto] = mutability
    shadowMap[proto] = proto
end

inherits(x::Symbol) = false
inherits(x::Expr) = x.head != :curly

"""
@protostruct(struct_ [, prefix_])

creates a struct that can have structure inherited from it and can inherit
structure.

additionally it creates an abstract type with a name given by the struct
definitions name and a prefix. The concrete type inherits from the abstract
type and anything which inherits the concrete types structure also inherits
behavior from the abstract type.

```Julia
julia> using StructuralInheritance

julia> @protostruct struct A{T}
           fieldFromA::T
       end
ProtoA

julia> @protostruct struct B{D} <: A{Complex{D}}
          fieldFromB::D
       end "SomeOtherPrefix"
SomeOtherPrefixB

julia> @protostruct struct C <: B{Int}
         fieldFromC
       end
ProtoC
```
"""
macro protostruct(struct_,prefix_ = "Proto",mutablilityOverride = false)
    protostruct(__module__,struct_,prefix_,mutablilityOverride)
end

function protostruct(__module__,struct_,prefix_,mutablilityOverride)
    try
        prefix = (prefix_ isa Union{String,Symbol}) ? string(prefix_) : string(__module__.eval(prefix_))

        if length(prefix) == 0
            throw("Prefix must have finite Length")
        end

        struct_ = macroexpand(__module__,struct_,recursive=true)::Expr

        mutability = struct_.args[1]::Bool
        newName,name,newStructLightName,lightname = newnames(struct_.args[2]::FieldType,__module__,prefix)

        newParameters,oldParameters = get2parameters(struct_.args[2]::FieldType)

        fields = extractfields(struct_)
        prototypeDefinition = abstracttype(name)
        structDefinition = rename(struct_,newName)
        SI = StructuralInheritance
        modulePath = Symbol[fullname(__module__)...]
        if !inherits(name)
            sanitize!(modulePath,fields,newParameters)
        else #inheritence case
            parentType = get(shadowMap,__module__.eval(deparametrize(name.args[2])),nothing)
            oldFields = get(fieldBacking,parentType ,FieldType[])
            oldMutability = get(mutabilityMap,parentType,nothing)

            if oldMutability != nothing && oldMutability != mutability
                if eval(mutablilityOverride) != true
                    throw("$(oldMutability ? "im" : "")mutable object"*
                    " inheriting from $(oldMutability ? "" : "im")mutable")
                end
            end

            assertcollisionfree(fields,oldFields)
            sanitize!(modulePath,fields,newParameters)
            oldFields = updateParameters(oldFields,
                                        get(parameterMap,parentType,FieldType[]),
                                         (newParameters,oldParameters),
                                         parentType,
                                         modulePath)
            fields = vcat(oldFields,fields)
            constructors = extractconstructors(struct_)
            structDefinition = replacefields(structDefinition,
                                             vcat(fields,constructors))
        end
        return esc(quote
            $prototypeDefinition
            $structDefinition
            function $SI.totuple(x::$(deparametrize(newStructLightName)))
                $(tupleexpander(:x,fields))
            end
            $SI.register($(deparametrize(newStructLightName)),
                         $(deparametrize(lightname)),
                         $(Meta.quot(fields)),
                         $(Meta.quot(newParameters)),
                         $mutability)
        end)

    catch e
        return :(throw($e))
    end
end

end

max_stars_repo_name: WschW/StrucralInheritence

max_issues_repo_name: WschW/StrucralInheritence

max_forks_repo_name: WschW/StrucralInheritence


text:
export rule

@rule KernelGCV(:z, Marginalisation) (q_y_x::MvNormalMeanCovariance, meta::KernelGCVMetadata) = begin
    dims = Int64(ndims(q_y_x) / 2)

    m_yx   = mean(q_y_x)
    cov_yx = cov(q_y_x)

    cov11 = @view cov_yx[1:dims,1:dims]
    cov12 = @view cov_yx[1:dims,dims+1:end]
    cov21 = @view cov_yx[dims+1:end,1:dims]
    cov22 = @view cov_yx[dims+1:end,dims+1:end]

    m1 = @view m_yx[1:dims]
    m2 = @view m_yx[dims+1:end]

    psi = cov11 + cov22 - cov12 - cov21 + (m1 - m2)*(m1 - m2)'

    kernelfunction = get_kernelfn(meta)

    logpdf = (z) -> begin
        gz = kernelfunction(z)
        return -0.5*(logdet(gz) + tr(cholinv(gz)*psi))
    end

    return FnWithApproximation(logpdf, get_approximation(meta))
end

@rule KernelGCV(:z, Marginalisation) (q_y_x::MvNormalMeanPrecision, meta::KernelGCVMetadata) = begin
    dims = Int64(ndims(q_y_x) / 2)

    m_yx   = mean(q_y_x)
    cov_yx = cov(q_y_x)

    cov11 = @view cov_yx[1:dims,1:dims]
    cov12 = @view cov_yx[1:dims,dims+1:end]
    cov21 = @view cov_yx[dims+1:end,1:dims]
    cov22 = @view cov_yx[dims+1:end,dims+1:end]

    m1 = @view m_yx[1:dims]
    m2 = @view m_yx[dims+1:end]

    psi = cov11 + cov22 - cov12 - cov21 + (m1 - m2)*(m1 - m2)'

    kernelfunction = get_kernelfn(meta)

    logpdf = (z) -> begin
        gz = kernelfunction(z)
        return -0.5*(logdet(gz) + tr(cholinv(gz)*psi))
    end

    return FnWithApproximation(logpdf, get_approximation(meta))
end


text:
export assembleStiffAndForce,assembleInternalForce,assembleMassMatrix!,tfAssembleInternalForce

@doc """
    Numerically assemble internal force vector, compute local internal force f_int
    from domain.state and then assemble to F_int
    - 'globdat': GlobalData
    - 'domain': Domain, finite element domain, for data structure
    - 'Δt':  Float64, current time step size
    Return F_int: Float64[neqs], internal force vector
"""->function assembleInternalForce(globdat::GlobalData, domain::Domain, Δt::Float64 = 0.0)
    Fint = zeros(Float64, domain.neqs)
    neles = domain.neles

    # Loop over the elements in the elementGroup
    for iele  = 1:neles
        element = domain.elements[iele]

      # Get the element nodes
        el_nodes = getNodes(element)

      # Get the element nodes
        el_eqns = getEqns(domain, iele)

        el_dofs = getDofs(domain, iele)

        el_state  = getState(domain, el_dofs)

        el_Dstate = getDstate(domain, el_dofs)

      # Get the element contribution by calling the specified action
        fint = getInternalForce(element, el_state, el_Dstate, Δt)

      # Assemble in the global array
        el_eqns_active = (el_eqns .>= 1)
        Fint[el_eqns[el_eqns_active]] += fint[el_eqns_active]
    end

    return Fint
end

@doc """
    Tensorflow version of assembling internal force vector, compute local internal
    force f_int and then assemble to F_int, which generates inverse problem automatically.

    - 'domain': Domain, finite element domain, for data structure
    - 'nn': Function strain -> stress, neural network constitutive law function
    - 'E_all': PyObject(Float64)[neles*nGauss, nstrains], neles*nGauss is the number of Gaussian quadrature points,
               nstrain is the number of strain components. All strains for the current time-step
    - 'DE_all': PyObject(Float64)[neles*nGauss, nstrains], neles*nGauss is the number of Gaussian quadrature points,
                nstrain is the number of strain components. All strains for the previous time-step
    - 'w∂E∂u_all': PyObject(Float64)[neles*nGauss, ndofs_per_element, nstrains], neles*nGauss is the number of Gaussian quadrature points,
                   ndofs_per_element is the number of freedoms per element, nstrain is the number of strain components.
                   Multiplication of the Gaussian weight and ∂E∂u^T for current time-step,

    - 'σ0_all': PyObject(Float64)[neles*nGauss, nstrains], neles*nGauss is the number of Gaussian quadrature points,
    nstrain is the number of strain components.  All stresses for the last time-step

    Return: internal force vector F_int, PyObject(Float64)[neqns] and the predicted stresses at the
            current step σ_all, PyObject(Float64)[neles*nGauss, nstrains]
"""->function tfAssembleInternalForce(domain::Domain, nn::Function, E_all::PyObject, DE_all::PyObject, w∂E∂u_all::PyObject, σ0_all::PyObject)
    neles = domain.neles
    nGauss = length(domain.elements[1].weights)
    neqns_per_elem = length(getEqns(domain, 1))
    nstrains = size(E_all, 2)

    @assert size(E_all) == (neles * nGauss, nstrains)
    @assert size(DE_all) == (neles * nGauss, nstrains)
    @assert size(σ0_all) == (neles * nGauss, nstrains)

  # el_eqns_all, the equation numbers related to the Gaussian points in the element, negative value means Drichlet boundary
    el_eqns_all = zeros(Int32, neles * nGauss, neqns_per_elem)
  # el_eqns_active_all = el_eqns_all > 0
    el_eqns_active_all = zeros(Bool, neles * nGauss, neqns_per_elem)

  # Loop over the elements in the elementGroup to construct el_eqns_active_all , el_eqns_all and w∂E∂u_all
    for iele  = 1:neles
    # Get the element nodes
        el_eqns = getEqns(domain, iele)

    # Assemble in the global array
        el_eqns_active_all[(iele - 1) * nGauss + 1:iele * nGauss,:] = repeat((el_eqns .>= 1)', nGauss, 1)
        el_eqns_all[(iele - 1) * nGauss + 1:iele * nGauss,:] = repeat(el_eqns', nGauss, 1)
    end

  # get stress at each Gaussian points
    σ_all = nn(E_all, DE_all, σ0_all)

  # compute fint at each Gaussian quadrature points, fints[igp] = w∂E∂u_all[igp] * σ_all[igp]
    fints = squeeze(tf.matmul(w∂E∂u_all, tf.expand_dims(σ_all, 2)))

  # call the Cpp assembler to construct Fint
    Fint = cpp_fint(fints, constant(el_eqns_all, dtype = Int32), constant(domain.neqs, dtype = Int32))

    return Fint, σ_all
end

# function tfAssembleInternalForce(domain::Domain, nn::Function, E_all::PyObject, DE_all::PyObject,
#      w∂E∂u_all::PyObject, σ0_all::PyObject, α::PyObject)
#   neles = domain.neles
#   nGauss = length(domain.elements[1].weights)
#   neqns_per_elem = length(getEqns(domain,1))
#   nstrains = size(E_all,2) #todo

#   @assert size(E_all)==(neles*nGauss, nstrains)
#   @assert size(DE_all)==(neles*nGauss, nstrains)
#   @assert size(σ0_all)==(neles*nGauss, nstrains)

#   # el_eqns_all, the equation numbers related to the Gaussian point, negative value means Drichlet boundary
#   el_eqns_all = zeros(Int32, neles*nGauss, neqns_per_elem)
#   # el_eqns_active_all = el_eqns_all > 0
#   el_eqns_active_all = zeros(Bool, neles*nGauss, neqns_per_elem)
#   # Loop over the elements in the elementGroup to construct el_eqns_active_all , el_eqns_all and w∂E∂u_all
#   for iele  = 1:neles
#     # Get the element nodes
#     el_eqns = getEqns(domain,iele)

#     # Assemble in the global array
#     el_eqns_active_all[(iele-1)*nGauss+1:iele*nGauss,:] = repeat((el_eqns .>= 1)', nGauss, 1)
#     el_eqns_all[(iele-1)*nGauss+1:iele*nGauss,:] = repeat(el_eqns', nGauss, 1)
#   end

#   # get stress at each Gaussian points
#   # @info "* ", E_all, DE_all, σ0_all
#   σ_all, α = nn(E_all, DE_all, σ0_all, α) # α is the internal variable
#   fints = squeeze(tf.matmul(w∂E∂u_all, tf.expand_dims(σ_all,2)))
#   Fint = cpp_fint(fints,constant(el_eqns_all, dtype=Int32),constant(domain.neqs, dtype=Int32))
#   return Fint, σ_all, α
# end

@doc """
    Numerically assemble internal force vector and stiffness matrix, compute local internal force f_int/
    stiffness matrix Slocal from domain.state and then assemble to F_int/Ksparse
    - 'globdat': GlobalData
    - 'domain': Domain, finite element domain, for data structure
    - 'Δt':  Float64, current time step size
    Return F_int: Float64[neqs], internal force vector
    Return Ksparse: Float64[neqs, neqns], sparse stiffness matrix
"""->function assembleStiffAndForce(globdat::GlobalData, domain::Domain, Δt::Float64 = 0.0)
    Fint = zeros(Float64, domain.neqs)
  # K = zeros(Float64, domain.neqs, domain.neqs)
    ii = Int64[]; jj = Int64[]; vv = Float64[]
    neles = domain.neles

  # Loop over the elements in the elementGroup
    for iele  = 1:neles
        element = domain.elements[iele]

    # Get the element nodes
        el_nodes = getNodes(element)

    # Get the element nodes equation numbers
        el_eqns = getEqns(domain, iele)

    # Get the element nodes dof numbers
        el_dofs = getDofs(domain, iele)

        el_state  = getState(domain, el_dofs)

        el_Dstate = getDstate(domain, el_dofs)

    # Get the element contribution by calling the specified action
        fint, stiff  = getStiffAndForce(element, el_state, el_Dstate, Δt)

    # Assemble in the global array
        el_eqns_active = el_eqns .>= 1
    # K[el_eqns[el_eqns_active], el_eqns[el_eqns_active]] += stiff[el_eqns_active,el_eqns_active]
        Slocal = stiff[el_eqns_active,el_eqns_active]
        Idx = el_eqns[el_eqns_active]
        for i = 1:length(Idx)
            for j = 1:length(Idx)
                push!(ii, Idx[i])
                push!(jj, Idx[j])
                push!(vv, Slocal[i,j])
            end
        end
        Fint[el_eqns[el_eqns_active]] += fint[el_eqns_active]
    end
    Ksparse = sparse(ii, jj, vv, domain.neqs, domain.neqs)

    return Fint, Ksparse
end

# function assembleStiffAndForce(globdat::GlobalData, domain::Domain, Δt::Float64 = 0.0)
#   Fint = zeros(Float64, domain.neqs)
#   K = zeros(Float64, domain.neqs, domain.neqs)
#   neles = domain.neles

#   # Loop over the elements in the elementGroup
#   for iele  = 1:neles
#     element = domain.elements[iele]

#     # Get the element nodes
#     el_nodes = getNodes(element)

#     # Get the element nodes
#     el_eqns = getEqns(domain,iele)

#     el_dofs = getDofs(domain,iele)

#     #@show "iele", iele, el_dofs

#     #@show "domain.state", iele, domain.state

#     el_state  = getState(domain,el_dofs)

#     el_Dstate = getDstate(domain,el_dofs)
#     # #@show "+++++", el_state, el_Dstate

#     # Get the element contribution by calling the specified action
#     #@info "ele id is ", iele
#     fint, stiff  = getStiffAndForce(element, el_state, el_Dstate, Δt)

#     # Assemble in the global array
#     el_eqns_active = el_eqns .>= 1
#     K[el_eqns[el_eqns_active], el_eqns[el_eqns_active]] += stiff[el_eqns_active,el_eqns_active]
#     Fint[el_eqns[el_eqns_active]] += fint[el_eqns_active]

#     # @info "Fint is ", Fint
#   end
#   return Fint, sparse(K)
# end

@doc """
    compute constant sparse mass matrix
    due to the time-dependent Dirichlet boundary condition
    mass matrix = M,    MID
                  MID'  MDD

    - 'globdat': GlobalData
    - 'domain': Domain, finite element domain, for data structure
    here M is Float64[neqns, neqns]
         MID is Float64[neqns, nd1]
         Mlumped is Float64[neqns]

    update M and MID and Mlumpe in globaldat
"""->function assembleMassMatrix!(globaldat::GlobalData, domain::Domain)
    Mlumped = zeros(Float64, domain.neqs)
    # M = zeros(Float64, domain.neqs, domain.neqs)
    iiM = Int64[]; jjM = Int64[]; vvM = Float64[]
    Mlumped = zeros(Float64, domain.neqs)
    neles = domain.neles

    nnodes, ndims = domain.nnodes, domain.ndims

    # MID = zeros(domain.neqs, sum(domain.EBC .== -2))
    iiMID = Int64[]; jjMID = Int64[]; vvMID = Float64[]

    # construct map from freedoms(first direction, second direction ...) to time-dependent Dirichlet freedoms
    dofs_to_EBCdofs = zeros(Int64, nnodes*ndims)
    if globaldat.EBC_func != nothing
        dof_id = 0
        for idof = 1:ndims
            for inode = 1:nnodes
                if (domain.EBC[inode, idof] == -2)
                    dof_id += 1
                    dofs_to_EBCdofs[inode + (idof - 1) * nnodes] = dof_id
                end
            end
        end
    end

    # Loop over the elements in the elementGroup
    for iele = 1:neles
        element = domain.elements[iele]

        # Get the element nodes
        el_nodes = getNodes(element)

        # Get the element nodes equation numbers
        el_eqns = getEqns(domain, iele)

        # Get the element nodes dof numbers
        el_dofs = getDofs(domain, iele)

        # Get the element contribution by calling the specified action
        lM, lumped = getMassMatrix(element)

        # Assemble in the global array

        el_eqns_active = (el_eqns .>= 1)

        # time-dependent Dirichlet boundary condition
        el_eqns_acc_active = (el_eqns .== -2)

        # M[el_eqns[el_eqns_active], el_eqns[el_eqns_active]] += lM[el_eqns_active, el_eqns_active]
        Idx = el_eqns[el_eqns_active]
        Mlocal = lM[el_eqns_active, el_eqns_active]
        for i = 1:length(Idx)
            for j = 1:length(Idx)
                push!(iiM, Idx[i])
                push!(jjM, Idx[j])
                push!(vvM, Mlocal[i,j])
            end
        end

        Mlumped[el_eqns[el_eqns_active]] += lumped[el_eqns_active]

        if globaldat.EBC_func != nothing
            Idx = el_eqns[el_eqns_active]
            Idy = dofs_to_EBCdofs[el_dofs[el_eqns_acc_active]]

            Mlocal = lM[el_eqns_active, el_eqns_acc_active]

            for i = 1:length(Idx)
                for j = 1:length(Idy)
                    push!(iiMID, Idx[i])
                    push!(jjMID, Idy[j])
                    push!(vvMID, Mlocal[i,j])
                end
            end
        end

    end

    globaldat.M = sparse(iiM, jjM, vvM, domain.neqs, domain.neqs)
    globaldat.Mlumped = Mlumped
    globaldat.MID = sparse(iiMID, jjMID, vvMID, domain.neqs, sum(domain.EBC .== -2))

end

# @doc """
#     compute constant mass matrix as sparse matrix
#     due to the time-dependent Dirichlet boundary condition
#     mass matrix = M,    MID
#                   MID'  MDD
#     save M and MID and lump(M)
# """->

# function assembleMassMatrix!(globaldat::GlobalData, domain::Domain)
#   Mlumped = zeros(Float64, domain.neqs)
#   M = zeros(Float64, domain.neqs, domain.neqs)
#   Mlumped = zeros(Float64, domain.neqs)
#   neles = domain.neles

#   MID = zeros(domain.neqs, sum(domain.EBC .== -2))
#   # Loop over the elements in the elementGroup
#   for iele = 1:neles
#       element = domain.elements[iele]

#       # Get the element nodes
#       el_nodes = getNodes(element)

#       # Get the element nodes
#       el_eqns = getEqns(domain,iele)

#       # Get the element contribution by calling the specified action
#       lM, lumped = getMassMatrix(element)

#       # Assemble in the global array

#       el_eqns_active = (el_eqns .>= 1)
#       el_eqns_acc_active = (el_eqns .== -2)

#       M[el_eqns[el_eqns_active], el_eqns[el_eqns_active]] += lM[el_eqns_active, el_eqns_active]
#       Mlumped[el_eqns[el_eqns_active]] += lumped[el_eqns_active]
#       MID[el_eqns[el_eqns_active], el_eqns[el_eqns_acc_active]] += lM[el_eqns_active, el_eqns_acc_active]

#   end

#   globaldat.M = sparse(M)
#   globaldat.Mlumped = sparse(Mlumped)
#   globaldat.MID = MID
# end


text:
"""
    load(fpath::String; gamma::Union{Nothing,Float64}=nothing, expand_paletted::Bool=false)
    load(s::IO; gamma::Union{Nothing,Float64}=nothing, expand_paletted::Bool=false)

Read a `.png` image from file at `fpath` or in IO stream `s`.
`gamma` can be used to override the automatic gamma correction, a value of 1.0
means no gamma correction.

The result will be an 8 bit (N0f8) image if the source bit depth is <= 8 bits, 16 bit (N0f16)
otherwise.

The number of channels (and transparency) of the source determines the color type of the output:
- 1 channel  -> Gray
- 2 channels -> GrayA
- 3 channels -> RGB
- 4 channels -> RGBA

When reading in simple paletted images, i.e. having a PLTE chunk and an 8 bit depth, the image will
be represented as an `IndirectArray` with `OffsetArray` `values` field. To always get back a plain
`Matrix` of colorants, use `expand_paletted=true`.
"""
function load(fpath::String; gamma::Union{Nothing,Float64}=nothing, expand_paletted::Bool=false)
    fp = open_png(fpath)
    png_ptr = create_read_struct()
    @debug "Load PNG File:" fpath png_ptr
    info_ptr = create_info_struct(png_ptr)
    png_init_io(png_ptr, fp)
    png_set_sig_bytes(png_ptr, PNG_BYTES_TO_CHECK)
    out = _load(png_ptr, info_ptr, gamma=gamma, expand_paletted=expand_paletted)
    close_png(fp)
    return out
end

maybe_lock(f, io::IO) = lock(f, io)
# IOStream doesn't support locking...
maybe_lock(f, io::IOStream) = f()
maybe_lock(f, io::IOBuffer) = f()
maybe_lock(f, io::Base64EncodePipe) = f()
maybe_lock(f, io::IOContext) = maybe_lock(f, io.io)

if VERSION < v"1.6.0-DEV.1652"
    _iswritable(x) = Base.iswritable(x)
    _iswritable(pipe::Base64EncodePipe) = Base.iswritable(pipe.io)
else
    _iswritable(x) = iswritable(x)
end

function load(s::IO; gamma::Union{Nothing,Float64}=nothing, expand_paletted::Bool=false)
    isreadable(s) || throw(ArgumentError("read failed, IOStream is not readable"))
    Base.eof(s) && throw(EOFError())

    png_ptr = create_read_struct()
    info_ptr = create_info_struct(png_ptr)

    maybe_lock(s) do
        if s isa IOBuffer
            png_set_read_fn(png_ptr, pointer_from_objref(s), readcallback_iobuffer_c[])
        else
            png_set_read_fn(png_ptr, s.handle, readcallback_c[])
        end
        # https://stackoverflow.com/questions/22564718/libpng-error-png-unsigned-integer-out-of-range
        png_set_sig_bytes(png_ptr, 0)
        return _load(png_ptr, info_ptr, gamma=gamma, expand_paletted=expand_paletted)
    end
end

function _readcallback(png_ptr::png_structp, data::png_bytep, length::png_size_t)::Cvoid
    a = png_get_io_ptr(png_ptr)
    ccall(:ios_readall, Csize_t, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), a, data, length)
    return
end

function _readcallback_iobuffer(png_ptr::png_structp, data::png_bytep, length::png_size_t)::Cvoid
    a = png_get_io_ptr(png_ptr)
    io = unsafe_pointer_to_objref(a)
    unsafe_read(io, Ptr{UInt8}(data), length)
    return
end

function _load(png_ptr, info_ptr; gamma::Union{Nothing,Float64}=nothing, expand_paletted::Bool=false)
    png_read_info(png_ptr, info_ptr)
    flags = png_get_valid(
        png_ptr, info_ptr,
        PNG_INFO_sRGB | PNG_INFO_gAMA | PNG_INFO_tRNS | PNG_INFO_cHRM | PNG_INFO_bKGD | PNG_INFO_PLTE
    )
    valid_sRGB = (flags & PNG_INFO_sRGB) != 0
    valid_gAMA = (flags & PNG_INFO_gAMA) != 0
    valid_tRNS = (flags & PNG_INFO_tRNS) != 0
    valid_cHRM = (flags & PNG_INFO_cHRM) != 0
    valid_bKGD = (flags & PNG_INFO_bKGD) != 0
    valid_PLTE = (flags & PNG_INFO_PLTE) != 0

    width = png_get_image_width(png_ptr, info_ptr)
    height = png_get_image_height(png_ptr, info_ptr)
    color_type = color_type_orig = png_get_color_type(png_ptr, info_ptr)
    bit_depth = bit_depth_orig = png_get_bit_depth(png_ptr, info_ptr)
    num_channels = png_get_channels(png_ptr, info_ptr)
    interlace_type = png_get_interlace_type(png_ptr, info_ptr)

    if valid_bKGD
        backgroundp = png_color_16p()
        if png_get_bKGD(png_ptr, info_ptr, Ptr{png_color_16p}(backgroundp)) != 0
            png_set_background(png_ptr, backgroundp, PNG_BACKGROUND_GAMMA_FILE, 1, 1.0)
        end
    end

    screen_gamma = PNG_DEFAULT_sRGB
    image_gamma = Ref{Cdouble}(-1.0)
    intent = Ref{Cint}(-1)
    if isnothing(gamma)
        if valid_sRGB
            if png_get_sRGB(png_ptr, info_ptr, intent) != 0
                png_set_gamma(png_ptr, screen_gamma, PNG_DEFAULT_sRGB)
            else
                if png_get_gAMA(png_ptr, info_ptr, image_gamma) == 0
                    image_gamma[] = 0.45455
                end
                png_set_gamma(png_ptr, screen_gamma, image_gamma[])
            end
        elseif valid_gAMA
            if png_get_gAMA(png_ptr, info_ptr, image_gamma) == 0
                image_gamma[] = 0.45455
            end
            png_set_gamma(png_ptr, screen_gamma, image_gamma[])
        end
    elseif gamma != 1
        image_gamma[] = gamma
        png_set_gamma(png_ptr, screen_gamma, image_gamma[])
    end

    read_as_paletted = !expand_paletted &&
        color_type == PNG_COLOR_TYPE_PALETTE &&
        bit_depth == 8 &&
        valid_PLTE

    if !read_as_paletted
        if color_type == PNG_COLOR_TYPE_PALETTE
            png_set_palette_to_rgb(png_ptr)
            color_type = PNG_COLOR_TYPE_RGB
        end

        if color_type == PNG_COLOR_TYPE_GRAY && bit_depth < 8
            png_set_expand_gray_1_2_4_to_8(png_ptr)
            png_set_packing(png_ptr)
            bit_depth = 8
        end

        if valid_tRNS
            png_set_tRNS_to_alpha(png_ptr)
            if color_type == PNG_COLOR_TYPE_GRAY || color_type == PNG_COLOR_TYPE_RGB
                color_type |= PNG_COLOR_MASK_ALPHA
            end
        end
        buffer_eltype = _buffer_color_type(color_type, bit_depth)
        bit_depth == 16 && png_set_swap(png_ptr)
    end
    n_passes = png_set_interlace_handling(png_ptr)
    png_read_update_info(png_ptr, info_ptr)

    # Gamma correction is applied to a palette after `png_read_update_info` is called
    if read_as_paletted
        palette_length = Ref{Cint}()
        # TODO: Figure out the lenght of paletted before calling png_get_PLTEs
        palette_buffer = Vector{RGB{N0f8}}(undef, PNG_MAX_PALETTE_LENGTH)
        png_get_PLTE(png_ptr, info_ptr, pointer_from_objref(palette_buffer), palette_length)
        palette = palette_buffer[1:palette_length[]]

        if valid_tRNS
            alpha_buffer = Vector{_AlphaBuffer}(undef, palette_length[])
            alphas_cnt = Ref{Cint}()
            png_get_tRNS(png_ptr, info_ptr, pointer_from_objref(alpha_buffer), alphas_cnt, C_NULL)
            if alphas_cnt[] > 1
                palette = map(x->RGBA(x[1], x[2].val), zip(palette, alpha_buffer))
            else
                # Seems that if there is only one transparency entry, it is applied to the first
                # color in palette while the rest of the colors are opaque
                α = alpha_buffer[1].val
                palette = map(x->RGBA(x[2], x[1] == 1 ? α : 1), enumerate(palette))
            end
        end
        buffer_eltype = UInt8
    end

    # We transpose to work around libpng expecting row-major arrays
    buffer = Array{buffer_eltype}(undef, width, height)

    @debug(
        "Read PNG info:",
        png_ptr,
        height,
        width,
        color_type_orig,
        color_type,
        bit_depth_orig,
        bit_depth,
        num_channels,
        interlace_type,
        gamma,
        image_gamma[],
        screen_gamma,
        intent[],
        read_as_paletted,
        n_passes,
        buffer_eltype,
        valid_sRGB,
        valid_gAMA,
        valid_tRNS,
        valid_cHRM,
        valid_bKGD,
        valid_PLTE,
        PNG_HEADER_VERSION_STRING
    )

    buffer = Base.invokelatest(_load!, buffer, png_ptr, info_ptr)
    if expand_paletted || color_type != PNG_COLOR_TYPE_PALETTE
        return buffer
    else
        # We got 0-based indices back from libpng and converting to 1-based could overflow UInt8.
        # Using UInt16 for index would cost us large part of the savings provided by IndirectArray.
        return IndirectArray(buffer, OffsetArray(palette, -1))
    end
end

function _load!(buffer::Matrix{T}, png_ptr, info_ptr) where T    # separate to support precompilation of permutedims
    png_read_image(png_ptr, map(pointer, eachcol(buffer)))
    png_read_end(png_ptr, info_ptr)
    png_destroy_read_struct(Ref{Ptr{Cvoid}}(png_ptr), Ref{Ptr{Cvoid}}(info_ptr), C_NULL)
    return permutedims(buffer, (2, 1))
end

function _buffer_color_type(color_type, bit_depth)
    bit_depth = Int(bit_depth)
    if color_type == PNG_COLOR_TYPE_GRAY
        colors_type = Gray{bit_depth > 8 ? Normed{UInt16,bit_depth} : Normed{UInt8,bit_depth}}
    elseif color_type == PNG_COLOR_TYPE_PALETTE
        colors_type = RGB{bit_depth == 16 ? N0f16 : N0f8}
    elseif color_type == PNG_COLOR_TYPE_RGB
        colors_type = RGB{bit_depth == 16 ? N0f16 : N0f8}
    elseif color_type == PNG_COLOR_TYPE_RGB_ALPHA
        colors_type = RGBA{bit_depth == 16 ? N0f16 : N0f8}
    elseif color_type == PNG_COLOR_TYPE_GRAY_ALPHA
        colors_type = GrayA{bit_depth > 8 ? Normed{UInt16,bit_depth} : Normed{UInt8,bit_depth}}
    else
        throw(error("Unknown color type: $color_type"))
    end
    return colors_type
end

### Write ##########################################################################################

const SupportedPaletteColor = Union{
    AbstractRGB{<:Union{N0f8,AbstractFloat}},
    TransparentRGB{T,<:Union{N0f8,AbstractFloat}} where T,
}

"""
    save(fpath::String, image::AbstractArray;
         compression_level::Integer=0, compression_strategy::Integer=3, filters::Integer=4)
    save(s::IO, image::AbstractArray;
         compression_level::Integer=0, compression_strategy::Integer=3, filters::Integer=4)

Writes `image` as a png to file at `fpath`, or to IO stream `s`.

## Arguments
- `compression_level`: 0 (`Z_NO_COMPRESSION`), 1 (`Z_BEST_SPEED`), ..., 9 (`Z_BEST_COMPRESSION`)
- `compression_strategy`: 0 (`Z_DEFAULT_STRATEGY`), 1 (`Z_FILTERED`), 2 (`Z_HUFFMAN_ONLY`), 3 (`Z_RLE`), 4 (`Z_FIXED`)
- `filters`: 0 (None), 1 (Sub), 2 (Up), 3 (Average), 4 (Paeth).

The saved image will have 16 bits of depth if the `image` has eltype that is based on `UInt16`,
8 bits otherwise.

The number of channels and element type of the `image` determines the color type of the
output:
- 0/1 channel Float / Integer / Normed or Gray eltype        -> `PNG_COLOR_TYPE_GRAY`
- 2 channels Float  / Integer / Normed or GrayA eltype       -> `PNG_COLOR_TYPE_GRAY_ALPHA`
- 3 channels Float  / Integer / Normed or RGB / BGR eltype   -> `PNG_COLOR_TYPE_RGB`
- 4 channels Float  / Integer / Normed or ARGB / ABGR eltype -> `PNG_COLOR_TYPE_RGB_ALPHA`

When `image` is an `IndirectArray` with up to 256 unique RGB colors, the result is encoded as a
"paletted image".

"""
function save(
    fpath::String,
    image::S;
    compression_level::Integer = Z_BEST_SPEED,
    compression_strategy::Integer = Z_RLE,
    filters::Integer = Int(PNG_FILTER_PAETH),
    palette::Union{Nothing,AbstractVector{<:Union{RGB{N0f8},RGBA{N0f8}}}} = nothing,
) where {
    T,
    S<:Union{AbstractMatrix,AbstractArray{T,3}}
}
    @assert Z_DEFAULT_STRATEGY <= compression_strategy <= Z_FIXED
    @assert Z_NO_COMPRESSION <= compression_level <= Z_BEST_COMPRESSION
    @assert 2 <= ndims(image) <= 3
    @assert size(image, 3) <= 4

    fp = ccall(:fopen, Ptr{Cvoid}, (Cstring, Cstring), fpath, "wb")
    fp == C_NULL && error("Could not open $(fpath) for writing")

    png_ptr = create_write_struct()
    @debug "Save PNG File:" fpath png_ptr
    info_ptr = create_info_struct(png_ptr)

    png_init_io(png_ptr, fp)

    _save(png_ptr, info_ptr, image,
        compression_level=compression_level,
        compression_strategy=compression_strategy,
        filters=filters,
        palette=palette)

    close_png(fp)
end
function save(
    s::IO,
    image::S;
    compression_level::Integer = Z_BEST_SPEED,
    compression_strategy::Integer = Z_RLE,
    filters::Integer = Int(PNG_FILTER_PAETH),
    palette::Union{Nothing,AbstractVector{<:Union{RGB{N0f8},RGBA{N0f8}}}} = nothing,
) where {
    T,
    S<:Union{AbstractMatrix,AbstractArray{T,3}}
}
    @assert Z_DEFAULT_STRATEGY <= compression_strategy <= Z_FIXED
    @assert Z_NO_COMPRESSION <= compression_level <= Z_BEST_COMPRESSION
    @assert 2 <= ndims(image) <= 3
    @assert size(image, 3) <= 4
    _iswritable(s) || throw(ArgumentError("write failed, IOStream is not writeable"))

    png_ptr = create_write_struct()
    info_ptr = create_info_struct(png_ptr)
    maybe_lock(s) do
        r = Ref{Any}(s)
        GC.@preserve r begin
            png_set_write_fn(png_ptr, r, writecallback_c[], C_NULL)

            _save(png_ptr, info_ptr, image,
                compression_level=compression_level,
                compression_strategy=compression_strategy,
                filters=filters,
                palette=palette)
        end
    end
end

function _save(png_ptr, info_ptr, image::S;
    compression_level::Integer = Z_BEST_SPEED,
    compression_strategy::Integer = Z_RLE,
    filters::Integer = Int(PNG_FILTER_PAETH),
    palette::Union{Nothing,AbstractVector{<:Union{RGB{N0f8},RGBA{N0f8}}}} = nothing,
) where {
    T,
    S<:Union{AbstractMatrix,AbstractArray{T,3}}
}
    image = _enforce_dense_cpu_array(image)
    height, width = size(image)[1:2]
    bit_depth = _get_bit_depth(image)
    color_type = _get_color_type(image)

    png_set_filter(png_ptr, PNG_FILTER_TYPE_BASE, UInt32(filters))
    png_set_compression_level(png_ptr, compression_level)
    png_set_compression_strategy(png_ptr, compression_strategy)

    if color_type == PNG_COLOR_TYPE_PALETTE
        # TODO: 1, 2, 4 bit-depth indices for palleted
        _png_check_paletted(image)
        palette = image.values
        color_count = length(palette)
        if eltype(palette) <: TransparentRGB
            png_set_PLTE(png_ptr, info_ptr, _standardize_palette(color.(palette)), color_count)
            png_set_tRNS(png_ptr, info_ptr, _palette_alpha(palette), color_count, C_NULL)
        else
            png_set_PLTE(png_ptr, info_ptr, _standardize_palette(palette), color_count)
        end
    else
        image_eltype = eltype(image)
        if (image_eltype <: BGR || image_eltype <: BGRA || image_eltype <: ABGR || image_eltype <: ARGB32)
            png_set_bgr(png_ptr)
        end

        if (image_eltype <: ABGR || image_eltype <: ARGB)
            png_set_swap_alpha(png_ptr)
        end

        if color_type == PNG_COLOR_TYPE_GRAY && bit_depth < 8
            png_set_packing(png_ptr)
            bit_depth = 8  # TODO: support 1, 2, 4 bit-depth gray images
        end
    end

    # gAMA and cHRM chunks should be always present for compatibility with older systems
    png_set_sRGB_gAMA_and_cHRM(png_ptr, info_ptr, PNG_sRGB_INTENT_PERCEPTUAL)

    @debug(
        "Write PNG info:",
        png_ptr,
        height,
        width,
        bit_depth,
        color_type,
        filters,
        compression_level,
        compression_strategy,
        palette,
        typeof(image)
    )

    # TODO: on error this throws an abort signal because we currently don't handle `longjmp`
    png_set_IHDR(
        png_ptr,
        info_ptr,
        width,
        height,
        bit_depth,
        color_type,
        PNG_INTERLACE_NONE,
        PNG_COMPRESSION_TYPE_BASE,
        PNG_FILTER_TYPE_BASE,
    )

    png_write_info(png_ptr, info_ptr)
    # Handles endianness for 16 bit, must be set after png_write_info
    bit_depth == 16 && png_set_swap(png_ptr)

    # We transpose to work around libpng expecting row-major arrays
    _write_image(permutedims(_prepare_buffer(image), (2, 1)), png_ptr, info_ptr)

    png_destroy_write_struct(Ref(png_ptr), Ref(info_ptr))
end

function _writecallback(png_ptr::png_structp, data::png_bytep, length::png_size_t)::Csize_t
    a = png_get_io_ptr(png_ptr)
    io = unsafe_load(Ptr{Any}(a))
    unsafe_write(io, data, length)
end

function _write_image(buf::AbstractArray{T,2}, png_ptr::Ptr{Cvoid}, info_ptr::Ptr{Cvoid}) where {T}
    ccall(
        (:png_write_image, libpng),
        Cvoid,
        (Ptr{Cvoid}, Ptr{Ptr{T}}),
        png_ptr,
        map(pointer, eachcol(buf)),
    )
    png_write_end(png_ptr, info_ptr)
end

function _png_check_paletted(image)
    palette = image.values
    color_count = length(palette)
    color_type = eltype(palette)

    ndims(image) != 2 && throw(ArgumentError("Only 2D `IndirectArrays` are supported"))
    color_count > 256 && throw(ArgumentError("Maximum size of `image.velues` is 256 colors"))
    if !(color_type <: SupportedPaletteColor)
        throw(ArgumentError(
            "Only 8-bit (transparent) RGB colors are supported for paletted images"
        ))
    end
end

_prepare_buffer(x::IndirectArray) = UInt8.(x.index .- first(axes(x.values, 1)))
_prepare_buffer(x::BitArray) = _prepare_buffer(collect(x))
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:Colorant{<:Normed}} = x
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:UInt8} =  reinterpret(Gray{N0f8}, x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:UInt16} = reinterpret(Gray{N0f16}, x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:Normed} = reinterpret(Gray{T}, x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:Gray{<:AbstractFloat}} =  Gray{N0f8}.(x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:GrayA{<:AbstractFloat}} = GrayA{N0f8}.(x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:RGB{<:AbstractFloat}} =   RGB{N0f8}.(x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:RGBA{<:AbstractFloat}} =  RGBA{N0f8}.(x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:BGR{<:AbstractFloat}} =   BGR{N0f8}.(x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:BGRA{<:AbstractFloat}} =  BGRA{N0f8}.(x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:ARGB{<:AbstractFloat}} =  ARGB{N0f8}.(x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:ABGR{<:AbstractFloat}} =  ABGR{N0f8}.(x)
_prepare_buffer(x::AbstractMatrix{<:T}) where {T<:Union{AbstractFloat,Bool}} = reinterpret(Gray{N0f8}, N0f8.(x))
_prepare_buffer(x::AbstractArray{T,3}) where {T<:Union{AbstractFloat,Bool}} = __prepare_buffer(N0f8.(x))
_prepare_buffer(x::AbstractArray{T,3}) where {T<:Union{UInt8,Int8}} = __prepare_buffer(reinterpret(N0f8, x))
_prepare_buffer(x::AbstractArray{T,3}) where {T<:Union{UInt16,Int16}} = __prepare_buffer(reinterpret(N0f16, x))
_prepare_buffer(x::AbstractArray{T,3}) where {T<:Normed} = __prepare_buffer(x)

function __prepare_buffer(x::AbstractArray{T,3}) where {T}
    nchannels = size(x, 3)
    if nchannels == 1
        ifelse(ndims(x) == 3, _prepare_buffer(dropdims(x, dims=3)), x)
    elseif nchannels == 2
        GrayA.(colorview(GrayA, view(x, :, :, 1), view(x, :, :, 2)))
    elseif nchannels == 3
        RGB.(colorview(RGB, view(x, :, :, 1), view(x, :, :, 2), view(x, :, :, 3)))
    elseif nchannels == 4
        RGBA.(colorview(RGBA, view(x, :, :, 1), view(x, :, :, 2), view(x, :, :, 3), view(x, :, :, 4)))
    else
        error("Unsupported number of channels $(nchannels) in input. Only <= 4 is expected.")
    end
end

# On performance: if the array type has efficient convert method to Array then this is
# almost a no-op
function _enforce_dense_cpu_array(img::AbstractArray)
    if Base.has_offset_axes(img)
        convert(Array, OffsetArrays.no_offset_view(img))
    else
        convert(Array, img)
    end
end
_enforce_dense_cpu_array(img::DenseArray) = img
_enforce_dense_cpu_array(img::OffsetArray) = _enforce_dense_cpu_array(parent(img))
_enforce_dense_cpu_array(img::IndirectArray) = img # PNGFiles has built-in support for this type

_get_bit_depth(img::BitArray) = 8  # TODO: write 1 bit-depth images
_get_bit_depth(img::AbstractArray{C}) where {C<:Colorant} = __get_bit_depth(eltype(C))
_get_bit_depth(img::AbstractArray{T}) where {T<:Normed} = __get_bit_depth(T)
# __get_bit_depth(::Type{Normed{T,1}}) where T = 1  # TODO: write 1 bit-depth images
# __get_bit_depth(::Type{Normed{T,2}}) where T = 2  # TODO: write 2 bit-depth images
# __get_bit_depth(::Type{Normed{T,4}}) where T = 4  # TODO: write 4 bit-depth images
__get_bit_depth(::Type{Normed{T,8}}) where T = 8
__get_bit_depth(::Type{Normed{T,16}}) where T = 16
__get_bit_depth(::Type{Normed{T,N}}) where {T,N} = ifelse(N <= 8, 8, 16)
__get_bit_depth(::Type{<:AbstractFloat}) = 8
_get_bit_depth(img::AbstractArray{T}) where {T<:AbstractFloat} = 8
_get_bit_depth(img::AbstractArray{<:Bool}) = 8  # TODO: write 1 bit-depth images
_get_bit_depth(img::AbstractArray{<:UInt8}) = 8
_get_bit_depth(img::AbstractArray{<:UInt16}) = 16

_get_color_type(x::AbstractArray{<:Gray}) = PNG_COLOR_TYPE_GRAY
_get_color_type(x::AbstractArray{<:GrayA}) = PNG_COLOR_TYPE_GRAY_ALPHA
_get_color_type(x::AbstractArray{<:AbstractRGB}) = PNG_COLOR_TYPE_RGB
_get_color_type(x::AbstractArray{<:AbstractARGB}) = PNG_COLOR_TYPE_RGBA
_get_color_type(x::AbstractArray{<:AbstractRGBA}) = PNG_COLOR_TYPE_RGBA
_get_color_type(x::IndirectArray) = PNG_COLOR_TYPE_PALETTE
function _get_color_type(
        x::AbstractArray{T, N}
    ) where {
        T<:Union{Normed,Unsigned,Bool,AbstractFloat},
        N
    }
    if N == 2
        return PNG_COLOR_TYPE_GRAY
    elseif N == 3
        d = size(x, 3)
        d == 1 && (return PNG_COLOR_TYPE_GRAY)
        d == 2 && (return PNG_COLOR_TYPE_GRAY_ALPHA)
        d == 3 && (return PNG_COLOR_TYPE_RGB)
        d == 4 && (return PNG_COLOR_TYPE_RGBA)
    end
    error("Number of dimensions $(N) in image not supported (only 2D or 3D Arrays are expected).")
end

_standardize_palette(p::AbstractVector{<:RGB}) = p
_standardize_palette(p::AbstractVector{<:AbstractRGB}) = RGB.(p)
_standardize_palette(p::AbstractVector{<:RGB{<:AbstractFloat}}) = RGB{N0f8}.(p)
_standardize_palette(p::AbstractVector{<:AbstractRGB{<:AbstractFloat}}) = RGB{N0f8}.(p)
_standardize_palette(p::OffsetArray) = _standardize_palette(parent(p))

_palette_alpha(p::OffsetArray) where {T} = _palette_alpha(collect(p))
_palette_alpha(p::AbstractVector{<:TransparentRGB{T,N0f8}}) where {T} = alpha.(p)
_palette_alpha(p::AbstractVector{<:TransparentRGB{T,<:AbstractFloat}}) where {T} = N0f8.(alpha.(p))


text:
export AshikhminShirleySingle, Lambert, LommelSeeliger

# Utility functions
costheta(w::Vector3) = w.z
sintheta2(w::Vector3) = max(0.0, 1.0 - costheta(w)^2)
sintheta(w::Vector3) = sqrt(sintheta2(w))
cosphi(w::Vector3) = costheta(w) ≈ 1.0 ? 1.0 : clamp(w.x / sintheta(w), -1.0, 1.0)
sinphi(w::Vector3) = costheta(w) ≈ 1.0 ? 0.0 : clamp(w.y / sintheta(w), -1.0, 1.0)

"""
    evaluate(B::BxDF, world_to_local::Transformation, w0::Vector3, w1::Vector3)

Evaluate the value of a BRDF given a local transformation and the incident and emergent
directions. The local transformation is one that rotates the surface normal to coincide
with the +Z direction.

"""
function evaluate(B::BxDF, world_to_local::Transformation, w0::Vector3, w1::Vector3)
    return evaluate(B, world_to_local(w0), world_to_local(w1))
end

######################################################################

"""
    Lambert{T<:Spectrum} <: BxDF

The Lambertian BRDF with spectrum `R`.

"""
struct Lambert{T<:Spectrum} <: BxDF
    R::T
end

"""
    evaluate(B::Lambert, w0::Vector3, w1::Vector3)

Evaluate the Lambertian BRDF for given direction vectors.

"""
evaluate(B::Lambert, w0::Vector3, w1::Vector3) = B.R ./ π

""""""
rho(B::Lambert, w0::Vector3) = B.R

"""
    generate_ray(B::Lambert, w0::Vector3) = nothing

Generate a new ray, weighed by the emergent light distribution of the Lambertian BRDF.

TODO: not implemented yet

"""
generate_ray(B::Lambert, w0::Vector3) = nothing

"""
    compute_pressure(B::Lambert, w0::Vector3, S::Spectrum)

Compute the radiation pressure on a surface

"""
function compute_pressure(B::Lambert, w0::Vector3, S::Spectrum)
    if costheta(w0) < 0.0
        return Vector3(0)
    end
    total_energy = integrate(B.R .* S)
    return -total_energy * (w0 + 2/3 * Z_AXIS)
end

######################################################################

"""
    LommelSeeliger{T<:Spectrum} <: BxDF

Lommel-Seeliger BRDF. Contains the surface spectrum and a phase function. The phase
function needs to be a one-parameter function defined on [0,π] and it's assumed to be
normalized so that its integral over the sphere equals 4π.

"""
struct LommelSeeliger{T<:Spectrum} <: BxDF
    R::T
    P::Function
end

"""
    evaluate(B::LommelSeeliger, w0::Vector3, w1::Vector3)

Evaluate Lommel-Seeliger BRDF for given direction vectors.

TODO: check physics

"""
function evaluate(B::LommelSeeliger, w0::Vector3, w1::Vector3)
    cos_alpha = clamp(dot(w0, w1), -1.0, 1.0)
    alpha = acos(cos_alpha)
    phase = B.P(alpha)
    return 0.25 * phase * B.R * max(0.0, 1.0 / (costheta(w0) + costheta(w1)))
end

SpectrumLike = Union{T,R} where T <: Spectrum where R <: Real

######################################################################

"""
    AshikhminShirleySingle{T<:Spectrum} <: BxDF

The Ashikhmin-Shirley BRDF where `R` is a spectrum, `d` is the diffuse component weigth
and `n` is the specular reflection width parameter. See Wetterer (2014).

"""
struct AshikhminShirleySingle{S<:SpectrumLike, T<:SpectrumLike} <: BxDF
    R::T
    specularity::S
    peak::Float64
end

"""
    Fresnel(R::Spectrum, costheta::Real, s::Real)

The Fresnel term of the Ashikhmin-Shirley BRDF (see Wetterer, 2014).
"""
@inline Fresnel(Rs::SpectrumLike, x::Real) = Rs .+ (1.0 .- Rs) .* (1.0 - x)^5

"""
    BlinnPhong(n::Real, hdn::Real)

The Blinn-Phong function used in the Ashikhmin-Shirley BRDF (see Wetterer, 2014).
"""
@inline BlinnPhong(n::SpectrumLike, hdn::Real) = (n .+ 1) .* hdn^n ./ 8π

@inline f5(x) = (one(x) - (one(x) - x/2)^5)

"""
    evaluate(B::AshikhminShirleySingle, w0::Vector3, w1::Vector3)

Evaluate Ashikhmin-Shirley BRDF for given direction vectors.

"""
function evaluate(B::AshikhminShirleySingle, w0::Vector3, w1::Vector3)
    h = normalize((w0 + w1) / 2.0)
    D = BlinnPhong(B.peak, h.z)
    F = Fresnel(B.R, costheta(w1))
    specular = F .* (D ./ (dot(h,w1) * max(costheta(w0), costheta(w1))))
    diffuse = B.R .* (1 .- B.specularity) .* (28/23π * f5(costheta(w0)) * f5(costheta(w1)))

    return diffuse .+ specular
end

######################################################################

"""
    SpecularDiffuse{T<:Spectrum} <: BxDF

Specular + Diffuse reflection. The fraction of diffusely reflected light is `fd`, the
fraction of specularly reflected light is `fs`. Absorption is `1 - fd - fs`.

TODO: Make all the parameters spectral.

"""
struct SpecularDiffuse{T<:Spectrum} <: BxDF
    R::T
    fd::Float64
    fs::Float64
end

"""
    evaluate(B::SpecularDiffuse, w0::Vector3, w1::Vector3)

Evaluate diffuse + specular BRDF for given direction vectors.

TODO: normalization is probably off...

"""
function evaluate(B::SpecularDiffuse, w0::Vector3, w1::Vector3)
    if normalize(w0 + w1) ≈ Z_AXIS
        return (B.fd / π + B.fs) .* B.R
    else
        return B.fd / π .* B.R
    end
end

"""
    compute_pressure(B::SpecularDiffuse, w0::Vector3)

Compute the radiation pressure on a surface.

"""
function compute_pressure(B::SpecularDiffuse, w0::Vector3, S::Spectrum)
    if costheta(w0) < 0.0
        return Vector3(0)
    end
    fa = 1.0 - B.fd - B.fs
    total_energy = integrate(B.R .* S)
    return -total_energy * ((fa + B.fd) * (w0 + 2/3 * Z_AXIS) + 2*B.fs*costheta(w0) * Z_AXIS)
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
module AlignmentAlgorithms

export nwscore,
       nwsetupmatrix,
       nwscorematrix
#export nwalign

include("needleman_wunsch.jl")

end


text:
#=

Functions related to capturing the past states of the notebook

We always deepcopy variables before saving them in the archive.
So, we use the routine can_copy to determine if deepcopy will work on them.
It is not perfect.

To capture all current variables in a dictionary, we use `main_to_dict`.

=#

"""
Go over every symbol in Main.
If `deepcopy` works on it, put it in the dict `d`.
Note that it doesn't capture functions, which is unfortunate.
"""
function main_to_dict()
    d = Dict{Any,Any}()
    for n in names(Main)
        # you would expect that every name in Main corresponds to a defined variable, but this assumption can fail.
        if isdefined(Main, n)
            val = @eval Main $n
            if !(objectid(val) ∈ DontSave)
                can, h = can_copy_and_hash(val)
                if can
                    copyval = haskey(VX.store, h) ? nothing : deepcopy(isa(val, Union{Expr,Symbol}) ? Meta.quot(val) : val)
                    d[n] = (h,copyval)
                end
            end
        end
    end

    return d
end

"""
Save the current variables and ans in VX (a Varchive).
"""
function save_state()
    # if this was spawned by TM.@thread, then its variables will be saved later.
    if !(IJulia.n ∈ spawned)
        di = main_to_dict()

        # get a version of ans we can copy, if it exists
        ans = IJulia.ans
        ansc = can_copy(ans) ? deepcopy(isa(ans,Union{Expr,Symbol}) ? Meta.quot(ans) : ans) : nothing

        put_state!(VX, IJulia.n, di, ansc)

        if Base.summarysize(VX) > Sys.total_memory()/3
            @warn "IJuliaTimeMachine state takes over 1/3 of system memory.  Consider IJuliaTimeMachine.clear_state()."
        end
        debug_mode && println(IJulia.orig_stdout[], "Saved state $(IJulia.n)")
    end
end

"""
    @past n

Return to the state after cell `n` was exectued.
"""
macro past(n_in)

    n = n_in
    if typeof(n_in) == Symbol
        n = @eval Main $(n_in)
    end

    if haskey(VX.past, n)
        local s = IJuliaTimeMachine.vars(n)
        local ans = VX.past[n].ans
        return quote
            ($([esc(x) for x in keys(s)]...),) =
                ($([esc(x) for x in values(s)]...),)
            $ans
        end
    else
        return quote
            error("State $($n) was not saved.")
        end
    end
end

"""
    @dict_to_main dict

Place all variables in dict into Main memory.
"""
macro dict_to_main(dict)
    s = @eval Main $(dict)
    return quote
        ($([esc(x) for x in keys(s)]...),) =
            ($([esc(x) for x in values(s)]...),)
        nothing
    end
end

"""
    clear_past(indices)

Empty all storage of the past.  Use to free up memory.
If indices is omitted, it clears all history.
"""
function clear_past()
    empty!(VX.store)
    empty!(VX.past)
    return nothing
end

function clear_past(indices)
    # first, mark all variables that are not in indices
    keepvar = Set()
    for (n, vs) in VX.past
        if !(n ∈ indices)
            for h in values(vs.vars)
                push!(keepvar,h)
            end
        end
    end

    # now, remove from store every unmarked variable
    for k in keys(VX.store)
        if !(k ∈ keepvar)
            delete!(VX.store, k)
        end
    end

    for k ∈ indices
        delete!(VX.past, k)
    end
end

vars(n::Int) = haskey(VX.past,n) ? vars(VX, n) : error("Cell $(n) was not saved.")

"""
Recall the answer from a cell
"""
ans(n::Int) = haskey(VX.past,n) ? VX.past[n].ans : error("Cell $(n) was not saved.")

max_stars_repo_name: danspielman/IJuliaTimeMachine

max_issues_repo_name: danspielman/IJuliaTimeMachine

max_forks_repo_name: danspielman/IJuliaTimeMachine


text:
# Copyright (c) 2019 Arpit Bhatia and contributors                               #src
#                                                                                #src
# Permission is hereby granted, free of charge, to any person obtaining a copy   #src
# of this software and associated documentation files (the "Software"), to deal  #src
# in the Software without restriction, including without limitation the rights   #src
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      #src
# copies of the Software, and to permit persons to whom the Software is          #src
# furnished to do so, subject to the following conditions:                       #src
#                                                                                #src
# The above copyright notice and this permission notice shall be included in all #src
# copies or substantial portions of the Software.                                #src
#                                                                                #src
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     #src
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       #src
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    #src
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         #src
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  #src
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  #src
# SOFTWARE.                                                                      #src

# # Working with Data Files

# **Originally Contributed by**: Arpit Bhatia

# In many cases we might need to read data available in an external file rather
# than type it into Julia ourselves.

# This tutorial is concerned with reading tabular data into Julia and using it
# for a JuMP model.

# We'll be reading data using the [DataFrames.jl](https://github.com/JuliaData/DataFrames.jl)
# package and some other packages specific to file types.

# The data are stored in the [`/docs/src/tutorials/Getting started/data`](https://github.com/jump-dev/JuMP.jl/tree/master/docs/src//tutorials/Getting started/data)
# directory of the JuMP source code.

const DATA_DIR = joinpath(@__DIR__, "data");

# !!! note
#     There are multiple ways to read the same kind of data into Julia. This
#     tutorial focuses on DataFrames.jl because it provides the ecosystem to
#     work with most of the required file types in a straightforward manner.

# ### DataFrames.jl

# The `DataFrames` package provides a set of tools for working with tabular
# data. It is available through the Julia package system.
# ```julia
# using Pkg
# Pkg.add("DataFrames")
# ```

import DataFrames

# ### What is a DataFrame?

# A DataFrame is a data structure like a table or spreadsheet. You can use it
# for storing and exploring a set of related data values. Think of it as a
# smarter array for holding tabular data.

# ## Reading Tabular Data into a DataFrame

# We will begin by reading data from different file formats into a DataFrame
# object.

# ### Excel Sheets

# Excel files can be read using the [XLSX.jl](https://github.com/filepenoris/XLSX.jl)
# package.

# ```julia
# Pkg.add("XLSX")
# ```

import XLSX

# To read a Excel file into a DataFrame, we use the following julia code. The
# first arguement to the `readtable` function is the file to be read and the
# second arguement is the name of the sheet.

excel_df = DataFrames.DataFrame(
    XLSX.readtable(joinpath(DATA_DIR, "SalesData.xlsx"), "SalesOrders")...
)

# ### CSV files

# CSV and other delimited text files can be read by the CSV.jl package.

# ```julia
# Pkg.add("CSV")
# ```

import CSV

# To read a CSV file into a DataFrame, we use the `CSV.read` function.

csv_df = CSV.read(joinpath(DATA_DIR, "StarWars.csv"), DataFrames.DataFrame)

# ### Other Delimited Files

# We can also use the `CSV.jl` package to read any other delimited text file
# format.

# By default, CSV.File will try to detect a file's delimiter from the first 10
# lines of the file.

# Candidate delimiters include `','`, `'\t'`, `' '`, `'|'`, `';'`, and `':'`. If
# it can't auto-detect the delimiter, it will assume `','`.

# Let's take the example of space separated data.

ss_df = CSV.read(joinpath(DATA_DIR, "Cereal.txt"), DataFrames.DataFrame)

# We can also specify the delimiter by passing the `delim` arguement.

delim_df = CSV.read(
    joinpath(DATA_DIR, "Soccer.txt"), DataFrames.DataFrame, delim = "::"
)

# Note that by default, are read-only. If we wish to make changes to the data
# read, we pass the `copycols = true` arguement in the function call.

ss_df = CSV.read(
    joinpath(DATA_DIR, "Cereal.txt"), DataFrames.DataFrame, copycols = true
)

# ## Working with DataFrames

# Now that we have read the required data into a DataFrame, let us look at some
# basic operations we can perform on it.

# ### Querying Basic Information

# The `size` function gets us the dimensions of the DataFrame.

DataFrames.size(ss_df)

# We can also us the `nrow` and `ncol` functions to get the number of rows and
# columns respectively.

DataFrames.nrow(ss_df), DataFrames.ncol(ss_df)

# The `describe` function gives basic summary statistics of data in a DataFrame.

DataFrames.describe(ss_df)

# Names of every column can be obtained by the `names` function.

DataFrames.names(ss_df)

# Corresponding data types are obtained using the broadcasted `eltype` function.

eltype.(ss_df)

# ### Accessing the Data

# Similar to regular arrays, we use numerical indexing to access elements of a
# DataFrame.

csv_df[1, 1]

# The following are different ways to access a column.

csv_df[!, 1]

#-

csv_df[!, :Name]

#-

csv_df.Name

#-

csv_df[:, 1] # Note that this creates a copy.

# The following are different ways to access a row.

csv_df[1:1, :]

#-

csv_df[1, :] # This produces a DataFrameRow.

# We can change the values just as we normally assign values.

# Assign a range to scalar.

excel_df[1:3, 5] .= 1

# Vector to equal length vector.

excel_df[4:6, 5] = [4, 5, 6]

# Subset of the DataFrame to another data frame of matching size.

excel_df[1:2, 6:7] =  DataFrames.DataFrame(
    [-2 -2; -2 -2], [Symbol("Unit Cost"), :Total]
)

#-

excel_df

# !!! tip
#     There are a lot more things which can be done with a DataFrame. Read the
#     [docs](https://juliadata.github.io/DataFrames.jl/stable/) for more
#     information.

# ## A Complete Modelling Example - Passport Problem

# Let's now apply what we have learnt to solve a real modelling problem.

# The [Passport Index Dataset](https://github.com/ilyankou/passport-index-dataset)
# lists travel visa requirements for 199 countries, in `.csv` format. Our task
# is to find out the minimum number of passports required to visit all
# countries.

# In this dataset, the first column represents a passport (=from) and each
# remaining column represents a foreign country (=to).

# The values in each cell are as follows:
# * 3 = visa-free travel
# * 2 = eTA is required
# * 1 = visa can be obtained on arrival
# * 0 = visa is required
# * -1 is for all instances where passport and destination are the same

# Our task is to find out the minimum number of passports needed to visit every
# country without requiring a visa.

# Thus, the values we are interested in are -1 and 3. Let us modify the data in
# the following manner:

passport_data = CSV.read(
    joinpath(DATA_DIR, "passport-index-matrix.csv"),
    DataFrames.DataFrame;
    copycols = true,
)

for i in 1:DataFrames.nrow(passport_data)
    for j in 2:DataFrames.ncol(passport_data)
        if passport_data[i, j] == -1 || passport_data[i, j] == 3
            passport_data[i, j] = 1
        else
            passport_data[i, j] = 0
        end
    end
end

# The values in the cells now represent:
# * 1 = no visa required for travel
# * 0 = visa required for travel

# Let us assossciate each passport with a decision variable $pass_{cntr}$ for
# each country. We want to minize the sum $\sum pass_{cntr}$ over all countries.

# Since we wish to visit all the countries, for every country, we should own at
# least one passport that lets us travel to that country visa free. For one
# destination, this can be mathematically represented as
# $\sum_{cntr \in world} passportdata_{cntr,dest} \cdot pass_{cntr} \geq 1$.

# Thus, we can represent this problem using the following model:

# ```math
# \begin{aligned}
# \min && \sum_{cntr \in World} pass_{cntr} \\
# \text{s.t.} && \sum_{cntr \in World} passportdata_{cntr,dest} \cdot pass_{cntr} \geq 1 && \forall dest \in World \\
# && pass_{cntr} \in \{0,1\} && \forall cntr \in World
# \end{aligned}
# ```

# We'll now solve the problem using JuMP.

using JuMP
import GLPK

# First, create the set of countries:

World = names(passport_data)[2:end]

# Then, create the model and intialize the decision variables:
model = Model(GLPK.Optimizer)
@variable(model, pass[cntr in World], Bin)

# Define the objective function

@objective(model, Min, sum(pass[cntr] for cntr in World))

#-

@constraint(model, [dest in World], passport_data[:, dest]' * pass >= 1)

# Now optimize!

optimize!(model)
println("Minimum number of passports needed: ", objective_value(model))

#-

optimal_passports = [cntr for cntr in World if value(pass[cntr]) > 0.5]
println("Countries:")
for p in optimal_passports
    println(" ", p)
end

# !!! note
#     We use `value(pass[i]) > 0.5` rather than `value(pass[i]) == 1` to avoid
#     excluding solutions like `pass[i] = 0.99999` that are "1" to some
#     tolerance.

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# open file and store content into a string
input = open("inputs//day1.txt", "r") do io
    read(io, String)
end

# split the string into a vector and parse it as integers
data = parse.(Int, split(input))

##### part 1 #####

# bit-sum of all the "true" when checking for positive in the derivation of data
result1 = sum(0 .< diff(data))

println("part1: ", result1)

##### part 2 #####

# generate the new data (gliding window of 3 terms sum)
data3 = zeros(Int, length(data) - 2)
for i = 1:length(data)-2
    data3[i] = sum(data[i:i+2])
end

# bit-sum of all the "true" when checking for positive in the derivation of data
result2 = sum(0 .< diff(data3))

println("part2: ", result2)


text:
using GeometricFlux
using GeometricFlux: neighbors, laplacian_matrix
using LightGraphs: SimpleGraph, SimpleDiGraph, add_edge!
using SimpleWeightedGraphs: SimpleWeightedGraph, SimpleWeightedDiGraph, add_edge!
using Zygote
using Test

using CUDAapi
if has_cuda()
    try
        using CuArrays
        @eval has_cuarrays() = true
    catch ex
        @warn "CUDA is installed, but CuArrays.jl fails to load" exception=(ex,catch_backtrace())
        @eval has_cuarrays() = false
    end
else
    has_cuarrays() = false
end

cuda_tests = [
    "cuda/scatter",
    "cuda/pool",
    "cuda/grad",
    "cuda/conv",
    "cuda/msgpass",
]

tests = [
    "layers/msgpass",
    "layers/conv",
    "layers/pool",
    "grad",
    "models",
    "linalg",
    "scatter",
    "graph/simplegraphs",
    "graph/weightedgraphs",
    "graph/utils",
    "utils",
]

if has_cuarrays()
    using Flux: gpu
    append!(tests, cuda_tests)
end

@testset "GeometricFlux" begin
    for t in tests
        include("$(t).jl")
    end
end


text:
using Test, YaoBase

@testset "test error & exceptions" begin
    include("error.jl")
    include("exceptions.jl")
end

@testset "test abstract interface" begin
    include("abstract_register.jl")
end

@testset "test utils" begin
    include("math.jl")
end