text:
# Functional implementation of https://egraphs-good.github.io/
# https://dl.acm.org/doi/10.1145/3434304

using DataStructures

"""
Abstract type representing an [`EGraph`](@ref) analysis,
attaching values from a join semi-lattice domain to
an EGraph
"""
const ClassMem = Dict{Int64,EClass}
const HashCons = Dict{ENode,Int64}
const Analyses = Set{Type{<:AbstractAnalysis}}
const SymbolCache = Dict{Any, Vector{Int64}}

"""
A concrete type representing an [`EGraph`].
See the [egg paper](https://dl.acm.org/doi/pdf/10.1145/3434304)
for implementation details
"""
mutable struct EGraph
    """stores the equality relations over e-class ids"""
    U::IntDisjointSets
    """map from eclass id to eclasses"""
    M::ClassMem             #
    H::HashCons             # hashcons
    """worklist for ammortized upwards merging"""
    dirty::Vector{Int64}
    root::Int64
    """A vector of analyses associated to the EGraph"""
    analyses::Analyses
    """
    a cache mapping function symbols to e-classes that
    contain e-nodes with that function symbol.
    """
    symcache::SymbolCache
end

EGraph() = EGraph(
    IntDisjointSets(0),
    ClassMem(),
    HashCons(),
    # ParentMem(),
    Vector{Int64}(),
    0,
    Analyses(),
    SymbolCache()
)

function EGraph(e)
    g = EGraph()
    rootclass = addexpr!(g, e)
    g.root = rootclass.id
    g
end

function canonicalize(g::EGraph, n::ENode)
    new_args = map(x -> find(g, x), n.args)
    typeof(n)(n.head, new_args, n.metadata)
end

function canonicalize!(g::EGraph, n::ENode)
    for i ∈ 1:ariety(n)
        n.args[i] = find(g, n.args[i])
    end
    return n
    # n.args = map(x -> find(g, x), n.args)
end

"""
Returns the canonical e-class id for a given e-class.
"""
find(g::EGraph, a::Int64)::Int64 = find_root!(g.U, a)
find(g::EGraph, a::EClass)::Int64 = find_root!(g.U, a.id)

geteclass(g::EGraph, a::Int64)::EClass = g.M[find(g, a)]
geteclass(g::EGraph, a::EClass)::Int64 = a

### Definition 2.3: canonicalization
# iscanonical(U::IntDisjointSets, n::Expr) = n == canonicalize(U, n)
iscanonical(g::EGraph, n::ENode) = n == canonicalize(g, n)
iscanonical(g::EGraph, e::EClass) = find(g, e.id) == e.id

function canonicalize!(g::EGraph, e::EClass)
    e.id = find(g, e.id)
end

"""
Inserts an e-node in an [`EGraph`](@ref)
"""
function add!(g::EGraph, n::ENode)::EClass
    @debug("adding ", n)

    n = canonicalize!(g, n)
    if haskey(g.H, n)
        return geteclass(g, find(g, g.H[n]))
    end
    @debug(n, " not found in H")

    id = push!(g.U) # create new singleton eclass

    for c_id ∈ n.args
        addparent!(g.M[c_id], n, id)
    end

    g.H[n] = id

    classdata = EClass(id, OrderedSet([n]), OrderedDict{ENode, Int64}())
    g.M[id] = classdata

    # cache the eclass for the symbol for faster matching
    sym = n.head
    if !haskey(g.symcache, sym)
        g.symcache[sym] = Int64[]
    end
    push!(g.symcache[sym], id)

    # make analyses for new enode
    for an ∈ g.analyses
        if !islazy(an)
            setdata!(classdata, an, make(an, g, n))
            modify!(an, g, id)
        end
    end

    return classdata
end

"""
Recursively traverse an [`Expr`](@ref) and insert terms into an
[`EGraph`](@ref). If `e` is not an [`Expr`](@ref), then directly
insert the literal into the [`EGraph`](@ref).
"""
function addexpr_rec!(g::EGraph, e)::EClass
    # e = preprocess(e)
    # println("========== $e ===========")
    if e isa EClass
        return e
    end

    if istree(e)
        args = getargs(e)
        n = length(args)
        class_ids = Vector{Int64}(undef, n)
        for i ∈ 1:n
            # println("child $child")
            @inbounds child = args[i]
            c_eclass = addexpr!(g, child)
            @inbounds class_ids[i] = c_eclass.id
        end
        node = ENode(e, class_ids)
        return add!(g, node)
    end

    return add!(g, ENode(e))
end

addexpr!(g::EGraph, e) = addexpr_rec!(g, preprocess(e))

function clean_enode!(g::EGraph, t::ENode, to::Int64)
    nt = canonicalize!(g, t)
    if nt != t
        delete!(g.H, t)
    end
    g.H[nt] = to
    return t
end

"""
Given an [`EGraph`](@ref) and two e-class ids, set
the two e-classes as equal.
"""
function Base.merge!(g::EGraph, a::Int64, b::Int64)::Int64
    id_a = find(g, a)
    id_b = find(g, b)
    id_a == id_b && return id_a
    to = union!(g.U, id_a, id_b)

    @debug "merging" id_a id_b

    from = (to == id_a) ? id_b : id_a

    push!(g.dirty, to)

    from_class = g.M[from]
    to_class = g.M[to]

    g.M[to] = union!(from_class, to_class)
    # g.M[from] = g.M[to]
    delete!(g.M, from)

    # mutable version
    # for an ∈ g.analyses
    #     if hasdata(from_class, an) && hasdata(to_class, an)
    #         from_data = getdata(from_class, an)
    #         to_data = getdata(to_class, an)
    #         setdata!(to_class, an, join(an, from_data, to_data))
    #     end
    # end

    return to
end

"""
This function restores invariants and executes
upwards merging in an [`EGraph`](@ref). See
the [egg paper](https://dl.acm.org/doi/pdf/10.1145/3434304)
for more details.
"""
function rebuild!(egraph::EGraph)
    while !isempty(egraph.dirty)
        todo = unique([find(egraph, id) for id ∈ egraph.dirty])
        empty!(egraph.dirty)
        for x ∈ todo
            repair!(egraph, x)
        end
    end

    for (sym, ids) ∈ egraph.symcache
        egraph.symcache[sym] = unique(ids .|> x -> find(egraph, x))
    end

    if egraph.root != 0
        egraph.root = find(egraph, egraph.root)
    end

    # for i ∈ 1:length(egraph.U)
    #     find_root!(egraph.U, i)
    # end
    # INVARIANTS ASSERTIONS
    # for (id, c) ∈  egraph.M
    # #     ecdata.nodes = map(n -> canonicalize(egraph.U, n), ecdata.nodes)
    #     for an ∈ egraph.analyses
    #         if haskey(an, id)
    #             @assert an[id] == mapreduce(x -> make(an, x), (x, y) -> join(an, x, y), c.nodes)
    #         end
    #     end
    #
    #     for n ∈ c
    #         # println(n)
    #         # println("canon = ", canonicalize(egraph, n))
    #         @assert egraph.H[canonicalize(egraph, n)] == find(egraph, id)
    #     end
    # end
end

function repair!(g::EGraph, id::Int64)
    id = find(g, id)
    ecdata = g.M[id]
    @debug "repairing " id

    for (p_enode, p_eclass) ∈ ecdata.parents
        clean_enode!(g, p_enode, find(g, p_eclass))
    end

    new_parents = OrderedDict{ENode,Int64}()

    for (p_enode, p_eclass) ∈ ecdata.parents
        p_enode = canonicalize!(g, p_enode)
        # deduplicate parents
        if haskey(new_parents, p_enode)
            @debug "merging classes" p_eclass (new_parents[p_enode])
            merge!(g, p_eclass, new_parents[p_enode])
        end
        new_parents[p_enode] = find(g, p_eclass)
    end
    ecdata.parents = new_parents
    @debug "updated parents " id g.parents[id]

    # ecdata.nodes = map(n -> canonicalize(g.U, n), ecdata.nodes)

    # Analysis invariant maintenance
    for an ∈ g.analyses
        hasdata(ecdata, an) && modify!(an, g, id)
        # modify!(an, id)
        # id = find(g, id)
        for (p_enode, p_id) ∈ ecdata.parents
            # p_eclass = find(g, p_eclass)
            p_eclass = g.M[p_id]
            if !islazy(an) && !hasdata(p_eclass, an)
                setdata!(p_eclass, an, make(an, g, p_enode))
            end
            if hasdata(p_eclass, an)
                p_data = getdata(p_eclass, an)

                new_data = join(an, p_data, make(an, g, p_enode))
                if new_data != p_data
                    setdata!(p_eclass, an, new_data)
                    push!(g.dirty, p_id)
                end
            end
        end
    end

    # ecdata.nodes = map(n -> canonicalize(g.U, n), ecdata.nodes)

end

"""
Recursive function that traverses an [`EGraph`](@ref) and
returns a vector of all reachable e-classes from a given e-class id.
"""
function reachable(g::EGraph, id::Int64)
    id = find(g, id)
    hist = Int64[id]
    todo = Int64[id]
    while !isempty(todo)
        curr = find(g, pop!(todo))
        for n ∈ g.M[curr]
            nn = canonicalize(g, n)
            # println("node in reachability is ", n)
            for c_id ∈ nn.args
                if c_id ∉ hist
                    push!(hist, c_id)
                    push!(todo, c_id)
                end
            end
        end
    end

    return hist
end


text:
using BranchAndPrune
using Test

import BranchAndPrune: BPLeaf, BPNode, BPTree, MissingImplementationError
import BranchAndPrune: discard_leaf!, newid, root

include("tree.jl")
include("interface.jl")

ext: jl


text:
using IProfile

@iprofile begin
function f1(x::Int)
    k = 5
    l = 12
    y = sum(1:1000)
    for j = 1:100000
        z = j^2
    end
    if z > 0
        z = z-1
    elseif z == 0
        z = 1000
    else
        z = z+1
    end
    return x+1
end
function f1(x::Float64)
    return x+2
end

function f1{T}(x::T)
    return x+5
end

unused_top_level_var = 2

# Test ternary operators with trivial return values. See issue #17.
function f2()
    rand()<0.5 ? 1 : 0
end

end #@iprofile begin

f1(7)
f1(7.3)
println("Less than 0.5: ", f2())

for i = 1:200
    f1(Int16(i))
end

@iprofile report
@iprofile clear


text:
module SPN

include("nodes.jl")
export Node, InnerNode, LeafNode, SumNode, ProdNode, IndicatorNode, GaussianNode
export connect!, setInput!, eval!, passDerivative!

include("network.jl")
export SumProductNetwork, computeDerivatives!
export numNodes, numSumNodes, numProdNodes, numLeafNodes

include("inference.jl")
export marginalInference!, conditionalInference!, mpeInference!

include("utils.jl")
export addExp, quantileMeans, dataLikelihood!

include("learn/parameter.jl")
export parameterLearnHardEM!, parameterLearnEM!, parameterLearnGD!

include("learn/poon.jl")
export structureLearnPoon

end # module


text:
using Plots, LaTeXStrings, Statistics

function Deposition(;len, tot_time)
    surf = [0 for i=1:len]
    AvgList = [0.0 for i=1:tot_time]
    for n in 2:tot_time+1
        randsurf = rand(1:len)
        for i in randsurf
            MAX = FindMax(surf, i, len)
            surf[i] = MAX
        end
        AvgList[n-1] = mean(surf)
    end
    return AvgList
end

function sides(n, L)
    if n == L
        return n-1 , 1
    elseif n == 1
        return L , n+1
    else
        return n-1, n+1
    end
