text:
"""
    add_overlap!(ata_model::AbstractModel; overlap_file = "overlap_matrix.csv", overlap_delim=";")

# Description

Add maximum overlap constraints to the `ATA.AbstractModel` as specified by the `overlap_file` which contains a `TxT` matrix defining the maximum number of shared items between each pair of tests.
Alternatively, the overlap matrix can be passed using the argument `overlap`.
It requires the [`start_ata`](#ATA.start_ata) step.

# Arguments

- **`ata_model::AbstractModel`** : Required. The model built with `start_ata()` and with settings loaded by [`start_ata`](#ATA.start_ata) function.
- **`overlap::Matrix{Float64} = Matrix{Float64}(undef, 0, 0)`**: Optional.
- **`overlap_file`** : Optional. Default: "overlap_matrix.csv". The path of the file containing the maximum overlap between each pair of tests in the form of a matrix with custom-separated values.
- **`overlap_delim`** : Optional. Default: ";". The custom-separator for the `overlap_file`.

"""
function add_overlap!(
    ata_model::AbstractModel;
    overlap::Matrix{Float64} = Matrix{Float64}(undef, 0, 0),
    overlap_file = "overlap_matrix.csv",
    overlap_delim = ";",
)

    #load overlap matrix
    if size(overlap, 1) > 0
        overlap_matrix = overlap
        success!(ata_model, "Overlap matrix loaded.")
    elseif isfile(overlap_file)
        try
            overlap_matrix = CSV.read(
                overlap_file,
                DataFrames.DataFrame,
                delim = overlap_delim,
                header = false,
            )
        catch e
            error!(ata_model, "Error in reading the overlap file.")
            return nothing
        end
        success!(ata_model, "Overlap file loaded.")
    else
        error!(
            ata_model,
            string(
                "Overlap file with name ",
                overlap_file,
                " does not exist. Provide a valid overlap matrix or a name of an existing file.",
            ),
        )
        return nothing
    end
    try
        T = ata_model.settings.T
        n_items = ata_model.settings.n_items

        if size(overlap_matrix, 1) > 0
            #ol_max = Vector{Vector{Int64}}(undef, T)
            # for t = 1:T
            # 	ol_max[t] = overlap_matrix[t, setdiff(collect(1:T), t)]
            # end
            overlap_matrix = Matrix{Float64}(overlap_matrix)[1:T, 1:T]
            DelimitedFiles.writedlm("opt/overlap.txt", overlap_matrix)
            ata_model.settings.ol_max = overlap_matrix
            ata_model.settings.to_apply[3] = true
        else
            ata_model.settings.ol_max =
                ones(ata_model.settings.n_items, ata_model.settings.n_items) .*
                ata_model.settings.n_items
            ata_model.settings.to_apply[3] = false
            success!(
                ata_model,
                string(
                    "No lines in ",
                    overlap_file,
                    ", overlap constraints are not applied.",
                ),
            )
        end
        JLD2.@save "opt/ata_model.jld2" ata_model
        success!(ata_model, "Maximum overlap constrained.")
    catch e
        error!(ata_model, string(sprint(showerror, e)))

    end
    return nothing
end


text:
#### Problem generation ####

m=[2,3,4,5,6]; # Problem size
I=10.*m; #number of costumers
J=m+1; # number of service facilities
K=3(m+1) #number of potential locations for defence facilities

function Ecualidian_distanse(x1::Float64,x2::Float64,y1::Float64,y2::Float64)
    dist=((x1-x2)^2+(y1-y2)^2)^(1/2)
    return dist
end

# customer locations
customer_location=zeros(I[1],2)
for i=1:I[1]
    customer_location[i,1]=1000*rand()
    customer_location[i,2]=1000*rand()
end

#service locations
service_location=zeros(J[1],2)
for j=1:J[1]
    service_location[j,1]=1000*rand()
    service_location[j,2]=1000*rand()
end

#potential defence positions
pdef_location=zeros(K[1],2)
for k=1:K[1]
    pdef_location[k,1]=1000*rand()
    pdef_location[k,2]=1000*rand()
end

distance_cus_ser=zeros(I[1],J[1])# distance between customer i and service facility j

for i=1:I[1]
    for j=1:J[1]
        dist=Ecualidian_distanse(customer_location[i,1],service_location[j,1],customer_location[i,2],service_location[j,2])
        distance_cus_ser[i,j]=dist
    end
end

distance_ser_pdef=zeros(J[1],K[1]) # distance between service facility j and potential position for defence facilities
for j=1:J[1]
    for k=1:K[1]
        dist=Ecualidian_distanse(service_location[j,1],pdef_location[k,1],service_location[j,2],pdef_location[k,2])
        distance_ser_pdef[j,k]=dist
    end
end

customer_demand=zeros(1,I[1]) # defining customer demand
for i=1:I[1]
    customer_demand[i]=rand(collect(5:5:100))
end

def_pos_cost=zeros(1,K[1]) # cost of stablishing a defence facility at position k
for k=1:K[1]
    def_pos_cost[k]=rand(collect(1000:1000:20000))
end

attacker_buddet=2; # attacker budget ### Must change based on  the paper
 r=60  #range of defence facilities

 #service facilities capacity
 ser_cap=zeros(1,J[1]);
 q_p_vec=zeros(1,J[1]);
for j=1:J[1]
    q_p_vec[j]=rand([1,m[1]+1])
end
su=(1+rand())*sum(customer_demand)
for j=1:J[1]
    ser_cap[j]=(su*q_p_vec[j])/sum(q_p_vec)
end
for j=1:J[1]
    a=rem(ser_cap[j],20)
    ser_cap[j]=ser_cap[j]-a
end

#defining probablity of defence
prob_matrix=zeros(size(distance_ser_pdef))
for j=1:J[1]
    p_l=0
    p_u=0
    for k=1:K[1]
        if 0<distance_ser_pdef[j,k]<=r
            p_u=0.95
        else
            p_u=0.15
        end
        if distance_ser_pdef[j,k]>r
            p_l=0.15
        else
            p_l=0.95
        end
        prob_matrix[j,k]=p_l+(((p_u-p_l)*(r-distance_ser_pdef[j,k]))/r)
    end
end

CS=2 #unit shipment cost of demand per unit distance
CO=500 # unit outsourcing cost of demand (independent of distance)

ext: jl

max_stars_repo_name: imanpanahi/MilProject

max_issues_repo_name: imanpanahi/MilProject

max_forks_repo_name: imanpanahi/MilProject


text:
## Residual Constrained Alternating Minimization
# Function Definitions

