text:
using UInt12Arrays
using Test

@testset "UInt12Arrays.jl" begin
    data = UInt8[0x21, 0x43, 0x65, 0xba, 0xdc, 0xfe]
    A = UInt12Vector(data)
    T = eltype(A)
    @test size(A) == (4,)
    @test A[1] == 0x0321
    @test A[2] == 0x0654
    @test A[3] == 0x0cba
    @test A[4] == 0x0fed
    @test T == UInt12Arrays.default_eltype
    @test typeof(A[1]) == T
    resize!(A, 5)
    @test size(A) == (5,)
    A[5] = 0x0987
    @test A[5] == 0x0987
    @test size(data) == (8,)
    push!(A, 0x777)
    @test A[end] == 0x777
    @test_throws InexactError A[end] = 0x7777
    reference_array = [0x321, 0x654, 0xcba, 0xfed, 0x987, 0x777]
    @test UInt12Vector{UInt16}(data, 6) == reference_array
    @test UInt12Vector{UInt12}(data, 6) == reference_array
    @test UInt12Array(data) == reference_array
    @test UInt12Array{UInt16, Vector{UInt8}, 1}(data, (6,)) == reference_array
    B = UInt12Vector(undef, 5)
    @test size(B) == (5,)
    @test size(B.data) == (8,)
    C = UInt12Vector(undef, 6)
    @test size(C) == (6,)
    @test size(C.data) == (9,)
    @test convert(Array{UInt16}, A) == A
    @test convert(UInt12Array, @view(A[1:end])) == A
    @test convert(UInt12Array, @view(A[1:3])) == A[1:3]
    @test convert(UInt12Array, @view(A[1:2])) == A[1:2]
    @test convert(UInt12Array, @view(A[1:1])) == A[1:1]
    @test_broken convert(UInt12Array, @view(A[2:end])) == A[2:end]
    @test convert(UInt12Array, @view(A[3:end])) == A[3:end]
    @test_broken convert(UInt12Array, @view(A[4:end])) == A[4:end]
    @test convert(Array{UInt16}, @view(A[1:end])) == A
    @test convert(Array{UInt16}, @view(A[1:3])) == A[1:3]
    @test convert(Array{UInt16}, @view(A[1:2])) == A[1:2]
    @test convert(Array{UInt16}, @view(A[1:1])) == A[1:1]
    @test convert(Array{UInt16}, @view(A[2:end])) == A[2:end]
    @test convert(Array{UInt16}, @view(A[3:end])) == A[3:end]
    @test convert(Array{UInt16}, @view(A[4:end])) == A[4:end]
    @test convert(Array, @view(A[1:1])) == A[1:1]
    @test convert(Array, @view(A[2:end])) == A[2:end]
    @test convert(Array, @view(A[3:end])) == A[3:end]
    @test convert(Array, @view(A[4:end])) == A[4:end]
end


text:
""" Graphviz drawing of categories, functors, diagrams, and related structures.
"""
module GraphvizCategories
export to_graphviz

using ...CategoricalAlgebra, ...Graphs, ..GraphvizGraphs
import ..Graphviz
import ..GraphvizGraphs: to_graphviz, to_graphviz_property_graph

function to_graphviz(F::FinFunctor; kw...)
  to_graphviz(to_graphviz_property_graph(F; kw...))
end

function to_graphviz_property_graph(F::FinFunctor; kw...)
  g = graph(dom(F))
  pg = to_graphviz_property_graph(g; kw...)
  for v in vertices(g)
    lᵥ = g[v, :vname]
    tᵥ = ob_map(F, v)
    set_vprop!(pg, v, :label, string(lᵥ, ":", tᵥ))
  end
  for e in edges(g)
    tₑ = hom_map(F, e)
    set_eprop!(pg, e, :label, string(tₑ))
  end
  pg
end

end

ext: jl


text:
function SharedFaceVolume(S₂::AbstractArray{T, 2},
            αs_S₁vertices_intermsof_S₂vertices::AbstractArray{Float64, 2},
            ordered_vertex_indices_S1::AbstractVector{Int},
            ordered_vertex_indices_S2::AbstractVector{Int}) where {T}

    # Indices of the nonshared vertices in simplex 1 (label it u)
    # and in simplex 2 (label it v)
    n_vertices  = size(S₂, 2)
    idx_u = ordered_vertex_indices_S1[n_vertices]
    idx_v = ordered_vertex_indices_S2[n_vertices]

    αᵤ = αs_S₁vertices_intermsof_S₂vertices[:, idx_u]
    last_α = αᵤ[idx_v]
    shared_vertex_inds_in_S2 = ordered_vertex_indices_S2[1:n_vertices-1]

    # If last coefficient is zero or negative, then the simplices intersect
	# in the shared face. So we just have to consider
	if last_α <= 0.0
		return 0.0
    else
        inds_positive_coeffs = findall(αᵤ .> 0)
        positive_coeffs = αᵤ[inds_positive_coeffs]

        # One can prove that the intersection in the nontrivial case (last_α > 0)
        # is a new simplex. This means that there will be only one extra
        # intersecting point. It is the convex combination of the points
        # with strictly positive convex expansion coefficients with those
        # coefficients (hence, the normalisation to get convexity).
        intpt = (S₂[:, inds_positive_coeffs] * positive_coeffs) ./ sum(positive_coeffs)

        return abs(det([ones(1, n_vertices);
                        [S₂[:, shared_vertex_inds_in_S2] intpt]]))
    end
end

function SharedFaceVertices(S₂::AbstractArray{T, 2},
            αs_S₁vertices_intermsof_S₂vertices::AbstractArray{Float64, 2},
            ordered_vertex_indices_S1::AbstractVector{Int},
            ordered_vertex_indices_S2::AbstractVector{Int}) where {T}

    # Indices of the nonshared vertices in simplex 1 (label it u)
    # and in simplex 2 (label it v)
    n_vertices  = size(S₂, 2)
    idx_u = ordered_vertex_indices_S1[n_vertices]
    idx_v = ordered_vertex_indices_S2[n_vertices]

    αᵤ = αs_S₁vertices_intermsof_S₂vertices[:, idx_u]
    last_α = αᵤ[idx_v]
    shared_vertex_inds_in_S2 = ordered_vertex_indices_S2[1:n_vertices-1]

    # If last coefficient is zero or negative, then the simplices intersect
	# in the shared face. So we just have to consider
	if last_α <= 0.0
		return Float64[]
    else
        inds_positive_coeffs = findall(αᵤ .> 0)
        positive_coeffs = αᵤ[inds_positive_coeffs]

        # One can prove that the intersection in the nontrivial case (last_α > 0)
        # is a new simplex. This means that there will be only one extra
        # intersecting point. It is the convex combination of the points
        # with strictly positive convex expansion coefficients with those
        # coefficients (hence, the normalisation to get convexity).
        intpt = (S₂[:, inds_positive_coeffs] * positive_coeffs) ./ sum(positive_coeffs)

        return [S₂[:, shared_vertex_inds_in_S2] intpt]
    end
end


text:
module BioMart

    using HTTP, CSV, LightXML, DataFrames, GenomicFeatures, Memoization

    include("config.jl")
    include("Database.jl")
    include("Dataset.jl")
    include("Query.jl")
    include("Filter.jl")
    include("Attribute.jl")

end # module


text:
@testset "scattered field" begin
    k0 = 0.1
    kin = 3k0
    θ_i = rand()*2π

    N = 1000
    shapes = [squircle(1, N);
            rounded_star(0.8, 0.2, 3, N);
            rounded_star(0.75, 0.25, 5, N);
            ellipse(0.2, 0.9, N)]
    for s in shapes
        sigma_mu = get_potentialPW(k0, kin, s, θ_i)
        dx = 2*norm(s.ft[1,:] - s.ft[2,:])
        #normal vector
        nvec = [s.dft[:,2] -s.dft[:,1]]./hypot.(s.dft[:,1],s.dft[:,2])
        p_in = s.ft - nvec*dx
        p_out = s.ft + nvec*dx
        u_in = scatteredfield(sigma_mu, kin, s, p_in)
        u_out = scatteredfield(sigma_mu, k0, s, p_out)
        #some error is allowed as this quadrature is imprecise close to the boundary
        u_inc = uinc(k0, s.ft, PlaneWave(θ_i))
        @test norm(u_inc + u_out - u_in)/norm(u_inc) < 1e-3
    end
end


text:
@testset "FIFOBuffer" begin

    f = HTTP.FIFOBuffer(0)
    @test read(f, Tuple{UInt8,Bool}) == (0x00, false)
    @test isempty(readavailable(f))

    f = HTTP.FIFOBuffer(1)
    @test read(f, Tuple{UInt8,Bool}) == (0x00, false)
    @test isempty(readavailable(f))

    @test write(f, 0x01) == 1
    @test write(f, 0x02) == 0
    @test read(f, Tuple{UInt8,Bool}) == (0x01, true)
    @test read(f, Tuple{UInt8,Bool}) == (0x00, false)
    @test isempty(readavailable(f))

    @test write(f, UInt8[0x01, 0x02]) == 1
    @test all(readavailable(f) .== UInt8[0x01])

    f = HTTP.FIFOBuffer(5)
    @test read(f, Tuple{UInt8,Bool}) == (0x00, false)
    @test isempty(readavailable(f))

    write(f, 0x01)
    @test read(f, Tuple{UInt8,Bool}) == (0x01, true)
    @test read(f, Tuple{UInt8,Bool}) == (0x00, false)
    @test isempty(readavailable(f))

    write(f, 0x01)
    write(f, 0x02)
    @test read(f, Tuple{UInt8,Bool}) == (0x01, true)
    @test read(f, Tuple{UInt8,Bool}) == (0x02, true)
    @test read(f, Tuple{UInt8,Bool}) == (0x00, false)
    @test isempty(readavailable(f))

    write(f, 0x01)
    write(f, 0x02)
    @test all(readavailable(f) .== UInt8[0x01, 0x02])
    @test read(f, Tuple{UInt8,Bool}) == (0x00, false)
    @test isempty(readavailable(f))

    write(f, 0x01)
    write(f, 0x02)
    write(f, 0x03)
    write(f, 0x04)
    write(f, 0x05)
    @test read(f, Tuple{UInt8,Bool}) == (0x01, true)
    @test read(f, Tuple{UInt8,Bool}) == (0x02, true)
    @test read(f, Tuple{UInt8,Bool}) == (0x03, true)
    @test read(f, Tuple{UInt8,Bool}) == (0x04, true)
    @test read(f, Tuple{UInt8,Bool}) == (0x05, true)
    @test read(f, Tuple{UInt8,Bool}) == (0x00, false)
    @test isempty(readavailable(f))

    write(f, 0x01)
    write(f, 0x02)
    write(f, 0x03)
    write(f, 0x04)
    write(f, 0x05)
    write(f, 0x06) == 0
    @test all(readavailable(f) .== UInt8[0x01, 0x02, 0x03, 0x04, 0x05])
    @test read(f, Tuple{UInt8,Bool}) == (0x00, false)
    @test isempty(readavailable(f))

    write(f, UInt8[])
    @test read(f, Tuple{UInt8,Bool}) == (0x00, false)
    @test isempty(readavailable(f))

    write(f, UInt8[0x01])
    @test read(f, Tuple{UInt8,Bool}) == (0x01, true)
    @test read(f, Tuple{UInt8,Bool}) == (0x00, false)
    @test isempty(readavailable(f))

    write(f, UInt8[0x01, 0x02])
    @test read(f, Tuple{UInt8,Bool}) == (0x01, true)
    @test read(f, Tuple{UInt8,Bool}) == (0x02, true)
    @test read(f, Tuple{UInt8,Bool}) == (0x00, false)
    @test isempty(readavailable(f))

    write(f, UInt8[0x01, 0x02])
    @test all(readavailable(f) .== UInt8[0x01, 0x02])
    @test read(f, Tuple{UInt8,Bool}) == (0x00, false)
    @test isempty(readavailable(f))

    write(f, UInt8[0x01, 0x02, 0x03, 0x04, 0x05])
    @test read(f, Tuple{UInt8,Bool}) == (0x01, true)
    @test read(f, Tuple{UInt8,Bool}) == (0x02, true)
    @test read(f, Tuple{UInt8,Bool}) == (0x03, true)
    @test read(f, Tuple{UInt8,Bool}) == (0x04, true)
    @test read(f, Tuple{UInt8,Bool}) == (0x05, true)
    @test read(f, Tuple{UInt8,Bool}) == (0x00, false)
    @test isempty(readavailable(f))

    write(f, UInt8[0x01, 0x02, 0x03, 0x04, 0x05])
    @test all(readavailable(f) .== UInt8[0x01, 0x02, 0x03, 0x04, 0x05])
    @test read(f, Tuple{UInt8,Bool}) == (0x00, false)
    @test isempty(readavailable(f))

    # overflow
    write(f, UInt8[0x01, 0x02, 0x03, 0x04, 0x05, 0x06]) == 5
    @test all(readavailable(f) .== UInt8[0x01, 0x02, 0x03, 0x04, 0x05])
    @test read(f, Tuple{UInt8,Bool}) == (0x00, false)
    @test isempty(readavailable(f))

    # condition notification
    # fill the buffer up
    @test write(f, UInt8[0x01, 0x02, 0x03, 0x04, 0x05]) == 5
    # a write task is started asynchronously, which means it will wait for a
    # notify from a read that there's space to write again
    tsk = @async begin
        @test write(f, UInt8[0x01, 0x02, 0x03, 0x04, 0x05]) == 5
    end
    # when data is read in the next readavailble, a notify is sent to f.cond
    # which wakes up `tsk` to write it's data
    @test all(readavailable(f) .== UInt8[0x01, 0x02, 0x03, 0x04, 0x05])
    # meanwhile, we start the next `readavailable` asynchronously, which will block
    # until `tsk` is done writing
    @sync @async @test all(readavailable(f) .== UInt8[0x01, 0x02, 0x03, 0x04, 0x05])
    @sync begin
        N = 100
        tsk1 = @async begin
            for i = 1:N
                @test write(f, UInt8[0x4a, 0x61, 0x63, 0x6f, 0x62]) == 5
            end
        end
        tsk2 = @async begin
            for i = 1:N
                @test all(readavailable(f) .== UInt8[0x4a, 0x61, 0x63, 0x6f, 0x62])
            end
        end
    end

    # buffer growing
    f = HTTP.FIFOBuffer(10)
    @test write(f, UInt8[0x01, 0x02, 0x03, 0x04, 0x05]) == 5
    @test write(f, UInt8[0x06, 0x07, 0x08, 0x09, 0x0a]) == 5
    @test all(readavailable(f) .== 0x01:0x0a)

    # read
    f = HTTP.FIFOBuffer(5)
    @test write(f, UInt8[0x01, 0x02, 0x03, 0x04, 0x05]) == 5
    @test all(read(f, 5) .== 0x01:0x05)
    @test write(f, UInt8[0x01, 0x02, 0x03, 0x04, 0x05]) == 5
    @test all(read(f, 6) .== 0x01:0x05)
    @test write(f, UInt8[0x01, 0x02, 0x03, 0x04, 0x05]) == 5
    @test isempty(read(f, 0))
    @test all(read(f, 2) .== 0x01:0x02)
    @test write(f, 0x01) == 1
    @test all(read(f, 5) .== UInt8[0x03, 0x04, 0x05, 0x01])
    @test write(f, UInt8[0x01, 0x02, 0x03, 0x04, 0x05]) == 5
    @test all(read(f, 2) .== 0x01:0x02)
    r = read(f, 3)
    @test all(r .== 0x03:0x05)

    f2 = HTTP.FIFOBuffer(f)
    @test f == f2

    f = HTTP.FIFOBuffer(5)
    @test isempty(read(f, 1))
    t = @async read(f, 1)
    write(f, 0x01)
    @test wait(t) == [0x01]

    @test write(f, [0x01, 0x02, 0x03, 0x04, 0x05]) == 5
    @test write(f, [0x01, 0x02]) == 0

    @test readavailable(f) == [0x01, 0x02, 0x03, 0x04, 0x05]

    # ensure we're in a wrap-around state
    f = HTTP.FIFOBuffer(5)
    @test write(f, [0x01, 0x02, 0x03]) == 3
    @test readavailable(f) == [0x01, 0x02, 0x03]
    @test write(f, [0x01, 0x02, 0x03, 0x04]) == 4
    @test f.f > f.l

    @test write(f, [0x05]) == 1
    @test readavailable(f) == [0x01, 0x02, 0x03, 0x04, 0x05]

    @test write(f, [0x01, 0x02, 0x03, 0x04]) == 4
    @test write(f, [0x05, 0x06]) == 1
    @test readavailable(f) == [0x01, 0x02, 0x03, 0x04, 0x05]

    # ensure that `read(..., ::Type{UInt8})` returns a `UInt8`
    # https://github.com/JuliaWeb/HTTP.jl/issues/41
    f = HTTP.FIFOBuffer(5)
    b = Array{UInt8}(3)
    @test write(f, [0x01, 0x02, 0x03, 0x04]) == 4
    @test readbytes!(f, b) == 3
    @test b == [0x01, 0x02, 0x03]
    @test read(f, UInt8) == 0x04
    @test_throws EOFError read(f, UInt8)

    # ensure we return eof == false if there are still bytes to be read
    f = HTTP.FIFOBuffer(5)
    write(f, [0x01, 0x02, 0x03, 0x04])
    close(f)
    @async begin
        @test !eof(f)
    end

    # Issue #45
    # Ensure that we don't encounter an EOF when reading before data is written
    f = HTTP.FIFOBuffer(5)
    bytes = [0x01, 0x02, 0x03, 0x04]
    @test !eof(f)
    @sync begin
        @async begin
            bytes_read = UInt8[]
            while !eof(f)
                push!(bytes_read, read(f, UInt8))
            end
            @test bytes_read == bytes
        end
        yield()
        write(f, bytes)
        close(f)
    end
    @test eof(f)
