text:
using LightXML
using JuMP, MultilinearOpt, BARON, Gurobi, CPLEX

mutable struct Feed
    id::String
    fc::Float64
    uc::Float64
    lb::Float64
    ub::Float64
    qualities::Dict
end

mutable struct Pool
    id::String
    fc::Float64
    uc::Float64
    cap::Float64
    qualities::Dict
end

mutable struct Product
    id::String
    ur::Float64
    lb::Float64
    ub::Float64
    qualities::Dict
end

mutable struct Connection
    source::String
    destination::String
    fc::Float64
    uc::Float64
end

const Quality = String

c( x::Feed) = x.uc
AL(x::Feed) = x.lb
AU(x::Feed) = x.ub

SS(x::Pool) = x.cap

d(x::Product) = x.ur
DU(x::Product) = x.ub
DL(x::Product) = x.lb

CC(x::Feed, y::Quality) = x.qualities[y]
PL(x::Product, y::Quality) = x.qualities[y][1]
PU(x::Product, y::Quality) = x.qualities[y][2]

# I = feeds, J = products, L = pools, K = qualities
function pooling_problem(I, J, L, K, conns)
    pairs = [(x.source,x.destination) for x in conns]
    Tx(i,l) = (i.id,l.id) in pairs
    Ty(l,j) = (l.id,j.id) in pairs
    Tz(i,j) = (i.id,j.id) in pairs

    m = Model(solver=CplexSolver())
    @variable(m, 0 ≤ q[i in I, l in L] ≤ (Tx(i,l) ? 1 : Inf))
    @variable(m, 0 ≤ y[l in L, j in J] ≤ min(SS(l), DU(j), sum(AU(i) for i in I if Tx(i,l))))
    @variable(m, 0 ≤ z[i in I, j in J] ≤ min(AU(i), DU(j)))

    # Standard Pooling Problem

    @objective(m, Min, sum(c(i)*q[i,l]*y[l,j] for i in I, l in L, j in J if Tx(i,l) && Ty(l,j)) - sum(d(j)*y[l,j] for l in L, j in J if Ty(l,j)) - sum((d(j)-c(i))*z[i,j] for i in I, j in J if Tz(i,j)))

    # Availability bounds
    for i in I
        if AL(i) > 0
            @constraint(m, sum(q[i,l]*y[l,j] for j in J, l in L if Tx(i,l) && Ty(l,j)) + sum(z[i,j] for j in J if Tz(i,j)) ≥ AL(i))
        end
        if AU(i) < Inf
            @constraint(m, sum(q[i,l]*y[l,j] for j in J, l in L if Tx(i,l) && Ty(l,j)) + sum(z[i,j] for j in J if Tz(i,j)) ≤ AU(i))
        end
    end

    # Pool Capacity ========================
    for l in L
        if SS(l) < Inf
            @constraint(m, sum(y[l,j] for j in J if Ty(l,j)) ≤ SS(l))
        end
    end

    # Product Demand =======================
    for j in J
        if DL(j) > 0
            @constraint(m, sum(y[l,j] for l in L if Ty(l,j)) + sum(z[i,j] for i in I if Tz(i,j)) ≥ DL(j))
        end
        if DU(j) < Inf
            @constraint(m, sum(y[l,j] for l in L if Ty(l,j)) + sum(z[i,j] for i in I if Tz(i,j)) ≤ DU(j))
        end
    end

    # Product Quality ======================
    for j in J, k in K
        if PL(j,k) > 0
            @constraint(m, PL(j,k)*(sum(y[l,j] for l in L if Ty(l,j)) + sum(z[i,j] for i in I if Tz(i,j))) ≤ sum(CC(i,k)*z[i,j] for i in I if Tz(i,j)) + sum(CC(i,k)*q[i,l]*y[l,j] for i in I, l in L if Tx(i,l) && Ty(l,j)))
        end
        if PU(j,k) < Inf
            @constraint(m, PU(j,k)*(sum(y[l,j] for l in L if Ty(l,j)) + sum(z[i,j] for i in I if Tz(i,j))) ≥ sum(CC(i,k)*z[i,j] for i in I if Tz(i,j)) + sum(CC(i,k)*q[i,l]*y[l,j] for i in I, l in L if Tx(i,l) && Ty(l,j)))
        end
    end

    # Simplex Constraint ===================
    for l in L
        @constraint(m, sum(q[i,l] for i in I if Tx(i,l)) == 1)
    end

    # PQ Cut ===============================
    for l in L, j in J
        @constraint(m, sum(q[i,l]*y[l,j] for i in I if Tx(i,l)) == y[l,j])
    end

    m
end

filebase = "/Users/huchette/Dropbox\ (MIT)/research/multilinear_relaxations/sample_problems/sample_problems/standard_lit_problems"

fp = open("/Users/huchette/Dropbox\ (MIT)/research/multilinear_relaxations/pooling_results.csv", "w")

for (roots, dirs, files) in walkdir(filebase), filename in files
    fi = joinpath(filebase,filename)
    xdoc = parse_file(fi)
    xroot = root(xdoc)

    feeds = Any[]
    pools = Any[]
    prods = Any[]
    conns = Any[]
    quals = Any[]
    for ce in child_elements(xroot)
        if name(ce) == "Feed"
            qualities = Dict(attribute(x,"quality_id") => parse(Float64, content(find_element(x,"Level"))) for x in get_elements_by_tagname(ce, "QualityLevel"))
            feed = Feed(attribute(ce, "id"),
                        parse(Float64, content(find_element(ce,"FixedCost"))),
                        parse(Float64, content(find_element(ce,"UnitCost"))),
                        parse(Float64, content(find_element(ce,"LowerBound"))),
                        parse(Float64, content(find_element(ce,"UpperBound"))),
                        qualities)
            push!(feeds, feed)
        elseif name(ce) == "Pool"
            qualities = Dict(attribute(x,"quality_id") => parse(Float64, content(find_element(x,"PercentReduction"))) for x in get_elements_by_tagname(ce, "QualityLevel"))
            pool = Pool(attribute(ce, "id"),
                        parse(Float64, content(find_element(ce,"FixedCost"))),
                        parse(Float64, content(find_element(ce,"UnitCost"))),
                        parse(Float64, content(find_element(ce,"Capacity"))),
                        qualities)
            push!(pools, pool)
        elseif name(ce) == "Product"
            qualities = Dict(attribute(x,"quality_id") => (parse(Float64, content(find_element(x,"LowerBound"))),parse(Float64, content(find_element(x,"UpperBound")))) for x in get_elements_by_tagname(ce, "QualityLevel"))
            product = Product(attribute(ce, "id"),
                           parse(Float64, content(find_element(ce,"UnitRevenue"))),
                           parse(Float64, content(find_element(ce,"LowerBound"))),
                           parse(Float64, content(find_element(ce,"UpperBound"))),
                           qualities)
            push!(prods, product)
        elseif name(ce) == "Connection"
            connection = Connection(attribute(ce, "source"),
                                    attribute(ce, "destination"),
                                    parse(Float64, content(find_element(ce,"FixedCost"))),
                                    parse(Float64, content(find_element(ce,"UnitCost"))))
            push!(conns, connection)
        elseif name(ce) == "Quality"
            push!(quals, attribute(ce, "id"))
        end
    end

    # MisenerLinear
    m1 = pooling_problem(feeds, prods, pools, quals, conns)
    relaxbilinear!(m1, method=:MisenerLinear, disc_level=33)
    t1 = @elapsed optimize!(m1)
    obj1 = objective_value(m1)
    println(fp, "$filename, MisenerLinear, $t1, $obj1")

    # MisenerLog1
    m2 = pooling_problem(feeds, prods, pools, quals, conns)
    relaxbilinear!(m2, method=:MisenerLog1, disc_level=33)
    t2 = @elapsed optimize!(m2)
    obj2 = objective_value(m2)
    println(fp, "$filename, MisenerLog1, $t2, $obj2")

    # MisenerLog2
    m3 = pooling_problem(feeds, prods, pools, quals, conns)
    relaxbilinear!(m3, method=:MisenerLog2, disc_level=33)
    t3 = @elapsed optimize!(m3)
    obj3 = objective_value(m3)
    println(fp, "$filename, MisenerLog2, $t3, $obj3")

    # Logarithmic 1D
    m4 = pooling_problem(feeds, prods, pools, quals, conns)
    relaxbilinear!(m4, method=:Logarithmic1D, disc_level=33)
    t4 = @elapsed optimize!(m4)
    obj4 = objective_value(m4)
    println(fp, "$filename, Logarithmic1D, $t4, $obj4")

    # Logarithmic 2D
    m5 = pooling_problem(feeds, prods, pools, quals, conns)
    relaxbilinear!(m5, method=:Logarithmic2D, disc_level=33)
    t5 = @elapsed optimize!(m5)
    obj5 = objective_value(m5)
    println(fp, "$filename, Logarithmic2D, $t5, $obj5")

    # ZigZag 1D
    m6 = pooling_problem(feeds, prods, pools, quals, conns)
    relaxbilinear!(m6, method=:ZigZag1D, disc_level=33)
    t6 = @elapsed optimize!(m6)
    obj6 = objective_value(m6)
    println(fp, "$filename, ZigZag1D, $t6, $obj6")

    # ZigZag 2D
    m7 = pooling_problem(feeds, prods, pools, quals, conns)
    relaxbilinear!(m7, method=:ZigZag2D, disc_level=33)
    t7 = @elapsed optimize!(m7)
    obj7 = objective_value(m7)
    println(fp, "$filename, ZigZag2D, $t7, $obj7")

    flush(fp)
end

close(fp)


text:
module MiniLogging

export get_logger, basic_config

include("Hierarchy.jl")
using .Hierarchy

const LogLevel = Int

const DEFINED_LEVELS = Dict(
    0 => :NOTSET,
    10 => :DEBUG,
    20 => :INFO,
    30 => :WARN,
    40 => :ERROR,
    50 => :CRITICAL,
)

for (value, symbol) in DEFINED_LEVELS
    @eval $symbol = $value
end

mutable struct Handler
    output::IO
    date_format::String
end

mutable struct Logger
    name::String
    level::LogLevel
    handlers::Vector{Handler}
end

function Base.show(io::IO, logger::MiniLogging.Logger)
    level_symbol = get(DEFINED_LEVELS, logger.level, nothing)
    if level_symbol == nothing
        level_str = string(logger.level)
    else
        level_str = string(level_symbol, ":", logger.level)
    end

    if (is_root(logger.name))
        print(io, "RootLogger($level_str)")
    else
        print(io,
            """Logger("$(logger.name)", $level_str)"""
        )
    end
end

Logger(name::String, level::LogLevel) = Logger(name, level, Handler[])

const TREE = Tree()
const ROOT = Logger("", WARN)
const LOGGERS = Dict{String, Logger}("" => ROOT)

is_root(name::String) = name == "" || name == "Main"

get_logger() = ROOT

"""
- `@assert get_logger("") == get_logger() == get_logger("Main")`
"""
function get_logger(name::String)::Logger
    if is_root(name)
        return get_logger()
    end

    if haskey(LOGGERS, name)
        return LOGGERS[name]
    end

    push!(TREE, name)
    logger = Logger(name, NOTSET)
    LOGGERS[name] = logger
    logger
end

get_logger(name) = get_logger(string(name))

is_not_set(logger::Logger) = logger.level == NOTSET

function get_effective_level(logger::Logger)::LogLevel
    logger_name = logger.name
    while !is_root(logger_name)
        if !is_not_set(logger)
            return logger.level
        end
        logger_name = parent_node(TREE, logger.name)
        logger = LOGGERS[logger_name]
    end
    # This is `ROOT`.
    return logger.level
end

is_enabled_for(logger::Logger, level::LogLevel) = level >= get_effective_level(logger)

has_handlers(logger::Logger) = !isempty(logger.handlers)

function get_effective_handlers(logger::Logger)::Vector{Handler}
    logger_name = logger.name
    while !is_root(logger_name)
        if has_handlers(logger)
            return logger.handlers
        end
        logger_name = parent_node(TREE, logger.name)
        logger = LOGGERS[logger_name]
    end
    # This is `ROOT`.
    return logger.handlers
end

function basic_config(level::LogLevel; date_format::String="%Y-%m-%d %H:%M:%S")
    ROOT.level = level
    handler = Handler(stderr, date_format)
    push!(ROOT.handlers, handler)
end

function basic_config(level::LogLevel, file_name::String; date_format::String="%Y-%m-%d %H:%M:%S", file_mode::String="a")
    ROOT.level = level
    f = open(file_name, file_mode)
    handler = Handler(f, date_format)
    push!(ROOT.handlers, handler)
end

write_log(output::T, color::Symbol, msg::AbstractString) where T <: IO = (print(output, msg); flush(output))
write_log(output::Base.TTY, color::Symbol, msg::AbstractString) = Base.printstyled(output, msg, color=color)

function _log(
        logger::Logger, level::LogLevel, color::Symbol,
        msg...
    )
    logger_name = is_root(logger.name) ? "Main" : logger.name
    level_symbol = get(DEFINED_LEVELS, level, nothing)
    if level_symbol == nothing
        level_str = string(level)
    else
        level_str = string(level_symbol)
    end

    for handler in get_effective_handlers(logger)
        t = Libc.strftime(handler.date_format, time())
        s = string(t, ":" , level_str, ":", logger_name, ":" , msg..., "\n")
        write_log(handler.output, color, s)
    end
end

function define_macro(macro_name::Symbol, level::LogLevel, color::Symbol)
    @eval macro $macro_name(logger, msg...)
        level = $level
        # This generates e.g. `:red`.
        color = $(Expr(:quote, color))
        msg = map(esc, msg)
        quote
            logger = $(esc(logger))
            if is_enabled_for(logger, $level)
                _log(logger, $level, $(Expr(:quote, color)), $(msg...))
            end
        end
    end
    callable_macro_name = Symbol("@", macro_name)
    @eval export $callable_macro_name
end

for (macro_name, level, color) in [
        (:debug, DEBUG, :cyan),
        (:info, INFO, :blue),
        (:warn, WARN, :magenta),
        (:error, ERROR, :red),
        (:critical, CRITICAL, :red)
    ]
    define_macro(macro_name, level, color)
end

function define_new_level(
        macro_name::Symbol,
        level::LogLevel,
        color::Symbol=":yellow"
    )
    DEFINED_LEVELS[level] = macro_name
    define_macro(macro_name, level, color)
end