end

function FindMax(surface, index_,L_surf)
    i1 , i2 = sides(index_, L_surf)
    maxlen = max(surface[i1],surface[index_] + 1,surface[i2])
    return maxlen
end
function Linear_fit(;Time, AvgList, tot_time)
    A = [hcat(Time[1:tot_time]) reshape(ones(tot_time), tot_time, 1)]
    b = reshape(AvgList[1:tot_time], tot_time, 1)
    line = (A \ b)
    return line
end

iternum = 1000
Parameters = Dict(:len => 300,
                    :tot_time => 50,)
allAvg = [ [0.0 for i in 1:Parameters[:tot_time]] for j = 1:iternum]
meanAvg = [0.0 for i in 1:Parameters[:tot_time]]
vars = [0.0 for i in 1:Parameters[:tot_time]]
for i in 1:iternum
    AvgList = Deposition(;Parameters...)
    allAvg[i] = AvgList
    meanAvg += AvgList
    print("\r$i")
end
meanAvg /= iternum
for i in 1:Parameters[:tot_time]
    vars[i] = std(hcat(allAvg...)[i,:])
end

Time = 0:Parameters[:tot_time]-1

Paraline = Dict(
                :Time => Time[1:end],
                :AvgList => meanAvg[1:end],
                :tot_time => 50)
Line = Linear_fit(;Paraline...)
X = Time[1]:Time[50]
Y = X .* Line[1] .+ Line[2]
Line[2]
plot(X,Y,c=:black,label = L"y = %$(round(Line[1],digits= 3))x + %$(round(Line[2],digits= 3))")
scatter!(Time[1:end], meanAvg[1:end],
    # xlims = (1, Parameters[:tot_time]),
    c = :steelblue,
    xlabel= L"Time",
    ylabel= L"H_{(t)}",
    title= L"~H_{(t)}-Time~\ (L = %$(Parameters[:len]))",
    label = L"Data\ point",
    yerror = vars)
    # legend = nothing)

savefig("../../Figs/Q2/H-t(L=$(Parameters[:len])).pdf")


