text:
################################################################################
#                   CORE DISPATCHVARIABLEREF METHOD EXTENSIONS
################################################################################
# Extend dispatch_variable_ref
function dispatch_variable_ref(
    model::InfiniteModel,
    index::SemiInfiniteVariableIndex
    )::SemiInfiniteVariableRef
    return SemiInfiniteVariableRef(model, index)
end

# Extend _add_data_object
function _add_data_object(
    model::InfiniteModel,
    object::VariableData{<:SemiInfiniteVariable}
    )::SemiInfiniteVariableIndex
    return MOIUC.add_item(model.semi_infinite_vars, object)
end

# Extend _data_dictionary (type based)
function _data_dictionary(
    model::InfiniteModel,
    ::Type{SemiInfiniteVariable}
    )::MOIUC.CleverDict{SemiInfiniteVariableIndex, VariableData{SemiInfiniteVariable{GeneralVariableRef}}}
    return model.semi_infinite_vars
end

# Extend _data_dictionary (reference based)
function _data_dictionary(
    vref::SemiInfiniteVariableRef
    )::MOIUC.CleverDict{SemiInfiniteVariableIndex, VariableData{SemiInfiniteVariable{GeneralVariableRef}}}
    return JuMP.owner_model(vref).semi_infinite_vars
end

# Extend _data_object
function _data_object(
    vref::SemiInfiniteVariableRef
    )::VariableData{SemiInfiniteVariable{GeneralVariableRef}}
    object = get(_data_dictionary(vref), JuMP.index(vref), nothing)
    if isnothing(object)
        error("Invalid point variable reference, cannot find ",
              "corresponding variable in the model. This is likely ",
              "caused by using the reference of a deleted variable.")
    end
    return object
end

"""
    internal_semi_infinite_variable(vref::SemiInfiniteVariableRef,
                                    key::Val{:my_ext_key})::SemiInfiniteVariable

Return the semi-infinite variable object of `vref` assuming it is an internal variable
made during measure expansion within an optimizer model. This will apply to
optimizer model extensions that utilize `add_measure_variable` in combination
with `expand_measure`.
"""
function internal_semi_infinite_variable end

# Extend _core_variable_object
function _core_variable_object(
    vref::SemiInfiniteVariableRef
    )::SemiInfiniteVariable{GeneralVariableRef}
    if !haskey(_data_dictionary(vref), JuMP.index(vref))
        model = JuMP.owner_model(vref)
        return internal_semi_infinite_variable(vref, Val(optimizer_model_key(model)))
    else
        return _data_object(vref).variable
    end
end

# Extend _object_numbers
function _object_numbers(vref::SemiInfiniteVariableRef)::Vector{Int}
    return _core_variable_object(vref).object_nums
end

# Extend _parameter_numbers
function _parameter_numbers(vref::SemiInfiniteVariableRef)::Vector{Int}
    return _core_variable_object(vref).parameter_nums
end

################################################################################
#                             DEFINITION METHODS
################################################################################
## Process the value arguments to be formatted correctly
## (avoid AffExpr array promotions)
# Not an AffExpr
function _process_value(v)
    return v
end

# An AffExpr
function _process_value(v::JuMP.GenericAffExpr)
    if isempty(v.terms)
        return v.constant
    elseif length(v.terms) == 1 && isone(collect(values(v.terms))[1]) &&
           iszero(v.constant)
        return first(keys(v.terms))
    else
        return v
    end
end

"""
    SemiInfinite{V, VT <: VectorTuple} <: InfOptVariableType

A `DataType` to assist in making semi-infinite variables. This can be passed as an
extra argument to `@variable` to make such a variable:
```julia
@variable(model, var_expr, SemiInfinite(inf_var, parameter_values...), kwargs...)
```
Here `parameter_values` must match the format of the infinite parameter
references associated with the infinite variable `inf_var` and can be comprised
of both real valued supports and/or infinite parameters.

**Fields**
- `infinite_variable_ref::V`
- `parameter_values::VT`: The infinite parameters and/or infinite
   parameter support values the variable will depend on.
"""
struct SemiInfinite{V, VT <: Collections.VectorTuple} <: InfOptVariableType
    infinite_variable_ref::V
    parameter_values::VT
    function SemiInfinite(ivref::V, vt::VT) where {V, VT <: Collections.VectorTuple}
        processed_vals = _process_value.(vt.values)
        if !isequal(processed_vals, vt.values)
            vt2 = Collections.VectorTuple(processed_vals, vt.ranges, vt.indices)
            return new{V, typeof(vt2)}(ivref, vt2)
        end
        return new{V, VT}(ivref, vt)
    end
end
function SemiInfinite(ivref, vals...)
    vt = Collections.VectorTuple(vals)
    return SemiInfinite(ivref, vt)
end

"""
    JuMP.build_variable(_error::Function, info::JuMP.VariableInfo,
                        var_type::SemiInfinite)::SemiInfiniteVariable{GeneralVariableRef}

Build and return a semi-infinite variable based on `info` and `var_type`. Errors
if the information stored in `var_type` is invalid. See [`SemiInfinite`](@ref)
for more information.

**Example**
```julia-repl
julia> y
y(t, x)

julia> info = VariableInfo(false, 0, false, 0, false, 0, true, 0, false, false);

julia> semi_inf_var = build_variable(error, info, SemiInfinite(y, t, 0));
```
"""
function JuMP.build_variable(
    _error::Function,
    info::JuMP.VariableInfo,
    var_type::SemiInfinite;
    extra_kwargs...
    )::SemiInfiniteVariable{GeneralVariableRef}
    # check for unneeded keywords
    for (kwarg, _) in extra_kwargs
        _error("Keyword argument $kwarg is not for use with semi-infinite ",
               "variables.")
    end
    # ensure that all the info was left alone (TODO relax this)
    if info.has_lb || info.has_ub || info.has_fix || info.has_start ||
       info.integer || info.binary
        _error("Semi-infinite variables do not currently support the ",
               "specification of variable information fields (e.g., lower ",
               "bounds, upper bounds, start values, etc). These are simply ",
               "inherited from the corresponding infinite variable.")
    end
    # check the infinite variable reference
    ivref = var_type.infinite_variable_ref
    if !(ivref isa GeneralVariableRef)
        _error("Expected an infinite variable/derivative reference dependency ",
               "of type `GeneralVariableRef`, but got an argument of type ",
               "$(typeof(ivref)).")
    end
    # check that the values are the same format as the infinite variable
    raw_params = var_type.parameter_values
    if !all(v isa Union{Real, GeneralVariableRef} for v in raw_params)
        _error("Unexpected inputs given, expected them to be parameter references ",
               "and real numbers.")
    elseif !Collections.same_structure(raw_params, raw_parameter_refs(ivref))
        _error("The parameter reference/value input format $(raw_params) does ",
               "not match that of the infinite variable $(ivref).")
    end
    # make the evaluation support dictionary
    inds = findall(p -> p isa Real, raw_params)
    eval_supports = Dict{Int, Float64}(i => raw_params[i] for i in inds)
    # dispatch to the main builder
    return JuMP.build_variable(_error, ivref, eval_supports)
end

"""
    JuMP.build_variable(_error::Function, ivref::GeneralVariableRef,
                        eval_supports::Dict{Int, Float64}; [check::Bool = true]
                        )::SemiInfiniteVariable{GeneralVariableRef}

Extend the `JuMP.build_variable` function to build a semi-infinite variable
based on the infinite variable/derivative/parameter function `ivref` with
reduction support `eval_supports`. Will check that input is appropriate if
`check = true`. Errors if `ivref` is not an infinite variable, `eval_supports`
violate infinite parameter domains, or if the support dimensions don't match the
infinite parameter dimensions of `ivref`. This is intended an internal method for
use in evaluating measures.
"""
function JuMP.build_variable(
    _error::Function,
    ivref::GeneralVariableRef,
    eval_supports::Dict{Int, Float64};
    check::Bool = true
    )::SemiInfiniteVariable{GeneralVariableRef}
    # check the inputs
    dvref = dispatch_variable_ref(ivref)
    if check && !(dvref isa Union{InfiniteVariableRef, DerivativeRef, ParameterFunctionRef})
         _error("Must specify an infinite variable/derivative/parameter function dependency.")
    elseif check && maximum(keys(eval_supports)) > length(parameter_list(dvref))
        _error("Support evaluation dictionary indices do not match the infinite " *
               "parameter dependencies of $(ivref).")
    end
    prefs = parameter_list(dvref)
    if check
        for (index, value) in eval_supports
            pref = prefs[index]
            if JuMP.has_lower_bound(pref) && !supports_in_domain(value, infinite_domain(pref))
                _error("Evaluation support violates infinite parameter domain(s).")
            end
        end
    end
    # get the parameter object numbers of the dependencies
    object_nums = Int[]
    for i in eachindex(prefs)
        if !haskey(eval_supports, i)
            union!(object_nums, _object_number(prefs[i]))
        end
    end
    # get the parameter numbers
    orig_nums = _parameter_numbers(ivref)
    param_nums = [orig_nums[i] for i in eachindex(orig_nums)
                  if !haskey(eval_supports, i)]
    # round the support values in accordance with the significant digits
    for (k, v) in eval_supports
        eval_supports[k] = round(v, sigdigits = significant_digits(prefs[k]))
    end
    # build the variable
    return SemiInfiniteVariable(ivref, eval_supports, param_nums, object_nums)
end

## Make dispatch functions to add semi-infinite variable based supports
# Independent parameters
function _add_semi_infinite_support(
    prefs::Vector{IndependentParameterRef},
    eval_supps::Dict{Int, Float64},
    counter::Int
    )::Int
    for pref in prefs
        if haskey(eval_supps, counter)
            add_supports(pref, eval_supps[counter], check = false,
                        label = UserDefined)
        end
        counter += 1
    end
    return counter
end

# Dependent parameters
function _add_semi_infinite_support(
    prefs::Vector{DependentParameterRef},
    eval_supps::Dict{Int, Float64},
    counter::Int
    )::Int
    len = length(prefs)
    if all(haskey(eval_supps, i) for i in counter:counter+len-1)
        supp = Matrix{Float64}(undef, len, 1)
        for i in eachindex(supp)
            supp[i] = eval_supps[counter + i - 1]
        end
        add_supports(prefs, supp, check = false, label = UserDefined)
    end
    return counter + len
end

# Update the parameter supports as needed
function _update_param_supports(
    ivref::Union{InfiniteVariableRef, DerivativeRef},
    eval_supps::Dict{Int, Float64}
    )::Nothing
    raw_prefs = raw_parameter_refs(ivref)
    counter = 1
    for i in 1:size(raw_prefs, 1)
        prefs = dispatch_variable_ref.(raw_prefs[i, :])
        counter = _add_semi_infinite_support(prefs, eval_supps, counter)
    end
end

# Update the parameter supports as needed
function _update_param_supports(
    ivref::ParameterFunctionRef,
    eval_supps::Dict{Int, Float64}
    )::Nothing
    return
end

"""
    JuMP.add_variable(model::InfiniteModel, var::SemiInfiniteVariable,
                      [name::String = ""])::GeneralVariableRef

Extend the `JuMP.add_variable` function to accomodate `InfiniteOpt`
semi-infinite variable types. Adds `var` to the infinite model `model` and
returns a [`GeneralVariableRef`](@ref). Primarily intended to be an internal
function used in evaluating measures.
"""
function JuMP.add_variable(
    model::InfiniteModel,
    var::SemiInfiniteVariable,
    name::String = "";
    add_support = true
    )::GeneralVariableRef
    ivref = var.infinite_variable_ref
    divref = dispatch_variable_ref(ivref)
    eval_supps = var.eval_supports
    JuMP.check_belongs_to_model(divref, model)
    existing_index = get(model.semi_lookup, (ivref, eval_supps), nothing)
    if isnothing(existing_index)
        data_object = VariableData(var, name)
        vindex = _add_data_object(model, data_object)
        push!(_semi_infinite_variable_dependencies(divref), vindex)
        if add_support
            _update_param_supports(divref, eval_supps)
        end
        model.semi_lookup[(ivref, eval_supps)] = vindex
    else
        vindex = existing_index
        if !isempty(name)
            model.semi_infinite_vars[vindex].name = name
        end
    end
    return _make_variable_ref(model, vindex)
end

################################################################################
#                          PARAMETER REFERENCE METHODS
################################################################################
"""
    infinite_variable_ref(vref::SemiInfiniteVariableRef)::GeneralVariableRef

Return the infinite variable/derivative/parameter function reference associated
with the semi-infinite variable `vref`.

**Example**
```julia-repl
julia> infinite_variable_ref(vref)
g(t, x)
```
"""
function infinite_variable_ref(vref::SemiInfiniteVariableRef)::GeneralVariableRef
    return _core_variable_object(vref).infinite_variable_ref
end

"""
    eval_supports(vref::SemiInfiniteVariableRef)::Dict{Int, Float64}

Return the evaluation supports associated with the semi-infinite variable
`vref`.

**Example**
```julia-repl
julia> eval_supports(vref)
Dict{Int64,Float64} with 1 entry:
  1 => 0.5
```
"""
function eval_supports(vref::SemiInfiniteVariableRef)::Dict{Int, Float64}
    return _core_variable_object(vref).eval_supports
end

