text:
cd(@__DIR__)
sbmlfile = joinpath("data", "reactionsystem_01.xml")
@parameters t, k1, c1
@variables s1(t), s2(t), s1s2(t)

kinetic_params = Dict{String,Tuple{Float64,String}}()

COMP1 = SBML.Compartment("c1", true, 3, 2.0, "nl", nothing, nothing)
SPECIES1 = SBML.Species("s1", "c1", false, nothing, nothing, (1.0, "substance"), nothing, true, nothing, nothing)  # Todo: Maybe not support units in initial_concentration?
SPECIES2 = SBML.Species("s2", "c1", false, nothing, nothing, nothing, (1.0, "substance/nl"), false, nothing, nothing)
KINETICMATH1 = SBML.MathIdent("k1")
KINETICMATH2 = SBML.MathApply("*", SBML.Math[
    SBML.MathIdent("k1"), SBML.MathIdent("s2")])
KINETICMATH3 = SBML.MathApply("-", SBML.Math[KINETICMATH2, KINETICMATH1])
REACTION1 = SBML.Reaction(Dict(), Dict("s1" => 1), kinetic_params, "", "", nothing, KINETICMATH1, false, nothing, nothing)
REACTION2 = SBML.Reaction(Dict("s2" => 1), Dict(), kinetic_params, "", "", nothing, KINETICMATH2, false, nothing, nothing)
REACTION3 = SBML.Reaction(Dict("s2" => 1), Dict(), kinetic_params, "", "", nothing, KINETICMATH3, true, nothing, nothing)
MODEL1 = SBML.Model(
    parameters = Dict("k1" => (1.0, "")),
    compartments = Dict("c1" => COMP1),
    species = Dict("s1" => SPECIES1),
    reactions = Dict("r1" => REACTION1),
)  # PL: For instance in the compartments dict, we may want to enforce that key and compartment.name are identical
MODEL2 = SBML.Model(
    parameters = Dict("k1" => (1.0, "")),
    compartments = Dict("c1" => COMP1),
    species = Dict("s2" => SPECIES2),
    reactions = Dict("r2" => REACTION2),
)
MODEL3 = SBML.Model(
    parameters = Dict("k1" => (1.0, "")),
    compartments = Dict("c1" => COMP1),
    species = Dict("s2" => SPECIES2),
    reactions = Dict("r3" => REACTION3),
)
# MODEL2 = SBML.Model(Dict("k1" => ("", 1.0)), Dict(), Dict("c1" => COMP1), Dict("s2" => SPECIES2), Dict("r2" => REACTION2), Dict(), Dict())
# MODEL3 = SBML.Model(Dict("k1" => ("", 1.0)), Dict(), Dict("c1" => COMP1), Dict("s2" => SPECIES2), Dict("r3" => REACTION3), Dict(), Dict())

# Test ReactionSystem constructor
rs = ReactionSystem(MODEL1)
@test isequal(Catalyst.get_eqs(rs), Catalyst.Reaction[Catalyst.Reaction(k1, nothing, [s1], nothing, [1.0])])
@test isequal(Catalyst.get_iv(rs), t)
@test isequal(Catalyst.get_states(rs), [s1])
@test isequal(Catalyst.get_ps(rs), [k1, c1])
@named rs = ReactionSystem(MODEL1)
isequal(nameof(rs), :rs)

rs = ReactionSystem(readSBML(sbmlfile))
@test isequal(Catalyst.get_eqs(rs), Catalyst.Reaction[Catalyst.Reaction(0.25c1 * k1, [s1, s2], [s1s2], [1.0, 1.0], [1.0])])
@test isequal(Catalyst.get_iv(rs), t)
@test isequal(Catalyst.get_states(rs), [s1, s1s2, s2])
@test isequal(Catalyst.get_ps(rs), [k1, c1])
@named rs = ReactionSystem(MODEL1)
isequal(nameof(rs), :rs)

rs = ReactionSystem(MODEL3)  # Contains reversible reaction
@test isequal(Catalyst.get_eqs(rs),
    Catalyst.Reaction[
        Catalyst.Reaction(0.5k1, [s2], nothing,
            [1.0], nothing),
        Catalyst.Reaction(k1, nothing, [s2],
            nothing, [1.0])])
@test isequal(Catalyst.get_iv(rs), t)
@test isequal(Catalyst.get_states(rs), [s2])
@test isequal(Catalyst.get_ps(rs), [k1, c1])

@test_nowarn convert(ModelingToolkit.ODESystem, rs)
@test_nowarn structural_simplify(convert(ModelingToolkit.ODESystem, rs))

# Test ODESystem constructor
odesys = ODESystem(MODEL1)
trueeqs = Equation[Differential(t)(s1)~k1]
@test isequal(Catalyst.get_eqs(odesys), trueeqs)
@test isequal(Catalyst.get_iv(odesys), t)
@test isequal(Catalyst.get_states(odesys), [s1])
@test isequal(Catalyst.get_ps(odesys), [k1, c1])
u0 = [s1 => 1.0]
par = [k1 => 1.0, c1 => 2.0]
@test isequal(ModelingToolkit.defaults(odesys), ModelingToolkit._merge(u0, par))  # PL: @Anand: for some reason this does not work with `Catalyst.get_default()`
@named odesys = ODESystem(MODEL1)
isequal(nameof(odesys), :odesys)
@test_nowarn structural_simplify(odesys)

odesys = ODESystem(readSBML(sbmlfile))
m = readSBML(sbmlfile)
trueeqs = Equation[Differential(t)(s1)~-0.25c1*k1*s1*s2,
    Differential(t)(s1s2)~0.25c1*k1*s1*s2,
    Differential(t)(s2)~-0.25c1*k1*s1*s2]
@test isequal(Catalyst.get_eqs(odesys), trueeqs)
@test isequal(Catalyst.get_iv(odesys), t)
@test isequal(Catalyst.get_states(odesys), [s1, s1s2, s2])
@test isequal(Catalyst.get_ps(odesys), [k1, c1])
u0 = [s1 => 2 * 1.0, s2 => 2 * 1.0, s1s2 => 2 * 1.0]
par = [k1 => 1.0, c1 => 2.0]
@test isequal(ModelingToolkit.defaults(odesys), ModelingToolkit._merge(u0, par))
@named odesys = ODESystem(MODEL1)
isequal(nameof(odesys), :odesys)

@test_nowarn ODEProblem(odesys, [], [0.0, 1.0], [])

# Test ODEProblem
oprob = ODEProblem(ODESystem(MODEL1), [], [0.0, 1.0], [])
sol = solve(oprob, Tsit5())
@test isapprox(sol.u, [[1.0], [2.0]])

@test_nowarn ODEProblem(ODESystem(readSBML(sbmlfile)), [], [0.0, 1.0], [])

# Test checksupport
# @test_nowarn SBMLToolkit.checksupport(MODEL1)
# sbc = SBML.Species("sbc", "c1", true, nothing, nothing, (1., "substance"), nothing, true)
# mod = SBML.Model(Dict("k1" => 1.), Dict(), Dict("c1" => COMP1), Dict("sbc" => sbc), Dict("r1" => REACTION1), Dict(), Dict())
# @test_logs (:warn, "Species sbc has `boundaryCondition` or is `constant`. This will lead to wrong results when simulating the `ReactionSystem`.") SBMLToolkit.checksupport(mod)

# Test checksupport
@test_nowarn SBMLToolkit.checksupport(sbmlfile)
@test_throws ErrorException SBMLToolkit.checksupport(joinpath("data", "unsupported.sbml"))

# Test _get_substitutions
truesubs = Dict(Num(Symbolics.variable(:c1; T = Real)) => c1,
    Num(Symbolics.variable(:k1; T = Real)) => k1,
    Num(Symbolics.variable(:s1; T = Real)) => s1)
subs = SBMLToolkit._get_substitutions(MODEL1)
@test isequal(subs, truesubs)

# Test mtk_reactions
reaction = SBMLToolkit.mtk_reactions(MODEL1)[1]
truereaction = Catalyst.Reaction(k1, nothing, [s1], nothing, [1])  # Todo: implement Sam's suggestion on mass action kinetics
@test isequal(reaction, truereaction)

km = SBML.MathTime("x")
reac = SBML.Reaction(Dict("s1" => 1), Dict(), kinetic_params, "", "", nothing, km, false, nothing, nothing)
mod = SBML.Model(
    parameters = Dict("k1" => (1.0, "")),
    compartments = Dict("c1" => COMP1),
    species = Dict("s1" => SPECIES1),
    reactions = Dict("r1" => reac),
)

@test isequal(Catalyst.DEFAULT_IV, SBMLToolkit.mtk_reactions(mod)[1].rate)

# test from_noncombinatoric
@test isequal(4k1, SBMLToolkit.from_noncombinatoric(k1, [2, 2], false))
@test isequal(k1, SBMLToolkit.from_noncombinatoric(k1, nothing, false))

# Test use_rate
@test isequal(SBMLToolkit.use_rate(k1 * s1, [s1], [1]), (k1, false))  # Case hOSU=true
@test isequal(SBMLToolkit.use_rate(k1 * s1 * s2 / (c1 + s2), [s1], [1]), (k1 * s1 * s2 / (c1 + s2), true))  # Case Michaelis-Menten kinetics

# Test getreagents
@test isequal((nothing, [s1], nothing, [1.0]), SBMLToolkit.getreagents(REACTION1.reactants, REACTION1.products, MODEL1))

constspec = SBML.Species("constspec", "c1", true, nothing, nothing, (1.0, "substance"), nothing, true, nothing, nothing)
ncs = SBMLToolkit.create_var("constspec", Catalyst.DEFAULT_IV)
kineticmath = SBML.MathApply("-", SBML.Math[
    SBML.MathApply("*", SBML.Math[
        SBML.MathIdent("k1"),
        SBML.MathIdent("constspec")]),
    SBML.MathIdent("k1")])
constreac = SBML.Reaction(Dict("constspec" => 1), Dict(), kinetic_params, "", "", nothing, kineticmath, false, nothing, nothing)

constmod = SBML.Model(
    parameters = Dict("k1" => (1.0, "")),
    compartments = Dict("c1" => COMP1),
    species = Dict("constspec" => constspec),
    reactions = Dict("r1" => constreac),
)
@test isequal(([ncs], [ncs], [1.0], [1.0]), SBMLToolkit.getreagents(constreac.reactants, constreac.products, constmod))
@test isequal((nothing, nothing, nothing, nothing), SBMLToolkit.getreagents(constreac.reactants, constreac.products, constmod; rev = true))

# Test getunidirectionalcomponents
km = SBML.MathApply("-", SBML.Math[KINETICMATH1, SBML.MathIdent("c1")])
sm = convert(Num, km)
kl = SBMLToolkit.getunidirectionalcomponents(sm)
@test isequal(kl, (k1, c1))

km = SBML.MathApply("-", SBML.Math[KINETICMATH1, KINETICMATH2])
sm = convert(Num, km)
fw, rv = SBMLToolkit.getunidirectionalcomponents(sm)
rv = substitute(rv, Dict(SBMLToolkit.create_var("s2") => SBMLToolkit.create_var("s2", Catalyst.DEFAULT_IV)))
@test isequal((fw, rv), (k1, k1 * s2))

km = SBML.MathIdent("s1s2")
sm1 = convert(Num, km)
sm2 = sm - sm1
@test_throws ErrorException SBMLToolkit.getunidirectionalcomponents(sm2)
@test_throws ErrorException SBMLToolkit.getunidirectionalcomponents(:k1)

# Test get_paramap
trueparamap = [k1 => 1.0, c1 => 2.0]
paramap = SBMLToolkit.get_paramap(MODEL1)
@test isequal(paramap, trueparamap)

# Test getmassaction
@test isequal(SBMLToolkit.getmassaction(k1 * s1, [s1], [1]), k1)  # Case hOSU=true
@test isequal(SBMLToolkit.getmassaction(k1 * s1 / c1, [s1], [1]), k1 / c1)  # Case hOSU=false
@test isequal(SBMLToolkit.getmassaction(k1 + c1, nothing, nothing), k1 + c1)  # Case zero order kineticLaw
@test isnan(SBMLToolkit.getmassaction(k1 * s1 * s2 / (c1 + s2), [s1], [1]))  # Case Michaelis-Menten kinetics
@test isnan(SBMLToolkit.getmassaction(k1 * s1 * Catalyst.DEFAULT_IV, [s1], [1]))  # Case kineticLaw with time

@test isnan(SBMLToolkit.getmassaction(k1 * s1 * s2, [s1], [1]))
@test isnan(SBMLToolkit.getmassaction(k1 + c1, [s1], [1]))
@test_throws ErrorException SBMLToolkit.getmassaction(k1, nothing, [1])

# default test
@test ModelingToolkit.defaults(m) == ModelingToolkit.defaults(ReactionSystem(m))


text:
using Test;
using LinearAlgebra;
using SparseArrays;
using ParamLevelSet;
using Random;
n = 10;
B = getSpMatBuilder(Int64,Float64,n, n, 30);

for k=1:10
	setNext!(B,k,k,1.0);
end
II = getSparseMatrix(B);
@test norm(II - SparseMatrixCSC(1.0I,n,n)) < 1e-14

B.V.*=100.0;
reset!(B);

for i=1:10
	JJ = randperm(n)[1:5];
	for j = 1:5
		setNext!(B,i,JJ[j],1.0);
	end