function AltMinChol{ET<:Number}(b::AbstractArray{ET}, r::Int, K::Int)
# Performs RCAM on a single matricised monochromatic slice
# Calls: minFrobQR

    # Generate L0
    m,n = size(b)
    L = collect(ET, randn(m,r))

    indi = (1:max(m,n))::UnitRange{Int64}
    local R

    # Begin Alterations
    for k = 1:K

        # Solve for R
     @time R = minFrobChol(L,b,indi)

        # Solve for L
     @time L = minFrobChol(R,b',indi)

    end

    return L, R
end
export AltMinChol

function minFrobChol{ET<:Number}(L::AbstractArray{ET}, b::AbstractArray{ET}, indi::UnitRange{Int64})
# Solves for each row of input
# Calls: minFrobRsub

    m = size(b,2)::Int
    r = size(L,2)::Int

    # Initialize R
    R = zeros(ET,m,r)

    for k =1:length(indi)

        if indi[k]<=m

             R[indi[k],:] = minFrobRsubChol(L, b, indi[k])

        end

    end

    R = conj(R)

    return R
end

function minFrobRsubChol{ET<:Number}(L::AbstractArray{ET},b::AbstractArray{ET},indi::Int)
# Solves a single row

    # Solve for the kth row of R
    bsub = b[:,indi]

    # Find the non zero rows of the kth column of b
    ind = find(x -> x != 0, b[:,indi])

    # Keep only the non zero
    L2 = L[ind,:]
    bsub = transpose(L2)*bsub[ind]

    # Chol fact
    L2 = L2'*L2
    Q = chol(L2)

    bsub2 = fwdsub(Q',bsub)

    v = backsub(Q,bsub2)

    return v

end

function backsub{ET<:Number}(A::AbstractArray{ET}, y::AbstractArray{ET})

    m = size(A,2)
    x = zeros(ET,m,1)

    x[m] = y[m]/A[m,m]

    for i = m-1:-1:1
        x[i] = (y[i]-dot(A[i,i+1:m], x[i+1:m]))/A[i,i]
    end

    return x
end

function fwdsub{ET<:Number}(A::AbstractArray{ET},y::AbstractArray{ET})

    m = size(A,2)
    x = zeros(ET,m,1)
    x[1] = y[1]/A[1,1]

    for i = 2:1:m
        x[i] = (y[i]-dot(A[i,i-1:-1:1], x[i-1:-1:1]))/A[i,i]
    end

    return x
end


text:
"""
ConfusionMatrix{C}

Confusion matrix with `C ≥ 2` classes. Rows correspond to predicted values
and columns to the ground truth.
"""
struct ConfusionMatrix{C}
    mat::Matrix
    labels::Vector{String}
end

"""
ConfusionMatrix(m, labels)

Instantiates a confusion matrix out of a square integer matrix `m`.
Rows are the predicted class, columns the ground truth. See also
the [wikipedia article](https://en.wikipedia.org/wiki/Confusion_matrix).
"""
function ConfusionMatrix(m::Matrix{Int}, labels::Vector{String})
    s = size(m)
    s[1] == s[2] || throw(ArgumentError("Expected a square matrix."))
    s[1] > 1 || throw(ArgumentError("Expected a matrix of size ≥ 2x2."))
    length(labels) == s[1] || throw(ArgumentError("As many labels as classes must be provided."))
    ConfusionMatrix{s[1]}(m, labels)
end

# allow to access cm[i,j] but not set (it's immutable)
Base.getindex(cm::ConfusionMatrix, inds...) = getindex(cm.mat, inds...)

"""
confusion_matrix(ŷ, y; rev=false)

Computes the confusion matrix given a predicted `ŷ` with categorical elements
and the actual `y`. Rows are the predicted class, columns the ground truth.
The ordering follows that of `levels(y)`.

## Keywords

* `rev=false`: in the binary case, this keyword allows to swap the ordering of
               classes.
* `perm=[]`:   in the general case, this keyword allows to specify a permutation
               re-ordering the classes.
* `warn=true`: whether to show a warning in case `y` does not have scientific type
               `OrderedFactor{2}` (see note below).

## Note

To decrease the risk of unexpected errors, if `y` does not have scientific type
`OrderedFactor{2}` (and so does not have a "natural ordering" negative-positive),
a warning is shown indicating the current order unless the user specifies,
explicitly either `rev` or `perm` in which case it's assumed the user is aware of
the class ordering.
"""
function confusion_matrix(ŷ::VC, y::VC;
                          rev::Union{Nothing,Bool}=nothing,
                          perm::Union{Nothing,Vector{<:Integer}}=nothing,
                          warn::Bool=true
                          ) where VC <: AbstractVector{<:CategoricalElement}
    check_dimensions(ŷ, y)
    levels_ = levels(y)
    nc = length(levels_)
    if rev !== nothing && rev && nc > 2
        throw(ArgumentError("Keyword `rev` can only be used in binary case."))
    end
    if perm !== nothing && !isempty(perm)
        length(perm) == nc || throw(ArgumentError("`perm` must be of length matching the number of classes."))
        Set(perm) == Set(collect(1:nc)) || throw(ArgumentError("`perm` must specify a valid permutation of the classes."))
    end
    # warning
    if rev === nothing && perm === nothing
        if warn &&
            if nc==2 && !(scitype_union(y) >: OrderedFactor{2})
                @warn "The classes are un-ordered,\n" *
                      "using: negative='$(levels_[1])' and positive='$(levels_[2])'.\n" *
                      "To suppress this warning, consider coercing to OrderedFactor."
            elseif !(scitype_union(y) >: OrderedFactor{nc})
                @warn "The classes are un-ordered,\n" *
                      "using order: $([l for l in levels_]).\n" *
                      "To suppress this warning, consider coercing to OrderedFactor."
            end
        end
        rev  = false
        perm = Int[]
    elseif rev !== nothing
        # rev takes precedence in binary case
        if rev
            perm = [2, 1]
        else
            perm = Int[]
        end
    end

    # No permutation
    if isempty(perm)
        cmat = zeros(Int, nc, nc)
        @inbounds for i in eachindex(y)
            cmat[int(ŷ[i]), int(y[i])] += 1
        end
        return ConfusionMatrix(cmat, string.(levels_))
    end

    # With permutation
    cmat = zeros(Int, nc, nc)
    iperm = invperm(perm)
    @inbounds for i in eachindex(y)
        cmat[iperm[int(ŷ[i])], iperm[int(y[i])]] += 1
    end
    return ConfusionMatrix(cmat, string.(levels_))
end

# synonym
confmat = confusion_matrix

# Machinery to display the confusion matrix in a non-confusing way
# (provided the REPL is wide enough)

splitw(w::Int) = (sp1 = div(w, 2); sp2 = w - sp1; (sp1, sp2))

function Base.show(stream::IO, m::MIME"text/plain", cm::ConfusionMatrix{C}) where C
    width    = displaysize(stream)[2]
    cw       = 13
    textlim  = 9
    totalwidth = cw * (C+1) + C + 2
    width < totalwidth && (show(stream, m, cm.mat); return)

    iob     = IOBuffer()
    wline   = s -> write(iob, s * "\n")
    splitcw = s -> (w = cw - length(s); splitw(w))
    cropw   = s -> length(s) > textlim ? s[1:prevind(s, textlim)] * "…" : s

    # 1.a top box
    " "^(cw+1) * "┌" * "─"^((cw + 1) * C - 1) * "┐" |> wline
    gt = "Ground Truth"
    w  = (cw + 1) * C - 1 - length(gt)
    sp1, sp2 = splitw(w)
    " "^(cw+1) * "│" * " "^sp1 * gt * " "^sp2 * "│" |> wline
    # 1.b separator
    "┌" * "─"^cw * "┼" * ("─"^cw * "┬")^(C-1) * "─"^cw * "┤" |> wline
    # 2.a description line
    pr = "Predicted"
    sp1, sp2 = splitcw(pr)
    partial = "│" * " "^sp1 * pr * " "^sp2 * "│"
    for c in 1:C
        # max = 10
        s = cm.labels[c] |> cropw
        sp1, sp2 = splitcw(s)
        partial *= " "^sp1 * s * " "^sp2 * "│"
    end
    partial |> wline
    # 2.b separating line
    "├" * "─"^cw * "┼" * ("─"^cw * "┼")^(C-1) * ("─"^cw * "┤") |> wline
    # 2.c line by line
    for c in 1:C
        # line
        s  = cm.labels[c] |> cropw
        sp1, sp2 = splitcw(s)
        partial = "│" * " "^sp1 * s * " "^sp2 * "│"
        for r in 1:C
            e = string(cm[c, r])
            sp1, sp2 = splitcw(e)
            partial *= " "^sp1 * e * " "^sp2 * "│"
        end
        partial |> wline
        # separator
        if c < C
            "├" * "─"^cw * "┼" * ("─"^cw * "┼")^(C-1) * ("─"^cw * "┤") |> wline
        end
    end
    # 2.d final line
    "└" * "─"^cw * "┴" * ("─"^cw * "┴")^(C-1) * ("─"^cw * "┘") |> wline
    write(stream, take!(iob))
end


text:
include("test_utils.jl")

A   = preallocated_LinearOperator(get_div_grad(32, 32, 32))
n   = size(A, 1)
b   = ones(n)
M   = nonallocating_opEye(n)
mem = 10

# without preconditioner and with Ap preallocated, CG needs 3 n-vectors: x, r, p
storage_cg(n) = 3 * n
storage_cg_bytes(n) = 8 * storage_cg(n)

expected_cg_bytes = storage_cg_bytes(n)
cg(A, b, M=M)  # warmup
actual_cg_bytes = @allocated cg(A, b, M=M)
@test actual_cg_bytes ≤ 1.1 * expected_cg_bytes

# without preconditioner and with Ap preallocated, DIOM needs:
# - 2 n-vectors: x, x_old
# - 2 (n*mem)-matrices: P, V
# - 1 mem-vector: L
# - 1 (mem+2)-vector: H
# - 1 mem-bitArray: p
storage_diom(mem, n) = (2 * n) + (2 * n * mem) + (mem) + (mem + 2) + (mem / 64)
storage_diom_bytes(mem, n) = 8 * storage_diom(mem, n)

expected_diom_bytes = storage_diom_bytes(mem, n)
diom(A, b, M=M, memory=mem)  # warmup
actual_diom_bytes = @allocated diom(A, b, M=M, memory=mem)
@test actual_diom_bytes ≤ 1.05 * expected_diom_bytes

# with Ap preallocated, CG-Lanczos needs 4 n-vectors: x, v, v_prev, p
storage_cg_lanczos(n) = 4 * n
storage_cg_lanczos_bytes(n) = 8 * storage_cg_lanczos(n)

expected_cg_lanczos_bytes = storage_cg_lanczos_bytes(n)
cg_lanczos(A, b)  # warmup
actual_cg_lanczos_bytes = @allocated cg_lanczos(A, b)
@test actual_cg_lanczos_bytes ≤ 1.1 * expected_cg_lanczos_bytes

# without preconditioner and with Ap preallocated, DQGMRES needs:
# - 1 n-vectors: x
# - 2 (n*mem)-matrices: P, V
# - 2 mem-vectors: c, s
# - 1 (mem+2)-vector: H
storage_dqgmres(mem, n) = (n) + (2 * n * mem) + (2 * mem) + (mem + 2)
storage_dqgmres_bytes(mem, n) = 8 * storage_dqgmres(mem, n)

expected_dqgmres_bytes = storage_dqgmres_bytes(mem, n)
dqgmres(A, b, M=M, memory=mem)  # warmup
actual_dqgmres_bytes = @allocated dqgmres(A, b, M=M, memory=mem)
@test actual_dqgmres_bytes ≤ 1.05 * expected_dqgmres_bytes

# without preconditioner and with Ap preallocated, CR needs 4 n-vectors: x, r, p, q
storage_cr(n) = 4 * n
storage_cr_bytes(n) = 8 * storage_cr(n)

expected_cr_bytes = storage_cr_bytes(n)
cr(A, b, M=M)  # warmup
actual_cr_bytes = @allocated cr(A, b, M=M)
@test actual_cr_bytes ≤ 1.1 * expected_cr_bytes

# without preconditioner and with Ap preallocated, CGS needs 5 n-vectors: x, r, u, p, q
storage_cgs(n) = 5 * n
storage_cgs_bytes(n) = 8 * storage_cgs(n)
expected_cgs_bytes = storage_cgs_bytes(n)
cgs(A, b, M=M)  # warmup
actual_cgs_bytes = @allocated cgs(A, b, M=M)
@test actual_cgs_bytes ≤ 1.1 * expected_cgs_bytes


text:
# abstract type Components end
# abstract type NonDispatchables <: Components end
abstract type NonDispatchables end

"Project information."
struct Project
    "lifetime (years)"
    lifetime
    "discount rate ∈ [0,1]"
    discount_rate
    "time step (h)"
    timestep
    # TODO dispatch_type?
end

"Diesel generator parameters."
struct DieselGenerator
    "Rated power (kW)"
    power_rated   # decision variable
    "Minimum load ratio ∈ [0,1]"
    minimum_load_ratio  # ever it is on, it will work at least `min_load_ratio` of the power_max
    # min_production = min_load_ratio * power_max   # TODO - maybe it's a internal variable
    "Fuel curve intercept coefficient (L/(h × kW))"
    F0
    "Fuel curve slope (L/(h × kW))"
    F1

    # economics
    "Fuel cost (currency unit/L)"
    fuel_cost
    "Investiment cost (currency unit/kW)"
    investment_cost
    "Operation and maintenance cost (currency unit/(kW.h))"
    om_cost
    "Replacement cost (currency unit/kW)"
    replacement_cost
    "Salvage cost (currency unit/kW)"
    salvage_cost
    "Lifetime (h)"
    lifetime
end

"Photovoltaic parameters."
struct Photovoltaic <: NonDispatchables
    "Rated power (kW)"
    power_rated   # decision variable
    "Derating factor ∈ [0,1]"
    derating_factor
    "Incident global solar radiation (kW/m²)"
    IT
    "Standard amount of global solar radiation (kW/m²)"
    IS

    # economics
    "Investiment cost (currency unit/kW)"
    investment_cost
    "Operation and maintenance cost (currency unit/kW)"
    om_cost
    "Replacement cost (currency unit/kW)"
    replacement_cost
    "Salvage cost (currency unit/kW)"
    salvage_cost
    "Lifetime (years)"
    lifetime

    # Photovoltaic(fPV, IT, IS, Y_PV) = new(fPV, IT, IS, Y_PV)
end

"Wind turbine parameters."
struct WindPower <: NonDispatchables
    "Rated power (kW)"
    power_rated
    "Cut-in speed (m/s)"
    U_cut_in
    "Cut-out speed (m/s)"
    U_cut_out
    "Rated speed (m/s)"
    U_rated
    "Wind speed at the measurement height (m/s)"
    Uanem
    "Hub height (m)"
    zhub
    "Measurement height (m)"
    zanem
    "Roughness length (m)"
    z0
    # TODO rho
    # TODO rho0

    # economics
    "Investiment cost (currency unit/kW)"
    investment_cost
    "Operation and maintenance cost (currency unit/kW)"
    om_cost
    "Replacement cost (currency unit/kW)"
    replacement_cost
    "Salvage cost (currency unit/kW)"
    salvage_cost
    "Lifetime (years)"
    lifetime
end

"Battery parameters."
struct Battery
    "Initial energy (kWh)"
    energy_initial
    "Rated energy capacity (kWh)"
    energy_max    # Eb_max
    "Minimum energy level (kWh)"
    energy_min    # Eb_min  TODO - it could be the minimum state of charge too
    "Maximum charge power ∈ ``\\mathbf{R}^-`` (kW)"
    power_min     # Pb_min - charge (negative)
    "Maximum discharge power (kW)"
    power_max     # Pb_max - discharge
    "Linear loss factor ∈ [0,1]"
    loss

    # economics
    "Investiment cost (currency unit/kWh)"
    investment_cost
    "Operation and maintenance cost (currency unit/kWh)"
    om_cost
    "Replacement cost (currency unit/kWh)"
    replacement_cost
    "Salvage cost (currency unit/kWh)"
    salvage_cost
    "Lifetime (years)"
    lifetime
    "Maximum number of cycles"
    lifetime_throughput  # max throughput
end

# Operation variables - Trajectory
struct OperVarsTraj
    # load
    #= "Net load at each time instant after using the renewables power (kW)"
    Pnl_req =#
    "Net load at each time instant after dispatch (kW)"
    power_net_load
    "Unmet load/Load shedding power at each time instant (kW)"
    power_shedding
    # diesel generator
    "Diesel generator power at each time instant (kW)"
    Pgen
    # battery
    "Battery energy at each time instant (kWh)"
    Ebatt
    "Battery power at each time instant (kW)"
    Pbatt
    "Maximum battery discharge power at time t (kW)"
    Pbatt_dmax
    "Maximum battery charge power at time t (kW)"
    Pbatt_cmax
    # renewables sources
    "Renewables curtailment power at each time instant (kW)"
    power_curtailment
end

# Operation variables - Aggregation
struct OperVarsAggr
    # load
    "Load energy served in one year (kWh)"
    energy_served
    "Maximum load shedding power (kW)"
    power_shedding_max
    "Maximum consecutive duration of load shedding (h)"
    shedding_duration_max
    "Load shedding energy in one year (kWh)"
    energy_shedding_total
    "Ratio between load shedding energy and total consumption in one year (%)"
    shedding_rate
    # diesel generator
    "Number of diesel generator operation hours in one year (h)"
    DG_operation_hours
    "Diesel consumption in one year (L)"
    fuel_consumption
    # battery
    "Number of completed battery cycles in one year"
    annual_throughput
    # renewables sources
    "Maximum renewables curtailment power (kW)"
    power_curtailment_max
    "Ratio between energy supplied by renewables and energy served in one year (%)"
    renewables_rate
end

# Total costs
struct TotalCosts
    # general
    "Levelized cost of electricity (currency unit)"
    lcoe
    "Cost of electricity (currency unit)"
    coe # annualized
    "Net present cost (currency unit)"
    npc
    "Present investment cost (currency unit)"
    total_investment_cost
    "Present replacement cost (currency unit)"
    total_replacement_cost
    "Present operation and maintenance cost (currency unit)"
    total_om_cost
    "Present salvage cost (currency unit)"
    total_salvage_cost

    # components
    "Generator's total present cost (currency unit)"
    DG_total_cost
    "Generator's present investment cost (currency unit)"
    DG_investment_cost
    "Generator's present replacement cost (currency unit)"
    DG_replacement_cost
    "Generator's present operation and maintenance cost (currency unit)"
    DG_om_cost
    "Generator's present salvage cost (currency unit)"
    DG_salvage_cost
    "Generator's present fuel cost (currency unit)"
    DG_fuel_cost

    "Battery's total present cost (currency unit)"
    BT_total_cost
    "Battery's present investment cost (currency unit)"
    BT_investment_cost
    "Battery's present replacement cost (currency unit)"
    BT_replacement_cost
    "Battery's present operation and maintenance cost (currency unit)"
    BT_om_cost
    "Battery's present salvage cost (currency unit)"
    BT_salvage_cost

    "Photovoltaic's total present cost (currency unit)"
    PV_total_cost
    "Photovoltaic's present investment cost (currency unit)"
    PV_investment_cost
    "Photovoltaic's present replacement cost (currency unit)"
    PV_replacement_cost
    "Photovoltaic's present operation and maintenance cost (currency unit)"
    PV_om_cost
    "Photovoltaic's present salvage cost (currency unit)"
    PV_salvage_cost

    "Wind turbine's total present cost (currency unit)"
    WT_total_cost
    "Wind turbine's present investment cost (currency unit)"
    WT_investment_cost
    "Wind turbine's present replacement cost (currency unit)"
    WT_replacement_cost
    "Wind turbine's present operation and maintenance cost (currency unit)"
    WT_om_cost
    "Wind turbine's present salvage cost (currency unit)"
    WT_salvage_cost
end

# Microgrid
struct Microgrid
    project::Project
    power_load
    dieselgenerator::DieselGenerator
    # photovoltaic::Photovoltaic
    # windpower::WindPower
    battery::Battery
    nondispatchables::Vector{NonDispatchables}
end


text:
module Finometrics
using Revise

#=if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end=#

#update steps
#cd [current project path]
#] activate .
#] update
#] resolve
#] gc
#] build
#] precompile

##################Dependencies

using  DataFrames, Distributions, StatsBase, GLM, CategoricalArrays,  Dates, LinearAlgebra, StatsModels, Formatting, NLopt

try
  using CUDA
  CUDA.allowscalar(false) #default
catch
  println("Note: CuArrays not installed")
  CUDA = Nothing #allows checks for CuArrays to work (should all return false)
end

#=function testfunc()
  local f::FormulaTerm

  f = @formula(x ~ y)

  println(f)
end

testfunc()=#

import Base: +, -, ==, >, <, ≥, ≤, isless, isequal
import Base: push!, length, iterate, broadcastable

#######################MACROS###################
#useful macro for conditionally running things in parallel
#BROKEN
#=macro mpar(cond, expr)
  quote
    if $(esc(cond))
        :($(Threads.@threads($expr)))
    else
        :($($expr))
    end
  end
end=#

######################Methods####################
export FMLM, #Regression methods
  project!, #regression methods (won't work)
  resid!,
  coeff!,
  homoskedasticΣ!,
  modifiedwhiteΣ!,
  whiteΣ!,
  neweywestΣ!,
  homoskedasticΣslow,
  modifiedwhiteΣslow,
  whiteΣslow,
  neweywestΣslow,
  FMExpr,
  FMSym,
  FMQR,
  βterm, #odd naming is due to conflicts with other packages
  R,
  R²,
  neweywestΣfunc,
  neweywestΣfuncslow,
  neweywestpanelΣ!,
  neweywestpanelΣslow!,
  neweywestpanelΣfunc,
  neweywestpanelΣslowfunc,
  clusteredΣ!,

  textable, #IO Mthods
  writetextable,
  array2string,
  num2str,

  skewnessStat, #FMstat methods
  kurtosisStat,
  JBStat,
  UnivariateTest,
  autoρ,
  ljungBox,
  test,
  pAutoρ,
  ARMA,
  ARMAFillYHat!,
  ARMALogLike,
  ARMAEstimate,
  generalizeARMA!,
  narrowARMA!,
  winsorizequantile,
  winsorizelevel,

  lagwithin!, #FMLagWithin
  differencewithin!,
  differencewithin2!,
  lagwithin2!,
  lagwithin2sorted!,
  lagwithin2sorted,
  leadwithin2!,

  FMΦInv, #FMNumerical
  FMΦ,
  inSorted,
  vbind,

  VAR, #FMVAR
  lagDF!,
  VAREstimate,
  VARPrediction,
  varΣ,
  propagateImpulse,

  FMSpecs, #FMSpecs
  computeFMLMresults!,

  NBool, #Exported types
  NDict,
  NString,
  NSymbol,
  NType,
  NInt,
  NFloat64,
  NDate,

  MString,
  MBool,
  MFloat64,
  MInt,
  MSymbol,
  MDate,
  ∞,

  CU32,
  CU64,

  DField,
  MDField,
  NDField,

  #yearmonth types
  YearQuarter,
  MYearQuarter,
  YearMonth,
  MYearMonth,
  eom,
  bom,
  eoq,
  boq,

  ##WARNING DEPRECIATED methods #WARNING
  texTable,
  writeTexTable,
  array2String,
  num2Str,
  vec2String,
  getCoeff!,
  getHomoskedΣ!,
  getModWhiteΣ!,
  getWhiteΣ!,
  getNeweyWest!,
  getModWhiteΣSlow,
  getHomoskedΣSlow,
  getWhiteΣSlow,
  getNeweyWestSlow,
  getTerm,
  getR,
  getR²,
  getNeweyWestFunc,
  getClustered!,
  winsorize!

##################Custom types

###FMReg Types
const FMExpr = Union{Symbol,Expr,Nothing}
const FMData = Union{Float64,Int,Date,Symbol, CategoricalValue, Missing, Nothing}#,
abstract type FMModel end

###Convenience Types for Export
const NBool = Union{Nothing, Bool}
const NSymbol = Union{Nothing, Symbol}
const NType = Union{Nothing, Type}
const NString = Union{Nothing, String}
const NInt = Union{Int, Nothing}
const NFloat64 = Union{Float64, Nothing}
const NDate = Union{Date, Nothing}

const MBool = Union{Bool, Missing}
const MInt = Union{Int, Missing}
const MFloat64 = Union{Float64, Missing}
const MSymbol = Union{Symbol, Missing}
const MString = Union{String, Missing}
const MDate = Union{Date, Missing}
const ∞ = Inf

const DField = Union{String, Symbol}
const MDField = Union{String, Symbol, Missing}
const NDField = Union{String, Symbol, Nothing}

###VAR Types
const Structure = Union{Nothing, Vector{Float64}}
const CoefficientErrors = Union{Nothing, Vector{Float64}}
const RegressionSS = Union{Nothing, Float64}
const RegressionN = Union{Nothing, Int}

################Constants
const spinlock = Threads.SpinLock()

#shared cosntants
const DEFAULT_STAR_LEGEND = "*p<0.1, **p<0.05, ***p<0.01"
const DEBUG_FMMOD = false
const DEFAULT_DECIMALS = 4

#Stat constants
const MAX_TIME = 30.
const F_TOL_ABS = 10. ^ -12.
const F_TOL_REL = 10. ^ -10.
const OPT_BOUND = 2.

#numerical constants
const sqrt2Inv = 1.0/(2.0^.5)
const sqrt2PiInv = 1.0/(2.0*π)^0.5
#  MFloat64, MInt, MBool, MSymbol}

const UPPER_PAD = "%This table was programatically generated using the TexTables package and my own code.\n"
const LOWER_PAD = "\n"

###################Files############

# include("$(pwd())/FMReg.jl")
include("FMReg.jl")
include("FMCov.jl")
include("FMIO.jl")
include("FMStat.jl")
include("FMLagWithin.jl")
include("FMVAR.jl")
include("FMNumerical.jl")
include("FMSpecs.jl")
include("FMYearQuarter.jl")
include("FMTest.jl") #comment this out when testing
include("FMRegBroken2SLS.jl")

#type is derived from a user-defined type
const MYearQuarter = Union{YearQuarter, Missing}
const MYearMonth = Union{YearMonth, Missing}

end

max_stars_repo_name: clintonTE/Finometrics

max_issues_repo_name: clintonTE/Finometrics

max_forks_repo_name: clintonTE/Finometrics


text:
module FED

using Printf: @printf, @sprintf
using SparseArrays: spzeros
using QuantumLattices: AbstractFockOperator, Operators, FOperator, CREATION, rank
using QuantumLattices.Prerequisites: Float
using QuantumLattices.Mathematics.VectorSpaces: VectorSpace

import SparseArrays: SparseMatrixCSC
import QuantumLattices: kind, id, dimension
import QuantumLattices.Mathematics.AlgebraOverFields: idtype

export FockBasis, GBasis, PBasis, SBasis, SPBasis

"""
    FockBasis{I<:Unsigned} <: VectorSpace{I}

Abatract type for all kinds of Fock basis.
"""
abstract type FockBasis{I<:Unsigned} <: VectorSpace{I} end
@inline kind(bs::FockBasis) = kind(typeof(bs))
@inline idtype(bs::FockBasis) = idtype(typeof(bs))
function Base.show(io::IO, bs::FockBasis)
    @printf io "%s:\n" repr(bs)
    for i = 1:dimension(bs)
        @printf io "  %s\n" string(bs[i], base=2)
    end
end
@inline Base.repr(bs::FockBasis) = @sprintf "%s(%s)" nameof(typeof(bs)) join(id(bs), ",")
@inline Base.issorted(bs::FockBasis) = true

"""
    GBasis{I}(nstate::Int) where I<:Unsigned

Generic fermionic/hard-core-bosonic basis, which uses no conserved quantities.
"""
struct GBasis{I<:Unsigned} <: FockBasis{I}
    nstate::Int
end
@inline kind(::Type{<:GBasis}) = :g
@inline idtype(::Type{<:GBasis}) = Int
@inline id(bs::GBasis) = bs.nstate
@inline dimension(bs::GBasis) = 2^bs.nstate
@inline Base.getindex(bs::GBasis, i::Int) = i-1
@inline Base.searchsortedfirst(bs::GBasis, i::Int) = i+1

"""
"""
struct PBasis{I<:Unsigned} <: FockBasis{I}
    nstate::Int
    nparticle::Int
    table::Vector{I}
end
@inline kind(::Type{<:PBasis}) = :p
@inline idtype(::Type{<:PBasis}) = Tuple{Int, Int}
@inline id(bs::PBasis) = (bs.nstate, bs.nparticle)
@inline dimension(bs::PBasis) = length(bs.table)
@inline Base.getindex(bs::PBasis, i::Int) = bs.table[i]

"""
"""
struct SBasis{I<:Unsigned} <: FockBasis{I}
    nstate::Int
    spinz::Rational{Int}
    table::Vector{I}
end
@inline kind(::Type{<:SBasis}) = :s
@inline idtype(::Type{<:SBasis}) = Tuple{Int,Rational{Int}}
@inline id(bs::SBasis) = (bs.nstate, bs.spinz)
@inline dimension(bs::SBasis) = length(bs.table)
@inline Base.getindex(bs::SBasis, i::Int) = bs.table[i]

"""
"""
struct SPBasis{I<:Unsigned} <: FockBasis{I}
    nstate::Int
    nparticle::Int
    spinz::Rational{Int}
    table::Vector{I}
end
@inline kind(::Type{<:SPBasis}) = :sp
@inline idtype(::Type{<:SPBasis}) = Tuple{Int, Int, Rational{Int}}
@inline id(bs::SPBasis) = (bs.nstate, bs.nparticle, bs.spinz)
@inline dimension(bs::SPBasis) = length(bs.table)
@inline Base.getindex(bs::SPBasis, i::Int) = bs.table[i]

"""
    SparseMatrixCSC(operator::AbstractFockOperator, basis::FockBasis, table) -> SparseMatrixCSC{valtype(operator), Int}
    SparseMatrixCSC(operator::AbstractFockOperator, basis₁::FockBasis, basis₂::FockBasis, table) -> SparseMatrixCSC{valtype(operator), Int}
    SparseMatrixCSC(opts::Operators, basis::FockBasis, table, dtype::Type{<:Number}=Complex{Float}) -> SparseMatrixCSC{dtype, Int}
    SparseMatrixCSC(opts::Operators, basis₁::FockBasis, basis₂::FockBasis, table, dtype::Type{<:Number}=Complex{Float}) -> SparseMatrixCSC{dtype, Int}

Get the CSC-formed sparse matrix representation of
1) an even-ranked operator,
2) an odd-ranked operator,
3) a set of even-ranked operators,
4) a set of odd-ranked operators.
"""
function SparseMatrixCSC(operator::AbstractFockOperator, basis::FockBasis, table)
    @assert rank(operator)%2==0 "SparseMatrixCSC error: two bases are needed for odd ranked Fock operators."
    ndata, data, indices, indptr = 1, zeros(valtype(operator), dimension(basis)), zeros(Int, dimension(basis)), zeros(Int, dimension(basis)+1)
    eye, temp, flag = one(basis|>eltype), zeros(basis|>eltype, rank(operator)+1), true
    seqs = NTuple{rank(operator), Int}(table[id(operator)[i]]-1 for i = 1:rank(operator))
    nambus = NTuple{rank(operator), Bool}(oid.index.nambu==CREATION for oid in id(operator))
    for i = 1:dimension(basis)
        flag = true
        indptr[i] = ndata
        temp[1] = basis[i]
        for j = 1:rank(operator)
            (temp[j]&eye<<seqs[j]>0)==nambus[j] && (flag = false; break)
            temp[j+1] = nambus[j] ? temp[j]|eye<<seqs[j] : temp[j]&~(eye<<seqs[j])
        end
        if flag
            nsign = 0
            isa(operator, FOperator) && for j = 1:rank(operator)
                for k = 0:seqs[j]-1
                    temp[j]&eye<<k>0 && (nsign += 1)
                end
            end
            indices[ndata] = searchsortedfirst(basis, temp[end])
            data[ndata] = operator.value*(-1)^nsign
            ndata += 1
        end
    end
    indptr[end] = ndata+1
    return SparseMatrixCSC(dimension(basis), dimension(basis), indptr, indices[1:ndata], data[1:ndata])
end
function SparseMatrixCSC(operator::AbstractFockOperator, basis₁::FockBasis, basis₂::FockBasis, table)
    @assert rank(operator)%2==1 "SparseMatrixCSC error: only one basis is needed for even ranked Fock operators."
    ndata, data, indices, indptr = 1, zeros(valtype(operator), dimension(basis₁)), zeros(Int, dimension(basis₁)), zeros(Int, dimension(basis₁)+1)
    eye, temp, flag = one(basis₁|>eltype), zeros(basis₁|>eltype, rank(operator)+1), true
    seqs = NTuple{rank(operator), Int}(table[id(operator)[i]]-1 for i = 1:rank(operator))
    nambus = NTuple{rank(operator), Bool}(oid.index.nambu==CREATION for oid in id(operator))
    for i = 1:dimension(basis₁)
        flag = true
        indptr[i] = ndata
        temp[1] = basis₁[i]
        for j = 1:rank(operator)
            (temp[j]&eye<<seqs[j]>0)==nambus[j] && (flag = false; break)
            temp[j+1] = nambus[j] ? temp[j]|eye<<seqs[j] : temp[j]&~(eye<<seqs[j])
        end
        if flag
            nsign = 0
            isa(operator,FOperator) && for j = 1:rank(operator)
                for k = 0:seqs[j]-1
                    temp[j]&eye<<k>0 && (nsign += 1)
                end
            end
            indices[ndata] = searchsortedfirst(basis₂, temp[end])
            data[ndata] = operator.value*(-1)^nsign
            ndata += 1
        end
    end
    indptr[end] = ndata+1
    return SparseMatrixCSC(dimension(basis₂), dimension(basis₁), indptr, indices[1:ndata], data[1:ndata])
end
function SparseMatrixCSC(opts::Operators, basis::FockBasis, table, dtype::Type{<:Number}=Complex{Float})
    result = spzeros(dtype, dimension(basis), dimension(basis))
    for opt in values(opts)
        result += SparseMatrixCSC(opt, basis, table)
    end
    return result
end
function SparseMatrixCSC(opts::Operators, basis₁::FockBasis, basis₂::FockBasis, table, dtype::Type{<:Number}=Complex{Float})
    result = spzeros(dtype, dimension(basis₂), dimension(basis₁))
    for opt in values(opts)
        result += SparseMatrixCSC(opt, basis₁, basis₂, table)
    end
    return result
end

end  # module


text:
T=9
tsname="joGaussian"
@testset "$tsname" begin
for t=2:T # start test loop
    m=4*t
    n=m+rand(-2:2:2)
    A=joGaussian(m,n)
    An=joGaussian(m,n,normalized=true)
    Ai=joGaussian(m,n,implicit=true)
    Ain=joGaussian(m,n,implicit=true,normalized=true)
    Ao=joGaussian(m,orthonormal=true)
    Aor=joGaussian(m,m+2,orthonormal=true)

    verbose && println("$tsname ($m)")
    @testset "$m x $m" begin
        @test isadjoint(A)[1]
        @test isadjoint(An)[1]
        @test isadjoint(Ai)[1]
        @test isadjoint(Ain)[1]
        @test isadjoint(Ao)[1]
        @test isadjoint(Aor)[1]
        @test islinear(A)[1]
        @test islinear(An)[1]
        @test islinear(Ai)[1]
        @test islinear(Ain)[1]
        @test islinear(Ao)[1]
        @test islinear(Aor)[1]
    end

end # end test loop
end


text: # c'est un commentaire # et ça aussi

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using Documenter, BioXP

makedocs(
    modules = [BioXP],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Harrison Smith",
    sitename = "BioXP.jl",
    pages = Any["index.md"]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/hbsmith/BioXP.jl.git",
)

max_stars_repo_name: hbsmith/BioXP

max_issues_repo_name: hbsmith/BioXP

max_forks_repo_name: hbsmith/BioXP


text:
using CompScienceMeshes
using Test

for T in [Float32, Float64]
    faces = readmesh(joinpath(@__DIR__, "assets/rectangle32.in"), T=T)
    edges = skeleton(faces,1)
    verts = skeleton(faces,0)
    bnd = boundary(faces)

    int_edges = submesh(interior_tpredicate(faces), edges)
    # pred = interior_tpredicate(int_edges)

    bnd_verts = skeleton(bnd,0)
    function pred(node)
        !(node in cells(bnd_verts))
    end
    int_verts = submesh(pred, verts)

    @test numcells(int_verts) == 9
    @test numcells(edges) - numcells(int_edges) == 16
end

max_stars_repo_name: krcools/CompScienceMeshes

max_issues_repo_name: krcools/CompScienceMeshes

max_forks_repo_name: krcools/CompScienceMeshes


text:
using GMM
using ModelsGenerators
using Ipopt
using NLopt

srand(1)
y, x, z = ModelsGenerators.sim_iv_d01(CP = 100)

g(theta) = z.*(y-x*theta)

W = inv(z'z);

W = eye(5);
out = gmm(g, [.1], W)

W = pinv(GMM.mfvcov(out, HC0())*100)
out = GMM.gmm(g, [.0], W)

coef(out)
sqrt(vcov(out))

sqrt(vcov(out, ParzenKernel(3))

W = eye(5);
out = GMM.gmm(g, [.1], W, solver = NLoptSolver(algorithm = :LD_LBFGS))

W = pinv(GMM.mfvcov(out, HC1()))
out = GMM.gmm(g, [.1], W, solver = NLoptSolver(algorithm = :LD_LBFGS))


text:
using Pkg.TOML
using Pkg.Artifacts
using RelocatableFolders

# Fix the NPM version for now
# Now that we're not distributing Jupyter stuff via NPM, I'd rather just check
# in compiled assets (maybe using Artifacts if that's relatively easy, otherwise
# just as files in git).
const WEBIO_VERSION = "0.8.15"

# This needs to be cleaned up, it was rushed to get both the build step and relocatable paths to work
path = @path normpath(joinpath(@__DIR__, "..", "packages"))
isdir(path) || mkpath(path)
const PACKAGES_PATH = @path path
path = normpath(joinpath(@__DIR__, "bundles"))
isdir(path) || mkpath(path)
const BUNDLES_PATH = @path path

function bundleurl(pkg::String, filename::String)
    return "https://unpkg.com/@webio/$(pkg)@$(WEBIO_VERSION)/dist/$(filename)"
end

path = joinpath(BUNDLES_PATH, "webio.bundle.js")
isfile(path) || mkpath(path)
const CORE_BUNDLE_PATH = path
const CORE_BUNDLE_URL = bundleurl("webio", "webio.bundle.js")

path = joinpath(BUNDLES_PATH, "generic-http.bundle.js")
isfile(path) || mkpath(path)
const GENERIC_HTTP_BUNDLE_PATH = @path path
const GENERIC_HTTP_BUNDLE_URL = bundleurl("generic-http-provider", "generic-http.bundle.js")

path = joinpath(BUNDLES_PATH, "mux.bundle.js")
isfile(path) || mkpath(path)
const MUX_BUNDLE_PATH = @path path
const MUX_BUNDLE_URL = bundleurl("mux-provider", "mux.bundle.js")

# Deprecated! Remove for WebIO version 1.0.0
path = joinpath(BUNDLES_PATH, "blink.bundle.js")
isfile(path) || mkpath(path)
const BLINK_BUNDLE_PATH = @path path
const BLINK_BUNDLE_URL = bundleurl("blink-provider", "blink.bundle.js")

function download_bundle(name::String, path::String, url::String)
    if !isfile(path)
        @info "Downloading WebIO $(name) bundle from unpkg..."
        download(url, path)
    end
end

# TODO: this is all an ugly hack to avoid trying to build JS when other packages (that use WebIO)
# are just trying to run their own tests. It desperately needs to be restructured.
function isci()
    return (
        get(ENV, "CI", "false") == "true" &&
        split(get(ENV, "GITHUB_REPOSITORY", "Foo/Bar.jl"), '/')[2] == "WebIO.jl"
    )
end

function download_js_bundles()
    if isci()
        # In CI, we always build the bundles from scratch.
        @info "Not downloading WebIO bundles in CI; building instead..."
        include("./_bundlejs.jl")
        return
    end

    mkpath(BUNDLES_PATH)
    bundle_artifact_path = artifact"web"
    for asset in readdir(bundle_artifact_path)
        @debug("Copying", bundle_artifact_path, asset, BUNDLES_PATH)
        cp(joinpath(bundle_artifact_path, asset), joinpath(BUNDLES_PATH, asset); force = true)
    end
    # These commands are probably still useful if you want to create the artifact.
    # download_bundle("core", CORE_BUNDLE_PATH, CORE_BUNDLE_URL)
    # download_bundle("generic-http", GENERIC_HTTP_BUNDLE_PATH, GENERIC_HTTP_BUNDLE_URL)
    # download_bundle("mux", MUX_BUNDLE_PATH, MUX_BUNDLE_URL)
    # download_bundle("blink", BLINK_BUNDLE_PATH, BLINK_BUNDLE_URL)
end


text:
function run_model!(; model::Agents.AgentBasedModel, steps::Int64)
    agent_df, model_df = Agents.run!(
        model, agent_step!, model_step!, steps,
        adata=[:pos, :socialnorm, :conscientiousness, :littering],
        mdata=[:new_trash, :trash_aggregated]
    )
    return agent_df, model_df
end

function run_experiment(; config::Config, replicates::Int64, rnd_seed::Int64)
    Random.seed!(rnd_seed)
    agent_df_list = DataFrames.DataFrame[]
    model_df_list = DataFrames.DataFrame[]
    networks = Dict()
    for rep in 1:replicates
        space = create_space(config)
        model = initialize_model(config, space)
        networks[rep] = deepcopy(model.space.graph)
        agent_df, model_df = run_model!(model=model, steps=config.n_iter)
        agent_df[!, :replicate] .= deepcopy(rep)
        model_df[!, :replicate] .= deepcopy(rep)
        push!(agent_df_list, deepcopy(agent_df))
        push!(model_df_list, deepcopy(model_df))
    end
    agent_df = reduce(vcat, agent_df_list)
    model_df = reduce(vcat, model_df_list)
    return networks, agent_df, model_df
end

function run_batch(; configs::Dict, replicates::Int64, batchname::String, compress::Bool)
    make_folders(batchname)
    path = joinpath("experiments", batchname)
    progress_bar = ProgressMeter.Progress(length(keys(configs)) * 2, 1, "Running simulations...")
    for config_key in keys(configs)
        networks, agent_df, model_df = run_experiment(config=configs[config_key], replicates=replicates, rnd_seed=0)
        adata_filename = "config_" * lpad(string(config_key), 2, "0") * "_adata.feather"
        mdata_filename = "config_" * lpad(string(config_key), 2, "0") * "_mdata.feather"
        ProgressMeter.next!(progress_bar)
        names_to_camelcase!(agent_df)
        names_to_camelcase!(model_df)
        write_dataframe(df=agent_df, filename=adata_filename, path=path)
        write_dataframe(df=model_df, filename=mdata_filename, path=path)
        write_networks(networks, path=path, config_key=config_key)
        ProgressMeter.next!(progress_bar)
    end
    println("Simulations done.")
    configs_to_dataframe(configs, path=path, write=true)
    if compress
        patterns = ["_graph.txt", "_mdata.feather", "_adata.feather"]
        archive_filenames = ["graphs.7z", "mdata.7z", "adata.7z"]
        for (pattern, archive_filename) in zip(patterns, archive_filenames)
            compress_data(path=path, pattern=pattern, archive_filename=archive_filename);
        end
    end
    println("Batch $batchname run successfully.")
    return "successful"
end


text:
@testsuite "random" AT->begin
    rng = if AT <: AbstractGPUArray
        GPUArrays.default_rng(AT)
    else
        Random.default_rng()
    end

    @testset "rand" begin  # uniform
        for T in (Int16, Int32, Int64,
                  Float16, Float32, Float64,
                  Complex{Float16}, Complex{Float32}, Complex{Float64},
                  Complex{Int32}, Complex{Int64}), d in (10, (10,10))
            A = AT{T}(undef, d)
            B = copy(A)
            rand!(rng, A)
            rand!(rng, B)
            @test Array(A) != Array(B)

            Random.seed!(rng)
            Random.seed!(rng, 1)
            rand!(rng, A)
            Random.seed!(rng, 1)
            rand!(rng, B)
            @test all(Array(A) .== Array(B))
        end

        A = AT{Bool}(undef, 1024)
        fill!(A, false)
        rand!(rng, A)
        @test true in Array(A)
        fill!(A, true)
        rand!(rng, A)
        @test false in Array(A)
    end

    @testset "randn" begin  # normally-distributed
        for T in (Float16, Float32, Float64), d in (2, (2,2))
            A = AT{T}(undef, d)
            B = copy(A)
            randn!(rng, A)
            randn!(rng, B)
            @test !any(Array(A) .== Array(B))

            Random.seed!(rng)
            Random.seed!(rng, 1)
            randn!(rng, A)
            Random.seed!(rng, 1)
            randn!(rng, B)
            @test all(Array(A) .== Array(B))
        end
    end
end


text:
using Test

using Bridge, StaticArrays, Distributions
using Statistics, Random, LinearAlgebra
POSSIBLE_PARAMS = [:regular, :simpleAlter, :complexAlter, :simpleConjug,
                   :complexConjug]
SRC_DIR = joinpath(Base.source_dir(), "..", "src")

include("test_ODE_solver_change_pt.jl")
include("test_blocking.jl")


text:
"""
    struct SymSparseMatrixCSR{T,Ti<:Integer} <: AbstractSparseMatrix{T,Ti}

Matrix type for storing symmetric sparse matrices in the
Compressed Sparse Row format. The standard way
of constructing SparseMatrixCSR is through the
[`symsparsecsr`](@ref) function.
"""
struct SymSparseMatrixCSR{Bi,T,Ti<:Integer} <: AbstractSparseMatrix{T,Ti}
    uppertrian :: SparseMatrixCSR{Bi,T,Ti}
end

function sparse_from_coo(M::Type{<:SymSparseMatrixCSR}, I,J,V,m,n)
  symsparsecsr(M, I,J,V,m,n)
end

# SparseMatrix interface implementation
@inline function push_coo!(::Type{<:SymSparseMatrixCSR},
        I::Vector,J::Vector,V::Vector,ik::Integer,jk::Integer,vk::Number)
    (ik>jk) && return nothing

    push!(I, ik)
    push!(J, jk)
    push!(V, vk)
    nothing
end

@inline function is_entry_stored(::Type{<:SymSparseMatrixCSR},i::Integer,j::Integer)
  (i<=j)
end

function finalize_coo!(T::Type{<:SymSparseMatrixCSR},
        I::Vector,J::Vector,V::Vector, m::Integer, n::Integer)
    m == n || throw(DimensionMismatch("matrix is not square: dimensions are ($m, $n)"))
    touched = zeros(Bool,m)
    for k in 1:length(I)
        Ik = I[k]
        Jk = J[k]
        if Ik == Jk
            touched[Ik] = true
        end
    end
    for k in 1:m
        if ! touched[k]
            push_coo!(T,I,J,V,k,k,zero(eltype(V)))
        end
    end
end

add_entry!(A::SymSparseMatrixCSR{Bi,Tv,Ti},v::Number,i::Integer,j::Integer,combine::Function=+) where {Bi,Tv,Ti<:Integer} =
        return i>j ? A : add_entry!(A.uppertrian,v,i,j,combine)

function copy_entries!(a::SymSparseMatrixCSR,b::SymSparseMatrixCSR)
  _copy_entries_sparse!(a,b)
end

# CompressedSparseMatrix interface implementation

hasrowmajororder(::Type{<:SymSparseMatrixCSR}) = true
hasrowmajororder(a::SymSparseMatrixCSR) = hasrowmajororder(SymSparseMatrixCSR)
hascolmajororder(::Type{<:SymSparseMatrixCSR}) = false
hascolmajororder(a::SymSparseMatrixCSR) = hascolmajororder(SymSparseMatrixCSR)
getptr(S::SymSparseMatrixCSR) = getptr(S.uppertrian)
getindices(S::SymSparseMatrixCSR) = colvals(S)
getindex(A::SymSparseMatrixCSR, x::Integer, y::Integer) = getindex(A.uppertrian,min(x,y),max(x,y))
nnz(S::SymSparseMatrixCSR) = nnz(S.uppertrian)
count(pred, S::SymSparseMatrixCSR) = count(pred, S.uppertrian)
nonzeros(S::SymSparseMatrixCSR) = nonzeros(S.uppertrian)
nzrange(S::SymSparseMatrixCSR, row::Integer) = nzrange(S.uppertrian, row)
findnz(S::SymSparseMatrixCSR) = findnz(S.uppertrian)
rowvals(S::SymSparseMatrixCSR) = rowvals(S.uppertrian)
colvals(S::SymSparseMatrixCSR) = colvals(S.uppertrian)

show(io::IO, S::SymSparseMatrixCSR) = show(io, S.uppertrian)

size(A::SymSparseMatrixCSR) = size(A.uppertrian)

*(A::SymSparseMatrixCSR, v::Vector) = (y = similar(v,A.uppertrian.n);mul!(y,A,v))

function symsparsecsr(T::Type{<:SymSparseMatrixCSR{Bi,Tv,Ti}},I::Vector{Ti},J::Vector{Ti},V::Vector{Tv},args...)  where {Bi,Tv,Ti}
    m = length(args)>0 ? args[1] : isempty(I) ? 0 : Int(maximum(I))
    n = length(args)>1 ? args[2] : m
    c = length(args)>2 ? args[3] : +
    m == n || throw(DimensionMismatch("matrix is not square: dimensions are ($m, $n)"))
    SymSparseMatrixCSR(sparsecsr(SparseMatrixCSR{Bi},I,J,V,m,n,c))
end

"""
    symsparsecsr(args...)

Same args than `sparse`.
"""
symsparsecsr(I,J,V,args...) =
           SymSparseMatrixCSR(sparsecsr(I,J,V,args...))
symsparsecsr(T::Type{<:SymSparseMatrixCSR},I::Vector{Ti},J::Vector{Ti},V::Vector{Tv},args...) where {Tv,Ti} =
           symsparsecsr(SparseMatrixCSR{1,Tv,Ti},I,J,V,args...)
symsparsecsr(T::Type{<:SymSparseMatrixCSR{Bi}},I::Vector{Ti},J::Vector{Ti},V::Vector{Tv},args...) where {Bi,Tv,Ti} =
           symsparsecsr(SparseMatrixCSR{Bi,Tv,Ti},I,J,V,args...)

function show(io::IO, ::MIME"text/plain", S::SymSparseMatrixCSR)
    xnnz = nnz(S)
    print(io, S.uppertrian.m, "×", S.uppertrian.n, " ",
              typeof(S), " with ", xnnz, " stored ",
              xnnz == 1 ? "entry" : "entries")
    if xnnz != 0
        print(io, ":")
        show(IOContext(io, :typeinfo => eltype(S)), S)
    end
end

function mul!(y::AbstractVector,A::SymSparseMatrixCSR,v::AbstractVector{T}) where {T}
    A.uppertrian.n == size(v, 1) || throw(DimensionMismatch())
    A.uppertrian.m == size(y, 1) || throw(DimensionMismatch())

    y .= zero(T)
    for row = 1:size(y, 1)
        @inbounds for nz in nzrange(A,row)
            col = A.uppertrian.colval[nz]-A.uppertrian.offset
            y[row] += A.uppertrian.nzval[nz]*v[col]
            row != col && (y[col] += A.uppertrian.nzval[nz]*v[row])
        end
    end
    return y
end

function convert(::Type{<:SymSparseMatrixCSR{Bi,Tv,Ti}}, x::SymSparseMatrixCSR{xBi,xTv,xTi}) where {Bi,Tv,Ti,xBi,xTv,xTi}
    if (Bi,Tv,Ti) == (xBi,xTv,xTi)
        return x
    else
        return SymSparseMatrixCSR(convert(SparseMatrixCSR{Bi,Tv,Ti}, x.uppertrian))
    end
end


text:
using EmpiricalCDFs

struct Pareto
    alpha::Float64
    x0::Float64
end

Pareto(alpha) = Pareto(alpha,1.0)
Pareto() = Pareto(1.0)

Base.rand(p::Pareto) = p.x0 * rand()^(-one(p.alpha)/p.alpha)

function Base.rand(p::Pareto, n::Int)
    a = Array(Float64, n)
    for i in 1:n
        a[i] = rand(p)
    end
    a
end

const odir = joinpath(Pkg.dir("EmpiricalCDFs"), "test")
const ofile = joinpath(odir, "paretocdf.bin")

function write_test_data()
    d = Pareto()
    cdf = EmpiricalCDF()
    append!(cdf, rand(d,10^5))
    sort!(cdf)
    write(ofile, cdf)
    nothing
end

import EmpiricalCDFs: mleKS, scanmle

# julia> read_test_data()
# EmpiricalCDFs.MLEKS{Float64}(2.0018144538762677,0.0010018144538762677,0.2488588894720618)

function read_test_data()
    cdf = readcdf(ofile)
    mleks = mleKS(cdf)
    println(mleks)
    scanmle(cdf,15)
end


text:
function nthhamming(n :: UInt64) # :: Tuple{UInt32, UInt32, UInt32}
    # take care of trivial cases too small for band size estimation to work...
    n < 1 && throw("nthhamming:  argument must be greater than zero!!!")
    n < 2 && return (0, 0, 0)
    n < 3 && return (1, 0, 0)

    # some constants...
    log2of2, log2of3, log2of5 = 1.0, log(2, 3), log(2, 5)
    fctr, crctn = 6.0 * log2of3 * log2of5, log(2, sqrt(30))
    log2est = (fctr * Float64(n))^(1.0 / 3.0) - crctn # log2 answer from WP formula
    log2hi = log2est + 1.0 / log2est; width = 2.0 / log2est # up to 2X higher/lower

    # some really really big constants representing the "roll-your-own" big logs...
    biglog2of2 = BigInt(1267650600228229401496703205376)
    biglog2of3 = BigInt(2009178665378409109047848542368)
    biglog2of5 = BigInt(2943393543170754072109742145491)

    # loop to find the count of regular numbers and band of possible candidates...
    count :: UInt64 = 0; band = Vector{Tuple{BigInt,Tuple{UInt32,UInt32,UInt32}}}()
    fiveslmt = UInt32(ceil(log2hi / log2of5)); fives :: UInt32 = 0
    while fives < fiveslmt
        log2p = log2hi - fives * log2of5
        threeslmt = UInt32(ceil(log2p / log2of3)); threes :: UInt32 = 0
        while threes < threeslmt
            log2q = log2p - threes * log2of3
            twos = UInt32(floor(log2q)); frac = log2q - twos; count += twos + 1
            if frac <= width
                biglog = biglog2of2 * twos + biglog2of3 * threes + biglog2of5 * fives
                push!(band, (biglog, (twos, threes, fives)))
            end
            threes += 1
        end
        fives += 1
    end

    # process the band found including checks for validity and range...
    n > count && throw("nthhamming:  band high estimate is too low!!!")
    ndx = count - n + 1
    ndx > length(band) && throw("nthhamming:  band width estimate is too narrow!!!")
    sort!(band, by=(tpl -> let (lg,_) = tpl; -lg end)) # sort in decending order

    # get and return the answer...
    _, tri = band[ndx]
    tri
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module SimpleMatch
export @match
using Compat

"""
  matches arguments against locally created function and returns match

Simple match implementation which uses a locally created function.
  This gives a very familiar syntax and is natural to use.

Example
```
a = 1
b = 2.0
matched = @match(a, b) do f
  f(::Int, ::AbstractFloat) = b + 3
  f(_...) = a + 9
end
# 5.0
```

We need to define this as macro, as it turns out to be relatively difficult
  to define a function within an anonymous do notation
"""
macro match(define_func!, args...)
  @assert define_func!.head == :(->)  "expecting anonymous function as first argument"
  @assert define_func!.args[1].head == :tuple  "expecting anonymous function as first argument"
  @assert length(define_func!.args[1].args) == 1  "expecting define_func! to take exactly 1 argument, found $(length(define_func!.args[1].args))"

  function_name = define_func!.args[1].args[1]  # single argument of function
  function_definition = define_func!.args[2]  # function body

  esc(quote
    let
      function $function_name end
      $function_definition
      $function_name($(args...))
    end
  end)
end

end # module


text:
"""
# Usage

    solve(FDDE::FDDEProblem, h, DelayPI())

Use explicit rectangular product integration algorithm to solve an FDDE problem.

### References

```tex
@article{2020,
   title={On initial conditions for fractional delay differential equations},
   ISSN={1007-5704},
   url={http://dx.doi.org/10.1016/j.cnsns.2020.105359},
   DOI={10.1016/j.cnsns.2020.105359},
   journal={Communications in Nonlinear Science and Numerical Simulation},
   author={Garrappa, Roberto and Kaslik, Eva},
   year={2020},
}
```
"""
struct DelayPI <: FractionalDiffEqAlgorithm end

function solve(FDDE::FDDEProblem, h, ::DelayPI)
    @unpack f, ϕ, α, τ, tspan = FDDE
    t0 = tspan[1]; T = tspan[2]
    N = ceil(Int, (T-t0)/h)
    t = t0 .+ h*collect(0:N)

    nn_al = collect(0:N).^α
    b = [0; nn_al[2:end].-nn_al[1:end-1]]/gamma(α+1)
    h_al = h^α

    y0 = ϕ(t0)
    y = zeros(N+1)

    g = zeros(N+1)
    y[1] = y0

    for n = 1:N
        tnm1 = t[n]
        if tnm1 <= τ
            y_nm1_tau = ϕ(tnm1-τ)
        else
            nm1_tau1 = floor(Int, (tnm1-τ)/h)
            nm1_tau2 = ceil(Int, (tnm1-τ)/h)
            if nm1_tau1 == nm1_tau2
                y_nm1_tau = y[nm1_tau1+1]
            else
                tt0 = t[nm1_tau1+1]
                tt1 = t[nm1_tau1+2]
                yy0 = y[nm1_tau1+1]
                yy1 = y[nm1_tau1+2]
                y_nm1_tau = ((tnm1-τ)-tt0)/(tt1-tt0)*yy1 + ((tnm1-τ)-tt1)/(tt0-tt1)*yy0
            end
        end

        g[n] = f(tnm1, y[n], y_nm1_tau)
        f_mem = 0
        for j = 0:n-1
            f_mem = f_mem + g[j+1]*b[n-j+1]
        end
        y[n+1] = y0 + h_al*f_mem
    end
    return y
end


text:
using YaoBase: @interface

"""
    parse_block(n, ex)

This function parse the julia object `ex` to a quantum block,
it defines the syntax of high level interfaces. `ex` can be
a function takes number of qubits `n` as input or it can be
a pair.
"""
function parse_block end

parse_block(n::Int, x::Function) = x(n)

function parse_block(n::Int, x::AbstractBlock{N}) where {N}
    n == N || throw(ArgumentError("number of qubits does not match: $x"))
    return x
end

parse_block(n::Int, x::Pair{Int,<:AbstractBlock{N}}) where {N} = x

"""
    prewalk(f, src::AbstractBlock)

Walk the tree and call `f` once the node is visited.
"""
@interface function prewalk(f::Base.Callable, src::AbstractBlock)
    out = f(src)
    for each in subblocks(src)
        prewalk(f, each)
    end
    return out
end

"""
    postwalk(f, src::AbstractBlock)

Walk the tree and call `f` after the children are visited.
"""
@interface function postwalk(f::Base.Callable, src::AbstractBlock)
    for each in subblocks(src)
        postwalk(f, each)
    end
    return f(src)
end

@interface blockfilter!(f, v::Vector, blk::AbstractBlock) = postwalk(x -> f(x) ? push!(v, x) : v, blk)

@interface blockfilter(f, blk) = blockfilter!(f, [], blk)

"""
    collect_blocks(block_type, root)

Return a [`ChainBlock`](@ref) with all block of `block_type` in root.
"""
@interface collect_blocks(::Type{T}, x::AbstractBlock) where {T<:AbstractBlock} =
    blockfilter!(x -> x isa T, T[], x)

#@interface expect(op::AbstractBlock, r::AbstractRegister) = r' * apply!(copy(r), op)

#expect(op::AbstractBlock, dm::DensityMatrix) = mapslices(x->sum(mat(op).*x)[], dm.state, dims=[1,2]) |> vec
expect(op::AbstractBlock, dm::DensityMatrix{1}) = sum(mat(op) .* dropdims(dm.state, dims = 3))

"""
    expect(op::AbstractBlock, reg::AbstractRegister{B}) -> Vector
    expect(op::AbstractBlock, dm::DensityMatrix{B}) -> Vector

expectation value of an operator.
"""
@interface function expect(op::AbstractBlock, dm::DensityMatrix{B}) where {B}
    mop = mat(op)
    [tr(view(dm.state, :, :, i) * mop) for i in 1:B]
end

expect(op::AbstractBlock, reg::AbstractRegister{1}) = reg' * apply!(copy(reg), op)

function expect(op::AbstractBlock, reg::AbstractRegister{B}) where {B}
    ket = apply!(copy(reg), op)
    C = conj!(reshape(ket.state, :, B))
    A = reshape(reg.state, :, B)
    dropdims(sum(A .* C, dims = 1), dims = 1) |> conj
end

function expect(op::Add, reg::AbstractRegister)
    sum(opi -> expect(opi, reg), op)
end

function expect(op::Scale, reg::AbstractRegister)
    factor(op) * expect(content(op), reg)
end

function expect(op, plan::Pair{<:AbstractRegister,<:AbstractBlock})
    expect(op, copy(plan.first) |> plan.second)
end

expect(op::Add, reg::AbstractRegister{1}) = invoke(expect, Tuple{Add,AbstractRegister}, op, reg)
expect(op::Scale, reg::AbstractRegister{1}) = invoke(expect, Tuple{Scale,AbstractRegister}, op, reg)

for FUNC in [:measure!, :measure_collapseto!, :measure_remove!, :measure]
    @eval function YaoBase.$FUNC(
        rng::AbstractRNG,
        op::AbstractBlock,
        reg::AbstractRegister,
        locs::AllLocs;
        kwargs...,
    ) where {B}
        $FUNC(rng::AbstractRNG, eigen!(mat(op) |> Matrix), reg, locs; kwargs...)
    end
end

# obtaining Dense Matrix of a block
LinearAlgebra.Matrix(blk::AbstractBlock) = Matrix(mat(blk))


text:
"""
    Client([; host="127.0.0.1", port=6379, database=0, password="", username="", ssl_config=nothing]) -> Client

Creates a Client instance connecting and authenticating to a Redis host, provide an `MbedTLS.SSLConfig`
(see `get_ssl_config`) for a secured Redis connection (SSL/TLS).

# Attributes
- `host::AbstractString`: Redis host.
- `port::Integer`: Redis port.
- `database::Integer`: Redis database index.
- `password::AbstractString`: Redis password if any.
- `username::AbstractString`: Redis username if any.
- `socket::Union{TCPSocket,MbedTLS.SSLContext}`: Socket used for sending and reveiving from Redis host.
- `lock::Base.AbstractLock`: Lock for atomic reads and writes from client socket.
- `ssl_config::Union{MbedTLS.SSLConfig,Nothing}`: Optional ssl config for secured redis connection.
- `is_subscribed::Bool`: Whether this Client is actively subscribed to any channels or patterns.
- `subscriptions::AbstractSet{<:AbstractString}`: Set of channels currently subscribed on.
- `psubscriptions::AbstractSet{<:AbstractString}`: Set of patterns currently psubscribed on.

# Note
- Connection parameters `host`, `port`, `database`, `password`, `username` will not change after
client istance is constructed, even with `SELECT` or `CONFIG SET` commands.

# Examples
Basic connection:
```julia-repl
julia> client = Client();

julia> set("key", "value"; client=client)
"OK"

julia> get("key"; client=client)
"value"

julia> execute(["DEL", "key"], client)
1
```

SSL/TLS connection:
```julia-repl
julia> ssl_config = get_ssl_config(ssl_certfile="redis.crt", ssl_keyfile="redis.key", ssl_ca_certs="ca.crt");

julia> client = Client(ssl_config=ssl_config);
```
"""
mutable struct Client
    host::AbstractString
    port::Integer
    database::Integer
    password::AbstractString
    username::AbstractString
    socket::Union{TCPSocket,MbedTLS.SSLContext}
    lock::Base.AbstractLock
    ssl_config::Union{MbedTLS.SSLConfig,Nothing}
    is_subscribed::Bool
    subscriptions::AbstractSet{<:AbstractString}
    psubscriptions::AbstractSet{<:AbstractString}
end

function Client(; host="127.0.0.1", port=6379, database=0, password="", username="", ssl_config=nothing)
    if isnothing(ssl_config)
        socket = connect(host, port)
    else
        socket = ssl_connect(host, port, ssl_config)
    end

    client = Client(
        host,
        port,
        database,
        password,
        username,
        socket,
        ReentrantLock(),
        ssl_config,
        false,
        Set{String}(),
        Set{String}()
    )

    !isempty(password * username) && auth(password, username; client=client)
    database != 0 && select(database; client=client)

    return client
end

"""
    get_ssl_config([; ssl_certfile=nothing, ssl_keyfile=nothing, ssl_ca_certs=nothing]) -> MbedTLS.SSLConfig

Loads ssl cert, key and ca cert files from provided directories into MbedTLS.SSLConfig object.

# Examples
```julia-repl
julia> ssl_config = get_ssl_config(ssl_certfile="redis.crt", ssl_keyfile="redis.key", ssl_ca_certs="ca.crt");
```
"""
function get_ssl_config(; ssl_certfile=nothing, ssl_keyfile=nothing, ssl_ca_certs=nothing)
    ssl_config = MbedTLS.SSLConfig(false)

    if !isnothing(ssl_certfile) && !isnothing(ssl_keyfile)
        cert = MbedTLS.crt_parse_file(ssl_certfile)
        key = MbedTLS.parse_keyfile(ssl_keyfile)
        MbedTLS.own_cert!(ssl_config, cert, key)
    end

    if !isnothing(ssl_ca_certs)
        ca_certs = MbedTLS.crt_parse_file(ssl_ca_certs)
        MbedTLS.ca_chain!(ssl_config, ca_certs)
    end

    return ssl_config
end

"""
    ssl_connect(host::AbstractString, port::Integer, ssl_config::MbedTLS.SSLConfig) -> MbedTLS.SSLContext

Connects to the redis host and port, returns a socket connection with ssl context.
"""
function ssl_connect(host::AbstractString, port::Integer, ssl_config::MbedTLS.SSLConfig)
    tcp = connect(host, port)
    io = MbedTLS.SSLContext()
    MbedTLS.setup!(io, ssl_config)
    MbedTLS.associate!(io, tcp)
    MbedTLS.hostname!(io, host)
    MbedTLS.handshake!(io)
    return io
end

"""
    GLOBAL_CLIENT = Ref{Client}()

Reference to a global Client object.
"""
const GLOBAL_CLIENT = Ref{Client}()

"""
    set_global_client(client::Client)
    set_global_client([; host="127.0.0.1", port=6379, database=0, password="", username="", ssl_config=nothing])

Sets a Client object as the `GLOBAL_CLIENT[]` instance.
"""
function set_global_client(client::Client)
    GLOBAL_CLIENT[] = client
end

function set_global_client(; host="127.0.0.1", port=6379, database=0, password="", username="", ssl_config=nothing)
    client = Client(; host=host, port=port, database=database, password=password, username=username,  ssl_config=ssl_config)
    set_global_client(client)
end

"""
    get_global_client() -> Client

Retrieves the `GLOBAL_CLIENT[]` instance, if unassigned then initialises it with default values
`host="127.0.0.1"`, `port=6379`, `database=0`, `password=""`, `username=""`.
"""
function get_global_client()
    if isassigned(GLOBAL_CLIENT)
        return GLOBAL_CLIENT[]
    else
        return set_global_client()
    end
end

"""
    copy(client::Client) -> Client

Creates a new Client instance, copying the connection parameters of the input.
"""
function Base.copy(client::Client)
    return Client(;
        host=client.host,
        port=client.port,
        database=client.database,
        password=client.password,
        username=client.username,
        ssl_config=client.ssl_config
    )
end

"""
    disconnect!(client::Client)

Closes the client socket connection, it will be rendered unusable.
"""
function disconnect!(client::Client)
    close(client.socket)
end

"""
    reconnect!(client::Client) -> Client

Reconnects the input client socket connection.
"""
function reconnect!(client::Client)
    disconnect!(client)

    if isnothing(client.ssl_config)
        new_socket = connect(client.host, client.port)
    else
        new_socket = ssl_connect(connect(client.host, client.port), client.host, client.ssl_config)
    end

    client.socket = new_socket
    !isempty(client.password * client.username) && auth(client.password, client.username; client=client)
    client.database != 0 && select(client.database; client=client)

    return client
end

"""
    flush!(client::Client)

Reads and discards any bytes that remain unread in the client socket.
"""
function flush!(client::Client)
    while bytesavailable(client.socket) > 0
        recv(client.socket)
    end
end

"""
    set_subscribed!(client::Client)

Marks the Client instance as subscribed, should not be used publicly.
"""
function set_subscribed!(client::Client)
    client.is_subscribed = true
end

"""
    set_unsubscribed!(client::Client)

Marks the Client instance as unsubscribed, should not be used publicly.
"""
function set_unsubscribed!(client::Client)
    client.is_subscribed = false
end

"""
    wait_until_subscribed(client::Client)

Blocks until client changes to a subscribed state.
"""
function wait_until_subscribed(client::Client)
    if !client.is_subscribed
        while !client.is_subscribed
            sleep(0.001)
        end
    end
end

"""
    wait_until_unsubscribed(client::Client)

Blocks until client changes to a unsubscribed state.
"""
function wait_until_unsubscribed(client::Client)
    if client.is_subscribed
        while client.is_subscribed
            sleep(0.001)
        end
    end
end

"""
    wait_until_channel_unsubscribed(client::Client[, channels...])

Blocks until client is unsubscribed from channel(s), leave empty to wait until unsubscribed from all channels.
"""
function wait_until_channel_unsubscribed(client::Client, channels...)
    if isempty(channels)
        while !isempty(client.subscriptions)
            sleep(0.001)
        end
    else
        while !isempty(intersect(client.subscriptions, Set{String}(channels)))
            sleep(0.001)
        end
    end
end

"""
    wait_until_pattern_unsubscribed(client::Client[, patterns...])

Blocks until client is unsubscribed from pattern(s), leave empty to wait until unsubscribed from all patterns.
"""
function wait_until_pattern_unsubscribed(client::Client, patterns...)
    if isempty(patterns)
        while !isempty(client.psubscriptions)
            sleep(0.001)
        end
    else
        while !isempty(intersect(client.psubscriptions, Set{String}(patterns)))
            sleep(0.001)
        end
    end
end


text:
export proxL1!, proxL1, normL1

"""
    proxL1!(x::Array{T}, λ::Float64; sparseTrafo::Trafo=nothing, kargs...) where T

performs soft-thresholding - i.e. proximal map for the Lasso problem.

# Arguments:
* `x::Array{T}`                 - Vector to apply proximal map to
* `λ::Float64`                  - regularization paramter
* `sparseTrafo::Trafo=nothing`  - sparsifying transform to apply
"""
function proxL1!(x::AbstractArray{Tc}, λ::T; sparseTrafo::Trafo=nothing, kargs...) where {T, Tc <: Union{T, Complex{T}}}
  ε = eps(T)

  if sparseTrafo != nothing
    z = sparseTrafo*x
    z .= max.((abs.(z).-λ),0) .* (z.+ε)./(abs.(z).+ε)
    x .= adjoint(sparseTrafo)*z
  else
    x .= max.((abs.(x).-λ),0) .* (x.+ε)./(abs.(x).+ε)
  end

  return x
end

"""
    normL1(x::Array{T}, λ::Float64; sparseTrafo::Trafo=nothing, kargs...) where T

returns the value of the L1-regularization term.
Arguments are the same as in `proxL1!`
"""
function normL1(x::T, λ::Float64; sparseTrafo::Trafo=nothing, kargs...) where T<:AbstractArray
  if sparseTrafo != nothing
    l1Norm = λ*norm(sparseTrafo*x,1)
  else
    l1Norm = λ*norm(x,1)
  end
  return l1Norm
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "NetCDFF"
version = v"4.5.3"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/Unidata/netcdf-fortran/archive/refs/tags/v4.5.3.tar.gz", "c6da30c2fe7e4e614c1dff4124e857afbd45355c6798353eccfa60c0702b495a")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/netcdf-fortran*/
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
if [[ "${target}" == *-mingw* ]]; then
    # Manually build shared library for Windows
    gfortran -shared -fPIC -o ${libdir}/libnetcdff.${dlext} -Wl,$(flagon --whole-archive) ${prefix}/lib/libnetcdff.a -Wl,$(flagon --no-whole-archive) -lnetcdf
fi
# Remove static libraries
rm ${prefix}/lib/*.a
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc = "glibc"),
    Platform("x86_64", "macos"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("x86_64", "windows"),
    Platform("i686", "windows"),
]
platforms = expand_gfortran_versions(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libnetcdff", :libnetcdff)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="NetCDF_jll", uuid="7243133f-43d8-5620-bbf4-c2c921802cf3"); compat="400.701.400 - 400.799"),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
    # `MbedTLS_jll` is an indirect dependency through NetCDF, we need to specify
    # a compatible build version for this to work.
    BuildDependency(PackageSpec(; name="MbedTLS_jll", version=v"2.24.0")),
]

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

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
include("..\\libs\\FunctionsRobot.jl")

function modifyChess!(r, N::Integer)
    back_path = BackPath(r, (Nord, West))
    R = NChessRobot(r, Coords(0, 0), N)
    marks!(R, (Sud, Ost))
    back!(R, back_path)
end


text:
mutable struct Tait{P <: QuonParam}
    g::PlanarMultigraph
    quon_params::Dict{Int, P}    # he_id => quon_param
    inputs::Vector{Int}     # list of input vertices
    outputs::Vector{Int}    # list of output vertices
    genuses::Set{Int}   # set of vertices on genuses
    locations::Dict{Int, Tuple{Float64, Float64}}   # v -> location
end

Base.show(io::IO, q::Tait) = print(io, "Tait graph with $(nv(q)) vertices ($(length(genuses(q))) genuses), $(nhe(q)) half edges, $(nf(q)) faces")

nv(q::Tait) = nv(q.g)
ne(q::Tait) = ne(q.g)
nf(q::Tait) = nf(q.g)
nhe(q::Tait) = nhe(q.g)
vertices(q::Tait) = vertices(q.g)
faces(q::Tait) = faces(q.g)
half_edges(q::Tait) = half_edges(q.g)
check_faces(q::Tait) = check_faces(q.g)
check_vertices(q::Tait) = check_vertices(q.g)
check_combinatorial_maps(q::Tait) = check_combinatorial_maps(q.g)
isolated_vertices(q::Tait) = isolated_vertices(q.g)

src(q::Tait, id) = src(q.g, id)
dst(q::Tait, id) = dst(q.g, id)
half_edge(q::Tait, id) = half_edge(q.g, id)
face(q::Tait, id) = face(q.g, id)
next(q::Tait, id) = next(q.g, id)
prev(q::Tait, id) = prev(q.g, id)
twin(q::Tait, id) = twin(q.g, id)
α(q::Tait, id) = α(q.g, id)
ϕ(q::Tait, id) = ϕ(q.g, id)
σ(q::Tait, id) = σ(q.g, id)
σ_inv(q::Tait, id) = σ_inv(q.g, id)
is_boundary(q::Tait, id) = is_boundary(q.g, id)
out_half_edge(q::Tait, id) = out_half_edge(q.g, id)
surrounding_half_edge(q::Tait, id) = surrounding_half_edge(q.g, id)
trace_face(q::Tait, id) = trace_face(q.g, id)
trace_vertex(q::Tait, id) = trace_vertex(q.g, id)
neighbors(q::Tait, id) = neighbors(q.g, id)
rem_face!(q::Tait, id) = rem_face!(q.g, id)
update_face!(q::Tait, id) = update_face!(q.g, id)
is_isolated(q::Tait, id) = is_isolated(q.g, id)
add_vertex!(q::Tait, id) = add_vertex!(q.g, id)
has_vertex(q::Tait, id) = has_vertex(q.g, id)
has_half_edge(q::Tait, id) = has_half_edge(q.g, id)
has_face(q::Tait, id) = has_face(q.g, id)

function add_edge!(q::Tait{P}, v1::Integer, v2::Integer, f::Integer, p::P) where P
    (new_he1, new_he2) = add_edge!(q.g, v1, v2, f)
    if !(0 in (new_he1, new_he2))
        q.quon_params[new_he1] = p
        q.quon_params[new_he2] = p
    end
    return (new_he1, new_he2)
end

function rem_vertex!(q::Tait, v::Integer; update::Bool = true)
    for he in trace_vertex(q, v)
        delete!(q.quon_params, twin(q, he))
        delete!(q.quon_params, he)
    end
    rem_vertex!(q.g, v; update = update)
    deleteat!(q.inputs, findall(isequal(v), q.inputs))
    deleteat!(q.outputs, findall(isequal(v), q.outputs))
    delete!(q.genuses, v)
    delete!(q.locations, v)

    return q
end

function rem_edge!(q::Tait, he_id::Integer; update::Bool = true)
    twin_id = twin(q, he_id)
    rem_edge!(q.g, he_id; update = update)
    delete!(q.quon_params, he_id)
    delete!(q.quon_params, twin_id)

    return q
end

function contract_edge!(q::Tait, he_id::Integer)
    twin_id = twin(q, he_id)
    (v1, v2) = contract_edge!(q.g, he_id)
    if v2 in q.genuses
        delete!(q.genuses, v2)
        push!(q.genuses, v1)
    end
    delete!(q.quon_params, he_id)
    delete!(q.quon_params, twin_id)
    (v1 != v2) && delete!(q.locations, v2)
    return q
end

function merge_graph!(A::Tait, B::Tait; vertical::Bool = true, delta::Float64 = 0.0)
    v_max_A = A.g.v_max
    he_max_A = A.g.he_max
    f_max_A = A.g.f_max
    merge_graph!(A.g, B.g)
    ys_A = [p[2] for p in values(A.locations)]
    y_max_A = maximum(ys_A)
    xs_A = [p[1] for p in values(A.locations)]
    x_max_A = maximum(xs_A)
    for (he_id, p) in B.quon_params
        A.quon_params[he_id + he_max_A] = p
    end
    for v in B.inputs
        push!(A.inputs, v + v_max_A)
    end
    for v in B.outputs
        push!(A.outputs, v + v_max_A)
    end
    for v in B.genuses
        push!(A.genuses, v + v_max_A)
    end
    if vertical
        for (v, p) in B.locations
            A.locations[v + v_max_A] = (p[1], p[2] + y_max_A + 1 - delta)
        end
    else
        for (v, p) in B.locations
            A.locations[v + v_max_A] = (p[1] + x_max_A + 1 - delta, p[2])
        end
    end

    return A
end

Base.copy(q::Tait) = Tait(
        copy(q.g), copy(q.quon_params), copy(q.inputs), copy(q.outputs),
        copy(q.genuses), copy(g.locations)
    )

quon_params(q::Tait) = q.quon_params
quon_param(q::Tait, he_id::Integer) = q.quon_params[he_id]
genuses(q::Tait) = sort!(collect(q.genuses))
is_genus(q::Tait, v::Integer) = (v in q.genuses)
function change_direction!(q::Tait, e_id::Integer)
    p = change_direction(q.quon_params[e_id])
    q.quon_params[e_id] = p
    q.quon_params[twin(q, e_id)] = p
    return q
end

"""
    is_open_vertex(q::Tait, v)

Returns `true` if `v` represent a special vertex for open edges.
"""
is_open_vertex(q::Tait, v::Integer) = (v in q.inputs) || (v in q.outputs)

"""
    is_open_half_edge(q::Tait, he)

Returns `true` if `he` is an open half edge.
"""
is_open_half_edge(q::Tait, he::Integer) = !haskey(q.quon_params, he)

function contract!(A::Tait, B::Tait, va::Vector{Int}, vb::Vector{Int})
    length(va) == length(vb) || error("Size of the input Tait graph and the output Tait graph mismatch")
    v_max_A = A.g.v_max
    merge_graph!(A, B; delta = 3.0)
    vb = vb .+ v_max_A
    contract_boundary_vertices!(A, va, vb)
    return A
end
contract!(A::Tait, B::Tait) = contract!(A, B, copy(A.outputs), copy(B.inputs))
contract!(A, B, C...) = contract!(contract!(A, B), C...)
merge_graph(A::Tait, B::Tait) = merge_graph!(copy(A), B)

function contract_boundary_vertices!(q::Tait, va::Vector{Int}, vb::Vector{Int})
    for (a, b) in zip(va, vb)
        out_a = trace_vertex(q, a)
        out_b = trace_vertex(q, b)
        ia = findfirst(he -> is_boundary(q, he), out_a)
        out_a = [out_a[ia+1:end]; out_a[1:ia]]
        ib = findfirst(he -> is_boundary(q, he), out_b)
        out_b = [out_b[ib+1:end]; out_b[1:ib]]
        reverse!(out_b)

        tp_a = [prev(q, twin(q, he)) for he in out_a]
        tp_b = [prev(q, twin(q, he)) for he in out_b]
        dst_a = [dst(q, he) for he in out_a]
        dst_b = [dst(q, he) for he in out_b]
        fs_a = [face(q, he) for he in out_a]
        fs_b = [face(q, he) for he in out_b]

        # update surrounding half edge
        ts_a = [twin(q, he) for he in out_a]
        ts_b = [twin(q, he) for he in out_b]
        ts = vcat(ts_a, ts_b)
        for he in ts
            f_he = face(q, he)
            he_f = surrounding_half_edge(q, f_he)
            while he_f in ts
                he_f = next(q, he_f)
            end
            q.g.f2he[f_he] = he_f
        end
        if dst_a[end] == dst_b[end] # this will create a self-loop
            # face(out_a[end]) should be 0
            he0 = surrounding_half_edge(q, 0)
            while (he0 in ts) || he0 == out_a[end]
                he0 = next(q, he0)
            end
            q.g.f2he[0] = he0
            q.g.f_max += 1
            new_f = q.g.f_max
            q.g.f2he[new_f] = out_a[end]
            q.g.he2f[out_a[end]] = new_f
        end

        # update faces
        for k in 1:(length(out_a)-1)
            if fs_a[k] != fs_b[k+1]
                for he in trace_face(q, fs_b[k+1])
                    q.g.he2f[he] = fs_a[k]
                end
                delete!(q.g.f2he, fs_b[k+1])
            end
        end

        # connecting 2 planar graphs
        for k in 1:length(out_a)
            t_a = ts_a[k]
            t_b = ts_b[k]
            q.g.twin[out_a[k]] = out_b[k]
            q.g.twin[out_b[k]] = out_a[k]
            q.g.next[tp_a[k]] = out_b[k]
            q.g.next[tp_b[k]] = out_a[k]
            q.g.half_edges[out_a[k]] = HalfEdge(dst_b[k], dst_a[k])
            q.g.half_edges[out_b[k]] = HalfEdge(dst_a[k], dst_b[k])
            q.g.v2he[dst_a[k]] = out_b[k]
            q.g.v2he[dst_b[k]] = out_a[k]
            delete!(q.g.he2f, t_a)
            delete!(q.g.he2f, t_b)
            delete!(q.g.next, t_a)
            delete!(q.g.next, t_b)
            delete!(q.g.twin, t_a)
            delete!(q.g.twin, t_b)
            delete!(q.g.half_edges, t_a)
            delete!(q.g.half_edges, t_b)
        end

        for v_rm in (a, b)
            delete!(q.g.v2he, v_rm)
            deleteat!(q.inputs, findall(isequal(v_rm), q.inputs))
            deleteat!(q.outputs, findall(isequal(v_rm), q.outputs))
            delete!(q.genuses, v_rm)
            delete!(q.locations, v_rm)
        end

        for k in 1:length(out_a)
            contract_edge!(q, out_a[k])
        end
    end

    # make sure the input and output vertices are sorted from left to right
    bds = Int[]
    for he in trace_face(q, 0)
        v = src(q, he)
        if v in q.inputs || v in q.outputs
            push!(bds, v)
        end
    end
    i_ids = findall(v -> (v in q.inputs), bds)
    o_ids = findall(v -> (v in q.outputs), bds)
    if !isempty(i_ids)
        # @show i_ids
        if i_ids[end] - i_ids[1] + 1 == length(i_ids)
            q.inputs[:] = bds[i_ids]
        else
            q.inputs[:] = [bds[(o_ids[end]+1):end]; bds[1:(o_ids[1]-1)]]
        end
    end
    if !isempty(o_ids)
        if o_ids[end] - o_ids[1] + 1 == length(o_ids)
            q.outputs[:] = bds[o_ids]
        else
            q.outputs[:] = [bds[(i_ids[end]+1):end]; bds[1:(i_ids[1]-1)]]
        end
        reverse!(q.outputs)
    end

    return q
end

function right_boundary(q::Tait)
    v_in = q.inputs[end]
    he_in = out_half_edge(q, v_in)
    @assert he_in != 0
    while !is_boundary(q, he_in)
        he_in = σ_inv(q, he_in)
    end
    he_bd = he_in
    right_bd = dst(q, he_bd)
    is_genus(q, right_bd) && return right_bd
end

function left_boundary(q::Tait)
    v_out = q.outputs[1]
    he_out = out_half_edge(q, v_out)
    @assert he_out != 0
    while !is_boundary(q, he_out)
        he_out = σ_inv(q, he_out)
    end
    he_bd = he_out
    left_bd = dst(q, he_bd)
    is_genus(q, left_bd) && return left_bd
end

function tensor_product!(A::Tait, B::Tait;
        right_bd_A::Integer = right_boundary(A),
        left_bd_B::Integer = left_boundary(B))
    v_max_A = A.g.v_max
    merge_graph!(A, B; vertical = false, delta = 1.0)
    merge_boundary_vertices!(A, right_bd_A, left_bd_B + v_max_A)
    return A
end

function merge_boundary_vertices!(q::Tait, v1::Integer, v2::Integer)
    v1_out = out_half_edge(q, v1)
    @assert v1_out != 0
    while !is_boundary(q, v1_out)
        v1_out = σ_inv(q, v1_out)
    end
    v1_in = twin(q, σ(q, v1_out))
    v2_out = out_half_edge(q, v2)
    @assert v2_out != 0
    while !is_boundary(q, v2_out)
        v2_out = σ_inv(q, v2_out)
    end
    v2_in = twin(q, σ(q, v2_out))
    for he_id in trace_vertex(q, v2)
        twin_id = twin(q, he_id)
        he = half_edge(q, he_id)
        q.g.half_edges[he_id] = HalfEdge(v1, he.dst)
        q.g.half_edges[twin_id] = HalfEdge(he.dst, v1)
    end
    q.g.next[v1_in] = v2_out
    q.g.next[v2_in] = v1_out
    delete!(q.g.v2he, v2)
    delete!(q.genuses, v2)
    delete!(q.locations, v2)
    return q
end


text:
mutable struct Entry{ET}
    value::ET
    isinverted::Bool

    function Entry{T}(dims...; static = true) where T
        static ? value = szeros(T,dims...) : value = zeros(T,dims...)
        new{typeof(value)}(value, false)
    end

    Entry(dims...; static = true) = Entry{Float64}(dims...; static = static)
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, entry::Entry)
    println(io, "Entry with value:")
    show(io, mime, entry.value)
end

function Base.zero(::Entry{ET}) where ET
    dims = [ET.parameters[1].parameters...]
    return Entry{ET.parameters[2]}(dims...)
end
function Base.zero(::Type{Entry{ET}}) where ET
    dims = [ET.parameters[1].parameters...]
    return Entry{ET.parameters[2]}(dims...)
end

function randomize!(entry::Entry)
    value = entry.value
    entry.value = randn(eltype(value), size(value))
end


text:
"""
Provides virtual interface for the gray-box environment and simulator.
"""
module GrayBox

using Random
using Distributions

export Simulation,
       Sample,
       Environment,
       EnvironmentSample,
       environment,
       transition!,
       pack,
       unpack

"""
    GrayBox.Simulation

Abstract base type for a gray-box simulation.
"""
abstract type Simulation end

"""
    GrayBox.Sample

Holds sampled value and log-probability.
"""
mutable struct Sample{T}
    value::T
    logprob::Real
end

"""
    GrayBox.Environment

Alias type for a dictionary of gray-box environment distributions.

     e.g., `Environment(:variable_name => Sampleable)`
"""
const Environment = Dict{Symbol, Sampleable}

"""
    GrayBox.EnvironmentSample

Alias type for a single environment sample.

    e.g., `EnvironmentSample(:variable_name => Sample(value, logprob))`
"""
const EnvironmentSample = Dict{Symbol, Sample}

"""
	GrayBox.State

Alias type for a state.
"""
const State = Union{Vector{Float32},Nothing}

"""
    environment(sim::GrayBox.Simulation)

Return all distributions used in the simulation environment.
"""
function environment(sim::Simulation)::Environment end

"""
    transition!(sim::Union{GrayBox.Simulation, GrayBox.Simulation})::Real

Given an input `sample::EnvironmentSample`, apply transition and return the transition log-probability.
"""
function transition!(sim::Simulation, sample::EnvironmentSample)::Real end

"""
    transition!(sim::GrayBox.Simulation)::Real

Apply a transition step, and return the transition log-probability (used with `ASTSeedAction`).
"""
function transition!(sim::Simulation)::Real end

"""
	state(sim::GrayBox.Simulation)::State

Get current state of simulation as a Vector{Float64}.
"""
function state(sim::Simulation)::State end

"""
	count_actions(sim::Simulation)::Int64

Get total dimension of action space, to be implemented by simulator.
"""
function count_actions(sim::Simulation)::Int64 end

"""
	pack(sim::GrayBox.Simulation, actions::Vector{<:Real})::EnvironmentSample
	unpack(sim::GrayBox.Simulation, actions::EnvironmentSample)::Vector{Float32}

Allows general form of actions, delegates translation to/from vector to simulator.
"""
function pack(sim::Simulation, actions::Vector{<:Real})::EnvironmentSample end
function unpack(sim::Simulation, actions::EnvironmentSample)::Vector{Float32} end

end # module GrayBox


text:
module Morsel

using HttpServer,
      HttpCommon,
      Meddle

export App,
       app,
       route,
       namespace,
       with,
       get,
       post,
       put,
       update,
       delete,
       start,
       urlparam,
       routeparam,
       param,
       unsafestring,

       # from HttpCommon
       GET,
       POST,
       PUT,
       UPDATE,
       DELETE,
       OPTIONS,
       HEAD,

       # from Routes
       match_route_handler

include("Routes.jl")

# This produces a dictionary that maps each type of request (GET, POST, etc.)
# to a `RoutingTable`, which is an alias to the `Tree` datatype specified in
# `Trees.jl`.
routing_tables() = (HttpMethodBitmask => RoutingTable)[method => RoutingTable()
                                            for method in HttpMethodBitmasks]

# An 'App' is simply a dictionary linking each HTTP method to a `RoutingTable`.
# The default constructor produces an empty `RoutingTable` for members of
# `HttpMethods`.
#
type App
    routes::Dict{HttpMethodBitmask, RoutingTable}
    state::Dict{Any,Any}
end
function app()
    App(routing_tables(), Dict{Any,Any}())
end

# This defines a route and adds it to the `app.routes` dictionary. As HTTP
# methods are bitmasked integers they can be combined using the bitwise or
# operator, e.g. `GET | POST` refers to a `GET` method and a `POST` method.
#
# Example:
#
#   function hello_world(req, res)
#       "Hello, world!"
#   end
#   route(hello_world, GET | POST, "/hello/world")
#
# Or using do syntax:
#
#   route(app, GET | POST, "/hello/world") do req, res
#       "Hello, world"
#   end
#
function route(handler::Function, app::App, methods::Int, path::String)
    prefix    = get(app.state, :routeprefix, "")
    withstack = get(app.state, :withstack, Midware[])
    handle    = handler
    if length(withstack) > 0
        stack  = middleware(withstack..., Midware( (req::MeddleRequest, res::Response) -> prepare_response(handler(req, res), req, res) ))
        handle = (req::MeddleRequest, res::Response) -> Meddle.handle(stack, req, res)
    end
    for method in HttpMethodBitmasks
        methods & method == method && register!(app.routes[method], prefix * path, handle)
    end
    app
end
route(a::App, m::Int, p::String, h::Function) = route(h, a, m, p)

function namespace(thunk::Function, app::App, prefix::String)
    beforeprefix = get(app.state, :routeprefix, "")
    app.state[:routeprefix] = beforeprefix * prefix
    thunk(app)
    app.state[:routeprefix] = beforeprefix
    app
end

namespace(thunk::Function, app::App, prefix::String, mid::Union(Midware,MidwareStack)...) = with((app) -> namespace(thunk, app, prefix), app, mid...)

function with(thunk::Function, app::App, stack::MidwareStack)
    withstack = get(app.state, :withstack, Midware[])
    beforelen = length(withstack)
    for mid in stack
      push!(withstack, mid)
    end
    app.state[:withstack] = withstack
    thunk(app)
    app.state[:withstack] = withstack[1:beforelen]
    app
end

with(thunk::Function, app::App, mid::Midware...) = with(thunk, app, middleware(mid...))

import Base.get

# These are shortcut functions for common calls to `route`.
# e.g `get` calls `route` with a `GET` as the method parameter.
#
get(h::Function, a::App, p::String)    = route(h, a, GET, p)
post(h::Function, a::App, p::String)   = route(h, a, POST, p)
put(h::Function, a::App, p::String)    = route(h, a, PUT, p)
update(h::Function, a::App, p::String) = route(h, a, UPDATE, p)
delete(h::Function, a::App, p::String) = route(h, a, DELETE, p)

sanitize(input::String) = replace(input,r"</?[^>]*>|</?|>","")
sanitize(x) = x

function validatedvalue(value::Any, validator::Function)
    value == nothing && return nothing
    if is(validator, string)
        value = sanitize(value)
    end
    validator(value)
end

function safelyaccess(req::MeddleRequest, stateKey::Symbol, valKey::Any, validator::Function)
    haskey(req.state, stateKey) ? validatedvalue(get(req.state[stateKey], valKey, nothing), validator) : nothing
end

# validator for getting unsafe ( raw ) input
#
unsafestring(input::String) = input

# Safe accessors for URL parameters, route parameters and POST data
#
function urlparam(req::MeddleRequest, key::String, validator::Function=string)
    safelyaccess(req, :url_params, key, validator)
end
function routeparam(req::MeddleRequest, key::String, validator::Function=string)
    safelyaccess(req, :route_params, key, validator)
end
function param(req::MeddleRequest, key::String, validator::Function=string)
    safelyaccess(req, :data, key, validator)
end
# support symbols...
function urlparam(req::MeddleRequest, key::Symbol, validator::Function=string)
    urlparam(req, string(key), validator)
end
function routeparam(req::MeddleRequest, key::Symbol, validator::Function=string)
    routeparam(req, string(key), validator)
end
function param(req::MeddleRequest, key::Symbol, validator::Function=string)
    param(req, string(key), validator)
end

# `prepare_response` simply sets the data field of the `Response` to the input
# string `s` and calls the middleware's `respond` function.
#
function prepare_response(data::String, req::MeddleRequest, res::Response)
    res.data = data
    respond(req, res)
end
function prepare_response(status::Int, req::MeddleRequest, res::Response)
    res.status = status
    respond(req, res)
end
function prepare_response(data::(Int, String), req::MeddleRequest, res::Response)
    res.status = data[1]
    res.data = data[2]
    respond(req, res)
end
prepare_response(r::Response, req::MeddleRequest, res::Response) = respond(req, r)

# `start` uses to `HttpServer.jl` and `Meddle.jl` packages to launch a webserver
# running `app` on the desired `port`.
#
# This is a blocking function.  Anything that appears after it in the source
# file will not run.
#
function start(app::App, port::Int)

    MorselApp = Midware() do req::MeddleRequest, res::Response
        path = String["/"]
        for comp in split(req.state[:resource], '/')
            !isempty(comp) && push!(path, comp)
        end
        methodizedRouteTable = app.routes[HttpMethodNameToBitmask[req.http_req.method]]
        handler, params = match_route_handler(methodizedRouteTable, path)
        for (k,v) in params
          req.params[symbol(k)] = v
        end
        if handler != nothing
            return prepare_response(handler(req, res), req, res)
        end
        respond(req, Response(404))
    end

    stack = middleware(DefaultHeaders, URLDecoder, CookieDecoder, BodyDecoder, MorselApp)
    http = HttpHandler((req, res) -> Meddle.handle(stack, MeddleRequest(req,Dict{Symbol,Any}(),Dict{Symbol,Any}()), res))
    http.events["listen"] = (port) -> println("Morsel is listening on $port...")

    server = Server(http)
    run(server, port)
end

end # module Morsel


text:
###
# Reflection tooling
##

if VERSION >= v"1.2.0-DEV.249"
    sptypes_from_meth_instance(mi) = Core.Compiler.sptypes_from_meth_instance(mi)
else
    sptypes_from_meth_instance(mi) = Core.Compiler.spvals_from_meth_instance(mi)
end

if VERSION >= v"1.2.0-DEV.320"
    const may_invoke_generator = Base.may_invoke_generator
else
    may_invoke_generator(meth, @nospecialize(atypes), sparams) = isdispatchtuple(atypes)
end

if VERSION < v"1.2.0-DEV.573"
    code_for_method(method, metharg, methsp, world, force=false) = Core.Compiler.code_for_method(method, metharg, methsp, world, force)
else
    code_for_method(method, metharg, methsp, world, force=false) = Core.Compiler.specialize_method(method, metharg, methsp, force)
end

function do_typeinf_slottypes(mi::Core.Compiler.MethodInstance, run_optimizer::Bool, params::Core.Compiler.Params)
    ccall(:jl_typeinf_begin, Cvoid, ())
    result = Core.Compiler.InferenceResult(mi)
    frame = Core.Compiler.InferenceState(result, false, params)
    frame === nothing && return (nothing, Any)
    if Compiler.typeinf(frame) && run_optimizer
        opt = Compiler.OptimizationState(frame)
        Compiler.optimize(opt, result.result)
        opt.src.inferred = true
    end
    ccall(:jl_typeinf_end, Cvoid, ())
    frame.inferred || return (nothing, Any)
    return (frame.src, result.result, frame.slottypes)
end

if :trace_inference_limits in fieldnames(Core.Compiler.Params)
    current_params() = Core.Compiler.CustomParams(ccall(:jl_get_world_counter, UInt, ()); trace_inference_limits=true)
else
    current_params() = Core.Compiler.Params(ccall(:jl_get_world_counter, UInt, ()))
end

function reflect(@nospecialize(F), @nospecialize(TT); optimize=true, params=current_params(), kwargs...)
    sig = Base.signature_type(F, TT)
    reflect(sig; optimize=true, params=params)
end

# TODO: deduplicate with callinfo(sig, rt, ref)
function reflect(@nospecialize(sig); optimize=true, params=current_params())
    methds = Base._methods_by_ftype(sig, -1, params.world)
    (methds === false || length(methds) < 1) && return nothing
    reflections = Reflection[]
    for x in methds
        atypes = x[1]
        sparams = x[2]
        meth = x[3]
        if isdefined(meth, :generator) && !may_invoke_generator(meth, atypes, sparams)
            continue
        end
        mi = code_for_method(meth, sig, x[2], params.world)
        if mi === nothing
            continue
        end
        ref = reflect(mi, optimize=optimize, params=params)
        push!(reflections, ref)
    end
    return reflections
end

function reflect(mi::Core.Compiler.MethodInstance; optimize=true, params=current_params())
    sptypes = sptypes_from_meth_instance(mi)
    (CI, rt, slottypes) = do_typeinf_slottypes(mi, optimize, params)

    ref = Reflection(CI, mi, slottypes, sptypes, params.world)
end


text:
# joSRtoCMO - converts shot record to common-midpoint offset

## helper module
module joSRtoCMO_etc
    using JOLI: jo_convert
    using LinearAlgebra
    # forward SRtoCMO
    function SRtoCMOfwd(x::Vector{vdt},nr::Int,ns::Int,RDT::DataType) where vdt<:Number
        x    = reshape(x,nr,ns);
        B    = zeros(RDT,nr,2*ns-1);
        for k = 1-ns:1:nr-1
        B[floor(Int,abs(k)/2)+1:floor(Int,abs(k)/2)+length(diag(x,k)),ns+k] = diag(x,k);
        end
        B = B[:];
        B=jo_convert(RDT,B,false);
        return B
    end
    # adjoint SRtoCMO
    function SRtoCMOadj(x::Vector{vdt},nr::Int,ns::Int,RDT::DataType) where vdt<:Number
        x = reshape(x,nr,2*ns-1);
        B = zeros(RDT,nr,ns);
        for k = 1-ns:1:nr-1
            B[diagind(B,k)]=vec(x[floor(Int,abs(k)/2)+1:floor(Int,abs(k)/2)+length(diag(B,k)),ns+k]);
        end
        B = B[:];
        B=jo_convert(RDT,B,false);
        return B
    end
end
using .joSRtoCMO_etc

export joSRtoCMO
"""
    julia> joSRtoCMO(M,N)

converter from shot record to common-midpoint offset record

# Signature

    joSRtoCMO(nr::Int,ns::Int;DDT::DataType=joComplex,RDT::DataType=DDT)

# Arguments

- M : ?
- N : ?

# Notes

- ?

# Examples

- ?

"""
function joSRtoCMO(nr::Int,ns::Int;DDT::DataType=joComplex,RDT::DataType=DDT)
    joLinearFunctionFwd_A((nr)*(2*ns-1),nr*ns,
        v1->joSRtoCMO_etc.SRtoCMOfwd(v1,nr,ns,RDT),
        v2->joSRtoCMO_etc.SRtoCMOadj(v2,nr,ns,DDT),
        DDT,RDT;
        name="joSRtoCMO")
end