text:
function _get_lb(model::Optimizer, index::MOI.VariableIndex)
    return cpo_java_numvar_getlb(model.inner, _info(model, index).variable)
end
function _get_ub(model::Optimizer, index::MOI.VariableIndex)
    return cpo_java_numvar_getub(model.inner, _info(model, index).variable)
end
function _has_lb(model::Optimizer, index::MOI.VariableIndex)
    return _get_lb(model, index) != -IloInfinity &&
           _get_lb(model, index) != Float64(IloMinInt)
end
function _has_ub(model::Optimizer, index::MOI.VariableIndex)
    return _get_ub(model, index) != IloInfinity &&
           _get_ub(model, index) != Float64(IloMaxInt)
end

function _bounds_to_set(model::Optimizer, index::MOI.VariableIndex)
    if _has_lb(model, index)
        if _has_ub(model, index)
            if _info(model, index).ub == _info(model, index).ub
                return MOI.EqualTo{typeof(_info(model, index).ub)}
            else
                return MOI.Interval{typeof(_info(model, index).ub)}
            end
        else
            return MOI.GreaterThan{typeof(_info(model, index).lb)}
        end
    else
        if _has_ub(model, index)
            return MOI.LessThan{typeof(_info(model, index).ub)}
        else
            # No bounds set!
            return nothing
        end
    end
end

function _assert_no_lb(model::Optimizer, i::MOI.VariableIndex, s)
    return _has_lb(model, i) &&
           throw(MOI.LowerBoundAlreadySet{_bounds_to_set(model, i), s}(i))
end
function _assert_no_ub(model::Optimizer, i::MOI.VariableIndex, s)
    return _has_ub(model, i) &&
           throw(MOI.UpperBoundAlreadySet{_bounds_to_set(model, i), s}(i))
end

function _set_lb(model::Optimizer, index::MOI.VariableIndex, lb::Real)
    lb = Float64(lb)
    info = _info(model, index)
    info.lb = lb
    cpo_java_numvar_setlb(model.inner, info.variable, lb)
    return
end
function _set_ub(model::Optimizer, index::MOI.VariableIndex, ub::Real)
    ub = Float64(ub)
    info = _info(model, index)
    info.ub = ub
    cpo_java_numvar_setub(model.inner, info.variable, ub)
    return
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{T}},
) where {T <: Real}
    index = MOI.VariableIndex(c.value)
    return MOI.is_valid(model, index) && _has_ub(model, index)
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{T}},
) where {T <: Real}
    index = MOI.VariableIndex(c.value)
    return MOI.is_valid(model, index) && _has_lb(model, index)
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.Interval{T}},
) where {T <: Real}
    index = MOI.VariableIndex(c.value)
    return MOI.is_valid(model, index) &&
           _has_lb(model, index) &&
           _has_ub(model, index)
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{T}},
) where {T <: Real}
    index = MOI.VariableIndex(c.value)
    return MOI.is_valid(model, index) &&
           _get_lb(model, index) == _get_ub(model, index)
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne},
)
    index = MOI.VariableIndex(c.value)
    if !MOI.is_valid(model, index)
        return false
    end

    info = _info(model, index)
    return info.type == BINARY || info.binary !== nothing
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.Integer},
)
    index = MOI.VariableIndex(c.value)
    if !MOI.is_valid(model, index)
        return false
    end

    info = _info(model, index)
    return info.type == INTEGER || info.integer !== nothing
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintFunction,
    c::MOI.ConstraintIndex{MOI.VariableIndex, <:Any},
)
    println(">> $c")
    MOI.throw_if_not_valid(model, c)
    return MOI.VariableIndex(MOI.VariableIndex(c.value))
end

function MOI.set(
    ::Optimizer,
    ::MOI.ConstraintFunction,
    ::MOI.ConstraintIndex{MOI.VariableIndex, S},
    ::MOI.VariableIndex,
) where {S}
    throw(MOI.SettingVariableIndexFunctionNotAllowed())
    return
end

function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{MOI.VariableIndex, S},
    s::S,
) where {S}
    if S <: MOI.Interval
        _set_ub(model, MOI.VariableIndex(c.value), s.upper, typeof(s.upper))
        _set_lb(model, MOI.VariableIndex(c.value), s.lower, typeof(s.lower))
    elseif S <: MOI.EqualTo
        _set_ub(model, MOI.VariableIndex(c.value), s.value, typeof(s.value))
        _set_lb(model, MOI.VariableIndex(c.value), s.value, typeof(s.value))
    elseif S <: MOI.GreaterThan
        _set_lb(model, MOI.VariableIndex(c.value), s.lower, typeof(s.lower))
    elseif S <: MOI.LessThan
        _set_ub(model, MOI.VariableIndex(c.value), s.upper, typeof(s.upper))
    end
    return
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.VariableIndex,
    s::MOI.LessThan{T},
) where {T <: Real}
    @assert _info(model, f).type == _type_to_variabletype(T)
    _assert_no_ub(model, f, s)
    _set_ub(model, f, s.upper)
    return MOI.ConstraintIndex{MOI.VariableIndex, typeof(s)}(f.value)
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.VariableIndex,
    s::MOI.GreaterThan{T},
) where {T <: Real}
    @assert _info(model, f).type == _type_to_variabletype(T)
    _assert_no_lb(model, f, s)
    _set_lb(model, f, s.lower)
    return MOI.ConstraintIndex{MOI.VariableIndex, typeof(s)}(f.value)
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.VariableIndex,
    s::MOI.EqualTo{T},
) where {T <: Real}
    @assert _info(model, f).type == _type_to_variabletype(T)
    _assert_no_lb(model, f, s)
    _assert_no_ub(model, f, s)
    _set_lb(model, f, s.value)
    _set_ub(model, f, s.value)
    return MOI.ConstraintIndex{MOI.VariableIndex, typeof(s)}(f.value)
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.VariableIndex,
    s::MOI.Interval{T},
) where {T <: Real}
    @assert _info(model, f).type == _type_to_variabletype(T)
    _assert_no_lb(model, f, s)
    _assert_no_ub(model, f, s)
    _set_lb(model, f, s.lower)
    _set_ub(model, f, s.upper)
    return MOI.ConstraintIndex{MOI.VariableIndex, typeof(s)}(f.value)
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.VariableIndex,
    ::MOI.ZeroOne,
)
    info = _info(model, f)
    @assert info.type == CONTINUOUS || info.type == INTEGER
    @assert info.binary === nothing

    # CPLEX Java API expects variables to be created with the right type.
    # Bad trick: create a binary variable, impose equality. Don't show this variable to MOI, though.
    bindex, cindex = MOI.add_constrained_variable(model, MOI.ZeroOne())
    eqcstr =
        cpo_java_eq(model.inner, info.variable, _info(model, bindex).variable)
    cpo_java_add(model.inner, eqcstr)
    info.binary = (bindex, cindex, eqcstr)

    return MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}(
        f.value,
    )
end

function MOI.add_constraint(
    model::Optimizer,
    f::MOI.VariableIndex,
    ::MOI.Integer,
)
    info = _info(model, f)
    @assert info.type == CONTINUOUS
    @assert info.integer === nothing

    # CPLEX Java API expects variables to be created with the right type.
    # Bad trick: create an integer variable, impose equality. Don't show this variable to MOI, though.
    bindex, cindex = MOI.add_constrained_variable(model, MOI.Integer())
    eqcstr =
        cpo_java_eq(model.inner, info.variable, _info(model, bindex).variable)
    cpo_java_add(model.inner, eqcstr)
    info.integer = (bindex, cindex, eqcstr)

    return MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}(
        f.value,
    )
end

# No similar function for a vector of variables, no way to do it more efficiently.

function MOI.delete(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{T}},
) where {T <: Real}
    MOI.throw_if_not_valid(model, c)
    _set_ub(
        model,
        MOI.VariableIndex(c.value),
        (T == Float64) ? IloInfinity : IloMaxInt,
    )
    _info(model, MOI.VariableIndex(c.value)).ub_name = ""
    model.name_to_constraint = nothing
    return
end

function MOI.delete(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{T}},
) where {T <: Real}
    MOI.throw_if_not_valid(model, c)
    _set_lb(
        model,
        MOI.VariableIndex(c.value),
        (T == Float64) ? -IloInfinity : IloMinInt,
    )
    _info(model, MOI.VariableIndex(c.value)).lb_name = ""
    model.name_to_constraint = nothing
    return
end

function MOI.delete(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{T}},
) where {T <: Real}
    MOI.throw_if_not_valid(model, c)
    _set_ub(
        model,
        MOI.VariableIndex(c.value),
        (T == Float64) ? IloInfinity : IloMaxInt,
    )
    _set_lb(
        model,
        MOI.VariableIndex(c.value),
        (T == Float64) ? -IloInfinity : IloMinInt,
    )
    _info(model, MOI.VariableIndex(c.value)).equalto_name = ""
    model.name_to_constraint = nothing
    return
end

function MOI.delete(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.Interval{T}},
) where {T <: Real}
    MOI.throw_if_not_valid(model, c)
    _set_ub(
        model,
        MOI.VariableIndex(c.value),
        (T == Float64) ? IloInfinity : IloMaxInt,
    )
    _set_lb(
        model,
        MOI.VariableIndex(c.value),
        (T == Float64) ? -IloInfinity : IloMinInt,
    )
    _info(model, MOI.VariableIndex(c.value)).interval_name = ""
    model.name_to_constraint = nothing
    return
end

function MOI.delete(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VariableIndex, S},
) where {S <: Union{MOI.ZeroOne, MOI.Integer}}
    MOI.throw_if_not_valid(model, c)
    info = _info(model, MOI.VariableIndex(c.value))
    @assert info.type == CONTINUOUS || info.type == INTEGER

    if S == MOI.ZeroOne
        @assert info.binary !== nothing
        cpo_java_remove(model.inner, info.binary[3])
        info.binary = nothing
    else
        @assert info.integer !== nothing
        cpo_java_remove(model.inner, info.integer[3])
        info.integer = nothing
    end

    return
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{T}},
) where {T <: Real}
    MOI.throw_if_not_valid(model, c)
    return MOI.GreaterThan(_get_lb(model, MOI.VariableIndex(c.value)))
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{T}},
) where {T <: Real}
    MOI.throw_if_not_valid(model, c)
    return MOI.LessThan(_get_ub(model, MOI.VariableIndex(c.value)))
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.EqualTo{T}},
) where {T <: Real}
    MOI.throw_if_not_valid(model, c)
    return MOI.EqualTo(_get_ub(model, MOI.VariableIndex(c.value)))
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.Interval{T}},
) where {T <: Real}
    MOI.throw_if_not_valid(model, c)
    return MOI.Interval(
        _get_lb(model, MOI.VariableIndex(c.value)),
        _get_ub(model, MOI.VariableIndex(c.value)),
    )
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne},
)
    MOI.throw_if_not_valid(model, c)
    return MOI.ZeroOne()
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{MOI.VariableIndex, MOI.Integer},
)
    MOI.throw_if_not_valid(model, c)
    return MOI.Integer()
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{F, S},
) where {F, S}
    return _info(model, c).set
end

# TODO: MOI.set? Remove the constraint and rebuild it?
# TODO: make variable integer/binary?
# TODO: constraint names here? They don't get passed to the solver.

# Implementations of _build_constraint that are mostly useful for reification.
# MOI.add_constraint rather works on the variable bounds than new constraints in this case.
function _build_constraint(
    model::Optimizer,
    f::MOI.VariableIndex,
    s::MOI.GreaterThan{T},
) where {T <: Real}
    return cpo_java_ge(model.inner, _parse(model, f), s.lower)
end

function _build_constraint(
    model::Optimizer,
    f::MOI.VariableIndex,
    s::MOI.LessThan{T},
) where {T <: Real}
    return cpo_java_le(model.inner, _parse(model, f), s.upper)
end

function _build_constraint(
    model::Optimizer,
    f::MOI.VariableIndex,
    s::MOI.EqualTo{T},
) where {T <: Real}
    return cpo_java_eq(model.inner, _parse(model, f), s.value)
end

function _build_constraint(
    model::Optimizer,
    f::MOI.VariableIndex,
    s::MOI.Interval{T},
) where {T <: Real}
    if s.lower == Inf
        return _build_constraint(model, f, MOI.LessThan(s.upper))
    elseif s.upper == Inf
        return _build_constraint(model, f, MOI.GreaterThan(s.lower))
    elseif s.lower == s.upper
        return _build_constraint(model, f, MOI.EqualTo(s.lower))
    end

    return cpo_java_range(model.inner, s.lower, _parse(model, f), s.upper)
end


text:
function fwd_grad_like(gp)
    dσ = ForwardDiff.gradient(x -> negloglikelihood(x, gp.dpars.G, gp.dpars.u, gp), gp.dpars.σ)
    dG = ForwardDiff.gradient(x -> negloglikelihood(gp.dpars.σ, x, gp.dpars.u, gp), gp.dpars.G)
    du = ForwardDiff.gradient(x -> negloglikelihood(gp.dpars.σ, gp.dpars.G, x, gp), gp.dpars.u)
    (dσ, dG, du)
end

function fit!(gp, its; fwd=true, ls_lr=2e-3, σ_lr=2e-3, show_like=false)

    opt_ls = Flux.NADAM(ls_lr)
    opt_σ = Flux.NADAM(σ_lr)

    for i in 1:its
    # Σ = Diagonal(vcat([gp.dpars.σ[i]^2 * ones(size(gp.data.X)[1]) for i in 1:gp.D]...))
    # B = VolterraGP.fill_K(gp.data.X, gp.data.X, gp) + Σ + 1e-5 * I
    # println(sort(eigvals(B)))
    # println(maximum(B' - /B))
        if fwd
            grads = fwd_grad_like(gp)

            Flux.Optimise.update!(opt_σ, gp.dpars.σ, grads[1])
            Flux.Optimise.update!(opt_ls, gp.dpars.G, grads[2])
            Flux.Optimise.update!(opt_ls, gp.dpars.u, grads[3])

        else #use reverse diff
            grads = Flux.gradient(Flux.params(gp.dpars.σ, gp.dpars.G, gp.dpars.u)) do
                negloglikelihood(gp)
            end

           for p in (gp.dpars.G, gp.dpars.u)
               Flux.Optimise.update!(opt_ls, p, grads[p])
           end

            Flux.Optimise.update!(opt_σ, gp.dpars.σ, grads[gp.dpars.σ])
        end

        println("it: ", i)

        # println("params: ", gp.dpars.G, gp.dpars.σ, gp.dpars.u)
        if show_like
            println(" negloglike:", negloglikelihood(gp))
        end

        if its ÷ 5 == 0
            plotgp(gp.data.X, gp)
        end
    end
