text:
"""
    const sorting_rules = Dict{Symbol, Vector{Symbol}}(
        :nty  => [:authors;:editors;:title;:date],
        :nyt  => [:authors;:editors;:date;:title],
        :y    => [:date]
    );

Implemented sorting rules for bibliography entry sorting.

See also [`sort_bibliography!`](@ref).
"""
const sorting_rules = Dict{Symbol,Vector{Symbol}}(
    :nty => [:authors; :editors; :title; :date],
    :nyt => [:authors; :editors; :date; :title],
    :y => [:date],
)

"""
    sort_bibliography!(
        bibliography::DataStructures.OrderedDict{String,Entry},
        sorting_rule::Symbol = :key
        )

Sorts the bibliography in place.

The sorting order can be set by specifying the `sorting_rule`. The sorting is
implemented via `isless()` functions. For detailed insight have a look at the
`isless()` implementation of Julia and `BibInternal.jl`.

Supported symbols for `sorting_rule` are:
- `:key` __(default)__: sort by bibliography keys e.g. BibTeX keys or `:id`
- the sorting rules defined in [`sorting_rules`](@ref)

!!! info "Note:"
    The sorting is not following explicitly bibliographic alphabetizing
    conventions. It follows standard comparator behaviour implied by the
    implemented `isless()` functions (string comparators).
"""
function sort_bibliography!(
    bibliography::DataStructures.OrderedDict{String,Entry}, sorting_rule::Symbol=:key
)
    # TODO: allow Union{Symbol,Vector{Symbol}} for a arbitrary custom sorting order
    #       this needs one additional type check and a check for allowed symbols in
    #       vector (check against fieldnames(Bibliography.BibInternal.Entry) )
    if sorting_rule == :key
        sort!(bibliography)
    elseif sorting_rule in keys(sorting_rules)
        sort!(
            bibliography;
            lt=(a, b) -> recursive_isless(a, b, sorting_rules[sorting_rule]),
            by=x -> bibliography[x],
        )
    else
        throw(ArgumentError("Unsupported sorting order!"))
    end
end

"""
    recursive_isless(a::Entry, b::Entry, fields::Tuple{Symbol},
                     depth::Int = 0)

Helper function for [`sort_bibliography!`](@ref).

This function allows recursive checking if `a < b` with descending importance.
The importance set for the comparison is defined by the argument `fields`. This
argument is a tuple consisting of symbols denoting the fields of the data type
`BibInternal.Entry`. The ordering implies the importance.

The `depth` argument is purely for iterating/recursive purposes.
"""
function recursive_isless(a::Entry, b::Entry, fields::Vector{Symbol}, depth::Int=1)::Bool
    depth > length(fields) && return false
    a_field = getfield(a, fields[depth])
    b_field = getfield(b, fields[depth])

    if a_field == b_field
        return recursive_isless(a, b, fields, depth + 1)
    end
    return a_field < b_field
end


text:
import NLPModels

export NLPOptiModel

#export obj, grad!, cons!, jprod!, jtprod!, proj!, prox!, objprox!

#
#  Given an AbstractNLPModel ´nlp´ representing the problem
#
#    optimize   f(x)
#    s.t.       xl <=   x  <= xu
#               cl <= c(x) <= cu
#
#  the call
#
#    prob = NLPOptiModel( nlp )
#
#  returns an NLPOptiModel <: AbstractOptiModel ´prob´ representing the same problem
#  in the form
#
#    optimize   f(x) + g(x)
#    s.t.       c(x) in S
#
#  where
#
#    g = indBox( xl, xu )
#    S = Box( cl, cu )
#
#  In addition to keeping `meta` as any OptiModel, an NLPOptiModel
#  also stores ´base´, the original problem ´nlp´.

# Type
mutable struct NLPOptiModel{M <: NLPModels.AbstractNLPModel} <: AbstractOptiModel
  meta::OptiModelMeta
  base::M
end

# Constructor
function NLPOptiModel( nlp::M; name=nlp.meta.name * "-nlpopti" ) where {M <: NLPModels.AbstractNLPModel}

  meta = OptiModelMeta(nlp.meta.nvar, nlp.meta.ncon, x0=nlp.meta.x0, y0=nlp.meta.y0, minimize=nlp.meta.minimize, name=name)

  return NLPOptiModel(meta, nlp)
end

##########################################################
# Necessary methods
##########################################################
# obj       : unchanged
# grad!     : unchanged
# cons!     : unchanged
# jprod!    : unchanged
# jtprod!   : unchanged
# proj!     : projection onto [cl, cu]
# prox!     : projection onto [xl, xu]
# objprox!  : projection onto [xl, xu], zero cost

function obj(prob::NLPOptiModel, x::AbstractVector)
  return NLPModels.obj(prob.base, x)
end

function grad!(prob::NLPOptiModel, x::AbstractVector, dfx::AbstractVector)
  return NLPModels.grad!(prob.base, x, dfx)
end

function cons!(prob::NLPOptiModel, x::AbstractVector, cx::AbstractVector)
  return NLPModels.cons!(prob.base, x, cx)
end

function jprod!(prob::NLPOptiModel, x::AbstractVector, v::AbstractVector, Jv::AbstractVector)
  return NLPModels.jprod!(prob.base, x, v, Jv)
end

function jtprod!(prob::NLPOptiModel, x::AbstractVector, v::AbstractVector, Jtv::AbstractVector)
  return NLPModels.jtprod!(prob.base, x, v, Jtv)
end

function proj!(prob::NLPOptiModel, cx::AbstractVector, px::AbstractVector)
  @lencheck prob.meta.ncon cx px
  px .= clamp.( cx, prob.base.meta.lcon, prob.base.meta.ucon)
  return nothing
end

function prox!(prob::NLPOptiModel, x::AbstractVector, a::Real, z::AbstractVector)
  @lencheck prob.meta.nvar x z
  z .= clamp.( x, prob.base.meta.lvar, prob.base.meta.uvar)
  return nothing
end

function objprox!(prob::NLPOptiModel, x::AbstractVector, a::Real, z::AbstractVector)
  @lencheck prob.meta.nvar x z
  z .= clamp.( x, prob.base.meta.lvar, prob.base.meta.uvar)
  return 0.0
end

#  for meth in (:obj, :grad!, :cons!, :jprod!, :jtprod!)
#    @eval NLPOptiMo.$meth(prob::NLPOptiModel, ) = NLPModels.$meth(prob.base, )
#  end

ext: jl


text:
using Images, ImageView, ImageSegmentation

# one object
img = load("sample-images/cat-3352842_640.jpg");
imshow(img);

seeds = [
    (CartesianIndex(220,250), 1), # object location
    (CartesianIndex(220,500), 2)  # background location
]

segments = seeded_region_growing(img, seeds)

imshow(map(i->segment_mean(segments,i), labels_map(segments)))

# multiple objects
img = load("sample-images/kittens-555822_640.jpg");
imshow(img);

seeds = [
    (CartesianIndex(130,90), 1),
    (CartesianIndex(130,200), 2),
    (CartesianIndex(130,300), 2),
    (CartesianIndex(130,420), 3),
    (CartesianIndex(30,50), 4)
]

segments = seeded_region_growing(img, seeds)

imshow(map(i->segment_mean(segments,i), labels_map(segments)))

# multiple objects
img = load("sample-images/bird-3183441_640.jpg");

seeds = [
    (CartesianIndex(240,120), 1),
    (CartesianIndex(295,70), 2),
    (CartesianIndex(319,40), 3),
    (CartesianIndex(90,300), 4),
    (CartesianIndex(295,325), 5),
    (CartesianIndex(76,135), 6)
]

segments = seeded_region_growing(img, seeds)

imshow(map(i->segment_mean(segments,i), labels_map(segments)))


text:
export Gaussian, prod!, convert

# Convert parameterizations
function convert(::Type{ProbabilityDistribution{V, GaussianMeanPrecision}}, dist::ProbabilityDistribution{V, GaussianMeanVariance}) where V<:VariateType
    w = cholinv(dist.params[:v])
    m = deepcopy(dist.params[:m])

    return ProbabilityDistribution(V, GaussianMeanPrecision, m=m, w=w)
end

function convert(::Type{ProbabilityDistribution{V, GaussianMeanPrecision}}, dist::ProbabilityDistribution{V, GaussianWeightedMeanPrecision}) where V<:VariateType
    w = deepcopy(dist.params[:w])
    m = cholinv(w)*dist.params[:xi]

    return ProbabilityDistribution(V, GaussianMeanPrecision, m=m, w=w)
end

function convert(::Type{ProbabilityDistribution{V, GaussianMeanVariance}}, dist::ProbabilityDistribution{V, GaussianMeanPrecision}) where V<:VariateType
    v = cholinv(dist.params[:w])
    m = deepcopy(dist.params[:m])

    return ProbabilityDistribution(V, GaussianMeanVariance, m=m, v=v)
end

function convert(::Type{ProbabilityDistribution{V, GaussianMeanVariance}}, dist::ProbabilityDistribution{V, GaussianWeightedMeanPrecision}) where V<:VariateType
    v = cholinv(dist.params[:w])
    m = v*dist.params[:xi]

    return ProbabilityDistribution(V, GaussianMeanVariance, m=m, v=v)
end

function convert(::Type{ProbabilityDistribution{V, GaussianWeightedMeanPrecision}}, dist::ProbabilityDistribution{V, GaussianMeanPrecision}) where V<:VariateType
    w = deepcopy(dist.params[:w])
    xi = w*dist.params[:m]

    return ProbabilityDistribution(V, GaussianWeightedMeanPrecision, xi=xi, w=w)
end

function convert(::Type{ProbabilityDistribution{V, GaussianWeightedMeanPrecision}}, dist::ProbabilityDistribution{V, GaussianMeanVariance}) where V<:VariateType
    w = cholinv(dist.params[:v])
    xi = w*dist.params[:m]

    return ProbabilityDistribution(V, GaussianWeightedMeanPrecision, xi=xi, w=w)
end

# Convert VariateTypes
convert(::Type{ProbabilityDistribution{Multivariate, GaussianMeanVariance}}, dist::ProbabilityDistribution{Univariate, GaussianMeanVariance}) =
    ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[dist.params[:m]], v=mat(dist.params[:v]))
convert(::Type{ProbabilityDistribution{Multivariate, GaussianMeanPrecision}}, dist::ProbabilityDistribution{Univariate, GaussianMeanPrecision}) =
    ProbabilityDistribution(Multivariate, GaussianMeanPrecision, m=[dist.params[:m]], w=mat(dist.params[:w]))
convert(::Type{ProbabilityDistribution{Multivariate, GaussianWeightedMeanPrecision}}, dist::ProbabilityDistribution{Univariate, GaussianWeightedMeanPrecision}) =
    ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[dist.params[:xi]], w=mat(dist.params[:w]))

