text:
### This will have profile tests for the ape budgets

@testset "APE budgets" begin
    sx,sy,sz,st = 100,100,100,10
    x = y = range(0.0,stop=10.0,length = sx);
    z = range(0.0,stop=500.0,length = sz);
    t = range(0.0,stop=50.0,length = st);
    dx = dy = x[2] - x[1];
    dz = z[2] - z[1];
    dt = t[2] - t[1];

    zero_profile = zeros(sx,sy,sz,st);
    one_profile = ones(sx,sy,sz,st);
    test_brunt = 1e-4*ones(sz,st)
    test_rho = ones(1,1,sz,st)
    test_surface = zeros(sx,sy,st)

    @test  getapebudget(zero_profile,zero_profile,zero_profile,zero_profile,test_brunt,
                    zero_profile,test_surface,zero_profile,test_rho,x,y,z,t,dx,dy,dz,dt,50.0)[3] ≈ zeros(st)

    ### If B is constant in space, advection should be zero
    B = 5ones(sx,sy,sz,st);
    for i in 1:st
        for j in 1:sy
            for k in 1:sx
                B[j,k,:,i] .= B[j,k,:,i] .* sin(2pi/st*i)
            end
        end
    end
    RAD_b = 2ones(sx,sy,sz,st);
    Fs = 2ones(sx,sy,st);
    U = 10ones(sx,sy,sz,st);
    V = 10ones(sx,sy,sz,st);
    W = 10ones(sx,sy,sz,st);
    N2 = 1e-4 .* ones(1,1,sz,st);
    rho0 = ones(1,1,sz,st);
    diabatic_other = get_diabatic_as_residual_buoyancy(B, RAD_b, Fs, U,V, W, N2, dx,dy, dz, dt);
    a = getapebudget(B,U,V,W,N2,RAD_b,Fs,diabatic_other, rho0, x,y,z,t,dx,dy,dz,dt,500)
    @test a[4] ≈ -(a[7] - a[8] - a[9] - a[10])

    ### Budget must be closed when we have Advection, Rad, and Surface fluxes
    B = 5ones(sx,sy,sz,st);
    for i in 1:st
        for j in 1:sy
            for k in 1:sx
                B[j,k,:,i] .= B[j,k,:,i] .* sin(2pi/st*i)*sin(2pi/sx*k)*sin(2pi/sy*j)
            end
        end
    end
    RAD_b = 2ones(sx,sy,sz,st);
    Fs = 2ones(sx,sy,st);
    U = 10rand(sx,sy,sz,st);
    V = 10rand(sx,sy,sz,st);
    W = 5rand(sx,sy,sz,st);
    N2 = 1e-4 .* ones(1,1,sz,st);
    rho0 = ones(1,1,sz,st);
    diabatic_other = get_diabatic_as_residual_buoyancy(B, RAD_b, Fs, U,V, W, N2, dx,dy, dz, dt);
    a = getapebudget(B,U,V,W,N2,RAD_b,Fs,diabatic_other, rho0, x,y,z,t,dx,dy,dz,dt,500)
    @test a[4] ≈ -(a[5] + a[6] + a[7] - a[8] - a[9] - a[10])

    ### Budget must be closed when we have Advection, Rad, and zero surface fluxes
    B = 5ones(sx,sy,sz,st);
    for i in 1:st
        for j in 1:sy
            for k in 1:sx
                B[j,k,:,i] .= B[j,k,:,i] .* sin(2pi/st*i)*sin(2pi/sx*k)*sin(2pi/sy*j)
            end
        end
    end
    RAD_b = 2ones(sx,sy,sz,st);
    Fs = zeros(sx,sy,st);
    U = 10rand(sx,sy,sz,st);
    V = 10rand(sx,sy,sz,st);
    W = 5rand(sx,sy,sz,st);
    N2 = 1e-4 .* ones(1,1,sz,st);
    rho0 = ones(1,1,sz,st);
    diabatic_other = get_diabatic_as_residual_buoyancy(B, RAD_b, Fs, U,V, W, N2, dx,dy, dz, dt);
    a = getapebudget(B,U,V,W,N2,RAD_b,Fs,diabatic_other, rho0, x,y,z,t,dx,dy,dz,dt,500)
    @test a[4] ≈ -(a[5] + a[6] + a[7] - a[8] - a[9] )

    ### Budget must be closed when we have Advection, zero Rad, and finite surface fluxes
    B = 5ones(sx,sy,sz,st);
    for i in 1:st
        for j in 1:sy
            for k in 1:sx
                B[j,k,:,i] .= B[j,k,:,i] .* sin(2pi/st*i)*sin(2pi/sx*k)*sin(2pi/sy*j)
            end
        end
    end
    RAD_b = zeros(sx,sy,sz,st);
    Fs = 2ones(sx,sy,st);
    U = 10rand(sx,sy,sz,st);
    V = 10rand(sx,sy,sz,st);
    W = 5rand(sx,sy,sz,st);
    N2 = 1e-4 .* ones(1,1,sz,st);
    rho0 = ones(1,1,sz,st);
    diabatic_other = get_diabatic_as_residual_buoyancy(B, RAD_b, Fs, U,V, W, N2, dx,dy, dz, dt);
    a = getapebudget(B,U,V,W,N2,RAD_b,Fs,diabatic_other, rho0, x,y,z,t,dx,dy,dz,dt,500)
    @test a[4] ≈ -(a[5] + a[6] + a[7] - a[9] - a[10])

   #  labels = ["mass" "KE" "APE" "∂tAPE" "AdvecU" "AdvecV" "Conversion" "Rad" "Convec" "Surf" "Res"]
   #  mycolors = distinguishable_colors(11) #NEEDS COLORS.JL

   #  fig = Figure(resolution = (1200, 700), backgroundcolor = RGBf0(0.98, 0.98, 0.98),title = "APE budget closure test")
   #  ax1 = fig[1,1] = Axis(fig)
   #  ylims!(ax1,extrema(reduce(vcat,a[4:11])))
   #  for i in 4:11
   #     lines!(ax1,a[i], label = labels[i], color = mycolors[i], linewidth = 5)
   #  end
   #  axislegend(ax1)
   # fig

end


text:
using Conda

Conda.add("thermo=0.1.40", channel="conda-forge")
# pyimport_conda("thermo=0.1.40", "thermo", "conda-forge")

ext: jl


text:
# Locally Linear Embedding (LLE)
# ------------------------
# Nonlinear dimensionality reduction by locally linear embedding,
# Roweis, S. & Saul, L., Science 290:2323 (2000)

using Logging

#### LLE type
struct LLE{T <: Real} <: AbstractDimensionalityReduction
    k::Int
    λ::AbstractVector{T}
    proj::Projection{T}

    LLE{T}(k::Int, λ::AbstractVector{T}, proj::Projection{T})  where T = new(k, λ, proj)
end

## properties
outdim(R::LLE) = size(R.proj, 1)
eigvals(R::LLE) = R.λ
neighbors(R::LLE) = R.k

## show
summary(io::IO, R::LLE) = print(io, "LLE(outdim = $(outdim(R)), neighbors = $(neighbors(R)))")

## interface functions
function fit(::Type{LLE}, X::AbstractMatrix{T};
             maxoutdim::Int=2, k::Int=12, tol::Real=1e-5,
             use_naive=false, use_eigen=false) where {T<:Real}
    # Construct NN graph
    D, E = find_nn(X, k)
    _, C = largest_component(SimpleWeightedGraph(adjmat(D,E)))
    X = X[:, C]
    n = length(C)

    # Correct indexes of neighbors if more then one connected component
    Ec = E
    if size(E,2) != n
        R = Dict(zip(C, collect(1:n)))
        Ec = zeros(Int,k,n)
        for i in 1 : n
            Ec[:,i] = map(j->get(R,j,C[i]), E[:,C[i]])
        end
    end

    if k > maxoutdim
        @warn("k > maxoutdim: regularization will be used")
    end

    # Reconstruct weights and compute embedding:
    # M = (I - w)'(I - w) = I - w'I - Iw + w'w
    M = spdiagm(0 => fill(one(T), n))
    Ones = fill(one(T), k, 1)
    for i in 1 : n
        J = Ec[:,i]
        Z = view(X, :, J) .- view(X, :, i)
        G = transpose(Z)*Z
        G += I * tol # regularize
        w = vec(G \ Ones)
        w ./= sum(w)
        ww = w*transpose(w)
        for (l, j) in enumerate(J)
            M[i,j] -= w[l]
            M[j,i] -= w[l]
            for (m, jj) in enumerate(J)
                M[j,jj] = ww[l,m]
            end
        end
    end

    if use_eigen
        λ, V = decompose(M, maxoutdim)
    else
        λ, V = specDecompose(M, tol, maxoutdim, use_naive)
    end
    return LLE{T}(k, λ, rmul!(transpose(V), sqrt(n)))
end

transform(R::LLE) = R.proj


text:
using jlmie
using Plots

# settings
# structure
nmat = 4  # refractive index of a sphere
radius = 100*1e-9

# environment
nenv = 1.00

# wavelenth
lbd0 = (400:1000)*1e-9; lbdp = lbd0.*1e9;

# other settings
nmax = -1  # -1: namx large enough (determined from x)

# main
m,x = jlmie_mx(nmat,radius,lbd0,nenv)
Qsca = jlmie_Qsca(m,x,nmax)
Qsca_ED = jlmie_Qsca_a_n(m,x,1)
Qsca_MD = jlmie_Qsca_b_n(m,x,1)
Qsca_EQ = jlmie_Qsca_a_n(m,x,2)
Qsca_MQ = jlmie_Qsca_b_n(m,x,2)
Qsca_EH = jlmie_Qsca_a_n(m,x,3)
Qsca_MH = jlmie_Qsca_b_n(m,x,3)
println("calculation completed")

# show results
# validation can be done by comparing to ITMO Mie Calculator
# https://physics.itmo.ru/en/mie
plt = plot(lbdp,Qsca,
    xlabel       = "Wavelength (nm)",
    ylabel       = "Scattering efficiency",
    # legend       = false,
    label        = "Total",
    # xlims        =(-3,3),
    # ylims        =(-3,3),
    # aspect_ratio =0.5,
    # title        ="TITLE",
    # linecolor    =:blue,
    # linewidth    =5,
    # linestyle    =:dot,
    # size         =(400,300),
) # https://qiita.com/I_ppp/items/dca3552affa6a672e4bd
plot!(lbdp,Qsca_ED,
    label = "ED",
    linestyle = :dash,
)
plot!(lbdp,Qsca_MD,
    label = "MD",
    linestyle = :dash,
)
plot!(lbdp,Qsca_EQ,
    label = "EQ",
    linestyle = :dash,
)
plot!(lbdp,Qsca_MQ,
    label = "MQ",
    linestyle = :dash,
)
plot!(lbdp,Qsca_EH,
    label = "EH",
    linestyle = :dash,
)
plot!(lbdp,Qsca_MH,
    label = "MH",
    linestyle = :dash,
)
display(plt)

# save
if true
    runningfilename = splitext(splitpath(@__FILE__)[end])[1]
    outfilename = ".\\result\\" * runningfilename
    savefig(plt,outfilename)
    println("Graph has been saved in result directory")
end

# find peaks
lbdp_peaks = zeros(6)
lbdp_peaks[1] = lbdp[findmax(Qsca_ED)[2]]
lbdp_peaks[2] = lbdp[findmax(Qsca_MD)[2]]
lbdp_peaks[3] = lbdp[findmax(Qsca_EQ)[2]]
lbdp_peaks[4] = lbdp[findmax(Qsca_MQ)[2]]
lbdp_peaks[5] = lbdp[findmax(Qsca_EH)[2]]
lbdp_peaks[6] = lbdp[findmax(Qsca_MH)[2]]

println("maximum of ED is at " * string(round(lbdp_peaks[1])))
println("maximum of MD is at " * string(round(lbdp_peaks[2])))
println("maximum of EQ is at " * string(round(lbdp_peaks[3])))
println("maximum of MQ is at " * string(round(lbdp_peaks[4])))
println("maximum of EH is at " * string(round(lbdp_peaks[5])))
println("maximum of MH is at " * string(round(lbdp_peaks[6])))

max_stars_repo_name: Hinamoooon/jlmie

max_issues_repo_name: Hinamoooon/jlmie

max_forks_repo_name: Hinamoooon/jlmie


text:
MOI.Utilities.@product_of_sets(
    Equalities,
    MOI.EqualTo{T},
)

MOI.Utilities.@product_of_sets(
    Inequalities,
    MOI.LessThan{T},
)

MOI.Utilities.@struct_of_constraints_by_set_types(
    EqualitiesOrInequalities,
    MOI.EqualTo{T},
    MOI.LessThan{T},
)

const QPForm{T} = MOI.Utilities.GenericModel{
    T,
    MOI.Utilities.ObjectiveContainer{T},
    MOI.Utilities.FreeVariables,
    EqualitiesOrInequalities{T}{
        MOI.Utilities.MatrixOfConstraints{
            T,
            MOI.Utilities.MutableSparseMatrixCSC{
                T,
                Int,
                MOI.Utilities.OneBasedIndexing,
            },
            MOI.Utilities.Hyperrectangle{T},
            Equalities{T},
        },
        MOI.Utilities.MatrixOfConstraints{
            T,
            MOI.Utilities.MutableSparseMatrixCSC{
                T,
                Int,
                MOI.Utilities.OneBasedIndexing,
            },
            MOI.Utilities.Hyperrectangle{T},
            Inequalities{T},
        },
    },
}