end

ext: jl


text:
include("fastpolys.jl")

using Main.FastPolys
using StaticArrays
using BenchmarkTools
using StaticPolynomials
using CommonSubexpressions

module NB5I

using Main.FastPolys
using StaticArrays
using StaticPolynomials
using CommonSubexpressions

const P1_1 = (1,2,3,4,5,6,7,8,9,10,)

const P2_1 = (1,1,1,2,2,3,1,1,1,2,3,4,2,2,3,4,3,4,5,5,6,5,5,6,7,6,7,8,8,9,)
const P2_2 = (2,3,4,3,4,4,5,6,7,5,6,7,8,9,8,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,)

const P3_1 = (1,2,3,4,5,6,7,8,9,10,)

const P4_1 = (1,1,1,2,1,1,1,2,2,3,4,3,4,2,3,4,5,5,6,7,)
const P4_2 = (2,2,3,3,5,5,6,5,5,6,7,6,7,8,8,9,6,8,8,9,)
const P4_3 = (3,4,4,4,6,7,7,8,9,8,9,10,10,9,10,10,7,9,10,10,)

const P5_1 = (1,2,3,4,5,6,7,8,9,10,)

const P6_1 = (1,1,2,3,4,)
const P6_2 = (2,5,5,6,7,)
const P6_3 = (3,6,8,8,9,)
const P6_4 = (4,7,9,10,10,)

const P7_1 = (1,2,3,4,5,6,7,8,9,10,)

const P8_1 = (1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,3,3,1,1,1,2,2,3,4,3,4,2,3,4,)
const P8_2 = (2,2,2,2,2,2,2,3,4,2,3,4,3,4,3,4,4,4,5,5,5,5,5,5,5,6,6,5,6,7,)
const P8_3 = (3,3,3,3,3,3,5,5,5,5,6,7,5,5,6,7,6,7,6,6,6,6,7,6,7,7,7,8,8,8,)
const P8_4 = (4,4,4,4,4,4,6,6,6,8,8,9,8,8,8,9,8,9,7,7,7,8,8,8,9,8,9,9,9,9,)
const P8_5 = (5,6,7,8,9,10,7,7,7,9,10,10,9,9,10,10,10,10,8,9,10,9,9,10,10,10,10,10,10,10,)

const P9_1 = (1,2,3,4,5,6,7,8,9,10,)

const P10_1 = (1,2,3,4,5,6,7,8,9,10,)

# Irreducible secondaries for NBody=5and deg=12
 # : definitions at the beginning of the file
const IS1_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS1_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS2_1 = (1,1,1,2,2,3,5,5,6,8,)
const IS2_2 = (2,3,4,3,4,4,6,7,7,9,)
const IS2_3 = (5,6,7,8,9,10,8,9,10,10,)

const IS3_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS3_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS4_1 = (1,1,1,2,2,3,1,1,1,2,3,4,2,2,3,4,3,4,5,5,6,5,5,6,7,6,7,8,8,9,)
const IS4_2 = (2,3,4,3,4,4,5,6,7,5,6,7,8,9,8,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,)

const IS5_1 = (1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,8,9,10,5,5,6,7,6,7,8,9,10,8,9,10,)
const IS5_2 = (2,3,4,1,1,1,3,4,2,2,4,3,1,1,1,2,2,3,6,7,5,5,7,6,5,5,6,9,8,8,)
const IS5_3 = (5,6,7,5,6,7,8,9,8,9,10,10,2,3,4,3,4,4,8,9,8,9,10,10,6,7,7,10,10,9,)

const IS6_1 = (1,1,1,2,2,3,4,3,4,2,3,4,1,1,1,2,2,3,4,3,4,2,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,6,7,5,6,7,8,9,8,9,10,10,)
const IS6_2 = (2,2,3,1,1,1,1,1,1,3,2,2,5,5,6,5,5,6,7,6,7,8,8,9,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,5,5,8,8,9,5,5,6,7,6,7,)
const IS6_3 = (3,4,4,3,4,2,2,4,3,4,4,3,6,7,7,8,9,8,9,10,10,9,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,7,7,6,9,10,10,9,8,10,10,8,9,)

const IS7_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,8,9,8,9,10,10,)
const IS7_2 = (3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,4,3,4,3,2,2,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,7,6,7,6,5,5,6,7,5,5,7,6,)
const IS7_3 = (5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,8,9,8,9,10,10,9,8,10,10,8,9,)

const IS8_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS8_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS9_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS9_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS10_1 = (1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,8,9,10,5,5,6,7,6,7,8,9,10,8,9,10,)
const IS10_2 = (2,3,4,1,1,1,3,4,2,2,4,3,1,1,1,2,2,3,6,7,5,5,7,6,5,5,6,9,8,8,)
const IS10_3 = (5,6,7,5,6,7,8,9,8,9,10,10,2,3,4,3,4,4,8,9,8,9,10,10,6,7,7,10,10,9,)

const IS11_1 = (1,1,1,2,2,3,1,1,1,2,3,4,2,2,3,4,3,4,5,5,6,5,5,6,7,6,7,8,8,9,)
const IS11_2 = (2,3,4,3,4,4,5,6,7,5,6,7,8,9,8,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,)
const IS11_3 = (5,6,7,8,9,10,2,3,4,1,1,1,3,4,2,2,4,3,8,9,10,6,7,5,5,7,6,10,9,8,)

const IS12_1 = (1,1,1,2,2,3,4,3,4,2,3,4,1,1,1,2,2,3,4,3,4,2,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,6,7,5,6,7,8,9,8,9,10,10,)
const IS12_2 = (2,2,3,1,1,1,1,1,1,3,2,2,5,5,6,5,5,6,7,6,7,8,8,9,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,5,5,8,8,9,5,5,6,7,6,7,)
const IS12_3 = (3,4,4,3,4,2,2,4,3,4,4,3,6,7,7,8,9,8,9,10,10,9,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,7,7,6,9,10,10,9,8,10,10,8,9,)

const IS13_1 = (1,1,1,1,1,1,2,2,3,2,2,3,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,5,5,6,7,6,7,8,8,9,5,5,6,5,5,6,7,6,7,8,8,9,)
const IS13_2 = (2,2,3,4,3,4,3,4,4,3,4,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,)
const IS13_3 = (3,4,2,2,4,3,1,1,1,4,3,2,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,1,1,1,2,2,3,4,3,4,2,3,4,7,6,5,9,8,10,10,8,9,5,6,7,)

const IS14_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,8,9,8,9,10,10,)
const IS14_2 = (3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,4,3,4,3,2,2,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,7,6,7,6,5,5,6,7,5,5,7,6,)
const IS14_3 = (5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,8,9,8,9,10,10,9,8,10,10,8,9,)

const IS15_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS15_2 = (2,2,2,2,3,3,1,1,1,1,1,1,1,1,1,1,1,1,3,3,2,2,2,2,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,3,3,2,2,2,2,3,3,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,6,5,5,5,5,6,7,5,5,7,6,5,5,5,5,6,6,5,5,6,7,6,7,)
const IS15_3 = (3,4,3,4,4,4,3,4,2,2,4,3,3,4,2,2,4,3,4,4,4,3,4,3,5,5,5,5,6,6,5,5,6,7,6,7,5,5,6,7,6,7,8,8,8,9,8,9,2,2,3,4,3,4,2,2,3,4,3,4,3,4,3,4,4,4,3,4,3,4,4,4,6,7,5,5,7,6,6,7,5,5,7,6,5,5,5,5,6,6,9,8,9,8,8,8,7,7,7,6,7,6,8,8,8,9,8,9,6,7,6,7,7,7,9,8,9,8,8,8,)
const IS15_4 = (5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,6,7,7,7,8,9,8,9,10,10,8,9,8,9,10,10,9,9,10,10,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,8,9,8,9,10,10,8,9,8,9,10,10,6,7,6,7,7,7,10,10,10,10,9,9,8,9,8,9,10,10,9,9,10,10,10,10,9,8,10,10,8,9,10,10,10,10,9,9,)

const IS16_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS16_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS17_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS17_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS18_1 = (1,1,1,2,2,3,1,1,1,2,3,4,2,2,3,4,3,4,5,5,6,5,5,6,7,6,7,8,8,9,)
const IS18_2 = (2,3,4,3,4,4,5,6,7,5,6,7,8,9,8,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,)

const IS19_1 = (1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,8,9,10,5,5,6,7,6,7,8,9,10,8,9,10,)
const IS19_2 = (2,3,4,1,1,1,3,4,2,2,4,3,1,1,1,2,2,3,6,7,5,5,7,6,5,5,6,9,8,8,)
const IS19_3 = (5,6,7,5,6,7,8,9,8,9,10,10,2,3,4,3,4,4,8,9,8,9,10,10,6,7,7,10,10,9,)

const IS20_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS20_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)
const IS20_3 = (5,6,7,5,6,7,8,9,8,9,10,10,2,3,4,1,1,1,3,4,2,2,4,3,2,3,4,1,1,1,3,4,2,2,4,3,8,9,8,9,10,10,6,7,5,5,7,6,6,7,5,5,7,6,10,10,9,8,9,8,)

const IS21_1 = (1,1,1,2,2,3,5,5,6,8,)
const IS21_2 = (2,3,4,3,4,4,6,7,7,9,)
const IS21_3 = (5,6,7,8,9,10,8,9,10,10,)

const IS22_1 = (1,1,1,2,2,3,4,3,4,2,3,4,1,1,1,2,2,3,4,3,4,2,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,6,7,5,6,7,8,9,8,9,10,10,)
const IS22_2 = (2,2,3,1,1,1,1,1,1,3,2,2,5,5,6,5,5,6,7,6,7,8,8,9,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,5,5,8,8,9,5,5,6,7,6,7,)
const IS22_3 = (3,4,4,3,4,2,2,4,3,4,4,3,6,7,7,8,9,8,9,10,10,9,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,7,7,6,9,10,10,9,8,10,10,8,9,)

const IS23_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS23_2 = (2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)
const IS23_3 = (3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,4,3,4,3,2,2,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,7,6,7,6,5,5,9,8,10,10,8,9,9,8,10,10,8,9,5,5,6,7,6,7,)

const IS24_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,8,9,8,9,10,10,)
const IS24_2 = (3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,4,3,4,3,2,2,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,7,6,7,6,5,5,6,7,5,5,7,6,)
const IS24_3 = (5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,8,9,8,9,10,10,9,8,10,10,8,9,)

const IS25_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS25_2 = (2,2,2,2,3,3,1,1,1,1,1,1,1,1,1,1,1,1,3,3,2,2,2,2,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,3,3,2,2,2,2,3,3,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,6,5,5,5,5,6,7,5,5,7,6,5,5,5,5,6,6,5,5,6,7,6,7,)
const IS25_3 = (3,4,3,4,4,4,3,4,2,2,4,3,3,4,2,2,4,3,4,4,4,3,4,3,5,5,5,5,6,6,5,5,6,7,6,7,5,5,6,7,6,7,8,8,8,9,8,9,2,2,3,4,3,4,2,2,3,4,3,4,3,4,3,4,4,4,3,4,3,4,4,4,6,7,5,5,7,6,6,7,5,5,7,6,5,5,5,5,6,6,9,8,9,8,8,8,7,7,7,6,7,6,8,8,8,9,8,9,6,7,6,7,7,7,9,8,9,8,8,8,)
const IS25_4 = (5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,6,7,7,7,8,9,8,9,10,10,8,9,8,9,10,10,9,9,10,10,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,8,9,8,9,10,10,8,9,8,9,10,10,6,7,6,7,7,7,10,10,10,10,9,9,8,9,8,9,10,10,9,9,10,10,10,10,9,8,10,10,8,9,10,10,10,10,9,9,)

const IS26_1 = (1,1,1,1,1,1,2,2,3,2,2,3,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,5,5,6,7,6,7,8,8,9,5,5,6,5,5,6,7,6,7,8,8,9,)
const IS26_2 = (2,2,3,4,3,4,3,4,4,3,4,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,)
const IS26_3 = (3,4,2,2,4,3,1,1,1,4,3,2,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,2,2,3,4,3,4,2,3,4,7,6,5,6,7,5,5,7,6,5,6,7,)
const IS26_4 = (5,5,6,7,6,7,8,9,10,8,9,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,8,9,10,6,7,5,5,7,6,10,9,8,8,9,10,9,8,10,10,8,9,10,9,8,)

const IS27_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS27_2 = (5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,8,9,8,9,10,10,6,7,5,5,7,6,6,7,5,5,7,6,10,10,9,8,9,8,8,9,8,9,10,10,6,7,5,5,7,6,6,7,5,5,7,6,10,10,9,8,9,8,)
const IS27_3 = (3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,4,3,4,3,2,2,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,7,6,7,6,5,5,9,8,10,10,8,9,9,8,10,10,8,9,5,5,6,7,6,7,)

const IS28_1 = (1,2,3,4,5,6,7,8,9,10,)

const IS29_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS29_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS30_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS30_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS31_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS31_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS32_1 = (1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,8,9,10,5,5,6,7,6,7,8,9,10,8,9,10,)
const IS32_2 = (2,3,4,1,1,1,3,4,2,2,4,3,1,1,1,2,2,3,6,7,5,5,7,6,5,5,6,9,8,8,)
const IS32_3 = (5,6,7,5,6,7,8,9,8,9,10,10,2,3,4,3,4,4,8,9,8,9,10,10,6,7,7,10,10,9,)

const IS33_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS33_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)
const IS33_3 = (5,6,7,5,6,7,8,9,8,9,10,10,2,3,4,1,1,1,3,4,2,2,4,3,2,3,4,1,1,1,3,4,2,2,4,3,8,9,8,9,10,10,6,7,5,5,7,6,6,7,5,5,7,6,10,10,9,8,9,8,)

const IS34_1 = (1,1,1,2,2,3,1,1,1,2,3,4,2,2,3,4,3,4,5,5,6,5,5,6,7,6,7,8,8,9,)
const IS34_2 = (2,3,4,3,4,4,5,6,7,5,6,7,8,9,8,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,)
const IS34_3 = (5,6,7,8,9,10,2,3,4,1,1,1,3,4,2,2,4,3,8,9,10,6,7,5,5,7,6,10,9,8,)

const IS35_1 = (1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,8,9,10,5,5,6,7,6,7,8,9,10,8,9,10,)
const IS35_2 = (2,3,4,1,1,1,3,4,2,2,4,3,1,1,1,2,2,3,6,7,5,5,7,6,5,5,6,9,8,8,)
const IS35_3 = (5,6,7,5,6,7,8,9,8,9,10,10,2,3,4,3,4,4,8,9,8,9,10,10,6,7,7,10,10,9,)

