text:
# see log/statistics.jl
function _reliable_level(ep::ErrorPropagator{T,N})::Int64 where {T,N}
    isempty(ep) && (return 1)                # results in NaN in std_error
    i = findlast(x -> x >= 32, ep.count)
    something(i, 1)
end

"""
    varN(ep::ErrorPropagator, i[, lvl])

Calculates the variance/N for the i-th argument of the error propagator at a
given binning level.
"""
function varN(ep::ErrorPropagator, i::Integer, lvl::Integer = _reliable_level(ep))
    n = ep.count[lvl]
    var(ep, i, lvl) / n
end

"""
    var(ep::ErrorPropagator, i[, lvl])

Calculates the variance for the i-th argument of the error propagator at a given
binning level.
"""
function var(
        ep::ErrorPropagator{T,N},
        i::Integer,
        lvl::Integer = _reliable_level(ep)
    ) where {N, T <: Real}

    n = ep.count[lvl]
    X = ep.sums1D[lvl][i]
    X2 = ep.sums2D[lvl][i, i]

    # lvl = 1 <=> original values
    # correct variance:
    # (∑ xᵢ^2) / (N-1) - (∑ xᵢ)(∑ xᵢ) / (N(N-1))
    X2 / (n - 1) - X^2 / (n*(n - 1))
end

function var(
        ep::ErrorPropagator{T,N},
        i::Integer,
        lvl::Integer = _reliable_level(ep)
    ) where {N, T <: Complex}

    n = ep.count[lvl]
    X = ep.sums1D[lvl][i]
    X2 = ep.sums2D[lvl][i, i]

    # lvl = 1 <=> original values
    (real(X2) + imag(X2)) / (n - 1) - (real(X)^2 + imag(X)^2) / (n*(n - 1))
end

function var(
        ep::ErrorPropagator{<: AbstractArray{T, D}, N},
        i::Integer,
        lvl::Integer = _reliable_level(ep)
    ) where {N, D, T <: Real}

    n = ep.count[lvl]
    X = ep.sums1D[lvl][i]
    X2 = ep.sums2D[lvl][i, i]

    @. X2 / (n - 1) - X^2 / (n*(n - 1))
end

function var(
        ep::ErrorPropagator{<: AbstractArray{T, D}, N},
        i::Integer,
        lvl::Integer = _reliable_level(ep)
    ) where {N, D, T <: Complex}

    n = ep.count[lvl]
    X = ep.sums1D[lvl][i]
    X2 = ep.sums2D[lvl][i, i]

    @. (real(X2) + imag(X2)) / (n - 1) - (real(X)^2 + imag(X)^2) / (n*(n - 1))
end

# NOTE works for all types of ErrorPropagators
"""
    mean(ep::ErrorPropagator, i[, lvl])

Calculates the mean for the i-th argument of the error propagator at a given
binning level.
"""
function mean(ep::ErrorPropagator, i::Integer, lvl::Integer = 1)
    ep.sums1D[lvl][i] / ep.count[lvl]
end

"""
    tau(ep::ErrorPropagator, i[, lvl])

Calculates the autocorrelation time tau for the i-th argument of the error
propagator at a given binning level.
"""
function tau(
        ep::ErrorPropagator{T,N},
        i::Integer,
        lvl::Integer = _reliable_level(ep)
    ) where {N , T <: Number}

    var_0 = varN(ep, i, 1)
    var_l = varN(ep, i, lvl)
    0.5 * (var_l / var_0 - 1)
end
function tau(
        ep::ErrorPropagator{T,N},
        i::Integer,
        lvl::Integer = _reliable_level(ep)
    ) where {N , T <: AbstractArray}

    var_0 = varN(ep, i, 1)
    var_l = varN(ep, i, lvl)
    @. 0.5 * (var_l / var_0 - 1)
end

"""
    std_error(ep::ErrorPropagator, i[, lvl])

Calculates the standard error of the mean for the i-th argument of the error
propagator at a given binning level.
"""
function std_error(
        ep::ErrorPropagator{T,N},
        i::Integer,
        lvl::Integer=_reliable_level(ep)
    ) where {N, T <: Number}

    sqrt(varN(ep, i, lvl))
end
function std_error(
        ep::ErrorPropagator{T,N},
        i::Integer,
        lvl::Integer=_reliable_level(ep)
    ) where {N, T <: AbstractArray}

    sqrt.(varN(ep, i, lvl))
end

# Generated functions
for name in [:varN, :var, :tau, :std_error]
    # generates functions fs(ep[, lvl]) = [f(ep, 1, lvl), .., f(ep, N_args, lvl)]
    @eval begin
        function $(Symbol(name, :s))(ep::ErrorPropagator, lvl=_reliable_level(ep))
            [$name(ep, i, lvl) for i in eachindex(ep.sums1D[1])]
        end
    end

    # generates functions
    # all_fs(ep) = [[f(ep, 1, lvl), .., f(ep, N_args, lvl)] for lvl in eachlvl]
    @eval begin
        function $(Symbol(:all_, name, :s))(ep::ErrorPropagator{T, N}) where {T, N}
            [$(Symbol(name, :s))(ep, lvl) for lvl in 1:N if ep.count[lvl] > 1]
        end
    end
end

# These should not default to _reliable_level, so keep them out of the
# function generation above
function means(ep::ErrorPropagator, lvl=1)
    [mean(ep, i, lvl) for i in eachindex(ep.sums1D[1])]
end
function all_means(ep::ErrorPropagator{T, N}) where {T, N}
    [means(ep, lvl) for lvl in 1:N if ep.count[lvl] > 1]
end

# Docs
@doc """
    varNs(ep::ErrorPropagator[, lvl])

Calculates the variance/N for each argument of the error propagator at a given
binning level.
""" varNs
@doc """
    vars(ep::ErrorPropagator[, lvl])

Calculates the variance for each argument of the error propagator at a given
binning level.
""" vars
@doc """
    means(ep::ErrorPropagator[, lvl])

Calculates the mean for each argument of the error propagator at a given
binning level.
""" means
@doc """
    std_errors(ep::ErrorPropagator[, lvl])

Calculates the standard error of the mean for each argument of the error
propagator at a given binning level.
""" std_errors

@doc """
    all_varNs(ep::ErrorPropagator)

Calculates the variance/N for each argument and binning level of the error
propagator. The result is indexed as `all_varNs(ep)[lvl][arg_idx]`.
""" all_varNs
@doc """
    all_vars(ep::ErrorPropagator)

Calculates the variance for each argument and binning level of the error
propagator. The result is indexed as `all_vars(ep)[lvl][arg_idx]`.
""" all_vars
@doc """
    all_means(ep::ErrorPropagator)

Calculates the mean for each argument and binning level of the error propagator.
The result is indexed as `all_means(ep)[lvl][arg_idx]`.
""" all_means
@doc """
    all_std_errors(ep::ErrorPropagator)

Calculates the standard error of the mean for each argument and binning level of
the error propagator. The result is indexed as `all_varNs(ep)[lvl][arg_idx]`.
""" all_std_errors

# Special Error Propagator functions

"""
    covmat(ep::ErrorPropagator[, lvl])

Returns the covariance matrix for a given level of the error propgator.
"""
function covmat(
        ep::ErrorPropagator{T, N}, lvl = _reliable_level(ep)
    ) where {N , T <: Number}

    invN = 1.0 / ep.count[lvl]
    invN1 = 1.0 / (ep.count[lvl] - 1)
    [
        (
            ep.sums2D[lvl][i, j] -
            ep.sums1D[lvl][i] * conj(ep.sums1D[lvl][j]) * invN
        ) * invN1
        for i in eachindex(ep.sums1D[lvl]), j in eachindex(ep.sums1D[lvl])
    ]
end
function covmat(
        ep::ErrorPropagator{T, N}, lvl = _reliable_level(ep)
    ) where {N , T <: AbstractArray}

    invN = 1.0 / ep.count[lvl]
    invN1 = 1.0 / (ep.count[lvl] - 1)
    [
        (
            ep.sums2D[lvl][i, j] -
            ep.sums1D[lvl][i] .* conj(ep.sums1D[lvl][j]) * invN
        ) * invN1
        for i in eachindex(ep.sums1D[lvl]), j in eachindex(ep.sums1D[lvl])
    ]
end

"""
    var(ep::ErrorPropagator, gradient[, lvl])

Gives the first-order variance estimate of a function `f` acting on the
arguments of the error propagator. `gradient` is either the gradient of `f` (a
function) or a vector `∇f(means(ep))`. To get an estimate mean value of `f`,
`mean(ep, f)` can be used.
"""
function var(
        ep::ErrorPropagator{T, N},
        gradient::Vector,
        lvl = _reliable_level(ep)
    ) where {T <: Real, N}

    result = 0.0
    invN = 1.0 / ep.count[lvl]
    # invNN1 = 1.0 / (ep.count[lvl] * (ep.count[lvl] - 1))
    invN1 = 1.0 / (ep.count[lvl] - 1)
    for i in eachindex(ep.sums1D[lvl])
        for j in eachindex(ep.sums1D[lvl])
            result += gradient[i] * gradient[j] * (
                ep.sums2D[lvl][i, j] -
                ep.sums1D[lvl][i] * ep.sums1D[lvl][j] * invN
            ) * invN1
        end
    end

    result
end
function var(
        ep::ErrorPropagator{T, N},
        gradient::Vector,
        lvl = _reliable_level(ep)
    ) where {T <: Complex, N}

    result = 0.0 * 0.0im
    invN = 1.0 / ep.count[lvl]
    # invNN1 = 1.0 / (ep.count[lvl] * (ep.count[lvl] - 1))
    invN1 = 1.0 / (ep.count[lvl] - 1)
    for i in eachindex(ep.sums1D[lvl])
        for j in eachindex(ep.sums1D[lvl])
            result += gradient[i] * conj(gradient[j]) * (
                ep.sums2D[lvl][i, j] -
                ep.sums1D[lvl][i] * conj(ep.sums1D[lvl][j]) * invN
            ) * invN1
        end
    end

    abs(result)
end

# Wrappers
function var(ep::ErrorPropagator, gradient::Function, lvl = _reliable_level(ep))
    grad = gradient(means(ep))
    if typeof(grad) <: AbstractArray
        return var(ep, grad, lvl)
    else
        return var(ep, [grad], lvl)
    end
end

"""
    std_error(ep::ErrorPropagator, gradient[, lvl])

Gives the first-order standard error estimate of a function `f` acting on the
arguments of the error propagator. `gradient` is either the gradient of `f` (a
function) or a vector `∇f(means(ep))`. To get an estimate mean value of `f`,
`mean(ep, f)` can be used.
"""
function std_error(ep::ErrorPropagator, gradient, lvl = _reliable_level(ep))
    sqrt(var(ep, gradient, lvl) / ep.count[lvl])
end

"""
    mean(ep, f[, lvl=1])

Returns an estimate for the mean value of `f`, where `f` is a function acting
on the sample pushed the error porpagator. `f` must be of the form `f(v)`, where
`v = [mean_arg1, mean_arg2, ..., mean_argN]` is a vector containing the averages
of each argument pushed to the error propagator.
"""
function mean(ep::ErrorPropagator, f::Function, lvl = 1)
    f(means(ep, lvl))
end


text:
using OrdinaryDiffEq, StaticArrays, Test

# More ODE event tests, cf. #201, #199, #198, #197
function test_callback_inplace(alg; kwargs...)
    f = (du, u, p, t) -> @. du = u
    cb = ContinuousCallback((u,t,int) -> u[1] - exp(1), terminate!)
    prob = ODEProblem(f, [5.0], (2.0, 0.0), callback=cb)
    sol = solve(prob, alg; kwargs...)
    sol.u[end][1] ≈ exp(1)
end

function test_callback_outofplace(alg; kwargs...)
    f = (u, p, t) -> copy(u)
    cb = ContinuousCallback((u,t,int) -> u[1] - exp(1), terminate!)
    prob = ODEProblem(f, [5.0], (2.0, 0.0), callback=cb)
    sol = solve(prob, alg; kwargs...)
    sol.u[end][1] ≈ exp(1)
end

function test_callback_scalar(alg; kwargs...)
    f = (u, p, t) -> u
    cb = ContinuousCallback((u,t,int) -> u - exp(1), terminate!)
    prob = ODEProblem(f, 5.0, (2.0, 0.0), callback=cb)
    sol = solve(prob, alg; kwargs...)
    sol.u[end] ≈ exp(1)
end

function test_callback_svector(alg; kwargs...)
    f = (u, p, t) -> u
    cb = ContinuousCallback((u,t,int) -> u[1] - exp(1), terminate!)
    prob = ODEProblem(f, SVector(5.0), (2.0, 0.0), callback=cb)
    sol = solve(prob, alg; kwargs...)
    sol.u[end][1] ≈ exp(1)
end

function test_callback_mvector(alg; kwargs...)
    f = (u, p, t) -> copy(u)
    cb = ContinuousCallback((u,t,int) -> u[1] - exp(1), terminate!)
    prob = ODEProblem(f, MVector(5.0), (2.0, 0.0), callback=cb)
    sol = solve(prob, alg; kwargs...)
    sol.u[end][1] ≈ exp(1)
end

