text:
    catch_phrase_words = (
        (  "habilidad", "acceso", "adaptador", "algoritmo", "alianza",
           "analista", "aplicación", "enfoque", "arquitectura",
           "archivo", "inteligencia artificial", "array", "actitud",
           "medición", "gestión presupuestaria", "capacidad", "desafío",
           "circuito", "colaboración", "complejidad", "concepto",
           "conglomeración", "contingencia", "núcleo", "fidelidad",
           "base de datos", "data-warehouse", "definición", "emulación",
           "codificar", "encriptar", "extranet", "firmware",
           "flexibilidad", "focus group", "previsión", "base de trabajo",
           "función", "funcionalidad", "Interfaz Gráfica", "groupware",
           "Interfaz gráfico de usuario", "hardware", "Soporte", "jerarquía",
           "conjunto", "implementación", "infraestructura", "iniciativa",
           "instalación", "conjunto de instrucciones", "interfaz",
           "intranet", "base del conocimiento", "red de area local",
           "aprovechar", "matrices", "metodologías", "middleware",
           "migración", "modelo", "moderador", "monitorizar",
           "arquitectura abierta", "sistema abierto", "orquestar",
           "paradigma", "paralelismo", "política", "portal",
           "estructura de precios", "proceso de mejora",
           "producto", "productividad", "proyecto", "proyección",
           "protocolo", "línea segura", "software", "solución",
           "estandardización", "estrategia", "estructura", "éxito",
           "superestructura", "soporte", "sinergia", "mediante",
           "marco de tiempo", "caja de herramientas", "utilización",
           "website", "fuerza de trabajo"
           ),
            (
            "24 horas", "24/7", "3ra generación", "4ta generación",
            "5ta generación", "6ta generación", "analizada",
            "asimétrica", "asíncrona", "monitorizada por red",
            "bidireccional", "bifurcada", "generada por el cliente",
            "cliente servidor", "coherente", "cohesiva", "compuesto",
            "sensible al contexto", "basado en el contexto",
            "basado en contenido", "dedicada",
            "generado por la demanda", "didactica", "direccional",
            "discreta", "dinámica", "potenciada", "acompasada",
            "ejecutiva", "explícita", "tolerante a fallos",
            "innovadora", "amplio ábanico", "global", "heurística",
            "alto nivel", "holística", "homogénea", "híbrida",
            "incremental", "intangible", "interactiva", "intermedia",
            "local", "logística", "maximizada", "metódica",
            "misión crítica", "móbil", "modular", "motivadora",
            "multimedia", "multiestado", "multitarea", "nacional",
            "basado en necesidades", "neutral", "nueva generación",
            "no-volátil", "orientado a objetos", "óptima", "optimizada",
            "radical", "tiempo real", "recíproca", "regional",
            "escalable", "secundaria", "orientada a soluciones",
            "estable", "estatica", "sistemática", "sistémica",
            "tangible", "terciaria", "transicional", "uniforme",
            "valor añadido", "vía web", "defectos cero", "tolerancia cero"
        ),
        (
            "adaptivo", "avanzado", "asimilado", "automatizado",
            "balanceado", "enfocado al negocio",
            "centralizado", "clonado", "compatible", "configurable",
            "multiplataforma", "enfocado al cliente", "personalizable",
            "descentralizado", "digitizado", "distribuido", "diverso",
            "mejorado", "en toda la empresa", "ergonómico", "exclusivo",
            "expandido", "extendido", "cara a cara", "enfocado",
            "de primera línea", "totalmente configurable",
            "basado en funcionalidad", "fundamental", "horizontal",
            "implementado", "innovador", "integrado", "intuitivo",
            "inverso", "administrado", "mandatorio", "monitoreado",
            "multicanal", "multilateral", "multi-capas", "en red",
            "basado en objetos", "de arquitectura abierta",
            "Open-source", "operativo", "optimizado", "opcional",
            "orgánico", "organizado", "perseverante", "persistente",
            "polarizado", "preventivo", "proactivo", "enfocado a ganancias",
            "programable", "progresivo", "llave pública",
            "enfocado a la calidad", "reactivo", "realineado",
            "re-contextualizado", "reducido", "con ingeniería inversa",
            "de tamaño adecuado", "robusto", "seguro", "compartible",
            "sincronizado", "orientado a equipos", "total",
            "universal", "actualizable", "centrado al usuario",
            "versátil", "virtual", "visionario",
        )
    )

    bsWords = (
        (
            "implementa", "utiliza", "integrata", "optimiza",
            "evoluciona", "transforma", "abraza", "habilia",
            "orquesta", "reinventa", "agrega", "mejora", "incentiviza",
            "modifica", "empondera", "monetiza", "fortalece",
            "facilita", "synergiza",  "crear marca", "crece",
            "sintetiza", "entrega", "mezcla", "incuba", "compromete",
            "maximiza", "inmediata", "visualiza", "inova",
            "escala", "libera", "maneja", "extiende", "revoluciona",
            "genera", "explota", "transición", "itera", "cultiva",
            "redefine", "recontextualiza",
        ),
        (
            "synergías", "paradigmas", "marcados", "socios",
            "infraestructuras", "plataformas", "iniciativas",
            "chanales", "communidades", "ROI", "soluciones",
            "portales", "nichos", "tecnologías", "contenido",
            "cadena de producción", "convergencia", "relaciones",
            "architecturas", "interfaces", "comercio electrónico",
            "sistemas", "ancho de banda", "modelos", "entregables",
            "usuarios", "esquemas", "redes", "aplicaciones", "métricas",
            "funcionalidades", "experiencias", "servicios web",
            "metodologías"
        ),
        (
            "valor agregado", "verticales", "proactivas", "robustas",
            "revolucionarias", "escalables", "de punta", "innovadoras",
            "intuitivas", "estratégicas", "e-business", "de misión crítica",
            "uno-a-uno", "24/7", "end-to-end", "globales", "B2B", "B2C",
            "granulares", "sin fricciones", "virtuales", "virales",
            "dinámicas", "24/365", "magnéticas", "listo para la web",
            "interactivas", "dot-com", "sexi", "en tiempo real",
            "eficientes", "front-end", "distribuidas", "extensibles",
            "llave en mano", "de clase mundial", "open-source",
            "plataforma cruzada", "de paquete", "empresariales",
            "integrado", "impacto total", "inalámbrica", "transparentes",
            "de siguiente generación", "lo último", "centrado al usuario",
            "visionarias", "personalizado", "ubicuas", "plug-and-play",
            "colaborativas", "holísticas", "ricas"
        ),
    )

    company_preffixes = ("Despacho", "Grupo", "Corporativo", "Club",
                         "Industrias", "Laboratorios", "Proyectos")

    company_suffixes = ("A.C.", "S.A.", "S.A. de C.V.", "S.C.",
        "S. R.L. de C.V.","e Hijos", "y Asociados")

    company_suffix()= company_suffixes[rand(1:length(company_suffixes))]

    company_prefix() = company_preffixes[rand(1:length(company_preffixes))]

     catch_phrase()= (resultado ="";
        for word_list in catch_phrase_words
            resultado *= ("$(word_list[rand(1:length(word_list))]) ")
        end;
        resultado)

     bs()=(resultado ="";
        for word_list in bsWords
            resultado *= ("$(word_list[rand(1:length(word_list))]) ")
        end;
        resultado)

    company()=formatscomp[rand(1:length(formatscomp))]

    formatscomp = (
        "$(last_name()) $(company_suffix())",
        "$(last_name())-$(last_name())",
        "$(company_prefix()) $(last_name())-$(last_name())",
        "$(company_prefix()) $(last_name()) y $(last_name())",
        "$(company_prefix()) $(last_name()), $(last_name()) y $(last_name())",
        "$(last_name())-$(last_name()) $(company_suffix())",
        "$(last_name()), $(last_name()) y $(last_name())",
        "$(last_name()) y $(last_name()) $(company_suffix())"
    )


text:
struct Trace{F,R<:TracedRNG}
    f::F
    ctask::Libtask.CTask
    rng::R
end

const Particle = Trace

function Trace(f, rng::TracedRNG)
    ctask = let f = f
        Libtask.CTask() do
            res = f(rng)
            Libtask.produce(nothing)
            return res
        end
    end

    # add backward reference
    newtrace = Trace(f, ctask, rng)
    addreference!(ctask.task, newtrace)

    return newtrace
end

function Trace(f, ctask::Libtask.CTask)
    return Trace(f, ctask, TracedRNG())
end

# Copy task
Base.copy(trace::Trace) = Trace(trace.f, copy(trace.ctask), deepcopy(trace.rng))

# step to the next observe statement and
# return the log probability of the transition (or nothing if done)
function advance!(t::Trace, isref::Bool)
    isref ? load_state!(t.rng) : save_state!(t.rng)
    inc_counter!(t.rng)

    # Move to next step
    return Libtask.consume(t.ctask)
end

# reset log probability
reset_logprob!(t::Trace) = nothing

reset_model(f) = deepcopy(f)
delete_retained!(f) = nothing

# Task copying version of fork for Trace.
function fork(trace::Trace, isref::Bool=false)
    newtrace = copy(trace)
    isref && delete_retained!(newtrace.f)

    # add backward reference
    addreference!(newtrace.ctask.task, newtrace)

    return newtrace
end

