text:
################################################################################
#                               VARIABLE DEFINITION
################################################################################
"""
    InfOptVariableType

An abstract `DataType` for variable type objects used to create `InfiniteOpt`
variables via `JuMP.@variable`.
"""
abstract type InfOptVariableType end

# Fallback
function _check_and_make_variable_ref(model::InfiniteModel, v::T, name) where {T}
    throw(ArgumentError("Invalid variable object type `$T`."))
end

# helper function for setting the info constraints (#TODO add test)
function _set_info_constraints(
    info::JuMP.VariableInfo,
    gvref::GeneralVariableRef,
    dvref::DispatchVariableRef
    )::Nothing
    model = JuMP.owner_model(gvref)
    if info.has_lb
        newset = MOI.GreaterThan(info.lower_bound)
        cref = JuMP.add_constraint(model, JuMP.ScalarConstraint(gvref, newset),
                                   is_info_constr = true)
        _set_lower_bound_index(dvref, JuMP.index(cref))
    end
    if info.has_ub
        newset = MOI.LessThan(info.upper_bound)
        cref = JuMP.add_constraint(model, JuMP.ScalarConstraint(gvref, newset),
                                   is_info_constr = true)
        _set_upper_bound_index(dvref, JuMP.index(cref))
    end
    if info.has_fix
        newset = MOI.EqualTo(info.fixed_value)
        cref = JuMP.add_constraint(model, JuMP.ScalarConstraint(gvref, newset),
                                   is_info_constr = true)
        _set_fix_index(dvref, JuMP.index(cref))
    end
    if info.binary
        cref = JuMP.add_constraint(model, JuMP.ScalarConstraint(gvref, MOI.ZeroOne()),
                                   is_info_constr = true)
        _set_binary_index(dvref, JuMP.index(cref))
    elseif info.integer
        cref = JuMP.add_constraint(model, JuMP.ScalarConstraint(gvref, MOI.Integer()),
                                   is_info_constr = true)
        _set_integer_index(dvref, JuMP.index(cref))
    end
    return
end

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

Extend the `JuMP.add_variable` function to accomodate `InfiniteOpt` variable
types. Adds a variable to an infinite model `model` and returns a
[`GeneralVariableRef`](@ref). Primarily intended to be an internal function of
the constructor macro `@variable`. However, it can be used in combination with
`JuMP.build_variable` to add variables to an infinite model object.
Errors if invalid parameters reference(s) or an invalid infinite variable
reference is included in `var`.

**Examples**
```julia-repl
julia> @infinite_parameter(m, t in [0, 10]);

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

julia> inf_var = build_variable(error, info, Infinite(t));

julia> ivref = add_variable(m, inf_var, "var_name")
var_name(t)

julia> pt_var = build_variable(error, info, Point(ivref, 0.5));

julia> pvref = add_variable(m, pt_var, "var_alias")
var_alias

julia> f_var = build_variable(error, info);

julia> fvref = add_variable(m, f_var, "var_name")
var_name
```
"""
function JuMP.add_variable(
    model::InfiniteModel,
    var::JuMP.AbstractVariable,
    name::String = "";
    kwargs...
    )::GeneralVariableRef
    dvref = _check_and_make_variable_ref(model, var, name; kwargs...)
    vindex = JuMP.index(dvref)
    gvref = _make_variable_ref(model, vindex)
    _set_info_constraints(var.info, gvref, dvref)
    model.name_to_var = nothing
    return gvref
end

# Add VariableConstrainedOnCreation (TODO change paradigm to leverage MOI.add_constrained_variable)
function JuMP.add_variable(
    model::InfiniteModel,
    variable::JuMP.VariableConstrainedOnCreation,
    name::String
    )::GeneralVariableRef
    vref = JuMP.add_variable(model, variable.scalar_variable, name)
    con = JuMP.ScalarConstraint(vref, variable.set)
    JuMP.add_constraint(model, con)
    return vref
end

# Add VariablesConstrainedOnCreation (TODO change paradigm to leverage MOI.add_constrained_variables)
function JuMP.add_variable(
    model::InfiniteModel,
    variables::JuMP.VariablesConstrainedOnCreation,
    names::AbstractArray{<:String}
    )
    shape = variables.shape
    vector_names = JuMP.vectorize(names, shape)
    vrefs = JuMP.add_variable.(model, variables.scalar_variables, vector_names)
    shaped_vrefs = JuMP.reshape_vector(vrefs, shape)
    shaped_set = JuMP.reshape_set(variables.set, shape)
    con = JuMP.build_constraint(error, shaped_vrefs, shaped_set)
    JuMP.add_constraint(model, con)
    return shaped_vrefs
end

################################################################################
#                            VARIABLE DEPENDENCIES
################################################################################
# Extend _measure_dependencies
function _measure_dependencies(vref::DecisionVariableRef)::Vector{MeasureIndex}
    return _data_object(vref).measure_indices
end

# Extend _constraint_dependencies
function _constraint_dependencies(vref::DecisionVariableRef)::Vector{InfOptConstraintIndex}
    return _data_object(vref).constraint_indices
end

"""
    used_by_measure(vref::DecisionVariableRef)::Bool

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

**Example**
```julia-repl
julia> used_by_measure(vref)
true
```
"""
function used_by_measure(vref::DecisionVariableRef)::Bool
    return !isempty(_measure_dependencies(vref))
end

"""
    used_by_constraint(vref::DecisionVariableRef)::Bool

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

**Example**
```julia-repl
julia> used_by_constraint(vref)
false
```
"""
function used_by_constraint(vref::DecisionVariableRef)::Bool
    return !isempty(_constraint_dependencies(vref))
end

"""
    used_by_objective(vref::DecisionVariableRef)::Bool

Return a `Bool` indicating if `vref` is used by the objective.

**Example**
```julia-repl
julia> used_by_objective(vref)
true
```
"""
function used_by_objective(vref::DecisionVariableRef)::Bool
    return _data_object(vref).in_objective
end

"""
    is_used(vref::DecisionVariableRef)::Bool

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

**Example**
```julia-repl
julia> is_used(vref)
true
```
"""
function is_used(vref::DecisionVariableRef)::Bool
    return used_by_measure(vref) || used_by_constraint(vref) || used_by_objective(vref)
end

################################################################################
#                                VARIABLE NAMING
################################################################################
"""
    JuMP.name(vref::DecisionVariableRef)::String

Extend `JuMP.name` to return the names of `InfiniteOpt` variables.

**Example**
```julia-repl
julia> name(vref)
"var_name"
```
"""
function JuMP.name(vref::DecisionVariableRef)::String
    object = get(_data_dictionary(vref), JuMP.index(vref), nothing)
    return object === nothing ? "" : object.name
end

"""
    JuMP.set_name(vref::DecisionVariableRef, name::String)::Nothing

Extend `JuMP.set_name` to set names of decision variables.

**Example**
```julia-repl
julia> set_name(vref, "var_name")

julia> name(vref)
"var_name"
```
"""
function JuMP.set_name(vref::DecisionVariableRef, name::String)::Nothing
    _data_object(vref).name = name
    JuMP.owner_model(vref).name_to_var = nothing
    return
end

# Make a variable reference
function _make_variable_ref(
    model::InfiniteModel,
    index::ObjectIndex
    )::GeneralVariableRef
    return GeneralVariableRef(model, index.value, typeof(index))
end

# Get the name_to_var Dictionary
function _var_name_dict(
    model::InfiniteModel
    )::Union{Nothing, Dict{String, ObjectIndex}}
    return model.name_to_var
end

# Update name_to_var
function _update_var_name_dict(
    model::InfiniteModel,
    var_dict::MOIUC.CleverDict
    )::Nothing
    name_dict = _var_name_dict(model)
    for (index, data_object) in var_dict
        var_name = data_object.name
        if haskey(name_dict, var_name)
            name_dict[var_name] = FiniteVariableIndex(-1) # dumby value
        else
            name_dict[var_name] = index
        end
    end
    model.name_to_var = name_dict
    return
end

"""
    JuMP.variable_by_name(model::InfiniteModel,
                          name::String)::Union{GeneralVariableRef, Nothing}

Extend `JuMP.variable_by_name` for `InfiniteModel` objects. Return the variable
reference assoociated with a variable name. Errors if multiple variables have the
same name. Returns nothing if no such name exists.

**Examples**
```julia-repl
julia> variable_by_name(m, "var_name")
var_name

julia> variable_by_name(m, "fake_name")

```
"""
function JuMP.variable_by_name(
    model::InfiniteModel,
    name::String
    )::Union{GeneralVariableRef, Nothing}
    if _var_name_dict(model) === nothing
        model.name_to_var = Dict{String, ObjectIndex}()
        _update_var_name_dict(model, model.infinite_vars)
        _update_var_name_dict(model, model.semi_infinite_vars)
        _update_var_name_dict(model, model.point_vars)
        _update_var_name_dict(model, model.finite_vars)
    end
    index = get(_var_name_dict(model), name, nothing)
    if index isa Nothing
        return nothing
    elseif index == FiniteVariableIndex(-1)
        error("Multiple variables have the name $name.")
    else
        return _make_variable_ref(model, index)
    end
end

################################################################################
#                              PARAMETER REFERENCES
################################################################################
# Extend parameter_refs for variables (this serves as a fallback for finite types)
function parameter_refs(vref::FiniteRef)::Tuple
    return ()
end

################################################################################
#                         VARIABLE OBJECT MODIFICATION
################################################################################
# Extend _set_core_variable_object
function _set_core_variable_object(
    vref::DecisionVariableRef,
    object::JuMP.AbstractVariable
    )::Nothing
    _data_object(vref).variable = object
    return
end

################################################################################
#                           VARIABLE INFO METHODS
################################################################################
# Get info
function _variable_info(vref::UserDecisionVariableRef)::JuMP.VariableInfo
    return _core_variable_object(vref).info
end

"""
    JuMP.has_lower_bound(vref::UserDecisionVariableRef)::Bool

Extend `JuMP.has_lower_bound` to return a `Bool` whether an `InfiniteOpt`
variable has a lower bound.

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

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

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

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

# Extend to return the index of the lower bound constraint associated with `vref`.
function _lower_bound_index(vref::UserDecisionVariableRef)::InfOptConstraintIndex
    if !JuMP.has_lower_bound(vref)
        error("Variable $(vref) does not have a lower bound.")
    end
    return _data_object(vref).lower_bound_index
end

# Extend to specify the index `cindex` of a lower bound constraint for `vref`.
function _set_lower_bound_index(
    vref::UserDecisionVariableRef,
    cindex::Union{InfOptConstraintIndex, Nothing}
    )::Nothing
    _data_object(vref).lower_bound_index = cindex
    return
end