"""
    raw_parameter_refs(vref::SemiInfiniteVariableRef)::VectorTuple

Return the raw [`VectorTuple`](@ref InfiniteOpt.Collections.VectorTuple) of the
parameter references that `vref` depends on. This is primarily an internal method
where [`parameter_refs`](@ref parameter_refs(vref::SemiInfiniteVariableRef))
is intended as the preferred user function.
"""
function raw_parameter_refs(vref::SemiInfiniteVariableRef)
    orig_prefs = raw_parameter_refs(infinite_variable_ref(vref))
    eval_supps = eval_supports(vref)
    delete_indices = [!haskey(eval_supps, i) for i in eachindex(orig_prefs)]
    return Collections.restricted_copy(orig_prefs, delete_indices)
end

"""
    parameter_refs(vref::SemiInfiniteVariableRef)::Tuple

Return the infinite parameter references associated with the semi-infinite variable
`vref`. This is formatted as a `Tuple` of containing the parameter references as
they were inputted to define the untranscripted infinite variable except, the
evaluated parameters are excluded.

**Example**
```julia-repl
julia> parameter_refs(vref)
(t, [x[1], x[2]])
```
"""
function parameter_refs(vref::SemiInfiniteVariableRef)
    return Tuple(raw_parameter_refs(vref))
end

"""
    parameter_list(vref::SemiInfiniteVariableRef)::Vector{GeneralVariableRef}

Return a vector of the parameter references that `vref` depends on. This is
primarily an internal method where [`parameter_refs`](@ref parameter_refs(vref::SemiInfiniteVariableRef))
is intended as the preferred user function.
"""
function parameter_list(vref::SemiInfiniteVariableRef)::Vector{GeneralVariableRef}
    orig_prefs = parameter_list(infinite_variable_ref(vref))
    eval_supps = eval_supports(vref)
    return [orig_prefs[i] for i in eachindex(orig_prefs) if !haskey(eval_supps, i)]
end

################################################################################
#                            VARIABLE DEPENDENCIES
################################################################################
# Extend _derivative_dependencies
function _derivative_dependencies(
    vref::SemiInfiniteVariableRef
    )::Vector{DerivativeIndex}
    return _data_object(vref).derivative_indices
end

"""
    used_by_derivative(vref::SemiInfiniteVariableRef)::Bool

Return a `Bool` indicating if `vref` is used by a derivative.

**Example**
```julia-repl
julia> used_by_derivative(vref)
true
```
"""
function used_by_derivative(vref::SemiInfiniteVariableRef)::Bool
    return !isempty(_derivative_dependencies(vref))
end

"""
    is_used(vref::SemiInfiniteVariableRef)::Bool

Return a `Bool` indicating if `vref` is used in the model.

**Example**
```julia-repl
julia> is_used(vref)
true
```
"""
function is_used(vref::SemiInfiniteVariableRef)::Bool
    if used_by_measure(vref) || used_by_constraint(vref)
        return true
    end
    for dindex in _derivative_dependencies(vref)
        if is_used(DerivativeRef(JuMP.owner_model(vref), dindex))
            return true
        end
    end
    return false
end

################################################################################
#                            VARIABLE INFO METHODS
################################################################################
"""
    JuMP.has_lower_bound(vref::SemiInfiniteVariableRef)::Bool

Extend `JuMP.has_lower_bound` to return a `Bool` whether the original
infinite variable of `vref` has a lower bound.

**Example**
```julia-repl
julia> has_lower_bound(vref)
true
```
"""
function JuMP.has_lower_bound(vref::SemiInfiniteVariableRef)::Bool
    return JuMP.has_lower_bound(infinite_variable_ref(vref))
end

"""
    JuMP.lower_bound(vref::SemiInfiniteVariableRef)::Float64

Extend `JuMP.lower_bound` to return the lower bound of the original
infinite variable of `vref`. Errors if `vref` doesn't have a lower bound.

**Example**
```julia-repl
julia> lower_bound(vref)
0.0
```
"""
function JuMP.lower_bound(vref::SemiInfiniteVariableRef)::Float64
    ivref = dispatch_variable_ref(infinite_variable_ref(vref))
    if !JuMP.has_lower_bound(ivref)
        error("Variable $(vref) does not have a lower bound.")
    end
    return JuMP.lower_bound(ivref)
end

# Extend to return the index of the lower bound constraint associated with the
# original infinite variable of `vref`.
function InfiniteOpt._lower_bound_index(
    vref::SemiInfiniteVariableRef
    )::InfOptConstraintIndex
    ivref = dispatch_variable_ref(infinite_variable_ref(vref))
    if !JuMP.has_lower_bound(ivref)
        error("Variable $(vref) does not have a lower bound.")
    end
    return InfiniteOpt._lower_bound_index(ivref)
end

"""
    JuMP.LowerBoundRef(vref::SemiInfiniteVariableRef)::InfOptConstraintRef

Extend `JuMP.LowerBoundRef` to extract a constraint reference for the
lower bound of the original infinite variable of `vref`.

**Example**
```julia-repl
julia> cref = LowerBoundRef(vref)
var >= 0.0
```
"""
function JuMP.LowerBoundRef(
    vref::SemiInfiniteVariableRef
    )::InfOptConstraintRef
    return JuMP.LowerBoundRef(infinite_variable_ref(vref))
end

"""
    JuMP.has_upper_bound(vref::SemiInfiniteVariableRef)::Bool

Extend `JuMP.has_upper_bound` to return a `Bool` whether the original
infinite variable of `vref` has an upper bound.

**Example**
```julia-repl
julia> has_upper_bound(vref)
true
```
"""
function JuMP.has_upper_bound(vref::SemiInfiniteVariableRef)::Bool
    return JuMP.has_upper_bound(infinite_variable_ref(vref))
end

"""
    JuMP.upper_bound(vref::SemiInfiniteVariableRef)::Float64

Extend `JuMP.upper_bound` to return the upper bound of the original
infinite variable of `vref`. Errors if `vref` doesn't have a upper bound.

**Example**
```julia-repl
julia> upper_bound(vref)
0.0
```
"""
function JuMP.upper_bound(vref::SemiInfiniteVariableRef)::Float64
    ivref = dispatch_variable_ref(infinite_variable_ref(vref))
    if !JuMP.has_upper_bound(ivref)
        error("Variable $(vref) does not have a upper bound.")
    end
    return JuMP.upper_bound(ivref)
end

# Extend to return the index of the upper bound constraint associated with the
# original infinite variable of `vref`.
function InfiniteOpt._upper_bound_index(
    vref::SemiInfiniteVariableRef
    )::InfOptConstraintIndex
    ivref = dispatch_variable_ref(infinite_variable_ref(vref))
    if !JuMP.has_upper_bound(ivref)
        error("Variable $(vref) does not have a upper bound.")
    end
    return InfiniteOpt._upper_bound_index(ivref)
end

"""
    JuMP.UpperBoundRef(vref::SemiInfiniteVariableRef)::InfOptConstraintRef

Extend `JuMP.UpperBoundRef` to extract a constraint reference for the
upper bound of the original infinite variable of `vref`.

**Example**
```julia-repl
julia> cref = UpperBoundRef(vref)
var <= 1.0
```
"""
function JuMP.UpperBoundRef(
    vref::SemiInfiniteVariableRef
    )::InfOptConstraintRef
    return JuMP.UpperBoundRef(infinite_variable_ref(vref))
end

"""
    JuMP.is_fixed(vref::SemiInfiniteVariableRef)::Bool

Extend `JuMP.is_fixed` to return `Bool` whether the original infinite
variable of `vref` is fixed.

**Example**
```julia-repl
julia> is_fixed(vref)
true
```
"""
function JuMP.is_fixed(vref::SemiInfiniteVariableRef)::Bool
    return JuMP.is_fixed(infinite_variable_ref(vref))
end

"""
    JuMP.fix_value(vref::SemiInfiniteVariableRef)::Float64

Extend `JuMP.fix_value` to return the fix value of the original infinite
variable of `vref`. Errors if variable is not fixed.

**Example**
```julia-repl
julia> fix_value(vref)
0.0
```
"""
function JuMP.fix_value(vref::SemiInfiniteVariableRef)::Float64
    ivref = dispatch_variable_ref(infinite_variable_ref(vref))
    if !JuMP.is_fixed(ivref)
        error("Variable $(vref) is not fixed.")
    end
    return JuMP.fix_value(ivref)
end

# Extend to return the index of the fix constraint associated with the original
# infinite variable of `vref`.
function InfiniteOpt._fix_index(
    vref::SemiInfiniteVariableRef
    )::InfOptConstraintIndex
    ivref = dispatch_variable_ref(infinite_variable_ref(vref))
    if !JuMP.is_fixed(ivref)
        error("Variable $(vref) is not fixed.")
    end
    return InfiniteOpt._fix_index(ivref)
end

"""
    JuMP.FixRef(vref::SemiInfiniteVariableRef)::InfOptConstraintRef

Extend `JuMP.FixRef` to return the constraint reference of the fix
constraint associated with the original infinite variable of `vref`. Errors
`vref` is not fixed.

**Examples**
```julia-repl
julia> cref = FixRef(vref)
var == 1.0
```
"""
function JuMP.FixRef(
    vref::SemiInfiniteVariableRef
    )::InfOptConstraintRef
    return JuMP.FixRef(infinite_variable_ref(vref))
end

# Fallback dispatch
function JuMP.start_value(vref::SemiInfiniteVariableRef)
    return JuMP.start_value(infinite_variable_ref(vref))
end

"""
    start_value_function(vref::SemiInfiniteVariableRef)::Union{Nothing, Function}

Return the function that is used to generate the start values of `vref` for
particular support values. Returns `nothing` if no start behavior has been
specified.

**Example**
```julia-repl
julia> start_value_func(vref)
my_func
```
"""
function start_value_function(vref::SemiInfiniteVariableRef)
    return start_value_function(infinite_variable_ref(vref))
end

"""
    JuMP.is_binary(vref::SemiInfiniteVariableRef)::Bool

Extend `JuMP.is_binary` to return `Bool` whether the original infinite
variable of `vref` is binary.

**Example**
```julia-repl
julia> is_binary(vref)
true
```
"""
function JuMP.is_binary(vref::SemiInfiniteVariableRef)::Bool
    return JuMP.is_binary(infinite_variable_ref(vref))
end

# Extend to return the index of the binary constraint associated with the
# original infinite variable of `vref`.
function InfiniteOpt._binary_index(
    vref::SemiInfiniteVariableRef
    )::InfOptConstraintIndex
    ivref = dispatch_variable_ref(infinite_variable_ref(vref))
    if !JuMP.is_binary(ivref)
        error("Variable $(vref) is not binary.")
    end
    return InfiniteOpt._binary_index(ivref)
end

"""
    JuMP.BinaryRef(vref::SemiInfiniteVariableRef)::InfOptConstraintRef

Extend `JuMP.BinaryRef` to return a constraint reference to the
constraint constrainting the original infinite variable of `vref` to be binary.
Errors if one does not exist.

**Example**
```julia-repl
julia> cref = BinaryRef(vref)
var binary
```
"""
function JuMP.BinaryRef(
    vref::SemiInfiniteVariableRef
    )::InfOptConstraintRef
    return JuMP.BinaryRef(infinite_variable_ref(vref))
end

"""
    JuMP.is_integer(vref::SemiInfiniteVariableRef)::Bool

Extend `JuMP.is_integer` to return `Bool` whether the original infinite
variable of `vref` is integer.

**Example**
```julia-repl
julia> is_integer(vref)
true
```
"""
function JuMP.is_integer(vref::SemiInfiniteVariableRef)::Bool
    return JuMP.is_integer(infinite_variable_ref(vref))
end

# Extend to return the index of the integer constraint associated with the
# original infinite variable of `vref`.
function InfiniteOpt._integer_index(
    vref::SemiInfiniteVariableRef
    )::InfOptConstraintIndex
    ivref = dispatch_variable_ref(infinite_variable_ref(vref))
    if !JuMP.is_integer(ivref)
        error("Variable $(vref) is not an integer.")
    end
    return InfiniteOpt._integer_index(ivref)
end

"""
    JuMP.IntegerRef(vref::SemiInfiniteVariableRef)::InfOptConstraintRef

Extend `JuMP.IntegerRef` to return a constraint reference to the
constraint constrainting the original infinite variable of `vref` to be integer.
Errors if one does not exist.

**Example**
```julia-repl
julia> cref = IntegerRef(vref)
var integer
```
"""
function JuMP.IntegerRef(
    vref::SemiInfiniteVariableRef
    )::InfOptConstraintRef
    return JuMP.IntegerRef(infinite_variable_ref(vref))
end

################################################################################
#                                  DELETION
################################################################################
# Extend _delete_variable_dependencies (for use with JuMP.delete)
function _delete_variable_dependencies(vref::SemiInfiniteVariableRef)::Nothing
    # remove mapping to infinite variable
    ivref = infinite_variable_ref(vref)
    filter!(e -> e != JuMP.index(vref), _semi_infinite_variable_dependencies(ivref))
    # delete associated derivative variables and mapping
    model = JuMP.owner_model(vref)
    for index in _derivative_dependencies(vref)
        JuMP.delete(model, dispatch_variable_ref(model, index))
    end
    # remove the lookup entry
    delete!(JuMP.owner_model(vref).semi_lookup, (ivref, eval_supports(vref)))
    return
