text:
import Base.show
import Base.==
import Base.isequal
import SymPy.solve
import SymPy.degree
import SymPy.subs

coeff(x::Sym, y::Sym) = x.coeff(y)
is_polynomial(x::Sym, y::Sym) = x.is_polynomial(y)
Poly(x::Sym, y::Sym) = sympy.Poly(x, y)

abstract type Recurrence end

struct CFiniteRecurrence <: Recurrence
    coeffs::Array{Sym}
    f::SymFunction
    n::Sym
    inhom::Sym
end

CFiniteRecurrence(coeffs::Array{Sym}, f::SymFunction, n::Sym) = CFiniteRecurrence(coeffs, f, n, Sym(0))
CFiniteRecurrence(coeffs::Array{Sym}, f::SymFunction, n::Symbol) = CFiniteRecurrence(coeffs, f, Sym(n))

function order(r::Recurrence)
    return length(r.coeffs) - 1
end

function relation(r::Recurrence)
    return hompart(r) + r.inhom
end

function Base.show(io::IO, r::Recurrence)
    print(io, "$(lhs(r)) = $(rhs(r))")
end

function hompart(r::Recurrence)
    return sum([c * r.f(r.n + (i - 1)) for (i, c) in enumerate(r.coeffs)])
end

function is_homogeneous(r::Recurrence)
    return simplify(r.inhom) == 0
end

function rhs(r::Recurrence)
    if length(r.coeffs) == 1
        return -r.inhom
    end
    simplify(-sum([c * r.f(r.n + (i - 1)) for (i, c) in enumerate(r.coeffs[1:end-1])]) - r.inhom)
end

function lhs(r::Recurrence)
    t = r.f(r.n + order(r))
    if length(r.coeffs) > 1
        t *= r.coeffs[end]
    end
    t
end

function homogeneous(r::Recurrence)
    n = r.n
    if is_homogeneous(r)
        return r
    elseif is_polynomial(r.inhom, r.n)
        d = convert(Int64, degree(Poly(r.inhom, n)))
        res = relation(r)
        for i in 1:d + 1
            res = (res |> subs(n, n + 1)) - res
        end
        res = simplify(res)

        idx = n + order(r) + d + 1
        c = coeff(res, r.f(idx))
        return eq2rec(simplify(res / c), r.f, r.n)
    else
        res = hompart(r) / r.inhom
        res = expand((res |> subs(n, n + 1)) - res)
        idx = n + order(r) + 1
        c = coeff(res, r.f(idx))
        return eq2rec(simplify(res / c), r.f, r.n)
    end
end

function rewrite(d::Dict{Sym,Sym})
    # TODO: successively apply rewrite rules do values of dict
    dd = Dict{Sym,Sym}()
    while dd != d
        dd = d
        for (k, v) in d
            for (g, h) in dd
                dd[g] = h |> subs(k, v)
            end
        end
    end
    dd
end

function closedform(orig::CFiniteRecurrence)
    # println("Original: ", orig)
    r = homogeneous(orig)
    # println("Homogeneous: ", r)

    shift = order(r) - order(orig)
    rh = rhs(orig)
    ord = order(orig)
    init = Dict{Sym,Sym}([(orig.f(i), SymPy.replace(rh, orig.n, i-ord)) for i in ord:shift+ord-1])
    # init = rewrite(init)

    rel = relation(r)
    # println("Homogeneous: ", homogeneous(r))
    w0  = Wild("w0")
    @syms lbd
    cpoly   = SymPy.replace(rel, r.f(r.n), 1)
    cpoly   = SymPy.replace(cpoly, r.f(r.n + w0), lbd^w0)
    # println("CPoly: ", cpoly |> simplify)
    # factors = factor_list(cpoly)
    proots   = roots(cpoly)
    # d = hcat([[uniquevar() * r.n^i, z^r.n] for (z, m) in roots for i in 0:m - 1])
    # println(d[:,2])
    ansatz = sum([sum([uniquevar() * r.n^i * z^r.n for i in 0:m - 1]) for (z, m) in proots])
    # println(ansatz)
    # println(free_symbols(ansatz(n)))
    unknowns = filter(e -> e != r.n, free_symbols(ansatz))
    system = [Eq(r.f(i), ansatz |> subs(r.n, i)) for i in 0:order(r) - 1]
    sol = solve(system, unknowns)
    sol = ansatz |> subs(sol)
    if !isempty(init)
        tmp = nothing
        while true
            tmp = (sol |> subs(init)) |> simplify
            if tmp == sol
                break
            end
            sol = tmp

        end
        sol = simplify(tmp)
    end

    exp = [z for (z, _) in proots]
    exp = filter(x -> x!=Sym(1), exp)
    push!(exp, Sym(1))
    coeff = exp_coeffs(sol, [z^r.n for z in exp])
    return CFiniteClosedForm(r.f, r.n, exp, coeff)
end

function exp_coeffs(expr::Sym, exp::Array{Sym,1})
    # assume if 1 in exp then it is at the end
    coeffs = Sym[]
    for ex in exp
        if ex == Sym(1)
            push!(coeffs, expr)
        else
            c = coeff(expr, simplify(ex))
            push!(coeffs, c)
            expr = simplify(expr - c*ex)
        end
    end
    coeffs
end

function rec_dependency(recs::Array{<: Recurrence,1})
    rec_dependency([], [], recs...)
end

function rec_dependency(indep, dep, rec::Recurrence, recs::Recurrence...)
    expr = rhs(rec)
    fns = symfunctions(expr)
    depvars = filter(x -> (x.func!=Sym(rec.f.x) && !in(Sym(0), x.args)), fns)

    if isempty(depvars)
        push!(indep, rec)
    else
        push!(dep, (rec, depvars))
    end
    rec_dependency(indep, dep, recs...)
end

function rec_dependency(indep, dep)
    indep, dep
end

function subs(cf::ClosedForm, r::CFiniteRecurrence)
    rel = relation(r)
    rhs = poly(cf)
    fns = symfunctions(rel)
    for fn in fns
        if fn.func == Sym(cf.f.x)
            w0 = Wild("w0")
            idx = match(cf.f(cf.n + w0), fn)[w0]
            sub = rhs |> SymPy.subs(cf.n, cf.n + idx)
            rel = rel |> SymPy.subs(cf.f(cf.n + idx), sub) |> simplify
        end
    end
    return eq2rec(rel, r.f, r.n)
end

function rec_solve(recsorig::Array{<: Recurrence,1})
    # println("Recurrences: ", recs)
    if length(recsorig) == 0
        return ClosedForm[]
    end

    varlist = [(r.f(r.n), r.f) for r in recsorig]
    recs, varlist = rec_simplify(Recurrence[], varlist, recsorig...)
    # println("Simplifed recurrences: ", recs)
    # println("Varlist: ", varlist)
    # recs = loop.body
    solved = ClosedForm[]
    solvable = true
    while solvable
        indep, dep = rec_dependency(recs)
        indepcnt = length(indep)
        depcnt = length(dep)

        if indepcnt == 0 && depcnt != 0
            error("Illegal coupling detected in loop body")
        elseif indepcnt == 0 && depcnt == 0
            break
        end

        recs = setdiff(recs, indep)
        for rec in indep
            cf = closedform(rec)

            # TODO: deal with the case when a recurrence is not solvable

            push!(solved, cf)
            recs = filter(x -> x!=rec, recs)

            recs = [subs(cf, r) for r in recs]
        end
        if isempty(recs)
            break
        end
    end

    inivars = []
    for cf in solved
        inivars = [inivars; initvars(cf)]
    end
    inivars = [(SymFunction(string(Sym(x.func))), x.args[1]) for x in inivars]
    inivars = filter(x -> x[2]>0, inivars)
    iniexpr = init_expr(recsorig, inivars)
    # println("Initial expression: ", iniexpr)

    for var in varlist
        # TODO: assume var is something like (t(n-1), x)
        fn = SymFunction(string(var[1].func))
        if fn != var[2]
            for cf in solved
                if cf.f == fn
                    sh = args(var[1])[1] - cf.n
                    tmp = shift(cf, sh)
                    tmp.f = var[2]
                    push!(solved, tmp)
                    break
                end
            end
        end
    end
    [subs!(cf, iniexpr...) for cf in solved]
end

function ==(f::SymFunction, g::SymFunction)
    isequal(Sym(f.x), Sym(g.x))
end

function init_expr(recs::Array{<:Recurrence,1}, vars::Array{Tuple{SymFunction,Sym},1})
    initrules = Pair[]
    for (var, idx) in vars
        rec = [r for r in recs if string(r.f) == string(var)][1]
        rhs = initial(rec, idx)
        while true
            newvars = [(SymFunction(string(Sym(t.func))), t.args[1]) for t in symfunctions(rhs) if t.args[1] > 0]
            # newvars = filter(x -> x[2]>0, newvars)
            if isempty(newvars)
                break
            end
            rules = init_expr(recs, newvars)
            rhs = (rhs |> subs(rules...)) |> simplify
        end
        push!(initrules, Pair(var(idx), rhs))
    end
    initrules
end

function initial(r::CFiniteRecurrence, n::Sym)
    if n < order(r)
        error("Cannot give expression for $(r.f(n))")
    end
    i = n - order(r)
    rhs(r) |> subs(r.n, i)
end

function initvars(r::ClosedForm)
    symfunctions(rhs(r))
end

#-------------------------------------------------------------------------------

function contains_func(expr::Sym, f::SymFunction)
    fns = symfunctions(expr)
    w0 = Wild("w0")
    args = [get(match(f(w0), fn), w0, nothing) for fn in fns]
    args = filter(x -> x!=nothing, args)
    !isempty(args)
end

function rec_subs(r::Recurrence, rec::Recurrence)
    # println("subs $(rec) in $(r)")
    expr = rhs(r)
    fns = symfunctions(expr)
    w0 = Wild("w0")
    args = [get(match(rec.f(w0), fn), w0, nothing) for fn in fns]
    args = filter(x -> x!=nothing, args)
    if isempty(args)
        return r
    end
    for arg in args
        diff = rec.n + 1 - arg
        sub = rhs(rec) |> subs(rec.n, rec.n - diff)
        expr = expr |> subs(rec.f(arg), sub)
    end
    # assume lhs of rec is always f(n+1) (just temporarily)
    # println("result: ", expr)
    eq2rec(r.f(r.n+1) - simplify(expr), r.f, r.n)
end

function rec_subs(expr::Sym, rec::Recurrence)
    # println("subs $(rec) in $(r)")
    # expr = rhs(r)
    fns = symfunctions(expr)
    w0 = Wild("w0")
    args = [get(match(rec.f(w0), fn), w0, nothing) for fn in fns]
    args = filter(x -> x!=nothing, args)
    for arg in args
        diff = rec.n + 1 - arg
        sub = rhs(rec) |> subs(rec.n, rec.n - diff)
        expr = expr |> subs(rec.f(arg), sub)
    end
    simplify(expr)
end

# function rec_simplify(recs::Array{<:Recurrence,1})
#     # processed = Recurrence[]
#     # for rec in recs
#     #     expr = rhs(rec)
#     #     # fns = symfunctions(expr)
#     #     # w0 = Wild("w0")
#     #     # args = [match(rec.f(w0), expr)[w0] for fn in fns]
#     #     if !contains_func(expr, rec.f) # not a real recurrence
#     #         processed = [rec_subs(r, rec) for r in processed]
#     #         recs = [rec_subs(r, rec) for r in recs]
#     #     else
#     #         push!(processed, rec)
#     #     end
#     # end
#     # processed
#     rec_simplify(Recurrence[], recs...)
# end

function rec_simplify(processed::Array{<:Recurrence}, varlist::Array{Tuple{Sym,SymFunction},1}, rec::Recurrence, recs::Recurrence...)
    expr = rhs(rec)
    # fns = symfunctions(expr)
    # w0 = Wild("w0")
    # args = [match(rec.f(w0), expr)[w0] for fn in fns]
    if !contains_func(expr, rec.f) # not a real recurrence
        varlist = [var_simplify(var, rec) for var in varlist]
        processed = [rec_subs(r, rec) for r in processed]
        recs = [rec_subs(r, rec) for r in recs]
    else
        push!(processed, rec)
    end
    rec_simplify(processed, varlist, recs...)
end

function rec_simplify(processed::Array{<:Recurrence}, varlist::Array{Tuple{Sym,SymFunction},1})
    processed, varlist
end

function var_simplify(var::Tuple{Sym,SymFunction}, rec::Recurrence)
    rec_subs(var[1], rec), var[2]
end


text:
#=
const sources = [""]

@compat if is_windows()
    flags = ["-fopenmp", "-Wall", "-O3", "-shared", "-march=native"]
    libname = "libmerlin.dll"
    cmd = `nvcc `
    run(cmd)
elseif is_apple()
    flags = ["-fPIC", "-Wall", "-O3", "-shared", "-march=native"]
    libname = "libmerlin.so"
    cmd = `$compiler $flags -o $libname $sources`
    println("Running $cmd")
    run(cmd)
elseif is_linux()
    flags = ["-fopenmp", "-fPIC", "-Wall", "-O3", "-shared", "-march=native"]
    libname = "libmerlin.so"
    cmd = `$compiler $flags -o $libname $sources`
    println("Running $cmd")
    run(cmd)