function prod!(
    x::ProbabilityDistribution{Univariate, F1},
    y::ProbabilityDistribution{Univariate, F2},
    z::ProbabilityDistribution{Univariate, GaussianWeightedMeanPrecision}=ProbabilityDistribution(Univariate, GaussianWeightedMeanPrecision, xi=0.0, w=1.0)) where {F1<:Gaussian, F2<:Gaussian}

    z.params[:xi] = unsafeWeightedMean(x) + unsafeWeightedMean(y)
    z.params[:w] = unsafePrecision(x) + unsafePrecision(y)

    return z
end

@symmetrical function prod!(
    x::ProbabilityDistribution{Univariate, F},
    y::ProbabilityDistribution{Univariate, PointMass},
    z::ProbabilityDistribution{Univariate, PointMass}=ProbabilityDistribution(Univariate, PointMass, m=0.0)) where F<:Gaussian

    z.params[:m] = y.params[:m]
    return z
end

function prod!(
    x::ProbabilityDistribution{Multivariate, F1},
    y::ProbabilityDistribution{Multivariate, F2},
    z::ProbabilityDistribution{Multivariate, GaussianWeightedMeanPrecision}=ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[NaN], w=transpose([NaN]))) where {F1<:Gaussian, F2<:Gaussian}

    z.params[:xi] = unsafeWeightedMean(x) + unsafeWeightedMean(y)
    z.params[:w] = unsafePrecision(x) + unsafePrecision(y)

    return z
end

@symmetrical function prod!(
    x::ProbabilityDistribution{Multivariate, F},
    y::ProbabilityDistribution{Multivariate, PointMass},
    z::ProbabilityDistribution{Multivariate, PointMass}=ProbabilityDistribution(Multivariate, PointMass, m=[NaN])) where F<:Gaussian

    z.params[:m] = deepcopy(y.params[:m])

    return z
end

function sample(dist::ProbabilityDistribution{Univariate, F}) where F<:Gaussian
    isProper(dist) || error("Cannot sample from improper distribution")
    (m,v) = unsafeMeanCov(dist)
    return sqrt(v)*randn() + m
end

function sample(dist::ProbabilityDistribution{Univariate, F}, n_samples::Int64) where F<:Gaussian
    isProper(dist) || error("Cannot sample from improper distribution")
    (m,v) = unsafeMeanCov(dist)

    return sqrt(v).*randn(n_samples) .+ m
end

function sample(dist::ProbabilityDistribution{Multivariate, F}) where F<:Gaussian
    isProper(dist) || error("Cannot sample from improper distribution")
    (m,V) = unsafeMeanCov(dist)
    return (cholesky(default_cholesky_mode, V)).U' *randn(dims(dist)) + m
end