println("inplace")
@test test_callback_inplace(BS3())
@test test_callback_inplace(BS5())
@test test_callback_inplace(SSPRK43())
@test test_callback_inplace(SSPRK432())
@test test_callback_inplace(SSPRK932())
@test test_callback_inplace(OwrenZen3())
@test test_callback_inplace(OwrenZen4())
@test test_callback_inplace(OwrenZen5())
@test test_callback_inplace(DP5())
@test test_callback_inplace(DP8())
@test test_callback_inplace(Feagin10())
@test test_callback_inplace(Feagin12())
@test test_callback_inplace(Feagin14())
@test test_callback_inplace(TanYam7())
@test test_callback_inplace(Tsit5())
@test test_callback_inplace(TsitPap8())
@test test_callback_inplace(Vern6())
@test test_callback_inplace(Vern7())
@test test_callback_inplace(Vern8())
@test test_callback_inplace(Vern9())
@test test_callback_inplace(Rosenbrock23())
@test test_callback_inplace(Rosenbrock32())
@test test_callback_inplace(AB3(); dt=0.1)
@test test_callback_inplace(ABM32(); dt=0.1)
@test test_callback_inplace(AB4(); dt=0.1)
@test test_callback_inplace(ABM43(); dt=0.1)
@test test_callback_inplace(AB5(); dt=0.1)
@test test_callback_inplace(ABM54(); dt=0.1)
@test test_callback_inplace(VCAB3())
@test test_callback_inplace(VCAB4())
@test test_callback_inplace(VCAB5())
@test test_callback_inplace(VCABM3())
@test test_callback_inplace(VCABM4())
@test test_callback_inplace(VCABM5())

println("outofplace")
@test test_callback_outofplace(BS3())
@test test_callback_outofplace(BS5())
@test test_callback_outofplace(SSPRK43())
@test test_callback_outofplace(SSPRK432())
@test test_callback_outofplace(SSPRK932())
@test test_callback_outofplace(OwrenZen3())
@test test_callback_outofplace(OwrenZen4())
@test test_callback_outofplace(OwrenZen5())
@test test_callback_outofplace(DP5())
@test test_callback_outofplace(DP8())
@test test_callback_outofplace(Feagin10())
@test test_callback_outofplace(Feagin12())
@test test_callback_outofplace(Feagin14())
@test test_callback_outofplace(TanYam7())
@test test_callback_outofplace(Tsit5())
@test test_callback_outofplace(TsitPap8())
@test test_callback_outofplace(Vern6())
@test test_callback_outofplace(Vern7())
@test test_callback_outofplace(Vern8())
@test test_callback_outofplace(Vern9())
@test test_callback_outofplace(Rosenbrock23())
@test test_callback_outofplace(Rosenbrock32())
@test test_callback_outofplace(AB3(); dt=0.1)
@test test_callback_outofplace(ABM32(); dt=0.1)
@test test_callback_outofplace(AB4(); dt=0.1)
@test test_callback_outofplace(ABM43(); dt=0.1)
@test test_callback_outofplace(AB5(); dt=0.1)
@test test_callback_outofplace(ABM54(); dt=0.1)
@test test_callback_outofplace(VCAB3())
@test test_callback_outofplace(VCAB4())
@test test_callback_outofplace(VCAB5())
@test test_callback_outofplace(VCABM3())
@test test_callback_outofplace(VCABM4())
@test test_callback_outofplace(VCABM5())

println("scalar")
@test test_callback_scalar(BS3())
@test test_callback_scalar(BS5())
@test test_callback_scalar(SSPRK43())
@test test_callback_scalar(SSPRK432())
@test test_callback_scalar(SSPRK932())
@test test_callback_scalar(OwrenZen3())
@test test_callback_scalar(OwrenZen4())
@test test_callback_scalar(OwrenZen5())
@test test_callback_scalar(DP5())
@test test_callback_scalar(DP8())
@test test_callback_scalar(Feagin10())
@test test_callback_scalar(Feagin12())
@test test_callback_scalar(Feagin14())
@test test_callback_scalar(TanYam7())
@test test_callback_scalar(Tsit5())
@test test_callback_scalar(TsitPap8())
@test test_callback_scalar(Vern6())
@test test_callback_scalar(Vern7())
@test test_callback_scalar(Vern8())
@test test_callback_scalar(Vern9())
@test test_callback_scalar(Rosenbrock23())
@test test_callback_scalar(Rosenbrock32())
@test test_callback_scalar(AB3(); dt=0.1)
@test test_callback_scalar(ABM32(); dt=0.1)
@test test_callback_scalar(AB4(); dt=0.1)
@test test_callback_scalar(ABM43(); dt=0.1)
@test test_callback_scalar(AB5(); dt=0.1)
@test test_callback_scalar(ABM54(); dt=0.1)
@test test_callback_scalar(VCAB3())
@test test_callback_scalar(VCAB4())
@test test_callback_scalar(VCAB5())
@test test_callback_scalar(VCABM3())
@test test_callback_scalar(VCABM4())
@test test_callback_scalar(VCABM5())

println("svector")
@test test_callback_svector(BS3())
@test test_callback_svector(BS5())
@test test_callback_svector(SSPRK43())
@test test_callback_svector(SSPRK432())
@test test_callback_svector(SSPRK932())
@test test_callback_svector(OwrenZen3())
@test test_callback_svector(OwrenZen4())
@test test_callback_svector(OwrenZen5())
@test test_callback_svector(DP5())
@test test_callback_svector(DP8())
@test test_callback_svector(Feagin10())
@test test_callback_svector(Feagin12())
@test test_callback_svector(Feagin14())
@test test_callback_svector(TanYam7())
@test test_callback_svector(Tsit5())
@test test_callback_svector(TsitPap8())
@test test_callback_svector(Vern6())
@test test_callback_svector(Vern7())
@test test_callback_svector(Vern8())
@test test_callback_svector(Vern9())
@test test_callback_svector(Rosenbrock23())
@test test_callback_svector(Rosenbrock32())
@test test_callback_svector(AB3(); dt=0.1)
@test test_callback_svector(ABM32(); dt=0.1)
@test test_callback_svector(AB4(); dt=0.1)
@test test_callback_svector(ABM43(); dt=0.1)
@test test_callback_svector(AB5(); dt=0.1)
@test test_callback_svector(ABM54(); dt=0.1)
@test test_callback_svector(VCAB3())
@test test_callback_svector(VCAB4())
@test test_callback_svector(VCAB5())
@test test_callback_svector(VCABM3())
@test test_callback_svector(VCABM4())
@test test_callback_svector(VCABM5())

print("mvector")
@test test_callback_mvector(BS3())
@test test_callback_mvector(BS5())
print(".")
@test test_callback_mvector(SSPRK43())
@test test_callback_mvector(SSPRK432())
@test test_callback_mvector(SSPRK932())
print(".")
@test test_callback_mvector(OwrenZen3())
@test test_callback_mvector(OwrenZen4())
print(".")
@test test_callback_mvector(OwrenZen5())
@test test_callback_mvector(DP5())
print(".")
@test test_callback_mvector(DP8())
@test test_callback_mvector(Feagin10())
print(".")
@test test_callback_mvector(Feagin12())
@test test_callback_mvector(Feagin14())
print(".")
@test test_callback_mvector(TanYam7())
@test test_callback_mvector(Tsit5())
print(".")
@test test_callback_mvector(TsitPap8())
@test test_callback_mvector(Vern6())
print(".")
@test test_callback_mvector(Vern7())
@test test_callback_mvector(Vern8())
print(".")
@test test_callback_mvector(Vern9())
@test test_callback_mvector(Rosenbrock23())
print(".")
@test test_callback_mvector(Rosenbrock32())
@test test_callback_mvector(AB3(); dt=0.1)
print(".")
@test test_callback_mvector(ABM32(); dt=0.1)
@test test_callback_mvector(AB4(); dt=0.1)
print(".")
@test test_callback_mvector(ABM43(); dt=0.1)
@test test_callback_mvector(AB5(); dt=0.1)
print(".")
@test test_callback_mvector(ABM54(); dt=0.1)
@test test_callback_mvector(VCAB3())
print(".")
@test test_callback_mvector(VCAB4())
@test test_callback_mvector(VCAB5())
print(".")
@test test_callback_mvector(VCABM3())
@test test_callback_mvector(VCABM4())
print(".")
@test test_callback_mvector(VCABM5())
println(".")

ext: jl


text:
functions = Dict(
    "atan2" => x -> "atan($x)",
    "add" => x -> "+($x)",
    "sub" => x -> "-($x)",
    "pos" => x -> "+($x)",
    "neg" => x -> "-($x)",
    "mul" => x -> "*($x)",
    "div" => x -> "/($x)",
    "convexHull" => x -> "hull($x)",
    "intersection" => x -> "intersect($x)",
    "interior" => x -> "isinterior($x)",
    "subset" => x -> "issubset($x)",
    "equal" => x -> "==($x)",
    "sqr" => x -> "$x^2",
    "pow" => x -> "^($x)",
    "rootn" => x -> "nthroot($x)",
    "sqrt" => x -> "sqrt($x)",
    "exp" => x -> "exp($x)",
    "exp2" => x -> "exp2($x)",
    "exp10" => x -> "exp10($x)",
    "log" => x -> "log($x)",
    "log2" => x -> "log2($x)",
    "log10" => x -> "log10($x)",
    "sin" => x -> "sin($x)",
    "cos" => x -> "cos($x)",
    "tan" => x -> "tan($x)",
    "cot" => x -> "cot($x)",
    "asin" => x -> "asin($x)",
    "acos" => x -> "acos($x)",
    "atan" => x -> "atan($x)",
    "acot" => x -> "acot($x)",
    "sinh" => x -> "sinh($x)",
    "cosh" => x -> "cosh($x)",
    "tanh" => x -> "tanh($x)",
    "coth" => x -> "coth($x)",
    "asinh" => x -> "asinh($x)",
    "acosh" => x -> "acosh($x)",
    "atanh" => x -> "atanh($x)",
    "acoth" => x -> "acoth($x)",
    "expm1" => x -> "expm1($x)",
    "logp1" => x -> "log1p($x)",
    "cbrt" => x -> "cbrt($x)",
    "csc" => x -> "csc($x)",
    "csch" => x -> "csch($x)",
    "hypot" => x -> "hypot($x)",
    "sec" => x -> "sec($x)",
    "sech" => x -> "sech($x)",
    "intervalPart" => x -> "interval_part($x)",
    "isEmpty" => x -> "isempty($x)",
    "isEntire" => x -> "isentire($x)",
    "isNaI" => x -> "isnai($x)",
    "less" => x -> "≤($x)",
    "strictLess" => x -> "<($x)",
    "precedes" => x -> "precedes($x)",
    "strictPrecedes" => x -> "strictprecedes($x)",
    "disjoint" => x -> "isdisjoint($x)",
    "newDec" => x -> "DecoratedInterval($x)",
    "setDec" => x -> "DecoratedInterval($x)",
    "decorationPart" => x -> "decoration($x)",
    "recip" => x -> "inv($x)",
    "fma" => x -> "fma($x)",
    "pown" => x -> "^($x)",
    "sign" => x -> "sign($x)",
    "ceil" => x -> "ceil($x)",
    "floor" => x -> "floor($x)",
    "trunc" => x -> "trunc($x)",
    "roundTiesToEven" => x -> "round($x)",
    "roundTiesToAway" => x -> "round($x, RoundNearestTiesAway)",
    "abs" => x -> "abs($x)",
    "min" => x -> "min($x)",
    "max" => x -> "max($x)",
    "inf" => x -> "inf($x)",
    "sup" => x -> "sup($x)",
    "mid" => x -> "mid($x)",
    "rad" => x -> "radius($x)",
    "midRad" => x -> "midpoint_radius($x)",
    "wid" => x -> "diam($x)",
    "mag" => x -> "mag($x)",
    "mig" => x -> "mig($x)",
    "overlap" => x -> "MISSING_overlap($x)",
    "isCommonInterval" => x -> "iscommon($x)",
    "isSingleton" => x -> "isthin($x)",
    "isMember" => x -> "in($x)",
    "cancelPlus" => x -> "cancelplus($x)",
    "cancelMinus" => x -> "cancelminus($x)",
    "sum_nearest" => x -> "sum($x)",
    "dot_nearest" => x -> "dot($x)",
    "sum_abs_nearest" => x -> "sum(abs.($x))",
    "sum_sqr_nearest" => x -> "sum($x.^2)",
    "sqrRev" => x -> "sqr_rev($x)[2]",
    "sqrRevBin" => x -> "sqr_rev($x)[2]",
    "absRev" => x -> "abs_rev($x)[2]",
    "absRevBin" => x -> "abs_rev($x)[2]",
    "pownRev" => x -> "power_rev($x)[2]",
    "pownRevBin" => x -> "power_rev($x)[2]",
    "sinRev" => x -> "sin_rev($x)[2]",
    "sinRevBin" => x -> "sin_rev($x)[2]",
    "cosRev" => x -> "cos_rev($x)[2]",
    "cosRevBin" => x -> "cos_rev($x)[2]",
    "tanRev" => x -> "tan_rev($x)[2]",
    "tanRevBin" => x -> "tan_rev($x)[2]",
    "coshRev" => x -> "cosh_rev($x)[2]",
    "coshRevBin" => x -> "cosh_rev($x)[2]",
    "mulRev" => x -> "mul_rev_IEEE1788($x)",
    "mulRevTen" => x -> "mul_rev_IEEE1788($x)",
    "mulRevToPair" => x -> "mul_rev_to_pair($x)",
    "powRev1" => x -> "pow_rev1($x)",
    "powRev2" => x -> "pow_rev2($x)"
)


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