end; # testset


text:
struct SimControls
    nu::Int64
    ctrl::MJMatrix
    normalized::Vector{mjtNum}
    ctrl_magnitude::MJMatrix

    function SimControls(data::MJData)
        mat = data.ctrl
        norm = zeros(data.nu)
        rng = data.model.actuator_ctrlrange

        new(data.nu, mat, norm, rng)
    end
end

function SetSimControls(ctrl::SimControls, norm)
    for i=1:ctrl.nu
        ctrl.normalized[i] = norm[i]
        mag = ctrl.ctrl_magnitude[i]
        adj = norm[i] * mag

        if (adj > mag)
            adj = mag
        elseif (adj < -mag)
            adj = -mag
        end

        ctrl.ctrl[i] = adj
    end

    return ctrl
end

function getitem(ctrl::SimControls, i)
    return ctrl.ctrl[i, 1]
end

function setitem!(ctrl::SimControls, val, i)
    ctrl.ctrl.mat[1, i] = val
end

mutable struct MJSimEnv{StateType, ControlType}
    model::MJModel
    data::MJData

    state::StateType
    controls::ControlType

    reward::Function
    testTerminal::Function

    t::Int
    terminal::Bool

    stateType
    controlType

    function MJSimEnv{StateType, ControlType}(m::MJModel, r::Function, tF::Function, s::StateType, c::ControlType) where {StateType, ControlType}
        d = MJData(m)
        mj_resetData(m, d)

        new(m, d, StateType(d), ControlType(d), r, tF, 0, false, StateType, ControlType)
    end
end

function GetState(env::MJSimEnv)
    env.state = env.stateType(env.data)
end

function GetReward(env::MJSimEnv, state=env.state)
    return env.reward(state)
end

function IsTerminalState(env::MJSimEnv, state=env.state)
    terminal = env.testTerminal(state)
    return terminal
end

#function SetControls(env::MJSimEnv, controls=env.controls)
#    for i=1:env.data.nu
#        env.data.ctrl[i] = controls.ctrl[i]
#    end
#end

#function Act(env::MJSimEnv, num_steps::Int = 3, controls=env.controls)
#    for i=1:num_steps
#        #had to do weird set controls each frame in v1.31, TODO?
#        mj_step1(env.model, env.data)
#        for i=1:env.data.nu
#            env.data.ctrl.mat[1, i] = controls.ctrl[i]
#        end
#        mj_step2(env.model, env.data)
#    end
#end

function ResetEnv(env::MJSimEnv)
    env.terminal = false
    env.t = 0
    mj_resetData(env.model, env.data)
    env.state = env.stateType(env.data)
    env.controls = env.controlType(env.data)
end


