text:
"""
This file contains necessary utility functions required for column generation algorithm.
"""
"""
our problem or model structure is as follows:
min z[1]+z[2]+z[3]
st.
    Ax+By <= 0

    z[1] >= y[1]*L[1]
    z[2] >= y[2]*L[2]
    z[3] >= y[3]*L[3]

    z[1] >= (y[1]-1)*U[1]+∇f1(̄x)(x-̄x)
    z[2] >= (y[2]-1)*U[2]+∇f2(̄x)(x-̄x)
    z[3] >= (y[3]-1)*U[3]+∇f3(̄x)(x-̄x)

    x = [x[1], x[2], x[3]], :Cont
    y = [y[1], y[2], y[3]], :Bin
    z = [z[1], z[2], z[3]], :Cont
    A = (mxn) matrix
    B = (mxn) matrix

First, we need to do dantzig-wolfe decomposition of the above problem .

The constraint set
"Ax+By <= 0"
acts as feasibility constraint and these constraints indirectly influence the objective function.

The constraint set
"z[1] >= y[1]*L[1]
 z[2] >= y[2]*L[2]
 z[3] >= y[3]*L[3]
 z[1] >= (y[1]-1)*U[1]+∇f1(̄x)(x-̄x)
 z[2] >= (y[2]-1)*U[2]+∇f2(̄x)(x-̄x)
 z[3] >= (y[3]-1)*U[3]+∇f3(̄x)(x-̄x)"
directly influence the objective function.

Now we can do Dantzig-Wolfe decomposition of the model as follows:

main problem:
min z[1]+z[2]+z[3]
(x,y) ∈ Q
st.
    z[1] >= y[1]*L[1]
    z[2] >= y[2]*L[2]
    z[3] >= y[3]*L[3]

    z[1] >= (y[1]-1)*U[1]+∇f1(̄x)(x-̄x)
    z[2] >= (y[2]-1)*U[2]+∇f2(̄x)(x-̄x)
    z[3] >= (y[3]-1)*U[3]+∇f3(̄x)(x-̄x)

where Q = {(x,y): Ax+By <= 0, x :Cont, y :Cont}.

The Dantzig-Wolfe relaxation of the main problem:

min z[1]+z[2]+z[3]
(λ,μ)
st.
    z[1] >= sum(u_p[1,i]*L[1]*λ[i] for i in Κ) + sum(u_d[1,j]*L[1]*μ[j] for j in Τ)
    z[2] >= sum(u_p[2,i]*L[2]*λ[i] for i in Κ) + sum(u_d[2,j]*L[2]*μ[j] for j in Τ)
    z[3] >= sum(u_p[3,i]*L[3]*λ[i] for i in Κ) + sum(u_d[3,j]*L[3]*μ[j] for j in Τ)

    z[1]+∇f1(̄x)*̄x+U[1] >= sum((u_p[1,i]*U[1]+∇f1(̄x)*v_p[i])*λ[i] for i in Κ)+sum((u_d[1,j]*U[1]+∇f1(̄x)*v_d[j])*μ[j] for j in Τ)
    z[2]+∇f2(̄x)*̄x+U[2] >= sum((u_p[2,i]*U[2]+∇f2(̄x)*v_p[i])*λ[i] for i in Κ)+sum((u_d[2,j]*U[2]+∇f2(̄x)*v_d[j])*μ[j] for j in Τ)
    z[3]+∇f3(̄x)*̄x+U[3] >= sum((u_p[3,i]*U[3]+∇f3(̄x)*v_p[i])*λ[i] for i in Κ)+sum((u_d[3,j]*U[3]+∇f3(̄x)*v_d[j])*μ[j] for j in Τ)

    sum(λ[i] for i in K) = 1
    λ[i] >= 0, ∀i ∈ K
    μ[j] >= 0, ∀j ∈ T

where (v_p[i], u_p[i]) ≡ (x,y) is the i^th extreme point of Q, and (v_d[j], u_d[j]) is the j^th extreme direction of Q.

Suppose the dual variable solution of the above main relaxed problem is (̄γ,̄π,̄σ).
The pricing subproblem:

max sum(L[i]*y[i]*̄γ[i] for i in 1:3)+sum((U[j]*y[j]+∇f[j](̄x)*x)*̄π[j] for j in 1:3)+̄σ
(x,y)
st. Ax+By <= 0
    x :Cont
    y :Cont
"""

"""
#What do i need to do
1. decompose milp model
2. build main problem and subproblem
3. how to keep track of constraint sets
4. how to reapetedly update main problem
5. we might include the set of constraints
"z[1] >= y[1]*L[1]
z[2] >= y[2]*L[2]
z[3] >= y[3]*L[3]"
in the subproblem to minimimze the size of the main problem.
6. and of course, every model would use JuMP environment
"""

"""
now consider building the pricing subproblem.
what we need?
1. the constraint set "Ax+By <= 0"
2. the obj function
"sum(L[i]*y[i]*̄γ[i] for i in 1:3)+sum((U[j]*y[j]+∇f[j](̄x)*x)*̄π[j] for j in 1:3)+̄σ"
we can update this function as JuMP expression and define obj function.
"""

"""
we need a data structure for recording data from column generation algorithm
"""

mutable struct lpSolutionObj
    solution  :: Vector{Float64}
    obj_value :: Float64
end

mutable struct dwmSolutionObj
    solution  :: Vector{Float64}
    obj_value :: Float64
end

mutable struct psmSolutionObj
    solution  :: Vector{Float64}
    obj_value :: Float64
end

mutable struct psmExtrmDirSolutionObj
    solution  :: Vector{Float64}
    obj_value :: Float64
end

#define cg data structure
mutable struct CGdata
    lp_model                         :: JuMP.Model
    dw_main_model                    :: JuMP.Model
    pricing_sub_model                :: JuMP.Model
    pricing_sub_extm_dir_model       :: JuMP.Model

    lp_x                             :: Vector{JuMP.VariableRef}
    dw_main_x                        #:: Vector{JuMP.VariableRef}
    pricing_sub_x                    :: Vector{JuMP.VariableRef}
    pricing_sub_extm_dir_x

    lp_num_var                       :: Int64
    dw_main_num_var                  :: Int64
    pricing_sub_num_var              :: Int64

    lp_num_constr                    :: Int64
    dw_main_num_constr               :: Int64
    pricing_sub_num_constr           :: Int64

    lp_le_constr                     :: Vector{Tuple{SAF, MOI.LessThan{Float64}}}
    lp_ge_constr                     :: Vector{Tuple{SAF, MOI.GreaterThan{Float64}}}
    lp_eq_constr                     :: Vector{Tuple{SAF, MOI.EqualTo{Float64}}}
    dw_main_le_constr                :: Vector{Tuple{SAF, MOI.LessThan{Float64}}}
    dw_main_ge_constr                :: Vector{Tuple{SAF, MOI.GreaterThan{Float64}}}
    dw_main_eq_constr                :: Vector{Tuple{SAF, MOI.EqualTo{Float64}}}
    pricing_sub_le_constr            :: Vector{Tuple{SAF, MOI.LessThan{Float64}}}
    pricing_sub_ge_constr            :: Vector{Tuple{SAF, MOI.GreaterThan{Float64}}}
    pricing_sub_eq_constr            :: Vector{Tuple{SAF, MOI.EqualTo{Float64}}}

    lp_obj_sense                     :: Symbol
    dw_main_obj_sense                :: Symbol
    pricing_sub_obj_sense            :: Symbol

    lp_l_var                         :: Vector{Float64}
    lp_u_var                         :: Vector{Float64}
    dw_main_l_var                    :: Vector{Float64}
    dw_main_u_var                    :: Vector{Float64}
    pricing_sub_l_var                :: Vector{Float64}
    pricing_sub_u_var                :: Vector{Float64}

    lp_solution                      :: Vector{Float64}
    dw_main_solution                 :: Vector{Float64}
    dw_main_dual_solution            :: Vector{Float64}
    pricing_sub_solution             :: Vector{Float64}
    pricing_sub_extr_dir_solution    :: Vector{Float64}

    cg_status                        :: Symbol
    cg_started                       :: Bool
    incumbent                        :: Vector{Float64}
    new_incumbent                    :: Bool
    total_time                       :: Float64
    obj_val                          :: Float64
    obj_bound                        :: Float64
    obj_gap                          :: Float64

    num_oa_iter                      :: Int64
    num_cg_iter                      :: Int64

    moi_bin_var                      :: Vector{Any}
    jump_bin_var                     :: Vector{Any}
    obj_terms_lower_bound            :: Vector{Float64}
    obj_terms_upper_bound            :: Vector{Float64}

    dw_main_constr_varidx
    dw_main_constr_coef
    dw_main_constr_type              :: Vector{Symbol}
    dw_main_constr_rhs               :: Vector{Float64}
    dw_main_constr_ref

    pricing_sub_disc2var_idx         :: Vector{Int64}
    pricing_sub_var_type             :: Vector{Symbol}

    λ_counter                        :: Int64
    μ_counter                        :: Int64
    extr_ptn_sol                     :: Vector{psmSolutionObj}
    extr_dir_sol                     :: Vector{psmExtrmDirSolutionObj}
    dw_main_sol                      :: Vector{dwmSolutionObj}

    psm_feasible                     :: Bool
    psm_extr_dir_feasible            :: Bool
    best_solution                    :: Vector{Float64}

    λ_val                            :: Vector{Float64}
    μ_val                            :: Vector{Float64}

    CGdata() = new()
end

function init_cg_data(cgd::CGdata, oad::OAdata)
    cgd.lp_num_var                       = 0
    cgd.dw_main_num_var                  = 0
    cgd.pricing_sub_num_var              = 0

    cgd.lp_num_constr                    = 0
    cgd.dw_main_num_constr               = 0
    cgd.pricing_sub_num_constr           = 0

    cgd.lp_le_constr                     = []
    cgd.lp_ge_constr                     = []
    cgd.lp_eq_constr                     = []
    cgd.dw_main_le_constr                = []
    cgd.dw_main_ge_constr                = []
    cgd.dw_main_eq_constr                = []
    cgd.pricing_sub_le_constr            = []
    cgd.pricing_sub_ge_constr            = []
    cgd.pricing_sub_eq_constr            = []

    cgd.lp_obj_sense                     = :Unknown
    cgd.dw_main_obj_sense                = :Unknown
    cgd.pricing_sub_obj_sense            = :Unknown

    cgd.lp_l_var                         = Float64[]
    cgd.lp_u_var                         = Float64[]
    cgd.dw_main_l_var                    = Float64[]
    cgd.dw_main_u_var                    = Float64[]
    cgd.pricing_sub_l_var                = Float64[]
    cgd.pricing_sub_u_var                = Float64[]

    cgd.lp_solution                      = Float64[]
    cgd.dw_main_solution                 = Float64[]
    cgd.dw_main_dual_solution            = Float64[]
    cgd.pricing_sub_solution             = Float64[]
    cgd.pricing_sub_extr_dir_solution    = Float64[]

    cgd.cg_status                        = :Unknown
    cgd.cg_started                       = false
    cgd.incumbent                        = Float64[]
    cgd.new_incumbent                    = false
    cgd.total_time                       = 0
    cgd.obj_val                          = Inf
    cgd.obj_bound                        = -Inf
    cgd.obj_gap                          = Inf

    cgd.num_oa_iter                      = oad.oa_iter
    cgd.num_cg_iter                      = 0

    cgd.moi_bin_var                      = oad.moi_bin_var
    cgd.jump_bin_var                     = []
    cgd.obj_terms_lower_bound            = oad.obj_terms_lower_bound
    cgd.obj_terms_upper_bound            = oad.obj_terms_upper_bound

    cgd.dw_main_constr_varidx            = oad.ref_linearized_constr_varidx
    cgd.dw_main_constr_coef              = oad.ref_linearized_constr_coef
    cgd.dw_main_constr_type              = oad.ref_linearized_constr_type
    cgd.dw_main_constr_rhs               = oad.ref_linearized_constr_rhs
    cgd.dw_main_constr_ref               = []

    cgd.pricing_sub_disc2var_idx         = []
    cgd.pricing_sub_var_type             = []

    cgd.λ_counter                        = 0
    cgd.μ_counter                        = 0
    cgd.extr_ptn_sol                     = []
    cgd.extr_dir_sol                     = []
    cgd.dw_main_sol                      = []
    cgd.dw_main_x                        = []

    cgd.psm_feasible                     = false
    cgd.psm_extr_dir_feasible            = false
    cgd.best_solution                    = []

    cgd.λ_val                            = []
    cgd.μ_val                            = []
end