const IS36_1 = (1,1,1,2,2,3,4,3,4,2,3,4,1,1,1,2,2,3,4,3,4,2,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,6,7,5,6,7,8,9,8,9,10,10,)
const IS36_2 = (2,2,3,1,1,1,1,1,1,3,2,2,5,5,6,5,5,6,7,6,7,8,8,9,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,5,5,8,8,9,5,5,6,7,6,7,)
const IS36_3 = (3,4,4,3,4,2,2,4,3,4,4,3,6,7,7,8,9,8,9,10,10,9,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,7,7,6,9,10,10,9,8,10,10,8,9,)

const IS37_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS37_2 = (2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)
const IS37_3 = (3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,4,3,4,3,2,2,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,7,6,7,6,5,5,9,8,10,10,8,9,9,8,10,10,8,9,5,5,6,7,6,7,)

const IS38_1 = (1,1,1,1,1,1,2,2,3,2,2,3,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,5,5,6,7,6,7,8,8,9,5,5,6,5,5,6,7,6,7,8,8,9,)
const IS38_2 = (2,2,3,4,3,4,3,4,4,3,4,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,6,7,7,8,9,8,9,10,10,9,10,10,)
const IS38_3 = (3,4,2,2,4,3,1,1,1,4,3,2,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,1,1,1,2,2,3,4,3,4,2,3,4,7,6,5,9,8,10,10,8,9,5,6,7,)

const IS39_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,8,9,8,9,10,10,)
const IS39_2 = (3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,4,3,4,3,2,2,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,7,6,7,6,5,5,6,7,5,5,7,6,)
const IS39_3 = (5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,8,9,8,9,10,10,9,8,10,10,8,9,)