@stochastic_model farmer begin
    @stage 1 begin
        @parameters begin
            Crops = [:wheat, :corn, :beets]
            Cost = Dict(:wheat=>150, :corn=>230, :beets=>260)
            Budget = 500
        end
        @decision(farmer, x[c in Crops] >= 0)
        @objective(farmer, Min, sum(Cost[c]*x[c] for c in Crops))
        @constraint(farmer, sum(x[c] for c in Crops) <= Budget)
    end
    @stage 2 begin
        @parameters begin
            Crops = [:wheat, :corn, :beets]
            Required = Dict(:wheat=>200, :corn=>240, :beets=>0)
            PurchasePrice = Dict(:wheat=>238, :corn=>210)
            SellPrice = Dict(:wheat=>170, :corn=>150, :beets=>36, :extra_beets=>10)
        end
        @uncertain ξ[c in Crops]
        @recourse(farmer, y[p in setdiff(Crops, [:beets])] >= 0)
        @recourse(farmer, w[s in Crops ∪ [:extra_beets]] >= 0)
        @objective(farmer, Min, sum(PurchasePrice[p] * y[p] for p in setdiff(Crops, [:beets]))
                   - sum(SellPrice[s] * w[s] for s in Crops ∪ [:extra_beets]))
        @constraint(farmer, minimum_requirement[p in setdiff(Crops, [:beets])],
            ξ[p] * x[p] + y[p] - w[p] >= Required[p])
        @constraint(farmer, minimum_requirement_beets,
            ξ[:beets] * x[:beets] - w[:beets] - w[:extra_beets] >= Required[:beets])
        @constraint(farmer, beets_quota, w[:beets] <= 6000)
    end
end

Crops = [:wheat, :corn, :beets]
ξ₁ = @scenario ξ[c in Crops] = [3.0, 3.6, 24.0] probability = 1/3
ξ₂ = @scenario ξ[c in Crops] = [2.5, 3.0, 20.0] probability = 1/3
ξ₃ = @scenario ξ[c in Crops] = [2.0, 2.4, 16.0] probability = 1/3

farmer_res = SPResult([170,80,250], Dict(1 => [0., 0., 310., 48, 6000, 0.], 2 => [0., 0., 225., 0., 5000, 0.], 3 => [0., 48., 140., 0., 4000, 0.]), -108390, -115405.56, 7015.56, 1150, -118600, -107240)
push!(problems, (farmer, [ξ₁,ξ₂,ξ₃], farmer_res,"Farmer"))


text:
using ConstrainedDynamics
using Test

@testset "ConstrainedDynamics.jl" begin
    # Write your own tests here.
end


text:
using Agents
using AgentsX
using Random
using InteractiveDynamics
using CairoMakie

include("LayeredAgent.jl")
include("LayeredAgent2.jl")

function initialise(step_function; numagents=100, seed=250)
    space = GridSpace((20, 20); periodic=false)
    properties = Dict(:step_function => step_function)
    rng = Random.MersenneTwister(seed)

    model = AgentsX.ABM(
        [LayeredAgent, LayeredAgent2], space;
        properties, rng
    )

    for n in 1:numagents/2
        agent = LayeredAgent(n, (1, 1),
                Social(rand(0:1)), Economic(Bool(rand(0:1))),
                Cognitive(Bool(rand(0:1)))
        )
        agent2 = LayeredAgent2(n*2, (1, 1),
                Economic(Bool(rand(0:1))),
                Cognitive(Bool(rand(0:1))),
                Physiological(Bool(rand(0:1)))
        )
        add_agent_single!(agent, model)
        add_agent_single!(agent2, model)
    end

    return model
end

model = initialise([move, think])

step!(model, agent_step!)


text:
module Laboratorio3

include("./sumBetweenElements.jl")
include("./freqTwo.jl")
include("./sumOfMult.jl")
include("./countInc.jl")
include("./countAAA.jl")

export sumBetweenElements, freqTwo, sumOfMult, countInc, countAAA

end #module


text:
## Element type
struct IIDMVolt <: AbstractNodeLabel
  busname::String
  busid::Int
  substationname::String
  nomV::Float64    # kV
  minV::Float64    # p.u.
  maxV::Float64    # p.u.
end

## 1. Variables creation
function create_vars!(element::T, bus::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}, scenario::String) where T <: IIDMVolt
  bus_vars[elemid] = Variable(variable_name("VOLT", bus, "", scenario), Complex)
  return
end

# function create_vars!(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}) where T <: MyType
#   return
# end

## 2. Power balance
# function Snodal(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}) where T <: MyType
#   return ["", Polynomial(), 0, 0]
# end

# function Sorig(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}) where T<:MyType
#   return Polynomial()
# end

# function Sdest(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}) where T<:MyType
#   return Polynomial()
# end

## 3. Constraints creation
function constraint(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}) where T <: IIDMVolt
  cstrname = "VOLTM_$busid"
  Vmin = element.minV
  Vmax = element.maxV
  nomV = element.nomV
  p = abs2(bus_vars["Volt"])
  return Dict{String, Constraint}(cstrname=> Vmin^2 << p << Vmax^2)
end

# function constraint(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}) where T <: MyType
#   return [[], [], [], []]
# end

# ## 4. Measure accessor
# function add_measure(element::T, busid::String, elemid::String, elem_formulation::Symbol, bus_vars::Dict{String, Variable}) where T <: IIDMVolt
#   return Dict("V_$(busid)" => abs2((1/element.nomV)*bus_vars["Volt"]) ==  element.V)
# end
#
# # function add_measure(element::T, link::Link, elemid::String, elem_formulation::Symbol, link_vars::Dict{String, Variable}) where T
# #   return Dict{String, Constraint}()
# # end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
struct PSO{T,S,fType} <: Optimizer

    # Optimization problem
    prob::Problem{fType,S}

    # Swarm of particles
    swarm::Swarm{T}

    # Settings
    initMethod::Symbol      # Only option now is :Uniform
    updateMethod::Symbol    # Only option now is :MATLAB

    # PSO specific parameters/options
    inertiaRange::Tuple{T,T}
    minNeighborFrac::T
    selfAdjustWeight::T
    socialAdjustWeight::T

    # Results data
    results::Results{T}

    function PSO{T}(prob::Problem{fType,S}, numParticles::Integer, initMethod::Symbol,
        updateMethod::Symbol, inertiaRange::Tuple{T,T}, minNeighborFrac::T,
        selfAdjustWeight::T, socialAdjustWeight::T) where {T,S,fType}

        # Initialize Swarm
        N = length(prob.LB)
        swarm = Swarm{T}(N, numParticles)

        # Initialize results
        results = Results{T}(undef, N)

        # Check that init and update methods exist
        if initMethod != :Uniform && initMethod != :LogisticsMap
            throw(ArgumentError("$initMethod is not a PSO initialization method."))
        end
        if updateMethod != :MATLAB
            throw(ArgumentError("$updateMethod is not a PSO update method."))
        end

        return new{T,S,fType}(prob, swarm, initMethod, updateMethod,
                        inertiaRange, minNeighborFrac, selfAdjustWeight,
                        socialAdjustWeight, results)
    end
end

# ===== Constructors

function PSO(prob::Problem{fType,S}; numParticles = 100, inertiaRange = (0.1, 1.1), minNeighborFrac = 0.25,
    selfAdjustWeight = 1.49, socialAdjustWeight = 1.49, initMethod::Symbol = :Uniform,
    updateMethod::Symbol = :MATLAB) where {S,fType}

    # Error checking
    length(inertiaRange) == 2 || throw(ArgumentError("inertiaRange must be of length 2."))
    minNeighborFrac > 0       || throw(ArgumentError("minNeighborFrac must be > 0."))

    # Type info
    T = typeof(1.0)
    nIRange = (T(inertiaRange[1]), T(inertiaRange[2]))

    # Call constructor
    return PSO{T}(prob, numParticles, initMethod, updateMethod, nIRange, T(minNeighborFrac),
        T(selfAdjustWeight), T(socialAdjustWeight))
end

# ===== Methods

function _optimize!(pso::PSO, opts::Options)

    # Initialize PSO
    initialize!(pso, opts)

    # Iterate PSO
    iterate!(pso, opts)

    # Display results and return
    if opts.display
        display(pso.results)
    end

    return pso.results
end

function initialize!(pso::PSO, opts::Options)

    # Initialize position and velocities
    if pso.initMethod == :Uniform
        uniformInitialization!(pso.swarm, pso.prob, opts)
    elseif pso.initMethod == :LogisticsMap
        logisticsMapInitialization!(pso.swarm, pso.prob, opts)
    else
        throw(ArgumentError("PSO initialization method is not implemented."))
    end

    # Evaluate Objective Function
    feval!(pso.swarm, pso.prob.f, opts; init = true)

    # Set swarm global best obj. func. value to Inf
    pso.swarm.b = Inf

    # Set global best
    setGlobalBest!(pso.swarm)

    # Initialize neighborhood size
    pso.swarm.n = max(2, floor(length(pso.swarm)*pso.minNeighborFrac))

    # Initial inertia
    if pso.inertiaRange[2] > 0
        pso.swarm.w = (pso.inertiaRange[2] > pso.inertiaRange[1]) ?
                            pso.inertiaRange[2] : pso.inertiaRange[1]
    else
        pso.swarm.w = (pso.inertiaRange[2] < pso.inertiaRange[1]) ?
                            pso.inertiaRange[2] : pso.inertiaRange[1]
    end

    # Initialize self and social adjustment
    pso.swarm.y₁ = pso.selfAdjustWeight
    pso.swarm.y₂ = pso.socialAdjustWeight

    # Print Status
    if opts.display
        printStatus(pso.swarm, 0.0, 0, 0)
    end

    return nothing
end

function iterate!(pso::PSO, opts::Options)

    # Initialize time and iteration counter
    t0  = time()
    stallT0 = t0
    iters = 0
    stallIters = 0
    fStall = Inf

    # Compute minimum neighborhood size
    minNeighborSize = max(2, floor(length(pso.swarm)*pso.minNeighborFrac))

    # Begin loop
    exitFlag = 0
    while exitFlag == 0

        # Update iteration counter
        iters += 1

        # Update swarm velocities
        updateVelocities!(pso.swarm)

        # Update positions
        step!(pso.swarm)

        # Enforce Bounds
        if !(all(isinf.(pso.prob.LB)) && all(isinf.(pso.prob.UB)))
            enforceBounds!(pso.swarm, pso.prob.LB, pso.prob.UB)
        end

        # Evaluate objective function
        feval!(pso.swarm, pso.prob.f, opts)

        # Update global best
        flag::Bool = setGlobalBest!(pso.swarm)

        # Update stall counter and neighborhood
        if flag
            pso.swarm.c = max(0, pso.swarm.c - 1)
            pso.swarm.n = minNeighborSize
        else
            pso.swarm.c += 1
            pso.swarm.n = min(pso.swarm.n + minNeighborSize,
                                length(pso.swarm) - 1)
        end

        # Update inertia
        if pso.swarm.c < 2
            pso.swarm.w *= 2.0
        elseif pso.swarm.c > 5
            pso.swarm.w /= 2.0
        end

        # Ensure inertia is in inertia range
        if  pso.swarm.w > pso.inertiaRange[2]
            pso.swarm.w = pso.inertiaRange[2]
        elseif pso.swarm.w < pso.inertiaRange[1]
            pso.swarm.w = pso.inertiaRange[1]
        end

        # Track stalling
        if fStall - pso.swarm.b > opts.funcTol
            fStall = pso.swarm.b
            stallIters = 0
            stallT0 = time()
        else
            stallIters += 1
        end

        # Stopping criteria
        if stallIters >= opts.maxStallIters
            exitFlag = 1
        elseif iters >= opts.maxIters
            exitFlag = 2
        elseif pso.swarm.b <= opts.objLimit
            exitFlag = 3
        elseif time() - stallT0 >= opts.maxStallTime
            exitFlag = 4
        elseif time() - t0 >= opts.maxTime
            exitFlag = 5
        end

        # Output Status
        if opts.display && iters % opts.displayInterval == 0
            printStatus(pso.swarm, time() - t0, iters, stallIters)
        end
    end

    # Set results
    setResults!(pso, iters, time() - t0, exitFlag)

    return nothing
end

function setResults!(pso::PSO, iters::Int, Δt::AbstractFloat, exitFlag::Int)
    pso.results.fbest = pso.swarm.b
    pso.results.xbest .= pso.swarm.d
    pso.results.iters = iters
    pso.results.time = Δt
    pso.results.exitFlag = exitFlag

    return nothing
end


text:
using Plots

using simulator.fresnelltools
using simulator.matrixcore
using simulator.materials.spesifics

save = false

n_1 = SiO2_core_Sellmeier
n_2 = Au
n_3 = Pd
n_4 = l -> 1

d2 = 10e-9
d2s = 0e-9:1e-9:40e-9
d3s = 0e-9:1e-9:60e-9

λs = 900e-9:1e-9:1250e-9
λs_broadband = 400e-9:1e-9:1800e-9
θd = 45
θds = 5:1:80

description = "Palladium on $(round(d2 * 1e9, digits=2)) nm gold"

fig_path = "../rapport/figures/thickness_example/"
fig_token = replace(description, " " => "_")

if save
    gr()
else
    plotly()
end
θ1 = θd * π / 180
θs = θds .* (π / 180)

