text:
a2q(A,X) = giac(:a2q,A,X)
abcuv(a,b,c) = giac(:abcuv,a,b,c)
abcuv(a,b,c,var) = giac(:abcuv,a,b,c,var)
adjoint_matrix(A) = giac(:adjoint_matrix, A)

cot(x::giac) = giac(:cot, x)
sec(x::giac) = giac(:sec, x)
csc(x::giac) = giac(:csc, x)
#sech(x::giac) = giac(:sech, x) # not available in Giac
#csch(x::giac) = giac(:csch, x) # not available in Giac
asec(x::giac) = giac(:asec, x)
acsc(x::giac) = giac(:acsc, x)
acot(x::giac) = giac(:acot, x)

sign(x::giac) = giac(:sign,x)
log(::Irrational{:e}, x::giac) = log(x)
log(base, x::giac)  = giac(:logb, x, base)
log10(x::giac) = log(10, x)
round(x::giac) = giac(:round, x)
ceil(x::giac) = giac(:ceil, x)
floor(x::giac) = giac(:floor, x)
trunc(x::giac) = giac(:trunc, x)
factorial(x::giac) = giac(:factorial, x)
binomial(x::giac, y::giac) = giac(:binomial, x, y)
binomial(x::giac, y) = giac(:binomial, x, y)
binomial(x, y::giac) = giac(:binomial, x, y)
expm1(x::giac) = giac(:EXPM1, x)
log1p(x::giac) = giac(:LNP1, x)
Ei(x::giac) = giac(:Ei, x)
sinint(x::giac) = giac(:Si, x)
cosint(x::giac) = giac(:Ci, x)
Psi(x::giac) = giac(:Psi, x)
Psi(x::giac, y) = giac(:Psi, x, y)
Psi(x, y::giac) = giac(:Psi, x, y)
digamma(x::giac) = Psi(x)
polygamma(n, x::giac) = Psi(x, n)
polygamma(n::giac, x) = Psi(x, n)

erf(x::giac) = giac(:erf, x)
erfc(x::giac) = giac(:erfc, x)
gamma(x::giac) = giac(:Gamma, x) #  Giac/Gamma -> Julia/gamma
beta(x::giac, y::giac) = giac(:Beta, x, y) # Giac/Beta -> Julia/beta
beta(x::giac, y) = giac(:Beta, x, y) # Giac/Beta -> Julia/beta
beta(x, y::giac) = giac(:Beta, x, y) # Giac/Beta -> Julia/beta
zeta(x::giac) = giac(:Zeta, x)
airyai(x::giac) = giac(:Airy_Ai, x) # Giac/Airy_Ai does not work as expected
airybi(x::giac) = giac(:Airy_Bi, x)  # Giac/Airy_Bi does not work as expected
besselj(nu,x::giac) = giac(:BesselJ, nu, x) # Giac/BesselJ -> Julia/besselj
bessely(nu,x::giac) = giac(:BesselY, nu, x) # Giac/BesselY -> Julia/bessely

evaluate(x) = x
evaluatef(x) = x
simplify(x) = x

evaluate(ex::giac) = giac(:eval, ex)
evaluatef(ex::giac) = giac(:evalf, ex)
evaluatef(ex::giac, n::Int) = giac(:evalf, ex, n)
evalf = evaluatef
simplify(ex::giac) = giac(:simplify, ex)
factor(ex::giac) = giac(:factor, ex)
expand(ex::giac) = giac(:expand, ex)
collect(ex::giac) = giac(:collect, ex)
#collect(exs::Array{giac_SYMB}) = giac(:collect, giac(exs))
collect(exs::Array{giac}) = giac(:collect, giac(exs))
factor(ex::giac, sr)= giac(:factor, ex, sr)
diff(ex::giac, vars...) = giac(:derive, ex, vars...)
diff(ex::giac, vars::Array) = giac(:derive, ex, vars)
partfrac(ex) = giac(:partfrac, ex)
partfrac(ex, var) = giac(:partfrac, ex, var)
equal(left, right) = giac(:equal, left, right)
⩦(a,b) = equal(a, b)
right(eq) = giac(:right, eq)
left(eq) = giac(:left, eq)
denominator(ex::giac) = giac(:denom, ex) #Giac/denom -> Julia/denominator
numerator(ex::giac) = giac(:numer, ex) #Giac/numer -> Julia/numerator

subst(ex, var, val) = giac(:subst, ex, var, val)
subst(ex, sp::Pair...) = subst(ex, [v for (v,w) in sp], [w for (v,w) in sp])

limit(ex, var, val) = giac(:limit, ex, var, val)
limit(ex, var_val::Pair) = giac(:limit, ex, var_val[1], var_val[2])
limit(ex, var, val, dir) = giac(:limit, ex, var, val, dir)
limit(ex, var_val::Pair, dir) = giac(:limit, ex, var_val[1], var_val[2], dir)
series(ex, var, val) = giac(:series, ex, var, val)
series(ex, var, val, order) = giac(:series, ex, var, val, order)
series(ex, var_val::Pair) = giac(:series, ex, var_val[1], var_val[2])
series(ex, var_val::Pair, order) = giac(:series, ex, var_val[1], var_val[2], order)
series(ex, var, val, order, dir) = giac(:series, ex, var, val, order, dir)
series(ex, var_val::Pair, order, dir) = giac(:series, ex, var_val[1], var_val[2], order, dir)
sum(ex::giac) = giac(:sum, ex)
sum(ex::giac, x::giac) = giac(:sum, ex, x)
sum(ex::giac, x::giac, a::Union{giac,Number}, b::Union{giac,Number}) = giac(:sum, ex, x, a, b)
sum(ex::giac, x::giac, a::Union{giac,Number}, b::Union{giac,Number}, step::Union{giac,Number}) =
    giac(:sum, ex, x, a, b, step)
curl(exs, vars) = giac(:curl, exs, vars)
divergence(exs, vars) = giac(:divergence, exs, vars)
grad(ex, vars) = giac(:grad, ex, vars)
hessian(ex, vars) = giac(:hessian, ex, vars)
solve(eqs, var) = giac(:solve, eqs, var)
preval(ex, a, b) = giac(:preval, ex, a, b)
preval(ex, a, b, var) = giac(:preval, ex, a, b, var)
sum_riemann(expr, vars) = giac(:sum_riemann, expr, vars)
taylor(ex, var) = giac(:taylor, ex, var)
taylor(ex, var, order) = giac(:taylor, ex, var, order)
taylor(ex, var_val::Pair) = giac(:taylor, ex, var_val[1], var_val[2])
taylor(ex, var_val::Pair, order) = giac(:taylor, ex, var_val[1], order, var_val[2])
divpc(p1, p2, order) = giac(:divpc, p1, p2, order)
solve(eq) = giac(:solve, eq)
cSolve(eq) = giac(:cSolve, eq)
cSolve(eqs, var) = giac(:cSolve, eqs, var)
zeros(ex::giac) = giac(:zeros, ex)
zeros(exs, var) = giac(:zeros, exs, var)
cZeros(eq::giac) = giac(:cZeros, eq)
cZeros(exs, var) = giac(:cZeros, exs, var)
fSolve(eqs, vars, guess) = giac(:fsolve, eqs, vars, guess)
fSolve(eqs, var_guess::Pair) = giac(:fsolve, eqs, var_guess[1], var_guess[2])
fSolve(eqs, vars) = giac(:fsolve, eqs, vars)
fSolve(eqs) = giac(:fsolve, eqs)
deSolve(eq, var) = giac(:desolve, eq, var)
deSolve(eq, timevar, var) = giac(:desolve, eq, timevar, var)
linsolve(eqs, vars) = giac(:linsolve, eqs, vars)

logcollect(ex) = giac(:lncollect, ex)
powexpand(ex) = giac(:powexpand, ex)
texpand(ex) = giac(:texpand, ex)
exp2pow(ex) = giac(:exp2pow, ex)
pow2exp(ex) = giac(:pow2exp, ex)
exp2trig(ex) = giac(:exp2trig, ex)
expexpand(ex) = giac(:expexpand, ex)
asin2acos(ex) = giac(:asin2acos, ex)
asin2atan(ex) = giac(:asin2atan, ex)
sin2costan(ex) = giac(:sin2costan, ex)
acos2asin(ex) = giac(:acos2asin, ex)
acos2atan(ex) = giac(:acos2atan, ex)
cos2sintan(ex) = giac(:cos2sintan, ex)
atan2asin(ex) = giac(:atan2asin, ex)
atan2acos(ex) = giac(:atan2acos, ex)
tan2sincos(ex) = giac(:tan2sincos, ex)
halftan(ex) = giac(:halftan, ex)
trigsin(ex) = giac(:trigsin, ex)
trigcos(ex) = giac(:trigcos, ex)
trigtan(ex) = giac(:trigtan, ex)
atrig2log(ex) = giac(:atrig2ln, ex)
tlin(ex) = giac(:tlin, ex)
tcollect(ex) = giac(:tcollect, ex)
trigexpand(ex) = giac(:trigexpand, ex)
trig2exp(ex) = giac(:trig2exp, ex)

gbasis(polys, vars) = giac(:gbasis, polys, vars)
greduce(poly, basis, vars) = giac(:greduce, poly, basis, vars)

integrate(ex::giac) = giac(:integrate, ex)
integrate(ex::giac, var::giac) = giac(:integrate, ex, var)
integrate(ex::giac, var::giac, a::Union{giac,Number}, b::Union{giac,Number}) =
    giac(:integrate, ex, var, a, b)

head(ex::giac) = giac(:sommet, ex)
args(ex::giac) = giac(:feuille, ex)

# to be continued ...


text:
export BlockShuffle, CycleShuffle, CircShift

#########################################################################
# BlockSuffle
#########################################################################
function get_uniform_blocklengths(L::Int, n::Int)
    # Compute block lengths
    N = floor(Int, L/n)
    R = L % n
    blocklengths = [N for i = 1:n]
    for i = 1:R
        blocklengths[i] += 1
    end
    return blocklengths
end

"""
    BlockShuffle(n::Int; shift = false)

A block shuffle surrogate constructed by dividing the time series
into `n` blocks of roughly equal width at random indices (end
blocks are wrapped around to the start of the time series).

If `shift` is `true`, then the input signal is circularly shifted by a
random number of steps prior to picking blocks.

Block shuffle surrogates roughly preserve short-range temporal properties
in the time series (e.g. correlations at lags less than the block length),
but break any long-term dynamical information (e.g. correlations beyond
the block length).

Hence, these surrogates can be used to test any null hypothesis aimed at
comparing short-range dynamical properties versus long-range dynamical
properties of the signal.
"""
struct BlockShuffle{I <: Integer, B <: Bool} <: Surrogate
    n::I
    shift::B

    function BlockShuffle(n::I; shift::B = false) where {I <: Integer, B <: Bool}
        return new{I, B}(n, shift)
    end

    # Split time series into 10 pieces by default. Shifting disabled by default.
    function BlockShuffle()
        return new{Int, Bool}(10, false)
    end
end

function surrogenerator(x::AbstractVector, bs::BlockShuffle, rng = Random.default_rng())
    bs.n < length(x) || error("The number of blocks exceeds number of available points")

    # The lengths of the blocks (one block will have a differing length if length of
    # time series in not a multiple of the number of blocks, so )
    blocklengths = get_uniform_blocklengths(length(x), bs.n)

    # The start index of each block.
    startinds = vcat(1, cumsum(blocklengths) .+ 1)

    # The data from which we will sample. This array may be circularly shifted.
    x_rotated = copy(x)

    # The order in which we will draw the blocks. Will be shuffled every time
    # a new surrogate is generated.
    draw_order = collect(1:bs.n)

    init = (
        blocklengths = blocklengths,
        startinds = startinds,
        x_rotated = x_rotated,
        draw_order = draw_order,
    )

    # The surrogate.
    s = similar(x)

    return SurrogateGenerator(bs, x, s, init, rng)
end

function (sg::SurrogateGenerator{<:BlockShuffle})()
    init_fields = (:blocklengths, :startinds, :x_rotated, :draw_order)
    blocklengths, startinds, x_rotated, draw_order = getfield.(Ref(sg.init), init_fields)
    x, s, n = sg.x, sg.s, sg.method.n

    # Circular shift, if desired
    if sg.method.shift
        circshift!(x_rotated, x, rand(sg.rng, 1:length(x)))
    end

    # Shuffle blocks
    shuffle!(sg.rng, draw_order)

    k = 1
    npts_sampled = 0
    for i in draw_order
        # The index of the data point which this block starts (in `x_rotated`)
        sᵢ = startinds[i]

        # The length of this block.
        l = blocklengths[i]

        # Indices of the block which is sampled.
        ix_from = sᵢ:(sᵢ + l - 1)

        # Indices in `s` into which this block will be placed.
        ix_into = (npts_sampled + 1):(npts_sampled + 1 + l - 1)

        # Do the assignment in-place behind a code barrier. This about
        # 3x as efficients as doing an elementwise assignment and
        # doesn't allocate.
        assign_block_to_surrogate!(s, x_rotated, ix_into, ix_from)

        npts_sampled += l
    end
    #@assert all(sort(s) .== sort(x_rotated))
    return s