const IS40_1 = (1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS40_2 = (2,2,2,2,3,3,1,1,1,1,1,1,1,1,1,1,1,1,3,3,2,2,2,2,2,2,3,4,3,4,1,1,1,1,1,1,3,4,2,2,4,3,3,4,2,2,4,3,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,3,3,2,2,2,2,3,3,1,1,1,1,1,1,2,2,3,4,3,4,2,2,3,4,3,4,2,2,3,4,3,4,6,6,5,5,5,5,6,7,5,5,7,6,5,5,5,5,6,6,5,5,6,7,6,7,)
const IS40_3 = (3,4,3,4,4,4,3,4,2,2,4,3,3,4,2,2,4,3,4,4,4,3,4,3,5,5,5,5,6,6,5,5,6,7,6,7,5,5,6,7,6,7,8,8,8,9,8,9,2,2,3,4,3,4,2,2,3,4,3,4,3,4,3,4,4,4,3,4,3,4,4,4,6,7,5,5,7,6,6,7,5,5,7,6,5,5,5,5,6,6,9,8,9,8,8,8,7,7,7,6,7,6,8,8,8,9,8,9,6,7,6,7,7,7,9,8,9,8,8,8,)
const IS40_4 = (5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,6,7,6,7,7,7,8,9,8,9,10,10,8,9,8,9,10,10,9,9,10,10,10,10,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,8,9,8,9,10,10,8,9,8,9,10,10,6,7,6,7,7,7,10,10,10,10,9,9,8,9,8,9,10,10,9,9,10,10,10,10,9,8,10,10,8,9,10,10,10,10,9,9,)

const IS41_1 = (1,2,3,4,5,6,7,8,9,10,)

const IS42_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS42_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS43_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS43_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS44_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS44_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

const IS45_1 = (1,2,3,4,5,6,7,8,9,10,)

const IS46_1 = (1,1,1,2,3,4,2,2,3,4,3,4,1,1,1,2,3,4,2,2,3,4,3,4,5,6,7,5,6,7,8,9,8,9,10,10,5,5,6,7,6,7,5,5,6,7,6,7,8,9,8,9,10,10,8,9,8,9,10,10,)
const IS46_2 = (2,3,4,1,1,1,3,4,2,2,4,3,5,6,7,5,6,7,8,9,8,9,10,10,1,1,1,2,3,4,2,2,3,4,3,4,6,7,5,5,7,6,8,9,8,9,10,10,5,5,6,7,6,7,9,8,10,10,8,9,)

# Primary invariants for NBody=5and deg=12
 # : definitions of the types at the beginning of the file
const P1 = Val((P1_1,))
const P2 = Val((P2_1,P2_2,))
const P3 = Val((P3_1,))
const P4 = Val((P4_1,P4_2,P4_3,))
const P5 = Val((P5_1,))
const P6 = Val((P6_1,P6_2,P6_3,P6_4,))
const P7 = Val((P7_1,))
const P8 = Val((P8_1,P8_2,P8_3,P8_4,P8_5,))
const P9 = Val((P9_1,))
const P10 = Val((P10_1,))
# Irreducible secondaries for NBody=5and deg=12
 # : definitions of the types at the beginning of the file
const IS1 = Val((IS1_1,IS1_2,))
const IS2 = Val((IS2_1,IS2_2,IS2_3,))
const IS3 = Val((IS3_1,IS3_2,))
const IS4 = Val((IS4_1,IS4_2,))
const IS5 = Val((IS5_1,IS5_2,IS5_3,))
const IS6 = Val((IS6_1,IS6_2,IS6_3,))
const IS7 = Val((IS7_1,IS7_2,IS7_3,))
const IS8 = Val((IS8_1,IS8_2,))
const IS9 = Val((IS9_1,IS9_2,))
const IS10 = Val((IS10_1,IS10_2,IS10_3,))
const IS11 = Val((IS11_1,IS11_2,IS11_3,))
const IS12 = Val((IS12_1,IS12_2,IS12_3,))
const IS13 = Val((IS13_1,IS13_2,IS13_3,))
const IS14 = Val((IS14_1,IS14_2,IS14_3,))
const IS15 = Val((IS15_1,IS15_2,IS15_3,IS15_4,))
const IS16 = Val((IS16_1,IS16_2,))
const IS17 = Val((IS17_1,IS17_2,))
const IS18 = Val((IS18_1,IS18_2,))
const IS19 = Val((IS19_1,IS19_2,IS19_3,))
const IS20 = Val((IS20_1,IS20_2,IS20_3,))
const IS21 = Val((IS21_1,IS21_2,IS21_3,))
const IS22 = Val((IS22_1,IS22_2,IS22_3,))
const IS23 = Val((IS23_1,IS23_2,IS23_3,))
const IS24 = Val((IS24_1,IS24_2,IS24_3,))
const IS25 = Val((IS25_1,IS25_2,IS25_3,IS25_4,))
const IS26 = Val((IS26_1,IS26_2,IS26_3,IS26_4,))
const IS27 = Val((IS27_1,IS27_2,IS27_3,))
const IS28 = Val((IS28_1,))
const IS29 = Val((IS29_1,IS29_2,))
const IS30 = Val((IS30_1,IS30_2,))
const IS31 = Val((IS31_1,IS31_2,))
const IS32 = Val((IS32_1,IS32_2,IS32_3,))
const IS33 = Val((IS33_1,IS33_2,IS33_3,))
const IS34 = Val((IS34_1,IS34_2,IS34_3,))
const IS35 = Val((IS35_1,IS35_2,IS35_3,))
const IS36 = Val((IS36_1,IS36_2,IS36_3,))
const IS37 = Val((IS37_1,IS37_2,IS37_3,))
const IS38 = Val((IS38_1,IS38_2,IS38_3,))
const IS39 = Val((IS39_1,IS39_2,IS39_3,))
const IS40 = Val((IS40_1,IS40_2,IS40_3,IS40_4,))
const IS41 = Val((IS41_1,))
const IS42 = Val((IS42_1,IS42_2,))
const IS43 = Val((IS43_1,IS43_2,))
const IS44 = Val((IS44_1,IS44_2,))
const IS45 = Val((IS45_1,))
const IS46 = Val((IS46_1,IS46_2,))

function invold(x1::SVector{10, T}) where {T}
   x2 = x1.*x1
   x3 = x2.*x1
   x4 = x3.*x1
   x5 = x4.*x1
   x6 = x5.*x1
   x7 = x6.*x1
   x8 = x7.*x1
   x9 = x8.*x1
   #------------------------------------------------
   # Primaries
   #------------------------------------------------

   # Primary invariants for NBody=5and deg=12
    # : what goes in the function for the evaluation
   P1 = fpoly((x1,) , NB5I.P1)
   P2 = fpoly((x1,x1,) , NB5I.P2)
   P3 = fpoly((x2,) , NB5I.P3)
   P4 = fpoly((x1,x1,x1,) , NB5I.P4)
   P5 = fpoly((x3,) , NB5I.P5)
   P6 = fpoly((x1,x1,x1,x1,) , NB5I.P6)
   P7 = fpoly((x4,) , NB5I.P7)
   P8 = fpoly((x1,x1,x1,x1,x1,) , NB5I.P8)
   P9 = fpoly((x5,) , NB5I.P9)
   P10 = fpoly((x6,) , NB5I.P10)

      #------------------------------------------------
      # Irreducible secondaries
      #------------------------------------------------

   # Irreducible secondaries for NBody=5and deg=12
    # : what goes in the function for the evaluation
   IS1 = fpoly((x2,x1,) , NB5I.IS1)
   IS2 = fpoly((x1,x1,x1,) , NB5I.IS2)
   IS3 = fpoly((x3,x1,) , NB5I.IS3)
   IS4 = fpoly((x2,x2,) , NB5I.IS4)
   IS5 = fpoly((x2,x1,x1,) , NB5I.IS5)
   IS6 = fpoly((x2,x1,x1,) , NB5I.IS6)
   IS7 = fpoly((x2,x1,x1,) , NB5I.IS7)
   IS8 = fpoly((x4,x1,) , NB5I.IS8)
   IS9 = fpoly((x3,x2,) , NB5I.IS9)
   IS10 = fpoly((x3,x1,x1,) , NB5I.IS10)
   IS11 = fpoly((x2,x2,x1,) , NB5I.IS11)
   IS12 = fpoly((x3,x1,x1,) , NB5I.IS12)
   IS13 = fpoly((x2,x2,x1,) , NB5I.IS13)
   IS14 = fpoly((x3,x1,x1,) , NB5I.IS14)
   IS15 = fpoly((x2,x1,x1,x1,) , NB5I.IS15)
   IS16 = fpoly((x5,x1,) , NB5I.IS16)
   IS17 = fpoly((x4,x2,) , NB5I.IS17)
   IS18 = fpoly((x3,x3,) , NB5I.IS18)
   IS19 = fpoly((x4,x1,x1,) , NB5I.IS19)
   IS20 = fpoly((x3,x2,x1,) , NB5I.IS20)
   IS21 = fpoly((x2,x2,x2,) , NB5I.IS21)
   IS22 = fpoly((x4,x1,x1,) , NB5I.IS22)
   IS23 = fpoly((x3,x2,x1,) , NB5I.IS23)
   IS24 = fpoly((x4,x1,x1,) , NB5I.IS24)
   IS25 = fpoly((x3,x1,x1,x1,) , NB5I.IS25)
   IS26 = fpoly((x2,x2,x1,x1,) , NB5I.IS26)
   IS27 = fpoly((x3,x2,x1,) , NB5I.IS27)
   IS28 = fpoly((x7,) , NB5I.IS28)
   IS29 = fpoly((x6,x1,) , NB5I.IS29)
   IS30 = fpoly((x5,x2,) , NB5I.IS30)
   IS31 = fpoly((x4,x3,) , NB5I.IS31)
   IS32 = fpoly((x5,x1,x1,) , NB5I.IS32)
   IS33 = fpoly((x4,x2,x1,) , NB5I.IS33)
   IS34 = fpoly((x3,x3,x1,) , NB5I.IS34)
   IS35 = fpoly((x3,x2,x2,) , NB5I.IS35)
   IS36 = fpoly((x5,x1,x1,) , NB5I.IS36)
   IS37 = fpoly((x4,x2,x1,) , NB5I.IS37)
   IS38 = fpoly((x3,x3,x1,) , NB5I.IS38)
   IS39 = fpoly((x5,x1,x1,) , NB5I.IS39)
   IS40 = fpoly((x4,x1,x1,x1,) , NB5I.IS40)
   IS41 = fpoly((x8,) , NB5I.IS41)
   IS42 = fpoly((x7,x1,) , NB5I.IS42)
   IS43 = fpoly((x6,x2,) , NB5I.IS43)
   IS44 = fpoly((x5,x3,) , NB5I.IS44)
   IS45 = fpoly((x9,) , NB5I.IS45)
   IS46 = fpoly((x8,x1,) , NB5I.IS46)

   return (@SVector [P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,
         IS1, IS2, IS3, IS4, IS5, IS6, IS7, IS8, IS9, IS10,
         IS11, IS12, IS13, IS14, IS15, IS16, IS17, IS18, IS19,
         IS20, IS21, IS22, IS23, IS24, IS25, IS26, IS27, IS28, IS29, IS30, IS31,
         IS32, IS33, IS34, IS35, IS36, IS37, IS38, IS39, IS40, IS41, IS42, IS43,
         IS44, IS45, IS46])
end

const PP1 = (P1_1,)
const PP2 = (P2_1,P2_2,)
const PP3 = (P3_1,)
const PP4 = (P4_1,P4_2,P4_3,)
const PP5 = (P5_1,)
const PP6 = (P6_1,P6_2,P6_3,P6_4,)
const PP7 = (P7_1,)
const PP8 = (P8_1,P8_2,P8_3,P8_4,P8_5,)
const PP9 = (P9_1,)
const PP10 = (P10_1,)

const ISS1 = (IS1_1,IS1_2,)
const ISS2 = (IS2_1,IS2_2,IS2_3,)
const ISS3 = (IS3_1,IS3_2,)
const ISS4 = (IS4_1,IS4_2,)
const ISS5 = (IS5_1,IS5_2,IS5_3,)
const ISS6 = (IS6_1,IS6_2,IS6_3,)
const ISS7 = (IS7_1,IS7_2,IS7_3,)
const ISS8 = (IS8_1,IS8_2,)
const ISS9 = (IS9_1,IS9_2,)
const ISS10 = (IS10_1,IS10_2,IS10_3,)
const ISS11 = (IS11_1,IS11_2,IS11_3,)
const ISS12 = (IS12_1,IS12_2,IS12_3,)
const ISS13 = (IS13_1,IS13_2,IS13_3,)
const ISS14 = (IS14_1,IS14_2,IS14_3,)
const ISS15 = (IS15_1,IS15_2,IS15_3,IS15_4,)
const ISS16 = (IS16_1,IS16_2,)
const ISS17 = (IS17_1,IS17_2,)
const ISS18 = (IS18_1,IS18_2,)
const ISS19 = (IS19_1,IS19_2,IS19_3,)
const ISS20 = (IS20_1,IS20_2,IS20_3,)
const ISS21 = (IS21_1,IS21_2,IS21_3,)
const ISS22 = (IS22_1,IS22_2,IS22_3,)
const ISS23 = (IS23_1,IS23_2,IS23_3,)
const ISS24 = (IS24_1,IS24_2,IS24_3,)
const ISS25 = (IS25_1,IS25_2,IS25_3,IS25_4,)
const ISS26 = (IS26_1,IS26_2,IS26_3,IS26_4,)
const ISS27 = (IS27_1,IS27_2,IS27_3,)
const ISS28 = (IS28_1,)
const ISS29 = (IS29_1,IS29_2,)
const ISS30 = (IS30_1,IS30_2,)
const ISS31 = (IS31_1,IS31_2,)
const ISS32 = (IS32_1,IS32_2,IS32_3,)
const ISS33 = (IS33_1,IS33_2,IS33_3,)
const ISS34 = (IS34_1,IS34_2,IS34_3,)
const ISS35 = (IS35_1,IS35_2,IS35_3,)
const ISS36 = (IS36_1,IS36_2,IS36_3,)
const ISS37 = (IS37_1,IS37_2,IS37_3,)
const ISS38 = (IS38_1,IS38_2,IS38_3,)
const ISS39 = (IS39_1,IS39_2,IS39_3,)
const ISS40 = (IS40_1,IS40_2,IS40_3,IS40_4,)
const ISS41 = (IS41_1,)
const ISS42 = (IS42_1,IS42_2,)
const ISS43 = (IS43_1,IS43_2,)
const ISS44 = (IS44_1,IS44_2,)
const ISS45 = (IS45_1,)
const ISS46 = (IS46_1,IS46_2,)

function _gen_sym(pows, vars)
   I = sortperm(vars)
   return Symbol(prod( "X$(v)_$(p)_x_"
                       for (v,p) in zip(vars[I], pows[I]) )[1:end-2])
end

function push_monomial!(exprs, syms, pows, vars)
   @assert length(vars) == length(pows)
   s = _gen_sym(pows, vars)
   if s in syms
      # if this monomial already exists, there is nothing to do
      # just return the generated symbol
      return s
   end
   # if the symbol does not yet exist we build it recursively
   # TODO: here, the recursion could be taken further and a proper search
   # over similar symbols could be performed, thigh might
   # speed up things further
   s1 = push_monomial!(exprs, syms, pows[1:end-1], vars[1:end-1])
   s2 = push_monomial!(exprs, syms, pows[end:end], vars[end:end])
   push!(exprs, :($(s) = $(s1)*$(s2)))
   push!(syms, s)
   return s
end

function push_polynomial!(exprs, syms, pows, vars, newsym)
   push!(exprs, :($(newsym) = 0.0))
   for n = 1:length(vars[1])
      vars_n = [ vars[m][n] for m = 1:length(pows) ]
      s = push_monomial!(exprs, syms, pows, vars_n)
      push!(exprs, :($(newsym) += $(s)))
   end
end

@generated function invnew(x1::SVector{10, T}) where {T}
   exprs = Expr[]
   syms = Symbol[]

   # start the symbol list by adding all powes x1[n]^p
   for n = 1:10
      s = _gen_sym([1], [n])
      push!(exprs, :($s = x1[$n]))
      push!(syms, s)
      for p = 2:9
         sold = _gen_sym([p-1], [n])
         snew = _gen_sym([p], [n])
         push!(exprs, :($(snew) = $(sold) * x1[$n]))
         push!(syms, snew)
      end
   end

   #------------------------------------------------
   # Primaries
   #------------------------------------------------

   # Primary invariants for NBody=5and deg=12
    # : what goes in the function for the evaluation
   push_polynomial!(exprs, syms, (1,) , NB5I.PP1, :p1)
   push_polynomial!(exprs, syms, (1,1,) , NB5I.PP2, :p2)
   push_polynomial!(exprs, syms, (2,) , NB5I.PP3, :p3)
   push_polynomial!(exprs, syms, (1,1,1,) , NB5I.PP4, :p4)
   push_polynomial!(exprs, syms, (3,) , NB5I.PP5, :p5)
   push_polynomial!(exprs, syms, (1,1,1,1,) , NB5I.PP6, :p6)
   push_polynomial!(exprs, syms, (4,) , NB5I.PP7, :p7)
   push_polynomial!(exprs, syms, (1,1,1,1,1,) , NB5I.PP8, :p8)
   push_polynomial!(exprs, syms, (5,) , NB5I.PP9, :p9)
   push_polynomial!(exprs, syms, (6,) , NB5I.PP10, :p10)

   push_polynomial!(exprs, syms, (2,1,) , NB5I.ISS1, :is1)
   push_polynomial!(exprs, syms, (1,1,1,) , NB5I.ISS2, :is2)
   push_polynomial!(exprs, syms, (3,1,) , NB5I.ISS3, :is3)
   push_polynomial!(exprs, syms, (2,2,) , NB5I.ISS4, :is4)
   push_polynomial!(exprs, syms, (2,1,1,) , NB5I.ISS5, :is5)
   push_polynomial!(exprs, syms, (2,1,1,) , NB5I.ISS6, :is6)
   push_polynomial!(exprs, syms, (2,1,1,) , NB5I.ISS7, :is7)
   push_polynomial!(exprs, syms, (4,1,) , NB5I.ISS8, :is8)
   push_polynomial!(exprs, syms, (3,2,) , NB5I.ISS9, :is9)
   push_polynomial!(exprs, syms, (3,1,1,) , NB5I.ISS10, :is10)
   push_polynomial!(exprs, syms, (2,2,1,) , NB5I.ISS11, :is11)
   push_polynomial!(exprs, syms, (3,1,1,) , NB5I.ISS12, :is12)
   push_polynomial!(exprs, syms, (2,2,1,) , NB5I.ISS13, :is13)
   push_polynomial!(exprs, syms, (3,1,1,) , NB5I.ISS14, :is14)
   push_polynomial!(exprs, syms, (2,1,1,1,) , NB5I.ISS15, :is15)
   push_polynomial!(exprs, syms, (5,1,) , NB5I.ISS16, :is16)
   push_polynomial!(exprs, syms, (4,2,) , NB5I.ISS17, :is17)
   push_polynomial!(exprs, syms, (3,3,) , NB5I.ISS18, :is18)
   push_polynomial!(exprs, syms, (4,1,1,) , NB5I.ISS19, :is19)
   push_polynomial!(exprs, syms, (3,2,1,) , NB5I.ISS20, :is20)
   push_polynomial!(exprs, syms, (2,2,2,) , NB5I.ISS21, :is21)
   push_polynomial!(exprs, syms, (4,1,1,) , NB5I.ISS22, :is22)
   push_polynomial!(exprs, syms, (3,2,1,) , NB5I.ISS23, :is23)
   push_polynomial!(exprs, syms, (4,1,1,) , NB5I.ISS24, :is24)
   push_polynomial!(exprs, syms, (3,1,1,1,) , NB5I.ISS25, :is25)
   push_polynomial!(exprs, syms, (2,2,1,1,) , NB5I.ISS26, :is26)
   push_polynomial!(exprs, syms, (3,2,1,) , NB5I.ISS27, :is27)
   push_polynomial!(exprs, syms, (7,) , NB5I.ISS28, :is28)
   push_polynomial!(exprs, syms, (6,1,) , NB5I.ISS29, :is29)
   push_polynomial!(exprs, syms, (5,2,) , NB5I.ISS30, :is30)
   push_polynomial!(exprs, syms, (4,3,) , NB5I.ISS31, :is31)
   push_polynomial!(exprs, syms, (5,1,1,) , NB5I.ISS32, :is32)
   push_polynomial!(exprs, syms, (4,2,1,) , NB5I.ISS33, :is33)
   push_polynomial!(exprs, syms, (3,3,1,) , NB5I.ISS34, :is34)
   push_polynomial!(exprs, syms, (3,2,2,) , NB5I.ISS35, :is35)
   push_polynomial!(exprs, syms, (5,1,1,) , NB5I.ISS36, :is36)
   push_polynomial!(exprs, syms, (4,2,1,) , NB5I.ISS37, :is37)
   push_polynomial!(exprs, syms, (3,3,1,) , NB5I.ISS38, :is38)
   push_polynomial!(exprs, syms, (5,1,1,) , NB5I.ISS39, :is39)
   push_polynomial!(exprs, syms, (4,1,1,1,) , NB5I.ISS40, :is40)
   push_polynomial!(exprs, syms, (8,) , NB5I.ISS41, :is41)
   push_polynomial!(exprs, syms, (7,1,) , NB5I.ISS42, :is42)
   push_polynomial!(exprs, syms, (6,2,) , NB5I.ISS43, :is43)
   push_polynomial!(exprs, syms, (5,3,) , NB5I.ISS44, :is44)
   push_polynomial!(exprs, syms, (9,) , NB5I.ISS45, :is45)
   push_polynomial!(exprs, syms, (8,1,) , NB5I.ISS46, :is46)

   push!(exprs, :((@SVector [p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,
         is1, is2, is3, is4, is5, is6, is7, is8, is9, is10,
         is11, is12, is13, is14,
         is15, is16, is17, is18, is19, is20, is21, is22, is23,
         is24, is25, is26, is27, is28, is29, is30, is31, is32,
         is33, is34, is35, is36, is37, is38, is39, is40, is41,
         is42, is43, is44, is45, is46])))
   codeblock = Expr(:block, exprs...)

   quote
      $(Expr(:meta, :inline))
      @inbounds $codeblock
   end
end

@generated function invnew!(out::MVector{56,T}, x1::SVector{10, T}) where {T}
   exprs = Expr[]
   syms = Symbol[]

   # start the symbol list by adding all powes x1[n]^p
   for n = 1:10
      s = _gen_sym([1], [n])
      push!(exprs, :($s = x1[$n]))
      push!(syms, s)
      for p = 2:9
         sold = _gen_sym([p-1], [n])
         snew = _gen_sym([p], [n])
         push!(exprs, :($(snew) = $(sold) * x1[$n]))
         push!(syms, snew)
      end
   end

   #------------------------------------------------
   # Primaries
   #------------------------------------------------

   # Primary invariants for NBody=5and deg=12
    # : what goes in the function for the evaluation
   push_polynomial!(exprs, syms, (1,) , NB5I.PP1, :(out[1]))
   push_polynomial!(exprs, syms, (1,1,) , NB5I.PP2, :(out[2]))
   push_polynomial!(exprs, syms, (2,) , NB5I.PP3, :(out[3]))
   push_polynomial!(exprs, syms, (1,1,1,) , NB5I.PP4, :(out[4]))
   push_polynomial!(exprs, syms, (3,) , NB5I.PP5, :(out[5]))
   push_polynomial!(exprs, syms, (1,1,1,1,) , NB5I.PP6, :(out[6]))
   push_polynomial!(exprs, syms, (4,) , NB5I.PP7, :(out[7]))
   push_polynomial!(exprs, syms, (1,1,1,1,1,) , NB5I.PP8, :(out[8]))
   push_polynomial!(exprs, syms, (5,) , NB5I.PP9, :(out[9]))
   push_polynomial!(exprs, syms, (6,) , NB5I.PP10, :(out[10]))

   push_polynomial!(exprs, syms, (2,1,) , NB5I.ISS1, :(out[11]))
   push_polynomial!(exprs, syms, (1,1,1,) , NB5I.ISS2, :(out[12]))
   push_polynomial!(exprs, syms, (3,1,) , NB5I.ISS3, :(out[13]))
   push_polynomial!(exprs, syms, (2,2,) , NB5I.ISS4, :(out[14]))
   push_polynomial!(exprs, syms, (2,1,1,) , NB5I.ISS5, :(out[15]))
   push_polynomial!(exprs, syms, (2,1,1,) , NB5I.ISS6, :(out[16]))
   push_polynomial!(exprs, syms, (2,1,1,) , NB5I.ISS7, :(out[17]))
   push_polynomial!(exprs, syms, (4,1,) , NB5I.ISS8, :(out[18]))
   push_polynomial!(exprs, syms, (3,2,) , NB5I.ISS9, :(out[19]))
   push_polynomial!(exprs, syms, (3,1,1,) , NB5I.ISS10, :(out[20]))
   push_polynomial!(exprs, syms, (2,2,1,) , NB5I.ISS11, :(out[21]))
   push_polynomial!(exprs, syms, (3,1,1,) , NB5I.ISS12, :(out[22]))
   push_polynomial!(exprs, syms, (2,2,1,) , NB5I.ISS13, :(out[23]))
   push_polynomial!(exprs, syms, (3,1,1,) , NB5I.ISS14, :(out[24]))
   push_polynomial!(exprs, syms, (2,1,1,1,) , NB5I.ISS15, :(out[25]))
   push_polynomial!(exprs, syms, (5,1,) , NB5I.ISS16, :(out[26]))
   push_polynomial!(exprs, syms, (4,2,) , NB5I.ISS17, :(out[27]))
   push_polynomial!(exprs, syms, (3,3,) , NB5I.ISS18, :(out[28]))
   push_polynomial!(exprs, syms, (4,1,1,) , NB5I.ISS19, :(out[29]))
   push_polynomial!(exprs, syms, (3,2,1,) , NB5I.ISS20, :(out[30]))
   push_polynomial!(exprs, syms, (2,2,2,) , NB5I.ISS21, :(out[31]))
   push_polynomial!(exprs, syms, (4,1,1,) , NB5I.ISS22, :(out[32]))
   push_polynomial!(exprs, syms, (3,2,1,) , NB5I.ISS23, :(out[33]))
   push_polynomial!(exprs, syms, (4,1,1,) , NB5I.ISS24, :(out[34]))
   push_polynomial!(exprs, syms, (3,1,1,1,) , NB5I.ISS25, :(out[35]))
   push_polynomial!(exprs, syms, (2,2,1,1,) , NB5I.ISS26, :(out[36]))
   push_polynomial!(exprs, syms, (3,2,1,) , NB5I.ISS27, :(out[37]))
   push_polynomial!(exprs, syms, (7,) , NB5I.ISS28, :(out[38]))
   push_polynomial!(exprs, syms, (6,1,) , NB5I.ISS29, :(out[39]))
   push_polynomial!(exprs, syms, (5,2,) , NB5I.ISS30, :(out[40]))
   push_polynomial!(exprs, syms, (4,3,) , NB5I.ISS31, :(out[41]))
   push_polynomial!(exprs, syms, (5,1,1,) , NB5I.ISS32, :(out[42]))
   push_polynomial!(exprs, syms, (4,2,1,) , NB5I.ISS33, :(out[43]))
   push_polynomial!(exprs, syms, (3,3,1,) , NB5I.ISS34, :(out[44]))
   push_polynomial!(exprs, syms, (3,2,2,) , NB5I.ISS35, :(out[45]))
   push_polynomial!(exprs, syms, (5,1,1,) , NB5I.ISS36, :(out[46]))
   push_polynomial!(exprs, syms, (4,2,1,) , NB5I.ISS37, :(out[47]))
   push_polynomial!(exprs, syms, (3,3,1,) , NB5I.ISS38, :(out[48]))
   push_polynomial!(exprs, syms, (5,1,1,) , NB5I.ISS39, :(out[49]))
   push_polynomial!(exprs, syms, (4,1,1,1,) , NB5I.ISS40, :(out[50]))
   push_polynomial!(exprs, syms, (8,) , NB5I.ISS41, :(out[51]))
   push_polynomial!(exprs, syms, (7,1,) , NB5I.ISS42, :(out[52]))
   push_polynomial!(exprs, syms, (6,2,) , NB5I.ISS43, :(out[53]))
   push_polynomial!(exprs, syms, (5,3,) , NB5I.ISS44, :(out[54]))
   push_polynomial!(exprs, syms, (9,) , NB5I.ISS45, :(out[55]))
   push_polynomial!(exprs, syms, (8,1,) , NB5I.ISS46, :(out[56]))

   push!(exprs, :(return out))
   codeblock = Expr(:block, exprs...)

   quote
      $(Expr(:meta, :inline))
      @inbounds $codeblock
   end
end

end

using LinearAlgebra

const A20 = [ @SVector rand(20) for n = 1:10 ]
const A10 = [ @SVector rand(10) for n = 1:10 ]
const A56 = [ @SVector rand(56) for n = 1:10 ]
const X = [ @SVector rand(10) for n = 1:10 ]
function runN(f, A)
   s = 0.0
   for n = 1:length(A)
      s += dot(A[n], f(X[n]))
   end
   return s
end

@info("Old Invariants Code")
@btime runN($(Main.NB5I.invold), $A56)
@btime runN($(Main.NB5I.invold), $A56)

x = @SVector rand(10)
Main.NB5I.invnew(x) ≈ Main.NB5I.invold(x)

@info("New Invariants Code")
@btime runN($(Main.NB5I.invnew), $A56)
@btime runN($(Main.NB5I.invnew), $A56)

J = @MVector zeros(56)
Main.NB5I.invnew!(J, x) ≈ Main.NB5I.invold(x)
invnewin = let J=@MVector zeros(56)
   x -> Main.NB5I.invnew!(J, x)
end

@info("Inplace Invariants Code")
@btime runN($invnewin, $A56)
@btime runN($invnewin, $A56)

##

#
# @btime runN($invnewin, $x, 10)
#
# @info("New Invariants Inline")
# @btime Main.NB5I.invnew!($J, $x)
# @btime Main.NB5I.invnew!($J, $x)

# Main.NB5I.invstp(x) ≈ Main.NB5I.invold(x)
# @info("StaticPolynomials")
# @btime Main.NB5I.invstp($x)
# @btime Main.NB5I.invstp($x)
#
# Main.NB5I.invstp2(x) ≈ Main.NB5I.invold(x)
# @btime Main.NB5I.invstp2($x)
# @btime Main.NB5I.invstp2($x)

# using Profile
# f(ff, x) = (for n = 1:10_000_000; ff(x); end)
# @profile f(Main.NB5I.invstp, x)
# Profile.print()

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
using JuliaFormatter
format(
    [
        "../../src",
        "../../test",
        "../../benchmark/run.jl",
    ],
    verbose=true,
)


text:
import OptimTests: initial_x, solution_optimum, optim_problem

initial_x(p::OptimizationProblem) = p.initial_x
solver_optimum(p::OptimizationProblem) = p.minimum

function cutest_fg!(nlp, x, g)
    grad!(nlp, x, g)
    obj(nlp, x)
end

function symmetrize!(h)
    for j = 1:size(h,2)
        for i = 1:j-1
            h[i,j] = h[j,i]
        end
    end
    h
end

function cutest_hess!(nlp, x, h)
    symmetrize!(copy!(h, hess(nlp, x)))
end

function optim_problem(op::OptimizationProblem,
                       F = real(zero(eltype(initial_x(op)))),
                       G = similar(initial_x(op)),
                       H = spzeros(eltype(initial_x(op)), size(initial_x(op))...),
                       ) # TODO: How do we choose sparse vs full for H?
    if op.istwicedifferentiable
        df = TwiceDifferentiable(objective(op), gradient(op),
                                 MVP.objective_gradient(op),
                                 MVP.hessian(op),
                                 initial_x(op),
                                 F, G, H)
    elseif op.isdifferentiable
        df = OnceDifferentiable(objective(op), gradient(op),
                                MVP.objective_gradient(op),
                                initial_x(op),
                                F, G)
    else
        error("Only implemented for differentiable problems.")
    end

    return df
end

"""
Look up a the best (known) minimum of a CUTEst model.

Returns (minimizer, minimum)
"""
function minimumlookup(nlp::CUTEstModel)
    # TODO: Create approximate minimum
    # TODO: Should we store the minimizer, or just give it back a NaN?
    # TODO: Should the NaN-minimizer be of size nlp.meta.nvar?
    return [NaN], NaN
end

function optimizationproblem(nlp::CUTEstModel)
    # TODO: Set up database of minimum values for CUTEst problems
    # It seems like all the non-constrained CUTEst problems are twice differentiable?
    if nlp.meta.ncon > 0
        warn("We currently only support unconstrained CUTEst models, but $(nlp.meta.name) has `ncon > 0`.")
    end
    minimizer, minimum = minimumlookup(nlp)
    op = OptimizationProblem("$(nlp.meta.name)-$(nlp.meta.nvar)",
                             x->obj(nlp, x),
                             (g, x)->copy!(g, grad(nlp, x)),
                             (g, x)->cutest_fg!(nlp, x, g),
                             (h, x)->cutest_hess!(nlp, x, h),
                             nothing,
                             nlp.meta.x0,
                             minimizer,
                             minimum,
                             true,
                             true)
end


text:
import LinearAlgebra: dot
# we could get away with only one additional allocation if x ≠ y, but
# this is actually worse for performance and does not allow for parallelism
function dot(x::AbstractVector, A::Inverse{<:Any, <:Union{<:Cholesky, <:CholeskyPivoted}}, y::AbstractVector)
	xp = copy(x)
	yp = x ≡ y ? xp : copy(y)
	dot!!(xp, A, yp)
end

# WARNING: overwrites x and y with intermediate results
function dot!!(x::AbstractVector, A::Inverse{<:Any, <:Cholesky}, y::AbstractVector)
	C = A.parent
	inverse_cholesky_dot!!(x, C, y)
end

# WARNING: overwrites x and y with intermediate results
function dot!!(x::AbstractVector, A::Inverse{<:Any, <:CholeskyPivoted}, y::AbstractVector)
	C = A.parent
	permute!(x, C.p)
	x ≡ y || permute!(y, C.p)
	inverse_cholesky_dot!!(x, C, y)
end

# helper function for ternary dot product with inverse cholesky
function inverse_cholesky_dot!!(x, C, y)
	if x ≡ y
		inverse_cholesky_norm!(C, y)
	else
		if length(x) < parallel_threshold
			sequential_inverse_cholesky_dot!!(x, C, y)
		else
			parallel_inverse_cholesky_dot!!(x, C, y)
		end
	end
end

@inline function sequential_inverse_cholesky_dot!!(x, C, y)
	L = C.U' # since getting L causes allocations
	Lx, Ly =  ldiv!(L, x), ldiv!(L, y)
	return dot(Lx, Ly) # sum(abs2, Ly) is ever so slightly faster
end
@inline function parallel_inverse_cholesky_dot!!(x, C, y)
	L = C.U' # since getting L causes allocations
	Lx = @spawn ldiv!(L, x)
	Ly = @spawn ldiv!(L, y)
	Lx, Ly = fetch(Lx), fetch(Ly)
	return dot(Lx, Ly) # sum(abs2, Ly) is ever so slightly faster
end
# helper function for ternary dot product with inverse cholesky
@inline function inverse_cholesky_norm!(C, y)
	L = C.U' # since getting L causes allocations
	Ly = ldiv!(L, y)
	sum(abs2, Ly)
end

# IDEA: could have non-allocating mul! for this
# advantage seems to be less pronounced than for dot, since with matrices,
# BLAS3 is called and already includes parallelism
function *(X, A::Inverse{<:Any, <:Union{<:Cholesky, <:CholeskyPivoted}}, Y)
	inverse_cholesky_mul(X, A.parent, Y)
end

# saving one matrix solve if X == Y'
function inverse_cholesky_mul(X, C, Y)
	if X ≡ Y'
		Y = copy(Y)
		if C isa CholeskyPivoted
			permute_rows!(Y, C.piv)
		end
		L = C.U'
		LY = ldiv!(L, Y)
		*(LY', LY)
	else
		*(X, C\Y)
	end
end

function diag(A::Inverse{<:Any, <:Cholesky})
	inverse_cholesky_diag(A.parent)
end
function diag(A::Inverse{<:Any, <:CholeskyPivoted})
	C = A.parent
	d = inverse_cholesky_diag(C)
	invpermute!(d, C.piv)
end
# computes the diagonal of the inverse of a Cholesky factorization
function inverse_cholesky_diag(C)
	issuccess(C) || throw("cannot invert unsuccessful factorization, C.info = $(C.info)")
	L = C.U'
	invL = inv(L)
	vec(sum(abs2, invL, dims = 1))
end

################## computes diagonal of ternary product ########################
# IDEA: could dispatch when diag(LazyMatrixProduct{}) is called
function diag_mul(X, A::Inverse{<:Any, <:Cholesky}, Y)
	C = A.parent
	inverse_cholesky_diag_mul(X, C, Y)
end
function diag_mul(X, A::Inverse{<:Any, <:CholeskyPivoted}, Y)
	C = A.parent
	CY = copy(Y)
	if X ≡ Y'
		permute_rows!(CY, C.piv)
		d = inverse_cholesky_diag_norm!(C, CY)
	else
		d = inverse_cholesky_diag_mul!(X, C, CY)
	end
end
# leaves type of C unspecified to work with any implementation of a Cholesky factorization
# as long as it has a C.U field and allows for a backslash solve, and issuccess
function inverse_cholesky_diag_mul(X, C, Y)
	issuccess(C) || throw("cannot invert unsuccessful factorization, C.info = $(C.info)")
	CY = copy(Y)
	if X ≡ Y'
		inverse_cholesky_diag_norm!(C, CY)
	else
		inverse_cholesky_diag_mul!(X, C, CY)
	end
end
# WARNING: overwrites Y
function inverse_cholesky_diag_norm!(C, Y)
	issuccess(C) || throw("cannot invert unsuccessful factorization, C.info = $(C.info)")
	L = C.U'
	ldiv!(L, Y)
	vec(sum(abs2, Y, dims = 1))
end
# WARNING: overwrites Y
function inverse_cholesky_diag_mul!(X, C, Y)
	issuccess(C) || throw("cannot invert unsuccessful factorization, C.info = $(C.info)")
	ldiv!(C, Y)
	diag(X*Y)
end

########################## used for CholeskyPivoted code #######################
function permute_rows!(X::AbstractMatrix, p)
	for x in eachcol(X) # permuting rows means we have to iterate through columns
		permute!(x, p)
	end
	return X
end
function permute_columns!(X::AbstractMatrix, p)
	for x in eachrow(X) # permuting columns means we have to iterate through rows
		permute!(x, p)
	end
	return X
end


text:
module COMTRADE

using DataFrames

export read_comtrade, @S_str

# package code goes here
include("datetime.jl")
include("cfg.jl")
include("main.jl")

end # module


text:
"""
    plot_results(
        ata_model::Union{MaximinModel, SoysterMaximinModel, DeJongMaximinModel, RobustMaximinModel, CCMaximinModel, MinimaxModel};
        plots_folder = "plots",
    )

# Description

Plot the ICFs and TIFs of the assembled tests.

# Arguments

- **`ata_model::AbstractModel`** : Required. The model built with `ATA` fuctions, `ata_model.design` matrix must be `IxT` or `nfsxT` if the items are grouped by friend sets.
- **`plots_folder`** : Optional. Default: "plots". The folder in which the output is stored.
"""
function plot_results(ata_model::Union{MaximinModel, SoysterMaximinModel, DeJongMaximinModel, RobustMaximinModel, CCMaximinModel, MinimaxModel}; plots_folder = "plots")
    if !isdir(plots_folder)
        mkdir(plots_folder)
    else
        success!(
            ata_model,
            string(
                "There is already a folder with this name, files in ",
                plots_folder,
                " will be overwritten.",
            ),
        )
    end
    T = ata_model.settings.T
    if size(ata_model.output.design, 1) > 0
        n_items = ata_model.settings.n_items
        length_max = [ata_model.constraints[t].length_max for t = 1:T]
        n_fs = size(ata_model.settings.fs.items, 1)
        if n_fs < ata_model.settings.n_items
            n = n_fs * T
            design = reshape(ata_model.output.design, n_fs, T)
            new_design = zeros(Float64, n_items, T)
            for t = 1:T
                new_design[
                    vcat(
                        ata_model.settings.fs.items[findall(
                            design[:, t] .== one(Float64),
                        )]...,
                    ),
                    t,
                ] .= one(Float64)
            end
            design = new_design
            DelimitedFiles.writedlm(string(plots_folder, "/designItemLevel.csv"), design)
        else
            n = n_items * T
            design = reshape(ata_model.output.design, n_items, T)
        end

        #TIF e ICF
        if ata_model.obj.name in
           ["maximin", "soyster_maximin", "de_jong_maximin", "cc_maximin", "minimax"]
            if isfile("sim_pool.csv")
                sim_pool = CSV.read("sim_pool.csv", DataFrames.DataFrame)
            else
                sim_pool = Float64[]
            end

            ThetasPlot = collect(range(-4, stop = 4, length = 101)) #nqp values in interval/r/n",
            IIFplot = item_info(
                ata_model.settings.irt.parameters,
                ThetasPlot,
                model = ata_model.settings.irt.model,
                parametrization = ata_model.settings.irt.parametrization,
                D = ata_model.settings.irt.D,
            )
            ICFplot = item_char(
                ata_model.settings.irt.parameters,
                ThetasPlot,
                model = ata_model.settings.irt.model,
                parametrization = ata_model.settings.irt.parametrization,
                D = ata_model.settings.irt.D,
            )[1][
                :,
                :,
                1,
            ]
            IIFf = Array{Float64,2}(undef, T, 101)
            ICFf = Array{Float64,2}(undef, T, 101)

            for t = 1:T
                IIFf[t, :] = [sum(IIFplot[findall(design[:, t] .> 0), i]) for i = 1:101]
                ICFf[t, :] = [sum(ICFplot[findall(design[:, t] .> 0), i]) for i = 1:101]
            end

            ATAPlot.plot_ATA(
                ata_model,
                IIFf,
                ICFf,
                design;
                sim_pool = sim_pool,
                plots_folder = plots_folder,
            )

            if ata_model.obj.name == "cc_maximin"
                ATAPlot.plot_ATA_CC(
                    ata_model,
                    IIFf,
                    ICFf,
                    design;
                    sim_pool = sim_pool,
                    plots_folder = plots_folder,
                )
            end

        end

    else
        println("No solution found")
    end

end


text:
"""
    int2binary(n::Int64, array=zeros(Int64,0))

Converts any integer to a binary representation as initiger array with 0, 1.

# Arguments
- `n::Int64`: Integer value which has to be converted
- `array=zeros(Int64,0)`: Empty array for return

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> ClassicAlgorithmsCollections.int2binary(131)
([1,0,0,0,0,0,1,1],1)
```
"""
function int2binary(n::Int64, array = zeros(Int64, 0))
    if n > 1
        int2binary(n >> 1, array)
    end
    n &= 1
    append!(array, n)
    return array, n
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
let
    for use_autodiff in (false, true)
        for (name, prob) in Optim.UnconstrainedProblems.examples
            if prob.isdifferentiable
                f_prob = prob.f
                res = Optim.optimize(f_prob, prob.initial_x, BFGS(), Optim.Options(autodiff = use_autodiff))
                @assert norm(Optim.minimizer(res) - prob.solutions) < 1e-2
            end
        end
    end
end


text:
module BadelineBoost

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Badeline Boost" => Ahorn.EntityPlacement(
        Maple.BadelineBoost
    ),
)