function system_single(λ::Float64, θ1::Float64, d2::Float64, d3::Float64)
    i_1_2_te, i_1_2_tm, θ2 = FresnellBoundrary(n_1(λ), n_2(λ), θ1)
    bulk_2 = FresnellSlab(n_2(λ), 2*π/λ, d2, θ2)
    i_2_3_te, i_2_3_tm, θ3 = FresnellBoundrary(n_2(λ), n_3(λ), θ2)
    bulk_3 = FresnellSlab(n_3(λ), 2*π/λ, d3, θ3)
    i_3_4_te, i_3_4_tm, _ = FresnellBoundrary(n_3(λ), n_4(λ), θ3)

    system_te = CascadeScattering([i_1_2_te, bulk_2, i_2_3_te, bulk_3, i_3_4_te])
    system_tm = CascadeScattering([i_1_2_tm, bulk_2, i_2_3_tm, bulk_3, i_3_4_tm])

    te_r = abs(system_te[2, 1])^2
    te_t = abs(system_te[1, 1])^2
    tm_r = abs(system_tm[2, 1])^2
    tm_t = abs(system_tm[1, 1])^2

    [te_r, tm_r, te_t, tm_t]
end

function get_min(λs, θ1::Float64, d2::Float64, d3::Float64)
    mins = [2.0, 2, 2, 2]
    λmins = [0.0, 0, 0, 0]

    for λ in λs
        spesifics = system_single(λ, θ1, d2, d3)
        for (i, min, spesific) in zip(1:4, mins, spesifics)
            if spesific <= min
                mins[i] = spesific
                λmins[i] = λ
            end
        end
    end

    mins, λmins
end

println("Settup complete, starting simulation")

# get_min_results = [get_min(λs, θ1, d2, d3_spesific) for d3_spesific in d3s]
# min_r_te = [item[1][1] for item in get_min_results]
# min_r_tm = [item[1][2] for item in get_min_results]
# min_λ_te = [item[2][1] for item in get_min_results]
# min_λ_tm = [item[2][2] for item in get_min_results]

## Plot plasmon minima for plasmon as a function of d3
# fig_r = plot(title="Minimum te and tm reflection for $description at θ = $θd degrees", legend=:topright, xlabel="thickness in nm", ylabel="minimum reflection coefficient", ticks=:native)
# plot!(fig_r, d3s .* 1e9, min_r_te, label="TE reflection")
# plot!(fig_r, d3s .* 1e9, min_r_tm, label="TM reflection")
# display(fig_r)
# if save savefig(fig_r, fig_path * "min_reflection_" * fig_token * ".pdf") end

## Plot frequency minima for plasmon as a function of d3
# fig_λ = plot(title="λ for minimum te and tm reflection for $description at θ = $θd degrees", legend=:topright, xlabel="thickness in nm", ylabel="λ at minimum reflection coefficient", ticks=:native)
# plot!(fig_λ, d3s .* 1e9, min_λ_te, label="TE reflection")
# plot!(fig_λ, d3s .* 1e9, min_λ_tm, label="TM reflection")
# display(fig_λ)
# if save savefig(fig_λ, fig_path * "lambda_min_reflection_" * fig_token * ".pdf") end

## Plot plasmon minima and λ for plasmon as a function of d3 for all d2s
# fig_r_all = plot(title="Minimum tm reflection for $description at θ = $θd degrees", legend=:topright, xlabel="thickness in nm", ylabel="minimum reflection coefficient", ticks=:native)
# fig_λ_all = plot(title="λ for minimum tm reflection for $description at θ = $θd degrees", legend=:topright, xlabel="thickness in nm", ylabel="λ at minimum reflection coefficient", ticks=:native)
# for d2_spesific in d2s
#     get_min_results_current = [get_min(λs, θ1, d2_spesific, d3_spesific) for d3_spesific in d3s]
#     min_r_tm_current = [item[1][2] for item in get_min_results_current]
#     min_λ_tm_current = [item[2][2] for item in get_min_results_current]

#     d2_label = round(d2_spesific * 1e9, digits=2)
#     plot!(fig_r_all, d3s .* 1e9, min_r_tm_current, label="d2 = $d2_label nm")
#     plot!(fig_λ_all, d3s .* 1e9, min_λ_tm_current, label="d2 = $d2_label nm")
# end
# display(fig_r_all)
# if save savefig(fig_r_all, fig_path * "min_reflection_all_" * fig_token * ".pdf") end
# display(fig_λ_all)
# if save savefig(fig_λ_all, fig_path * "lambda_min_reflection_all_" * fig_token * ".pdf") end

## Plot plasmon minima and λ for plasmon as a function of d3 for all θds
# fig_r_all = plot(title="Minimum tm reflection for $description at θ = $θd degrees", legend=:topright, xlabel="thickness in nm", ylabel="minimum reflection coefficient", ticks=:native)
# fig_λ_all = plot(title="λ for minimum tm reflection for $description at θ = $θd degrees", legend=:topright, xlabel="thickness in nm", ylabel="λ at minimum reflection coefficient", ticks=:native)
# for θ1_spesific in θs
#     get_min_results_current = [get_min(λs_broadband, θ1_spesific, d2, d3_spesific) for d3_spesific in d3s]
#     min_r_tm_current = [item[1][2] for item in get_min_results_current]
#     min_λ_tm_current = [item[2][2] for item in get_min_results_current]

#     θ_label = round(θ1_spesific * 180 / π, digits=2)
#     plot!(fig_r_all, d3s .* 1e9, min_r_tm_current, label="θ = $θ_label")
#     plot!(fig_λ_all, d3s .* 1e9, min_λ_tm_current, label="θ = $θ_label")
# end
# display(fig_r_all)
# if save savefig(fig_r_all, fig_path * "min_reflection_angles_" * fig_token * ".pdf") end
# display(fig_λ_all)
# if save savefig(fig_λ_all, fig_path * "lambda_min_reflection_angles_" * fig_token * ".pdf") end

## Plot plasmon spectral response for all thicknesses
# fig_all = plot(title="Wavelength respons for all thicknesses for $description at θ = $θd", legend=:outerright, xlabel="wavelength in nm", ylabel="reflection coefficient", ticks=:native)
# for d in d3s
#     systemdata = reduce(vcat, transpose.(system_single.(λs, θ1, d2, d)))
#     d = round(d * 1e9, digits=2)
#     plot!(fig_all, λs .* 1e9, systemdata[:, 1], label="TE, d = $d")
#     plot!(fig_all, λs .* 1e9, systemdata[:, 2], label="TM, d = $d")
# end
# display(fig_all)
# if save savefig(fig_all, fig_path * "responses_" * fig_token * ".pdf") end

println("Finished without errors")


text:
Serialization.deserialize(p::AbstractPickle, file::AbstractString) = open(file, "r") do io
  Serialization.deserialize(p,  io)
end
Serialization.deserialize(p::AbstractPickle, io::IO) = load(p, io)

loads(s; proto=DEFAULT_PROTO) = loads(Pickler(proto), s)
loads(p::AbstractPickle, s) = load(p, IOBuffer(s))

load(file::AbstractString; proto=DEFAULT_PROTO) = open(f->load(f; proto=proto), file)
load(io::IO; proto=DEFAULT_PROTO) = load(Pickler(proto), io)
load(p::AbstractPickle, file::AbstractString) = open(f->load(p, f), file)
function load(p::AbstractPickle, io::IO)
  while !eof(io)
    opcode = read(io, OpCode)
    run!(p, opcode, io)

    isequal(OpCodes.STOP)(opcode) && return pop!(p.stack)
  end
end

function run!(p::AbstractPickle, op::OpCode, io::IO)
  argf = OpCodes.argument(op)
  arg = isnothing(argf) ? nothing : argf(io)

  execute!(p, Val(op), arg)
end

for op in :(INT, LONG, LONG1, LONG4,
            STRING, BINSTRING, SHORT_BINSTRING,
            BYTEARRAY8,
            UNICODE, SHORT_BINUNICODE, BINUNICODE, BINUNICODE8,
            FLOAT, BINFLOAT).args
  @eval execute!(p::AbstractPickle, ::Val{OpCodes.$op}, arg) = push!(p.stack, arg)
end

for op in :(BINBYTES, SHORT_BINBYTES, BINBYTES8).args
  @eval execute!(p::AbstractPickle, ::Val{OpCodes.$op}, arg) = push!(p.stack, codeunits(String(arg)))
end

for op in :(BININT, BININT1, BININT2).args
  @eval execute!(p::AbstractPickle, ::Val{OpCodes.$op}, arg) = push!(p.stack, Int(arg))
end

for (op, arg) in zip(
  :(NONE, NEWTRUE, NEWFALSE, EMPTY_LIST, EMPTY_TUPLE, EMPTY_DICT, EMPTY_SET).args,
  :(nothing, true, false, [], (), Dict(), Set()).args
)
  @eval execute!(p::AbstractPickle, ::Val{OpCodes.$op}, arg) = push!(p.stack, $arg)
end

"""
  similar to `Expr` but mutable

used to hold unconstructable object. require mutable for pickle memo work.
"""
mutable struct Defer
  head::Symbol
  args::Vector{Any}

  Defer(t::Tuple) = new(t...)
  Defer(head, @nospecialize(args...)) = new(head, collect(Any, args))
end

function show(io::IO, def::Defer)
  if isempty(def.args)
    print(io, "Defer(:$(def.head))")
  else
    print(io, "Defer(:$(def.head), ")
    join(io, def.args, ", ")
    print(io, ')')
  end
end

function wrap!(def::Defer, tag, @nospecialize(args...))
  old_def = Defer((def.head, def.args))
  def.head = tag
  def.args = [old_def, args...]
  def
end

isdefer(::Defer) = true
isdefer(c::Union{Array, Tuple, Set}) = any(isdefer, c)
isdefer(d::Dict) = any(isdefer, values(d)) || any(isdefer, keys(d))
isdefer(x) = false

macro deferexecutef(op, tag, popf, pushf)
  return quote
  end
end

_setkey!(obj, rpair) = setindex!(obj, rpair...)
_setkeys!(obj, pairs) = foreach(i->setindex!(obj, pairs[i+1], pairs[i]), 1:2:length(pairs))
_addvalue(obj, values) = foreach(x->push!(obj, x), values)

for (op, tag, popf, pushf) in (
  (:APPEND,   :(:append),   pop!,    push!),
  (:APPENDS,  :(:appends),  unmark!, append!),
  (:SETITEM,  :(:setitem),  pop2!,   _setkey!),
  (:SETITEMS, :(:setitems), unmark!, _setkeys!),
  (:ADDITEMS, :(:additems), unmark!, _addvalue))
  @eval function execute!(p::AbstractPickle, ::Val{OpCodes.$op}, @nospecialize(arg))
    value = $popf(p.stack)
    obj = first(p.stack)
    if obj isa Defer
      if obj.head == $tag
        push!(obj.args, value)
      else
        wrap!(obj, $tag, value)
      end
    else
      $pushf(obj, value)
    end
  end
end

execute!(p::AbstractPickle, ::Val{OpCodes.LIST}, arg) = push!(p.stack, unmark!(p.stack))
execute!(p::AbstractPickle, ::Val{OpCodes.TUPLE}, arg) = push!(p.stack, Tuple(unmark!(p.stack)))

execute!(p::AbstractPickle, ::Val{OpCodes.TUPLE1}, arg) = push!(p.stack, tuple(pop!(p.stack)))
execute!(p::AbstractPickle, ::Val{OpCodes.TUPLE2}, arg) = push!(p.stack, reverse(pop2!(p.stack)))

execute!(p::AbstractPickle, ::Val{OpCodes.TUPLE3}, arg) = push!(p.stack, reverse(pop3!(p.stack)))

function execute!(p::AbstractPickle, ::Val{OpCodes.DICT}, arg)
  pairs = unmark!(p.stack)
  push!(p.stack, Dict(pairs[i]=>pairs[i+1] for i = 1:2:length(pairs)))
end

execute!(p::AbstractPickle, ::Val{OpCodes.POP}, arg) =
  isempty(p.stack) ? unmark!(p.stack) : pop!(p.stack)
execute!(p::AbstractPickle, ::Val{OpCodes.DUP}, arg) = push!(p.stack, first(p.stack))
execute!(p::AbstractPickle, ::Val{OpCodes.MARK}, arg) = mark!(p.stack)
execute!(p::AbstractPickle, ::Val{OpCodes.POP_MARK}, arg) = unmark!(p.stack)

for op in :(GET, BINGET, LONG_BINGET).args
  @eval execute!(p::AbstractPickle, ::Val{OpCodes.$op}, arg) =
    push!(p.stack, p.memo[Int(arg)])
end

for op in :(PUT, BINPUT, LONG_BINPUT).args
  @eval execute!(p::AbstractPickle, ::Val{OpCodes.$op}, arg) =
    setindex!(p.memo, first(p.stack), Int(arg))
end

execute!(p::AbstractPickle, ::Val{OpCodes.MEMOIZE}, arg) =
  setindex!(p.memo, first(p.stack), length(p.memo))

function execute!(p::AbstractPickle, ::Val{OpCodes.FROZENSET}, arg)
  fz = lookup(p.mt, "__main__", "frozenset")
  obj = isnothing(fz) ? Defer(gensym(:frozenset)) : fz()
  push!(p.stack, obj)
end

function execute!(p::AbstractPickle, ::Val{OpCodes.GLOBAL}, arg)
  func = lookup(p.mt, arg...)
  obj = isnothing(func) ? Defer(Symbol(join(arg, '.'))) : func
  push!(p.stack, obj)