function sample(dist::ProbabilityDistribution{Multivariate, F}, n_samples::Int64) where F<:Gaussian
    isProper(dist) || error("Cannot sample from improper distribution")
    (m,V) = unsafeMeanCov(dist)
    U = (cholesky(default_cholesky_mode, V)).U
    d = dims(dist)

    return [U' *randn(d) + m for i in 1:n_samples]
end

# Entropy functional
function differentialEntropy(dist::ProbabilityDistribution{Univariate, F}) where F<:Gaussian
    return  0.5*log(unsafeCov(dist)) +
            0.5*log(2*pi) +
            0.5
end

function differentialEntropy(dist::ProbabilityDistribution{Multivariate, F}) where F<:Gaussian
    return  0.5*logdet(unsafeCov(dist)) +
            (dims(dist)/2)*log(2*pi) +
            (dims(dist)/2)
end


text:
using NeXLMatrixCorrection

struct DBKRatio
    database::SQLite.DB
    pkey::Int
    campaign::Int
    spectrum::Int
    primary::CharXRay
    lines::Vector{CharXRay}
    mode::String
    measured::Dict{Symbol,Any}
    reference::Dict{Symbol,Any}
    kratio::UncertainValue
end

function Base.show(io::IO, kr::DBKRatio)
    print( io, "$(krname(kr)): $(kr.lines)] = $(kr.kratio)")
end

function krname(kr::DBKRatio)
    nm = haskey(kr.measured, :Composition) ? name(kr.measured[:Composition]) : "Unspecified"
    return "K[($nm @ $(0.001*kr.measured[:BeamEnergy]) keV)/($(name(kr.reference[:Composition])) @ $(0.001*kr.reference[:BeamEnergy]) keV)"
end

function Base.read(db::SQLite.DB, ::Type{DBKRatio}, pkey::Int)::DBKRatio
    stmt = SQLite.Stmt(db, "SELECT * FROM KRATIO WHERE PKEY=?;")
    q = DBInterface.execute(stmt, (pkey,))
    SQLite.done(q) && error("There is no k-ratio with pkey = $(pkey).")
    r = SQLite.Row(q)
    primary = CharXRay(r[:ELEMENT], Transition(SubShell(r[:INNER]), SubShell(r[:OUTER])))
    lines = map(s -> parse(CharXRay, s), split(r[:LINES], ","))
    meas = Dict{Symbol,Any}(:BeamEnergy => r[:MEASE0], :TakeOffAngle => r[:MEASTOA])
    if r[:MEASURED] != -1
        meas[:Composition] = read(db, Material, r[:MEASURED])
    end
    ref = Dict{Symbol,Any}(:BeamEnergy => r[:REFE0], :TakeOffAngle => r[:REFTOA], :Composition => read(db, Material, r[:REFERENCE]))
    kr = uv(r[:KRATIO], r[:DKRATIO])
    return DBKRatio(db, pkey, r[:CAMPAIGN], r[:SPECPKEY], primary, lines, r[:MODE], meas, ref, kr)
end

function Base.findall(db::SQLite.DB, ::Type{DBKRatio}; material::Union{String,Nothing}=nothing, campaign::Union{Int,Nothing}=nothing, elm::Union{Element,Nothing}=nothing, mink::Float64=0.1)::Vector{DBKRatio}
    bs, args = "KRATIO >= ?", [ mink, ]
    if !isnothing(material)
        matkey = find(db, Material, material)
        @assert matkey!=-1 "Unable to find the material `$material`."
        bs *= " AND MEASURED = ?"
        push!(args, matkey)
    end
    if !isnothing(campaign)
        bs *= " AND CAMPAIGN = ?"
        push!(args, campaign)
    end
    if !isnothing(elm)
        bs *= " AND ELEMENT=?"
        push!(args, z(elm))
    end
    stmt = SQLite.Stmt(db, "SELECT PKEY FROM KRATIO WHERE $bs;")
    q = DBInterface.execute(stmt, args)
    return SQLite.done(q) ? [] : [read(db, DBKRatio, r[:PKEY]) for r in q]
end

function NeXLUncertainties.asa(
    ::Type{DataFrame}, #
    krs::AbstractVector{DBKRatio}; #
    withComputedKs::Bool = false, #
    mc::Type{<:MatrixCorrection} = XPP,
    fc::Type{<:FluorescenceCorrection} = ReedFluorescence,
    cc::Type{<:CoatingCorrection} = Coating
)
    fm, sm, cm = Union{Float64,Missing}, Union{String,Missing}, Union{Material,Missing}
    fs, lines, mease0, meastoa, meascomp = Int[], Vector{CharXRay}[], Float64[], Float64[], cm[]
    refe0, reftoa, refcomp, krv, dkrv, cks, ratio = Float64[], Float64[], cm[], Float64[], Float64[], fm[], fm[]
    for kr in krs
        push!(fs, kr.campaign)
        push!(lines, kr.lines)
        meas = kr.measured
        push!(mease0, get(meas, :BeamEnergy, missing))
        push!(meastoa, get(meas, :TakeOffAngle, missing))
        push!(meascomp, get(meas, :Composition, missing))
        ref = kr.reference
        push!(refe0, get(ref, :BeamEnergy, missing))
        push!(reftoa, get(ref, :TakeOffAngle, missing))
        push!(refcomp, get(ref, :Composition, missing))
        push!(krv, value(kr.kratio))
        push!(dkrv, σ(kr.kratio))
        if withComputedKs
            elm = element(kr.lines[1])
            if any(ismissing.((meascomp[end], mease0[end], meastoa[end], refcomp[end], refe0[end], reftoa[end]))) ||
               (NeXLCore.nonneg(meascomp[end], elm) < 1.0e-6) ||
               (NeXLCore.nonneg(refcomp[end], elm) < 1.0e-6) || (energy(inner(kr.primary))>0.95*min(mease0[end],refe0[end]))
                push!(cks, missing)
                push!(ratio, missing)
            else
                br = [ kr.primary ]
                zs = zafcorrection(mc, fc, cc, meascomp[end], br, mease0[end])
                zr = zafcorrection(mc, fc, cc, refcomp[end], br, refe0[end])
                k =
                    gZAFc(zs, zr, meastoa[end], reftoa[end] ) * NeXLCore.nonneg(meascomp[end], elm) /
                    NeXLCore.nonneg(refcomp[end], elm)
                push!(cks, k)
                push!(ratio, value(kr.kratio) / k)
            end
        end
    end
    res = DataFrame(
        Campaign = fs,
        Lines = lines,
        Cmeas = meascomp,
        E0meas = mease0,
        TOAmeas = meastoa,
        Cref = refcomp,
        E0ref = refe0,
        TOAref = reftoa,
        K = krv,
        ΔK = dkrv,
    )
    if withComputedKs
        res[:, :Kxpp] = cks
        res[:, :Ratio] = ratio
    end
    return res
end

NeXLUncertainties.asa(::Type{KRatio}, dbkr::DBKRatio)::KRatio =
    KRatio(dbkr.lines, dbkr.measured, dbkr.reference, dbkr.reference[:Composition], dbkr.kratio)

function Base.write(
    db::SQLite.DB,
    ::Type{DBKRatio},
    campaign::DBCampaign,
    meas::Spectrum,
    meascomp::Union{Material,Missing},
    res::FilterFitResult,
)
    stmt1 = SQLite.Stmt(
        db,
        "INSERT INTO KRATIO(CAMPAIGN, SPECPKEY, ELEMENT, INNER, OUTER, MODE, MEASURED, MEASNAME, MEASE0, MEASTOA, " *
        "REFERENCE, REFNAME, REFE0, REFTOA, PRINCIPAL, LINES, KRATIO, DKRATIO) VALUES (" *
        "?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?, ?);",
    )
    meascompidx = ismissing(meascomp) ? -1 : write(db, meascomp)
    measname = ismissing(meascomp) ? "Unknown" : name(meascomp)
    for lbl in filter(l -> (l isa CharXRayLabel) && (value(res[l]) > 0.0), labels(res))
        # @show lbl, value(res[lbl])
        ref, br = spectrum(lbl), brightest(lbl.xrays)
        refcomp = ref[:Composition]
        refcompidx = write(db, refcomp)
        r = DBInterface.execute(
            stmt1,
            (
                campaign.pkey, #
                meas[:PKEY],
                z(element(br)),
                inner(br).subshell.index,
                outer(br).subshell.index,
                "EDX", #
                meascompidx,
                measname,
                meas[:BeamEnergy],
                meas[:TakeOffAngle], #
                refcompidx,
                refcomp.name,
                ref[:BeamEnergy],
                ref[:TakeOffAngle], #
                repr(br),
                join(repr.(lbl.xrays), ","), #
                value(res[lbl]),
                uncertainty(res[lbl]),
            ),
        )
    end
end

"""
    NeXLMatrixCorrection.quantify(#
        krs::AbstractVector{DBKRatio};
        strip::AbstractVector{Element} = Element[],
        mc::Type{<:MatrixCorrection} = XPP,
        fc::Type{<:FluorescenceCorrection} = ReedFluorescence,
        cc::Type{<:CoatingCorrection} = Coating,
        kro::KRatioOptimizer = SimpleKRatioOptimizer(1.5),
        unmeasured::UnmeasuredElementRule = NullUnmeasuredRule(),
    )::Vector{IterationResult}

Quantify a collection of `DBKRatio`.
"""
function NeXLMatrixCorrection.quantify(#
    krs::AbstractVector{DBKRatio};
    strip::AbstractVector{Element} = Element[],
    mc::Type{<:MatrixCorrection} = XPP,
    fc::Type{<:FluorescenceCorrection} = ReedFluorescence,
    cc::Type{<:CoatingCorrection} = Coating,
    kro::KRatioOptimizer = SimpleKRatioOptimizer(1.5),
    unmeasured::UnmeasuredElementRule = NullUnmeasuredRule(),
)::Vector{IterationResult}
    iter = Iteration(mc, fc, cc, unmeasured = unmeasured)
    map(unique(kr.spectrum for kr in krs)) do spec
        skrs = asa.(KRatio, filter(kr->kr.spectrum==spec, krs))
        okrs = optimizeks(kro, filter(kr -> !(element(kr) in strip), skrs))
        quantify(iter, label("Unknown[$spec]"), okrs)
    end
end


text:
radius(::Hydrogen) = 0.2
radius(::Carbon) = 0.3
radius(::Nitrogen) = 8
radius(::Oxygen) = 9

function writetoblenderscript(mol::Molecule, filename)
    lines = String[]
    push!(lines, "import bpy", "import numpy")

    for element in ("Hydrogen", "Carbon", "Nitrogen", "Oxygen")
        push!(lines, """$(lowercase(element))_material = bpy.data.materials.new(name='$(element)Material')""")
    end

    for atom in mol.atoms
        x, y, z = position(atom)
        push!(lines, """bpy.ops.mesh.primitive_uv_sphere_add(radius=$(radius(atom)), location=($x, $y, $z))""")
    end

    for bond in mol.bonds
        x, y, z = (position(bond.atom2) .+ position(bond.atom1)) ./ 2
        δ = position(bond.atom2) .- position(bond.atom1)
        d = norm(δ)
        #z += d / 2
        θ = π/2 - atan(δ[3], sqrt(δ[1]^2 + δ[2]^2))
        ϕ = atan(δ[2], δ[1])

        push!(lines, """bpy.ops.mesh.primitive_cylinder_add(depth=$(norm(δ)), radius=0.1, location=($x, $y, $z), rotation=(0.0, $θ, $ϕ))""")
    end

    script = join(lines, "\n")
    write(filename, script)
end


text:
using FactCheck

include("LabelMe.jl")
include("XML.jl")

FactCheck.exitstatus()


text:
export Intlist
export is_empty, top, next, contains, push!, pop!, pop_all!, clear!

import Base.push!, Base.pop!, Base.contains

const NIL = Int64(-1)
const EMPTY = Int64(0)

mutable struct Intlist
    capacity::Int64
    size::Int64

    values::Vector{Int64}
    ip::Int64

    Intlist(capacity::Int64) = new(capacity, 0, fill(EMPTY, capacity), NIL)
end

function is_empty(il::Intlist)::Bool
    return 0 == il.size
end

function top(il::Intlist)::Int64
    return il.ip
end

@inbounds function next(il::Intlist, key::Int64)
    if NIL == key
        return key
    else
        return il.values[key]
    end
end

@inbounds function Base.contains(il::Intlist, key::Int64)::Bool
    return key > 0 && key <= il.capacity && EMPTY != il.values[key]
end

@inbounds function Base.:push!(il::Intlist, key::Int64)
    @assert key > 0 && key <= il.capacity

    if contains(il, key)
        return
    end

    il.values[key] = il.ip
    il.ip = key

    il.size += 1

    return
end

@inbounds function Base.:pop!(il::Intlist)::Int64
    @assert il.size > 0

    local key = il.ip
    il.ip = il.values[il.ip]
    il.values[key] = EMPTY
    il.size -= 1

    return key
end

@inbounds function pop_all!(il::Intlist)::Vector{Int64}
    local n = il.size
    local values = zeros(Int64, n)

    for i = 1:n
        values[i] = pop!(il)
    end

    return values
end

@inbounds function clear!(il::Intlist)
    while il.size > 0
        pop!(il)
    end

    return
end


text:
using Literate
using Printf

repo_src = joinpath(@__DIR__,"..","src")
notebooks_dir = joinpath(@__DIR__,"..","notebooks")

files = [
  "Poisson equation"=>"poisson.jl",
  "Code validation"=>"validation.jl",
  "Linear elasticity"=>"elasticity.jl",
  "p-Laplacian"=>"p_laplacian.jl",
  "Hyper-elasticity"=>"hyperelasticity.jl",
  "Poisson equation (with DG)"=> "dg_discretization.jl",
  "Darcy equation (with RT)"=>"darcy.jl",
  "Incompressible Navier-Stokes"=>"inc_navier_stokes.jl",
  "Stokes equation" => "stokes.jl",
  "Isotropic damage model" => "isotropic_damage.jl",
  "Fluid-Structure Interaction"=>"fsi_tutorial.jl",
  "Low-level API Poisson equation"=>"poisson_dev_fe.jl"]

Sys.rm(notebooks_dir;recursive=true,force=true)
for (i,(title,filename)) in enumerate(files)
  notebook_prefix = string("t",@sprintf "%03d_" i)
  notebook = string(notebook_prefix,splitext(filename)[1])
  notebook_title = string("# # Tutorial ", i, ": ", title)
  function preprocess_notebook(content)
    return string(notebook_title, "\n\n", content)
  end
  Literate.notebook(joinpath(repo_src,filename), notebooks_dir; name=notebook, preprocess=preprocess_notebook, documenter=false, execute=false)
end

deps_jl = joinpath(@__DIR__, "deps.jl")

if isfile(deps_jl)
  rm(deps_jl)
end

open(deps_jl,"w") do f
  println(f, "# This file is automatically generated")
  println(f, "# Do not edit")
  println(f)
  println(f, :(const files = $files))
end

max_stars_repo_name: WenjieYao/Tutorials

max_issues_repo_name: WenjieYao/Tutorials

max_forks_repo_name: WenjieYao/Tutorials


text:
function process(group::Group; config = Dict())
    for child in group.children
        #child.config[:active] && process(child, config = config)
        process(child, config = config)
    end
    group
end

function process(doc::Document; config = Dict())
    !doc.config[:active] && return

    abs_src = joinpath(doc.root, doc.rel_path)
    config = deepcopy(config)

    haskey(config, "preprocess") && @warn "Overriding preprocess function in config"

    config = Literate.create_configuration(abs_src, user_config = config, user_kwargs = ())

    pre = x -> preprocess(x, doc, config = config)
    config["preprocess"] = pre

    builds = doc.config[:builds]
    paths = map(p -> joinpath(STAGING_DIR, p), PATHS)

    #copy_deps(doc) # TODO not needed with File type

    if doc.kind === :documenter && :markdown in builds
        abs_dst = joinpath(STAGING_DIR, doc.rel_path)
        mkpath(dirname(abs_dst))
        content = pre(read(abs_src, String))
        open(io -> write(io, content), abs_dst, "w")
    else
        if :markdown in builds
            abs_dst = joinpath(STAGING_DIR, doc.rel_path)
            mkpath(dirname(abs_dst))
            Literate.markdown(abs_src, dirname(abs_dst), config = config)
        end
        if :script in builds
            abs_dst = joinpath(paths.script, doc.rel_path)
            mkpath(dirname(abs_dst))
            Literate.script(abs_src, dirname(abs_dst), config = config)
        end
        if :notebook in builds
            abs_dst = joinpath(paths.notebook, doc.rel_path)
            mkpath(dirname(abs_dst))
            Literate.notebook(abs_src, dirname(abs_dst), config = config)
        end
    end
end

function process(file::File; config = Dict())
    abs_dst = joinpath(STAGING_DIR, file.rel_path)
    mkpath(dirname(abs_dst))
    cp(joinpath(file.root, file.rel_path), abs_dst)
end

function copy_deps(doc::Document)
    if !isempty(doc.config[:deps])
        for rel_path in doc.config[:deps]
            rel_path = joinpath(dirname(doc.rel_path), rel_path) # relative to doc.root
            abs_src = joinpath(doc.root, rel_path)
            abs_dst = joinpath(STAGING_DIR, rel_path)
            mkpath(dirname(abs_dst))
            cp(abs_src, abs_dst)
        end
    end
    nothing
end

function build_pages(group::Group)
    pages = filter(x -> !isnothing(x), map(_build_pages, group.children))
    sort!(pages, by = p -> p.first)
    pages = map(p -> p.second, pages)
end

function _build_pages(group::Group)
    pages = []
    for child in group.children
        page = _build_pages(child)
        page !== nothing && push!(pages, page)
    end
    sort!(pages, by = p -> p.first)
    pages = map(p -> p.second, pages)
    group.config[:weight] => group.config[:short_title] => pages
end

function _build_pages(doc::Document)
    if doc.config[:active] && (:markdown in doc.config[:builds])
        path = first(splitext(doc.rel_path)) * ".md"
        page = doc.config[:short_title] => path
        page = doc.config[:hide] ? hide(page) : page
        return doc.config[:weight] => page
    else
        return nothing
    end
end

_build_pages(::File) = nothing

function indented_println(xs...; indent = 0)
    for _ = 1:(Base.indent_width*indent)
        print(' ')
    end
    println(xs...)
end

function print_pages(index, indent = 0)
    for (title, page_or_section) in index
        if page_or_section isa String
            indented_println(title, " => ", page_or_section, indent = indent)
        else
            indented_println(title, indent = indent)
            print_pages(page_or_section, indent + 1)
        end
    end
end

function execute_block(block::String)
    m = Module(gensym())
    # eval(expr) is available in the REPL (i.e. Main) so we emulate that for the sandbox
    Core.eval(m, :(eval(x) = Core.eval($m, x)))
    # modules created with Module() does not have include defined
    # abspath is needed since this will call `include_relative`
    Core.eval(m, :(include(x) = Base.include($m, abspath(x))))

    # r is the result
    # status = (true|false)
    # _: backtrace
    # str combined stdout, stderr output
    r, status, _, str = Literate.Documenter.withoutput() do
        include_string(m, block)
    end
    if !status
        error("""
             $(sprint(showerror, r))
             when executing the following code block:

             $block
             """)
    end
    m
end

islocalbuild() = get(ENV, "GITHUB_ACTIONS", nothing) != "true"


text:
using DelimitedFiles

function read_vector(doc, n)
    vector=[]
    data = readdlm(doc, skipstart=n)
    for i in data[1,:]
        append!(vector, i)
    end
    return(n+1, vector)
end

function read_matrix(P, Q, doc, n)
    matrix = Matrix{Int8}(undef, P, Q)
    data = readdlm(doc, skipstart=n)
    for i = 1:P
        for j = 1:Q
            matrix[i,j] = data[i,j]
        end
    end
    return(n+P, matrix)
end

function read_data(doc)
    open(doc, "r") do f
        n = 11
        C = parse(Int, readline(f)) #number of containers
        P = parse(Int, readline(f)) #number of positions
        CP = parse(Int, readline(f)) #number of compatible combinations
        S = parse(Int, readline(f)) #number of transport vehicles
        Q = parse(Int, readline(f)) #number of quay cranes
        J = parse(Int, readline(f)) #number of bays
        tt = parse(Int, readline(f)) #quay crane travel time per bay
        d = parse(Int, readline(f)) #clearance, in number of bays
        mt = parse(Int, readline(f)) #maximum time
        beta = parse(Int, readline(f)) #beta
        EFT = parse(Int, readline(f)) #Expected Finishing Time

        n, ci = read_vector(doc, n) #container class vector
        n, pj = read_vector(doc, n) #position class vector

        n, qj = read_matrix(P,Q,doc,n) #quay cranes position matrix
        n, sj = read_matrix(P,S,doc,n) #transfer vehicle position matrix
        n, prejj = read_matrix(P,P,doc,n) #position precedence matrix
        n, cpij = read_matrix(CP, 3, doc,n) #combination c, p, time(cp)
        return(C,P,CP,S,Q,J,tt,d,mt,beta,EFT,ci,pj,qj,sj,prejj,cpij)
    end
end

doc = "./data/60C_10Type_LessDense_2QC_6TV_mod.txt"
C,P,CP,S,Q,J,tt,d,mt,beta,EFT,ci,pj,qj,sj,prejj,cpij = read_data(doc)

max_stars_repo_licenses: Zlib MIT

max_issues_repo_licenses: Zlib MIT

max_forks_repo_licenses: Zlib MIT


text:
# this file tests the rk4 function
# this uses rk4 to integrate the function t^4 + t^3 + t^2 +t + 1
# thus the function supplied to rk4 is d/dt of that function
# and q_vec_rk[1] is initially function(t=0) = 1
function true_f(t)
#   return t + 1
  return t^4 + t^3 + t^2 + t + 1
end

function test_f(q_vec, res_vec, opts, t)
#   return res_vec[1] = 1
  res_vec[1] = 4*t^3 + 3*t^2 + 2*t + 1
end

pre_func = (a...) -> 3
function post_func(q_vec, res_vec, opts; calc_norm=true)
  return norm(res_vec)
end

"""
  Use rk4 to solve an 4th order polynoial ODE, for which it is exact
"""
function test_rk4()
    @testset "----- testing rk4 -----" begin
    opts_rk = Dict{Any, Any}("physics" => "Euler", "run_type" => 1, "smb_name" => "", "output_freq" => 1, "write_vis" => false, "use_itermax" => false, "numBC" => 0)
    read_input(opts_rk)

    q_vec_rk = [1.0]
    res_vec_rk = [0.0]
    ctx = (q_vec_rk, res_vec_rk)

    delta_t = 0.1
    t_max = 1.0

    t = rk4(test_f, delta_t, t_max, q_vec_rk, res_vec_rk, pre_func, post_func, ctx, opts_rk, real_time=true)

    f_approx = q_vec_rk[1]
    f_exact = true_f(t)
    err = abs(f_approx - f_exact)

    @test isapprox( t, t_max) atol=1e-14
    @test  err  < 1e-14# machine precision
  end  # end facts block

  @testset "----- testing lserk54 -----" begin
    opts_rk = Dict{Any, Any}("physics" => "Euler", "run_type" => 30, "smb_name" => "", "output_freq" => 1, "write_vis" => false, "use_itermax" => false, "numBC" => 0)
    read_input(opts_rk)

    q_vec_rk = [1.0]
    res_vec_rk = [0.0]
    ctx = (q_vec_rk, res_vec_rk)

    delta_t = 0.1
    t_max = 1.0

    t = lserk54(test_f, delta_t, t_max, q_vec_rk, res_vec_rk, pre_func, post_func, ctx, opts_rk, real_time=true)

    f_approx = q_vec_rk[1]
    f_exact = true_f(t)
    err = abs(f_approx - f_exact)

    @test  err  < 1e-14# machine precision

  end

  return nothing
end

#test_rk4()
add_func1!(EulerTests, test_rk4, [TAG_NLSOLVERS, TAG_SHORTTEST])


text:
"""
    rowrange(n, m, i)

Utility function that returns a range based on an iteration index `i`,
the number of elements `n` and number of ranges `m`.
"""
@inline rowrange(n::Int, m::Int, i::Int) = begin
    rs = floor(Int, n/m)  # row step
    rr = rs*(i-1)+1 : rs*i    # row range
    return rr
end

"""
    quantized_eltype(k)

Determinies the minimum `Unsigned` type that can hold the value `k`.
"""
quantized_eltype(k) = begin
    b = clamp(log2(k), 1, MAX_BITS)
    minbits = 1024
    local mintype
    for (nbits, utype) in BITS_TO_TYPE
        if b <= nbits && nbits < minbits
            minbits = nbits
            mintype = utype
        end
    end
    if minbits > MAX_BITS
        @error "Number is too large to fit in $MAX_BITS bits."
    end
    return mintype
end

"""
    quantization_type(method)

Returns a quantization type based on the quantization method to be employed.
"""
quantization_type(method::Symbol) = begin
    ifelse(method == :rvq, AdditiveQuantization(), OrthogonalQuantization())
end


text:
module Zlib

using Compat

import Base: read, read!, readuntil, readbytes!, write, close, eof

export compress, decompress, crc32

const Z_NO_FLUSH      = 0
const Z_PARTIAL_FLUSH = 1
const Z_SYNC_FLUSH    = 2
const Z_FULL_FLUSH    = 3
const Z_FINISH        = 4
const Z_BLOCK         = 5
const Z_TREES         = 6

const Z_OK            = 0
const Z_STREAM_END    = 1
const Z_NEED_DICT     = 2
const ZERRNO          = -1
const Z_STREAM_ERROR  = -2
const Z_DATA_ERROR    = -3
const Z_MEM_ERROR     = -4
const Z_BUF_ERROR     = -5
const Z_VERSION_ERROR = -6

@unix_only const libz = "libz"
@windows_only const libz = "zlib1"

# The zlib z_stream structure.
type z_stream
    next_in::Ptr{UInt8}
    avail_in::Cuint
    total_in::Culong

    next_out::Ptr{UInt8}
    avail_out::Cuint
    total_out::Culong

    msg::Ptr{UInt8}
    state::Ptr{Void}

    zalloc::Ptr{Void}
    zfree::Ptr{Void}
    opaque::Ptr{Void}

    data_type::Cint
    adler::Culong
    reserved::Culong

    function z_stream()
        strm = new()
        strm.next_in   = C_NULL
        strm.avail_in  = 0
        strm.total_in  = 0
        strm.next_out  = C_NULL
        strm.avail_out = 0
        strm.total_out = 0
        strm.msg       = C_NULL
        strm.state     = C_NULL
        strm.zalloc    = C_NULL
        strm.zfree     = C_NULL
        strm.opaque    = C_NULL
        strm.data_type = 0
        strm.adler     = 0
        strm.reserved  = 0
        strm
    end
end

type gz_header
    text::Cint          # true if compressed data believed to be text
    time::Culong        # modification time
    xflags::Cint        # extra flags (not used when writing a gzip file)
    os::Cint            # operating system
    extra::Ptr{UInt8}   # pointer to extra field or Z_NULL if none
    extra_len::Cuint    # extra field length (valid if extra != Z_NULL)
    extra_max::Cuint    # space at extra (only when reading header)
    name::Ptr{UInt8}    # pointer to zero-terminated file name or Z_NULL
    name_max::Cuint     # space at name (only when reading header)
    comment::Ptr{UInt8} # pointer to zero-terminated comment or Z_NULL
    comm_max::Cuint     # space at comment (only when reading header)
    hcrc::Cint          # true if there was or will be a header crc
    done::Cint          # true when done reading gzip header (not used
                        # when writing a gzip file)
    gz_header() = new(0,0,0,0,0,0,0,0,0,0,0,0,0)
end

function zlib_version()
    ccall((:zlibVersion, libz), Ptr{UInt8}, ())
end

type Writer <: IO
    strm::z_stream
    io::IO
    closed::Bool

    Writer(strm::z_stream, io::IO, closed::Bool) =
        (w = new(strm, io, closed); finalizer(w, close); w)
end

function Writer(io::IO, level::Integer, gzip::Bool=false, raw::Bool=false)
    if !(1 <= level <= 9)
        error("Invalid zlib compression level.")
    end

    strm = z_stream()
    ret = ccall((:deflateInit2_, libz),
                Int32, (Ptr{z_stream}, Cint, Cint, Cint, Cint, Cint, Ptr{UInt8}, Int32),
                &strm, level, 8, raw? -15 : 15+gzip*16, 8, 0, zlib_version(), sizeof(z_stream))

    if ret != Z_OK
        error("Error initializing zlib deflate stream.")
    end

    if gzip && false
        hdr = gz_header()
        ret = ccall((:deflateSetHeader, libz),
            Cint, (Ptr{z_stream}, Ptr{gz_header}),
            &strm, &hdr)
        if ret != Z_OK
            error("Error setting gzip stream header.")
        end
    end

    Writer(strm, io, false)
end

Writer(io::IO, gzip::Bool=false, raw::Bool=false) = Writer(io, 9, gzip, raw)

function write(w::Writer, p::Ptr, nb::Integer)
    w.strm.next_in = p
    w.strm.avail_in = nb
    outbuf = Array(UInt8, 1024)

    while true
        w.strm.avail_out = length(outbuf)
        w.strm.next_out = pointer(outbuf)

        ret = ccall((:deflate, libz),
                    Int32, (Ptr{z_stream}, Int32),
                    &w.strm, Z_NO_FLUSH)
        if ret != Z_OK
            error("Error in zlib deflate stream ($(ret)).")
        end

        n = length(outbuf) - w.strm.avail_out
        if n > 0 && write(w.io, outbuf[1:n]) != n
            error("short write")
        end
        if w.strm.avail_out != 0
            break
        end
    end
    nb
end

# If this is not provided, Base.IO write methods will write
# arrays one element at a time.
function write{T}(w::Writer, a::Array{T})
    if isbits(T)
        write(w, pointer(a), length(a)*sizeof(T))
    else
        invoke(write, (IO, Array), w, a)
    end
end

# Copied from Julia base/io.jl
function write{T,N,A<:Array}(w::Writer, a::SubArray{T,N,A})
    if !isbits(T) || stride(a,1)!=1
        return invoke(write, (Any, AbstractArray), s, a)
    end
    colsz = size(a,1)*sizeof(T)
    if N<=1
        return write(s, pointer(a, 1), colsz)
    else
        # WARNING: cartesianmap(f,dims) is deprecated, use for idx = CartesianRange(dims)
        #     f(idx.I...)

        if VERSION >= v"0.4.0-"
            for idx in CartesianRange(tuple(1, size(a)[2:end]...))
                write(w, pointer(a, idx.I), colsz)
            end
        else
            cartesianmap((idxs...)->write(w, pointer(a, idxs), colsz),
                        tuple(1, size(a)[2:end]...))
        end

        return colsz*Base.trailingsize(a,2)
    end
end

function write(w::Writer, b::UInt8)
    write(w, UInt8[b])
end

function close(w::Writer)
    if w.closed
        return
    end
    w.closed = true

    # flush zlib buffer using Z_FINISH
    inbuf = Array(UInt8, 0)
    w.strm.next_in = pointer(inbuf)
    w.strm.avail_in = 0
    outbuf = Array(UInt8, 1024)
    ret = Z_OK
    while ret != Z_STREAM_END
        w.strm.avail_out = length(outbuf)
        w.strm.next_out = pointer(outbuf)
        ret = ccall((:deflate, libz),
                    Int32, (Ptr{z_stream}, Int32),
                    &w.strm, Z_FINISH)
        if ret != Z_OK && ret != Z_STREAM_END
            error("Error in zlib deflate stream ($(ret)).")
        end
        n = length(outbuf) - w.strm.avail_out
        if n > 0 && write(w.io, outbuf[1:n]) != n
            error("short write")
        end
    end

    ret = ccall((:deflateEnd, libz), Int32, (Ptr{z_stream},), &w.strm)
    if ret == Z_STREAM_ERROR
        error("Error: zlib deflate stream was prematurely freed.")
    end
end

function compress(input::Vector{UInt8}, level::Integer, gzip::Bool=false, raw::Bool=false)
    b = IOBuffer()
    w = Writer(b, level, gzip, raw)
    write(w, input)
    close(w)
    takebuf_array(b)
end

function compress(input::AbstractString, level::Integer, gzip::Bool=false, raw::Bool=false)
    compress(convert(Vector{UInt8}, input), level, gzip, raw)
end

compress(input::Vector{UInt8}, gzip::Bool=false, raw::Bool=false) = compress(input, 9, gzip, raw)
compress(input::AbstractString, gzip::Bool=false, raw::Bool=false) = compress(input, 9, gzip, raw)

type Reader <: IO
    strm::z_stream
    io::IO
    buf::IOBuffer
    closed::Bool
    bufsize::Int
    stream_end::Bool

    Reader(strm::z_stream, io::IO, buf::IOBuffer, closed::Bool, bufsize::Int) =
        (r = new(strm, io, buf, closed, bufsize, false); finalizer(r, close); r)
end

function Reader(io::IO, raw::Bool=false; bufsize::Int=4096)
    strm = z_stream()
    ret = ccall((:inflateInit2_, libz),
                Int32, (Ptr{z_stream}, Cint, Ptr{UInt8}, Int32),
                &strm, raw? -15 : 47, zlib_version(), sizeof(z_stream))
    if ret != Z_OK
        error("Error initializing zlib inflate stream.")
    end

    Reader(strm, io, PipeBuffer(), false, bufsize)
end

# Fill up the buffer with at least minlen bytes of uncompressed data,
# unless we have already reached EOF.
function fillbuf(r::Reader, minlen::Integer)
    ret = Z_OK
    while nb_available(r.buf) < minlen && !eof(r.io) && ret != Z_STREAM_END
        input = read(r.io, UInt8, min(nb_available(r.io), r.bufsize))
        r.strm.next_in = pointer(input)
        r.strm.avail_in = length(input)
        #outbuf = Array(UInt8, r.bufsize)

        while true
            #r.strm.next_out = outbuf
            #r.strm.avail_out = length(outbuf)
            (r.strm.next_out, r.strm.avail_out) = Base.alloc_request(r.buf, convert(UInt, r.bufsize))
            actual_bufsize_out = r.strm.avail_out
            ret = ccall((:inflate, libz),
                        Int32, (Ptr{z_stream}, Int32),
                        &r.strm, Z_NO_FLUSH)
            if ret == Z_DATA_ERROR
                error("Error: input is not zlib compressed data: $(bytestring(r.strm.msg))")
            elseif ret != Z_OK && ret != Z_STREAM_END && ret != Z_BUF_ERROR
                error("Error in zlib inflate stream ($(ret)).")
            end
            if (nbytes = actual_bufsize_out - r.strm.avail_out) > 0
                #write(r.buf, pointer(outbuf), nbytes)
                # TODO: the last two parameters are not used by notify_filled()
                # and can be removed if Julia PR #4484 is merged
                Base.notify_filled(r.buf, convert(Int, nbytes), C_NULL, convert(UInt, 0))
            end
            if r.strm.avail_out != 0
                break
            end
        end
    end

    if ret == Z_STREAM_END
        r.stream_end = true
    end

    nb_available(r.buf)
end

# This is to fix the ambiguity with Base.read!
function read!(r::Reader, a::Vector{UInt8})
    nb = length(a)
    if fillbuf(r, nb) < nb
        throw(EOFError())
    end
    read!(r.buf, a)
    a
end

function read!{T}(r::Reader, a::Array{T})
    if isbits(T)
        nb = length(a)*sizeof(T)
        if fillbuf(r, nb) < nb
            throw(EOFError())
        end
        read!(r.buf, a)
    else
        invoke(read!, (IO, Array), r, a)
    end
    a
end

# This function needs to be fast because other read calls use it.
function read(r::Reader, ::Type{UInt8})
    if nb_available(r.buf) < 1 && fillbuf(r, 1) < 1
        throw(EOFError())
    end
    read(r.buf, UInt8)
end

# This is faster than using the generic implementation in Base. We use
# it (indirectly) for decompress below.
readbytes!(r::Reader, b::AbstractArray{UInt8}, nb=length(b)) =
    readbytes!(r.buf, b, fillbuf(r, nb))

function readuntil(r::Reader, delim::UInt8)
    nb = search(r.buf, delim)
    while nb == 0
        offset = nb_available(r.buf)
        fillbuf(r, offset+r.bufsize)
        if nb_available(r.buf) == nb
            break
        end
        # TODO: add offset here when https://github.com/JuliaLang/julia/pull/4485
        # is merged
        nb = search(r.buf, delim) #, offset)
    end
    if nb == 0;  nb == nb_available(r.buf); end
    read!(r.buf, Array(UInt8, nb))
end

function close(r::Reader)
    if r.closed
        return
    end
    r.closed = true

    ret = ccall((:inflateEnd, libz), Int32, (Ptr{z_stream},), &r.strm)
    if ret == Z_STREAM_ERROR
        error("Error: zlib inflate stream was prematurely freed.")
    end
end

function eof(r::Reader)
    # Detecting EOF is somewhat tricky: we might not have reached
    # EOF in r.io but decompressing the remaining data might
    # yield no uncompressed data. So, make sure we can get at least
    # one more byte of decompressed data before we say we haven't
    # reached EOF yet.
    nb_available(r.buf) == 0 && eof(r.io)
end

function decompress(input::Vector{UInt8}, raw::Bool=false)
    r = Reader(IOBuffer(input), raw)
    b = readbytes(r)
    if !r.stream_end
        error("Error: zlib compressed data is incomplete or truncated")
    end
    close(r)
    b
end

decompress(input::AbstractString, raw::Bool=false) = decompress(convert(Vector{UInt8}, input), raw)

function crc32(data::Vector{UInt8}, crc::Integer=0)
    convert(UInt32, (ccall((:crc32, libz),
                 Culong, (Culong, Ptr{UInt8}, Cuint),
                 crc, data, length(data))))
end

crc32(data::AbstractString, crc::Integer=0) = crc32(convert(Vector{UInt8}, data), crc)

end # module


text:
include("common/constructor_validations.jl")
include("common/device_constructor_utils.jl")
include("thermalgeneration_constructor.jl")
include("hydrogeneration_constructor.jl")
include("branch_constructor.jl")
include("renewablegeneration_constructor.jl")
include("load_constructor.jl")
include("storage_constructor.jl")


text:
using Logging, Test, HSM

struct RootStateTestState <: HSM.AbstractHsmState
    state_info::HSM.HsmStateInfo

    RootStateTestState(parent_state) = new(HSM.HsmStateInfo(parent_state))
end

"""
    basic_state_machine()

Builds and returns the following state machine:

```plantuml
@startuml
state A {
    state B
    state C {
        state D
        state E
    }
}
@enduml
```
"""
function root_state_state_machine()
    A = RootStateTestState(nothing)
    B = RootStateTestState(A)
    C = RootStateTestState(A)
    D = RootStateTestState(C)
    E = RootStateTestState(C)

    return A, B, C, D, E
end

@testset "root_state()" begin
    A, B, C, D, E = root_state_state_machine()

    @test HSM.root_state(A) == A
    @test HSM.root_state(B) == A
    @test HSM.root_state(C) == A
    @test HSM.root_state(D) == A
    @test HSM.root_state(E) == A
end


text:
__precompile__()

module HeapedVectors

import Base: getindex, length, push!, isempty,
        pop!, filter!, popfirst!

export HeapedVector

import ..StrategyBase:filter_elements!
using ..StrategyBase

struct HeapedVector{T, F<:Function} <: Strategy
    data::Vector{T}
    by::F
    function HeapedVector(v::Vector{T}, by::F) where {T, F}
        new{T, F}(heaping(v, by), by)
    end
end

HeapedVector(data::Vector{T}) where {T} = HeapedVector(data, identity)

function heaping(v, by)
    ar = typeof(v[1])[]
    for i = 1:length(v)
        insert!(ar, i, v[i])
        floatup!(ar, length(ar), by)
    end
    return ar
end

function floatup!(ar, index, by)
    par = convert(Int, floor(index/2))
    if index <= 1
        return ar
    end
    if by(ar[index]) < by(ar[par])
       ar[par], ar[index] = ar[index], ar[par]
    end
    floatup!(ar, par, by)
end

function push!(v::HeapedVector{T}, x::T) where {T}
    insert!(v.data, length(v.data)+1, x)
    floatup!(v.data, length(v.data), v.by)
    return v
end

isempty(v::HeapedVector) = isempty(v.data)

function popfirst!(v::HeapedVector{T}) where {T}
    if length(v.data) == 0
         return
    end

    if length(v.data) > 2
        v.data[length(v.data)], v.data[1] = v.data[1], v.data[length(v.data)]
    	minm = pop!(v.data)
    	bubbledown!(v::HeapedVector{T}, 1)

    elseif length(v.data) == 2
        v.data[length(v.data)], v.data[1] = v.data[1], v.data[length(v.data)]
    	minm = pop!(v.data)
    else
    	minm = pop!(v.data)
    end
    return minm
end

function bubbledown!(v::HeapedVector{T}, index) where{T}
    left = index*2
    right = index*2+1
    smallest = index

    if length(v.data)+1 > left && v.by(v.data[smallest]) > v.by(v.data[left])
        smallest = left
    end

    if length(v.data)+1 > right && v.by(v.data[smallest]) > v.by(v.data[right])
	   smallest = right
    end

    if smallest != index
       v.data[index], v.data[smallest] = v.data[smallest], v.data[index]
	   bubbledown!(v, smallest)
    end
end

function filter_elements!(A::HeapedVector{T}, x::T) where{T}
    func(y) = A.by(y) <= A.by(x)
    filter!(func, A.data)

    if length(A.data) == 0
        return A
    end

    heaping(A.data, A.by)
    return A

end

end


text:
using Ipopt
using DelimitedFiles
using PyPlot
import ForwardDiff

using CSV
using DataFrames
# using YAML

# set up boundary constraint wrapper function
function boundary_wrapper(x, params)
    # include relevant globals
    params.boundary_vertices
    params.boundary_normals

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    # get and return boundary distances
    return ff.ray_trace_boundary(boundary_vertices, boundary_normals, turbine_x, turbine_y)
end

# set up spacing constraint wrapper function
function spacing_wrapper(x, params)
    # include relevant globals
    params.rotor_diameter

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    # get and return spacing distances
    return 2.0*rotor_diameter[1] .- ff.turbine_spacing(turbine_x,turbine_y)
end

# set up objective wrapper function
function aep_wrapper(x, params)
    # include relevant globals
    params.turbine_z
    params.rotor_diameter
    params.hub_height
    params.turbine_yaw
    params.ct_model
    params.generator_efficiency
    params.cut_in_speed
    params.cut_out_speed
    params.rated_speed
    params.rated_power
    params.windresource
    params.power_models
    params.model_set
    params.rotor_points_y
    params.rotor_points_z
    params.obj_scale

    # get number of turbines
    nturbines = Int(length(x)/2)

    # extract x and y locations of turbines from design variables vector
    turbine_x = x[1:nturbines]
    turbine_y = x[nturbines+1:end]

    # calculate AEP
    AEP = obj_scale*ff.calculate_aep(turbine_x, turbine_y, turbine_z, rotor_diameter,
                hub_height, turbine_yaw, ct_model, generator_efficiency, cut_in_speed,
                cut_out_speed, rated_speed, rated_power, windresource, power_models, model_set,
                rotor_sample_points_y=rotor_points_y,rotor_sample_points_z=rotor_points_z)

    # return the objective as an array
    return [AEP]
end

# objective function
function obj(x)
    global funcalls_AEP
    AEP = -aep_wrapper(x)[1]
    append!(funcalls_AEP, -AEP)
    return AEP
end

# constraint function
function con(x, g)
    # calculate spacing constraint value
    spacing_con = spacing_wrapper(x)

    # calculate boundary constraint
    boundary_con = boundary_wrapper(x)

    # combine spacing and boundary constraint values
    g[:] = [spacing_con; boundary_con]
end

# objective gradient function
function obj_grad(x, grad_f)
    grad_f[:] = -ForwardDiff.jacobian(aep_wrapper,x)
end

# constraint gradients function
function con_grad(x, mode, rows, cols, values)
    if mode == :Structure
        # report the sparcity structure of the jacobian
        for i = 1:prob.m
            rows[(i-1)*prob.n+1:i*prob.n] = i*ones(Int,prob.n)
            cols[(i-1)*prob.n+1:i*prob.n] = 1:prob.n
        end
        return rows, cols
    else
        # calculate spacing constraint jacobian
        ds_dx = ForwardDiff.jacobian(spacing_wrapper, x)

        # calculate boundary constraint jacobian
        db_dx = ForwardDiff.jacobian(boundary_wrapper, x)

        # combine constraint jacobians into an overall array of constraint derivatives
        for i = 1:prob.m
            for j = 1:prob.n
                values[(i-1)*prob.n+j] = [ds_dx; db_dx][i, j]
            end
        end
    end
    return values
end

# this function will be called at each iteration of the optimization
function intermediate(
    alg_mod::Int,
    iter_count::Int,
    obj_value::Float64,
    inf_pr::Float64, inf_du::Float64,
    mu::Float64, d_norm::Float64,
    regularization_size::Float64,
    alpha_du::Float64, alpha_pr::Float64,
    ls_trials::Int)
    global iter_AEP
    append!(iter_AEP, -obj_value)

    return true
end

# set globals for iteration history
iter_AEP = zeros(Float64, 0)
funcalls_AEP = zeros(Float64, 0)
global iter_AEP
global funcalls_AEP

# import model set with wind farm and related details
include("./model_sets/model_set_7_ieacs3.jl")

# scale objective to be between 0 and 1
obj_scale = 1E-12

# set wind farm boundary parameters
boundary_vertices = [10363.8 6490.3; 9449.7 1602.2; 9387.0 1056.6; 9365.1 625.5; 9360.8 360.2; 9361.5 126.9; 9361.3 137.1; 7997.6 1457.9; 6098.3 3297.5;
    8450.3 6455.3; 8505.4 6422.3; 9133.0 6127.4; 9332.8 6072.6; 9544.2 6087.1; 9739.0 6171.2; 9894.9 6316.9; 10071.8 6552.5; 10106.9 6611.1]
boundary_normals = [0.9829601758936983 -0.1838186405319916; 0.9934614633172962 -0.11416795042154541; 0.9987121579438882 -0.050734855622757584;
    0.9998686751666075 -0.01620593781838486; 0.9999954987444023 0.0030004151269687495; -0.9998078216567232 -0.019604074934516894; -0.6957179389375846 -0.718315076718037;
    -0.6957275377423737 -0.7183057797532565; -0.8019887481131871 0.5973391397688945; 0.5138086803485797 0.8579047965820281; 0.4252760929807897 0.905063668886888;
    0.2645057513093967 0.9643841078762402; -0.0684295708121141 0.9976559496331737; -0.39636379138742883 0.9180935381958544; -0.6828023205475376 0.7306031693435896;
    -0.7996740386176392 0.6004343694034798; -0.8578802011411015 0.5138497450520954; 0.42552559023380465 0.9049463918134445]

# set globals for use in wrapper functions
struct params_struct{}
    model_set
    rotor_points_y
    rotor_points_z
    turbine_z
    ambient_ti
    rotor_diameter
    boundary_vertices
    boundary_normals
    obj_scale
    hub_height
    turbine_yaw
    ct_model
    generator_efficiency
    cut_in_speed
    cut_out_speed
    rated_speed
    rated_power
    windresource
    power_models
end

params = params_struct(model_set, rotor_points_y, rotor_points_z, turbine_z, ambient_ti,
    rotor_diameter, boundary_vertices, boundary_normals, obj_scale, hub_height, turbine_yaw,
    ct_model, generator_efficiency, cut_in_speed, cut_out_speed, rated_speed, rated_power,
    windresource, power_models)

# initialize design variable array
x_initial = [copy(turbine_x);copy(turbine_y)]

# calculate initial state AEPs
state_aeps = ff.calculate_state_aeps(turbine_x, turbine_y, turbine_z, rotor_diameter,
                hub_height, turbine_yaw, ct_model, generator_efficiency, cut_in_speed,
                cut_out_speed, rated_speed, rated_power, windresource, power_models, model_set;
                rotor_sample_points_y=[0.0], rotor_sample_points_z=[0.0], hours_per_year=365.0*24.0)

dir_aep = zeros(20)
for i in 1:20
    for j in 1:20
        dir_aep[i] += state_aeps[(i-1)*20 + j]
    end
end

# add initial turbine location to plot
plt.clf()
for i = 1:length(turbine_x)
    plt.gcf().gca().add_artist(plt.Circle((turbine_x[i],turbine_y[i]), rotor_diameter[1]/2.0, fill=false,color="C0"))
end

# get number of design variables
n_designvariables = length(x_initial)

# get number of constraints
function numberofspacingconstraints(nturb)
    # calculates number of spacing constraints needed for given number of turbines
    ncon = 0
    for i = 1:nturb-1; ncon += i; end
    return ncon
end
n_spacingconstraints = numberofspacingconstraints(nturbines)
n_boundaryconstraints = length(boundary_wrapper(x_initial, params))
n_constraints = n_spacingconstraints + n_boundaryconstraints

# set general lower and upper bounds for design variables
lb = ones(n_designvariables) * -Inf
ub = ones(n_designvariables) * Inf

# set lower and upper bounds for constraints
lb_g = ones(n_constraints) * -Inf
ub_g = zeros(n_constraints)

# create the problem
prob = createProblem(n_designvariables, lb, ub, n_constraints, lb_g, ub_g, n_designvariables*n_constraints, 0,
    obj, con, obj_grad, con_grad)
addOption(prob, "hessian_approximation", "limited-memory")
prob.x = x_initial

# set function to be ran at every iteration
setIntermediateCallback(prob, intermediate)

# generate wrapper function surrogates
spacing_wrapper(x) = spacing_wrapper(x, params)
aep_wrapper(x) = aep_wrapper(x, params)
boundary_wrapper(x) = boundary_wrapper(x, params)

# run and time optimization
t1 = time()
status = solveProblem(prob)
t2 = time()
clkt = t2-t1
xopt = prob.x
fopt = prob.obj_val
info = Ipopt.ApplicationReturnStatus[status]

# print initial objective value
println("nturbines: ", nturbines)
println("rotor diameter: ", rotor_diameter[1])
println("starting AEP value (MWh): ", aep_wrapper(x_initial, params)[1]*1E5)

# print optimization results
println("Finished in : ", clkt, " (s)")
println("info: ", info)
println("end AEP value (MWh): ", aep_wrapper(xopt)[1]*1E5)

# extract final turbine locations
turbine_x = copy(xopt[1:nturbines])
turbine_y = copy(xopt[nturbines+1:end])

# add final turbine locations to plot
for i = 1:length(turbine_x)
    plt.gcf().gca().add_artist(plt.Circle((turbine_x[i],turbine_y[i]), rotor_diameter[1]/2.0, fill=false,color="C1", linestyle="--"))
end

# add wind farm boundary to plot
plt.gcf().gca().plot([boundary_vertices[:,1];boundary_vertices[1,1]],[boundary_vertices[:,2];boundary_vertices[1,2]], color="C2")

# set up and show turbine layout plot
axis("equal")
xlim(5500, 11000)
ylim(-500, 7000)
plt.show()
savefig("ipoptresults_turbinelayout.png")

# create obj value history plots
plt.clf()
plt.gcf().gca().plot(1:length(funcalls_AEP), funcalls_AEP)
title("Objective Value History")
xlabel("Function Call")
ylabel("AEP (TWh)")
savefig("ipoptresults_plot_AEP_history_by_funcall.png")

plt.clf()
plt.gcf().gca().plot(1:length(iter_AEP), iter_AEP)
title("Objective Value at Each Iteration")
xlabel("Iteration")
ylabel("AEP (TWh)")
savefig("ipoptresults_plot_AEP_history_by_iter.png")

# write text file with optimization summary info
io = open("ipoptresults_summary.txt", "a")
write(io, "number of turbines: ")
writedlm(io, nturbines)
write(io, "rotor diameter: ")
writedlm(io, rotor_diameter[1])
write(io, "starting AEP value (MWh): ")
writedlm(io, aep_wrapper(x_initial)[1]*1E5)
write(io, "Finished in : ")
writedlm(io, clkt)
write(io, "end AEP value (MWh): ")
writedlm(io, aep_wrapper(xopt)[1]*1E5)
write(io, "\noptimal turbine coordinates:\n")
writedlm(io, [turbine_x turbine_y])
write(io, "\neach function call value\n")
writedlm(io, funcalls_AEP)
close(io)

# write results to csv files
dataforcsv_xopt = DataFrame(xopt = xopt)
dataforcsv_funceval = DataFrame(function_value = funcalls_AEP)
dataforcsv_funceval_at_iteration = DataFrame(function_value = iter_AEP)
CSV.write("cs3_xopt.csv", dataforcsv_xopt)
CSV.write("cs3_functionvalue_log.csv", dataforcsv_funceval)
CSV.write("cs3_functionvalue_at_iteration_log.csv", dataforcsv_funceval_at_iteration)

# # write results to a YAML file (not complete)
# dataforyaml = [
#     "optimal_turbine_coordinates" => [turbine_x turbine_y]
#     "function_evaluation_log" => funcalls_AEP
#     "function_evaluation_at_iterations_log" => iter_AEP
#     ]
# YAML.write_file("cs3_results.yml", dataforyaml)


text:
using Distances
using Clustering
using DelimitedFiles
using PdfFraud

# set up the paths
reprpath = "vae_runs/2"
reprpath = joinpath(PdfFraud.experiment_path, reprpath, "representations")

# collect the files and encodings
encfs = filter(x->x!="index.csv",readdir(reprpath));
encodings = map(x->readdlm(x, ','), joinpath.(reprpath, encfs));

# get labels
_, labels, _, _ = PdfFraud.load_pdf_data();
linds = map(x->x!=nothing, labels)
Y = Int.(labels[linds]);

# compute distances for one encoding and use it to compute silhouttes
enc = encodings[11];
X = Array(enc[linds,:]');
silhouettes_encs(data, labels) = silhouettes(labels, pairwise(Euclidean(), data, dims=2))
silhs = silhouettes_encs(X, Y);


text:
@testset "ESCB2" begin
  n = 3
  ε = .1
  reward = Distribution[Bernoulli(.5 + ((i == j) ? ε : 0.)) for i in 1:n, j in 1:n]

  if ! is_travis
    @testset "Exact" begin
      instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingLPSolver(Gurobi.Optimizer))

      @testset "Algorithm" begin
        Random.seed!(1)
        n_rounds = 20
        s = simulate(instance, ESCB2(ESCB2Exact()), n_rounds)
        @test s.round == n_rounds
        @test sum(values(s.arm_counts)) <= n_rounds * n
        @test s.regret ≈ 3.0
        @test s.reward ≈ 32.0
        @test s.regret + s.reward ≈ 35.0

        for arm in keys(s.arm_counts)
          @test s.arm_average_reward[arm] ≈ s.arm_reward[arm] / s.arm_counts[arm]
        end
      end

      @testset "Trace" begin
        n = 3
        ε = .1
        reward = Distribution[Bernoulli(.5 + ((i == j) ? ε : 0.)) for i in 1:n, j in 1:n]
        instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingLPSolver(Gurobi.Optimizer))

        n_rounds = 20
        test_policy_trace(instance, ESCB2(ESCB2Exact()), n_rounds, n_rounds * n, 3.0, 32.0)
      end
    end
  end

  @testset "Greedy" begin
    instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingMunkresSolver())

    @testset "Algorithm" begin
      Random.seed!(1)
      n_rounds = 20
      s = simulate(instance, ESCB2(ESCB2Greedy()), n_rounds)
      @test s.round == n_rounds
      @test sum(values(s.arm_counts)) == n_rounds * n
      @test s.regret ≈ 3.0
      @test s.reward ≈ 32.0
      @test s.regret + s.reward ≈ 35.0

      for arm in keys(s.arm_counts)
        @test s.arm_average_reward[arm] ≈ s.arm_reward[arm] / s.arm_counts[arm]
      end
    end

    @testset "Trace" begin
      n = 3
      ε = .1
      reward = Distribution[Bernoulli(.5 + ((i == j) ? ε : 0.)) for i in 1:n, j in 1:n]
      instance = UncorrelatedPerfectBipartiteMatching(reward, PerfectBipartiteMatchingMunkresSolver())

      n_rounds = 20
      test_policy_trace(instance, ESCB2(ESCB2Greedy()), n_rounds, n_rounds * n, 3.0, 32.0)
    end
  end

  @testset "Budgeted (all at once)" begin
    # Not available for many instances, for now…
    m = 2
    r = Distribution[Bernoulli(.5 + ((i % 2 == 0) ? ε : 0.)) for i in 1:n]
    i = MSet(r, m, MSetAlgosSolver())

    @testset "Algorithm" begin
      Random.seed!(1)
      n_rounds = 20
      s = simulate(i, ESCB2(ESCB2Budgeted(.1, true)), n_rounds)
      @test s.round == n_rounds
      @test sum(values(s.arm_counts)) == n_rounds * m
      @test s.regret ≈ 1.5
      @test s.reward ≈ 26.0
      @test s.regret + s.reward ≈ 27.5

      for arm in keys(s.arm_counts)
        @test s.arm_average_reward[arm] ≈ s.arm_reward[arm] / s.arm_counts[arm]
      end
    end

    @testset "Trace" begin
      n_rounds = 20
      test_policy_trace(i, ESCB2(ESCB2Budgeted(.1, true)), n_rounds, n_rounds * m, 1.5, 26.0)
    end
  end

  @testset "Budgeted (one by one)" begin
    # Not available for many instances, for now…
    m = 2
    r = Distribution[Bernoulli(.5 + ((i % 2 == 0) ? ε : 0.)) for i in 1:n]
    i = MSet(r, m, MSetAlgosSolver())

    @testset "Algorithm" begin
      Random.seed!(1)
      n_rounds = 20
      s = simulate(i, ESCB2(ESCB2Budgeted(.1, false)), n_rounds)
      @test s.round == n_rounds
      @test sum(values(s.arm_counts)) <= n_rounds * m
      @test s.regret ≈ 1.5
      @test s.reward ≈ 26.0
      @test s.regret + s.reward ≈ 27.5

      for arm in keys(s.arm_counts)
        @test s.arm_average_reward[arm] ≈ s.arm_reward[arm] / s.arm_counts[arm]
      end
    end

    @testset "Trace" begin
      n_rounds = 20
      test_policy_trace(i, ESCB2(ESCB2Budgeted(.1, false)), n_rounds, n_rounds * m, 1.5, 26.0)
    end
  end

  @testset "optimise_linear_sqrtlinear" begin
    # Test data coming from 1000 rounds, with ε <= 1/sqrt(1000) ≈ 0.03.

    m = 2
    distr = Distribution[Bernoulli(.5) for i in 1:5]
    i = MSet(distr, m, MSetAlgosSolver())
    ilp = ! is_travis && MSet(distr, m, MSetLPSolver(Gurobi.Optimizer))
    w = Dict(1 => 0.16666666666666666,2 => 0.8041237113402062,3 => 0.16666666666666666,4 => 0.7921348314606742,5 => 0.16666666666666666)
    s2 = Dict(1 => 1.0818164611145353,2 => 0.10037472319619399,3 => 1.0818164611145353,4 => 0.10939717022506537,5 => 1.0818164611145353)
    f_n = log(1000) / 2
    exact_obj = 3.0669981573660285

    algo = ESCB2Budgeted(.01, true)
    sol, rd = optimise_linear_sqrtlinear(i, algo, w, s2, f_n, 1000, with_trace=true)
    @test escb2_index(w, s2, f_n, sol) ≈ exact_obj
    @test rd.best_objective ≈ exact_obj

    if ! is_travis
      algo = ESCB2Exact()
      sol, rd = optimise_linear_sqrtlinear(ilp, algo, w, s2, f_n, 1000, with_trace=true)
      @test escb2_index(w, s2, f_n, sol) ≈ exact_obj
      @test rd.best_objective ≈ exact_obj
    end

    m = 25
    distr = Distribution[Bernoulli(.5) for i in 1:50]
    i = MSet(distr, m, MSetAlgosSolver())
    ilp = ! is_travis && MSet(distr, m, MSetLPSolver(Gurobi.Optimizer))
    w = Dict(1 => 0.24242424242424243,2 => 0.7777777777777778,3 => 0.35714285714285715,4 => 0.8048048048048048,5 => 0.21875,6 => 0.8308605341246291,7 => 0.13793103448275862,8 => 0.7685185185185185,9 => 0.21875,10 => 0.7689393939393939,11 => 0.1111111111111111,12 => 0.7692307692307693,13 => 0.24242424242424243,14 => 0.7877813504823151,15 => 0.04,16 => 0.7509433962264151,17 => 0.24242424242424243,18 => 0.827485380116959,19 => 0.07692307692307693,20 => 0.7813765182186235,21 => 0.2647058823529412,22 => 0.7592592592592593,23 => 0.13793103448275862,24 => 0.8245614035087719,25 => 0.1111111111111111,26 => 0.8054711246200608,27 => 0.16666666666666666,28 => 0.7746913580246914,29 => 0.1935483870967742,30 => 0.8033333333333333,31 => 0.2972972972972973,32 => 0.8230088495575221,33 => 0.21875,34 => 0.7781350482315113,35 => 0.24242424242424243,36 => 0.8011695906432749,37 => 0.13793103448275862,38 => 0.8176470588235294,39 => 0.1935483870967742,40 => 0.8159509202453987,41 => 0.07692307692307693,42 => 0.8147058823529412,43 => 0.2972972972972973,44 => 0.7466666666666667,45 => 0.21875,46 => 0.7704918032786885,47 => 0.16666666666666666,48 => 0.7647058823529411,49 => 0.21875,50 => 0.8224852071005917)
    s2 = Dict(1 => 5.434941671158229,2 => 0.6227537331535471,3 => 4.270311313052894,4 => 0.5385978232679326,5 => 5.604783598381924,6 => 0.5322049707662361,7 => 6.184588798214537,8 => 0.5535588739142641,9 => 5.604783598381924,10 => 0.6793677088947786,11 => 6.642706486971169,12 => 0.6898195198008522,13 => 5.434941671158229,14 => 0.5766979908302944,15 => 7.1741230059288625,16 => 0.6768040571631002,17 => 5.434941671158229,18 => 0.5244241963398292,19 => 6.898195198008522,20 => 0.7261258103166864,21 => 5.275090445535929,22 => 0.5535588739142641,23 => 6.184588798214537,24 => 0.5244241963398292,25 => 6.642706486971169,26 => 0.5451461250705822,27 => 5.978435838274052,28 => 0.5535588739142641,29 => 5.78558306929747,30 => 0.5978435838274052,31 => 4.847380409411394,32 => 0.5290651184313321,33 => 5.604783598381924,34 => 0.5766979908302944,35 => 5.434941671158229,36 => 0.5244241963398292,37 => 6.184588798214537,38 => 0.5275090445535928,39 => 5.78558306929747,40 => 0.5501628072031336,41 => 6.898195198008522,42 => 0.5275090445535928,43 => 4.847380409411394,44 => 0.5978435838274052,45 => 5.604783598381924,46 => 0.7350535866730392,47 => 5.978435838274052,48 => 0.5552726784774661,49 => 5.604783598381924,50 => 0.5306303998468094)
    f_n = log(1000) / 2
    exact_obj = 29.29724069607937

    ε = 0.014581
    algo = ESCB2Budgeted(ε, true)
    sol, rd = optimise_linear_sqrtlinear(i, algo, w, s2, f_n, 1000, with_trace=true)
    @test escb2_index(w, s2, f_n, sol) ≈ exact_obj atol=ε
    @test rd.best_objective ≈ exact_obj atol=ε

    if ! is_travis
      algo = ESCB2Exact()
      sol, rd = optimise_linear_sqrtlinear(ilp, algo, w, s2, f_n, with_trace=true)
      @test escb2_index(w, s2, f_n, sol) ≈ exact_obj
      @test rd.best_objective ≈ exact_obj
    end
  end
