text:
############################################
################ Structures ################
############################################

function glb(types...)
    # If a lower bound is in the types then it is greatest
    # else we just return Unknown for now
    for type in types
        all(type <: t_ for t_ in types) && return type
    end
    return Unknown
end

function input_target_scitypes(models, metalearner)
    # The target scitype is defined as the greatest lower bound of the
    # metalearner and the base models in the library
    all_tg_scitypes = [target_scitype(m) for m in models]
    tg_scitype = glb(target_scitype(metalearner), all_tg_scitypes...)
    # The input scitype is defined as the greatest lower bound of the
    # base models in the library
    inp_scitype = glb([input_scitype(m) for m in models]...)

    return inp_scitype, tg_scitype
end

mutable struct DeterministicStack{modelnames, inp_scitype, tg_scitype} <: DeterministicComposite
   models::Vector{Supervised}
   metalearner::Deterministic
   resampling
   function DeterministicStack(modelnames, models, metalearner, resampling)
        inp_scitype, tg_scitype = input_target_scitypes(models, metalearner)
        return new{modelnames, inp_scitype, tg_scitype}(models, metalearner, resampling)
   end
end

mutable struct ProbabilisticStack{modelnames, inp_scitype, tg_scitype} <: ProbabilisticComposite
    models::Vector{Supervised}
    metalearner::Probabilistic
    resampling
    function ProbabilisticStack(modelnames, models, metalearner, resampling)
        inp_scitype, tg_scitype = input_target_scitypes(models, metalearner)
        return new{modelnames, inp_scitype, tg_scitype}(models, metalearner, resampling)
    end
 end

const Stack{modelnames, inp_scitype, tg_scitype} =
    Union{DeterministicStack{modelnames, inp_scitype, tg_scitype},
            ProbabilisticStack{modelnames, inp_scitype, tg_scitype}}