end

function execute!(p::AbstractPickle, ::Val{OpCodes.STACK_GLOBAL}, arg)
  name = pop!(p.stack)
  scope = pop!(p.stack)
  func = lookup(p.mt, scope, name)
  obj = isnothing(func) ? Defer(Symbol(join((scope, name), '.'))) : func
  push!(p.stack, obj)
end

function execute!(p::AbstractPickle, ::Val{OpCodes.REDUCE}, arg)
  args = pop!(p.stack)
  f = first(p.stack)
  res = f isa Defer ? Defer(:reduce, f, args...) : f(args...)

  updatefirst!(p.stack, res)
end

execute!(p::AbstractPickle, ::Val{OpCodes.PROTO}, arg) = @assert protocol(p) >= arg "imcompatible protocol version: require version >= $arg"

# FRAMEING is ignored, but can be added if we found performance is bounded by io
for op in :(STOP, FRAME).args
  @eval execute!(p::AbstractPickle, ::Val{OpCodes.$op}, arg) = nothing
end

function execute!(p::AbstractPickle, ::Val{OpCodes.PERSID}, arg)
  f = lookup(p.mt, "persistent_load")
  obj = isnothing(f) ? Defer(:persistent_load, arg) : f(arg)
  push!(p.stack, obj)
end

function execute!(p::AbstractPickle, ::Val{OpCodes.BINPERSID}, arg)
  f = lookup(p.mt, "persistent_load")
  pid = pop!(p.stack)
  obj = isnothing(f) ? Defer(:persistent_load, pid) : f(pid)
  push!(p.stack, obj)
end

function execute!(p::AbstractPickle, ::Val{OpCodes.INST}, arg)
  f = lookup(p.mt, arg...)
  args = unmark!(p.stack)
  res = isnothing(f) ? Defer(:inst, f, args...) : f(args...)
  push!(p.stack, res)
end

function execute!(p::AbstractPickle, ::Val{OpCodes.OBJ}, arg)
  args = unmark!(p.stack)
  cls = popfirst!(args)
  res = cls isa Defer ? Defer(:obj, cls, args...) : cls(args...)
  push!(p.stack, res)
end

function execute!(p::AbstractPickle, ::Val{OpCodes.NEWOBJ}, arg)
  args, cls = pop2!(p.stack)
  res = cls isa Defer ? Defer(:newobj, cls, args...) : cls(args...)
  push!(p.stack, res)
end

function execute!(p::AbstractPickle, ::Val{OpCodes.NEWOBJ_EX}, arg)
  kwargs, args, cls = pop3!(p.stack)
  if cls isa Defer
    res = Defer(:newobj_ex, cls, kwargs, args...)
  else
    res = cls(args...; kwargs...)
  end
  push!(p.stack, res)
end

function objtypename(T)
  namestr = string(Base.typename(typeof(T)))
  if startswith(namestr, "typename(")
    namestr = namestr[10:end-1]
  end
  return namestr
end

function execute!(p::AbstractPickle, ::Val{OpCodes.BUILD}, arg)
  args = pop!(p.stack)
  obj = first(p.stack)
  if obj isa Defer
    wrap!(obj, :build, args)
  else
    build = lookup(p.mt, "__build__", objtypename(obj)) # hack for dispatch
    if isnothing(build)
      newobj = Defer(:build, obj, args)
    else
      newobj = build(obj, args)
    end
    hasref(p.memo, obj) &&
      maybeupdate!(p.memo, obj, newobj)
    maybeupdatefirst!(p.stack, newobj)
  end
end

# PROTOCOL 5 not supported
# execute!(p::AbstractPickle, ::Val{OpCodes.NEXT_BUFFER}, arg) = nothing
# execute!(p::AbstractPickle, ::Val{OpCodes.READONLY_BUFFER}, arg) = nothing

# Extension not supported, fire an issue if you need it.
# execute!(p::AbstractPickle, ::Val{OpCodes.EXT1}, arg) = read_uint1
# execute!(p::AbstractPickle, ::Val{OpCodes.EXT2}, arg) = read_uint2
# execute!(p::AbstractPickle, ::Val{OpCodes.EXT4}, arg) = read_int4


text:
###Test a Standard Metropolis-Hastings MCMC for a sine target model

###Innitialize variables

#MCMC iterations and burnin iterations
niterations = 100
nburnin = 10

#Generalized MH with multiple proposals per iteration
nproposals = 50

#Time points to simulate the sine-cosine model
timepoints = linspace(0.0,10.0,100)

###Model parameter values (coefficients of sine and cosine)
alower = 2.5
areal = 3.0
aupper = 3.5

###The variance of the normal noise on the data
fcov = [0.01]

###Create a sine-cosine with measurement data and uniform priors on the parameters
fparas = parameters([:a],[alower],[aupper],[areal])
fdata = data(:array,timepoints,sin(3.0*timepoints))
fnoise = noise(:gaussian,fcov)
fmodel = model(:target,fparas,fdata,fnoise,(t,p)->sin(p[1]*t);name="SineTestModel")
show(fmodel)

###Create different samplers
fnparas = numparas(fmodel)
fsampler1 = sampler(:mh,:normal,0.1,eye(1))
fsampler2 = sampler(:adaptive,0.1,fnparas)

###Create different tuners for the samplers
ftuner1 = tuner(:scale,10,0.5,:erf)
ftuner2 = tuner(:monitor,10)

###Create a Generalized Metropolis-Hastings runner
p = policy(:mh,nproposals)
r = runner(p,niterations,nproposals;numburnin=nburnin)
show(r)

###Run the MCMC with MH sampler
c1 = run!(r,fmodel,fsampler1,ftuner1)
show(c1)

###Run the MCMC with Adaptive MH sampler
c2 = run!(r,fmodel,fsampler2,ftuner2)
show(c2)


text:
using Test
using QuantumOptics
using LinearAlgebra

@testset "timeevolution_twolevel" begin