end


text:
using Pkg
Pkg.activate(pwd())

using Merly

server = Merly.app()

@page "/" ""

@page "/user/:id>" "{{id}}"

@route POST "/user" begin
    res.body = ""
end

server.start(config = Dict("host" => "0.0.0.0", "port" => 3000), verbose = true)

max_stars_repo_name: paulcsmith/web-frameworks

max_issues_repo_name: paulcsmith/web-frameworks

max_forks_repo_name: paulcsmith/web-frameworks


text:
mutable struct HSDE{T1<:ProximableFunction,T2<:DualConeProduct}
    indAffine::T1
    indCones::T2
    n::Int64
end

function HSDE(model::FOSMathProgModel; direct=false)
    m,n = size(model.A)

    S1 = if direct
            # Uses direct solver from ProximalOperators
            Q = [spzeros(n,n) model.A'      model.c;
                -model.A      spzeros(m,m)  model.b;
                -model.c'     -model.b'     0      ]
            IndAffine([sparse(Q) -I], zeros(size(Q,1)))
        else
            Q = HSDEMatrixQ(model.A, model.b, model.c)
            # Q = [spzeros(n,n) model.A'      model.c;
            #     -model.A      spzeros(m,m)  model.b;
            #     -model.c'     -model.b'     0      ]
            # Using CG on KKT system
            AffinePlusLinear(Q, zeros(size(Q,1)), zeros(size(Q,1)), 1, decreasing_accuracy=true)
        end
    S2 = DualConeProduct(model.K1,model.K2)
    m,n = S2.m, S2.n
    status_generator = (mo, checki, eps, verbose, debug) ->
        HSDEStatus(m, n, 0, mo, :Continue, checki, eps, verbose, false, direct, time_ns(), model.init_duration, debug)
    return HSDE(S1, S2, 2*(size(model.A,1)+size(model.A,2)+1)), status_generator