end


text:
#=
gams:
- Julia version: 1.3.1
- Author: Berk
- Date: 2020-12-10
=#

"""Returns all constants from GAMS Model. """
function constants(gams::Dict{String, Any})
    consts = Dict{String, Any}(); #
    for p in ("parameters", "tables")
        if haskey(gams, p)
            for (k, v) in gams[p]
                v = v isa Ref ? v[] : v
                consts[GAMSFiles.getname(k)] = v;
            end
        end
    end
    return consts
end

""" Moves GAMS constants into the global scope.
    NOTE: This is only for debugging purposes, could cause bugs.
"""
function constants_to_global(gams::Dict{String, Any})
    consts = constants(gams)
    for (key, value) in consts
        eval(Meta.parse("$(key) = $(value)"))
    end
end

ext: jl


text:
using Combinatorics: combinations
using Random

"""
    total_combinations(arity, n_way)

Given an array of the number of values for each parameter and the
level of coverage, return the number of n-tuples required for complete
coverage.
"""
function total_combinations(arity, n_way)
    param_cnt = length(arity)
    sum(prod(arity[key_set]) for key_set in combinations(1:param_cnt, n_way))
end

"""
    next_multiplicative!(values, arity)

Given a set of values for parameters, this returns the next possible
value. When it reaches the end, it cycles back to the start. The first
value is all ones. The last value equals the arity.
"""
function next_multiplicative!(values, arity)
    carry = 1
    for slot_idx in length(values):-1:1
        values[slot_idx] += carry
        if arity[slot_idx] < values[slot_idx]
            values[slot_idx] = 1
            carry = 1
        else
            carry = 0
        end
    end
end

"""
    all_combinations(arity, n_way)

This represents possible coverage as a matrix, one column per parameter,
zero if not used. `arity` is a list of the number of values for each parameter,
and `n_way` is the order of the combinations, most commonly 2-way.
"""
function all_combinations(arity, n_way)
    v_cnt = length(arity)
    indices = collect(combinations(1:v_cnt, n_way))
    combinations_cnt = total_combinations(arity, n_way)

    coverage = zeros(Int, combinations_cnt, v_cnt)
    idx = 1
    for indices_idx in 1:size(indices, 1)
        offset = indices[indices_idx]
        sub_arity = arity[offset]
        sub_cnt = prod(sub_arity)
        values = copy(sub_arity)
        for sub_idx in 1:sub_cnt
            next_multiplicative!(values, sub_arity)
            coverage[idx, offset] = values
            idx += 1
        end
    end
    coverage
end

"""
The wayness is a dictionary from an integer, the wayness, to a set
of lists of indices that should have that wayness together. The `base_wayness`
is the `n_way` for the rest of the variables. It should be less than
the other sets of waynesses.
"""
function multi_way_coverage(arity, wayness, base_wayness)
    orders = sort(collect(keys(wayness)), rev = true)
    @assert minimum(orders) > base_wayness
    @assert maximum(orders) <= length(arity)

    param_cnt = length(arity)
    order_combos = Any[]
    for order in orders
        # The parameter set is a list of parameter indices.
        for param_set in wayness[order]
            high_combos = all_combinations(arity[param_set], order)
            widened = zeros(Int, size(high_combos, 1), param_cnt)
            widened[:, param_set] = high_combos
            push!(order_combos, widened)
        end
        # We could remove duplicates of higher orders.
    end
    push!(order_combos, all_combinations(arity, base_wayness))
    vcat(order_combos...)
end

function combination_histogram(allc, arity)
    width = maximum(arity)
    hist = zeros(Int, length(arity), width)
    for row_idx in 1:size(allc, 1)
        for col_idx in 1:width
            if allc[row_idx, col_idx] > 0
                hist[col_idx, allc[row_idx, col_idx]] += 1
            end
        end
    end
    hist
end

function most_to_cover(allc, row_cnt)
   argmax(vec(sum(allc[1:row_cnt, :] .!= 0, dims = 1)))
end

"""
    coverage_by_parameter(allc, row_cnt)

Given a coverage array and a the number of rows to search,
find how many times a nonzero appears in each column.
"""
function coverage_by_parameter(allc, row_cnt)
    vec(sum(allc[1:row_cnt, :] .!= 0, dims = 1))
end

"""
    coverage_by_value(allc, row_cnt, arity, param_idx)

Given a particular parameter, look at that column of the coverage matrix
and return a histogram of how many times each value appears in
that column.
"""
function coverage_by_value(allc, row_cnt, arity, param_idx)
    hist = zeros(Int, arity[param_idx] + 1)
    for row_idx in 1:row_cnt
        hist[allc[row_idx, param_idx] + 1] += 1
    end
    hist[2:end]
end

function most_common_value(allc, row_cnt, arity, param_idx)
    hist = zeros(Int, arity[param_idx] + 1)
    for row_idx in 1:row_cnt
        hist[allc[row_idx, param_idx] + 1] += 1
    end
    argmax(hist[2:end])
end

"""
    most_matches_existing(allc, row_cnt, arity, existing, param_idx)

There is a coverage array, and there is a list of existing choices for the
values of parameters. This asks, given another parameter, how many matches
would each of its values produce, were it chosen. The return value is a
vector where the first entry is the number of matches for value = 1, the second
entry is number of matches for value = 2, and so on.
"""
function most_matches_existing(allc, row_cnt, arity, existing, param_idx)
    @assert existing[param_idx] == 0
    param_cnt = length(arity)
    # The match count can only be as large as n_way - 1 because the
    # column with the parameter is non-zero and will be zero in existing.
    # If n-way is 3 and we only know one parameter so far, that's another limit
    # to the possible match size.
    max_known = sum(existing .!= 0)
    # params_known = min(sum(existing .!= 0), n_way - 1)
    hist = zeros(Int, arity[param_idx])
    for row_idx in 1:row_cnt
        # The given parameter column is part of this match.
        if allc[row_idx, param_idx] != 0
            match_cnt = 0
            n_way = 0
            for match_idx in 1:param_cnt
                if allc[row_idx, match_idx] != 0
                    n_way += 1
                end
                if existing[match_idx] != 0 && allc[row_idx, match_idx] == existing[match_idx]
                    match_cnt += 1
                end
            end
            if match_cnt == min(max_known, n_way)
                hist[allc[row_idx, param_idx]] += 1
            end
        end
    end
    hist
end

combination_number(n, m) = prod(n:-1:(n-m+1)) ÷ factorial(m)

function pairs_in_entry(entry, n_way)
    n = length(entry)
    ans = zeros(Int, combination_number(n, n_way), n)
    row_idx = 1
    for param_idx in combinations(1:length(entry), n_way)
        for col_idx in 1:n_way
            ans[row_idx, param_idx[col_idx]] = entry[param_idx[col_idx]]
        end
        row_idx += 1
    end
    ans
end

"""
    tuples_in_trials(trials, n_way)

Find the number of distinct tuples in a set of trials. This uses a different method
to count and store those tuples, as a check on other methods. It returns a dictionary
where each key is the index of non-zero parameters and the value is the set of all
ways those values appear.
"""
function tuples_in_trials(trials, n_way)
    param_cnt = length(trials[1])
    param_combo = [Tuple(x) for x in combinations(1:param_cnt, n_way)]
    # The dictionary with keys = parameter indices and dictionary
    # value is set of tuples of parameter values.
    seen = Dict(x => Set{NTuple{n_way,Int}}() for x in param_combo)
    for trial_idx in 1:length(trials)
        entry = trials[trial_idx]
        for params in param_combo
            push!(seen[params], Tuple(entry[[x for x in params]]))
        end
    end
    seen
end

function coverage_by_tuple(trials, n_way)
    accumulated = tuples_in_trials(trials, n_way)
    sum([length(aset) for aset in values(accumulated)])
end

"""
    add_coverage!(allc, row_cnt, entry)

The coverage matrix, `allc`, has `row_cnt` entries that are considered
uncovered, and the rest have been covered already. This adds a new entry,
which is a complete set of parameter choices. For every parameter that's
covered, this moves those to the end. It returns a new `row_cnt` which is
the number of initial uncovered rows.
"""
function add_coverage!(allc, row_cnt, entry)
    param_cnt = length(entry)

    covers = zeros(Int, row_cnt)
    cover_cnt = 0
    # Find matches before reordering them to the end.
    for row_idx in 1:row_cnt
        match_cnt = 0
        n_way = 0
        for match_idx in 1:param_cnt
            if allc[row_idx, match_idx] != 0
                n_way += 1
            end
            if allc[row_idx, match_idx] == entry[match_idx]
                match_cnt += 1
            end
        end
        if match_cnt == n_way
            cover_cnt += 1
            covers[cover_cnt] = row_idx
        end
    end
    # Given the matches, we can swap them to the end of the matrix.
    # Work from the end in case a match is near row_cnt.
    for cover_idx in cover_cnt:-1:1
        if row_cnt > 1
            save = allc[row_cnt, :]
            allc[row_cnt, :] = allc[covers[cover_idx], :]
            allc[covers[cover_idx], :] = save
            row_cnt -= 1
        else
            row_cnt -= 1
        end
    end
    row_cnt
end

"""
    match_score(allc, row_cnt, entry)

If we were to add this entry to the trials, how many uncovered
n-tuples would it now cover? `allc` is the matrix of n-tuples.
`row_cnt` is the first set of rows, representing uncovered tuples.
`n_way` is the length of each tuple. Entry is a set of putative
paramter values. Returns an integer number of newly-covered tuples.

This is a place to look into alternative algorithms. We could weight
the match score by increasing the score for greater wayness.
"""
function match_score(allc, row_cnt, entry)
    param_cnt = length(entry)
    cover_cnt = 0
    for row_idx in 1:row_cnt
        param_match_cnt = 0
        n_way = 0
        for match_idx in 1:param_cnt
            if allc[row_idx, match_idx] == entry[match_idx]
                param_match_cnt += 1
            end
            if allc[row_idx, match_idx] != 0
                n_way += 1
            end
        end
        if param_match_cnt == n_way
            cover_cnt += 1
        end
    end
    cover_cnt
end

"""
    argmin_rand(rng, v)

Given a vector, find the index of the smallest value. If more than
one value is the smallest, then randomly choose among the smallest
values.
"""
function argmin_rand(rng, v)
    small = typemax(v[1])
    small_extra_cnt = 0
    small_idx = -1
    for i in 1:length(v)
        if v[i] < small
            small = v[i]
            small_extra_cnt = 0
            small_idx = i
        elseif v[i] == small
            small_extra_cnt += 1
        # else not the smallest.
        end
    end
    if small_extra_cnt == 0
        return small_idx
    else
        which = rand(rng, 1:(small_extra_cnt + 1))
        which_cnt = 1
        for s_idx in small_idx:length(v)
            if v[s_idx] == small
                if which_cnt == which
                    return s_idx
                end
                which_cnt += 1
            end
        end
    end
    return 0
end

function n_way_coverage(arity, n_way, M, rng)
    param_cnt = length(arity)
    allc = all_combinations(arity, n_way)
    # Every combination is nonzero.
    @assert sum(sum(allc, dims = 2) == 0) == 0
    remain = size(allc, 1)
    maximum_match_score = combination_number(param_cnt, n_way)
    param_cnt = length(arity)

    # Array of arrays.
    coverage = Array{Array{Int64,1},1}()

    trials = zeros(Int, M, param_cnt)
    trial_scores = zeros(Int, M)
    params = zeros(Int, param_cnt)
    entry = zeros(Int, param_cnt)

    loop_idx = 1
    while remain > 0
        params[:] = 1:param_cnt
        param_coverage = coverage_by_parameter(allc, remain)
        params[1] = argmin_rand(rng, -param_coverage)
        params[params[1]] = 1
        for trial_idx in 1:M
            params[2:end] = shuffle(rng, params[2:end])
            entry[:] .= 0
            candidate_params = coverage_by_value(allc, remain, arity, params[1])
            entry[params[1]] = argmin_rand(rng, -candidate_params)
            for p_idx in 2:param_cnt
                candidate_values = most_matches_existing(allc, remain, arity, entry, params[p_idx])
                entry[params[p_idx]] = argmin_rand(rng, -candidate_values)
            end
            score = match_score(allc, remain, entry)
            trial_scores[trial_idx] = score
            trials[trial_idx, :] = entry
        end
        chosen_idx = argmin_rand(rng, -trial_scores)
        maximum_score = trial_scores[chosen_idx]
        if maximum_score > 0
            chosen_trial = trials[chosen_idx, :]
            remain = add_coverage!(allc, remain, chosen_trial)
            push!(coverage, chosen_trial)
        else
            println("error because nothing was covered")
        end
        loop_idx += 1
    end
    coverage
end