end

function assign_block_to_surrogate!(s, x_rotated, ix_into, ix_from)
    @assert length(ix_into) == length(ix_from)
    n = length(ix_into)
    @inbounds for k = 1:n
        ito = ix_into[k]
        ifr = ix_from[k]
        s[ito] = x_rotated[ifr]
    end
end

#########################################################################
# CycleShuffle
#########################################################################
"""
    CycleShuffle(n::Int = 7, σ = 0.5)

Cycle shuffled surrogates[^Theiler1994] that identify successive local peaks in the data and shuffle the
cycles in-between the peaks. Similar to [`BlockShuffle`](@ref), but here
the "blocks" are defined as follows:
1. The timeseries is smoothened via convolution with a Gaussian (`DSP.gaussian(n, σ)`).
2. Local maxima of the smoothened signal define the peaks, and thus the blocks in between them.
3. The first and last index of timeseries can never be peaks and thus signals that
   should have peaks very close to start or end of the timeseries may not perform well. In addition,
   points before the first or after the last peak are never shuffled.
3. The defined blocks are randomly shuffled as in [`BlockShuffle`](@ref).

CSS are used to test the null hypothesis that the signal is generated by a periodic
oscillator with no dynamical correlation between cycles,
i.e. the evolution of cycles is not deterministic.

See also [`PseudoPeriodic`](@ref).

[^Theiler1994]: J. Theiler, On the evidence for low-dimensional chaos in an epileptic electroencephalogram, [Phys. Lett. A 196](https://doi.org/10.1016/0375-9601(94)00856-K)
"""
struct CycleShuffle{T <: AbstractFloat} <: Surrogate
    n::Int
    σ::T
end
CycleShuffle(n = 7, σ = 0.5) = CycleShuffle{typeof(σ)}(n, σ)

function surrogenerator(x::AbstractVector, cs::CycleShuffle, rng = Random.default_rng())
    n, N = cs.n, length(x)
    g = DSP.gaussian(n, cs.σ)
    smooth = DSP.conv(x, g)
    r = length(smooth) - N
    smooth = iseven(r) ? smooth[r÷2+1:end-r÷2] : smooth[r÷2+1:end-r÷2-1]
    peaks = findall(i -> smooth[i-1] < smooth[i] && smooth[i] > smooth[i+1], 2:N-1)
    blocks = [collect(peaks[i]:peaks[i+1]-1) for i in 1:length(peaks)-1]
    init =  (blocks = blocks, peak1 = peaks[1])
    s = copy(x)
    SurrogateGenerator(cs, x, s, init, rng)
end

function (sg::SurrogateGenerator{<:CycleShuffle})()
    blocks, peak1 = sg.init
    x = sg.x
    s = sg.s
    shuffle!(sg.rng, blocks)
    i = peak1
    for b in blocks
        s[(0:length(b)-1) .+ i] .= @view x[b]
        i += length(b)
    end
    return s
end

#########################################################################
# Timeshift
#########################################################################
random_shift(n::Integer, rng) = n
random_shift(n::AbstractVector{<:Integer}, rng) = rand(rng, n)

"""
    CircShift(n)
Surrogates that are circularly shifted versions of the original timeseries.

`n` can be an integer (the surrogate is the original time series shifted
by `n` indices), or any vector of integers, which which means that each
surrogate is shifted by an integer selected randomly among the entries in `n`.
"""
struct CircShift{N} <: Surrogate
    n::N
end

function surrogenerator(x, sd::CircShift, rng = Random.default_rng())
    return SurrogateGenerator(sd, x, similar(x), nothing, rng)
end

function (sg::SurrogateGenerator{<:CircShift})()
    x, s = sg.x, sg.s
    shift = random_shift(sg.method.n, sg.rng)
    circshift!(s, x, shift)

    return s
end

ext: jl


text:
using StochasticProcesses

bm = BrownianMotion()
t = linspace(0, 1, 10000)
sim(bm, t, 100)

println("sim(BrownianMotion) 10000x10000:")
for i in 1:5
  @time sim(bm, t, 10000)
end


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

@with_kw mutable struct QuasiGradientData{T <: AbstractFloat}
    Q::T = 1e10
    master_objective::AffineDecisionFunction{T} = zero(AffineDecisionFunction{T})
    no_objective::Bool = false
    iterations::Int = 1
end

@with_kw mutable struct QuasiGradientParameters{T <: AbstractFloat}
    debug::Bool = false
    time_limit::T = Inf
    log::Bool = true
    keep::Bool = true
    offset::Int = 0
    indent::Int = 0
end

"""
    QuasiGradientAlgorithm

Functor object for the quasi-gradient algorithm.

...
# Algorithm parameters
- `log::Bool = true`: Specifices if quasi-gradient procedure should be logged on standard output or not.
...
"""
struct QuasiGradientAlgorithm{T <: AbstractFloat,
                              A <: AbstractVector,
                              ST <: StageDecompositionStructure,
                              M <: MOI.AbstractOptimizer,
                              E <: AbstractQuasiGradientExecution,
                              P <: AbstractProximal,
                              S <: AbstractStep,
                              C <: AbstractTerminationCriterion,
                              Pr <: AbstractProgress} <: AbstractQuasiGradient
    structure::ST
    num_subproblems::Int
    data::QuasiGradientData{T}
    parameters::QuasiGradientParameters{T}

    # Master
    master::M
    decisions::DecisionMap
    x::A
    ξ::A
    c::A
    gradient::A
    Q_history::A

    # Execution
    execution::E

    # Prox
    prox::P

    # Step
    step::S

    # Termination
    criterion::C

    progress::Pr

    function QuasiGradientAlgorithm(structure::StageDecompositionStructure,
                                    x₀::AbstractVector,
                                    _execution::AbstractExecution,
                                    subproblems::AbstractSubProblemState,
                                    _prox::AbstractProx,
                                    step::AbstractStepSize,
                                    _criterion::AbstractTermination; kw...)
        # Sanity checks
        length(x₀) != num_decisions(structure, 1) && error("Incorrect length of starting guess, has ", length(x₀), " should be ", num_decisions(structure, 1))
        n = num_subproblems(structure, 2)
        n == 0 && error("No subproblems in stochastic program. Cannot run L-shaped procedure.")
        # Float types
        T = promote_type(eltype(x₀), Float32)
        x₀_ = convert(AbstractVector{T}, copy(x₀))
        A = typeof(x₀_)
        # Structure
        ST = typeof(structure)
        M = typeof(backend(structure.first_stage))
        # Prox
        prox = _prox(structure, x₀_, T)
        P = typeof(prox)
        # Step
        stepsize = step(T)
        S = typeof(stepsize)
        # Termination
        criterion = _criterion(T)
        C = typeof(criterion)
        # Execution policy
        execution = _execution(structure, x₀_, subproblems, T)
        E = typeof(execution)
        # Algorithm parameters
        params = QuasiGradientParameters{T}(; kw...)
        # Progress
        progress = Progress(criterion, "$(indentstr(params.indent))Quasi-gradient Progress ")
        Pr = typeof(progress)

        quasigradient = new{T,A,ST,M,E,P,S,C,Pr}(structure,
                                                 n,
                                                 QuasiGradientData{T}(),
                                                 params,
                                                 backend(structure.first_stage),
                                                 structure.decisions[1],
                                                 x₀_,
                                                 copy(x₀_),
                                                 zero(x₀_),
                                                 zero(x₀_),
                                                 A(),
                                                 execution,
                                                 prox,
                                                 stepsize,
                                                 criterion,
                                                 progress)
        # Initialize solver
        initialize!(quasigradient)
        return quasigradient
    end
end

function show(io::IO, quasigradient::QuasiGradientAlgorithm)
    println(io, typeof(quasigradient).name.name)
    println(io, "State:")
    show(io, quasigradient.data)
    println(io, "Parameters:")
    show(io, quasigradient.parameters)
end

function show(io::IO, ::MIME"text/plain", quasigradient::QuasiGradientAlgorithm)
    show(io, quasigradient)
end

function (quasigradient::QuasiGradientAlgorithm)()
    # Reset timer
    quasigradient.progress.tinit = quasigradient.progress.tlast = time()
    # Start workers (if any)
    start_workers!(quasigradient)
    # Start procedure
    while true
        status = iterate!(quasigradient)
        if status !== nothing
            close_workers!(quasigradient)
            return status
        end
    end
end


text:
module SurveyAnalysis
include("SurveyDesign.jl")

end


text:
using Test
using DataFrames, CSV
using MultiDimEquations
using IndexedTables

println("*** Start testing of MultiDimEquations...")

# *** NEW TEST
# Checking empty variable definition
(exp,consumption)   = defVars(["product","regions","years"],[String,String,Int64],n=2)
(exp2,consumption2) = defVars([3,2,5],n=2,missingValue=missing)
consumption10 = defVars(["product","regions","years"],[String,String,Int64])
consumption20 = defVars([3,2,5])
@test consumption10 == consumption
@test Base.size(consumption20) == Base.size(consumption2)
exp["banana","Canada",2010] = 2
consumption["banana","Canada",2010] = 4
@test exp["banana","Canada",2010] == 2 # testing they are different memory locations
@test exp["banana","Canada",2010] + 1 == 3
consumption2[2,1,5] = 2
@test consumption2[2,1,5] + 1 == 3
exp2[2,1,5] = 10
@test consumption2[2,1,5] == 2

# *** NEW TEST
# checking loading variables from long-format dataframe
df = CSV.read(IOBuffer("""
reg prod var value
us banana production 10
us banana transfCoef    0.6
us apples production    7
us apples transfCoef    0.7
us apples trValues      5
us juice production     missing
us juice transfCoef     missing
us juice trValues       missing
eu banana production    5
eu banana transfCoef    0.7
eu apples production    8
eu apples transfCoef    0.8
eu apples trValues 4
eu juice production missing
eu juice transfCoef missing
eu juice trValues missing
"""), DataFrame, delim=" ", ignorerepeated=true, copycols=true, missingstring="missing") # missing eu banana trValues      1 and us banana trValues      2

prodDf      = df[df.var .== "production",:]
prodSingle  = defLoadVar(prodDf,["reg","prod"], valueCol="value",sparse=true)
prodSingle2 = defLoadVar(prodDf,["reg","prod"], valueCol="value",sparse=false)

(production,transfCoef,trValues)     = defLoadVars(["production","transfCoef","trValues"], df,["reg","prod"], varNameCol="var", valueCol="value",sparse=true)
(production2,transfCoef2,trValues2)  = defLoadVars(["production","transfCoef","trValues"], df,["reg","prod"], varNameCol="var", valueCol="value",sparse=false,fullKeys=true)
trValues3  = defLoadVars(["trValues"], df,["reg","prod"], varNameCol="var", valueCol="value",sparse=false,fullKeys=false)
@test Base.size(trValues2) == (2,3)
@test Base.size(trValues3) == (2,2)

reg = unique(df.reg)
products = unique(df.prod)

@test trValues["eu","apples"] == 4
@test trValues2[findfirst(isequal("eu"),reg),findfirst(isequal("apples"),products)] == 4

# *** NEW TEST
# checking the @meq macro

primPr    = products[1:2]
secPr     = [products[3]]
primPrIdx = [1,2]
secPrIdx  = [3]
nReg      = length(reg)

# Refilling the intentionally left out values...
trValues["eu","banana"] = 1
trValues["us","banana"] = 2
trValues2[findfirst(isequal("eu"),reg),findfirst(isequal("banana"),products)] = 1
trValues2[findfirst(isequal("us"),reg),findfirst(isequal("banana"),products)] = 2

# Checking further computations..
[production[r,sp]   =  sum(trValues[r,pp] * transfCoef[r,pp]  for pp in primPr) for r in reg, sp in secPr]
a = deepcopy(production)
@meq production[r in reg, sp in secPr]   = sum(trValues[r,pp] * transfCoef[r,pp]  for pp in primPr)
b = deepcopy(production)
[production2[r,sp]   =  sum(trValues2[r,pp] * transfCoef2[r,pp]  for pp in primPrIdx) for r in eachindex(reg), sp in secPrIdx]
c = deepcopy(production2)
@meq production2[r in eachindex(reg), sp in secPrIdx]   =  sum(trValues2[r,pp] * transfCoef2[r,pp]  for pp in primPrIdx)
d = deepcopy(production2)