end

# TODO Unclear why this is nessesary/differenet
getinitialvalue(model::FOSMathProgModel) = HSDE_getinitialvalue(model)

getinitialvalue(model::FOSMathProgModel, alg, data) =
    HSDE_getinitialvalue(model)

populate_solution(model::FOSMathProgModel, alg, data, x, status) =
    HSDE_populatesolution(model, x, status)

function HSDE_getinitialvalue(model::FOSMathProgModel)
    m, n = size(model.A)
    l = m+n+1
    x = zeros(2*l)
    x[l]  = 1.0
    x[2l] = 1.0
    return x
end

function HSDE_populatesolution(model::FOSMathProgModel, x, status::HSDEStatus)
    m, n = size(model.A)
    l = m+n+1
    @assert length(x) == 2l
    τ = x[l]
    κ = x[2l]
    #TODO Eveluate status
    endstatus = status.status
    if endstatus == :Continue
        endstatus = :Indeterminate
    end
    solution = Solution(x[1:n]/τ, x[n+1:n+m]/τ, x[l+n+1:l+n+m]/τ, endstatus)
end


text:
# Playing Cards: Unit Tests
using Test
using PlayingCards

@testset "card.jl" begin
    a = Card(1, 1)
    b = Card(7, ♡)
    c = Card(jack, ♠)
    d = Card(Seven, 2)
    e = Card(Five, clubs)
    f = Card(five, Clubs)
    g = 10♡
    h = 4Diamonds
    i = 6clubs

    @test (a > b) == true
    @test (d < c) == true
    @test (e === f) == true
    @test name(g) === "Ten of Hearts"
    @test name(h) === "Four of Diamonds"
    @test name(i) === "Six of Clubs"