Ahorn.nodeLimits(entity::Maple.BadelineBoost) = 0, -1

sprite = "objects/badelineboost/idle00.png"

function Ahorn.selection(entity::Maple.BadelineBoost)
    nodes = get(entity.data, "nodes", ())
    x, y = Ahorn.position(entity)

    res = Ahorn.Rectangle[Ahorn.getSpriteRectangle(sprite, x, y)]

    for node in nodes
        nx, ny = Int.(node)

        push!(res, Ahorn.getSpriteRectangle(sprite, nx, ny))
    end

    return res
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.BadelineBoost)
    px, py = Ahorn.position(entity)

    for node in get(entity.data, "nodes", ())
        nx, ny = Int.(node)

        theta = atan(py - ny, px - nx)
        Ahorn.drawArrow(ctx, px, py, nx + cos(theta) * 8, ny + sin(theta) * 8, Ahorn.colors.selection_selected_fc, headLength=6)
        Ahorn.drawSprite(ctx, sprite, nx, ny)

        px, py = nx, ny
    end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.BadelineBoost, room::Maple.Room)
    x, y = Ahorn.position(entity)
    Ahorn.drawSprite(ctx, sprite, x, y)
end

end


text:
# video at: https://vimeo.com/190052649

addprocs(7)
using RoME, IncrementalInference, Gadfly, Colors, KernelDensityEstimate
# for p in procs()
# # @everywhere begin
# remotecall_fetch(p, ()->using RoME)
# remotecall_fetch(p, ()->using IncrementalInference)
# remotecall_fetch(p, ()->using Gadfly)
# remotecall_fetch(p, ()->using Colors)
# # end
# end