end

ext: jl


text:
using transD_GP
##
ntimesperdecade = 2
T = 10 .^(-3:1/ntimesperdecade:6)
f = 1 ./T
## model 1 conductor above
ρ = [10, 100]
z = [0, 5000]
## make a synthetic
F = transD_GP.MT1DInversion.create_synthetic(ρ, z, f, rseed=11)
transD_GP.get_misfit(F.d_log10_ρ, F.d_phase_deg, F.σ_log10_ρ, F.σ_phase_deg, f, ρ, z)


text:
# Adaptation of BayesLogit C code by Nicholas Polson, James Scott, Jesse Windle, 2012-2019

const __PI = 3.14159265358979323846264338327950288419716939937510582097494459230781640628
const _TERMS = 100

include("invert_Y.jl")
include("polyagamma_devroye.jl")
include("polyagamma_sp.jl")
include("polyagamma_alt.jl")
include("polyagamma_normal.jl")

struct PolyaGamma{T1<:Real,T2<:Real} <: Distributions.ContinuousUnivariateDistribution
    h::T1
    z::T2
    PolyaGamma{T1,T2}(h, z) where {T1,T2} = new{T1,T2}(h, z)
end

function PolyaGamma(h::T1, z::T2; check_args = true) where {T1<:Real,T2<:Real}
    check_args && Distributions.@check_args(PolyaGamma, h > zero(h))
    return PolyaGamma{T1,T2}(h, z)
end

#### Outer constructors
PolyaGamma(h::T) where {T<:Real} = PolyaGamma(h, one(T))
PolyaGamma() = PolyaGamma(1, 1.0, check_args = false)

#@distr_support PolyaGamma -Inf Inf #TODO

#### Parameters

scale(d::PolyaGamma) = d.h
tilting(d::PolyaGamma) = d.z

params(d::PolyaGamma) = (d.h, d.z)
partype(::PolyaGamma{T1,T2}) where {T1,T2} = (T1, T2)

#### Statistics

mean(d::PolyaGamma) = d.h / (2 * d.z) * tanh(d.z / 2) # ord.h / (2 * d.z) * ((_exp_c(d.z) - 1) / (1 + _exp_c(d.z)))

var(d::PolyaGamma) = d.h / (4 * d.z^3) * (sinh(d.z) - d.z) * sech(d.z / 2)^2

## DEPRECATED IN BAYESLOGIT
## Calculate coefficient n in density of PG(1.0, 0.0), i.e. J* from Devroye.
##------------------------------------------------------------------------------
function _a_coef(n, x)
    if (x > _TRUNC)
        return __PI * (n + 0.5) * _exp_c(-(n + 0.5)^2 * __PI^2 * x / 2)
    else
        return (2 / __PI / x)^1.5 * __PI * (n + 0.5) * _exp_c(-2 * (n + 0.5)^2 / x)
    end
end

## DEPRECATED IN BAYESLOGIT
function _jacobi_logpdf(z::Float64, x::Float64; ntrunc::Int = _TERMS)
    v = mapreduce(n -> (iseven(n) ? 1 : -1) * _a_coef(n, x), +, 0:ntrunc)
    return _log_cosh_c(z) - x * z^2 / 2 + _log_c(v)
end

"""
    The log coefficients of the infinite sum for the density of PG(b, 0).
    See [^Polson], section 2.3.

    [^Polson]     See [^Polson], section 2.3.

"""
function _pg_logcoef(x, b, n)
    loggamma(n + b) - loggamma(n + 1) - loggamma(b) + _log_c(2n + b) -
    _log_c(2π * x^3) / 2 - (2n + b)^2 / 8x
end

"""
   log density of the PG(b, 0) distribution.
   See [^Polson], section 2.3.

   [^Polson]: Polson, Nicholas G., Scott, James G., & Windle, J., (2013). Bayesian Inference for Logistic Models Using P olya–Gamma Latent Variables, Journal of the American Statistical Association, 108:504, 1339–1349, DOI: 10.1080/01621459.2013.
"""
function _pg0_logpdf(x::Float64, b::Float64; ntrunc::Int = _TERMS)
    v = zero(x)
    for n = 0:ntrunc
        v += (iseven(n) ? 1 : -1) * _exp_c(_pg_logcoef(x, b, n))
    end
    return (b - 1) * _log_c(2) + _log_c(abs(v))
end

"""
    log density of the PG(b, c) distribution.
    See Polson et al. 2013, section 2.2 and equation (5).

    [^Polson]: Polson, Nicholas G., Scott, James G., & Windle, J., (2013). Bayesian Inference for Logistic Models Using P olya–Gamma Latent Variables, Journal of the American Statistical Association, 108:504, 1339–1349, DOI: 10.1080/01621459.2013.
"""
function _pg_logpdf(b::Float64, c::Float64, x::Float64; ntrunc::Int = _TERMS)
    b * _log_cosh_c(c / 2) - x * c^2 / 2 + _pg0_logpdf(x, b; ntrunc = ntrunc)
end

function logpdf(d::PolyaGamma, x::Float64; ntrunc::Int = _TERMS)
    if d.h == 1
        return _jacobi_logpdf(d.z / 2, 4 * x; ntrunc = ntrunc) + _log_c(4)
    else
        return _pg_logpdf(d.h, d.z, x; ntrunc = ntrunc)
    end
end

function Distributions.pdf(d::PolyaGamma, x::Real; ntrunc::Int = _TERMS)
    _exp_c(logpdf(d, x; ntrunc = ntrunc))
end

function Distributions.rand(rng::Distributions.AbstractRNG, d::PolyaGamma)
    if d.h == 1.0
        return Distributions.rand(rng, PolyaGammaDevRoye1Sampler(float(d.h), d.z))
    elseif d.h == 2.0
        return Distributions.rand(rng, PolyaGammaDevRoyeSampler(float(d.h), d.z))
    elseif d.h <= 13.0
        return Distributions.rand(rng, PolyaGammaAltSampler(float(d.h), d.z, 200))
    elseif d.h <= 170.0
        return Distributions.rand(rng, PolyaGammaSPSampler(float(d.h), d.z, 100))
    else
        m = pg_m1(d.h, d.z)
        v = pg_m2(d.h, d.z) - m * m
        return Distributions.rand(rng, Distributions.Normal(m, sqrt(v)))
    end
end


text:
module DEACpopulation
using Random;
const functionNames = ["randPop",
                       "flatPop",
                       "gaussPop",
                       "gaussPop2",
                       "gaussOne",
                       "gaussTwo",
                       "gaussThree"]
#Generate population functions
function randPop(rng::MersenneTwister,genomeSize::Int64,populationSize::Int64,omega::Array{Float64,1})
    rand(rng,genomeSize,populationSize)
end

function flatPop(rng::MersenneTwister,genomeSize::Int64,populationSize::Int64,omega::Array{Float64,1})
    rand(rng) .* ones(genomeSize,populationSize)
end

function gaussPop(rng::MersenneTwister,genomeSize::Int64,populationSize::Int64,omega::Array{Float64,1})
    p = (rand(rng,populationSize));
    mu = (rand(rng,populationSize)*maximum(omega));
    alpha = (rand(rng,populationSize)*maximum(omega)/4);
    (p'.*exp.(-((omega .- mu').^2)./(alpha'.^2)/2)./alpha'/sqrt(2*pi))
end

function gaussPop2(rng::MersenneTwister,genomeSize::Int64,populationSize::Int64,omega::Array{Float64,1})
    p = ones(populationSize);
    mu = Array{Float64}(undef,populationSize);
    for i = 0:populationSize-1
        mu[i+1]=omega[(i%genomeSize)+1];
    end
    alpha = (rand(rng::MersenneTwister,populationSize));
    (p'.*exp.(-((omega .- mu').^2)./(alpha'.^2)/2)./alpha'/sqrt(2*pi));
end

function gaussOne(rng::MersenneTwister,genomeSize::Int64,populationSize::Int64,omega::Array{Float64,1})
    gaussPop(rng,genomeSize,populationSize,omega)
end

function gaussTwo(rng::MersenneTwister,genomeSize::Int64,populationSize::Int64,omega::Array{Float64,1})
    gaussPop(rng,genomeSize,populationSize,omega) + gaussPop(rng,genomeSize,populationSize,omega)
end

function gaussThree(rng::MersenneTwister,genomeSize::Int64,populationSize::Int64,omega::Array{Float64,1})
    gaussPop(rng,genomeSize,populationSize,omega) + gaussPop(rng,genomeSize,populationSize,omega) + gaussPop(rng,genomeSize,populationSize,omega)
end

end


text:
#
## Constant params
const MOI_OPTIMAL_STATUSES = [
    MOI.OPTIMAL,
    MOI.ALMOST_OPTIMAL,
    MOI.LOCALLY_SOLVED,
    MOI.ALMOST_LOCALLY_SOLVED,
]

@enum(TargetDevice,
    CPU,
    CUDADevice,
    Mixed,
)

#
#
# Algorithmic parameters
#
"""
    AlgParams

Specifies ProxAL's algorithmic parameters.

| Parameter | Description | Default value |
| :--- | :--- | :--- |
| `decompCtgs::Bool` | if true: decompose across contingencies (along with time) | false
| `jacobi::Bool` |     if true: do Jacobi updates, else do Gauss-Siedel updates | true
| `iterlim::Int` |     maximum number of ProxAL iterations | 100
| `nlpiterlim::Int` |  maximum number of NLP subproblem iterations | 100
| `tol::Float64` |     tolerance used for ProxAL termination | 1.0e-4
| `zero::Float64` |    tolerance below which is regarded as zero | 1.0e-8
| `θ_t::Float64` |     see [Formulation](@ref) | 1.0
| `θ_c::Float64` |     see [Formulation](@ref) | 1.0
| `ρ_t::Float64` |     AL penalty weight for ramp constraints | 1.0
| `ρ_c::Float64` |     AL penalty weight for ctgs constraints | 1.0
| `updateρ_t::Bool` |  if true: dynamically update `ρ_t` | true
| `updateρ_c::Bool` |  if true: dynamically update `ρ_c` | true
| `τ::Float64`       | Proximal weight parameter | 3.0
| `updateτ::Bool` |    if true: dynamically update `τ` | true
| `verbose::Int` |     level of output: 0 (none), 1 (stdout) | 0
| `mode::Symbol` |     computation mode `∈ [:nondecomposed, :coldstart, :lyapunov_bound]` | `:coldstart`
| `optimizer::Any` |   NLP solver | `nothing`
| `gpu_optimizer::Any` | GPU-compatible NLP solver | `nothing`
| `nr_tol::Float64`    | Tolerance of the Newton-Raphson algorithm (used only in `ExaPFBackend()` model) | 1e-10
| `init_opt::Bool` |   if true: initialize block OPFs with base OPF solution | false
| `device::TargetDevice` | Target device to deport the resolution of the optimization problem | CPU
"""
Base.@kwdef mutable struct AlgParams
    decompCtgs::Bool = false # decompose contingencies (along with time)
    jacobi::Bool     = true  # if true: do jacobi, else do gauss-siedel
    iterlim::Int     = 100   # maximum number of ADMM iterations
    nlpiterlim::Int  = 100   # maximum number of NLP subproblem iterations
    tol::Float64     = 1e-3  # tolerance used for ADMM termination
    zero::Float64    = 1e-8  # tolerance below which is regarded as zero
    θ_t::Float64     = 1.0   # weight_quadratic_penalty_time
    θ_c::Float64     = 1.0   # weight_quadratic_penalty_ctgs
    ρ_t::Float64     = 1.0   # AL parameter for ramp constraints
    ρ_c::Float64     = 1.0   # AL parameter for ctgs constraints
    updateρ_t::Bool  = true  # Dynamically update ρ for ramp constraints
    updateρ_c::Bool  = true  # Dynamically update ρ for ctgs constraints
    τ::Float64       = 3.0   # Proximal coefficient
    updateτ::Bool    = true  # Dynamically update τ
    verbose::Int     = 0     # level of output: 0 (none), 1 (stdout), 2 (+plots), 3 (+outfiles)
    mode::Symbol            = :coldstart     # computation mode [:nondecomposed, :coldstart, :lyapunov_bound]
    optimizer::Any          = nothing        # NLP solver for fullmodel and subproblems
    gpu_optimizer::Any      = nothing        # GPU-compatible NLP solver for fullmodel and subproblems
    nr_tol::Float64         = 1e-10          # Tolerance of the Newton-Raphson algorithm (for ExaBlockBackend backend)
    init_opf::Bool          = false
    device::TargetDevice    = CPU
    verbose_inner::Int      = 0
    tron_rho_pq::Float64    = 4e2
    tron_rho_pa::Float64    = 4e4
    tron_scale::Float64     = 1e-4
    tron_inner_iterlim::Int = 800
    tron_outer_iterlim::Int = 20
    tron_outer_eps::Float64 = 1e-4
end

"""
    ModelInfo

Specifies the ACOPF model structure.

| Parameter | Description | Default value |
| :--- | :--- | :--- |
| `num_time_periods::Int` | number of time periods | 1
| `num_ctgs::Int` | number of line contingencies | 0
| `load_scale::Float64` | load multiplication factor | 1.0
| `ramp_scale::Float64` | multiply this with ``p_{g}^{max}`` to get generator ramping ``r_g`` | 1.0
| `corr_scale::Float64` | multiply this with ``r_g`` to get generator ramping for corrective control | 0.1
| `obj_scale::Float64` | objective multiplication factor | 1.0e-3
| `allow_obj_gencost::Bool` | model generator cost | true
| `allow_constr_infeas::Bool` | allow constraint infeasibility | false
| `allow_line_limits::Bool` | allow line flow limits | true
| `weight_constr_infeas::Float64` | quadratic penalty weight for constraint infeasibilities | 1.0
| `weight_freq_ctrl::Float64` | quadratic penalty weight for frequency violations | 1.0
| `weight_ctgs::Float64` | linear weight of contingency objective function | 1.0
| `case_name::String` | name of case file | ""
| `savefile::String` | name of save file | ""
| `time_link_constr_type::Symbol` | `∈ [:penalty, :equality, :inequality]` see [Formulation](@ref) | `:penalty`
| `ctgs_link_constr_type::Symbol` | `∈ [:frequency_penalty, :frequency_equality, :preventive_penalty, :preventive_equality, :corrective_penalty, :corrective_equality, :corrective_inequality]`, see [Formulation](@ref) | `:frequency_penalty`
"""
Base.@kwdef mutable struct ModelInfo
    num_time_periods::Int = 1
    num_ctgs::Int = 0
    load_scale::Float64 = 1.0
    ramp_scale::Float64 = 1.0
    corr_scale::Float64 = 0.1
    obj_scale::Float64 = 1e-3
    allow_obj_gencost::Bool = true
    allow_constr_infeas::Bool = false
    allow_line_limits::Bool = true
    weight_constr_infeas::Float64 = 1.0
    weight_freq_ctrl::Float64 = 1.0
    weight_ctgs::Float64 = 1.0
    case_name::String = ""
    savefile::String = ""
    time_link_constr_type::Symbol = :penalty # ∈ [:penalty, :equality, :inequality]
    ctgs_link_constr_type::Symbol = :frequency_penalty # ∈ [:frequency_penalty, :frequency_equality, :preventive_penalty, :preventive_equality, :corrective_penalty, :corrective_equality, :corrective_inequality]
