text:
"Holds the tableau of an Specialised Partitioned Additive Runge-Kutta method for Variational systems."
struct TableauSLRK{DT <: Number} <: AbstractTableau{DT}
    name::Symbol
    o::Int
    s::Int
    r::Int

    q::Tableau{DT}
    p::Tableau{DT}

    q̃::Tableau{DT}
    p̃::Tableau{DT}

    ω::Matrix{DT}
    d::Vector{DT}

    function TableauSLRK(name::Symbol, o::Int, s::Int,
                        q::Tableau{DT}, p::Tableau{DT},
                        q̃::Tableau{DT}, p̃::Tableau{DT},
                        ω::Matrix{DT}, d=DT[]) where {DT}

        @assert s > 0 "Number of stages s must be > 0"

        @assert s==q.s==p.s==q̃.s==p̃.s
        @assert size(ω,1)==s
        @assert size(ω,2)==s+1

        @assert length(d)==0 || length(d)==s

        new{DT}(name, o, s, s, q, p, q̃, p̃, ω, d)
    end
end

"Parameters for right-hand side function of Specialised Partitioned Additive Runge-Kutta methods for Variational systems."
const ParametersSLRK = AbstractParametersSPARK{:slrk}

@doc raw"""
Specialised Partitioned Additive Runge-Kutta integrator for degenerate
variational systems with projection on secondary constraint.

This integrator solves the following system of equations for the internal stages,
```math
\begin{aligned}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a^1_{ij} V_{n,j} + h \sum \limits_{j=1}^{\sigma} a^2_{ij} \tilde{\Lambda}_{n,j} , & i &= 1, ..., s , \\
P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \bar{a}^1_{ij} F^1_{n,j} + h \sum \limits_{j=1}^{s} \bar{a}^2_{ij} F^2_{n,j} + h \sum \limits_{j=1}^{\sigma} \bar{a}^3_{ij} \tilde{F}^3_{n,j} , & i &= 1, ..., s , \\
0 &= \Phi_{n,i} , & i &= 1, ..., s ,
\end{aligned}
```
the projective stages
```math
\begin{aligned}
\tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \alpha^1_{ij} \tilde{V}_{n,j} + h \sum \limits_{j=1}^{\sigma} \alpha^2_{ij} \tilde{\Lambda}_{n,j} , & i &= 1, ..., \sigma , \\
\tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \bar{\alpha}^1_{ij} F^1_{n,j} + h \sum \limits_{j=1}^{s} \bar{\alpha}^2_{ij} \tilde{F}^2_{n,j} + h \sum \limits_{j=1}^{\sigma} \bar{\alpha}^3_{ij} \tilde{F}^3_{n,j} , & i &= 1, ..., \sigma , \\
0 &= \tilde{\Phi}_{n,i} , & i &= 1, ..., \sigma , \\
0 &= \sum \limits_{j=1}^{\sigma} \omega_{ij} \tilde{\Psi}_{n,i} , & i &= 1, ..., \sigma-1 ,
\end{aligned}
```
and update rule
```math
\begin{aligned}
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b^1_{i} V_{n,i} + h \sum \limits_{i=1}^{\sigma} b^2_{i} \tilde{\Lambda}_{n,i} , \\
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b^1_{i} F^1_{n,i} + h \sum \limits_{i=1}^{s} b^2_{i} F^2_{n,i} + h \sum \limits_{i=1}^{\sigma} b^3_{i} \tilde{F}^3_{n,i} , \\
0 &= \phi (q_{n+1}, p_{n+1}) ,
\end{aligned}
```
with definitions
```math
\begin{aligned}
F^1_{n,i} &= \nabla H (Q_{n,i}) , & i &= 1, ..., s , \\
F^2_{n,i} &= \nabla \vartheta (Q_{n,i}) \cdot V_{n,i} , & i &= 1, ..., s , \\
\tilde{F}^3_{n,i} &= \nabla \phi (\tilde{Q}_{n,i}, \tilde{P}_{n,i}) \cdot \tilde{\Lambda}_{n,i} , & i &= 1, ..., \sigma , \\
\Phi_{n,i} &= \phi (Q_{n,i}, P_{n,i}) = P_{n,i} - \vartheta (Q_{n,i}) , & i &= 1, ..., s , \\
\tilde{\Phi}_{n,i} &= \phi (\tilde{Q}_{n,i}, \tilde{P}_{n,i}) = \tilde{P}_{n,i} - \vartheta (\tilde{Q}_{n,i}) , & i &= 1, ..., \sigma , \\
\tilde{\Psi}_{n,i} &= \psi (\tilde{Q}_{n,i}, \tilde{V}_{n,i}, \tilde{P}_{n,i}, \tilde{F}_{n,i}) = \tilde{F}_{n,i} - \tilde{V}_{n,i} \cdot \nabla \vartheta (\tilde{Q}_{n,i}) , & i &= 1, ..., \sigma ,
\end{aligned}
```
so that
```math
\begin{aligned}
F^1_{n,i} + F^2_{n,i} &= \frac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
\phi (Q_{n,i}, P_{n,i}) &= P_{n,i} - \frac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
\psi (\tilde{Q}_{n,i}, \tilde{V}_{n,i}, \tilde{P}_{n,i}, \tilde{F}_{n,i}) &= \tilde{F}_{n,i} - \tilde{V}_{n,i} \cdot \nabla \frac{\partial L}{\partial v} (\tilde{Q}_{n,i}, \tilde{V}_{n,i}) , & i &= 1, ..., \sigma .
\end{aligned}
```
"""
struct IntegratorSLRK{DT, TT, D, S, PT <: ParametersSLRK{DT,TT,D,S,S},
                                    ST <: NonlinearSolver{DT},
                                    IT <: InitialGuessIODE{TT}} <: AbstractIntegratorVSPARK{DT,TT,D,S,S}
    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    function IntegratorSLRK(params::ParametersSLRK{DT,TT,D,S}, solver::ST, iguess::IT, caches) where {DT,TT,D,S,ST,IT}
        new{DT, TT, D, S, typeof(params), ST, IT}(params, solver, iguess, caches)
    end

    function IntegratorSLRK{DT,D}(equations::NamedTuple, tableau::TableauSLRK{TT}, Δt::TT) where {DT,TT,D,ST}
        # get number of stages
        S = tableau.s

        N = 4*D*S

        if isdefined(tableau, :d) && length(tableau.d) > 0
            N += D
        end

        # create params
        params = ParametersSLRK{DT,D}(equations, tableau, Δt)

        # create cache dict
        caches = CacheDict(params)

        # create solver
        solver = create_nonlinear_solver(DT, N, params, caches)

        # create initial guess
        iguess = InitialGuessIODE(get_config(:ig_extrapolation), equations[:v̄], equations[:f̄], Δt)

        # create integrator
        IntegratorSLRK(params, solver, iguess, caches)
    end

    function IntegratorSLRK(equation::LDAE{DT}, tableau::TableauSLRK, Δt; kwargs...) where {DT}
        IntegratorSLRK{DT, ndims(equation)}(get_functions(equation), tableau, Δt; kwargs...)
    end
end

GeometricBase.nconstraints(::IntegratorSLRK{DT,TT,D,S}) where {DT,TT,D,S} = D

function initial_guess!(int::IntegratorSLRK{DT}, sol::AtomicSolutionPDAE{DT},
                        cache::IntegratorCacheSPARK{DT}=int.caches[DT]) where {DT}
    for i in 1:nstages(int)
        evaluate!(int.iguess, sol.q̄, sol.p̄, sol.v̄, sol.f̄,
                              sol.q, sol.p, sol.v, sol.f,
                              cache.q̃, cache.p̃, cache.ṽ, cache.f̃,
                              tableau(int).q̃.c[i], tableau(int).p̃.c[i])

        for k in 1:ndims(int)
            int.solver.x[4*(ndims(int)*(i-1)+k-1)+1] = (cache.q̃[k] - sol.q[k])/timestep(int)
            int.solver.x[4*(ndims(int)*(i-1)+k-1)+2] = (cache.p̃[k] - sol.p[k])/timestep(int)
            int.solver.x[4*(ndims(int)*(i-1)+k-1)+3] = cache.ṽ[k]
            int.solver.x[4*(ndims(int)*(i-1)+k-1)+4] = 0
        end
    end

    if isdefined(tableau(int), :d) && length(tableau(int).d) > 0
        for k in 1:ndims(int)
            int.solver.x[4*ndims(int)*pstages(int)+k] = 0
        end
    end
end

function compute_stages!(x::Vector{ST}, cache::IntegratorCacheSPARK{ST,D,S},
                                        params::ParametersSLRK{DT,TT,D,S}) where {ST,DT,TT,D,S}
    local t::TT

    for i in 1:S
        for k in 1:D
            # copy y to Y, Z and Λ
            cache.Yp[i][k] = x[4*(D*(i-1)+k-1)+1]
            cache.Zp[i][k] = x[4*(D*(i-1)+k-1)+2]
            cache.Vp[i][k] = x[4*(D*(i-1)+k-1)+3]
            cache.Λp[i][k] = x[4*(D*(i-1)+k-1)+4]

            # compute Q and P
            cache.Qp[i][k] = params.q[k] + params.Δt * cache.Yp[i][k]
            cache.Pp[i][k] = params.p[k] + params.Δt * cache.Zp[i][k]
        end

        # compute f(X)
        t = params.t + params.Δt * params.tab.p.c[i]
        params.equs[:f](t, cache.Qp[i], cache.Vp[i], cache.Fp[i])
        params.equs[:g](t, cache.Qp[i], cache.Λp[i], cache.Gp[i])
        params.equs[:ϕ](t, cache.Qp[i], cache.Pp[i], cache.Φp[i])
        params.equs[:ψ](t, cache.Qp[i], cache.Pp[i], cache.Vp[i], cache.Fp[i], cache.Ψp[i])
    end

    if length(params.tab.d) > 0
        for k in 1:D
            cache.μ[k] = x[4*D*S+k]
        end
    end

    # compute q and p
    cache.q̃ .= params.q
    cache.p̃ .= params.p
    for i in 1:S
        cache.q̃ .+= params.Δt .* params.tab.q.b[i] .* cache.Vp[i]
        cache.q̃ .+= params.Δt .* params.tab.q̃.b[i] .* cache.Λp[i]
        cache.p̃ .+= params.Δt .* params.tab.p.b[i] .* cache.Fp[i]
        cache.p̃ .+= params.Δt .* params.tab.p̃.b[i] .* cache.Gp[i]
    end

    # compute ϕ(q,p)
    t = params.t + params.Δt
    params.equs[:ϕ](t, cache.q̃, cache.p̃, cache.ϕ̃)
end

"Compute stages of specialised partitioned additive Runge-Kutta methods for variational systems."
function Integrators.function_stages!(y::Vector{ST}, b::Vector{ST}, params::ParametersSLRK{DT,TT,D,S},
                                      caches::CacheDict) where {ST,DT,TT,D,S}

    # get cache for internal stages
    cache = caches[ST]

    compute_stages!(y, cache, params)

    # compute b = - [(Y-AV-AU), (Z-AF-AG), Φ, ωΨ]
    for i in 1:S
        for k in 1:D
            b[4*(D*(i-1)+k-1)+1] = - cache.Yp[i][k]
            b[4*(D*(i-1)+k-1)+2] = - cache.Zp[i][k]
            b[4*(D*(i-1)+k-1)+3] = - cache.Φp[i][k]
            b[4*(D*(i-1)+k-1)+4] = params.tab.ω[i,S+1] * cache.ϕ̃[k]
            for j in 1:S
                b[4*(D*(i-1)+k-1)+1] += params.tab.q.a[i,j] * cache.Vp[j][k]
                b[4*(D*(i-1)+k-1)+1] += params.tab.q̃.a[i,j] * cache.Λp[j][k]
                b[4*(D*(i-1)+k-1)+2] += params.tab.p.a[i,j] * cache.Fp[j][k]
                b[4*(D*(i-1)+k-1)+2] += params.tab.p̃.a[i,j] * cache.Gp[j][k]
                b[4*(D*(i-1)+k-1)+4] += params.tab.ω[i,j]   * cache.Ψp[j][k]
            end
        end
    end

    if length(params.tab.d) > 0
        for i in 1:S
            for k in 1:D
                b[4*(D*(i-1)+k-1)+2] += cache.μ[k] * params.tab.d[i] / params.tab.p.b[i]
                b[4*(D*(i-1)+k-1)+3] += cache.μ[k] * params.tab.d[i] / params.tab.p.b[i]
            end
        end

        for k in 1:D
            b[4*D*S+k] = 0
            for i in 1:S
                b[4*D*S+k] -= cache.Vp[i][k] * params.tab.d[i]
            end
        end
    end
end

function update_solution!(int::IntegratorSLRK{DT,TT}, sol::AtomicSolutionPDAE{DT,TT},
                          cache::IntegratorCacheSPARK{DT}=int.caches[DT]) where {DT,TT}
    # compute final update
    update_solution!(sol.q, cache.Vp, int.params.tab.q.b, timestep(int))
    update_solution!(sol.q, cache.Λp, int.params.tab.q̃.b, timestep(int))
    update_solution!(sol.p, cache.Fp, int.params.tab.p.b, timestep(int))
    update_solution!(sol.p, cache.Gp, int.params.tab.p̃.b, timestep(int))
end


text:
# Zygote loads ForwardDiff, so this adjoint will autmatically be loaded together
# with `using Zygote`.