"""
    JuMP.set_lower_bound(vref::UserDecisionVariableRef, lower::Real)::Nothing

Extend `JuMP.set_lower_bound` to specify the lower bound of an `InfiniteOpt`
variable `vref`. Errors if `vref` is fixed.

**Example**
```julia-repl
julia> set_lower_bound(vref, -1)

julia> lower_bound(vref)
-1.0
```
"""
function JuMP.set_lower_bound(
    vref::UserDecisionVariableRef,
    lower::Real
    )::Nothing
    newset = MOI.GreaterThan(convert(Float64, lower))
    gvref = _make_variable_ref(JuMP.owner_model(vref), JuMP.index(vref))
    new_constr = JuMP.ScalarConstraint(gvref, newset)
    model = JuMP.owner_model(vref)
    if JuMP.has_lower_bound(vref)
        cindex = _lower_bound_index(vref)
        cref = _make_constraint_ref(model, cindex)
        _set_core_constraint_object(cref, new_constr)
        set_optimizer_model_ready(model, false)
    else
        @assert !JuMP.is_fixed(vref) "$vref is fixed, cannot set lower bound."
        cref = JuMP.add_constraint(model, new_constr, is_info_constr = true)
        _set_lower_bound_index(vref, JuMP.index(cref))
    end
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(true, convert(Float64, lower),
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

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

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

**Example**
```julia-repl
var ≥ 0.0
```
"""
function JuMP.LowerBoundRef(vref::UserDecisionVariableRef)::InfOptConstraintRef
    cindex = _lower_bound_index(vref)
    model = JuMP.owner_model(vref)
    return _make_constraint_ref(model, cindex)
end

"""
    JuMP.delete_lower_bound(vref::UserDecisionVariableRef)::Nothing

Extend `JuMP.delete_lower_bound` to delete lower bound of `vref`. Errors if it
doesn't have a lower bound.

**Example**
```julia-repl
julia> delete_lower_bound(vref)

julia> has_lower_bound(vref)
false
```
"""
function JuMP.delete_lower_bound(vref::UserDecisionVariableRef)::Nothing
    JuMP.delete(JuMP.owner_model(vref), JuMP.LowerBoundRef(vref))
    _set_lower_bound_index(vref, nothing)
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(false, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

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

Extend `JuMP.has_upper_bound` to return a `Bool` whether an `InfiniteOpt`
variable has an upper bound.

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

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

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

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

# Extend to return the index of the upper bound constraint associated with `vref`.
function _upper_bound_index(vref::UserDecisionVariableRef)::InfOptConstraintIndex
    if !JuMP.has_upper_bound(vref)
        error("Variable $(vref) does not have a upper bound.")
    end
    return _data_object(vref).upper_bound_index
end

# Extend to specify the index `cindex` of a upper bound constraint for `vref`
function _set_upper_bound_index(
    vref::UserDecisionVariableRef,
    cindex::Union{InfOptConstraintIndex, Nothing}
    )::Nothing
    _data_object(vref).upper_bound_index = cindex
    return
end

"""
    JuMP.set_upper_bound(vref::UserDecisionVariableRef, upper::Real)::Nothing

Extend `JuMP.set_upper_bound` to specify the upper bound of an `InfiniteOpt`
variable `vref`. Errors if `vref` is fixed.

**Example**
```julia-repl
julia> set_upper_bound(vref, 1)

julia> upper_bound(vref)
1.0
```
"""
function JuMP.set_upper_bound(
    vref::UserDecisionVariableRef,
    upper::Real
    )::Nothing
    newset = MOI.LessThan(convert(Float64, upper))
    gvref = _make_variable_ref(JuMP.owner_model(vref), JuMP.index(vref))
    new_constr = JuMP.ScalarConstraint(gvref, newset)
    model = JuMP.owner_model(vref)
    if JuMP.has_upper_bound(vref)
        cindex = _upper_bound_index(vref)
        cref = _make_constraint_ref(model, cindex)
        _set_core_constraint_object(cref, new_constr)
        set_optimizer_model_ready(model, false)
    else
        @assert !JuMP.is_fixed(vref) "$vref is fixed, cannot set upper bound."
        cref = JuMP.add_constraint(model, new_constr, is_info_constr = true)
        _set_upper_bound_index(vref, JuMP.index(cref))
    end
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  true, convert(Float64, upper),
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

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

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

**Example**
```julia-repl
julia> cref = UpperBoundRef(vref)
var ≤ 1.0
```
"""
function JuMP.UpperBoundRef(vref::UserDecisionVariableRef)::InfOptConstraintRef
    cindex = _upper_bound_index(vref)
    model = JuMP.owner_model(vref)
    return _make_constraint_ref(model, cindex)
end

"""
    JuMP.delete_upper_bound(vref::UserDecisionVariableRef)::Nothing

Extend `JuMP.delete_upper_bound` to delete the upper bound of `vref`. Errors if
it doesn't have an upper bound.

**Example**
```julia-repl
julia> delete_upper_bound(vref)

julia> has_upper_bound(vref)
false
```
"""
function JuMP.delete_upper_bound(vref::UserDecisionVariableRef)::Nothing
    JuMP.delete(JuMP.owner_model(vref), JuMP.UpperBoundRef(vref))
    _set_upper_bound_index(vref, nothing)
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  false, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

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

Extend `JuMP.is_fixed` to return `Bool` whether an `InfiniteOpt` variable is
fixed.

**Example**
```julia-repl
julia> is_fixed(vref)
true
```
"""
JuMP.is_fixed(vref::UserDecisionVariableRef)::Bool = _variable_info(vref).has_fix

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

Extend `JuMP.fix_value` to return the fix value of an `InfiniteOpt` variable.
Errors if variable is not fixed.

**Example**
```julia-repl
julia> fix_value(vref)
0.0
```
"""
function JuMP.fix_value(vref::UserDecisionVariableRef)::Float64
    if !JuMP.is_fixed(vref)
        error("Variable $(vref) is not fixed.")
    end
    return _variable_info(vref).fixed_value
end

# Extend to return the index of the fix constraint associated with `vref`.
function _fix_index(vref::UserDecisionVariableRef)::InfOptConstraintIndex
    if !JuMP.is_fixed(vref)
        error("Variable $(vref) is not fixed.")
    end
    return _data_object(vref).fix_index
end

# Extend to set the index of the fix constraintassociated with `vref`.
function _set_fix_index(
    vref::UserDecisionVariableRef,
    cindex::Union{InfOptConstraintIndex, Nothing}
    )::Nothing
    _data_object(vref).fix_index = cindex
    return
end

"""
    JuMP.fix(vref::UserDecisionVariableRef, value::Real;
             force::Bool = false)::Nothing

Extend `JuMP.fix` to fix the value of an `InfiniteOpt` variable. Errors if
variable has a lower and/or an upper bound(s) unless `force = true`.

**Examples**
```julia-repl
julia> fix(vref, 3)

julia> fix_value(vref)
3.0

julia> fix(vref2, 2, force = true)

julia> fix_value(vref2)
2.0
```
"""
function JuMP.fix(
    vref::UserDecisionVariableRef,
    value::Real;
    force::Bool = false
    )::Nothing
    new_set = MOI.EqualTo(convert(Float64, value))
    model = JuMP.owner_model(vref)
    gvref = _make_variable_ref(JuMP.owner_model(vref), JuMP.index(vref))
    new_constr = JuMP.ScalarConstraint(gvref, new_set)
    if JuMP.is_fixed(vref)  # Update existing fixing constraint.
        cindex = _fix_index(vref)
        cref = _make_constraint_ref(model, cindex)
        _set_core_constraint_object(cref, new_constr)
        set_optimizer_model_ready(model, false)
    else  # Add a new fixing constraint.
        if  JuMP.has_upper_bound(vref) ||  JuMP.has_lower_bound(vref)
            if !force
                error("Unable to fix $(vref) to $(value) because it has " *
                      "existing variable bounds. Consider calling " *
                      "`JuMP.fix(variable, value; force=true)` which will " *
                      "delete existing bounds before fixing the variable.")
            end
            if  JuMP.has_upper_bound(vref)
                 JuMP.delete_upper_bound(vref)
            end
            if  JuMP.has_lower_bound(vref)
                 JuMP.delete_lower_bound(vref)
            end
        end
        cref = JuMP.add_constraint(model, new_constr, is_info_constr = true)
        _set_fix_index(vref, JuMP.index(cref))
    end
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(false, info.lower_bound,
                                                  false, info.upper_bound,
                                                  true, convert(Float64, value),
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

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

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

**Examples**
```julia-repl
julia> cref = FixRef(vref)
var = 1.0
```
"""
function JuMP.FixRef(vref::UserDecisionVariableRef)::InfOptConstraintRef
    cindex = _fix_index(vref)
    model = JuMP.owner_model(vref)
    return _make_constraint_ref(model, cindex)
end

"""
    JuMP.unfix(vref::UserDecisionVariableRef)::Nothing

Extend `JuMP.unfix` to unfix `vref`. Errors if it is not fixed.

**Example**
```julia-repl
julia> unfix(vref)

julia> is_fixed(vref)
false
```
"""
function JuMP.unfix(vref::UserDecisionVariableRef)::Nothing
    JuMP.delete(JuMP.owner_model(vref), JuMP.FixRef(vref))
    _set_fix_index(vref, nothing)
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  false, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, info.integer))
    return
end

"""
    JuMP.start_value(vref::UserDecisionVariableRef)::Union{Nothing, Float64}

Extend `JuMP.start_value` to return starting value of `InfiniteOpt` variable if
it has one. Returns `nothing` otherwise.

**Example**
```julia-repl
julia> start_value(vref)
0.0
```
"""
function JuMP.start_value(vref::UserDecisionVariableRef)::Union{Nothing, Float64}
    if _variable_info(vref).has_start
        return _variable_info(vref).start
    else
        return
    end
end

"""
    JuMP.set_start_value(vref::UserDecisionVariableRef, value::Real)::Nothing

Extend `JuMP.set_start_value` to specify the start value of `InfiniteOpt`
variables.

**Example**
```julia-repl
julia> set_start_value(vref, 1)

julia> start_value(vref)
1.0
```
"""
function JuMP.set_start_value(
    vref::UserDecisionVariableRef,
    value::Real
    )::Nothing
    info = _variable_info(vref)
    set_optimizer_model_ready(JuMP.owner_model(vref), false)
    _update_variable_info(vref,
                          JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                            info.has_ub, info.upper_bound,
                                            info.has_fix, info.fixed_value,
                                            true, Float64(value),
                                            info.binary, info.integer))
    return
end

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

Extend `JuMP.is_binary` to return `Bool` whether an `InfiniteOpt` variable is
binary.

**Example**
```julia-repl
true
```
"""
JuMP.is_binary(vref::UserDecisionVariableRef)::Bool = _variable_info(vref).binary

# Extend to return the index of the binary constraint associated with `vref`.
function _binary_index(vref::UserDecisionVariableRef)::InfOptConstraintIndex
    if !JuMP.is_binary(vref)
        error("Variable $(vref) is not binary.")
    end
    return _data_object(vref).zero_one_index
end

# Extend to specify the index of the binary constraint associated with `vref`.
function _set_binary_index(
    vref::UserDecisionVariableRef,
    cindex::Union{InfOptConstraintIndex, Nothing}
    )::Nothing
    _data_object(vref).zero_one_index = cindex
    return
end

"""
    JuMP.set_binary(vref::UserDecisionVariableRef)::Nothing

Extend `JuMP.set_binary` to specify an `InfiniteOpt` variable as a binary
variable. Errors if `vref` is an integer variable.

**Example**
```julia-repl
julia> set_binary(vref)

julia> is_binary(vref)
true
```
"""
function JuMP.set_binary(vref::UserDecisionVariableRef)::Nothing
    if JuMP.is_binary(vref)
        return
    elseif JuMP.is_integer(vref)
        error("Cannot set the variable_ref $(vref) to binary as it " *
              "is already integer.")
    end
    gvref = _make_variable_ref(JuMP.owner_model(vref), JuMP.index(vref))
    cref = JuMP.add_constraint(JuMP.owner_model(vref),
                               JuMP.ScalarConstraint(gvref, MOI.ZeroOne()),
                               is_info_constr = true)
    _set_binary_index(vref, JuMP.index(cref))
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  true, info.integer))
    return
end

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

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

**Example**
```julia-repl
julia> cref = BinaryRef(vref)
var binary
```
"""
function JuMP.BinaryRef(vref::UserDecisionVariableRef)::InfOptConstraintRef
    cindex = _binary_index(vref)
    model = JuMP.owner_model(vref)
    return _make_constraint_ref(model, cindex)
end

"""
    JuMP.unset_binary(vref::UserDecisionVariableRef)::Nothing

Extend `JuMP.unset_binary` to unset `vref` as a binary variable. Errors if it is
not binary.

```julia-repl
julia> unset_binary(vref)

julia> is_binary(vref)
false
```
"""
function JuMP.unset_binary(vref::UserDecisionVariableRef)::Nothing
    JuMP.delete(JuMP.owner_model(vref), JuMP.BinaryRef(vref))
    _set_binary_index(vref, nothing)
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  false, info.integer))
    return
end

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

Extend `JuMP.is_integer` to return `Bool` whether an `InfiniteOpt` variable is
integer.

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

# Extend to return the index of the integer constraintassociated with `vref`.
function _integer_index(vref::UserDecisionVariableRef)::InfOptConstraintIndex
    if !JuMP.is_integer(vref)
        error("Variable $(vref) is not an integer.")
    end
    return _data_object(vref).integrality_index
end

# Extend to specify the index of the integer constraint associated with `vref`.
function _set_integer_index(
    vref::UserDecisionVariableRef,
    cindex::Union{InfOptConstraintIndex, Nothing}
    )::Nothing
    _data_object(vref).integrality_index = cindex
    return
end

"""
    JuMP.set_integer(vref::UserDecisionVariableRef)::Nothing

Extend `JuMP.set_integer` to specify an `InfiniteOpt` variable as a integer
variable. Errors if `vref` is an binary variable.

