text:
export LPHRepresentation

# No copy since I do not modify anything and a copy is done when building a polyhedron
mutable struct LPHRepresentation{T, MT<:AbstractMatrix{T}} <: MixedHRep{T}
    # lb <= Ax <= ub
    # l <= x <= u
    A::MT
    l::AbstractVector{T}
    u::AbstractVector{T}
    colleqs::BitSet
    colgeqs::BitSet
    coleqs::BitSet
    lb::AbstractVector{T}
    ub::AbstractVector{T}
    rowleqs::BitSet
    rowgeqs::BitSet
    roweqs::BitSet

    function LPHRepresentation{T, MT}(A::MT, l::AbstractVector{T}, u::AbstractVector{T}, lb::AbstractVector{T}, ub::AbstractVector{T}) where {T, MT<:AbstractMatrix{T}}
        if length(l) != length(u) || size(A, 2) != length(l)
            throw(DimensionMismatch("The length of l and u must be equal to the number of rows of A"))
        end
        if length(lb) != length(ub) || size(A, 1) != length(lb)
            throw(DimensionMismatch("The length of lb and ub must be equal to the number of columns of A"))
        end
        colleqs = BitSet()
        colgeqs = BitSet()
        coleqs = BitSet()
        for i in 1:size(A, 2)
            leq = u[i] < typemax(T)
            geq = l[i] > typemin(T)
            if leq && geq && _isapprox(l[i], u[i])
                push!(coleqs, i)
            else
                if leq
                    push!(colleqs, i)
                end
                if geq
                    push!(colgeqs, i)
                end
            end
        end
        rowleqs = BitSet()
        rowgeqs = BitSet()
        roweqs = BitSet()
        for i in 1:size(A, 1)
            leq = ub[i] < typemax(T)
            geq = lb[i] > typemin(T)
            if leq && geq && _isapprox(lb[i], ub[i])
                push!(roweqs, i)
            else
                if leq
                    push!(rowleqs, i)
                end
                if geq
                    push!(rowgeqs, i)
                end
            end
        end
        new{T, typeof(A)}(A, l, u, colleqs, colgeqs, coleqs, lb, ub, rowleqs, rowgeqs, roweqs)
    end
end
FullDim(rep::LPHRepresentation) = size(rep.A, 2)

LPHRepresentation(A::AbstractMatrix{T}, l::AbstractVector{T}, u::AbstractVector{T}, lb::AbstractVector{T}, ub::AbstractVector{T}) where {T <: Real} = LPHRepresentation{T, typeof(A)}(A, l, u, lb, ub)
function LPHRepresentation(A::AbstractMatrix, l::AbstractVector, u::AbstractVector, lb::AbstractVector, ub::AbstractVector)
    T = promote_type(eltype(A), eltype(l), eltype(u), eltype(lb), eltype(ub))
    AT = AbstractMatrix{T}(A)
    LPHRepresentation{T, typeof(AT)}(AT, AbstractVector{T}(l), AbstractVector{T}(u), AbstractVector{T}(lb), AbstractVector{T}(ub))
end

hvectortype(::Type{LPHRepresentation{T, MT}}) where {T, MT} = vectortype(MT)
similar_type(::Type{LPHRepresentation{S, MT}}, ::FullDim, ::Type{T}) where {S, T, MT} = LPHRepresentation{T, similar_type(MT, T)}
fulltype(::Type{LPHRepresentation{T, MT}}) where {T, MT} = LPHRepresentation{T, MT}

LPHRepresentation(h::HRep{T}) where {T} = LPHRepresentation{T}(h)
LPHRepresentation{T}(h::HRep) where {T} = convert(LPHRepresentation{T, hmatrixtype(typeof(h), T)}, h)

#function LPHRepresentation{T}(it::HRepIterator{T}) where {T}
#    A = Matrix{T}(length(it), N)
#    lb = Vector{T}(length(it))
#    ub = Vector{T}(length(it))
#    MPB.HighLevelInterface.warn_no_inf(T)
#    l = fill(typemin(T), N)
#    u = fill(typemax(T), N)
#    for (i, h) in enumerate(it)
#        A[i,:] = h.a
#        ub[i] = h.β
#        if islin(h)
#            lb[i] = ub[i]
#        else
#            lb[i] = typemin(T)
#        end
#    end
#    LPHRepresentation{T}(A, l, u, lb, ub)
#end
function LPHRepresentation{T, MT}(d::FullDim,
                                  hyperplanes::ElemIt{<:HyperPlane{T}},
                                  halfspaces::ElemIt{<:HalfSpace{T}}) where {T, MT}
    nhyperplane = length(hyperplanes)
    nhrep = nhyperplane + length(halfspaces)
    N = fulldim(d)
    A = emptymatrix(MT, nhrep, N)
    lb = Vector{T}(undef, nhrep)
    ub = Vector{T}(undef, nhrep)
    MPB.HighLevelInterface.warn_no_inf(T)
    l = fill(typemin(T), N)
    u = fill(typemax(T), N)
    for (i, h) in enumerate(hyperplanes)
        A[i,:] = h.a
        lb[i] = h.β
        ub[i] = h.β
    end
    for (i, h) in enumerate(halfspaces)
        A[nhyperplane+i,:] = h.a
        lb[nhyperplane+i] = typemin(T)
        ub[nhyperplane+i] = h.β
    end
    LPHRepresentation{T, MT}(A, l, u, lb, ub)
end

Base.copy(lp::LPHRepresentation{T, MT}) where {T, MT} = LPHRepresentation{T, MT}(copy(lp.A), copy(lp.l), copy(lp.u), copy(lp.colleqs), copy(lp.colgeqs), copy(lp.coleqs), copy(lp.lb), copy(lp.ub), copy(lp.rowleqs), copy(lp.rowgeqs), copy(lp.roweqs))

Base.length(idxs::Indices{T, <:HyperPlane{T}, LPHRepresentation{T}}) where {T} = length(idxs.rep.coleqs) + length(idxs.rep.roweqs)
Base.length(idxs::Indices{T, <:HalfSpace{T}, LPHRepresentation{T}}) where {T} = length(idxs.rep.colleqs) + length(idxs.rep.colgeqs) + length(idxs.rep.rowleqs) + length(idxs.rep.rowgeqs)

# m rows (constraints) and n columns (variables)
# state : colrow, i, lgeq
# * colrow :
#   - 1 : columns (variables)
#   - 2 : rows (constraints)
#   - 3 : done
# * i : index of row or column
# * lgeq :
#   - 1 : ⟨a, x⟩ ≤ β
#   - 2 : ⟨a, x⟩ ≥ β
#   - 3 : ⟨a, x⟩ = β

# HyperPlane index j :
#   1 <= j <= n   : column j
# n+1 <= j <= n+m : row    j-n
function _index2state(lp, idx::HyperPlaneIndex)
    m = size(lp.A, 1)
    n = size(lp.A, 2)
    j = idx.value
    @assert j >= 1
    if j <= n
        1, j, 3
    elseif j <= n+m
        2, j-n, 3
    else
        3, 0, 0
    end
end

# HalfSpace index k :
# first bit :
# - 0 : lgeq = 1
# - 1 : lgeq = 2
# remaining bits : j
#   1 <= j+1 <= n   : column j
# n+1 <= j+1 <= n+m : row    j-n
function _index2state(lp, idx::HalfSpaceIndex)
    m = size(lp.A, 1)
    n = size(lp.A, 2)
    k = idx.value-1
    lgeq = (k & 1) + 1
    j = (k >> 1) + 1 # drop first bit
    @assert j >= 1
    if j <= n
        1, j, lgeq
    elseif j <= n+m
        2, j-n, lgeq
    else
        3, 0, 0
    end
end

function Base.isvalid(lp::LPHRepresentation{T}, idx::HIndex{T}) where {T}
    colrow, i, lgeq = _index2state(lp, idx)
    if colrow == 1
        lgeqs = (lp.colleqs, lp.colgeqs, lp.coleqs)
    else
        lgeqs = (lp.rowleqs, lp.rowgeqs, lp.roweqs)
    end
    1 <= colrow <= 2 && i in lgeqs[lgeq]
end

function getaβ(lp::LPHRepresentation{T}, idx::HIndex{T}) where {T}
    colrow, i, lgeq = _index2state(lp, idx)
    if colrow == 1
        a = origin(hvectortype(typeof(lp)), fulldim(lp))
        a[i] = lgeq == 2 ? -one(T) : one(T)
        β = lgeq == 2 ? -lp.l[i] : lp.u[i]
    else
        @assert colrow == 2
        a = lgeq == 2 ? -lp.A[i,:] : lp.A[i,:]
        β = lgeq == 2 ? -lp.lb[i] : lp.ub[i]
    end
    a, β
end
done(idxs::HIndices{T, <:LPHRepresentation{T}}, idx::HIndex{T}) where {T} = _index2state(idxs.rep, idx)[1] == 3
Base.get(lp::LPHRepresentation{T}, idx::HIndex{T}) where {T} = valuetype(idx)(getaβ(lp, idx)...)

dualtype(::Type{<:LPHRepresentation{T}}, ::Type{AT}) where {T, AT} = dualtype(Intersection{T, AT, Int}, AT)
dualfullspace(h::LPHRepresentation, d::FullDim, ::Type{T}, ::Type{AT}) where {T, AT} = dualfullspace(Intersection{T, AT, Int}, d, T, AT)


text:
using Flux

const FluxLayer = Union{
    Chain, Dense, Maxout, typeof(RNN), typeof(LSTM), typeof(GRU), Conv, ConvTranspose, MaxPool, MeanPool,
    DepthwiseConv, Dropout, AlphaDropout, LayerNorm, BatchNorm, InstanceNorm, GroupNorm,
    Hook
}

"""
    AttrType

Enum for layer field types.

- `ParamAttr`

    Field with trainable parameters, e.g. `Dense.W`
- `ChildAttr`

    Field that contains children layers, e.g. `Chain.layers`
- `SettingAttr`

    Fields with hyperparameters that define the layer structure, e.g. `Conv.pad`
"""
@enum AttrType begin
    ParamAttr
    SettingAttr
    ChildAttr
end

"""
    attrtype(layer::L, name::Val)::AttrType

Dispatch on a `layer`, and a field name (wrapped in a Val)
to determine what `AttrType` the field is.

If there is no custom method defined, for a layer type and field name,
`attrtype` will additionally dispatch on the field's value to give a
reasonable default.

By default the layers will be classified as following:

- `ParamAttr`

    Fields of type `AbstractArray{T,<:AbstractFloat}` are treated as trainable
    parameters

- `ChildAttr`

    Fields that are a Flux layer are treated as children.

- `SettingAttr`

    Any other

To customize for your type, define a method like

```julia
attrtype(::MyLayer, ::Val{:myfieldname}) = ChildAttr
```
"""
attrtype(::T, ::Val, ::Any) where T = SettingAttr
attrtype(::T, ::Val, ::AbstractArray) where T = ParamAttr
attrtype(::T, ::Val, ::FluxLayer) where T = ChildAttr
attrtype(::T, ::Val, ::NTuple{N, FluxLayer}) where {N, T} = ChildAttr

# dispatch on
attrtype(layer::Any, field::Val) = attrtype(layer, field, getfield(layer, getval(field)))
attrtype(::Chain, ::Val{:layers}) = ChildAttr
attrtype(::Hook, ::Val{:layer}) = ChildAttr

getval(::Val{T}) where T = T

ext: jl


text:
#---------------------------------------------------------
# # [Shepp-Logan Phantoms](@id 07-shepp)
#---------------------------------------------------------