mutable struct QPDiff <: DiffModel
    # storage for problem data in matrix form
    model::QPForm{Float64}
    # includes maps from matrix indices to problem data held in `optimizer`
    # also includes KKT matrices
    # also includes the solution
    gradient_cache::Union{Nothing,QPCache}

    # caches for sensitivity output
    # result from solving KKT/residualmap linear systems
    # this allows keeping the same `gradient_cache`
    # if only sensitivy input changes
    forw_grad_cache::Union{Nothing,QPForwBackCache}
    back_grad_cache::Union{Nothing,QPForwBackCache}

    # sensitivity input cache using MOI like sparse format
    input_cache::DiffInputCache

    x::Vector{Float64} # Primal
    λ::Vector{Float64} # Dual of inequalities
    ν::Vector{Float64} # Dual of equalities
end
function QPDiff()
    return QPDiff(QPForm{Float64}(), nothing, nothing, nothing, DiffInputCache(), Float64[], Float64[], Float64[])
end

function MOI.empty!(model::QPDiff)
    MOI.empty!(model.model)
    model.gradient_cache = nothing
    model.forw_grad_cache = nothing
    model.back_grad_cache = nothing
    empty!(model.input_cache)
    empty!(model.x)
    empty!(model.λ)
    empty!(model.ν)
    return
end

const EQ = MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},MOI.EqualTo{Float64}}
const LE = MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},MOI.LessThan{Float64}}

function MOI.set(model::QPDiff, ::MOI.ConstraintPrimalStart, ci::MOI.ConstraintIndex, value)
    MOI.throw_if_not_valid(model, ci)
    return # Ignored
end

# We want to stay consistent with the variable `ν` defined in (3) of
# Left hand side of eq. (6) in https://arxiv.org/pdf/1703.00443.pdf
# However, in eq. (6), they put it in the lagrangian as
# `+ ν ⋅ (Az - b)`
# while in MOI, we put it as
# `- ν ⋅ (Az - b)`
# so the we should reverse the sign if we want to use the same equations
# as in the paper.
function MOI.set(model::QPDiff, ::MOI.ConstraintDualStart, ci::EQ, value)
    MOI.throw_if_not_valid(model, ci)
    _enlarge_set(model.ν, MOI.Utilities.rows(_equalities(model), ci), -value)
end

function MOI.set(model::QPDiff, ::MOI.ConstraintDualStart, ci::LE, value)
    MOI.throw_if_not_valid(model, ci)
    _enlarge_set(model.λ, MOI.Utilities.rows(_inequalities(model), ci), -value)
end