**Example**
```julia-repl
julia> set_integer(vref)

julia> is_integer(vref)
true
```
"""
function JuMP.set_integer(vref::UserDecisionVariableRef)::Nothing
    if JuMP.is_integer(vref)
        return
    elseif JuMP.is_binary(vref)
        error("Cannot set the variable_ref $(vref) to integer as it " *
              "is already binary.")
    end
    gvref = _make_variable_ref(JuMP.owner_model(vref), JuMP.index(vref))
    cref = JuMP.add_constraint(JuMP.owner_model(vref),
                               JuMP.ScalarConstraint(gvref, MOI.Integer()),
                               is_info_constr = true)
    _set_integer_index(vref, JuMP.index(cref))
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, true))
    return
end

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

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

**Example**
```julia-repl
julia> cref = IntegerRef(vref)
var integer
```
"""
function JuMP.IntegerRef(vref::UserDecisionVariableRef)::InfOptConstraintRef
    cindex = _integer_index(vref)
    model = JuMP.owner_model(vref)
    return _make_constraint_ref(model, cindex)
end

"""
    JuMP.unset_integer(vref::UserDecisionVariableRef)::Nothing

Extend `JuMP.unset_integer` to unset `vref` as an integer variable. Errors if it
is not an integer variable.

```julia-repl
julia> unset_integer(vref)

julia> is_integer(vref)
false
```
"""
function JuMP.unset_integer(vref::UserDecisionVariableRef)::Nothing
    JuMP.delete(JuMP.owner_model(vref), JuMP.IntegerRef(vref))
    _set_integer_index(vref, nothing)
    info = _variable_info(vref)
    _update_variable_info(vref, JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                                  info.has_ub, info.upper_bound,
                                                  info.has_fix, info.fixed_value,
                                                  info.has_start, info.start,
                                                  info.binary, false))
    return
end

################################################################################
#                          MODEL VARIABLE QUERIES
################################################################################
# Define finite variable struct for queries
struct FiniteVariable end

"""
    JuMP.num_variables(model::InfiniteModel, [type])::Int

Extend `JuMP.num_variables` to return the number of `InfiniteOpt` variables
assigned to `model`. By default, the total number of infinite, semi-infinite,
point, and finite variables is returned. The amount of a particular type is
obtained by specifying the concrete variable type via `type`. Type options
include:
 - `InfiniteVariable`: all infinite variables
 - `SemiInfiniteVariable`: all semi-infinite variables
 - `PointVariable`: all point variables
 - `FiniteVariable`: all finite variables

**Example**
```julia-repl
julia> num_variables(model)
3

julia> num_variables(model, InfiniteVariable)
2
```
"""
function JuMP.num_variables(model::InfiniteModel)::Int
    num_vars = JuMP.num_variables(model, InfiniteVariable)
    num_vars += JuMP.num_variables(model, SemiInfiniteVariable)
    num_vars += JuMP.num_variables(model, PointVariable)
    num_vars += JuMP.num_variables(model, FiniteVariable)
    return num_vars
end

# Particular variable types
function JuMP.num_variables(model::InfiniteModel,
                            type::Type{C})::Int where {C}
    return length(_data_dictionary(model, type))
end

"""
    JuMP.all_variables(model::InfiniteModel, [type])::Vector{GeneralVariableRef}

Extend `JuMP.all_variables`] to return a list of all the variable references
associated with `model`. By default, all of the infinite, semi-infinite, point,
and finite variables is returned. Those of a particular type is obtained by
specifying the concrete variable type via `type`. Type options include:
 - `InfiniteVariable`: all infinite variables
 - `SemiInfiniteVariable`: all semi-infinite variables
 - `PointVariable`: all point variables
 - `FiniteVariable`: all finite variables

**Examples**
```julia-repl
julia> all_variables(model)
4-element Array{GeneralVariableRef,1}:
 y(t)
 w(t, x)
 y(0)
 z

julia> all_variables(model, PointVariable)
1-element Array{GeneralVariableRef,1}:
 y(0)
```
"""
function JuMP.all_variables(model::InfiniteModel)::Vector{GeneralVariableRef}
    vrefs_list = JuMP.all_variables(model, InfiniteVariable)
    append!(vrefs_list, JuMP.all_variables(model, SemiInfiniteVariable))
    append!(vrefs_list, JuMP.all_variables(model, PointVariable))
    append!(vrefs_list, JuMP.all_variables(model, FiniteVariable))
    return vrefs_list
end

# Particular variable types
function JuMP.all_variables(model::InfiniteModel,
                            type::Type{C}
                            )::Vector{GeneralVariableRef} where {C}
    vrefs_list = Vector{GeneralVariableRef}(undef, JuMP.num_variables(model, type))
    for (i, (index, _)) in enumerate(_data_dictionary(model, type))
        vrefs_list[i] = _make_variable_ref(model, index)
    end
    return vrefs_list
end

################################################################################
#                                 DELETION
################################################################################
# Helper function to delete the info constraints
function _delete_info_constraints(vref::UserDecisionVariableRef)::Nothing
    # remove variable info constraints associated with vref
    if JuMP.has_lower_bound(vref)
        JuMP.delete_lower_bound(vref)
    end
    if JuMP.has_upper_bound(vref)
        JuMP.delete_upper_bound(vref)
    end
    if JuMP.is_fixed(vref)
        JuMP.unfix(vref)
    end
    if JuMP.is_binary(vref)
        JuMP.unset_binary(vref)
    elseif JuMP.is_integer(vref)
        JuMP.unset_integer(vref)
    end
    return
end

"""
    JuMP.delete(model::InfiniteModel, vref::DecisionVariableRef)::Nothing

Extend `JuMP.delete` to delete `InfiniteOpt` variables and their dependencies.
Errors if variable is invalid, meaning it has already been deleted or it belongs
to another model.

**Example**
```julia-repl
julia> print(model)
Min measure(g(t)*t) + z
Subject to
 z ≥ 0.0
 g(t) + z ≥ 42.0, ∀ t ∈ [0, 6]
 g(0.5) = 0

julia> delete(model, g)

julia> print(model)
Min measure(t) + z
Subject to
 z ≥ 0.0
 z ≥ 42.0
```
"""
function JuMP.delete(model::InfiniteModel, vref::DecisionVariableRef)::Nothing
    @assert JuMP.is_valid(model, vref) "Variable is invalid."
    # update the optimizer model status
    if is_used(vref)
        set_optimizer_model_ready(model, false)
    end
    # delete attributes specific to the variable type
    _delete_variable_dependencies(vref)
    gvref = _make_variable_ref(model, JuMP.index(vref))
    # remove from measures if used
    for mindex in _measure_dependencies(vref)
        mref = dispatch_variable_ref(model, mindex)
        func = measure_function(mref)
        data = measure_data(mref)
        if func isa GeneralVariableRef
            new_func = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
            new_meas = Measure(new_func, data, Int[], Int[], true)
        else
            _remove_variable(func, gvref)
            new_meas = build_measure(func, data)
        end
        _set_core_variable_object(mref, new_meas)
    end
    # remove from constraints if used
    for cindex in _constraint_dependencies(vref)
        cref = _make_constraint_ref(model, cindex)
        func = JuMP.jump_function(JuMP.constraint_object(cref))
        if func isa GeneralVariableRef
            set = JuMP.moi_set(JuMP.constraint_object(cref))
            new_func = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
            new_constr = JuMP.ScalarConstraint(new_func, set)
            _set_core_constraint_object(cref, new_constr)
            empty!(_object_numbers(cref))
        elseif func isa AbstractArray{GeneralVariableRef}
            JuMP.delete(model, cref)
        else
            _remove_variable(func, gvref)
            # update the object numbers if vref is infinite
            if vref isa Union{InfiniteVariableRef, SemiInfiniteVariableRef}
                _data_object(cref).object_nums = sort(_object_numbers(func))
            end
        end
    end
    # remove from objective if vref is in it
    if used_by_objective(vref)
        if JuMP.objective_function(model) isa GeneralVariableRef
            new_func = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
            JuMP.set_objective_function(model, new_func)
            JuMP.set_objective_sense(model, MOI.FEASIBILITY_SENSE)
        else
            _remove_variable(JuMP.objective_function(model), gvref)
        end
    end
    # delete the variable information
    _delete_data_object(vref)
    return
end

################################################################################
#                           INTEGRALITY RELAXATION
################################################################################
"""
    JuMP.relax_integrality(model::InfiniteModel)::Function

Modifies `model` to "relax" all binary and integrality constraints on
variables. Specifically,

- Binary constraints are deleted, and variable bounds are tightened if
  necessary to ensure the variable is constrained to the interval ``[0, 1]``.
- Integrality constraints are deleted without modifying variable bounds.
- All other constraints are ignored (left in place). This includes discrete
  constraints like SOS and indicator constraints.

Returns a function that can be called without any arguments to restore the
original model. The behavior of this function is undefined if additional
changes are made to the affected variables in the meantime.

**Example**
```julia-repl
julia> undo_relax = relax_integrality(model);

julia> print(model)
Min x + ∫{t ∈ [0, 10]}(y(t))
Subject to
 x ≥ 0.0
 y(t) ≥ 1.0
 x ≤ 1.0
 y(t) ≤ 10.0

julia> undo_relax()

julia> print(model)
Min x + ∫{t ∈ [0, 10]}(y(t))
Subject to
 y(t) ≥ 1.0
 y(t) ≤ 10.0
 y(t) integer
 x binary
```
"""
function JuMP.relax_integrality(model::InfiniteModel)
    # TODO ensure variables are not semi-continous/integer
    # get the variable info
    info_pre_relaxation = map(v -> (dispatch_variable_ref(v),
                              _variable_info(dispatch_variable_ref(v))),
                              JuMP.all_variables(model))
    # relax the appropriate variables
    for (vref, info) in info_pre_relaxation
        if info.integer
            JuMP.unset_integer(vref)
        elseif info.binary
            JuMP.unset_binary(vref)
            if isnan(info.lower_bound)
                JuMP.set_lower_bound(vref, 0.0)
            else
                JuMP.set_lower_bound(vref, max(0.0, info.lower_bound))
            end
            if isnan(info.upper_bound)
                JuMP.set_upper_bound(vref, 1.0)
            else
                JuMP.set_upper_bound(vref, min(1.0, info.upper_bound))
            end
        end
    end
    function unrelax()::Nothing
        for (vref, info) in info_pre_relaxation
            if info.integer
                JuMP.set_integer(vref)
            elseif info.binary
                JuMP.set_binary(vref)
                if info.has_lb
                    JuMP.set_lower_bound(vref, info.lower_bound)
                else
                    JuMP.delete_lower_bound(vref)
                end
                if info.has_ub
                    JuMP.set_upper_bound(vref, info.upper_bound)
                else
                    JuMP.delete_upper_bound(vref)
                end
            end
        end
        return
    end
    return unrelax
end


text:
using OctreeBH
using StaticArrays
using Test

const BOXSIZE = 1.0

const N = 3 #spatial dimension
const T = Float64

Npart = 1000 #number of particles

@testset "gravtree" begin

#randomly distributing particles
X = [@SVector rand(N) for _ in 1:Npart]

topnode_length = @SVector(ones(N)) * BOXSIZE  #actual length of tree

center = topnode_length .* 0.5

#individual smoothing length (only relevant for scatter ngbs)
hsml = zeros(Npart)

mass = ones(Npart) #particle mass

part = [Data{N,T}(SVector(X[i]), i, hsml[i], mass[i]) for i in eachindex(X)]

#build the tree
tree = buildtree(part, center, topnode_length);

#search center
#x = @SVector rand(N)

boxsizes = @SVector(ones(N)) * BOXSIZE  #for periodic B.C.

#softening = BOXSIZE / Npart^(1/3.)
softening = 0.0
ANGLE = 0.0

pot = zeros(Npart)
pot_bf = zeros(Npart)
for i in eachindex(X)
    ga = GravTreeGather{N,T}()
    gravity_treewalk!(ga,X[i],tree,ANGLE,softening,boxsizes)
    pot[i] = ga.pot

    ga_bf = GravTreeGather{N,T}()
    gravity_bruteforce!(ga_bf,X[i],X,mass,boxsizes)
    pot_bf[i] = ga_bf.pot
end
@test all(pot .≈ pot_bf)

end

ext: jl


text:
function PipeCoefsOver(wn_dict::Dict{Any,Any}, Q_base::Dict{String,Array{Float64,1}})
    pipes = wn_dict["pipes"]
    aPipe = Dict{String, Array{Float64,1}}()
    bPipe = Dict{String, Array{Float64,1}}()

    for (key, pipe) in pipes
        name = pipe["name"]
        Q_base_current = Q_base[name]
        Q_power_current = Q_base_current.^1.852
        aPipe_current = Vector{Union{Nothing,Float64}}(nothing, length(Q_base_current)-1)
        bPipe_current = Vector{Union{Nothing,Float64}}(nothing, length(Q_base_current)-1)
        for i in range(1, length(aPipe_current))
            aPipe_current[i] = (Q_power_current[i+1]-Q_power_current[i])/(Q_base_current[i+1]-Q_base_current[i])
            bPipe_current[i] =Q_power_current[i] - aPipe_current[i]*Q_base_current[i]
        end
        aPipe[name] = aPipe_current
        bPipe[name] = bPipe_current
    end
    return(Q_base, aPipe, bPipe)