"""
    Stack(;metalearner=nothing, resampling=CV(), name1=model1, name2=model2, ...)

Implements the two-layer generalized stack algorithm introduced by
[Wolpert
(1992)](https://www.sciencedirect.com/science/article/abs/pii/S0893608005800231)
and generalized by [Van der Laan et al
(2007)](https://biostats.bepress.com/ucbbiostat/paper222/). Returns an
instance of type `ProbablisiticStack` or `DeterministicStack`,
depending on the prediction type of `metalearner`.

When training a machine bound to such an instance:

- The data is split into training/validation sets according to the
  specified `resampling` strategy.

- Each base model `model1`, `model2`, ... is trained on each training
  subset and outputs predictions on the corresponding validation
  sets. The multi-fold predictions are spliced together into a
  so-called out-of-sample prediction for each model.

- The adjudicating model, `metalearner`, is subsequently trained on
  the out-of-sample predictions to learn the best combination of base
  model predictions.

- Each base model is retrained on all supplied data for purposes of
  passing on new production data onto the adjudicator for making new
  predictions

### Arguments

- `metalearner::Supervised`: The model that will optimize the desired
  criterion based on its internals.  For instance, a LinearRegression
  model will optimize the squared error.

- `resampling`: The resampling strategy used
  to prepare out-of-sample predictions of the base learners.
  It can be a user-defined strategy, the only
  caveat being that it should have a `nfolds` attribute.

- `name1=model1, name2=model2, ...`: the `Supervised` model instances
  to be used as base learners.  The provided names become properties
  of the instance created to allow hyper-parameter access

### Example

The following code defines a `DeterministicStack` instance for
learning a `Continuous` target, and demonstrates that:

- Base models can be `Probabilistic` models even if the stack
  itself is `Deterministic` (`predict_mean` is applied in such cases).

- As an alternative to hyperparameter optimization, one can stack
  multiple copies of given model, mutating the hyper-parameter used in
  each copy.

```julia
using MLJ

DecisionTreeRegressor = @load DecisionTreeRegressor pkg=DecisionTree
EvoTreeRegressor = @load EvoTreeRegressor
XGBoostRegressor = @load XGBoostRegressor
KNNRegressor = @load KNNRegressor pkg=NearestNeighborModels
LinearRegressor = @load LinearRegressor pkg=MLJLinearModels

X, y = make_regression(500, 5)

stack = Stack(;metalearner=LinearRegressor(),
                resampling=CV(),
                constant=ConstantRegressor(),
                tree_2=DecisionTreeRegressor(max_depth=2),
                tree_3=DecisionTreeRegressor(max_depth=3),
                evo=EvoTreeRegressor(),
                knn=KNNRegressor(),
                xgb=XGBoostRegressor())

mach = machine(stack, X, y)
evaluate!(mach; resampling=Holdout(), measure=rmse)

```

"""
function Stack(;metalearner=nothing, resampling=CV(), named_models...)
    metalearner === nothing &&
        throw(ArgumentError("No metalearner specified. Use Stack(metalearner=...)"))

    nt = NamedTuple(named_models)
    modelnames = keys(nt)
    models = collect(nt)

    if metalearner isa Deterministic
        stack =  DeterministicStack(modelnames, models, metalearner, resampling)
    elseif metalearner isa Probabilistic
        stack = ProbabilisticStack(modelnames, models, metalearner, resampling)
    else
        throw(ArgumentError("The metalearner should be a subtype
                    of $(Union{Deterministic, Probabilistic})"))
    end
    # Issuing clean! statement
    message = MMI.clean!(stack)
    isempty(message) || @warn message

    # Warning if either input_scitype/target_scitype is
    # Unknown at construction time
    params = typeof(stack).parameters
    params[end-1] == Unknown && @warn "Could not infer input_scitype of the stack"
    params[end] == Unknown && @warn "Could not infer target_scitype of the stack"

    return stack
end

function MMI.clean!(stack::Stack{modelnames, inp_scitype, tg_scitype}) where {modelnames,inp_scitype,tg_scitype}
    # We only carry checks and don't try to correct the arguments here
    message = ""
    # Checking target_scitype and input_scitype have not been changed from the original stack
    glb_inp_scitype, glb_tg_scitype = input_target_scitypes(getfield(stack, :models), stack.metalearner)
    glb_inp_scitype == inp_scitype ||
            throw(DomainError(inp_scitype, "The newly inferred input_scitype of the stack doesn't
            match its original one. You have probably changed one of the base models or the metalearner
            to a non compatible type."))
    glb_tg_scitype == tg_scitype ||
            throw(DomainError(tg_scitype, "The newly inferred target_scitype of the stack doesn't
            match its original one. You have probably changed one of the base model or the metalearner
            to a non compatible type."))
    # Checking the target scitype is consistent with either Probabilistic/Deterministic Stack
    target_scitype(stack.metalearner) <: Union{AbstractArray{<:Continuous}, AbstractArray{<:Finite}} ||
        throw(ArgumentError("The metalearner should have target_scitype:
                $(Union{AbstractArray{<:Continuous}, AbstractArray{<:Finite}})"))

    return message
end

Base.propertynames(::Stack{modelnames}) where modelnames = tuple(:resampling, :metalearner, modelnames...)

function Base.getproperty(stack::Stack{modelnames}, name::Symbol) where modelnames
    name === :metalearner && return getfield(stack, :metalearner)
    name === :resampling && return getfield(stack, :resampling)
    models = getfield(stack, :models)
    for j in eachindex(modelnames)
        name === modelnames[j] && return models[j]
    end
    error("type Stack has no property $name")
end

function Base.setproperty!(stack::Stack{modelnames}, _name::Symbol, val) where modelnames
    _name === :metalearner && return setfield!(stack, :metalearner, val)
    _name === :resampling && return setfield!(stack, :resampling, val)
    idx = findfirst(==(_name), modelnames)
    idx isa Nothing || return getfield(stack, :models)[idx] = val
    error("type Stack has no property $name")
end

MMI.target_scitype(::Type{<:Stack{modelnames, input_scitype, target_scitype}}) where
    {modelnames, input_scitype, target_scitype} = target_scitype

MMI.input_scitype(::Type{<:Stack{modelnames, input_scitype, target_scitype}}) where
    {modelnames, input_scitype, target_scitype} = input_scitype

MLJBase.load_path(::Type{<:ProbabilisticStack}) = "MLJBase.ProbabilisticStack"
MLJBase.load_path(::Type{<:DeterministicStack}) = "MLJBase.DeterministicStack"
MLJBase.package_name(::Type{<:Stack}) = "MLJBase"
MLJBase.package_uuid(::Type{<:Stack}) = "a7f614a8-145f-11e9-1d2a-a57a1082229d"
MLJBase.package_url(::Type{<:Stack}) = "https://github.com/alan-turing-institute/MLJBase.jl"
MLJBase.package_license(::Type{<:Stack}) = "MIT"

###########################################################
################# Node operations Methods #################
###########################################################

getfolds(y::AbstractNode, cv::CV, n::Int) =
    source(train_test_pairs(cv, 1:n))

getfolds(y::AbstractNode, cv::StratifiedCV, n::Int) =
    node(YY->train_test_pairs(cv, 1:n, YY), y)

trainrows(X::AbstractNode, folds::AbstractNode, nfold) =
    node((XX, ff) -> selectrows(XX, ff[nfold][1]), X, folds)

testrows(X::AbstractNode, folds::AbstractNode, nfold) =
    node((XX, ff) -> selectrows(XX, ff[nfold][2]), X, folds)

pre_judge_transform(ŷ::Node, ::Type{<:Probabilistic}, ::Type{<:AbstractArray{<:Finite}}) =
    node(ŷ -> pdf(ŷ, levels(first(ŷ))), ŷ)

pre_judge_transform(ŷ::Node, ::Type{<:Probabilistic}, ::Type{<:AbstractArray{<:Continuous}}) =
    node(ŷ->mean.(ŷ), ŷ)

pre_judge_transform(ŷ::Node, ::Type{<:Deterministic}, ::Type{<:AbstractArray{<:Continuous}}) =
    ŷ

function oos_set(m::Stack, folds::AbstractNode, Xs::Source, ys::Source)
    Zval = []
    yval = []
    # Loop over the cross validation folds to build a training set for the metalearner.
    for nfold in 1:m.resampling.nfolds
        Xtrain = trainrows(Xs, folds, nfold)
        ytrain = trainrows(ys, folds, nfold)
        Xtest = testrows(Xs, folds, nfold)
        ytest = testrows(ys, folds, nfold)

        # Train each model on the train fold and predict on the validation fold
        # predictions are subsequently used as an input to the metalearner
        Zfold = []
        for model in getfield(m, :models)
            mach = machine(model, Xtrain, ytrain)
            ypred = predict(mach, Xtest)
            # Dispatch the computation of the expected mean based on
            # the model type and target_scytype
            ypred = pre_judge_transform(ypred, typeof(model), target_scitype(model))
            push!(Zfold, ypred)
        end

        Zfold = hcat(Zfold...)

        push!(Zval, Zfold)
        push!(yval, ytest)
    end

    Zval = MLJBase.table(vcat(Zval...))
    yval = vcat(yval...)

    Zval, yval
end

#######################################
################# Fit #################
#######################################
"""
    fit(m::Stack, verbosity::Int, X, y)
"""
function fit(m::Stack, verbosity::Int, X, y)
    n = nrows(y)

    Xs = source(X)
    ys = source(y)

    folds = getfolds(ys, m.resampling, n)

    Zval, yval = oos_set(m, folds, Xs, ys)

    metamach = machine(m.metalearner, Zval, yval)

    # Each model is retrained on the original full training set
    Zpred = []
    for model in getfield(m, :models)
        mach = machine(model, Xs, ys)
        ypred = predict(mach, Xs)
        ypred = pre_judge_transform(ypred, typeof(model), target_scitype(model))
        push!(Zpred, ypred)
    end

    Zpred = MLJBase.table(hcat(Zpred...))
    ŷ = predict(metamach, Zpred)

    # We can infer the Surrogate by two calls to supertype
    mach = machine(supertype(supertype(typeof(m)))(), Xs, ys; predict=ŷ)

    return!(mach, m, verbosity)

end


text:
using Resonance
using CairoMakie

md = CSV.read("data/wrangled.csv", DataFrame)

##

reads = CSV.read("data/read_counts.csv", DataFrame)
(fig, ax, hs) = hist(reads.count, axis=(
    xlabel="Total reads (single)",
    ylabel="# of samples"
))

save("slides/assets/kneaddatacounts.png", fig)

##

@rsubset!(md, !ismissing(:subject),
                    !ismissing(:sample),
                    )
@rsubset!(md, !startswith(:sample, "C"),
                    !startswith(:sample, "z"),
                    )

spec = metaphlan_profiles(filter(f-> contains(f, "profile") && any(s-> contains(f, s), md.sample),
    readdir("/grace/echo/analysis/biobakery3/links/metaphlan", join=true)), :species)

sns = map(samplenames(spec)) do s
    replace(s, r"_S\d+_profile"=>"")
end

idx = unique([findfirst(==(x), sns) for x in sns])

spec = CommunityProfile(abundances(spec)[:, idx], features(spec), MicrobiomeSample.(sns[idx]))

set!(spec, md)

for sampl in samples(spec)
    sn = name(sampl)
    set!(spec, sn, :richness, count(!=(0), abundances(spec[:, sampl])))
end

(fig, ax, hs) = hist([get(spec, sn, :richness) for sn in samplenames(spec)],
    axis=(
        xlabel="Richness (# of taxa)",
        ylabel="# of samples"
    ))

save("slides/assets/echo_richness.png", fig)
fig

##

fig = Figure()
ax = Axis(fig[1,1], ylabel="# of samples", xlabel = "Richness (# of taxa)")

idx = [any(col-> get(spec, s, col, false), [:age0to3mo, :age3to6mo, :age6to12mo]) for s in samplenames(spec)]
hist!(ax, [get(spec, sn, :richness, missing) for sn in samplenames(spec)][Not(idx)], label="Age > 12mo")
hist!(ax, [get(spec, sn, :richness, missing) for sn in samplenames(spec)][idx], label="Age < 12mo")
axislegend(ax)

save("slides/assets/echo_richness_byage.png", fig)
fig

##

yidx = [any(col-> get(spec, s, col, false), [:age0to3mo]) for s in samplenames(spec)]
oidx = [any(col-> get(spec, s, col, false), [:age12moplus]) for s in samplenames(spec)]

ratios = map(featurenames(spec)) do feature
    log(mean(abundances(spec[feature,oidx])) / mean(abundances(spec[feature, yidx])))
end

onlyyoung = count(isnan, ratios)
onlyold = count(isinf, ratios)

filter!(x-> !isinf(x) && !isnan(x), ratios)

#

fig = Figure(resolution=(1200, 800))
yngax = Axis(fig[1,1], title="Only in <3 mo", ylabel="# of species")
bothax = Axis(fig[1, 2:3], title="Present in both", xlabel="Log ratio abundances")
oldax = Axis(fig[1,4], title="Only in >12 mo")

barplot!(yngax, [1], [onlyyoung])
hist!(bothax, ratios)
barplot!(oldax, [1], [onlyold])

hidexdecorations!.([yngax, oldax])
bothax.ygridvisible = true
oldax.ygridvisible = true

tightlimits!(yngax, Bottom())
tightlimits!(bothax, Left(), Right(), Bottom())
tightlimits!(oldax, Bottom())

linkyaxes!(yngax, bothax, oldax)
save("slides/assets/age_taxon_ratios.png", fig)
fig
##

metabs = CSV.read("data/metabolites.csv", DataFrame)

##

humann_join("/grace/echo/analysis/biobakery3/links/humann/genefamilies",
            "/grace/echo/analysis/biobakery3/links/humann/all_genefamilies.tsv";
            file_name="genefamilies")

func = humann_profiles("/grace/echo/analysis/biobakery3/links/humann/all_genefamilies.tsv", stratified=false)

sns = map(samplenames(func)) do s
    replace(s, r"_S\d+_genefamilies"=>"")
end

idx = unique([findfirst(==(x), sns) for x in sns])

func = CommunityProfile(abundances(func)[:, idx], features(func), sns[idx])

set!(func, md)

for sampl in samples(func)
    sn = name(sampl)
    set!(func, sn, :richness, count(!=(0), abundances(func[:, sampl])))
end

##

(fig, ax, hs) = hist([get(func, sn, :richness) for sn in samplenames(func)],
    axis=(
        xlabel="Richness (# of gene functions)",
        ylabel="# of samples"
    ))

save("slides/assets/echo_richness_func.png", fig)
fig

##

fig = Figure()
ax = Axis(fig[1,1], ylabel="# of samples", xlabel = "Richness (# of gene functions)")

idx = [any(col-> get(func, s, col, false), [:age0to3mo, :age3to6mo, :age6to12mo]) for s in samplenames(func)]
hist!(ax, [get(func, sn, :richness, missing) for sn in samplenames(func)][Not(idx)], label="Age > 12mo")
hist!(ax, [get(func, sn, :richness, missing) for sn in samplenames(func)][idx], label="Age < 12mo")
axislegend(ax)

save("slides/assets/echo_richness__func_byage.png", fig)
fig

ext: jl

max_stars_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_issues_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_forks_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -6.7172365869297215 <= q <= 3.8657529446129524)

                     add_NL_constraint(m, :(swish(0.10971879992509148 + -0.31290714990830937*$(x[1]) + 0.781487131205628*$(x[2])) + swish(-0.5708088154317625 + 0.5853291463969379*$(x[1]) + -0.9903406373931247*$(x[2])) + swish(0.5809119298414265 + 0.5024686961828384*$(x[1]) + -0.4096226237713658*$(x[2])) + swish(-0.6435921083509291 + 0.52943737688555*$(x[1]) + 0.2727617072630939*$(x[2])) + swish(-0.9019716271422107 + 0.4041768804112231*$(x[1]) + 0.5029634163532655*$(x[2])) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
include("common.jl")

# Produce the data needed for the plots of the paper

###
### Cohen-Bergstresser
###

function cb_gamma_versus_Ecut_Double64(;n_bands=8)
    reset_timer!(to)

    h5name = "cb_gamma_Ecut_$(n_bands)bands_Double64_fast_gap.hdf5"
    isfile(h5name) && error("Remove $h5name first")
    for Ecut in append!(collect(3:19), collect(20:2:70))
        println("#\n#-- Ecut = $Ecut\n#")
        tol = 1e-28
        Ecut < 35 && (tol = 1e-25)
        Ecut < 22 && (tol = 5e-17)

        fast = false
        (Ecut > 55) && (fast = true)

        data = run_cohen_bergstresser(Ecut=Ecut, T=Double64, tol=tol, only_gamma=true,
                                      n_bands=n_bands)
        errors1 = estimate_discretisation_error(data, fast=fast)
        store_results(h5name, merge(data, errors1))
    end

    display(to)
end

function cb_gamma_versus_Ecut_Float64(;n_bands=8)
    reset_timer!(to)

    h5name = "cb_gamma_Ecut_$(n_bands)bands_Float64.hdf5"
    isfile(h5name) && error("Remove $h5name first")
    for Ecut in collect(3:70)
        println("#\n#-- Ecut = $Ecut\n#")
        tol = 1e-10
        Ecut > 12 && (tol = 1e-12)

        data = run_cohen_bergstresser(Ecut=Ecut, T=Float64, tol=tol, only_gamma=true,
                                      n_bands=n_bands)
        errors1 = estimate_discretisation_error(data, fast=false)
        errors2 = estimate_arithmetic_error(data)
        store_results(h5name, merge(data, errors1, errors2))
    end

    display(to)
end

function cb_bands(;Ecut=10, T=Float64, tol=6)
    reset_timer!(to)

    h5name = "cb_bands.hdf5"
    isfile(h5name) && error("Remove $h5name first")

    tol = 10.0^(-tol)
    data = run_cohen_bergstresser(Ecut=Ecut, T=T, tol=tol, kline_density=7)
    errors1 = estimate_discretisation_error(data, fast=true)
    errors2 = estimate_arithmetic_error(data)
    store_results(h5name, merge(data, errors1, errors2))

    display(to)
end

###
### Linear silicon
###

function si_gamma_versus_Ecut_full_bands_Float64(;maxEcut=36)
    reset_timer!(to)
    h5name = "si_gamma_Ecut_full_bands_Float64.hdf5"
    isfile(h5name) && error("Remove $h5name first")

    for Ecut in collect(6:2:maxEcut)
        println("#\n#-- Ecut = $Ecut\n#")
        tol = 1e-10
        Ecut > 45 && (tol = 5e-12)

        data = run_linear_silicon(Ecut=Ecut, T=Float64, tol=tol, only_gamma=true)
        errors1 = estimate_discretisation_error(data, full_block=true)
        store_results(h5name, merge(data, errors1))
    end

    display(to)
end

function si_gamma_versus_Ecut_Float64(;n_bands=25, maxEcut=100)
    reset_timer!(to)

    h5name = "si_gamma_Ecut_$(n_bands)bands_Float64.hdf5"
    isfile(h5name) && error("Remove $h5name first")
    for Ecut in collect(6:2:maxEcut)
        println("#\n#-- Ecut = $Ecut\n#")
        tol = 1e-10
        Ecut > 45 && (tol = 5e-12)

        data = run_linear_silicon(Ecut=Ecut, T=Float64, tol=tol, only_gamma=true,
                                  n_bands=n_bands)
        errors1 = estimate_discretisation_error(data, fast=false)
        store_results(h5name, merge(data, errors1))
    end

    display(to)
end
si_gamma_versus_Ecut_25_bands_Float64() = si_gamma_versus_Ecut_Float64(n_bands=25, maxEcut=100)
si_gamma_versus_Ecut_50_bands_Float64() = si_gamma_versus_Ecut_Float64(n_bands=50, maxEcut=60)

function si_gamma_residual_versus_Ecut2(;Ecut=20)
    reset_timer!(to)

    h5name = "si_gamma_residual_versus_Ecut2.hdf5"
    isfile(h5name) && error("Remove $h5name first")

    # Determine minimal Ecut2:
    basis = PlaneWaveBasis(model_linear_silicon(), Ecut)
    minEcut2 = ceil(minimal_Ecut2(basis))

    for Ecut2 in collect(minEcut2:100)
        sEcut2 = @sprintf "%.3f" Ecut2
        println("#\n#-- Ecut2 = $sEcut2\n#")
        data = run_linear_silicon(Ecut=Ecut, T=Float64, tol=1e-10, only_gamma=true)
        errors1 = estimate_discretisation_error(data, fast=false, Ecut2=Ecut2)
        store_results(h5name, merge(data, errors1), subdir=sEcut2)
    end

    display(to)
end

function si_bands(;Ecut=42, T=Float64, tol=6, n_bands=35)
    reset_timer!(to)

    h5name = "si_bands.hdf5"
    isfile(h5name) && error("Remove $h5name first")

    tol = 10.0^(-tol)
    data = run_linear_silicon(Ecut=Ecut, T=T, tol=tol, kline_density=7, n_bands=n_bands)
    errors1 = estimate_discretisation_error(data, fast=true)
    errors2 = estimate_arithmetic_error(data, n_bands=10)  # Otherwise too expensive
    store_results(h5name, merge(data, errors1, errors2), n_bands=10)

    display(to)
end

function main()
    functions = [
        cb_gamma_versus_Ecut_Double64,
        cb_gamma_versus_Ecut_Float64,
        cb_bands,
        #
        si_gamma_residual_versus_Ecut2,
        si_gamma_versus_Ecut_25_bands_Float64,
        si_gamma_versus_Ecut_50_bands_Float64,
        si_gamma_versus_Ecut_full_bands_Float64,
        si_bands,
    ]

    for fun in functions
        println("#\n#-- Running $fun\n#")
        fun()
        println("\n#\n#" * "-" ^ 70 * "\n#\n")
    end
end

(abspath(PROGRAM_FILE) == @__FILE__) && main()


text:
using Coverage
# defaults to src/; alternatively, supply the folder name as argument
coverage = process_folder()
# Get total coverage for all Julia files
covered_lines, total_lines = get_summary(coverage)
# Or process a single file
@show get_summary(process_file("src/Q.jl"))
@show get_summary(process_file("src/_k.jl"))


text:
abstract type AbstractBinomialModelState end

struct StateVariables
    s::OptimizingIR.ImmutableVariable
    u::OptimizingIR.ImmutableVariable
    d::OptimizingIR.ImmutableVariable
    ln_u::OptimizingIR.ImmutableVariable
    ln_d::OptimizingIR.ImmutableVariable
end

struct BinomialModelState{T<:InterestRates.DailyDatesRange} <: AbstractBinomialModelState
    dates_range::T
    current_time_step::Int # starts at n = 1
    current_state_of_nature::Int # starts at i = 1, 1 <= i <= n
    state_variables::StateVariables
end

struct BinomialModelInitialState <: AbstractBinomialModelState
end

function initial_state(ctx::CompilerContext{B}) where {B<:BinomialModelDaily}
    return BinomialModelInitialState()
end

"Returns tuple var_u, var_d, var_ln_u, var_ln_d"
function volatility_match(ctx::CompilerContext{B}, dates_range::InterestRates.DailyDatesRange, state::BinomialModelInitialState) where {B<:BinomialModelDaily}
    pricing_model = get_pricing_model(ctx)
    σ = lower!(ctx, Volatility(pricing_model.underlying), state)
    Δt = lower!(ctx, Konst(InterestRates.yearfractionvalue(dates_range)), state)
    u_value = add_instruction!(ctx, OP_CRR_VOLATILITY_MATCH_UP_FACTOR, σ, Δt) # u = exp(σ*sqrt(Δt))
    d_value = add_instruction!(ctx, OP_BIN_SCALAR_DIV, lower!(ctx, Konst(1.0), state), u_value) # d = 1 / u
    ln_u_value = add_instruction!(ctx, OP_UNARY_SCALAR_LOG, u_value) # ln(u)
    ln_d_value = add_instruction!(ctx, OP_UNARY_SCALAR_SUB, ln_u_value) # ln(d) = -ln(u)
    return bind_local_immutable_variable!(ctx, :u, u_value), bind_local_immutable_variable!(ctx, :d, d_value), bind_local_immutable_variable!(ctx, :ln_u, ln_u_value), bind_local_immutable_variable!(ctx, :ln_d, ln_d_value)
end

function lower!(ctx::CompilerContext{B}, c::WhenAt, initial_state::BinomialModelInitialState) :: OptimizingIR.ImmutableValue where {B<:BinomialModelDaily}
    return lower_binomial_tree_algorithm!(ctx, c, initial_state)
end

function lower!(ctx::CompilerContext{B}, c::Anytime, initial_state::BinomialModelInitialState) :: OptimizingIR.ImmutableValue where {B<:BinomialModelDaily}
    return lower_binomial_tree_algorithm!(ctx, c, initial_state)
end

@inline function lower_binomial_tree_algorithm!(ctx::CompilerContext{B}, c::C, initial_state::BinomialModelInitialState) :: OptimizingIR.ImmutableValue where {B<:BinomialModelDaily, C<:Union{WhenAt,Anytime}}
    startdate = get_pricing_date(ctx)
    enddate = C <: WhenAt ? c.maturity : c.expiry
    model = get_pricing_model(ctx)
    dates_range = InterestRates.DailyDatesRange(startdate, enddate, model.numeraire_daycount_convention)
    number_of_steps = length(dates_range)

    # initial underlying price
    var_s = bind_local_immutable_variable!(ctx, :s, lower!(ctx, Unit(model.underlying), initial_state))

    # up and down underlying price factors
    var_u, var_d, var_ln_u, var_ln_d = volatility_match(ctx, dates_range, initial_state)

    state_variables = StateVariables(var_s, var_u, var_d, var_ln_u, var_ln_d)

    # evaluate contract on maturity for all states
    # at maturity, there are n states of nature
    # where n is the number of steps in the dates_range
    opt_prices_at_maturity = [ lower!(ctx, c.c, BinomialModelState(dates_range, number_of_steps, i, state_variables)) for i in 1:number_of_steps ]

    current_state_prices = opt_prices_at_maturity
    for time_step in (number_of_steps-1):-1:1
        next_state_prices = current_state_prices
        current_state_prices = Vector{OptimizingIR.ImmutableValue}(undef, time_step) # could be optimized to avoid allocations

        q = lower_risk_neutral_probability!(ctx, dates_range, time_step, state_variables)
        df = lower_numeraire_discountfactor_forward!(ctx, dates_range, time_step)

        # for each state of nature
        for i in 1:time_step
            current_state = BinomialModelState(dates_range, time_step, i, state_variables)
            current_state_prices[i] = lower_current_state_price!(ctx, c, current_state, df, q, next_state_prices[i+1], next_state_prices[i])
        end
    end

    return current_state_prices[1]
end

@inline function lower_current_state_price!(ctx::CompilerContext, c::WhenAt, current_state::BinomialModelState, df, q, next_state_price_up, next_state_price_down)
    # european-style
    return add_instruction!(ctx, OP_BINOMIAL_RISK_NEUTRAL_PRICING_STEP, df, q, next_state_price_up, next_state_price_down)
end

@inline function lower_current_state_price!(ctx::CompilerContext, c::Anytime, current_state::BinomialModelState, df, q, next_state_price_up, next_state_price_down)
    # american-style
    early_exercise_value = lower!(ctx, c.c, current_state)
    hold_value = add_instruction!(ctx, OP_BINOMIAL_RISK_NEUTRAL_PRICING_STEP, df, q, next_state_price_up, next_state_price_down)
    return add_instruction!(ctx, OP_BIN_SCALAR_MAX, hold_value, early_exercise_value)
end

@inline function lower_numeraire_discountfactor_forward!(ctx::CompilerContext{B}, dates_range::InterestRates.DailyDatesRange, current_time_step::Integer) :: OptimizingIR.ImmutableValue where {B<:BinomialModelDaily}
    model = get_pricing_model(ctx)

    fwd_start_date = dates_range[current_time_step]
    fwd_end_date = dates_range[current_time_step + 1]

    numeraire_carry_symbol = get_riskfree_curve_symbol(model.static_model, ctx.attr["riskfree_curves"], SpotCurrency(model.static_model.functional_currency))
    df_numeraire_carry = DiscountFactorForward(numeraire_carry_symbol, fwd_start_date, fwd_end_date)
    return lower!(clone_context_with_model(ctx, model.static_model), df_numeraire_carry, get_pricing_date(ctx))
end

@inline function lower_risk_neutral_probability!(ctx::CompilerContext{B}, dates_range::InterestRates.DailyDatesRange, current_time_step::Integer, state_variables::StateVariables) :: OptimizingIR.ImmutableValue where {B<:BinomialModelDaily}

    model = get_pricing_model(ctx)

    fwd_start_date = dates_range[current_time_step]
    fwd_end_date = dates_range[current_time_step + 1]

    numeraire_carry_symbol = get_riskfree_curve_symbol(model.static_model, ctx.attr["riskfree_curves"], SpotCurrency(model.static_model.functional_currency))

    # For Stocks, the carry is the dividend yield
    underlying_carry_symbol = get_riskfree_curve_symbol(model.static_model, ctx.attr["riskfree_curves"], model.underlying)

    # exp(r-q)Δt = df_underlying_carry / df_numeraire_carry
    df_underlying_carry = DiscountFactorForward(underlying_carry_symbol, fwd_start_date, fwd_end_date)
    df_numeraire_carry = DiscountFactorForward(numeraire_carry_symbol, fwd_start_date, fwd_end_date)
    ert = lower!(clone_context_with_model(ctx, model.static_model), df_underlying_carry / df_numeraire_carry, get_pricing_date(ctx))

    d = state_variables.d
    u = state_variables.u

    return add_instruction!(ctx, OP_BINOMIAL_RISK_NEUTRAL_PROBABILITY, d, ert, u)
end

function lower!(ctx::CompilerContext{B}, c::Unit{Stock}, state::BinomialModelState) :: OptimizingIR.ImmutableValue where {B<:BinomialModelDaily}
    model = get_pricing_model(ctx)
    @assert c.o == model.underlying "Expected $(model.underlying) as underlying. Got. $(c.o)."
    s = state.state_variables.s
    ln_u = state.state_variables.ln_u
    ln_d = state.state_variables.ln_d
    n = lower!(ctx, Konst(state.current_time_step), state)
    i = lower!(ctx, Konst(state.current_state_of_nature), state)
    return add_instruction!(ctx, OP_BINOMIAL_UNDERLYING_PRICE_AT, s, n, i, ln_u, ln_d)
end

function lower!(ctx::CompilerContext{B}, v::Volatility, ::BinomialModelInitialState) :: OptimizingIR.ImmutableValue where {B<:BinomialModelDaily}
    return add_input_riskfactor!(ctx, v)
end

function lower!(ctx::CompilerContext{B}, unit::Unit{Stock}, ::BinomialModelInitialState) where {B<:BinomialModelDaily}
    # gets the stock price at the pricing_date, using the StaticHedgingModel
    model = get_pricing_model(ctx)
    pricing_date = get_pricing_date(ctx)
    static_context = clone_context_with_model(ctx, model.static_model)
    return lower!(static_context, unit, pricing_date)
end

function lower!(ctx::CompilerContext{BinomialModelDaily{C,U,T}}, c::Unit{SpotCurrency{C}}, s::AbstractBinomialModelState) :: OptimizingIR.ImmutableValue where {C<:Currencies.Currency, U, T}
    # Same currency -> no exch.
    # Acquisition at state `s`, possibly a future date.
    return lower!(ctx, Konst(1.0), s)
end


text:
## Y matrix function

"""
    Ymatrix(distributions::Vector{<:SimpleSDMLayer}; transform::Bool=false, observed::Bool=false)

Creates a matrix Y (a site-by-species community data table) from a vector of
layers with species distributions. Sites (rows) are ordered as in the layer
grids and species (columns) are ordered as in the vector of layers.

In the matrix Y, values for sites where at least one species is present are
either 0.0 for absence and 1.0 for presence. On the other hand, values for sites
where no species are present are all `nothing`. These sites can be removed with
`observed=true` (default is `false`).

Setting `transform=true` (default is `false`) will apply the Hellinger
transformation to the Y matrix using the function `decostand` from the package
`vegan` in *R*. However, as this step is only needed to compute LCBD scores, it
is recommended to keep the default and use the `transform` argument from the
`lcbd` function instead.
"""
function Ymatrix(
    distributions::Vector{<:SimpleSDMLayer}; transform::Bool=false, observed::Bool=false
)
    # Get distributions as vectors
    distributions_vec = [vec(d.grid) for d in distributions]
    # Create matrix Y by combining distribution vectors
    Y = hcat(distributions_vec...)

    # Get indices of sites with observations
    inds_obs = _indsobs(Y)
    # Create matrix Yobs with observed sites only, replace nothings by zeros
    Yobs = _Yobs(Y, inds_obs)
    # Apply Hellinger transformation (using vegan in R)
    if transform
        Yobs = _Ytransf(Yobs)
    end
    # Replace values in original matrix
    Y[inds_obs, :] = Yobs

    if observed
        return Yobs
    else
        return Y
    end
end

function _indsobs(Y)
    # Verify if sites have observations
    sites_obs = [any(y -> !isnothing(y) && y > 0, yrow) for yrow in eachrow(Y)]
    # Get indices of sites with observations
    inds_obs = findall(sites_obs)
    return inds_obs
end

function _indsnotobs(Y)
    # Verify if sites have observations
    sites_obs = [any(y -> !isnothing(y) && y > 0, yrow) for yrow in eachrow(Y)]
    # Get indices of sites without observations
    inds_notobs = findall(.!sites_obs)
    return inds_notobs
end

function _Yobs(Y, inds_obs)
    # Create matrix Yobs with observed sites only
    Yobs = Y[inds_obs, :]
    # Replace nothings by zeros for observed sites (~true absences)
    replace!(Yobs, nothing => 0.0)
    return Yobs
end
_Yobs(Y) = _Yobs(Y, _indsobs(Y))

function _Ytransf(Yobs)
    # Remove type Nothing from Array (weird effects with RCall)
    Yobs = Array{Float32}(Yobs)
    # Apply Hellinger transformation (using vegan in R)
    @rput Yobs
    begin
        R"""
            Ytransf <- vegan::decostand(Yobs, "hel")
        """
    end
    @rget Ytransf
    return Ytransf
end

## Richness functions

"""
    richness(Y::Matrix, layer::SimpleSDMLayer)

Computes the species richness for the sites in the community matrix `Y` and
returns a `SimpleSDMResponse` with the same dimensions and coordinates and
dimensions as `layer`.
"""
function richness(Y::Matrix, layer::SimpleSDMLayer)
    @assert size(Y, 1) == length(layer.grid) "Y and layer must have the same number of sites"
    # Create necessary Y elements
    inds_obs = _indsobs(Y)
    Yobs = _Yobs(Y, inds_obs)
    # Create empty empty vector
    sums = fill(nothing, size(Y, 1)) |> Array{Union{Nothing,Float32}}
    # Get number of species per observed site
    sums_obs = map(sum, eachrow(Yobs))
    sums[inds_obs] = sums_obs
    # Reshape to grid format
    sums = reshape(sums, size(layer)) |> Array
    ## Create SimpleSDMLayer
    return SimpleSDMResponse(sums, layer)
end

"""
    gamma(Y::Matrix)

Returns the gamma diversity in the region covered by the matrix `Y`, i.e. the
number of species with a least one presence in `Y`. Note that this will differ
from the number of columns when some species have no observations in the region.
"""
function gamma(Y::Matrix)
    inds_obs = _indsobs(Y)
    Yobs = _Yobs(Y, inds_obs)
    sp_counts = sum(Yobs; dims=1)
    return sum(sp_counts .> 0)
end

## Beta diversity function

"""
    lcbd(Y::Matrix, layer::SimpleSDMLayer; transform::Bool=true, relative::Bool=true)

Computes the LCBD (local contributions to beta diversity) scores for the sites
in the matrix `Y` and returns a `SimpleSDMResponse` with the same dimensions and
coordinates as `layer`.

Setting `transform=true` (the default) will apply the Hellinger transformation
to the Y matrix using the function `decostand` from the package `vegan` in *R*,
which is appropriate for presence-absence data.

Setting `relative=true` (default is `false`) will rescale the LCBD scores as
relative to the maximum score (whose value will be one).
"""
function lcbd(Y::Matrix, layer::SimpleSDMLayer; transform::Bool=true, relative::Bool=false)
    @assert size(Y, 1) == length(layer.grid) "Y and layer must have the same number of sites"

    # Create necessary Y elements
    inds_obs = _indsobs(Y)
    Yobs = _Yobs(Y, inds_obs)

    # Apply hellinger transformation (if requested)
    if transform
        Yobs = _Ytransf(Yobs)
    end
    # Compute beta diversity statistics
    BDstats = betadiv(Yobs)

    # Extract LCBD values
    LCBDvals = BDstats.LCBDi
    # Scale LCBDi values to maximum value
    if relative
        LCBDvals = LCBDvals ./ maximum(LCBDvals)
    end

    # Create empty grid
    LCBDgrid = fill(nothing, size(layer)) |> Array{Union{Nothing,Float32}}
    # Fill-in grid
    LCBDgrid[inds_obs] = LCBDvals
    # Create SimpleSDMLayer with LCBD values
    LCBDlayer = SimpleSDMResponse(LCBDgrid, layer)

    return LCBDlayer
end

"""
    beta_total(Y::Matrix; transform::Bool=true)

Returns the total beta diversity from the matrix `Y`, i.e. the unbiased &
comparable estimator of variance in `Y`.
"""
function beta_total(Y::Matrix; transform::Bool=true)
    # Create necessary Y elements
    inds_obs = _indsobs(Y)
    Yobs = _Yobs(Y, inds_obs)

    # Apply hellinger transformation
    if transform
        Yobs = _Ytransf(Yobs)
    end
    # Compute beta diversity statistics
    BDstats = betadiv(Yobs)

    return BDstats.BDtotal
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
using Altro
using TrajectoryOptimization
using Test
using Altro: ALTROSolver
using RobotDynamics
const RD = RobotDynamics
const TO = TrajectoryOptimization
if !isdefined(Main,:TEST_TIME)
    TEST_TIME = true
end
TEST_TIME = false
v = true
ci = haskey(ENV, "CI")

## Double Integrator
v && println("Double Integrator")
solver = ALTROSolver(Problems.DoubleIntegrator()...)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 < 1
@test max_violation(solver) < 1e-6
@test iterations(solver) ∈ [8,9] # 8
@test solver.stats.gradient[end] < 1e-8
@test status(solver) == Altro.SOLVE_SUCCEEDED

## Pendulum
v && println("Pendulum")
solver = ALTROSolver(Problems.Pendulum()..., verbose=2)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 < 2
@test max_violation(solver) < 1e-6
@test iterations(solver) == 18 # 17
@test solver.stats.gradient[end] < 1e-1
@test status(solver) == Altro.SOLVE_SUCCEEDED

## Cartpole
v && println("Cartpole")
solver = ALTROSolver(Problems.Cartpole()..., save_S=true, verbose=2)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 <  10
@test max_violation(solver) < 1e-6
@test iterations(solver) == 41 # 40
@test solver.stats.gradient[end] < 1e-1
@test status(solver) == Altro.SOLVE_SUCCEEDED
@test Altro.usestatic(Altro.get_ilqr(solver)) == true

##
solver = ALTROSolver(Problems.Cartpole()..., projected_newton=false, use_static=Val(false))
b = benchmark_solve!(solver)
iters = iterations(solver)
J = cost(solver)
if !Sys.iswindows() && VERSION > v"1.5"
    @test b.allocs == 0
end
@test Altro.usestatic(Altro.get_ilqr(solver)) == false
@test get_trajectory(solver) isa SampledTrajectory{Any,Any}

# Use Static arrays
solver = ALTROSolver(Problems.Cartpole()..., projected_newton=false, use_static=Val(true))
b = benchmark_solve!(solver)
@test iterations(solver) == iters
@test cost(solver) ≈ J
if !Sys.iswindows() && VERSION > v"1.5"
    @test b.allocs == 0
end

solver = Altro.iLQRSolver(Problems.Cartpole()...)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time /1e6 < 10
!Sys.iswindows() && (@test b.allocs == 0)
@test iterations(solver) == 47
@test status(solver) == Altro.SOLVE_SUCCEEDED

## Acrobot
if !ci
    v && println("Acrobot")
    solver = ALTROSolver(Problems.Acrobot()...)
    b = benchmark_solve!(solver)
    TEST_TIME && @test minimum(b).time / 1e6 < 15
    @test max_violation(solver) < 1e-6
    @test iterations(solver) == 50 # 50
    @test solver.stats.gradient[end] < 1e-1
    @test status(solver) == Altro.SOLVE_SUCCEEDED
end

## Parallel Park
if !ci
    v && println("Parallel Park")
    solver = ALTROSolver(Problems.DubinsCar(:parallel_park)...)
    b =  benchmark_solve!(solver)
    TEST_TIME && @test minimum(b).time /1e6 < 10
    @test max_violation(solver) < 1e-6
    @test iterations(solver) == 13 # 13
    @test solver.stats.gradient[end] < 1e-1
    @test status(solver) == Altro.SOLVE_SUCCEEDED

    ## Three Obstacles
    solver = ALTROSolver(Problems.DubinsCar(:three_obstacles)...)
    b = benchmark_solve!(solver)  # TODO: figure out why the projected newton allocates here
    TEST_TIME && @test minimum(b).time /1e6 < 6
    @test max_violation(solver) < 1e-6
    @test iterations(solver) == 20 # 20
    @test solver.stats.gradient[end] < 1e-1  # 1e-2
    @test status(solver) == Altro.SOLVE_SUCCEEDED

    solver = ALTROSolver(Problems.DubinsCar(:three_obstacles)..., projected_newton=false)
    @test solver.opts.projected_newton == false
    @test solver.stats.gradient[end] < 1e-1
    b = benchmark_solve!(solver)
    if !Sys.iswindows()  # not sure why this fails on Windows?
        @test b.allocs == 0
        @test status(solver) == Altro.SOLVE_SUCCEEDED
    end
end

## Escape
v && println("Escape")
solver = ALTROSolver(Problems.DubinsCar(:escape)..., infeasible=true, R_inf=0.1)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 < 35  # was 25
@test max_violation(solver) < 1e-5
@test iterations(solver) == 14 # 13
@test solver.stats.gradient[end] < 1e-1
@test status(solver) == Altro.SOLVE_SUCCEEDED

## Zig-zag
v && println("Quadrotor")
solver = ALTROSolver(Problems.Quadrotor(:zigzag)...)
b = benchmark_solve!(solver)
TEST_TIME && @test minimum(b).time / 1e6 < 60
@test max_violation(solver) < 1e-6
@test iterations(solver) == 14 # 16
@test solver.stats.gradient[end] < 2e-2
@test status(solver) == Altro.SOLVE_SUCCEEDED

if !ci
    v && println("Quadrotor (no pn)")
    solver = ALTROSolver(Problems.Quadrotor(:zigzag)..., projected_newton=false)
    b = benchmark_solve!(solver)
    @test iterations(solver) - 60 <= 2 # 60
    @test status(solver) == Altro.SOLVE_SUCCEEDED
    @test solver.stats.gradient[end] < 0.3
    if !Sys.iswindows() && VERSION > v"1.5"  # not sure why this fails on Windows?
        b.allocs == 0
    end

    # # Test infeasible Quadrotor (note that this allocates for the static-bp)
    # v && println("Quadrotor (infeasible)")
    # solver = ALTROSolver(Problems.Quadrotor(:zigzag)..., projected_newton=false,
    #     infeasible=true, static_bp=false, constraint_tolerance=1e-4, verbose=2)
    # b = benchmark_solve!(solver)
    # @test iterations(solver) == 25
    # @test max_violation(solver) < 1e-4
    # @test solver.stats.gradient[end] < 1e-4
end

# Barrell Roll
if !ci
    v && println("Barrell Roll")
    solver = ALTROSolver(Problems.YakProblems(costfun=:QuatLQR, termcon=:quatvec)..., use_static=Val(true))
    b = benchmark_solve!(solver)
    TEST_TIME && @test minimum(b).time / 1e6 < 100
    @test max_violation(solver) < 1e-6
    @test iterations(solver) == 20 # 20
    @test solver.stats.gradient[end] < 1e0  # 1e-3
    @test status(solver) == Altro.SOLVE_SUCCEEDED
end


text:
function generateValues( attribute::Attribute, nVals::Int )

    pVals = attribute.initValueWeights ./ sum( attribute.initValueWeights )
    return attribute.initValues[rand( attribute.initRNG, Categorical( pVals ),
        nVals )]

end  # function generateValues( attribute, nVals )


text:
using MagSim
using Plots;
gr();
using DataFrames
using CSV
using LaTeXStrings

# %%

βs = range(1e-3, stop = 1, length = 200);
Ls = [4, 8, 16, 32, 64];
prefix = "2020-02-29-simulating-magnetism-part-3";

@time sims = simulate_Ising_parallel(metropolis!, Ls, βs);
df = sim_dict_to_df(sims);
CSV.write(prefix * ".csv", df);

# %%

function plot_data(df::DataFrame)
    plots = Dict()
    Os = [:e, :c, :m, :χ]
    labels = [L"e", L"c", L"|m|", L"\chi"]
    Ls = unique(df.L)

    for (i, O) in enumerate(Os)
        O ∈ [:e, :m] ? legend_bool = true : legend_bool = false
        p = plot(
            dpi = 200,
            xaxis = (latexstring("\\beta")),
            yaxis = (labels[i]),
            legend = legend_bool,
        )

        colors = [:firebrick1, :orange, :forestgreen, :cyan3, :dodgerblue]
        for (i, L) in enumerate(Ls)
            df_ = df[(df.L.==L), :]
            x = df_[:, :β]
            y = df_[:, O]
            plot!(
                x,
                y,
                marker = (:star4, colors[i], 4, stroke(0)),
                line = (colors[i], 2),
                label = latexstring("L = $(string(L))"),
            )
        end

        plots[O] = p
    end

    return plots
end

plots = plot_data(df);
plot_e_c = plot(plots[:e], plots[:c], layout = (2, 1), size = (600, 500));
savefig(plot_e_c, "../images/" * prefix * "/ising_metropolis_e_c.png")
plot_m_χ = plot(plots[:m], plots[:χ], layout = (2, 1), size = (600, 500));
savefig(plot_m_χ, "../images/" * prefix * "/ising_metropolis_m_chi.png")


text:
type CompositeQueryResultItem
	# A result item part of an overall set of result of a CompositeQuery
	id::UUID.Uuid
	typelabel::String
	item::Container
	containedobject::Any
	hassource::Bool
	sourceresultitem::CompositeQueryResultItem
	storedvalues::Dict{String, Any}
	isnew::Bool

	function CompositeQueryResultItem(newitem::Container, source::CompositeQueryResultItem)
		# Constructs a ResultItem given an item (Vertex or Edge) and a previous result
		r = new()

		r.item = newitem
		r.id = newitem.id
		r.typelabel = newitem.typelabel
		if isdefined(newitem,:containedobject)
			r.containedobject = newitem.containedobject
		end
		r.hassource = true
		r.sourceresultitem = source
		r.storedvalues = copy(source.storedvalues)
		r.isnew = true

		return r
	end

	function CompositeQueryResultItem(newitem::Container)
		# Constructs a ResultItem given a container (Vertex or Edge)
		r = new()

		r.item = newitem
		r.id = newitem.id
		r.typelabel = newitem.typelabel
		r.containedobject = UnspecifiedValue
		r.hassource = false
		r.storedvalues = Dict{String,Any}()
		r.isnew = true

		return r
	end
end

function get(resultitem::CompositeQueryResultItem, name::String, default::Any=UnspecifiedValue)
	# get the value of a specified property, returning a default value if there is no other property value associated
	return get(resultitem.item,name,default)
end

function getstored(resultitem::CompositeQueryResultItem, name::String, default::Any=UnspecifiedValue)
	# retrieve a stored value
	return Base.get(resultitem.storedvalues,name,default)
end

function store(resultitem::CompositeQueryResultItem, name::String, value::Any)
	# store a value on the as part of the result item
	resultitem.storedvalues[name] = value
end


text:
export iauPdp
"""
p-vector inner (=scalar=dot) product.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   a      double[3]     first p-vector
   b      double[3]     second p-vector

Returned (function value):
         double        a . b

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauPdp(double a[3], double b[3])
function iauPdp(a::AbstractVector{<:Real}, b::AbstractVector{<:Real})
   return ccall((:iauPdp, libsofa_c), Cdouble,
       (Ptr{Cdouble}, Ptr{Cdouble}),
       convert(Array{Float64, 1}, a),
       convert(Array{Float64, 1}, b))
end


text:
# Compute C = A * B, using a cache-obvious multiplication algorithm

function add_matmul_rec!(m,n,p, i0,j0,k0, C,A,B)
    if m+n+p <= 64   # base case: naive matmult for sufficiently large matrices
        for i = 1:m
            for k = 1:p
                c = zero(eltype(C))
                for j = 1:n
                    @inbounds c += A[i0+i,j0+j] * B[j0+j,k0+k]
                end
                @inbounds C[i0+i,k0+k] += c
            end
        end
    else
        m2 = m ÷ 2; n2 = n ÷ 2; p2 = p ÷ 2
        add_matmul_rec!(m2, n2, p2, i0, j0, k0, C, A, B)

        add_matmul_rec!(m-m2, n2, p2, i0+m2, j0, k0, C, A, B)
        add_matmul_rec!(m2, n-n2, p2, i0, j0+n2, k0, C, A, B)
        add_matmul_rec!(m2, n2, p-p2, i0, j0, k0+p2, C, A, B)

        add_matmul_rec!(m-m2, n-n2, p2, i0+m2, j0+n2, k0, C, A, B)
        add_matmul_rec!(m2, n-n2, p-p2, i0, j0+n2, k0+p2, C, A, B)
        add_matmul_rec!(m-m2, n2, p-p2, i0+m2, j0, k0+p2, C, A, B)

        add_matmul_rec!(m-m2, n-n2, p-p2, i0+m2, j0+n2, k0+p2, C, A, B)
    end
    return C
end

function matmul_rec!(C, A, B)
    m,n = size(A)
    n,p = size(B)
    size(C) == (m,p) || error("incorrect dimensions ", size(C), " ≠ $m × $p")
    fill!(C, 0)
    return add_matmul_rec!(m,n,p, 0,0,0, C,A,B)
end

matmul_rec(A, B) = matmul_rec!(Array{promote_type(eltype(A), eltype(B))}(undef,size(A,1), size(B,2)),A, B)

# Example:
# A = rand(50,60)
# B = rand(60,70)
# matmul_rec(A,B)

max_stars_repo_path: Code/Numerical Analysis & Machine Learning/Matrix Multiplication/CacheObliviousMatrixMultiplication.jl

max_stars_repo_name: BambooFlower/Math-Scripts

max_issues_repo_name: BambooFlower/Math-Scripts

max_forks_repo_name: BambooFlower/Math-Scripts


text:
"""
  IsingJ

  ver 0.2.2 Released Sep 2020
  This code is distributed under the constitution
  of MIT LICENSE.

  Log of ising.jl

  2019/04/24  Released by PearCandy                 """

module IsingJ
    include("./Params.jl")
    include("./MCsampling.jl")
    include("./Initial.jl")
    include("./Smoothing_mag.jl")
    using ArgParse
    using Dates
    using Printf
    using Random
    using .Params
    using .MCsampling
    using .Initial
    using .Smoothing_mag
    export main

    # main
    function main()
        parser= ArgParseSettings("Julia code for solving 2D Ising model") # parser
        @add_arg_table parser begin
            "--lower_temp", "-l"
            arg_type = Float64
            default  = 0.1
            help     = "lowest temperature"
            "--upper_temp", "-u"
            arg_type = Float64
            default  = 5.0
            help     = "highest temperature"
            "--step_num", "-n"
            arg_type = Int
            default  = 30
            help     = "step number of temperature"
            "--sweep_num", "-s"
            arg_type = Int
            default  = 2000
            help     = "montecarlo sweep number"
            "--warm_up", "-w"
            arg_type = Int
            default  = 2000
            help     = "warm up steps before sampling"
            "--exch_j", "-j"
            arg_type = Float64
            default  = 1.0
            help     = "exchange interaction"
            "--ext_field", "-e"
            arg_type = Float64
            default  = 0.0
            help     = "magnetic field"
            "--site_size", "-N"
            arg_type = Int
            default  = 64
            help     = "simulation size"
        end
            args = parse_args(parser)

    """ Arguments
        N_site    : number of site
        j0_val    : temperature
        low_t     : minimum temperature
        high_t    : maximum temperature
        ext       : external field
        inum_t    : number of temperature steps
        sweep_num : sweep number
        warm_up   : number of "warming-up" iteration  """

        N_site    = args["site_size"]
        j0_val    = args["exch_j"]
        low_t     = args["lower_temp"]
        up_t      = args["upper_temp"]
        ext_val   = args["ext_field"]
        inum_t    = args["step_num"]
        sweep_num = args["sweep_num"]
        warm_up   = args["warm_up"]
        cnst      = 2.26918  # Tc for square lattice

        # set variables
        site_size=64
        exch_j=1.0
        lower_temp=0.1
        upper_temp=5.0
        ext_field=0.0
        step_num=30
        sweep_num=2000
        warm_up=2000

        params = Params.Variables(site_size,exch_j,lower_temp,upper_temp,
                                  ext_field,step_num,sweep_num,warm_up)
        N_site    = params.site_size
        j0_val    = params.exch_j
        low_t     = params.lower_temp
        up_t      = params.upper_temp
        ext_val   = params.ext_field
        inum_t    = params.step_num
        sweep_num = params.sweep_num
        warm_up   = params.warm_up

        # initialization of the physical quantities
        spin      = zeros(N_site,N_site)
        mag       = 0
        mag_ave   = 0
        dt        = (up_t-low_t)/float(inum_t)
        _mag      = []

        Random.seed!(100)
        spin_ini = cold_start(N_site)  # Initialization, random_start(N_site)
        td=Dates.today()               # start motecalro simulation
        println(" ")
        println(td)
        println(" ")
        println("     I s i n g J")
        println("     Version 0.2.0, Released July 2019")
        println(" ")
        println("     Copyright (C) PearCandy")
        println(" ")
        println("     Usage: julia Ising.jl -h")
        println("---------------------------------------")
        open("log.dat","w") do io
            println(io,"#temperature magnetization energy Cv Chi Tc")
            for i in 1:inum_t       # main routine (montecalro sampling)
                _temp = low_t + dt*float(i)
                spin,mag,eng,cv,chi = mcsampling(spin_ini,sweep_num,warm_up,N_site,j0_val,ext_val,_temp)
                @printf(io, " %+06f"  ,_temp)
                @printf(io, " %+06f"  , mag)
                @printf(io, " %+06f"  , eng)
                @printf(io, " %+06f"  , cv)
                @printf(io, " %+06f"  , chi)
                @printf(io, " %+06f\n", cnst)
                append!(_mag,abs(mag))

                if i%10==0          # log
                    println(i,"/",inum_t,"steps")
                else
                end
            end
        end

        output_T_vs_mag(low_t,dt,_mag)   # smoothing T vs magnetization

        #-- plot (using python) --
        #    pyplot=pyimport("matplotlib.pyplot")
        #    numpy=pyimport("numpy")
        #    x=[i for i in 1:N_site]
        #    y=[i for i in 1:N_site]
        #    xx,yy=numpy.meshgrid(x,y)
        #    pyplot.pcolormesh(xx,yy,spin,cmap="bwr")
        #    pyplot.axes().set_aspect("equal")
        #    pyplot.show()

    end
end # module
#=-----------------------------------------------------------=#
#if occursin(PROGRAM_FILE,@__FILE__)
#    using .IsingJ
#    main()
#end


text:
using Random, JuMP, JuDGE, Test

include("solvers/setup_gurobi.jl")
#include("solvers/setup_cplex.jl")
#include("solvers/setup_coin.jl")
#include("solvers/setup_glpk.jl")

# this test is a test based off of a presentation by Prof. Andy Philpott
function knapsack_fixed()
   mytree = narytree(2,2)
   function invest_cost(node)
      if node == get_node(mytree,[1])
         180.0
      elseif node == get_node(mytree,[1,1])
         50.0
      elseif node == get_node(mytree,[1,2])
         60.0
      elseif node == get_node(mytree,[1,1,1])
         40.0
      elseif node == get_node(mytree,[1,1,2])
         60.0
      elseif node == get_node(mytree,[1,2,1])
         10.0
      elseif node == get_node(mytree,[1,2,2])
         10.0
      end
   end

   function item_volume(node)
      if node == get_node(mytree,[1])
         [6, 2, 1, 1, 1]
      elseif node == get_node(mytree,[1,1])
         [8, 2, 2, 2, 1]
      elseif node == get_node(mytree,[1,2])
         [8, 1, 1, 1, 3]
      elseif node == get_node(mytree,[1,1,1])
         [4, 4, 3, 1, 2]
      elseif node == get_node(mytree,[1,1,2])
         [1, 3, 1, 1, 2]
      elseif node == get_node(mytree,[1,2,1])
         [7, 3, 1, 1, 1]
      elseif node == get_node(mytree,[1,2,2])
         [2, 5, 2, 1, 2]
      end
   end

   function item_reward(node)
      if node == get_node(mytree,[1])
         [60, 20, 10, 15, 10]
      elseif node == get_node(mytree,[1,1])
         [8, 10, 20, 20, 10]
      elseif node == get_node(mytree,[1,2])
         [8, 10, 15, 10, 30]
      elseif node == get_node(mytree,[1,1,1])
         [40, 40, 35, 10, 20]
      elseif node == get_node(mytree,[1,1,2])
         [15, 35, 15, 15, 20]
      elseif node == get_node(mytree,[1,2,1])
         [70, 30, 15, 15, 10]
      elseif node == get_node(mytree,[1,2,2])
         [25, 50, 25, 15, 20]
      end
   end

   ### with judge
   function sub_problems(node)
      model = Model(JuDGE_SP_Solver)
      @expansion(model, bag, Bin)
      @capitalcosts(model, bag*invest_cost(node))
      @variable(model, y[1:5], Bin)
      @constraint(model, BagExtension, sum(y[i]*item_volume(node)[i] for i in 1:5) <= 3 + 4 * bag)
      @objective(model, Min, sum(-item_reward(node)[i] * y[i] for i in 1:5))
      return model
   end

   judy = JuDGEModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_MP_Solver)
   JuDGE.solve(judy)

   println("Objective: "*string(JuDGE.get_objval(judy)))
   JuDGE.print_expansions(judy,onlynonzero=false)

   println("Re-solved Objective: " * string(resolve_subproblems(judy)))

   deteq = DetEqModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_DE_Solver)
   JuDGE.solve(deteq)
   println("Deterministic Equivalent Objective: " * string(JuDGE.get_objval(deteq)))

   return objective_value(judy.master_problem)
end

# this test is a test based off of the test in the old JuDGE
function knapsack_random()
   Random.seed!(1)
   # how many investments?
   numinvest = 2;

   # number of items to pick from in the knapsack?
   numitems = 20

   # size of tree?
   degree = 3
   height = 3

   totalnodes = Int64((degree^(height+1) - 1)/(degree-1))

   investcost = zeros(totalnodes,numinvest)
   for i = 1:totalnodes
      investcost[i,:] = (rand(numinvest)*2  + 2*[2.0,3.5])*(1-((i-1)/(totalnodes*1.2)))
   end

   investvol = [40,50]
   initialcap = 80

   itemvolume = zeros(totalnodes,numitems)
   for i = 1:totalnodes
      itemvolume[i,:] = ((rand(numitems) .- 0.5)*2)*2 + collect(range(4,22,length = numitems))
   end

   itemcost = zeros(totalnodes,numitems)
   for i = 1:totalnodes
      itemcost[i,:] = ((rand(numitems) .- 0.5)*2)*0.5 + collect(range(0.5,1,length = numitems))
   end

   mytree = narytree(height,degree)

   nodes = collect(mytree)
   function data(node, input)
      input[findall(x -> x == node, nodes)[1], :]
   end

   function sub_problems(node)
      model = Model(JuDGE_SP_Solver)
      @expansion(model, bag[1:numinvest], Bin)
      @capitalcosts(model, sum(data(node,investcost)[i] * bag[i] for i in  1:numinvest))
      @variable(model, y[1:numitems], Bin)
      @constraint(model, BagExtension ,sum( y[i]*data(node,itemvolume)[i] for i in 1:numitems) <= initialcap + sum(bag[i]*investvol[i] for i in 1:numinvest))
      @objective(model, Min, sum(-data(node,itemcost)[i] * y[i] for i in 1:numitems))
      return model
   end

   function format_output(s::Symbol,values)
      if s==:bag
         return sum(values[i]*investvol[i] for i in 1:numinvest)
      end
      return nothing
   end

   judy = JuDGEModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_MP_Solver)
   JuDGE.solve(judy)

   println("Objective: "*string(JuDGE.get_objval(judy)))
   JuDGE.print_expansions(judy,format=format_output)

   println("Re-solved Objective: " * string(resolve_subproblems(judy)))

   deteq = DetEqModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_DE_Solver)
   JuDGE.solve(deteq)
   println("Deterministic Equivalent Objective: " * string(JuDGE.get_objval(deteq)))

   return objective_value(judy.master_problem)
end

function knapsack_branch_and_price()
   Random.seed!(50)
   # how many investments?
   numinvest = 5;

   # number of items to pick from in the knapsack?
   numitems = 10

   # size of tree?
   degree = 3
   height = 4

   totalnodes = Int64((degree^(height+1) - 1)/(degree-1))

   investcost = zeros(totalnodes,numinvest)
   for i = 1:totalnodes
      investcost[i,:] = ([1,1.8,3.5,6.8,13.5])*(1-((i-1)/(totalnodes*1.2)))
   end

   investvol = [1,2,4,8,16]
   initialcap = 0

   itemvolume = zeros(totalnodes,numitems)
   for i = 1:totalnodes
      itemvolume[i,:] = ((rand(numitems))*2) + collect(range(4,22,length = numitems))
   end

   itemcost = zeros(totalnodes,numitems)
   for i = 1:totalnodes
      itemcost[i,:] = ((rand(numitems) .- 0.5)*2)*2# + collect(range(0.5,1,length = numitems))
   end

   mytree = narytree(height,degree)

   nodes = collect(mytree)
   function data(node, input)
      input[findall(x -> x == node, nodes)[1], :]
   end

   function sub_problems(node)
      model = Model(JuDGE_SP_Solver)
      @expansion(model, bag[1:numinvest], Bin)
      @capitalcosts(model, sum(data(node,investcost)[i] * bag[i] for i in  1:numinvest))
      @variable(model, y[1:numitems], Bin)
      @constraint(model, BagExtension ,sum( y[i]*data(node,itemvolume)[i] for i in 1:numitems) <= initialcap + sum(bag[i]*investvol[i] for i in 1:numinvest))
      @objective(model, Min, sum(-data(node,itemcost)[i] * y[i] for i in 1:numitems))
      return model
   end

   function format_output(s::Symbol,values)
      if s==:bag
         return sum(values[i]*investvol[i] for i in 1:numinvest)
      end
      return nothing
   end

   judy = JuDGEModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_MP_Solver)

   best=JuDGE.branch_and_price(judy,termination=Termination(rlx_abstol=10^-7,inttol=10^-6))

   println("Objective: "*string(JuDGE.get_objval(best)))
   JuDGE.print_expansions(best,format=format_output)

   deteq = DetEqModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_DE_Solver)
   JuDGE.solve(deteq)
   println("Deterministic Equivalent Objective: " * string(JuDGE.get_objval(deteq)))

   return objective_value(best.master_problem)
end

function knapsack_risk_averse()
   Random.seed!(50)
   # how many investments?
   numinvest = 5;

   # number of items to pick from in the knapsack?
   numitems = 10

   # size of tree?
   degree = 3
   height = 4

   totalnodes = Int64((degree^(height+1) - 1)/(degree-1))

   investcost = zeros(totalnodes,numinvest)
   for i = 1:totalnodes
      investcost[i,:] = ([1,1.8,3.5,6.8,13.5])*(1-((i-1)/(totalnodes*1.2)))
   end

   investvol = [1,2,4,8,16]
   initialcap = 0

   itemvolume = zeros(totalnodes,numitems)
   for i = 1:totalnodes
      itemvolume[i,:] = ((rand(numitems))*2) + collect(range(4,22,length = numitems))
   end

   itemcost = zeros(totalnodes,numitems)
   for i = 1:totalnodes
      itemcost[i,:] = ((rand(numitems) .- 0.5)*2)*2# + collect(range(0.5,1,length = numitems))
   end

   mytree = narytree(height,degree)

   nodes = collect(mytree)
   function data(node, input)
      input[findall(x -> x == node, nodes)[1], :]
   end

   function sub_problems(node)
      model = Model(JuDGE_SP_Solver)
      @expansion(model, bag[1:numinvest], Bin)
      @capitalcosts(model, sum(data(node,investcost)[i] * bag[i] for i in  1:numinvest))
      @variable(model, y[1:numitems], Bin)
      @constraint(model, BagExtension ,sum( y[i]*data(node,itemvolume)[i] for i in 1:numitems) <= initialcap + sum(bag[i]*investvol[i] for i in 1:numinvest))
      @objective(model, Min, sum(-data(node,itemcost)[i] * y[i] for i in 1:numitems))
      return model
   end

   function format_output(s::Symbol,values)
      if s==:bag
         return sum(values[i]*investvol[i] for i in 1:numinvest)
      end
      return nothing
   end

   judy = JuDGEModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_MP_Solver, risk=Risk(0.5,0.05))

   best=JuDGE.branch_and_price(judy,termination=Termination(rlx_abstol=10^-6,inttol=10^-6))
   println("Objective: "*string(JuDGE.get_objval(best)))
   JuDGE.print_expansions(best,format=format_output)

   println("Re-solved Objective: " * string(resolve_subproblems(best)))

   deteq = DetEqModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_DE_Solver, risk=Risk(0.5,0.05))
   JuDGE.solve(deteq)
   println("Deterministic Equivalent Objective: " * string(JuDGE.get_objval(deteq)))

   return objective_value(best.master_problem)
end

function knapsack_delayed_investment(;CVaR=RiskNeutral())
   Random.seed!(100)
   # how many investments?
   numinvest = 2;

   # number of items to pick from in the knapsack?
   numitems = 25

   # size of tree?
   degree = 3
   height = 3

   totalnodes = Int64((degree^(height+1) - 1)/(degree-1))

   investcost = zeros(totalnodes,numinvest)
   for i = 1:totalnodes
      investcost[i,:] = (rand(numinvest)*2  + 2*[2.0,3.5])*(1-((i)/(totalnodes)))
   end

   # investvol = [40,45,50,70]
   investvol = [40,50]
   initialcap = 80

   itemvolume = zeros(totalnodes,numitems)
   for i = 1:totalnodes
      itemvolume[i,:] = ((rand(numitems) .- 0.5)*2)*2 + collect(range(4,22,length = numitems))
   end

   itemcost = zeros(totalnodes,numitems)
   for i = 1:totalnodes
      itemcost[i,:] = ((rand(numitems) .- 0.5)*2)*0.5 + collect(range(0.5,1,length = numitems))
   end

   mytree = narytree(height,degree)

   nodes = collect(mytree)
   function data(node, input)
      input[findall(x -> x == node, nodes)[1], :]
   end

   function sub_problems(node)
      model = Model(JuDGE_SP_Solver)
      @expansion(model, bag[1:numinvest], Bin, lag=1)
      @capitalcosts(model, sum(data(node,investcost)[i] * bag[i] for i in  1:numinvest))
      @variable(model, y[1:numitems], Bin)
      @constraint(model, BagExtension ,sum( y[i]*data(node,itemvolume)[i] for i in 1:numitems) <= initialcap + sum(bag[i]*investvol[i] for i in 1:numinvest))
      @objective(model, Min, sum(-data(node,itemcost)[i] * y[i] for i in 1:numitems))
      return model
   end

   function format_output(s::Symbol,values)
      if s==:bag_bought
         return sum(values[i]*investvol[i] for i in 1:numinvest)
      elseif s==:bag_arrived
         return 0.0
      end
      return nothing
   end

   judy = JuDGEModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_MP_Solver, risk=CVaR)

   judy=JuDGE.branch_and_price(judy,termination=Termination(rlx_abstol=10^-6,inttol=10^-6))
   println("Objective: "*string(JuDGE.get_objval(judy, risk=CVaR)))
   JuDGE.print_expansions(judy,format=format_output)

   println("Re-solved Objective: " * string(resolve_subproblems(judy)))

   deteq = DetEqModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_DE_Solver, risk=CVaR)
   JuDGE.solve(deteq)
   println("Deterministic Equivalent Objective: " * string(JuDGE.get_objval(deteq, risk=CVaR)))
   return objective_value(judy.master_problem)
end

function knapsack_shutdown()
   mytree = narytree(2,2)

   divest_revenue=Dict{AbstractTree,Float64}()
   divest_revenue[get_node(mytree,[1])]=60.0
   divest_revenue[get_node(mytree,[1,1])]=58.0
   divest_revenue[get_node(mytree,[1,2])]=35.0
   divest_revenue[get_node(mytree,[1,1,1])]=40.0
   divest_revenue[get_node(mytree,[1,1,2])]=35.0
   divest_revenue[get_node(mytree,[1,2,1])]=25.0
   divest_revenue[get_node(mytree,[1,2,2])]=15.0

   item_volume=Dict{AbstractTree,Array{Float64,1}}()
   item_volume[get_node(mytree,[1])]=[6, 2, 1, 1, 1]
   item_volume[get_node(mytree,[1,1])]=[8, 2, 2, 2, 1]
   item_volume[get_node(mytree,[1,2])]=[8, 1, 1, 1, 3]
   item_volume[get_node(mytree,[1,1,1])]=[4, 4, 3, 1, 2]
   item_volume[get_node(mytree,[1,1,2])]=[1, 3, 1, 1, 2]
   item_volume[get_node(mytree,[1,2,1])]=[7, 3, 1, 1, 1]
   item_volume[get_node(mytree,[1,2,2])]=[2, 5, 2, 1, 2]

   item_reward=Dict{AbstractTree,Array{Float64,1}}()
   item_reward[get_node(mytree,[1])]=[60, 20, 10, 15, 10]
   item_reward[get_node(mytree,[1,1])]=[8, 10, 20, 20, 10]
   item_reward[get_node(mytree,[1,2])]=[8, 10, 15, 10, 30]
   item_reward[get_node(mytree,[1,1,1])]=[40, 40, 35, 10, 20]
   item_reward[get_node(mytree,[1,1,2])]=[15, 35, 15, 15, 20]
   item_reward[get_node(mytree,[1,2,1])]=[70, 30, 15, 15, 10]
   item_reward[get_node(mytree,[1,2,2])]=[25, 50, 25, 15, 20]

   ### with judge
   function sub_problems(node)
      model = Model(JuDGE_SP_Solver)
      @shutdown(model, bag, Bin)
      @capitalcosts(model, -bag*divest_revenue[node])
      @variable(model, y[1:5], Bin)
      @constraint(model, BagExtension, sum(y[i]*item_volume[node][i] for i in 1:5) <= 4 - 2 * bag)
      @objective(model, Min, sum(-item_reward[node][i] * y[i] for i in 1:5))
      return model
   end

   function format_output(s::Symbol,value)
      if s==:bag
         return -2.0*value
      end
      return nothing
   end

   judy = JuDGEModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_MP_Solver)
   JuDGE.solve(judy)

   println("Objective: "*string(JuDGE.get_objval(judy)))
   JuDGE.print_expansions(judy,format=format_output)

   println("Re-solved Objective: " * string(resolve_subproblems(judy)))

   deteq = DetEqModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_DE_Solver)
   JuDGE.solve(deteq)
   println("Deterministic Equivalent Objective: " * string(JuDGE.get_objval(deteq)))

   return objective_value(judy.master_problem)
end

function knapsack_budget()
   mytree = narytree(2,2)

   invest_cost = Dict( zip( collect(mytree,order=:breadth), [15, 8, 8, 4, 4, 4, 4]) )

   item_volume = Dict( zip( collect(mytree,order=:breadth), [ [4, 3, 3, 1, 2],
                                                            [5, 3, 4, 2, 1],
                                                            [5, 4, 2, 7, 2],
                                                            [5, 4, 1, 8, 2],
                                                            [3, 1, 5, 6, 3],
                                                            [2, 5, 8, 4, 6],
                                                            [7, 5, 4, 2, 3] ]) )

   item_reward = Dict( zip( collect(mytree,order=:breadth), [ [32, 9, 9, 4, 8],
                                                            [30, 12, 40, 10, 9],
                                                            [20, 28, 12, 42, 12],
                                                            [40, 28, 9, 24, 10],
                                                            [15, 7, 20, 48, 12],
                                                             [10, 30, 54, 32, 30],
                                                            [32, 25, 24, 14, 24] ]) )

   num_items=5
   num_invest=6
   initial_volume = 6
   invest_volume = [2,2,2,3,3,3]

   ### with judge
   function sub_problems(node)
      sp = Model(JuDGE_SP_Solver)
      @expansion(sp, invest[1:num_invest], Bin)
      @capitalcosts(sp, sum(invest[i]*invest_volume[i] for i=1:num_invest)*invest_cost[node])
      @variable(sp, y[1:num_items], Bin)
      @constraint(sp, BagExtension, sum(y[i]*item_volume[node][i] for i in 1:num_items) <= initial_volume + sum(invest_volume[i] * invest[i] for i in 1:num_invest))
      @objective(sp, Min, sum(-item_reward[node][i] * y[i] for i in 1:num_items))
      return sp
   end

   function format_output(s::Symbol,value)
      if s==:invest
         return sum(invest_volume[i]*value[i] for i in 1:num_invest)
      end
      return nothing
   end

   function budget(model,tree)
      for node in collect(tree)
         @constraint(model,sum(invest_cost[node]*invest_volume[i]*invest[node][i] for i in 1:num_invest)<=40)
      end
   end

   judy = JuDGEModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_MP_Solver, sideconstraints=budget)
   #JuDGE.solve(judy)
   judy = JuDGE.branch_and_price(judy)
   println("Objective: "*string(JuDGE.get_objval(judy)))
   JuDGE.print_expansions(judy, format=format_output)

   println("Re-solved Objective: " * string(resolve_subproblems(judy)))

   deteq = DetEqModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_DE_Solver, sideconstraints=budget)
   JuDGE.solve(deteq)
   println("Deterministic Equivalent Objective: " * string(JuDGE.get_objval(deteq)))

   return objective_value(judy.master_problem)
end

@test knapsack_fixed() ≈ -131.25 atol = 1e-3
@test knapsack_random() ≈ -34.749 atol = 1e-3
@test knapsack_branch_and_price() ≈ -0.69456 atol = 1e-4
@test knapsack_risk_averse() ≈ -0.27292 atol = 1e-4
@test knapsack_delayed_investment() ≈ -34.058 atol = 1e-3
@test knapsack_delayed_investment(CVaR=Risk(0.95,0.05)) ≈ -31.344 atol = 1e-3
@test knapsack_shutdown() ≈ -145.25 atol = 1e-3
@test knapsack_budget() ≈ -159.0 atol = 1e-3


text:
# Here I include some maps between the experimental
# data ids and the model ids

## ------------------------------------------------------------------
# maps between Kayser2005 https://doi.org/10.1099/mic.0.27481-0.
# data and the model
function load_mets_map()
    mets_map = Dict()
    mets_map["GLC"] = "glc_DASH_D_e"
    mets_map["SUCC"] = "succ_e"
    mets_map["FORM"] = "for_e"
    mets_map["THM"] = "thm_e"
    mets_map["NH4"] = "nh4_e"
    mets_map["CIT"] = "cit_e"
    mets_map["CO2"] = "co2_e"
    mets_map["O2"] = "o2_e"
    mets_map["AC"] = "ac_e"
    mets_map["PYR"] = "pyr_e"
    mets_map["LAC"] = "lac_DASH_D_e"
    mets_map["MAL"] = "mal_DASH_D_e";
    for (k, v) in mets_map
        mets_map[v] = k
    end
    return mets_map
end

# EX_glc_LPAREN_e_RPAREN_
function load_rxns_map()
    rxns_map = Dict()
    rxns_map["D"] = "BiomassEcoli"
    rxns_map["AC"] = "EX_ac_LPAREN_e_RPAREN_"
    rxns_map["NH4"] = "EX_nh4_LPAREN_e_RPAREN_"
    rxns_map["GLC"] = "EX_glc_LPAREN_e_RPAREN_"
    rxns_map["THM"] = "EX_thm_LPAREN_e_RPAREN_"
    rxns_map["CO2"] = "EX_co2_LPAREN_e_RPAREN_"
    rxns_map["FORM"] = "EX_for_LPAREN_e_RPAREN_"
    rxns_map["CIT"] = "EX_cit_LPAREN_e_RPAREN_"
    rxns_map["SUCC"] = "EX_succ_LPAREN_e_RPAREN_"
    rxns_map["O2"] = "EX_o2_LPAREN_e_RPAREN_"
    rxns_map["MAL"] = "EX_mal_L_LPAREN_e_RPAREN_"
    rxns_map["LAC"] = "EX_lac_D_LPAREN_e_RPAREN_"
    rxns_map["PYR"] = "EX_pyr_LPAREN_e_RPAREN_"
    for (k, v) in rxns_map
        rxns_map[v] = k
    end
    return rxns_map
end

## ------------------------------------------------------------------
# enzymatic costs
# from Beg, (2007) https://doi.org/10.1073/pnas.0609845104.

# A map between model ids and the reactions reported in Beg2007
function load_beg_rxns_map()
    beg_rxns_map = Dict(
        "carbonyl reductase (NADPH)" => ["P5CR"],
        "alcohol dehydrogenase (NADP+)" => ["ALCD19"],
        "quinate/shikimate dehydrogenase" => ["SHK3Dr"],
        "malate dehydrogenase (decarboxylating)" => ["MDH","MDH2", "MDH3"],
        "3alpha-hydroxysteroid dehydrogenase (B-specific)" => [""],
        "2-hydroxy-3-oxopropionate reductase" => [""],
        "glucose dehydrogenase (acceptor)" => ["G6PDH2r", "UDPGD"],
        "cellobiose dehydrogenase (acceptor)" => [""],
        "peroxidase" => [""],
        "catechol 2,3-dioxygenase" => [""],
        "arachidonate 8-lipoxygenase" => [""],
        "calcidiol 1-monooxygenase" => [""],
        "nitric-oxide synthase" => [""],
        "phenylalanine 4-monooxygenase" => [""],
        "tryptophan 5-monooxygenase" => [""],
        "Carboxylate reductase" => ["P5CR"],
        "arsenate reductase (donor)" => [""],
        "biliverdin reductase" => [""],
        "15-oxoprostaglandin 13-oxidase" => [""],
        "coproporphyrinogen oxidase" => ["CPPPGO","PPPGO"],
        "long-chain-acyl-CoA dehydrogenase" => [""],
        "butyryl-CoA dehydrogenase" => [""],
        "acyl-CoA dehydrogenase" => [""],
        "L-amino-acid oxidase" => ["ASPO3","ASPO4","ASPO5","ASPO6"],
        "amine oxidase (flavin-containing)" => ["PYAM5PO"],
        "methylenetetrahydrofolate reductase [NAD(P)H]" => ["MTHFR2"],
        "formyltetrahydrofolate dehydrogenase" => ["MTHFD"],
        "sarcosine oxidase" => [""],
        "nitrate reductase (NADH)" => [""],
        "nitrite reductase (NO-forming)" => [""],
        "nitrate reductase" => ["NO3R2"],
        "trypanothione-disulfide reductase" => [""],
        "glutathione-disulfide reductase" => [""],
        "thioredoxin-disulfide reductase" => [""],
        "thiol oxidase" => [""],
        "nitrate reductase (cytochrome)" => [""],
        "aspartate carbamoyltransferase" => ["ASPCT"],
        "serine O-acetyltransferase" => ["SERAT"],
        "protein-glutamine gamma-glutamyltransferase" => [""],
        "gamma-glutamyltransferase" => ["CRNBTCT"],
        "citrate (Si)-synthase" => ["CS"],
        "kaempferol 3-O-galactosyltransferase" => [""],
        "NAD+ ADP-ribosyltransferase" => ["NNDMBRT"],
        "di-trans,poly-cis-decaprenylcistransferase" => ["ACGAMT"],
        "cystathionine gamma-synthase" => [""],
        "adenosine kinase" => ["ADNK1"],
        "glycerate kinase" => ["GLYCK"],
        "galactokinase" => ["GALKr"],
        "[pyruvate dehydrogenase (acetyl-transferring)] kinase" => [""],
        "guanylate kinase" => ["GK1"],
        "FMN adenylyltransferase" => ["FMNAT"],
        "tRNA adenylyltransferase" => [""],
        "aryl sulfotransferase" => [""],
        "aminoacyl-tRNA hydrolase" => [""],
        "carboxymethylenebutenolidase" => [""],
        "ubiquitin thiolesterase" => [""],
        "fructose-bisphosphatase" => ["FBP","FBA"],
        "[phosphorylase] phosphatase" => [""],
        "phosphoglycolate phosphatase" => ["PGLYCP"],
        "protein-tyrosine-phosphatase" => [""],
        "inositol-polyphosphate 5-phosphatase" => ["MI1PP"],
        "3',5'-cyclic-GMP phosphodiesterase" => [""],
        "beta-glucosidase" => ["MLTG1","MLTG2","MLTG3","MLTG4","MLTG5"],
        "beta-glucuronidase" => [""],
        "glucosylceramidase" => [""],
        "cyclomaltodextrinase" => [""],
        "alpha-N-arabinofuranosidase" => [""],
        "purine nucleosidase" => ["AMPN","AHCYSNS","CMPN","MTAN","NMNN"],
        "rRNA N-glycosylase" => [""],
        "NAD+ nucleosidase" => [""],
        "Xaa-Pro aminopeptidase" => [""],
        "dipeptidyl-peptidase I" => [""],
        "peptidyl-dipeptidase A" => [""],
        "coagulation factor Xa" => [""],
        "t-Plasminogen activator" => [""],
        "cathepsin B" => [""],
        "envelysin" => [""],
        "amidase" => ["","GSPMDA","NMNDA","NNAM"],
        "formamidase" => [""],
        "arginase" => [""],
        "guanidinoacetase" => ["GUAD"],
        "apyrase" => [""],
        "phloretin hydrolase" => [""],
        "Orotidine-5'-phosphate decarboxylase" => ["OMPDC"],
        "4-Hydroxybenzoate decarboxylase" => ["OPHBDC"],
        "Threonine aldolase" => ["THRAr"],
        "enoyl-CoA hydratase" => [""],
        "Uroporphyrinogen-III synthase" => ["UPP3S"],
        "dihydroxy-acid dehydratase" => ["DHAD1","DHAD2"],
        "pectin lyase" => [""],
        "DNA-(apurinic or apyrimidinic site) lyase" => [""],
        "lactoylglutathione lyase" => ["LGTHL"],
        "guanylate cyclase" => [""],
        "dTDP-4-dehydrorhamnose 3,5-epimerase" => ["TDPDRE"],
        "UDP-glucose 4-epimerase" => ["UDPG4E"],
        "Triose-phosphate isomerase" => ["TPI"],
        "steroid DELTA-isomerase" => [""],
        "dodecenoyl-CoA isomerase" => [""],
        "Glutamate-1-semialdehyde 2,1-aminomutase" => [""],
        "Chalcone isomerase" => [""],
        "Chloromuconate cycloisomerase" => [""],
        "Tyrosine-tRNA ligase" => [""],
        "Threonine-tRNA ligase" => [""],
        "Isoleucine-tRNA ligase" => [""],
        "Lysine-tRNA ligase" => [""],
        "formate-tetrahydrofolate ligase" => [""],
        "Adenylosuccinate synthase" => ["ADSS"],
        "DNA ligase (NAD+)" => [""]
    )

    # inverse relation
    for (beg_id, model_ids) in beg_rxns_map
        for model_id in model_ids
            beg_ids = get!(beg_rxns_map, model_id, [])
            push!(beg_ids, beg_id)
        end
    end
    return beg_rxns_map
end

## ------------------------------------------------------------------
# base model exch met map
# A quick way to get exchages from mets and te other way around
load_exch_met_map() = lprocdir(iJR904, "exch_met_map.bson"; verbose = false)

## ------------------------------------------------------------------
# The intakes bounds of the network are determined by the
# medium concentration in the Chemostat model (see Cossios paper)
# This is a base medium for modeling

function load_base_intake_info()
    return Dict(
        "EX_glc_LPAREN_e_RPAREN_" => Dict("c"=> maximum(Fd.val(:cGLC)), "lb"=> -ABS_MAX_BOUND),

        "EX_nh4_LPAREN_e_RPAREN_" => Dict("c"=> MAX_CONC, "lb"=> -ABS_MAX_BOUND),
        "EX_o2_LPAREN_e_RPAREN_"  => Dict("c"=> MAX_CONC, "lb"=> -ABS_MAX_BOUND),
        "EX_pi_LPAREN_e_RPAREN_"  => Dict("c"=> MAX_CONC, "lb"=> -ABS_MAX_BOUND),
        "EX_so4_LPAREN_e_RPAREN_" => Dict("c"=> MAX_CONC, "lb"=> -ABS_MAX_BOUND),
    )
end

function intake_info(exp)
    return Dict(
        "EX_glc_LPAREN_e_RPAREN_" => Dict("c"=> Fd.val(:cGLC, exp), "lb"=> -ABS_MAX_BOUND),

        "EX_nh4_LPAREN_e_RPAREN_" => Dict("c"=> MAX_CONC, "lb"=> -ABS_MAX_BOUND),
        "EX_o2_LPAREN_e_RPAREN_"  => Dict("c"=> MAX_CONC, "lb"=> -ABS_MAX_BOUND),
        "EX_pi_LPAREN_e_RPAREN_"  => Dict("c"=> MAX_CONC, "lb"=> -ABS_MAX_BOUND),
        "EX_so4_LPAREN_e_RPAREN_" => Dict("c"=> MAX_CONC, "lb"=> -ABS_MAX_BOUND),
    )
end

function load_krebs_iders()
    krebs_iders = ["SUCD1", "SUCOAS", "AKGDH", "ICDHyr",
        "ACONT", "CS", "MDH", "FUM", "MALS", "ICL"
    ]
end

function load_inners_idermap()
    kreps_idermap = Dict(
        "SUCD1" => ["SUCD1i"],
        "SUCOAS" => ["SUCOAS_fwd", "SUCOAS_bkwd"],
        "AKGDH" => ["AKGDH"],
        "ICDHyr" => ["ICDHyr_fwd", "ICDHyr_bkwd"],
        "ACONT" => ["ACONT_bkwd", "ACONT_fwd"],
        "CS" => ["CS"],
        "MDH" => ["MDH_fwd", "MDH_bkwd"],
        "FUM" => ["FUM_fwd", "FUM_bkwd"],
        "MALS" => ["MALS"],
        "ICL" => ["ICL"]
    )
end


text:
#############################################################################
#############################################################################
#
# This file contains units tests for polynomial operations
#
#############################################################################
#############################################################################

"""
Test product of polynomials.
"""
function prod_test_poly(;N::Int = 10, N_prods::Int = 20, seed::Int = 0)
    Random.seed!(seed)
    for _ in 1:N
        p1 = rand(Polynomial)
        p2 = rand(Polynomial)
        prod = p1*p2
        @assert leading(prod) == leading(p1)*leading(p2)
    end

    for _ in 1:N
        p_base = Polynomial(Term(1,0))
        for _ in 1:N_prods
            p = rand(Polynomial)
            prod = p_base*p
            @assert leading(prod) == leading(p_base)*leading(p)
            p_base = prod
        end
    end
    println("prod_test_poly - PASSED")
end

"""
Test product of polynomials modulo P.
"""
function prod_test_poly_mod_p(;N::Int = 10, N_prods::Int = 20, seed::Int = 0)
    Random.seed!(seed)
    sample_primes_array = [3,5,7,11,13,17,19]
    for _ in 1:N
        rand_prime = rand(sample_primes_array)
        p1 = rand(Polynomial)
        p2 = rand(Polynomial)
        p1mod = PolynomialModP(p1, rand_prime)
        p2mod = PolynomialModP(p2, rand_prime)
        prod = p1mod*p2mod
        @assert mod(leading(prod), rand_prime) == leading(mod(p1mod.terms*p2mod.terms, rand_prime))
    end

    for _ in 1:N
        rand_prime = rand(sample_primes_array)
        p_base = PolynomialModP(Polynomial(Term(1,0)), rand_prime)
        for _ in 1:N_prods
            p = PolynomialModP(rand(Polynomial), rand_prime)
            prod = p_base*p
            prod = PolynomialModP(prod, rand_prime)
            @assert leading(prod.terms) == leading(mod(p_base.terms*p.terms, rand_prime))
            p_base = prod
        end
    end
    println("prod_test_poly_mod_p - PASSED")
end

"""
Test power of polynomial/integer mod P
"""
function power_test_poly(;N::Int = 10, seed::Int = 0)
    sample_primes_array = [3,5,7,11,13,17,19]
    Random.seed!(seed)
    x = x_poly()
    print("Proceeding with power tests ~ [")
    for i in 1:N
        exponent = rand(10:20)
        rand_prime = rand(sample_primes_array)
        p1 = rand(Polynomial)
        out = 1
        for i in 1:exponent
            out = mod(out*p1, rand_prime)
        end
        if ^(p1, exponent, rand_prime) - out != 0
            println("--- TEST FAILED ---")
            break
        end
        print("-")
    end
    print("]\n")
    println("power_test_poly - COMPLETED\n")
end

"""
Test power of polynomial mod P
"""
function power_test_poly_modp(;N::Int = 10, seed::Int = 0)
    sample_primes_array = [3,5,7,11,13,17,19]
    Random.seed!(seed)
    x = x_poly()
    print("Proceeding with power tests using PolynomialModP data type ~ [")
    for i in 1:N
        exponent = rand(10:20)
        rand_prime = rand(sample_primes_array)
        p1 = rand(Polynomial)
        p1mod = PolynomialModP(p1, rand_prime)
        out = 1
        for i in 1:exponent
            out = mod(out*p1mod.terms, p1mod.prime)
        end
        if ^(p1mod, exponent) - out != 0
            println("--- TEST FAILED ---")
            break
        end
        print("-")
    end
    print("]\n")
    println("power_test_poly_modp - COMPLETED\n")
end

"""
Benchmark multiplication
"""
function prod_benchmark(;N::Int = 50, seed::Int = 0)
    println("Begin prod_benchmark: Demonstrating performance of polynomial multiplication algorithm")
    Random.seed!(seed)
    p1 = 0
    p2 = 0
    x = x_poly()
    for i in 1:N
        p1 += rand(1:N)*x^(i)
        p2 += rand(1:N)*x^(i)
        if i%10 == 0
            print("Time taken to multiply two polynomials with $i terms:")
            @time prod = p1*p2
            @assert leading(prod) == leading(p1)*leading(p2)
        end
    end
    println("prod_benchmark COMPLETED\n")
end

function prod_modP_benchmark(;N::Int = 50, seed::Int = 0)
    println("Begin prod_modP_benchmark: Demonstrating performance of polynomial multiplication algorithm in modulo P")
    Random.seed!(seed)
    sample_primes_array = [3,5,7,11,13,17,19]
    p1 = 0
    p2 = 0
    x = x_poly()
    for i in 1:N
        rand_prime = rand(sample_primes_array)
        p1 += rand(1:N)*x^(i)
        p2 += rand(1:N)*x^(i)
        p1mod = PolynomialModP(p1, rand_prime)
        p2mod = PolynomialModP(p2, rand_prime)
        if i%10 == 0
            print("Time taken to multiply two polynomials with $i terms:")
            @time prod = p1mod*p2mod
            @assert leading(prod) == leading(mod((mod(p1mod.terms, rand_prime) * mod(p2mod.terms, rand_prime)), rand_prime))
        end
    end
    println("prod_modP_benchmark COMPLETED\n")
end

function crt_benchmark(;N::Int = 50, seed::Int = 0)
    println("Begin crt_benchmark: Demonstrating performance of polynomial multiplication using the Chinese Remainder Theorem")
    Random.seed!(seed)
    sample_primes_array = [3,5,7,11,13,17,19]
    p1 = 0
    p2 = 0
    x = x_poly()
    for i in 1:N
        x = x_poly()
        rand_prime = rand(sample_primes_array)
        p1 += rand(1:N)*x^(i)
        p2 += rand(1:N)*x^(i)
        if i%10 == 0
            print("Time taken to multiply two polynomials with $i terms:")
            @time prod = mult_poly_with_crt(p1, p2)
            @assert prod == p1 * p2 "Error: CRT does not net p1*p2"
        end
    end
    println("crt_benchmark COMPLETED\n")
end

"""
Test derivative of polynomials (as well as product).
"""
function prod_derivative_test_poly(;N::Int = 10^3,  seed::Int = 0)
    Random.seed!(seed)
    for _ in 1:N
        p1 = rand(Polynomial)
        p2 = rand(Polynomial)
        p1d = derivative(p1)
        p2d = derivative(p2)
        @assert (p1d*p2) + (p1*p2d) == derivative(p1*p2)
    end
    println("prod_derivative_test_poly - PASSED")
end

"""
Test derivative of polynomials modulo P (as well as product).
"""
function prod_derivative_test_poly_mod_p(;N::Int = 10^3,  seed::Int = 0)
    Random.seed!(seed)
    sample_primes_array = [3,5,7,11,13,17,19]
    for _ in 1:N
        rand_prime = rand(sample_primes_array)
        p1 = rand(Polynomial)
        p2 = rand(Polynomial)
        p1mod = PolynomialModP(p1, rand_prime)
        p2mod = PolynomialModP(p2, rand_prime)
        p1d = derivative(p1)
        p2d = derivative(p2)
        p1mod_d = PolynomialModP(p1d, rand_prime)
        p2mod_d = PolynomialModP(p2d, rand_prime)
        @assert mod((p1d*p2) + (p1*p2d), rand_prime) == mod(derivative(p1mod.terms*p2mod.terms), rand_prime)
    end
    println("prod_derivative_test_poly_mod_p - PASSED")
end

"""
Test division of polynomials modulo p.
"""
## REMOVE
function division_test_poly(;prime::Int = 101, N::Int = 10^4, seed::Int = 0)
    Random.seed!(seed)
    for _ in 1:N
        p1 = rand(Polynomial)
        p2 = rand(Polynomial)
        p_prod = p1*p2
        q, r = Polynomial(), Polynomial()
        try
            q, r = divide(p_prod, p2)(prime)
            if (q, r) == (nothing,nothing)
                println("Unlucky prime: $p1 is reduced to $(p1 % prime) modulo $prime")
                continue
            end
        catch e
            if typeof(e) == DivideError
                @assert mod(p2, prime) == 0
            else
                throw(e)
            end
        end
        @assert iszero( mod(q*p2+r - p_prod, prime) )
    end
    println("division_test_poly - PASSED")
end

"""
Test the extended euclid algorithm for polynomials modulo p.
"""
function ext_euclid_test_poly(;prime::Int=101, N::Int = 10^3, seed::Int = 0)
    Random.seed!(seed)
    for _ in 1:N
        p1 = rand(Polynomial)
        p2 = rand(Polynomial)
        g, s, t = extended_euclid_alg(p1, p2, prime)
        @assert mod(s*p1 + t*p2 - g, prime) == 0
    end
    println("ext_euclid_test_poly - PASSED")
end

"""
Test division of polynomials modulo p using PolynomialModP data type.
"""
function division_test_poly_mod_p(;prime::Int = 101, N::Int = 10^3, seed::Int = 0)
    Random.seed!(seed)
    for _ in 1:N
        p1 = rand(Polynomial)
        p2 = rand(Polynomial)
        p1mod = PolynomialModP(p1, prime)
        p2mod = PolynomialModP(p2, prime)
        p_prod = PolynomialModP(p1mod.terms*p2mod.terms, prime)
        q, r = Polynomial(), Polynomial()
        try
            q, r = divide(p_prod, p2mod)
            if (q, r) == (nothing,nothing)
                println("Unlucky prime: $p1 is reduced to $(p1 % prime) modulo $prime")
                continue
            end
        catch e
            if typeof(e) == DivideError
                @assert mod(p2mod.terms, prime) == 0
            else
                throw(e)
            end
        end
        @assert iszero(mod(q*p2mod.terms+r - p_prod.terms, prime))
    end
    println("division_test_poly_mod_p - PASSED")
end

"""
Test the extended euclid algorithm for polynomials modulo p.
"""
function ext_euclid_test_poly_mod_p(;prime::Int=101, N::Int = 10^3, seed::Int = 0)
    Random.seed!(seed)
    for _ in 1:N
        p1 = rand(Polynomial)
        p2 = rand(Polynomial)
        p1mod = PolynomialModP(p1, prime)
        p2mod = PolynomialModP(p2, prime)
        g, s, t = extended_euclid_alg(p1mod, p2mod)
        @assert mod(s*p1mod.terms + t*p2mod.terms - g, prime) == 0
    end
    println("ext_euclid_test_poly_mod_p - PASSED")
end