else
    throw("Unknown OS.")
end
=#

ext: jl


text:
using Documenter, LatinHypercubeSampling, Random

makedocs(;
        modules = [LatinHypercubeSampling],
        format = Documenter.HTML(   assets = ["assets/favicon.ico"],
                                prettyurls = get(ENV, "CI", nothing) == "true"),
        sitename = "LatinHypercubeSampling.jl",
        strict = true,
        clean = true,
        checkdocs = :none,
        pages = [
                "Home" => "index.md",
                "Manual" => [
                        "man/lhcoptim.md",
                        "man/sublhcoptim.md",
                        "man/categorical.md",
                ]
        ]
)

deploydocs(
    repo = "github.com/MrUrq/LatinHypercubeSampling.jl.git"
)

max_stars_repo_name: MrUrq/LatinHypercubeSampling

max_issues_repo_name: MrUrq/LatinHypercubeSampling

max_forks_repo_name: MrUrq/LatinHypercubeSampling


text:
############################################################
# FactCheck.jl
# A testing framework for Julia
# http://github.com/JuliaLang/FactCheck.jl
# MIT Licensed
############################################################
# Assertion helpers
# - not
# - exactly
# - roughly
# - anyof
# - less_than
# - less_than_or_equal
# - greater_than_or_equal
############################################################

# not: logical not for values and functions
not(x) = isa(x, Function) ? (y) -> !x(y) :
                            (y) -> x != y

# exactly: tests object/function equality (i.e. ===)
exactly(x) = (y) -> x === y

# approx/roughly: Comparing numbers approximately
roughly(x::Number, atol) = (y::Number) -> isapprox(y, x, atol=atol)
roughly(x::Number; kvtols...) = (y::Number) -> isapprox(y, x; kvtols...)

roughly(A::AbstractArray, atol) = (B::AbstractArray) -> begin
    size(A) != size(B) && return false
    return isapprox(A, B, atol=atol)
end
roughly(A::AbstractArray; kvtols...) = (B::AbstractArray) -> begin
    size(A) != size(B) && return false
    return isapprox(A, B; kvtols...)
end

# anyof: match any of the arguments
anyof(x...) = y -> any(arg->(isa(arg,Function) ? arg(y)::Bool : (y==arg)), x)

# less_than: Comparing two numbers
less_than(compared) = (compare) -> compare < compared

# less_than_or_equal: Comparing two numbers
less_than_or_equal(compared) = (compare) -> compare <= compared

# greater_than: Comparing two numbers
greater_than(compared) = (compare) -> compare > compared

# greater_than_or_equal: Comparing two numbers
greater_than_or_equal(compared) = (compare) -> compare >= compared


text:
module EvaluateOdds

include("./gans.jl")
include("./vaes.jl")
include("./micegen.jl")
import JLD
using Statistics
import Random
using Distributed
import BoltzmannMachines
using SplitApplyCombine
const BMs = BoltzmannMachines

using LinearAlgebra

function logoddsratio(x)
   N, dim = size(x)
   OR = LowerTriangular(ones(dim,dim))
   for i = 1:dim
      for j = 1:(i-1)
         a = N - sum(min.(x[:,i] .+ x[:,j], 1))      # ij_00
         diff = x[:,i] .- x[:,j]
         b = sum(diff .== -1)                        # ij_01
         c = sum(diff .== 1)                         # ij_10
         d = sum(x[:,i] .* x[:,j])                   # ij_11
         OR[i,j] = max(a,0.5)*max(d,0.5)/max(b,0.5)/max(c,0.5)
      end
   end
   return LowerTriangular(log.(OR))
end

function oddsdelta(x1, x2)
    mean((logoddsratio(x1) .- logoddsratio(x2)).^2)
end

function my_fitdbm(x::Matrix{Float64}, epochs::Int)
   dbm = BMs.fitdbm(x;
         learningrate = 0.001,
         pretraining = [
               BMs.TrainLayer(nhidden = size(x, 2));
               BMs.TrainLayer(nhidden = 10)],
         epochspretraining = epochs,
         epochs = 20)
end

function getdata(datasetindex)
   JLD.load("data/$(datasetindex)_50.jld")["$(datasetindex)_50"]
end

const nvariables = 50
const ntest = 100
const nsamples = 500

function validationdataset(i)
   getdata(i)[1001:2000, 1:nvariables]
end

const datasetids = 1:30
const validationoddsratios = Dict([(i, logoddsratio(validationdataset(i))) for i in datasetids])
const training_datasets = Dict([(i, getdata(i)[1:nsamples, 1:nvariables]) for i in datasetids])
const test_datasets = Dict([(i, getdata(i)[(nsamples+1):(nsamples + ntest), 1:nvariables])
      for i in datasetids])

# Collect the data sets from the different sites and evaluate the results
function evaluatexgen(fitgenresults)
   xgens = map(r -> r[:xgen], fitgenresults)
   xgens = filter(x -> length(size(x)) == 2, xgens)
   xgen = vcat(xgens...)
   if isempty(xgen)
      delta = Inf
   else
      validationodds = validationoddsratios[fitgenresults[1][:datasetindex]]
      delta = mean((logoddsratio(xgen) .- validationodds).^2)
   end

   oddsdeltatrain = minimum(map(r -> r[:oddsdeltatrain], fitgenresults))

   # The different sites have been trained with the same model, the same seed,
   # and the same number of epochs, so it suffices to take the value of the
   # first array element.
   (seed = fitgenresults[1][:seed],
         datasetindex = fitgenresults[1][:datasetindex],
         modeltype = fitgenresults[1][:modeltype],
         epochs = fitgenresults[1][:epoch],
         nsites = fitgenresults[1][:nsites],
         oddsdeltatrain = oddsdeltatrain,
         oddsdelta = delta)
end

function evaluateoddsratios(modeldefinitions, seeds)
   nsites = [1,2,5,20]

   args = Iterators.product(
      seeds,
      datasetids,
      nsites,
      modeldefinitions)

   pmap(EvaluateOdds.evaluateodds, args; retry_delays = ExponentialBackOff(n = 5))
end

function evaluateodds((seed, datasetindex, nsites, modeldefinition))
   evaluatexgen(map(datasetpart ->
      fitgenerate(seed, datasetindex, nsites, datasetpart, modeldefinition),
         1:nsites))
end

function _initvae()
   VAEs.initvae(nvariables = nvariables, nhidden = nvariables, nz = 10)
end

function _trainvae!(model, x, epochs)
   VAEs.trainvae!(model, x; epochs = epochs, learningrate = 0.002)
end

function _initgan()
   GANs.initgan(nvariables)
end

function _traingan!(model, x, epochs)
   GANs.traingan!(model, x; int_epochs = epochs)
end

function extract_trainingdata(datasetindex, nsites, datasetpart)
   xtrain = training_datasets[datasetindex]
   npersite = floor(Int, nsamples / nsites)
   xpart = xtrain[((datasetpart-1)*npersite+1):(datasetpart*npersite),:]
end

function fitgenerate(
      seed,
      datasetindex,
      nsites,
      datasetpart,
      modeldefinition)

   function evaluation((model_evaluation, xgen, oddsdeltatrain, epoch))
      (seed = seed, datasetindex = datasetindex,
         model_evaluation = model_evaluation,
         nsites = nsites, datasetpart = datasetpart,
         modeltype = modeldefinition["model"],
         epoch = epoch,
         xgen = xgen,
         oddsdeltatrain = oddsdeltatrain)
   end

   Random.seed!(seed)

   xtest = test_datasets[datasetindex]
   xpart = extract_trainingdata(datasetindex, nsites, datasetpart)

   ngensamples = size(xpart, 1)

   if modeldefinition["model"] == "dbm"
      return evaluation(evaldbm(modeldefinition, xpart, xtest))

   elseif modeldefinition["model"] == "mice"
      return evaluation(evalmice(modeldefinition, xpart, xtest))

   elseif modeldefinition["model"] == "gan"
      return evaluation(eval_multepochs(modeldefinition, xpart, xtest;
         initfun = _initgan,
         trainfun! = _traingan!,
         samplefun = model -> GANs.samples(model, ngensamples)))

   elseif modeldefinition["model"] == "vae"
      return evaluation(eval_multepochs(modeldefinition, xpart, xtest;
         initfun = _initvae,
         trainfun! = _trainvae!,
         samplefun = model -> VAEs.samples(model, ngensamples)))

   elseif modeldefinition["model"] == "im"
      return evaluation(evalim(modeldefinition, xpart, xtest))

   else
      error("unknown model")
   end
end

function evaldbm(modeldefinition, xpart, xtest)
   ngensamples = size(xpart, 1)
   epochs = modeldefinition["epochs"]
   dbm = my_fitdbm(xpart, epochs)
   xgen = BMs.samples(dbm, ngensamples)
   model_evaluation = oddsdelta(xgen, xtest)
   oddsdeltatrain = oddsdelta(xpart, xgen)
   model_evaluation, xgen, oddsdeltatrain, epochs
end

function evalmice(modeldefinition, xpart, xtest)
   ngensamples = size(xpart, 1)
   model_evaluation = Inf
   oddsdeltatrain = Inf
   xgen = []
   try
      model = MICE.fitmice(xpart)
      xgen = MICE.genmice(model, ngensamples)
      model_evaluation = oddsdelta(xgen, xtest)
      oddsdeltatrain = oddsdelta(xpart, xgen)
   catch ex
   end
   epoch = 0 # not relevant for mice
   model_evaluation, xgen, oddsdeltatrain, epoch
end

function evalim(modeldefinition, xpart, xtest)
   ngensamples = size(xpart, 1)
   # odds ratio of independent variables is 1
   varmeans = mean(xpart, dims = 1)
   xgen = float.(rand(Float64, size(xpart)) .< varmeans)
   model_evaluation = oddsdelta(xgen, xtest)
   oddsdeltatrain = oddsdelta(xpart, xgen)
   epoch = 0 # not relevant for im
   model_evaluation, xgen, oddsdeltatrain, epoch
end

function eval_multepochs(modeldefinition, xpart, xtest;
      initfun::Function, trainfun!::Function, samplefun::Function)

   epochs = vcat(0, modeldefinition["epochs"])
   bestoddsdelta = Inf
   bestxgen = zeros(1,1)
   bestepoch = 0

   model = initfun()

   for i in 2:length(epochs)
      trainfun!(model, xpart, epochs[i] - epochs[i-1])
      # let sampling not change the random numbers for model training
      rng = copy(Random.default_rng());
      xgen = samplefun(model)
      copy!(Random.default_rng(), rng);
      newoddsdelta = oddsdelta(xgen, xtest)
      if newoddsdelta < bestoddsdelta
         bestoddsdelta = newoddsdelta
         bestxgen = xgen
         bestepoch = epochs[i]
      end
   end
   model_evaluation = bestoddsdelta
   xgen = bestxgen
   oddsdeltatrain = oddsdelta(xpart, xgen)
   model_evaluation, xgen, oddsdeltatrain, bestepoch
end

end

max_stars_repo_name: stefan-m-lenz/dist-gen-comp

max_issues_repo_name: stefan-m-lenz/dist-gen-comp

max_forks_repo_name: stefan-m-lenz/dist-gen-comp


text:
@testset "icd9" begin
    icd9_add_dots    = DiagnosisClassification.icd9_add_dots
    icd9_remove_dots = DiagnosisClassification.icd9_remove_dots

    @testset "icd9_add_dots" begin
        test_cases = [
            # Input is too short for a dot
            "003"    => "003",
            "V01"    => "V01",
            "E000"   => "E000",

            # Input contains a dot
            "003.2"  => "003.2",
            "003.20" => "003.20",
            "V01.7"  => "V01.7",
            "V01.71" => "V01.71",
            "E000.0" => "E000.0",

            # Input does not contain a dot
            "0032"   => "003.2",
            "00320"  => "003.20",
            "V017"   => "V01.7",
            "V0171"  => "V01.71",
            "E0000"  => "E000.0",
        ]
        for test_case in test_cases
            input           = test_case[1]
            expected_output = test_case[2]
            @test icd9_add_dots(input) == expected_output
        end
        @test_throws ErrorException icd9_add_dots("")
    end

    @testset "icd9_remove_dots" begin
        test_cases = [
            # Input is too short for a dot
            "003"    => "003",
            "V01"    => "V01",
            "E000"   => "E000",

            # Input contains a dot
            "003.2"  => "0032",
            "003.20" => "00320",
            "V01.7"  => "V017",
            "V01.71" => "V0171",
            "E000.0" => "E0000",

            # Input does not contain a dot
            "0032"   => "0032",
            "00320"  => "00320",
            "V017"   => "V017",
            "V0171"  => "V0171",
            "E0000"  => "E0000",
        ]
        for test_case in test_cases
            input           = test_case[1]
            expected_output = test_case[2]
            @test icd9_remove_dots(input) == expected_output
        end
        @test_throws ErrorException icd9_remove_dots("")
    end

    @testset "Roundtrip icd9_add_dots" begin
        test_cases = [
            # Input is too short for a dot
            "003"    => "003",
            "V01"    => "V01",
            "E000"   => "E000",

            # Input contains a dot
            "003.2"  => "0032",
            "003.20" => "00320",
            "V01.7"  => "V017",
            "V01.71" => "V0171",
            "E000.0" => "E0000",

            # Input does not contain a dot
            "0032"   => "0032",
            "00320"  => "00320",
            "V017"   => "V017",
            "V0171"  => "V0171",
            "E0000"  => "E0000",
        ]
        for test_case in test_cases
            input           = test_case[1]
            expected_output = test_case[2]
            @test icd9_remove_dots(icd9_add_dots(input)) == expected_output
        end
    end

    @testset "Roundtrip icd9_remove_dots" begin
        test_cases = [
            # Input is too short for a dot
            "003"    => "003",
            "V01"    => "V01",
            "E000"   => "E000",

            # Input contains a dot
            "003.2"  => "003.2",
            "003.20" => "003.20",
            "V01.7"  => "V01.7",
            "V01.71" => "V01.71",
            "E000.0" => "E000.0",

            # Input does not contain a dot
            "0032"   => "003.2",
            "00320"  => "003.20",
            "V017"   => "V01.7",
            "V0171"  => "V01.71",
            "E0000"  => "E000.0",
        ]
        for test_case in test_cases
            input           = test_case[1]
            expected_output = test_case[2]
            @test icd9_add_dots(icd9_remove_dots(input)) == expected_output
        end
    end