function n_way_coverage_init(arity, n_way, seed, M, rng)
    param_cnt = length(arity)
    allc = all_combinations(arity, n_way)
    # Every combination is nonzero.
    @assert sum(sum(allc, dims = 2) == 0) == 0
    remain = size(allc, 1)
    maximum_match_score = combination_number(param_cnt, n_way)
    param_cnt = length(arity)

    # Array of arrays.
    coverage = Array{Array{Int64,1},1}()

    # Don't save the seeds to the coverage array because the user knows
    # what they are.
    for seed_idx in 1:size(seed, 1)
        remain = add_coverage!(allc, remain, seed[seed_idx, :])
    end

    trials = zeros(Int, M, param_cnt)
    trial_scores = zeros(Int, M)
    params = zeros(Int, param_cnt)
    entry = zeros(Int, param_cnt)

    exhausted_by_filter = false
    loop_idx = 1
    while remain > 0 || exhausted_by_filter
        params[:] = 1:param_cnt
        param_coverage = coverage_by_parameter(allc, remain)
        params[1] = argmin_rand(rng, -param_coverage)
        params[params[1]] = 1
        for trial_idx in 1:M
            params[2:end] = shuffle(rng, params[2:end])
            entry[:] .= 0
            candidate_params = coverage_by_value(allc, remain, arity, params[1])
            entry[params[1]] = argmin_rand(rng, -candidate_params)
            for p_idx in 2:param_cnt
                candidate_values = most_matches_existing(allc, remain, arity, entry, params[p_idx])
                entry[params[p_idx]] = argmin_rand(rng, -candidate_values)
            end
            score = match_score(allc, remain, entry)
            trial_scores[trial_idx] = score
            trials[trial_idx, :] = entry
        end
        chosen_idx = argmin_rand(rng, -trial_scores)
        maximum_score = trial_scores[chosen_idx]
        if maximum_score > 0
            chosen_trial = trials[chosen_idx, :]
            remain = add_coverage!(allc, remain, chosen_trial)
            push!(coverage, chosen_trial)
        else
            exhausted_by_filter = true
        end
        loop_idx += 1
    end
    coverage
end

function remove_combinations(allc, disallow)
    allow_cnt = 0
    allowed = zeros(Int, size(allc, 1))
    for i in 1:size(allc, 1)
        if !disallow(allc[i, :])
            allow_cnt += 1
            allowed[allow_cnt] = i
        # else disallowed
        end
    end
    allc[allowed[1:allow_cnt], :]
end

function arguments_to_arity(arguments)
    [length(x) for x in arguments]
end

function indices_to_arguments(arguments, indices)
    Tuple(arguments[i[1]][i[2]] for i in zip(1:length(indices), indices))
end

function n_way_coverage_filter(arity, n_way, disallow, seed, M, rng)
    param_cnt = length(arity)
    allc = all_combinations(arity, n_way)
    # Every combination is nonzero.
    @assert sum(sum(allc, dims = 2) == 0) == 0
    before_disallow = size(allc, 1)
    allc = remove_combinations(allc, disallow)
    remain = size(allc, 1)
    maximum_match_score = combination_number(param_cnt, n_way)
    param_cnt = length(arity)

    # Array of arrays.
    coverage = Array{Array{Int64,1},1}()

    # Don't save the seeds to the coverage array because the user knows
    # what they are.
    for seed_idx in 1:size(seed, 1)
        remain = add_coverage!(allc, remain, seed[seed_idx, :])
    end

    trials = zeros(Int, M, param_cnt)
    trial_scores = zeros(Int, M)
    params = zeros(Int, param_cnt)
    entry = zeros(Int, param_cnt)

    loop_idx = 1
    while remain > 0
        params[:] = 1:param_cnt
        param_coverage = coverage_by_parameter(allc, remain)
        params[1] = argmin_rand(rng, -param_coverage)
        params[params[1]] = 1
        trial_cnt = 0
        while trial_cnt < M
            params[2:end] = shuffle(rng, params[2:end])
            entry[:] .= 0
            candidate_params = coverage_by_value(allc, remain, arity, params[1])
            entry[params[1]] = argmin_rand(rng, -candidate_params)
            for p_idx in 2:param_cnt
                candidate_values = most_matches_existing(allc, remain, arity, entry, params[p_idx])
                entry[params[p_idx]] = argmin_rand(rng, -candidate_values)
            end
            if !disallow(entry)
                trial_cnt += 1
                score = match_score(allc, remain, entry)
                trial_scores[trial_cnt] = score
                trials[trial_cnt, :] = entry
            end
        end
        chosen_idx = argmin_rand(rng, -trial_scores[1:trial_cnt])
        chosen_trial = trials[chosen_idx, :]
        maximum_score = trial_scores[chosen_idx]
        if 0 < maximum_score
            remain = add_coverage!(allc, remain, chosen_trial)
            push!(coverage, chosen_trial)
        # else Failure to cover can happen for a bad draw in the shuffle.
        end
        loop_idx += 1
    end
    coverage
end

function n_way_coverage_multi(arity, allc, disallow, seed, M, rng)
    param_cnt = length(arity)
    remain = size(allc, 1)
    before_disallow = size(allc, 1)
    allc = remove_combinations(allc, disallow)

    # Array of arrays.
    coverage = Array{Array{Int64,1},1}()

    # Don't save the seeds to the coverage array because the user knows
    # what they are.
    for seed_idx in 1:size(seed, 1)
        remain = add_coverage!(allc, remain, seed[seed_idx, :])
    end

    trials = zeros(Int, M, param_cnt)
    trial_scores = zeros(Int, M)
    params = zeros(Int, param_cnt)
    entry = zeros(Int, param_cnt)

    loop_idx = 1
    while remain > 0
        params[:] = 1:param_cnt
        param_coverage = coverage_by_parameter(allc, remain)
        params[1] = argmin_rand(rng, -param_coverage)
        params[params[1]] = 1
        trial_cnt = 0
        while trial_cnt < M
            params[2:end] = shuffle(rng, params[2:end])
            entry[:] .= 0
            candidate_params = coverage_by_value(allc, remain, arity, params[1])
            entry[params[1]] = argmin_rand(rng, -candidate_params)
            for p_idx in 2:param_cnt
                candidate_values = most_matches_existing(allc, remain, arity, entry, params[p_idx])
                entry[params[p_idx]] = argmin_rand(rng, -candidate_values)
            end
            if !disallow(entry)
                trial_cnt += 1
                score = match_score(allc, remain, entry)
                trial_scores[trial_cnt] = score
                trials[trial_cnt, :] = entry
            end
        end
        chosen_idx = argmin_rand(rng, -trial_scores[1:trial_cnt])
        chosen_trial = trials[chosen_idx, :]
        maximum_score = trial_scores[chosen_idx]
        if 0 < maximum_score
            remain = add_coverage!(allc, remain, chosen_trial)
            push!(coverage, chosen_trial)
        # else Failure to cover can happen for a bad draw in the shuffle.
        end
        loop_idx += 1
    end
    coverage
end

max_stars_repo_name: adolgert/test-scientific-code

max_issues_repo_name: adolgert/test-scientific-code

max_forks_repo_name: adolgert/test-scientific-code


text:
using CSV
using DataFrames
using DelimitedFiles
using SnpArrays
using Random

"""
Generate numeric matrix of `x` where columns are groups.

# Inputs
- `x`: A `SnpArray` that contain original snps and its knockoffs
- `groups`: A `Vector{Int}`. `groups[i]` is the group membership of snp `i`
- `knockoff` A `BitVector`. `knockoff[i] == true` means snp `i` is the knockoff

# Outputs
-`X`: numeric matrix of `x` where columns are groups
"""
function make_X_by_group(x::SnpArray, groups::Vector{Int}, knockoff::BitVector)
    n, p = size(x)
    length(groups) == p || error("Number of SNPs in `x` has different length than the `groups` vector")
    X = zeros(n, length(unique(groups)))
    snp = zeros(n)
    # loop over snp by snp
    @inbounds for i in 1:p
        copyto!(snp, @view(x[:, i]))
        knockoff[i] && continue # skip knockoffs if requested
        X[:, groups[i]] .+= snp
    end
    return X
end

"""
Generate numeric matrix of `x` where columns are groups and their knockoffs.

# Outputs
- `X`: numeric matrix of `x` where columns are groups and their knockoffs
- `original_list`: `original_list[i]` is the column of `X` that stores the `i`th group
- `knockoff_list`: `knockoff_list[i]` is the column of `X` that stores the knockoff
    for the `i`th group.
"""
function make_Xko_by_group(x::SnpArray, groups::Vector{Int}, knockoff::BitVector)
    n, p = size(x)
    length(groups) == p || error("Number of SNPs in `x` has different length than the `groups` vector")
    P = length(unique(groups))
    X = zeros(n, 2P)
    snp = zeros(n)
    # loop over snp by snp, adding knockoffs to index 2g and originals to 2g-1
    @inbounds for i in 1:p
        copyto!(snp, @view(x[:, i]))
        col = knockoff[i] ? 2groups[i] : 2groups[i] - 1
        X[:, col] .+= snp
    end
    # permute knockoffs and originals with probability 0.5
    original_list = [2i - 1 for i in 1:P]
    knockoff_list = [2i for i in 1:P]
    for j in 1:P
        if rand() < 0.5
            col1, col2 = 2j - 1, 2j
            @inbounds for k in 1:n
                X[k, col1], X[k, col2] = X[k, col2], X[k, col1]
            end
            original_list[j], knockoff_list[j] = knockoff_list[j], original_list[j]
        end
    end
    return X, original_list, knockoff_list
end

#
# import key and data
#
chr = 10
keyfile = "/scratch/users/bbchu/ukb/groups/Radj20_K50_s0/ukb_gen_chr$chr.key"
df = CSV.read(keyfile, DataFrame)
groups = df[!, :Group]
knockoff = convert(BitVector, df[!, :Knockoff])
x = SnpArray("/scratch/users/bbchu/ukb/groups/Radj20_K50_s0/ukb_gen_chr$chr.bed")

#
# generate grouped X (without knockoff) and save
#
# @time X = make_X_by_group(x, groups, knockoff)
# @time writedlm("chr$chr.noKnockoffs.data.txt", X)

#
# generate grouped X (includes knockoff) and save
#
# Random.seed!(2021)
# @time Xko, original_list, knockoff_list = make_Xko_by_group(x, groups, knockoff)
# @time writedlm("chr$chr.data.txt", Xko)
# @time writedlm("chr$chr.original.txt", Xko)
# @time writedlm("chr$chr.knockoff.txt", Xko)


text:
type Body
  mass::Float64
  pos::Vector{Float64}
  vel::Vector{Float64}

  id::Integer

  function Body(;mass=0.0, pos=Vector{Float64}, vel=Vector{Float64}, id=0)
    # accel_hist_len = 4
    # accel_hist = zeros(Float64, NDIMS, accel_hist_len)

    new(mass, pos, vel, id)
  end
end

function Body(d::Dict{String, Any})
  b = Body(
    mass = d["mass"],
    pos = d["pos"],
    vel = d["vel"]
  )
  b.id = if haskey(d, "id")
    d["id"]
  else
    0
  end
  b
end

######## IO ###########

import Base.show
function show(io::IO, a::Body)
  print(io, "mass:  ", a.mass, "\n")
  print(io, "pos:   ", join(a.pos, ", "), "\n")
  print(io, "vel:   ", join(a.vel, ", "), "\n")
end

# function pp(io::IO, a::Body, nb::NBodySystem)
#   acc = accel(a, nb)
#   show(io, a)
#   print(io, "accel: ", join(acc, ", "), "\n")
# end
#
# function pp(a::Body, nb::NBodySystem)
#   pp(STDOUT, a, nb)
# end

# import Base.print
# function print(io::IO, a::Body)
#   @printf io "%24.16e\n" a.mass
#
#   map(Array[a.pos, a.vel]) do b
#     map(x-> @printf(io,"%24.16e",x), b)
#     print(io, "\n")
#   end
# end

function read_body()
  mass = parse(Float64, readline())
  pos = map(x -> parse(Float64, x), split(readline()))
  vel = map(x -> parse(Float64, x), split(readline()))
  return Body(mass, pos, vel)
end

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

function accel(body::Body, bodies::Array{Body}, soft_len)
  acc = zeros(body.vel)
  for a in bodies
    if a != body
      r = a.pos - body.pos
      r2 = dot(r,r) + soft_len*soft_len
      r3 = r2*sqrt(r2)
      acc += r*(a.mass/r3)
    end
  end
  acc
end

function jerk(body::Body, bodies::Array{Body}, soft_len)
  j = zeros(body.vel)
  for a in bodies
    if a != body
      r = a.pos - body.pos
      r2 = dot(r,r) + soft_len*soft_len
      r3 = r2*sqrt(r2)
      v = a.vel - body.vel
      j += (v-r*(3*dot(r,v)/r2))*(a.mass/r3)
    end
  end
  j
end

function collision_time_scale(body::Body, bodies::Array{Body})
  time_scale_sq = Inf
  for b in bodies
    if b != body
      r = b.pos - body.pos
      v = b.vel - body.vel

      r2 = dot(r, r)
      v2 = dot(v, v)

      estimate_sq = r2 / v2

      if time_scale_sq > estimate_sq
        time_scale_sq = estimate_sq
      end

      a = (body.mass + b.mass) / r2
      estimate_sq = sqrt(r2) / a

      if time_scale_sq > estimate_sq
        time_scale_sq = estimate_sq
      end
    end
  end

  sqrt(time_scale_sq)
end

####### Energy #########

function kin_energy(body::Body)
  return 0.5 * body.mass * dot(body.vel, body.vel)
end

function pot_energy(body::Body, bodies, soft_len)
  p = 0
  for a in bodies
    if a != body
      r = a.pos - body.pos
      p += -body.mass*a.mass/sqrt(dot(r,r) + soft_len*soft_len)
    end
  end
  return p
end