@everywhere using RoME, RoMEPlotting
@everywhere using IncrementalInference
@everywhere using Gadfly
@everywhere using Colors
@everywhere using KernelDensityEstimate
# end

# , Gadfly

# using CloudGraphs, Neo4j
# include("BlandAuthDB.jl")
#
# configuration = CloudGraphs.CloudGraphConfiguration(dbaddress, 7474, dbusr, dbpwd, mongoaddress, 27017, false, "", "");
# cloudGraph = connect(configuration);
# # register types of interest in CloudGraphs
# registerGeneralVariableTypes!(cloudGraph)
# IncrementalInference.setCloudDataLayerAPI!()

# togglePrtStbLines()

GTp = Dict{String, Vector{Float64}}()
GTl = Dict{String, Vector{Float64}}()

GTp["l100"] = [0.0;0]
GTp["l101"] = [50.0;0]
GTp["l102"] = [100.0;0]
GTp["l103"] = [100.0;50.0]
GTp["l104"] = [100.0;100.0]
GTp["l105"] = [50.0;100.0]
GTp["l106"] = [0.0;100.0]
GTp["l107"] = [0.0;50.0]
GTp["l108"] = [0.0;-50.0]
GTp["l109"] = [0.0;-100.0]
GTp["l110"] = [50.0;-100.0]
GTp["l111"] = [100.0;-100.0]
GTp["l112"] = [100.0;-50.0]

GTl["l1"] = [10.0;30]
GTl["l2"] = [30.0;-30]
GTl["l3"] = [80.0;40]
GTl["l4"] = [120.0;-50]

cursor = [0.0;0]

function landmsInRange(GTl::Dict{String, Vector{Float64}}, cur::Vector{Float64};
			lim::Float64=100.0)
	inrange = Dict{String, Float64}()
	for gtl in GTl
		rho = norm(gtl[2]-cur)
		if rho < lim
			inrange[gtl[1]] = rho
		end
	end
	return inrange
end

function isInFG!(fgl::FactorGraph, lbl::Symbol; N=100, ready=0)
	v = nothing
	if !haskey(fgl.IDs, lbl)
		init = 300*randn(2,N)
		v = addNode!(fgl, lbl, init, diagm([1000.0;1000.0]), N=N, ready=ready)
	else
		v = getVert(fgl, lbl)
	end
	v
end

function addLandmsOnPose!(fgl::FactorGraph, pose::Graphs.ExVertex, GTl::Dict{String, Float64};
			ready=0,N=100)
	for gtl in GTl
		println("addLandmsOnPose! -- adding $(gtl[1])")
		v = isInFG!(fgl, Symbol(gtl[1]), N=N,ready=ready)
		# add the constraint
		ppr = Point2DPoint2DRange([gtl[2]], 2.0, [1.0])
		addFactor!(fgl, [pose;v], ppr, ready=ready)
	end
	nothing
end

function addNewPose!(fgl::FactorGraph, from::Symbol, lbl::Symbol, GTp;
			ready=0, N=N)

	init = 300*randn(2,N)

	v = addNode!(fgl, lbl, init, diagm([1000.0;1000.0]), N=N, ready=ready)
	rhoZ = norm(GTp[string(lbl)]-GTp[string(from)])
	ppr = Point2DPoint2DRange([rhoZ], 3.0, [1.0])
	f = addFactor!(fgl, [from,lbl], ppr, ready=ready)
  # f = addFactor!(fgl, [getVert(fgl, from);v], ppr, ready=ready)
  initializeNode!(fgl,lbl)
	# pts = evalFactor2(fgl, f, v.index)
	# setVal!(v, pts)
	# updateFullVert!(fgl,v)
	getVert(fgl, lbl)
end

function drive(fgl::FactorGraph, GTp, GTl, from, to; N=100)
	v = addNewPose!(fgl, from, to, GTp, N=N)
	# v = getVert(fgl, to)
	addLandmsOnPose!(fgl, v, landmsInRange(GTl, GTp[string(to)], lim=120.0), N=N )
  println("added landmark")
	writeGraphPdf(fgl)
	nothing
end

function batchsolve(fgl::FactorGraph; N::Int=100)
	tree = wipeBuildNewTree!(fgl, drawpdf=true)
	inferOverTree!(fgl, tree, N=N)
	nothing
end

@everywhere begin
function drawQuadLandms(fgl; file="", w=30cm, h=20cm,
			h11 = Gadfly.context(),
			h12 = Gadfly.context())
	h21 = haskey(fgl.IDs,"l3") ? drawMarginalContour(fgl,"l3") : Gadfly.context()
	h22 = haskey(fgl.IDs,"l4") ? drawMarginalContour(fgl,"l4") : Gadfly.context()
	hh = vstack(hstack(h11,h12),hstack(h21,h22))
	if length(file) > 0
		Gadfly.draw(PNG(file,w,h),hh)
	end
	hh
end

function drawLandmMargOver(fgl::FactorGraph, lbl::String,
		file,
		gtlayers;
		w=30cm, h=20cm)

	hover = Gadfly.context()
	if haskey(fgl.IDs,lbl)
		hover =  drawMarginalContour(fgl,lbl)
		# for gtl in gtlayers
		# 	push!(hover.layers, gtl)
		# end
	end

	if length(file) > 0
		Gadfly.draw(PNG(file,w,h), hover)
	end
	hover
end

function drawGroundTruth(GTp, orderp, GTl=nothing, orderl=[]; drawranges=true, interp=false, t=0)
	LAYERS = Gadfly.Layer[]
	X = Float64[]
	Y = Float64[]
	for o in orderp
		push!(X, GTp[o][1])
		push!(Y, GTp[o][2])
	end
	if interp
		X[end] = t*GTp[orderp[length(orderp)-1]][1]+(1-t)*GTp[orderp[end]][1]
		Y[end] = t*GTp[orderp[length(orderp)-1]][2]+(1-t)*GTp[orderp[end]][2]
	end
	push!(LAYERS, Gadfly.layer(x=X, y=Y, Geom.point)[1])
	for i in 2:(length(orderp))
		push!(LAYERS, Gadfly.layer(x=[X[i-1];X[i]], y=[Y[i-1];Y[i]], Geom.path())[1])
	end
	Xlr = Float64[]
	Ylr = Float64[]
	Lblr = String[]
	Xla = Float64[]
	Yla = Float64[]
	Lbla = String[]
	for o in orderl
		if o == "l1" || o == "l2"
			push!(Lblr, o)
			push!(Xlr, GTl[o][1])
			push!(Ylr, GTl[o][2])
		else
			push!(Lbla, o)
			push!(Xla, GTl[o][1])
			push!(Yla, GTl[o][2])
		end
	end
	if length(orderl) > 0
		push!(LAYERS, Gadfly.layer(x=Xlr, y=Ylr, Geom.point, Theme(default_color=colorant"red", point_size=4pt)  )[1]) # label=Lblr, Geom.label
		push!(LAYERS, Gadfly.layer(x=Xla, y=Yla, Geom.point, Theme(default_color=colorant"cyan", point_size=4pt) )[1]) # label=Lbla, Geom.label
		if drawranges
			idx = !interp ? length(X) : length(X)-1
			gtlcur = landmsInRange(GTl, [X[idx];Y[idx]])
			for lbl in keys(gtlcur)
				push!(LAYERS, Gadfly.layer(x=[X[idx];GTl[lbl][1]],y=[Y[idx];GTl[lbl][2]],Geom.path(), Theme(default_color=colorant"magenta"))[1])
			end
		end
	end
	Gadfly.plot(LAYERS,
		Guide.title("Ground Truth"),
		Coord.Cartesian(xmin=-150,xmax=150,ymin=-150, ymax=150))
end

function drawConvenience1(fgl, ii, fr, folderloc, h12, h12a, interp)
	drawQuadLandms(fgl,file="$(folderloc)/overlms$(fr).png",h12=h12);
	draw(PNG("$(folderloc)/gt$(fr).png",30cm,20cm),h12)
	ii -= interp ? 1 : 0
	lstPosePl = RoME.drawMarginalContour(fgl,"l$(100+ii)")
	ii += interp ? 1 : 0
	draw(PNG("$(folderloc)/lst$(fr).png",30cm,20cm),lstPosePl)
	push!(h12a.layers, lstPosePl.layers[1])
	draw(PNG("$(folderloc)/lstOver$(fr).png",30cm,20cm),h12a)
	nothing
end

function drawConvenience2(fgl, fr, folderloc, h12b)
	ppl = drawLandms(fgl,showmm=true,lbls=false,from=100)
	ppl.coord = Coord.Cartesian(xmin=-150,xmax=150,ymin=-150,ymax=150)

	@sync begin
		allposesfile2 = "$(folderloc)/allPoses$(fr).png"
		@async draw(PNG(allposesfile2,30cm,20cm),ppl)

		# for l in ppl.layers
		templyrs = union(h12b.layers, ppl.layers)
		h12b.layers = templyrs
		# end

		allposesfile = "$(folderloc)/allPosesOver$(fr).png"
		draw(PNG(allposesfile,30cm,20cm),h12b)
	end
	nothing
end

end #everywhere

function drawAllVidImages(GTp, GTl, fgl, ii, fr; drawranges=true, interp=false, t=0)
	poselbls = ["l$(100+i)" for i in 0:ii]
	lmlbls = ["l$(i)" for i in 1:4]
	folderloc = "/home/dehann/irosVid/new"
	rr = Dict{Int,Future}()

	h12 = drawGroundTruth(GTp, poselbls, GTl, lmlbls, drawranges=drawranges, interp=interp, t=t)
	# overlay ground truth and landmark estimation
	rr[1] = @spawn drawLandmMargOver(fgl, "l3", "$(folderloc)/lm3Over$(fr).png", h12.layers)
	rr[2] = @spawn drawLandmMargOver(fgl, "l4", "$(folderloc)/lm4Over$(fr).png", h12.layers)

	h12a = drawGroundTruth(GTp, poselbls, GTl, lmlbls, drawranges=drawranges, interp=interp, t=t)
	rr[3] = @spawn drawConvenience1(fgl, ii, fr, folderloc, h12, h12a, interp)

	h12b = drawGroundTruth(GTp, poselbls, GTl, lmlbls, drawranges=drawranges, interp=interp, t=t)
	rr[4] = @spawn drawConvenience2(fgl, fr, folderloc, h12b)

	for r in rr
		fetch(r[2])
	end

	# adding kld plots
  nothing
end

function draw30AllFast(GTp, GTl, fgl, th, offset, start, ending;
				interp=true, drawranges=true)
	# @sync begin
	# 	@async [drawAllVidImages(GTp, GTl, fgl, th, i+offset, interp=interp, drawranges=drawranges, t=(30-i)/30.0) for i in start:minimum([10;ending])];
	# 	@async [drawAllVidImages(GTp, GTl, fgl, th, i+offset, interp=interp, drawranges=drawranges, t=(30-i)/30.0) for i in 11:minimum([20;ending])];
	# 	@async [drawAllVidImages(GTp, GTl, fgl, th, i+offset, interp=interp, drawranges=drawranges, t=(30-i)/30.0) for i in 21:ending];
	# end

  # doesnt work properly yet, redrawing old
	@sync begin
		for i in start:ending
			@async drawAllVidImages(GTp, GTl, fgl, th, i+offset, interp=interp, drawranges=drawranges, t=(30-i)/30.0)
		end
	end

	nothing
