text:
## A small driver file to run large sample tests combine
#run like this
#julia -p 4 -L testCLTHarness.jl testPOAPLargeSample.jl
#pass arguments for things via ARGS[1] is numRun per batch

spath = "./Results/POAPLargeSample_"
param_path = "./Results/param_portExp_mtn2.csv"
S_grid = collect(1:5:50)
n = 2^17

numRuns = parse(Int, ARGS[1])
dist_type = ARGS[2]

start_time = time_ns()
a = @spawn test_POAPLargeSample(spath, param_path, numRuns, n, S_grid, 8675309, dist_type, 1.0, 100., 1.)
b = @spawn test_POAPLargeSample(spath, param_path, numRuns, n, S_grid, 5164174290, dist_type, 1.0, 100., 1.)
c = @spawn test_POAPLargeSample(spath, param_path, numRuns, n, S_grid, 5167462266, dist_type, 1.0, 100., 1.)
d = @spawn test_POAPLargeSample(spath, param_path, numRuns, n, S_grid, 123456, dist_type, 1.0, 100., 1.)

######
file_a = fetch(a)
file_b = fetch(b)
file_c = fetch(c)
file_d = fetch(d)

time_stamp = (time_ns() - start_time) * 1e-9

##read everyone in, throw away a line
data, header = readdlm(file_a, ',', header=true)

data_t = readdlm(file_b, ',', skipstart=1)
data_t[:, 1] .+= numRuns
data = vcat(data, data_t)

data_t = readdlm(file_c, ',', skipstart=1)
data_t[:, 1] .+= 2numRuns
data = vcat(data, data_t)

data_t = readdlm(file_d, ',', skipstart=1)
data_t[:, 1] .+= 3numRuns
data = vcat(data, data_t)

#strip the name of file_a to make the numbers better
println("Filea \t", file_a)
println("spath \t", spath)

#indx = search(file_a, spath)[end] + 1
indx =  findfirst(spath, file_a)[end] + 1

f = open("$(spath)_$(file_a[indx:end])_full_$(4numRuns).csv", "w")
writedlm(f,  header, ',')
writedlm(f,  data, ',')
close(f)

println("Num Paths: \t $(numRuns) \t Time:", time_stamp )

max_stars_repo_path: Normal Model/testPOAPLargeSample.jl

max_issues_repo_path: Normal Model/testPOAPLargeSample.jl

max_forks_repo_path: Normal Model/testPOAPLargeSample.jl


text:
module RestrictedTriangulationsTests

using Test
using Gridap.ReferenceFEs
using Gridap.Arrays
using Gridap.Geometry

domain = (0,1,0,1)
partition = (10,10)
oldgrid = CartesianGrid(domain,partition)

cell_to_oldcell = collect(1:34)
trian = RestrictedTriangulation(oldgrid,cell_to_oldcell)
test_triangulation(trian)

n_cells = num_cells(trian)

@test restrict(collect(1:n_cells),trian) == cell_to_oldcell
@test reindex(collect(1:n_cells),trian) == cell_to_oldcell
@test get_cell_id(trian) == cell_to_oldcell

oldcell_to_mask = fill(false,num_cells(oldgrid))
oldcell_to_mask[cell_to_oldcell] .= true

trian = RestrictedTriangulation(oldgrid,oldcell_to_mask)
test_triangulation(trian)

@test restrict(collect(1:n_cells),trian) == cell_to_oldcell
@test reindex(collect(1:n_cells),trian) == cell_to_oldcell
@test get_cell_id(trian) == cell_to_oldcell

oldmodel = CartesianDiscreteModel(domain,partition)

trian = Triangulation(oldmodel,oldcell_to_mask)
@test isa(trian,RestrictedTriangulation)
test_triangulation(trian)

@test restrict(collect(1:n_cells),trian) == cell_to_oldcell
@test reindex(collect(1:n_cells),trian) == cell_to_oldcell
@test get_cell_id(trian) == cell_to_oldcell

trian = Triangulation(oldmodel,cell_to_oldcell)
@test isa(trian,RestrictedTriangulation)
test_triangulation(trian)

@test restrict(collect(1:n_cells),trian) == cell_to_oldcell
@test reindex(collect(1:n_cells),trian) == cell_to_oldcell
@test get_cell_id(trian) == cell_to_oldcell

end # module

ext: jl


text:
module QJuliaIntrinsics

import QJuliaRegisters

#create function/type alias
#SSE
m128d   = QJuliaRegisters.double2
m128    = QJuliaRegisters.float4
#AVX/AVX2
m256d   = QJuliaRegisters.double4
m256    = QJuliaRegisters.float8
#AVX3
m512d   = QJuliaRegisters.double8
m512    = QJuliaRegisters.float16

m256dfma(a,b,c) = ccall("llvm.fma.v4f64", llvmcall, m256d, (m256d, m256d, m256d), a, b, c)
m256fma(a,b,c ) = ccall("llvm.fma.v8f32", llvmcall, m256, (m256, m256, m256), a, b, c)

m512dfma(a,b,c) = ccall("llvm.fma.v8f64", llvmcall, m512d, (m512d, m512d, m512d), a, b, c)
m512fma(a,b,c ) = ccall("llvm.fma.v16f32", llvmcall, m512, (m512, m512, m512), a, b, c)

prefetch(a, i1, i2, i3) = ccall("llvm.prefetch", llvmcall, Cvoid, (Ptr{Cchar}, Int32, Int32, Int32), a, i1, i2, i3)

prefetchT0(a) = mmprefetch(a, Int32(0), Int32(3), Int(1))
prefetchT1(a) = mmprefetch(a, Int32(0), Int32(2), Int(1))
prefetchT2(a) = mmprefetch(a, Int32(0), Int32(1), Int(1))

@inline function mm_mul(x::T, y::T) where T <:AbstractFloat
   return x * y
end

@inline function mm_mul(xmm::m128, ymm::m128)
  (VecElement(xmm[1].value * ymm[1].value),
   VecElement(xmm[2].value * ymm[2].value),
   VecElement(xmm[3].value * ymm[3].value),
   VecElement(xmm[4].value * ymm[4].value))
end

@inline function mm_mul(xmm::m128d, ymm::m128d)
  (VecElement(xmm[1].value * ymm[1].value),
   VecElement(xmm[2].value * ymm[2].value))
end

@inline function mm_mul(xmm::m256, ymm::m256)
  (VecElement(xmm[1].value * ymm[1].value),
   VecElement(xmm[2].value * ymm[2].value),
   VecElement(xmm[3].value * ymm[3].value),
   VecElement(xmm[4].value * ymm[4].value),
   VecElement(xmm[5].value * ymm[5].value),
   VecElement(xmm[6].value * ymm[6].value),
   VecElement(xmm[7].value * ymm[7].value),
   VecElement(xmm[8].value * ymm[8].value))
end

@inline function mm_mul(xmm::m256d, ymm::m256d)
  (VecElement(xmm[1].value * ymm[1].value),
   VecElement(xmm[2].value * ymm[2].value),
   VecElement(xmm[3].value * ymm[3].value),
   VecElement(xmm[4].value * ymm[4].value))
end

@inline function mm_add(x::T, y::T) where T <:AbstractFloat
   return x + y
end

@inline function mm_add(xmm::m128, ymm::m128)
  (VecElement(xmm[1].value + ymm[1].value),
   VecElement(xmm[2].value + ymm[2].value),
   VecElement(xmm[3].value + ymm[3].value),
   VecElement(xmm[4].value + ymm[4].value))
end

@inline function mm_add(xmm::m128d, ymm::m128d)
  (VecElement(xmm[1].value + ymm[1].value),
   VecElement(xmm[2].value + ymm[2].value))
end

@inline function mm_add(xmm::m256, ymm::m256)
  (VecElement(xmm[1].value + ymm[1].value),
   VecElement(xmm[2].value + ymm[2].value),
   VecElement(xmm[3].value + ymm[3].value),
   VecElement(xmm[4].value + ymm[4].value),
   VecElement(xmm[5].value + ymm[5].value),
   VecElement(xmm[6].value + ymm[6].value),
   VecElement(xmm[7].value + ymm[7].value),
   VecElement(xmm[8].value + ymm[8].value))
end

@inline function mm_add(xmm::m256d, ymm::m256d)
  (VecElement(xmm[1].value + ymm[1].value),
   VecElement(xmm[2].value + ymm[2].value),
   VecElement(xmm[3].value + ymm[3].value),
   VecElement(xmm[4].value + ymm[4].value))
end

@inline function mm_sub(x::T, y::T) where T <:AbstractFloat
   return x - y
end

@inline function mm_sub(xmm::m128, ymm::m128)
  (VecElement(xmm[1].value - ymm[1].value),
   VecElement(xmm[2].value - ymm[2].value),
   VecElement(xmm[3].value - ymm[3].value),
   VecElement(xmm[4].value - ymm[4].value))
end

@inline function mm_sub(xmm::m128d, ymm::m128d)
  (VecElement(xmm[1].value - ymm[1].value),
   VecElement(xmm[2].value - ymm[2].value))
end

@inline function mm_sub(xmm::m256, ymm::m256)
  (VecElement(xmm[1].value - ymm[1].value),
   VecElement(xmm[2].value - ymm[2].value),
   VecElement(xmm[3].value - ymm[3].value),
   VecElement(xmm[4].value - ymm[4].value),
   VecElement(xmm[5].value - ymm[5].value),
   VecElement(xmm[6].value - ymm[6].value),
   VecElement(xmm[7].value - ymm[7].value),
   VecElement(xmm[8].value - ymm[8].value))
end

@inline function mm_sub(xmm::m256d, ymm::m256d)
  (VecElement(xmm[1].value - ymm[1].value),
   VecElement(xmm[2].value - ymm[2].value),
   VecElement(xmm[3].value - ymm[3].value),
   VecElement(xmm[4].value - ymm[4].value))
end

@inline function mm_mad(a::T, x::T, y::T) where T <:AbstractFloat
   return a * x + y
end

@inline function mm_mad(a::T, xmm::m128, ymm::m128) where T <: AbstractFloat
  (VecElement(a * xmm[1].value + ymm[1].value),
   VecElement(a * xmm[2].value + ymm[2].value),
   VecElement(a * xmm[3].value + ymm[3].value),
   VecElement(a * xmm[4].value + ymm[4].value))
end

@inline function mm_mad(a::T, xmm::m128d, ymm::m128d) where T <: AbstractFloat
  (VecElement(a * xmm[1].value + ymm[1].value),
   VecElement(a * xmm[2].value + ymm[2].value))
end

@inline function mm_mad(a::T, xmm::m256, ymm::m256) where T <: AbstractFloat
   zmm = m256((a, a, a, a, a, a, a, a))
   return m256fma(zmm,xmm,ymm)
end

@inline function mm_mad(a::T, xmm::m256d, ymm::m256d) where T <: AbstractFloat
   zmm = m256d((a, a, a, a))
   return m256dfma(zmm,xmm,ymm)
end

end #QJuliaIntrinsics

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
#lmm.jl

struct LMMLogMsg
    type::Symbol
    msg::String