function decompose_lp_model(model::JuMP.Model, optimizer::MOI.AbstractOptimizer, op::HybridProblem, oad::OAdata, cgd::CGdata)
    cgd.lp_model = model
    #extract all the constraints from the lp model and save it as cgd lb constraints
    extract_lp_constraints(model, cgd)
    #set user defined lp constraints as the constraints for pricing subproblem
    cgd.pricing_sub_le_constr = optimizer.linear_le_constraints
    cgd.pricing_sub_ge_constr = optimizer.linear_ge_constraints
    cgd.pricing_sub_eq_constr = optimizer.linear_eq_constraints
    #set obj sense
    #obj_sense       = JuMP.objective_sense(model)
    (JuMP.objective_sense(model) == MOI.MAX_SENSE) && (cgd.lp_obj_sense = :MAX_SENSE)
    (JuMP.objective_sense(model) == MOI.MIN_SENSE) && (cgd.lp_obj_sense = :MIN_SENSE)
    cgd.dw_main_obj_sense  = cgd.lp_obj_sense
    (cgd.dw_main_obj_sense == :MIN_SENSE) && (cgd.pricing_sub_obj_sense = :MIN_SENSE)
    (cgd.dw_main_obj_sense == :MAX_SENSE) && (cgd.pricing_sub_obj_sense = :MAX_SENSE)
    #set num var
    cgd.lp_num_var              = JuMP.num_variables(model)
    cgd.pricing_sub_num_var     = length(optimizer.variable_info)
    cgd.pricing_sub_var_type    = oad.ref_var_type[1:cgd.pricing_sub_num_var]
    #set num constr
    cgd.lp_num_constr           = length(cgd.lp_le_constr)+length(cgd.lp_ge_constr)+length(cgd.lp_eq_constr)
    cgd.pricing_sub_num_constr  = length(cgd.pricing_sub_le_constr)+length(cgd.pricing_sub_ge_constr)+length(cgd.pricing_sub_eq_constr)
    cgd.dw_main_num_constr      = cgd.lp_num_constr-cgd.pricing_sub_num_constr+1
    #set var bound
    lp_var = JuMP.all_variables(model)
    for i in 1:cgd.lp_num_var
        JuMP.has_lower_bound(lp_var[i]) && push!(cgd.lp_l_var, JuMP.lower_bound(lp_var[i]))
        !JuMP.has_lower_bound(lp_var[i]) && (op.var_type[i] == :Bin) && push!(cgd.lp_l_var, 0.0)
        JuMP.has_upper_bound(lp_var[i]) && push!(cgd.lp_u_var, JuMP.upper_bound(lp_var[i]))
        !JuMP.has_upper_bound(lp_var[i]) && (op.var_type[i] == :Bin) && push!(cgd.lp_u_var, 1.0)
    end
    #get data from original problem
    cgd.pricing_sub_l_var = op.l_var
    cgd.pricing_sub_u_var = op.u_var
    cgd.pricing_sub_disc2var_idx = op.disc2var_idx
    #decompose the lp model
    construct_pricing_sub_model(cgd)
    construct_dw_main_model(cgd)
end

"""
This function generates initial Dantzig-Wolfe relaxed main problem.
"""
function construct_dw_main_model(cgd::CGdata)
    #initial vals
    nλ = 100 #max 1000 extreme points
    nμ = 100 #max 1000 extreme directions
    nbin         = length(cgd.jump_bin_var)
    #num_constr   = cgd.dw_main_num_constr
    constr_type  = cgd.dw_main_constr_type
    rhs          = cgd.dw_main_constr_rhs
    nitr         = cgd.num_oa_iter
    num_constr   = nbin+nbin*nitr+1
    #the dw relaxed main problem
    dwm = Model()
    @variable(dwm, λ[1:nλ] >= 0)
    @variable(dwm, μ[1:nμ] >= 0)
    @variable(dwm, t[1:nbin])
    #define constr containers
    constr_set1 = []
    constr_set2 = []
    constr_set3 = []
    #define dwd constr for "t >= yL" set of constr
    for i in 1:nbin
        ep = JuMP.@expression(dwm, 0)
        for k in 1:nλ
            cp = 0
            ep = ep+JuMP.@expression(dwm, cp*λ[k])
        end
        ed = JuMP.@expression(dwm, 0)
        for l in 1:nμ
            cd = 0
            ed = ed+JuMP.@expression(dwm, cd*μ[l])
        end
        con=@constraint(dwm, -ep-ed+t[i] >= 0.0)
        push!(constr_set1, con)
    end
    #define dwd constr for "t >= (y-1)U+∇f(̄x)(x-̄x)" set of constr
    indx = 0
    for i in 1:nitr
        for j in 1:nbin
            indx += 1
            ep = JuMP.@expression(dwm, 0)
            for k in 1:nλ
                cp = 0
                ep = ep+JuMP.@expression(dwm, cp*λ[k])
            end
            ed = JuMP.@expression(dwm, 0)
            for l in 1:nμ
                cd = 0
                ed = ed+JuMP.@expression(dwm, cd*μ[l])
            end
            if constr_type[indx] == :(<=)
                con=@constraint(dwm, ep+ed-t[j] <= rhs[indx])
                push!(constr_set2, con)
            elseif constr_type[indx] == :(>=)
                con=@constraint(dwm, ep+ed+t[j] >= rhs[indx])
                push!(constr_set2, con)
            end
        end
    end
    #define convexity constr for "sum(λ[i]) == 1" constr
    con=@constraint(dwm, sum(λ[k]*0 for k in 1:nλ) == 1.0)
    push!(constr_set3, con)
    constr_ref = [constr_set1, constr_set2, constr_set3]
    #define obj func
    (cgd.dw_main_obj_sense == :MIN_SENSE) && @objective(dwm, Min, sum(t[i] for i in 1:nbin))
    (cgd.dw_main_obj_sense == :MAX_SENSE) && @objective(dwm, Max, sum(t[i] for i in 1:nbin))
    cgd.dw_main_constr_ref = constr_ref
    cgd.dw_main_model = dwm
    cgd.dw_main_x = [λ, μ, t]
end

"""
This function constructs the dw pricing sub model.
"""
function construct_pricing_sub_model(cgd)
    lb = cgd.pricing_sub_l_var
    ub = cgd.pricing_sub_u_var
    psm = Model()
    @variable(psm, lb[i] <= x[i=1:cgd.pricing_sub_num_var] <= ub[i])
    backend = JuMP.backend(psm)
    llc = cgd.pricing_sub_le_constr
    lgc = cgd.pricing_sub_ge_constr
    lec = cgd.pricing_sub_eq_constr
    for constr_type in [llc, lgc, lec]
        for constr in constr_type
            MOI.add_constraint(backend, constr[1], constr[2])
        end
    end
    #convert moi variable references to jump variables
    bin_vars = []
    siz = length(cgd.moi_bin_var)
    for i in 1:siz
        push!(bin_vars, JuMP.VariableRef(psm, cgd.moi_bin_var[i].args[2]))
    end
    cgd.jump_bin_var = bin_vars
    cgd.pricing_sub_model = psm
    cgd.pricing_sub_x = x
    #set obj function
    #max sum(L[i]*y[i]*̄γ[i] for i in 1:3)+sum((U[j]*y[j]+∇f[j](̄x)*x)*̄π[j] for j in 1:3)+̄σ
end

"""
This function extracts constraints from the input model of the cg algorithm and saves them in the cg data structure
"""
function extract_lp_constraints(model::JuMP.Model, cgd::CGdata)
    #list of different types of constraints
    constr_types = JuMP.list_of_constraint_types(model)
    #number of different constraint types
    siz = length(constr_types)
    #reserve constraints for each type
    for i in 1:siz
        num_constr = JuMP.num_constraints(model, constr_types[i][1], constr_types[i][2])
        constrs    = JuMP.all_constraints(model, constr_types[i][1], constr_types[i][2])
        for j in 1:num_constr
            con_obj = JuMP.constraint_object(constrs[j])
            func    = JuMP.moi_function(con_obj.func)
            set     = con_obj.set
            (constr_types[i][2] == MOI.EqualTo{Float64}) && push!(cgd.lp_eq_constr, (func, set))
            (constr_types[i][2] == MOI.LessThan{Float64}) && push!(cgd.lp_le_constr, (func, set))
            (constr_types[i][2] == MOI.GreaterThan{Float64}) && push!(cgd.lp_ge_constr, (func, set))
        end
    end
end

"""
This function constructs pricing sub model for extreme directions.
"""
function construct_pricing_sub_model_for_extreme_directions(cgd::CGdata)
    #define model
    psmed = Model()
    @variable(psmed, x[1:cgd.pricing_sub_num_var]) #moi references are created according to the order of variable creation
    @variable(psmed, 0 <= y[1:cgd.pricing_sub_num_var] <= 1)
    @variable(psmed, 0 <= z[1:cgd.pricing_sub_num_var] <= 1)
    for i in 1:cgd.pricing_sub_num_var
        if cgd.pricing_sub_var_type[i] == :Bin
            JuMP.set_lower_bound(x[i], 0)
            JuMP.set_upper_bound(x[i], 1)
        else
            JuMP.set_lower_bound(x[i], -1)
            JuMP.set_upper_bound(x[i], 1)
        end
    end
    #need a function to modify the input set of constraints to the set required for ext direction
    constr_types = JuMP.list_of_constraint_types(cgd.pricing_sub_model)
    #number of different constraint types
    len = length(constr_types)
    #define psmed constraints
    for i in 1:len
        num_constr = JuMP.num_constraints(cgd.pricing_sub_model, constr_types[i][1], constr_types[i][2])
        constrs    = JuMP.all_constraints(cgd.pricing_sub_model, constr_types[i][1], constr_types[i][2])
        for j in 1:num_constr
            con_obj = JuMP.constraint_object(constrs[j])
            func    = JuMP.moi_function(con_obj.func)
            exp = JuMP.@expression(psmed, 0)
            for term in func.terms
                exp = exp+JuMP.@expression(psmed, term.coefficient*(y[term.variable_index.value]-z[term.variable_index.value]))
            end
            (constr_types[i][2] == MOI.LessThan{Float64}) && @constraint(psmed, exp <= 0)
            (constr_types[i][2] == MOI.GreaterThan{Float64}) && @constraint(psmed, -1*exp <= 0)
            if (constr_types[i][2] == MOI.EqualTo{Float64})
                @constraint(psmed, exp <= 0)
                @constraint(psmed, -1*exp <= 0)#need some clarification here, cuz this could be con_obj.set.value
            end
        end
    end
    for i in 1:cgd.pricing_sub_num_var
        @constraint(psmed, x[i] == y[i]-z[i])
    end
    cgd.pricing_sub_extm_dir_model = psmed
    cgd.pricing_sub_extm_dir_x     = [x, y, z]
end

max_stars_repo_name: MdMahbubarRahman/Duran

max_issues_repo_name: MdMahbubarRahman/Duran

max_forks_repo_name: MdMahbubarRahman/Duran


text:
module Chirps

using HDF5, ..Fct

thresholds = [
  5   1.5
  6   4.0
  7   5.9
  8   7.3
  9   8.6
 10   9.7
 11  10.8
 12  11.8
 13  12.8
 14  13.8
 15  14.8
 16  15.8
 17  16.9
]

L=size(thresholds)[1]

offsets = accumulate(+,vcat([1],[1+2^(Int64(thresholds[i,1])-1) for i in 1:L]))

#even-degree Chebyshev coefficients of linear chirps
chirps = Vector{Complex{Float64}}(undef,offsets[end])

degrees = [2 ^ Int64(p) for p in thresholds[:,1]]
rates = 2.0 .^ thresholds[:,2]

function getchirp(i)# 1 <= i <= L
    view(chirps,offsets[i]:offsets[i+1]-1)
end

function store(file)
    h=h5open(file,"w") #overwrite
    for i in 1:L
        N = degrees[i]
        a = BigFloat(rates[i])
        s = Fct.chebsample(x->exp(a*im*x^2),2N,T=Complex{BigFloat})
        g = getchirp(i)
        h[string(N)] = g[:] = Complex{Float64}.(Fct.chebcoeffs(s)[1:2:N+1])
    end
    close(h)
end

function load(file)
    h=h5open(file,"r")
    for i in 1:L
        g = getchirp(i)
        g[:] = read(h, string(degrees[i]))
    end
    close(h)
end

end #module

ext: jl


text:
using DrWatson
@quickactivate :FractalDimension # uses DynamicalSystems, PyPlot
include(srcdir("style.jl"))

N = Int(1e5)

qH = 1
qC = 2

system = :henonheiles_chaotic
data = :embed_system
embedding_ds = 2:7
Cmethod = "standard"

# Calculate values for H, C
eHs, eCs, Hs, Cs = [Vector{Float64}[] for i in 1:4]
for d in embedding_ds
    params = @strdict N d qH qC data system
    if Cmethod ≠ "standard"
        params["Cmethod"] = Cmethod
    end

    params["z"] = -1

    # This is the main call that calculates everything
    output, s = produce_or_load(
        datadir("main"), params, make_C_H;
        prefix = string(data), suffix = "jld2", force = false,
        ignores = ["data"], tag = false,
    )
    @unpack eH, eC, H, C = output
    push!(eHs, eH); push!(Hs, H); push!(eCs, eC); push!(Cs, C)
end

# Do the actual plot
labels = ["\$d=$(d)\$" for d in embedding_ds]
legendtitle = "delay embedding (Hénon-Heiles, Δ ≈ 3)"

fig, axs = mainplot(
    Hs, Cs, eHs, eCs, labels, legendtitle;
    qH, qC, tol = 0.25,
)

wsave(plotsdir("paper", "embedding_$(system)"), fig)

max_stars_repo_name: Datseris/FractalDimension

max_issues_repo_name: Datseris/FractalDimension

max_forks_repo_name: Datseris/FractalDimension


text:
using LinearAlgebra
using Adapt: adapt

import CUDAdrv
import CUDAdrv: CuPtr, CU_NULL

@testset "Memory" begin
  CuArrays.alloc(0)

  @test (CuArrays.@allocated CuArray{Int32}(undef,1)) == 4

  ret, out = @grab_output CuArrays.@time CuArray{Int32}(undef, 1)
  @test isa(ret, CuArray{Int32})
  @test occursin("1 GPU allocation: 4 bytes", out)

  ret, out = @grab_output CuArrays.@time Base.unsafe_wrap(CuArray, CuPtr{Int32}(12345678), (2, 3))
  @test isa(ret, CuArray{Int32})
  @test !occursin("GPU allocation", out)
end