text:
@testset "LP Solver" begin
    @testset "Issue 83 TimeLimit" begin
        glpk_optimizer =
            optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => GLPK.GLP_MSG_OFF)
        model = Model(optimizer_with_attributes(
            CS.Optimizer,
            "lp_optimizer" => glpk_optimizer,
            "logging" => [],
            "time_limit" => 0.01,
        ))

        # Variables
        @variable(model, inclusion[h = 1:3], Bin)
        @variable(model, 0 <= allocations[h = 1:3, a = 1:3] <= 1, Int)
        @variable(model, 0 <= days[h = 1:3, a = 1:3] <= 5, Int)

        # Constraints
        @constraint(
            model,
            must_include[h = 1:3],
            sum(allocations[h, a] for a in 1:3) <= inclusion[h]
        )
        # at least n
        @constraint(model, min_hospitals, sum(inclusion[h] for h in 1:3) >= 3)
        # every h must be allocated at most one a
        @constraint(model, must_visit[h = 1:3], sum(allocations[h, a] for a in 1:3) <= 1)
        # every allocated h must have fewer than 5 days of visits per week
        @constraint(
            model,
            max_visits[h = 1:3],
            sum(days[h, a] for a in 1:3) <= 5 * inclusion[h]
        )

        @objective(model, Max, sum(days[h, a] * 5 for h in 1:3, a in 1:3))
        optimize!(model)
        @test JuMP.termination_status(model) == MOI.TIME_LIMIT ||
              JuMP.termination_status(model) == MOI.OPTIMAL
        if JuMP.termination_status(model) == MOI.OPTIMAL
            @test JuMP.solve_time(model) <= 0.1
        end
    end

    @testset "Issue 83 max_bt_steps" begin
        glpk_optimizer =
            optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => GLPK.GLP_MSG_OFF)
        model = Model(optimizer_with_attributes(
            CS.Optimizer,
            "lp_optimizer" => glpk_optimizer,
            "logging" => [],
            "max_bt_steps" => 2,
        ))

        # Variables
        @variable(model, inclusion[h = 1:3], Bin)
        @variable(model, 0 <= allocations[h = 1:3, a = 1:3] <= 1, Int)
        @variable(model, 0 <= days[h = 1:3, a = 1:3] <= 5, Int)

        # Constraints
        @constraint(
            model,
            must_include[h = 1:3],
            sum(allocations[h, a] for a in 1:3) <= inclusion[h]
        )
        # at least n
        @constraint(model, min_hospitals, sum(inclusion[h] for h in 1:3) >= 3)
        # every h must be allocated at most one a
        @constraint(model, must_visit[h = 1:3], sum(allocations[h, a] for a in 1:3) <= 1)
        # every allocated h must have fewer than 5 days of visits per week
        @constraint(
            model,
            max_visits[h = 1:3],
            sum(days[h, a] for a in 1:3) <= 5 * inclusion[h]
        )

        @objective(model, Max, sum(days[h, a] * 5 for h in 1:3, a in 1:3))
        optimize!(model)
        if JuMP.termination_status(model) != MOI.OTHER_LIMIT
            @test JuMP.objective_value(model) ≈ 75
        end
    end

    @testset "Issue 83" begin
        glpk_optimizer =
            optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => GLPK.GLP_MSG_OFF)
        model = Model(optimizer_with_attributes(
            CS.Optimizer,
            "lp_optimizer" => glpk_optimizer,
            "logging" => [],
        ))

        # Variables
        @variable(model, inclusion[h = 1:3], Bin)
        @variable(model, 0 <= allocations[h = 1:3, a = 1:3] <= 1, Int)
        @variable(model, 0 <= days[h = 1:3, a = 1:3] <= 5, Int)

        # Constraints
        @constraint(
            model,
            must_include[h = 1:3],
            sum(allocations[h, a] for a in 1:3) <= inclusion[h]
        )
        # at least n
        @constraint(model, min_hospitals, sum(inclusion[h] for h in 1:3) >= 3)
        # every h must be allocated at most one a
        @constraint(model, must_visit[h = 1:3], sum(allocations[h, a] for a in 1:3) <= 1)
        # every allocated h must have fewer than 5 days of visits per week
        @constraint(
            model,
            max_visits[h = 1:3],
            sum(days[h, a] for a in 1:3) <= 5 * inclusion[h]
        )

        @objective(model, Max, sum(days[h, a] * 5 for h in 1:3, a in 1:3))
        optimize!(model)
        @test JuMP.objective_value(model) ≈ 75

        com = CS.get_inner_model(model)
    end

    @testset "Combine lp with table constraint" begin
        cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
        model = Model(optimizer_with_attributes(
            CS.Optimizer,
            "lp_optimizer" => cbc_optimizer,
            "logging" => [],
        ))

        # Variables
        @variable(model, 1 <= x[1:10] <= 15, Int)

        # Constraints
        @constraint(model, sum(x[1:5]) >= 10)
        @constraint(model, sum(x[6:10]) < 16)

        table_left = [
            1 2 3 4 5 # sum 16
            1 2 3 4 6 #     17
            4 5 5 3 4 #     21
        ]

        table_right = [
            1 2 3 1 1 # sum 8
            1 3 2 2 1 #     9
            1 1 1 1 4 #     8
        ]

        @constraint(model, x[1:5] in CS.TableSet(table_left))
        @constraint(model, x[6:10] in CS.TableSet(table_right))

        @objective(model, Max, sum(x))
        optimize!(model)

        @test JuMP.objective_value(model) ≈ 30
        @test sum(JuMP.value.(x[6:10])) ≈ 9
        @test sum(JuMP.value.(x[1:5])) ≈ 21
    end

    @testset "LP with bigger table constraint" begin
        cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
        model = Model(optimizer_with_attributes(
            CS.Optimizer,
            "lp_optimizer" => cbc_optimizer,
            "logging" => [],
            "keep_logs" => true,
        ))

        # Variables
        @variable(model, 1 <= x[1:10] <= 15, Int)

        # Constraints
        @constraint(model, sum(x[1:4]) > 9)
        @constraint(model, sum(x[5:10]) <= 15)

        table_left = Array{Int,2}(undef, (0, 5))
        for r in permutations(1:10, 5)
            s = sum(r)
            if 10 < s < 20 || 28 < s < 47
                table_left = vcat(table_left, r')
            end
        end

        table_right = [
            1 2 3 1 1 # sum 8
            1 3 3 2 1 #    10
            1 1 3 2 1 #     8
            1 1 1 2 4 #     9
            4 5 5 3 4 #     21
        ]

        @constraint(model, x[1:5] in CS.TableSet(table_left))
        @constraint(model, x[6:10] in CS.TableSet(table_right))

        @objective(model, Max, sum(x))
        optimize!(model)

        @test JuMP.termination_status(model) == MOI.OPTIMAL
        @test JuMP.objective_value(model) ≈ 49
        @test sum(JuMP.value.(x[1:4])) >= 10
        @test sum(JuMP.value.(x[5:10])) <= 15
        @test allunique(JuMP.value.(x[1:5]))
    end

    @testset "Combine lp with all different + not equal" begin
        cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
        model = Model(optimizer_with_attributes(
            CS.Optimizer,
            "lp_optimizer" => cbc_optimizer,
            "logging" => [],
            "traverse_strategy" => :DBFS,
        ))

        # Variables
        @variable(model, 1 <= x[1:10] <= 15, Int)

        # Constraints
        @constraint(model, sum(x[1:5]) >= 10)
        @constraint(model, sum(x[6:10]) <= 15)
        # disallow x[1] = 11 and x[2] == 12
        @constraint(model, 1.5 * x[1] + 2 * x[2] != 40.5)
        @constraint(model, x in CS.AllDifferentSet())

        @objective(model, Max, sum(x))
        optimize!(model)
        # possible solution 11+12+13+14+15  + 1+2+3+4+5
        # only works fast if the all different bound works
        @test JuMP.objective_value(model) ≈ 80
        @test sum(JuMP.value.(x[6:10])) ≈ 15
        @test !(1.5 * JuMP.value(x[1]) + 2 * JuMP.value(x[2]) ≈ 40.5)
    end

    @testset "Combine lp with all different + not equal + DFS" begin
        cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
        model = Model(optimizer_with_attributes(
            CS.Optimizer,
            "lp_optimizer" => cbc_optimizer,
            "logging" => [],
            "traverse_strategy" => :DFS,
            "branch_split" => :Biggest,
            "keep_logs" => true,
        ))

        # Variables
        @variable(model, 1 <= x[1:10] <= 15, Int)

        # Constraints
        @constraint(model, sum(x[1:5]) >= 10)
        @constraint(model, sum(x[6:10]) <= 15)
        # disallow x[1] = 11 and x[2] == 12
        @constraint(model, 1.5 * x[1] + 2 * x[2] != 40.5)
        @constraint(model, x in CS.AllDifferentSet())

        @objective(model, Max, sum(x))
        optimize!(model)
        # possible solution 11+12+13+14+15  + 1+2+3+4+5
        # only works fast if the all different bound works
        @test JuMP.objective_value(model) ≈ 80
        @test sum(JuMP.value.(x[6:10])) ≈ 15
        @test !(1.5 * JuMP.value(x[1]) + 2 * JuMP.value(x[2]) ≈ 40.5)

        com = CS.get_inner_model(model)
        @test is_solved(com)
        @test general_tree_test(com)
    end

    @testset "Combine lp with all different + not equal + DFS + Split Half" begin
        cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
        model = Model(optimizer_with_attributes(
            CS.Optimizer,
            "lp_optimizer" => cbc_optimizer,
            # "logging" => [],
            "traverse_strategy" => :DFS,
            "branch_split" => :InHalf,
        ))

        # Variables
        @variable(model, 1 <= x[1:10] <= 15, Int)

        # Constraints
        @constraint(model, sum(x[1:5]) >= 10)
        @constraint(model, sum(x[6:10]) <= 15)
        # disallow x[1] = 11 and x[2] == 12
        @constraint(model, 1.5 * x[1] + 2 * x[2] != 40.5)
        @constraint(model, x in CS.AllDifferentSet())

        @objective(model, Max, sum(x))
        optimize!(model)
        # possible solution 11+12+13+14+15  + 1+2+3+4+5
        # only works fast if the all different bound works
        @test JuMP.objective_value(model) ≈ 80
        @test sum(JuMP.value.(x[6:10])) ≈ 15
        @test !(1.5 * JuMP.value(x[1]) + 2 * JuMP.value(x[2]) ≈ 40.5)
    end
end


text:
module linalg_functions
    include("gram_schmidt.jl")
end


text:
#vulnerability curve
Pfun(ψ::T,ψ₅₀::T,b::T) where {T<:Float64} = (1/2)^(ψ/ψ₅₀)^b

#Invers vulnerability curve
Pfunᵢₙᵥ(Pval::T,ψ₅₀::T,b::T) where {T<:Float64} = ψ₅₀*(log(Pval)/log(0.5))^(1/b)

#Effective saturation
CalcSₑ(θₛ::T;θₛₐₜ::T=0.41,θᵣ::T=0.006) where {T<:Float64} = (θₛ-θᵣ)/(θₛₐₜ-θᵣ)

#Soil water content (volumetric) to soil water potential (MPa) (Water retention curve)
function θₛ2ψₛ(θₛ::T;θₛₐₜ::T=0.41,θᵣ::T=0.006,λ::T=1.0,ψₐ::T=-0.097706) where {T<:Float64}
    Sₑ = CalcSₑ(θₛ;θₛₐₜ=θₛₐₜ,θᵣ=θᵣ)

    return ψₐ*Sₑ^(-1/λ)
end

#Relative Soil conductance (ratio betwen actual and maximal)
function Re_Kₛᵣfun(θₛ::T;θₛₐₜ::T=0.41,θᵣ::T=0.006,p::T=4.66) where {T<:Float64}
    Sₑ = CalcSₑ(θₛ;θₛₐₜ=θₛₐₜ,θᵣ=θᵣ)
    return Sₑ^p
end

#Calculate canopy conductance
function Calc_K_cost(gₛ::T,model::CCPHStruct;limit_up::T=1.0,limit_lo::T=0.12) where {T<:Float64}
    ψ₅₀,b,Kₓₗ₀,g,ρ_H2O,θₛ = model.hydPar.ψ₅₀,model.hydPar.b,model.hydPar.Kₓₗ₀,model.cons.g,model.cons.ρ_H2O,model.env.θₛ

    #Calculate tree height
    H = model.treesize.H
    #Calculate transpiration
    E = 1.6*gₛ*model.env.VPD/model.env.P
    #Caluclate soil water potential
    ψₛ = θₛ2ψₛ(θₛ)
    #Calcualate soil conductance
    re_Kₛᵣ = Re_Kₛᵣfun(θₛ)
    Kₛᵣ =   model.hydPar.Kₛᵣ₀*re_Kₛᵣ

    #Soil water potential adjusted for gravitational pressure (MPa)
    ψₛ_g = ψₛ-H*ρ_H2O*g*10^-6

    #Calculate canopy conductance
    try
        K_cost = bisection(x->x-Pfun(ψₛ_g-E*(2*Kₓₗ₀*x+Kₛᵣ)/(2*Kₓₗ₀*x*Kₛᵣ),ψ₅₀,b),limit_lo, limit_up) #lower than 0.12 resultts in hydraulic failure
        Kₓₗ = Kₓₗ₀*K_cost
        ψ_c = ψₛ_g-E*(Kₓₗ+Kₛᵣ)/(Kₓₗ*Kₛᵣ)
        return K_cost, Kₓₗ, ψ_c
    catch
        error("Could not find a feasable canopy conductance")
    end
end

#Calculate stomatal (gₛ) for a given K_cost (Pval)
function Calc_K_costᵢₙᵥ(Pval::Float64,model::CCPHStruct)
    ψ₅₀,b,Kₓₗ₀,g,ρ_H2O,θₛ = model.hydPar.ψ₅₀,model.hydPar.b,model.hydPar.Kₓₗ₀,model.cons.g,model.cons.ρ_H2O,model.env.θₛ

    #Calc target ψ
    ψ_target = Pfunᵢₙᵥ(Pval,ψ₅₀,b)
    #Calc taget Kₓₗ
    Kₓₗ = Kₓₗ₀*Pval
    #Calculate tree height
    H = model.treesize.H
    #Caluclate soil water potential
    ψₛ = θₛ2ψₛ(θₛ)
    #Calcualate soil conductance
    re_Kₛᵣ = Re_Kₛᵣfun(θₛ)
    Kₛᵣ =   model.hydPar.Kₛᵣ₀*re_Kₛᵣ
    #Soil water potential adjusted for gravitational pressure (MPa)
    ψₛ_g = ψₛ-H*ρ_H2O*g*10^-6
    #Calculate target E
    E_target = -(ψ_target-ψₛ_g)*(2*Kₓₗ*Kₛᵣ)/(2*Kₓₗ+Kₛᵣ)
    #Calcualte target gₛ
    gₛ_target = E_target*model.env.P/(1.6*model.env.VPD)

    return gₛ_target
end


text:
using Statistics, StatsBase
using PyCall
xr = pyimport("xarray")
using Blosc
using Elefridge
using JLD
using ZfpCompression
using TranscodingStreams,CodecZstd

ZstdCompressorL3 = ZstdCompressor(level=3)
TranscodingStreams.initialize(ZstdCompressorL3)

ZstdCompressorL10 = ZstdCompressor(level=10)
TranscodingStreams.initialize(ZstdCompressorL10)

ZstdCompressorL22 = ZstdCompressor(level=22)
TranscodingStreams.initialize(ZstdCompressorL22)

path = "/Users/milan/cams/unstructured/"
filelist = filter(x->endswith(x,".grib"),readdir(path))

n = length(filelist)
varnames = fill("",n)
methods = ["ZstdL3","ZstdL10","ZstdL22"]
stats = ["cf","meanerror","absmedian","abs90%","absmax","decmedian","dec90%","decmax"]

E = load("/Users/milan/cams/error/round_zstd_all_unstr.jld")
E = E["E"]

D = load("/Users/milan/cams/entropy/keepbits_98.jld")

## RUN
function rs_zfp(r::Int)
    if r > 6
        rd = 6
    elseif r > 2
        rd = 5
    else
        rd = 4
    end
    return r+rd
end

function statscalc(X0::Array{T,N},Xa::Array{T,N}) where {T,N}
    Xamean = mean(Xa)
    X0mean = mean(X0)
    meanerror = abs(X0mean-Xamean)/X0mean       # normalised

    X0amean = mean(abs.(X0))

    abserr = abs.(X0-Xa)/X0amean                # normalised
    abserr = sort(vec(abserr))
    abs50 = T(quantile(abserr,.5,sorted=true))
    abs90 = T(quantile(abserr,.9,sorted=true))
    absmax = abserr[end]

    # exception cases for decimal error
    Xa0 = Xa .== zero(T)
    X00 = X0 .== zero(T)

    Xas = signbit.(Xa)
    X0s = signbit.(X0)

    decerr = abs.(log10.(abs.(X0./Xa)))

    decerr[Xas .⊻ X0s] .= Inf        # signs don't match
    decerr[Xa0 .⊻ X00] .= Inf        # one is zero the other one not
    decerr[Xa0 .& X00] .= zero(T)    # both zero, set decimal error to zero

    decerr = sort(vec(decerr))
    dec50 = T(quantile(decerr,.5,sorted=true))
    dec90 = T(quantile(decerr,.9,sorted=true))
    decmax = decerr[end]

    return [meanerror,abs50,abs90,absmax,dec50,dec90,decmax]
end

for (i,file) in enumerate(filelist)
    varname = split(split(file,"0_")[end],".")[1]
    varnames[i] = varname

    println("---")
    println("Reading $varname")
    X = xr.open_dataarray(joinpath(path,file),engine="cfgrib").data

    # transpose array
    # use only half the globe
    X = permutedims(X[:,1:(348528÷2)],[2,1])  # unravels as longitude, latitude, level

    # number of bits of information
    r = D["infbits"][D["varnames"] .== varname][1]-9
    r = max(r,2)    # don't use less than two

    println("$r real information bits.")

    Xr = round(X,r)

    ## Round+Zstd
    E[i,1,2:end] = statscalc(X,Xr)
    E[i,2,2:end] = E[i,1,2:end]         # copy across for Zstd
    E[i,3,2:end] = E[i,1,2:end]

    Xr8 = unsafe_wrap(Array, Ptr{UInt8}(pointer(Xr)), sizeof(Xr))

    E[i,1,1] = sizeof(X)/sizeof(transcode(ZstdCompressorL3,Xr8))
    println("ZstdL3 finished.")
    E[i,2,1] = sizeof(X)/sizeof(transcode(ZstdCompressorL10,Xr8))
    println("ZstdL10 finished.")
    E[i,3,1] = sizeof(X)/sizeof(transcode(ZstdCompressorL22,Xr8))
    println("ZstdL22 finished.")

    println(E[i,:,1])
    @save "/Users/milan/cams/error/round_zstd_all_unstr.jld" varnames methods stats E
end


text:
# TODO: sub_pixel_smoothing in 2-dim case, take care of smoothing along edges
# TODO: fix Bravais in 1-dimension

################################################################################
###############   SUB PIXEL SMOOTHING    #######################################
################################################################################
function construct_εFr(bnd::Boundary, dis::Discretization, sys::System)
    x = dis.X
    y = dis.Y
    domains = Array{Int}(undef,size(x)...)
    which_domains!(domains, x, y, bnd, sys)
    r = domains

    xb, yb = Array{Float64}(undef,size(x)...), Array{Float64}(undef,size(y)...)
    bravais_coordinates_unit_cell!(xb,yb,x,y,domains,sys)

    ε = Array{ComplexF64}(undef,dis.N[1],dis.N[2])
    F = Array{Float64}(undef,dis.N[1],dis.N[2])
    for i ∈ eachindex(r)
        j = r[i]
        ε[i] = sys.ε_by_region[j](x[i],y[i],sys.params_by_region[j])
        F[i] = sys.F_by_region[j](x[i],y[i],sys.params_by_region[j])
    end
    return ε, F, r
end

"""
    ε, regions = sub_pixel_smoothing(bnd, dis, sys; display=false)
"""
function sub_pixel_smoothing!(bnd::Boundary, dis::Discretization{Cartesian}, sys::System, ε, F, r; display::Bool=false)
    x = dis.X
    y = dis.Y
    sub_pixel_num = dis.sub_pixel_num
    if dis.N[2] == 1

    elseif dis.N[1] == 1

    else
        sub_x = Array{Float64}(undef, sub_pixel_num, sub_pixel_num)
        sub_y = Array{Float64}(undef, sub_pixel_num, sub_pixel_num)
        xb = Array{Float64}(undef, sub_pixel_num, sub_pixel_num)
        yb = Array{Float64}(undef, sub_pixel_num, sub_pixel_num)
        domains = Array{Int}(undef, sub_pixel_num, sub_pixel_num)
        ε_temp = Array{ComplexF64}(undef,sub_pixel_num, sub_pixel_num)
        F_temp = Array{Float64}(undef,sub_pixel_num, sub_pixel_num)
        if display
            pg = Progress((size(x,1)-2)*(size(y,2)-2), PROGRESS_UPDATE_TIME::Float64, "sub-pixel smoothing ")
        end
        for h ∈ CartesianIndices(x)
            i,j = h[1],h[2]
            if all((1,1).<(i,j).<size(x))
                nearestNeighborFlag = r[i,j]!==r[i,j+1] || r[i,j]!==r[i,j-1] || r[i,j]!==r[i+1,j] || r[i,j]!==r[i-1,j]
                nextNearestNeighborFlag = r[i,j]!==r[i+1,j+1] || r[i,j]!==r[i-1,j-1] || r[i,j]!==r[i+1,j-1] || r[i,j]!==r[i-1,j+1]
                if nearestNeighborFlag || nextNearestNeighborFlag
                    x_min = (x[i,j]+x[i-1,j])/2
                    x_max = (x[i,j]+x[i+1,j])/2
                    y_min = (y[i,j]+y[i,j-1])/2
                    y_max = (y[i,j]+y[i,j+1])/2
                    for k ∈ CartesianIndices(sub_x)
                        sub_x[k] = x_min + (k[1]-1)*(x_max-x_min)/(sub_pixel_num-1)
                        sub_y[k] = y_min + (k[2]-1)*(y_max-y_min)/(sub_pixel_num-1)
                    end
                    which_domains!(domains, sub_x, sub_y, bnd, sys)
                    bravais_coordinates_unit_cell!(xb, yb, sub_x, sub_y, domains, sys)
                    for l ∈ eachindex(xb)
                        d = domains[l]
                        ε_temp[l] = sys.ε_by_region[d](sub_x[d],sub_y[d],sys.params_by_region[d])
                        F_temp[l] = sys.F_by_region[d](sub_x[d],sub_y[d],sys.params_by_region[d])
                    end
                    ε[i,j] = mean(ε_temp)
                    F[i,j] = mean(F_temp)
                end
                if display
                    next!(pg)
                end
            end
        end
    end
    return nothing
end
function sub_pixel_smoothing!(bnd::Boundary, dis::Discretization{Polar}, sys::System, ε, F, r; display::Bool=false)
    R = broadcast((a,b)->a,dis.x...)
    Θ = broadcast((a,b)->b,dis.x...)
    sub_pixel_num = dis.sub_pixel_num
    if dis.N[2] == 1

    else
        sub_x = Array{Float64}(undef, sub_pixel_num, sub_pixel_num)
        sub_y = Array{Float64}(undef, sub_pixel_num, sub_pixel_num)
        domains = Array{Int}(undef, sub_pixel_num, sub_pixel_num)
        ε_temp = Array{ComplexF64}(undef,sub_pixel_num, sub_pixel_num)
        F_temp = Array{Float64}(undef,sub_pixel_num, sub_pixel_num)
        if display
            pg = Progress((size(x,1)-2)*(size(y,2)-2), PROGRESS_UPDATE_TIME::Float64, "sub-pixel smoothing ")
        end
        for h ∈ CartesianIndices(r)
            i,j = h[1],h[2]
            if all((1,1).<(i,j).<size(r))
                nearestNeighborFlag = r[i,j]!==r[i,j+1] || r[i,j]!==r[i,j-1] || r[i,j]!==r[i+1,j] || r[i,j]!==r[i-1,j]
                nextNearestNeighborFlag = r[i,j]!==r[i+1,j+1] || r[i,j]!==r[i-1,j-1] || r[i,j]!==r[i+1,j-1] || r[i,j]!==r[i-1,j+1]
                if nearestNeighborFlag || nextNearestNeighborFlag
                    r_min = (R[i,j]+R[i-1,j])/2
                    r_max = (R[i,j]+R[i+1,j])/2
                    θ_min = (Θ[i,j]+Θ[i,j-1])/2
                    θ_max = (Θ[i,j]+Θ[i,j+1])/2
                    for kj ∈ 1:sub_pixel_num
                        sub_θ = θ_min + (kj-1)*(θ_max-θ_min)/(sub_pixel_num-1)
                        sc = sincos(sub_θ)
                        for ki ∈ 1:sub_pixel_num
                            sub_r = r_min + (ki-1)*(r_max-r_min)/(sub_pixel_num-1)
                            sub_x[ki,kj] = sub_r*sc[2]
                            sub_y[ki,kj] = sub_r*sc[1]
                        end
                    end
                    which_domains!(domains, sub_x, sub_y, bnd, sys)
                    for l ∈ eachindex(domains)
                        d = domains[l]
                        ε_temp[l] = sys.ε_by_region[d](sub_x[d],sub_y[d],sys.params_by_region[d])
                        F_temp[l] = sys.F_by_region[d](sub_x[d],sub_y[d],sys.params_by_region[d])
                    end
                    ε[i,j] = mean(ε_temp)
                    F[i,j] = mean(F_temp)
                end
                if display
                    next!(pg)
                end
            end
        end
    end
    return nothing
end

function which_domains(x,y,bnd,sys)
    domains = zeros(Int,size(x)...)
    which_domains!(domains,x,y,bnd,sys)
    return domains
end
function which_domains!(domains,x,y,bnd,sys) # need to add boundary
    for i ∈ eachindex(x)
        domains[i] = which_domain(x[i],y[i],bnd,sys)
    end
    return nothing
end

"""
    which_domain(x, y, bnd, sys)
"""
function which_domain(x, y, bnd::Boundary, sys::System)
    domain = findfirst(map(z->z.is_in_domain(x,y),sys.domains))
    @assert !isnothing(domain) "domain not found for (x,y)=($x,$y). has a background domain been defined?"
    return domain
end

################################################################################
###############  STANDARD DOMAINS  #######################################
################################################################################
"""
    left_domain(x, y, domain_index, bnd, sys)
"""
left_domain(x,y,idx::Int,bnd::Boundary,sys::System) =
    sys.domains[idx].which_asymptote == :left && x<bnd.∂Ω_tr[1][1]

"""
    right_domain(x, y, domain_index, bnd, sys)
"""
function right_domain(x,y,idx::Int,bnd::Boundary,sys::System)
    return sys.domains[idx].which_asymptote == :right && x>bnd.∂Ω_tr[1][2]
end

"""
    bottom_domain(x, y, domain_index, bnd, sys)
"""
function bottom_domain(x,y,idx::Int,bnd::Boundary,sys::System)
    return sys.domains[idx].which_asymptote == :bottom && y<bnd.∂Ω_tr[2][1]
end

"""
    top_domain(x, y, domain_index, bnd, sys)
"""
function top_domain(x,y,idx::Int,bnd::Boundary,sys::System)
    return sys.domains[idx].which_asymptote == :top && y>bnd.∂Ω_tr[2][2]
end


text:
# Titulos constituintes dos indices S&P500, DJIA, NIKKEY225, FTSE100, DAX30, CAC40, PSI20 e EURO STOXX 50, em 31/03/2021

# The S&P 500 stock market index, maintained by S&P Dow Jones Indices, comprises 505 common stocks issued by 500 large-cap companies and traded on American stock exchanges (including the 30 companies that compose the Dow Jones Industrial Average), and covers about 80 percent of the American equity market by capitalization. Although called the S&P 500, the index contains 505 stocks because it includes two share classes of stock from 5 of its component companies.

const SP500n = ["A", "AAL", "AAP", "AAPL", "ABBV", "ABC", "ABMD", "ABT", "ACN", "ADBE", "ADI", "ADM", "ADP", "ADSK", "AEE", "AEP", "AES", "AFL", "AIG", "AIZ", "AJG", "AKAM", "ALB", "ALGN", "ALK", "ALL", "ALLE", "ALXN", "AMAT", "AMCR", "AMD", "AME", "AMGN", "AMP", "AMT", "AMZN", "ANET", "ANSS", "ANTM", "AON", "AOS", "APA", "APD", "APH", "APTV", "ARE", "ATO", "ATVI", "AVB", "AVGO", "AVY", "AWK", "AXP", "AZO", "BA", "BAC", "BAX", "BBY", "BDX", "BEN", "BF-B", "BIIB", "BIO", "BK", "BKNG", "BKR", "BLK", "BLL", "BMY", "BR", "BRK-B", "BSX", "BWA", "BXP", "C", "CAG", "CAH", "CARR", "CAT", "CB", "CBOE", "CBRE", "CCI", "CCL", "CDNS", "CDW", "CE", "CERN", "CF", "CFG", "CHD", "CHRW", "CHTR", "CI", "CINF", "CL", "CLX", "CMA", "CMCSA", "CME", "CMG", "CMI", "CMS", "CNC", "CNP", "COF", "COG", "COO", "COP", "COST", "CPB", "CPRT", "CRM", "CSCO", "CSX", "CTAS", "CTLT", "CTSH", "CTVA", "CTXS", "CVS", "CVX", "CZR", "D", "DAL", "DD", "DE", "DFS", "DG", "DGX", "DHI", "DHR", "DIS", "DISCA", "DISCK", "DISH", "DLR", "DLTR", "DOV", "DOW", "DPZ", "DRE", "DRI", "DTE", "DUK", "DVA", "DVN", "DXC", "DXCM", "EA", "EBAY", "ECL", "ED", "EFX", "EIX", "EL", "EMN", "EMR", "ENPH", "EOG", "EQIX", "EQR", "ES", "ESS", "ETN", "ETR", "ETSY", "EVRG", "EW", "EXC", "EXPD", "EXPE", "EXR", "F", "FANG", "FAST", "FB", "FBHS", "FCX", "FDX", "FE", "FFIV", "FIS", "FISV", "FITB", "FLIR", "FLT", "FMC", "FOX", "FOXA", "FRC", "FRT", "FTNT", "FTV", "GD", "GE", "GILD", "GIS", "GL", "GLW", "GM", "GNRC", "GOOG", "GOOGL", "GPC", "GPN", "GPS", "GRMN", "GS", "GWW", "HAL", "HAS", "HBAN", "HBI", "HCA", "HD", "HES", "HFC", "HIG", "HII", "HLT", "HOLX", "HON", "HPE", "HPQ", "HRL", "HSIC", "HST", "HSY", "HUM", "HWM", "IBM", "ICE", "IDXX", "IEX", "IFF", "ILMN", "INCY", "INFO", "INTC", "INTU", "IP", "IPG", "IPGP", "IQV", "IR", "IRM", "ISRG", "IT", "ITW", "IVZ", "J", "JBHT", "JCI", "JKHY", "JNJ", "JNPR", "JPM", "K", "KEY", "KEYS", "KHC", "KIM", "KLAC", "KMB", "KMI", "KMX", "KO", "KR", "KSU", "L", "LB", "LDOS", "LEG", "LEN", "LH", "LHX", "LIN", "LKQ", "LLY", "LMT", "LNC", "LNT", "LOW", "LRCX", "LUMN", "LUV", "LVS", "LW", "LYB", "LYV", "MA", "MAA", "MAR", "MAS", "MCD", "MCHP", "MCK", "MCO", "MDLZ", "MDT", "MET", "MGM", "MHK", "MKC", "MKTX", "MLM", "MMC", "MMM", "MNST", "MO", "MOS", "MPC", "MPWR", "MRK", "MRO", "MS", "MSCI", "MSFT", "MSI", "MTB", "MTD", "MU", "MXIM", "NCLH", "NDAQ", "NEE", "NEM", "NFLX", "NI", "NKE", "NLOK", "NLSN", "NOC", "NOV", "NOW", "NRG", "NSC", "NTAP", "NTRS", "NUE", "NVDA", "NVR", "NWL", "NWS", "NWSA", "NXPI", "O", "ODFL", "OKE", "OMC", "ORCL", "ORLY", "OTIS", "OXY", "PAYC", "PAYX", "PBCT", "PCAR", "PEAK", "PEG", "PENN", "PEP", "PFE", "PFG", "PG", "PGR", "PH", "PHM", "PKG", "PKI", "PLD", "PM", "PNC", "PNR", "PNW", "POOL", "PPG", "PPL", "PRGO", "PRU", "PSA", "PSX", "PVH", "PWR", "PXD", "PYPL", "QCOM", "QRVO", "RCL", "RE", "REG", "REGN", "RF", "RHI", "RJF", "RL", "RMD", "ROK", "ROL", "ROP", "ROST", "RSG", "RTX", "SBAC", "SBUX", "SCHW", "SEE", "SHW", "SIVB", "SJM", "SLB", "SNA", "SNPS", "SO", "SPG", "SPGI", "SRE", "STE", "STT", "STX", "STZ", "SWK", "SWKS", "SYF", "SYK", "SYY", "T", "TAP", "TDG", "TDY", "TEL", "TER", "TFC", "TFX", "TGT", "TJX", "TMO", "TMUS", "TPR", "TRMB", "TROW", "TRV", "TSCO", "TSLA", "TSN", "TT", "TTWO", "TWTR", "TXN", "TXT", "TYL", "UA", "UAA", "UAL", "UDR", "UHS", "ULTA", "UNH", "UNM", "UNP", "UPS", "URI", "USB", "V", "VAR", "VFC", "VIAC", "VLO", "VMC", "VNO", "VRSK", "VRSN", "VRTX", "VTR", "VTRS", "VZ", "WAB", "WAT", "WBA", "WDC", "WEC", "WELL", "WFC", "WHR", "WLTW", "WM", "WMB", "WMT", "WRB", "WRK", "WST", "WU", "WY", "WYNN", "XEL", "XLNX", "XOM", "XRAY", "XYL", "YUM", "ZBH", "ZBRA", "ZION", "ZTS", "^GSPC"]

const DJ30n = [:AAPL, :AMGN, :AXP, :BA, :CAT, :CRM, :CSCO, :CVX, :DD, :DIS, :GS, :HD, :HON, :IBM, :INTC, :JNJ, :JPM, :KO, :MCD, :MMM, :MRK,:MSFT, :NKE, :PG, :TRV, :UNH, :V,  :VZ, :WBA, :WMT, "^DJI", "^GSPC"]

const NIKKEI225n = ["6857.T", "8267.T", "5201.T", "2802.T", "6770.T", "6113.T", "9202.T", "8304.T", "2502.T", "3407.T", "4503.T", "7832.T", "5108.T", "7751.T", "6952.T", "9022.T", "9502.T", "4519.T", "7762.T", "1721.T", "7186.T", "8253.T", "4751.T", "7912.T", "8750.T", "4568.T", "6367.T", "1925.T", "8601.T", "2432.T", "4061.T", "6902.T", "4324.T", "4631.T", "5714.T", "9020.T", "6361.T", "4523.T", "5020.T", "6954.T", "9983.T", "6504.T", "4901.T", "5803.T", "6702.T", "8354.T", "5801.T", "6674.T", "1808.T", "7205.T", "6305.T", "7004.T", "6501.T", "7267.T", "5019.T", "7013.T", "1605.T", "3099.T", "7202.T", "8001.T", "3086.T", "8697.T", "6178.T", "2914.T", "5411.T", "1963.T", "6473.T", "1812.T", "4452.T", "7012.T", "9107.T", "9433.T", "9008.T", "9009.T", "2801.T", "2503.T", "5406.T", "6301.T", "9766.T", "4902.T", "6326.T", "3405.T", "6971.T", "4151.T", "2413.T", "8002.T", "1333.T", "8252.T", "8628.T", "7261.T", "2269.T", "6479.T", "4188.T", "8058.T", "6503.T", "8802.T", "7011.T", "9301.T", "5711.T", "7211.T", "8306.T", "8031.T", "4183.T", "7003.T", "8801.T", "5706.T", "9104.T", "8411.T", "8725.T", "6701.T", "3659.T", "5333.T", "2282.T", "2871.T", "7731.T", "5214.T", "9062.T", "5703.T", "3863.T", "5202.T", "5401.T", "1332.T", "9432.T", "9101.T", "4021.T", "7201.T", "2002.T", "3105.T", "6988.T", "8604.T", "6471.T", "6472.T", "9613.T", "1802.T", "9007.T", "3861.T", "6703.T", "6103.T", "7733.T", "6645.T", "9532.T", "4578.T", "5541.T", "6752.T", "4755.T", "6098.T", "8308.T", "7752.T", "2501.T", "7735.T", "9735.T", "6724.T", "1928.T", "3382.T", "6753.T", "1803.T", "4063.T", "8303.T", "4507.T", "4911.T", "4004.T", "9412.T", "9434.T", "9984.T", "2768.T", "8630.T", "6758.T", "7270.T", "3436.T", "4005.T", "8053.T", "4506.T", "5802.T", "6302.T", "5713.T", "8316.T", "8309.T", "5232.T", "8830.T", "7269.T", "8795.T", "5233.T", "1801.T", "6976.T", "2531.T", "8233.T", "4502.T", "6762.T", "3401.T", "4543.T", "8331.T", "5631.T", "9503.T", "8355.T", "5101.T", "9001.T", "9602.T", "5707.T", "5301.T", "8766.T", "4043.T", "9501.T", "8035.T", "9531.T", "8804.T", "9005.T", "3289.T", "7911.T", "3402.T", "4042.T", "5332.T", "5901.T", "3101.T", "7203.T", "8015.T", "4704.T", "4208.T", "3103.T", "9021.T", "7951.T", "7272.T", "9064.T", "6506.T", "6841.T", "4689.T", "^N225"]

const FTSE100n = ["III.L", "ADM.L", "AAL.L", "ANTO.L", "AHT.L", "ABF.L", "AZN.L", "AUTO.L", "AVST.L", "AVV.L", "AV.L", "BME.L", "BA.L", "BARC.L", "BDEV.L", "BKG.L", "BHP.L", "BP.L", "BATS.L", "BLND.L", "BT-A.L", "BNZL.L", "BRBY.L", "CCH.L", "CPG.L", "CRH.L", "CRDA.L", "DCC.L", "DGE.L", "ENT.L", "EVR.L", "EXPN.L", "FERG.L", "FLTR.L", "FRES.L", "GSK.L", "GLEN.L", "HLMA.L", "HL.L", "HIK.L", "HSBA.L", "IHG.L", "IMB.L", "INF.L", "ICP.L", "IAG.L", "ITRK.L", "JD.L", "JMAT.L", "JET.L", "KGF.L", "LAND.L", "LGEN.L", "LLOY.L", "LSEG.L", "MNG.L", "MRO.L", "MNDI.L", "NG.L", "NWG.L", "NXT.L", "OCDO.L", "PSON.L", "PSH.L", "PSN.L", "PHNX.L", "POLY.L", "PRU.L", "RB.L", "REL.L", "RSW.L", "RTO.L", "RMV.L", "RIO.L", "RR.L", "RDSA.L", "RSA.L", "SGE.L", "SBRY.L", "SDR.L", "SMT.L", "SGRO.L", "SVT.L", "SMDS.L", "SMIN.L", "SN.L", "SKG.L", "SPX.L", "SSE.L", "STAN.L", "SLA.L", "STJ.L", "TW.L", "TSCO.L", "ULVR.L", "UU.L", "VOD.L", "WEIR.L", "WTB.L", "WPP.L", "^FTSE"]

const DAX30n = ["ADS.DE", "ALV.DE", "BAS.DE", "BAYN.DE", "BEI.DE", "BMW.DE", "CON.DE", "1COV.DE", "DAI.DE", "DHER.DE", "DBK.DE", "DB1.DE", "DPW.DE", "DTE.DE", "DWNI.DE", "EOAN.DE", "FRE.DE", "FME.DE", "HEI.DE", "HEN3.DE", "IFX.DE", "LIN.DE", "MRK.DE", "MTX.DE", "MUV2.DE", "RWE.DE", "SAP.DE", "SIE.DE", "VOW3.DE", "VNA.DE", "^GDAXI"]

const CAC40n = ["AI.PA", "AIR.PA", "ALO.PA", "MT.AS", "ATO.PA", "CS.PA", "BNP.PA", "EN.PA", "CAP.PA", "CA.PA", "ACA.PA", "BN.PA", "DSY.PA", "ENGI.PA", "EL.PA", "RMS.PA", "KER.PA", "OR.PA", "LR.PA", "MC.PA", "ML.PA", "ORA.PA", "RI.PA", "PUB.PA", "RNO.PA", "SAF.PA", "SGO.PA", "SAN.PA", "SU.PA", "GLE.PA", "STLA.PA", "STM.PA", "TEP.PA", "HO.PA", "FP.PA", "URW.AS", "VIE.PA", "DG.PA", "VIV.PA", "WLN.PA", "^FCHI"]

const FTSEMIBn = ["A2A.MI", "AMP.MI", "ATL.MI", "AZM.MI", "BGN.MI", "BAMI.MI", "BPE.MI", "BRE.MI", "BZU.MI", "CPR.MI", "CNHI.MI", "DIA.MI", "ENEL.MI", "ENI.MI", "EXO.MI", "RACE.MI", "FBK.MI", "G.MI", "ISP.MI", "IG.MI", "JUVE.MI", "LDO.MI", "MB.MI", "MONC.MI", "PIRC.MI", "PST.MI", "PRY.MI", "REC.MI", "SPM.MI", "SFER.MI", "SRG.MI", "STLA.MI", "STM.MI", "TIT.MI", "TEN.MI", "TRN.MI", "UBI.MI", "UCG.MI", "UNI.MI", "US.MI", "FTSEMIB.MI"]

const PSI20n = ["ALTR.LS", "BCP.LS", "COR.LS", "CTT.LS", "EDP.LS", "EDPR.LS", "GALP.LS", "IBS.LS", "JMT.LS", "EGL.LS", "NBA.LS", "NOS.LS", "NVG.LS", "PHR.LS", "RAM.LS", "RENE.LS", "SEM.LS", "SON.LS", "PSI20.LS"]

const EUROSTOXX50n = ["ADS.DE", "AD.AS", "AI.PA", "AIR.PA", "ALV.DE", "ABI.BR", "ASML.AS", "AMS.MC", "CS.PA", "BBVA.MC", "SAN.MC", "BAS.DE", "BAYN.DE", "BMW.DE", "BNP.PA", "CRG.IR", "SGO.PA", "DAI.DE", "DPW.DE", "DTE.DE", "ENEL.MI", "ENGI.PA", "ENI.MI", "EOAN.DE", "EL.PA", "FRE.DE", "BN.PA", "IBE.MC", "ITX.MC", "INGA.AS", "ISP.MI", "LIN.DE", "OR.PA", "MC.PA", "MUV2.DE", "NOKIA.HE", "ORA.PA", "PHIA.AS", "SAF.PA", "SAN.PA", "SAP.DE", "SU.PA", "SIE.DE", "GLE.PA", "TEF.MC", "FP.PA", "URW.AS", "UNA.AS", "DG.PA", "VIV.PA", "VOW.DE"]

const ASX50n = ["A2M.AX", "AGL.AX", "AMC.AX", "ANZ.AX", "APA.AX", "ALD.AX", "ALL.AX", "ASX.AX", "AZJ.AX", "BHP.AX", "BXB.AX", "COH.AX", "COL.AX", "CBA.AX", "CPU.AX", "CSL.AX", "DXS.AX", "FMG.AX", "GMG.AX", "GPT.AX", "IAG.AX", "JHX.AX", "LLC.AX", "MQG.AX", "MPL.AX", "MGR.AX", "NAB.AX", "NCM.AX", "OSH.AX", "ORI.AX", "ORG.AX", "QAN.AX", "QBE.AX", "RHC.AX", "RIO.AX", "STO.AX", "SCG.AX", "SHL.AX", "S32.AX", "SGP.AX", "SUN.AX", "SYD.AX", "TLS.AX", "TCL.AX", "TWE.AX", "VCX.AX", "WES.AX", "WBC.AX", "WPL.AX", "WOW.AX", "^AFLI"]

const HKSEn = ["0005.HK", "0011.HK", "0388.HK", "0939.HK", "1299.HK", "1398.HK", "2318.HK", "2388.HK", "2628.HK", "3328.HK", "3988.HK", "0002.HK", "0003.HK", "0006.HK", "1038.HK", "0012.HK", "0016.HK", "0017.HK", "0101.HK", "0688.HK", "0823.HK", "0960.HK", "1109.HK", "1113.HK", "1997.HK", "2007.HK", "0001.HK", "0027.HK", "0066.HK", "0175.HK", "0241.HK", "0267.HK", "0288.HK", "0386.HK", "0669.HK", "0700.HK", "0762.HK", "0857.HK", "0883.HK", "0941.HK", "1044.HK", "1093.HK", "1177.HK", "1810.HK", "1876.HK", "1928.HK", "2018.HK", "2020.HK", "2269.HK", "2313.HK", "2319.HK", "2382.HK", "3690.HK", "6862.HK", "9988.HK", "^HSI"]

const SMIn = ["NESN.SW", "NOVN.SW", "ROG.SW", "ZURN.SW", "UBSG.SW", "ABBN.SW", "CFR.SW", "LONN.SW", "SIKA.SW", "ALC.SW", "GIVN.SW", "LHN.SW", "SCMN.SW", "CSGN.SW", "PGHN.SW", "SREN.SW", "GEBN.SW", "SGSN.SW", "SLHN.SW", "UHR.SW"]

const FXn = ["EURUSD=X", "EURGBP=X", "EURJPY=X", "EURCHF=X", "EURSEK=X", "EURDKK=X", "EURCAD=X", "EURAUD=X", "EURKRW=X", "HKD=X"]

"""
    yahoo(simbolos, data_inicial, data_final, intervalo = "1d")

Obter dados de Yahoo.Finance

## Argumentos
* `simbolos`: simbolos de mercado, e.g "AAPL", "^GSPC", ["IBM", "FB"]
* `data_inicial`: Data inicial do periodo, e.g. Date(2019,12,20) ou DateTime(2019,12,20,8,30,0)
* `data_final`: Data final do periodo, e.g. Date(2020,12,20) ou DateTime(2020,12,20,8,30,0)
* `intervalo`: Intervalo de amostragem, e.g. "1d" (default), "1wk", "1mo", "3mo"
"""
function yahoo(symbol, date1 = Date(1900,1,1), date2 = Date(Dates.now()), interval::String = "1d")
    if isa(symbol, String) == true
        nsymb = 1
    elseif isa(symbol, Symbol) == true
        nsymb = 1
    else
        nsymb = length(symbol)
    end

    date1 = DateTime(date1)
    date2 = DateTime(date2)
    if date1 > date2
        date1, date2 = date2, date1
    end
    date2 = date2 +  Day(1)
    from = string(round(Int64, datetime2unix(date1)))
    to = string(round(Int64, datetime2unix(date2)))
    host = rand(["query1", "query2"])
    for i in 1:nsymb
        if nsymb == 1 && (isa(symbol, String) == true || isa(symbol, Symbol) == true)
            symb = String(symbol)
        else
            symb = String(symbol[i])
        end
        println(symb)
        url = "https://$host.finance.yahoo.com/v7/finance/chart/$symb?&interval=$interval&period1=$from&period2=$to"
        response = HTTP.get(url, cookies = true)
        body = JSON.parse(String(response.body))["chart"]["result"][1]
        # values = body["indicators"]["quote"][1]
        a = replace!(body["indicators"]["adjclose"][1]["adjclose"], nothing => missing)
        new = DataFrame(
            Time = Dates.Date.(unix2datetime.(body["timestamp"])),
            # Open = convert(Vector{Float64},values["open"]),
            # High =  convert(Vector{Float64},values["high"]),
            # Low = convert(Vector{Float64},values["low"]),
            # Close = convert(Vector{Float64},values["close"]),
            #AdjClose = convert(Vector{Float64},body["indicators"]["adjclose"][1]["adjclose"]),
            AdjClose = convert(Vector{Union{Missing, Float64}},a),
            # Volume = convert(Vector{Int64},values["volume"])
        )

        DataFrames.rename!(new,[:Data, Symbol(symb)])
        #DataFrames.deleterows!(new, isnothing.(new).Close)
        if i == 1
            global old = new
        else
            old = outerjoin(old, new, on = :Data)
        end
    end
    return old
end
yahoo(s::Vector{Symbol}, date1, date2, interval) = yahoo(String.(s), date1, date2, interval)
yahoo(s::Vector{Symbol}, date1, date2) = yahoo(String.(s), date1, date2, "1d")
yahoo(s::Vector{Symbol}, date1) = yahoo(String.(s), date1, Date(Dates.now()), "1d")
yahoo(s::Symbol, date1, date2, interval) = yahoo(String(s), date1, date2, interval)
yahoo(s::Symbol, date1, date2) = yahoo(String(s), date1, date2, "1d")
yahoo(s::Symbol, date1) = yahoo(String(s), date1, Date(Dates.now()), "1d")

"""
    retornos(precos::DataFrame, tipo_retorno)

Calcular retornos (simples, ou logaritmicos) a partir de tabela de precos, com datas na primeira coluna

## Argumentos
* `precos`: tabela de precos, com datas na primeira coluna
* `tipo_retorno`: :simples, :log
"""
function retornos(df::DataFrame, tipo_retorno::Symbol)
    c = Matrix(df)

    cot = convert(Matrix{Union{Missing, Float64}},c[:,2:end])
    dt  = convert(Vector{Date}, c[:,1])
    nret = size(c,1) - 1
    ncol = size(c,2) - 1

    ret   = Array{Union{Missing, Float64}}(undef, nret, ncol)
    datas = Array{Date}(undef,nret)
    if tipo_retorno == :log
        for i = 1:nret
            ret[i,:] = log.(cot[i+1,:]./cot[i,:])
            datas[i] = dt[i+1]
        end
    elseif tipo_retorno == :simples
        for i = 1:nret
            ret[i,:] = cot[i+1,:]./cot[i,:] .- 1.0
            datas[i] = dt[i+1]
        end
    end
    d = DataFrame(Data = datas)
    nomes = names(df)
    popfirst!(nomes)
    r = DataFrame(ret, Symbol.(nomes))

    return  hcat(d, r)
end

function ETL(retornos, alfa, V)
    nobs = length(retornos)
    corte = quantile(retornos, alfa)
    excedimentos = 0
    n_excedimentos = 0
    for i in 1:nobs
        if retornos[i] <= corte
            excedimentos += retornos[i]
            n_excedimentos += 1
        end
    end
    return - (excedimentos/n_excedimentos - mean(retornos)) * V
end

function EWMA(retornos, lambda)
    m = length(retornos)
    media = mean(retornos)
    sigma2 = 0
    for k in 1:m
        sigma2 += (1-lambda)*lambda^(k-1) * (retornos[end-k+1] - media)^2
    end
    return sqrt(sigma2)
end

function semaforo(retornos, alfa)
    nobs = lenght(retornos)
    excep = zeros(nobs)
    sigma = std(retornos)
    RaR = - sigma * quantile(Normal(0,1), alfa)
    for i in 1:nobs
        if retornos[i] <= -RaR
            excep[i] = 1
        end
    end
    return sum(excep)
end

function cvm(Sigma)
    n = size(Sigma,1)
    modelo = Model(Ipopt.Optimizer)
    @variable(modelo, w[1:n])
    @objective(modelo, Min, w' * Sigma * w)
    @constraint(modelo, sum(w) == 1)
    @constraint(modelo, w .>= 0)
    optimize!(modelo)
    return value.(w)
end

function fe(mu, Sigma, mu_k)
    n = size(Sigma,1)
    modelo = Model(Ipopt.Optimizer)
    @variable(modelo, w[1:n])
    @objective(modelo, Min, w' * Sigma * w)
    @constraint(modelo, sum(w) == 1)
    @constraint(modelo, w .>= 0)
    @constraint(modelo, mu' * w >= mu_k)
    optimize!(modelo)
    return value.(w)
  end

function gfe(mu, Sigma, z_alfa, T, V_0, ncarteiras = 10)
    n = length(mu)
    modelo = Model(Ipopt.Optimizer)
    set_silent(modelo)
    @variable(modelo, w[1:n])
    @objective(modelo, Min, w' * Sigma * w)
    @constraint(modelo, sum(w) == 1)
    @constraint(modelo, w .>= 0)
    optimize!(modelo)
    w_CVM = value.(w)
    mu_CVM = mu' * w_CVM
    VaR_k = zeros(ncarteiras)
    mu_k = zeros(ncarteiras)
    for i in 0:ncarteiras-1
        modelo = Model(Ipopt.Optimizer)
        set_silent(modelo)
        @variable(modelo, w[1:n])
        @objective(modelo, Min, w' * Sigma * w)
        @constraint(modelo, sum(w) == 1)
        @constraint(modelo, w .>= 0)
        mu_max = maximum(mu)
        incremento = (mu_max-mu_CVM)/(ncarteiras-1)
        mu_k[i+1] = mu_CVM + incremento * i
        @constraint(modelo, mu' * w >= mu_k[i+1])
        optimize!(modelo)
        w_k = value.(w)
        VaR_k[i+1] = -z_alfa * sqrt(T) * sqrt(w_k' * Sigma * w_k) * V_0
    end
    aVaR = -z_alfa * sqrt(T) * sqrt.(diag(Sigma)) * V_0
    fig = plot(VaR_k,mu_k, xlabel = "VaR da Carteira (alfa = $(alfa))", ylabel = "Valor Esperado do Retorno da Carteira", label = "Fronteira Eficiente", xlim = (0, maximum(aVaR) * 1.1), ylim = (0, maximum(mu)*1.1), legend = :bottomright)
    fig = scatter!(aVaR, mu, label = "Ativos")
    VaR_CVM = -z_alfa * sqrt(T) * sqrt(w_CVM' * Sigma * w_CVM) * V_0
    fig = scatter!([VaR_CVM], [mu_CVM], label = "Carteira VaR Minimo")
    return fig
end

function alocar(mu, Sigma, lista, ncarteiras = 10)
    modelo = Model(Ipopt.Optimizer)
    set_silent(modelo)
    n = length(mu)
    @variable(modelo, w[1:n])
    @objective(modelo, Min, w' * Sigma * w)
    @constraint(modelo, sum(w) == 1)
    @constraint(modelo, w .>= 0)
    optimize!(modelo)
    w_CVM = value.(w)
    mu_k = zeros(ncarteiras)
    w_k = zeros(n,ncarteiras)
    mu_CVM = mu' * w_CVM
    for i in 0:ncarteiras-1
        modelo = Model(Ipopt.Optimizer)
        set_silent(modelo)
        @variable(modelo, w[1:n])
        @objective(modelo, Min, w' * Sigma * w)
        @constraint(modelo, sum(w) == 1)
        @constraint(modelo, w .>= 0)
        mu_max = maximum(mu)
        incremento = (mu_max-mu_CVM)/(ncarteiras-1)
        mu_k[i+1] = mu_CVM + incremento * i
        @constraint(modelo, mu' * w >= mu_k[i+1])
        optimize!(modelo)
        w_k[:,i+1] = value.(w)
    end
    println(w_k)
    ticklabel = "P" .* string.(collect(1:ncarteiras))
    ticklabel[1] = "CVM"
    ticklabel[ncarteiras] = "CRM"
    fig = groupedbar(w_k', bar_position = :stack, bar_width=1.0, xlabel = "Carteiras", xticks=(1:ncarteiras, ticklabel), ylabel = "Pesos", label = permutedims(String.(lista)), legend = :bottomright)
    return fig
end

#=
function installedx()
    @warn "Pkg.installed() is deprecated"
    deps = dependencies()
    installs = Dict{String, VersionNumber}()
    for (uuid, dep) in deps
        dep.is_direct_dep || continue
        dep.version === nothing && continue
        installs[dep.name] = dep.version
    end
    return installs
end
=#


text:
############################################################
# joSA-distribute/gather constructors            ###########
############################################################

############################################################
# joSAdistribute/gather constructors #######################

export joSAdistribute
"""
    julia> joSAdistribute(m [,parts]; [DT])
    julia> joSAdistribute(wpool, m [,parts]; [DT])
    julia> joSAdistribute(m, nvc [,parts]; [DT])
    julia> joSAdistribute(wpool, m, nvc [,parts]; [DT])
    julia> joSAdistribute(dst)
    julia> joSAdistribute(SA)

defines operator to distribute serial vector into DistributedArrays' vector

# Signature

    joSAdistribute(m::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(),m);
        kwargs...) where INT<:Integer
    joSAdistribute(wpool::WorkerPool,m::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(wpool),m);
        DT::DataType=joFloat,gclean::Bool=false) where INT<:Integer
    joSAdistribute(m::Integer,nvc::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(),nvc);
        kwargs...) where INT<:Integer
    joSAdistribute(wpool::WorkerPool,m::Integer,nvc::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(wpool),nvc);
        DT::DataType=joFloat,gclean::Bool=false) where INT<:Integer
    joSAdistribute(dst::joPAsetup;gclean::Bool=false)
    joSAdistribute(A::joSAdistributedLinearOperator;gclean::Bool=false)

# Arguments

- `m`: length of the vector
- `nvc`: # of columns in multi-vector - if given then multi-vector is distributed over 2nd dimension
- `parts`: custom partitioning of distributed dimension
- `wpool`: custom WorkerPool
- `DT`: DataType for joPAsetup
- `dst`: joPAsetup
- `SA`: joSAdistributedLinearOperator
- `glcean`: clean SArray after gathering

# Notes

- no type conversions are attempted at the moment (i.e. DT is used as for any other JOLI operator)
- adjoint/transpose of the joSAdistribute will gather distributed vector into serial vector

# Examples

- `joSAdistribute(5)`: distribute vector of lenght 5 into default WorkerPool
- `joSAdistribute(5,2)`: distribute multi-vector of lenght 5 with 2 columns into default WorkerPool
- `joSAdistribute(5)'`: gather vector of lenght 5
- `joSAdistribute(5,2)'`: gather multi-vector of lenght 5 with 2 columns

"""
function joSAdistribute(wpool::WorkerPool,m::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(wpool),m);
        DT::DataType=joFloat,gclean::Bool=false) where INT<:Integer

    length(parts)==nworkers(wpool) || throw(joPAsetupException("joSAdistribute: lenght(parts) does not much nworkers()"))
    sum(parts)==m || throw(joPAsetupException("joSAdistribute: sum(parts) does not much m"))
    dst=joPAsetup(wpool,(m,),1,parts=parts)
    return joSAdistribute{DT,DT,1}("joSAdistributeV",m,m,1,
        v1->scatter(v1,dst),
        v2->sdata(v2),
        v3->sdata(v3),
        v4->scatter(v4,dst),
        @joNF, @joNF, @joNF, @joNF,
        dst,gclean)
end
joSAdistribute(m::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(),m);
        kwargs...) where INT<:Integer = joSAdistribute(WorkerPool(workers()),m,parts;kwargs...)

function joSAdistribute(wpool::WorkerPool,m::Integer,nvc::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(wpool),nvc);
        DT::DataType=joFloat,gclean::Bool=false) where INT<:Integer

    length(parts)==nworkers(wpool) || throw(joPAsetupException("joSAdistribute: lenght(parts) does not much nworkers()"))
    sum(parts)==nvc || throw(joPAsetupException("joSAdistribute: sum(parts) does not much nvc"))
    dst=joPAsetup(wpool,(m,nvc),2,parts=parts)
    return joSAdistribute{DT,DT,2}("joSAdistributeMV:$nvc",m,m,nvc,
        v1->scatter(v1,dst),
        v2->sdata(v2),
        v3->sdata(v3),
        v4->scatter(v4,dst),
        @joNF, @joNF, @joNF, @joNF,
        dst,gclean)
end
joSAdistribute(m::Integer,nvc::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(),nvc);
        kwargs...) where INT<:Integer = joSAdistribute(WorkerPool(workers()),m,nvc,parts;kwargs...)
function joSAdistribute(dst::joPAsetup;gclean::Bool=false)
    m=dst.dims[1]
    nvc=dst.dims[2]
    DT=dst.DT
    return joSAdistribute{DT,DT,2}("joSAdistributeMV:$nvc",m,m,nvc,
        v1->scatter(v1,dst),
        v2->sdata(v2),
        v3->sdata(v3),
        v4->scatter(v4,dst),
        @joNF, @joNF, @joNF, @joNF,
        dst,gclean)
end
joSAdistribute(A::joSAdistributedLinearOperator;gclean::Bool=false) = joSAdistribute(A.PAs_in,gclean=gclean)

export joSAgather
"""
    julia> joSAgather(m [,parts]; [DT])
    julia> joSAgather(wpool, m [,parts]; [DT])
    julia> joSAgather(m, nvc [,parts]; [DT])
    julia> joSAgather(wpool, m, nvc [,parts]; [DT])
    julia> joSAgather(dst)
    julia> joSAgather(SA)

defines operator to gather DistributedArrays' vector into serial vector

# Signature

    joSAgather(m::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(),m);
        kwargs...) where INT<:Integer
    joSAgather(wpool::WorkerPool,m::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(wpool),m);
        DT::DataType=joFloat,gclean::Bool=false) where INT<:Integer
    joSAgather(m::Integer,nvc::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(),nvc);
        kwargs...) where INT<:Integer
    joSAgather(wpool::WorkerPool,m::Integer,nvc::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(wpool),nvc);
        DT::DataType=joFloat,gclean::Bool=false) where INT<:Integer
    joSAgather(dst::joPAsetup;gclean::Bool=false)
    joSAgather(A::joSAdistributedLinearOperator;gclean::Bool=false)

# Arguments

- `m`: length of the vector
- `nvc`: # of columns in multi-vector - if given then multi-vector is distributed over 2nd dimension
- `parts`: custom partitioning of distributed dimension
- `wpool`: custom WorkerPool
- `DT`: DataType for joPAsetup
- `glcean`: clean SArray after gathering
- `dst`: joPAsetup
- `SA`: joSAdistributedLinearOperator

# Notes

- no type conversions are attempted at the moment (i.e. DT is used as for any other JOLI operator)
- adjoint/transpose of the joSAgather will distribute serial vector into DistributedArrays' vector

# Examples

- `joSAgather(5)`: gather vector of lenght 5
- `joSAgather(5,2)`: gather multi-vector of lenght 5 with 2 columns
- `joSAgather(5)'`: distribute vector of lenght 5 into default WorkerPool
- `joSAgather(5,2)'`: distribute multi-vector of lenght 5 with 2 columns into default WorkerPool

"""
function joSAgather(wpool::WorkerPool,m::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(wpool),m);
        DT::DataType=joFloat,gclean::Bool=false) where INT<:Integer

    length(parts)==nworkers(wpool) || throw(joPAsetupException("joSAgather: lenght(parts) does not much nworkers()"))
    sum(parts)==m || throw(joPAsetupException("joSAgather: sum(parts) does not much m"))
    dst=joPAsetup(wpool,(m,),1,parts=parts)
    return joSAgather{DT,DT,1}("joSAgatherV",m,m,1,
        v1->sdata(v1),
        v2->scatter(v2,dst),
        v3->scatter(v3,dst),
        v4->sdata(v4),
        @joNF, @joNF, @joNF, @joNF,
        dst,gclean)
end
joSAgather(m::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(),m);
        kwargs...) where INT<:Integer = joSAgather(WorkerPool(workers()),m,parts;kwargs...)

function joSAgather(wpool::WorkerPool,m::Integer,nvc::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(wpool),nvc);
        DT::DataType=joFloat,gclean::Bool=false) where INT<:Integer

    length(parts)==nworkers(wpool) || throw(joPAsetupException("joSAgather: lenght(parts) does not much nworkers()"))
    sum(parts)==nvc || throw(joPAsetupException("joSAgather: sum(parts) does not much nvc"))
    dst=joPAsetup(wpool,(m,nvc),2,parts=parts)
    return joSAgather{DT,DT,2}("joSAgatherMV:$nvc",m,m,nvc,
        v1->sdata(v1),
        v2->scatter(v2,dst),
        v3->scatter(v3,dst),
        v4->sdata(v4),
        @joNF, @joNF, @joNF, @joNF,
        dst,gclean)
end
joSAgather(m::Integer,nvc::Integer,
        parts::Vector{INT}=joPAsetup_etc.balanced_partition(nworkers(),nvc);
        kwargs...) where INT<:Integer = joSAgather(WorkerPool(workers()),m,nvc,parts;kwargs...)
function joSAgather(dst::joPAsetup;gclean::Bool=false)
    m=dst.dims[1]
    nvc=dst.dims[2]
    DT=dst.DT
    return joSAgather{DT,DT,2}("joSAgatherMV:$nvc",m,m,nvc,
        v1->sdata(v1),
        v2->scatter(v2,dst),
        v3->scatter(v3,dst),
        v4->sdata(v4),
        @joNF, @joNF, @joNF, @joNF,
        dst,gclean)
end
joSAgather(A::joSAdistributedLinearOperator;gclean::Bool=false) = joSAgather(A.PAs_out,gclean=gclean)


text:
# extensions of merge methods

# Home of all extended merge! methods
merge(S::EventTraceData; v::Int64=KW.v) = (U = deepcopy(S); merge!(U, v=v); return U)
merge!(S::EventTraceData, U::EventTraceData; v::Int64=KW.v) = ([append!(getfield(S, f), getfield(U, f)) for f in tracefields]; S.n += U.n; merge!(S; v=v))
merge!(S::EventTraceData, C::EventChannel; v::Int64=KW.v) = merge!(S, EventTraceData(C), v=v)

function merge(A::Array{EventTraceData,1}; v::Int64=KW.v)
  L::Int64 = length(A)
  n = sum([A[i].n for i = 1:L])
  T = EventTraceData(n)
  [setfield!(T, f, vcat([getfield(A[i],f) for i = 1:L]...)) for f in SeisIO.datafields]
  merge!(T, v=v)
  return T
end
merge(S::EventTraceData, U::EventTraceData; v::Int64=KW.v) = merge(Array{EventTraceData,1}([S,U]), v=v)
merge(S::EventTraceData, C::EventChannel; v::Int64=KW.v) = merge(S, EventTraceData(C), v=v)
merge(C::EventChannel, S::EventTraceData; v::Int64=KW.v) = merge(EventTraceData(C), S, v=v)
merge(C::EventChannel, D::EventChannel; v::Int64=KW.v) = (S = EventTraceData(C,D); merge!(S, v=v); return S)

+(S::EventTraceData, C::EventChannel) = +(S, EventTraceData(C))
+(C::EventChannel, S::EventTraceData) = +(S, EventTraceData(C))
+(C::EventChannel, D::EventChannel) = +(EventTraceData(C), EventTraceData(D))

# Multiplication
# distributivity: (S1+S2)*S3) == (S1*S3 + S2*S3)
*(S::EventTraceData, U::EventTraceData) = merge(Array{EventTraceData,1}([S,U]))
*(S::EventTraceData, C::EventChannel) = merge(S, EventTraceData(C))
function *(C::EventChannel, D::EventChannel)
  s1 = deepcopy(C)
  s2 = deepcopy(D)
  S = merge(EventTraceData(s1),EventTraceData(s2))
  return S
end


text:
using Documenter, PolygonMinkowskiSum

makedocs(
    modules = [PolygonMinkowskiSum],
    format = Documenter.HTML(),
    checkdocs = :exports,
    sitename = "PolygonMinkowskiSum.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/goretkin/PolygonMinkowskiSum.jl.git",
)


text:
include("../src/UGM.jl")

include("testSmallDemo.jl")


text:
@testset "Polygons" begin

@testset "construction" begin
    points = Point2f0[
        (0,6),
        (0,0),
        (3,0),
        (4,1),
        (6,1),
        (8,0),
        (12,0),
        (13,2),
        (8,2),
        (8,4),
        (11,4),
        (11,6),
        (6,6),
        (4,3),
        (2,6)
    ]
    mesh = GLNormalMesh(points)
    faces = polygon2faces(points, Triangle{Int})
    @test eltype(faces) == Triangle{Int}
    @test facetype(mesh) == GLTriangle

end

@testset "area-2d" begin
    points = Point2f0[
        (0,0),
        (1,0),
        (1,1),
        (0,1)
    ]
    @test area(points) ≈ 1f0
    @test area(reverse(points)) ≈ -1f0
end

@testset "area-2d-nonconvex" begin
    points = Point2f0[
        (0,0),
        (1,0),
        (0.5,0.5),
        (1,1),
        (0,1),
        (0.5,0.5)
    ]
    @test area(points) ≈ 0.5f0
    @test area(reverse(points)) ≈ -0.5f0
end

@testset "area-3d" begin
    points = Point3f0[
        (0,0,0),
        (1,0,0),
        (1,1,0),
        (0,1,0)
    ]
    @test area(points) ≈ 1f0
    @test area(reverse(points)) ≈ 1f0
end

@testset "area-3d-nonconvex" begin
    points = Point3f0[
        (0,0,0),
        (1,0,0),
        (0.5,0.5,0),
        (1,1,0),
        (0,1,0),
        (0.5,0.5,0)
    ]
    @test area(points) ≈ 0.5f0
    @test area(reverse(points)) ≈ 0.5f0
end

end


text:
"""
Author: Wayne Wang
Email: wayneyw@umich.edu
Last modified: 08/26/2021
"""

module TributaryPCA

using MPI
using LinearAlgebra
using Distributions
using Random
using SparseArrays
using FileIO

include("utils_mpi.jl")
include("oja_mpi.jl")
include("oja.jl")

end

max_stars_repo_name: ayanbiswas/PRISM

max_issues_repo_name: ayanbiswas/PRISM

max_forks_repo_name: ayanbiswas/PRISM


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Private functions and macros.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Apply an alignment to a LaTeX table cell.
function _latex_apply_cell_alignment(str, alignment, j, num_printed_cols,
                                     show_row_number, vlines, left_vline,
                                     mid_vline, right_vline)

    a = _latex_alignment(alignment)
    aux_j = show_row_number ? j+1 : j

    # We only need to add left vertical line if it is the first
    # column.
    lvline = (0 ∈ vlines) && (aux_j-1 == 0) ? left_vline : ""

    # For the right vertical line, we must check if it is a mid line
    # or right line.
    if aux_j ∈ vlines
        rvline = (j == num_printed_cols) ? right_vline : mid_vline
    else
        rvline = ""
    end

    # Wrap the data into the multicolumn environment.
    return _latex_envs(str, "multicolumn{1}{$(lvline)$(a)$(rvline)}")
end

# Convert the alignment symbol into a LaTeX alignment string.
function _latex_alignment(s::Symbol)
    if (s == :l) || (s == :L)
        return "l"
    elseif (s == :c) || (s == :C)
        return "c"
    elseif (s == :r) || (s == :R)
        return "r"
    else
        error("Invalid LaTeX alignment symbol: $s.")
    end
end

# Get the LaTeX table description (alignment and vertical columns).
function _latex_table_desc(id_cols::AbstractVector,
                           alignment::Vector{Symbol},
                           show_row_number::Bool,
                           row_number_alignment::Symbol,
                           vlines::AbstractVector,
                           left_vline::AbstractString,
                           mid_vline::AbstractString,
                           right_vline::AbstractString)

    Δc = show_row_number ? 1 : 0

    str = "{"

    0 ∈ vlines && (str *= left_vline)

    # Add the alignment information of the row number column if required.
    Δc = 0
    if show_row_number
        Δc = 1
        str *= _latex_alignment(row_number_alignment)
        1 ∈ vlines && (str *= mid_vline)
    end

    # Process the alignment of all the other columns.
    for i = 1:length(id_cols)
        str *= _latex_alignment(alignment[id_cols[i]])

        if i+Δc ∈ vlines
            if i != length(id_cols)
                str *= mid_vline
            else
                str *= right_vline
            end
        end
    end

    str *= "}"

    return str
end

# Wrap the `text` into LaTeX environment(s).
_latex_envs(text::AbstractString, envs::Vector{String}) = _latex_envs(text, envs, length(envs))
_latex_envs(text::AbstractString, env::String) = "\\" * string(env) * "{" * text * "}"

function _latex_envs(text::AbstractString, envs::Vector{String}, i::Int)
    @inbounds if i > 0
        str = _latex_envs(text, envs[i])
        return _latex_envs(str, envs, i -= 1)
    end

    return text
end

# This was adapted from Julia `escape_string` function. In case of LaTeX, we
# should not escape the sequence `\\`, i.e. the backslash that is already
# escaped.
function _str_latex_escaped(io::IO, s::AbstractString, esc="")
    a = Iterators.Stateful(s)
    for c in a
        if c in esc
            print(io, '\\', c)
        elseif isascii(c)
            c == '\0'          ? print(io, escape_nul(peek(a))) :
            c == '\e'          ? print(io, "\\e") :
          # c == '\\'          ? print(io, "\\\\") :
            '\a' <= c <= '\r'  ? print(io, '\\', "abtnvfr"[Int(c)-6]) :
            isprint(c)         ? print(io, c) :
                                 print(io, "\\x", string(UInt32(c), base = 16, pad = 2))
        elseif !isoverlong(c) && !ismalformed(c)
            isprint(c)         ? print(io, c) :
            c <= '\x7f'        ? print(io, "\\x", string(UInt32(c), base = 16, pad = 2)) :
            c <= '\uffff'      ? print(io, "\\u", string(UInt32(c), base = 16, pad = need_full_hex(peek(a)) ? 4 : 2)) :
                                 print(io, "\\U", string(UInt32(c), base = 16, pad = need_full_hex(peek(a)) ? 8 : 4))
        else # malformed or overlong
            u = bswap(reinterpret(UInt32, c))
            while true
                print(io, "\\x", string(u % UInt8, base = 16, pad = 2))
                (u >>= 8) == 0 && break
            end
        end
    end
end

_str_latex_escaped(s::AbstractString, esc="") =
    sprint(_str_latex_escaped, s, esc, sizehint=lastindex(s))


text:
using MixedModels
using Test

using MixedModels: dataset

const gfms = Dict(
    :cbpp => [@formula((incid/hsz) ~ 1 + period + (1|herd))],
    :contra => [@formula(use ~ 1+age+abs2(age)+urban+livch+(1|urbdist))],
    :grouseticks => [@formula(ticks ~ 1+year+ch+ (1|index) + (1|brood) + (1|location))],
    :verbagg => [@formula(r2 ~ 1+anger+gender+btype+situ+(1|subj)+(1|item))],
)

@testset "contra" begin
    contra = dataset(:contra)
    gm0 = fit(MixedModel, only(gfms[:contra]), contra, Bernoulli(), fast=true);
    @test gm0.lowerbd == zeros(1)
    @test isapprox(gm0.θ, [0.5720734451352923], atol=0.001)
    @test isapprox(deviance(gm0), 2361.657188518064, atol=0.001)
    gm1 = fit(MixedModel, only(gfms[:contra]), contra, Bernoulli());
    @test isapprox(gm1.θ, [0.573054], atol=0.005)
    @test lowerbd(gm1) == vcat(fill(-Inf, 7), 0.)
    @test isapprox(deviance(gm1), 2361.54575, rtol=0.00001)
    @test isapprox(loglikelihood(gm1), -1180.77288, rtol=0.00001)
    @test dof(gm0) == length(gm0.β) + length(gm0.θ)
    @test nobs(gm0) == 1934
    fit!(gm0, fast=true, nAGQ=7)
    @test isapprox(deviance(gm0), 2360.9838, atol=0.001)
    gm1 = fit(MixedModel, only(gfms[:contra]), contra, Bernoulli(), nAGQ=7)
    @test isapprox(deviance(gm1), 2360.8760, atol=0.001)
    @test gm1.β == gm1.beta
    @test gm1.θ == gm1.theta
    @test isnan(gm1.σ)
    @test length(gm1.y) == size(gm1.X, 1)
    @test :θ in propertynames(gm0)
    @testset "GLMM rePCA" begin
        @test length(MixedModels.PCA(gm0)) == 1
        @test length(MixedModels.rePCA(gm0)) == 1
        @test length(gm0.rePCA) == 1
    end
    # the next three values are not well defined in the optimization
    #@test isapprox(logdet(gm1), 75.7217, atol=0.1)
    #@test isapprox(sum(abs2, gm1.u[1]), 48.4747, atol=0.1)
    #@test isapprox(sum(gm1.resp.devresid), 2237.349, atol=0.1)
    show(IOBuffer(), gm1)
    show(IOBuffer(), BlockDescription(gm0))
end

@testset "cbpp" begin
    cbpp = dataset(:cbpp)
    gm2 = fit(MixedModel, only(gfms[:cbpp]), cbpp, Binomial(), wts=float(cbpp.hsz))
    @test deviance(gm2,true) ≈ 100.09585619892968 atol=0.0001
    @test sum(abs2, gm2.u[1]) ≈ 9.723054788538546 atol=0.0001
    @test logdet(gm2) ≈ 16.90105378801136 atol=0.0001
    @test isapprox(sum(gm2.resp.devresid), 73.47174762237978, atol=0.001)
    @test isapprox(loglikelihood(gm2), -92.02628186840045, atol=0.001)
    @test isnan(sdest(gm2))
    @test varest(gm2) == 1
end

@testset "verbagg" begin
    gm3 = fit(MixedModel, only(gfms[:verbagg]), dataset(:verbagg), Bernoulli())
    @test deviance(gm3) ≈ 8151.40 rtol=1e-5
    @test lowerbd(gm3) == vcat(fill(-Inf, 6), zeros(2))
    @test fitted(gm3) == predict(gm3)
    # these two values are not well defined at the optimum
    @test isapprox(sum(x -> sum(abs2, x), gm3.u), 273.29646346940785, rtol=1e-3)
    @test sum(gm3.resp.devresid) ≈ 7156.550941446312 rtol=1e-4
end

@testset "grouseticks" begin
    center(v::AbstractVector) = v .- (sum(v) / length(v))
    grouseticks = dataset(:grouseticks)
    grouseticks.ch = center(grouseticks.height)
    gm4 = fit(MixedModel, only(gfms[:grouseticks]), grouseticks, Poisson(), fast=true)  # fails in pirls! with fast=false
    @test isapprox(deviance(gm4), 851.4046, atol=0.001)
    # these two values are not well defined at the optimum
    #@test isapprox(sum(x -> sum(abs2, x), gm4.u), 196.8695297987013, atol=0.1)
    #@test isapprox(sum(gm4.resp.devresid), 220.92685781326136, atol=0.1)
end

@testset "goldstein" begin # from a 2020-04-22 msg by Ben Goldstein to R-SIG-Mixed-Models
    goldstein =
        categorical!(
            DataFrame(
                group = repeat(1:10, outer=10),
                y = [
                    83, 3, 8, 78, 901, 21, 4, 1, 1, 39,
                    82, 3, 2, 82, 874, 18, 5, 1, 3, 50,
                    87, 7, 3, 67, 914, 18, 0, 1, 1, 38,
                    86, 13, 5, 65, 913, 13, 2, 0, 0, 48,
                    90, 5, 5, 71, 886, 19, 3, 0, 2, 32,
                    96, 1, 1, 87, 860, 21, 3, 0, 1, 54,
                    83, 2, 4, 70, 874, 19, 5, 0, 4, 36,
                    100, 11, 3, 71, 950, 21, 6, 0, 1, 40,
                    89, 5, 5, 73, 859, 29, 3, 0, 2, 38,
                    78, 13, 6, 100, 852, 24, 5, 0, 1, 39
                    ],
                ),
            :group,
        )
    gform = @formula(y ~ 1 + (1|group))
    m1 = fit(MixedModel, gform, goldstein, Poisson())
    @test deviance(m1) ≈ 193.5587302384811 rtol=1.e-5
    @test only(m1.β) ≈ 4.192196439077657 atol=1.e-5
    @test only(m1.θ) ≈ 1.838245201739852 atol=1.e-5
    m11 = fit(MixedModel, gform, goldstein, Poisson(), nAGQ=11)
    @test deviance(m11) ≈ 193.51028088736842 rtol=1.e-5
    @test only(m11.β) ≈ 4.192196439077657 atol=1.e-5
    @test only(m11.θ) ≈ 1.838245201739852 atol=1.e-5
end


text:
using Interpolations, Test

@testset "MixedTests" begin
    A = rand(6,5)
    knots = (range(1, stop=size(A,1), length=size(A,1)), range(1, stop=size(A,2), length=size(A,2)))
    itp = @inferred(interpolate(knots, A, (Gridded(Linear()),NoInterp())))
    @inferred(itp(2, 2))
    @inferred(itp(CartesianIndex((2,2))))
    for j = 2:size(A,2)-1, i = 2:size(A,1)-1
        @test itp(i,j) ≈ A[i,j]
        @test itp(CartesianIndex(i,j)) ≈ A[i,j]
    end
    @inferred(itp(knots...)) ≈ A

    A = rand(8,20)
    knots = ([x^2 for x = 1:8], [0.2y for y = 1:20])
    itp = interpolate(knots, A, Gridded(Linear()))
    @test itp(4,1.2) ≈ A[2,6]

    @test_throws  ErrorException interpolate(knots, A, (Gridded(Linear()),NoInterp()))
end


text:
using Statistics
using LinearAlgebra

coord = [0.0 0.0
         0.34 0.76
         0.98 0.23]

ends = [1 2 0.01
        2 3 0.01]

A,xc,yc,Ix,Iy,Ixy,theta,I1,I2,J,xs,ys,Cw,B1,B2,wn = cutwp_prop2(coord,ends)


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests of the Two Body orbit propagator.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 2-4. Solving Kepler's Problem [1, p. 94-95].
#
#   Initial position:
#
#       r0_ijk = + 1131.340 i - 2282.343 j + 6672.423 k [km]
#       v0_ijk = - 5.64305  i + 4.30333  j + 2.42879  k [km/s]
#
#    After 40 min., the solution of Kepler's problem leads to the following
#    position:
#
#       rf_ijk = - 4219.7527 i + 4363.0292 j - 3958.7666 k [km]
#       vf_ijk = - 3.689866  i - 1.916735  j - 6.112511  k [km/s]
#
#   Using the SatelliteToolbox Two Body orbit propagator, the following result
#   was obtained:
#
#       julia> orb = rv_to_kepler(1131340, -2282343, 6672423, -5643.05, 4303.33, 2428.79)
#       julia> orbp = init_orbit_propagator(Val{:twobody}, orb)
#       julia> (o,r,v) = step!(orbp,40*60)
#       julia> @printf("rf_ijk = %+f i %+f j %+f k", r[1]/1000, r[2]/1000, r[3]/1000)
#           rf_ijk = -4219.752738 i +4363.029177 j -3958.766617 k
#       julia> @printf("vf_ijk = %+f i %+f j %+f k", v[1]/1000, v[2]/1000, v[3]/1000)
#           vf_ijk = +3.689866 i -1.916735 j -6.112511 k
#
#   Where it is possible to see that the solution is exactly the same as
#   presented in [1].
#
################################################################################

@testset "Two Body orbit propagator" begin
    aux     = rv_to_kepler(1131340., -2282343., 6672423., -5643.05, 4303.33, 2428.79)
    orb     = Orbit(DatetoJD(1986,6,19,18,35,0), aux.a, aux.e, aux.i, aux.Ω,
                    aux.ω, aux.f)
    orbp    = init_orbit_propagator(Val{:twobody}, orb)
    (o,r,v) = step!(orbp,40*60)

    # Testing position.
    @test r[1]/1000 ≈ -4219.7527 atol=1e-3
    @test r[2]/1000 ≈ +4363.0292 atol=1e-3
    @test r[3]/1000 ≈ -3958.7666 atol=1e-3

    # Testing velocity.
    @test v[1]/1000 ≈ +3.689866 atol=1e-6
    @test v[2]/1000 ≈ -1.916735 atol=1e-6
    @test v[3]/1000 ≈ -6.112511 atol=1e-6

    # Test using the function `propagate!`
    # ==========================================================================

    (o,r,v) = propagate!(orbp, 40*60)

    # Testing position.
    @test r[1]/1000 ≈ -4219.7527 atol=1e-3
    @test r[2]/1000 ≈ +4363.0292 atol=1e-3
    @test r[3]/1000 ≈ -3958.7666 atol=1e-3

    # Testing velocity.
    @test v[1]/1000 ≈ +3.689866 atol=1e-6
    @test v[2]/1000 ≈ -1.916735 atol=1e-6
    @test v[3]/1000 ≈ -6.112511 atol=1e-6

    # Test using the function `propagate_to_epoch!`
    # ==========================================================================

    #                                        1986/6/19 18:35:00 + 40 min.
    (o,r,v) = propagate_to_epoch!(orbp, DatetoJD(1986,6,19,19,15,0))

    # Testing position.
    @test r[1]/1000 ≈ -4219.7527 atol=1e-3
    @test r[2]/1000 ≈ +4363.0292 atol=1e-3
    @test r[3]/1000 ≈ -3958.7666 atol=1e-3

    # Testing velocity.
    @test v[1]/1000 ≈ +3.689866 atol=1e-6
    @test v[2]/1000 ≈ -1.916735 atol=1e-6
    @test v[3]/1000 ≈ -6.112511 atol=1e-6
end


text:
using AdventOfCodeSolutions
using Test

function input(puzzle::Puzzle{2020, 2, n}) where n
    io = openInput(puzzle)
    return split(read(io, String), '\n', keepempty=false)
end

function parsePassword(s)
    range, char, password = split(s, ' ')
    left, right = parse.(Int, split(range, '-'))

    return (left:right, char[1], password)
end

function isValidPart1(s)
    range, char, password = parsePassword(s)

    return count(==(char), password) in range
end

solve(::Puzzle{2020, 2, 1}, input) = count(isValidPart1, input)

function isValidPart2(s)
    range, char, password = parsePassword(s)

    return xor(
        password[range.start] == char,
        password[range.stop] == char
    )
end

solve(::Puzzle{2020, 2, 2}, input) = count(isValidPart2, input)

@testset "Day 2" begin
    @test(solve(Puzzle(2, 1), input(Puzzle(2, 1))) == 556)
    @test(solve(Puzzle(2, 2), input(Puzzle(2, 2))) == 605)
end


text:
using Gadfly

Gadfly.set_default_plot_format(:png)

#Captures figures
function Base.display(report::Report, m::MIME"image/png", p::Plot)
    chunk = report.cur_chunk

    full_name, rel_name = get_figname(report, chunk)

    docformat = formats[report.formatdict[:doctype]]

    #Add to results for term chunks and store otherwise
    if chunk[:term]
      chunk[:figure] = [rel_name]

      if report.term_state == :text
        report.cur_result *= "\n" * report.formatdict[:codeend] * "\n"
      end

      report.cur_result *= formatfigures(chunk, docformat)
      report.term_state = :fig
      chunk[:figure] = String[]
    else
      push!(report.figures, rel_name)
    end

    report.fignum += 1

    w = chunk[:fig_width]inch
    h = chunk[:fig_height]inch
    format = chunk[:fig_ext]
    dpi = chunk[:dpi]

    #This is probably not the correct way to handle different formats, but it works.
    if format == ".png"
        try
          draw(PNG(full_name, w, h, dpi=dpi), p)
        catch
          draw(PNG(full_name, w, h), p) #Compose < 0.3.1, Gadfly < 0.3.1
        end
    elseif format == ".pdf"
        draw(PDF(full_name, w, h), p)
    elseif format == ".ps"
        draw(PS(full_name, w, h), p)
    elseif format == ".svg"
        draw(SVG(full_name, w, h), p)
    elseif format == ".js.svg"
        draw(SVGJS(full_name, w, h), p)
    else:
        warn("Can't save figure. Unsupported format")
    end
end


text:
let
    struct TigerPOMDPTestFixture <: POMDP{Bool, Int, Bool} end
    POMDPs.states(::TigerPOMDPTestFixture) = (true, false)
    POMDPs.stateindex(::TigerPOMDPTestFixture, s) = Int(s) + 1
    POMDPs.actions(m::TigerPOMDPTestFixture) = 0:2
    POMDPs.actionindex(m::TigerPOMDPTestFixture, s::Int) = s+1
    POMDPs.observations(::TigerPOMDPTestFixture) = (true, false)
    POMDPs.obsindex(::TigerPOMDPTestFixture, o) = Int(o) + 1

    pomdp = TigerPOMDPTestFixture()

    @test ordered_states(pomdp) == [false, true]
    @test ordered_observations(pomdp) == [false, true]
    @test ordered_actions(pomdp) == [0,1,2]
end

struct TM <: POMDP{Int, Int, Int} end
POMDPs.states(::TM) = [1,3]
POMDPs.stateindex(::TM, s::Int) = s

@test_throws ErrorException ordered_states(TM())


text:
push!(LOAD_PATH,"../src/")

using Documenter
using DataIO

makedocs(
    sitename = "DataIO.jl",
    authors = "Tobias Frilling",
    format = Documenter.HTML(),
    modules = [DataIO],
    pages = [
        "Home" => "index.md",
        "File Formats" => [
            "Data `(*.lrn)`" => "lrn.md",
            "Matrix `(*.umx)`" => "umx.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/ckafi/DataIO.jl"
)

max_stars_repo_name: ckafi/DataIO

max_issues_repo_name: ckafi/DataIO

max_forks_repo_name: ckafi/DataIO


text:
module GPUMatFac

include("model.jl")
include("losses.jl")
include("regularizer.jl")
include("fit.jl")
include("impute.jl")
include("transform.jl")
include("model_io.jl")

BLAS.set_num_threads(1)

end # module


text:
struct MonteCarlo <: AbstractMonteCarlo
    n::Integer
    MonteCarlo(n) = n > 0 ? new(n) : error("n must be greater than zero")
end

struct SobolSampling <: AbstractQuasiMonteCarlo
    n::Integer
    SobolSampling(n) = n > 0 ? new(n) : error("n must be greater than zero")
end

struct HaltonSampling <: AbstractQuasiMonteCarlo
    n::Integer
    HaltonSampling(n) = n > 0 ? new(n) : error("n must be greater than zero")
end

struct LatinHypercubeSampling <: AbstractQuasiMonteCarlo
    n::Integer
    LatinHypercubeSampling(n) = n > 0 ? new(n) : error("n must be greater than zero")
end

struct LatticeRuleSampling <: AbstractQuasiMonteCarlo
    n::Integer
    LatticeRuleSampling(n) = n > 0 ? new(n) : error("n must be greater than zero")
end

function sample(inputs::Array{<:UQInput}, sim::MonteCarlo)
    return sample(inputs, sim.n)
end

function sample(inputs::Array{<:UQInput}, sim::AbstractQuasiMonteCarlo)
    random_inputs = filter(i -> isa(i, RandomUQInput), inputs)
    deterministic_inputs = filter(i -> isa(i, DeterministicUQInput), inputs)

    n_rv = count_rvs(random_inputs)

    u = qmc_samples(sim, n_rv)

    samples = quantile.(Normal(), u)
    samples = DataFrame(names(random_inputs) .=> eachcol(samples))

    if !isempty(deterministic_inputs)
        samples = hcat(samples, sample(deterministic_inputs, sim.n))
    end

    to_physical_space!(inputs, samples)

    return samples
end

function qmc_samples(sim::SobolSampling, rvs::Integer)
    return transpose(QuasiMonteCarlo.sample(sim.n, zeros(rvs), ones(rvs), SobolSample()))
end

function qmc_samples(sim::HaltonSampling, rvs::Integer)
    coprimes = [prime(i) for i in 1:rvs]
    return transpose(
        QuasiMonteCarlo.sample(sim.n, zeros(rvs), ones(rvs), LowDiscrepancySample(coprimes))
    )
end

function qmc_samples(sim::LatinHypercubeSampling, rvs::Integer)
    return transpose(
        QuasiMonteCarlo.sample(sim.n, zeros(rvs), ones(rvs), LatinHypercubeSample())
    )
end

function qmc_samples(sim::LatticeRuleSampling, rvs::Integer)
    return transpose(
        QuasiMonteCarlo.sample(sim.n, zeros(rvs), ones(rvs), LatticeRuleSample())
    )
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
# ---
# title: 1472. Design Browser History
# id: problem1472
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Design
# link: <https://leetcode.com/problems/design-browser-history/description/>
# hidden: true
# ---
#
# You have a **browser** of one tab where you start on the `homepage` and you
# can visit another `url`, get back in the history number of `steps` or move
# forward in the history number of `steps`.
#
# Implement the `BrowserHistory` class:
#
#   * `BrowserHistory(string homepage)` Initializes the object with the `homepage` of the browser.
#   * `void visit(string url)` Visits `url` from the current page. It clears up all the forward history.
#   * `string back(int steps)` Move `steps` back in history. If you can only return `x` steps in the history and `steps > x`, you will return only `x` steps. Return the current `url` after moving back in history **at most** `steps`.
#   * `string forward(int steps)` Move `steps` forward in history. If you can only forward `x` steps in the history and `steps > x`, you will forward only `x` steps. Return the current `url` after forwarding in history **at most** `steps`.
#
#
#
# **Example:**
#
#
#
#     Input:
#     ["BrowserHistory","visit","visit","visit","back","back","forward","visit","forward","back","back"]
#     [["leetcode.com"],["google.com"],["facebook.com"],["youtube.com"],[1],[1],[1],["linkedin.com"],[2],[2],[7]]
#     Output:
#     [null,null,null,null,"facebook.com","google.com","facebook.com",null,"linkedin.com","google.com","leetcode.com"]
#
#     Explanation:
#     BrowserHistory browserHistory = new BrowserHistory("leetcode.com");
#     browserHistory.visit("google.com");       // You are in "leetcode.com". Visit "google.com"
#     browserHistory.visit("facebook.com");     // You are in "google.com". Visit "facebook.com"
#     browserHistory.visit("youtube.com");      // You are in "facebook.com". Visit "youtube.com"
#     browserHistory.back(1);                   // You are in "youtube.com", move back to "facebook.com" return "facebook.com"
#     browserHistory.back(1);                   // You are in "facebook.com", move back to "google.com" return "google.com"
#     browserHistory.forward(1);                // You are in "google.com", move forward to "facebook.com" return "facebook.com"
#     browserHistory.visit("linkedin.com");     // You are in "facebook.com". Visit "linkedin.com"
#     browserHistory.forward(2);                // You are in "linkedin.com", you cannot move forward any steps.
#     browserHistory.back(2);                   // You are in "linkedin.com", move back two steps to "facebook.com" then to "google.com". return "google.com"
#     browserHistory.back(7);                   // You are in "google.com", you can move back only one step to "leetcode.com". return "leetcode.com"
#
#
#
#
# **Constraints:**
#
#   * `1 <= homepage.length <= 20`
#   * `1 <= url.length <= 20`
#   * `1 <= steps <= 100`
#   * `homepage` and `url` consist of  '.' or lower case English letters.
#   * At most `5000` calls will be made to `visit`, `back`, and `forward`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@testset "hyperparameters" begin
    empty!(HYPERPARAMETERS)

    withenv("SM_HP_HPA" => "1", "SM_HP_HPB" => "2", "OTHER_HPC" => "3.0", "OTHER_HPD" => "hi") do
        @test HYPERPARAMETERS == Dict()

        hp = hyperparam(Float64, :hpa)
        @test hp == 1.0
        @test hp isa Float64
        @test HYPERPARAMETERS == Dict(:hpa => 1.0)
        @test HYPERPARAMETERS[:hpa] isa Float64

        hp = hyperparam(String, :hpb)
        @test hp == "2"
        @test HYPERPARAMETERS == Dict(:hpa => 1.0, :hpb => "2")

        # Log if type changes
        @test_log LOGGER "notice" "type of HYPERPARAMETERS[:hpa]" hyperparam(Int, :hpa)
        # No log if value is unchanged
        @test_nolog LOGGER "notice" "HYPERPARAMETERS[:hpa]" hyperparam(Int, :hpa)
        hp = HYPERPARAMETERS[:hpa]
        @test hp == 1
        @test hp isa Int
        @test HYPERPARAMETERS == Dict(:hpa => 1, :hpb => "2")

        hp = hyperparam(String, :hpc; prefix="OTHER_")
        @test hp == "3.0"
        @test HYPERPARAMETERS == Dict(:hpa => 1, :hpb => "2", :hpc => "3.0")
        # Log if old value != new value
        @test_log LOGGER "notice" "Overwriting HYPERPARAMETERS[:hpc]" hyperparam(Float64, :hpc; prefix="OTHER_")
        @test HYPERPARAMETERS == Dict(:hpa => 1, :hpb => "2", :hpc => 3.0)

        # No Key
        @test_throws KeyError hyperparam(:hpa; prefix="OTHER_")
        # Parse error
        @test_throws ArgumentError hyperparam(Float64, :hpd; prefix="OTHER_")
        @test HYPERPARAMETERS == Dict(:hpa => 1, :hpb => "2", :hpc => 3)

        empty!(HYPERPARAMETERS)
        @test HYPERPARAMETERS == Dict()

        hp = hyperparams(:hpc; prefix="OTHER_")
        @test hp == (hpc=3,)
        @test hp.hpc isa Float64
        @test HYPERPARAMETERS == Dict(:hpc => 3.0)
        @test HYPERPARAMETERS[:hpc] isa Float64

        @test_throws ArgumentError hyperparams(:hpa, :hpb; types=[Int, String, String])
        @test HYPERPARAMETERS == Dict(:hpc => 3)

        hp = hyperparams(:hpa, :hpb; types=[Int, String])
        @test hp == (hpa=1, hpb="2")
        @test hp.hpa isa Int
        @test HYPERPARAMETERS == Dict(:hpa => 1, :hpb => "2", :hpc => 3)
        @test HYPERPARAMETERS[:hpa] isa Int

        save_hyperparam(:hpd, 4, prefix="OTHER_")
        @test HYPERPARAMETERS[:hpd] isa Int
        @test HYPERPARAMETERS[:hpd] == 4
        @test ENV["OTHER_HPD"] == "4"
    end

    @testset "report_hyperparameters" begin
        # Ensure hyperparameters are what we expect
        empty!(HYPERPARAMETERS)

        withenv("SM_HP_HPA" => "1", "SM_HP_HPB" => "2", "SM_HP_HPC" => "3") do
            # we use hpa and hpb, but we want hpc to just be scoped from environment
            hyperparam(Float64, :hpa)
            hyperparam(String, :hpb)
            # It will have to guess the type of hpc

            mktmpdir() do dir
                @test_log LOGGER "info" "hyperparameters: hpb=2" report_hyperparameters(dir)

                contents = read(joinpath(dir, "hyperparameters.json"), String)
                @test occursin("\"hpa\": 1.0", contents)
                @test occursin("\"hpb\": \"2\"", contents)
                @test occursin("\"hpc\": 3", contents)
            end
        end
    end

    @testset "_parse_hyper: guessing type correctly" begin
        _parse_hyper = Hyperparameters._parse_hyper

        @test _parse_hyper("true") isa Bool
        @test _parse_hyper("true") == true
        @test _parse_hyper("false") isa Bool
        @test _parse_hyper("false") == false

        @test _parse_hyper("1") isa Int
        @test _parse_hyper("1") == 1
        @test _parse_hyper("-1") isa Int
        @test _parse_hyper("-1") == -1

        @test _parse_hyper("2.0") isa Float64
        @test _parse_hyper("2.0") == 2.0
        @test _parse_hyper("-2.0") isa Float64
        @test _parse_hyper("-2.0") == -2.0

        @test _parse_hyper("three") isa String
        @test _parse_hyper("three") == "three"

        # if not matching julia lowercase convention get strings:
        @test _parse_hyper("TRUE") isa String
        @test _parse_hyper("True") isa String
    end
end


text:
# options.jl: Wrapper for internal DBoptlist objects, with Julia support

#TODO: define
# - convert(Ptr{DBoptlist}, silooptions)
# - convert(silooptions, Ptr{DBoptlist})
# - define push!/pop!/setindex!/getindex for silooptions