end


text:
using Test
using Bobby

function test_fen(fen, depth, result)
    b = Bobby.fen_to_bitboard(fen)
    pt = Bobby.perft(b, depth, b.player_color)
    if pt.nodes == result
        return true
    else
        println(pt.nodes)
        return false
    end
end

@testset "parser.jl" begin
    include("test_parser.jl")
end

@testset "converters.jl" begin
    include("test_converters.jl")
end

@testset "nights.jl" begin
    include("test_nights.jl")
end

@testset "king.jl" begin
    include("test_king.jl")
end

@testset "rooks.jl" begin
    include("test_rooks.jl")
end

@testset "bishops.jl" begin
    include("test_bishops.jl")
end

@testset "queen.jl" begin
    include("test_queen.jl")
end

@testset "pawns.jl" begin
    include("test_pawns.jl")
end

@testset "check.jl" begin
    include("test_check.jl")
end

# @testset "perft.jl" begin
#     include("test_perft.jl")
# end


text:
global const MODE_H = Int8(1)
global const MODE_A = Int8(-1)
global const MODE_HA = Int8(0)

global const GRB_ENV = Gurobi.Env()
global const GRB_ENV_2 = Gurobi.Env()
global const GRB_ENV_3 = Gurobi.Env()

struct IdxSet
    binrepr::Int64
end

struct RR2
    nteams::Int
    nslots::Int
    isphased::Bool
end

abstract type Constraint end

struct CAcon <: Constraint
    teams1::IdxSet
    teams2::IdxSet
    slots::IdxSet
    ub::Int
    mode::Int
    pen::Int
end

struct BRcon <: Constraint
    teams::IdxSet
    slots::IdxSet
    ub::Int
    mode::Int
    pen::Int
end

struct GAcon <: Constraint
    slots::IdxSet
    lb::Int
    ub::Int
    meetings::Vector{Tuple{Int,Int}}
    pen::Int
end

struct FAcon <: Constraint
    teams::IdxSet
    slots::IdxSet
    ub::Int
    pen::Int
end

struct SEcon <: Constraint
    teams::IdxSet
    lb::Int
    pen::Int
end

struct Instance
    rr2::RR2
    constrCA::Vector{CAcon}
    constrBR::Vector{BRcon}
    constrGA::Vector{GAcon}
    constrFA::Vector{FAcon}
    constrSE::Vector{SEcon}
end

struct TabuSearchOptimizer
    maxiter::Int
    tabulength::Int
    infeaslength::Int
    rr2::RR2
    constraints::Vector{Constraint}
    conentries::BitArray{3}
end

const Move = Tuple{Symbol, Int64, Int64, Vararg{Int64, N} where N}

mutable struct SearchState
    it::Int # iteration count
    hpen::Int
    spen::Int # current value
    sched::Matrix{Int}
    optval::Int # best known value
    optsched::Matrix{Int} # best known solution
    lastimprov::Int # last iteration in which global improvement of obj function was fou
    offset::Vector{Int}
    moves::Vector{Move}
    tabu::Dict{Move,Int}
    blocked::Dict{Move,Int}
    pen::Vector{Float64} # penalty for each constraint
end

function IdxSet(A::Vector{T}) where {T<:Integer}
    binrepr = 0
    for i in A
        if binrepr & (1 << (i - 1)) !== 0
            throw(ArgumentError("Invalid Array: has duplicates"))
        end
        if i < 1 || i > 64
            throw(ArgumentError("Invalid Array: has values not in range 1:64"))
        end
        binrepr += 2^(i - 1)
    end
    return IdxSet(binrepr)
end

Base.eltype(::IdxSet) = Int
@inline Base.in(n::Int64, s::IdxSet) = s.binrepr & (1 << (n - 1)) !== 0
@inline Base.length(s::IdxSet) = count_ones(s.binrepr)

function Base.iterate(s::IdxSet, state = 1)
    trunc = s.binrepr >>> (state - 1)
    trunc === 0 && return nothing
    c = trailing_zeros(trunc)
    return state + c, state + c + 1
end

show(io::IO, x::IdxSet) = join(string.(i for i in x), ' ')

function Instance(rr2::RR2, cons::Vector{Constraint})
    D = Dict(CAcon => 1, BRcon => 2, GAcon => 3, FAcon => 4, SEcon => 5)
    groupedcons = (
        Vector{CAcon}(),
        Vector{BRcon}(),
        Vector{GAcon}(),
        Vector{FAcon}(),
        Vector{SEcon}(),
    )
    for c in cons
        push!(groupedcons[D[typeof(c)]], c)
    end
    return Instance(rr2, groupedcons...)
end

struct PGBACut
    ub::Int
    assignments::Vector{Tuple{Int,Int}}
end