# PG requires keeping all randomness for the reference particle
# Create new task and copy randomness
function forkr(trace::Trace)
    newf = reset_model(trace.f)
    Random123.set_counter!(trace.rng, 1)

    ctask = let f = trace.ctask.task.code
        Libtask.CTask() do
            res = f()(trace.rng)
            Libtask.produce(nothing)
            return res
        end
    end

    # add backward reference
    newtrace = Trace(newf, ctask, trace.rng)
    addreference!(ctask.task, newtrace)

    return newtrace
end

# create a backward reference in task_local_storage
function addreference!(task::Task, trace::Trace)
    if task.storage === nothing
        task.storage = IdDict()
    end
    task.storage[:__trace] = trace

    return task
end

current_trace() = current_task().storage[:__trace]

"""
Data structure for particle filters
- `effectiveSampleSize(pc :: ParticleContainer)`: Return the effective sample size of the particles in `pc`
"""
mutable struct ParticleContainer{T<:Particle,R<:TracedRNG}
    "Particles."
    vals::Vector{T}
    "Unnormalized logarithmic weights."
    logWs::Vector{Float64}
    "Traced RNG to replay the resampling step"
    rng::R
end

function ParticleContainer(particles::Vector{<:Particle})
    return ParticleContainer(particles, zeros(length(particles)), TracedRNG())
end

function ParticleContainer(particles::Vector{<:Particle}, r::TracedRNG)
    return ParticleContainer(particles, zeros(length(particles)), r)
end

Base.collect(pc::ParticleContainer) = pc.vals
Base.length(pc::ParticleContainer) = length(pc.vals)
Base.@propagate_inbounds Base.getindex(pc::ParticleContainer, i::Int) = pc.vals[i]

function Base.rand(rng::Random.AbstractRNG, pc::ParticleContainer)
    index = randcat(rng, getweights(pc))
    return pc[index]
end

# registers a new x-particle in the container
function Base.push!(pc::ParticleContainer, p::Particle)
    push!(pc.vals, p)
    push!(pc.logWs, 0.0)
    return pc
end

# clones a theta-particle
function Base.copy(pc::ParticleContainer)
    # fork particles
    vals = eltype(pc.vals)[fork(p) for p in pc.vals]

    # copy weights
    logWs = copy(pc.logWs)

    # Copy rng and states
    rng = copy(pc.rng)

    return ParticleContainer(vals, logWs, rng)
end

"""
    reset_logweights!(pc::ParticleContainer)

Reset all unnormalized logarithmic weights to zero.
"""
function reset_logweights!(pc::ParticleContainer)
    fill!(pc.logWs, 0.0)
    return pc
end

"""
    increase_logweight!(pc::ParticleContainer, i::Int, x)

Increase the unnormalized logarithmic weight of the `i`th particle with `x`.
"""
function increase_logweight!(pc::ParticleContainer, i, logw)
    pc.logWs[i] += logw
    return pc
end

"""
    getweights(pc::ParticleContainer)

Compute the normalized weights of the particles.
"""
getweights(pc::ParticleContainer) = StatsFuns.softmax(pc.logWs)

"""
    getweight(pc::ParticleContainer, i)

Compute the normalized weight of the `i`th particle.
"""
getweight(pc::ParticleContainer, i) = exp(pc.logWs[i] - logZ(pc))

"""
    logZ(pc::ParticleContainer)

Return the logarithm of the normalizing constant of the unnormalized logarithmic weights.
"""
logZ(pc::ParticleContainer) = StatsFuns.logsumexp(pc.logWs)

"""
    effectiveSampleSize(pc::ParticleContainer)

Compute the effective sample size ``1 / ∑ wᵢ²``, where ``wᵢ```are the normalized weights.
"""
function effectiveSampleSize(pc::ParticleContainer)
    Ws = getweights(pc)
    return inv(sum(abs2, Ws))
end

"""
    update_keys!(pc::ParticleContainer)

Create new unique keys for the particles in the ParticleContainer
"""
function update_keys!(pc::ParticleContainer, ref::Union{Particle,Nothing}=nothing)
    # Update keys to new particle ids
    nparticles = length(pc)
    n = ref === nothing ? nparticles : nparticles - 1
    for i in 1:n
        pi = pc.vals[i]
        k = split(pi.rng.rng.key)
        Random.seed!(pi.rng, k[1])
    end
    return nothing
end

"""
    resample_propagate!(rng, pc::ParticleContainer[, randcat = resample_systematic,
                        ref = nothing; weights = getweights(pc)])

Resample and propagate the particles in `pc`.

Function `randcat` is used for sampling ancestor indices from the categorical distribution
of the particle `weights`. For Particle Gibbs sampling, one can provide a reference particle
`ref` that is ensured to survive the resampling step.
"""
function resample_propagate!(
    rng::Random.AbstractRNG,
    pc::ParticleContainer,
    randcat=resample_systematic,
    ref::Union{Particle,Nothing}=nothing;
    weights=getweights(pc),
)
    # check that weights are not NaN
    @assert !any(isnan, weights)

    # sample ancestor indices
    n = length(pc)
    nresamples = ref === nothing ? n : n - 1
    indx = randcat(rng, weights, nresamples)

    # count number of children for each particle
    num_children = zeros(Int, n)
    @inbounds for i in indx
        num_children[i] += 1
    end

    # fork particles
    particles = collect(pc)
    children = similar(particles)
    j = 0
    @inbounds for i in 1:n
        ni = num_children[i]

        if ni > 0
            # fork first child
            pi = particles[i]
            isref = pi === ref
            p = isref ? fork(pi, isref) : pi
            nseeds = isref ? ni - 1 : ni

            seeds = split(p.rng.rng.key, nseeds)
            !isref && Random.seed!(p.rng, seeds[1])

            children[j += 1] = p
            # fork additional children
            for k in 2:ni
                part = fork(p, isref)
                Random.seed!(part.rng, seeds[k])
                children[j += 1] = part
            end
        end
    end

    if ref !== nothing
        # Insert the retained particle. This is based on the replaying trick for efficiency
        # reasons. If we implement PG using task copying, we need to store Nx * T particles!
        @inbounds children[n] = ref
    end

    # replace particles and log weights in the container with new particles and weights
    pc.vals = children
    reset_logweights!(pc)

    return pc
end

function resample_propagate!(
    rng::Random.AbstractRNG,
    pc::ParticleContainer,
    resampler::ResampleWithESSThreshold,
    ref::Union{Particle,Nothing}=nothing;
    weights=getweights(pc),
)
    # Compute the effective sample size ``1 / ∑ wᵢ²`` with normalized weights ``wᵢ``
    ess = inv(sum(abs2, weights))

    if ess ≤ resampler.threshold * length(pc)
        resample_propagate!(rng, pc, resampler.resampler, ref; weights=weights)
    else
        update_keys!(pc, ref)
    end

    return pc
end

"""
    reweight!(pc::ParticleContainer)

Check if the final time step is reached, and otherwise reweight the particles by
considering the next observation.
"""
function reweight!(pc::ParticleContainer, ref::Union{Particle,Nothing}=nothing)
    n = length(pc)

    particles = collect(pc)
    numdone = 0
    for i in 1:n
        p = particles[i]

        # Obtain ``\\log p(yₜ | y₁, …, yₜ₋₁, x₁, …, xₜ, θ₁, …, θₜ)``, or `nothing` if the
        # the execution of the model is finished.
        # Here ``yᵢ`` are observations, ``xᵢ`` variables of the particle filter, and
        # ``θᵢ`` are variables of other samplers.
        isref = p === ref
        score = advance!(p, isref)

        if score === nothing
            numdone += 1
        else
            # Increase the unnormalized logarithmic weights.
            increase_logweight!(pc, i, score)

            # Reset the accumulator of the log probability in the model so that we can
            # accumulate log probabilities of variables of other samplers until the next
            # observation.
            reset_logprob!(p)
        end
    end

    # Check if all particles are propagated to the final time point.
    numdone == n && return true

    # The posterior for models with random number of observations is not well-defined.
    if numdone != 0
        error(
            "mis-aligned execution traces: # particles = ",
            n,
            " # completed trajectories = ",
            numdone,
            ". Please make sure the number of observations is NOT random.",
        )
    end

    return false
end

"""
    sweep!(rng, pc::ParticleContainer, resampler)

Perform a particle sweep and return an unbiased estimate of the log evidence.

The resampling steps use the given `resampler`.

# Reference

Del Moral, P., Doucet, A., & Jasra, A. (2006). Sequential monte carlo samplers.
Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(3), 411-436.
"""
function sweep!(
    rng::Random.AbstractRNG,
    pc::ParticleContainer,
    resampler,
    ref::Union{Particle,Nothing}=nothing,
)
    # Initial step:
    # Resample and propagate particles.
    resample_propagate!(rng, pc, resampler, ref)

    # Compute the current normalizing constant ``Z₀`` of the unnormalized logarithmic
    # weights.
    # Usually it is equal to the number of particles in the beginning but this
    # implementation covers also the unlikely case of a particle container that is
    # initialized with non-zero logarithmic weights.
    logZ0 = logZ(pc)

    # Reweight the particles by including the first observation ``y₁``.
    isdone = reweight!(pc, ref)

    # Compute the normalizing constant ``Z₁`` after reweighting.
    logZ1 = logZ(pc)

    # Compute the estimate of the log evidence ``\\log p(y₁)``.
    logevidence = logZ1 - logZ0

    # For observations ``y₂, …, yₜ``:
    while !isdone
        # Resample and propagate particles.
        resample_propagate!(rng, pc, resampler, ref)

        # Compute the current normalizing constant ``Z₀`` of the unnormalized logarithmic
        # weights.
        logZ0 = logZ(pc)

        # Reweight the particles by including the next observation ``yₜ``.
        isdone = reweight!(pc, ref)

        # Compute the normalizing constant ``Z₁`` after reweighting.
        logZ1 = logZ(pc)

        # Compute the estimate of the log evidence ``\\log p(y₁, …, yₜ)``.
        logevidence += logZ1 - logZ0
    end

    return logevidence