end


text:
module ExceptionTests

include("mpiexec.jl")
run_mpi_driver(procs=8,file="driver_exception.jl")

end # module


text:
function _debias1(Ψx, Ψy, θ, ω, θ_ind::Int; refit=true)
    μx = vec(mean(Ψx, dims=2))
    if refit == false
        r = rhat(θ, Ψy)
        θ1 = θ[θ_ind]
        for l in findall(!iszero, ω)
            θ1 += ω[l] * ( μx[l] - mean( r .* Ψy[l, :] ) )
        end
        return θ1
    else
        supp = _find_supp(θ_ind, θ)
        θk = KLIEP(Ψx[supp, :], Ψy[supp, :], CD_KLIEP())

        r = rhat(θk, Ψy[supp, :])
        θ1 = θk[end]
        for l in findall(!iszero, ω)
            θ1 += ω[l] * ( μx[l] - mean( r .* Ψy[l, :] ) )
        end
        return θ1
    end
end

function _debias1(Ψx, Ψy, θ, Hinv, θ_ind::Union{Vector{Int},UnitRange}; refit=true)
    θ1 = Vector{Float64}(undef, length(θ_ind))
    μx = vec(mean(Ψx, dims=2))
    if refit == false
        for k in 1:length(θ_ind)
            r = rhat(θ, Ψy)
            θ1[k] = θk[θ_ind[k]]
            for l in findall(!iszero, Hinv[k])
                θ1[k] += Hinv[k][l] * ( μx[l] - mean( r .* Ψy[l, :] ) )
            end
        end
        return θ1
    else
        for k in 1:length(θ_ind)
            supp = _find_supp(θ_ind[k], θ)
            θk = KLIEP(Ψx[supp, :], Ψy[supp, :], CD_KLIEP())

            r = rhat(θk, Ψy[supp, :])
            θ1[k] = θk[end]
            for l in findall(!iszero, Hinv[k])
                θ1[k] += Hinv[k][l] * ( μx[l] - mean( r .* Ψy[l, :] ) )
            end
        end
        return θ1
    end
end

function _debias2(Ψx, Ψy, θ, ω, θ_ind::Int)
    supp = _find_supp(θ_ind, θ, ω)
    θk = KLIEP(Ψx[supp,:], Ψy[supp,:], CD_KLIEP())
    θk[end]
end

function _debias2(Ψx, Ψy, θ, Hinv, θ_ind::Union{Vector{Int},UnitRange})
    θ2 = Vector{Float64}(undef, length(θ_ind))
    for k in 1:length(θ_ind)
        supp = _find_supp(θ_ind[k], θ, Hinv[k])
        θk = KLIEP(Ψx[supp, :], Ψy[supp, :], CD_KLIEP())
        θ2[k] = θk[end]
    end
    θ2
end


text:
"""
$(TYPEDSIGNATURES)

A heuristic estimator for the mu and sigma values of a GMM.
"""
function adaptive_parzen_normal(
    mus::Vector{Float64}, prior_mu::Float64, prior_sigma::Float64, config::Config
)::NTuple{3, Vector{Float64}}

    if prior_sigma <= 0
        throw(DimensionMismatch("prior_sigma: $(prior_sigma) is less than or equal to 0"))
    end

    srtd_mus = []
    sigma = []
    prior_pos = 1
    # sortperm must be used here because order is applied to unsorted_weights below
    order = sortperm(mus)

    if isempty(mus)
        srtd_mus = [prior_mu]
        sigma = [prior_sigma]
    elseif length(mus) == 1
        if prior_mu < first(mus)
            srtd_mus = [prior_mu, first(mus)]
            sigma = [prior_sigma, prior_sigma * 0.5]
        else
            prior_pos = 2
            srtd_mus = [first(mus), prior_mu]
            sigma = [prior_sigma * 0.5, prior_sigma]
        end
    elseif length(mus) >= 2
        # Create new_mus, which is sorted, and in which the prior has been inserted
        prior_pos = searchsortedfirst(mus[order], prior_mu)
        srtd_mus = mus[order]
        splice!(srtd_mus, prior_pos:(prior_pos - 1), prior_mu)
        sigma = zero(srtd_mus)
        sigma[2:end-1] = max.(
            srtd_mus[2:end-1] - srtd_mus[1:end-2],
            srtd_mus[3:end] - srtd_mus[2:end-1]
        )
        lsigma = srtd_mus[2] - srtd_mus[1]
        usigma = srtd_mus[end] - srtd_mus[end-1]
        sigma[1] = lsigma
        sigma[end] = usigma
    end

    if config.linear_forgetting < length(mus)
        unsorted_weights = ForgettingWeights.forgetting_weights(
            length(mus), config.linear_forgetting
        )
        if length(unsorted_weights) + 1 != length(srtd_mus)
            throw(DimensionMismatch(string(
                "length(unsorted_weights) + 1: ", length(unsorted_weights) + 1,
                " doesn't equal length(srtd_mus): ", length(srtd_mus)
            )))
        end
        sorted_weights = unsorted_weights[order]
        splice!(sorted_weights, prior_pos:(prior_pos - 1), config.prior_weight)
    else
        sorted_weights = ones(length(srtd_mus))
        sorted_weights[prior_pos] = config.prior_weight
    end

    # Magic formula
    maxsigma = prior_sigma
    minsigma = prior_sigma / min(100.0, (1.0 + length(srtd_mus)))
    clamp!(sigma, minsigma, maxsigma)
    sigma[prior_pos] = prior_sigma

    if all(sigma .<= 0)
        throw(DimensionMismatch(string(
            "minimum(sigma): ", minimum(sigma),
            " minsigma: ", minsigma,
            " maxsigma: ", maxsigma
        )))
    end

    sorted_weights /= sum(sorted_weights)

    return sorted_weights, srtd_mus, sigma
end


text:
using MonteCarlo

################################################################################
### Occupation + Pairing Correlation
################################################################################

run_idx = 1

begin
    if run_idx == 1
        betas = (2.0, 5.0, 7.0)
        lattice = TriangularLattice(4)
        mus = vcat(-2.0, -1.5, -1.25:0.05:-1.0, -0.8:0.2:0.8, 0.9:0.05:1.25, 1.5, 2.0)
    elseif run_idx == 2
        betas = (2.0, 5.0, 7.0)
        lattice = TriangularLattice(6)
        mus = vcat(-2.0:0.25:-0.25, -0.1:0.1:1.1, 1.25, 1.5, 2.0)
    elseif run_idx == 3
        betas = (2.0, 5.0, 6.0, 7.0)
        lattice = TriangularLattice(8)
        mus = vcat(-2.0:0.5:-0.5, -0.1:0.1:1.1, 1.25, 1.5, 2.0)
    else
        error("Invalid run_idx = $run_idx")
    end

    dqmcs = []
    counter = 0
    N = length(mus) * length(betas)
    @time for beta in betas, mu in mus
        counter += 1
        print("\r[", lpad("$counter", 2), "/$N]")
        m = HubbardModelAttractive(l = lattice, t = 1.0, U = 4.0, mu = mu)
        dqmc = DQMC(
            m, beta = beta, delta_tau = 0.125, safe_mult = 8,
            thermalization = 1000, sweeps = 1000, measure_rate = 1,
            recorder = Discarder()
        )
        dqmc[:occ] = occupation(dqmc, m)
        dqmc[:PC] = pairing_correlation(dqmc, m, kernel = MonteCarlo.pc_kernel)
        run!(dqmc, verbose = false)

        # for simplicity we just keep the whole simulation around
        push!(dqmcs, dqmc)
    end

    N = length(mus)
    occs = []
    Δoccs = []
    pcs = []
    Δpcs = []

    for i in 0:length(betas)-1
        # Measurements are saved in a LogBinner from BinningAnalysis by default.
        # Taking the mean (std_error) of a LogBinner will return the Monte Carlo
        # average (error). Occupation measurements happen per site, so we need
        # another mean afterwards.
        _occs = [2 * mean(mean(dqmcs[N*i + j][:occ])) for j in 1:N]
        doccs = [2 * mean(std_error(dqmcs[N*i + j][:occ])) for j in 1:N]
        push!(occs, _occs)
        push!(Δoccs, doccs)

        # pairing correlations are saved in a partially processed state - a 3D matrix
        # where each index corresponds to vectors between sites
        # y_{i, j, k} = ∑_x ⟨c_{x, ↑} c_{x+j, ↓} c_{x+i+k, ↓}^† c_{x+i, ↑}^†
        # The vectors corresponding to the indices i, j, k are returned by
        # directions(lattice(dqmc)). To compute the pairing correlation of a certain
        # symmetry, we need to apply the weights corresponding vector indices j, k.
        # For s-wave symmetry these weights are always (1, 0, ..., 0) (only vector 0).
        # To match the paper the index i should just be summed over. This is
        # equivalent to a q=0 Fourier transform.
        _pcs = [sum(mean(dqmcs[N*i + j][:PC])[:, 1, 1]) for j in 1:N]
        dpcs = [sum(std_error(dqmcs[N*i + j][:PC])[:, 1, 1]) for j in 1:N]
        push!(pcs, _pcs)
        push!(Δpcs, dpcs)
    end
end

using CairoMakie, FileIO, Colors

