text:
# The following is largely inspired from JuMP/test/JuMPExtension.jl

@enum Level BOTH LOWER_ONLY UPPER_ONLY DUAL_OF_LOWER

abstract type AbstractBilevelModel <: JuMP.AbstractModel end

Base.broadcastable(model::AbstractBilevelModel) = Ref(model)

mutable struct ConstraintInfo{T<:Union{Float64, Vector{Float64}}}
    start::T
    upper::T
    lower::T
    function ConstraintInfo{Float64}()
        new(NaN, NaN, NaN)
    end
    function ConstraintInfo{Vector{Float64}}(N::Integer)
        new(fill(NaN, N), fill(NaN, N), fill(NaN, N))
    end
end

mutable struct VariableInfo{T<:Union{Float64, Vector{Float64}}}
    upper::T
    lower::T
    function VariableInfo{Float64}()
        new(NaN, NaN)
    end
    function VariableInfo{Vector{Float64}}(N::Integer)
        new(fill(NaN, N), fill(NaN, N))
    end
end

mutable struct BilevelModel <: AbstractBilevelModel
    # Structured data
    upper::JuMP.AbstractModel
    lower::JuMP.AbstractModel

    # Model data
    nextvaridx::Int                                 # Next variable index is nextvaridx+1
    variables::Dict{Int, JuMP.AbstractVariable}     # Map varidx -> variable
    varnames::Dict{Int, String}                     # Map varidx -> name
    var_level::Dict{Int, Level}
    var_upper::Dict{Int, JuMP.AbstractVariableRef}
    var_lower::Dict{Int, JuMP.AbstractVariableRef}
    var_info::Dict{Int, VariableInfo}

    # upper level decisions that are "parameters" of the second level
    upper_to_lower_link::Dict{JuMP.AbstractVariableRef, JuMP.AbstractVariableRef}
    # lower level decisions that are input to upper
    lower_to_upper_link::Dict{JuMP.AbstractVariableRef, JuMP.AbstractVariableRef}
    # lower level decisions that are input to upper
    upper_var_to_lower_ctr_link::Dict{JuMP.AbstractVariableRef, JuMP.ConstraintRef}
    # joint link
    link::Dict{JuMP.AbstractVariableRef, JuMP.AbstractVariableRef}

    nextconidx::Int                                 # Next constraint index is nextconidx+1
    constraints::Dict{Int, JuMP.AbstractConstraint} # Map conidx -> variable
    connames::Dict{Int, String}                     # Map varidx -> name
    ctr_level::Dict{Int, Level}
    ctr_upper::Dict{Int, JuMP.ConstraintRef}
    ctr_lower::Dict{Int, JuMP.ConstraintRef}
    ctr_info::Dict{Int, ConstraintInfo}

    upper_objective_sense::MOI.OptimizationSense
    upper_objective_function::JuMP.AbstractJuMPScalar

    lower_objective_sense::MOI.OptimizationSense
    lower_objective_function::JuMP.AbstractJuMPScalar

    # solution data
    solver#::MOI.ModelLike
    upper_to_sblm
    lower_to_sblm
    lower_dual_to_sblm
    sblm_to_solver
    lower_primal_dual_map

    objdict::Dict{Symbol, Any}    # Same that JuMP.Model's field `objdict`

    function BilevelModel()

        model = new(
            JuMP.Model(),
            JuMP.Model(),

            # var
            0, Dict{Int, JuMP.AbstractVariable}(),   Dict{Int, String}(),    # Model Variables
            Dict{Int, Level}(), Dict{Int, JuMP.AbstractVariable}(), Dict{Int, JuMP.AbstractVariable}(),
            Dict{Int, VariableInfo}(),
            # links
            Dict{JuMP.AbstractVariable, JuMP.AbstractVariable}(), Dict{JuMP.AbstractVariable, JuMP.AbstractVariable}(),
            Dict{JuMP.AbstractVariable, JuMP.ConstraintRef}(),
            Dict{JuMP.AbstractVariable, JuMP.AbstractVariable}(),
            #ctr
            0, Dict{Int, JuMP.AbstractConstraint}(), Dict{Int, String}(),    # Model Constraints
            Dict{Int, Level}(), Dict{Int, JuMP.AbstractConstraint}(), Dict{Int, JuMP.AbstractConstraint}(),
            Dict{Int, ConstraintInfo}(),
            #obj
            MOI.FEASIBILITY_SENSE,
            zero(JuMP.GenericAffExpr{Float64, BilevelVariableRef}), # Model objective
            MOI.FEASIBILITY_SENSE,
            zero(JuMP.GenericAffExpr{Float64, BilevelVariableRef}), # Model objective

            nothing,
            nothing,
            nothing,
            nothing,
            nothing,
            nothing,
            Dict{Symbol, Any}(),
            )

        return model
    end
end

abstract type InnerBilevelModel <: AbstractBilevelModel end
struct UpperModel <: InnerBilevelModel
    m::BilevelModel
end
Upper(m::BilevelModel) = UpperModel(m)
UpperToLower(m::BilevelModel) = UpperModel(m)
struct LowerModel <: InnerBilevelModel
    m::BilevelModel
end
Lower(m::BilevelModel) = LowerModel(m)
LowerToUpper(m::BilevelModel) = LowerModel(m)
bilevel_model(m::InnerBilevelModel) = m.m
mylevel_model(m::UpperModel) = bilevel_model(m).upper
mylevel_model(m::LowerModel) = bilevel_model(m).lower
level(m::LowerModel) = LOWER_ONLY
level(m::UpperModel) = UPPER_ONLY
mylevel_ctr_list(m::LowerModel) = bilevel_model(m).ctr_lower
mylevel_ctr_list(m::UpperModel) = bilevel_model(m).ctr_upper
mylevel_var_list(m::LowerModel) = bilevel_model(m).var_lower
mylevel_var_list(m::UpperModel) = bilevel_model(m).var_upper

# obj

mylevel_obj_sense(m::LowerModel) = bilevel_model(m).lower_objective_sense
mylevel_obj_function(m::LowerModel) = bilevel_model(m).lower_objective_function
mylevel_obj_sense(m::UpperModel) = bilevel_model(m).upper_objective_sense
mylevel_obj_function(m::UpperModel) = bilevel_model(m).upper_objective_function

set_mylevel_obj_sense(m::LowerModel, val) = bilevel_model(m).lower_objective_sense = val
set_mylevel_obj_function(m::LowerModel, val) = bilevel_model(m).lower_objective_function = val
set_mylevel_obj_sense(m::UpperModel, val) = bilevel_model(m).upper_objective_sense = val
set_mylevel_obj_function(m::UpperModel, val) = bilevel_model(m).upper_objective_function = val

function set_link!(m::UpperModel, upper::JuMP.AbstractVariableRef, lower::JuMP.AbstractVariableRef)
    bilevel_model(m).upper_to_lower_link[upper] = lower
    bilevel_model(m).link[upper] = lower
    nothing
end
function set_link!(m::LowerModel, upper::JuMP.AbstractVariableRef, lower::JuMP.AbstractVariableRef)
    bilevel_model(m).lower_to_upper_link[lower] = upper
    bilevel_model(m).link[upper] = lower
    nothing
end

# Models to deal with variables that are not exchanged between models
abstract type SingleBilevelModel <: AbstractBilevelModel end
struct UpperOnlyModel <: SingleBilevelModel
    m::BilevelModel
end
UpperOnly(m::BilevelModel) = UpperOnlyModel(m)
struct LowerOnlyModel <: SingleBilevelModel
    m::BilevelModel
end
LowerOnly(m::BilevelModel) = LowerOnlyModel(m)

bilevel_model(m::SingleBilevelModel) = m.m
mylevel_model(m::UpperOnlyModel) = bilevel_model(m).upper
mylevel_model(m::LowerOnlyModel) = bilevel_model(m).lower
level(m::LowerOnlyModel) = LOWER_ONLY
level(m::UpperOnlyModel) = UPPER_ONLY
mylevel_var_list(m::LowerOnlyModel) = bilevel_model(m).var_lower
mylevel_var_list(m::UpperOnlyModel) = bilevel_model(m).var_upper

in_upper(l::Level) = l == BOTH || l == UPPER_ONLY || l == DUAL_OF_LOWER
in_lower(l::Level) = l == BOTH || l == LOWER_ONLY

#### Model ####

# Variables
struct BilevelVariableRef <: JuMP.AbstractVariableRef
    model::BilevelModel # `model` owning the variable
    idx::Int       # Index in `model.variables`
    level::Level
end
mylevel(v::BilevelVariableRef) = v.level
in_level(v::BilevelVariableRef, level::Level) = (
    v.level === BOTH ||
    v.level === level ||
    (v.level === DUAL_OF_LOWER && level === UPPER_ONLY))

in_upper(v::BilevelVariableRef) = in_upper(mylevel(v))
in_lower(v::BilevelVariableRef) = in_lower(mylevel(v))
upper_ref(v::BilevelVariableRef) = v.model.var_upper[v.idx]
lower_ref(v::BilevelVariableRef) = v.model.var_lower[v.idx]

function bound_ref(v::BilevelVariableRef)
    if mylevel(v) == LOWER_ONLY
        return lower_ref(v)
    elseif mylevel(v) == UPPER_ONLY || mylevel(v) == DUAL_OF_LOWER
        return upper_ref(v)
    elseif mylevel(v) == BOTH
        m = v.model
        i = v.idx
        upper_ref = m.var_upper[i]
        # this checks actual owner:
        if haskey(m.upper_to_lower_link, upper_ref)
            return upper_ref
        else
            return m.var_lower[i]
        end
    else
        error("Unknown level")
    end
end

function solver_ref(v::BilevelVariableRef)
    m = v.model
    if mylevel(v) == LOWER_ONLY
        return m.sblm_to_solver[
            m.lower_to_sblm[JuMP.index(lower_ref(v))]]
    else
        return m.sblm_to_solver[
            m.upper_to_sblm[JuMP.index(upper_ref(v))]]
    end
end

Base.broadcastable(v::BilevelVariableRef) = Ref(v)
Base.copy(v::BilevelVariableRef) = v
# TODO
# Base.copy(v::BilevelVariableRef, new_model::BilevelModel) = BilevelVariableRef(new_model, v.idx, v.level) # also requires a map
Base.:(==)(v::BilevelVariableRef, w::BilevelVariableRef) =
    v.model === w.model && v.idx == w.idx && v.level == w.level
JuMP.owner_model(v::BilevelVariableRef) = v.model
JuMP.isequal_canonical(v::BilevelVariableRef, w::BilevelVariableRef) = v == w
JuMP.variable_type(::AbstractBilevelModel) = BilevelVariableRef
# add in BOTH levels
function JuMP.add_variable(inner::InnerBilevelModel, v::JuMP.AbstractVariable, name::String="")
    m = bilevel_model(inner)
    m.nextvaridx += 1
    vref = BilevelVariableRef(m, m.nextvaridx, BOTH)
    # break info so that bounds go to correct level
    var_upper, var_lower = split_variable(inner, v)
    v_upper = JuMP.add_variable(m.upper, var_upper, name)
    m.var_upper[vref.idx] = v_upper
    v_lower = JuMP.add_variable(m.lower, var_lower, name)
    m.var_lower[vref.idx] = v_lower
    m.var_level[vref.idx] = BOTH
    set_link!(inner, v_upper, v_lower)
    m.variables[vref.idx] = v # save complete data
    JuMP.set_name(vref, name)
    m.var_info[vref.idx] = empty_info(v)
    vref
end
function JuMP.add_variable(single::SingleBilevelModel, v::JuMP.AbstractVariable, name::String="")
    m = bilevel_model(single)
    m.nextvaridx += 1
    vref = BilevelVariableRef(m, m.nextvaridx, level(single))
    v_level = JuMP.add_variable(mylevel_model(single), v, name)
    mylevel_var_list(single)[vref.idx] = v_level
    m.var_level[vref.idx] = level(single)
    m.variables[vref.idx] = v
    JuMP.set_name(vref, name)
    m.var_info[vref.idx] = empty_info(v)
    vref