end

ext: jl


text:
export add_ite, add_ite_const, add_eval_const
export add_leq
export add_comp

function add_ite(mgr::Ptr{Manager}, f::Ptr{Node}, g::Ptr{Node}, h::Ptr{Node})
    result = Cudd_addIte(mgr, f, g, h)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end

function add_ite_const(mgr::Ptr{Manager}, f::Ptr{Node}, g::Ptr{Node}, h::Ptr{Node})
    result = Cudd_addIteConstant(mgr, f, g, h)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end

function add_eval_const(mgr::Ptr{Manager}, f::Ptr{Node}, g::Ptr{Node})
    result = Cudd_addEvalConst(mgr, f, g)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end

function add_leq(mgr::Ptr{Manager}, f::Ptr{Node}, g::Ptr{Node})
    result = Cudd_addLeq(mgr, f, g)
    return result
end

function add_comp(mgr::Ptr{Manager}, f::Ptr{Node})
    result = Cudd_addCmpl(mgr, f)
    if result == C_NULL
        throw(OutOfMemoryError())
    end
    return result
end


text:
#=
Interpolation Methods

=#

using OffsetArrays: OffsetVector

"""
private function, find index in array where x would be inserted for interpolation
between xvec[i] and xvec[i+1]
"""
function findindex(xvec, x)

    n = length(xvec)
    i = searchsortedlast(xvec, x)

    # this version allows extrapolation
    if i == 0
        i = 1
    elseif i == n
        i = n - 1
    end

    # this version prevents extrapolation
    # if i == 0
    #     throw(DomainError(x, "x falls below range of provided spline data"))
    # elseif i == n
    #     if x == xvec[n]
    #         i = n - 1
    #     else
    #         throw(DomainError(x, "x falls above range of provided spline data"))
    #     end
    # end

    return i
end

## ------ Akima Interpolation  ---------

# struct used internally
struct Akima{TX, TY, TCoeff}

    xdata::TX
    ydata::TY
    p0::Vector{TCoeff}
    p1::Vector{TCoeff}
    p2::Vector{TCoeff}
    p3::Vector{TCoeff}

end

"""
    Akima(xdata, ydata, delta_x=0.0)

Creates an akima spline at node points: `xdata`, `ydata`.  This is a 1D spline that avoids
overshooting issues common with many other polynomial splines resulting in a more
natural curve.  It also only depends on local points (`i-2`...`i+2`) allow for more efficient
computation.  `delta_x` is the half width of a smoothing interval used for the absolute
value function.  Set `delta_x=0` to recover the original akima spline.  The smoothing
is only useful if you want to differentiate xdata and ydata.  In many case the nodal points
are fixed so this is not needed.  Returns an akima spline object (Akima struct).
This function, only performs construction of the spline, not evaluation.
This is useful if you want to evaluate the same mesh at multiple different conditions.
A convenience method exists below to perform both in one shot.
"""
function Akima(xdata::AbstractVector{TFX}, ydata::AbstractVector{TFY}, delta_x=0.0) where {TFX, TFY}

    # setup
    eps = 1e-30
    n = length(xdata)
    TCoeff = promote_type(TFX, TFY)

    # compute segment slopes
    m = OffsetVector(zeros(TCoeff, n+3), -1:n+1)
    for i = 1:n-1
        m[i] = (ydata[i+1] - ydata[i]) / (xdata[i+1] - xdata[i])
    end

    # estimation for end points
    m[0] = 2.0*m[1] - m[2]
    m[-1] = 2.0*m[0] - m[1]
    m[n] = 2.0*m[n-1] - m[n-2]
    m[n+1] = 2.0*m[n] - m[n-1]

    # slope at points
    t = zeros(TCoeff, n)
    for i = 1:n
        m1 = m[i-2]
        m2 = m[i-1]
        m3 = m[i]
        m4 = m[i+1]
        w1 = abs_smooth(m4 - m3, delta_x)
        w2 = abs_smooth(m2 - m1, delta_x)
        if (w1 < eps && w2 < eps)
            t[i] = 0.5*(m2 + m3)  # special case to avoid divide by zero
        else
            t[i] = (w1*m2 + w2*m3) / (w1 + w2)
        end
    end

    # polynomial cofficients
    p0 = zeros(TCoeff, n-1)
    p1 = zeros(TCoeff, n-1)
    p2 = zeros(TCoeff, n-1)
    p3 = zeros(TCoeff, n-1)
    for i = 1:n-1
        dx = xdata[i+1] - xdata[i]
        t1 = t[i]
        t2 = t[i+1]
        p0[i] = ydata[i]
        p1[i] = t1
        p2[i] = (3.0*m[i] - 2.0*t1 - t2)/dx
        p3[i] = (t1 + t2 - 2.0*m[i])/dx^2
    end

    return Akima(xdata, ydata, p0, p1, p2, p3)
end

function (spline::Akima)(x::Number)

    j = findindex(spline.xdata, x)

    # evaluate polynomial
    dx = x - spline.xdata[j]
    y = spline.p0[j] + spline.p1[j]*dx + spline.p2[j]*dx^2 + spline.p3[j]*dx^3

    return y
end

(spline::Akima)(x::AbstractVector) = spline.(x)

"""
    akima(x, y, xpt)

A convenience method to perform construction and evaluation of the spline in one step.
See docstring for Akima for more details.

**Arguments**
- `x, y::Vector{Float}`: the node points
- `xpt::Vector{Float} or ::Float`: the evaluation point(s)

**Returns**
- `ypt::Vector{Float} or ::Float`: interpolated value(s) at xpt using akima spline.
"""
akima(x, y, xpt) = Akima(x, y)(xpt)

"""
    derivative(spline, x)

Computes the derivative of an Akima spline at x.

**Arguments**
- `spline::Akima}`: an Akima spline
- `x::Float`: the evaluation point(s)

**Returns**
- `dydx::Float`: derivative at x using akima spline.
"""
function derivative(spline::Akima, x)

    j = findindex(spline.xdata, x)

    # evaluate polynomial
    dx = x - spline.xdata[j]
    dydx = spline.p1[j] + 2*spline.p2[j]*dx + 3*spline.p3[j]*dx^2

    return dydx
end

"""
    gradient(spline, x)

Computes the gradient of a Akima spline at x.

**Arguments**
- `spline::Akima}`: an Akima spline
- `x::Vector{Float}`: the evaluation point(s)

**Returns**
- `dydx::Vector{Float}`: gradient at x using akima spline.
"""
gradient(spline::Akima, x) = derivative.(Ref(spline), x)

# ------ Linear Interpolation  ---------

"""
    linear(xdata, ydata, x::Number)

Linear interpolation.

**Arguments**
- `xdata::Vector{Float64}`: x data used in constructing interpolation
- `ydata::Vector{Float64}`: y data used in constructing interpolation
- `x::Float64`: point to evaluate spline at

**Returns**
- `y::Float64`: value at x using linear interpolation
"""
function linear(xdata, ydata, x::Number)

    i = findindex(xdata, x)

    # lienar interpolation
    eta = (x - xdata[i]) / (xdata[i+1] - xdata[i])
    y = ydata[i] + eta*(ydata[i+1] - ydata[i])

    return y
end

"""
    linear(xdata, ydata, x::AbstractVector)

Convenience function to perform linear interpolation at multiple points.
"""
linear(xdata, ydata, x::AbstractVector) = linear.(Ref(xdata), Ref(ydata), x)

"""
derivative of linear interpolation at `x::Number`
"""
function derivative(xdata, ydata, x)

    i = findindex(xdata, x)
    dydx = (ydata[i+1] - ydata[i]) / (xdata[i+1] - xdata[i])

    return dydx
end

"""
gradient of linear interpolation at `x::Vector`
"""
gradient(xdata, ydata, x) = derivative.(Ref(xdata), Ref(ydata), x)

# ------- higher order recursive 1D interpolation ------

"""
   interp2d(interp1d, xdata, ydata, fdata, xpt, ypt)

2D interpolation using recursive 1D interpolation.  This approach is likely less efficient than a more direct
2D interpolation method, especially one you can create separate creation from evaluation,
but it is generalizable to any spline approach and any dimension.

**Arguments**
- `interp1d`: any spline function of form: ypt = interp1d(xdata, ydata, xpt) where data are the known
    data(node) points and pt are the points where you want to evaluate the spline at.
- `xdata::Vector{Float}`, `ydata::Vector{Float}`: Define the 2D grid
- `fdata::Matrix{Float}`: where fdata[i, j] is the function value at xdata[i], ydata[j]
- `xpt::Vector{Float}`, `ypt::Vector{Float}`: the locations where you want to evaluate the spline

**Returns**
- `fhat::Matrix{Float}`: where fhat[i, j] is the estimate function value at xpt[i], ypt[j]
"""
function interp2d(interp1d, xdata, ydata, fdata, xpt, ypt)

    ny = length(ydata)
    nxpt = length(xpt)
    nypt = length(ypt)

    R = promote_type(eltype(xpt), eltype(ypt))

    yinterp = Array{R}(undef, ny, nxpt)
    output = Array{R}(undef, nxpt, nypt)

    for i = 1:ny
        yinterp[i, :] .= interp1d(xdata, fdata[:, i], xpt)
    end
    for i = 1:nxpt
        output[i, :] .= interp1d(ydata, yinterp[:, i], ypt)
    end

    return output