@test a["us","juice"] ≈ b["us","juice"] ≈ c[findfirst(isequal("us"),reg),findfirst(isequal("juice"),products)] ≈ d[findfirst(isequal("us"),reg),findfirst(isequal("juice"),products)] ≈ 4.7

consumption  = defVars(["reg","prod"],[String,String])
consumption2 = defVars([2,3])

@meq consumption[r in reg, pp in primPr] = production[r,pp] - trValues[r,pp]
@meq consumption[r in reg, sp in secPr]  = production[r, sp]
totalConsumption = sum(consumption[r,p] for r in reg, p in products)
@test totalConsumption == 26.6

@meq consumption2[r in eachindex(reg), pp in primPrIdx] = production2[r,pp] - trValues2[r,pp]
@meq consumption2[r in eachindex(reg), sp in secPrIdx]  = production2[r, sp]
totalConsumption2 = sum(consumption2[r,p] for r in eachindex(reg), p in eachindex(products))
@test totalConsumption2 == 26.6

# *** NEW TEST
# checking getSafe()

@test getSafe(production,("us","banana"))      == 10
@test getSafe(production,("us","oranges"),0.0) == 0.0

# *** NEW TEST
# checking toDataFrame()

content  = [["banana","banana","apple","apple","orange"],["us",missing,"us","eu","us"],[1.1,2.2,3.3,missing,5.5]]
dimNames = ["item","region"]
t        = NDSparse(content...,names=Symbol.(dimNames))
df       = toDataFrame(t)
@test df[((df.item .== "banana") .* ismissing.(df.region)),"value"][1] == 2.2

max_forks_repo_name: sylvaticus/MultiDimEquations


text:
module PIPEs

using ExprRules, StatsBase, LinearAlgebra
using ..PPTs

using ExprOptimization: ExprOptAlgorithm, ExprOptResult
import ExprOptimization: optimize

export PPT, PIPE

"""
    PIPE

Probabilistic Incremental Program Evolution. Example parameters from paper are indicated in parentheses)
# Arguments:
- `ppt_params::PPT`: parameters for PPT  (e.g., [0.8, 0.2])
- `pop_size::Int`: population size
- `iterations::Int`: number of iterations
- `p_elitist::Float64`: elitist update probability (e.g., 0.2)
- `c::Float64`: learning rate multiplier (e.g., 0.1)
- `α::Float64`: learning rate (e.g., 0.05)
- `ϵ::Float64`: fitness constant (e.g., 1)
- `p_mutation::Float64`: mutation probability (e.g., 0.2)
- `β::Float64`: mutation rate (e.g., 0.6)
- `p_threshold::Float64`: prune threshold (e.g., 0.999)
- `max_depth::Int`: maximum depth of derivation tree
"""
struct PIPE <: ExprOptAlgorithm
    ppt_params::PPT
    pop_size::Int
    iterations::Int
    p_elitist::Float64
    c::Float64
    α::Float64
    ϵ::Float64
    p_mutation::Float64
    β::Float64
    p_threshold::Float64
    max_depth::Int
end

"""
    optimize(p::PIPE, grammar::Grammar, typ::Symbol, loss::Function; kwargs...)

Expression tree optimization using the PIPE algorithm with parameters p, grammar 'grammar', start symbol typ, and loss function 'loss'.  Loss function has the form: los::Float64=loss(node::RuleNode, grammar::Grammar).
"""
function optimize(p::PIPE, grammar::Grammar, typ::Symbol, loss::Function; kwargs...)
    pipe(p, grammar, typ, loss; kwargs...)
end

"""
    pipe(p::PIPE, grammar::Grammar, typ::Symbol, loss::Function)

Probabilistic Incremental Program Evolution (PIPE) optimization algorithm with parameters p, grammar 'grammar', start symbol typ, and loss function 'loss'.  Loss function has the form: los::Float64=loss(node::RuleNode, grammar::Grammar).

Reference: R. Salustowicz and J. Schmidhuber, "Probabilistic Incremental Program Evolution",
    Evolutionary Computation, vol. 5, no. 2, pp. 123-141, 1997.
"""
function pipe(p::PIPE, grammar::Grammar, typ::Symbol, loss::Function)
    iseval(grammar) && error("PIPE does not support _() functions in the grammar")

    best_tree, best_loss = RuleNode(0), Inf
    pp = p.ppt_params
    pop = Vector{RuleNode}(undef,p.pop_size)
    losses = Vector{Float64}(undef,p.pop_size)

    ppt = PPTNode(p.ppt_params, grammar)
    for i = 1:p.iterations
        if rand() < p.p_elitist #elitist learning
            update!(p, ppt, grammar, best_tree, best_loss, best_loss)
        else #generational learning
            for j = 1:p.pop_size
                pop[j] = rand(pp, ppt, grammar, typ)
            end
            best_tree, best_loss = evaluate!(loss, grammar, pop, losses, best_tree, best_loss)
            update!(p, ppt, grammar, pop[1], losses[1], best_loss)
            mutate!(ppt, grammar, pop[1], p.p_mutation, p.β)
            prune!(ppt, grammar, p.p_threshold)
        end
    end
    ExprOptResult(best_tree, best_loss, get_executable(best_tree, grammar), nothing)
end

"""
    p_target(pp::PPT, ppt::PPTNode, grammar::Grammar, x_best::RuleNode, y_best::Float64,
                  y_elite::Float64, α::Float64, ϵ::Float64)

Compute the target probability p_target.  See PIPE paper for description of equation.
"""
function p_target(pp::PPT, ppt::PPTNode, grammar::Grammar, x_best::RuleNode, y_best::Float64,
                  y_elite::Float64, α::Float64, ϵ::Float64)
    p_best = probability(pp, ppt, grammar, x_best)
    return p_best + (1-p_best)*α*(ϵ - y_elite)/(ϵ - y_best)
end

"""
    update!(p::PIPE, ppt::PPTNode, grammar::Grammar, x::RuleNode, y::Float64,
                 y_elite::Float64)

Update the ppt probabilities toward individual x with loss y given elite loss y_elite.
"""
function update!(p::PIPE, ppt::PPTNode, grammar::Grammar, x::RuleNode, y::Float64,
                 y_elite::Float64)
    pp, c, α, ϵ = p.ppt_params, p.c, p.α, p.ϵ
    p_targ = p_target(pp, ppt, grammar, x, y, y_elite, α, ϵ)
    while probability(pp, ppt, grammar, x) < p_targ
        _update!(ppt, grammar, x, c, α)
    end
    return ppt
end
function _update!(ppt::PPTNode, grammar::Grammar, x::RuleNode, c::Float64, α::Float64)
    typ = return_type(grammar, x)
    i = something(findfirst(isequal(x.ind),grammar[typ]), 0)
    p = ppt.ps[typ]
    p[i] += c*α*(1-p[i])
    psum = sum(p)
    for j in 1 : length(p)
        if j != i
            p[j] *= (1- (1-psum)/(p[j]-psum))
        end
    end
    for (pptchild,xchild) in zip(ppt.children, x.children)
        _update!(pptchild, grammar, xchild, c, α)
    end
    return ppt
end

"""
    mutate!(ppt::PPTNode, grammar::Grammar, x_best::RuleNode, p_mutation::Float64, β::Float64;
        sqrtlen::Float64=sqrt(length(x_best)))

Mutate ppt node towards individual x_best using mutation probability p_mutation, mutation rate β using the
square-root length criteria.
"""
function mutate!(ppt::PPTNode, grammar::Grammar, x_best::RuleNode, p_mutation::Float64, β::Float64;
    sqrtlen::Float64=sqrt(length(x_best)))

    typ = return_type(grammar, x_best)
    p = ppt.ps[typ]
    prob = p_mutation/(length(p)*sqrtlen)
    for i in 1 : length(p)
        if rand() < prob
            p[i] += β*(1-p[i])
        end
    end
    normalize!(p, 1)
    for (pptchild,xchild) in zip(ppt.children, x_best.children)
        mutate!(pptchild, grammar, xchild, p_mutation, β,
                sqrtlen=sqrtlen)
    end
    ppt
end

"""
    evaluate!(loss::Function, grammar::Grammar, pop::Vector{RuleNode}, losses::Vector{Float64},
        best_tree::RuleNode, best_loss::Float64)

Evaluate the loss function for population and sort.  Update the globally best tree.
"""
function evaluate!(loss::Function, grammar::Grammar, pop::Vector{RuleNode}, losses::Vector{Float64},
                   best_tree::RuleNode, best_loss::Float64)
    Threads.@threads for i in 1:length(losses)
        losses[i] = loss(pop[i], grammar)
    end
    perm = sortperm(losses)
    pop[:], losses[:] = pop[perm], losses[perm]
    if losses[1] < best_loss
        best_tree, best_loss = pop[1], losses[1]
    end
    (best_tree, best_loss)
end

end #module


text:
module FlagColors
using Colors, ColorSchemes
export flagcolors

include("../resources/flagcolors.jl")

end


text:
using MINLPTests, JuMP, SCIP, Test

const OPTIMIZER = JuMP.with_optimizer(SCIP.Optimizer, display_verblevel=0)
const OBJTOL = 1e-4
const PRIMALTOL = 1e-3
const DUALTOL = NaN  # to disable the query

MINLPTests.test_directory(
    "nlp", OPTIMIZER,
    objective_tol=OBJTOL, primal_tol=PRIMALTOL, dual_tol=DUALTOL,
    termination_target=MINLPTests.TERMINATION_TARGET_GLOBAL,
    primal_target=MINLPTests.PRIMAL_TARGET_GLOBAL,
    include=[
        "005_010",
        "007_010",
    ])

MINLPTests.test_directory(
    "nlp-cvx", OPTIMIZER,
    objective_tol=OBJTOL, primal_tol=PRIMALTOL, dual_tol=DUALTOL,
    termination_target=MINLPTests.TERMINATION_TARGET_GLOBAL,
    primal_target=MINLPTests.PRIMAL_TARGET_GLOBAL,
    include=[
        "001_010",
        "002_010",
        "101_010",
        "101_012",
        "102_010",
        "102_011",
        "102_012",
        "103_010",
        "103_011",
        "103_012",
        "103_013",
        "103_014",
        "104_010",
        "105_010",
        "105_011",
        "105_012",
        "105_013",
        "201_010",
        "201_011",
        "202_010",
        "202_011",
        "202_012",
        "202_013",
        "202_014",
        "203_010",
        "204_010",
        "205_010",
    ])

MINLPTests.test_directory(
    "nlp-mi", OPTIMIZER,
    objective_tol=OBJTOL, primal_tol=PRIMALTOL, dual_tol=DUALTOL,
    termination_target=MINLPTests.TERMINATION_TARGET_GLOBAL,
    primal_target=MINLPTests.PRIMAL_TARGET_GLOBAL,
    include=[
        "005_010",
        "007_010",
        "007_020",
])


text:
include("../common.jl")
build_julia(ARGS, v"1.5.4")

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
# --------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120--------130
# vector potential
function getAt(t::TF, pulse::Pulse{At, TF}) where TF
	@unpack tL, tR, T, A, ω, ϕ = pulse
	if t < tL || tR < t
		return 0
	end
	return A*get_env(t, pulse)*sin(ω*t + ϕ)
end
function getAt(t::TVF, pulse::Pulse{At}) where {TVF <: AbstractVector{<:AbstractFloat}}
	At = zeros(eltype(t), size(t))
	for it in 1:length(t)
		At[it] = getAt(t[it], pulse)
	end
	return At
end
function getAt(t::TF, laser::Vector{Pulse{At, TF}}) where {TF <: AbstractFloat}
	At = zero(t)
	for ip in 1:length(laser)
		At += getAt(t, laser[ip])
	end
	return At
end
function getAt(t::TVF, laser::Vector{Pulse{At}}) where {TVF <: AbstractVector{<:AbstractFloat}}
	At = zeros(eltype(t), size(t))
	for it in 1:length(t)
		At[it] = getAt(t[it], laser)
	end
	return At
end
# electric field
# Et = -dA/dt
function getEt(t::TF, pulse::Pulse{At})where TF <: AbstractFloat
	@unpack tL, tR, T = pulse
	if t < tL || tR < t
		return 0
	end
	dt = T*1e-5
	Et =  (1/60dt)*(getAt(t + 3dt, pulse) - getAt(t - 3dt, pulse)) -
			(3/20dt)*(getAt(t + 2dt, pulse) - getAt(t - 2dt, pulse)) +
			(3/4dt)*(getAt(t + dt, pulse) - getAt(t - dt, pulse))
	return -Et
end
function getEt(t::TVF, pulse::Pulse{At}) where {TVF <: AbstractVector{<:AbstractFloat}}
	Et = zeros(eltype(t), size(t))
	for it in 1:length(t)
		Et[it] = getEt(t[it], pulse)
	end
	return Et