text:
function update_aux!(edmf::EDMF_PrognosticTKE{N_up}, gm, grid, state, Case, param_set, TS) where {N_up}
    #####
    ##### Unpack common variables
    #####
    kc_surf = kc_surface(grid)
    kf_surf = kf_surface(grid)
    kc_toa = kc_top_of_atmos(grid)
    up = edmf.UpdVar
    en = edmf.EnvVar
    en_thermo = edmf.EnvThermo
    ρ0_f = face_ref_state(state).ρ0
    p0_c = center_ref_state(state).p0
    ρ0_c = center_ref_state(state).ρ0
    α0_c = center_ref_state(state).α0
    g = CPP.grav(param_set)
    c_m = CPEDMF.c_m(param_set)
    KM = center_aux_turbconv(state).KM
    KH = center_aux_turbconv(state).KH
    surface = Case.Sur
    obukhov_length = surface.obukhov_length
    FT = eltype(grid)
    prog_gm = center_prog_grid_mean(state)
    prog_gm_f = face_prog_grid_mean(state)
    aux_up = center_aux_updrafts(state)
    aux_up_f = face_aux_updrafts(state)
    aux_en = center_aux_environment(state)
    aux_en_f = face_aux_environment(state)
    aux_gm = center_aux_grid_mean(state)
    aux_tc_f = face_aux_turbconv(state)
    aux_tc = center_aux_turbconv(state)
    aux_bulk = center_aux_bulk(state)
    prog_en = center_prog_environment(state)
    aux_en_2m = center_aux_environment_2m(state)
    prog_up = center_prog_updrafts(state)
    prog_up_f = face_prog_updrafts(state)
    aux_en_unsat = aux_en.unsat
    aux_en_sat = aux_en.sat
    m_entr_detr = aux_tc.ϕ_temporary
    ∇m_entr_detr = aux_tc.ψ_temporary
    w_up_c = aux_tc.w_up_c
    w_en_c = aux_tc.w_en_c
    wvec = CC.Geometry.WVector
    Δt = TS.dt
    max_area = edmf.max_area
    sde_model = edmf.sde_model

    #####
    ##### center variables
    #####

    @inbounds for k in real_center_indices(grid)
        #####
        ##### Set primitive variables
        #####
        @inbounds for i in 1:N_up
            if is_surface_center(grid, k)
                if prog_up[i].ρarea[k] / ρ0_c[k] >= edmf.minimum_area
                    aux_up[i].θ_liq_ice[k] = edmf.h_surface_bc[i]
                    aux_up[i].q_tot[k] = edmf.qt_surface_bc[i]
                    aux_up[i].area[k] = edmf.area_surface_bc[i]
                else
                    aux_up[i].θ_liq_ice[k] = prog_gm.θ_liq_ice[k]
                    aux_up[i].q_tot[k] = prog_gm.q_tot[k]
                end
            else
                if prog_up[i].ρarea[k] / ρ0_c[k] >= edmf.minimum_area
                    aux_up[i].θ_liq_ice[k] = prog_up[i].ρaθ_liq_ice[k] / prog_up[i].ρarea[k]
                    aux_up[i].q_tot[k] = prog_up[i].ρaq_tot[k] / prog_up[i].ρarea[k]
                    aux_up[i].area[k] = prog_up[i].ρarea[k] / ρ0_c[k]
                else
                    aux_up[i].θ_liq_ice[k] = prog_gm.θ_liq_ice[k]
                    aux_up[i].q_tot[k] = prog_gm.q_tot[k]
                    aux_up[i].area[k] = 0
                end
            end
        end

        #####
        ##### compute bulk
        #####
        aux_bulk.q_tot[k] = 0
        aux_bulk.θ_liq_ice[k] = 0
        aux_bulk.area[k] = sum(i -> aux_up[i].area[k], 1:N_up)
        if aux_bulk.area[k] > 0
            @inbounds for i in 1:N_up
                aux_bulk.q_tot[k] += aux_up[i].area[k] * aux_up[i].q_tot[k] / aux_bulk.area[k]
                aux_bulk.θ_liq_ice[k] += aux_up[i].area[k] * aux_up[i].θ_liq_ice[k] / aux_bulk.area[k]
            end
        else
            aux_bulk.q_tot[k] = prog_gm.q_tot[k]
            aux_bulk.θ_liq_ice[k] = prog_gm.θ_liq_ice[k]
        end
        aux_en.area[k] = 1 - aux_bulk.area[k]
        aux_en.tke[k] = prog_en.ρatke[k] / (ρ0_c[k] * aux_en.area[k])
        aux_en.Hvar[k] = prog_en.ρaHvar[k] / (ρ0_c[k] * aux_en.area[k])
        aux_en.QTvar[k] = prog_en.ρaQTvar[k] / (ρ0_c[k] * aux_en.area[k])
        aux_en.HQTcov[k] = prog_en.ρaHQTcov[k] / (ρ0_c[k] * aux_en.area[k])

        #####
        ##### decompose_environment
        #####
        a_bulk_c = aux_bulk.area[k]
        val1 = 1 / (1 - a_bulk_c)
        val2 = a_bulk_c * val1
        aux_en.q_tot[k] = max(val1 * prog_gm.q_tot[k] - val2 * aux_bulk.q_tot[k], 0) #Yair - this is here to prevent negative QT
        aux_en.θ_liq_ice[k] = val1 * prog_gm.θ_liq_ice[k] - val2 * aux_bulk.θ_liq_ice[k]

        #####
        ##### saturation_adjustment and buoyancy
        #####
        ts_gm = thermo_state_pθq(param_set, p0_c[k], prog_gm.θ_liq_ice[k], prog_gm.q_tot[k])
        aux_tc.θ_virt[k] = TD.virtual_pottemp(ts_gm)

        ts_en = thermo_state_pθq(param_set, p0_c[k], aux_en.θ_liq_ice[k], aux_en.q_tot[k])

        aux_en.T[k] = TD.air_temperature(ts_en)
        aux_en.θ_virt[k] = TD.virtual_pottemp(ts_en)
        aux_en.θ_dry[k] = TD.dry_pottemp(ts_en)
        aux_en.q_liq[k] = TD.liquid_specific_humidity(ts_en)
        aux_en.q_ice[k] = TD.ice_specific_humidity(ts_en)
        rho = TD.air_density(ts_en)
        aux_en.buoy[k] = buoyancy_c(param_set, ρ0_c[k], rho)

        update_sat_unsat(en_thermo, state, k, ts_en)
        aux_en.RH[k] = TD.relative_humidity(ts_en)

        @inbounds for i in 1:N_up
            if aux_up[i].area[k] < edmf.minimum_area && k > kc_surf && aux_up[i].area[k - 1] > 0.0
                qt = aux_up[i].q_tot[k - 1]
                h = aux_up[i].θ_liq_ice[k - 1]
                ts_up = thermo_state_pθq(param_set, p0_c[k], h, qt)
            else
                ts_up = thermo_state_pθq(param_set, p0_c[k], aux_up[i].θ_liq_ice[k], aux_up[i].q_tot[k])
            end
            aux_up[i].q_liq[k] = TD.liquid_specific_humidity(ts_up)
            aux_up[i].q_ice[k] = TD.ice_specific_humidity(ts_up)
            aux_up[i].T[k] = TD.air_temperature(ts_up)
            ρ = TD.air_density(ts_up)
            aux_up[i].buoy[k] = buoyancy_c(param_set, ρ0_c[k], ρ)
            aux_up[i].RH[k] = TD.relative_humidity(ts_up)
        end
        aux_gm.buoy[k] = (1.0 - aux_bulk.area[k]) * aux_en.buoy[k]
        @inbounds for i in 1:N_up
            aux_gm.buoy[k] += aux_up[i].area[k] * aux_up[i].buoy[k]
        end
        @inbounds for i in 1:N_up
            aux_up[i].buoy[k] -= aux_gm.buoy[k]
        end
        aux_en.buoy[k] -= aux_gm.buoy[k]

        #####
        ##### compute bulk thermodynamics
        #####
        aux_bulk.q_liq[k] = 0
        aux_bulk.q_ice[k] = 0
        aux_bulk.T[k] = 0
        aux_bulk.RH[k] = 0
        aux_bulk.buoy[k] = 0
        if a_bulk_c > 0
            @inbounds for i in 1:N_up
                aux_bulk.q_liq[k] += aux_up[i].area[k] * aux_up[i].q_liq[k] / a_bulk_c
                aux_bulk.q_ice[k] += aux_up[i].area[k] * aux_up[i].q_ice[k] / a_bulk_c
                aux_bulk.T[k] += aux_up[i].area[k] * aux_up[i].T[k] / a_bulk_c
                aux_bulk.RH[k] += aux_up[i].area[k] * aux_up[i].RH[k] / a_bulk_c
                aux_bulk.buoy[k] += aux_up[i].area[k] * aux_up[i].buoy[k] / a_bulk_c
            end
        else
            aux_bulk.RH[k] = aux_en.RH[k]
            aux_bulk.T[k] = aux_en.T[k]
        end

        #####
        ##### update_GMV_diagnostics
        #####
        aux_gm.q_liq[k] = (aux_bulk.area[k] * aux_bulk.q_liq[k] + (1 - aux_bulk.area[k]) * aux_en.q_liq[k])
        aux_gm.q_ice[k] = (aux_bulk.area[k] * aux_bulk.q_ice[k] + (1 - aux_bulk.area[k]) * aux_en.q_ice[k])
        aux_gm.T[k] = (aux_bulk.area[k] * aux_bulk.T[k] + (1 - aux_bulk.area[k]) * aux_en.T[k])
        aux_gm.buoy[k] = (aux_bulk.area[k] * aux_bulk.buoy[k] + (1 - aux_bulk.area[k]) * aux_en.buoy[k])

        has_condensate = TD.has_condensate(aux_bulk.q_liq[k] + aux_bulk.q_ice[k])
        aux_bulk.cloud_fraction[k] = if has_condensate && a_bulk_c > 1e-3
            1
        else
            0
        end
    end

    #####
    ##### face variables: diagnose primitive, diagnose env and compute bulk
    #####
    # TODO: figure out why `ifelse` is allocating
    @inbounds for i in 1:N_up
        a_up_bcs = (; bottom = CCO.SetValue(edmf.area_surface_bc[i]), top = CCO.Extrapolate())
        If = CCO.InterpolateC2F(; a_up_bcs...)
        a_min = edmf.minimum_area
        a_up = aux_up[i].area
        @. aux_up_f[i].w = ifelse(If(a_up) >= a_min, max(prog_up_f[i].ρaw / (ρ0_f * If(a_up)), 0), FT(0))
    end
    @inbounds for i in 1:N_up
        aux_up_f[i].w[kf_surf] = edmf.w_surface_bc[i]
    end

    parent(aux_tc_f.bulk.w) .= 0
    a_bulk_bcs = (; bottom = CCO.SetValue(sum(edmf.area_surface_bc)), top = CCO.Extrapolate())
    Ifb = CCO.InterpolateC2F(; a_bulk_bcs...)
    @inbounds for i in 1:N_up
        a_up = aux_up[i].area
        a_up_bcs = (; bottom = CCO.SetValue(edmf.area_surface_bc[i]), top = CCO.Extrapolate())
        Ifu = CCO.InterpolateC2F(; a_up_bcs...)
        @. aux_tc_f.bulk.w += ifelse(Ifb(aux_bulk.area) > 0, Ifu(a_up) * aux_up_f[i].w / Ifb(aux_bulk.area), FT(0))
    end
    # Assuming gm.W = 0!
    @. aux_en_f.w = -Ifb(aux_bulk.area) / (1 - Ifb(aux_bulk.area)) * aux_tc_f.bulk.w

    #####
    #####  diagnose_GMV_moments
    #####
    get_GMV_CoVar(edmf, grid, state, :Hvar, :θ_liq_ice)
    get_GMV_CoVar(edmf, grid, state, :QTvar, :q_tot)
    get_GMV_CoVar(edmf, grid, state, :HQTcov, :θ_liq_ice, :q_tot)

    # TODO - use this inversion in free_convection_windspeed and not compute zi twice
    edmf.zi = get_inversion(grid, state, param_set, surface.Ri_bulk_crit)

    update_surface(Case, grid, state, gm, TS, param_set)
    update_forcing(Case, grid, state, gm, TS, param_set)
    update_radiation(Case, grid, state, gm, TS, param_set)

    compute_pressure_plume_spacing(edmf, param_set)

    #####
    ##### compute_updraft_closures
    #####

    Ic = CCO.InterpolateF2C()
    ∇c = CCO.DivergenceF2C()
    LB = CCO.LeftBiasedC2F(; bottom = CCO.SetValue(FT(0)))

    @inbounds for i in 1:N_up
        # compute ∇m at cell centers
        a_up = aux_up[i].area
        w_up = aux_up_f[i].w
        w_en = aux_en_f.w
        w_gm = prog_gm_f.w
        @. m_entr_detr = a_up * (Ic(w_up) - Ic(w_gm))
        @. ∇m_entr_detr = ∇c(wvec(LB(m_entr_detr)))
        @. w_up_c = Ic(w_up)
        @. w_en_c = Ic(w_en)
        @inbounds for k in real_center_indices(grid)
            # entrainment

            q_cond_up = TD.condensate(TD.PhasePartition(aux_up[i].q_tot[k], aux_up[i].q_liq[k], aux_up[i].q_ice[k]))
            q_cond_en = TD.condensate(TD.PhasePartition(aux_en.q_tot[k], aux_en.q_liq[k], aux_en.q_ice[k]))

            if aux_up[i].area[k] > 0.0

                εδ_model_vars = GeneralizedEntr(;
                    q_cond_up = q_cond_up,
                    q_cond_en = q_cond_en,
                    w_up = w_up_c[k],
                    w_en = w_en_c[k],
                    b_up = aux_up[i].buoy[k],
                    b_en = aux_en.buoy[k],
                    tke = aux_en.tke[k],
                    dMdz = ∇m_entr_detr[k],
                    M = m_entr_detr[k],
                    a_up = aux_up[i].area[k],
                    a_en = aux_en.area[k],
                    R_up = edmf.pressure_plume_spacing[i],
                    RH_up = aux_up[i].RH[k],
                    RH_en = aux_en.RH[k],
                    max_area = max_area,
                    zc_i = grid.zc[k].z,
                    Δt = Δt,
                )

                er = entr_detr(param_set, εδ_model_vars, edmf.entr_closure)
                aux_up[i].entr_sc[k] = er.ε_dyn
                aux_up[i].detr_sc[k] = er.δ_dyn
                # stochastic closure
                stoch_ε = stochastic_closure(param_set, sde_model, Entrainment())
                stoch_δ = stochastic_closure(param_set, sde_model, Detrainment())
                aux_up[i].entr_sc[k] *= stoch_ε
                aux_up[i].detr_sc[k] *= stoch_δ

                aux_up[i].frac_turb_entr[k] = er.ε_turb
                aux_up[i].horiz_K_eddy[k] = er.K_ε
            else
                aux_up[i].entr_sc[k] = 0.0
                aux_up[i].detr_sc[k] = 0.0
                aux_up[i].frac_turb_entr[k] = 0.0
                aux_up[i].horiz_K_eddy[k] = 0.0
            end
        end
    end

    @inbounds for i in 1:N_up

        # pressure
        b_up = aux_up[i].buoy
        a_up = aux_up[i].area
        w_up = aux_up_f[i].w
        w_en = aux_en_f.w
        updraft_top = up.updraft_top[i]
        asp_ratio = 1.0

        b_bcs = (; bottom = CCO.SetValue(b_up[kc_surf]), top = CCO.SetValue(b_up[kc_toa]))
        a_bcs = (; bottom = CCO.SetValue(edmf.area_surface_bc[i]), top = CCO.SetValue(FT(0)))
        w_bcs = (; bottom = CCO.SetValue(wvec(FT(0))), top = CCO.SetValue(wvec(FT(0))))
        bcs = (; a_up = a_bcs, b_up = b_bcs, w_up = w_bcs)

        nh_press_buoy = aux_up_f[i].nh_pressure_b
        nh_press_adv = aux_up_f[i].nh_pressure_adv
        nh_press_drag = aux_up_f[i].nh_pressure_drag

        nh_press_buoy .= nh_pressure_buoy(param_set, a_up, b_up, ρ0_f, asp_ratio, bcs)
        nh_press_adv .= nh_pressure_adv(param_set, updraft_top, a_up, ρ0_f, w_up, bcs)
        nh_press_drag .= nh_pressure_drag(param_set, updraft_top, a_up, ρ0_f, w_up, w_en, bcs)
        @. aux_up_f[i].nh_pressure = nh_press_buoy + nh_press_adv + nh_press_drag
    end

    #####
    ##### compute_eddy_diffusivities_tke
    #####

    # Subdomain exchange term
    Ic = CCO.InterpolateF2C()
    b_exch = center_aux_turbconv(state).b_exch
    parent(b_exch) .= 0
    a_en = aux_en.area
    w_en = aux_en_f.w
    tke_en = aux_en.tke
    @inbounds for i in 1:N_up
        a_up = aux_up[i].area
        w_up = aux_up_f[i].w
        δ_dyn = aux_up[i].detr_sc
        ε_turb = aux_up[i].frac_turb_entr
        @. b_exch +=
            a_up * Ic(w_up) * δ_dyn / a_en * (1 / 2 * (Ic(w_up) - Ic(w_en))^2 - tke_en) -
            a_up * Ic(w_up) * (Ic(w_up) - Ic(w_en)) * ε_turb * Ic(w_en) / a_en
    end

    Shear² = center_aux_turbconv(state).Shear²
    ∂qt∂z = center_aux_turbconv(state).∂qt∂z
    ∂θl∂z = center_aux_turbconv(state).∂θl∂z
    ∂θv∂z = center_aux_turbconv(state).∂θv∂z
    ∂qt∂z_sat = center_aux_turbconv(state).∂qt∂z_sat
    ∂θl∂z_sat = center_aux_turbconv(state).∂θl∂z_sat
    ∂θv∂z_unsat = center_aux_turbconv(state).∂θv∂z_unsat

    ∇0_bcs = (; bottom = CCO.Extrapolate(), top = CCO.Extrapolate())
    If0 = CCO.InterpolateC2F(; ∇0_bcs...)
    u_gm = prog_gm.u
    v_gm = prog_gm.v
    w_en = aux_en_f.w
    # compute shear
    @. Shear² = (∇c(wvec(If0(u_gm))))^2 + (∇c(wvec(If0(v_gm))))^2 + (∇c(wvec(w_en)))^2

    q_tot_en = aux_en.q_tot
    θ_liq_ice_en = aux_en.θ_liq_ice
    θ_virt_en = aux_en.θ_virt
    @. ∂qt∂z = ∇c(wvec(If0(q_tot_en)))
    @. ∂θl∂z = ∇c(wvec(If0(θ_liq_ice_en)))
    @. ∂θv∂z = ∇c(wvec(If0(θ_virt_en)))

    # Second order approximation: Use dry and cloudy environmental fields.
    cf = aux_en.cloud_fraction
    shm = copy(cf)
    parent(shm) .= shrink_mask(vec(cf))

    # Since NaN*0 ≠ 0, we need to conditionally replace
    # our gradients by their default values.
    @. ∂qt∂z_sat = ∇c(wvec(If0(aux_en_sat.q_tot)))
    @. ∂θl∂z_sat = ∇c(wvec(If0(aux_en_sat.θ_liq_ice)))
    @. ∂θv∂z_unsat = ∇c(wvec(If0(aux_en_unsat.θ_virt)))
    for k in real_center_indices(grid)
        if shm[k] == 0
            ∂qt∂z_sat[k] = ∂qt∂z[k]
            ∂θl∂z_sat[k] = ∂θl∂z[k]
            ∂θv∂z_unsat[k] = ∂θv∂z[k]
        end
    end

    @inbounds for k in real_center_indices(grid)

        # buoyancy_gradients
        if edmf.bg_closure == BuoyGradMean()
            # First order approximation: Use environmental mean fields.
            ts_en = TD.PhaseEquil_pTq(param_set, p0_c[k], aux_en.T[k], aux_en.q_tot[k])
            bg_kwargs = (;
                t_sat = aux_en.T[k],
                qv_sat = TD.vapor_specific_humidity(ts_en),
                qt_sat = aux_en.q_tot[k],
                θ_sat = aux_en.θ_dry[k],
                θ_liq_ice_sat = aux_en.θ_liq_ice[k],
                ∂θv∂z_unsat = ∂θv∂z[k],
                ∂qt∂z_sat = ∂qt∂z[k],
                ∂θl∂z_sat = ∂θl∂z[k],
                p0 = p0_c[k],
                en_cld_frac = aux_en.cloud_fraction[k],
                alpha0 = α0_c[k],
            )
            bg_model = EnvBuoyGrad(edmf.bg_closure; bg_kwargs...)

        elseif edmf.bg_closure == BuoyGradQuadratures()

            bg_kwargs = (;
                t_sat = aux_en_sat.T[k],
                qv_sat = aux_en_sat.q_vap[k],
                qt_sat = aux_en_sat.q_tot[k],
                θ_sat = aux_en_sat.θ_dry[k],
                θ_liq_ice_sat = aux_en_sat.θ_liq_ice[k],
                ∂θv∂z_unsat = ∂θv∂z_unsat[k],
                ∂qt∂z_sat = ∂qt∂z_sat[k],
                ∂θl∂z_sat = ∂θl∂z_sat[k],
                p0 = p0_c[k],
                en_cld_frac = aux_en.cloud_fraction[k],
                alpha0 = α0_c[k],
            )
            bg_model = EnvBuoyGrad(edmf.bg_closure; bg_kwargs...)
        end
        bg = buoyancy_gradients(param_set, bg_model)

        # Limiting stratification scale (Deardorff, 1976)
        # compute ∇Ri and Pr
        ∇_Ri = gradient_Richardson_number(bg.∂b∂z, Shear²[k], eps(0.0))
        aux_tc.prandtl_nvec[k] = turbulent_Prandtl_number(param_set, obukhov_length, ∇_Ri)

        ml_model = MinDisspLen(;
            z = grid.zc[k].z,
            obukhov_length = obukhov_length,
            tke_surf = aux_en.tke[kc_surf],
            ustar = surface.ustar,
            Pr = aux_tc.prandtl_nvec[k],
            p0 = p0_c[k],
            ∇b = bg,
            Shear² = Shear²[k],
            tke = aux_en.tke[k],
            b_exch = b_exch[k],
        )

        ml = mixing_length(param_set, ml_model)
        aux_tc.mls[k] = ml.min_len_ind
        aux_tc.mixing_length[k] = ml.mixing_length
        aux_tc.ml_ratio[k] = ml.ml_ratio

        KM[k] = c_m * ml.mixing_length * sqrt(max(aux_en.tke[k], 0.0))
        KH[k] = KM[k] / aux_tc.prandtl_nvec[k]

        aux_en_2m.tke.buoy[k] = -aux_en.area[k] * ρ0_c[k] * KH[k] * bg.∂b∂z
    end

    #####
    ##### compute covarinaces tendendies
    #####
    tke_press = aux_en_2m.tke.press
    w_en = aux_en_f.w
    parent(tke_press) .= 0
    @inbounds for i in 1:N_up
        w_up = aux_up_f[i].w
        nh_press = aux_up_f[i].nh_pressure
        @. tke_press += (Ic(w_en) - Ic(w_up)) * Ic(nh_press)
    end

    compute_covariance_entr(edmf, grid, state, Val(:tke), Val(:w))
    compute_covariance_entr(edmf, grid, state, Val(:Hvar), Val(:θ_liq_ice))
    compute_covariance_entr(edmf, grid, state, Val(:QTvar), Val(:q_tot))
    compute_covariance_entr(edmf, grid, state, Val(:HQTcov), Val(:θ_liq_ice), Val(:q_tot))
    compute_covariance_shear(edmf, grid, state, gm, Val(:tke), Val(:w))
    compute_covariance_shear(edmf, grid, state, gm, Val(:Hvar), Val(:θ_liq_ice))
    compute_covariance_shear(edmf, grid, state, gm, Val(:QTvar), Val(:q_tot))
    compute_covariance_shear(edmf, grid, state, gm, Val(:HQTcov), Val(:θ_liq_ice), Val(:q_tot))
    compute_covariance_dissipation(edmf, grid, state, :tke, param_set)
    compute_covariance_dissipation(edmf, grid, state, :Hvar, param_set)
    compute_covariance_dissipation(edmf, grid, state, :QTvar, param_set)
    compute_covariance_dissipation(edmf, grid, state, :HQTcov, param_set)

    # TODO defined again in compute_covariance_shear and compute_covaraince
    @inbounds for k in real_center_indices(grid)
        aux_en_2m.tke.rain_src[k] = 0
        aux_en_2m.Hvar.rain_src[k] = ρ0_c[k] * aux_en.area[k] * 2 * aux_en.Hvar_rain_dt[k]
        aux_en_2m.QTvar.rain_src[k] = ρ0_c[k] * aux_en.area[k] * 2 * aux_en.QTvar_rain_dt[k]
        aux_en_2m.HQTcov.rain_src[k] = ρ0_c[k] * aux_en.area[k] * aux_en.HQTcov_rain_dt[k]
    end

    get_GMV_CoVar(edmf, grid, state, :tke, :w)

    compute_diffusive_fluxes(edmf, grid, state, gm, Case, TS, param_set)

    # TODO: use dispatch
    if edmf.Precip.precipitation_model == "clima_1m"
        # helper to calculate the rain velocity
        # TODO: assuming gm.W = 0
        # TODO: verify translation
        term_vel_rain = aux_tc.term_vel_rain
        term_vel_snow = aux_tc.term_vel_snow
        prog_pr = center_prog_precipitation(state)

        @inbounds for k in real_center_indices(grid)
            term_vel_rain[k] = CM1.terminal_velocity(param_set, CM1.RainType(), ρ0_c[k], prog_pr.q_rai[k])
            term_vel_snow[k] = CM1.terminal_velocity(param_set, CM1.SnowType(), ρ0_c[k], prog_pr.q_sno[k])
        end
    end

    @inbounds for i in 1:N_up
        up.updraft_top[i] = 0.0

        @inbounds for k in real_center_indices(grid)
            if aux_up[i].area[k] > 1e-3
                up.updraft_top[i] = max(up.updraft_top[i], grid.zc[k])
            end
        end
    end

