text:
"""
    Submodel{iip, F, NS, NI, NP, JR, JS, JI, JP, GT, IS, II, IP, OS, OI, OP}

Base type for submodels in `Aeroelasticity.jl`.
"""
struct Submodel{iip, F, NS, NI, NP, C, JR, JS, JI, JP, GT, IS, II, IP, OS, OI, OP}
    # residual function
    f::F
    # dimensions
    nx::NS
    ny::NI
    np::NP
    # constants
    constants::C
    # jacobians/gradients
    ratejac::JR
    statejac::JS
    inputjac::JI
    paramjac::JP
    tgrad::GT
    # input/output functions
    setstate::IS
    setinput::II
    setparam::IP
    sepstate::OS
    sepinput::OI
    sepparam::OP
end

"""
    Submodel{iip}(f, nx, ny, np; kwargs...)

Define a custom submodel for use with `Aeroelasticity.jl`.  Note that `iip` should be set to
`true` if `f` is defined in-place and `false` otherwise.

# Arguments
 - `f`: Residual function (`resid = f(dx, x, y, p, t)` or `f(resid, dx, x, y, p, t)`)
 - `nx`: Number of states.
 - `ny`: Number of inputs.
 - `np`: Number of parameters.

# Keyword Arguments
 - `ratejac = Nonlinear()`:  Rate jacobian description (of type [`AbstractJacobian`](@ref))
 - `statejac = Nonlinear()`: State jacobian description (of type [`AbstractJacobian`](@ref))
 - `inputjac = Nonlinear()`: Input jacobian description (of type [`AbstractJacobian`](@ref))
 - `paramjac = Nonlinear()`: Parameter jacobian description (of type [`AbstractJacobian`](@ref))
 - `tgrad = Nonlinear()`: Time gradient description (of type [`AbstractGradient`](@ref))
 - `setstate = (states; x) -> states .= x`: Function for setting rates/states using a
    custom set of keyword arguments.
 - `setinput = (inputs; y) -> inputs .= y`: Function for setting inputs using a custom set
    of keyword arguments.
 - `setparam = (params; p) -> params .= p`: Function for setting parameters using a custom
    set of keyword arguments.
 - `sepstate = (x) -> (x = x,)`: Function for separating rates/states into named variables.
 - `sepinput = (y) -> (y = y,)`: Function for separating inputs into named variables.
 - `sepparam = (p) -> (p = p,)`: Function for separating parameters into named variables.
"""
function Submodel{iip}(f, nx, ny, np;
    constants = (),
    ratejac = Nonlinear(),
    statejac = Nonlinear(),
    inputjac = Nonlinear(),
    paramjac = Nonlinear(),
    tgrad = Nonlinear(),
    setstate = (states; x) -> states .= x,
    setinput = (inputs; y) -> inputs .= y,
    setparam = (params; p) -> params .= p,
    sepstate = (x) -> (x = x,),
    sepinput = (y) -> (y = y,),
    sepparam = (p) -> (p = p,)) where iip

    # add jacobian definitions, if necessary
    ratejac = add_rate_jacobian(ratejac, f, iip, nx, ny, np)
    statejac = add_state_jacobian(statejac, f, iip, nx, ny, np)
    inputjac = add_input_jacobian(inputjac, f, iip, nx, ny, np)
    paramjac = add_parameter_jacobian(paramjac, f, iip, nx, ny, np)
    tgrad = add_time_gradient(tgrad, f, iip, nx, ny, np)

    return Submodel{iip,
        typeof(f), typeof(nx), typeof(ny), typeof(np), typeof(constants),
        typeof(ratejac), typeof(statejac), typeof(inputjac), typeof(paramjac), typeof(tgrad),
        typeof(setstate), typeof(setinput), typeof(setparam),
        typeof(sepstate), typeof(sepinput), typeof(sepparam)}(
            f, nx, ny, np, constants,
            ratejac, statejac, inputjac, paramjac, tgrad,
            setstate, setinput, setparam,
            sepstate, sepinput, sepparam)
end

"""
    Submodel(np; kwargs...)

Define a submodel with no state variables (see [`Submodel`](@ref))

# Arguments
 - `np`: Number of parameters.

# Keyword Arguments
 - `constants = ()`: Constants associated with this model
 - `setparam = (params; p) -> params .= p`: Function for setting parameters using a custom
    set of keyword arguments.
 - `sepparam = (p) -> (p = p,)`: Function for separating parameters into named variables.
"""
function Submodel(np;
    constants = (),
    setparam = (params; p) -> params .= p,
    sepparam = (p) -> (p = p,))

    return Submodel{false}((dx, x, y, p, t) -> Float64[], Val(0), Val(0), np;
        constants = constants,
        ratejac = Empty(),
        statejac = Empty(),
        inputjac = Empty(),
        paramjac = Empty(),
        tgrad = Empty(),
        setstate = (x) -> x,
        setinput = (y) -> y,
        setparam = setparam,
        sepstate = (x) -> (),
        sepinput = (y) -> (),
        sepparam = sepparam)
end

"""
    Coupling{iip, G, NX, NY, NP, NPC, JR, JS, JP, GT, IP, OP}

Base type for coupling models in Aeroelasticity.jl.
"""
struct Coupling{iip, G, NX, NY, NP, NPC, JR, JS, JP, GT, IP, OP}
    # coupling function
    g::G
    # dimensions
    nx::NX
    ny::NY
    np::NP
    npc::NPC
    # jacobians/gradients
    ratejac::JR
    statejac::JS
    paramjac::JP
    tgrad::GT
    # input/output functions
    setparam::IP
    sepparam::OP
end

"""
    Coupling{iip}(g, nx, ny, np, npc; kwargs...)

Define a custom coupling model for use with Aeroelasticity.jl.  Note that `iip` should
be set to `true` if `g` is defined in-place and `false` otherwise.

# Arguments
 - `g`: Coupling function (`y = g(dx, x, p, t)` or `g(y, dx, x, p, t)`)
 - `nx`: Number of states.
 - `ny`: Number of inputs.
 - `np`: Number of parameters (including parameters introduced by the coupling).
 - `npc`: Number of additional parameters introduced by the coupling.

# Keyword Arguments
 - `ratejac = Nonlinear()`:  Rate jacobian description (of type [`AbstractJacobian`](@ref))
 - `statejac = Nonlinear()`: State jacobian description (of type [`AbstractJacobian`](@ref))
 - `paramjac = Nonlinear()`: Parameter jacobian description (of type [`AbstractJacobian`](@ref))
 - `tgrad = Nonlinear()`: Time gradient description (of type [`AbstractGradient`](@ref))
 - `setparam = (params; p) -> params .= p`: Function for setting parameters introduced by
    the coupling using a custom set of keyword arguments.
 - `sepparam = (p) -> (p = p,)`: Function for separating parameters introduced by the
    coupling into named variables.
"""
function Coupling{iip}(g, nx, ny, np, npc;
    ratejac = Nonlinear(),
    statejac = Nonlinear(),
    paramjac = Nonlinear(),
    tgrad = Nonlinear(),
    setparam = (params; p) -> params .= p,
    sepparam = (p) -> (p = p,)) where {iip}

    # add jacobian definitions, if necessary
    ratejac = add_coupling_rate_jacobian(ratejac, g, iip, nx, ny, np)
    statejac = add_coupling_state_jacobian(statejac, g, iip, nx, ny, np)
    paramjac = add_coupling_parameter_jacobian(paramjac, g, iip, nx, ny, np)
    tgrad = add_coupling_time_gradient(tgrad, g, iip, nx, ny, np)

    return Coupling{iip,
        typeof(g), typeof(nx), typeof(ny), typeof(np), typeof(npc),
        typeof(ratejac), typeof(statejac), typeof(paramjac), typeof(tgrad),
        typeof(setparam), typeof(sepparam)}(
            g, nx, ny, np, npc,
            ratejac, statejac, paramjac, tgrad,
            setparam, sepparam)
end

"""
    CoupledModel{M,C}

Base type for coupled models in Aeroelasticity.jl.
"""
struct CoupledModel{M,C}
    submodels::M
    coupling::C
end

"""
    CoupledModel(submodels, coupling)

Define a coupled model given a tuple of submodels (of type [`Submodel`](@ref)) and a
coupling model (of type [`Coupling`](@ref)).
"""
CoupledModel(submodels, coupling)

# --- Inplaceness Trait --- #

abstract type InPlaceness end
struct InPlace <: InPlaceness end
struct OutOfPlace <: InPlaceness end

inplaceness(::Submodel{true}) = InPlace()
inplaceness(::Submodel{false}) = OutOfPlace()

inplaceness(::Coupling{true}) = InPlace()
inplaceness(::Coupling{false}) = OutOfPlace()

function inplaceness(model::CoupledModel)
    if any(isinplace.(model.submodels)) || isinplace(model.coupling)
        return InPlace()
    else
        return OutOfPlace()
    end
end

isinplace(x) = isinplace(inplaceness(x))
isinplace(::InPlace) = true
isinplace(::OutOfPlace) = false


text:
include("Types.jl")
include("Network.jl")
include("Base.jl")
include("Files.jl")
include("Compute.jl")

ext: jl

max_stars_repo_name: Pooksoft/BinomialEquityAssetPriceModel

max_issues_repo_name: Pooksoft/BinomialEquityAssetPriceModel

max_forks_repo_name: Pooksoft/BinomialEquityAssetPriceModel


text:
Random.seed!(43)

X = randn(100, 20)
Xup = rand(25, 20)
l = size(Xup, 1) + 1
Xprim = vcat(X, Xup)[l:end,:]
@testset "data updat" begin
  @test dataupdat(X, Xup) ≈ Xprim
end

@testset "moment updates" begin
  x = ones(6, 2)
  y = 2*ones(2,2)
  M3 = moment(x, 3)
  M4 = moment(x, 4)
  M3up = momentupdat(M3, x, y)
  @testset "simple test" begin
    Mup = moment(dataupdat(x,y),3)
    @test Array(Mup) ≈ Array(M3up)
    x0 = Array{Float64}(undef, (0,0))
    y0 = Array{Float64}(undef, (0,0))
    @test M3 == momentupdat(M3, x0, y0)
  end
  @testset "moment array" begin
    Ma = momentarray(x, 4, 2)
    @test Array(Ma[3]) ≈ Array(M3)
    @test Array(Ma[4]) ≈ Array(M4)
    MM = momentupdat(Ma, x, y)
    @test Array(M3up) ≈Array(MM[3])
  end
end

@testset "moment exceptions" begin
  x = ones(10,4);
  y = 2*ones(5,3);
  m = moment(x, 3);
  @test_throws DimensionMismatch momentupdat(m, x, y)
  y = 2*ones(5,4)
  @test_throws DimensionMismatch momentupdat(m, x[:, 1:3], y)
  y = 2*ones(15,4)
  @test_throws BoundsError momentupdat(m, x, y)
end

@testset "moments to cumulants" begin
  m1 = moment(X, 1)
  m2 = moment(X, 2)
  m3 = moment(X, 3)
  m4 = moment(X, 4)
  m5 = moment(X, 5)
  c = cumulants(X, 5)
  m2c = [m1, m2, m3, m4, m5]
  @testset "moms2cums!" begin
    moms2cums!(m2c)
    @test Array(c[1]) ≈ Array(m2c[1])
    @test Array(c[2]) ≈ Array(m2c[2])
    @test Array(c[3]) ≈ Array(m2c[3])
    @test Array(c[4]) ≈ Array(m2c[4])
    @test Array(c[5]) ≈ Array(m2c[5])
  end
  @testset "cums2moms" begin
    mm = cums2moms(c);
    @test Array(mm[1]) ≈ Array(m1)
    @test Array(mm[2]) ≈ Array(m2)
    @test Array(mm[3]) ≈ Array(m3)
    @test Array(mm[4]) ≈ Array(m4)
    @test Array(mm[5]) ≈ Array(m5)
  end
end

@testset "simple cumulants update" begin
  x = ones(10, 3)
  y = 2*ones(2,3)
  s = DataMoments(x, 6, 2)
  c2 = cumulantsupdate!(s, y)
  xx = dataupdat(x,y)
  c3 = cumulants(xx, 6)
  @test Array(c2[1]) ≈ Array(c3[1])
  @test Array(c2[2]) ≈ Array(c3[2])
  @test Array(c2[3]) ≈ Array(c3[3])
  @test Array(c2[4]) ≈ Array(c3[4])
  @test Array(c2[5]) ≈ Array(c3[5])
  @test Array(c2[6]) ≈ Array(c3[6])
end

@testset "cumulants updates larger data" begin
    c = cumulants(X)
    Xp = dataupdat(X, Xup)
    s = DataMoments(X, 4, 4)
    C = cumulantsupdate!(s, Xup)
    @test s.X ≈ Xp
    @test s.b == 4
    @test s.d == 4
    @test Array(s.M[4]) ≈ Array(moment(Xp, 4))
    CC = cumulants(Xp, 4)
    @test Array(C[3]) ≈ Array(CC[3])
    @test Array(C[4]) ≈ Array(CC[4])
end

@testset "multiprocessing cumulants update" begin
  addprocs(2)
  eval(Expr(:toplevel, :(@everywhere using CumulantsUpdates)))
    Xp = dataupdat(X, Xup)
    s = DataMoments(X, 4, 4)
    C = cumulantsupdate!(s, Xup)
    @test s.X ≈ Xp
    @test s.b == 4
    @test s.d == 4
    @test Array(s.M[4]) ≈ Array(moment(Xp, 4))
    CC = cumulants(Xp, 4)
    @test Array(C[3]) ≈ Array(CC[3])
    @test Array(C[4]) ≈ Array(CC[4])
end

@testset "cumulants update exceptions" begin
  x = ones(10,4);
  y = 2*ones(5,3);
  s = DataMoments(x, 4, 2)
  @test_throws DimensionMismatch cumulantsupdate!(s, y)
  y = 2*ones(15,4)
  @test_throws BoundsError cumulantsupdate!(s, y)
  s1 = DataMoments(x[:, 1:3], 4, 2)
  y = 2*ones(5,4)
  @test_throws UndefVarError cumulantsupdat(s1, y)
end

println(VERSION)
if VERSION <= v"1.6.1"
  @testset "save and load" begin
    x = ones(10,4);
    s = DataMoments(x, 4, 2)
    @test savedm(s, "/tmp/cumdata.jld2") == nothing
    #println(load("/tmp/cumdata.jld2"))
    s1 = loaddm("/tmp/cumdata.jld2")
    @test s1.X == s.X
    @test s1.d == s.d
    @test s1.b == s.b
    @test Array(s1.M[4]) == Array(s.M[4])
  end
end


text:
export MeanCov, MeanVar
import Base.iterate

#=
Online statistics

m   # mean(x[1:n])
m2  #  sum of squares of differences from the current mean, ``\textstyle\sum_{i=1}^n (x_i - \bar x_n)^2``
n   # number of iterations

Implemented as "dependent iterator" where next has an additional argument
=#

_diag(x) = diag(x)
_diag(x::Number) = x

"""
    mcstart(x) -> state

Create state for random chain online statitics. The entries/value of `x` are ignored
"""
mcstart(yy::Array{T}) where {T} = (zeros(T, size(yy))/1, zeros(typeof(outer(zero(T))), size(yy))/1, 0)
mcstart(y::T) where {T<:Number} = (zero(T)/one(T), zero(T)/one(T), 0)

"""
    mcnext(state, x) -> state

Update random chain online statistics when new chain value `x` was
observed. Return new `state`.
"""
function mcnext(mc, x)
    m, m2, n = mc
    delta = x - m
    m = m + delta/(n + 1)
    delta2 = x - m
    m2 = m2 + delta .* delta2
    m, m2, n + 1
end

function mcnext(mc,  x::Vector{<:AbstractArray}) # fix me: use covariance
    m, m2, n = mc
    delta = x - m
    m = m + delta/(n+1)
    m2 = m2 + map(outer, delta, (x - m))
    m, m2, n + 1
end

function mcnext!(mc,  x::Vector{<:AbstractArray})
    m, m2, n = mc
    for i in 1:length(m2)
        delta = x[i] - m[i]
        m[i]  += (delta)/(n+1)
        m2[i] += outer(delta, x[i]-m[i])
    end
    m, m2, n + 1
end

"""
    mcmeanband(mc)

Compute marginal confidence interval for the chain *mean* using normal approximation
"""
function mcbandmean(mc)
    m, m2, k = mc
    Q = sqrt(2.)*erfinv(0.95)

    ste = eltype(m)[sqrt.(_diag(v)*(1/(k - 1))) for v in m2 ]*sqrt(1/k)
    m - Q*ste, m + Q*ste
end

"""
    mcband(mc)

Compute marginal 95% coverage interval for the chain from normal approximation.
"""
function mcband(mc)
    m, m2, k = mc
    Q = sqrt(2.)*erfinv(0.95)

    std = eltype(m)[sqrt.(_diag(v) * (1/(k - 1))) for v in m2]
    m-Q*std, m+Q*std
end

"""
    mcstats(mc)

Compute mean and covariance estimates.
"""
function mcstats(mc)
    m, m2, k = mc
    cov = m2/(k - 1)
    m, cov
end

"""
    mcmarginalstats(mcstates) -> mean, std

Compute `mean`` and marginal standard deviations `std` for 2d plots.
"""
function mcmarginalstats(states)
    xx, vv = Bridge.mcstats(states[1])
    Xmean = copy(xx)
    Xstd = map(x->sqrt.(diag(x)), vv)
    for i in 2:length(states)
        pop!(Xmean); pop!(Xstd)
        xx, vv = Bridge.mcstats(states[i])
        append!(Xmean, xx)
        append!(Xstd, map(x->sqrt.(diag(x)), vv))
    end
    Xmean, Xstd
end

"""
    MeanCov(itr)

Iterator interface for online mean and covariance estimates.
Iterates are triples `mean, λ, cov/λ`

Example:

```
c = Channel{Vector}(1)
m = Bridge.MeanCov(c)
put!(c, rand(5))
u = iterate(m)
put!(c, rand(5))
u = iterate(m, u[2])
close(c)
u[1][1]
```

```
m = Bridge.MeanCov(Channel{Vector{Float64}}(1))
u = register!(m, rand(5))
u = register!(m, rand(5), u)
close(m)
u[1][1]
```

"""
struct MeanCov{T}
    iter::T
end
import Base: length, IteratorSize
length(x::MeanCov) = length(x.iter)
IteratorSize(::Type{MeanCov{T}}) where {T} = IteratorSize(T) isa Base.HasShape ?  Base.HasLength() : IteratorSize(T)

function register!(m::MeanCov{<:Channel}, x)
    put!(m.iter, x)
    iterate(m)
end
function register!(m::MeanCov{<:Channel}, x, u)
    put!(m.iter, x)
    iterate(m, u[2])
end
Base.close(m::MeanCov{<:Channel}) = close(m.iter)

@inline function iterate(mc::MeanCov)
    u = iterate(mc.iter)
    if u === nothing
        return nothing
    end
    x, state = u
    n = 0
    delta = copy(x)
    m = delta/(n+1)

    m2 = outer(delta, x - m)
    (m, 1/n, m2), (m, m2, delta, n + 1, state)
end

iterate(mc::MeanCov, mcstate) = iterate_(mc, ismutable(typeof(mcstate[1])), mcstate...)

function iterate_(mc::MeanCov, ::Val{false}, m, m2, delta, n, state )
    u = iterate(mc.iter, state)
    if u === nothing
        return nothing
    end
    x, state = u
    delta = x - m
    m  += delta/(n+1)
    m2 += outer(delta, x - m)
    (m, 1/n, m2), (m, m2, delta, n + 1, state)
end

function iterate_(mc::MeanCov, ::Val{true}, m, m2, delta, n, state )
    u = iterate(mc.iter, state)
    if u === nothing
        return nothing
    end
    x, state = u

    for i in eachindex(m)
        delta[i] = x[i] - m[i]
        m[i] += (delta[i])/(n+1)
    end
    for i in eachindex(m)
        for j in eachindex(m)
            m2[i,j] += outer(delta[i], x[j]-m[j])
        end
    end
    (m, 1/n, m2), (m, m2, delta, n + 1, state)
end

struct MeanVar{T}
    iter::T
end
import Base: length, IteratorSize
length(x::MeanVar) = length(x.iter)
IteratorSize(::Type{MeanVar{T}}) where {T} = IteratorSize(T) isa Base.HasShape ?  Base.HasLength() : IteratorSize(T)

function register!(m::MeanVar{<:Channel}, x)
    put!(m.iter, x)
    iterate(m)
end
function register!(m::MeanVar{<:Channel}, x, u)
    put!(m.iter, x)
    iterate(m, u[2])
end
Base.close(m::MeanVar{<:Channel}) = close(m.iter)

@inline function iterate(mc::MeanVar)
    u = iterate(mc.iter)
    if u === nothing
        return nothing
    end
    x, state = u
    n = 0
    delta = copy(x)
    m = delta/(n+1)

    m2 =(delta).*(x - m)
    (m, 1/n, m2), (m, m2, delta, n + 1, state)
end

iterate(mc::MeanVar, mcstate) = iterate_(mc, ismutable(typeof(mcstate[1])), mcstate...)

function iterate_(mc::MeanVar, ::Val{false}, m, m2, delta, n, state )
    u = iterate(mc.iter, state)
    if u === nothing
        return nothing
    end
    x, state = u
    delta = x - m
    m  += delta/(n+1)
    m2 += delta .* (x - m)
    (m, 1/n, m2), (m, m2, delta, n + 1, state)
end

function iterate_(mc::MeanVar, ::Val{true}, m, m2, delta, n, state )
    u = iterate(mc.iter, state)
    if u === nothing
        return nothing
    end
    x, state = u

    for i in eachindex(m)
        delta[i] = x[i] - m[i]
        m[i] += (delta[i])/(n+1)
        m2[i] += outer(delta[i], x[i] - m[i])
    end
    (m, 1/n, m2), (m, m2, delta, n + 1, state)
end

export OnlineStat
import Base: push!
import Statistics: mean, cov, std

"""

```
stats = map(OnlineStat, (x, θ, η))

map(push!, stats, (x, θ, η))

mean.(stats)
cov.(stats)
```
"""
mutable struct OnlineStat{S}
    state::S
    function OnlineStat(x)
        m = Bridge.MeanVar(Channel{typeof(x)}(1))
        u = Bridge.register!(m, x)
        new{typeof((m,u...))}((m,u...))
    end
end
function push!(o::OnlineStat, x)
    @assert o.state != nothing # cannot happen for now
    m, val, s = o.state
    u = Bridge.register!(m, x, (val, s))
    o.state = (o.state[1], u[1], u[2])
end
mean(o::OnlineStat) = o.state[2][1]
std(o::OnlineStat) = sqrt.(o.state[2][2]*o.state[2][3])
#=
S = OnlineStat(ones(5))
for i in 2:10
    push!(S, i*ones(5))