end
function getEt(t::TF, laser::Vector{Pulse{At}}) where {TF <: AbstractFloat}
	Et = zero(t)
	for ip in 1:length(laser)
		Et += getEt(t, laser[ip])
	end
	return Et
end
function getEt(t::TVF, laser::Vector{Pulse{At}}) where {TVF <: AbstractVector{<:AbstractFloat}}
	Et = zeros(eltype(t), size(t))
	for it in 1:length(t)
		Et[it] = getEt(t[it], laser)
	end
	return Et
end


text:
module BARON

using MathProgBase
using MathProgBase.SolverInterface
using Compat

export BaronSolver
struct BaronSolver <: AbstractMathProgSolver
    options
end

BaronSolver(;kwargs...) = BaronSolver(kwargs)

const baron_exec = ENV["BARON_EXEC"]

mutable struct BaronMathProgModel <: AbstractNonlinearModel
    options

    xˡ::Vector{Float64}
    xᵘ::Vector{Float64}
    gˡ::Vector{Float64}
    gᵘ::Vector{Float64}

    nvar::Int
    ncon::Int

    obj
    constrs

    vartypes::Vector{Symbol}

    v_names::Vector{String}
    c_names::Vector{String}

    sense::Symbol

    x₀::Vector{Float64}

    probfile::String
    sumfile::String
    resfile::String

    objval::Float64
    dualbound::Float64
    solution::Vector{Float64}
    walltime::Float64
    status::Symbol

    d::AbstractNLPEvaluator

    function BaronMathProgModel(;options...)
        dir = tempdir()
        push!(options, (:ResName, joinpath(dir, "res.lst")))
        push!(options, (:TimName, joinpath(dir, "tim.lst")))
        push!(options, (:SumName, joinpath(dir, "sum.lst")))
        new(options,
        zeros(0),
        zeros(0),
        zeros(0),
        zeros(0),
        0,
        0,
        :(0),
        Expr[],
        Symbol[],
        String[],
        String[],
        :Min,
        zeros(0),
        "",
         "",
        "",
        NaN,
        NaN,
        [NaN],
        0.0,
        :NotSolved)
    end
end

MathProgBase.NonlinearModel(s::BaronSolver) = BaronMathProgModel(;s.options...)
MathProgBase.LinearQuadraticModel(s::BaronSolver) = MathProgBase.NonlinearToLPQPBridge(MathProgBase.NonlinearModel(s))

verify_support(c) = c

function verify_support(c::Expr)
    if c.head == :comparison
        map(verify_support, c.args)
        return c
    end
    if c.head == :call
        if c.args[1] in (:+, :-, :*, :/, :exp, :log)
            return c
        elseif c.args[1] in (:<=, :>=, :(==))
            map(verify_support, c.args[2:end])
            return c
        elseif c.args[1] == :^
            @assert isa(c.args[2], Real) || isa(c.args[3], Real)
            return c
        else # TODO: do automatic transformation for x^y, |x|
            error("Unsupported expression $c")
        end
    end
    return c
end

function MathProgBase.loadproblem!(m::BaronMathProgModel,
    nvar, ncon, xˡ, xᵘ, gˡ, gᵘ, sense,
    d::MathProgBase.AbstractNLPEvaluator)

    @assert nvar == length(xˡ) == length(xᵘ)
    @assert ncon == length(gˡ) == length(gᵘ)

    m.xˡ, m.xᵘ = xˡ, xᵘ
    m.gˡ, m.gᵘ = gˡ, gᵘ
    m.sense = sense
    m.nvar, m.ncon = nvar, ncon
    m.d = d

    m.v_names = ["x$i" for i in 1:nvar]
    m.c_names = ["e$i" for i in 1:ncon]

    MathProgBase.initialize(d, [:ExprGraph])

    m.obj = verify_support(MathProgBase.obj_expr(d))
    m.vartypes = fill(:Cont, nvar)

    m.constrs = map(1:ncon) do c
        verify_support(MathProgBase.constr_expr(d,c))
    end

    dir = tempdir()
    m.probfile = joinpath(dir, "baron_problem.bar")
    m.sumfile  = joinpath(dir, "sum.lst")
    m.resfile  = joinpath(dir, "res.lst")
    m
end

function MathProgBase.setvartype!(m::BaronMathProgModel, cat::Vector{Symbol})
    @assert all(x-> (x in (:Cont,:Bin,:Int)), cat)
    m.vartypes = cat
end

function print_var_definitions(m, fp, header, condition)
    idx = filter(condition, 1:m.nvar)
    if !isempty(idx)
        println(fp, header, join([m.v_names[i] for i in idx], ", "), ";")
    end
end

# to_str(x::Int) = string(x)
# to_str(x) = (@show(x); string(float(x)))
to_str(x) = string(x)

function to_str(c::Expr)
    if c.head == :comparison
        if length(c.args) == 3
            return join([to_str(c.args[1]), c.args[2], c.args[3]], " ")
        elseif length(c.args) == 5
            return join([c.args[1], c.args[2], to_str(c.args[3]),
                         c.args[4], c.args[5]], " ")
        end
    elseif c.head == :call
        if c.args[1] in (:<=,:>=,:(==))
            if length(c.args) == 3
                return join([to_str(c.args[2]), to_str(c.args[1]), to_str(c.args[3])], " ")
            elseif length(c.args) == 5
                return join([to_str(c.args[1]), to_str(c.args[2]), to_str(c.args[3]), to_str(c.args[4]), to_str(c.args[5])], " ")
            end
        elseif c.args[1] in (:+,:-,:*,:/,:^)
            if all(d->isa(d, Real), c.args[2:end]) # handle unary case
                return string(eval(c))
            elseif c.args[1] == :- && length(c.args) == 2
                return string("(-$(to_str(c.args[2])))")
            else
                return string("(", join([to_str(d) for d in c.args[2:end]], string(c.args[1])), ")")
            end
        elseif c.args[1] in (:exp,:log)
            if isa(c.args[2], Real)
                return string(eval(c))
            else
                return string(c.args[1], "( ", to_str(c.args[2]), " )")
            end
        end
    elseif c.head == :ref
        if c.args[1] == :x
            @assert isa(c.args[2], Int)
            return "x$(c.args[2])"
        else
            error("Unrecognized reference expression $c")
        end
    end
end

function write_bar_file(m::BaronMathProgModel)
    fp = open(m.probfile, "w")

    # First: process any options
    println(fp, "OPTIONS{")
    for (opt,setting) in m.options
        if isa(setting, AbstractString) # wrap it in quotes
             println(fp, "$opt: $('"')$setting$('"');")
        else
            println(fp, "$opt: $setting;")
        end
    end
    println(fp, "}")
    println(fp)

    # Next, define variables
    print_var_definitions(m, fp, "BINARY_VARIABLES ",   v->(m.vartypes[v]==:Bin))
    print_var_definitions(m, fp, "INTEGER_VARIABLES ",  v->(m.vartypes[v]==:Int))
    print_var_definitions(m, fp, "POSITIVE_VARIABLES ", v->(m.vartypes[v]==:Cont && m.xˡ[v]==0))
    print_var_definitions(m, fp, "VARIABLE ",           v->(m.vartypes[v]==:Cont && m.xˡ[v]!=0))
    println(fp)

    # Print variable bounds
    if any(c->!isinf(c), m.xˡ)
        println(fp, "LOWER_BOUNDS{")
        for (i,l) in enumerate(m.xˡ)
            if !isinf(l)
                println(fp, "$(m.v_names[i]): $l;")
            end
        end
        println(fp, "}")
        println(fp)
    end
    if any(c->!isinf(c), m.xᵘ)
        println(fp, "UPPER_BOUNDS{")
        for (i,u) in enumerate(m.xᵘ)
            if !isinf(u)
                println(fp, "$(m.v_names[i]): $u;")
            end
        end
        println(fp, "}")
        println(fp)
    end

    # Now let's declare the equations
    if !isempty(m.constrs)
        println(fp, "EQUATIONS ", join(m.c_names, ", "), ";")
        for (i,c) in enumerate(m.constrs)
            str = to_str(c)
            println(fp, "$(m.c_names[i]): $str;")
        end
        println(fp)
    end

    # Now let's do the objective
    print(fp, "OBJ: ")
    print(fp, m.sense == :Min ? "minimize " : "maximize ")
    print(fp, to_str(m.obj))
    println(fp, ";")
    println(fp)

    if !isempty(m.x₀)
        println(fp, "STARTING_POINT{")
        for (i,v) in enumerate(m.x₀)
            println(fp, "$(m.v_names[i]): $v;")
        end
        println(fp, "}")
    end
    close(fp)
end

const normal_status = [
    ["***", "Normal", "completion", "***"]
]

const userlimits_status = [
    ["***", "Max.", "allowable", "nodes", "in", "memory", "reached", "***"],
    ["***", "Max.", "allowable", "BaR", "iterations", "reached", "***"],
    ["***", "Max.", "allowable", "CPU", "time", "exceeded", "***"],
    ["***", "Max.", "allowable", "time", "exceeded", "***"],
    ["***", "Problem", "is", "numerically", "sensitive", "***"],
    ["***", "Heuristic", "termination", "***"],
    ["***", "Globality", "is", "therefore", "not", "guaranteed", "***"]
]

const infeasible_status = [
    ["***", "Problem", "is", "infeasible", "***"],
    ["***", "No", "feasible", "solution", "was", "found", "***"],
    ["***", "Infeasibility", "is", "therefore", "not", "guaranteed", "***"],
    ["***", "Insufficient", "Memory", "for", "Data", "structures", "***"],
]

const unbounded_status = [
    ["***", "Problem", "is", "unbounded", "***"],
    ["***", "User" ,"did", "not", "provide", "appropriate", "variable", "bounds", "***"]
]

const numerical_status = [
    ["***", "Problem", "is", "numerically", "sensitive", "***"]
]

const interrupt_status = [
    ["***", "Search", "interrupted", "by", "user", "***"]
]

const error_status = [
    ["***", "A", "potentially", "catastrophic", "access", "violation", "just", "took", "place", "***"]
]

function read_results(m::BaronMathProgModel)

    # First, we read the summary file to get the solution status
    stat_code = []
    fp = open(m.sumfile, "r")
    stat = :Undefined
    t = -1.0
    n = -99
    while true
        line = readline(fp)
        spl = split(chomp(line))
        if !isempty(spl) && spl[1] == "***"
            push!(stat_code, spl)
        elseif length(spl)>=3 && spl[1:3] == ["Wall", "clock", "time:"]
            t = parse(Float64,match(r"\d+.\d+", line).match)
        elseif length(spl)>=3 && spl[1:3] == ["Best", "solution", "found"]
            n = parse(Int,match(r"-?\d+", line).match)
        # Grab dual bound if solved during presolve (printed directly to summary file)
        elseif length(spl)>=3 && spl[1:3] == ["Lower", "bound", "is"]
            m.dualbound = parse(Float64, spl[4])
        # Grab dual bound if branching (need to get it from solver update log)
        elseif spl == ["Iteration", "Open", "nodes", "Time", "(s)", "Lower", "bound", "Upper", "bound"]
            while true
                line = readline(fp)
                spl = split(chomp(line))
                if isempty(spl)
                    break
                end
                # Lowerbound is 4th column in table, but log line might include * for heuristic solution
                m.dualbound = parse(Float64, spl[end-1])
            end
        end
        eof(fp) && break
    end
    close(fp)

    t < 0.0 && warn("No solution time is found in sum.lst")
    n == -99 && error("No solution node information found sum.lst")

    # Now navigate to the right problem status by looking at the main status
    if stat_code[1] in userlimits_status
        m.status = :UserLimit
        n == -3 && (m.status = :Infeasible)
    elseif stat_code[1] in normal_status
        m.status = :Optimal
        n == -3 && (m.status = :Infeasible)
    elseif stat_code[1] in infeasible_status
        m.status = :Infeasible
    elseif stat_code[1] in unbounded_status
        m.status = :Unbounded
    elseif stat_code[1] in numerical_status
        m.status = :NumericalSensitive
        n == -3 && (m.status = :Infeasible)
    elseif stat_code[1] in interrupt_status
        m.status = :Interrupted
        n == -3 && (m.status = :Infeasible)
    elseif stat_code[1] in error_status
        m.status = :Error
    end

    m.walltime = t # Track the time

    # Next, we read the results file to get the solution
    x = fill(NaN, m.nvar)
    m.objval = NaN
    if n != -3 # parse as long as there exist a solution
        fp = open(m.resfile, "r")
        while true
            startswith(readline(fp), "The best solution found") && break
            eof(fp) && error()
        end
        readline(fp)
        readline(fp)

        while true
            line = chomp(readline(fp))
            parts = split(line)
            isempty(parts) && break
            mt = match(r"\d+", parts[1])
            mt == nothing && error("Cannot find appropriate variable index from $(parts[1])")
            v_idx = parse(Int,mt.match)
            v_val = float(parts[3])
            x[v_idx] = v_val
        end
        m.solution = x
        line = readline(fp)
        val = match(r"[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?", chomp(line))
        if val != nothing
            m.objval = float(val.match)
        end
    end
    nothing