end
function MOI.delete!(m::AbstractBilevelModel, vref::BilevelVariableRef)
    error("No deletion on bilevel models")
    delete!(m.variables, vref.idx)
    delete!(m.varnames, vref.idx)
end
MOI.is_valid(m::BilevelModel, vref::BilevelVariableRef) = vref.idx in keys(m.variables)
JuMP.num_variables(m::BilevelModel) = length(m.variables)
JuMP.num_variables(m::InnerBilevelModel) = JuMP.num_variables(bilevel_model(m))
function empty_info(::JuMP.AbstractVariable)
    return VariableInfo{Float64}()
end

"""
Split variable because actual owner of the variable should be the one holding the bounds.
"""
function split_variable(::UpperModel, v::JuMP.AbstractVariable)
    var_upper = v
    var_lower = JuMP.ScalarVariable(JuMP.VariableInfo(
        false, NaN,
        false, NaN,
        false, NaN,
        v.info.has_start, v.info.start,
        false, false))
    return var_upper, var_lower
end
function split_variable(::LowerModel, v::JuMP.AbstractVariable)
    var_lower = v
    var_upper = JuMP.ScalarVariable(JuMP.VariableInfo(
        false, NaN,
        false, NaN,
        false, NaN,
        v.info.has_start, v.info.start,
        false, false))
    return var_upper, var_lower
end

# Internal function
variable_info(vref::BilevelVariableRef) = vref.model.variables[vref.idx].info
function update_variable_info(vref::BilevelVariableRef, info::JuMP.VariableInfo)
    vref.model.variables[vref.idx] = JuMP.ScalarVariable(info)
end