basis = SpinBasis(1//2)

# Random 2 level Hamiltonian
a1 = 0.5
a2 = 1.9
c = 1.3
d = -4.7

data = [a1 c-1im*d; c+1im*d a2]
H = DenseOperator(basis, data)

a = (a1 + a2)/2
b = (a1 - a2)/2
r = [c d b]

sigma_r = c*sigmax(basis) + d*sigmay(basis) + b*sigmaz(basis)

U(t) = exp(-1im*a*t)*(cos(norm(r)*t)*one(basis) - 1im*sin(norm(r)*t)*sigma_r/norm(r))

# Random initial state
psi0 = randstate(basis)
T = [0:0.5:1;]

f(t, psi::Ket) = @test 1e-5 > norm(psi - U(t)*psi0)
timeevolution.schroedinger(T, psi0, H; fout=f)
timeevolution.mcwf(T, psi0, H, []; fout=f)

f(t, rho::DenseOpType) = @test 1e-5 > tracedistance(rho, dm(U(t)*psi0))
timeevolution.master(T, psi0, H, []; fout=f)

end # testset


text:
using Printf
using Dates

# const WindThreshold=800
# const windSpaceMin = 500#meters
# const windIdStartPos = 0#meters
# const aboveThresholdPercent= 50

# const WindingChangeTimeMaxInterval = 2 #seconds

# const windAlarmSwitchONMinTIME = 1*60/2
# const windAlarmSwitchOFFMinTIME = 1*60/2

const baseDatetime ="20170211"*"000000"
include("../../common/windLogLst.jl")

# println(Base.@__MODULE__ )
include("../../common/common.jl")
# curpath = pwd()

# jlPath =joinpath(curpath ,"../Jl-Aux/glbCfg/WindIDCfg.jl")

# getCfg( Base.@__MODULE__,jlPath)
# println("\n",windAlarmSwitchOFFMinTIME,"\n")
mutable struct Analsys

    SegID1::UInt
    SegID2::UInt

    AvgValue::UInt

    L0Status::Bool
    L0time::UInt

    L1Status::Bool
    L1time::UInt

    L2Status::Bool
    L2time::UInt

    final_Status::Bool
    final_time::UInt # the new Alarm Start Time
    Last_finalTime::UInt ## Last Alarm start TIme
    Switched ::Bool
    function Analsys()
        new()
    end
end

function l0_handle(aobj::Analsys ,dt::String, reatime_eventOn::Bool ,avg::UInt )
    aobj.AvgValue= avg
    aobj.L0Status = reatime_eventOn
    curtime =  DT2UInt( dt )
    if aobj.L0time == UInt( 0 )
        # println( "-----",aobj.L0time ," DT:",dt )
        aobj.L0time = curtime
        aobj.L1time = curtime
        aobj.L2time = curtime
        aobj.final_time = curtime
        # exit(-1)
    else
        # println( "---Not Zero  --",aobj.L0time ," DT:",dt )
        aobj.L0time = curtime
        # exit(1)
    end
#    println( curtime ," ,,", aobj.final_time)
end

function l1_handle(aobj::Analsys )
    if aobj.L1Status  != aobj.L0Status
        aobj.L1time = aobj.L0time
        aobj.L1Status  = aobj.L0Status
    end
end

function l2_handle(aobj::Analsys )
    L1_persistTime =  aobj.L0time -aobj.L1time
    if L1_persistTime > WindingChangeTimeMaxInterval
        if  aobj.L2Status != aobj.L1Status
            aobj.L2Status = aobj.L1Status
            aobj.L2time  = aobj.L1time
        end
    end
end

function final_handle(aobj::Analsys )
    L2_persistTime =  aobj.L0time - aobj.L2time
    if   true == aobj.L2Status  && L2_persistTime > windAlarmSwitchONMinTIME
        if aobj.final_Status != aobj.L2Status
            aobj.final_Status = aobj.L2Status
            aobj.Last_finalTime = aobj.final_time
            aobj.final_time = aobj.L2time
            aobj.Switched =true
            return
        end
    end

    if  false == aobj.L2Status  && L2_persistTime > windAlarmSwitchOFFMinTIME
        if aobj.final_Status != aobj.L2Status
            aobj.final_Status = aobj.L2Status
            aobj.Last_finalTime = aobj.final_time
            aobj.final_time = aobj.L2time
            aobj.Switched = true
            return
        end
    end
    aobj.Switched = false
end

function init(  aObj::Analsys)
    # aObj.chnID =0
    aObj.SegID1=aObj.SegID2 =0x00;
    # dtValue=0;
    aObj.AvgValue =0;
    aObj.L0Status=false
    aObj.L0time=0

    aObj.L1Status=false
    aObj.L1time=0

    aObj.L2Status=false
    aObj.L2time=0

    aObj.final_Status=false
    aObj.final_time=0

    Switched =false
end

# global SpaceSegs = []#::Array{UInt,1}
# global WindAlarmSegs = []#::Array{Analsys,1  }

function init!(SpaceSegs,mertersPerSeg::Float64 ,beginID::UInt16,endID::UInt16 )
    println( "windSpaceMin=$windSpaceMin")
    step = floor(Int,windSpaceMin/mertersPerSeg)
    start = floor(Int,windIdStartPos/mertersPerSeg) +beginID
    maxSegID = min(  floor(Int,50*1000/mertersPerSeg)  , endID)

    for pos = start : step : maxSegID
        append!(SpaceSegs,[ pos,min(pos+step, maxSegID) ] )
    end
    #  println(start , " :",step," :",maxSegID ," ---:",beginID, " ",endID)
end

function init!(WindAlarmSegs,SpaceSegs )
    for idx = 1: 2 : length(SpaceSegs )
        aobj = Analsys()
        init(  aobj )
        aobj.SegID1 = SpaceSegs[idx]
        aobj.SegID2 = SpaceSegs[idx+1]
        # println( ":",SpaceSegs[idx]," --  ",SpaceSegs[idx+1]  )##############
        append!(WindAlarmSegs,[aobj])
    end
end

function isin( es1::UInt ,es2::UInt  ,a )
    len1 = length(a)
    for i = 1: len1
        aobj = a[i]
        as1 = aobj.SegID1
        as2 = aobj.SegID2
        inrng = ( es2 >= as1 && es2 <= as2 ) ||
            ( es1 >= as1 && es1 <= as2 ) ||
            (es1 < as1 && es2 > as2 )
        if inrng == true
            isWind = aobj.final_Status  ### 获取大风的状态，final_Status
            return isWind
        end
    end
    println( es1," ~ ",es2," not found winding on ")
    return false
end

function isinWindRng(es1::UInt ,es2::UInt  )
    if  es1 >= CalChn2Obj.startID || es2 >= CalChn2Obj.startID
        # println("Chnid =2 ",CalChn2Obj.startID)
        return isin( es1,es2, CalChn2Obj.WindAlarmSegs )
    else
        # println("Chnid =1")
        return isin( es1,es2, CalChn1Obj.WindAlarmSegs )
    end
end
function isinWindRng_1(es1::UInt ,es2::UInt  )
    inwind = isin( es1,es2, CalChn1Obj.WindAlarmSegs )
    if inwind == nothing  ## 未发现该区段
         inwind = isin( es1,es2, CalChn2Obj.WindAlarmSegs )
         if inwind == nothing
                 return false
         else
                return false
         end

    else
        return inwind
    end
end

################################################

function showAlarmWindingSwitch(winding::Analsys )
    # return
    if winding.Switched == true
        persisttime = winding.L0time -  winding.final_time
        if winding.final_Status == 1
            windStr = "1"
        else
            windStr = "0"
        end
        dtstr = UInt2DTStr( winding.L0time )
        info1 = @sprintf( "%llu,ID=%04u:%04u,EventFirsttime：%llu,L0time:%s, Winding: %s\r\n",
        winding.L0time,winding.SegID1,winding.SegID2,
           winding.final_time, dtstr, windStr)
        fileLog_wind.LogEvent(info1);

    end
end

const title1 = "curtime,SegRange,Winding,StartTime,EndTime,LastingTime,CurDateTime\r\n"

function showAlarmWinding(winding::Analsys )
    # return
    persisttime =winding.L0time -  winding.final_time
    @printf( "From ID:%04u,To ID:%04u,time:%s,L0time：%s,persistTIme:%llu,Winding:%u\n",
    winding.SegID1,winding.SegID2,
        UInt2DTStr(winding.final_time), UInt2DTStr( winding.L0time ), persisttime,winding.final_Status)
end

function showAlarmWinding(winding::Analsys )
    # return
    persisttime =winding.L0time -  winding.final_time
    @printf( "From ID:%04u,To ID:%04u,time:%s,L0time：%s,persistTIme:%llu,Winding:%u\n",
    winding.SegID1,winding.SegID2,
        UInt2DTStr(winding.final_time), UInt2DTStr( winding.L0time ), persisttime,winding.final_Status)
end

function showAlarmWindingL1(winding::Analsys )
    # return
    persisttime =winding.L0time -  winding.L1time
    @printf( "From ID:%04u,To ID:%04u,time:%s,persistTIme:%llu,Winding:%u\n",
    winding.SegID1,winding.SegID2,
        UInt2DTStr(winding.L1time), persisttime,winding.L1Status)
end

function showAlarmWindingL2(winding::Analsys )
    # return
    persisttime =winding.L0time -  winding.L2time
    @printf( "From ID:%04u,To ID:%04u,time:%s,persistTIme:%llu,Winding:%u\n",
    winding.SegID1,winding.SegID2,
        UInt2DTStr(winding.L2time), persisttime,winding.L2Status)
end
function updateWinding(aobj::Analsys,curDTStr::String,reatime_eventOn::Bool ,curAvgValue::UInt  )
    l0_handle( aobj ,curDTStr, reatime_eventOn ,curAvgValue )
    l1_handle( aobj   )
    l2_handle( aobj   )
    final_handle( aobj )
    showAlarmWindingSwitch(aobj )
    return
    # showAlarmWinding( aobj )
    # println( winding.isWind ,"---$(winding.status)-----------$(winding.lastOntime)---------------")
end

function FeatureAnal(a ,windingObject ,dt)
    cnt =0
    val =0
    segNum = length( a )
    for val0 in a
        val+=  val0
        if val0  > WindThreshold
            cnt +=1
            end
    end
    # avg=0
    # ratio=10
    # try
        avg = floor( UInt,val / segNum )
        ratio = floor( UInt32,cnt / (1*segNum) *100 )
    # return

    # print(  UInt(avg) ," ,",ratio ,"---")
    if (ratio > aboveThresholdPercent)
        updateWinding(windingObject , dt , true ,UInt(avg) )
        # return true
    else
        updateWinding(windingObject ,dt , false ,UInt(avg)  )
        # return false;
    end

end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

@testset "test_triangle_barycenter" begin
     triangle1 = Triangle()
     @test Raytracing.triangle_barycenter(triangle1) ≈ Point(√3/6, 0.0, 0.0)

     P1 = Point(1.0, 1.0, 1.0)
     P2 = Point(3.0, -1.0, 1.0)
     P3 = Point(2.0, 0.0, 4.0)
     triangle2 = Triangle(P1,P2,P3)

     @test Raytracing.triangle_barycenter(triangle2) ≈ Point(2.0, 0.0, 2.0)
end

@testset "test_Hit" begin
     triangle = Triangle()
     barycenter =  Raytracing.triangle_barycenter(triangle)
     ray1 = Ray(barycenter + 2.0*VEC_Z, -VEC_Z)

     intersection1 = ray_intersection(triangle, ray1)

     @test typeof(intersection1) == HitRecord
     @test HitRecord(
          Raytracing.triangle_barycenter(triangle),
          Normal(0.0, 0.0, 1.0),
          Vec2d(1/3, 1/3),
          2.0,
          ray1
          ) ≈ intersection1
end

@testset "test_noHit" begin
     triangle = Triangle()

     ray1 = Ray(Point(0, 0, 1), VEC_X)
     intersection1 = ray_intersection(triangle, ray1)
     @test isnothing(intersection1)

     ray2 = Ray(Point(1, 2, -10), VEC_X+VEC_Y)
     intersection2 = ray_intersection(triangle, ray2)
     @test isnothing(intersection2)

     ray3 = Ray(Point(10, 0, 1), -VEC_Z)
     intersection3 = ray_intersection(triangle, ray3)
     @test isnothing(intersection3)
end

@testset "test_general_triangle" begin
     P1 = Point(1.0, 0.5, 0.0)
     P2 = Point(1.0, -0.5, 0.0)
     P3 = Point(1.0, 0.0, √3/2)
     triangle = Triangle(SVector{3}(P1,P2, P3))

     ray1 = Ray(Point(0, 0,  √3/6), VEC_X)
     intersection1 = ray_intersection(triangle, ray1)
     @test typeof(intersection1) == HitRecord
     @test HitRecord(
          Raytracing.triangle_barycenter(triangle),
          Normal(-1.0, 0.0, 0.0),
          Vec2d(1/3, 1/3),
          1.0,
          ray1
          ) ≈ intersection1

     # Check if the triangle failed to move by trying to hit the basic shape
     @test isnothing( ray_intersection(triangle, Ray( Point(0.5, 0, 1), -VEC_Z ) ) )
end

@testset "test_Normals" begin
     P1 = Point(1.0, 1.0, 1.0)
     P2 = Point(3.0, -1.0, 1.0)
     P3 = Point(2.0, 0.0, √6+1)
     triangle = Triangle(SVector{3}(P1,P2,P3))

     ray1 = Ray(Point(0, 0, (2+√6)/3), VEC_X)
     intersection1 = ray_intersection(triangle, ray1)
     @test typeof(intersection1) == HitRecord
     @test intersection1.normal ≈ Normal(-√2/2, -√2/2, 0.0)

     ray2 = Ray(Point(2.0, 1.0, (2+√6)/3), -VEC_Y)
     intersection2 = ray_intersection(triangle, ray2)
     @test typeof(intersection2) == HitRecord
     @test intersection2.normal ≈ Normal(√2/2, √2/2, 0.0)
end

@testset "test_UV_Coordinates" begin
     triangle = Triangle()

     barycenter = Raytracing.triangle_barycenter(triangle)
     @test Raytracing.triangle_point_to_uv(triangle, barycenter) ≈ Vec2d(1/3, 1/3)
     ray1 = Ray(Point(barycenter.x, 0.0, 1.0), -VEC_Z)
     @test ray_intersection(triangle, ray1).surface_point ≈ Vec2d(1/3, 1/3)

     P2 = Point(√3/3, 0.0, 0.0)
     @test Raytracing.triangle_point_to_uv(triangle, P2) ≈ Vec2d(1/6, 1/6)
     ray2 = Ray(Point(P2.x, 0.0, 1.0), -VEC_Z)
     @test ray_intersection(triangle, ray2).surface_point ≈ Vec2d(1/6, 1/6)

     P3 = Point(√3/12, 1/4, 0.0)
     @test Raytracing.triangle_point_to_uv(triangle, P3) ≈ Vec2d(2/3, 1/6)
     ray2 = Ray(Point(P3.x, P3.y, 1.0), -VEC_Z)
     @test ray_intersection(triangle, ray2).surface_point ≈ Vec2d(2/3, 1/6)

     P4 = Point(0.0, 0.5, 0.0)
     @test Raytracing.triangle_point_to_uv(triangle, P4) ≈ Vec2d(1.0, 0.0)
end


text:
@everywhere push!(LOAD_PATH, "../")
@everywhere push!(LOAD_PATH, "./")

include("swimmerfunctions.jl")

####################################### Experiment Parameters / HyperParameters

model_file = "$(modeldir)/swimmer.xml"

# simulator parameters
T = 500
numT = 32
skip = 5
niter = 30

# solver parameters
gamma = 0.995
gae = 0.98

fullFIM = true

norm_step_size = 0.05

cg_iter= 12
cg_reg = 1e-6
cg_tol = 1e-10

#### deterministic randomization
srand(12345) # random seed, set globally

dtype = Float64

my_mjsys       = mjw.load_model(model_file, skip, "normal")

if METHOD == :NPG ######################################### NPG startup
   my_policy      = Policy.GLP{dtype}(my_mjsys.ns, my_mjsys.nu) # inputs: n, m
   #my_policy      = Policy.NN{dtype}(my_mjsys.ns, my_mjsys.nu, 32, 2) # inputs: n, m

   pgmodel        = PolicyGradModel(copy(my_policy.theta),
                                    my_mjsys,
                                    my_policy,
                                    T, numT, 10, niter)

   # NPG baseline aka value function approximation
   baseline    = Baseline.Quadratic{Float64}(1e-5, my_mjsys.ns, T, numT)
   #baseline    = Baseline.NN{Float64}(my_mjsys.ns, T, numT, 128)

   pg_specs   = NPGStrategy{dtype}(pgmodel,
                                   baseline,
                                   fullFIM,
                                   norm_step_size,
                                   gamma, gae,
                                   cg_iter, cg_reg, cg_tol)

elseif METHOD == :MPPI ######################################### MPPI startup
   mpcT = 100
   mppi = MPPIModel(eye(my_mjsys.nu)*0.1, # sigma
                    0.3,                  # lambda
                    my_mjsys,
                    mpcT, # horizon for mpc; set to T for non-MPC mode
                    T,
                    numT;
                    theta = randn(my_mjsys.nu, mpcT),
                    s0 = [my_mjsys.d.qpos; my_mjsys.d.qvel], # set initial state here
                   )

   mppi_specs = MPPIStrategy(mppi)
else
   error("WHAT YOU DOING")
end


text:
using UKrig
using Test

@testset "UKrig.jl" begin
   #FIXME ... add more examples
end


text:
export RainbowLearner

using Flux
using Zygote
using StatsBase
using Random

mutable struct RainbowLearner{
    Tq<:AbstractApproximator,
    Tt<:AbstractApproximator,
    Tf,
    Ts,
    Tss<:Union{Int,Nothing},
    R<:AbstractRNG,
} <: AbstractLearner
    approximator::Tq
    target_approximator::Tt
    loss_func::Tf
    Vₘₐₓ::Float32
    Vₘᵢₙ::Float32
    n_actions::Int
    n_atoms::Int
    support::Ts
    stack_size::Tss
    delta_z::Float32
    γ::Float32
    batch_size::Int
    update_horizon::Int
    min_replay_history::Int
    update_freq::Int
    target_update_freq::Int
    update_step::Int
    default_priority::Float64
    rng::R
end

function RainbowLearner(;
    approximator,
    target_approximator,
    loss_func,
    Vₘₐₓ,
    Vₘᵢₙ,
    n_actions,
    n_atoms = 51,
    support = collect(range(Float32(-Vₘₐₓ), Float32(Vₘₐₓ), length = n_atoms)),
    stack_size = 4,
    delta_z = Float32(support[2] - support[1]),
    γ = 0.99,
    batch_size = 32,
    update_horizon = 1,
    min_replay_history = 32,
    update_freq = 1,
    target_update_freq = 500,
    update_step = 0,
    default_priority = 100.0,
    seed = nothing,
)
    copyto!(approximator, target_approximator)  # force sync
    support = send_to_device(device(approximator), support)
    rng = MersenneTwister(seed)
    RainbowLearner(
        approximator,
        target_approximator,
        loss_func,
        Vₘₐₓ,
        Vₘᵢₙ,
        n_actions,
        n_atoms,
        support,
        stack_size,
        delta_z,
        γ,
        batch_size,
        update_horizon,
        min_replay_history,
        update_freq,
        target_update_freq,
        update_step,
        default_priority,
        rng,
    )
end

function (learner::RainbowLearner)(obs)
    state = send_to_device(device(learner.approximator), get_state(obs))
    logits = batch_estimate(learner.approximator, state)
    q = learner.support .* softmax(reshape(logits, :, learner.n_actions))
    # probs = vec(sum(q, dims=1)) .+ legal_action
    vec(sum(q, dims = 1)) |> send_to_host
end

function RLBase.update!(learner::RainbowLearner, batch)
    learner.update_step += 1
    learner.update_step % learner.update_freq == 0 || return nothing

    Q, Qₜ, γ, loss_func, n_atoms, n_actions, support, delta_z, update_horizon, batch_size =
        learner.approximator,
        learner.target_approximator,
        learner.γ,
        learner.loss_func,
        learner.n_atoms,
        learner.n_actions,
        learner.support,
        learner.delta_z,
        learner.update_horizon,
        learner.batch_size

    states, rewards, terminals, next_states = map(
        x -> send_to_device(device(Q), x),
        (batch.states, batch.rewards, batch.terminals, batch.next_states),
    )
    actions = CartesianIndex.(batch.actions, 1:batch_size)
    target_support =
        reshape(rewards, 1, :) .+
        (reshape(support, :, 1) * reshape((γ^update_horizon) .* (1 .- terminals), 1, :))

    updated_priorities = send_to_device(device(Q), Vector{Float32}())

    next_logits = batch_estimate(Qₜ, next_states)
    next_probs = reshape(softmax(reshape(next_logits, n_atoms, :)), n_atoms, n_actions, :)
    next_q = reshape(sum(support .* next_probs, dims = 1), n_actions, :)
    # next_q_argmax = argmax(cpu(next_q .+ next_legal_actions), dims=1)
    next_prob_select = select_best_probs(next_probs, next_q)

    target_distribution = project_distribution(
        target_support,
        next_prob_select,
        support,
        delta_z,
        learner.Vₘᵢₙ,
        learner.Vₘₐₓ,
    )

    gs = gradient(Flux.params(Q)) do
        logits = reshape(batch_estimate(Q, states), n_atoms, n_actions, :)
        select_logits = logits[:, actions]
        batch_losses = loss_func(select_logits, target_distribution)

        updated_priorities = vec(clamp.(sqrt.(batch_losses .+ 1f-10), 1.f0, 1.f2))
        target_priorities = 1.0f0 ./ sqrt.(updated_priorities .+ 1f-10)
        normalized_target_priorities = target_priorities ./ maximum(target_priorities)

        mean(Zygote.dropgrad(normalized_target_priorities) .* batch_losses)
    end

    update!(Q, gs)

    if learner.update_step % learner.target_update_freq == 0
        copyto!(Qₜ, Q)
    end

    updated_priorities |> send_to_host
end

@inline function select_best_probs(probs, q)
    q_argmax = argmax(q, dims = 1)
    prob_select = @inbounds probs[:, q_argmax] # !!! without @inbounds it would be really slow
    reshape(prob_select, :, length(q_argmax))
end

function project_distribution(supports, weights, target_support, delta_z, vmin, vmax)
    batch_size, n_atoms = size(supports, 2), length(target_support)
    clampped_support = clamp.(supports, vmin, vmax)
    tiled_support = reshape(
        repeat(clampped_support; outer = (n_atoms, 1)),
        n_atoms,
        n_atoms,
        batch_size,
    )

    projection =
        clamp.(
            1 .- abs.(tiled_support .- reshape(target_support, 1, :)) ./ delta_z,
            0,
            1,
        ) .* reshape(weights, n_atoms, 1, batch_size)
    reshape(sum(projection, dims = 1), n_atoms, batch_size)
end

function RLBase.extract_experience(t::AbstractTrajectory, learner::RainbowLearner)
    s = learner.stack_size
    h = learner.update_horizon
    n = learner.batch_size
    γ = learner.γ
    if length(t) > learner.min_replay_history
        # 1. sample indices based on priority
        inds = Vector{Int}(undef, n)
        valid_ind_range = isnothing(s) ? (1:(length(t)-h)) : (s:(1:(length(t)-h)))
        for i in 1:n
            ind, p = sample(learner.rng, get_trace(t, :priority))
            while ind ∉ valid_ind_range
                ind, p = sample(learner.rng, get_trace(t, :priority))
            end
            inds[i] = ind
        end

        # 2. extract SARTS
        states = consecutive_view(get_trace(t, :state), inds; n_stack = s)
        actions = consecutive_view(get_trace(t, :action), inds)
        next_states = consecutive_view(get_trace(t, :state), inds .+ h; n_stack = s)
        consecutive_rewards = consecutive_view(get_trace(t, :reward), inds; n_horizon = h)
        consecutive_terminals =
            consecutive_view(get_trace(t, :terminal), inds; n_horizon = h)
        rewards, terminals = zeros(Float32, n), fill(false, n)

        # 3. make sure that we only consider experiences in current episode for rewards and terminals
        for i in 1:n
            m = findfirst(view(consecutive_terminals, :, i))
            if isnothing(m)
                terminals[i] = false
                rewards[i] = discount_rewards_reduced(view(consecutive_rewards, :, i), γ)
            else
                terminals[i] = true
                rewards[i] = discount_rewards_reduced(view(consecutive_rewards, 1:m, i), γ)
            end
        end
        inds,
        (
            states = states,
            actions = actions,
            rewards = rewards,
            terminals = terminals,
            next_states = next_states,
        )
    else
        nothing
    end
end

function RLBase.update!(p::QBasedPolicy{<:RainbowLearner}, t::AbstractTrajectory)
    indexed_experience = extract_experience(t, p)
    if !isnothing(indexed_experience)
        inds, experience = indexed_experience
        priorities = update!(p.learner, experience)
        if !isnothing(priorities)
            get_trace(t, :priority)[inds] .= priorities
        end
    end
end

function (
    agent::Agent{<:QBasedPolicy{<:RainbowLearner},<:CircularCompactPSARTSATrajectory}
)(
    ::PostActStage,
    obs,
)
    push!(
        agent.trajectory;
        reward = get_reward(obs),
        terminal = get_terminal(obs),
        priority = agent.policy.learner.default_priority,
    )
    nothing
end


text:
"""
    ScalarBatchAccumulator(net, state, batch_size)

A ScalarBatchAccumulator is used to evaluate the contribution to
local observables ⟨σ|Ô|σ'⟩ψ(σ')/ψ(σ) , but by computing the neural
network ψ(σ) in batches of size `batch_size`. This is essential
to extrract a speedup when using GPUs.

This is an internal implementation detail of NeuralQuantum, and should
not be relied upon.

Once constructed, a ScalarBatchAccumulator is supposed to be used as follows:
- if `isfull(sba) == true` you should not push new elements to it (an error
will be throw otherwise)
- data is pushed as `sba(⟨σ|Ô|σ'⟩, σ', ψ(σ))`.
The configuration should be passed as a vector (if ket) or as a tuple
of two vectors (if density matrix).
"""
mutable struct AccumulatorObsScalar{A,B,C,D} <: AbstractObservableAccumulator
    res::A

    ψ_σ::B       # Buffers alls the <σ|ψ> of the denominator
    mel_buf::C      # ⟨σ|Ô|σ'⟩ in the buffer

    logψ_acc::D
end

function AccumulatorObsScalar(net::NeuralNetwork, hilb, v, batch_sz)
    bnet      = cached(net, batch_sz)
    CT        = Complex{real(out_type(net))}
    v         = state(hilb, bnet)

    acc       = AccumulatorLogPsi(bnet, v, batch_sz)

    w         = trainable_first(net)
    RT        = real(eltype(w))

    res       = zero(eltype(out_similar(bnet)))
    mel_buf   = similar(w, CT, 1, batch_sz)

    return AccumulatorObsScalar(
        res,
        res,
        mel_buf, acc)
end

@forward AccumulatorObsScalar.logψ_acc Base.length, Base.count,
    isfull

accum(a::AccumulatorObsScalar) = a.logψ_acc

function init!(c::AccumulatorObsScalar, σ, ψ_σ)
    c.ψ_σ     = ψ_σ
    c.mel_buf .= 0
    c.res     = 0

    init!(accum(c), σ)
    return nothing
end

# This resets the value at the end of a process_accumulatr
function reset!(c::AccumulatorObsScalar, σ)
    c.mel_buf .= 0

    init!(accum(c), σ)
    return nothing
end

function (c::AccumulatorObsScalar)(mel::Number, cngs_l, cngs_r, v)
    isfull(c) && (process_accumulator!(c); reset!(c, v))

    # If the matrix element is zero, don't do anything
    mel == 0.0 && return acc

    n_cngs_l = isnothing(cngs_l) ? 0 : length(cngs_l)
    n_cngs_r = isnothing(cngs_r) ? 0 : length(cngs_r)

    # If there are no changes, just sum it
    if n_cngs_l == 0 && n_cngs_r == 0
        c.res += mel
    else
        accum(c)(cngs_l, cngs_r)
        @inbounds c.mel_buf[count(c)] = mel
    end

    return c
end

function (c::AccumulatorObsScalar)(mel::Number, cngs, v)
    isfull(c) && (process_accumulator!(c); reset!(c, v))

    # If the matrix element is zero, don't do anything
    mel == 0.0 && return c

    n_cngs = isnothing(cngs) ? 0 : length(cngs)

    # If there are no changes, just sum it
    if n_cngs == 0
        c.res += mel
    else
        accum(c)(cngs)
        @inbounds c.mel_buf[count(c)] = mel
    end

    return c
end

"""
    process_accumulator!(c)

Processes all states stored in the accumulator, by computing their
relative local contribution.

It is safe to call this even if the accumulator is not full. In this
case all data beyond the count should be disregarded as it was
not initialized.

The output will be returned. You should not assume ownership of
the output, as it is preallocated and will be used for further
computations of the accumulator.
"""
function process_accumulator!(c::AccumulatorObsScalar)
    count(c) == 0 && return nothing

    ψ_σp = process_accumulator!(accum(c))

    c.mel_buf .*= exp.(ψ_σp .- c.ψ_σ)
    c.res += sum(c.mel_buf)

    return nothing
end

function finalize!(c::AccumulatorObsScalar)
    process_accumulator!(c)
    return c.res
end


text:
# import Particles package
using Particles

# parameters
Δt = 0.1    # time step size
nt = 200    # number of time steps
np = 10000  # number of particles
nx = 16     # number of grid points

# random initial conditions
x₀ = randn(np)
v₀ = randn(np)

# shift x₀ to the interval [0,1]
xmax = ceil(maximum(abs.(x₀)))
x₀ .+= xmax
x₀ ./= 2*xmax

# concatenate x₀ and v₀
z₀ = [ collect(z) for z in zip(x₀, v₀) ]

# vector field
function lorentz_force!(ż, z, p::PoissonSolver)
    ż[1] = z[2]
    ż[2] = eval_field(p, z[1])
end

# create Poisson solver
p = PoissonSolver{eltype(x₀)}(nx)
solve!(p, x₀)

# create an Equation instance with callbacks
equ = Equation((ż, z) -> lorentz_force!(ż, z, p), z₀; f_preproc = z -> solve!(p, reduce(hcat, z)[1,:]))

# create Simulation instance
sim = Simulation(equ, ExplicitEuler, Δt, nt)

# run simulation
z = sim()

# load Plots package
using Plots

# compute plot ranges
vmax = ceil(maximum(abs.(v₀)))
xlim = (0, 1)
vlim = (-vmax, +vmax)

# plot initial condition
z₀plot = reduce(hcat, z₀)
scatter(mod.(z₀plot[1,:], 1), z₀plot[2,:],
        marker = 3,
        xlim = xlim,
        ylim = vlim,
        title = "Vlasov-Poisson",
        legend = false,
        size = (800, 600)
)

# save figure to file
savefig("vlasov_poisson_z₀.png")

# create animation
# anim = @animate for n in 0:nt
#     scatter(mod.(z[1,:,n], 1), z[2,:,n],
#         marker = 3,
#         xlim = xlim,
#         ylim = vlim,
#         title = "Vlasov-Poisson",
#         legend = false,
#         size = (800, 600)
#     )
# end

# save animation to file
# gif(anim, "vlasov_poisson_anim.gif", fps=10)

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
struct ResCorrs{T<:AbstractFloat} <: AbstractVector{Array{T,2}}
    L_mats::Array{Array{T,2},1}

    num_params::Integer
    num_layers::Integer
end

function ResCorrs(vec::Array{T,1}, num_params::Integer, num_layers::Integer) where T<:Real
    ResCorrs{T}(
        [reshape(vec[((div(l*(l+1),2)-1)*num_params^2+1):((div((l+1)*(l+2),2)-1)*
        num_params^2)],num_params, (l+1)*num_params) for l in 1:num_layers],
        num_params,
        num_layers)
end

Base.vec(rescorrs::ResCorrs{T}) where T <: Real =
    Array{T}(vcat([vec(l) for l in rescorrs.L_mats]...))
Base.print(rescorrs::ResCorrs) = print(vec(rescorrs))
Base.length(::Type{ResCorrs}, num_params, num_layers) =
    num_params*num_params*(div((num_layers+2)*(num_layers+1),2)-1)
Base.size(rescorrs::ResCorrs; vec_size=false) = ifelse(vec_size,
    length(ResCorrs, rescorrs.num_params, rescorrs.num_layers),
    rescorrs.num_layers)

Base.getindex(rescorrs::ResCorrs{T}, idxs::Integer...) where T <: Real =
    rescorrs.L_mats[idxs...]


text:
function solve(
    problem::LeastSquaresProblem,
    x,
    method::Union{Adam,AdaMax,LineSearch,ParticleSwarm,TrustRegion},
    options::LeastSquaresOptions,
)
    td = OptimizationProblem(problem.residuals, problem.bounds)
    res = solve(td, x, method, OptimizationOptions())
end

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
""" Julia program to Check if a number is a Palindrome Number or not.
A palindrome number is a number such that the reverse of the number is equal to the number itself
"""

function check_palindrome(num)
    temp = num
    rev = 0
    while (temp > 0)
        rem = temp % 10
        rev = (rev * 10) + rem
        temp = temp ÷ 10;
    end
    if (num == rev)
        return true
    else
        return false
    end
end

print("Enter the number: ")
num = readline()
num = parse(Int, num)
res = check_palindrome(abs(num))
if res
    println("The given number $num is a Palindrome Number.")
else
    println("The given number $num is not a Palindrome Number.")
end

"""
Time Complexity: O(log(num)), where 'num' is the given number
Space Complexity: O(1)

SAMPLE INPUT AND OUTPUT

SAMPLE 1
Enter the number: 121
The given number 121 is a Palindrome Number.

SAMPLE 2
Enter the number: 1234
The given number 1234 is not a Palindrome Number.

"""


text:
# we write a function that returns a random vector from [0,1]^d with exponential density
function randExp(d, lam)

	# we generate a uniformly random vector
	y = rand(d)

	# we transform (elementwise) to random numbers in [exp(-lambda),1]
	y = exp(-lam) .+ (1-exp(-lam))*y

	# we now get numbers with exponential density
	x = -log.(y)/lam
	return x
end

max_stars_repo_path: 1a/Monte Carlo integration/randExp.jl

max_issues_repo_path: 1a/Monte Carlo integration/randExp.jl

max_forks_repo_path: 1a/Monte Carlo integration/randExp.jl


text:
"""
    MaybeHotMatrix{T, U, V} <: AbstractMatrix{V}

A matrix-like structure for representing one-hot encoded variables.
Like `Flux.OneHotMatrix` but supports `missing` values.

Construct with the [`maybehotbatch`](@ref) function.

See also: [`MaybeHotVector`](@ref), [`maybehot`](@ref).
"""
struct MaybeHotMatrix{T, U, V} <: AbstractMatrix{V}
    I::Vector{T}
    l::U
    MaybeHotMatrix(I::Vector{T}, l::U) where {T <: Integer, U <: Integer} = new{T, U, Bool}(I, l)
    MaybeHotMatrix(I::Vector{Missing}, l::U) where U <: Integer = new{Missing, U, Missing}(I, l)
    function MaybeHotMatrix(I::Vector{T}, l::U) where {T <: Maybe{Integer}, U <: Integer}
        new{T, U, Maybe{Bool}}(I, l)
    end

    function MaybeHotMatrix{T, U, V}(I, l) where {T <: Maybe{Integer}, U <: Integer, V <: Maybe{Bool}}
        new{T, U, V}(convert(Vector{T}, I), convert(U, l))
    end
end

MaybeHotMatrix(i::Integer, l::Integer) = MaybeHotMatrix([i], l)
MaybeHotMatrix(x::MaybeHotVector) = MaybeHotMatrix([x.i], x.l)

Base.size(X::MaybeHotMatrix) = (X.l, length(X.I))
Base.length(X::MaybeHotMatrix) = X.l * length(X.I)

Base.getindex(X::MaybeHotMatrix, idcs...) = (@boundscheck checkbounds(X, idcs...); _getindex(X, idcs...))
_getindex(X::MaybeHotMatrix, i::Union{Integer, AbstractVector}, j::Integer) = X.I[j] .== i
_getindex(X::MaybeHotMatrix, i::Integer, ::Colon) = X.I .== i
_getindex(X::MaybeHotMatrix, idcs::CartesianIndex{2}) = _getindex(X, Tuple(idcs)...)
_getindex(X::MaybeHotMatrix, ::Colon, i::Integer) = MaybeHotVector(X.I[i], X.l)
_getindex(X::MaybeHotMatrix, ::Colon, i::AbstractArray) = MaybeHotMatrix(X.I[i], X.l)
_getindex(X::MaybeHotMatrix, ::Colon, ::Colon) = MaybeHotMatrix(copy(X.I), X.l)

function Base.convert(::Type{<:MaybeHotMatrix{T, U}}, X::MaybeHotMatrix) where {T, U}
    MaybeHotMatrix(convert(Vector{T}, X.I), convert(U, X.l))
end

function Base.promote_rule(::Type{MaybeHotMatrix{T, U, V}},
        ::Type{MaybeHotMatrix{A, B, C}}) where {T, A, U, B, V, C}
    MaybeHotMatrix{promote_type(T, A), promote_type(U, B), promote_type(V, C)}
end

function _check_l(Xs::AbstractVecOrTuple{MaybeHotMatrix})
    l = Xs[1].l
    if any(!isequal(l), (X.l for X in Xs))
        DimensionMismatch(
            "Number of rows of `MaybeHotMatrix`s to hcat must correspond"
        ) |> throw
    end
    l
end

Base.hcat(Xs::T...) where T <: MaybeHotMatrix = _typed_hcat(T, Xs)
Base.hcat(Xs::MaybeHotMatrix...) = _typed_hcat(_promote_types(Xs...), Xs)
function _typed_hcat(::Type{T}, Xs::Tuple{Vararg{MaybeHotMatrix}}) where T <: MaybeHotMatrix
    T(vcat([X.I for X in Xs]...), _check_l(Xs))
end

Base.reduce(::typeof(hcat), Xs::Vector{<:MaybeHotMatrix}) = _typed_hcat(mapreduce(typeof, promote_type, Xs), Xs)
function _typed_hcat(::Type{T}, Xs::AbstractVector{<:MaybeHotMatrix}) where T <: MaybeHotMatrix
    T(reduce(vcat, [X.I for X in Xs]), _check_l(Xs))
end

A::AbstractMatrix * B::MaybeHotMatrix = (_check_mul(A, B); _mul(A, B))
Zygote.@adjoint A::AbstractMatrix * B::MaybeHotMatrix = (_check_mul(A, B); Zygote.pullback(_mul, A, B))

_mul(A::AbstractMatrix, B::MaybeHotMatrix{Missing}) = fill(missing, size(A, 1), size(B, 2))
_mul(A::AbstractMatrix, B::MaybeHotMatrix{<:Integer}) = A[:, B.I]
function _mul(A::AbstractMatrix, B::MaybeHotMatrix)
    C = zeros(Maybe{eltype(A)}, size(A, 1), size(B, 2))
    @inbounds for (k,i) in enumerate(B.I)
        C[:, k] = A * MaybeHotVector(i, B.l)
    end
    C
end

Flux.onehotbatch(X::MaybeHotMatrix{<:Integer}) = Flux.onehotbatch(X.I, 1:X.l)

"""
    maybehotbatch(ls, labels)

Return a [`MaybeHotMatrix`](@ref) in which each column corresponds to one element of `ls`
containing `1` at its first occurence in `labels` with all other elements set to `0`.

# Examples
```jldoctest
julia> maybehotbatch([:c, :a], [:a, :b, :c])
3×2 MaybeHotMatrix with eltype Bool:
 ⋅  1
 ⋅  ⋅
 1  ⋅

julia> maybehotbatch([missing, 2], 1:3)
3×2 MaybeHotMatrix with eltype Union{Missing, Bool}:
 missing    ⋅
 missing   true
 missing    ⋅
```

See also: [`maybehot`](@ref), [`MaybeHotMatrix`](@ref), [`MaybeHotVector`](@ref).
"""
maybehotbatch(L, labels) = MaybeHotMatrix([maybehot(l, labels).i for l in L], length(labels))

Base.hash(X::MaybeHotMatrix, h::UInt) where {T, U, V, W} = hash((X.I, X.l), h)
(X1::MaybeHotMatrix == X2::MaybeHotMatrix) = X1.I == X2.I && X1.l == X2.l
Base.isequal(X1::MaybeHotMatrix, X2::MaybeHotMatrix) = isequal(X1.I, X2.I) && X1.l == X2.l


text:
function cubetest(mdl::ModelInit,method_cube::String="CPU")

    mdl = init_model(mdl)
    mdl.Utils = execute_cube(mdl, method_cube)

    v1 = [1,1,1,1,1]
    v2 = [1,1,1,1,2]
    v3 = Int.(floor.( length.( mdl.GridDef.Grids ) / 2 ))
    v4 = push!(length.( mdl.GridDef.Grids[1:( length(mdl.GridDef.Grids) - 1 ) ] ),1 )
    v5 = length.( mdl.GridDef.Grids )
    test_val = Vector{Vector{Int}}([v1, v2, v3, v4, v5])

    out = Vector{Float64}()
    v = [v1,v2,v3,v4,v5]
    for v = test_val

        push!(out, mdl.Utils.f(v, mdl))

    end

    return( mdl.Utils.f.(v, [mdl]) )
    # return(out)

end


text:
#
# Initial author: Martin Otter, DLR-SR (first version: Jan. 14, 2017)
# License: MIT (expat)

"""
   module Media

Test module to evaluate how complex media could be implemented in (pure) Julia and then used in Modia.
"""
module Media

# Definition of medium states
abstract type MediumState end

mutable struct State_pT <: MediumState
   p::Float64
   T::Float64
end

mutable struct State_ph <: MediumState
   p::Float64
   h::Float64
end

# Definition of generic Medium and its functions
abstract type AbstractMedium end
density(               m::AbstractMedium,state::MediumState) = error("... function not defined for $m")
specificEnthalpy(      m::AbstractMedium,state::MediumState) = error("... function not defined for $m")
specificInternalEnergy(m::AbstractMedium,state::MediumState) = error("... function not defined for $m")

setState_pT(           m::AbstractMedium,p,T)                = error("... function not defined for $m")
setState_ph(           m::AbstractMedium,p,h)                = error("... function not defined for $m")

end

"""
    module SimpleWater

Simple medium model of water
"""
module SimpleWater
   import Media

   # Constants of medium that have the same value for every medium instance
   const cp_const = 4184
   const cv_const = 4184
   const d_const  = 995.586
   const T0       = 273.15

   # Constants of medium specific to a medium instance
   struct Medium <: Media.AbstractMedium
      h_offset::Float64
      Medium(;h_offset=0.0) = new(h_offset)
   end

   # Functions of medium
   Media.density(               m::Medium, state::Media.State_pT) = d_const
   Media.specificEnthalpy(      m::Medium, state::Media.State_pT) = cp_const*(state.T - T0) + m.h_offset
   Media.specificInternalEnergy(m::Medium, state::Media.State_pT) = cv_const*(state.T - T0)

   Media.setState_pT(m::Medium, p, T)::Media.State_pT = Media.State_pT(p,T)
   Media.setState_ph(m::Medium, p, h)::Media.State_pT = Media.State_pT(p,T0+(h-m.h_offset)/cp_const)
end

"""
    module Test_Media

Test media models
"""
module Test_Media

import Media
import SimpleWater

#=
@model FluidPort begin
  medium
  p=Float()
  m_flow=Float(flow=true)
end
=#

# Set a medium at one connector
medium = SimpleWater.Medium(h_offset=10.0)

# Propagate medium in connector

# Inside a model
state = Media.setState_ph(medium, 1e5, 300.0)   # p,h from connector
d     = Media.density(medium,state)
h     = Media.specificEnthalpy(medium,state)

println("d = ", d)
println("h = ", h)

end


text:
# module MarkovChain

# import ..JuliaRL

import Random
import Base.copy

module MarkovChainParams

const RIGHT = 1
const LEFT = 2

const ACTIONS = [RIGHT, LEFT]

end

mutable struct MarkovChain <: AbstractEnvironment
    agent_pos::Int64
    size::Int64
    MarkovChain(size) = new(size/2, size)
    MarkovChain(agent_pos, size) = new(agent_pos, size)
end

JuliaRL.is_terminal(env::MarkovChain) = (env.agent_pos == env.size) || (env.agent_pos == 1)
is_terminal(env::MarkovChain, state) = (state == env.size) || (state == 1)
JuliaRL.get_reward(env::MarkovChain) = if (env.agent_pos == env.size) 1.0 else 0.0 end
get_reward(env::MarkovChain, state) = if (state == env.size) 1.0 else 0.0 end
JuliaRL.get_state(env::MarkovChain) = env.agent_pos
Base.size(env::MarkovChain) = env.size
JuliaRL.get_actions(env::MarkovChain) = MarkovChainParams.ACTIONS
num_actions(env::MarkovChain) = 2
get_states(env::MarkovChain)=collect(1:env.size)

function JuliaRL.reset!(env::MarkovChain; rng=Random.GLOBAL_RNG, kwargs...)
    env.agent_pos = rand(rng, 2:(env.size - 1))
end

function JuliaRL.environment_step!(env::MarkovChain, action; rng=Random.GLOBAL_RNG, kwargs...)
    mcp = MarkovChainParams
    @assert !JuliaRL.is_terminal(env)
    if action == mcp.RIGHT
        env.agent_pos -= 1
    elseif action == mcp.LEFT
        env.agent_pos += 1
    else
        @error "Action invalid $(action)"
    end
end

function _step(env::MarkovChain, action, state; rng=Random.GLOBAL_RNG, kwargs...)
    mcp = MarkovChainParams
    @assert !JuliaRL.is_terminal(env)
    agent_pos = state
    if action == mcp.RIGHT
        agent_pos -= 1
    elseif action == mcp.LEFT
        agent_pos += 1
    else
        @error "Action invalid $(action)"
    end
    return agent_pos, get_reward(env, agent_pos), is_terminal(env, agent_pos)
end

function DynamicProgramming(env::MarkovChain, pi, gamma, thresh=1e-20)
    v = zeros(size(env))
    # v[size(env)] = 1
    delta = 100000
    while thresh < delta
        delta = 0
        for i = 2:9
            _v = v[i]
            r_left = 0
            r_right = if (i == size(env)-1) 1 else 0 end
            v[i] = pi[1]*(r_left + gamma*v[i-1]) + pi[2]*(r_right + gamma*v[i+1])
            delta = max(delta, v[i] - _v)
        end
    end
    return v
end

function Base.show(io::IO, env::MarkovChain)
    println("MARKOV CHAIN: $(get_state(env)), $(get_reward(env)), $(is_terminal(env))")
    model = fill("0", env.size)
    model[1] = "1"
    println(prod(model.*" "))
    model = fill("-", env.size)
    model[env.agent_pos] = "a"
    println(prod(model.*" "))
end