text:
import Random: shuffle

struct Node{T}
    value::NTuple{2,T}
    value_id::Int
    left::Int
    right::Int
end

struct KDTree{T}
    nodes::Vector{Node{T}}
end

function KDTree(values::Vector{NTuple{2,T}}) where {T}
    n = length(values)
    nodes = [Node((0.0, 0.0), 0, 0, 0) for _ in 1:n]
    node_id = 1
    nodes[node_id] = Node(values[node_id], node_id, 0, 0)
    tree = KDTree(nodes)
    next_node = node_id + 1
    for value_id in 2:n
        next_node = insert!(tree, values[value_id], value_id, node_id, next_node, 1)
    end
    return tree
end

function insert!(tree, value, value_id, node_id, next_node, dim)
    root = tree.nodes[node_id]
    if value[dim] > root.value[dim]
        if root.right == 0
            tree.nodes[next_node] = Node(value, value_id, 0, 0)
            tree.nodes[node_id] = Node(root.value, root.value_id, root.left, next_node)
            next_node += 1
        else
            next_node = insert!(tree, value, value_id, root.right, next_node, dim % 2 + 1)
        end
    else
        if root.left == 0
            tree.nodes[next_node] = Node(value, value_id, 0, 0)
            tree.nodes[node_id] = Node(root.value, root.value_id, next_node, root.right)
            next_node += 1
        else
            next_node = insert!(tree, value, value_id, root.left, next_node, dim % 2 + 1)
        end
    end
    return next_node
end

function spatial_order(tree)
    n = length(tree.nodes)
    order = zeros(Int, n)
    spatial_order!(order, tree, 1, 1)
    return order
end

function spatial_order!(order, tree, index, node_id)
    root = tree.nodes[node_id]
    if root.left != 0
        index = spatial_order!(order, tree, index, root.left)
    end
    order[root.value_id] = index
    index += 1
    if root.right != 0
        index = spatial_order!(order, tree, index, root.right)
    end
    return index
end

function biased_random_order(vertices)
    n = length(vertices)
    m = Int(ceil(log(2, n))) # Number of insertion rounds
    vertices = shuffle(vertices)
    order = spatial_order(KDTree(vertices))
    reverse = true
    start = 1
    for i in 1:m - 1
        count = Int(max(1, floor(n/2^(m + 1 - i))))
        stop = start + count - 1
        order_round = @view order[start:stop]
        perm = sortperm(order_round; rev=reverse)
        vertices[start:stop] .= @view vertices[start:stop][perm]
        reverse = !reverse
        start = stop + 1
    end
    order_round = @view order[start:n]
    perm = sortperm(order_round; rev=reverse)
    vertices[start:n] .= @view vertices[start:n][perm]
    return vertices
end

max_stars_repo_name: PaulMatlashewski/DelaunayTriangulations

max_issues_repo_name: PaulMatlashewski/DelaunayTriangulations

max_forks_repo_name: PaulMatlashewski/DelaunayTriangulations


text:
# script that outputs some statistics of the model fits

using DataFrames, Gadfly, Colors, MAT, CSV
import Cairo, Fontconfig

include("common.jl")
include("models.jl")
include("behaviorstats.jl")
include("simmodel.jl")
include("modelfitting.jl")

# settings (when manually evaluating the objective function)
const blockreps = 10
const simtrials = 10^5
const burnin = 1000
const nrats = 1 # do not divide by nrats = 4
const δt = 0.0001
const t = Theme(minor_label_font_size=6pt,
                major_label_font_size=6pt,
                grid_line_width=0pt, grid_color=RGB(1.0, 1.0, 1.0))
const fitpath = "fits"
const datapath = "../fitdata/$(fitpath)"
const figpath = "../fitdata/$(fitpath)"

# models
const modeldict = Dict("opt" => OptimalModel,
                       "optlapse" => OptimalLapseModel,
                       "optcollapse" => OptimalColLapseModel,
                       "optcolbiaslapse" => OptimalColBiasLapseModel,
                       "adflapse" => ADFLapseModel,
                       "adfcollapse" => ADFColLapseModel,
                       "adfbiaslapse" => ADFBiasLapseModel,
                       "adfcolbiaslapse" => ADFColBiasLapseModel,
                       "randomcolbiaslapse" => RandomColBiasLapseModel,
                       "randomcollapse" => RandomColLapseModel,
                       "gammalapse" => GammaLapseModel,
                       "gammacollapse" => GammaColLapseModel,
                       "gammabiaslapse" => GammaBiasLapseModel,
                       "gammacolbiaslapse" => GammaColBiasLapseModel,
                       "deltalapse" => DeltaLapseModel,
                       "deltacollapse" => DeltaColLapseModel,
                       "deltabiaslapse" => DeltaBiasLapseModel,
                       "deltacolbiaslapse" => DeltaColBiasLapseModel,
                       "deltaexpcolbiaslapse" => DeltaExpColBiasLapseModel,
                       "lcalapse" => LCALapseModel,
                       "lcacollapse" => LCAColLapseModel)
# separate list to preserve order that keys(modeldict) wouldn't have
# const modelnames = ["opt", "optlapse", "optcollapse",
#                     "adflapse", "adfcollapse", "adfbiaslapse", "adfcolbiaslapse",
#                     "gammalapse", "gammacollapse", "gammabiaslapse", "gammacolbiaslapse",
#                     "deltalapse", "deltacollapse", "deltabiaslapse", "deltacolbiaslapse",
#                     "deltaexpcolbiaslapse"]
const modelnames = ["opt", "optlapse", "optcollapse", "optcolbiaslapse",
                    "adflapse", "adfcollapse", "adfbiaslapse", "adfcolbiaslapse",
                    "randomcolbiaslapse", "randomcollapse",
                    "gammalapse", "gammacollapse", "gammabiaslapse", "gammacolbiaslapse",
                    "deltaexpcolbiaslapse", "lcalapse", "lcacollapse"]
# const modelnames = ["opt", "optlapse", "optcollapse",
#                     "adflapse", "adfcollapse", "adfbiaslapse", "adfcolbiaslapse",
#                     "gammalapse", "gammacollapse", "gammabiaslapse", "gammacolbiaslapse",
#                     "deltaexpcolbiaslapse"]
const basemodel = "adfcolbiaslapse"  # base model for model comparison
@assert issubset(Set(modelnames), Set(keys(modeldict))) # unordered comparison

# objective functions
# second tuple component is trials, 1 = ident, 2 = categ, 3 = itlvd
const objfndict = Dict("psychchron1" => (PsychChron1ObjFn, [1, 2]),
                       "psychchron1ident" => (PsychChron1IdentObjFn, [1]),
                       "psychchron1categ" => (PsychChron1CategObjFn, [2]),
                       "psychchron2" => (PsychChron2ObjFn, [1, 2]),
                       "psychchron3" => (PsychChron3ObjFn, [1, 2]),
                       "psychchronseq1" => (PsychChronSeq1ObjFn, [1, 2]),
                       "psychchronseq2" => (PsychChronSeq2ObjFn, [1, 2]),
                       "psychchronseq3" => (PsychChronSeq3ObjFn, [1, 2]),
                       "psychchron1itlvd" => (PsychChron1ItlvdObjFn, [3]))
#const objfns = ["psychchron1", "psychchron1ident", "psychchron2", "psychchron3",
#                "psychchronseq1", "psychchronseq2", "psychchronseq3"]
const objfns = ["psychchron1", "psychchronseq1"]
@assert issubset(Set(objfns), Set(keys(objfndict)))

# model fit criteria, converted to log-evidence in log(10)
BIC(avgllh, ϕn, trials) = (trials * avgllh - 0.5ϕn * log(trials)) / log(10)
AIC(avgllh, ϕn, trials) = (trials * avgllh - ϕn) / log(10)
AICc(avgllh, ϕn, trials) = AIC(avgllh, ϕn, trials) - (ϕn * (ϕn+1) / (trials-ϕn-1)) / log(10)

# returns number of trials in categ/ident dataset as vector
datatrials(ratid::Int=0) = [size(d, 1) / nrats
    for d in (load_data(:categ, ratid), load_data(:ident, ratid))]
datatrials(objfn::AbstractString, ratid::Int=0) = sum(datatrials(ratid)[objfndict[objfn][2]])

# returns dataframe with best-fitting parameters and avg. log-likelihood
function loadmodelfits(modelname, objfnname, ratid::Int=0)

    fitfilename = ratid == 0 ? "$(datapath)/fit_$(modelname)_$(objfnname).csv" :
                               "$(datapath)/fit_$(modelname)_$(objfnname)_rat$(ratid).csv"
    if !isfile(fitfilename)
        error("Could not find fit file $fitfilename")
    end
    d = CSV.read(fitfilename)
    # identify max-llh parameters
    maxllh, maxi = findmax(d[:llh])
    if maxllh != d[end,:llh]
        warning("Last element in $fitfilename is not best-fitting")
    end
    return d[maxi,:]
end