end

MathProgBase.setwarmstart!(m::BaronMathProgModel, v::Vector{Float64}) = (m.x₀ = v)

function MathProgBase.optimize!(m::BaronMathProgModel)
    write_bar_file(m)
    run(`$baron_exec $(m.probfile)`)

    read_results(m)
end

MathProgBase.status(m::BaronMathProgModel) = m.status

MathProgBase.numvar(m::BaronMathProgModel) = m.nvar
MathProgBase.getsolution(m::BaronMathProgModel) = m.solution
MathProgBase.getobjval(m::BaronMathProgModel) = m.objval
MathProgBase.getobjbound(m::BaronMathProgModel) = m.dualbound
MathProgBase.getsolvetime(m::BaronMathProgModel) = m.walltime

end


text:
abstract type AbstractSymmetricConnectivityModel end

@inline (mdl::AbstractSymmetricConnectivityModel)(n_vertices, n_pins, n_local_nets, n_remote_nets, k) = mdl(n_vertices, n_pins, n_local_nets, n_remote_nets)

struct AffineSymmetricConnectivityModel{Tv} <: AbstractSymmetricConnectivityModel
    α::Tv
    β_vertex::Tv
    β_pin::Tv
    β_local_net::Tv
    β_remote_net::Tv
end

function AffineSymmetricConnectivityModel(; α = false, β_vertex = false, β_pin = false, β_local_net = false, β_remote_net = false)
    AffineSymmetricConnectivityModel(promote(α, β_vertex, β_pin, β_local_net, β_remote_net)...)
end

@inline cost_type(::Type{AffineSymmetricConnectivityModel{Tv}}) where {Tv} = Tv

(mdl::AffineSymmetricConnectivityModel)(n_vertices, n_pins, n_local_nets, n_remote_nets) = mdl.α + n_vertices * mdl.β_vertex + n_pins * mdl.β_pin + n_local_nets * mdl.β_local_net + n_remote_nets * mdl.β_remote_net

struct SymmetricConnectivityOracle{Ti, Net, DiaNet, Mdl} <: AbstractOracleCost{Mdl}
    pos::Vector{Ti}
    net::Net
    dianet::DiaNet
    mdl::Mdl
end

oracle_model(ocl::SymmetricConnectivityOracle) = ocl.mdl

function oracle_stripe(hint::AbstractHint, mdl::AbstractSymmetricConnectivityModel, A::SparseMatrixCSC{Tv, Ti}; net=nothing, adj_A=nothing, kwargs...) where {Tv, Ti}
    @inbounds begin
        m, n = size(A)
        @assert m == n
        N = nnz(A)
        pos = A.colptr

        if net === nothing
            net = netcount(hint, A)
        end

        dianet = dianetcount(hint, A)

        return SymmetricConnectivityOracle(pos, net, dianet, mdl)
    end
end

@inline function (cst::SymmetricConnectivityOracle{Ti, Mdl})(j::Ti, j′::Ti, k...) where {Ti, Mdl}
    @inbounds begin
        w = cst.pos[j′] - cst.pos[j]
        d = cst.net(j, j′)
        r = cst.dianet(j, j′) - (j′ - j)
        l = d - r
        return cst.mdl(j′ - j, w, l, r, k...)
    end
end

@inline function (stp::Step{Ocl})(_j, _j′, _k...) where {Ti, Mdl, Ocl <: SymmetricConnectivityOracle{Ti, Mdl}}
    @inbounds begin
        cst = stp.ocl
        j = destep(_j)
        j′ = destep(_j′)
        k = maptuple(destep, _k...)
        w = cst.pos[j′] - cst.pos[j]
        d = Step(cst.net)(_j, _j′)
        r = Step(cst.dianet)(_j, _j′) - (j′ - j)
        l = d - r
        return cst.mdl(j′ - j, w, l, r, k...)
    end
end


text:
module lmp_data

using Statistics
using Printf

include("Constants.jl")
#=
export Const_k_b, Const_n_a, Const_density_wat,
    Const_kcal2j, Const_kcalm2j, Const_kcalm2t,
    Const_g2kg, Const_kg2g, Const_gm2g, Const_gm2kg,
    Const_an2m, Const_an2nm, Const_nm2m, Const_cm2an, Const_dm2m, Const_m2dm,
    Const_fs2s, Const_ps2s, Const_ns2s
=#

include("useful_funcs.jl")
export dist, str2float, diag, copy_array, leastsq

include("read_dump.jl")
export read_dump, read_dump_prim, split_info!, split_atom, momentum!, mass!, time_step!, coord_scl!

include("box.jl")
export genr_box_diag, genr_box_inv, PBC!

include("diffusion.jl")
export MSD!, diffusion

include("RDF.jl")
export RDF, RDF_double, num_coordination

end # module


text:
abstract GrammarToken

type GrammarLiteral <: GrammarToken
  value::AbstractString
end
type GrammarSymbol <: GrammarToken
  value::AbstractString
end

type Grammar
  root::AbstractString
  rules::Dict{AbstractString, Vector{Vector{GrammarToken}}}

  function Grammar(root::AbstractString, rules::Dict{AbstractString, Vector{AbstractString}})
    p_rules = Dict{AbstractString, Vector{Vector{GrammarToken}}}()
    for (sym, entries) in rules
      p_rules[sym] = [parse_rule_entry(entry) for entry in  entries]
    end
    new(root, p_rules)
  end
end

function derive(g::Grammar, s::Vector{Int}, max_wraps::Int)
  derivation = ""

  q = GrammarToken[GrammarSymbol(g.root)]
  i = 1
  len = length(s)
  wraps = 0

  while !isempty(q)
    token = shift!(q)

    if isa(token, GrammarLiteral)
      derivation = derivation * token.value
    else
      options = g.rules[token.value]
      if length(options) == 1
        prepend!(q, options[1])
      else
        if i >= len
          if wraps == max_wraps
            error("Genotype exhausted")
          end
          i = 1
          wraps += 1
        end
        prepend!(q, options[(s[i] % length(options)) + 1])
        i += 1
      end
    end
  end
  return derivation
end

function parse_rule_entry(entry::AbstractString)
  tokens = GrammarToken[]
  while !isempty(entry)
    left, tag, entry = partition(entry, r"<[^\<\>]+>")

    if !isempty(left)
      push!(tokens, GrammarLiteral(left))
    end
    if !isempty(tag)
      push!(tokens, GrammarSymbol(tag[2:end-1]))
    end
  end
  return tokens
end


text:
using StanSample, Test

binom_model = "
// Inferring a Rate
data {
  int<lower=1> n;
  int<lower=0> k;
}
parameters {
  real<lower=0,upper=1> theta;
  real<lower=0,upper=1> thetaprior;
}
model {
  // Prior Distribution for Rate Theta
  theta ~ beta(1, 1);
  thetaprior ~ beta(1, 1);

  // Observed Counts
  k ~ binomial(n, theta);
}
generated quantities {
  int<lower=0> postpredk;
  int<lower=0> priorpredk;

  postpredk <- binomial_rng(n, theta);
  priorpredk <- binomial_rng(n, thetaprior);
}
"

sm = SampleModel("binomial", binom_model)

binom_data = Dict("n" => 10, "k" => 5)

rc = stan_sample(sm, data=binom_data)

if success(rc)
  # Fetch cmdstan summary data frame.
  df = read_summary(sm)
  @test df[df.parameters .== :theta, :mean][1] ≈ 0.5 rtol=0.1

end


text:
using EmDee
using Test
using CUDA
import Chemfiles

# function test_cells()
#     N = 1000
#     L = 1.0
#     cutoff = 0.2
#     x = CUDA.rand(Float32, 3, N)
#     y = x .+ 0.01
#     cells_x = Cells(x, L, cutoff)
#     cells_y = Cells(y, L, cutoff)
#     update_cells!(cells_x, y, L)
#     return all(Array(cells_x.index) .== Array(cells_y.index)) &&
#            all(Array(cells_x.population) .== Array(cells_y.population))
# end

function test_compute_nonbonded(xyz_file, L, cutoff, switch)
    frame = read(Chemfiles.Trajectory(xyz_file))
    xyz_data = Chemfiles.positions(frame)
    positions = CUDA.cu(xyz_data)
    N = size(xyz_data, 2)
    model = LennardJonesModel(cutoff, switch)
    atoms = CUDA.cu(fill(LennardJonesAtom(1, 1), N))

    forces_ref = CUDA.zeros(Float32, 3, N)
    energies_ref = CUDA.zeros(N)
    virials_ref = CUDA.zeros(N)
    naively_compute_nonbonded!(forces_ref, energies_ref, virials_ref, positions, L, model, atoms)

    tiles = nonbonded_computation_tiles(N)
    forces = CUDA.zeros(Float32, 3, N)
    energies = CUDA.zeros(N)
    virials = CUDA.zeros(N)
    compute_nonbonded!(forces, energies, virials, positions, L,
                       tiles, model, atoms, Val(FORCES | ENERGIES | VIRIALS))

    return maximum(forces .- forces_ref) < 1.0f-4 &&
           maximum(energies .- energies_ref) < 1.0f-4 &&
           maximum(virials .- virials_ref) < 1.0f-4
end

function test_system(pdb_file, xml_file)
    force_field = ForceField(xml_file)
    system = System(pdb_file, force_field)
    topology = Chemfiles.Topology(system)
    return length(system) == 1519 && Chemfiles.count_residues(topology) == 500
end

@testset "EmDee.jl[CPU]" begin
    @test test_system(joinpath(@__DIR__, "data", "dibenzo-p-dioxin-in-water.").*["pdb", "xml"]...)
end

if CUDA.functional()
    @testset "EmDee.jl[GPU]" begin
        # @test test_cells()
        @test test_compute_nonbonded(joinpath(@__DIR__, "data", "lj_sample.xyz"), 10, 3, 2.5)
    end
end


text:
@testset "html partition" begin
    parts = """
        ABC1 <!-- DEF --> ABC2 <script> GHI </script> ABC3 {{ JKL }} ABC4
        """ |> FP.html_partition
    for (i, c) in enumerate(("ABC1", "DEF", "ABC2", "> GHI", "ABC3", "JKL", "ABC4"))
        @test isapproxstr(FP.content(parts[i]), c)
    end
end


text:
#****************************************************************#
# VeysiADN 12 May 2021
# Advanced Robotics Homework 5 Codes
#****************************************************************#
using LinearAlgebra
const lAlgebra = LinearAlgebra
#****************************************************************#
#Computes exponential coordinate ‘so3mat’
#(skew symmetric matrix) of the rotation matrix R