end

function PipeCoefsUnder(wn_dict::Dict{Any,Any}, Q_base::Dict{String,Array{Float64,1}})
    pipes = wn_dict["pipes"]
    aPipe = Dict{String, Array{Float64,1}}()
    bPipe = Dict{String, Array{Float64,1}}()
    for (key, pipe) in pipes
        name = pipe["name"]
        Q_base_current = Q_base[name]
        Q_power_current = Array{Union{Nothing,Float64}}(nothing,length(Q_base_current))
        aPipe_current = Array{Union{Nothing,Float64}}(nothing,length(Q_base_current)) #note # of slopes = # of points unlike over estimate
        bPipe_current = Array{Union{Nothing,Float64}}(nothing,length(Q_base_current))
        for i = 1:length(Q_base_current)
            Q_power_current[i] = Q_base_current[i]^1.852
            bPipe_current[i] = 1.852 * Q_base_current[i]^0.852 #slope of tangent line
            aPipe_current[i] = Q_power_current[i] - bPipe_current[i]*Q_base_current[i]
        end
        aPipe[name] = aPipe_current
        bPipe[name] = bPipe_current
    end
    return(Q_base, aPipe, bPipe)
end

function Q_base_pipe(wn_dict::Dict{Any,Any}, Qmin::Dict{String,Float64}, Qmax::Dict{String,Float64}, n::Int64, Q_lb::Float64, log_density::Union{Int64,Float64}, density_boost::Union{Int64, Float64}, lin_density::Union{Float64,Int64}, linlog_cut::Union{Float64,Int64},dH_critical::Float64)
    pipes = wn_dict["pipes"]
    Q_base = Dict{String, Array{Float64,1}}()
    for (key, pipe) in pipes
        name = pipe["name"]
        Q_bar = maximum([abs(Qmin[name]), abs(Qmax[name])])
        if Q_bar <= Q_lb
            Q_base[name] = [0, Q_lb]
        else
            HW_multiplier = (1/pipe["roughness"])^1.852 *10.67 *pipe["length"]*(1/pipe["diameter"])^4.87
            if Q_bar ^1.852*HW_multiplier >= dH_critical
                log_density_current = log_density*density_boost
                log_density_current = Int(floor(log10(Q_bar*linlog_cut/Q_lb)*log_density_current))+1
                lin_density_current = lin_density* density_boost
            else
                log_density_current = log_density_current
                log_density_current = Int(floor(log10(Q_bar *linlog_cut/Q_lb)*log_density_current)) +1
                lin_density_current = lin_density
            end
            Q_base_log = exp10.(range(log10(Q_lb), log10(Q_bar*linlog_cut), step = log_density_current))
            Q_base_lin = exp10.(range(Q_bar*linlog_cut, Q_bar, step = lin_density))
            Q_base[name] = [0; Q_base_log; Q_base_lin[1:length(Q_base_lin)]]
        end
        Q_base_current = Q_base[name]
        length(Q_base[name]) > n+1 ? error("Please Choose another n > $length(Q_base[name])") : nothing
    end
    return(Q_base)
end


text:
# X and Z components in fT
# conducting 100 Ωm
## model
zfixed   = [-1e6,   0, 40, 60]
rho      = [1e12,   100, 100, 100]

ross_primary = [25.7847     -18.7313]
ross_secondary = [3.9416     -5.8586
    1.8035     -3.6031
    1.1267     -2.6355
    0.6896     -1.8802
   0.38484     -1.2491
   0.20371    -0.79275
   0.10128    -0.47781
  0.049214    -0.28167
  0.023819    -0.16467
  0.011302    -0.09436
0.0051838   -0.052416
0.0022996    -0.02822
0.0009899   -0.014747
0.00041294  -0.0075244
0.00015465  -0.0036212]

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
create(c::Client, ::Type{Invite}, ch::DiscordChannel; kwargs...) = create_channel_invite(c, ch.id; kwargs...)

retrieve(c::Client, ::Type{Invite}, g::AbstractGuild) = get_guild_invites(c, g.id)
retrieve(c::Client, ::Type{Invite}, ch::DiscordChannel) = get_channel_invites(c, ch.id)
retrieve(c::Client, ::Type{Invite}, invite::AbstractString; kwargs...) = get_invite(c, invite; kwargs...)

delete(c::Client, i::Invite) = delete_invite(c, i.code)


text:
using Pkg
Pkg.activate(".")

using Reproduce
using Logging

const save_loc = "markov_chain_results"
const exp_file = "experiment/markov_chain.jl"
const exp_module_name = :MarkovChainExperiment
const exp_func_name = :main_experiment
# const alphas = collect(0.0:0.1:12)
const u_0 = 10 .^ (0:0.25:3)
const alphas = [0.0, 0.001, 0.1, 1.0]
const policies = ["easy", "hard", "hardest"]

function make_arguments(args::Dict)
    new_args=["--policy", args["policy"]]
    return new_args
end

function main()

    as = ArgParseSettings()
    @add_arg_table as begin
        "numworkers"
        arg_type=Int64
        default=1
        "--numruns"
        arg_type=Int64
        default=100
        "--saveloc"
        arg_type=String
        default=save_loc
        "--numjobs"
        action=:store_true
    end
    parsed = parse_args(as)
    num_workers = parsed["numworkers"]

    arg_dict = Dict([
        # "outhorde"=>["onestep", "chain"],
        "policy"=>policies
    ])
    arg_list = ["policy"]
    alg_list = ["--wisrupam", "--normalize_eta"]
    static_args = [alg_list; ["--exp_loc", parsed["saveloc"], "--numruns", string(parsed["numruns"]), "--numiter", "1000", "--alphas", string.(alphas)..., "--init_u", string.(u_0)...]]
    args_iterator = ArgIterator(arg_dict, static_args; arg_list=arg_list, make_args=make_arguments)

    if parsed["numjobs"]
        @info "This experiment has $(length(collect(args_iterator))) jobs."
        println(collect(args_iterator)[num_workers])
        exit(0)
    end

    # experiment = Experiment(save_loc)
    experiment = Experiment(save_loc,
                            exp_file,
                            exp_module_name,
                            exp_func_name,
                            args_iterator)

    create_experiment_dir(experiment)
    add_experiment(experiment; settings_dir="settings")
    ret = job(experiment; num_workers=num_workers)
    post_experiment(experiment, ret)
end

main()


text:
# Inspired from `SetProg/test/Test/square.jl`

using LinearAlgebra, Test

using SwitchOnSafety
using Polyhedra
const Sets = SetProg.Sets
using MultivariatePolynomials

function invariant_square_test(
    factory, variable::SetProg.AbstractVariable, set_test; kws...)
    □ = polyhedron(HalfSpace([1, 0], 1.0) ∩ HalfSpace([-1, 0], 1) ∩ HalfSpace([0, 1], 1) ∩ HalfSpace([0, -1], 1))
    A = [0.0 -1.0
         1.0  0.0]
    system = ConstrainedLinearDiscreteSystem(A, □)

    set = invariant_set(system, factory, variable; verbose=0, kws...)
    set_test(set)
end

const atol = 1e-4
const rtol = 1e-4

@testset "Ellipsoid homogeneous with $factory" for factory in sdp_factories
    invariant_square_test(
        factory, Ellipsoid(symmetric=true, dimension=2),
        ◯ -> begin
            @test ◯ isa Sets.Polar{Float64, Sets.EllipsoidAtOrigin{Float64}}
            @test Sets.polar(◯).Q ≈ Symmetric([1.0 0.0; 0.0 1.0]) atol=atol rtol=rtol
        end,
        volume_heuristic = nth_root)
end

@testset "Convex quadratic homogeneous with $factory" for factory in sdp_factories
    # cholesky miss condition :: not positive definite :: info = 27 :: line 785 in sdpa_linear.cpp
    issdpa(factory) && continue
    invariant_square_test(
        factory, PolySet(degree=2, convex=true, symmetric=true),
        ◯ -> begin
            @test ◯ isa Sets.Polar{Float64, Sets.ConvexPolynomialSublevelSetAtOrigin{Float64}}
            x, y = Sets.space_variables(◯)
            ◯_polar = Sets.polar(◯)
            @test ◯_polar.p ≈ x^2 + y^2 atol=atol rtol=rtol
        end,
        volume_heuristic = nth_root)
end


text:
struct DisjointPacker{Mtd, Mtd′}
    mtd::Mtd
    mtd′::Mtd′
end

function pack_plaid(A::SparseMatrixCSC, method::DisjointPacker; kwargs...)
    Φ = pack_stripe(A, method.mtd; kwargs...)
    Π = pack_stripe(PermutedDimsArray(A, (2, 1)), method.mtd′, Φ; kwargs...)
    return (Π, Φ)
end

struct AlternatingPacker{Mtds}
    mtds::Mtds
end

AlternatingPacker(mtds...) = AlternatingPacker{typeof(mtds)}(mtds)

function pack_plaid(A::SparseMatrixCSC, method::AlternatingPacker; adj_A = nothing, kwargs...)
    if adj_A === nothing
        adj_A = adjointpattern(A)
    end
    Φ = pack_stripe(A, method.mtds[1]; adj_A=adj_A, kwargs...)
    Π = pack_stripe(adj_A, method.mtds[2], Φ; adj_A=A, kwargs...)
    for (i, mtd) in enumerate(method.mtds[3:end])
        if isodd(i)
            Φ = pack_stripe(A, mtd, Π; adj_A=adj_A, kwargs...)
        else
            Π = pack_stripe(adj_A, mtd, Φ; adj_A=A, kwargs...)
        end
    end
    return (Π, Φ)
end

struct SymmetricPacker{Mtds}
    mtds::Mtds
end

SymmetricPacker(mtds...) = SymmetricPacker{typeof(mtds)}(mtds)

function pack_plaid(A::SparseMatrixCSC, method::SymmetricPacker; adj_A = nothing, kwargs...)
    if adj_A === nothing
        adj_A = adjointpattern(A)
    end
    Π = pack_stripe(A, method.mtds[1]; adj_A=adj_A, kwargs...)
    for (i, mtd) in enumerate(method.mtds[2:end])
        if isodd(i)
            Π = pack_stripe(A, mtd, Π; adj_A=adj_A, kwargs...)
        else
            Π = pack_stripe(adj_A, mtd, Π; adj_A=A, kwargs...)
        end
    end
    return (Π, Π)
end


text:
module RandomProjectionTrees

using StatsBase: sample
using LinearAlgebra: adjoint, norm

include("tree.jl")
include("split.jl")
include("utils.jl")

export RandomProjectionTree, search_rptree, RandomProjectionForest, approx_knn

end # module


text:
#=
## Generate `PowerDynamics.jl`-Node with `BlockSystems`

