text:
# This file contains some popular informats
# This is experimental currently

Base.@propagate_inbounds function strip_location(x, lo, hi)
    for i in lo:hi
        if x.data[i] != 0x20
            lo = i
            break
        end
    end
    for i in hi:-1:lo
        if x.data[i] != 0x20
            hi = i
            break
        end
    end
    lo, hi
end

# In general any function defined as informat must have these specifications:
# * it must take three positional arguments, x, lo, hi, where x is a custom structure and x.data is a vector of UInt8
# * function can change the values of x.data but only within lo:hi
# * function must do the operations in place and return nothing

function COMMA!(x, lo, hi)
    cnt = hi
    for i in hi:-1:lo
        if !(x.data[i] in (UInt8(','), UInt8('$'), UInt8('£')))
            x.data[cnt] = x.data[i]
            cnt -= 1
            cnt < lo && break
        end
    end
    # for i in lo:cnt
    #     x.data[i] = 0x20
    # end
    fill!(view(x.data, lo:cnt), 0x20)
    nothing
end
function COMMAX!(x, lo, hi)
    cnt = hi
    for i in hi:-1:lo
        if !(x.data[i] in (UInt8('.'))) #  TODO we should take care of('€') it is not UInt8
            if x.data[i] == UInt8(',')
                x.data[cnt] = UInt8('.')
            else
                x.data[cnt] = x.data[i]
            end
            cnt -= 1
            cnt < lo && break
        end
    end
    # for i in lo:cnt
    #     x.data[i] = 0x20
    # end
    fill!(view(x.data, lo:cnt), 0x20)
    nothing
end

# treats NA,na,Na,nA,... as missing value for numeric columns
function NUM_NA!(x, lo, hi)
    lo, hi = strip_location(x, lo, hi)
    flag = false
    if hi-lo+1 == 2
        if x.data[lo] in (UInt8('N'), UInt8('n')) && x.data[hi] in (UInt8('A'), UInt8('a'))
            flag = true
        end
    end

    if flag
        fill!(view(x.data, lo:hi), 0x20)
    end
    nothing
end

function CHAR_NA!(x, lo, hi)
    lo, hi = strip_location(x, lo, hi)
    flag = false
    if hi-lo+1 == 2
        if x.data[lo] in (UInt8('N'), UInt8('n')) && x.data[hi] in (UInt8('A'), UInt8('a'))
            flag = true
        end
    end

    if flag
        # for i in lo:hi
        #     x.data[i] = 0x20
        # end
        # x.data[lo] = UInt8('.')
        fill!(view(x.data, lo:hi), 0x20)
    end
    nothing
end

Base.@propagate_inbounds function BOOL!(x, lo, hi)
    lo, hi = strip_location(x, lo, hi)
    if length(lo:hi) == 1
        if x.data[lo] in (0x54, 0x74, 0x31)
            x.data[lo] = 0x31
        elseif x.data[lo] in (0x46, 0x66, 0x30)
            x.data[lo] = 0x30
        end
    elseif length(lo:hi) == 4
        if x.data[lo] in (0x54, 0x74) && x.data[lo+1] in (0x52,0x72) && x.data[lo+2] in (0x55, 0x75) && x.data[lo+3] in (0x45, 0x65)
            x.data[lo] = 0x31
            # x.data[lo+1] = 0x20
            # x.data[lo+2] = 0x20
            # x.data[lo+3] = 0x20
            fill!(view(x.data, lo+1:lo+3), 0x20)
        end
    elseif length(lo:hi) == 5
        if x.data[lo] in (0x46, 0x66) && x.data[lo+1] in (0x41,0x61) && x.data[lo+2] in (0x4c, 0x6c) && x.data[lo+3] in (0x53, 0x73) && x.data[lo+4] in (0x45, 0x65)
            x.data[lo] = 0x30
            # x.data[lo+1] = 0x20
            # x.data[lo+2] = 0x20
            # x.data[lo+3] = 0x20
            # x.data[lo+4] = 0x20
            fill!(view(x.data, lo+1:lo+4), 0x20)
        end
    end
    nothing
end


text:
function (fn::GlobalSchedulerFunction)(thunk)
    promise = Channel{GenericTask}(1)
    function bg_waiter()  # aka bgw
        @record(:bgw_begin)
        try
            local task = Schedulers.spawn(thunk)
            trywait(task)
            @record(:bgw_done, task = taskid(task))
            @check isfinished(task)
            put!(promise, task)
        catch err
            @record(:bgw_error, exception = (err, catch_backtrace()))
            @error(
                "Unexpected error in background waiter",
                exception = (err, catch_backtrace())
            )
        finally
            close(promise)
        end
        @record(:bgw_end)
    end
    Schedulers.spawn(bg_waiter, schedulerof(fn))
    task = take!(promise)
    @check isfinished(task)
    return fetch(task)
end

const DEFAULT_WORKSTEALING_SCHEDULER_HANDLE =
    Ref{Union{WorkStealingScheduler,Nothing}}(nothing)
const DEFAULT_PRIORITIZED_SCHEDULER_HANDLE =
    Ref{Union{PrioritizedScheduler,Nothing}}(nothing)

schedulerof(::typeof(Schedulers.global_workstealing)) =
    DEFAULT_WORKSTEALING_SCHEDULER_HANDLE[]::WorkStealingScheduler
schedulerof(::typeof(Schedulers.global_prioritized)) =
    DEFAULT_PRIORITIZED_SCHEDULER_HANDLE[]::PrioritizedScheduler

close_if_something(x) = close(something(x))
close_if_something(::Nothing) = nothing

function init_global_schedulers()
    close_if_something(DEFAULT_WORKSTEALING_SCHEDULER_HANDLE[])
    close_if_something(DEFAULT_PRIORITIZED_SCHEDULER_HANDLE[])
    DEFAULT_WORKSTEALING_SCHEDULER_HANDLE[] = start(Schedulers.workstealing())
    DEFAULT_PRIORITIZED_SCHEDULER_HANDLE[] = start(Schedulers.prioritized())
    atexit() do
        close_if_something(DEFAULT_WORKSTEALING_SCHEDULER_HANDLE[])
        close_if_something(DEFAULT_PRIORITIZED_SCHEDULER_HANDLE[])
        DEFAULT_WORKSTEALING_SCHEDULER_HANDLE[] = nothing
        DEFAULT_PRIORITIZED_SCHEDULER_HANDLE[] = nothing
    end
end

ext: jl


text:
type Model
    env::Env # Cplex environment
    lp::Ptr{Void} # Cplex problem (lp)
    has_int::Bool # problem has integer variables?
    has_qc::Bool # problem has quadratic constraints?
    callback::Any
    terminator::Vector{Cint}
end

function Model(env::Env, lp::Ptr{Void})
    model = Model(env, lp, false, false, nothing, Cint[0])
    finalizer(model, free_problem)
    set_terminate(model)
    model
end

function Model(env::Env, name::String)
    @assert is_valid(env)
    stat = Array(Cint, 1)
    tmp = @cpx_ccall(createprob, Ptr{Void}, (Ptr{Void}, Ptr{Cint}, Ptr{Cchar}), env.ptr, stat, name)
    if tmp == C_NULL
        throw(CplexError(env, stat))
    end
    return Model(env, tmp)
end

# internal function that wraps finalizer for model and environment together
function _Model(env::Env)
    @assert is_valid(env)
    stat = Array(Cint, 1)
    tmp = @cpx_ccall(createprob, Ptr{Void}, (Ptr{Void}, Ptr{Cint}, Ptr{Cchar}), env.ptr, stat, "CPLEX.jl")
    if tmp == C_NULL
        throw(CplexError(model.env, stat))
    end
    model = Model(env, tmp, false, false, nothing, Cint[0])
    finalizer(model, model::Model -> begin
            tmp = model.env
            free_problem(model)
            close_CPLEX(tmp)
        end)
    set_terminate(model)
    return model
end

function read_model(model::Model, filename::String)
    stat = @cpx_ccall(readcopyprob, Cint, (Ptr{Void}, Ptr{Void}, Ptr{Cchar}, Ptr{Cchar}), model.env.ptr, model.lp, filename, C_NULL)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
end

function write_model(model::Model, filename::String)
    if endswith(filename,".mps")
        filetype = "MPS"
    elseif endswith(filename,".lp")
        filetype = "LP"
    else
        error("Unrecognized file extension: $filename (Only .mps and .lp are supported)")
    end
    stat = @cpx_ccall(writeprob, Cint, (Ptr{Void}, Ptr{Void}, Ptr{Cchar}, Ptr{Cchar}), model.env.ptr, model.lp, filename, filetype)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
end

## TODO: deep copy model, reset model

function get_sense(model::Model)
    sense_int = @cpx_ccall(getobjsen, Cint, (
                           Ptr{Void},
                           Ptr{Void},
                           ),
                           model.env.ptr, model.lp)
    if sense_int == 1
        sense = :Min
    elseif sense_int == -1
        sense = :Max
    else
        error("CPLEX: problem object or environment does not exist")
    end
    return sense
end

function set_sense!(model::Model, sense)
    if sense == :Min
        @cpx_ccall(chgobjsen, Void, (Ptr{Void}, Ptr{Void}, Cint), model.env.ptr, model.lp, 1)
    elseif sense == :Max
        @cpx_ccall(chgobjsen, Void, (Ptr{Void}, Ptr{Void}, Cint), model.env.ptr, model.lp, -1)
    else
        error("Unrecognized objective sense $sense")
    end
end

function get_obj(model::Model)
    nvars = num_var(model)
    obj = Array(Cdouble, nvars)
    stat = @cpx_ccall(getobj, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Ptr{Cdouble},
                      Cint,
                      Cint
                      ),
                      model.env.ptr, model.lp, obj, 0, nvars-1)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
    return obj
end

const type_map = Dict(
     0 => :LP,
     1 => :MILP,
     3 => :MILP,
     5 => :QP,
     7 => :MIQP,
     8 => :MIQP,
    10 => :QCP,
    11 => :MIQCP
)

function get_prob_type(model::Model)
  ret = @cpx_ccall(getprobtype, Cint, (
                   Ptr{Void},
                   Ptr{Void}),
                   model.env.ptr, model.lp)
  ret == -1 && error("No problem of environment")
  return type_map[Int(ret)]
end

function set_obj!(model::Model, c::Vector)
    nvars = num_var(model)
    stat = @cpx_ccall(chgobj, Cint, (
                        Ptr{Void},
                        Ptr{Void},
                        Cint,
                        Ptr{Cint},
                        Ptr{Cdouble}
                        ),
                        model.env.ptr, model.lp, nvars, Cint[0:nvars-1;], float(c))
    if stat != 0
        throw(CplexError(model.env, stat))
    end
end

set_warm_start!(model::Model, x::Vector{Float64}) = set_warm_start!(model, Cint[1:length(x);], x)

function set_warm_start!(model::Model, indx::IVec, val::FVec)
    stat = @cpx_ccall(addmipstarts, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Cint,
                      Cint,
                      Ptr{Cint},
                      Ptr{Cint},
                      Ptr{Cdouble},
                      Ptr{Cint},
                      Ptr{Ptr{Cchar}}
                      ),
                      model.env.ptr, model.lp, 1, length(indx), Cint[0], indx.-1, val, Cint[CPX_MIPSTART_AUTO], C_NULL)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
end

function free_problem(model::Model)
    tmp = Ptr{Void}[model.lp]
    stat = @cpx_ccall(freeprob, Cint, (Ptr{Void}, Ptr{Void}), model.env.ptr, tmp)
    if stat != 0
        throw(CplexError(model.env, stat))
    end
end

function close_CPLEX(env::Env)
    tmp = Ptr{Void}[env.ptr]
    stat = @cpx_ccall(closeCPLEX, Cint, (Ptr{Void},), tmp)
    if stat != 0
        throw(CplexError(env, stat))
    end
end

function set_terminate(model::Model)
    stat = @cpx_ccall(setterminate, Cint, (Ptr{Void},Ptr{Cint}), model.env.ptr, model.terminator)
    if stat != 0
        throw(CplexError(env, stat))
    end
end

terminate(model::Model) = (model.terminator[1] = 1)


text:
function LinearConstraint(arg0::RealVector, arg1::Relationship, arg2::jdouble)
    return LinearConstraint((RealVector, Relationship, jdouble), arg0, arg1, arg2)
end