end


text:
using Test
using SparseArrays
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

@testset "bool.jl file Tests" begin

    n,m,p = 1,1,1
    V,(VV,EV,FV,CV) = Lar.cuboidGrid([n,m,p],true)
    cube = V,FV,EV

    threecubes = Lar.Struct([ cube,
        Lar.t(.3,.4,.25), Lar.r(pi/5,0,0), Lar.r(0,0,pi/12), cube,
        Lar.t(-.2,.4,-.2), Lar.r(0,pi/5,0), Lar.r(0,pi/12,0), cube ])

    V,FV,EV = Lar.struct2lar(threecubes)

        @testset "spaceindex Tests" begin
            @test typeof((V,FV)) == Tuple{Array{Float64,2},Array{Array{Int64,1},1}}
    		@test Lar.spaceindex([.5,.5,.5])((V,FV)) == [5,6,11,12,9]
			@test typeof(Lar.spaceindex([.5,.5,.5])((V,FV))) == Array{Int64,1}
            @test  Lar.spaceindex([0,0,0.])((V,FV)) == [1,19,3,5]
            @test  Lar.spaceindex([0,0,1.])((V,FV)) == [1,19,3,5]
        end

        @testset "rayintersection Tests" begin
            V,(VV,EV,FV,CV) = Lar.simplex(3,true);
    		@test Lar.rayintersection([.333,.333,0])(V,FV,4) == [0.333,0.333,0.3340000000000001]
			@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
        end

        @testset "planemap Tests" begin
    		@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
        end

        @testset "settestpoints Tests" begin
    		@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
        end

        @testset "testinternalpoint Tests" begin
    		@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
        end

        @testset "getinternalpoint Tests" begin
    		@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
        end

        @testset "chainbasis2solids Tests" begin
    		@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
        end

        @testset "internalpoints Tests" begin
    		@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
        end

        @testset "bool3d Tests" begin
    		@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
			@test bbbb == cccc
        end

end


text:
@sumProductRule(:node_type     => Beta,
                :outbound_type => Message{Beta},
                :inbound_types => (Nothing, Message{PointMass}, Message{PointMass}),
                :name          => SPBetaOutNMM)

@sumProductRule(:node_type     => Beta,
                :outbound_type => Message{Function},
                :inbound_types => (Message, Nothing, Message),
                :name          => SPBetaMNM)

@sumProductRule(:node_type     => Beta,
                :outbound_type => Message{Function},
                :inbound_types => (Message, Message, Nothing),
                :name          => SPBetaMMN)

@naiveVariationalRule(:node_type     => Beta,
                      :outbound_type => Message{Beta},
                      :inbound_types => (Nothing, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VBBetaOut)

@naiveVariationalRule(:node_type     => Beta,
                      :outbound_type => Message{Function},
                      :inbound_types => (ProbabilityDistribution, Nothing, ProbabilityDistribution),
                      :name          => VBBetaA)

@naiveVariationalRule(:node_type     => Beta,
                      :outbound_type => Message{Function},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, Nothing),
                      :name          => VBBetaB)

mutable struct SPBetaOutMCNMM <: SumProductRule{Beta} end
outboundType(::Type{SPBetaOutMCNMM}) = Message{SampleList}
function isApplicable(::Type{SPBetaOutMCNMM}, input_types::Vector{Type})
  nothing_inputs = 0
  point_inputs = 0

  for input_type in input_types
      if input_type == Nothing
          nothing_inputs += 1
      elseif matches(input_type, Message{PointMass})
          point_inputs += 1
      end
  end

  return (nothing_inputs == 1) && (point_inputs != 2)
end


text:
if basename(pwd()) == "aoc"
    cd("2015/16")
end

ForgottenInt{T} = Union{Int,Missing}

struct Compound
    children::ForgottenInt
    cats::ForgottenInt
    samoyeds::ForgottenInt
    pomeranians::ForgottenInt
    akitas::ForgottenInt
    vizslas::ForgottenInt
    goldfish::ForgottenInt
    trees::ForgottenInt
    cars::ForgottenInt
    perfumes::ForgottenInt
end
Compound(;
    children::ForgottenInt = missing,
    cats::ForgottenInt = missing,
    samoyeds::ForgottenInt = missing,
    pomeranians::ForgottenInt = missing,
    akitas::ForgottenInt = missing,
    vizslas::ForgottenInt = missing,
    goldfish::ForgottenInt = missing,
    trees::ForgottenInt = missing,
    cars::ForgottenInt = missing,
    perfumes::ForgottenInt = missing
) = Compound(
    children,
    cats,
    samoyeds,
    pomeranians,
    akitas,
    vizslas,
    goldfish,
    trees,
    cars,
    perfumes
)

function Compound(str::AbstractString)
    Compound(; (
        Symbol(m[1]) => parse(Int, m[2]) for m in eachmatch(r"(\S+): (\d+)", str)
    )...)
end

loadsues() = [Compound(line[findfirst(':', line):end]) for line in eachline("input.txt")]

tickertape = Compound("""
    children: 3
    cats: 7
    samoyeds: 2
    pomeranians: 3
    akitas: 0
    vizslas: 0
    goldfish: 5
    trees: 3
    cars: 2
    perfumes: 1
""")

function comparefields(a, b, fields, fn = (==))
    all(
        begin
            af = getfield(a, f)
            bf = getfield(b, f)
            ismissing(af) || !ismissing(bf) && fn(af, bf)
        end
        for f in fields
    )
end

Base.issubset(a::Compound, b::Compound) = comparefields(a, b, fieldnames(Compound))
Base.issubset(c::Compound) = Base.Fix2(issubset, c)

# part 1
findfirst(issubset(tickertape), sues)

function isweirdsubset(a::Compound, b::Compound)
    comparefields(a, b, (:children, :samoyeds, :akitas, :vizslas, :cars, :perfumes)) &&
        comparefields(a, b, (:cats, :trees), (>)) &&
        comparefields(a, b, (:pomeranians, :goldfish), (<))
end
isweirdsubset(c::Compound) = Base.Fix2(isweirdsubset, c)