text:
include("../src/linalg_module.jl")

module lutest
using ..linalg
using BenchmarkTools
using LinearAlgebra
using PyPlot
PyPlot.pygui(true)

function test_linsolve_lu()

    #=
    A = [
        1   4   -1;
        3   0   5;
        2   2   1.0
        ]
    b = [10.0, 2, 3]
    =#

    A = [
        3e-6    2       1;
        2       2       2;
        1       2       -1
        ]
    b = [3+3e-6, 6, 2]

    n_dim = 1000
    A = rand(n_dim, n_dim)
    b = rand(n_dim)

    x_lu = linalg.linsolve_lu(A, b)
    @time x_theo = A \ b

    println("x_lu:")
    # println(x_lu)
    println( "residuum norm x_lu:  $(norm( A * x_lu - b )) " )

    println("x_theo:")
    # println(x_theo)
    println( "residuum norm x_theo:  $(norm( A * x_theo - b )) " )

end

function test_lu_decomposition()
    A = [
        1   4   -1;
        3   0   5;
        2   2   1.0
        ]
    b = [10.0, 2, 3]

    (L, U, p) = linalg.lu_decomp_partialpivot(A)

    println("L:")
    show(stdout, MIME("text/plain"), L)
    println("\n")

    println("U:")
    show(stdout, MIME("text/plain"), U)
    println("\n")

    println("p:")
    show(stdout, MIME("text/plain"), p)
    println("\n")

    return (L, U, p)
end

function benchmark_linsolve_lu(A, b)
    b_ = @benchmark linalg.linsolve_lu($A, $b)

    println("=== benchmark_linsolve_lu ===")
    show(
            stdout,
            MIME("text/plain"),
            b_ )
    println("\n\n")
end

function profiler_linsolve_lu()

    N = 3_000
    A = rand(N, N)
    A = diagm( diag(A) ) .+ A
    b = rand(N)

    x = linalg.linsolve_lu(A, b)

end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
struct LSQOutput
    β::Vector{Float64}
    ζ::Vector{Float64}
    Σβ::Matrix{Float64}
    Σζ::Matrix{Float64}
    Σβζ::Matrix{Float64}
    χ²::Float64
    dof::Int
    βindex::Dict{Symbol, Union{Int, UnitRange{Int}}}
    zindex::Dict{Symbol, Union{Int, UnitRange{Int}}}
end

_names_ordered(index) = map(x->x[1], sort(collect(index), by=x->x[2]))

"""Return names of measured parameters `z`."""
names_z(output::LSQOutput) = _names_ordered(output.zindex)

"""Return names of unknown parameters `β`."""
names_β(output::LSQOutput) = _names_ordered(output.βindex)

"""Return standard uncertainty of unknown parameters `β`."""
uncertainty_β(output::LSQOutput) = [sqrt(output.Σβ[i, i]) for i in 1:length(output.β)]

"""Return covariance matrix of unknown parameters `β`."""
function covariance_β(output::LSQOutput)
    names = names_β(output)
    return NamedArray(output.Σβ, (names, names))
end

"""Return correlation matrix of unknown parameters `β`."""
function correlation_β(output::LSQOutput)
    cov = covariance_β(output)
    N = size(cov)[1]
    corr = similar(cov)
    for i in 1:N
        for j in 1:N
            corr[i, j] = cov[i, j]/√(cov[i, i]*cov[j, j])
        end
    end
    return corr
end

"""Display covariance matrix of unknown parameters `β`."""
function print_cov_β(output::LSQOutput)
    pretty_table(
        covariance_β(output);
        row_names=names_β(output),
        header=names_β(output),
        title="covariance of β parameters"
    )
end

"""Display correlation matrix of unknown parameters `β`."""
function print_corr_β(output::LSQOutput)
    pretty_table(
        correlation_β(output);
        row_names=names_β(output),
        header=names_β(output),
        title="correlation of β parameters"
    )
end


text:
module SimpleTrees

using ..Factories: Inference, Argument, Parameter, FunctionFactory, TypeFactory
using ..Factories: addfields!, addparams!, addargs!, addwhereparams!, extendbody!
using ..Factories: MixEscaped, Escaped, UnEscaped, rawexpr
using ..TypeTraits: efficientoperations

export simpletreedepth, simpletreewidth
export AbstractSimpleTree
export root, parent, children
export addnode!, deletenode!
export isleaf, level
export ancestor, descendants, siblings, leaves
export subtree, move!
export SimpleTreeCore, @simpletree, SimpleTree

abstract type SimpleTreeIteration end
struct SimpleTreeDepth <: SimpleTreeIteration end
struct SimpleTreeWidth <: SimpleTreeIteration end

"""
    simpletreedepth

Indicate that the iteration over a tree is depth-first.
"""
const simpletreedepth = SimpleTreeDepth()

"""
    simpletreewidth

Indicate that the iteration over a tree is width-first.
"""
const simpletreewidth = SimpleTreeWidth()

"""
    AbstractSimpleTree{Node, Data}

Abstract type for all concrete trees.
"""
abstract type AbstractSimpleTree{N, D} end

"""
    eltype(tree::AbstractSimpleTree)
    eltype(::Type{<:AbstractSimpleTree{N, D}}) where {N, D}

Get the eltype of a tree.
"""
Base.eltype(tree::AbstractSimpleTree) = eltype(typeof(tree))
Base.eltype(::Type{<:AbstractSimpleTree{N, D}}) where {N, D} = Pair{N, D}

"""
    root(tree::AbstractSimpleTree) -> Union{keytype(tree), Nothing}

Get a tree's root node.
"""
root(tree::AbstractSimpleTree) = tree.TREECORE.root

"""
    haskey(tree::AbstractSimpleTree{N}, node::N) where N -> Bool

Check whether a node is in a tree.
"""
Base.haskey(tree::AbstractSimpleTree{N}, node::N) where N = haskey(tree.TREECORE.contents, node)

"""
    length(tree::AbstractSimpleTree) -> Int

Get the number of a tree's nodes.
"""
Base.length(tree::AbstractSimpleTree) = length(tree.TREECORE.contents)

"""
    parent(tree::AbstractSimpleTree{N}, node::N, superparent::Union{N, Nothing}=nothing) where N -> Union{N, Nothing}

Get the parent of a tree's node. When `node` is the tree's root, return `superparent`.
"""
parent(tree::AbstractSimpleTree{N}, node::N, superparent::Union{N, Nothing}=nothing) where {N} = (node == root(tree) ? superparent : tree.TREECORE.parent[node])

"""
    children(tree::AbstractSimpleTree) -> Vector{keytype(tree)}
    children(tree::AbstractSimpleTree, ::Nothing) -> Vector{keytype(tree)}
    children(tree::AbstractSimpleTree{N}, node::N) where N -> Vector{N}

Get the children of a tree's node.
"""
children(tree::AbstractSimpleTree) = children(tree, nothing)
children(tree::AbstractSimpleTree, ::Nothing) = (root(tree) === nothing) ? error("children error: empty tree!") : [root(tree)]
children(tree::AbstractSimpleTree{N}, node::N) where N = tree.TREECORE.children[node]

"""
    addnode!(tree::AbstractSimpleTree{N}, node::N) where N} -> typeof(tree)
    addnode!(tree::AbstractSimpleTree{N}, ::Nothing, node::N) where N -> typeof(tree)
    addnode!(tree::AbstractSimpleTree{N}, parent::N, node::N) where N -> typeof(tree)

Update the structure of a tree by adding a node. When the parent is `nothing`, the input tree must be empty and the input node becomes the tree's root.
"""
addnode!(tree::AbstractSimpleTree{N}, node::N) where {N} = addnode!(tree, nothing, node)
function addnode!(tree::AbstractSimpleTree{N}, ::Nothing, node::N) where N
    @assert root(tree) === nothing "addnode! error: not empty tree."
    tree.TREECORE.root = node
    tree.TREECORE.children[node] = N[]
    return tree
end
function addnode!(tree::AbstractSimpleTree{N}, parent::N, node::N) where N
    @assert haskey(tree, parent) "addnode! error: parent($parent) not in tree."
    @assert !haskey(tree, node) "addnode! error: node($node) already in tree."
    push!(tree.TREECORE.children[parent], node)
    tree.TREECORE.parent[node] = parent
    tree.TREECORE.children[node] = N[]
    return tree
end

"""
    deletenode!(tree::AbstractSimpleTree{N}, node::N) where N -> typeof(tree)

Update the structure of a tree by deleting a node.
"""
function deletenode!(tree::AbstractSimpleTree{N}, node::N) where N
    @assert haskey(tree, node) "deletenode! error: node($node) not in tree."
    if node == root(tree)
        tree.TREECORE.root = nothing
    else
        pnode = parent(tree, node)
        haskey(tree.TREECORE.children, pnode) && filter!(!=(node), tree.TREECORE.children[pnode])
    end
    pop!(tree.TREECORE.contents, node)
    pop!(tree.TREECORE.parent, node, nothing)
    pop!(tree.TREECORE.children, node)
    return tree
end

"""
    getindex(tree::AbstractSimpleTree{N}, node::N) where N -> N

Get the data of a tree's node.
"""
Base.getindex(tree::AbstractSimpleTree{N}, node::N) where {N} = tree.TREECORE.contents[node]

"""
    setindex!(tree::AbstractSimpleTree{N, D}, data::D, node::N) where {N, D}

Set the data of a tree's node.
"""
Base.setindex!(tree::AbstractSimpleTree{N, D}, data::D, node::N) where {N, D} = (tree.TREECORE.contents[node] = data)

"""
    empty(tree::AbstractSimpleTree)

Construct an empty tree of the same type with the input one.
"""
@generated function Base.empty(tree::AbstractSimpleTree{N, D}) where {N, D}
    contents = Expr[:(getfield(tree, $name)) for name in QuoteNode.(fieldnames(tree)) if name != QuoteNode(:TREECORE)]
    return :(($tree)($(contents...), SimpleTreeCore{$N, $D}()))
end

"""
    keytype(tree::AbstractSimpleTree)
    keytype(::Type{<:AbstractSimpleTree{N}}) where N

Get a tree's node type.
"""
Base.keytype(tree::AbstractSimpleTree) = keytype(typeof(tree))
Base.keytype(::Type{<:AbstractSimpleTree{N}}) where N = N

"""
    valtype(tree::AbstractSimpleTree)
    valtype(::Type{<:AbstractSimpleTree{N, D} where N}) where D=D

Get a tree's data type.
"""
Base.valtype(tree::AbstractSimpleTree) = valtype(typeof(tree))
Base.valtype(::Type{<:AbstractSimpleTree{N, D} where N}) where D = D

"""
    ==(t1::T, t2::T) where T<:AbstractSimpleTree -> Bool

Overloaded equivalent operator.
"""
Base.:(==)(t1::T, t2::T) where T <: AbstractSimpleTree = ==(efficientoperations, t1, t2)

"""
    isequal(t1::T, t2::T) where T<:AbstractSimpleTree -> Bool

Overloaded equivalent operator.
"""
Base.isequal(t1::T, t2::T) where T <: AbstractSimpleTree = isequal(efficientoperations, t1, t2)

"""
    keys(tree::AbstractSimpleTree{N}, ::SimpleTreeDepth, node::Union{N, Nothing}=root(tree)) where N
    keys(tree::AbstractSimpleTree{N}, ::SimpleTreeWidth, node::Union{N, Nothing}=root(tree)) where N

Iterate over a tree's nodes starting from a certain `node` by depth first search or width first search.
"""
function Base.keys(tree::AbstractSimpleTree{N}, ti::SimpleTreeIteration, node::Union{N, Nothing}=root(tree)) where N
    TreeKeys{typeof(ti), N, valtype(tree), typeof(tree)}(tree, node)
end
struct TreeKeys{P<:SimpleTreeIteration, N, D, T<:AbstractSimpleTree{N, D}}
    tree::T
    node::Union{N, Nothing}
end
Base.eltype(::Type{<:TreeKeys{<:SimpleTreeIteration, N, D, <:AbstractSimpleTree{N, D}} where D}) where {N} = N
Base.IteratorSize(::Type{<:TreeKeys}) = Base.SizeUnknown()
Base.iterate(tk::TreeKeys) = (root(tk.tree) === nothing) ? nothing : (tk.node, copy(children(tk.tree, tk.node)))
function Base.iterate(tk::TreeKeys{SimpleTreeDepth, N}, state::Vector{N}) where N
    (length(state) == 0) ? nothing : (node = popfirst!(state); prepend!(state, children(tk.tree, node)); (node, state))
end
function Base.iterate(tk::TreeKeys{SimpleTreeWidth, N}, state::Vector{N}) where N
    (length(state) == 0) ? nothing : (node = popfirst!(state); append!(state, children(tk.tree, node)); (node, state))
end

"""
    values(tree::AbstractSimpleTree{N}, ::SimpleTreeDepth, node::Union{N, Nothing}=root(tree)) where N
    values(tree::AbstractSimpleTree{N}, ::SimpleTreeWidth, node::Union{N, Nothing}=root(tree)) where N

Iterate over a tree's data starting from a certain `node` by depth first search or width first search.
"""
Base.values(tree::AbstractSimpleTree{N}, ti::SimpleTreeIteration, node::Union{N, Nothing}=root(tree)) where {N} = (tree[key] for key in keys(tree, ti, node))