#=
This page illustrates the Shepp-Logan phantoms in the Julia package
[`ImagePhantoms`](https://github.com/JuliaImageRecon/ImagePhantoms.jl).

This page was generated from a single Julia file:
[07-shepp.jl](@__REPO_ROOT_URL__/07-shepp.jl).
=#

#md # In any such Julia documentation,
#md # you can access the source code
#md # using the "Edit on GitHub" link in the top right.

#md # The corresponding notebook can be viewed in
#md # [nbviewer](http://nbviewer.jupyter.org/) here:
#md # [`07-shepp.ipynb`](@__NBVIEWER_ROOT_URL__/07-shepp.ipynb),
#md # and opened in [binder](https://mybinder.org/) here:
#md # [`07-shepp.ipynb`](@__BINDER_ROOT_URL__/07-shepp.ipynb).

# ### Setup

# Packages needed here.

using ImagePhantoms
using ImageGeoms: ImageGeom, axesf
using MIRTjim: jim, prompt

# The following line is helpful when running this example.jl file as a script;
# this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? jim(:prompt, true) : prompt(:draw);

# ### Overview

# There are several variations of the Shepp-Logan phantom available.

# ### CT version

#=
Original version from:
Larry A Shepp, Benjamin F Logan,
"The Fourier reconstruction of a head section,"
IEEE Transactions on Nuclear Science, 21(3):21-42, June 1974.
[doi](http://doi.org/10.1109/TNS.1974.6499235)

This (default) version has low soft tissue contrast,
so it usually should be displayed
with a narrow window the using `clim` option:
=#

image1 = shepp_logan(256) # CT version by default
jim(image1, "SheppLogan (original CT version)", clim=(0.95, 1.05))

# ### Over-sampling

#=
When generating ellipse phantoms,
it is generally preferable
to "over-sample" the ellipse values within each pixel
to account for partial volume effect.
A factor of 3 over-sampling (along both axes) typically suffices,
so this factor is the default for the `shepp_logan` method.
Here is how the phantom image looks without over-sampling:
=#

image1o = shepp_logan(256; oversample=1)
jim(image1o, "No over sampling", clim=(0.95,1.05))

#=
Note that boundaries of the interior ellipses look smoother
in the original version with the default 3× over-sampling.
Most of the remaining examples use the recommended default over-sampling.
=#

# ### Toft version

#=
This version is from:
Toft, Peter Aundal & Sørensen, John Aasted
"The Radon transform-theory and implementation,"
Technical University of Denmark (DTU), 1996. Page 201.
[pdf](https://files.openpdfs.org/1ra51GP6gJO.pdf)
=#

image2 = shepp_logan(256, SheppLoganToft())
jim(image2, "SheppLoganToft")

# ### Emission tomography version

#=
This version has low intensity for the skull
because typical PET/SPECT radiotracers
do not accumulate in bone regions.
It is probably also useful for MRI,
because typical MRI scans
have low signal from bone.
=#

image3 = shepp_logan(256, SheppLoganEmis())
jim(image3, "SheppLoganEmis")

# ### BrainWeb version

#=
This version was inspired by the
[BrainWeb](https://brainweb.bic.mni.mcgill.ca)
phantoms that have integer indices
for each of the different regions.
It should not be used directly,
but rather one should assign meaningful intensity values
to each of the integer indices.
=#

image4 = shepp_logan(256, SheppLoganBrainWeb())
jim(image4, "SheppLoganBrainWeb")

# For the BrainWeb version, there is no over-sampling by default,
# to preserve the integer indices.

# ### Disjoint middle ellipses

# Sometimes it can be more convenient
# to have the middle ellipses be non-overlapping:

ob = ellipse_parameters(SheppLoganBrainWeb(), disjoint=true)
ob[:,end] = 1:10
ob = Ellipse(ob)
x = LinRange(-0.4, 0.4, 206)
y = LinRange(-0.5, 0.5, 256)
oversample = 3
image5 = phantom(x, y, ob, oversample)
jim(x, y, image5, "Disjoint")

# ### Comedy version

image6 = shepp_logan(256, SouthPark(); fovs=(1,1))
jim(image6, "SouthPark")


text:
"""
    get_raw_moments(sol::EnsembleSolution, order::Int; naive::Bool=true, b::Int=2)

Given an `EnsembleSolution` of [DifferentialEquations ensemble simulation]
(https://diffeq.sciml.ai/stable/features/ensemble/#Performing-an-Ensemble-Simulation),
return a Dictionary of raw moments computed up to the specified `order` at each time step.

# Notes
- For example, the dictionary key `(2,0,1)` maps to an array containing the values of
  the raw moment ``μ_{201}`` at each time step.
- It is assumed that the time steps are all at the same time point for all trajectories
  (i.e., fixed `dt` used by the integrator or values were saved using `saveat`, as discussed [here]
  (https://diffeq.sciml.ai/stable/features/ensemble/#Time-steps-vs-time-points)).
- Moments are computed using [Cumulants.jl](https://github.com/iitis/Cumulants.jl) internally.
  The naive algorithm of moment tensor calculations (`naive=true`) is usually faster for small
  systems but the [proposed novel algorithm](https://arxiv.org/pdf/1701.05420.pdf) (`naive=false`)
  *should* be more efficient in case of many marginal variables. The [block size]
  (https://github.com/iitis/Cumulants.jl#block-size) `b` can also be specified for the
  novel algorithm and may have a significant effect on its performance.
- Only useful if higher order moments are needed: DifferentialEquations
  has a number of far more efficient and flexible ensemble statistics functions for
  means, variances and correlations, see [this tutorial]
  (https://diffeq.sciml.ai/stable/features/ensemble/#Analyzing-an-Ensemble-Experiment)
  for more details.
"""
function get_raw_moments(sol::EnsembleSolution, order::Int; naive::Bool=true, b::Int=2)

    # TODO: improve performance (it's rather poor now)
    # TODO: consider dropping the Dict and using DiffEqArrays

    if naive
        f_alg = naivemoment
        b = ()
    else
        f_alg = moment
    end

    no_t_pts = length(sol.u[1])
    N = length(sol.u[1].u[1])

    iter_all = construct_iter_all(N, order)
    iter_order = [filter(x -> sum(x) == m, iter_all) for m in 1:order]
    iter_μ = vcat(iter_order...)

    keys = Dict([iter => vcat([fill(i, n) for (i,n) in enumerate(iter)]...) for iter in iter_μ])
    μ = Dict([iter => Array{Float64}(undef, no_t_pts) for iter in iter_μ])

    for t_pt in 1:no_t_pts

        tslice = componentwise_vectors_timestep(sol, t_pt)
        tslice = hcat(tslice...) / 1

        for m in 1:order
            data = Array(f_alg(tslice, m, b...))
            for iter in iter_order[m]
                μ[iter][t_pt] = data[keys[iter]...]
            end
        end

    end

    μ

end

"""
    get_central_moments(sol::EnsembleSolution, order::Int; naive::Bool=true, b::Int=2)

Given an `EnsembleSolution` of [DifferentialEquations ensemble simulation]
(https://diffeq.sciml.ai/stable/features/ensemble/#Performing-an-Ensemble-Simulation),
return a Dictionary of central moment estimates computed up to the specified `order`
at each time step. See the notes of [`get_raw_moments`](@ref) function for more information.
"""
function get_central_moments(sol::EnsembleSolution, order::Int; naive::Bool=true, b::Int=2)

    if naive
        f_alg = naivemoment
        b = ()
    else
        f_alg = moment
    end

    no_t_pts = length(sol.u[1])
    N = length(sol.u[1].u[1])

    iter_all = construct_iter_all(N, order)
    iter_order = [filter(x-> sum(x) == m, iter_all) for m in 1:order]
    iter_M = filter(x -> sum(x) > 1, iter_all)

    keys = Dict([iter => vcat([fill(i, n) for (i,n) in enumerate(iter)]...) for iter in iter_all])
    M = Dict([iter => Array{Float64}(undef, no_t_pts) for iter in iter_M])

    μ = Dict()
    μ[Tuple(fill(0, N))] = 1.

    for t_pt in 1:no_t_pts
        tslice = componentwise_vectors_timestep(sol, t_pt)
        tslice = hcat(tslice...) / 1

        for m in 1:order
            data = Array(f_alg(tslice, m, b...))
            for iter in iter_order[m]
                μ[iter] = data[keys[iter]...]
            end
        end
        M_temp = raw_to_central_moments(N, order, μ)
        for iter in iter_M
            M[iter][t_pt] = M_temp[iter]
        end
    end

    M

end

"""
    get_cumulants(sol::EnsembleSolution, order::Int; naive::Bool=true, b::Int=2)

Given an `EnsembleSolution` of [DifferentialEquations ensemble simulation]
(https://diffeq.sciml.ai/stable/features/ensemble/#Performing-an-Ensemble-Simulation),
return a Dictionary of cumulant estimates computed up to the specified `order`
at each time step. See the notes of [`get_raw_moments`](@ref) function for more information.
"""
function get_cumulants(sol::EnsembleSolution, order::Int; naive::Bool=true, b::Int=2)

    no_t_pts = length(sol.u[1])
    N = length(sol.u[1].u[1])

    iter_all = construct_iter_all(N, order)
    iter_order = [filter(x-> sum(x) == m, iter_all) for m in 1:order]
    iter_κ = vcat(iter_order...)

    keys = Dict([iter => vcat([fill(i, n) for (i,n) in enumerate(iter)]...) for iter in iter_κ])
    κ = Dict([iter => Array{Float64}(undef, no_t_pts) for iter in iter_κ])

    if naive

        for t_pt in 1:no_t_pts

            tslice = componentwise_vectors_timestep(sol, t_pt)
            tslice = hcat(tslice...) / 1

            for m in 1:order
                data = Array(naivecumulant(tslice, m))
                for iter in iter_order[m]
                    κ[iter][t_pt] = data[keys[iter]...]
                end
            end

        end

    else

        for t_pt in 1:no_t_pts

            tslice = componentwise_vectors_timestep(sol, t_pt)
            tslice = hcat(tslice...) / 1
            κ_tensor = cumulants(tslice, order, b)

            for m in 1:order

                data = Array(κ_tensor[m])

                for iter in iter_order[m]
                    κ[iter][t_pt] = data[keys[iter]...]
                end

            end

        end

    end

    κ

end

"""
    get_moments_FSP(sol::ODESolution, order::Int, moment_type::String)

Given an `ODESolution` obtained using [FiniteStateProjection.jl](https://github.com/kaandocal/FiniteStateProjection.jl),
return a Dictionary of moments computed up to the specified `order` at each time step.
Here, `moment_type` specifies the type of moments to be computed: available options are `raw`,
`central` or `cumulant`.

# Notes
- The `ODESolution` represents the time-evolution of the probability density function
  that is the solution of the Chemical Master Equation approximated using Finite State Projection
  algorithms. See the documentation of [FiniteStateProjection.jl](https://github.com/kaandocal/FiniteStateProjection.jl)
  for more information.
"""
function get_moments_FSP(sol::ODESolution, order::Int, moment_type::String)

    @assert any(moment_type .== ["raw", "central", "cumulant"]) "Moment type "*moment_type*" does not exist"

    state_space = size(sol)[1:end-1]
    N = length(state_space)
    no_t_pts = length(sol.u)

    iter_moments = MomentClosure.construct_iter_all(N, order)[2:end]
    moments = Dict([iter => Array{Float64}(undef, no_t_pts) for iter in iter_moments])

    iter_state = Iterators.product((0:i-1 for i in state_space)...)

    if moment_type == "raw"
        for t_pt in 1:no_t_pts
            tslice = sol[t_pt]
            for μ_ind in iter_moments
                moments[μ_ind][t_pt] = sum(prod(n_state.^μ_ind) * tslice[(n_state.+1)...] for n_state in iter_state)
            end
        end
    else
        μ = Dict{Tuple{fill(Int, N)...}, Float64}()
        μ[Tuple(fill(0, N))] = 1.

        for t_pt in 1:no_t_pts
            tslice = sol[t_pt]
            for μ_ind in iter_moments
                μ[μ_ind] = sum(prod(n_state.^μ_ind) * tslice[(n_state.+1)...] for n_state in iter_state)
            end

            if moment_type == "central"
                moment_temp = MomentClosure.raw_to_central_moments(N, order, μ)
            else
                moment_temp = MomentClosure.cumulants_to_raw_moments(N, order, μ)
            end

            for iter in iter_moments
                moments[iter][t_pt] = moment_temp[iter]
            end
        end
    end

    moments

end

"""
    deterministic_IC(u₀::Array{T, 1}, eqs::MomentEquations) where T<:Real

Given an array of initial molecule numbers and the corresponding moment equations,
return a mapping of each moment to its initial value under deterministic initial conditions.

# Notes
- The means are set to initial molecule numbers (as they take the values specified in
  `u₀` with probability one). The higher order raw moments are products of the corresponding
  powers of the means whereas the higher order central moments are simply zero.
- The ordering of `u₀` elements must be consistent with the ordering of species
  in the corresponding reaction system (can be checked with the `speciesmap` function).
- As higher-order moment functions under log-normal, gamma, derivative matching and
  the conditional closures involve moments raised to negative powers, setting initial
  molecule numbers of certain species to *zeros* will result in NaN errors when solving
  the ODEs (the specifics depend on the system at hand).
"""
function deterministic_IC(u₀::Array{T, 1}, eqs::MomentEquations) where T<:Real

    if eqs isa ClosedMomentEquations
        sys = eqs.open_eqs
    else
        sys = eqs
    end

    odes = eqs.odes
    N = sys.N
    if N != length(u₀)
        error("length of the passed IC vector and the number of species in the system are inconsistent")
    end

    μ_map = [sys.μ[iter] => u₀[i] for (i, iter) in enumerate(sys.iter_1)]

    no_states = length(odes.states)
    if typeof(sys) == CentralMomentEquations
        moment_map = [odes.states[i] => 0.0 for i in N+1:no_states]
    else
        reverse_μ = Dict(μ => iter for (iter, μ) in sys.μ)
        moment_map = [odes.states[i] => prod(u₀ .^ reverse_μ[odes.states[i]]) for i in N+1:no_states]
    end

    vcat(μ_map, moment_map)

end

"""
    format_moment_eqs(eqs::MomentEquations)

Given `MomentEquations`, return an array of formatted strings representing
each ODE that can be visualised using [Latexify](https://github.com/korsbo/Latexify.jl).
Although Latexify can be applied directly on [`ModelingToolkit.ODESystem`]
(https://mtk.sciml.ai/stable/systems/ODESystem/) (accessed by `eqs.odes`),
`format_moment_eqs` makes it visually more appealing: simplifies the symbolic
expressions further (expanding all terms), makes the time-dependence of moment
variables implicit (removes all `(t)`) and removes all trailing zeros (`2.0` → `2`).
"""
function format_moment_eqs(eqs::MomentEquations)

    odes = eqs.odes
    exprs  = []
    for i in 1:size(odes.eqs)[1]
        key = odes.states[i]
        eq = odes.eqs[i].rhs
        expr = "d"*string(key)*"/dt = "*string(eq)
        expr = replace(expr, "(t)"=>"")
        expr = replace(expr, ".0"=>"")
        push!(exprs, expr)
    end
    exprs

end

"""
    format_closure(eqs::ClosedMomentEquations)

Given [`ClosedMomentEquations`](@ref), return an array of formatted strings
representing closure functions for each higher order moment. The symbolic
expressions are formatted similarly to [`format_moment_eqs`](@ref) and can
be visualised using [Latexify](https://github.com/korsbo/Latexify.jl).

# Notes
- `include_all` argument specifies whether all higher-order moments (`include_all=true`)
  or only those encountered in the moment ODEs specifically (`include_all=false`,
  the default option) will be returned.
"""
function format_closure(eqs::ClosedMomentEquations; format_all::Bool=false)
    closure = eqs.closure
    exprs = []

    if format_all
        iter = keys(closure)
    else
        iter = setdiff(eqs.open_eqs.odes.states, eqs.odes.states)
    end

    for i in iter
        eq = closure[i]
        expr = string(i)*" = "*string(eq)
        expr = replace(expr, "(t)"=>"")
        expr = replace(expr, ".0"=>"")
        push!(exprs, expr)
    end
    exprs
end

@latexrecipe function f(eqs::MomentEquations, type=:equations; print_all=false)

    env --> :align
    starred --> true
    cdot --> false

    if type == :equations
        return format_moment_eqs(eqs)
    elseif type == :closure
        return format_closure(eqs, format_all=print_all)
    else
        error("supported arguments are only `:equations` or `:closure`")
    end

end


text:
# Julia wrapper for header: /usr/include/libavformat/avformat.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_get_packet,
    av_append_packet,
    av_stream_get_r_frame_rate,
    av_stream_set_r_frame_rate,
    av_stream_get_recommended_encoder_configuration,
    av_stream_set_recommended_encoder_configuration,
    av_stream_get_parser,
    av_stream_get_end_pts,
    av_format_get_probe_score,
    av_format_get_video_codec,
    av_format_set_video_codec,
    av_format_get_audio_codec,
    av_format_set_audio_codec,
    av_format_get_subtitle_codec,
    av_format_set_subtitle_codec,
    av_format_get_data_codec,
    av_format_set_data_codec,
    av_format_get_metadata_header_padding,
    av_format_set_metadata_header_padding,
    av_format_get_opaque,
    av_format_set_opaque,
    av_format_get_control_message_cb,
    av_format_set_control_message_cb,
    av_format_get_open_cb,
    av_format_set_open_cb,
    av_format_inject_global_side_data,
    av_fmt_ctx_get_duration_estimation_method,
    avformat_version,
    avformat_configuration,
    avformat_license,
    av_register_all,
    av_register_input_format,
    av_register_output_format,
    avformat_network_init,
    avformat_network_deinit,
    av_iformat_next,
    av_oformat_next,
    av_muxer_iterate,
    av_demuxer_iterate,
    avformat_alloc_context,
    avformat_free_context,
    avformat_get_class,
    avformat_new_stream,
    av_stream_add_side_data,
    av_stream_new_side_data,
    av_stream_get_side_data,
    av_new_program,
    avformat_alloc_output_context2,
    av_find_input_format,
    av_probe_input_format,
    av_probe_input_format2,
    av_probe_input_format3,
    av_probe_input_buffer2,
    av_probe_input_buffer,
    avformat_open_input,
    av_demuxer_open,
    avformat_find_stream_info,
    av_find_program_from_stream,
    av_program_add_stream_index,
    av_find_best_stream,
    av_read_frame,
    av_seek_frame,
    avformat_seek_file,
    avformat_flush,
    av_read_play,
    av_read_pause,
    avformat_close_input,
    avformat_write_header,
    avformat_init_output,
    av_write_frame,
    av_interleaved_write_frame,
    av_write_uncoded_frame,
    av_interleaved_write_uncoded_frame,
    av_write_uncoded_frame_query,
    av_write_trailer,
    av_guess_format,
    av_guess_codec,
    av_get_output_timestamp,
    av_hex_dump,
    av_hex_dump_log,
    av_pkt_dump2,
    av_pkt_dump_log2,
    av_codec_get_id,
    av_codec_get_tag,
    av_codec_get_tag2,
    av_find_default_stream_index,
    av_index_search_timestamp,
    av_add_index_entry,
    av_url_split,
    av_dump_format,
    av_get_frame_filename2,
    av_get_frame_filename,
    av_filename_number_test,
    av_sdp_create,
    av_match_ext,
    avformat_query_codec,
    avformat_get_riff_video_tags,
    avformat_get_riff_audio_tags,
    avformat_get_mov_video_tags,
    avformat_get_mov_audio_tags,
    av_guess_sample_aspect_ratio,
    av_guess_frame_rate,
    avformat_match_stream_specifier,
    avformat_queue_attached_pictures,
    av_apply_bitstream_filters,
    avformat_transfer_internal_stream_timing_info,
    av_stream_get_codec_timebase

function av_get_packet(s, pkt, size::Integer)
    ccall((:av_get_packet, libavformat), Cint, (Ptr{AVIOContext}, Ptr{AVPacket}, Cint), s, pkt, size)
end

function av_append_packet(s, pkt, size::Integer)
    ccall((:av_append_packet, libavformat), Cint, (Ptr{AVIOContext}, Ptr{AVPacket}, Cint), s, pkt, size)
end

function av_stream_get_r_frame_rate(s)
    ccall((:av_stream_get_r_frame_rate, libavformat), AVRational, (Ptr{AVStream},), s)
end

function av_stream_set_r_frame_rate(s, r::AVRational)
    ccall((:av_stream_set_r_frame_rate, libavformat), Cvoid, (Ptr{AVStream}, AVRational), s, r)
end

function av_stream_get_recommended_encoder_configuration(s)
    ccall((:av_stream_get_recommended_encoder_configuration, libavformat), Cstring, (Ptr{AVStream},), s)
end

function av_stream_set_recommended_encoder_configuration(s, configuration)
    ccall((:av_stream_set_recommended_encoder_configuration, libavformat), Cvoid, (Ptr{AVStream}, Cstring), s, configuration)
end

function av_stream_get_parser(s)
    ccall((:av_stream_get_parser, libavformat), Ptr{AVCodecParserContext}, (Ptr{AVStream},), s)
end

function av_stream_get_end_pts(st)
    ccall((:av_stream_get_end_pts, libavformat), Int64, (Ptr{AVStream},), st)
end

function av_format_get_probe_score(s)
    ccall((:av_format_get_probe_score, libavformat), Cint, (Ptr{AVFormatContext},), s)
end

function av_format_get_video_codec(s)
    ccall((:av_format_get_video_codec, libavformat), Ptr{AVCodec}, (Ptr{AVFormatContext},), s)
end

function av_format_set_video_codec(s, c)
    ccall((:av_format_set_video_codec, libavformat), Cvoid, (Ptr{AVFormatContext}, Ptr{AVCodec}), s, c)
end

function av_format_get_audio_codec(s)
    ccall((:av_format_get_audio_codec, libavformat), Ptr{AVCodec}, (Ptr{AVFormatContext},), s)
end

function av_format_set_audio_codec(s, c)
    ccall((:av_format_set_audio_codec, libavformat), Cvoid, (Ptr{AVFormatContext}, Ptr{AVCodec}), s, c)
end

function av_format_get_subtitle_codec(s)
    ccall((:av_format_get_subtitle_codec, libavformat), Ptr{AVCodec}, (Ptr{AVFormatContext},), s)
end

function av_format_set_subtitle_codec(s, c)
    ccall((:av_format_set_subtitle_codec, libavformat), Cvoid, (Ptr{AVFormatContext}, Ptr{AVCodec}), s, c)
end

function av_format_get_data_codec(s)
    ccall((:av_format_get_data_codec, libavformat), Ptr{AVCodec}, (Ptr{AVFormatContext},), s)
end

function av_format_set_data_codec(s, c)
    ccall((:av_format_set_data_codec, libavformat), Cvoid, (Ptr{AVFormatContext}, Ptr{AVCodec}), s, c)
end

function av_format_get_metadata_header_padding(s)
    ccall((:av_format_get_metadata_header_padding, libavformat), Cint, (Ptr{AVFormatContext},), s)
end

function av_format_set_metadata_header_padding(s, c::Integer)
    ccall((:av_format_set_metadata_header_padding, libavformat), Cvoid, (Ptr{AVFormatContext}, Cint), s, c)
end

function av_format_get_opaque(s)
    ccall((:av_format_get_opaque, libavformat), Ptr{Cvoid}, (Ptr{AVFormatContext},), s)
end

function av_format_set_opaque(s, opaque)
    ccall((:av_format_set_opaque, libavformat), Cvoid, (Ptr{AVFormatContext}, Ptr{Cvoid}), s, opaque)
end

function av_format_get_control_message_cb(s)
    ccall((:av_format_get_control_message_cb, libavformat), av_format_control_message, (Ptr{AVFormatContext},), s)
end

function av_format_set_control_message_cb(s, callback::av_format_control_message)
    ccall((:av_format_set_control_message_cb, libavformat), Cvoid, (Ptr{AVFormatContext}, av_format_control_message), s, callback)
end

function av_format_get_open_cb(s)
    ccall((:av_format_get_open_cb, libavformat), AVOpenCallback, (Ptr{AVFormatContext},), s)
end

function av_format_set_open_cb(s, callback::AVOpenCallback)
    ccall((:av_format_set_open_cb, libavformat), Cvoid, (Ptr{AVFormatContext}, AVOpenCallback), s, callback)
end

function av_format_inject_global_side_data(s)
    ccall((:av_format_inject_global_side_data, libavformat), Cvoid, (Ptr{AVFormatContext},), s)
end

function av_fmt_ctx_get_duration_estimation_method(ctx)
    ccall((:av_fmt_ctx_get_duration_estimation_method, libavformat), AVDurationEstimationMethod, (Ptr{AVFormatContext},), ctx)
end

function avformat_version()
    ccall((:avformat_version, libavformat), UInt32, ())
end

function avformat_configuration()
    ccall((:avformat_configuration, libavformat), Cstring, ())
end

function avformat_license()
    ccall((:avformat_license, libavformat), Cstring, ())
end

function av_register_all()
    ccall((:av_register_all, libavformat), Cvoid, ())
end

function av_register_input_format(format)
    ccall((:av_register_input_format, libavformat), Cvoid, (Ptr{AVInputFormat},), format)
end

function av_register_output_format(format)
    ccall((:av_register_output_format, libavformat), Cvoid, (Ptr{AVOutputFormat},), format)
end

function avformat_network_init()
    ccall((:avformat_network_init, libavformat), Cint, ())
end

function avformat_network_deinit()
    ccall((:avformat_network_deinit, libavformat), Cint, ())
end

function av_iformat_next(f)
    ccall((:av_iformat_next, libavformat), Ptr{AVInputFormat}, (Ptr{AVInputFormat},), f)
end

function av_oformat_next(f)
    ccall((:av_oformat_next, libavformat), Ptr{AVOutputFormat}, (Ptr{AVOutputFormat},), f)
end

function av_muxer_iterate(opaque)
    ccall((:av_muxer_iterate, libavformat), Ptr{AVOutputFormat}, (Ptr{Ptr{Cvoid}},), opaque)
end

function av_demuxer_iterate(opaque)
    ccall((:av_demuxer_iterate, libavformat), Ptr{AVInputFormat}, (Ptr{Ptr{Cvoid}},), opaque)
end

function avformat_alloc_context()
    ccall((:avformat_alloc_context, libavformat), Ptr{AVFormatContext}, ())
end

function avformat_free_context(s)
    ccall((:avformat_free_context, libavformat), Cvoid, (Ptr{AVFormatContext},), s)
end

function avformat_get_class()
    ccall((:avformat_get_class, libavformat), Ptr{AVClass}, ())
end

function avformat_new_stream(s, c)
    ccall((:avformat_new_stream, libavformat), Ptr{AVStream}, (Ptr{AVFormatContext}, Ptr{AVCodec}), s, c)
end

function av_stream_add_side_data(st, _type::AVPacketSideDataType, data, size::Csize_t)
    ccall((:av_stream_add_side_data, libavformat), Cint, (Ptr{AVStream}, AVPacketSideDataType, Ptr{UInt8}, Csize_t), st, _type, data, size)
end

function av_stream_new_side_data(stream, _type::AVPacketSideDataType, size::Integer)
    ccall((:av_stream_new_side_data, libavformat), Ptr{UInt8}, (Ptr{AVStream}, AVPacketSideDataType, Cint), stream, _type, size)
end

function av_stream_get_side_data(stream, _type::AVPacketSideDataType, size)
    ccall((:av_stream_get_side_data, libavformat), Ptr{UInt8}, (Ptr{AVStream}, AVPacketSideDataType, Ptr{Cint}), stream, _type, size)
end

function av_new_program(s, id::Integer)
    ccall((:av_new_program, libavformat), Ptr{AVProgram}, (Ptr{AVFormatContext}, Cint), s, id)
end

function avformat_alloc_output_context2(ctx, oformat, format_name, filename)
    ccall((:avformat_alloc_output_context2, libavformat), Cint, (Ptr{Ptr{AVFormatContext}}, Ptr{AVOutputFormat}, Cstring, Cstring), ctx, oformat, format_name, filename)
end

function av_find_input_format(short_name)
    ccall((:av_find_input_format, libavformat), Ptr{AVInputFormat}, (Cstring,), short_name)
end

function av_probe_input_format(pd, is_opened::Integer)
    ccall((:av_probe_input_format, libavformat), Ptr{AVInputFormat}, (Ptr{AVProbeData}, Cint), pd, is_opened)
end

function av_probe_input_format2(pd, is_opened::Integer, score_max)
    ccall((:av_probe_input_format2, libavformat), Ptr{AVInputFormat}, (Ptr{AVProbeData}, Cint, Ptr{Cint}), pd, is_opened, score_max)
end

function av_probe_input_format3(pd, is_opened::Integer, score_ret)
    ccall((:av_probe_input_format3, libavformat), Ptr{AVInputFormat}, (Ptr{AVProbeData}, Cint, Ptr{Cint}), pd, is_opened, score_ret)
end

function av_probe_input_buffer2(pb, fmt, url, logctx, offset::Integer, max_probe_size::Integer)
    ccall((:av_probe_input_buffer2, libavformat), Cint, (Ptr{AVIOContext}, Ptr{Ptr{AVInputFormat}}, Cstring, Ptr{Cvoid}, UInt32, UInt32), pb, fmt, url, logctx, offset, max_probe_size)
end

function av_probe_input_buffer(pb, fmt, url, logctx, offset::Integer, max_probe_size::Integer)
    ccall((:av_probe_input_buffer, libavformat), Cint, (Ptr{AVIOContext}, Ptr{Ptr{AVInputFormat}}, Cstring, Ptr{Cvoid}, UInt32, UInt32), pb, fmt, url, logctx, offset, max_probe_size)
end

function avformat_open_input(ps, url, fmt, options)
    ccall((:avformat_open_input, libavformat), Cint, (Ptr{Ptr{AVFormatContext}}, Cstring, Ptr{AVInputFormat}, Ref{Ptr{AVDictionary}}), ps, url, fmt, options)
end

function av_demuxer_open(ic)
    ccall((:av_demuxer_open, libavformat), Cint, (Ptr{AVFormatContext},), ic)
end

function avformat_find_stream_info(ic, options)
    ccall((:avformat_find_stream_info, libavformat), Cint, (Ptr{AVFormatContext}, Ptr{Ptr{AVDictionary}}), ic, options)
end

function av_find_program_from_stream(ic, last, s::Integer)
    ccall((:av_find_program_from_stream, libavformat), Ptr{AVProgram}, (Ptr{AVFormatContext}, Ptr{AVProgram}, Cint), ic, last, s)
end

function av_program_add_stream_index(ac, progid::Integer, idx::Integer)
    ccall((:av_program_add_stream_index, libavformat), Cvoid, (Ptr{AVFormatContext}, Cint, UInt32), ac, progid, idx)
end

function av_find_best_stream(ic, _type::AVMediaType, wanted_stream_nb::Integer, related_stream::Integer, decoder_ret, flags::Integer)
    ccall((:av_find_best_stream, libavformat), Cint, (Ptr{AVFormatContext}, AVMediaType, Cint, Cint, Ptr{Ptr{AVCodec}}, Cint), ic, _type, wanted_stream_nb, related_stream, decoder_ret, flags)
end

function av_read_frame(s, pkt)
    ccall((:av_read_frame, libavformat), Cint, (Ptr{AVFormatContext}, Ptr{AVPacket}), s, pkt)
end

function av_seek_frame(s, stream_index::Integer, timestamp::Int64, flags::Integer)
    ccall((:av_seek_frame, libavformat), Cint, (Ptr{AVFormatContext}, Cint, Int64, Cint), s, stream_index, timestamp, flags)
end

function avformat_seek_file(s, stream_index::Integer, min_ts::Int64, ts::Int64, max_ts::Int64, flags::Integer)
    ccall((:avformat_seek_file, libavformat), Cint, (Ptr{AVFormatContext}, Cint, Int64, Int64, Int64, Cint), s, stream_index, min_ts, ts, max_ts, flags)
end

function avformat_flush(s)
    ccall((:avformat_flush, libavformat), Cint, (Ptr{AVFormatContext},), s)
end

function av_read_play(s)
    ccall((:av_read_play, libavformat), Cint, (Ptr{AVFormatContext},), s)
end

function av_read_pause(s)
    ccall((:av_read_pause, libavformat), Cint, (Ptr{AVFormatContext},), s)
end

function avformat_close_input(s)
    ccall((:avformat_close_input, libavformat), Cvoid, (Ptr{Ptr{AVFormatContext}},), s)
end

function avformat_write_header(s, options)
    ccall((:avformat_write_header, libavformat), Cint, (Ptr{AVFormatContext}, Ref{Ptr{AVDictionary}}), s, options)
end

function avformat_init_output(s, options)
    ccall((:avformat_init_output, libavformat), Cint, (Ptr{AVFormatContext}, Ref{Ptr{AVDictionary}}), s, options)
end

function av_write_frame(s, pkt)
    ccall((:av_write_frame, libavformat), Cint, (Ptr{AVFormatContext}, Ptr{AVPacket}), s, pkt)
end

function av_interleaved_write_frame(s, pkt)
    ccall((:av_interleaved_write_frame, libavformat), Cint, (Ptr{AVFormatContext}, Ptr{AVPacket}), s, pkt)
end

function av_write_uncoded_frame(s, stream_index::Integer, frame)
    ccall((:av_write_uncoded_frame, libavformat), Cint, (Ptr{AVFormatContext}, Cint, Ptr{AVFrame}), s, stream_index, frame)
end

function av_interleaved_write_uncoded_frame(s, stream_index::Integer, frame)
    ccall((:av_interleaved_write_uncoded_frame, libavformat), Cint, (Ptr{AVFormatContext}, Cint, Ptr{AVFrame}), s, stream_index, frame)
end

function av_write_uncoded_frame_query(s, stream_index::Integer)
    ccall((:av_write_uncoded_frame_query, libavformat), Cint, (Ptr{AVFormatContext}, Cint), s, stream_index)
end

function av_write_trailer(s)
    ccall((:av_write_trailer, libavformat), Cint, (Ptr{AVFormatContext},), s)
end

function av_guess_format(short_name, filename, mime_type)
    ccall((:av_guess_format, libavformat), Ptr{AVOutputFormat}, (Cstring, Cstring, Cstring), short_name, filename, mime_type)
end

function av_guess_codec(fmt, short_name, filename, mime_type, _type::AVMediaType)
    ccall((:av_guess_codec, libavformat), AVCodecID, (Ptr{AVOutputFormat}, Cstring, Cstring, Cstring, AVMediaType), fmt, short_name, filename, mime_type, _type)
end

function av_get_output_timestamp(s, stream::Integer, dts, wall)
    ccall((:av_get_output_timestamp, libavformat), Cint, (Ptr{AVFormatContext}, Cint, Ptr{Int64}, Ptr{Int64}), s, stream, dts, wall)
end

function av_hex_dump(f, buf, size::Integer)
    ccall((:av_hex_dump, libavformat), Cvoid, (Ptr{Cvoid}, Ptr{UInt8}, Cint), f, buf, size)
end

function av_hex_dump_log(avcl, level::Integer, buf, size::Integer)
    ccall((:av_hex_dump_log, libavformat), Cvoid, (Ptr{Cvoid}, Cint, Ptr{UInt8}, Cint), avcl, level, buf, size)
end

function av_pkt_dump2(f, pkt, dump_payload::Integer, st)
    ccall((:av_pkt_dump2, libavformat), Cvoid, (Ptr{Cvoid}, Ptr{AVPacket}, Cint, Ptr{AVStream}), f, pkt, dump_payload, st)
end

function av_pkt_dump_log2(avcl, level::Integer, pkt, dump_payload::Integer, st)
    ccall((:av_pkt_dump_log2, libavformat), Cvoid, (Ptr{Cvoid}, Cint, Ptr{AVPacket}, Cint, Ptr{AVStream}), avcl, level, pkt, dump_payload, st)
end

function av_codec_get_id(tags, tag::Integer)
    ccall((:av_codec_get_id, libavformat), AVCodecID, (Ptr{Ptr{AVCodecTag}}, UInt32), tags, tag)
end

function av_codec_get_tag(tags, id::AVCodecID)
    ccall((:av_codec_get_tag, libavformat), UInt32, (Ptr{Ptr{AVCodecTag}}, AVCodecID), tags, id)
end

function av_codec_get_tag2(tags, id::AVCodecID, tag)
    ccall((:av_codec_get_tag2, libavformat), Cint, (Ptr{Ptr{AVCodecTag}}, AVCodecID, Ptr{UInt32}), tags, id, tag)
end

function av_find_default_stream_index(s)
    ccall((:av_find_default_stream_index, libavformat), Cint, (Ptr{AVFormatContext},), s)
end

function av_index_search_timestamp(st, timestamp::Int64, flags::Integer)
    ccall((:av_index_search_timestamp, libavformat), Cint, (Ptr{AVStream}, Int64, Cint), st, timestamp, flags)
end

function av_add_index_entry(st, pos::Int64, timestamp::Int64, size::Integer, distance::Integer, flags::Integer)
    ccall((:av_add_index_entry, libavformat), Cint, (Ptr{AVStream}, Int64, Int64, Cint, Cint, Cint), st, pos, timestamp, size, distance, flags)
end

function av_url_split(proto, proto_size::Integer, authorization, authorization_size::Integer, hostname, hostname_size::Integer, port_ptr, path, path_size::Integer, url)
    ccall((:av_url_split, libavformat), Cvoid, (Cstring, Cint, Cstring, Cint, Cstring, Cint, Ptr{Cint}, Cstring, Cint, Cstring), proto, proto_size, authorization, authorization_size, hostname, hostname_size, port_ptr, path, path_size, url)
end

function av_dump_format(ic, index::Integer, url, is_output::Integer)
    ccall((:av_dump_format, libavformat), Cvoid, (Ptr{AVFormatContext}, Cint, Cstring, Cint), ic, index, url, is_output)
end

function av_get_frame_filename2(buf, buf_size::Integer, path, number::Integer, flags::Integer)
    ccall((:av_get_frame_filename2, libavformat), Cint, (Cstring, Cint, Cstring, Cint, Cint), buf, buf_size, path, number, flags)
end

function av_get_frame_filename(buf, buf_size::Integer, path, number::Integer)
    ccall((:av_get_frame_filename, libavformat), Cint, (Cstring, Cint, Cstring, Cint), buf, buf_size, path, number)
end

function av_filename_number_test(filename)
    ccall((:av_filename_number_test, libavformat), Cint, (Cstring,), filename)
end

function av_sdp_create(ac, n_files::Integer, buf, size::Integer)
    ccall((:av_sdp_create, libavformat), Cint, (Ptr{Ptr{AVFormatContext}}, Cint, Cstring, Cint), ac, n_files, buf, size)
end

function av_match_ext(filename, extensions)
    ccall((:av_match_ext, libavformat), Cint, (Cstring, Cstring), filename, extensions)
end

function avformat_query_codec(ofmt, codec_id::AVCodecID, std_compliance::Integer)
    ccall((:avformat_query_codec, libavformat), Cint, (Ptr{AVOutputFormat}, AVCodecID, Cint), ofmt, codec_id, std_compliance)
end

function avformat_get_riff_video_tags()
    ccall((:avformat_get_riff_video_tags, libavformat), Ptr{AVCodecTag}, ())
end

function avformat_get_riff_audio_tags()
    ccall((:avformat_get_riff_audio_tags, libavformat), Ptr{AVCodecTag}, ())
end

function avformat_get_mov_video_tags()
    ccall((:avformat_get_mov_video_tags, libavformat), Ptr{AVCodecTag}, ())
end

function avformat_get_mov_audio_tags()
    ccall((:avformat_get_mov_audio_tags, libavformat), Ptr{AVCodecTag}, ())
end

function av_guess_sample_aspect_ratio(format, stream, frame)
    ccall((:av_guess_sample_aspect_ratio, libavformat), AVRational, (Ptr{AVFormatContext}, Ptr{AVStream}, Ptr{AVFrame}), format, stream, frame)
end

function av_guess_frame_rate(ctx, stream, frame)
    ccall((:av_guess_frame_rate, libavformat), AVRational, (Ptr{AVFormatContext}, Ptr{AVStream}, Ptr{AVFrame}), ctx, stream, frame)
end

function avformat_match_stream_specifier(s, st, spec)
    ccall((:avformat_match_stream_specifier, libavformat), Cint, (Ptr{AVFormatContext}, Ptr{AVStream}, Cstring), s, st, spec)
end

function avformat_queue_attached_pictures(s)
    ccall((:avformat_queue_attached_pictures, libavformat), Cint, (Ptr{AVFormatContext},), s)
end

function av_apply_bitstream_filters(codec, pkt, bsfc)
    ccall((:av_apply_bitstream_filters, libavformat), Cint, (Ptr{AVCodecContext}, Ptr{AVPacket}, Ptr{AVBitStreamFilterContext}), codec, pkt, bsfc)
end

function avformat_transfer_internal_stream_timing_info(ofmt, ost, ist, copy_tb::AVTimebaseSource)
    ccall((:avformat_transfer_internal_stream_timing_info, libavformat), Cint, (Ptr{AVOutputFormat}, Ptr{AVStream}, Ptr{AVStream}, AVTimebaseSource), ofmt, ost, ist, copy_tb)
end

function av_stream_get_codec_timebase(st)
    ccall((:av_stream_get_codec_timebase, libavformat), AVRational, (Ptr{AVStream},), st)
end


text:
using SortingAlgorithms

function main(points)
  try
    ch_points = []
    ch_points = convex_hull!(points)

    output = [[],[],[],[]]

    for ch_point in ch_points
      push!(output[1], ch_point[1])
      push!(output[2], ch_point[2])

      deleteat!(points, findfirst(points, ch_point))
    end

    for point in points
      push!(output[3], point[1])
      push!(output[4], point[2])
    end

    return output
  catch e
    if isa(e, MethodError)
      return "Please provide array of tuples where each element is Int64"
    else
      return e
    end
  end
end

function convex_hull!(points::Array{Tuple{Int64,Int64},1})
  if size(points, 1) < 3
    return points
  end

  # Sort all points
  points = sort!(union(points), alg = HeapSort)

  # Arrays for top and bottom parts of convex hull
  calculated_ch_top = []
  calculated_ch_bottom = []

  # Build top part of convex hull
  build_convex_hull_part!(calculated_ch_top, points, 1:1:length(points))

  # Build bottom part of convex hull
  build_convex_hull_part!(calculated_ch_bottom, points, length(points):-1:1)

  # Remove duplicated points:
  # End of top part == beginning of the bottom part
  pop!(calculated_ch_top)

  # End of bottom part == beginning of the top part
  pop!(calculated_ch_bottom)

  # Top + Bottom parts == convex hull
  convex_hull = [calculated_ch_top; calculated_ch_bottom]

  return convex_hull
end

function counter_clockwise(p1, p2, p3)
  return ((p2[1] - p1[1]) * (p3[2] - p1[2])) - ((p2[2] - p1[2]) * (p3[1] - p1[1]))
end

function generate_points(number_of_points, range_min, range_max)
  points = [(0,0)]
  pop!(points)
  for i = 1:number_of_points
    generated_point = (rand(range_min:range_max), rand(range_min:range_max))
    push!(points, generated_point)
  end
  return points
end

function build_convex_hull_part!(hull_points, points_to_check, iterator_range)
  for i = iterator_range
    point = points_to_check[i]

    while length(hull_points) > 1 && counter_clockwise(hull_points[end-1], hull_points[end], point) >= 0
      pop!(hull_points)
    end

    push!(hull_points, point)
  end
end

max_stars_repo_name: MrShemek/ConvexHullTDD

max_issues_repo_name: MrShemek/ConvexHullTDD

max_forks_repo_name: MrShemek/ConvexHullTDD


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Htop_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Htop")
JLLWrappers.@generate_main_file("Htop", UUID("3902149a-7cde-5d36-8de2-61bc48a40c47"))
end  # module Htop_jll


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to the general functions using Euler angle and axis.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/angleaxis.jl
# ========================

# Functions: inv
# --------------

@testset "General functions of Euler angle and axis: inv" begin
    # Float64
    # ==========================================================================
    T = Float64
    av = EulerAngleAxis(deg2rad(T(20)), T[sqrt(2) / 2, 0, sqrt(2) / 2])
    iav = inv(av)
    @test eltype(iav) === T
    @test iav.a ≈ deg2rad(T(20))
    @test iav.v ≈ T[-sqrt(2) / 2, 0, -sqrt(2)/2]

    av = EulerAngleAxis(deg2rad(T(-20)), T[sqrt(2) / 2, 0, sqrt(2) / 2])
    iav = inv(av)
    @test eltype(iav) === T
    @test iav.a ≈ deg2rad(T(20))
    @test iav.v ≈ T[sqrt(2) / 2, 0, sqrt(2)/2]

    # Float32
    # ==========================================================================
    T = Float32
    av = EulerAngleAxis(deg2rad(T(20)), T[sqrt(2) / 2, 0, sqrt(2) / 2])
    iav = inv(av)
    @test eltype(iav) === T
    @test iav.a ≈ deg2rad(T(20))
    @test iav.v ≈ T[-sqrt(2) / 2, 0, -sqrt(2)/2]

    av = EulerAngleAxis(deg2rad(T(-20)), T[sqrt(2) / 2, 0, sqrt(2) / 2])
    iav = inv(av)
    @test eltype(iav) === T
    @test iav.a ≈ deg2rad(T(20))
    @test iav.v ≈ T[sqrt(2) / 2, 0, sqrt(2)/2]
end

# Functions: show
# ---------------

@testset "General functions of Euler angle and axis: show" begin
    buf = IOBuffer()
    io = IOContext(buf)
    av = EulerAngleAxis(π / 3, [sqrt(3) / 3, -sqrt(3) / 3, sqrt(3) / 3])

    # Extended printing.
    show(io, MIME"text/plain"(), av)
    expected = """
        EulerAngleAxis{Float64}:
          Euler angle : 1.0472 rad  (60.0°)
          Euler axis  : [0.57735, -0.57735, 0.57735]"""
    @test String(take!(io.io)) == expected

    # Comapct printing.
    show(io, av)
    expected = """
        EulerAngleAxis{Float64}: θ = 1.0472 rad, v = [0.57735, -0.57735, 0.57735]"""
    @test String(take!(io.io)) == expected

    # Colors.
    io = IOContext(buf, :color => true)
    show(io, MIME"text/plain"(), av)
    expected = """
        EulerAngleAxis{Float64}:
        \e[32;1m  Euler angle : \e[0m1.0472 rad  (60.0°)
        \e[33;1m  Euler axis  : \e[0m[0.57735, -0.57735, 0.57735]"""
    @test String(take!(io.io)) == expected
end


text:
# calculation of steady adjoint

using PETSc2

"""
  Use an existing [`LinearSolver`](@ref) to compute the adjoint.
  A preconditioner and linear operator can be obtained from
  [`NonlinearSolver.getNewtonPCandLO`](@ref).

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * ls: a LinearSolver object
   * functionalData: an [`AbstractFunctional`](@ref) to compute the adjoint of

  **Inputs/Outputs**

   * adjoint_vec: Array{Tsol, 1} of length mesh.numDof, (overwritten)

  **Keyword Arguments**

   * recalc_jac: recalc the linear operator inside `ls`, default false
   * recalc_pc: recalc the preconditioner inside `ls`, default false
   * start_comm: do parallel communication before recomputing the pc and/or jac,
                 default false
"""
function calcAdjoint(mesh::AbstractDGMesh{Tmsh},
                  sbp::AbstractOperator, eqn::AbstractSolutionData{Tsol, Tres}, opts,
                  ls::LinearSolver, functionalData::AbstractFunctional,
                  adjoint_vec::Array{Tsol,1}; recalc_jac=false,
                  recalc_pc=false, start_comm=false) where {Tmsh, Tsol, Tres}

  if !opts["need_adjoint"]
    error("""must specify opts["need_adjoint"] for adjoint computation""")
  end

  # recalc operators if requested
  ctx_residual = (evalResidual,)
  if recalc_jac && recalc_pc
    calcPCandLO(ls, mesh, sbp, eqn, opts, ctx_residual, 0.0, start_comm=start_comm)
  elseif recalc_jac
    calcLinearOperator(ls, mesh, sbp, eqn, opts, ctx_residual, 0.0, start_comm=start_comm)
  elseif recalc_pc
    calcPC(ls, mesh, sbp, eqn, opts, ctx_residual, 0.0, start_comm=start_comm)
  end

  #TODO: store this somewhere?
  func_deriv = zeros(Tsol, mesh.numDof)

  # 3D array into which func_deriv_arr gets interpolated
  func_deriv_arr = zeros(eqn.q)

  # Calculate df/dq_bndry on edges where the functional is calculated and put
  # it back in func_deriv_arr
  evalFunctionalDeriv_q(mesh, sbp, eqn, opts, functionalData, func_deriv_arr)

  # Assemble func_deriv
  array3DTo1D(mesh, sbp, eqn, opts, func_deriv_arr, func_deriv)
  scale!(func_deriv, -1.0)

  # do transpose solve
  _adjoint_vec = zeros(real(Tsol), length(adjoint_vec))
  linearSolveTranspose(ls, real(func_deriv), _adjoint_vec)
  copy!(adjoint_vec, _adjoint_vec)

  # Output/Visualization options for Adjoint
  if opts["write_adjoint"]
    outname = string("adjoint_vec_", mesh.myrank,".dat")
    f = open(outname, "w")
    for i = 1:length(adjoint_vec)
      println(f, real(adjoint_vec[i]))
    end
    close(f)
  end

  if opts["write_adjoint_vis"]
    saveSolutionToMesh(mesh, adjoint_vec)
    fname = "adjoint_field"
    writeVisFiles(mesh, fname)
  end

  return nothing
end

"""
  Convenience method for solving the adjoint is you don't already have a
  [`LinearSolver`](@ref).  It creates the default one used for Newton's method.
  Note that this can be bad if globalization is used because the globalization
  might be added to the linear operator used for the adjoint solve

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * functionalData: functional to compute the adjoint for

  **Inputs/Outputs**

   * adjoint_vec: vector, same shape and element type as `eqn.q_vec`, to be
                  overwritten with the adjoint

  **Keyword Arguments**

   * start_comm: if true, start parallel communication before computing the
                 linear operator for the adjoint solve, default true
"""
function calcAdjoint(mesh::AbstractDGMesh{Tmsh},
                  sbp::AbstractOperator, eqn::AbstractSolutionData{Tsol, Tres}, opts,
                  functionalData::AbstractFunctional,
                  adjoint_vec::Array{Tsol,1}; start_comm=true) where {Tmsh, Tsol, Tres}

  ls = createLinearSolver(mesh, sbp, eqn, opts)
  calcAdjoint(mesh, sbp, eqn, opts, ls, functionalData, adjoint_vec;
              recalc_pc=true, recalc_jac=true, start_comm=start_comm)

  return nothing
end


text:
#####
##### Model
#####
crc3, parameters, environment, alphas, stepsymbol, cl_stall, nx, nu = MT.build_crc3(; nwings)
crc3_model = MT.AircraftDynamics(crc3, parameters, environment, stepsymbol, cl_stall, nx, nu)

#####
##### re-nondimensionalize
#####
function alpha_2_CF_CM(crc3_model, alpha; vinf = 1.0, rotors_on = false, U = [0.0, 0.0]* 2*pi/60.0, frame = MT.Wind())
    # build state
    X = MT.level_flight_state(alpha, vinf)

    # get forces/moments
    force, M_y = MT.aerodynamics!(crc3_model, X, U; rotors_on)

    # get normalization parameters
    Sref = crc3_model.aircraft.wing_system.system.reference[1].S
    b = crc3_model.aircraft.wing_system.system.reference[1].b
    c = crc3_model.aircraft.wing_system.system.reference[1].c
    rho = crc3_model.environment.ρ
    qinf = 0.5 * rho * vinf^2

    # get force/moment coefficients
    CF = force ./ qinf ./ Sref
    CM = M_y / qinf / Sref / c

    return CF, CM
end

#####
##### test!
#####
function test_cl_alpha(crc3_model; alphas = range(-20.0, stop=25.0, length=26) .* pi/180, vinf = 1.0, rotors_on = false, U = [0.0, 0.0], frame = MT.Wind(), save_figures = true, clear_figure = true)
    # initialize vars
    CFs = zeros(2,length(alphas))
    CMs = zeros(length(alphas))

    # build lift curves
    for (i,alpha) in enumerate(alphas)
        CF, CM = alpha_2_CF_CM(crc3_model, alpha; vinf, rotors_on, U, frame)
        CFs[:,i] .= CF
        CMs[i] = CM
    end

    fig_cf = PP.figure("test_CF")
    if clear_figure; fig_cf.clear(); end
    fig_cf.add_subplot(2, 1, 1, ylabel = MT.LS.L"C_L")
    fig_cf.add_subplot(2, 1, 2, ylabel = MT.LS.L"C_D", xlabel = MT.LS.L"\alpha[^\circ]")
    axs_cf = fig_cf.get_axes()
    axs_cf[2].plot(alphas .* 180/pi, CFs[1,:])
    axs_cf[1].plot(alphas .* 180/pi, CFs[2,:])

    fig_cm = PP.figure("test_CM")
    if clear_figure; fig_cm.clear(); end
    fig_cm.add_subplot(1,1,1, ylabel = MT.LS.L"C_M", xlabel = MT.LS.L"\alpha[^\circ]")
    axs_cm = fig_cm.get_axes()[1]
    axs_cm.plot(alphas .* 180/pi, CMs)

    if save_figures
        fig_cf.savefig(joinpath(MT.topdirectory, "data", "plots", MT.TODAY, "test_aerodynamics_cf_$(nwings)_wing.png"))
        fig_cm.savefig(joinpath(MT.topdirectory, "data", "plots", MT.TODAY, "test_aerodynamics_cm_$(nwings)_wing.png"))
    end

    return alphas, CFs, CMs
end

_ = test_cl_alpha(crc3_model; U=[0.0])

max_stars_repo_name: byuflowlab/MultifidelityTrajectories

max_issues_repo_name: byuflowlab/MultifidelityTrajectories

max_forks_repo_name: byuflowlab/MultifidelityTrajectories


text:
using Ramnet.Models: train!, predict
using Ramnet.Models.Nodes

using Ramnet.Utils: stack

@testset "Nodes" begin
    # Standard Node
    all_active   = ones(Bool, 9)
    all_inactive = zeros(Bool, 9)

    std_n = DictNode()

    train!(std_n, all_active)

    @test predict(std_n, all_active) == 1
    @test predict(std_n, all_inactive) == 0
    @test predict(std_n, stack(all_active, all_inactive)) == [1, 0]

    # Accumulator Node
    acc_n = AccNode()

    train!(acc_n, all_active)

    @test predict(acc_n, all_active) == true
    @test predict(acc_n, all_active; b=1) == false

    train!(acc_n, all_active)

    @test predict(acc_n, all_active) == true
    @test predict(acc_n, all_active; b=1) == true

    train!(acc_n, all_inactive)

    @test predict(acc_n, stack(all_active, all_inactive)) == Bool[1, 1]
    @test predict(acc_n, stack(all_active, all_inactive); b=1) == Bool[1, 0]
end

max_stars_repo_name: RafaelFK/ramnet

max_issues_repo_name: RafaelFK/ramnet

max_forks_repo_name: RafaelFK/ramnet


text:
using BlackBoxOptim: general_stop_reason

@testset "general_stop_reason" begin
    problem = BlackBoxOptim.example_problems["Sphere"]

    @testset "Within fitness tolerance" begin
        res = bboptimize(problem; NumDimensions = 1, TraceMode = :silent,
            FitnessTolerance = 1e-1)
        @test general_stop_reason(res) == "Within fitness tolerance of optimum"
    end

    @testset "Delta fitness below tolerance" begin
        res = bboptimize(problem; NumDimensions = 1, TraceMode = :silent,
            MinDeltaFitnessTolerance = 1e-1, FitnessTolerance = 1e-10)
        @test general_stop_reason(res) == "Delta fitness below tolerance"
    end
end


text:
struct SilvisModel{T}<:AbstractLESModel
    c::T
    cp::T
    Δ²::T
    tau::SymTrTenField{T,3,2,false}
    reduction::Vector{T}
end

function SilvisModel(c::T,cp::T,Δ::Real,dim::NTuple{3,Integer}) where {T<:Real}
    data = SymTrTenField{T}(dim,(LX,LY,LZ))
    #fill!(data,0)
    reduction = zeros(THR ? Threads.nthreads() : 1)
    return SilvisModel{T}(c,cp,Δ^2, data,reduction)
end

is_Silvis(a::Union{<:SilvisModel,Type{<:SilvisModel}}) = true
@inline @par is_Silvis(s::Type{T}) where {T<:@par(AbstractSimulation)} = is_Silvis(LESModelType)
@inline is_Silvis(s::T) where {T<:AbstractSimulation} = is_Silvis(T)
is_Silvis(a) = false

msg(a::SilvisModel) = "\nLES model: Silvis\nS Constant: $(a.c)\nP Constant: $(a.cp)\nFilter Width: $(sqrt(a.Δ²))\n"

fvs(S,w) = norm(S⋅w)/(fsqrt(S:S)*norm(w))

@inline function Silvis_eddy_viscosity(S::SymTen,f::Real,c::Real,Δ2::Real)
    return c*c*Δ2*f*f*f*fsqrt(S:S)/2
end

@inline function Silvis_P_coeff(cp::Real,f::Real,Δ2::Real)
    return cp*Δ2*f*f*f*f/4
end


text:
import Compat: Libc, @compat

function test_numerics()
    PostgresType = PostgreSQL.PostgresType
    values = @compat Any[Int16(4), Int32(4), Int64(4), Float32(4), Float64(4)]

    types = [PostgresType{:int2}, PostgresType{:int4}, PostgresType{:int8},
        PostgresType{:float4}, PostgresType{:float8}]

    p = convert(Ptr{UInt8}, Libc.malloc(8))
    try
        for i = 1:length(values)
            p = PostgreSQL.storestring!(p, string(values[i]))
            data = PostgreSQL.jldata(types[i], p)
            testsameness(data, values[i])
        end
    finally
        Libc.free(p)
    end

    p = Libc.malloc(8)
    try
        for i = 1:length(values)
            p = PostgreSQL.pgdata(types[i], convert(Ptr{UInt8}, p), values[i])
            data = PostgreSQL.jldata(types[i], p)
            testsameness(data, values[i])
        end
    finally
        Libc.free(p)
    end
end

function test_strings()
    PGType = PostgreSQL.PostgresType
    for typ in [PGType{:varchar}, PGType{:text}, PGType{:bpchar}]
        for str in Any["foobar", "fooba\u211D"]
            p = PostgreSQL.pgdata(typ, convert(Ptr{UInt8}, C_NULL), str)
            try
                data = PostgreSQL.jldata(typ, p)
                @test typeof(str) == typeof(data)
                @test str == data
            finally
                Libc.free(p)
            end
        end
    end
end

function test_bytea()
    typ = PostgreSQL.PostgresType{:bytea}
    bin = (UInt8)[0x01, 0x03, 0x42, 0xab, 0xff]
    p = PostgreSQL.pgdata(typ, convert(Ptr{UInt8}, C_NULL), bin)
    try
        data = PostgreSQL.jldata(typ, p)
        @test typeof(bin) == typeof(data)
        @test bin == data
    finally
        Libc.free(p)
    end
end

function test_json()
    PGType = PostgreSQL.PostgresType
    for typ in Any[PGType{:json}, PGType{:jsonb}]
        dict1 = @compat Dict("bobr dobr" => [1, 2, 3], "foo" => 3.0)
        p = PostgreSQL.pgdata(typ, convert(Ptr{UInt8}, C_NULL), dict1)
        try
            dict2 = PostgreSQL.jldata(typ, p)
            @test dict1 == dict2
        finally
            Libc.free(p)
        end
    end
end

test_numerics()
test_strings()
test_bytea()
test_json()

max_stars_repo_name: beffaCo/updatePostgresql

max_forks_repo_name: beffaCo/updatePostgresql


text:
using ARDESPOT, ParticleFilters, BeliefUpdaters, POMDPSimulators
using POMCPOW, MCTS, POMDPPolicies

include("pomdps.jl")
include("visibility_polygons.jl")

# This function overrides an ARDESPOT default behavior to
# make state-dependent actions spaces. Native support for
# that funcationality does not exist yet apparently.
actions(p::SkyNet, b::ARDESPOT.ScenarioBelief) = actions(p, b.scenarios[1][2])
function actions(p::SkyNet, b::POMCPOW.POWTreeObsNode)
    @show belief(b)
    actions(p, rand(belief(b)))
end
function MCTS.next_action(gen::RandomActionGenerator, p::SkyNet, s, snode::AbstractStateNode)
    rand(actions(p, rand(s)))
end
function POMDPPolicies.action(policy::RandomPolicy, b::Nothing)
    return rand(policy.rng, actions(policy.problem))
end

n_parts = 10000

if !isdefined(Main, :p) || !(p isa SkyNet)

    @eval begin
        R = 3
        n = 2
        println("Building SkyNet with R = $R, n = $n")
        n_parts_prev = n_parts
        p = SkyNet(R = R,
                   n_agents = n,
                   partition_points =[CartesianIndex((2, 3)),
                                       CartesianIndex((5,13)),
                                       # CartesianIndex((13,7))
                                       ])

        println("initializing everything else with $n_parts particles")
        ag0 = Agent.(p.partition_points, p.partitions)

        b0 = [SNState(ag0, [rand_pos(graph(p))]) for i in 1:n_parts]
        s0 = rand(b0)

        ## MAYBE MAKE OUR OWN
        belief_updater = SIRParticleFilter(p, n_parts)
    end
end

K = 1000
D = 100
solver = DESPOTSolver(T_max = 10.0,
                      K = K,
                      # max_trials = 10000,
                      # tree_in_info = true,
                      D = D,
                      bounds = IndependentBounds(-1e6, 1e6, check_terminal = true),
                      random_source = MemorizingSource(K, D, Random.GLOBAL_RNG, min_reserve=1000))
# solver = POMCPOWSolver() # has hecka keywords
policy = solve(solver, p)

println("simulating...")
history = simulate(HistoryRecorder(max_steps = 2), p, policy, belief_updater, b0, s0)

for (s, a, r) in zip(history.state_hist, history.action_hist, history.reward_hist)
    println("Adversary at:\t", adversaries(s)[1])
    println("Action taken:\t", a)
    println("Total reward:\t", r)
end

# for (s, a, o, r) in stepthrough(p, policy, "saor", max_steps = 2)
    # println("Adversary at:\t", adversaries(s))
    # println("Action taken:\t", a)
    # println("Total reward:\t", r)
# end


text:
export cudaDataType, libraryPropertyType

@enum libraryPropertyType::Cint begin
    MAJOR_VERSION = 0
    MINOR_VERSION
    PATCH_LEVEL
end

@enum cudaDataType::Cint begin
    R_16F = 2  # `CUDA_R_16F`, real as a half
    C_16F = 6  # `CUDA_C_16F`, complex as a pair of half numbers
    R_32F = 0  # `CUDA_R_32F`, real as a float
    C_32F = 4  # `CUDA_C_32F`, complex as a pair of float numbers
    R_64F = 1  # `CUDA_R_64F`, real as a double
    C_64F = 5  # `CUDA_C_64F`, complex as a pair of double numbers
    R_8I  = 3  # `CUDA_R_8I`,  real as a signed char
    C_8I  = 7  # `CUDA_C_8I`,  complex as a pair of signed char numbers
    R_8U       # `CUDA_R_8U`,  real as a unsigned char
    C_8U       # `CUDA_C_8U`,  complex as a pair of unsigned char numbers
    R_32I      # `CUDA_R_32I`, real as a signed int
    C_32I      # `CUDA_C_32I`, complex as a pair of signed int numbers
    R_32U      # `CUDA_R_32U`, real as a unsigned int
    C_32U      # `CUDA_C_32U`, complex as a pair of unsigned int numbers
end

function cudaDataType(T::DataType)
    if T == Float32
        return R_32F
    elseif T == ComplexF32
        return C_32F
    elseif T == Float16
        return R_16F
    elseif T == ComplexF16
        return C_16F
    elseif T == Float64
        return R_64F
    elseif T == ComplexF64
        return C_64F
    elseif T == Int8
        return R_8I
    elseif T == Complex{Int8}
        return C_8I
    elseif T == Int32
        return R_32I
    elseif T == Complex{Int32}
        return C_32I
    elseif T == UInt8
        return R_8U
    elseif T == Complex{UInt8}
        return C_8U
    elseif T == UInt32
        return R_32U
    elseif T == Complex{UInt32}
        return C_32U
    else
        throw(ArgumentError("cudaDataType equivalent for input type $T does not exist!"))
    end
end


text:
using Animations
using Colors
using Javis
using NBodySimulator
using StaticArrays

m1 = MassBody(SVector(-0.995492, 0.0, 0.0), SVector(-0.347902, -0.53393, 0.0), 1.0)
m2 = MassBody(SVector(0.995492, 0.0, 0.0), SVector(-0.347902, -0.53393, 0.0), 1.0)
m3 = MassBody(SVector(0.0, 0.0, 0.0), SVector(0.695804, 1.067860, 0.0), 1.0)

G = 1
system = GravitationalSystem([m1, m2, m3], G)

tspan = (0.0, 2pi);
simulation = NBodySimulation(system, tspan)

sim_result = run_simulation(simulation,saveat = 0.01)
sim_data = sim_result.solution.u

# extract the data
planet1_pos = [(sim_data[x][1, 1], sim_data[x][2, 1]) for x in 1:length(sim_data)]
planet2_pos = [(sim_data[x][1, 2], sim_data[x][2, 2]) for x in 1:length(sim_data)]
planet3_pos = [(sim_data[x][1, 3], sim_data[x][2, 3]) for x in 1:length(sim_data)]

# scaling factor
scaling_factor = 210
#Convert to Points
planet1_point = Point.((scaling_factor .* getindex.(planet1_pos, 1)), ((scaling_factor .* getindex.(planet1_pos, 2))))
planet2_point = Point.((scaling_factor .* getindex.(planet2_pos, 1)), ((scaling_factor .* getindex.(planet2_pos, 2))))
planet3_point = Point.((scaling_factor .* getindex.(planet3_pos, 1)), ((scaling_factor .* getindex.(planet3_pos, 2))))

#Julia logo colors
logocolors = Colors.JULIA_LOGO_COLORS
red = logocolors.red
green = logocolors.green
purple = logocolors.purple
blue = logocolors.blue

function ground(args...)
    background("black")
    sethue(blue)
end

function planet(p, r, color)
    sethue(color)
    return circle(p, r, :fill)
end

function make_animation(red, green, purple)
    nframes = 200

    vid = Video(500, 500)
    Background(1:nframes, ground)

    planet1_anim = Animation(
        collect(range(0, stop = 1, length = length(planet1_point))),
        planet1_point,
        linear()
    )

    planet2_anim = Animation(
        collect(range(0, stop = 1, length = length(planet2_point))),
        planet2_point,
        linear()
    )

    planet3_anim = Animation(
        collect(range(0, stop = 1, length = length(planet3_point))),
        planet3_point,
        linear()
    )

    planet1 = Object(1:nframes, (args...) -> planet(O, 20.0, red))
    planet2 = Object(1:nframes, (args...) -> planet(O, 20.0, purple))
    planet3 = Object(1:nframes, (args...) -> planet(O, 20.0, green))

    act!(planet1, Action(1:nframes, planet1_anim, translate()))
    act!(planet2, Action(1:nframes, planet2_anim, translate()))
    act!(planet3, Action(1:nframes, planet3_anim, translate()))

    render(vid; pathname = "three_body.gif")
end

make_animation(red, green, purple)

max_stars_repo_name: Ved-Mahajan/graphic-dev

max_issues_repo_name: Ved-Mahajan/graphic-dev

max_forks_repo_name: Ved-Mahajan/graphic-dev


text:
for (i,f) ∈ enumerate((:second,:third,:fourth,:fifth,:sixth,:seventh,:eighth,:ninth,:tenth,:eleventh,:twelfth,:thirteenth))
  @eval @inline $f(x) = @inbounds getindex(x, $(i+1))
end


text:
"""
    write_mtg(file, mtg; kwargs...)
    write_mtg(file, mtg, classes, description, features)

Write an mtg file to disk.

# Arguments

- `file::String`: The path to the MTG file to write.
- `mtg`: the mtg
- `classes`: the classes section
- `description`: the description section
- `features`: the features section

# Note

kwargs can be used to give zero, one or two of the classes, description and features
instead of all. In this case the missing ones are recomputed using [`get_classes`](@ref),
[`get_features`](@ref) or [`get_description`](@ref).

# Examples

```julia
file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","simple_plant.mtg")
mtg = read_mtg(file)
write_mtg("test.mtg",mtg)
```
"""
function write_mtg(file, mtg; kwargs...)
    kwargs = (; kwargs...)

    if !haskey(kwargs, :classes)
        classes = get_classes(mtg)
    end

    if !haskey(kwargs, :description)
        description = nothing
    end

    if !haskey(kwargs, :features)
        features = get_features(mtg)
    end

    write_mtg(file, mtg, classes, description, features)
end

function write_mtg(file, mtg, classes, description, features)
    @info "Writing mtg to $file"
    open(file, "w") do io
        # Code section:
        writedlm(io, ["CODE:" "FORM-A"])

        # Classes section:
        writedlm(io, [""])
        writedlm(io, ["CLASSES:"])
        writedlm(io, reshape(names(classes), (1, :)))
        writedlm(io, eachrow(classes))

        # Description section:
        writedlm(io, [""])
        writedlm(io, ["DESCRIPTION:"])

        # Description is optional
        if description !== nothing
            # Reformat the RIGHT column to match how it is written in an MTG
            right = fill("", size(description)[1])

            for i in 1:length(description.RIGHT)
                right = join(description.RIGHT[i], ",")
            end
            description[!, :RIGHT] .= right

            writedlm(io, reshape(names(description), (1, :)))

            writedlm(io, eachrow(description))
        else
            writedlm(io, ["LEFT" "RIGHT" "RELTYPE" "MAX"])
        end

        # Features section:
        writedlm(io, [""])
        writedlm(io, ["FEATURES:"])
        writedlm(io, reshape(names(features), (1, :)))
        writedlm(io, eachrow(features))

        # MTG section:
        writedlm(io, [""])
        writedlm(io, ["MTG:"])
        mtg_df, mtg_colnames = paste_node_mtg(mtg, features)
        writedlm(io, reshape(mtg_colnames, (1, :)), quotes=false)
        for i in 1:length(mtg_df["mtg_print"])
            writedlm(io, reshape([mtg_df[k][i] for k in keys(mtg_df)], (1, :)), quotes=false)
        end
    end
end

function paste_node_mtg(mtg, features)
    @mutate_mtg!(
        mtg,
        lead = MultiScaleTreeGraph.get_leading_tabs(node),
        mtg_print = MultiScaleTreeGraph.paste_mtg_node(node),
        mtg_refer = MultiScaleTreeGraph.get_reference(node)
    )

    attributes = Dict_attrs(mtg, ["mtg_print", features.NAME..., "lead", "mtg_refer"])
    max_tabs = maximum(attributes["lead"])

    # Build the "ENTITY-CODE" column with necessary "^", leading and trailing tabs
    attributes["mtg_print"] = string.(
        # Add the leading tabulations:
        repeat.("\t", attributes["lead"]),
        # Add the "^" keyword before mtg print in case we refer to the column above:
        attributes["mtg_refer"],
        # Add the mtg printing (e.g. "/Axis0"):
        attributes["mtg_print"],
        # Add the trailing tabulations:
        repeat.("\t", max_tabs .- attributes["lead"])
    )
    # Remove the lead and mtg_refer columns now that we used it:
    pop!(attributes, "lead")
    pop!(attributes, "mtg_refer")

    # Replacing all nothing values by tabulations:
    for (key, val) in attributes
        replace!(val, nothing => "")
    end
    # Renaming first column and adding tabs:
    mtg_colnames = collect(keys(attributes))
    mtg_colnames[1] = string("ENTITY-CODE", repeat("\t", max_tabs))
    return attributes, mtg_colnames
end

"""
    paste_mtg_node(node)

Parse the mtg node as it should appear in the mtg file.
"""
function paste_mtg_node(node)
    index = node.MTG.index === nothing ? "" : string(node.MTG.index)
    node.MTG.link * node.MTG.symbol * index
end

"""
    get_leading_tabs(node)

Get the number of tabulation the node should have when writting it to a file based on the
topology of its parent.
"""
function get_leading_tabs(node)
    if isroot(node)
        return 0
    else
        node.MTG.link == "+" ? node.parent[:lead] + 1 : node.parent[:lead]
    end
end

"""
    get_reference(node)

Get the preceding "^" keyword if needed, *i.e.* in case we refer to the parent node in the
same mtg file column.
"""
function get_reference(node)
    if isroot(node)
        return ""
    else
        node.MTG.link == "+" ? "" : "^"
    end
end

function Dict_attrs(mtg, attrs)
    df = OrderedDict()
    for var in attrs
        push!(df, var => descendants(mtg, var, self=true))
    end
    df
end


text:
function Base.sort!(dt::DataTable; cols=Any[], alg=nothing,
                    lt=isless, by=identity, rev=false, order=Forward)
    if !(isa(by, Function) || eltype(by) <: Function)
        msg = "'by' must be a Function or a vector of Functions. Perhaps you wanted 'cols'."
        throw(ArgumentError(msg))
    end
    ord = ordering(dt, cols, lt, by, rev, order)
    _alg = Sort.defalg(dt, ord; alg=alg, cols=cols)
    sort!(dt, _alg, ord)
end

function Base.sort!(dt::DataTable, a::Base.Sort.Algorithm, o::Base.Sort.Ordering)
    p = sortperm(dt, a, o)
    pp = similar(p)
    c = columns(dt)

    for (i,col) in enumerate(c)
        # Check if this column has been sorted already
        if any(j -> c[j]===col, 1:i-1)
            continue
        end

        copy!(pp,p)
        Base.permute!!(col, pp)
    end
    dt
end


text:
"""
An abstract type to define custom spacing distributions.
"""
abstract type AbstractSpacing end

struct Cosine  <: AbstractSpacing end
struct Sine    <: AbstractSpacing end
struct Uniform <: AbstractSpacing end

spacing(n, :: Uniform, flip = false) = uniform_spacing(0., 1., n + 1)
spacing(n, :: Cosine, flip = false)  = cosine_spacing(0.5, 1., n + 1)
spacing(n, :: Sine, flip = false) = sine_spacing(0., 1., (n + 1) * ifelse(flip, -1, 1))
spacing(n, distribution :: AbstractSpacing = Cosine(), flip = false) = spacing(n, distribution, flip)

Base.reverse(q :: AbstractSpacing) = q
Base.Broadcast.broadcastable(q::AbstractSpacing) = Ref(q)

function chop_sections(set1, set2, n :: Integer, distribution = Cosine(); flip = false)
    space = spacing(n, distribution, flip)
    @views [ weighted_vector.(set1, set2, μ) for μ ∈ space ][1:end-1]
end

chop_coordinates(coords, n, spacing = Cosine(), flip = false) = @views [ reduce(vcat, chop_sections.(coords[1:end-1], coords[2:end], n, spacing; flip = flip)); [ coords[end] ] ]

chop_spans(xyzs, div, spacing = Cosine(), flip = false) = permutedims(combinedimsview(map(xyz -> chop_coordinates(xyz, div, spacing, flip), eachrow(xyzs))))

chop_chords(xyzs, div, spacing = Cosine()) = combinedimsview(map(xyz -> chop_coordinates(xyz, div, spacing), eachcol(xyzs)))

chop_wing(xyzs, span_num, chord_num; span_spacing = Cosine(), chord_spacing = Cosine(), flip = false) = chop_chords(chop_spans(xyzs, span_num, span_spacing, flip), chord_num, chord_spacing)

transform_coordinates(xyz, twist, section) = eachrow(xyz * RotY(-twist)') .+ Ref(section)

function chop_spanwise_sections(scaled_foils, twisties, leading_xyz, span_num, spacings, flip = false)
    # Reverse direction if left side
    if flip
        scaled_foils = reverse(scaled_foils)
        twisties     = reverse(twisties)
    end

    # Rotate and translate coordinates
    foil_coords = combinedimsview(transform_coordinates.(scaled_foils, twisties, leading_xyz))

    # Chop up spanwise sections
    chop_spans(foil_coords, span_num, spacings, flip)
end


text:
# using Pkg; for p in ("Knet","ArgParse"); haskey(Pkg.installed(),p) || Pkg.add(p); end

"""

julia bilstm-tagger-withchar.jl # to use with default options on CPU
julia bilstm-tagger-withchar.jl --usegpu # to use with default options on GPU
julia bilstm-tagger-withchar.jl -h # to see all options with default values

This example implements a named entity tagger built on top of a BiLSTM
neural network similar to the model defined in 'Bidirectional LSTM-CRF Models
for Sequence Tagging', Zhiheng Huang, Wei Xu, Kai Yu, arXiv technical report
1508.01991, 2015. Originally, this model implemented for dynet-benchmarks.
This model generates UNK words' embeddings by using an another BiLSTM network
which takes characters as input.

* Paper url: https://arxiv.org/pdf/1508.01991.pdf
* DyNet report: https://arxiv.org/abs/1701.03980
* Benchmark repo: https://github.com/neulab/dynet-benchmark

"""
module CharTagger
using Knet, CUDA, ArgParse, Dates, Random, Printf

include(Knet.dir("data","wikiner.jl"))
const F = Float32
t00 = now()

function main(args)
    s = ArgParseSettings()
    s.description = "Bidirectional LSTM Tagger (with chars) in Knet"

    @add_arg_table s begin
        ("--usegpu"; action=:store_true; help="use GPU or not")
        ("--cembed"; arg_type=Int; default=20; help="char embedding size")
        ("--wembed"; arg_type=Int; default=128; help="word embedding size")
        ("--hidden"; arg_type=Int; default=50; help="BiLSTM hidden size")
        ("--mlp"; arg_type=Int; default=32; help="MLP size")
        ("--timeout"; arg_type=Int; help="max timeout (in seconds)"; default=600)
        ("--epochs"; arg_type=Int; default=100; help="epochs")
        ("--minoccur"; arg_type=Int; default=6; help="word min occurence limit")
        ("--report"; arg_type=Int; default=500; help="report period in iters")
        ("--valid"; arg_type=Int; default=10000; help="valid period in iters")
        ("--seed"; arg_type=Int; default=-1; help="random seed")
    end

    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)
    o[:seed] > 0 && Knet.setseed(o[:seed])
    atype = o[:atype] = !o[:usegpu] ? Array{Float32} : KnetArray{Float32}
    datadir = abspath(joinpath(@__DIR__, "../data/tags"))
    datadir = isdir(datadir) ? datadir : WIKINER_DIR

    # load WikiNER data
    data = WikiNERData(datadir, o[:minoccur])

    # build model
    w, srnns = initweights(
        atype, o[:hidden], length(data.w2i), data.ntags, data.nchars,
        o[:wembed], o[:cembed], o[:mlp], o[:usegpu])
    opt = optimizers(w, Adam)

    # train bilstm tagger
    nwords = data.nwords; ntags = data.ntags; nchars = data.nchars
    println("nwords=$nwords, ntags=$ntags, nchars=$nchars"); flush(stdout)
    println("startup time: ", Int((now()-t00).value)*0.001); flush(stdout)
    t0 = now()
    all_time = dev_time = all_tagged = this_tagged = this_loss = 0
    o[:timeout] = o[:timeout] <= 0 ? Inf : o[:timeout]
    for epoch = 1:o[:epochs]
        shuffle!(data.trn)
        for k = 1:length(data.trn)
            iter = (epoch-1)*length(data.trn) + k
            if o[:report] > 0 && iter % o[:report] == 0
                @printf("%f\n", this_loss/this_tagged); flush(stdout)
                all_tagged += this_tagged
                this_loss = this_tagged = 0
                all_time = Int((now()-t0).value)*0.001
            end

            if all_time > o[:timeout] || o[:valid] > 0 && iter % o[:valid] == 0
                dev_start = now()
                good_sent = bad_sent = good = bad = 0.0
                for sent in data.dev
                    input = make_input(sent, data.w2i, data.c2i)
                    nwords = length(sent)
                    ypred,_ = predict(w, input, srnns)
                    ypred = map(
                        x->data.i2t[x], mapslices(argmax,Array(ypred),dims=1))
                    ygold = map(x -> x[2], sent)
                    same = true
                    for (y1,y2) in zip(ypred, ygold)
                        if y1 == y2
                            good += 1
                        else
                            bad += 1
                            same = false
                        end
                    end
                    if same
                        good_sent += 1
                    else
                        bad_sent += 1
                    end
                end
                dev_time += Int((now()-dev_start).value)*0.001
                train_time = Int((now()-t0).value)*0.001-dev_time

                @printf(
                    "tag_acc=%.4f, sent_acc=%.4f, time=%.4f, word_per_sec=%.4f\n",
                    good/(good+bad), good_sent/(good_sent+bad_sent), train_time,
                    all_tagged/train_time); flush(stdout)
                all_time > o[:timeout] && return
            end

            # train with instance
            input = make_input(data.trn[k],data.w2i,data.c2i)
            out = make_output(data.trn[k],data.t2i)
            batch_loss = train!(w,input,out,srnns,opt)
            this_loss += batch_loss
            this_tagged += length(data.trn[k])
        end
        @printf("epoch %d finished\n", epoch-1); flush(stdout)
    end
end

function make_input(sample, w2i, c2i)
    seq, is_word = Any[], Bool[]
    words = map(x->x[1], sample)
    for word in words
        push!(is_word, haskey(w2i, word) && word != UNK)
        if is_word[end]
            push!(seq, w2i[word])
        else
            chars = [PAD; codeunits(word); PAD]
            push!(seq, convert(Array{Int32}, map(c->c2i[c], chars)))
        end
    end

    # construct rare word block - it's for efficiency
    rare_words = filter(x->isa(x,Array), seq)
    sort!(rare_words, by=x->length(x), rev=true)
    longest = length(rare_words) != 0 ? length(rare_words[1]) : 0
    rblock = Int32[]
    batchsizes = zeros(Int32, longest)
    for t = 1:longest
        for i = 1:length(rare_words)
            length(rare_words[i]) < t && break
            push!(rblock, rare_words[i][t])
            batchsizes[t] += 1
        end
    end

    cinds = findall(is_word)
    rinds = findall(.!is_word)
    cwords = seq[is_word]
    cwords = reshape(cwords, 1, length(cwords))
    cwords = convert(Array{Int32}, cwords)
    rwords = rblock

    return (cwords, cinds, rwords, rinds, vec(batchsizes))
end

function make_output(sample,t2i)
    map(s->t2i[s[2]], sample)
end

function initweights(
    atype, hidden, words, tags, chars, wembed, cembed, mlp, usegpu, winit=0.01)
    w = Array{Any}(undef,8)
    _birnninit(x,y) = rnninit(x,y; bidirectional=true, usegpu=usegpu)

    # init rnns
    srnn1, wrnn1 = _birnninit(wembed, hidden)
    w[1] = wrnn1
    srnn2, wrnn2 = _birnninit(cembed, div(wembed,2))
    w[2] = wrnn2

    # weight/bias params for MLP network
    w[3] = convert(atype, winit*randn(mlp, 2*hidden))
    w[4] = convert(atype, zeros(mlp, 1))
    w[5] = convert(atype, winit*randn(tags, mlp))
    w[6] = convert(atype, winit*randn(tags, 1))

    # word/char embeddings
    w[7] = convert(atype, winit*randn(wembed, words))
    w[8] = convert(atype, winit*randn(cembed, chars))
    return w, [srnn1, srnn2]
end

# loss function
function loss(w, input, ygold, srnns)
    py, _ = predict(w,input,srnns)
    return nll(py,ygold)
end

lossgradient = gradloss(loss)

function predict(w,input,srnns)
    x = encoder(w,input,srnns[2])
    x = reshape(x, size(x,1), 1, size(x,2))
    r = srnns[1]; wr = w[1]
    wmlp, bmlp = w[3], w[4]
    wy, by = w[5], w[6]
    y, hy, cy = rnnforw(r,wr,x)
    y2 = reshape(y,size(y,1),size(y,2)*size(y,3))
    y3 = wmlp * y2 .+ bmlp
    return wy*y3.+by, hy, cy
end

# encoder - it generates embeddings
function encoder(w,input,srnn)
    # expand input tuple
    cwords, cinds, rwords, rinds, bs = input

    # common words' embedding
    cembed = w[end-1][:,cwords]
    cembed = reshape(cembed, size(cembed,1), size(cembed)[end])
    length(rinds) == 0 && return cembed

    r = srnn; wr = w[2]
    c0 = w[end][:,rwords]
    y, hy, cy = rnnforw(r,wr,c0; hy=true, cy=true, batchSizes=bs)
    r0 = permutedims(hy, (3,1,2))
    rembed = reshape(r0, size(r0,1)*size(r0,2), size(r0,3))

    e0 = hcat(cembed, rembed)
    e1 = e0[:,[cinds...,rinds...]]
    return e1
end

function train!(w,input,ygold,srnns,opt)
    gloss, lossval = lossgradient(w, input, ygold, srnns)
    for k = 1:length(w)
        gloss[k] != nothing && update!(w[k], gloss[k], opt[k])
    end
    return lossval*(length(input[2])+length(input[4]))
end

splitdir(PROGRAM_FILE)[end] == "bilstm-tagger-withchar.jl" && main(ARGS)
end # module


text:
const SDL = ShaleDrillingLikelihood.SDLParameters

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing

    revs = (
        DrillingRevenue{Unconstrained, NoTrend,       GathProcess, Learn, WithRoyalty},
        DrillingRevenue{Unconstrained, TimeTrend,     GathProcess, Learn, WithRoyalty},
        DrillingRevenue{Unconstrained, TimeFE,        GathProcess, Learn, WithRoyalty},
        DrillingRevenue{Unconstrained, TimeTrendDgt0, GathProcess, Learn, WithRoyalty},
    )

    costs = (
        DrillingCost_TimeFE,
        DrillingCost_TimeFE_rigrate,
        DrillingCost_TimeFE_costdiffs,
        DrillingCost_TimeFE_rig_costdiffs,
        DrillingCost_DoubleTimeFE,
    )

    scraps = (
        ScrapValue_Zero,
        ScrapValue_Constant,
        ScrapValue_Price,
    )

    exts = (
        ExtensionCost_Constant,
    )

    for R in revs
        for C in costs
            for S in scraps
                for E in exts
                    DR = DrillReward{R,C,E,S}
                    @assert precompile(Tuple{Type{DataDrillPrimitive}, DR, String})
                    @assert precompile(Tuple{Type{DataProduce},        DR, String})
                    @assert precompile(Tuple{Type{DataRoyalty},        DR, String})
                    @assert precompile(Tuple{Type{DynamicDrillModel},  DR, Float64,LeasedProblem,Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},UnitRange{Float64}},SparseArrays.SparseMatrixCSC{Float64,Int64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Bool})
                    @assert precompile(Tuple{Core.kwftype(typeof(SDL.GridTransition)),NamedTuple{(:minp,),Tuple{Float64}},typeof(SDL.GridTransition),DataDrillPrimitive{DR, ExogTimeVars{2,Tuple{Float64,Float64},StepRange{Dates.Date,Dates.Month}},Tuple{Float64,Float64},Int64,LeasedProblem},Float64,Int64})
                end
            end
        end
    end

    @assert precompile(Tuple{typeof(SDL.parse_commandline)})
    @assert precompile(Tuple{typeof(SDL.PsiSpace),Int64})
    @assert precompile(Tuple{typeof(grad_simloglik!),SubArray{Float64,1,SharedArrays.SharedArray{Float64,2},Tuple{Array{Int64,1},Int64},  false}, ObservationGroup{DataProduce{ProductionModel,Float64},Int64},Array{Float64,1},SimulationDraws{Float64,1,SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true}}})
    @assert precompile(Tuple{typeof(simloglik!),     SubArray{Float64,1,SharedArrays.SharedArray{Float64,2},Tuple{UnitRange{Int64},Int64},true},  ObservationGroup{DataRoyalty{RoyaltyModel,Int64,Array{Float64,1},Float64},Int64},SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true},SimulationDraws{Float64,1,SubArray{Float64,1,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true}},Bool})
    @assert precompile(Tuple{typeof(update_reo!),Array{Float64,1}})
    @assert precompile(Tuple{Core.kwftype(typeof(start_up_workers)),NamedTuple{(:nprocs,),Tuple{Int64}},typeof(start_up_workers),Base.EnvDict})

end


text:
using JuMP
using SeaPearl

struct Edge
    vertex1     :: Int
    vertex2     :: Int
end

struct InputData
    edges               :: Array{Edge}
    numberOfEdges       :: Int
    numberOfVertices    :: Int
end

struct OutputData
    numberOfColors      :: Int
    edgeColors          :: Array{Int}
    optimality          :: Bool
end

include("IOmanager.jl")

function outputFromSeaPearl(sol::SeaPearl.Solution; optimality=false)
    numberOfColors = 0
    edgeColors = Int[]

    for i in 1:length(keys(sol))
        color = sol[string(i)]
        if !(color in edgeColors)
            numberOfColors += 1
        end
        push!(edgeColors, color)
    end

    return OutputData(numberOfColors, edgeColors, optimality)
end

struct GraphColoringVariableSelection  <: SeaPearl.AbstractVariableSelection{true}
    sortedPermutation
    degrees
end
function (vs::GraphColoringVariableSelection)(model::SeaPearl.CPModel)
    sortedPermutation, degrees = vs.sortedPermutation, vs.degrees
    maxDegree = 0
    toReturn = nothing
    for i in sortedPermutation
        if !SeaPearl.isbound(model.variables[string(i)])
            if isnothing(toReturn)
                toReturn = model.variables[string(i)]
                maxDegree = degrees[i]
            end
            if degrees[i] < maxDegree
                return toReturn
            end

            if length(model.variables[string(i)].domain) < length(toReturn.domain)
                toReturn = model.variables[string(i)]
            end
        end
    end
    return toReturn
end

function solve_coloring_MOI(input_file; benchmark=false)

    # use input data to fill the model
    input = getInputData(input_file)

    model = SeaPearl.Optimizer()

    for i in 1:input.numberOfVertices
        MOI.add_constrained_variable(model, MOI.Interval(1, input.numberOfVertices))
    end

    degrees = zeros(Int, input.numberOfVertices)

    for e in input.edges
        MOI.add_constraint(model, MOI.VectorOfVariables([MOI.VariableIndex(e.vertex1), MOI.VariableIndex(e.vertex2)]), SeaPearl.VariablesEquality(false))
        degrees[e.vertex1] += 1
        degrees[e.vertex2] += 1
    end

    sortedPermutation = sortperm(degrees; rev=true)

    # define the heuristic used for variable selection
    variableheuristic = GraphColoringVariableSelection(sortedPermutation, degrees)
    MOI.set(model, SeaPearl.MOIVariableSelectionAttribute(), variableheuristic)

    MOI.set(model, SeaPearl.VariableSelection(), variableheuristic)

    solution = MOI.optimize!(model)

    output = outputFromSeaPearl(solution)
    printSolution(output)
end

function solve_coloring_JuMP(input_file; benchmark=false)

    # use input data to fill the model
    input = getInputData(input_file)

    model = Model(SeaPearl.Optimizer)

    @variable(model, 1 <= x[1:input.numberOfVertices] <= input.numberOfVertices)

    degrees = zeros(Int, input.numberOfVertices)
    for e in input.edges
        @constraint(model, [x[e.vertex1], x[e.vertex2]] in SeaPearl.NotEqualSet())
        #update degrees
        degrees[e.vertex1] += 1
        degrees[e.vertex2] += 1
    end

    @variable(model, 1 <= y <= input.numberOfVertices)
    for i in 1:input.numberOfVertices
        @constraint(model, x[i] <= y)
    end
    @objective(model, Min, y)

    sortedPermutation = sortperm(degrees; rev=true)

    # define the heuristic used for variable selection
    variableheuristic = GraphColoringVariableSelection(sortedPermutation, degrees)
    MOI.set(model, SeaPearl.MOIVariableSelectionAttribute(), variableheuristic)

    optimize!(model)
    status = MOI.get(model, MOI.TerminationStatus())

    # output = outputFromSeaPearl(solution)
    # printSolution(output)
    println(model)
    println(status)
    println(has_values(model))
    println(value.(x))
    println(value(y))

end

max_stars_repo_name: corail-research/SeaPearlZoo

max_issues_repo_name: corail-research/SeaPearlZoo

max_forks_repo_name: corail-research/SeaPearlZoo