function MatrixLog3(R)
    acosinput = (lAlgebra.tr(R) - 1) / 2
    if acosinput >= 1
        return zeros(3, 3)
    else
        theta = acos(acosinput)
        return 1 / 2 / sin(theta) * (R - R')
    end
end
#****************************************************************#
# Computes exponential coordinate ‘se3mat’
# (twist matrix) of the rigid body transformation matrix T
function MatrixLog6(T)
    if( size(T)!=(4,4) )
        println("!Warning! : Transformation matrix should be 4x4.")
        return
    end
    R = T[1:3,1:3]
    t = T[1:3,4]
    omg_h = MatrixLog3(R)
    w = [omg_h[3,2];omg_h[1,3];omg_h[2,1]]
    if omg_h == zeros(3, 3)
        return vcat(hcat(zeros(3, 3), T[1:3, 4]), [0 0 0 0])
    else
        theta = acos((lAlgebra.tr(R) - 1) / 2)
        return vcat(hcat(omg_h,inv((I-R)*omg_h + w*w'*theta)*t),
                          [0 0 0 0])
    end
end
#****************************************************************#

max_stars_repo_name: veysiadn/AdvancedRoboticsJulia

max_issues_repo_name: veysiadn/AdvancedRoboticsJulia

max_forks_repo_name: veysiadn/AdvancedRoboticsJulia


text:
const TransparentColorN{N,C<:Color,T} = TransparentColor{C,T,N}

struct Hue{T<:Real}
    angle::T
    isgray::Bool
    Hue(angle::T, isgray::Bool=false) where {T <: Real} = new{T}(angle, isgray)
    Hue(c::C) where {T, C<:Union{HSV{T}, HSL{T}, HSI{T},
                                 AHSV{T}, AHSL{T}, AHSI{T}, HSVA{T}, HSLA{T}, HSIA{T}}} =
        new{T}(c.h, c.s == zero(T))
    Hue(c::C) where {T, C<:Union{LCHab{T}, LCHuv{T},
                                 ALCHab{T}, ALCHuv{T}, LCHabA{T}, LCHuvA{T}}} =
        new{T}(c.h, c.c == zero(T))
end

const SeparableBlendMode = Union{
    BlendMode{:normal},
    BlendMode{:multiply},
    BlendMode{:screen},
    BlendMode{:overlay},
    BlendMode{:darken},
    BlendMode{:lighten},
    BlendMode{Symbol("color-dodge")},
    BlendMode{Symbol("color-burn")},
    BlendMode{Symbol("hard-light")},
    BlendMode{Symbol("soft-light")},
    BlendMode{:difference},
    BlendMode{:exclusion},
}

const NonSeparableBlendMode = Union{
    BlendMode{:hue},
    BlendMode{:saturation},
    BlendMode{:color},
    BlendMode{:luminosity},
}

const OpaqueInSourceOutOperation = Union{
    typeof(CompositeCopy),
    typeof(CompositeSourceOver),
    typeof(CompositeSourceIn),
    typeof(CompositeSourceAtop),
}

const OpaqueInDestinationOutOperation = Union{
    typeof(CompositeDestination),
    typeof(CompositeDestinationOver),
    typeof(CompositeDestinationIn),
    typeof(CompositeDestinationAtop),
}


text:
using DeepRL
using Flux
using Random
using ProgressMeter
# using Plots
using TensorBoardLogger
using Logging
using LinearAlgebra
using BSON: @save
# using Distributions

function main_experiment(seed,
                         num_frames,
                         save_loc,
                         checkin_step;
                         gamename="breakout",
                         prog_meter_offset=0)

    max_episode_length = 27000 # From Dopamine.

    save_loc = joinpath(save_loc, "run_$(seed)")
    if !isdir(save_loc)
        mkpath(save_loc)
    end

    model_save_loc = joinpath(save_loc, "models")
    if !isdir(model_save_loc)
        mkdir(model_save_loc)
    end

    res_save_file = joinpath(save_loc, "results.bson")

    Random.seed!(Random.GLOBAL_RNG, seed)

    env = DeepRL.DopamineALEAtari(gamename, rand(UInt16))
    agent = DeepRL.DopamineDQNBaseline(env; minimal_action_set=true)

    @info "Number Actions: $(length(agent.acting_policy.action_set))"

    # env = DeepRL.RevisitingALEAtari(gamename, rand(UInt16))
    # agent = DeepRL.RevisitingALEDQNBaseline(env)

    total_rews = Array{Int,1}()
    steps = Array{Int,1}()

    front = ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇']
    p = ProgressMeter.Progress(
        num_frames;
        dt=1,
        desc="Step: ",
        barglyphs=ProgressMeter.BarGlyphs('|','█',front,' ','|'),
        barlen=Int64(floor(100/length(front))),
        offset=prog_meter_offset)

    start_time = time()

    lg=TBLogger(joinpath(dirname(save_loc), "tensorboard_logs/run_$(seed)"), min_level=Logging.Info)

    eps = 0
    total_steps = 0
    prev_log_steps = 0
    with_logger(lg) do
        while sum(steps) < num_frames
            cur_step = sum(steps)
            episode_cut_off = min(max_episode_length, num_frames - sum(steps))
            tr, stp =
                run_episode!(env, agent, episode_cut_off) do (s, a, s′, r)
                    next!(p, showvalues=[(:step, cur_step), (:episode, eps)])
                    if (cur_step) % checkin_step == 0
                        model = cpu(agent.model)
                        total_time = time() - start_time
                        @save model_save_loc*"/step_$(cur_step).bson" model total_rews steps total_time
                    end
                    cur_step+=1
                end
            push!(total_rews, tr)
            push!(steps, stp)

            @info "Train" returns = tr steps = stp
            eps += 1
        end

    end
    end_time = time()
    close(env)

    model = cpu(agent.model)
    total_time = end_time - start_time
    @save res_save_file model total_rews steps total_time

end


text:
module MaxRectangle
export max_rectangle

using Parameters
using StatsBase
using ProgressMeter
using PyCall

plt = pyimport("matplotlib.pyplot")
patches = pyimport("matplotlib.patches")
cmap = pyimport("matplotlib.cm")
colors = pyimport("matplotlib.colors")

@with_kw mutable struct Cell
    dim::Int64
    rect_list::Array{Vector{Int64}} = []
    dist_to_zero::Vector{Int64} = [0 for __ in 1:dim]
    max_rect::Vector{Int64} = [0 for __ in 1:dim]
    color::NTuple{4, Float64} = (1.0, 1.0, 1.0, 0.0)
    point_list::Array{Vector{Float64}} = []
    point_list_rects::Array{Array{Vector{Float64}}} = []
    point_list_max_rect::Array{Vector{Float64}} = []
end

function Base.show(io::IO, object::Cell)
    println(io, "$(object.dist_to_zero)\n$(object.max_rect)")
end

function find_candidates(neighbors, grid, grid_idx, cellsize)
    """
    Routine to find candidates for the maximum rectangle for a cell
    based on rectangles of neighboring cells.
    """
    rect_candidates = Array{Int64, 2}(
        undef,
        sum(neighbor->length(neighbor.rect_list), neighbors),
        length(neighbors)
    )
    rect_candidates_points = []
    candidate_idx = 1

    # Neighbors correspond to neighboring dimensions in decreasing order,
    # e.g. z, y, x in the three dimensional case.
    for (i_neighbor, neighbor) in enumerate(neighbors)
        # i_neighbor iterates grid indices
        # dim iterates dimensions
        dim = length(neighbors) - i_neighbor + 1
        for (rect_idx, rect) in enumerate(neighbor.rect_list)
            # Construct candidate based on neighbor and current cell
            local_rect::Vector{Int64} = [
                (
                    i_dim == dim ?
                    # At the dimension of the neighboring rectangle, the size just increments by 1
                    rect[i_dim] + 1 :
                    # Otherwise, we have to check how far the candidate can be spanned
                    # in the direction of the corresponding dimension
                    min(grid[grid_idx].dist_to_zero[i_dim], rect[i_dim])
                )
                # Dimensional iteration
                for i_dim in 1:length(neighbors)
            ]
            # Points of candidate include all points of the regarded neighboring rectangle
            local_points = neighbor.point_list_rects[rect_idx]
            # Starting index of traceback procedure to collect all additional points
            tb_start = CartesianIndex(
                [
                    (
                        i_index == i_neighbor ?
                        # At the neighboring index, the traceback index equals the grid index
                        grid_idx[i_index] :
                        # The other dimensions have to be iterated.
                        # Starting from the grid index, we go back as far as possible,
                        # limited by dist_to_zero and the neighboring rectangle
                        grid_idx[i_index] - min(
                            grid[grid_idx].dist_to_zero[length(neighbors) - i_index + 1],
                            rect[length(neighbors) - i_index + 1]
                        ) + 1
                    )
                    # Grid index iteration
                    for i_index in 1:length(neighbors)
                ]...
            )
            local_points = vcat(
                local_points,
                [grid[tb_idx].point_list for tb_idx in tb_start:grid_idx]...
            )

            push!(rect_candidates_points, local_points)
            rect_candidates[candidate_idx, :] = local_rect
            candidate_idx += 1
        end
    end
    return rect_candidates, rect_candidates_points
end

function max_rectangle(data::Array{Float64}, cellsize::Float64, dim::Int64)
    """
    Multidimensional implementation of algorithm to find the maximum rectangle in a grid.
    The data is expected to be normalized.
    """
    nb_cells = convert(Int64, ceil(1 / cellsize))

    # Build grid based on nb_cells and dim
    grid = [
        Cell(dim=dim) for __ in CartesianIndices(Base.OneTo.(Tuple([nb_cells for __ in 1:dim])))
    ]

    # Store points in corresponding grid cells
    for point_idx in 1:size(data, 1)
        grid_idx = CartesianIndex(
            # List comprehension of grid coordinates
            [
                convert(Int64, floor(data[point_idx, dim_idx] / cellsize)) + 1
                for dim_idx in reverse(1:dim)
            ]... # Splat list components for usage as arguments for constructing CartesianIndex
        )
        push!(grid[grid_idx].point_list, data[point_idx, :])
    end

    global_max = 0
    progress = Progress(nb_cells ^ dim, 1)
    grid_indices = CartesianIndices(grid)
    for grid_idx in grid_indices
        if !isempty(grid[grid_idx].point_list)
            # Admittedly brainfuck-ish
            grid[grid_idx].dist_to_zero = [
                # Construct each element of dist_to_zero independently
                # by using the dist_to_zero of previous cells in each dimension.
                # Note the reversed order of storage of the coordinates in grid,
                # e.g., grid[z, y, x] in three dimensions.
                grid[
                    CartesianIndex(
                        [
                            idx_dim == dim_idx ? max(idx_value - 1, 1) : idx_value
                            for (idx_dim, idx_value) in enumerate(Tuple(grid_idx))
                        ]...
                    )
                ].dist_to_zero[dim - dim_idx + 1] + 1
                for dim_idx in reverse(1:dim)
            ]
            row = grid_idx[1]
            col = grid_idx[2]
            # Check if all of the dimensional components of the index are 1 except for one.
            # Then, max_rect and point_list just accumulate along the dimension which is not 1.
            if count(x->x == 1, Tuple(grid_idx)) >= dim - 1
                # dist_to_zero is already calculated
                grid[grid_idx].max_rect = grid[grid_idx].dist_to_zero
                push!(grid[grid_idx].rect_list, grid[grid_idx].dist_to_zero)
                point_list = grid[grid_idx].point_list
                if grid[grid_idx].dist_to_zero != [1 for __ in 1:dim]
                    # point_list is the accumulation of point_list at current cell and point_list
                    # of previous cell along the dimension whose index component is not 1
                    point_list = vcat(
                        point_list,
                        grid[
                            # max is evaluated for each index component independently.
                            # Due to only being one component differing from 1,
                            # the index is always reduced along the desired dimension.
                            max(grid_idx - oneunit(grid_idx), first(grid_indices))
                        ].point_list_max_rect
                    )
                end
                push!(grid[grid_idx].point_list_rects, point_list)
                grid[grid_idx].point_list_max_rect = point_list
            else
                # If there is more than one dimension differing from 1, the current cell is
                # surrounded by other cells.
                # The max_rect of the current cell is estimated as a combination of a rect from
                # the rect_list of non-overlapping rects of a neighboring cell
                # and the possibilities, which the inclusion of the current rect offers.
                neighboring_cells = Array{Cell, 1}(undef, dim)
                for neighbor_idx in 1:dim
                    neighboring_cells[neighbor_idx] = grid[
                        # Reduce the dimensional components of grid_idx by 1 one ofter the other
                        # to get neighboring cells.
                        max(
                            grid_idx - CartesianIndex(
                                (neighbor_idx == dim_idx ? 1 : 0 for dim_idx in 1:dim)...
                            ),
                            CartesianIndex((1 for __ in 1:dim)...)
                        )
                    ]
                end
                if all(neighbor->isempty(neighbor.rect_list), neighboring_cells)
                    # If rect_list of all neighbors are empty, the current cell is the start
                    # of a new rect.
                    grid[grid_idx].max_rect = [1 for __ in 1:dim]
                    push!(grid[grid_idx].rect_list, [1 for __ in 1:dim])
                    push!(grid[grid_idx].point_list_rects, grid[grid_idx].point_list)
                    grid[grid_idx].point_list_max_rect = grid[grid_idx].point_list
                else
                    max_rect = [0 for __ in 1:dim]
                    point_list_max_rect = []
                    # Find all rectangle candidates using all neighboring cells
                    candidates, candidate_points = find_candidates(
                        neighboring_cells,
                        grid,
                        grid_idx,
                        cellsize
                    )

                    for candidate_idx in 1:size(candidates, 1)
                        candidate = candidates[candidate_idx, :]
                        if prod(candidate) > prod(max_rect)
                            max_rect = candidate
                            point_list_max_rect = candidate_points[candidate_idx]
                        end
                        # Only rectangles, which are not completely overlapped by other
                        # rectangles should be remembered.
                        # So the overlapping between the candidate and the rectangles of
                        # the current cell has to be checked.
                        (
                            # If rect_list is empty, the candidate can be appended to it,
                            # since it will be the only member.
                            isempty(grid[grid_idx].rect_list) ?
                            append_rect = true :
                            append_rect = false
                        )
                        not_overlapped = ones(length(grid[grid_idx].rect_list))
                        for rect_idx in 1:size(grid[grid_idx].rect_list, 1)
                            rect = grid[grid_idx].rect_list[rect_idx]
                            # If a lager rectangle if found,
                            # no other members of rect_list have to be checked
                            if all(rect .>= candidate)
                                append_rect = false
                                break
                            # If the candidate is larger than a member,
                            # the member has to be removed from rect_list
                            elseif all(rect .<= candidate)
                                not_overlapped[rect_idx] = 0
                            end
                            append_rect = true
                        end
                        if any(not_overlapped .!= 1)
                            grid[grid_idx].rect_list = [
                                grid[grid_idx].rect_list[rect_idx]
                                for rect_idx in 1:size(not_overlapped, 1)
                                if not_overlapped[rect_idx] == 1
                            ]
                            grid[grid_idx].point_list_rects = [
                                grid[grid_idx].point_list_rects[rect_idx]
                                for rect_idx in 1:size(not_overlapped, 1)
                                if not_overlapped[rect_idx] == 1
                            ]
                        end
                        if append_rect
                            push!(
                                grid[grid_idx].point_list_rects,
                                candidate_points[candidate_idx]
                            )
                            push!(grid[grid_idx].rect_list, candidate)
                        end
                    end
                    grid[grid_idx].max_rect = max_rect
                    grid[grid_idx].point_list_max_rect = point_list_max_rect
                end
            end
            # Adjust the global maximal product of rectangle components for the colorscale limit
            rect_size = prod(grid[grid_idx].max_rect)
            if rect_size > global_max
                global_max = rect_size
            end
        end
        next!(progress)
    end
    colormap = cmap.get_cmap("viridis")
    normalize = colors.Normalize(vmin=0, vmax=global_max)
    for grid_idx in CartesianIndices(grid)
        rect_size = prod(grid[grid_idx].max_rect)
        if rect_size != 0
            grid[grid_idx].color = colormap(normalize(rect_size))
        end
    end

    return grid
end

end

max_stars_repo_name: effaeff/distribution-sampling-jl

max_issues_repo_name: effaeff/distribution-sampling-jl

max_forks_repo_name: effaeff/distribution-sampling-jl


text:
import LinearAlgebra: wrapperop

for AT in [:AdjOrTransAbsHalfIntMatrix, :AdjOrTransAbsHalfIntVector]
	for IT in [:Integer, :Real]
		@eval function Base.isassigned(A::$AT, I::$IT...)
			checkbounds(Bool, A, I...) || return false
		    K = trimtoN(I,Val(2))
		    isassigned(parent(A), reverse(K)...)
		end
	end
	for IT in [:Int, :Real]
		@eval @propagate_inbounds function Base.getindex(A::$AT, i::$IT, j::$IT)
			@boundscheck checkbounds(A, i, j)
			wrapperop(A)(A.parent[j,i])
		end

		@eval @propagate_inbounds function Base.setindex!(A::$AT, val, i::$IT, j::$IT)
			@boundscheck checkbounds(A, i, j)
			A.parent[j,i] = wrapperop(A)(val)
			A
		end
	end
end

# AbstractArray interface, additional definitions to retain wrapper over vectors where appropriate
for DT in [:Int, :Real]
	@eval @propagate_inbounds Base.getindex(v::AdjOrTransAbsHalfIntVector, ::Colon, is::AbstractArray{$DT}) = wrapperop(v)(v.parent[is])
end
@propagate_inbounds Base.getindex(v::AdjOrTransAbsHalfIntVector, ::Colon, ::Colon) = wrapperop(v)(copy(v.parent))

@propagate_inbounds function Base.getindex(a::AdjOrTransAbsHalfIntVecOrMat, i::Int)
	@boundscheck checkbounds(a, i)
	c = CartesianIndicesHalfInt(axes(a))
	inds = Tuple(c[i])
	a[inds...]
end

@propagate_inbounds function Base.setindex!(a::AdjOrTransAbsHalfIntVecOrMat, val, i::Int)
	@boundscheck checkbounds(a, i)
	c = CartesianIndicesHalfInt(axes(a))
	inds = Tuple(c[i])
	a[inds...] = val
	a
end

function LinearAlgebra.diagind(h::AbstractHalfIntegerMatrix, k::Integer=0)
	diagind(size(h,1),size(h,1),k)
end
# These methods fall back to the parent by default, assuming that parent(x) != x
# Subtypes should specialize them if custom behaviour is desired
for f in [:isdiag, :isbanded, :issymmetric, :ishermitian, :isposdef, :det]
	@eval LinearAlgebra.$f(h::AbstractHalfIntegerMatrix) = $f(parent(h))
end

for f in [:istriu, :istril]
	@eval LinearAlgebra.$f(h::AbstractHalfIntegerMatrix, k::Integer=0) = $f(parent(h), k)
end
LinearAlgebra.diag(h::AbstractHalfIntegerMatrix, k::Integer=0) = parent(h)[diagind(h,k)]

function LinearAlgebra.inv(h::HalfIntArray)
	HalfIntArray(inv(parent(h)), h.offsets)
end
function LinearAlgebra.inv(s::SpinMatrix)
	SpinMatrix(inv(parent(s)), s.j)
end

function Base.:(*)(h1::HalfIntMatrix, h2::HalfIntMatrix)
	axes(h1,2) == axes(h2,1) ||
	throw(DimensionMismatch("A has dimension $(axes(h1)) but B has dimension $(axes(h2))"))

	y1, y2 = map(unwraphalfint,(h1,h2))
	HalfIntArray(y1 * y2, h1.offsets)
end
function Base.:(*)(S1::SpinMatrix, S2::SpinMatrix)
	S1.j == S2.j ||
	throw(DimensionMismatch("A has dimension $(axes(S1)) but B has dimension $(axes(S2))"))

	y1, y2 = map(unwraphalfint,(S1,S2))
	SpinMatrix(y1 * y2, S1.j)
end


text:
export preprocess_split

using mlpack._Internal.io

import mlpack_jll
const preprocess_splitLibrary = mlpack_jll.libmlpack_julia_preprocess_split

# Call the C binding of the mlpack preprocess_split binding.
function preprocess_split_mlpackMain()
  success = ccall((:preprocess_split, preprocess_splitLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module preprocess_split_internal
  import ..preprocess_splitLibrary

end # module

"""
    preprocess_split(input; [input_labels, no_shuffle, seed, test_ratio, verbose])

This utility takes a dataset and optionally labels and splits them into a
training set and a test set. Before the split, the points in the dataset are
randomly reordered. The percentage of the dataset to be used as the test set can
be specified with the `test_ratio` parameter; the default is 0.2 (20%).

The output training and test matrices may be saved with the `training` and
`test` output parameters.

Optionally, labels can be also be split along with the data by specifying the
`input_labels` parameter.  Splitting labels works the same way as splitting the
data. The output training and test labels may be saved with the
`training_labels` and `test_labels` output parameters, respectively.

So, a simple example where we want to split the dataset `X` into `X_train` and
`X_test` with 60% of the data in the training set and 40% of the dataset in the
test set, we could run

```julia
julia> using CSV
julia> X = CSV.read("X.csv")
julia> X_test, _, X_train, _ = preprocess_split(X; test_ratio=0.4)
```

Also by default the dataset is shuffled and split; you can provide the
`no_shuffle` option to avoid shuffling the data; an example to avoid shuffling
of data is:

```julia
julia> using CSV
julia> X = CSV.read("X.csv")
julia> X_test, _, X_train, _ = preprocess_split(X; no_shuffle=1,
            test_ratio=0.4)
```

If we had a dataset `X` and associated labels `y`, and we wanted to split these
into `X_train`, `y_train`, `X_test`, and `y_test`, with 30% of the data in the
test set, we could run

```julia
julia> using CSV
julia> X = CSV.read("X.csv")
julia> y = CSV.read("y.csv"; type=Int)
julia> X_test, y_test, X_train, y_train = preprocess_split(X;
            input_labels=y, test_ratio=0.3)
```

# Arguments

 - `input::Array{Float64, 2}`: Matrix containing data.
 - `input_labels::Array{Int, 2}`: Matrix containing labels.
 - `no_shuffle::Bool`: Avoid shuffling and splitting the data.  Default
      value `false`.

 - `seed::Int`: Random seed (0 for std::time(NULL)).  Default value `0`.

 - `test_ratio::Float64`: Ratio of test set; if not set,the ratio defaults
      to 0.2  Default value `0.2`.

 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `test::Array{Float64, 2}`: Matrix to save test data to.
 - `test_labels::Array{Int, 2}`: Matrix to save test labels to.
 - `training::Array{Float64, 2}`: Matrix to save training data to.
 - `training_labels::Array{Int, 2}`: Matrix to save train labels to.

"""
function preprocess_split(input;
                          input_labels = missing,
                          no_shuffle::Union{Bool, Missing} = missing,
                          seed::Union{Int, Missing} = missing,
                          test_ratio::Union{Float64, Missing} = missing,
                          verbose::Union{Bool, Missing} = missing,
                          points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, preprocess_splitLibrary), Nothing, ());

  IORestoreSettings("Split Data")

  # Process each input argument before calling mlpackMain().
  IOSetParamMat("input", input, points_are_rows)
  if !ismissing(input_labels)
    IOSetParamUMat("input_labels", input_labels, points_are_rows)
  end
  if !ismissing(no_shuffle)
    IOSetParam("no_shuffle", convert(Bool, no_shuffle))
  end
  if !ismissing(seed)
    IOSetParam("seed", convert(Int, seed))
  end
  if !ismissing(test_ratio)
    IOSetParam("test_ratio", convert(Float64, test_ratio))
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("test")
  IOSetPassed("test_labels")
  IOSetPassed("training")
  IOSetPassed("training_labels")
  # Call the program.
  preprocess_split_mlpackMain()

  return IOGetParamMat("test", points_are_rows),
         IOGetParamUMat("test_labels", points_are_rows),
         IOGetParamMat("training", points_are_rows),
         IOGetParamUMat("training_labels", points_are_rows)
end


text:
#=
Copyright 2020 Ali Mert Ceylan

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
=#

module AutoDiff

using NaNMath
using SpecialFunctions
using DiffRules

include("core.jl")
include("interface.jl")

export
    BD,
    convert,
    promote_rule,
    ndims,
    size,
    getindex,
    setindex!,
    length,
    similar,
    broadcasted,
    +,
    -,
    *,
    /,
    ^,
    sum,
    mean,
    abs2,
    transpose,
    adjoint,
    inv,
    det,
    params,
    @differentiate

end # module


text:
using StaticArrays

struct Cylinder{T}
    midpoint::SVector{3,T}
    rot::SMatrix{3,3,T,9}
    r::T
    h::T
end

struct CylinderGenerator{Fr,Fh}
    r::Fr
    h::Fh
end

function (g::CylinderGenerator)(midpoint, rot)
    Cylinder(midpoint, rot, g.r(), g.h())
end

volume(shape::Cylinder) = π * shape.r^2 * h

function set_midpoint(shape::Cylinder, midpoint)
    return Cylinder(midpoint, t.rot, t.r, t.h)
end

function bounding_box(shape::Cylinder)
    lower = (shape.midpoint + shape.rot * SA[0, 0, -shape.h/2]) .- shape.r
    upper = (shape.midpoint + shape.rot * SA[0, 0,  shape.h/2]) .+ shape.r
    return lower, upper
end

function inside(shape::Cylinder, p0)
    p = shape.rot' * (p0 - shape.midpoint)
    return -shape.h/2 < p[3] < shape.h/2 && shape.r^2 < (p[1]^2 + p[2]^2)
end


text:
module CoupledSystems

using LinearAlgebra
using SparseArrays, SparseDiffTools
using FiniteDiff, ForwardDiff, ReverseDiff, DiffResults
using NLsolve

import Base.setindex, Base.front, Base.tail

export @var
export NamedVar
export name
export value
export combine
export combine!
export separate

export AbstractDiffMethod
export AbstractAD, ForwardAD, ReverseAD
export AbstractFD, ForwardFD, CentralFD, ComplexFD
export Forward, Reverse
export Direct, Adjoint

export AbstractSparsityPattern, DensePattern, SparsePattern

export AbstractComponent
export AbstractExplicitComponent, ExplicitComponent, ExplicitSystem
export AbstractImplicitComponent, ImplicitComponent, ImplicitSystem

export AbstractSolver, Newton

export outputs, outputs!, outputs!!, outputs!!!
export jacobian, jacobian!, jacobian!!, jacobian!!!
export outputs_and_jacobian, outputs_and_jacobian!, outputs_and_jacobian!!, outputs_and_jacobian!!!
export residuals, residuals!, residuals!!, residuals!!!
export residual_input_jacobian, residual_input_jacobian!, residual_input_jacobian!!, residual_input_jacobian!!!
export residual_output_jacobian, residual_output_jacobian!, residual_output_jacobian!!, residual_output_jacobian!!!
export residuals_and_input_jacobian, residuals_and_input_jacobian!, residuals_and_input_jacobian!!, residuals_and_input_jacobian!!!
export residuals_and_output_jacobian, residuals_and_output_jacobian!, residuals_and_output_jacobian!!, residuals_and_output_jacobian!!!
export residuals_and_jacobians, residuals_and_jacobians!, residuals_and_jacobians!!, residuals_and_jacobians!!!

# named variables
include("variables.jl")

# functions for computing derivatives
include("derivatives.jl")

# definition of components and calling functions
include("components.jl")

# functions for finding jacobian matrix sparsity
include("sparsity.jl")

# constructors for explicit and implicit components
include("constructors.jl")

# collection of built-in solvers for converting implicit components to explicit components
include("solvers.jl")

end


text:
using Base.Test
using Elemental

m, n = 10, 10

for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
    for mat in (Elemental.Matrix, Elemental.DistMatrix)
        A = mat(T)
        Elemental.gaussian!(A, m, n)

        @test norm(A, 1) ≈ norm(T[A[i,j] for i = 1:m, j = 1:n], 1)
        @test norm(A) ≈ norm(T[A[i,j] for i = 1:m, j = 1:n])
        @test norm(A, Inf) ≈ norm(T[A[i,j] for i = 1:m, j = 1:n], Inf)
        @test countnz(A) ≈ countnz(T[A[i,j] for i = 1:m, j = 1:n])
    end
end


text:
## chernoff.jl
##
## The code below is intended to go with the Distributions package of Julia.
## It was written by Joris Pinkse, joris@psu.edu, on December 12, 2017.  Caveat emptor.
##
## It computes pdf, cdf, moments, quantiles, and random numbers for the Chernoff distribution.
##
## Most of the symbols have the same meaning as in the Groeneboom and Wellner paper.
##
## Random numbers are drawn using a Ziggurat algorithm.  To obtain draws in the tails, the
## algorithm reverts to quantiles, which is slow.
"""
    Chernoff()

The *Chernoff distribution* is the distribution of the random variable
```math
\\argmax_{t \\in (-\\infty,\\infty)} ( G(t) - t^2 ),
```
where ``G`` is standard two--sided Brownian motion.

The distribution arises as the limit distribution of various cube-root-n consistent estimators,
including the isotonic regression estimator of Brunk, the isotonic density estimator of Grenander,
the least median of squares estimator of Rousseeuw, and the maximum score estimator of Manski.

For theoretical results, see e.g. Kim and Pollard, Annals of Statistics, 1990.  The code for the
computation of pdf and cdf is based on the algorithm described in Groeneboom and Wellner,
Journal of Computational and Graphical Statistics, 2001.

```julia
Chernoff()
pdf(Chernoff(),x::Real)
cdf(Chernoff(),x::Real)
logpdf(Chernoff(),x::Real)
survivor(Chernoff(),x::Real)
mean(Chernoff())
var(Chernoff())
skewness(Chernoff())
kurtosis(Chernoff())
kurtosis(Chernoff(), excess::Bool)
mode(Chernoff())
entropy(Chernoff())
rand(Chernoff())
rand(rng, Chernoff()
cdf(Chernoff(),-x)              #For tail probabilities, use this instead of 1-cdf(Chernoff(),x)
```
"""
struct Chernoff <: ContinuousUnivariateDistribution
end

module ChernoffComputations
    import QuadGK.quadgk
    # The following arrays of constants have been precomputed to speed up computation.
    # The first two correspond to the arrays a and b in the Groeneboom and Wellner article.
    # The array airyai_roots contains roots of the airyai functions (atilde in the paper).
    # Finally airyai_prime contains contains airyaiprime evaluated at the corresponding atildes
    const a=[
         0.14583333333333334
        -0.0014105902777777765
        -2.045269786155239e-5
         1.7621771897199738e-6
        -5.7227597511834636e-8
         1.2738111577329235e-9
        -2.1640824986157608e-11
         2.87925759635786e-13
        -2.904496582565578e-15
         1.8213982724416428e-17
         4.318809086319947e-20
        -3.459634830348895e-21
         6.541295360751684e-23
        -8.92275524974091e-25
         1.0102911018761945e-26
        -9.976759882474874e-29
        ]
    const b=[
         0.6666666666666666
         0.021164021164021163
        -0.0007893341226674574
         1.9520978781876193e-5
        -3.184888134149933e-7
         2.4964953625552273e-9
         3.6216439304006735e-11
        -1.874438299727797e-12
         4.393984966701305e-14
        -7.57237085924526e-16
         1.063018673013022e-17
        -1.26451831871265e-19
         1.2954000902764143e-21
        -1.1437790369170514e-23
         8.543413293195206e-26
        -5.05879944592705e-28
        ]
    const airyai_roots=[
        -2.338107410459767
        -4.08794944413097
        -5.520559828095551
        -6.786708090071759
        -7.944133587120853
        -9.02265085334098
        -10.040174341558085
        -11.008524303733262
        -11.936015563236262
        -12.828776752865757
        -13.691489035210719
        -14.527829951775335
        -15.340755135977997
        -16.132685156945772
        -16.90563399742994
        -17.66130010569706
        ]
    const airyai_prime=[
         0.7012108227206906
        -0.8031113696548534
         0.865204025894141
        -0.9108507370495821
         0.9473357094415571
        -0.9779228085695176
         1.0043701226603126
        -1.0277386888207862
         1.0487206485881893
        -1.0677938591574279
         1.0853028313507
        -1.1015045702774968
         1.1165961779326556
        -1.1307323104931881
         1.1440366732735523
        -1.156609849116565
        ]

    const cuberoottwo = cbrt(2.0)
    const sqrthalfpi = sqrt(0.5*pi)
    const sqrttwopi = sqrt(2.0*pi)

    function p(y::Real)
        if iszero(y)
            return -sqrt(0.5*pi)
        end

        (y > 0) || throw(DomainError(y, "argument must be positive"))

        cnsty = y^(-1.5)
        if (y <= 1.0)
            return sum([(b[k]*cnsty - a[k]*sqrthalfpi)*y^(3*k) for k=1:length(a)])-sqrthalfpi
        else
            return sum([exp(cuberoottwo*a*y) for a in airyai_roots]) * 2 * sqrttwopi * exp(-y*y*y/6) - cnsty
        end
    end

    function g(x::Real)
        function g_one(y::Real)
            return p(y) * exp(-0.5*y*(2*x+y)*(2*x+y))
        end
        function g_two(y::Real)
            z = 2*x+y*y
            return (z*y*y + 0.5 * z*z) * exp(-0.5*y*y*z*z)
        end
        if (x <= -1.0)
            return cuberoottwo*cuberoottwo * exp(2*x*x*x/3.0) * sum([exp(-cuberoottwo*airyai_roots[k]*x) / airyai_prime[k] for k=1:length(airyai_roots)])
        else
            return 2*x - (quadgk(g_one, 0.0, Inf)[1] - 4*quadgk(g_two, 0.0, Inf)[1]) / sqrttwopi   # should perhaps combine integrals
        end
    end

    _pdf(x::Real) = g(x)*g(-x)*0.5
    _cdf(x::Real) = (x < 0.0) ? _cdfbar(-x) : 0.5 + quadgk(_pdf,0.0,x)[1]
    _cdfbar(x::Real) = (x < 0.0) ? _cdf(x) : quadgk(_pdf, x, Inf)[1]
end

pdf(d::Chernoff, x::Real) = ChernoffComputations._pdf(x)
logpdf(d::Chernoff, x::Real) = log(ChernoffComputations.g(x))+log(ChernoffComputations.g(-x))+log(0.5)
cdf(d::Chernoff, x::Real) = ChernoffComputations._cdf(x)

function quantile(d::Chernoff, tau::Real)
    # some commonly used quantiles were precomputed
    precomputedquants=[
        0.0 -Inf;
        0.01 -1.171534341573129;
        0.025 -0.9981810946684274;
        0.05 -0.8450811886357725;
        0.1 -0.6642351964332931;
        0.2 -0.43982766604886553;
        0.25 -0.353308035220429;
        0.3 -0.2751512847290148;
        0.4 -0.13319637678583637;
        0.5 0.0;
        0.6 0.13319637678583637;
        0.7 0.2751512847290147;
        0.75 0.353308035220429;
        0.8 0.4398276660488655;
        0.9 0.6642351964332931;
        0.95 0.8450811886357724;
        0.975 0.9981810946684272;
        0.99 1.17153434157313;
        1.0 Inf
        ]
    (0.0 <= tau && tau <= 1.0) || throw(DomainError(tau, "illegal value of tau"))
    present = searchsortedfirst(precomputedquants[:, 1], tau)
    if present <= size(precomputedquants, 1)
        if tau == precomputedquants[present, 1]
            return precomputedquants[present, 2]
        end
    end

    # one good approximation of the quantiles can be computed using Normal(0.0, stdapprox) with stdapprox = 0.52
    stdapprox = 0.52
    dnorm = Normal(0.0, 1.0)
    if tau < 0.001
        return -newton(x -> tau - ChernoffComputations._cdfbar(x), ChernoffComputations._pdf, quantile(dnorm, 1.0 - tau)*stdapprox)

    end
    if tau > 0.999
        return newton(x -> 1.0 - tau - ChernoffComputations._cdfbar(x), ChernoffComputations._pdf, quantile(dnorm, tau)*stdapprox)
    end
    return newton(x -> ChernoffComputations._cdf(x) - tau, ChernoffComputations._pdf, quantile(dnorm, tau)*stdapprox)   # should consider replacing x-> construct for speed
end

minimum(d::Chernoff) = -Inf
maximum(d::Chernoff) = Inf
insupport(d::Chernoff, x::Real) = isnan(x) ? false : true

mean(d::Chernoff) = 0.0
var(d::Chernoff) = 0.26355964132470455
modes(d::Chernoff) = [0.0]
mode(d::Chernoff) = 0.0
skewness(d::Chernoff) = 0.0
kurtosis(d::Chernoff) = -0.16172525511461888
kurtosis(d::Chernoff, excess::Bool) = kurtosis(d) + (excess ? 0.0 : 3.0)
entropy(d::Chernoff) = -0.7515605300273104

### Random number generation
rand(d::Chernoff) = rand(GLOBAL_RNG, d)
function rand(rng::AbstractRNG, d::Chernoff)                 # Ziggurat random number generator --- slow in the tails
    # constants needed for the Ziggurat algorithm
    A = 0.03248227216266608
    x = [
        1.4765521793744492
        1.3583996502410562
        1.2788224934376338
        1.2167121025431031
        1.164660153310361
        1.1191528874523227
        1.0782281238946987
        1.0406685077375248
        1.0056599129954287
        0.9726255909850547
        0.9411372703351518
        0.910863725882819
        0.8815390956471935
        0.8529422519848634
        0.8248826278765808
        0.7971898990526088
        0.769705944382039
        0.7422780374708061
        0.7147524811697309
        0.6869679724643997
        0.6587479056872714
        0.6298905501492661
        0.6001554584431008
        0.5692432986584453
        0.5367639126935895
        0.5021821750911241
        0.4647187417226889
        0.42314920361072667
        0.37533885097957154
        0.31692143952775814
        0.2358977457249061
        1.0218214689661219e-7
       ]
    y=[
        0.02016386420423385
        0.042162593823411566
        0.06607475557706186
        0.09147489698007219
        0.11817165794330549
        0.14606157249935905
        0.17508555351826158
        0.20521115629454248
        0.2364240467790298
        0.26872350680813245
        0.3021199879039766
        0.3366338395880132
        0.37229479658922043
        0.4091420247190392
        0.4472246406136998
        0.48660269400571066
        0.5273486595545326
        0.5695495447104147
        0.6133097942704644
        0.6587552778727138
        0.7060388099508409
        0.7553479187650504
        0.8069160398536568
        0.8610391369707446
        0.9181013320624011
        0.97861633948841
        1.0432985878313983
        1.113195213997703
        1.1899583790598025
        1.2764995726449935
        1.3789927085182485
        1.516689116183566
        ]
    n = length(x)
    i = rand(rng, 0:n-1)
    r = (2.0*rand(rng)-1) * ((i>0) ? x[i] : A/y[1])
    rabs = abs(r)
    if rabs < x[i+1]
        return r
    end
    s = (i>0) ? (y[i]+rand(rng)*(y[i+1]-y[i])) : rand(rng)*y[1]
    if s < 2.0*ChernoffComputations._pdf(rabs)
        return r
    end
    if i > 0
        return rand(rng, d)
    end
    F0 = ChernoffComputations._cdf(A/y[1])
    tau = 2.0*rand(rng)-1 # ~ U(-1,1)
    tauabs = abs(tau)
    return quantile(d, tauabs + (1-tauabs)*F0) * sign(tau)
end