@testset "Array" begin
  xs = CuArray{Int}(undef, 2, 3)
  @test collect(CuArray([1 2; 3 4])) == [1 2; 3 4]
  @test collect(cu[1, 2, 3]) == [1, 2, 3]
  @test collect(cu([1, 2, 3])) == [1, 2, 3]
  @test testf(vec, rand(5,3))
  @test cu(1:3) === 1:3
  @test Base.elsize(xs) == sizeof(Int)
  @test CuArray{Int, 2}(xs) === xs

  @test_throws ArgumentError Base.unsafe_convert(Ptr{Int}, xs)
  @test_throws ArgumentError Base.unsafe_convert(Ptr{Float32}, xs)

  # Check that allowscalar works
  @test_throws ErrorException xs[1]
  @test_throws ErrorException xs[1] = 1

  # unsafe_wrap
  @test Base.unsafe_wrap(CuArray, CU_NULL, 1; own=false).pooled == false
  @test Base.unsafe_wrap(CuArray, CU_NULL, 2)                == CuArray{Nothing,1}(CU_NULL, (2,),  false)
  @test Base.unsafe_wrap(CuArray{Nothing}, CU_NULL, 2)       == CuArray{Nothing,1}(CU_NULL, (2,),  false)
  @test Base.unsafe_wrap(CuArray{Nothing,1}, CU_NULL, 2)     == CuArray{Nothing,1}(CU_NULL, (2,),  false)
  @test Base.unsafe_wrap(CuArray, CU_NULL, (1,2))            == CuArray{Nothing,2}(CU_NULL, (1,2), false)
  @test Base.unsafe_wrap(CuArray{Nothing}, CU_NULL, (1,2))   == CuArray{Nothing,2}(CU_NULL, (1,2), false)
  @test Base.unsafe_wrap(CuArray{Nothing,2}, CU_NULL, (1,2)) == CuArray{Nothing,2}(CU_NULL, (1,2), false)

  @test collect(CuArrays.zeros(2, 2)) == zeros(Float32, 2, 2)
  @test collect(CuArrays.ones(2, 2)) == ones(Float32, 2, 2)

  @test collect(CuArrays.fill(0, 2, 2)) == zeros(Float32, 2, 2)
  @test collect(CuArrays.fill(1, 2, 2)) == ones(Float32, 2, 2)
end

@testset "Adapt" begin
  A = rand(Float32, 3, 3)
  dA = CuArray(A)
  @test adapt(Array, dA) ≈ A
  @test adapt(CuArray, A) ≈ dA
end

@testset "Broadcast" begin
  @test testf((x)       -> fill!(x, 1),  rand(3,3))
  @test testf((x, y)    -> map(+, x, y), rand(2, 3), rand(2, 3))
  @test testf((x)       -> sin.(x),      rand(2, 3))
  @test testf((x)       -> log.(x) .+ 1, rand(2, 3))
  @test testf((x)       -> 2x,           rand(2, 3))
  @test testf((x)       -> x .^ 0,      rand(2, 3))
  @test testf((x)       -> x .^ 1,      rand(2, 3))
  @test testf((x)       -> x .^ 2,      rand(2, 3))
  @test testf((x)       -> x .^ 3,      rand(2, 3))
  @test testf((x)       -> x .^ 5,      rand(2, 3))
  @test testf((x)       -> (z = Int32(5); x .^ z),      rand(2, 3))
  @test testf((x)       -> (z = Float64(π); x .^ z),      rand(2, 3))
  @test testf((x)       -> (z = Float32(π); x .^ z),      rand(Float32, 2, 3))
  @test testf((x, y)    -> x .+ y,       rand(2, 3), rand(1, 3))
  @test testf((z, x, y) -> z .= x .+ y,  rand(2, 3), rand(2, 3), rand(2))
  @test (CuArray{Ptr{Cvoid}}(undef, 1) .= C_NULL) == CuArray([C_NULL])
  @test CuArray([1,2,3]) .+ CuArray([1.0,2.0,3.0]) == CuArray([2,4,6])

  @eval struct Whatever{T}
      x::Int
  end
  @test Array(Whatever{Int}.(CuArray([1]))) == Whatever{Int}.([1])
end

@testset "Cufunc" begin
  gelu1(x) = oftype(x, 0.5) * x * (1 + tanh(oftype(x, √(2/π))*(x + oftype(x, 0.044715) * x^3)))
  sig(x) = one(x) / (one(x) + exp(-x))
  f(x) = gelu1(log(x)) * sig(x) * tanh(x)
  g(x) = x^7 - 2 * x^f(x^2) + 3

  CuArrays.@cufunc gelu1(x) = oftype(x, 0.5) * x * (1 + tanh(oftype(x, √(2/π))*(x + oftype(x, 0.044715) * x^3)))
  CuArrays.@cufunc sig(x) = one(x) / (one(x) + exp(-x))
  CuArrays.@cufunc f(x) = gelu1(log(x)) * sig(x) * tanh(x)
  CuArrays.@cufunc g(x) = x^7 - 2 * x^f(x^2) + 3

  @test :gelu1 ∈ CuArrays.cufuncs()
  @test :sig ∈ CuArrays.cufuncs()
  @test :f ∈ CuArrays.cufuncs()
  @test :g ∈ CuArrays.cufuncs()
  @test testf((x)  -> gelu1.(x), rand(3,3))
  @test testf((x)  -> sig.(x),  rand(3,3))
  @test testf((x)  -> f.(x),    rand(3,3))
  @test testf((x)  -> g.(x),    rand(3,3))
end

# https://github.com/JuliaGPU/CUDAnative.jl/issues/223
@testset "Ref Broadcast" begin
  foobar(idx, A) = A[idx]
  @test CuArray([42]) == foobar.(CuArray([1]), Base.RefValue(CuArray([42])))
end

@testset "Broadcast Fix" begin
  @test testf(x -> log.(x), rand(3,3))
  @test testf((x,xs) -> log.(x.+xs), Ref(1), rand(3,3))

  if CuArrays.has_cudnn()
    using NNlib

    @test testf(x -> logσ.(x), rand(5))
  end
end

@testset "Reduce" begin
  @test testf(x -> sum(x, dims=1), rand(2, 3))
  @test testf(x -> sum(x, dims=2), rand(2, 3))
  @test testf(x -> sum(x -> x^2, x, dims=1), rand(2, 3))
  @test testf(x -> prod(x, dims=2), rand(2, 3))

  @test testf(x -> sum(x), rand(2, 3))
  @test testf(x -> prod(x), rand(2, 3))

  @test testf(x -> minimum(x), rand(2, 3))
  @test testf(x -> minimum(x, dims=2), rand(2, 3))
  @test testf(x -> minimum(x, dims=(2, 3)), rand(2, 3, 4))

  @test testf(x -> maximum(x), rand(2, 3))
  @test testf(x -> maximum(x, dims=2), rand(2, 3))
  @test testf(x -> maximum(x, dims=(2, 3)), rand(2, 3, 4))

  myreducer(x1, x2) = x1 + x2 # bypass optimisations for sum()
  @test testf(x -> reduce(myreducer, x, dims=(2, 3), init=0.0), rand(2, 3, 4))
  @test testf(x -> reduce(myreducer, x, init=0.0), rand(2, 3))
  @test testf(x -> reduce(myreducer, x, dims=2, init=42.0), rand(2, 3))

  ex = ErrorException("Please supply a neutral element for &. E.g: mapreduce(f, &, A; init = 1)")
  @test_throws ex mapreduce(t -> t > 0.5, &, cu(rand(2, 3)))
  @test testf(x -> mapreduce(t -> t > 0.5, &, x, init=true), rand(2, 3))

  ex = UndefKeywordError(:init)
  cub = map(t -> t > 0.5, cu(rand(2, 3)))
  @test_throws ex reduce(|, cub)
  @test testf(x -> reduce(|, x, init=false), map(t -> t > 0.5, cu(rand(2, 3))))
end

@testset "0D" begin
  x = CuArray{Float64}(undef)
  x .= 1
  @test collect(x)[] == 1
  if VERSION >= v"1.3-"
    # broken test that throws
    # https://github.com/JuliaGPU/GPUArrays.jl/issues/204
    @test_throws CUDAnative.InvalidIRError x /= 2
  else
    x /= 2
    @test collect(x)[] == 0.5
  end
end

@testset "SubArray" begin
  @test testf(rand(5)) do x
    y = x[2:4]
    y .= 1
    x
  end

  @test testf(rand(5)) do x
    y = view(x, 2:4)
    y .= 1
    x
  end

  @test testf(x->view(x, :, 1:4, 3), rand(Float32, 5, 4, 3))

  let x = cu(rand(Float32, 5, 4, 3))
    @test_throws BoundsError view(x, :, :, 1:10)

    # Contiguous views should return new CuArray
    @test view(x, :, 1, 2) isa CuVector{Float32}
    @test view(x, 1:4, 1, 2) isa CuVector{Float32}
    @test view(x, :, 1:4, 3) isa CuMatrix{Float32}
    @test view(x, :, :, 1) isa CuMatrix{Float32}
    @test view(x, :, :, :) isa CuArray{Float32,3}
    @test view(x, :) isa CuVector{Float32}
    @test view(x, 1:3) isa CuVector{Float32}

    # Non-contiguous views should fall back to base's SubArray
    @test view(x, 1:3, 1:3, 3) isa SubArray
    @test view(x, 1, :, 3) isa SubArray
    @test view(x, 1, 1:4, 3) isa SubArray
    @test view(x, :, 1, 1:3) isa SubArray
    @test view(x, :, 1:2:4, 1) isa SubArray
    @test view(x, 1:2:5, 1, 1) isa SubArray
  end

  # non-contiguous copyto!
  let x = CuArrays.rand(4, 4)
    y = view(x, 2:3, 2:3)

    # to gpu
    gpu = CuArray{eltype(y)}(undef, size(y))
    copyto!(gpu, y)
    @test Array(gpu) == Array(y)

    # to cpu
    cpu = Array{eltype(y)}(undef, size(y))
    copyto!(cpu, y)
    @test cpu == Array(y)
  end

  # bug in parentindices conversion
  let x = CuArray{Int}(undef, 1, 1)
    x[1,:] .= 42
    @test Array(x)[1,1] == 42
  end

  # bug in copyto!
  ## needless N type parameter
  @test testf((x,y)->copyto!(y, selectdim(x, 2, 1)), ones(2,2,2), zeros(2,2))
  ## inability to copyto! smaller destination
  @test testf((x,y)->copyto!(y, selectdim(x, 2, 1)), ones(2,2,2), zeros(3,3))
end

@testset "reshape" begin
  A = [1 2 3 4
       5 6 7 8]
  gA = reshape(CuArray(A),1,8)
  _A = reshape(A,1,8)
  _gA = Array(gA)
  @test all(_A .== _gA)
  A = [1,2,3,4]
  gA = reshape(CuArray(A),4)
end

@testset "$f! with diagonal $d" for (f, f!) in ((triu, triu!), (tril, tril!)),
                                          d in -2:2
  A = randn(10, 10)
  @test f(A, d) == Array(f!(CuArray(A), d))
end

@testset "Utilities" begin
  t = Base.@elapsed ret = CuArrays.@sync begin
    # TODO: do something that takes a while on the GPU
    #       (need to wrap clock64 in CUDAnative for that)
    42
  end
  @test t >= 0
  @test ret == 42
end

@testset "accumulate" begin
  for n in (0, 1, 2, 3, 10, 10_000, 16384, 16384+1) # small, large, odd & even, pow2 and not
    @test testf(x->accumulate(+, x), rand(n))
  end

  # multidimensional
  for (sizes, dims) in ((2,) => 2,
                        (3,4,5) => 2,
                        (1, 70, 50, 20) => 3)
    @test testf(x->accumulate(+, x; dims=dims), rand(Int, sizes))
  end

  # using initializer
  for (sizes, dims) in ((2,) => 2,
                        (3,4,5) => 2,
                        (1, 70, 50, 20) => 3)
    @test testf(x->accumulate(+, x; dims=dims, init=100.), rand(Int, sizes))
  end

  # in place
  @test testf(x->(accumulate!(+, x, copy(x)); x), rand(2))

  # specialized
  @test testf(cumsum, rand(2))
  @test testf(cumprod, rand(2))
end

@testset "logical indexing" begin
  @test CuArray{Int}(undef, 2)[CuArray{Bool}(undef, 2)] isa CuArray
  @test CuArray{Int}(undef, 2, 2)[CuArray{Bool}(undef, 2, 2)] isa CuArray
  @test CuArray{Int}(undef, 2, 2, 2)[CuArray{Bool}(undef, 2, 2, 2)] isa CuArray

  @test CuArray{Int}(undef, 2)[Array{Bool}(undef, 2)] isa CuArray
  @test CuArray{Int}(undef, 2, 2)[Array{Bool}(undef, 2, 2)] isa CuArray
  @test CuArray{Int}(undef, 2, 2, 2)[Array{Bool}(undef, 2, 2, 2)] isa CuArray

  @test testf((x,y)->x[y], rand(2), rand(Bool, 2))
  @test testf((x,y)->x[y], rand(2, 2), rand(Bool, 2, 2))
  @test testf((x,y)->x[y], rand(2, 2, 2), rand(Bool, 2, 2, 2))

  @test testf(x -> x[x .> 0.5], rand(2))
  @test testf(x -> x[x .> 0.5], rand(2,2))
  @test testf(x -> x[x .> 0.5], rand(2,2,2))

  @test testf(x -> filter(y->y .> 0.5, x), rand(2))
  @test testf(x -> filter(y->y .> 0.5, x), rand(2,2))
  @test testf(x -> filter(y->y .> 0.5, x), rand(2,2,2))