# part 2
findfirst(isweirdsubset(tickertape), sues)

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
function rbc()
    ∞ = Inf
    @variables α, β, ρ, δ, σ
    @variables t::Integer, k(..), z(..), c(..), q(..)
    x = [k, z]
    y = [c, q]
    p = [α, β, ρ, δ, σ]

    H = [1 / c(t) - (β / c(t + 1)) * (α * exp(z(t + 1)) * k(t + 1)^(α - 1) + (1 - δ)),
         c(t) + k(t + 1) - (1 - δ) * k(t) - q(t), q(t) - exp(z(t)) * k(t)^α,
         z(t + 1) - ρ * z(t)]

    steady_states = [k(∞) ~ (((1 / β) - 1 + δ) / α)^(1 / (α - 1)), z(∞) ~ 0,
                     c(∞) ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1)) -
                            δ * (((1 / β) - 1 + δ) / α)^(1 / (α - 1)),
                     q(∞) ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1))]

    steady_states_iv = [k(∞) ~ (((1 / β) - 1 + δ) / α)^(1 / (α - 1)), z(∞) ~ 0,
                        c(∞) ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1)) -
                               δ * (((1 / β) - 1 + δ) / α)^(1 / (α - 1)),
                        q(∞) ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1))]

    n_ϵ = 1
    n_x = length(x)
    n_y = length(y)
    n_p = length(p)
    Γ = reshape([σ], n_ϵ, n_ϵ)
    η = reshape([0; -1], n_x, n_ϵ) # η is n_x * n_ϵ matrix

    return H, (; t, x, y, p, steady_states, steady_states_iv, Γ, η), "rbc"
end

function rbc_multiple_shocks()
    H, nt = rbc()
    σ = nt.Γ[1, 1]
    Γ = [σ 0; σ σ]
    η = [0 -1; 0 -1]
    return H, merge(nt, (; Γ, η)), "rbc_multiple_shocks"
end

function rbc_reorder_ss()
    H, nt = rbc()
    steady_states = reverse(nt.steady_states)
    return H, merge(nt, (; steady_states)), "rbc_reorder_ss"
end

function rbc_solve_steady_state()
    H, nt = rbc()
    steady_states = nothing
    return H, merge(nt, (; steady_states)), "rbc_solve_steady_state"
end

function rbc_solve_steady_state_different_iv()
    H, nt = rbc()

    ∞ = Inf
    @variables α, β, ρ, δ, σ
    @variables t::Integer, k(..), z(..), c(..), q(..)
    x = [k, z]
    y = [c, q]
    p = [α, β, ρ, δ, σ]

    steady_states_iv = [
        k(∞) ~ 0.1,# i.e., not (((1 / β) - 1 + δ) / α)^(1 / (α - 1)),
        z(∞) ~ 0.01,
        c(∞) ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1)) -
                δ * (((1 / β) - 1 + δ) / α)^(1 / (α - 1)),
        q(∞) ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1)),
    ]

    steady_states = nothing

    return H, merge(nt, (; steady_states, steady_states_iv)),
           "rbc_solve_steady_state_different_iv"
end

# Add observation equation with a single source of measurement error
function rbc_observables()
    H, nt = rbc()

    n_x = length(nt.x)
    n_y = length(nt.y)
    n_z = 2 # number of observables
    Q = zeros(n_z, n_x + n_y)
    Q[1, 1] = 1.0
    Q[2, 3] = 1.0

    @variables Ω_1
    Ω = [Ω_1, Ω_1]
    p = [nt.p; Ω_1]
    return H, merge(nt, (; p, Ω, Q)), "rbc_observables"
end

# Add observation equation with two parameters in the source of measurement error
function rbc_observables_separate_variance()
    H, nt = rbc()

    n_x = length(nt.x)
    n_y = length(nt.y)
    n_z = 2 # number of observables
    Q = zeros(n_z, n_x + n_y)
    Q[1, 1] = 1.0
    Q[2, 3] = 1.0

    @variables Ω_1, Ω_2
    Ω = [Ω_1, Ω_2]
    p = [nt.p; Ω_1; Ω_2]
    return H, merge(nt, (; p, Ω, Q)), "rbc_observables_separate_variance"
end


text:
"""
    SimulationLog(; value_dict::Dict{Symbol,Vector}=Dict{Symbol,Vector}(), is_active::Bool=false)

A log of saved variables from a simulation. Variables can be accessed with dot notation or the `getproperty`
function.
"""
Base.@kwdef mutable struct SimulationLog
    value_dict::Dict{Symbol, Array} = Dict{Symbol, Array}()
    is_active::Bool = false
end

function Base.show(io::IO, log::SimulationLog)
    println(io, "SimulationLog with signals:")
    for (name, val) in value_dict(log)
        println(io, "  $name :: $(typeof(val))")
    end
end

# Field access
"""
    value_dict(log::SimulationLog)

Get the value dictionary of a `SimulationLog`.
"""
value_dict(log::SimulationLog) = getfield(log, :value_dict)

"""
    is_active(log::SimulationLog)

See if simulation is currently active (able to be logged to).
"""
is_active(log::SimulationLog) = getfield(log, :is_active)

# Setting fields
"""
    activate!(log::SimulationLog=GLOBAL_LOG)

Activate a `SimulationLog`. If no log is given, activate the `GLOBAL_LOG`.
"""
activate!(log::SimulationLog=GLOBAL_LOG) = setfield!(log, :is_active, true)

"""
    deactivate!(log::SimulationLog=GLOBAL_LOG)

Deactivate a `SimulationLog`. If no log is given, deactivate the `GLOBAL_LOG`.
"""
deactivate!(log::SimulationLog=GLOBAL_LOG) = setfield!(log, :is_active, false)

"""
    reset!(log::SimulationLog=GLOBAL_LOG)

Reset the value dictionary of a `SimulationLog`. If no log is given, the `GLOBAL_LOG` will
be reset
"""
reset!(log::SimulationLog=GLOBAL_LOG) = setfield!(log, :value_dict, Dict{Symbol, Array}())

# Value access
Base.keys(log::SimulationLog) = keys(value_dict(log))
Base.propertynames(log::SimulationLog) = Tuple(keys(value_dict(log)))

Base.getproperty(log::SimulationLog, s::Symbol) = value_dict(log)[s]
function Base.setproperty!(log::SimulationLog, s::Symbol, val)
    val_dict = value_dict(log)
    val_dict[s] = val
end

Base.getindex(log::SimulationLog, s::Symbol) = getproperty(log, s)
Base.setindex!(log::SimulationLog, val, s::Symbol) = setproperty!(log, s, val)

"""
    struct Logged{T,N,A,S<:AbstractTimeseriesSolution{T,N,A}} <: AbstractTimeseriesSolution{T,N,A}

    Logged(sol::AbstractTimeseriesSolution, log::SimulationLog)

Logged differential equation solution. All properties of the underlying solution
can be accessed with the `getproperty` or 'dot' accessing. The inner `SimulationLog` can
be accessed with `.log`. A `Logged` solution can be used the same way as its underlying
solution (i.e. indexed with `sol[i]` and interpolated with `sol(t)`)
"""
struct Logged{T,N,A,S<:AbstractTimeseriesSolution{T,N,A}} <: AbstractTimeseriesSolution{T,N,A}
    sol::S
    log::SimulationLog
end

function Base.show(io::IO, mime::MIME"text/plain", sol::Logged)
    show(io, mime, sol.sol)
    println(io, "\nWith logged values:")
    for (name, val) in value_dict(sol.log)
        println(io, " $name: $(join(string.(size(val)), "×"))-element $(typeof(val))")
    end
end

(sol::Logged)(args...; kwargs...) = sol.sol(args...; kwargs...)

Base.getindex(sol::Logged, i::Int) = getindex(sol.sol, i)
Base.getindex(sol::Logged, i::Colon) = getindex(sol.sol, i)

@inline Base.getproperty(sol::Logged, s::Symbol) = getproperty(sol, Val(s))
@inline Base.getproperty(sol::Logged, ::Val{:log}) = getfield(sol, :log)
@inline Base.getproperty(sol::Logged, ::Val{:sol}) = getfield(sol, :sol)
@inline Base.getproperty(sol::Logged, ::Val{s}) where {s} = getproperty(sol.sol, s)

Base.propertynames(sol::Logged) = (:log, propertynames(sol.sol)...)

max_stars_repo_name: jonniedie/SimulationLogs

max_issues_repo_name: jonniedie/SimulationLogs

max_forks_repo_name: jonniedie/SimulationLogs


text:
D = [
	-5  -9  -9   7   0;
	 7   0   4   0   9;
	 8   0   5   4  -6
]

# when data is representad as a vector
Y = [
    -1;
     0;
     5
]

X = matching_pursuit(Y, D)
@test norm(Y-D*X) < 1e-6

# when data is representad as a matrix
Y = [
    -1  2;
     0  3;
     5 -5
]

X = matching_pursuit(Y, D)
@test norm(Y-D*X) < 1e-6

# when max_iter < 0
@test_throws ArgumentError matching_pursuit(Y, D; max_iter = 0)
@test_throws ArgumentError matching_pursuit(Y, D; max_iter = -1)

# when tolerance < 0
@test_throws ArgumentError matching_pursuit(Y, D, tolerance = 0.)
@test_throws ArgumentError matching_pursuit(Y, D; tolerance = -1.0)

# when dimensions of data and atoms do not match

Y = [
     2;
    -5
]

@test_throws ArgumentError matching_pursuit(Y, D)


text:
function sais{T}(text::Vector{T}, sa::Vector{Int}, fs::Int, n::Int, k::Int, isbwt::Bool)
    pidx = 0
    flags = 0
    if k <= 256
        C = zeros(Int, k)
        if k <= fs
            B = offset_to_array(sa, n+fs-k+1)
            flags = 1
        else
            B = zeros(Int, k)
            flags = 3
        end
    elseif k <= fs
        C = offset_to_array(sa, n+fs-k+1)
        if k <= fs-k
            B = offset_to_array(sa, n+fs-2k+1)
            flags = 0
        elseif k <= 1024
            B = zeros(Int, k)
            flags = 2
        else
            B = C
            flags = 8
        end
    else
        C = B = zeros(Int, k)
        flags = 4 | 8
    end

    # stage 1: reduce the problem by at least 1/2
    # sort all the LMS-substrings
    setcounts!(text, C, n, k)
    getbuckets(C, B, k, true)
    for i = 1:n
        sa[i] = 0
    end
    b = -1
    i = j = n
    m = 0
    c0 = c1 = text[n]
    i -= 1
    while 1 <= i && ((c0 = text[i]) >= c1)
        c1 = c0
        i -= 1
    end
    while 1 <= i
        c1 = c0
        i -= 1
        while 1 <= i && ((c0 = text[i]) <= c1)
            c1 = c0
            i -= 1
        end
        if 1 <= i
            0 <= b && (sa[b+1] = j)
            b = (B[c1+1] -= 1)
            j = i-1
            m += 1
            c1 = c0
            i -= 1
            while 1 <= i && ((c0 = text[i]) >= c1)
                c1 = c0
                i -= 1
            end
        end
    end
    if 1 < m
        LMSsort(text, sa, C, B, n, k)
        name = LMSpostproc(text, sa, n, m)
    elseif m == 1
        sa[b+1] = j + 1
        name = 1
    else
        name = 0
    end

    # stage 2
    if name < m
        newfs = n + fs - 2m
        if flags & (1 | 4 | 8) == 0
            if (k + name) <= newfs
                newfs -= k
            else
                flags |= 8
            end
        end
        j = 2m + newfs
        for i = (m + (n >> 1)):-1:(m+1)
            if sa[i] != 0
                sa[j] = sa[i] - 1
                j -= 1
            end
        end
        RA = offset_to_array(sa, m+newfs+1)
        sais(RA, sa, newfs, m, name, false)

        i = n
        j = 2m
        c0 = c1 = text[n]
        while 1 <= (i -= 1) && ((c0 = text[i]) >= c1)
            c1 = c0
        end
        while 1 <= i
            c1 = c0
            while 1 <= (i -= 1) && ((c0 = text[i]) <= c1)
                c1 = c0
            end
            if 1 <= i
                sa[j] = i
                j -= 1
                c1 = c0
                while 1 <= (i -= 1) && ((c0 = text[i]) >= c1)
                    c1 = c0
                end
            end
        end
        for i = 1:m
            sa[i] = sa[m+sa[i]+1]
        end
        if flags & 4 != 0
            C = B = zeros(Int, k)
        end
        if flags & 2 != 0
            B = zeros(Int, k)
        end
    end

    # stage 3
    flags & 8 != 0 && setcounts!(text, C, n, k)
    if 1 < m
        getbuckets(C, B, k, true)
        i = m - 1
        j = n
        p = sa[m]
        c1 = text[p+1]
        while true
            c0 = c1
            q = B[c0+1]
            while q < j
                j -= 1
                sa[j+1] = 0
            end
            while true
                j -= 1
                sa[j+1] = p
                i -= 1
                i < 0 && break
                p = sa[i+1]
                c1 = text[p+1]
                c1 != c0 && break
            end
            i < 0 && break
        end
        while 0 < j
            j -= 1
            sa[j+1] = 0
        end
    end
    #if (isbwt == false) { induceSA(T, SA, C, B, n, k); }
    #else { pidx = computeBWT(T, SA, C, B, n, k); }
    isbwt == false || throw("isbwt is true.")
    induceSA(text, sa, C, B, n, k)
    sa+1