# FIXME: This is inefficient, replace with the commented code below once Zygote supports it.
@adjoint function poissonbinomial_pdf_fft(x::AbstractArray{T}) where T<:Real
    fft = poissonbinomial_pdf_fft(x)
    return  fft, Δ -> begin
        ((ForwardDiff.jacobian(x -> poissonbinomial_pdf_fft(x), x)::Matrix{T})' * Δ,)
    end
end

# The code below doesn't work because of bugs in Zygote. The above is inefficient.
#=
ZygoteRules.@adjoint function poissonbinomial_pdf_fft(x::AbstractArray{<:Real})
    return ZygoteRules.pullback(poissonbinomial_pdf_fft_zygote, x)
end
function poissonbinomial_pdf_fft_zygote(p::AbstractArray{T}) where {T <: Real}
    n = length(p)
    ω = 2 * one(T) / (n + 1)
    lmax = ceil(Int, n/2)
    x1 = [one(T)/(n + 1)]
    x_lmaxp1 = map(1:lmax) do l
        logz = zero(T)
        argz = zero(T)
        for j=1:n
            zjl = 1 - p[j] + p[j] * cospi(ω*l) + im * p[j] * sinpi(ω * l)
            logz += log(abs(zjl))
            argz += atan(imag(zjl), real(zjl))
        end
        dl = exp(logz)
        return dl * cos(argz) / (n + 1) + dl * sin(argz) * im / (n + 1)
    end
    x_lmaxp2_end = [conj(x[l + 1]) for l in lmax:-1:1 if n + 1 - l > l]
    x = vcat(x1; x_lmaxp1, x_lmaxp2_end)
    y = [sum(x[j] * cis(-π * float(T)(2 * mod(j * k, n)) / n) for j in 1:n) for k in 1:n]
    return max.(0, real.(y))
end
function poissonbinomial_pdf_fft_zygote2(p::AbstractArray{T}) where {T <: Real}
    n = length(p)
    ω = 2 * one(T) / (n + 1)
    x = Vector{Complex{T}}(undef, n+1)
    lmax = ceil(Int, n/2)
    x[1] = one(T)/(n + 1)
    for l=1:lmax
        logz = zero(T)
        argz = zero(T)
        for j=1:n
            zjl = 1 - p[j] + p[j] * cospi(ω*l) + im * p[j] * sinpi(ω * l)
            logz += log(abs(zjl))
            argz += atan(imag(zjl), real(zjl))
        end
        dl = exp(logz)
        x[l + 1] = dl * cos(argz) / (n + 1) + dl * sin(argz) * im / (n + 1)
        if n + 1 - l > l
            x[n + 1 - l + 1] = conj(x[l + 1])
        end
    end
    max.(0, real.(_dft_zygote(copy(x))))
end
function _dft_zygote(x::Vector{T}) where T
    n = length(x)
    y = Zygote.Buffer(zeros(complex(float(T)), n))
    @inbounds for j = 0:n-1, k = 0:n-1
        y[k+1] += x[j+1] * cis(-π * float(T)(2 * mod(j * k, n)) / n)
    end
    return copy(y)
end
=#

ext: jl


text:
const ln2pi = log(2*pi);

"""
    norm2(x::Array{Float64,1})

Computes ||x||^2 more precisely than norm(x)^2.
"""
function norm2(x::Array{Float64,1})
    return sum(x.^2)
end

"""
    norm2(X::Array{Float64,2})

Computes the square of the Frobenius matrix norm ||X||^2_FRO more precisely than vecnorm(A,2)^2.
"""
function norm2(X::Array{Float64,2})
    return sum(X.^2)
end

"""
   delta(new, old)

Computes delta for convergence control for a given variable.
"""
function delta(new, old)
    return norm(old - new)/norm(old)
end

"""
    scaleY(Y::Array{Float64,2})

Scales down a 2 dimensional array so it has approx. standard normal distribution.
"""
function scaleY(Y::Array{Float64,2})
    L, M = size(Y)
    mu = mean(Y,2);
    sigma = var(Y,2);

    # if there are NaN, then sigma is zero for a given row ->
    # the scaled down row is also zero
    # but we treat this more economically by setting the denominator for a given row to one
    # also, deal with numerical zeroes
    den = sigma
    den[abs(den) .<= 1e-15] = 1.0
    den[den .== 0.0] = 1.0
    den = repmat(sqrt(den),1, M)
    nom = Y - repmat(mu, 1, M)
    nom[abs(nom) .<= 1e-8] = 0.0
    Y = nom./den
    return Y
end

"""
    scaleY(Y::Array{Float32,2})

Scales down a 2 dimensional array so it has approx. standard normal distribution.
"""
function scaleY(Y::Array{Float32,2})
    Y64 = convert(Array{Float64,2}, Y)
    Y64 = scaleY(Y64)
    return Y64
    #return convert(Array{Float32,2}, Y64)
end

"""
    preprocess(Y::Array{Float64,2}, lambda::Float64; verb = false)

Scales Y to be normal, deletes unnecessary rows and then multiplies it by lambda
for better algorithm convergence.
"""
function preprocess(Y::Array{Float64,2}, lambda::Float64; verb = false)
    sY = scaleY(Y)
    rowsums = sum(abs(sY), 2)
    L, M = size(Y)
    used_rows = 1:L
    used_rows = used_rows[rowsums .>= 1e-5]
    # there are either zeroes or the variance is so small that the row carries no information

    if verb
        println("Original problem size: $L rows, $(L-size(used_rows)[1]) rows not relevant and are not used.")
    end

    return lambda*sY[used_rows,:]
end

"""
    preprocess(Y::Array{Float32,2}, lambda::Float64; verb = false)

Scales Y to be normal, deletes unnecessary rows and then multiplies it by lambda
for better algorithm convergence.
"""
function preprocess(Y::Array{Float32,2}, lambda::Float64; verb = false)
    Y64 = convert(Array{Float64,2}, Y)
    return preprocess(Y64, lambda, verb = verb)
end

"""
    traceXTY(X, Y)

Computes tr(X^TY) more effectively than trace(X'*Y).
"""
function traceXTY(X::Array{Float64,2}, Y::Array{Float64,2})
    return(sum(X.*Y))
end

"""
    normalEntropy(Sigma)

Computes the entropy of a multivariate normal distribution with covariance Sigma.
"""
function normalEntropy(Sigma::Array{Float64,2})
    m, n = size(Sigma)
    if n != m
        error("Sigma must be square!")
    end
    d = det(Sigma)
    # numerical issues for too small determinants
    if d < eps(0.0) # smallest nonzero float
        d = eps(0.0)
    end

    return m/2 + m/2*ln2pi + 1/2*log(d)
end

"""
    normalEntropy(diagSigma)

Computes the entropy of a multivariate normal distribution with
diagonal of covariance diagSigma.
"""
function normalEntropy(diagSigma::Array{Float64,1})
    n = size(diagSigma)[1]

    return n/2 + n/2*ln2pi + 1/2*sum(log(diagSigma))
end

"""
    gammaEntropy(a, b)

Computes the entropy of a gamma distribution with shape a and rate b.
"""
function gammaEntropy(a::Float64, b::Float64)
    return a + log(b) + lgamma(a) + (1-a)*digamma(a)
end

"""
    gammaELn(a, b)

Computes E[ln(x)] where x follows gamma distribution with shape a and rate b.
"""
function gammaELn(a::Float64, b::Float64)
    return digamma(a) - log(b)
end


text:
using Plots
using Interpolations

pyplot()

x = [0.0, 0.2, 0.5, 0.6, 0.9, 1.0]
y = [10.0, 0.0, -5.0, 10.0, -8.0, -2.0]
grid = 0.0:0.002:1.0

itypes = [LinearMonotonicInterpolation(),
    FiniteDifferenceMonotonicInterpolation(),
    CardinalMonotonicInterpolation(0.0),
    CardinalMonotonicInterpolation(0.5),
    CardinalMonotonicInterpolation(1.0),
    FritschCarlsonMonotonicInterpolation(),
    FritschButlandMonotonicInterpolation(),
    SteffenMonotonicInterpolation()
    ]

scatter(x, y, label = "data")
for itype in itypes
    itp = interpolate(x, y, itype)
    plot!([x for x in grid], [itp(x) for x in grid], label = string(itype))
end
gui()


text:
#include("../ReadWrite/Header.jl")

"""
**SeisPatch**
*Creation of overlapping 5d patches from a 5d volume*
**IN**
* in: input filename
* out: prefix for output filenames
* style="sxsygxgy"
* min_isx=0
* max_isx=0
* min_isy=0
* max_isy=0
* min_igx=0
* max_igx=0
* min_igy=0
* max_igy=0
* min_imx=0
* max_imx=0
* min_imy=0
* max_imy=0
* min_ihx=0
* max_ihx=0
* min_ihy=0
* max_ihy=0
* min_ih=0
* max_ih=0
* min_iaz=0
* max_iaz=0
* it_WL=9e9
* it_WO=0
* ix1_WL=9e9
* ix1_WO=0
* ix2_WL=9e9
* ix2_WO=0
* ix3_WL=9e9
* ix3_WO=0
* ix4_WL=9e9
* ix4_WO=0
**OUT**
*Credits: A. Stanton, 2015*
"""

function SeisPatch(in::String,out::String;style="sxsygxgy",min_isx=0,max_isx=0,min_isy=0,max_isy=0,min_igx=0,max_igx=0,min_igy=0,max_igy=0,min_imx=0,max_imx=0,min_imy=0,max_imy=0,min_ihx=0,max_ihx=0,min_ihy=0,max_ihy=0,min_ih=0,max_ih=0,min_iaz=0,max_iaz=0,it_WL=9e9,it_WO=0,ix1_WL=9e9,ix1_WO=0,ix2_WL=9e9,ix2_WO=0,ix3_WL=9e9,ix3_WO=0,ix4_WL=9e9,ix4_WO=0)

	if (style == "sxsygxgy")
		key = ["t","isx","isy","igx","igy"]
		min_ix1 = min_isx
		max_ix1 = max_isx
		min_ix2 = min_isy
		max_ix2 = max_isy
		min_ix3 = min_igx
		max_ix3 = max_igx
		min_ix4 = min_igy
		max_ix4 = max_igy
	elseif (style=="mxmyhxhy")
		key = ["t","imx","imy","ihx","ihy"]
		min_ix1 = min_imx
		max_ix1 = max_imx
		min_ix2 = min_imy
		max_ix2 = max_imy
		min_ix3 = min_ihx
		max_ix3 = max_ihx
		min_ix4 = min_ihy
		max_ix4 = max_ihy
	elseif (style=="mxmyhaz")
		key = ["t","imx","imy","ih","iaz"]
		min_ix1 = min_imx
		max_ix1 = max_imx
		min_ix2 = min_imy
		max_ix2 = max_imy
		min_ix3 = min_ih
		max_ix3 = max_ih
		min_ix4 = min_iaz
		max_ix4 = max_iaz
	elseif (style=="sxsyhxhy")
		key = ["t","isx","isy","ihx","ihy"]
		min_ix1 = min_isx
		max_ix1 = max_isx
		min_ix2 = min_isy
		max_ix2 = max_isy
		min_ix3 = min_ihx
		max_ix3 = max_ihx
		min_ix4 = min_ihy
		max_ix4 = max_ihy
	elseif (style=="gxgyhxhy")
		key = ["t","igx","igy","ihx","ihy"]
		min_ix1 = min_igx
		max_ix1 = max_igx
		min_ix2 = min_igy
		max_ix2 = max_igy
		min_ix3 = min_ihx
		max_ix3 = max_ihx
		min_ix4 = min_ihy
		max_ix4 = max_ihy
	elseif (style=="sxsyhaz")
		key = ["t","isx","isy","ih","iaz"]
		min_ix1 = min_isx
		max_ix1 = max_isx
		min_ix2 = min_isy
		max_ix2 = max_isy
		min_ix3 = min_ih
		max_ix3 = max_ih
		min_ix4 = min_iaz
		max_ix4 = max_iaz
	elseif (style=="gxgyhaz")
		key = ["t","igx","igy","ih","iaz"]
		min_ix1 = min_igx
		max_ix1 = max_igx
		min_ix2 = min_igy
		max_ix2 = max_igy
		min_ix3 = min_ih
		max_ix3 = max_ih
		min_ix4 = min_iaz
		max_ix4 = max_iaz
	else
		error("style not defined.")
	end
	nx1 = max_ix1 - min_ix1 + 1
	nx2 = max_ix2 - min_ix2 + 1
	nx3 = max_ix3 - min_ix3 + 1
	nx4 = max_ix4 - min_ix4 + 1

	filename_data = ParseDataName(in)
	filename_headers = ParseHeaderName(in)
	extent = ReadTextHeader(in)

	nt = extent.n1
	dt = extent.d1
	ot = extent.o1

	it_WL  = it_WL  > nt  ? nt  : it_WL
	ix1_WL = ix1_WL > nx1 ? nx1 : ix1_WL
	ix2_WL = ix2_WL > nx2 ? nx2 : ix2_WL
	ix3_WL = ix3_WL > nx3 ? nx3 : ix3_WL
	ix4_WL = ix4_WL > nx4 ? nx4 : ix4_WL

	tmax = ot + dt*nt
	it_NW = Int(floor(nt/(it_WL-it_WO)))
	ix1_NW = Int(floor(nx1/(ix1_WL-ix1_WO)))
	ix2_NW = Int(floor(nx2/(ix2_WL-ix2_WO)))
	ix3_NW = Int(floor(nx3/(ix3_WL-ix3_WO)))
	ix4_NW = Int(floor(nx4/(ix4_WL-ix4_WO)))

	if (ot + dt*(it_NW-1)*(it_WL-it_WO) + dt*it_WL < tmax)
		it_NW += 1
	end
	if (min_ix1 + (ix1_NW-1)*(ix1_WL-ix1_WO) + ix1_WL < max_ix1)
		ix1_NW += 1
	end
	if (min_ix2 + (ix2_NW-1)*(ix2_WL-ix2_WO) + ix2_WL < max_ix2)
		ix2_NW += 1
	end
	if (min_ix3 + (ix3_NW-1)*(ix3_WL-ix3_WO) + ix3_WL < max_ix3)
		ix3_NW += 1
	end
	if (min_ix4 + (ix4_NW-1)*(ix4_WL-ix4_WO) + ix4_WL < max_ix4)
		ix4_NW += 1
	end

	NW=it_NW*ix1_NW*ix2_NW*ix3_NW*ix4_NW
	patch_list = Patch[]
	patch_names = String[]
	for it_W = 1 : it_NW
		mint = ot + dt*(it_W-1)*(it_WL-it_WO)
		maxt = mint + dt*(it_WL - 1)
		if (maxt >= tmax)
			maxt = tmax
		end
		for ix1_W = 1 : ix1_NW
			minx1 = min_ix1 + (ix1_W-1)*(ix1_WL-ix1_WO)
			maxx1 = minx1 + ix1_WL - 1
			if (maxx1 >= max_ix1)
				maxx1 =  max_ix1
			end
			for ix2_W = 1 : ix2_NW
				minx2 = min_ix2 + (ix2_W-1)*(ix2_WL-ix2_WO)
				maxx2 = minx2 + ix2_WL - 1
				if (maxx2 >= max_ix2)
					maxx2 = max_ix2
				end
				for ix3_W = 1 : ix3_NW
					minx3 = min_ix3 + (ix3_W-1)*(ix3_WL-ix3_WO)
					maxx3 = minx3 + ix3_WL - 1
					if (maxx3 >= max_ix3)
						maxx3 = max_ix3
					end
					for ix4_W = 1 : ix4_NW
						minx4 = min_ix4 + (ix4_W-1)*(ix4_WL-ix4_WO)
						maxx4 = minx4 + ix4_WL - 1
						if (maxx4 >= max_ix4)
							maxx4 = max_ix4
						end
						patch_name = join([out "_" it_W "_" ix1_W "_" ix2_W "_" ix3_W "_" ix4_W])
						minval=[mint minx1 minx2 minx3 minx4]
						maxval=[maxt maxx1 maxx2 maxx3 maxx4]
						patch_names = push!(patch_names,patch_name)
						push!(patch_list,Patch(in,patch_name,key,mint,maxt,minx1,maxx1,minx2,maxx2,minx3,maxx3,minx4,maxx4))
					end

				end
			end
		end
	end
	pmap(grab_patch,patch_list)
	return patch_names
end

type Patch
	in
	name
	key
	mint
	maxt
	minx1
	maxx1
	minx2
	maxx2
	minx3
	maxx3
	minx4
	maxx4
end

function grab_patch(patch)

	minval=[patch.mint patch.minx1 patch.minx2 patch.minx3 patch.minx4]
	maxval=[patch.maxt patch.maxx1 patch.maxx2 patch.maxx3 patch.maxx4]
	SeisWindow(patch.in,patch.name,key=patch.key,minval=minval,maxval=maxval)

end


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

using Compat

# const CUDNN_MAJOR = 4
# const CUDNN_MINOR = 0
# const CUDNN_PATCHLEVEL = 4
# const CUDNN_VERSION = CUDNN_MAJOR * 1000 + CUDNN_MINOR * 100 + CUDNN_PATCHLEVEL
const CUDNN_DIM_MAX = 8
const CUDNN_LRN_MIN_N = 1
const CUDNN_LRN_MAX_N = 16
const CUDNN_LRN_MIN_K = 1.0e-5
const CUDNN_LRN_MIN_BETA = 0.01
const CUDNN_BN_MIN_EPSILON = 1.0e-5

type cudnnContext
end

typealias cudnnHandle_t Ptr{cudnnContext}

# begin enum cudnnStatus_t
typealias cudnnStatus_t UInt32
const CUDNN_STATUS_SUCCESS = (UInt32)(0)
const CUDNN_STATUS_NOT_INITIALIZED = (UInt32)(1)
const CUDNN_STATUS_ALLOC_FAILED = (UInt32)(2)
const CUDNN_STATUS_BAD_PARAM = (UInt32)(3)
const CUDNN_STATUS_INTERNAL_ERROR = (UInt32)(4)
const CUDNN_STATUS_INVALID_VALUE = (UInt32)(5)
const CUDNN_STATUS_ARCH_MISMATCH = (UInt32)(6)
const CUDNN_STATUS_MAPPING_ERROR = (UInt32)(7)
const CUDNN_STATUS_EXECUTION_FAILED = (UInt32)(8)
const CUDNN_STATUS_NOT_SUPPORTED = (UInt32)(9)
const CUDNN_STATUS_LICENSE_ERROR = (UInt32)(10)
# end enum cudnnStatus_t

type cudnnTensorStruct
end

typealias cudnnTensorDescriptor_t Ptr{cudnnTensorStruct}

type cudnnConvolutionStruct
end

typealias cudnnConvolutionDescriptor_t Ptr{cudnnConvolutionStruct}

type cudnnPoolingStruct
end

typealias cudnnPoolingDescriptor_t Ptr{cudnnPoolingStruct}

type cudnnFilterStruct
end

typealias cudnnFilterDescriptor_t Ptr{cudnnFilterStruct}

type cudnnLRNStruct
end

typealias cudnnLRNDescriptor_t Ptr{cudnnLRNStruct}

type cudnnActivationStruct
end

typealias cudnnActivationDescriptor_t Ptr{cudnnActivationStruct}

# begin enum cudnnDataType_t
typealias cudnnDataType_t UInt32
const CUDNN_DATA_FLOAT = (UInt32)(0)
const CUDNN_DATA_DOUBLE = (UInt32)(1)
const CUDNN_DATA_HALF = (UInt32)(2)
# end enum cudnnDataType_t

# begin enum cudnnNanPropagation_t
typealias cudnnNanPropagation_t UInt32
const CUDNN_NOT_PROPAGATE_NAN = (UInt32)(0)
const CUDNN_PROPAGATE_NAN = (UInt32)(1)
# end enum cudnnNanPropagation_t

# begin enum cudnnTensorFormat_t
typealias cudnnTensorFormat_t UInt32
const CUDNN_TENSOR_NCHW = (UInt32)(0)
const CUDNN_TENSOR_NHWC = (UInt32)(1)
# end enum cudnnTensorFormat_t

# begin enum cudnnAddMode_t
typealias cudnnAddMode_t UInt32
const CUDNN_ADD_IMAGE = (UInt32)(0)
const CUDNN_ADD_SAME_HW = (UInt32)(0)
const CUDNN_ADD_FEATURE_MAP = (UInt32)(1)
const CUDNN_ADD_SAME_CHW = (UInt32)(1)
const CUDNN_ADD_SAME_C = (UInt32)(2)
const CUDNN_ADD_FULL_TENSOR = (UInt32)(3)
# end enum cudnnAddMode_t

# begin enum cudnnConvolutionMode_t
typealias cudnnConvolutionMode_t UInt32
const CUDNN_CONVOLUTION = (UInt32)(0)
const CUDNN_CROSS_CORRELATION = (UInt32)(1)
# end enum cudnnConvolutionMode_t

# begin enum cudnnConvolutionFwdPreference_t
typealias cudnnConvolutionFwdPreference_t UInt32
const CUDNN_CONVOLUTION_FWD_NO_WORKSPACE = (UInt32)(0)
const CUDNN_CONVOLUTION_FWD_PREFER_FASTEST = (UInt32)(1)
const CUDNN_CONVOLUTION_FWD_SPECIFY_WORKSPACE_LIMIT = (UInt32)(2)
# end enum cudnnConvolutionFwdPreference_t

# begin enum cudnnConvolutionFwdAlgo_t
typealias cudnnConvolutionFwdAlgo_t UInt32
const CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_GEMM = (UInt32)(0)
const CUDNN_CONVOLUTION_FWD_ALGO_IMPLICIT_PRECOMP_GEMM = (UInt32)(1)
const CUDNN_CONVOLUTION_FWD_ALGO_GEMM = (UInt32)(2)
const CUDNN_CONVOLUTION_FWD_ALGO_DIRECT = (UInt32)(3)
const CUDNN_CONVOLUTION_FWD_ALGO_FFT = (UInt32)(4)
const CUDNN_CONVOLUTION_FWD_ALGO_FFT_TILING = (UInt32)(5)
# end enum cudnnConvolutionFwdAlgo_t

type cudnnConvolutionFwdAlgoPerf_t
    algo::cudnnConvolutionFwdAlgo_t
    status::cudnnStatus_t
    time::Cfloat
    memory::Cint
end

# begin enum cudnnConvolutionBwdFilterPreference_t
typealias cudnnConvolutionBwdFilterPreference_t UInt32
const CUDNN_CONVOLUTION_BWD_FILTER_NO_WORKSPACE = (UInt32)(0)
const CUDNN_CONVOLUTION_BWD_FILTER_PREFER_FASTEST = (UInt32)(1)
const CUDNN_CONVOLUTION_BWD_FILTER_SPECIFY_WORKSPACE_LIMIT = (UInt32)(2)
# end enum cudnnConvolutionBwdFilterPreference_t

# begin enum cudnnConvolutionBwdFilterAlgo_t
typealias cudnnConvolutionBwdFilterAlgo_t UInt32
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_0 = (UInt32)(0)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_1 = (UInt32)(1)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_FFT = (UInt32)(2)
const CUDNN_CONVOLUTION_BWD_FILTER_ALGO_3 = (UInt32)(3)
# end enum cudnnConvolutionBwdFilterAlgo_t

type cudnnConvolutionBwdFilterAlgoPerf_t
    algo::cudnnConvolutionBwdFilterAlgo_t
    status::cudnnStatus_t
    time::Cfloat
    memory::Cint
end

# begin enum cudnnConvolutionBwdDataPreference_t
typealias cudnnConvolutionBwdDataPreference_t UInt32
const CUDNN_CONVOLUTION_BWD_DATA_NO_WORKSPACE = (UInt32)(0)
const CUDNN_CONVOLUTION_BWD_DATA_PREFER_FASTEST = (UInt32)(1)
const CUDNN_CONVOLUTION_BWD_DATA_SPECIFY_WORKSPACE_LIMIT = (UInt32)(2)
# end enum cudnnConvolutionBwdDataPreference_t

# begin enum cudnnConvolutionBwdDataAlgo_t
typealias cudnnConvolutionBwdDataAlgo_t UInt32
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_0 = (UInt32)(0)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_1 = (UInt32)(1)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_FFT = (UInt32)(2)
const CUDNN_CONVOLUTION_BWD_DATA_ALGO_FFT_TILING = (UInt32)(3)
# end enum cudnnConvolutionBwdDataAlgo_t

type cudnnConvolutionBwdDataAlgoPerf_t
    algo::cudnnConvolutionBwdDataAlgo_t
    status::cudnnStatus_t
    time::Cfloat
    memory::Cint
end

# begin enum cudnnSoftmaxAlgorithm_t
typealias cudnnSoftmaxAlgorithm_t UInt32
const CUDNN_SOFTMAX_FAST = (UInt32)(0)
const CUDNN_SOFTMAX_ACCURATE = (UInt32)(1)
const CUDNN_SOFTMAX_LOG = (UInt32)(2)
# end enum cudnnSoftmaxAlgorithm_t

# begin enum cudnnSoftmaxMode_t
typealias cudnnSoftmaxMode_t UInt32
const CUDNN_SOFTMAX_MODE_INSTANCE = (UInt32)(0)
const CUDNN_SOFTMAX_MODE_CHANNEL = (UInt32)(1)
# end enum cudnnSoftmaxMode_t

# begin enum cudnnPoolingMode_t
typealias cudnnPoolingMode_t UInt32
const CUDNN_POOLING_MAX = (UInt32)(0)
const CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING = (UInt32)(1)
const CUDNN_POOLING_AVERAGE_COUNT_EXCLUDE_PADDING = (UInt32)(2)
# end enum cudnnPoolingMode_t

# begin enum cudnnActivationMode_t
typealias cudnnActivationMode_t UInt32
const CUDNN_ACTIVATION_SIGMOID = (UInt32)(0)
const CUDNN_ACTIVATION_RELU = (UInt32)(1)
const CUDNN_ACTIVATION_TANH = (UInt32)(2)
const CUDNN_ACTIVATION_CLIPPED_RELU = (UInt32)(3)
# end enum cudnnActivationMode_t

# begin enum cudnnLRNMode_t
typealias cudnnLRNMode_t UInt32
const CUDNN_LRN_CROSS_CHANNEL_DIM1 = (UInt32)(0)
# end enum cudnnLRNMode_t

# begin enum cudnnDivNormMode_t
typealias cudnnDivNormMode_t UInt32
const CUDNN_DIVNORM_PRECOMPUTED_MEANS = (UInt32)(0)
# end enum cudnnDivNormMode_t

# begin enum cudnnBatchNormMode_t
typealias cudnnBatchNormMode_t UInt32
const CUDNN_BATCHNORM_PER_ACTIVATION = (UInt32)(0)
const CUDNN_BATCHNORM_SPATIAL = (UInt32)(1)
# end enum cudnnBatchNormMode_t


text:
__precompile__()

module InteractUIkit

using Reexport
@reexport using InteractBase
using WebIO

import InteractBase:
    toggle,
    getclass

export UIkit

struct UIkit<:InteractBase.WidgetTheme; end

const uikit_icons_min_js = joinpath(@__DIR__, "..", "assets", "uikit-icons.min.js")
const uikit_min_css = joinpath(@__DIR__, "..", "assets", "uikit.min.css")
const uikit_min_js = joinpath(@__DIR__, "..", "assets", "uikit.min.js")
const toggle_switch_css = joinpath(@__DIR__, "..", "assets", "toggle-switch.css")
const font_awesome = isdefined(InteractBase, :font_awesome) ? InteractBase.font_awesome : joinpath(Pkg.dir("InteractBase"), "assets", "all.js")

InteractBase.libraries(::UIkit) = [
    uikit_icons_min_js,
    uikit_min_css,
    uikit_min_js,
    toggle_switch_css,
    font_awesome
]

function __init__()
    InteractBase.settheme!(UIkit())
    nothing
end

include("classes.jl")
include("widgets.jl")

end # module


text:
struct LinearDiffusionObservationModel{S, F} <: ObservationModel{Vector{S}, Vector{S}, ContinuousTime}
    observation_matrix::F
    LinearDiffusionObservationModel(C::AbstractMatrix{S}) where S = new{S, typeof(C)}(C)
end

#####################
### BASIC METHODS ###
#####################

drift(model::LinearDiffusionObservationModel)     = model.observation_matrix
state_dim(model::LinearDiffusionObservationModel) = size(drift(model), 2)
obs_dim(model::LinearDiffusionObservationModel)   = size(drift(model), 1)
noise_dim(model::LinearDiffusionObservationModel) = obs_dim(model)

function observation_function(obs_model::LinearDiffusionObservationModel)
    C = drift(obs_model)
    h(x) = C*x
    return h
end

function Base.show(io::IO, ::MIME"text/plain", model::LinearDiffusionObservationModel)
    print(io, "Linear diffusion process model for the observation
    type of hidden state:                   ", state_dim(model),"-dimensional vector
    type of observation:                    ", obs_dim(model),"-dimensional vector
    number of independent Brownian motions: ", noise_dim(model))
end

function Base.show(io::IO, model::LinearDiffusionObservationModel)
    print(io, obs_dim(model),"-dimensional linear diffusion with ", noise_dim(model), "-dimensional Brownian motion")
end

function Base.show(io::IO, ::MIME"text/plain", model::Type{LinearDiffusionObservationModel{S,T}}) where {S,T}
    print(io, "LinearDiffusionObservationModel{", S, ",", T, "}")
end

function Base.show(io::IO, model::Type{LinearDiffusionObservationModel{S,T}}) where {S,T}
    print(io, "LinearDiffusionObservationModel")
end

######################
### TIME EVOLUTION ###
######################

function (model::LinearDiffusionObservationModel)(x::AbstractVector{T}, dt) where T
    dV = randn(T, noise_dim(model))
    return drift(model) * x * dt + dV * sqrt(dt)
end

function (model::LinearDiffusionObservationModel)(x::AbstractMatrix{T}, dt) where T
    dV = randn(T, noise_dim(model), size(x, 2))
    return drift(model) * x * dt + dV * sqrt(dt)
end

# convert linear model into general model
function DiffusionObservationModel(ob_mod::LinearDiffusionObservationModel)
    h = observation_function(ob_mod)
    return DiffusionObservationModel{eltype(state_type(ob_mod)), eltype(obs_type(ob_mod)), typeof(h)}(state_dim(ob_mod), obs_dim(ob_mod), h)
end

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
# Advent of Code, Day 12
#
# Find the number of paths from start to end in a graph, allowing repeat
# visits for "large caves" (capitalized), but either no (part 1) or at
# most 2 (part 2) visits for "large caves" (lower case).
#
# AK, 12 Dec 2021

part1 = false       # Set to true for Part 2
paths = []          # List of solutions found
f = "input.txt"     # Input file name

# Try all possible trajectories from the current location
function traverse(path::Vector{String}, node::String)

    # Add this node to a copy of the list of visited nodes
    path = copy(path)
    push!(path, node)

    # If this path ends with "end", add it to solutions
    if node == "end"
        push!(paths, join(path, ','))
        return
    end

    # Otherwise, recursively try each potential next node
    for dest in dests[node]
        if legal(path, dest)
            traverse(path, dest)
        end
    end

end

# Determine if it is legal to visit from the current path to the given location
function legal(path::Vector{String}, node::String)

    # For part 1, this is just uppercase, or not already visited
    if part1
        return !(isSmall(node) && node in path)
    end

    # Can't go back to start
    if node == "start"
        return false
    end

    # Can't allow if any duplicate small caves
    if dupes(filter(isSmall, path))
        return !(isSmall(node) && node in path)
    end

    # Otherwise, okay
    return true

end

# Check whether there are any duplicates in a list
function dupes(l::Vector{String})::Bool
    ctr = Dict{String,Int}()
    for s in l
        if haskey(ctr, s)
            return true
        end
        ctr[s] = 1
    end
    return false
end

# Check whether a name is all lower case
function isSmall(s::String)::Bool
    return s == lowercase(s)
end

# Read input file, accumulate dictionary of possible next destination
# for any node, in either direction
dests = Dict{String,Vector{String}}()
lines = filter(s -> length(s) > 0, split(read(f, String), "\n"))
for l in lines
    s, d = split(l, '-')
    if ! haskey(dests, s)
        dests[s] = []
    end
    if ! haskey(dests, d)
        dests[d] = []
    end
    push!(dests[s], d)
    push!(dests[d], s)
end

# Start traversing from the start node
emptyPath = Vector{String}()
traverse(emptyPath, "start")

# Show number of paths found
println("Number of paths: ", length(paths))


text:
using DataDeps

register(DataDep(
  "monks-problems-1",
  "http://archive.ics.uci.edu/ml/datasets/MONK%27s+Problems",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/monks-problems/monks-1.train",
  "d9c328b6f0e095ab10f7bc8468f1d2141f3d08be73e665c3bf0a62994c922606",
  post_fetch_method=(path -> begin
    # First char of each row is whitespace, remove
    write(path, join(strip.(readlines(path)), "\n"))

    UCIData.process_dataset(path,
        target_index=1,
        feature_indices=2:7,
        id_indices=8,
        categoric_indices=2:7,
        delim=' ',
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
using SafeTestsets, Test, Pkg

const GROUP = get(ENV, "GROUP", "All")

function activate_downstream_env()
    Pkg.activate("downstream")
    Pkg.develop(PackageSpec(path=dirname(@__DIR__)))
    Pkg.instantiate()
end

if GROUP == "All" || GROUP == "Core"
    @safetestset "Macro Test" begin include("macro.jl") end
    @safetestset "Differentiation Test" begin include("diff.jl") end
    @safetestset "Is Linear or Affine Test" begin include("islinear_affine.jl") end
    @safetestset "Overloading Test" begin include("overloads.jl") end
    @safetestset "Build Function Test" begin include("build_function.jl") end
    @safetestset "Build Function Array Test" begin include("build_function_arrayofarray.jl") end
    @safetestset "Build Targets Test" begin include("build_targets.jl") end
    @safetestset "Latexify Test" begin include("latexify.jl") end
end

if GROUP == "Downstream"
    activate_downstream_env()
    #@time @safetestset "ParameterizedFunctions MATLABDiffEq Regression Test" begin include("downstream/ParameterizedFunctions_MATLAB.jl") end
end


text:
const MINVAL = 1.0e-15
const PI = 3.141592653589793
const MAXVAL = 1.0e10
const MINMU = 1.0e-5
const MINIMP = 0.0001
const MAXIMP = 0.9999
const MAXCONPAIR = 50
const MAXVFS = 200
const MAXVFSNAME = 100

const NEQDATA = 7
const NDYN = 10
const NGAIN = 10
const NBIAS = 10
const NREF = 2
const NIMP = 5
const NSOLVER = 1000

@enum mjtDisableBit::Cint begin      # disable default feature bitflags
    DSBL_CONSTRAINT   = 1<<0         # entire constraint solver
    DSBL_EQUALITY     = 1<<1         # equality constraints
    DSBL_FRICTIONLOSS = 1<<2         # joint and tendon frictionloss constraints
    DSBL_LIMIT        = 1<<3         # joint and tendon limit constraints
    DSBL_CONTACT      = 1<<4         # contact constraints
    DSBL_PASSIVE      = 1<<5         # passive forces
    DSBL_GRAVITY      = 1<<6         # gravitational forces
    DSBL_CLAMPCTRL    = 1<<7         # clamp control to specified range
    DSBL_WARMSTART    = 1<<8         # warmstart constraint solver
    DSBL_FILTERPARENT = 1<<9         # remove collisions with parent body
    DSBL_ACTUATION    = 1<<10        # apply actuation forces
    DSBL_REFSAFE      = 1<<11        # integrator safety: make ref[0]>=2*timestep

    NDISABLE          = 12           # number of disable flags
end

@enum mjtEnableBit::Cint begin       # enable optional feature bitflags
    ENBL_OVERRIDE     = 1<<0         # override contact parameters
    ENBL_ENERGY       = 1<<1         # energy computation
    ENBL_FWDINV       = 1<<2         # record solver statistics
    ENBL_SENSORNOISE  = 1<<3         # add noise to sensor data

    # NENABLE           = 4          # number of enable flags
end
const NENABLE = 4

@enum mjtJoint::Cint begin           # type of degree of freedom
    JNT_FREE          = 0            # global position and orientation (quat)       (7)
    JNT_BALL                         # orientation (quat) relative to parent        (4)
    JNT_SLIDE                        # sliding distance along body-fixed axis       (1)
    JNT_HINGE                        # rotation angle (rad) around body-fixed axis  (1)
end

@enum mjtGeom::Cint begin            # type of geometric shape
                                     # regular geom types
    GEOM_PLANE        = 0            # plane
    GEOM_HFIELD                      # height field
    GEOM_SPHERE                      # sphere
    GEOM_CAPSULE                     # capsule
    GEOM_ELLIPSOID                   # ellipsoid
    GEOM_CYLINDER                    # cylinder
    GEOM_BOX                         # box
    GEOM_MESH                        # mesh

    NGEOMTYPES                       # number of regular geom types

                                     # rendering-only geom types: not used in mjModel not counted in mjNGEOMTYPES
    GEOM_ARROW        = 100          # arrow
    GEOM_ARROW1                      # arrow without wedges
    GEOM_ARROW2                      # arrow in both directions
    GEOM_LINE                        # line
    GEOM_SKIN                        # skin
    GEOM_LABEL                       # text label

    GEOM_NONE         = 1001         # missing geom type
end

@enum mjtCamLight::Cint begin        # tracking mode for camera and light
    CAMLIGHT_FIXED    = 0            # pos and rot fixed in body
    CAMLIGHT_TRACK                   # pos tracks body, rot fixed in global
    CAMLIGHT_TRACKCOM                # pos tracks subtree com, rot fixed in body
    CAMLIGHT_TARGETBODY              # pos fixed in body, rot tracks target body
    CAMLIGHT_TARGETBODYCOM           # pos fixed in body rot tracks target subtree com
end

@enum mjtTexture::Cint begin         # type of texture
    TEXTURE_2D        = 0            # 2d texture, suitable for planes and hfields
    TEXTURE_CUBE                     # cube texture, suitable for all other geom types
    TEXTURE_SKYBOX                   # cube texture used as skybox
end

@enum mjtIntegrator::Cint begin      # integrator mode
    INT_EULER         = 0            # semi-implicit Euler
    INT_RK4                          # 4th-order Runge Kutta
end

@enum mjtCollision::Cint begin       # collision mode for selecting geom pairs
    COL_ALL           = 0            # test precomputed and dynamic pairs
    COL_PAIR                         # test predefined pairs only
    COL_DYNAMIC                      # test dynamic pairs only
end

@enum mjtCone::Cint begin            # type of friction cone
    CONE_PYRAMIDAL     = 0           # pyramidal
    CONE_ELLIPTIC                    # elliptic
end

@enum mjtJacobian::Cint begin        # type of constraint Jacobian
    JAC_DENSE          = 0           # dense
    JAC_SPARSE                       # sparse
    JAC_AUTO                         # dense if nv<60 sparse otherwise
end

@enum mjtSolver::Cint begin          # constraint solver algorithm
    SOL_PGS            = 0           # PGS    (dual)
    SOL_CG                           # CG     (primal)
    SOL_NEWTON                       # Newton (primal)
end

@enum mjtEq::Cint begin              # type of equality constraint
    EQ_CONNECT        = 0            # connect two bodies at a point (ball joint)
    EQ_WELD                          # fix relative position and orientation of two bodies
    EQ_JOINT                         # couple the values of two scalar joints with cubic
    EQ_TENDON                        # couple the lengths of two tendons with cubic
    EQ_DISTANCE                      # fix the contact distance betweent two geoms
end

@enum mjtWrap::Cint begin            # type of tendon wrap object
    WRAP_NONE         = 0            # null object
    WRAP_JOINT                       # constant moment arm
    WRAP_PULLEY                      # pulley used to split tendon
    WRAP_SITE                        # pass through site
    WRAP_SPHERE                      # wrap around sphere
    WRAP_CYLINDER                    # wrap around (infinite) cylinder
end

@enum mjtTrn::Cint begin             # type of actuator transmission
    TRN_JOINT         = 0            # force on joint
    TRN_JOINTINPARENT                # force on joint, expressed in parent frame
    TRN_SLIDERCRANK                  # force via slider-crank linkage
    TRN_TENDON                       # force on tendon
    TRN_SITE                         # force on site

    TRN_UNDEFINED     = 1000         # undefined transmission type
end

@enum mjtDyn::Cint begin             # type of actuator dynamics
    DYN_NONE          = 0            # no internal dynamics; ctrl specifies force
    DYN_INTEGRATOR                   # integrator: da/dt = u
    DYN_FILTER                       # linear filter: da/dt = (u-a) / tau
    DYN_MUSCLE                       # piece-wise linear filter with two time constants
    DYN_USER                         # user-defined dynamics type
end

@enum mjtGain::Cint begin            # type of actuator gain
    GAIN_FIXED        = 0            # fixed gain
    GAIN_MUSCLE                      # muscle FLV curve computed by mju_muscleGain()
    GAIN_USER                        # user-defined gain type
end

@enum mjtBias::Cint begin            # type of actuator bias
    BIAS_NONE         = 0            # no bias
    BIAS_AFFINE                      # const + kp*length + kv*velocity
    BIAS_MUSCLE                      # muscle passive force computed by mju_muscleBias()
    BIAS_USER                        # user-defined bias type
end

@enum mjtObj::Cint begin             # type of MujoCo object
    OBJ_UNKNOWN       = 0            # unknown object type
    OBJ_BODY                         # body
    OBJ_XBODY                        # body, used to access regular frame instead of i-frame
    OBJ_JOINT                        # joint
    OBJ_DOF                          # dof
    OBJ_GEOM                         # geom
    OBJ_SITE                         # site
    OBJ_CAMERA                       # camera
    OBJ_LIGHT                        # light
    OBJ_MESH                         # mesh
    OBJ_SKIN                         # skin
    OBJ_HFIELD                       # heightfield
    OBJ_TEXTURE                      # texture
    OBJ_MATERIAL                     # material for rendering
    OBJ_PAIR                         # geom pair to include
    OBJ_EXCLUDE                      # body pair to exclude
    OBJ_EQUALITY                     # equality constraint
    OBJ_TENDON                       # tendon
    OBJ_ACTUATOR                     # actuator
    OBJ_SENSOR                       # sensor
    OBJ_NUMERIC                      # numeric
    OBJ_TEXT                         # text
    OBJ_TUPLE                        # tuple
    OBJ_KEY                          # keyframe
end

@enum mjtConstraint::Cint begin      # type of constraint
    CNSTR_EQUALITY    = 0            # equality constraint
    CNSTR_FRICTION_DOF               # dof friction
    CNSTR_FRICTION_TENDON            # tendon friction
    CNSTR_LIMIT_JOINT                # joint limit
    CNSTR_LIMIT_TENDON               # tendon limit
    CNSTR_CONTACT_FRICTIONLESS       # frictionless contact
    CNSTR_CONTACT_PYRAMIDAL          # frictional contact, pyramidal friction cone
    CNSTR_CONTACT_ELLIPTIC           # frictional contact elliptic friction cone
end

@enum mjtConstraintState::Cint begin # constraint state
    CNSTRSTATE_SATISFIED = 0         # constraint satisfied, zero cost (limit, contact)
    CNSTRSTATE_QUADRATIC             # quadratic cost (equality, friction, limit, contact)
    CNSTRSTATE_LINEARNEG             # linear cost, negative side (friction)
    CNSTRSTATE_LINEARPOS             # linear cost, positive side (friction)
    CNSTRSTATE_CONE                  # squared distance to cone cost (elliptic contact)
end

@enum mjtSensor::Cint begin          # type of sensor
                                     # common robotic sensors attached to a site
    SENS_TOUCH        = 0            # scalar contact normal forces summed over sensor zone
    SENS_ACCELEROMETER               # 3D linear acceleration, in local frame
    SENS_VELOCIMETER                 # 3D linear velocity, in local frame
    SENS_GYRO                        # 3D angular velocity, in local frame
    SENS_FORCE                       # 3D force between site's body and its parent body
    SENS_TORQUE                      # 3D torque between site's body and its parent body
    SENS_MAGNETOMETER                # 3D magnetometer
    SENS_RANGEFINDER                 # scalar distance to nearest geom or site along z-axis

                                     # sensors related to scalar joints tendons, actuators
    SENS_JOINTPOS                    # scalar joint position (hinge and slide only)
    SENS_JOINTVEL                    # scalar joint velocity (hinge and slide only)
    SENS_TENDONPOS                   # scalar tendon position
    SENS_TENDONVEL                   # scalar tendon velocity
    SENS_ACTUATORPOS                 # scalar actuator position
    SENS_ACTUATORVEL                 # scalar actuator velocity
    SENS_ACTUATORFRC                 # scalar actuator force

                                     # sensors related to ball joints
    SENS_BALLQUAT                    # 4D ball joint quaterion
    SENS_BALLANGVEL                  # 3D ball joint angular velocity

                                     # joint and tendon limit sensors in constraint space
    SENS_JOINTLIMITPOS               # joint limit distance-margin
    SENS_JOINTLIMITVEL               # joint limit velocity
    SENS_JOINTLIMITFRC               # joint limit force
    SENS_TENDONLIMITPOS              # tendon limit distance-margin
    SENS_TENDONLIMITVEL              # tendon limit velocity
    SENS_TENDONLIMITFRC              # tendon limit force

                                     # sensors attached to an object with spatial frame: (x)body geom, site, camera
    SENS_FRAMEPOS                    # 3D position
    SENS_FRAMEQUAT                   # 4D unit quaternion orientation
    SENS_FRAMEXAXIS                  # 3D unit vector: x-axis of object's frame
    SENS_FRAMEYAXIS                  # 3D unit vector: y-axis of object's frame
    SENS_FRAMEZAXIS                  # 3D unit vector: z-axis of object's frame
    SENS_FRAMELINVEL                 # 3D linear velocity
    SENS_FRAMEANGVEL                 # 3D angular velocity
    SENS_FRAMELINACC                 # 3D linear acceleration
    SENS_FRAMEANGACC                 # 3D angular acceleration

                                     # sensors related to kinematic subtrees; attached to a body (which is the subtree root)
    SENS_SUBTREECOM                  # 3D center of mass of subtree
    SENS_SUBTREELINVEL               # 3D linear velocity of subtree
    SENS_SUBTREEANGMOM               # 3D angular momentum of subtree

                                     # user-defined sensor
    SENS_USER                        # sensor data provided by mjcb_sensor callback
end

@enum mjtStage::Cint begin           # computation stage
    STAGE_NONE        = 0            # no computations
    STAGE_POS                        # position-dependent computations
    STAGE_VEL                        # velocity-dependent computations
    STAGE_ACC                        # acceleration/force-dependent computations
end

@enum mjtDataType::Cint begin        # data type for sensors
    DATATYPE_REAL     = 0            # real values, no constraints
    DATATYPE_POSITIVE                # positive values; 0 or negative: inactive
    DATATYPE_AXIS                    # 3D unit vector
    DATATYPE_QUATERNION              # unit quaternion
end

@enum mjtLRMode::Cint begin          # mode for actuator length range computation
    LRMODE_NONE   = 0                # do not process any actuators
    LRMODE_MUSCLE                    # process muscle actuators
    LRMODE_MUSCLEUSER                # process muscle and user actuators
    LRMODE_ALL                       # process all actuators
end

#------------------------------ mjLROpt ------------------------------------------------

mutable struct mjLROpt # options for mj_setLengthRange()
    # flags
    mode::Cint         # which actuators to process (mjtLRMode)
    useexisting::Cint  # use existing length range if available
    uselimit::Cint     # use joint and tendon limits if available

    # algorithm parameters
    accel::mjtNum      # target acceleration used to compute force
    maxforce::mjtNum   # maximum force 0: no limit
    timeconst::mjtNum  # time constant for velocity reduction min 0.01
    timestep::mjtNum   # simulation timestep 0: use mjOption.timestep
    inttotal::mjtNum   # total simulation time interval
    inteval::mjtNum    # evaluation time interval (at the end)
    tolrange::mjtNum   # convergence tolerance (relative to range)

    mjLROpt() = new()
end

mutable struct mjVFS
    nfile::Cint
    filename::SVector{MAXVFS, SVector{MAXVFSNAME, UInt8}}
    filesize::SVector{MAXVFS, Cint}
    filedata::SVector{MAXVFS, Ptr{Cvoid}}

    mjVFS() = new()
end

struct mjOption
    timestep::mjtNum
    apirate::mjtNum

    impratio::mjtNum
    tolerance::mjtNum
    noslip_tolerance::mjtNum
    mpr_tolerance::mjtNum

    gravity::SVector{3, mjtNum}
    wind::SVector{3, mjtNum}
    magnetic::SVector{3, mjtNum}
    density::mjtNum
    viscosity::mjtNum

    o_margin::mjtNum
    o_solref::SVector{NREF, mjtNum}
    o_solimp::SVector{NIMP, mjtNum}

    integrator::Cint
    collision::Cint
    cone::Cint
    jacobian::Cint
    solver::Cint
    iterations::Cint
    noslip_iterations::Cint
    mpr_iterations::Cint
    disableflags::Cint
    enableflags::Cint

    mjOption() = new()
end

struct _global
    fovy::Cfloat
    ipd::Cfloat
    linewidth::Cfloat
    glow::Cfloat
    offwidth::Cint
    offheight::Cint
end

struct quality
    shadowsize::Cint
    offsamples::Cint
    numslices::Cint
    numstacks::Cint
    numquads::Cint
end

struct headlight
    ambient::SVector{3, Cfloat}
    diffuse::SVector{3, Cfloat}
    specular::SVector{3, Cfloat}
    active::Cint
end

struct map
    stiffness::Cfloat
    stiffnessrot::Cfloat
    force::Cfloat
    torque::Cfloat
    alpha::Cfloat
    fogstart::Cfloat
    fogend::Cfloat
    znear::Cfloat
    zfar::Cfloat
    haze::Cfloat
    shadowclip::Cfloat
    shadowscale::Cfloat
    actuatortendon::Cfloat
end

struct scale
    forcewidth::Cfloat
    contactwidth::Cfloat
    contactheight::Cfloat
    connect::Cfloat
    com::Cfloat
    camera::Cfloat
    light::Cfloat
    selectpoint::Cfloat
    jointlength::Cfloat
    jointwidth::Cfloat
    actuatorlength::Cfloat
    actuatorwidth::Cfloat
    framelength::Cfloat
    framewidth::Cfloat
    constraint::Cfloat
    slidercrank::Cfloat
end

struct rgba
    fog::SVector{4, Cfloat}
    haze::SVector{4, Cfloat}
    force::SVector{4, Cfloat}
    inertia::SVector{4, Cfloat}
    joint::SVector{4, Cfloat}
    actuator::SVector{4, Cfloat}
    actuatornegative::SVector{4, Cfloat}
    actuatorpositive::SVector{4, Cfloat}
    com::SVector{4, Cfloat}
    camera::SVector{4, Cfloat}
    light::SVector{4, Cfloat}
    selectpoint::SVector{4, Cfloat}
    connect::SVector{4, Cfloat}
    contactpoint::SVector{4, Cfloat}
    contactforce::SVector{4, Cfloat}
    contactfriction::SVector{4, Cfloat}
    contacttorque::SVector{4, Cfloat}
    contactgap::SVector{4, Cfloat}
    rangefinder::SVector{4, Cfloat}
    constraint::SVector{4, Cfloat}
    slidercrank::SVector{4, Cfloat}
    crankbroken::SVector{4, Cfloat}
end

struct mjVisual
    _global::_global
    quality::quality
    headlight::headlight
    map::map
    scale::scale
    rgba::rgba

    mjVisual() = new()
end

struct mjStatistic
    meaninertia::mjtNum
    meanmass::mjtNum
    meansize::mjtNum
    extent::mjtNum
    center::SVector{3, mjtNum}
end

mutable struct mjModel
    nq::Cint
    nv::Cint
    nu::Cint
    na::Cint
    nbody::Cint
    njnt::Cint
    ngeom::Cint
    nsite::Cint
    ncam::Cint
    nlight::Cint
    nmesh::Cint
    nmeshvert::Cint
    nmeshtexvert::Cint
    nmeshface::Cint
    nmeshgraph::Cint
    nskin::Cint
    nskinvert::Cint
    nskintexvert::Cint
    nskinface::Cint
    nskinbone::Cint
    nskinbonevert::Cint
    nhfield::Cint
    nhfielddata::Cint
    ntex::Cint
    ntexdata::Cint
    nmat::Cint
    npair::Cint
    nexclude::Cint
    neq::Cint
    ntendon::Cint
    nwrap::Cint
    nsensor::Cint
    nnumeric::Cint
    nnumericdata::Cint
    ntext::Cint
    ntextdata::Cint
    ntuple::Cint
    ntupledata::Cint
    nkey::Cint
    nuser_body::Cint
    nuser_jnt::Cint
    nuser_geom::Cint
    nuser_site::Cint
    nuser_cam::Cint
    nuser_tendon::Cint
    nuser_actuator::Cint
    nuser_sensor::Cint
    nnames::Cint

    nM::Cint
    nemax::Cint
    njmax::Cint
    nconmax::Cint
    nstack::Cint
    nuserdata::Cint
    nmocap::Cint
    nsensordata::Cint

    nbuffer::Cint

    opt::mjOption
    vis::mjVisual
    stat::mjStatistic

    buffer::Ptr{Cvoid}

    qpos0::Ptr{mjtNum}
    qpos_spring::Ptr{mjtNum}

    body_parentid::Ptr{Cint}
    body_rootid::Ptr{Cint}
    body_weldid::Ptr{Cint}
    body_mocapid::Ptr{Cint}
    body_jntnum::Ptr{Cint}
    body_jntadr::Ptr{Cint}
    body_dofnum::Ptr{Cint}
    body_dofadr::Ptr{Cint}
    body_geomnum::Ptr{Cint}
    body_geomadr::Ptr{Cint}
    body_simple::Ptr{mjtByte}
    body_sameframe::Ptr{mjtByte}
    body_pos::Ptr{mjtNum}
    body_quat::Ptr{mjtNum}
    body_ipos::Ptr{mjtNum}
    body_iquat::Ptr{mjtNum}
    body_mass::Ptr{mjtNum}
    body_subtreemass::Ptr{mjtNum}
    body_inertia::Ptr{mjtNum}
    body_invweight0::Ptr{mjtNum}
    body_user::Ptr{mjtNum}

    jnt_type::Ptr{Cint}
    jnt_qposadr::Ptr{Cint}
    jnt_dofadr::Ptr{Cint}
    jnt_bodyid::Ptr{Cint}
    jnt_group::Ptr{Cint}
    jnt_limited::Ptr{mjtByte}
    jnt_solref::Ptr{mjtNum}
    jnt_solimp::Ptr{mjtNum}
    jnt_pos::Ptr{mjtNum}
    jnt_axis::Ptr{mjtNum}
    jnt_stiffness::Ptr{mjtNum}
    jnt_range::Ptr{mjtNum}
    jnt_margin::Ptr{mjtNum}
    jnt_user::Ptr{mjtNum}

    dof_bodyid::Ptr{Cint}
    dof_jntid::Ptr{Cint}
    dof_parentid::Ptr{Cint}
    dof_Madr::Ptr{Cint}
    dof_simplenum::Ptr{Cint}
    dof_solref::Ptr{mjtNum}
    dof_solimp::Ptr{mjtNum}
    dof_frictionloss::Ptr{mjtNum}
    dof_armature::Ptr{mjtNum}
    dof_damping::Ptr{mjtNum}
    dof_invweight0::Ptr{mjtNum}
    dof_M0::Ptr{mjtNum}

    geom_type::Ptr{Cint}
    geom_contype::Ptr{Cint}
    geom_conaffinity::Ptr{Cint}
    geom_condim::Ptr{Cint}
    geom_bodyid::Ptr{Cint}
    geom_dataid::Ptr{Cint}
    geom_matid::Ptr{Cint}
    geom_group::Ptr{Cint}
    geom_priority::Ptr{Cint}
    geom_sameframe::Ptr{mjtByte}
    geom_solmix::Ptr{mjtNum}
    geom_solref::Ptr{mjtNum}
    geom_solimp::Ptr{mjtNum}
    geom_size::Ptr{mjtNum}
    geom_rbound::Ptr{mjtNum}
    geom_pos::Ptr{mjtNum}
    geom_quat::Ptr{mjtNum}
    geom_friction::Ptr{mjtNum}
    geom_margin::Ptr{mjtNum}
    geom_gap::Ptr{mjtNum}
    geom_user::Ptr{mjtNum}
    geom_rgba::Ptr{Cfloat}

    site_type::Ptr{Cint}
    site_bodyid::Ptr{Cint}
    site_matid::Ptr{Cint}
    site_group::Ptr{Cint}
    site_sameframe::Ptr{mjtByte}
    site_size::Ptr{mjtNum}
    site_pos::Ptr{mjtNum}
    site_quat::Ptr{mjtNum}
    site_user::Ptr{mjtNum}
    site_rgba::Ptr{Cfloat}

    cam_mode::Ptr{Cint}
    cam_bodyid::Ptr{Cint}
    cam_targetbodyid::Ptr{Cint}
    cam_pos::Ptr{mjtNum}
    cam_quat::Ptr{mjtNum}
    cam_poscom0::Ptr{mjtNum}
    cam_pos0::Ptr{mjtNum}
    cam_mat0::Ptr{mjtNum}
    cam_fovy::Ptr{mjtNum}
    cam_ipd::Ptr{mjtNum}
    cam_user::Ptr{mjtNum}

    light_mode::Ptr{Cint}
    light_bodyid::Ptr{Cint}
    light_targetbodyid::Ptr{Cint}
    light_directional::Ptr{mjtByte}
    light_castshadow::Ptr{mjtByte}
    light_active::Ptr{mjtByte}
    light_pos::Ptr{mjtNum}
    light_dir::Ptr{mjtNum}
    light_poscom0::Ptr{mjtNum}
    light_pos0::Ptr{mjtNum}
    light_dir0::Ptr{mjtNum}
    light_attenuation::Ptr{Cfloat}
    light_cutoff::Ptr{Cfloat}
    light_exponent::Ptr{Cfloat}
    light_ambient::Ptr{Cfloat}
    light_diffuse::Ptr{Cfloat}
    light_specular::Ptr{Cfloat}

    mesh_vertadr::Ptr{Cint}
    mesh_vertnum::Ptr{Cint}
    mesh_texcoordadr::Ptr{Cint}
    mesh_faceadr::Ptr{Cint}
    mesh_facenum::Ptr{Cint}
    mesh_graphadr::Ptr{Cint}
    mesh_vert::Ptr{Cfloat}
    mesh_normal::Ptr{Cfloat}
    mesh_texcoord::Ptr{Cfloat}
    mesh_face::Ptr{Cint}
    mesh_graph::Ptr{Cint}

    skin_matid::Ptr{Cint}
    skin_rgba::Ptr{Cfloat}
    skin_inflate::Ptr{Cfloat}
    skin_vertadr::Ptr{Cint}
    skin_vertnum::Ptr{Cint}
    skin_texcoordadr::Ptr{Cint}
    skin_faceadr::Ptr{Cint}
    skin_facenum::Ptr{Cint}
    skin_boneadr::Ptr{Cint}
    skin_bonenum::Ptr{Cint}
    skin_vert::Ptr{Cfloat}
    skin_texcoord::Ptr{Cfloat}
    skin_face::Ptr{Cint}
    skin_bonevertadr::Ptr{Cint}
    skin_bonevertnum::Ptr{Cint}
    skin_bonebindpos::Ptr{Cfloat}
    skin_bonebindquat::Ptr{Cfloat}
    skin_bonebodyid::Ptr{Cint}
    skin_bonevertid::Ptr{Cint}
    skin_bonevertweight::Ptr{Cfloat}

    hfield_size::Ptr{mjtNum}
    hfield_nrow::Ptr{Cint}
    hfield_ncol::Ptr{Cint}
    hfield_adr::Ptr{Cint}
    hfield_data::Ptr{Cfloat}

    tex_type::Ptr{Cint}
    tex_height::Ptr{Cint}
    tex_width::Ptr{Cint}
    tex_adr::Ptr{Cint}
    tex_rgb::Ptr{mjtByte}

    mat_texid::Ptr{Cint}
    mat_texuniform::Ptr{mjtByte}
    mat_texrepeat::Ptr{Cfloat}
    mat_emission::Ptr{Cfloat}
    mat_specular::Ptr{Cfloat}
    mat_shininess::Ptr{Cfloat}
    mat_reflectance::Ptr{Cfloat}
    mat_rgba::Ptr{Cfloat}

    pair_dim::Ptr{Cint}
    pair_geom1::Ptr{Cint}
    pair_geom2::Ptr{Cint}
    pair_signature::Ptr{Cint}
    pair_solref::Ptr{mjtNum}
    pair_solimp::Ptr{mjtNum}
    pair_margin::Ptr{mjtNum}
    pair_gap::Ptr{mjtNum}
    pair_friction::Ptr{mjtNum}

    exclude_signature::Ptr{Cint}

    eq_type::Ptr{Cint}
    eq_obj1id::Ptr{Cint}
    eq_obj2id::Ptr{Cint}
    eq_active::Ptr{mjtByte}
    eq_solref::Ptr{mjtNum}
    eq_solimp::Ptr{mjtNum}
    eq_data::Ptr{mjtNum}

    tendon_adr::Ptr{Cint}
    tendon_num::Ptr{Cint}
    tendon_matid::Ptr{Cint}
    tendon_group::Ptr{Cint}
    tendon_limited::Ptr{mjtByte}
    tendon_width::Ptr{mjtNum}
    tendon_solref_lim::Ptr{mjtNum}
    tendon_solimp_lim::Ptr{mjtNum}
    tendon_solref_fri::Ptr{mjtNum}
    tendon_solimp_fri::Ptr{mjtNum}
    tendon_range::Ptr{mjtNum}
    tendon_margin::Ptr{mjtNum}
    tendon_stiffness::Ptr{mjtNum}
    tendon_damping::Ptr{mjtNum}
    tendon_frictionloss::Ptr{mjtNum}
    tendon_lengthspring::Ptr{mjtNum}
    tendon_length0::Ptr{mjtNum}
    tendon_invweight0::Ptr{mjtNum}
    tendon_user::Ptr{mjtNum}
    tendon_rgba::Ptr{Cfloat}

    wrap_type::Ptr{Cint}
    wrap_objid::Ptr{Cint}
    wrap_prm::Ptr{mjtNum}

    actuator_trntype::Ptr{Cint}
    actuator_dyntype::Ptr{Cint}
    actuator_gaintype::Ptr{Cint}
    actuator_biastype::Ptr{Cint}
    actuator_trnid::Ptr{Cint}
    actuator_group::Ptr{Cint}
    actuator_ctrllimited::Ptr{mjtByte}
    actuator_forcelimited::Ptr{mjtByte}
    actuator_dynprm::Ptr{mjtNum}
    actuator_gainprm::Ptr{mjtNum}
    actuator_biasprm::Ptr{mjtNum}
    actuator_ctrlrange::Ptr{mjtNum}
    actuator_forcerange::Ptr{mjtNum}
    actuator_gear::Ptr{mjtNum}
    actuator_cranklength::Ptr{mjtNum}
    actuator_acc0::Ptr{mjtNum}
    actuator_length0::Ptr{mjtNum}
    actuator_lengthrange::Ptr{mjtNum}
    actuator_user::Ptr{mjtNum}

    sensor_type::Ptr{Cint}
    sensor_datatype::Ptr{Cint}
    sensor_needstage::Ptr{Cint}
    sensor_objtype::Ptr{Cint}
    sensor_objid::Ptr{Cint}
    sensor_dim::Ptr{Cint}
    sensor_adr::Ptr{Cint}
    sensor_cutoff::Ptr{mjtNum}
    sensor_noise::Ptr{mjtNum}
    sensor_user::Ptr{mjtNum}

    numeric_adr::Ptr{Cint}
    numeric_size::Ptr{Cint}
    numeric_data::Ptr{mjtNum}

    text_adr::Ptr{Cint}
    text_size::Ptr{Cint}
    text_data::Ptr{UInt8}

    tuple_adr::Ptr{Cint}
    tuple_size::Ptr{Cint}
    tuple_objtype::Ptr{Cint}
    tuple_objid::Ptr{Cint}
    tuple_objprm::Ptr{mjtNum}

    key_time::Ptr{mjtNum}
    key_qpos::Ptr{mjtNum}
    key_qvel::Ptr{mjtNum}
    key_act::Ptr{mjtNum}

    name_bodyadr::Ptr{Cint}
    name_jntadr::Ptr{Cint}
    name_geomadr::Ptr{Cint}
    name_siteadr::Ptr{Cint}
    name_camadr::Ptr{Cint}
    name_lightadr::Ptr{Cint}
    name_meshadr::Ptr{Cint}
    name_skinadr::Ptr{Cint}
    name_hfieldadr::Ptr{Cint}
    name_texadr::Ptr{Cint}
    name_matadr::Ptr{Cint}
    name_pairadr::Ptr{Cint}
    name_excludeadr::Ptr{Cint}
    name_eqadr::Ptr{Cint}
    name_tendonadr::Ptr{Cint}
    name_actuatoradr::Ptr{Cint}
    name_sensoradr::Ptr{Cint}
    name_numericadr::Ptr{Cint}
    name_textadr::Ptr{Cint}
    name_tupleadr::Ptr{Cint}
    name_keyadr::Ptr{Cint}
    names::Ptr{UInt8}
end


text:
"""
    BetaML

The Beta Machine Learning toolkit
https://github.com/sylvaticus/BetaML.jl

For documentation, please look at the individual modules:

- `?BetaML.Perceptron`: Linear and kernel classifiers
- `?BetaML.Nn`: Neural networks
- `?BetaML.Clusters`: Clustering algorithms and collaborative filering using clusters
- `?BetaML.Utils`: Various utility functions (scale, one-hot, distances, kernels,..)

While the code is organised in different sub-modules, all objects are re-exported at the BetaML root level, hence the functionality of this package can be accessed by simply `using BetaML` and then employing the required function directly.

"""
module BetaML

import MLJModelInterface
const MMI = MLJModelInterface

using ForceImport, Reexport

include("Api.jl")        # Shared names across modules
include("Utils/Utils.jl")      # Utility function
include("Nn/Nn.jl")         # Neural Networks
include("Perceptron/Perceptron.jl") # Perceptron-like algorithms
include("Trees/Trees.jl")      # Decision Trees and ensembles (Random Forests)
include("Clustering/Clustering.jl") # Clustering algorithms

# "Merging" of the modules...
@force    using .Api
@reexport using .Api
@force    using .Utils
@reexport using .Utils
@force    using .Nn
@reexport using .Nn
@force    using .Perceptron
@reexport using .Perceptron
@force    using .Trees
@reexport using .Trees
@force    using .Clustering
@reexport using .Clustering

# ------------------------------------------------------------------------------
#MLJ interface...
const MLJ_PERCEPTRON_MODELS = (PerceptronClassifier, KernelPerceptronClassifier, PegasosClassifier)
const MLJ_TREES_MODELS      = (DecisionTreeClassifier, DecisionTreeRegressor, RandomForestClassifier, RandomForestRegressor)
const MLJ_CLUSTERING_MODELS = (KMeans, KMedoids, GMMClusterer, MissingImputator)
const MLJ_INTERFACED_MODELS = (MLJ_PERCEPTRON_MODELS..., MLJ_TREES_MODELS..., MLJ_CLUSTERING_MODELS...)

function __init__()
    MMI.metadata_pkg.(MLJ_INTERFACED_MODELS,
        name       = "BetaML",
        uuid       = "024491cd-cc6b-443e-8034-08ea7eb7db2b", # see your Project.toml
        url        = "https://github.com/sylvaticus/BetaML.jl",  # URL to your package repo
        julia      = true,          # is it written entirely in Julia?
        license    = "MIT",       # your package license
        is_wrapper = false,    # does it wrap around some other package?
    )
end

end # module


text:
"""
    μϵF_induced_by_ϵjs(pts, ϵF, ϵjs)

## Returns
The measures ``\\{\\mu(b_f)\\}`` of the bins ``b_f \\in B_{\\epsilon_{F}}``
(``B_{\\epsilon_{F}}`` is a partition given by the binning scheme `ϵF`)
induced by the measures of the bins ``b_j \\in B_{\\epsilon_{j}}``
(``B_{\\epsilon_{j}}`` is a partition given by the binning scheme `ϵⱼ`)
for several different ``B_{\\epsilon_{j}}`` constructed from the
binning schemes `ϵⱼs`

This is essentially averaging the invariant measures obtained at different
resolutions, expressing the average at the final resolution given by
the binning scheme `ϵF`.

## Explanation

Consider a cloud of `points` and two separate binning schemes, `ϵF` and `ϵⱼ`.
The binning schemes dictate how to create a rectangular partition covering the
points.

Create the partition ``B_{\\epsilon_{F}}`` from `ϵF` and a separate partition
``B_{\\epsilon_{j}}`` from `ϵⱼ`. Then, compute an invariant measure over each
partition separately.

The boxes in ``B_{\\epsilon_{F}}`` are have sizes different from
the boxes in ``B_{\\epsilon_{j}}``. Thus, the invariant measure associated
with each of the partitions convey something about the dynamics at different
scales. We may want to know how the dynamics exposed by
``B_{\\epsilon_{j}}`` is expressed in terms of another coarse-graining
``B_{\\epsilon_{F}}``.

This function computes the measure of the bins ``b_f \\in B_{\\epsilon_{F}}``
induced by the bins ``b_j \\in B_{\\epsilon_{j}}``. That is, for each ``b_f``,
find the set of bins
``\\{b_k \\in B_{\\epsilon_{j}} : b_k \\cap b_f \\neq \\emptyset \\}``.

Extending this to multiple partitions (many different resolutions),
we can decide on a final rectangular grid (resolution) and compute the measure
over that grid as an average over multiple separate partitions (finer, coarser, or both,
depending on the specific application).

## Arguments
- **`points`**: The points representing observations from the system from
    which to compute an invariant measure.
- **`ϵF`**: Edge lengths along each coordinate axis in the partition ``B_{\\epsilon_{F}}``.
- **`ϵⱼ`**: Edge lengths along each coordinate axis in the partition ``B_{\\epsilon_{j}}``.
"""
function μϵF_induced_by_ϵjs(points, ϵF, ϵjs)
    n_partitions = length(ϵjs)

    μ_across_binsizes = Vector{Vector{Float64}}(undef, n_partitions)

    for i=1:n_partitions
        μ_across_binsizes[i] = μϵF_induced_by_ϵj(points, ϵF, ϵjs[i])
    end

    # Gather in a matrix and transpose. Now the i-th column
    # and a-th row give the measure indiced for the i-th
    # state of the target partition induced by the a-th
    # bin size ϵₐ.
    μ = transpose(hcat(μ_across_binsizes...,))
end

export μϵF_induced_by_ϵjs

"""
    average_measure(points, ϵF, ϵjs)

Express the invariant measure of a rectangular box covering of `points`, at
the resolution given by the binning scheme `ϵⱼ` as the invariant measure at
a separate resolution given by the binning scheme `ϵF`.

See the documentation for [`μϵF_induced_by_ϵjs`](@ref) for more details.

## Arguments
- **`points`**: The points representing observations from the system from
    which to compute an invariant measure.
- **`ϵF`**: Edge lengths along each coordinate axis in the partition ``B_{\\epsilon_{F}}``.
- **`ϵjs`**: Edge lengths along each coordinate axis in the (multiple) partitions ``B_{\\epsilon_{j}}``.
"""
function average_measure(points, ϵF, ϵjs)
    n_partitions = length(ϵjs)

    measures = μϵF_induced_by_ϵjs(points, ϵF, ϵjs)

    # Compute the average measure for the final partition
    n_states_μF = size(measures, 2)
    avg_measure = zeros(Float64, n_states_μF)
    for i in 1:n_states_μF
        avg_measure[i] = (1/n_partitions) * sum(measures[:, i])
    end

    return avg_measure
end

export average_measure


text:
using Distributions
using ZeroInflation
using Test
import Turing: BernoulliLogit
import Turing: BinomialLogit

@testset "ZeroInflation.jl" begin
    @testset "constructing ZeroInflated" begin
        d = zeroinflated(Normal(), 0.25)
        @test d isa ZeroInflation.ZeroInflated
        @test d.mix isa Bernoulli
        @test d.mix.p == 0.25
        @test d.original isa Normal

        d = zeroinflated(Normal(), BernoulliLogit(1.0))
        @test d.mix isa BinomialLogit
        @test d.mix.logitp == 1.0
    end

    @testset "accessing ZeroInflated" begin
        d = zeroinflated(Poisson(2.0), 0.5)
        @test params(d) == (2.0, 0.5)
        @test minimum(d) == 0
        @test maximum(d) == Inf

        d = zeroinflated(Normal(), 0.1)
        @test params(d) == (0.0, 1.0, 0.1)
        @test minimum(d) == -Inf
        @test maximum(d) == Inf

        d = zeroinflated(Uniform(-0.5, 0.5), 0.1)
        @test params(d) == (-0.5, 0.5, 0.1)
        @test minimum(d) == -0.5
        @test maximum(d) == 0.5

        d = zeroinflated(Uniform(-2.0, -1.0), 0.1)
        @test minimum(d) == -2.0
        @test maximum(d) == 0.0
    end

    @testset "rand" begin
        @test rand(zeroinflated(Normal(), 1.0)) == 0.0
        @test rand(zeroinflated(Normal(), 0.0)) != 0.0

        @test rand(zeroinflated(Normal(), 0.0)) isa AbstractFloat
        @test rand(zeroinflated(Poisson(1.0), 0.0)) isa Integer
    end

    @testset "logpdf" begin
        d = zeroinflated(Poisson(2.0), 0.2)
        samp = rand(d, 1000)
        pointwise_loglik = logpdf.(d, samp)
        @test length(pointwise_loglik) == length(samp)
        @test isapprox(sum(pointwise_loglik), logpdf(d, samp))
    end
end


text:
include("dfs.jl")
include("ilds.jl")
include("rbs.jl")

"""
    initroot!(toCall::Stack{Function}, ::F, model::CPModel, variableHeuristic::AbstractVariableSelection, valueSelection::ValueSelection) where F <: SearchStrategy
Initialisation function that fill the toCall Stack according to a certain strategy.
    """
function initroot!(toCall::Stack{Function}, ::F, model::CPModel, variableHeuristic::AbstractVariableSelection, valueSelection::ValueSelection) where F <: SearchStrategy
    throw(ErrorException("Search Strategy $(F) (initroot! function  ) not implemented."))
end

"""
search!(model::CPModel, strategy::S, variableHeuristic::AbstractVariableSelection, valueSelection::ValueSelection=BasicHeuristic(); out_solver::Bool=false) where S <: SearchStrategy
Perform a search following a specific strategy in the `model` using `variableHeuristic` to choose which domain will be changed
at each branching and using `valueSelection` to choose how the branching will be done.

"""
function search!(model::CPModel, strategy::S, variableHeuristic::AbstractVariableSelection, valueSelection::ValueSelection=BasicHeuristic(); out_solver::Bool=false) where S <: SearchStrategy
    tic()
    # create env and get first observation
    valueSelection(InitializingPhase, model)

    toCall = Stack{Function}()
    # Starting at the root node with an empty stack
    currentStatus = initroot!(toCall, strategy,  model, variableHeuristic, valueSelection)

    while !isempty(toCall)
        # Stop right away if reached a limit
        if currentStatus == :NodeLimitStop || currentStatus == :SolutionLimitStop || currentStatus == :TimeLimitStop || (out_solver && (currentStatus in [:Infeasible, :FoundSolution]))
            break
        end

        if currentStatus != :SavingState
            # set reward and metrics
            valueSelection(StepPhase, model, currentStatus)
        end

        currentProcedure = pop!(toCall)
        currentStatus = currentProcedure(model)
    end
    # set final reward and last observation
    model.statistics.numberOfSolutions=sum(map(x->!isnothing(x),model.statistics.solutions))
    valueSelection(EndingPhase, model, currentStatus)

    if currentStatus == :NodeLimitStop || currentStatus == :SolutionLimitStop || currentStatus == :TimeLimitStop || (out_solver & (currentStatus in [:Infeasible, :FoundSolution]))
        return currentStatus
    end

    if isa(strategy, DFSearch) && !all(map(x->isnothing(x),model.statistics.solutions)) == 1       #Only the DFS search can give the optimality certificate
        return :Optimal
    elseif !all(map(x->isnothing(x),model.statistics.solutions)) == 1
        return :NonOptimal
    end

    return :Infeasible
end


text:
module ImPlot

using CImGui
import CImGui: ImVec2, ImVec4

for i in instances(CImGui.ImGuiCond_)
    @eval import CImGui.LibCImGui: $(Symbol(i))
end

include("libcimplot.jl")

using .LibCImPlot
export IMPLOT_AUTO, IMPLOT_AUTO_COL
export ImPlotPoint, ImPlotRange, ImPlotLimits

const ImPlotData = Union{Float32, Float64, Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64}
const IMPLOT_ENUMS = [ImPlotFlags_, ImPlotAxisFlags_, ImPlotCol_, ImPlotStyleVar_, ImPlotMarker_,
          ImPlotColormap_, ImPlotLocation_, ImPlotOrientation_, ImPlotYAxis_]

# Export plot flags
for i in IMPLOT_ENUMS
    for j in instances(i)
        @eval export $(Symbol(j))
    end
end

function BeginPlot(title_id::String, x_label = C_NULL, y_label = C_NULL, size::ImVec2 = ImVec2(-1,0);
                    flags = ImPlotFlags_None,
                    x_flags = ImPlotAxisFlags_None,
                    y_flags = ImPlotAxisFlags_None,
                    y2_flags = ImPlotAxisFlags_None,
                    y3_flags = ImPlotAxisFlags_None)::Bool

    LibCImPlot.BeginPlot(title_id, x_label, y_label, size,
                         flags, x_flags, y_flags, y2_flags, y3_flags)
end

include("constructors.jl")
include("lines.jl")
include("stairs.jl")
include("shaded.jl")
include("scatter.jl")
include("heatmap.jl")
include("digital.jl")
include("barchart.jl")
include("piechart.jl")
include("errorbars.jl")
include("stems.jl")
include("other.jl")
include("util.jl")
include("color.jl")
include("styling.jl")

end # module


text:
vis = Visualizer()
open(vis)

# Create Roadway
roadway_opts = MergingRoadwayOptions(lane_width=0.26)
roadway = build_roadway(roadway_opts)

# Create players
T = Float64
p = 2
model = DoubleIntegratorGame(p=p)

players = Vector{Player{T}}(undef, p)
player_opts = HighwayPlayerOptions{model.ni[1],model.mi[1],T}(r_cost=0.22)

players[1] = Player(model, roadway.lane[3],
    Q=Diagonal(SVector{model.ni[1],T}(0.0, 0.4, 0.4, 0.05)),
    R=Diagonal(SVector{model.mi[1],T}(0.1, 0.1)),
    x0=VehicleState(0.15, -0.1, 0.0, 0.7),
    xf=VehicleState(1.9, -0.1, 0.0, 0.7),
    μ=200.0,
    r_col=0.00,
    r_cost=0.16,
    opts=player_opts)
players[2] = Player(model, roadway.lane[6],
    Q=Diagonal(SVector{model.ni[2],T}(0.0, 0.4, 0.4, 0.05)),
    R=Diagonal(SVector{model.mi[2],T}(0.1, 0.1)),
    x0=roadway.lane[5].centerline(0.15, 0.7),
    xf=roadway.lane[5].centerline(1.9, 0.7),
    μ=200.0,
    r_col=0.00,
    r_cost=0.16,
    opts=player_opts)

# Create Scenario
sce = Scenario(model, roadway, players)

# Create AutoProblem
N = 50
dt = 0.05
solver_opts = Options()
solver_opts.inner_print = true
solver_opts.outer_print = true
solver_opts.reg_0 = 1e-5
prob = GameProblem(N, dt, sce, solver_opts)

init_traj!(prob.pdtraj, f=rand, amplitude=1e-8)
solver_opts.shift = 0
@time newton_solve!(prob)

plot!(model, prob.pdtraj.pr)
plot!(prob.stats)

# Visualize trajectory
set_scenario!(vis, sce)
set_env!(vis, VehicleState(2.0, 0.0, 0.0, 0.0))
set_camera_birdseye!(vis, height=7.3)
build_waypoint!(vis, sce.player, N)

@time set_state!(vis, model, sce, Algames.state(prob.pdtraj.pr[10]))
@time set_traj!(vis, model, sce, prob.pdtraj.pr)
@time set_waypoint_traj!(vis, model, sce, prob.pdtraj.pr)

function vizz(prob)
    init_traj = Vector{PrimalDualTraj}()
    sol_traj = Vector{PrimalDualTraj}()

    prob.opts.shift = 2^10
    prob.opts.f_init = rand
    prob.opts.amplitude_init = 1e-1
    solver_opts.reg_0 = 1e-4
    for k = 1:10
        prob.opts.seed = k
        prob.opts.ϵ_opt = 1e-6
        prob.opts.ϵ_dyn = 1e-6
        prob.opts.ϵ_con = 1e-6
        prob.opts.ϵ_sta = 1e-6

        prob.opts.outer_iter = 0
        @time newton_solve!(prob)
        init = deepcopy(prob.pdtraj)
        #
        prob.opts.outer_iter = 20
        @time newton_solve!(prob)
        sol = deepcopy(prob.pdtraj)

        opt = mean(abs.(prob.stats.opt_vio[end].vio))
        sta = mean(abs.(prob.stats.sta_vio[end].vio))
        dyn = mean(abs.(prob.stats.dyn_vio[end].vio))
        con = mean(abs.(prob.stats.con_vio[end].vio))
        opts = prob.opts
        if opt < opts.ϵ_opt && sta < opts.ϵ_sta && dyn < opts.ϵ_dyn && con < opts.ϵ_con
            # plot!(prob.stats)
            # plot!(prob.model, init.pr)
            plot!(prob.model, sol.pr)
            @time set_waypoint_traj!(vis, model, sce, init.pr)
            # sleep(1.0)
            # plot!(prob.model, prob.pdtraj.pr)
            @time set_traj!(vis, model, sce, prob.pdtraj.pr)
            build_waypoint!(vis, sce.player, N, key=k)
            @time set_waypoint_traj!(vis, model, sce, sol.pr, key=k)
            sleep(4.0)
            # push!(init_traj, init)
            push!(sol_traj, sol)
        end
    end
    return init_traj, sol_traj
end

# init_traj, sol_traj = vizz(prob)
# init_traj
# sol_traj

init_traj!(prob.pdtraj; x0=prob.x0, f=prob.opts.f_init, amplitude=prob.opts.amplitude_init, s=prob.opts.shift)
Algames.rollout!(Algames.RK3, prob.model, prob.pdtraj.pr)
set_waypoint_traj!(vis, model, sce, prob.pdtraj.pr)
plot!(model, prob.pdtraj.pr)

function figgg(prob::Algames.GameProblem)
    prob.opts.shift = 2^10
    prob.opts.α_dual = 0.0
    prob.opts.αx_dual = zeros(p)
    prob.opts.f_init = x -> (rand(x) .- 0.5)
    prob.opts.amplitude_init = 1.0

    prob.opts.ϵ_opt = 1e-6
    prob.opts.ϵ_dyn = 1e-6
    prob.opts.ϵ_con = 1e-6
    prob.opts.ϵ_sta = 1e-6

    init = []
    sol = []
    plt = plot()
    plt2 = plot()
    for k = 1:10
        prob.opts.seed = k
        prob.opts.outer_iter = 0
        @time newton_solve!(prob)
        init_traj = deepcopy(prob.pdtraj)

        prob.opts.outer_iter = 20
        @time newton_solve!(prob)
        sol_traj = deepcopy(prob.pdtraj)
        plot!(prob.stats)

        opt = mean(abs.(prob.stats.opt_vio[end].vio))
        sta = mean(abs.(prob.stats.sta_vio[end].vio))
        dyn = mean(abs.(prob.stats.dyn_vio[end].vio))
        con = mean(abs.(prob.stats.con_vio[end].vio))
        opts = prob.opts
        if opt < opts.ϵ_opt && sta < opts.ϵ_sta && dyn < opts.ϵ_dyn && con < opts.ϵ_con
            plot!(model, init_traj.pr, plt=plt)
            plot!(model, sol_traj.pr, plt=plt2)
            push!(sol, sol_traj)
            push!(init, init_traj)
        end
    end
    display(plt)
    return nothing
end

figgg(prob)

plot!


text:
# All of this functionality is utilised only in the AD tests. Can be safely ignored if
# you're concerned with understanding how TemporalGPs works.

using ChainRulesCore: backing

# Functionality to test my testing functionality.
are_harmonised(a::Any, b::AbstractZero) = true
are_harmonised(a::AbstractZero, b::Any) = true
are_harmonised(a::AbstractZero, b::AbstractZero) = true

are_harmonised(a::Number, b::Number) = true

function are_harmonised(a::AbstractArray, b::AbstractArray)
    return all(ab -> are_harmonised(ab...), zip(a, b))
end

are_harmonised(a::Tuple, b::Tuple) = all(ab -> are_harmonised(ab...), zip(a, b))

function are_harmonised(a::Composite{<:Any, <:Tuple}, b::Composite{<:Any, <:Tuple})
    return all(ab -> are_harmonised(ab...), zip(a, b))
end

function are_harmonised(
    a::Composite{<:Any, <:NamedTuple},
    b::Composite{<:Any, <:NamedTuple},
)
    return all(
        name -> are_harmonised(getproperty(a, name), getproperty(b, name)),
        union(fieldnames(typeof(a)), fieldnames(typeof(b))),
    )
end

# Functionality to make it possible to compare different kinds of differentials. It's not
# entirely clear how much sense this makes mathematically, but it seems to work in a
# practical sense at the minute.
harmonise(a::Any, b::AbstractZero) = (a, b)
harmonise(a::AbstractZero, b::Any) = (a, b)
harmonise(a::AbstractZero, b::AbstractZero) = (a, b)

# Resolve ambiguity.
harmonise(a::AbstractZero, b::Composite{<:Any, <:NamedTuple}) = (a, b)

harmonise(a::Number, b::Number) = (a, b)

function harmonise(a::Tuple, b::Tuple)
    vals = map(harmonise, a, b)
    return first.(vals), last.(vals)
end
function harmonise(a::AbstractArray, b::AbstractArray)
    vals = map(harmonise, a, b)
    return first.(vals), last.(vals)
end

function harmonise(a::Composite{<:Any, <:Tuple}, b::Composite{<:Any, <:Tuple})
    vals = map(harmonise, backing(a), backing(b))
    return (Composite{Any}(first.(vals)...), Composite{Any}(last.(vals)...))
end

harmonise(a::Composite{<:Any, <:Tuple}, b::Tuple) = harmonise(a, Composite{Any}(b...))

harmonise(a::Tuple, b::Composite{<:Any, <:Tuple}) = harmonise(Composite{Any}(a...), b)

function harmonise(
    a::Composite{<:Any, <:NamedTuple{names}},
    b::Composite{<:Any, <:NamedTuple{names}},
) where {names}
    vals = map(harmonise, values(backing(a)), values(backing(b)))
    a_harmonised = Composite{Any}(; NamedTuple{names}(first.(vals))...)
    b_harmonised = Composite{Any}(; NamedTuple{names}(last.(vals))...)
    return (a_harmonised, b_harmonised)
end

function harmonise(a::Composite{<:Any, <:NamedTuple}, b::Composite{<:Any, <:NamedTuple})

    # Compute names missing / present in each data structure.
    a_names = propertynames(backing(a))
    b_names = propertynames(backing(b))
    mutual_names = intersect(a_names, b_names)
    all_names = (union(a_names, b_names)..., )
    a_missing_names = setdiff(all_names, a_names)
    b_missing_names = setdiff(all_names, b_names)

    # Construct `Composite`s with the same names.
    a_vals = map(name -> name ∈ a_names ? getproperty(a, name) : Zero(), all_names)
    b_vals = map(name -> name ∈ b_names ? getproperty(b, name) : Zero(), all_names)
    a_unioned_names = Composite{Any}(; NamedTuple{all_names}(a_vals)...)
    b_unioned_names = Composite{Any}(; NamedTuple{all_names}(b_vals)...)

    # Harmonise those composites.
    return harmonise(a_unioned_names, b_unioned_names)
end

function harmonise(a::Composite{<:Any, <:NamedTuple}, b)
    b_names = fieldnames(typeof(b))
    vals = map(name -> getfield(b, name), b_names)
    return harmonise(
        a, Composite{Any}(; NamedTuple{b_names}(vals)...),
    )
end

harmonise(a::Composite{<:Any, <:NamedTuple}, b::AbstractZero) = (a, b)

harmonise(a, b::Composite{<:Any, <:NamedTuple}) = reverse(harmonise(b, a))

# Special-cased handling for `Adjoint`s. Due to our usual AD setup, a differential for an
# Adjoint can be represented either by a matrix or a `Composite`. Both ought to `to_vec` to
# the same thing though, so this should be fine for now, if a little unsatisfactory.
function harmonise(a::Adjoint, b::Composite{<:Adjoint, <:NamedTuple})
    return Composite{Any}(parent=parent(a)), b
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/MortarContact2D.jl/blob/master/LICENSE

"""
    2D mortar contact mechanics for JuliaFEM.

# Problem types

- `Mortar2D` to couple non-conforming meshes (mesh tying).
- `MortarContact2D` to add contact constraints.

"""
module MortarContact2D

using FEMBase, SparseArrays, LinearAlgebra, Statistics

include("mortar2d.jl")
include("contact2d.jl")
export Mortar2D, Contact2D

end


text:
using JuMP, MAT

function load_sedumi_file(filename)
   	file = matopen(filename)
   	if exists(file, "A")
      		A = read(file, "A");
   	elseif exists(file, "At")
      		A = read(file, "At")';
   	else
      		@assert false "The MATLAB file must contain either A or At as a variable"
   	end
   	b = Vector(read(file, "b")[:]);
   	c = Vector(read(file, "c")[:]);
   	K = read(file, "K")
   	close(file)

	# Cast Int entries
   	for key in ["f", "l"]
      		if haskey(K, key)
         			K[key] = Int(K[key])
      		else
         			K[key] = 0
      		end
   	end
	# Cast Vector{Int} entries
   	for key in ["q", "r", "s"]
      		if haskey(K, key)
         			@assert isa(K[key], AbstractArray) || isa(K[key], Number)

         			if isa(K[key], AbstractArray)
            				K[key] = Vector{Int}(K[key][:])
         			else # isa(K[key], Number)
            				K[key] = [Int(K[key])]
         			end
         			K[key] = K[key][K[key] .> 0]
      		else
         			K[key] = Int[]
      		end
   	end

   	return A, b, c, K
end

function solve_sedumi_jump(A, b, c, K, solver; kwargs...)
	# Solve Problems in Sedumi form, i.e.
	# minimize    c'x
	# subject to  Ax = b
	#             x \in K
	# See https://github.com/sqlp/sedumi/blob/master/sedumi.m#L49-L92

   	model = Model(with_optimizer(solver; kwargs...))

	# Fisrt handle all variables except semidefinite matrices
   	constraint_types = [MOI.Reals;
		MOI.Nonnegatives;
		fill(MOI.SecondOrderCone, length(K["q"]));
		fill(MOI.RotatedSecondOrderCone, length(K["r"]))]
   	dimensions = [K["f"]; K["l"]; K["q"]; K["r"]]

   	x = Vector{VariableRef}(undef, 0)
   	for (dimension, constraint_type) in zip(dimensions, constraint_types)
      		x_ = @variable(model, [1:dimension])
      		if constraint_type != MOI.Reals # Some solvers don't accept MOI.Reals
         			@constraint(model, x_ in constraint_type(dimension))
      		end
      		append!(x, x_)
   	end

	# Now handle semidefinite matrices
   	for dimension in K["s"]
      		x_ = @variable(model, [1:dimension, 1:dimension], PSD)
      		append!(x, reshape(x_, length(x_)))
   	end

   	@constraint(model, A * x .== b)
   	@objective(model, Min, dot(c, x))

   	JuMP.optimize!(model)
    return value.(x), JuMP.objective_value(model), JuMP.termination_status(model)
end

function solve_sedumi_jump_dual(A, b, c, K, solver; kwargs...)
	# Solve the dual of problem described in Sedumi form, i.e.
	# maximize    b'y
	# subject to  c - A'y \in K*
	# See https://github.com/sqlp/sedumi/blob/master/sedumi.m#L49-L92

   	model = Model(with_optimizer(solver; kwargs...))

	# Fisrt handle all variables except semidefinite matrices
   	constraint_types = [MOI.Zeros;
		MOI.Nonnegatives;
		fill(MOI.SecondOrderCone, length(K["q"]));
		fill(MOI.RotatedSecondOrderCone, length(K["r"]))]
   	dimensions = [K["f"]; K["l"]; K["q"]; K["r"]]

   	@variable(model, y[1:length(b)])
   	start_idx = 1
   	for (dimension, constraint_type) in zip(dimensions, constraint_types)
      		indices = start_idx:start_idx + dimension - 1
      		dimension > 0 && @constraint(model, c[indices] - A[:, indices]' * y in constraint_type(dimension))
      		start_idx += dimension
   	end
   	for dimension in K["s"]
      		indices = Vector(start_idx:start_idx + dimension^2 - 1)[upper_triangular_indices(dimension)]
      		@constraint(model, c[indices] - A[:, indices]' * y in MOI.PositiveSemidefiniteConeTriangle(dimension))
      		start_idx += dimension^2
   	end

   	@objective(model, Max, dot(b, y))

   	JuMP.optimize!(model)
    return value.(y), JuMP.objective_value(model), JuMP.termination_status(model)
end

function upper_triangular_indices(n)
   	indices = Vector{Int}(undef, Int(n * (n + 1) / 2))
   	index = 0
   	counter = 0
   	for j = 1:n, i = 1:n
      		index += 1
      		if i <= j
         			counter += 1
         			indices[counter] = index
      		end
   	end
   	return indices
end