"""
    pairs(tree::AbstractSimpleTree{N}, ::SimpleTreeDepth, node::Union{N, Nothing}=root(tree)) where N
    pairs(tree::AbstractSimpleTree{N}, ::SimpleTreeWidth, node::Union{N, Nothing}=root(tree)) where N

Iterate over a tree's (node, data) pairs starting from a certain `node` by depth first search or width first search.
"""
Base.pairs(tree::AbstractSimpleTree{N}, ti::SimpleTreeIteration, node::Union{N, Nothing}=root(tree)) where {N} = ((key, tree[key]) for key in keys(tree, ti, node))

"""
    isleaf(tree::AbstractSimpleTree{N}, node::N) where N -> Bool

Judge whether a tree's node is a leaf (a node without children) or not.
"""
isleaf(tree::AbstractSimpleTree{N}, node::N) where {N} = length(children(tree, node)) == 0

"""
    level(tree::AbstractSimpleTree{N}, node::N) where N -> Int

Get the level of tree's node.
"""
function level(tree::AbstractSimpleTree{N}, node::N) where N
    result = 1
    while node != root(tree)
        result += 1
        node = parent(tree, node)
    end
    return result
end

"""
    ancestor(tree::AbstractSimpleTree{N}, node::N, generation::Int=1) where N -> N

Get the ancestor of a tree's node of the n-th generation.
"""
function ancestor(tree::AbstractSimpleTree{N}, node::N, generation::Int=1) where N
    @assert generation >= 0 "ancestor error: generation($generation) must be non-negative."
    result = node
    for i = 1:generation
        result = parent(tree, result)
    end
    result
end

"""
    descendants(tree::AbstractSimpleTree{N}, node::N, generation::Int=1) where N -> Vector{N}

Get the descendants of a tree's node of the nth generation.
"""
function descendants(tree::AbstractSimpleTree{N}, node::N, generation::Int=1) where N
    @assert generation >= 0 "descendants error: generation($generation) must be non-negative."
    result = N[node]
    for i = 1:generation
        result = vcat((children(tree, node) for node in result)...)
    end
    result
end

"""
    siblings(tree::AbstractSimpleTree{N}, node::N) where N -> Vector{N}

Get the siblings (other nodes sharing the same parent) of a tree's node.
"""
siblings(tree::AbstractSimpleTree{N}, node::N) where N = filter(!=(node), children(tree, parent(tree, node)))

"""
    leaves(tree::AbstractSimpleTree) -> Vector{keytype(tree)}

Get a tree's leaves.
"""
leaves(tree::AbstractSimpleTree) = keytype(tree)[node for node in keys(tree, simpletreedepth) if isleaf(tree, node)]

"""
    push!(tree::AbstractSimpleTree{N, D}, node::N, data::D) where {N, D} -> typeof(tree)
    push!(tree::AbstractSimpleTree{N, D}, parent::Union{N, Nothing}, node::N, data::D) where {N, D} -> typeof(tree)

Push a new node to a tree. When `parent` is `nothing`, this function set the root node of an empty tree.
"""
Base.push!(tree::AbstractSimpleTree{N, D}, node::N, data::D) where {N, D} = push!(tree, nothing, node, data)
function Base.push!(tree::AbstractSimpleTree{N, D}, parent::Union{N, Nothing}, node::N, data::D) where {N, D}
    addnode!(tree, parent, node)
    tree[node] = data
    return tree
end

"""
    append!(tree::AbstractSimpleTree{N, D}, subtree::AbstractSimpleTree{N, D}) where {N, D} -> typeof(tree)
    append!(tree::AbstractSimpleTree{N, D}, node::Union{N, Nothing}, subtree::AbstractSimpleTree{N, D}) where {N, D} -> typeof(tree)

Append a subtree to a tree.
"""
Base.append!(tree::AbstractSimpleTree{N, D}, subtree::AbstractSimpleTree{N, D}) where {N, D} = append!(tree, nothing, subtree)
function Base.append!(tree::AbstractSimpleTree{N, D}, node::Union{N, Nothing}, subtree::AbstractSimpleTree{N, D}) where {N, D}
    for (key, value) in pairs(subtree, simpletreewidth)
        push!(tree, parent(subtree, key, node), key, value)
    end
    return tree
end

"""
    delete!(tree::AbstractSimpleTree{N}, node::N) where N -> typeof(tree)

Delete a node and all its descendants from a tree.
"""
function Base.delete!(tree::AbstractSimpleTree{N}, node::N) where N
    for key in collect(N, keys(tree, simpletreedepth, node))
        deletenode!(tree, key)
    end
    return tree
end

"""
    empty!(tree::AbstractSimpleTree) -> typeof(tree)

Empty a tree.
"""
Base.empty!(tree::AbstractSimpleTree) = delete!(tree, root(tree))

"""
    subtree(tree::AbstractSimpleTree{N}, node::N) where N -> typeof(tree)

Get a subtree whose root is `node`.
"""
function subtree(tree::AbstractSimpleTree{N}, node::N) where N
    result = empty(tree)
    for (i, (key, value)) in enumerate(pairs(tree, simpletreedepth, node))
        push!(result, (i == 1) ? nothing : parent(tree, key), key, value)
    end
    return result
end

"""
    move!(tree::AbstractSimpleTree{N}, node::N, parent::N) where N -> typeof(tree)

Move a subtree to a new position.
"""
function move!(tree::AbstractSimpleTree{N}, node::N, parent::N) where N
    sub = subtree(tree, node)
    delete!(tree, node)
    append!(tree, parent, sub)
    return tree
end

"""
    SimpleTreeCore()

The core of a tree.
"""
mutable struct SimpleTreeCore{N, D}
    root::Union{N, Nothing}
    contents::Dict{N, D}
    parent::Dict{N, N}
    children::Dict{N, Vector{N}}
    SimpleTreeCore{N, D}() where {N, D} = new{N, D}(nothing, Dict{N, D}(), Dict{N, N}(), Dict{N, Vector{N}}())
end

"""
    ==(tc1::TC, tc2::TC) where TC<:SimpleTreeCore -> Bool
    isequal(tc1::TC, tc2::TC) where TC<:SimpleTreeCore -> Bool

Overloaded equivalent operator.
"""
Base.:(==)(tc1::TC, tc2::TC) where {TC<:SimpleTreeCore} = ==(efficientoperations, tc1, tc2)
Base.isequal(tc1::TC, tc2::TC) where {TC<:SimpleTreeCore} = isequal(efficientoperations, tc1, tc2)

"""
    @simpletree structdef treeparams::Union{Expr, Nothing}=nothing

Decorate a "raw" struct to be a subtype of `AbstractSimpleTree`.
!!! note
    1. A "raw" struct means:
       - It has no explicit supertype;
       - It has no inner constructor;
       - It has no attribute `:TREECORE`.
    2. The keytype and valtype can be assigned by the argument `treeparams` in the form `{keytype, valtype}`.
       - When the formal argument names of keytype and valtype are not assigned, they can be automatically generated by the functioin `gensym`.
         For example, the struct `Tree` after the decoration by the following code
         ```julia
         @simpletree(struct Tree{N} info::Vector{N} end)
         ```
         will have three type parameters.
       - When the keytype or valtype is a definite type, use symbol `:(::)` before that type.
         For example, the struct `Tree` decorated by the following code
         ```julia
         @simpletree(struct Tree{N} info::Vector{N} end, {::String, D<:Number})
         ```
         only has two type parameters `N` and `D<:Number`.
       - When the formal argument names of keytype and valtype overlap with those of the raw struct type parameters, the duplicates will be considered as the same.
         For example, the decorated struct `Tree` by the following code
         ```julia
         @simpletree(struct Tree{D} info::Vector{D} end, {::String, D<:Number})
         ```
         has only one type parameter `D<:Number`, where the `D` in the `info::Vector{D}` is the same `D` with that in the decorated attribute `TREECORE::SimpleTreeCore{String, D}`.
       - When the formal argument names of keytype and valtype have no intersection with those of the raw struct type parameters,
         the type parameters of the decorated struct will be just extended by keytype and valtype.
         For example, the decorated struct `Tree` by the following code
         ```julia
         @simpletree(struct Tree{T} info::Vector{T} end, {N<:Tuple, D<:Number})
         ```
         has 3 type parameters, `T`, `N<:Tuple` and `D<:Number`.
"""
macro simpletree(structdef, treeparams::Union{Expr, Nothing}=nothing)
    tf = TypeFactory(structdef)
    fieldnames = [field.name for field in tf.fields]
    paramnames = [param.name for param in tf.params]
    @assert tf.supertype == Inference(:Any) "@simpletree error: no explicit supertype except `Any` is allowed."
    @assert length(tf.constructors) == 0 "@simpletree error: no inner constructor is allowed."
    @assert :TREECORE ∉ fieldnames "@simpletree error: :TREECORE is a reserved attribute name."
    if treeparams === nothing
        treeparamnames = (gensym(), gensym())
        treeparams = [Parameter(treeparamnames[1]), Parameter(treeparamnames[2])]
        truetreeparamnames = treeparamnames
        escapednames = ()
    else
        @assert (treeparams.head == :braces) && (length(treeparams.args) == 2) "@simpletree error: not supported treeparams."
        treeparams = [Parameter(arg) for arg in treeparams.args]
        (treeparams[1].name === nothing) && (treeparams[1].head ≠ :(::)) && (treeparams[1].name = gensym())
        (treeparams[2].name === nothing) && (treeparams[2].head ≠ :(::)) && (treeparams[2].name = gensym())
        treeparamnames = tuple(((tp.name === nothing) ? tp(rawexpr) : tp.name for tp in treeparams)...)
        truetreeparamnames = tuple((tp.name for tp in treeparams if tp.name ≠ nothing)...)
        escapednames = vcat((names(tp(rawexpr)) for tp in treeparams if tp.name === nothing)...)
        filter!(tp->(tp.name ∉ paramnames) && (tp.name ≠ nothing), treeparams)
    end
    tf.supertype = Inference(:(AbstractSimpleTree{$(treeparamnames...)}))
    append!(paramnames, (tp.name for tp in treeparams))
    addfields!(tf, :(TREECORE::SimpleTreeCore{$(treeparamnames...)}))
    addparams!(tf, treeparams...)
    outer = FunctionFactory(name=tf.name)
    addparams!(outer, truetreeparamnames...)
    addargs!(outer, (Argument(name=field.name, type=field.type) for field in tf.fields[1:end-1])...)
    addwhereparams!(outer, tf.params...)
    extendbody!(outer, :($(tf.name)($(fieldnames...), SimpleTreeCore{$(treeparamnames...)}())))
    sbtreedef = tf(MixEscaped(Escaped(tf.name, escapednames...), UnEscaped(paramnames..., :AbstractSimpleTree, :SimpleTreeCore)))
    outerdef = outer(MixEscaped(Escaped(tf.name, escapednames...), UnEscaped(paramnames...)))
    return Expr(:block, :(Base.@__doc__($sbtreedef)), outerdef)
end

"""
    SimpleTree{N, D}() where {N, D}

The minimum tree structure that implements all the default tree methods.
"""
struct SimpleTree{N, D} <: AbstractSimpleTree{N, D}
    TREECORE::SimpleTreeCore{N, D}
end
SimpleTree{N, D}() where {N, D} = SimpleTree(SimpleTreeCore{N, D}())

end #module


text:
# conversions between ROS and _PCL

@info "Caesar._PCL is loading tools related to RobotOS.jl."
@debug "Caesar._PCL RobotOS.jl requires user to first run as @rostypegen (on messages below) in `Main` before loading Caesar."
@debug "  @rosimport sensor_msgs.msg: PointCloud2"

using ..RobotOS

@rosimport sensor_msgs.msg: PointCloud2

Base.convert(::Type{UInt64}, rost::RobotOS.Time) = UInt64(rost.secs)*1000000 + trunc(UInt64, rost.nsecs*1e-3)

# Really odd constructor, strong assumption that user FIRST ran @rostypegen in Main BEFORE loading Caesar
# https://docs.ros.org/en/hydro/api/pcl_conversions/html/pcl__conversions_8h_source.html#l00208
function PCLPointCloud2(msg::Main.sensor_msgs.msg.PointCloud2)
  header = Header(;stamp  = msg.header.stamp,
                  seq     = msg.header.seq,
                  frame_id= msg.header.frame_id )
  #

  # all PointField elements
  pfs = PointField[PointField(;name=pf_.name,offset=pf_.offset,datatype=pf_.datatype,count=pf_.count) for pf_ in msg.fields]

  endian = msg.is_bigendian ? _PCL_ENDIAN_BIG_BYTE : _PCL_ENDIAN_LITTLE_WORD
  pc2 = PCLPointCloud2(;header,
                        height     = msg.height,
                        width      = msg.width,
                        fields     = pfs,
                        data       = msg.data,
                        is_bigendian= endian,
                        point_step = msg.point_step,
                        row_step   = msg.row_step,
                        is_dense   = UInt8(msg.is_dense) )
  #
end

#


text:
abstract Covariance
import Distances

abstract Kernel

function flatten_hyp(m)
	lik_hyp = Float64[]
	for n in names(m.lik)
		push!(lik_hyp, m.lik.(n)...)
	end
	cov_hyp = Float64[]
	for n in names(m.k)
		push!(cov_hyp, m.k.(n)...)
	end
	vcat(lik_hyp, cov_hyp)