We want to model the
[`VSIVoltagePT1`](https://juliaenergy.github.io/PowerDynamics.jl/dev/node_types/#PowerDynamics.VSIVoltagePT1)
with the help of `BlockSystems`.

We start by defining some general stuff...
=#
using BlockSystems
using ModelingToolkit

@parameters t
D = Differential(t)
nothing #hide

#=
### Defining the common Blocks
We need some common parts such as filters or integrators. The actual Blocks will
be defined by using these 'blueprints'. Stuff like this should go to a library of
components at some point.

=#
# #### low pass filter
@parameters τ input(t)
@variables filtered(t)

lpf = IOBlock([D(filtered) ~ 1/τ * (- filtered + input)],
              [input], [filtered])

# #### voltage source

@parameters ω(t) v(t) τ
@variables u_i(t) u_r(t) A(t)

## explicit algebraic equation for A will be reduced at connect
voltage_source = IOBlock([A ~ 1/τ * (v/√(u_i^2 + u_r^2) - 1),
                          D(u_r) ~ -ω * u_i + A*u_r,
                          D(u_i) ~  ω * u_r + A*u_i],
                         [ω, v], [u_i, u_r])

# #### Droop control

@parameters K u_ref x_ref x(t)
@variables u(t)

droop_control = IOBlock([
    u ~ - K * (x - x_ref) + u_ref # output is the droop voltage v
    ], [x], [u])

#=
After defining the blueprints we can create the Blocks of the system
based on them:
=#

p_filter = IOBlock(lpf, name = :p_filter)
q_filter = IOBlock(lpf, name = :q_filter)
p_droop = IOBlock(droop_control, name = :p_droop)
q_droop = IOBlock(droop_control, name = :q_droop)
v_source = IOBlock(voltage_source, name = :v_source)

nothing #hide
#=
We can put the blocks together to form this system:

```
     +----------+  +-------------+
 P --| p_filter |--|   p_droop   |   +----------+
     |    τ     |  | u_ref, xref |---|          |
     +----------+  +-------------+  ω| v_source |-- u_r
                                     |    τ     |
     +----------+  +-------------+  v|          |-- u_i
 Q --| q_filter |--|   q_droop   |---|          |
     |    τ     |  | u_ref, xref |   +----------+
     +----------+  +-------------+

```
=#
gfi = IOSystem([p_filter.filtered => p_droop.x,
                q_filter.filtered => q_droop.x,
                p_droop.u => v_source.ω,
                q_droop.u => v_source.v],
               [p_filter, q_filter, p_droop, q_droop, v_source],
               name = :GridForming,
               namespace_map = [p_filter.input => :P_in,
                                q_filter.input => :Q_in,
                                p_filter.filtered => :p_filtered,
                                q_filter.filtered => :q_filtered,
                                ## parameter names which match VSIVoltagePT1
                                v_source.τ => :τ_v, # time constant voltage control delay
                                p_filter.τ => :τ_P, # time constant active power measurement
                                q_filter.τ => :τ_Q, # time constant reactive power measurement
                                p_droop.K  => :K_P, # droop constant frequency droop
                                q_droop.K  => :K_Q, # droop constant voltage droop
                                q_droop.u_ref => :V_r, # reference/ desired voltage
                                p_droop.u_ref => :ω_r, # reference/ desired frequency
                                p_droop.x_ref => :P, # active (real) power infeed
                                q_droop.x_ref => :Q], # reactive (imag) power infeed                .
               outputs = [v_source.u_i, v_source.u_r])

#=
There is still one problem: our connected system has the inputs `P` and `Q`. In order to construct
IONodes the IOSystems have to be systems which convert some complex current to complex voltage:
```
       +---------+
i_r -->|   ???   |--> u_r
i_i -->|         |--> u_i
       +---------+
```

We can achieve this by defining a nother block which converts `(i, u) ↦ (P, Q)`

```
               +----------+  +-------------+
     +-----+   | p_filter |--|   p_droop   |   +----------+
i_r--|     |-P-|   p_τ    |  |   P, ω_r    |---|          |
i_i--| pow |   +----------+  +-------------+  ω| v_source |---+--u_r
     |     |                                   |    τ     |   |
 +---|     |   +----------+  +-------------+  v|          |-+-|--u_i
 | +-|     |-Q-| q_filter |--|   q_droop   |---|          | | |
 | | +-----+   |   q_τ    |  |   Q, v_r    |   +----------+ | |
 | |           +----------+  +-------------+                | |
 | +--------------------------------------------------------+ |
 +-----------------------------------------------------------+

```
=#
@parameters u_i(t) u_r(t) i_i(t) i_r(t)
@variables P_in(t) Q_in(t)
pow = IOBlock([P_in ~ u_r*i_r + u_i*i_i,
               Q_in ~ u_i*i_r - u_r*i_i],
              [u_i, u_r, i_i, i_r], [P_in, Q_in], name=:pow)

gfi2 = IOSystem([gfi.u_i => pow.u_i,
                 gfi.u_r => pow.u_r,
                 pow.P_in => gfi.P_in,
                 pow.Q_in => gfi.Q_in],
                [gfi, pow], outputs=[gfi.u_i, gfi.u_r], name=:gfi_i_to_u)
nothing #hide

# and the system can be reduced to a single IOBlock
connected = connect_system(gfi2)
nothing #hide

#=
The following part should live inside `PowerDynamics.jl` at some point. It defines an
`IONode` which will be constructed from an `IOBlock` as a new `AbstractNode` type.
=#
using PowerDynamics
using NetworkDynamics: ODEVertex
struct IONode{T} <: AbstractNode
    block::IOBlock
    generated::T
    parameters::Vector{Float64}
end

function IONode(blk::IOBlock, parameters::Dict)
    ## BlockSpec: blk must by of type (i_r, i_i) ↦ (u_r, u_i)
    spec = BlockSpec([:i_r, :i_i], [:u_r, :u_i])
    @assert spec(blk) "Block has to follow PowerDynamics i/o conventions!"
    ## TODO check parameters
    gen = generate_io_function(blk,
                               f_states=[blk.u_r, blk.u_i],
                               f_inputs=[blk.i_r, blk.i_i],
                               f_params=keys(parameters), warn=false);
    IONode(blk, gen, collect(values(parameters)))
end

import PowerDynamics.construct_vertex
function construct_vertex(ion::IONode)
    gen = ion.generated
    function rhs!(dx, x, e_s, e_d, _p, t)
        i = total_current(e_s, e_d)
        gen.f_ip(dx, x, (real(i), imag(i)), ion.parameters, t)
    end
    ODEVertex(f = rhs!, dim = length(gen.states), mass_matrix = gen.massm, sym = Symbol.(gen.states))
end

import PowerDynamics.symbolsof
PowerDynamics.symbolsof(ionode::IONode) = Symbol.(ionode.generated.states)
PowerDynamics.dimension(ionode::IONode) = length(ionode.generated.states)

#=
Let's try it out:
The parameters have to be provided as a Dict. This dict should contain all internal parameters of the block.
The list of parameters can be retrieved from the block:
=#
@info "connected system:" connected.iparams

#=
The parameter dict can be defined in several ways, I think the first option is the most convenient
because it is a pure `Symbol` which does not depend on any variables in the scope. The namespace separator
can be typed as `\_+<tab>`.
```
:τ_P => 1.0
:gfi_i_to_u₊τ_P => 1.0
connected.τ_P => 1.0
```
=#
para = Dict(
    :τ_v => 0.005,    # time constant voltage control delay
    :τ_P => 0.5,      # time constant active power measurement
    :τ_Q => 0.5,      # time constant reactive power measurement
    :K_P => 0.396,     # droop constant frequency droop
    :K_Q => 0.198,     # droop constant voltage droop
    :V_r => 1.0,   # reference/ desired voltage
    :P   => 0.303, # active (real) power infeed
    :Q   => 0.126, # reactive (imag) power infeed                .
    :ω_r => 0.0)   # refrence/ desired frequency
nothing #hide

# Now let's test whether this works in PD...
using PowerDynamics: SlackAlgebraic, PiModelLine,PowerGrid,find_operationpoint
using OrderedCollections: OrderedDict

buses=OrderedDict(
    "bus1"=> IONode(connected, para),
    "bus2"=> SlackAlgebraic(U=1))

branches=OrderedDict(
    "branch1"=> PiModelLine(from= "bus1", to = "bus2",y=4.999131600798035-1im*15.263086523179553, y_shunt_km=0.0528/2, y_shunt_mk=0.0528/2))

powergrid = PowerGrid(buses,branches)
operationpoint = find_operationpoint(powergrid)
timespan= (0.0,0.1)
nothing #hide

# simulating a voltage perturbation at node 1
fault1 = ChangeInitialConditions(node="bus1", var=:u_r, f=Inc(0.2))
solution1 = simulate(fault1, powergrid, operationpoint, timespan)
using Plots
plot(solution1.dqsol)

#=
## Compare against `VSIVoltagePT1`
We  want to compare the newly defined node against the `VSIVoltagePT1` node for random data.
Both nodes have slightly different states:
 - `symbolsof(node_bs) == [:u_r, :u_i, :p_filtered, :q_filtered]`
 - `symbolsof(node_pd) == [:u_r, :u_i, :ω, :q_m]`
We need adapt the initial state for the different third variable:
```math
\omega = \omega_r - K_P\cdot(p_{\mathrm{filtered}} - P)\\
\dot\omega = - K_P \cdot \dot{p}_{\mathrm{filtered}}
```
=#
using Test
@testset "Compare IONode vs. PD Node" begin
    for i in 1:100
        para = Dict(
            :τ_v => rand(), # time constant voltage control delay
            :τ_P => rand(), # time constant active power measurement
            :τ_Q => rand(), # time constant reactive power measurement
            :K_P => rand(), # droop constant frequency droop
            :K_Q => rand(), # droop constant voltage droop
            :V_r => rand(), # reference/ desired voltage
            :P   => rand(), # active (real) power infeed
            :Q   => rand(), # reactive (imag) power infeed                .
            :ω_r => 0.0)    # refrence/ desired frequency

        ## create IONode
        node_bs = IONode(connected, para)
        f_bs = construct_vertex(node_bs).f

        ## create PDNode
        ## the PD node does not have the explicit ω_r parameter
        para_pd = delete!(copy(para), :ω_r)
        nt = NamedTuple{Tuple(keys(para_pd))}(values(para_pd))
        node_pd = VSIVoltagePT1(; nt...)
        f_pd = construct_vertex(node_pd).f

        ## create fake "edge data", 4 incoming, 4 outgooing with 4 states each
        es = [randn(4) for i in 1:4]
        ed = [randn(4) for i in 1:4]

        ## select random time
        t = rand()

        ## chose random initial state and account for initial ω in PD node
        x_bs = randn(4)
        x_pd = copy(x_bs)
        x_pd[3] = - para[:K_P] * (x_bs[3] - para[:P])

        ## create arrays for the results
        dx_bs = similar(x_bs)
        dx_pd = similar(x_pd)

        ## call both functions
        f_bs(dx_bs, x_bs, es, ed, nothing, t)
        f_pd(dx_pd, x_pd, es, ed, nothing, t)

        ## compare results
        ## we have to correct variable 3 of bs implementation to match dω
        @test dx_bs[1] ≈ dx_pd[1]                # u_r
        @test dx_bs[2] ≈ dx_pd[2]                # u_i
        @test - para[:K_P] * dx_bs[3] ≈ dx_pd[3] # ω
        @test dx_bs[4] ≈ dx_pd[4]                # q_filtered
    end
end
nothing #hide

#=
it works 🎉

## Benchmark
We can also run a quick benchmark of both node functions:
=#
para = Dict(:τ_v=>rand(),:τ_P=>rand(), :τ_Q=>rand(),
            :K_P=>rand(), :K_Q=>rand(), :V_r=>rand(),
            :P=>rand(), :Q=>rand(), :ω_r=>0.0)

node_bs = IONode(connected, para)
f_bs = construct_vertex(node_bs).f
para_pd = delete!(copy(para), :ω_r)
nt = NamedTuple{Tuple(keys(para_pd))}(values(para_pd))
node_pd = VSIVoltagePT1(; nt...)
f_pd = construct_vertex(node_pd).f

es = [randn(4) for i in 1:4]
ed = [randn(4) for i in 1:4]
t = rand()

## chose random initial state and account for initial ω in PD node
x_bs = randn(4)
x_pd = copy(x_bs)
x_pd[3] = - para[:K_P] * (x_bs[3] - para[:P])
dx = similar(x_bs)

using BenchmarkTools
@btime $f_bs($dx, $x_bs, $es, $ed, nothing, $t)
@btime $f_pd($dx, $x_pd, $es, $ed, nothing, $t)
# it seems like the `IONode` is even a bit faster.


text:
using Test

# Category
##########

A, B = Ob(FreeCategory, :A), Ob(FreeCategory, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)

# Expression types
@test isa(A, FreeCategory.Ob) && isa(f, FreeCategory.Hom)
@test isa(A, ObExpr) && isa(f, HomExpr)
@test isa(A, CategoryExpr) && isa(f, CategoryExpr)

# Domains and codomains
@test dom(f) == A
@test codom(f) == B
@test dom(compose(f,g)) == A
@test codom(compose(f,g)) == A
@test_throws SyntaxDomainError compose(f,f)

# Associativity
@test compose(compose(f,g),f) == compose(f,compose(g,f))

# Extra syntax
@test compose(f,g,f) == compose(compose(f,g),f)
@test compose([f,g,f]) == compose(compose(f,g),f)
@test f⋅g == compose(f,g)
@test f⋅g⋅f⋅g == compose(f,g,f,g)

# String format
@test string(A) == "A"
@test string(f) == "f"
@test string(compose(f,g)) == "compose(f,g)"
@test string(compose(f,g,f)) == "compose(f,g,f)"
@test string(Ob(FreeCategory, nothing)) != ""

# S-expressions
@test sexpr(A) == ":A"
@test sexpr(f) == ":f"
@test sexpr(compose(f,g)) == "(compose :f :g)"
@test sexpr(compose(f,g,f)) == "(compose :f :g :f)"

# Infix notation (Unicode)
@test unicode(A) == "A"
@test unicode(A, all=true) == "A"
@test unicode(f) == "f"
@test unicode(f, all=true) == "f: A → B"
@test unicode(id(A)) == "id{A}"
@test unicode(compose(f,g)) == "f⋅g"

# Infix notation (LaTeX)
@test latex(A) == "A"
@test latex(A, all=true) == raw"$A$"
@test latex(f) == "f"
@test latex(f, all=true) == raw"$f : A \to B$"
@test latex(id(A)) == "\\mathrm{id}_{A}"
@test latex(compose(f,g)) == "f \\cdot g"

@test latex(Ob(FreeCategory, "x")) == "x"
@test latex(Ob(FreeCategory, "sin")) == "\\mathrm{sin}"
@test latex(Ob(FreeCategory, "\\alpha")) == "\\alpha"

# (Co)presheaf
##############

A, B = Ob(FreePresheaf, :A, :B)
f = Hom(:f, A, B)
x = El(:x, B)
@test ob(x) == B
@test ob(coact(f, x)) == A
@test f ⋅ x == coact(f, x)

A, B = Ob(FreeCopresheaf, :A, :B)
f = Hom(:f, A, B)
x = El(:x, A)
@test ob(x) == A
@test ob(act(x, f)) == B
@test x ⋅ f == act(x, f)

# Infix notation (Unicode)
@test unicode(x) == "x"
@test unicode(x, all=true) == "x: A"
@test unicode(act(x, f)) == "x⋅f"
@test unicode(act(x, f), all=true) == "x⋅f: B"

# Infix notation (LaTeX)
@test latex(x) == "x"
@test latex(x, all=true) == raw"$x : A$"
@test latex(act(x, f)) == raw"x \cdot f"
@test latex(act(x, f), all=true) == raw"$x \cdot f : B$"


text:
@symbol_func function perimeter(cur_reactor::AbstractReactor)
  cur_perimeter = 2.0

  cur_perimeter *= cur_reactor.pi

  cur_perimeter *= a(cur_reactor)

  cur_perimeter *= sqrt(
    g_p(cur_reactor)
  )

  cur_perimeter
end


text:
using Pkg
Pkg.activate(".")
Pkg.add(url="https://github.com/AlgebraicJulia/CombinatorialSpaces.jl")
Pkg.add(url="https://github.com/bosonbaas/Catlab.jl", rev="rem_fix")
Pkg.develop(url="../..")
Pkg.instantiate()


text:
"""
    CoverPoly(p1::Pillar, p2::Pillar)

A polygon defined by the two pillars `p1` and `p2` where `p1` and `p2`
point to points ``x ∈ E_α(f)``.
"""
struct CoverPoly
    p1::Pillar
    p2::Pillar
end

struct CoveringPolygon
    p1_id::Int
    p2_id::Int
    order::SVector{2, Int}
    left_id::Int
    right_id::Int
end

mutable struct Covering
    pillars::Vector{Pillar}
    polygons::Dict{Int, CoveringPolygon}
    id_counter::Int # id of latest inserted polygon, increment this to get unique ids
end
function Covering(pillars::Vector{Pillar}, polygons::Vector{CoveringPolygon})
    id_counter = length(polygons)
    Covering(pillars, Dict(enumerate(polygons)), id_counter)
end

Covering(S::Spine2D; kwargs...) = Covering(S, S.amoeba_approximation.grid; kwargs...)
Covering(S::Spine2D, grid::Grid2D; kwargs...) = Covering(S, limits(grid)...; kwargs...)
function Covering(S::Spine2D, xmin::Real, xmax::Real, ymin::Real, ymax::Real)
    handle_halfray_intersections=true
    spine = S.curve

    # we have to ensue that the points of the complement complements are inside the
    # window. So we enlarge it if necessary
    CCs = components_complement(S)
    lims = covering_limits(CCs, xmin, xmax, ymin, ymax)
    intersect_vertices, vertices, linesegments = vertices_linesegments(S.curve, lims)

    order_pillars_dict = make_order_pillars_dict(CCs, vertices, linesegments)
    sorted_order_pillars = sort_order_pillars(order_pillars_dict, CCs, intersect_vertices)

    # We need to take care of the extra covering at the intersections
    # due to our sorting these are origins of the first and the last pillar
    if handle_halfray_intersections
        for (cc, pillars) in zip(CCs, sorted_order_pillars)
            if cc.bounded
                continue
            end
            p_start = intersection_cover(pillars[1], pillars[2], xmin, xmax, ymin, ymax)
            pushfirst!(pillars, p_start)
            p_end = intersection_cover(pillars[end], pillars[end-1], xmin, xmax, ymin, ymax)
            push!(pillars, p_end)
        end
    end
    # sorted_order_pillars is now a vector where each entry is a vector containing
    # all pillars to a given order (sorted in positive direction).
    # We can now construct the polygonal cover

    k = 0
    pillar_vec = Vector{Pillar}()
    polygons = Vector{CoveringPolygon}()
    poly_id = 0
    pillar_id = 0
    for (cc, pillars) in zip(CCs, sorted_order_pillars)
        n = length(pillars)
        if cc.bounded
            j = n
            i = 1
            while i ≤ n
                p1, p2 = pillar_id + j, pillar_id + i
                left_id = poly_id + j
                if i == n
                    right_id = poly_id + 1
                else
                    right_id = poly_id + i + 1
                end
                poly = CoveringPolygon(p1, p2, cc.order, left_id, right_id)
                push!(polygons, poly)
                j = i
                i += 1
            end
            poly_id += n
        else
            for i=1:n-1
                p1, p2 = pillar_id + i, pillar_id + i + 1
                left_id = i == 1 ? 0 : poly_id + i - 1
                right_id = i == n - 1 ? 0 : poly_id + i + 1
                poly = CoveringPolygon(p1, p2, cc.order, left_id, right_id)
                push!(polygons, poly)
            end
            poly_id += n - 1
        end
        pillar_id += n
        append!(pillar_vec, pillars)
    end
    Covering(pillar_vec, polygons)
end

function covering_limits(CCs, xmin, xmax, ymin, ymax)
    x_complement_points = first.(point.(CCs))
    y_complement_points = last.(point.(CCs))
    xmin = min(xmin, minimum(x_complement_points))
    xmax = max(xmax, maximum(x_complement_points))
    ymin = min(ymin, minimum(y_complement_points))
    ymax = max(ymax, maximum(y_complement_points))

    return xmin, xmax, ymin, ymax
end

function vertices_linesegments(spine, limits)
    linesegments = Vector{SVector{2, Point}}()
    all_vertices = Vector{SVector{2,Float64}}()
    xmin, xmax, ymin, ymax = limits
    for (origin, direction) in halfrays(spine)
        intersec, _, _ = intersection_halfray_limits(origin, direction, xmin, xmax, ymin, ymax)
        push!(linesegments, SVector(origin, intersec))
        push!(all_vertices, intersec)
    end

    intersect_vertices = copy(all_vertices)
    append!(linesegments, segments(spine))
    append!(all_vertices, vertices(spine))

    intersect_vertices, all_vertices, linesegments
end
function make_order_pillars_dict(CCs, vertices, linesegments)
    order_pillars_dict = map(cc -> (cc.order, Vector{Pillar}()), CCs) |> Dict
    for a in vertices
        for cc in CCs
            p = cc.x
            Δap = p - a
            nointersection = all(linesegments) do l
                !linesegments_intersect_pp_pd(l, (a, Δap))
            end
            # we found a proper segment
            if nointersection
                push!(order_pillars_dict[cc.order], (a, Δap))
            end
        end
    end
    return order_pillars_dict
end

function update_neighbor(p::CoveringPolygon, nb_id::Int, left_or_right::Symbol)
    if left_or_right == :left
        CoveringPolygon(p.p1_id, p.p2_id, p.order, nb_id, p.right_id)
    elseif left_or_right == :right
        CoveringPolygon(p.p1_id, p.p2_id, p.order, p.left_id, nb_id)
    else
        throw(error("Expected `:left` or `:right`, got `:$(left_or_right)`"))
    end
end

function cartesian_polygon(covering::Covering, cp::CoveringPolygon)
    cartesian_polygon(cp, covering.pillars)
end
function cartesian_polygon(cp::CoveringPolygon, pillars::Vector{Pillar})
    p1, p2 = pillars[cp.p1_id], pillars[cp.p2_id]
    (p2[1], sum(p2), sum(p1), p1[1])
end

area(covering::Covering, cp::CoveringPolygon) = area(cp, covering.pillars)
area(cp::CoveringPolygon, pillars) = polyarea(cartesian_polygon(cp, pillars))

polygons(covering::Covering) = covering.polygons
pillar(covering::Covering, id::Int) = covering.pillars[id]
polygon(covering::Covering, id::Int) = covering.polygons[id]

function _set!(covering, id::Int, pillar::Pillar)
    covering.pillars[id] = pillar
    covering
end
function _set!(covering, id::Int, poly::CoveringPolygon)
    covering.polygons[id] = poly
    covering
end

"""
    new_poly_id!(covering)

Create a new poly id.
"""
function new_poly_id!(covering::Covering)
    covering.id_counter += 1
    covering.id_counter
end
"""
    add_pillar!(covering, pillar)

Add a new pillar to the covering and return its id.
"""
function add_pillar!(covering::Covering, pillar::Pillar)
    push!(covering.pillars, pillar)
    length(covering.pillars)
end

"""
    add_polygon!(covering, id, poly)

Add a polygon with the given `id` to the covering.
"""
function add_polygon!(covering::Covering, id::Int, poly::CoveringPolygon)
    _set!(covering, id, poly)
    covering
end

"""
    remove_polygon!(covering, id, poly)

Remove a polygon with the given `id` from the covering.
"""
function remove_polygon!(covering::Covering, id::Int)
    delete!(covering.polygons, id)
    covering
end

"""
    update_polygon_neighbor!(covering, id, new_neighbor_id, left_right)

Update the left or right neighbor of the polygon with the given `id` to be
`new_neighbor_id`.
"""
function update_polygon_neighbor!(covering::Covering, id::Int, new_neighbor_id::Int, left_right::Symbol)
    poly = covering.polygons[id]
    covering.polygons[id] = update_neighbor(poly, new_neighbor_id, left_right)
    covering
end

function cartesian_polygons(covering::Covering)
    map(cp -> cartesian_polygon(cp, covering.pillars), values(covering.polygons))
end

function sort_order_pillars(order_pillars_dict, CCs, intersect_vertices)
    map(CCs) do cc
        pillars = order_pillars_dict[cc.order]
        vertices = first.(pillars)

        hull_indices = convexhull_vertices(vertices)
        hull = vertices[hull_indices]
        dirs = (last.(pillars))[hull_indices]
        sorted_pillars = collect(zip(hull, dirs))
        if cc.bounded
            return sorted_pillars
        else
            n = length(hull)
            if any(v -> v == hull[n], intersect_vertices) &&
               any(v -> v == hull[1], intersect_vertices)
               return sorted_pillars
            end

            for i=1:n-1
                if any(v -> v == hull[i], intersect_vertices) &&
                   any(v -> v == hull[i+1], intersect_vertices)
                   return [sorted_pillars[i+1:end]; sorted_pillars[1:i]]
                end
            end
        end
    end
end

function intersection_cover(intersection_pillar, node_pillar, xmin, xmax, ymin, ymax)
    intersection, Δintersection_p = intersection_pillar
    node, Δnode_p = node_pillar
    direction = intersection - node
    # we recompute the intersection to get the canidates for the pillar direciton
    _, d1, d2 = intersection_halfray_limits(node, 2direction, xmin, xmax, ymin, ymax)

    # we need to have a proper line segment to apply `pillars_onsameside`
    if pillars_onsameside((node, Δnode_p), (intersection, d1))
        dir = d1
    else
        dir = d2
    end
    # d is the projection of `dir` onto the line parallel to `direction`
    # through `p`
    d = (Δintersection_p × direction) / (dir × direction) * dir
    # # we maybe have to scale d, this can happen for
    # λ = scale_to_fit(intersection, d, xmin, xmax, ymin, ymax)
    (intersection, d)
end

"""
    scalepillars!(f, covering::Covering)

Scale each pillar of `covering` by a factor `t` returned by the function `f`.
`f` has as an input the origin and direction of the pillar.
"""
function scalepillars!(f, covering::Covering)
    N = length(covering.pillars)
    for k=1:N
        origin, dir = pillar(covering, k)
        t = f(origin, dir)
        covering.pillars[k] = (origin, t * dir)
    end
    covering
end

function split!(covering::Covering, polygon_id::Int, new_pillar::Pillar)
    pillar_id = add_pillar!(covering, new_pillar)
    poly = polygon(covering, polygon_id)

    new_poly_id_left = new_poly_id!(covering)
    new_poly_id_right = new_poly_id!(covering)

    # create new polygons
    left = CoveringPolygon(
        poly.p1_id, pillar_id, poly.order, poly.left_id, new_poly_id_right)
    right = CoveringPolygon(
        pillar_id, poly.p2_id, poly.order, new_poly_id_left, poly.right_id)

    # now we have to update the left and right neighbours of p
    if poly.left_id != 0
        update_polygon_neighbor!(covering, poly.left_id, new_poly_id_left, :right)
    end
    if poly.right_id != 0
        update_polygon_neighbor!(covering, poly.right_id, new_poly_id_right, :left)
    end

    remove_polygon!(covering, polygon_id)
    add_polygon!(covering, new_poly_id_left, left)
    add_polygon!(covering, new_poly_id_right, right)

    (new_poly_id_left, new_poly_id_right)
end

# AREA
function intersection_point(covering::Covering, poly::CoveringPolygon, ε)
    l1, l2 = approximation_lines(covering, poly, ε)
    lineintersection(l1, l2)
end
function intersection_point(c::Covering, poly_id::Int, ε)
    intersection_point(c, polygon(covering, poly_id), ε)
end

function approximation_lines(covering::Covering, poly::CoveringPolygon, ε)
    if poly.left_id != 0 && poly.right_id != 0
        left = polygon(covering, poly.left_id)
        right = polygon(covering, poly.right_id)

        # we have to account for the possible ε error
        a1, v1 = pillar(covering, poly.p1_id)
        a2, v2 = pillar(covering, poly.p2_id)

        b1 = a1 + (1 - ε / norm(v1)) * v1
        b2 = a2 + (1 - ε / norm(v2)) * v2

        l1 = (b1, sum(pillar(covering, left.p1_id)))
        l2 = (b2, sum(pillar(covering, right.p2_id)))

        return l1, l2
    elseif poly.left_id == 0
        # We have one of these special polygons
        right = polygon(covering, poly.right_id)

        # a1 == a2 and (a1 + t*v1) is part of the boundary
        a1, v1 = pillar(covering, poly.p1_id)
        a2, v2 = pillar(covering, poly.p2_id)

        b2 = a2 + (1 - ε / norm(v2)) * v2

        l1 = (a1, a1 + v1)
        l2 = (b2, sum(pillar(covering, right.p2_id)))

        return l1, l2
    else
        #poly.right_id == 0
        # We have one of these special polygons
        left = polygon(covering, poly.left_id)

        # a1 == a2 and (a2 + t*v2) is part of the boundary
        a1, v1 = pillar(covering, poly.p1_id)
        a2, v2 = pillar(covering, poly.p2_id)

        b1 = a1 + (1 - ε / norm(v1)) * v1

        l2 = (a2, a2 + v2)
        l1 = (b1, sum(pillar(covering, left.p1_id)))

        return l1, l2
    end
end

function approximated_local_error(covering::Covering, poly_id::Int, ε)
    approximated_local_error(covering, polygon(covering, poly_id), ε)
end
function approximated_local_error(covering::Covering, p::CoveringPolygon, ε)
    z = intersection_point(covering, p, ε)
    poly = cartesian_polygon(covering, p)
    if iscontained(poly, z)
        if p.left_id == 0
            x1 = sum(pillar(covering, p.p1_id))
            a2, v2 = pillar(covering, p.p2_id)
            y2 = a2 + (1 - ε / norm(v2)) * v2
            x2 = a2 + v2
            polyarea(z, y2, x2, x1)
        elseif p.right_id == 0
            a1, v1 = pillar(covering, p.p1_id)
            x1 = a1 + v1
            x2 = sum(pillar(covering, p.p2_id))

            y1 = a1 + (1 - ε / norm(v1)) * v1

            polyarea(z, x2, x1, y1)
        else
            a1, v1 = pillar(covering, p.p1_id)
            a2, v2 = pillar(covering, p.p2_id)
            y1 = a1 + (1 - ε / norm(v1)) * v1
            y2 = a2 + (1 - ε / norm(v2)) * v2
            x1, x2 = a1 + v1, a2 + v2
            return polyarea(z, y2, x2, x1, y1)
        end
    else
        return polyarea(poly)
    end
end

"""
    side(A::Point, B::Point, p::Point)

`side(A,B,p)>0` if `p` lies on one side of the line through `A` and `B` and `side(A,B,p)<0`
if `p` lies on te other side and `side(A,B,p) == 0` if `p` lies on the line.
"""
side(A, B, p) = (p[1] - A[1]) * (B[2] - A[2]) - (p[2] - A[2]) * (B[1] - A[1])

"""
    pillars_onsameside(p1::Pillar, p2::Pillar)

Check if two pillars point to the same side as defined by the line between the origins of
`p1` and `p2`.
"""
function pillars_onsameside(p1::Pillar, p2::Pillar)
     side(p1[1], p2[1], p1[1] + p1[2]) * side(p1[1], p2[1], p2[1] + p2[2]) >= 0
end

"""
    scale_to_fit(origin, d, xmin, xmax, ymin, ymax)

Return a scaling factor ``λ`` such that `origin + λ * d` is on the boundary of the window
determined by `xmin`, `xmax`, `ymin` and `ymax`.
"""
function scale_to_fit(origin::Point, d::Direction, xmin, xmax, ymin, ymax)
        x0, y0 = origin
        dx, dy = d

        if dx ≥ 0 && dy ≥ 0
            λx = (xmax - x0) / dx
            λy = (ymax - y0) / dy
        elseif dx < 0 && dy ≥ 0
            λx = (xmin - x0) / dx
            λy = (ymax - y0) / dy
        elseif dx ≥ 0 && dy < 0
            λx = (xmax - x0) / dx
            λy = (ymin - y0) / dy
        else
            λx = (xmin - x0) / dx
            λy = (ymin - y0) / dy
        end

        min(λx, λy)
end

"""
    linesegments_intersect(q::Point, s::Direction, p::Point, r::Direction)

Determines whether the line segments `(q, q + s)` and `(p, p + r)` intersect
"""
function linesegments_intersect(q::Point, s::Direction, p::Point, r::Direction)
    # algorithm taken from
    # https://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect
    unum = (q - p) × r
    tnum = (q - p) × s
    denom = r × s

    # colinear
    if unum ≈ 0.0 && denom ≈ 0
        norm2 = r ⋅ r
        t0 = (q - p) ⋅ r / norm2
        t1 = t0 + s ⋅ r / norm2
        if s ⋅ r < 0
            return 0 ≤ t0 ≤ 1 || 0 ≤ t1 ≤ 1 || t0 ≤ 0 && 1 ≤ t1
        else
            return 0 ≤ t0 ≤ 1 || 0 ≤ t1 ≤ 1 || t1 ≤ 0 && 1 ≤ t0
        end
    elseif denom ≈ 0
        return false
    else
        t = tnum / denom
        u = unum / denom

        # we only count proper intersections, not if endpoints coincide
        return 0 < t < 1 && 0 < u < 1
    end
end
@inline function linesegments_intersect(a::Tuple{Point, Direction}, b::Tuple{Point, Direction})
    linesegments_intersect(a[1], a[2], b[1], b[2])
end

function linesegments_intersect_pp_pd(a::SVector{2, Point}, b::Tuple{Point, Direction})
    linesegments_intersect(a[1], a[2] - a[1], b[1], b[2])
end

@recipe function plot(covering::Covering)
    fillalpha--> 0.4
    linecolor --> :dodgerblue
    fillcolor --> :dodgerblue
    linewidth --> 1.0
    legend --> false
    for (p1, p2, p3, p4) in cartesian_polygons(covering)
        @series begin
            st := :shape
            x1, y1 = p1
            x2, y2 = p2
            x3, y3 = p3
            x4, y4 = p4
            x5, y5 = p1

            [x1, x2, x3, x4, x5], [y1, y2, y3, y4, y5]
        end
    end
end


text:
using NeuralQuantum, Random, Plots

Random.seed!(Random.GLOBAL_RNG, 1234)

# Parameters of the transverse field ising model
N = 20
h = 1.0
J = 1.0

# Constructs the Hilbert space for N 1//2 spins.
hilb = HomogeneousSpin(N, 1//2)

# Builds the hamiltonian
H = LocalOperator(hilb)
for i=1:N
    global H  -= h * sigmax(hilb, i)
    global H  += J * sigmaz(hilb, i) * sigmaz(hilb, mod(i, N)+1)
end

net  = RBM(Float32, N, 1, af_logcosh)
init_random_pars!(net, sigma=0.01)

sampl = MetropolisSampler(LocalRule(), 125, N, burn=100)
algo  = SR(ϵ=(0.1), algorithm=sr_cg, precision=1e-3)

is = BatchedSampler(net, sampl, H, algo; batch_sz=8)

optimizer = Optimisers.Descent(0.1)

Evalues = Float64[];
Eerr = Float64[];
for i=1:300
    ldata, prec = sample!(is)
    ob = compute_observables(is)

    println("$i - $ldata")

    push!(Evalues, real(ldata.mean))
    push!(Eerr, ldata.error)
    grad = precondition!(prec, algo, i)
    Optimisers.update!(optimizer, net, grad)
end

plot(Evalues, yerr=Eerr)

# N=20 (thanks netket)
exact = -1.274549484318e00 * 20
hline!([exact, exact])


text:
const data = parse.(Int, split("219748365", ""))

function destination(cups::Vector{Int}, current::Int, n::Int=9)
    for i in 1:4
        a = mod1(current - i, n)
        a ∉ cups && return a
    end
end

function move!(allcups::Vector{Int}, current_ind::Int)
    current = allcups[current_ind]
    allcups = circshift(allcups, -current_ind)
    cups = splice!(allcups, 1:3)
    next_val = allcups[1]
    dest = findfirst(==(destination(cups, current)), allcups)
    allcups = prepend!(circshift(allcups, -dest), cups)
    return (allcups, findfirst(==(next_val), allcups))
end

function play(allcups::Vector{Int}, n::Int=100)
    dest = 1
    for i in 1:n
        i % 1000 == 0 && println("$(i / n) %")
        allcups, dest = move!(allcups, dest)
    end
    return allcups
end

function result(allcups::Vector{Int})
    one = findfirst(==(1), allcups)
    return join(circshift(allcups, -one)[1:8])
end

println("Part 1: $(result(play(data)))")

function tonextsdict(data::Vector{Int}, p1::Bool=false)
    nexts = Dict{Int,Int}()
    for (d, n) in zip(data, data[2:end])
        nexts[d] = n
    end
    if p1
        nexts[data[end]] = data[1]
        return nexts
    end
    nexts[data[end]] = 10
    for i in 10:999_999
        nexts[i] = i + 1
    end
    nexts[1_000_000] = data[1]
    return nexts
end

function move2!(nexts::Dict{Int,Int}, current::Int)
    n1 = nexts[current]
    n2 = nexts[n1]
    n3 = nexts[n2]
    nextcurrent = nexts[n3]
    dest = destination([n1, n2, n3], current, length(nexts))
    # insert here
    finalnext = nexts[dest]
    nexts[current] = nextcurrent
    nexts[dest] = n1
    nexts[n3] = finalnext
    return nexts, nextcurrent
end

function play2!(nexts::Dict{Int,Int}, start::Int, n::Int=10_000_000)
    current = start
    for i in 1:n
        nexts, current = move2!(nexts, current)
    end
    return nexts
end

function part2(data)
    nexts = tonextsdict(data)
    play2!(nexts, data[1])
    a = nexts[1]
    b = nexts[a]
    return a * b
end

println("Part 2: $(part2(data))")


text:
"""
An AbstractYield is an object which can be used as an argument to:

- zero-coupon spot rates via [`zero`](@ref)
- discount factor via [`discount`](@ref)
- accumulation factor via [`accumulation`](@ref)

It can be be constructed via:

- zero rate curve with [`Zero`](@ref)
- forward rate curve with [`Forward`](@ref)
- par rate curve with [`Par`](@ref)
- typical OIS curve with [`OIS`](@ref)
- typical constant maturity treasury (CMT) curve with [`CMT`](@ref)
"""
abstract type AbstractYield end

# make interest curve broadcastable so that you can broadcast over multiple`time`s in `interest_rate`
Base.Broadcast.broadcastable(ic::T) where {T<:AbstractYield} = Ref(ic)

struct YieldCurve{T,U,V} <: AbstractYield
    rates::T
    maturities::U
    zero::V # function time -> continuous zero rate
end


text:
# Example - BVAR
using Plots, CSV, StatsBase
include(joinpath(Pkg.dir("VectorAutoregressions"),"src","bvar.jl"))
plotly()

y = CSV.read(joinpath(Pkg.dir("VectorAutoregressions"),"test","bvar_data.csv"), header = false)
y = convert(Array,y[:,1:3])
const λ = 0.1:0.1:1
const τ = 10*(0.1:0.1:1)
const ε = 0.0001
const p = 4
const H = 10
const reps = 10000
const burnin = 8000
const max_try = 1000
const update = 1000
prior = Hyperparameter(λ,τ,ε,p,H,reps,burnin,max_try,update)

mForecast = fit_bvar(y,prior)
pFore = plot(layout = grid(1,3));
for j in 1:size(mForecast,3)
plot!(pFore,size(y,1)-H:size(y,1)+H, [y[size(y,1)-H:end,j]; median(mForecast[:,:,j],1)'], legend = false,subplot = j);
plot!(pFore,size(y,1)+1:size(y,1)+H, [[percentile(mForecast[:,i,j],16) for i in 1:H] median(mForecast[:,:,j],1)' [percentile(mForecast[:,i,j],84) for i in 1:H]], color = "red", line = [:dash :solid :dash], subplot = j)
end

gui(pFore)


text:
using Random
using LinearAlgebra
using Statistics
using Optim
using DataFrames
using DataFramesMeta
using CSV
using HTTP
using GLM

# read in function to create state transitions for dynamic model
include("create_grids.jl")

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 1: reshaping the data
#:::::::::::::::::::::::::::::::::::::::::::::::::::
# load in the data

url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdataBeta0.csv"
df = CSV.read(HTTP.get(url).body)

# create bus id variable
df = @transform(df, bus_id = 1:size(df,1))

#---------------------------------------------------
# reshape from wide to long (must do this twice be-
# cause DataFrames.stack() requires doing it one
# variable at a time)
#---------------------------------------------------
# first reshape the decision variable
dfy = @select(df, :bus_id,:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20,:RouteUsage,:Branded)
dfy_long = DataFrames.stack(dfy, Not([:bus_id,:RouteUsage,:Branded]))
rename!(dfy_long, :value => :Y)
dfy_long = @transform(dfy_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfy_long, Not(:variable))

# next reshape the odometer variable
dfx = @select(df, :bus_id,:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20)
dfx_long = DataFrames.stack(dfx, Not([:bus_id]))
rename!(dfx_long, :value => :Odometer)
dfx_long = @transform(dfx_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfx_long, Not(:variable))

# join reshaped df's back together
df_long = leftjoin(dfy_long, dfx_long, on = [:bus_id,:time])
sort!(df_long,[:bus_id,:time])

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 2: estimate a static version of the model
#:::::::::::::::::::::::::::::::::::::::::::::::::::

alpha_hat_glm = glm(@formula(Y ~ RouteUsage + Branded), df_long, Binomial(), LogitLink())
println(alpha_hat_glm)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 3a: read in data for dynamic model
#:::::::::::::::::::::::::::::::::::::::::::::::::::

url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdata.csv"
df = CSV.read(HTTP.get(url).body)

Y = Matrix(df[:,[:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20]])

Odo = Matrix(df[:,[:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20]])

X = Matrix(df[:,[:Xst1,:Xst2,:Xst3,:Xst4,:Xst5,:Xst6,:Xst7,:Xst8,:Xst9,:Xst10,:Xst11,:Xst12,:Xst13,:Xst14,:Xst15,:Xst16,:Xst17,:Xst18,:Xst19,:Xst20]])

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 3b: generate state transition matrices
#:::::::::::::::::::::::::::::::::::::::::::::::::::
zval,zbin,xval,xbin,xtran = create_grids()

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 3c: Compute the future value terms
#:::::::::::::::::::::::::::::::::::::::::::::::::::

f_v=[zeros(size(xtran, 1)), zeros(2), zeros(size(X, 2)+1)]


text:
# run gen_op(ohlc) before calling this macro
# note that Open() will be defined in REPL as TimeSeriesQueries.Open()

# @select ohlc where Open
macro select(s1::Symbol, s2::Symbol, s3::Symbol)
    Expr(:call, s3, esc(:($s1)))
end

# @select ohlc when monthname "October"
macro select(s1::Symbol, s2::Symbol, s3::Symbol, str::String)
    Expr(:call, :when, esc(:($s1)), s3, str)
end

# @select ohlc when year 2000
macro select(s1::Symbol, s2::Symbol, s3::Symbol, num::Int)
    Expr(:call, :when, esc(:($s1)), s3, num)
end

# @select ohlc where Open when monthname "January"
macro select(s1::Symbol, s2::Symbol, s3::Symbol, s4::Symbol, s5::Symbol,  str::String)
    chain(s3, Expr(:call, :when, esc(:($s1)), s5, str))
end

# @select ohlc where Open when year 2000
macro select(s1::Symbol, s2::Symbol, s3::Symbol, s4::Symbol, s5::Symbol,  num::Int)
    chain(s3, Expr(:call, :when, esc(:($s1)), s5, num))
end

# @select ohlc where Open .> 100 (or .> Close)
macro select(s1::Symbol, s2::Symbol, ex::Expr)
    # a check of comparison is another columne or an Int/Float
    typeof(eval(ex.args[3])) == Int || typeof(eval(ex.args[3])) == Float64 ?
    # limit of two scenarios
    chain(ex.args[1], Expr(:call, ex.args[2], esc(:($s1))), ex.args[3]) :
    chain(ex.args[1], Expr(:call, ex.args[2], esc(:($s1))), Expr(:call, ex.args[3], esc(:($s1))))
end

# @select ohlc where Open .> 100 (or .> Close) when dayname "Friday"
macro select(s1::Symbol, s2::Symbol, ex::Expr, s3::Symbol, s4::Symbol, str::String)
    # a check of comparison is another columne or an Int/Float
    typeof(eval(ex.args[3])) == Int || typeof(eval(ex.args[3])) == Float64 ?
    # limit of two scenarios
    chain(:when, Expr(:call, ex.args[1], Expr(:call, ex.args[2], esc(:($s1))), ex.args[3]), s4, str) :
    chain(:when, Expr(:call, ex.args[1], Expr(:call, ex.args[2], esc(:($s1))), Expr(:call, ex.args[3], esc(:($s1)))), s4, str)
end

# julia> @select ohlc where Open .> Close (or .> Close )when month 2
macro select(s1::Symbol, s2::Symbol, ex::Expr, s3::Symbol, s4::Symbol, num::Int)
    # a check of comparison is another columne or an Int/Float
    typeof(eval(ex.args[3])) == Int || typeof(eval(ex.args[3])) == Float64 ?
    # limit of two scenarios
    chain(:when, Expr(:call, ex.args[1], Expr(:call, ex.args[2], esc(:($s1))), ex.args[3]), s4, num) :
    chain(:when, Expr(:call, ex.args[1], Expr(:call, ex.args[2], esc(:($s1))), Expr(:call, ex.args[3], esc(:($s1)))), s4, num)
end


text:
a =
println(a)


text:
using Transit, Test
function run_test(path)
    println("Running $path")
    include(path)
end

# run_test("tset.jl")
run_test("decoder.jl")
run_test("encoder.jl")
run_test("roundtrip.jl")
# run_test("exemplar.jl")


text:
# -----------------------------------------------------------
#
# Using boundary terms and integration (summation) by parts
# to construct more accurate derivative matrices for time
# evolution using Lax-Friedricks-type fluxes
#
# All derivatives are computed here in pos basis first
# Then conjugated into hier basis
#
# -----------------------------------------------------------

# Efficiency criticality: MED
# Most computations time is taken up by functions called from
# other script files

# Accuracy criticality: LOW
# The accuracy is limited only by functions in other scripts
# and ODE.jl

# Construct the time evolution matrix given a laplacian and
# initial data given f0, v0
function wave_data{A<:AbstractArray, T<:Real}(laplac::A,
							f0coeffs::Array{T,1}, v0coeffs::Array{T,1})
	len = length(f0coeffs)
	rows = rowvals(laplac)
	vals = nonzeros(laplac)

	I = Int[]
	J = Int[]
	V = T[]
	for col = 1:len
		for i in nzrange(laplac, col)
			row = rows[i]
			val = vals[i]
			push!(I,row+len)
			push!(J,col)
			push!(V,val)
		end
	end
	for i = 1:len
		push!(I, i)
		push!(J, i+len)
		push!(V, 1.0)
	end
	RHS = sparse(I, J, V, 2*len, 2*len, +)
	y0 = Array{T}([i<=len ? f0coeffs[i] : v0coeffs[i-len] for i in 1:2*len])
	return RHS, y0
end

# The reason we have a wave evolution in 1D is to test
# that the standard "position" and multiresolution "heirarchical"
# bases agree for wave PDE evolution
function wave_evolve_1D(k::Int, max_level::Int,
							f0::Function, v0::Function,
							time0::Real, time1::Real; basis="hier", order = "45",
							kwargs...)
	if basis == "pos"
		f0coeffs = pos_vcoeffs_DG(k, max_level, f0)
		v0coeffs = pos_vcoeffs_DG(k, max_level, v0)
	elseif basis == "hier"
		f0coeffs = vcoeffs_DG(1, k, max_level, f0)
		v0coeffs = vcoeffs_DG(1, k, max_level, v0)
	elseif basis == "nodal"
		throw(MethodError("Nodal basis wave evolution not yet implemented"))
	elseif basis == "point"
		m2p = modal2points_1D(k, max_level)
		f0coeffs = m2p*vcoeffs_DG(1, k, max_level, f0)
		v0coeffs = m2p*vcoeffs_DG(1, k, max_level, v0)
	else
		throw(ArgumentError(:basis))
	end

	D_op = periodic_DLF_matrix(k, max_level; basis=basis)
	laplac= *(D_op,D_op)

	RHS, y0 = wave_data(laplac, f0coeffs, v0coeffs)
	if order == "45"
		soln = ode45((t,x)->*(RHS,x), y0, [time0,time1]; kwargs...)
	elseif order == "78"
		soln = ode78((t,x)->*(RHS,x), y0, [time0,time1]; kwargs...)
	else
		throw(ArgumentError(:order))
	end
	return soln
end

function norm_squared{T<:Real}(coeffs::Array{T}; basis="hier")

	if basis=="hier" || basis=="pos"
		return sum(i^2 for i in coeffs)
	else
		throw(MethodError("Nodal and point basis norms not yet implemented"))
	end
end

function energy_func_1D{T<:Real}(k, level,
						soln::Tuple{Array{T,1},Array{Array{T,1},1}};
						basis = "hier")
	len			= length(soln[1])
	num_coeffs	= Int(round(length(soln[2][1])/2))
	times		= copy(soln[1])
	energies	= Array{T}(len)
	D_op		= periodic_DLF_matrix(k, level; basis=basis)

	for i in 1:len
		ux   = *(D_op, soln[2][i][1:num_coeffs])
		udot = soln[2][i][num_coeffs+1:end]
		energies[i] = norm_squared(ux) + norm_squared(udot)

	end
	return (times, energies)
end

function wave_evolve(D::Int, k::Int, n::Int,
							  f0::Function, v0::Function,
							  time0::Real, time1::Real;
							  order="45", scheme="sparse", kwargs...)

	f0coeffs = vcoeffs_DG(D, k, n, f0; scheme=scheme)
	v0coeffs = vcoeffs_DG(D, k, n, v0; scheme=scheme)
	laplac   = laplacian_matrix(D, k, n; scheme=scheme)

	RHS, y0 = wave_data(laplac, f0coeffs, v0coeffs)
	if order == "45"
		soln = ode45((t,x)->*(RHS,x), y0, [time0,time1]; kwargs...)
	elseif order == "78"
		soln = ode78((t,x)->*(RHS,x), y0, [time0,time1]; kwargs...)
	else
		throw(ArgumentError(:order))
	end
	return soln
end

function energy_func{T<:Real}(D::Int, k::Int, n::Int,
					 soln::Tuple{Array{T,1},Array{Array{T,1},1}};
					 scheme = "sparse")

	len 		= length(soln[1])
	num_coeffs	= Int(round(length(soln[2][1])/2))
	times 		= copy(soln[1])
	energies	= Array{T}(len)
	D_ops		= grad_matrix(D, k, n; scheme=scheme)

	for i in 1:len
		ux   = [*(D_ops[j], soln[2][i][1:num_coeffs]) for j in 1:D]
		udot = soln[2][i][num_coeffs+1:end]
		energies[i] = sum([norm_squared(ux[j]) for j in 1:D])+norm_squared(udot)
	end
	return (times, energies)
end