end

function layerCircle(;cent=[0.0;0.0], radius=1.0, c="deepskyblue", N=200)
	TH = linspace(0,2pi,200)
	X = real(radius*exp(TH*im))+cent[1]
	Y = imag(radius*exp(TH*im))+cent[2]
	layer(x=X,y=Y, Geom.path(), Theme(default_color=parse(Colorant,c)))
end
function plotCircle(;cent=[0.0;0.0], radius=1.0, c="deepskyblue", N=200, drawcent=false)
	PL = []
	push!(PL, layerCircle(cent=cent, radius=radius, c=c, N=N) )
	!drawcent ? nothing : push!(PL, layer(x=[cent[1]],y=[cent[2]], Geom.point, Theme(default_color=parse(Colorant,c))))
	Gadfly.plot(PL...)
end

function drawFirstPoseIllustration(GTl, GTp)
	l1, l2 = GTl["l1"], GTl["l2"]
	p1, p2 = GTp["l100"], GTp["l101"]
	rho1, rho2 = norm(p1-l1), norm(p1-l2)

	PL = Gadfly.Layer[]

	PL = union(PL, plotCircle(cent=[0.0;0],radius=3.0, c="blue", drawcent=false).layers)
	PL = union(PL, plotCircle(cent=[36.0;12],radius=3.0, c="blue", drawcent=false).layers)

	PL = union(PL, plotCircle(cent=l1,radius=rho1, c="red", drawcent=true).layers)
	PL = union(PL, plotCircle(cent=l2,radius=rho2, c="red", drawcent=true).layers)

	Gadfly.plot(PL,
		Guide.title("Parametric illustration"),
		Coord.Cartesian(xmin=-150,xmax=150,ymin=-150, ymax=150) )
end

function drawSecondPoseIllustration(GTl, GTp)
	l1, l2 = GTl["l1"], GTl["l2"]
	p1, p2 = GTp["l100"], GTp["l101"]
	rho1, rho2 = norm(p1-l1), norm(p1-l2)

	PL = Gadfly.Layer[]

	PL = union(PL, plotCircle(cent=[0.0;0],radius=norm(p1-p2), c="blue", drawcent=true).layers)
	PL = union(PL, plotCircle(cent=[36.0;12],radius=norm(p1-p2), c="blue", drawcent=true).layers)

	PL = union(PL, plotCircle(cent=l1,radius=rho1, c="gray").layers)
	PL = union(PL, plotCircle(cent=l2,radius=rho2, c="gray").layers)

	Gadfly.plot(PL,
		Guide.title("Parametric illustration"),
		Coord.Cartesian(xmin=-150,xmax=150,ymin=-150, ymax=150) )
end

function drawFirstL3Illustration(GTl, GTp)
	l1, l2, l3 = GTl["l1"], GTl["l2"], GTl["l3"]
	p1, p2 = GTp["l100"], GTp["l101"]
	rho1, rho2 = norm(p1-l1), norm(p1-l2)

	PL = Gadfly.Layer[]

	PL = union(PL, plotCircle(cent=[0.0;0],radius=norm(p1-l3), c="cyan", drawcent=true).layers)
	PL = union(PL, plotCircle(cent=[36.0;12],radius=norm(p1-l3), c="cyan", drawcent=true).layers)

	PL = union(PL, plotCircle(cent=l1,radius=rho1, c="gray").layers)
	PL = union(PL, plotCircle(cent=l2,radius=rho2, c="gray").layers)

	push!(PL, layer(x=[l3[1]], y=[l3[2]], Geom.point, Theme(default_color=colorant"cyan", point_size=4pt)  )[1] )

	Gadfly.plot(PL,
		Guide.title("Parametric illustration"),
		Coord.Cartesian(xmin=-150,xmax=150,ymin=-150, ymax=150) )
end

function drawIllustrations(GTl, GTp, folderloc)
	pl = drawGroundTruth(GTp, String[], GTl, ["l$(i)" for i in 1:4], drawranges=false, interp=false, t=0)
	filename = "$(folderloc)/gtPos0.png"
	draw(PNG(filename,30cm,20cm),pl)

	# first pose location
	pl2 = drawFirstPoseIllustration(GTl, GTp)
	filename2 = "$(folderloc)/gtArc0.png"
	draw(PNG(filename2,30cm,20cm),pl2)

	pl3 = drawSecondPoseIllustration(GTl, GTp)
	filename3 = "$(folderloc)/gtPose2Pred.png"
	draw(PNG(filename3,30cm,20cm),pl3)

	pl4 = drawFirstL3Illustration(GTl, GTp)
	filename4 = "$(folderloc)/gtLm3Pred.png"
	draw(PNG(filename4,30cm,20cm),pl4)

	nothing
end

function evaluateAccuracy(fgl::FactorGraph, GTp)
  # get last pose
  xx,ll = ls(fg)
  maxPosition = 0
  for l in ll
    val = parse(Int,string(l)[2:end])
    maxPosition = val >= 100 && val > maxPosition ? val : maxPosition
  end
  vsym = Symbol("l$(maxPosition)")
  pX = getVertKDE(fg, vsym)
  truePos = GTp[string(vsym)]

  postlikeli = evaluateDualTree(pX, reshape(truePos,2,1))[1]
  MaxPointError = truePos - getKDEMax(pX, N=1000)
  return postlikeli, norm(MaxPointError)
end

function evaluateAllAccuracy(fgl::FactorGraph, GTp, GTl)
  xx,ll = ls(fg)
  LK = Float64[]
  PE = Float64[]
  for l in ll
    pX = getVertKDE(fg, l)
    truePos = haskey(GTp, string(l)) ? GTp[string(l)] : GTl[string(l)]
    push!(LK, evaluateDualTree(pX, reshape(truePos,2,1))[1])
    MaxPointError = truePos - getKDEMax(pX, N=1000)
    push!(PE, norm(MaxPointError))
  end
  return LK, PE
end

folderloc2 = "/home/dehann/irosVid/new"
drawIllustrations(GTl, GTp, folderloc2)

# N = 100

fg = initfg()

for N in collect(25:25:25)

fg = initfg()
# fg.sessionname="SESSranges"
# fg.cg = cloudGraph

# Some starting position
init = 300*randn(2,N)
v1 = addNode!(fg, :l100, init, diagm([1000.0;1000.0]), N=N, ready=0)

lmv1 = landmsInRange(GTl, GTp["l100"])
addLandmsOnPose!(fg, v1, lmv1, N=N )

# must pin landmarks for guage
pp2 = PriorPoint2D(GTl["l1"], diagm([1.0;1.0]), [1.0])
f = addFactor!(fg,[:l1], pp2)
# f = addFactor!(fg,[getVert(fg,:l1)], pp2)
pp2 = PriorPoint2D(GTl["l2"], diagm([1.0;1.0]), [1.0])
f = addFactor!(fg, [:l2], pp2)

writeGraphPdf(fg)

batchsolve(fg, N=N)
lk, pe = evaluateAccuracy(fg, GTp)
fid = open("/home/dehann/Videos/slamedonutLstPoseLks_$(N).txt","w")
write(fid, "$(lk)\n")
close(fid)
fid = open("/home/dehann/Videos/slamedonutLstPoseErr_$(N).txt","w")
write(fid, "$(pe)\n")
close(fid)

drawAllVidImages(GTp, GTl, fg, 0, 0);
draw30AllFast(GTp, GTl, fg, 1, 0, 1, 29)
# @sync begin
# @async [drawAllVidImages(GTp, GTl, fg, 1, i+0, interp=true, t=(30-i)/30.0) for i in 1:10];
# @async [drawAllVidImages(GTp, GTl, fg, 1, i+0, interp=true, t=(30-i)/30.0) for i in 11:20];
# @async [drawAllVidImages(GTp, GTl, fg, 1, i+0, interp=true, t=(30-i)/30.0) for i in 21:29];
# end
## drawAllVidImages(GTp, GTl, fg, 0, 29, interp=false)

addNewPose!(fg, :l100, :l101, GTp, N=N)
initializeNode!(fg,:l101)
writeGraphPdf(fg)
batchsolve(fg, N=N)
lk, pe = evaluateAccuracy(fg, GTp)
fid = open("/home/dehann/Videos/slamedonutLstPoseLks_$(N).txt","a")
write(fid, "$(lk)\n")
close(fid)
fid = open("/home/dehann/Videos/slamedonutLstPoseErr_$(N).txt","a")
write(fid, "$(pe)\n")
close(fid)

draw30AllFast(GTp, GTl, fg, 1, 30, 0, 14, drawranges=false, interp=false)
# [drawAllVidImages(GTp, GTl, fg, 1, 30+i, drawranges=false) for i in 0:14]

addLandmsOnPose!(fg, getVert(fg, :l101),
								landmsInRange(GTl, GTp["l101"]), N=N )
writeGraphPdf(fg)
batchsolve(fg, N=N)
lk, pe = evaluateAccuracy(fg, GTp)
fid = open("/home/dehann/Videos/slamedonutLstPoseLks_$(N).txt","a")
write(fid, "$(lk)\n")
close(fid)
fid = open("/home/dehann/Videos/slamedonutLstPoseErr_$(N).txt","a")
write(fid, "$(pe)\n")
close(fid)

draw30AllFast(GTp, GTl, fg, 1, 45, 0, 14, interp=false)
# [drawAllVidImages(GTp, GTl, fg, 1, 45+i) for i in 0:14]

fgd = deepcopy(fg)
drive(fg, GTp, GTl, :l101, :l102, N=N)
initializeNode!(fg,:l102)
writeGraphPdf(fg)

batchsolve(fg, N=N)
lk, pe = evaluateAccuracy(fg, GTp)
fid = open("/home/dehann/Videos/slamedonutLstPoseLks_$(N).txt","a")
write(fid, "$(lk)\n")
close(fid)
fid = open("/home/dehann/Videos/slamedonutLstPoseErr_$(N).txt","a")
write(fid, "$(pe)\n")
close(fid)

draw30AllFast(GTp, GTl, fgd, 2, 60, 0, 29)
# [drawAllVidImages(GTp, GTl, fg, 2, i+60, interp=true, t=(30-i)/30.0) for i in 0:29]

# initializeNode!(fg,:l102)

fgd = deepcopy(fg)
drive(fg, GTp, GTl, :l102, :l103, N=N)
batchsolve(fg, N=N)
lk, pe = evaluateAccuracy(fg, GTp)
fid = open("/home/dehann/Videos/slamedonutLstPoseLks_$(N).txt","a")
write(fid, "$(lk)\n")
close(fid)
fid = open("/home/dehann/Videos/slamedonutLstPoseErr_$(N).txt","a")
write(fid, "$(pe)\n")
close(fid)
# drawLandms(fg,showmm=true)
# h12 = drawGroundTruth(GTp, [:l$(00+i)" for i in 0:3], GTl, [:l$()" for i in 1:4])
# drawQuadLandms(fg,file="irosVid/test4.png",h12=h12);
# draw(PNG("irosVid/gt4.png",30cm,20cm),h12)
# lstPosePl = RoME.drawMarginalContour(fg,:l103)
# draw(PNG("irosVid/lst4.png",30cm,20cm),lstPosePl)
# push!(h12.layers, lstPosePl.layers[1])
# draw(PNG("irosVid/lstOver4.png",30cm,20cm),h12)
# drawAllVidImages(GTp, GTl, fg, 3, 4)

draw30AllFast(GTp, GTl, fgd, 3, 90, 0, 29)
# [drawAllVidImages(GTp, GTl, fg, 3, i+90, interp=true, t=(30-i)/30.0) for i in 0:29]

fgd = deepcopy(fg)
drive(fg, GTp, GTl, :l103, :l104, N=N)
# addLandmsOnPose!(fg, getVert(fg, :l104),
# 								landmsInRange(GTl, GTp["l104"]), N=N )
# writeGraphPdf(fg)
batchsolve(fg, N=N)
lk, pe = evaluateAccuracy(fg, GTp)
fid = open("/home/dehann/Videos/slamedonutLstPoseLks_$(N).txt","a")
write(fid, "$(lk)\n")
close(fid)
fid = open("/home/dehann/Videos/slamedonutLstPoseErr_$(N).txt","a")
write(fid, "$(pe)\n")
close(fid)

draw30AllFast(GTp, GTl, fgd, 4, 120, 0, 29)
# [drawAllVidImages(GTp, GTl, fg, 4, i+120, interp=true, t=(30-i)/30.0) for i in 0:29]

fgd = deepcopy(fg)
drive(fg, GTp, GTl, :l104, :l105, N=N)
# addLandmsOnPose!(fg, getVert(fg, :l104),
# 								landmsInRange(GTl, GTp["l104"]), N=N )
writeGraphPdf(fg)

batchsolve(fg, N=N)
lk, pe = evaluateAccuracy(fg, GTp)
fid = open("/home/dehann/Videos/slamedonutLstPoseLks_$(N).txt","a")
write(fid, "$(lk)\n")
close(fid)
fid = open("/home/dehann/Videos/slamedonutLstPoseErr_$(N).txt","a")
write(fid, "$(pe)\n")
close(fid)

draw30AllFast(GTp, GTl, fgd, 5, 150, 0, 29)
# [drawAllVidImages(GTp, GTl, fg, 5, i+150, interp=true, t=(30-i)/30.0) for i in 0:29]

fgd = deepcopy(fg)
drive(fg, GTp, GTl, :l105, :l106, N=N)
batchsolve(fg, N=N)
lk, pe = evaluateAccuracy(fg, GTp)
fid = open("/home/dehann/Videos/slamedonutLstPoseLks_$(N).txt","a")
write(fid, "$(lk)\n")
close(fid)
fid = open("/home/dehann/Videos/slamedonutLstPoseErr_$(N).txt","a")
write(fid, "$(pe)\n")
close(fid)

# drawAllVidImages(GTp, GTl, fgd, 6, 7)
draw30AllFast(GTp, GTl, fgd, 6, 180, 0, 29)
# [drawAllVidImages(GTp, GTl, fg, 6, i+180, interp=true, t=(30-i)/30.0) for i in 0:29]

fgd = deepcopy(fg)
drive(fg, GTp, GTl, :l106, :l107, N=N)
batchsolve(fg, N=N)
lk, pe = evaluateAccuracy(fg, GTp)
fid = open("/home/dehann/Videos/slamedonutLstPoseLks_$(N).txt","a")
write(fid, "$(lk)\n")
close(fid)
fid = open("/home/dehann/Videos/slamedonutLstPoseErr_$(N).txt","a")
write(fid, "$(pe)\n")
close(fid)

