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,
}

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

mutable struct VariableInfo{T<:Union{Float64, Vector{Float64}}}
    upper::T
    lower::T
    function VariableInfo()
        new{Float64}(NaN, NaN)
    end
    function VariableInfo(N::Integer)
        new{Vector{Float64}}(fill(NaN, N), fill(NaN, N))
    end
end
function VariableInfo(_info::ConstraintInfo{T}) where T
    if isa(_info.upper, Number)
        info = VariableInfo()
        info.lower = _info.lower
        info.upper = _info.upper
    else
        info = VariableInfo(length(_info.upper))
        info.lower .= _info.lower
        info.upper .= _info.upper
    end
    return info
end

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

abstract type AbstractBilevelSolverMode{T} end

mutable struct NoMode{T} <: AbstractBilevelSolverMode{T}
end

mutable struct SOS1Mode{T} <: AbstractBilevelSolverMode{T}
    function SOS1Mode()
        return new{Float64}()
    end
end

mutable struct ComplementMode{T} <: AbstractBilevelSolverMode{T}
    with_slack::Bool
    function ComplementMode(;with_slack = false)
        return new{Float64}(with_slack)
    end
end

mutable struct ProductMode{T} <: AbstractBilevelSolverMode{T}
    epsilon::T
    with_slack::Bool
    function ProductMode(eps::T=zero(Float64); with_slack = false) where T
        return new{Float64}(eps, with_slack)
    end
end

@enum IndicatorSetting ZERO_ONE ZERO_ZERO ONE_ONE

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

struct ComplementBoundCache
    # internal usage
    upper::Dict{VI, VariableInfo}
    lower::Dict{VI, VariableInfo}
    ldual::Dict{CI, ConstraintInfo}
    # full map
    map::Dict{VI, VariableInfo}
    function ComplementBoundCache()
        return new(
            Dict{VI, VariableInfo}(),
            Dict{VI, VariableInfo}(),
            Dict{CI, ConstraintInfo}(),
            Dict{VI, VariableInfo}(),
        )
    end
end

abstract type AbstractBoundedMode{T} <: AbstractBilevelSolverMode{T} end

mutable struct FortunyAmatMcCarlMode{T} <: AbstractBoundedMode{T}
    with_slack::Bool
    primal_big_M::Float64
    dual_big_M::Float64
    cache::ComplementBoundCache
    function FortunyAmatMcCarlMode(;with_slack = false,
        primal_big_M = Inf, dual_big_M = Inf)
        return new{Float64}(
            with_slack,
            primal_big_M,
            dual_big_M,
            ComplementBoundCache()
        )
    end
end

mutable struct MixedMode{T} <: AbstractBoundedMode{T}
    default::AbstractBilevelSolverMode{T}
    constraint_mode_map_c::Dict{CI, AbstractBilevelSolverMode{T}}
    constraint_mode_map_v::Dict{VI, AbstractBilevelSolverMode{T}}
    cache::ComplementBoundCache
    function MixedMode(;default = SOS1Mode())
        return new{Float64}(
            default,
            Dict{CI, AbstractBilevelSolverMode{Float64}}(),
            Dict{VI, AbstractBilevelSolverMode{Float64}}(),
            ComplementBoundCache()
        )
    end
end

function reset!(::AbstractBilevelSolverMode)
    return nothing
end
function reset!(mode::AbstractBoundedMode)
    mode.cache = ComplementBoundCache()
    return nothing
end

function appush!(col, element::AbstractVector)
    append!(col, element)
    return nothing
end
function appush!(col, element)
    push!(col, element)
    return nothing
end

function build_full_map!(mode,
    upper_idxmap, lower_idxmap, lower_dual_idxmap, lower_primal_dual_map)
    return nothing
end
function build_full_map!(mode::AbstractBoundedMode,
        upper_idxmap, lower_idxmap, lower_dual_idxmap, lower_primal_dual_map)
    _build_bound_map!(mode.cache,
        upper_idxmap, lower_idxmap, lower_dual_idxmap, lower_primal_dual_map)
    return nothing
end
function _build_bound_map!(mode::ComplementBoundCache,
    upper_idxmap, lower_idxmap, lower_dual_idxmap, lower_primal_dual_map)
    empty!(mode.map)
    for (k,v) in mode.upper
        mode.map[upper_idxmap[k]] = v
    end
    for (k,v) in mode.lower
        mode.map[lower_idxmap[k]] = v
    end
    for (k,v) in mode.ldual
        vec = lower_primal_dual_map.primal_con_dual_var[k]#[1] # TODO check this scalar
        for var in vec
            mode.map[lower_dual_idxmap[var]] = VariableInfo(v)
        end
    end
    return nothing
end

mutable struct StrongDualityMode{T} <: AbstractBilevelSolverMode{T}
    inequality::Bool
    epsilon::T
    function StrongDualityMode(eps::T=zero(Float64); inequality = true) where T
        return new{T}(inequality, eps)
    end
end

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

function accept_vector_set(mode::AbstractBilevelSolverMode{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(::MixedMode{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.copy(MOI.get(primal_model, MOI.ConstraintFunction(), ci))::F
    set = MOI.copy(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, ci, 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.copy(MOI.get(primal_model, MOI.ConstraintFunction(), ci))::F
    set = MOI.copy(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, ci, func, set_with_zero(set), map[ci][1])
    return con
end

function set_with_zero(set::S) where {S<:SCALAR_SETS}
    return S(0.0)
end
function set_with_zero(set)
    return MOI.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}();
    copy_names::Bool = false,
    pass_start::Bool = false
    )

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

    #=
        Initialize Lower DUAL level model
    =#
    # dualize the second level
    dual_problem = Dualization.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)
    if copy_names
        pass_names(m, upper, upper_idxmap)
    end

    #=
        Pass Lower level model
    =#

    handle_lower_objective_sense(lower)

    # 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}())
        # deepcopy and 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)
    if copy_names
        pass_names(m, lower, lower_idxmap)
    end

    #=
        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)
    if copy_names
        pass_names(m, lower_dual, lower_dual_idxmap)
    end

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

    # build map bound map for FortunyAmatMcCarlMode
    build_full_map!(mode,
    upper_idxmap, lower_idxmap, lower_dual_idxmap, lower_primal_dual_map)

    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, copy_names, pass_start)
            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::StrongDualityMode{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)

    if !mode.inequality
        c = MOIU.normalize_and_add_constraint(m, func, MOI.EqualTo(zero(T)))
        MOI.set(m, MOI.ConstraintName(), c, "lower_strong_duality")
        return CI[c]
    else
        func_up = MOIU.operate(-, T, func, mode.epsilon)
        c_up = MOIU.normalize_and_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 = MOIU.normalize_and_add_constraint(m, func_lo, MOI.GreaterThan(zero(T)))
        MOI.set(m, MOI.ConstraintName(), c_lo, "lower_strong_duality_lo")

        return CI[c_up, c_lo]
    end
end

function add_complement(mode::MixedMode{T}, m, comp::Complement,
        idxmap_primal, idxmap_dual, copy_names::Bool, pass_start::Bool) where T
    _mode = get_mode(mode, comp.constraint, idxmap_primal)
    accept_vector_set(_mode, comp)
    add_complement(_mode, m, comp,
        idxmap_primal, idxmap_dual, copy_names, pass_start)
end

function get_mode(mode::MixedMode{T}, ci::CI{F,S}, map) where {
    T,
    F<:MOI.SingleVariable,
    S<:Union{MOI.EqualTo{T}, MOI.LessThan{T}, MOI.GreaterThan{T}}
}
    # key = map[VI(ci.value)]
    key = VI(ci.value)
    if haskey(mode.constraint_mode_map_v, key)
        return mode.constraint_mode_map_v[key]
    else
        return mode.default
    end
end
function get_mode(mode::MixedMode{T}, ci::CI{F,S}, map) where T where {F,S}
    key = ci#map[ci]
    # @show ci, key, map
    # @show mode.constraint_mode_map_c
    if haskey(mode.constraint_mode_map_c, key)
        return mode.constraint_mode_map_c[key]
    else
        return mode.default
    end
end

function add_complement(mode::ComplementMode{T}, m, comp::Complement,
    idxmap_primal, idxmap_dual, copy_names::Bool, pass_start::Bool) where T
    f = comp.func_w_cte
    s = comp.set_w_zero
    v = comp.variable
    out_var = VI[]
    out_ctr = CI[]

    with_slack = mode.with_slack

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

    dual = idxmap_dual[v]

    if with_slack
        slack, slack_in_set = MOI.add_constrained_variable(m, s)
        new_f = MOIU.operate(-, T, f_dest, MOI.SingleVariable(slack))
        equality = MOIU.normalize_and_add_constraint(m, new_f, MOI.EqualTo(zero(T)))

        if pass_start
            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
        end

        c = MOI.add_constraint(m,
            MOI.VectorOfVariables([slack, dual]),
            MOI.Complements(1))
        if copy_names
            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(), c, "compl_complWslk_($(nm))")
        end

        appush!(out_var, slack)
        appush!(out_ctr, slack_in_set)
        appush!(out_ctr, equality)
        appush!(out_ctr, c)
    else
        new_f = MOIU.operate(vcat, T, f_dest, MOI.SingleVariable(dual))

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

        if copy_names
            nm = MOI.get(m, MOI.VariableName(), dual)
            MOI.set(m, MOI.ConstraintName(), c, "compl_compl_($(nm))")
        end

        appush!(out_ctr, c)
    end

    return out_var, out_ctr
end

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

    if comp.is_vec
        error("Vector constraint is not supported by SOS1 mode")
    end

    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]))

    if copy_names
        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))")
    end

    return slack, slack_in_set, equality, c1
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)

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

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

    out_var = VI[]
    out_ctr = CI[]

    eps = mode.epsilon
    with_slack = mode.with_slack

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

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

    if with_slack
        slack, slack_in_set = if comp.is_vec
            MOI.add_constrained_variables(m, s)
        else
            MOI.add_constrained_variable(m, s)
        end
        new_f = MOIU.operate(-, T, f_dest, only_variable_functions(slack))
        if comp.is_vec
            equality = MOIU.normalize_and_add_constraint(m, new_f, MOI.Zeros(length(slack)))
        else
            equality = MOIU.normalize_and_add_constraint(m, new_f, MOI.EqualTo(zero(T)))
        end

        prod_f = MOIU.operate(dot, T, only_variable_functions(slack), only_variable_functions(dual))

        prod_f1 = MOIU.operate(-, T, prod_f, eps)
        c1 = MOIU.normalize_and_add_constraint(m,
            prod_f1,
            MOI.LessThan{Float64}(0.0))
        if comp.is_vec
            prod_f2 = MOIU.operate(+, T, prod_f, eps)
            c2 = MOIU.normalize_and_add_constraint(m,
                prod_f2,
                MOI.GreaterThan{Float64}(0.0))
        end

        appush!(out_var, slack)
        appush!(out_ctr, slack_in_set)
        appush!(out_ctr, equality)
        appush!(out_ctr, c1)
        if comp.is_vec
            appush!(out_ctr, c2)
        end

        if pass_start
            val = MOIU.eval_variables(
                x-> nothing_to_nan(MOI.get(m, MOI.VariablePrimalStart(), x)), f_dest)
            if comp.is_vec
                for i in eachindex(val)
                    if !isnan(val[i])
                        MOI.set(m, MOI.VariablePrimalStart(), slack[i], val[i])
                    end
                end
            else
                if !isnan(val)
                    MOI.set(m, MOI.VariablePrimalStart(), slack, val)
                end
            end
        end

        if copy_names
            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))")
            if comp.is_vec
                MOI.set(m, MOI.ConstraintName(), c1, "compl_prodWslk2_($(nm))")
            end
        end
    else
        new_f = MOIU.operate(dot, T, f_dest, only_variable_functions(dual))
        new_f1 = MOIU.operate(-, T, new_f, eps)
        c1 = MOIU.normalize_and_add_constraint(m,
            new_f1,
            MOI.LessThan{T}(0.0))
        if comp.is_vec # conic
            new_f2 = MOIU.operate(+, T, new_f, eps)
            c2 = MOIU.normalize_and_add_constraint(m,
                new_f2,
                MOI.GreaterThan{T}(0.0))
        end

        # TODO(?): if eps == 0 then add equality

        appush!(out_ctr, c1)
        if comp.is_vec
            appush!(out_ctr, c2)
        end

        if copy_names
            nm = if comp.is_vec
                MOI.get.(m, MOI.VariableName(), dual)
            else
                MOI.get(m, MOI.VariableName(), dual)
            end
            MOI.set(m, MOI.ConstraintName(), c1, "compl_prod_($(nm))")
            if comp.is_vec
                MOI.set(m, MOI.ConstraintName(), c2, "compl_prod2_($(nm))")
            end
        end
    end
    return out_var, out_ctr
end

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

    method = mode.mode

    is_tight = false
    has_start = false

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

    dual = idxmap_dual[v]

    if comp.is_vec
        error("Vector constraint is (currently) not supported by indicator mode")
    end

    if copy_names
        nm = if comp.is_vec
            MOI.get.(m, MOI.VariableName(), dual)
        else
            MOI.get(m, MOI.VariableName(), dual)
        end
    end

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

    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())
        if copy_names
            MOI.set(m, MOI.VariableName(), vb2, "compl_bin2_($(nm))")
        end

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

    pre_f1, pre_s1 = MOIU.normalize_constant(f_dest, MOI.EqualTo(zero(T)))
    f1 = MOIU.operate(vcat, T, SVF(vb1), pre_f1)
    f2 = MOIU.operate(vcat, T, SVF(vb2), SVF(dual))

    if pass_start
        val = MOIU.eval_variables(
            x-> nothing_to_nan(MOI.get(m, MOI.VariablePrimalStart(), x)), f_dest)
        if !isnan(val)
            is_tight = abs(val) < 1e-8
            has_start = true
        end
    end

    if method == ONE_ONE
        s1 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(pre_s1)
        s2 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.EqualTo(zero(T)))
        if pass_start && has_start
            MOI.set(m, MOI.VariablePrimalStart(), vb1, ifelse(is_tight, 1.0, 0.0))
            MOI.set(m, MOI.VariablePrimalStart(), vb2, ifelse(is_tight, 0.0, 1.0))
        end
    elseif method == ZERO_ZERO
        s1 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}(pre_s1)
        s2 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}(MOI.EqualTo(zero(T)))
        if pass_start && has_start
            MOI.set(m, MOI.VariablePrimalStart(), vb1, ifelse(is_tight, 0.0, 1.0))
            MOI.set(m, MOI.VariablePrimalStart(), vb2, ifelse(is_tight, 1.0, 0.0))
        end
    else
        s1 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(pre_s1)
        s2 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}(MOI.EqualTo(zero(T)))
        if pass_start && has_start
            MOI.set(m, MOI.VariablePrimalStart(), vb1, ifelse(is_tight, 1.0, 0.0))
        end
    end

    # do not MOIU.normalize_and_add_constraint because are vector functions
    c1 = MOI.add_constraint(m, to_vector_affine(f1), s1)
    c2 = MOI.add_constraint(m, to_vector_affine(f2), s2)

    if copy_names
        MOI.set(m, MOI.ConstraintName(), c1, "compl_ind1_($(nm))")
        MOI.set(m, MOI.ConstraintName(), c2, "compl_ind2_($(nm))")
    end
    return 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

function get_bounds(var, map, fallback_bound = Inf)
    if haskey(map, var)
        info = map[var]
        # TODO deal with precision and performance
        lower = ifelse(info.lower != -Inf, info.lower, -fallback_bound)
        upper = ifelse(info.upper != +Inf, info.upper, +fallback_bound)
        return IntervalArithmetic.interval(lower, upper)
    elseif 0.0 <= fallback_bound <= Inf
        return IntervalArithmetic.interval(-fallback_bound, fallback_bound)
    else
        error("variable $var has no finite bounds defined")
    end
end

function set_bound(inter::IntervalArithmetic.Interval, ::LT{T}) where T
    return inter.hi
end
function set_bound(inter::IntervalArithmetic.Interval, ::GT{T}) where T
    return inter.lo
end

function add_complement(mode::FortunyAmatMcCarlMode{T}, m, comp::Complement,
    idxmap_primal, idxmap_dual, copy_names::Bool, pass_start::Bool) where T

    f = comp.func_w_cte
    s = comp.set_w_zero
    v = comp.variable

    is_tight = false
    has_start = false

    if mode.with_slack
        slack, slack_in_set = MOI.add_constrained_variable(m, s)
    end
    f_dest = MOIU.map_indices.(Ref(idxmap_primal), f)

    f_bounds = MOIU.eval_variables(vi -> get_bounds(vi, mode.cache.map, mode.primal_big_M), f_dest)

    if pass_start
        val = MOIU.eval_variables(
            x-> nothing_to_nan(MOI.get(m, MOI.VariablePrimalStart(), x)), f_dest)
        if !isnan(val)
            is_tight = abs(val) < 1e-8
            has_start = true
        end
    end

    if mode.with_slack
        new_f = MOIU.operate(-, T, f_dest, MOI.SingleVariable(slack))
        equality = MOIU.normalize_and_add_constraint(m, new_f, MOI.EqualTo(zero(T)))
        if pass_start && has_start
            MOI.set(m, MOI.VariablePrimalStart(), slack, val)
        end
    end

    dual = idxmap_dual[v]
    v_bounds = get_bounds(dual, mode.cache.map, mode.dual_big_M)

    bin = MOI.add_variable(m)
    if pass_start && has_start && is_tight
        MOI.set(m, MOI.VariablePrimalStart(), bin, 1.0)
        MOI.set(m, MOI.VariablePrimalStart(), dual, 0.0)
    else
        MOI.set(m, MOI.VariablePrimalStart(), bin, 0.0)
    end

    s1 = flip_set(s)
    s2 = flip_set(s)

    Ms = set_bound(f_bounds, s1)
    Mv = set_bound(v_bounds, s2)

    if isnan(Ms) || abs(Ms) >= Inf || isnan(Mv) || abs(Mv) >= Inf
        error("It was not possible to automatically compute bounds"*
            " for a complementarity pair, please add the arguments"*
            " primal_big_M and dual_big_M to FortunyAmatMcCarlMode")
    end

    if mode.with_slack
        f1 = MOI.ScalarAffineFunction{T}(
            MOI.ScalarAffineTerm{T}.(
                [one(T), -Ms], [slack, bin]
            ),
            0.0
        )
    else
        push!(f_dest.terms, MOI.ScalarAffineTerm{T}(-Ms, bin))
        f1 = f_dest
    end

    f2 = MOI.ScalarAffineFunction{T}(
        MOI.ScalarAffineTerm{T}.(
            [one(T), Mv], [dual, bin]
        ),
        -Mv
    )

    c1 = MOIU.normalize_and_add_constraint(m, f1, s2)
    c2 = MOIU.normalize_and_add_constraint(m, f2, s2)
    c3 = MOI.add_constraint(m, MOI.SingleVariable(bin), MOI.ZeroOne())

    if copy_names
        nm = MOI.get(m, MOI.VariableName(), dual)
        if mode.with_slack
            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))")
        end
        MOI.set(m, MOI.VariableName(), bin, "bin_($(nm))")
        MOI.set(m, MOI.ConstraintName(), c1, "compl_fa_sl_($(nm))")
        MOI.set(m, MOI.ConstraintName(), c2, "compl_fa_dl_($(nm))")
        MOI.set(m, MOI.ConstraintName(), c2, "compl_fa_bn_($(nm))")
    end

    # if mode.with_slack
    #     return slack, slack_in_set, equality, c1
    # else
    # end
    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 handle_lower_objective_sense(lower::MOI.ModelLike)
    lower_objective_sense = MOI.get(lower, MOI.ObjectiveSense())
    if lower_objective_sense == MOI.FEASIBILITY_SENSE
        throw(ErrorException("Lower level models with objective_sense: " *
                            lower_objective_sense *
                            " are not supported."))
    end
    return
end


text:
# push!(LOAD_PATH,"../src/")
using Documenter, Collatz
# ~ Frob the operational context of jldoctest blocks with a using <this>
DocMeta.setdocmeta!(Collatz, :DocTestSetup, :(using Collatz); recursive=true)
@time makedocs(
    sitename="Collatz",
    modules=[Collatz],
    checkdocs = :exports,
    strict = true,
    format = Documenter.HTML(collapselevel = 1),
    pages = [
        "Introduction" => "index.md",
        "Functions" => "functions.md",
        "Extra" => [
            "Internal Modules" => "internal_modules.md",
            "Constants" => "constants.md",
            "Internal Functions" => "internal_functions.md"
        ]
    ],
)

ext: jl

max_stars_repo_name: Skenvy/Collatz

max_issues_repo_name: Skenvy/Collatz

max_forks_repo_name: Skenvy/Collatz


text:
using MINLPTests

solver = JuMP.optimizer_with_attributes(EAGO.Optimizer, "relative_tolerance" => 1E-9)

minlp_nlp_exclude = String[
    "001_010", # Unbounded box, check solution bad if not gradient-based....
    "002_010", # Unbounded box
    "004_010", # Unbounded box
    "004_011", # Unbounded box
    "005_010", # Unbounded box
    "006_010", # Unbounded box
    "007_010", # Unbounded box
    "008_010", # Unbounded box
    "008_011", # Unbounded box
    "005_011"  # \ operator not in JuMP
]
MINLPTests.test_nlp(solver, exclude = minlp_nlp_exclude,
                            objective_tol = 1E-3,
                            primal_tol = 1E-3,
                            dual_tol = NaN,
                            termination_target = MINLPTests.TERMINATION_TARGET_GLOBAL,
                            primal_target = MINLPTests.PRIMAL_TARGET_GLOBAL)

minlp_nlp_cvx_exclude = String[
    "001_010",
    "001_011", # convex quadratic objective... (linear unbounded...)
    "002_010",
    "002_011", # unbounded linear problem & convex quadratic objective
    "101_010", # convex quadratic nl constraints...
    "101_011",
    "101_012",
    "102_010",
    "102_011",
    "102_012",
    "102_013",
    "102_014",
    "103_010",
    "103_011",
    "103_012",
    "103_013",
    "103_014",
    "104_010",
    "105_010",
    "105_011",
    "105_012",
    "105_013",
    "106_010", # simple bounded domain
    "106_011", #
    "107_010",
    "107_011",
    "107_012",
    "108_010",
    "108_011",
    "108_012",
    "108_013",
    "109_010",
    "109_011",
    "109_012",
    "110_010",
    "110_011",
    "110_012",
    "201_010",
    "201_011",
    "202_010",
    "202_011",
    "202_012",
    "202_013",
    "202_014",
    "203_010",
    "204_010",
    "205_010",
    "206_010",
    "210_010",
    "210_011",
    "210_012",
    "501_010",
    "501_011"
]
MINLPTests.test_nlp_cvx(solver, exclude = minlp_nlp_cvx_exclude,
                                objective_tol = 1E-3,
                                primal_tol = 1E-3,
                                dual_tol = NaN,
                                termination_target = MINLPTests.TERMINATION_TARGET_GLOBAL,
                                primal_target = MINLPTests.PRIMAL_TARGET_GLOBAL)

minlp_nlp_mi_exclude = String[
    "001_010",  # no box constraints
    "002_010",

    #"003_010",   # TODO: Fix 003_010 - 003_016
    "003_011",  # FAIL
    "003_012",  # FAIL
    "003_013",  # FAIL
    "003_014",  # FAIL Never converges...
    "003_015",  #FAIL
    "003_016",

    "004_010",
    "004_011",
    "004_012",
    "005_010",
    "005_011",   # \ operator not in JuMP
    "006_010",
    "007_010",
    #"007_020"    # no way of specifying
]
MINLPTests.test_nlp_mi(solver, exclude = minlp_nlp_mi_exclude,
objective_tol = 1E-3, primal_tol = 1E-3,
termination_target = MINLPTests.TERMINATION_TARGET_GLOBAL,
primal_target = MINLPTests.PRIMAL_TARGET_GLOBAL)

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
using DataStructures

# Parse function and type documentation from julia source files.

# Metadata parsed from a declarations's preceeding comment.
type DeclarationComment
    decltype::String
    description::String
    args::Union(OrderedDict, Nothing)
    sections::Dict

    function DeclarationComment(decltype, description)
        new(decltype, description, nothing, Dict())
    end
end

# pattern to extract comments immediately preceeding declarations:
const decl_comment_pat =
    r"
    ((?:\h*\#[^\n]*\n)+)
    \h*(function|macro|immutable|abstract|type|const)\s+([A-Za-z_][A-Za-z0-9_?!]*)
    "xm

# pattern to extract comments immediately preceeding
# assignment style function declarations.
const compact_decl_comment_pat =
    r"
    ((?:\h*\#[^\n]*\n)+)              # comments
    \h*([A-Za-z_][A-Za-z0-9_?!]*)     # function name
    \h*((?:\{(?:[^{}]++|(?-1))*+\})?) # type parameters (balanced braces)
    \h*(\((?:[^()]++|(?-1))*+\))      # balanced parenthesis
    \h*=[^=]                          # assignment and not ==
    "xm

# pattern to strip leading whitespace and '#' characters from comments.
const comment_strip_pat = r"\h*\#+"

# Extract comments immediately preceeding declarations.
#
# Args:
#   input: Text which is expected to be julia source code.
#
# Returns:
#   A three-tuple of the form (decltype, name, comment), where
#   decltype is one of "function", "immutable", "type".
#
function extract_declaration_comments(input::String)
    mats = {}
    comment_strip(txt) = replace(txt, comment_strip_pat, "")

    for mat in eachmatch(decl_comment_pat, input)
        push!(mats, (mat.captures[2], mat.captures[3], comment_strip(mat.captures[1])))
    end

    for mat in eachmatch(compact_decl_comment_pat, input)
        push!(mats, ("function", mat.captures[2], comment_strip(mat.captures[1])))
    end
    mats
end

# Harvest declaration comments from a package
#
# Args:
#   package: Name of a currently installed package.
#
# Returns:
#   A dictionary mapping declared identifiers to DeclarationComment objects.
#
function harvest(package::String)
    srcdir = joinpath(Pkg.dir(package), "src")
    filenames = {}
    for filename in walkdir(srcdir)
        if match(r"\.jl$", filename) != nothing
            push!(filenames, filename)
        end
    end
    harvest(filenames)
end

# Extract and parse decleration-preceeding comments from a set of files.
#
# Args:
#   filename: Names of julia source files to harvest from.
#
# Returns:
#   A dictionary mapping declared identifiers to DeclarationComment objects.
#
function harvest(filenames::Vector)
    declarations = Dict()
    for filename in filenames
        for (decltype, name, comment) in extract_declaration_comments(readall(filename))
            if decltype == "macro"
                name = string("@", name)
            end
            declarations[name] = parse_comment(decltype, comment)
        end
    end
    declarations
end

# Return a substring of input that is of equal or greater indentation than the
# line starting at i.
function get_indented_block(input::String, i::Integer=1)
    output = IOBuffer()
    n = length(input)
    block_indent = 0
    while i < n && (input[i] == ' ' || input[i] == '\t')
        block_indent += 1
        write(output, input[i])
        i = nextind(input, i)
    end

    j = search(input, '\n', i)
    if j == 0
        j = n
    end

    write(output, input[i:j])
    i = nextind(input, j)

    while i < n
        j = i
        line_indent = 0
        while j < n && (input[j] == ' ' || input[j] == '\t')
            line_indent += 1
            j = nextind(input, j)
        end

        k = search(input, '\n', j)
        if k == 0
            k = n
        end

        if line_indent < block_indent && j < k
            break
        end

        write(output, input[i:k])
        i = nextind(input, k)
    end
    takebuf_string(output)
end

# TODO: We should not just take fixed section names. Anything should be allowed.
# pattern used to find fileds in a function declaration comment
const func_field_pat = r"^\h*(Args|Returns|Modifies|Throws)\s*:\h*\r?\n"im

# Parse a function declaration comment
#
# Args:
#   input: Comment preceeding a funciton declaration, stripped of the leading
#          '#' characters.
#
# Returns:
#   A DeclarationComment object.
#
function parse_comment(decltype, input::String)
    mat = match(func_field_pat, input)
    if mat == nothing
        return DeclarationComment(decltype, strip(input))
    end

    metadata = DeclarationComment(decltype, input[1:mat.offset-1])
    while mat != nothing
        typ = mat.captures[1]
        if typ == "Args"
            args = get_indented_block(input, mat.offset + length(mat.match))
            offset = mat.offset + length(args)
            metadata.args = parse_comment_args(args)
        else
            section = strip(get_indented_block(input, mat.offset + length(mat.match)))
            metadata.sections[typ] = section
            offset = mat.offset + length(section)
        end

        mat = match(func_field_pat, input, offset)
    end

    metadata
end

# pattern to match argument names/descriptions.
const arg_desc_pat = r"^(\h*)([\w_][\w\d_\!]*(?:\.\.\.)?)\h*:\h*(.*)\r?"m

# Parse the "Args" section of a function declaration comment.
#
# Args:
#   input: Everything under the "Args" block.
#
# Returns:
#   A dictionary mapping arugment names to their descriptions.
#
function parse_comment_args(input::String)
    args = OrderedDict()
    mat = match(arg_desc_pat, input)
    while mat != nothing
        indent = length(mat.captures[1])
        name = mat.captures[2]
        desc = IOBuffer()
        write(desc, mat.captures[3])

        offset = mat.offset + length(mat.match)
        if offset < length(input) && input[offset] == '\n'
            offset = nextind(input, offset)
            next_line_indent = 0
            i = offset
            while i < length(input) && (input[i] == ' ' || input[i] == '\t')
                next_line_indent += 1
                i = nextind(input, i)
            end
            if next_line_indent > indent
                desc_rest = get_indented_block(input, offset)
                write(desc, desc_rest)
                offset += length(desc_rest)
            end
        end

        args[name] = takebuf_string(desc)
        mat = match(arg_desc_pat, input, offset)
    end

    args
end

# Generate markdown for each declaration comment.
function generate_declaration_markdown(declarations::Dict)
    declaration_markdown = Dict()
    for (name, comment) in declarations
        out = IOBuffer()

        # TODO: in pandoc > 1.10 we can apply a class directly to a header
        println(out, "
")

        println(out, comment.description)

        if comment.args != nothing
            println(out, "\n#### Args")
            for (arg_name, arg_desc) in comment.args
                @printf(out, "  * `%s`: %s\n", arg_name, arg_desc)
            end
        end

        for (section_name, section_content) in comment.sections
            println(out, "\n#### ", section_name)
            println(out, section_content)
        end

        println(out, "
")

        declaration_markdown[name] = takebuf_string(out)
    end
    declaration_markdown
end

# Insert function declaration comments into a markdown document.
function expand_declaration_docs(input::String, declaration_markdown::Dict)
    Mustache.render(input, declaration_markdown)
end


text:
function test_softmax_loss_layer(backend::Backend, tensor_dim, use_weights::Bool, T, eps)
  println("-- Testing SoftmaxLossLayer on $(typeof(backend)){$T} $(use_weights ? "(with weights)" : "")...")

  if use_weights
    op_dim = tensor_dim-1
  else
    op_dim = max(abs(rand(Int)) % tensor_dim, 1)
  end
  dims = rand(6:11, tensor_dim)
  dims_label = copy(dims); dims_label[op_dim] = 1
  dims = tuple(dims...)
  dims_label = tuple(dims_label...)
  println("    > $dims (operate on dimension $op_dim)")

  input = rand(T, dims) + convert(T, 0.01)

  input_blob = make_blob(backend, input)
  diff_blob = make_blob(backend, T, size(input))

  if use_weights
    weights = rand(T, dims[1:end-1]) + convert(T, 0.1)
    weights = weights .* (dims[op_dim] ./ sum(weights,op_dim))
  else
    weights = []
  end

  label = abs.(rand(Int, dims_label)) .% dims[op_dim]
  label = convert(Array{T}, label)
  label_blob = make_blob(backend, label)

  inputs = Blob[input_blob, label_blob]

  loss_weight = abs(rand(T))

  layer = SoftmaxLossLayer(bottoms=[:pred, :labels], weights=weights, normalize=:local, dim=op_dim, weight=loss_weight)
  state = setup(backend, layer, inputs, Blob[diff_blob])

  println("    > Forward")
  forward(backend, state, inputs)

  new_shape = [size(weights)..., 1]
  rep_shape = ones(Int, tensor_dim); rep_shape[end] = dims[end]
  weights = repeat(reshape(weights, new_shape...), inner=rep_shape)

  expected_loss = 0
  expected_grad = zeros(T, size(input))

  dim_pre, dim_prob, dim_post = split_dims(input, op_dim)
  canonical_input = reshape(input, dim_pre, dim_prob, dim_post)
  canonical_grad = reshape(expected_grad, dim_pre, dim_prob, dim_post)
  if !isempty(weights)
    weights = reshape(weights, (dim_pre, dim_prob, dim_post))
  end
  label = reshape(label, dim_pre, 1, dim_post)
  for i = 1:dim_pre
    for j = 1:dim_post
      pred = exp.(canonical_input[i,:,j])
      pred /= sum(pred)
      if isempty(weights)
        canonical_grad[i,:,j] = pred
        canonical_grad[i,round(Int, label[i,1,j])+1,j] -= 1
        expected_loss += -log(pred[round(Int, label[i,1,j])+1])
      else
        y = round(Int, label[i,1,j])+1
        canonical_grad[i,:,j] = pred .* weights[i,y,j]
        canonical_grad[i,y,j] -= weights[i,y,j]
        expected_loss += -log(pred[y]) * weights[i,y,j]
      end
    end
  end
  scale = dims[op_dim] / prod(dims)
  expected_loss *= scale * loss_weight
  expected_grad *= scale * loss_weight
  expected_grad = reshape(expected_grad, size(input))

  @test -eps < state.loss - expected_loss < eps

  println("    > Backward")
  backward(backend, state, inputs, Blob[diff_blob])
  grad = zeros(T, size(input))
  copy!(grad, diff_blob)

  @test all(-eps .< grad - expected_grad .< eps)

  shutdown(backend, state)
end
function test_softmax_loss_layer(backend::Backend, use_weights::Bool, T, eps)
  for i in [2,4,5]
    test_softmax_loss_layer(backend, i, use_weights, T, eps)
  end
end
function test_softmax_loss_layer(backend::Backend, T, eps)
  test_softmax_loss_layer(backend, false, T, eps)
  test_softmax_loss_layer(backend, true, T, eps)
end
function test_softmax_loss_layer(backend::Backend)
  test_softmax_loss_layer(backend, Float64, 1e-5)
  test_softmax_loss_layer(backend, Float32, 1e-3)
end

if test_cpu
  test_softmax_loss_layer(backend_cpu)
end
if test_gpu
  test_softmax_loss_layer(backend_gpu)
end


text:
using Test, LinearAlgebra
using Yao
using QuAlgorithmZoo
using Yao.Blocks

witness_vec = zeros(1<<4)
witness_vec[[0, 3, 5, 10, 12, 15].+1] .= 1
witness_op = matrixgate(Diagonal(-witness_vec))

nbit = 4
depth = 4
cnots = [1=>2, 3=>4, 2=>3, 4=>1]
gen = random_diff_circuit(nbit, depth, cnots, mode=:Merged) |> autodiff(:QC)

dispatch!(gen, :random)
loss_func = () -> expect(witness_op, zero_state(nbit) |> gen)

function train(gen::AbstractBlock, g_learning_rate::Real, niter::Int)
    diff_blocks = collect(gen, AbstractDiff)
    for i in 1:niter
        ggrad = opdiff.(()->zero_state(nbit)|>gen, diff_blocks, Ref(witness_op))
        dispatch!(-, gen, ggrad.*g_learning_rate)
        if i%5==1 println("Step $i, loss = $(loss_func())") end
    end
end

train(gen, 0.1, 50)

using PyPlot
pygui(true)
pl = apply!(zero_state(nbit), gen)|>probs
bar(0:1<<nbit-1, pl)

gcf()
display(gcf())


text:
# auxiliary/axpby.jl
#
# Simple wrapper for the operation of computing α * x + β * y, and two
# special singleton types (immutables) to remove any overhead of multiplication
# by one or addition by zero.

immutable Zero <: Integer
end
immutable One <: Integer
end

const _zero = Zero()
const _one = One()

axpby(α::One,  x, β::One,  y) = x+y
axpby(α::Zero, x, β::One,  y) = y
axpby(α::One,  x, β::Zero, y) = x
axpby(α::Zero, x, β::Zero, y) = zero(y)

axpby(α::One,  x, β,       y) = x+β*y
axpby(α::Zero, x, β,       y) = β*y
axpby(α,       x, β::Zero, y) = α*x
axpby(α,       x, β::One,  y) = α*x+y

axpby(α,       x, β,       y) = α*x+β*y


text:
#=
day4:
- Julia version: 1.5.2
- Author: Paul.Mealor
- Date: 2021-12-04
=#

using AoC, Test

using Base.Iterators: flatten

function createGame(lines)
    numbers = split(lines[1],",") .|> n->parse(Int, n)
    lines = lines[2:end]
    cards = Array{Int}(undef, 5,5,0)
    for line in 1:6:length(lines)
        card = flatten(lines[line+1:line+5] .|> strip .|> (l -> split(l, r" +")))
        card = parse.(Int, card) |> nn -> reshape(nn, 5,5)
        cards = cat(cards,card; dims=3)
    end
    return (numbers, cards)
end

function calcResult(number, card, ticks)
    return number * sum(filter(i->!ticks[i], eachindex(ticks)) .|> i->card[i])
end

function isWinner(card, ticks)
    for j in 1:size(card, 2)
        if reduce(&, ticks[:,j])
            return true
        end
    end
    for j in 1:size(card, 1)
        if reduce(&, ticks[j,:])
            return true
        end
    end
    return false
end

function bingoGame(lines)
    (numbers,cards) = createGame(lines)

    ticks = zeros(Bool, size(cards))

    cardsInPlay = [1:size(cards,3)...]
    results = []

    for number in numbers
        for i in [cardsInPlay...]
            card=view(cards, :,:,i)
            cardTicks=view(ticks, :,:,i)
            for j in eachindex(card)
                if card[j] == number
                    cardTicks[j] = true
                end
            end

            if isWinner(card, cardTicks)
                push!(results, (number, card, cardTicks))
                setdiff!(cardsInPlay, [i])
            end
        end
    end
    return results
end

part1(lines) = bingoGame(lines)[1] |> result -> calcResult(result...)
part2(lines) = bingoGame(lines)[end] |> result -> calcResult(result...)

@test part1(exampleLines(4,1)) == 4512
@test part2(exampleLines(4,1)) == 1924

lines(4) |> ll -> @time part1(ll) |> show
lines(4) |> ll -> @time part2(ll) |> show


text:
module Integrate

using Geotherm.Geometry: Point2D

export runge_kutta

# This is a helper function and should not be exported.
function runge_kutta_iter(p::Point2D, f, h = 0.01)
    x, y = p.x, p.y
    k1 = h * f(x, y)
    k2 = h * f(x + h / 2, y + k1 / 2)
    k3 = h * f(x + h / 2, y + k2 / 2)
    k4 = h * f(x + h, y + k3)
    return Point2D(x + h, y + (k1 + 2 * k2 + 2 * k3 + k4) / 6)
end

"""
    runge_kutta(p0::Point2D, f, h=0.01, n=1000)

Do a Runge–Kutta fourth-order integration of function `f` with starting point `p0`, step `n`, and step-size `h`.
"""
function runge_kutta(p0::Point2D, f, h = 0.01, n = 1000)
    trace = Point2D[p0]
    for i in 1:(n-1)
        p_next = runge_kutta_iter(trace[i], f, h)
        push!(trace, p_next)
    end
    return trace
end

end


text:
using FFMPEG
using Test

text_execute(f) = try
    f()
    return true
catch e
    @warn "can't execute" exception=e
    return false
end

@testset "FFMPEG.jl" begin
    FFMPEG.versioninfo()

    # Test run and parse output
    out = FFMPEG.exe(`-version`, collect=true)
    @test occursin("ffmpeg version ",out[1])

    out = FFMPEG.exe(`-version`, command=FFMPEG.ffprobe, collect=true)
    @test occursin("ffprobe version ",out[1])

    # Test different invokation methods
    @test text_execute(() -> FFMPEG.exe("-version"))
    @test text_execute(() -> FFMPEG.exe(`-version`))
    @test text_execute(() -> FFMPEG.exe(`-version`, collect=true))
    @test text_execute(() -> FFMPEG.ffmpeg_exe(`-version`))
    @test text_execute(() -> FFMPEG.ffprobe_exe(`-version`))
    @test text_execute(() -> ffmpeg`-version`)
    @test text_execute(() -> ffprobe`-version`)
    @test text_execute(() -> @ffmpeg_env run(`$ffmpeg -version`))
end


text:
using LinearAlgebra: real
using PyPlot, SpecialFunctions, LinearAlgebra
using Peacock
polarisation = TM
eye(M) = Matrix{ComplexF64}(I, M, M)
diag_R(mat) = mat[diagind(mat)]

function mask(M)
    mat = zeros(M,M)
    for i = -(M-1):M-1
        mat = mat + diagm(i=>(-1)^(i)*ones(M-abs(i)))
    end
    return mat
end

r = 0.3
epf(x, y) = x^2 + y^2 < r^2 ? 9 : 1
muf(x,y) = 1
a1 = [1, 0]; a2 = [0, 1]
d1 = 0.001; d2 = 0.001
#geometry = Geometry(epf,muf,a1,a2,d1,d2)
#solver = Solver(geometry,11,11) #Solver里的Kx，Ky和下面的反过来了，相同的时候还是有问题，应该是a1，a2把方向给固定了

m = 5
M = (2*m+1)^2

kx = [i for i in -m:m]
kx = repeat(kx,1, 2*m+1)
ky = kx'
K = hcat(kx[:],ky[:])#为何没有*2pi
Kp = K
#K = hcat(diag_R(solver.Ky),diag_R(solver.Kx))/2/pi
#Kp = K

Kxp = [i for i in Kp[:,1]]
Kxp = repeat(Kxp,1,length(Kp[:,1]))
Kx = copy(Kxp)
Kxp = Kxp'

Kyp = [i for i in Kp[:,2]]
Kyp = repeat(Kyp,1,length(Kp[:,1]))
Ky = copy(Kyp)
Kyp = Kyp'

function geo(kx,ky)
    k=sqrt.(kx.^2 .+ ky.^2)
    eb = 1
    eab=8
    r=0.2
    ep=eab*r*besselj.(1,2*pi*k*r)./k.*cos.(pi*kx)
    ep[k.==0] .= (eb+eab*pi*r^2)
    return ep
end

#ep0 = solver.epc #不一样的原因是Peacock里Kx-Kxp对应的是qs，Ky-Kyp对应的是ps
#ep2 = geo(Ky-Kyp,Kx-Kxp)
epc = geo(Kx-Kxp,Ky-Kyp)
#epc = ep0.*mask(size(ep0)[1])
muc = 1

ii = 1
kym = (-1:0.01:1)*0.5
for omega = 0.4
kTD = []
kTDi = []
kTA = []
for ky in kym

    LHS11 = zeros(ComplexF64, M, M)
    LHS12 = eye(M)
    Kp = K + ones(M) * [0 ky]

    if polarisation == TE
        LHS21 = epc * (omega^2*muc*eye(M) .- inv(epc).*(Kp*transpose(Kp)))
        LHS22 = -epc*(inv(epc).*(K[:,1]*ones(1,M)+ones(M)*K[:,1]'))
        LHS = [LHS11 LHS12; LHS21 LHS22]
        RHS = epc
        label = L"H_z"

    elseif polarisation == TM
        LHS21 = (omega^2*epc - diagm(K[:,1].^2 + (ky .+ K[:,2]).^2))
        LHS22 = diagm(-2*K[:,1])
        LHS = [LHS11 LHS12; LHS21 LHS22]
        RHS = muc

        label = L"E_z"
    end
    kT, VT = eigen(LHS,sortby=x -> -abs(x)) #TODO注意这里的排序可能不同, 排序对结果不影响
    #kT, VT = eigen(LHS)
    VT = VT[1:M, :]
    #归一化，不用归一化也行
    VT = VT * diagm(vec(1 ./sqrt.(sum(abs.(VT).^2,dims=1))))
    maxVT,idVT=findmax(abs.(VT), dims=1)
    VT=VT*diagm(vec(maxVT./VT[idVT]))

    dT=1/omega*diag_R(VT'*inv(RHS)*(diagm(K[:,1])*VT+VT*diagm(kT)))

    i0= (abs.(real.(dT)) .> abs.(imag.(dT))).*((real.(dT)+imag.(dT)) .< 0) + (abs.(real.(dT)) .< abs.(imag.(dT))).*(imag.(kT).< 0)
    i0 = Bool.(i0)

    i3 = []
    for i4 = 0:1
        i0 = .!i0
        i1= 1:2*M
        i1 = i1[i0]
        i2 = sortperm(abs.(real.(kT[i0])))
        buf = sort(abs.(real.(kT[i0])))
        i1=i1[i2]
        i2 = (abs.([imag.(kT[i1]); Inf] - [Inf; imag.(kT[i1])] ) .> 1e-10) .| (abs.([imag.(kT[i1]); Inf]) .<= 1e-10) .| (abs.(abs.([real.(kT[i1]); Inf]) .- abs.([Inf; real.(kT[i1])])) .> 1e-10) .| (abs.(abs.([real.(kT[i1]); Inf]) .+ abs.([Inf; real.(kT[i1])]) .- 1) .>= 0.1)
        i1=i1[i2[1:(length(i2)-1)]]
        i3=[i3; i1[1:2m+1]]
    end
    kT=kT[i3]
    dT=dT[i3]
    VT=VT[:,i3]

    buffer=copy(kT)
    i1 = abs.(imag.(buffer)) .< 1e-6
    buffer1=real.(buffer)
    buffer1[.!i1] .= NaN #剔除所有虚部绝对值大于1e-5的值
    if length(buffer1) < length(kT)
        buffer1[length(buffer1)+1:lenm] .= NaN
    end
#    push!(kTD, buffer1')
	kTD = [kTD; buffer1[:]] #所有纯实数解
    buffer[i1] .= NaN
    if length(buffer) < length(kT)
        buffer[length(buffer)+1:lenm] .= NaN
    end
#    push!(kTDi, buffer')
	kTDi = [kTDi; buffer[:]] #所有复数解
    kTA = [kTA; kT[:]]
end
kTD = reshape(kTD,Int(length(kTD)/length(kym)),length(kym))
kTDi = reshape(kTDi,Int(length(kTDi)/length(kym)),length(kym))
kTA = reshape(kTA,Int(length(kTA)/length(kym)),length(kym))
kTD = transpose(kTD)
kTDi = transpose(kTDi)
kTA = transpose(kTA)
figure()
plot(real(kTA[:]),(kym*ones(1,2*(2m+1)))[:],".")
plot(imag(kTA[:]),(kym*ones(1,2*(2m+1)))[:],".")
xlabel(L"k_x(2\pi/a)")
ylabel(L"k_y(2\pi/a)")
title("ω=$omega(2πc/a)")
legend([L"Re(k_{x})",L"Im(k_x)"])
grid()
xlim([-1,1])
savefig("r_0p3_eps_9_omega_$omega.svg")
close()
end
#kTA
#plot(kTD,".")


text:
using Distributions
using Base.Test

const randi = Distributions.randi
n = 1_000_000

x = Int[randi(10) for i = 1:n]
@test min(x) == 1
@test max(x) == 10

x = Int[randi(3, 12) for i = 1:n]
@test min(x) == 3
@test max(x) == 12


text:
# The second type is ignored. For `T = VariableRef`, it should be `Float64` anyway
# since JuMP only supports `Float64`.
_promote_sum(T::Type, ::Type=Float64) = MA.promote_operation(+, T, T)
# `+` is not defined between `MOI.VariableIndex`.
_promote_sum(::Type{MOI.VariableIndex}, T::Type=Float64) = MOI.ScalarAffineFunction{T}

_promote_add_mul(T::Type) = MA.promote_operation(MA.add_mul, T, T, T)
_promote_add_mul(::Type{MOI.VariableIndex}) = MOI.ScalarQuadraticFunction{Float64}
function MP.polynomial(p::GramMatrix{MOI.VariableIndex, B, U}) where {B, U}
    Q = convert(Vector{U}, p.Q.Q)
    return MP.polynomial(GramMatrix(SymMatrix(Q, p.Q.n), p.basis))
end
#function MP.polynomial(p::GramMatrix{F}) where {F <: MOI.AbstractFunction}
#    MP.polynomial(p, MOIU.promote_operation(+, Float64, F, F))
#end

function primal_value(model, p::GramMatrix{MOI.VariableIndex})
    Q = MOI.get(model, MOI.VariablePrimal(), p.Q.Q)
    return GramMatrix(SymMatrix(Q, p.Q.n), p.basis)
end

_complex(::Type{T}, ::Type) where {T} = T
_complex(::Type{T}, ::Type{SumOfSquares.COI.HermitianPositiveSemidefiniteConeTriangle}) where {T} = Complex{T}

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
using Test
import LinearAlgebra.norm
@testset "Example: ManyModelsforWaterWaves" begin
    include("../examples/ManyModelsforWaterWaves.jl")

    Diff,=Integrate(1;μ=.1,ϵ=.1,p=2,N=2^8,L=10,T=1,dt=0.01,
            dealias=0,iterate=true,precond=true,method=1,maxiter=10,
            name=nothing)

    @test 1e-3< Diff[1] <2e-3   # SGN
    @test 1e-4< Diff[2] <4e-4   # WGN
    @test 1e-5< Diff[3] <4e-5   # IK2

    Diff,=Integrate(2;μ=.1,ϵ=.1,p=2,N=2^8,L=10,T=1,dt=0.01,
            dealias=0,iterate=true,precond=true,method=2,maxiter=10,
            name=nothing)

    @test 2e-3< Diff[1] <4e-3   # SGN
    @test 4e-4< Diff[2] <8e-4   # WGN
    @test 1e-4< Diff[3] <2e-4   # IK2

end

@testset "Example: Solitary Wave Whitham and KdV" begin
    include("../examples/SolitaryWaveWhitham.jl")

    a=PlotSolitaryWaveKdV(3)
    @test 4e-11 < a < 5e-11

    b=PlotSolitaryWaveWhitham(1.17)
    @test 0.0674 < b <  0.0675

end

@testset "Example: Solitary Wave Whitham-Boussinesq" begin
    include("../examples/SolitaryWaveWhithamBoussinesq.jl")

    b=IntegrateSolitaryWaveWhithamBoussinesq(N=2^6,dt=0.01/1.05)
    @test 9.9e-11 < b < 10e-11

    a=PlotSolitaryWaveWhithamBoussinesq(;c=1.1,α=1,L=20,N=2^9,μ=0.1,ϵ=0.1)
    @test 0.0316 < a < 0.0317

end

@testset "Example: Spectral Stability Green Naghdi" begin
    include("../examples/SpectralStabilityGreenNaghdi.jl")
    param = ell(0.3,0.005,0.75)
    σ,param,η,u,v,mesh=figspecCW(param;P=1,N=2^4)

    @test param == (h₀ = 0.3, h₁ = 0.3021875, h₂ = 0.305, a₀ = 0.3, a₁ = 0.0050000000000000044, H₀ = 0.3034497099913216, c = 0.166283361314354, λ = 5.189093852987355, m = 0.7499999999999991, κ = 0.3682704215595702)
    @test mesh == Mesh((N=2^4,L=param.λ))
    @test all((norm(σ),norm(η),norm(u),norm(v)) .≈ (92.75180153770384, 2.7862039840304336, 1.5268134085348848, 1.526814071924427))

    σ,η,u,v,mesh=figspecSW(1.2;L=10*π,N=2^6,M=10)
    @test mesh == Mesh((N=2^6,L=10*π))
    @test all((norm(σ),norm(η),norm(u),norm(v)) .≈ (95.64108213520367, 0.7411137440019041, 0.6632854772327135, 0.7064480663393101))

end

@testset "Study: Rectified WW2" begin
    include("../examples/StudyRectifiedWW2.jl")
    problem,blowup_time,blowup,error_energy=IntegrateWW2(init=1,μ=1,ϵ=0.1,L=20,N=2^6,T=1,dt = 0.01,dealias=0,δ=0,m=-1)
    @test (blowup_time,blowup) == (1.,false)
    @test 9e-7 < error_energy < 10e-7

    problem,blowup_time,blowup,error_energy=IntegrateWW2(init=1,μ=1,ϵ=10,L=20,N=2^14,T=0.1,dt = 0.01,dealias=0,δ=0,m=-1)
    @test blowup == true || blowup_time in (0.06,0.07)

end

@testset "Study: Whitham-Green-Naghdi" begin
    include("../examples/StudyWhithamGreenNaghdi.jl")

    norms = (13.564148081919193, 13.564148081919193, 12.01160049359638)

    η,u,v,mesh = PlotSolitaryWaveWGN1(c=2,N=2^9,L=7*π,verbose=false,name=nothing)
    @test all((norm(η),norm(η),norm(v)) .≈ norms )
    @test mesh == Mesh((N=2^9,L=7*π))

    η,u,v,mesh = PlotSolitaryWaveWGN2(c=2,N=2^9,L=7*π,verbose=false,name=nothing)
    @test all((norm(η),norm(η),norm(v)) .≈ norms )
    @test mesh == Mesh((N=2^9,L=7*π))

    η,u,v,mesh = PlotSolitaryWaveWGN3(c=2,N=2^9,L=7*π,verbose=false,name=nothing)
    @test all((norm(η),norm(η),norm(v)) .≈ norms )
    @test mesh == Mesh((N=2^9,L=7*π))

    Jac,Jacstar,FFT,IFFT = PlotJacobianWGN(;c=2,L=10*π,N=2^6,SGN=false,verbose=false,name=nothing)
    @test norm(Jac) .≈ 98.38287473770893 || norm(Jacstar) .≈ 73.20464101226274
    @test norm(FFT) == 2^6 || norm(IFFT) == 1

    pb = IntegrateSolitaryWaveWGN(;SGN=true,c=2,N=2^8,L=8*π,T=1,dt=1/100,name=nothing)
    η0,u0,v0=SolitaryWaveSerreGreenNaghdi((c=2,N=2^8,L=8*π,ϵ=1,μ=1); x₀ = 2)
    η,v=solution(pb)
    @test norm(η-η0)/norm(η) < 2e-7 || norm(v-v0)/norm(v) < 4e-8

    pb = StabilitySolitaryWaveWGN(;p=2,c=2,N=2^6,L=5*π,T=1,dt=10/10^3,SGN=true,precond=true,iterate=true,dealias=0,name=nothing)
    η,v=solution(pb)
    @test norm(η) ≈ 5.720544559359206 || norm(v) ≈ 5.080826080626654

    pb = IntegrateWGN(2;δ=0.1,N=2^6,L=3*π,x₀=-3,T= 1,dt = 10/10^3,SGN=false,dealias=0,iterate=true,precond=true,name=nothing)
    η,v=solution(pb)
    @test norm(η) ≈ 1.0808748139030515 || norm(v) ≈ 0.9929455087249658

end


text:
using BellScenario
using Test

@testset "BellScenario.jl" begin
    # Write your own tests here.
end


text:
################################################################################
#
#   SQUARE LATTICE
#
################################################################################

# Implementation
# - implementation 1
# - labels <: Any
function getUnitcellSquare(
            unitcell_type  :: Type{U},
            implementation :: Val{1}
        ) :: U where {LS,LB,S<:AbstractSite{LS,2},B<:AbstractBond{LB,2},U<:AbstractUnitcell{S,B}}

    # return a new Unitcell
    return newUnitcell(
        # Type of the unitcell
        U,
        # lattice vectors
        Vector{Float64}[
            Float64[1, 0],
            Float64[0, 1]
        ],
        # sites
        S[
            newSite(S, Float64[0,0], getDefaultLabelN(LS,1))
        ],
        # bonds
        B[
            newBond(B, 1,1, getDefaultLabel(LB), (+1,0)),
            newBond(B, 1,1, getDefaultLabel(LB), (-1,0)),
            newBond(B, 1,1, getDefaultLabel(LB), (0,+1)),
            newBond(B, 1,1, getDefaultLabel(LB), (0,-1))
        ]
    )
end

# Implementation
# - implementation 3
# - labels <: Any
function getUnitcellSquare(
            unitcell_type  :: Type{U},
            implementation :: Val{3}
        ) :: U where {LS,LB,S<:AbstractSite{LS,2},B<:AbstractBond{LB,2},U<:AbstractUnitcell{S,B}}

    # return a new Unitcell
    return newUnitcell(
        # Type of the unitcell
        U,
        # lattice vectors
        Vector{Float64}[
            Float64[1, 0],
            Float64[0, 1]
        ],
        # sites
        S[
            newSite(S, Float64[0.0,0.0], getDefaultLabelA(LS)),
            newSite(S, Float64[0.5,0.0], getDefaultLabelB(LS)),
            newSite(S, Float64[0.5,0.5], getDefaultLabelA(LS)),
            newSite(S, Float64[0.0,0.5], getDefaultLabelB(LS))
        ],
        # bonds
        B[
            newBond(B, 1,2, getDefaultLabelY(LB), (-1,0)),
            newBond(B, 1,4, getDefaultLabelY(LB), (0,-1)),
            newBond(B, 1,2, getDefaultLabelX(LB), (0,0)),
            newBond(B, 1,4, getDefaultLabelX(LB), (0,0)),

            newBond(B, 2,1, getDefaultLabelY(LB), (1,0)),
            newBond(B, 2,3, getDefaultLabelX(LB), (0,0)),
            newBond(B, 2,1, getDefaultLabelX(LB), (0,0)),
            newBond(B, 2,3, getDefaultLabelY(LB), (0,-1)),

            newBond(B, 3,4, getDefaultLabelX(LB), (0,0)),
            newBond(B, 3,2, getDefaultLabelX(LB), (0,0)),
            newBond(B, 3,4, getDefaultLabelY(LB), (1,0)),
            newBond(B, 3,2, getDefaultLabelY(LB), (0,1)),

            newBond(B, 4,3, getDefaultLabelX(LB), (0,0)),
            newBond(B, 4,1, getDefaultLabelY(LB), (0,1)),
            newBond(B, 4,3, getDefaultLabelY(LB), (-1,0)),
            newBond(B, 4,1, getDefaultLabelX(LB), (0,0)),
        ]
    )
end

# Implementation
# - implementation -3 (same as 3 but with pi flux)
# - labels <: Any
function getUnitcellSquare(
            unitcell_type  :: Type{U},
            implementation :: Val{-3}
        ) :: U where {LS,LB,S<:AbstractSite{LS,2},B<:AbstractBond{LB,2},U<:AbstractUnitcell{S,B}}

    # return a new Unitcell
    return newUnitcell(
        # Type of the unitcell
        U,
        # lattice vectors
        Vector{Float64}[
            Float64[1, 0],
            Float64[0, 1]
        ],
        # sites
        S[
            newSite(S, Float64[0.0,0.0], getDefaultLabelA(LS)),
            newSite(S, Float64[0.5,0.0], getDefaultLabelB(LS)),
            newSite(S, Float64[0.5,0.5], getDefaultLabelA(LS)),
            newSite(S, Float64[0.0,0.5], getDefaultLabelB(LS))
        ],
        # bonds
        B[
            newBond(B, 1,2, getDefaultLabelY(LB), (-1,0)),
            newBond(B, 1,4, getDefaultLabelY(LB), (0,-1)),
            newBond(B, 1,2, getDefaultLabelX(LB), (0,0)),
            newBond(B, 1,4, getDefaultLabelX(LB), (0,0)),

            newBond(B, 2,1, getDefaultLabelY(LB), (1,0)),
            newBond(B, 2,3, flipLabel(getDefaultLabelX(LB)), (0,0)),
            newBond(B, 2,1, getDefaultLabelX(LB), (0,0)),
            newBond(B, 2,3, flipLabel(getDefaultLabelY(LB)), (0,-1)),

            newBond(B, 3,4, getDefaultLabelX(LB), (0,0)),
            newBond(B, 3,2, flipLabel(getDefaultLabelX(LB)), (0,0)),
            newBond(B, 3,4, getDefaultLabelY(LB), (1,0)),
            newBond(B, 3,2, flipLabel(getDefaultLabelY(LB)), (0,1)),

            newBond(B, 4,3, getDefaultLabelX(LB), (0,0)),
            newBond(B, 4,1, getDefaultLabelY(LB), (0,1)),
            newBond(B, 4,3, getDefaultLabelY(LB), (-1,0)),
            newBond(B, 4,1, getDefaultLabelX(LB), (0,0)),
        ]
    )
end


text:
using BinDeps

@BinDeps.setup

if is_linux()
    hdf5 = library_dependency("libhdf5", aliases = ["libhdf5", "libhdf5_serial", "libhdf5_serial.so.10", "libhdf5_openmpi", "libhdf5_mpich"])
    provides(AptGet, "hdf5-tools", hdf5)
    provides(Pacman, "hdf5", hdf5)
    provides(Yum, "hdf5", hdf5)
end

if is_windows()
    using WinRPM
    hdf5 = library_dependency("libhdf5")
    provides(WinRPM.RPM, "hdf5", hdf5, os = :Windows )
end

if is_apple()
    using Homebrew
    hdf5 = library_dependency("libhdf5")
    provides(Homebrew.HB, "homebrew/science/hdf5", hdf5, os = :Darwin )
end

provides(Sources, URI("https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.10/hdf5-1.10.0-patch1/src/hdf5-1.10.0-patch1.tar.gz"), hdf5)
provides(BuildProcess, Autotools(libtarget = joinpath("src", "libhdf5.la")), hdf5)

@BinDeps.install Dict(:libhdf5 => :libhdf5)


text:
module VX
using CEnum

include("libopenvx_enums.jl")
include("libopenvx_constants.jl")
include("libopenvx_compat.jl")
include("libopenvx_types.jl")
include("libopenvx_api.jl")

end


text:
abstract type SVOptMethod end

struct mini <: SVOptMethod end

function (svhc::mini)(f, a, b; ϵ = 1e-5)
	ρ = (3-sqrt(5))/2
	w1 = a + (1.0 - ρ)*(b - a)
	w2 = a + ρ*(b - a)
	wL = min(w1, w2)
	wP = max(w1, w2)
	a1 = a
	b1 = b
	while abs(b1-a1) > ϵ
		if f(wL) < f(wP)
			b1 = wP
			wP = wL
			wL = b1 - ρ*(b1 - a1)
		else
			a1 = wL
			wL = wP
			wP = a1 + ρ*(b1 - a1)
		end
		wP1, wL1 = wP, wL
		wL = min(wP1, wL1)
		wP = max(wP1, wL1)
	end

return (a1+b1)/2.0

end

max_stars_repo_name: madzias/Julia

max_issues_repo_name: madzias/Julia

max_forks_repo_name: madzias/Julia


text:
using Test, SparseArrays, Missings

@testset "Missings" begin
    x = Missings.replace([1, 2, missing, 4], 3)
    @test eltype(x) === Int
    @test length(x) == 4
    @test size(x) == (4,)
    @test collect(x) == collect(1:4)
    @test collect(x) isa Vector{Int}
    x = Missings.replace([1, 2, missing, 4], 3.0)
    @test eltype(x) === Int
    @test length(x) == 4
    @test size(x) == (4,)
    @test collect(x) == collect(1:4)
    @test collect(x) isa Vector{Int}
    x = Missings.replace([1 2; missing 4], 3)
    @test eltype(x) === Int
    @test length(x) == 4
    @test size(x) == (2, 2)
    @test collect(x) == [1 2; 3 4]
    @test collect(x) isa Matrix{Int}
    x = Missings.replace((v for v in [missing, 1, missing, 2, 4]), 0)
    @test length(x) == 5
    @test size(x) == (5,)
    @test eltype(x) === Any
    @test collect(x) == [0, 1, 0, 2, 4]
    @test collect(x) isa Vector{Int}

    x = Missings.fail([1, 2, 3, 4])
    @test eltype(x) === Int
    @test length(x) == 4
    @test size(x) == (4,)
    @test collect(x) == [1, 2, 3, 4]
    @test collect(x) isa Vector{Int}
    x = Missings.fail([1 2; 3 4])
    @test eltype(x) === Int
    @test length(x) == 4
    @test size(x) == (2, 2)
    @test collect(x) == [1 2; 3 4]
    @test collect(x) isa Matrix{Int}
    @test_throws MissingException collect(Missings.fail([1, 2, missing, 4]))
    x = Missings.fail(v for v in [1, 2, 4])
    @test eltype(x) === Any
    @test length(x) == 3
    @test size(x) == (3,)
    @test collect(x) == [1, 2, 4]
    @test collect(x) isa Vector{Int}

    x = skipmissing([1, 2, missing, 4])
    @test eltype(x) === Int
    @test collect(x) == [1, 2, 4]
    @test collect(x) isa Vector{Int}
    x = skipmissing([1  2; missing 4])
    @test eltype(x) === Int
    @test collect(x) == [1, 2, 4]
    @test collect(x) isa Vector{Int}
    x = collect(skipmissing([missing]))
    @test eltype(x) === Union{}
    @test isempty(collect(x))
    @test collect(x) isa Vector{Union{}}
    x = collect(skipmissing(Union{Int, Missing}[]))
    @test eltype(x) === Int
    @test isempty(collect(x))
    @test collect(x) isa Vector{Int}
    x = skipmissing([missing, missing, 1, 2, missing, 4, missing, missing])
    @test eltype(x) === Int
    @test collect(x) == [1, 2, 4]
    @test collect(x) isa Vector{Int}
    x = skipmissing(v for v in [missing, 1, missing, 2, 4])
    @test eltype(x) === Any
    @test collect(x) == [1, 2, 4]
    @test collect(x) isa Vector{Int}

    @test levels(1:1) == levels([1]) == levels([1, missing]) == levels([missing, 1]) == [1]
    @test levels(2:-1:1) == levels([2, 1]) == levels([2, missing, 1]) == [1, 2]
    @test levels([missing, "a", "c", missing, "b"]) == ["a", "b", "c"]
    @test levels([Complex(0, 1), Complex(1, 0), missing]) == [Complex(0, 1), Complex(1, 0)]
    @test levels(sparse([0 3 2])) == [0, 2, 3]
    @test typeof(levels([1])) === typeof(levels([1, missing])) === Vector{Int}
    @test typeof(levels(["a"])) === typeof(levels(["a", missing])) === Vector{String}
    @test typeof(levels(sparse([1]))) === Vector{Int}
    @test isempty(levels([missing]))
    @test isempty(levels([]))

    @test Missings.T(Union{Int, Missing}) == Int
    @test Missings.T(Any) == Any
    @test Missings.T(Missing) == Union{}
    @test Missings.T(Union{Array{Int}, Missing}) == Array{Int}

    @test isequal(missings(1), [missing])
    @test isequal(missings(Int, 1), [missing])
    @test missings(Int, 1) isa Vector{Union{Int, Missing}}
    @test isequal(missings(Union{Int, Missing}, 1, 2), [missing missing])
    @test missings(Union{Int, Missing}, 1, 2) isa Matrix{Union{Int, Missing}}
    @test Union{Int, Missing}[1,2,3] == (Union{Int, Missing})[1,2,3]
    x = missings(Int, (1, 2))
    @test isa(x, Matrix{Union{Int, Missing}})
    @test isequal(x, [missing missing])
    x = missings((1, 2))
    @test isa(x, Matrix{Missing})
    @test isequal(x, [missing missing])

    @test allowmissing([1]) == [1]
    @test allowmissing([1]) isa AbstractVector{Union{Int, Missing}}
    @test allowmissing(Any[:a]) == [:a]
    @test allowmissing(Any[:a]) isa AbstractVector{Any}
    @test isequal(allowmissing([1, missing]), [1, missing])
    @test allowmissing([1, missing]) isa AbstractVector{Union{Int, Missing}}
    @test isequal(allowmissing([missing]), [missing])
    @test allowmissing([missing]) isa AbstractVector{Missing}

    @test allowmissing([1 1]) == [1 1]
    @test allowmissing([1 1]) isa AbstractArray{Union{Int, Missing}, 2}
    @test allowmissing([:a 1]) == [:a 1]
    @test allowmissing([:a 1]) isa AbstractArray{Any, 2}
    @test isequal(allowmissing([1 missing]), [1 missing])
    @test allowmissing([1 missing]) isa AbstractArray{Union{Int, Missing}, 2}
    @test isequal(allowmissing([missing missing]), [missing missing])
    @test allowmissing([missing missing]) isa AbstractArray{Missing, 2}

    @test disallowmissing(Union{Int, Missing}[1]) == [1]
    @test disallowmissing(Union{Int, Missing}[1]) isa AbstractVector{Int}
    @test disallowmissing([1]) == [1]
    @test disallowmissing([1]) isa AbstractVector{Int}
    @test disallowmissing(Any[:a]) == [:a]
    @test disallowmissing(Any[:a]) isa AbstractVector{Any}
    @test_throws MethodError disallowmissing([1, missing])
    @test_throws MethodError disallowmissing([missing])

    @test disallowmissing(Union{Int, Missing}[1 1]) == [1 1]
    @test disallowmissing(Union{Int, Missing}[1 1]) isa AbstractArray{Int, 2}
    @test disallowmissing([1 1]) == [1 1]
    @test disallowmissing([1 1]) isa AbstractArray{Int, 2}
    @test disallowmissing([:a 1]) == [:a 1]
    @test disallowmissing([:a 1]) isa AbstractArray{Any, 2}
    @test_throws MethodError disallowmissing([1 missing])
    @test_throws MethodError disallowmissing([missing missing])

    # Lifting
    @test passmissing(sqrt)(4) == 2.0
    @test isequal(passmissing(sqrt)(missing), missing)
    @test isequal(passmissing(sqrt).([missing, 4]), [missing, 2.0])
    @test passmissing((x,y)->"$x $y")(1, 2) == "1 2"
    @test isequal(passmissing((x,y)->"$x $y")(missing), missing)
    @test_throws ErrorException passmissing(string)(missing, base=2)
end


text:
using Thrift

function metadata(path)
    io = open(path)
    sz = filesize(io)
    seek(io, sz - SZ_PAR_MAGIC - SZ_FOOTER)

    # read footer size as little endian signed Int32
    meta_len = read(io, Int32)
    datasize = sz - meta_len - 2SZ_PAR_MAGIC - SZ_FOOTER
    seek(io, SZ_PAR_MAGIC + datasize)
    filemetadata = read_thrift(io, PAR2.FileMetaData)
end


text:
#########################
### Forest fire model ###
#########################
using Agents
using Random

mutable struct Tree <: AbstractAgent
  id::Int
  pos::Tuple{Int, Int}
  status::Int  # 1: green, 2: burning, 3: burned
end

# we can put the model initiation in a function
function model_initiation(;d, griddims, seed)
  Random.seed!(seed)

  space = GridSpace(griddims, moore = true)

  properties = Dict(:d => d)
  forest = ABM(Tree, space; properties=properties, scheduler=random_activation)

  # create and add trees to each node with probability d, which determines the density of the forest
  for node in 1:nv(forest)
    pp = rand()
    if pp <= forest.properties[:d]
      # Set all trees in the first column on fire.
      if model[node].pos[1] == 1
        tree = Tree(node, (1,1), 2)
      else
        tree = Tree(node, (1,1), 1)
      end
      add_agent!(tree, node, forest)
    end
  end
  return forest
end

function tree_step!(Tree, forest)
  if Tree.status == 2
    Tree.status = 3
    neighbor_cells = node_neighbors(tree, forest)
    for cell in neighbor_cells
      treeid = get_node_contents(cell, forest)
      if length(treeid) != 0 # the cell is not empty
        treen = forest.agents[treeid[1]]
        if treen.status == 1
          treen.status = 2
        end
      end
    end
  end
end


text:
struct StatelessSoftmax{T} <: AbstractStatelessPolicy where {T}
    θ::T

    function StatelessSoftmax(::Type{T}, num_actions::Int) where {T}
        θ = zeros(T, num_actions)
        return new{typeof(θ)}(θ)
    end

    function StatelessSoftmax(num_actions::Int) where {T}
        return StatelessSoftmax(Float64, num_actions)
    end
end

function params(π::StatelessSoftmax)
    return (π.θ,)
end

function paramsvec(π::StatelessSoftmax)
    return copy(vec(π.θ))
end

function paramsfromvec!(π::StatelessSoftmax, θ)
    @. π.θ = θ
    return π
end

function rrule(::typeof(paramsfromvec!), π::StatelessSoftmax, θ)
    @. π.θ = θ
    function paramsfromvec_statelesssoftmax!(ȳ)
        dθ = ȳ.θ
        return NoTangent(), ȳ, dθ
    end
    return π, paramsfromvec_statelesssoftmax!
end

struct LinearSoftmax{T} <: AbstractPolicy where {T}
    θ::T

    function LinearSoftmax(::Type{T}, num_features::Int, num_actions::Int) where {T}
        θ = zeros(T, (num_features, num_actions))
        return new{typeof(θ)}(θ)
    end

    function LinearSoftmax(num_features::Int, num_actions::Int)
        return LinearSoftmax(Float64, num_features, num_actions)
    end
end

function params(π::LinearSoftmax)
    return (π.θ,)
end

function paramsvec(π::LinearSoftmax)
    return copy(vec(π.θ))
end

function paramsfromvec!(π::LinearSoftmax, θ)
    π.θ .= reshape(θ, size(π.θ))
    return π
end

function rrule(::typeof(paramsfromvec!), π::LinearSoftmax, θ)
    π.θ .= reshape(θ, size(π.θ))
    function paramsfromvec_softmax!(ȳ)
        dθ = vec(ȳ.θ)
        return NoTangent(), ȳ, dθ
    end
    return π, paramsfromvec_softmax!
end

struct SoftmaxBuffer{TA, TP, TS} <: Any
    action::TA
    p::TP
    ψ::TS

    function SoftmaxBuffer(π::StatelessSoftmax)
        T = eltype(π.θ)
        n = length(π.θ)
        p = zeros(T, n)
        ψ = (zero(π.θ),)
        action = zeros(Int, 1)
        return new{typeof(action),typeof(p), typeof(ψ)}(action, p, ψ)
    end

    function SoftmaxBuffer(π::LinearSoftmax)
        T = eltype(π.θ)
        n = size(π.θ, 2)
        p = zeros(T, n)
        ψ = (zero(π.θ),)
        action = zeros(Int, 1)
        return new{typeof(action),typeof(p),typeof(ψ)}(action, p, ψ)
    end
end

function softmax!(out::AbstractVector, x::AbstractVector)
    @. out = x
    softmax!(out)
    return nothing
end

function softmax!(x::AbstractVector)
    maxx = maximum(x)
    x .-= maxx
    @. x = exp(x)
    x ./= sum(x)
    return nothing
end

function rrule(::typeof(softmax!), x::AbstractVector)
    softmax!(x)
    function softmax_pullback(ȳ)
        b = dot(ȳ, x)
        gx = @. (ȳ - b) * x
        return NoTangent(), gx
    end
    return x, softmax_pullback
end

function softmax(x::Vector{T}) where {T}
    maxx = maximum(x)
    y = x .- maxx
    @. y = exp(y)
    tot = sum(y)
    z = @. y / tot
    return z
end

function rrule(::typeof(softmax), x::Vector{T}) where {T}
    y = softmax(x)
    function softmax_pullback(ȳ)
        b = dot(ȳ, y)
        gx = @. (ȳ - b) * y
        return NoTangent(), gx
    end
    return y, softmax_pullback
end

function sample_discrete(p)
    n = length(p)
    i = 1
    c = p[1]
    u = rand()
    while c < u && i < n
        c += p[i += 1]
    end

    return i
end

function (π::StatelessSoftmax)()
    p = softmax(π.θ)
    return Categorical(p)
end

function (π::StatelessSoftmax)(buff::SoftmaxBuffer)
    softmax!(buff.p, π.θ)
    return Categorical(buff.p)
end

function (π::StatelessSoftmax)(s)
    return π()
end

function (π::StatelessSoftmax)(buff::SoftmaxBuffer, s)
    return π(buff)
end

function logpdf(π::StatelessSoftmax, a)
    p = similar(π.θ)
    softmax!(p, π.θ)
    return log(p[a])
end

function rrule(::typeof(logpdf), π::StatelessSoftmax, a)
    p = similar(π.θ)
    softmax!(p, π.θ)
    logp = log(p[a])
    T = typeof(logp)
    function logpdf_stateless_softmax_pullback(Ȳ)
        p *= -T(1.0)
        p[a] += T(1.0)
        p *= Ȳ
        dθ = Tangent{typeof(π)}(θ=p)
        return NoTangent(), dθ, NoTangent()
    end
    return logp, logpdf_stateless_softmax_pullback
end

function logpdf!(buff::SoftmaxBuffer, π::StatelessSoftmax, a)
    softmax!(buff.p, π.θ)
    return log(buff.p[a])
end

function logpdf!(buff::SoftmaxBuffer, π::StatelessSoftmax, s, a)
    return logpdf!(buff, π, a)
end

function grad_logpdf!(ψ, π::StatelessSoftmax, a)
    softmax!(ψ[1], π.θ)
    logp = log(ψ[1][a])
    T = eltype(π.θ)
    @. ψ[1] *= -T(1.0)
    ψ[1][a] += T(1.0)
    return logp, ψ
end

function grad_logpdf!(buff::SoftmaxBuffer, π::StatelessSoftmax, a)
    return grad_logpdf!(buff.ψ, π, a)
end

function grad_logpdf!(buff::SoftmaxBuffer, π::StatelessSoftmax, s, a)
    return grad_logpdf!(buff.ψ, π, a)
end

function sample_with_trace!(ψ, action, π::StatelessSoftmax)
    softmax!(ψ[1], π.θ)
    a = sample_discrete(ψ[1])
    action[1] = a

    logp = log(ψ[1][a])
    T = typeof(logp)
    @. ψ[1] *= -T(1.0)
    ψ[1][a] += T(1.0)
    return a, logp, ψ
end

function sample_with_trace!(buff::SoftmaxBuffer, π::StatelessSoftmax)
    return sample_with_trace!(buff.ψ, buff.action, π)
end

function sample_with_trace!(ψ, action, π::StatelessSoftmax, s)
    return sample_with_trace!(ψ, action, π)
end

function sample_with_trace!(buff::SoftmaxBuffer, π::StatelessSoftmax, s)
    return sample_with_trace!(buff.ψ, buff.action, π)
end

function (π::LinearSoftmax)(s)
    T = eltype(π.θ)
    p = zeros(T, size(π.θ,2))
    mul!(p, π.θ', s)
    softmax!(p)
    return Categorical(p)
end

function rrule(π::LinearSoftmax, s)
    T = eltype(π.θ)
    p = zeros(T, size(π.θ,2))
    mul!(p, π.θ', s)
    softmax!(p)
    d = Categorical(p)
    function linear_softmax_dist_pullback(ȳ)
        b = dot(ȳ.p, p)
        dp = @. (ȳ.p - b) * p
        ψ = zero(π.θ)
        G = ψ'
        for i in 1:length(dp)
            @. G[i, :] = s * dp[i]
        end
        ds = zero(s)
        mul!(ds, π.θ, dp)
        dθ = Tangent{typeof(π)}(θ=ψ)
        return dθ, ds
    end
    return d, linear_softmax_dist_pullback
end

function (π::LinearSoftmax)(buff::SoftmaxBuffer, s)
    mul!(buff.p, π.θ', s)
    softmax!(buff.p)
    return Categorical(buff.p)
end

function logpdf(π::LinearSoftmax, s, a)
    return logpdf_softmax(π.θ, s, a)
end

function logpdf_softmax(θ, s, a)
    T = eltype(θ)
    p = zeros(T, size(θ,2))
    mul!(p, θ', s)
    softmax!(p)
    return log(p[a])
end

function rrule(::typeof(logpdf_softmax), θ, s, a)
    T = eltype(θ)
    p = zeros(T, size(θ,2))
    mul!(p, θ', s)
    softmax!(p)
    logp = log(p[a])

    function logpdf_linear_softmax_pullback(ȳ)
        p *= -T(1.0)
        p[a] += T(1.0)
        p *= ȳ
        ψ = zero(θ)
        G = ψ'
        for i in 1:length(p)
            @. G[i, :] = s * p[i]
        end
        ds = zero(s)
        mul!(ds, θ, p)
        return NoTangent(), ψ, ds, NoTangent()
    end
    return logp, logpdf_linear_softmax_pullback
end

function logpdf_softmax(θ, s::AbstractMatrix, a)
    T = eltype(θ)
    num_a = size(θ,2)
    n = size(s,2)
    p = zeros(T, (num_a, n))
    mul!(p, θ', s)
    logps = zeros(T, (1,n))
    for i in 1:n
        p_i = @view p[:, i]
        softmax!(p_i)
        logps[1,i] = log(p_i[a[i]])
    end
    return logps
end

function rrule(::typeof(logpdf_softmax), θ, s::AbstractMatrix, a)
    T = eltype(θ)
    num_a = size(θ,2)
    n = size(s,2)
    p = zeros(T, (num_a, n))
    mul!(p, θ', s)
    logps = zeros(T, (1,n))
    for i in 1:n
        p_i = @view p[:, i]
        softmax!(p_i)
        logps[1,i] = log(p_i[a[1,i]])
    end

    function logpdf_linear_softmax_batch_pullback(ȳ)
        ψ = zero(θ)
        G = ψ'
        for i in 1:n
            p[:, i] *= -one(T)
            p[a[1, i], i] += one(T)
            p[:, i] *= ȳ[1, i]
            for j in 1:num_a
                x = @view s[:, i]
                @. G[j, :] += x * p[j,i]
            end
        end
        ds = zero(s)
        mul!(ds, θ, p)
        return NoTangent(), ψ, ds, ZeroTangent()
    end
    return logps, logpdf_linear_softmax_batch_pullback
end

function logpdf!(buff::SoftmaxBuffer, π::LinearSoftmax, s, a)
    mul!(buff.p, π.θ', s)
    softmax!(buff.p)
    return log(buff.p[a])
end

function grad_logpdf!(ψ, π::LinearSoftmax, s, a)
    T = eltype(π.θ)
    p = zeros(T, size(π.θ,2))
    mul!(p, π.θ', s)
    softmax!(p)
    logpa = log(p[a])
    T = typeof(logpa)
    @. p *= -T(1.0)
    p[a] += T(1.0)
    G = ψ[1]'
    for i in 1:length(p)
        @. G[i, :] = s * p[i]
    end
    return logpa, ψ
end

function grad_logpdf!(buff::SoftmaxBuffer, π::LinearSoftmax, s, a)
    mul!(buff.p, π.θ', s)
    softmax!(buff.p)
    logpa = log(buff.p[a])
    T = typeof(logpa)
    @. buff.p *= -T(1.0)
    buff.p[a] += T(1.0)
    G = buff.ψ[1]'
    for i in 1:length(buff.p)
        @. G[i, :] = s * buff.p[i]
    end
    return logpa, buff.ψ
end

function sample_with_trace!(ψ, action, π::LinearSoftmax, s)
    T = eltype(π.θ)
    p = zeros(T, size(π.θ,2))
    mul!(p, π.θ', s)
    softmax!(p)
    action[1] = sample_discrete(p)
    logp = log(p[action[1]])
    T = typeof(logp)
    @. p *= -T(1.0)

    p[action[1]] += T(1.0)
    G = ψ[1]'
    for i in 1:length(p)
        @. G[i, :] = s * p[i]
    end
    return action[1], logp, ψ
end

function sample_with_trace!(buff::SoftmaxBuffer, π::LinearSoftmax, s)
    mul!(buff.p, π.θ', s)
    softmax!(buff.p)
    a = sample_discrete(buff.p)
    buff.action[1] = a
    logp = log(buff.p[a])
    T = typeof(logp)
    @. buff.p *= -T(1.0)

    buff.p[a] += T(1.0)
    G = buff.ψ[1]'
    for i in 1:length(buff.p)
        @. G[i, :] = s * buff.p[i]
    end
    return a, logp, buff.ψ
end


text:
function Jdist!(J::Vector{Float64},x::Vector{Float64},C::Float64,R0::Float64,
    R1::Float64,Ws::Float64,W1::Float64,rho::Float64,beta::Float64,c0::Float64,
    a1::Float64,a2::Float64,h::Float64)
    J[1] = -a1/(4*a2*Ws^4)*(0.125*(W1-Ws*x[1])+c0) - (-3*W1^4/(8192*Ws^4) +
        W1^3*x[1]^2/(2048*Ws^3) - W1^3*c0/(128*Ws^4) + 3*W1^2*x[1]^2/(4096*Ws^2) -
        3*W1^2*c0^2/(64*Ws^4) - 3*W1*x[1]^3/(2048*Ws) + 3*W1*x[1]^2*c0/(128*Ws^2) -
        3*W1*x[1]*c0^2/(32*Ws^3) + 5*x[1]^4/8192 - x[1]^3*c0/(64*Ws) +
        9*x[1]^2*c0^2/(64*Ws^2) - x[1]*c0^3/(2*Ws^3) + c0^4/(2*Ws^4));
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
"""
    BoundaryFunction{B, X1, X2, F}

A wrapper for user-defined boundary condition functions.
"""
struct BoundaryFunction{B, X1, X2, F} <: Function
    func :: F

    """
        BoundaryFunction{B, X1, X2}(func)

    A wrapper for user-defined boundary condition functions on the
    boundary specified by symbol `B` and at location `(X1, X2)`.

    Example
    =======
    julia> using Oceananigans: BoundaryCondition, BoundaryFunction, Flux, Cell

    julia> top_tracer_flux = BoundaryFunction{:z, Cell, Cell}((x, y, t) -> cos(2π*x) * cos(t))
    (::BoundaryFunction{:z,Cell,Cell,getfield(Main, Symbol("##7#8"))}) (generic function with 1 method)

    julia> top_tracer_bc = BoundaryCondition(Flux, top_tracer_flux);
    """
    function BoundaryFunction{B, X1, X2}(func) where {B, X1, X2}
        B ∈ (:x, :y, :z) || throw(ArgumentError("The boundary B at which the BoundaryFunction is
                                                to be applied must be either :x, :y, or :z."))
        new{B, X1, X2, typeof(func)}(func)
    end
end


text:
function apply_boundary_conditions_on_z_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :infinite, :fixed}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[:, :, 1, rbi] .= grid_boundary_factors[1] .* view(rbpot, :, :, 2, rbi) # infinity boundaries in r
    nothing
end
function apply_boundary_conditions_on_z_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :infinite, :infinite}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[:, :, 1, rbi] .= grid_boundary_factors[1] .* view(rbpot, :, :, 2, rbi)
    rbpot[:, :, end, rbi] .= grid_boundary_factors[2] .* view(rbpot, :, :, size(rbpot, 3) - 1, rbi)
    nothing
end
function apply_boundary_conditions_on_z_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :reflecting, :fixed}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[:, :, 1, rbi] .= view(rbpot, :, :, 3, rbi) # infinity boundaries in r
    nothing
end
function apply_boundary_conditions_on_z_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :fixed, :infinite}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[:, :, end, rbi] .= grid_boundary_factors[2] .* view(rbpot, :, :, size(rbpot, 3) - 1, rbi) # infinity boundaries in r
    nothing
end
function apply_boundary_conditions_on_z_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :fixed, :reflecting}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[:, :, end, rbi] .= view(rbpot, :, :, size(rbpot, 3) - 2, rbi) # infinity boundaries in r
    nothing
end
function apply_boundary_conditions_on_z_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :reflecting, :reflecting}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[:, :, 1, rbi] .= view(rbpot, :, :, 3, rbi)
    rbpot[:, :, end, rbi] .= view(rbpot, :, :, size(rbpot, 3) - 2, rbi)
    nothing
end
function apply_boundary_conditions_on_z_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :fixed, :fixed}, int::Interval{:closed, :closed, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    nothing
end
function apply_boundary_conditions_on_z_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :periodic, :periodic}, int::Interval{:closed, :open, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[:, :,   1, rbi] .= view(rbpot,  :, :, size(rbpot, 3) - 1, rbi)
    rbpot[:, :, end, rbi] .= view(rbpot,  :, :,       2, rbi)
    nothing
end

function apply_boundary_conditions_on_x_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :fixed, :fixed}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :infinite, :infinite}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[  1, :, :, rbi] .= grid_boundary_factors[1] .* view(rbpot, 2, :, :, rbi)  # hmm, this is probably not fully correct since this is the red-black dimension
    rbpot[end, :, :, rbi] .= grid_boundary_factors[2] .* view(rbpot, size(rbpot, 1) - 1, :, :, rbi)  # anyhow its just an approximation
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :fixed, :infinite}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[end, :, :, rbi] .= grid_boundary_factors[2] .* view(rbpot, size(rbpot, 1) - 1, :, :, rbi)    # infinity boundaries in x (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :infinite, :fixed}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[  1, :, :, rbi] .= grid_boundary_factors[1] .* view(rbpot, 2, :, :, rbi)  # infinity boundaries in x (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :infinite, :reflecting}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[  1, :, :, rbi] .= grid_boundary_factors[1] .* view(rbpot, 2, :, :, rbi)  # infinity boundaries in x (only ±1 because this is the compressed (redblack) dimension)
    rbpot[end, :, :, rbi] .= view(rbpot, size(rbpot, 1) - 1, :, :, rbi) # (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :reflecting, :infinite}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[1, :, :, rbi]   .= view(rbpot, 2, :, :, rbi) #  (only ±1 because this is the compressed (redblack) dimension)
    rbpot[end, :, :, rbi] .= grid_boundary_factors[2] .* view(rbpot, size(rbpot, 1) - 1, :, :, rbi)    # infinity boundaries in x (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :reflecting, :reflecting}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[1, :, :, rbi]   .= view(rbpot, 2, :, :, rbi) #  (only ±1 because this is the compressed (redblack) dimension)
    rbpot[end, :, :, rbi] .= view(rbpot, size(rbpot, 1) - 1, :, :, rbi) # (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :fixed, :reflecting}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[end, :, :, rbi] .= view(rbpot, size(rbpot, 1) - 1, :, :, rbi) # (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :reflecting, :fixed}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[1, :, :, rbi]   .= view(rbpot, 2, :, :, rbi) #  (only ±1 because this is the compressed (redblack) dimension)
    nothing
end
function apply_boundary_conditions_on_x_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :periodic, :periodic}, int::Interval{:closed, :open, T},
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[  1, :, :, rbi] .= view(rbpot, size(rbpot, 1), :, :, rbi)
    rbpot[end, :, :, rbi] .= view(rbpot, 1, :, :, rbi)
    nothing
end

function apply_boundary_conditions_on_y_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :fixed, :fixed}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :infinite, :infinite}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[ :, 1, :, rbi]   .= grid_boundary_factors[1] .* view(rbpot, :, 2, :, rbi)
    rbpot[ :, end, :, rbi] .= grid_boundary_factors[2] .* view(rbpot, :, size(rbpot, 2) - 1, :, rbi)
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :fixed, :infinite}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[ :, end, :, rbi] .= grid_boundary_factors[2] .* view(rbpot, :, size(rbpot, 2) - 1, :, rbi)
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :infinite, :fixed}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[ :, 1, :, rbi]   .= grid_boundary_factors[1] .* view(rbpot, :, 2, :, rbi)
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :reflecting, :reflecting}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[ :, 1, :, rbi]   .= view(rbpot, :, 3, :, rbi)
    rbpot[ :, end, :, rbi] .= view(rbpot, :, size(rbpot, 2) - 2, :, rbi)
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :fixed, :reflecting}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[ :, end, :, rbi] .= view(rbpot, :, size(rbpot, 2) - 2, :, rbi)
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :reflecting, :fixed}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[ :, 1, :, rbi]   .= view(rbpot, :, 3, :, rbi)
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :infinite, :reflecting}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[ :, 1, :, rbi]   .= grid_boundary_factors[1] .* view(rbpot, :, 2, :, rbi)
    rbpot[ :, end, :, rbi] .= view(rbpot, :, size(rbpot, 2) - 2, :, rbi)
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :reflecting, :infinite}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[ :, 1, :, rbi]   .= view(rbpot, :, 3, :, rbi)
    rbpot[ :, end, :, rbi] .= grid_boundary_factors[2] .* view(rbpot, :, size(rbpot, 2) - 1, :, rbi)
    nothing
end
function apply_boundary_conditions_on_y_axis!(  rbpot::AbstractArray{T, 4}, rbi::Int, ax::DiscreteAxis{T, :periodic, :periodic}, int::Interval,
                                                grid_boundary_factors::NTuple{2, T})::Nothing where {T}
    rbpot[ :, 1, :, rbi]   .= view(rbpot, :, size(rbpot, 2) - 1, :, rbi)
    rbpot[ :, end, :, rbi] .= view(rbpot, :, 2, :, rbi)
    nothing
end

function apply_boundary_conditions!(pcs::PotentialCalculationSetup{T, Cartesian}, update_even_points::Val{even_points}, only2d::Val{only_2d}) where {T, even_points, only_2d}
    rbi::Int = even_points ? rb_even::Int : rb_odd::Int
    apply_boundary_conditions_on_x_axis!( pcs.potential, rbi, pcs.grid.axes[1], pcs.grid.axes[1].interval, pcs.grid_boundary_factors[1])
    apply_boundary_conditions_on_y_axis!( pcs.potential, rbi, pcs.grid.axes[2], pcs.grid.axes[2].interval, pcs.grid_boundary_factors[2])
    apply_boundary_conditions_on_z_axis!( pcs.potential, rbi, pcs.grid.axes[3], pcs.grid.axes[3].interval, pcs.grid_boundary_factors[3])
    nothing
end


text:
using FlicsFlow
using Test

@testset "FlicsFlow.jl" begin
    # Write your tests here.
	casedir = "../examples/elbow"
	@test readopenfoammesh(casedir).nodes[1].index == 1
end


text:
export Cheb_dir_AD

function Cheb_dir_AD(nlp            :: AbstractNLPModel,
                     nlp_at_x       :: NLPAtX,
                     x              :: AbstractVector,
                     Nwtdirection   :: Function)

    dₙ = Nwtdirection(nlp_at_x.Hx, nlp_at_x.gx, verbose = false)

    ∇³fdₙdₙ = ∇f³xuv(nlp, x, dₙ, dₙ)

    dCp = Nwtdirection(nlp_at_x.Hx, ∇³fdₙdₙ, verbose = false)

    dC = dₙ + 0.5 * dCp

    return dC
end