end

@testset "deck.jl" begin
    deck = Deck()

    @test length(deck) == 52
    @test deck[end] == Card(King, ♠)
    @test length(deck[1:10]) == 10
    @test deal(deck) == Card(ace, Clubs)
end

@testset "hand.jl" begin
    deck = Deck()
    hand = Hand(cards = deal(deck, 5))

    @test length(deck) == 47
    @test length(hand) == 5
    @test name(hand[3]) == "Ace of Hearts"
    @test length(hand[1:3]) == 3
end

max_stars_repo_name: Joshua-Robison/PlayingCards

max_issues_repo_name: Joshua-Robison/PlayingCards

max_forks_repo_name: Joshua-Robison/PlayingCards


text:
function _build_species_table(model_dictionary::Dict{String,Any})::DataFrame

    # initialize -
    species_table = DataFrame(symbol=String[], type=Symbol[], compartment=String[], sequence=Union{Missing,String,FASTA.Record}[])

    try

        # ok, so lets build an intermediate representation that is a bunch of DataFrames -
        list_of_species_dictionaries = model_dictionary["list_of_model_species"]
        for species_dictionary in list_of_species_dictionaries

            # grab -
            local_data_row = Union{String,Symbol,Missing,FASTA.Record}[]
            push!(local_data_row, species_dictionary["symbol"])
            push!(local_data_row, Symbol(species_dictionary["type"]))
            push!(local_data_row, species_dictionary["compartment"])

            # Lets load the sequence -
            seq_path = species_dictionary["sequence"]
            if (seq_path == "")
                push!(local_data_row, missing)
            else
                open(FASTA.Reader, seq_path) do reader
                    for record in reader
                        push!(local_data_row, record)
                    end
                end
            end

            # push into data frame -
            push!(species_table, tuple(local_data_row...))
        end

        # return -
        return species_table
    catch error
        rethrow(error)
    end