end

@testset "generic fallbacks" begin
    a = rand(Int8, 3, 3)
    b = rand(Int8, 3, 3)
    d_a = CuArray{Int8}(a)
    d_b = CuArray{Int8}(b)
    d_c = d_a*d_b
    @test collect(d_c) == a*b
    a = rand(Complex{Int8}, 3, 3)
    b = rand(Complex{Int8}, 3, 3)
    d_a = CuArray{Complex{Int8}}(a)
    d_b = CuArray{Complex{Int8}}(b)
    d_c = d_a'*d_b
    @test collect(d_c) == a'*b
    d_c = d_a*d_b'
    @test collect(d_c) == a*b'
    d_c = d_a'*d_b'
    @test collect(d_c) == a'*b'
    d_c = transpose(d_a)*d_b'
    @test collect(d_c) == transpose(a)*b'
    d_c = d_a'*transpose(d_b)
    @test collect(d_c) == a'*transpose(b)
    d_c = transpose(d_a)*d_b
    @test collect(d_c) == transpose(a)*b
    d_c = d_a*transpose(d_b)
    @test collect(d_c) == a*transpose(b)
    d_c = transpose(d_a)*transpose(d_b)
    @test collect(d_c) == transpose(a)*transpose(b)
    d_c = rmul!(copy(d_a), Complex{Int8}(2, 2))
    @test collect(d_c) == a*Complex{Int8}(2, 2)
    d_c = lmul!(Complex{Int8}(2, 2), copy(d_a))
    @test collect(d_c) == Complex{Int8}(2, 2)*a
end

@testset "reverse" begin
    # 1-d out-of-place
    @test testf(x->reverse(x), rand(1000))
    @test testf(x->reverse(x, 10), rand(1000))
    @test testf(x->reverse(x, 10, 90), rand(1000))

    # 1-d in-place
    @test testf(x->reverse!(x), rand(1000))
    @test testf(x->reverse!(x, 10), rand(1000))
    @test testf(x->reverse!(x, 10, 90), rand(1000))

    # n-d out-of-place
    for shape in ([1, 2, 4, 3], [4, 2], [5], [2^5, 2^5, 2^5]),
        dim in 1:length(shape)
      @test testf(x->reverse(x; dims=dim), rand(shape...))

      cpu = rand(shape...)
      gpu = CuArray(cpu)
      reverse!(gpu; dims=dim)
      @test Array(gpu) == reverse(cpu; dims=dim)
    end
end

@testset "permutedims" begin
    @test testf(x->permutedims(x, [1, 2]), rand(4, 4))

    inds = rand(1:100, 150, 150)
    @test testf(x->permutedims(view(x, inds, :), (3, 2, 1)), rand(100, 100))
end

@testset "findall" begin
    # 1D
    @test testf(x->findall(x), rand(Bool, 100))
    @test testf(x->findall(y->y>0.5, x), rand(100))

    # ND
    let x = rand(Bool, 10, 10)
      @test findall(x) == Array(findall(CuArray(x)))
    end
    let x = rand(10, 10)
      @test findall(y->y>0.5, x) == Array(findall(y->y>0.5, CuArray(x)))
    end
end

@testset "findfirst" begin
    # 1D
    @test testf(x->findfirst(x), rand(Bool, 100))
    @test testf(x->findfirst(y->y>0.5, x), rand(100))

    # ND
    let x = rand(Bool, 10, 10)
      @test findfirst(x) == findfirst(CuArray(x))
    end
    let x = rand(10, 10)
      @test findfirst(y->y>0.5, x) == findfirst(y->y>0.5, CuArray(x))
    end
end

@testset "findmin" begin
  let x = rand(Float32, 100)
      @test findmin(x) == findmin(CuArray(x))
  end
  let x = rand(Float32, 10, 10)
      @test findmin(x) == findmin(CuArray(x))
  end
end

@testset "findmax" begin
  let x = rand(Float32, 100)
      @test findmax(x) == findmax(CuArray(x))
  end
  let x = rand(Float32, 10, 10)
      @test findmax(x) == findmax(CuArray(x))
  end
end


text:
# include("solvers/taboo.jl")

"""
    main_test(args)

Cette action est destinée à la mise au point de code (nouveau solveurs,
nouveaux voisinages, ...)

Une fois la fonctionnalité opérationnel, il faut en faire en test unitaire
(i.e. en créant un nouveau fichier dans tests/text_xxx.jl)
"""
function main_test(args)
    # p1 = Plane("p2", 1, 100, 150, 200, 1.0, 2.0)
    # p2 = Plane("p22", 1, 100, 150, 200, 1.0)

    if Args.get(:infile) == "NO_INFILE"
        Args.set(:infile, "data/alpx/01.alpx")
    end
    println("="^70)
    println("Début de l'action test")
    inst = Instance(Args.get(:infile))

    sol = Solution(inst)
    println(to_s(sol)) # OK
    solve!(sol)
    println(to_s(sol)) # OK
    Random.shuffle!(sol.planes)
    println(to_s(sol)) # OK
    solve!(sol)
    println(to_s(sol)) # OK

    println("Fin de l'action test")
end

main_test(Args.args)


text:
struct AngularVelocity{T} <: FieldVector{3, T}
    x1::T
    x2::T
    x3::T
end
AngularVelocity(x1, x2, x3) = AngularVelocity(Base.promote(x1, x2, x3)...)

@inline (::Type{R})(ω::AngularVelocity) where {R<:RotationVec} = R(ω...)
@inline (::Type{R})(ω::AngularVelocity) where {R<:Rotation{3}} = R(RotationVec(ω))

@inline (::Type{T})(r::RotationVec) where {T<:AngularVelocity} = T(r.sx, r.sy, r.sz)
@inline (::Type{T})(r::Rotation{3}) where {T<:AngularVelocity} = T(RotationVec(r))

groupexp(ω::AngularVelocity) = RotationVec(ω)

grouplog(r::Rotation{3}) = AngularVelocity(r)

Rotations.rotation_angle(ω::AngularVelocity) = norm(ω)


text:
function tcm_estimateRanking_kcv_weighted_ass(N, CS_list, nfolds, assortment_folds, assortments, transaction_counts, numInitializations, time_limit, perTol)

	unique_fold_ids = unique(assortment_folds)

	KL_total_holdout_fold = zeros(length(CS_list), nfolds)
	noinf_count_holdout_fold = zeros(length(CS_list), nfolds)

	start_time = time()

	for fold_ind in 1:nfolds
		fold = unique_fold_ids[fold_ind]

		# trainfolds_assortments, trainfolds_transaction_counts = tcm_transactionsToCounts(N, trainfolds_transactions)
		# testfold_assortments, testfold_transaction_counts = tcm_transactionsToCounts(N, testfold_transactions)

		trainfolds_assortments = assortments[ find(assortment_folds .!= fold), :]
		testfold_assortments = assortments[ find(assortment_folds .== fold), :]

		trainfolds_transaction_counts = transaction_counts[ find(assortment_folds .!= fold), :]
		testfold_transaction_counts = transaction_counts[ find(assortment_folds .== fold), :]

		for CSind in 1:length(CS_list)
			println("CS size ", CS_list[CSind], ", fold ", fold)
			CS_size = CS_list[CSind]

			ranking_lambda, orderings, ranking_loglik, ranking_elapsed_time = tcm_estimateRanking_CS(N, trainfolds_assortments, trainfolds_transaction_counts, time_limit, perTol, CS_size)
			ranking_K = length(ranking_lambda)

			testfold_ranking_predict, testfold_ranking_predict_by_class = tcm_predictRanking(N, ranking_K, ranking_lambda, orderings, testfold_assortments)

			testfold_KL, testfold_noinf_count, testfold_KL_avg = tcm_evaluateKL(N, testfold_ranking_predict, testfold_transaction_counts)

			KL_total_holdout_fold[CSind, fold_ind] = testfold_KL #testfold_KL / sum(testfold_transaction_counts)
			noinf_count_holdout_fold[CSind, fold_ind] = testfold_noinf_count
		end
	end

	end_time = time()
	elapsed_time = end_time - start_time

	KL_avg_by_CS = sum(KL_total_holdout_fold, 2) ./ sum(noinf_count_holdout_fold, 2)  * 100

	best_CSind = indmin(KL_avg_by_CS)
	best_CS_size = CS_list[best_CSind]
	best_KL_avg_cv = KL_avg_by_CS[best_CSind]

	@show KL_avg_by_CS

	return best_CS_size, best_KL_avg_cv, elapsed_time
end

function tcm_estimateRanking_kcv_weighted(N, CS_list, nfolds, transaction_folds, transactions, time_limit, perTol)
	unique_fold_ids = unique(transaction_folds)

	KL_total_holdout_fold = zeros(length(CS_list), nfolds)
	noinf_count_holdout_fold = zeros(length(CS_list), nfolds)

	start_time = time()

	for fold_ind in 1:nfolds
		fold = unique_fold_ids[fold_ind]

		@show fold

		@show unique_fold_ids

		println("Pre-subsetting ")

		#@show transaction_folds
		#temp = find(transaction_folds .!= fold)
		#@show temp
		@show size(transaction_folds)
		@show size(transactions)
		@show fold

		temp = find(transaction_folds .!= fold)
		@show temp[1:5]
		trainfolds_transactions = transactions[ find(transaction_folds .!= fold), :]
		println("Post-subsetting")
		testfold_transactions = transactions[ find(transaction_folds .== fold), :]

		trainfolds_assortments, trainfolds_transaction_counts = tcm_transactionsToCounts(N, trainfolds_transactions)
		testfold_assortments, testfold_transaction_counts = tcm_transactionsToCounts(N, testfold_transactions)

		for i in 1:length(CS_list)
			println("CS size ", CS_list[i], ", fold ", fold)
			input_CS_size = CS_list[i]

			lambda, orderings, loglik, elapsed_time = tcm_estimateRanking_CS(N, trainfolds_assortments, trainfolds_transaction_counts, time_limit, perTol, input_CS_size)
			K = length(lambda)
			testfold_ranking_predict, testfold_ranking_predict_by_class = tcm_predictRanking(N, K, lambda, orderings, testfold_assortments)

			testfold_KL, testfold_noinf_count, testfold_KL_avg = tcm_evaluateKL(N, testfold_ranking_predict, testfold_transaction_counts)

			KL_total_holdout_fold[i, fold_ind] = testfold_KL #testfold_KL / sum(testfold_transaction_counts)
			noinf_count_holdout_fold[i, fold_ind] = testfold_noinf_count
		end
	end

	end_time = time()
	elapsed_time = end_time - start_time

	KL_avg_by_CS_size = sum(KL_total_holdout_fold, 2) ./ sum(noinf_count_holdout_fold, 2)  * 100

	@show KL_total_holdout_fold
	@show noinf_count_holdout_fold
	@show KL_avg_by_CS_size

	best_i = indmin(KL_avg_by_CS_size)
	best_CS_size = CS_list[best_i]
	best_KL_avg_cv = KL_avg_by_CS_size[best_i]

	return best_CS_size, best_KL_avg_cv, elapsed_time

end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
using Test, AlgoritmiStruttureDati.Laboratorio3
a = [4, 3, 2, 1]
b = [10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0]

@testset "countInc.jl" begin
    @test countInc(a) == 3
    @test countInc(b) == 10
end


text:
const GLOBAL_DEBUG = 0

import Gallium

macro if_debug(expr)
    return quote
        if GLOBAL_DEBUG > 0 && isinteractive()
            $(esc(expr))
        end
    end