end
sais(text::Vector) = sais(text, Array{Int}(length(text)), 0, length(text), 65536, false)

function setcounts!{T}(text::Vector{T}, C::Vector{Int}, n::Int, k::Int)
    for i = 1:k
        C[i] = 0
    end
    for i = 1:n
        C[text[i]+1] += 1
    end
end

function getbuckets(C::Vector{Int}, B::Vector{Int}, k::Int, isend::Bool)
    sum = 0
    if isend
        for i = 1:k
            sum += C[i]
            B[i] = sum
        end
    else
        for i = 1:k
            sum += C[i]
            B[i] = sum - C[i]
        end
    end
end

offset_to_array(x::Vector, i::Int) = unsafe_wrap(Array, pointer(x,i), length(x)-i+1)

function LMSsort{T<:Integer}(text::Vector{T}, sa, C, B, n::Int, k::Int)
    C == B && setcounts!(text, C, n, k)
    getbuckets(C, B, k, false)
    j = n - 1
    c1 = text[j+1]
    b = B[c1+1]
    j -= 1
    sa[b+1] = text[j+1] < c1 ? ~j : j
    b += 1
    for i = 1:n
        if 0 < (j = sa[i])
            if (c0 = text[j+1]) != c1
                B[c1+1] = b
                c1 = c0
                b = B[c1+1]
            end
            j -= 1
            sa[b+1] = text[j+1] < c1 ? ~j : j
            b += 1
            sa[i] = 0
        elseif j < 0
            sa[i] = ~j
        end
    end
    C == B && setcounts!(text, C, n, k)
    getbuckets(C, B, k, true)
    c1 = 0
    b = B[c1+1]
    for i = n:-1:1
        if 0 < (j = sa[i])
            c0 = text[j+1]
            if c0 != c1
                B[c1+1] = b
                c1 = c0
                b = B[c1+1]
            end
            j -= 1
            b -= 1
            sa[b+1] = text[j+1] > c1 ? ~(j+1) : j
            sa[i] = 0
        end
    end
end

function LMSpostproc{T}(text::Vector{T}, sa::Vector{Int}, n::Int, m::Int)
    i = 1
    while (p = sa[i]) < 0
        sa[i] = ~p
        i += 1
    end
    if i-1 < m
        j = i
        i += 1
        while true
            if (p = sa[i]) < 0
                sa[j] = ~p
                j += 1
                sa[i] = 0
                j-1 == m && break
            end
            i += 1
        end
    end

    i = j = n
    c0 = c1 = text[n]
    while 1 <= (i -= 1) && ((c0 = text[i]) >= c1)
        c1 = c0
    end
    while 1 <= i
        c1 = c0
        while 1 <= (i -= 1) && ((c0 = text[i]) <= c1)
            c1 = c0
        end
        if 1 <= i
            sa[m + (i >> 1) + 1] = j - i
            j = i + 1
            c1 = c0
            while 1 <= (i -= 1) && ((c0 = text[i]) >= c1)
                c1 = c0
            end
        end
    end
    name = 0
    q = n
    qlen = 0
    for i = 1:m
        p = sa[i]
        plen = sa[m + (p >> 1) + 1]
        diff = true
        if plen == qlen && (q + plen < n)
            j = 0
            while j < plen && text[p+j+1] == text[q+j+1]
                j += 1
            end
            j == plen && (diff = false)
        end
        if diff != false
            name += 1
            q = p
            qlen = plen
        end
        sa[m + (p >> 1) + 1] = name
    end
    return name
end

function induceSA{T}(text::Vector{T}, sa::Vector{Int}, C::Vector{Int}, B::Vector{Int}, n::Int, k::Int)
    C == B && setcounts!(text, C, n, k)
    getbuckets(C, B, k, false)
    j = n - 1
    c1 = text[j+1]
    b = B[c1+1]
    sa[b+1] = 0 < j && text[j] < c1 ? ~j : j
    b += 1
    for i = 1:n
        j = sa[i]
        sa[i] = ~j
        if 0 < j
            j -= 1
            if (c0 = text[j+1]) != c1
                B[c1+1] = b
                c1 = c0
                b = B[c1+1]
            end
            sa[b+1] = 0 < j && text[j] < c1 ? ~j : j
            b += 1
        end
    end
    C == B && setcounts!(text, C, n, k)
    getbuckets(C, B, k, true)
    c1 = 0
    b = B[c1+1]
    for i = n:-1:1
        if 0 < (j = sa[i])
            j -= 1
            if (c0 = text[j+1]) != c1
                B[c1+1] = b
                c1 = c0
                b = B[c1+1]
            end
            b -= 1
            sa[b+1] = j == 0 || text[j] > c1 ? ~j : j
        else
            sa[i] = ~j
        end
    end
end


text:
using Pkg
Pkg.activate(@__DIR__)

using ProgressiveHedging
import JuMP
import Ipopt

function two_stage_model(scenario_id::ScenarioID)

    model = JuMP.Model(()->Ipopt.Optimizer())
    JuMP.set_optimizer_attribute(model, "print_level", 0)
    JuMP.set_optimizer_attribute(model, "tol", 1e-12)
    JuMP.set_optimizer_attribute(model, "acceptable_tol", 1e-12)

    scen = value(scenario_id)

    ref = JuMP.@variable(model, x >= 0.0)
    stage1 = [ref]

    ref = JuMP.@variable(model, y >= 0.0)
    stage2 = [ref]

    b_s = scen == 0 ? 11.0 : 4.0
	c_s = scen == 0 ? 0.5 : 10.0

    JuMP.@constraint(model, x + y == b_s)

    JuMP.@objective(model, Min, 1.0*x + c_s*y)

    return JuMPSubproblem(model,
	                  scenario_id,
                          Dict(stid(1) => stage1,
			       stid(2) => stage2)
                          )
end

scen_tree = two_stage_tree(2)
(niter, abs_res, rel_res, obj, soln_df, phd) = solve(scen_tree,
                                                     two_stage_model,
                                                     ScalarPenaltyParameter(1.0)
                                                     )
@show niter
@show abs_res
@show rel_res
@show obj
@show soln_df

dual_df = retrieve_w(phd)
@show dual_df

raw_values_df = retrieve_no_hats(phd)
@show raw_values_df


text: ### v0.8 Deprecations


text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_lowered)" --startup-file=no %s 2>&1 | FileCheck %s

sle_int(x, y) = Base.sle_int(x, y)
emit(sle_int, Int64, Int64)

# CHECK: module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.sle_int), Int64, Int64}"(%arg0: !jlir<"typeof(Main.sle_int)">, %arg1: i64, %arg2: i64) -> i1 attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     %0 = cmpi sle, %arg1, %arg2 : i64
# CHECK-NEXT:     return %0 : i1
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK:   llvm.func @"Tuple{typeof(Main.sle_int), Int64, Int64}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: i64, %arg2: i64) -> i1 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.icmp "sle" %arg1, %arg2 : i64
# CHECK-NEXT:     llvm.return %0 : i1
# CHECK-NEXT:   }
# CHECK-NEXT:   llvm.func @"_mlir_ciface_Tuple{typeof(Main.sle_int), Int64, Int64}"(%arg0: !llvm.ptr<struct<"struct_jl_value_type", opaque>>, %arg1: i64, %arg2: i64) -> i1 attributes {llvm.emit_c_interface, sym_visibility = "nested"} {
# CHECK-NEXT:     %0 = llvm.call @"Tuple{typeof(Main.sle_int), Int64, Int64}"(%arg0, %arg1, %arg2) : (!llvm.ptr<struct<"struct_jl_value_type", opaque>>, i64, i64) -> i1
# CHECK-NEXT:     llvm.return %0 : i1
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: femtomc/brutus

max_issues_repo_name: femtomc/brutus

max_forks_repo_name: femtomc/brutus


text:
# using Pandas
using Dates

"""
Append table or row to CSV file.

julia> using LarSurf
julia> LarSurf.add_to_csv("file.csv", (:first=>[4, 5], :second=>["a", "b"]))
"""
function add_to_csv(df::Tuple, data_file::AbstractString)
    dfo = Dict(df)
    add_to_csv(dfo, data_file)
end

function add_to_csv(df::Dict, data_file::AbstractString)
    # wrap values into array
    one_row = false
    for (key, val) in df
        if !isa(val, Array)
            one_row = true
            break
        end
    end
    if one_row
        add_row_to_csv(df, data_file)
    else
        for i=1:length(df[collect(keys(df))[1]])
            new_row_dict = Dict()
            for (key, val) in df
                new_row_dict[key] = val[i]
            end
            add_row_to_csv(new_row_dict, data_file)

        end
    end
end

# function add_to_csv(df::Dict, data_file::AbstractString)
#     # wrap values into array
#     wrap_with_array = false
#     for (key, val) in df
#         if !isa(val, Array)
#             wrap_with_array = true
#             break
#         end
#     end
#
#     new_dict = Dict()
#     for (key, val) in df
#         if wrap_with_array
#             new_dict[key] = [val]
#         else
#             new_dict[key] = val
#         end
#     end
#
#
#     dfo = Pandas.DataFrame(new_dict)
#     add_to_csv(dfo, data_file)
#
# end

# function add_to_csv(df::Pandas.DataFrame, data_file::AbstractString)
#
#     if isfile(data_file)
#         df0 = read_csv(data_file)
#         df = concat((df0, df); ignore_index=true, sort=false)
#     end
#     to_csv(df, data_file; index=false) #, sep=";")
# end