text:
C1s = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1])]
C2s = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("C_2^1", [-1.0 -1.2246467991473532e-16 0.0; 1.2246467991473532e-16 -1.0 0.0; 0.0 0.0 1.0])]
C3s = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("C_3^1", [-0.4999999999999998 -0.8660254037844387 0.0; 0.8660254037844387 -0.4999999999999998 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_3^2", [-0.5000000000000003 0.8660254037844384 0.0; -0.8660254037844384 -0.5000000000000003 0.0; 0.0 0.0 1.0])]
C4s = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("C_4^1", [6.123233995736766e-17 -1.0 0.0; 1.0 6.123233995736766e-17 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2^1", [-1.0 -1.2246467991473532e-16 0.0; 1.2246467991473532e-16 -1.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_4^3", [-1.8369701987210297e-16 1.0 0.0; -1.0 -1.8369701987210297e-16 0.0; 0.0 0.0 1.0])]
C5s = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("C_5^1", [0.30901699437494745 -0.9510565162951535 0.0; 0.9510565162951535 0.30901699437494745 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_5^2", [-0.8090169943749473 -0.5877852522924731 0.0; 0.5877852522924731 -0.8090169943749473 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_5^3", [-0.8090169943749475 0.587785252292473 0.0; -0.587785252292473 -0.8090169943749475 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_5^4", [0.30901699437494734 0.9510565162951535 0.0; -0.9510565162951535 0.30901699437494734 0.0; 0.0 0.0 1.0])]
C6s = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("C_6^1", [0.5000000000000001 -0.8660254037844386 0.0; 0.8660254037844386 0.5000000000000001 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_3^1", [-0.4999999999999998 -0.8660254037844388 0.0; 0.8660254037844388 -0.4999999999999998 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2^1", [-1.0 -3.885780586188048e-16 0.0; 3.885780586188048e-16 -1.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_3^2", [-0.5000000000000006 0.8660254037844385 0.0; -0.8660254037844385 -0.5000000000000006 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_6^5", [0.49999999999999944 0.8660254037844392 0.0; -0.8660254037844392 0.49999999999999944 0.0; 0.0 0.0 1.0])]

C1cn = ["A"]
C1ct = [1.0]
C2cn = ["A", "B"]
C2ct = [1.0 1.0; 1.0 -1.0]
C3cn = ["A", "E"]
C3ct = [1.0 1.0 1.0; 2.0 -0.9999999999999996 -1.0000000000000009]
C4cn = ["A", "B", "E"]
C4ct = [1.0 1.0 1.0 1.0; 1.0 -1.0 1.0 -1.0; 2.0 1.2246467991473532e-16 -2.0 -3.6739403974420594e-16]
C5cn = ["A", "E1", "E2"]
C5ct =  [1.0   1.0        1.0        1.0        1.0;
         2.0   0.6180339887498949   -1.6180339887498947   -1.618033988749895    0.6180339887498945;
         2.0  -1.6180339887498947    0.6180339887498945    0.6180339887498955  -1.6180339887498953]
C6cn = ["A", "B", "E1", "E2"]
C6ct = [1.0   1.0   1.0   1.0   1.0   1.0;
        1.0  -1.0   1.0  -1.0   1.0  -1.0;
        2.0   1.0  -1.0  -2.0  -1.0   1.0;
        2.0  -1.0  -1.0   2.0  -1.0  -1.0]


text:
import JSON2

function validate_serialization(sys::System)
    #path, io = mktemp()
    # For some reason files aren't getting deleted when written to /tmp. Using current dir.
    path = "test_system_serialization.json"
    io = open(path, "w")
    @info "Serializing to $path"
    sys_ext = get_ext(sys)
    sys_ext["data"] = 5
    ext_test_bus_name = ""
    try
        IS.prepare_for_serialization!(sys.data, path)
        bus = collect(PSY.get_components(PSY.Bus, sys))[1]
        ext_test_bus_name = PSY.get_name(bus)
        ext = PSY.get_ext(bus)
        ext["test_field"] = 1
        to_json(io, sys)
    catch
        close(io)
        rm(path)
        rethrow()
    end
    close(io)

    ts_file = nothing
    try
        ts_file = open(path) do file
            JSON2.read(file).data.time_series_storage_file
        end
        sys2 = System(path)
        sys_ext2 = get_ext(sys2)
        sys_ext2["data"] != 5 && return false
        bus = PSY.get_component(PSY.Bus, sys2, ext_test_bus_name)
        ext = PSY.get_ext(bus)
        ext["test_field"] != 1 && return false
        return IS.compare_values(sys, sys2)
    finally
        @debug "delete temp file" path
        rm(path)
        rm(ts_file)
    end
end

@testset "Test JSON serialization of CDM data" begin
    sys = create_rts_system()
    @test validate_serialization(sys)
    text = JSON2.write(sys)
    @test length(text) > 0
end

@testset "Test JSON serialization of matpower data" begin
    sys = System(PowerSystems.PowerModelsData(joinpath(MATPOWER_DIR, "case5_re.m")))

    # Add a Probabilistic forecast to get coverage serializing it.
    bus = Bus(nothing)
    bus.name = "Bus1234"
    add_component!(sys, bus)
    tg = RenewableFix(nothing)
    tg.bus = bus
    add_component!(sys, tg)
    tProbabilisticForecast =
        PSY.Probabilistic("scalingfactor", Hour(1), DateTime("01-01-01"), [0.5, 0.5], 24)
    add_forecast!(sys, tg, tProbabilisticForecast)

    @test validate_serialization(sys)
end

@testset "Test JSON serialization of ACTIVSg2000 data" begin
    path = joinpath(DATA_DIR, "ACTIVSg2000", "ACTIVSg2000.m")
    sys = System(PowerSystems.PowerModelsData(path))
    validate_serialization(sys)
end


text:
@testset "Bounding boxes" begin
  # sphere
  @test boundingbox(Sphere(P2(0,0),T(1))) == Box(P2(-1,-1),P2(1,1))
  @test boundingbox(Sphere(P2(1,1),T(1))) == Box(P2(0,0),P2(2,2))

  # point set
  @test boundingbox(PointSet(T[0 1 2; 0 2 1])) == Box(P2(0,0), P2(2,2))
  @test boundingbox(PointSet(T[1 2; 2 1])) == Box(P2(1,1), P2(2,2))

  # cartesian grid
  @test boundingbox(CartesianGrid{T}(10,10)) == Box(P2(0,0), P2(10,10))
  @test boundingbox(CartesianGrid{T}(100,200)) == Box(P2(0,0), P2(100,200))
  @test boundingbox(CartesianGrid((10,10), T.((1,1)), T.((1,1)))) == Box(P2(1,1), P2(11,11))
end


text:
# Make a macro that will create matrix R and vector r for restricting the VAR
# The user will write beta[1,2,1] - beta[2,1,3] = 10 (beta[from, to, lag]) and through
# regular expressions this will get parsed into a vector of the matrix
# macro restrict(ex) gamaa restricts the coefficients

function β(from, to, lag)
	col = fill(0.0, (vars^2)*lags + vars*ntyp)
	println(1 + vars*ntyp + (lag-1)*(vars^2) + (from-1)*vars + (to-1))
	col[1 + vars*ntyp + (lag-1)*(vars^2) + (from-1)*vars + (to-1)] = 1.0
	return col
end

function betaa(from, to, lag)
	return β(from, to, lag)
end

function γ(eq)
	col = fill(0.0, (vars^2)*lags + vars*ntyp)
	col[eq] = 1.0
	return col
end

function gamaa(eq)
	return γ(eq)
end

macro restrictions(ex)
	ind = [2:2:length(ex.args)]
	R = apply(Expr, prepend!([:($(ex.args[i].args[1])) for i=ind], [:call, :hcat]));
	r = [:($(ex.args[i].args[2].args[2])) for i=ind];
	return :(tuple($R, $r))
end

# vars, lags, and ntyp have to be defined
# vars = e.vars
# ntyp = e.ntyp
# lags = e.lags

# @restrictions begin
#     gamaa(3)=1
# end

# @restrictions begin
#     gamaa(4) + betaa(1,2,3) + betaa(3,2,1) = 10
# 	betaa(3,4,2) - betaa(3,1,3) = 5
# end

# Heh, hard
# function makeR(R::Matrix{Float64})
# 	base = eye(size(R)[1])

# 	if numberInColumn==1
# 		base

# 	end
# end

# function testRestrictions(R::Matrix{Float64})
# 	@assume rank(R)==apply(max, size(mat)) "The constraints are not independent. One can be expressed as linear combination of other."
# 	# Alright, this is making me nuts :-D
# 	@assume all(mapslices((x) -> sum(map((y) -> y==0.0, x)), getR(mat), [2]) .>= apply(min, size(R))) "Your constants yield the model overidentified. For example, beta(1,2,3) + beta(2,3,1) = 0 and beta(1,2,3) + beta(2,2,1) = 0 imply only equality of beta(2,2,1) - beta(2,3,1) = 0 and does not depend on beta(1,2,3)."
# end


text:
module NNetTest

using OnlineAI, FactCheck

srand(1)

function xor_data()
    inputs = [0 0; 0 1; 1 0; 1 1]
    targets = float(sum(inputs,2) .== 1)

    # all sets are the same
    inputs = inputs .- mean(inputs,1)
    DataPoints(inputs, targets)
end

function testxor{LAYER}(::Type{LAYER};
        hiddenLayerNodes = [4],
        hiddenMapping = SoftsignMapping(),
        finalMapping = SigmoidMapping(),
        params = NetParams(),
        solverParams = SolverParams(maxiter=100000),
        inputTransformer = IdentityTransformer(),
        doPretrain = false)

    # all xor inputs and results
    inputs = [0 0; 0 1; 1 0; 1 1]
    targets = float(sum(inputs,2) .== 1)

    # all sets are the same
    inputs = inputs .- mean(inputs,1)
    data = DataPoints(inputs, targets)
    sampler = SimpleSampler(data)

    # hiddenLayerNodes = [2]
    net = buildRegressionNet(
        LAYER,
        ncols(inputs),
        ncols(targets),
        hiddenLayerNodes;
        hiddenMapping = hiddenMapping,
        finalMapping = finalMapping,
        params = params,
        solverParams = solverParams,
        inputTransformer = inputTransformer
    )
    show(net)

    if doPretrain
        pretrain(net, sampler, sampler)
    end

    stats = solve!(net, sampler, sampler)

    output = vec(predict(net, inputs))
    for (o, d) in zip(output, data)
        println("Result: input=$(d.x) target=$(d.y) output=$o")
    end

    net, output, stats
end

for LAYER in (Layer, NormalizedLayer)
    facts("NeuralNet{$LAYER}") do

        atol = 0.05
        solverParams = SolverParams(
            maxiter=50000,
            erroriter=10000,
            minerror=1e-3,
            plotiter=-1,
            plotfields=Symbol[:x, :xhat, :β, :α, :δy, :y, :w, :b, :δΣ, :Σ, :a]
        )

        updater = AdaMaxUpdater()
        net, output, stats = testxor(
            LAYER,
            params=NetParams(updater=updater, mloss=L2DistLoss()),
            solverParams=solverParams,
            doPretrain=false
        )
        @fact output --> roughly([0., 1., 1., 0.], atol=atol)

        net, output, stats = testxor(
            LAYER,
            params=NetParams(updater=AdagradUpdater(), mloss=L2DistLoss()),
            solverParams=solverParams
        )
        @fact output --> roughly([0., 1., 1., 0.], atol=atol)

        net, output, stats = testxor(
            LAYER,
            params=NetParams(updater=AdaMaxUpdater(), mloss=CrossentropyLoss()),
            finalMapping=SigmoidMapping(),
            solverParams=solverParams
        )
        @fact output --> roughly([0., 1., 1., 0.], atol=atol)

        solverParams.maxiter = 300000
        net, output, stats = testxor(
            LAYER,
            params=NetParams(updater=AdaMaxUpdater(), mloss=CrossentropyLoss(),
            dropout=Dropout(1.0,0.9)),
            finalMapping=SigmoidMapping(),
            solverParams=solverParams,
            hiddenLayerNodes = [6,6,6]
        )
        @fact output --> roughly([0., 1., 1., 0.], atol=atol)

    end # facts
end

function test_pretrain(; solve=true, pretr=true, netparams=NetParams(), kwargs...)
    data = xor_data()
    sampler = StratifiedSampler(data)

    # nin = 1; f = x->x[1:1]
    nin = 2; f = nop
    net = buildRegressionNet(
        Layer, nin,1,[2]; params=netparams,
        solverParams=SolverParams(maxiter=10000),
        inputTransformer=f
    )
    if pretr
        pretrain(net, sampler, sampler; kwargs...)
    end

    if solve
        solve!(net, sampler, sampler)
    end

    net, Float64[predict(net,d.x)[1] for d in data]
end

end # module
nn = NNetTest


text:
using Solc
using JSON

installSolc("latest")


text:
## Run this from the hitting_times_model directory
model_str = (@__DIR__)*"/absorbing_model/"
pth = model_str*"stationary_distribution/"
include(model_str*"model_def.jl")
mkpath(pth*"data")
mkpath(pth*"figs")

include((@__DIR__)*"default_params.jl")
# playing with meta programming, sorry its hard to read...
for (c_o,o) in enumerate(orders)
    for (c_m,m) in enumerate((:dg,:order1,:qbdrap))
        @eval $(Symbol(m,"_d0",o)),$(Symbol(m,"_B",o)),$(Symbol(m,"_dq",o)) = $(Symbol(m,"_model"))($o)
        @eval $(Symbol(m,"_B",o))[:,1] .= 1.0
        @eval bvec = zeros(1,size($(Symbol(m,"_B",o)),1))
        bvec[1] = 1.0
        @eval $(Symbol(m,"_π",o)) = SFMDistribution(bvec/Matrix($(Symbol(m,"_B",o))),$(Symbol(m,"_dq",o)))
        @eval $(Symbol(m,"_pdf",o)) = pdf($(Symbol(m,"_π",o)))
        @eval $(Symbol(m,"_cdf",o)) = cdf($(Symbol(m,"_π",o)))
    end
end

p1 = plot(layout = (3,length(orders)))
q1 = plot(layout = (3,length(orders)))

std_plot!(p,fun,args...; kwargs...) = plot!(
    p,
    fun,
    -eps(), model.b+eps();
    # ylim=(-0.1,2.1),
    xticks=false,
    yticks=false,
    grid=false,
    label=false,
    kwargs...,
)

_pdf, _cdf = sfm.stationary_distribution_x(model)
true_pdf(x,i) = _pdf(x)[i]
true_cdf(x,i) = _cdf(x)[i]

L1_cdf_errors = DataFrame(DG = Float64[], Order_1 = Float64[], QBDRAP = Float64[])
ks_errors = DataFrame(DG = Float64[], Order_1 = Float64[], QBDRAP = Float64[])
L1_pdf_errors = DataFrame(DG = Float64[], Order_1 = Float64[], QBDRAP = Float64[])
L2_pdf_errors = DataFrame(DG = Float64[], Order_1 = Float64[], QBDRAP = Float64[])

L1_cdf_error_row = zeros(1,3)
ks_error_row = zeros(1,3)
L1_pdf_error_row = zeros(1,3)
L2_pdf_error_row = zeros(1,3)

# what a truely terrible piece of code!
for (c_o,o) in enumerate(orders)
    for (c_m,m) in enumerate((:dg,:order1,:qbdrap))
        for i in 1:2
            # @eval pm = sum($(Symbol(m,"_π",o))[[1:2;end-1:end]])
            @eval std_plot!($(Symbol("p",1)).layout.grid[$c_m,$c_o],x->$(Symbol(m,"_pdf",o))(x,$i))
            @eval std_plot!($(Symbol("p",1)).layout.grid[$c_m,$c_o],x->true_pdf(x,$i))
            @eval std_plot!($(Symbol("q",1)).layout.grid[$c_m,$c_o],x->$(Symbol(m,"_cdf",o))(x,$i))
            @eval std_plot!($(Symbol("q",1)).layout.grid[$c_m,$c_o],x->true_cdf(x,$i))

        end
        for plt in (:p,:q)
            if c_o==1
                @eval plot!($(Symbol(plt,1)).layout.grid[$c_m,$c_o],ylabel=$(uppercase(string(m))))
                if plt==:p
                    @eval plot!($(Symbol(plt,1)).layout.grid[$c_m,$c_o],yticks=0:0.002:0.006)
                else
                    @eval plot!($(Symbol(plt,1)).layout.grid[$c_m,$c_o],yticks=0:0.01:0.05)
                end
            end
            if c_m==1
                @eval plot!($(Symbol(plt,1)).layout.grid[$c_m,$c_o],title=$(string(o)))
            end
            if c_m==3
                @eval plot!($(Symbol(plt,1)).layout.grid[$c_m,$c_o],xticks=0:5:10)
            end
        end
        @eval (L1_cdf_error_row[$c_m] = log10(DiscretisedFluidQueues.Lp(x->$(Symbol(m,"_cdf",o))(x,1),x->true_cdf(x,1),range(-eps(),10+eps(),length=n_err_evals)) +
                                                        DiscretisedFluidQueues.Lp(x->$(Symbol(m,"_cdf",o))(x,2),x->true_cdf(x,2),range(-eps(),10+eps(),length=n_err_evals))))
        @eval (ks_error_row[$c_m] = log10(DiscretisedFluidQueues.kolmogorov_smirnov(x->$(Symbol(m,"_cdf",o))(x,1),x->true_cdf(x,1),range(-eps(),10+eps(),length=n_err_evals)) +
                                                        DiscretisedFluidQueues.kolmogorov_smirnov(x->$(Symbol(m,"_cdf",o))(x,2),x->true_cdf(x,2),range(-eps(),10+eps(),length=n_err_evals))))
        @eval (L1_pdf_error_row[$c_m] = log10(DiscretisedFluidQueues.Lp(x->$(Symbol(m,"_pdf",o))(x,1),x->true_pdf(x,1),range(-eps(),10+eps(),length=n_err_evals)) +
                                                        DiscretisedFluidQueues.Lp(x->$(Symbol(m,"_pdf",o))(x,2),x->true_pdf(x,2),range(-eps(),10+eps(),length=n_err_evals))))
        @eval (L2_pdf_error_row[$c_m] = log10(DiscretisedFluidQueues.Lp(x->$(Symbol(m,"_pdf",o))(x,1),x->true_pdf(x,1),range(-eps(),10+eps(),length=n_err_evals),2) +
                                                        DiscretisedFluidQueues.Lp(x->$(Symbol(m,"_pdf",o))(x,2),x->true_pdf(x,2),range(-eps(),10+eps(),length=n_err_evals),2)))
    end
    push!(L1_cdf_errors,L1_cdf_error_row)
    push!(ks_errors,ks_error_row)
    push!(L1_pdf_errors,L1_pdf_error_row)
    push!(L2_pdf_errors,L2_pdf_error_row)
end
file = pth*"figs/"
savefig(p1,file*"pdfs.svg")
savefig(q1,file*"cdfs.svg")
titles = [
    "L¹ error between the true CDF and approximations";
    "KS error between the true CDF and approximations";
    "L¹ error between the true PDF and approximations";
    "L² error between the true PDF and approximations";
]
for (c_err,err) in enumerate((:L1_cdf_errors,:ks_errors,:L1_pdf_errors,:L2_pdf_errors))
    file = pth*"data/"*string(err)
    @eval CSV.write($file*".csv",$err)
    @eval q = plot()
    df = @eval $err
    for m in names(df)
        @eval plot!(q,orders,$(err)[:,$(string(m))],label=$m)
    end
    @eval plot!(q,xlabel="Order", ylabel="log₁₀ Error",
        title=titles[$c_err],
        legend=:outertopright)
    #display(q)
    # #display((@__DIR__)*)
    file = pth*"figs/"*string(err)
    @eval savefig(q,$file*".svg")
end


text:
@params mutable struct AssembleK{T} <: AbstractFunction{T}
    problem::StiffnessTopOptProblem
    K::AbstractMatrix{T}
    global_dofs::AbstractVector{<:Integer} # preallocated dof vector for a cell
end

function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::AssembleK)
    return println("TopOpt global linear stiffness matrix assembly function")
end

function AssembleK(problem::StiffnessTopOptProblem)
    dh = problem.ch.dh
    k = ndofs_per_cell(dh)
    global_dofs = zeros(Int, k)
    return AssembleK(problem, initialize_K(problem), global_dofs)
end

"""
    assembleK(Kes) = (Σ_ei x_e * K_ei)

Forward-pass function call.
"""
function (ak::AssembleK{T})(Kes::AbstractVector{<:AbstractMatrix{T}}) where {T}
    @unpack problem, K, global_dofs = ak
    dh = problem.ch.dh
    if K isa Symmetric
        K.data.nzval .= 0
        assembler = Ferrite.AssemblerSparsityPattern(K.data, T[], Int[], Int[])
    else
        K.nzval .= 0
        assembler = Ferrite.AssemblerSparsityPattern(K, T[], Int[], Int[])
    end
    Ke = zeros(T, size(Kes[1]))
    TK = eltype(Kes)
    for (i, _) in enumerate(CellIterator(dh))
        celldofs!(global_dofs, dh, i)
        Ke = TK isa Symmetric ? Kes[i].data : Kes[i]
        Ferrite.assemble!(assembler, global_dofs, Ke)
    end
    return copy(K)
end

"""
    ChainRulesCore.rrule(ak::AssembleK{T}, Kes)

`rrule` for autodiff.

Let's consider composite function `g(F(...))`, where
`F` can be a struct-valued, vector-valued, or matrix-valued function.
In the case here, `F = AssembleK`. Then `rrule` wants us to find `g`'s derivative
w.r.t each *output* of `F`, given `g`'s derivative w.r.t. each *input* of `F`.
Here, `F: K_e -> K_ij = sum_e K_e_ij`. Then `dK_ij/dK_e_ij = 1`.
And we know `Delta_ij = dg/dK_ij`.

And our goal for `rrule` is to go from `dg/dKij` to `dg/dK_e_ij`,
which has the same structure as the input `K_e`.

    dg/dK_e_ij = sum_i'j' dg/dK_i'j' dK_i'j'/dK_e_ij
    		   = dg/dK_ij dK_ij/dK_e_ij
    		   # i, j above are global indices
    		   # i, j below are local indices
    		   = Delta[global_dofs[i], global_dofs[j]]

    (dK_i'j'/dK_e_ij = 0 unless i' == i, j' == j, (i, j) in e)

which can be shortened as:

    dg/dK_e = Delta[global_dofs, global_dofs]
"""
function ChainRulesCore.rrule(
    ak::AssembleK{T}, Kes::AbstractVector{<:AbstractMatrix{T}}
) where {T}
    @unpack problem, K, global_dofs = ak
    dh = problem.ch.dh
    # * forward-pass
    K = ak(Kes)
    n_dofs = length(global_dofs)
    function assembleK_pullback(Δ)
        ΔKes = [zeros(T, n_dofs, n_dofs) for _ in 1:getncells(dh.grid)]
        for (ci, _) in enumerate(CellIterator(dh))
            celldofs!(global_dofs, dh, ci)
            ΔKes[ci] = Δ[global_dofs, global_dofs]
        end
        return Tangent{typeof(ak)}(; problem=NoTangent(), K=Δ, global_dofs=NoTangent()),
        ΔKes
    end
    return K, assembleK_pullback
end


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

using Documenter, SimpleSDMLayers
using GBIF
using Statistics

makedocs(
    sitename = "Simple SDM Layers",
    modules = [SimpleSDMLayers],
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Types" => "man/types.md",
            "Overloads" => "man/overloads.md",
            "Other operations" => "man/operations.md",
            "Data" => "man/data.md",
            "Plotting" => "man/plotting.md"
        ],
        "Examples" => [
            "Temperature data" => "examples/temperature.md",
            "GBIF integration" => "examples/gbif.md",
            "Importing raster data" => "examples/import.md",
            "Sliding window analysis" => "examples/slidingwindow.md",
            "Landcover data" => "examples/landcover.md",
            "Landcover consensus" => "examples/consensus.md"
        ]
    ]
)

run(`find . -type f -size +40M -delete`)

deploydocs(
    repo = "github.com/EcoJulia/SimpleSDMLayers.jl.git",
    push_preview = true
)


text:
# exports
export create, Cone, show

mutable struct Cone <: Shape
    radius::Float64
    length::Float64
end

function Base.show(io::IO, x::Cone)
    println(io, "Cone")
    println(io, "Radius: $(x.radius)")
    println(io, "Length: $(x.length)")
end

function create(c::Cone; resolution=nothing, rand_=0.0, type::Int64=1)
    radius = c.radius
    length_ = c.length
    x, v, y, vol, type_ = unpack(create(Disk(radius, length_), resolution=resolution, rand_=rand_))
    X = @view x[1:2, :]
    z = @view x[3:3, :]

    mask = vec(sum(X.^2, dims=1) .<= ((length_ .- z)./length_.*radius).^2)
    mesh = x[:, mask]
    return Dict(:x => mesh, :v => zeros(size(mesh)), :y => copy(mesh), :volume => vol[mask], :type => type*ones(Int64, sum(mask)))
end

max_stars_repo_name: sajidmannan/PDMesh

max_issues_repo_name: sajidmannan/PDMesh

max_forks_repo_name: sajidmannan/PDMesh


text:
@testset "Linear systems" begin
    @testset "LinearUnivariateTimeInvariant" begin
        y = rand(100)
        Z = ones(3)
        d = rand()
        T = Matrix(Diagonal(ones(3, 3)))
        c = zeros(3)
        R = [
            1.0 0.0
            0.0 1.0
            0.0 0.0
        ]
        H = rand()
        Q = Matrix(Diagonal(rand(2, 2)))
        # Test if it throws an error with an invalid system matrices
        @test_throws MethodError LinearUnivariateTimeInvariant{Float64}(
            y, Z, T, R, d, c[1], H, Q
        )
        @test_throws MethodError LinearUnivariateTimeInvariant{Float64}(
            y, Z, T, R, d, c, H[:, :], Q
        )
        @test_throws MethodError LinearUnivariateTimeInvariant{Float64}(
            y, Z, T, R[1], d, c, H, Q
        )
        @test_throws DimensionMismatch LinearUnivariateTimeInvariant{Float64}(
            y, Z, T[:, 1:2], R, d, c, H, Q
        )
        @test_throws DimensionMismatch LinearUnivariateTimeInvariant{Float64}(
            y, Z[1:2], T, R, d, c, H, Q
        )
        @test_throws DimensionMismatch LinearUnivariateTimeInvariant{Float64}(
            y, Z, T, R, d, c, H, Q[:, 1:1]
        )

        system = StateSpaceModels.LinearUnivariateTimeInvariant{Float64}(
            y, Z, T, R, d, c, H, Q
        )

        n_obs_sim = 200
        initial_state = [0.0; 0.0; 0.7]
        sim = StateSpaceModels.simulate(system, initial_state, n_obs_sim)
        @test length(sim) == 200
        @test maximum(sim) <= 1e5 # This should not explode

        sim, alpha = StateSpaceModels.simulate(
            system, initial_state, n_obs_sim; return_simulated_states=true
        )
        @test size(alpha) == (200, 3)
        # The third state has 0.0 variance and initial_state 0.7 so this should not change
        @test all(alpha[:, 3] .== 0.7)
    end

    @testset "LinearUnivariateTimeVariant" begin
        # TODO
    end

    @testset "LinearMultivariateTimeInvariant" begin
        # TODO
    end

    @testset "LinearMultivariateTimeVariant" begin
        # TODO
    end
end


text:
abstract type DynamicComponent <: DeviceParameter end
abstract type DynamicInjection <: Device end
supports_services(d::DynamicInjection) = false

function get_dynamic_components(device::T) where {T <: DynamicInjection}
    components_set = Set()
    for field in fieldnames(T)
        if fieldtype(T, field) <: DynamicComponent
            push!(components_set, getfield(device, field))
        end
    end
    return components_set
end


text:
@testset "Maps and Reductions" begin
    @testset "sum" begin
        @testset "Vector" begin
            M = 3
            frule_test(sum, (randn(M), randn(M)))
            rrule_test(sum, randn(), (randn(M), randn(M)))
        end
        @testset "Matrix" begin
            M, N = 3, 4
            frule_test(sum, (randn(M, N), randn(M, N)))
            rrule_test(sum, randn(), (randn(M, N), randn(M, N)))
        end
        @testset "Array{T, 3}" begin
            M, N, P = 3, 7, 11
            frule_test(sum, (randn(M, N, P), randn(M, N, P)))
            rrule_test(sum, randn(), (randn(M, N, P), randn(M, N, P)))
        end
        @testset "keyword arguments" begin
            n = 4
            X = randn(n, n+1)
            y, pullback = rrule(sum, X; dims=2)
            ȳ = randn(size(y))
            _, x̄_ad = pullback(ȳ)
            x̄_fd = only(j′vp(central_fdm(5, 1), x->sum(x, dims=2), ȳ, X))
            @test x̄_ad ≈ x̄_fd atol=1e-9 rtol=1e-9
        end
    end  # sum
end


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/libim_image_/include/im_image.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@static Sys.iswindows() ? (const libim_image_ = "im") : (const libim_image_ = "libim")  # Name of IM shared lib.

struct imImage
    width::Cint
    height::Cint
    color_space::Cint
    data_type::Cint
    has_alpha::Cint
    depth::Cint
    line_size::Cint
    plane_size::Cint
    size::Cint
    count::Cint
    data::Ptr{Ptr{Cvoid}}
    palette::Ptr{Clong}
    palette_count::Cint
    attrib_table::Ptr{Cvoid}
end

function imImageCreate(width::Cint, height::Cint, color_space::Cint, data_type::Cint)
  ccall( (:imImageCreate, libim_image_), Ptr{imImage}, (Cint, Cint, Cint, Cint), width, height, color_space, data_type)
end

function imImageInit(width::Cint, height::Cint, color_mode::Cint, data_type::Cint, data_buffer::Ptr{Cvoid}, palette::Ptr{Clong}, palette_count::Cint)
  ccall( (:imImageInit, libim_image_), Ptr{imImage}, (Cint, Cint, Cint, Cint, Ptr{Cvoid}, Ptr{Clong}, Cint), width, height, color_mode, data_type, data_buffer, palette, palette_count)
end

function imImageCreateBased(image::Ptr{imImage}, width::Cint, height::Cint, color_space::Cint, data_type::Cint)
  ccall( (:imImageCreateBased, libim_image_), Ptr{imImage}, (Ptr{imImage}, Cint, Cint, Cint, Cint), image, width, height, color_space, data_type)
end

function imImageDestroy(image::Ptr{imImage})
  ccall( (:imImageDestroy, libim_image_), Cvoid, (Ptr{imImage},), image)
end

function imImageAddAlpha(image::Ptr{imImage})
  ccall( (:imImageAddAlpha, libim_image_), Cvoid, (Ptr{imImage},), image)
end

function imImageSetAlpha(image::Ptr{imImage}, alpha::Cfloat)
  ccall( (:imImageSetAlpha, libim_image_), Cvoid, (Ptr{imImage}, Cfloat), image, alpha)
end

function imImageRemoveAlpha(image::Ptr{imImage})
  ccall( (:imImageRemoveAlpha, libim_image_), Cvoid, (Ptr{imImage},), image)
end

function imImageReshape(image::Ptr{imImage}, width::Cint, height::Cint)
  ccall( (:imImageReshape, libim_image_), Cvoid, (Ptr{imImage}, Cint, Cint), image, width, height)
end

function imImageCopy(src_image::Ptr{imImage}, dst_image::Ptr{imImage})
  ccall( (:imImageCopy, libim_image_), Cvoid, (Ptr{imImage}, Ptr{imImage}), src_image, dst_image)
end

function imImageCopyData(src_image::Ptr{imImage}, dst_image::Ptr{imImage})
  ccall( (:imImageCopyData, libim_image_), Cvoid, (Ptr{imImage}, Ptr{imImage}), src_image, dst_image)
end

function imImageCopyAttributes(src_image::Ptr{imImage}, dst_image::Ptr{imImage})
  ccall( (:imImageCopyAttributes, libim_image_), Cvoid, (Ptr{imImage}, Ptr{imImage}), src_image, dst_image)
end

function imImageMergeAttributes(src_image::Ptr{imImage}, dst_image::Ptr{imImage})
  ccall( (:imImageMergeAttributes, libim_image_), Cvoid, (Ptr{imImage}, Ptr{imImage}), src_image, dst_image)
end

function imImageCopyPlane(src_image::Ptr{imImage}, src_plane::Cint, dst_image::Ptr{imImage}, dst_plane::Cint)
  ccall( (:imImageCopyPlane, libim_image_), Cvoid, (Ptr{imImage}, Cint, Ptr{imImage}, Cint), src_image, src_plane, dst_image, dst_plane)
end

function imImageDuplicate(image::Ptr{imImage})
  ccall( (:imImageDuplicate, libim_image_), Ptr{imImage}, (Ptr{imImage},), image)
end

function imImageClone(image::Ptr{imImage})
  ccall( (:imImageClone, libim_image_), Ptr{imImage}, (Ptr{imImage},), image)
end

function imImageSetAttribute(image::Ptr{imImage}, attrib::String, data_type::Cint, count::Cint, data::Ptr{Cvoid})
  ccall( (:imImageSetAttribute, libim_image_), Cvoid, (Ptr{imImage}, Ptr{UInt8}, Cint, Cint, Ptr{Cvoid}), image, attrib, data_type, count, data)
end

function imImageGetAttribute(image::Ptr{imImage}, attrib::String, data_type::Ptr{Cint}, count::Ptr{Cint})
  ccall( (:imImageGetAttribute, libim_image_), Ptr{Cvoid}, (Ptr{imImage}, Ptr{UInt8}, Ptr{Cint}, Ptr{Cint}), image, attrib, data_type, count)
end

function imImageGetAttribInteger(image::Ptr{imImage}, attrib::String, index::Integer)
  ccall( (:imImageGetAttribInteger, libim_image_), Ptr{Cvoid}, (Ptr{imImage}, Ptr{UInt8}, Ptr{Cint}), image, attrib, index)
end

function imImageGetAttributeList(image::Ptr{imImage}, attrib::Ptr{String}, attrib_count::Ptr{Cint})
  ccall( (:imImageGetAttributeList, libim_image_), Cvoid, (Ptr{imImage}, Ptr{String}, Ptr{Cint}), image, attrib, attrib_count)
end

function imImageGetAttribReal(image::Ptr{imImage}, attrib::String, index::Integer)
  ccall( (:imImageGetAttribReal, libim_image_), Ptr{Cvoid}, (Ptr{imImage}, Ptr{UInt8}, Ptr{Cint}), image, attrib, index)
end

function imImageGetAttribString(image::Ptr{imImage}, attrib::Ptr{UInt8})
  ccall( (:imImageGetAttribString, libim_image_), Ptr{UInt8}, (Ptr{imImage}, Ptr{UInt8}), image, attrib)
end

function imImageGetAttributeList(image::Ptr{imImage}, attrib::Ptr{Ptr{UInt8}}, attrib_count::Ptr{Int})
  ccall( (:imImageGetAttributeList, libim_image_), Cvoid, (Ptr{imImage}, Ptr{Ptr{UInt8}}, Cint), image, attrib, attrib_count)
end

function imImageClear(image::Ptr{imImage})
  ccall( (:imImageClear, libim_image_), Cvoid, (Ptr{imImage},), image)
end

function imImageIsBitmap(image::Ptr{imImage})
  ccall( (:imImageIsBitmap, libim_image_), Cint, (Ptr{imImage},), image)
end

function imImageSetPalette(image::Ptr{imImage}, palette::Ptr{Clong}, palette_count::Cint)
  ccall( (:imImageSetPalette, libim_image_), Cvoid, (Ptr{imImage}, Ptr{Clong}, Cint), image, palette, palette_count)
end

function imImageMatchSize(image1::Ptr{imImage}, image2::Ptr{imImage})
  ccall( (:imImageMatchSize, libim_image_), Cint, (Ptr{imImage}, Ptr{imImage}), image1, image2)
end

function imImageMatchColor(image1::Ptr{imImage}, image2::Ptr{imImage})
  ccall( (:imImageMatchColor, libim_image_), Cint, (Ptr{imImage}, Ptr{imImage}), image1, image2)
end

function imImageMatchDataType(image1::Ptr{imImage}, image2::Ptr{imImage})
  ccall( (:imImageMatchDataType, libim_image_), Cint, (Ptr{imImage}, Ptr{imImage}), image1, image2)
end

function imImageMatchColorSpace(image1::Ptr{imImage}, image2::Ptr{imImage})
  ccall( (:imImageMatchColorSpace, libim_image_), Cint, (Ptr{imImage}, Ptr{imImage}), image1, image2)
end

function imImageMatch(image1::Ptr{imImage}, image2::Ptr{imImage})
  ccall( (:imImageMatch, libim_image_), Cint, (Ptr{imImage}, Ptr{imImage}), image1, image2)
end

function imImageSetMap(image::Ptr{imImage})
  ccall( (:imImageSetMap, libim_image_), Cvoid, (Ptr{imImage},), image)
end

function imImageSetBinary(image::Ptr{imImage})
  ccall( (:imImageSetBinary, libim_image_), Cvoid, (Ptr{imImage},), image)
end

function imImageSetGray(image::Ptr{imImage})
  ccall( (:imImageSetGray, libim_image_), Cvoid, (Ptr{imImage},), image)
end

function imImageMakeBinary(image::Ptr{imImage})
  ccall( (:imImageMakeBinary, libim_image_), Cvoid, (Ptr{imImage},), image)
end

function imImageMakeGray(image::Ptr{imImage})
  ccall( (:imImageMakeGray, libim_image_), Cvoid, (Ptr{imImage},), image)
end

function imFileLoadImage(ifile::Ptr{Cint}, index::Cint, error::Ptr{Cint})
  ccall( (:imFileLoadImage, libim_image_), Ptr{imImage}, (Ptr{Cint}, Cint, Ptr{Cint}), ifile, index, error)
end

function imFileLoadImageFrame(ifile::Ptr{Cint}, index::Cint, image::Ptr{imImage}, error::Ptr{Cint})
  ccall( (:imFileLoadImageFrame, libim_image_), Cvoid, (Ptr{Cint}, Cint, Ptr{imImage}, Ptr{Cint}), ifile, index, image, error)
end

function imFileLoadBitmap(ifile::Ptr{Cint}, index::Cint, error::Ptr{Cint})
  ccall( (:imFileLoadBitmap, libim_image_), Ptr{imImage}, (Ptr{Cint}, Cint, Ptr{Cint}), ifile, index, error)
end

function imFileLoadImageRegion(ifile::Ptr{Cint}, index::Cint, bitmap::Cint, error::Ptr{Cint}, xmin::Cint, xmax::Cint, ymin::Cint, ymax::Cint, width::Cint, height::Cint)
  ccall( (:imFileLoadImageRegion, libim_image_), Ptr{imImage}, (Ptr{Cint}, Cint, Cint, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint, Cint), ifile, index, bitmap, error, xmin, xmax, ymin, ymax, width, height)
end

function imFileLoadBitmapFrame(ifile::Ptr{Cint}, index::Cint, image::Ptr{imImage}, error::Ptr{Cint})
  ccall( (:imFileLoadBitmapFrame, libim_image_), Cvoid, (Ptr{Cint}, Cint, Ptr{imImage}, Ptr{Cint}), ifile, index, image, error)
end

function imFileSaveImage(ifile::Ptr{Cint}, image::Ptr{imImage})
  ccall( (:imFileSaveImage, libim_image_), Cint, (Ptr{Cint}, Ptr{imImage}), ifile, image)
end

function imFileImageLoad(file_name::String, index::Cint, error::Ptr{Cint})
  ccall( (:imFileImageLoad, libim_image_), Ptr{imImage}, (Ptr{UInt8}, Cint, Ptr{Cint}), file_name, index, error)
end

function imFileImageLoadBitmap(file_name::String, index::Cint, error::Ptr{Cint})
  ccall( (:imFileImageLoadBitmap, libim_image_), Ptr{imImage}, (Ptr{UInt8}, Cint, Ptr{Cint}), file_name, index, error)
end

function imFileImageLoadRegion(file_name::String, index::Cint, bitmap::Cint, error::Ptr{Cint}, xmin::Cint, xmax::Cint, ymin::Cint, ymax::Cint, width::Cint, height::Cint)
  ccall( (:imFileImageLoadRegion, libim_image_), Ptr{imImage}, (Ptr{UInt8}, Cint, Cint, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint, Cint), file_name, index, bitmap, error, xmin, xmax, ymin, ymax, width, height)
end

function imFileImageSave(file_name::String, format::String, image::Ptr{imImage})
  ccall( (:imFileImageSave, libim_image_), Cint, (Ptr{UInt8}, Ptr{UInt8}, Ptr{imImage}), file_name, format, image)
end


text:
using Random

include("../tlwe.jl")

function tlwe_test()
    n = 630
    α = 2e-15
    μ = 50000000
    rng = RandomDevice()

    num_try_key = 10000
    num_try_whole = 100

    for i in 1:num_try_whole
        key = LweKey(rng, n)
        for j in 1:num_try_key
            message = rand(rng, Bool)
            ciphertext = lwe_encrypt_bool(rng, message, key, α, μ)
            decrypted_text = lwe_decrypt_bool(ciphertext, key)
            if message != decrypted_text
                return false
            end
        end
    end

    return true
end

tlwe_test()

max_stars_repo_name: bbtatt/BBTTFHE

max_issues_repo_name: bbtatt/BBTTFHE

max_forks_repo_name: bbtatt/BBTTFHE


text:
module Math

export Vector3, Point, len, unit, dot, randominunitsphere, squaredlength, cross, lenhost, unithost, Transform
export identitytransform, translation, rotation, permute, scale

using CUDAnative
using CthulhuVision.Random

struct Vector3
    x::Float32
    y::Float32
    z::Float32
end

@inline Base.:+(a::Vector3, b::Vector3) :: Vector3 = Vector3(a.x + b.x, a.y + b.y, a.z + b.z)
@inline Base.:-(a::Vector3, b::Vector3) :: Vector3 = Vector3(a.x - b.x, a.y - b.y, a.z - b.z)
@inline Base.:-(a::Vector3) :: Vector3 = Vector3(-a.x, -a.y, -a.z)
@inline Base.:*(t::Float32, v::Vector3) :: Vector3 = Vector3(t*v.x, t*v.y, t*v.z)
@inline Base.:/(a::Vector3, t::Float32) :: Vector3 = Vector3(a.x / t, a.y / t, a.z / t)
@inline len(a::Vector3) = CUDAnative.sqrt(a.x*a.x + a.y*a.y + a.z*a.z)
@inline lenhost(a::Vector3) = sqrt(a.x*a.x + a.y*a.y + a.z*a.z)
@inline squaredlength(a::Vector3) :: Float32 = a.x*a.x + a.y*a.y + a.z*a.z

@inline dot(a::Vector3, b::Vector3) :: Float32 = a.x * b.x + a.y * b.y + a.z * b.z

@inline function cross(a::Vector3, b::Vector3) :: Vector3
    Vector3(
          a.y * b.z - a.z * b.y,
        -(a.x * b.z - a.z * b.x),
          a.x * b.y - a.y * b.x
    )
end

@inline function unit(a::Vector3) :: Vector3
    l = len(a)
    a / l
end

@inline function unithost(a::Vector3) :: Vector3
    l = lenhost(a)
    a / l
end

@inline function byindex(v::Vector3, index) :: Float32
    if index == 0
        v.x
    elseif index == 1
        v.y
    else
        v.z
    end
end

@inline permute(v::Vector3, x, y, z) :: Vector3 = Vector3(byindex(v, x), byindex(v, y), byindex(v, z))

@inline function randominunitsphere(rng::UniformRNG) :: Vector3
    p = Vector3(1.0f0, 1.0f0, 1.0f0)

    while squaredlength(p) >= 1.0f0
        p = 2.0f0 * Vector3(next(rng), next(rng), next(rng)) - Vector3(1f0, 1f0, 1f0)
    end

    p
end

#
# Point
#

struct Point
    x::Float32
    y::Float32
    z::Float32
end

@inline function byindex(p::Point, index) :: Float32
    if index == 0
        p.x
    elseif index == 1
        p.y
    else
        p.z
    end
end

@inline Base.:-(p::Point, v::Vector3) :: Point = Point(p.x - v.x, p.y - v.y, p.z - v.z)
@inline Base.:-(a::Point, b::Point) :: Vector3 = Vector3(a.x - b.x, a.y - b.y, a.z - b.z)
@inline Base.:+(a::Point, b::Vector3) :: Point = Point(a.x + b.x, a.y + b.y, a.z + b.z)
@inline Base.abs(p::Point) :: Point = Point(abs(p.x), abs(p.y), abs(p.z))
@inline permute(p::Point, x, y, z) :: Point = Point(byindex(p, x), byindex(p, y), byindex(p, z))
@inline dot(a::Vector3, b::Point) :: Float32 = a.x * b.x + a.y * b.y + a.z * b.z

#
# Transformation matrix 4x4
#

struct Transform
    e11::Float32
    e12::Float32
    e13::Float32
    e14::Float32
    e21::Float32
    e22::Float32
    e23::Float32
    e24::Float32
    e31::Float32
    e32::Float32
    e33::Float32
    e34::Float32
    e41::Float32
    e42::Float32
    e43::Float32
    e44::Float32
end

@inline function Base.:*(a::Transform, b::Transform) :: Transform
    e11 = a.e11 * b.e11 + a.e12 * b.e21 + a.e13 * b.e31 + a.e14 * b.e41
    e12 = a.e11 * b.e12 + a.e12 * b.e22 + a.e13 * b.e32 + a.e14 * b.e42
    e13 = a.e11 * b.e13 + a.e12 * b.e23 + a.e13 * b.e33 + a.e14 * b.e43
    e14 = a.e11 * b.e14 + a.e12 * b.e24 + a.e13 * b.e34 + a.e14 * b.e44

    e21 = a.e21 * b.e11 + a.e22 * b.e21 + a.e23 * b.e31 + a.e24 * b.e41
    e22 = a.e21 * b.e12 + a.e22 * b.e22 + a.e23 * b.e32 + a.e24 * b.e42
    e23 = a.e21 * b.e13 + a.e22 * b.e23 + a.e23 * b.e33 + a.e24 * b.e43
    e24 = a.e21 * b.e14 + a.e22 * b.e24 + a.e23 * b.e34 + a.e24 * b.e44

    e31 = a.e31 * b.e11 + a.e32 * b.e21 + a.e33 * b.e31 + a.e34 * b.e41
    e32 = a.e31 * b.e12 + a.e32 * b.e22 + a.e33 * b.e32 + a.e34 * b.e42
    e33 = a.e31 * b.e13 + a.e32 * b.e23 + a.e33 * b.e33 + a.e34 * b.e43
    e34 = a.e31 * b.e14 + a.e32 * b.e24 + a.e33 * b.e34 + a.e34 * b.e44

    e41 = a.e41 * b.e11 + a.e42 * b.e21 + a.e43 * b.e31 + a.e44 * b.e41
    e42 = a.e41 * b.e12 + a.e42 * b.e22 + a.e43 * b.e32 + a.e44 * b.e42
    e43 = a.e41 * b.e13 + a.e42 * b.e23 + a.e43 * b.e33 + a.e44 * b.e43
    e44 = a.e41 * b.e14 + a.e42 * b.e24 + a.e43 * b.e34 + a.e44 * b.e44

    Transform(
        e11, e12, e13, e14,
        e21, e22, e23, e24,
        e31, e32, e33, e34,
        e41, e42, e43, e44,
    )
end

@inline function Base.:*(a::Transform, b::Vector3) :: Vector3
    w = 1f0
    x = a.e11 * b.x + a.e12 * b.y + a.e13 * b.z + a.e14 * w
    y = a.e21 * b.x + a.e22 * b.y + a.e23 * b.z + a.e24 * w
    z = a.e31 * b.x + a.e32 * b.y + a.e33 * b.z + a.e34 * w

    Vector3(x, y, z)
end

@inline function Base.:*(a::Transform, b::Point) :: Point
    w = 1f0
    x = a.e11 * b.x + a.e12 * b.y + a.e13 * b.z + a.e14 * w
    y = a.e21 * b.x + a.e22 * b.y + a.e23 * b.z + a.e24 * w
    z = a.e31 * b.x + a.e32 * b.y + a.e33 * b.z + a.e34 * w

    Point(x, y, z)
end

########
# Host #
########

function identitytransform() :: Transform
    Transform(
        1f0, 0f0, 0f0, 0f0,
        0f0, 1f0, 0f0, 0f0,
        0f0, 0f0, 1f0, 0f0,
        0f0, 0f0, 0f0, 1f0,
    )
end

function translation(v::Vector3) :: Transform
    Transform(
        1f0, 0f0, 0f0, v.x,
        0f0, 1f0, 0f0, v.y,
        0f0, 0f0, 1f0, v.z,
        0f0, 0f0, 0f0, 0f0,
    )
end

function rotation(θ::Float32, axis::Vector3) :: Transform
    u = unithost(axis)

    e11 = cos(θ) + u.x*u.x*(1f0 - cos(θ))
    e12 = u.x * u.y * (1f0 - cos(θ)) - u.z * sin(θ)
    e13 = u.x * u.z * (1f0 - cos(θ)) + u.y * sin(θ)
    e14 = 0f0

    e21 = u.y * u.x * (1f0 - cos(θ)) + u.z * sin(θ)
    e22 = cos(θ) + u.y * u.y * (1f0 - cos(θ))
    e23 = u.y * u.z * (1f0 - cos(θ)) - u.x * sin(θ)
    e24 = 0f0

    e31 = u.z * u.x * (1f0 - cos(θ)) - u.y * sin(θ)
    e32 = u.z * u.y * (1f0 - cos(θ)) + u.x * sin(θ)
    e33 = cos(θ) + u.z * u.z * (1f0 - cos(θ))
    e34 = 0f0

    e41 = 0f0
    e42 = 0f0
    e43 = 0f0
    e44 = 1f0

    Transform(
        e11, e12, e13, e14,
        e21, e22, e23, e24,
        e31, e32, e33, e34,
        e41, e42, e43, e44,
    )
end

function scale(v::Vector3) :: Transform
    Transform(
        v.x, 0f0, 0f0, 0f0,
        0f0, v.y, 0f0, 0f0,
        0f0, 0f0, v.z, 0f0,
        0f0, 0f0, 0f0, 1f0,
    )
end

end


text:
function visualize!(vis, model,
    q_body, q_f1, q_f2, q_f3, q_f4;
    q_body_ref = q_body,
    ref_tl = 0.5,
    r = 0.025,
    l_torso = 0.19,
    w_torso = 0.1025,
    Δt = 0.1)

	default_background!(vis)

	setobject!(vis["torso"],
    	Rect(Vec(-l_torso, -w_torso, -0.025),Vec(2.0 * l_torso, 2.0 * w_torso, 0.05)),
    	MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

    setobject!(vis["torso_ref"],
    	Rect(Vec(-l_torso, -w_torso, -0.025),Vec(2.0 * l_torso, 2.0 * w_torso, 0.05)),
    	MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, ref_tl)))

	feet1 = setobject!(vis["feet1"], Sphere(Point3f0(0),
        convert(Float32, r)),
        MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

	feet2 = setobject!(vis["feet2"], Sphere(Point3f0(0),
		convert(Float32, r)),
		MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

	feet3 = setobject!(vis["feet3"], Sphere(Point3f0(0),
		convert(Float32, r)),
		MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

	feet4 = setobject!(vis["feet4"], Sphere(Point3f0(0),
        convert(Float32, r)),
        MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

	anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

	T = length(q_body)
	p_shift = [0.0, 0.0, r]

	for t = 1:T
		MeshCat.atframe(anim, t) do
			rot = MRP(q_body[t][4:6]...)
            rot_ref = MRP(q_body_ref[t][4:6]...)

			p_torso = q_body[t][1:3] + p_shift
            p_torso_ref = q_body_ref[t][1:3] + p_shift
			p_foot1 = q_f1[t] + p_shift
			p_foot2 = q_f2[t] + p_shift
			p_foot3 = q_f3[t] + p_shift
			p_foot4 = q_f4[t] + p_shift

			settransform!(vis["torso"], compose(Translation(p_torso), LinearMap(rot)))
            settransform!(vis["torso_ref"], compose(Translation(p_torso_ref), LinearMap(rot_ref)))
			settransform!(vis["feet1"], Translation(p_foot1))
			settransform!(vis["feet2"], Translation(p_foot2))
			settransform!(vis["feet3"], Translation(p_foot3))
			settransform!(vis["feet4"], Translation(p_foot4))
		end
	end

	# settransform!(vis["/Cameras/default"],
	#     compose(Translation(0.0, 0.0, -1.0), LinearMap(RotZ(-pi / 2.0))))

	MeshCat.setanimation!(vis, anim)
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using MichiganAutonomousVehicles
using NLOptControl
using VehicleModels
using Parameters
using DataFrames
using PrettyPlots, Plots

case_name = "s1"
c = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/planner/","RTPP",".yaml")))
c["vehicle"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/vehicle/","hmmwv",".yaml")))
c["goal"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/case/",case_name,".yaml")))["goal"]
c["X0"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/case/",case_name,".yaml")))["X0"]
c["obstacle"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/case/",case_name,".yaml")))["obstacle"]
setConfig(c,"misc";(:N=>15),(:model=>:ThreeDOFv3),(:solver=>:KNITRO),(:integrationScheme=>:trapezoidal))
setConfig(c,"X0";(:ux=>5.0))# the type needs to be a float
# otherwise, this will pop up: WARNING: Ipopt finished with status Invalid_Number_Detected
fixYAML(c)   # fix messed up data types
n = initializeAutonomousControl(c);
simMPC!(n;updateFunction=updateAutoParams!,checkFunction=checkCrash)

function mpcSim(n,idx)
  mpcPlot(n,idx)
end
# (:mpc_lines=>[(6.5,:blue,:solid)])
# (:mpc_markers =>(:circle,:blueviolet,0.0,0.0)),
a = 2
plotSettings(;(:X0p=>true), (:mpc_markers =>(:circle,:blueviolet,0.0,1.)),(:polyPts=>true),(:simulate=>true),(:plant=>true),(:plantOnly=>false),(:size=>(a*900,a*600)),(:format=>"png"));
anim = @animate for idx in 1:5#length(n.r.ocp.dfs)
  mpcSim(n,idx)
end
cd(n.r.resultsDir)
  gif(anim,"mainSim.gif",fps=Int(ceil(1/n.mpc.v.tex)));
  run(`ffmpeg -f gif -i mainSim.gif RESULT.mp4`)
cd(n.r.mainDir)

savePlantData!(n)
saveOptData(n)


text:
@testset "Reification" begin
    mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Int}()))
    config = MOIT.Config()

    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> (MOIU.mock_optimize!(mock, [0, 1])),
        (mock::MOIU.MockOptimizer) -> (MOIU.mock_optimize!(mock, [0, 1])),
    )
    COIDT.reificationtest(mock, config)
end


text:
module PandorasBoxPandorasBox

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/pandorasBox" PandorasBox(x::Integer, y::Integer, completeChapter::Bool=false, dialog::String="")

const placements = Ahorn.PlacementDict(
    "Pandora's Box (Pandora's Box)" => Ahorn.EntityPlacement(
        PandorasBox
    )
)

sprite = "objects/pandorasBox/pandorasBox/box_idle0"

function Ahorn.selection(entity::PandorasBox)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y, jx=0.5, jy=1.0)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::PandorasBox, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0, jx=0.5, jy=1.0)

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
module Ex5_10
using FaultDetectionTools
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

# Example 5.10 - Solution of an EFDIP
println("Example 5.10")

# enter output and fault vector dimensions
p = 3; mf = 3;
# generate random dimensions for system order and input vectors
nu = Int(floor(1+4*rand())); mu = Int(floor(1+4*rand()));
nd = Int(floor(1+4*rand())); md = Int(floor(1+4*rand()));

# define random Gu(s) and Gd(s) with triplex sensor redundancy
# and Gf(s) for triplex sensor faults
Gu = ones(3,1)*rss(nu,1,mu); # enter Gu(s) in state-space form
Gd = ones(3,1)*rss(nd,1,md); # enter Gd(s) in state-space form
Gf = eye(3);                 # enter Gf(s) for sensor faults
atol = 1.e-7;                # tolerance for rank tests

# build model with faults
#sysf = [Gu Gd Gf]
sysf = fdimodset([Gu Gd Gf], c = 1:mu, d = mu.+(1:md), f = (mu+md).+(1:mf))

S = [ 0 1 1; 1 0 1; 1 1 0] .> 0;  # enter structure matrix

## Step 1) of Procedure EFDI
# compute Q1, the left nullspace of [Gu Gd;I 0], and  Rf1 = Q1*[Gf;0]
# QRf contains Q1 in QRf[:,1:p+mu] and Rf1 in QRf[:,p+mu+1:p+mu+mf]
QRf, info = glnull([sysf.sys; eye(mu,mu+md+mf)],mf; atol)

## Step 2): of Procedure EFDI
# initialize overall filter Q and fault detection system Rf
nb = size(S,1);      # number of filters
Qt = similar(Vector{typeof(QRf)},nb)
Rft = similar(Vector{typeof(QRf)},nb)
for i = 1:nb
    # Step 2.1): define  Rf11 = Rf1[:,indd], Rf12 = Rf1[:,indf]
    indd = Vector(1:mf)[S[i,:] .== false]
    indf = Vector(1:mf)[S[i,:] .== true]
    # pack [ Rf11 Rf12 [Q1 Rf1] ]
    sysc = fdimodset(QRf, d = (p+mu) .+ indd, f = (p+mu) .+ indf, aux = Vector(1:p+mu+mf))
    # Step 2.2): apply  Procedure EFD to {0,Rf11,Rf12,[Q1 Rf1]}
    # to determine a least order Q1i such that Q1i*Rf11 = 0
    # QRfauxi contains: [Q1i*Rf12 Q1i*Q1 Q1i*Rf1]
    _, QRfauxi, infoi = efdsyn(sysc; rdim = 1, atol)
    # extract [Q1i*Q1 Q1i*Rf1]
    QRfi = QRfauxi.sys[:,QRfauxi.aux]
    # Step 2.3): extract Q[i] = Q1i*Q1 and Rf[i] = Q1i*Rf1
    Qt[i] = QRfi[:,1:p+mu]
    Rft[i] = QRfi[:,(p+mu).+(1:mf)]
end

# normalize Q and Rf to match example
scale = sign.([ Rft[1].D[1,2], Rft[2].D[1,3], Rft[3].D[1,1]])
Q = FDIFilter(scale .* Qt, p, mu)
R = FDIFilterIF(scale .* Rft, 0, 0, mf)
println("gap = $(fdif2ngap(R,S)[1])")
end # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
# Plotting functionalities from Matplotlib.

using PyPlot

export plot, pcolormesh, pcolormeshslice, vdfslice, streamplot, quiver, plotmesh

@static if matplotlib.__version__ >= "3.3"
   matplotlib.rc("image", cmap="turbo") # set default colormap
end

matplotlib.rc("font", size=14)
matplotlib.rc("xtick", labelsize=10)
matplotlib.rc("ytick", labelsize=10)

"""
    plot(meta, var, ax=nothing; kwargs)

Plot `var` from `meta` of 1D VLSV data. If `ax===nothing`, plot on the current active axes.
"""
function PyPlot.plot(meta::MetaVLSV, var, ax=nothing; kwargs...)
   data = readvariable(meta, var)

   x = LinRange(meta.coordmin[1], meta.coordmax[1], meta.ncells[1])

   if isnothing(ax) ax = plt.gca() end

   ax.plot(x, data; kwargs...)
end

"""
    streamplot(meta, var, ax=nothing; comp="xy", axisunit=RE, kwargs...)

Wrapper over Matplotlib's streamplot function. The `comp` option can take a subset of "xyz"
in any order. `axisunit` can be chosen from `RE, SI`.
The keyword arguments can be any valid Matplotlib arguments into streamplot.

# Optional arguments
- `comp`: a subset of "xyz" in any order.
- `axisunit`: chosen from RE and SI.
"""
function PyPlot.streamplot(meta::MetaVLSV, var::AbstractString, ax=nothing;
   comp="xy", axisunit=RE, kwargs...)

   X, Y, v1, v2 = set_vector(meta, var, comp, axisunit)

   if isnothing(ax) ax = plt.gca() end

   streamplot(X, Y, v1, v2; kwargs...)
end

"""
    quiver(meta, var, ax=nothing; comp="xy", axisunit=RE, stride=10, kwargs...)

Wrapper over Matplotlib's quiver function. If `ax===nothing`, plot on the current active
axes. The `comp` option can take a subset of "xyz" in any order. `axisunit` can be chosen
from `RE, SI`.
The keyword arguments can be any valid Matplotlib arguments into quiver.

# Optional arguments
- `comp`: a subset of "xyz" in any order.
- `axisunit`: chosen from RE and SI.
- `stride::Integer`: arrow strides in number of cells.
"""
function PyPlot.quiver(meta::MetaVLSV, var::AbstractString, ax=nothing;
   comp="xy", axisunit::AxisUnit=RE, stride::Integer=10, kwargs...)

   X, Y, v1, v2 = set_vector(meta, var, comp, axisunit)

   if isnothing(ax) ax = plt.gca() end

   Xq  = @view X[1:stride:end, 1:stride:end]
   Yq  = @view Y[1:stride:end, 1:stride:end]
   v1q = @view v1[1:stride:end, 1:stride:end]
   v2q = @view v2[1:stride:end, 1:stride:end]

   ax.quiver(Xq, Yq, v1q, v2q; kwargs...)
end

function set_vector(meta::MetaVLSV, var, comp, axisunit::AxisUnit)
   @unpack ncells, coordmin, coordmax = meta
   if occursin("x", comp)
      v1_ = 1
      if occursin("y", comp)
         v2_ = 2
         sizes = [ncells[1], ncells[2]]
         plotrange = [coordmin[1], coordmax[1], coordmin[2], coordmax[2]]
      else
         v2_ = 3
         sizes = [ncells[1], ncells[3]]
         plotrange = [coordmin[1], coordmax[1], coordmin[3], coordmax[3]]
      end
   else
      v1_, v2_ = 2, 3
      sizes = [ncells[2], ncells[3]]
      plotrange = [coordmin[2], coordmax[2], coordmin[3], coordmax[3]]
   end

   data = readvariable(meta, var)

   if !startswith(var, "fg_") # vlasov grid
      @assert ndims(data) == 2 && size(data,1) == 3 "Vector data required!"
      data = reshape(data, 3, sizes[1], sizes[2])
   end

   v1 = data[v1_,:,:]
   v2 = data[v2_,:,:]

   x, y = get_axis(axisunit, plotrange, sizes)

   # meshgrid: note the array ordering difference between Julia and Python!
   X = [i for _ in y, i in x]
   Y = [j for j in y, _ in x]

   X, Y, v1', v2'
end

"""
    pcolormesh(meta::MetaVLSV, var::AbstractString, ax=nothing;
       op=:mag, axisunit=RE, colorscale=Linear, vmin=-Inf, vmax=Inf, addcolorbar=true,
       kwargs...)

Plot a variable using pseudocolor from 2D VLSV data.
If `ax` is provided, then it will plot on that axes.
If 3D or AMR grid detected, it will pass arguments to [`pcolormeshslice`](@ref).

# Optional arguments
- `op::Symbol`: the component of a vector, chosen from `:mag, :x, :y, :z, :1, :2, :3`.
- `axisunit::AxisUnit`: the unit of axis ∈ `RE, SI`.
- `colorscale::ColorScale`: Linear, Log, or SymLog.
- `vmin::Float`: minimum data range. Set to maximum of data if not specified.
- `vmax::Float`: maximum data range. Set to minimum of data if not specified.
- `addcolorbar::Bool`: whether to add a colorbar to the colormesh.

`pcolormesh(meta, var)`

`pcolormesh(meta, var, axisunit=SI)`

`pcolormesh(data, func, colorscale=Log)`
"""
function PyPlot.pcolormesh(meta::MetaVLSV, var::AbstractString, ax=nothing; op=:mag,
   axisunit::AxisUnit=RE, colorscale::ColorScale=Linear, addcolorbar=true,
   vmin=-Inf, vmax=Inf, kwargs...)

   if ndims(meta) == 3 || meta.maxamr > 0
      # check if origin and normal exist in kwargs
      normal = haskey(kwargs, :normal) ? kwargs.data.normal : :y
      origin = haskey(kwargs, :origin) ? kwargs.data.origin : 0.0
      kwargs = Base.structdiff(kwargs.data, (normal = normal, origin = origin))
      c = pcolormeshslice(meta, var, ax; op, axisunit, colorscale, addcolorbar, vmin, vmax,
         normal, origin, kwargs...)
      return c
   end

   pArgs = set_args(meta, var, axisunit)

   x, y = get_axis(pArgs)
   data = prep2d(meta, var, op)'

   if var in ("fg_b", "fg_e", "vg_b_vol", "vg_e_vol") || endswith(var, "vg_v")
      rho_ = findfirst(endswith("rho"), meta.variable)
      if !isnothing(rho_)
         rho = readvariable(meta, meta.variable[rho_])
         rho = reshape(rho, pArgs.sizes[1], pArgs.sizes[2])
         mask = findall(==(0.0), rho')

         if ndims(data) == 2
            @inbounds data[mask] .= NaN
         else
            ind = CartesianIndices((pArgs.sizes[2], pArgs.sizes[1]))
            for m in mask
               @inbounds data[:, ind[m][1], ind[m][2]] .= NaN
            end
         end
      end
   end

   cnorm, cticks = set_colorbar(colorscale, vmin, vmax, data)

   if isnothing(ax) ax = plt.gca() end

   if colorscale != SymLog
      c = ax.pcolormesh(x, y, data; norm=cnorm, shading="nearest", kwargs...)
   else
      c = ax.pcolormesh(x, y, data; norm=cnorm, cmap=matplotlib.cm.RdBu, shading="nearest",
         kwargs...)
   end

   set_plot(c, ax, pArgs, cticks, addcolorbar)

   return c
end

"""
    pcolormeshslice(meta, var, ax=nothing; kwargs...)

Plot pseudocolor var on a 2D slice of 3D vlsv data.
If `ax` is provided, then it will plot on that axes.
It would be easier to call [`pcolormesh`](@ref), since it auto-detects dimension.

# Optional arguments
- `op::Symbol`: the component of a vector, chosen from `:mag, :x, :y, :z, :1, :2, :3`.
- `origin::Float`: center of slice plane in the normal direction.
- `normal::Symbol`: the normal direction of cut plane, chosen from `:x, :y, :z`.
- `axisunit::AxisUnit`: the unit of axis ∈ `RE, SI`.
- `colorscale::ColorScale`: color scale for data ∈ (`Linear`, `Log`, `SymLog`)
- `vmin::Real`: minimum data range. Set to maximum of data if not specified.
- `vmax::Real`: maximum data range. Set to minimum of data if not specified.
- `addcolorbar::Bool`: whether to add a colorbar to the colormesh.

`pcolormeshslice(meta, var)`

`pcolormeshslice(meta, var, op=:z, origin=1.0, normal=:x)`

`pcolormeshslice(data, func, colorscale=Log)`
"""
function pcolormeshslice(meta::MetaVLSV, var::AbstractString, ax=nothing; op::Symbol=:mag,
   origin=0.0, normal::Symbol=:y, axisunit::AxisUnit=RE, colorscale::ColorScale=Linear,
   addcolorbar=true, vmin::Real=-Inf, vmax::Real=Inf, kwargs...)

   pArgs = set_args(meta, var, axisunit; normal, origin)

   data = prep2dslice(meta, var, normal, op, pArgs)'
   x, y = get_axis(pArgs)

   cnorm, cticks = set_colorbar(colorscale, vmin, vmax, data)

   if isnothing(ax) ax = plt.gca() end

   c = ax.pcolormesh(x, y, data; norm=cnorm, shading="nearest", kwargs...)

   set_plot(c, ax, pArgs, cticks, addcolorbar)

   c
end

"Set colorbar norm and ticks in a given range `v1` to `v2` for `data` in `colorscale`."
function set_colorbar(colorscale::ColorScale, v1, v2, data=[1.0])
   vmin, vmax = set_lim(v1, v2, data, colorscale)
   if colorscale == Linear
      levels = matplotlib.ticker.MaxNLocator(nbins=255).tick_values(vmin, vmax)
      cnorm = matplotlib.colors.BoundaryNorm(levels, ncolors=256, clip=true)
      ticks = matplotlib.ticker.LinearLocator(numticks=9)
   elseif colorscale == Log # logarithmic
      cnorm = matplotlib.colors.LogNorm(;vmin, vmax)
      ticks = matplotlib.ticker.LogLocator(base=10, subs=collect(0:9))
   else # symmetric log
      linthresh = 1.0
      logstep = 1

      logthresh = floor(Int, log10(linthresh))
      minlog = ceil(Int, log10(-vmin))
      maxlog = ceil(Int, log10(vmax))

      cnorm = matplotlib.colors.SymLogNorm(;linthresh, linscale=0.03, vmin, vmax, base=10)
      ticks = [ [-(10.0^x) for x in minlog:-logstep:logthresh]..., 0.0,
         [10.0^x for x in logthresh+1:logstep:maxlog]..., ]
   end

   cnorm, ticks
end

"Configure customized plot."
function set_plot(c, ax, pArgs::PlotArgs, cticks, addcolorbar)
   @unpack str_title, strx, stry, cb_title = pArgs

   if addcolorbar
      cb = colorbar(c; ax, ticks=cticks, fraction=0.04, pad=0.02)
      !isempty(cb_title) && cb.ax.set_ylabel(cb_title)
      cb.ax.tick_params(direction="in")
   end

   ax.set_title(str_title, fontweight="bold")
   ax.set_xlabel(strx)
   ax.set_ylabel(stry)
   ax.set_aspect("equal")

   # Set border line widths
   for loc in ("left", "bottom", "right", "top")
      edge = get(ax.spines, loc, nothing)
      edge.set_linewidth(2.0)
   end

   ax.xaxis.set_tick_params(width=2.0, length=3)
   ax.yaxis.set_tick_params(width=2.0, length=3)
   return
end

"""
    vdfslice(meta, location, ax=nothing; kwargs...)

Plot the 2D slice cut of phase space distribution function at `location` within velocity
range `limits`. If `ax===nothing`, plot on the current active axes.
# Optional arguments
- `unit::AxisUnit`: location unit in `SI`, `RE`.
- `unitv::String`: velocity unit in ("km/s", "m/s").
- `limits::Vector{Real}`: velocity space range given in [xmin, xmax, ymin, ymax].
- `slicetype`: symbol for choosing the slice type from :xy, :xz, :yz, :bperp, :bpar, :bpar1.
- `center`: symbol for setting the reference frame from :bulk, :peak.
- `vslicethick`: setting the velocity space slice thickness in the normal direction. If set
to 0, the whole distribution along the normal direction is projected onto a plane. Currently
this is only meaningful when `center` is set such that a range near the bulk/peak normal
velocity is selected!
- `fmin, fmax`: minimum and maximum VDF values for plotting.
- `weight::Symbol`: choosing distribution weights from phase space density or particle flux
between `:particle` and `:flux`.
- `flimit`: minimum VDF threshold for plotting.
- `kwargs...`: any valid keyword argument for hist2d.
"""
function vdfslice(meta::MetaVLSV, location, ax=nothing; limits=[-Inf, Inf, -Inf, Inf],
   verbose=false, species="proton", fmin=-Inf, fmax=Inf, unit=SI, unitv="km/s",
   slicetype=:nothing, vslicethick=0.0, center=:nothing, weight=:particle, flimit=-1.0,
   kwargs...)

   v1, v2, r1, r2, fweight, strx, stry, str_title =
      prep_vdf(meta, location;
         species, unit, unitv, slicetype, vslicethick, center, weight, flimit, verbose)

   isinf(fmin) && (fmin = minimum(fweight))
   isinf(fmax) && (fmax = maximum(fweight))

   verbose && @info "Active f range is $fmin, $fmax"

   if isnothing(ax) ax = plt.gca() end

   cnorm = matplotlib.colors.LogNorm(vmin=fmin, vmax=fmax)

   h = ax.hist2d(v1, v2, bins=(r1, r2), weights=fweight, norm=cnorm)

   ax.set_title(str_title, fontweight="bold")
   ax.set_xlabel(strx, weight="black")
   ax.set_ylabel(stry, weight="black")
   ax.set_aspect("equal")
   ax.grid(color="grey", linestyle="-")
   ax.tick_params(direction="in")

   cb = colorbar(h[4]; ax, fraction=0.046, pad=0.02)
   cb.ax.tick_params(direction="in")
   cb_title = cb.ax.set_ylabel("f(v)")

   if slicetype in (:bperp, :bpar, :bpar1)
      # Draw vector of magnetic field direction
   end
   plt.tight_layout()

   h[4] # h[1] is 2D data, h[2] is x axis, h[3] is y axis
end

"""
    plotmesh(meta; projection="3d", origin=0.0, marker="+", kwargs...)

Plot mesh cell centers from axis view `projection`. `projection` should be either "3d", "x",
"y" or "z". `origin` is center of projection plane in the normal direction.
"""
function plotmesh(meta::MetaVLSV, ax=nothing; projection="3d", origin=0.0, marker="+",
   kwargs...)
   @unpack coordmin, coordmax, cellid = meta
   if projection in ("x", "y", "z")
      dirp, dir1, dir2 =
         if projection == "x"
            1, 2, 3
         elseif projection == "y"
            2, 1, 3
         else # "z"
            3, 1, 2
         end
      sliceoffset = origin - coordmin[dirp]
      ids, _ = getslicecell(meta, sliceoffset, dirp, coordmin[dirp], coordmax[dirp])
   else # 3D
      ids = cellid
   end

   centers = [zeros(SVector{3, Float32}) for _ in ids]
   for (i, id) in enumerate(ids)
      @inbounds centers[i] = getcellcoordinates(meta, id)
   end

   if isnothing(ax) ax = plt.gca() end

   if projection in ("x", "y", "z")
      x1 = getindex.(centers, dir1)
      x2 = getindex.(centers, dir2)
      s = ax.scatter(x1, x2; marker, kwargs...)
   else
      x1 = getindex.(centers, 1)
      x2 = getindex.(centers, 2)
      x3 = getindex.(centers, 3)
      s = ax.scatter(x1, x2, x3; marker, kwargs...)
   end
   s
end


text:
@testset "Transforms" begin
  # compositions
  for c in [Composition(1,2,3),
            Composition(0,2,3),
            Composition(1,2,0)]
    @test alrinv(alr(c)) == c
    @test clrinv(clr(c)) == c
    @test ilrinv(ilr(c)) == c
  end

  # convenience methods for heap-allocated vectors
  @test alrinv([1,2,3]) isa Composition{4}
  @test clrinv([1,2,3]) isa Composition{3}
  @test ilrinv([1,2,3]) isa Composition{4}

  # TableTransform.jl API
  t = (a=[1.,0.,1.], b=[2.,2.,2.], c=[3.,3.,0.])
  n, c = apply(ALR(), t)
  talr = revert(ALR(), n, c)
  n, c = apply(CLR(), t)
  tclr = revert(CLR(), n, c)
  n, c = apply(ILR(), t)
  tilr = revert(ILR(), n, c)
  @test Tables.matrix(talr) ≈ Tables.matrix(tclr)
  @test Tables.matrix(tclr) ≈ Tables.matrix(tilr)
  @test Tables.matrix(talr) ≈ Tables.matrix(tilr)

  n, c = apply(ALR(:b), t)
  talr = revert(ALR(:b), n, c)
  @test n    |> Tables.columnnames |> collect == [:a, :c]
  @test talr |> Tables.columnnames |> collect == [:a, :b, :c]

  n, c = apply(ILR(:b), t)
  tilr = revert(ILR(:b), n, c)
  @test n    |> Tables.columnnames |> collect == [:a, :c]
  @test tilr |> Tables.columnnames |> collect == [:a, :b, :c]

  # Tests for Closure
  t = (a=[2.,66.,0.], b=[4.,22.,2.], c=[4.,12.,98.])
  n, c = apply(Closure(), t)
  tcls = revert(Closure(), n, c)
  @test Tables.matrix(n) ≈ [0.2 0.4 0.4; 0.66 0.22 0.12; 0.00 0.02 0.98;]
  @test Tables.matrix(tcls) ≈ Tables.matrix(t)

  # Tests for Remainder
  t = (a=[2.,66.,0.], b=[4.,22.,2.], c=[4.,12.,98.])
  n, c = apply(Remainder(), t)
  trem = revert(Remainder(), n, c)
  Xt = Tables.matrix(t)
  Xn = Tables.matrix(n)
  @test Xn[:, 1:end-1] == Xt
  @test all(x -> 0 ≤ x ≤ c, Xn[:, end])
  @test n    |> Tables.columnnames |> collect == [:a, :b, :c, :remainder]
  @test trem |> Tables.columnnames |> collect == [:a, :b, :c]

  t = (a=[1.,10.,0.], b=[1.,5.,0.], c=[4.,2.,1.])
  n = reapply(Remainder(), t, c)
  Xn = Tables.matrix(n)
  @test all(x -> 0 ≤ x ≤ c, Xn[:, end])

  t = (a=[1.,10.,0.], b=[1.,5.,0.], remainder=[4.,2.,1.])
  names = t |> Remainder() |> Tables.columnnames |> collect
  @test names == [:a, :b, :remainder, :remainder_]

  t = (a=[1.,10.,0.], b=[1.,5.,0.], remainder=[4.,2.,1.])
  n1, c1 = apply(Remainder(), t)
  n2 = reapply(Remainder(), t, c1)
  @test n1 == n2
end


text:
using ArcadeLearningEnvironment
using Test

# TODO: find a public domain ROM and test on that

ale = ALE_new()

# This value is set in ale.cfg
@test getFloat(ale, "repeat_action_probability") == 0.0

ALE_del(ale)


text:
include("../src/Iteration.jl")
include("testing.jl")

importall Iteration

print("testing iteration ")

@assert_equal({[1, 2], [3, 4]}, [e for e = slicerator([1, 2, 3, 4], 2)])
@assert_equal({[1, 2], [3, 4], [5, 6]}, [e for e = slicerator([1, 2, 3, 4, 5, 6], 2)])
@assert_equal({[1, 2, 3], [4, 5, 6]}, [e for e = slicerator([1, 2, 3, 4, 5, 6], 3)])

@assert_equal({[3, 4], [5, 6]}, [e for e = offset(1, slicerator([1, 2, 3, 4, 5, 6], 2))])
@assert_equal({[5, 6]}, [e for e = offset(2, slicerator([1, 2, 3, 4, 5, 6], 2))])
@assert_equal({}, [e for e = offset(3, slicerator([1, 2, 3, 4, 5, 6], 2))])
@assert_equal({}, [e for e = offset(4, slicerator([1, 2, 3, 4, 5, 6], 2))])

println()
println()


text:
function _second_order_sensitivities(loss,prob,alg,sensealg::ForwardDiffOverAdjoint,
                                     args...;kwargs...)
   ForwardDiff.jacobian(prob.p) do p
     x = Zygote.gradient(p) do _p
       loss(solve(prob,alg,args...;p=_p,sensealg=sensealg.adjalg,kwargs...))
     end
     first(x)
   end
end

function _second_order_sensitivity_product(loss,v,prob,alg,sensealg::ForwardDiffOverAdjoint,
                                           args...;kwargs...)

   θ = ForwardDiff.Dual.(prob.p,v)
   _loss = p -> loss(solve(prob,alg,args...;p=p,sensealg=sensealg.adjalg,kwargs...))
   getindex.(ForwardDiff.partials.(Zygote.gradient(_loss,θ)[1]),1)
end