end
mean(S)
cov(S)
cov(repeat(collect(1:10)', outer=5)') == cov(S)
=#


text:
@testset "Substitution" begin

    substitution = AgentBasedModels.substitution

    abm = @agent(1,
            [l1,l2]::Local,
            [g1,g2]::Global,
            [ga1,ga2]::GlobalArray,
            [id1,id2]::Identity,
            [m1,m2]::Medium)
    s = SimulationFree(abm,box=[(:x,0,1)],medium=[MediumFlat("Dirichlet",10)])
    p = AgentBasedModels.Program_(abm,s)

    @test substitution(:(x,x),p) == :(localV[ic1_, 1],localV[ic1_, 1])
    @test substitution(:(l1,l1),p) == :(localV[ic1_, 2],localV[ic1_, 2])
    @test substitution(:(l2,l2),p) == :(localV[ic1_, 3],localV[ic1_, 3])
    @test substitution(:(g1,g1),p) == :(globalV[1],globalV[1])
    @test substitution(:(g2,g2),p) == :(globalV[2],globalV[2])
    @test substitution(:(ga1,ga1),p) == :(ga1,ga1)
    @test substitution(:(ga2,ga2),p) == :(ga2,ga2)
    @test substitution(:(id,id),p) == :(identityV[ic1_,1],identityV[ic1_,1])
    @test substitution(:(id1,id1),p) == :(identityV[ic1_,2],identityV[ic1_,2])
    @test substitution(:(id2,id2),p) == :(identityV[ic1_,3],identityV[ic1_,3])
    @test substitution(:(m1,m1),p) == :(mediumV[ic1_,1],mediumV[ic1_,1])
    @test substitution(:(m2,m2),p) == :(mediumV[ic1_,2],mediumV[ic1_,2])

    @test substitution(:(x_new,x_new),p,post="_new",addName="Copy") == :(localVCopy[ic1_, 1],localVCopy[ic1_, 1])
    @test substitution(:(l1_new,l1_new),p,post="_new",addName="Copy") == :(localVCopy[ic1_, 2],localVCopy[ic1_, 2])
    @test substitution(:(l2_new,l2_new),p,post="_new",addName="Copy") == :(localVCopy[ic1_, 3],localVCopy[ic1_, 3])
    @test substitution(:(g1_new,g1_new),p,post="_new",addName="Copy") == :(globalVCopy[1],globalVCopy[1])
    @test substitution(:(g2_new,g2_new),p,post="_new",addName="Copy") == :(globalVCopy[2],globalVCopy[2])
    @test substitution(:(ga1_new,ga1_new),p,post="_new",addName="Copy") == :(ga1Copy,ga1Copy)
    @test substitution(:(ga2_new,ga2_new),p,post="_new",addName="Copy") == :(ga2Copy,ga2Copy)
    @test substitution(:(id_new,id_new),p,post="_new",addName="Copy") == :(identityVCopy[ic1_,1],identityVCopy[ic1_,1])
    @test substitution(:(id1_new,id1_new),p,post="_new",addName="Copy") == :(identityVCopy[ic1_,2],identityVCopy[ic1_,2])
    @test substitution(:(id2_new,id2_new),p,post="_new",addName="Copy") == :(identityVCopy[ic1_,3],identityVCopy[ic1_,3])
    @test substitution(:(m1_new,m1_new),p,post="_new",addName="Copy") == :(mediumVCopy[ic1_,1],mediumVCopy[ic1_,1])
    @test substitution(:(m2_new,m2_new),p,post="_new",addName="Copy") == :(mediumVCopy[ic1_,2],mediumVCopy[ic1_,2])

    @test substitution(:(∇(x)),p, op=:∇,index=[[-1,0,0],[0,0,0],[+1,0,0]],factor=[1,1,1],opF=:(X/2)) ==
        :(((localV[ic1_ + -1, 1] + localV[ic1_, 1]) + localV[ic1_ + 1, 1]) / 2)

    @test substitution(:(∇(x)),p, op=:∇,nargs=1) == :(∇(x))
    @test substitution(:(∇(x)),p, op=:∇,nargs=1,addArgs=[:dt]) == :(∇(x,dt))
    @test substitution(:(∇(x)),p, op=:∇,nargs=1,opF=:g) == :(g(x))

end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
#
# Execution machinery
#

function price_input_vector(program::AbstractPricer, scenario::Scenario) :: Vector{Float64}
    return price_input_vector(program.input_riskfactors, scenario, program.currency)
end

# A NonCashRiskFactor has no Currency information, therefore we don't need to apply `exch`
function parse_risk_factor_value(scenario::Scenario, risk_factor::NonCashRiskFactor, ::TO_CURRENCY) :: Float64 where {TO_CURRENCY<:Currencies.Currency}
    return scenario[risk_factor]
end

function parse_risk_factor_value(scenario::Scenario, risk_factor::SpotCurrency{FROM_CURRENCY}, to::TO_CURRENCY) :: Float64 where {FROM_CURRENCY<:Currencies.Currency, TO_CURRENCY<:Currencies.Currency}
    exch(scenario, FROM_CURRENCY(), to)
end

function parse_risk_factor_value(scenario::Scenario, risk_factor::CashRiskFactor, to::Currencies.Currency) :: Float64
    exch(scenario, scenario[risk_factor], to)
end

function price_input_vector(
        input_symbols::OptimizingIR.LookupTable{RiskFactor},
        scenario::Scenario,
        functional_currency::Currencies.Currency) :: Vector{Float64}

    if isempty(input_symbols)
        # there are no risk factors
        return Vector{Float64}()
    end

    input_vector = zeros(Float64, length(input_symbols))

    for risk_factor in input_symbols
        if !isa(risk_factor, SpotCurrency)
            # For currencies, `exch` checks wether currency prices are available.
            # For any other risk factor, the Scenario must provide values for Risk Factors.
            # See `parse_risk_factor_value`.
            @assert haskey(scenario, risk_factor) "Scenario has no value for Risk Factor $risk_factor."
        end

        @inbounds input_vector[ OptimizingIR.indexof(input_symbols, risk_factor) ] = parse_risk_factor_value(scenario, risk_factor, functional_currency)
    end

    return input_vector
end

# converts `program.pricing_function`
# to something that accepts a vector as input
# so that ForwardDiff can do its thing
# Regarding `invokelatest`, see https://felipenoris.github.io/OptimizingIR.jl/dev/#World-Age-Problem-1
function price_closure(program::Pricer) :: Function
    if isa(program.pricing_function, Function)
        x -> Base.invokelatest(program.pricing_function, x)
    else
        x -> program.pricing_function(x)
    end
end

# Regarding `invokelatest`, see https://felipenoris.github.io/OptimizingIR.jl/dev/#World-Age-Problem-1
function price_closure(program::CashflowPricer) :: Function
    # the CashflowPricer returns the contract price along with
    # the prices for each cashflow.

    if isempty(program.output_index_to_cashflow_type)
        # there are no events for this contract
        if isa(program.pricing_function, Function)
            return x -> Base.invokelatest(program.pricing_function, x)
        else
            return x -> program.pricing_function(x)
        end
    else
        if isa(program.pricing_function, Function)
            return x -> Base.invokelatest(program.pricing_function, x)[program.price_output_index]
        else
            return x -> program.pricing_function(x)[program.price_output_index]
        end
    end
end

function price(program::AbstractPricer, scenario::Scenario) :: Number
    x = price_input_vector(program, scenario)
    f = price_closure(program)
    return f(x)
end

function eachcashflow(p::CashflowPricer, scenario::Scenario) :: Vector{Cashflow}

    if isempty(p.output_index_to_cashflow_type)
        return Vector{Cashflow}()
    else

        cashflows = Vector{Cashflow}(undef, length(p.output_index_to_cashflow_type))

        x = price_input_vector(p, scenario)
        result = p.pricing_function(x)

        i = 1
        for (index, cftype) in p.output_index_to_cashflow_type
            val = result[index]
            cf = Cashflow(cftype.event, cftype.maturity, val, p.currency)
            cashflows[i] = cf
            i += 1
        end
    end

    return cashflows
end

function exposures(::DeltaNormalExposuresMethod, program::AbstractPricer, scenario::Scenario) :: ExposureResult
    result = ExposureResult()
    first_order_exposures!(result, program, scenario)
    return result
end

function exposures(::DeltaGammaApproxExposuresMethod, program::AbstractPricer, scenario::Scenario) :: ExposureResult
    result = ExposureResult()
    first_order_exposures!(result, program, scenario)
    second_order_exposures!(result, program, scenario)
    return result
end

function first_order_exposures!(result::ExposureResult, program::AbstractPricer, scenario::Scenario)

    x = price_input_vector(program, scenario)

    if isempty(x)
        return result
    end

    f_grad = args -> ForwardDiff.gradient(price_closure(program), args)
    exposures_vector = x .* f_grad(x)

    # sanity check
    _risk_factors = riskfactors(program)
    @assert length(exposures_vector) == length(_risk_factors) "Size of exposures vector ($(length(exposures_vector))) should equal the size of inputs vector ($(length(_risk_factors)))."

    for i in 1:length(exposures_vector)
        # consider as zero exposure values below eps(1_000_000_000_000.0)
        if abs(exposures_vector[i]) > 0.0001220703125
            @inbounds result[_risk_factors[i]] = exposures_vector[i]
        end
    end

    nothing
end

"""
Returns second order exposures approximation as:

``\\frac{1}{2} \\frac{ \\partial^2 P }{ \\partial S^2 } S^2``

where `P` is the pricing function,
and `S` is a Risk Factor.
"""
function second_order_exposures!(result::ExposureResult, program::AbstractPricer, scenario::Scenario)

    x = price_input_vector(program, scenario)

    if isempty(x)
        return result
    end

    f_hessian = args -> ForwardDiff.hessian(price_closure(program), args)
    exposures_vector = 0.5 * (LinearAlgebra.diag(f_hessian(x)) .* (x.^2))

    # sanity check
    _risk_factors = riskfactors(program)
    @assert length(exposures_vector) == length(_risk_factors) "Size of exposures vector ($(length(exposures_vector))) should equal the size of inputs vector ($(length(_risk_factors)))."

    for i in 1:length(exposures_vector)
        # consider as zero exposure values below eps(1_000_000_000_000.0)
        if abs(exposures_vector[i]) > 0.0001220703125
            @inbounds result[SecondOrderRiskFactor(_risk_factors[i])] = exposures_vector[i]
        end
    end

    nothing
end


text:
using Nanomsg
using Base.Test

# jl_nn_errno_check("Transform")

# println("Symbol tree:", CSymbols.J_NN_ALL)

@assert(isbits(_NNSymbolProperties))

#println("Symbol tree:", CSymbols.J_NN_ALL)

addr = "tcp://127.0.0.1:7789"

pub = Socket(CSymbols.AF_SP, CSymbols.NN_BUS)
sub = Socket(CSymbols.AF_SP, CSymbols.NN_BUS)

bind(pub, addr)
connect(sub, addr)

function loop(msg, sender, recvr)
	send(sender, msg)
	if isa(msg, AbstractString)
		result = recv(recvr, AbstractString, CSymbols.NN_NO_FLAG)
		@assert(msg == result)
	else
		result = recv(recvr, CSymbols.NN_NO_FLAG)
		@assert(msg == result)
	end
end

data = ones(UInt8, 11)

bench_string() = loop("TestMessage", pub, sub)
bench_string_alt() = loop("\u2200 x \u2203 y", pub, sub)
bench_array() = loop(data, pub, sub)

bench_string()

#using Benchmark

#println(compare([bench_string, bench_string_alt, bench_array], 100))

#@time bench_string()
#@time bench_string_alt()
#@time bench_array()


text:
#Amino acid...
const AA_DICT = Dict(
    'A' => 1, 'C' => 2, 'D' => 3, 'E' => 4, 'F' => 5,
    'G' => 6, 'H' => 7, 'I' => 8, 'K' => 9, 'L' => 10,
    'M' => 11, 'N' => 12, 'P' => 13, 'Q' => 14, 'R' => 15,
    'S' => 16, 'T' => 17, 'V' => 18, 'W' => 19, 'Y' => 20,
    '*' => 21, 'X' => 22
    )

#possibly add ambig handling for aas...

#Nucleotide...
const NT_DICT = Dict('A' => 1, 'C' => 2, 'G' => 3, 'T' => 4)

#disambiguate!
const IUPACbool = Dict{Char,Array{Bool,1}}(Dict())
IUPACbool['A']=[true,false,false,false]
IUPACbool['C']=[false,true,false,false]
IUPACbool['G']=[false,false,true,false]
IUPACbool['T']=[false,false,false,true]
IUPACbool['U']=[false,false,false,true]
IUPACbool['R']=[true,false,true,false]
IUPACbool['Y']=[false,true,false,true]
IUPACbool['S']=[false,true,true,false]
IUPACbool['W']=[true,false,false,true]
IUPACbool['K']=[false,false,true,true]
IUPACbool['M']=[true,true,false,false]
IUPACbool['B']=[false,true,true,true]
IUPACbool['D']=[true,false,true,true]
IUPACbool['H']=[true,true,false,true]
IUPACbool['V']=[true,true,true,false]
IUPACbool['N']=[true,true,true,true];

function resolve_base(c)
    uc = uppercase(c)
    if uc in keys(IUPACbool)
        return ['A','C','G','T'][rand((1:4)[IUPACbool[uc]])]
    else
        return uc
    end
end

function resolve_seq(s)
    return join(resolve_base.(collect(s)))
end

#encode string
function string2encoding(seq::String, lookup_dic::Dict{Char,Int}; missing_chars = Set{Char}())
    vec = Int64[]
    for c in seq
        if !haskey(lookup_dic, c)
            push!(missing_chars, c)
            push!(vec, 0)
        else
            push!(vec, lookup_dic[c])
        end
    end
    return vec, missing_chars
end

max_stars_repo_name: alecpnkw/SeqUMAP

max_issues_repo_name: alecpnkw/SeqUMAP

max_forks_repo_name: alecpnkw/SeqUMAP


text:
@inline function safe_load(ptr::Ref{T}) where T
    dest = Ref{T}()
    dest_ptr = Base.unsafe_convert(Ptr{T}, dest)
    src_ptr = Base.unsafe_convert(Ptr{T}, ptr)
    if probe_read_kernel(dest_ptr, sizeof(T), src_ptr) == 0
        return dest[]
    else
        return nothing
    end
end

@inline elemptr(ptr::Ptr{T}, ::Val{field}) where {T,field} =
    reinterpret(Ptr{API._fieldtype(T, Val(field))}, ptr + API.offsetof(T, Val(field)))
@inline arrptr(ptr::Ptr{T}, idx) where T =
    ptr + (sizeof(T)*(idx-1))
@inline function elemptr(ptr::Cptr{T}, ::Val{field}) where {T,field}
    S = API._fieldtype(typeof(ptr), Val(field))
    if S <: Cptr
        # unsafe, need to use helper
        ref = ZeroInitRef(S)
        ref_ptr = Base.unsafe_convert(Ptr{S}, ref)
        if probe_read_kernel(ref_ptr, sizeof(S), reinterpret(Ptr{S}, ptr)) != 0
            return nothing
        end
        return ref[]
    else
        return getproperty(ptr, field)
    end
end
@inline function elemptr(ptr::Cptr{T}) where T
    # unsafe, need to use helper
    ref = ZeroInitRef(T)
    ref_ptr = Base.unsafe_convert(Ptr{T}, ref)
    if probe_read_kernel(ref_ptr, sizeof(T), reinterpret(Ptr{T}, ptr)) != 0
        return nothing
    end
    return ref[]
end
@inline arrptr(ptr::Cptr{T}, idx) where T =
    getindex(ptr, idx)
@inline elemptr(::Nothing, field) = nothing
@inline arrptr(::Nothing, field) = nothing
@inline _toptr(ptr::Ptr) = ptr
@inline _toptr(ref::Base.RefValue{T}) where T = Base.unsafe_convert(Ptr{T}, ref)
@inline _toptr(ptr::Cptr) = ptr

"Returns a pointer to the final element referenced in `ex`, where the
parent object of `ex` is a `Ptr` or `Ref`."
macro elemptr(ex)
    @assert Meta.isexpr(ex, :(.)) || Meta.isexpr(ex, :ref) "@elemptr expects a struct field or array access"
    rootex = Expr(:block, Expr(:block), nothing)
    refex = rootex
    while true
        if Meta.isexpr(ex, :(.))
            refex.args[2] = :($elemptr(nothing, $(Val(ex.args[2].value))))
        elseif Meta.isexpr(ex, :ref)
            if length(ex.args) == 2
                refex.args[2] = :($arrptr(nothing, $(esc(ex.args[2]))))
            else
                refex.args[2] = :($elemptr(nothing))
            end
        else
            break
        end
        refex = refex.args[2]
        ex = ex.args[1]
    end
    refex.args[2] = :($_toptr($(esc(ex))))
    rootex
end


text:
"""
    GadgetConfig(snapprops, subprops, particle_types)

Returns a GADGET Configuration struct with the snap and sub file properties and the mapping
from particle type to the index in the GADGET files.
"""
struct GadgetConfig
    snapprops::Dict{String}
    subprops::Dict{String}
    particle_types::Dict{Symbol,Int}
end

const _snapprops_default, _subprops_default, _particle_types_default = let
    mass = (h -> 1e10 / h.h0, u"Msun")
    pos = (h -> 1 / (h.h0 * (h.z + 1)), u"kpc")
    vel = (h -> 1 / sqrt(h.z + 1), u"km/s")

    snapprops = Dict("POS" => pos, "VEL" => vel, "MASS" => mass)

    subprops = Dict(
        "MTOT" => mass,
        "GPOS" => pos,
        "MVIR" => mass,
        "RVIR" => pos,
        "MSUB" => mass,
        "SPOS" => pos,
        "SVEL" => (h -> 1, u"km/s"), # physical units
    )

    particle_types = Dict(:gas => 0, :dm => 1, :stars => 4, :bh => 5)

    snapprops, subprops, particle_types
end

"""
    GadgetConfig(; snapprops, subprops, particle_types)

Convenience method for creating a `GadgetConfig` where the properties of the snap or sub file are
not specified.

The default values are as follows:

```julia
using Unitful, UnitfulAstro

mass = (h -> 1e10 / h.h0, u"Msun")
pos = (h -> 1 / (h.h0 * (h.z + 1)), u"kpc")
vel = (h -> 1 / sqrt(h.z + 1), u"km/s")

snapprops = Dict("POS" => pos, "VEL" => vel, "MASS" => mass)

subprops = Dict(
    "MTOT" => mass,
    "GPOS" => pos,
    "MVIR" => mass,
    "RVIR" => pos,
    "MSUB" => mass,
    "SPOS" => pos,
    "SVEL" => (h -> 1, u"km/s"), # physical units
)

particle_types = Dict(:gas => 0, :dm => 1, :stars => 4, :bh => 5)
```
"""
function GadgetConfig(;
    snapprops=_snapprops_default,
    subprops=_subprops_default,
    particle_types=_particle_types_default,
)
    GadgetConfig(snapprops, subprops, particle_types)
end

function Base.show(io::IO, ::MIME"text/plain", obj::GadgetConfig)
    print(io, "GadgetConfig")
end


text:
module TestZygote

include("preamble.jl")

using Zygote

f0 = WeightedAdd(1, 2)
f1 = WeightedAdd([1, 10], 2)

@test gradient((x, y) -> sum(f0.(x, y)), [1, 10], 100) == ([1, 1], 4)
@test gradient((x, y) -> sum(f1.(x, y)), [1, 10], 100) == ([1, 10], 4)

c0 = poly([1, 10], [2, 20])
c1 = AddCall(poly([1, 10], 0), poly(0, [2, 20]))
@test gradient(x -> sum(c0.(x)), [1, 2]) == ([2, 20],)
@test gradient(x -> sum(c1.(x)), [1, 2]) == ([2, 20],)

end  # module


text:
# Interface to NCBI Entrez utilities
# http://www.ncbi.nlm.nih.gov/books/NBK25497/
# Date: May 6, 2016
# Authors: Isabel Restrepo, Paul Stey
# BCBI - Brown University
# Version: Julia 0.4.5
using Requests
using LightXML
using HttpCommon

using  XMLconvert

# Helper function to build the url and open a handle to it
# Uses HTTP POST instead of GET for long queries
function open_entrez(cgi, params, post=false)
    if get(params, "tool", "") == ""
        params["tool"] = "BioJulia"
    end
    if get(params, "email", "") == ""
        error("Email address is required to search Entrez")
    end
    #open a handle to Entrez
    if post
        response = Requests.post(cgi; data=params)

    else
        response = Requests.get(cgi; query=params)
    end

    println("NCBI Response: ", STATUS_CODES[response.status])

    if response.status != 200
        error("Bad response from NCBI ENTREZ")
    end

    # convert binary xml data to ascii
    return ASCIIString(response.data)
end

"""
    esearch(search_dict)

Request list of UIDs matiching a query - see [NCBI Entrez:Esearch](http://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.ESearch)

###Arguments

* `search_dic::Dict`: dictionary specifying search criteria

###Output

* `::ASCIIString`: XML response from NCBI

###Example

```julia
search_dic = Dict("db"=>"pubmed", "term" => search_term,
"retstart" => 0, "retmax"=>5, "tool" =>"BioJulia",
"email" => "email")
esearch_response = esearch(search_dic)
```

###Note

* email must be a valid email address (otherwise pubmed will block you)
* search_term corresponds to a valid [PubMed Search](http://www.ncbi.nlm.nih.gov/pubmed/advanced). It may contain one or more filtering criteria using AND/OR.
For instance:

`search_term = (asthma[MeSH Terms]) AND ("2001/01/29"[Date - Publication] : "2010"[Date - Publication])`.

"""
function esearch(search_dic)
    # Seach Entrez database
    cgi = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi"
    variables = Dict("db"=>"db")
    variables = merge(variables, search_dic)
    # return variables
    return open_entrez(cgi, variables, false)
end

"""
    efetch(fetch_dic, id_list)

Retrieve data records from a list of UIDs - see [NCBI Entrez: EFetch](http://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.EFetch)

###Arguments

* `fetch_dic::Dict` - Dictionary specifying fetch criteria
* `id_list::Array` - List of ids e.g those embedded in response from esearch

###Results

* `::ASCIIString` - XML response from NCBI

###Example

```julia
# get the list of ids
if !haskey(esearch_dict, "IdList")
  error("Error: IdList not found")
end

ids = []

for id_node in esearch_dict["IdList"][1]["Id"]
  push!(ids, id_node)
end

# define the fetch dictionary
fetch_dic = Dict("db"=>"pubmed","tool" =>"BioJulia",
"email" => email, "retmode" => "xml", "rettype"=>"null")

# fetch
efetch_response = efetch(fetch_dic, ids)
```

"""
function efetch(fetch_dic, id_list)
    post = false
    cgi = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi"
    # NCBI prefers an HTTP POST instead of an HTTP GET if there are
    # more than about 200 IDs
    if length(id_list) > 200
        post = true
    end
    id_list = join(id_list,",")
    variables = Dict("id"=>id_list)
    variables = merge(variables, fetch_dic)
    return open_entrez(cgi, variables,  post)
end

"""
    elink(elink_dict)

Lists, checks or returns UIDs linked to an input list of UIDs in the same or
different Entrez database. For more info see
[NCBI Entrez:ELink](http://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.ELink)

###Arguments

* `elink_dict::Dict`: dictionary specifying elink inputs as credentials, ids...

###Output

* `::ASCIIString`: XML response from NCBI

###Example

```julia
pmid = "19304878"
elink_dict = Dict("dbfrom" =>"pubmed", "id" => pmid,
                  "linkname" => "pubmed_pubmed", "email"=>email)
elink_response = elink(elink_dict)
```
"""

function elink(elink_dict)
    cgi = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/elink.fcgi"
    return open_entrez(cgi, elink_dict)
end

"""
    esummary(esummary_dict)

Return document summaries for a list of input UIDs. For more info see
[NCBI Entrez:ESummary](http://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.ESummary)

###Arguments

* `esummary_dict::Dict`: dictionary specifying esummary inputs as credentials, ids...

###Output

* `::ASCIIString`: XML response from NCBI

###Example

```julia
pmid = "30367"
esummary_dict = Dict("db" =>"pubmed", "id" => pmid, "email"=>email)
esummary_response = esummary(esummary_dict)
```
"""
function esummary(esummary_dict)
    cgi = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi"
    return open_entrez(cgi, esummary_dict)
end

"""
    eparse(response::ASCIIString)

Converts NCBI XML response into a Julia dictionary

"""
function eparse(ncbi_response::ASCIIString)
    xdoc = parse_string(ncbi_response)
    # get the root element
    xroot = root(xdoc)  # an instance of XMLElement
    # get all child nodes and append to dictionary
    node_element = xroot
    #convert attributes to elements
    attributes_to_elements!(xroot)
    #convert to a dictionary
    dict = xml2dict(xroot)
    return dict
end

"""
    eparse_from_file(xml_file::ASCIIString)

Converts NCBI XML (previously saved) file into a Julia dictionary

"""
function eparse_from_file(xml_file::ASCIIString)
    xdoc = parse_file(xml_file)
    # get the root element
    xroot = root(xdoc)  # an instance of XMLElement
    # get all child nodes and append to dictionary
    node_element = xroot
    #convert attributes to elements
    attributes_to_elements!(xroot)
    #convert to a dictionary
    dict = xml2dict(xroot)
    return dict
end

# """
#     save_efetch(efetch_dict, db_path)
#
# Save the results (dictionary) of an entrez fetch to a SQLite database.
#
# ####Note:
#
# It is best to assure the databese file does not exist. If the file
# path corresponds to an exixting database, the system
# attempts to use that database, which must contain the correct tables.
#
# ###Example
#
# ```julia
# #SQLite example
# db_backend("SQLite")
# db_config = Dict(:db_path=>"test_db.slqite", :overwrite=>true)
# db = save_efetch(efetch_dict, db_config)
# ```
#
# """
# function save_efetch(efetch_dict, db_config)
#     #init database with its structure only if file doesn't exist
#     db = DB.init_database(db_config)
#
#     if !haskey(efetch_dict, "PubmedArticle")
#         println("Error: Could not save to DB key:PubmedArticleSet not found")
#         return
#     end
#     articles = efetch_dict["PubmedArticle"]
#
#     #articles should be an array
#     if !isa(articles, Array{Any, 1})
#         println("Error: Could not save to DB articles should be in an Array")
#         return
#     end
#
#     println("Saving " , length(articles) ,  " articles to database")
#
#     for article in articles
#
#         if !haskey(article,"MedlineCitation")
#             println("Error: Could not save to DB key:MedlineCitation not found")
#             return
#         end
#
#         pmid = DB.NULL[_db_backend[1]];
#         title = DB.NULL[_db_backend[1]];
#         pubYear = DB.NULL[_db_backend[1]];
#
#
#         # PMID is used as primary key - therefore it must be present
#         if haskey(article["MedlineCitation"][1],"PMID")
#             pmid = article["MedlineCitation"][1]["PMID"][1]["PMID"][1]
#         else
#             println("Error: Could not save to DB key:PMID not found - cannot be NULL")
#             return
#         end
#
#         # Retrieve basic article info
#         if haskey(article["MedlineCitation"][1],"Article")
#             if haskey(article["MedlineCitation"][1]["Article"][1], "ArticleTitle")
#                 title = article["MedlineCitation"][1]["Article"][1]["ArticleTitle"][1]
#             end
#             if haskey(article["MedlineCitation"][1]["Article"][1], "ArticleDate")
#                 if haskey(article["MedlineCitation"][1]["Article"][1]["ArticleDate"][1], "Year")
#                     pubYear = article["MedlineCitation"][1]["Article"][1]["ArticleDate"][1]["Year"][1]
#                 end
#             else  #series of attempts to pull a publication year from alternative xml elements
#                 try
#                     pubYear = article["MedlineCitation"][1]["Article"][1]["Journal"][1]["JournalIssue"][1]["PubDate"][1]["Year"][1]
#                 catch
#                     try
#                         pubYear = article["MedlineCitation"][1]["Article"][1]["Journal"][1]["JournalIssue"][1]["PubDate"][1]["MedlineDate"][1]
#                         pubYear = parse(Int64, pubYear[1:4])
#                     catch
#                         println("Warning: No date found")
#                     end
#                 end
#             end
#
#             # Save article data
#             DB.insert_row(db, "article", Dict(:pmid => pmid,
#             :title=>title,
#             :pubYear=>pubYear))
#
#             # insert all authors
#             forename = DB.NULL[_db_backend[1]]
#             lastname = DB.NULL[_db_backend[1]]
#             if haskey(article["MedlineCitation"][1]["Article"][1], "AuthorList")
#                 authors = article["MedlineCitation"][1]["Article"][1]["AuthorList"][1]["Author"]
#                 for author in authors
#
#                     if author["ValidYN"][1] == "N"
#                         continue
#                     end
#
#                     if haskey(author, "ForeName")
#                         forename = author["ForeName"][1]
#                     else
#                         forname = "UNKNOWN"
#                     end
#
#                     if haskey(author, "LastName")
#                         lastname = author["LastName"][1]
#                     else
#                         println("Skipping Author: ", author)
#                         continue
#                     end
#
#                     # Authors must be unique - insert only if it doesn't exist
#                     # DB.exists(db, "author", )
#
#                     # Save author data
#                     author_id = DB.insert_row(db, "author",
#                     Dict(:id => DB.NULL[_db_backend[1]],
#                     :forename => forename,
#                     :lastname => lastname))
#
#                     if (author_id >= 0 )
#                         DB.insert_row(db, "author2article",
#                         Dict(:aid =>author_id, :pmid => pmid))
#                     end
#
#                 end
#             end
#
#             # Save related "keywords" of MESH Descriptors
#             if haskey(article["MedlineCitation"][1], "MeshHeadingList")
#                 if haskey(article["MedlineCitation"][1]["MeshHeadingList"][1], "MeshHeading")
#                     mesh_headings = article["MedlineCitation"][1]["MeshHeadingList"][1]["MeshHeading"]
#                     for heading in mesh_headings
#
#                         if !haskey(heading,"DescriptorName")
#                             println("Error: MeshHeading must have DescriptorName")
#                             return
#                         end
#
#                         #save descriptor
#                         descriptor_name = heading["DescriptorName"][1]["DescriptorName"][1]
#                         descriptor_name = normalize_string(descriptor_name, casefold=true)
#
#                         did = heading["DescriptorName"][1]["UI"][1]
#                         did_int = parse(Int64, did[2:end])  #remove preceding D
#
#                         DB.insert_row(db, "mesh_descriptor",
#                         Dict(:id=>did_int, :name=>descriptor_name))
#
#                         heading["DescriptorName"][1]["MajorTopicYN"][1] == "Y" ? dmjr = 1 : dmjr = 0
#
#                         #save the qualifiers
#                         if haskey(heading,"QualifierName")
#                             qualifiers = heading["QualifierName"]
#                             for qual in qualifiers
#                                 qualifier_name = qual["QualifierName"][1]
#                                 qualifier_name = normalize_string(qualifier_name, casefold=true)
#
#                                 qid = qual["UI"][1]
#                                 qid_int = parse(Int64, qid[2:end])  #remove preceding Q
#
#                                 DB.insert_row(db, "mesh_qualifier",
#                                 Dict(:id=>qid_int, :name=>qualifier_name) )
#
#                                 qual["MajorTopicYN"][1] == "Y" ? qmjr = 1 : qmjr = 0
#
#                                 #save the heading related to this paper
#                                 DB.insert_row(db, "mesh_heading",
#                                 Dict(:id=>DB.NULL[_db_backend[1]], :pmid=> pmid, :did=>did_int,
#                                 :qid=>qid_int, :dmjr=>dmjr, :qmjr=>qmjr) )
#                             end
#                         else
#                             #save the heading related to this paper
#                             DB.insert_row(db, "mesh_heading",
#                             Dict(:id=>DB.NULL[_db_backend[1]], :pmid=> pmid, :did=>did_int,
#                             :qid=>DB.NULL[_db_backend[1]], :dmjr=>dmjr, :qmjr=>DB.NULL[_db_backend[1]]) )
#                         end
#                     end
#                 end
#             end
#
#         end
#
#     end
#
#     return db
#
# end


text:
using DSGE, Test, Dates
using JLD2, OrderedCollections, HDF5, Nullables

path = dirname(@__FILE__)
save_output = true

# Initialize model object
m = AnSchorfheide(testing = true)
dir = joinpath(saveroot(m), "output_data", "an_schorfheide", "ss0")

output_vars = [:histstates, :histobs, :histpseudo, :histshocks,
                   :forecastobs, :shockdecobs, :dettrendobs, :trendobs, :irfobs]

@testset "Test main forecast I/O functions and metadata writing" begin
    # get_forecast_input_file
    @test get_forecast_input_file(m, :init) == ""
    @test get_forecast_input_file(m, :init_draw_shocks) == ""
    @test get_forecast_input_file(m, :prior) == ""
    @test get_forecast_input_file(m, :mode) == rawpath(m, "estimate", "paramsmode.h5")
    @test get_forecast_input_file(m, :mode_draw_shocks) == rawpath(m, "estimate", "paramsmode.h5")
    m <= Setting(:sampling_method, :invalid_method)
    @test_throws ErrorException get_forecast_input_file(m, :full) == rawpath(m, "estimate", "smcsave.h5")
    m <= Setting(:sampling_method, :SMC)
    @test get_forecast_input_file(m, :full) == rawpath(m, "estimate", "smcsave.h5")
    m <= Setting(:sampling_method, :MH)
    @test get_forecast_input_file(m, :full) == rawpath(m, "estimate", "mhsave.h5")
    global overrides = forecast_input_file_overrides(m)
    overrides[:full] = "invalidpath"
    @test_throws ErrorException get_forecast_input_file(m, :full)
    overrides[:mode] = tempname()
    @test_throws ErrorException get_forecast_input_file(m, :mode)
    overrides[:full] = mktemp()[1]
    @test get_forecast_input_file(m, :full) == overrides[:full]
    @test get_forecast_input_file(m, :subset) == overrides[:full]
    overrides[:subset] = mktemp()[1]
    @test get_forecast_input_file(m, :subset) == overrides[:subset]

    # get_forecast_filename
    @test get_forecast_filename(m, :mode, :none, :histstates) == joinpath(dir, "forecast", "raw", "histstates_cond=none_para=mode_test.jld2")
    @test get_forecast_filename(m, :mode, :none, :histstates; pathfcn = workpath) == joinpath(dir, "forecast", "work", "histstates_cond=none_para=mode_test.jld2")
    @test get_forecast_filename(m, :mode, :none, :histstates; forecast_string = "test") == joinpath(dir, "forecast", "raw", "histstates_cond=none_fcid=test_para=mode_test.jld2")
    @test get_forecast_filename(m, :mode, :none, :histstates; fileformat = :h5) == joinpath(dir, "forecast", "raw", "histstates_cond=none_para=mode_test.h5")

    # get_forecast_filestring_addl
    @test DSGE.get_forecast_filestring_addl(:mode, :none) == ["para=mode", "cond=none"]
    @test DSGE.get_forecast_filestring_addl(:mode, :none; forecast_string = "test") == ["para=mode", "cond=none", "fcid=test"]
    @test_throws ErrorException DSGE.get_forecast_filestring_addl(:subset, :none)

    # get_forecast_output_files
    dict = Dict{Symbol, String}()

    for var in output_vars
        dict[var] = get_forecast_filename(m, :mode, :none, var)
    end
    @test get_forecast_output_files(m, :mode, :none, output_vars) == dict

    # write_forecast_outputs
    if haskey(ENV, "FRED_API_KEY") || isfile(joinpath(ENV["HOME"],".freddatarc"))
        global df = load_data(m, verbose = :none)
        overrides = forecast_input_file_overrides(m)
        overrides[:full] = "$(path)/../reference/smcsave_.h5"
        m <= Setting(:sampling_method, :SMC)
        m <= Setting(:forecast_block_size, 100)
        block_inds, block_inds_thin = DSGE.forecast_block_inds(m, :full)
        forecast_output_files = DSGE.get_forecast_output_files(m, :full, :none, [:histstates, :histobs, :forecaststates, :forecastobs])
        # Test that write_forecast_outputs and combine_raw_forecast_output_and_metadata run without deprecration
        for block in 1:20
            forecast_output = load("$(path)/../reference/write_forecast_outputs_$(block)_in.jld2", "forecast_output")
            DSGE.write_forecast_outputs(m, :full, [:histstates, :histobs, :forecaststates, :forecastobs],forecast_output_files, forecast_output, df = df, block_number = Nullables.Nullable(block), block_inds = block_inds_thin[block])
        end
        DSGE.combine_raw_forecast_output_and_metadata(m, forecast_output_files, verbose = :none)

        # write_forecast_metadata
        for var in output_vars
            JLD2.jldopen(dict[var], "w") do file
                DSGE.write_forecast_metadata(m, file, var)
            end
        end
        JLD2.jldopen(dict[:histstates], "r") do file
            dates = read(file, "date_indices")
            @test dates[date_mainsample_start(m)] == 1
            @test dates[date_mainsample_end(m)] == length(dates)
            @test read(file, "state_indices") == merge(m.endogenous_states, m.endogenous_states_augmented)
            @test all(x -> x == Symbol("identity"), values(read(file, "state_revtransforms")))
        end
        JLD2.jldopen(dict[:histobs], "r") do file
            @test read(file, "observable_indices") == m.observables
            @test !all(x -> x == Symbol("identity"), values(read(file, "observable_revtransforms")))
        end
        JLD2.jldopen(dict[:histpseudo], "r") do file
            @test read(file, "pseudoobservable_indices") == m.pseudo_observables
            @test !all(x -> x == Symbol("identity"), values(read(file, "pseudoobservable_revtransforms")))
        end
        JLD2.jldopen(dict[:histshocks], "r") do file
            @test read(file, "shock_indices") == m.exogenous_shocks
            @test all(x -> x == Symbol("identity"), values(read(file, "shock_revtransforms")))
        end
        JLD2.jldopen(dict[:histshocks], "r") do file
            @test read(file, "shock_indices") == m.exogenous_shocks
            @test all(x -> x == Symbol("identity"), values(read(file, "shock_revtransforms")))
        end
        JLD2.jldopen(dict[:forecastobs], "r") do file
            dates = read(file, "date_indices")
            @test dates[date_forecast_start(m)] == 1
            @test dates[date_forecast_end(m)] == length(dates)
        end
        for class in [:shockdec, :dettrend, :trend]
            JLD2.jldopen(dict[Symbol(class, "obs")], "r") do file
                dates = read(file, "date_indices")
                @test dates[date_mainsample_start(m)] == 1
                @test dates[date_forecast_end(m)] == length(dates)
            end
        end
        JLD2.jldopen(dict[:irfobs], "r") do file
            @test !haskey(file, "date_indices")
            @test all(x -> x == Symbol("identity"), values(read(file, "observable_revtransforms")))
        end
    else
        @warn "Skipping some tests because FRED_API_KEY not present"
    end
end


text:
push!(LOAD_PATH, "../LF_common/")
using m_LF3d

include("../LF_common/sparse_LF3d.jl")
include("build_nabla2.jl")

function test_main( NN::Array{Int64} )

  println("NN = ", NN)

  Nx = NN[1]
  Ny = NN[2]
  Nz = NN[3]
  Npoints = Nx*Ny*Nz

  AA = [0.0, 0.0, 0.0]
  BB = [16.0, 16.0, 16.0]

  # Initialize LF
  LF = init_LF3d_c( NN, AA, BB, verbose=true )

  @time nabla2 = build_nabla2( LF )

  @time nabla2_v2 = get_Laplacian3d_kron( LF )

  #println(nabla2)
  #println(nabla2_v2)

  #println( full(nabla2) - full(nabla2_v2) )

  #println( nabla2 - nabla2_v2 )

  srand(123)
  x = rand(Npoints)

  # Compare multiplication result
  Lx = nabla2 * x
  Lx_v2 = nabla2_v2 * x

  println( sum(Lx-Lx_v2) )

end

test_main( [30, 30, 30] )
#test_main( [3,2,4] )


text:
# import Glob

type PlotAnimation
  t::Int
  widget::Widget
  dir::ASCIIString
end

function animation(widget::Widget, dir::ASCIIString)

  # remove existing pngs
  dir = string(dir, "/")
  #map(rm, Glob.glob("qtanim*.png", dir))

  PlotAnimation(1, widget, dir)
end

function saveframe(animation::PlotAnimation)
  savepng(animation.widget, string(animation.dir, @sprintf("qtanim%04d.png", animation.t)))
  animation.t += 1
end

makegif(animation::PlotAnimation) = run(`ffmpeg -framerate 20 -i $(animation.dir)/qtanim%04d.png -y $(animation.dir)/out.gif`)


text:
"""
    combine(uvals::Vector{AbstractUncertainValue}; n = 10000*length(uvals),
        bw::Union{Nothing, Real} = nothing)

Combine multiple uncertain values into a single uncertain value. This is
done by resampling each uncertain value in `uvals`, `n` times  each,
then pooling these draws together. Finally, a kernel density estimate to the final
distribution is computed over those draws.

The KDE bandwidth is controlled by `bw`. By default, `bw = nothing`; in this case,
the bandwidth is determined using the `KernelDensity.default_bandwidth` function.

!!! tip

    For very wide, close-to-normal distributions, the default bandwidth may work well.
    If you're combining very peaked distributions or discrete populations, however,
    you may want to lower the bandwidth significantly.

# Example

```julia
v1 = UncertainValue(Normal, 1, 0.3)
v2 = UncertainValue(Normal, 0.8, 0.4)
v3 = UncertainValue([rand() for i = 1:3], [0.3, 0.3, 0.4])
v4 = UncertainValue(Normal, 3.7, 0.8)
uvals = [v1, v2, v3, v4];

combine(uvals)
combine(uvals, n = 20000) # adjust number of total draws
```
"""
function combine(uvals::Vector{AbstractUncertainValue}; n = 10000*length(uvals),
        bw::Union{Nothing, Real} = nothing)
    N = length(uvals)
    draws = zeros(Float64, N*n)

    for (i, uval) in enumerate(uvals)
        draws[(i-1)*n+1:i*n] = resample(uval, n)
    end

    return UncertainValue(UnivariateKDE, draws,
        bandwidth = bw isa Real ? bw : default_bandwidth(draws))
end

"""
    combine(uvals::Vector{AbstractUncertainValue}, weights::ProbabilityWeights;
        n = 10000*length(uvals),
        bw::Union{Nothing, Real} = nothing)

Combine multiple uncertain values into a single uncertain value. This is
done by resampling each uncertain value in `uvals` proportionally to the provided
relative analytic `weights` indicating their relative importance (these are normalised by
default, so don't need to sum to 1), then pooling these draws together. Finally, a kernel
density estimate to the final distribution is computed over the `n` total draws.

Providing `ProbabilityWeights` leads to the exact same behaviour as for `AnalyticWeights`,
but may be more appropriote when, for example, weights have been determined
quantitatively.

The KDE bandwidth is controlled by `bw`. By default, `bw = nothing`; in this case,
the bandwidth is determined using the `KernelDensity.default_bandwidth` function.

!!! tip

    For very wide, close-to-normal distributions, the default bandwidth may work well.
    If you're combining very peaked distributions or discrete populations, however,
    you may want to lower the bandwidth significantly.

# Example

```julia
v1 = UncertainValue(Normal, 1, 0.3)
v2 = UncertainValue(Normal, 0.8, 0.4)
v3 = UncertainValue([rand() for i = 1:3], [0.3, 0.3, 0.4])
v4 = UncertainValue(Normal, 3.7, 0.8)
uvals = [v1, v2, v3, v4];

# Two difference syntax options
combine(uvals, ProbabilityWeights([0.2, 0.1, 0.3, 0.2]))
combine(uvals, pweights([0.2, 0.1, 0.3, 0.2]), n = 20000) # adjust number of total draws
```

"""
function combine(uvals::Vector{AbstractUncertainValue}, weights::ProbabilityWeights;
        n = 10000*length(uvals),
        bw::Union{Nothing, Real} = nothing)

    length(uvals) == length(weights) ? "Number of values != number of weights" : nothing

    # Scale the number of draws of each value according to their relative probability
    L = length(uvals)
    wts = weights ./ weights.sum
    Ns = [ceil(Int, n*wts[i]) for i = 1:L]

    N = sum(Ns)

    draws = zeros(Float64, N)
    for (i, uval) in enumerate(uvals)
        if i == 1
            draws[1:sum(Ns[1:i])] = resample(uval, Ns[i])
        else
            draws[sum(Ns[1:(i-1)])+1:sum(Ns[1:i])] = resample(uval, Ns[i])
        end
    end

    return UncertainValue(UnivariateKDE, draws,
        bandwidth = bw isa Real ? bw : default_bandwidth(draws))
end

"""
    combine(uvals::Vector{AbstractUncertainValue}, weights::AnalyticWeights;
        n = 10000*length(uvals),
        bw::Union{Nothing, Real} = nothing)

Combine multiple uncertain values into a single uncertain value. This is
done by resampling each uncertain value in `uvals` proportionally to the provided
relative probability `weights` (these are normalised by default, so don't need
to sum to 1), then pooling these draws together. Finally, a kernel density
estimate to the final distribution is computed over the `n` total draws.

Providing `AnalyticWeights` leads to the exact same behaviour as for `ProbabilityWeights`,
but may be more appropriote when relative importance weights are assigned subjectively,
and not based on quantitative evidence.

The KDE bandwidth is controlled by `bw`. By default, `bw = nothing`; in this case,
the bandwidth is determined using the `KernelDensity.default_bandwidth` function.

!!! tip

    For very wide, close-to-normal distributions, the default bandwidth may work well.
    If you're combining very peaked distributions or discrete populations, however,
    you may want to lower the bandwidth significantly.

# Example

```julia
v1 = UncertainValue(Normal, 1, 0.3)
v2 = UncertainValue(Normal, 0.8, 0.4)
v3 = UncertainValue([rand() for i = 1:3], [0.3, 0.3, 0.4])
v4 = UncertainValue(Normal, 3.7, 0.8)
uvals = [v1, v2, v3, v4];

# Two difference syntax options
combine(uvals, AnalyticWeights([0.2, 0.1, 0.3, 0.2]))
combine(uvals, aweights([0.2, 0.1, 0.3, 0.2]), n = 20000) # adjust number of total draws
```
"""
function combine(uvals::Vector{AbstractUncertainValue}, weights::AnalyticWeights;
    n = 10000*length(uvals),
    bw::Union{Nothing, Real} = nothing)

    length(uvals) == length(weights) ? "Number of values != number of weights" : nothing

    # Scale the number of draws of each value according to their relative probability
    L = length(uvals)
    wts = weights ./ weights.sum
    Ns = [ceil(Int, n*wts[i]) for i = 1:L]

    N = sum(Ns)

    draws = zeros(Float64, N)
    for (i, uval) in enumerate(uvals)
        if i == 1
            draws[1:sum(Ns[1:i])] = resample(uval, Ns[i])
        else
            draws[sum(Ns[1:(i-1)])+1:sum(Ns[1:i])] = resample(uval, Ns[i])
        end
    end

    return UncertainValue(UnivariateKDE, draws,
        bandwidth = bw isa Real ? bw : default_bandwidth(draws))
end

"""
    combine(uvals::Vector{AbstractUncertainValue}, weights::Weights;
        n = 10000*length(uvals),
        bw::Union{Nothing, Real} = nothing)

Combine multiple uncertain values into a single uncertain value. This is
done by resampling each uncertain value in `uvals` proportionally to the provided `weights`
(these are normalised by default, so don't need to sum to 1), then pooling these draws
together. Finally, a kernel density estimate to the final distribution is computed over
the `n` total draws.

Providing `Weights` leads to the exact same behaviour as for `ProbabilityWeights` and
`AnalyticalWeights`.

The KDE bandwidth is controlled by `bw`. By default, `bw = nothing`; in this case,
the bandwidth is determined using the `KernelDensity.default_bandwidth` function.

!!! tip

    For very wide, close-to-normal distributions, the default bandwidth may work well.
    If you're combining very peaked distributions or discrete populations, however,
    you may want to lower the bandwidth significantly.

# Example

```julia
v1 = UncertainValue(Normal, 1, 0.3)
v2 = UncertainValue(Normal, 0.8, 0.4)
v3 = UncertainValue([rand() for i = 1:3], [0.3, 0.3, 0.4])
v4 = UncertainValue(Normal, 3.7, 0.8)
uvals = [v1, v2, v3, v4];

# Two difference syntax options
combine(uvals, Weights([0.2, 0.1, 0.3, 0.2]))
combine(uvals, weights([0.2, 0.1, 0.3, 0.2]), n = 20000) # adjust number of total draws
```
"""
function combine(uvals::Vector{AbstractUncertainValue}, weights::Weights;
        n = 10000*length(uvals),
        bw::Union{Nothing, Real} = nothing)

    length(uvals) == length(weights) ? "Number of values != number of weights" : nothing

    # Scale the number of draws of each value according to their relative probability
    L = length(uvals)
    wts = weights ./ weights.sum
    Ns = [ceil(Int, n*wts[i]) for i = 1:L]

    N = sum(Ns)

    draws = zeros(Float64, N)
    for (i, uval) in enumerate(uvals)
        if i == 1
            draws[1:sum(Ns[1:i])] = resample(uval, Ns[i])
        else
            draws[sum(Ns[1:(i-1)])+1:sum(Ns[1:i])] = resample(uval, Ns[i])
        end
    end

    return UncertainValue(UnivariateKDE, draws,
        bandwidth = bw isa Real ? bw : default_bandwidth(draws))

end

"""
    combine(uvals::Vector{AbstractUncertainValue}, weights::FrequencyWeights;
        bw::Union{Nothing, Real} = nothing)

Combine multiple uncertain values into a single uncertain value. This is
done by resampling each uncertain value in `uvals` according to their relative
frequencies (the absolute number of draws provided by `weights`). Finally, a kernel density
estimate to the final distribution is computed over the `sum(weights)` total draws.

The KDE bandwidth is controlled by `bw`. By default, `bw = nothing`; in this case,
the bandwidth is determined using the `KernelDensity.default_bandwidth` function.

!!! tip

    For very wide and close-to-normal distributions, the default bandwidth may work well.
    If you're combining very peaked distributions or discrete populations, however,
    you may want to lower the bandwidth significantly.

# Example

```julia
v1 = UncertainValue(Normal, 1, 0.3)
v2 = UncertainValue(Normal, 0.8, 0.4)
v3 = UncertainValue([rand() for i = 1:3], [0.3, 0.3, 0.4])
v4 = UncertainValue(Normal, 3.7, 0.8)
uvals = [v1, v2, v3, v4];

# Two difference syntax options
combine(uvals, FrequencyWeights([100, 500, 343, 7000]))
combine(uvals, pweights([1410, 550, 223, 801]))
```
"""
function combine(uvals::Vector{AbstractUncertainValue}, weights::FrequencyWeights;
        bw::Union{Nothing, Real} = nothing)

    length(uvals) == length(weights) ? "Number of values != number of weights" : nothing

    # Scale the number of draws of each value according to their relative probability
    Ns = weights
    N = sum(Ns)

    draws = zeros(Float64, N)
    for (i, uval) in enumerate(uvals)
        if i == 1
            draws[1:sum(Ns[1:i])] = resample(uval, Ns[i])
        else
            draws[sum(Ns[1:(i-1)])+1:sum(Ns[1:i])] = resample(uval, Ns[i])
        end
    end

    return UncertainValue(UnivariateKDE, draws,
        bandwidth = bw isa Real ? bw : default_bandwidth(draws))
end

export combine,
    ProbabilityWeights, pweights,
    Weights, weights,
    AnalyticWeights, aweights,
    FrequencyWeights, fweights


text:
# This file contains additional helper functions for computing the steady state.

# Discretize a state space dimension of I gridpoints
# Where agridparam is the "bending" coefficient of the grid
# i.e. grid_param = 1 implies a uniform grid
# And grid_min/grid_max are the start and endpoints of the grid
@inline function construct_asset_grid(I::Int64, grid_param::Int64, grid_min::Float64, grid_max::Float64)
    a  = collect(range(0, stop=1, length=I))
    a  = a .^ (1 / grid_param)
    a  = grid_min .+ (grid_max - grid_min) * a
    return a
end

# Provide the income and asset grids
@inline function initialize_diff_grids(a::Vector{Float64}, I::Int64, J::Int64)
    daf    = similar(a) # forward difference for a
    dab    = similar(a) # backward difference for a
    adelta = similar(a) # size of differences

    for i=1:I
        # Create a grid of lengths of overlapping intervals in a dimension.
        # The purpose is generally to compute Riemann integrals
        # by taking midpoint Riemann sums and dividing by two to adjust for
        # or average out the added Lebesgue measure given by using overlapping intervals.
        daf[i] = (i==I) ? a[end] - a[end-1] : a[i+1] - a[i]
        dab[i] = (i==1) ? a[2] - a[1]       : a[i] - a[i-1]

        # Ex:    adelta[2] =    (a_3 - a_1) / 2; if we integrate f, then
        # f[2] * adelta[2] = f_2(a_3 - a_1) / 2
        if i==1
            adelta[1]   = 0.5 * daf[1]
        elseif i==I
            adelta[end] = 0.5 * daf[end-1]
        else
            adelta[i]   = 0.5 * (daf[i-1] + daf[i])
        end
    end
    azdelta = repeat(adelta, J)

    return daf, dab, azdelta
end

# P is the Markov transition matrix
# n_income_states is the number of income states (the number of discrete states in the distribution)
# iter_num is the number of iterations willing to be accepted for convergence
@inline function compute_stationary_income_distribution(P::Matrix{Float64}, n_income_states::Int64;
                                                iter_num::Int64 = 50)
    Pt = P'
    g_z = fill(1/n_income_states, n_income_states)
    for n = 1:iter_num
        g_z_new = (speye(n_income_states) - Pt * 1000)\g_z
        diff    = maximum(abs.(g_z_new - g_z))
        if diff < 1e-5
            break
        end
        g_z = g_z_new
    end
    return g_z
end

# initial_ygrid is the unscaled income values for each income state
# income_distr is the stationary income distribution
# meanlabeff is the mean labor efficiency value that scales z such that
# The expected value of income with respect to income_distr is meanlabeff
# n_gridpoints should be the same as the number of grid points in the asset process
@inline function construct_labor_income_grid(initial_ygrid::Vector{Float64},
                                     income_distr::Vector{Float64},
                                     meanlabeff::Float64, n_gridpoints::Int64)
    z       = exp.(initial_ygrid)
    z_bar   = dot(z, income_distr)
    z  	    = (meanlabeff/z_bar) .* z
    zz      = ones(n_gridpoints, 1) * z'
    return zz
end

@inline function construct_household_problem_functions(V::Matrix{S}, w::T,
                                                       coefrra::R, frisch::R,
                                                       labtax::R,
                                                       labdisutil::R) where {R<:AbstractFloat,T<:Real,S<:Number}

    @inline function util(c::U, h::U) where {U<:Number}
        f(x::U) = coefrra == 1.0 ? log(x) : x^(1-coefrra) / (1-coefrra)
        return f(c) - labdisutil * (h ^ (1 + 1/frisch)/(1 + 1/frisch))
    end

    @inline income(h::U, z::Float64, profshare::V, lumptransfer::V, r::T,
                   a::Float64) where {T<:Number,U<:Number,V<:Number} = h * z * w * (1 - labtax) + lumptransfer + profshare + r * a

    @inline labor(z::U, val::V) where {U<:Number,V<:Number} = (z * w * (1 - labtax) * val / labdisutil) ^ frisch

    return util, income, labor
end

@inline function construct_initial_diff_matrices(V::Matrix{T},
                                                 Vaf::Matrix{T}, Vab::Matrix{T},
                                                 income::Function, labor::Function,
                                                 h::Matrix{U}, h0::Matrix{U},
                                                 zz::Matrix{S}, profshare::Matrix{T},
                                                 lumptransfer::T,
                                                 amax::S, amin::S, coefrra::S, r::S,
                                                 daf::Vector{S}, dab::Vector{S},
                                                 maxhours::S) where {S<:Number,T<:Number,U<:Number}
    I,J = size(V)
    cf  = similar(V)
    hf  = similar(V)
    cb  = similar(V)
    hb  = similar(V)

    for j=1:J, i=1:I
        if i==I
            Vaf[end, j] = income(h[end, j], zz[end, j], profshare[end, j], lumptransfer, r, amax) ^ (-coefrra)
            Vab[i, j]   = (V[i, j] - V[i-1, j]) / dab[i]
        elseif i==1
            Vaf[i, j]   = (V[i+1, j] - V[i, j]) / daf[i]
            Vab[1, j]   = income(h0[1, j], zz[1, j], profshare[1, j], lumptransfer, r, amin) ^ (-coefrra)
        else
            Vaf[i, j]   = (V[i+1, j] - V[i, j]) / daf[i]
            Vab[i, j]   = (V[i, j] - V[i-1, j]) / dab[i]
        end
        cf[i,j] = Vaf[i,j] ^ (-1 / coefrra)
        cb[i,j] = Vab[i,j] ^ (-1 / coefrra)

        hf[i,j] = min(norm(labor(zz[i,j], Vaf[i,j])), maxhours)
        hb[i,j] = min(norm(labor(zz[i,j], Vab[i,j])), maxhours)
    end

    return Vaf, Vab, cf, hf, cb, hb
end

# For initialization
@inline function calculate_ss_equil_vars(zz::Matrix{Float64}, m_ss::Float64, meanlabeff::Float64,
                                 lumptransferpc::Float64, govbondtarget::Float64)

    N_ss         = complex(1/3) # steady state hours: so that quarterly GDP = 1 in s.s
    Y_ss         = complex(1.)
    B_ss         = govbondtarget * Y_ss
    profit_ss    = complex((1 - m_ss) * Y_ss)
    profshare    = zz / meanlabeff * profit_ss
    lumptransfer = complex(lumptransferpc * Y_ss)

    return N_ss, Y_ss, B_ss, profit_ss, profshare, lumptransfer
end

@inline function calculate_ss_equil_vars(zz::Matrix{Float64},
                                 h::Matrix{ComplexF64},
                                 g::Vector{ComplexF64}, azdelta::Vector{Float64},
                                 aa::Matrix{Float64}, m_ss::Float64,
                                 meanlabeff::Float64, lumptransferpc::Float64,
                                 govbondtarget::Float64)

    # equilibrium objects
    Y_ss = N_ss  = sum(vec(zz) .* vec(h) .* g .* azdelta)
    Y_ss         = N_ss
    B_ss         = sum(g .* vec(aa) .* azdelta)
    profit_ss    = (1 - m_ss) * Y_ss
    profshare    = zz / meanlabeff * profit_ss
    lumptransfer = lumptransferpc * Y_ss
    bond_err     = B_ss / Y_ss - govbondtarget

    return N_ss, Y_ss, B_ss, profit_ss, profshare, lumptransfer, bond_err
end

@inline function hours_iteration(income::Function, labor::Function,
                                 zz::Matrix{Float64},
                                 profshare::Matrix{T},
                                 lumptransfer::T,
                                 aa::Matrix{Float64},
                                 coefrra::Float64, r::S,
                                 cf::Matrix{T}, hf::Matrix{T},
                                 cb::Matrix{T}, hb::Matrix{T},
                                 c0::Matrix{T}, h0::Matrix{U},
                                 maxhours::Float64,
                                 niter_hours::Int64) where {S<:Number,T<:Number,U<:Number}
    I, J = size(zz)
    for ih = 1:niter_hours
        for j=1:J, i=1:I
            if i==I
                cf[end, j] = income(hf[end, j], zz[end, j], profshare[end, j], lumptransfer, r, aa[end, j])
                hf[end, j] = labor(zz[end, j], cf[end, j] ^ (-coefrra))
                hf[end, j] = min(norm(hf[end, j]), maxhours)
            elseif i==1
                cb[1, j] = income(hb[1, j], zz[1, j], profshare[1, j], lumptransfer, r, aa[1, j])
                hb[1, j] = labor(zz[1, j], cb[1, j] ^ (-coefrra))
                hb[1, j] = min(norm(hb[1, j]), maxhours)
            end
            c0[i, j] = income(h0[i, j], zz[i, j], profshare[i, j], lumptransfer, r, aa[i, j])
            h0[i, j] = labor(zz[i, j], c0[i, j]^(-coefrra))
            h0[i, j] = min(norm(h0[i, j]), maxhours)
        end
    end
    return cf, hf, cb, hb, c0, h0
end

# For compute_steady_state; choose upwinding direction
@inline function upwind(ρ::Float64, V::Matrix{T}, args...;
                        Δ_HJB::Float64 = 1e6) where {T<:Number}
    A, u, h, c, s = upwind(args...)
    I, J = size(u)
    B    = (1 / Δ_HJB + ρ) * speye(T, I*J) - A
    b    = reshape(u, I*J) + reshape(V, I*J) / Δ_HJB
    V    = reshape(B \ b, I, J)
    return V, A, u, h, c, s
end

# For equilibrium_conditions; choose upwinding direction.
@inline function upwind(util::Function,
                        A_switch::SparseMatrixCSC{S, Int64},
                        cf::Matrix{T}, cb::Matrix{T},
                        c0::Matrix{T}, hf::Matrix{T},
                        hb::Matrix{T}, h0::Matrix{T},
                        sf::Matrix{T}, sb::Matrix{T},
                        Vaf::Matrix{T}, Vab::Matrix{T},
                        daf::Vector{Float64},
                        dab::Vector{Float64}) where {T<:Number, S<:Number}

    I,J = size(sb)
    h   = similar(sb)
    c   = similar(sb)
    s   = similar(sb)
    u   = similar(sb)
    X   = similar(sb)
    Z   = similar(sb)
    Y   = similar(sb)

    for i in eachindex(sb)

        Vf = (cf[i] > 0) * (util(cf[i], hf[i]) + sf[i] * Vaf[i]) + (cf[i] <= 0) * (-1e12)
        Vb = (cb[i] > 0) * (util(cb[i], hb[i]) + sb[i] * Vab[i]) + (cb[i] <= 0) * (-1e12)
        V0 = (c0[i] > 0) * util(c0[i], h0[i]) + (c0[i] <= 0) * (-1e12)

        Iunique = (sb[i] < 0) * (1 - (sf[i] > 0)) + (1 - (sb[i] < 0)) * (sf[i] > 0)
        Iboth = (sb[i] < 0) * (sf[i] > 0)

        Ib = Iunique * (sb[i] < 0) * (Vb > V0) + Iboth * (Vb == max(max(Vb, Vf), V0))
        If = Iunique * (sf[i] > 0) * (Vf > V0) + Iboth * (Vf == max(max(Vb, Vf), V0))
        I0 = 1 - Ib - If

        h[i] = hf[i] * If + hb[i] * Ib + h0[i] * I0
        c[i] = cf[i] * If + cb[i] * Ib + c0[i] * I0
        s[i] = sf[i] * If + sb[i] * Ib
        u[i] = util(c[i], h[i])

        # Construct A matrix
        X[i] = -Ib * sb[i] / dab[((i-1) % I) + 1]
        Z[i] =  If * sf[i] / daf[((i-1) % I) + 1]
        Y[i] = -Z[i] - X[i]
    end

    # R: Pretty sure the indexing for this ought be the same as in the KrusellSmith model
    X[1,:] .= T == ComplexF64 ? complex(0.) : 0.
    Z[I,:] .= T == ComplexF64 ? complex(0.) : 0.

    A = spdiagm(-1 => reshape(X,I*J)[2:I*J], 0 => reshape(Y,I*J), 1 => reshape(Z,I*J)[1:I*J-1]) + A_switch

    return A, u, h, c, s
end

# Using the market clearing condition on bonds to determine whether or not
# an equilibrium has been reached
@inline function check_bond_market_clearing(bond_err::ComplexF64, crit_S::Float64,
                                            r::Float64, r_min::Float64, r_max::Float64,
                                            r_ρ::Float64, ρ_min::Float64, ρ_max::Float64,
                                            iter_r::Bool, iter_ρ::Bool)
    clearing_condition = false
    # Using the market clearing condition on bonds to determine whether or not
    # an equilibrium has been reached
    if abs(bond_err) > crit_S
        if bond_err > 0
            if iter_r
                r_max  = r
                r      = 0.5 * (r + r_min)
            elseif iter_ρ
                ρ_min = r_ρ
                r_ρ   = 0.5 * (r_ρ + ρ_max)
            end
        else
            if iter_r
                r_min  = r
                r      = 0.5 * (r + r_max)
            elseif iter_ρ
                ρ_max = r_ρ
                r_ρ   = 0.5 * (r_ρ + ρ_min)
            end
        end
    else
        clearing_condition = true
    end
    return r, r_min, r_max, r_ρ, ρ_min, ρ_max, clearing_condition
end


text:
# Test how well neutral evolution works with different parameter settings.

function run_test_neutral_evolution( repetitions::Int64, numinputs::Int64, numinstructions_rng::IntRange, numregisters_rng::IntRange, max_steps;
      cartesian::Bool=:false, csvfile::String="" )
  if nprocs() > 2
    process_repetitions = Int(floor(repetitions/(nprocs()-1)))
    println("process_repetitions: ",process_repetitions)
    df_list = pmap( x->test_neutral_evolution( process_repetitions, numinputs, numinstructions_rng, numregisters_rng, max_steps, cartesian=cartesian ),
        collect(1:(nprocs()-1)) )
    result_df = df_list[1]
    for i = 2:length(df_list)
      result_df = average_compatible_dataframes( result_df, df_list[i] )
    end
  else
    result_df = test_neutral_evolution( repetitions, numinputs, numinstructions_rng, numregisters_rng, max_steps, cartesian=cartesian  )
  end
  if length(csvfile) > 0
    open( csvfile, "w" ) do f
      hostname = chomp(open("/etc/hostname") do f read(f,String) end)
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", default_funcs(numinputs))
      println(f,"# cartesian: ",cartesian)
      println(f,"# repetitions: ",repetitions)
      println(f,"# max_steps: ",max_steps)
      CSV.write( f, result_df, append=true, writeheader=true )
    end
  end
  result_df
end

function test_neutral_evolution( repetitions::Int64, p::Parameters, max_steps::Int64;
      cartesian::Bool=:false, print_steps::Bool=:false )
  test_neutral_evolution( repetitions, p.numinputs, p.numinteriors, p.numlevelsback, max_steps,
      cartesian=cartesian, print_steps=print_steps )
end

function test_neutral_evolution( repetitions::Int64, numinputs::Int64, numinstructions_rng::IntRange, numregisters_rng::IntRange,
      max_steps::Int64; cartesian::Bool=:false, print_steps::Bool=:false )
  df = DataFrame()
  df.numinputs = Int64[]
  if cartesian
    df.numgates = Int64[]
    df.levelsback = Int64[]
  else
    df.numinstructions = Int64[]
    df.numregisterss = Int64[]
  end
  df.steps = Float64[]
  df.fail_fract = Float64[]
  funcs = lin_funcs(numinputs)
  for ni in numinstructions_rng
    for nr in numregisters_rng
      println("nr: ",nr)
      p = Parameters(numinputs,1,ni,nr)
      sum_steps = 0
      sum_fail = 0
      for i = 1:repetitions
        gl = randgoal(p)
        if cartesian
          c = random_chromosome( p )
        else
          c = rand_lcircuit(p,funcs)
        end
        (c,step) = neutral_evolution( c, gl, max_steps, print_steps=print_steps )
        if repetitions > 1000 && i % 100 == 0
          println("i: ",i,"  gl: ",@sprintf("0x%x",gl[1]),"  ni: ",ni,"  nr: ",nr,"  step: ",step)
        elseif repetitions <= 1000
          println("i: ",i,"  gl: ",@sprintf("0x%x",gl[1]),"  ni: ",ni,"  nr: ",nr,"  step: ",step)
        end
        sum_steps += step
        sum_fail += step==max_steps ? 1 : 0
      end
      df_row = ( numinputs, ni, nr, sum_steps/repetitions, sum_fail/repetitions )
      push!(df,df_row)
    end
  end
  df
end

function average_compatible_dataframes( df1::DataFrame, df2::DataFrame )
  @assert names(df1) == names(df2)
  @assert size(df1) == size(df2)
  result_df = DataFrame()
  for i = 1:size(df2)[2]
    if df1[:,i] == df2[:,i]
      insertcols!(result_df,i,names(df1)[i]=>df1[:,i])
    else
      insertcols!(result_df,i,names(df1)[i]=>(df1[:,i]+df2[:,i])./2.0)
    end
  end
  result_df
end


text:
"
Generator represents the way that a population is going to be generated.

# Fields
-`_popSize::UInt32`: size of the population, number of individuals.
-`_generateOneByOne::Bool`: specifies wether the generation method generates
    a single individual each run or the whole population at a time.
-`_method::Function`: method used for generating the individuals.
-`_varArgs::Array{Any}`: arguments necessary for the generation method.
"
mutable struct Generator
    _popSize::UInt32
    _generateOneByOne::Bool
    _method::Function
    _varArgs::Array{Any}
end # struct

"""
    getPopSize(gen::Generator)::UInt32

Returns the population size assigned.
"""
getPopSize(gen::Generator)::UInt32 = gen._popSize
# function

"""
    getMethod(gen::Generator)::Function

Returns the method for generate the population.
"""
getMethod(gen::Generator)::Function = gen._method
# function

"""
    isGeneratedOneByOne(gen::Generator)::Bool

Returns whether the method generates individual one by one or all individuals in a single call
"""
isGeneratedOneByOne(gen::Generator)::Bool = gen._generateOneByOne
#function

"""
    getFunctionArgs(gen::Generator)::Array

Obtains the aditional arguments associated to initialization method.
"""
getFunctionArgs(gen::Generator)::Array = gen._varArgs
#function

"""
    genPopulation_(generator::Generator, experimentInfo::ExperimentInfo)::Array{Individual}

Generates a population of individuals according to the method given by the Generator.
"""
function genPopulation_(generator::Generator, experimentInfo::ExperimentInfo)::Array{Individual}

    popSize = getPopSize(generator)
    population = Array{Individual}(undef, popSize)
    rng = getRNG(experimentInfo)

    method = getMethod(generator)
    functionArgs = getFunctionArgs(generator)

    if getIndividualType(experimentInfo) <: GAGenotype
        if isGeneratedOneByOne(generator)
            for i = 1:popSize
                genotype = method(rng, functionArgs...)
                population[i] = Individual(genotype)
            end
        else
            populationGenotype = method(popSize, rng, functionArgs...)
            map!(x->Individual(x), population, populationGenotype)
        end
    else
        if isGeneratedOneByOne(generator)
            for i = 1:popSize
                genotype = method(experimentInfo._GPExperimentInfo, rng, functionArgs...)
                population[i] = Individual(genotype)
            end
        else
            populationGenotype = method(experimentInfo._GPExperimentInfo, popSize, rng, functionArgs...)
            map!(x->Individual(x), population, populationGenotype)
        end
    end

    return population
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
module GeometryPrimitives
using StaticArrays, LinearAlgebra
using Statistics: mean
using Makie

export Shape, Shape1, Shape2, Shape3
export level, surfpt_nearby, normal, bounds, translate

abstract type Shape{N,N²} end # a solid geometric shape in N dimensions (N² = N*N is needed in some shapes, e.g., Cuboid)
const Shape1 = Shape{1,1}
const Shape2 = Shape{2,4}
const Shape3 = Shape{3,9}

Base.ndims(o::Shape{N}) where {N} = N

# The following functions return Any due to the limitations of Julia's dispatch system.
# Therefore, always call them with return type assertions.  See
# https://discourse.julialang.org/t/extending-base-in-type-stably/5341/12
# https://github.com/JuliaLang/julia/issues/23210
level(x::AbstractVector{<:Real}, o::Shape{N}) where {N} = level(SVector{N}(x), o)
Base.in(x::AbstractVector{<:Real}, o::Shape{N}) where {N} = level(x,o) ≤ 0
surfpt_nearby(x::AbstractVector{<:Real}, o::Shape{N}) where {N} = surfpt_nearby(SVector{N}(x), o)
normal(x::AbstractVector{<:Real}, o::Shape) = surfpt_nearby(x, o)[2]  # outward direction even for x inside o
translate(s::Shape{N}, ∆::AbstractVector{<:Real}) where {N} = translate(s, SVector{N}(∆))
translate(s::Shape{N}, ∆::SVector{N,<:Real}) where {N} = (s2 = deepcopy(s); s2.c += ∆; s2)  # default implementation

function orthoaxes(n::SVector{3,<:Real})
    u_temp = abs(n[3]) < abs(n[1]) ? SVector(0,0,1) : SVector(1,0,0)
    v = normalize(n × u_temp)
    u = v × n

    return u, v
end

function orthoaxes(n::SVector{N,<:Real}) where {N}
    u_temp = abs(n[3]) < abs(n[1]) ? SVector(0,0,1) : SVector(1,0,0)
    v = normalize(n × u_temp)
    u = v × n

    return u, v
end

include("hyper/hyper.jl")
include("planar/planar.jl")
include("prism/prism.jl")
include("util/util.jl")

end # module


text:
"""
Return path to the notebook starting with "/docs/src/".
The "/docs/src/" assumption is required since we don't know the module.
"""
function _relative_notebook_path(bopts::BuildOptions, in_path::String)::Union{String,Nothing}
    absolute_path = joinpath(bopts.dir, in_path)
    pos = findfirst("/docs/src/", absolute_path)
    if !isnothing(pos)
        relative_path = absolute_path[first(pos) + 1:end]
        return string(relative_path)
    else
        return nothing
    end
end

function _editurl_text(
        bopts::BuildOptions,
        github_repo::String,
        branch::String,
        in_path::String
    )
    path = _relative_notebook_path(bopts, in_path)
    url = "https://github.com/$github_repo/blob/$branch/$path"
    return """
        ```@meta
        EditURL = "$url"
        ```
        """
end

function _editurl_text(bopts::BuildOptions, in_path::String)::String
    github_repo = get(ENV, "GITHUB_REPOSITORY", "")
    if github_repo != ""
        # Only when we're in GitHub's CI, it is possible to figure out the path to the file.
        # "refs/heads/$(branchname)" for branch, "refs/tags/$(tagname)" for tags.
        # Thanks to Documenter for this comment.
        github_ref = get(ENV, "GITHUB_REF", "")
        if github_ref != ""
            prefix = "refs/heads/"
            if contains(github_ref, prefix)
                branch = github_ref[length(prefix) + 1:end]
                return _editurl_text(bopts, github_repo, branch, in_path)
            else
                # We're on a tag so maybe that's a triggered deploy for stable docs?
                branch = "main"
                return _editurl_text(bopts, github_repo, branch, in_path)
            end
        else
            @warn "github_ref was empty which is unexpected"
            return ""
        end
    else
        return ""
    end
end


text:
@checked struct LLVMError
    ref::API.LLVMErrorRef
end
Base.unsafe_convert(::Type{API.LLVMErrorRef}, err::LLVMError) = err.ref

"""
    consume(err::LLVMError)

Consumes the error without handling it.
"""
function consume(err::LLVMError)
    API.LLVMConsumeError(err)
end

# NOTE: handles the error, calling consume on `err` is illegal
function Base.convert(::Type{LLVMException}, err::LLVMError)
    errmsg = API.LLVMGetErrorMessage(err)
    local msg
    try
        msg = Base.unsafe_string(errmsg)
    finally
        API.LLVMDisposeErrorMessage(errmsg)
    end
    LLVMException(msg)
end

macro check(expr)
    quote
        err_ref = $(esc(expr))
        if err_ref != C_NULL
            err = LLVMError(err_ref)
            throw(convert(LLVMException, err))
        end
    end
end


text:
include("switchedsystems.jl")
include("sosdemo2.jl")


text:
using Catlab
using LibPQ, Tables

include(joinpath(@__DIR__, "FLS.jl"))

conn = LibPQ.Connection("dbname=models")
result = execute(conn, "SELECT typname FROM pg_type WHERE oid = 16")
data = columntable(result)

"""Data to specify a table"""
struct Tab
  name::String
  attrs::Vector{Tuple{String, String, Bool}}
  fks::Vector{Union{Pair{String,String},String}}
end

"""SQL create table statement"""
function add_tab(t::Tab, conn::LibPQ.Connection)::Nothing
  attrs = ["$k $v NOT NULL " * (u ? "UNIQUE" : "") for (k, v, u) in t.attrs]
  fks = [let (key, tgt) = k isa String ? ("$(k)_id", k) : k;
          ("$key INTEGER NOT NULL",
           "FOREIGN KEY ($key) REFERENCES $tgt ($(tgt)_id)") end
       for k in t.fks]
  fk1, fk2 = isempty(fks) ? ([],[]) : map(collect, zip(fks...)) # all attr decs, then all constraints
  parentunique = t.name == "PARENT" ? ["UNIQUE (parent, child)"] : []
  stmt = """CREATE TABLE IF NOT EXISTS $(t.name)
    \t($(t.name)_id SERIAL PRIMARY KEY,
    \t$(join(vcat(attrs,fk1, fk2, parentunique), ",\n\t")))"""
  execute(conn, stmt)
  return nothing
end

tabs = [Tab("FLS", [("hash", "DECIMAL(20,0)", true),
                    ("name", "TEXT", true),
                    ("jdump","TEXT", true)], []),
    Tab("Premodel", [("hash", "DECIMAL(20,0)", true),
                     ("jdump", "TEXT", true),
                     ("size", "INTEGER", false),
                     ("fired", "BOOLEAN", false)], ["FLS"]), # ought be boolean
    Tab("Model", [("hash", "DECIMAL(20,0)", true),
                  ("jdump", "TEXT", true)], ["Premodel"]),
    Tab("Init", [("cards", "TEXT", true)], ["Premodel"]),
    Tab("Parent", [("delta", "TEXT", false)], ["parent"=>"Premodel", "child" => "Premodel"])]

function init_db(; reset::Bool=false)::LibPQ.Connection
  conn = LibPQ.Connection("dbname=models")

  if reset
    for t in reverse(tabs)
      execute(conn, "DROP TABLE IF EXISTS $(t.name)")
    end
    for t in tabs
      add_tab(t, conn)
    end
  end
  return conn
end

"""
Given initial cardinalities for every nonlimit object, log the premodel they
correspond to.
"""
function init_premodel(db::LibPQ.Connection, F::FLS, cards::Vector{Int})::Int
  syms = sort(collect(realobs(F)))
  length(syms) == length(cards) || error("bad length: cards $cards syms $syms")
  d = Dict(zip(syms, cards))
  dstr = string(sort(collect(d)))
  r = columntable(execute(db, "SELECT Premodel_id FROM Init WHERE cards=\$1",
                              [dstr]))
  if !isempty(r[:premodel_id])
    return r[:premodel_id][1]
  end

  i = add_premodel(db, F, initrel(F, d))
  execute(db, "INSERT INTO Init (cards, Premodel_id) VALUES (\$1,\$2)", [dstr,i])
  return i
end

"""Add FLS (if not already in DB) and return id"""
function add_fls(db::LibPQ.Connection, F::FLS)::Int
  ghash = canonical_hash(to_graph(F.schema))
  # check if already in db
  r = columntable(execute(db, "SELECT FLS_id FROM FLS WHERE hash=\$1",[ghash]))
  if !(isempty(r[:fls_id]))
    return r[:fls_id][1]
  end

  z = columntable(execute(db,
    "INSERT INTO FLS (hash, name, jdump) VALUES (\$1,\$2,\$3) RETURNING FLS_id",
    [ghash,F.name,to_json(F)]))[:fls_id][1]
  return z
end

"""Get FLS by primary key id"""
function get_fls(db::LibPQ.Connection, i::Int)::FLS
  z = columntable(execute(db, "SELECT jdump FROM FLS WHERE FLS_id=\$1",[i]))
  return fls_from_json(z[:jdump][1])
end

"""Add premodel (if it doesn't exist) """
function add_premodel(db::LibPQ.Connection, F::FLS, m::StructACSet,
                      parent::Union{Int, Nothing}=nothing,
                      csd::Union{ChaseStepData, Nothing}=nothing,
                      chash::Union{UInt64, Nothing}=nothing)::Int
  chash = chash === nothing ? canonical_hash(m) : chash
  # check if already in db
  r = columntable(execute(db,
    "SELECT Premodel_id FROM Premodel WHERE hash=\$1",[chash]))[:premodel_id]
  if isempty(r)
    fid = add_fls(db, F)
    r = columntable(execute(db,
      """INSERT INTO Premodel (hash, jdump, FLS_id, fired, size)
         VALUES (\$1,\$2, \$3, false, \$4) RETURNING Premodel_id""",
            [chash,generate_json_acset(m), fid, relsize(F, m)]))[:premodel_id]
  end

  if parent isa Int # Add parents if any
    x =  columntable(execute(db,
      "SELECT Parent_id FROM Parent WHERE parent=\$1 AND child=\$2",
      [parent, r[1]]))[:parent_id]
    if isempty(x)
      execute(db, "INSERT INTO Parent(parent, child, delta) VALUES(\$1, \$2, \$3)", [parent, r[1], JSON.json(csd_to_dict(csd))])
    end
  end
  return r[1]
end

"""Add a premodel that is a model (if not already there). Return id"""
function add_model(db::LibPQ.Connection, F::FLS, relm::StructACSet,
                   parent::Union{Int, Nothing}=nothing,
                   csd::Union{ChaseStepData, Nothing}=nothing,
                   relm_hsh::Union{UInt64, Nothing}=nothing)::Int
  m, partial = crel_to_cset(F, relm) # fail if nonfunctional
  !partial || error("adding a partial model")
  chash = canonical_hash(m)
  pid = add_premodel(db, F, relm, parent, csd, relm_hsh)
  # check if already in db
  r = columntable(execute(db, "SELECT Model_id FROM Model WHERE hash=\$1",
                              [chash]))[:model_id]
  if !(isempty(r))
    return r[1]
  end

  z = columntable(execute(db, """INSERT INTO Model (hash, jdump, Premodel_id)
                                 VALUES (\$1,\$2,\$3) RETURNING Model_id""",
                          [chash, generate_json_acset(m), pid]))[:model_id]
  return z[1]
end

function get_premodel_id(db::LibPQ.Connection, crel::StructACSet)::Union{Nothing, Int}
  z = columntable(execute(db, """SELECT Premodel_id FROM Premodel
      WHERE hash=\$1""",[canonical_hash(crel)]))[:premodel_id]
  return (isempty(z) ? nothing : z[1])
end

"""Get premodel by primary key ID"""
function get_premodel(db::LibPQ.Connection, i::Int)::Pair{FLS, StructACSet}
  z = columntable(execute(db, """SELECT FLS.jdump AS f, Premodel.jdump as p
    FROM Premodel JOIN FLS USING (FLS_id) WHERE Premodel_id=\$1""",[i]))
  fls = fls_from_json(z[:f][1])
  crel = parse_json_acset(fls.crel, z[:p][1])
  return fls => crel
end

"""Get Model by primary key ID"""
function get_model(db::LibPQ.Connection, i::Int)::Pair{FLS, StructACSet}
  z = columntable(execute(db, """SELECT FLS.jdump AS f, Model.jdump as p
    FROM Model JOIN Premodel USING (Premodel_id)
               JOIN FLS USING (FLS_id) WHERE Model_id=\$1""",[i]))
  fls = fls_from_json(z[:f][1])
  cset = parse_json_acset(fls.cset, z[:p][1])
  return fls => cset
end

"""All premodel hashes for an FLS that have been chased"""
function get_seen(db::LibPQ.Connection, fls::FLS)::Set{UInt64}
  fid = add_fls(db, fls)
  z = columntable(execute(db,
    "SELECT hash FROM Premodel WHERE fired AND FLS_id=\$1",[fid]))
  return Set(map(UInt64, z[:hash]))
end