end
@test sum(B.V[:]) == 50;
II = getSparseMatrix(B);
IIT = getSparseMatrixTransposed(B);
@test norm(II - IIT') < 1e-14

ext: jl


text:
baremodule BenchmarkConfigSweeps

# Config specifires
function env end
function nthreads end
function julia end

# Sweep execution
function run end

# Data loading
function load end
function simpletable end
function flattable end

module Internal

using ..BenchmarkConfigSweeps: BenchmarkConfigSweeps

import JSON3
import StructTypes
import Tables
using Accessors: @set
using ArgCheck: @argcheck
using AutoHashEquals: @auto_hash_equals
using BenchmarkTools: BenchmarkTools, BenchmarkGroup

include("utils.jl")
include("configurations.jl")
include("run.jl")
include("info.jl")
include("data.jl")

end  # module Internal

Internal.define_docstring()

end  # baremodule BenchmarkConfigSweeps


text:
@testset "simple broadcasting" begin

    b = rand(2)
    c = rand(3)
    bc = rand(2,3)

    ## compare to ordinary broadcasting
    Z = zeros(3,2)
    A = bc' .+ b' .* c .^ 3
    @cast B[j,i] := bc[i,j] .+ b[i] .* c[j]^3
    @cast Z[j,i] = bc[i,j] .+ b[i] .* c[j]^3
    @test A ==B == Z

    Z1 = zeros(3,2,1)
    @cast Z1[j,i,1] = bc[i,j] .+ b[i] .* c[j]^3
    @test A == Z1[:,:,1]

    A = @. exp(-bc + 2*c') / sqrt(b)
    @cast B[i,j] := exp(-bc[i,j] + 2*c[j]) / sqrt(b[i])
    @test A == B

    ## test permutedims in orient function
    bcdef = rand(2,3,4,5,6);
    @cast A[d,e,c,f,b] := 10 + bcdef[b,c,d,e,f];
    @test A[3,4,2,5,1] == 10 + bcdef[1,2,3,4,5]

    ## check that f(A) is allowed as array name
    exp_all(A) = similar(A) .= exp.(A)
    A = exp_all(bc) .* sqrt.(b .+ 3)
    @cast B[i,j] := exp_all(bc)[i,j] * sqrt(b[i] + 3)
    @test A == B

    ## check that log(2) doesn't get broadcast
    global log_n = 0
    two=2.0
    log_count(x) = (global log_n += 1; log(x))
    A = log.(bc) ./ log_count(2)
    @cast B[i,j] := log(bc[i,j]) / log_count(2)
    @cast C[i,j] := log(bc[i,j]) / log_count(two)
    @test A == B == C
    @test log_n == 3

end
@testset "broadcast reduction" begin

    b = rand(2)
    c = rand(3)
    bc = rand(2,3)

    Z = zeros(3)
    A = sum( bc' .+ b' .* c ; dims=2) |> vec
    @reduce B[j] := sum(i) bc[i,j] .+ b[i] .* c[j]
    @reduce Z[j] = sum(i)  bc[i,j] .+ b[i] .* c[j]
    @test A ==B == Z

    ## complete reduction
    @reduce S[] := sum(i,j) bc[j,i] + 1
    T = @reduce sum(i,j) bc[j,i] + 1
    @test S[] ≈ sum(bc .+ 1) ≈ T

    ## with output shaping / fixing
    bcde = rand(2,3,4,5);
    @reduce A[d⊗c,_,b] := sum(e) bcde[-b,c,d,e]
    B = similar(A);
    @reduce B[d⊗c,_,b] = sum(e) bcde[-b,c,d,e]
    @test A ≈ B

    @reduce A[_,d⊗b,_] := sum(e,c) bcde[b,c,d,e]
    B = similar(A);
    @reduce B[_,d⊗b,_] = sum(c,e) bcde[b,c,d,e]
    @test A ≈ B

end
@testset "complex" begin

    C = rand(3,3) .+ im .* rand(3,3)
    R = rand(1:33, 3,3)

    @cast A[i,j] := R[i,j] + C[j,i]'
    @test A == R .+ C'

    @cast A[i,j] = R[j,i] + C[i,j]' + 10im
    @test A == R' .+ conj.(C) .+ 10im

    @cast B[i,j] := R[i,j] + C'[i,j]
    @test B == R .+ C'

end
@testset "diag" begin

    B = rand(4,5)
    @test TensorCast.diagview(B) == TensorCast.diag(B) # really LinearAlgebra!

    R = rand(2,2)
    @cast V[i] := R[i,i]
    @cast V2[i] := R[i,i] lazy=false

    @test V[1] == R[1,1]
    @test V[2] == R[2,2]
    @test V == V2

    # R4 = rand(1:10, 4)
    # @cast W[i] := R4[i⊗i]^2 # diag(reshape(R4, (sz_i, sz_i))) with sz_i = (:)
    # @test W[1] := R4[1]^2
    # @test W[2] := R4[4]^2

    M = [repeat([10i+j],2) for i=1:3, j=1:3]
    @cast A[i,k] := M[i,i][k]

    @test A[1,1] == A[1,2] == 11
    @test A[3,1] == A[3,2] == 33

end
@testset "mapslices etc" begin

    M = rand(-10:10, 3,4)
    g(x) = x ./ sum(abs,x)
    @cast N[i,j] := g(M[:,j])[i]
    N2 = mapslices(g, M, dims=1)
    @test N == N2

    # V = [ rand(1:99, 2) for _=1:5, _=1:1 ]
    # @cast W[i,j] := (v->v .// 2)(V[j,_])[i] # TODO
    # W2 = hcat(V...) .// 2
    # @test W == W2

    fun(m::AbstractMatrix, a=0) = vec(sum(m,dims=1)) .+ a
    X = rand(2,3,4)
    Y = rand(4)
    @cast Z[i,j] := fun(X[:,i,:],42)[j] + Y[j]^2
    Z2 = dropdims(sum(X, dims=1), dims=1) .+ 42 .+ transpose(Y).^2
    @test Z ≈ Z2

    A = rand(2:2:10, 3,2,4)
    ff(v::AbstractVector) = v .^2
    B = rand(-1000:1000:1000, 2,1,4)
    CM = rand(2,3)
    C = eachcol(CM) |> collect
    @cast Y[i,j,k] := ff(A[j,:,k])[i] + B[i,_,k] + C[j][i]
    Y2 = permutedims(A, (2,1,3)) .^2 .+ B .+ CM
    @test Y ≈ Y2

end


text:
using Roots
using Statistics
using FileIO
using Parameters
include("random_matrices.jl")
include("network_sims_gpu.jl")

resultsdir = "results2"

etas = collect(-1.0f0:0.1f0:1.0f0)
ind = parse(Int64, ARGS[1])

N = 3000
eta = etas[ind]
s = 0.0f0

J_seeds = collect(5000:5019)
h0_seed = 6000

h0 = generate_nv(N, Float32, seed=h0_seed)

wait = 200.0f0
t_bound = 2000.0f0
dt = 0.01f0

startstep = round(Int, wait / dt) + 1

@unpack target_avg_var = load("vanilla_target_avg_var.jld2")

# Difference between the average variance of the network ensemble and target
function f(g)
    network_avg_vars = Vector{Float32}(undef, length(J_seeds))
    for (i, J_seed) in enumerate(J_seeds)
        J = generate_J(N, eta, seed=J_seed)
        A = sim_activity_gpu(J, g, h0; s=s, dt=dt, t_bound=t_bound,
                minimum_norm=1e-90)
        if size(A, 2) < startstep
            avg_var = 0
        else
            avg_var = mean(A[:,startstep:end].^2)
        end
       network_avg_vars[i] = avg_var
    end
    return mean(network_avg_vars) - target_avg_var
end

if eta == -1.0f0
    lo_bound = 10f0
else
    lo_bound = 1 / (1 + eta)
end
hi_bound = 20f0
@time g = find_zero(f, (lo_bound, hi_bound), Roots.A42(), verbose=true)

out = Dict{String,Any}()
@pack! out = N, g, eta, s, J_seeds, h0_seed, wait, t_bound, dt, target_avg_var
save(joinpath(resultsdir, "Rootfinding_E$(ind)_N=$N,eta=$eta,s=$s.jld2"), out)


text:
function read_species_basics(yaml)
    n_species = length(yaml["phases"][1]["species"])
    n_reactions = length(yaml["reactions"])
    species_names = yaml["phases"][1]["species"]
    elements = yaml["phases"][1]["elements"]
    n_elements = length(elements)

    ele_matrix = zeros(n_elements, n_species)

    _species_names =
        [yaml["species"][i]["name"] for i = 1:length(yaml["species"])]

    for (i, species) in enumerate(species_names)
        spec = yaml["species"][findfirst(x -> x == species, _species_names)]
        for j = 1:n_elements
            if haskey(spec["composition"], elements[j])
                ele_matrix[j, i] = spec["composition"][elements[j]]
            end
        end
    end

    return n_species, n_reactions, species_names, elements, n_elements, ele_matrix
end
"""
    CreateSolution(mech)

Reaction mechanism is interepreted here. Part of the infomation are read in
from the yaml file, pary of them are from the pre-processed .npz file from
ReacTorch and Cantera
test for math enviroment
"""
function CreateSolution(mech)
    yaml = YAML.load_file(mech)

    #### Basics
    n_species, n_reactions, species_names,
    elements, n_elements, ele_matrix = read_species_basics(yaml)

    #### Thermo
    if yaml["phases"][1]["thermo"] == "ideal-gas" # switch to work with Cantera standard
        thermo = IdealGasThermo(yaml)
    else
        constructorThermo = Symbol(yaml["phases"][1]["thermo"], :Thermo)
        thermo = @eval($constructorThermo)(yaml)
    end

    #### Kinetic data

    npz = npzread("$mech.npz")
    MW = npz["molecular_weights"]
    efficiencies_coeffs_full = npz["efficiencies_coeffs"]
    product_stoich_coeffs = sparse(npz["product_stoich_coeffs"])
    reactant_stoich_coeffs = sparse(npz["reactant_stoich_coeffs"])
    reactant_orders = sparse(npz["reactant_orders"])
    is_reversible = npz["is_reversible"]
    Arrhenius_coeffs = npz["Arrhenius_coeffs"]
    if haskey(npz, "Arrhenius_A0")
        Arrhenius_A0 = npz["Arrhenius_A0"]
        Arrhenius_b0 = npz["Arrhenius_b0"]
        Arrhenius_Ea0 = npz["Arrhenius_Ea0"]
    else
        Arrhenius_A0 = []
        Arrhenius_b0 = []
        Arrhenius_Ea0 = []
    end

    if haskey(npz, "Troe_A")
        Troe_A = npz["Troe_A"]
        Troe_T1 = npz["Troe_T1"]
        Troe_T2 = npz["Troe_T2"]
        Troe_T3 = npz["Troe_T3"]
    else
        Troe_A = []
        Troe_T1 = []
        Troe_T2 = []
        Troe_T3 = []
    end
    Arrhenius_0 = hcat(Arrhenius_A0, Arrhenius_b0, Arrhenius_Ea0)
    Troe_ = hcat(Troe_A, Troe_T1, Troe_T2, Troe_T3)

    index_three_body = []
    index_falloff = []
    index_falloff_Troe = []
    j = 1
    for i = 1:n_reactions
        reaction = yaml["reactions"][i]
        if haskey(reaction, "type")
            if reaction["type"] == "three-body"
                push!(index_three_body, i)
            end
            if reaction["type"] == "falloff"
                push!(index_falloff, i)
                if haskey(reaction, "Troe")
                    push!(index_falloff_Troe, j)
                    j = j + 1
                else
                    push!(index_falloff_Troe, -1)
                end
            end
        end
    end

    i_reactant = []
    i_product = []
    for i = 1:n_reactions
        push!(i_reactant, findall(reactant_orders[:, i] .> 0.01))
        push!(i_product, findall(product_stoich_coeffs[:, i] .> 0.01))
    end

    vk = sparse(product_stoich_coeffs - reactant_stoich_coeffs)
    vk_sum = sum(vk, dims=1)[1, :]

    for i in 1:n_reactions
        if !((i in index_three_body) | (i in index_falloff))
            efficiencies_coeffs_full[:, i] .= 0.0
        end
    end
    efficiencies_coeffs = sparse(efficiencies_coeffs_full)

    reaction = Reaction(
        product_stoich_coeffs,
        reactant_stoich_coeffs,
        reactant_orders,
        is_reversible,
        Arrhenius_coeffs,
        Arrhenius_0,
        Troe_,
        index_three_body,
        index_falloff,
        index_falloff_Troe,
        efficiencies_coeffs,
        i_reactant,
        i_product,
        n_reactions,
        vk,
        vk_sum,
    )

    #### Transport data

    if haskey(npz, "species_viscosities_poly")
        species_viscosities_poly = npz["species_viscosities_poly"]
        thermal_conductivity_poly = npz["thermal_conductivity_poly"]
        binary_diff_coeffs_poly = npz["binary_diff_coeffs_poly"]
        poly_order = size(species_viscosities_poly)[1]
    else
        species_viscosities_poly = zeros(2, 2)
        thermal_conductivity_poly = zeros(2, 2)
        binary_diff_coeffs_poly = zeros(2, 2)
        poly_order = 6
    end

    trans = Transport(poly_order,
                      species_viscosities_poly,
                      thermal_conductivity_poly,
                      binary_diff_coeffs_poly)

    gas = Solution(
        n_species,
        n_reactions,
        MW,
        species_names,
        elements,
        ele_matrix,
        thermo,
        trans,
        reaction,
    )
    return gas
end

export CreateSolution


text:
module DA

export etkf

using Statistics
using LinearAlgebra

"""
Ensemble transform Kalman filter (ETKF)
"""
function etkf(;E, R_inv, inflation, H, y)
    D, m = size(E)

    x_m = mean(E, dims=2)
    X = (E .- x_m)/sqrt(m - 1)

    X = sqrt(inflation)*X

    y_m = H(x_m)
    Y = (vcat([H(E[:, i]) for i=1:m]...) .- y_m)'/sqrt(m - 1)
    Ω = inv(Symmetric(I + Y'*R_inv*Y))
    w = Ω*Y'*R_inv*(y - y_m)'

    E = x_m .+ X*(w .+ sqrt(m - 1)*sqrt(Ω))

    return E
end

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAutoscalingV2beta2MetricIdentifier <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    selector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: selector

    function IoK8sApiAutoscalingV2beta2MetricIdentifier(;name=nothing, selector=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta2MetricIdentifier, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiAutoscalingV2beta2MetricIdentifier, Symbol("selector"), selector)
        setfield!(o, Symbol("selector"), selector)
        o
    end
end # type IoK8sApiAutoscalingV2beta2MetricIdentifier

const _property_map_IoK8sApiAutoscalingV2beta2MetricIdentifier = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("selector")=>Symbol("selector"))
const _property_types_IoK8sApiAutoscalingV2beta2MetricIdentifier = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("selector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta2MetricIdentifier }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta2MetricIdentifier))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta2MetricIdentifier }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAutoscalingV2beta2MetricIdentifier[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta2MetricIdentifier }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta2MetricIdentifier[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta2MetricIdentifier)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAutoscalingV2beta2MetricIdentifier }, name::Symbol, val)
end


text:
# Aer

function get_basic_backends()
    qiskit.BasicAer.backends()
end

function get_basic_backend(b)
    qiskit.BasicAer.get_backend(b)
end


text:
@everywhere using AntennaExpts
@everywhere using PIFAProblem
using PIFA_GE
using ExprSearch
using ExprSearch.GE

function ExprSearch.GE.evaluate!(p::GEESParams, pop::GEPopulation, result::GEESResult,
    problem::PIFAAnt, default_expr)
    #first convert all derivtrees to exprs
    for ind in pop
        if isnull(ind.fitness)
            try
                ind.expr = get_expr(ind.derivtree)
            catch e
                if !isa(e, IncompleteException)
                    rethrow(e)
                end
                ind.expr = default_expr
                ind.fitness = Nullable{Float64}(realmax(Float64))
            end
        end
    end
    asyncdriver(pop, result)
end

function asyncdriver(pop::GEPopulation, result::GEESResult)
    np = nprocs()  # determine the number of processes available
    n = length(pop)
    i = 1
    function nextidx()
        while i <= n
            idx = i
            i += 1
            if isnull(pop[idx].fitness)
                return idx
            end
        end
        n+1
    end
    @sync begin
        for p=1:np
            if p != myid() || np == 1
                @async begin
                    #println("myid=", myid())
                    while true
                        idx = nextidx()
                        if idx > n
                            break
                        end
                        #println("idx=", idx, ", p=", p)
                        fitness = remotecall_fetch(worker_expr_fitness, p, pop[idx].expr)
                        pop[idx].fitness = Nullable{Float64}(fitness)
                        result.totalevals += 1
                        if fitness < result.fitness
                            result.fitness = fitness
                            copy!(result.tree, pop[idx].derivtree)
                            result.best_at_eval = result.totalevals
                            result.expr = pop[idx].expr
                        end
                    end
                end
            end
        end
    end
end

@everywhere problem = YagiAnt()
@everywhere worker_expr_fitness(expr) = begin
    #println("myid=", myid())
    expr_fitness(problem, expr)
end

pifa_ge(; seed=1, pop_size=20000, tournament_size=4000, iterations=20, top_keep=0.0)


text:
"""
    ProjectFile(; version=v"0.1.0")

Creates a `Project.toml`.

## Keyword Arguments
- `version::VersionNumber`: The initial version of created packages.
"""
@plugin struct ProjectFile <: Plugin
    version::VersionNumber = v"0.1.0"
end

# Other plugins like Tests will modify this file.
priority(::ProjectFile, ::typeof(hook)) = typemax(Int) - 5

function hook(p::ProjectFile, t::Template, pkg_dir::AbstractString)
    toml = Dict(
        "name" => basename(pkg_dir),
        "uuid" => string(@mock uuid4()),
        "authors" => t.authors,
        "version" => string(p.version),
        "compat" => Dict("julia" => compat_version(t.julia)),
    )
    write_project(joinpath(pkg_dir, "Project.toml"), toml)
end

# Taken from:
# https://github.com/JuliaLang/Pkg.jl/blob/v1.7.0/src/project.jl#L175-L177

function project_key_order(key::String)
    _project_key_order = ["name", "uuid", "keywords", "license", "desc", "deps", "compat"]
    return something(findfirst(x -> x == key, _project_key_order), length(_project_key_order) + 1)
end

write_project(path::AbstractString, dict) =
    open(io -> write_project(io, dict), path; write = true)
write_project(io::IO, dict) =
    TOML.print(io, dict, sorted = true, by = key -> (project_key_order(key), key))

"""
    compat_version(v::VersionNumber) -> String

Format a `VersionNumber` to exclude trailing zero components.
"""
function compat_version(v::VersionNumber)
    return if v.patch == 0 && v.minor == 0
        "$(v.major)"
    elseif v.patch == 0
        "$(v.major).$(v.minor)"
    else
        "$(v.major).$(v.minor).$(v.patch)"
    end
end


text:
@testset "land_sea_mask" begin
    lon,lat,data = GeoDatasets.landseamask(;resolution='c',grid=5)
    @test size(data,1) == length(lon)
    @test size(data,2) == length(lat)

    @test_throws ErrorException GeoDatasets.landseamask(;resolution='c',grid=-999)
    @test_throws ErrorException GeoDatasets.landseamask(;resolution='g',grid=5)

    mask = GeoDatasets.LandSeaMask()
    @test GeoDatasets.is_land(mask, 52.200475, 0.1138)
    @test !GeoDatasets.is_land(mask, 56.615841, 3.39206)
end

max_stars_repo_name: JuliaGeo/GeoDatasets

max_issues_repo_name: JuliaGeo/GeoDatasets

max_forks_repo_name: JuliaGeo/GeoDatasets


text:
export ## Types
       CompositeQSystem,
       Control,
       Dissipation,
       Field,
       Interaction,
       ParametricInteraction,
       QSystem,
       ## Methods
       dim,
       label,
       lowering,
       number,
       raising,
       strength,
       X

abstract QSystem

label(q::QSystem) = q.label
dim(q::QSystem) = q.dim
raising(q::QSystem) = diagm(sqrt(1:(dim(q)-1)), -1)
lowering(q::QSystem) = diagm(sqrt(1:(dim(q)-1)), 1)
number(q::QSystem) = raising(q) * lowering(q)
X(q::QSystem) = raising(q) + lowering(q)
hamiltonian(q::QSystem, t) = hamiltonian(q)

abstract Control
label(c::Control) = c.label

abstract Interaction
strength(i::Interaction) = i.strength
hamiltonian(i::Interaction, t) = hamiltonian(i)

abstract ParametricInteraction <: Interaction

abstract Dissipation

# [todo] - Should COmpositeQSystem <: QSystem ?
type CompositeQSystem
    # [feature] - Use something like OrderedDict for component enumeration
    subSystems::Vector{QSystem}
    interactions::Vector{Interaction}
    parametericInteractions::Vector{ParametricInteraction}
    subSystemExpansions::Vector{Vector{Vector{Int}}}
    interactionExpansions::Vector{Vector{Vector{Int}}}
    dissipators::Vector{Dissipation}
end


text:
using Documenter, DataStreams

makedocs(
    modules = [DataStreams],
)

deploydocs(
    deps = Deps.pip("mkdocs", "mkdocs-material", "python-markdown-math"),
    repo = "github.com/JuliaData/DataStreams.jl.git"
)


text:
using Base.Iterators
using LinearAlgebra

function readInput(filename::String)::Matrix{Int}
    open(filename) do file
        return transpose(hcat([parse.(Int, split(line, "")) for line in readlines(file)]...))
    end
end

function flash(mat::Matrix{Int}, flashing_at::Int)::Tuple{Matrix{Int}, Matrix{Bool}}
    w, h = size(mat)
    is_flashing = falses(size(mat))
    s_old, s_new = -1, sum(is_flashing)
    while s_old != s_new
        for i = 1:h
            for j = 1:w
                if mat[i, j] > flashing_at && !is_flashing[i, j]
                    y_min = max(1, j - 1)
                    y_max = min(w, j + 1)
                    x_min = max(1, i - 1)
                    x_max = min(h, i + 1)
                    mat[x_min:x_max, y_min:y_max] .+= 1

                    is_flashing[i, j] = true
                end
            end
        end
        s_old, s_new = s_new, sum(is_flashing)
    end

    return mat, is_flashing
end

function part1(mat::Matrix{Int}; steps::Int = 100, flashing_at::Int = 9)::Int
    n_flashes = 0
    for _ in 1:steps
        mat, _ = flash(mat .+ 1, flashing_at)
        n_flashes += sum(mat .> flashing_at)
        mat[mat .> flashing_at] .= 0
    end

    return n_flashes
end

function part2(mat::Matrix{Int}; flashing_at::Int = 9)::Int
    for i in countfrom(1)
        mat, is_flashing = flash(mat .+ 1, flashing_at)
        if all(is_flashing)
            return i
        end

        mat[mat .> flashing_at] .= 0
    end
end

mat = readInput("input.txt")

println(part1(mat, steps=100))
println(part2(mat))


text:
"""
        InvariantDistribution(dist, nonzero_inds)

    Invariant visitiation frequencies over a partitioned state space. Fields are
    `dist` (the probability distribution) and `nonzero_inds` (the indices of `dist`
    with nonzero entries).
"""
struct InvariantDistribution <: AbstractInvariantDistribution
    dist::Vector{Float64} # Distribution over the simplices
    nonzero_inds::Vector{Int} # indices of nonzero entries
end

Base.length(iv::InvariantDistribution) = length(iv.dist)
Base.size(iv::InvariantDistribution) = size(iv.dist)

get_distribution(i::InvariantDistribution) = i.dist

export
InvariantDistribution,
get_distribution


text:
# Legit.jl -- Convert French law in LEGI format to Git & Markdown
# By: Emmanuel Raviart <emmanuel.raviart@data.gouv.fr>
#
# Copyright (C) 2015 Etalab
# https://github.com/etalab/Legit.jl
#
# The Legit.jl package is licensed under the MIT "Expat" License.

const number_by_latin_extension = @compat Dict{String, String}(
  "bis" => "0000020",
  "ter" => "0000030",
  "quater" => "0000040",
  "quinquies" => "0000050",
  "sexies" => "0000060",
  "septies" => "0000070",
  "octies" => "0000080",
  "nonies" => "0000090",
  "decies" => "0000100",
  "undecies" => "0000110",
  "duodecies" => "0000120",
  "terdecies" => "0000130",
  "quaterdecies" => "0000140",
  "quindecies" => "0000150",
  "quinquedecies" => "0000150",
  "sexdecies" => "0000160",
  "septdecies" => "0000170",
  "octodecies" => "0000180",
  "novodecies" => "0000190",
  "vicies" => "0000200",
  "unvicies" => "0000210",
  "duovicies" => "0000220",
  "tervicies" => "0000230",
  "quatervicies" => "0000240",
  "quinvicies" => "0000250",
  "sexvicies" => "0000260",
  "septvicies" => "0000270",
  "octovicies" => "0000280",
  "novovicies" => "0000290",
  "tricies" => "0000300",
  "untricies" => "0000310",
  "duotricies" => "0000320",
  "tertricies" => "0000330",
  "quatertricies" => "0000340",
  "quintricies" => "0000350",
  "sextricies" => "0000360",
  "septtricies" => "0000370",
  "octotricies" => "0000380",
  "novotricies" => "0000390",
)

const number_by_slug = @compat Dict{String, String}(
  "premier" => "0000010",
  "premiere" => "0000010",
  "deuxieme" => "0000020",
  "troisieme" => "0000030",
  "quatrieme" => "0000040",
  "cinquieme" => "0000050",
  "sixieme" => "0000060",
  "septieme" => "0000070",
  "huitieme" => "0000080",
  "neuvieme" => "0000090",
  "dixieme" => "0000100",
)

abstract Node
abstract AbstractTableOfContent <: Node

type RootNode <: Node
  title::String
end

@compat type Article <: Node
  container::AbstractTableOfContent
  start_date::Union(Date, Nothing)
  stop_date::Union(Date, Nothing)
  dict::Dict  # Dict{String, Any}
  next_version::Nullable{Article}  # Next version of the same article (may have the same ID)

  function Article(container::AbstractTableOfContent, start_date::Union(Date, Nothing), stop_date::Union(Date, Nothing),
      dict::Dict, next_version::Nullable{Article})
    article = new(container, start_date, stop_date, dict, next_version)
    finalizer(article, free!)
    return article
  end
end

Article(container::AbstractTableOfContent, start_date::Union(Date, Nothing), stop_date::Union(Date, Nothing),
  dict::Dict) = @compat Article(container, start_date, stop_date, dict, Nullable{Article}())

type Document <: AbstractTableOfContent
  container::Node  # TODO: SimpleNode or RootNode or Section
  texte_version::Dict  # Dict{String, Any}
  textelr::Dict  # Dict{String, Any}

  function Document(container::Node, texte_version::Dict, textelr::Dict)
    document = new(container, texte_version, textelr)
    finalizer(document, free!)
    return document
  end
end

type NonArticle <: Node
  container::AbstractTableOfContent
  kind::String
  title::String
  content::XMLElement
end

type Changed
  articles::Array{Union(Article, NonArticle)}
  deleted_articles::Array{Article}

  Changed() = new(Article[], Article[])
end

type Section <: Node
  short_title::String
  sortable_title::String
  title::String
  child_by_name::Dict{String, Node}

  Section(short_title::String, title::String) = new(short_title, node_sortable_title(short_title), title,
    @compat Dict{String, Node}())
end

Section(title::String) = Section(title, title)

Section() = Section("", "")

type SimpleNode <: Node
  container::Node
  title::String
end

type TableOfContent <: AbstractTableOfContent
  container::AbstractTableOfContent
  start_date::Union(Date, Nothing)
  stop_date::Union(Date, Nothing)
  dict::Dict  # Dict{String, Any}

  function TableOfContent(container::AbstractTableOfContent, start_date::Union(Date, Nothing),
      stop_date::Union(Date, Nothing), dict::Dict)
    table_of_content = new(container, start_date, stop_date, dict)
    finalizer(table_of_content, free!)
    return table_of_content
  end
end

type UnparsedSection <: Node
  short_title::String
  sortable_title::String

  UnparsedSection(short_title::String) = new(short_title, node_sortable_title(short_title))
end

function commonmark(article::Article, mode::String; depth::Int = 1)
  blocks = String[
    "#" ^ depth,
    " ",
    node_title(article),
    "\n\n",
  ]

  content = commonmark(article.dict["BLOC_TEXTUEL"]["CONTENU"])
  content = join(map(strip, split(content, '\n')), '\n')
  while searchindex(content, "\n\n\n") > 0
    content = replace(content, "\n\n\n", "\n\n")
  end
  push!(blocks, strip(content))
  push!(blocks, "\n")

  nota = get(article.dict, "NOTA", nothing)
  if nota !== nothing
    push!(blocks, "\n## Nota\n\n")
    content = commonmark(nota["CONTENU"])
    content = join(map(strip, split(content, '\n')), '\n')
    while searchindex(content, "\n\n\n") > 0
      content = replace(content, "\n\n\n", "\n\n")
    end
    push!(blocks, strip(content))
    push!(blocks, "\n")
  end

  return join(blocks)
end

function commonmark(non_article::NonArticle, mode::String; depth::Int = 1)
  blocks = String[]
  if !isempty(non_article.title)
    push!(blocks,
      "#" ^ depth,
      " ",
      non_article.title,
      "\n\n",
    )
  end
  content = commonmark(non_article.content)
  content = join(map(strip, split(content, '\n')), '\n')
  while searchindex(content, "\n\n\n") > 0
    content = replace(content, "\n\n\n", "\n\n")
  end
  push!(blocks, strip(content))
  push!(blocks, "\n\n")
  return join(blocks)
end

commonmark(section::Section, mode::String; depth::Int = 1) = string(
  "#" ^ depth,
  " ",
  node_title(section),
  "\n",
  commonmark_children(section, mode; depth = depth),
)

function commonmark(xhtml_element::XMLElement; depth::Int = 1)
  blocks = String[]
  for xhtml_node in child_nodes(xhtml_element)
    if is_textnode(xhtml_node)
      push!(blocks, content(xhtml_node))
    elseif is_elementnode(xhtml_node)
      xhtml_child = XMLElement(xhtml_node)
      child_name = name(xhtml_child)
      if child_name == "blockquote"
        push!(blocks, "\n")
        child_text = commonmark(xhtml_child, depth = depth)
        push!(blocks, join(map(line -> string("> ", strip(line)), split(strip(child_text), '\n')), '\n'))
        push!(blocks, "\n")
      elseif child_name == "br"
        push!(blocks, "\n\n")
      elseif child_name in ("abbr", "acronym", "code", "div", "font", "hr", "ol", "sub", "sup", "table", "u", "ul")
        push!(blocks, string(xhtml_child))
      elseif child_name in ("em", "i")
        content_commonmark = strip(commonmark(xhtml_child, depth = depth))
        if !isempty(content_commonmark)
          push!(blocks, "_")
          push!(blocks, content_commonmark)
          push!(blocks, "_")
        end
      elseif child_name in ("h1", "h2", "h3", "h4", "h5")
        content_commonmark = strip(commonmark(xhtml_child, depth = depth))
        if !isempty(content_commonmark)
          push!(blocks, "#" ^ (parseint(child_name[2]) + 1))
          push!(blocks, content_commonmark)
          push!(blocks, "\n\n")
        end
      elseif child_name == "NOTES"
        content_commonmark = strip(commonmark(xhtml_child, depth = depth))
        if !isempty(content_commonmark)
          push!(blocks, content_commonmark)
          push!(blocks, "\n\n")
        end
      elseif child_name == "p"
        push!(blocks, "\n\n")
        content_commonmark = rstrip(commonmark(xhtml_child, depth = depth))
        if !isempty(content_commonmark)
          push!(blocks, content_commonmark)
          push!(blocks, "\n\n")
        end
      elseif child_name in ("b", "strong")
        content_commonmark = strip(commonmark(xhtml_child, depth = depth))
        if !isempty(content_commonmark)
          push!(blocks, "**")
          push!(blocks, content_commonmark)
          push!(blocks, "**")
        end
      else
        error("Unexpected XHTML element $child_name in:\n$(string(xhtml_element)).")
      end
    end
  end
  return join(blocks)
end

commonmark_children(article::Article, mode::String; depth::Int = 1, link_prefix::String = "") = ""

function commonmark_children(section::Section, mode::String; depth::Int = 1, link_prefix::String = "")
  blocks = String[
    "\n",
  ]
  children_infos = [
    (node_sortable_title(child), name, child)
    for (name, child) in section.child_by_name
  ]
  sort!(children_infos)
  if mode == "single-page"
    for (index, (sortable_number, name, child)) in enumerate(children_infos)
      if index > 1
        push!(blocks, "\n")
      end
      push!(blocks, commonmark(child, mode; depth = depth + 1))
    end
  else
    indent = "  " ^ (depth - 1)
    for (sortable_number, name, child) in children_infos
      push!(blocks, "$(indent)- [$(node_short_title(child))]($link_prefix$name)\n")
      if mode == "deep"
        push!(blocks, commonmark_children(child, mode; depth = depth + 1, link_prefix = string(link_prefix, name, '/')))
      end
    end
  end
  return join(blocks)
end

free!(article::Article) = free!(article.dict)

function free!(dict::Dict)
  for (key, value) in dict
    if key == "CONTENU"
      # LighXML.jl needs explicit garbage collecting.
      dict[key] = nothing
      free(value)
    else
      free!(value)
    end
  end
end

function free!(document::Document)
  free!(document.texte_version)
  free!(document.textelr)
end

function free!(array::Array)
  for value in array
    free!(value)
  end
end

free!(table_of_content::TableOfContent) = free!(table_of_content.dict)

function free!(value)
end

function link_articles(articles_by_id)
  for (article_id, same_id_articles) in articles_by_id
    sort!(same_id_articles, by = article -> min_date(node_start_date(article), node_stop_date(article)))
  end
  for (article_id, same_id_articles) in articles_by_id
    for (article_index, article) in enumerate(same_id_articles)
      previous_article_with_same_id = article_index == 1 ? nothing : same_id_articles[article_index - 1]
      versions = article.dict["VERSIONS"]["VERSION"]
      version_index = findfirst(version -> version["LIEN_ART"]["@id"] == article_id, versions)
      if version_index > 1
        previous_version_article = nothing
        previous_version_article_index = 0
        previous_version_index = version_index - 1
        while previous_version_index > 0
          previous_version_id = versions[previous_version_index]["LIEN_ART"]["@id"]
          previous_version_articles = get(articles_by_id, previous_version_id, Article[])
          # Note: When previous_version_articles is empty, the article has an empty date interval. Skip it.
          if !isempty(previous_version_articles)
            previous_version_article_index = findlast(previous_version_articles) do previous_version_article
              return min_date(node_start_date(previous_version_article), node_stop_date(previous_version_article)) <
                min_date(node_start_date(article), node_stop_date(article))
            end
            if previous_version_article_index > 0
              previous_version_article = previous_version_articles[previous_version_article_index]
            end
            break
          end
          previous_version_index -= 1
        end
        previous_article = previous_article_with_same_id === nothing ?
          previous_version_article :
          previous_version_article !== nothing &&
              min_date(node_start_date(previous_article_with_same_id), node_stop_date(previous_article_with_same_id)) <
              min_date(node_start_date(previous_version_article), node_stop_date(previous_version_article)) ?
            previous_version_article :
            previous_article_with_same_id
      else
        previous_article = previous_article_with_same_id
      end
      if previous_article !== nothing
        if isnull(previous_article.next_version)
          previous_article.next_version = Nullable(article)
        else
          article_start_date = node_start_date(article)
          next_article = get(previous_article.next_version)
          next_article_start_date = node_start_date(next_article)
          warn(string(
            "Previous article ",
            node_id(previous_article),
            '@',
            node_start_date(previous_article),
            '-',
            node_stop_date(previous_article),
            " of ",
            article_id,
            '@',
            article_start_date,
            '-',
            node_stop_date(article),
            " already has a next version ",
            node_id(next_article),
            '@',
            next_article_start_date,
            '-',
            node_stop_date(next_article),
          ))
          if article_start_date < next_article_start_date
            previous_article.next_version = Nullable(article)
            article.next_version = Nullable(next_article)
          elseif article_start_date == next_article_start_date
            # The same article number appears twice at the same date. Don't even try to repair.
          else
            @assert isnull(next_article.next_version)
            next_article.next_version = Nullable(article)
          end
        end
      end
    end
  end
end

function load_article(dir::String, id::String)
  article_file_path = joinpath(dir, "article", id[1:4], id[5:8], id[9:10], id[11:12], id[13:14], id[15:16], id[17:18],
    id * ".xml")
  if !ispath(article_file_path)
    warn("Missing article file $article_file_path.")
    return @compat Dict{String, Any}()
  end
  article_xml_document = parse_file(article_file_path)
  return Convertible(parse_xml_element(root(article_xml_document))) |> pipe(
    element_to_article,
    require,
  ) |> to_value
  # free(article_xml_document)
end

max_date(left::Date, right::Date) = max(left, right)

max_date(left::Date, ::Nothing) = left

max_date(::Nothing, right::Date) = right

max_date(::Nothing, ::Nothing) = nothing

min_date(left::Date, right::Date) = min(left, right)

min_date(left::Date, ::Nothing) = left

min_date(::Nothing, right::Date) = right

min_date(::Nothing, ::Nothing) = nothing

function node_dir_name(simple_node::SimpleNode)
  dir_name = slugify(node_short_title(simple_node))
  if isempty(dir_name)
    dir_name = "sans-titre"
  end
  return dir_name
end

node_dir_name(section::Section) = section.dir_name

function node_dir_name(table_of_content::Document)
  dir_name = slugify(node_short_title(table_of_content))
  if isempty(dir_name)
    dir_name = "sans-titre"
  end
  return dir_name
end

function node_dir_name(table_of_content::TableOfContent)
  dir_name = slugify(node_number_and_simple_title(node_short_title(table_of_content))[2])
  if isempty(dir_name)
    dir_name = "sans-titre"
  end
  return dir_name
end

node_filename(article::Article) = string("article-", slugify(node_number(article)), ".md")

node_filename(non_article::NonArticle) = string(non_article.kind, ".md")

node_filename(node::Node) = node_dir_name(node) * ".md"

node_git_dir(article::Article) = node_git_dir(article.container)

node_git_dir(non_article::NonArticle) = node_git_dir(non_article.container)

node_git_dir(::RootNode) = ""

function node_git_dir(node::Node)
  container_git_dir = node_git_dir(node.container)
  dir_name = node_dir_name(node)
  return isempty(container_git_dir) ? dir_name : string(container_git_dir, '/', dir_name)
end

node_git_file_path(node::Node) = string(node_git_dir(node), '/', node_filename(node))

node_id(article::Article) = article.dict["META"]["META_COMMUN"]["ID"]

node_id(document::Document) = document.texte_version["META"]["META_COMMUN"]["ID"]

node_id(table_of_content::TableOfContent) = table_of_content.dict["ID"]

node_id(node_dict::Dict) = node_dict["META"]["META_COMMUN"]["ID"]

node_name(table_of_content::AbstractTableOfContent) = node_dir_name(table_of_content)

node_name(article::Article) = node_filename(article)

node_name(non_article::NonArticle) = node_filename(non_article)

node_name(simple_node::SimpleNode) = node_dir_name(simple_node)

node_number(table_of_content::AbstractTableOfContent) = node_number(node_short_title(table_of_content))

node_number(article::Article) = get(article.dict["META"]["META_SPEC"]["META_ARTICLE"], "NUM", "")

node_number(section::Section) = node_number(node_short_title(section))

node_number(short_title::String) = node_number_and_simple_title(short_title)[1]

node_number_and_simple_title(table_of_content::AbstractTableOfContent) = node_number_and_simple_title(node_short_title(
  table_of_content))

function node_number_and_simple_title(short_title::String)
  number_fragments = String[]
  simple_title_fragments = String[]
  for fragment in split(strip(short_title))
    fragment_lower = lowercase(fragment)
    if fragment_lower == "n°"
      continue
    end
    if startswith(fragment_lower, "n°")
      fragment = fragment[endof("n°"):end]
      fragment_lower = fragment_lower[endof("n°"):end]
    end
    slug = slugify(fragment_lower)
    if slug in ("chapitre", "livre", "paragraphe", "partie", "section", "sous-paragraphe", "sous-section",
        "sous-sous-paragraphe", "titre")
      push!(simple_title_fragments, fragment)
    elseif slug in ("annexe", "legislative", "preliminaire", "reglementaire", "rubrique", "sommaire", "suite",
        "tableau")
      # Partie législative, partie réglementaire, chapître préliminaire
      push!(number_fragments, fragment)
      push!(simple_title_fragments, fragment)
    elseif isdigit(fragment) || slug in ("ier", "unique") || ismatch(r"^[ivxlcdm]+$", fragment_lower) ||
        slug in keys(number_by_latin_extension) || slug in keys(number_by_slug) ||
        length(slug) <= 3 && all(letter -> 'a' <= letter <= 'z', slug) &&
          !(slug in ("de", "des", "du", "en", "la", "le", "les")) ||
        2 <= length(slug) <= 5 && 'a' <= slug[1] <= 'z' && isdigit(slug[2 : end]) ||
        3 <= length(slug) <= 6 && all(letter -> 'a' <= letter <= 'z', slug[1 : 2]) && isdigit(slug[3 : end])
      push!(number_fragments, fragment)
      push!(simple_title_fragments, fragment)
    elseif isempty(number_fragments)
      push!(simple_title_fragments, fragment)
    else
      break
    end
  end
  @assert !isempty(simple_title_fragments) "Empty simplification of title for: $title."
  return join(number_fragments, ' '), join(simple_title_fragments, ' ')
end

function node_short_title(document::Document)
  short_title = join(
    split(document.texte_version["META"]["META_SPEC"]["META_TEXTE_VERSION"]["TITRE"]), ' ')
  nor = get(document.texte_version["META"]["META_SPEC"]["META_TEXTE_CHRONICLE"], "NOR", "")
  if !isempty(nor)
    short_title = string(short_title, " (", nor, ")")
  end
  return short_title
end

node_short_title(non_article::NonArticle) = @compat Dict{String, String}(
  "notes" => "Nota",
  "signataires" => "Signataires",
  "visas" => "Visas",
)[non_article.kind]

node_short_title(section::Section) = section.short_title

node_short_title(section::UnparsedSection) = section.short_title

node_short_title(node::Node) = node_title(node)

node_sortable_title(article::Article) = node_sortable_title(node_number(article), node_short_title(article))

node_sortable_title(non_article::NonArticle) = @compat Dict{String, String}(
  "notes" => "9999998",
  "signataires" => "9999997",
  "visas" => "0000002",
)[non_article.kind]

node_sortable_title(section::Section) = section.sortable_title

node_sortable_title(simple_node::SimpleNode) = node_short_title(simple_node)

node_sortable_title(section::UnparsedSection) = section.sortable_title

node_sortable_title(table_of_content::AbstractTableOfContent) = node_sortable_title(node_short_title(table_of_content))

node_sortable_title(short_title::String) = isempty(short_title) ?
  "" :
  node_sortable_title(node_number_and_simple_title(short_title)...)

function node_sortable_title(number::String, simple_title::String)
  if isempty(number)
    return slugify(simple_title)
  end
  number_fragments = String[]
  slug = slugify(number)
  slug = replace(slug, "-a-l-article-", "-")
  for fragment in split(slug, '-')
    if isdigit(fragment)
      @assert length(fragment) <= 6 "Fragment is too long: $fragment"
      push!(number_fragments, string("000000", fragment, '0')[end - 6 : end])
    elseif fragment == "preliminaire"
      push!(number_fragments, "0000005")
    elseif fragment in ("ier", "legislative", "unique")
      push!(number_fragments, "0000010")
    elseif fragment in ("reglementaire", "suite")
      push!(number_fragments, "0000020")
    elseif fragment == "rubrique"
      push!(number_fragments, "7000000")
    elseif fragment == "sommaire"
      push!(number_fragments, "8000000")
    elseif fragment == "annexe"
      push!(number_fragments, "9000000")
    elseif fragment == "tableau"
      push!(number_fragments, "9500000")
    elseif ismatch(r"^[ivxlcdm]+$", fragment)
      value = 0
      for letter in fragment
        digit = [
          'i' => 1,
          'v' => 5,
          'x' => 10,
          'l' => 50,
          'c' => 100,
          'd' => 500,
          'm' => 1000,
        ][letter]
        if digit > value
          value = digit - value
        else
          value += digit
        end
      end
      @assert value < 1000000
      push!(number_fragments, string("000000", value, '0')[end - 6 : end])
    else
      number = get(number_by_latin_extension, fragment, "")
      if !isempty(number)
        push!(number_fragments, number)
      else
        number = get(number_by_slug, fragment, "")
        if !isempty(number)
          push!(number_fragments, number)
        elseif length(fragment) <= 3 && all(letter -> 'a' <= letter <= 'z', fragment) &&
            !(fragment in ("de", "des", "du", "en", "la", "le", "les"))
          push!(number_fragments, fragment)
        elseif 2 <= length(fragment) <= 5 && 'a' <= fragment[1] <= 'z' && isdigit(fragment[2 : end])
          push!(number_fragments, fragment[1 : 1])
          push!(number_fragments, string("000000", fragment[2 : end], '0')[end - 6 : end])
        elseif 3 <= length(fragment) <= 6 && all(letter -> 'a' <= letter <= 'z', fragment[1 : 2]) &&
            isdigit(fragment[3 : end])
          push!(number_fragments, fragment[1 : 2])
          push!(number_fragments, string("000000", fragment[3 : end], '0')[end - 6 : end])
        else
          push!(number_fragments, fragment)
        end
      end
    end
  end
  return join(number_fragments, '-')
end

node_start_date(article::Article) = article.start_date

node_start_date(document::Document) = get(document.texte_version["META"]["META_SPEC"]["META_TEXTE_VERSION"],
  "DATE_DEBUT", nothing)

node_start_date(table_of_content::TableOfContent) = table_of_content.start_date

function node_stop_date(node::Union(Article, TableOfContent))
  stop_date = node.stop_date
  if stop_date !== nothing && stop_date < node.start_date
    # May occur when ETAT = MODIFIE_MORT_NE.
    stop_date = node.start_date
  end
  return stop_date
end

node_stop_date(document::Document) = get(document.texte_version["META"]["META_SPEC"]["META_TEXTE_VERSION"],
  "DATE_FIN", nothing)

node_structure(table_of_content::Document) = table_of_content.textelr["STRUCT"]

node_structure(table_of_content::TableOfContent) = table_of_content.dict["STRUCTURE_TA"]

node_title(article::Article) = string("Article ", node_number(article))

node_title(simple_node::SimpleNode) = simple_node.title

node_title(section::Section) = section.title

node_title(document::Document) = join(
  split(document.texte_version["META"]["META_SPEC"]["META_TEXTE_VERSION"]["TITREFULL"]), ' ')

node_title(table_of_content::TableOfContent) = table_of_content.dict["TITRE_TA"]

parse_section_commonmark(repository::GitRepo, ::Nothing, section_short_title::String) = Section(section_short_title)

function parse_section_commonmark(repository::GitRepo, entry::GitTreeEntry, section_short_title::String)
  oid = Oid(entry)
  blob = lookup_blob(repository, oid)
  lines = split(text(blob), '\n')
  if !isempty(lines)
    title_line = lines[1]
    @assert startswith(title_line, "# ")
    section_title = title_line[3 : end]
    section = Section(isempty(section_short_title) ? section_title : section_short_title, section_title)
    if length(lines) >= 2
      @assert isempty(lines[2])
      for list_item_line in lines[3 : end]
        list_item_line_match = match(r"^- \[(?P<short_title>.+)\]\((?P<slug>.+)\)$", list_item_line)
        if list_item_line_match !== nothing
          section.child_by_name[list_item_line_match.captures[2]] = UnparsedSection(list_item_line_match.captures[1])
        end
      end
    end
    return section
  end
  return Section(section_short_title)
end

function parse_structure(table_of_content::AbstractTableOfContent, articles_by_id::Dict{String, Vector{Article}},
    changed_by_message_by_date::Dict{Date, Dict{String, Changed}}, dir::String)
  structure = node_structure(table_of_content)
  table_of_content_start_date = node_start_date(table_of_content)
  if table_of_content_start_date === nothing
    return
  end
  table_of_content_stop_date = node_stop_date(table_of_content)

  children = TableOfContent[]
  for lien_section_ta in get(structure, "LIEN_SECTION_TA", Dict{String, Any}[])
    lien_start_date = get(lien_section_ta, "@debut", nothing)
    if lien_start_date === nothing
      continue
    end
    if table_of_content_stop_date !== nothing && lien_start_date >= table_of_content_stop_date
      continue
    end
    lien_stop_date = get(lien_section_ta, "@fin", nothing)
    if lien_stop_date !== nothing && (lien_stop_date <= table_of_content_start_date ||
        lien_stop_date <= lien_start_date)
      continue
    end

    section_ta = nothing
    section_ta_file_path = joinpath(dir, "section_ta" * lien_section_ta["@url"])
    if !ispath(section_ta_file_path)
      warn("Missing SECTION_TA file $section_ta_file_path.")
      continue
    end
    try
      section_ta_xml_document = parse_file(section_ta_file_path)
      section_ta = Convertible(parse_xml_element(root(section_ta_xml_document))) |> pipe(
        element_to_section_ta,
        require,
      ) |> to_value
      # free(section_ta_xml_document)
    catch
      warn("An exception occured in file $section_ta_file_path.")
      rethrow()
    end
    child = TableOfContent(table_of_content, max_date(table_of_content_start_date, lien_start_date),
      min_date(table_of_content_stop_date, lien_stop_date), section_ta)
    # When the start date of the table of content is in conflict with the stop date of its previous version, consider
    # that the start date is right and correct the stop date.
    if !isempty(children)
      previous_child = children[end]
      if node_dir_name(previous_child) == node_dir_name(child) &&
          (previous_child.stop_date === nothing || previous_child.stop_date > child.start_date)
        previous_child.stop_date = child.start_date
      end
    end
    push!(children, child)
  end
  for child in children
    if child.stop_date !== nothing && child.start_date >= child.stop_date
      continue
    end
    parse_structure(child, articles_by_id, changed_by_message_by_date, dir)
  end

  child_articles = Article[]
  for lien_article in get(structure, "LIEN_ART", Dict{String, Any}[])
    lien_start_date = get(lien_article, "@debut", nothing)
    if lien_start_date === nothing
      continue
    end
    if table_of_content_stop_date !== nothing && lien_start_date >= table_of_content_stop_date
      continue
    end
    lien_stop_date = get(lien_article, "@fin", nothing)
    if lien_stop_date !== nothing && (lien_stop_date <= table_of_content_start_date ||
        lien_stop_date <= lien_start_date)
      continue
    end

    article_id = lien_article["@id"]
    article_dict = load_article(dir, article_id)
    if isempty(article_dict)
      continue
    end
    article = Article(table_of_content, max_date(table_of_content_start_date, lien_start_date),
      min_date(table_of_content_stop_date, lien_stop_date), article_dict)
    # When the start date of the article is in conflict with the stop date of its previous version, consider
    # that the start date is right and correct the stop date.
    if !isempty(child_articles)
      previous_article = child_articles[end]
      if node_filename(previous_article) == node_filename(article) &&
          (previous_article.stop_date === nothing || previous_article.stop_date > article.start_date)
        previous_article.stop_date = article.start_date
      end
    end
    push!(child_articles, article)
  end
  for article in child_articles
    if article.stop_date !== nothing && article.start_date >= article.stop_date
      continue
    end
    article_dict = article.dict
    article_id = article_dict["META"]["META_COMMUN"]["ID"]
    meta_article = article_dict["META"]["META_SPEC"]["META_ARTICLE"]
    same_id_articles = get!(articles_by_id, article_id) do
      return Article[]
    end
    push!(same_id_articles, article)
    try
      start_messages = String[]
      stop_messages = String[]
      for lien in get(article_dict["LIENS"], "LIEN", Dict{String, Any}[])
        if get(lien, "@datesignatexte", nothing) === nothing
          continue
        end

        if lien["@sens"] == "cible" && lien["@typelien"] in ("CREE", "DEPLACE", "MODIFIE") ||
            lien["@sens"] == "source" && lien["@typelien"] in ("MODIFICATION", "TRANSPOSITION")
          if meta_article["DATE_DEBUT"] <= lien["@datesignatexte"]
            info("Unexpected date $(lien["@datesignatexte"]) after DATE_DEBUT article $(meta_article["DATE_DEBUT"]) " *
              "in $article_id for: $lien. Ignoring link...")
          else
            message = split(lien["^text"], " - ")[1]
            if !(message in start_messages)
              push!(start_messages, message)
            end
          end
        elseif lien["@sens"] == "source" && lien["@typelien"] in ("ABROGATION", "DISJONCTION", "PEREMPTION",
            "SUBSTITUTION", "TRANSFERT") ||
            lien["@sens"] == "cible" && lien["@typelien"] in ("ABROGE", "DISJOINT", "PERIME", "TRANSFERE")
          stop_date = get(meta_article, "DATE_FIN", nothing)
          if stop_date !== nothing
            if stop_date <= lien["@datesignatexte"]
              info("Unexpected date $(lien["@datesignatexte"]) after DATE_FIN article $(stop_date) in " *
                "$article_id for: $lien. Ignoring link...")
            else
              message = split(lien["^text"], " - ")[1]
              if !(message in stop_messages)
                push!(stop_messages, message)
              end
            end
          end
        end
      end

      creation_date = node_start_date(article)
      @assert creation_date !== nothing
      if isempty(start_messages)
        push!(start_messages, "Modifications d'origine indéterminée")
      end
      changed_by_message = get!(changed_by_message_by_date, creation_date) do
        return @compat Dict{String, Changed}()
      end
      changed = get!(changed_by_message, join(start_messages, ", ", " et ")) do
        return Changed()
      end
      push!(changed.articles, article)

      deletion_date = node_stop_date(article)
      if deletion_date !== nothing
        if isempty(stop_messages)
          push!(stop_messages, "Suppressions d'origine indéterminée")
        end
        changed_by_message = get!(changed_by_message_by_date, deletion_date) do
          return @compat Dict{String, Changed}()
        end
        changed = get!(changed_by_message, join(stop_messages, ", ", " et ")) do
          return Changed()
        end
        push!(changed.deleted_articles, article)
      end
    catch
      warn("An exception occured in $(node_filename(article)) [$(get(meta_article, "ETAT", "inconnu"))]: $article_id.")
      rethrow()
    end
  end
end

function parse_xml_element(xml_element::XMLElement)
  element = @compat Dict{String, Any}()
  for attribute in attributes(xml_element)
    element[string('@', name(attribute))] = value(attribute)
  end
  previous = element
  for xml_node in child_nodes(xml_element)
    if is_textnode(xml_node)
      if previous === element
        element["^text"] = get(element, "^text", "") * content(xml_node)
      else
        previous["^tail"] = get(element, "^tail", "") * content(xml_node)
      end
    elseif is_elementnode(xml_node)
      xml_child = XMLElement(xml_node)
      child_name = name(xml_child)
      if child_name == "CONTENU"
        @assert !(child_name in element)
        element[child_name] = xml_child
      else
        child = parse_xml_element(xml_child)
        same_children = get!(element, child_name) do
          return Dict{String, Any}[]
        end
        push!(same_children, child)
        previous = child
      end
    end
  end
  return element
end

function slugify(string::String; separator::Char = '-', transform::Function = lowercase)
  simplified = replace(replace(string, "N°", "no "), "n°", "no ")
  return Slugify.slugify(simplified, separator = separator, transform = transform)
end


text:
function translate(x)
  A0 = [0 0 0 0
        x[1] 0 0 0
        x[2] x[3] 0 0
        x[4] x[5] x[6] 0]

  A1 = [0 0 0 0
        x[7] 0 0 0
        x[8] x[9] 0 0
        x[10] x[11] x[12] 0]

  B0 = [0 0 0 0
        x[13] 0 0 0
        x[14] x[15] 0 0
        x[16] x[17] x[18] 0]

  B1 = [0 0 0 0
        x[19] 0 0 0
        x[20] x[21] 0 0
        x[22] x[23] x[24] 0]

  α = [x[25]; x[26]; x[27]; x[28]]

  β1 = [x[29]; x[30]; x[31]; x[32]]
  β2 = [x[33]; x[34]; x[35]; x[36]]
  β3 = [x[37]; x[38]; x[39]; x[40]]
  β4 = [x[41]; x[42]; x[43]; x[44]]
  return A0,A1,B0,B1,α,β1,β2,β3,β4
end


text:
# ---
# title: "`mgradient`"
# id: op_mgradient
# cover: assets/mgradient.png
# ---

# There are three commonly used morphological gradient definitions: the beucher gradident
# `dilate(img) - erode(img)`, the external half-gradient `dilate(img) - img`, and the
# internal half-gradident `img - erode(img)`.

using ImageMorphology
using TestImages
using ImageBase
using ImageShow

img = restrict(testimage_dip3e("fig0940")) # rice
img01 = @. Gray(img > 0.5)
mosaic(img, img01; nrow=1)

# ## Beucher gradient

# The default mode is beucher gradient defined as `dilate(img) - erode(img)`.

mosaic(mgradient(img), mgradient(img01); nrow=1)

# There are two nice properties of the beucher gradient:
#
# - it is always positive if the structuring element is symmetric
# - it is a self-complementary operator -- `mgradient(img) == mgradient(complement.(img))`

all(mgradient(rand(32, 32)) .> 0) # always positive

#-

A = rand(32, 32)
mgradient(A) == mgradient(1 .- A) # self-complementary

# ## Half gradient

# If we look closely, we can see that the beucher gradient boundary is always two-pixel wide
# -- one pixel per each side.

A = falses(7, 7)
A[3:5, 3:5] .= true
A

#-

Int.(mgradient(A))

# This is because `dilate` extends the boundary for one pixel, while `erode` shrinks for one
# pixel. This observation introduces so-called _half-gradients_.
#
# The _external half-gradient_, also known as _half-gradient by dilation_, is defined as
# `dilate(img) - img`:

Int.(mgradient(A; mode=:external)) # external boundary

# The _internal half-gradient_, also known as _half-gradient by erosion_, is defined as
# `img - erode(img)`:

Int.(mgradient(A; mode=:internal)) # internal boundary

# Note that the result of internal gradient is inside the original image, and that of
# external gradient is outside the original image. Also, external gradient and internal
# are complementary to each other:

A = rand(32, 32)
mgradient(A; mode=:external) == mgradient(1 .- A; mode=:internal)

# ## Thick gradient

# The gradient thickness is controlled by the SE size, or the `r` parameter. When `r>1`, it
# is usually called _thick gradient_.

mosaic([mgradient(img01; r) for r in 1:3]; nrow=1)

# ## See also

# There are many useful gradient definitions in image processing. For instance, ImageBase.jl
# provides the finite-difference version of gradients `fdiff` and `fgradient`.
# ImageFiltering.jl provides `sobel`, `prewitt` and other gradient filters.

fout = abs.(ImageBase.FiniteDiff.fdiff(img; dims=1))
mout = mgradient(img; dims=1)
mosaic(fout, mout; nrow=1)

# For a comprehensive and more accurate documentation, please check the [`mgradient`](@ref)
# reference page.

# save cover image #src
using FileIO #src
mkpath("assets") #src
img = Gray.(restrict(testimage("blobs"))) #src
save("assets/mgradient.png", mgradient(img; r=2)) #src


text:
"""
    calculate_stored_energy(sim::Simulation; consider_multiplicity::Bool = true)

Calculates and returns the energy stored in the [`ElectricField`](@ref) of a
[`SolidStateDetector`](@ref) in a given [`Simulation`](@ref) in units of J.

## Arguments
* `sim::Simulation{T}`: [`Simulation`](@ref) with `sim.detector` for which the stored energy is calculated.

## Keywords
* `consider_multiplicity::Bool = true`: Whether symmetries of the system should be taken into account.
    For example, in case of true coaxial detector center around the origin and calculated on a cartesian grid
    with the `x-axis` going from `[0, x_max]` and the `y-axis` going from `[0, y_max]` the multiplicity is 4
    and, if `consider_multiplicity == true`, the returned value is already multiplied by 4.
"""
function calculate_stored_energy(sim::Simulation; consider_multiplicity::Bool = true)
    @assert !ismissing(sim.electric_potential) "Electric potential has not been calculated yet. Please run `calculate_electric_potential!(sim)` first."
    calculate_stored_energy(sim.electric_potential, sim.ϵ_r; consider_multiplicity)
end

function calculate_stored_energy(ep::ElectricPotential{T,3,CS}, ϵ::DielectricDistribution{T,3,CS}; consider_multiplicity::Bool = true) where {T <: SSDFloat, CS}
    cylindrical = CS == Cylindrical
    phi_2D = cylindrical && size(ep, 2) == 1
    ep3d = phi_2D ? get_2π_potential(ep, n_points_in_φ = 2) : _get_closed_potential(ep)
    grid = ep3d.grid
    W::T = 0
    for i3 in 1:size(grid, 3)-1
        for i2 in 1:size(grid, 2)-1
            for i1 in 1:size(grid, 1)-1
                w1, w2, w3 = voxel_widths(grid, i1, i2, i3)
                dV = voxel_volume(grid, i1, i2, i3, w1, w2, w3)

                _ϵ = ϵ.data[i1 + 1, i2 + 1, i3 + 1]

                ep000 = ep3d.data[i1    , i2    , i3    ]
                ep100 = ep3d.data[i1 + 1, i2    , i3    ]
                ep010 = ep3d.data[i1    , i2 + 1, i3    ]
                ep110 = ep3d.data[i1 + 1, i2 + 1, i3    ]
                ep001 = ep3d.data[i1    , i2    , i3 + 1]
                ep101 = ep3d.data[i1 + 1, i2    , i3 + 1]
                ep011 = ep3d.data[i1    , i2 + 1, i3 + 1]
                ep111 = ep3d.data[i1 + 1, i2 + 1, i3 + 1]

                efv1 = ( (ep100 - ep000) + (ep110 - ep010) + (ep101 - ep001) + (ep111 - ep011) ) / (4 * w1)
                efv2 = if cylindrical
                    _w2 = (grid[2].ticks[i2 + 1] - grid[2].ticks[i2])
                    if i1 == 1
                        ((ep110 - ep100)/(_w2*grid[1].ticks[i1+1]) +
                        (ep111 - ep101)/(_w2*grid[1].ticks[i1+1])) / 2
                    else
                        ((ep010 - ep000)/(_w2*grid[1].ticks[i1]) +
                         (ep110 - ep100)/(_w2*grid[1].ticks[i1+1]) +
                         (ep011 - ep001)/(_w2*grid[1].ticks[i1]) +
                         (ep111 - ep101)/(_w2*grid[1].ticks[i1+1])) / 4
                    end
                else
                    ( (ep010 - ep000) + (ep110 - ep100) + (ep011 - ep001) + (ep111 - ep101) ) / (4 * w2)
                end
                efv3 = ( (ep001 - ep000) + (ep101 - ep100) + (ep011 - ep010) + (ep111 - ep110) ) / (4 * w3)
                W += sum((efv1, efv2, efv3).^2) * dV * _ϵ
            end
        end
    end
    E = W * ϵ0 / 2 * u"J"
    phi_2D && (E *= 2)
    return consider_multiplicity ? E * multiplicity(grid) : E
end


text:
module Contribution

using ..TypeDBClient

const g = TypeDBClient

export entity_set_owns

function entity_set_owns(thing_type::Type{<:g.AbstractThingType},
                entity::AbstractString,
                attribute_type::AbstractString,
                transaction::g.AbstractCoreTransaction,
                is_key = false,
                overriden::g.Optional{<:AbstractString} = nothing)

    cm = ConceptManager(transaction)
    loc_entity = nothing
    rem_entitiy = nothing
    loc_attribute = nothing
     try
         loc_entity = get(cm, thing_type, entity)
         rem_entitiy = g.as_remote(loc_entity, transaction)
         loc_attribute = get(cm, AttributeType, attribute_type)

         if overriden !== nothing
             overriden_attr = get(cm, AttributeType, overriden)
         else
             overriden_attr = nothing
         end

         @assert loc_entity !== nothing &&
                 rem_entitiy !== nothing &&
                 loc_attribute !== nothing

         set_owns(rem_entitiy, loc_attribute, is_key, overriden_attr)
     catch ex
         if loc_entity !== nothing
             @debug loc_entity
         else
             @debug "get entity_type not successfull. Please check your input"
         end
         if loc_attribute !== nothing
             @debug loc_attribute
         else
             @debug "get attribute not successfull. Please check your input"
         end
         rethrow(ex)
     end
     return nothing
 end

 function entity_set_owns(entity::AbstractString,
            attribute_type::AbstractString,
            transaction,
            is_key = false,
            overriden::g.Optional{String} = nothing)

     entity_set_owns(EntityType, entity, attribute_type, transaction, is_key, overriden)
 end

end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export divider

"""
    divider(;kwargs...)

A Divider component.
Horizontal line with optional label or vertical divider. For more information, see: https://mantine.dev/core/divider/
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `className` (String; optional): Often used with CSS to style elements with common properties
- `color` (a value equal to: "dark", "gray", "red", "pink", "grape", "violet", "indigo", "blue", "cyan", "teal", "green", "lime", "yellow", "orange"; optional): Line color from theme, defaults to gray in light color scheme and to dark in dark color scheme
- `label` (String; optional): Adds text after line in horizontal orientation
- `labelPosition` (a value equal to: "right", "left", "center"; optional): Label position
- `orientation` (a value equal to: "horizontal", "vertical"; optional): Line orientation
- `size` (a value equal to: "xs", "sm", "md", "lg", "xl" | Real; optional): Sets height in horizontal orientation and width in vertical
- `style` (Dict; optional): Inline style override
- `variant` (a value equal to: "dashed", "dotted", "solid"; optional): Divider borderStyle
"""
function divider(; kwargs...)
        available_props = Symbol[:id, :className, :color, :label, :labelPosition, :orientation, :size, :style, :variant]
        wild_props = Symbol[]
        return Component("divider", "Divider", "dash_mantine_components", available_props, wild_props; kwargs...)
end


text:
using EasyDataAugmentation
using Test

@testset "EasyDataAugmentation.jl" begin
    sentence = "Hello world!"
    word = "world"
    @test EasyDataAugmentation.substitute_word(sentence, word) == "Hello !"

    input_file = "../data/sst2_train_500.txt"
    output_file = "../data/output.txt"
    @test EasyDataAugmentation.easy_data_augmentation(input_file, output_file) == output_file
    @test isfile(output_file)

    sentence = "Try' cleanInNG This` SENTENCE"
    clean_sentence = "try cleaninng this sentence"
    @test EasyDataAugmentation.clean_text_sentence(sentence) == clean_sentence

    sentence = "Try' cleanInNG This`            SENTENCE"
    clean_sentence = "try cleaninng this sentence"
    @test EasyDataAugmentation.clean_text_sentence(sentence) == clean_sentence

    sentence = "Try' cle@nInNG!        This` SENTENCE \n That SENTENCE"
    clean_sentence = "try cleninng this sentence that sentence"
    @test EasyDataAugmentation.clean_text_sentence(sentence) == clean_sentence

    sentence = " clean THIS SENTENCE @·"
    clean_sentence = "clean this sentence"
    @test EasyDataAugmentation.clean_text_sentence(sentence) == clean_sentence

end


text:
using DataFrames

version = 3
othp = 0.00
telp = 0.15
stp = 0.00
infile = "acs.txt"
outfile = "acs_dirty$version.txt"
f = open(outfile, "w")

df = readtable(infile, separator=',', header=false)

iST = 1
iTEL = 14
d = ncol(df) # 37
n = nrow(df) # 671153

# fill in p
p = othp * ones(d)
p[iTEL] = telp
p[iST] = stp

for i in 1:n
    row = Array{Any}(df[i,:]) # sketchy to allow us to print ""
    for j in 1:d
        if rand() < p[j]
            row[j] = ""
        end
    end

    println(f, join(row, ","))
end

close(f)

max_stars_repo_name: jfeser/ImputeDB

max_issues_repo_name: jfeser/ImputeDB

max_forks_repo_name: jfeser/ImputeDB


text:
using DataFrames

macro c(ret_type, func, arg_types, lib)
    local args_in = Any[ symbol(string('a',x)) for x in 1:length(arg_types.args) ]
    quote
        $(esc(func))($(args_in...)) = ccall( ($(string(func)), $(Expr(:quote, lib)) ),
                                            $ret_type, $arg_types, $(args_in...) )
    end
end

#stdlib stuff
@c Clong strtol (Ptr{UInt8}, Ptr{Ptr{UInt8}}, Cint) libc
@c Clonglong strtoll (Ptr{UInt8}, Ptr{Ptr{UInt8}}, Cint) libc
@c Cfloat strtof (Ptr{UInt8}, Ptr{Ptr{UInt8}}) libc
@c Cdouble strtod (Ptr{UInt8}, Ptr{Ptr{UInt8}}) libc

abstract AbstractPostgresType{T}

type PostgresType{T} <: AbstractPostgresType{T}
    name::Symbol
    naval::T
end

type PostgresEnumType{T <: UTF8String} <: AbstractPostgresType{T}
    name::Symbol
    naval::T
    enumvals::Set{T}
end

type PostgresDomainType{T} <: AbstractPostgresType{T}
    name::Symbol
    basetype::PostgresType{T}
end

function Base.show(io::IO, t::PostgresType)
        print(io, "$(t.name) -> $(typeof(t.naval))")
end

function Base.show(io::IO, t::PostgresEnumType)
        print(io, "$(t.name) ∈ $(t.enumvals)")
end

function Base.show(io::IO, t::PostgresDomainType)
        print(io, "($(t.name) <: $(t.basetype.name)) -> $(typeof(t.basetype.naval))")
end

naval(t::PostgresType) = t.naval
naval(t::PostgresEnumType) = t.naval
naval(t::PostgresDomainType) = t.basetype.naval

###########################################################################
### base parsers

function unsafe_parse(t::AbstractPostgresType, value::UTF8String)
    value
end

function unsafe_parse(t::AbstractPostgresType, ptr::Ptr{UInt8})
    unsafe_parse(t, utf8(ptr))
end

function unsafe_parse(t::PostgresDomainType, value::UTF8String)
    unsafe_parse(t.basetype, value)
end

###########################################################################
### c parsers

function unsafe_parse{T <: Int16}(::PostgresType{T}, ptr::Ptr{UInt8})
    convert(Int16, strtol(ptr, C_NULL, 10))
end

function unsafe_parse{T <: Int32}(::PostgresType{T}, ptr::Ptr{UInt8})
    strtol(ptr, C_NULL, 10)
end

function unsafe_parse{T <: Int64}(::PostgresType{T}, ptr::Ptr{UInt8})
    strtoll(ptr, C_NULL, 10)
end

function unsafe_parse{T <: Float32}(::PostgresType{T}, ptr::Ptr{UInt8})
    strtof(ptr, C_NULL)
end

function unsafe_parse{T <: Float64}(::PostgresType{T}, ptr::Ptr{UInt8})
    strtod(ptr, C_NULL)
end

###########################################################################
### julia parsers

#should not need with the c parsers
#function unsafe_parse{T <: Real}(::PostgresType{T}, value::UTF8String)
#    parse(T, value)
#end

function unsafe_parse{T <: Char}(t::PostgresType{T}, value::UTF8String)
    value == "" ? Char('\0') : Char(value[1])
end

function unsafe_parse{T <: BigFloat}(::PostgresType{T}, value::UTF8String)
    parse(BigFloat, value)
end

function unsafe_parse{T <: Date}(::PostgresType{T}, value::UTF8String)
    y, m, d = split(value, "-")
    Date(parse(Int, y), parse(Int, m), parse(Int, d))
end

function unsafe_parse{T <: Bool}(::PostgresType{T}, value::UTF8String)
    value == "t"
end

function unsafe_parse{T <: Vector{UInt8}}(::PostgresType{T}, value::UTF8String)
    hex2bytes(value[3:end])
end

function unsafe_parse{T <: BitVector}(::PostgresType{T}, value::UTF8String)
    BitVector([c=='1' for c in value])
end

#function unsafe_parse{T <: DateTime}(::PostgresType{T}, value::UTF8String)
#    DateTime(value[1:end-3], "y-m-d H:M:S.s")
#end

###########################################################################
### value converters

type PostgresValue{T}
    pgtype::AbstractPostgresType
    value::UTF8String
end
function Base.show(io::IO, val::PostgresValue)
    print(io, val.value)
end
function Base.showall(io::IO, val::PostgresValue)
    esc = val.pgtype.name == :bytea ? "e" : ""
    print(io, "$esc'$(val.value)'::$(val.pgtype.name)")
end

PostgresValue(val) =
    PostgresValue(base_types_byname[:text]utf8, utf8(string(val)))

PostgresValue{T <: AbstractString}(val::T) =
    PostgresValue{T}(base_types_byname[:text], utf8(string(val)))

PostgresValue{T <: Int16}(val::T) =
    PostgresValue{T}(base_types_byname[:int2], utf8(string(val)))

PostgresValue{T <: Int32}(val::T) =
    PostgresValue{T}(base_types_byname[:int4], utf8(string(val)))

PostgresValue{T <: Int64}(val::T) =
    PostgresValue{T}(base_types_byname[:int8], utf8(string(val)))

PostgresValue{T <: Float32}(val::T) =
    PostgresValue{T}(base_types_byname[:float4], utf8(string(val)))

PostgresValue{T <: Float64}(val::T) =
    PostgresValue{T}(base_types_byname[:float8], utf8(string(val)))

PostgresValue{T <: BigFloat}(val::T) =
    PostgresValue{T}(base_types_byname[:numeric], utf8(string(val)))

PostgresValue{T <: Date}(val::T) =
    PostgresValue{T}(base_types_byname[:date], utf8(string(val)))

PostgresValue{T <: Bool}(val::T) =
    PostgresValue{T}(base_types_byname[:bool], utf8(string(val)))

function PostgresValue{T <: Vector{UInt8}}(val::T)
    v = "\\\\x" * join(map(x -> repr(x)[3:end], val))
    PostgresValue{T}(base_types_byname[:bytea], v)
end

function PostgresValue{T <: BitVector}(val::T)
    v = join(map(x -> x ? "1" : "0", val))
    PostgresValue{T}(base_types_byname[:varbit], v)
end

function pgtypes(df::DataFrame)
    a = AbstractPostgresType[]
    for i in 1:size(df)[2]
        if issubtype(typeof(df[i]), PooledDataArray)
            val = df[i].pool[1]
        else
            val = df[i].data[1]
        end
        push!(a, PostgresValue(val).pgtype)
    end
    a
end

###########################################################################
### instances

# Oids have to be unique against the _whole_ database,
# while the name of a type only has to be unique per namespace.

# This is everything the manual documents.
# There is more but they are obsolete or are not used for data columns.

#XXX Do not use NaN for naval because (NaN == NaN) == false.
#    this will fail the unit tests for round trips
#    (there is big rant on NaN from the developers.)

base_types = Dict(

    # default if we cannot find it in the Dict
    0            =>  PostgresType{UTF8String}(:jlunknown, UTF8String("∅")),

    16           =>  PostgresType{Bool}(:bool, false),
    17           =>  PostgresType{Vector{UInt8}}(:bytea, UInt8[0]),

    #### BITS
    1560         =>  PostgresType{BitVector}(:bit, BitVector([false])),
    1562         =>  PostgresType{BitVector}(:varbit, BitVector([false])),

    #### NUMBERS
    20           =>  PostgresType{Int64}(:int8, 0),
    21           =>  PostgresType{Int16}(:int2, 0),
    23           =>  PostgresType{Int32}(:int4, 0),
    700          =>  PostgresType{Float32}(:float4, 0),
    701          =>  PostgresType{Float64}(:float8, 0),
    1700         =>  PostgresType{BigFloat}(:numeric, 0),
    # good luck :) depends on locale.
    # Currency package looks pretty sweet for this.
    #790          =>  PostgresType{Float64}(:money, 0),

    #### STRINGS
    # varchar and bpchar (blank padded char) are realy subsets of text.
    # PG does not really have a character type like julia.
    25           =>  PostgresType{UTF8String}(:text, UTF8String("∅")),
    1043         =>  PostgresType{UTF8String}(:varchar, UTF8String("∅")),
    1042         =>  PostgresType{UTF8String}(:bpchar, UTF8String("∅")),
    705          =>  PostgresType{UTF8String}(:unknown, UTF8String("∅")),
    # probably deprecated. Cast it and it always comes out bpchar.
    #18           =>  PostgresType{UTF8String}(:char, UTF8String("∅")),

    #### STRUCTURED TEXT
    # I do not find these necessary to support parsers
    # as text is text so do what you will with it.
    # I know there is already json support in julia
    # but what are you going to do with a DataTable
    # full a json or xml?
    114          =>  PostgresType{UTF8String}(:json, UTF8String("{}")),
    3802         =>  PostgresType{UTF8String}(:jsonb, UTF8String("{}")),
    ### PG not built with xml support on my build
    #142  xml     =>  PostgresType{UTF8String}(:xml, UTF8String("")),

    #### TIME
    1082         =>  PostgresType{Date}(:date, Date()),
    # julia only uses 5 decimals of precision for the seconds :(
    #1114         =>  PostgresType{DateTime}(:timestamp, DateTime()),
    #1083 │ time
    #1184 │ timestamptz
    #1266 │ timetz
    #1186 │ interval

    #### OID (INTERNAL) TYPES
    # You never hear of these unless your a Postgres geek.
    #24           =>  PostgresType{Int32}(:regproc, 0),
    #26           =>  PostgresType{Int32}(:oid, 0),
    #2202         =>  PostgresType{Int32}(:regprocedure, 0),
    #2203         =>  PostgresType{Int32}(:regoper, 0),
    #2204         =>  PostgresType{Int32}(:regoperator, 0),
    #2205         =>  PostgresType{Int32}(:regclass, 0),
    #2206         =>  PostgresType{Int32}(:regtype, 0),
    #3734         =>  PostgresType{Int32}(:regconfig, 0),
    #3769         =>  PostgresType{Int32}(:regdictionary, 0),
    #194          =>  PostgresType{UTF8String}(:pg_node_tree, UTF8String("∅")),
    #19           =>  PostgresType{UTF8String}(:name, UTF8String("∅")),
    #3220 │ pg_lsn ;  LSN (Log Sequence Number)

    #### GEOM
    #600 │ point
    #601 │ lseg
    #602 │ path
    #603 │ box
    #604 │ polygon
    #628 │ line
    #718 │ circle

    #### NETWORK
    #650 │ cidr
    #869 │ inet
    #829 │ macaddr; strangely a base user defined type

    #### RANGE
    #3904 │ int4range
    #3906 │ numrange
    #3908 │ tsrange
    #3910 │ tstzrange
    #3912 │ daterange
    #3926 │ int8range

    #MISC.
    #3614 │ tsvector ; text search
    #3615 │ tsquery  ; text serach
    #2950 │ uuid; Universally Unique Identifiers (UUID)

    # Pseudo-Types are abstract so you should
    # should never see their oid's from libpq fetches
)

base_types_byname = Dict([(t[2].name, t[2]) for t in base_types])

max_stars_repo_name: NCarson/Postgres

max_issues_repo_name: NCarson/Postgres

max_forks_repo_name: NCarson/Postgres


text:
# ----------------------------------------------------------------
# Total CO₂ Emissions (accounting for oxidation of CH₄).
# ----------------------------------------------------------------

@defcomp total_co2_emissions begin

	exogenous_CO₂_emissions = Parameter(index=[time]) # Annual carbon dioxide emissions (GtC yr⁻¹).
	oxidized_CH₄_to_CO₂     = Parameter(index=[time]) # Methane of fossil fule ogirin that has been oxidized to carbon dioxide.

	total_CO₂_emissions     = Variable(index=[time])  # Total CO₂ emissions (GtC yr⁻¹).

    function run_timestep(p, v, d, t)

    	# Calculate total CO₂ emissions as the sum of exogenous sources and CO₂ from oxidized CH₄.
    	v.total_CO₂_emissions[t] = p.exogenous_CO₂_emissions[t] + p.oxidized_CH₄_to_CO₂[t]
    end
end


text:
module Infos
using JSON
export Info

# map datatype of python to Julia
const DATATYPE_MAP = Dict{String, DataType}(
    "bool"      => Bool,
    "uint8"     => UInt8,
    "uint16"    => UInt16,
    "uint32"    => UInt32,
    "uint64"    => UInt64,
    "float32"   => Float32,
    "float64"   => Float64
)

module InfoScales

using BigArrays.Codings

export InfoScale

const ENCODING_MAP = Dict{String,Any}(
    # note that the raw encoding in cloud storage will be automatically gzip encoded!
    "raw"       => GzipCoding,
    "jpeg"      => JPEGCoding,
    "blosclz"   => BlosclzCoding,
    "gzip"      => GzipCoding,
    "zstd"      => ZstdCoding
)

mutable struct InfoScale
    key         ::Symbol
    chunkSizes  ::Vector{NTuple{3,Int}}
    encoding    ::DataType
    resolution  ::NTuple{3,Float64}
    volumeSize  ::NTuple{3,Int}
    voxelOffset ::CartesianIndex{3}
end

"""
    InfoScale(; key::Symbol=:4_4_40,
                    chunkSizes::Vector{NTuple{3,Int}}=[(64,64,64)],
                    encoding::DataType=GzipCoding,
                    resolution::NTuple{3,Float64}=map(x->Float64(parse(x)),
                                                      split(string(key), "_")),
                    volumeSize::NTuple{3,Int}=(1024,1024,1024),
                    voxelOffset::CartesianIndex{3}=CartesianIndex{3}(0,0,0))

"""
function InfoScale(; key::Symbol=Symbol("4_4_40"),
                    chunkSizes::Vector{NTuple{3,Int}}=[(64,64,64)],
                    encoding::DataType=GzipCoding,
                    resolution::NTuple{3,Float64}=Tuple(map(x->Float64(Meta.parse(x)),
                                                            split(string(key), "_"))),
                    volumeSize::NTuple{3,Int}=(1024,1024,1024),
                    voxelOffset::CartesianIndex{3}=CartesianIndex{3}(0,0,0))
    InfoScale(key, chunkSizes, encoding, resolution, volumeSize, voxelOffset)
end

function InfoScale(d::Dict{Symbol, Any})
    chunkSizes = map(x->tuple(Int.(x)...), d[:chunk_sizes])
    encoding = ENCODING_MAP[ d[:encoding] ]
    key = d[:key] |> Symbol
    resolution = tuple(Float64.(d[:resolution])...)
    volumeSize = tuple(Int.(d[:size])...)
    voxelOffset = CartesianIndex{3}(Int.(d[:voxel_offset])...)
    InfoScale(key, chunkSizes, encoding, resolution, volumeSize, voxelOffset)
end

function Base.show(io::IO, self::InfoScale)
    show(io, get_key(self))
#    println("chunk size:    ", get_chunk_size(self))
#    println("encoding:      ", get_encoding(self))
#    println("resolution:    ", get_resolution(self))
#    println("volume size:   ", get_volume_size(self))
#    println("voxel offset:  ", get_voxel_offset(self))
end

function Base.Dict(self::InfoScale)
    d = Dict{Symbol,Any}()
    d[:key] = get_key(self)
    d[:chunk_sizes] = [[get_chunk_size(self)...]]
    for (k,v) in ENCODING_MAP
        if v == get_encoding(self)
            d[:encoding] = string(k)
        end
    end
    d[:size] = [get_volume_size(self)...]
    d[:voxel_offset] = [Tuple(get_voxel_offset(self))...]
    return d
end

@inline function Base.string(self::InfoScale)
    d = Dict(self)
    JSON.json(d)
end

@inline function get_key(self::InfoScale)
    self.key
end

@inline function get_chunk_size(self::InfoScale)
    @assert length(self.chunkSizes) == 1
    self.chunkSizes[1]
end

@inline function set_chunk_size!(self::InfoScale, chunkSize::NTuple{3,Int})
    @assert length(self.chunkSizes) == 1
    self.chunkSizes[1]=chunkSize
end

@inline function get_offset(self::InfoScale)
    self.voxelOffset
end
@inline function set_offset(self::InfoScale, offset::CartesianIndex{N}) where N
    self.voxelOffset = offset
end

@inline function get_encoding(self::InfoScale) self.encoding end
@inline function set_encoding!(self::InfoScale, encoding::DataType)
    self.encoding = encoding
end

"""
    set_encoding!(self::InfoScale, encoding::Symbol)

the encoding map is:
$(ENCODING_MAP)
"""
@inline function set_encoding!(self::InfoScale, encoding::String)
    self.encoding = ENCODING_MAP[encoding]
end

@inline function get_resolution(self::InfoScale) self.resolution end
@inline function set_resolution!(self::InfoScale, resolution::NTuple{3,T}) where T
    self.resolution = map(Float64, resolution)
end

@inline function get_volume_size(self::InfoScale) self.volumeSize end
@inline function set_volume_size!(self::InfoScale, volumeSize::NTuple{3,Int})
    self.volumeSize=volumeSize
end

@inline function get_voxel_offset(self::InfoScale) self.voxelOffset end
@inline function set_voxel_offset!(self::InfoScale, voxelOffset::CartesianIndex{3})
    self.voxelOffset = voxelOffset
end

function get_properties(self::InfoScale)
    chunkSize = get_chunk_size(self)
    encoding = get_encoding(self)
    resolution = get_resolution(self)
    voxelOffset = get_voxel_offset(self)
    volumeSize = get_volume_size(self)
    return chunkSize, encoding, resolution, voxelOffset, volumeSize
end

"""
    generate_next_mip(self::InfoScale)

only downsample the images in XY plane by 2 times.
meaning increase the resolution by 2x in X and Y axis
the chunk size will remain the same
the encoding will also be the same
"""
function generate_next_mip(self::InfoScale)
    resolution = (self.resolution[1]*2.0, self.resolution[2]*2.0, self.resolution[3])
    key = Symbol("$(round(Int,resolution[1]))_$(round(Int,resolution[2]))_$(round(Int,resolution[3]))")
    chunkSize = get_chunk_size(self)
    encoding = get_encoding(self)
    volumeSize = map(div, get_volume_size(self), (2,2,1))
    voxelOffset = CartesianIndex(map(div, Tuple(get_voxel_offset(self)), (2,2,1)))

    return InfoScale(key, [chunkSize], encoding, resolution, volumeSize, voxelOffset)
end

end # end of InfoScales module

using .InfoScales

mutable struct Info{T}
    mesh        ::String
    numChannels ::Int
    scales      ::Vector{InfoScale}
    skeletons   ::String
    layerType   ::Symbol
    function Info(dataType::DataType, mesh::String,
                    numChannels::Int, scales::Vector{InfoScale},
                    skeletons::String, layerType::Symbol)
        new{dataType}(mesh, numChannels, scales, skeletons, layerType)
    end
end

"""
    Info(;
        dataType::DataType=UInt8,
        mesh::String="",
        numChannels::Int=1,
        scales::Vector{InfoScale}=[InfoScale()],
        skeletons::String="",
        layerType::Symbol=:image)

layerType: the layer type defined in neuroglancer precomputed format: {image, segmentation}
"""
function Info(;
                dataType::DataType=UInt8,
                mesh::String="",
                numChannels::Int=1,
                scales::Vector{InfoScale}=[InfoScale()],
                skeletons::String="",
                layerType::Symbol=:image,
                volumeSize::NTuple{3,Int}=(1024,1024,1024),
                chunkSize::NTuple{3,Int}=(64,64,64),
                voxelOffset::CartesianIndex{3}=CartesianIndex{3}((0,0,0)),
                numMip::Int=1)
    infoScale = scales[1]
    InfoScales.set_volume_size!(infoScale, volumeSize)
    InfoScales.set_chunk_size!(infoScale, chunkSize)
    InfoScales.set_voxel_offset!(infoScale, voxelOffset)

    for mip in 2:numMip
        infoScale = InfoScales.generate_next_mip(infoScale)
        push!(scales, infoScale)
    end

    Info(dataType, mesh, numChannels, scales, skeletons, layerType)
end

function Info(d::Dict{Symbol,Any})
    dataType = DATATYPE_MAP[ d[:data_type] ]
    mesh = get(d, :mesh, "")
    numChannels = d[:num_channels]
    scales = map(InfoScale, d[:scales])
    skeletons = get(d, :skeletons, "")
    layerType = Symbol( d[:type] )
    Info(dataType, mesh, numChannels, scales, skeletons, layerType)
end

@inline function Info(str::AbstractString)
    d = JSON.parse(str, dicttype=Dict{Symbol, Any})
    Info(d)
end

@inline function Info(data::Vector{UInt8})
    Info(String(data))
end

function Base.show(self::Info{T}) where T
    println("\ndata type:   ", T)
    println("mesh:          ", get_mesh(self))
    println("num of channel:", get_num_channels(self))
    println("scales:        ", get_scales(self))
    println("skeletons:     ", get_skeletons(self))
    println("layer type:    ", get_layer_type(self))
end

"""
    Base.Dict(self::Info)

the transformation follows JSON format
"""
function Base.Dict(self::Info{T}) where T
    d = Dict{Symbol, Any}()
    for (k,v) in DATATYPE_MAP
        if v == T
            d[:data_type] = string(k)
        end
    end

    if !isempty(get_mesh(self))
        d[:mesh] = get_mesh(self)
    end

    d[:num_channels] = get_num_channels(self)
    d[:scales] = map(Dict, get_scales(self))
    d[:skeletons] = string(get_skeletons(self))
    d[:type] = string(get_layer_type(self))
    return d
end

"""
    Base.string(self::Info)

"""
@inline function Base.string(self::Info)
    d = Dict(self)
    JSON.json(d)
end

############ get the properties ##########
@inline function get_mesh(self::Info) self.mesh end
@inline function set_mesh!(self::Info, mesh::String) self.mesh=mesh end

@inline function get_key(self::Info, mip::Integer=1)
    InfoScales.get_key( get_scales(self)[mip] )
end

@inline function get_chunk_size(self::Info, mip::Integer=1)
    InfoScales.get_chunk_size( self.scales[mip] )
end
@inline function set_chunk_size!(self::Info, chunkSize::NTuple{3,Int}, mip::Integer)
    InfoScales.set_chunk_size!( self.scales[mip], chunkSize)
end
@inline function set_chunk_size!(self::Info, chunkSize::NTuple{3,Int})
    for infoScale in self.scales
        InfoScales.set_chunk_size!( infoScale )
    end
end

@inline function get_offset(self::Info, mip::Integer=1)
    InfoScales.get_offset( self.scales[mip] )
end
@inline function set_offset!(self::Info{T}, offset::CartesianIndex{N}, mip::Integer) where {T,N}
    InfoScales.set_offset!( self.scales[mip], offset)
end

@inline function get_volume_size(self::Info, mip::Integer=1)
    InfoScales.get_volume_size(self.scales[mip])
end
@inline function set_volume_size!(self::Info, volumeSize::NTuple{3,Int}, mip::Integer=1)
    InfoScales.set_volume_size!(self.scales[mip], volumeSize)
end

@inline function get_encoding(self::Info, mip::Integer=1)
    InfoScales.get_encoding( get_scales(self)[mip] )
end
@inline function set_encoding!(self::Info, encoding::Union{String, DataType})
    for scale in get_scales(self)
        InfoScales.set_encoding!(scale, encoding)
    end
    nothing
end

@inline function get_num_channels(self::Info) self.numChannels end
@inline function set_num_channels!(self::Info, numChannels::Int) self.numChannels=numChannels end

@inline function get_data_type(info::Info{T}) where T
    T
end

@inline function Base.size(self::Info)
    get_volume_size(self, 0)
end

function Base.ndims(info::Info)
    numChannels = get_num_channels(info)
    if numChannels == 1
        return 3
    else
        return 4
    end
end

@inline function get_scales(self::Info) self.scales end
@inline function set_scales!(self::Info, scales::Vector{InfoScale}) self.scales=scales end

@inline function get_skeletons(self::Info) self.skeletons end
@inline function set_skeletons!(self::Info, skeletons::String) self.skeletons=skeletons end

@inline function get_layer_type(self::Info) self.layerType end
@inline function set_layer_type(self::Info, layerType::Symbol) self.layerType=layerType end

"""
    get_properties_in_mip_level(self::Info, mip::Integer=1)
"""
function get_properties_in_mip_level(self::Info, mip::Integer=1)
    numChannels = get_num_channels(self)
    infoScale = self.scales[mip]
    chunkSize, encoding, resolution, voxelOffset, volumeSize =
                                        InfoScales.get_properties(infoScale)
    if numChannels > 1
        chunkSize = (chunkSize..., numChannels)
        volumeSize = (volumeSize..., numChannels)
        voxelOffset = CartesianIndex(Tuple(voxelOffset)..., 0)
    end
    return chunkSize, encoding, resolution, voxelOffset, volumeSize
end

"""
    get_properties_in_mip_level(self::Info, key::Symbol)
"""
function get_properties_in_mip_level(self::Info, key::Symbol)
    for (i, infoScale) in self.scales |> enumerate
        if key == InfoScales.get_key(infoScale)
            return get_properties_in_mip_level(self, i)
        end
    end
    @error "did not find any corresponding mip level with key: " key
end

## TODO: some functions to create info file

end # end of module