end

function inflate_hyp(m, hyp_vec)
	liks = hyp_vec[1:length(m.lik)]
	covs = hyp_vec[length(m.lik) + 1:end]
	typeof(m.k)(covs), typeof(m.lik)(liks)
end

type SqExpIso <: Kernel
	log_l
	log_s
end
SqExpIso(v::Array) = SqExpIso(v[1], v[2])

function kernel{T <: FloatingPoint}(h::SqExpIso, X1::Array{T, 2},
						X2::Array{T, 2})
	l2 = exp(2h.log_l)
	s2 = exp(2h.log_s)
	s2 * exp(-Distances.pairwise(Distances.SqEuclidean(), X1', X2') / (2l2))
end
kernel{T <: FloatingPoint}(h::SqExpIso, X::Array{T, 2}) = kernel(h, X, X)

type SqExpIsoPeriodic <: Kernel
	log_l
	log_s
	log_period
end
SqExpIsoPeriodic(v::Array) = SqExpIsoPeriodic(v[1], v[2], v[3])

function kernel{T <: FloatingPoint}(h::SqExpIsoPeriodic, X1::Array{T, 2},
						X2::Array{T, 2})
	l2 = exp(2h.log_l)
	s2 = exp(2h.log_s)
	period1 = exp(h.log_period)
	s2 *
	exp(-Distances.pairwise(Distances.SqEuclidean(),
		cos(2pi * X1 / period1)', cos(2pi * X2 / period1)') / (2l2))
end
kernel{T <: FloatingPoint}(h::SqExpIsoPeriodic, X::Array{T, 2}) = kernel(h, X, X)

type SqExpIso2Periodic <: Kernel
	log_l
	log_s
	log_period1
	log_period2
end
SqExpIso2Periodic(v::Array) = SqExpIso2Periodic(v[1], v[2], v[3], v[4])

function kernel{T <: FloatingPoint}(h::SqExpIso2Periodic, X1::Array{T, 2},
						X2::Array{T, 2})
	l2 = exp(2h.log_l)
	s2 = exp(2h.log_s)
	period1 = exp(h.log_period1)
	period2 = exp(h.log_period2)
	s2 *
	(exp(-Distances.pairwise(Distances.SqEuclidean(),
		cos(2pi * X1 / period1)', cos(2pi * X2 / period1)') / (2l2)) .+
	exp(-Distances.pairwise(Distances.SqEuclidean(),
		cos(2pi * X1 / period2)', cos(2pi * X2 / period2)') / (2l2)))
end
kernel{T <: FloatingPoint}(h::SqExpIso2Periodic, X::Array{T, 2}) = kernel(h, X, X)

type SqExpIso3Periodic <: Kernel
	log_l
	log_s
	log_period1
	log_period2
	log_period3
end
SqExpIso3Periodic(v::Array) = SqExpIso(v[1], v[2], v[3], v[4], v[5])

function kernel{T <: FloatingPoint}(h::SqExpIso3Periodic, X1::Array{T, 2},
						X2::Array{T, 2})
	l2 = exp(2h.log_l)
	s2 = exp(2h.log_s)
	period1 = exp(h.log_period1)
	period2 = exp(h.log_period2)
	period3 = exp(h.log_period3)
	s2 *
	(exp(-Distances.pairwise(Distances.SqEuclidean(),
		cos(2pi * X1 / period1)', cos(2pi * X2 / period1)') / (2l2)) .+
	exp(-Distances.pairwise(Distances.SqEuclidean(),
		cos(2pi * X1 / period2)', cos(2pi * X2 / period2)') / (2l2)) .+
	exp(-Distances.pairwise(Distances.SqEuclidean(),
		cos(2pi * X1 / period3)', cos(2pi * X2 / period3)') / (2l2)))
end
kernel{T <: FloatingPoint}(h::SqExpIso3Periodic, X::Array{T, 2}) = kernel(h, X, X)

type SqExpIso4Periodic <: Kernel
	log_l
	log_s
	log_period1
	log_period2
	log_period3
	log_period4
end
SqExpIso4Periodic(v::Array) = SqExpIso(v[1], v[2], v[3], v[4],
	v[5], v[6], v[7])

function kernel{T <: FloatingPoint}(h::SqExpIso4Periodic, X1::Array{T, 2},
						X2::Array{T, 2})
	l2 = exp(2h.log_l)
	s2 = exp(2h.log_s)
	period1 = exp(h.log_period1)
	period2 = exp(h.log_period2)
	period3 = exp(h.log_period3)
	period4 = exp(h.log_period4)
	s2 *
	(exp(-Distances.pairwise(Distances.SqEuclidean(),
		cos(2pi * X1 / period1)', cos(2pi * X2 / period1)') / (2l2)) .+
	exp(-Distances.pairwise(Distances.SqEuclidean(),
		cos(2pi * X1 / period2)', cos(2pi * X2 / period2)') / (2l2)) .+
	exp(-Distances.pairwise(Distances.SqEuclidean(),
		cos(2pi * X1 / period3)', cos(2pi * X2 / period3)') / (2l2)) .+
	exp(-Distances.pairwise(Distances.SqEuclidean(),
		cos(2pi * X1 / period4)', cos(2pi * X2 / period4)') / (2l2)))
end
kernel{T <: FloatingPoint}(h::SqExpIso4Periodic, X::Array{T, 2}) = kernel(h, X, X)

type MaternIso32 <: Kernel
	log_l
	log_s
end
MaternIso32(v::Array) = MaternIso32(v[1], v[2])
function kernel{T <: FloatingPoint}(h::MaternIso32, X1::Array{T, 2},
						X2::Array{T, 2})
	l = exp(h.log_l)
	s = exp(h.log_s)
	s3dol = sqrt(3) * Distances.pairwise(Distances.Euclidean(), X1', X2') / l
	s * (1 .+ s3dol) .* exp(-s3dol)
end
kernel{T <: FloatingPoint}(h::MaternIso32, X::Array{T, 2}) = kernel(h, X, X)

type MaternARD32 <: Kernel
	log_ls::Array{Float64,}
	log_s
end
MaternARD32(v::Array) = MaternARD32(v[1:end-1], v[end])
function kernel{T <: FloatingPoint}(h::MaternARD32, X1::Array{T, 2},
						X2::Array{T, 2})
	ls2 = 2exp(2h.log_ls[:])
	s = exp(h.log_s)
	s3dol = sqrt(3) * Distances.pairwise(Distances.WeightedEuclidean(1./ls2), X1', X2')
	s * (1 .+ s3dol) .* exp(-s3dol)
end
kernel{T <: FloatingPoint}(h::MaternARD32, X::Array{T, 2}) = kernel(h, X, X)


text:
abstract type AbstractProx <: AbstractPolicy end

# function prox!(prox::AbstractProx,step::Real,xprev::AbstractVector,gcurr::AbstractVector,xcurr::AbstractVector)
#     error("No defined prox! function for prox policy ", typeof(prox))
# end

function prox!() end

function prox_wrapper(step::Cdouble, xprev_b::Ptr{Cdouble},
                      xprev_e::Ptr{Cdouble}, gcurr_b::Ptr{Cdouble}, xcurr_b::Ptr{Cdouble},
                      prox_data::Ptr{Cvoid})::Ptr{Cdouble}
    prox_policy = unsafe_pointer_to_objref(prox_data)::AbstractProx
    ptrdiff = Int(xprev_e - xprev_b)
    N = divrem(ptrdiff, sizeof(Cdouble))[1]
    xprev = unsafe_wrap(Array, xprev_b, N)
    gcurr = unsafe_wrap(Array, gcurr_b, N)
    xcurr = unsafe_wrap(Array, xcurr_b, N)
    prox!(prox_policy, step, xprev, gcurr, xcurr)
    return xcurr_b + ptrdiff
end

module Prox

using LinearAlgebra
using Parameters
using POLO: AbstractProx, AbstractPolicyParameters
import POLO: initialize!, prox!

function initialize!(policy::AbstractProx,x₀::Vector{Float64})
    nothing
end

include("none.jl")
include("affine.jl")
include("box.jl")
include("halfspace.jl")
include("l1norm.jl")
include("l2ball.jl")
include("orthant.jl")

end


text:
"""
    DependencyTree

A natural language sentence annotated with dependency structure.
"""
struct DependencyTree{T<:Dependency} <: AbstractGraph{Int}
    graph::SimpleDiGraph
    tokens::Vector{T}
    mwts::Vector{MultiWordToken}
    emptytokens::Vector{EmptyToken}
    root::Int

    function DependencyTree(graph, tokens, mwts, emptytokens, root; check=true, kwargs...)
        g = new{eltype(tokens)}(graph, tokens, mwts, emptytokens, root)
        check && check_tree(g; kwargs...)
        return g
    end
end

"""
   DependencyTree(T, tokens; add_id=false; kwargs...)

Make a dependency tree with tokens of type `T` from `tokens`.
"""
function DependencyTree(T::Type{<:Dependency}, tokens; add_id=false, kwargs...)
    A, mwts, es = T[], MultiWordToken[], EmptyToken[]
    for (id, token) in enumerate(tokens)
        add_id ? dependency = T(id, token...) : dependency = T(token...)
        push!(A, dependency)
    end
    DependencyTree(A; mwts=mwts, emptytokens=es, kwargs...)
end

function DependencyTree(tokens::Vector{<:Dependency}; mwts=MultiWordToken[], emptytokens=EmptyToken[], kwargs...)
    rt = 0
    graph = SimpleDiGraph(length(tokens))
    for dep in tokens
        isroot(dep) && continue
        i, h = id(dep), head(dep)
        iszero(h) && (rt = i)
        add_edge!(graph, h, i) # arrows point from head to dependent
    end
    return DependencyTree(graph, tokens, mwts, emptytokens, rt; kwargs...)
end

function DependencyTree{T}(lines::AbstractVector{S}; add_id=false, kwargs...) where {T<:Dependency,S<:AbstractString}
    A, mwts, emptytokens = T[], MultiWordToken[], EmptyToken[]
    for (i, line) in enumerate(lines)
        try
            push!(A, T(line, add_id=add_id))
        catch err
            if isa(err, MultiWordTokenError)
                push!(mwts, MultiWordToken(line))
            elseif isa(err, EmptyTokenError)
                push!(emptytokens, EmptyToken(line))
            else
                throw(err)
            end
        end
    end
    DependencyTree(A; mwts=mwts, emptytokens=emptytokens, kwargs...)
end

"""
    DependencyTree{T}(sentence::String; add_id=false, kwargs...)

Read a dependency tree from `sentence`.
"""
function DependencyTree{T}(sentence::String; add_id=false, kwargs...) where T
    ls = String.(filter(x -> x != "", split(strip(sentence), "\n")))
    DependencyTree{T}(ls; add_id=false, kwargs...)
end

"""
    check_tree(g, check_single_head=true, check_has_root=true, check_projective=false)

Ensure the well-formedness of the dependency graph `g`, throwing an
error if g is not well-formed.
"""
function check_tree(tree::DependencyTree; check_single_head=true, check_has_root=true,
                        check_projective=false)
    check_has_root && iszero(tree.root) && throw(RootlessGraphError(tree))
    if check_single_head
        !has1root(tree) && throw(MultipleRootsError(tree))
        !is_weakly_connected(tree.graph) && throw(GraphConnectivityError(tree,"not weakly connected!"))
    end
    check_projective && !isprojective(tree) && throw(NonProjectiveGraphError(tree))
    for i in 1:length(tree)
        n_inc = length(inneighbors(tree.graph, i))
        # root node and its dependency on predicate are
        # represented implicitly, so 0 is expected here
        n_inc == 0 && head(tree, i) == 0 ? continue :
        n_inc != 1 && throw(GraphConnectivityError(tree, "node $i should have exactly 1 incoming connection (has $n_inc)"))
    end
    return nothing
end

has1root(tree::DependencyTree) = count(iszero ∘ head, tokens(tree)) == 1

isprojective(tree::DependencyTree) =
    # For every arc (i,l,j) there is a directed path from i to every
    # word k such that min(i,j) < k < max(i,j)
    all(isprojective(tree, src(e), dst(e)) for e in edges(tree.graph))

isprojective(g::DependencyTree, head::Int, dep::Int) =
    all(k -> has_path(g.graph, head, k), min(head, dep):max(head, dep))

dependents(g::DependencyTree, id::Int) =
    iszero(id) ? [g.root] : outneighbors(g.graph, id)
deprel(g::DependencyTree, id::Int) = deprel(g[id])
deptype(g::DependencyTree) = eltype(g.tokens)
form(g::DependencyTree, id::Int) = form(g[id])
has_arc(g::DependencyTree, h::Int, d::Int) = head(g[d]) == h
has_dependency(g::DependencyTree, h::Int, d::Int) = head(g[d]) == h
head(g::DependencyTree, id::Int) = head(g[id])
root(g::DependencyTree) = g[0]
token(g::DependencyTree, id) = g[id]
tokens(g::DependencyTree) = g.tokens
tokens(g::DependencyTree, is) = [token(g, i) for i in is if 0 <= i <= length(g)]

function leftdeps(g::DependencyTree, dep::Dependency)
    i = id(dep)
    g[filter(d -> d < i, dependents(g, id(dep)))]