function JuMP.has_lower_bound(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        return !isnan(get_dual_lower_bound(get_constrain_ref(vref)))
    end
    return variable_info(vref).has_lb
end
function JuMP.lower_bound(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        return get_dual_lower_bound(get_constrain_ref(vref))
    end
    @assert !JuMP.is_fixed(vref)
    variable_info(vref).lower_bound
end
function JuMP.set_lower_bound(vref::BilevelVariableRef, lower)
    if mylevel(vref) == DUAL_OF_LOWER
        set_dual_lower_bound(get_constrain_ref(vref), lower)
        return vref.model.variables[vref.idx]
    end
    info = variable_info(vref)
    JuMP.set_lower_bound(bound_ref(vref), lower)
    update_variable_info(vref,
                         JuMP.VariableInfo(true, lower,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end
function JuMP.delete_lower_bound(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        set_dual_lower_bound(get_constrain_ref(vref), NaN)
        return vref.model.variables[vref.idx]
    end
    info = variable_info(vref)
    JuMP.delete_lower_bound(bound_ref(vref))
    update_variable_info(vref,
                         JuMP.VariableInfo(false, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end
function JuMP.has_upper_bound(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        return !isnan(get_dual_upper_bound(get_constrain_ref(vref)))
    end
    return variable_info(vref).has_ub
end
function JuMP.upper_bound(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        return get_dual_upper_bound(get_constrain_ref(vref))
    end
    @assert !JuMP.is_fixed(vref)
    variable_info(vref).upper_bound
end
function JuMP.set_upper_bound(vref::BilevelVariableRef, upper)
    if mylevel(vref) == DUAL_OF_LOWER
        set_dual_upper_bound(get_constrain_ref(vref), upper)
        return vref.model.variables[vref.idx]
    end
    info = variable_info(vref)
    JuMP.set_upper_bound(bound_ref(vref), upper)
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           true, upper,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end
function JuMP.delete_upper_bound(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        set_dual_upper_bound(get_constrain_ref(vref), NaN)
        return vref.model.variables[vref.idx]
    end
    info = variable_info(vref)
    JuMP.delete_upper_bound(bound_ref(vref))
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           false, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end
JuMP.is_fixed(vref::BilevelVariableRef) = variable_info(vref).has_fix
JuMP.fix_value(vref::BilevelVariableRef) = variable_info(vref).fixed_value
function JuMP.fix(vref::BilevelVariableRef, value)
    if mylevel(vref) == DUAL_OF_LOWER
        error("Dual variable cannot be fixed.")
    end
    info = variable_info(vref)
    JuMP.fix(bound_ref(vref), value)
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           true, value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end
function JuMP.unfix(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        error("Dual variable cannot be fixed.")
    end
    info = variable_info(vref)
    JuMP.unfix(bound_ref(vref))
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           false, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, info.integer))
end
function JuMP.start_value(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        return get_dual_start(get_constrain_ref(vref))
    end
    variable_info(vref).start
end
function JuMP.set_start_value(vref::BilevelVariableRef, start)
    if mylevel(vref) == DUAL_OF_LOWER
        set_dual_start(get_constrain_ref(vref), start)
        return vref.model.variables[vref.idx]
    end
    info = variable_info(vref)
    in_upper(vref) && JuMP.set_start_value(upper_ref(vref), start)
    in_lower(vref) && JuMP.set_start_value(lower_ref(vref), start)
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           true, start,
                                           info.binary, info.integer))
end
JuMP.is_binary(vref::BilevelVariableRef) = variable_info(vref).binary
function JuMP.set_binary(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        error("Dual variable cannot be binary.")
    end
    @assert !JuMP.is_integer(vref)
    info = variable_info(vref)
    JuMP.set_binary(bound_ref(vref))
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           true, info.integer))
end
function JuMP.unset_binary(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        error("Dual variable cannot be binary.")
    end
    info = variable_info(vref)
    JuMP.unset_binary(bound_ref(vref))
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           false, info.integer))
end
JuMP.is_integer(vref::BilevelVariableRef) = variable_info(vref).integer
function JuMP.set_integer(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        error("Dual variable cannot be integer.")
    end
    @assert !JuMP.is_binary(vref)
    info = variable_info(vref)
    JuMP.set_integer(bound_ref(vref))
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, true))
end
function JuMP.unset_integer(vref::BilevelVariableRef)
    if mylevel(vref) == DUAL_OF_LOWER
        error("Dual variable cannot be integer.")
    end
    info = variable_info(vref)
    JuMP.unset_integer(bound_ref(vref))
    update_variable_info(vref,
                         JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                           info.has_ub, info.upper_bound,
                                           info.has_fix, info.fixed_value,
                                           info.has_start, info.start,
                                           info.binary, false))
end

# Constraints
struct BilevelConstraintRef
    model::BilevelModel # `model` owning the constraint
    idx::Int       # Index in `model.constraints`
    level::Level
end
JuMP.constraint_type(::AbstractBilevelModel) = BilevelConstraintRef
my_level(cref::BilevelConstraintRef) = cref.level
function JuMP.add_constraint(m::BilevelModel, c::JuMP.AbstractConstraint, name::String="")
    error(
        "Can't add constraint directly to the bilevel model `m`, "*
        "attach the constraint to the upper or lower model "*
        "with @constraint(Upper(m), ...) or @constraint(Lower(m), ...)")
end
# function constraint_object(con_ref::ConstraintRef{Model, _MOICON{FuncType, SetType}}) where
#     {FuncType <: MOI.AbstractScalarFunction, SetType <: MOI.AbstractScalarSet}
#     model = con_ref.model
#     f = MOI.get(model, MOI.ConstraintFunction(), con_ref)::FuncType
#     s = MOI.get(model, MOI.ConstraintSet(), con_ref)::SetType
#     return ScalarConstraint(jump_function(model, f), s)
# end
# JuMP.add_constraint(m::UpperModel, c::JuMP.VectorConstraint, name::String="") =
#     error("no vec ctr")
function JuMP.add_constraint(m::InnerBilevelModel, c::Union{JuMP.ScalarConstraint{F,S},JuMP.VectorConstraint{F,S}}, name::String="") where {F,S}
    blm = bilevel_model(m)
    blm.nextconidx += 1
    cref = BilevelConstraintRef(blm, blm.nextconidx, level(m))
    func = JuMP.jump_function(c)
    # @show c
    # @show func
    level_func = replace_variables(func, bilevel_model(m), mylevel_model(m), mylevel_var_list(m), level(m))
    level_c = JuMP.build_constraint(error, level_func, c.set)
    level_cref = JuMP.add_constraint(mylevel_model(m), level_c, name)
    blm.ctr_level[cref.idx] = level(m)
    mylevel_ctr_list(m)[cref.idx] = level_cref
    blm.constraints[cref.idx] = c
    blm.ctr_info[cref.idx] = empty_info(c)
    JuMP.set_name(cref, name)
    cref
end
function MOI.delete!(m::AbstractBilevelModel, cref::BilevelConstraintRef)
    error("can't delete")
    delete!(m.constraints, cref.idx)
    delete!(m.connames, cref.idx)
end
MOI.is_valid(m::BilevelModel, cref::BilevelConstraintRef) = cref.idx in keys(m.constraints)
MOI.is_valid(m::InnerBilevelModel, cref::BilevelConstraintRef) =
    MOI.is_valid(bilevel_model(m), cref) && bilevel_model(m).ctr_level[cref.idx] == level(m)
function JuMP.constraint_object(cref::BilevelConstraintRef, F::Type, S::Type)
    c = cref.model.constraints[cref.idx]
    # `TypeError` should be thrown is `F` and `S` are not correct
    # This is needed for the tests in `constraints.jl`
    c.func::F
    c.set::S
    c
end
function empty_info(c::JuMP.ScalarConstraint{F,S}) where {F,S}
    return ConstraintInfo{Float64}()
end
function empty_info(c::JuMP.VectorConstraint{F,S}) where {F,S}
    return ConstraintInfo{Vector{Float64}}(MOI.dimension(c.set))
end

function set_dual_start(cref::BilevelConstraintRef, value::T) where T<:Number
    cref.model.ctr_info[cref.idx].start = value
end
function set_dual_start(cref::BilevelConstraintRef, value::T) where T<:Vector{S} where S
    array = cref.model.ctr_info[cref.idx].start
    @assert length(array) == length(value)
    copyto!(array, value)
end
function get_dual_start(cref::BilevelConstraintRef)
    cref.model.ctr_info[cref.idx].start
end
function set_dual_upper_bound(cref::BilevelConstraintRef, value::T) where T<:Number
    cref.model.ctr_info[cref.idx].upper = value
end
function set_dual_upper_bound(cref::BilevelConstraintRef, value::T) where T<:Vector{S} where S
    array = cref.model.ctr_info[cref.idx].upper
    @assert length(array) == length(value)
    copyto!(array, value)
end
function get_dual_upper_bound(cref::BilevelConstraintRef)
    cref.model.ctr_info[cref.idx].upper
end
function set_dual_lower_bound(cref::BilevelConstraintRef, value::T) where T<:Number
    cref.model.ctr_info[cref.idx].lower = value
end
function set_dual_lower_bound(cref::BilevelConstraintRef, value::T) where T<:Vector{S} where S
    array = cref.model.ctr_info[cref.idx].lower
    @assert length(array) == length(value)
    copyto!(array, value)
end
function get_dual_lower_bound(cref::BilevelConstraintRef)
    cref.model.ctr_info[cref.idx].lower
end

function set_primal_upper_bound_hint(vref::BilevelVariableRef, value::T) where T<:Number
    vref.model.var_info[vref.idx].upper = value
end
function get_primal_upper_bound_hint(vref::BilevelVariableRef)
    vref.model.var_info[vref.idx].upper
end
function set_primal_lower_bound_hint(vref::BilevelVariableRef, value::T) where T<:Number
    vref.model.var_info[vref.idx].lower = value
end
function get_primal_lower_bound_hint(vref::BilevelVariableRef)
    vref.model.var_info[vref.idx].lower
end

function JuMP.value(cref::BilevelConstraintRef)
    if my_level(cref) == BilevelJuMP.LOWER_ONLY
        # Constraint index on the lower model
        con_lower_idx = cref.model.ctr_lower[cref.idx].index
        # Single bilevel model constraint associated with the lower level constraint
        con_sblm_idx = cref.model.lower_to_sblm[con_lower_idx]
    else
        # Constraint index on the lower model
        con_upper_idx = cref.model.ctr_upper[cref.idx].index
        # Single bilevel model constraint associated with the lower level constraint
        con_sblm_idx = cref.model.upper_to_sblm[con_upper_idx]
    end
    # Solver constraint associated with the single bilevel model constraint
    con_solver_idx = cref.model.sblm_to_solver[con_sblm_idx]
    return MOI.get(cref.model.solver, MOI.ConstraintPrimal(), con_solver_idx)
end
# variables again (duals)
# code for using dual variables associated with lower level constraints
# in the upper level

struct DualOf
    ci::BilevelConstraintRef
end
struct DualVariableInfo
    info::JuMP.VariableInfo
    ci::BilevelConstraintRef
end
function JuMP.build_variable(
    _error::Function,
    info::JuMP.VariableInfo,
    dual_of::DualOf;
    extra_kw_args...,
)
    for (kwarg, _) in extra_kw_args
        _error("Unrecognized keyword argument $kwarg")
    end

    if info.has_lb
        set_dual_lower_bound(dual_of.ci, info.lower_bound)
        # info.has_lb = false
        # info.lower_bound = NaN
    end
    if info.has_ub
        set_dual_upper_bound(dual_of.ci, info.upper_bound)
        # info.has_ub = false
        # info.upper_bound = NaN
    end
    info.has_fix   && _error("Dual variable does not support fixing")
    if info.has_start
        set_dual_start(dual_of.ci, info.start)
        # info.has_start = false
        # info.start = NaN
    end
    info.binary    && _error("Dual variable cannot be binary")
    info.integer   && _error("Dual variable cannot be integer")

    info = JuMP.VariableInfo(false, NaN, false, NaN,
                             false, NaN, false, NaN,
                             false, false)

    if dual_of.ci.level != LOWER_ONLY
        error("Variables can only be tied to LOWER level constraints, got $(dual_of.ci.level) level")
    end

    return DualVariableInfo(
        info,
        dual_of.ci
    )
end
function JuMP.add_variable(inner::UpperModel, dual_info::DualVariableInfo, name::String="")
    # TODO vector version
    m = bilevel_model(inner)
    m.nextvaridx += 1
    vref = BilevelVariableRef(m, m.nextvaridx, DUAL_OF_LOWER)
    v_upper = JuMP.add_variable(m.upper, JuMP.ScalarVariable(dual_info.info), name)
    m.var_upper[vref.idx] = v_upper
    m.var_level[vref.idx] = DUAL_OF_LOWER
    m.upper_var_to_lower_ctr_link[v_upper] = m.ctr_lower[dual_info.ci.idx] # TODO improve this
    m.variables[vref.idx] = JuMP.ScalarVariable(dual_info.info)
    JuMP.set_name(vref, name)
    m.var_info[vref.idx] = empty_info(dual_info)
    vref
end
function empty_info(::DualVariableInfo)
    return VariableInfo{Float64}()
end

function get_constrain_ref(vref::BilevelVariableRef)
    model = vref.model
    ctr_ref = model.upper_var_to_lower_ctr_link[model.var_upper[vref.idx]]
    idx = -1
    for (ind, ref) in model.ctr_lower
        if ref == ctr_ref
            idx = ind
        end
    end
    @assert idx != -1
    return BilevelConstraintRef(model, idx, DUAL_OF_LOWER)
end

# Objective
function JuMP.set_objective(m::InnerBilevelModel, sense::MOI.OptimizationSense,
                            f::JuMP.AbstractJuMPScalar)
    set_mylevel_obj_sense(m, sense)
    set_mylevel_obj_function(m, f)
    level_f = replace_variables(f, bilevel_model(m), mylevel_model(m), mylevel_var_list(m), level(m))
    JuMP.set_objective(mylevel_model(m), sense, level_f)
end
JuMP.objective_sense(m::InnerBilevelModel) = mylevel_obj_sense(m)
JuMP.objective_function_type(m::InnerBilevelModel) = typeof(mylevel_obj_function(m))
JuMP.objective_function(m::InnerBilevelModel) = mylevel_obj_function(m)
function JuMP.objective_function(m::InnerBilevelModel, FT::Type)
    mylevel_obj_function(m) isa FT || error("The objective function is not of type $FT")
    mylevel_obj_function(m)
end

# todo remove
JuMP.objective_sense(m::AbstractBilevelModel) = MOI.FEASIBILITY_SENSE
# end todo remove
JuMP.num_variables(m::AbstractBilevelModel) = JuMP.num_variables(bilevel_model(m))
JuMP.show_constraints_summary(::Any, ::AbstractBilevelModel) = "no summary"
JuMP.show_backend_summary(::Any, ::AbstractBilevelModel) = "no summary"
JuMP.object_dictionary(m::BilevelModel) = m.objdict
JuMP.object_dictionary(m::AbstractBilevelModel) = JuMP.object_dictionary(bilevel_model(m))
JuMP.show_objective_function_summary(::IO, ::AbstractBilevelModel) = "no summary"
# TODO
# function JuMP.constraints_string(print_mode, model::MyModel)

bileve_obj_error() = error("There is no objective for BilevelModel use Upper(.) and Lower(.)")

function JuMP.set_objective(m::BilevelModel, sense::MOI.OptimizationSense,
    f::JuMP.AbstractJuMPScalar)
    bileve_obj_error()
end
JuMP.objective_sense(m::BilevelModel) = JuMP.objective_sense(m.upper)#bileve_obj_error()
JuMP.objective_function_type(model::BilevelModel) = bileve_obj_error()
JuMP.objective_function(model::BilevelModel) = bileve_obj_error()
function JuMP.objective_function(model::BilevelModel, FT::Type)
    bileve_obj_error()
end

# Names
JuMP.name(vref::BilevelVariableRef) = vref.model.varnames[vref.idx]
function JuMP.set_name(vref::BilevelVariableRef, name::String)
    vref.model.varnames[vref.idx] = name
end
JuMP.name(cref::BilevelConstraintRef) = cref.model.connames[cref.idx]
function JuMP.set_name(cref::BilevelConstraintRef, name::String)
    cref.model.connames[cref.idx] = name
end
# TODO
# function JuMP.variable_by_name(model::MyModel, name::String)
# see jump extensions
# function JuMP.constraint_by_name(model::MyModel, name::String)

# replace variables
function replace_variables(var::BilevelVariableRef,
    model::BilevelModel,
    inner::JuMP.AbstractModel,
    variable_map::Dict{Int, V},
    level::Level) where {V<:JuMP.AbstractVariableRef}
    if var.model === model && in_level(var, level)
        return variable_map[var.idx]
    elseif var.model === model
        error("Variable $(var) belonging Only to $(var.level) level, was added in the $(level) level.")
    else
        error("A BilevelModel cannot have expression using variables of a BilevelModel different from itself")
    end
end
function replace_variables(aff::JuMP.GenericAffExpr{C, BilevelVariableRef},
    model::BilevelModel,
    inner::JuMP.AbstractModel,
    variable_map::Dict{Int, V},
    level::Level) where {C,V<:JuMP.AbstractVariableRef}
    result = JuMP.GenericAffExpr{C, JuMP.VariableRef}(0.0)#zero(aff)
    result.constant = aff.constant
    for (coef, var) in JuMP.linear_terms(aff)
        JuMP.add_to_expression!(result,
        coef,
        replace_variables(var, model, model, variable_map, level))
    end
    return result
end
function replace_variables(quad::JuMP.GenericQuadExpr{C, BilevelVariableRef},
    model::BilevelModel,
    inner::JuMP.AbstractModel,
    variable_map::Dict{Int, V},
    level::Level) where {C,V<:JuMP.AbstractVariableRef}
    aff = replace_variables(quad.aff, model, model, variable_map, level)
    quadv = JuMP.GenericQuadExpr{C, JuMP.VariableRef}(aff)
    for (coef, var1, var2) in JuMP.quad_terms(quad)
        JuMP.add_to_expression!(quadv,
        coef,
        replace_variables(var1, model, model, variable_map, level),
        replace_variables(var2, model, model, variable_map, level))
    end
    return quadv
end
replace_variables(funcs::Vector, args...) = map(f -> replace_variables(f, args...), funcs)

function print_lp(m, name)
    dest = MOI.FileFormats.Model(format = MOI.FileFormats.FORMAT_MOF)
    MOI.copy_to(dest, m)
    MOI.write_to_file(dest, name)
end

JuMP.optimize!(::T) where {T<:AbstractBilevelModel} =
    error("cant solve a model of type: $T ")
function JuMP.optimize!(model::BilevelModel, optimizer, mode::BilevelSolverMode{T} = SOS1Mode();
    lower_prob = "", upper_prob = "", bilevel_prob = "", solver_prob = "") where T

    if !MOI.is_empty(optimizer)
        error("An empty optimizer must be provided")
    end

    upper = JuMP.backend(model.upper)
    lower = JuMP.backend(model.lower)

    if length(lower_prob) > 0
        print_lp(lower, lower_prob)
    end
    if length(upper_prob) > 0
        print_lp(upper, upper_prob)
    end

    moi_upper = JuMP.index.(
        collect(values(model.upper_to_lower_link)))
    moi_link = JuMP.index(model.link)
    moi_link2 = index2(model.upper_var_to_lower_ctr_link)

    single_blm, upper_to_sblm, lower_to_sblm, lower_primal_dual_map, lower_dual_to_sblm =
    build_bilevel(upper, lower, moi_link, moi_upper, mode, moi_link2)

    # pass lower level dual variables info (start, upper, lower)
    for (idx, info) in model.ctr_info
        if haskey(model.ctr_lower, idx)
            ctr = model.ctr_lower[idx]
            pre_duals = lower_primal_dual_map.primal_con_dual_var[JuMP.index(ctr)] # vector
            duals = map(x->lower_dual_to_sblm[x], pre_duals)
            pass_dual_info(single_blm, duals, info)
        end
    end
    for (idx, info) in model.var_info
        if haskey(model.var_lower, idx)
            var = lower_to_sblm[JuMP.index(model.var_lower[idx])]
        elseif haskey(model.var_upper, idx)
            var = upper_to_sblm[JuMP.index(model.var_upper[idx])]
        else
            continue
        end
        pass_primal_info(single_blm, var, info)
    end

    if length(bilevel_prob) > 0
        print_lp(single_blm, bilevel_prob)
    end
    # print_lp(single_blm, "bilevel_orig.mof")

    solver = optimizer#MOI.Bridges.full_bridge_optimizer(optimizer, Float64)
    sblm_to_solver = MOI.copy_to(solver, single_blm, copy_names = true)
    # print_lp(solver, "bilevel_bridge.mof")
    # print_lp(solver.model, "bilevel_cache.mof")

    if length(solver_prob) > 0
        print_lp(solver, solver_prob)
    end

    MOI.optimize!(solver)

    model.solver  = solver#::MOI.ModelLike
    model.upper_to_sblm = upper_to_sblm
    model.lower_to_sblm = lower_to_sblm
    model.lower_dual_to_sblm = lower_dual_to_sblm
    model.lower_primal_dual_map = lower_primal_dual_map
    model.sblm_to_solver = sblm_to_solver

    return nothing
end
function pass_dual_info(single_blm, dual, info::ConstraintInfo{Float64})
    if !isnan(info.start)
        MOI.set(single_blm, MOI.VariablePrimalStart(), dual[], info.start)
    end
    if !isnan(info.upper) &&
        !MOI.is_valid(single_blm, CI{SVF,LT{Float64}}(dual[].value))
        MOI.add_constraint(single_blm,
            SVF(dual[]), LT{Float64}(info.upper))
    end
    if !isnan(info.lower) &&
        !MOI.is_valid(single_blm, CI{SVF,GT{Float64}}(dual[].value))
        MOI.add_constraint(single_blm,
            SVF(dual[]), GT{Float64}(info.lower))
    end
    return
end
function pass_dual_info(single_blm, dual, info::ConstraintInfo{Vector{Float64}})
    for i in eachindex(dual)
        if !isnan(info.start[i])
            MOI.set(single_blm, MOI.VariablePrimalStart(), dual[i], info.start[i])
        end
        if !isnan(info.upper[i]) &&
            !MOI.is_valid(single_blm, CI{SVF,LT{Float64}}(dual[i].value))
            MOI.add_constraint(single_blm,
                SVF(dual[i]), LT{Float64}(info.upper[i]))
        end
        if !isnan(info.lower[i]) &&
            !MOI.is_valid(single_blm, CI{SVF,GT{Float64}}(dual[i].value))
            MOI.add_constraint(single_blm,
                SVF(dual[i]), MOI.GreaterThan{Float64}(info.lower[i]))
        end
    end
    return
end

function pass_primal_info(single_blm, primal, info::VariableInfo{Float64})
    if !isnan(info.upper) &&
        !MOI.is_valid(single_blm, CI{SVF,LT{Float64}}(primal.value))
        MOI.add_constraint(single_blm,
            SVF(primal), LT{Float64}(info.upper))
    end
    if !isnan(info.lower) &&
        !MOI.is_valid(single_blm, CI{SVF,GT{Float64}}(primal.value))
        MOI.add_constraint(single_blm,
            SVF(primal), GT{Float64}(info.lower))
    end
    return
end

function JuMP.index(d::Dict)
    ret = Dict{VI,VI}()
    # sizehint!(ret, length(d))
    for (k,v) in d
        ret[JuMP.index(k)] = JuMP.index(v)
    end
    return ret
end
function index2(d::Dict)
    ret = Dict{VI,CI}()
    # sizehint!(ret, length(d))
    for (k,v) in d
        ret[JuMP.index(k)] = JuMP.index(v)
    end
    return ret
end

function JuMP.value(v::BilevelVariableRef)::Float64
    m = owner_model(v)
    solver = m.solver
    ref = solver_ref(v)
    return MOI.get(solver, MOI.VariablePrimal(), ref)
end

function JuMP.dual(cref::BilevelConstraintRef)
    # Right now this code assumes there is no possibility for vectorized constraints
    if my_level(cref) == BilevelJuMP.LOWER_ONLY
        # Constraint index on the lower model
        con_lower_idx = cref.model.ctr_lower[cref.idx].index
        # Dual variable associated with constraint index
        model_var_idxs = cref.model.lower_primal_dual_map.primal_con_dual_var[con_lower_idx]
        # Single bilevel model variable associated with the dual variable
        sblm_var_idxs = MOI.VariableIndex[]
        for vi in model_var_idxs
            push!(sblm_var_idxs, cref.model.lower_dual_to_sblm[vi])
        end
        # Solver variable associated withe the sblm model
        solver_var_idxs = MOI.VariableIndex[]
        for vi in sblm_var_idxs
            push!(solver_var_idxs, cref.model.sblm_to_solver[vi])
        end
        return MOI.get(cref.model.solver, MOI.VariablePrimal(), solver_var_idxs)
    else
        error("Dual solutions of upper level constraints are not available")
    end
end
function JuMP.primal_status(model::BilevelModel)
    return MOI.get(model.solver, MOI.PrimalStatus())
end
JuMP.dual_status(model::BilevelModel) = error("dual status cant be queried for BilevelModel")
function JuMP.termination_status(model::BilevelModel)
    return MOI.get(model.solver, MOI.TerminationStatus())
end
function JuMP.raw_status(model::BilevelModel)
    return MOI.get(model.solver, MOI.RawStatusString())
end
function JuMP.objective_value(model::BilevelModel)
    return MOI.get(model.solver, MOI.ObjectiveValue())
end
function lower_objective_value(model::BilevelModel)
    # Create a dict with lower variables in the lower objective
    vals = Dict()
    for v in model.lower_objective_function.terms.keys
        vals[MOI.VariableIndex(v.idx)] = JuMP.value(v)
    end
    # Evaluate the lower objective expression
    return MOIU.eval_variables(vi -> vals[vi],
        model.lower.moi_backend.model_cache.model.objective)
end


text:
function HelpMarker(desc)
    CImGui.TextDisabled("(?)")
    if (CImGui.IsItemHovered())
        CImGui.BeginTooltip()
        CImGui.PushTextWrapPos(CImGui.GetFontSize() * 35.0)
        CImGui.TextUnformatted(desc)
        CImGui.PopTextWrapPos()
        CImGui.EndTooltip()
    end
end

ext: jl


text:
file_path_csv = "../../../../../../externals/Core.Math.Data/data/HolisticWare.Core.Sample.Data/CSV/Pejcic_318.csv"

#------------------------------------------------------------------------------------
# low level file API
# https://en.wikibooks.org/wiki/Introducing_Julia/Working_with_text_files
file_handle = open(file_path_csv)

totaltime, totallines = open(file_path_csv) do file_handle_1
    linecounter = 0
    timetaken = @elapsed for l in eachline(file_handle_1)
        linecounter += 1
    end
    (timetaken, linecounter)
end

close(file_handle)

open(file_path_csv) do file_handle_1
    line = 1
    while !eof(file_handle_1)
      x = readline(file_handle_1)
      println("$line $x")
      line += 1
    end
  end

  open(file_path_csv) do file_handle_1
    for i in enumerate(eachline(file_handle_1))
        println(i[1], ": ", i[2])
      end
  end

  data_as_string = read(file_path_csv, String)

#------------------------------------------------------------------------------------


text:
using Documenter, HTLogger

makedocs(sitename="HTLogger Documentation")

deploydocs( repo = "github.com/MLackner/HTLogger.jl" )


text:
# Unfortunately, considering only horizontal and vertical lines doesn't give
# you the full picture; you need to also consider diagonal lines.
#
# Because of the limits of the hydrothermal vent mapping system, the lines in
# your list will only ever be horizontal, vertical, or a diagonal line at
# exactly 45 degrees. In other words:
#
# An entry like 1,1 -> 3,3 covers points 1,1, 2,2, and 3,3.
# An entry like 9,7 -> 7,9 covers points 9,7, 8,8, and 7,9.
#
# Considering all lines from the above example would now produce the following
# diagram:
#
# 1.1....11.
# .111...2..
# ..2.1.111.
# ...1.2.2..
# .112313211
# ...1.2....
# ..1...1...
# .1.....1..
# 1.......1.
# 222111....
#
# You still need to determine the number of points where at least two lines
# overlap. In the above example, this is still anywhere in the diagram with a 2
# or larger - now a total of 12 points.
#
# Consider all of the lines. At how many points do at least two lines overlap?

map = zeros(1000, 1000)
for l in eachline("input.txt")
  # parse the line
  (x1, y1, x2, y2) = parse.(Int, split(l, r",| -> "))

  # add one to each point
  if x1 == x2 || y1 == y2
    # make sure x1<x2, and y1<y2 - it simplifies the loop below
    x1 > x2 && ((x1,x2) = (x2,x1))
    y1 > y2 && ((y1,y2) = (y2,y1))

    for x in x1:x2, y in y1:y2
      map[x, y] += 1
    end
  else
    # diagonal
    for i in 0:abs(x2-x1)
      map[x1+copysign(i, x2-x1), y1+copysign(i, y2-y1)] += 1
    end
  end
end

# count points where at least two lines overlap
println("Result: ", count(>(1), map))


text:
module EGraphs

include("../docstrings.jl")

using RuntimeGeneratedFunctions
RuntimeGeneratedFunctions.init(@__MODULE__)

import ..Rule
import ..getrhsfun

using ..TermInterface
using ..Util

include("enode.jl")
export ENode

include("abstractanalysis.jl")
export AbstractAnalysis

include("eclass.jl")
export EClass
export hasdata
export getdata
export setdata!

include("egg.jl")
export find
export geteclass
export ariety
export EGraph
export merge!
export addexpr!
export addanalysis!
export rebuild!

include("analysis.jl")
export analyze!

include("ematch.jl")
include("Schedulers/Schedulers.jl")
export Schedulers
using .Schedulers

include("saturation_report.jl")
include("saturation_params.jl")
export SaturationParams
include("saturation.jl")
export saturate!
include("equality.jl")
export areequal
export @areequal
export @areequalg

include("extraction.jl")
export extract!
export ExtractionAnalysis
export astsize
export astsize_inv
export @extract

end


text:
using CUDA;
#using JSON;

#dicttxt = JSON.parsefile("config_example.json")

#x2 = dicttxt["x"]
#y2 = dicttxt["y"]
N = 256*2
x_d = CUDA.fill(convert(Float64,3), N)  # a vector stored on the GPU filled with 1.0 (Float32)
y_d = CUDA.fill(convert(Float64,2), N)
z_d = CUDA.fill(convert(Float64,1), N)# a vector stored on the GPU filled with 2.0

function myadd(x::Float32, y::Float32)

    return x+2*y
end

function myadd(x::Float64, y::Float64)
    return x+y
end

using Test
function gpu_add3!(y, x, z)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    z[index] = myadd(y[index], x[index])

    #@cuprintln(blockIdx().x, "   ", threadIdx().x)

    if threadIdx().x == 1
        @cuprintln("y=", y[index], "  x=", x[index], "  z=",  z[index])
    end

    return
end

numblocks = ceil(Int, N/256)

@cuda threads=256 blocks=numblocks gpu_add3!(y_d, x_d,z_d);
#synchronize()
display(z_d)


text:
@testset "008" begin
    @test tdpe.e008(4) == 5832
end

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
include("dataset.jl")
include("timespans.jl")
include("serialization.jl")
include("printing.jl")
include(joinpath(dirname(@__DIR__), "examples", "flac.jl"))
include(joinpath(dirname(@__DIR__), "examples", "tour.jl"))


text:
"""
  Resources

Server resource abstractions.
"""
module Resources

using HTTP
using HTTP: Response, Request, Router

using ..HttpErrors: HttpError, bad_request, internal_server_error, method_not_allowed
using ..HttpMethods: HttpMethod
using ..HttpMethods

export DeserializeError, Resource, deserialize, process, serialize, handle, handler

"""
    Resource Interface

Abstract type respresenting a resource.

# Implementation
Given `S <: HttpMethod`, an abstract type `T <: Resource` which implements this
interface *must* implement the following methods:
- `deserialize(::S, ::T, req::Request)::U`
- `process(::S, ::T, data::U)::V`
- `serialize(::S, ::T, data::V)::Response`
(mathematically, `U`, `V` are intended to represent type varables bound consistently
across all three methods, although this is informal notation here, and there are no
formal enforcements of this principle).

Implementing the first three methods will enable users to use the function
`handle(::S, ::T, req::Request)` which serially invokes `deserialize`, `process`,
and then `serialize`, wrapping each in a small layer of extra logic.

The function given by `handler(::S, ::T)` serially invokes `deserialize`, `process`,
and then `serialize`.

Note that Julia will not throw a compile-time error if the `Resource` interface is
only partially implemented; it is up to the programmer to ensure that resources are
completely written if they want to make use of `handler`.

This issue may be solved in the future by an resource interface DSL, i.e. macro.

If one wants to register a resource and method with a `Rest.Middleware.Routers.Router`
instance, then one needs to additionally implement `path(::T)`, which specifies the
path at which a resource resides.
"""
abstract type Resource end

# DeserializeError

"""
    DeserializeError

Error to throw during deserialization.
"""
struct DeserializeError <: Exception end

# path

"""
    path(::Resource)

Specifiy the path at which a resource resides.
"""
function path end

# deserialize

"""
    deserialize(method::HttpMethod, resource::Resource, req::Request)

Deserialize an API request.

The result will be piped into `x -> API.process(method, resource, x)` for processing and
then subsequently into `x -> API.serialize(method, resource, x)` for serialization into a
response which will be returned to the client.

Throwing `API.DeserializeError` inside of `deserialize` will interrupt request handling
and send a 400-code response to the client. Any other error will send a 500-code
response to the client.
"""
function deserialize end

# process

"""
    process(method::HttpMethod, resource::Resource, data)

Process data in the handling an API request.

The result will be piped into `x -> API.serialize(method, resource, x)` for serialization into a
response which will be returned to the client.

Any error thrown in `process` will interrupt request handling and send a 500-code
response to the client.
"""
function process end

# serialize

"""
    serialize(method::HttpMethod, resource::Resource, data)::HTTP.Response

Serialize an object into an HTTP response.

Any error thrown in `serialize` will interrupt request handling and send a 500-code
response to the client.
"""
function serialize end

# handler

"""
    handle(method::HttpMethod, resource::Resource, req::Request)

Handle an `HTTP.Request` for `resource` using `method`.
"""
function handle(method::HttpMethod, resource::Resource, req::Request)
  req_data = try
    deserialize(method, resource, req)
  catch err
    # TODO: Determine if DeserializeError is useful
    """
    if isa(err, DeserializeError)
      bad_request(err.msg)
    elseif isa(err, HttpError)
      rethrow(err)
    else
      sprint(showerror, err) |> internal_server_error
    end
    """
    if isa(err, HttpError)
      rethrow(err)
    else
      sprint(showerror, err) |> bad_request
    end
  end

  proc_result = try
    process(method, resource, req_data)
  catch err
    if isa(err, HttpError)
      rethrow(err)
    else
      sprint(showerror, err) |> internal_server_error
    end
  end

  try
    response = serialize(method, resource, proc_result)

    if !isa(response, Response)
      internal_server_error(
        "Serialization did not produce an `HTTP.Response`; check method implementation"
      )
    end

    return response
  catch err
    if isa(err, HttpError)
      rethrow(err)
    else
      sprint(showerror, err) |> internal_server_error
    end
  end
end

"""
    handler(::HttpMethod, ::Resource)

Produce a function which handles an `HTTP.Request` for a resource, given a method.
"""
handler(method::HttpMethod, resource::Resource) = (req -> handle(method, resource, req))

end # module Api


text:
using Plots
using Random
Random.seed!(1)

include_ddp()

# Model
include_model("double_integrator")

p = 2 * model.n + 2 * model.m + 2 * model.n

function f(model::DoubleIntegratorContinuous, x, u, w)
    [x[2]; u[1]]
end

function fd(model::DoubleIntegratorContinuous{Midpoint, FixedTime}, x, u, w, h, t)
	if t == 1
		θ = u[1 .+ (1:p)]
    	return [x + h * f(model, x + 0.5 * h * f(model, x, view(u, 1:1), w), view(u, 1:1), w); θ]
	else
		θ = x[2 .+ (1:p)]
		return [view(x, 1:2) + h * f(model, view(x, 1:2) + 0.5 * h * f(model, view(x, 1:2), u, w), u, w); θ]
	end
end

# Time
T = 11
h = 0.1

model = DoubleIntegratorContinuous{Midpoint, FixedTime}(2, 1, 0)
n = [t == 1 ? model.n : model.n + p for t = 1:T]
m = [t == 1 ? model.m + p : model.m for t = 1:T]

# Initial conditions, controls, disturbances
x1 = [0.0; 0.0]
xT = [1.0; 0.0]

Q10 = ones(model.n)
q10 = randn(model.n)
R10 = ones(model.m)
r10 = randn(model.m)
QT0 = ones(model.n)
qT0 = randn(model.n)
θ0 = [Q10; q10; R10; r10; QT0; qT0]
# θ0 = [1.0]
ū = [t == 1 ? [0.01 * randn(1); θ0] : 0.01 * randn(1) for t = 1:T-1]
w = [zeros(model.d) for t = 1:T-1]

# Rollout
x̄ = rollout(model, x1, ū, w, h, T)

# Objective

obj = StageCosts([NonlinearCost() for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		t == 1 ? (θ = u[1 .+ (1:p)]) : (θ = x[2 .+ (1:p)])

		Q = Diagonal(θ[1:model.n])
		q = θ[model.n .+ (1:model.n)]
		R = Diagonal(θ[2 * model.n .+ (1:model.m)])
		r = θ[2 * model.n + model.m .+ (1:model.m)]

        return 1.0e-5 * u' * u #+ view(x, 1:model.n)' * Q * view(x, 1:model.n) + q' * view(x, 1:model.n) + view(u, 1:model.m)' * R * view(u, 1:model.m) + r' * view(u, 1:model.m)
    elseif t == T
		θ = x[2 .+ (1:p)]

		Q = Diagonal(θ[2 * model.n + 2 * model.m .+ (1:model.n)])
		q = θ[2 * model.n + 2 * model.m + model.n .+ (1:model.n)]

        return view(x, 1:model.n)' * Q * view(x, 1:model.n) + q' * view(x, 1:model.n)
    else
        return 0.0
    end
end

g(obj, x̄[1], ū[1], 1)
g(obj, x̄[2], ū[2], 2)
g(obj, x̄[T], nothing, T)
objective(obj, x̄, ū)

# Constraints
p_con = [t == 1 ? 0 : (t == T ? model.n : 0) for t = 1:T]
info_1 = Dict()#):x1l => [0.0; -1.0], :x1u => [0.0; 1.0], :inequality => (1:2 * model.n))
info_M = Dict()#:xM => xM)
info_t = Dict()
info_T = Dict(:xT => xT)
con_set = [StageConstraint(p_con[t], t < T ? (t == 1 ? info_1 : info_t) : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	# if t == 1
	# 	c[1:2] .= view(u, 2:3) - cons.con[t].info[:x1u]
	# 	c[3:4] = cons.con[t].info[:x1l] - view(u, 2:3)
	# end
	# if t == Tm
	# 	c[1:2] .= view(x, 1:2) - cons.con[t].info[:xM]
	# end
	if t == T
		c[1:2] .= view(x, 1:2) - cons.con[t].info[:xT]
	end
end

prob = problem_data(model, obj, con_set, copy(x̄), copy(ū), w, h, T, n = n, m = m)
prob.m_data.obj
# objective_derivatives!(prob.m_data.obj, prob.m_data)

t = T
gx(z) = g(obj, z, t == T ? nothing : ū[t], t)
gu(z) = g(obj, x̄[t], z, t)
gz(z) = g(obj, z[1:n[t]], z[n[t] .+ (1:m[t])], t)

data = prob.m_data
ForwardDiff.gradient!(data.obj_deriv.gx[t], gx, x̄[t])
ForwardDiff.gradient!(data.obj_deriv.gu[t], gu, ū[t])
# ForwardDiff.hessian!(data.obj_deriv.gxx[t], gx, x̄[t])
using Symbolics
@variables x_sym[1:n[t]]

gx_sym = g(obj, x_sym, nothing, t)
gx_sym = simplify.(gx_sym)
dd_gx_sym = Symbolics.hessian(gx_sym, x_sym, simplify = true)

dd_gx_expr = Symbolics.build_function(dd_gx_sym, x_sym)[1]
dd_gx = eval(dd_gx_expr)
show(dd_gx(x̄[t]))

ForwardDiff.hessian!(data.obj_deriv.guu[t], gu, ū[t])
data.obj_deriv.gux[t] .= ForwardDiff.hessian(gz,
	[x̄[t]; ū[t]])[n[t] .+ (1:m[t]), 1:n[t]]

# Solve
@time constrained_ddp_solve!(prob,
    max_iter = 1000, max_al_iter = 10,
	ρ_init = 1.0, ρ_scale = 10.0)

# # Solve
# @time ddp_solve!(prob,
#     max_iter = 100, verbose = true)

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)

# x[1] = u[1][2:3]
# x̄[1] = ū[1][2:3]
# norm(x[1] - x[T][1:2])

# Visualize
using Plots
plot(hcat([xt[1:2] for xt in x]...)[:, 1:T]', color = :black, label = "")
# plot(hcat(u..., u[end])', linetype = :steppost)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
# An illustration of invariance in Julia Sets.

abstract Animal
type Cat <: Animal name end
type Dog <: Animal name end

dogs = Set([Dog("rex"), Dog("spot")])
@assert typeof(dogs) == Set{Dog}      # Julia infers this
@assert isa(dogs, Set{Dog})
@assert !isa(dogs, Set{Animal})       # Surprised?

@assert (try
  local animals::Set{Animal} = dogs   # Not allowed!
catch
  1
end) == 1


text:
using Scapin.Grid

@testset "Grid module" begin
    @testset "Connectivity" begin
        @testset "Connectivity, 2d" begin
            """
            1┬────────┬────────┬────────┐
             │        │        │        │
             │ (1, 4) │ (2, 4) │ (3, 4) │
             │        │        │        │
            4┼────────┼────────┼────────│
             │        │        │        │
             │ (1, 3) │ (2, 3) │ (3, 3) │
             │        │        │        │
            3┼────────┼────────┼────────│
             │        │        │        │
             │ (1, 2) │ (2, 2) │ (3, 2) │
             │        │        │        │
            2┼────────┼────────┼────────│
             │        │        │        │
             │ (1, 1) │ (2, 1) │ (3, 1) │
             │        │        │        │
            1┼────────┼────────┼────────┤
             1        2        3        1
            """
            N = (3, 4)
            d = size(N, 1)
            dims = Tuple(1:n for n in N)
            linear = LinearIndices(dims)
            @test cell_vertices(CartesianIndex(1, 1), N) ==
                  map(CartesianIndex, [(1, 1) (1, 2); (2, 1) (2, 2)])
            @test cell_vertices(CartesianIndex(2, 1), N) ==
                  map(CartesianIndex, [(2, 1) (2, 2); (3, 1) (3, 2)])
            @test cell_vertices(CartesianIndex(3, 1), N) ==
                  map(CartesianIndex, [(3, 1) (3, 2); (1, 1) (1, 2)])
            @test cell_vertices(CartesianIndex(1, 2), N) ==
                  map(CartesianIndex, [(1, 2) (1, 3); (2, 2) (2, 3)])
            @test cell_vertices(CartesianIndex(2, 2), N) ==
                  map(CartesianIndex, [(2, 2) (2, 3); (3, 2) (3, 3)])
            @test cell_vertices(CartesianIndex(3, 2), N) ==
                  map(CartesianIndex, [(3, 2) (3, 3); (1, 2) (1, 3)])
            @test cell_vertices(CartesianIndex(1, 3), N) ==
                  map(CartesianIndex, [(1, 3) (1, 4); (2, 3) (2, 4)])
            @test cell_vertices(CartesianIndex(2, 3), N) ==
                  map(CartesianIndex, [(2, 3) (2, 4); (3, 3) (3, 4)])
            @test cell_vertices(CartesianIndex(3, 3), N) ==
                  map(CartesianIndex, [(3, 3) (3, 4); (1, 3) (1, 4)])
            @test cell_vertices(CartesianIndex(1, 4), N) ==
                  map(CartesianIndex, [(1, 4) (1, 1); (2, 4) (2, 1)])
            @test cell_vertices(CartesianIndex(2, 4), N) ==
                  map(CartesianIndex, [(2, 4) (2, 1); (3, 4) (3, 1)])
            @test cell_vertices(CartesianIndex(3, 4), N) ==
                  map(CartesianIndex, [(3, 4) (3, 1); (1, 4) (1, 1)])
        end

        @testset "Connectivity, 3d" begin
            """
            1───────────2───────────3───────────4
            │           │           │           │
            │ (1, 4, 1) │ (2, 4, 1) │ (3, 4, 1) │
            │           │           │           │
            10─────────11───────────12─────────10
            │           │           │           │
            │ (1, 3, 1) │ (2, 3, 1) │ (3, 3, 1) │
            │           │           │           │
            7───────────8───────────9───────────7
            │           │           │           │
            │ (1, 2, 1) │ (2, 2, 1) │ (3, 2, 1) │
            │           │           │           │
            4───────────5───────────6───────────4
            │           │           │           │
            │ (1, 1, 1) │ (2, 1, 1) │ (3, 1, 1) │
            │           │           │           │
            1───────────2───────────3───────────1

            13─────────14───────────15─────────13
            │           │           │           │
            │ (1, 4, 2) │ (2, 4, 2) │ (3, 4, 2) │
            │           │           │           │
            22─────────23───────────24─────────22
            │           │           │           │
            │ (1, 3, 2) │ (2, 3, 2) │ (3, 3, 2) │
            │           │           │           │
            19─────────20───────────21─────────19
            │           │           │           │
            │ (1, 2, 2) │ (2, 2, 2) │ (3, 2, 2) │
            │           │           │           │
            16─────────17───────────18─────────16
            │           │           │           │
            │ (1, 1, 2) │ (2, 1, 2) │ (3, 1, 2) │
            │           │           │           │
            13─────────14───────────15─────────13

            25─────────26───────────27─────────25
            │           │           │           │
            │ (1, 4, 3) │ (2, 4, 3) │ (3, 4, 3) │
            │           │           │           │
            34─────────35───────────36─────────34
            │           │           │           │
            │ (1, 3, 3) │ (2, 3, 3) │ (3, 3, 3) │
            │           │           │           │
            31─────────32───────────33─────────31
            │           │           │           │
            │ (1, 2, 3) │ (2, 2, 3) │ (3, 2, 3) │
            │           │           │           │
            28─────────29───────────30─────────28
            │           │           │           │
            │ (1, 1, 3) │ (2, 1, 3) │ (3, 1, 3) │
            │           │           │           │
            25─────────26───────────27─────────25

            37─────────38───────────39─────────37
            │           │           │           │
            │ (1, 4, 4) │ (2, 4, 4) │ (3, 4, 4) │
            │           │           │           │
            46─────────47───────────48─────────46
            │           │           │           │
            │ (1, 3, 4) │ (2, 3, 4) │ (3, 3, 4) │
            │           │           │           │
            43─────────44───────────45─────────43
            │           │           │           │
            │ (1, 2, 4) │ (2, 2, 4) │ (3, 2, 4) │
            │           │           │           │
            40─────────41───────────42─────────40
            │           │           │           │
            │ (1, 1, 4) │ (2, 1, 4) │ (3, 1, 4) │
            │           │           │           │
            37─────────38───────────39─────────37

            49─────────50───────────51─────────49
            │           │           │           │
            │ (1, 4, 5) │ (2, 4, 5) │ (3, 4, 5) │
            │           │           │           │
            58─────────59───────────60─────────58
            │           │           │           │
            │ (1, 3, 5) │ (2, 3, 5) │ (3, 3, 5) │
            │           │           │           │
            55─────────56───────────57─────────55
            │           │           │           │
            │ (1, 2, 5) │ (2, 2, 5) │ (3, 2, 5) │
            │           │           │           │
            52─────────53───────────54─────────52
            │           │           │           │
            │ (1, 1, 5) │ (2, 1, 5) │ (3, 1, 5) │
            │           │           │           │
            49─────────50───────────51─────────49

            1───────────2───────────3───────────4
            │           │           │           │
            │           │           │           │
            │           │           │           │
            10─────────11───────────12─────────10
            │           │           │           │
            │           │           │           │
            │           │           │           │
            7───────────8───────────9───────────7
            │           │           │           │
            │           │           │           │
            │           │           │           │
            4───────────5───────────6───────────4
            │           │           │           │
            │           │           │           │
            │           │           │           │
            1───────────2───────────3───────────1

            """
            N = (3, 4, 5)
            d = size(N, 1)
            dims = Tuple(1:n for n in N)
            linear = LinearIndices(dims)
            transform(a) = reshape(map(CartesianIndex, a), 2, 2, 2)
            do_test(cell, vertices) =
                cell_vertices(CartesianIndex(cell), N) == transform(vertices)

            @test do_test(
                (1, 1, 1),
                [
                    (1, 1, 1),
                    (2, 1, 1),
                    (1, 2, 1),
                    (2, 2, 1),
                    (1, 1, 2),
                    (2, 1, 2),
                    (1, 2, 2),
                    (2, 2, 2),
                ],
            )
            @test do_test(
                (2, 1, 1),
                [
                    (2, 1, 1),
                    (3, 1, 1),
                    (2, 2, 1),
                    (3, 2, 1),
                    (2, 1, 2),
                    (3, 1, 2),
                    (2, 2, 2),
                    (3, 2, 2),
                ],
            )
            @test do_test(
                (3, 1, 1),
                [
                    (3, 1, 1),
                    (1, 1, 1),
                    (3, 2, 1),
                    (1, 2, 1),
                    (3, 1, 2),
                    (1, 1, 2),
                    (3, 2, 2),
                    (1, 2, 2),
                ],
            )
            @test do_test(
                (1, 2, 1),
                [
                    (1, 2, 1),
                    (2, 2, 1),
                    (1, 3, 1),
                    (2, 3, 1),
                    (1, 2, 2),
                    (2, 2, 2),
                    (1, 3, 2),
                    (2, 3, 2),
                ],
            )
            @test do_test(
                (2, 2, 1),
                [
                    (2, 2, 1),
                    (3, 2, 1),
                    (2, 3, 1),
                    (3, 3, 1),
                    (2, 2, 2),
                    (3, 2, 2),
                    (2, 3, 2),
                    (3, 3, 2),
                ],
            )
            @test do_test(
                (3, 2, 1),
                [
                    (3, 2, 1),
                    (1, 2, 1),
                    (3, 3, 1),
                    (1, 3, 1),
                    (3, 2, 2),
                    (1, 2, 2),
                    (3, 3, 2),
                    (1, 3, 2),
                ],
            )
            @test do_test(
                (1, 3, 1),
                [
                    (1, 3, 1),
                    (2, 3, 1),
                    (1, 4, 1),
                    (2, 4, 1),
                    (1, 3, 2),
                    (2, 3, 2),
                    (1, 4, 2),
                    (2, 4, 2),
                ],
            )
            @test do_test(
                (2, 3, 1),
                [
                    (2, 3, 1),
                    (3, 3, 1),
                    (2, 4, 1),
                    (3, 4, 1),
                    (2, 3, 2),
                    (3, 3, 2),
                    (2, 4, 2),
                    (3, 4, 2),
                ],
            )
            @test do_test(
                (3, 3, 1),
                [
                    (3, 3, 1),
                    (1, 3, 1),
                    (3, 4, 1),
                    (1, 4, 1),
                    (3, 3, 2),
                    (1, 3, 2),
                    (3, 4, 2),
                    (1, 4, 2),
                ],
            )
            @test do_test(
                (1, 4, 1),
                [
                    (1, 4, 1),
                    (2, 4, 1),
                    (1, 1, 1),
                    (2, 1, 1),
                    (1, 4, 2),
                    (2, 4, 2),
                    (1, 1, 2),
                    (2, 1, 2),
                ],
            )
            @test do_test(
                (2, 4, 1),
                [
                    (2, 4, 1),
                    (3, 4, 1),
                    (2, 1, 1),
                    (3, 1, 1),
                    (2, 4, 2),
                    (3, 4, 2),
                    (2, 1, 2),
                    (3, 1, 2),
                ],
            )
            @test do_test(
                (3, 4, 1),
                [
                    (3, 4, 1),
                    (1, 4, 1),
                    (3, 1, 1),
                    (1, 1, 1),
                    (3, 4, 2),
                    (1, 4, 2),
                    (3, 1, 2),
                    (1, 1, 2),
                ],
            )

            @test do_test(
                (1, 1, 2),
                [
                    (1, 1, 2),
                    (2, 1, 2),
                    (1, 2, 2),
                    (2, 2, 2),
                    (1, 1, 3),
                    (2, 1, 3),
                    (1, 2, 3),
                    (2, 2, 3),
                ],
            )
            @test do_test(
                (2, 1, 2),
                [
                    (2, 1, 2),
                    (3, 1, 2),
                    (2, 2, 2),
                    (3, 2, 2),
                    (2, 1, 3),
                    (3, 1, 3),
                    (2, 2, 3),
                    (3, 2, 3),
                ],
            )
            @test do_test(
                (3, 1, 2),
                [
                    (3, 1, 2),
                    (1, 1, 2),
                    (3, 2, 2),
                    (1, 2, 2),
                    (3, 1, 3),
                    (1, 1, 3),
                    (3, 2, 3),
                    (1, 2, 3),
                ],
            )
            @test do_test(
                (1, 2, 2),
                [
                    (1, 2, 2),
                    (2, 2, 2),
                    (1, 3, 2),
                    (2, 3, 2),
                    (1, 2, 3),
                    (2, 2, 3),
                    (1, 3, 3),
                    (2, 3, 3),
                ],
            )
            @test do_test(
                (2, 2, 2),
                [
                    (2, 2, 2),
                    (3, 2, 2),
                    (2, 3, 2),
                    (3, 3, 2),
                    (2, 2, 3),
                    (3, 2, 3),
                    (2, 3, 3),
                    (3, 3, 3),
                ],
            )
            @test do_test(
                (3, 2, 2),
                [
                    (3, 2, 2),
                    (1, 2, 2),
                    (3, 3, 2),
                    (1, 3, 2),
                    (3, 2, 3),
                    (1, 2, 3),
                    (3, 3, 3),
                    (1, 3, 3),
                ],
            )
            @test do_test(
                (1, 3, 2),
                [
                    (1, 3, 2),
                    (2, 3, 2),
                    (1, 4, 2),
                    (2, 4, 2),
                    (1, 3, 3),
                    (2, 3, 3),
                    (1, 4, 3),
                    (2, 4, 3),
                ],
            )
            @test do_test(
                (2, 3, 2),
                [
                    (2, 3, 2),
                    (3, 3, 2),
                    (2, 4, 2),
                    (3, 4, 2),
                    (2, 3, 3),
                    (3, 3, 3),
                    (2, 4, 3),
                    (3, 4, 3),
                ],
            )
            @test do_test(
                (3, 3, 2),
                [
                    (3, 3, 2),
                    (1, 3, 2),
                    (3, 4, 2),
                    (1, 4, 2),
                    (3, 3, 3),
                    (1, 3, 3),
                    (3, 4, 3),
                    (1, 4, 3),
                ],
            )
            @test do_test(
                (1, 4, 2),
                [
                    (1, 4, 2),
                    (2, 4, 2),
                    (1, 1, 2),
                    (2, 1, 2),
                    (1, 4, 3),
                    (2, 4, 3),
                    (1, 1, 3),
                    (2, 1, 3),
                ],
            )
            @test do_test(
                (2, 4, 2),
                [
                    (2, 4, 2),
                    (3, 4, 2),
                    (2, 1, 2),
                    (3, 1, 2),
                    (2, 4, 3),
                    (3, 4, 3),
                    (2, 1, 3),
                    (3, 1, 3),
                ],
            )
            @test do_test(
                (3, 4, 2),
                [
                    (3, 4, 2),
                    (1, 4, 2),
                    (3, 1, 2),
                    (1, 1, 2),
                    (3, 4, 3),
                    (1, 4, 3),
                    (3, 1, 3),
                    (1, 1, 3),
                ],
            )

            @test do_test(
                (1, 1, 3),
                [
                    (1, 1, 3),
                    (2, 1, 3),
                    (1, 2, 3),
                    (2, 2, 3),
                    (1, 1, 4),
                    (2, 1, 4),
                    (1, 2, 4),
                    (2, 2, 4),
                ],
            )
            @test do_test(
                (2, 1, 3),
                [
                    (2, 1, 3),
                    (3, 1, 3),
                    (2, 2, 3),
                    (3, 2, 3),
                    (2, 1, 4),
                    (3, 1, 4),
                    (2, 2, 4),
                    (3, 2, 4),
                ],
            )
            @test do_test(
                (3, 1, 3),
                [
                    (3, 1, 3),
                    (1, 1, 3),
                    (3, 2, 3),
                    (1, 2, 3),
                    (3, 1, 4),
                    (1, 1, 4),
                    (3, 2, 4),
                    (1, 2, 4),
                ],
            )
            @test do_test(
                (1, 2, 3),
                [
                    (1, 2, 3),
                    (2, 2, 3),
                    (1, 3, 3),
                    (2, 3, 3),
                    (1, 2, 4),
                    (2, 2, 4),
                    (1, 3, 4),
                    (2, 3, 4),
                ],
            )
            @test do_test(
                (2, 2, 3),
                [
                    (2, 2, 3),
                    (3, 2, 3),
                    (2, 3, 3),
                    (3, 3, 3),
                    (2, 2, 4),
                    (3, 2, 4),
                    (2, 3, 4),
                    (3, 3, 4),
                ],
            )
            @test do_test(
                (3, 2, 3),
                [
                    (3, 2, 3),
                    (1, 2, 3),
                    (3, 3, 3),
                    (1, 3, 3),
                    (3, 2, 4),
                    (1, 2, 4),
                    (3, 3, 4),
                    (1, 3, 4),
                ],
            )
            @test do_test(
                (1, 3, 3),
                [
                    (1, 3, 3),
                    (2, 3, 3),
                    (1, 4, 3),
                    (2, 4, 3),
                    (1, 3, 4),
                    (2, 3, 4),
                    (1, 4, 4),
                    (2, 4, 4),
                ],
            )
            @test do_test(
                (2, 3, 3),
                [
                    (2, 3, 3),
                    (3, 3, 3),
                    (2, 4, 3),
                    (3, 4, 3),
                    (2, 3, 4),
                    (3, 3, 4),
                    (2, 4, 4),
                    (3, 4, 4),
                ],
            )
            @test do_test(
                (3, 3, 3),
                [
                    (3, 3, 3),
                    (1, 3, 3),
                    (3, 4, 3),
                    (1, 4, 3),
                    (3, 3, 4),
                    (1, 3, 4),
                    (3, 4, 4),
                    (1, 4, 4),
                ],
            )
            @test do_test(
                (1, 4, 3),
                [
                    (1, 4, 3),
                    (2, 4, 3),
                    (1, 1, 3),
                    (2, 1, 3),
                    (1, 4, 4),
                    (2, 4, 4),
                    (1, 1, 4),
                    (2, 1, 4),
                ],
            )
            @test do_test(
                (2, 4, 3),
                [
                    (2, 4, 3),
                    (3, 4, 3),
                    (2, 1, 3),
                    (3, 1, 3),
                    (2, 4, 4),
                    (3, 4, 4),
                    (2, 1, 4),
                    (3, 1, 4),
                ],
            )
            @test do_test(
                (3, 4, 3),
                [
                    (3, 4, 3),
                    (1, 4, 3),
                    (3, 1, 3),
                    (1, 1, 3),
                    (3, 4, 4),
                    (1, 4, 4),
                    (3, 1, 4),
                    (1, 1, 4),
                ],
            )
            @test do_test(
                (1, 1, 4),
                [
                    (1, 1, 4),
                    (2, 1, 4),
                    (1, 2, 4),
                    (2, 2, 4),
                    (1, 1, 5),
                    (2, 1, 5),
                    (1, 2, 5),
                    (2, 2, 5),
                ],
            )
            @test do_test(
                (2, 1, 4),
                [
                    (2, 1, 4),
                    (3, 1, 4),
                    (2, 2, 4),
                    (3, 2, 4),
                    (2, 1, 5),
                    (3, 1, 5),
                    (2, 2, 5),
                    (3, 2, 5),
                ],
            )
            @test do_test(
                (3, 1, 4),
                [
                    (3, 1, 4),
                    (1, 1, 4),
                    (3, 2, 4),
                    (1, 2, 4),
                    (3, 1, 5),
                    (1, 1, 5),
                    (3, 2, 5),
                    (1, 2, 5),
                ],
            )
            @test do_test(
                (1, 2, 4),
                [
                    (1, 2, 4),
                    (2, 2, 4),
                    (1, 3, 4),
                    (2, 3, 4),
                    (1, 2, 5),
                    (2, 2, 5),
                    (1, 3, 5),
                    (2, 3, 5),
                ],
            )
            @test do_test(
                (2, 2, 4),
                [
                    (2, 2, 4),
                    (3, 2, 4),
                    (2, 3, 4),
                    (3, 3, 4),
                    (2, 2, 5),
                    (3, 2, 5),
                    (2, 3, 5),
                    (3, 3, 5),
                ],
            )
            @test do_test(
                (3, 2, 4),
                [
                    (3, 2, 4),
                    (1, 2, 4),
                    (3, 3, 4),
                    (1, 3, 4),
                    (3, 2, 5),
                    (1, 2, 5),
                    (3, 3, 5),
                    (1, 3, 5),
                ],
            )
            @test do_test(
                (1, 3, 4),
                [
                    (1, 3, 4),
                    (2, 3, 4),
                    (1, 4, 4),
                    (2, 4, 4),
                    (1, 3, 5),
                    (2, 3, 5),
                    (1, 4, 5),
                    (2, 4, 5),
                ],
            )
            @test do_test(
                (2, 3, 4),
                [
                    (2, 3, 4),
                    (3, 3, 4),
                    (2, 4, 4),
                    (3, 4, 4),
                    (2, 3, 5),
                    (3, 3, 5),
                    (2, 4, 5),
                    (3, 4, 5),
                ],
            )
            @test do_test(
                (3, 3, 4),
                [
                    (3, 3, 4),
                    (1, 3, 4),
                    (3, 4, 4),
                    (1, 4, 4),
                    (3, 3, 5),
                    (1, 3, 5),
                    (3, 4, 5),
                    (1, 4, 5),
                ],
            )
            @test do_test(
                (1, 4, 4),
                [
                    (1, 4, 4),
                    (2, 4, 4),
                    (1, 1, 4),
                    (2, 1, 4),
                    (1, 4, 5),
                    (2, 4, 5),
                    (1, 1, 5),
                    (2, 1, 5),
                ],
            )
            @test do_test(
                (2, 4, 4),
                [
                    (2, 4, 4),
                    (3, 4, 4),
                    (2, 1, 4),
                    (3, 1, 4),
                    (2, 4, 5),
                    (3, 4, 5),
                    (2, 1, 5),
                    (3, 1, 5),
                ],
            )
            @test do_test(
                (3, 4, 4),
                [
                    (3, 4, 4),
                    (1, 4, 4),
                    (3, 1, 4),
                    (1, 1, 4),
                    (3, 4, 5),
                    (1, 4, 5),
                    (3, 1, 5),
                    (1, 1, 5),
                ],
            )
            @test do_test(
                (1, 1, 5),
                [
                    (1, 1, 5),
                    (2, 1, 5),
                    (1, 2, 5),
                    (2, 2, 5),
                    (1, 1, 1),
                    (2, 1, 1),
                    (1, 2, 1),
                    (2, 2, 1),
                ],
            )
            @test do_test(
                (2, 1, 5),
                [
                    (2, 1, 5),
                    (3, 1, 5),
                    (2, 2, 5),
                    (3, 2, 5),
                    (2, 1, 1),
                    (3, 1, 1),
                    (2, 2, 1),
                    (3, 2, 1),
                ],
            )
            @test do_test(
                (3, 1, 5),
                [
                    (3, 1, 5),
                    (1, 1, 5),
                    (3, 2, 5),
                    (1, 2, 5),
                    (3, 1, 1),
                    (1, 1, 1),
                    (3, 2, 1),
                    (1, 2, 1),
                ],
            )
            @test do_test(
                (1, 2, 5),
                [
                    (1, 2, 5),
                    (2, 2, 5),
                    (1, 3, 5),
                    (2, 3, 5),
                    (1, 2, 1),
                    (2, 2, 1),
                    (1, 3, 1),
                    (2, 3, 1),
                ],
            )
            @test do_test(
                (2, 2, 5),
                [
                    (2, 2, 5),
                    (3, 2, 5),
                    (2, 3, 5),
                    (3, 3, 5),
                    (2, 2, 1),
                    (3, 2, 1),
                    (2, 3, 1),
                    (3, 3, 1),
                ],
            )
            @test do_test(
                (3, 2, 5),
                [
                    (3, 2, 5),
                    (1, 2, 5),
                    (3, 3, 5),
                    (1, 3, 5),
                    (3, 2, 1),
                    (1, 2, 1),
                    (3, 3, 1),
                    (1, 3, 1),
                ],
            )
            @test do_test(
                (1, 3, 5),
                [
                    (1, 3, 5),
                    (2, 3, 5),
                    (1, 4, 5),
                    (2, 4, 5),
                    (1, 3, 1),
                    (2, 3, 1),
                    (1, 4, 1),
                    (2, 4, 1),
                ],
            )
            @test do_test(
                (2, 3, 5),
                [
                    (2, 3, 5),
                    (3, 3, 5),
                    (2, 4, 5),
                    (3, 4, 5),
                    (2, 3, 1),
                    (3, 3, 1),
                    (2, 4, 1),
                    (3, 4, 1),
                ],
            )
            @test do_test(
                (3, 3, 5),
                [
                    (3, 3, 5),
                    (1, 3, 5),
                    (3, 4, 5),
                    (1, 4, 5),
                    (3, 3, 1),
                    (1, 3, 1),
                    (3, 4, 1),
                    (1, 4, 1),
                ],
            )
            @test do_test(
                (1, 4, 5),
                [
                    (1, 4, 5),
                    (2, 4, 5),
                    (1, 1, 5),
                    (2, 1, 5),
                    (1, 4, 1),
                    (2, 4, 1),
                    (1, 1, 1),
                    (2, 1, 1),
                ],
            )
            @test do_test(
                (2, 4, 5),
                [
                    (2, 4, 5),
                    (3, 4, 5),
                    (2, 1, 5),
                    (3, 1, 5),
                    (2, 4, 1),
                    (3, 4, 1),
                    (2, 1, 1),
                    (3, 1, 1),
                ],
            )
            @test do_test(
                (3, 4, 5),
                [
                    (3, 4, 5),
                    (1, 4, 5),
                    (3, 1, 5),
                    (1, 1, 5),
                    (3, 4, 1),
                    (1, 4, 1),
                    (3, 1, 1),
                    (1, 1, 1),
                ],
            )
        end
    end
end


text:
function [X,J,dXdx]=JacobiCart(ksi,eta,F,Grid)
n1=size(ksi,1);
n2=size(eta,2);
X=zeros(n1,n2,3);
X(:,:,1)=0.25*(F.P(1,1).*(1-ksi)*(1-eta)+...
  F.P(1,2).*(1+ksi)*(1-eta)+...
  F.P(1,3).*(1+ksi)*(1+eta)+...
  F.P(1,4).*(1-ksi)*(1+eta));
X(:,:,2)=0.25*(F.P(2,1).*(1-ksi)*(1-eta)+...
  F.P(2,2).*(1+ksi)*(1-eta)+...
  F.P(2,3).*(1+ksi)*(1+eta)+...
  F.P(2,4).*(1-ksi)*(1+eta));
X(:,:,3)=0;

dXdx=zeros(n1,n2,2,2);

dXdx(:,:,1,1)=0.25*(-F.P(1,1)+F.P(1,2)+...
  F.P(1,3)-F.P(1,4)+...
  ( F.P(1,1)-F.P(1,2)+...
  F.P(1,3)-F.P(1,4))*(ones(n1,1)*eta));
dXdx(:,:,1,2)=0.25*(-F.P(1,1)-F.P(1,2)+...
  F.P(1,3)+F.P(1,4)+...
  ( F.P(1,1)-F.P(1,2)+...
  F.P(1,3)-F.P(1,4))*(ksi*ones(1,n2)));
dXdx(:,:,2,1)=0.25*(-F.P(2,1)+F.P(2,2)+...
  F.P(2,3)-F.P(2,4)+...
  ( F.P(2,1)-F.P(2,2)+...
  F.P(2,3)-F.P(2,4))*(ones(n1,1)*eta));
dXdx(:,:,2,2)=0.25*(-F.P(2,1)-F.P(2,2)+...
  F.P(2,3)+F.P(2,4)+...
  ( F.P(2,1)-F.P(2,2)+...
  F.P(2,3)-F.P(2,4))*(ksi*ones(1,n2)));
J=dXdx(:,:,1,1).*dXdx(:,:,2,2)-dXdx(:,:,1,2).*dXdx(:,:,2,1);
end


text:
using DynamicalSystems
using Plots
plotly()

# Simulate a sine wave using DynamicalSystems, compare accuracy and timing to builtin sine and use this to construct Process type containing (only) relevant metadata

X0 = [0.0]
T = 2000.0
μ = 1.0 # Hz

function dynamicalSine(X, μ, t)
    dX = μ.*cos(μ.*t) # dx/dt
    return SVector{1}(dX)
end
function dynamicalSineJ(X, μ, t)
    J = @SMatrix [-μ^2*sin(μ.*t)]
end

# (x, μ, t) -> sin(μ.*t)
ds = ContinuousDynamicalSystem(dynamicalSine, X0, μ, dynamicalSineJ) # Easy system defined anonymously
data = trajectory(ds, T; dt=0.01)

plot(data[:, 1], seriestype=:path)
plot!(title="Sine")


text:
#
# step.jl --
#
# Implement the STEP (Select The Easiest Point) method for global univariate
# optimization [1].
#
# [1] S. Swarzberg, G. Seront & H. Bersini, "S.T.E.P.: the easiest way to
#     optimize a function," in IEEE World Congress on Computational
#     Intelligence, Proceedings of the First IEEE Conference on Evolutionary
#     Computation,vol. 1, pp. 519-524 (1994).
#
# -----------------------------------------------------------------------------
#
# This file is part of OptimPackNextGen.jl which is licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2015-2018, Éric Thiébaut.
# <https://github.com/emmt/OptimPackNextGen.jl>.
#

# FIXME: use priority queues
# FIXME: fix estimation of precision (compared to Yorick version)
# FIXME: make sure no other minima exist when precision test is satisfied

module Step

# Use the same floating point type for scalars as in OptimPackNextGen.
import OptimPackNextGen.Float

using Printf

"""
# Cyclic singly linked list
"""
mutable struct Item{T}
    next::Item{T}
    data::T
    (::Type{Item{T}})(next::Item{T}, data::T) where {T} = new{T}(next, data)
    function (::Type{Item{T}})(data::T) where {T}
        newitem = new{T}()
        newitem.next = newitem
        newitem.data = data
        return newitem
    end
end
Item(data::T) where {T} = Item{T}(data)

"""
Append a new item after a given item of a cyclic singly linked list.
"""
function append!(item::Item{T}, data::T) where {T}
    newitem = Item{T}(item.next, data)
    item.next = newitem
    return newitem
end

mutable struct NodeData
    x::Float # position
    y::Float # function value
    q::Float # "quality" factor
end

# Default absolute and relative tolerances:
const TOL = (floatmin(Float), sqrt(eps(Float)))

@inline sqrtdifmin(lvl,val) = sqrt(val - lvl)
@inline sqrtdifmax(lvl,val) = sqrt(lvl - val)

for (func, cmp, incr, wgt) in ((:minimize, <, -, :sqrtdifmin),
                               (:maximize, >, +, :sqrtdifmax))
    @eval begin

        function $func(f::Function, a::Float, b::Float;
                       maxeval::Int=100000,
                       tol::NTuple{2,Float}=TOL,
                       alpha::Float=0.0, beta::Float=0.0,
                       verb::Bool=false,
                       printer::Function=default_printer, output::IO=stdout)

            maxeval ≥ 2 || error("parameter `maxeval` must be at least 2")
            tol[1] ≥ 0 || error("absolute tolerance `tol[1]` must be nonnegative")
            0 ≤ tol[2] ≤ 1 || error("relative tolerance `tol[2]` must be in [0,1]")
            alpha ≥ 0 || error("parameter `alpha` must be nonnegative")
            beta ≥ 0 || error("parameter `beta` must be nonnegative")

            if a > b
                (a, b) = (b, a)
            elseif a == b
                return (a, f(a), zero(Float), 1)
            end

            fa = f(a)
            fb = f(b)
            list = Item(NodeData(a, fa, 0))
            append!(list, NodeData(b, fb, 0))
            evaluations::Int = 2

            if $cmp(fa, fb)
                xbest = a
                fbest = fa
            else
                xbest = b
                fbest = fb
            end
            xtol::Float = (b - a)/2
            rehash::Bool = true
            verb && printer(output, evaluations, xbest, fbest, xtol)
            c = xbest
            while true
                if xtol ≤ hypot(tol[1], xbest*tol[2])
                    break
                end
                if evaluations ≥ maxeval
                    @warn "too many evaluations"
                    break
                end

                # Find where to split.
                split = list
                if rehash
                    qmin = Float(Inf)
                    t = hypot(alpha*fbest, beta)
                    level = $incr(fbest, t)
                    n2 = list
                    x2 = n2.data.x
                    w2 = $wgt(level, n2.data.y)
                    while true
                        n1 = n2
                        n2 = n2.next
                        x1 = x2
                        x2 = n2.data.x
                        x2 > x1 || break
                        w1 = w2
                        w2 = $wgt(level, n2.data.y)
                        q = (w1 + w2)/(x2 - x1)
                        n1.data.q = q
                        if q < qmin
                            qmin = q
                            split = n1
                        end
                    end
                    rehash = false
                else
                    qmin = split.data.q
                    n1 = list
                    x1 = n1.data.x
                    while true
                        n2 = n1.next
                        x2 = n2.data.x
                        x2 > x1 || break
                        q = n1.data.q
                        if q < qmin
                            qmin = q
                            split = n1
                        end
                        n1, x1 = n2, x2
                    end
                end

                # Split the chosen interval.
                evaluations += 1
                x0 = split.data.x
                x1 = split.next.data.x
                c = (x0 + x1)/2
                fc = f(c)
                if $cmp(fc, fbest)
                    xbest = c
                    fbest = fc
                    xtol = (x1 - x0)/2
                    rehash = true
                    verb && printer(output, evaluations, xbest, fbest, xtol)
                end
                if rehash
                    # All Q factors have to be recomputed.
                    q = zero(Float)
                else
                    # Compute the Q factors for the split interval and for the
                    # new one.
                    w1 = $wgt(level, split.data.y)
                    w2 = $wgt(level, fc)
                    w3 = $wgt(level, split.next.data.y)
                    e = (x1 - x0)/2
                    split.data.q = (w1 + w2)/e
                    q = (w2 + w3)/e
                end
                append!(split, NodeData(c, fc, q))
            end

            # Return best point found so far.
            (xbest, fbest, xtol, evaluations)
        end

        function $func(f::Function, a::Real, b::Real;
                       maxeval::Integer=10000,
                       tol::NTuple{2,Real}=TOL,
                       alpha::Real=0.0, beta::Real=0.0,
                       kwds...)
            $func(f, Float(a), Float(b);
                  maxeval = Int(maxeval),
                  tol = (Float(tol[1]), Float(tol[2])),
                  alpha = Float(alpha),
                  beta = Float(beta),
                  kwds...)
        end

    end
end

"""
# Find a global minimum or maximum

    Step.minimize(f, a, b) -> (xbest, fbest, xtol, n)
    Step.maximize(f, a, b) -> (xbest, fbest, xtol, n)

finds a global minimum (resp. maximum) of `f(x)` in the interval `[a,b]` and
returns its position `xbest`, the corresponding function value `fbest =
f(xbest)`, the uncertainty `xtol` and the number of function evaluations needed
to find it.

The algorithm is based on the STEP method described in:

> Swarzberg, S., Seront, G. & Bersini, H., "S.T.E.P.: the easiest way to
> optimize a function," in IEEE World Congress on Computational Intelligence,
> Proceedings of the First IEEE Conference on Evolutionary Computation, vol. 1,
> pp. 519-524 (1994).

The following optional keywords can be used:

* `printer` can be set with a user defined function to print iteration
  information, its signature is:

      printer(io::IO, iter::Integer, eval::Integer, rejects::Integer,
              f::Real, gnorm::Real, stp::Real)

  where `io` is the output stream, `iter` the iteration number (`iter = 0` for
  the starting point), `eval` is the number of calls to `fg!`, `rejects` is the
  number of times the computed direction was rejected, `f` and `gnorm` are the
  value of the function and norm of the gradient at the current point, `stp` is
  the length of the step to the current point.

* `output` specifies the output stream for printing information (`stdout` is
  used by default).

""" minimize

@doc @doc(minimize) maximize

function default_printer(eval::Int, xm::Float, fm::Float, prec::Float)
    default_printer(stdout, eval, xm, fm, prec)
end

function default_printer(io::IO, eval::Int, xm::Float, fm::Float, prec::Float)
    if eval < 3
        @printf(io, "# %s%s\n# %s%s\n",
                "EVALS              X                      F(X)        ",
                "       PREC",
                "------------------------------------------------------",
                "-------------")
    end
    @printf(io, "%7d  %23.15e  %23.15e  %10.2e\n", eval, xm, fm, prec)
end

# Simple parabola.  To be minimized over [-1,2].
testParabola(x) = x*x

# Brent's 5th function.  To be minimized over [-10,10].
testBrent5(x) = (x - sin(x))*exp(-x*x)

# Michalewicz's 1st function.  To be minimized over [-1,2].
testMichalewicz1(x) = x*sin(10.0*x)

# Michalewicz's 2nd function.  To be maximized over [0,pi].
function testMichalewicz2(x)
    s = 0.0
    a = sin(x)
    b = x*x/pi
    for i in 1:10
        s += a*(sin(b*i)^20)
    end
    return s
end

function runtests()
    println("\n# Simple parabola:")
    (xbest, fbest, xtol, n) = minimize(testParabola, -1, 2, verb=true,
                                       maxeval=100, alpha=0, beta=0)
    println("x = $xbest ± $xtol, f(x) = $fbest, n = $n")

    println("\n# Brent's 5th function:")
    (xbest, fbest, xtol, n) = minimize(testBrent5, -10, 10, verb=true,
                                       maxeval=1000)
    println("x = $xbest ± $xtol, f(x) = $fbest, n = $n")

    println("\n# Michalewicz's 1st function:")
    (xbest, fbest, xtol, n) = minimize(testMichalewicz1, -1, 2, verb=true,
                                       maxeval=1000)
    println("x = $xbest ± $xtol, f(x) = $fbest, n = $n")

    println("\n# Michalewicz's 2nd function:")
    (xbest, fbest, xtol, n) = maximize(testMichalewicz2, 0, pi, verb=true,
                                       maxeval=1000, tol=(1e-12,0))
    println("x = $xbest ± $xtol, f(x) = $fbest, n = $n")
end

end # module


text: issorted("abc") # true

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode