text:
const SVF = MOI.SingleVariable
const VVF = MOI.VectorOfVariables
const SAF{T} = MOI.ScalarAffineFunction{T}
const VAF{T} = MOI.VectorAffineFunction{T}

const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

const GT = MOI.GreaterThan
const LT = MOI.LessThan
const ET = MOI.EqualTo

const SCALAR_SETS = Union{
    MOI.GreaterThan{Float64},
    MOI.LessThan{Float64},
    MOI.EqualTo{Float64},
}

const VECTOR_SETS = Union{
    MOI.SecondOrderCone,
    MOI.RotatedSecondOrderCone,
}

# abstract type AsbtractBilevelOptimizer end
# struct SOS1Optimizer{O} <: AsbtractBilevelOptimizer
#     solver::O
#     # options
# end
# function SOS1Optimizer(solver::O) where O
#     return SOS1Optimizer{O}(solver)
# end

struct Complement#{M1 <: MOI.ModelLike, M2 <: MOI.ModelLike, F, S}
    is_vec
    # primal::M1
    func_w_cte#::F
    set_w_zero#::S
    # dual::M2
    variable#::VI
    # var_set#::S2
end

abstract type BilevelSolverMode{T} end

mutable struct SOS1Mode{T} <: BilevelSolverMode{T}
    epsilon::T
    function SOS1Mode()
        return new{Float64}(zero(Float64))
    end
end

mutable struct PositiveSOS1Mode{T} <: BilevelSolverMode{T}
    epsilon::T
    function PositiveSOS1Mode()
        return new{Float64}(zero(Float64))
    end
end

mutable struct ComplementMode{T} <: BilevelSolverMode{T}
    epsilon::T
    function ComplementMode()
        return new{Float64}(zero(Float64))
    end
end

mutable struct ComplementWithSlackMode{T} <: BilevelSolverMode{T}
    epsilon::T
    function ComplementWithSlackMode()
        return new{Float64}(zero(Float64))
    end
end

mutable struct ProductMode{T} <: BilevelSolverMode{T}
    epsilon::T
    function ProductMode()
        return new{Float64}(zero(Float64))
    end
    function ProductMode{T}(eps::T) where T
        return new{Float64}(eps)
    end
end

mutable struct ProductWithSlackMode{T} <: BilevelSolverMode{T}
    epsilon::T
    function ProductWithSlackMode()
        return new{Float64}(zero(Float64))
    end
    function ProductWithSlackMode{T}(eps::T) where T
        return new{Float64}(eps)
    end
end

@enum IndicatorSetting ZERO_ONE ZERO_ZERO ONE_ONE

mutable struct IndicatorMode{T} <: BilevelSolverMode{T}
    mode::IndicatorSetting
    function IndicatorMode()
        return new{Float64}(ONE_ONE)
    end
end

abstract type StrongDualityMode{T} <: BilevelSolverMode{T} end

mutable struct StrongDualityEqualityMode{T} <: StrongDualityMode{T}
    function StrongDualityEqualityMode()
        return new{Float64}()
    end
end
mutable struct StrongDualityInequalityMode{T} <: StrongDualityMode{T}
    epsilon::T
    function StrongDualityInequalityMode()
        return new{Float64}(zero(Float64))
    end
    function StrongDualityInequalityMode{T}(eps::T) where T
        return new{Float64}(eps)
    end
end

ignore_dual_objective(::BilevelSolverMode{T}) where T = true
ignore_dual_objective(::StrongDualityMode{T}) where T = false

function accept_vector_set(mode::BilevelSolverMode{T}, con::Complement) where T
    if con.is_vec
        error("Set $(typeof(con.set_w_zero)) is not accepted when solution method is $(typeof(mode))")
    end
    return nothing
end
accept_vector_set(::ProductMode{T}, ::Complement) where T = nothing
accept_vector_set(::ProductWithSlackMode{T}, ::Complement) where T = nothing

function get_canonical_complements(primal_model, primal_dual_map)
    map = primal_dual_map.primal_con_dual_var
    out = Complement[]
    for ci in keys(map)
        con = get_canonical_complement(primal_model, map, ci)
        push!(out, con)
    end
    return out
end
function get_canonical_complement(primal_model, map,
    ci::CI{F,S}) where {F, S<:VECTOR_SETS}
    T = Float64
    func = MOI.get(primal_model, MOI.ConstraintFunction(), ci)::F
    set = MOI.get(primal_model, MOI.ConstraintSet(), ci)::S
    dim = MOI.dimension(set)
    # vector sets have no constant
    # for i in 1:dim
    #     func.constant[i] = Dualization.set_dot(i, set, T) *
    #         Dualization.get_scalar_term(primal_model, i, ci)
    # end
    # todo - set dot on function
    con = Complement(true, func, set_with_zero(set), map[ci])
    return con
end
function get_canonical_complement(primal_model, map,
    ci::CI{F,S}) where {F, S<:SCALAR_SETS}
    T = Float64
    func = MOI.get(primal_model, MOI.ConstraintFunction(), ci)::F
    set = MOI.get(primal_model, MOI.ConstraintSet(), ci)::S
    constant = Dualization.set_dot(1, set, T) *
        Dualization.get_scalar_term(primal_model, 1, ci)
    if F == MOI.SingleVariable
        func = MOIU.operate(+, T, func, constant)
    else
        func.constant = constant
    end
    # todo - set dot on function
    con = Complement(false, func, set_with_zero(set), map[ci][1])
    return con
end

function set_with_zero(set::S) where {S<:SCALAR_SETS} where T
    return S(0.0)
end
function set_with_zero(set)
    return copy(set)
end

function build_bilevel(
    upper::MOI.ModelLike, lower::MOI.ModelLike,
    link::Dict{VI,VI}, upper_variables::Vector{VI},
    mode,
    upper_var_lower_ctr::Dict{VI,CI} = Dict{VI,CI}()
    )

    # Start with an empty problem
    moi_mode = MOIU.AUTOMATIC
    # m = MOIU.CachingOptimizer(MOIU.Model{Float64}(), moi_mode)
    m = MOIU.CachingOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()), moi_mode)

    # add the first level
    copy_names = false

    #=
        Initialize Lower DUAL level model
    =#
    # dualize the second level
    dual_problem = dualize(lower,
        dual_names = DualNames("dual_","dual_"),
        variable_parameters = upper_variables,
        ignore_objective = ignore_dual_objective(mode))
    lower_dual = dual_problem.dual_model
    lower_primal_dual_map = dual_problem.primal_dual_map

    #=
        Pass Upper level model
    =#

    # key are from src, value are from dest
    upper_idxmap = MOIU.default_copy_to(m, upper, copy_names)
    pass_names(m, upper, upper_idxmap)

    #=
        Pass Lower level model
    =#

    # cache and delete lower objective
    if !ignore_dual_objective(mode)
        # get primal obj
        tp_primal_obj = MOI.get(lower, MOI.ObjectiveFunctionType())
        @assert tp_primal_obj !== nothing
        lower_primal_obj = MOI.get(lower, MOI.ObjectiveFunction{tp_primal_obj}())
        # delete dual obj
        # MOI.set(lower, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0))
    end

    # initialize map to lower level model
    lower_idxmap = MOIU.IndexMap()
    for (upper_key, lower_val) in link
        lower_idxmap[lower_val] = upper_idxmap[upper_key]
    end

    # append the second level primal
    append_to(m, lower, lower_idxmap, copy_names, allow_single_bounds = true)
    pass_names(m, lower, lower_idxmap)

    #=
        Pass Dual of Lower level model
    =#

    # initialize map to lower level model
    if !ignore_dual_objective(mode)
        # get dual obj
        tp_dual_obj = MOI.get(lower_dual, MOI.ObjectiveFunctionType())
        @assert tp_dual_obj !== nothing
        lower_dual_obj = MOI.get(lower_dual, MOI.ObjectiveFunction{tp_dual_obj}())
        # delete dual obj
        # MOI.set(lower_dual, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0))
    end

    # initialize map to lower level model
    lower_dual_idxmap = MOIU.IndexMap()
    # for QP's there are dual variable that are tied to:
    # primal variables
    for (lower_primal_var_key, lower_dual_quad_slack_val) in lower_primal_dual_map.primal_var_dual_quad_slack
        lower_dual_idxmap[lower_dual_quad_slack_val] = lower_idxmap[lower_primal_var_key]
    end
    # and to upper level variable which are lower level parameters
    for (lower_primal_param_key, lower_dual_param_val) in lower_primal_dual_map.primal_parameter
        lower_dual_idxmap[lower_dual_param_val] = lower_idxmap[lower_primal_param_key]
    end
    # Dual variables might appear in the upper level
    for (upper_var, lower_con) in upper_var_lower_ctr
        var = lower_primal_dual_map.primal_con_dual_var[lower_con][1] # TODO check this scalar
        lower_dual_idxmap[var] = upper_idxmap[upper_var]
    end

    # append the second level dual
    append_to(m, lower_dual, lower_dual_idxmap, copy_names)
    pass_names(m, lower_dual, lower_dual_idxmap)

    #=
        Additional Optimiality conditions (to complete the KKT)
    =#

    if ignore_dual_objective(mode)
        # complementary slackness
        comps = get_canonical_complements(lower, lower_primal_dual_map)
        for comp in comps
            if !is_equality(comp.set_w_zero)
                accept_vector_set(mode, comp)
                add_complement(mode, m, comp, lower_idxmap, lower_dual_idxmap)
            else
                # println("eq in complement")
            end
        end
    else
        # strong duality
        lower_dual_obj
        lower_primal_obj
        add_strong_duality(mode, m, lower_primal_obj, lower_dual_obj, lower_idxmap, lower_dual_idxmap)
    end

    return m, upper_idxmap, lower_idxmap, lower_primal_dual_map, lower_dual_idxmap
end

function add_strong_duality(mode::StrongDualityEqualityMode{T}, m, primal_obj, dual_obj,
    idxmap_primal, idxmap_dual) where T

    primal = MOIU.map_indices.(Ref(idxmap_primal), primal_obj)
    dual   = MOIU.map_indices.(Ref(idxmap_dual), dual_obj)

    func = MOIU.operate(-, T, primal, dual)

    c = MOI.add_constraint(m, func, MOI.EqualTo(zero(T)))

    MOI.set(m, MOI.ConstraintName(), c, "lower_strong_duality")

    return c
end
function add_strong_duality(mode::StrongDualityInequalityMode{T}, m, primal_obj, dual_obj,
    idxmap_primal, idxmap_dual) where T

    primal = MOIU.map_indices.(Ref(idxmap_primal), primal_obj)
    dual   = MOIU.map_indices.(Ref(idxmap_dual), dual_obj)

    func = MOIU.operate(-, T, primal, dual)

    func_up = MOIU.operate(-, T, func, mode.epsilon)
    c_up = MOI.add_constraint(m, func_up, MOI.LessThan(zero(T)))
    MOI.set(m, MOI.ConstraintName(), c_up, "lower_strong_duality_up")

    func_lo = MOIU.operate(+, T, func, mode.epsilon)
    c_lo = MOI.add_constraint(m, func_lo, MOI.GreaterThan(zero(T)))
    MOI.set(m, MOI.ConstraintName(), c_lo, "lower_strong_duality_lo")

    return c_up, c_lo
end

function add_complement(mode::ComplementMode{T}, m, comp::Complement, idxmap_primal, idxmap_dual) where T
    f = comp.func_w_cte
    s = comp.set_w_zero
    v = comp.variable

    f_dest = MOIU.map_indices.(Ref(idxmap_primal), f)

    dual = idxmap_dual[v]

    new_f = MOIU.operate(vcat, T, f_dest, MOI.SingleVariable(dual))

    c1 = MOI.add_constraint(m,
        new_f,
        MOI.Complements(1))

    nm = MOI.get(m, MOI.VariableName(), dual)
    MOI.set(m, MOI.ConstraintName(), c1, "compl_compl_($(nm))")

    return c1
end

function add_complement(mode::ComplementWithSlackMode{T}, m, comp::Complement, idxmap_primal, idxmap_dual) where T
    f = comp.func_w_cte
    s = comp.set_w_zero
    v = comp.variable

    slack, slack_in_set = MOI.add_constrained_variable(m, s)
    f_dest = MOIU.map_indices.(Ref(idxmap_primal), f)
    new_f = MOIU.operate(-, T, f_dest, MOI.SingleVariable(slack))
    equality = MOIU.normalize_and_add_constraint(m, new_f, MOI.EqualTo(zero(T)))

    # add start value to slack
    val = MOIU.eval_variables(
        x-> nothing_to_nan(MOI.get(m, MOI.VariablePrimalStart(), x)), f_dest)
    if !isnan(val)
        MOI.set(m, MOI.VariablePrimalStart(), slack, val)
    end

    dual = idxmap_dual[v]
    c1 = MOI.add_constraint(m,
        MOI.VectorOfVariables([slack, dual]),
        MOI.Complements(1))

    nm = MOI.get(m, MOI.VariableName(), dual)
    MOI.set(m, MOI.VariableName(), slack, "slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), slack_in_set, "bound_slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), equality, "eq_slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), c1, "compl_complWslk_($(nm))")

    return slack, slack_in_set, equality, c1