"""
Add description of segmentation to dictionary
"""
function segmentation_description_to_dict!(dct, segmentation, data_name="data_")
    n_zeros = sum(segmentation.==0)
    n_non_zeros = sum(segmentation.!=0)
    dct[data_name * "size_1"] = size(segmentation)[1]
    dct[data_name * "size_2"] = size(segmentation)[2]
    if length(size(segmentation)) > 2
        dct[data_name * "size_3"] = size(segmentation)[3]
    end
    dct[data_name * "zeros_number"] = n_zeros
    dct[data_name * "non_zeros_number"] = n_non_zeros
    dct[data_name * "voxel_number"] = prod(size(segmentation))
    return dct
end

"""
experiment is string used to tag experiment
"""
function timed_and_more_to_dict!(dct, timed; experiment=nothing, datetime=true)
    dct["elapsed"] = timed[2]
    dct["alocated"] = timed[3]
    if experiment != nothing
        dct["experiment"] = experiment
    end
    if datetime
        dct[prefix * "datetime"] = Dates.format(Dates.now(), "yyyy-mm-dd HH:MM:SS")
    end
    return dct
end

"""
experiment is string used to tag experiment
"""
function timed_to_dict!(dct, timed, prefix="")
    dct[prefix * "elapsed"] = timed[2]
    dct[prefix * "alocated"] = timed[3]
    return dct
end

function datetime_to_dict!(dct, key="datetime")
    dct["datetime"] = Dates.format(Dates.now(), "yyyy-mm-dd HH:MM:SS")
    return dct
end

function size_to_dict!(dct, data_size, data_name="")
    dct[data_name * "size_1"] = data_size[1]
    dct[data_name * "size_2"] = data_size[2]
    dct[data_name * "size_3"] = data_size[3]
    return dct
end


text:
module LotkaVolterra2dProblem

    using GeometricIntegrators.Equations

    export lotka_volterra_2d_ode,  lotka_volterra_2d_pdae,
           lotka_volterra_2d_iode, lotka_volterra_2d_idae,
           lotka_volterra_2d_vode, lotka_volterra_2d_vdae,
           lotka_volterra_2d_dg, lotka_volterra_2d_vdae_slrk

    Δt = 0.01
    nt = 10

    const A1=1.0
    const A2=1.0
    const B1=1.0
    const B2=2.0

    const X0=1.0
    const Y0=1.0

    function ϑ₁(t, q)
        q[2] + log(q[2]) / q[1]
    end

    function ϑ₂(t, q)
        q[1]
    end

    function v₁(t, q)
        q[1] * (A2*q[2] - B2)
    end

    function v₂(t, q)
        q[2] * (B1 - A1*q[1])
    end

    const q₀=[X0, Y0]
    const v₀=[v₁(0, q₀), v₂(0, q₀)]

    function dϑ₁d₁(t, q)
        - log(q[2]) / q[1]^2
    end

    function dϑ₁d₂(t, q)
        1 + 1 / (q[1] * q[2])
    end

    function dϑ₂d₁(t, q)
        one(eltype(q))
    end

    function dϑ₂d₂(t, q)
        zero(eltype(q))
    end

    function d²ϑ₁d₁d₁(t, q)
        + 2 * log(q[2]) / q[1]^3
    end

    function d²ϑ₁d₁d₂(t, q)
        - 1 / (q[1]^2 * q[2])
    end

    function d²ϑ₁d₂d₁(t, q)
        - 1 / (q[1]^2 * q[2])
    end

    function d²ϑ₁d₂d₂(t, q)
        - 1 / (q[1] * q[2]^2)
    end

    function d²ϑ₂d₁d₁(t, q)
        zero(eltype(q))
    end

    function d²ϑ₂d₁d₂(t, q)
        zero(eltype(q))
    end

    function d²ϑ₂d₂d₁(t, q)
        zero(eltype(q))
    end

    function d²ϑ₂d₂d₂(t, q)
        zero(eltype(q))
    end

    function ϑ(t::Number, q::AbstractVector, Θ::AbstractVector)
        Θ[1] = ϑ₁(t,q)
        Θ[2] = ϑ₂(t,q)
        nothing
    end

    function ϑ(t::Number, q::AbstractVector)
        [ϑ₁(t,q), ϑ₂(t,q)]
    end

    function ω(t, q, Ω)
        Ω[1,1] = 0
        Ω[1,2] = dϑ₁d₂(t,q) - dϑ₂d₁(t,q)

        Ω[2,1] = dϑ₂d₁(t,q) - dϑ₁d₂(t,q)
        Ω[2,2] = 0

        nothing
    end

    function f₁(t, q, v)
        dϑ₁d₁(t,q) * v[1] + dϑ₂d₁(t,q) * v[2]
    end

    function f₂(t, q, v)
        dϑ₁d₂(t,q) * v[1] + dϑ₂d₂(t,q) * v[2]
    end

    function g₁(t, q, v)
        dϑ₁d₁(t,q) * v[1] + dϑ₁d₂(t,q) * v[2]
    end

    function g₂(t, q, v)
        dϑ₂d₁(t,q) * v[1] + dϑ₂d₂(t,q) * v[2]
    end

    function g̅₁(t, q, v)
        d²ϑ₁d₁d₁(t,q) * q[1] * v[1] + d²ϑ₁d₂d₁(t,q) * q[1] * v[2] + d²ϑ₂d₁d₁(t,q) * q[2] * v[1] + d²ϑ₂d₂d₁(t,q) * q[2] * v[2]
    end

    function g̅₂(t, q, v)
        d²ϑ₁d₁d₂(t,q) * q[1] * v[1] + d²ϑ₁d₂d₂(t,q) * q[1] * v[2] + d²ϑ₂d₁d₂(t,q) * q[2] * v[1] + d²ϑ₂d₂d₂(t,q) * q[2] * v[2]
    end

    function hamiltonian(t, q)
        A1*q[1] + A2*q[2] - B1*log(q[1]) - B2*log(q[2])
    end

    function dHd₁(t, q)
        A1 - B1 / q[1]
    end

    function dHd₂(t, q)
        A2 - B2 / q[2]
    end

    function lotka_volterra_2d_dH(t, q, dH)
        dH[1] = dHd₁(t, q)
        dH[2] = dHd₂(t, q)
        nothing
    end

    function lotka_volterra_2d_ϑ(t, q, Θ)
        ϑ(t, q, Θ)
    end

    function lotka_volterra_2d_ϑ(t, q, v, Θ)
        ϑ(t, q, Θ)
    end

    function lotka_volterra_2d_ϑ(κ, t, q, v, Θ)
        Θ[1] = (1-κ) * ϑ₁(t,q) - κ * f₁(t,q,q)
        Θ[2] = (1-κ) * ϑ₂(t,q) - κ * f₂(t,q,q)
        nothing
    end

    function lotka_volterra_2d_ω(t, q, Ω)
        ω(t, q, Ω)
    end

    function lotka_volterra_2d_v(t, q, v)
        v[1] = v₁(t, q)
        v[2] = v₂(t, q)
        nothing
    end

    function lotka_volterra_2d_v(t, q, p, v)
        lotka_volterra_2d_v(t, q, v)
    end

    function lotka_volterra_2d_v_ham(t, q, p, v)
        v[1] = 0
        v[2] = 0
        nothing
    end

    function lotka_volterra_2d_f(t::Real, q::Vector, v::Vector, f::Vector)
        f[1] = f₁(t,q,v) - dHd₁(t,q)
        f[2] = f₂(t,q,v) - dHd₂(t,q)
        nothing
    end

    function lotka_volterra_2d_f(κ::Real, t::Real, q::Vector, v::Vector, f::Vector)
        f[1] = (1-κ) * f₁(t,q,v) - κ * (g₁(t,q,v) + g̅₁(t,q,v)) - dHd₁(t,q)
        f[2] = (1-κ) * f₂(t,q,v) - κ * (g₂(t,q,v) + g̅₂(t,q,v)) - dHd₂(t,q)
        nothing
    end

    function lotka_volterra_2d_f_ham(t::Real, q::Vector, f::Vector)
        f[1] = - dHd₁(t,q)
        f[2] = - dHd₂(t,q)
        nothing
    end

    function lotka_volterra_2d_f_ham(t::Real, q::Vector, v::Vector, f::Vector)
        lotka_volterra_2d_f_ham(t, q, f)
    end

    function lotka_volterra_2d_g(t::Real, q::Vector, v::Vector, g::Vector)
        g[1] = f₁(t,q,v)
        g[2] = f₂(t,q,v)
        nothing
    end

    function lotka_volterra_2d_g(t::Real, q::Vector, p::Vector, v::Vector, g::Vector)
        lotka_volterra_2d_g(t, q, v, g)
    end

    function lotka_volterra_2d_g(κ::Real, t::Real, q::Vector, v::Vector, g::Vector)
        g[1] = (1-κ) * f₁(t,q,v) - κ * (g₁(t,q,v) + g̅₁(t,q,v))
        g[2] = (1-κ) * f₂(t,q,v) - κ * (g₂(t,q,v) + g̅₂(t,q,v))
        nothing
    end

    # function lotka_volterra_2d_g(κ::Real, t::Real, q::Vector, v::Vector, g::Vector)
    #     g[1] = (1-κ) * g₁(t,q,v) - κ * g̅₁(t,q,v) - κ * f₁(t,q,v)
    #     g[2] = (1-κ) * g₂(t,q,v) - κ * g̅₂(t,q,v) - κ * f₂(t,q,v)
    #     nothing
    # end

    function lotka_volterra_2d_g̅(t::Real, q::Vector, v::Vector, g::Vector)
        g[1] = g₁(t,q,v)
        g[2] = g₂(t,q,v)
        nothing
    end

    function lotka_volterra_2d_g̅(t::Real, q::Vector, p::Vector, v::Vector, g::Vector)
        lotka_volterra_2d_g̅(t, q, v, g)
    end

    function lotka_volterra_2d_u(t, q, v, u)
        u[1] = v[1]
        u[2] = v[2]
        nothing
    end

    function lotka_volterra_2d_u(t, q, p, v, u)
        lotka_volterra_2d_u(t, q, v, u)
    end

    function lotka_volterra_2d_ϕ(t, q, p, ϕ)
        ϕ[1] = p[1] - ϑ₁(t,q)
        ϕ[2] = p[2] - ϑ₂(t,q)
        nothing
    end

    function lotka_volterra_2d_ψ(t, q, p, v, f, ψ)
        ψ[1] = f[1] - g₁(t,q,v)
        ψ[2] = f[2] - g₂(t,q,v)
        nothing
    end

    function lotka_volterra_2d_ode(q₀=q₀)
        ODE(lotka_volterra_2d_v, q₀; h=hamiltonian)
    end

    function lotka_volterra_2d_iode(q₀=q₀, p₀=ϑ(0, q₀))
        IODE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f,
             lotka_volterra_2d_g, q₀, p₀;
             h=hamiltonian, v=lotka_volterra_2d_v)
    end

    function lotka_volterra_2d_vode(q₀=q₀, p₀=ϑ(0, q₀))
        VODE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f,
             lotka_volterra_2d_g, q₀, p₀;
             v=lotka_volterra_2d_v, Ω=lotka_volterra_2d_ω, ∇H=lotka_volterra_2d_dH)
    end

    function lotka_volterra_2d_idae(q₀=q₀, p₀=ϑ(0, q₀), λ₀=zero(q₀))
        IDAE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f,
             lotka_volterra_2d_u, lotka_volterra_2d_g,
             lotka_volterra_2d_ϕ, q₀, p₀, λ₀;
             v=lotka_volterra_2d_v)
    end

    function lotka_volterra_2d_pdae(q₀=q₀, p₀=ϑ(0, q₀), λ₀=zero(q₀))
        PDAE(lotka_volterra_2d_v_ham, lotka_volterra_2d_f_ham,
             lotka_volterra_2d_u, lotka_volterra_2d_g,
             lotka_volterra_2d_ϕ, q₀, p₀, λ₀)
    end

    function lotka_volterra_2d_vdae(q₀=q₀, p₀=ϑ(0, q₀), λ₀=zero(q₀))
        VDAE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f_ham,
             lotka_volterra_2d_g, lotka_volterra_2d_g̅,
             lotka_volterra_2d_ϕ, lotka_volterra_2d_ψ,
             q₀, p₀, λ₀; v=lotka_volterra_2d_v)
    end

    function lotka_volterra_2d_vdae_slrk(q₀=q₀, p₀=ϑ(0, q₀), λ₀=zero(q₀))
        VDAE(lotka_volterra_2d_ϑ, lotka_volterra_2d_f,
             lotka_volterra_2d_g, lotka_volterra_2d_g̅,
             lotka_volterra_2d_ϕ, lotka_volterra_2d_ψ,
             q₀, p₀, λ₀; v=lotka_volterra_2d_v)
    end

    function lotka_volterra_2d_dg(q₀=q₀; κ=0)
        lotka_volterra_2d_ϑ_κ(t, q, v, p) = lotka_volterra_2d_ϑ(κ, t, q, v, p)
        lotka_volterra_2d_f_κ(t, q, v, f) = lotka_volterra_2d_f(κ, t, q, v, f)
        lotka_volterra_2d_g_κ(t, q, λ, g) = lotka_volterra_2d_g(κ, t, q, λ, g)

        IODE(lotka_volterra_2d_ϑ_κ, lotka_volterra_2d_f_κ,
             lotka_volterra_2d_g_κ, q₀, p₀;
             v=lotka_volterra_2d_v)
    end