end

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
const empty_vector_of_array = VectorOfArray(Array{Vector{Float64}}(Float64[]))
Base.@kwdef struct model_states{vec1D,vec2D,R1}
    """
    If you want to add anything to this struct, you must also check/modify set_vars!`.
    Otherwise, it may not work as intended
    """
    # Matrices (vectors in space and time)
    Y::vec2D = ( vec2D == VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
    YP::vec2D = ( vec2D == VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
    c_e::vec2D = ( vec2D == VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
    c_s_avg::vec2D = ( vec2D == VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
    T::vec2D = ( vec2D == VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
    film::vec2D = ( vec2D == VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
    Q::vec2D = ( vec2D == VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
    j::vec2D = ( vec2D == VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
    j_s::vec2D = ( vec2D == VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
    Φ_e::vec2D = ( vec2D == VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
    Φ_s::vec2D = ( vec2D == VectorOfArray{Float64,2,Array{Array{Float64,1},1}} ) ? copy(empty_vector_of_array) : nothing
    # Vectors (vectors in time, not space)
    I::vec1D = ( vec1D == Array{Float64,1} ) ? Float64[] : nothing
    t::vec1D = ( vec1D == Array{Float64,1} ) ? Float64[] : nothing
    V::vec1D = ( vec1D == Array{Float64,1} ) ? Float64[] : nothing
    P::vec1D = ( vec1D == Array{Float64,1} ) ? Float64[] : nothing
    SOC::vec1D = ( vec1D == Array{Float64,1} ) ? Float64[] : nothing
    SOH::vec1D = ( vec1D == Array{Float64,1} ) ? Float64[] : nothing
    # Culmination of all the runs
    results::R1 = R1 == Array{run_results,1} ? run_results[] : nothing
end

@inline model_states_logic(output::Symbol) = model_states_logic((output,))
@inline model_states_logic() = model_states_logic(())

eval(quote
@inline function model_states_logic(outputs::T) where T<:Tuple
    outputs_tot = $(fieldnames(model_states)[1:end-1])

    use_all = :all ∈ outputs

    x = @inbounds (use_all || field ∈ outputs for field in outputs_tot)

    return model_states{Bool,Bool,T}(x..., outputs), outputs
end
end)


text:
using Parameters

export hamiltonian, u, ω, ϑ, ϑ₁, ϑ₂, ϑ₃, ϑ₄, dϑ, β₁, β₂, β₃, B, B₁, B₂, B₃, b₁, b₂, b₃, dH

@inline u(t,q) = q[4]

ϑ₁(t,q) = A₁(t,q) + u(t,q) * b₁(t,q)
ϑ₂(t,q) = A₂(t,q) + u(t,q) * b₂(t,q)
ϑ₃(t,q) = A₃(t,q) + u(t,q) * b₃(t,q)
ϑ₄(t,q) = zero(eltype(q))

dϑ₁dx₁(t,q) = dA₁dx₁(t,q) + u(t,q) * db₁dx₁(t,q)
dϑ₁dx₂(t,q) = dA₁dx₂(t,q) + u(t,q) * db₁dx₂(t,q)
dϑ₁dx₃(t,q) = dA₁dx₃(t,q) + u(t,q) * db₁dx₃(t,q)
dϑ₁dx₄(t,q) = b₁(t,q)

dϑ₂dx₁(t,q) = dA₂dx₁(t,q) + u(t,q) * db₂dx₁(t,q)
dϑ₂dx₂(t,q) = dA₂dx₂(t,q) + u(t,q) * db₂dx₂(t,q)
dϑ₂dx₃(t,q) = dA₂dx₃(t,q) + u(t,q) * db₂dx₃(t,q)
dϑ₂dx₄(t,q) = b₂(t,q)

dϑ₃dx₁(t,q) = dA₃dx₁(t,q) + u(t,q) * db₃dx₁(t,q)
dϑ₃dx₂(t,q) = dA₃dx₂(t,q) + u(t,q) * db₃dx₂(t,q)
dϑ₃dx₃(t,q) = dA₃dx₃(t,q) + u(t,q) * db₃dx₃(t,q)
dϑ₃dx₄(t,q) = b₃(t,q)

dϑ₄dx₁(t,q) = zero(eltype(q))
dϑ₄dx₂(t,q) = zero(eltype(q))
dϑ₄dx₃(t,q) = zero(eltype(q))
dϑ₄dx₄(t,q) = zero(eltype(q))

function ϑ(t::Number, q::AbstractVector, Θ::AbstractVector)
    Θ[1] = ϑ₁(t,q)
    Θ[2] = ϑ₂(t,q)
    Θ[3] = ϑ₃(t,q)
    Θ[4] = ϑ₄(t,q)
    nothing
end

function ϑ(t::Number, q::AbstractVector, k::Int)
    if k == 1
        ϑ₁(t, q)
    elseif k == 2
        ϑ₂(t, q)
    elseif k == 3
        ϑ₃(t, q)
    elseif k == 4
        ϑ₄(t, q)
    else
        nothing
    end
end

function ω(t, q, Ω)
    Ω[1,1] = 0
    Ω[1,2] = dϑ₁dx₂(t,q) - dϑ₂dx₁(t,q)
    Ω[1,3] = dϑ₁dx₃(t,q) - dϑ₃dx₁(t,q)
    Ω[1,4] = dϑ₁dx₄(t,q) - dϑ₄dx₁(t,q)

    Ω[2,1] = dϑ₂dx₁(t,q) - dϑ₁dx₂(t,q)
    Ω[2,2] = 0
    Ω[2,3] = dϑ₂dx₃(t,q) - dϑ₃dx₂(t,q)
    Ω[2,4] = dϑ₂dx₄(t,q) - dϑ₄dx₂(t,q)

    Ω[3,1] = dϑ₃dx₁(t,q) - dϑ₁dx₃(t,q)
    Ω[3,2] = dϑ₃dx₂(t,q) - dϑ₂dx₃(t,q)
    Ω[3,3] = 0
    Ω[3,4] = dϑ₃dx₄(t,q) - dϑ₄dx₃(t,q)

    Ω[4,1] = dϑ₄dx₁(t,q) - dϑ₁dx₄(t,q)
    Ω[4,2] = dϑ₄dx₂(t,q) - dϑ₂dx₄(t,q)
    Ω[4,3] = dϑ₄dx₃(t,q) - dϑ₃dx₄(t,q)
    Ω[4,4] = 0

    nothing
end

function D²ϑ₁(t, q, D²ϑ)
    D²ϑ[1,1] = d²A₁dx₁dx₁(t,q) + u(t,q) * d²b₁dx₁dx₁(t,q)
    D²ϑ[1,2] = d²A₁dx₁dx₂(t,q) + u(t,q) * d²b₁dx₁dx₂(t,q)
    D²ϑ[1,3] = d²A₁dx₁dx₃(t,q) + u(t,q) * d²b₁dx₁dx₃(t,q)
    D²ϑ[1,4] = db₁dx₁(t,q)

    D²ϑ[2,1] = d²A₁dx₂dx₁(t,q) + u(t,q) * d²b₁dx₂dx₁(t,q)
    D²ϑ[2,2] = d²A₁dx₂dx₂(t,q) + u(t,q) * d²b₁dx₂dx₂(t,q)
    D²ϑ[2,3] = d²A₁dx₂dx₃(t,q) + u(t,q) * d²b₁dx₂dx₃(t,q)
    D²ϑ[2,4] = db₁dx₂(t,q)

    D²ϑ[3,1] = d²A₁dx₃dx₁(t,q) + u(t,q) * d²b₁dx₃dx₁(t,q)
    D²ϑ[3,2] = d²A₁dx₃dx₂(t,q) + u(t,q) * d²b₁dx₃dx₂(t,q)
    D²ϑ[3,3] = d²A₁dx₃dx₃(t,q) + u(t,q) * d²b₁dx₃dx₃(t,q)
    D²ϑ[3,4] = db₁dx₃(t,q)

    D²ϑ[4,1] = db₁dx₁(t,q)
    D²ϑ[4,2] = db₁dx₂(t,q)
    D²ϑ[4,3] = db₁dx₃(t,q)
    D²ϑ[4,4] = 0

    nothing
end

function D²ϑ₂(t, q, D²ϑ)
    D²ϑ[1,1] = d²A₂dx₁dx₁(t,q) + u(t,q) * d²b₂dx₁dx₁(t,q)
    D²ϑ[1,2] = d²A₂dx₁dx₂(t,q) + u(t,q) * d²b₂dx₁dx₂(t,q)
    D²ϑ[1,3] = d²A₂dx₁dx₃(t,q) + u(t,q) * d²b₂dx₁dx₃(t,q)
    D²ϑ[1,4] = db₂dx₁(t,q)

    D²ϑ[2,1] = d²A₂dx₂dx₁(t,q) + u(t,q) * d²b₂dx₂dx₁(t,q)
    D²ϑ[2,2] = d²A₂dx₂dx₂(t,q) + u(t,q) * d²b₂dx₂dx₂(t,q)
    D²ϑ[2,3] = d²A₂dx₂dx₃(t,q) + u(t,q) * d²b₂dx₂dx₃(t,q)
    D²ϑ[2,4] = db₂dx₂(t,q)

    D²ϑ[3,1] = d²A₂dx₃dx₁(t,q) + u(t,q) * d²b₂dx₃dx₁(t,q)
    D²ϑ[3,2] = d²A₂dx₃dx₂(t,q) + u(t,q) * d²b₂dx₃dx₂(t,q)
    D²ϑ[3,3] = d²A₂dx₃dx₃(t,q) + u(t,q) * d²b₂dx₃dx₃(t,q)
    D²ϑ[3,4] = db₂dx₃(t,q)

    D²ϑ[4,1] = db₂dx₁(t,q)
    D²ϑ[4,2] = db₂dx₂(t,q)
    D²ϑ[4,3] = db₂dx₃(t,q)
    D²ϑ[4,4] = 0

    nothing
end

function D²ϑ₃(t, q, D²ϑ)
    D²ϑ[1,1] = d²A₃dx₁dx₁(t,q) + u(t,q) * d²b₃dx₁dx₁(t,q)
    D²ϑ[1,2] = d²A₃dx₁dx₂(t,q) + u(t,q) * d²b₃dx₁dx₂(t,q)
    D²ϑ[1,3] = d²A₃dx₁dx₃(t,q) + u(t,q) * d²b₃dx₁dx₃(t,q)
    D²ϑ[1,4] = db₃dx₁(t,q)

    D²ϑ[2,1] = d²A₃dx₂dx₁(t,q) + u(t,q) * d²b₃dx₂dx₁(t,q)
    D²ϑ[2,2] = d²A₃dx₂dx₂(t,q) + u(t,q) * d²b₃dx₂dx₂(t,q)
    D²ϑ[2,3] = d²A₃dx₂dx₃(t,q) + u(t,q) * d²b₃dx₂dx₃(t,q)
    D²ϑ[2,4] = db₃dx₂(t,q)

    D²ϑ[3,1] = d²A₃dx₃dx₁(t,q) + u(t,q) * d²b₃dx₃dx₁(t,q)
    D²ϑ[3,2] = d²A₃dx₃dx₂(t,q) + u(t,q) * d²b₃dx₃dx₂(t,q)
    D²ϑ[3,3] = d²A₃dx₃dx₃(t,q) + u(t,q) * d²b₃dx₃dx₃(t,q)
    D²ϑ[3,4] = db₃dx₃(t,q)

    D²ϑ[4,1] = db₃dx₁(t,q)
    D²ϑ[4,2] = db₃dx₂(t,q)
    D²ϑ[4,3] = db₃dx₃(t,q)
    D²ϑ[4,4] = 0

    nothing
end

function D²ϑ₄(t, q, D²ϑ)
    D²ϑ .= 0
    nothing
end

function D²ϑd₁(t, q, D²ϑ)
    D²ϑ[1,1] = d²A₁dx₁dx₁(t,q) + u(t,q) * d²b₁dx₁dx₁(t,q)
    D²ϑ[1,2] = d²A₁dx₁dx₂(t,q) + u(t,q) * d²b₁dx₁dx₂(t,q)
    D²ϑ[1,3] = d²A₁dx₁dx₃(t,q) + u(t,q) * d²b₁dx₁dx₃(t,q)
    D²ϑ[1,4] = db₁dx₁(t,q)

    D²ϑ[2,1] = d²A₂dx₁dx₁(t,q) + u(t,q) * d²b₂dx₁dx₁(t,q)
    D²ϑ[2,2] = d²A₂dx₁dx₂(t,q) + u(t,q) * d²b₂dx₁dx₂(t,q)
    D²ϑ[2,3] = d²A₂dx₁dx₃(t,q) + u(t,q) * d²b₂dx₁dx₃(t,q)
    D²ϑ[2,4] = db₂dx₁(t,q)

    D²ϑ[3,1] = d²A₃dx₁dx₁(t,q) + u(t,q) * d²b₃dx₁dx₁(t,q)
    D²ϑ[3,2] = d²A₃dx₁dx₂(t,q) + u(t,q) * d²b₃dx₁dx₂(t,q)
    D²ϑ[3,3] = d²A₃dx₁dx₃(t,q) + u(t,q) * d²b₃dx₁dx₃(t,q)
    D²ϑ[3,4] = db₃dx₁(t,q)

    D²ϑ[4,1] = 0
    D²ϑ[4,2] = 0
    D²ϑ[4,3] = 0
    D²ϑ[4,4] = 0

    nothing
end

function D²ϑd₂(t, q, D²ϑ)
    D²ϑ[2,1] = d²A₁dx₂dx₁(t,q) + u(t,q) * d²b₁dx₂dx₁(t,q)
    D²ϑ[2,2] = d²A₁dx₂dx₂(t,q) + u(t,q) * d²b₁dx₂dx₂(t,q)
    D²ϑ[2,3] = d²A₁dx₂dx₃(t,q) + u(t,q) * d²b₁dx₂dx₃(t,q)
    D²ϑ[2,4] = db₁dx₂(t,q)

    D²ϑ[2,1] = d²A₂dx₂dx₁(t,q) + u(t,q) * d²b₂dx₂dx₁(t,q)
    D²ϑ[2,2] = d²A₂dx₂dx₂(t,q) + u(t,q) * d²b₂dx₂dx₂(t,q)
    D²ϑ[2,3] = d²A₂dx₂dx₃(t,q) + u(t,q) * d²b₂dx₂dx₃(t,q)
    D²ϑ[2,4] = db₂dx₂(t,q)

    D²ϑ[3,1] = d²A₂dx₃dx₁(t,q) + u(t,q) * d²b₂dx₃dx₁(t,q)
    D²ϑ[3,2] = d²A₂dx₃dx₂(t,q) + u(t,q) * d²b₂dx₃dx₂(t,q)
    D²ϑ[3,3] = d²A₂dx₃dx₃(t,q) + u(t,q) * d²b₂dx₃dx₃(t,q)
    D²ϑ[3,4] = db₂dx₃(t,q)

    D²ϑ[4,1] = 0
    D²ϑ[4,2] = 0
    D²ϑ[4,3] = 0
    D²ϑ[4,4] = 0

    nothing
end

function D²ϑd₃(t, q, D²ϑ)
    D²ϑ[2,1] = d²A₁dx₃dx₁(t,q) + u(t,q) * d²b₁dx₃dx₁(t,q)
    D²ϑ[2,2] = d²A₁dx₃dx₂(t,q) + u(t,q) * d²b₁dx₃dx₂(t,q)
    D²ϑ[2,3] = d²A₁dx₃dx₃(t,q) + u(t,q) * d²b₁dx₃dx₃(t,q)
    D²ϑ[2,4] = db₁dx₃(t,q)

    D²ϑ[2,1] = d²A₃dx₂dx₁(t,q) + u(t,q) * d²b₃dx₂dx₁(t,q)
    D²ϑ[2,2] = d²A₃dx₂dx₂(t,q) + u(t,q) * d²b₃dx₂dx₂(t,q)
    D²ϑ[2,3] = d²A₃dx₂dx₃(t,q) + u(t,q) * d²b₃dx₂dx₃(t,q)
    D²ϑ[2,4] = db₃dx₂(t,q)

    D²ϑ[3,1] = d²A₃dx₃dx₁(t,q) + u(t,q) * d²b₃dx₃dx₁(t,q)
    D²ϑ[3,2] = d²A₃dx₃dx₂(t,q) + u(t,q) * d²b₃dx₃dx₂(t,q)
    D²ϑ[3,3] = d²A₃dx₃dx₃(t,q) + u(t,q) * d²b₃dx₃dx₃(t,q)
    D²ϑ[3,4] = db₃dx₃(t,q)

    D²ϑ[4,1] = 0
    D²ϑ[4,2] = 0
    D²ϑ[4,3] = 0
    D²ϑ[4,4] = 0

    nothing
end

function D²ϑd₄(t, q, D²ϑ)
    D²ϑ[1,1] = db₁dx₁(t,q)
    D²ϑ[1,2] = db₁dx₂(t,q)
    D²ϑ[1,3] = db₁dx₃(t,q)
    D²ϑ[1,4] = 0

    D²ϑ[2,1] = db₂dx₁(t,q)
    D²ϑ[2,2] = db₂dx₂(t,q)
    D²ϑ[2,3] = db₂dx₃(t,q)
    D²ϑ[2,4] = 0

    D²ϑ[3,1] = db₃dx₁(t,q)
    D²ϑ[3,2] = db₃dx₂(t,q)
    D²ϑ[3,3] = db₃dx₃(t,q)
    D²ϑ[3,4] = 0

    D²ϑ[4,1] = 0
    D²ϑ[4,2] = 0
    D²ϑ[4,3] = 0
    D²ϑ[4,4] = 0

    nothing
end

function dϑ(t, q, dϑ)
    dϑ[1,1] = dϑ₁dx₁(t,q)
    dϑ[1,2] = dϑ₁dx₂(t,q)
    dϑ[1,3] = dϑ₁dx₃(t,q)
    dϑ[1,4] = dϑ₁dx₄(t,q)

    dϑ[2,1] = dϑ₂dx₁(t,q)
    dϑ[2,2] = dϑ₂dx₂(t,q)
    dϑ[2,3] = dϑ₂dx₃(t,q)
    dϑ[2,4] = dϑ₂dx₄(t,q)

    dϑ[3,1] = dϑ₃dx₁(t,q)
    dϑ[3,2] = dϑ₃dx₂(t,q)
    dϑ[3,3] = dϑ₃dx₃(t,q)
    dϑ[3,4] = dϑ₃dx₄(t,q)

    dϑ[4,1] = dϑ₄dx₁(t,q)
    dϑ[4,2] = dϑ₄dx₂(t,q)
    dϑ[4,3] = dϑ₄dx₃(t,q)
    dϑ[4,4] = dϑ₄dx₄(t,q)

    nothing
end

β₁(t,q) = dϑ₃dx₂(t,q) - dϑ₂dx₃(t,q)
β₂(t,q) = dϑ₁dx₃(t,q) - dϑ₃dx₁(t,q)
β₃(t,q) = dϑ₂dx₁(t,q) - dϑ₁dx₂(t,q)

# function β(t,q)
#    return sqrt(β1(t,q)^2 + β2(t,q)^2 + β3(t,q)^2)
# end

function hamiltonian(t, q, params)
    @unpack μ = params
    0.5 * u(t,q)^2 + μ*B(t,q)
end

dHdx₁(t, q, μ) = μ * dBdx₁(t,q)
dHdx₂(t, q, μ) = μ * dBdx₂(t,q)
dHdx₃(t, q, μ) = μ * dBdx₃(t,q)
dHdx₄(t, q, μ) = u(t,q)

function dH(t, q, dH, params)
    @unpack μ = params
    dH[1] = dHdx₁(t, q, μ)
    dH[2] = dHdx₂(t, q, μ)
    dH[3] = dHdx₃(t, q, μ)
    dH[4] = dHdx₄(t, q, μ)
    nothing
end

f₁(t,q,v) = dϑ₁dx₁(t,q) * v[1] + dϑ₂dx₁(t,q) * v[2] + dϑ₃dx₁(t,q) * v[3] + dϑ₄dx₁(t,q) * v[4]
f₂(t,q,v) = dϑ₁dx₂(t,q) * v[1] + dϑ₂dx₂(t,q) * v[2] + dϑ₃dx₂(t,q) * v[3] + dϑ₄dx₂(t,q) * v[4]
f₃(t,q,v) = dϑ₁dx₃(t,q) * v[1] + dϑ₂dx₃(t,q) * v[2] + dϑ₃dx₃(t,q) * v[3] + dϑ₄dx₃(t,q) * v[4]
f₄(t,q,v) = dϑ₁dx₄(t,q) * v[1] + dϑ₂dx₄(t,q) * v[2] + dϑ₃dx₄(t,q) * v[3] + dϑ₄dx₄(t,q) * v[4]

g₁(t,q,v) = dϑ₁dx₁(t,q) * v[1] + dϑ₁dx₂(t,q) * v[2] + dϑ₁dx₃(t,q) * v[3] + dϑ₁dx₄(t,q) * v[4]
g₂(t,q,v) = dϑ₂dx₁(t,q) * v[1] + dϑ₂dx₂(t,q) * v[2] + dϑ₂dx₃(t,q) * v[3] + dϑ₂dx₄(t,q) * v[4]
g₃(t,q,v) = dϑ₃dx₁(t,q) * v[1] + dϑ₃dx₂(t,q) * v[2] + dϑ₃dx₃(t,q) * v[3] + dϑ₃dx₄(t,q) * v[4]
g₄(t,q,v) = dϑ₄dx₁(t,q) * v[1] + dϑ₄dx₂(t,q) * v[2] + dϑ₄dx₃(t,q) * v[3] + dϑ₄dx₄(t,q) * v[4]

function g̅₁(t, q, v)
    D²ϑ = zeros(eltype(q), length(q), length(v))
    D²ϑd₁(t, q, D²ϑ)
    return transpose(q) * (D²ϑ * v)
end

function g̅₂(t, q, v)
    D²ϑ = zeros(eltype(q), length(q), length(v))
    D²ϑd₂(t, q, D²ϑ)
    return transpose(q) * (D²ϑ * v)
end

function g̅₃(t, q, v)
    D²ϑ = zeros(eltype(q), length(q), length(v))
    D²ϑd₃(t, q, D²ϑ)
    return transpose(q) * (D²ϑ * v)
end

function g̅₄(t, q, v)
    D²ϑ = zeros(eltype(q), length(q), length(v))
    D²ϑd₄(t, q, D²ϑ)
    return transpose(q) * (D²ϑ * v)
end

function guiding_center_4d_v(t, q::Vector{DT}, v::Vector{DT}, params) where {DT}
    @unpack μ = params

    local lB₁ = B₁(t,q)
    local lB₂ = B₂(t,q)
    local lB₃ = B₃(t,q)

    local lβ₁ = β₁(t,q)
    local lβ₂ = β₂(t,q)
    local lβ₃ = β₃(t,q)
    local lβ  = lβ₁ * dϑ₁dx₄(t,q) + lβ₂ * dϑ₂dx₄(t,q) + lβ₃ * dϑ₃dx₄(t,q)

    local ∇₁B = dBdx₁(t,q)
    local ∇₂B = dBdx₂(t,q)
    local ∇₃B = dBdx₃(t,q)

    v[1] = ( u(t,q) * lβ₁ - μ * ( ∇₂B * dϑ₃dx₄(t,q) - ∇₃B * dϑ₂dx₄(t,q) ) ) / lβ
    v[2] = ( u(t,q) * lβ₂ - μ * ( ∇₃B * dϑ₁dx₄(t,q) - ∇₁B * dϑ₃dx₄(t,q) ) ) / lβ
    v[3] = ( u(t,q) * lβ₃ - μ * ( ∇₁B * dϑ₂dx₄(t,q) - ∇₂B * dϑ₁dx₄(t,q) ) ) / lβ
    v[4] = - μ * ( ∇₁B * lβ₁ + ∇₂B * lβ₂ + ∇₃B * lβ₃ ) / lβ

    nothing
end

function guiding_center_4d_v(t, q::Vector{DT}, p::Vector{DT}, v::Vector{DT}, params) where {DT}
    guiding_center_4d_v(t, q, v, params)
end

# function guiding_center_4d_ϑ(t, q::Vector{DT}, Θ::Vector{DT}) where {DT}
#     ϑ(t, q, Θ)
# end

function guiding_center_4d_ϑ(t, q::Vector{DT}, v::Vector{DT}, Θ::Vector{DT}, params) where {DT}
    ϑ(t, q, Θ)
end

function guiding_center_4d_ϑ(t, q::Vector{DT}, v::Vector{DT}, Θ::Vector{DT}, params, κ) where {DT}
    Θ[1] = (1-κ) * ϑ₁(t,q) - κ * f₁(t,q,q)
    Θ[2] = (1-κ) * ϑ₂(t,q) - κ * f₂(t,q,q)
    Θ[3] = (1-κ) * ϑ₃(t,q) - κ * f₃(t,q,q)
    Θ[4] = (1-κ) * ϑ₄(t,q) - κ * f₄(t,q,q)
    nothing
end

function guiding_center_4d_f(t, q::Vector{DT}, v::Vector{DT}, f::Vector{DT}, params) where {DT}
    @unpack μ = params
    f[1] = f₁(t,q,v) - dHdx₁(t,q,μ)
    f[2] = f₂(t,q,v) - dHdx₂(t,q,μ)
    f[3] = f₃(t,q,v) - dHdx₃(t,q,μ)
    f[4] = f₄(t,q,v) - dHdx₄(t,q,μ)
    nothing
end

function guiding_center_4d_f(t, q::Vector{DT}, v::Vector{DT}, f::Vector{DT}, params, κ) where {DT}
    @unpack μ = params
    f[1] = (1-κ) * f₁(t,q,v) - κ * (g₁(t,q,v) + g̅₁(t,q,v)) - dHdx₁(t,q,μ)
    f[2] = (1-κ) * f₂(t,q,v) - κ * (g₂(t,q,v) + g̅₂(t,q,v)) - dHdx₂(t,q,μ)
    f[3] = (1-κ) * f₃(t,q,v) - κ * (g₃(t,q,v) + g̅₃(t,q,v)) - dHdx₃(t,q,μ)
    f[4] = (1-κ) * f₄(t,q,v) - κ * (g₄(t,q,v) + g̅₄(t,q,v)) - dHdx₄(t,q,μ)
    nothing
end

function guiding_center_4d_g(t::Number, q::Vector, λ::Vector, g::Vector, params)
    g[1] = f₁(t,q,λ)
    g[2] = f₂(t,q,λ)
    g[3] = f₃(t,q,λ)
    g[4] = f₄(t,q,λ)
    nothing
end

# function guiding_center_4d_g(t, q::Vector{DT}, λ::Vector{DT}, g::Vector{DT}, params, κ::DT=0) where {DT}
#     g[1] = g₁(t,q,λ)
#     g[2] = g₂(t,q,λ)
#     g[3] = g₃(t,q,λ)
#     g[4] = g₄(t,q,λ)
#     nothing
# end

function guiding_center_4d_g(t, q::Vector{DT}, λ::Vector{DT}, g::Vector{DT}, params, κ::DT) where {DT}
    g[1] = (1-κ) * f₁(t,q,λ) - κ * (g₁(t,q,λ) + g̅₁(t,q,λ))
    g[2] = (1-κ) * f₂(t,q,λ) - κ * (g₂(t,q,λ) + g̅₂(t,q,λ))
    g[3] = (1-κ) * f₃(t,q,λ) - κ * (g₃(t,q,λ) + g̅₃(t,q,λ))
    g[4] = (1-κ) * f₄(t,q,λ) - κ * (g₄(t,q,λ) + g̅₄(t,q,λ))
    nothing
end

# function guiding_center_4d_g(t, q::Vector{DT}, λ::Vector{DT}, g::Vector{DT}, params, κ::DT=0) where {DT}
#     g[1] = (1-κ) * g₁(t,q,λ) - κ * (f₁(t,q,λ) + g̅₁(t,q,λ))
#     g[2] = (1-κ) * g₂(t,q,λ) - κ * (f₂(t,q,λ) + g̅₂(t,q,λ))
#     g[3] = (1-κ) * g₃(t,q,λ) - κ * (f₃(t,q,λ) + g̅₃(t,q,λ))
#     g[4] = (1-κ) * g₄(t,q,λ) - κ * (f₄(t,q,λ) + g̅₄(t,q,λ))
#     nothing
# end

function guiding_center_4d_dH(t::DT, q::Vector{DT}, ∇H::Vector{DT}, params) where {DT}
    dH(t, q, ∇H, params)
end

function guiding_center_4d_ω(t::DT, q::Vector{DT}, Ω::Matrix{DT}, params) where {DT}
    ω(t, q, Ω)
end

function guiding_center_4d_λ(t::DT, q::Vector{DT}, μ, λ::Vector{DT}, Ω::Matrix{DT}, dh::Vector{DT}) where {DT}
    dH(t, q, μ, dh)
    ω(t, q, Ω)
    λ .= inv(Ω) * dh
    nothing
end

function guiding_center_4d_λ(t::DT, q::Vector{DT}, μ, λ::Vector{DT}) where {DT}
    D  = size(q, 1)
    guiding_center_4d_λ(t, q, μ, λ, zeros(DT,D,D), zeros(DT,D))
end

function guiding_center_4d_pᵢ(qᵢ, tᵢ=0)
    pᵢ = zero(qᵢ)

    if ndims(qᵢ) == 1
        ϑ(tᵢ, qᵢ, pᵢ)
    else
        for i in 1:size(qᵢ,2)
            ϑ(tᵢ, (@view qᵢ[:,i]), (@view pᵢ[:,i]))
        end
    end
    pᵢ
end

function guiding_center_4d_λᵢ(qᵢ::Array{DT}, μ, Δt::DT=1, tᵢ::DT=0) where {DT}
    D  = size(qᵢ, 1)
    λᵢ = zero(qᵢ)
    Ω  = zeros(DT, D, D)
    dh = zeros(DT, D)

    if ndims(qᵢ) == 1
        guiding_center_4d_λ(qᵢ, μ, λᵢ, Ω, dh, Δt, tᵢ)
    else
        for i in 1:size(qᵢ,2)
            guiding_center_4d_λ((@view qᵢ[:,i]), μ,  (@view λᵢ[:,i]), Ω, dh, Δt, tᵢ)
        end
    end

    return λᵢ
end


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

struct StaticPool <: AbstractThreadPool
    tids :: Vector{Int}
end

#############################
# Constructors
#############################

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

The main StaticPool object.
"""
function StaticPool(init_thrd::Integer=1, nthrds::Integer=Threads.nthreads())
    thrd0 = min(init_thrd, Threads.nthreads())
    thrd1 = min(init_thrd+nthrds-1, Threads.nthreads())
    return StaticPool(thrd0:thrd1)
end

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

function tmap(fn, pool::StaticPool, itr)
    data = collect(itr)
    applicable(fn, first(data)) || error("function can't be applied to iterator contents")
    N = length(data)
    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
            for i in start(ind):finish(ind)
                @inbounds result[i] = fn(Base.unsafe_getindex(data, i))
            end
        end
    end

    tasks = Channel(Inf)
    for tid in 1:Threads.nthreads()
        ind = findfirst(t->t==tid, pool.tids)
        if !isnothing(ind)
            task = @tspawnat tid _fn(ind)
            put!(tasks, task)
        end
    end
    _sync_end(tasks)

    return result
end

# From Julia PR #39007
function _schedule_result(t, rs)
    schedule(Task(() -> begin
        ex = nothing
        try
            wait(t)
        catch e
            ex = e
        finally
            put!(rs, ex)
        end
    end))
end

# Modified from Julia PR #39007 to duplicate error handling of tmap
function _sync_end(c::Channel{Any})
    try
        n = 0
        isready(c) || return
        while true
            t = take!(c)
            if t isa Exception
                while isready(c)                       # Clean the channel
                    take!(c)
                end
                throw(t)
            elseif isnothing(t)                        # Successful result from monitor.
                n -= 1                                 # Leave if nothing remains.
                n == 0 && !isready(c) && break
            else                                       # Not a monitor result - must be
                n += 1                                 # a new waitable. Create a new
                _schedule_result(t, c)                 # monitor and return to Channel
            end
        end
    finally
        close(c)
    end
    nothing
end


text:
const libllvmgeneral = find_library(["libllvmgeneral",], ["/home/jake/.julia/v0.3/LLVM/deps/usr/lib"])
@assert(libllvmgeneral != "", """Failed to find required library libllvmgeneral.
                            Try re-running the package script using Pkg.build("LLVM")""")
const libllvm = find_library(["libLLVM-3.6.0svn"])
@assert(libllvm != "", "Failed to find llvm library libLLVM-3.6.0svn")


text:
using Documenter
using Skewlia

makedocs(
    sitename = "Skewlia",
    format = Documenter.HTML(),
    modules = [Skewlia]
)

deploydocs(
    repo = "github.com/sysread/Skewlia.git",
)

max_stars_repo_name: sysread/Skewlia

max_issues_repo_name: sysread/Skewlia

max_forks_repo_name: sysread/Skewlia


text:
# Faster than Polynomials.jl with less memory allocation + single-point type
# stability; adapted from Octave and
# https://github.com/JuliaMath/Polynomials.jl/blob/master/src/Polynomials.jl
function poly(x::Array{T,1}) where T <: Number
  n = length(x)
  y = zeros(T, n+1)
  y[1] = one(T)
  for j = 1:n
    y[2:j+1] .-= x[j].*y[1:j]
  end
  return y
end

function polyval(p::Array{T1,1}, x::T2) where {T1 <: Number, T2 <: Number}
  y = T2(p[1])
  for i = 2:lastindex(p)
    y = p[i] .+ x*y
  end
  return y
end

function polyval(p::Array{T1,1}, x::Array{T2,1}) where {T1 <: Number, T2 <: Number}
  y = ones(T2, length(x)) .* p[1]
  for i = 2:length(p)
    broadcast!(*, y, y, x)
    broadcast!(+, y, y, p[i])
  end
  return y
end

function polyfit(x::Array{T1,1}, y::Array{T2,1}, n::Integer=1) where {T1 <: Real, T2 <: Real}
  nx = length(x)
  nx == length(y) || error("SeisIO.polyfit requires length(t) == length(x)")
  -1 < n < nx || throw(DomainError)
  A = Array{T2, 2}(undef, length(x), n+1)
  A[:,n+1] .= one(T2)
  for i = n:-1:1
      A[:,i] .= A[:,i+1] .* x
  end
  return A \ y
end

# Convert to Float64 or use improved sum
function linreg(t::Array{Float64,1}, x::AbstractArray{T,1}) where T
  n = length(t)
  st = sum(t)
  sx = sum(x)
  stt = dot(t,t)
  stx = dot(t,x)
  sxx = dot(x,x)
  d = n*stt - st*st
  b = (stt*sx - st*stx)/d
  a = (n*stx - st*sx)/d
  return T[a,b]
end

# p = linreg(x, fs)
function linreg(x::AbstractArray{T,1}, dt::Float64) where T
  n = length(x)
  t = (1:n)*dt
  st = sum(t)
  sx = sum(x)
  stt = dot(t,t)
  stx = dot(t,x)
  sxx = dot(x,x)
  d = n*stt - st*st
  b = (stt*sx - st*stx)/d
  a = (n*stx - st*sx)/d
  return T[a,b]
end

#= CHANGELOG for this file
2019-08-08
* poly, polyval, and polyfit should now always output powers in descending
order, i.e., p^n ... p^0
* BUG fixed: polyval(p::Array{T,1}, x::T) used power ordering of Polynomials.jl;
corrected to be consistent with other routines
* polyfit now allows order n=0 and takes any Integer for n; n=0 returns the mean
2019-08-19
* added linreg for low-memory linear regression; identical to SAC detrend
2019-09-03
* bug fix for rare situation where float precision led to length(t) != length(x)
=#


text:
using Test
using Schemata

using CSV
using DataFrames
using Dates
using TimeZones

const CATEGORICAL = true
const REQUIRED    = true
const UNIQUE      = true

include("test_inmemory_tables.jl")
include("test_ondisk_tables.jl")


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -37.14009556780016 <= q <= 48.365016418887166)

                     add_NL_constraint(m, :(gelu(0.24902435571631587 + -0.02312729681609138*gelu(-0.46870412880299783 + 0.1700038462852933*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.44738483682757924*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.5112184633122538*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.17714072857386487*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9529511507705823*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.5418072450548097*gelu(0.5019192867350393 + 0.2997113620022893*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.5721484689942011*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.4073624372378828*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + -0.07630883358086304*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.3103910161241421*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.9789829058919253*gelu(-0.652369829242545 + -0.998242239199818*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.041041478696029365*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.4348093907280841*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.7773011179238121*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7537599741964742*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.13446361488468206*gelu(-0.40141742757451926 + 0.7289092837988314*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + 0.16077592157977172*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.6378484913808045*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9326240829024273*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9638370498324353*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.3069132295129271*gelu(0.17680227033398044 + 0.8293506390949057*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.09356004658531791*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.16358082594109513*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9029174038261001*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7318887196209252*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3])))) + gelu(-0.017552134571214406 + -0.41004797399293524*gelu(-0.46870412880299783 + 0.1700038462852933*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.44738483682757924*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.5112184633122538*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.17714072857386487*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9529511507705823*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.47019448283071075*gelu(0.5019192867350393 + 0.2997113620022893*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.5721484689942011*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.4073624372378828*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + -0.07630883358086304*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.3103910161241421*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.35795867979675133*gelu(-0.652369829242545 + -0.998242239199818*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.041041478696029365*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.4348093907280841*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.7773011179238121*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7537599741964742*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.44715886126593407*gelu(-0.40141742757451926 + 0.7289092837988314*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + 0.16077592157977172*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.6378484913808045*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9326240829024273*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9638370498324353*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.22638221307132733*gelu(0.17680227033398044 + 0.8293506390949057*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.09356004658531791*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.16358082594109513*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9029174038261001*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7318887196209252*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3])))) + gelu(-0.4029592795849051 + -0.6197606146665708*gelu(-0.46870412880299783 + 0.1700038462852933*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.44738483682757924*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.5112184633122538*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.17714072857386487*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9529511507705823*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.8464388391219866*gelu(0.5019192867350393 + 0.2997113620022893*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.5721484689942011*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.4073624372378828*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + -0.07630883358086304*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.3103910161241421*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.6262663939568953*gelu(-0.652369829242545 + -0.998242239199818*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.041041478696029365*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.4348093907280841*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.7773011179238121*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7537599741964742*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.06090941414894013*gelu(-0.40141742757451926 + 0.7289092837988314*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + 0.16077592157977172*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.6378484913808045*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9326240829024273*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9638370498324353*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.6699649674430836*gelu(0.17680227033398044 + 0.8293506390949057*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.09356004658531791*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.16358082594109513*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9029174038261001*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7318887196209252*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3])))) + gelu(0.6579858782021573 + -0.21911632564812766*gelu(-0.46870412880299783 + 0.1700038462852933*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.44738483682757924*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.5112184633122538*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.17714072857386487*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9529511507705823*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.7292607357990395*gelu(0.5019192867350393 + 0.2997113620022893*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.5721484689942011*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.4073624372378828*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + -0.07630883358086304*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.3103910161241421*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.23882038501163416*gelu(-0.652369829242545 + -0.998242239199818*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.041041478696029365*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.4348093907280841*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.7773011179238121*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7537599741964742*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.6855544316424598*gelu(-0.40141742757451926 + 0.7289092837988314*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + 0.16077592157977172*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.6378484913808045*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9326240829024273*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9638370498324353*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.2868272308560216*gelu(0.17680227033398044 + 0.8293506390949057*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.09356004658531791*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.16358082594109513*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9029174038261001*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7318887196209252*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3])))) + gelu(-0.10147744083820065 + -0.7175110849623545*gelu(-0.46870412880299783 + 0.1700038462852933*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.44738483682757924*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.5112184633122538*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.17714072857386487*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9529511507705823*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.9583278917345623*gelu(0.5019192867350393 + 0.2997113620022893*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.5721484689942011*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.4073624372378828*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + -0.07630883358086304*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.3103910161241421*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.7313133539296679*gelu(-0.652369829242545 + -0.998242239199818*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.041041478696029365*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.4348093907280841*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.7773011179238121*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7537599741964742*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + 0.008349755673621395*gelu(-0.40141742757451926 + 0.7289092837988314*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + 0.16077592157977172*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + 0.6378484913808045*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9326240829024273*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.9638370498324353*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3]))) + -0.13663535205769417*gelu(0.17680227033398044 + 0.8293506390949057*gelu(-0.07115429328131828 + -0.06282459839957566*$(x[1]) + 0.4662931752239805*$(x[2]) + -0.47993632100881634*$(x[3])) + -0.09356004658531791*gelu(-0.42367185224307535 + 0.829552620135646*$(x[1]) + -0.7267035823407642*$(x[2]) + 0.6637528273154039*$(x[3])) + -0.16358082594109513*gelu(0.6408115001973855 + -0.9252249601617657*$(x[1]) + 0.6902522661180681*$(x[2]) + 0.6531628217744303*$(x[3])) + 0.9029174038261001*gelu(0.3778204937005931 + 0.21019890124389518*$(x[1]) + -0.5904566970281198*$(x[2]) + -0.06968695857557394*$(x[3])) + 0.7318887196209252*gelu(0.7716365722893217 + 0.12636272698834716*$(x[1]) + 0.15699019471654063*$(x[2]) + -0.5007466869555728*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Indexing
# ========
#
# Indexing methods for biological sequences.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# Getting and setting elements in a biological sequence.

function Base.checkbounds(seq::BioSequence, range::UnitRange)
    if 1 ≤ range.start && range.stop ≤ endof(seq)
        return true
    end
    throw(BoundsError(seq, range))
end

function Base.checkbounds(seq::BioSequence, locs::AbstractVector{Bool})
    if length(seq) == length(locs)
        return true
    end
    throw(BoundsError(seq, locs))
end

function Base.checkbounds(seq::BioSequence, locs::AbstractVector)
    for i in locs
        checkbounds(seq, i)
    end
    return true
end

function checkdimension(from::Integer, to::Integer)
    if from == to
        return true
    end
    throw(DimensionMismatch(string(
        "attempt to assign ",
        from, " elements to ",
        to,   " elements")))
end

function checkdimension(seq::BioSequence, locs::AbstractVector)
    return checkdimension(length(seq), length(locs))
end

function checkdimension(seq::BioSequence, locs::AbstractVector{Bool})
    return checkdimension(length(seq), sum(locs))
end

# creates a bit mask for given number of bits `n`
mask(n::Integer) = (UInt64(1) << n) - 1
mask{A<:Alphabet}(::Type{A}) = mask(bitsof(A))

# assumes `i` is positive and `bitsof(A)` is a power of 2
@inline function bitindex{A}(seq::BioSequence{A}, i::Integer)
    return BitIndex((i + first(seq.part) - 2) << trailing_zeros(bitsof(A)))
end

@inline function inbounds_getindex{A}(seq::BioSequence{A}, i::Integer)
    j = bitindex(seq, i)
    @inbounds return decode(A, (seq.data[index(j)] >> offset(j)) & mask(A))
end

Base.getindex(seq::BioSequence, part::UnitRange) = BioSequence(seq, part)
Base.view(seq::BioSequence, part::UnitRange) = BioSequence(seq, part)

function Base.setindex!{A,T<:Integer}(seq::BioSequence{A},
                                      other::BioSequence{A},
                                      locs::AbstractVector{T})
    checkbounds(seq, locs)
    checkdimension(other, locs)
    orphan!(seq)
    for (i, x) in zip(locs, other)
        unsafe_setindex!(seq, x, i)
    end
    return seq
end

function Base.setindex!{A,T<:Integer}(seq::BioSequence{A},
                                      other::BioSequence{A},
                                      locs::UnitRange{T})
    checkbounds(seq, locs)
    checkdimension(other, locs)
    return copy!(seq, locs.start, other, 1)
end

function Base.setindex!{A}(seq::BioSequence{A},
                           other::BioSequence{A},
                           locs::AbstractVector{Bool})
    checkbounds(seq, locs)
    checkdimension(other, locs)
    orphan!(seq)
    i = j = 0
    while (i = findnext(locs, i + 1)) > 0
        unsafe_setindex!(seq, other[j+=1], i)
    end
    return seq
end

function Base.setindex!{A}(seq::BioSequence{A}, other::BioSequence{A}, ::Colon)
    return setindex!(seq, other, 1:endof(seq))
end

function Base.setindex!(seq::BioSequence, x, i::Integer)
    checkbounds(seq, i)
    orphan!(seq)
    return unsafe_setindex!(seq, x, i)
end

function Base.setindex!{A,T<:Integer}(seq::BioSequence{A}, x, locs::AbstractVector{T})
    checkbounds(seq, locs)
    bin = enc64(seq, x)
    orphan!(seq)
    for i in locs
        encoded_setindex!(seq, bin, i)
    end
    return seq
end

function Base.setindex!{A}(seq::BioSequence{A}, x, locs::AbstractVector{Bool})
    checkbounds(seq, locs)
    bin = enc64(seq, x)
    orphan!(seq)
    i = j = 0
    while (i = findnext(locs, i + 1)) > 0
        encoded_setindex!(seq, bin, i)
    end
    return seq
end

Base.setindex!{A}(seq::BioSequence{A}, x, ::Colon) = setindex!(seq, x, 1:endof(seq))

# this is "unsafe" because of no bounds check and no orphan! call
@inline function unsafe_setindex!{A}(seq::BioSequence{A}, x, i::Integer)
    bin = enc64(seq, x)
    return encoded_setindex!(seq, bin, i)
end

@inline function encoded_setindex!{A}(seq::BioSequence{A}, bin::UInt64, i::Integer)
    j, r = bitindex(seq, i)
    @inbounds seq.data[j] = (bin << r) | (seq.data[j] & ~(mask(A) << r))
    return seq
end


text:
using ElectricFields
using Unitful

using PyPlot
using Jagot.plotting
plot_style("ggplot")

function savedocfig(name,dir="figures")
    fig = gcf()
    filename = joinpath(@__DIR__, "src", dir, "$(name).svg")
    savefig(filename,
            transparent=false,
            facecolor=fig.get_facecolor())
    close(fig)
    if isfile(filename)
        println("Saved $(name) to $(filename)")
    else
        @warn "Saving $(name) to $(filename) failed"
    end
end

function savedfigure(fun::Function, name, args...; kwargs...)
    cfigure(name, args...; kwargs...) do
        fun()
    end
    savedocfig(name)
end

function index_example()
    @field(IR) do
        I₀ = 1e14u"W/cm^2"
        λ = 800.0u"nm"
        τ = 6.2u"fs"
        σmax = 6.0
    end
    t = timeaxis(IR)
    tplot = 24.3e-3t
    Fv = field_amplitude(IR, t)
    Fenv = field_envelope.(IR, t)
    Av = vector_potential(IR, t)
    Iv = instantaneous_intensity.(IR, t)
    Ienv = intensity.(IR, t)

    savedfigure("index_example") do
        csubplot(311, nox=true) do
            plot(tplot, Av)
            ylabel(L"$A(t)$ [au]")
        end
        csubplot(312, nox=true) do
            plot(tplot, Fv)
            plot(tplot, Fenv)
            ylabel(L"$F(t)$ [au]")
        end
        csubplot(313) do
            plot(tplot, Iv)
            plot(tplot, Ienv)
            ylabel(L"$I(t)$ [au]")
        end
        xlabel(L"$t$ [fs]")
    end
end

function index_polarized_example()
    @field(A) do
        I₀ = 1.0
        ω = 1.0
        cycles = 6.0
        env = :cos²
        ξ = 1.0
    end

    @field(B) do
        I₀ = 1.0
        ω = 2.0
        cycles = 6.0
        env = :cos²
        ξ = -1.0
    end

    F = A + delay(B, 3/2π)

    t = timeaxis(F)
    Fv = field_amplitude(F, t)

    savedfigure("index_polarized_example", figsize=(8,10)) do
        csubplot(211) do
            plot(t, Fv[:,1], label=L"F_x(t)")
            plot(t, Fv[:,3], label=L"F_z(t)")
            legend()
            xlabel(L"$t$ [jiffies]")
        end
        csubplot(212, projection="3d") do
            plot3D(t, Fv[:,1], Fv[:,3])
            xlabel(L"$t$ [jiffies]")
            ylabel(L"F_x(t)")
            gca().set_zlabel(L"F_z(t)")
        end
    end
end

macro echo(expr)
    println(expr)
    :(@time $expr)
end

@info "Documentation plots"
mkpath("docs/src/figures")
@echo index_example()
@echo index_polarized_example()


text:
module GenomeLoader

using BioSequences, GenomicFeatures, FASTX

export load_genome, get_seq

struct GenomeRecord{T}
	name::String
	metadata::String
	seq::T
end

function convert_N_to_gap!(seq)
    for i = 1:length(seq)
        (seq[i] == DNA_N) && (seq[i] = DNA_Gap)
    end
end

function load_genome(fastafile; verbose=true)
	indexfile = fastafile*".fai"
	!isfile(indexfile) && error("Index: $indexfile not found")
	reader = FASTA.Reader(open(fastafile, "r"), index=indexfile)

    genome = Dict{String, GenomeRecord{LongSequence{DNAAlphabet{4}}}}()

    for record in reader
        id = FASTA.identifier(record)
        verbose && println("Loading $id...")
        seq = FASTA.sequence(record)
        convert_N_to_gap!(seq)
        genome[id] = GenomeRecord(id, "", seq)
    end
	genome
end

@inline function get_seq(genome, chr, range)
    ((range.start < 1) || (range.stop .> length(genome[chr].seq))) && return genome[chr].seq[1:0]
    genome[chr].seq[range]
end

@inline function get_seq(genome, chr, range, strand::String)
    seq = get_seq(genome, chr, range)
    (strand == "-") && reverse_complement!(seq)
    seq
end

@inline function get_seq(genome, chr, range, strand::Strand)
    seq = get_seq(genome, chr, range)
    (strand == STRAND_NEG) && reverse_complement!(seq)
    seq
end

end # module


text:
# uses variables from setup_tests.jl

STDNORMAL = Normal()
G(x::Real) = cdf(Normal(0, σ), x)
μ(a::Real) = σ * (pdf(STDNORMAL, quantile(STDNORMAL, 1-a)) - a * quantile(STDNORMAL, 1-a))

### Verifications: endogenous divorce model ###

"""
Element-by-element calculation of steady state flow equations
```julia
mres = ℓ_m - um .* (1 .+ λ .* ((α ./ (δ * (1 .- α) .+ ψ_m .+ ψ_f')) * uf))
fres = ℓ_f - uf .* (1 .+ λ .* ((α ./ (δ * (1 .- α) .+ ψ_m .+ ψ_f'))' * um))
```
"""
function sse_stoch(M::SearchMatch)

	mres = similar(M.ℓ_m)
	fres = similar(M.ℓ_f)

	for i in 1:length(M.ℓ_m)
		mres[i] = M.ℓ_m[i] - M.u_m[i] * (1 + M.λ * dot(M.α[i,:] ./
		          (M.δ * (1 .- M.α[i,:]) .+ M.ψ_m[i] .+ M.ψ_f), M.u_f))
	end
	for j in 1:length(M.ℓ_f)
		fres[j] = M.ℓ_f[j] - M.u_f[j] * (1 + M.λ * dot(M.α[:,j] ./
		          (M.δ * (1 .- M.α[:,j]) .+ M.ψ_f[j] .+ M.ψ_m), M.u_m))
	end
	return mres, fres
end

"""
Element-by-element calculation of value function equations
```julia

v_m = 0.5*λ * (μα * u_f)
v_f = 0.5*λ * (μα' * u_m)
```
where `μα = μ.(α) ./ (r + δ + ψ_m .+ ψ_f')`
"""
function vf_stoch(M::SearchMatch)
	mres = similar(M.v_m)
	fres = similar(M.v_f)

	for i in 1:length(M.v_m)
		mres[i] = 2*M.v_m[i] - M.λ * dot(μ.(M.α[i,:]) ./
		                             (M.r + M.δ + M.ψ_m[i] .+ M.ψ_f), M.u_f)
	end
	for j in 1:length(M.v_f)
		fres[j] = 2*M.v_f[j] - M.λ * dot(μ.(M.α[:,j]) ./
		                             (M.r + M.δ + M.ψ_f[j] .+ M.ψ_m), M.u_m)
	end

	return mres, fres
end

function surplus_stoch(M::SearchMatch)
	s = similar(M.h)

	for i in 1:length(M.v_m), j in 1:length(M.v_f)
		s[i,j] = M.h[i,j] - M.v_m[i] - M.v_f[j] + M.δ * μ(M.α[i,j]) /
		                               (M.r + M.δ + M.ψ_m[i] .+ M.ψ_f[j])
	end

	return s
end

### Tests: endogenous divorce model ###

# symmetric case

@test rsym.v_m ≈ rsym.v_f #expected symmetry of values
@test rsym.u_m ≈ rsym.u_f #expected symmetry of singles
@test rsym.α ≈ rsym.α' #expected symmetry of matching

# asymmetric case

# check convergence
rmsse, rfsse = sse_stoch(rasym)
rmvf, rfvf = vf_stoch(rasym)

# valid solution
@test rmsse ≈ zero(rmsse) atol=1e-7 #market equilibrium: single men did not converge
@test rfsse ≈ zero(rfsse) atol=1e-7 #market equilibrium: single women did not converge
# convergence seems to stall around 2e-5...
@test maximum(abs, rmvf) ≈ 0 atol=3e-5 #matching equilibrium: man values did not converge
@test maximum(abs, rfvf) ≈ 0 atol=1e-7 #matching equilibrium: woman values did not converge
@test rasym.α ≈ 1 .- G.(-surplus_stoch(rasym)) atol=1e-5 #matching equilibrium: match probabilities α did not converge

# sex ratio effects on singles
@test all(rsym.u_f .≤ rasym.u_f) #expected more singlehood with fewer available men


text:
using RenderingGeometry

"""

    2D Component Geometry Definitions

Use for cycling through tests common to certain types.
"""
const pair_types = (Vector2i, Vector2s, Vector2d, Point2i, Point2s, Point2d)
const pair_symbols = (:Vector2i, :Vector2s, :Vector2d, :Point2i, :Point2s, :Point2d)
pair_type_iterable = zip(pair_types, pair_symbols)

const pair_int_types = (pair_types[1], pair_types[3])
const pair_int_symbols = (pair_symbols[1], pair_symbols[3])
pair_int_iterable = zip(pair_int_types, pair_int_symbols)

const pair_float_types = (pair_types[2:3]..., pair_types[5:end]...)
const pair_float_symbols = (pair_symbols[2:3]..., pair_symbols[5:end]...)
pair_float_iterable = zip(pair_float_types, pair_float_symbols)

const pair_vector_types = pair_types[1:3]
const pair_vector_symbols = pair_symbols[1:3]
pair_vector_iterable = zip(pair_vector_types, pair_vector_symbols)

const pair_point_types = pair_types[4:end]
const pair_point_symbols = pair_symbols[4:end]
pair_point_iterable = zip(pair_point_types, pair_point_symbols)

IntegerPair = Union{Vector2i, Point2i}
FloatPair = Union{Vector2s, Point2s}
DoublePair = Union{Vector2d, Point2d}

VectorPair = Union{Vector2i, Vector2s, Vector2d}
PointPair = Union{Point2i, Point2s, Point2d}

function generate_set(type::Type)
    vcat([type(), type(rand(-5000:1:5000, 1)...)], [type(rand(-5000:1:5000, 2)...) for i in 1:2])
end