end
leftdeps(g::DependencyTree, i::Int) = filter(d -> d < i, dependents(g, i))

function rightdeps(g::DependencyTree, dep::Dependency)
    i = id(dep)
    g[filter(d -> d > i, dependents(g, id(dep)))]
end
rightdeps(g::DependencyTree, i::Int) = filter(d -> d > i, dependents(g, i))

leftmostdep(g::DependencyTree, args...) = leftmostdep(g.tokens, args...)
rightmostdep(g::DependencyTree, args...) = rightmostdep(g.tokens, args...)

toconllu(g::DependencyTree) = string(toconllu.(g)...)

function Base.show(io::IO, g::DependencyTree)
    T = "DependencyTree{$(eltype(g))}\n"
    token(t) = "$(id(t))\t$(form(t))\t$(deprel(t))\t$(head(t))"
    print(io,T,join([token(t) for t in g],"\n"))
end
function Base.show(io::IO, g::DependencyTree{UntypedDependency})
    T = "DependencyTree{$(eltype(g))}\n"
    token(t) = "$(id(t))\t$(form(t))\t$(head(t))"
    print(io,T,join([token(t) for t in g],"\n"))
end

==(g1::DependencyTree, g2::DependencyTree) = all(g1.tokens .== g2.tokens)
Base.eltype(g::DependencyTree) = eltype(g.tokens)
Base.getindex(g::DependencyTree, i) = i == 0 ? root(eltype(g)) : g.tokens[i]
Base.iterate(g::DependencyTree, state=1) = iterate(g.tokens, state)
Base.length(g::DependencyTree) = length(g.tokens)


text:
import Base.==, Base.<=, Base.>=, Base.<, Base.>

### Linear equality constraint
mutable struct EqConstraint <: Constraint
    head::Symbol
    id_hash::UInt64
    lhs::AbstractExpr
    rhs::AbstractExpr
    size::Tuple{Int,Int}
    dual::ValueOrNothing

    function EqConstraint(lhs::AbstractExpr, rhs::AbstractExpr)
        if lhs.size == rhs.size || lhs.size == (1, 1)
            sz = rhs.size
        elseif rhs.size == (1, 1)
            sz = lhs.size
        else
            error(
                "Cannot create equality constraint between expressions of size $(lhs.size) and $(rhs.size)",
            )
        end
        id_hash = hash((lhs, rhs, :(==)))
        return new(:(==), id_hash, lhs, rhs, sz, nothing)
    end
end

function vexity(c::EqConstraint)
    vex = vexity(c.lhs) + (-vexity(c.rhs))
    # You can't have equality constraints with concave/convex expressions
    if vex == ConvexVexity() || vex == ConcaveVexity()
        vex = NotDcp()
    end
    return vex
end

function conic_form!(c::EqConstraint, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, c)
        if !(sign(c.lhs) == ComplexSign() || sign(c.rhs) == ComplexSign())
            expr = c.lhs - c.rhs
            objective = conic_form!(expr, unique_conic_forms)
            new_constraint =
                ConicConstr([objective], :Zero, [c.size[1] * c.size[2]])
            unique_conic_forms.conic_constr_to_constr[new_constraint] = c
        else
            real_expr = real(c.lhs - c.rhs)
            imag_expr = imag(c.lhs - c.rhs)
            real_objective = conic_form!(real_expr, unique_conic_forms)
            imag_objective = conic_form!(imag_expr, unique_conic_forms)
            new_constraint = ConicConstr(
                [real_objective, imag_objective],
                :Zero,
                [c.size[1] * c.size[2], c.size[1] * c.size[2]],
            )
            unique_conic_forms.conic_constr_to_constr[new_constraint] = c
        end
        cache_conic_form!(unique_conic_forms, c, new_constraint)
    end
    return get_conic_form(unique_conic_forms, c)
end

==(lhs::AbstractExpr, rhs::AbstractExpr) = EqConstraint(lhs, rhs)
==(lhs::AbstractExpr, rhs::Value) = ==(lhs, Constant(rhs))
==(lhs::Value, rhs::AbstractExpr) = ==(Constant(lhs), rhs)

### Linear inequality constraints
mutable struct LtConstraint <: Constraint
    head::Symbol
    id_hash::UInt64
    lhs::AbstractExpr
    rhs::AbstractExpr
    size::Tuple{Int,Int}
    dual::ValueOrNothing

    function LtConstraint(lhs::AbstractExpr, rhs::AbstractExpr)
        if sign(lhs) == ComplexSign() || sign(rhs) == ComplexSign()
            error(
                "Cannot create inequality constraint between expressions of sign $(sign(lhs)) and $(sign(rhs))",
            )
        else
            if lhs.size == rhs.size || lhs.size == (1, 1)
                sz = rhs.size
            elseif rhs.size == (1, 1)
                sz = lhs.size
            else
                error(
                    "Cannot create inequality constraint between expressions of size $(lhs.size) and $(rhs.size)",
                )
            end
        end
        id_hash = hash((lhs, rhs, :(<=)))
        return new(:(<=), id_hash, lhs, rhs, sz, nothing)
    end
end

function vexity(c::LtConstraint)
    vex = vexity(c.lhs) + (-vexity(c.rhs))
    if vex == ConcaveVexity()
        vex = NotDcp()
    end
    return vex
end

function conic_form!(c::LtConstraint, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, c)
        expr = c.rhs - c.lhs
        objective = conic_form!(expr, unique_conic_forms)
        new_constraint =
            ConicConstr([objective], :NonNeg, [c.size[1] * c.size[2]])
        unique_conic_forms.conic_constr_to_constr[new_constraint] = c
        cache_conic_form!(unique_conic_forms, c, new_constraint)
    end
    return get_conic_form(unique_conic_forms, c)
end

<=(lhs::AbstractExpr, rhs::AbstractExpr) = LtConstraint(lhs, rhs)
<=(lhs::AbstractExpr, rhs::Value) = <=(lhs, Constant(rhs))
<=(lhs::Value, rhs::AbstractExpr) = <=(Constant(lhs), rhs)
<(lhs::AbstractExpr, rhs::AbstractExpr) = LtConstraint(lhs, rhs)
<(lhs::AbstractExpr, rhs::Value) = <=(lhs, Constant(rhs))
<(lhs::Value, rhs::AbstractExpr) = <=(Constant(lhs), rhs)

mutable struct GtConstraint <: Constraint
    head::Symbol
    id_hash::UInt64
    lhs::AbstractExpr
    rhs::AbstractExpr
    size::Tuple{Int,Int}
    dual::ValueOrNothing

    function GtConstraint(lhs::AbstractExpr, rhs::AbstractExpr)
        if sign(lhs) == ComplexSign() || sign(rhs) == ComplexSign()
            error(
                "Cannot create inequality constraint between expressions of sign $(sign(lhs)) and $(sign(rhs))",
            )
        else
            if lhs.size == rhs.size || lhs.size == (1, 1)
                sz = rhs.size
            elseif rhs.size == (1, 1)
                sz = lhs.size
            else
                error(
                    "Cannot create inequality constraint between expressions of size $(lhs.size) and $(rhs.size)",
                )
            end
        end
        id_hash = hash((lhs, rhs, :(>=)))
        return new(:(>=), id_hash, lhs, rhs, sz, nothing)
    end
end

function vexity(c::GtConstraint)
    vex = -vexity(c.lhs) + (vexity(c.rhs))
    if vex == ConcaveVexity()
        vex = NotDcp()
    end
    return vex
end

function conic_form!(c::GtConstraint, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, c)
        expr = c.lhs - c.rhs
        objective = conic_form!(expr, unique_conic_forms)
        new_constraint =
            ConicConstr([objective], :NonNeg, [c.size[1] * c.size[2]])
        unique_conic_forms.conic_constr_to_constr[new_constraint] = c
        cache_conic_form!(unique_conic_forms, c, new_constraint)
    end
    return get_conic_form(unique_conic_forms, c)
end

>=(lhs::AbstractExpr, rhs::AbstractExpr) = GtConstraint(lhs, rhs)
>=(lhs::AbstractExpr, rhs::Value) = >=(lhs, Constant(rhs))
>=(lhs::Value, rhs::AbstractExpr) = >=(Constant(lhs), rhs)
>(lhs::AbstractExpr, rhs::AbstractExpr) = GtConstraint(lhs, rhs)
>(lhs::AbstractExpr, rhs::Value) = >=(lhs, Constant(rhs))
>(lhs::Value, rhs::AbstractExpr) = >=(Constant(lhs), rhs)

function +(
    constraints_one::Array{<:Constraint},
    constraints_two::Array{<:Constraint},
)
    constraints = append!(Constraint[], constraints_one)
    return append!(constraints, constraints_two)
end
function +(constraint_one::Constraint, constraint_two::Constraint)
    return [constraint_one] + [constraint_two]
end
function +(constraint_one::Constraint, constraints_two::Array{<:Constraint})
    return [constraint_one] + constraints_two
end
function +(constraints_one::Array{<:Constraint}, constraint_two::Constraint)
    return constraints_one + [constraint_two]
end


text:
"""

    X = L96(Float32)

Integrates the Lorenz1996 model. Standard parameters

    L96(::Type{T}=Float64,              # number format for RHS
    ::Type{Tprog}=T;                    # number format for prognostic variables
    N::Int=10_000,                      # number of time steps
    n::Int=36,                          # number of variables
    X::Array{Float64,1}=zeros(36),      # initial conditions
    α::Real=1.0,                        # friction parameter
    F::Float64=8.0,                     # forcing constant
    s::Float64=1.0,                     # scaling
    η::Float64=0.01,                    # strength of initial perturbation at X[1] if no X provided
    Δt::Float64=0.01,                   # time step
    scheme::String="RK4"                # time integration scheme
    output::Bool=false                  # return every time step?

# Examples
```jldoc
julia> X1 = L96(Float32);
julia> X2 = L96(Float32,n=6,Δt=0.005);
julia> X3 = L96(Float16,Float32);
```
"""
function L96(::Type{T},                         # number format for RHS
            ::Type{Tprog};                      # number format for prognostic variables
            N::Int=10_000,                      # number of time steps
            n::Int=36,                          # number of variables
            X::Vector=zeros(36),                # initial conditions
            α::Real=1.0,                        # friction parameter
            F::Real=8.0,                        # forcing constant
            s::Real=1.0,                        # scaling
            η::Real=0.01,                       # strength of initial perturbation at X[1] if no X provided
            Δt::Real=0.1,                       # time step
            scheme::String="RK4",               # time integration scheme
            output::Bool=false                  # return every time step?
            ) where {T<:AbstractFloat,Tprog<:AbstractFloat}

            # if both n and X are specified then they should match!
            if n != 36 && X != zeros(36) && length(X) != n
                throw(error("n=$n and length(X)=$(length(X)) do not match."))
            end

            # add the forcing for equilibrium initial conditions
            if X == zeros(36)   # the default
                X = zeros(n) .+ F
                X[1] += η # add small perturbation on the first variable
            end

            if scheme == "RK4"
                return RKn(T,Tprog,N,X,α,F,s,Δt,output,RKo=4)
            elseif scheme == "RK3"
                return RKn(T,Tprog,N,X,α,F,s,Δt,output,RKo=3)
            elseif scheme == "RK2"
                return RKn(T,Tprog,N,X,α,F,s,Δt,output,RKo=2)
            else
                throw(error("Scheme $scheme not implemented."))
            end
end

"""Calls L96 if no type is specified. Float64 for T and Tprog."""
function L96(;kwargs...)                 # if no type specified
    L96(Float64,Float64;kwargs...)       # use Float64 for everything
end

"""Calls L96 if only one type T is specified. T for RHS and prognostic variables."""
function L96(::Type{T};kwargs...) where {T<:AbstractFloat}
    L96(T,T;kwargs...)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Telescope"
version = v"0.1.2"

# Collection of sources required to build this package
sources = [
    ArchiveSource("https://github.com/jhigginbotham64/Telescope/archive/refs/tags/v$(version).tar.gz",
              "025a24a5dfe4c4d5910ae0cad7664ec2c2fc21cce8ffaae56cc6b42300249de6"),
]

# Bash recipe for building across all platforms
script = raw"""
export CFLAGS="-I${includedir}"
export CXXFLAGS="-I${includedir}"
cd $WORKSPACE/srcdir/Telescope*
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -S .. \
    -B .
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; exclude=p->!(Sys.islinux(p) && nbits(p) == 64))

platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    FileProduct("include/telescope.h", :telescope_h),
    LibraryProduct("libtelescope", :libtelescope),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("OpenGLMathematics_jll"),
    Dependency("SDL2_jll"),
    Dependency("SDL2_image_jll"),
    Dependency("SDL2_mixer_jll"),
    Dependency("SDL2_net_jll"),
    Dependency("SDL2_ttf_jll"),
    Dependency("Shaderc_jll"),
    BuildDependency("Vulkan_Headers_jll"),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"8.1.0")

max_stars_repo_name: cmhamel/Yggdrasil

max_issues_repo_name: cmhamel/Yggdrasil

max_forks_repo_name: cmhamel/Yggdrasil


text:
include(pwd() * "/test/equation_test.jl")
include(pwd() * "/test/symbolic_utils_tests.jl")
include(pwd() * "/test/preliminary_tests.jl")
include(pwd() * "/test/to_expr_test.jl")
include(pwd() * "/test/domain_test.jl")