end


text:
using LightXML

struct AnimatedTile
    name::String
    path::String

    origX::Number
    origY::Number

    posX::Number
    posY::Number

    delay::Number
end

struct Mask
    mask::Array{Int, 2}
end

# Convert mask string into a 3x3 matrix
function Mask(s::String)
    data = replace(replace(s, "-" => ""), "x" => "2")

    return Mask(transpose(reshape(parse.(Int, (collect(data))), (3, 3))))
end

function sortScore(m::Mask)
    return length(filter(v -> v == 2, m.mask))
end

# Get quads for a given mask
function convertTileString(s::String)
    parts = split(s, ";")
    res = Tuple{Integer, Integer}[]

    for part in parts
        x, y = parse.(Int, split(part, ","))

        push!(res, (x, y))
    end

    return res
end

function loadAnimatedTilesXML(fn::String)
    xdoc = parse_file(fn)
    xroot = root(xdoc)

    animatedTiles = AnimatedTile[]

    for elem in child_elements(xroot)
        push!(animatedTiles, AnimatedTile(
            attribute(elem, "name"),
            attribute(elem, "path"),

            parseNumber(attribute(elem, "origX")),
            parseNumber(attribute(elem, "origY")),

            parseNumber(attribute(elem, "posX")),
            parseNumber(attribute(elem, "posY")),

            parseNumber(attribute(elem, "delay"))
        ))
    end

    return animatedTiles
end

function loadTilesetXML(fn::String)
    xdoc = parse_file(fn)
    xroot = root(xdoc)

    paths = Dict{Char, String}()
    masks = Dict{Char, Any}()
    padding = Dict{Char, Array{Tuple{Integer, Integer}, 1}}()
    center = Dict{Char, Array{Tuple{Integer, Integer}, 1}}()
    ignores = Dict{Char, Array{Char, 1}}()

    orderedElements = sort(collect(child_elements(xroot)), by=set -> attribute(set, "copy") !== nothing)

    for set in orderedElements
        id = attribute(set, "id")[1] # Char doesn't work on 1 length string?
        path = attribute(set, "path")
        copyTileset = attribute(set, "copy")
        ignore = attribute(set, "ignores")

        paths[id] = "tilesets/$path"

        if ignore !== nothing
            ignores[id] = collect(ignore)
        end

        if copyTileset !== nothing
            padding[id] = deepcopy(padding[copyTileset[1]])
            center[id] = deepcopy(center[copyTileset[1]])
            masks[id] = deepcopy(masks[copyTileset[1]])
        end

        currMasks = Tuple{Mask, Array{Tuple{Integer, Integer}, 1}, String}[]

        for c in child_elements(set)
            mask = attribute(c, "mask")
            tilesString = attribute(c, "tiles")

            if mask == "padding"
                padding[id] = convertTileString(tilesString)

            elseif mask == "center"
                center[id] = convertTileString(tilesString)

            else
                sprites = attribute(c, "sprites")
                push!(currMasks, (
                    Mask(mask),
                    convertTileString(tilesString),
                    isa(sprites, String) ? sprites : ""
                ))
            end
        end

        if !isempty(currMasks)
            masks[id] = currMasks
            sort!(masks[id], by=m -> sortScore(m[1]), rev=false)
        end
    end

    return paths, masks, padding, center, ignores
end

struct TilerMeta
    paths::Dict{Char, String}
    masks::Dict{Char, Any}
    paddings::Dict{Char, Array{Tuple{Integer, Integer}, 1}}
    centers::Dict{Char, Array{Tuple{Integer, Integer}, 1}}
    ignores::Dict{Char, Array{Char, 1}}
end

TilerMeta(s::String) = TilerMeta(loadTilesetXML(s)...)

fgTilerMeta = nothing
bgTilerMeta = nothing
animatedTilesMeta = nothing

function loadFgTilerMeta(side::Union{Side, Nothing}, fn::String, loadDefault::Bool=false)
    path = joinpath(storageDirectory, "XML", "ForegroundTiles.xml")

    if side !== nothing && !loadDefault
        meta = get(Dict{String, Any}, side.data, "meta")
        metaPath = get(meta, "ForegroundTiles", "")

        hasRoot, modRoot = getModRoot(fn)
        xmlPath = joinpath(modRoot, metaPath)

        if !isempty(metaPath) && hasRoot
            path = xmlPath
        end
    end

    try
        global fgTilerMeta = TilerMeta(path)

    catch e
        if !loadDefault
            println(Base.stderr, "Failed to load custom ForegroundTiles XML")
            println(Base.stderr, e)

            loadFgTilerMeta(side, fn, true)
        end
    end
end

function loadBgTilerMeta(side::Union{Side, Nothing}, fn::String, loadDefault::Bool=false)
    path = joinpath(storageDirectory, "XML", "BackgroundTiles.xml")

    if side !== nothing && !loadDefault
        meta = get(Dict{String, Any}, side.data, "meta")
        metaPath = get(meta, "BackgroundTiles", "")

        hasRoot, modRoot = getModRoot(fn)
        xmlPath = joinpath(modRoot, metaPath)

        if !isempty(metaPath) && hasRoot
            path = xmlPath
        end
    end

    try
        global bgTilerMeta = TilerMeta(path)

    catch e
        if !loadDefault
            println(Base.stderr, "Failed to load custom BackgroundTiles XML")
            println(Base.stderr, e)

            loadBgTilerMeta(side, fn, true)
        end
    end
end

function loadAnimatedTilesMeta(side::Union{Side, Nothing}, fn::String, loadDefault::Bool=false)
    path = joinpath(storageDirectory, "XML", "AnimatedTiles.xml")

    if side !== nothing && !loadDefault
        meta = get(Dict{String, Any}, side.data, "meta")
        metaPath = get(meta, "AnimatedTiles", "")

        hasRoot, modRoot = getModRoot(fn)
        xmlPath = joinpath(modRoot, metaPath)

        if !isempty(metaPath) && hasRoot
            path = xmlPath
        end
    end

    try
        global animatedTilesMeta = loadAnimatedTilesXML(path)

    catch e
        if !loadDefault
            println(Base.stderr, "Failed to load custom AnimatedTiles XML")
            println(Base.stderr, e)

            loadAnimatedTilesMeta(side, fn, true)
        end
    end
end

function loadXMLMeta()
    side = loadedState.side
    filename = loadedState.filename

    loadFgTilerMeta(side, filename)
    loadBgTilerMeta(side, filename)
    loadAnimatedTilesMeta(side, filename)
end

function getTile(tiles::Tiles, x::Integer, y::Integer)
    return get(tiles.data, (y, x), ' ')
end

function checkPadding(tiles::Tiles, x::Integer, y::Integer)
    # Checks for '0' even though getTile might return ' '
    # This is due to quirks/special flags in the actuall autotiler
    return getTile(tiles, x - 2, y) == '0' || getTile(tiles, x + 2, y) == '0' || getTile(tiles, x, y - 2) == '0' || getTile(tiles, x, y + 2) == '0'
end

# Return the "mask" value for a tile
function checkTile(value::Char, target::Char, ignore::Array{Char, 1}=Char[])
    return !(target == '0' || target in ignore || ('*' in ignore && value != target))
end

function checkMask(adj::Array{Bool, 2}, mask::Mask)
    mask = mask.mask
    for i in 1:9
        if mask[i] != 2
            if adj[i] != mask[i]
                return false
            end
        end
    end

    return true
end

function getMaskQuads(x::Integer, y::Integer, tiles::Tiles, meta::TilerMeta)
    value = tiles.data[y, x]

    masks = get(meta.masks, value, Tuple{Mask, Array{Tuple{Integer, Integer}, 1}, String}[])
    ignore = get(meta.ignores, value, Char[])

    adjacent = get(tiles.data, (y - 1:y + 1, x - 1:x + 1), value)
    adjacent = Bool[checkTile(value, target, ignore) for target in adjacent]

    for data in masks
        mask, quads, sprites = data

        if checkMask(adjacent, mask)
            return quads, sprites
        end
    end

    if checkPadding(tiles, x, y)
        return get(meta.paddings, value, Tuple{Integer, Integer}[(0, 0)]), ""

    else
        return get(meta.centers, value, Tuple{Integer, Integer}[(0, 0)]), ""
    end

    return Tuple{Integer, Integer}[(5, 12)], ""
end


text:
using LazyGraphs
using Documenter

DocMeta.setdocmeta!(LazyGraphs, :DocTestSetup, :(using LazyGraphs); recursive=true)

makedocs(;
    modules=[LazyGraphs],
    authors="Cédric Simal, University of Namur",
    repo="https://github.com/csimal/LazyGraphs.jl/blob/{commit}{path}#{line}",
    sitename="LazyGraphs.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://csimal.github.io/LazyGraphs.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/csimal/LazyGraphs.jl",
    devbranch="main",
)


text:
function LinearObjectiveFunction(arg0::RealVector, arg1::jdouble)
    return LinearObjectiveFunction((RealVector, jdouble), arg0, arg1)
end

function LinearObjectiveFunction(arg0::Vector{jdouble}, arg1::jdouble)
    return LinearObjectiveFunction((Vector{jdouble}, jdouble), arg0, arg1)
end

function equals(obj::LinearObjectiveFunction, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_coefficients(obj::LinearObjectiveFunction)
    return jcall(obj, "getCoefficients", RealVector, ())
end

function get_constant_term(obj::LinearObjectiveFunction)
    return jcall(obj, "getConstantTerm", jdouble, ())
end

function hash_code(obj::LinearObjectiveFunction)
    return jcall(obj, "hashCode", jint, ())
end

function value(obj::LinearObjectiveFunction, arg0::RealVector)
    return jcall(obj, "value", jdouble, (RealVector,), arg0)
end

function value(obj::LinearObjectiveFunction, arg0::Vector{jdouble})
    return jcall(obj, "value", jdouble, (Vector{jdouble},), arg0)
end