end
"""
    LMM(model, data; contrasts=Dict{Symbol,Any}(),  random::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, repeated::Union{Nothing, VarEffect} = nothing)

Make Linear-Mixed Model object.

`model`: is a fixed-effect model (`@formula`)

`data`: tabular data

`contrasts`: contrasts for fixed factors

`random`: vector of random effects or single random effect

`repeated`: is a repeated effect (only one)
"""
struct LMM{T} <: MetidaModel
    model::FormulaTerm
    mf::ModelFrame
    mm::ModelMatrix
    covstr::CovStructure{T}
    data::LMMData{T}
    nfixed::Int
    rankx::Int
    result::ModelResult
    maxvcbl::Int
    log::Vector{LMMLogMsg}

    function LMM(model, data; contrasts=Dict{Symbol,Any}(),  random::Union{Nothing, VarEffect, Vector{VarEffect}} = nothing, repeated::Union{Nothing, VarEffect} = nothing)
        #need check responce - Float
        if repeated === nothing && random === nothing
            error("No effects specified!")
        end

        tv = termvars(model)
        if !isnothing(random)
            union!(tv, termvars(random))
        end
        if !isnothing(repeated)
            union!(tv, termvars(repeated))
        end
        data, data_ = StatsModels.missing_omit(NamedTuple{tuple(tv...)}(Tables.columntable(data)))

        lmmlog = Vector{LMMLogMsg}(undef, 0)
        mf     = ModelFrame(model, data; contrasts = contrasts)
        mm     = ModelMatrix(mf)
        nfixed = nterms(mf)
        if repeated === nothing
            repeated = VarEffect(Metida.@covstr(1|1), Metida.ScaledIdentity())
        end
        if random === nothing
            random = VarEffect(Metida.@covstr(0|1), Metida.RZero())
        end
        if !isa(random, Vector) random = [random] end

        if repeated.covtype.s == :SI && !isa(repeated.model, ConstantTerm)
            lmmlog!(lmmlog, 1, LMMLogMsg(:WARN, "Repeated effect not a constant, but covariance type is SI. "))
        end
        lmmdata = LMMData(mm.m, response(mf))
        covstr = CovStructure(random, repeated, data)
        rankx =  rank(mm.m)
        if rankx != size(mm.m, 2)
            lmmlog!(lmmlog, 1, LMMLogMsg(:WARN, "Fixed-effect matrix not full-rank!"))
        end
        new{eltype(mm.m)}(model, mf, mm, covstr, lmmdata, nfixed, rankx, ModelResult(), findmax(length, covstr.vcovblock)[1], lmmlog)
    end
end
"""
    lcontrast(lmm::LMM, i::Int)

L-contrast matrix for `i` fixed effect.
"""
function lcontrast(lmm::LMM, i::Int)
    n = nterms(lmm.mf)
    if i > n || n < 1 error("Factor number out of range 1-$(n)") end
    inds = findall(x -> x==i, lmm.mm.assign)
    if typeof(lmm.mf.f.rhs.terms[i]) <: CategoricalTerm
        mxc   = zeros(size(lmm.mf.f.rhs.terms[i].contrasts.matrix, 1), size(lmm.mm.m, 2))
        mxcv  = view(mxc, :, inds)
        mxcv .= lmm.mf.f.rhs.terms[i].contrasts.matrix
        mx    = zeros(size(lmm.mf.f.rhs.terms[i].contrasts.matrix, 1) - 1, size(lmm.mm.m, 2))
        for i = 2:size(lmm.mf.f.rhs.terms[i].contrasts.matrix, 1)
            mx[i-1, :] .= mxc[i, :] - mxc[1, :]
        end
    else
        mx = zeros(length(inds), size(lmm.mm.m, 2))
        for i = 1:length(inds)
            mx[i, inds[i]] = 1
        end
    end
    mx
end
################################################################################
"""
    thetalength(lmm::LMM)

Length of theta vector.
"""
function thetalength(lmm)
    lmm.covstr.tl
end

"""
    coefn(lmm)

Coef number.
"""
function coefn(lmm)
    length(lmm.result.beta)
end

"""
    theta(lmm::LMM)

Return theta vector.
"""
function theta(lmm::LMM)
    copy(lmm.result.theta)
end
"""
    rankx(lmm::LMM)

Return rank of `X` matrix.
"""
function rankx(lmm::LMM)
    Int(lmm.rankx)
end
################################################################################
function lmmlog!(io, lmmlog::Vector{LMMLogMsg}, verbose, vmsg)
    if verbose == 1
        push!(lmmlog, vmsg)
    elseif verbose == 2
        println(io, vmsg)
        push!(lmmlog, vmsg)
    elseif verbose == 3
        if vmsg.type == :ERROR println(io, vmsg) end
        push!(lmmlog, vmsg)
    end
end
function lmmlog!(lmmlog::Vector{LMMLogMsg}, verbose, vmsg)
    lmmlog!(stdout, lmmlog, verbose, vmsg)
end
function lmmlog!(io, lmm::LMM, verbose, vmsg)
    lmmlog!(io, lmm.log, verbose, vmsg)
end
#MetidaNLopt use this
function lmmlog!(lmm::LMM, verbose, vmsg)
    lmmlog!(stdout, lmm, verbose, vmsg)
end
function lmmlog!(lmm::LMM, vmsg)
    lmmlog!(stdout, lmm, 1, vmsg)
end
################################################################################

function Base.show(io::IO, lmm::LMM)
    println(io, "Linear Mixed Model: ", lmm.model)
    for i = 1:length(lmm.covstr.random)
        println(io, "Random $i: ")
        if lmm.covstr.random[i].covtype.s == :ZERO
            println(io, "   No")
            continue
        end
        println(io, "    Model: $(lmm.covstr.random[i].model === nothing ? "nothing" : string(lmm.covstr.random[i].model, "|", lmm.covstr.random[i].subj))")
        println(io, "    Type: $(lmm.covstr.random[i].covtype.s) ($(lmm.covstr.t[i])), Subjects: $(lmm.covstr.sn[i])")
        #println(io, "   Coefnames: $(coefnames(lmm.covstr.schema[i]))")
    end
    println(io, "Repeated: ")
    println(io, "    Model: $(lmm.covstr.repeated.model === nothing ? "nothing" : string(lmm.covstr.repeated.model, "|", lmm.covstr.repeated.subj))")
    println(io, "    Type: $(lmm.covstr.repeated.covtype.s) ($(lmm.covstr.t[end]))")
    println(io, "    Blocks: $(length(lmm.covstr.vcovblock)), Maximum block size: $(lmm.maxvcbl)")
    #println(io, "")
    if lmm.result.fit
        print(io, "Status: ")
        printresult(io, lmm.result.optim)
        if any(x-> x.type == :ERROR, lmm.log)
            printstyled(io, "  See error(s) in log. Final results can be wrong!\n"; color = :red)
        elseif any(x-> x.type == :WARN, lmm.log)
            printstyled(io, " See warnings in log.\n"; color = :yellow)
        else
            println(io, " (No Errors)")
        end
        #Optim.converged(lmm.result.optim) ? printstyled(io, "converged \n"; color = :green) : printstyled(io, "not converged \n"; color = :red)
        #if length(lmm.log) > 0  printstyled(io, "Warnings! See lmm.log \n"; color = :yellow) end
        ##println(io, "")
        println(io, "    -2 logREML: ", round(lmm.result.reml, sigdigits = 6), "    BIC: ", round(bic(lmm), sigdigits = 6))
        println(io, "")
        println(io, "    Fixed-effects parameters:")
        #println(io, "")
        #chl = '─'
        z = lmm.result.beta ./ lmm.result.se
        mx  = hcat(coefnames(lmm.mf), round.(lmm.result.beta, sigdigits = 6), round.(lmm.result.se, sigdigits = 6), round.(z, sigdigits = 6), round.(ccdf.(Chisq(1), abs2.(z)), sigdigits = 6))
        mx  = vcat(["Name" "Estimate" "SE" "z" "Pr(>|z|)"], mx)
        printmatrix(io, mx)
        println(io, "")
        println(io, "    Variance components:")
        #println(io, "")
        println(io, "    θ vector: ", round.(lmm.result.theta, sigdigits = 6))
        #println(io, "")

        mx = hcat(Matrix{Any}(undef, lmm.covstr.tl, 1), lmm.covstr.rcnames, lmm.covstr.ct, round.(lmm.result.theta, sigdigits = 6))

        for i = 1:length(lmm.covstr.random)
            if lmm.covstr.random[i].covtype.s != :ZERO
                view(mx, lmm.covstr.tr[i], 1) .= "Random $i"
            end
        end
        view(mx, lmm.covstr.tr[end], 1) .= "Residual"
        for i = 1:lmm.covstr.tl
            if mx[i, 3] == :var mx[i, 4] = round.(mx[i, 4]^2, sigdigits = 6) end
        end
        printmatrix(io, mx)
    else
        println(io, "Not fitted.")
    end
end

function printresult(io, res::T) where T <: Optim.MultivariateOptimizationResults
    Optim.converged(res) ? printstyled(io, "converged"; color = :green) : printstyled(io, "not converged"; color = :red)
end
function printresult(io, res)
    if res[3] == :FTOL_REACHED || res[3] == :XTOL_REACHED || res[3] == :SUCCESS
        printstyled(io, "converged ($(string(res[3])))"; color = :green)
    else
        printstyled(io, "not converged ($(string(res[3])))"; color = :red)
    end
end

function Base.show(io::IO, lmmlog::LMMLogMsg)
    if lmmlog.type == :INFO
        printstyled(io, "  INFO  : "; color = :blue)
        println(io, lmmlog.msg)
    elseif lmmlog.type == :WARN
        printstyled(io, "  WARN  : "; color = :yellow)
        println(io, lmmlog.msg)
    elseif lmmlog.type == :ERROR
        printstyled(io, "  ERROR : "; color = :red)
        println(io, lmmlog.msg)
    end
end
"""
    getlog(lmm::LMM)

Return fitting log.
"""
function getlog(lmm::LMM)
    lmm.log
end


text:
import Libxc: evaluate!, Functional

include("lda_x.jl")
include("lda_c_vwn.jl")

# This file extends the evaluate! function from Libxc.jl for cases where the Array type
# is not a plain Julia array and the Floating point type is not Float64
function Libxc.evaluate!(func::Functional, ::Val{:lda}, rho::AbstractArray;
                         zk=nothing, vrho=nothing, v2rho2=nothing)
    func.n_spin == 1  || error("Fallback functionals only for $(func.n_spin) == 1")
    isnothing(v2rho2) || error("Fallback functionals only for 0-th and 1-st derivative")

    zk = reshape(zk, size(rho))
    func.identifier == :lda_x     && return     lda_x!(rho, E=zk, Vρ=vrho)
    func.identifier == :lda_c_vwn && return lda_c_vwn!(rho, E=zk, Vρ=vrho)

    error("Fallback functional for $(string(func.identifier)) not implemented.")
end


text:
abstract type CallbackExecutor end

struct LinearExecutor <: CallbackExecutor end

function handle(executor::LinearExecutor, event::FitEvent, phase, learner)
    idxs = Zygote.ignore() do
        topologicalsort(learner.callbacks.graph)
    end
    for i in idxs
        _on(event, phase, learner.callbacks.cbs[i], learner)
    end
end

function topologicalsort(g::AbstractGraph)
    nodes = Int[]
    unvisited = trues(nv(g))
    while any(unvisited)
        v = findfirst(unvisited)
        dfs(g, v) do v
            push!(nodes, v)
            unvisited[v] = false
        end
    end
    return nodes
end

function dfs(f, g, v, visited = falses(nv(g)))
    for v_ in inneighbors(g, v)
        if !visited[v_]
            dfs(f, g, v_, visited)
        end
    end
    visited[v] = true
    return f(v)
end


text:
########################### Thermal Generation Models ######################################
abstract type AbstractThermalFormulation <: AbstractDeviceFormulation end
abstract type AbstractThermalDispatchFormulation <: AbstractThermalFormulation end
abstract type AbstractThermalUnitCommitment <: AbstractThermalFormulation end
struct ThermalBasicUnitCommitment <: AbstractThermalUnitCommitment end
struct ThermalStandardUnitCommitment <: AbstractThermalUnitCommitment end
struct ThermalDispatch <: AbstractThermalDispatchFormulation end
struct ThermalRampLimited <: AbstractThermalDispatchFormulation end
struct ThermalDispatchNoMin <: AbstractThermalDispatchFormulation end