function _gradient_cache(model::QPDiff)
    if model.gradient_cache !== nothing
        return model.gradient_cache
    end
    A = convert(SparseMatrixCSC{Float64, Int}, _equalities(model).coefficients)
    b = _equalities(model).constants.upper
    G = convert(SparseMatrixCSC{Float64, Int}, _inequalities(model).coefficients)
    h = _inequalities(model).constants.upper

    nz = size(A, 2)
    # TODO ideally, the objective should be stored in matrix form in `model.model`
    objective_function = MOI.get(model.model, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
    sparse_array_obj = sparse_array_representation(objective_function, nz)
    Q = sparse_array_obj.quadratic_terms
    q = sparse_array_obj.affine_terms

    LHS = create_LHS_matrix(model.x, model.λ, Q, G, h, A)

    model.gradient_cache = QPCache(LHS)

    return model.gradient_cache
end

function _equalities(model::QPDiff)
    return MOI.Utilities.constraints(
        model.model.constraints,
        MOI.ScalarAffineFunction{Float64},
        MOI.EqualTo{Float64},
    )
end
function _inequalities(model::QPDiff)
    return MOI.Utilities.constraints(
        model.model.constraints,
        MOI.ScalarAffineFunction{Float64},
        MOI.LessThan{Float64},
    )
end

# TODO: create test functions for the methods

# """
#     Left hand side of eqn(6) in https://arxiv.org/pdf/1703.00443.pdf
# """
# function create_LHS_matrix(z, λ, Q, G, h, A=nothing)
#     if A == nothing || size(A)[1] == 0
#         return [Q                G';
#                 Diagonal(λ) * G    Diagonal(G * z - h)]
#     else
#         @assert size(A)[2] == size(G)[2]
#         p, n = size(A)
#         m    = size(G)[1]
#         return [Q                  G'                    A';
#                 Diagonal(λ) * G    Diagonal(G * z - h)   zeros(m, p);
#                 A                  zeros(p, m)           zeros(p, p)]
#     end
# end

"""
    create_LHS_matrix(z, λ, Q, G, h, A=nothing)

Inverse matrix specified on RHS of eqn(7) in https://arxiv.org/pdf/1703.00443.pdf

Helper method while calling `_backward_quad`
"""
function create_LHS_matrix(z, λ, Q, G, h, A=nothing)::AbstractMatrix{Float64}
    if (A === nothing || size(A)[1] == 0) && (G === nothing || size(G)[1] == 0)
        return Q
    elseif A === nothing || size(A)[1] == 0
        return [Q         G' * Diagonal(λ);
                G         Diagonal(G * z - h)]
    elseif G === nothing || size(G)[1] == 0
        p, n = size(A)
        return [Q         A';
                A         spzeros(p, p)]
    else
        p, n  = size(A)
        m, n2 = size(G)
        if n != n2
            throw(DimensionError("Sizes of $A and $G do not match"))
        end
        return [Q         G' * Diagonal(λ)       A';
                G         Diagonal(G * z - h)    spzeros(m, p);
                A         spzeros(p, m)          spzeros(p, p)]
    end
end
# TODO: this is the transpose, check back for usage

# """
#     Right hand side of eqn(6) in https://arxiv.org/pdf/1703.00443.pdf
# """
# function create_RHS_matrix(z, dQ, dq, λ, dG, dh, ν=nothing, dA=nothing, db=nothing)
#     if dA == nothing || size(dA)[1] == 0
#         return -[dQ * z + dq + dG' * λ      ;
#                  Diagonal(λ) * (dG * z - dh)]
#     else
#         return -[dQ * z + dq + dG' * λ + dA' * ν;
#                  Diagonal(λ) * (dG * z - dh)    ;
#                  dA * z - db                    ]
#     end
# end

const _QP_SET_TYPES = Union{
    MOI.GreaterThan{Float64},
    MOI.LessThan{Float64},
    MOI.EqualTo{Float64},
    # MOI.Interval{Float64},
}

const _QP_FUNCTION_TYPES = Union{
    MOI.VariableIndex,
    MOI.ScalarAffineFunction{Float64},
}

_qp_supported(::Type{F}, ::Type{S}) where {F <: _QP_FUNCTION_TYPES, S <: _QP_SET_TYPES} = true
_qp_supported(::Type{F}, ::Type{S}) where {F, S} = false
function _qp_supported(model::MOI.AbstractOptimizer)
    return all(FS -> _qp_supported(FS...), MOI.get(model, MOI.ListOfConstraintTypesPresent()))
end

function MOI.get(model::QPDiff, ::ForwardOutVariablePrimal, vi::VI)
    return model.forw_grad_cache.dz[vi.value]
end

function _get_db(model::QPDiff, ci::LE)
    i = ci.value
    # dh = -Diagonal(λ) * dλ
    return model.λ[i] * model.back_grad_cache.dλ[i]
end
function _get_db(model::QPDiff, ci::EQ)
    return model.back_grad_cache.dν[ci.value]
end

"""
    backward(model::QPDiff)

Method to differentiate optimal solution `z` and return
product of jacobian matrices (`dz / dQ`, `dz / dq`, etc) with
the backward pass vector `dl / dz`

The method computes the product of
1. jacobian of problem solution `z*` with respect to
    problem parameters set with the [`BackwardInVariablePrimal`](@ref)
2. a backward pass vector `dl / dz`, where `l` can be a loss function

Note that this method *does not returns* the actual jacobians.

For more info refer eqn(7) and eqn(8) of https://arxiv.org/pdf/1703.00443.pdf
"""
function backward(model::QPDiff)
    gradient_cache = _gradient_cache(model)
    LHS = gradient_cache.lhs

    neq = length(model.ν)
    nineq = length(model.λ)

    dl_dz = zeros(length(model.x))
    for (vi, value) in model.input_cache.dx
        dl_dz[vi.value] = value
    end

    RHS = [dl_dz; zeros(nineq + neq)]

    nv = length(model.x)
    Q = view(LHS, 1:nv, 1:nv)
    partial_grads = if norm(Q) ≈ 0
        -lsqr(LHS, RHS)
    else
        -LHS \ RHS
    end

    dz = partial_grads[1:nv]
    dλ = partial_grads[nv+1:nv+nineq]
    dν = partial_grads[nv+nineq+1:end]

    model.back_grad_cache = QPForwBackCache(dz, dλ, dν)
    return
    # dQ = 0.5 * (dz * z' + z * dz')
    # dq = dz
    # dG = Diagonal(λ) * (dλ * z' + λ * dz') # was: Diagonal(λ) * dλ * z' - λ * dz')
    # dh = -Diagonal(λ) * dλ
    # dA = dν * z'+ ν * dz' # was: dν * z' - ν * dz'
    # db = -dν
    # todo, check MOI signs for dA and dG
end

_linsolve(A, b) = A \ b
# See https://github.com/JuliaLang/julia/issues/32668
_linsolve(A, b::SparseVector) = A \ Vector(b)

# Just a hack that will be removed once we use `MOIU.MatrixOfConstraints`
struct _QPSets end
MOI.Utilities.rows(::_QPSets, ci::MOI.ConstraintIndex) = ci.value

"""
    forward(model::QPDiff)
"""
function forward(model::QPDiff)
    gradient_cache = _gradient_cache(model)
    LHS = gradient_cache.lhs

    z = model.x
    nv = length(z)

    objective_function = _convert(MOI.ScalarQuadraticFunction{Float64}, model.input_cache.objective)
    sparse_array_obj = sparse_array_representation(objective_function, nv)
    dQ = sparse_array_obj.quadratic_terms
    dq = sparse_array_obj.affine_terms

    # The user sets the constraint function in the sense `func`-in-`set` while
    # `db` and `dh` corresponds to the tangents of the set constants. Therefore,
    # we should multiply the constant by `-1`. For `GreaterThan`, we needed to
    # multiply by `-1` to transform it to `LessThan` so it cancels out.
    db = zero(model.ν)
    _fill(isequal(MOI.EqualTo{Float64}), (::Type{MOI.EqualTo{Float64}}) -> true, gradient_cache, model.input_cache, _QPSets(), db)
    dh = zero(model.λ)
    _fill(!isequal(MOI.EqualTo{Float64}), !isequal(MOI.GreaterThan{Float64}), gradient_cache, model.input_cache, _QPSets(), dh)

    dAi = zeros(Int, 0)
    dAj = zeros(Int, 0)
    dAv = zeros(Float64, 0)
    _fill(isequal(MOI.EqualTo{Float64}), isequal(MOI.GreaterThan{Float64}), gradient_cache, model.input_cache, _QPSets(), dAi, dAj, dAv)
    dA = sparse(dAi, dAj, dAv, length(model.ν), nv)

    dGi = zeros(Int, 0)
    dGj = zeros(Int, 0)
    dGv = zeros(Float64, 0)
    _fill(!isequal(MOI.EqualTo{Float64}), isequal(MOI.GreaterThan{Float64}), gradient_cache, model.input_cache, _QPSets(), dGi, dGj, dGv)
    dG = sparse(dGi, dGj, dGv, length(model.λ), nv)

    λ = model.λ
    ν = model.ν
    RHS = [
        dQ * z + dq + dG' * λ + dA' * ν
        λ .* (dG * z) - λ .* dh
        dA * z - db
    ]

    Q = view(LHS, 1:nv, 1:nv)
    partial_grads = if norm(Q) ≈ 0
        -lsqr(LHS', RHS)
    else
        -_linsolve(LHS', RHS)
    end

    dz = partial_grads[1:nv]
    dλ = partial_grads[nv+1:nv+length(λ)]
    dν = partial_grads[nv+length(λ)+1:end]

    model.forw_grad_cache = QPForwBackCache(dz, dλ, dν)
    return
end

function MOI.get(model::QPDiff, ::BackwardOutObjective)
    ∇z = model.back_grad_cache.dz
    z = model.x
    # `∇z * z' + z * ∇z'` doesn't work, see
    # https://github.com/JuliaArrays/LazyArrays.jl/issues/178
    dQ = LazyArrays.@~ (∇z .* z' + z .* ∇z') / 2.0
    return MatrixScalarQuadraticFunction(
        VectorScalarAffineFunction(∇z, 0.0),
        dQ,
    )
end

# quadratic matrix indexes are split by type either == or <=
function _get_dA(model::QPDiff, ci::EQ)
    i = ci.value
    dz = model.back_grad_cache.dz
    dν = model.back_grad_cache.dν
    return lazy_combination(+, dν[i], model.x, model.ν[i], dz)
end
function _get_dA(model::QPDiff, ci::LE)
    i = ci.value
    dz = model.back_grad_cache.dz
    dλ = model.back_grad_cache.dλ
    l = model.λ[i]
    return lazy_combination(+, l * dλ[i], model.x, l, dz)
end


text:
using ZigZagBoomerang
using StaticArrays
using LinearAlgebra
using SparseArrays
using Random
using Test
using Statistics

include("blockchol.jl")

function tofull(Γ::AbstractMatrix)
    n = size(Γ, 1)
    d = size(Γ[1,1],1)
    Γfull = Matrix(Γ)
    [Γfull[i[2],j[2]][i[1],j[1]] for i in vec(CartesianIndices((d,n))), j in vec(CartesianIndices((d,n)))]
    #[Γfull[i[1],j[1]][i[2],j[2]] for i in vec(CartesianIndices((n,d))), j in vec(CartesianIndices((n,d)))]
end
function tofull(x::Vector)
    n = size(x, 1)
    d = size(x[1],1)
    [x[i[2]][i[1]] for i in vec(CartesianIndices((d,n)))]
end
@testset "Vector of SVector" begin
    Random.seed!(2)

    d = 2
    𝕏 = SArray{Tuple{d},Float64,1,d}

    n = 5

    I_nd = Diagonal(fill(SMatrix{d,d}(1.0I), n))

    Γ0 = sparse(Tridiagonal([SMatrix{d,d}([0.0 -0.4;-0.0 -0.0]) for i in 1:n-1], fill(SMatrix{d,d}([1.0 -0.4;-0.4 1.0]), n), [SMatrix{d,d}([0.0 -0.0;-0.4 -0.0]) for i in 1:n-1]))
    Γ = Γ0*Γ0

#    ∇ϕ(x, i, Γ) = sum(Γ[i,j]*x[j] for j in 1:n)
    ∇ϕ(x, i, Γ) = ZigZagBoomerang.idot(Γ, i, x)
    t0 = 0.0
    x0 = randn(𝕏, n)
    θ0 = [randn(𝕏) for i in 1:n]

    μ = 0*x0
    c = [50.0 for i in 1:n]
    σ = [SMatrix{d,d}(1.0I) for i in 1:n]
    Z = ZigZag(Γ, μ, σ; λref=0.02, ρ=0.0) # need refreshments!
    T = 1000.0

    @time trace, (tT, xT, θT), (acc, num) = spdmp(∇ϕ, t0, x0, θ0, T, c, Z, Γ)

    xs = last.(collect(discretize(trace, 2.0)))
    L = lchol(Matrix(Γ))
    Σ = cholinverse!(L, Matrix(I_nd))
    @test_broken mean(norm.(cov(xs) - Σ)) < 3/sqrt(T)
end

@testset "SVector" begin
    Random.seed!(1)

    d = 5
    Γ = sparse(SymTridiagonal(1.0ones(d), -0.4ones(d-1)))
    ∇ϕ!(y, x::T,  Γ) where {T} = T(Γ*x)::SVector

    t0 = 0.0
    x0 = @SVector randn(5)
    θ0 = @SVector ones(Float64, 5)

    μ = 0*x0
    c = 100.0
    σ = [SMatrix{d,d}(1.0I) for i in 1:n]
    BP = BouncyParticle(Γ, x0*0, 0.5, 0., nothing, MMatrix{d,d}(sparse(cholesky(Symmetric(Γ)).L)))
    T = 800.0

    @time trace, (tT, xT, θT), (acc, num) = pdmp(∇ϕ!, t0, x0, θ0, T, c, BP, Γ)
    xs = last.(collect(discretize(trace, 0.01)))
    @test mean(abs.(cov(xs) - inv(Matrix(Γ)))) < 3/sqrt(T)
end


text:
# Generate center coordinates for specific index
"""
    st_coords(ga; mid)

`st_coords` is for regular GeoArray.

# Return
- `X,Y`: A matrix with the dimension of [NX, NY]
"""
function st_coords(ga::GeoArray; mid::Vector{Int} = [1, 1])
    x, y = st_dim(ga; mid)
    meshgrid(x, y)
end

"""
    st_coordinates(ga)

`st_coordinates` is similar as `st_dim`, but for irregular GeoArray.

# Return
Matrix{StaticArrays.SVector{2, Float64}}, e.g.:
```julia
[110.004, 42.9958]  [110.004, 42.9875]
[110.013, 42.9958]  [110.013, 42.9875]
```
"""
function st_coordinates(ga::GeoArray, p::SVector{2, Int})
    ga.f(p .- 0.5)
end
st_coordinates(ga::GeoArray, p::Vector{Int}) = st_coordinates(ga, SVector{2}(p))

function st_coordinates(ga::GeoArray)
    # (ui, uj) = size(ga)[1:2]
    # ci = [st_coordinates(ga, SVector{2}(i,j)) for i in 1:ui, j in 1:uj]
    nrow, ncol = size(ga)[1:2]
    X = zeros(nrow, ncol)
    Y = zeros(nrow, ncol)
    @inbounds for i = 1:nrow, j = 1:ncol
        ci = st_coordinates(ga, SVector{2}(i,j))
        X[i, j] = ci[1]
        Y[i, j] = ci[2]
    end
    X, Y
end

function st_coordinates(ga::GeoArray, dim::Symbol)
    if is_rotated(ga)
        error("This method cannot be used for a rotated GeoArray")
    end
    if dim==:x
        ui = size(ga,1)
        ci = [st_coordinates(ga, SVector{2}(i,1))[1] for i in 1:ui]
    elseif dim==:y
        uj = size(ga,2)
        ci = [st_coordinates(ga, SVector{2}(1,j))[2] for j in 1:uj]
    else
        error("Use :x or :y as second argument")
    end
    return ci
end

"""
    st_coordinates!(ga, x::AbstractUnitRange, y::AbstractUnitRange)

Set AffineMap of `GeoArray` by specifying the *center coordinates* for each x, y
dimension by a `UnitRange`.
"""
function st_coordinates!(ga, x::AbstractUnitRange, y::AbstractUnitRange)
    size(ga)[1:2] != (length(x), length(y)) && error("Size of `GeoArray` $(size(ga)) does not match size of (x,y): $((length(x),length(y))). Note that this function takes *center coordinates*.")
    ga.f = unitrange_to_affine(x, y)
    ga
end

## for missing  ----------------------------------------------------------------
function st_coordinatesnotmissing(ga::GeoArray)
    (ui, uj) = size(ga)[1:2]
    ci = [st_coordinates(ga, SVector{2}(i,j)) for i in 1:ui, j in 1:uj if ~ismissing(ga.A[i, j])]
end

function st_coordinatesmissing(ga::GeoArray)
    (ui, uj) = size(ga)[1:2]
    ci = [st_coordinates(ga, SVector{2}(i,j)) for i in 1:ui, j in 1:uj if ismissing(ga.A[i, j])]
end

function indexmissing(ga::GeoArray)
    (ui, uj) = size(ga)[1:2]
    ci = [[i,j] for i in 1:ui, j in 1:uj if ismissing(ga.A[i, j])]
end

export st_coordinates, st_coordinates!


text:
arg_dict = Dict{Any, Any}(
"physics" => "Euler",
"run_type" => 5,
"jac_type" => 3,
"jac_method" => 2,
"newton_globalize_euler" => true,
"order" => 1,
"real_time" => false,
"variable_type" => :conservative,
"smb_name" => "SRCMESHES/2D_airfoil_np2.smb",
"dmg_name" => "SRCMESHES/2D_airfoil.dmg",
"IC_name" => "ICFreeStream",
"numBC" => 2,
"BC1" => [8],
"BC1_name" => "FreeStreamBC",
"BC2" => [5],
"BC2_name" => "noPenetrationBC",
"aoa" => 2.0, # -2.0,
"Ma" => 0.5,# default working = 0.3
"use_DG" => true,
"Flux_name" => "RoeFlux",
"operator_type" => "SBPOmega",
"exact_visualization" => true,
"t_max" => 500000.0,
"delta_t" => 0.0001,
"res_tol" => 1e-10,
"res_abstol" => 1e-10,
"res_reltol" => 1e-9,
"step_tol" => 1e-10,
"krylov_abstol" => 1e-16,
"krylov_reltol" => 1e-6,
"itermax" => 5000000000, # For RK4

# Integral contributions
"addVolumeIntegrals" => true,
"addFaceIntegrals" => true,
"addStabilization" => true,
"addBoundaryIntegrals" => true,

"writeq" => false,
"write_edge_vertnums" => false,
"write_face_vertnums" => false,
"write_qic" => false,
"writeboundary" => false,
"write_res" => false,
"print_cond" => false,
"write_counts" => false,
"write_vis" => false,
"output_freq" => 100000,
"solve" => true,
"do_postproc" => true,
"calc_functional" => false,
"force_solution_complex" => true,

"functional_name1" => "lift",
"functional_bcs1" => [2],
"num_functionals" => 1,
"need_adjoint" => true, # Necessary for reverse mode

"exact_soln_func" => "ICFreeStream",
)


text:
Random.seed!(1)

nums = [1,2,2,3,3,3,4,4,4,4,4]
lett = ["a","b","b","c","c","c","d","d","d","d","d"]

@testset "oversample" begin
    @test typeof(oversample) <: Function
    @test oversample === MLDataPattern.oversample

    @testset "questionable results to unusual parameters" begin
        @test_throws MethodError oversample(x->x+1, 1)
        @test_throws MethodError oversample(uppercase, "test")
        @test_throws MethodError oversample((1,1:3))
        @test_throws MethodError oversample((1:3,1))
        @test_throws MethodError oversample(identity, 3.)
    end

    @testset "nobs mismatch" begin
        @test_throws DimensionMismatch oversample((1:3,1:4))
        @test_throws DimensionMismatch oversample((1:4,1:3))
        @test_throws DimensionMismatch oversample((X,Yt))
        @test_throws DimensionMismatch oversample((y,Yt))
        @test_throws DimensionMismatch oversample((X,y), true, (ObsDim.First(),))
        @test_throws DimensionMismatch oversample(identity, (X,y), true, (ObsDim.First(),))
    end

    for vals in (nums, lett)
        @testset "Vector of $(eltype(vals)); shuffle=true" begin
            for res in (
                    @inferred(oversample(vals)),
                    @inferred(oversample(vals, 1)),
                    @inferred(oversample(vals, 1, ObsDim.First())),
                    @inferred(oversample(vals, 1, true)),
                    @inferred(oversample(vals, 1, true, ObsDim.First())),
                    @inferred(oversample(vals, true)),
                    @inferred(oversample(vals, true, ObsDim.First())),
                    @inferred(oversample(identity, vals)),
                    @inferred(oversample(identity, vals, true)),
                    @inferred(oversample(identity, vals, 1)),
                    @inferred(oversample(identity, vals, 1, true)),
                    @inferred(oversample(identity, vals, true, ObsDim.First())),
                    @inferred(oversample(identity, vals, 1, true, ObsDim.First())),
                    oversample(vals, fraction=1, shuffle=true, obsdim=1),
                    oversample(identity, vals, fraction=1, shuffle=true, obsdim=1)
                )
                @test typeof(res) <: SubArray{eltype(vals),1}
                @test length(res) == 20
                lm = labelfreq(res)
                @test nlabel(lm) == 4
                @test all(lbl ∈ label(vals) for lbl in label(res))
                @test all(val == 5 for val in values(lm))
                # little hack to test that order is random.
                @test any(res != oversample(vals) for i = 1:5)
            end
        end

        @testset "Vector of $(eltype(vals)); shuffle=false" begin
            for res in (
                    @inferred(oversample(vals, false)),
                    @inferred(oversample(vals, 1, false)),
                    @inferred(oversample(vals, false, ObsDim.First())),
                    @inferred(oversample(vals, 1, false, ObsDim.First())),
                    @inferred(oversample(identity, vals, false)),
                    @inferred(oversample(identity, vals, 1, false)),
                    @inferred(oversample(identity, vals, false, ObsDim.First())),
                    @inferred(oversample(identity, vals, 1, false, ObsDim.First())),
                    oversample(vals, shuffle=false, obsdim=1),
                    oversample(vals, fraction=1, shuffle=false, obsdim=1),
                    oversample(identity, vals, fraction=1, shuffle=false, obsdim=1)
                )
                @test typeof(res) <: SubArray{eltype(vals),1}
                @test length(res) == 20
                lm = labelfreq(res)
                @test nlabel(lm) == 4
                @test all(lbl ∈ label(vals) for lbl in label(res))
                @test all(val == 5 for val in values(lm))
                @test res[1:length(vals)] == vals
                # little hack to test that order is deterministic.
                @test all(res == oversample(vals,false) for i = 1:3)
            end
        end
    end

    @testset "Basic" begin
        #n_src = 200
        src = [1, 2,2,2, 3,3,3,3,3, 4,4,4,4,4,4,4,4]
        oversampled = oversample(src)
        @test all(counts(oversampled).==counts(oversampled)[1])
        @test all(x ∈ oversampled for x in unique(src))
        @test length(oversampled) > length(src)
        oversampled = oversample(src, 0.5)
        @test countmap(oversampled, alg = :dict) == Dict(4=>8,2=>4,3=>5,1=>4)
        for oversampled in (
                @inferred(oversample(src, 0.2)),
                @inferred(oversample(src, 0.2, true)),
                @inferred(oversample(src, 0.2, false)),
                oversample(src, fraction=0.2),
                @inferred(oversample(identity, src, 0.2)),
                @inferred(oversample(identity, src, 0.2, false)),
                oversample(identity, src, fraction=0.2),
            )
            @test all(counts(oversampled) .>= 2)
            @test any(counts(oversampled) .<= 3)
        end
    end

    @testset "ObsView" begin
        n_src = 200
        src = rand([1, 2,2, 3,3,3, 4,4,4,4], n_src)
        src2 = obsview(src)
        @test typeof(src2) <: ObsView
        oversampled2 = getobs.(@inferred(oversample(src2)))
        @test all(counts(oversampled2).==counts(oversampled2)[1])
        @test all(x ∈ oversampled2 for x in unique(src))
        @test length(oversampled2) > n_src
    end

    @testset "Advanced" begin
        n_src = 200
        data = rand(n_src, 50) #50 features
        lbs = rand([1, 2,2, 3,3,3, 4,4,4,4], n_src)

        od = ObsDim.First(), ObsDim.First()

        data_os, lbls_os = oversample((data, lbs); obsdim=od)
        @test all(counts(lbls_os).==counts(lbls_os)[1])
        @test all(x ∈ lbls_os for x in unique(lbls_os))
        @test nobs(data_os, ObsDim.First()) == nobs(lbls_os)
        @test nobs(lbls_os) > n_src
    end

    @testset "MultiFactor Label" begin
        n_factors = 4
        n_observations = 200

        src = rand([1, 2,2, 3,3,3, 4,4,4,4], (n_factors, n_observations))
        oversampled = oversample(src)

        src_cnts = labelfreq(obsview(src))
        os_cnts = labelfreq(obsview(oversampled))

        @test Set(keys(os_cnts))==Set(keys(src_cnts))
        @test size(oversampled,2) > n_observations
        @test all(cnt == first(os_cnts)[2] for (kk, cnt) in os_cnts)
    end

    @testset "MultiFactor Label with fun" begin
        n_observations = 200
        src = rand([1, 2,2, 3,3,3, 4,4,4,4], (2, n_observations))
        sampled = oversample(x->x[1]>x[2], src)
        @assert sum(src[1,:].>src[2,:])!=n_observations//2

        @test size(sampled,2) > n_observations
        @test sum(sampled[1,:].>sampled[2,:]) == sum(sampled[1,:].<=sampled[2,:])
    end
end

@testset "undersample" begin
    @test typeof(undersample) <: Function
    @test undersample === MLDataPattern.undersample

    @testset "questionable results to unusual parameters" begin
        @test_throws MethodError undersample(x->x+1, 1)
        @test_throws MethodError undersample(uppercase, "test")
        @test_throws MethodError undersample((1,1:3))
        @test_throws MethodError undersample((1:3,1))
        @test_throws MethodError undersample(identity, 3.)
    end

    @testset "nobs mismatch" begin
        @test_throws DimensionMismatch undersample((1:3,1:4))
        @test_throws DimensionMismatch undersample((1:4,1:3))
        @test_throws DimensionMismatch undersample((X,Yt))
        @test_throws DimensionMismatch undersample((y,Yt))
        @test_throws DimensionMismatch undersample((X,y), true, (ObsDim.First(),))
        @test_throws DimensionMismatch undersample(identity, (X,y), true, (ObsDim.First(),))
    end

    for vals in (nums, lett)
        @testset "Vector of $(eltype(vals)); shuffle=true" begin
            for res in (
                    @inferred(undersample(vals, true)),
                    @inferred(undersample(vals, true, ObsDim.First())),
                    @inferred(undersample(identity, vals, true)),
                    @inferred(undersample(identity, vals, true, ObsDim.First())),
                    undersample(vals, shuffle=true, obsdim=1),
                    undersample(identity, vals, shuffle=true, obsdim=1)
                )
                @test typeof(res) <: SubArray{eltype(vals),1}
                @test length(res) == 4
                lm = labelfreq(res)
                @test nlabel(lm) == 4
                @test all(lbl ∈ label(vals) for lbl in label(res))
                @test all(val == 1 for val in values(lm))
                # little hack to test that order is random.
                @test any(res != undersample(vals,true) for i = 1:5)
            end
        end

        @testset "Vector of $(eltype(vals)); shuffle=false" begin
            for res in (
                    @inferred(undersample(vals)),
                    @inferred(undersample(vals, false)),
                    @inferred(undersample(vals, false, ObsDim.First())),
                    @inferred(undersample(identity, vals)),
                    @inferred(undersample(identity, vals, false)),
                    @inferred(undersample(identity, vals, false, ObsDim.First())),
                    undersample(vals, shuffle=false, obsdim=1),
                    undersample(identity, vals, shuffle=false, obsdim=1)
                )
                @test typeof(res) <: SubArray{eltype(vals),1}
                @test length(res) == 4
                lm = labelfreq(res)
                @test nlabel(lm) == 4
                @test all(lbl ∈ label(vals) for lbl in label(res))
                @test all(val == 1 for val in values(lm))
                @test res == sort(unique(vals))
                # little hack to test that order is deterministic.
                @test all(res == undersample(vals) for i = 1:3)
            end
        end
    end

    @testset "Basic" begin
        n_src = 2000
        src = rand([1, 2,2, 3,3,3, 4,4,4,4], n_src)
        sampled = undersample(src)
        @test all(counts(sampled).==counts(sampled)[1])
        @test all( x ∈ sampled for x in unique(src))
        @test length(sampled) < n_src
    end

    @testset "ObsView" begin
        n_src = 200
        src = rand([1, 2,2, 3,3,3, 4,4,4,4], n_src)
        src2 = obsview(src)
        @test typeof(src2) <: ObsView
        sampled = getobs.(@inferred(undersample(src2)))
        @test all(counts(sampled).==counts(sampled)[1])
        @test all( x ∈ sampled for x in unique(src))
        @test length(sampled) < n_src
    end

    @testset "Advanced" begin
        n_src = 200
        data = rand(n_src, 50) #50 features
        lbs = rand([1, 2,2, 3,3,3, 4,4,4,4], n_src)

        od = ObsDim.First(), ObsDim.First()

        data_os, lbls_os = undersample((data, lbs); obsdim=od)
        @test all(counts(lbls_os).==counts(lbls_os)[1])
        @test all( x ∈ lbls_os for x in unique(lbls_os))
        @test nobs(data_os, ObsDim.First()) == nobs(lbls_os)
        @test nobs(lbls_os) < n_src
    end

    @testset "MultiFactor Label" begin
        n_factors = 4
        n_observations = 2000

        src = rand([1, 2,2, 3,3,3, 4,4,4,4], (n_factors, n_observations))
        sampled = undersample(src)

        src_cnts = labelfreq(obsview(src))
        os_cnts =  labelfreq(obsview(sampled))

        @test Set(keys(os_cnts))==Set(keys(src_cnts))
        @test size(sampled,2) < n_observations

        first_os_count = first(os_cnts)[2]
        @test all(cnt == first_os_count for (kk, cnt) in os_cnts)
    end

    @testset "MultiFactor Label with fun" begin
        n_observations = 200

        src = rand([1, 2,2, 3,3,3, 4,4,4,4], (2, n_observations))
        sampled = undersample(x->x[1]>x[2], src)
        @assert sum(src[1,:].>src[2,:])!=n_observations//2

        @test size(sampled,2) < n_observations
        @test sum(sampled[1,:].>sampled[2,:]) == sum(sampled[1,:].<=sampled[2,:])
    end
end


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

@testset "ilu0" begin
    @testset for elty in [Float32,Float64,Complex64,Complex128]
        A = rand(elty,m,m)
        @testset "csr" begin
            d_A = CudaSparseMatrixCSR(sparse(A))
            info = CUSPARSE.sv_analysis('N','G','U',d_A,'O')
            d_B = CUSPARSE.ilu0('N',d_A,info,'O')
            h_B = to_host(d_B)
            Alu = lufact(full(A),Val{false})
            Ac = sparse(Alu[:L]*Alu[:U])
            h_A = ctranspose(h_B) * h_B
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
        @testset "csc" begin
            d_A = CudaSparseMatrixCSC(sparse(A))
            info = CUSPARSE.sv_analysis('N','G','U',d_A,'O')
            d_B = CUSPARSE.ilu0('N',d_A,info,'O')
            h_B = to_host(d_B)
            Alu = lufact(full(A),Val{false})
            Ac = sparse(Alu[:L]*Alu[:U])
            h_A = ctranspose(h_B) * h_B
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
    end
end

@testset "ilu02" begin
    @testset for elty in [Float32,Float64,Complex64,Complex128]
        @testset "csr" begin
            A = rand(elty,m,m)
            A += transpose(A)
            A += m * eye(elty,m)
            d_A = CudaSparseMatrixCSR(sparse(A))
            d_B = CUSPARSE.ilu02(d_A,'O')
            h_A = to_host(d_B)
            Alu = lufact(full(A),Val{false})
            Ac = sparse(Alu[:L]*Alu[:U])
            h_A = ctranspose(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
        @testset "csc" begin
            A = rand(elty,m,m)
            A += transpose(A)
            A += m * eye(elty,m)
            d_A = CudaSparseMatrixCSC(sparse(A))
            d_B = CUSPARSE.ilu02(d_A,'O')
            h_A = to_host(d_B)
            Alu = lufact(full(A),Val{false})
            Ac = sparse(Alu[:L]*Alu[:U])
            h_A = ctranspose(h_A) * h_A
            @test h_A.rowval ≈ Ac.rowval
            @test reduce(&, isfinite.(h_A.nzval))
        end
    end
end


text:
module Brutus

function __init__()
    ccall((:brutus_init, "libbrutus"), Cvoid, (Any,), @__MODULE__)
end

@enum DumpOption::UInt8 begin
    DumpIRCode        = 0
    DumpTranslated    = 1
    DumpCanonicalized = 2
    DumpLoweredToStd  = 4
    DumpLoweredToLLVM = 8
    DumpTranslateToLLVM = 16
end

function find_invokes(IR)
    callees = Core.MethodInstance[]
    for stmt in IR.stmts
        if stmt isa Expr
            if stmt.head == :invoke
                mi = stmt.args[1]
                push!(callees, mi)
            end
        end
    end
    return callees
end

# Emit MLIR IR to stdout
function emit(@nospecialize(ft), @nospecialize(tt);
              emit_fptr::Bool=true,
              dump_options::Vector{DumpOption}=DumpOption[])
    name = (ft <: Function) ? nameof(ft.instance) : nameof(ft)

    # get first method instance matching signature
    entry_mi = get_methodinstance(Tuple{ft, tt...})
    IR, rt = code_ircode(entry_mi)

    if DumpIRCode in dump_options
        println("return type: ", rt)
        println("IRCode:\n")
        println(IR)
    end

    worklist = [IR]
    methods = Dict{Core.MethodInstance, Tuple{Core.Compiler.IRCode, Any}}(
        entry_mi => (IR, rt)
    )

    while !isempty(worklist)
        code = pop!(worklist)
        callees = find_invokes(code)
        for callee in callees
            if !haskey(methods, callee)
                _code, _rt = code_ircode(callee)

                methods[callee] = (_code, _rt)
                push!(worklist, _code)
            end
        end
    end

    # generate LLVM bitcode and load it
    dump_flags = reduce(|, map(UInt8, dump_options), init=0)
    fptr = ccall((:brutus_codegen, "libbrutus"),
                 Ptr{Nothing},
                 (Any, Any, Cuchar, Cuchar),
                 methods, entry_mi, emit_fptr, dump_flags)
    if convert(Int, fptr) == 0
        return nothing
    end

    return fptr, rt
end

###
# Call Interface
###

struct MemrefDescriptor{T, N}
    allocated :: Ptr{T}
    aligned   :: Ptr{T}
    offset    :: Int64
    size      :: NTuple{N, Int64}
    strides   :: NTuple{N, Int64}
end

function Base.convert(::Type{MemrefDescriptor{T, N}}, arr::Array{T, N}) where {T, N}
    allocated = pointer(arr)
    aligned   = pointer(arr)
    offset    = Int64(0)
    size      = Base.size(arr)
    strides   = (1, cumprod(size)[1:end-1]...)

    # colum-major to row-major
    size = reverse(size)
    strides = reverse(strides)

    MemrefDescriptor{T, N}(allocated, aligned, offset, size, strides)
end

struct Thunk{F, RT, TT}
    f::F
    ptr::Ptr{Cvoid}
end

function thunk(f::F, tt::TT=Tuple{}) where {F<:Base.Callable, TT<:Type}
    # TODO: Use GPUCompiler.cached_compilation
    result = emit(F, (tt.parameters...,))
    if result === nothing
        error("failed to emit function")
    end
    fptr, rt, = result
    Thunk{F, rt, tt}(f, fptr)
end

# Need to pass struct as pointer, to match cifacme ABI
abi(::Type{<:Array{T, N}}) where {T, N} = Ref{MemrefDescriptor{T, N}}
function abi(T::DataType)
    if isprimitivetype(T)
        return T
    else
        return Any
    end
end

@inline (thunk::Thunk)(args...) = __call(thunk, args)

@generated function __call(thunk::Thunk{f, RT, TT}, args::TT) where {f, RT, TT}
    nargs = map(abi, (args.parameters...,))
    _args = (:(args[$i]) for i in 1:length(nargs))

    # Insert function type up-front
    nargs = (Any, nargs...)

    expr = quote
        ccall(thunk.ptr, $(abi(RT)), ($(nargs...),), thunk.f, $(_args...))::RT
    end
    return expr
end

function call(f::F, args...) where F
    TT = Tuple{map(Core.Typeof, args)...}
    return thunk(f, TT)(args...)
end

function get_methodinstance(@nospecialize(sig);
                            world=Base.get_world_counter(),
                            interp=Core.Compiler.NativeInterpreter(world))
    ms = Base._methods_by_ftype(sig, 1, Base.get_world_counter())
    @assert length(ms) == 1
    m = ms[1]
    mi = ccall(:jl_specializations_get_linfo,
               Ref{Core.MethodInstance}, (Any, Any, Any),
               m[3], m[1], m[2])
    return mi
end

function code_ircode_by_signature(@nospecialize(sig);
                                  world=Base.get_world_counter(),
                                  interp=Core.Compiler.NativeInterpreter(world))
    return [code_ircode(ccall(:jl_specializations_get_linfo,
                              Ref{Core.MethodInstance},
                              (Any, Any, Any),
                              data[3], data[1], data[2]);
                        world=world, interp=interp)
            for data in Base._methods_by_ftype(sig, -1, world)]
end

function code_ircode(@nospecialize(f), @nospecialize(types=Tuple);
                     world=Base.get_world_counter(),
                     interp=Core.Compiler.NativeInterpreter(world))
    return [code_ircode(mi; world=world, interp=interp)
            for mi in Base.method_instances(f, types, world)]
end

function code_ircode(mi::Core.Compiler.MethodInstance;
                     world=Base.get_world_counter(),
                     interp=Core.Compiler.NativeInterpreter(world))
    ccall(:jl_typeinf_begin, Cvoid, ())
    result = Core.Compiler.InferenceResult(mi)
    frame = Core.Compiler.InferenceState(result, false, interp)
    frame === nothing && return nothing
    if Core.Compiler.typeinf(interp, frame)
        opt_params = Core.Compiler.OptimizationParams(interp)
        opt = Core.Compiler.OptimizationState(frame, opt_params, interp)
        ir = Core.Compiler.run_passes(opt.src, opt.nargs - 1, opt)
        opt.src.inferred = true
    end
    ccall(:jl_typeinf_end, Cvoid, ())
    frame.inferred || return nothing
    # TODO(yhls): Fix this upstream
    resize!(ir.argtypes, opt.nargs)
    return ir => Core.Compiler.widenconst(result.result)
end

include("reflection.jl")

end # module

max_stars_repo_name: femtomc/brutus

max_issues_repo_name: femtomc/brutus

max_forks_repo_name: femtomc/brutus


text:
#
# Copyright (c) 2021 Chernobylll

# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMI
# using Plots TODO: can be deleted

pathToFMU = joinpath(dirname(@__FILE__), "../model/Dymola/2020x/IO.fmu")

myFMU = fmiLoad(pathToFMU)
fmiInstantiate!(myFMU; loggingOn=true)

fmiSetupExperiment(myFMU, 0.0)

fmiEnterInitializationMode(myFMU)

fmiGetString(myFMU, "p_string")

rndReal = 100 * rand()
rndInteger = round(Integer, 100 * rand())
rndBoolean = rand() > 0.5
rndString = "Not random!"

fmiSetReal(myFMU, "p_real", rndReal)
display("$rndReal == $(fmiGetReal(myFMU, "p_real"))")

fmiSetInteger(myFMU, "p_integer", rndInteger)
display("$rndInteger == $(fmiGetInteger(myFMU, "p_integer"))")

fmiSetBoolean(myFMU, "p_boolean", rndBoolean)
display("$rndBoolean == $(fmiGetBoolean(myFMU, "p_boolean"))")

fmiSetString(myFMU, "p_string", rndString)
display("$rndString == $(fmiGetString(myFMU, "p_string"))")

fmiExitInitializationMode(myFMU)

fmiUnload(myFMU)

max_stars_repo_name: Chernobylll/fmi-codes-and-variants

max_issues_repo_name: Chernobylll/fmi-codes-and-variants

max_forks_repo_name: Chernobylll/fmi-codes-and-variants


text:
import Base: merge

###### merge ####################

#function merge{T,N}(ta1::TimeArray{T,N}, ta2::TimeArray{T,N}; method="inner")
function merge{T}(ta1::TimeArray{T}, ta2::TimeArray{T}; colnames = [""], method="inner",)
    tstamp = [Date(1,1,1):Year(1):Date(length(ta1),1,1)]
    n      = 1
      for i in 1:length(ta1)
        for j in 1:length(ta2)
          if ta1.timestamp[i] == ta2.timestamp[j]
            tstamp[n] = ta1.timestamp[i]
            n+=1
          end
        end
      end
    tstamp = tstamp[1:n-1] # trim down the length, if necessary

    val1 = ta1[tstamp].values
    val2 = ta2[tstamp].values
    vals = hcat(val1, val2)

    #if length(colnames) != 2
    if length(colnames) < 2
    cnames = copy(ta1.colnames) # otherwise ta1 gets contaminated
      for m in 1:length(ta2.colnames)
        push!(cnames, ta2.colnames[m])
      end
     else cnames = colnames
     end

    TimeArray(tstamp, vals, cnames)
end

# collapse ######################

function collapse{T,N}(ta::TimeArray{T,N}, f::Function; period::Function=week)

  w = [period(ta.timestamp[t]) for t in 1:length(ta)] # get weekly id from entire array
  z = Int[]; j = 1
  for i=1:length(ta) - 1 # create unique period ID array
    if w[i] != w[i+1]
      push!(z, j)
      j = j+1
    else
      push!(z,j)
    end
  end

  # account for last row
  w[length(ta)]  ==  w[length(ta)-1] ? # is the last row the same period as 2nd to last row?
  push!(z, z[length(z)]) :
  push!(z, z[length(z)] + 1)

  # pre-allocate timestamp and value arrays
  tstamps = [Date(1,1,1):Year(1):Date(maximum(z),1,1)]
  vals    = zeros(maximum(z)) # number of unique periods
  #replace their values except for the last row
  for i = 1:maximum(z)-1  # iterate over period ID groupings
    temp       = ta[findfirst(z .== i):findfirst(z .== i+1)-1] # period's worth that will be squished
    tstamps[i] = temp[length(temp)].timestamp[1]
    vals[i]    = f(temp.values)
  end
  # and once again account for the last temp that isn't looped on above
  lasttemp = ta[findfirst(z .== maximum(z)):length(ta)]
  lasttempindex = lasttemp[length(lasttemp)].timestamp
  lasttempvalue = f(ta.values)

  # complete the tstamp and vals arrays with the last value
  tstamps[length(tstamps)] = lasttempindex[1]
  vals[length(vals)]       = lasttempvalue

  TimeArray(tstamps, vals, ta.colnames)
end


text:
@testset "VREP constructor" begin
    N = Float64

    # dimension unknown => VPolytope
    @test setrep(VREP(δ=0.1)) == VPolytope{N, Vector{N}}

    # dimension is two => VPolygon with static arrays
    @test setrep(VREP(δ=0.1, dim=2)) == VPolygon{N, SVector{2, N}}

    # dimension is two => VPolygon with static arrays
    @test setrep(VREP(δ=0.1, dim=2, static=true)) == VPolygon{N, SVector{2, N}}

    # dimension is two, no static
    @test setrep(VREP(δ=0.1, dim=2, static=false)) == VPolygon{N, Vector{N}}

    # dimension is two, use static arrays by default
    @test setrep(VREP(δ=0.1, dim=5)) == VPolytope{N, SVector{5, N}}
end

@testset "VREP static arrays conversions" begin
    N = Float64

    V2 = convert(VPolygon, BallInf(zeros(2), 1.0))
    V3 = convert(VPolytope, BallInf(zeros(3), 1.0))

    # no-op
    @test RA._reconvert(V2, Val(false), Missing) == V2
    @test RA._reconvert(V3, Val(false), Missing) == V3

    # convert to static arrays
    @test RA._reconvert(V2, Val(true), Val(2)) isa VPolygon{Float64, SVector{2, N}}
    @test RA._reconvert(V3, Val(true), Val(3)) isa VPolytope{Float64, SVector{3, N}}
end

@testset "VREP for 2D system" begin
    A  = [0   1.
          -1  0]
    X0 = VPolygon([[1.0, 0.0], [1.2, 0.0], [1.1, 0.2]])
    prob = @ivp(x' = A*x, x(0) ∈ X0)
    sol = solve(prob, tspan=(0.0, 2.0), alg=VREP(δ=1e-2))
    @test setrep(sol) == VPolygon{Float64, Vector{Float64}}

    # statically sized array
    sol = solve(prob, tspan=(0.0, 2.0), alg=VREP(δ=1e-2, static=true, dim=2))
    @test setrep(sol) == VPolygon{Float64, SVector{2, Float64}}
end

@testset "VREP for 5D system" begin
    A  = [1. 0  0 -5
          4 -2  4 -3
         -4  0  0  1
          5 -2  2  3.]

    X0_12 = VPolygon([[1.0, 0.0], [1.2, 0.0], [1.1, 0.2]])
    X0_34 = BallInf(zeros(2), 0.2)
    X0 = VPolytope(vertices_list(X0_12 × X0_34))

    # default Polyhedra backend works, but shows warning messages of this sort:
    # glp_simplex: unable to recover undefined or non-optimal solution
    #prob = @ivp(x' = A*x, x(0) ∈ X0)
    #sol = solve(prob, tspan=(0.0, 1.0), alg=VREP(δ=1e-3, static=true, dim=4))

    # change the backend
    prob = @ivp(x' = A*x, x(0) ∈ X0)
    sol = solve(prob, tspan=(0.0, 1.0), alg=VREP(δ=1e-3, static=true, dim=4, backend=CDDLib.Library()))
end

@testset "VREP for sdof with distinct approximations" begin
    prob, _ = harmonic_oscillator()
    tmax = 2.0

    a = solve(prob, tspan=(0.0, tmax), alg=VREP(δ=0.1, approx_model=Forward(setops=:concrete)));
    b = solve(prob, tspan=(0.0, tmax), alg=VREP(δ=0.1, approx_model=StepIntersect(setops=:concrete)));
    c = solve(prob, tspan=(0.0, tmax), alg=VREP(δ=0.1, approx_model=CorrectionHull(exp=:base)));
end


text:
using Test
using TightBindingLattice

@testset "iscompatible" begin
    # hypercube and operation
    @testset "orthocube-operation" begin
        #hypercube = HypercubicLattice([2 0; 0 2])
        orthocube = OrthoCube([2 0; 0 2])
        @test_throws DimensionMismatch iscompatible(orthocube, TranslationOperation([1, 0, 0]))
        @test iscompatible(orthocube, TranslationOperation([1, 0]))

        @test_throws DimensionMismatch iscompatible(orthocube, PointOperation([1 0 0; 0 1 0; 0 0 0]))
        @test iscompatible(orthocube, PointOperation([0 1; 1 0]))
        @test !iscompatible(orthocube, PointOperation([0 2; 1 0]))
        @test !iscompatible(OrthoCube([3 0; 0 2]), PointOperation([0 1; 1 0]))
    end

    # hypercube and symmetry
    @testset "hypercube-symmetry" begin
        orthocube = OrthoCube([2 0; 0 2])
        @test  iscompatible(orthocube, TranslationSymmetry([2 0; 0 2]))
        @test  iscompatible(orthocube, TranslationSymmetry([2 2; 0 2]))
        @test !iscompatible(orthocube, TranslationSymmetry([1 0; 0 1]))
        @test !iscompatible(orthocube, TranslationSymmetry([2 0; 0 1]))

        psym = project(PointSymmetryDatabase.find("4mm"), [1 0 0; 0 1 0])
        @test  iscompatible(orthocube, psym)
        @test  iscompatible(OrthoCube([2 2; 0 2]), psym)
        @test !iscompatible(OrthoCube([3 0; 0 2]), psym)
    end

    # translation symmetry and point operation
    @testset "translation-point" begin
        @test  iscompatible(TranslationSymmetry([2 0; 0 2]), PointOperation([0 1; 1 0]))
        @test  iscompatible(TranslationSymmetry([2 2; 0 2]), PointOperation([0 1; 1 0]))
        @test !iscompatible(TranslationSymmetry([3 0; 0 2]), PointOperation([0 1; 1 0]))

        psym = project(PointSymmetryDatabase.find("4mm"), [1 0 0; 0 1 0])
        @test  iscompatible(TranslationSymmetry([2 0; 0 2]), psym)
        @test  iscompatible(TranslationSymmetry([2 2; 0 2]), psym)
        @test !iscompatible(TranslationSymmetry([3 0; 0 2]), psym)

        @test  iscompatible(psym, TranslationSymmetry([2 0; 0 2]))
        @test  iscompatible(psym, TranslationSymmetry([2 2; 0 2]))
        @test !iscompatible(psym, TranslationSymmetry([3 0; 0 2]))
    end

    @testset "translation-irrep-point-operation" begin
        # shape matching, momentum check
        tsym = TranslationSymmetry([3 0; 0 3])
        for irrep_index in 1:num_irreps(tsym)
            k = tsym.fractional_momenta[irrep_index]
            @test  iscompatible(tsym, irrep_index, PointOperation([0 1; 1  0])) == (k[1] == k[2])
            @test  iscompatible(tsym, irrep_index, PointOperation([1 0; 0 -1])) == (k[2] == 0 || k[2] == 1//2)
        end
        tsym = TranslationSymmetry([3 3; 0 3])
        for irrep_index in 1:num_irreps(tsym)
            k = tsym.fractional_momenta[irrep_index]
            @test  iscompatible(tsym, irrep_index, PointOperation([0 1; 1  0])) == (k[1] == k[2])
            @test  iscompatible(tsym, irrep_index, PointOperation([1 0; 0 -1])) == (k[2] == 0 || k[2] == 1//2)
        end

        # shape mismatch
        tsym = TranslationSymmetry([3 0; 0 2])
        for irrep_index in 1:num_irreps(tsym)
            k = tsym.fractional_momenta[irrep_index]
            @test !iscompatible(tsym, irrep_index, PointOperation([0 1; 1 0])) # not compatible with tsym
            @test  iscompatible(tsym, irrep_index, PointOperation([1 0; 0 -1])) == (k[2] == 0 || k[2] == 1//2)
        end

        psym = project(PointSymmetryDatabase.find("4mm"), [1 0 0; 0 1 0])
        @test  iscompatible(TranslationSymmetry([2 0; 0 2]), psym)
        @test  iscompatible(TranslationSymmetry([2 2; 0 2]), psym)
        @test !iscompatible(TranslationSymmetry([3 0; 0 2]), psym)
    end

    @testset "translation-irrep-point-symmetry" begin
        tsym = TranslationSymmetry([3 0; 0 3])
        psym = project(PointSymmetryDatabase.find("4mm"), [1 0 0; 0 1 0])
        for irrep_index in 1:num_irreps(tsym)
            k = tsym.fractional_momenta[irrep_index]
            @test iscompatible(tsym, irrep_index, psym) == (k in [[0//1, 0//1], [1//2, 1//2]])
            # 4mm only at Gamma and M points
        end

        tsym = TranslationSymmetry([4 0; 0 4])
        psym = project(PointSymmetryDatabase.find("4mm"), [1 0 0; 0 1 0])
        for irrep_index in 1:num_irreps(tsym)
            k = tsym.fractional_momenta[irrep_index]
            @test iscompatible(tsym, irrep_index, psym) == (k in [[0//1, 0//1], [1//2, 1//2]])
            # 4mm only at Gamma and M points
        end

        tsym = TranslationSymmetry([4 4; 0 4])
        psym = project(PointSymmetryDatabase.find("4mm"), [1 0 0; 0 1 0])
        for irrep_index in 1:num_irreps(tsym)
            k = tsym.fractional_momenta[irrep_index]
            @test iscompatible(tsym, irrep_index, psym) == (k in [[0//1, 0//1], [1//2, 1//2]])
            # 4mm only at Gamma and M points
        end

        # shape mismatch
        tsym = TranslationSymmetry([4 2; 0 4])
        psym = project(PointSymmetryDatabase.find("4mm"), [1 0 0; 0 1 0])
        for irrep_index in 1:num_irreps(tsym)
            k = tsym.fractional_momenta[irrep_index]
            @test !iscompatible(tsym, irrep_index, psym)
            # 4mm already not compatible with lattice shape
        end
    end

    @testset "lattice-symmetry" begin
        unitcell = make_unitcell([1.0 0.0; 0.0 1.0]; SiteType=String)
        addsite!(unitcell, "X", FractCoord([0, 0], [0.5, 0.0]))

        lattice = make_lattice(unitcell, [2 0; 0 2])
        @test_throws DimensionMismatch iscompatible(lattice, TranslationOperation([1,2,3]))
        for i in -5:5, j in -5:5
            @test iscompatible(lattice, TranslationOperation([i, j]))
        end
        @test iscompatible(lattice, PointOperation([1 0; 0 -1]))
        @test !iscompatible(lattice, PointOperation([0 1; 1 0]))

        @test iscompatible(lattice, TranslationSymmetry([2 0; 0 2]))
        @test iscompatible(lattice, TranslationSymmetry([2 2; 0 2]))
        # ^- This is a choice. Whether to require equality, or just compatibility

        psym = project(PointSymmetryDatabase.get(7), [1 0 0; 0 1 0])
        @test iscompatible(lattice, psym)

        psym = project(PointSymmetryDatabase.get(13), [1 0 0; 0 1 0])
        @test !iscompatible(lattice, psym)
    end
end


text:
using Flux: onehotbatch
using BSON
using GZip
using Tar
using ProgressMeter
using DelimitedFiles

const TRAINING_OUT_DIR = "train"
const TEXT_DIR = "train_txt"
const PHONES = split("h#	q	eh	dx	iy	r	ey	ix	tcl	sh	ow	z	s	hh	aw	m	t	er	l	w	aa	hv	ae	dcl	y	axr	d	kcl	k	ux	ng	gcl	g	ao	epi	ih	p	ay	v	n	f	jh	ax	en	oy	dh	pcl	ah	bcl	el	zh	uw	pau	b	uh	th	ax-h	em	ch	nx	eng")
phn2num = Dict(phone=>i for (i, phone) in enumerate(PHONES))
phn2num["sil"] = 1

"""
    createData(data_dir, out_dir)

Extracts data from files in `data_dir` and saves results in `out_dir`.
"""
function createData(out_dir)

  println("Downloading data...")
  download("https://era.library.ualberta.ca/items/d0dd87c9-8091-467d-8c30-bd23a2a1fb55/download/192e48bf-ffbd-4cb9-9e42-3cde0312994d", "mald.tar.gz")

  print("Decompressing data... ")

  local s
  GZip.open("mald.tar.gz") do fh
    s = read(fh)
  end

  open("mald.tar", "w") do w
    write(w, s)
  end
  rm("mald.tar.gz")

  println("DONE")
  print("Extracting data... ")
  Tar.extract("mald.tar", "extracted_mald")
  rm("mald.tar")
  println("DONE")

  data_dir = "extracted_mald/sR_bson_subset"

  @showprogress for fname in readdir(data_dir)

    BSON.@load joinpath(data_dir, fname) mfccs labs
    x = mfccs
    labels = [phn2num[x] for x in vec(labs)]
    class_nums = collect(1:62)

    y = Int.(onehotbatch(labels, class_nums))
    BSON.@save joinpath(out_dir, fname) x y
  end
  rm("extracted_mald", recursive=true)
end

"""
    createTextData(from_dir, out_dir)

Reads in BSON data from `from_dir` and writes it out as a tab-delimited file in `out_dir`.
"""
function createTextData(from_dir, out_dir)

  @showprogress for fname in readdir(from_dir)

    fname = joinpath(from_dir, fname)
    BSON.@load fname x y

    fname = basename(fname)
    fname = joinpath(out_dir, fname)
    writedlm(replace(fname, ".bson" => ".txt"), x, '\t')
    writedlm(replace(fname, ".bson" => "_labs.txt"), y, '\t')
  end
end

if ! isdir(TRAINING_OUT_DIR) mkdir(TRAINING_OUT_DIR) end
createData(TRAINING_OUT_DIR)

println("Converting to text")
if ! isdir(TEXT_DIR) mkdir(TEXT_DIR) end
createTextData(TRAINING_OUT_DIR, TEXT_DIR)

max_stars_repo_name: maetshju/flux-ctc-grad

max_issues_repo_name: maetshju/flux-ctc-grad

max_forks_repo_name: maetshju/flux-ctc-grad


text:
#JOB FILEPATHS
#GFF3 feature database, FASTA genome and index paths
danio_gff_path = "/bench/PhD/seq/GRCz11/Danio_rerio.GRCz11.94.gff3"
danio_genome_path = "/bench/PhD/seq/GRCz11/GCA_000002035.4_GRCz11_genomic.fna"
danio_gen_index_path = "/bench/PhD/seq/GRCz11/GCA_000002035.4_GRCz11_genomic.fna.fai"
#sample record and hmm serialisation output path
sample_output = "/bench/PhD/NGS_binaries/BBM/refinement_samples"
#non registered libs

#GENERAL SETUP
@info "Loading libraries..."
using Distributed,Serialization

#JOB CONSTANTS
#CONSTANTS FOR GENOMIC SAMPLING
const sample_set_length = Int64(16e6)
const sample_window_min = 10
const sample_window_max = 3000
const perigenic_pad = 500
const partitions = 3 #exonic, periexonic, intragenic

#Setup sampling workers
@info "Spawning workers..."
pool_size = partitions                          #number of workers to use
worker_pool = addprocs(pool_size) # add processes up to the worker pool size + 1 control process
@everywhere using BioBackgroundModels, Random
@everywhere Random.seed!(1)

#GET SEQUENCE OBSERVATIONS TO TRAIN AND TEST MODELS - PARTITIONING GENOME INTO EXONIC, PERIEXONIC, INTRAGENIC SEQUENCE
@info "Building observation db..."
#generate the genomic sample dataframe
if isfile(sample_output) #if the sample DB has already been built for the current project, terminate
    @error "Existing sample dataframe at $sample_output"
else #otherwise, build it from scratch
    @info "Setting up sampling jobs..."
    #setup worker channels, input gets the genome partitions
    channels = setup_sample_jobs(danio_genome_path, danio_gen_index_path, danio_gff_path, sample_set_length, sample_window_min, sample_window_max, perigenic_pad)

    @info "Sampling genome..."
    sample_record_dfs=execute_sample_jobs(channels, worker_pool)

    @info "Serializing sample dataframes..."
    serialize(sample_output, sample_record_dfs) #write the dataframes to file
end

rmprocs(worker_pool)

@info "Done sampling jobs!"


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

@testset "Clebsch-Gordan" begin

@info("Testset Clebsch-Gordan")

#---

using PyCall, Test, ACE, ACE.SphericalHarmonics, JuLIP.Testing, StaticArrays
using JuLIP: evaluate
using ACE.SphericalHarmonics: index_y
using ACE.RPI.Rotations3D: ClebschGordan

#---

try
	sympy = pyimport("sympy")
	spin = pyimport("sympy.physics.quantum.spin")

	pycg(j1, m1, j2, m2, j3, m3, T=Float64) =
      		spin.CG(j1, m1, j2, m2, j3, m3).doit().evalf().__float__()

	cg = ClebschGordan()

	@info("compare implementation against `sympy`")
	let ntest = 0
   		while ntest < 200
      			j1 = rand(0:10)
      			j2 = rand(0:10)
      			J = rand(abs(j1-j2):min(10,(j1+j2)))
      			m1 = rand(-j1:j1)
      			for m2 = -j2:j2
         			M = m1+m2
         			if abs(M) <= J
            				ntest += 1
            				print_tf(@test cg(j1,m1, j2,m2, J,M) ≈ pycg(j1,m1, j2,m2, J,M))
            				# print_tf(@test clebschgordan(j1,m1, j2,m2, J,M) ≈ )
         			end
      			end
   		end
	end
	println()

	##

	@info("Checking the SphH expansion in terms of CG coeffs")
	# expansion coefficients of a product of two spherical harmonics in terms a
	# single spherical harmonic
	# see e.g. https://en.wikipedia.org/wiki/Clebsch–Gordan_coefficients
	# this is the magic formula that we need, on which everything else is based
	for ntest = 1:200
   		# two random Ylm  ...
   		l1, l2 = rand(1:10), rand(1:10)
   		m1, m2 = rand(-l1:l1), rand(-l2:l2)
   		# ... evaluated at random spherical coordinates
   		θ = rand() * π
   		φ = (rand()-0.5) * 2*π
   		R = SVector( cos(φ)*sin(θ), sin(φ)*sin(θ), cos(θ) )
   		# evaluate all relevant Ylms (up to l1 + l2)
   		Ylm = evaluate(SHBasis(l1+l2), R)
   		# evaluate the product p = Y_l1_m1 * Y_l2_m2
   		p = Ylm[index_y(l1,  m1)] * Ylm[index_y(l2,m2)]
   		# and its expansion in terms of CG coeffs
   		p2 = 0.0
   		M = m1 + m2  # all other coeffs are zero

   		for L = abs(M):(l1+l2)
      			p2 += sqrt( (2*l1+1)*(2*l2+1) / (4 * π * (2*L+1)) ) *
            		cg(l1,  0, l2,  0, L, 0) *
            		cg(l1, m1, l2, m2, L, M) *
            		Ylm[index_y(L, M)]
   		end
   		print_tf((@test (p ≈ p2) || (abs(p-p2) < 1e-15)))
	end
	println()

catch e   # try importing some python packages
	@warn "sympy didn't import, No Clebsch-Gordan tests are running"
end

#---

end # @testset


text:
#module phit

include("roundPerRound.jl")
#using .roundPerRound

#end


text:
y2 = @. sin(x^2) * exp(-x/10)
plot!(x, y2, col="blue", lwidth=0.05, label="plot 2")

xlabel("x-axis")
ylabel("y-axis")
xticks([-pi/2, 0, pi/2], ["π/2", "0", "π/2"])
ylim(-1.5, 1.5)
yticks(-1:0.25:1)

legend()


text:
# ft_setvidpid
export ft_setvidpid

function ft_setvidpid(vid::Integer, pid::Integer)
  @assert vid>=0
  @assert pid>=0
  ft_status = ccall((:FT_SetVIDPID, d2xx),
                     Cuint,
                     (Cuint, Cuint),
                     vid, pid)
  checkstatus(ft_status)
  return nothing
end


text:
# Implementation of sensitivities for unary linalg optimisations.
# Implementation of sensitivities for binary linalg optimisations.
const A = ∇Array
const S = ∇Scalar
const AS = Union{∇Scalar, ∇Array}

# Sensitivities for the Kronecker product:
import LinearAlgebra: kron
@explicit_intercepts kron Tuple{A, A}

# The allocating versions simply allocate and then call the in-place versions.
∇(::typeof(kron), ::Type{Arg{1}}, p, Y::A, Ȳ::A, A::A, B::A) =
    ∇(zeroslike(A), kron, Arg{1}, p, Y, Ȳ, A, B)
∇(::typeof(kron), ::Type{Arg{2}}, p, Y::A, Ȳ::A, A::A, B::A) =
    ∇(zeroslike(B), kron, Arg{2}, p, Y, Ȳ, A, B)

function ∇(Ā::A, ::typeof(kron), ::Type{Arg{1}}, p, Y::A, Ȳ::A, A::A, B::A)
    (I, J), (K, L), m = size(A), size(B), length(Y)
    for j = reverse(1:J), l = reverse(1:L), i = reverse(1:I)
        aij, āij = A[i, j], Ā[i, j]
        for k = reverse(1:K)
            āij += Ȳ[m] * B[k, l]
            m -= 1
        end
        Ā[i, j] = āij
    end
    return Ā
end
function ∇(B̄::A, ::typeof(kron), ::Type{Arg{2}}, p, Y::A, Ȳ::A, A::A, B::A)
    (I, J), (K, L), m = size(A), size(B), length(Y)
    for j = reverse(1:J), l = reverse(1:L), i = reverse(1:I)
        aij = A[i, j]
        for k = reverse(1:K)
            B̄[k, l] += Ȳ[m] * aij
            m -= 1
        end
    end
    return B̄
end

@explicit_intercepts Base.:+ Tuple{A, UniformScaling}
∇(::typeof(+), ::Type{Arg{1}}, p, Y::∇Array, Ȳ::∇Array, A::∇Array, B::UniformScaling) = Ȳ

@explicit_intercepts Base.:+ Tuple{UniformScaling, A}
∇(::typeof(+), ::Type{Arg{2}}, p, Y::∇Array, Ȳ::∇Array, A::UniformScaling, B::∇Array) = Ȳ

# `copy` materializes `Adjoint` and `Transpose` wrappers but can be called on anything
import Base: copy
@explicit_intercepts copy Tuple{Any}
∇(::typeof(copy), ::Type{Arg{1}}, p, Y, Ȳ, A) = copy(Ȳ)


text:
using TestingNewCodingStyle
using Statistics: mean, var
using StatsBase: fit, normalize, Histogram
using ProgressMeter
using Plots, FileIO
# using Revise

using Distributed
addprocs(3)
# addprocs(4, exeflags=`--threads 3`, enable_threaded_blas=true)

# machine3 = [("pc3@192.168.15.7",5)]
# addprocs(machine3;
#     multiplex=true,
#     exeflags=`--threads 2`,
#     enable_threaded_blas=true, dir="/home/pc3")

# @everywhere push!(LOAD_PATH, pwd())
# @everywhere using TestingNewCodes

fullPATH = homedir()*"/Documentos/TestingNewCodingStyle"
for p in workers()
    fetch(@spawnat p cd(fullPATH))
    fetch(@spawnat p push!(LOAD_PATH, pwd()))
end
@everywhere using TestingNewCodingStyle

@everywhere function get_intensity_sample(atoms, s, Δ, sensors)
    # atoms = Cube(N, kL)
    laser = Gaussian_3D(estimate_waist(atoms), s, Δ)
    simulation = ScalarProblem(atoms, laser)

    βₙ = get_steady_state(simulation)
    intensities = get_scattered_intensity(simulation, βₙ, sensors)
    GC.gc()
    return intensities
end

## data that I will use
ρλ³_range = range(5,45,length=5)
maxRepAvailable = 10
kL = 32.4;

# each density needs different repetitions to convergence
low_density_rep = maxRepAvailable
high_density_rep = maxRepAvailable/2
rep_range = round.(Int, range(low_density_rep,high_density_rep, length=length(ρλ³_range)))

s = 1e-6
nSensors = 360
sensors = ring_on_space(num_pts=nSensors, kR=1.5kL, θ=5π/12)
Δ_range = range(0,2,length=3)

variances = zeros(length(ρλ³_range), length(Δ_range))
for (iρ, ρλ³) in enumerate(ρλ³_range)
    ρ_over_k₀³ = ρλ³ / (2π)^3
    N = floor(Int, ρ_over_k₀³ * kL^3)

    atoms_saved = load(fullPATH*"/saved/Cubes_kL=$(kL), rho=$(ρλ³), maxRepAvailable=$(maxRepAvailable).jld2", "atoms_saved")
    maxRep = rep_range[iρ]

    for (iΔ, Δ) in enumerate(Δ_range)
        # println(iρ, " -- ", iΔ)
        println("density: $(iρ)/$(length(ρλ³_range)) -- detuning: $(iΔ)/$(length(Δ_range))" )
        # @showprogress
        many_intensities = pmap(1:maxRep) do rep
            intensities = get_intensity_sample(atoms_saved[rep], s, Δ, sensors)
        end

        all_intensities = Float64[]
        for i in many_intensities
            for j = 1:length(i)
                push!(all_intensities, i[j] )
            end
        end
        all_intensities_over_mean =  all_intensities./ mean(all_intensities)

        variances[iρ, iΔ] = var(all_intensities_over_mean)
    end
    # heatmap(Δ_range, ρλ³_range, variances,
    #     c=:jet, xlabel="Δ/Γ", ylabel="ρλ³",
    #     title="s=$(s), kL=$(kL)"
    # ) |> display
end

# save("heatmap_variances.jld2", Dict("variances" => variances))

using FileIO
variances = load("src/heatmap_variances.jld2", "variances")

ρλ³_range = range(5,45,length=5)
Δ_range = range(-1,2,length=5)

using Plots

#Δ_range, ρλ³_range,
heatmap( variances,
        c=:jet, xlabel="Δ/Γ", ylabel="ρλ³",
        title=""
    )
# savefig("florent_middle_resolution.png")


text:
##############################################################################
##
## Matrix constructors
##
##############################################################################

for (f, basef) in ((:dataeye, :eye), )
    @eval begin
        ($f)(n::Int) = DataArray(($basef)(n), falses(n, n))
        ($f)(n::Int, p::Int) = DataArray(($basef)(n, p), falses(n, p))
    end
end
for (f, basef) in ((:datadiagm, :diagm), )
    @eval begin
        ($f)(vals::Vector) = DataArray(($basef)(vals), falses(length(vals), length(vals)))
    end
end

##############################################################################
##
## Extract the matrix diagonal
##
##############################################################################

function diag{T}(dm::DataMatrix{T})
    return DataArray(diag(dm.data), diag(dm.na))
end

##############################################################################
##
## Size information
##
##############################################################################

nrow(dm::DataMatrix) = size(dm, 1)
ncol(dm::DataMatrix) = size(dm, 2)


text:
using BinaryBuilder

name = "PROPACK"
version = v"0.1"

sources = [
    ArchiveSource("https://github.com/optimizers/PROPACK/archive/v1.0.tar.gz",
                  "0d029a4c2cdcdb9b18a4fae77593a562f79406c3f79839ee948782b37974a10e")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/PROPACK-*/

OPENBLAS="-lopenblas"
FFLAGS=(-xf77-cpp-input)
if [[ ${nbits} == 64 ]] && [[ ${target} != aarch64* ]]; then
  OPENBLAS="-lopenblas64_"
  FFLAGS+=(-fdefault-integer-8 -fno-align-commons)

  for sym in caxpy cdotc cdotu ccopy cgemv clarnv clascl cscal csscal daxpy dbdsdc dbdsqr dcopy ddot dgemm dgemv dlamch dlapy2 dlarnv dlartg dlascl dnrm2 dznrm2 drot dscal lsame saxpy sbdsdc sbdsqr scnrm2 scopy sdot sgemm sgemv slamch slapy2 slarnv slartg slascl snrm2 srot sscal zaxpy zcopy zdotc zdotu zdscal zgemv zlarnv zlascl zscal
  do
    FFLAGS+=("-D${sym}=${sym}_64")
  done
fi

make SLIB=${dlext} FFLAG=\"${FFLAGS[*]}\" BLAS=${OPENBLAS}
cp complex8/libcpropack.${dlext} complex16/libzpropack.${dlext} single/libspropack.${dlext} double/libdpropack.${dlext} ${libdir}/
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

platforms = expand_gfortran_versions(platforms)

products = [
    LibraryProduct("libcpropack", :libcpropack),
    LibraryProduct("libzpropack", :libzpropack),
    LibraryProduct("libspropack", :libspropack),
    LibraryProduct("libdpropack", :libdpropack),
]

dependencies = [
    Dependency("OpenBLAS_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: kinesis_video_signaling
using AWS.Compat
using AWS.UUIDs

"""
    GetIceServerConfig()

Gets the Interactive Connectivity Establishment (ICE) server configuration information, including URIs, username, and password which can be used to configure the WebRTC connection. The ICE component uses this configuration information to setup the WebRTC connection, including authenticating with the Traversal Using Relays around NAT (TURN) relay server.  TURN is a protocol that is used to improve the connectivity of peer-to-peer applications. By providing a cloud-based relay service, TURN ensures that a connection can be established even when one or more peers are incapable of a direct peer-to-peer connection. For more information, see A REST API For Access To TURN Services.  You can invoke this API to establish a fallback mechanism in case either of the peers is unable to establish a direct peer-to-peer connection over a signaling channel. You must specify either a signaling channel ARN or the client ID in order to invoke this API.

# Required Parameters
- `ChannelARN`: The ARN of the signaling channel to be used for the peer-to-peer connection between configured peers.

# Optional Parameters
- `ClientId`: Unique identifier for the viewer. Must be unique within the signaling channel.
- `Service`: Specifies the desired service. Currently, TURN is the only valid value.
- `Username`: An optional user ID to be associated with the credentials.
"""
get_ice_server_config(ChannelARN; aws_config::AWSConfig=global_aws_config()) = kinesis_video_signaling("POST", "/v1/get-ice-server-config", Dict{String, Any}("ChannelARN"=>ChannelARN); aws_config=aws_config)
get_ice_server_config(ChannelARN, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_video_signaling("POST", "/v1/get-ice-server-config", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ChannelARN"=>ChannelARN), args)); aws_config=aws_config)

"""
    SendAlexaOfferToMaster()

This API allows you to connect WebRTC-enabled devices with Alexa display devices. When invoked, it sends the Alexa Session Description Protocol (SDP) offer to the master peer. The offer is delivered as soon as the master is connected to the specified signaling channel. This API returns the SDP answer from the connected master. If the master is not connected to the signaling channel, redelivery requests are made until the message expires.

# Required Parameters
- `ChannelARN`: The ARN of the signaling channel by which Alexa and the master peer communicate.
- `MessagePayload`: The base64-encoded SDP offer content.
- `SenderClientId`: The unique identifier for the sender client.

"""
send_alexa_offer_to_master(ChannelARN, MessagePayload, SenderClientId; aws_config::AWSConfig=global_aws_config()) = kinesis_video_signaling("POST", "/v1/send-alexa-offer-to-master", Dict{String, Any}("ChannelARN"=>ChannelARN, "MessagePayload"=>MessagePayload, "SenderClientId"=>SenderClientId); aws_config=aws_config)
send_alexa_offer_to_master(ChannelARN, MessagePayload, SenderClientId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = kinesis_video_signaling("POST", "/v1/send-alexa-offer-to-master", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ChannelARN"=>ChannelARN, "MessagePayload"=>MessagePayload, "SenderClientId"=>SenderClientId), args)); aws_config=aws_config)


text:
using XLSX, JSON
using DataFrames
using Statistics: mean
include("../src/JOH.jl")
include("MDMKP_id_df.jl")

function vec_to_DF(vec, problem_id)
	dataframe = JOH.Matheur.Model_DF()
	for row in vec
		push!(dataframe, row)
	end

	#add problem id to this dataframe
	id_col = repeat([problem_id], nrow(dataframe))
	dataframe[!, "problem_id"] = id_col

	dataframe
end

function read_problem_dir(folder, problem_id)
	filepath = joinpath(folder, problem_id)
	results = []

	phase = 0
	while true
		phase += 1

		model_path = joinpath(filepath, "$phase", "model.mps")
		results_path = joinpath(filepath, "$phase", "results.json")
		sol_path = joinpath(filepath, "$phase", "start_sol.json")
		try
			results_file = open(results_path, "r")
			phase_results = read(results_file, String)
			close(results_file)

			push!(results, phase_results)
		catch SystemError
			break
		end
	end

	vec_to_DF(JSON.parse.(results), problem_id)
end

function summarize_df(df)
	su = DataFrame()
	su[!, "id"] = [parse(Int, df.problem_id[1])]
	highest_rtol = -1.0
	lowest_gap = -1.0
	for phase in 1:nrow(df)
		su[!, "p$(phase)_obj"] = [-1*df.objective[phase]]
		su[!, "p$(phase)_gap"] = [df.gap[phase]]
		su[!, "p$(phase)_dettime"] = [df.dettime[phase]]
		su[!, "p$(phase)_solve_time"] = [df.solve_time[phase]]
		su[!, "p$(phase)_elapsed_time"] = [df.elapsed_time[phase]]
		su[!, "p$(phase)_rtol"] = [df.rtol[phase]]
		su[!, "p$(phase)_term"] = [df.term_stat[phase]]

		highest_rtol = df.rtol[phase]
		lowest_gap = df.gap[phase]
	end
	su[!, "total_solve_time"] = [sum(df.solve_time)]
	su[!, "total_elapsed_time"] = [sum(df.elapsed_time)]
	su[!, "total_dettime"] = [sum(df.dettime)]
	su[!, "highest_rtol"] = [highest_rtol]
	su[!, "lowest_gap"] = [lowest_gap]
	su
end

function read_result_dir(folder)
	results = []
	summaries = []
	for problem in readdir(folder)
		df = read_problem_dir(folder, problem)
		println(df)
		push!(results, df)
		push!(summaries, summarize_df(df))
	end
	vcat(results...), vcat(summaries..., cols=:union)
end

function join_summary_problem_tables(summary_df, problem_df)
	new_results = innerjoin(problem_df, summary_df, on="id")
	new_results
end

function category_section_table(summary_df, category)
	filtered = filter(row->row.category==category, summary_df)
	filtered
end

function load_summary(folder)
	#make complete and summary dataframes
	complete, summary = read_result_dir(folder)
	problem_df = MDMKP_id_dataframe()
	complete, join_summary_problem_tables(summary, problem_df)
end

function case_dataset_pivot(summary, category)
	f = category_section_table(summary, category)
	groups = groupby(f, [:case, :dataset])
	vasko_summary = combine(groups, nrow => :n_problems,
		:total_elapsed_time => mean => :aver_run_time,
		:highest_rtol => mean => :aver_tol,
		:lowest_gap => mean => :aver_gap)

	XLSX.writetable(
		"category_$(category)_grouped_stats.xlsx",
		results =(
			collect(DataFrames.eachcol(vasko_summary)),
			DataFrames.names(vasko_summary) ))
end

function write_excel_decision_tree_results(; folder="results/decision_tree/")
	complete, summary = load_summary(folder)

	XLSX.writetable(
		"test(decision_tree_results).xlsx",
		results =(
			collect(DataFrames.eachcol(summary)),
			DataFrames.names(summary) ),
		complete_results =(
			collect(DataFrames.eachcol(complete)),
			DataFrames.names(complete) ))
end

function add_column(df, col, data)
	replace(col, " "=>"_")
	println(df, col, data)
	df[!, Symbol(col)] = repeat([data], nrow(df))
	df
end

function load_song_results(; folder="results/full_song_mdmkp")

	load_df(file) = load_summary(file)[2]
	load_method_df(method) = add_column(load_df(method), "method",
		last(splitpath(method)))

	methods = readdir(folder, join=true)
	methods = filter(x->isdir(x), methods)
	dataframes = load_method_df.(methods)

    data = vcat(dataframes..., cols=:union)
	data
end

df = load_song_results(folder="results/decision_tree/D")
XLSX.writetable("decision_tree_D.xlsx", collect(DataFrames.eachcol(df)), DataFrames.names(df))

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
using SOtoZulip
using DBInterface

include("../configuration.jl")

db = getdb(SODB)

########################################
# UP
########################################

queries = [
"""
ALTER TABLE questions
ADD COLUMN deleted INTEGER NOT NULL DEFAULT 0
""",
"""
ALTER TABLE answers
ADD COLUMN deleted INTEGER NOT NULL DEFAULT 0
"""
]

DBInterface.execute(db, "BEGIN TRANSACTION")
try
    DBInterface.execute.(Ref(db), queries)
    DBInterface.execute(db, "COMMIT TRANSACTION")
catch
    DBInterface.execute(db, "ROLLBACK")
end


text:
using Serialization, Statistics, DataFrames
import GEAS:MaCS, macs_2_base, macs_args, qtl_breeding, sim_QTL

function jenko_fork(fork, qtl, dic, nedit, esires, df, class)
    for t in nedit
        cc = "$class, $t edit"
        @info cc
        dic[:nEdit] = t         # how many (top) QTN to edit
        arg = (; dic...)
        base = deserialize(fork)
        ped, _ = qtl_breeding(base, arg, qtl, eSires = esires)
        tdf = combine(ped, :tbv => mean => :mtbv)
        tdf.sel .= cc
        append!(df, tdf)
    end
end

"""
    function jenko2015(; nid = 1000, nchr = 10, nsnp = 1000, maf = .05, dir = "dat/jenko2015")
---
## Description
This is to repeat Jenko 2015.
"""
function jenko2015(; nid = 1000, nchr = 10, nsnp = 1000, maf = .05, dir = "dat/jenko2015")
    @info "Simulation of a base with `MaCS`"
    par = macs_args(nid)
    tmp = joinpath(dir, "macs")
    isdir(tmp) || mkpath(tmp)

    # MaCS(par.cattle, nchr, dir = tmp)
    # base = macs_2_base(10, nsnp, dir = tmp)
    # serialize("$dir/base.ser", base)
    base = deserialize("$dir/base.ser")
    qtl = sim_QTL(base, nchr * nsnp)[1]

    dic = Dict(
        :nSire=> 25,
        :nDam => 500,
        :nSib => 2,             # such that selection is only on sires
        :nG8n => 21,
        :h²   => .5)            # dummy, only to utilize GEAS

    @info "Genomic selection on generation -20 → 0"

    arg = (; dic...)
    # Note, orginal work was just using TBV. No evaluation method envolved.
    ped, snp = qtl_breeding(base, arg, qtl)
    df = combine(ped, :tbv => mean => :mtbv)
    df.g8n[:] = arg.nG8n-1 : -1 : 0
    df.sel .= "pre-fork"        # 0 for generation -20 → 0

    # create new base, and insert editing option
    base.hap[:, :] = snp    # the new base, only to replace genotypes
    serialize("$dir/fork.ser", base) # this base will go through many different scenarios

    @info "Generation 1 → 20"
    dic[:nG8n] = 20

    @info "Continue selection with no editing"
    arg = (; dic...)
    ped, _ = qtl_breeding(base, arg, qtl)
    tdf = combine(ped, :tbv => mean => :mtbv)
    tdf.sel .= "GS only"     # 1 for continued normal selection from 0
    append!(df, tdf)

    # Editing the top 5 bulls, 25, 50, 100 QTN edits
    nedit = [25, 50, 100]
    esires = [ones(Int, 5); zeros(Int, 20)]
    jenko_fork("$dir/fork.ser", qtl, dic, nedit, esires, df, "top 5 sires")

    # Editing the top 10 bulls, 1, 5, 10, 20 QTN edits
    nedit = [1, 5, 10, 20]
    esires = [ones(Int, 10); zeros(Int, 15)]
    jenko_fork("$dir/fork.ser", qtl, dic, nedit, esires, df, "top 10 sires")

    # Editing the bottom 10 bulls, 1, 5, 10, 20 QTN edits
    esires = [zeros(Int, 15); ones(Int, 10)]
    jenko_fork("$dir/fork.ser", qtl, dic, nedit, esires, df, "bottom 10 sires")

    @info "Editing the all 25 bulls, 1, 5, 10 20 edits"
    esires = ones(Int, 25)
    jenko_fork("$dir/fork.ser", qtl, dic, nedit, esires, df, "all sires")

    df
end

function mytest(dic)
    dic[:tmp] => [1, 2, 3]
end


text:
using Test
using BenchmarkTools
include("../src/Optimizers.jl")

@testset "Least squares - simple" begin
    A = 1.0 .+ rand(3,3)
    b = rand(3)
    _optimize(A,b,LSTSQ())
end

@testset "Least squares - simple HD" begin
    A = 1.0 .+ rand(3,3)
    b = rand(2,3)
    _optimize(A,b,LSTSQ())
end

@testset "Thresholded Least squares - matrix" begin
    A = 1.0 .+ rand(3,3)
    b = rand(2,3)
    _optimize(A,b,STLSQ(1.0))
end


text:
using Documenter
using Jute

makedocs(
    modules = [Jute],
    format = Documenter.HTML(prettyurls=false),
    sitename = "Jute.jl",
    authors = "Bogdan Opanchuk",
    pages = [
        "Home" => "index.md",
        "Manual" => "manual.md",
        "Public API" => "public.md",
        "Internals" => "internals.md",
        "Version history" => "history.md",
    ],
)

deploydocs(
    repo = "github.com/fjarri/Jute.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
# Here we go through an example using the Lorenz equations
using Random, Distributions, LinearAlgebra
const plotit = false
Random.seed!(1234)

# Initialize arrays for timestepping
ṡ = [[0.0] for i = 1:3]
s = [[randn() + 10] for i = 1:3]
s̃ = copy(s)

function lorenz!(ṡ, s; ρ = 28, β = 8 / 3, σ = 10)
    ẋ, ẏ, ż = ṡ
    x, y, z = s
    @. ẋ = σ * (-x + y)
    @. ẏ = -y + (ρ - z) * x
    @. ż = -β * z + x * y
    return nothing
end

function closure_lorenz!(ρ, β, σ)
    function rhs!(ṡ, s)
        lorenz!(ṡ, s, ρ = ρ, β = β, σ = σ)
    end
end

function step!(s, ṡ, s̃, rhs!, Δt)
    rhs!(ṡ, s)
    @. s̃ = s + Δt * ṡ
    @. s = s + Δt * 0.5 * ṡ
    rhs!(ṡ, s̃)
    @. s = s + Δt * 0.5 * ṡ
    return nothing
end

# timestep and record solution
timeseries = []
push!(timeseries, [(s...)...])
Δt = 0.005
N = 2^19 - 1

for i = 1:N
    step!(s, ṡ, s̃, lorenz!, Δt)
    push!(timeseries, [(s...)...])
end

x = [timeseries[i][1] for i in eachindex(timeseries)]
y = [timeseries[i][2] for i in eachindex(timeseries)]
z = [timeseries[i][3] for i in eachindex(timeseries)]
#=
if plotit
using GLMakie, Printf, LaTeXStrings
fig = Figure(resolution = (1600, 1300))

averaging_slider = Slider(fig, range = 1:16, startvalue = 1)
i = averaging_slider.value
x̅ₜ = @lift(mean(reshape(x, (2^$i, 2^(19 - $i))), dims = 1)[:])
time_string = @lift("Averaging Interval T = " * @sprintf("%0.3f", $i * Δt))

ax = Axis(fig[1, 1],
    xlabel = L" \frac{1}{\tau+T}\int_{\tau}^T x(t) dt ", xlabelsize = 40,
    ylabel = "pdf", ylabelsize = 40,
    xticklabelsize = 30, yticklabelsize = 30)

@lift ylims!(ax, (0, 2 / sqrt(2 * π * var($x̅ₜ))))
xlims!(ax, (-23, 23))

d1 = density!(ax, x̅ₜ, color = (:purple, 0.15), strokewidth = 1)

fig[2, 1] = vgrid!(
    Label(fig, time_string, width = nothing, textsize = 30),
    averaging_slider,
)

display(fig)

frames = 16
fps = 5
record(fig, pwd() * "/lorenz_ekixample.mp4"; framerate = fps) do io
    for ii = 1:frames
        i[] = ii
        sleep(1 / fps)
        recordframe!(io)
    end
end
end
=#
##
function closure_forward_map(i)
    function forward_map(C)
        ρ, σ, β = C
        rhs! = closure_lorenz!(ρ, β, σ)
        Δt = 0.005
        N = 2^16 - 1

        Random.seed!(1234)
        # Initialize arrays for timestepping
        ṡ = [[0.0] for i = 1:3]
        s = [[randn() + 10] for i = 1:3]
        s̃ = copy(s)
        timeseries = []
        push!(timeseries, [(s...)...])

        for i = 1:N
            step!(s, ṡ, s̃, rhs!, Δt)
            push!(timeseries, [(s...)...])
        end

        x = [timeseries[i][1] for i in eachindex(timeseries)]
        y = [timeseries[i][2] for i in eachindex(timeseries)]
        z = [timeseries[i][3] for i in eachindex(timeseries)]

        # use the Nusselt number as the observation
        Nuₜ = mean(reshape(z .* β, (2^i, 2^(16 - i))), dims = 1)[:]
        z² = mean(reshape(z .* z, (2^i, 2^(16 - i))), dims = 1)[:]
        z³ = mean(reshape(z .* z .* z, (2^i, 2^(16 - i))), dims = 1)[:]
        x²z = mean(reshape(x .* x .* z, (2^i, 2^(16 - i))), dims = 1)[:]
        return [Nuₜ, z²,]
    end
end

𝒢 = closure_forward_map(12)
C = [28, 10, 8 / 3]

forward_map(C) = mean.(𝒢(C))

y = 𝒢(C)
y̅ = mean.(y)
Γ = Diagonal(var.(y))
J = length(C) * 10
ξ = MvNormal(zeros(length(y)), Γ)
# Guassianize

prior = MvNormal([28, 10, 8 / 3], Diagonal([3.0, 3.0, 1.0]))
prior = MvNormal([20, 5, 4], Diagonal([5.0, 4.0, 1.0]))
# Construct empirical prior
u = [rand(prior) for i = 1:J]
u₀ = copy(u)

timeseries = []
push!(timeseries, copy(u))

N = 4
h = 1 / N
⊗(a, b) = a * b'
for i = 1:N
    println("loop ", i)
    u̅ = mean(u)
    G = forward_map.(u) # error handling needs to go here
    G̅ = mean(G)

    # define covariances
    Cᵘᵖ = (u[1] - u̅) ⊗ (G[1] - G̅)
    Cᵖᵖ = (G[1] - G̅) ⊗ (G[1] - G̅)
    for j = 2:J
        Cᵘᵖ += (u[j] - u̅) ⊗ (G[j] - G̅)
        Cᵖᵖ += (G[j] - G̅) ⊗ (G[j] - G̅)
    end
    Cᵘᵖ *= 1 / (J - 1)
    Cᵖᵖ *= 1 / (J - 1)

    # ensemblize the data
    y = [y̅ + rand(ξ) for i = 1:J]
    r = y - G

    # update
    if length(Γ) > 1
        Cᵖᵖ_factorized = cholesky(Symmetric(Cᵖᵖ + 1 / h * Γ))
    else
        Cᵖᵖ_factorized = Cᵖᵖ + 1 / h * Γ
    end
    for j = 1:J
        u[j] += Cᵘᵖ * (Cᵖᵖ_factorized \ r[j])
    end
    push!(timeseries, copy(u))
end

residual = C - mean(timeseries[end])

println("The relative error of ρ is ", abs(residual[1]) / C[1])
println("The relative error of σ is ", abs(residual[2]) / C[2])
println("The relative error of β is ", abs(residual[3]) / C[3])

max_stars_repo_name: bischtob/SEC

max_issues_repo_name: bischtob/SEC

max_forks_repo_name: bischtob/SEC


text:
using CliOptions

spec = CliOptionSpec(
    HelpOption(),
    # FlagOption("-0", "--null";
    #            help = "Use NUL character as delimiter."),
    Positional("command"; default = "echo",
               help = "The command to execute for each items."),
    Positional("arguments"; default = String[], multiple = true,
               help = "The arguments to be passed to the commands."),
    use_double_dash = true,
)
options = parse_args(spec)

for line in eachline(stdin)
    cmd = Cmd([options.command, options.arguments..., line])
    run(cmd)
end