function LinearConstraint(arg0::RealVector, arg1::jdouble, arg2::Relationship, arg3::RealVector, arg4::jdouble)
    return LinearConstraint((RealVector, jdouble, Relationship, RealVector, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function LinearConstraint(arg0::Vector{jdouble}, arg1::Relationship, arg2::jdouble)
    return LinearConstraint((Vector{jdouble}, Relationship, jdouble), arg0, arg1, arg2)
end

function LinearConstraint(arg0::Vector{jdouble}, arg1::jdouble, arg2::Relationship, arg3::Vector{jdouble}, arg4::jdouble)
    return LinearConstraint((Vector{jdouble}, jdouble, Relationship, Vector{jdouble}, jdouble), arg0, arg1, arg2, arg3, arg4)
end

function equals(obj::LinearConstraint, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_coefficients(obj::LinearConstraint)
    return jcall(obj, "getCoefficients", RealVector, ())
end

function get_relationship(obj::LinearConstraint)
    return jcall(obj, "getRelationship", Relationship, ())
end

function get_value(obj::LinearConstraint)
    return jcall(obj, "getValue", jdouble, ())
end

function hash_code(obj::LinearConstraint)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
pd1 = ProbabilityDistribution(["head","tails"],[0.5,0.5])
pd2 = ProbabilityDistribution(["head","tails"],[0.3,0.7])

a1 = DBayesianNode(:hulu, pd1)
a2 = DBayesianNode(:bulu, pd2)
b3 = CBayesianNode(:esben, x->x^0.5)

e1 = BayesianEdge(a1,a2)
e2 = BayesianEdge(a1,b3)

@test edge_index(e1) == 0
@test edge_index(e2) == 0
@test source(e1) == a1
@test target(e1) == a2
@test target(e2) == b3
@test_throws e2 = BayesianEdge(b3,a1)

k1 = DBayesianNode(:bjarke, pd1)
k2 = DBayesianNode(:thomas, pd2)
k3 = DBayesianNode(:esben, pd1)

h1 = BayesianEdge(k1,k2)
h2 = BayesianEdge(k1,k2)
h3 = BayesianEdge(k2,k1)
h4 = BayesianEdge(k1,k3)

@test h1 == h1
@test h1 == h2
@test h1 != h3
@test h3 != h4
@test h1 != h4

h1.index = 4
h2.index = 3
h3.index = 2
h4.index = 1

@test h1 == h1
@test h1 != h2
@test h1 != h3
@test h3 != h4
@test h1 != h4


text:
using POMDPFiles
using POMDPModels
using Test

pomdp = TigerPOMDP()
write(stdout, pomdp)


text: # Asian Option pricing with Monte Carlo simulation.


text:
@test_nowarn GibbsTypePriors.expected_number_of_clusters_Dirichlet(500, 60.)
@test_nowarn GibbsTypePriors.expected_number_of_clusters_Dirichlet(500, 60., 50)
@test_nowarn GibbsTypePriors.expected_number_of_cluster_2PD(10, 0., 0.4, 20)
@test_nowarn GibbsTypePriors.expected_number_of_cluster_2PD(10, 0., 0.4)
@test_nowarn GibbsTypePriors.expected_number_of_cluster_stable(10, 0.4, 20)
@test_nowarn GibbsTypePriors.expected_number_of_cluster_stable(10, 0.4)
@test_nowarn GibbsTypePriors.expected_number_of_clusters_PD_direct(100, 10., 0.1)
@test_nowarn GibbsTypePriors.expected_number_of_clusters_PD_direct(100, 0.01, 0.9)
@test_nowarn GibbsTypePriors.variance_number_of_clusters_PD_direct(100, 0.01, 0.9)
@test_nowarn GibbsTypePriors.variance_number_of_clusters_PD_direct(1000, 0.01, 0.1)

@test_nowarn GibbsTypePriors.E_Dirichlet(500, 60.)
@test_nowarn GibbsTypePriors.E_2PD(10, 0., 0.4)
@test_nowarn GibbsTypePriors.E_PY(10, 0., 0.4)
@test_nowarn GibbsTypePriors.E_stable(10, 0.4)

@testset "Expectation DP" begin
    @test_nowarn GibbsTypePriors.E_Dirichlet(500, 60.)
    @test_nowarn GibbsTypePriors.E_Dirichlet(100, 0.01)
    @test_nowarn GibbsTypePriors.expected_number_of_clusters_Dirichlet(500, 60.)
    @test_nowarn GibbsTypePriors.expected_number_of_clusters_Dirichlet_direct(1000, 60.)
    @test_nowarn GibbsTypePriors.expected_number_of_clusters_Dirichlet_direct(100, 0.01)
    @test_nowarn GibbsTypePriors.expected_number_of_clusters_Dirichlet_Multinomial_direct(100, 100, 0.01)
    @test_nowarn GibbsTypePriors.expected_number_of_clusters_Dirichlet_Multinomial_direct(1000, 1000, 0.01)
    @test_nowarn GibbsTypePriors.expected_number_of_clusters_Dirichlet_Multinomial(100, 100, 0.01)
    @test_nowarn GibbsTypePriors.expected_number_of_clusters_Dirichlet_Multinomial(100, 100, 0.1)
        @testset "Comparison between methods" begin
        @test Float64(GibbsTypePriors.expected_number_of_clusters_DPM(100, 100, 1.0)) ≈ Float64(GibbsTypePriors.expected_number_of_clusters_Dirichlet_Multinomial(100, 100, 1.0))
        @test_nowarn GibbsTypePriors.expected_number_of_clusters_Dirichlet(500, 60., 150) < @test_nowarn GibbsTypePriors.expected_number_of_clusters_Dirichlet(500, 60.)
        @test Float64(GibbsTypePriors.expected_number_of_clusters_Dirichlet(500, 60.)) ≈  Float64(GibbsTypePriors.E_Dirichlet(500, 60.))
        @test Float64(GibbsTypePriors.expected_number_of_clusters_Dirichlet_direct(100, 6.)) ≈  Float64(GibbsTypePriors.E_Dirichlet(100, 6.))
        @test Float64(GibbsTypePriors.expected_number_of_clusters_Dirichlet_Multinomial_direct(500,500, 6.)) ==  Float64(GibbsTypePriors.expected_number_of_clusters_Dirichlet_Multinomial(500,500,  6.))
        @test Float64(GibbsTypePriors.expected_number_of_clusters_Dirichlet_Multinomial_direct(500,100, 6.)) <  Float64(GibbsTypePriors.E_Dirichlet(500,  6.))
        @test Float64(GibbsTypePriors.expected_number_of_clusters_Dirichlet_Multinomial(500,300, 6.)) <  Float64(GibbsTypePriors.E_Dirichlet(500,  6.))
    end
end

@testset "Expectation PY" begin
    @test_nowarn GibbsTypePriors.E_PY(500, 60., 0.1)
    @test_nowarn GibbsTypePriors.E_PY(100, 60., 0.9)
    @test_nowarn GibbsTypePriors.expected_number_of_clusters_PD_direct(100, 60., 0.9)
    @test_nowarn GibbsTypePriors.E_PY(500, 0.1, 0.99)
    @test_nowarn GibbsTypePriors.expected_number_of_clusters_PYM(50,50,  0.1, 0.5)
    @test_nowarn GibbsTypePriors.expected_number_of_clusters_PYM(50,50,  0.1, 0.5, 50)
    @testset "Comparison between methods" begin
        @test Float64(GibbsTypePriors.expected_number_of_clusters_PYM(100,100,  60., 0.1))<  Float64(GibbsTypePriors.expected_number_of_clusters_PD_direct(100, 60., 0.1))
        @test Float64(GibbsTypePriors.E_PY(100, 60., 0.1)) ≈  Float64(GibbsTypePriors.expected_number_of_clusters_PD_direct(100, 60., 0.1))
        @test Float64(GibbsTypePriors.expected_number_of_cluster_2PD(10, 0., 0.4)) == Float64(GibbsTypePriors.E_PY(10, 0., 0.4))
        @test Float64(GibbsTypePriors.E_PY(100, 6., 0.9)) ≈  Float64(GibbsTypePriors.expected_number_of_clusters_PD_direct(100, 6., 0.9))
        @test Float64(GibbsTypePriors.E_2PD(100, 6., 0.9, 50)) < Float64(GibbsTypePriors.E_PY(100, 6.0, 0.9))
        @test Float64(GibbsTypePriors.E_PY(100, 60., 0.1)) ≈  Float64(GibbsTypePriors.expected_number_of_clusters_PD_direct(100, 60., 0.1))

    end
end

@testset "Expectation NGG" begin
    @test_nowarn GibbsTypePriors.expected_number_of_clusters_NGG(50, 60., 0.1, 50)
    @test_nowarn GibbsTypePriors.expected_number_of_clusters_NGG(50, 60., 0.1)
    @test_nowarn GibbsTypePriors.expected_number_of_clusters_NGGM(50, 50, 60., 0.1,10)
    @testset "Comparison between methods" begin
        @test Float64(GibbsTypePriors.expected_number_of_clusters_NGGM(30, 50, 60., 0.1))< Float64(GibbsTypePriors.expected_number_of_clusters_NGG(30, 60., 0.1))
        @test Float64(GibbsTypePriors.expected_number_of_clusters_NGG(20, 60., 0.1, 10))< Float64(GibbsTypePriors.expected_number_of_clusters_NGG(20, 60., 0.1))
        @test Float64(GibbsTypePriors.E_PY(100, 60., 0.1)) ≈  Float64(GibbsTypePriors.expected_number_of_clusters_PD_direct(100, 60., 0.1))
        @test Float64(GibbsTypePriors.expected_number_of_cluster_2PD(10, 0., 0.4)) == Float64(GibbsTypePriors.E_PY(10, 0., 0.4))
        @test Float64(GibbsTypePriors.E_PY(100, 6., 0.9)) ≈  Float64(GibbsTypePriors.expected_number_of_clusters_PD_direct(100, 6., 0.9))
        @test Float64(GibbsTypePriors.E_2PD(100, 6., 0.9, 50)) < Float64(GibbsTypePriors.E_PY(100, 6.0, 0.9))

    end
end

@testset "Variance for Gibbs-type priors" begin
    @testset "Variance DP" begin
        @test_nowarn GibbsTypePriors.V_Dirichlet(100, 60.)
        @test_nowarn GibbsTypePriors.variance_number_of_clusters_Dirichlet(100, 0.6)
        @test_nowarn GibbsTypePriors.variance_number_of_clusters_Dirichlet_direct(1000, 60.)
        @test_nowarn GibbsTypePriors.variance_number_of_clusters_Dirichlet_direct(1000, 0.001)
        @test_nowarn GibbsTypePriors.variance_number_of_clusters_DPM(100, 100, 0.1)
        @test_nowarn GibbsTypePriors.variance_number_of_clusters_Dirichlet(100, 0.1)
        @testset "Comparison" begin
            @test Float64(GibbsTypePriors.variance_number_of_clusters_DPM(100, 100, 0.1)) < Float64(GibbsTypePriors.V_Dirichlet(100, 0.1))
            @test Float64(GibbsTypePriors.variance_number_of_clusters_Dirichlet(100, 60.)) ≈ Float64(GibbsTypePriors.V_Dirichlet(100, 60.))
            @test Float64(GibbsTypePriors.variance_number_of_clusters_Dirichlet_direct(100, 60.)) ≈ Float64(GibbsTypePriors.V_Dirichlet(100, 60.))
        end
    end
    @testset "Variance PY" begin
        @test_nowarn GibbsTypePriors.V_PY(500, 60., 0.1)
        @test_nowarn GibbsTypePriors.V_PY(100, 60., 0.9)
        @test_nowarn GibbsTypePriors.V_PY(1000, 60., 0.9)
        @test_nowarn GibbsTypePriors.V_PY(100, 60., 0.5)
        @test_nowarn GibbsTypePriors.variance_number_of_clusters_PD(20, 5., 0.9)
        @test_nowarn GibbsTypePriors.variance_number_of_clusters_PD_direct(20, 5., 0.9)
        @test_nowarn GibbsTypePriors.variance_number_of_clusters_PD(20, 5., 0.9)
        @test_nowarn GibbsTypePriors.variance_number_of_clusters_PYM(20,50,  5., 0.9)
        @testset "Comparison between methods" begin
            @test Float64(GibbsTypePriors.variance_number_of_clusters_PD(100, 0.5, 0.5)) ≈ Float64(GibbsTypePriors.V_PY(100, 0.5, 0.5))
            @test Float64(GibbsTypePriors.variance_number_of_clusters_PD(100, 0.5, 0.5)) ≈ Float64(GibbsTypePriors.variance_number_of_clusters_PD_direct(100, 0.5, 0.5))
        end
    end
    @testset "Variance NGG" begin
        @test_nowarn GibbsTypePriors.variance_number_of_clusters_NGG(50, 60., 0.1)
        @test_nowarn GibbsTypePriors.variance_number_of_clusters_NGG(10, 6., 0.8)
        @test_nowarn GibbsTypePriors.variance_number_of_clusters_NGGM(10, 10, 6., 0.8)
    end
end

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


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

struct LoggedStaticPool <: AbstractThreadPool
    tids :: Vector{Int}
    lck  :: ReentrantLock
    recs :: Vector{ThreadRecord}
    t0   :: UInt64
    log  :: ThreadLog
end

#############################
# Constructors / Finalizer
#############################

"""
    LoggedStaticPool(init_thrd=1, nthrds=Threads.nthreads())

The main LoggedStaticPool object.
"""
function LoggedStaticPool(init_thrd::Integer=1, nthrds::Integer=Threads.nthreads())
    thrd0 = min(init_thrd, Threads.nthreads())
    thrd1 = min(thrd0+nthrds-1, Threads.nthreads())
    return LoggedStaticPool(thrd0:thrd1, ReentrantLock(), ThreadRecord[], time_ns(), ThreadLog())
end

function Base.close(pool::LoggedStaticPool)
    _recordstolog!(pool.log, pool.recs, pool. t0)
end

#############################
# ThreadPool Interface
#############################

function tmap(fn, pool::LoggedStaticPool, itr)
    data = collect(itr)
    applicable(fn, first(data)) || error("function can't be applied to iterator contents")
    N = length(data)
    sizehint!(pool.recs, N)
    result = Array{_detect_type(fn, data), ndims(data)}(undef, size(data))
    nthrds = length(pool.tids)
    njobs = div(N,nthrds)
    remjobs = N % nthrds

    len(ind) = max(0, njobs + (nthrds-ind+1 <= remjobs ? 1 : 0))
    finish(ind) = sum([len(x) for x in 1:ind])
    start(ind) = finish(ind)-len(ind)+1

    _fn(ind) = begin
        if finish(ind) > 0
            tid = Threads.threadid()
            for i in start(ind):finish(ind)
                lock(pool.lck)
                push!(pool.recs, (i, tid, true, time_ns()))
                unlock(pool.lck)
                @inbounds result[i] = fn(Base.unsafe_getindex(data, i))
                lock(pool.lck)
                push!(pool.recs, (i, tid, false, time_ns()))
                unlock(pool.lck)
            end
        end
    end

    Threads.@threads for tid in 1:Threads.nthreads()
        ind = findfirst(t->t==tid, pool.tids)
        isnothing(ind) || _fn(ind)
    end

    return result
end


text:
####
#### Defines list of developer documents
####

dev_docs = Any[
    "Coding style" => "DevDocs/CodeStyle.md",
    "Acceptable Unicode" => "DevDocs/AcceptableUnicode.md",
    "Model variable list" => "DevDocs/ModelVariableList.md",
    "Model output" => "DevDocs/ModelOutput.md",
    "Diagnostic variable list" => "DevDocs/DiagnosticVariableList.md",
    "Custom System Image" => "DevDocs/SystemImage.md",
]


text:
# interfacing with Random standard library

using Random

using GPUArrays

export rand_logn!, rand_poisson!

# the interface is split in two levels:
# - functions that extend the Random standard library, and take an RNG as first argument,
#   will only ever dispatch to CURAND and as a result are limited in the types they support.
# - functions that take an array will dispatch to either CURAND or GPUArrays
# - non-exported functions are provided for constructing GPU arrays from only an eltype

mutable struct RNG <: Random.AbstractRNG
    handle::curandGenerator_t
    ctx::CuContext
    typ::Int

    function RNG(typ=CURAND_RNG_PSEUDO_DEFAULT)
        handle_ref = Ref{curandGenerator_t}()
        @allocates curandCreateGenerator(handle_ref, typ)

        obj = new(handle_ref[], context(), typ)
        finalizer(unsafe_destroy!, obj)
        return obj
    end
end

function unsafe_destroy!(rng::RNG)
    if CUDAdrv.isvalid(rng.ctx)
        curandDestroyGenerator(rng)
    end
end

Base.unsafe_convert(::Type{curandGenerator_t}, rng::RNG) = rng.handle

## seeding

function Random.seed!(rng::RNG, seed=Base.rand(UInt64), offset=0)
    curandSetPseudoRandomGeneratorSeed(rng, seed)
    curandSetGeneratorOffset(rng, offset)
    @allocates curandGenerateSeeds(rng)
    return
end

# convenience function that seeds both generators
function seed!(seed=Base.rand(UInt64))
    Random.seed!(generator(), seed)
    Random.seed!(GPUArrays.global_rng(CuArray, device()), seed)
end

## in-place

# uniform
const UniformType = Union{Type{Float32},Type{Float64}}
const UniformArray = CuArray{<:Union{Float32,Float64}}
function Random.rand!(rng::RNG, A::CuArray{Float32})
    @allocates curandGenerateUniform(rng, A, length(A))
    return A
end
function Random.rand!(rng::RNG, A::CuArray{Float64})
    @allocates curandGenerateUniformDouble(rng, A, length(A))
    return A
end

# some functions need pow2 lengths: use a padded array and copy back to the original one
function inplace_pow2(A, f)
    len = length(A)
    if len > 1 && ispow2(len)
        f(A)
    else
        padlen = max(2, nextpow(2, len))
        B = similar(A, padlen)
        f(B)
        copyto!(A, 1, B, 1, len)
        CuArrays.unsafe_free!(B)
    end
    A
end

# normal
const NormalType = Union{Type{Float32},Type{Float64}}
const NormalArray = CuArray{<:Union{Float32,Float64}}
function Random.randn!(rng::RNG, A::CuArray{Float32}; mean=0, stddev=1)
    inplace_pow2(A, B->@allocates(curandGenerateNormal(rng, B, length(B), mean, stddev)))
    return A
end
function Random.randn!(rng::RNG, A::CuArray{Float64}; mean=0, stddev=1)
    inplace_pow2(A, B->@allocates(curandGenerateNormalDouble(rng, B, length(B), mean, stddev)))
    return A
end

# log-normal
const LognormalType = Union{Type{Float32},Type{Float64}}
const LognormalArray = CuArray{<:Union{Float32,Float64}}
function rand_logn!(rng::RNG, A::CuArray{Float32}; mean=0, stddev=1)
    inplace_pow2(A, B->@allocates(curandGenerateLogNormal(rng, B, length(B), mean, stddev)))
    return A
end
function rand_logn!(rng::RNG, A::CuArray{Float64}; mean=0, stddev=1)
    inplace_pow2(A, B->@allocates(curandGenerateLogNormalDouble(rng, B, length(B), mean, stddev)))
    return A
end

# poisson
const PoissonType = Union{Type{Cuint}}
const PoissonArray = CuArray{Cuint}
function rand_poisson!(rng::RNG, A::CuArray{Cuint}; lambda=1)
    @allocates curandGeneratePoisson(rng, A, length(A), lambda)
    return A
end

## out of place

# some functions need pow2 lengths: construct a compatible array and return part of it
function outofplace_pow2(shape, ctor, f)
    len = prod(shape)
    if ispow2(len)
        A = ctor(shape)
        f(A)
    else
        padlen = max(2, nextpow(2, len))
        A = ctor(padlen)
        f(A)
        B = reshape(A[1:len], shape)
        return B
    end
end

# arrays
Random.rand(rng::RNG, T::UniformType, dims::Dims) =
    Random.rand!(rng, CuArray{T}(undef, dims))
Random.randn(rng::RNG, T::NormalType, dims::Dims; kwargs...) =
    outofplace_pow2(dims, shape->CuArray{T}(undef, dims), A->randn!(rng, A; kwargs...))
rand_logn(rng::RNG, T::LognormalType, dims::Dims; kwargs...) =
    outofplace_pow2(dims, shape->CuArray{T}(undef, dims), A->rand_logn!(rng, A; kwargs...))
rand_poisson(rng::RNG, T::PoissonType, dims::Dims; kwargs...) =
    rand_poisson!(rng, CuArray{T}(undef, dims); kwargs...)

# specify default types
Random.rand(rng::RNG, dims::Dims; kwargs...) = Random.rand(rng, Float32, dims; kwargs...)
Random.randn(rng::RNG, dims::Dims; kwargs...) = Random.randn(rng, Float32, dims; kwargs...)
rand_logn(rng::RNG, dims::Dims; kwargs...) = rand_logn(rng, Float32, dims; kwargs...)
rand_poisson(rng::RNG, dims::Dims; kwargs...) = rand_poisson(rng, Cuint, dims; kwargs...)

# support all dimension specifications
Random.rand(rng::RNG, T::UniformType, dim1::Integer, dims::Integer...) =
    Random.rand(rng, T, Dims((dim1, dims...)))
Random.randn(rng::RNG, T::NormalType, dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn(rng, T, Dims((dim1, dims...)); kwargs...)
rand_logn(rng::RNG, T::LognormalType, dim1::Integer, dims::Integer...; kwargs...) =
    rand_logn(rng, T, Dims((dim1, dims...)); kwargs...)
rand_poisson(rng::RNG, T::PoissonType, dim1::Integer, dims::Integer...; kwargs...) =
    rand_poisson(rng, T, Dims((dim1, dims...)); kwargs...)

## functions that dispatch to either CURAND or GPUArrays

# CURAND in-place
Random.rand!(A::UniformArray) = Random.rand!(generator(), A)
Random.randn!(A::NormalArray; kwargs...) = Random.randn!(generator(), A; kwargs...)
rand_logn!(A::LognormalArray; kwargs...) = rand_logn!(generator(), A; kwargs...)
rand_poisson!(A::PoissonArray; kwargs...) = rand_poisson!(generator(), A; kwargs...)

# CURAND out-of-place
rand(T::UniformType, dims::Dims) = Random.rand(generator(), T, dims)
randn(T::NormalType, dims::Dims; kwargs...) = Random.randn(generator(), T, dims; kwargs...)
rand_logn(T::LognormalType, dims::Dims; kwargs...) = rand_logn(generator(), T, dims; kwargs...)
rand_poisson(T::PoissonType, dims::Dims; kwargs...) = rand_poisson(generator(), T, dims; kwargs...)

# support all dimension specifications
rand(T::UniformType, dim1::Integer, dims::Integer...) =
    Random.rand(generator(), T, Dims((dim1, dims...)))
randn(T::NormalType, dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn(generator(), T, Dims((dim1, dims...)); kwargs...)
rand_logn(T::LognormalType, dim1::Integer, dims::Integer...; kwargs...) =
    rand_logn(generator(), T, Dims((dim1, dims...)); kwargs...)
rand_poisson(T::PoissonType, dim1::Integer, dims::Integer...; kwargs...) =
    rand_poisson(generator(), T, Dims((dim1, dims...)); kwargs...)

# GPUArrays in-place
Random.rand!(A::CuArray) = Random.rand!(GPUArrays.global_rng(A), A)
Random.randn!(A::CuArray; kwargs...) =
    error("CuArrays does not support generating normally-distrubyted random numbers of type $(eltype(A))")
rand_logn!(A::CuArray; kwargs...) =
    error("CuArrays does not support generating lognormally-distributed random numbers of type $(eltype(A))")
rand_poisson!(A::CuArray; kwargs...) =
    error("CuArrays does not support generating Poisson-distributed random numbers of type $(eltype(A))")

# GPUArrays out-of-place
rand(T::Type, dims::Dims) = Random.rand!(CuArray{T}(undef, dims...))
randn(T::Type, dims::Dims; kwargs...) = Random.randn!(CuArray{T}(undef, dims...); kwargs...)
rand_logn(T::Type, dims::Dims; kwargs...) = rand_logn!(CuArray{T}(undef, dims...); kwargs...)
rand_poisson(T::Type, dims::Dims; kwargs...) = rand_poisson!(CuArray{T}(undef, dims...); kwargs...)

# support all dimension specifications
rand(T::Type, dim1::Integer, dims::Integer...) =
    Random.rand!(CuArray{T}(undef, dim1, dims...))
randn(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn!(CuArray{T}(undef, dim1, dims...); kwargs...)
rand_logn(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    rand_logn!(CuArray{T}(undef, dim1, dims...); kwargs...)
rand_poisson(T::Type, dim1::Integer, dims::Integer...; kwargs...) =
    rand_poisson!(CuArray{T}(undef, dim1, dims...); kwargs...)

# untyped out-of-place
rand(dim1::Integer, dims::Integer...) =
    Random.rand(generator(), Dims((dim1, dims...)))
randn(dim1::Integer, dims::Integer...; kwargs...) =
    Random.randn(generator(), Dims((dim1, dims...)); kwargs...)
rand_logn(dim1::Integer, dims::Integer...; kwargs...) =
    rand_logn(generator(), Dims((dim1, dims...)); kwargs...)
rand_poisson(dim1::Integer, dims::Integer...; kwargs...) =
    rand_poisson(generator(), Dims((dim1, dims...)); kwargs...)


text:
using Flux

"""
    init!(a, initializer)

Initialize an array `a` with `initializer`.

Example

```julia
layer = Conv((3, 3), 3 => 16)
init!(layer.weight, init_kaiming_normal)
```
"""
init!(a::AbstractArray{T}, initializer) where T = copy!(a, initializer(T, size(a)...))

# initializers

"""
    init_kaiming_uniform(T, dims...; gain=sqrt(2))

Kaiming uniform initializer
"""
function init_kaiming_uniform(T, dims...; gain=sqrt(2))
  fan_in = length(dims) <= 2 ? dims[end] : div(*(dims...), dims[end])
  bound = sqrt(3.0) * gain / sqrt(fan_in) |> Float32
  return rand(Float32, dims...) .* (2bound) .+ -bound
end
init_kaiming_uniform(dims...; gain = sqrt(2)) = init_kaiming_uniform(Float32, dims...; gain = gain)

"""
    init_normal_uniform(T, dims...; gain=sqrt(2))

Kaiming normal initializer
"""
function init_kaiming_normal(T, dims...; gain=sqrt(2))
  fan_in = length(dims) <= 2 ? dims[end] : div(*(dims...), dims[end])
  std = gain / sqrt(fan_in) |> T
  return randn(T, dims...) .* std
end
init_kaiming_normal(dims...; gain = sqrt(2)) = init_kaiming_normal(Float32, dims...; gain = gain)

init_zeros = Flux.zeros

"""
    Initialization(layertype, field, initializer)

Describes an initialization for field `field` of all layers of type
`layertype` to be initialized with `initializer`.

Use [`initmodel!`](@ref) to apply to a model.

## Example

`Initialization(Conv, :bias, init_zeros)` -> initialize the bias of all
Conv layers with zeros.
"""
struct Initialization
    layertype::Type
    field::Symbol
    initializer
end

"""
    initmodel!(model, inits::Vector{Initialization})

Applies `inits` to `model`.

See [`Initialization`](@ref).

## Example

```julia
model = Chain(Dense(10, 10), Dense(10, 10))
initmodel!(model, [
    Initialization(Dense, :W, init_kaiming_normal),
    Initialization(Dense, :bias, init_zeros)
])
```
"""
function initmodel!(model, inits::AbstractVector{Initialization})
    for layer in IterLayers(model)
        for init in inits
            if (layer isa init.layertype) && hasfield(typeof(layer), init.field)
                init!(getfield(layer, init.field), init.initializer)
            end
        end
    end
    return model
end


text:
export Locations, CtrlLocations, AbstractLocations, LocationError, merge_locations

const LocationStorageTypes = Union{Int,NTuple{N,Int} where N,UnitRange{Int},Vector{Int}}

abstract type AbstractLocations end

"""
    merge_locations(locations...)

Construct a new `Locations` by merging two or more existing locations.
"""
merge_locations(x::AbstractLocations, y::AbstractLocations, locations::AbstractLocations...) =
    merge_locations(merge_locations(x, y), locations...)

merge_locations(x::AbstractLocations) = x

"""
    Locations <: AbstractLocations

Type to annotate locations in quantum circuit.

    Locations(x)

Create a `Locations` object from a raw location statement. Valid storage types are:

- `Int`: single position
- `NTuple{N, Int}`: a list of locations
- `UnitRange{Int}`: contiguous locations

Other types will be converted to the storage type via `Tuple`.
"""
struct Locations{T<:LocationStorageTypes} <: AbstractLocations
    storage::T

    Base.@pure function Locations(x::T) where {T<:LocationStorageTypes}
        new{T}(x)
    end
end

# skip it if x is a location
Locations(x::Locations) = x
Locations(xs...) = Locations(xs)
Locations(xs::NTuple{N,<:Locations}) where {N} = merge_locations(xs...)
Locations(x::NTuple{N,T}) where {N,T} = throw(LocationError("expect Int, got $T"))

Base.@propagate_inbounds Base.getindex(l::Locations, idx...) = Locations(getindex(l.storage, idx...))
Base.length(l::Locations) = length(l.storage)
Base.iterate(l::Locations) = iterate(l.storage)
Base.iterate(l::Locations, st) = iterate(l.storage, st)
Base.eltype(::Type{T}) where {T<:Locations} = Int
Base.eltype(x::Locations) = Int
Base.Tuple(x::Locations) = (x.storage...,)

struct LocationError <: Exception
    msg::String
end

function merge_locations(l1::Locations, l2::Locations)
    Locations((l1.storage..., l2.storage...))
end

# location mapping
# TODO: preserve sign when indexing
# TODO: provide a @inlocation macro via Expr(:meta, :inlocation, true) so when we compile to Julia functions
#       we can use unsafe_mapping directly
@inline unsafe_mapping(parent::Locations{Int}, sub::Locations{Int}) = parent
@inline unsafe_mapping(parent::Locations{Int}, sub::Locations{NTuple{N,Int}}) where {N} = parent
@inline unsafe_mapping(parent::Locations{Int}, sub::Locations{UnitRange{Int}}) = parent
@inline unsafe_mapping(parent::Locations{NTuple{N,Int}}, sub::Locations{Int}) where {N} =
    Locations(@inbounds parent[sub.storage])
@inline unsafe_mapping(parent::Locations{NTuple{N,Int}}, sub::Locations{NTuple{M,Int}}) where {N,M} =
    Locations(map(x -> @inbounds(parent[x]), sub.storage))
@inline unsafe_mapping(parent::Locations{NTuple{N,Int}}, sub::Locations{UnitRange{Int}}) where {N} =
    Locations(@inbounds parent[sub.storage])
@inline unsafe_mapping(parent::Locations{UnitRange{Int}}, sub::Locations{Int}) =
    Locations(@inbounds parent[sub.storage])
@inline unsafe_mapping(parent::Locations{UnitRange{Int}}, sub::Locations{NTuple{N,Int}}) where {N} =
    Locations(map(x -> @inbounds(parent[x]), sub.storage))
@inline unsafe_mapping(parent::Locations{UnitRange{Int}}, sub::Locations{UnitRange{Int}}) =
    Locations(@inbounds parent[sub.storage])

@inline map_error(parent, sub) = throw(LocationError("got $sub in parent space $parent"))

@noinline function map_check(parent, sub)
    map_check_nothrow(parent, sub) || map_error(parent, sub)
end

function map_check_nothrow(parent::Locations{Int}, sub::Locations{Int})
    sub.storage == 1
end

function map_check_nothrow(parent::Locations{Int}, sub::Locations{Tuple{Int}})
    sub.storage[1] == 1
end

function map_check_nothrow(parent::Locations{Int}, sub::Locations{NTuple{N,Int}}) where {N}
    false
end

function map_check_nothrow(parent::Locations{Int}, sub::Locations{UnitRange{Int}})
    (length(sub) == 1) && (sub.storage.start == 1)
end

function map_check_nothrow(parent::Locations{NTuple{N,Int}}, sub::Locations{Int}) where {N}
    1 <= sub.storage <= N
end

function map_check_nothrow(
    parent::Locations{NTuple{N,Int}},
    sub::Locations{NTuple{M,Int}},
) where {N,M}
    all(x -> (1 <= x <= N), sub.storage)
end

function map_check_nothrow(parent::Locations{NTuple{N,Int}}, sub::Locations{UnitRange{Int}}) where {N}
    (1 <= sub.storage.start) && (sub.storage.stop <= N)
end

function map_check(parent::Locations{UnitRange{Int}}, sub::Locations{Int})
    1 <= sub.storage <= length(parent)
end

function map_check(parent::Locations{UnitRange{Int}}, sub::Locations{NTuple{N,Int}}) where {N}
    all(x -> (1 <= x <= length(parent)), sub.storage)
end

function map_check(parent::Locations{UnitRange{Int}}, sub::Locations{UnitRange{Int}})
    (1 <= sub.storage.start) && (sub.storage.stop <= length(parent))
end

# comparing
function Base.:(==)(l1::Locations, l2::Locations)
    length(l1) == length(l2) || return false
    return all(l1.storage .== l2.storage)
end

struct CtrlFlags{L,N}
    data::NTuple{N,UInt64}
end

function flags(bits::Vararg{UInt8,L}) where {L}
    N = Base.num_bit_chunks(L)
    data = ntuple(N) do k
        x = UInt64(0)
        p = UInt64(0)
        for i in (64k-63):min(64k, L)
            x += bits[i] << p
            p += UInt64(1)
        end
        return x
    end
    return CtrlFlags{L,N}(data)
end

function _default_flags(L::Int)
    N = Base.num_bit_chunks(L)
    return CtrlFlags{L,N}(ntuple(x -> UInt64(0), N))
end

function Base.all(flags::CtrlFlags)
    return all(iszero, flags.data)
end

function Base.getindex(flags::CtrlFlags, idx::Int)
    i1, i2 = Base.get_chunks_id(idx)
    u = UInt64(1) << i2
    r = (flags.data[i1] & u) == 0
    return r
end

function merge_flags(a::CtrlFlags{LA}, b::CtrlFlags{LB}) where {LA,LB}
    L = LA + LB
    N = Base.num_bit_chunks(L)
    last_chunk_len = rem(LA, 64)

    if iszero(last_chunk_len)
        return CtrlFlags{L,N}((a.data..., b.data...))
    else
        return flags(ntuple(k -> UInt8(!a[k]), LA)..., ntuple(k -> UInt8(!b[k]), LB)...)
    end
end

function Base.show(io::IO, x::CtrlFlags{L,N}) where {L,N}
    print(io, "CtrlFlags(")
    p = L - 1
    for d in x.data
        for k in 0:min(63, p)
            print(io, d >> k & UInt64(1))
        end
        p -= 64
    end
    print(io, ")")
end

# CtrlLocations
struct CtrlLocations{T<:LocationStorageTypes,L,N} <: AbstractLocations
    storage::Locations{T}
    flags::CtrlFlags{L,N}

    Base.@pure CtrlLocations(storage::Locations{T}, flags::CtrlFlags{L,N}) where {T,L,N} =
        new{T,L,N}(storage, flags)
end

# skip itself
CtrlLocations(x::CtrlLocations) = x
CtrlLocations(x::Locations) = CtrlLocations(x, _default_flags(length(x)))
CtrlLocations(x::LocationStorageTypes, configs::NTuple{L,UInt8}) where {L} =
    CtrlLocations(Locations(x), flags(configs...))
CtrlLocations(xs...) = CtrlLocations(Locations(xs...))

Base.length(l::CtrlLocations) = length(l.storage)
Base.iterate(l::CtrlLocations) = iterate(l.storage)
Base.iterate(l::CtrlLocations, st) = iterate(l.storage, st)

function Base.show(io::IO, x::Locations)
    print(io, "Locations(")
    print_locations(io, x)
    print(io, ")")
end

function Base.show(io::IO, x::CtrlLocations)
    print(io, "CtrlLocations(")
    print_locations(io, x)
    print(io, ")")
end

function print_locations(io::IO, x::Locations)
    if x.storage isa AbstractRange
        return printstyled(io, x.storage; color = :light_blue)
    end

    nlocations = length(x)
    for i in 1:nlocations
        printstyled(io, x.storage[i]; color = :light_blue)

        if i != nlocations
            print(io, ", ")
        end
    end
end

function print_locations(io::IO, x::CtrlLocations)
    if all(x.flags)
        print_locations(io, x.storage)
    else
        nlocations = length(x)
        for i in 1:nlocations
            l = x.storage.storage[i]
            if x.flags[i]
                @show "pass"
                printstyled(io, l; color = :light_blue)
            else
                printstyled(io, "!", l; color = :light_blue)
            end

            if i != nlocations
                print(io, ", ")
            end
        end
    end
end

function merge_locations(l1::CtrlLocations, l2::CtrlLocations)
    CtrlLocations(merge_locations(l1.storage, l2.storage), merge_flags(l1.flags, l2.flags))
end

# NOTE: CtrlLocations can not be mapped by Locations
@inline unsafe_mapping(parent::Locations, sub::CtrlLocations) =
    CtrlLocations(unsafe_mapping(parent, sub.storage), sub.flags)
map_check_nothrow(parent::Locations, sub::CtrlLocations) = map_check_nothrow(parent, sub.storage)

Base.:(==)(l1::CtrlLocations, l2::CtrlLocations) =
    (l1.storage == l2.storage) && (l1.flags == l2.flags)

# parent location has to be Locations since CtrlLocations can't be parent
@inline function Base.getindex(parent::Locations, sub::AbstractLocations)
    map_check(parent, sub)
    return unsafe_mapping(parent, sub)
end


text:
module TestsUKBB

using Test
using BGEN
using CSV
using DataFrames
using TMLEEpistasis
using CategoricalArrays
using TOML

include("helper_fns.jl")

@testset "Test read_bgen function" begin
    # This file as an associated .bgi
    bgen_file = BGEN.datadir("example.8bits.bgen")
    b = TMLEEpistasis.read_bgen(bgen_file)
    @test b.idx isa BGEN.Index

    # This file as no associated .bgi
    bgen_file = BGEN.datadir("example.9bits.bgen")
    b = TMLEEpistasis.read_bgen(bgen_file)
    @test b.idx === nothing
end

@testset "Test samples_genotype" begin
    probabilities = [0.3 0.2 0.9;
                     0.5 0.2 0.05;
                     0.2 0.6 0.05]
    variant_genotypes = ["AA", "AT", "TT"]

    threshold = 0.9
    genotypes = TMLEEpistasis.samples_genotype(
        probabilities,
        variant_genotypes,
        threshold)
    @test genotypes[1] === genotypes[2] === missing
    @test genotypes[3] == "AA"

    threshold = 0.55
    genotypes = TMLEEpistasis.samples_genotype(
        probabilities,
        variant_genotypes,
        threshold)
    @test genotypes[1]  === missing
    @test genotypes[2] == "TT"
    @test genotypes[3] == "AA"

end

@testset "Test UKBBGenotypes function" begin
    build_query_file(threshold=0.95)
    queries = TMLEEpistasis.parse_queries(queryfile)
    genotypes = TMLEEpistasis.UKBBGenotypes(queryfile, queries)
    # I only look at the first 10 rows
    # SAMPLE_ID
    @test genotypes[1:9, "SAMPLE_ID"] == ["sample_00$i" for i in 1:9]
    # RSID_10
    @test genotypes[1:10, "RSID_10"] == repeat(["AG"], 10)
    # RSID_100
    @test all(genotypes[1:10, "RSID_100"] .=== ["AG", "AA", "AG", missing, "AG", "AG", missing, "AG", "GG", "AG"])

    rm(queryfile)
end

@testset "Test variant_genotypes" begin
    b = TMLEEpistasis.read_bgen(BGEN.datadir("example.8bits.bgen"))
    queries = [
        "query_1" => (RSID_10=["AA", "GG"], RSID_100=["AA", "GG"]),
        "query_2" => (RSID_10=["AA", "GG"], RSID_100=["AA", "GA"])
    ]
    # Defaults to "AG"
    v = variant_by_rsid(b, "RSID_10")
    @test TMLEEpistasis.variant_genotypes(v, queries) == ["AA", "AG", "GG"]
    # "GA" is specified by user so returned
    v = variant_by_rsid(b, "RSID_100")
    @test TMLEEpistasis.variant_genotypes(v, queries) == ["AA", "GA", "GG"]

end

@testset "Test preprocess" begin
    confounders = CSV.File(confoundersfile) |> DataFrame
    n = size(confounders)[1]
    genotypes = DataFrame(
        SAMPLE_ID = confounders.SAMPLE_ID,
        RSID_1    = vcat(Vector{Missing}(missing, 10), collect(1:n-10)),
        RSID_2    = vcat(collect(1:n-10), Vector{Missing}(missing, 10))
        )

    # Continuous phenotypes
    phenotypes = CSV.File(phenotypefile, select=["eid", "continuous_phenotype"]) |> DataFrame
    T, W, y = TMLEEpistasis.preprocess(genotypes, confounders, phenotypes;verbosity=0)

    @test y == phenotypes.continuous_phenotype[11:n-10]

    @test size(T) == (480, 2)
    @test size(W) == (480, 2)
    @test size(y) == (480,)

    @test names(T) == ["RSID_1", "RSID_2"]
    @test names(W) == ["PC1", "PC2"]

    @test T.RSID_1 == collect(1:n-20)
    @test T.RSID_1 isa CategoricalArray
    @test T.RSID_2 == collect(11:n-10)
    @test T.RSID_2 isa CategoricalArray

    # Binary phenotypes
    phenotypes = CSV.File(phenotypefile, select=["eid", "categorical_phenotype"]) |> DataFrame
    T, W, y = TMLEEpistasis.preprocess(genotypes, confounders, phenotypes;
                verbosity=0)

    @test y isa CategoricalArray{Bool,1,UInt32}
end

@testset "Test PhenotypeTMLEEpistasis with continuous target" begin
    # Only one continuous phenotype
    estimatorfile = joinpath("config", "tmle_config.toml")
    build_query_file()
    parsed_args = Dict(
        "phenotypes" => phenotypefile,
        "confounders" => confoundersfile,
        "queries" => queryfile,
        "estimator" => estimatorfile,
        "output" => "cont_results.csv",
        "phenotypes-list" => "data/phen_list_1.csv",
        "verbosity" => 0
    )

    UKBBVariantRun(parsed_args)

    results = CSV.File(parsed_args["output"]) |> DataFrame
    @test results.QUERYNAME == ["QUERY_1", "QUERY_2"]
    @test names(results) == ["PHENOTYPE", "QUERYNAME", "QUERYSTRING", "ESTIMATE", "PVALUE", "LOWER_BOUND", "UPPER_BOUND", "STD_ERROR", "QSTACK_COEFS", "NROWS", "TCOUNTS"]
    @test size(results) == (2, 11)

    # Clean
    rm(parsed_args["output"])
    rm(parsed_args["queries"])
end

@testset "Test PhenotypeTMLEEpistasis with categorical target and recovery mode" begin
    estimatorfile = joinpath("config", "tmle_config.toml")
    build_query_file()
    parsed_args = Dict(
        "phenotypes" => phenotypefile,
        "confounders" => confoundersfile,
        "queries" => queryfile,
        "estimator" => estimatorfile,
        "output" => "full_results.csv",
        "phenotypes-list" => "data/phen_list_2.csv",
        "verbosity" => 0
    )

    initial_results = DataFrame(PHENOTYPE=[:continuous_phenotype, :continuous_phenotype],
                    QUERYNAME=[:QUERY_DONE, :QUERY_DONE],
                    QUERYSTRING=["QUERYSTRING_1", "QUERYSTRING_2"],
                    ESTIMATE=[1., 2.],
                    PVALUE=[1., 2.],
                    LOWER_BOUND=[1., 2.],
                    UPPER_BOUND=[1., 2.],
                    STD_ERROR=[1., 2.],
                    QSTACK_COEFS=["TOTO", "TOTO"],
                    NROWS=[0, 0],
                    TCOUNTS=["TOTO", "TOTO"]
                    )
    CSV.write(parsed_args["output"], initial_results)

    UKBBVariantRun(parsed_args)

    results = CSV.File(parsed_args["output"]) |> DataFrame

    @test results[!, :NROWS] == [0, 0, 477, 477]
    @test results[!, :TCOUNTS] == ["TOTO", "TOTO", "RSID_10 & RSID_100: AG AA 65 | AG AG 248 | AG GG 145 | GG AA 7 | GG AG 4 | GG GG 8", "RSID_10 & RSID_100: AG AA 65 | AG AG 248 | AG GG 145 | GG AA 7 | GG AG 4 | GG GG 8"]
    @test results.QUERYNAME == ["QUERY_DONE", "QUERY_DONE", "QUERY_1", "QUERY_2"]
    @test names(results) == ["PHENOTYPE", "QUERYNAME", "QUERYSTRING", "ESTIMATE", "PVALUE", "LOWER_BOUND", "UPPER_BOUND", "STD_ERROR", "QSTACK_COEFS", "NROWS", "TCOUNTS"]
    @test size(results) == (4, 11)

    # Clean
    rm(parsed_args["output"])
    rm(parsed_args["queries"])
end

@testset "Test PhenotypeCrossValidation with categorical target and recovery mode" begin
    estimatorfile = joinpath("config", "tmle_config.toml")
    build_query_file()
    parsed_args = Dict(
        "phenotypes" => phenotypefile,
        "confounders" => confoundersfile,
        "queries" => queryfile,
        "estimator" => estimatorfile,
        "output" => "full_results.csv",
        "phenotypes-list" => "data/phen_list_2.csv",
        "verbosity" => 0
    )

    initial_results = DataFrame(PHENOTYPE=[:continuous_phenotype],
                    VARIANTS=["nothing"],
                    Q_RESULTSTRING=["nothing"],
                    G_RESULTSTRING=["nothing"],
                    )
    CSV.write(parsed_args["output"], initial_results)

    UKBBVariantRun(parsed_args, run_fn=TMLEEpistasis.PhenotypeCrossValidation)

    results = CSV.File(parsed_args["output"]) |> DataFrame

    @test results[2, "VARIANTS"] == "RSID_10 && RSID_100"
    # Check Q
    Qres = sort(split.(split(results[2, :Q_RESULTSTRING], " | "), ": "))
    @test Qres[1][1] == "HALClassifier_1"
    @test Qres[2][1] == "InteractionLMClassifier_1"
    @test Qres[3][1] == "XGBoostClassifier_1"
    for i in 1:3
        m_string, std_string = split(Qres[i][2], " ")
        parse(Float64, split(m_string, "=")[2])
        parse(Float64, split(std_string, "=")[2])
    end
    # Check G
    Gres = sort(split.(split(results[2, :G_RESULTSTRING], " | "), ": "))
    @test Gres[1][1] == "LogisticClassifier_1"
    @test Gres[2][1] == "XGBoostClassifier_1"
    for i in 1:2
        m_string, std_string = split(Gres[i][2], " ")
        parse(Float64, split(m_string, "=")[2])
        parse(Float64, split(std_string, "=")[2])
    end

    # Clean
    rm(parsed_args["output"])
    rm(parsed_args["queries"])
end

@testset "Test PhenotypeCrossValidation with continuous target" begin
    # Only one continuous phenotype
    estimatorfile = joinpath("config", "tmle_config.toml")
    build_query_file()
    parsed_args = Dict(
        "phenotypes" => phenotypefile,
        "confounders" => confoundersfile,
        "queries" => queryfile,
        "estimator" => estimatorfile,
        "output" => "cont_results.csv",
        "phenotypes-list" => "data/phen_list_1.csv",
        "verbosity" => 0
    )

    UKBBVariantRun(parsed_args, run_fn=TMLEEpistasis.PhenotypeCrossValidation)

    results = CSV.File(parsed_args["output"]) |> DataFrame

    @test results[1, "VARIANTS"] == "RSID_10 && RSID_100"
    # Check Q
    Qres = sort(split.(split(results[1, :Q_RESULTSTRING], " | "), ": "))
    @test Qres[1][1] == "HALRegressor_1"
    @test Qres[2][1] == "InteractionLMRegressor_1"
    @test Qres[3][1] == "XGBoostRegressor_1"
    @test Qres[4][1] == "XGBoostRegressor_2"
    for i in 1:4
        m_string, std_string = split(Qres[i][2], " ")
        parse(Float64, split(m_string, "=")[2])
        parse(Float64, split(std_string, "=")[2])
    end
    # Check G
    Gres = sort(split.(split(results[1, :G_RESULTSTRING], " | "), ": "))
    @test Gres[1][1] == "LogisticClassifier_1"
    @test Gres[2][1] == "XGBoostClassifier_1"
    for i in 1:2
        m_string, std_string = split(Gres[i][2], " ")
        parse(Float64, split(m_string, "=")[2])
        parse(Float64, split(std_string, "=")[2])
    end

    # Clean
    rm(parsed_args["output"])
    rm(parsed_args["queries"])
end

end;

true


text:
using OffsetArrays
using SpecialFunctions: gamma, erfcx
using PyPlot
using Polynomials
using Printf
using LinearAlgebra

Γ(z) = gamma(z)

function pade(α::T, β::T, m::Integer, n::Integer) where T <: AbstractFloat
    if mod(m+n,2) == 0
        throw(ArgumentError("m+n must be odd"))
    end
    r = div(m+n-1,2)
    a = OffsetArray{T}(undef, 0:m-1)
    b = OffsetArray{T}(undef, 0:n-1)
    pow = one(T)
    for k = 0:m-1
        a[k] = pow / Γ(β+k*α)
        pow = -pow
    end
    b[0] = zero(T)
    pow = one(T)
    for k = 1:n-1
        pow = -pow
        b[k] = pow * sinpi(k*α-β) * Γ(1+k*α-β) / π
    end
    A = zeros(T, 2r-1, 2r-1)
    rhs = zeros(T, 2r-1)
    if m ≥ r+1
        for k = 1:r-1
            A[k,k] = one(T)
            for j = 1:k
                A[k,r-1+j] = -a[k-j]
            end
            rhs[k] = a[k]
        end
        for k = r:m-1
            for j = 1:r
                A[k,r-1+j] = -a[k-j]
            end
            rhs[k] = a[k]
        end
        for k = r-n+1:r-1
            row = k+r
            A[row,k] = one(T)
            for j = k+1:r
                col = j + r - 1
                A[row,col] = -b[j-k]
            end
        end
    else
        for k = 1:m-1
            A[k,k] = one(T)
            for j = 1:k
                col = r - 1 + j
                A[k,col] = -a[k-j]
            end
            rhs[k] = a[k]
        end
        for k = 1:n-r-1
            row = k + m - 1
            for j = 1:r
                col = r - 1 + j
                A[row,col] = -b[j+k]
            end
            rhs[row] = b[k]
        end
        for j = 1:r
            col = r - 1 + j
            A[r,col] = -b[j]
        end
        rhs[r] = -a[0]
        for k = 1:r-1
            row = r + k
            A[row,k] = one(T)
            for j = k+1:r
                col = r - 1 + j
                A[row,col] = -b[j-k]
            end
        end
    end
    coeffs = A \ rhs
    p = OffsetArray{T}(undef, 0:r-1)
    q = OffsetArray{T}(undef, 0:r)
    p[0] = a[0]
    for k = 1:r-1
        p[k] = coeffs[k]
    end
    q[0] = one(T)
    for k = 1:r
        q[k] = coeffs[k+r-1]
    end
    return p, q, A, rhs
end

#T = BigFloat
T = Float64
r = 3
α = 1 / parse(T, "2")
β = parse(T, "1")
x = range(zero(T), parse(T, "10"), length=201)
y = erfcx.(x)
figure(1)
for m = r:5
    n = 2r + 1 - m
    @printf("%3d  %3d", m, n)
    p, q, A, rhs = pade(α, β, m, n)
    @printf("  %5.2f\n", cond(A))
    P = Polynomial(p[0:r-1])
    Q = Polynomial(q[0:r])
    err = P.(x) ./ Q.(x) - y
    plot(x, err, label="m, n = $m, $n")
end
legend()
grid(true)

figure(2)
#for (r, m, n) in [(4, 5, 4), (8, 9, 8), (12, 13, 12), (16, 17, 16)]
for (r, m, n) in [(4, 5, 4), (8, 9, 8), (12, 13, 12)]
    p, q, A, rhs = pade(α, β, m, n)
    println("cond(A) = ", cond(A))
    display(p)
    display(q)
    P = Polynomial(p[0:r-1])
    Q = Polynomial(q[0:r])
    err = P.(x) ./ Q.(x) - y
    semilogy(x, abs.(err), label="m, n = $m, $n")
end
legend()
grid(true)


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

using Markdown
using InteractiveUtils

# ╔═╡ f643a980-5a13-11eb-1016-5f8ebc2b8370
using Plots

# ╔═╡ 13a311d0-5a11-11eb-1377-cd0b5a5ed7de
html"
The Lotka-Volterra Equations
"

# ╔═╡ 55724410-5a1f-11eb-2e9a-813c4e69c2b9
html"
Creating struct LotkaVolterra representing quantities in the equation:
"

# ╔═╡ 4f8f5730-5a11-11eb-1965-0d923a0baf77
@Base.kwdef mutable struct LotkaVolterra
	🐰::Float64	# Prey population
	🐺::Float64	# Predator population
	α::Float64 	# Prey growth factor per year
	β::Float64 	# Prey death factor per year
	δ::Float64	# Predator growth factor per year
	γ::Float64 	# Predator death factor per year
end

# ╔═╡ 839c59c0-5a1f-11eb-24c4-c7a3de982ca9
html"
Time step of 3x a month. The choice is arbitrary.
"

# ╔═╡ 13e4e050-5a12-11eb-1764-cf2089784c84
τ = 1/36

# ╔═╡ a1bc2480-5a1f-11eb-1238-cd6e072da32a
html"
The actual Lotka-Volterra differential equation:
"

# ╔═╡ 2bb5867e-5a12-11eb-135c-cf6728140506
function step!(l::LotkaVolterra)
	l.🐰 += τ*(l.α*l.🐰 - l.β*l.🐰*l.🐺)
	l.🐺 += τ*(l.δ*l.🐰*l.🐺 - l.γ*l.🐺)
end

# ╔═╡ b5717610-5a1f-11eb-03b8-5783009fa91f
html"
Setting things up for plotting:
"

# ╔═╡ 2ac5e370-5a1f-11eb-1966-e3f612893e56
begin
maxsteps = 6000
maxplots = 200
plotstep = maxsteps ÷ maxplots
xplot = zeros(maxplots)
yplot = zeros(maxplots)
tplot = zeros(maxplots)
end

# ╔═╡ c845dd80-5a1f-11eb-0ae1-41ce8fdd7be7
html"
Using real information about baboon/cheetah populations over time, we finally numerically solve the equations:
"

# ╔═╡ b144a830-5a12-11eb-3c23-9338902fef1e
let
lv = LotkaVolterra(🐰=10.0, 🐺=10.0, α=1.1, β=0.4, δ=0.4, γ=0.1)
iplot = 1
for step ∈ 1:maxsteps
	if step%plotstep == 1
		xplot[iplot] = lv.🐰
		yplot[iplot] = lv.🐺
		tplot[iplot] = τ*step
		iplot += 1
	end
	step!(lv)
end
end

# ╔═╡ fdc4e570-5a13-11eb-126a-df201e5b7370
let
plot(tplot, xplot, label="Baboon Population", lw=2)
plot!(tplot, yplot, label="Cheetah Population", lw=2)
plot!(title="Lotka-Volterra Equations", xlabel="Time (years)", ylabel="Population")
end

# ╔═╡ fd1f6030-5a1f-11eb-3c28-7d0a32456864
html"
Now, quickly, a less extreme example:
"

# ╔═╡ bcc7adad-03a4-49a9-9135-286fa3fc61f5
begin
steps = 720
plots = 180
pstep = steps ÷ plots
populations = 0.6:0.1:1.8
🐰 = zeros(length(populations), plots)
🐺 = zeros(length(populations), plots)
t = zeros(plots)
end

# ╔═╡ c19bd9c7-819b-437f-a628-30b2897d9958
let
for (i, pop) ∈ enumerate(populations)
	lk = LotkaVolterra(🐰=pop, 🐺=pop, α=2/3, β=4/3, δ=1, γ=1)
	iplot = 1
	for step ∈ 1:steps
		if step%pstep == 1
			🐰[i, iplot] = lk.🐰
			🐺[i, iplot] = lk.🐺
			t[iplot] = τ*step
			iplot += 1
		end
		step!(lk)
	end
end
end

# ╔═╡ 55b5d660-5a17-11eb-3db6-8394b4373f1b
let
plot(t, 🐰[1, :], label="Prey Population", lw=4)
plot!(t, 🐺[1, :], label="Predator Population", lw=4)
plot!(title="Lotka-Volterra Model",
		xlabel="Time (Years)", ylabel="Population (Thousands)")
end

# ╔═╡ 8950aca0-5a19-11eb-0e00-6f0f90f52a3f
let
plot(title = "Lotka-Volterra Phase Portrait")
xlabel!("Prey Population (Thousands)")
ylabel!("Predator Population (Thousands)")
for (i, pop) ∈ enumerate(populations)
	plot!(🐰[i, 1:100], 🐺[i, 1:100], label="x₀=y₀=$(pop)")
end
xₛ = 1/1 			# point of stability: xₛ = γ/δ = 1/1
yₛ = (2/3)/(4/3) 	# point of stability: yₛ = α/β = (2/3)/(4/3)
scatter!([xₛ], [yₛ], label="Equlibrium ($xₛ, $yₛ)")
end

# ╔═╡ Cell order:
# ╟─13a311d0-5a11-11eb-1377-cd0b5a5ed7de
# ╟─55724410-5a1f-11eb-2e9a-813c4e69c2b9
# ╠═4f8f5730-5a11-11eb-1965-0d923a0baf77
# ╟─839c59c0-5a1f-11eb-24c4-c7a3de982ca9
# ╠═13e4e050-5a12-11eb-1764-cf2089784c84
# ╟─a1bc2480-5a1f-11eb-1238-cd6e072da32a
# ╠═2bb5867e-5a12-11eb-135c-cf6728140506
# ╟─b5717610-5a1f-11eb-03b8-5783009fa91f
# ╠═2ac5e370-5a1f-11eb-1966-e3f612893e56
# ╟─c845dd80-5a1f-11eb-0ae1-41ce8fdd7be7
# ╠═b144a830-5a12-11eb-3c23-9338902fef1e
# ╠═f643a980-5a13-11eb-1016-5f8ebc2b8370
# ╠═fdc4e570-5a13-11eb-126a-df201e5b7370
# ╟─fd1f6030-5a1f-11eb-3c28-7d0a32456864
# ╠═bcc7adad-03a4-49a9-9135-286fa3fc61f5
# ╠═c19bd9c7-819b-437f-a628-30b2897d9958
# ╠═55b5d660-5a17-11eb-3db6-8394b4373f1b
# ╠═8950aca0-5a19-11eb-0e00-6f0f90f52a3f

max_stars_repo_name: benharmonics/GeneralPhysics

max_issues_repo_name: benharmonics/GeneralPhysics

max_forks_repo_name: benharmonics/GeneralPhysics


text:
##
using DelimitedFiles, Logging, Pkg, Test
cd(@__DIR__())
pkg"activate ./Intcode"
using Intcode
logfile = "intcode.log"
function do_with_logger(f::Function, level = Logging.Info, logfile=logfile)
    open(logfile, "w+") do io
        with_logger(ConsoleLogger(io, level)) do
            f(io)
        end
        flush(io)
    end
end

##
try
    pkg"test Intcode"
catch e
    @error e
end

##
init_data = readdlm("input.txt", ',', Int) |> vec;

##
do_with_logger() do io
    ic = IntcodeMachine(init_data)
    run_intcode!(ic)
    res = fetch(ic)
    isnothing(res) && @error "TEST program failed"
    @info something(res)
    !isnothing(io) && flush(io)
end

##
include(joinpath("Intcode", "test", "test_cmp_jmp.jl"))

##
do_with_logger(test_cmp, Logging.Debug)

##
do_with_logger(test_jmp, Logging.Debug)

##
do_with_logger(test_larger, Logging.Debug)


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

@testset "get travel mode" begin
	sim = initSim("data/cities/small/1/sim_config.xml") # see travel.csv file in same folder
	travel = sim.travel # shorthand

	# first test getTravelSetsStartTimesIndex(), which getTravelMode!() relies on
	f(t::Float) = JEMSS.getTravelSetsStartTimesIndex(travel, t)
	result = true
	travelSetsStartTimes = [0, 0.2, 0.5, 0.8]
	for t = 0:0.01:1.0
		expectedValue = findlast(t .>= travelSetsStartTimes)
		result &= f(t) == expectedValue
	end
	@test result

	# test getTravelMode!()
	result = true
	f(p::Priority, t::Float) = JEMSS.getTravelMode!(travel, p, sim.time; startTime = t).index
	expectedValues = Dict(
		(highPriority, 0.0) => 1,
		(medPriority, 0.0) => 2,
		(lowPriority, 0.0) => 3,
		(highPriority, 1.0) => 2,
		(medPriority, 1.0) => 3,
		(lowPriority, 1.0) => 4,
	)
	for priority in priorities, t = (0.0, 1.0)
		result &= f(priority, t) == expectedValues[(priority, t)]
	end
	@test result

	# check that getTravelMode!() updates travel.recentSetsStartTimesIndex
	function f(t::Float)
		JEMSS.getTravelMode!(travel, highPriority, t)
		return travel.recentSetsStartTimesIndex
	end
	result = true
	travelSetsStartTimes = [0, 0.2, 0.5, 0.8]
	for t = 0:0.01:1.0
		expectedValue = findlast(t .>= travelSetsStartTimes)
		result &= f(t) == expectedValue
	end
	@test result
end


text:
# slt_header.jl
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 0.7
# 09.12.2016
# Last Edit: 18.04.2019

"""# slt_header

Load an SLT header

`slt_header(f,AnalogIn,freq)`\n
* **f**::String = File name
* **AnalogIn**::Int = Number of analog signals
* **freq**::Number or String = Sample frequency"""
function slt_header(f::String,AnalogIn::Int,freq::String)
	return slt_header(f,AnalogIn,parse(Float64,freq))
end # slt_header(f::String,AnalogIn::Int,freq::String)

function slt_header(f::String,AnalogIn::Int,freq::Number)
	# Prepare Output Variables
	bpr = Int # Bytes Per Record
	eddymeasver = [] # Eddy Meas Version
	t0 = DateTime # Initial Timestamp
	ch = [] # Channels
	bm = [] # Bit Masks
	l = [] # Number of Rows of Data
	stpos = [] # Starting Position of Data

	fid = open(f,"r")
	try
		# Header
		bpr = Int(read!(fid,Array{Int8}(undef,1))[1]) # Bytes Per Record
		eddymeasver = read!(fid,Array{Int8}(undef,1))
		dom = Int(read!(fid,Array{Int8}(undef,1))[1]) # Day of Month
		m = Int(read!(fid,Array{Int8}(undef,1))[1]) # Month number
		yr = 100*Int(read!(fid,Array{Int8}(undef,1))[1]) + Int(read!(fid,Array{Int8}(undef,1))[1]) # Year
		h = Int(read!(fid,Array{Int8}(undef,1))[1]) # Hour
		minut = Int(read!(fid,Array{Int8}(undef,1))[1]) # Minute
		t0 = DateTime(yr,m,dom,h,minut,0) # Initial file time

		fs = stat(f).size # File size in bytes
		l = (fs - 8 - 2*AnalogIn)/bpr # Number of rows of data

		# Bit Masks and Channels
		bm = read!(fid,Array{Int8}(undef,2*AnalogIn))
		ch = bm[2:2:end] # Channels
		bm = bm[1:2:end] # Bit masks

		stpos = position(fid) # Position of the first data point
	catch
		println("Error loading file: " * f)
	end
	close(fid)

	# Format Output
	output = DataFrame(FileName = f,T0 = t0,Line_Count = l,Start_Pos = stpos,Channels = collect(Array[ch]),Bit_Mask = collect(Array[bm]),EddyMeas_Version = eddymeasver,BytesPerRecord = bpr)

	return output
end # slt_header(f::String,AnalogIn::Int,freq::Number)

max_stars_repo_name: JeremyRueffer/ClimateDataIO

max_issues_repo_name: JeremyRueffer/ClimateDataIO


text:
typealias Weight Matrix{Float32}
typealias Bias   Vector{Float32}
typealias Signal Matrix{Float32}

abstract Optimizer
abstract Layer
abstract Filter
abstract RecurrentLayer <: Layer
abstract ActivationFunc

immutable HiddenLayer <: Layer
    batch_size::Int64
    W::Weight
    ΔW::Weight
    b::Bias
    Δb::Bias
    ΔE::Signal
    X::Signal
    U::Signal
    Z::Signal
    Zd::Signal
    activation::ActivationFunc
    post_filters::Vector{Filter}

    function HiddenLayer(batch_size, x_size, h_size; activation=LeakyRelu(), post_filters=Filter[])
        w_in    = randn(h_size, x_size)  * 0.1
        dw_in   = zeros(h_size, x_size)
        b       = zeros(h_size)
        db      = zeros(h_size)
        de      = zeros(h_size, batch_size)
        x       = zeros(x_size, batch_size)
        u       = zeros(h_size, batch_size)
        z       = zeros(h_size, batch_size)
        dz      = zeros(h_size, batch_size)
        new(batch_size, w_in, dw_in, b, db, de, x, u, z, dz, activation, post_filters)
    end
end

immutable OutputLayer <: Layer
    batch_size::Int64
    W::Weight
    ΔW::Weight
    b::Bias
    Δb::Bias
    ΔE::Signal
    X::Signal
    U::Signal
    Z::Signal
    Zd::Signal
    activation::ActivationFunc
    loss::Function
    post_filters::Vector{Filter}

    function OutputLayer(batch_size, x_size, y_size; activation=IdentityAct(), loss=quadratic_cost, post_filters=Filter[])
        w_in    = randn(y_size, x_size)  * 0.05
        dw_in   = zeros(y_size, x_size)
        b       = zeros(y_size)
        db      = zeros(y_size)
        de      = zeros(y_size, batch_size)
        x       = zeros(x_size, batch_size)
        u       = zeros(y_size, batch_size)
        z       = zeros(y_size, batch_size)
        dz      = zeros(y_size, batch_size)
        new(batch_size, w_in, dw_in, b, db, de, x, u, z, dz, activation, loss, post_filters)
    end
end

immutable SimpleRecurrentLayer <: RecurrentLayer
    batch_size::Int64
    W::Weight
    ΔW::Weight
    Wh::Weight
    ΔWh::Weight
    b::Bias
    Δb::Bias
    ΔE::Signal
    X::Signal
    U::Signal
    Uh::Signal
    Z::Signal
    Zd::Signal
    activation::ActivationFunc
    post_filters::Vector{Filter}

    function SimpleRecurrentLayer(batch_size, x_size, h_size; activation=LeakyRelu(), post_filters=Filter[])
        w_in  = randn(h_size, x_size)  * 0.1
        dw_in = zeros(h_size, x_size)
        w_h   = eye(h_size) * 0.98
        dw_h  = zeros(h_size, h_size)
        b     = zeros(h_size)
        db    = zeros(h_size)
        ΔE   = zeros(h_size, batch_size)
        X     = zeros(x_size, batch_size)
        U     = zeros(h_size, batch_size)
        Uh    = zeros(h_size, batch_size)
        Z     = zeros(h_size, batch_size)
        Zd    = zeros(h_size, batch_size)
        layer = new(batch_size, w_in, dw_in, w_h, dw_h, b, db, ΔE, X, U, Uh, Z, Zd, activation, post_filters)
        reset!(layer, zeros(Float32, x_size, batch_size); force=true)
        layer
    end
end

immutable Network
    batch_size::Int64
    layers::Array{Layer}
    optimizers::Array{Optimizer}
    rnn_optimizers::Array{Optimizer}

    function Network(batch_size)
        new(batch_size, Layer[], Optimizer[], Optimizer[])
    end
end

function push!(net::Network, layer::Layer, optimizer::Optimizer)
    Base.push!(net.layers,         layer)
    Base.push!(net.optimizers,     optimizer)
    Base.push!(net.rnn_optimizers, optimizer)
    true
end

function push!(net::Network, layer::Layer, optimizer::Optimizer, rnn_optimizer::Optimizer)
    Base.push!(net.layers,         layer)
    Base.push!(net.optimizers,     optimizer)
    Base.push!(net.rnn_optimizers, rnn_optimizer)
    true
end

function reset!(l::Layer, nextX::Signal)
    # nothing to do
end

function reset!(l::SimpleRecurrentLayer, nextX::Signal; force=false)
    x_size     = size(l.W, 2)
    h_size     = size(l.W, 1)
    batch_size = size(nextX, 2)

    if force || batch_size != size(l.X,2)
        copy!(l.ΔE, zeros(h_size, batch_size))
        copy!(l.X  , zeros(x_size, batch_size))
        copy!(l.U  , zeros(h_size, batch_size))
        copy!(l.Uh , zeros(h_size, batch_size))
        copy!(l.Z  , zeros(h_size, batch_size))
        copy!(l.Zd , zeros(h_size, batch_size))
    end
end


text:
export addbmm, addbmm!, bmm!, bmm

function bmm!(tA::Char, tB::Char, result::Array{T, 2}, alpha::T, batch1::Array{T, 3}, batch2::Array{T, 3}) where T
    @assert size(batch1, 3) == size(batch2, 3) "batch size does not match"

    batch = size(batch1, 3)
    for i = 1:batch
        @inbounds vbatch1 = view(batch1, :, :, i)
        @inbounds vbatch2 = view(batch2, :, :, i)
        BLAS.gemm!(tA, tB, alpha, vbatch1, vbatch2, one(T), result)
    end
    return result
end

bmm!(tA::Char, tB::Char, result::Array{T, 2}, batch1::Array{T, 3}, batch2::Array{T, 3}) where T =
    bmm!(tA, tB, result, one(T), batch1, batch2)
bmm(tA::Char, tB::Char, alpha::T, batch1::Array{T, 3}, batch2::Array{T, 3}) where T =
    bmm!(tA, tB, zeros(size(batch1, 1), size(batch2, 2)), alpha, batch1, batch2)
bmm(tA::Char, tB::Char, batch1::Array{T, 3}, batch2::Array{T, 3}) where T =
    bmm(tA, tB, one(T), batch1, batch2)

function addbmm!(tA::Char, tB::Char, result::Array{T, 2}, beta::T, mat::Array{T, 2}, alpha::T, batch1::Array{T, 3}, batch2::Array{T, 3}) where T
    bmm!(tA, tB, result, alpha, batch1, batch2)
    result += beta * mat
    return result
end

addbmm!(tA::Char, tB::Char, result::Array{T, 2}, mat::Array{T, 2}, batch1::Array{T, 3}, batch2::Array{T, 3}) where T =
    addbmm!(tA, tB, result, one(T), mat, one(T), batch1, batch2)

addbmm(tA::Char, tB::Char, beta::T, mat::Array{T, 2}, alpha::T, batch1::Array{T, 3}, batch2::Array{T, 3}) where T =
    addbmm!(tA, tB, zeros(size(batch1, 1), size(batch2, 2)), beta, mat, alpha, batch1, batch2)
addbmm(tA::Char, tB::Char, mat::Array{T, 2}, batch1::Array{T, 3}, batch2::Array{T, 3}) where T =
    addbmm!(tA, tB, one(T), mat, one(T), batch1, batch2)


text:
"""
Adds pose nodes to graph with a prior on Z, pitch, and roll.
"""
function handle_poses!(slam::SyncrSLAM,
                       msg::pose_node_t)
    id = msg.id
    println("[Caesar.jl] Received pose msg for x$(id)")

    mean = msg.mean
    covar = msg.covar
    t = [mean[1], mean[2], mean[3]]
    qw = mean[4]
    qxyz = [mean[5], mean[6], mean[7]]
    q = Quaternion(qw,qxyz) # why not a (w,x,y,z) constructor?
    pose = SE3(t,q)
    euler = Euler(q)

    node_label = Symbol("x$(id)")
    varRequest = VariableRequest(node_label, "Pose3", nothing, ["POSE"])
    resp = addVariable(slam.syncrconf, slam.robotId, slam.sessionId, varRequest)

    if id == 0
        println("[Caesar.jl] First pose")
        # this is the first msg, and it does not carry odometry, but the prior on the first node.

        # add 6dof prior
        initPosePrior = RoME.PriorPose3( MvNormal( veeEuler(pose), diagm([covar...]) ) )
        packedPrior = convert(RoME.PackedPriorPose3, initPosePrior)

        # 3. Build the factor request (again, we can make this way easier and transparent once it's stable)
        fctBody = FactorBody(string(typeof(initPosePrior)), string(typeof(packedPrior)), "JSON", JSON.json(packedPrior))
        fctRequest = FactorRequest([node_label], fctBody, false, false)
        resp = addFactor(slam.syncrconf, slam.robotId, slam.sessionId, fctRequest)

        # set robot parameters in the first pose, this will become a separate node in the future
        println("[Caesar.jl] Setting robot parameters")
        setRobotParameters!(slam)
    end
end

"""
Handle ZPR priors on poses.
"""
function handle_priors!(slam::SyncrSLAM,
                         msg::prior_zpr_t)

    id = msg.id
    println("[Caesar.jl] Adding prior on RPZ to x$(id)")

    node_label = Symbol("x$(id)")
    # 1. Build up the prior
    z = msg.z
    pitch = msg.pitch
    roll = msg.roll
    var_z = msg.var_z
    var_pitch = msg.var_pitch
    var_roll = msg.var_roll
    rp_dist = MvNormal( [roll;pitch], diagm([var_roll, var_pitch]))
    z_dist = Normal(z, var_z)
    prior_rpz = RoME.PartialPriorRollPitchZ(rp_dist, z_dist)

    # 2. Pack the prior (we can automate this step soon, but for now it's hand cranking)
    packed_prior_rpz = convert(RoME.PackedPartialPriorRollPitchZ, prior_rpz)

    # 3. Build the factor request (again, we can make this way easier and transparent once it's stable)
    fctBody = FactorBody(string(typeof(prior_rpz)), string(typeof(packed_prior_rpz)), "JSON", JSON.json(packed_prior_rpz))
    fctRequest = FactorRequest([node_label], fctBody, false, false)
    resp = addFactor(slam.syncrconf, slam.robotId, slam.sessionId, fctRequest)
end

"""
Handle partial x, y, and heading odometry constraints between Pose3 variables.
"""
function handle_partials!(slam::SyncrSLAM,
                         msg::pose_pose_nh_t)
    println(" --- Handling odometry change...")
    @show msg

    origin_id = msg.node_1_id
    destination_id = msg.node_2_id
    origin_label = Symbol("x$(origin_id)")
    destination_label = Symbol("x$(destination_id)")

    println("[Caesar.jl] Adding XYH odometry constraint between(x$(origin_id), x$(destination_id))")

    mean = msg.mean
    covar = msg.covar
    t = [mean[1], mean[2], mean[3]]
    qw = mean[4]
    qxyz = [mean[5], mean[6], mean[7]]
    q = Quaternion(qw,qxyz) # why not a (w,x,y,z) constructor?
    pose = SE3(t,q)
    euler = Euler(q)
    yaw = euler.Y

    xyh_dist = MvNormal([t[1], t[2], yaw], diagm([covar[1], covar[2], covar[3]]))
    xyh_factor = RoME.PartialPose3XYYaw(xyh_dist)
    # 2. Pack the prior (we can automate this step soon, but for now it's hand cranking)
    packed_xyh_factor = convert(RoME.PackedPartialPose3XYYaw, xyh_factor)

    # 3. Build the factor request (again, we can make this way easier and transparent once it's stable)
    fctBody = FactorBody(string(typeof(xyh_factor)), string(typeof(packed_xyh_factor)), "JSON", JSON.json(packed_xyh_factor))
    fctRequest = FactorRequest([origin_label; destination_label], fctBody, false, false)
    resp = addFactor(slam.syncrconf, slam.robotId, slam.sessionId, fctRequest)
end

"""
Handle loop closure proposals with chance of being a null hypothesis likelihood.
"""
function handle_loops!(slam::SyncrSLAM,
                       msg::pose_pose_xyh_nh_t)

    println(" --- Handling loop...")
    # return nothing

    origin_id = msg.node_1_id
    destination_id = msg.node_2_id
    origin_label = Symbol("x$(origin_id)")
    destination_label = Symbol("x$(destination_id)")

    delta_x = msg.delta_x
    delta_y = msg.delta_y
    delta_yaw = msg.delta_yaw

    var_x = msg.var_x
    var_y = msg.var_y
    var_yaw = msg.var_yaw
    confidence = msg.confidence

    println("[Caesar.jl] Adding XYH-NH loop closure constraint between (x$(origin_id), x$(destination_id))")
    xyh_dist = MvNormal([delta_x, delta_y, delta_yaw], diagm([var_x, var_y, var_yaw]))
    xyh_factor = RoME.PartialPose3XYYawNH(xyh_dist ,[1.0-confidence, confidence]) # change to NH
    # 2. Pack the prior (we can automate this step soon, but for now it's hand cranking)
    packed_xyh_factor = convert(RoME.PackedPartialPose3XYYawNH, xyh_factor)

    # 3. Build the factor request (again, we can make this way easier and transparent once it's stable)
    fctBody = FactorBody(string(typeof(xyh_factor)), string(typeof(packed_xyh_factor)), "JSON", JSON.json(packed_xyh_factor))
    fctRequest = FactorRequest([origin_label; destination_label], fctBody, false, false)
    resp = addFactor(slam.syncrconf, slam.robotId, slam.sessionId, fctRequest)
end

"""
Callback for caesar_point_cloud_t msgs. Adds point cloud to SLAM_Client
"""
function handle_clouds!(slam::SyncrSLAM,
                        msg::point_cloud_t)
    id = msg.id
    last_pose = "x$(id)"
    println("[Caesar.jl] Got cloud $id")

    points = JSON.json(msg.points)
    colors = JSON.json(msg.colors)

    node = getNode(slam.syncrconf, slam.robotId, slam.sessionId, last_pose)
    pointCloud = BigDataElementRequest("PointCloud", "Mongo", "Pointcloud from HAUV", points, "application/binary")
    colors = BigDataElementRequest("Colors", "Mongo", "Colors from HAUV", colors, "application/binary")
    addDataElement(slam.syncrconf, slam.robotId, slam.sessionId, node.id, pointCloud)
    addDataElement(slam.syncrconf, slam.robotId, slam.sessionId, node.id, colors)
    # serialized_point_cloud = BSONObject(Dict("pointcloud" => points))
    # appendvertbigdata!(slaml, vert, "BSONpointcloud", string(serialized_point_cloud).data)
    # serialized_colors = BSONObject(Dict("colors" => colors))
    # appendvertbigdata!(slaml, vert, "BSONcolors", string(serialized_colors).data)
end


text:
#jl #! format: off
# # Large Scale Production Cost Modeling with [PowerSimulations.jl](https://github.com/NREL-SIIP/PowerSimulations.jl)

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# This example shows a basic PCM simulation using the system data assembled in the
# [US-System example](https://nbviewer.jupyter.org/github/NREL-SIIP/SIIPExamples.jl/blob/master/notebook/2_PowerSystems_examples/08_US_system.ipynb).

# ### Dependencies
using SIIPExamples
using PowerSystems
using PowerSimulations
using PowerGraphics
using Logging
using Dates

pkgpath = pkgdir(SIIPExamples)
PSI = PowerSimulations
plotlyjs()

# ### Optimization packages
# You can use the free solvers as in one of the other PowerSimulations.jl examples, but on
# large problems it's useful to have a solver with better performance. I'll use the Xpress
# solver since I have a license, but others such as Gurobi or CPLEX work well too.
using Xpress
solver = optimizer_with_attributes(Xpress.Optimizer, "MIPRELSTOP" => 0.1, "OUTPUTLOG" => 1)

# ### Load the US `System`.
# If you have run the
# [US-System example](https://nbviewer.jupyter.org/github/NREL-SIIP/SIIPExamples.jl/blob/master/notebook/2_PowerSystems_examples/08_US-System.ipynb), the data will
# be serialized in the json and H5 format, so we can efficiently deserialize it:

sys = System(joinpath(pkgpath, "US-System", "SIIP", "sys.json"))
horizon = 24;
interval = Dates.Hour(24);
transform_single_time_series!(sys, horizon, interval);

# ### Selecting flow limited lines
# Since PowerSimulations will apply constraints by component type (e.g. Line), we need to
# change the component type of the lines on which we want to enforce flow limits. So, let's
# change the device type of certain branches from Line to MonitoredLine differentiate
# treatment when we build the model. Here, we can select inter-regional lines, or lines
# above a voltage threshold.

for line in get_components(Line, sys)
    if (get_base_voltage(get_from(get_arc(line))) >= 230.0) &&
       (get_base_voltage(get_to(get_arc(line))) >= 230.0)
        #if get_area(get_from(get_arc(line))) != get_area(get_to(get_arc(line)))
        @info "Changing $(get_name(line)) to MonitoredLine"
        convert_component!(MonitoredLine, line, sys)
    end
end

# ### Create a `template`
# Now we can create a `template` that applies an unbounded formulation to `Line`s and the standard
# flow limited formulation to `MonitoredLine`s.
template = ProblemTemplate(DCPPowerModel)
set_device_model!(template, Line, StaticBranchUnbounded)
set_device_model!(template, TapTransformer, StaticBranchUnbounded)
set_device_model!(template, MonitoredLine, StaticBranch)
set_device_model!(template, ThermalStandard, ThermalStandardUnitCommitment)
set_device_model!(template, RenewableDispatch, RenewableFullDispatch)
set_device_model!(template, PowerLoad, StaticPowerLoad)
set_device_model!(template, HydroDispatch, FixedOutput)

# ### Build and execute single step problem
op_problem = DecisionModel(template, sys; optimizer = solver, horizon = 24)

build!(op_problem, output_dir = mktempdir(), console_level = Logging.Info)

solve!(op_problem)

# ### Analyze results
res = ProblemResults(op_problem)
plot_fuel(res)

# ## Sequential Simulation
# In addition to defining the formulation template, sequential simulations require
# definitions for how information flows between problems.
sim_folder = mkpath(joinpath(pkgpath, "Texas-sim"))
models = SimulationModels(
    decision_models = [
        DecisionModel(
            template,
            sys,
            name = "UC",
            optimizer = solver,
            system_to_file = false,
        ),
    ],
)
DA_sequence =
    SimulationSequence(models = models, ini_cond_chronology = IntraProblemChronology())

# ### Define and build a simulation
sim = Simulation(
    name = "Texas-test",
    steps = 3,
    models = models,
    sequence = DA_sequence,
    simulation_folder = "Texas-sim",
)

build!(
    sim,
    console_level = Logging.Info,
    file_level = Logging.Debug,
    recorders = [:simulation],
)

# ### Execute the simulation
#nb execute!(sim)

# ### Load and analyze results
#nb results = SimulationResults(sim);
#nb uc_results = get_decision_problem_results(results, "UC");

#nb plot_fuel(uc_results)

max_stars_repo_name: NREL-SIIP/Examples

max_issues_repo_name: NREL-SIIP/Examples

max_forks_repo_name: NREL-SIIP/Examples


text:
# Julia wrapper for header: /usr/include/expat.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@c None XML_SetElementDeclHandler (XML_Parser, XML_ElementDeclHandler) libexpat
@c None XML_SetAttlistDeclHandler (XML_Parser, XML_AttlistDeclHandler) libexpat
@c None XML_SetXmlDeclHandler (XML_Parser, XML_XmlDeclHandler) libexpat
@c XML_Parser XML_ParserCreate (Ptr{XML_Char},) libexpat
@c XML_Parser XML_ParserCreateNS (Ptr{XML_Char}, XML_Char) libexpat
@c XML_Parser XML_ParserCreate_MM (Ptr{XML_Char}, Ptr{XML_Memory_Handling_Suite}, Ptr{XML_Char}) libexpat
@c XML_Bool XML_ParserReset (XML_Parser, Ptr{XML_Char}) libexpat
@c None XML_SetEntityDeclHandler (XML_Parser, XML_EntityDeclHandler) libexpat
@c None XML_SetElementHandler (XML_Parser, XML_StartElementHandler, XML_EndElementHandler) libexpat
@c None XML_SetStartElementHandler (XML_Parser, XML_StartElementHandler) libexpat
@c None XML_SetEndElementHandler (XML_Parser, XML_EndElementHandler) libexpat
@c None XML_SetCharacterDataHandler (XML_Parser, XML_CharacterDataHandler) libexpat
@c None XML_SetProcessingInstructionHandler (XML_Parser, XML_ProcessingInstructionHandler) libexpat
@c None XML_SetCommentHandler (XML_Parser, XML_CommentHandler) libexpat
@c None XML_SetCdataSectionHandler (XML_Parser, XML_StartCdataSectionHandler, XML_EndCdataSectionHandler) libexpat
@c None XML_SetStartCdataSectionHandler (XML_Parser, XML_StartCdataSectionHandler) libexpat
@c None XML_SetEndCdataSectionHandler (XML_Parser, XML_EndCdataSectionHandler) libexpat
@c None XML_SetDefaultHandler (XML_Parser, XML_DefaultHandler) libexpat
@c None XML_SetDefaultHandlerExpand (XML_Parser, XML_DefaultHandler) libexpat
@c None XML_SetDoctypeDeclHandler (XML_Parser, XML_StartDoctypeDeclHandler, XML_EndDoctypeDeclHandler) libexpat
@c None XML_SetStartDoctypeDeclHandler (XML_Parser, XML_StartDoctypeDeclHandler) libexpat
@c None XML_SetEndDoctypeDeclHandler (XML_Parser, XML_EndDoctypeDeclHandler) libexpat
@c None XML_SetUnparsedEntityDeclHandler (XML_Parser, XML_UnparsedEntityDeclHandler) libexpat
@c None XML_SetNotationDeclHandler (XML_Parser, XML_NotationDeclHandler) libexpat
@c None XML_SetNamespaceDeclHandler (XML_Parser, XML_StartNamespaceDeclHandler, XML_EndNamespaceDeclHandler) libexpat
@c None XML_SetStartNamespaceDeclHandler (XML_Parser, XML_StartNamespaceDeclHandler) libexpat
@c None XML_SetEndNamespaceDeclHandler (XML_Parser, XML_EndNamespaceDeclHandler) libexpat
@c None XML_SetNotStandaloneHandler (XML_Parser, XML_NotStandaloneHandler) libexpat
@c None XML_SetExternalEntityRefHandler (XML_Parser, XML_ExternalEntityRefHandler) libexpat
@c None XML_SetExternalEntityRefHandlerArg (XML_Parser, Ptr{Void}) libexpat
@c None XML_SetSkippedEntityHandler (XML_Parser, XML_SkippedEntityHandler) libexpat
@c None XML_SetUnknownEncodingHandler (XML_Parser, XML_UnknownEncodingHandler, Ptr{Void}) libexpat
@c None XML_DefaultCurrent (XML_Parser,) libexpat
@c None XML_SetReturnNSTriplet (XML_Parser, Int32) libexpat
@c None XML_SetUserData (XML_Parser, Ptr{Void}) libexpat
@c Int32 XML_SetEncoding (XML_Parser, Ptr{XML_Char}) libexpat
@c None XML_UseParserAsHandlerArg (XML_Parser,) libexpat
@c Int32 XML_UseForeignDTD (XML_Parser, XML_Bool) libexpat
@c Int32 XML_SetBase (XML_Parser, Ptr{XML_Char}) libexpat
@c Ptr{XML_Char} XML_GetBase (XML_Parser,) libexpat
@c Int32 XML_GetSpecifiedAttributeCount (XML_Parser,) libexpat
@c Int32 XML_GetIdAttributeIndex (XML_Parser,) libexpat
@c Int32 XML_Parse (XML_Parser, Ptr{Uint8}, Int32, Int32) libexpat
@c Ptr{Void} XML_GetBuffer (XML_Parser, Int32) libexpat
@c Int32 XML_ParseBuffer (XML_Parser, Int32, Int32) libexpat
@c Int32 XML_StopParser (XML_Parser, XML_Bool) libexpat
@c Int32 XML_ResumeParser (XML_Parser,) libexpat
@c None XML_GetParsingStatus (XML_Parser, Ptr{XML_ParsingStatus}) libexpat
@c XML_Parser XML_ExternalEntityParserCreate (XML_Parser, Ptr{XML_Char}, Ptr{XML_Char}) libexpat
@c Int32 XML_SetParamEntityParsing (XML_Parser, Void) libexpat
@c Int32 XML_SetHashSalt (XML_Parser, Uint32) libexpat
@c Int32 XML_GetErrorCode (XML_Parser,) libexpat
@c XML_Size XML_GetCurrentLineNumber (XML_Parser,) libexpat
@c XML_Size XML_GetCurrentColumnNumber (XML_Parser,) libexpat
@c XML_Index XML_GetCurrentByteIndex (XML_Parser,) libexpat
@c Int32 XML_GetCurrentByteCount (XML_Parser,) libexpat
@c Ptr{Uint8} XML_GetInputContext (XML_Parser, Ptr{Int32}, Ptr{Int32}) libexpat
@c None XML_FreeContentModel (XML_Parser, Ptr{XML_Content}) libexpat
@c Ptr{Void} XML_MemMalloc (XML_Parser, size_t) libexpat
@c Ptr{Void} XML_MemRealloc (XML_Parser, Ptr{Void}, size_t) libexpat
@c None XML_MemFree (XML_Parser, Ptr{Void}) libexpat
@c None XML_ParserFree (XML_Parser,) libexpat
@c Ptr{XML_LChar} XML_ErrorString (Void,) libexpat
@c Ptr{XML_LChar} XML_ExpatVersion () libexpat
@c XML_Expat_Version XML_ExpatVersionInfo () libexpat
@c Ptr{XML_Feature} XML_GetFeatureList () libexpat


text:
### All possible fields
### If it initialisation lacks some fields then it might be resolved with simpler constructors
### Similar type could be made for RestructurisationResults
struct StabilisationResults
    zc
    vinit
    Finit
    vres
    Fres
end

function Base.show(io::IO, r::StabilisationResults)
    println("Results of Stabilisation Algorithm $(typeof(r.zc))")
    println(" * vinit: $(r.vinit)")
    println(" * Finit: $(r.Finit)")
    println(" * vres: $(r.vres)")
    println(" * Fres: $(r.Fres)")
end

# High level interface. If you trust enough ;)
function stabilise(points,faces,v;method=:Zinchenko2013,vp=zeros(size(points)...))

    n = Array{Float64}(undef,size(points)...)
    NormalVectors!(n,points,faces,i -> FaceVRing(i,faces))

    ### Projecting previous velocity
    for xkey in 1:size(points,2)
        P = eye(3) - n[:,xkey]*n[:,xkey]'
        v[:,xkey] = n[:,xkey]*dot(n[:,xkey],v[:,xkey]) + P*vp[:,xkey]
    end

    if method==:Zinchenko2013
        zc = Zinchenko2013(points,faces,n)
    elseif method==:Zinchenko1997
        zc = Zinchenko1997(i -> VertexVRing(i,faces))
    elseif method==:Erdmanis2016
        zc = Erdmanis2016()
    else
        zc = method   ### For construction outside the usual
    end

    vres = copy(v)
    Finit = F(points,faces,v,zc)
    stabilise!(points,faces,n,vres,zc)
    Fres = F(points,faces,vres,zc)
    res = StabilisationResults(zc,v,Finit,vres,Fres)
    return res
end

# function stabilise(points,faces,v;method=:Zinchenko2013)

#     n = Array(Float64,size(points)...)
#     for i in 1:size(points,2)
#         iter = FaceVRing(i,faces)
#         n[:,i] = vnormal(i,points,faces,iter)
#     end

#     stabilise(points,faces,n,v;method=method)
# end