end

# == PUBLIC METHODS BELOW HERE ==================================================================== #

"""
    read_model_document(model::VLJuliaModelObject) -> VLResult

This documentation is going to be awesome.
"""
function read_model_document(model::VLJuliaModelObject)::VLResult

    # initialize -
    intermediate_representation_dictionary = Dict{String,Any}()

    try

        # where is the main file located?
        path_to_model_file = model.path_to_model_file

        # load the model JSON file -
        json_model_dictionary = JSON.parsefile(path_to_model_file)

        # get model species, transscription and translation models -
        model_species_table = _build_species_table(json_model_dictionary)
        list_of_transcription_models = json_model_dictionary["list_of_transcription_models"]
        list_of_translation_models = json_model_dictionary["list_of_translation_models"]
        system_type = json_model_dictionary["system"]["system_type"]

        # how many species do we have?
        (number_of_species, _) = size(model_species_table)

        # setup up the system dimension -
        intermediate_representation_dictionary["number_of_species"] = number_of_species
        intermediate_representation_dictionary["number_of_transcription_models"] = length(list_of_transcription_models)
        intermediate_representation_dictionary["number_of_translation_models"] = length(list_of_translation_models)
        intermediate_representation_dictionary["system_type"] = system_type

        # Build a species table -
        intermediate_representation_dictionary["model_species_table"] = model_species_table

        # Grab the system dictionary -
        intermediate_representation_dictionary["system_dictionary"] = json_model_dictionary["system"]

        # Grab the transcription models -
        intermediate_representation_dictionary["list_of_transcription_models"] = list_of_transcription_models

        # Grab the translation models -
        intermediate_representation_dictionary["list_of_translation_models"] = list_of_translation_models

        # return -
        return VLResult(intermediate_representation_dictionary)
    catch error
        return VLResult{Exception}(error)
    end
