text:
@testset "pow problems" begin
    @testset "pow1_test" begin
    #=
    primal
        max z
    s.t.
        x^0.9 * y^(0.1) >= |z| (i.e (x, y, z) are in the 3d power cone with a=0.9) :w_3, w_4, w_5
        x == 2 :w_2
        y == 1 :w_2

   dual
        min -w_2 - 2w_1
    s.t.
        w_1 + w_3 == 0
        w_2 + w_4 == 0
        w_5 == -1
        (w_3, w_4, w_5) ∈ DualPowerCone
   =#
       primal_model = pow1_test()
       dual_model, primal_dual_map = dual_model_and_map(primal_model)

       @test MOI.get(dual_model, MOI.NumberOfVariables()) == 5
       list_of_cons =  MOI.get(dual_model, MOI.ListOfConstraints())
       @test Set(list_of_cons) == Set([
           (SAF{Float64}, MOI.EqualTo{Float64})
           (VVF, MOI.DualPowerCone{Float64})
       ])
       @test MOI.get(dual_model, MOI.NumberOfConstraints{VVF, MOI.DualPowerCone{Float64}}()) == 1
       @test MOI.get(dual_model, MOI.NumberOfConstraints{SAF{Float64}, MOI.EqualTo{Float64}}()) == 3
       obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
       @test obj_type == SAF{Float64}
       obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
       @test MOI.get(dual_model, MOI.ObjectiveSense()) == MOI.MIN_SENSE
       @test MOI.constant(obj) == 0.0
       @test MOI.coefficient.(obj.terms) == [-1.0; -2.0]

       eq_con1_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
       eq_con1_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
       @test MOI.coefficient.(eq_con1_fun.terms) == [1.0; 1.0]
       @test MOI.constant.(eq_con1_fun) == 0.0
       @test MOI.constant(eq_con1_set) == 0.0
       eq_con2_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
       eq_con2_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
       @test MOI.coefficient.(eq_con2_fun.terms) == [1.0; 1.0]
       @test MOI.constant.(eq_con2_fun) == 0.0
       @test MOI.constant(eq_con2_set) == 0.0
       eq_con3_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(4))
       eq_con3_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(4))
       @test MOI.coefficient.(eq_con3_fun.terms) == [1.0]
       @test MOI.constant.(eq_con3_fun) == 0.0
       @test MOI.constant(eq_con3_set) == -1.0

       dual_exp_con = MOI.get(dual_model, MOI.ConstraintFunction(), CI{VVF, MOI.DualPowerCone}(1))
       @test dual_exp_con.variables == VI.(3:5)

       primal_con_dual_var = primal_dual_map.primal_con_dual_var
       @test primal_con_dual_var[CI{SAF{Float64}, MOI.EqualTo{Float64}}(2)] == [VI(1)]
       @test primal_con_dual_var[CI{SAF{Float64}, MOI.EqualTo{Float64}}(3)] == [VI(2)]
       @test primal_con_dual_var[CI{VVF, MOI.PowerCone{Float64}}(1)] == VI.(3:5)

       primal_var_dual_con = primal_dual_map.primal_var_dual_con
       @test primal_var_dual_con[VI(1)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(2)
       @test primal_var_dual_con[VI(2)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(3)
       @test primal_var_dual_con[VI(3)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(4)
   end

   @testset "pow2_test" begin
    #=
    primal
        max z
    s.t.
        x^0.9 * y^(0.1) >= |z| (i.e (x, y, z) are in the 3d power cone with a=0.9) :w_3, w_4, w_5
        x == 2 :w_2
        y == 1 :w_2

   dual
        min -w_2 - 2w_1
    s.t.
        w_1 + w_3 == 0
        w_2 + w_4 == 0
        w_5 == -1
        (w_3, w_4, w_5) ∈ DualPowerCone
   =#
       primal_model = pow2_test()
       dual_model, primal_dual_map = dual_model_and_map(primal_model)

       @test MOI.get(dual_model, MOI.NumberOfVariables()) == 5
       list_of_cons =  MOI.get(dual_model, MOI.ListOfConstraints())
       @test Set(list_of_cons) == Set([
           (SAF{Float64}, MOI.EqualTo{Float64})
           (VVF, MOI.DualPowerCone{Float64})
       ])
       @test MOI.get(dual_model, MOI.NumberOfConstraints{VVF, MOI.DualPowerCone{Float64}}()) == 1
       @test MOI.get(dual_model, MOI.NumberOfConstraints{SAF{Float64}, MOI.EqualTo{Float64}}()) == 3
       obj_type = MOI.get(dual_model, MOI.ObjectiveFunctionType())
       @test obj_type == SAF{Float64}
       obj = MOI.get(dual_model, MOI.ObjectiveFunction{obj_type}())
       @test MOI.get(dual_model, MOI.ObjectiveSense()) == MOI.MIN_SENSE
       @test MOI.constant(obj) == 0.0
       @test MOI.coefficient.(obj.terms) == [-1.0; -2.0]

       eq_con1_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
       eq_con1_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(2))
       @test MOI.coefficient.(eq_con1_fun.terms) == [1.0; 1.0]
       @test MOI.constant.(eq_con1_fun) == 0.0
       @test MOI.constant(eq_con1_set) == 0.0
       eq_con2_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
       eq_con2_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(3))
       @test MOI.coefficient.(eq_con2_fun.terms) == [1.0; 1.0]
       @test MOI.constant.(eq_con2_fun) == 0.0
       @test MOI.constant(eq_con2_set) == 0.0
       eq_con3_fun = MOI.get(dual_model, MOI.ConstraintFunction(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(4))
       eq_con3_set = MOI.get(dual_model, MOI.ConstraintSet(), CI{SAF{Float64}, MOI.EqualTo{Float64}}(4))
       @test MOI.coefficient.(eq_con3_fun.terms) == [1.0]
       @test MOI.constant.(eq_con3_fun) == 0.0
       @test MOI.constant(eq_con3_set) == -1.0

       dual_exp_con = MOI.get(dual_model, MOI.ConstraintFunction(), CI{VVF, MOI.DualPowerCone}(1))
       @test dual_exp_con.variables == VI.(3:5)

       primal_con_dual_var = primal_dual_map.primal_con_dual_var
       @test primal_con_dual_var[CI{SAF{Float64}, MOI.EqualTo{Float64}}(2)] == [VI(1)]
       @test primal_con_dual_var[CI{SAF{Float64}, MOI.EqualTo{Float64}}(3)] == [VI(2)]
       @test primal_con_dual_var[CI{VAF{Float64}, MOI.PowerCone{Float64}}(1)] == VI.(3:5)

       primal_var_dual_con = primal_dual_map.primal_var_dual_con
       @test primal_var_dual_con[VI(1)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(2)
       @test primal_var_dual_con[VI(2)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(3)
       @test primal_var_dual_con[VI(3)] == CI{SAF{Float64}, MOI.EqualTo{Float64}}(4)
   end
end

ext: jl


text:
module PlutoStaticHTML

using Base64: base64encode
using Pkg:
    Types.Context,
    Types.UUID,
    Operations
using Pluto:
    Cell,
    CellOutput,
    Notebook,
    PkgCompat.dependencies,
    PlutoRunner,
    ServerSession,
    SessionActions,
    WorkspaceManager,
    generate_html,
    load_notebook_nobackup,
    update_dependency_cache!,
    update_run!,
    update_save_run!

include("module_doc.jl")
include("context.jl")
include("html.jl")
include("build.jl")

export notebook2html, run_notebook!
export parallel_build!

end # module


text:
"""
    index(dgt::AbstractVector{T}; base::Integer=2) where T <: Integer

Convert a digits to the integer index using the relation:

`number = ∑ᵢ bits[i] * base^(L-i) + 1`

The summaton is evaluate using the efficieint polynomial evaluation method.
"""
@inline function index(dgt::AbstractVector{T}; base::Integer=2) where T <: Integer
    N = zero(T)
    for i = 1:length(dgt)
        N *= base
        N += dgt[i]
    end
    N + one(T)
end

"""
    index(dgt::AbstractVector{T}, sites::AbstractVector{<:Integer}; base::Integer=2) where T <: Integer

Convert a sub-digits (subarray of `dgt`) to the integer index.
"""
@inline function index(dgt::AbstractVector{T}, sites::AbstractVector{<:Integer}; base::Integer=2) where T <: Integer
    N = zero(T)
    for i in sites
        N *= base
        N += dgt[i]
    end
    N + one(T)
end

"""
    change!(dgt::AbstractVector{T}, ind::Integer; base::Integer=2) where T <: Integer

Change the digits to that with the target index.
This method is the inverse of `index`.
"""
@inline function change!(dgt::AbstractVector{T}, ind::Integer; base::Integer=2) where T <: Integer
    N = convert(T, ind) - one(T)
    for i = length(dgt):-1:1
        N, dgt[i] = divrem(N, base)
    end
end

"""
    change!(dgt::AbstractVector{T}, sites::AbstractVector{<:Integer}, ind::Integer; base::Integer=2) where T <: Integer

Change the sub-digits to that with the target index.
This method is the inverse of `index`.
"""
@inline function change!(dgt::AbstractVector{T}, sites::AbstractVector{<:Integer}, ind::Integer; base::Integer=2) where T <: Integer
    N = convert(T, ind) - one(T)
    for i = length(sites):-1:1
        N, dgt[sites[i]] = divrem(N, base)
    end
end

"""
    binary_search(list::AbstractVector{<:Integer}, i<:Integer)

Return the position of i in a sorted list using binary search algorithm.
"""
function binary_search(list::AbstractVector{<:Integer}, i::Integer)
    l::Int = 1
    r::Int = length(list)
    c::Int = (l + r) ÷ 2
    while true
        t = list[c]
        (i < t) ? (r = c - 1) : (i > t) ? (l = c + 1) : break
        (l > r) ? (c = 0; break) : (c = (l + r) ÷ 2)
    end
    c
end

"""
    selectindex(f, L, rg; base=2, alloc=1000)

Select the legal indices for a basis.

Inputs:
-------
- `f`    : Function for digits that tells whether a digits is valid.
- `L`    : Length of the system.
- `rg`   : Range of interation.
- `base` : (Optional) Base, `base=2` by default.
- `alloc`: (Optional) Pre-allocation memory for the list of indices, `alloc=1000` by default.

Outputs:
--------
- `I`: List of indices in a basis.
"""
function selectindex(
    f,
    L::Integer,
    rg::UnitRange;
    base::Integer=2,
    alloc::Integer=1000
)
    dgt = zeros(Int64, L)
    I = Int[]
    sizehint!(I, alloc)
    for i in rg
        change!(dgt, i, base=base)
        f(dgt) && append!(I, i)
    end
    I
end

"""
    selectindexnorm(f, L, rg; base=2, alloc=1000)

Select the legal indices and corresponding norm for a basis.

Inputs:
-------
- `f`    : Function for digits that tells whether a digits is valid as well as its normalization.
- `L`    : Length of the system.
- `rg`   : Range of interation.
- `base` : (Optional) Base, `base=2` by default.
- `alloc`: (Optional) Pre-allocation memory for the list of indices, `alloc=1000` by default.

Outputs:
--------
- `I`: List of indices in a basis.
- `R`: List of normalization for each states.
"""
function selectindexnorm(
    f,
    L::Integer,
    rg::UnitRange;
    base::Integer=2,
    alloc::Integer=1000
)
    dgt = zeros(Int64, L)
    I, R = Int[], Float64[]
    sizehint!(I, alloc)
    sizehint!(R, alloc)
    for i in rg
        change!(dgt, i, base=base)
        Q, N = f(dgt, i)
        if Q
            append!(I, i)
            append!(R, N)
        end
    end
    I, R
end

"""
    dividerange(maxnum::Integer, nthreads::Integer)

Divide the interation range equally according to the number of threads.

Inputs:
-------
- `maxnum`  : Maximum number of range.
- `nthreads`: Number of threads.

Outputs:
--------
- `list`: List of ranges.
"""
function dividerange(maxnum::Integer, nthreads::Integer)
    list = Vector{UnitRange{Int64}}(undef, nthreads)
    eachthreads, left = divrem(maxnum, nthreads)
    start = 1
    for i = 1:left
        stop  = start + eachthreads
        list[i] = start:stop
        start = stop+1
    end
    for i = left+1:nthreads-1
        stop  = start + eachthreads - 1
        list[i] = start:stop
        start = stop+1
    end
    list[nthreads] = start:maxnum
    list
end

"""
    selectindex_threaded(f, L, rg; base=2, alloc=1000)

Select the legal indices with multi-threads.

Inputs:
-------
- `f`    : Function for digits that tells whether a digits is valid as well as its normalization.
- `L`    : Length of the system.
- `base` : (Optional) Base, `base=2` by default.
- `alloc`: (Optional) Pre-allocation memory for the list of indices, `alloc=1000` by default.

Outputs:
--------
- `I`: List of indices in a basis.
"""
function selectindex_threaded(
    f,
    L::Integer;
    base::Integer=2,
    alloc::Integer=1000
)
    nt = Threads.nthreads()
    ni = dividerange(base^L, nt)
    nI = Vector{Vector{Int}}(undef, nt)
    Threads.@threads for ti in 1:nt
        nI[ti] = selectindex(f, L, ni[ti], base=base, alloc=alloc)
    end
    I = vcat(nI...)
    I
end

"""
    selectindexnorm_threaded(f, L, rg; base=2, alloc=1000)

Select the legal indices and corresponding norm for a basis with multi-threads.

Inputs:
-------
- `f`    : Function for digits that tells whether a digits is valid as well as its normalization.
- `L`    : Length of the system.
- `base` : (Optional) Base, `base=2` by default.
- `alloc`: (Optional) Pre-allocation memory for the list of indices, `alloc=1000` by default.

Outputs:
--------
- `I`: List of indices in a basis.
- `R`: List of normalization for each states.
"""
function selectindexnorm_threaded(
    f,
    L::Integer;
    base::Integer=2,
    alloc::Integer=1000
)
    nt = Threads.nthreads()
    ni = dividerange(base^L, nt)
    nI = Vector{Vector{Int}}(undef, nt)
    nR = Vector{Vector{Float64}}(undef, nt)
    Threads.@threads for ti in 1:nt
        nI[ti], nR[ti] = selectindexnorm(f, L, ni[ti], base=base, alloc=alloc)
    end
    I, R = vcat(nI...), vcat(nR...)
    I, R
end

# Level statistics
export gapratio, meangapratio
function gapratio(E::AbstractVector{<:Real})
    dE = diff(E)
    r = zeros(length(dE)-1)
    for i = 1:length(r)
        if dE[i] < dE[i+1]
            r[i] = dE[i]/dE[i+1]
        elseif dE[i] > dE[i+1]
            r[i] = dE[i+1]/dE[i]
        else
            r[i] = 1.0
        end
    end
    r
end

meangapratio(E::AbstractVector{<:Real}) = sum(gapratio(E)) / (length(E) - 2)


text:
JSON.lower(ms::TopNMetricSpec) = non_nothing_dict(ms)

"""
    Numeric(metric::String)

Numeric topN metric spec.
"""
struct Numeric <: TopNMetricSpec
    type::String
    metric::String
    Numeric(metric) = new("metric", metric)
end

"""
    Dimension()
    Dimension(; ordering=nothing, previousStop=nothing)

Dimension topN metric spec.

ordering and previousStop should be `String`s if provided.
"""
struct Dimension <: TopNMetricSpec
    type::String
    ordering
    previousStop
    function Dimension(; ordering=nothing, previousStop=nothing)
        ordering === nothing || (ordering = lowercase(ordering)) ∈ ["lexicographic", "alphanumeric", "numeric", "strlen"] || error("Invalid ordering")
        nothing_or_type(previousStop, String)
        new("dimension", ordering, previousStop)
    end
end

"""
    Inverted(metric::TopNMetricSpec)

Inverting topN metric spec.
"""
struct Inverted <: TopNMetricSpec
    type::String
    metric::TopNMetricSpec
    Inverted(metric) = new("inverted", metric)
end

"""
    TopN(dataSource::DataSource, intervals::Vector{Interval}, granularity::Granularity,
        dimension::DimensionSpec, threshold::Uint64, metric::TopNMetricSpec; <keyword arguments>)

Returns a Druid.TopNResult which is compatible with the Tables.jl interface.

A method with all arguments as keyword arguments is also provided.

# Arguments
- aggregations::Vector{<:Aggregator} = nothing
- postAggregations::Vector{<:PostAggregator} = nothing
- filter::Filter = nothing
- virtualColumns::Vector{<:VirtualColumn} = nothing
- context::Dict = nothing
"""
mutable struct TopN <: Query
    queryType::String
    dataSource::DataSource
    intervals::Vector{<:Interval}
    granularity::Granularity
    dimension::DimensionSpec
    threshold::UInt64
    metric::TopNMetricSpec
    aggregations
    postAggregations
    filter
    virtualColumns
    context
    function TopN(
        dataSource, intervals, granularity, dimension, threshold, metric;
        aggregations=nothing, postAggregations=nothing, filter=nothing, virtualColumns=nothing, context=nothing
    )
        nothing_or_type(aggregations, Vector{<:Aggregator})
        nothing_or_type(postAggregations, Vector{<:PostAggregator})
        nothing_or_type(filter, Filter)
        nothing_or_type(virtualColumns, Vector{<:VirtualColumn})
        nothing_or_type(context, Dict)
        if !isa(metric, Numeric) && !(isa(metric, Inverted) && isa(metric.metric, Numeric))
            if !(aggregations === nothing) || !(postAggregations === nothing)
                error("Aggregations and post aggregations can only be specified for Numeric metric specs")
            end
        end
        new(
            "topN", dataSource, intervals, granularity, dimension, threshold, metric,
            aggregations, postAggregations, filter, virtualColumns, context
        )
    end
end
TopN(
    ; dataSource, intervals, granularity, dimension, threshold, metric,
    aggregations=nothing, postAggregations=nothing, filter=nothing, virtualColumns=nothing, context=nothing
) = TopN(
    dataSource, intervals, granularity, dimension, threshold, metric;
    aggregations, postAggregations, filter, virtualColumns, context
)

function execute(client::Client, query::TopN; pretty=false)
    res = convert(Vector{Dict}, JSON.parse(execute_native_query(client, query; pretty)))
    prefixes = ["top" * string(i) * '_' for i in 1:query.threshold]
    if length(res) != 0
        names = vec([Symbol(p * k) for (p, k) ∈ Iterators.product(prefixes, keys(res[1]["result"][1]))])
        push!(names, :timestamp)
    else
        names = Symbol[]
    end
    TopNResult(names, res)
end

struct TopNResult <: QueryResult
    names::Vector{Symbol}
    inner_array::Vector{Dict}
end

names(tr::TopNResult) = getfield(tr, :names)

Base.summary(io::IO, tr::TopNResult) = print(io, string(length(tr)) * "-element " * string(typeof(tr)))
Base.getindex(tr::TopNResult, i::Int) = (i <= length(tr) || throw(BoundsError(tr, i))) && TopNRow(i, tr)

Tables.rowaccess(::TopNResult) = true
Tables.rows(tr::TopNResult) = tr

Base.eltype(::TopNResult) = TopNRow
Base.length(tr::TopNResult) = length(getfield(tr, :inner_array))
Base.iterate(tr::TopNResult, state=1) = state > length(tr) ? nothing : (TopNRow(state, tr), state + 1)

function Tables.schema(tr::TopNResult)
    if length(tr) != 0
        types = [typeof(cl) for cl in tr[1]]
    else
        types = DataType[]
    end
    Tables.Schema(names(tr), types)
end

struct TopNRow <: Tables.AbstractRow
    row::Int
    source::TopNResult
end

function Tables.getcolumn(tr::TopNRow, name::Symbol)
    if name == :timestamp
        getfield(getfield(tr, :source), :inner_array)[getfield(tr, :row)][string(name)]
    else
        pattern = r"^top(\d)_(.*)$"
        m = match(pattern, string(name))
        getfield(getfield(tr, :source), :inner_array)[getfield(tr, :row)]["result"][parse(Int, m.captures[1])][m.captures[2]]
    end
end

Tables.getcolumn(tr::TopNRow, i::Int) = Tables.getcolumn(tr, names(getfield(tr, :source))[i])
Tables.columnnames(tr::TopNRow) = names(getfield(tr, :source))


text:
module DB

using ...DBUtils
using ..Entrez
using SQLite
using MySQL
using DataStreams, DataFrames
using NullableArrays

export init_pubmed_db_mysql,
       init_pubmed_db_sqlite,
       get_value,
       all_pmids,
       get_article_mesh,
       db_insert!

get_value{T}(val::Nullable{T}) = get(val)
get_value(val)= val
get_value{T}(val_array::Array{T}) = val_array
get_value{T}(val_array::NullableArray{T, 1}) = val_array.values

function init_pubmed_db_mysql(config)

    println("Initializing MySQL Database")

    #intput dictionary must have the following keys
    if haskey(config, :host) && haskey(config, :dbname) &&
       haskey(config, :username) && haskey(config, :pswd) &&
       haskey(config, :overwrite)

       mysql_code=nothing
       try
           filename = Pkg.dir() * "/BioMedQuery/src/Entrez/create_pubmed_db.sql"
           println(filename)
           f = open(filename, "r")
           mysql_code = readall(f)
           close(f)
       catch
           error("Could not read create_entrez_db.sql")
       end

       db = DBUtils.init_mysql_database(host = config[:host], dbname =config[:dbname],
       username = config[:username], pswd= config[:pswd],
       overwrite = config[:overwrite], mysql_code = mysql_code)

       return db
   end
end

function init_pubmed_db_sqlite(config)
    println("Initializing SQLite Database")
    if haskey(config, :db_path) && haskey(config, :overwrite)
        db = init_pubmed_db_sqlite(config[:db_path], config[:overwrite])
        return db
    else
        println("Error with following configuration:")
        println(config)
        println("Must contain: db_path")
        error("Improper configuration for entrez_sqlite:init_database")
    end
end

# Creates a database with all necessary tables to store
# Entrez related searches. All tables are empty at this point
# If a database existis at the given path - an error is ruturned an the user
# is asked whether he intended to clean the existing file
function init_pubmed_db_sqlite(path::ASCIIString, overwrite=false)

    if isfile(path)
        if overwrite
            rm(path)
        else
            println("Database found. Returning existing database.")
            return SQLite.DB(path)
        end
    end

    #Create database file
    db = SQLite.DB(path)

    #Create tables to store
    SQLite.query(db, "CREATE TABLE
    article(pmid INTEGER NOT NULL PRIMARY KEY,
    title TEXT,
    pubYear INTEGER)")

    SQLite.query(db, "CREATE TABLE
    author(id INTEGER PRIMARY KEY AUTOINCREMENT,
    forename TEXT,
    lastname TEXT NOT NULL,
    CONSTRAINT unq UNIQUE(forename,  lastname) )")

    SQLite.query(db, "CREATE TABLE
    author2article(aid INTEGER, pmid INTEGER,
    FOREIGN KEY(aid) REFERENCES author(id),
    FOREIGN KEY(pmid) REFERENCES article(pmid),
    PRIMARY KEY(aid, pmid) )")

    #--------------------------
    # MeshHeading Tables
    #--------------------------

    #Descriptor
    #The id corresponds to the DUI of mesh library
    #Adding a "D" at the beginning of the id, allows for
    #lookup in the mesh browerser
    # https://www.nlm.nih.gov/mesh/MBrowser.html
    SQLite.query(db, "CREATE TABLE
    mesh_descriptor(id INTEGER NOT NULL PRIMARY KEY ,
                    name TEXT UNIQUE )")

    #Qualifier
    SQLite.query(db, "CREATE TABLE
    mesh_qualifier(id INTEGER NOT NULL PRIMARY KEY ,
                   name TEXT UNIQUE )")

    #Heading
    SQLite.query(db, "CREATE TABLE
    mesh_heading(id INTEGER PRIMARY KEY AUTOINCREMENT,
                 pmid INTEGER, did INTEGER, qid INTEGER,
                 dmjr TEXT, qmjr TEXT,
                 FOREIGN KEY(pmid) REFERENCES article(pmid),
                 FOREIGN KEY(did) REFERENCES mesh_descriptor(id),
                 FOREIGN KEY(qid) REFERENCES mesh_qualifier(id),
                 CONSTRAINT unq UNIQUE(pmid, did, qid) )")

    return db

end

"""
    all_pmids(db)
Return all PMIDs stored in the *article* table of the input database
"""
function all_pmids(db)
    query = db_query(db, "SELECT pmid FROM article;")
    return get_value(query[1])
end

"""
    all_mesh(db)
Return all PMIDs stored in the *article* table of the input database
"""
function all_mesh(db)
    query = db_query(db, "SELECT name FROM mesh_descriptor;")
    return get_value(query[1])
end

"""
    get_article_mesh(db, pmid)
Get the all mesh-descriptors associated with a give article
"""
function get_article_mesh(db, pmid::Integer)

    query_string = "SELECT md.name
                      FROM mesh_heading as mh,
                           mesh_descriptor as md
                     WHERE mh.did = md.id
                      AND mh.pmid = $pmid"
    query  = db_query(db, query_string)
    #return data array
    return get_value(query.columns[1])

end

function db_insert!(db, article::PubMedArticle, verbose=false)
    #------- PMID - TITLE - YEAR
    isnull(article.pmid) && error("NULL PMID")

    # Save article data
    insert_row!(db, "article", Dict(:pmid =>article.pmid.value,
                                    :title=>get(article.title, ""),
                                    :pubYear=>get(article.year, 0)), verbose)

    #------- AUTHORS
    for au in article.authors
        if isnull(au[:LastName])
           println("Skipping author, null field: ", au)
           continue
        end

        author_id = insert_row!(db, "author",
        Dict(:id => nothing,
             :forename => get(au[:ForeName], "Unknown"),
             :lastname => get(au[:LastName], nothing)), verbose)

        if author_id < 0
            sel = db_select(db, ["id"], "author",
            Dict(:forename => get(au[:ForeName], "Unknown"),
                 :lastname => au[:LastName].value))
            if length(sel[1]) > 0
                author_id = get_value(sel[1][1])
                if verbose
                    println("Author already in db: ", au)
                end
                insert_row!(db, "author2article",
                Dict(:aid =>author_id, :pmid => article.pmid.value), verbose)
            else
                error("Can't save nor find Author: ", au)
            end
        end
    end
end

function db_insert!(db, pmid::Int64, mesh_heading_list::MeshHeadingList, verbose=false)
    for heading in mesh_heading_list

        did_int = heading.descriptor_id.value
        descriptor_name = heading.descriptor_name.value
        dmjr = get(heading.descriptor_mjr, nothing)

        #Save Descriptor
        insert_row!(db, "mesh_descriptor",
        Dict(:id=>did_int,
             :name=>descriptor_name),
             verbose)

        if isempty(heading.qualifier_id)
            #Save Headings
            insert_row!(db, "mesh_heading",
            Dict(:id=>nothing,
                 :pmid=> pmid,
                 :did=>did_int,
                 :qid=>nothing,
                 :dmjr=>nothing, :qmjr=>nothing), verbose )
        else

            for i=1:length(heading.qualifier_id)
                qid_int = get(heading.qualifier_id[i], -1)
                qualifier_name = get(heading.qualifier_name[i], nothing)
                qmjr = get(heading.qualifier_mjr[i], nothing)

                #Save Qualifiers`
                insert_row!(db, "mesh_qualifier",
                Dict(:id=>qid_int,
                     :name=>qualifier_name),
                     verbose )

                #Save Headings
                insert_row!(db, "mesh_heading",
                Dict(:id=>nothing,
                     :pmid=> pmid,
                     :did=>did_int,
                     :qid=>qid_int,
                     :dmjr=>dmjr, :qmjr=>qmjr), verbose )
            end
        end

    end
end

end #module


text:
using QuadraticRelaxations

using Mosek, MosekTools

using Dates

const ROOT_DIR = joinpath(@__DIR__, "..")
include(joinpath(ROOT_DIR, "local_config.jl"))
include(joinpath(ROOT_DIR, "solver_factories.jl"))
include(joinpath(ROOT_DIR, "box_qp/util.jl"))
include(joinpath(ROOT_DIR, "box_qp/justitia_interface.jl"))

const TIME_LIMIT = 10 * 60.0
const NUM_THREADS = 1

const NUM_LAYERS = 3
const EPIGRAPH = false

const SHIFT = QuadraticRelaxations.SemidefiniteShift(Mosek.Optimizer)
const SOC_LOWER_BOUND = false

const FAMILIES = ("basic", "extended", "extended2")

const INSTANCE_DIR = joinpath(ROOT_DIR, "box_qp", "data")

const CONFIG = FactoryConfig(TIME_LIMIT, NUM_THREADS, CPLEX_PATH)

INSTANCES = Dict{String,BoxQPInstance}()
for family in FAMILIES, instance_name in readdir(joinpath(INSTANCE_DIR, family))
    INSTANCES["$family-$instance_name"] =
        _create_instance(family, instance_name, primal_factory(CONFIG), EPIGRAPH)
end

APPROACHES = Dict{String,BoxQPApproach}()
for (t_name, transformation) in [
    "direct" => identity,
    "direct-shift" =>
        reformulate_with(QuadraticRelaxations.ExactFormulation(), SHIFT, SOC_LOWER_BOUND),
]
    APPROACHES[string("grb3", "-", t_name)] =
        BoxQPApproach(gurobi_direct_mipfocus_3_factory, transformation)
end

for (r_name, relaxation) in [
    "nn" => QuadraticRelaxations.NeuralNetRelaxation(NUM_LAYERS),
    "cda" => QuadraticRelaxations.CDARelaxation(NUM_LAYERS),
    "standard-nmdt" => QuadraticRelaxations.StandardNMDT(NUM_LAYERS),
    "tightened-nmdt" => QuadraticRelaxations.TightenedNMDT(NUM_LAYERS),
]
    APPROACHES[string("grb3", "-", r_name)] = BoxQPApproach(
        gurobi_direct_mipfocus_3_factory,
        reformulate_with(relaxation, SHIFT, SOC_LOWER_BOUND),
    )
end


text:
export global_queue

# Context management and global state
function global_queue(dev::MtlDevice)
    get!(task_local_storage(), (:MtlCommandQueue, dev)) do
        MtlCommandQueue(dev)
    end
end

function MTL.device()
    get!(task_local_storage(), :MtlDevice) do
        MtlDevice(1)
    end
end


text:
using BinaryBuilder

name = "Zstd"
version = v"1.4.5"

sources = [
    ArchiveSource("https://github.com/facebook/zstd/releases/download/v$version/zstd-$version.tar.gz",
                  "98e91c7c6bf162bf90e4e70fdbc41a8188b9fa8de5ad840c401198014406ce9e"),
]

script = raw"""
cd ${WORKSPACE}/srcdir/zstd-*/
mkdir build-zstd && cd build-zstd
if [[ "${target}" == i686-linux-musl ]]; then
    # We can't run executables for i686-linux-musl in the BB environment
    sed -i "s/needs_exe_wrapper = false/needs_exe_wrapper = true/" "${MESON_TARGET_TOOLCHAIN}"
elif [[ "${target}" == aarch64-linux-gnu ]]; then
    # Work around https://github.com/facebook/zstd/issues/1872
    #
    #   ../programs/util.c:76:29: error: ‘stat_t’ has no member named ‘st_mtim’
    #            timebuf[1] = statbuf->st_mtim;
    #                                ^
    sed -i "s/c_args = \[\]/c_args = ['-D_POSIX_C_SOURCE']/" "${MESON_TARGET_TOOLCHAIN}"
fi
meson --cross-file="${MESON_TARGET_TOOLCHAIN}" ../build/meson/
ninja -j${nproc}
ninja install
"""

platforms = supported_platforms()

products = [
    LibraryProduct("libzstd", :libzstd),
    ExecutableProduct("zstd", :zstd),
    ExecutableProduct("zstdmt", :zstdmt),
]

dependencies = Dependency[]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using CSV
using DataFrames
using Plots
using JLD2
using DiscretisedFluidQueues
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/plot_utils.jl")
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/error_metrics.jl")

@halfwidth_plot_defaults()

plotlyjs()
function p()
ks_data = CSV.read((@__DIR__)*"/../data/meshs_ks_func_count_4.csv",DataFrame)
ks_data.DG_lin_lim = CSV.read((@__DIR__)*"/../data2/meshs_ks_func_count_4.csv",DataFrame).DG_limiter

plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot,:dash]
markerstyles_vec = [:cross,:diamond,:circle,:dot,:diamond]
for (c,col) in enumerate(names(ks_data))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),ks_data[:,col],
        label=col,linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2)
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="KS error - CDF")
# plot!(legend=(0.225,0.625))
plot!()
@add_lines!(ks_data,("Unif","QBDRAP"),@__DIR__)
savefig((@__DIR__)*"/meshs_ks_error_formatted.pdf")
end
p()

function p()
l2_pdf_data = CSV.read((@__DIR__)*"/../data/meshs_l2_pdf_func_count_4.csv",DataFrame)
l2_pdf_data.DG_lin_lim = CSV.read((@__DIR__)*"/../data2/meshs_l2_pdf_func_count_4.csv",DataFrame).DG_limiter

plot()
linestyles_vec = [:solid,:dash,:dashdot, :dot,:dash]
markerstyles_vec = [:cross,:diamond,:circle,:cross,:diamond]
for (c,col) in enumerate(names(l2_pdf_data))
    os = collect(1:2:21)
    (c==5)&&(os.+=1)
    plot!(log10.(os),l2_pdf_data[:,col],
        label=col,linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2)
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="L² error - PDF")
# plot!(legend=(0.225,0.625))
plot!()
@add_lines!(l2_pdf_data,("Unif","QBDRAP"),@__DIR__)
savefig((@__DIR__)*"/meshs_l2_pdf_error_formatted.pdf")
end
p()

function p()
@fullwidth_plot_defaults()
dir = @__DIR__()
os = 1:2:7
plot(layout=(5,length(os)))
linestyles_vec = [:solid,:solid,:solid,:solid]
markerstyles_vec = [:cross,:diamond,:circle]
xs = range(2+eps(),7-eps(),length=200)
for (c,o) in enumerate(os)
    model = BoundedFluidQueue(fill(0.0,2,2),[1.0; 0.0],10.0)

    mesh = DGMesh(0.0:10.0,o)
    dq = DiscretisedFluidQueue(model,mesh)
    jldopen(dir*"/../../../coeffs_matrix.jld2") do f
        coeffs = f["coeff_matrix"][4,c,1]
        dt = SFMDistribution(coeffs,dq)
        plot!(xs,pdf(dt).(xs,1),subplot=c,label=(c==1 ? "DG" : false),
            linestyle=linestyles_vec[1],ylims = (-0.1,1.1),  yticks=false,
            # marker=markerstyles_vec[1],
            linewidth=2,xticks=false)
    end
    plot!(subplot=c, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylabel=(c==1 ? "DG" : false), grid=false,
        title = string("Dim. ",o))

    mesh = DGMesh(0.0:(1.0./o):10.0,1)
    dq = DiscretisedFluidQueue(model,mesh)
    jldopen(dir*"/../../../coeffs_matrix.jld2") do f
        coeffs = f["coeff_matrix"][4,c,2]
        dt = SFMDistribution(coeffs,dq)
        plot!(xs,pdf(dt).(xs,1),subplot=c+4,label=(c==1 ? "Unif." : false),
            linestyle=linestyles_vec[2],ylims = (-0.1,1.1),  yticks=false,
            # marker=markerstyles_vec[2],
            linewidth=2,xticks=false)
    end
    (c==1)&&plot!(subplot=c+4, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylabel=(c==1 ? "Unif." : false), grid=false)

    mesh = DGMesh(0.0:10.0,o)
    dq = DiscretisedFluidQueue(model,mesh)
    jldopen(dir*"/../../../coeffs_matrix.jld2") do f
        coeffs = f["coeff_matrix"][4,c,3]
        dt = SFMDistribution(coeffs,dq)
        plot!(xs,pdf(dt).(xs,1),subplot=c+8,label=(c==1 ? "DG (limiter)" : false),
            linestyle=linestyles_vec[3],ylims = (-0.1,1.1),  yticks=false,
            # marker=markerstyles_vec[3],
            linewidth=2,xticks=false)
    end
    (c==1)&&plot!(subplot=c+8, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylims = (-0.1,1.1),
        ylabel=(c==1 ? "DG-limit" : false), grid=false)

    mesh = FRAPMesh(0.0:10.0,o)
    dq = DiscretisedFluidQueue(model,mesh)
    jldopen(dir*"/../../../coeffs_matrix.jld2") do f
        coeffs = f["coeff_matrix"][4,c,4]
        dt = SFMDistribution(coeffs,dq)
        plot!(xs,pdf(dt).(xs,1),subplot=c+12,label=(c==1 ? "Q-R" : false),
            linestyle=linestyles_vec[3],ylims = (-0.1,1.1), yticks=false,
            # marker=markerstyles_vec[3],
            linewidth=2,xticks=false, grid=false)
    end
    (c==1)&&plot!(subplot=c+12, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylims = (-0.1,1.1),
        ylabel=(c==1 ? "Q-R" : false), grid=false)
    # plot!(xs,x->2*x*(x<0.5) + Float64(x>=0.5),subplot=c,label=(c==1 ? "truth" : false),
    #     linestyle=:dot,
    #     # marker=markerstyles_vec[3],
    #     color=:grey,
    #     linewidth=2)

    # plot!(subplot=c, yticks=(c==1 ? (0.1:0.3:1.2) : false),
    #     ylims = (-0.1,1.1), xticks = 0.0:3:10,
    #     ylabel=(c==1 ? "probability" : false), grid=false,
    #     title = string("Dim. ",o),
    #     legend=:outertopright)
    plot!(legend=false)

    model = BoundedFluidQueue([0 0; 1.0 -1.0],[1.0; 0.0],10.0)

    mesh = DGMesh(0.0:(2/(o+1)):10.0,2)
    dq = DiscretisedFluidQueue(model,mesh)
    thepdf0(z,i) = 1.0*(z<1.0)*(i==1)*(z>=0.0)*(z<6.0)
    d0 = SFMDistribution(thepdf0,dq)
    dt = integrate_time(d0,build_full_generator(dq),4.0, StableRK4(0.005); limiter=GeneralisedMUSCL)

    plot!(xs,pdf(dt).(xs,1),subplot=c+16,label=(c==1 ? "DG lin lim" : false),
        linestyle=linestyles_vec[3],ylims = (-0.1,1.1), yticks=false,
        # marker=markerstyles_vec[3],
        linewidth=2,xticks=0:2:10, grid=false)
    (c==1)&&plot!(subplot=c+16, yticks=(c==1 ? (0.1:0.4:1.2) : false),
        ylims = (-0.1,1.1),
        ylabel=(c==1 ? "DG lin lim" : false), grid=false)
    plot!(legend=false)
end
plot!()
savefig((@__DIR__)*"/pdfs_formatted.pdf")
end
p()

function p()
L1_cell_probs_errors = DataFrame(DG = Float64[], Unif = Float64[], DG_limiter = Float64[], QBDRAP = Float64[], DG_lin_lim = Float64[])
truth = zeros(10)
truth[5] = 1.0
truth_pm = zeros(3)
for (c_o,o) in enumerate(1:2:21)
    row = zeros(5)
    for (c_m,m) in enumerate(names(L1_cell_probs_errors))
        (c_m<5)&&jldopen((@__DIR__)*"/../../../coeffs_matrix.jld2") do f
            model = BoundedFluidQueue(fill(0.0,2,2),[1.0; 0.0],10.0)

            mesh = (c_m<4) ? DGMesh(0.0:10.0,o) : FRAPMesh(0.0:10.0,o)
            dq = DiscretisedFluidQueue(model,mesh)

            coeffs = f["coeff_matrix"][4,c_o,c_m]
            dt = SFMDistribution(coeffs,dq)

            x_vals = range(0.5,9.5;length=10)

            approx = sum(cell_probs(dt).(x_vals,(1:2)'),dims=2)
            approx_pm = [coeffs[1:N₋(model)];coeffs[end-N₊(model)+1:end]]

            row[c_m] = log10(sum(abs.(approx-truth)) + sum(abs.(approx_pm-truth_pm)))
        end
        (c_m==5)&&begin
            model = BoundedFluidQueue([0 0; 1.0 -1.0],[1.0; 0.0],10.0)

            mesh = DGMesh(0.0:(2/(o+1)):10.0,2)
            dq = DiscretisedFluidQueue(model,mesh)
            thepdf0(z,i) = 1.0*(z<1.0)*(i==1)*(z>=0.0)*(z<6.0)
            d0 = SFMDistribution(thepdf0,dq)
            dt = integrate_time(d0,build_full_generator(dq),4.0, StableRK4(0.005); limiter=GeneralisedMUSCL)
            coeffs = dt.coeffs

            x_vals = (mesh.nodes[1:end-1]+mesh.nodes[2:end])./2

            approx = sum(reshape(cell_probs(dt).(x_vals,1),Int((o+1)/2),10),dims=1)[:]
            approx_pm = [coeffs[1:N₋(model)];coeffs[end-N₊(model)+1:end]]

            row[c_m] = log10(sum(abs.(approx-truth)) + sum(abs.(approx_pm-truth_pm)))
        end
    end
    push!(L1_cell_probs_errors,row)
end
@halfwidth_plot_defaults()
file = (@__DIR__)*"/../data/meshs_l1_cell_"*"func_count_4"
CSV.write(file*".csv",L1_cell_probs_errors)
q = plot()
linestyles_vec = [:solid,:dash,:dashdot,:dot,:dash]
markerstyles_vec = [:cross,:diamond,:circle,:dot,:diamond]
for (c,reconstruction) in enumerate(names(L1_cell_probs_errors))
    os = collect(1:2:21)
    (c==5)&&(os .+= 1)
    plot!(log10.(os),L1_cell_probs_errors[:,reconstruction],label=reconstruction,
        linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],
        linewidth=2)
end
plot!(xlabel="Dimension", ylabel="Error",
    title="Error between cell probabilities",)
    # legend=(0.225,0.625))
error_ticks!(plot!())
@add_lines!(L1_cell_probs_errors,("Unif","QBDRAP"),@__DIR__)
savefig((@__DIR__)*"/L1_cell_probs.pdf")
end
p()
# os = 1:2:7
# plot(layout=(1,length(os)))
# linestyles_vec = [:solid,:dash,:dashdot]
# markerstyles_vec = [:cross,:diamond,:circle]
# xs = range(eps(),1-eps(),length=200)
# for (c,o) in enumerate(os)
#     model = BoundedFluidQueue(fill(0.0,1,1),[1.0],1.0)

#     mesh = DGMesh(0.0:1.0,o)
#     dq = DiscretisedFluidQueue(model,mesh)
#     d0 = SFMDistribution((x,i)->2*Float64(x<=0.5),dq,TrapezoidRule;fun_evals=10_001)
#     plot!(xs,pdf(d0).(xs,1),subplot=c,label=(c==1 ? "DG" : false),
#         linestyle=linestyles_vec[1],
#         # marker=markerstyles_vec[1],
#         linewidth=2)

#     mesh = DGMesh(0.0:(1.0./o):1.0,1)
#     dq = DiscretisedFluidQueue(model,mesh)
#     d0 = SFMDistribution((x,i)->2*Float64(x<=0.5),dq,TrapezoidRule;fun_evals=10_001)
#     plot!(xs,pdf(d0).(xs,1),subplot=c,label=(c==1 ? "Unif." : false),
#         linestyle=linestyles_vec[2],
#         # marker=markerstyles_vec[2],
#         linewidth=2)

#     mesh = FRAPMesh(0.0:1.0,o)
#     dq = DiscretisedFluidQueue(model,mesh)
#     d0 = SFMDistribution((x,i)->2*Float64(x<=0.5),dq,TrapezoidRule;fun_evals=10_001)
#     plot!(xs,pdf(d0).(xs,1),subplot=c,label=(c==1 ? "QBD-RAP" : false),
#         linestyle=linestyles_vec[3],
#         # marker=markerstyles_vec[3],
#         linewidth=2)

#     plot!(xs,x->2*Float64(x<=0.5),subplot=c,label=(c==1 ? "truth" : false),
#         linestyle=:dot,
#         # marker=markerstyles_vec[3],
#         color=:grey,
#         linewidth=2)

#     plot!(subplot=c, yticks=(c==1 ? (-2:2:6) : false),
#         ylims = (-3,7.2), xticks = 0.25:0.5:0.75,
#         ylabel=(c==1 ? "probability" : false), grid=false,
#         title = string("Dim. ",o),
#         legend=(0.1,0.8))
# end
# plot!()
# savefig((@__DIR__)*"/pdfs_formatted.pdf")


text:
using Test
import CryptoGroups
import CryptoGroups: PGroup, specialize, Enc, Dec, <|, generator, PGroup, ECGroup

function elgamal_test(g)

    sk = 5
    pk = g^sk
    r = 3
    m = g^5
    r2 = 6

    enc = Enc(pk, g)
    dec = Dec(sk)

    @test dec(enc(m, r)) == m
    @test enc(enc(m, r), r2) == enc(m, r + r2)

    ### Shuffle generation

    sk = 5
    pk = g^sk

    enc = Enc(pk, g)

    m_vec = [g^4, g^2, g^3]
    e_vec = enc.(m_vec, [2, 3, 7])

    ### The shuffling
    r_vec = Int[4, 2, 3]

    e_enc = enc.(e_vec, r_vec)
    ψ = sortperm(e_enc)
    sort!(e_enc)

    @test sort(dec.(e_enc)) == sort(m_vec)

    #m_vec = [g, g^2, g^3]
    #e_vec = enc.(m_vec, 1)
end

let
    q = 11
    p = 2*q + 1

    G = PGroup(p, q)
    g = G(3)

    elgamal_test(g)
end

import CryptoGroups
import CryptoGroups: specialize, <|, generator, PGroup, ECGroup, Specs

let
    spec = Specs.MODP_1024
    G = specialize(PGroup, spec)
    g = G <| generator(spec)

    elgamal_test(g)
end

let
    spec = Specs.Curve_P_256

    G = specialize(ECGroup, spec; name = :P_256)
    g = G <| generator(spec)

    elgamal_test(g)
end

let
    spec = Specs.Curve_B_163_PB
    G = specialize(ECGroup, spec; name = :B_163_PB)
    g = G <| generator(spec)

    elgamal_test(g)
end


text:
import Base.push!

# Methods for altering arrays containing nodes (children) depending on the chosen Loss-function
# FNN, CCM & Prediction

"""
    push!(children::Union{Array{Node,1},Nothing}, n::EmbeddingPars, Γ::AbstractLoss, current_node::AbstractTreeElement)

Adds new `children`/nodes to the tree below `current_node` with embedding parameters `n`, according to the loss `Γ`.
"""
function push!(children::Array{Node,1}, n::EmbeddingPars, Γ::AbstractLoss, current_node::AbstractTreeElement)
    Base.push!(children, Node(n, [get_τs(current_node); τ(n)], [get_ts(current_node); t(n)], nothing))
end

# L-function is computed as increments of the last value, that's why here it has to be added to the total in this function
function push!(children::Array{Node,1}, n::EmbeddingPars, Γ::L_statistic, current_node::Root)
    Base.push!(children, Node(EmbeddingPars(τ=τ(n),t=t(n),L=(current_node.Lmin+L(n)),temp=temp(n)), [get_τs(current_node); τ(n)], [get_ts(current_node); t(n)], nothing))
end

function push!(children::Array{Node,1}, n::EmbeddingPars, Γ::L_statistic, current_node::Node)
    Base.push!(children, Node(EmbeddingPars(τ=τ(n),t=t(n),L=(L(current_node)+L(n)), temp=temp(n)), [get_τs(current_node); τ(n)], [get_ts(current_node); t(n)], nothing))
end

"""
    init_embedding_params(Γ::AbstractLoss, N::Int)

Return the initial embedding parameters and loss function value, based on the chosen loss function. Every new loss should get a new function, otherwise the default (0, 1, 99999, nothing) is returned.
"""
function init_embedding_params(Γ::AbstractLoss, N::Int)
    return [EmbeddingPars(τ=0, t=1, L=99999f0)]
end
function init_embedding_params(Γ::FNN_statistic, N::Int)
    return [EmbeddingPars(τ=0, t=1, L=1f0)]
end
function init_embedding_params(Γ::L_statistic, N::Int)
    return [EmbeddingPars(τ=0, t=1, L=0f0)]
end
function init_embedding_params(Γ::CCM_ρ, N::Int)
    return [EmbeddingPars(τ=0, t=1, L=0f0)]
end

"""
    get_potential_delays(optimalg::AbstractMCDTSOptimGoal, Ys::Dataset, τs, w::Int, τ_vals,
                    ts_vals, L_old ; kwargs...]) → τ_pot, ts_pot, L_pot, flag, temps

    Computes an vector of potential embedding parameters: the potential delay `τ_pot` and
    time series values `ts_pot`, which would each result in a potential Loss-statistic value
    `L_pot`, by using an embedding method specified in `optimalg` [^Kraemer2022] (see [`MCDTSOptimGoal`](@ref))
    and for a range of possible delay values `τs`. The input dataset `Ys` can be
    multivariate. `w` is the Theiler window (neighbors in time with index `w` close
    to the point, that are excluded from being true neighbors. `w=0` means to
    exclude only the point itself, and no temporal neighbors. In case of multivariate
    time series input choose `w` as the maximum of all `wᵢ's`. `τ_vals` and `ts_vals`
    describe the embedding up to the current embedding cycle.

    ## Keyword arguments
    * See [`mcdts_embedding`](@ref) for a list of all keywords.
"""
function get_potential_delays(optimalg::AbstractMCDTSOptimGoal, Yss::Dataset{D, T},
                τs, w::Int, τ_vals, ts_vals, L_old; kwargs...) where {D, T}

    Ys = DelayEmbeddings.standardize(Yss)

    # compute actual embedding trajectory Y_act
    if typeof(optimalg.Γ) == Prediction_error
        Y_act = genembed(Ys, τ_vals .* (-1), ts_vals) # ensure causality for forecasts
    else
        Y_act = genembed(Ys, τ_vals, ts_vals)
    end

    # compute potential delay values with corresponding time series values and
    # Loss-values
    embedding_pars = embedding_cycle(optimalg, Y_act, Ys, τs, w, τ_vals, ts_vals; kwargs...)

    if isempty(embedding_pars)
        flag = true
        return EmbeddingPars[], flag
    end

    embedding_pars, converge = get_embedding_params_according_to_loss(optimalg.Γ,
                                            embedding_pars, L_old)

    return embedding_pars, converge
end

"""
    get_embedding_params_according_to_loss(Γ::AbstractLoss, τ_pot, ts_popt, L_pot, L_old)

    Helper function for [`get_potential_delays`](@ref). Computes the potential
    delay-, time series- and according Loss-values with respect to the actual loss
    in the current embedding cycle.
"""
function get_embedding_params_according_to_loss(Γ::AbstractLoss, embedding_pars::Vector{EmbeddingPars}, L_old)
    threshold = Γ.threshold
    L_pot = L.(embedding_pars)
    if (minimum(L_pot) ≥ L_old)
        return EmbeddingPars[], true
    elseif (minimum(L_pot) ≤ threshold)
        ind = L_pot .< L_old
        return embedding_pars[ind], true
    else
        ind = L_pot .< L_old
        return embedding_pars[ind], false
    end
end

"""
    Perform a potential embedding cycle from the multi- or univariate Dataset `Ys`.
    Return the possible delays `τ_pot`, the associated time series `ts_pot` and
    the corresponding L-statistic-values, `L_pot` for each peak, i.e. for each
    (`τ_pot`, `ts_pot`) pair. If `FNN=true`, `L_pot` stores the corresponding
    fnn-statistic-values.
"""
function embedding_cycle(optimalg::AbstractMCDTSOptimGoal, Y_act, Ys, τs,
                                                    w, τ_vals, ts_vals; kwargs...)

    # Compute Delay-pre-selection method according to `optimalg.Λ`
    delay_pre_selection_statistic = get_delay_statistic(optimalg.Λ, Ys, τs, w, τ_vals, ts_vals; kwargs... )

    # update τ_vals, ts_vals, Ls, ε★s
    embedding_params = pick_possible_embedding_params(optimalg.Γ, optimalg.Λ, delay_pre_selection_statistic, Y_act, Ys, τs, w, τ_vals, ts_vals; kwargs...)

    return embedding_params
end

"""
    Compute all possible τ-values (and according time series numbers) and their
    corresponding Loss-statistics for the input delay_pre_selection_statistic `dps`.
"""
function pick_possible_embedding_params(Γ::AbstractLoss, Λ::AbstractDelayPreselection, dps, Y_act::Dataset{D, T}, Ys, τs, w::Int, τ_vals, ts_vals; kwargs...) where {D, T}

    embedding_pars = EmbeddingPars[]
    for ts = 1:size(Ys,2)
        # compute loss and its corresponding index w.r.t `delay_pre_selection_statistic`

        # zero-padding of dps in order to also cover τ=0 (important for the multivariate case)
        L_trials, max_idx, temp = compute_loss(Γ, Λ, vec([0; dps[:,ts]]), Y_act, Ys, τs, w, ts, τ_vals, ts_vals; kwargs...)
        if isempty(max_idx)
            tt = max_idx
        else
            tt = τs[max_idx.-1]
            if typeof(tt)==Int
                tt = [tt]
            end
        end
        for i_trial ∈ 1:length(L_trials)
            embedding_par = EmbeddingPars(τ=tt[i_trial], t=ts, L=L_trials[i_trial], temp=temp)
            Base.push!(embedding_pars, embedding_par)
        end

    end

    return embedding_pars
end

"""
    Compute the loss of a given delay-preselection statistic `dps` and the loss
    determined by `optimalg.Γ`.
"""
compute_loss


text:
#! /usr/bin/env julia

using Benchmark
using MeshSlicer
using DataFrames

outfile = Pkg.dir("MeshSlicer")*"/perf/slice1_benchmarks.csv"
example_dir = Pkg.dir("MeshSlicer")*"/examples/"

function slice1()
    mesh = PolygonMesh(example_dir*"cell_lamp.stl");
    layers = 1000;

    # find slice heights
    step = (mesh.bounds.zmax - mesh.bounds.zmin)/layers

    heights = [mesh.bounds.zmin:step:mesh.bounds.zmax]
    a = MeshSlice(mesh,heights)
end

out = benchmark(slice1, "Slicing", "cell lamp", 10)

if isfile(outfile)
    previous = readtable(outfile)
    out = vcat(previous, out)
end
writetable(outfile, out)


text:
include("../src/dependencies.jl")
# include("src/dependencies.jl") # This line is just here to test the plot function call at the end of the file from the Julia REPL.

lambda_range = [i for i in 0.2:0.2:1.2]

# Simulation recorded data
sojurn_times = Float64[]
system_job_totals, total_orbiting_jobs = Int64[], Int64[]

# Processed simulation data for project plots.
mean_system_job_totals = zeros(Float64, length(lambda_range))
proportion_orbiting_jobs = zeros(Float64, length(lambda_range))
processed_sojurn_times = []

# function plot_data_callbacks(time::Float64, state::NetworkState, event::TimedEvent)
#     state.number_in_system_decreased && push!(sojurn_times, time - event.arrival_time)
#     push!(total_orbiting_jobs, state.orbiting_jobs)
#     push!(system_job_totals, state.number_in_system)
#     if state.number_in_system_decreased
#         println(state.queues)
#         println(state.number_in_system, "\n")
#     end
#     return nothing
# end

prev_time = 0
prev_number_in_system = 0
number_in_system_integral = 0.0

function plot_data_callbacks(time::Float64, state::NetworkState, event::TimedEvent)
    state.number_in_system_decreased && push!(sojurn_times, time - event.arrival_time)
    global number_in_system_integral += (time - prev_time) * prev_number_in_system
    global prev_number_in_system = state.number_in_system
    global prev_time = time
    return nothing
end

"""
Call simulate function below for lambdas in lambda_range
"""

current_scenario = NetworkParameters(  L=5,
                                gamma_shape = 3.0,
                                λ = 1.0,
                                η = 4.0,
                                μ_vector = collect(5:-1:1),
                                P = [0   0.5 0.5 0   0;
                                     0   0   0   1   0;
                                     0   0   0   0   1;
                                     0.5 0   0   0   0;
                                     0.2 0.2 0.2 0.2 0.2],
                                Q = [0 0 0 0 0;
                                     1 0 0 0 0;
                                     1 0 0 0 0;
                                     1 0 0 0 0;
                                     1 0 0 0 0],
                                p_e = AnalyticWeights([0.2, 0.2, 0, 0, 0.6]),
                                K = [-1, -1, 10, 10, 10])

Random.seed!(0)
@time comp = 0
@time simulate(NetworkState(12, zeros(Int8, current_scenario.L), 0), TimedEvent(ArrivalEvent(),0.0, 0, 0.0), current_scenario, max_time = 100.0, callback = plot_data_callbacks)
# fill(Int8(4),3)
# scenario_plots(processed_sojurn_times, mean_system_job_totals, proportion_orbiting_jobs, 999, lambda_range)


text:
### General constraint Jacobians match known solutions
pI = 6
n = 3
m = 3

function cI(cdot,x,u)
    cdot[1] = x[1]^3 + u[1]^2
    cdot[2] = x[2]*u[2]
    cdot[3] = x[3]
    cdot[4] = u[1]^2
    cdot[5] = u[2]^3
    cdot[6] = u[3]
end

c_jac = TrajectoryOptimization.generate_general_constraint_jacobian(cI,pI,n,m)

x = [1;2;3]
u = [4;5;6]

A = zeros(6,3)
B = zeros(6,3)

# cx, cu = c_jac(x,u)
c_jac(A,B,x,u)
cx_known = [3 0 0; 0 5 0; 0 0 1; 0 0 0; 0 0 0; 0 0 0]
cu_known = [8 0 0; 0 2 0; 0 0 0; 8 0 0; 0 75 0; 0 0 1]

@test all(A .== cx_known)
@test all(B .== cu_known)
###

# ### Custom equality constraint on quadrotor quaternion state: sqrt(q1^2 + q2^2 + q3^2 + q4^2) == 1
# opts = TrajectoryOptimization.SolverOptions()
# opts.verbose = false
# opts.constraint_tolerance = 1e-3
# opts.cost_tolerance_intermediate = 1e-3
# opts.cost_tolerance = 1e-3
# ######################
#
# ### Set up model, objective, solver ###
# # Model
# n = 13 # states (quadrotor w/ quaternions)
# m = 4 # controls
# dt = 0.05
# model, obj_uncon = TrajectoryOptimization.Dynamics.quadrotor
#
# # -control limits
# u_min = -50.0
# u_max = 50.0
#
# # -constraint that quaternion should be unit
# function cE(cdot,x,u)
#     cdot[1] = sqrt(x[4]^2 + x[5]^2 + x[6]^2 + x[7]^2) - 1.0
# end
#
# obj_con = TrajectoryOptimization.ConstrainedObjective(obj_uncon, u_min=u_min, u_max=u_max, cE=cE)#,cI=cI)
# # Solver
# # - Initial control and state trajectories
# solver = TrajectoryOptimization.Solver(model,obj_con,integration=:rk4,dt=dt,opts=opts)
# U = 10.0*ones(solver.model.m, solver.N)
#
# ##################
#
# ### Solve ###
# results, stats = TrajectoryOptimization.solve(solver,U)
# #############
# @test stats["c_max"][end] < opts.constraint_tolerance


text:
import MuJoCo
using MuJoCo

include("simenvi.jl")

function invstep(sim::MJSimEnv, x_prev, x_curr, x_next)

    # separate out x`
    q_prev = x_prev[...,:sim.data.nq]
    q_curr = x_curr[...,:sim.data.nq]
    q_next = x_next[...,:sim.data.nq]

    # compute time derivatives
    qpos = q_curr
    qvel = (q_next - q_curr) / sim.dt
    qacc = (q_prev - 2.0 * q_curr + q_next) / (sim.dt * sim.dt)

    assert qpos.shape[-1] == sim.model.nq
    assert qvel.shape[-1] == sim.model.nv
    assert qacc.shape[-1] == sim.model.nv
    assert qpos.shape[:-1] == qvel.shape[:-1] and qvel.shape[:-1] == qacc.shape[:-1], \
        'All inputs must have the same batch size, but:\nQPOS=%s\nQVEL=%s\nU=%s' \
        % (qpos.shape[:-1], qvel.shape[:-1], qacc.shape[:-1])
    batch_shape = list(qpos.shape[:-1])

    state = np.empty(batch_shape, dtype=sim.dtype)

    for flat_i in range(np.prod(batch_shape)):
        i = np.unravel_index(flat_i, batch_shape)
        sim.model.data.qpos = qpos[i]
        sim.model.data.qvel = qvel[i]
        sim.model.data.qacc = qacc[i]

        # perform inverse dynamics step
        #mjlib.mj_inverse(sim.model.ptr, sim.model.data.ptr)
        # recreate mj_inverse:
        mjlib.mj_invPosition(sim.model.ptr, sim.model.data.ptr)
        mjlib.mj_invVelocity(sim.model.ptr, sim.model.data.ptr)
        #mjlib.mj_invConstraint(sim.model.ptr, sim.model.data.ptr)
        qfrc_inverse = sim.model.data.qfrc_inverse
        mjlib.mj_rne(sim.model.ptr, sim.model.data.ptr, 1, qfrc_inverse.ctypes.data_as(ctypes.POINTER(ctypes.c_double)))
        # # override forces
        # efc_force = sim.model.data.efc_force
        # qfrc_constraint = sim.model.data.qfrc_constraint
        # mjlib.mj_mulJacTVec(sim.model.ptr, sim.model.data.ptr, \
        #     qfrc_constraint.ctypes.data_as(ctypes.POINTER(ctypes.c_double)), \
        #     efc_force.ctypes.data_as(ctypes.POINTER(ctypes.c_double)))
        # sim.model.data.qfrc_constraint = qfrc_constraint

        sim.model.data.qfrc_inverse = qfrc_inverse \
            + sim.model.dof_armature * sim.model.data.qacc \
            - sim.model.data.qfrc_passive# - sim.model.data.qfrc_constraint

        # record results
        for result_field, source_field, source_idx in sim.fields:
            val = getattr(sim.model.data, source_field)
            if source_idx is not None: val = val[source_idx]
            state[result_field][i] = val.flatten()

    # convert the numpy structured array type back into a dict so it has .items() etc.
    # this is fine since it just creates views and does not copy.
    state_dict = {k: state[k] for k in state.dtype.names}
    # add custom fields to dictionary
    state_dict["frc_root"] = state_dict["qfrc_inverse"][...,sim.unactuated_dofs()]
    state_dict["frc_actuator"] = state_dict["qfrc_inverse"][...,sim.actuated_dofs()]
    # penetration
    state_dict["efc_penetration"] = -np.minimum(state_dict["efc_pos"],0.0)

    return state_dict
end

function dynamics_linearization(inverse_dynamics, x_means, eps)
    f_x = []
    for i=1:length(x_means)
        xi_mean = x_means[i]
        T = size(xi_mean)[1]
        dimx = size(xi_mean)[2]

        #calculate finite difference samples
        x_samples = Vector(x_samples)
        for j=1:length(x_means)
            x_samples[j] = repeat(x_means[j], outer=[1, 1, dimx * 2])
        end
        for k=1:dimx
            x_samples[i][k, :, k] -= eps
            x_samples[i][k, :, k + dimx] += eps
        end

        #compute inverse_dynamics
        f_x_samples = inverse_dynamics(x_samples)

        #calculate finite differences from samples
        f_xi = {}
        keys = keys(f_x_samples)
        for j=1:length(f_x_samples)
            v_samples = get(f_x_samples, keys[j], nothing)
            dimv = shape(v_samples)[3]
            v_x = zeros([dimx, dimv, T])
            for k=1:dimx
                v_x[k, :, :] = (v_samples[:, :, k + dimx] - v_samples[:. :, k]) / 2 / eps
            end
            push(f_xi, keys[j], v_x)
        end
        push(f_x, f_xi)
    end

    return (inverse_dynamics(x_means), f_x)
end

function create_cost_quadratic_element(cost)
    w = cost["w"]
    r = get(cost, "r", nothing)
    c = w * VectorVectorProduct(r, r) / 2
    c_x = []
    c_xx = []
    r_x = get(cost, "r_x", nothing)
    for i=1:length(r_x)
        c_xi = w * MatrixVectorProduct(r_x[i], r)
        push!(c_x, c_xi)

        c_xix = []
        for j=1:length(r_x)
            c_xixj = w * MatrixMatrixProduct(r_x[i], r_x[j])
            push!(c_xix, c_xixj)
        end
        push!(c_xx, c_xix)
    end
    return (c, c_x, c_xx)
end

function create_cost_quadratic(costlist)
    C, C_x, C_xx = create_cost_quadratic_element(costlist[1])
    for i=2:length(costlist)
        c, c_x, c_xx = create_cost_quadratic_element(costlist[i])
        C += c
        for i=1:length(c_x)
            C_x[i] += c_x[i]

            for j=1:length(c_xx[i])
                C_xx[i][j] += c_xx[i][j]
            end
        end
    end

    return (C, C_x, C_xx)
end

function direct_solve(c, c_x, c_xx, lambda)
    T = length(c)

    #assemble sparse trajectory Hessian and dense gradient
    l_x = repeat(nothing, outer=[T])
    l_xx = repeat(nothing, outer=[T, T])

    for t=1:T
        for i=1:3
            ti = t + i - 1
            if (ti < 0 || ti > (T-1))
                continue
            end
            if l_x[ti] == nothing
                l_x[ti] = 0.0
            end
            l_x[ti] += c_x[i][:, t]
            for j=1:3
                tj = t + j - 1
                if (tj < 0 || tj > (T-1))
                    continue
                end
                if l_xx[tj, ti] == nothing
                    l_xx[tj, ti] = 0.0
                end
                l_xx[tj, ti] += c_xx[i][j][:,:, t]
            end
        end
    end

end

function direct_trajectory_optimization(qpos_init, qvel_init, T, sim::MJSimEnv, inverse_dynamics,
    num_iterations, X=nothing)
    force_init = zeros()
    if X == nothing
        force_init = zeros(size(sim.data.njmax))
        X = repeat(vcat(qpos_init, qvel_init), outer=[1, T]) + randn(size(X)[1], size(X)[2]) * .1
    end

    f = nothing
    f_x = nothing

    #Perform each iteration
    for iteration=1:num_iterations
        #collect inverse_dynamics inputs
        x_prev = vcat(X[:, [1]], X[:, 2:T]))
        x_curr = X
        x_next = vcat(X[:, 2:], X[:, [T]])

        #linearize inverse_dynamics
        f, f_x = dynamics_linearization(inverse_dynamics, [x_prev, x_curr, x_next], 1e-6)

        #calculate cost
        cost = []

        #get cost lineariation
        c, c_x, c_xx = create_cost_quadratic(cost)

        #update trajectory parameters
        dX = direct_solve(c, c_x, c_xx, .1)
        X += reshape(dX, size(X))
    end

    return f
end


text:
module LibGL

import DataStructures

using ModernGL

using CoreExtended
using FileManager

include("GLLists.jl")
include("GLDebugControl.jl")
include("GLExtendedFunctions.jl")
include("StorageManager.jl")
include("ShaderManager.jl")
include("GLTexture.jl")

using JLGEngine.GraphicsManager
using .GLLists
using .GLExtendedFunctions

function __init__()
end

lastMessage = (0, 0, 0, 0, 0)

#hasDuplicates(msg) = (for i=1:length(lastMessage) if msg[i] != lastMessage[i] return false end end; true)
hasDuplicates(msg) = (for (i,j) in zip(msg,lastMessage) if i != j return false end end; true)

function openglerrorcallback(
                source::GLenum, typ::GLenum,
                id::GLuint, severity::GLenum,
                len::GLsizei, message::Ptr{GLchar},
                userParam::Ptr{Void}
            )

		msg = (source, typ, id, severity, len)
		if hasDuplicates(msg) return end # ignore duplicates
		global lastMessage = msg

		#source = GL_DEBUG_SOURCE_API, GL_DEBUG_SOURCE_WINDOW_SYSTEM_, GL_DEBUG_SOURCE_SHADER_COMPILER, GL_DEBUG_SOURCE_THIRD_PARTY, GL_DEBUG_SOURCE_APPLICATION, GL_DEBUG_SOURCE_OTHER, GL_DONT_CARE
		#typ = GL_DEBUG_TYPE_ERROR, GL_DEBUG_TYPE_DEPRECATED_BEHAVIOR, GL_DEBUG_TYPE_UNDEFINED_BEHAVIOR, GL_DEBUG_TYPE_PORTABILITY, GL_DEBUG_TYPE_PERFORMANCE, GL_DEBUG_TYPE_MARKER, GL_DEBUG_TYPE_PUSH_GROUP, GL_DEBUG_TYPE_POP_GROUP, or GL_DEBUG_TYPE_OTHER, GL_DONT_CARE
		#severity = GL_DEBUG_SEVERITY_LOW, GL_DEBUG_SEVERITY_MEDIUM, or GL_DEBUG_SEVERITY_HIGH, GL_DONT_CARE
    errormessage = 	"\n"*
                    " __OPENGL________________________________________________________\n"*
                    "|\n"*
                    #"| type: $(GLENUM(typ).name) :: id: $(GLENUM(id).name)\n"*
										#"| source: $(GLENUM(source).name) :: severity: $(GLENUM(severity).name)\n"*
										(userParam != C_NULL ? "| UserParam : NOT NULL\n" : "")*
                    (len > 0 ? "| "*ascii(unsafe_string(message, len))*"\n" : "")*
										(haskey(LIST_ERROR,id) ? "| "*LIST_ERROR[id]*"\n" : "")*
                    "|________________________________________________________________\n"
    #if typ == GL_DEBUG_TYPE_ERROR
				warn(errormessage)
    #end
    nothing
end

const _openglerrorcallback = cfunction(openglerrorcallback, Void,
                                        (GLenum, GLenum,
                                        GLuint, GLenum,
                                        GLsizei, Ptr{GLchar},
                                        Ptr{Void}))

function SetDebugging(debug::Bool)
  if debug SetErrorCallBack() end
end

function SetErrorCallBack()
  @static if is_apple()
			warn("OpenGL debug message callback not available on osx")
			return
  end
	flags = getValByRef((ref) -> glGetIntegerv(GL_CONTEXT_FLAGS, ref), GLint(0))
  if (flags & GL_CONTEXT_FLAG_DEBUG_BIT) != 0
		glEnable(GL_DEBUG_OUTPUT)
		glEnable(GL_DEBUG_OUTPUT_SYNCHRONOUS)
		glDebugMessageCallbackARB(_openglerrorcallback, C_NULL)
	end
end

function init()
	SetErrorCallBack()
end

GetVersion(key::Symbol) = haskey(LIST_INFO,key) ? string("OpenGL ",glString(LIST_INFO[key])) : ""

const PREFIX = :(GL_)
const MOD = current_module()
function getMode(key::Symbol, value::Any)
	global MOD
	if key == :POLYGONMODE
		k=eval(Symbol(PREFIX,value[1]))
		v=eval(Symbol(PREFIX,value[2]))
		:($MOD.glPolygonMode($k,$v);)
	else
		:()
	end
end

#function createShader(source::FileSource, range=range(1,stat(source.path).size))
#	(dir,typ,ext) = splitFileName(source.path)
#	Shader(get(SHADER_FILE_TYPE, typ, :FRAGMENT),FileSourcePart(source,range))
#end

function resetBuffers()	StorageManager.unbindStorages() end

function resetRegisters()
	#StorageManager.unbindStorages()
	ShaderManager.reset()
end

#function linkTo(obj::AbstractGraphicsReference)
#	if obj == nothing return end
#	if isa(obj,AbstractGraphicsShaderProgram) ShaderManager.linkTo(obj) end
#end

end #LibGL

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
@testset "3 stages serial with $solver" for solver in lp_solvers
    include("prob5.2_data.jl")

    numScen = 2

    x = Matrix{Vector{JuMP.Variable}}(undef, 2,2)
    v = Matrix{Vector{JuMP.Variable}}(undef, 2,2)
    y = Matrix{Matrix{JuMP.Variable}}(undef, 2,2)
    models = Vector{JuMP.Model}(undef, 3)

    model = StructuredModel(num_scenarios=numScen)

    x[1,1] = @variable(model, [1:n], lowerbound=0)
    v[1,1] = @variable(model, [1:n], lowerbound=0)
    @constraints model begin
        x[1,1] .== v[1,1]
    end
    @objective(model, Min, dot(ic, v[1]))

    models[1] = model

    for s in 1:numScen
        model = StructuredModel(parent=models[1], prob=p2[s], same_children_as=(s == 1 ? nothing : models[2]), id=s, num_scenarios=numScen)
        y[1,s] = @variable(model, [1:n, 1:m], lowerbound=0)
        x[2,s] = @variable(model, [1:n], lowerbound=0)
        v[2,s] = @variable(model, [1:n], lowerbound=0)
        @constraints model begin
            x[2,s] .== x[1,1] + v[2,s]
            demand[j=1:m], sum(y[1,s][:,j]) == D2[j,s]
            ylim[i=1:n], sum(y[1,s][i,:]) <= x[1,1][i]
        end
        @objective(model, Min, dot(ic, v[2,s]) + dot(C, y[1,s] * T))
        if s == 1
            models[2] = model
        end
    end
    for s in 1:numScen
        model = StructuredModel(parent=models[2], prob=p2[s], same_children_as=(s == 1 ? nothing : models[3]), id=s)
        y[2,s] = @variable(model, [1:n, 1:m], lowerbound=0)
        @constraints model begin
            demand[j=1:m], sum(y[2,s][:,j]) == D2[j,s]
            ylim[i=1:n], sum(y[2,s][i,:]) <= x[2,1][i]
        end
        @objective(model, Min, dot(C, y[2,s] * T))
        if s == 1
            models[3] = model
        end
    end

    if !isdefined(@__MODULE__, :testniter_3stages)
        include("prob5.2_3stages_niter.jl")
    end
    fulltest(models[1], 3, 406712.49, [2986,0,7329,854], 402593.71614, 17319.095064, testniter_3stages, solver)
end


text:
include("cuextract.jl")

function investigate_function(io::IO, oh, dbgs, didata, name)
    sects = sections_for_oh(oh)[2]
    text = first(filter(x->contains(lowercase(bytestring(sectionname(x))),"text"),sects))
    for c in didata
        for sp in c.SPs
            if name != sp.linkage_name
                continue
            end
            rangelists = map(sp.variables) do loc
                if loc == 0
                    return []
                elseif loc == (-1 % UInt)
                    return UnitRange{UInt}[0:typemax(UInt)]
                else
                    return rangelist(dbgs,loc)
                end
            end
            println(io, "Dump for ", sp.name)
            range = sp.low_pc:sp.high_pc
            seek(oh,sectionoffset(text))
            data = readbytes(oh,sectionsize(text))
            disassemble2(data, range; io=io, rangelists = rangelists, rloffset = c.low_pc, offset = sectionoffset(text))
        end
    end
end

function investigate2(oh, debugoh=oh;
        dumpsource = false, dumpdies = false, dumploctables = false, fromproc = false)
    dorelocate = true

    # Relocate the object file
    if dorelocate
        new_buf = copy(oh.io)
        seekstart(new_buf)
        ELF.relocate!(new_buf,oh)
        seekstart(new_buf)
        oh = readmeta(new_buf)
    end

    dbgs = debugsections(debugoh);
    strtab = load_strtab(dbgs.debug_str)

    if dumpdies
        dt = DIETrees(dbgs)
        it = start(dt)
        while !done(dt, it)
            tree, it = next(dt, it)
            show(STDOUT, tree; strtab = strtab)
        end
    end

    didata = process_cus(dbgs, strtab)

#=
    llvmf = map(FindFunctionNamed,linkage_names)
    juliaf = map(x->x == C_NULL ? nothing : getLam(x),llvmf)
    if dumpsource
        for i = 1:length(llvmf)
            print(Base.uncompressed_ast(juliaf[i]).args[3])
            icxx" $(llvmf[i])->getParent()->dump(); "
        end
    end
    text = first(filter(x->contains(lowercase(bytestring(sectionname(x))),"text"),sects))
    rls = rangelists(dbgs)
    for i = 1:length(juliaf)
        println("Debug Info for ",names[i])
        if juliaf[i] == nothing
            println("AST not found")
        else
            ast = Base.uncompressed_ast(juliaf[i])
            println("Parameters:")
            params = ast.args[1]
            for var in params
                println(" - ",var.args[1])
            end
            lvars = ast.args[2][1]
            if !isempty(lvars)
                println("Local Variables:")
                for var in lvars
                    println(" - ",var[1])
                end
            end
        end
        range = ranges[i]
        @show range
        this_rls = collect(filter(x->!isempty(x),map(rls) do rl
            collect(filter(rl) do entry
                first(range) <= first(entry) && last(entry) <= last(range)
            end)
        end))
        if fromproc
            disassemble2(sectionaddress(text),sectionsize(text), range; rangelists=this_rls)
        else
            seek(oh,sectionoffset(text))
            data = readbytes(oh,sectionsize(text))
            disassemble2(data, range; rangelists=this_rls, offset = sectionoffset(text))
        end
    end
    if dumploctables
        seek(oh, sectionoffset(dbgs.debug_loc))
        while position(oh) < sectionoffset(dbgs.debug_loc)+sectionsize(dbgs.debug_loc)
            show(read(oh, DWARF.LocationList{UInt64}))
        end
    end
    (llvmf, juliaf)
=#
end


text:
"""
Interop types for CUDA native.
"""
immutable CUarray{T,N}
    ptr::Ptr{T}
    dims::N
    strides::N
end

CUDA.cubox(x::Union{CuArray,CuSubArray}) = CUDA.cubox(CUarray(x))
CUDA.cubox(t::NTuple{1,Int}) = CUDA.cubox(Cint1(t[1]))
CUDA.cubox(t::NTuple{2,Int}) = CUDA.cubox(Cint2(t[1],t[2]))
CUDA.cubox(t::NTuple{3,Int}) = CUDA.cubox(Cint3(t[1],t[2],t[3]))
CUDA.cubox(t::NTuple{4,Int}) = CUDA.cubox(Cint4(t[1],t[2],t[3],t[4]))
CUDA.cubox(t::NTuple{5,Int}) = CUDA.cubox(Cint5(t[1],t[2],t[3],t[4],t[5]))

CUarray(x::Union{CuArray,CuSubArray}) = CUarray(pointer(x), cint(size(x)), cint(strides(x)))

function CUarray3(x::CuArray, dim::Int)
    dims = Cint[1, size(x,dim), 1]
    for i = 1:dim-1
        dims[1] *= size(x, i)
    end
    for i = dim+1:ndims(x)
        dims[3] *= size(x, i)
    end
    cdims = Cint3(dims[1],dims[2],dims[3])
    cstrides = Cint3(Cint(1), cdims.i1, cdims.i1*cdims.i2)
    CUarray(pointer(x), cdims, cstrides)
end

immutable Cint1
    i1::Cint
end

immutable Cint2
    i1::Cint
    i2::Cint
end

immutable Cint3
    i1::Cint
    i2::Cint
    i3::Cint
end

immutable Cint4
    i1::Cint
    i2::Cint
    i3::Cint
    i4::Cint
end

immutable Cint5
    i1::Cint
    i2::Cint
    i3::Cint
    i4::Cint
    i5::Cint
end

cint(t::NTuple{1,Int}) = Cint1(t[1])
cint(t::NTuple{2,Int}) = Cint2(t[1],t[2])
cint(t::NTuple{3,Int}) = Cint3(t[1],t[2],t[3])
cint(t::NTuple{4,Int}) = Cint4(t[1],t[2],t[3],t[4])
cint(t::NTuple{5,Int}) = Cint5(t[1],t[2],t[3],t[4],t[5])
cint(i::Int...) = cint(i)


text:
using ApproxFunOrthogonalPolynomials, ApproxFunRational,
 ApproxFunFourier, ApproxFunBase, ApproxFun, AbstractIterativeSolvers, Plots

### Scalar case
L = 1.
g = x -> .2exp(-x^2)
G1 = Fun(zai(g),OscLaurent(1.,L),100) + Fun(zai(g),OscLaurent(0.,L),100)
h = x -> .2/(1 + x^2)
H2 = Fun(zai(h),OscLaurent(-1.,L),100)
𝓒 = Cauchy(-1)
𝓢 = Cauchy(1)

G=zeros(Fun,2,2)
G[1,1] = G1
G[2,2] = G1
G[1,2] = H2
G[2,1] = H2
G = copy(Fun(G))

One = pad(Fun(1.,OscLaurent(0.0,L)),40)
b = G*Fun([One;0*One])

op = x -> condense(x - G*(𝓒*x))
out = GMRES(op,b,⋅,1e-3,100, condense)
u = sum([out[2][i]*out[1][i] for i=1:length(out[2]) ])

id = Array([1. 0.; 0. 1.])
𝓢 = Cauchy(1)
𝓒 = Cauchy(-1)
c = id*[1.0; 0.]

id + G(.1)

jump = x -> (𝓢*u)(x) + c - (id + G(x))*((𝓒*u)(x)+c)
jump(.11)


text:
################################################################################

"""

    extractSignalBin(edfDf::DataFrame, params::Dict)

# Description
Use `extractSignalBin` on EDF file per channel from shell arguments. Returns a dictionary with channel names as keys.

See also: [`extractFFT`](@ref)
"""
function extractSignalBin(edfDf::DataFrame, params::Dict)
  if haskey(params, "window-size") && haskey(params, "bin-overlap")
    return extractSignalBin(edfDf, binSize = params["window-size"], binOverlap = params["bin-overlap"])
  else
    @error "Variables are not defined in dictionary"
  end
end

################################################################################

"""

    extractSignalBin(channel::Array;
    binSize::T, binOverlap::T) where T <: Number

# Description
Bin channel signal.

See also: [`extractFFT`](@ref)
"""
function extractSignalBin(channel::Array; binSize::T, binOverlap::T) where T <: Number
  # define variables
  stepSize = floor(Int64, binSize / binOverlap)
  signalSteps = 1:stepSize:length(channel)
  signalAr = Array{Float64}(
    undef,
    length(signalSteps),
    binSize
  )

  # iterate over signal bins
  for stepIx ∈ 1:length(signalSteps)
    signalBoundry = signalSteps[stepIx] + binSize - 1

    # extract signals
    if signalBoundry <= length(channel)
      signalAr[stepIx, :] = channel[signalSteps[stepIx]:signalBoundry]

    # adjust last bin
    elseif signalBoundry > length(channel)
      signalAr[stepIx, :] = [
        channel[signalSteps[stepIx]:end];
        (signalBoundry - length(channel) |> abs |> zeros)
      ]
    end

  end
  return signalAr
end

################################################################################

"""

    extractSignalBin(edfDf::DataFrames.DataFrame;
    binSize::T, binOverlap::T) where T <: Number

# Description
Use `extractSignalBin` on EDF file per channel. Returns a dictionary with channel names as keys.

See also: [`extractFFT`](@ref)
"""
function extractSignalBin(edfDf::DataFrames.DataFrame; binSize::T, binOverlap::T) where T <: Number
  @info("Binning channels signals...")
  channelDc = Dict()
  signalAr = begin
    stepSize = floor(Int64, binSize / binOverlap)
    signalSteps = 1:stepSize:size(edfDf, 1)
    Array{Float64}(
      undef,
      1,
      binSize,
      length(signalSteps)
    )
  end

  # iterate on dataframe channels
  for (ψ, ε) ∈ enumerate(names(edfDf))
    tmpAr = extractSignalBin(edfDf[:, ψ], binSize = binSize, binOverlap = binOverlap)
    for β ∈ 1:size(tmpAr, 1)
      signalAr[:, :, β] = tmpAr[β, :]
    end
    channelDc[ε] = copy(signalAr)
  end
  return channelDc
end

################################################################################

max_stars_repo_name: DanielRivasMD/MindReader

max_issues_repo_name: DanielRivasMD/MindReader

max_forks_repo_name: DanielRivasMD/MindReader


text:
"""
    diffci(x1, n1, x2, n2; level = 0.95, method = :default)

Proportion difference (x1 / n1 - x2 / n2) confidence interval.

* 'method'

- `:mn` | `:default`- Miettinen & Nurminen;  Miettinen, O. and Nurminen, M. (1985), Comparative analysis of two rates. Statist. Med., 4: 213-226. doi:10.1002/sim.4780040211;
- `:wald` - Wald CI without CC;
- `:waldcc` - Wald CI with CC;
- `:nhs`  - Newcombes Hybrid (wilson) Score interval; Newcombe RG (1998), Interval Estimation for the Difference Between Independent Proportions: Comparison of Eleven Methods. Statistics in Medicine 17, 873-890;
- `:nhscc` - Newcombes Hybrid Score CC; Newcombe (1998);
- `:ac` -  Agresti-Caffo interval; Agresti A, Caffo B., “Simple and effective confidence intervals for proportions and differences of proportions result from adding two successes and two failures”, American Statistician 54: 280–288 (2000);
- `:ha` - Hauck-Andersen; Hauck, W. W., & Anderson, S. (1986). A Comparison of Large-Sample Confidence Interval Methods for the Difference of Two Binomial Probabilities. The American Statistician, 40(4), 318–322. doi:10.1080/00031305.1986.10475426 ;
- `:mnmee` - Mee maximum likelihood method; Mee RW (1984) Confidence bounds for the difference between two probabilities,Biometrics40:1175-1176
- `:mover` - Method of variance estimates recovery;
- `:jeffrey` - Brown, Li's Jeffreys.

# Reference:

* Brown, L.D., Cai, T.T., and DasGupta, A. Interval estimation for a binomial proportion. Statistical Science, 16(2):101–117, 2001.
* Farrington, C. P. and Manning, G. (1990), “Test Statistics and Sample Size Formulae for Comparative Binomial Trials with Null Hypothesis of Non-zero Risk Difference or Non-unity Relative Risk,” Statistics in Medicine, 9, 1447–1454
* Li HQ, Tang ML, Wong WK. Confidence intervals for ratio of two Poisson rates using the methodof variance estimates recovery. Computational Statistics 2014; 29(3-4):869-889
* Brown, L., Cai, T., & DasGupta, A. (2003). INTERVAL ESTIMATION IN EXPONENTIAL FAMILIES. Statistica Sinica, 13(1), 19-49.
"""
function diffci(x1, n1, x2, n2; level = 0.95, method = :default)

    alpha    = 1 - level
    if method == :mn || method == :default
        ci_diff_mn(x1, n1, x2, n2, alpha)
    elseif method == :wald
        ci_diff_wald(x1, n1, x2, n2, alpha)
    elseif method == :waldcc
        ci_diff_wald_cc(x1, n1, x2, n2, alpha)
    elseif method == :nhs
        ci_diff_nhs(x1, n1, x2, n2, alpha)
    elseif method == :nhscc
        ci_diff_nhs_cc(x1, n1, x2, n2, alpha)
    elseif method == :ac
        ci_diff_ac(x1, n1, x2, n2, alpha)
    elseif method == :ha
        ci_diff_ha(x1, n1, x2, n2, alpha)
    elseif method == :mnmee
        ci_diff_mnmee(x1, n1, x2, n2, alpha)
    elseif method == :mover
        ci_diff_mover(x1, n1, x2, n2, alpha)
    elseif method == :jeffrey
        ci_diff_jeffrey(x1, n1, x2, n2, alpha)
    else
        throw(ArgumentError("unknown ci method=$(method)"))
    end
end
"""
    diffci(contab::ConTab; level = 0.95, method = :default)

ConTab 2X2:

A | B
--|--
C | D

Difference: A / (A + B) - C / (C + D)
"""
function diffci(contab::ConTab; level = 0.95, method = :default)
    if !(size(contab.tab, 1) == size(contab.tab, 2) == 2) throw(ArgumentError("CI only for 2 X 2 tables.")) end
    alpha    = 1 - level
    x1 = contab.tab[1,1]
    n1 = x1 + contab.tab[1,2]
    x2 = contab.tab[2,1]
    n2 = x2 + contab.tab[2,2]
    diffci(x1, n1, x2, n2; level = level, method = method)
end
"""
    orci(contab::ConTab; level = 0.95, method = :default)

- `:mn`
- `:woolf`
- `:awoolf`
- `:mover`
"""
function orci(contab::ConTab; level = 0.95, method = :default)
    if !(size(contab.tab, 1) == size(contab.tab, 2) == 2) throw(ArgumentError("CI only for 2 X 2 tables.")) end
    alpha    = 1 - level
    x1 = contab.tab[1,1]
    n1 = x1 + contab.tab[1,2]
    x2 = contab.tab[2,1]
    n2 = x2 + contab.tab[2,2]
    if method == :mn || method == :default
        ci_or_mn(x1, n1, x2, n2, alpha)
    elseif method == :woolf
        ci_or_woolf(x1, n1, x2, n2, alpha)
    elseif method == :awoolf
        ci_or_awoolf(x1, n1, x2, n2, alpha)
    elseif method == :mover
        ci_or_mover(x1, n1, x2, n2, alpha)
    else
        throw(ArgumentError("unknown ci method=$(method)"))
    end
end
"""
    rrci(contab::ConTab; level = 0.95, method = :default)

- `:mn`
- `:cli`
- `:li`
- `:mover`
"""
function rrci(contab::ConTab; level = 0.95, method = :default)
    if !(size(contab.tab, 1) == size(contab.tab, 2) == 2) throw(ArgumentError("CI only for 2 X 2 tables.")) end
    alpha    = 1 - level
    x1 = contab.tab[1,1]
    n1 = x1 + contab.tab[1,2]
    x2 = contab.tab[2,1]
    n2 = x2 + contab.tab[2,2]
    if method == :mn || method == :default
        ci_rr_mn(x1, n1, x2, n2, alpha)
    elseif method == :cli
        ci_rr_cli(x1, n1, x2, n2, alpha)
    elseif method == :li
        ci_rr_li(x1, n1, x2, n2, alpha)
    elseif method == :mover
        ci_rr_mover(x1, n1, x2, n2, alpha)
    else
        throw(ArgumentError("unknown ci method=$(method)"))
    end
end

"""
    propci(x::Int, n::Int; level = 0.95, method = :default)

`method`:

- `:wilson` | `:default` - Wilson's confidence interval (CI) for a single proportion (wilson score);
- `:wilsoncc` - Wilson's CI with continuity correction (CC);
- `:cp` - Clopper-Pearson exact CI;
- `:blaker` - Blaker exact CI for discrete distributions;
- `:soc` - SOC: Second-Order corrected CI;
- `:arc` - Arcsine CI;
- `:wald` - Wald CI without CC;
- `:waldcc` - Wald CI with CC;
- `:ac`
- `:jeffrey`
"""
function propci(x::Int, n::Int; level = 0.95, method = :default)
    if  x > n throw(ArgumentError("x > n")) end
    alpha    = 1 - level
    if method == :wilson || method == :default
        fx = ci_prop_wilson
    elseif method==:wilsoncc
        fx = ci_prop_wilson_cc
    elseif method==:cp
        fx = ci_prop_cp
    elseif method==:blaker
        fx = ci_prop_blaker
    elseif method==:soc
        fx = ci_prop_soc
    elseif method==:arc
        fx = ci_prop_arc
    elseif method==:wald
        fx = ci_prop_wald
    elseif method==:waldcc
        fx = ci_prop_wald_cc
    elseif method==:ac
        fx = ci_prop_ac
    elseif method==:jeffrey
        fx = ci_prop_jeffrey
    else
        throw(ArgumentError("unknown method!"))
    end
    fx(x, n, alpha)
end
"""
    propci(contab::ConTab; level = 0.95, method = :default)
"""
function propci(contab::ConTab; level = 0.95, method = :default)
    alpha    = 1 - level
    if  size(contab.tab, 2) != 2 throw(ArgumentError("CI only for N X 2 tables.")) end
    if size(contab, 1) > 1
        v = Vector{Tuple{Float64, Float64}}(undef, size(contab, 1))
        for i = 1:size(contab, 1)
            x = contab.tab[i,1]
            n = x + contab.tab[i,2]
            println(x, " : ", n)
            v[i] = propci(x, n; level = level, method = method)
        end
        return Tuple(v)
    else
        x = contab.tab[1,1]
        n = x + contab.tab[1,2]
        return propci(x, n; level = level, method = method)
    end
end

"""
    propci(prop::Proportion; level = 0.95, method = :default)
"""
function propci(prop::Proportion; level = 0.95, method = :default)
    propci(prop.x, prop.n; level = level, method = method)
end

"""
    mpropci(contab::ConTab; level = 0.95, method = :default)

Multinomial proportions confidence interval.

`method`:

- `goodman` | `default` Goodman, L.A. (1965). On Simultaneous Confidence Intervals for Multinomial Proportions. Technometrics 7: 247-254.

"""
function mpropci(contab::ConTab; level = 0.95, method = :default)
    alpha    = 1 - level
    if  size(contab.tab, 2) <= 2 throw(ArgumentError("CI only for N X M tables where M > 2")) end

    if method == :goodman || method == :default
        fx = ci_prop_goodman
    else
        throw(ArgumentError("unknown method!"))
    end
    if size(contab, 1) > 1
        v = Vector{Vector{Tuple{Float64, Float64}}}(undef, size(contab, 1))
        for i = 1:size(contab, 1)
            v[i] = fx(view(contab.tab, i, :), alpha)
        end
        return v
    else
        return fx(view(contab.tab, 1, :), alpha)
    end
end

################################################################################
# Proportion difference CI
################################################################################
# Method of Mee 1984 with Miettinen and Nurminen modification n / (n - 1) Newcombe 1998
# Score intervals for the difference of two binomial proportions
@inline function mle_diff(p1, n1, p2, n2, δ)
    if p1 - p2 - δ == 0 return 0.0 end
    θ = n2 / n1
    a = 1 + θ
    b = -(1 + θ + p1 + θ * p2 + δ * (θ + 2))
    c = δ^2 + δ * (2 * p1 + θ + 1) + p1 + θ * p2
    d = -p1 * δ * (1 + δ)
    v = (b / a / 3)^3 - b * c/(6 * a * a) + d / 2 / a
    u = sign(v) * sqrt((b / 3 / a)^2 - c / 3 / a)
    w = (pi + acos(v / u^3)) / 3
    p1n = 2 * u * cos(w) - b / 3 /a
    p2n = p1n - δ
    return p1n, p2n
end
@inline function mn_diff_z_val(p1, n1, p2, n2, est, δ)
    p1n, p2n = mle_diff(p1, n1, p2, n2, δ)
    return (est - δ)^2 / ((n1 + n2) / (n1 + n2 - 1) * (p1n * (1 - p1n) / n1 + p2n * (1 - p2n) / n2))
end
function ci_diff_mn(x1, n1, x2, n2, alpha; atol::Float64 = 1E-8)
    lcis, ucis = ci_diff_nhs_cc(x1, n1, x2, n2, alpha)
    p1       = x1 / n1
    p2       = x2 / n2
    est      = p1 - p2
    z        = quantile(Chisq(1), 1 - alpha)
    fmnd(x)  = mn_diff_z_val(p1, n1, p2, n2, est, x) - z
    if fmnd(lcis) * fmnd(est - eps()) < 0.0
        ll = lcis
        lu = est - eps()
    else
        ll = -1.0 + eps()
        lu = lcis
    end
    if fmnd(ucis) * fmnd(est + eps()) < 0.0
        ul = est + eps()
        uu = ucis
    else
        ul = ucis
        uu = 1.0 - eps()
    end
    lci = find_zero(fmnd, (ll, lu), atol = atol)
    uci = find_zero(fmnd, (ul, uu), atol = atol)
    return lci, uci
end
# Wald
function ci_diff_wald(x1, n1, x2, n2, alpha)
    p1       = x1 / n1
    p2       = x2 / n2
    est      = p1 - p2
    z        = quantile(Normal(), 1 - alpha / 2)
    se       = sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)
    return est - z * se, est + z * se
end
# Wald continuity correction
function ci_diff_wald_cc(x1, n1, x2, n2, alpha)
    p1       = x1 / n1
    p2       = x2 / n2
    est      = p1 - p2
    z        = quantile(Normal(), 1 - alpha / 2)
    se       = sqrt(p1 * (1 - p1) / n1 + p2 * (1 - p2) / n2)
    cc       = (1 / n1 + 1 / n2) / 2
    return est - z * se - cc, est + z * se + cc
end
# Newcombes Hybrid (wilson) Score interval for the difference of proportions
# Newcombe 1998
function ci_diff_nhs(x1, n1, x2, n2, alpha)
    p1       = x1 / n1
    p2       = x2 / n2
    est      = p1 - p2
    z        = quantile(Normal(), 1 - alpha / 2)
    lci1, uci1 = ci_prop_wilson(x1, n1, alpha)
    lci2, uci2 = ci_prop_wilson(x2, n2, alpha)
    return  est - z * sqrt(lci1 * (1 - lci1)/n1 + uci2 * (1 - uci2) / n2), est + z * sqrt(uci1 * (1 - uci1) / n1 + lci2 * (1 - lci2) / n2)
end
# Newcombes Hybrid Score continuity correction
function ci_diff_nhs_cc(x1, n1, x2, n2, alpha)
    p1       = x1 / n1
    p2       = x2 / n2
    est      = p1 - p2
    z        = quantile(Normal(), 1 - alpha / 2)
    lci1, uci1 = ci_prop_wilson_cc(x1, n1, alpha)
    lci2, uci2 = ci_prop_wilson_cc(x2, n2, alpha)
    return  est - sqrt((p1 - lci1)^2 + (uci2 - p2)^2), est + sqrt((uci1 - p1)^2 + (p2 - lci2)^2)
end
# Agresti-Caffo interval for the difference of proportions
# Agresti A, Caffo B., “Simple and effective confidence intervals for proportions and differences of proportions result from adding two successes and two failures”, American Statistician 54: 280–288 (2000)
function ci_diff_ac(x1, n1, x2, n2, alpha)
    z        = quantile(Normal(), 1 - alpha / 2)
    p1I      = (x1 + 1) / (n1 + 2)
    p2I      = (x2 + 1) / (n2 + 2)
    n1I      = n1 + 2
    n2I      = n2 + 2
    est      = p1I - p2I
    se       = sqrt(p1I * (1 - p1I) / n1I + p2I * (1 - p2I) / n2I)
    return  est - z * se, est + z * se
end
# Hauck-Andersen
# Hauck, W. W., & Anderson, S. (1986). A Comparison of Large-Sample Confidence Interval Methods for the Difference of Two Binomial Probabilities. The American Statistician, 40(4), 318–322. doi:10.1080/00031305.1986.10475426

function ci_diff_ha(x1, n1, x2, n2, alpha)
    p1       = x1 / n1
    p2       = x2 / n2
    est      = p1 - p2
    z        = quantile(Normal(), 1 - alpha / 2)
    se       = sqrt(p1 * (1 - p1) / (n1 - 1) + p2 * (1 - p2) / (n2 - 1))
    cc       = 1 / min(n1, n2)
    return est - z * se - cc, est + z * se + cc
end
# Mee RW (1984) Confidence bounds for the difference between two probabilities,Biometrics40:1175-1176
# MN - no correction
@inline function mnmee_z_val(p1, n1, p2, n2, est, δ)
    p1n, p2n = mle_diff(p1, n1, p2, n2, δ)
    return (est - δ)^2 / (p1n * (1 - p1n) / n1 + p2n * (1 - p2n) / n2)
end
function ci_diff_mnmee(x1, n1, x2, n2, alpha; atol::Float64 = 1E-8)
    lcis, ucis = ci_diff_nhs_cc(x1, n1, x2, n2, alpha)
    p1       = x1 / n1
    p2       = x2 / n2
    est      = p1 - p2
    z        = quantile(Chisq(1), 1 - alpha)
    fmnd(x)  = mnmee_z_val(p1, n1, p2, n2, est, x) - z
    if fmnd(lcis) * fmnd(est - eps()) < 0.0
        ll = lcis
        lu = est - eps()
    else
        ll = -1.0 + eps()
        lu = lcis
    end
    if fmnd(ucis) * fmnd(est + eps()) < 0.0
        ul = est + eps()
        uu = ucis
    else
        ul = ucis
        uu = 1.0 - eps()
    end
    lci = find_zero(fmnd, (ll, lu), atol = atol)
    uci = find_zero(fmnd, (ul, uu), atol = atol)
    return  lci, uci
end
# Brown, Li's Jeffreys
function ci_diff_jeffreys(x1, n1, x2, n2, alpha)
    p1   = (x1 +  1/2) / (n1 + 1)
    p2   = (x2 +  1/2) / (n2 + 1)
    se   = sqrt(p1*(1 - p1) / n1 + p2 * (1 - p2) / n2)
    z    = quantile(Normal(), 1 - alpha / 2)
    est  = p1 - p2
    return  max(-1.0, est - z * se), min(1.0, est + z * se)
end
# Method of variance estimates recovery
function ci_diff_mover(x1, n1, x2, n2, alpha)
    p1       = x1 / n1
    p2       = x2 / n2
    est      = p1 - p2
    Z        = quantile(Normal(), 1 - alpha / 2)
    lci1, uci1   = ci_prop_wilson(x1, n1, alpha)
    lci2, uci2   = ci_prop_wilson(x2, n2, alpha)
    lci      = est - sqrt((p1 - lci1)^2 + (uci2 - p2)^2)
    uci      = est + sqrt((uci1 - p1)^2 + (p2 - lci2)^2)
    return lci, uci
end
################################################################################
# Odd ratio CI
################################################################################
@inline function mle_or(φ, x1, n1, x2, n2)
    a  = n2 * (φ-1)
    b  = φ * n1 + n2 - (x1 + x2) * (φ - 1)
    c  = -(x1 + x2)
    p2 = (-b + sqrt(b * b - 4 * a * c)) / a / 2
    p1 = p2 * φ/(1 + p2 * (φ - 1))
    return p1, p2
end
@inline function mle_or_z_val(φ, x1, n1, x2, n2)
    p1 = x1 / n1
    pmle1, pmle2 = mle_or(φ, x1, n1, x2, n2)
    return (n1 * (p1 - pmle1))^2 * (1 / (n1 * pmle1 * (1 - pmle1)) + 1/(n2 * pmle2 * (1 - pmle2))) / ((n1 + n2)/(n1 + n2 - 1))
end
################################################################################
# Miettinen O. S., Nurminen M. (1985) Comparative analysis of two rates.Statistics in Medicine4,213–226
# MN Score
function ci_or_mn(x1, n1, x2, n2, alpha; atol::Float64 = 1E-8)
    z        = quantile(Chisq(1), 1 - alpha)
    fmnor(x) = mle_or_z_val(x, x1, n1, x2, n2) - z
    if (x1== 0 && x2 == 0) || (x1 == n1 && x2 == n2)
        return 0.0, Inf
    elseif x1==0 || x2 == n2
        return 0.0, find_zero(fmnor, atol, atol = atol)
    elseif x1 == n1 || x2 == 0
        return find_zero(fmnor, atol, atol = atol), Inf
    else
        est  = (x1 / (n1 - x1)) / (x2 / (n2 - x2))
        return find_zero(fmnor, atol, atol = atol), find_zero(fmnor, est + atol, atol = atol)
    end
end
# Woolf logit
# Woolf, B. (1955). On estimating the relation between blood group and disease. Annals of human genetics, 19(4):251-253.
function ci_or_woolf(x1, n1, x2, n2, alpha)
    xa        = x1
    xb        = n1 - x1
    xc        = x2
    xd        = n2 - x2
    est       = xa*xd/xc/xb
    estI      = log(est)
    se        = sqrt(1/xa + 1/xb + 1/xc + 1/xd)
    z         = quantile(Normal(), 1 - alpha / 2)
    return exp(estI - z * se), exp(estI + z * se)
end
# Adjusted Woolf interval (Gart adjusted logit) Lawson, R (2005):Smallsample confidence intervals for the odds ratio.  Communication in Statistics Simulation andComputation, 33, 1095-1113.
# Gart, J. J. (1966). Alternative analyses of contingency tables. Journal of the Royal Statistical Society. Series B (Methodological), 28:164-179.
function ci_or_awoolf(x1, n1, x2, n2, alpha)
    xa        = x1 +  1/2
    xb        = n1 - x1 +  1/2
    xc        = x2 + 1/2
    xd        = n2 - x2 +  1/2
    est       = xa*xd/xc/xb
    estI      = log(est)
    se        = sqrt(1/xa + 1/xb + 1/xc + 1/xd)
    z         = quantile(Normal(), 1 - alpha/2)
    return exp(estI - z * se), exp(estI + z * se)
end
# Method of variance estimates recovery
# Donner, A. and Zou, G. (2012). Closed-form confidence intervals for functions of the normal mean and standard deviation. Statistical Methods in Medical Research, 21(4):347-359.
function ci_or_mover(x1, n1, x2, n2, alpha)
    p1       = (x1/(n1-x1))
    p2       = (x2/(n2-x2))
    est      = p1/p2
    z        = quantile(Normal(), 1 - alpha / 2)
    lci1, uci1   = ci_prop_wilson(x1, n1, alpha)
    lci2, uci2   = ci_prop_wilson(x2, n2, alpha)
    vl1      = lci1/(1 - lci1)
    vu1      = uci1/(1 - uci1)
    vl2      = lci2/(1 - lci2)
    vu2      = uci2/(1 - uci2)
    lci      = (p1*p2-sqrt((p1*p2)^2 - vl1*vu2*(2*p1-vl1)*(2*p2-vu2)))/(vu2*(2*p2 - vu2))
    uci      = (p1*p2+sqrt((p1*p2)^2 - vu1*vl2*(2*p1-vu1)*(2*p2-vl2)))/(vl2*(2*p2 - vl2))
    return lci, uci
end
################################################################################
# Risk ratio CI
################################################################################
@inline function mle_rr(φ, x1, n1, x2, n2)
    a = (n1 + n2) * φ
    b = -(φ * (x1 + n2) + x2 + n1)
    c = x1 + x2
    p2 = (-b - sqrt(b * b - 4 * a * c)) / 2 / a
    p1 = p2 * φ
    return p1, p2
end
@inline function mle_rr_z_val(φ, x1, n1, x2, n2)
    p1 = x1 / n1
    p2 = x2 / n2
    pmle1, pmle2 = mle_rr(φ, x1, n1, x2, n2)
    return (p1 - φ * p2)^2 / ((pmle1 * (1 - pmle1) / n1 + φ * φ * pmle2 * (1 - pmle2) / n2) * ((n1 + n2 - 1) / (n1 + n2)))
end
################################################################################
# Miettinen-Nurminen Score interval
# Miettinen, O. and Nurminen, M. (1985), Comparative analysis of two rates. Statist. Med., 4: 213-226. doi:10.1002/sim.4780040211
function ci_rr_mn(x1, n1, x2, n2, alpha; atol::Float64 = 1E-8)
    z        = quantile(Chisq(1), 1 - alpha)
    fmnrr(x) = mle_rr_z_val(x, x1, n1, x2, n2) - z
    if (x1==0 && x2==0) || (x1==n1 && x2==n2)
        return  0.0, Inf
    elseif x1==0 || x2==n2
        return 0.0, find_zero(fmnrr, atol, atol=atol)
    elseif x1==n1 || x2 == 0
        return find_zero(fmnrr, atol, atol=atol), Inf
    else
        est = (x1 / n1) / (x2 / n2)
        return find_zero(fmnrr, atol, atol=atol), find_zero(fmnrr, est + atol, atol=atol)
    end
end
# Crude log interval
# Gart, JJand Nam, J (1988): Approximate interval estimation of the ratio of binomial parameters: Areview and corrections for skewness. Biometrics 44, 323-338.
function ci_rr_cli(x1, n1, x2, n2, alpha)
    x1I       = x1 + 1/2
    x2I       = x2 + 1/2
    n1I       = n1 + 1/2
    n2I       = n2 + 1/2
    estI      = log((x1I / n1I) / (x2I / n2I))
    se        = sqrt(1 / x2I + 1 / x1I - 1 / n2I - 1 / n1I)
    est       = (x1 / n1) / (x2 / n2)
    z         =  quantile(Normal(), 1 - alpha / 2)
    return exp(estI - z * se), exp(estI + z * se)
end
# Katz D, Baptista J, Azen SP and Pike MC. Obtaining confidence intervals for the risk ratio in cohort studies. Biometrics 1978; 34: 469–474
function ci_rr_li(x1, n1, x2, n2, alpha)
    est       = (x1 / n1) / (x2 / n2)
    estI      = log(est)
    se        = sqrt(1 / x2 + 1 / x1 - 1 / n2 - 1 / n1)
    z         = quantile(Normal(), 1 - alpha / 2)
    return exp(estI - z * se), exp(estI + z * se)
end
# Method of variance estimates recovery (Donner, Zou, 2012)
function ci_rr_mover(x1, n1, x2, n2, alpha)
    p1       = x1 / n1
    p2       = x2 / n2
    est      = p1 / p2
    Z        = quantile(Normal(), 1 - alpha / 2)
    lci1, uci1   = ci_prop_wilson(x1, n1, alpha)
    lci2, uci2   = ci_prop_wilson(x2, n2, alpha)
    lci      = (p1 * p2 - sqrt((p1 * p2)^2 - lci1 * uci2 * (2 * p1 - lci1) * (2 * p2 - uci2))) / (uci2 * (2 * p2 - uci2))
    uci      = (p1 * p2 + sqrt((p1 * p2)^2 - uci1 * lci2 * (2 * p1 - uci1) * (2 * p2 - lci2))) / (lci2 * (2 * p2 - lci2))
    return  lci, uci
end
################################################################################
# Prorotion CI
################################################################################
# Wilson’s confidence interval for a single proportion, wilson score
# Wilson, E.B. (1927) Probable inference, the law of succession, and statistical inferenceJ. Amer.Stat. Assoc22, 209–212
function ci_prop_wilson(x, n, alpha)
    z   = abs(quantile(Normal(), 1 - alpha / 2))
    p   = x / n
    d   = 1 + (z^2) / n
    se  = sqrt((p * (1 - p) + (z^2) / (4 * n)) / n) / d
    est = (p + (z^2) / (2 * n)) / d
    return est - z * se, est + z * se
end
# Wilson CC
# Newcombe, R. G. (1998). "Two-sided confidence intervals for the single proportion: comparison of seven methods". Statistics in Medicine. 17 (8): 857–872. doi:10.1002/(SICI)1097-0258(19980430)17:8<857::AID-SIM777>3.0.CO;2-E. PMID 959561
function ci_prop_wilson_cc(x, n, alpha)
    z = abs(quantile(Normal(), 1 - alpha / 2))
    p = x / n
    l = (2*n*p+z*z-1-z*sqrt(z*z-2-1/n+4*p*(n*(1-p)+1)))/2/(n+z*z)
    u = (2*n*p+z*z+1+z*sqrt(z*z+2-1/n+4*p*(n*(1-p)-1)))/2/(n+z*z)
    return min(p, l), max(p, u)
end
#Clopper-Pearson exatct CI
#Clopper, C. and Pearson, E.S. (1934) The use of confidence or fiducial limits illustrated in the caseof the binomial.Biometrika26, 404–413.
function ci_prop_cp(x, n, alpha)
    if x == 0
        ll = 0.0
        ul = 1 - (alpha / 2)^(1 / n)
    elseif x == n
        ul = 1.0
        ll = (alpha / 2)^(1 / n)
    else
        ll = 1/(1+(n-x+1)/(x*quantile(FDist(2*x, 2*(n-x+1)), alpha/2)))
        ul = 1/(1+(n-x) / ((x+1)*quantile(FDist(2*(x+1), 2*(n-x)), 1-alpha/2)))
    end
    return ll, ul
end
# Blaker CI
# Blaker, H. (2000). Confidence curves and improved exact confidence intervals for discrete distributions,Canadian Journal of Statistics28 (4), 783–798
function ci_prop_blaker(x, n, alpha; atol::Float64 = 1E-8)
    lower = 0; upper = 1;
    fx(p) =  acceptbin(x, n, p) - alpha
    if n != 0
        lower = quantile(Beta(x, n-x+1), alpha/2)
        lower = find_zero(fx, lower, atol=atol)
    end
    if x != n
        upper = quantile(Beta(x+1, n-x), 1-alpha/2)
        upper = find_zero(fx, upper, atol=atol)
    end
    return lower,upper
end
@inline function acceptbin(x, n, p)
    BIN = Binomial(n,p)
    p1 = 1-cdf(BIN,x-1)
    p2 =   cdf(BIN,x)
    a1 = p1 + cdf(BIN, quantile(BIN,p1)-1)
    a2 = p2+1-cdf(BIN, quantile(BIN,1-p2))
    return min(a1,a2)
end
# SOC  Second-Order corrected
# T. Tony Cai One-sided confdence intervals in discrete distributions doi:10.1016/j.jspi.2004.01.00
function ci_prop_soc(x, n, alpha)
    p  = x/n
    k  = quantile(Normal(), 1-alpha/2)
    k2 = k^2
    η  = k2/3+1/6
    γ1 = -(k2*13/18+17/18)
    γ2 = k2/18+7/36
    m  = (x+η)/(n+2*η)
    b  = k*sqrt(p*(1-p)+(γ1*p*(1-p)+γ2)/n)/sqrt(n)
    return m-b, m+b
end
# Arcsine
function ci_prop_arc(x, n, alpha)
    q = quantile(Normal(), 1-alpha/2)
    p = x/n
    z = q/(2*sqrt(n))
    return sin(asin(sqrt(p))-z)^2, sin(asin(sqrt(p))+z)^2
end
# Wald CI
function ci_prop_wald(x, n, alpha)
    p=x/n
    b = quantile(Normal(), 1-alpha/2)*sqrt(p*(1-p)/n)
    return p-b, p+b
end
# Wald CI CC
function ci_prop_wald_cc(x::Int, n::Int, alpha::Real)
    p=x/n
    b = quantile(Normal(), 1-alpha/2)*sqrt(p*(1-p)/n)
    cc = 0.5 / n
    return p-b-cc, p+b+cc
end
# Agresti-Coull
function ci_prop_ac(x, n, alpha)
    z = quantile(Normal(), 1 - alpha / 2)
    n = n + z^2
    est = (x + z ^ 2 / 2) / n
    se = sqrt(est * (1 - est) / n)
    return (est-z*se, est+z*se)
end

# Jeffreys interval
function ci_prop_jeffrey(x, n, alpha)
    (quantile(Beta(x + 1/2, n - x + 1/2), alpha/2), quantile(Beta(x + 1/2, n - x + 1/2), 1-alpha/2))
end

################################################################################
# Goodman CI
################################################################################
# Goodman, L.A. (1965). On Simultaneous Confidence Intervals for Multinomial Proportions. Technometrics 7: 247-254.
# https://blogs.sas.com/content/iml/2017/02/15/confidence-intervals-multinomial-proportions.html
# https://rdrr.io/cran/CoinMinD/man/GM.html
function ci_prop_goodman(v, alpha::T2) where T2
    k   = length(v)
    s   = sum(v)
    p   = v ./ s
    chi = quantile(Chisq(one(Int)), one(T2) - alpha/k)
    ci  = Vector{Tuple{Float64, Float64}}(undef, k)
    d   = 2 * (chi + s)
     @inbounds @simd for i = 1:k
        ci[i] = ( (chi + 2v[i] - sqrt(chi*chi + 4v[i]*chi*(1.0 - v[i]/s))) / d , (chi + 2v[i] + sqrt(chi*chi + 4v[i]*chi*(1.0 - v[i]/s))) / d )
    end
    ci
end