# prints the fit parameters for all objective functions for a given model and given rats
function printfitparams(modelname, ratids::Vector{Int}=[0])
    # use first objective function to get parameter names
    d = loadmodelfits(modelname, objfns[1], ratids[1])
    ϕnames = names(d)
    # print parameters
    objfnlen = maximum([length(s) for s in objfns])
    ratidlen = maximum([floor(Int, ratid > 0 ? log10(ratid) : 0.0)+1
                        for ratid in ratids]) + 3
    ratidstrs = [ratid > 0 ? "rat$(ratid)" : "all" for ratid in ratids]
    println(repeat("-", objfnlen + ratidlen + 1), "  ",
            join([@sprintf("%7s", string(ϕi)) for ϕi in ϕnames], "  "))
    # print fits
    for objfn in objfns
        for ratidx in 1:length(ratids)
            d = loadmodelfits(modelname, objfn, ratids[ratidx])
            if ratidx == 1
                print(repeat(" ", objfnlen - length(objfn)), objfn)
            else
                print(repeat(" ", objfnlen))
            end
            println(repeat(" ", ratidlen - length(ratidstrs[ratidx]) + 1),
                    ratidstrs[ratidx], "  ",
                    join([@sprintf("%7.2f", x) for x in Vector(d)], "  "))
        end
    end
end

# prints the fit quality for all models for the given objective function
function printfitquality(objfn, ratids::Vector{Int}=Int[])
    # remove 0's from ratids's vector, as they are default
    fitmeasures = ["BIC", "AIC", "AICc"]
    ratids = ratids[ratids .!= 0]
    ratnum = length(ratids)
    modelnum = length(modelnames)
    ratidlen = maximum([floor(Int, log10(ratid))+1 for ratid in ratids]) + 3
    ratidstrs = ["rat$(ratid)" for ratid in ratids]
    modellen = max(maximum([length(s) for s in modelnames]), 3) # include "all" here
    println(repeat("-", modellen + ratidlen + 1), "  ",
            join([@sprintf("%10s      ", s) for s in ["avg. llh", "BIC", "AIC", "AICc"]], "  "))
    fitstats = Array{Float64}(undef, ratnum, modelnum, length(fitmeasures)+1)
    fitstatsall = Array{Float64}(undef, modelnum, length(fitmeasures))
    for modelidx = 1:modelnum
        # fits across all rats
        modelname = modelnames[modelidx]
        d = loadmodelfits(modelname, objfn)
        avgllh = d[:llh]
        ϕn = length(d) - 1
        trials = datatrials(objfn)
        fitq = [avgllh; [f(avgllh, ϕn, trials) for f in [BIC, AIC, AICc]]]
        println(repeat(" ", modellen - length(modelname)), modelname, "  ",
                repeat(" ", ratidlen - 3), "all ",
                join([@sprintf("%10.2f      ", x) for x in fitq], "  "))
        fitstatsall[modelidx, :] = fitq[2:end]
        # fits for individual rats
        if length(ratids) == 0
            continue
        end
        for ratidx in 1:ratnum
            ratid = ratids[ratidx]
            d = loadmodelfits(modelname, objfn, ratid)
            avgllh = d[:llh]
            ϕn = length(d) - 1
            trials = datatrials(objfn, ratid)
            fitstats[ratidx,modelidx,:] = [avgllh; [f(avgllh, ϕn, trials) for f in [BIC, AIC, AICc]]]
            println(repeat(" ", modellen + 1), " ",
                    repeat(" ", ratidlen - length(ratidstrs[ratidx])), ratidstrs[ratidx], " ",
                    join([@sprintf("%10.2f      ", x) for x in fitstats[ratidx,modelidx,:]], "  "))
        end
        # averages across rats
        println(repeat(" ", modellen + 1), " ", repeat(" ", ratidlen - 3), "avg ",
                join([@sprintf("%10.2f±%6.2f", mean(fitstats[:,modelidx,i]), √(var(fitstats[:,modelidx,i]) / ratnum))
                      for i in 1:ratnum], " "))
    end
    return fitstats[:,:,2:end], fitstatsall
end