end
# == PUBLIC METHODS ABOVE HERE ==================================================================== #


text:
using Charlotte
using WebAssembly
using Base.Test

# @testset "Import" begin

@noinline function mathfun(x)
    # x + sin(x)
    2x
end
@noinline function mathfun1(x)
    4x
end
function mathfun2(x, y)
    mathfun(3x) + mathfun1(y)
end

m = wasm_module([mathfun => Tuple{Float64},
                 mathfun2 => Tuple{Float64, Float64}])

function docos(x)
    ccall((:jscos, "imports"), Float64, (Float64,), x)
end
m1 = wasm_module([docos => Tuple{Float64}])

function docos2(x)
    ccall((:jscos, "imports"), Float64, (Float64,Float64), x, 33.3)
end
m2 = wasm_module([docos2 => Tuple{Float64}])

# BROKEN stuff to try to test memory
# const s = "hello"
# function stringtest()
#     # length(s)
#     s[1]
# end
# ms = wasm_module([stringtest => Tuple{}])
# const a = [1,2,3]
# function arraytest()
#     a[3]
# end
# ma = wasm_module([arraytest => Tuple{}])

# ms = wasm_module([stringtest => Tuple{}])

## Better UI:
# m = @wasm begin
#     mathfun(Float64)
#     mathfun2(Float64, Float64)
#     # Enter more functions here...
# end

# write("test.wat", m)

# end

WA = WebAssembly

relu(x) = ifelse(x < 0, 0, x)
result = @code_wasm relu(1)
expected = WA.Func(Symbol("#relu_Int64"), [WA.i64], [WA.i64], [], WA.Block([WA.Const(0), WA.Local(0), WA.Local(0), WA.Const(0), WA.Op(WA.i64, :lt_s), WA.Select(), WA.Return()]))
@test result.params == expected.params
@test result.returns == expected.returns
@test result.locals == expected.locals
@test result.body.body == expected.body.body