end

function add_complement(mode::SOS1Mode{T}, m, comp::Complement, idxmap_primal, idxmap_dual) where T
    f = comp.func_w_cte
    s = comp.set_w_zero
    v = comp.variable

    slack, slack_in_set = MOI.add_constrained_variable(m, s)
    f_dest = MOIU.map_indices.(Ref(idxmap_primal), f)
    new_f = MOIU.operate(-, T, f_dest, MOI.SingleVariable(slack))
    equality = MOIU.normalize_and_add_constraint(m, new_f, MOI.EqualTo(zero(T)))

    dual = idxmap_dual[v]
    c1 = MOI.add_constraint(m,
        MOI.VectorOfVariables([slack, dual]),
        MOI.SOS1([1.0, 2.0]))

    nm = MOI.get(m, MOI.VariableName(), dual)
    MOI.set(m, MOI.VariableName(), slack, "slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), slack_in_set, "bound_slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), equality, "eq_slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), c1, "compl_sos1_($(nm))")

    return slack, slack_in_set, equality, c1
end

function add_complement(mode::PositiveSOS1Mode{T}, m, comp::Complement, idxmap_primal, idxmap_dual) where T
    f = comp.func_w_cte
    s = comp.set_w_zero
    v = comp.variable

    f_dest = MOIU.map_indices.(Ref(idxmap_primal), f)

    if typeof(s) <: MOI.LessThan # 0
        # requires flipping
        # flipped slack
        slack, slack_in_set = MOI.add_constrained_variable(m, MOI.GreaterThan{T}(0.0))
        new_f = MOIU.operate(+, T, f_dest, MOI.SingleVariable(slack))
        # flipped dual
        real_dual = idxmap_dual[v]
        dual, dual_in_set = MOI.add_constrained_variable(m, MOI.GreaterThan{T}(0.0))
        # dual + real_dual == 0
        opposite = MOIU.normalize_and_add_constraint(m,
            MOI.ScalarAffineFunction(
                [MOI.ScalarAffineTerm(1.0, real_dual),
                 MOI.ScalarAffineTerm(1.0, dual)],
                0.0),
            MOI.EqualTo(zero(T)))

        nm = MOI.get(m, MOI.VariableName(), dual)
        MOI.set(m, MOI.VariableName(), slack, "flip_dual_($(nm))")
        MOI.set(m, MOI.ConstraintName(), slack_in_set, "flip_dual_in_set_($(nm))")
        MOI.set(m, MOI.ConstraintName(), opposite, "flip_dual_eq_($(nm))")
    elseif typeof(s) <: MOI.GreaterThan # 0
        slack, slack_in_set = MOI.add_constrained_variable(m, s)
        new_f = MOIU.operate(-, T, f_dest, MOI.SingleVariable(slack))
        dual = idxmap_dual[v]
    else
        error("Unexpected set type: $s, while building complment constraints.")
    end

    equality = MOIU.normalize_and_add_constraint(m, new_f, MOI.EqualTo(zero(T)))

    c1 = MOI.add_constraint(m,
        MOI.VectorOfVariables([slack, dual]),
        MOI.SOS1([1.0, 2.0]))

    nm = MOI.get(m, MOI.VariableName(), dual)
    MOI.set(m, MOI.VariableName(), slack, "slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), slack_in_set, "bound_slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), equality, "eq_slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), c1, "compl_sos1_($(nm))")

    return slack, slack_in_set, equality, c1
end

function flip_set(set::MOI.LessThan{T}) where T
    return MOI.GreaterThan{T}(0.0)
end
function flip_set(set::MOI.GreaterThan{T}) where T
    return MOI.LessThan{T}(0.0)
end
is_equality(set::S) where {S<:MOI.AbstractSet} = false
is_equality(set::MOI.EqualTo{T}) where T = true
is_equality(set::MOI.Zeros) = true

only_variable_functions(v::MOI.VariableIndex) = MOI.SingleVariable(v)
only_variable_functions(v::Vector{MOI.VariableIndex}) = MOI.VectorOfVariables(v)

function add_complement(mode::ProductMode{T}, m, comp::Complement, idxmap_primal, idxmap_dual) where T
    f = comp.func_w_cte
    s = comp.set_w_zero
    v = comp.variable

    eps = mode.epsilon

    f_dest = MOIU.map_indices(x->idxmap_primal[x], f)

    dual = comp.is_vec ? map(x->idxmap_dual[x], v) : idxmap_dual[v]

    new_f = MOIU.operate(dot, T, f_dest, only_variable_functions(dual))
    new_f1 = MOIU.operate(-, T, new_f, eps)
    c1 = MOI.add_constraint(m,
        new_f1,
        MOI.LessThan{T}(0.0))
    if comp.is_vec # conic
        new_f2 = MOIU.operate(+, T, new_f, eps)
        c2 = MOI.add_constraint(m,
            new_f2,
            MOI.GreaterThan{T}(0.0))
    end

    # c1 = MOI.add_constraint(m,
    #     new_f,
    #     MOI.EqualTo(zero(T)))

    nm = comp.is_vec ? MOI.get.(m, MOI.VariableName(), dual) : MOI.get(m, MOI.VariableName(), dual)

    MOI.set(m, MOI.ConstraintName(), c1, "compl_prod_($(nm))")
    if comp.is_vec
        MOI.set(m, MOI.ConstraintName(), c2, "compl_prod2_($(nm))")
    end

    return c1
end

nothing_to_nan(val) = ifelse(val === nothing, NaN, val)

function add_complement(mode::ProductWithSlackMode{T}, m, comp::Complement, idxmap_primal, idxmap_dual) where T
    f = comp.func_w_cte
    s = comp.set_w_zero
    v = comp.variable

    eps = mode.epsilon

    slack, slack_in_set = MOI.add_constrained_variable(m, s)
    f_dest = MOIU.map_indices.(Ref(idxmap_primal), f)
    new_f = MOIU.operate(-, T, f_dest, MOI.SingleVariable(slack))
    equality = MOIU.normalize_and_add_constraint(m, new_f, MOI.EqualTo(zero(T)))

    # add start value to slack
    val = MOIU.eval_variables(
        x-> nothing_to_nan(MOI.get(m, MOI.VariablePrimalStart(), x)), f_dest)
    if !isnan(val)
        MOI.set(m, MOI.VariablePrimalStart(), slack, val)
    end

    dual = idxmap_dual[v]

    prod_f = MOIU.operate(*, T, MOI.SingleVariable(slack), MOI.SingleVariable(dual))

    prod_f1 = MOIU.operate(-, T, prod_f, eps)
    c1 = MOI.add_constraint(m,
        prod_f1,
        MOI.LessThan{Float64}(0.0))
    if comp.is_vec # conic
        prod_f2 = MOIU.operate(+, T, prod_f, eps)
        c2 = MOI.add_constraint(m,
            prod_f2,
            MOI.GreaterThan{Float64}(0.0))
        MOI.set(m, MOI.ConstraintName(), c1, "compl_prodWslk2_($(nm))")
    end

    # c1 = MOI.add_constraint(m,
    #     prod_f1,
    #     MOI.EqualTo(zero(T)))

    # TODO dont add slack if single variable
    # TODO use dot product

    nm = MOI.get(m, MOI.VariableName(), dual)
    MOI.set(m, MOI.VariableName(), slack, "slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), slack_in_set, "bound_slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), equality, "eq_slk_($(nm))")
    MOI.set(m, MOI.ConstraintName(), c1, "compl_prodWslk_($(nm))")

    return slack, slack_in_set, equality, c1
end

function add_complement(mode::IndicatorMode{T}, m, comp::Complement, idxmap_primal, idxmap_dual) where T
    f = comp.func_w_cte
    s = comp.set_w_zero
    v = comp.variable

    method = mode.mode

    f_dest = MOIU.map_indices(x->idxmap_primal[x], f)

    dual = idxmap_dual[v]

    if comp.is_vec
        @show s
        error("vector constraint not supported")
    end

    nm = comp.is_vec ? MOI.get.(m, MOI.VariableName(), dual) : MOI.get(m, MOI.VariableName(), dual)

    vb1 = MOI.add_variable(m)
    MOI.set(m, MOI.VariableName(), vb1, "compl_bin1_($(nm))")

    cb1 = MOI.add_constraint(m, SVF(vb1), MOI.ZeroOne())
    if method == ONE_ONE || method == ZERO_ZERO
        # second binary
        vb2 = MOI.add_variable(m)
        cb2 = MOI.add_constraint(m, SVF(vb2), MOI.ZeroOne())
        MOI.set(m, MOI.VariableName(), vb2, "compl_bin2_($(nm))")

        # z1 + z2 == 1
        fb = MOIU.operate(+, T, SVF(vb1), SVF(vb2))
        cb = MOI.add_constraint(m, fb, MOI.EqualTo{T}(one(T)))
        MOI.set(m, MOI.ConstraintName(), cb, "compl_sum_bin_($(nm))")
    else
        vb2 = vb1
    end

    f1 = MOIU.operate(vcat, T, SVF(vb1), f_dest)
    f2 = MOIU.operate(vcat, T, SVF(vb2), SVF(dual))

    if method == ONE_ONE
        s1 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.EqualTo(zero(T)))
        s2 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.EqualTo(zero(T)))
    elseif method == ZERO_ZERO
        s1 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}(MOI.EqualTo(zero(T)))
        s2 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}(MOI.EqualTo(zero(T)))
    else
        s1 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.EqualTo(zero(T)))
        s2 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}(MOI.EqualTo(zero(T)))
    end

    c1 = MOI.add_constraint(m, to_vector_affine(f1), s1)
    c2 = MOI.add_constraint(m, to_vector_affine(f2), s2)

    MOI.set(m, MOI.ConstraintName(), c1, "compl_ind1_($(nm))")
    MOI.set(m, MOI.ConstraintName(), c2, "compl_ind2_($(nm))")

    return c1
end

function to_vector_affine(f::MOI.VectorAffineFunction{T}) where T
    return f
end
function to_vector_affine(f::MOI.VectorOfVariables)
    return MOI.VectorAffineFunction{Float64}(f)
end

function pass_names(dest, src, map)
    for vi in MOI.get(src, MOI.ListOfVariableIndices())
        name = MOI.get(src, MOI.VariableName(), vi)
        if name != ""
            MOI.set(dest, MOI.VariableName(), map[vi], name)
        end
    end
    for (F,S) in MOI.get(src, MOI.ListOfConstraints())
        for con in MOI.get(src, MOI.ListOfConstraintIndices{F,S}())
            name = MOI.get(src, MOI.ConstraintName(), con)
            if name != ""
                MOI.set(dest, MOI.ConstraintName(), map[con], name)
            end
        end
    end
end

function append_to(dest::MOI.ModelLike, src::MOI.ModelLike, idxmap, copy_names::Bool; allow_single_bounds::Bool = true)
    # MOI.empty!(dest)

    # idxmap = MOIU.IndexMap()

    vis_src = MOI.get(src, MOI.ListOfVariableIndices())
    constraint_types = MOI.get(src, MOI.ListOfConstraints())
    single_variable_types = [S for (F, S) in constraint_types
                             if F == MOI.SingleVariable && allow_single_bounds]
    vector_of_variables_types = [S for (F, S) in constraint_types
                                 if F == MOI.VectorOfVariables]

    # The `NLPBlock` assumes that the order of variables does not change (#849)
    if MOI.NLPBlock() in MOI.get(src, MOI.ListOfModelAttributesSet())
        error("not nlp for now")
        vector_of_variables_not_added = [
            MOI.get(src, MOI.ListOfConstraintIndices{MOI.VectorOfVariables, S}())
            for S in vector_of_variables_types
        ]
        single_variable_not_added = [
            MOI.get(src, MOI.ListOfConstraintIndices{MOI.SingleVariable, S}())
            for S in single_variable_types
        ]
    else
        vector_of_variables_not_added = [
            MOIU.copy_vector_of_variables(dest, src, idxmap, S)
            for S in vector_of_variables_types
        ]
        single_variable_not_added = [
            MOIU.copy_single_variable(dest, src, idxmap, S)
            for S in single_variable_types
        ]
    end

    # MOIU.copy_free_variables(dest, idxmap, vis_src, MOI.add_variables)
    for vi in vis_src
        if !haskey(idxmap.varmap, vi)
            var = MOI.add_variable(dest)
            idxmap.varmap[vi] = var
        end
    end

    # Copy variable attributes
    MOIU.pass_attributes(dest, src, copy_names, idxmap, vis_src)

    # Copy model attributes
    # attention HERE to no pass objective functions!
    # pass_attributes(dest, src, copy_names, idxmap)

    # Copy constraints
    MOIU.pass_constraints(dest, src, copy_names, idxmap,
                     single_variable_types, single_variable_not_added,
                     vector_of_variables_types, vector_of_variables_not_added)

    return idxmap
end

using LinearAlgebra

# scalar
function MOIU.promote_operation(::typeof(LinearAlgebra.dot), ::Type{T},
    ::Type{<:Union{MOI.SingleVariable, MOI.ScalarAffineFunction{T}}},
    ::Type{T}
    ) where T
    MOI.ScalarAffineFunction{T}
end
function MOIU.promote_operation(::typeof(LinearAlgebra.dot), ::Type{T},
    ::Type{T},
    ::Type{<:Union{MOI.SingleVariable, MOI.ScalarAffineFunction{T}}}
    ) where T
    MOI.ScalarAffineFunction{T}
end
function MOIU.promote_operation(::typeof(LinearAlgebra.dot), ::Type{T},
    ::Type{<:Union{MOI.SingleVariable, MOI.ScalarAffineFunction{T}}},
    ::Type{<:Union{MOI.SingleVariable, MOI.ScalarAffineFunction{T}}}
    ) where T
    MOI.ScalarQuadraticFunction{T}
end
function MOIU.promote_operation(::typeof(LinearAlgebra.dot), ::Type{T},
    ::Type{MOI.ScalarQuadraticFunction{T}},
    ::Type{T}
    ) where T
    MOI.ScalarQuadraticFunction{T}
end
function MOIU.promote_operation(::typeof(LinearAlgebra.dot), ::Type{T},
    ::Type{T},
    ::Type{MOI.ScalarQuadraticFunction{T}}
    ) where T
    MOI.ScalarQuadraticFunction{T}
end
# flip
function MOIU.operate(::typeof(LinearAlgebra.dot), ::Type{T},
    f::Union{
        MOI.SingleVariable,
        MOI.ScalarAffineFunction{T},
        MOI.ScalarQuadraticFunction{T}
        },
    α::T) where T
    return MOIU.operate(LinearAlgebra.dot, T, α, f)
end
# pass to *
function MOIU.operate(::typeof(LinearAlgebra.dot), ::Type{T},
    f::Union{
        T,
        MOI.SingleVariable,
        MOI.ScalarAffineFunction{T}
        },
    g::Union{
        MOI.SingleVariable,
        MOI.ScalarAffineFunction{T}
        }
    ) where T
    return MOIU.operate(*, T, f, g)
end
function MOIU.operate(::typeof(LinearAlgebra.dot), ::Type{T},
    α::T,
    f::MOI.ScalarQuadraticFunction{T}
    ) where T
    return MOIU.operate(*, T, f, α)
end

# vector
function MOIU.promote_operation(::typeof(LinearAlgebra.dot), ::Type{T},
    ::Type{<:Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}}},
    ::Type{Vector{T}}
    ) where T
    MOI.VectorAffineFunction{T}
end
function MOIU.promote_operation(::typeof(LinearAlgebra.dot), ::Type{T},
    ::Type{Vector{T}},
    ::Type{<:Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}}}
    ) where T
    MOI.VectorAffineFunction{T}
end
function MOIU.promote_operation(::typeof(LinearAlgebra.dot), ::Type{T},
    ::Type{<:Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}}},
    ::Type{<:Union{MOI.VectorOfVariables, MOI.VectorAffineFunction{T}}}
    ) where T
    MOI.VectorQuadraticFunction{T}
end
function MOIU.promote_operation(::typeof(LinearAlgebra.dot), ::Type{T},
    ::Type{MOI.VectorQuadraticFunction{T}},
    ::Type{Vector{T}}
    ) where T
    MOI.VectorQuadraticFunction{T}
end
function MOIU.promote_operation(::typeof(LinearAlgebra.dot), ::Type{T},
    ::Type{Vector{T}},
    ::Type{MOI.VectorQuadraticFunction{T}}
    ) where T
    MOI.VectorQuadraticFunction{T}
end
# flip
function MOIU.operate(::typeof(LinearAlgebra.dot), ::Type{T},
    f::Union{
        MOI.VectorOfVariables,
        MOI.VectorAffineFunction{T},
        MOI.VectorQuadraticFunction{T}
        },
    α::Vector{T}) where T
    return MOIU.operate(LinearAlgebra.dot, T, α, f)
end
# pass to _operate(LinearAlgebra.dot, ...)
function MOIU.operate(::typeof(LinearAlgebra.dot), ::Type{T},
    f::Union{
        Vector{T},
        MOI.VectorOfVariables,
        MOI.VectorAffineFunction{T}
        },
    g::Union{
        MOI.VectorOfVariables,
        MOI.VectorAffineFunction{T}
        }
    ) where T
    return _operate(LinearAlgebra.dot, T, f, g)
end
function MOIU.operate(::typeof(LinearAlgebra.dot), ::Type{T},
    α::T,
    f::MOI.VectorQuadraticFunction{T}
    ) where T
    return _operate(LinearAlgebra.dot, T, f, α)
end
function _operate(::typeof(LinearAlgebra.dot), ::Type{T},
    f::Union{
        Vector{T},
        MOI.VectorOfVariables,
        MOI.VectorAffineFunction{T},
        MOI.VectorQuadraticFunction{T}
        },
    g::Union{
        MOI.VectorOfVariables,
        MOI.VectorAffineFunction{T},
        MOI.VectorQuadraticFunction{T}
    }) where T

    dim = MOI.output_dimension(g)
    if MOI.output_dimension(f) != dim
        throw(DimensionMismatch("f and g are of different MOI.output_dimension's!"))
    end

    fs = MOIU.scalarize(f)
    gs = MOIU.scalarize(g)

    out = MOIU.operate(*, T, fs[1], gs[1])
    for i in 2:dim
        MOIU.operate!(+, T, out, MOIU.operate(*, T, fs[i], gs[i]))
    end

    return out
end
MOIU.scalarize(v::Vector{T}) where T<:Number = v
MOI.output_dimension(v::Vector{T}) where T<:Number = length(v)


text:
using JSON
using CSV
using Dates
using Dates: DateFormat
using NamedArrays

using LazySets
# using Colors
using Combinatorics
using IOCapture

using JuMP

ext: jl


text:
include("lm_getindex.jl")
include("lm-hcat.jl")


text:
using SphericalBesselFunctions
using Test

include("reference.jl")

@testset "SphericalBesselFunctions.jl" begin
    @testset "Accuracy" begin
        nx = 1001
        x = 10 .^ range(-1, stop=4, length=nx)
        nℓ = 105
        j, j′, y, y′ = @time bessels(x, nℓ)

        sel = 1:min(nℓ,104)

        nℓref = min(nℓ,length(sel))
        jref = zeros(nx, nℓref)
        j′ref = zeros(nx, nℓref)
        yref = zeros(nx, nℓref)
        y′ref = zeros(nx, nℓref)

        @time for (i,x) in enumerate(x)
            bessels_ref!(view(jref, i, :),
                         view(j′ref, i, :),
                         view(yref, i, :),
                         view(y′ref, i, :),
                         x)
        end

        j = j[:,sel]
        j′ = j′[:,sel]
        y = y[:,sel]
        y′ = y′[:,sel]
        jref = jref[:,sel]
        j′ref = j′ref[:,sel]
        yref = yref[:,sel]
        y′ref = y′ref[:,sel]

        @test j ≈ jref
        @test j′ ≈ j′ref

        # The irregular functions disagree wildly for x < ℓ, so we
        # mask that out
        mask = x .< (0:nℓ-2)'
        y[mask] .= 0
        y′[mask] .= 0
        yref[mask] .= 0
        y′ref[mask] .= 0

        # Just to make sure we did not mask out the whole matrix...
        @test any(!iszero, y)
        @test any(!iszero, y′)

        @test y ≈ yref
        @test y′ ≈ y′ref
    end
end


text:
using TensorFlow
using Distributions

# Generate some synthetic data
x = randn(100, 50)
w = randn(50, 10)
y_prob = exp.(x*w)
y_prob ./= sum(y_prob,2)

function draw(probs)
    y = zeros(size(probs))
    for i in 1:size(probs, 1)
        idx = rand(Categorical(probs[i, :]))
        y[i, idx] = 1
    end
    return y
end

y = draw(y_prob)

# Build the model
sess = Session(Graph())

X = placeholder(Float64, shape=[-1, 50])
Y_obs = placeholder(Float64, shape=[-1, 10])

variable_scope("logisitic_model"; initializer=Normal(0, .001)) do
    global W = get_variable("W", [50, 10], Float64)
    global B = get_variable("B", [10], Float64)
end

Y=nn.softmax(X*W + B)

Loss = -reduce_sum(log(Y).*Y_obs)
optimizer = train.AdamOptimizer()
minimize_op = train.minimize(optimizer, Loss)
saver = train.Saver()

# Run training
run(sess, global_variables_initializer())
checkpoint_path = mktempdir()
info("Checkpoint files saved in $checkpoint_path")
for epoch in 1:100
    cur_loss, _ = run(sess, [Loss, minimize_op], Dict(X=>x, Y_obs=>y))
    println(@sprintf("Current loss is %.2f.", cur_loss))
    train.save(saver, sess, joinpath(checkpoint_path, "logistic"), global_step=epoch)
end

max_stars_repo_name: octonion/tensorflow

max_issues_repo_name: octonion/tensorflow

max_forks_repo_name: octonion/tensorflow


text:
function macroexpand(o::Org, m::Macro)
    if m.name in ("title", "author", "email") &&
        haskey(o.keywords, m.name)
        join(o.keywords[m.name], ' ')
    elseif m.name == "keyword" &&
        haskey(o.keywords, m.name)
        join(o.keywords[m.name], ' ')
        # TODO support DATE, Base.Libc.strftime will probaby be useful
        # TODO support n macro ... somehow
    elseif m.name in keys(o.macros)
        o.macros[m.name](m.arguments)
    elseif m.name == "results" && length(m.arguments) == 1
        "=$(m.arguments[1])="
    end
end


text:
import DataStructures: OrderedDict
import StaticArrays: MVector
import Serialization
import SHA

################
# NSitem
################

abstract type AbstNSitem end

struct NScst_item{T} <: AbstNSitem
    obj::T
end

mutable struct NSnoncst_item{T} <: AbstNSitem
    obj::T
end

Base.copy(x::Wild.AbstNSitem) = typeof(x)(x.obj)

################
# NSX{X}
################

struct NSX{X} <: AbstNS
    __dict::OrderedDict{Symbol, AbstNSitem}
    __fix_lck::MVector{2, Bool}

#    NSX{X}() where X = new{X}(#= __dict    =# OrderedDict{Symbol, AbstNSitem}(),
#                              #= __fix_lck =# MVector{2, Bool}(false, false))
    NSX{X}(; ka...) where X =
        begin
            g = new(#= __dict    =# OrderedDict{Symbol, AbstNSitem}(),
                    #= __fix_lck =# MVector{2, Bool}(false, false))
            setprp(g, k, v) =
                begin
                    s = string(k)
                    if length(s) > 3 && s[4] == '_'
                        length(s) == 4 &&
                            Base.error("$(k) can't be used for property.")
                        pre = Symbol(s[1:3])
                        if pre in [:cst, :dfn, :prp, :req, :mth, :fnc]
                            setprp(Base.getproperty(g, pre), Symbol(s[5:end]), v)
                        else
                            Base.setproperty!(g, k, v)
                        end
                    else
                        Base.setproperty!(g, k, v)
                    end
                end

            for (k,v) in pairs(ka)
                setprp(g, k, v)
            end
            g
        end;
end

const NS = NSX{nothing}

abstract type __NSFlgCodeMode end

################
# NS
################

Base.setproperty!(ns::AbstNS, atr::Symbol, x) =
    begin
        hasfield(typeof(ns), atr) && (return Base.setfield!(ns, atr, x))

        atr == :exe && (return x(ns))

        haskey(_NSdict0, atr) &&
            Base.error("""'$(atr)' can't be used for property""")

        d = ns.__dict

        if haskey(d, atr)
            ns._fixed && Base.error("this NS is fixed!")

            o = d[atr].obj
            isa(o, NSPrp) && (o.fnc(ns, x); return)

            isa(d[atr], NScst_item) && Base.error("""'$(atr)' is const.""")
        else
            ns._lcked && Base.error("this NS is locked!")
        end

        y = isa(x, AbstNSitem) ? x : NSnoncst_item(x)

        d[atr] = (isa(y.obj, NSDfn)
                  ? (isa(y, NScst_item) ? NScst_item : NSnoncst_item)(y.obj.fnc(ns))
                  : copy(y))
    end

Base.haskey(o::AbstNS, key::Symbol) = key ∈ o._keys

Base.propertynames(ns::AbstNS, private::Bool=false) =
    (private
     ? tuple(Base.keys(ns.__dict)...,
             Base.keys(_NSdict0)...,
             Base.fieldnames(typeof(ns))...)
     : tuple(Base.keys(ns.__dict)...,
             Base.keys(_NSdict0)...))

Base.hasproperty(ns::AbstNS, atr::Symbol) =
    Base.hasfield(typeof(ns), atr) ||
    haskey(_NSdict0, atr) ||
    haskey(ns.__dict, atr)

Base.getproperty(ns::AbstNS, atr::Symbol) =
    begin
        Base.hasfield(typeof(ns), atr) && (return Base.getfield(ns, atr))

        haskey(_NSdict0, atr) && (return _NSdict0[atr](ns))

        d = ns.__dict

        haskey(d, atr) ||
            error("""this NS does not have a property named '$(atr)'.""")

        x = d[atr].obj;
        isa(x, Union{NSPrp, NSMth}) && (return x(ns))
        isa(x, NSFnc) && (return x.fnc)
        isa(x, NSReq) &&
            (y = x(ns);
             d[atr] = (isa(d[atr], NScst_item) ? NScst_item : NSnoncst_item)(y);
             return y)
        return x
    end

################
# __NSX_CodeMode
################

__NSX_CodeMode_CodeType = Vector{NamedTuple{(:atr, :obj),Tuple{Symbol,Any}}}
struct __NSX_CodeMode <: AbstNS
    __code::__NSX_CodeMode_CodeType
    __instances
    __parallel::NSnoncst_item{Bool}
    __NSX_CodeMode() =
        new(#= __code      =# __NSX_CodeMode_CodeType(),
            #= __instances =# [],
            #= __parallel  =# NSnoncst_item{Bool}(false))
end

__divNn(N::Integer, n::Integer) =
begin
    k1 = (let v = fill(N ÷ n, n); v[1:N%n] .+= 1; accumulate(+, v); end)
    k0 = (let v = similar(k1); v[1] = 1; @. v[2:end] = k1[1:end-1] + 1; v; end)
    [i0:i1 for (i0,i1) in zip(k0,k1)]
end

Base.setproperty!(ns::__NSX_CodeMode, atr::Symbol, x) =
    begin
        hasfield(typeof(ns), atr) && (return Base.setfield!(ns, atr, x))

        haskey(_NSdict0, atr) &&
            Base.error("""'$(atr)' can't be used for property""")

        y = atr == :exe ? x : isa(x, AbstNSitem) ? x : NSnoncst_item(x)

        push!(ns.__code, NamedTuple{(:atr, :obj), Tuple{Symbol, Any}}((atr, y)))

        inst = ns.__instances[1]
        if ns.__parallel.obj
            @sync for r in __divNn(length(inst), Threads.nthreads())
                Threads.@spawn foreach(i -> Base.setproperty!(i.o, atr, x), inst[r])
            end
        else
            foreach(i -> Base.setproperty!(i.o, atr, x), inst)
        end
    end

Base.haskey(ns::__NSX_CodeMode, key::Symbol) = key ∈ propertynames(ns)

Base.propertynames(ns::__NSX_CodeMode, private::Bool=false) =
    tuple(Base.keys(_NSdict0)..., Base.fieldnames(typeof(ns))...)

Base.hasproperty(ns::__NSX_CodeMode, atr::Symbol) =
    Base.hasfield(typeof(ns), atr) || haskey(_NSdict0, atr)

Base.getproperty(ns::__NSX_CodeMode, atr::Symbol) =
    begin
        Base.hasfield(typeof(ns), atr) && (return Base.getfield(ns, atr))
        haskey(_NSdict0, atr) && (return _NSdict0[atr](ns))
        error("""This NS does not have a property named '$(atr)'.""")
    end;

################
# __NSClsInstance{X}
################

struct __NSClsInstance{X} <: AbstNS
    __dict::OrderedDict{Symbol, AbstNSitem}
    __fix_lck::MVector{2, Bool}
    cls::NS

    __NSClsInstance{X}(cls) where X =
        new{X}(#= __dict    =# OrderedDict{Symbol, AbstNSitem}(),
               #= __fix_lck =# MVector{2, Bool}(false, false),
               #= cls       =# cls)
end

Base.getproperty(nsi::__NSClsInstance, atr::Symbol) =
    begin
        Base.hasfield(typeof(nsi), atr) && (return Base.getfield(nsi, atr))

        haskey(_NSdict0, atr) && (return _NSdict0[atr](nsi))

        d = nsi.__dict

        if haskey(d, atr)
            x = d[atr].obj;
            isa(x, Union{NSTagFunc{:prp}, NSTagFunc{:mth}}) && (return x(nsi))
            isa(x, NSTagFunc{:fnc}) && (return x.fnc)
            isa(x, NSTagFunc{:req}) &&
                (y = x(nsi);
                 d[atr] = (isa(d[atr], NScst_item) ? NScst_item : NSnoncst_item)(y);
                 return y)
            return x
        else
            haskey(nsi.cls, atr) && (return Base.getproperty(nsi.cls, atr))
            error("""This NS does not have a property named '$(atr)'.""")
        end
    end

################
# >>, >>>
################

Base.:>>( g::AbstNS, h::AbstNS) = h.import(g)
Base.:>>>(g::AbstNS, h::AbstNS) = h.deepimport(g)

################
# genNSX, ns
################

# genNSX() = NSX{Symbol("NS_", string(bytes2hex(SHA.sha256(string(time_ns())))))}
genNSX() = NSX{gensym()}
genNSX(X) = NSX{X}

nsx() = genNSX()()
nsx(X) = genNSX(X)()


text:
"""
Data pre- / post-processing for sequential problems. Takes a ProfileData object and prepares it for use in GP.
"""

"""
`get_predictors_targets(vavg::Array, problem::Sequential_T)`
# Description
    Returns x and y, the scaled predictor and target pairs from which to extract the training and verification data sets for "T" profiles.

    model( predictor ) -> target
         model( T[i] ) -> (T[i+1]-T[i])/Δt ≈ ∂t(T)

# Arguments
- `vavg`: (Array)                Nt-length array of D-length vectors. Data from which to extract x and y, the predictors and corresponding predictions.
- `problem`: (Sequential_dT)     Sequential_T object associated with the data (output of get_problem)

# *--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*
# | Sequential_dT                              |
# |                                            |
# |   predictor       target                   |
# |   T[i] --model--> (T[i+1]-T[i])/Δt ≈ ∂t(T) |
# |                                            |
# |   model output                             |
# |   model(T[i])                              |
# |                                            |
# |   prediction                               |
# |   predicted T[i+1] = model(T[i])Δt + T[i]  |
# *--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*

"""
function get_predictors_targets(vavg::Array, problem::Sequential_dT)

    # scale according to initial temperature profile vavg[1]
    scaling = min_max_scaling(vavg[1][end]-vavg[1][1], minimum(vavg[1]))
    vavg = [scale(vec, scaling) for vec in vavg]

    # store the scaling as an attribute in the problem for later use
    problem.scaling = scaling

    # calculate predictors and targets from the scaled data
    predictors = vavg[1:end-1] # T[i] for i = 1, ..., Nt-1
    targets = (vavg[2:end] - predictors) / problem.Δt # (T[i+1]-T[i])/Δt'

    return (predictors, targets)
end

"""
`postprocess_prediction(predictor, model_output, problem::Sequential_dT)`

# Description
Takes in a scaled predictor, T[i], the scaled model output on the predictor, model(T[i]), and a Sequential_T object.
Returns the unscaled prediction (predicted temperature profile), T[i+1], computed from T[i] and model(T[i]) by

    predicted T[i+1] = model(T[i])Δt + T[i]

# Arguments
- `predictor`: (Array)            T[i], the scaled predictor for a temperature profile
- `model_output`: (Array)         model(T[i]), the scaled model output for a temperature profile
- `problem`: (Sequential_dT)
"""
function postprocess_prediction(predictor, model_output, problem::Sequential_dT)

    prediction = model_output * problem.Δt + predictor
    prediction = [unscale(vec, problem.scaling) for vec in prediction]

    return prediction

end

"""
`get_predictors_targets(vavg::Array, problem::Sequential_T)`

# Description
    Returns x and y, the predictors and target predictions from which to extract the training and verification data for "T" profiles.

    model( predictors ) -> targets
          model( T[i] ) -> T[i+1]

# Arguments
- `vavg`: (Array)               Nt-length array of D-length vectors. Data from which to extract x and y, the predictors and corresponding predictions.
- `problem`: (Sequential_T)     Sequential_T object associated with the data (output of get_problem)

# *--*--*--*--*--*--*--*--*--*--*
# | Sequential_T                |
# |   T[i] --model--> T[i+1]    |
# *--*--*--*--*--*--*--*--*--*--*
"""
function get_predictors_targets(vavg, problem::Sequential_T)

    # scale according to initial temperature profile vavg[1]
    scaling = min_max_scaling(vavg[1][end]-vavg[1][1], minimum(vavg[1]))
    vavg = [scale(vec, scaling) for vec in vavg]

    # store the scaling as an attribute in the problem for later use
    problem.scaling = scaling # this mutates the input

    # get predictors and targets
    predictors = vavg[1:end-1] # T[i] for i = 1, ..., Nt-1
    targets = vavg[2:end] # T[i+1]

    return (predictors, targets)
end

"""
`postprocess_prediction(predictor, prediction, problem::Sequential_T)`

# Description
Takes in a scaled predictor, T[i], the scaled GP prediction on the predictor, G(T[i]), and a Sequential_T object.
Returns the predicted temperature profile, T[i+1], computed from T[i] and G(T[i]) by

           prediction = model( predictor )
     predicted T[i+1] = model( T[i] )

# Arguments
- `predictor`: (Array)            T[i], the scaled predictor for a temperature profile
- `prediction`: (Array)           model(T[i), the scaled prediction for a temperature profile
- `problem`: (Sequential_T)
"""
function postprocess_prediction(predictor, prediction, problem::Sequential_T)
    # println("prediction$([unscale(vec, problem.scaling) for vec in prediction])")

    return [unscale(vec, problem.scaling) for vec in prediction]
end

"""
`get_predictors_targets(vavg, problem::Sequential_wT)`

# Description
    Returns x and y, the predictors and targets from which to extract the training and verification data for "wT" profiles.

     model( predictor ) -> target
         model( wT[i] ) -> wT[i+1]

# Arguments
- `vavg`: (Array)                  Nt-length array of D-length vectors. Data from which to extract x and y, the predictors and corresponding predictions.
- `problem`: (Sequential_wT)       Sequential_wT object associated with the data (output of get_problem)

# *--*--*--*--*--*--*--*--*--*--*
# | Sequential_wT               |
# |   wT[i] --model--> wT[i+1]  |
# *--*--*--*--*--*--*--*--*--*--*
"""
function get_predictors_targets(vavg, problem::Sequential_wT)

    # scale according to initial temperature profile vavg[1]
    # scaling = min_max_scaling(vavg[1][end]-vavg[1][1], minimum(vavg[1]))
    # vavg = [scale(vec, scaling) for vec in vavg]

    # store the scaling as an attribute in the problem for later use
    # problem.scaling = scaling # this mutates the input

    # get predictors and predictions
    predictors = vavg[1:end-1] # wT[i] for i = 1, ..., Nt-1
    targets = vavg[2:end] # wT[i+1]

    return (predictors, targets)
end

"""
`postprocess_prediction(predictor, prediction, problem::Sequential_wT)`

# Description
Takes in a scaled predictor, wT[i], the scaled model prediction on the predictor, model(wT[i]), and a Sequential_wT object.
Returns the temperature profile, T[i+1], computed from model(T[i]) by

           prediction = model(predictor)
    predicted wT[i+1] = model(wT[i])

# Arguments
- `predictor`: (Array)                wT[i], the predictor for a wT profile
- `prediction`: (Array)               predicted wT[i+1], the model prediction for a wT profile
- `problem`: (Sequential_wT)          Sequential_wT object associated with the data
"""
function postprocess_prediction(predictor, prediction, problem::Sequential_wT)
    # prediction = [unscale(vec, problem.scaling) for vec in prediction]
    return prediction
end

"""
get_predictors_targets(vavg::Array, problem::Residual_KPP)
# Description
    Returns x and y, the scaled predictors and target predictions from which to extract the training and verification data for temperature profiles.
    Scales the predictors and targets using min-max scaling based on the initial temperature profile from the les simulation.

    model( predictors ) = targets
    model( KPP(T[i])  ) = T[i] - KPP(T[i])

# Arguments
- `vavg`: (Array)               Nt-length array of D-length vectors. Data from which to extract x and y, the predictors and corresponding predictions.
-  `problem`: (Residual_KPP)    Residual_KPP object associated with the data (output of get_problem)

# *--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*
# | Residual_KPP                                     |
# |                                                  |
# |   predictor            target                    |
# |   KPP(T[i]) --model--> T[i] - KPP(T[i])          |
# |                                                  |
# |   model output                                   |
# |   model(KPP(T[i]))                               |
# |                                                  |
# |   predicted T[i] = model(KPP(T[i])) + KPP(T[i])  |
# *--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*

# *--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*
# | Residual_TKE                                     |
# |                                                  |
# |   predictor            target                    |
# |   TKE(T[i]) --model--> T[i] - TKE(T[i])          |
# |                                                  |
# |   model output                                   |
# |   model(TKE(T[i]))                               |
# |                                                  |
# |   predicted T[i] = model(TKE(T[i])) + TKE(T[i])  |
# *--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*--*
"""
function get_predictors_targets(vavg::Array, problem::Union{Sequential_KPP, Sequential_TKE})

    # scale according to initial temperature profile vavg[1]
    scaling = min_max_scaling(vavg[1][end]-vavg[1][1], minimum(vavg[1]))
    vavg = [scale(vec, scaling) for vec in vavg]

    # store the scaling as an attribute in the problem for later use
    problem.scaling = scaling

    # scale physics data with same scaling
    physics_data = [scale(vec, scaling) for vec in problem.physics_data] # kpp(i; T[i-1])

    targets = (vavg .- physics_data) / problem.Δt # residual
    return (vavg, targets)
end

"""
# Description
Takes in a predictor, T[i], the GP prediction on the predictor, G(T[i]), and a Residual_KPP object.
Returns the predicted temperature profile, T[i+1], computed from T[i] and G(T[i]) by

           prediction = model(predictor) # residual
     predicted T[i+1] = model(KPP(T[i])) + KPP(T[i])

# Arguments
- `predictor`: (Array)            T[i], the predictor for a temperature profile
- `prediction`: (Array)           model(T[i), the prediction for a temperature profile
- `problem`: (Residual_KPP)
"""
function postprocess_prediction(predictor, prediction, problem::Union{Sequential_KPP, Sequential_TKE})

    #unscale predictor
    predictor = [unscale(vec, problem.scaling) for vec in predictor] # kpp data

    #unscale prediction
    prediction = [(vec * problem.scaling.Δv) for vec in prediction] # residual

    return predictor .+ prediction * problem.Δt
end

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
module BoysFunction

using QuadGK
using SpecialFunctions

include("functions.jl")
include("recursions.jl")

export boys_gamma,
        boys_quadgk,
        boys_rec_largex,
        downward_recursion

# Write your package code here.

end


text:
using RVM
using Base.Test

@test hello("Tump") == "Hello, Tump!"


text:
# __precompile__()

module LightDarkLPDM

# importall POMDPs
using POMDPs
using StaticArrays
using Combinatorics
using Distributions
# using Plots
# using POMDPToolbox
using Parameters # for @with_kw
using Printf
using Random
# using ParticleFilters # for AbstractParticleBelief
using LPDM

include("1D/lightdark1d.jl")
include("2D/lightdark2d.jl")
# include("lightdark2dtarget.jl")
# include("lightdark2dfilter.jl")
# include("lightdark1ddespot.jl")
include("1D/lightdark1dlpdm.jl")
# include("lightdark2ddespot.jl")
include("2D/lightdark2dlpdm.jl")

# include("lightdark2dvis.jl")
export
    AbstractLD1,
    AbstractLD2,
    LightDark1D,
    LightDark2D,
    LightDark2DTarget,
    # LightDark1DDespot,
    LightDark1DLpdm,
    # LightDark2DDespot,
    LightDark2DLpdm,
    LightDark2DKalman,
    Normal1D,
    SymmetricNormal2D,
    Vec2,
    obs_std,
    rand,
    state_distribution,
    generate_sor_det,
    LD1State,
    LD1Action,
    LD1Obs,
    LD2State,
    LD2Action,
    LD2Obs
end # module

max_stars_repo_name: ebalaban/LightDarkLPDM

max_issues_repo_name: tomerarnon/LightDarkLPDM

max_forks_repo_name: tomerarnon/LightDarkLPDM


text:
# Autogenerated wrapper script for ODEInterface_jll for i686-w64-mingw32-libgfortran5
export libbvp_m_proxy, libbvpsol, libbvpsol_i32, libcolnew, libcolnew_i32, libddeabm, libddeabm_i32, libddebdf, libddebdf_i32, libdop853, libdop853_i32, libdopri5, libdopri5_i32, libodex, libodex_i32, libradau, libradau5, libradau5_i32, libradau_i32, librodas, librodas_i32, libseulex, libseulex_i32

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("ODEInterface")
JLLWrappers.@declare_library_product(libbvp_m_proxy, "bvp_m_proxy.dll")
JLLWrappers.@declare_library_product(libbvpsol, "bvpsol.dll")
JLLWrappers.@declare_library_product(libbvpsol_i32, "bvpsol_i32.dll")
JLLWrappers.@declare_library_product(libcolnew, "colnew.dll")
JLLWrappers.@declare_library_product(libcolnew_i32, "colnew_i32.dll")
JLLWrappers.@declare_library_product(libddeabm, "ddeabm.dll")
JLLWrappers.@declare_library_product(libddeabm_i32, "ddeabm_i32.dll")
JLLWrappers.@declare_library_product(libddebdf, "ddebdf.dll")
JLLWrappers.@declare_library_product(libddebdf_i32, "ddebdf_i32.dll")
JLLWrappers.@declare_library_product(libdop853, "dop853.dll")
JLLWrappers.@declare_library_product(libdop853_i32, "dop853_i32.dll")
JLLWrappers.@declare_library_product(libdopri5, "dopri5.dll")
JLLWrappers.@declare_library_product(libdopri5_i32, "dopri5_i32.dll")
JLLWrappers.@declare_library_product(libodex, "odex.dll")
JLLWrappers.@declare_library_product(libodex_i32, "odex_i32.dll")
JLLWrappers.@declare_library_product(libradau, "radau.dll")
JLLWrappers.@declare_library_product(libradau5, "radau5.dll")
JLLWrappers.@declare_library_product(libradau5_i32, "radau5_i32.dll")
JLLWrappers.@declare_library_product(libradau_i32, "radau_i32.dll")
JLLWrappers.@declare_library_product(librodas, "rodas.dll")
JLLWrappers.@declare_library_product(librodas_i32, "rodas_i32.dll")
JLLWrappers.@declare_library_product(libseulex, "seulex.dll")
JLLWrappers.@declare_library_product(libseulex_i32, "seulex_i32.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libbvp_m_proxy,
        "bin\\bvp_m_proxy.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libbvpsol,
        "bin\\bvpsol.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libbvpsol_i32,
        "bin\\bvpsol_i32.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcolnew,
        "bin\\colnew.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcolnew_i32,
        "bin\\colnew_i32.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libddeabm,
        "bin\\ddeabm.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libddeabm_i32,
        "bin\\ddeabm_i32.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libddebdf,
        "bin\\ddebdf.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libddebdf_i32,
        "bin\\ddebdf_i32.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdop853,
        "bin\\dop853.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdop853_i32,
        "bin\\dop853_i32.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdopri5,
        "bin\\dopri5.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdopri5_i32,
        "bin\\dopri5_i32.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libodex,
        "bin\\odex.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libodex_i32,
        "bin\\odex_i32.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libradau,
        "bin\\radau.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libradau5,
        "bin\\radau5.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libradau5_i32,
        "bin\\radau5_i32.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libradau_i32,
        "bin\\radau_i32.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        librodas,
        "bin\\rodas.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        librodas_i32,
        "bin\\rodas_i32.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libseulex,
        "bin\\seulex.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libseulex_i32,
        "bin\\seulex_i32.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# This file was generated, do not modify it. # hide
predict(mach, rows=71:73)

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
push!(LOAD_PATH, "../src/")

DOCUMENTER_DEBUG=true

using Documenter, HnswAnn

makedocs(
    format = Documenter.HTML(prettyurls = false),
    sitename = "HnswAnn.jl",
    pages = Any[
        "Introduction" => "INTRO.md",
        "API documentation" => "api.md",
        "Internals" => "internals.md"
    ]
)


text:
# Gradient Bandits
"""
    Gradient Bandit

    Based on: Sec 2.8, R. S. Sutton, A. G. Barto, and A. B. Book, Reinforcement Learning : An Introduction, Second Edi. The MIt Press, 2017.
"""
mutable struct GradientBandit <: BanditAlgorithmBase
    noOfArms::Int64
    noOfSteps::Int64
    lastPlayedArm::Int64

    α::Float64
    avgValue::Float64
    preference::Vector{Float64}
    pDist::Categorical

    function GradientBandit( noOfArms::Integer, α::Real )
        new( noOfArms,
             0,
             0,
             α,
             0,
             zeros( Float64, noOfArms ),
             Categorical( 1/noOfArms*ones(noOfArms) ) )
    end
end

function get_arm_index( agent::GradientBandit )
    agent.lastPlayedArm = rand( agent.pDist )
    return agent.lastPlayedArm
end

function update_reward!( agent::GradientBandit, r::Real )
    agent.avgValue = (agent.avgValue*agent.noOfSteps+r)/(agent.noOfSteps+1)
    agent.noOfSteps += 1

    # Update for all arms
    agent.preference -= (agent.α*(r-agent.avgValue).*agent.pDist.p)
    # Perform update for selected arm
    agent.preference[agent.lastPlayedArm] += agent.α*(r-agent.avgValue)

    agent.pDist = Categorical( exp.(agent.preference)/sum(exp.(agent.preference)) )

    nothing
end

function reset!( agent::GradientBandit )
    agent.noOfSteps     = 0
    agent.lastPlayedArm = 0

    agent.avgValue      = 0
    agent.preference    = zeros( Float64, agent.noOfArms )
    agent.pDist         = Categorical( 1/agent.noOfArms*ones(agent.noOfArms) )

    nothing
end

function info_str( agent::GradientBandit, latex::Bool )
    if latex
        return @sprintf( "GradientBandit (\$\\alpha=%4.3f\$)", agent.α )
    else
        return @sprintf( "GradientBandit (α=%4.3f)", agent.α )
    end
end


text:
# Generate Passwords in Julia
using ProgressBars
using Random

# WARNING: Do not use this code to generate actual passwords!
function generate_passwords()
    num_passwords = parse(Int64, Base.prompt("How many passwords do you want to generate?"))
    password_length = parse(Int64, Base.prompt("How long should each password be?"))

    # Create an empty vector / array
    password_holder = []

    # Generate a progress bar to show how close we are to being done
    for i in ProgressBar(1:num_passwords)
        # Add the new password into the password holder
        push!(password_holder, randstring(password_length))
        sleep(0.2) # Manually slowdown the generation of passwords
    end

    # Only show the passwords if there are less than 100
    if length(password_holder) <= 100
        # Loop through each password one by one
        for password in password_holder
            print("\n", password)
        end
    end
end


text:
module SparseTimeSeries
using JSON3

export
    Event,
    Segment,
    EventSeries,
    timestamp,
    value,
    align,
    cumtime,
    fill_forward_event,
    fill_forward_value,
    fuse,
    segments,
    select,
    splice,
    next,
    previous

include("helpers.jl")
include("event_types.jl")
include("segment_types.jl")
include("series_types.jl")

JSON3.StructType(::Type{<:Event}) = JSON3.Struct()
JSON3.StructType(::Type{<:Segment}) = JSON3.Struct()
JSON3.StructType(::Type{<:EventSeries}) = JSON3.Struct()

end # module


text:
# FIXME `ListOfVariableIndices` is not modified by bridges
#       see https://github.com/JuliaOpt/MathOptInterface.jl/issues/693
#    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
#                                                                [[1.0], [0.0], [0.0]],
#                                                                [1.0]),
#                                  (mock) -> MOIU.mock_optimize!(mock,
#                                                                [[1.0], [0.0], [0.0]],
#                                                                [1.0]),
#                                  (mock) -> MOIU.mock_optimize!(mock,
#                                                                [[0.0], [0.0], [1.0], [0.0]],
#                                                                [2.0]),
#                                  (mock) -> MOIU.mock_optimize!(mock,
#                                                                [[0.0], [0.0], [2.0], [0.0]],
#                                                                [2.0]),
#                                  (mock) -> MOIU.mock_optimize!(mock,
#                                                                [[0.0], [1.0], [0.0], [0.0]],
#                                                                [54.4045, 1.0]),
#                                  (mock) -> MOIU.mock_optimize!(mock,
#                                                                [[0.0], [2.0], [0.0]],
#                                                                [61.2652, 1.0]),
#                                  (mock) -> MOIU.mock_optimize!(mock,
#                                                                [[0.0], [0.0], [2.0]],
#                                                                [79.6537, 2.0]),
#                                  (mock) -> MOIU.mock_optimize!(mock,
#                                                                [[0.0], [1.0], [1.0], [0.0]],
#                                                                [30.9277, 1.5, -0.5]))
#    MOIT.linear1test(bridged, config)
MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                            [[1.0], [0.0], [0.0]],
                                                            [1.0]))
MOIT.linear2test(bridged, config)
MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                            [[3.0], [0.0]],
                                                            [-1.0]),
                              (mock) -> MOIU.mock_optimize!(mock,
                                                            [[0.0], [3.0]],
                                                            [-0.0]))
MOIT.linear3test(bridged, config)
MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                            [[0.0], [0.0]],
                                                            [-0.0]),
                              (mock) -> MOIU.mock_optimize!(mock,
                                                            [[0.0], [0.0]],
                                                            [-0.0]),
                              (mock) -> MOIU.mock_optimize!(mock,
                                                            [[0.0], [0.0]],
                                                            [-0.0]))
MOIT.linear4test(bridged, config)
MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                            [[1.3333], [1.3333], [0.0], [0.0]],
                                                            [0.3333, 0.3333]),
                              (mock) -> MOIU.mock_optimize!(mock,
                                                            [[2.0], [0.0], [0.0], [2.0]],
                                                            [0.5, -0.0]),
                              (mock) -> MOIU.mock_optimize!(mock,
                                                            [[4.0], [0.0], [0.0]],
                                                            [1.0]),
                              (mock) -> MOIU.mock_optimize!(mock,
                                                            [[2.0], [0.0]],
                                                            [0.5]))
MOIT.linear5test(bridged, config)
MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                            [[0.0], [0.0], [0.0 0.0; 0.0 0.0], [0.0 0.0; 0.0 0.0]],
                                                            [-1.0, 1.0]),
                              (mock) -> MOIU.mock_optimize!(mock,
                                                            [[0.0], [0.0], [168.953 0.0; 0.0 68.9528], [107.78 0.0; 0.0 107.78]],
                                                            [-1.0, 1.0]),
                              (mock) -> MOIU.mock_optimize!(mock,
                                                            [[0.0], [0.0], [250.152 0.0; 0.0 150.152], [150.152 0.0; 0.0 250.152]],
                                                            [-1.0, 1.0]))
MOIT.linear6test(bridged, config)
MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                            [[0.0], [0.0], [0.0 0.0; 0.0 0.0], [0.0 0.0; 0.0 0.0]],
                                                            [-1.0, 1.0]),
                              (mock) -> MOIU.mock_optimize!(mock,
                                                            [[0.0], [0.0], [168.953 0.0; 0.0 68.9528], [107.78 0.0; 0.0 107.78]],
                                                            [-1.0, 1.0]),
                              (mock) -> MOIU.mock_optimize!(mock,
                                                            [[0.0], [0.0], [250.152 0.0; 0.0 150.152], [150.152 0.0; 0.0 250.152]],
                                                            [-1.0, 1.0]))
MOIT.linear7test(bridged, config)

MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                            MOI.INFEASIBLE,
                                                            tuple(),
                                                            [1.0]))
MOIT.linear8atest(bridged, config)
MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                            MOI.DUAL_INFEASIBLE,
                                                            (MOI.INFEASIBILITY_CERTIFICATE, [[0.7709], [0.2291], [0.3126]])))
MOIT.linear8btest(bridged, config)
MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                            MOI.DUAL_INFEASIBLE,
                                                            (MOI.INFEASIBILITY_CERTIFICATE, [[0.5], [0.5]])))
MOIT.linear8ctest(bridged, config)
MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                            [[29.0909], [36.3636], [4.5455], [0.0], [0.0001]],
                                                            [-0.0, 11.3636, 0.8636]))
MOIT.linear9test(bridged, config)
# FIXME 5 <= x + y <= 10 is bridged into x + y - z == 0, 5 <= z <= 10 and
# then it tries to add two SingleVariable constraints to `z`. We should drop
# support for MOI.LessThan once variables bridges works.
#    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
#                                                                [[5.0], [5.0], [5.0], [0.0]],
#                                                                [-0.0, 1.0]),
#                                  (mock) -> MOIU.mock_optimize!(mock,
#                                                                [[2.5], [2.5], [0.0], [5.0]],
#                                                                [-1.0, 0.0]),
#                                  (mock) -> MOIU.mock_optimize!(mock,
#                                                                [[1.0], [1.0], [0.0], [10.0]],
#                                                                [-1.0, -0.0]),
#                                  (mock) -> MOIU.mock_optimize!(mock,
#                                                                [[6.0], [6.0], [10.0], [0.0]],
#                                                                [-0.0, 1.0]))
#    MOIT.linear10test(bridged, config)
MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                            [[1.0], [0.0], [3.2439 0.0; 0.0 2.2439], [3.2439 0.0; 0.0 2.2439]],
                                                            [0.0, -1.0]),
                              (mock) -> MOIU.mock_optimize!(mock,
                                                            [[0.0], [1.0], [1.3765 0.0; 0.0 0.8765], [1.3765 0.0; 0.0 0.8765]],
                                                            [-1.0, 0.0]))
MOIT.linear11test(bridged, config)
MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                            MOI.INFEASIBLE,
                                                            tuple(),
                                                            [1.0, 3.0]))
MOIT.linear12test(bridged, config)
MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
                                                            [[24.0266], [27.6758 0.0; 0.0 22.6705], [27.6758 0.0; 0.0 22.6705]],
                                                            [0.0, 0.0]))
MOIT.linear13test(bridged, config)
# FIXME z >= 0.0 followed by z <= 1.0. We need to drop support for
#       SingleVariable-in-LessThan but we need variable bridge otherwise,
#       it creates a slack
#    MOIU.set_mock_optimize!(mock, (mock) -> MOIU.mock_optimize!(mock,
#                                                                [[0.0], [0.5], [1.0], [0.0], [0.0]],
#                                                                [2.0, 1.0]),
#                                  (mock) -> MOIU.mock_optimize!(mock,
#                                                                [[1.0], [0.0]],
#                                                                [1.0]))
#    MOIT.linear14test(bridged, config)


text:
module Random

# TODO: rename rand_radial -> rand???

import LinearAlgebra: norm
import ACE: rand_radial, scaling,
            PositionState, ACEBasis,
            rand_sphere, rand_rot, rand_refl, rand_O3,
            B1pComponent

using Random: shuffle

using StaticArrays: @SMatrix, SVector

export rand_nhd, rand_config, rand_sym, randcoeffs, randcombine, rand_vec3

# -------------------------------------------
#   random neighbourhoods and  configurations

rand_radial(Rn::B1pComponent) =
         Rn.meta["rin"] + rand() * (Rn.meta["rcut"] - Rn.meta["rin"])

rand_vec3(Rn::B1pComponent) = rand_radial(Rn) * rand_sphere()

Base.rand(::Type{TX}, basis::ACEBasis) where {TX <: PositionState} =
         TX( (rr = rand_vec3(basis),) )

Base.rand(T::Type{TX}, basis::ACEBasis, N::Integer) where {TX <: PositionState} =
         [ rand(T, basis) for _=1:N ]

function rand_rot(Xs::AbstractVector)
   Q = rand_rot()
   return [ Q * X for X in Xs ]
end

function rand_refl(Xs::AbstractVector)
   σ = rand_refl()
   return [ σ * X for X in Xs ]
end

# --------------------------------------------------------------
# random operations on neighbourhoods, mostly for testing

# TODO: rewrite for generic ACE

# function rand_perm(Rs, Zs)
#    @assert length(Rs) == length(Zs)
#    p = shuffle(1:length(Rs))
#    return Rs[p], Zs[p]
# end
#
#     random potentials / random basis / random potentials ???

# # TODO: we have an issue with eltypes that needs to be fixed!!!
#
# function randcoeffs(basis; diff = 2)
#    ww = scaling(basis, diff)
#    c = 2 * (rand(rfltype(basis), length(basis)) .- 0.5) ./ ww
#    return c / norm(c)
# end
#
# randcombine(basis; diff = 2) =
#    combine(basis, randcoeffs(basis; diff = diff))

# # move to utility???
# function rand(::Type{ACE.RPI.RPIBasis}; kwargs...)
#
# end

end


text:
using DashBootstrapComponents, DashHtmlComponents

card = dbc_card([
    dbc_cardheader(
        dbc_tabs(
            [
                dbc_tab(label="Tab 1", tab_id="tab-1"),
                dbc_tab(label="Tab 2", tab_id="tab-2"),
            ],
            id="card-tabs",
            card=true,
            active_tab="tab-1",
        ),
    ),
    dbc_cardbody(html_p(id="card-content", className="card-text")),
]);

callback!(
    app,
    Output("card-content", "children"),
    Input("card-tabs", "active_tab"),
) do active_tab
    return "This is tab $active_tab"
end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
using Base: isexported

"""
    stdlibobject(name::Symbol)

If `name` is an export of Base or Core, return `x` where `x` is
the object itself. Otherwise, return `nothing`.
"""
function stdlibobject(ex::Symbol)
    if isexported(Base, ex) && isdefined(Base, ex)
        getfield(Base, ex)
    elseif isexported(Core, ex) && isdefined(Core, ex)
        getfield(Core, ex)
    else
        nothing
    end
end

"""
    parsetype(ctx::LintContext, ex::Expr) :: Type

Obtain a supertype of the type represented by `ex`.
"""
function parsetype(ctx::LintContext, ex)
    obj = abstract_eval(ctx, ex)
    if obj !== nothing && isa(obj, Type)
        obj
    elseif isexpr(ex, :curly)
        obj = abstract_eval(ctx, ex.args[1])
        if obj !== nothing && isa(obj, Type) && obj !== Union
            obj
        else
            Any
        end
    else
        Any
    end
end

"""Returns true if ex is `:` or `Colon()`"""
function iscolon(ex)::Bool
    return ex == :(:) || ex == :(Colon())
end

function guesstype(ex::Symbol, ctx::LintContext)
    result = lookup(ctx, ex)
    if result == nothing
        Any  # conservative guess
    else
        result.typeactual
    end
end

# Tag definitions
abstract type ExpressionTag end
struct TupleTag <: ExpressionTag end
struct ColonTag <: ExpressionTag end
struct BlockTag <: ExpressionTag end
struct ReturnTag <: ExpressionTag end
struct MacroCallTag <: ExpressionTag end
struct CallTag <: ExpressionTag end
struct RefTag <: ExpressionTag end
struct ComparisonTag <: ExpressionTag end
struct CurlyTag <: ExpressionTag end
struct TernaryIfTag <: ExpressionTag end
struct LambdaTag <: ExpressionTag end
struct AssignTag <: ExpressionTag end
struct DefaultTag <: ExpressionTag end

# condition → tag
# TODO(felipe): Use Match.jl
function get_tag_per_condition(ex::Expr)
    if isexpr(ex, :tuple)
        return TupleTag
    end
    if isexpr(ex, :(::)) && length(ex.args) == 2
        return ColonTag
    end
    if isexpr(ex, :block)
        return BlockTag
    end
    if isexpr(ex, :return)
        return ReturnTag
    end
    if isexpr(ex, :call)
        return CallTag
    end
    if isexpr(ex, :macrocall)
        return MacroCallTag
    end
    if isexpr(ex, :ref) # it could be a ref a[b] or an array Int[1,2,3], Vector{Int}[]
        return RefTag
    end
    if isexpr(ex, :comparison)
        return ComparisonTag
    end
    if isexpr(ex, :curly)
        return CurlyTag
    end
    # simple if statement e.g. test ? 0 : 1
    if isexpr(ex, :if) && 2 ≤ length(ex.args) ≤ 3
        return TernaryIfTag
    end

    if isexpr(ex, :(->))
        return LambdaTag
    end
    if isexpr(ex, :(=))
        return AssignTag
    end

    return DefaultTag
end

# Specializations for each type below
function guesstype(ex::Expr, ctx::LintContext, ::Type{TupleTag})::Type
    ts = Type[]
    for a in ex.args
        push!(ts, guesstype(a, ctx))
    end
    return Tuple{ts...}
end

function guesstype(ex::Expr, ctx::LintContext, ::Type{ColonTag})::Type
    parsetype(ctx, ex.args[2])
end

function guesstype(ex::Expr, ctx::LintContext, ::Type{BlockTag})::Type
    isempty(ex.args) ? Nothing : guesstype(ex.args[end], ctx)
end

function guesstype(ex::Expr, ctx::LintContext, ::Type{ReturnTag})::Type
    guesstype(ex.args[1], ctx)
end

function guesstype(ex::Expr, ctx::LintContext, ::Type{CallTag})::Type
    fn = ex.args[1]
    if any(x -> isexpr(x, :kw) || isexpr(x, :(...)), ex.args[2:end])
        # TODO: smarter way to deal with kw/vararg
        return Any
    end
    argtypes = map(x -> guesstype(x, ctx), ex.args[2:end])

    # infer return types of Base functions
    obj = abstract_eval(ctx, fn)
    type_argtypes = [isa(t, Type) ? t : Any for t in argtypes]
    if obj !== nothing
        inferred = StaticTypeAnalysis.infertype(obj, type_argtypes)
        if inferred ≠ Any
            return inferred
        end
    end
    return Any
end

function guesstype(ex::Expr, ctx::LintContext, ::Type{MacroCallTag})::Type
    if ex.args[1] == Symbol("@sprintf")
        return String
    elseif ex.args[1] == Symbol("@compat")
        return guesstype(ex.args[2], ctx)
    end
end

function guesstype(ex::Expr, ctx::LintContext, ::Type{RefTag})::Type
        if isexpr(ex.args[1], :curly) # must be a datatype, right?
            elt = abstract_eval(ctx, ex.args[1])
            if elt !== nothing && isa(elt, Type)
                return Vector{elt}
            else
                return Vector
            end
        end

        if isa(ex.args[1], Symbol)
            what = registersymboluse(ex.args[1], ctx)
            if what <: Type
                elt = abstract_eval(ctx, ex.args[1])
                if elt !== nothing && isa(elt, Type)
                    return Vector{elt}
                else
                    return Vector
                end
            end
        end
        # not symbol, or symbol but it refers to a variable
        partyp = guesstype(ex.args[1], ctx)
        if isa(partyp, Symbol)
            # we are in a context of a constructor of a new type, so it's
            # difficult to figure out the content
            return Any
        elseif partyp <: AbstractArray && !(partyp <: AbstractRange)
            eletyp = StaticTypeAnalysis.eltype(partyp)
            nd=0
            try
                nd = ndims(partyp) # This may throw if we couldn't infer the dimension
            catch
                return Any
            end
            dim_diff = nd - (length(ex.args)-1)
            if dim_diff < 0
                if nd == 0 && ex.args[2] == 1 # ok to do A[1] for a 0-dimensional array
                    return eletyp
                else
                    msg(ctx, :E436, ex, "more indices than dimensions")
                    return Any
                end
            end
            colon_arguments = filter(iscolon, ex.args[2:end])
            dim_diff+=length(colon_arguments)

            if dim_diff != 0
                return Array{eletyp, dim_diff} # is this strictly right?
            else
                return eletyp
            end
        else
            argtypes = [guesstype(x, ctx) for x in ex.args]
            type_argtypes = [isa(t, Type) ? t : Any for t in argtypes]
            inferred = StaticTypeAnalysis.infertype(getindex, type_argtypes)
            if ctx.versionreachable(VERSION) && inferred == Union{}
                indtypes = if length(type_argtypes) == 1
                    "no indices"
                else
                    string("index types ", join(type_argtypes[2:end], ", "))
                end
                msg(ctx, :E522, ex,
                    string("indexing $(type_argtypes[1]) with ",
                           indtypes,
                           " is not supported"))
            end
            return inferred
        end
end

function guesstype(ex::Expr, ctx::LintContext, ::Type{ComparisonTag})::Type
    Bool
end

function guesstype(ex::Expr, ctx::LintContext, ::Type{CurlyTag})::Type
    Type
end

function guesstype(ex::Expr, ctx::LintContext, ::Type{TernaryIfTag})::Type
    tt = guesstype(ex.args[2], ctx)
    ft = if length(ex.args) == 3
        guesstype(ex.args[3], ctx)
    else
        Nothing
    end
    if tt == ft
        # we need this case because tt and ft might be symbols
        return tt
    elseif isa(tt, Type) && isa(ft, Type)
        return Union{tt, ft}
    else
        return Any
    end
end

function guesstype(ex::Expr, ctx::LintContext, ::Type{LambdaTag})::Type
    Function
end

function guesstype(ex::Expr, ctx::LintContext, ::Type{DefaultTag})::Type
    Any
end

function guesstype(ex::Expr, ctx::LintContext, ::Type{AssignTag})::Type
    return guesstype(ex.args[2], ctx)
end

# `guesstype` dispatcher (resolves the tag)
function guesstype(ex::Expr, ctx::LintContext)::Type
    ex = ExpressionUtils.expand_trivial_calls(ex)
    return guesstype(ex, ctx, get_tag_per_condition(ex))
end

guesstype(ex, ctx::LintContext) = lexicaltypeof(ex)


text:
export iauEpb2jd
"""
Besselian Epoch to Julian Date.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   epb      double    Besselian Epoch (e.g. 1957.3)

Returned:
   djm0     double    MJD zero-point: always 2400000.5
   djm      double    Modified Julian Date

Note:

   The Julian Date is returned in two pieces, in the usual SOFA
   manner, which is designed to preserve time resolution.  The
   Julian Date is available as a single number by adding djm0 and
   djm.

Reference:

   Lieske, J.H., 1979, Astron.Astrophys. 73, 282.

This revision:  2013 August 13

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauEpb2jd(double epb, double *djm0, double *djm)

function iauEpb2jd(epb::Real)
   ref_djm0 = Ref{Float64}(0.0)
   ref_djm  = Ref{Float64}(0.0)

   ccall((:iauEpb2jd, libsofa_c), Cvoid,
         (Cdouble, Ref{Cdouble}, Ref{Cdouble}),
         convert(Float64, epb),
         ref_djm0, ref_djm)

   return ref_djm0[], ref_djm[]
end


text:
function add_parameters(ps_m::CanonicalModel,
                        data::Matrix,
                        param_name::Symbol,
                        axs...)

    ps_m.parameters[param_name] = JuMP.Containers.DenseAxisArray{PJ.ParameterRef}(undef, axs...)

    Cidx = CartesianIndices(length.(axs))

    for idx in Cidx
        ps_m.parameters[param_name].data[idx] = PJ.add_parameter(ps_m.JuMPmodel,data[idx])
    end

    return

end

function add_parameters(ps_m::CanonicalModel,
                        ts_data::Vector{Tuple{String,Int64, Vector{Float64}}},
                        time_range::UnitRange{Int64},
                        param_name::Symbol,
                        expression::Symbol)

    ps_m.parameters[param_name] = JuMP.Containers.DenseAxisArray{PJ.ParameterRef}(undef, [r[1] for r in ts_data], time_range)

    for t in time_range, r in ts_data
        ps_m.parameters[param_name][r[1], t] = PJ.add_parameter(ps_m.JuMPmodel, r[3][t]);
        _add_to_expression!(ps_m.expressions[expression], r[2], t, ps_m.parameters[param_name][r[1], t])
    end

    return

end


text:
using Test: @test, @testset, @inferred, @test_throws
import RestrictedBoltzmannMachines as RBMs

@testset "obs" begin
    X, Y = randn(10, 4), randn(5, 4)
    @test (@inferred RBMs._nobs(X, Y)) == 4
    @test_throws AssertionError RBMs._nobs(X, randn(5, 3))
    @test (@inferred RBMs._getobs([2, 3], X, Y)) == (X[:, [2, 3]], Y[:, [2, 3]])

    X, Y, Z = randn(10, 4), nothing, randn(5, 4)
    @test (@inferred RBMs._nobs(X, Y, Z)) == 4
    @test_throws AssertionError RBMs._nobs(X, randn(5, 3), nothing)
    @test (@inferred RBMs._getobs([2, 3], X, Y, Z)) == (X[:, [2, 3]], nothing, Z[:, [2, 3]])
end

@testset "minibatches" begin
    @test RBMs.minibatch_count(10; batchsize = 3) == 4
    @test RBMs.minibatches(7; batchsize=3, shuffle=false) == [[1,2,3], [4,5,6], [7,1,2]]
    @test RBMs.minibatches(7; batchsize=5, shuffle=false) == [[1,2,3,4,5], [6,7,1,2,3]]
    @test length(RBMs.minibatches(7; batchsize=3)) == RBMs.minibatch_count(7; batchsize=3)

    X, Y, Z = randn(10, 4), nothing, randn(5, 4)
    batches = RBMs.minibatches(X, Y, Z; batchsize=2, shuffle=false)
    @test length(batches) == RBMs.minibatch_count(4; batchsize=2) == 2
    @test length(batches) == RBMs.minibatch_count(X, Y, Z; batchsize=2)
    @test batches[1] == (X[:,1:2], nothing, Z[:,1:2])
    @test batches[2] == (X[:,3:4], nothing, Z[:,3:4])

    X = randn(7, 10)
    batches = RBMs.minibatches(X, X, copy(X); batchsize=3, shuffle=true)
    for (x1, x2, x3) in batches
        @test size(x1) == (7, 3)
        @test x1 == x2 == x3
    end

    batch = only(RBMs.minibatches([1,2,3], [1,2,3]'; batchsize=3))
    p = sortperm(batch[1])
    @test batch[1][p] == batch[2]'[p] == [1,2,3]
end

@testset "training_epochs" begin
    nsamples = 5421
    nupdates = 10403
    batchsize = 5
    nepochs = RBMs.training_epochs(; nsamples, nupdates, batchsize)
    @test RBMs.minibatch_count(nsamples; batchsize) * (nepochs - 1) ≤ nupdates
    @test RBMs.minibatch_count(nsamples; batchsize) * (nepochs + 1) ≥ nupdates
end


text:
module HexapoleDataFit

    include("types.jl")
    include("ellipse.jl")
    include("IO.jl")
    include("ellipse_fitting.jl")
    include("plot.jl")
    include("transformations.jl")
    include("kernels.jl")

end # module


text:
# example taken from Aragón, J. L., R. A. Barrio, T. E. Woolley, R. E. Baker, and P. K. Maini. “Nonlinear Effects on Turing Patterns: Time Oscillations and Chaos.” Physical Review E 86, no. 2 (August 8, 2012): 026201. https://doi.org/10.1103/PhysRevE.86.026201.
using Revise
	using DiffEqOperators, ForwardDiff, DifferentialEquations, Sundials
	using BifurcationKit, LinearAlgebra, Plots, SparseArrays, Parameters, Setfield
	const BK = BifurcationKit

norminf(x) = norm(x, Inf)
f(u, v, p) = p.η * (      u + p.a * v - p.C * u * v - u * v^2)
g(u, v, p) = p.η * (p.H * u + p.b * v + p.C * u * v + u * v^2)

function Laplacian(N, lx, bc = :Dirichlet)
	hx = 2lx/N
	D2x = CenteredDifference(2, 2, hx, N)
	if bc == :Neumann
		Qx = Neumann0BC(hx)
	elseif bc == :Dirichlet
		Qx = Dirichlet0BC(typeof(hx))
	elseif bc == :Periodic
		Qx = PeriodicBC(typeof(hx))
	end
	D2xsp = sparse(D2x * Qx)[1] |> sparse
end

@views function NL!(dest, u, p, t = 0.)
	N = div(length(u), 2)
	u1 =  u[1:N]
	u2 =  u[N+1:end]

	dest[1:N]     .= f.(u1, u2, Ref(p))
	dest[N+1:end] .= g.(u1, u2, Ref(p))

	return dest
end

function Fbr!(f, u, p, t = 0.)
	mul!(f, p.Δ, u)
	f .= f .+ NL(u, p)
end

NL(u, p) = NL!(similar(u), u, p)
Fbr(x, p, t = 0.) = Fbr!(similar(x), x, p)
dNL(x, p, dx) = ForwardDiff.derivative(t -> NL(x .+ t .* dx, p), 0.)

function dFbr(x, p, dx)
	f = similar(dx)
	mul!(f, p.Δ, dx)

	nl = dNL(x, p, dx)
	f .= f .+ nl
end

Jbr(x, p) = sparse(ForwardDiff.jacobian(x -> Fbr(x, p), x))
jet = BK.getJet(Fbr, Jbr)

####################################################################################################
N = 100
	n = 2N
	lx = 3pi /2
	X = LinRange(-lx,lx, N)
	Δ = Laplacian(N, lx, :Neumann)
	D = 0.08
	par_br = (η = 1.0, a = -1., b = -3/2., H = 3.0, D = D, C = -0.6, Δ = blockdiag(D*Δ, Δ))

	u0 = cos.(2X)
	solc0 = vcat(u0, u0)
####################################################################################################
# eigls = DefaultEig()
eigls = EigArpack(0.5, :LM)
optnewton = NewtonPar(eigsolver = eigls, verbose=true, maxIter = 3200, tol=1e-9)

out, = @time newton(Fbr, Jbr, solc0, par_br, optnewton, normN = norminf)
	plot();plot!(X,out[1:N]);plot!(X,solc0[1:N], label = "sol0",line=:dash)

optcont = ContinuationPar(dsmax = 0.051, ds = -0.001, pMin = -1.8, detectBifurcation = 3, nev = 21, plotEveryStep = 50, newtonOptions = optnewton, maxSteps = 370, nInversion = 10, maxBisectionSteps = 25)

	br, = @time continuation(Fbr, Jbr, solc0, (@set par_br.C = -0.2), (@lens _.C), (@set optcont.newtonOptions.verbose = false);
		plot = true, verbosity = 2,
		recordFromSolution = (x, p) -> norm(x, Inf),
		plotSolution = (x, p; kwargs...) -> plot!(x[1:end÷2];label="",ylabel ="u", kwargs...))

plot(br)
####################################################################################################
# branching from Hopf bp using aBS-Trapezoid
opt_po = NewtonPar(tol = 1e-10, verbose = true, maxIter = 120)

eig = EigKrylovKit(tol= 1e-10, x₀ = rand(2N), verbose = 2, dim = 40)
eig = DefaultEig()
optcontpo = ContinuationPar(dsmin = 0.001, dsmax = 0.015, ds= 0.01, pMin = -1.8, maxSteps = 60, newtonOptions = (@set opt_po.eigsolver = eig), nev = 25, precisionStability = 1e-7, detectBifurcation = 3, dsminBisection = 1e-6)
M = 100
	probFD = PeriodicOrbitTrapProblem(M = M)
	br_po, = @time continuation(
		# arguments for branch switching from the first
		# Hopf bifurcation point
		jet..., br, 1,
		# arguments for continuation
		optcontpo, probFD;
		# OPTIONAL parameters
		# we want to jump on the new branch at phopf + δp
		# ampfactor is a factor to increase the amplitude of the guess
		verbosity = 2,
		plot = true,
		jacobianPO = :FullSparseInplace,
		plotSolution = (x, p;kwargs...) ->  (heatmap!(reshape(x[1:end-1], 2*N, M)'; ylabel="time", color=:viridis, kwargs...);plot!(br, subplot=1)),
		recordFromSolution = (u, p) -> maximum(u[1:end-1]),#BK.maximumPOTrap(u, N, M; ratio = 2),
		normC = norminf)

plot(br, br_po, label = "")
####################################################################################################
# branching from PD using aBS
br_po_pd, = @time BK.continuation(
		# arguments for branch switching from the first
		# Hopf bifurcation point
		br_po, 1, setproperties(br_po.contparams; detectBifurcation = 3, plotEveryStep = 1, ds = 0.01);
		# OPTIONAL parameters
		# we want to jump on the new branch at phopf + δp
		# ampfactor is a factor to increase the amplitude of the guess
		ampfactor = 0.9, δp = -0.01,
		verbosity = 3,
		plot = true,
		# jacobianPO = :FullSparseInplace,
		jacobianPO = :BorderedSparseInplace,
		plotSolution = (x, p;kwargs...) ->  (heatmap!(reshape(x[1:end-1], 2*N, M)'; ylabel="time", color=:viridis, kwargs...);plot!(br_po, subplot=1)),
		recordFromSolution = (u, p) -> maximum(u[1:end-1]),#BK.maximumPOTrap(u, N, M; ratio = 2),
		normC = norminf)

plot(br, br_po, br_po_pd, label = "")
####################################################################################################
# Period doubling
ind_pd = 1
vec_pd = geteigenvector(eig,
		br_po.eig[br_po.specialpoint[ind_pd].idx][2],
		br_po.specialpoint[ind_pd].ind_ev)

orbitguess_f = br_po.specialpoint[1].x .+ 0.1 * real.(vec_pd)
orbitguess_f[end] *= 2
#
# outpo_pd, _, _ = @time BK.newton(
# 		poTrap(br_po.bifpoint[1].param),
# 		orbitguess_f, opt_po, :FullLU;
# 		normN = norminf)
# 	printstyled(color=:red, "--> T = ", outpo_pd[end], ", amplitude = ", BK.amplitude(outpo_pd, 2N, M),"\n")
# 	BK.plotPeriodicPOTrap(outpo_pd, N, M)
####################################################################################################
# shooting
par_br_hopf = @set par_br.C = -0.86
f1 = DiffEqArrayOperator(par_br.Δ)
f2 = NL!
prob_sp = SplitODEProblem(f1, f2, solc0, (0.0, 280.0), par_br_hopf)

sol = @time solve(prob_sp, ETDRK2(krylov=true); abstol=1e-14, reltol=1e-14, dt = 0.1, progress = true)

prob_ode = ODEProblem(Fbr, solc0, (0.0, 280.0), par_br_hopf)
sol = @time solve(prob_ode, Rodas4P2(); abstol=1e-14, reltol=1e-7, dt = 0.1, progress = true)
orbitsection = Array(sol[:,[end]])
# orbitsection = orbitguess[:, 1]

initpo = vcat(vec(orbitsection), 3.)

BK.plotPeriodicShooting(initpo[1:end-1], 1);title!("")

probSh = ShootingProblem(prob_sp, ETDRK2(krylov=true), [sol(280.0)]; abstol=1e-14, reltol=1e-14, dt = 0.1, parallel = true)
# probSh = ShootingProblem(prob_ode, Rodas4P(), [sol(280.0)]; abstol=1e-10, reltol=1e-4, parallel = true)

probSh(initpo, par_br_hopf)

ls = GMRESIterativeSolvers(reltol = 1e-7, N = length(initpo), maxiter = 50, verbose = false)
	# ls = GMRESKrylovKit(verbose = 0, dim = 200, atol = 1e-9, rtol = 1e-5)
	optn = NewtonPar(verbose = true, tol = 1e-9,  maxIter = 20, linsolver = ls)
	# deflationOp = BK.DeflationOperator(2 (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [outpo])
	outposh, _, flag = @time newton(probSh, initpo, par_br_hopf, optn;
		callbackN = (x, f, J, res, iteration; kw...) -> (@show x[end];true),
		normN = norminf)
	flag && printstyled(color=:red, "--> T = ", outposh[end], ", amplitude = ", BK.getAmplitude(probSh, outposh, par_br_hopf; ratio = 2),"\n")

plot(initpo[1:end-1], label = "Init guess")
	plot!(outposh[1:end-1], label = "sol")

eig = EigKrylovKit(tol= 1e-12, x₀ = rand(2N), verbose = 2, dim = 40)
# eig = DefaultEig()
optcontpo = ContinuationPar(dsmin = 0.0001, dsmax = 0.01, ds= -0.005, pMin = -1.8, maxSteps = 50, newtonOptions = (@set optn.eigsolver = eig), nev = 10, precisionStability = 1e-2, detectBifurcation = 3)
	br_po_sh, = @time continuation(probSh, outposh, par_br_hopf, (@lens _.C), optcontpo;
		verbosity = 3,	plot = true,
		linearAlgo = MatrixFreeBLS(@set ls.N = probSh.M*n+2),
		finaliseSolution = (z, tau, step, contResult; kw...) ->
			(BK.haseigenvalues(contResult) && Base.display(contResult.eig[end].eigenvals) ;true),
		plotSolution = (x, p; kwargs...) -> BK.plotPeriodicShooting!(x[1:end-1], 1; kwargs...),
		recordFromSolution = (u, p) -> BK.getMaximum(probSh, u, (@set par_br_hopf.C = p.p); ratio = 2), normC = norminf)

# branches = [br_po_sh]
# push!(branches, br_po_sh)
# plot(branches...)

plot(br_po_sh, br, label = "")

####################################################################################################
# shooting Period Doubling
par_br_pd = @set par_br.C = -1.32
f1 = DiffEqArrayOperator(par_br.Δ)
f2 = NL!
prob_sp = SplitODEProblem(f1, f2, solc0, (0.0, 300.0), par_br_pd; abstol=1e-14, reltol=1e-14, dt = 0.1)
# solution close to the PD point.

solpd = @time solve(prob_sp, ETDRK2(krylov=true))
	# heatmap(sol.t, X, sol[1:N,:], color=:viridis, xlim=(20,280.0))

plot(solpd.t, solpd[N÷2, :], xlim=(290,296.2))

solpd.t[end-100]

orbitsectionpd = Array(solpd[:,end-100])
initpo_pd = vcat(vec(orbitsectionpd), 6.2)
BK.plotPeriodicShooting(initpo_pd[1:end-1], 1);title!("")

ls = GMRESIterativeSolvers(reltol = 1e-7, N = length(initpo_pd), maxiter = 50, verbose = false)
	# ls = GMRESKrylovKit(verbose = 0, dim = 200, atol = 1e-9, rtol = 1e-5)
	optn = NewtonPar(verbose = true, tol = 1e-9,  maxIter = 120, linsolver = ls)
	# deflationOp = BK.DeflationOperator(2 (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [outpo])
	outposh_pd, _, flag = @time newton(probSh, initpo_pd, par_br_pd, optn;
		callback = (x, f, J, res, iteration, itlinear, options; kwargs...) -> (@show x[end];true),
		normN = norminf)
	flag && printstyled(color=:red, "--> T = ", outposh_pd[end], ", amplitude = ", BK.getAmplitude(probSh, outposh_pd, (@set par_br.C = -0.86); ratio = 2),"\n")

	plot(initpo[1:end-1], label = "Init guess")
	plot!(outposh_pd[1:end-1], label = "sol")

optcontpo = ContinuationPar(dsmin = 0.0001, dsmax = 0.005, ds= -0.001, pMin = -1.8, maxSteps = 500, newtonOptions = (@set optn.eigsolver = eig), nev = 10, precisionStability = 1e-2, detectBifurcation = 0)
	br_po_sh_pd, = @time continuation(
		probSh, outposh_pd, par_br_pd, (@lens _.C),
		optcontpo; verbosity = 3,
		plot = true,
		linearAlgo = MatrixFreeBLS(@set ls.N = probSh.M*n+2),
		finaliseSolution = (z, tau, step, contResult; k...) ->
			(Base.display(contResult.eig[end].eigenvals) ;println("--> T = ", z.u[end]);true),
		plotSolution = (x, p; kwargs...) -> (BK.plotPeriodicShooting!(x[1:end-1], 1; kwargs...); plot!(br_po_sh; subplot=1, legend=false)),
		recordFromSolution = (u, p; k...) -> BK.getMaximum(probSh, u, (@set par_br_pd.C = p); ratio = 2), normC = norminf)

plot(br_po_sh_pd, br, label = "");title!("")
####################################################################################################
# branching from Hopf bp using aBS- Shooting
ls = GMRESIterativeSolvers(reltol = 1e-7, maxiter = 50, verbose = false)
eig = EigKrylovKit(tol= 1e-10, x₀ = rand(2N), verbose = 2, dim = 40)
eig = DefaultEig()

opt_po = NewtonPar(tol = 1e-10, verbose = true, maxIter = 120, linsolver  = ls)
optcontpo = ContinuationPar(dsmin = 0.0001, dsmax = 0.01, ds= -0.005, pMin = -1.8, maxSteps = 50, newtonOptions = (@set opt_po.eigsolver = eig), nev = 20, precisionStability = 1e-2, detectBifurcation = 3)

probPO = ShootingProblem(1, prob_sp, ETDRK2(krylov=true); abstol=1e-14, reltol=1e-14, dt = 0.1)
	br_po, = @time continuation(
		# arguments for branch switching from the first
		# Hopf bifurcation point
		jet..., br, 1,
		# arguments for continuation
		optcontpo, probPO;
		# OPTIONAL parameters
		# we want to jump on the new branch at phopf + δp
		# ampfactor is a factor to increase the amplitude of the guess
		δp = 0.005,
		verbosity = 3,
		plot = true,
		linearAlgo = MatrixFreeBLS(@set ls.N = probPO.M*n+2),
		finaliseSolution = (z, tau, step, contResult; kw...) ->
			(BK.haseigenvalues(contResult) && Base.display(contResult.eig[end].eigenvals) ;true),
		plotSolution = (x, p; kwargs...) -> (BK.plotPeriodicShooting!(x[1:end-1], 1; kwargs...);plot!(br, subplot=1)),
		recordFromSolution = (u, p) -> BK.getMaximum(probPO, u, (@set par_br.C = p.p); ratio = 2), normC = norminf)

plot(br, br_po, label = "")

# aBS from PD
@set! br_po.contparams.newtonOptions.tol = 1e-7
br_po_pd, = BK.continuation(br_po, 1, setproperties(br_po.contparams, detectBifurcation = 0, maxSteps = 5, ds = 0.01, plotEveryStep = 1);
	verbosity = 3, plot = true,
	ampfactor = .2, δp = -0.01,
	usedeflation = false,
	# for aBS from period doubling, we double the sections
	linearAlgo = MatrixFreeBLS(@set ls.N = 2probPO.M*n+2),
	plotSolution = (x, p; kwargs...) -> begin
		outt = BK.getPeriodicOrbit(p.prob, x, (@set par_br_hopf.C = p.p))
		heatmap!(outt[:,:]'; color = :viridis, subplot = 3)
		plot!(br_po; legend=false, subplot=1)
	end,
	recordFromSolution = (u, p) -> (BK.getMaximum(p.prob, u, (@set par_br_hopf.C = p.p); ratio = 2)), normC = norminf
	)

plot(br_po, br_po_pd, legend=false)
####################################################################################################
# aBS Poincare Shooting

br_po.contparams.newtonOptions.linsolver.solver.N
br_po_pd.contparams.newtonOptions.linsolver.solver.N

####################################################################################################