########################### Active Dispatch Variables ######################################
"""
This function add the variables for power generation output to the model
"""
function activepower_variables!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.ThermalGen}
    if get_warm_start(psi_container.settings)
        initial_value = d -> PSY.get_activepower(d)
    else
        initial_value = nothing
    end
    add_variable(
        psi_container,
        devices,
        variable_name(ACTIVE_POWER, T),
        false,
        :nodal_balance_active;
        ub_value = d -> PSY.get_activepowerlimits(d).max,
        lb_value = d -> PSY.get_activepowerlimits(d).min,
        init_value = initial_value,
    )
    return
end

"""
This function add the variables for power generation output to the model
"""
function reactivepower_variables!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.ThermalGen}
    if get_warm_start(psi_container.settings)
        initial_value = d -> PSY.get_activepower(d)
    else
        initial_value = nothing
    end
    add_variable(
        psi_container,
        devices,
        variable_name(REACTIVE_POWER, T),
        false,
        :nodal_balance_reactive;
        ub_value = d -> PSY.get_reactivepowerlimits(d).max,
        lb_value = d -> PSY.get_reactivepowerlimits(d).min,
        init_value = initial_value,
    )
    return
end

"""
This function add the variables for power generation commitment to the model
"""
function commitment_variables!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
) where {T <: PSY.ThermalGen}
    time_steps = model_time_steps(psi_container)
    if get_warm_start(psi_container.settings)
        initial_value = d -> (PSY.get_activepower(d) > 0 ? 1.0 : 0.0)
    else
        initial_value = nothing
    end

    add_variable(psi_container, devices, variable_name(ON, T), true)
    var_names = (variable_name(START, T), variable_name(STOP, T))
    for v in var_names
        add_variable(psi_container, devices, v, true; init_value = initial_value)
    end

    return
end

activepower_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{<:PSY.ThermalGen},
    model::DeviceModel{<:PSY.ThermalGen, <:AbstractThermalFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::SemiContinuousFF,
) = nothing

activepower_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{<:PSY.ThermalGen},
    model::DeviceModel{<:PSY.ThermalGen, ThermalDispatchNoMin},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::SemiContinuousFF,
) = nothing

"""
This function adds the active power limits of generators when there are no CommitmentVariables
"""
function activepower_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, <:AbstractThermalDispatchFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Nothing,
) where {T <: PSY.ThermalGen}
    constraint_data = Vector{DeviceRange}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        name = PSY.get_name(d)
        limits = PSY.get_activepowerlimits(d)
        range_data = DeviceRange(name, limits)
        add_device_services!(range_data, d, model)
        constraint_data[ix] = range_data
    end
    device_range(
        psi_container,
        constraint_data,
        constraint_name(ACTIVE_RANGE, T),
        variable_name(ACTIVE_POWER, T),
    )
    return
end

"""
This function adds the active power limits of generators when there are CommitmentVariables
"""
function activepower_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, <:AbstractThermalFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Nothing,
) where {T <: PSY.ThermalGen}
    constraint_data = Vector{DeviceRange}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        limits = PSY.get_activepowerlimits(d)
        name = PSY.get_name(d)
        range_data = DeviceRange(name, limits)
        add_device_services!(range_data, d, model)
        constraint_data[ix] = range_data
    end
    device_semicontinuousrange(
        psi_container,
        constraint_data,
        constraint_name(ACTIVE_RANGE, T),
        variable_name(ACTIVE_POWER, T),
        variable_name(ON, T),
    )
    return
end

"""
This function adds the active power limits of generators when there are
    no CommitmentVariables
"""
function activepower_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, ThermalDispatchNoMin},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Nothing,
) where {T <: PSY.ThermalGen}
    constraint_data = Vector{DeviceRange}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        limits = (min = 0.0, max = PSY.get_activepowerlimits(d).max)
        name = PSY.get_name(d)
        range_data = DeviceRange(name, limits)
        add_device_services!(range_data, d, model)
        constraint_data[ix] = range_data
    end

    var_key = variable_name(ACTIVE_POWER, T)
    variable = get_variable(psi_container, var_key)
    # If the variable was a lower bound != 0, not removing the LB can cause infeasibilities
    for v in variable
        if JuMP.has_lower_bound(v)
            JuMP.set_lower_bound(v, 0.0)
        end
    end

    device_range(
        psi_container,
        constraint_data,
        constraint_name(ACTIVE_RANGE, T),
        variable_name(ACTIVE_POWER, T),
    )
    return
end

"""
This function adds the reactive  power limits of generators when there are CommitmentVariables
"""
function reactivepower_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, <:AbstractThermalDispatchFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.ThermalGen}
    constraint_data = Vector{DeviceRange}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        name = PSY.get_name(d)
        limits = PSY.get_reactivepowerlimits(d)
        range_data = DeviceRange(name, limits)
        #add_device_services!(range_data, d, model)
        # Uncomment when we implement reactive power services
        constraint_data[ix] = range_data
    end

    device_range(
        psi_container,
        constraint_data,
        constraint_name(REACTIVE_RANGE, T),
        variable_name(REACTIVE_POWER, T),
    )
    return
end

"""
This function adds the reactive power limits of generators when there CommitmentVariables
"""
function reactivepower_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, <:AbstractThermalFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.ThermalGen}
    constraint_data = Vector{DeviceRange}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        limits = PSY.get_reactivepowerlimits(d)
        name = PSY.get_name(d)
        range_data = DeviceRange(name, limits)
        #add_device_services!(range_data, d, model)
        # Uncomment when we implement reactive power services
        constraint_data[ix] = range_data
    end

    device_semicontinuousrange(
        psi_container,
        constraint_data,
        constraint_name(REACTIVE_RANGE, T),
        variable_name(REACTIVE_POWER, T),
        variable_name(ON, T),
    )
    return
end

### Constraints for Thermal Generation without commitment variables ####
"""
This function adds the Commitment Status constraint when there are CommitmentVariables
"""
function commitment_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, D},
    ::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.ThermalGen, D <: AbstractThermalFormulation, S <: PM.AbstractPowerModel}
    device_commitment(
        psi_container,
        get_initial_conditions(psi_container, ICKey(DeviceStatus, T)),
        constraint_name(COMMITMENT, T),
        (variable_name(START, T), variable_name(STOP, T), variable_name(ON, T)),
    )
    return
end

########################## Make initial Conditions for a Model #############################
function initial_conditions!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    device_formulation::Type{D},
) where {T <: PSY.ThermalGen, D <: AbstractThermalUnitCommitment}
    status_init(psi_container.initial_conditions, devices)
    output_init(psi_container.initial_conditions, devices)
    duration_init(psi_container.initial_conditions, devices)
    return
end

function initial_conditions!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    device_formulation::Type{ThermalBasicUnitCommitment},
) where {T <: PSY.ThermalGen}
    status_init(psi_container.initial_conditions, devices)
    output_init(psi_container.initial_conditions, devices)
    return
end

function initial_conditions!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    device_formulation::Type{D},
) where {T <: PSY.ThermalGen, D <: AbstractThermalDispatchFormulation}
    output_init(psi_container.initial_conditions, devices)
    return
end

########################### Ramp/Rate of Change constraints ################################
"""
This function gets the data for the generators
"""
function _get_data_for_rocc(
    initial_conditions::Vector{InitialCondition},
    resolution::Dates.TimePeriod,
)
    if resolution > Dates.Minute(1)
        minutes_per_period = Dates.value(Dates.Minute(resolution))
    else
        minutes_per_period = Dates.value(Dates.Second(resolution)) / 60
    end
    lenght_devices = length(initial_conditions)
    ini_conds = Vector{InitialCondition}(undef, lenght_devices)
    ramp_params = Vector{UpDown}(undef, lenght_devices)
    minmax_params = Vector{MinMax}(undef, lenght_devices)
    idx = 0
    for ic in initial_conditions
        g = ic.device
        name = PSY.get_name(g)
        non_binding_up = false
        non_binding_down = false
        ramplimits = PSY.get_ramplimits(g)
        rating = PSY.get_rating(g)
        if !isnothing(ramplimits)
            p_lims = PSY.get_activepowerlimits(g)
            max_rate = abs(p_lims.min - p_lims.max) / minutes_per_period
            if (ramplimits.up * rating >= max_rate) & (ramplimits.down * rating >= max_rate)
                @debug "Generator $(name) has a nonbinding ramp limits. Constraints Skipped"
                continue
            else
                idx += 1
            end
            ini_conds[idx] = ic
            ramp_params[idx] = (
                up = ramplimits.up * rating * minutes_per_period,
                down = ramplimits.down * rating * minutes_per_period,
            )
            minmax_params[idx] = p_lims
        end
    end
    if idx < lenght_devices
        deleteat!(ini_conds, (idx + 1):lenght_devices)
        deleteat!(ramp_params, (idx + 1):lenght_devices)
        deleteat!(minmax_params, (idx + 1):lenght_devices)
    end
    return ini_conds, ramp_params, minmax_params
end

"""
This function adds the ramping limits of generators when there are CommitmentVariables
"""
function ramp_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, D},
    system_formulation::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.ThermalGen, D <: AbstractThermalFormulation, S <: PM.AbstractPowerModel}
    time_steps = model_time_steps(psi_container)
    resolution = model_resolution(psi_container)
    key = ICKey(DevicePower, T)
    initial_conditions = get_initial_conditions(psi_container, key)
    rate_data = _get_data_for_rocc(initial_conditions, resolution)
    ini_conds, ramp_params, minmax_params =
        _get_data_for_rocc(initial_conditions, resolution)
    if !isempty(ini_conds)
        # Here goes the reactive power ramp limits when versions for AC and DC are added
        device_mixedinteger_rateofchange(
            psi_container,
            (ramp_params, minmax_params),
            ini_conds,
            constraint_name(RAMP, T),
            (
                variable_name(ACTIVE_POWER, T),
                variable_name(START, T),
                variable_name(STOP, T),
            ),
        )
    else
        @warn "Data doesn't contain generators with ramp limits, consider adjusting your formulation"
    end
    return
end

function ramp_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, D},
    system_formulation::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {
    T <: PSY.ThermalGen,
    D <: AbstractThermalDispatchFormulation,
    S <: PM.AbstractPowerModel,
}
    time_steps = model_time_steps(psi_container)
    resolution = model_resolution(psi_container)
    initial_conditions = get_initial_conditions(psi_container, ICKey(DevicePower, T))
    rate_data = _get_data_for_rocc(initial_conditions, resolution)
    ini_conds, ramp_params, minmax_params =
        _get_data_for_rocc(initial_conditions, resolution)
    if !isempty(ini_conds)
        # Here goes the reactive power ramp limits when versions for AC and DC are added
        device_linear_rateofchange(
            psi_container,
            ramp_params,
            ini_conds,
            constraint_name(RAMP, T),
            variable_name(ACTIVE_POWER, T),
        )
    else
        @warn "Data doesn't contain generators with ramp limits, consider adjusting your formulation"
    end
    return
end

########################### time duration constraints ######################################
"""
If the fraction of hours that a generator has a duration constraint is less than
the fraction of hours that a single time_step represents then it is not binding.
"""
function _get_data_for_tdc(
    initial_conditions_on::Vector{InitialCondition},
    initial_conditions_off::Vector{InitialCondition},
    resolution::Dates.TimePeriod,
)

    steps_per_hour = 60 / Dates.value(Dates.Minute(resolution))
    fraction_of_hour = 1 / steps_per_hour
    lenght_devices_on = length(initial_conditions_on)
    lenght_devices_off = length(initial_conditions_off)
    @assert lenght_devices_off == lenght_devices_on
    time_params = Vector{UpDown}(undef, lenght_devices_on)
    ini_conds = Matrix{InitialCondition}(undef, lenght_devices_on, 2)
    idx = 0
    for (ix, ic) in enumerate(initial_conditions_on)
        g = ic.device
        @assert g == initial_conditions_off[ix].device
        non_binding_up = false
        non_binding_down = false
        timelimits = PSY.get_timelimits(g)
        name = PSY.get_name(g)
        if !isnothing(timelimits)
            if (timelimits.up <= fraction_of_hour) & (timelimits.down <= fraction_of_hour)
                @debug "Generator $(name) has a nonbinding time limits. Constraints Skipped"
            else
                idx += 1
            end
            ini_conds[idx, 1] = ic
            ini_conds[idx, 2] = initial_conditions_off[ix]
            up_val = round(timelimits.up * steps_per_hour, RoundUp)
            down_val = round(timelimits.down * steps_per_hour, RoundUp)
            time_params[idx] = time_params[idx] = (up = up_val, down = down_val)
        end
    end
    if idx < lenght_devices_on
        ini_conds = ini_conds[1:idx, :]
        deleteat!(time_params, (idx + 1):lenght_devices_on)
    end
    return ini_conds, time_params
end

function time_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, D},
    system_formulation::Type{S},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.ThermalGen, D <: AbstractThermalFormulation, S <: PM.AbstractPowerModel}
    parameters = model_has_parameters(psi_container)
    resolution = model_resolution(psi_container)
    initial_conditions_on = get_initial_conditions(psi_container, ICKey(TimeDurationON, T))
    initial_conditions_off =
        get_initial_conditions(psi_container, ICKey(TimeDurationOFF, T))
    ini_conds, time_params =
        _get_data_for_tdc(initial_conditions_on, initial_conditions_off, resolution)
    if !(isempty(ini_conds))
        if parameters
            device_duration_parameters(
                psi_container,
                time_params,
                ini_conds,
                constraint_name(DURATION, T),
                (variable_name(ON, T), variable_name(START, T), variable_name(STOP, T)),
            )
        else
            device_duration_retrospective(
                psi_container,
                time_params,
                ini_conds,
                constraint_name(DURATION, T),
                (variable_name(ON, T), variable_name(START, T), variable_name(STOP, T)),
            )
        end
    else
        @warn "Data doesn't contain generators with time-up/down limits, consider adjusting your formulation"
    end
    return
end

########################### Cost Function Calls#############################################
function cost_function(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::Type{<:AbstractThermalDispatchFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.ThermalGen}
    add_to_cost(psi_container, devices, variable_name(ACTIVE_POWER, T), :variable)
    return
end

function cost_function(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::Type{ThermalDispatchNoMin},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.ThermalGen}
    resolution = model_resolution(psi_container)
    dt = Dates.value(Dates.Minute(resolution)) / 60
    variable = get_variable(psi_container, variable_name(ACTIVE_POWER, T))

    # uses the same cost function whenever there is NO PWL
    function _ps_cost(d::PSY.ThermalGen, cost_component::PSY.VariableCost)
        return ps_cost(
            psi_container,
            variable_name(ACTIVE_POWER, T),
            PSY.get_name(d),
            cost_component,
            dt,
            1.0,
        )
    end

    # This function modified the PWL cost data when present
    function _ps_cost(
        d::PSY.ThermalGen,
        cost_component::PSY.VariableCost{Vector{NTuple{2, Float64}}},
    )
        gen_cost = JuMP.GenericAffExpr{Float64, _variable_type(psi_container)}()
        if !haskey(psi_container.variables, :PWL_cost_vars)
            time_steps = model_time_steps(psi_container)
            container = add_var_container!(
                psi_container,
                :PWL_cost_vars,
                PSY.get_name(d),
                time_steps,
                1:length(cost_component);
                sparse = true,
            )
        else
            container = get_variable(psi_container, :PWL_cost_vars)
        end
        for (t, var) in enumerate(variable[PSY.get_name(d), :])
            pwlvars = JuMP.@variable(
                psi_container.JuMPmodel,
                [i = 1:length(cost_component)],
                base_name = "{$(variable)}_{pwl}",
                start = 0.0,
                lower_bound = 0.0,
                upper_bound = PSY.get_breakpoint_upperbounds(cost_component)[i]
            )
            slopes = PSY.get_slopes(cost_component)
            first_pair = cost_component.cost[1]
            if slopes[1] != 0.0
                slopes[1] =
                    (
                        first_pair[1]^2 - slopes[1] * first_pair[2] +
                        COST_EPSILON * first_pair[2]
                    ) / (first_pair[1] * first_pair[2])
            end
            if slopes[1] < 0 || slopes[1] <= slopes[2]
                throw(IS.ConflictingInputsError("The PWL cost data provided for generator $(PSY.get_name(d)) is not compatible with a No Min Cost."))
            end
            for (ix, pwlvar) in enumerate(pwlvars)
                JuMP.add_to_expression!(gen_cost, slopes[ix] * pwlvar)
                container[(PSY.get_name(d), t, ix)] = pwlvar
            end

            c = JuMP.@constraint(
                psi_container.JuMPmodel,
                variable == sum([pwlvar for (ix, pwlvar) in enumerate(pwlvars)])
            )
            JuMP.add_to_expression!(gen_cost, c)
        end
        return sign * gen_cost * d
    end

    for d in devices
        cost_component = PSY.get_variable(PSY.get_op_cost(d))
        cost_expression = _ps_cost(d, cost_component)
        T_ce = typeof(cost_expression)
        T_cf = typeof(psi_container.cost_function)
        if T_cf <: JuMP.GenericAffExpr && T_ce <: JuMP.GenericQuadExpr
            psi_container.cost_function += cost_expression
        else
            JuMP.add_to_expression!(psi_container.cost_function, cost_expression)
        end
    end

    return
end

function cost_function(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::Type{<:AbstractThermalFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
) where {T <: PSY.ThermalGen}
    #Variable Cost component
    add_to_cost(psi_container, devices, variable_name(ACTIVE_POWER, T), :variable)
    #Commitment Cost Components
    add_to_cost(psi_container, devices, variable_name(START, T), :startup)
    add_to_cost(psi_container, devices, variable_name(STOP, T), :shutdn)
    add_to_cost(psi_container, devices, variable_name(ON, T), :fixed)
    return
end


text:
# Packages
using Random
using LinearAlgebra
using Statistics
using Printf

# outter functions
include("../utils/byrow.jl")
include("../utils/dare.jl")
include("../utils/dlyap.jl")
include("../utils/simula01.jl")
include("../ACQR/ACQR_ksmoother.jl")
include("../ACQR/ACQR_ksmoother_s.jl")

# inner functions
include("ABCDQR_simula.jl")
include("ABCDQR_kfilter.jl")
include("ABCDQR_kfilter_s.jl")
include("ABCDQR_em_starting.jl")
include("ABCDQR_em.jl")
include("ABCDQR_em_s.jl")


text:
# ---
# title: 746. Min Cost Climbing Stairs
# id: problem746
# author: zhwang
# date: 2022-01-19
# difficulty: Easy
# categories: Array, Dynamic Programming
# link: <https://leetcode.com/problems/min-cost-climbing-stairs/description/>
# hidden: true
# ---
#
# On a staircase, the `i`-th step has some non-negative cost `cost[i]` assigned
# (0 indexed).
#
# Once you pay the cost, you can either climb one or two steps. You need to find
# minimum cost to reach the top of the floor, and you can either start from the
# step with index 0, or the step with index 1.
#
# **Example 1:**
#
#
#
#     Input: cost = [10, 15, 20]
#     Output: 15
#     Explanation: Cheapest is start on cost[1], pay that cost and go to the top.
#
#
# **Example 2:**
#
#
#
#     Input: cost = [1, 100, 1, 1, 1, 100, 1, 1, 100, 1]
#     Output: 6
#     Explanation: Cheapest is start on cost[0], and only step on 1s, skipping cost[3].
#
#
# **Note:**
#
#   1. `cost` will have a length in the range `[2, 1000]`.
#   2. Every `cost[i]` will be an integer in the range `[0, 999]`.
#
#
## @lc code=start
using LeetCode

## Dynamic Programming
function min_cost_climbing_stairs(cost::Vector{Int})::Int
    first, second = cost[1], cost[2]
    for c in @view cost[3:end]
        first, second = second, min(first, second) + c
    end
    return min(first, second)
end

## @lc code=end


text:
module day21

parseline(line) = Set.(split.(split(rstrip(replace(line, "," => ""), ')'), " (contains "), " "))

function safeingredients(info, ingredients, allergens)
    allergmap = Dict(a => [i for i in info if a ∈ i[2]] for a in allergens)
    candidates = Dict(a => intersect(map(ing -> ing[1], i)...) for (a, i) in allergmap)
    setdiff(Set(ingredients), union(values(candidates)...))
end

info = parseline.(readlines("day21_input.txt"))
# info = [[Set(["mxmxvkd", "kfcds", "sqjhc", "nhms"]), Set(["dairy", "fish"])],
#         [Set(["trh", "fvjkl", "sbzzf", "mxmxvkd"]), Set(["dairy"])],
#         [Set(["sqjhc", "fvjkl"]), Set(["soy"])],
#         [Set(["sqjhc", "mxmxvkd", "sbzzf"]), Set(["fish"])]]
ingredients = union(map(x -> x[1], info)...)
allergens = union(map(x -> x[2], info)...)
safe = safeingredients(info, ingredients, allergens)
println("Part 1: count = $(sum(map(item -> length(safe ∩ item[1]), info)))")

# Part 2
function knownallergens(info, allergens)
    allergmap = Dict(a => [i for i in info if a ∈ i[2]] for a in allergens)
    candidates = Dict(a => intersect(map(ing -> ing[1], i)...) for (a, i) in allergmap)
    known = Dict()
    while !isempty(candidates)
        merge!(known, filter(x -> length(last(x)) == 1, candidates))
        filter!(x -> first(x) ∉ keys(known), candidates)
        candidates = Dict(a => setdiff(ilist, union(values(known)...)) for (a, ilist) in candidates)
    end
    Dict(a => only(i) for (a, i) in known)
end

known = knownallergens(info, allergens)
println("Part 2: $(join([known[a] for a in sort(collect(allergens))], ','))")

end

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
##
using Colors

function t_colors()
    # RGBA(0.3, 0, 1, 1)
    RGBA(0.3, 1, 1, 1)
end
export t_colors
##
function echo_fin()
    outrs(`echo fin`)
end
export echo_fin
##


text:
using BSON
using RCall

R"load('Ricker_single.rda')"
R"ls()"

@rget y X y_valid X_valid Y_valid tprime

bson("Ricker_single.bson", y = deepcopy(y),
    X = deepcopy(X), y_valid = deepcopy(y_valid),
    X_valid = deepcopy(X_valid), Y_valid = deepcopy(Y_valid),
    tprime = deepcopy(tprime))


text:
# test suite that can be used for all packages inheriting from GPUArrays

module TestSuite

export supported_eltypes

using GPUArrays

using LinearAlgebra
using Random
using Test

using FFTW
using FillArrays

convert_array(f, x) = f(x)
convert_array(f, x::Base.RefValue) = x[]

function compare(f, AT::Type{<:AbstractGPUArray}, xs...; kwargs...)
    cpu_in = convert_array.(copy, xs)
    gpu_in = convert_array.(AT, xs)
    cpu_out = f(cpu_in...; kwargs...)
    gpu_out = f(gpu_in...; kwargs...)
    collect(cpu_out) ≈ collect(gpu_out)
end

function supported_eltypes()
    (Float32, Float64, Int32, Int64, ComplexF32, ComplexF64)
end

include("testsuite/construction.jl")
include("testsuite/gpuinterface.jl")
include("testsuite/indexing.jl")
include("testsuite/io.jl")
include("testsuite/base.jl")
include("testsuite/vector.jl")
include("testsuite/mapreduce.jl")
include("testsuite/broadcasting.jl")
include("testsuite/linalg.jl")
include("testsuite/fft.jl")
include("testsuite/blas.jl")
include("testsuite/random.jl")

"""
Runs the entire GPUArrays test suite on array type `AT`
"""
function test(AT::Type{<:AbstractGPUArray})
    TestSuite.test_construction(AT)
    TestSuite.test_gpuinterface(AT)
    TestSuite.test_indexing(AT)
    TestSuite.test_io(AT)
    TestSuite.test_base(AT)
    TestSuite.test_mapreduce(AT)
    TestSuite.test_broadcasting(AT)
    TestSuite.test_linalg(AT)
    TestSuite.test_fft(AT)
    TestSuite.test_blas(AT)
    TestSuite.test_random(AT)
end

end


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

"""
jltest/EnhancedTestSet.jl extend the types and methods that support unit testing.

Acknowledgements
----------------
* Much of the core functionality of the `EnhancedTestSet` type (and associated methods)
  were taken directly from the TestSetExtensions package developed by Spencer Russell and
  his collaborators (https://github.com/ssfrr/TestSetExtensions.jl). `EnhancedTestSet`
  builds upon the foundation of `ExtendedTestSet` to add support for fail-fast
  functionality.
"""

# --- Exports

export EnhancedTestSet, EnhancedTestSetException

# --- Imports

# Standard library
using Test: Test
using Test: AbstractTestSet, DefaultTestSet, FallbackTestSet
using Test: Result, Fail, Error, Pass

# External Packages
using DeepDiffs
using Distributed

# --- Types

"""
    struct EnhancedTestSet{T<:AbstractTestSet} <: AbstractTestSet

Extension of the `AbstracctTestSet` type that provides the following functionality:

* display diffs (when available) for comparison test failures;

* support fail-fast (i.e., stop testing at first failure).
"""
struct EnhancedTestSet{T<:AbstractTestSet} <: AbstractTestSet
    wrapped::T

    """
    Two inner constructors

    * one for subtypes of AbstractTestSet that possess a constructor with a `desc` argument

    * one for subtypes of AbstractTestSet that do not possess a constructor with a `desc`
      argument
    """
    EnhancedTestSet{T}(desc) where {T} = new(T(desc))

    EnhancedTestSet{FallbackTestSet}(desc) = new(FallbackTestSet())
end

struct EnhancedTestSetException <: Exception
    msg::String
end

# --- Functions/Methods

"""
    EnhancedTestSet(description::AbstractString; <keyword arguments>)

Construct an EnhancedTestSet with the specified `description`.

# Keyword Arguments

* `wrap::Type{<:AbstractTestSet}`: test set type to wrap. Default: `DefaultTestSet`
"""
function EnhancedTestSet(
    description::AbstractString; wrap::Type{<:AbstractTestSet}=DefaultTestSet
)
    return EnhancedTestSet{wrap}(description)
end

function Test.record(ts::EnhancedTestSet{T}, res::Fail) where {T}
    println("\n=====================================================")
    Test.record(ts.wrapped, res)

    return nothing
end

function Test.record(ts::EnhancedTestSet{DefaultTestSet}, res::Fail)
    if Distributed.myid() == 1
        println("\n=====================================================")
        printstyled(ts.wrapped.description, ": "; color=:white)

        if res.test_type === :test
            try
                test_expr = if isa(res.data, Expr)
                    res.data
                elseif isa(res.data, String)
                    Meta.parse(res.data)
                end

                if test_expr.head === :call && test_expr.args[1] === Symbol("==")
                    test_expr_diff =
                        if isa(test_expr.args[2], String) && isa(test_expr.args[3], String)
                            deepdiff(test_expr.args[2], test_expr.args[3])
                        elseif test_expr.args[2].head === :vect &&
                            test_expr.args[3].head === :vect
                            deepdiff(test_expr.args[2].args, test_expr.args[3].args)
                        elseif test_expr.args[2].head === :call &&
                            test_expr.args[3].head === :call &&
                            test_expr.args[2].args[1].head === :curly &&
                            test_expr.args[3].args[1].head === :curly
                            deepdiff(
                                Base.eval(test_expr.args[2].args),
                                Base.eval(test_expr.args[3].args),
                            )
                        end

                    if !isa(test_expr_diff, DeepDiffs.SimpleDiff)
                        # The test was an comparison between things we can diff,
                        # so display the diff
                        printstyled("Test Failed\n"; bold=true, color=Base.error_color())
                        println("  Expression: ", res.orig_expr)
                        printstyled("\n  Diff:\n"; color=Base.info_color())
                        println(test_expr_diff)
                    else
                        # Fallback to the default printing if there is no diff.
                        println(res)
                    end
                else
                    # Fallback to the default printing for non-equality comparisons
                    println(res)
                end
            catch ex
                println(res)
            end
        else
            # Fallback to the default printing if test_type != :test
            println(res)
        end

        Base.show_backtrace(stdout, Test.scrub_backtrace(backtrace()))
        println("\n=====================================================")
    end
    push!(ts.wrapped.results, res)
    return res, backtrace()
end

function Test.record(ts::EnhancedTestSet{T}, res::Error) where {T}
    # Ignore errors generated from failed FallbackTestSet
    if occursin(r"^(Test.)*FallbackTestSetException", res.value) || (
        occursin(r"^(TestTools.jltest.)*EnhancedTestSetException", res.value) &&
        occursin("FallbackTestSetException occurred", res.value)
    )
        throw(EnhancedTestSetException("FallbackTestSetException occurred"))
    end

    println("\n=====================================================")
    Test.record(ts.wrapped, res)
    println("=====================================================")
    return nothing
end

function Test.record(ts::EnhancedTestSet{T}, res::Pass) where {T}
    printstyled("."; color=:green)
    Test.record(ts.wrapped, res)
    return res
end

Test.record(ts::EnhancedTestSet{T}, res) where {T} = Test.record(ts.wrapped, res)

# When recording DefaultTestSet results to an EnhancedTestSet{FallbackTestSet},
# throw an exception if there are any failures or errors in the DefaultTestSet.
#
# Note: this method is only needed for backward compatibility with Julia<=1.3
function Test.record(ts::EnhancedTestSet{FallbackTestSet}, res::DefaultTestSet)
    # Check for failures and errors
    passes, fails, errors, broken, _, _, _, _ = Test.get_test_counts(res)
    if (fails > 0) || (errors > 0)
        throw(
            EnhancedTestSetException(
                "Failure or error occurred in DefaultTestSet nested within FallbackTestSet."
            ),
        )
    end

    return res
end

function Test.finish(ts::EnhancedTestSet{T}) where {T}
    Test.get_testset_depth() == 0 && print("\n\n")
    Test.finish(ts.wrapped)
    return nothing
end


text:
module OverconfidentTraders

using Distributions,
      Parameters,
      Random,
      Statistics,
      Roots,
      LinearAlgebra

using QuantEcon: qnwnorm
# Include Types
include("Types.jl")

# Include Functions
include("Functions.jl")
end # module


text:
module Evolutionary
    using Random, LinearAlgebra, Statistics
    using Base: @kwdef
    using UnPack: @unpack
    using StackViews
    using NLSolversBase: NLSolversBase, AbstractObjective, ConstraintBounds,
                         AbstractConstraints, nconstraints_x, nconstraints
    import NLSolversBase: f_calls, value, value!, value!!
    import Base: show, copy, minimum, summary, identity, getproperty, rand, getindex, length, copyto!, setindex!

    export AbstractStrategy, strategy, mutationwrapper,
           IsotropicStrategy, AnisotropicStrategy, NoStrategy,
           isfeasible, BoxConstraints, apply!, penalty, penalty!, bounds,
           PenaltyConstraints, WorstFitnessConstraints, MixedTypePenaltyConstraints,
           EvolutionaryObjective, ismultiobjective,
           # ES mutations
           gaussian, cauchy,
           # GA mutations
           flip, bitinversion, uniform, BGA, inversion, insertion, swap2, scramble,
           shifting, PM, MIPM, PLM,
           # ES recombinations
           average, marriage,
           # GA recombinations
           SPX, TPX, UX, SHFX,
           DC, AX, WAX, IC, LC, HX, LX, MILX, SBX,
           PMX, OX1, CX, OX2, POS,
           # GA selections
           ranklinear, uniformranking, roulette, rouletteinv, sus, susinv,
           tournament, truncation,
           # DE selections
           random, permutation, randomoffset, best,
           # DE recombinations
           BINX, EXPX,
           # GP exports
           Terminal, subtree, point, hoist, shrink,
           # Optimization methods
           ES, CMAES, GA, DE, TreeGP, NSGA2,
           # re-export
           value, value!, value!!, f_calls

    # optimize API
    include("api/types.jl")
    include("api/objective.jl")
    include("api/results.jl")
    include("api/utilities.jl")
    include("api/constraints.jl")
    include("api/optimize.jl")
    include("api/expressions.jl")
    include("api/moea.jl")

    # Evolution Strategy
    include("es.jl")
    include("cmaes.jl")

    # Genetic Algorithms
    include("ga.jl")
    include("nsga2.jl")

    # Differential Evolution
    include("de.jl")

    # Genetic Programming
    include("api/protected.jl")
    include("gp.jl")

    # ES & GA recombination functions
    include("recombinations.jl")

    # ES & GA mutation functions
    include("mutations.jl")

    # GA selection functions
    include("selections.jl")

    @deprecate uniform(v1, v2) UX(v1, v2)
    @deprecate uniformbin BINX
    @deprecate exponential EXPX
    @deprecate singlepoint SPX
    @deprecate twopoint TPX
    @deprecate domainrange BGA
    @deprecate waverage WAX
    @deprecate intermediate IC
    @deprecate line LC
    @deprecate discrete DC

end


text:
include("../core/blocks.jl")
include("../../configs.jl")

"""
Upsampling and Pyramid Network for the RetinaFace.
"""
struct FPN o6; o2; o3; o4; o5; merge4; merge3; merge2; scale_cnt; end

function FPN(;in_size=2048, dtype=Array{Float32s}, scale_cnt=5)
    return FPN(
        ConvBn(3, 3, in_size, 256, init=xavier, stride=2, bias=false, padding=1, dtype=dtype),

        ConvBn(1, 1, 256, 256, f=relu, bias=false, dtype=dtype),
        ConvBn(1, 1, 512, 256, f=relu, bias=false, dtype=dtype),
        ConvBn(1, 1, 1024, 256, f=relu, bias=false, dtype=dtype),
        ConvBn(1, 1, 2048, 256, f=relu, bias=false, dtype=dtype),

        ConvBn(3, 3, 256, 256, f=relu, bias=false, dtype=dtype, padding=1),
        ConvBn(3, 3, 256, 256, f=relu, bias=false, dtype=dtype, padding=1),
        ConvBn(3, 3, 256, 256, f=relu, bias=false, dtype=dtype, padding=1),
        scale_cnt
    )
end

function (fpn::FPN)(xs; train=true)
    c2, c3, c4, c5 = xs
    # extra conv layer
    if fpn.scale_cnt == 5
        p6 = fpn.o6(c5, train=train)
        # refucing the dim. to 256 and upsampling
        p5 = fpn.o5(c5, train=train)
        p4 = fpn.o4(c4, train=train) + unpool(p5)
        p4 = fpn.merge4(p4, train=train)
        p3 = fpn.o3(c3, train=train) + unpool(p4)
        p3 = fpn.merge3(p3, train=train)
        p2 = fpn.o2(c2, train=train) + unpool(p3)
        p2 = fpn.merge2(p2, train=train)
        return [p2, p3, p4, p5, p6]

    elseif fpn.scale_cnt == 3
        p5 = fpn.o5(c5, train=train)
        p4 = fpn.o4(c4, train=train) + unpool(p5)
        p4 = fpn.merge4(p4, train=train)
        p3 = fpn.o3(c3, train=train) + unpool(p4)
        p3 = fpn.merge3(p3, train=train)
        return [p3, p4, p5]
    end
    return nothing
end

max_stars_repo_name: barisbatuhan/retinaface

max_issues_repo_name: barisbatuhan/retinaface

max_forks_repo_name: barisbatuhan/retinaface


text:
function iscompatible(system::SystemNamelist, cell_parameters::CellParametersCard)
    ibrav, celldm = system.ibrav, system.celldm
    if iszero(ibrav)
        if optionof(cell_parameters) in ("bohr", "angstrom")
            return all(iszero, celldm)
        else  # "alat"
            return !iszero(first(celldm))  # first(celldm) != 0
        end
    else
        return false
    end
end # function iscompatible
iscompatible(x::CellParametersCard, y::SystemNamelist) = iscompatible(y, x)


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     cell_metadata_json: true
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.11.3
#   kernelspec:
#     display_name: Julia 1.6.0
#     language: julia
#     name: julia-1.6
# ---

# +
using ArgoData, DataFrames, CSV, NCDatasets, Plots

fil="ArgoToMITprof.yml"
meta=DownloadArgo.mitprof_interp_setup(fil)
greylist=DataFrame(CSV.File(meta["dirIn"]*"../ar_greylist.txt"));
# -

meta

f=1
println(meta["dirIn"]*meta["fileInList"][f])
argo_data=Dataset(meta["dirIn"]*meta["fileInList"][f])
haskey(argo_data.dim,"N_PROF") ? np=argo_data.dim["N_PROF"] : np=NaN

m=1
prof=DownloadArgo.GetOneProfile(argo_data,m)

k=findall((!ismissing).(prof["T"]))[200]
prof["p"][k]-401.799987792968724

scatter(prof["S"],prof["T"])

# + {"cell_style": "split"}
#scatter(prof["T"],-prof["p"])

# + {"cell_style": "split"}
#scatter(prof["S"],-prof["p"])
# +
lonlatISbad=false
(prof["lat"]<-90.0)|(prof["lat"]>90.0) ? lonlatISbad=true : nothing
(prof["lon"]<-180.0)|(prof["lon"]>360.0) ? lonlatISbad=true : nothing

#if needed then reset lon,lat after issuing a warning
lonlatISbad==true ? println("warning: out of range lon/lat was reset to 0.0,-89.99") : nothing
lonlatISbad ? (prof["lon"],prof["lat"])=(0.0,-89.99) : nothing

#if needed then fix longitude range to 0-360
(~lonlatISbad)&(prof["lon"]>180.0) ? prof["lon"]-=360.0 : nothing

# +
#if needed then convert pressure to depth
(~meta["inclZ"])&(~lonlatISbad) ? DownloadArgo.prof_PtoZ!(prof,meta) : nothing
println(prof[meta["var_out"][1]][200]-398.625084513574966)

#if needed then convert T to potential temperature θ
meta["TPOTfromTINSITU"] ? DownloadArgo.prof_TtoΘ!(prof,meta) : nothing
println(prof["T"][200]-13.80094224720384374)

T_step1=prof["T"]; S_step1=prof["S"]; D_step1=prof["depth"];

#interpolate to standard depth levels
DownloadArgo.prof_interp!(prof,meta)

# + {"cell_style": "center"}
scatter(T_step1,-D_step1,title="temperature")
scatter!(prof["T"],-meta["z_std"])

# + {"cell_style": "center"}
scatter(S_step1,-D_step1,title="salinity")
scatter!(prof["S"],-meta["z_std"])
# -


text:
module EverestMusicLayer

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Music Layer (Everest)" => Ahorn.EntityPlacement(
        Maple.EverestMusicLayerTrigger,
        "rectangle"
    )
)

end


text:
export Random

"""
    Random(;args)

Randomly generated initial data, based on provided (optional arguments) :
- `L` is the typical wavelength (default is `L=1`)
- `s` is the (real) Sobolev index regularity (default is `s=∞`)
- `λ` is the length of spatial localization (default is none)

The initial data `(η,v)` are generated through randomly chosen Fourier coefficients,
multiplied with weigth `w=10^(-|k|L/(2π))` if `s=∞`, or `w=1/(1+9(|k|L/(2π))^(s+1/2))` otherwise.
If `λ` is provided, the function in spatial variables is multiplied by `exp(-|x/λ|^2)`,
and in any case normalized to have maximum absolute value 1.

"""
struct Random <: InitialData

    η
    v

    function Random(;L=1,s=Inf,λ=nothing)

        function generate( x )

            k = Mesh( x ).k
            if s == Inf
                w = 10 .^(-abs.(k*L/(2*π)))
            else
                w = 1 ./( 1 .+ 9*abs.(k*L/(2*π)).^(s+1/2) )
            end
            if λ == nothing
                φ = zero(x).+1
            else
                φ = exp.(-abs.((x/λ).^2))
            end

            θ = 2*π*rand(Float64,length(x))
            r = rand(Float64,length(x))
            Fourier = r.*exp.(-1im.*θ).*w
            Physic = real.(ifft(Fourier)).*φ
            return Physic./maximum(abs.(Physic))
        end
        η( x ) = generate( x )
        v( x ) = generate( x )

    	new( η,v )

    end

end


text:
include("loadtrace.jl")

type dpaResults
  diff_of_means::Array
  res_dir::String
end

#cipher type will determine the sbox to be used for the selectiong function (if at all needed)
#usually the first round is chosen for the DPA attack, however one may extend the attack
#for the other round by extending the if conditionals, by default the first round attack
#is implemented which will work for many SPN as the first round initial step is common for
#many SPN which is addition(xor) with key followed by substituion
#target_bit in the following represents the bit which is to be targeted in the selection function
function dpa(attack_round, target_bit, cipher_type, trace_data)
  unit_size = trace_data.unit_size
  num_unit = trace_data.num_unit
  plaintexts = trace_data.plaintexts
  ciphertexts = trace_data.ciphertexts
  traces = trace_data.traces
  res_dir = trace_data.res_dir
  unit_in_pt = div(length(plaintexts[1]),num_unit)
  keymax = 2^unit_size-1
  sbox = S_boxes[cipher_type]

  n_sample = length(plaintexts)
  diff_vecs_units = []

  if attack_round == "first"
    for i in 1:num_unit
      diff_vecs = Array{Float64}[]
      for subkey in 0:keymax
         marker = zeros(Int,n_sample)
         for p in 1:n_sample
           pt = plaintexts[p]
           ptval = parse(Int, pt[(i-1)*unit_in_pt+1:i*unit_in_pt], 16)
           xor_val = ptval$subkey
           sval = sbox[xor_val+1]
           if bin(sval,unit_size)[target_bit] == '1'
             marker[i] = 1
           end
         end
         #grouping the traces together depending upon
         #which class they belong to 'one' or 'zero'
         diff = mean(traces[find(marker .== 1),:],1)' - mean(traces[find(marker .== 0),:],1)'
         push!(diff_vecs,diff)
       end
       push!(diff_vecs_units,diff_vecs)
     end
   end

   x = dpaResults(diff_vecs_units,res_dir)
 end

max_stars_repo_name: udion/zypher

max_issues_repo_name: udion/zypher

max_forks_repo_name: udion/zypher


text:
@testset "parameterestimation.jl" begin
    include(joinpath("parameterestimation", "bayesian_test.jl"))
    include(joinpath("parameterestimation", "maximumlikelihood_test.jl"))
    include(joinpath("parameterestimation", "probabilityweightedmoment_test.jl"))

end


text:
using Test
using MPI: mpiexec

# run tests on Travis CI in parallel
const TRIXI_TEST = get(ENV, "TRIXI_TEST", "all")
const TRIXI_MPI_NPROCS = clamp(Sys.CPU_THREADS, 2, 3)
const TRIXI_NTHREADS   = clamp(Sys.CPU_THREADS, 2, 3)