draw30AllFast(GTp, GTl, fgd, 7, 210, 0, 29)
# [drawAllVidImages(GTp, GTl, fg, 7, i+210, interp=true, t=(30-i)/30.0) for i in 0:29]

fgd = deepcopy(fg)
drive(fg, GTp, GTl, :l107, :l108, N=N)
# addLandmsOnPose!(fg, getVert(fg, :l108),
# 								landmsInRange(GTl, GTp["l108"]), N=N )
writeGraphPdf(fg)

batchsolve(fg, N=N)
lk, pe = evaluateAccuracy(fg, GTp)
fid = open("/home/dehann/Videos/slamedonutLstPoseLks_$(N).txt","a")
write(fid, "$(lk)\n")
close(fid)
fid = open("/home/dehann/Videos/slamedonutLstPoseErr_$(N).txt","a")
write(fid, "$(pe)\n")
close(fid)

draw30AllFast(GTp, GTl, fgd, 8, 240, 0, 29)
# [drawAllVidImages(GTp, GTl, fg, 8, i+240, interp=true, t=(30-i)/30.0) for i in 0:29]

fgd = deepcopy(fg)
drive(fg, GTp, GTl, :l108, :l109, N=N)
# addLandmsOnPose!(fg, getVert(fg, :l109),
# 								landmsInRange(GTl, GTp["l109"]), N=N )
writeGraphPdf(fg)

batchsolve(fg, N=N)
lk, pe = evaluateAccuracy(fg, GTp)
fid = open("/home/dehann/Videos/slamedonutLstPoseLks_$(N).txt","a")
write(fid, "$(lk)\n")
close(fid)
fid = open("/home/dehann/Videos/slamedonutLstPoseErr_$(N).txt","a")
write(fid, "$(pe)\n")
close(fid)

draw30AllFast(GTp, GTl, fgd, 9, 270, 0, 29)
# [drawAllVidImages(GTp, GTl, fg, 9, i+270, interp=true, t=(30-i)/30.0) for i in 0:29]

fgd = deepcopy(fg)
drive(fg, GTp, GTl, :l109, :l110, N=N)
# addLandmsOnPose!(fg, getVert(fg, :l110),
# 								landmsInRange(GTl, GTp["l110"]), N=N )
writeGraphPdf(fg)

batchsolve(fg, N=N)
lk, pe = evaluateAccuracy(fg, GTp)
fid = open("/home/dehann/Videos/slamedonutLstPoseLks_$(N).txt","a")
write(fid, "$(lk)\n")
close(fid)
fid = open("/home/dehann/Videos/slamedonutLstPoseErr_$(N).txt","a")
write(fid, "$(pe)\n")
close(fid)

draw30AllFast(GTp, GTl, fgd, 10, 300, 0, 29)
# [drawAllVidImages(GTp, GTl, fg, 10, i+300, interp=true, t=(30-i)/30.0) for i in 0:29]

fgd = deepcopy(fg)
drive(fg, GTp, GTl, :l110, :l111, N=N)
# addLandmsOnPose!(fg, getVert(fg, :l111),
# 								landmsInRange(GTl, GTp["l111"]), N=N )
writeGraphPdf(fg)

batchsolve(fg, N=N)
lk, pe = evaluateAccuracy(fg, GTp)
fid = open("/home/dehann/Videos/slamedonutLstPoseLks_$(N).txt","a")
write(fid, "$(lk)\n")
close(fid)
fid = open("/home/dehann/Videos/slamedonutLstPoseErr_$(N).txt","a")
write(fid, "$(pe)\n")
close(fid)

draw30AllFast(GTp, GTl, fgd, 11, 330, 0, 29)
# [drawAllVidImages(GTp, GTl, fg, 11, i+330, interp=true, t=(30-i)/30.0) for i in 0:29]

fgd = deepcopy(fg)
drive(fg, GTp, GTl, :l111, :l112, N=N)

tsec = @elapsed batchsolve(fg, N=N)
fid = open("/home/dehann/Videos/slamedonutBatchTime_$(N).txt","w")
write(fid, "$(tsec)")
close(fid)

lk, pe = evaluateAccuracy(fg, GTp)
fid = open("/home/dehann/Videos/slamedonutLstPoseLks_$(N).txt","a")
write(fid, "$(lk)\n")
close(fid)
fid = open("/home/dehann/Videos/slamedonutLstPoseErr_$(N).txt","a")
write(fid, "$(pe)\n")
close(fid)

LK, PE = evaluateAllAccuracy(fg, GTp, GTl)

fid = open("/home/dehann/Videos/slamedonutAllPE_$(N).txt","w")
for pe in PE
  write(fid, "$(pe)\n")
end
close(fid)
fid = open("/home/dehann/Videos/slamedonutAllLK_$(N).txt","w")
for lk in LK
  write(fid, "$(lk)\n")
end
close(fid)

writeGraphPdf(fg)
wipeBuildNewTree!(fg, drawpdf=true)

draw30AllFast(GTp, GTl, fgd, 12, 360, 0, 29)
# [drawAllVidImages(GTp, GTl, fg, 12, i+360, interp=true, t=(30-i)/30.0) for i in 0:29]

draw30AllFast(GTp, GTl, fg, 12, 390, 0, 29, interp=false)
# [drawAllVidImages(GTp, GTl, fg, 12, i+390) for i in 0:29]

# run(`convert -delay 100 irosVid/test*.png irosVid/anim.gif`)

run(`ffmpeg -y -i $(folderloc2)/lstOver%d.png -threads 8 -vcodec libx264 -s 1920x1080 -b:v 2M /home/dehann/Videos/slamedonutLastPose_$(N).mp4`)
# run(`ffmpeg -y -i $(folderloc2)/lm3Over%d.png -vcodec libx264 -s 1920x1080 -b:v 2M /home/dehann/Videos/slamedonutLm3_$(N).mp4`)
# run(`ffmpeg -y -i $(folderloc2)/lm4Over%d.png -vcodec libx264 -s 1920x1080 -b:v 2M /home/dehann/Videos/slamedonutLm4_$(N).mp4`)

end

#


text:
@testset "Product with non-commutative monomials and variables" begin
    Mod.@ncpolyvar x y z
    m1 = x * (y * x)
    @test variables(m1)[1] == x
    @test variables(m1)[2] == y
    @test variables(m1)[3] == x
    @test all(exponents(m1) .== 1)
    m2 = y * (y * x)
    @test variables(m2)[1] == y
    @test variables(m2)[2] == x
    @test collect(exponents(m2)) == [2, 1]
    for m3 in [(x*x)*(y^2*x), (x*x*y)*(y*x), x*(x*y^2*x), (x^2*y^2)*x]
        @test variables(m3)[1] == x
        @test variables(m3)[2] == y
        @test variables(m3)[3] == x
        @test collect(exponents(m3)) == [2, 2, 1]
    end
    mv = monovec([x, y, x*y])
    @test x * mv == [x^2, x*y, x^2*y]
    @test mv * x == [x^2, y*x, x*y*x]
    @test y^2 * mv == [y^2*x, y^3, y^2*x*y]
    @test mv * y^2 == [x*y^2, y^3, x*y^3]
    m4 = y*x^2
    @test variables(m4)[1] == y
    @test variables(m4)[2] == x
    @test collect(exponents(m4)) == [1, 2]

    @test (x*y) * (3x^2 + 2x*y + 4y^2) == 3x*y*x^2 + 2x*y*x*y + 4x*y^3

    P = [2 3 4;
         3 4 5;
         4 5 6]
    Q = [4 3 5;
         3 2 4
         5 4 6]
    p = @inferred polynomial(P, [x*y, x^2, y^2])
    @test coefficients(p) == [4, 3, 5, 4, 3, 2, 6, 5, 4]
    @test monomials(p) == [x^4, x^3*y, x^2*y^2, x*y^3, x*y*x^2, x*y*x*y, y^4, y^2*x^2, y^2*x*y]
    p = @inferred polynomial(Q, monovec([x*y, x^2, y^2]))
    @test coefficients(p) == [4, 3, 5, 4, 3, 2, 6, 5, 4]
    @test monomials(p) == [x^4, x^3*y, x^2*y^2, x*y^3, x*y*x^2, x*y*x*y, y^4, y^2*x^2, y^2*x*y]
    q = @inferred polynomial(zeros(Float64, 0, 0), emptymonovec(x))
    @test q isa AbstractPolynomialLike{Float64}
    @test q == 0
end
@testset "Noncommutative quadratic" begin
    Mod.@ncpolyvar x[1:2]
    Q = Hermitian([1 2 + 3im; 2 - 3im 4])
    p = 1x[1]^2 + (2 + 3im) * x[1] * x[2] + (2 - 3im) * x[2] * x[1] + 4x[2]^2
    @test polynomial(Q, x) == p
    @test polynomial(Q, monovec(x)) == p
end


text:
"""
    GoldenSection <: UnivariateAlgorithm

The golden section search algorithm optimize a univariate function by iteratively shrinking
the initial bracket around the minimum with the golden ratio as an approximation to the
Fibonacci sequence.
"""
struct GoldenSection <: UnivariateAlgorithm end

function _optimize(f, lower::T, upper::T, alg::GoldenSection;
        reltol, abstol, maxiter) where {T}
    p = φ - 1
    x = T(p * lower + (1 - p) * upper)
    yx = f(x)
    converged = false
    iter = 1
    while true
        x_tol = reltol * abs(x) + abstol
        if abs(upper - lower) < 2x_tol
            converged = true
            break
        end
        iter == maxiter && break
        iter += 1
        new_x = T(p * upper + (1 - p) * lower)
        new_yx = f(new_x)
        if yx > new_yx
            lower, x, yx = x, new_x, new_yx
        else
            lower, upper = new_x, lower
        end
    end
    return Solution(converged, iter, x, yx)
end


text:
# Julia wrapper for header: /usr/include/scip/cons_symresack.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConshdlrSymresack(scip)
    ccall((:SCIPincludeConshdlrSymresack, libscip), SCIP_RETCODE, (Ptr{SCIP_},), scip)
end

function SCIPcreateSymbreakCons(scip, cons, name, perm, vars, nvars, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateSymbreakCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{Cint}, Ptr{Ptr{SCIP_VAR}}, Cint, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, perm, vars, nvars, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsSymresack(scip, cons, name, perm, vars, nvars, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
    ccall((:SCIPcreateConsSymresack, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{Cint}, Ptr{Ptr{SCIP_VAR}}, Cint, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32, UInt32), scip, cons, name, perm, vars, nvars, initial, separate, enforce, check, propagate, _local, modifiable, dynamic, removable, stickingatnode)
end

function SCIPcreateConsBasicSymresack(scip, cons, name, perm, vars, nvars)
    ccall((:SCIPcreateConsBasicSymresack, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONS}}, Cstring, Ptr{Cint}, Ptr{Ptr{SCIP_VAR}}, Cint), scip, cons, name, perm, vars, nvars)
end


text:
function update!(t::DAQTask)
    #synchronize Julia DAQTask to current properties
    #of the associated TaskHandle in C land
end

function Base.push!(t::DAQTask, chan::DAQChannel; kwargs...)
end

function Base.push!(t::DAQTask, dev::DAQDevice; kwargs...)
end

function start(task::DAQTask)
    DAQmx.StartTask(task.handle) |> catch_error
    return println("Running $(task.name)...")
end

function stop(task::DAQTask)
    DAQmx.StopTask(task.handle) |> catch_error
    if length(task.conditions) > 0
        while length(task.conditions) > 0
            x = pop!(task.conditions)
            ccall(:uv_async_send, Nothing, (Ptr{Nothing},), x)
            close(x)
        end
        println("Stopped $(task.name) and closed async conditions...")
    else
        println("Stopped $(task.name)...")
    end
    return
end

function clear(task::DAQTask)
    DAQmx.ClearTask(task.handle) |> catch_error
    task = nothing
    return
end

function isrunning(task::DAQTask)
    istaskdone = Ref{UInt32}()
    DAQmx.IsTaskDone(task.handle, istaskdone) |> catch_error
    return !(Bool(istaskdone[]))
end

function refresh!(task::DAQTask)
    DAQmx.ClearTask(task.handle) |> catch_error
    # rebuild task from fields of DAQTask and replace TaskHandle
    return
end


text:
if !isdefined(Main, :CustomStencil)
    include("../module/CustomStencil.jl")
    using Main.CustomStencil
end
include("misc.jl")
radius = 4
l = 2*radius + 1
coefs = zeros(Float64, (l,l,l))

for i = 1:radius
    global radius
    posi = radius + i + 1
    negi = radius -i + 1
    mid = radius + 1
    coefs[posi,mid,mid] = 1
    coefs[negi,mid,mid] = 1
    coefs[mid, posi, mid] = 1
    coefs[mid, negi, mid] = 1
    coefs[mid, mid, posi] = 1
    coefs[mid, mid, negi] = 1

end
mid = radius + 1
coefs[mid,mid,mid] = 1
#
# star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1, 4, c[i]*(
#             D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
#             D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))

st_inst = CustomStencil.NewStencilInstance(coefs)

nx = 5
ny = 5
nz = 5
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

t_steps = 5
## Actual time steps
data = CreateData(dx,dy,dz,radius)
dout = ApplyStencil(st_inst, data, t_steps)
#heatmap(dout1[:,:,Int(round(dz/2))])
o = Array(dout)

function cpu_star_stencil(data, radius,coefs)
    out = zeros(Float32,size(data))
    for i = (radius+1):(size(data,1)-radius)
        for j = (radius+1):(size(data,2)-radius)
            for k = (radius+1):(size(data,3)-radius)
                c = coefs[1]*data[i,j,k]
                for r = 1:radius
                    c += coefs[r+1]*(data[i + r,j,k] + data[i - r,j,k] +
                    data[i,j + r,k] + data[i,j - r,k] +
                    data[i,j,k + r] + data[i,j,k - r])
                end
                println(c)
                out[i,j,k] = c
            end
        end
    end
    return out
end

c_coefs = Float64.([1 1 1 1 1])
cout = cpu_star_stencil(data, radius, c_coefs)
sum(cout)
for i = 2:t_steps
    global cout
    cout = cpu_star_stencil(cout, radius, c_coefs)
end

@assert(sum(abs.(cout .- o))./(reduce(*, size(data))) < 0.001)
mean_err = sum(abs.(cout .- o))./(reduce(*, size(data)))

using Plots
shift = 4
heatmap(cout[:,:,Int(round(dz/2))+shift])
heatmap(o[:,:,Int(round(dz/2))+shift])
sum(cout)
heatmap(o[:,:,Int(round(dz/2))+shift] .- cout[:,:,Int(round(dz/2))+shift])

diff =  abs.(cout .- o)#./(reduce(*, size(data)))
findmax(diff)