end

"""
    interp3d(interp1d, xdata, ydata, zdata, fdata, xpt, ypt, zpt)

Same as interp2d, except in three dimension.
"""
function interp3d(interp1d, xdata, ydata, zdata, fdata, xpt, ypt, zpt)

    nz = length(zdata)
    nxpt = length(xpt)
    nypt = length(ypt)
    nzpt = length(zpt)

    R = promote_type(eltype(xpt), eltype(ypt), eltype(zpt))
    zinterp = Array{R}(undef, nz, nxpt, nypt)
    output = Array{R}(undef, nxpt, nypt, nzpt)

    for i = 1:nz
        zinterp[i, :, :] .= interp2d(interp1d, xdata, ydata, fdata[:, :, i], xpt, ypt)
    end
    for j = 1:nypt
        for i = 1:nxpt
            output[i, j, :] .= interp1d(zdata, zinterp[:, i, j], zpt)
        end
    end

    return output
end

"""
    interp4d(interp1d, xdata, ydata, zdata, fdata, xpt, ypt, zpt)

Same as interp3d, except in four dimensions.
"""
function interp4d(interp1d, xdata, ydata, zdata, tdata, fdata, xpt, ypt, zpt, tpt)

    nt = length(tdata)
    nxpt = length(xpt)
    nypt = length(ypt)
    nzpt = length(zpt)
    ntpt = length(tpt)

    R = promote_type(eltype(xpt), eltype(ypt), eltype(zpt), eltype(tpt))
    tinterp = Array{R}(undef, nt, nxpt, nypt, nzpt)
    output = Array{R}(undef, nxpt, nypt, nzpt, ntpt)

    for i = 1:nt
        tinterp[i, :, :, :] .= interp3d(interp1d, xdata, ydata, zdata, fdata[:, :, :, i], xpt, ypt, zpt)
    end
    for k = 1:nzpt
        for j = 1:nypt
            for i = 1:nxpt
                output[i, j, k, :] .= interp1d(tdata, tinterp[:, i, j, k], tpt)
            end
        end
    end

    return output
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -47.221875796417265 <= q <= 49.11122836315772)

                     add_NL_constraint(m, :(log(1 + exp(0.42790162605340454 + 0.9948146339173967*log(1 + exp(0.30429674417798 + 0.01326158862793747*log(1 + exp(0.7020479917236249 + 0.5157446936496526*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.8162922098274223*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6507571114046486*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + 0.3903970834590602*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.8506753528903843*log(1 + exp(0.7594455380381384 + 0.9088497207479063*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + 0.34867162336420376*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + -0.586942758624089*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.05285512797069236*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.7699605184902198*log(1 + exp(0.5333438081093944 + -0.9884706408297048*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.9213430875393112*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6584439947423988*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.39397409418232554*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.9481609482474158*log(1 + exp(-0.8688273038537231 + -0.5762764679031651*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.4115583086286181*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.07598378122288318*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.5005567670062852*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))))) + 0.5690258167368487*log(1 + exp(0.4003333530904172 + 0.4565361301167439*log(1 + exp(0.7020479917236249 + 0.5157446936496526*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.8162922098274223*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6507571114046486*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + 0.3903970834590602*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + -0.41218949132242644*log(1 + exp(0.7594455380381384 + 0.9088497207479063*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + 0.34867162336420376*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + -0.586942758624089*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.05285512797069236*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.6719961553691296*log(1 + exp(0.5333438081093944 + -0.9884706408297048*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.9213430875393112*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6584439947423988*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.39397409418232554*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.9454422589116072*log(1 + exp(-0.8688273038537231 + -0.5762764679031651*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.4115583086286181*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.07598378122288318*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.5005567670062852*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))))) + -0.632632285898699*log(1 + exp(-0.4740061987673023 + -0.21503151130953757*log(1 + exp(0.7020479917236249 + 0.5157446936496526*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.8162922098274223*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6507571114046486*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + 0.3903970834590602*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.7641417558802628*log(1 + exp(0.7594455380381384 + 0.9088497207479063*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + 0.34867162336420376*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + -0.586942758624089*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.05285512797069236*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.8796988576372375*log(1 + exp(0.5333438081093944 + -0.9884706408297048*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.9213430875393112*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6584439947423988*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.39397409418232554*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.6466469730144713*log(1 + exp(-0.8688273038537231 + -0.5762764679031651*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.4115583086286181*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.07598378122288318*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.5005567670062852*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))))) + -0.0799939753544896*log(1 + exp(-0.9533597279323933 + 0.019939841305541606*log(1 + exp(0.7020479917236249 + 0.5157446936496526*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.8162922098274223*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6507571114046486*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + 0.3903970834590602*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + -0.8134733780415981*log(1 + exp(0.7594455380381384 + 0.9088497207479063*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + 0.34867162336420376*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + -0.586942758624089*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.05285512797069236*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + -0.10348468441851777*log(1 + exp(0.5333438081093944 + -0.9884706408297048*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.9213430875393112*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6584439947423988*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.39397409418232554*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.008542712959329624*log(1 + exp(-0.8688273038537231 + -0.5762764679031651*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.4115583086286181*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.07598378122288318*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.5005567670062852*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))))))) + log(1 + exp(-0.31282642357163937 + 0.400776902433146*log(1 + exp(0.30429674417798 + 0.01326158862793747*log(1 + exp(0.7020479917236249 + 0.5157446936496526*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.8162922098274223*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6507571114046486*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + 0.3903970834590602*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.8506753528903843*log(1 + exp(0.7594455380381384 + 0.9088497207479063*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + 0.34867162336420376*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + -0.586942758624089*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.05285512797069236*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.7699605184902198*log(1 + exp(0.5333438081093944 + -0.9884706408297048*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.9213430875393112*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6584439947423988*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.39397409418232554*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.9481609482474158*log(1 + exp(-0.8688273038537231 + -0.5762764679031651*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.4115583086286181*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.07598378122288318*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.5005567670062852*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))))) + -0.8700113611274283*log(1 + exp(0.4003333530904172 + 0.4565361301167439*log(1 + exp(0.7020479917236249 + 0.5157446936496526*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.8162922098274223*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6507571114046486*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + 0.3903970834590602*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + -0.41218949132242644*log(1 + exp(0.7594455380381384 + 0.9088497207479063*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + 0.34867162336420376*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + -0.586942758624089*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.05285512797069236*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.6719961553691296*log(1 + exp(0.5333438081093944 + -0.9884706408297048*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.9213430875393112*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6584439947423988*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.39397409418232554*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.9454422589116072*log(1 + exp(-0.8688273038537231 + -0.5762764679031651*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.4115583086286181*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.07598378122288318*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.5005567670062852*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))))) + -0.8173535730312884*log(1 + exp(-0.4740061987673023 + -0.21503151130953757*log(1 + exp(0.7020479917236249 + 0.5157446936496526*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.8162922098274223*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6507571114046486*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + 0.3903970834590602*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.7641417558802628*log(1 + exp(0.7594455380381384 + 0.9088497207479063*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + 0.34867162336420376*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + -0.586942758624089*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.05285512797069236*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.8796988576372375*log(1 + exp(0.5333438081093944 + -0.9884706408297048*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.9213430875393112*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6584439947423988*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.39397409418232554*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.6466469730144713*log(1 + exp(-0.8688273038537231 + -0.5762764679031651*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.4115583086286181*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.07598378122288318*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.5005567670062852*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))))) + -0.10383649475667989*log(1 + exp(-0.9533597279323933 + 0.019939841305541606*log(1 + exp(0.7020479917236249 + 0.5157446936496526*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.8162922098274223*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6507571114046486*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + 0.3903970834590602*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + -0.8134733780415981*log(1 + exp(0.7594455380381384 + 0.9088497207479063*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + 0.34867162336420376*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + -0.586942758624089*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.05285512797069236*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + -0.10348468441851777*log(1 + exp(0.5333438081093944 + -0.9884706408297048*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.9213430875393112*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6584439947423988*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.39397409418232554*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.008542712959329624*log(1 + exp(-0.8688273038537231 + -0.5762764679031651*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.4115583086286181*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.07598378122288318*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.5005567670062852*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))))))) + log(1 + exp(0.07748802701261592 + -0.48287575350847955*log(1 + exp(0.30429674417798 + 0.01326158862793747*log(1 + exp(0.7020479917236249 + 0.5157446936496526*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.8162922098274223*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6507571114046486*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + 0.3903970834590602*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.8506753528903843*log(1 + exp(0.7594455380381384 + 0.9088497207479063*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + 0.34867162336420376*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + -0.586942758624089*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.05285512797069236*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.7699605184902198*log(1 + exp(0.5333438081093944 + -0.9884706408297048*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.9213430875393112*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6584439947423988*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.39397409418232554*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.9481609482474158*log(1 + exp(-0.8688273038537231 + -0.5762764679031651*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.4115583086286181*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.07598378122288318*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.5005567670062852*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))))) + 0.9032398979921337*log(1 + exp(0.4003333530904172 + 0.4565361301167439*log(1 + exp(0.7020479917236249 + 0.5157446936496526*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.8162922098274223*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6507571114046486*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + 0.3903970834590602*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + -0.41218949132242644*log(1 + exp(0.7594455380381384 + 0.9088497207479063*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + 0.34867162336420376*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + -0.586942758624089*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.05285512797069236*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.6719961553691296*log(1 + exp(0.5333438081093944 + -0.9884706408297048*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.9213430875393112*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6584439947423988*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.39397409418232554*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.9454422589116072*log(1 + exp(-0.8688273038537231 + -0.5762764679031651*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.4115583086286181*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.07598378122288318*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.5005567670062852*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))))) + 0.8154202752057489*log(1 + exp(-0.4740061987673023 + -0.21503151130953757*log(1 + exp(0.7020479917236249 + 0.5157446936496526*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.8162922098274223*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6507571114046486*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + 0.3903970834590602*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.7641417558802628*log(1 + exp(0.7594455380381384 + 0.9088497207479063*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + 0.34867162336420376*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + -0.586942758624089*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.05285512797069236*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.8796988576372375*log(1 + exp(0.5333438081093944 + -0.9884706408297048*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.9213430875393112*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6584439947423988*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.39397409418232554*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.6466469730144713*log(1 + exp(-0.8688273038537231 + -0.5762764679031651*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.4115583086286181*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.07598378122288318*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.5005567670062852*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))))) + 0.8885160528665201*log(1 + exp(-0.9533597279323933 + 0.019939841305541606*log(1 + exp(0.7020479917236249 + 0.5157446936496526*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.8162922098274223*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6507571114046486*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + 0.3903970834590602*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + -0.8134733780415981*log(1 + exp(0.7594455380381384 + 0.9088497207479063*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + 0.34867162336420376*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + -0.586942758624089*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.05285512797069236*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + -0.10348468441851777*log(1 + exp(0.5333438081093944 + -0.9884706408297048*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.9213430875393112*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6584439947423988*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.39397409418232554*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.008542712959329624*log(1 + exp(-0.8688273038537231 + -0.5762764679031651*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.4115583086286181*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.07598378122288318*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.5005567670062852*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))))))) + log(1 + exp(0.6389491685818425 + -0.9024035564379238*log(1 + exp(0.30429674417798 + 0.01326158862793747*log(1 + exp(0.7020479917236249 + 0.5157446936496526*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.8162922098274223*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6507571114046486*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + 0.3903970834590602*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.8506753528903843*log(1 + exp(0.7594455380381384 + 0.9088497207479063*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + 0.34867162336420376*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + -0.586942758624089*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.05285512797069236*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.7699605184902198*log(1 + exp(0.5333438081093944 + -0.9884706408297048*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.9213430875393112*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6584439947423988*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.39397409418232554*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.9481609482474158*log(1 + exp(-0.8688273038537231 + -0.5762764679031651*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.4115583086286181*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.07598378122288318*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.5005567670062852*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))))) + -0.8688708754259227*log(1 + exp(0.4003333530904172 + 0.4565361301167439*log(1 + exp(0.7020479917236249 + 0.5157446936496526*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.8162922098274223*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6507571114046486*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + 0.3903970834590602*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + -0.41218949132242644*log(1 + exp(0.7594455380381384 + 0.9088497207479063*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + 0.34867162336420376*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + -0.586942758624089*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.05285512797069236*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.6719961553691296*log(1 + exp(0.5333438081093944 + -0.9884706408297048*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.9213430875393112*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6584439947423988*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.39397409418232554*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.9454422589116072*log(1 + exp(-0.8688273038537231 + -0.5762764679031651*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.4115583086286181*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.07598378122288318*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.5005567670062852*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))))) + -0.9399903753655909*log(1 + exp(-0.4740061987673023 + -0.21503151130953757*log(1 + exp(0.7020479917236249 + 0.5157446936496526*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.8162922098274223*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6507571114046486*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + 0.3903970834590602*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.7641417558802628*log(1 + exp(0.7594455380381384 + 0.9088497207479063*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + 0.34867162336420376*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + -0.586942758624089*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.05285512797069236*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.8796988576372375*log(1 + exp(0.5333438081093944 + -0.9884706408297048*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.9213430875393112*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6584439947423988*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.39397409418232554*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.6466469730144713*log(1 + exp(-0.8688273038537231 + -0.5762764679031651*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.4115583086286181*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.07598378122288318*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.5005567670062852*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))))) + 0.5656846187923494*log(1 + exp(-0.9533597279323933 + 0.019939841305541606*log(1 + exp(0.7020479917236249 + 0.5157446936496526*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.8162922098274223*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6507571114046486*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + 0.3903970834590602*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + -0.8134733780415981*log(1 + exp(0.7594455380381384 + 0.9088497207479063*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + 0.34867162336420376*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + -0.586942758624089*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.05285512797069236*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + -0.10348468441851777*log(1 + exp(0.5333438081093944 + -0.9884706408297048*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.9213430875393112*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.6584439947423988*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.39397409418232554*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))) + 0.008542712959329624*log(1 + exp(-0.8688273038537231 + -0.5762764679031651*log(1 + exp(-0.990780872931377 + 0.7098282678057655*$(x[1]) + -0.9123645360112271*$(x[2]))) + -0.4115583086286181*log(1 + exp(-0.26001558800187974 + -0.2578347550825866*$(x[1]) + -0.8440952761970419*$(x[2]))) + 0.07598378122288318*log(1 + exp(0.6359235954642828 + 0.0835591752979501*$(x[1]) + -0.4064923799931912*$(x[2]))) + -0.5005567670062852*log(1 + exp(0.919513762198326 + -0.48206842436164976*$(x[1]) + 0.2512150299638187*$(x[2]))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module DifferentialForms

using Reexport

# The order of these include statements matters
include("Defs.jl")
include("Forms.jl")
include("Multivectors.jl")

# @reexport using .Defs
using .Defs
export unit, hodge, ⋆, bitsign, sort_perm
@reexport using .Forms
@reexport using .Multivectors

end


text:
"""
    DiGraph{V,E}()
    DiGraph{V}()
    DiGraph()

Construct a static directed graph with node type `V` and edge type `E`. `DiGraph()` constructs a static directed graph with integer nodes and edges.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> g = DiGraph()
DiGraph 0 nodes, 0 edges

julia> add_edge!(g, 1, 2)
Node(1)->Node(2)

julia> add_edge!(g, 2, 3)
Node(2)->Node(3)

julia> adjacency_matrix(g)
3×3 Array{Float64,2}:
 0.0  1.0  0.0
 0.0  0.0  1.0
 0.0  0.0  0.0

julia> add_node!(g, 4)
Node(4)

julia> nodes(g)
4-element Array{EvolvingGraphs.Node{Int64},1}:
 Node(1)
 Node(2)
 Node(3)
 Node(4)
```
"""
mutable struct DiGraph{V<:AbstractNode, E<:AbstractEdge, IncList, VK} <: AbstractStaticGraph{V, E}
    nodes::Vector{V}
    edges::Vector{E}
    forward_ilist::IncList # forward incidence list
    backward_ilist::IncList # backward incidence list
    indexof::Dict{VK,Int}
end
DiGraph{V, E}() where V where E = DiGraph(V[], E[], [], [], Dict{eltype(V), Int}())
DiGraph{V}() where V = DiGraph{V, Edge{V}}()
DiGraph() = DiGraph{Node{Int}, Edge{Node{Int}}}()
is_directed(g::DiGraph) = true

nodes(g::DiGraph) = g.nodes
num_nodes(g::DiGraph) = length(g.nodes)
edges(g::DiGraph) = g.edges
num_edges(g::DiGraph) = length(g.edges)

function add_bunch_of_edges!(g::DiGraph, ebunch)
    for e in ebunch
        if length(e) == 2
            i,j = e
            add_node!(g,i)
            add_node!(g,j)
            add_edge!(g,i,j)
        else
            error("Each edge in ebunch must have 2 elemnts: source and target")
        end
    end
    g
end

"""
    out_edges(g, v)

Return the outward edges of node `v` in a static graph `g`, where `v` can be a node or a node key.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> g = DiGraph()
DiGraph 0 nodes, 0 edges

julia> add_bunch_of_edges!(g, [(1,2), (2,3), (3,4)])
DiGraph 4 nodes, 3 edges

julia> out_edges(g, 1)
1-element Array{EvolvingGraphs.Edge{EvolvingGraphs.Node{Int64}},1}:
 Node(1)->Node(2)
```
"""
out_edges{V}(g::DiGraph{V}, v::V) = g.forward_ilist[node_index(v)]
out_edges{V}(g::DiGraph{V}, key::eltype(V)) = (v = Node(g.indexof[key],key); out_edges(g, v))

"""
    out_degree(g, v)

Return the number of outward edges of node `v` in static graph `g`, where `v` can be a node of a node key.
"""
out_degree{V}(g::DiGraph{V}, v::V) = length(out_edges(g, v))
out_degree{V}(g::DiGraph{V}, key::eltype(V)) = (v = Node(g.indexof[key],key); out_degree(g, v))

"""
    in_edges(g, v)

Return the inward edges of node `v` in static graph `g`.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> g = DiGraph()
DiGraph 0 nodes, 0 edges

julia> add_bunch_of_edges!(g, [(2,1), (3,1), (4,1)])
DiGraph 4 nodes, 3 edges

julia> in_edges(g,1)
3-element Array{EvolvingGraphs.Edge{EvolvingGraphs.Node{Int64}},1}:
 Node(2)->Node(1)
 Node(3)->Node(1)
 Node(4)->Node(1)

julia> in_degree(g,1)
3
```
"""
in_edges{V}(g::DiGraph{V}, v::V) = g.backward_ilist[node_index(v)]
in_edges{V}(g::DiGraph{V}, key::eltype(V)) = (v = Node(g.indexof[key],key); in_edges(g, v))

"""
    in_degree(g, v)

Return the number of inward edges of node `v` in static graph `g`.
"""
in_degree{V}(g::DiGraph{V}, v::V) = length(in_edges(g, v))
in_degree{V}(g::DiGraph{V}, key::eltype(V)) = (v = Node(g.indexof[key],key); in_degree(g, v))

function add_node!{K,E}(g::DiGraph{Node{K},E}, key::K; return_graph::Bool = false)
    index = get(g.indexof, key, 0)
    if index != 0
        v = Node(index,key)
        return return_graph? (g, v) : v
    end

    v = Node{K}(g, key)
    push!(g.nodes, v)
    push!(g.forward_ilist, E[])
    push!(g.backward_ilist, E[])
    g.indexof[key] = node_index(v)
    return_graph ? (g,v) : v
end

function add_edge!{V, E}(g::DiGraph{V ,E}, u::V, v::V, e::E)
    ui = node_index(u)
    vi = node_index(v)
    if !(e in g.forward_ilist[ui])
        push!(g.forward_ilist[ui], e)
    end
    if !(e in g.backward_ilist[vi])
        push!(g.backward_ilist[vi], e)
    end
    if !(e in g.edges)
        push!(g.edges, e)
    end
    e
end

function add_edge!{V,E}(g::DiGraph{V,E}, u::V, v::V)
    e = E(u, v)
    add_edge!(g, u, v, e)
end

function add_edge!{K}(g::DiGraph{Node{K}}, u_key::K, v_key::K)
    g, u = add_node!(g, u_key, return_graph = true)
    g, v = add_node!(g, v_key, return_graph = true)

    add_edge!(g, u, v)
end

"""
    aggregate_graph(g)

Convert an evolving graph `g` to a static graph by aggregate the edges to the same timestamp.

# Example

```jldoctest
julia> using EvolvingGraphs

julia> eg = EvolvingGraph()
Directed EvolvingGraph 0 nodes, 0 static edges, 0 timestamps

julia> add_bunch_of_edges!(eg, [(1,2,3), (2,3,4), (1,2,1)])
Directed EvolvingGraph 3 nodes, 3 static edges, 3 timestamps

julia> aggregate_graph(eg)
DiGraph 3 nodes, 2 edges
```
"""
function aggregate_graph{V}(eg::AbstractEvolvingGraph{V})
    g = DiGraph{V,Edge{V}}()
    for e in edges(eg)
        add_edge!(g, node_key(source(e)), node_key(target(e)))
    end
    g
end

find_node{V}(g::DiGraph{V}, v::V) = v in g.nodes ? v : false
find_node(g::DiGraph, key) = get(g.indexof, key, false)

function adjacency_matrix(g::DiGraph)
    n = num_nodes(g)
    es = edges(g)
    A = zeros(Float64, n, n)
    for e in es
        i = node_index(source(e))
        j = node_index(target(e))
        A[i,j] = 1.
    end
    return A
end

function sparse_adjacency_matrix(g::DiGraph)
    es = edges(g)
    n = num_nodes(g)
    is = Int[]
    js = Int[]

    for e in es
        i = node_index(source(e))
        j = node_index(target(e))
        push!(is, i)
        push!(js, j)
    end
    vs = ones(Float64, length(is))
    return sparse(is, js, vs, n, n)
end


text:
"""
Run some analyses on the wine dataset
"""

import DelimitedFiles

const wine_data_url = "https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data"
const wine_data_file = download(wine_data_url)

const raw_wine_data = DelimitedFiles.readdlm(wine_data_file,',', Float64)
const wine_quant = raw_wine_data[:,2:end]
const wine_labels = Int.(raw_wine_data[:,1])

Plots.scatter(wine_quant[:,1], wine_quant[:,2], group = wine_labels)
Plots.title!("Basic wine scatter plot")
Plots.savefig(joinpath(base_img, "basic_wine_scatter.pdf"))

# Fit a MLE on two first variables, second log-transformed
const res_normal = fit_mle(MvNormal, [wine_quant[:,1] log.(wine_quant[:,2])]')

# function for contour plot
function cont_func(x1, x2)
    pdf(res_normal, [x1,x2])
end

const wine_kde = KernelDensity.kde((wine_quant[:,1], log.(wine_quant[:,2])))

# contour plot for KDE
function cont_kde(x1, x2)
   pdf(wine_kde, x1, x2)
end

# local scope not to add plots to global one
let
    p = Plots.contour(11.0:0.05:15.0, -0.5:0.05:2.5, cont_func)
    Plots.scatter!(p, wine_quant[:,1], log.(wine_quant[:,2]), label = "Data points")
    Plots.title!(p, "Wine scatter plot & Gaussian maximum likelihood estimation")
    Plots.savefig(p, joinpath(base_img, "gaussian_contour_mle.pdf"))

    p = Plots.contour(11.0:0.05:15.0, -0.5:0.05:2.5, cont_kde)
    Plots.scatter!(p, wine_quant[:,1], log.(wine_quant[:,2]), group = wine_labels)
    Plots.title!(p, "Wine scatter plot & Kernel Density Estimation")
    Plots.savefig(p, joinpath(base_img, "gaussian_contour_kde.pdf"))
end

# fit a (Normal × LogNormal) on x1 × x2

function build_product_distribution(p)
    return Product([
        Normal(p[1], p[2]),
        LogNormal(p[3], p[4]),
    ])
end

function loglike(p)
    d = build_product_distribution(p)
    return loglikelihood(d.v[1], wine_quant[:,1]) + loglikelihood(d.v[2], wine_quant[:,2] .- 0.73)
end

∇L(p) = ForwardDiff.gradient(loglike, p)

let
    Random.seed!(42)
    # random start
    p = [10.0 + 3.0 * rand(), rand()+1, 2.0 + 3.0*rand(), rand()+1]
    iter = 1
    maxiter = 5000
    mean_paths = (Vector{Float64}(undef, maxiter), Vector{Float64}(undef, maxiter))
    while iter <= maxiter && sum(abs.(∇L(p))) >= 10^-6
        mean_paths[1][iter] = p[1]
        mean_paths[2][iter] = p[3]
        p = p + 0.05 * inv(iter+5) * ∇L(p)
        p[2] = p[2] < 0 ? -p[2] : p[2]
        p[4] = p[4] < 0 ? -p[4] : p[4]
        @info sum(abs.(∇L(p)))
        iter += 1
    end
    @info iter
    @info p
    d = build_product_distribution(p)
    cont_product(x1, x2) = pdf(d, [x1, x2])

    plt = Plots.contour(11.0:0.05:15.0, -0.5:0.05:2.5, cont_product, contours = false, legend = nothing)
    Plots.scatter!(plt, wine_quant[:,1], wine_quant[:,2] .- 0.73, legend = nothing)
    Plots.title!(plt, "Data points & final contour")
    Plots.xticks!(plt, 11.0:1.0:15.0)
    Plots.xaxis!(plt, (10.9,15.0))
    Plots.yaxis!(plt, (-0.5,5.2))
    plt1 = Plots.plot(collect(11.0:0.05:14.9), [pdf(d.v[1], x) for x in 11.0:0.05:14.9], legend = nothing, width=5, color = :blue)
    Plots.title!(plt1, "Normal marginal")
    Plots.xticks!(plt1, 11.0:1.0:15.0)
    Plots.xaxis!(plt1, (10.9,15.0))
    plt2 = Plots.plot([pdf(d.v[2], x) for x in 0.0:0.05:5.0], collect(0.0:0.05:5.0), legend = nothing, width=5, color = :red)
    Plots.yaxis!(plt2, (-0.5,5.2))
    Plots.title!(plt2, "Lognormal marginal")
    pfinal = Plots.plot(plt, plt2, plt1, Plots.plot(framestyle = :none), layout = (2,2), size = (1200,800))
    Plots.savefig(pfinal, joinpath(base_img, "wine_product_dist.pdf"))

    pltconv = Plots.scatter(wine_quant[:,1], wine_quant[:,2] .- 0.73, label = "Data points", markersize = 3)
    Plots.plot!(pltconv, mean_paths[1][1:iter-1], exp.(mean_paths[2][1:iter-1]), label = "", linestyle = :dot, linecolor = :black, alpha = 0.7)
    Plots.scatter!(pltconv, mean_paths[1][1:iter-1], exp.(mean_paths[2][1:iter-1]), label = "Mean iterations", marker = :cross, linecolor = :black, alpha = 0.8)
    Plots.scatter!(pltconv, [p[1]], [exp(p[3])], label = "Final mean", color = :red)
    Plots.title!(pltconv, "Likelihood maximization - $(iter-1) iterations")
    Plots.savefig(pltconv, joinpath(base_img, "wine_product_convergence.pdf"))
end

# reduced block for benchmark
function bench_func()
    # random start
    p = [10.0 + 3.0 * rand(), rand()+1, 2.0 + 3.0*rand(), rand()+1]
    iter = 1
    maxiter = 5000
    while iter <= maxiter && sum(abs.(∇L(p))) >= 10^-6
        p = p + 0.05 * inv(iter+5) * ∇L(p)
        p[2] = p[2] < 0 ? -p[2] : p[2]
        p[4] = p[4] < 0 ? -p[4] : p[4]
        iter += 1
    end
    return build_product_distribution(p)
end

let
    X = [wine_quant[:,1] log.(wine_quant[:,2])]
    n = size(X, 1)
    Nk = 3
    for maxiter in (0, 20, 22, 25, 30, 50, 60, 65, 67, 70, 80)
        (dists, prior, Z, l, iter) = expectation_maximization(MvNormal, X, Nk, maxiter)
        contour_mixture(x1, x2) = maximum([pdf(dists[k], [x1,x2]) for k in 1:Nk])
        colors = map(1:n) do i
            (_, j) = findmax(Z[i,:])
            j == 1 ? :red  :
            j == 2 ? :blue :
            :green
        end
        plt = Plots.contour(11.0:0.05:15.0, -0.5:0.05:2.5, contour_mixture, contours = false, label = "", size = (600,400))
        Plots.scatter!(plt, X[:,1],X[:,2], label = "", color = colors, markershape = :cross)
        Plots.scatter!(plt, [mean(d)[1] for d in dists], [mean(d)[2] for d in dists], color = :black, label = "Centroids")
        Plots.title!("Estimated mixture, $maxiter iterations")
        Plots.xlabel!("x1")
        Plots.ylabel!("x2")
        Plots.savefig(joinpath(base_img, "wine_mixture_iter$(maxiter)_$iter.pdf"))
    end
end


text:
#!/bin/bash
#=
JULIA="${JULIA:-julia}"
JULIA_CMD="${JULIA_CMD:-$JULIA --color=yes --startup-file=no}"
# below gets the directory name of the script, even if there is a symlink involved
# from https://stackoverflow.com/questions/59895/get-the-source-directory-of-a-bash-script-from-within-the-script-itself
SOURCE="${BASH_SOURCE[0]}"
while [ -h "$SOURCE" ]; do # resolve $SOURCE until the file is no longer a symlink
  DIR="$( cd -P "$( dirname "$SOURCE" )" >/dev/null 2>&1 && pwd )"
  SOURCE="$(readlink "$SOURCE")"
  [[ $SOURCE != /* ]] && SOURCE="$DIR/$SOURCE" # if $SOURCE was a relative symlink, we need to resolve it relative to the path where the symlink file was located
done
DIR="$( cd -P "$( dirname "$SOURCE" )" >/dev/null 2>&1 && pwd )"
export JULIA_PROJECT=$DIR/..
export JULIA_LOAD_PATH=@:@stdlib  # exclude default environment
exec $JULIA_CMD -e 'include(popfirst!(ARGS))' "$SOURCE" "$@"
=#
using ArgParse
using ARMA
using HDF5
using Printf
using Pope.NoiseAnalysis

delete!(ENV, "PYTHONPATH")
using PyCall, PyPlot
pdf = pyimport("matplotlib.backends.backend_pdf")

function parse_commandline()
    s = ArgParseSettings()
    s.description="""Make summary plots of HDF5 noise files (produced by
noise_analysis.jl)."""

    @add_arg_table s begin
        "--outputfile", "-o"
            help = "store the results in OUTPUTFILE (a PDF file) instead of the inferred file or files"
            arg_type = String
        "--replaceoutput", "-r"
            help = "delete and replace any existing output files (default: false)"
            action = :store_true
        "hdf5files"
            help = "1 or more HDF5 files to analyze"
            arg_type = String
            action = :store_arg
            nargs = '+'
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()

    inputs = parsed_args["hdf5files"]
    if parsed_args["outputfile"] == nothing
        outputs = [splitext(input)[1]*".pdf" for input in inputs]
    else
        base,ext = splitext(parsed_args["outputfile"])
        outputs = [base*"_$(i)"*ext for i in 1:length(inputs)]
    end
    for (input,output) in zip(inputs, outputs)
        @show input, output
        pdffile = pdf.PdfPages(output)
        try
            noiseplots(input, pdffile)
        finally
            pdffile.close()
        end
        @show output
    end
end

function noiseplots(input::AbstractString, pdffile)
    NVERT = 4
    MAXCORRSAMPLES = 50
    h5open(input, "r") do hfile
        channels = [parse(Int, k) for k in names(hfile)]
        sort!(channels)

        fig = figure(1, figsize=(11.5, 8))
        allfigs = [fig]

        for (i,chan) in enumerate(channels)
            g = hfile["$(chan)/noise"]
            noise = NoiseAnalysis.hdf5load(g)
            p = noise.model.p
            q = noise.model.q
            Nf = length(noise.psd)
            N = length(noise.autocorr)
            println("Loaded noise data for $(chan)")

            vert = (i-1)%NVERT
            isbottom = i%NVERT == 0
            NHOR = 3
            if N > 2MAXCORRSAMPLES
                NHOR = 5
            end

            ax = subplot(NVERT,NHOR,1+vert*NHOR)
            freq = noise.freqstep*collect(0:Nf-1)
            dfreq = noise.freqstep
            psd_model = ARMA.model_psd(noise.model, Nf) / freq[end]
            loglog(freq[3:end], noise.psd[3:end], "r", freq[2:end], psd_model[2:end], "b")
            title("Chan $(chan) noise PSD")
            isbottom && xlabel("Frequency (Hz)")
            ylabel("PSD (arbs\$^2\$/Hz)")
            grid(true)

            ax = subplot(NVERT,NHOR,2+vert*NHOR)
            ac_model = ARMA.model_covariance(noise.model, N)
            plot(noise.autocorr, "r", ac_model, "b")
            plot([0,N], [0,0], color="gray")
            plot([0],noise.autocorr[1], "ro", [0],ac_model[1], "bo")
            title("Acorr, model($p,$q)")
            isbottom && xlabel("Samples")
            var = noise.autocorr[1]
            for j=1:p
                a = abs(noise.model.expampls[j])
                b = noise.model.expbases[j]
                phi = atan(imag(b), real(b))
                timeconst = -1.0/log(abs(b))
                msg = @sprintf("%.3f*exp(-t/%.3f)", a, timeconst)
                if abs(phi*N) > .01 && imag(b)<0
                    msg = @sprintf("period = %.2f", -2π/phi)
                else
                end
                text(N, var*(1-0.09*j), msg, ha="right", fontsize="small")
            end

            ax = subplot(NVERT,NHOR,3+vert*NHOR)
            plot(noise.autocorr-ac_model, "g")
            plot([0,N], [0,0], color="gray")
            title("Resid (data-model)")
            isbottom && xlabel("Samples")

            if N > 2MAXCORRSAMPLES
                ax = subplot(NVERT,NHOR,4+vert*NHOR)
                ac_model = ARMA.model_covariance(noise.model, MAXCORRSAMPLES)
                plot([0,MAXCORRSAMPLES], [0,0], color="gray")
                plot(noise.autocorr[1:MAXCORRSAMPLES], "r", ac_model, "b")
                plot([0],noise.autocorr[1], "ro", [0],ac_model[1], "bo")
                title("Acorr zoom")
                isbottom && xlabel("Samples")

                ax = subplot(NVERT,NHOR,5+vert*NHOR)
                plot([0,MAXCORRSAMPLES], [0,0], color="gray")
                plot(noise.autocorr[1:MAXCORRSAMPLES]-ac_model, "g")
                title("Resid zoom")
                isbottom && xlabel("Samples")
            end

            if vert == NVERT-1 || chan == channels[end]
                tight_layout()
                pdffile.savefig(fig)
                println("Rendered a page")
                plt.close(fig)
                fig = figure(1, figsize=(11.5, 8))
            end
        end
    end
end

@time main()


text:
# Copyright 2015 Mitchell Kember. Subject to the MIT License.

using Gadfly
using Base.Test

# Plots the distribution of permutations for shuffling n cards with the custom
# shuffling function, using the given number of bins in the histogram. Returns
# the Gadfly plot object.
function distribution(n, samples, bins, shuffle)
  perms = Array(Int, samples)
  for i = 1:samples
    perms[i] = perm_index(shuffle([0:n-1]))
  end
  plot(
    x=perms, Geom.histogram(bincount=bins),
    yintercept=[samples/bins], Geom.hline(color="yellowgreen"),
    Guide.title(""), Guide.xlabel(""), Guide.ylabel(""),
    Scale.y_continuous(minvalue=0)
  )
end

# Saves the distribution as an SVG file with the given name.
function save_svg(dist, name)
  draw(SVG(name, 6inch, 4inch), dist)
end

# Returns a random number between a and b inclusive.
function next_rand(a, b)
  rand(a:b)
end

# Shuffles xs in-place using a bad algorithm.
function simple_shuffle!(xs)
  n = length(xs)
  for _ = 1:n
    a = next_rand(1, n)
    b = next_rand(1, n)
    xs[a], xs[b] = xs[b], xs[a]
  end
  xs
end

# Shuffles xs in-place using a good algorithm.
function knuth_shuffle!(xs)
  n = length(xs)
  for i = 1:n
    j = next_rand(i, n)
    xs[i], xs[j] = xs[j], xs[i]
  end
  xs
end

# Identifies the permutation of [0:n-1] that xs corresponds to, where n is the
# length of xs. Returns zero for ascending order and n!-1 for desceding order.
# Does not check that xs is a permutation.
function perm_index(xs)
  if isempty(xs)
    return 0
  end
  n = length(xs)
  available = ones(xs)
  index = 0
  f = factorial(n-1)
  for i = 1:n-1
    val = xs[i]
    pos = 0
    for j = 1:val
      pos += available[j]
    end
    index += pos * f
    available[val+1] = 0
    f /= n - i
  end
  index
end

# Verifies permutation indices for permutations of [0] through [0:3].
function test_perm_index()
  order = Any[
    [0],
    [0 1; 1 0],
    [0 1 2; 0 2 1; 1 0 2; 1 2 0; 2 0 1; 2 1 0],
    [0 1 2 3; 0 1 3 2; 0 2 1 3; 0 2 3 1; 0 3 1 2; 0 3 2 1;
     1 0 2 3; 1 0 3 2; 1 2 0 3; 1 2 3 0; 1 3 0 2; 1 3 2 0;
     2 0 1 2; 2 0 2 1; 2 1 0 3; 2 1 3 0; 2 3 0 1; 2 3 1 0;
     3 0 1 2; 3 0 2 1; 3 1 0 2; 3 1 2 0; 3 2 0 1; 3 2 1 0]
  ]
  for ps in order
    for i = 1:size(ps, 1)
      @test perm_index(ps[i,:]) == i - 1
    end
  end
end

# Plots the distributions for the simple shuffle and for the Knuth shuffle, and
# saves each to an SVG file.
function main()
  const n = 20
  const samples = 1_000_000
  const bins = 200
  save_svg(distribution(n, samples, bins, simple_shuffle!), "simple.svg")
  save_svg(distribution(n, samples, bins, knuth_shuffle!), "knuth.svg")
end


text:
function f(x::Vector)
  (100.0 - x[1])^2 + x[2]^2
end

function rosenbrock(x::Vector)
  (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

initial_x = [0.0, 0.0]

results = Optim.nelder_mead(f, initial_x)

@assert results.converged
@assert norm(results.minimum - [100.0, 0.0]) < 0.01
@assert length(results.trace.states) == 0

results = Optim.nelder_mead(rosenbrock, initial_x)

@assert results.converged
@assert norm(results.minimum - [1.0, 1.0]) < 0.01
@assert length(results.trace.states) == 0


text:
# Deprecated

function drawLineBetween!(vis,
        cgl::CloudGraph,
        session::AbstractString,
        fr::Symbol,
        to::Symbol;
        scale=0.01,
        name::Symbol=:edges,
        subname::Union{Nothing,Symbol}=nothing,
        color=RGBA(0,1.0,0,0.5)  )
  #

  cv1 = getCloudVert(cgl, session, sym=fr)
  v1 = cloudVertex2ExVertex(cv1)
  cv2 = getCloudVert(cgl, session, sym=to)
  v2 = cloudVertex2ExVertex(cv2)

  drawLineBetween!(vis,session,v1,v2,scale=scale,name=name,subname=subname,color=color   )
  nothing
end

function drawAllBinaryFactorEdges!(vis,
                                   cgl::CloudGraph,
                                   session::AbstractString;
                                   scale=0.01   )
  #

  sloth = findAllBinaryFactors(cgl, session)

  @showprogress 1 "Drawing all binary edges..." for (teeth, toe) in sloth
    color = pointToColor(toe[3])
    drawLineBetween!(vis, cgl, session, toe[1], toe[2], subname=toe[3], scale=scale, color=color)
  end
  nothing
end

function fetchdrawdepthcloudbycvid!(vis,
                                    cloudGraph::CloudGraphs.CloudGraph,
                                    cvid::Int,
                                    vsym::Symbol,
                                    poseswithdepth::Dict,
                                    param::Dict,
                                    sesssym::Symbol;
                                    depthcolormaps=Dict(),
                                    # imshape=(480,640),
                                    wTb::CoordinateTransformations.AbstractAffineMap=
                                            Translation(0,0,0.0) ∘ LinearMap(
                                            CoordinateTransformations.Quat(1.0, 0, 0, 0))   )
                                    # bTc::CoordinateTransformations.AbstractAffineMap=
                                    #       Translation(0,0,0.6) ∘ LinearMap(
                                    #       CoordinateTransformations.Quat(0.5, -0.5, 0.5, -0.5))  )
  #

  if !haskey(poseswithdepth, vsym)
    cache = Dict()

    # fetch copy of big data from CloudGraphs
    cv = CloudGraphs.get_vertex(cloudGraph, cvid, true )

    # depthcolormaps = could be one or more or these options
    # 0, 1 or 2+ color maps per pointcloud
    #  ("none" => "depthframe_image", "none" => "pointcloud")
    # or
    #  ("keyframe_rgb" => "depthframe_image",
    #  "keyframe_segnet" => "depthframe_image")
    # or
    #  ("colors" => "pointcloud")
    for (va, ke) in depthcolormaps
      # prep the detph pointcloud
      cachepointclouds!(cache, cv, ke, param)

      if haskey(cache, ke)
        rgb = Array{Colorant,2}()
        if va[1:4] != "none"
          rgb = retrievePointcloudColorInfo!(cv, va)
        end

        # add color information to the pointcloud
        pointcloud = prepcolordepthcloud!( cvid, cache[ke], rgb=rgb )

        # publish the pointcloud data to Director viewer
        if pointcloud != nothing
          drawpointcloud!(vis, poseswithdepth, vsym, pointcloud, va, param, sesssym, wTb=wTb)
        end
      end
    end
  end
  nothing
end

function fetchdrawposebycvid!(vis,
      cloudGraph::CloudGraphs.CloudGraph,
      cvid::Int,
      poseswithdepth::Dict,
      param::Dict;
      session::AbstractString="",
      depthcolormaps=Dict()  )
      # imshape=(480,640),
      # bTc::CoordinateTransformations.AbstractAffineMap=
      #       Translation(0,0,0.6) ∘ LinearMap(
      #       CoordinateTransformations.Quat(0.5, -0.5, 0.5, -0.5))  )
  #

  # skip big data elements at first
  cv = CloudGraphs.get_vertex(cloudGraph, cvid, false )
  vert = cloudVertex2ExVertex(cv)

  # extract and draw new poses
  wTb = drawpose!(vis, vert, session=session )
  # also draw pose points from variable marginal belief approximation KDE
  drawposepoints!(vis, vert, session=session )

  # also update any parallel transform paths, previous and new
  updateparallelposes!(vis, vert, poseswithdepth, wTb=wTb)

  # check if we can draw depth pointclouds, and add new ones to parallel transform paths
  fetchdrawdepthcloudbycvid!(vis,
        cloudGraph,
        cvid,
        Symbol(vert.label),
        poseswithdepth,
        param,
        Symbol(session),
        depthcolormaps=depthcolormaps,
        # imshape=imshape,
        wTb=wTb  )

  sleep(0.005)
  nothing
end

# dbcoll,
function drawdbsession(vis,
                       cloudGraph::CloudGraphs.CloudGraph,
                       addrdict,
                       param::Dict,
                       poseswithdepth::Dict )
                       # bTc::CoordinateTransformations.AbstractAffineMap=
                       #       Translation(0,0,0.6) ∘ LinearMap(
                       #       CoordinateTransformations.Quat(0.5, -0.5, 0.5, -0.5) )    )
  #

  @show session = addrdict["session"]
  sesssym = Symbol(session)
  DRAWDEPTH = addrdict["draw depth"]=="y" # not going to support just yet
  DRAWEDGES = addrdict["draw edges"]=="y" # not going to support just yet

  # fg = Caesar.initfg(sessionname=addrdict["session"], cloudgraph=cloudGraph)
  println("Fetching pose IDs to be drawn...")
  IDs = getExVertexNeoIDs(cloudGraph.neo4j.connection, label="POSE", session=session, reqbackendset=false);
  landmIDs = getExVertexNeoIDs(cloudGraph.neo4j.connection, label="LANDMARK", session=session, reqbackendset=false);

  @showprogress 1 "Drawing LANDMARK IDs..." for (vid,cvid) in landmIDs
    cv = nothing
    # skip big data elements
    cv = CloudGraphs.get_vertex(cloudGraph, cvid, false )
    vert = cloudVertex2ExVertex(cv)
    x = Symbol(vert.label)

    # vert = getVert(fg, x, api=localapi)
    drawpoint!(vis, vert, session=session)
    # drawposepoints!(vis, vert, session=session )
  end

  depthcolormaps = !DRAWDEPTH  ? Dict() : Dict("keyframe_rgb" => "depthframe_image", "keyframe_segnet" => "depthframe_image", "BSONcolors" => "BSONpointcloud", "none"=>"BSONpointcloud")

  @showprogress 1 "Drawing POSE IDs..." for (vid,cvid) in IDs
    fetchdrawposebycvid!(vis,
          cloudGraph,
          cvid,
          poseswithdepth,
          param,
          session=session,
          depthcolormaps=depthcolormaps  )
  end

  if DRAWEDGES
    println("Going to draw edges...")
    drawAllBinaryFactorEdges!(vis, cloudGraph, session)
  end
  nothing
end


text:
using Documenter, SphericalBesselFunctions

isdefined(Main, :NOPLOTS) && NOPLOTS || include("plots.jl")

makedocs(;
    modules=[SphericalBesselFunctions],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/jagot/SphericalBesselFunctions.jl/blob/{commit}{path}#L{line}",
    sitename="SphericalBesselFunctions.jl",
    authors="Stefanos Carlström <stefanos.carlstrom@gmail.com>"
)

deploydocs(;
    repo="github.com/jagot/SphericalBesselFunctions.jl",
)