begin

    fig = Figure(resolution = (800, 800))
    top = Axis(fig[1, 1])
    bot = Axis(fig[2, 1])

    # References
    p = joinpath(pkgdir(MonteCarlo), "docs/src/examples/assets/triangular")
    if run_idx == 1
        top_ref = FileIO.load(joinpath(p, "tri_Hub_ref1_1.png"))
        bot_ref = FileIO.load(joinpath(p, "tri_Hub_ref1_2.png"))
    elseif run_idx == 2
        top_ref = FileIO.load(joinpath(p, "tri_Hub_ref2_1.png"))
        bot_ref = FileIO.load(joinpath(p, "tri_Hub_ref2_2.png"))
    elseif run_idx == 3
        top_ref = FileIO.load(joinpath(p, "tri_Hub_ref3_1.png"))
        bot_ref = FileIO.load(joinpath(p, "tri_Hub_ref3_2.png"))
    end
    ip = image!(top, -2..2, 0..2, top_ref'[:, end:-1:1])
    translate!(ip, 0, 0, -1)
    if run_idx == 1
        ip = image!(bot, -2..2, 0..2.5, bot_ref'[:, end:-1:1], transparency=true)
    elseif run_idx == 2
        ip = image!(bot, -2..2, 0..3.5, bot_ref'[:, end:-1:1], transparency=true)
    elseif run_idx == 3
        ip = image!(bot, -2..2, 0..6.3, bot_ref'[:, end:-1:1], transparency=true)
    end
    translate!(ip, 0, 0, -1)

    c = HSV(250, 0.6, 1)
    markers = length(betas) == 3 ? ('■', '□', 'o') : ('■', '□', '△', 'o')
    for (i, (ys, dys)) in enumerate(zip(occs, Δoccs))
        band!(top, mus, ys .- dys, ys .+ dys, color = (:red, 0.3))
        lines!(top, mus, ys, color = (c, 0.5), linewidth=2)
        scatter!(top, mus, ys, color = c, marker = markers[i])
    end

    labels = map(beta -> "β = $beta", collect(betas))
    axislegend(top, top.scene.plots[5:3:end], labels, position = :rb)

    for (i, (ys, dys)) in enumerate(zip(pcs, Δpcs))
        band!(bot, mus, ys .- dys, ys .+ dys, color = (:red, 0.3), transparency=true)
        lines!(bot, mus, ys, color = (c, 0.5), linewidth=2)
        scatter!(bot, mus, ys, color = c, marker = markers[i])
    end

    xlims!(top, -2 , 2)
    ylims!(top, 0 , 2)
    xlims!(bot, -2 , 2)
    run_idx == 1 && ylims!(bot, 0 , 2.5)
    run_idx == 2 && ylims!(bot, 0 , 3.5)
    run_idx == 3 && ylims!(bot, 0 , 6.3)

    display(fig)
    run_idx == 1 && CairoMakie.save(joinpath(p, "fig1_comparison.png"), fig)
    run_idx == 2 && CairoMakie.save(joinpath(p, "fig2_comparison.png"), fig)
    run_idx == 3 && CairoMakie.save(joinpath(p, "fig3_comparison.png"), fig)
end

################################################################################
### Charge Density
################################################################################

using MonteCarlo, LinearAlgebra

Ls = (6, 8)
betas = [8.0]
dqmcs = []

function my_kernel(mc, m, ij::NTuple{2}, G::AbstractArray)
    i, j = ij
    4 * (I[j, i] - G[j, i]) * G[i, j]
end

counter = 0
@time for beta in betas, L in Ls
    counter += 1
    print("\r[", lpad("$counter", 2), "/$(length(betas) * length(Ls))]")
    lattice = TriangularLattice(L)

    m = HubbardModelAttractive(l = lattice, t = 1.0, U = 4.0, mu = 0)
    dqmc = DQMC(
        m, beta = beta, delta_tau = 0.125, safe_mult = 8,
        thermalization = 1000, sweeps = 1000, measure_rate = 1,
        recorder = Discarder()
    )
    dqmc[:CDC] = MonteCarlo.Measurement(dqmc, m, Greens, EachSitePairByDistance(), my_kernel)
    run!(dqmc, verbose = false)

    # for simplicity we just keep the whole simulation around
    push!(dqmcs, dqmc)
end

qs = vcat(
    range(Float64[0, 0],   Float64[pi, 0], length=10),
    range(Float64[pi, 0],  Float64[pi, pi], length=10),
    range(Float64[pi, pi], Float64[0, 0], length=10),
)
ys = map(dqmcs) do dqmc
    MonteCarlo.fourier_transform(
        qs, directions(MonteCarlo.lattice(dqmc)), mean(dqmc[:CDC])
    ) |> real
end

using CairoMakie, FileIO, Colors

fig = Figure(resolution = (800, 800))
ax = Axis(fig[1, 1])

# References
p = joinpath(pkgdir(MonteCarlo), "docs/src/examples/assets/triangular")
ref = FileIO.load(joinpath(p, "tri_Hubbard_CDSF.png"))
ip = image!(ax, 0..3, 0..1.5, ref'[:, end:-1:1])
translate!(ip, 0, 0, -1)

c = HSV(250, 0.6, 1)
for (i, ys) in enumerate(ys)
    xs = vcat(range(0, 1, length=10), range(1, 2, length=10), range(2, 3, length=10))
    Makie.lines!(ax, xs, ys, color = (c, 0.5), linewidth=2)
    scatter!(ax, xs, ys, color = c, marker = ('■', 'o', 'x', '+', '*')[i])
end

ax.xticks[] = ([0,1,2,3], ["(0,0)", "(π, 0)", "(π,π)", "(0,0)"])
xlims!(ax, 0, 3)
ylims!(ax, 0 , 1.5)

fig
CairoMakie.save(joinpath(p, "fig6_comparison.png"), fig)

################################################################################
### Spin Susceptibility
################################################################################

using MonteCarlo, LinearAlgebra

Ls = (6,)
betas = [1.0, 2.0, 4.0, 5.0, 6.0, 7.0, 8.0]
dqmcs = []

counter = 0
@time for beta in betas, L in Ls
    counter += 1
    print("\r[", lpad("$counter", 2), "/$(length(betas))]")
    lattice = TriangularLattice(L)

    m = HubbardModelAttractive(l = lattice, t = 1.0, U = 4.0, mu = 0)
    dqmc = DQMC(
        m, beta = beta, delta_tau = 0.125, safe_mult = 8,
        thermalization = 1000, sweeps = 1000, measure_rate = 1,
        recorder = Discarder()
    )
    dqmc[:SDS] = spin_density_susceptibility(dqmc, m, :z)
    run!(dqmc, verbose = false)

    # for simplicity we just keep the whole simulation around
    push!(dqmcs, dqmc)
end

xs = map(dqmc -> 1.0 / dqmc.parameters.beta, dqmcs)
ys = map(dqmc -> (dqmc[:SDS] |> mean |> sum), dqmcs)
dys = map(dqmc -> sqrt((dqmc[:SDS] |> std_error).^2 |> sum), dqmcs)

using CairoMakie, FileIO, Colors

fig = Figure(resolution = (800, 800))
ax = Axis(fig[1, 1])

# References
p = joinpath(pkgdir(MonteCarlo), "docs/src/examples/assets/triangular")
ref = FileIO.load(joinpath(p, "tri_Hubbard_SDS.png"))
ip = image!(ax, 0..1, 0..0.20, ref'[:, end:-1:1])
translate!(ip, 0, 0, -1)

c = HSV(250, 0.6, 1)
band!(ax, xs, ys .- dys, ys .+ dys, color = (:red, 0.3))
Makie.lines!(ax, xs, ys, color = (c, 0.5), linewidth=2)
scatter!(ax, xs, ys, color = c, marker = ('■', 'o', 'x', '+', '*')[1])

xlims!(ax, 0, 1)
ylims!(ax, 0 , 0.20)

fig
CairoMakie.save(joinpath(p, "fig7_comparison.png"), fig)


text:
using DifferentialEquations
using PyCall
using PyPlot

function tusharabmodels()
    # Quick plots of Tushar's a and b
    v = collect(LinRange(1e-15, 1e1, 100000))
    amin = 0.045
    ahighdiff =  0.0319 - amin
    δ = 0.7 # m/s
    aplot = @. amin + ahighdiff * 0.5 * (1 + tanh((log10(v) - log10(5e-3)) / δ))
    bmin = 0.085
    bhighdiff =  0.2811 - bmin
    ##  0.5 * (1 + tanh((log10(v)-np.log10(5e-3)) / 7e-1))
    bplot = @. bmin + bhighdiff * 0.5 * (1 + tanh((log10(v) - log10(5e-3)) / δ))
    plot(v, aplot, "-r", label="a")
    plot(v, bplot, "-b", label="b")
    xscale("log")
    xlabel(L"\log_{10} v")
    ylabel(L"a, b")
    legend()
end

function aginglaw(v, θ, dc)
    return 1 - θ * v / dc
end

function sliplaw(v, θ, dc)
    return -v * θ / dc * log(v * θ / dc)
end

function calcdadv(v)
    amin = 0.045
    ahighdiff =  0.0319 - amin
    δ = 0.7 # m/s
    # a = amin + ahighdiff * 0.5 * (1 + tanh((log10(v) - log10(5e-3)) / δ))
    dadv = ahighdiff * 0.5 * (sech(( -log10(v) + log10(5e-3)) / (δ * log(10))  ))^2 / (δ * v * log(10))
    return dadv
end

function calcdbdv(v)
    bmin = 0.085
    bhighdiff =  0.2811 - bmin
    δ = 0.7 # m/s
    #  bmin + bhighdiff * 0.5 * (1 + tanh((log10(v) - log10(5e-3)) / δ))
    dbdv = bhighdiff * 0.5 * (sech(( -log10(v) + log10(5e-3)) / (δ * log(10))  ))^2 / (δ * v * log(10))
    return dbdv
end

function calcdvθclassic!(du, u, p, t)
    dc, η, σn, a, b, μ, Vp, L, ρ, statelaw = p
    θ = u[1]
    v = u[2]
    du[1] = statelaw(v, θ, dc)
    du[2] = 1 / (η / σn + a / v) * (μ * (Vp - v) / (L * σn) - b * du[1] / θ)
    return nothing
end

function calcdvθab!(du, u, p, t)
    dc, η, σn, μ, vp, L, ρ, statelaw, fstar, vstar, θstar = p
    θ = u[1]
    v = u[2]
    a = u[3]
    b = u[4]
    dθ = du[1]
    dv = du[2]
    du[1] = statelaw(v, θ, dc)
    da = calcdadv(v) * dv
    db = calcdbdv(v) * dv
    numclassic = (μ * (vp - v) / (L * σn) - b * du[1] / θ)
    numab = -da * log(abs(v) / vstar) - db * log(abs(θ) / θstar)
    numσn = 0
    denom = (η / σn + a / v)
    du[2] = (numclassic + numab + numσn) / denom
    du[3] = da
    du[4] = db
    return nothing
end

function plottimeseries(sol, siay, titlelabel)
    t1 = [x / siay for x in sol.t]
    θ1 = [x[1] for x in sol.u]
    v1 = [x[2] for x in sol.u]

    figure(figsize = (12, 6))
    subplot(2, 2, 1)
    plot(t1, θ1, "-b", linewidth = 0.5)
    yscale("log")
    ylabel(L"\theta")

    subplot(2, 2, 2)
    plot(1:1:length(t1), θ1, "-b", linewidth = 0.5)
    yscale("log")
    ylabel(L"\theta")

    subplot(2, 2, 3)
    plot(t1, v1, "-b", linewidth = 0.5)
    yscale("log")
    xlabel("time (years)")
    ylabel("v (m/s)")

    subplot(2, 2, 4)
    plot(1:1:length(t1), v1, "-b", linewidth = 0.5)
    yscale("log")
    xlabel("time (step)")
    ylabel("v (m/s)")
    suptitle(titlelabel)

    return nothing
end

function sliding()
    # Model parameters
    siay = 365.25 * 24 * 60 * 60
    tspan = (0.0, siay * 330.0) # This will run
    # tspan = (0.0, siay * 302.0) # This will run but fail for the a, b case
    μ = 3e10
    ν = 0.25
    ρ = 2700.0
    η = μ / (2.0 * sqrt(μ / ρ))
    L = 60 * 1e3
    a = 0.015
    b = 0.02
    vp = 1e-9
    σn = 30e6
    dc = 0.1
    abstol = 1e-4
    reltol = 1e-4

    # Parameters need for d(a, b σn)/dt
    fstar = 0.6 # ???
    vstar = vp # ???
    θstar = 1e9 # ???

    # Show Tushar scalings for a and b
    tusharabmodels()

    # Time integrate - classic
    icsclassic = [1e8; vp / 1000]
    pclassic = (dc, η, σn, a, b, μ, vp, L, ρ, aginglaw)
    probclassic = ODEProblem(calcdvθclassic!, icsclassic, tspan, pclassic)
    solclassic = solve(probclassic, RK4(), abstol = abstol, reltol = reltol)
    plottimeseries(solclassic, siay, "RSF (classic)")

    # Time integrate - a, b can evolve
    icsab = [1e8 ; vp / 1000 ; a ; b]
    pab = (dc, η, σn, μ, vp, L, ρ, aginglaw, fstar, vstar, θstar)
    probab = ODEProblem(calcdvθab!, icsab, tspan, pab)
    solab = solve(probab, RK4(), abstol = abstol, reltol = reltol)
    plottimeseries(solab, siay, "RSF (ab)")

    return nothing
end
sliding()


text:
#=
day18:
- Julia version: 1.7.0
- Author: Paul.Mealor
- Date: 2021-12-18
=#

using AoC, Test

struct SfNumber
    left::Union{Int,SfNumber}
    right::Union{Int,SfNumber}
end
Base.show(io::IO, n::SfNumber) = print(io, "[",n.left,",",n.right,"]")

function sfparse(line)
    if !isnothing(match(r"^\d+$", line))
        return parse(Int, line)
    end
    depth = 0
    for i in 2:length(line)-1
        c = line[i]
        if c == '['
            depth += 1
        elseif c == ']'
            depth -= 1
        elseif c == ','
            if depth == 0
                return SfNumber(sfparse(line[2:i-1]), sfparse(line[i+1:length(line)-1]))
            end
        end
    end
    throw(ArgumentError(line))
end

@test sfparse("[[1,2],[[3,4],5]]") == SfNumber(SfNumber(1,2),SfNumber(SfNumber(3,4),5))
@test sfparse("[[1,2],[[15,4],5]]") == SfNumber(SfNumber(1,2),SfNumber(SfNumber(15,4),5))

explode(n::Int, depth=1) = (0,0,n)
function explode(n::SfNumber, depth=1)
    if depth > 4
        return (n.left, n.right, 0)
    end

    (leftFragment, rightFragment, left) = explode(n.left, depth+1)
    if left != n.left
        right = applyExplodeRight(n.right, rightFragment)
        return (leftFragment, 0, SfNumber(left, right))
    end

    (leftFragment, rightFragment, right) = explode(n.right, depth+1)
    if right != n.right
        left = applyExplodeLeft(n.left, leftFragment)
        return (0, rightFragment, SfNumber(left, right))
    end

    return (0,0,n)
end

applyExplodeRight(n::Int, fragment) = n+fragment
applyExplodeRight(n::SfNumber, fragment) = SfNumber(applyExplodeRight(n.left, fragment), n.right)

applyExplodeLeft(n::Int, fragment) = n+fragment
applyExplodeLeft(n::SfNumber, fragment) = SfNumber(n.left, applyExplodeLeft(n.right, fragment))

@test explode(sfparse("[[[[[9,8],1],2],3],4]"))[3] == sfparse("[[[[0,9],2],3],4]")
@test explode(sfparse("[7,[6,[5,[4,[3,2]]]]]"))[3] == sfparse("[7,[6,[5,[7,0]]]]")
@test explode(sfparse("[[6,[5,[4,[3,2]]]],1]"))[3] == sfparse("[[6,[5,[7,0]]],3]")
@test explode(sfparse("[[3,[2,[1,[7,3]]]],[6,[5,[4,[3,2]]]]]"))[3] == sfparse("[[3,[2,[8,0]]],[9,[5,[4,[3,2]]]]]")
@test explode(sfparse("[[3,[2,[8,0]]],[9,[5,[4,[3,2]]]]]"))[3] == sfparse("[[3,[2,[8,0]]],[9,[5,[7,0]]]]")

split(n::Int) = n > 9 ? SfNumber(n ÷ 2, (n+1) ÷ 2) : n
function split(n::SfNumber)
    left = split(n.left)
    if left != n.left
        return SfNumber(left, n.right)
    end

    right = split(n.right)
    if right != n.right
        return SfNumber(n.left, right)
    end

    return n
end

@test split(9) == 9
@test split(10) == SfNumber(5,5)
@test split(11) == SfNumber(5,6)
@test split(12) == SfNumber(6,6)

function sfreduce(n)
    while true
        (_,_,newValue) = explode(n)
        if newValue == n
            newValue = split(n)
            if newValue == n
                return n
            end
        end
        n = newValue
    end
end

@test sfreduce(sfparse("[[[[[4,3],4],4],[7,[[8,4],9]]],[1,1]]")) == sfparse("[[[[0,7],4],[[7,8],[6,0]]],[8,1]]")

sfsum(l,r) = sfreduce(SfNumber(l,r))
sfsum(args...) = sfsum(sfsum(args[1], args[2]), args[3:end]...)
@test sfsum(sfparse("[[[0,[4,5]],[0,0]],[[[4,5],[2,6]],[9,5]]]"), sfparse("[7,[[[3,7],[4,3]],[[6,3],[8,8]]]]")) == sfparse("[[[[4,0],[5,4]],[[7,7],[6,0]]],[[8,[7,7]],[[7,9],[5,0]]]]")
@test sfsum(sfparse("[[[[7,7],[7,7]],[[8,7],[8,7]]],[[[7,0],[7,7]],9]]"), sfparse("[[[[4,2],2],6],[8,7]]")) == sfparse("[[[[8,7],[7,7]],[[8,6],[7,7]]],[[[0,7],[6,6]],[8,7]]]")

magnitude(n::Int) = n
magnitude(n::SfNumber) = 3*magnitude(n.left) + 2*magnitude(n.right)
@test magnitude(sfparse("[[1,2],[[3,4],5]]")) == 143
@test magnitude(sfparse("[[[[0,7],4],[[7,8],[6,0]]],[8,1]]")) == 1384
@test magnitude(sfparse("[[[[1,1],[2,2]],[3,3]],[4,4]]")) == 445
@test magnitude(sfparse("[[[[3,0],[5,3]],[4,4]],[5,5]]")) == 791
@test magnitude(sfparse("[[[[5,0],[7,4]],[5,5]],[6,6]]")) == 1137
@test magnitude(sfparse("[[[[8,7],[7,7]],[[8,6],[7,7]]],[[[0,7],[6,6]],[8,7]]]")) == 3488

part1(lines) = magnitude(sfsum(sfparse.(lines)...))
@test part1(exampleLines(18,1)) == 4140

function part2(lines)
    lines = sfparse.(lines)
    biggest = 0
    for i in eachindex(lines)
        for j in eachindex(lines)
            if i != j
                value = magnitude(sfsum(lines[i],lines[j]))
                if value > biggest
                    biggest = value
                end
            end
        end
    end
    return biggest
end

@test part2(exampleLines(18,1)) == 3993

@show lines(18) |> ll -> @time part1(ll)
@show lines(18) |> ll -> @time part2(ll)


text:
export sets

const sets = [:GRN, :M19, :DOM, :RIX, :XLN, :HOU, :AKH, :AER, :KLD, :EMN, :SOI, :OGW, :BFZ]
const set_counter = Dict( (s,0.0) for s in sets )

max_stars_repo_name: asieev/ArenaSim

max_issues_repo_name: asieev/ArenaSim

max_forks_repo_name: asieev/ArenaSim


text:
struct MyFunctor
    a::Float64
end
(mf::MyFunctor)(b) = mf.a + b

@testset "ignore_derivatives.jl" begin
    @testset "function" begin
        f() = return 4.0

        y, ẏ = frule((1.0,), ignore_derivatives, f)
        @test y == f()
        @test ẏ == NoTangent()

        y, pb = rrule(ignore_derivatives, f)
        @test y == f()
        @test pb(1.0) == (NoTangent(), NoTangent())
    end

    @testset "argument" begin
        arg = 2.1

        y, ẏ = frule((1.0,), ignore_derivatives, arg)
        @test y == arg
        @test ẏ == NoTangent()

        y, pb = rrule(ignore_derivatives, arg)
        @test y == arg
        @test pb(1.0) == (NoTangent(), NoTangent())
    end

    @testset "functor" begin
        mf = MyFunctor(1.0)

        # as an argument
        y, ẏ = frule((1.0,), ignore_derivatives, mf)
        @test y == mf
        @test ẏ == NoTangent()

        y, pb = rrule(ignore_derivatives, mf)
        @test y == mf
        @test pb(1.0) == (NoTangent(), NoTangent())

        # when called
        y, ẏ = frule((1.0,), ignore_derivatives, () -> mf(3.0))
        @test y == mf(3.0)
        @test ẏ == NoTangent()

        y, pb = rrule(ignore_derivatives, () -> mf(3.0))
        @test y == mf(3.0)
        @test pb(1.0) == (NoTangent(), NoTangent())
    end
end


text:
"""
    ContactConstraint{T} <: Constraint{T}

    constraint containing information for contact node.

    id: unique identifying number
    name: unique identifying name
    model: type of contact model: ImpactContact, LinearContact, NonlinearContact
    parent_id: identifying number of Body experiencing contact
    child_id: always 0
    impulses: contact impulses applied to Body
    impulses_dual: dual contact impulses, used by solver to enforce correct contact behaviors
"""
mutable struct ContactConstraint{T,N,Nc,Cs,N½} <: Constraint{T,N}
    # ID
    id::Int64
    name::Symbol

    # contact model
    model::Cs

    # neighbor IDs
    parent_id::Int
    child_id::Int

    # variables
    impulses::Vector{SVector{N½,T}}
    impulses_dual::Vector{SVector{N½,T}}

    function ContactConstraint(data;
        name::Symbol=Symbol("contact_" * randstring(4)))

        model, parent_id, child_id = data
        T = typeof(model).parameters[1]

        N = length(model)
        N½ = Int64(N/2)

        impulses = [neutral_vector(model) for i = 1:2]
        impulses_dual = [neutral_vector(model) for i = 1:2]
        new{T,N,1,typeof(model),N½}(getGlobalID(), name, model, parent_id, child_id, impulses, impulses_dual)
    end
end

# function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, constraint::ContactConstraint)
#     summary(io, constraint)
#     println(io, "")
#     println(io, "id:            "*string(constraint.id))
#     println(io, "name:          "*string(constraint.name))
#     println(io, "model:         "*string(constraint.model))
#     println(io, "parent_id:     "*string(constraint.parent_id))
#     println(io, "child_id:      "*string(constraint.child_id))
#     println(io, "impulses:      "*string(constraint.impulses))
#     println(io, "impulses_dual: "*string(constraint.impulses_dual))
# end

"""
    contact_constraint(bodies::Vector{Body})

    generate ContactConstraints for each Body in list

    normal: surface normal for each contact point
    friction coefficient: value of coefficient of friction for each contact point (optional for ImpactContact)
    contact_origins: the offset with respect to the center of Body for each contact point (optional)
    contact_radius: radius for each contact (optional)
    contact_type: :nonlinear, :linear, :impact
"""
function contact_constraint(bodies::Vector{Body{T}},
        normal::AbstractVector{<:AbstractVector};
        friction_coefficient::AbstractVector{T}=ones(length(normal)),
        contact_origins::AbstractVector=[szeros(T, 3) for i=1:length(normal)],
        contact_radius::AbstractVector=[0.0 for i=1:length(normal)],
        names::Vector{Symbol}=[Symbol("contact_" * randstring(4)) for i = 1:length(normal)],
        contact_type::Symbol=:nonlinear) where T

    n = length(normal)
    @assert n == length(bodies) == length(normal) == length(friction_coefficient) == length(contact_origins) == length(contact_radius)
    contacts = Vector{ContactConstraint}()
    for i = 1:n
        contact = contact_constraint(bodies[i], normal[i],
            friction_coefficient=friction_coefficient[i],
            contact_origin=contact_origins[i],
            contact_radius=contact_radius[i],
            name=names[i],
            contact_type=contact_type)
        push!(contacts, contact)
    end
    contacts = [contacts...] # vector typing
    return contacts
end

function contact_constraint(body::Body{T},
        normal::AbstractVector{<:AbstractVector};
        friction_coefficient::AbstractVector{T}=ones(length(normal)),
        contact_origins::AbstractVector=[szeros(T, 3) for i=1:length(normal)],
        contact_radius::AbstractVector=[0.0 for i=1:length(normal)],
        names::Vector{Symbol}=[Symbol("contact_" * randstring(4)) for i = 1:length(normal)],
        contact_type::Symbol=:nonlinear) where T
    n = length(normal)
    @assert n == length(normal) == length(friction_coefficient) == length(contact_origins) == length(contact_radius)
    return contact_constraint(fill(body, n), normal,
        friction_coefficient=friction_coefficient,
        contact_origins=contact_origins,
        contact_radius=contact_radius,
        names=names,
        contact_type=contact_type)
end

function contact_constraint(body::Body{T},
        normal::AbstractVector{T};
        friction_coefficient::T=1.0,
        contact_origin::AbstractVector{T}=szeros(T, 3),
        contact_radius::T=0.0,
        name::Symbol=Symbol("contact_" * randstring(4)),
        contact_type::Symbol=:nonlinear) where T

    if contact_type == :nonlinear
        model = NonlinearContact(body, normal, friction_coefficient,
            contact_origin=contact_origin,
            contact_radius=contact_radius)
    elseif contact_type == :linear
        model = LinearContact(body, normal, friction_coefficient,
            contact_origin=contact_origin,
            contact_radius=contact_radius)
    elseif contact_type == :impact
        model = ImpactContact(body, normal,
            contact_origin=contact_origin,
            contact_radius=contact_radius)
    else
        @warn "unknown contact_type"
    end
    contacts = ContactConstraint((model, body.id, 0);
        name=name)
    return contacts
end


text:
using Thetafunc

println("On teste le package Theta")

Thetafunc_success=true
r=0.1

#theta(x::Vector,r::Float64)
#Les fonctions thetas doivent etre nulle en 0, <=1 pour x>=0 et <0 pour x<0
if Thetafunc.theta(0.0,r)!=0 && Thetafunc.theta(1.0+1e6,r)>=1 && Thetafunc(-1e-6,r)>=0
 println("Error in the definition of theta")
 plot(x->Thetafunc.theta(x,r),-0.1,2.0)
 #draw(SVG("output.svg",6inch,3inch),plot(x->Thetamod.theta(x,r),-0.1,2.0))
 Thetafunc_success=false
end

#dtheta(x::Vector,r::Float64)
#la dérivée des fonctions theta doit être croissante
if Thetafunc.theta(1.0+1e6,r)>=0 && (Thetafunc.theta(1.0,r)-Thetafunc.theta(2.0,r)>0)
 println("Error in the definition of dtheta")
 plot(x->Thetafunc.dtheta(x,r),-0.1,2.0)
 Thetafunc_success=false
end

#invtheta(x::Vector,r::Float64)
#Tester aussi les cas interdits ici pas défini pour x>1 !!!
if norm(Thetafunc.invtheta(Thetafunc.theta(1e-6,r),r)-1e-6)>precision || norm(Thetafunc.invtheta(Thetafunc.theta(-1e-6,r),r)+1e-6)>precision
 println("Error in the definition of invtheta")
 println(norm(Thetafunc.invtheta(Thetafunc.theta(1e-6,r),r)-1e-6))
 println(norm(Thetafunc.invtheta(Thetafunc.theta(-1e-6,r),r)+1e-6))
 plot(x->Thetafunc.invtheta(Thetafunc.theta(x,r),r),-0.1,1.0)
 Thetafunc_success=false
end

#AlphaThetaMax(x::Float64,dx::Float64,y::Float64,dy::Float64,r::Float64,s::Float64,t::Float64)
#AlphaThetaMax : résoud l'équation en alpha où psi(x,r,s,t)=s+t*theta(x,r)
r=1.0;s=1.0;t=0.0 #résoud l'éq.

x=1.0;y=2.0;dx=1.0;dy=-1.0;
#doit rendre (0,1)
if !isequal(Thetafunc.AlphaThetaMax(x,dx,y,dy,r,s,t),[0.0,1.0])
 println(Thetafunc.AlphaThetaMax(x,dx,y,dy,r,s,t))
 Thetafunc_success=false
end

x=0.0;y=2.0;dx=1.0;dy=-2.0;
#solution : (1,0.5)
if !isequal(Thetafunc.AlphaThetaMax(x,dx,y,dy,r,s,t),[1.0,0.5])
 println(Thetafunc.AlphaThetaMax(x,dx,y,dy,r,s,t))
 Thetafunc_success=false
end

r=1.0;s=1.0;t=1.0
#on tape la contrainte :
x=1.0;y=2.0;dx=1.0;dy=0.0;
if !(minimum(Thetafunc.AlphaThetaMax(x,dx,y,dy,r,s,t))>=0.0) || !(maximum(Thetafunc.AlphaThetaMax(x,dx,y,dy,r,s,t))==Inf)
 println(Thetafunc.AlphaThetaMax(x,dx,y,dy,r,s,t))
 Thetafunc_success=false
end

x=2.0;y=1.0;dx=-1.0;dy=2.0;
if !(minimum(Thetafunc.AlphaThetaMax(x,dx,y,dy,r,s,t))>=0.0) || !(maximum(Thetafunc.AlphaThetaMax(x,dx,y,dy,r,s,t))<Inf)
 println(Thetafunc.AlphaThetaMax(x,dx,y,dy,r,s,t))
 Thetafunc_success=false
end

#on est sous (s,s) et donc on ne tape pas la contrainte virtuelle
x=0.0;y=0.0;dx=0.0;dy=-2.0;
if !isequal(maximum(Thetafunc.AlphaThetaMax(x,dx,y,dy,r,s,t)),Inf)
 println(Thetafunc.AlphaThetaMax(x,dx,y,dy,r,s,t))
 Thetafunc_success=false
end

#Bilan : donner une sortie bilan
if Thetafunc_success==true
 println("Theta.jl passes the test !")
else
 println("Theta.jl contains some errors")
end


text:
module TargetFunction

immutable ObjFun
    name::ASCIIString
    f_x::Function
    g_x::Function
    about_alpha::Function
    about_alpha_prime::Function
end

function f_x(x::Array{Float64, 1})
    f = 10.0 * (x[2] - x[1]*x[1])*(x[2] - x[1]*x[1]) + (1 - x[1])*(1 - x[1])
    return f
end

function g_x(x::Array{Float64, 1})
    g1 = -40.0 * (x[1]*x[2] - x[1]*x[1]*x[1]) - 2.0 + 2.0 * x[1]
    g2 = 20.0 * (x[2] - x[1]*x[1])
    return [g1; g2]
end

function about_alpha(x::Array{Float64, 1}, s::Array{Float64, 1})
    function along_s(alpha::Float64)
        return f_x(x + s*alpha)
    end
    return along_s
end

function about_alpha_prime(x::Array{Float64, 1}, s::Array{Float64, 1})
    function along_s(alpha::Float64)
        return g_x(x + alpha * s)' * s
    end
    return along_s
end

rosenbrock = ObjFun("Rosenbrock", f_x, g_x, about_alpha, about_alpha_prime)

end

max_stars_repo_name: wyli/gradient

max_issues_repo_name: wyli/gradient

max_forks_repo_name: wyli/gradient


text:
module KuramotoSivashinksy

using FFTW

export ksintegrateNaive, ksintegrateTuned

function ksintegrateNaive(u, Lx, dt, Nt, nplot);
    Nx = length(u)                  # number of gridpoints
    kx = vcat(0:Nx/2-1, 0, -Nx/2+1:-1)  # integer wavenumbers: exp(2*pi*i*kx*x/L)
    alpha = 2*pi*kx/Lx              # real wavenumbers:    exp(i*alpha*x)
    D = 1im*alpha                   # D = d/dx operator in Fourier space
    L = alpha.^2 - alpha.^4         # linear operator -D^2 - D^4 in Fourier space
    G = -0.5*D                      # -1/2 D operator in Fourier space
    Nplot = round(Int64, Nt/nplot)+1  # total number of saved time steps

    x = collect(0:Nx-1)*Lx/Nx
    t = collect(0:Nplot)*dt*nplot
    U = zeros(Nplot, Nx)

    # some convenience variables
    dt2  = dt/2
    dt32 = 3*dt/2
    A =  ones(Nx) + dt2*L
    B = (ones(Nx) - dt2*L).^(-1)

    Nn  = G.*fft(u.*u); # -1/2 d/dx(u^2) = -u u_x, collocation calculation
    Nn1 = Nn;

    U[1,:] = u; # save initial value u to matrix U
    np = 2;     # counter for saved data

    # transform data to spectral coeffs
    u  = fft(u);

    # timestepping loop
    for n = 0:Nt-1
        Nn1 = Nn;                       # shift N^{n-1} <- N^n
        Nn  = G.*fft(real(ifft(u)).^2); # compute N^n = -u u_x

        u = B .* (A .* u + dt32*Nn - dt2*Nn1); # CNAB2 formula

        if mod(n, nplot) == 0
            U[np,:] = real(ifft(u))
            np += 1
        end
    end
    U,x,t
end

function ksintegrateTuned(u, Lx, dt, Nt);
    u = (1+0im)*u                       # force u to be complex
    Nx = length(u)                      # number of gridpoints
    kx = vcat(0:Nx/2-1, 0:0, -Nx/2+1:-1)# integer wavenumbers: exp(2*pi*kx*x/L)
    alpha = 2*pi*kx/Lx                  # real wavenumbers:    exp(alpha*x)
    D = 1im*alpha                       # spectral D = d/dx operator
    L = alpha.^2 - alpha.^4             # spectral L = -D^2 - D^4 operator
    G = -0.5*D                          # spectral -1/2 D operator, to eval -u u_x = 1/2 d/dx u^2

    # convenience variables
    dt2  = dt/2
    dt32 = 3*dt/2
    A =  ones(Nx) + dt2*L
    B = (ones(Nx) - dt2*L).^(-1)

    # compute in-place FFTW plans
    FFT! = plan_fft!(u, flags=FFTW.ESTIMATE)
    IFFT! = plan_ifft!(u, flags=FFTW.ESTIMATE)

    # compute uf == Fourier coeffs of u and Nnf == Fourier coeffs of -u u_x
    # FFT!(u);
    Nn  = G.*fft(u.^2); # Nnf == -1/2 d/dx (u^2) = -u u_x, spectral
    Nn1 = copy(Nn);     # use Nnf1 = Nnf at first time step
    FFT!*u;

    # timestepping loop, many vector ops unrolled to eliminate temporary vectors
    for n = 0:Nt

        for i = 1:length(Nn)
            @inbounds Nn1[i] = Nn[i];
            @inbounds Nn[i] = u[i];
        end

        IFFT!*Nn; # in-place FFT

        for i = 1:length(Nn)
            @fastmath @inbounds Nn[i] = Nn[i]*Nn[i];
        end

        FFT!*Nn;

        for i = 1:length(Nn)
            @fastmath @inbounds Nn[i] = G[i]*Nn[i];
        end

        for i = 1:length(u)
            @fastmath @inbounds u[i] = B[i]* (A[i] * u[i] + dt32*Nn[i] - dt2*Nn1[i]);
        end

    end
    u = real(ifft(u))
end

end # module


text:
### A Pluto.jl notebook ###
# v0.18.0

using Markdown
using InteractiveUtils

# ╔═╡ aff9313b-2fd8-46ca-b399-47fa21b324c5
begin
    import Pkg
    Pkg.activate(Base.current_project())
    Pkg.instantiate()
end

# ╔═╡ bcca4c48-89be-11ec-3c89-3189fd524f08
md"""
# Template-matching in Julia using Pluto.jl

## Introduction

Julia is a general programming language suited for scientific computing applications. Moreover, Julia is a dynamic, high-level language: executing a program written in Julia does not require the user to compile it beforehand, as for C or Fortran, and its syntax is similar to MATLAB or Python. Nonetheless, its performance is close to ones of static languages.

Pluto is a Julia library for writing notebooks similar (but better in many respect) to Jupyter notebooks. In this notebook I will work out an example and very briefly introduce Julia and Pluto. The goal is to template match some signals from acoustic measurements of hydraulic fracture growth performed under triaxial confinement in laboratory (more informations [here](https://www.researchgate.net/profile/Seyyedmaalek-Momeni/publication/352780078_Combining_active_and_passive_acoustic_methods_to_image_hydraulic_fracture_growth_in_laboratory_experiments/links/60d87837a6fdccb745ea29dd/Combining-active-and-passive-acoustic-methods-to-image-hydraulic-fracture-growth-in-laboratory-experiments.pdf)).
"""

# ╔═╡ 11ba1299-05e6-46c1-8716-ac06d5fc46a5
md"""
## Julia and Pluto.jl

### Pluto.jl

A Pluto notebook is a sequence of cells. Each cell contain some Julia code, which can be executed. In Julia, every statement is an expression and hence evaluates to a value. Pluto converts that value to an HTML element and shows it. What you are reading is the output of a cell containing just a string, which returns the string itself [^1]. Pluto is able to parse the content of cells and detect dependencies between one cell and another. When the user evaluates a cell, for example changing the value of a variable, Pluto will re-evaluate all the cells that depend on that one and in the right order.

Notice that Pluto requires one statement per cell, hence if you want to put more than one statement in a single cell you need to wrap them between the keywords `begin` and `end`, which are used to chain statements. The resulting block will evaluate to the value of the last statement.

For example

```julia
begin
	x = 1
	y = x + 1
	"hello world"
	2 * y
end
```

evaluates to 4.

[^1]: Actually, the value that is rendered as the piece of text you are reading is not a string, but an object whose type is defined in the Markdown standard library. You can create these object prefixing the string literals with `md`. In this way, the content will be formatted using Markdown. Otherwise you can use plain strings.
"""

# ╔═╡ 6dca0719-eeb3-42d5-8585-f89a77e579b3
md"""
### Defining Functions in Julia

Julia does not have methods attached to objects and accessible using the dot notation like in Python or C++. The main way of organizing code in Julia is by writing functions, and eventually grouping them into modules.

Function are defined using the `function` keyword. For example, the following will multiply the first two arguments and add the product to the third.

```julia
function f(x, y, z)
	x * y + z
end
```

Th `return` keyword is implied in the last statement of a function definition, but can be supplyied if one wants to.

Notice that the first time a function is called, Julia compiles on the fly that function. This is called just in time compilation or JIT. Hence the execution time of a function call change drastically between the first and subsequent calls. This is something to keep in mind when writing Julia code. There are some subtleties related to JIT compilation and the types of the arguments of a function, but we will talk about that in the next section.
"""

# ╔═╡ 364a442a-17ef-455d-9608-aecd3126db59
md"""
### Arrays and other types in Julia

In theory, one might write Julia code without thinking or caring about types, like in Python. However, it is practically impossibile to write and debug efficient code without some knowledge of the topic. One of the most important things to know is that the compiler can produce optimized machine code only for type-stable functions, where a function is called type-stabile if the type of its output is determined only from the types of its arguments (and not from their values).

The most used type in scientific computing is probably the array, i.e. a data type that represent indexed objects of the same type and contiguous in memory. Julia has a rich type system that includes arrays. In particular, arrays are an example of parametric types in Julia, i.e. types that depends on other types. The arguments of a dependent type are surrounded by braces, so for example `Array{Float64, 1}` is the type of one-dimensional arrays of 64-bit floating points numbers. To produce an exemplar of this type, you have to call its constructor. The constructor is a function with the same name of the type. In the case of one-dimensional array (or vector) constructors, they take two arguments, the value with which filling the array and the number of elements, so for example `Array{Float64, 1}(1.0, 10)` will be a ten element vector filled with the number 1.0. Array indexing in Julia is similar to Python with one important exceptions, in Julia one starts to count from 1 instead of 0, as you would do in Fortran or Matlab. Also, the last element is indicated with the keyword `end` and you can apply a function that acts on scalars element-wise on arrays posfixing a dot to its name in the function call, as in the function definition below (really, the dot notation is used for broadcasting, a powerful concept that should be familiar to Numpy user).
"""

# ╔═╡ Cell order:
# ╠═aff9313b-2fd8-46ca-b399-47fa21b324c5
# ╟─bcca4c48-89be-11ec-3c89-3189fd524f08
# ╟─11ba1299-05e6-46c1-8716-ac06d5fc46a5
# ╟─6dca0719-eeb3-42d5-8585-f89a77e579b3
# ╟─364a442a-17ef-455d-9608-aecd3126db59


text:
__precompile__()

module DeepDream

using Flux, Metalhead, Images, FileIO
using CuArrays

export load_image, generate_image, save_image,
	   load_model, deepdream, dream_batch,
	   load_guide_image, writevideo, recurdream

include("utils.jl")
include("dream.jl")

end # module


text:
# An exception type thrown when a tracker object is not found
type ChangeTrackerNotDefinedException <: Exception end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, BinaryBuilderBase, Pkg

include("../../fancy_toys.jl")

# The version of this JLL is decoupled from the upstream version.
# Whenever we package a new upstream release, we initially map its
# version X.Y.Z to X00.Y00.Z00 (i.e., multiply each component by 100).
# So for example version 2.6.3 would become 200.600.300.
#
# Moreover, all our packages using this JLL use `~` in their compat ranges.
#
# Together, this allows us to increment the patch level of the JLL for minor tweaks.
# If a rebuild of the JLL is needed which keeps the upstream version identical
# but breaks ABI compatibility for any reason, we can increment the minor version
# e.g. go from 200.600.300 to 200.601.300.
# To package prerelease versions, we can also adjust the minor version; e.g. we may
# map a prerelease of 2.7.0 to 200.690.000.
#
# There is currently no plan to change the major version, except when upstream itself
# changes its major version. It simply seemed sensible to apply the same transformation
# to all components.

name = "GAP"
base_version = v"400.1190.200"
upstream_version = v"4.12.0-dev"

# Collection of sources required to complete build
sources = [
    # snapshot of GAP master branch leading up to GAP 4.12:
    GitSource("https://github.com/gap-system/gap.git", "3c581ce4d4f9fb5cf68f963d2b71dd8b9f6cb349"),
#    ArchiveSource("https://github.com/gap-system/gap/releases/download/v$(upstream_version)/gap-$(upstream_version)-core.tar.gz",
#                  "2b6e2ed90fcae4deb347284136427105361123ac96d30d699db7e97d094685ce"),
    DirectorySource("../bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/gap*

for f in ${WORKSPACE}/srcdir/patches/*.patch; do
    atomic_patch -p1 ${f}
done

# run autogen.sh if compiling from it source and/or if configure was patched
./autogen.sh

# provide some generated code
cp ${WORKSPACE}/srcdir/generated/c_*.c src/

# compile GAP
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --with-gmp=${prefix} \
    --with-readline=${prefix} \
    --with-zlib=${prefix} \
    --with-gc=julia \
    --with-julia
make -j${nproc}

# install GAP binaries
make install-bin install-headers install-libgap

# also install config.h
cp build/config.h ${prefix}/include/gap

# the license
install_license LICENSE

# get rid of *.la files, they just cause trouble
rm ${prefix}/lib/*.la

# get rid of the wrapper shell script, which is useless for us
mv ${prefix}/bin/gap.real ${prefix}/bin/gap

# install gac and sysinfo.gap
mkdir -p ${prefix}/share/gap/
cp gac sysinfo.gap ${prefix}/share/gap/

# We deliberately do NOT install the GAP library, documentation, etc. because
# they are identical across all platforms; instead, we use another platform
# independent artifact to ship them to the user.
"""

function configure(julia_version, libjulia_version)
    # These are the platforms we will build for by default, unless further
    # platforms are passed in on the command line
    platforms = supported_platforms(; experimental=(julia_version > v"1.6"))

    # we only care about 64bit builds
    filter!(p -> nbits(p) == 64, platforms)

    # Windows is not supported
    filter!(!Sys.iswindows, platforms)

    # adjust the JLL version
    global version = VersionNumber(base_version.major, base_version.minor, base_version.patch + julia_version.minor)

    if julia_version >= v"1.6"
        # add julia_version to platform tuple

    #= disabled for now, until we have way to build versions of the JLL against
       Julia dev versions

        foreach(platforms) do p
            BinaryPlatforms.add_tag!(p.tags, "julia_version", string(julia_version))
        end
=#
    end

    return platforms
end

# The products that we will ensure are always built
products = [
    ExecutableProduct("gap", :gap),
    LibraryProduct("libgap", :libgap),
]

const init_block = """

    sym = dlsym(libgap_handle, :GAP_InitJuliaMemoryInterface)
    ccall(sym, Nothing, (Any, Ptr{Nothing}), @__MODULE__, C_NULL)
"""

max_stars_repo_name: singularitti/Yggdrasil

max_issues_repo_name: singularitti/Yggdrasil

max_forks_repo_name: singularitti/Yggdrasil


text:
function test_header(hdStr :: String)
    test_divider();
    println("Testing  $hdStr")
    return nothing
end

function test_divider()
    println("-------------------")
    return nothing
end

# -------------


text:
highlighter_exists(config::CthulhuConfig) =
    Sys.which(config.highlighter.exec[1]) !== nothing

function highlight(io, x, lexer, config::CthulhuConfig)
    _print = endswith(x, '\n') ? print : println
    config.enable_highlighter || return _print(io, x)
    if lexer == "llvm"
        InteractiveUtils.print_llvm(io, x)
    elseif lexer == "asm"
        InteractiveUtils.print_native(io, x)
    else
        if !highlighter_exists(config)
            @warn "Highlighter command $(config.highlighter.exec[1]) does not exist."
            return _print(io, x)
        end
        cmd = `$(config.highlighter) $lexer`
        open(pipeline(cmd; stdout=io, stderr=stderr), "w") do io
            _print(io, x)
        end
    end
end

function cthulhu_llvm(io::IO, mi, optimize, debuginfo, interp::CthulhuInterpreter, config::CthulhuConfig,
                      dump_module = false)
    dump = InteractiveUtils._dump_function_linfo_llvm(
        mi, get_world_counter(interp),
        #=wrapper=# false, #=strip_ir_metadata=# true,
        dump_module,
        optimize, debuginfo != DInfo.none ? :source : :none,
        Base.CodegenParams())
    highlight(io, dump, "llvm", config)
end

function cthulhu_native(io::IO, mi, optimize, debuginfo, interp::CthulhuInterpreter, config::CthulhuConfig)
    dump = InteractiveUtils._dump_function_linfo_native(
        mi, get_world_counter(interp),
        #=wrapper=# false, #=syntax=# config.asm_syntax,
        debuginfo != DInfo.none ? :source : :none,
        #=binary=# false)
    highlight(io, dump, "asm", config)
end

function cthulhu_ast(io::IO, mi, optimize, debuginfo, ::CthulhuInterpreter, config::CthulhuConfig)
    meth = mi.def::Method
    ast = definition(Expr, meth)
    if ast!==nothing
        if !config.pretty_ast
            dump(io, ast; maxdepth=typemax(Int))
        else
            show(io, ast)
            # Meta.show_sexpr(io, ast)
            # Could even highlight the above as some kind-of LISP
        end
    else
        @warn "Could not retrieve AST of $meth. AST display requires Revise.jl to be loaded."
    end
end

function cthulhu_source(io::IO, mi, optimize, debuginfo, ::CthulhuInterpreter, config::CthulhuConfig)
    meth = mi.def::Method
    def = definition(String, meth)
    if isnothing(def)
        return @warn "couldn't retrieve source of $meth"
    end
    src, line = def
    highlight(io, src, "julia", config)
end

using Base.IRShow: IRShow, _stmt, _type, should_print_ssa_type, IRShowConfig, show_ir

const __debuginfo = merge(IRShow.__debuginfo, Dict(
    :compact => src -> src isa CodeInfo ? __debuginfo[:source](src)
                                        : IRShow.inline_linfo_printer(src)
))

function is_type_unstable(code::Union{IRCode, CodeInfo}, idx::Int, used::BitSet)
    stmt = _stmt(code, idx)
    type = _type(code, idx)
    should_print_ssa_type(stmt) || return false
    # `used` only contains used SSA values and ignores slots
    in_use = in(idx, used) || Meta.isexpr(stmt, :(=))
    return in_use && is_type_unstable(type)
end
is_type_unstable(@nospecialize(type)) = type isa Type && (!Base.isdispatchelem(type) || type == Core.Box)

cthulhu_warntype(args...; kwargs...) = cthulhu_warntype(stdout::IO, args...; kwargs...)
function cthulhu_warntype(io::IO, debuginfo::AnyDebugInfo,
    src::Union{CodeInfo,IRCode}, @nospecialize(rt), mi::Union{Nothing,MethodInstance}=nothing;
    hide_type_stable::Bool=false, inline_cost::Bool=false, interp::CthulhuInterpreter=CthulhuInterpreter())
    if inline_cost
        isa(mi, MethodInstance) || error("Need a MethodInstance to show inlining costs. Call `cthulhu_typed` directly instead.")
    end
    cthulhu_typed(io, debuginfo, src, rt, mi; iswarn=true, hide_type_stable, inline_cost, interp)
    return nothing
end

# # for API consistency with the others
# function cthulhu_typed(io::IO, mi::MethodInstance, optimize, debuginfo, params, config::CthulhuConfig)
#     interp = mkinterp(mi)
#     (; src, rt, infos, slottypes) = lookup(interp, mi, optimize)
#     ci = Cthulhu.preprocess_ci!(src, mi, optimize, config)
#     cthulhu_typed(io, debuginfo, src, rt, mi)
# end

cthulhu_typed(io::IO, debuginfo::DebugInfo, args...; kwargs...) =
    cthulhu_typed(io, Symbol(debuginfo), args...; kwargs...)
function cthulhu_typed(io::IO, debuginfo::Symbol,
    src::Union{CodeInfo,IRCode}, @nospecialize(rt), mi::Union{Nothing,MethodInstance};
    iswarn::Bool=false, hide_type_stable::Bool=false,
    remarks::Union{Nothing,Remarks}=nothing, inline_cost::Bool=false,
    type_annotations::Bool=true, interp::CthulhuInterpreter=CthulhuInterpreter())

    debuginfo = IRShow.debuginfo(debuginfo)
    lineprinter = __debuginfo[debuginfo]
    rettype = ignorelimited(rt)
    iolim = IOContext(io, :limit=>true)
    lambda_io = iolim

    if isa(src, Core.CodeInfo)
        # we're working on pre-optimization state, need to ignore `LimitedAccuracy`
        src = copy(src)
        src.ssavaluetypes = Base.mapany(ignorelimited, src.ssavaluetypes::Vector{Any})
        src.rettype = ignorelimited(src.rettype)

        if src.slotnames !== nothing
            slotnames = Base.sourceinfo_slotnames(src)
            lambda_io = IOContext(lambda_io, :SOURCE_SLOTNAMES => slotnames)
            show_variables(io, src, slotnames)
        end
    end

    if iswarn
        print(io, "Body")
        InteractiveUtils.warntype_type_printer(iolim, rettype, true)
        println(io)
    else
        isa(mi, MethodInstance) || throw("`mi::MethodInstance` is required")
        print(io, "│ ─ ")
        println(iolim, Callsite(-1, MICallInfo(mi, rettype), :invoke))
    end

    # preprinter configuration
    if src isa IRCode && inline_cost
        isa(mi, MethodInstance) || throw("`mi::MethodInstance` is required")
        code = src isa IRCode ? src.stmts.inst : src.code
        cst = Vector{Int}(undef, length(code))
        params = Core.Compiler.OptimizationParams(interp)
        maxcost = Core.Compiler.statement_costs!(cst, code, src, Any[mi.sparam_vals...], false, params)
        nd = ndigits(maxcost)
        _lineprinter = lineprinter(src)
        function preprinter(io, linestart, idx)
            str = idx > 0 ? lpad(cst[idx], nd+1) : " "^(nd+1)
            str = sprint(io -> Base.printstyled(io, str; color=:green); context=:color=>true)
            if debuginfo === :source
                str *= " "
                linestart *= " "^(nd+2)
            end
            return str * _lineprinter(io, linestart, idx)
        end
    else
        preprinter = lineprinter(src)
    end
    # postprinter configuration
    _postprinter = if type_annotations
        iswarn ? InteractiveUtils.warntype_type_printer : IRShow.default_expr_type_printer
    else
        Returns(nothing)
    end
    if !isnothing(remarks)
        function postprinter(io::IO, @nospecialize(typ), used::Bool)
            _postprinter(io, typ, used)
            haskey(io, :idx) || return
            idx = io[:idx]::Int
            for (pc, remark) in remarks
                if pc == idx
                    printstyled(io, ' ', remark; color=:light_black)
                end
            end
        end
    else
        postprinter = _postprinter
    end

    should_print_stmt = hide_type_stable ? is_type_unstable : Returns(true)
    bb_color = (src isa IRCode && debuginfo === :compact) ? :normal : :light_black

    irshow_config = IRShowConfig(preprinter, postprinter; should_print_stmt, bb_color)

    show_ir(lambda_io, src, irshow_config)
    return nothing
end

function show_variables(io, src, slotnames)
    println(io, "Variables")
    slottypes = src.slottypes
    for i = 1:length(slotnames)
        print(io, "  ", slotnames[i])
        if isa(slottypes, Vector{Any})
            InteractiveUtils.warntype_type_printer(io, slottypes[i], true)
        end
        println(io)
    end
    println(io)
end

# These are standard code views that don't need any special handling,
# This namedtuple maps toggle::Symbol to function
const CODEVIEWS = (;
    # typed=cthulhu_typed,
    llvm=cthulhu_llvm,
    native=cthulhu_native,
    ast=cthulhu_ast,
    source=cthulhu_source,
)

"""
    Cthulhu.Bookmark

A `Cthulhu.Bookmark` remembers a method marked by `b` key during a descent.
It can be used with the following functions:

* `descend(::Bookmark)`, `descend_code_typed(::Bookmark)`,
  `descend_code_warntype(::Bookmark)`: continue the descent.
* `code_typed(::Bookmark)`, `code_warntype([::IO,] ::Bookmark)`: show typed IR
* `code_llvm([::IO,] ::Bookmark)`: pretty-print LLVM IR
* `code_native([::IO,] ::Bookmark)`: pretty-print native code
"""
struct Bookmark
    mi::MethodInstance
    interp::CthulhuInterpreter
end

"""
    Cthulhu.BOOKMARKS :: Vector{Bookmark}

During a descent, methods can be "bookmarked" by pressing `b` key.  It
pushes a [`Cthulhu.Bookmark`](@ref) into `Cthulhu.BOOKMARKS`.  This can be
used to, e.g., continue descending by `descend(Cthulhu.BOOKMARKS[end])`.
See [`Cthulhu.Bookmark`](@ref) for other usages.
"""
const BOOKMARKS = Bookmark[]

# Turn off `optimize` and `debuginfo` for default `show` so that the
# output is smaller.
function Base.show(
    io::IO, ::MIME"text/plain", b::Bookmark;
    optimize::Bool=false, debuginfo::AnyDebugInfo=:none, iswarn::Bool=false, hide_type_stable::Bool=false)
    world = get_world_counter(b.interp)
    CI, rt = InteractiveUtils.code_typed(b; optimize)
    if get(io, :typeinfo, Any) === Bookmark  # a hack to check if in Vector etc.
        print(io, Callsite(-1, MICallInfo(b.mi, rt), :invoke))
        print(io, " (world: ", world, ")")
        return
    end
    println(io, "Cthulhu.Bookmark (world: ", world, ")")
    cthulhu_typed(io, debuginfo, CI, rt, b.mi; iswarn, hide_type_stable, b.interp)
end

function InteractiveUtils.code_typed(b::Bookmark; optimize::Bool=true)
    interp = b.interp
    mi = b.mi
    (; src, rt, codeinf) = lookup(interp, mi, optimize)
    src = preprocess_ci!(src, mi, optimize, CONFIG)
    if src isa IRCode
        nargs = Int((mi.def::Method).nargs) - 1
        Core.Compiler.replace_code_newstyle!(codeinf, src, nargs)
    end
    return codeinf => rt
end

InteractiveUtils.code_warntype(b::Bookmark; kw...) =
    InteractiveUtils.code_warntype(stdout::IO, b; kw...)
function InteractiveUtils.code_warntype(
    io::IO,
    b::Bookmark;
    debuginfo::AnyDebugInfo = :source,
    hide_type_stable::Bool = true,
    kw...,
)
    CI, rt = InteractiveUtils.code_typed(b; kw...)
    cthulhu_warntype(io, debuginfo, CI, rt, b.mi; hide_type_stable, b.interp)
end

InteractiveUtils.code_llvm(b::Bookmark) = InteractiveUtils.code_llvm(stdout::IO, b)
InteractiveUtils.code_llvm(
    io::IO,
    b::Bookmark;
    optimize = true,
    debuginfo = :source,
    dump_module = false,
    config = CONFIG,
) = cthulhu_llvm(
    io,
    b.mi,
    optimize,
    debuginfo == :source,
    b.interp,
    config,
    dump_module,
)

InteractiveUtils.code_native(b::Bookmark; kw...) =
    InteractiveUtils.code_native(stdout::IO, b; kw...)
InteractiveUtils.code_native(
    io::IO,
    b::Bookmark;
    optimize = true,
    debuginfo = :source,
    config = CONFIG,
) = cthulhu_native(io, b.mi, optimize, debuginfo == :source, b.interp, config)


text:
"""
A mutable collection of `Monetary` values of various currencies.

Like regular `Monetary` values, `Basket` values support basic arithmetic
operations, with both other baskets and with raw monetary values. The
constructor for `Basket` accepts either a monetary value or a vector of
monetary values.

```julia
basket = Basket([1USD, 5EUR])        # \$1 + 5€
basket += 2EUR                       # \$1 + 7€
basket *= 2                          # \$2 + 14€
```

`Basket` values also support iteration. The iteration order of currencies is
undefined. For instance, the following will print `1.00 USD` and `5.00 EUR`, in
some order. Iteration will skip over any currencies with zero weight in the
basket.

```julia
for amt in Basket([1USD, 5EUR])
    println(amt)
end
```

One can also mutate the basket with index notation, or with `push!`.

```julia
basket = Basket([1USD, 2EUR])
basket[:USD] = 3USD  # Basket([3USD, 2EUR])
push!(basket, 10GBP) # Basket([3USD, 2EUR, 10GBP])
```
"""
struct Basket <: AbstractMonetary
    table::SMDict
    Basket(ms::Union{AbstractArray,Tuple}) = new(buildtable(ms))
end

# basket outer constructor and convert-constructor
Basket() = Basket(())
(::Type{T})(x) where T <: Basket = convert(T, x)

# access methods
Base.length(b::Basket) = length(b.table)
Base.haskey(b::Basket, k) = haskey(b.table, k) && !iszero(b.table[k])
Base.getindex(b::Basket, T::Symbol) = get(b.table, T, zero(Monetary{T}))
Base.iterate(b::Basket, s...) = iterate(values(b.table), s...)

# arithmetic methods (for static & dynamic baskets)
- b::Basket = Basket([-x for x in b])
b::Basket + c::Basket = Basket((b, c))
b::Basket * k::Real   = Basket([x * k for x in b])
k::Real   * b::Basket = Basket([k * x for x in b])
b::Basket / k::Real   = Basket([x / k for x in b])

# methods for dynamic baskets
function Base.setindex!(b::Basket, m::Monetary, k::Symbol)
    if currency(m) ≠ k
        throw(ArgumentError("Monetary value type does not match currency"))
    end
    b.table[k] = m
    deleteifzero!(b.table, k)
    b
end

# somewhat strange because Baskets behave like collections and like numbers...
# maybe add! is a better name
Base.push!(b::Basket, m::Monetary) = (b[currency(m)] += m; b)
function Base.push!(b::Basket, c::Basket)
    for m in c
        push!(b, m)
    end
    b
end

# other methods (eltype, iszero, zero, ==)
          iszero(b::Basket) = isempty(b)
  Base.zero(::Type{Basket}) = Basket()
   Base.one(::Type{Basket}) = 1
Base.eltype(::Type{Basket}) = Monetary


text:
mutable struct RegisteredFunction
    s         :: Symbol
    dimension :: Integer
    f         :: Function
    gradf     :: Union{Nothing, Function}
    grad2f    :: Union{Nothing, Function}
    autodiff  :: Bool
end

# Options for the solver (more details like defaults in solver.jl)
mutable struct SolverOptions
    nl_solver                           :: Any# needs to be set
    log_levels                          :: Vector{Symbol}
    silent                              :: Bool
    atol                                :: Float64
    num_resolve_root_relaxation         :: Int64
    branch_strategy                     :: Symbol
    gain_mu                             :: Float64
    strong_branching_perc               :: Float64
    strong_branching_nsteps             :: Int64
    strong_branching_time_limit  :: Float64
    strong_restart                      :: Bool
    reliability_branching_threshold     :: Int64
    reliability_branching_perc          :: Float64
    incumbent_constr                    :: Bool
    obj_epsilon                         :: Float64
    time_limit                          :: Float64
    mip_gap                             :: Float64
    best_obj_stop                       :: Float64
    solution_limit                      :: Int64
    all_solutions                       :: Bool
    list_of_solutions                   :: Bool
    processors                          :: Int64
    traverse_strategy                   :: Symbol
    feasibility_pump                    :: Bool
    feasibility_pump_time_limit         :: Float64
    feasibility_pump_tolerance_counter  :: Int64
    tabu_list_length                    :: Int64
    num_resolve_nlp_feasibility_pump    :: Int64
    mip_solver                          :: Any
    allow_almost_solved                 :: Bool
    allow_almost_solved_integral        :: Bool
    registered_functions                :: Union{Nothing, Vector{RegisteredFunction}}
    seed                                :: Int64

    # only for testing
    force_parallel                      :: Bool
    debug                               :: Bool
    debug_write                         :: Bool
    debug_file_path                     :: String

    fixed_gain_mu                       :: Bool
end

mutable struct SolutionObj
    solution    :: Vector{Float64}
    objval      :: Float64
end

# Juniper MOI struct

mutable struct JuniperProblem
    nl_solver           :: Any
    nl_solver_options   :: Vector{Pair}

    model               :: JuMP.Model

    relaxation_status   :: MOI.TerminationStatusCode
    relaxation_objval   :: Float64
    relaxation_solution :: Vector{Float64}

    status              :: MOI.TerminationStatusCode
    objval              :: Float64
    best_bound          :: Float64

    x                   :: Vector{JuMP.VariableRef}
    primal_start        :: Vector{Real}
    num_constr          :: Int64
    num_nl_constr       :: Int64
    num_q_constr        :: Int64
    num_l_constr        :: Int64
    num_var             :: Int64
    l_var               :: Vector{Float64}
    u_var               :: Vector{Float64}

    has_nl_objective    :: Bool
    nlp_evaluator       :: MOI.AbstractNLPEvaluator

    objective           :: Union{SVF, SAF, SQF, Nothing}

    disc2var_idx        :: Vector{Int64}
    var2disc_idx        :: Vector{Int64}

    var_type            :: Vector{Symbol}
    obj_sense           :: Symbol
    num_disc_var        :: Int64

    solution            :: Vector{Float64}

    soltime             :: Float64
    options             :: SolverOptions
    solutions           :: Vector{SolutionObj}
    nsolutions          :: Int64

    mip_solver          :: Any
    mip_solver_options  :: Vector{Pair}

    relaxation_time     :: Float64
    start_time          :: Float64

    # Info
    nintvars            :: Int64
    nbinvars            :: Int64
    nnodes              :: Int64
    ncuts               :: Int64
    nbranches           :: Int64
    nlevels             :: Int64

    fpump_info          :: Dict{Symbol,Float64}

    # debug
    debugDict           :: Dict{Symbol,Any}

    JuniperProblem() = new()
end

###########################################################################
########################## FPump ##########################################
###########################################################################
mutable struct TabuList
    sols      :: Vector{Vector{Float64}}
    length    :: Int64
    pointer   :: Int64

    TabuList() = new()
end

############################################################################
######################## Tree structure ####################################
############################################################################
mutable struct BnBNode
    idx                 :: Int64
    level               :: Int64
    l_var               :: Vector{Float64}
    u_var               :: Vector{Float64}
    solution            :: Vector{Float64}
    var_idx             :: Int64
    state               :: Symbol
    relaxation_state    :: MOI.TerminationStatusCode
    best_bound          :: Float64
    path                :: Vector{String}   # list of parent hashes
    hash                :: String
end

mutable struct Incumbent
    objval      :: Float64
    solution    :: Vector{Float64}
    only_almost :: Bool
end

mutable struct GainObj
    minus           :: Vector{Float64} # gain of objective per variable on left node
    plus            :: Vector{Float64} # gain of objective per variable on right node
    minus_counter   :: Vector{Int64} # obj_gain_m / obj_gain_mc => average gain on left node
    plus_counter    :: Vector{Int64} # obj_gain_p / obj_gain_pc => average gain on right node
    # counter of how often one child is infeasible (>= 0) subtract 1 if both Optimal
    # can be negative temporately in step_obj but not in tree.obj_gain_
    inf_counter     :: Vector{Int64}
end

mutable struct BnBTreeObj
    m               :: Juniper.JuniperProblem
    incumbent       :: Incumbent
    limit           :: Symbol
    obj_gain        :: GainObj
    disc2var_idx    :: Vector{Int64}
    var2disc_idx    :: Vector{Int64}
    options         :: Juniper.SolverOptions
    obj_fac         :: Int64 # factor for objective 1 if max -1 if min
    start_time      :: Float64
    nsolutions      :: Int64
    branch_nodes    :: Vector{BnBNode}
    best_bound      :: Float64
    global_solver   :: Bool

    BnBTreeObj() = new()
end

mutable struct StrongBranchStep
    var_idx             :: Int64
    l_relaxation_state  :: MOI.TerminationStatusCode
    r_relaxation_state  :: MOI.TerminationStatusCode
    init_restart        :: Bool

    StrongBranchStep() = new()
end

# the object holds information for the current step
mutable struct StepObj
    node                :: BnBNode # current branch node
    var_idx             :: Int64   # variable to branch on
    state               :: Symbol  # if infeasible => break (might be set by strong branching)
    nrestarts           :: Int64
    gain_gap            :: Float64
    obj_gain            :: GainObj
    idx_time            :: Float64
    node_idx_time       :: Float64
    upd_gains_time      :: Float64
    node_branch_time    :: Float64
    branch_time         :: Float64
    integral            :: Vector{BnBNode}
    branch              :: Vector{BnBNode}
    l_nd                :: BnBNode
    r_nd                :: BnBNode
    counter             :: Int64
    upd_gains           :: Symbol
    strong_disc_vars    :: Vector{Int64}
    branch_strategy     :: Symbol # the strategy that is used in this step
    strong_branching    :: Union{Nothing,Vector{StrongBranchStep}} # used only if debug=true and branch_strategy of this step is :Strong

    StepObj() = new()
end

mutable struct TimeObj
    solve_leaves_get_idx :: Float64
    solve_leaves_branch :: Float64
    branch :: Float64
    get_idx :: Float64
    upd_gains :: Float64
end