# returns the absolute and relative fit quality plots for all models for the
# given statistics
function genfitqualityplots(fitstats, fitstatsall, ratids::Vector{Int})
    fitmeasurenames = ["BIC", "AIC", "AICc"]
    ratnum = length(ratids)
    modelnum = length(modelnames)
    @assert(ratnum > 0)
    @assert(size(fitstats, 1) == ratnum)
    @assert(size(fitstats, 2) == modelnum)
    @assert(size(fitstats, 3) == length(fitmeasurenames))
    @assert(size(fitstatsall, 1) == modelnum)
    @assert(size(fitstatsall, 2) == length(fitmeasurenames))
    # separate plot per fit measure
    ps = Array{Any}(undef, length(fitmeasurenames))
    for fitmeasureidx = 1:length(fitmeasurenames)
        substats = fitstats[:,:,fitmeasureidx]
        fitmeasurename = fitmeasurenames[fitmeasureidx]
        # average statistics
        fitmean = mean(substats, dims=1)
        fitse = sqrt.(var(substats, dims=1) / ratnum)
        ps[fitmeasureidx] = plot(
            layer(x=1:modelnum, y=fitstatsall[:,fitmeasureidx] / ratnum,
                  Geom.point, Theme(t, default_color=RGB(0.8, 0.0, 0.0))),
            layer(x=1:modelnum, y=fitmean, ymin=fitmean.-fitse, ymax=fitmean.+fitse,
                  Geom.point, Geom.errorbar),
            layer(x=vec(repeat(collect(1:modelnum)', ratnum)), y=vec(substats),
                  Geom.point, Stat.x_jitter(range=0.4), Theme(t, default_color=RGB(0.5, 0.5, 0.5))),
            layer(yintercept=[0.0], Geom.hline(color="black")),
            Guide.xticks(ticks=collect(1:modelnum), orientation=:vertical),
            Scale.x_discrete(labels=(i->modelnames[i])),
            Guide.xlabel(""),
            Guide.ylabel("log model evidence (using $(fitmeasurename); log10)", orientation=:vertical),
            Coord.cartesian(ymax=0.0))
    end
    p_abs = vstack(ps...)
    # plot comparison to base model
    basemodelidx = findall(modelnames .== basemodel)
    fitdiff = fitstats[:, basemodelidx, :] .- fitstats[:, modelnames .!= basemodel, :]
    fitdiffall = fitstatsall[basemodelidx, :] .- fitstatsall[modelnames .!= basemodel, :]
    diffmodels = modelnames[modelnames .!= basemodel]
    ps = Array{Any}(undef, length(fitmeasurenames))
    for fitmeasureidx = 1:length(fitmeasurenames)
        substats = fitdiff[:,:,fitmeasureidx]
        fitmeasurename = fitmeasurenames[fitmeasureidx]
        # average statistics
        fitmean = mean(substats, dims=1)
        fitse = sqrt.(var(substats, dims=1) / ratnum)
        ps[fitmeasureidx] = plot(
            layer(x=1:(modelnum-1), y=fitdiffall[:,fitmeasureidx] / ratnum,
                  Geom.point, Theme(t, default_color=RGB(0.8, 0.0, 0.0))),
            layer(x=1:(modelnum-1), y=fitmean, ymin=fitmean.-fitse, ymax=fitmean.+fitse,
                  Geom.point, Geom.errorbar),
            layer(x=vec(repeat(collect(1:(modelnum-1))', ratnum)), y=vec(substats),
                  Geom.point, Stat.x_jitter(range=0.4), Theme(t, default_color=RGB(0.5, 0.5, 0.5))),
            layer(yintercept=[0.0], Geom.hline(color="black")),
            Guide.xticks(ticks=collect(1:(modelnum-1)), orientation=:vertical),
            Scale.x_discrete(labels=(i->diffmodels[i])),
            Guide.xlabel(""),
            Guide.ylabel("log model evidence (using $(fitmeasurename); log10)\n$(basemodel) vs. others",
                         orientation=:vertical))
    end
    p_rel = vstack(ps...)
    return p_abs, p_rel
end

# plots the fit quality for all models for the given objective function
function plotfitquality(objfn, ratids::Vector{Int})
    fitmeasures = [BIC, AIC, AICc]
    ratnum = length(ratids)
    @assert(ratnum > 0)
    # collect model fit statistics
    modelnum = length(modelnames)
    fitstats = Array{Float64}(undef, ratnum, modelnum, length(fitmeasures))
    fitstatsall = Array{Float64}(undef, modelnum, length(fitmeasures))
    for modelidx = 1:modelnum
        # first global, then for individual rats
        d = loadmodelfits(modelnames[modelidx], objfn)
        fitstatsall[modelidx, :] = [
            fitmeasure(d[:llh], length(d)-1, datatrials(objfn))
            for fitmeasure in fitmeasures]
        for ratidx = 1:ratnum
            d = loadmodelfits(modelnames[modelidx], objfn, ratids[ratidx])
            fitstats[ratidx, modelidx, :] = [
                fitmeasure(d[:llh], length(d)-1, datatrials(objfn, ratids[ratidx]))
                for fitmeasure in fitmeasures]
        end
    end
    # generate and save plots
    p_abs, p_rel = genfitqualityplots(fitstats, fitstatsall, ratids)
    draw(PDF("$figpath/fitqual_$(objfn).pdf", 8.5inch, 11inch), p_abs)
    println("Plot writting to $figpath/fitqual_$(objfn).pdf")
    draw(PDF("$figpath/fitqual_$(objfn)_relativ.pdf", 8.5inch, 11inch), p_rel)
    println("Plot written to $figpath/fitqual_$(objfn)_relativ.pdf")
end

# returns average log-likelihood for one objective function, using another
# objective function for the parameters. It returns this average
# log-likelihood and the number of trials, and the number of parameters
function avgcrossllh(paramobjfn, evalobjfn, modelname, ratid::Int=0)
    # trial sequence to fit data on
    d_ident = load_data(:ident, ratid)
    d_categ = load_data(:categ, ratid)
    # instantiate model and get model parameters
    m = modeldict[modelname]()
    ϕ = Vector(loadmodelfits(modelname, paramobjfn)[ϕnames(m)])
    ϕn = length(ϕ)
    # evaluate objective function
    evalobjfntype = objfndict[evalobjfn][1]
    evalobjfntrials = datatrials(evalobjfn, ratid)
    Random.seed!(0)
    evidence, ids = catevidence(ident_cs, categ_cs, 10, 10^5)
    o = evalobjfntype(m, evidence, ids, burnin, d_ident, d_categ, nrats)
    avgllh = llh(o, ϕ)
    return avgllh, evalobjfntrials, ϕn
end

# prints the cross fit quality, using parameters for the paramobjfn and
# evaluating it on the evalobjfn. It returns an array of fitmeasures per rat
# and model, and one per model for the meta rat (fits across all rats)
function printcrossfitquality(paramobjfn, evalobjfn, ratids::Vector{Int})
    fitmeasures = [(BIC, "BIC"), (AIC, "AIC"), (AICc, "AICc")]
    ratids = ratids[ratids .!= 0]
    ratnum = length(ratids)
    ratidlen = maximum([floor(Int, log10(ratid))+1 for ratid in ratids]) + 3
    ratidstrs = ["rat$(ratid)" for ratid in ratids]
    @assert(ratnum > 0)
    # collect model fit statistics
    modelnum = length(modelnames)
    modellen = max(maximum([length(s) for s in modelnames]), 3) # include "all" here
    fitstats = Array{Float64}(undef, ratnum, modelnum, length(fitmeasures)+1)
    fitstatsall = Array{Float64}(undef, modelnum, length(fitmeasures))
    for modelidx = 1:modelnum
        modelname = modelnames[modelidx]
        # fits across all rats
        avgllh, evalobjfntrials, ϕn = avgcrossllh(paramobjfn, evalobjfn, modelname)
        fitq = [avgllh; [f(avgllh, ϕn, evalobjfntrials) for f in [BIC, AIC, AICc]]]
        println(repeat(" ", modellen - length(modelname)), modelname, "  ",
                repeat(" ", ratidlen - 3), "all ",
                join([@sprintf("%10.2f      ", x) for x in fitq], "  "))
        fitstatsall[modelidx, :] = fitq[2:end]
        # fits for individual rats
        if length(ratids) == 0
            continue
        end
        avgllhs = Array{Float64}(undef, ratnum)
        for ratidx in 1:ratnum
            ratid = ratids[ratidx]
            avgllh, evalobjfntrials, ϕn = avgcrossllh(paramobjfn, evalobjfn, modelname, ratid)
            fitq = [avgllh; [f(avgllh, ϕn, evalobjfntrials) for f in [BIC, AIC, AICc]]]
            println(repeat(" ", modellen + 1), " ",
                    repeat(" ", ratidlen - length(ratidstrs[ratidx])), ratidstrs[ratidx], " ",
                    join([@sprintf("%10.2f      ", x) for x in fitq], "  "))
            fitstats[ratidx,modelidx,:] = fitq
        end
        # averages across rats
        println(repeat(" ", modellen + 1), " ", repeat(" ", ratidlen - 3), "avg ",
                join([@sprintf("%10.2f±%6.2f", mean(fitstats[:,modelidx,i]),
                               √(var(fitstats[:,modelidx,i]) / ratnum))
                      for i in 1:ratnum], " "))
    end
    # don't return the avgllh fit for individual rats
    return fitstats[:,:,2:end], fitstatsall
end

# writes the fit quality measures to a MAT file
function writefitquality(fitstats, fitstatsall, ratids::Vector{Int}, fname)
    fitmeasurenames = ["BIC", "AIC", "AICc"]
    ratnum = length(ratids)
    modelnum = length(modelnames)
    @assert(ratnum > 0)
    @assert(size(fitstats, 1) == ratnum)
    @assert(size(fitstats, 2) == modelnum)
    @assert(size(fitstats, 3) == length(fitmeasurenames))
    @assert(size(fitstatsall, 1) == modelnum)
    @assert(size(fitstatsall, 2) == length(fitmeasurenames))
    # write matrices to MAT file
    d = Dict(
        "fitstats" => fitstats,
        "fitstatsall" => fitstatsall,
        "fitmeasures" => fitmeasurenames,
        "modelnames" => modelnames,
        "ratids" => ratids
        )
    matwrite(fname, d)
end

# print model parameters
# for modelname in modelnames
#     println("Model $modelname")
#     printfitparams(modelname, [0, 1, 2, 3, 4])
#     println()
# end
# println()

# print model fit statistics
for objfn in objfns
    println("Objective function $objfn")
    fitstats, fitstatsall = printfitquality(objfn, [0, 1, 2, 3, 4])
    fname = "$datapath/fitqual_$(objfn).mat"
    writefitquality(fitstats, fitstatsall, [1, 2, 3, 4], fname)
    println("Data written to $fname")
    println()
end

# plot model fit stats
for objfn in objfns
    println("Generating plots for objective function $objfn")
    plotfitquality(objfn, [1, 2, 3, 4])
end

# # print model fit statistics when evaluated with psychchronseq1 objective
println("Fitted on psychchron1, evaluated on psychchronseq1")
fitstats, fitstatsall = printcrossfitquality("psychchron1", "psychchronseq1", [1, 2, 3, 4])
fname = "$datapath/fitqual_psychchron1_to_psychchronseq1.mat"
writefitquality(fitstats, fitstatsall, [1, 2, 3, 4], fname)
println("Data written to $fname")
println()
println("Generating plots")
p_abs, p_rel = genfitqualityplots(fitstats, fitstatsall, [1, 2, 3, 4])
draw(PDF("$figpath/fitqual_psychchron1_to_psychchronseq1.pdf", 8.5inch, 11inch), p_abs)
println("Plot written to $figpath/fitqual_psychchron1_to_psychchronseq1.pdf")
draw(PDF("$figpath/fitqual_psychchron1_to_psychchronseq1_relativ.pdf", 8.5inch, 11inch), p_rel)
println("Plot written to $figpath/fitqual_psychchron1_to_psychchronseq1_relativ.pdf")

ext: jl


text: println("Hiiiiiiii from julia ", VERSION)


text:
using DataFrames
using SparseArrays

export IndexedDF, getData, getCount, removeSamples, getValues, valueMean
export FastIDF

mutable struct IndexedDF
  df::DataFrame
  index::Vector{Vector{Vector{Int64}}}

  function IndexedDF(df::DataFrame, dims::Vector{Int64})
    ## indexing all columns D - 1 columns (integers)
    index = [ [Int64[] for j in 1:i] for i in dims ]
    for i in 1:size(df, 1)
      for mode in 1:length(dims)
        j = df[i, mode]
        push!(index[mode][j], i)
      end
    end
    new(df, index)
  end
end

IndexedDF(df::DataFrame, dims::Tuple) = IndexedDF(df, Int64[i for i in dims])
IndexedDF(df::DataFrame) = IndexedDF(df, Int64[maximum(df[:,i]) for i in 1 : size(df,2)-1])

valueMean(idf::IndexedDF) = mean(idf.df[:,end])
import Base.size
size(idf::IndexedDF) = tuple( [length(i) for i in idf.index]... )
size(idf::IndexedDF, i::Integer) = length(idf.index[i])

nnz(idf::IndexedDF) = size(idf.df, 1)

function removeSamples(idf::IndexedDF, samples)
  df = idf.df[ setdiff(1:size(idf.df, 1), samples), :]
  return IndexedDF(df, size(idf))
end

getValues(idf::IndexedDF) = convert(Array, idf.df[:, end])
getMode(idf::IndexedDF, mode::Integer) = idf.df[:, mode]
getData(idf::IndexedDF, mode::Integer, i::Integer)  = idf.df[ idf.index[mode][i], :]
getCount(idf::IndexedDF, mode::Integer, i::Integer) = length( idf.index[mode][i] )
getI(idf::IndexedDF, mode::Integer, i::Integer)     = idf.index[mode][i]

## FastIDF used in sampling of latent variables
mutable struct FastIDF{Ti,Tv}
  ids::Matrix{Ti}
  values::Vector{Tv}
  index::Vector{Vector{Vector{Int64}}}
end

FastIDF(idf::IndexedDF) = FastIDF(convert(Matrix, idf.df[:,1:end-1]), convert(Array, idf.df[:,end]), idf.index)
function FastIDF(df::DataFrame, dims::Vector{Int64})
  ## indexing all columns D - 1 columns (integers)
  index = [ [Int64[] for j in 1:i] for i in dims ]
  for i in 1:size(df, 1)
    for mode in 1:length(dims)
      j = df[i, mode]
      push!(index[mode][j], i)
    end
  end
  ix = convert(Matrix, df[:, 1:end-1])
  v  = convert(Array, df[:, end])
  return FastIDF(ix, v, index)
end

function getData(f::FastIDF, mode::Integer, i::Integer)
  id = f.index[mode][i]
  return f.ids[id,:], f.values[id]
end

import Base.size
size(f::FastIDF) = tuple( [length(i) for i in f.index]... )
size(f::FastIDF, i::Integer) = length(f.index[i])
nnz(f::FastIDF) = size(f.values, 1)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -47.221875796417265 <= q <= 49.11122836315772)

                     add_NL_constraint(m, :(gelu(0.42790162605340454 + 0.9948146339173967*gelu(0.30429674417798 + 0.01326158862793747*gelu(0.7020479917236249 + 0.5157446936496526*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.8506753528903843*gelu(0.7594455380381384 + 0.9088497207479063*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.7699605184902198*gelu(0.5333438081093944 + -0.9884706408297048*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.9481609482474158*gelu(-0.8688273038537231 + -0.5762764679031651*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + 0.5690258167368487*gelu(0.4003333530904172 + 0.4565361301167439*gelu(0.7020479917236249 + 0.5157446936496526*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.41218949132242644*gelu(0.7594455380381384 + 0.9088497207479063*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.6719961553691296*gelu(0.5333438081093944 + -0.9884706408297048*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.9454422589116072*gelu(-0.8688273038537231 + -0.5762764679031651*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + -0.632632285898699*gelu(-0.4740061987673023 + -0.21503151130953757*gelu(0.7020479917236249 + 0.5157446936496526*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.7641417558802628*gelu(0.7594455380381384 + 0.9088497207479063*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.8796988576372375*gelu(0.5333438081093944 + -0.9884706408297048*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.6466469730144713*gelu(-0.8688273038537231 + -0.5762764679031651*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + -0.0799939753544896*gelu(-0.9533597279323933 + 0.019939841305541606*gelu(0.7020479917236249 + 0.5157446936496526*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.8134733780415981*gelu(0.7594455380381384 + 0.9088497207479063*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.10348468441851777*gelu(0.5333438081093944 + -0.9884706408297048*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.008542712959329624*gelu(-0.8688273038537231 + -0.5762764679031651*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + gelu(-0.31282642357163937 + 0.400776902433146*gelu(0.30429674417798 + 0.01326158862793747*gelu(0.7020479917236249 + 0.5157446936496526*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.8506753528903843*gelu(0.7594455380381384 + 0.9088497207479063*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.7699605184902198*gelu(0.5333438081093944 + -0.9884706408297048*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.9481609482474158*gelu(-0.8688273038537231 + -0.5762764679031651*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + -0.8700113611274283*gelu(0.4003333530904172 + 0.4565361301167439*gelu(0.7020479917236249 + 0.5157446936496526*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.41218949132242644*gelu(0.7594455380381384 + 0.9088497207479063*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.6719961553691296*gelu(0.5333438081093944 + -0.9884706408297048*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.9454422589116072*gelu(-0.8688273038537231 + -0.5762764679031651*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + -0.8173535730312884*gelu(-0.4740061987673023 + -0.21503151130953757*gelu(0.7020479917236249 + 0.5157446936496526*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.7641417558802628*gelu(0.7594455380381384 + 0.9088497207479063*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.8796988576372375*gelu(0.5333438081093944 + -0.9884706408297048*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.6466469730144713*gelu(-0.8688273038537231 + -0.5762764679031651*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + -0.10383649475667989*gelu(-0.9533597279323933 + 0.019939841305541606*gelu(0.7020479917236249 + 0.5157446936496526*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.8134733780415981*gelu(0.7594455380381384 + 0.9088497207479063*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.10348468441851777*gelu(0.5333438081093944 + -0.9884706408297048*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.008542712959329624*gelu(-0.8688273038537231 + -0.5762764679031651*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + gelu(0.07748802701261592 + -0.48287575350847955*gelu(0.30429674417798 + 0.01326158862793747*gelu(0.7020479917236249 + 0.5157446936496526*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.8506753528903843*gelu(0.7594455380381384 + 0.9088497207479063*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.7699605184902198*gelu(0.5333438081093944 + -0.9884706408297048*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.9481609482474158*gelu(-0.8688273038537231 + -0.5762764679031651*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + 0.9032398979921337*gelu(0.4003333530904172 + 0.4565361301167439*gelu(0.7020479917236249 + 0.5157446936496526*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.41218949132242644*gelu(0.7594455380381384 + 0.9088497207479063*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.6719961553691296*gelu(0.5333438081093944 + -0.9884706408297048*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.9454422589116072*gelu(-0.8688273038537231 + -0.5762764679031651*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + 0.8154202752057489*gelu(-0.4740061987673023 + -0.21503151130953757*gelu(0.7020479917236249 + 0.5157446936496526*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.7641417558802628*gelu(0.7594455380381384 + 0.9088497207479063*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.8796988576372375*gelu(0.5333438081093944 + -0.9884706408297048*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.6466469730144713*gelu(-0.8688273038537231 + -0.5762764679031651*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + 0.8885160528665201*gelu(-0.9533597279323933 + 0.019939841305541606*gelu(0.7020479917236249 + 0.5157446936496526*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.8134733780415981*gelu(0.7594455380381384 + 0.9088497207479063*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.10348468441851777*gelu(0.5333438081093944 + -0.9884706408297048*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.008542712959329624*gelu(-0.8688273038537231 + -0.5762764679031651*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + gelu(0.6389491685818425 + -0.9024035564379238*gelu(0.30429674417798 + 0.01326158862793747*gelu(0.7020479917236249 + 0.5157446936496526*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.8506753528903843*gelu(0.7594455380381384 + 0.9088497207479063*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.7699605184902198*gelu(0.5333438081093944 + -0.9884706408297048*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.9481609482474158*gelu(-0.8688273038537231 + -0.5762764679031651*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + -0.8688708754259227*gelu(0.4003333530904172 + 0.4565361301167439*gelu(0.7020479917236249 + 0.5157446936496526*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.41218949132242644*gelu(0.7594455380381384 + 0.9088497207479063*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.6719961553691296*gelu(0.5333438081093944 + -0.9884706408297048*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.9454422589116072*gelu(-0.8688273038537231 + -0.5762764679031651*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + -0.9399903753655909*gelu(-0.4740061987673023 + -0.21503151130953757*gelu(0.7020479917236249 + 0.5157446936496526*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.7641417558802628*gelu(0.7594455380381384 + 0.9088497207479063*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.8796988576372375*gelu(0.5333438081093944 + -0.9884706408297048*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.6466469730144713*gelu(-0.8688273038537231 + -0.5762764679031651*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2])))) + 0.5656846187923494*gelu(-0.9533597279323933 + 0.019939841305541606*gelu(0.7020479917236249 + 0.5157446936496526*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.8162922098274223*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6507571114046486*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + 0.3903970834590602*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.8134733780415981*gelu(0.7594455380381384 + 0.9088497207479063*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + 0.34867162336420376*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + -0.586942758624089*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.05285512797069236*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + -0.10348468441851777*gelu(0.5333438081093944 + -0.9884706408297048*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.9213430875393112*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.6584439947423988*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.39397409418232554*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))) + 0.008542712959329624*gelu(-0.8688273038537231 + -0.5762764679031651*gelu(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2])) + -0.4115583086286181*gelu(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2])) + 0.07598378122288318*gelu(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2])) + -0.5005567670062852*gelu(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(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:
using DSPopt

include("farmer_models.jl")

m = farmer_stochastic()
DSPopt.classifyConstrs(m)

@testset "Parent model" begin
    @test length(m.variables) == 3
    @test length(m.constraints) == 1
    start, index, value, rlbd, rubd, obj, clbd, cubd, ctype, cname, nqrows, linnzcnt, quadnzcnt, rhs, sense, linstart, linind, linval, quadstart, quadrow, quadcol, quadval = DSPopt.get_model_data(m)
    @test length(start) == length(m.constraints) + 1
    @test start[end] == length(index)
    @test start == [0, 3]
    @test index == [0, 1, 2]
    @test value == [1., 1., 1.]
    @test rlbd == [-Inf]
    @test rubd == [Budget]
    @test obj == Cost
    @test clbd == zeros(3)
    @test cubd == zeros(3) .+ Inf
    @test ctype == "III"
end

@testset "Child model $i" for (i,subm) in m.children
    @test i > 0
    @test m.probability[i] == probability[i]
    @test length(subm.variables) == 6
    @test length(subm.constraints) == 4
    start, index, value, rlbd, rubd, obj, clbd, cubd, ctype, cname, nqrows, linnzcnt, quadnzcnt, rhs, sense, linstart, linind, linval, quadstart, quadrow, quadcol, quadval = DSPopt.get_model_data(subm, i)
    @test length(start) == length(subm.constraints) + 1
    @test start[end] == length(index)
    @test start == [0, 3, 6, 9, 10]
    @test index == [0, 3, 5, 1, 4, 6, 2, 7, 8, 7]
    @test value == [Yield[i,1], 1., -1., Yield[i,2], 1., -1., Yield[i,3], -1., -1., 1.]
    @test rlbd == [Minreq; -Inf]
    @test rubd == [Inf, Inf, Inf, 6000]
    @test obj == [Purchase; -Sell]
    @test clbd == zeros(6)
    @test cubd == zeros(6) .+ Inf
    @test ctype == "CCCCCC"
end

@testset "optimize!: $j" for j in instances(DSPopt.Methods)
    dsp.is_stochastic = true
    status = DSPopt.optimize!(m, is_stochastic = true, solve_type = j)
    @test DSPopt.termination_status(m) == MOI.OPTIMAL
    if dsp.solve_type in [DSPopt.DW, DSPopt.ExtensiveForm]
        @test isapprox(objective_value(m), -108390.)
    else
        @test objective_value(m) >= -108390.
    end
    @test isapprox(dual_objective_value(m), -108390., atol=0.1)

    x = m[:x]
    primsol = value()
    dualsol = dual()
    if dsp.solve_type == DSPopt.Dual
        for k = 0:3
            @test primsol[k] != []
        end
        @test dualsol != []
    else
        @test isapprox(primsol[0], [170.0, 80.0, 250.0])
        if dsp.solve_type != DSPopt.Benders
            @test isapprox(primsol[1], [0.0, 0.0, 310.0, 48.0, 6000.0, 0.0])
            @test isapprox(primsol[2], [0.0, 0.0, 225.0, 0.0, 5000.0, 0.0])
            @test isapprox(primsol[3], [0.0, 48.0, 140.0, 0.0, 4000.0, 0.0])
        end
        @test dualsol == []
        @test isapprox(value(x[1]), 170.0)
        @test isapprox(value(x[2]), 80.0)
        @test isapprox(value(x[3]), 250.0)
    end
    DSPopt.freeSolver(dsp)
    @testset "freeModel" begin
        DSPopt.freeModel(dsp)
        @test dsp.p != C_NULL
        @test length(dsp.numRows) == 0
        @test length(dsp.numCols) == 0
        @test isnan(dsp.primVal)
        @test isnan(dsp.dualVal)
        @test length(dsp.colVal) == 0
        @test length(dsp.rowVal) == 0
        @test dsp.nblocks == -1
        @test dsp.block_ids == []
        @test dsp.is_stochastic == false
        @test dsp.solve_type == DSPopt.Dual
        @test length(dsp.quadConstrs) == 0
        @test length(dsp.linConstrs) == 0
    end
end


text:
function EcksteinHechlerPropagator(arg0::Orbit, arg1::AttitudeProvider, arg10::PropagationType, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return EcksteinHechlerPropagator((Orbit, AttitudeProvider, PropagationType, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function EcksteinHechlerPropagator(arg0::Orbit, arg1::AttitudeProvider, arg2::UnnormalizedSphericalHarmonicsProvider)
    return EcksteinHechlerPropagator((Orbit, AttitudeProvider, UnnormalizedSphericalHarmonicsProvider), arg0, arg1, arg2)
end

function EcksteinHechlerPropagator(arg0::Orbit, arg1::AttitudeProvider, arg2::jdouble, arg3::UnnormalizedSphericalHarmonicsProvider)
    return EcksteinHechlerPropagator((Orbit, AttitudeProvider, jdouble, UnnormalizedSphericalHarmonicsProvider), arg0, arg1, arg2, arg3)
end

function EcksteinHechlerPropagator(arg0::Orbit, arg1::AttitudeProvider, arg2::jdouble, arg3::UnnormalizedSphericalHarmonicsProvider, arg4::PropagationType)
    return EcksteinHechlerPropagator((Orbit, AttitudeProvider, jdouble, UnnormalizedSphericalHarmonicsProvider, PropagationType), arg0, arg1, arg2, arg3, arg4)
end

function EcksteinHechlerPropagator(arg0::Orbit, arg1::AttitudeProvider, arg2::jdouble, arg3::UnnormalizedSphericalHarmonicsProvider, arg4::UnnormalizedSphericalHarmonicsProvider_UnnormalizedSphericalHarmonics)
    return EcksteinHechlerPropagator((Orbit, AttitudeProvider, jdouble, UnnormalizedSphericalHarmonicsProvider, UnnormalizedSphericalHarmonicsProvider_UnnormalizedSphericalHarmonics), arg0, arg1, arg2, arg3, arg4)
end

function EcksteinHechlerPropagator(arg0::Orbit, arg1::AttitudeProvider, arg2::jdouble, arg3::UnnormalizedSphericalHarmonicsProvider, arg4::UnnormalizedSphericalHarmonicsProvider_UnnormalizedSphericalHarmonics, arg5::PropagationType)
    return EcksteinHechlerPropagator((Orbit, AttitudeProvider, jdouble, UnnormalizedSphericalHarmonicsProvider, UnnormalizedSphericalHarmonicsProvider_UnnormalizedSphericalHarmonics, PropagationType), arg0, arg1, arg2, arg3, arg4, arg5)
end

function EcksteinHechlerPropagator(arg0::Orbit, arg1::AttitudeProvider, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble)
    return EcksteinHechlerPropagator((Orbit, AttitudeProvider, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8)
end

function EcksteinHechlerPropagator(arg0::Orbit, arg1::AttitudeProvider, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return EcksteinHechlerPropagator((Orbit, AttitudeProvider, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function EcksteinHechlerPropagator(arg0::Orbit, arg1::UnnormalizedSphericalHarmonicsProvider)
    return EcksteinHechlerPropagator((Orbit, UnnormalizedSphericalHarmonicsProvider), arg0, arg1)
end

function EcksteinHechlerPropagator(arg0::Orbit, arg1::UnnormalizedSphericalHarmonicsProvider, arg2::PropagationType)
    return EcksteinHechlerPropagator((Orbit, UnnormalizedSphericalHarmonicsProvider, PropagationType), arg0, arg1, arg2)
end

function EcksteinHechlerPropagator(arg0::Orbit, arg1::jdouble, arg2::UnnormalizedSphericalHarmonicsProvider)
    return EcksteinHechlerPropagator((Orbit, jdouble, UnnormalizedSphericalHarmonicsProvider), arg0, arg1, arg2)
end

function EcksteinHechlerPropagator(arg0::Orbit, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble)
    return EcksteinHechlerPropagator((Orbit, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function EcksteinHechlerPropagator(arg0::Orbit, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble)
    return EcksteinHechlerPropagator((Orbit, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8)
end

function propagate_orbit(obj::EcksteinHechlerPropagator, arg0::AbsoluteDate)
    return jcall(obj, "propagateOrbit", CartesianOrbit, (AbsoluteDate,), arg0)
end

function reset_initial_state(obj::EcksteinHechlerPropagator, arg0::SpacecraftState)
    return jcall(obj, "resetInitialState", void, (SpacecraftState,), arg0)
end

function reset_initial_state(obj::EcksteinHechlerPropagator, arg0::SpacecraftState, arg1::PropagationType)
    return jcall(obj, "resetInitialState", void, (SpacecraftState, PropagationType), arg0, arg1)
end


text:
export iauBi00
"""
Frame bias components of IAU 2000 precession-nutation models (part
of MHB2000 with additions).

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

Status:  canonical model.

Returned:
    dpsibi,depsbi  double  longitude and obliquity corrections
    dra            double  the ICRS RA of the J2000.0 mean equinox

Notes:

    1. The frame bias corrections in longitude and obliquity (radians)
    are required in order to correct for the offset between the GCRS
    pole and the mean J2000.0 pole.  They define, with respect to the
    GCRS frame, a J2000.0 mean pole that is consistent with the rest
    of the IAU 2000A precession-nutation model.

    2. In addition to the displacement of the pole, the complete
    description of the frame bias requires also an offset in right
    ascension.  This is not part of the IAU 2000A model, and is from
    Chapront et al. (2002).  It is returned in radians.

    3. This is a supplemented implementation of one aspect of the IAU
    2000A nutation model, formally adopted by the IAU General
    Assembly in 2000, namely MHB2000 (Mathews et al. 2002).

References:

    Chapront, J., Chapront-Touze, M. & Francou, G., Astron.
    Astrophys., 387, 700, 2002.

    Mathews, P.M., Herring, T.A., Buffet, B.A., "Modeling of nutation
    and precession   New nutation series for nonrigid Earth and
    insights into the Earth's interior", J.Geophys.Res., 107, B4,
    2002.  The MHB2000 code itself was obtained on 9th September 2002
    from ftp://maia.usno.navy.mil/conv2000/chapter5/IAU2000A.

This revision:  2013 June 18

SOFA release 2018-01-30

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

# void iauBi00(double *dpsibi, double *depsbi, double *dra)

function iauBi00()
    # Allocate return value
    ref_dpsibi = Ref{Float64}(0.0)
    ref_depsbi = Ref{Float64}(0.0)
    ref_dra    = Ref{Float64}(0.0)

    ccall((:iauBi00, libsofa_c), Cvoid,
        (Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
        ref_dpsibi, ref_depsbi, ref_dra)

    return ref_dpsibi[], ref_depsbi[], ref_dra[]
end


text:
module arg_handler

const doc = """This is the Julia vobs program
I can help you with:
Convert VFLD to SQLite (cmd: vobs_to_sqlite)

Usage:
  vobs.jl [vobs_to_sqlite] [options]
  vobs.jl -h | --help

Options:
  -h --help                     Show this screen.
  --starttime=<YYYY-MM-DD-HH>   start time.
  --endtime=<YYYY-MM-DD-HH>     end time.
  --indir=<str>                 Input directory [default: ~/].

vobs_to_sqlite:
  --sqlite-file=<str>           Output to files [default: ~/out.db]

Development shortcuts:
    julia --project=~/git/vobs/ vobs_util.jl vobs_to_sqlite --starttime=2013-10-01-00 --endtime=2013-10-02-00 --sqlite-file=/home/kah/git/vobs/scr/test.db --indir=/home/kah/tmp/vobs

Notes:
    Assumes the naming of vobs files follows a structure like so:
    vobs2013103123 (must end with 10 date characters)

    It is recommended not to cover more than one month with "vobs_to_sqlite" command.
"""

const avail_commands = ["vobs_to_sqlite"]

function main_args(args)
    help_msg(args)
    cmds = make_dict(args)
    cmd_message = command_msg(cmds)
    return cmd_message
end

function help_msg(args)
    """Checks if just the doc string should be printed"""
    if "-h" in args || "--help" in args
        print(doc)
        exit(0)
    end
end

function make_dict(args)
    """Makes a dictionary of input arguments"""
    cmds = Dict()

    for k in args
        if startswith(k,"--")       # Identifier for an option
            opt = split(k,"=")
            cmds[opt[1]] = opt[2]
        elseif startswith(k,"-")    # Identifier for a flag
            cmds[k] = "flag"
        else                        # Identifier for a command
            cmds[k] = "command"
        end
    end
    return cmds
end

function command_msg(cmds)
    """Checks which command we are supposed to run"""

    for cmd in avail_commands

        # command: vobs_to_sqlite
        if cmd == "vobs_to_sqlite" && haskey(cmds,cmd) # Determine which command to run
            cmd_message = check_options(cmds,cmd)
            return cmd_message # We only allow one command at a time for now
        elseif cmd == "date_naming" && haskey(cmds,cmd)
            cmd_message = check_options(cmds,cmd)
            return cmd_message
        end

    end
end

function check_options(cmds, cmd)
     """Further check for necessary options for this specific command"""
    if cmd == "vobs_to_sqlite"

        starttime = key_check("--starttime", cmds)
        endtime = key_check("--endtime", cmds)

        if starttime === missing || endtime === missing
            raise_missing_option("--starttime and --endtime")
        end

        indir = key_check("--indir", cmds)
        indir===missing ? indir = "~/" : nothing

        sqlitefile = key_check("--sqlite-file", cmds)
        sqlitefile===missing ? sqlitefile = "~/out.db" : nothing

        cmd_message = String.((cmd, starttime, endtime, indir, sqlitefile))

        return cmd_message
    else
        print("ERROR. No such option: $cmd")
        exit(0)
    end

end

function key_check(key,dict)
    """Checks if a key exist in a dict and if so, return the value"""
    if haskey(dict,key)
        value = dict[key]
        return value
    else
        return missing
    end
end

function raise_missing_option(option)
    println(option*" necessary for this command")
    exit(0)
end

end


text:
# ---------------------------------------------------------------------------- #
#
#   meshF2D.jl
#
#   Type for 2D frame meshes
#   Inherits from the abstract "meshF" type
#
#   intrico
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    MeshF2D

MeshF2D type:
Type for 2D meshes consisting of triangles. It holds the nodes and connectivity
information.
"""
struct MeshF2D <: MeshF

    dim::Int64              # Dimension of the problem

    n::Int64                # Number of nodes

    p::Matrix{Float64}      # Nodal locations
    t::Matrix{Int64}        # Triangle - node connectivity
    t2f::Matrix{Int64}      # Triangle - face connectivity
    f::Matrix{Int64}        # Face - node/triangle connectivity
    fb::Matrix{Int64}       # Boundary face info
    nodes::Array{Float64,3} # Nodes on which solution is evaluated -- needed to
                          #  compute Jacobian on face
    n2e::Vector{Vector{Int64}} # Node to edge connectivity

    boundtags::Vector{String} # boundary tags

end

"""
    MeshF( mesh::divido.Mesh2D )

Constructor that generates frame from `divido` mesh structure in `mesh`.
"""
function MeshF( mesh::divido.Mesh2D )

    n2e = genNEconnec2D( mesh.f, mesh.n )

    MeshF2D( 2, mesh.n, mesh.p, mesh.t, mesh.t2f, mesh.f, mesh.fb, mesh.nodes, n2e, mesh.boundtags )

end

"""
    MeshF2D( name::String )

Constructor for frame from mesh written in `name`.
"""
function MeshF2D( name::String )

    if name[end-3:end] == ".su2"
        (p_, t_, bel_, tags_) = divido.readSU2_2D( name )
    elseif name[end-3:end] == ".msh" # BAMG
        (p_, t_, bel_, tags_ ) = divido.readBAMG( name )
    elseif name[end-4:end] == ".mesh" # FEFLOA
        (p_, t_, bel_, tags_ ) = divido.readFEFLOA_2D( name )
    else
        error("MeshF2D: Unknown mesh type")
    end

    porder_ = divido.P1()

    (f_, t2f_, nodes_, ploc_, tloc_, fb_) = divido.genmesh( porder_, p_, t_, bel_ )

    n_ = size( p_, 1 )

    n2e_ = genNEconnec2D( f_, n_ )

    MeshF2D( 2, n_, p_, t_, t2f_, f_, fb_, nodes_, n2e_, tags_ )

end

"""
    genNEconnec2D( e::Matrix{Int64}, n::Int64 )

Generates node to edge connectivity given edge connectivity (`e`).
"""
function genNEconnec2D( e::Matrix{Int64}, n::Int64 )

    n2e = fill( fill(0,0), n )

    for jj in 1:size(e,1)

        ii1 = e[jj,1]
        ii2 = e[jj,2]

        if isempty( n2e[ii1] .== jj )
            append!(n2e[ii1], jj)
        end
        if isempty( n2e[ii2] .== jj )
            append!(n2e[ii2], jj)
        end

    end

    for kk in 1:n
        n2e[kk] = sort( n2e[kk] )
    end

    return n2e

end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =========================================================================================
"""
    struct Job{T<:IO}
        event::Condition
        done::Ref{Bool}
        file::T
        buffer::Vector{UInt8}
    end

Encapsulates all the data needed to execute a task that reads the output from Tesseract and
provides it to the Julia application.

__Values:__

| Name   | Description
| :----- | :----------
| event  | The event that gets set when there is more data or the done flag is set.
| done   | The flag that indicates when the tesseract process has finished.
| file   | The file that is being populated by tesseract.
| buffer | The buffer to reuse when reading the data from the file.

__Constructors:__

    Job(event, done, file)

Create a new instance of the object.  The buffer is initialized to the size used by the
reader.
"""
struct Job{T<:IO}
    event::Condition
    done::Ref{Bool}
    file::T
    buffer::Vector{UInt8}
    Job(event::Condition, done::Ref{Bool}, file::T) where T = new{T}(event, done, file, Vector{UInt8}(undef, 4096))
end

# =========================================================================================
"""
    start_job(
        processor::Function,
        event::Condition,
        done::Ref{Bool},
        path::AbstractString,
        filename::AbstractString,
        transcode::Bool = false
    )::Task

Create a job and start a task to execute it.

__Arguments:__

| T | Name      | Default | Description
|:--| :-------- | :------ | :----------
| R | processor |         | The function to call on the task to process the job.
| R | event     |         | The event that will be signalled when there is data.
| R | done      |         | The flag that indicates the job is done.
| R | path      |         | The directory created for the job.
| R | filename  |         | The name of the file Tesseract will be writing to.
| O | transcode | `false` | Should the input stream be converted into UTF8 from Latin-1?
"""
function start_job(
            processor::Function,
            event::Condition,
            done::Ref{Bool},
            path::AbstractString,
            filename::AbstractString,
            transcode::Bool = false
        )::Task
    local file = open(joinpath(path, filename); read = true)
    local job

    if transcode == true
        job = Job(event, done, Latin1(file))
    else
        job = Job(event, done, file)
    end

    return @async begin
        try
            # -----------------------------------------------------------------------------
            # Run the task.
            processor(job)
        finally
            # -----------------------------------------------------------------------------
            # Make sure everything is cleaned up.
            isopen(file) && close(file)
            ispath(path) && rm(path; recursive = true)
        end
    end
end

# =========================================================================================
"""
    job_read(
        j::Job
    )::Bool

Read data from the file, waiting if there is no data to read.  If there is no more data
and the job is done returns `false`.

__Arguments:__

| T | Name | Default    | Description
|:--| :--- | :--------- | :----------
| R | j    |            | The job to read the data for.
"""
function job_read(
            j::Job
        )::Bool
    # -------------------------------------------------------------------------------------
    # Make sure the buffer is correctly sized.
    resize!(j.buffer, 4096)

    # -------------------------------------------------------------------------------------
    # Repeat while we've read nothing.
    local read = 0

    while read == 0

        # ---------------------------------------------------------------------------------
        # Read any available bytes.
        read = readbytes!(j.file, j.buffer)

        # ---------------------------------------------------------------------------------
        # If no bytes where read and the task is done, return false.
        if read == 0 && j.done[] == true
            return false
        end

        # ---------------------------------------------------------------------------------
        # If we're read nothing, wait for the event to trigger before we try again.
        if read == 0
            wait(j.event)
        end

    end

    # ---------------------------------------------------------------------------------
    # Resize the buffer down to the bytes read.
    resize!(j.buffer, read)
    return true
end

# =========================================================================================
"""
    job_create_byte_array(
        j::Job,
        out::TessOutput{Vector{UInt8}}
    )::Nothing

Read the file specified by the job and convert the whole file into a byte array.

__Arguments:__

| T | Name | Default    | Description
|:--| :--- | :--------- | :----------
| R | j    |            | The job this function is being called for.
| R | out  |            | The object to pass the byte array to.
"""
function job_create_byte_array(
            j::Job,
            out::TessOutput{Vector{UInt8}}
        )::Nothing
    local data = Vector{UInt8}()

    while job_read(j) == true
        append!(data, j.buffer)
    end

    out.result = data

    nothing
end

# =========================================================================================
"""
    job_create_string(
        j::Job,
        out::TessOutput{String}
    )::Nothing

Read the file specified by the job and convert the whole file into a string.

__Arguments:__

| T | Name | Default    | Description
|:--| :--- | :--------- | :----------
| R | j    |            | The job this function is being called for.
| R | out  |            | The object to pass the string to.
"""
function job_create_string(
            j::Job,
            out::TessOutput{String}
        )::Nothing
    local data = Vector{UInt8}()

    while job_read(j) == true
        append!(data, j.buffer)
    end

    out.result = String(data)

    nothing
end

# =========================================================================================
"""
    job_dispatch(
        j::Job,
        dispatch::Function,
        separator::String
    )::Nothing

Read the file specified by the job, break the file into lines and call the dispatch
function so they can be processed.

__Arguments:__

| T | Name      | Default | Description
|:--| :-------- | :------ | :----------
| R | j         |         | The job this function is being called for.
| R | dispatch  |         | The function to pass the lines to.
| R | separator |         | The string that Tesseract will be using to separate pages.
"""
function job_dispatch(
            j::Job,
            dispatch::Function,
            separator::String
        )::Nothing
    local sepbytes = Vector{UInt8}(separator)
    local seplen   = length(separator) + 1
    local line = Vector{UInt8}()

    function page_separator(arr)::Bool

        if length(sepbytes) == 0 || length(arr) < length(sepbytes)
            return false
        end
        for i in 1:length(sepbytes)
            @inbounds if sepbytes[i] != arr[i]
                return false
            end
        end
        return true
    end

    function notify(arr)
        if page_separator(arr)
            dispatch(separator)
            if length(arr) >= seplen
                dispatch(String(arr[seplen:length(arr)]))
            end
        else
            dispatch(String(arr))
        end
    end

    # -------------------------------------------------------------------------------------
    # Read the text and process it.
    while job_read(j) == true
        local s = 1

        # ---------------------------------------------------------------------------------
        # Walk the buffer and extract all the lines to pass to the notify function.
        for (i, v) in enumerate(j.buffer)
            if v == 0x0a && length(line) > 0
                append!(line, view(j.buffer, s:i))
                notify(line)
                s = i + 1
            elseif v == 0x0a
                notify(j.buffer[s:i])
                s = i + 1
            end
        end

        # ---------------------------------------------------------------------------------
        # Copy the remainder of the buffer to the line.
        if s <= length(j.buffer)
            append!(line, view(j.buffer, s:length(j.buffer)))
        end
    end

    # -------------------------------------------------------------------------------------
    # Tell the client about the last line.
    if length(line) > 0
        notify(line)
    end

    nothing
end

# =========================================================================================
"""
    job_copy(
        j::Job,
        output::IO
    )::Nothing

Read the contents of the file specified by the job and write it the output file.

__Arguments:__

| T | Name   | Default | Description
|:--| :----- | :------ | :----------
| R | j      |         | The job this function is being called for.
| R | output |         | The function to pass the lines to.

__Details:__

This method is currently used for transcoding the output from from one encoding to another.
"""
function job_copy(
            j::Job,
            output::IO
        )::Nothing
    # -------------------------------------------------------------------------------------
    # Read the text and process it.
    while job_read(j) == true
        write(output, j.buffer)
    end

    nothing
end


text:
#!/usr/bin/env julia

# Draw some unicodes to an image.
# https://github.com/GiovineItalia/Compose.jl/pull/360#issuecomment-539283765

using Compose
import Cairo, Fontconfig

c = compose(context(), text(0.5, 0.5, "𝑎 𝑏 𝑐 𝑑 𝑒 𝑓", hcenter, vcenter), font("Segoe UI"), fontsize(20pt) ) # Note: the string is composed with \ita \itb etc.

imgs = [
    PDF("unicode.pdf")
    SVG("unicode.svg")
    PNG("unicode.png")
]

draw.(imgs, [c])


text:
using CheckedArithmeticCore
using Test

struct MyType end
struct MySafeType end
Base.convert(::Type{MySafeType}, ::MyType) = MySafeType()
CheckedArithmeticCore.safearg_type(::Type{MyType}) = MySafeType
CheckedArithmeticCore.accumulatortype(::typeof(+), ::Type{MyType}) = MyType
CheckedArithmeticCore.accumulatortype(::typeof(+), ::Type{MySafeType}) = MySafeType
CheckedArithmeticCore.accumulatortype(::typeof(*), ::Type{MyType}) = MySafeType
CheckedArithmeticCore.accumulatortype(::typeof(*), ::Type{MySafeType}) = MySafeType
Base.promote_rule(::Type{MyType}, ::Type{MySafeType}) = MySafeType

# fallback
@test safearg(MyType()) === MySafeType()

# fallback
@test safeconvert(UInt16, 0x12) === 0x0012

# fallback
@test accumulatortype(MyType) === MySafeType
@test accumulatortype(MyType, MyType) === MySafeType
@test accumulatortype(+, MyType) === MyType
@test accumulatortype(*, MyType) === MySafeType
@test accumulatortype(+, MyType, MySafeType) === MySafeType
@test accumulatortype(*, MySafeType, MyType) === MySafeType

# acc
@test acc(MyType()) === MySafeType()
@test acc(+, MyType()) === MyType()


text:
using Documenter
using Literate
using Plots
ENV["GKSwstype"] = "nul"  # to avoid GKS error on docs build ("nul" for headless systems)

using LongwaveModePropagator
using LongwaveModePropagator.Fields

#==
Generate examples
==#

const ROOT_DIR = abspath(@__DIR__, "..")
const EXAMPLES_DIR = joinpath(ROOT_DIR, "examples")
const BUILD_DIR = joinpath(ROOT_DIR, "docs", "build")
const OUTPUT_DIR = joinpath(ROOT_DIR, "docs", "src", "generated")

# repo_root_url is only assigned with CI, so we have to specify it for local builds
# can't use abspath on Windows (no colons) so we calculate relpath
const REPO_ROOT_URL = relpath(ROOT_DIR, joinpath(BUILD_DIR, "generated"))

examples = [
    "basic.jl",
    "io.jl",
    "meshgrid.jl",
    "meshgrid2.jl",
    "integratedreflection.jl",
    "wavefieldintegration.jl",
    "magneticfield.jl",
    "interpretinghpbeta.jl",
    "multiplespecies.jl",
    "interpolatingfunctions.jl",
    "ground.jl"
]

for example in examples
    example_filepath = joinpath(EXAMPLES_DIR, example)
    if get(ENV, "CI", nothing) == "true"
        Literate.markdown(example_filepath, OUTPUT_DIR, documenter=true)
    else
        # local
        Literate.markdown(example_filepath, OUTPUT_DIR, documenter=true,
                          repo_root_url=REPO_ROOT_URL)
    end
end

#==
Organize page hierarchies
==#

example_pages = Any[
    "generated/basic.md",
    "generated/io.md",
    "generated/meshgrid.md",
    "generated/meshgrid2.md",
    "generated/integratedreflection.md",
    "generated/wavefieldintegration.md",
    "generated/magneticfield.md",
    "generated/interpretinghpbeta.md",
    "generated/multiplespecies.md",
    "generated/interpolatingfunctions.md",
    "generated/ground.md"
]

library_pages = Any[
    "lib/public.md",
    "lib/internals.md"
]

pages = Any[
    "Home" => "index.md",
    "Examples" => example_pages,
    "Library" => library_pages,
    "References" => "references.md"
]

#==
Build and deploy docs
==#

format = Documenter.HTML(
    collapselevel = 1,
    prettyurls=get(ENV, "CI", nothing)=="true"
)

makedocs(
    sitename = "LongwaveModePropagator",
    authors = "Forrest Gasdia",
    format = format,
    pages = pages,
    modules = [LongwaveModePropagator]
)

deploydocs(
    repo = "github.com/fgasdia/LongwaveModePropagator.jl.git",
    devbranch = "main",
)


text:
model, obj_uncon = TrajectoryOptimization.Dynamics.pendulum

# -Constrained objective
obj = ConstrainedObjective(obj_uncon, u_min=-2, u_max=2, x_min=[-20;-20], x_max=[20;20])

# Solver
solver = Solver(model,obj,integration=:rk4,N=51)
U0 = ones(solver.model.m,solver1.N)
X0 = line_trajectory(solver)

results = init_results(solver,X0,U0)
rollout!(results,solver)
update_jacobians!(results, solver)
update_constraints!(results,solver)
_backwardpass!(results,solver)

results.bp

bp = results.bp
n,m,N = get_sizes(solver)
n̄,nn = get_num_states(solver)
m̄,mm = get_num_controls(solver)
p,pI,pE = get_num_constraints(solver)
p_N,pI_N,pE_N = get_num_terminal_constraints(solver)

nm = nn + mm
Nz = nn*N + mm*(N-1)
Np = p*(N-1) + p_N
∇²L = zeros(Nz,Nz)
∇c = zeros(Np,Nz)

for k = 1:N
    if k < N
        idx = ((k-1)*nm + 1):k*nm
        ∇²L[idx,idx] = [bp.Qxx[k] bp.Qux[k]'; bp.Qux[k] bp.Quu[k]]

        idx2 = ((k-1)*p + 1):k*p
        ∇c[idx2,idx] = [results.Cx[k] results.Cu[k]]
    else
        idx = ((k-1)*nm + 1):Nz
        ∇²L[idx,idx] = results.S[N]

        idx2 = ((k-1)*p + 1):Np
        ∇c[idx2,idx] = results.Cx[N]
    end
end

∇²L
∇c

C = vcat(results.C...)
λ = vcat(results.λ...)

results.λ
λ[1] = 10.
active_set = vcat(results.active_set...)

∇c̄ = ∇c[active_set,:]

λ[active_set] += (∇c̄*(∇²L\∇c̄'))\C[active_set]

# update the results
for k = 1:N
    if k != N
        idx_pI = pI
        idx = (k-1)*p+1:k*p
    else
        idx_pI = pI_N
        idx = (k-1)*p+1:Np
    end
    results.λ[k] = max.(solver.opts.dual_min, min.(solver.opts.dual_max, λ[idx]))
    results.λ[k][1:idx_pI] = max.(0.0,results.λ[k][1:idx_pI])
end


text:
"""
    AbstractInitializer

The abstract base class for all initializers.

To define a new initializer, it is
enough to derive a new type, and implement one or more of the following methods:

    _init_weight(self :: AbstractInitializer, name :: Base.Symbol, array :: NDArray)
    _init_bias(self :: AbstractInitializer, name :: Base.Symbol, array :: NDArray)
    _init_gamma(self :: AbstractInitializer, name :: Base.Symbol, array :: NDArray)
    _init_beta(self :: AbstractInitializer, name :: Base.Symbol, array :: NDArray)

Or, if full behavior customization is needed, override the following function

    init(self :: AbstractInitializer, name :: Base.Symbol, array :: NDArray)
"""
abstract AbstractInitializer

function init{T<:AbstractInitializer}(self :: T, name :: Base.Symbol, array :: NDArray)
  strname = string(name)
  if endswith(strname, "bias")
    _init_bias(self, name, array)
  elseif endswith(strname, "gamma")
    _init_gamma(self, name, array)
  elseif endswith(strname, "beta")
    _init_beta(self, name, array)
  elseif endswith(strname, "weight")
    _init_weight(self, name, array)
  elseif endswith(strname, "moving_mean")
    _init_zero(self, name, array)
  elseif endswith(strname, "moving_var")
    _init_zero(self, name, array)
  else
    _init_default(self, name, array)
  end
end

function _init_bias(self :: AbstractInitializer, name :: Base.Symbol, array :: NDArray)
  array[:] = 0
end
function _init_gamma(self :: AbstractInitializer, name :: Base.Symbol, array :: NDArray)
  array[:] = 1
end
function _init_beta(self :: AbstractInitializer, name :: Base.Symbol, array :: NDArray)
  array[:] = 0
end
function _init_zero(self :: AbstractInitializer, name :: Base.Symbol, array :: NDArray)
  array[:] = 0
end

function _init_default(self :: AbstractInitializer, name :: Base.Symbol, array :: NDArray)
  error("Do not know how to init $name")
end

"""
    UniformInitializer

Initialize weights according to a uniform distribution within the provided scale.
"""
immutable UniformInitializer <: AbstractInitializer
  scale :: AbstractFloat
end
"""
    UniformInitializer(scale=0.07)

Construct a `UniformInitializer` with the specified scale.
"""
UniformInitializer() = UniformInitializer(0.07)

function _init_weight(self :: UniformInitializer, name :: Base.Symbol, array :: NDArray)
  rand!(-self.scale, self.scale, array)
end

"""
    NormalInitializer

Initialize weights according to a univariate Gaussian distribution.
"""
immutable NormalInitializer <: AbstractInitializer
  μ :: AbstractFloat
  σ :: AbstractFloat
end
"""
    NormalIninitializer(; mu=0, sigma=0.01)

Construct a `NormalInitializer` with mean `mu` and variance `sigma`.
"""
NormalInitializer(; mu=0, sigma=0.01) = NormalInitializer(mu, sigma)

function _init_weight(self :: NormalInitializer, name :: Base.Symbol, array :: NDArray)
  randn!(self.μ, self.σ, array)
end

"""
    XavierInitializer

The initializer documented in the paper [Bengio and Glorot 2010]: *Understanding
the difficulty of training deep feedforward neuralnetworks*.

There are several different version of the XavierInitializer used in the wild.
The general idea is that the variance of the initialization distribution is controlled
by the dimensionality of the input and output. As a distribution one can either choose
a normal distribution with μ = 0 and σ² or a uniform distribution from -σ to σ.

Several different ways of calculating the variance are given in the literature or are
used by various libraries.

* [Bengio and Glorot 2010]: `mx.XavierInitializer(distribution = mx.xv_uniform, regularization = mx.xv_avg, magnitude = 1)`
* [K. He, X. Zhang, S. Ren, and J. Sun 2015]: `mx.XavierInitializer(distribution = mx.xv_gaussian, regularization = mx.xv_in, magnitude = 2)`
* caffe_avg: `mx.XavierInitializer(distribution = mx.xv_uniform, regularization = mx.xv_avg, magnitude = 3)`
"""

@enum XavierDistribution xv_uniform xv_normal
@enum XavierRegularization xv_avg xv_in xv_out

immutable XavierInitializer <: AbstractInitializer
  distribution :: XavierDistribution
  regularization :: XavierRegularization
  magnitude :: Float64
end

XavierInitializer(; distribution = xv_uniform, regularization = xv_avg, magnitude = 3.0) = XavierInitializer(distribution, regularization, magnitude)

function _init_weight(self :: XavierInitializer, name :: Base.Symbol, array :: NDArray)
  dims    = size(array)
  fan_in  = prod(dims[2:end])
  fan_out = dims[1]

  if self.regularization == xv_avg
    factor = (fan_in + fan_out) / 2
  elseif self.regularization == xv_in
    factor = fan_in
  elseif self.regularization == xv_out
    factor = fan_out
  end

  σ = √(self.magnitude / factor)

  if self.distribution == xv_uniform
    rand!(-σ, σ, array)
  elseif self.distribution == xv_normal
    randn!(0.0, σ, array)
  end
end


text:
using Gtk.ShortNames

function keypresswindow()

    # This code creates  the Gtk widgets on the screen.
    txt = "Type Y or N"
    win = Window("Keypress Test", 250, 30) |> (Frame() |> ((vbox = Box(:v)) |> (lab = Label(txt))))

    # this is the keystroke processing code, a function and a callback for the function.
    function keycall(w, event)
        ch = Char(event.keyval)
        set_gtk_property!(lab,:label, ch in('n','N','y','Y') ? "You hit the $ch key." : txt)
    end
    Gtk.signal_connect(keycall, win, "key-press-event")

    # this code sets up a proper exit when the widow is closed.
    c = Condition()
    endit(w) = notify(c)
    Gtk.signal_connect(endit, win, :destroy)
    Gtk.showall(win)
    wait(c)
end

keypresswindow()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
name = "Julia"
age = 9