@time @testset "Trixi.jl tests" begin
  # This is placed first since tests error out otherwise if `TRIXI_TEST == "all"`,
  # at least on some systems.
  @time if TRIXI_TEST == "all" || TRIXI_TEST == "2d_mpi"
    # Do a dummy `@test true`:
    # If the process errors out the testset would error out as well,
    # cf. https://github.com/JuliaParallel/MPI.jl/pull/391
    @test true

    # Based on `runtests.jl` from `MPI.jl` and `PencilArrays.jl`
    # On Julia v1.5 and before, precompilation is strictly serial and any attempt
    # to use it in parallel will result in race conditions and probably errors.
    # Hence, the additional flag `--compiled-modules=no` is required for Julia
    # versions older than v1.6.
    mpiexec() do cmd
      run(`$cmd -n $TRIXI_MPI_NPROCS $(Base.julia_cmd()) --threads=1 --check-bounds=yes test_examples_2d_parallel.jl`)
    end
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "2d_threaded"
    # Do a dummy `@test true`:
    # If the process errors out the testset would error out as well,
    # cf. https://github.com/JuliaParallel/MPI.jl/pull/391
    @test true

    run(`$(Base.julia_cmd()) --threads=$TRIXI_NTHREADS --check-bounds=yes --code-coverage=none test_examples_2d_parallel.jl`)
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "1d"
    include("test_examples_1d.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "2d_part1"
    include("test_examples_2d_part1.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "2d_part2"
    include("test_examples_2d_part2.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "2d_part3"
    include("test_examples_2d_part3.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "3d_part1"
    include("test_examples_3d_part1.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "3d_part2"
    include("test_examples_3d_part2.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "misc"
    include("test_unit.jl")
    include("test_special_elixirs.jl")
    include("test_visualization.jl")
  end

  @time if TRIXI_TEST == "all" || TRIXI_TEST == "paper-self-gravitating-gas-dynamics"
    include("test_paper-self-gravitating-gas-dynamics.jl")
  end
end


text:
#
# Correctness Tests
#

using Optim

my_tests = [#"test/bfgs.jl",
            "test/curve_fit.jl",
            "test/gradient_descent.jl",
            "test/grid_search.jl",
            "test/l_bfgs.jl",
            "test/levenberg_marquardt.jl",
            "test/newton.jl",
            "test/nelder_mead.jl",
            "test/optimize.jl",
            "test/simulated_annealing.jl"]

println("Running tests:")

for my_test in my_tests
    println(" * $(my_test)")
    include(my_test)
end


text:
mutable struct Assembly{F, L<:LeadTime, C <: Tuple}
    production_cost::F
    capacity::Float64
    leadtime::L
    components::C
    requirements::IdDict{Any, Float64}
    pull_orders::IdDict{Any,Float64}
    batchsize_log::Vector{Float64}
    cost_log::Vector{Float64}
    name::String
end

function Assembly(production_cost, components::Pair{<:Any, <:Number}...; leadtime = LeadTime(0,0), capacity = Inf, name = "assembly")
    @assert hasmethod(production_cost, Tuple{Assembly}) "production_cost must have a method with (::Assembly) signature"
    @assert length(components) >= 1 "Assembly must have at least one component"
    comps = tuple([first(pair) for pair in components]...)
    Assembly(production_cost, Float64(capacity), leadtime, comps, IdDict{Any, Float64}(components), IdDict{Any,Float64}(), zeros(0), zeros(0), name)
end

function Assembly(fixed_order_cost::Number, variable_order_cost::Number, components::Pair{<:Any, <:Number}...; leadtime = LeadTime(0,0), capacity = Inf, name = "assembly")
    Assembly(FixedLinearOrderCost(fixed_order_cost, variable_order_cost), components..., leadtime = leadtime, capacity = capacity, name = name)
end

RLBase.state(a::Assembly) = state(a.leadtime)
state_size(a::Assembly) = state_size(a.leadtime)
print_state(a::Assembly) = print_state(a.leadtime)

function pull!(ass::Assembly, quantity::Number, issuer)
    @assert 0 == length(ass.pull_orders) "assembly received multiple orders"
    push!(ass.pull_orders, issuer => min(quantity, ass.capacity))
    destination, quantity = first(ass.pull_orders)
    for (component, required) in ass.requirements
        pull!(component, quantity*required, ass)
    end
    nothing
end

function Base.push!(ass::Assembly, quantity, source)
    ass.pull_orders[source] = quantity
    nothing
end

function dispatch!(ass::Assembly)
    destination = only(setdiff(keys(ass.pull_orders), (c.inventory for c in ass.components)))
    quantity = ass.pull_orders[destination]
    for (component, required) in ass.requirements
        quantity = min(quantity, ass.pull_orders[component.inventory]/required)
    end
    for (component, required) in ass.requirements
        excess = pop!(ass.pull_orders, component.inventory) - required*quantity
        push!(component.inventory, excess, ass)
    end
    ass.pull_orders[destination] = quantity
    push!(ass.batchsize_log, quantity)
    push!(ass.leadtime, quantity, destination)
    dispatch!(ass.leadtime)
    nothing
end

function reward!(ass::Assembly)
    cost = ass.production_cost(ass)
    push!(ass.cost_log, cost)
    empty!(ass.pull_orders)

    return -cost + reward!(ass.leadtime)
end

inventory_position(ass::Assembly) = inventory_position(ass.leadtime)

function RLBase.reset!(a::Assembly)
    empty!(a.cost_log)
    empty!(a.batchsize_log)
    reset!(a.leadtime)
end

children(ass::Assembly) = ass.components

#cost functions defined in supplier.jl
function (f::FixedLinearOrderCost)(ass::Assembly)
    q = sum(values(ass.pull_orders))
    return (q > 0)*(f.K + f.c*q)
end

(f::LinearOrderCost)(ass::Assembly) = f.c*sum(values(ass.pull_orders))

Base.show(io::IO, ass::Assembly{F}) where {F} = print(io, "Assembly($(ass.name), $F, $(ass.leadtime.leadtime))")

max_stars_repo_path: src/elementary blocks/assembly.jl

max_issues_repo_path: src/elementary blocks/assembly.jl

max_forks_repo_path: src/elementary blocks/assembly.jl


text:
@testset "Test Spline Basis Evaluation" begin

    # constuct specific case that we can compute by hand

    n = 5
    a = 0
    b = 1
    params = SplineParams(range(a,stop=b,length=n), 0, 1)

    x = rand(10000)

    basestr = @inferred BasisMatrices.evalbase(params, x)
    base = Array(basestr)

    # vector of basis functions

    function linBspline(u)
        if u >= 0 && u < 0.25
            return [1-4*u,4*u,0,0,0]

        elseif u >= 0.25 && u < 0.5
            return [0,2*(1 - 2*u),4*u-1,0,0]

        elseif u >= 0.5 && u < 0.75
            return [0,0,3 - 4*u,2*(2*u - 1),0]

        else
            return [0,0,0,4*(1 - u),4*u-3]
        end
    end

    manualbase = Array{Float64}(undef, length(x), n)

    for i in 1:length(x)
        manualbase[i, :] = linBspline(x[i])
    end

    @testset "test evalbase with linear B spline" begin
        @test all(manualbase .== base)
    end

    # test stuff that isn't implemented
    params2 = SplineParams(10, 0, 1, 2)
    @test_throws ErrorException evalbase(SplineSparse, params, nodes(params), -1)
    @test_throws ErrorException evalbase(params, nodes(params), -1)

    @test BasisMatrices.family(SplineParams) == Spline
    @test BasisMatrices.family_name(SplineParams) == "Spline"
    @test issparse(SplineParams)

    @test BasisMatrices.family(params) == Spline
    @test BasisMatrices.family_name(params) == "Spline"
    @test issparse(params)

    @testset "non Float64 eltypes" begin
        p = SplineParams(Float32[0.1, 0.2, 0.3, 0.4, 0.5], 0, 3)
        @test Float32 == eltype(@inferred BasisMatrices.evalbase(p, Float32[0.25, 0.35]))
        @test Float32 == eltype(@inferred BasisMatrices.evalbase(p, Float32[0.25, 0.35], 1))
        @test Float64 == eltype(@inferred BasisMatrices.evalbase(p, [0.25, 0.35]))
    end
end


text:
using ITensors
using ITensorInfiniteMPS
using Test
using Random

#Currently, VUMPS cannot give the right result as the subspace expansion is too small
#This is meant to test the generalized translation rules
@time @testset "vumpsmpo_fqhe" begin
  Random.seed!(1234)

  function fermions_space_shift(pos; p=1, q=1, conserve_momentum=true, momentum_shift=1)
    if !conserve_momentum
      return [QN("Nf", -p) => 1, QN("Nf", q - p) => 1]
    else
      return [
        QN(("Nf", -p), ("NfMom", -p * pos + momentum_shift)) => 1,
        QN(("Nf", q - p), ("NfMom", (q - p) * pos + momentum_shift)) => 1,
      ]
    end
  end

  function initstate(n)
    if mod(n, 3) == 1
      return 2
    else
      return 1
    end
  end

  # VUMPS arguments
  cutoff = 1e-8
  maxdim = 3
  tol = 1e-8
  maxiter = 50
  outer_iters = 1

  #### Test 1/3
  model = Model"fqhe_2b_pot"()
  model_kwargs = (Ly=3.0, Vs=[1.0], prec=1e-6) #parent Hamiltonian of Laughlin 13
  p = 1
  q = 3
  conserve_momentum = true
  momentum_shift = 1

  #ψ1 = InfMPS(s, n->[1, 2, 2, 1, 1, 1][n]);

  @testset "VUMPS/TDVP with: multisite_update_alg = $multisite_update_alg, conserve_qns = $conserve_qns, nsites = $nsites" for multisite_update_alg in
                                                                                [
      "sequential"
    ],
    conserve_qns in [true],
    nsites in [6],
    time_step in [-Inf]

    vumps_kwargs = (
      multisite_update_alg=multisite_update_alg,
      tol=tol,
      maxiter=maxiter,
      outputlevel=0,
      time_step=time_step,
    )
    subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)

    function fermion_momentum_translator(i::Index, n::Integer; N=nsites)
      ts = tags(i)
      translated_ts = ITensorInfiniteMPS.translatecelltags(ts, n)
      new_i = replacetags(i, ts => translated_ts)
      for j in 1:length(new_i.space)
        ch = new_i.space[j][1][1].val
        mom = new_i.space[j][1][2].val
        new_i.space[j] = Pair(
          QN(("Nf", ch), ("NfMom", mom + n * N * ch)), new_i.space[j][2]
        )
      end
      return new_i
    end

    fermionic_space = [
      fermions_space_shift(
        x; p=p, q=q, conserve_momentum=conserve_momentum, momentum_shift=momentum_shift
      ) for x in 1:nsites
    ]
    s = infsiteinds(
      "Fermion", nsites; space=fermionic_space, translator=fermion_momentum_translator
    )
    ψ = InfMPS(s, initstate)

    Hmpo = InfiniteMPOMatrix(model, s; model_kwargs...)
    # Alternate steps of running VUMPS and increasing the bond dimension
    ψ = tdvp(Hmpo, ψ; vumps_kwargs...)
    for _ in 1:outer_iters
      println("Subspace expansion")
      ψ = @time subspace_expansion(ψ, Hmpo; subspace_expansion_kwargs...)
      println("TDVP")
      ψ = @time tdvp(Hmpo, ψ; vumps_kwargs...)
    end
  end
end


text:
function ades!(jj::Int64, input::SparseVector{Int64}, fock::SparseVector{Int64})

    # fock is a fock state, a vector having integers
    # its last element is square of coefficient
    fock .= input
    focksize = length(fock)

    fock[jj] = fock[jj]-1
    fock[focksize] = input[focksize]*(fock[jj]+1)

    if fock[jj] < 0
       fock[focksize] = 0
    end

    # return fock

end

max_stars_repo_name: leios/SOCED

max_issues_repo_name: leios/SOCED

max_forks_repo_name: leios/SOCED


text:
using CosmoMMF
using Documenter

DocMeta.setdocmeta!(CosmoMMF, :DocTestSetup, :(using CosmoMMF); recursive=true)

makedocs(;
    modules=[CosmoMMF],
    authors="James Sunseri",
    repo="https://github.com/James11222/CosmoMMF.jl/blob/{commit}{path}#{line}",
    sitename="CosmoMMF.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://James11222.github.io/CosmoMMF.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/James11222/CosmoMMF.jl",
    devbranch = "main"
)


text:
export add, ∑

"""
    add(x::Int64, y::Int64)

Returns x + y

# Examples

```jldoctest
julia> x, y = 1, 3; add(x, y)
4

julia> @assert "add(x::Int64, y::Int64)" == which_sig(add, x, y)
```
"""
add(x::Int64, y::Int64) = x + y

"""
    add(x::T, y::T) where {T<:Number}

Returns x + y

# Examples

```jldoctest
julia> x, y = 1+3im, 3+5im; add(x, y)
4 + 8im

julia> @assert "add(x::T, y::T) where T<:Number" == which_sig(add, x, y)
```
"""
add(x::T, y::T) where {T<:Number} = x + y

"""
    add(x::T, y::T) where {T<:AbstractString}

Returns x*y

# Examples

```jldoctest
julia> x, y = "my", "string"; add(x, y)
"mystring"

julia> @assert "add(x::T, y::T) where T<:AbstractString" == which_sig(add, x, y)
```
"""
add(x::T, y::T) where {T<:AbstractString} = x*y

"""
    add(x, y)

Returns x + y

# Examples

```jldoctest
julia> x, y = Vector([1,2,3]), Vector([4,5,6]); add(x, y)
3-element Array{Int64,1}:
 5
 7
 9

julia> @assert "add(x, y)" == which_sig(add, x, y)
```
"""
add(x,y) = x + y

@doc raw"""
    ∑(x::Int64...)
Adds all parameters:

```math
\sum_i^n x_1, x_2, \dots x_n \forall x_k \in \mathbb{N}
```
# Examples

```jldoctest
julia> x1, x2, x3, x4, x5 = 1:5; ∑(x1,x2,x3,x4,x5)
15

julia> @assert "∑(x::Int64...)" == which_sig(∑, x1, x2, x3, x4, x5)
```
"""
∑(x::Int64...) = sum(x)

@doc raw"""
    ∑(x::Int64[])
Adds all parameters:

# Examples

```jldoctest
julia> x = Vector(1:5); ∑(x)
15

julia> @assert "∑(x::AbstractArray{Int64,1})" == which_sig(∑, x)

```
"""
∑(x::AbstractArray{Int64,1}) = ∑(x...)

max_stars_repo_name: binnisb/JuliaLearn

max_issues_repo_name: binnisb/JuliaLearn

max_forks_repo_name: binnisb/JuliaLearn


text:
@testset "C" begin
    @testset "ccifrm" begin
        @test ccifrm(2, 3000) == (13000, "ITRF93", 399)
        @test ccifrm(9999, 3000) === nothing
    end
    @testset "cgv2el" begin
        vec1 = [1.0, 1.0, 1.0]
        vec2 = [1.0, -1.0, 1.0]
        cent = [-1.0, 1.0, -1.0]
        ellipse = cgv2el(cent, vec1, vec2)
        @test center(ellipse) ≈ [-1.0, 1.0, -1.0]
        @test semi_major(ellipse) ≈ [sqrt(2.0), 0.0, sqrt(2.0)]
        @test semi_minor(ellipse) ≈ [0.0, sqrt(2.0), 0.0]
    end
    @testset "chbder" begin
        cp = [1.0, 3.0, 0.5, 1.0, 0.5, -1.0, 1.0]
        x2s = [0.5, 3.0]
        dpdxs = chbder(cp, x2s, 1.0, 3)
        expected = [-0.340878, 0.382716, 4.288066, -1.514403]
        @test dpdxs ≈ expected rtol=1e-6
    end
    @testset "cidfrm" begin
        @test cidfrm(399) == (10013, "IAU_EARTH")
        @test cidfrm(999999) === nothing
    end
    @testset "ckcov" begin
        try
            sclk = path(CASSINI, :sclk)
            ck = path(CASSINI, :ck)
            furnsh(sclk)
            ckid = ckobj(ck)[1]
            @test length(ckid) == 1
            cover = ckcov(ck, ckid, false, "INTERVAL", 0.0, "SCLK")
            @test cover[1:2] ≈ [267832537952.000000, 267839247264.000000]
            @test cover[3:4] ≈ [267839256480.000000, 267867970464.000000]
            @test cover[5:6] ≈ [267868006304.000000, 267876773792.000000]
        finally
            kclear()
        end
    end
    @testset "ckgp" begin
        try
            ck = path(CASSINI, :ck)
            furnsh(ck, path(CASSINI, :sclk), path(CASSINI, :ik),
                   path(CASSINI, :fk), path(CASSINI, :pck))
            ckid = ckobj(ck)[1]
            cover = ckcov(ck, ckid, false, "INTERVAL", 0.0, "SCLK")
            cmat, clkout = ckgp(ckid, cover[1], 256, "J2000")
            expected_cmat = [0.5064665782997639365 -0.75794210739897316387 0.41111478554891744963;
                             -0.42372128242505308071 0.19647683351734512858 0.88422685364733510927;
                             -0.7509672961490383436 -0.6220294331642198804 -0.22164725216433822652]
            @test cmat ≈ expected_cmat
            @test clkout == 267832537952.0
        finally
            kclear()
        end
    end
    @testset "ckgpav" begin
        try
            ck = path(CASSINI, :ck)
            furnsh(ck, path(CASSINI, :sclk), path(CASSINI, :ik),
                   path(CASSINI, :fk), path(CASSINI, :pck))
            ckid = ckobj(ck)[1]
            cover = ckcov(ck, ckid, false, "INTERVAL", 0.0, "SCLK")
            cmat, av, clkout = ckgpav(ckid, cover[1], 256, "J2000")
            expected_cmat = [0.5064665782997639365 -0.75794210739897316387 0.41111478554891744963;
                             -0.42372128242505308071 0.19647683351734512858 0.88422685364733510927;
                             -0.7509672961490383436 -0.6220294331642198804 -0.22164725216433822652]
            expected_av = [-0.00231258422150853885, -0.00190333614370416515, -0.00069657429072504716]
            @test cmat ≈ expected_cmat
            @test av ≈ expected_av
            @test clkout == 267832537952.0
        finally
            kclear()
        end
    end
    @testset "cklpf/ckupf" begin
        file = tempname()
        try
            ifname = "Test CK type 1 segment created by cklpf"
            handle = ckopn(file, ifname, 10)
            ckw01(handle, 1.1, 4.1, -77701, "J2000", "Test type 1 CK segment",
                  [1.1, 4.1], [[1.0, 1.0, 1.0, 1.0], [2.0, 2.0, 2.0, 2.0]],
                  [[0.0, 0.0, 1.0], [0.0, 0.0, 2.0]])
            ckcls(handle)
            kclear()
            handle = cklpf(file)
            ckupf(handle)
            ckcls(handle)
            @test isfile(file)
        finally
            kclear()
        end
    end
    @testset "ckw01" begin
        file = tempname()
        try
            handle = ckopn(file)
            init_size = filesize(file)
            quats = [ones(4), ones(4)]
            avvs = [ones(3), ones(3)]
            sclkdp = [1.0, 2.0]
            ckw01(handle, 1.0, 2.0, -7701, "J2000", "GLL SCAN PLT - DATA TYPE 1", sclkdp, quats, avvs)
            ckcls(handle)
            end_size = filesize(file)
            @test end_size != init_size
        finally
            kclear()
        end
    end
    @testset "ckw02" begin
        try
            ck2 = tempname()
            inst = -77702
            maxrec = 201
            secpertick = 0.001
            segid = "test type 2 ck segment"
            ifname = "test ck type 2 segment created by cspice_ckw02"
            ncomch = 0
            ref = "j2000"
            spacing_ticks = 10.0
            spacing_secs = spacing_ticks * secpertick
            rate = 0.01
            handle = ckopn(ck2, ifname, ncomch)
            init_size = filesize(ck2)
            quats = fill([0.0, 0.0, 0.0, 1.0], maxrec)
            av = fill(zeros(3), maxrec)
            rates = fill(secpertick, maxrec)
            sclkdp = collect(1:maxrec) .* spacing_ticks
            sclkdp .+= 1000.0
            starts = sclkdp
            stops = sclkdp .+ (0.8 * spacing_ticks)
            begtime = sclkdp[1]
            endtime = sclkdp[end]
            ckw02(handle, begtime, endtime, inst, ref, segid, starts, stops, quats, av, rates)
            ckcls(handle)
            end_size = filesize(ck2)
            @test end_size != init_size
        finally
            kclear()
        end
    end
    @testset "ckw03" begin
        try
            ck3 = tempname()
            inst = -77702
            ref = "j2000"
            maxrec = 201
            secpertick = 0.001
            segid = "test type 3 ck segment"
            ifname = "test ck type 3 segment created by cspice_ckw03"
            spacing_ticks = 10.0
            spacing_secs = spacing_ticks * secpertick
            rate = 0.01
            handle = ckopn(ck3, ifname, 0)
            init_size = filesize(ck3)
            quats = fill([0.0, 0.0, 0.0, 1.0], maxrec)
            sclkdp = collect(1:maxrec) .* spacing_ticks
            sclkdp .+= 1000.0
            starts = sclkdp[1:2:end-1]
            begtime = sclkdp[1]
            endtime = sclkdp[end]
            ckw03(handle, begtime, endtime, inst, ref, segid, sclkdp, quats, starts)
            ckcls(handle)
            end_size = filesize(ck3)
            @test end_size != init_size
        finally
            kclear()
        end
    end
    @testset "ckw05" begin
        try
            ck5 = tempname()
            epochs = collect(0.0:1.0)
            inst = [-41000, -41001, -41002, -41003]
            segid = "ck type 05 test segment"
            avflag = true
            type0data = [[ 9.999e-1, -1.530e-4, -8.047e-5, -4.691e-4, 0.0, 0.0, 0.0, 0.0],
                         [9.999e-1, -4.592e-4, -2.414e-4, -1.407e-3,
                          -7.921e-10, -1.616e-7, -8.499e-8,  -4.954e-7]]
            type1data = [[ 9.999e-1, -1.530e-4, -8.047e-5, -4.691e-4],
                         [9.999e-1, -4.592e-4, -2.414e-4, -1.407e-3]]
            type2data = [[0.959, -0.00015309, -8.0476e-5, -0.00046913,
                          0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                         [0.959, -0.00045928, -0.00024143, -0.0014073,
                          -7.921e-10, -1.616e-7, -8.499e-8, -4.954e-7,
                          3.234e-7, 1.7e-7, 9.91e-7, 3.234e-7, 1.7e-9, 9.91e-9]]
            type3data = [[0.959, -0.00015309, -8.0476e-05, -0.00046913, 0.0, 0.0, 0.0],
                         [0.959, -0.00045928, -0.00024143, -0.0014073, 3.234e-7, 1.7e-7, 9.91e-7]]

            handle = ckopn(ck5, " ", 0)
            init_size = filesize(ck5)
            # test subtype 0
            ckw05(handle, 0, 15, epochs[1], epochs[end], inst[1], "j2000", avflag, segid, epochs,
                  type0data, 1000.0, 1, epochs)
            # test subtype 1
            ckw05(handle, 1, 15, epochs[1], epochs[end], inst[2], "j2000", avflag, segid, epochs,
                  type1data, 1000.0, 1, epochs)
            # test subtype 2
            ckw05(handle, 2, 15, epochs[1], epochs[end], inst[3], "j2000", avflag, segid, epochs,
                  type2data, 1000.0, 1, epochs)
            # test subtype 3
            ckw05(handle, 3, 15, epochs[1], epochs[end], inst[4], "j2000", avflag, segid, epochs,
                  type3data, 1000.0, 1, epochs)
            ckcls(handle)

            end_size = filesize(ck5)
            @test end_size != init_size
            furnsh(ck5)
            cmat, av, clk = ckgpav(-41000, epochs[1 ]+ 0.5, 1.0, "j2000")
            @test clk ≈ 0.5
        finally
            kclear()
        end
    end
    @testset "clight" begin
        @test clight() == 299792.458
    end
    @testset "clpool" begin
        pdpool("test", [1.0, 2.0, 3.0])
        clpool()
        @test gdpool("test") === nothing
    end
    @testset "cmprss" begin
        s = "1,,2,,,3,,,,"
        @test cmprss(',', 1, s) == "1,2,3,"
        @test cmprss(",", 2, s) == "1,,2,,3,,"
        @test cmprss(",", 3, s) == "1,,2,,,3,,,"
    end
    @testset "cnmfrm" begin
        @test cnmfrm("norbert") === nothing
        code, name = cnmfrm("IO")
        @test code == 10023
        @test name == "IAU_IO"
    end
    @testset "conics" begin
        try
            furnsh(path(CORE, :lsk), path(CORE, :spk), path(CORE, :gm_pck))
            et = str2et("Dec 25, 2007")
            state, ltime = spkezr("Moon", et, "J2000", "NONE", "EARTH")
            mu = bodvrd("EARTH", "GM")[1]
            elts = oscelt(state, et, mu)
            later = et + 7.0 * spd()
            later_state = conics(elts, later)
            state, ltime = spkezr("Moon", later, "J2000", "NONE", "EARTH")
            pert = later_state .- state
            expected_pert = [-7.48885583081946242601e+03,
                            3.97608014470621128567e+02,
                            1.95744667259379639290e+02,
                            -3.61527427787390887026e-02,
                            -1.27926899069508159812e-03,
                            -2.01458906615054056388e-03]
            @test expected_pert ≈ pert
        finally
            kclear()
        end
    end
    @testset "convrt" begin
        @test convrt(300.0, "statute_miles", "km") == 482.8032
        @test convrt(1.0, "parsecs", "lightyears") ≈ 3.2615638
        @test_throws SpiceError convrt(1.0, "parsecs", "hours")
        @test_throws SpiceError convrt(1.0, "parsecs", "norbert")
    end
    @testset "cpos/cposr" begin
        string = "BOB, JOHN, TED, AND MARTIN...."
        @test cpos(string, " ,", 1) == 4
        @test cpos(string, " ,", 5) == 5
        @test cpos(string, " ,", 6) == 10
        @test cpos(string, " ,", 11) == 11
        @test cpos(string, " ,", 12) == 15
        @test cpos(string, " ,", 16) == 16
        @test cpos(string, " ,", 17) == 20
        @test cpos(string, " ,", 21) == -1
        @test cpos(string, " ,", -112) == 4
        @test cpos(string, " ,", -1) == 4
        @test cpos(string, " ,", 1230) == -1

        @test cposr(string, " ,", 30) == 20
        @test cposr(string, " ,", 26) == 20
        @test cposr(string, " ,", 19) == 16
        @test cposr(string, " ,", 15) == 15
        @test cposr(string, " ,", 14) == 11
        @test cposr(string, " ,", 10) == 10
        @test cposr(string, " ,", 9) == 5
        @test cposr(string, " ,", 4) == 4
        @test cposr(string, " ,", 3) == -1
        @test cposr(string, " ,", 231) == 20
        @test cposr(string, " ,", 31) == 20
        @test cposr(string, " ,", -1) == -1
        @test cposr(string, " ,", -10) == -1
    end
    @testset "cvpool" begin
        try
            pdpool("TEST_VAR_CVPOOL", [-646.0])
            swpool("TEST_CVPOOL", ["TEST_VAR_CVPOOL"])
            pdpool("TEST_VAR_CVPOOL", [565.0])
            updated = cvpool("TEST_CVPOOL")
            @test updated
            clpool()
        finally
            kclear()
        end
    end
    @testset "cyllat" begin
        @test collect(cyllat(1.0, deg2rad(180.0), -1.0)) ≈ [sqrt(2), π, -π/4.0]
    end
    @testset "cylrec" begin
        @test cylrec(0.0, deg2rad(33.0), 0.0) ≈ [0.0, 0.0, 0.0]
    end
    @testset "cylsph" begin
        @test collect(cylsph(1.0, deg2rad(180.0), 1.0)) ≈ [1.4142, deg2rad(45.0), deg2rad(180.0)] rtol=1e-4
    end
end