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

function dispatch_variable_ref(model::InfiniteModel,
                               index::FiniteParameterIndex
                               )::FiniteParameterRef
    return FiniteParameterRef(model, index)
end

# Extend _add_data_object
function _add_data_object(model::InfiniteModel,
                          object::ScalarParameterData{<:IndependentParameter}
                          )::IndependentParameterIndex
    index =  MOIUC.add_item(model.independent_params, object)
    push!(model.param_object_indices, index)
    return index
end

function _add_data_object(model::InfiniteModel,
                          object::ScalarParameterData{<:FiniteParameter}
                          )::FiniteParameterIndex
    return MOIUC.add_item(model.finite_params, object)
end

# Extend _data_dictionary (type based)
function _data_dictionary(model::InfiniteModel,
    ::Type{IndependentParameter})::MOIUC.CleverDict
    return model.independent_params
end

function _data_dictionary(model::InfiniteModel,
    ::Type{FiniteParameter})::MOIUC.CleverDict
    return model.finite_params
end

# Extend _data_dictionary (ref based)
function _data_dictionary(pref::IndependentParameterRef)::MOIUC.CleverDict
    return JuMP.owner_model(pref).independent_params
end

function _data_dictionary(pref::FiniteParameterRef)::MOIUC.CleverDict
    return JuMP.owner_model(pref).finite_params
end

# Extend _data_object
function _data_object(pref::ScalarParameterRef)::AbstractDataObject
    object = get(_data_dictionary(pref), JuMP.index(pref), nothing)
    object === nothing && error("Invalid scalar parameter reference, cannot find " *
                           "corresponding parameter in the model. This is likely " *
                           "caused by using the reference of a deleted parameter.")
    return object
end

################################################################################
#                             CORE OBJECT METHODS
################################################################################
# Extend _core_variable_object for IndependentParameterRefs
function _core_variable_object(pref::IndependentParameterRef)::IndependentParameter
    return _data_object(pref).parameter
end

# Extend _core_variable_object for FiniteParameterRefs
function _core_variable_object(pref::FiniteParameterRef)::FiniteParameter
    return _data_object(pref).parameter
end

# Extend _parameter_number
function _parameter_number(pref::IndependentParameterRef)::Int
    return _data_object(pref).parameter_num
end

# Extend _parameter_numbers
function _parameter_numbers(pref::IndependentParameterRef)::Vector{Int}
    return [_parameter_number(pref)]
end

# Extend _object_number
function _object_number(pref::IndependentParameterRef)::Int
    return _data_object(pref).object_num
end

# Extend _object_numbers
function _object_numbers(pref::IndependentParameterRef)::Vector{Int}
    return [_object_number(pref)]
end

## Set helper methods for adapting data_objects with parametric changes
# No change needed
function _adaptive_data_update(pref::ScalarParameterRef, param::P,
                               data::ScalarParameterData{P})::Nothing where {P <: ScalarParameter}
    data.parameter = param
    return
end

# Reconstruction is necessary
function _adaptive_data_update(pref::ScalarParameterRef, param::P1,
                               data::ScalarParameterData{P2})::Nothing  where {P1, P2}
    new_data = ScalarParameterData(param, data.object_num, data.parameter_num,
                                   data.name, data.parameter_func_indices,
                                   data.infinite_var_indices,
                                   data.derivative_indices, data.measure_indices,
                                   data.constraint_indices, data.in_objective,
                                   data.generative_measures,
                                   data.has_internal_supports,
                                   data.has_generative_supports,
                                   data.has_deriv_constrs)
    _data_dictionary(pref)[JuMP.index(pref)] = new_data
    return
end

# Extend _set_core_variable_object for ScalarParameterRefs
function _set_core_variable_object(pref::ScalarParameterRef,
                                   param::ScalarParameter)::Nothing
    _adaptive_data_update(pref, param, _data_object(pref))
    return
end

################################################################################
#                           MACRO DEFINITION HELPERS
################################################################################
# Internal structure for building InfOptParameters
mutable struct _ParameterInfoExpr
    has_lb::Bool
    lower_bound::Any
    has_ub::Bool
    upper_bound::Any
    has_dist::Bool
    distribution::Any
    has_set::Bool
    set::Any
end

# Default constructor
function _ParameterInfoExpr(; lower_bound = NaN, upper_bound = NaN,
                            distribution = NaN, set = NaN)
    # isnan(::Expr) is not defined so we need to do !== NaN
    return _ParameterInfoExpr(lower_bound !== NaN, lower_bound,
                              upper_bound !== NaN, upper_bound,
                              distribution !== NaN, distribution,
                              set !== NaN, set)
end

# Internal function for use in processing valid key word arguments
function _is_set_keyword(kw::Expr)
    return kw.args[1] in [:set, :lower_bound, :upper_bound, :distribution]
end

# Extend to assist in building InfOptParameters
function _set_lower_bound_or_error(_error::Function,
                                   info::_ParameterInfoExpr, lower)::Nothing
    info.has_lb && _error("Cannot specify parameter lower_bound twice")
    info.has_dist && _error("Cannot specify parameter lower_bound and " *
                            "distribution")
    info.has_set && _error("Cannot specify parameter lower_bound and set")
    info.has_lb = true
    info.lower_bound = lower
    return
end

# Extend to assist in building InfOptParameters
function _set_upper_bound_or_error(_error::Function,
                                    info::_ParameterInfoExpr, upper)::Nothing
    info.has_ub && _error("Cannot specify parameter upper_bound twice")
    info.has_dist && _error("Cannot specify parameter upper_bound and " *
                            "distribution")
    info.has_set && _error("Cannot specify parameter upper_bound and set")
    info.has_ub = true
    info.upper_bound = upper
    return
end

# Extend to assist in building InfOptParameters
function _dist_or_error(_error::Function, info::_ParameterInfoExpr, dist)::Nothing
    info.has_dist && _error("Cannot specify parameter distribution twice")
    (info.has_lb || info.has_ub) && _error("Cannot specify parameter " *
                                           "distribution and upper/lower bounds")
    info.has_set && _error("Cannot specify parameter distribution and set")
    info.has_dist = true
    info.distribution = dist
    return
end

# Extend to assist in building InfOptParameters
function _set_or_error(_error::Function, info::_ParameterInfoExpr, set)::Nothing
    info.has_set && _error("Cannot specify variable fixed value twice")
    (info.has_lb || info.has_ub) && _error("Cannot specify parameter set and " *
                                           "upper/lower bounds")
    info.has_dist && _error("Cannot specify parameter set and distribution")
    info.has_set = true
    info.set = set
    return
end

# Construct an expression to build an infinite set (use with @infinite_macro)
function _constructor_set(_error::Function, info::_ParameterInfoExpr)
    if (info.has_lb || info.has_ub) && !(info.has_lb && info.has_ub)
        _error("Must specify both an upper bound and a lower bound")
    elseif info.has_lb
        check = :(isa($(info.lower_bound), Real))
        return :($(check) ? IntervalSet($(info.lower_bound), $(info.upper_bound)) : error("Bounds must be a number."))
    elseif info.has_dist
        check = :(isa($(info.distribution), Distributions.UnivariateDistribution))
        return :($(check) ? UniDistributionSet($(info.distribution)) : error("Distribution must be a Distributions.UnivariateDistribution."))
    elseif info.has_set
        check1 = :(isa($(info.set), InfiniteScalarSet))
        check2 = :(isa($(info.set), Distributions.UnivariateDistribution))
        return :($(check1) ? $(info.set) : ($(check2) ? UniDistributionSet($(info.set)) : error("Set must be a subtype of InfiniteScalarSet.")))
    else
        _error("Must specify upper/lower bounds, a distribution, or a set")
    end
end

################################################################################
#                            PARAMETER DEFINITION
################################################################################
# Define the default derivative evaluation method
const DefaultDerivativeMethod = FiniteDifference()

# Check that supports don't violate the set bounds
function _check_supports_in_bounds(_error::Function,
                                   supports::Union{<:Real, Vector{<:Real}},
                                   set::AbstractInfiniteSet)::Nothing
    if !supports_in_set(supports, set)
        _error("Supports violate the set domain bounds.")
    end
    return
end

"""
    build_parameter(_error::Function, set::InfiniteScalarSet;
                    [num_supports::Int = 0,
                    supports::Union{Real, Vector{<:Real}} = Real[],
                    sig_digits::Int = DefaultSigDigits,
                    derivative_method::AbstractDerivativeMethod = DefaultDerivativeMethod]
                    )::IndependentParameter

Returns a [`IndependentParameter`](@ref) given the appropriate information.
This is analagous to `JuMP.build_variable`. Errors if supports violate the
bounds associated with `set`. This is meant to primarily serve as a
helper method for [`@independent_parameter`](@ref). Here `derivative_method`
specifies the numerical evalution method that will be applied to derivatives that
are taken with respect to this infinite parameter.

**Example**
```julia-repl
julia> build_parameter(error, IntervalSet(0, 3), supports = Vector(0:3))
IndependentParameter{IntervalSet,FiniteDifference{Backward},NoGenerativeSupports}([0, 3], DataStructures.SortedDict(0.0 => Set([UserDefined]),1.0 => Set([UserDefined]),2.0 => Set([UserDefined]),3.0 => Set([UserDefined])), 12, FiniteDifference{Backward}(Backward(), true), NoGenerativeSupports())
```
"""
function build_parameter(_error::Function,
    set::InfiniteScalarSet;
    num_supports::Int = 0,
    supports::Union{Real, Vector{<:Real}} = Real[],
    sig_digits::Int = DefaultSigDigits,
    derivative_method::AbstractDerivativeMethod = DefaultDerivativeMethod,
    extra_kw_args...
    )::IndependentParameter
    for (kwarg, _) in extra_kw_args
        _error("Unrecognized keyword argument $kwarg")
    end
    label = UserDefined
    length_supports = length(supports)
    if !isempty(supports)
        supports = round.(supports, sigdigits = sig_digits)
        _check_supports_in_bounds(_error, supports, set)
        num_supports == 0 || @warn("Ignoring num_supports since supports is not empty.")
    elseif num_supports != 0
        supports, label = generate_support_values(set, num_supports = num_supports,
                                                  sig_digits = sig_digits)
    end
    supports_dict = DataStructures.SortedDict{Float64, Set{DataType}}(
                                            i => Set([label]) for i in supports)
    if length_supports != 0 && (length(supports_dict) != length_supports)
        @warn("Support points are not unique, eliminating redundant points.")
    end
    return IndependentParameter(set, supports_dict, sig_digits, derivative_method,
                                generative_support_info(derivative_method))
end

"""
    build_parameter(_error::Function, value::Real)::FiniteParameter

Returns a [`FiniteParameter`](@ref) given the appropriate information.
This is analagous to `JuMP.build_variable`. This is meant to primarily serve as
a helper method for [`@finite_parameter`](@ref).

**Example**
```jldoctest; setup = :(using InfiniteOpt)
julia> build_finite_parameter(error, 1)
FiniteParameter(1.0)
```
"""
function build_parameter(_error::Function, value::Real;
                         extra_kw_args...)::FiniteParameter
    for (kwarg, _) in extra_kw_args
        _error("Unrecognized keyword argument $kwarg")
    end
    return FiniteParameter(value)
end

"""
    add_parameter(model::InfiniteModel,
                  p::Union{IndependentParameterRef, FiniteParameterRef},
                  [name::String = ""])::GeneralVariableRef

Returns a [`GeneralVariableRef`](@ref) associated with the parameter `p` that is added
to `model`. This adds a parameter to the model in a manner similar to
`JuMP.add_variable`. This can be used to add parameters with the use of
[`@infinite_parameter`](@ref). [`build_parameter`](@ref) should be used to
construct `p`.

**Example**
```jldoctest; setup = :(using InfiniteOpt; model = InfiniteModel())
julia> p = build_independent_parameter(error, IntervalSet(0, 3), supports = Vector(0:3));

julia> param_ref = add_parameter(model, p, "name")
name
```
"""
function add_parameter(model::InfiniteModel, p::IndependentParameter,
                       name::String = "")::GeneralVariableRef
    obj_num = length(_param_object_indices(model)) + 1
    param_num = model.last_param_num += 1
    data_object = ScalarParameterData(p, obj_num, param_num, name)
    obj_index = _add_data_object(model, data_object)
    return GeneralVariableRef(model, obj_index.value, typeof(obj_index))
end

# FiniteParameter
function add_parameter(model::InfiniteModel, p::FiniteParameter,
                       name::String = "")::GeneralVariableRef
    data_object = ScalarParameterData(p, -1, -1, name)
    obj_index = _add_data_object(model, data_object)
    return GeneralVariableRef(model, obj_index.value, typeof(obj_index))
end

################################################################################
#                           PARAMETER DEPENDENCIES
################################################################################
# Extend _infinite_variable_dependencies
function _infinite_variable_dependencies(pref::ScalarParameterRef
    )::Vector{InfiniteVariableIndex}
    return _data_object(pref).infinite_var_indices
end

# Extend _parameter_function_dependencies
function _parameter_function_dependencies(pref::ScalarParameterRef
    )::Vector{ParameterFunctionIndex}
    return _data_object(pref).parameter_func_indices
end

# Extend _derivative_dependencies
function _derivative_dependencies(pref::ScalarParameterRef
    )::Vector{DerivativeIndex}
    return _data_object(pref).derivative_indices
end

# Extend _measure_dependencies
function _measure_dependencies(pref::ScalarParameterRef
    )::Vector{MeasureIndex}
    return _data_object(pref).measure_indices
end

# Extend _constraint_dependencies
function _constraint_dependencies(pref::ScalarParameterRef
    )::Vector{ConstraintIndex}
    return _data_object(pref).constraint_indices
end

# Extend _generative_measures
function _generative_measures(pref::ScalarParameterRef
    )::Vector{MeasureIndex}
    return _data_object(pref).generative_measures
end

################################################################################
#                             USED_BY FUNCTIONS
################################################################################
"""
    used_by_infinite_variable(pref::IndependentParameterRef)::Bool

Return true if `pref` is used by an infinite variable or false otherwise.

**Example**
```julia-repl
julia> used_by_infinite_variable(t)
true
```
"""
function used_by_infinite_variable(pref::IndependentParameterRef)::Bool
    return !isempty(_infinite_variable_dependencies(pref))
end

# FiniteParameter
used_by_infinite_variable(pref::FiniteParameterRef)::Bool = false

"""
    used_by_parameter_function(pref::IndependentParameterRef)::Bool

Return true if `pref` is used by an infinite parameter function or false otherwise.

**Example**
```julia-repl
julia> used_by_parameter_function(t)
false
```
"""
function used_by_parameter_function(pref::IndependentParameterRef)::Bool
    return !isempty(_parameter_function_dependencies(pref))
end

# FiniteParameter
used_by_parameter_function(pref::FiniteParameterRef)::Bool = false

"""
    used_by_measure(pref::Union{IndependentParameterRef, FiniteParameterRef})::Bool

Return true if `pref` is used by a measure or false otherwise.

**Example**
```julia-repl
julia> used_by_measure(t)
false
```
"""
function used_by_measure(pref::ScalarParameterRef)::Bool
    return !isempty(_measure_dependencies(pref))
end

"""
    used_by_constraint(pref::Union{IndependentParameterRef, FiniteParameterRef})::Bool

Return true if `pref` is used by a constraint or false otherwise.

**Example**
```julia-repl
julia> used_by_constraint(t)
true
```
"""
function used_by_constraint(pref::ScalarParameterRef)::Bool
    return !isempty(_constraint_dependencies(pref))
end

"""
    used_by_objective(pref::FiniteParameterRef)::Bool

Return true if `pref` is used by the objective function.

**Example**
```julia-repl
```
"""
function used_by_objective(pref::FiniteParameterRef)::Bool
    return _data_object(pref).in_objective
end

# IndependentParameter
used_by_objective(::IndependentParameterRef)::Bool = false

"""
    used_by_derivative(pref::IndependentParameterRef)::Bool

Return true if `pref` is used by a derivative or false otherwise.

**Example**
```julia-repl
julia> used_by_derivative(t)
false
```
"""
function used_by_derivative(pref::IndependentParameterRef)::Bool
    return !isempty(_derivative_dependencies(pref))
end

# FiniteParameter
used_by_derivative(::FiniteParameterRef)::Bool = false

"""
    is_used(pref::Union{IndependentParameterRef, FiniteParameterRef})::Bool

Return true if `pref` is used in the model or false otherwise.

**Example**
```julia-repl
julia> is_used(t)
true
```
"""
function is_used(pref::ScalarParameterRef)::Bool
    return used_by_measure(pref) || used_by_constraint(pref) ||
           used_by_infinite_variable(pref) || used_by_objective(pref) ||
           used_by_derivative(pref) || used_by_parameter_function(pref)
end

################################################################################
#                              NAME METHODS
################################################################################

"""
    JuMP.name(pref::Union{IndependentParameterRef, FiniteParameterRef})::String

Extend the [`JuMP.name`](@ref JuMP.name(::JuMP.VariableRef)) function to
accomodate infinite parameters. Returns the name string associated with `pref`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @independent_parameter(model, t in [0, 1]))
julia> name(t)
"t"
```
"""
function JuMP.name(pref::ScalarParameterRef)::String
    object = get(_data_dictionary(pref), JuMP.index(pref), nothing)
    return object === nothing ? "" : object.name
end

"""
    JuMP.set_name(pref::ScalarParameterRef, name::String)

Extend the [`JuMP.set_name`](@ref JuMP.set_name(::JuMP.VariableRef, ::String))
function to accomodate infinite parameters. Set a new base name to be associated
with `pref`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @independent_parameter(model, t in [0, 1]))
julia> set_name(t, "time")

julia> name(t)
"time"
```
"""
function JuMP.set_name(pref::ScalarParameterRef, name::String)::Nothing
    _data_object(pref).name = name
    JuMP.owner_model(pref).name_to_param = nothing
    return
end

# Make a parameter reference
function _make_parameter_ref(model::InfiniteModel,
                             index::AbstractInfOptIndex)::GeneralVariableRef
    return GeneralVariableRef(model, MOIUC.key_to_index(index), typeof(index))
end

function _make_parameter_ref(model::InfiniteModel,
                             index::DependentParameterIndex
                             )::GeneralVariableRef
    return GeneralVariableRef(model, MOIUC.key_to_index(index.object_index),
                              typeof(index), index.param_index)
end

# Get the name_to_param Dictionary
function _param_name_dict(model::InfiniteModel
    )::Union{Dict{String, AbstractInfOptIndex}, Nothing}
    return model.name_to_param
end

# Update name_to_param
function _update_param_name_dict(model::InfiniteModel,
    param_dict::MOIUC.CleverDict{K, V}
    )::Nothing where {K, V <: ScalarParameterData}
    name_dict = _param_name_dict(model)
    for (index, data_object) in param_dict
        param_name = data_object.name
        if haskey(name_dict, param_name)
            # IndependentParameterIndex(-1) is a special value that means
            # this string does not map to a unique variable name.
            name_dict[param_name] = IndependentParameterIndex(-1)
        else
            name_dict[param_name] = index
        end
    end
    model.name_to_param = name_dict
    return
end

function _update_param_name_dict(model::InfiniteModel,
    param_dict::MOIUC.CleverDict{K, V}
    )::Nothing where {K, V <: MultiParameterData}
    name_dict = _param_name_dict(model)
    for (index, data_object) in param_dict
        param_nums = data_object.parameter_nums
        for i in eachindex(param_nums)
            name = data_object.names[i]
            if haskey(name_dict, name)
                # IndependentParameterIndex(-1) is a special value that means
                # this string does not map to a unique variable name.
                name_dict[name] = IndependentParameterIndex(-1)
            else
                individual_index = DependentParameterIndex(index, i)
                name_dict[name] = individual_index
            end
         end
    end
    model.name_to_param = name_dict
    return
end

# TODO: this is not changing the name registered in model
"""
    parameter_by_name(model::InfiniteModel,
                      name::String)::Union{GeneralVariableRef, Nothing}

Return the parameter reference assoociated with a parameter name. Errors if
multiple parameters have the same name. Returns nothing if no such name exists.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @independent_parameter(model, t in [0, 1], supports = [0, 1]))
julia> parameter_by_name(model, "t")
t
```
"""
function parameter_by_name(model::InfiniteModel,
                           name::String)::Union{GeneralVariableRef, Nothing}
    if _param_name_dict(model) === nothing
        model.name_to_param = Dict{String, AbstractInfOptIndex}()
        _update_param_name_dict(model, model.independent_params)
        _update_param_name_dict(model, model.dependent_params)
        _update_param_name_dict(model, model.finite_params)
    end
    index = get(_param_name_dict(model), name, nothing)
    if index isa Nothing
        return nothing
    elseif index == IndependentParameterIndex(-1)
        error("Multiple variables have the name $name.")
    else
        return _make_parameter_ref(model, index)
    end
end

################################################################################
#                       GENERATIVE SUPPORT FUNCTIONS
################################################################################
# Extend copy for NoGenerativeSupports
function Base.copy(d::NoGenerativeSupports)::NoGenerativeSupports
    return NoGenerativeSupports()
end

# Extend copy for UniformGenerativeInfo
function Base.copy(d::UniformGenerativeInfo)::UniformGenerativeInfo
    return UniformGenerativeInfo(copy(d.support_basis), d.label)
end

"""
    support_label(info::AbstractGenerativeInfo)::DataType

Return the support label to be associated with generative supports produced in
accordance with `info`. This is intended an internal method that should be
extended for user defined types of [`AbstractGenerativeInfo`](@ref).
"""
function support_label(info::AbstractGenerativeInfo)
    error("`support_label` not defined for generative support info type " *
          "$(typeof(info)).")
end

# UniformGenerativeInfo
function support_label(info::UniformGenerativeInfo)::DataType
    return info.label
end

# NoGenerativeSupports
function support_label(info::NoGenerativeSupports)::DataType
    return _NoLabel
end

"""
    generative_support_info(pref::IndependentParameterRef)::AbstractGenerativeInfo

Return the generative support information associated with `pref`.
"""
function generative_support_info(pref::IndependentParameterRef)::AbstractGenerativeInfo
    return _core_variable_object(pref).generative_supp_info
end

"""
    has_generative_supports(pref::IndependentParameterRef)::Bool

Return whether generative supports have been added to `pref` in accordance
with its generative support info.
"""
function has_generative_supports(pref::IndependentParameterRef)::Bool
    return _data_object(pref).has_generative_supports
end

# Specify if a parameter has generative supports
function _set_has_generative_supports(pref::IndependentParameterRef,
                                      status::Bool)::Nothing
    _data_object(pref).has_generative_supports = status
    return
end

# Reset (remove) the generative supports if needed
function _reset_generative_supports(pref::IndependentParameterRef)::Nothing
    if has_generative_supports(pref)
        label = support_label(generative_support_info(pref))
        delete_supports(pref, label = label) # this also calls _set_has_generative_supports
    end
    return
end

# Specify the generative_support_info
function _set_generative_support_info(pref::IndependentParameterRef,
    info::AbstractGenerativeInfo)::Nothing
    sig_digits = significant_digits(pref)
    method = derivative_method(pref)
    set = _parameter_set(pref)
    supps = _parameter_supports(pref)
    new_param = IndependentParameter(set, supps, sig_digits, method, info)
    _reset_generative_supports(pref)
    _set_core_variable_object(pref, new_param)
    if is_used(pref)
        set_optimizer_model_ready(JuMP.owner_model(pref), false)
    end
    return
end

"""
    make_generative_supports(info::AbstractGenerativeInfo,
                             pref::IndependentParameterRef,
                             existing_supps::Vector{Float64}
                             )::Vector{Float64}

Generate the generative supports for `pref` in accordance with `info` and the
`existing_supps` that `pref` has. The returned supports should not include
`existing_supps`. This is intended as internal method to enable
[`add_generative_supports`](@ref) and should be extended for any user defined
`info` types that are created to enable new measure and/or derivative evaluation
techniques that require the creation of generative supports.
"""
function make_generative_supports(info::AbstractGenerativeInfo, pref, supps)
    error("`make_generative_supports` is not defined for generative support " *
          "info of type $(typeof(info)).")
end

# UniformGenerativeInfo
function make_generative_supports(info::UniformGenerativeInfo,
    pref, supps)::Vector{Float64}
    # collect the preliminaries
    basis = info.support_basis
    num_internal = length(basis)
    num_existing = length(supps)
    num_existing <= 1 && error("$(pref) does not have enough supports for " *
                                "creating generative supports.")
    internal_nodes = Vector{Float64}(undef, num_internal * (num_existing - 1))
    # generate the internal node supports
    for i in Iterators.take(eachindex(supps), num_existing - 1)
        lb = supps[i]
        ub = supps[i+1]
        internal_nodes[(i-1)*num_internal+1:i*num_internal] = basis * (ub - lb) .+ lb
    end
    return internal_nodes
end

## Define internal dispatch methods for adding generative supports
# AbstractGenerativeInfo
function _add_generative_supports(pref, info::AbstractGenerativeInfo)::Nothing
    if !has_generative_supports(pref)
        existing_supps = supports(pref, label = All)
        supps = make_generative_supports(info, pref, existing_supps)
        add_supports(pref, supps, label = support_label(info))
        _set_has_generative_supports(pref, true)
    end
    return
end

# NoGenerativeSupports
function _add_generative_supports(pref, info::NoGenerativeSupports)::Nothing
    return
end

"""
    add_generative_supports(pref::IndependentParameterRef)::Nothing

Create generative supports for `pref` if needed in accordance with its
generative support info using [`make_generative_supports`](@ref) and add them to
`pref`. This is intended as an internal function, but can be useful user defined
optimizer model extensions that utlize our support system.
"""
function add_generative_supports(pref::IndependentParameterRef)::Nothing
    info = generative_support_info(pref)
    _add_generative_supports(pref, info)
    return
end

################################################################################
#                        DERIVATIVE METHOD FUNCTIONS
################################################################################
# Determine if any derivatives have derivative constraints
function has_derivative_constraints(pref::IndependentParameterRef)::Bool
    return _data_object(pref).has_deriv_constrs
end

# Make update function for whether it has derivative supports
function _set_has_derivative_constraints(pref::IndependentParameterRef,
                                         status::Bool)::Nothing
    _data_object(pref).has_deriv_constrs = status
    return
end

"""
    derivative_method(pref::IndependentParameterRef)::AbstractDerivativeMethod

Returns the numerical derivative evaluation method employed with `pref` when it
is used as an operator parameter in a derivative.

**Example**
```julia-repl
julia> derivative_method(pref)
FiniteDifference(Backward, true)
```
"""
function derivative_method(pref::IndependentParameterRef)::AbstractDerivativeMethod
    return _core_variable_object(pref).derivative_method
end

# Make method to reset derivative constraints (supports are handled separately)
function _reset_derivative_constraints(pref::Union{IndependentParameterRef,
                                                   DependentParameterRef})::Nothing
    if has_derivative_constraints(pref)
        @warn("Support/method changes will invalidate existing derivative evaluation " *
              "constraints that have been added to the InfiniteModel. Thus, " *
              "these are being deleted.")
        for idx in _derivative_dependencies(pref)
            delete_derivative_constraints(DerivativeRef(JuMP.owner_model(pref), idx))
        end
        _set_has_derivative_constraints(pref, false)
    end
    return
end

"""
    set_derivative_method(pref::IndependentParameterRef,
                          method::AbstractDerivativeMethod)::Nothing

Specfies the desired derivative evaluation method `method` for derivatives that are
taken with respect to `pref`. Any internal supports exclusively associated with
the previous method will be deleted. Also, if any derivatives were evaluated
manually, the associated derivative evaluation constraints will be deleted. Errors
if new derivative method generates supports that are incompatible with existing
measures.

**Example**
```julia-repl
julia> set_derivative_method(d, OrthogonalCollocation(2))

```
"""
function set_derivative_method(pref::IndependentParameterRef,
    method::NonGenerativeDerivativeMethod
    )::Nothing
    old_param = _core_variable_object(pref)
    set = _parameter_set(pref)
    supps = _parameter_supports(pref)
    sig_figs = significant_digits(pref)
    if isempty(_generative_measures(pref))
        _reset_generative_supports(pref)
        new_param = IndependentParameter(set, supps, sig_figs, method,
                                         NoGenerativeSupports())
    else
        info = generative_support_info(pref)
        new_param = IndependentParameter(set, supps, sig_figs, method, info)
    end
    _reset_derivative_constraints(pref)
    _set_core_variable_object(pref, new_param)
    if is_used(pref)
        set_optimizer_model_ready(JuMP.owner_model(pref), false)
    end
    return
end

# GenerativeDerivativeMethod
function set_derivative_method(pref::IndependentParameterRef,
    method::GenerativeDerivativeMethod
    )::Nothing
    new_info = generative_support_info(method)
    old_info = generative_support_info(pref)
    if !isempty(_generative_measures(pref)) && new_info != old_info
        error("Generative derivative method conflicts with existing generative " *
              "measures.")
    end
    old_param = _core_variable_object(pref)
    set = _parameter_set(pref)
    supps = _parameter_supports(pref)
    sig_figs = significant_digits(pref)
    new_param = IndependentParameter(set, supps, sig_figs, method, new_info)
    _reset_derivative_constraints(pref)
    _reset_generative_supports(pref)
    _set_core_variable_object(pref, new_param)
    if is_used(pref)
        set_optimizer_model_ready(JuMP.owner_model(pref), false)
    end
    return
end

################################################################################
#                               SET FUNCTIONS
################################################################################
# Internal functions
function _parameter_set(pref::IndependentParameterRef)::InfiniteScalarSet
    return _core_variable_object(pref).set
end
function _update_parameter_set(pref::IndependentParameterRef,
                               set::AbstractInfiniteSet)::Nothing
    # old supports will always be discarded
    sig_digits = significant_digits(pref)
    method = derivative_method(pref)
    info = generative_support_info(pref)
    new_param = IndependentParameter(set, DataStructures.SortedDict{Float64, Set{DataType}}(),
                                     sig_digits, method, info)
    _set_core_variable_object(pref, new_param)
    _reset_derivative_constraints(pref)
    _set_has_generative_supports(pref, false)
    _set_has_internal_supports(pref, false)
    if is_used(pref)
        set_optimizer_model_ready(JuMP.owner_model(pref), false)
    end
    return
end

"""
    infinite_set(pref::IndependentParameterRef)::InfiniteScalarSet

Return the infinite set associated with `pref`.

**Example**
```julia-repl
julia> infinite_set(t)
[0, 1]
```
"""
function infinite_set(pref::IndependentParameterRef)::InfiniteScalarSet
    return _parameter_set(pref)
end

"""
    set_infinite_set(pref::IndependentParameterRef,
                     set::InfiniteScalarSet)::Nothing

Reset the infinite set of `pref` with another `InfiniteScalarSet`. An error will
be thrown if `pref` is being used by some measure.

**Example**
```julia-repl
julia> set_infinite_set(t, IntervalSet(0, 2))

julia> infinite_set(t)
[0, 2]
```
"""
function set_infinite_set(pref::IndependentParameterRef,
                          set::InfiniteScalarSet)::Nothing
    if used_by_measure(pref)
        error("$pref is used by a measure so changing its " *
              "infinite set is not allowed.")
    end
    _update_parameter_set(pref, set)
    return
end

################################################################################
#                        LOWER/UPPER BOUND FUNCTIONS
################################################################################
"""
    JuMP.has_lower_bound(pref::IndependentParameterRef)::Bool

Extend the `JuMP.has_lower_bound` function to accomodate infinite parameters.
Return true if the set associated with `pref` has a defined lower bound or if a
lower bound can be found. Extensions with user-defined infinite set types
should extend `JuMP.has_lower_bound(set::NewType)`.

**Example**
```julia-repl
julia> has_lower_bound(t)
true
```
"""
function JuMP.has_lower_bound(pref::IndependentParameterRef)::Bool
    set = _parameter_set(pref)
    return JuMP.has_lower_bound(set)
end

"""
    JuMP.lower_bound(pref::IndependentParameterRef)::Real

Extend the `JuMP.lower_bound` function to accomodate infinite parameters.
Returns the lower bound associated with the infinite set. Errors if such a bound
is not well-defined.

**Example**
```julia-repl
julia> lower_bound(t)
0.0
```
"""
function JuMP.lower_bound(pref::IndependentParameterRef)::Real
    set = _parameter_set(pref)
    if !JuMP.has_lower_bound(pref)
        error("Parameter $(pref) does not have a lower bound.")
    end
    return JuMP.lower_bound(set)
end

"""
    JuMP.set_lower_bound(pref::IndependentParameterRef, lower::Real)::Nothing

Extend the `JuMP.set_lower_bound` function to accomodate infinite parameters.
Updates the infinite set lower bound if such an operation is supported. Set
extensions that seek to employ this should extend
`JuMP.set_lower_bound(set::NewType, lower::Number)`.

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

julia> lower_bound(t)
-1.0
```
"""
function JuMP.set_lower_bound(pref::IndependentParameterRef, lower::Real)::Nothing
    set = _parameter_set(pref)
    new_set = JuMP.set_lower_bound(set, lower)
    _update_parameter_set(pref, new_set)
    return
end

"""
    JuMP.has_upper_bound(pref::IndependentParameterRef)::Bool

Extend the `JuMP.has_upper_bound` function to accomodate infinite parameters.
Return true if the set associated with `pref` has a defined upper bound or if a
upper bound can be found. Extensions with user-defined sets should extend
`JuMP.has_upper_bound(set::NewType)`.

**Example**
```julia-repl
julia> has_upper_bound(t)
true
```
"""
function JuMP.has_upper_bound(pref::IndependentParameterRef)::Bool
    set = _parameter_set(pref)
    return JuMP.has_upper_bound(set)
end

"""
    JuMP.upper_bound(pref::IndependentParameterRef)::Real

Extend the `JuMP.upper_bound` function to accomodate infinite parameters.
Returns the upper bound associated with the infinite set. Errors if such a bound
is not well-defined. Extensions with user-defined set types should extend
`JuMP.has_upper_bound(set::NewType)` and `JuMP.upper_bound(set::NewType)` if
appropriate.

**Example**
```julia-repl
julia> upper_bound(t)
1.0
```
"""
function JuMP.upper_bound(pref::IndependentParameterRef)::Real
    set = _parameter_set(pref)
    if !JuMP.has_upper_bound(pref)
        error("Parameter $(pref) does not have a upper bound.")
    end
    return JuMP.upper_bound(set)
end

"""
    JuMP.set_upper_bound(pref::IndependentParameterRef, lower::Real)::Nothing

Extend the `JuMP.set_upper_bound` function to accomodate infinite parameters.
Updates the infinite set upper bound if and only if it is an IntervalSet. Errors
otherwise. Extensions with user-defined infinite sets should extend
`JuMP.set_upper_bound(set::NewType, upper::Number)` if appropriate.

**Example**
```julia-repl
julia> set_upper_bound(t, 2)

julia> upper_bound(t)
2.0
```
"""
function JuMP.set_upper_bound(pref::IndependentParameterRef, upper::Real)::Nothing
    set = _parameter_set(pref)
    new_set = JuMP.set_upper_bound(set, upper)
    _update_parameter_set(pref, new_set)
    return
end

################################################################################
#                               SUPPORT FUNCTIONS
################################################################################
# Internal functions
function _parameter_supports(pref::IndependentParameterRef)
    return _core_variable_object(pref).supports
end
function _parameter_support_values(pref::IndependentParameterRef)::Vector{Float64}
    return collect(keys(_parameter_supports(pref)))
end
function _update_parameter_supports(pref::IndependentParameterRef,
    supports::DataStructures.SortedDict{Float64, Set{DataType}})::Nothing
    set = _parameter_set(pref)
    method = derivative_method(pref)
    sig_figs = significant_digits(pref)
    info = generative_support_info(pref)
    new_param = IndependentParameter(set, supports, sig_figs, method, info)
    _set_core_variable_object(pref, new_param)
    _reset_derivative_constraints(pref)
    _set_has_generative_supports(pref, false)
    if is_used(pref)
        set_optimizer_model_ready(JuMP.owner_model(pref), false)
    end
    return
end

"""
    has_internal_supports(pref::Union{IndependentParameterRef, DependentParameterRef})::Bool

Indicate if `pref` has internal supports that will be hidden from the user by
default.
"""
function has_internal_supports(
    pref::Union{IndependentParameterRef, DependentParameterRef}
    )::Bool
    return _data_object(pref).has_internal_supports
end

# update has internal supports
function _set_has_internal_supports(
    pref::Union{IndependentParameterRef, DependentParameterRef},
    status::Bool
    )::Nothing
    _data_object(pref).has_internal_supports = status
    return
end

"""
    significant_digits(pref::IndependentParameterRef)::Int

Return the number of significant digits enforced on the supports of `pref`.

**Example**
```julia-repl
julia> significant_digits(t)
12
```
"""
function significant_digits(pref::IndependentParameterRef)::Int
    return _core_variable_object(pref).sig_digits
end

"""
    num_supports(pref::IndependentParameterRef;
                 [label::Type{<:AbstractSupportLabel} = PublicLabel])::Int

Return the number of support points associated with `pref`. By default, only the
number of public supports are counted. The full amount can be determined by setting
`label = All`. Moreover, the amount of labels that satisfy `label` is obtained
using an [`AbstractSupportLabel`](@ref).

**Example**
```julia-repl
julia> num_supports(t)
2
```
"""
function num_supports(pref::IndependentParameterRef;
                      label::Type{<:AbstractSupportLabel} = PublicLabel)::Int
    supports_dict = _parameter_supports(pref)
    if label == All || (!has_internal_supports(pref) && label == PublicLabel)
        return length(supports_dict)
    else
        return count(p -> any(v -> v <: label, p[2]), supports_dict)
    end
end

"""
    has_supports(pref::IndependentParameterRef)::Bool

Return true if `pref` has supports or false otherwise.

**Example**
```julia-repl
julia> has_supports(t)
true
```
"""
has_supports(pref::IndependentParameterRef)::Bool = !isempty(_parameter_supports(pref))

"""
    supports(pref::IndependentParameterRef;
             [label::Type{<:AbstractSupportLabel} = PublicLabel])::Vector{Float64}

Return the support points associated with `pref`. Errors if there are no
supports. Users can query just support points generated by a certain method
using the keyword argument `label`. By default, the function returns all public
support points regardless of the associated label. The full collection is given by setting
`label = All`. Moreover, the amount of labels that satisfy `label` is obtained
using an [`AbstractSupportLabel`](@ref).

**Example**
```julia-repl
julia> supports(t)
2-element Array{Float64,1}:
 0.0
 1.0
```
"""
function supports(pref::IndependentParameterRef;
                  label::Type{<:AbstractSupportLabel} = PublicLabel)::Vector{Float64}
    if label == All || (!has_internal_supports(pref) && label == PublicLabel)
        return _parameter_support_values(pref)
    else
        return findall(x -> any(v -> v <: label, x), _parameter_supports(pref))
    end
end

# Return a matrix os supports when given a vector of IndependentParameterRefs (for measures)
function supports(prefs::Vector{IndependentParameterRef};
                  label::Type{<:AbstractSupportLabel} = PublicLabel,
                  use_combinatorics::Bool = true)::Matrix{Float64}
    # generate the support matrix considering all the unique combinations
    if use_combinatorics
        supp_list = Tuple(supports(p, label = label) for p in prefs)
        inds = CartesianIndices(ntuple(i -> 1:length(supp_list[i]), length(prefs)))
        supps = Matrix{Float64}(undef, length(prefs), length(inds))
        for (k, idx) in enumerate(inds)
            supps[:, k] = [supp_list[i][j] for (i, j) in enumerate(idx.I)]
        end
        return supps
    # generate the support matrix while negating the unique combinations
    else
        num_supps = num_supports(first(prefs), label = label)
        trans_supps = Matrix{Float64}(undef, num_supps, length(prefs))
        for i in eachindex(prefs)
            supp = supports(prefs[i], label = label)
            if length(supp) != num_supps
                error("Cannot simultaneously query the supports of multiple " *
                      "independent parameters if the support dimensions do not match " *
                      "while ignoring the combinatorics. Try setting `use_combinatorics = true`.")
            else
                @inbounds trans_supps[:, i] = supp
            end
        end
        return permutedims(trans_supps)
    end
end

"""
    set_supports(pref::IndependentParameterRef, supports::Vector{<:Real};
                 [force::Bool = false,
                 label::Type{<:AbstractSupportLabel} = UserDefined]
                 )::Nothing

Specify the support points for `pref`. Errors if the supports violate the bounds
associated with the infinite set. Warns if the points are not unique. If `force`
this will overwrite exisiting supports otherwise it will error if there are
existing supports.

**Example**
```julia-repl
julia> set_supports(t, [0, 1])

julia> supports(t)
2-element Array{Int64,1}:
 0
 1
```
"""
function set_supports(pref::IndependentParameterRef, supports::Vector{<:Real};
                      force::Bool = false,
                      label::Type{<:AbstractSupportLabel} = UserDefined
                      )::Nothing
    if has_supports(pref) && !force
        error("Unable set supports for $pref since it already has supports." *
              " Consider using `add_supports` or use set `force = true` to " *
              "overwrite the existing supports.")
    end
    set = _parameter_set(pref)
    supports = round.(supports, sigdigits = significant_digits(pref))
    _check_supports_in_bounds(error, supports, set)
    supports_dict = DataStructures.SortedDict{Float64, Set{DataType}}(
                                            i => Set([label]) for i in supports)
    if length(supports_dict) != length(supports)
        @warn("Support points are not unique, eliminating redundant points.")
    end
    _update_parameter_supports(pref, supports_dict)
    _set_has_internal_supports(pref, label <: InternalLabel)
    return
end

"""
    add_supports(pref::IndependentParameterRef,
                 supports::Union{Real, Vector{<:Real}};
                 [label::Type{<:AbstractSupportLabel} = UserDefined])::Nothing

Add additional support points for `pref` with identifying label `label`.

**Example**
```julia-repl
julia> add_supports(t, 0.5)

julia> supports(t)
3-element Array{Float64,1}:
 0.0
 0.5
 1.0

julia> add_supports(t, [0.25, 1])

julia> supports(t)
4-element Array{Float64,1}:
 0.0
 0.25
 0.5
 1.0
```
"""
function add_supports(pref::IndependentParameterRef,
                      supports::Union{Real, Vector{<:Real}};
                      label::Type{<:AbstractSupportLabel} = UserDefined,
                      check::Bool = true)::Nothing
    set = infinite_set(pref)
    supports = round.(supports, sigdigits = significant_digits(pref))
    check && _check_supports_in_bounds(error, supports, set)
    supports_dict = _parameter_supports(pref)
    for s in supports
        if haskey(supports_dict, s)
            push!(supports_dict[s], label)
        else
            supports_dict[s] = Set([label])
        end
    end
    _reset_derivative_constraints(pref)
    _reset_generative_supports(pref)
    if label <: InternalLabel
        _set_has_internal_supports(pref, true)
    end
    if is_used(pref)
        set_optimizer_model_ready(JuMP.owner_model(pref), false)
    end
    return
end

"""
    delete_supports(pref::IndependentParameterRef;
                    [label::Type{<:AbstractSupportLabel} = All])::Nothing

Delete the support points for `pref`. If `label != All` then delete `label` and
any supports that solely depend on it.

**Example**
```julia-repl
julia> delete_supports(t)

julia> supports(t)
ERROR: Parameter t does not have supports.
```
"""
function delete_supports(pref::IndependentParameterRef;
                         label::Type{<:AbstractSupportLabel} = All)::Nothing
    supp_dict = _parameter_supports(pref)
    if has_derivative_constraints(pref)
        @warn("Deleting supports invalidated derivative evaluations. Thus, these " *
              "are being deleted as well.")
        for idx in _derivative_dependencies(pref)
            delete_derivative_constraints(DerivativeRef(JuMP.owner_model(pref), idx))
        end
        _set_has_derivative_constraints(pref, false)
    end
    if label == All
        if used_by_measure(pref)
            error("Cannot delete the supports of $pref since it is used by " *
                  "a measure.")
        end
        empty!(supp_dict)
        _set_has_generative_supports(pref, false)
        _set_has_internal_supports(pref, false)
    else
        if has_generative_supports(pref) && support_label(generative_support_info(pref)) != label
            label = Union{label, support_label(generative_support_info(pref))}
        end
        _set_has_generative_supports(pref, false)
        filter!(p -> !all(v -> v <: label, p[2]), supp_dict)
        for (k, v) in supp_dict
            filter!(l -> !(l <: label), v)
        end
        if has_internal_supports(pref) && num_supports(pref, label = InternalLabel) == 0
            _set_has_internal_supports(pref, false)
        end
    end
    if is_used(pref)
        set_optimizer_model_ready(JuMP.owner_model(pref), false)
    end
    return
end

# Make dispatch for an array of parameters
function delete_supports(prefs::AbstractArray{<:IndependentParameterRef};
                         label::Type{<:AbstractSupportLabel} = All)::Nothing
    delete_supports.(prefs, label = label)
    return
end

"""
    parameter_value(pref::FiniteParameterRef)::Float64

Return the value of a finite parameter reference `pref`. Errors if it is
an infinite parameter.

**Example**
```julia-repl
julia> value(cost)
42.0
```
"""
function parameter_value(pref::FiniteParameterRef)::Real
    return _core_variable_object(pref).value
end

"""
    JuMP.set_value(pref::FiniteParameterRef, value::Real)::Nothing

Set the value of `pref` so long as it is a finite parameter. Errors if it is
an infinite parameter.

**Example**
```julia-repl
julia> set_value(cost, 27)

julia> value(cost)
27.0
```
"""
function JuMP.set_value(pref::FiniteParameterRef, value::Real)::Nothing
    _data_object(pref).parameter = FiniteParameter(value)
    if is_used(pref)
        set_optimizer_model_ready(JuMP.owner_model(pref), false)
    end
    return
end

"""
    fill_in_supports!(pref::IndependentParameterRef;
                      [num_supports::Int = DefaultNumSupports])::Nothing

Automatically generate support points for a particular independent parameter `pref`.
Generating `num_supports` for the parameter. The supports are generated uniformly
if the underlying infinite set is an `IntervalSet` or they are generating randomly
accordingly to the distribution if the set is a `UniDistributionSet`.
Will add nothing if there are supports
and `modify = false`. Extensions that use user defined set types should extend
[`generate_and_add_supports!`](@ref) and/or [`generate_support_values`](@ref)
as needed. Errors if the infinite set type is not recognized.

**Example**
```julia-repl
julia> fill_in_supports!(x, num_supports = 4)

julia> supports(x)
4-element Array{Number,1}:
 0.0
 0.333
 0.667
 1.0

```
"""
function fill_in_supports!(pref::IndependentParameterRef;
                           num_supports::Int = DefaultNumSupports,
                           modify::Bool = true)::Nothing
    set = infinite_set(pref)
    current_amount = length(_parameter_supports(pref))
    if (modify || current_amount == 0) && current_amount < num_supports
        generate_and_add_supports!(pref, set,
                                   num_supports = num_supports - current_amount,
                                   adding_extra = (current_amount > 0))
    end
    return
end

"""
    generate_and_add_supports!(pref::IndependentParameterRef,
                               set::AbstractInfiniteSet,
                               [method::Type{<:AbstractSupportLabel}];
                               [num_supports::Int = DefaultNumSupports])::Nothing

Generate supports for independent parameter `pref` via [`generate_support_values`](@ref)
and add them to `pref`. This is intended as an extendable internal method for
[`fill_in_supports!`](@ref fill_in_supports!(::IndependentParameterRef)).
Most extensions that empoy user-defined infinite sets can typically enable this
by extending [`generate_support_values`](@ref). Errors if the infinite set type
is not recognized.
"""
function generate_and_add_supports!(pref::IndependentParameterRef,
                                    set::AbstractInfiniteSet;
                                    num_supports::Int = DefaultNumSupports,
                                    adding_extra::Bool = false)::Nothing
    sig_digits = significant_digits(pref)
    if isa(set, IntervalSet) && adding_extra
        supports, label = generate_support_values(set, MCSample,
                                                  num_supports = num_supports,
                                                  sig_digits = sig_digits)
    else
        supports, label = generate_supports(set,
                                            num_supports = num_supports,
                                            sig_digits = sig_digits)
    end
    add_supports(pref, supports, label = label)
    return
end

# Dispatch with method
function generate_and_add_supports!(pref::IndependentParameterRef,
                                    set::AbstractInfiniteSet,
                                    method::Type{<:AbstractSupportLabel};
                                    num_supports::Int = DefaultNumSupports,
                                    adding_extra::Bool = false)::Nothing
    sig_digits = significant_digits(pref)
    supports, label = generate_supports(set, method,
                                        num_supports = num_supports,
                                        sig_digits = sig_digits)
    add_supports(pref, supports, label = label)
    return
end

################################################################################
#                               DELETE FUNCTIONS
################################################################################
# Check if parameter is used by measure data and error if it is to prevent bad
# deleting behavior
function _check_param_in_data(pref::GeneralVariableRef,
                              data::AbstractMeasureData)::Nothing
    prefs = parameter_refs(data)
    if (pref == prefs || pref in prefs)
        error("Unable to delete `$pref` since it is used to evaluate measures.")
    end
    return
end

# Update a reduced variable associated with an infinite variable whose parameter
# was removed
function _update_reduced_variable(vref::ReducedVariableRef,
                                  delete_indices::UnitRange{Int})::Nothing
    eval_supps = eval_supports(vref)
    new_supports = Dict{Int, Float64}()
    num_indices = length(delete_indices)
    for (index, support) in eval_supps
        if index < first(delete_indices)
            new_supports[index] = support
        elseif index > last(delete_indices)
            new_supports[index - num_indices] = support
        end
    end
    obj_nums = _object_numbers(vref)
    param_nums = _parameter_numbers(vref)
    new_var = ReducedVariable(infinite_variable_ref(vref), new_supports,
                              param_nums, obj_nums)
    _set_core_variable_object(vref, new_var)
    # TODO account for possibility that this becomes a point variable
    # TODO account for possibility that this becomes an infinite variable
    return
end

# Update the dependent measures
function _update_measures(model::InfiniteModel,
                          pref::GeneralVariableRef)::Nothing
    for mindex in _measure_dependencies(pref)
        mref = dispatch_variable_ref(model, mindex)
        func = measure_function(mref)
        if func isa GeneralVariableRef
            data = measure_data(mref)
            new_func = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
            new_meas = Measure(new_func, data, Int[], Int[], true)
            _set_core_variable_object(mref, new_meas)
        else
            _remove_variable(func, pref)
        end
    end
    return
end

# Update the dependent constraints
function _update_constraints(model::InfiniteModel,
                             pref::GeneralVariableRef)::Nothing
    for cindex in _constraint_dependencies(pref)
        cref = _temp_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))
        else
            _remove_variable(func, pref)
        end
    end
    return
end

# Remove given object/parameter number and update the list
function _update_number_list(nums::Vector{Int}, list::Vector{Int})::Nothing
    filter!(e -> !(e in nums), list)
    max_num = maximum(nums)
    for i in eachindex(list)
        if list[i] > max_num
            list[i] -= length(nums)
        end
    end
    return
end

# Update the model with the removed parameter/object numbers
function _update_model_numbers(model::InfiniteModel, obj_num::Int,
                               param_nums::Vector{Int})::Nothing
    # update the independent parameters
    for (_, object) in _data_dictionary(model, IndependentParameter)
        if object.object_num > obj_num
            object.object_num -= 1
            object.parameter_num -= length(param_nums)
        end
    end
    # update the dependent parameters
    for (_, object) in _data_dictionary(model, DependentParameters)
        if object.object_num > obj_num
            object.object_num -= 1
            object.parameter_nums = object.parameter_nums .- length(param_nums)
        end
    end
    # update the infinite parameter functions
    for (_, object) in model.param_functions
        _update_number_list([obj_num], object.func.object_nums)
        _update_number_list(param_nums, object.func.parameter_nums)
    end
    # update the infinite variables
    for vref in JuMP.all_variables(model, InfiniteVariable)
        _update_number_list([obj_num], _object_numbers(vref))
        _update_number_list(param_nums, _parameter_numbers(vref))
    end
    # update the reduced variables
    for vref in JuMP.all_variables(model, ReducedVariable)
        _update_number_list([obj_num], _object_numbers(vref))
        _update_number_list(param_nums, _parameter_numbers(vref))
    end
    # update the measures
    for mref in all_measures(model)
        _update_number_list([obj_num], _object_numbers(mref))
        _update_number_list(param_nums, _parameter_numbers(mref))
    end
    # update the constraints
    for (_, object) in model.constraints
        _update_number_list([obj_num], object.object_nums)
    end
    # update the central info
    deleteat!(_param_object_indices(model), obj_num)
    model.last_param_num -= length(param_nums)
    return
end

"""
    JuMP.delete(model::InfiniteModel, pref::ScalarParameterRef)::Nothing

Extend [`JuMP.delete`](@ref JuMP.delete(::JuMP.Model, ::JuMP.VariableRef)) to delete
scalar parameters and their dependencies. All variables, constraints, and
measure functions that depend on `pref` are updated to exclude it. Errors if the
parameter is contained in an `AbstractMeasureData` datatype that is employed by
a measure since the measure becomes invalid otherwise. Thus, measures that
contain this dependency must be deleted first. Note that
[`parameter_refs`](@ref parameter_refs(::AbstractMeasureData)) needs to be
extended to allow deletion of parameters when custom `AbstractMeasureData`
datatypes are used. Note that any dependent infinite variables will have their
start values reset via [`reset_start_value_function`](@ref).

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

julia> delete(model, x)

julia> print(model)
Min measure(g(t)*t) + z
Subject to
 g(t) + z ≥ 42.0, ∀ t ∈ [0, 6]
 g(0.5) = 0
```
"""
function JuMP.delete(model::InfiniteModel, pref::IndependentParameterRef)::Nothing
    @assert JuMP.is_valid(model, pref) "Parameter reference is invalid."
    gvref = _make_parameter_ref(JuMP.owner_model(pref), JuMP.index(pref))
    # ensure deletion is okay (pref isn't used by measure data)
    for mindex in _measure_dependencies(pref)
        data = measure_data(dispatch_variable_ref(model, mindex))
        _check_param_in_data(gvref, data)
    end
    # ensure pref is not used by a parameter function
    if used_by_parameter_function(pref)
        error("Cannot delete `$pref` since it is used by an infinite parameter " *
              "function.")
    end
    # update optimizer model status
    if is_used(pref)
        set_optimizer_model_ready(model, false)
    end
    # delete dependence of measures on pref
    _update_measures(model, gvref)
    # update infinite variables that depend on pref
    for vindex in _infinite_variable_dependencies(pref)
        # remove the parameter dependence
        vref = InfiniteVariableRef(model, vindex)
        prefs = raw_parameter_refs(vref)
        delete_index = findfirst(isequal(gvref), prefs)
        deleteat!(prefs, delete_index)
        reset_start_value_function(vref)
        # update any point variables that depend on vref accordingly
        for pindex in _point_variable_dependencies(vref)
            pvref = PointVariableRef(model, pindex)
            deleteat!(raw_parameter_values(pvref), delete_index)
        end
        # update any reduced variables that depend on vref accordingly
        for rindex in _reduced_variable_dependencies(vref)
            rvref = ReducedVariableRef(model, rindex)
            _update_reduced_variable(rvref, delete_index:delete_index)
        end
    end
    # delete any derivatives that use pref
    for index in _derivative_dependencies(pref)
        JuMP.delete(model, dispatch_variable_ref(model, index))
    end
    # update constraints in mapping to remove the parameter
    _update_constraints(model, gvref)
    # delete parameter information stored in model
    obj_num = _object_number(pref)
    param_nums = _parameter_numbers(pref)
    _delete_data_object(pref)
    # update the object numbers and parameter numbers
    _update_model_numbers(model, obj_num, param_nums)
    return
end

# FiniteParameterRef
function JuMP.delete(model::InfiniteModel, pref::FiniteParameterRef)::Nothing
    @assert JuMP.is_valid(model, pref) "Parameter reference is invalid."
    # update optimizer model status
    if is_used(pref)
        set_optimizer_model_ready(model, false)
    end
    gvref = _make_parameter_ref(model, JuMP.index(pref))
    # delete dependence of measures on pref
    _update_measures(model, gvref)
    # update constraints in mapping to remove the parameter
    _update_constraints(model, gvref)
    # update the objective if necessary
    if used_by_objective(pref)
        func = JuMP.objective_function(model)
        if func 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(func, gvref)
        end
    end
    # delete parameter information stored in model
    _delete_data_object(pref)
    return
end


text:
function answer1()
    h = 0
    v = 0
    open("input") do file
        for line in eachline(file)
            order = split(line)
            if order[1] == "forward"
                h += parse(Int, order[2])
            elseif order[1] == "up"
                v -= parse(Int, order[2])
            else
                v += parse(Int, order[2])
            end
        end
    end
    h, v, h*v
end
@show answer1()

function answer2()
    h = 0
    v = 0
    aim = 0
    open("input") do file
        for line in eachline(file)
            order = split(line)
            if order[1] == "down"
                aim += parse(Int, order[2])
            elseif order[1] == "up"
                aim -= parse(Int, order[2])
            else
                h += parse(Int, order[2])
                v += aim * parse(Int, order[2])
            end
        end
    end
    h, v, h*v
end
@show answer2()

ext: jl


text:
# Navigate an object through a series of obstacles, while minimizing the jerk of the trajectory.
# The domain is broken down into a polyhedral covering, and the time interval into subintervals.
# In each subinterval, we enforce the object is contained entirely within one of the polyhedral
# subregions, which enforces that we do not run into any obstacles.
#
# The trajectory is described by polynomials p_x and p_y, which describe the x and y position,
# respectively, as a function of time. This is a polynomial optimization problem, which we relax
# into a mixed-integer sum of squares problem. The problem is infeasible for polynomials of degree
# 4 or less, and feasible for degree 5.
#
# Example written and by Joey Huchette.

# source: https://github.com/JuliaOpt/Pajarito.jl/blob/master/examples/collisionavoidance.jl

using JuMP, PolyJuMP, SumOfSquares, MultivariatePolynomials
using Pajarito

# Specify Pajarito solver

mip_solver_drives = true
rel_gap = 1e-4

# using Cbc
# mip_solver = CbcSolver()

using CPLEX
mip_solver = CplexSolver(
    CPX_PARAM_SCRIND=(mip_solver_drives ? 1 : 0),
    CPX_PARAM_EPINT=1e-9,
    CPX_PARAM_EPRHS=1e-9,
    CPX_PARAM_EPGAP=(mip_solver_drives ? 1e-5 : 1e-9)
)

# using Mosek
# conic_solver = MosekSolver(LOG=0)
# conic_solver = MosekSolver(MSK_DPAR_INTPNT_CO_TOL_PFEAS=1e-6,MSK_DPAR_INTPNT_CO_TOL_DFEAS=1e-6,MSK_DPAR_INTPNT_CO_TOL_REL_GAP=1e-5,MSK_DPAR_INTPNT_TOL_INFEAS=1e-8,MSK_DPAR_INTPNT_CO_TOL_MU_RED=1e-6)

micp_solver = PajaritoSolver(
    mip_solver_drives=mip_solver_drives,
    log_level=3,
    rel_gap=rel_gap,
    mip_solver=mip_solver,
    # cont_solver=conic_solver,
    solve_relax=false,
    solve_subp=false,
    prim_cuts_only=true
)

# Create model

struct Box
    xl::Float64
    xu::Float64
    yl::Float64
    yu::Float64
end

boxes = Box[
    Box(0.0,1.0,0.0,0.3),
    Box(0.8,1.7,0.1,0.3),
    Box(1.4,1.9,0.2,0.4),
    Box(1.0,1.7,0.3,0.5),
    Box(0.5,1.4,0.4,0.6),
    Box(0.0,1.0,0.5,0.7),
    Box(0.2,1.0,0.6,0.8),
    Box(0.5,1.3,0.7,0.9),
    Box(1.0,2.0,0.7,1.0)
]

N = 8 # Number of trajectory pieces
d = 2 # dimension of space
r = 5 # dimension of polynomial trajectories
M = 2 # number of horizontal segments

domain = Box(0,M,0,1)

X₀ = Dict(:x=>0, :y=>0)
X₀′ = Dict(:x=>1, :y=>0)
X₀′′ = Dict(:x=>0, :y=>0)
X₁ = Dict(:x=>2, :y=>1)

T = linspace(0, 1, N+1)

Tmin = minimum(T)
Tmax = maximum(T)

model = SOSModel(solver=micp_solver)

@polyvar(t)
Z = monomials([t], 0:r)

@variable(model, H[1:N,boxes], Bin)

Mxl, Mxu, Myl, Myu = domain.xl, domain.xu, domain.yl, domain.yu
p = Dict()
for j in 1:N
    @constraint(model, sum(H[j,box] for box in boxes) == 1)

    p[(:x,j)] = @polyvariable(model, _, Z)
    p[(:y,j)] = @polyvariable(model, _, Z)
    for box in boxes
        xl, xu, yl, yu = box.xl, box.xu, box.yl, box.yu
        @polyconstraint(model, p[(:x,j)] >= Mxl + (xl-Mxl)*H[j,box], domain = (t >= T[j] && t <= T[j+1]))
        @polyconstraint(model, p[(:x,j)] <= Mxu + (xu-Mxu)*H[j,box], domain = (t >= T[j] && t <= T[j+1]))
        @polyconstraint(model, p[(:y,j)] >= Myl + (yl-Myl)*H[j,box], domain = (t >= T[j] && t <= T[j+1]))
        @polyconstraint(model, p[(:y,j)] <= Myu + (yu-Myu)*H[j,box], domain = (t >= T[j] && t <= T[j+1]))
    end
end

for axis in (:x,:y)
    @constraint(model, p[(axis,1)]([Tmin], [t]) == X₀[axis])
    @constraint(model, differentiate(p[(axis,1)], t)([Tmin], [t]) == X₀′[axis])
    @constraint(model, differentiate(p[(axis,1)], t, 2)([Tmin], [t]) == X₀′′[axis])

    for j in 1:N-1
        @constraint(model, p[(axis,j)]([T[j+1]], [t]) == p[(axis,j+1)]([T[j+1]], [t]))
        @constraint(model, differentiate(p[(axis,j)], t )([T[j+1]], [t]) == differentiate(p[(axis,j+1)], t)([T[j+1]], [t]))
        @constraint(model, differentiate(p[(axis,j)], t, 2)([T[j+1]], [t]) == differentiate(p[(axis,j+1)], t, 2)([T[j+1]], [t]))
    end

    @constraint(model, p[(axis,N)]([Tmax], [t]) == X₁[axis])
end

@variable(model, γ[keys(p)] ≥ 0)
for (key,val) in p
    @constraint(model, γ[key] ≥ norm(differentiate(val, t, 3)))
end
@objective(model, Min, sum(γ))

solve(model)

# Interpret solution: eval_poly(r) gives the trajectory location at time r

PP = Dict(key => getvalue(p[key]) for key in keys(p))
HH = getvalue(H)

function eval_poly(r)
    for i in 1:N
        if T[i] <= r <= T[i+1]
            return PP[(:x,i)]([r], [t]), PP[(:y,i)]([r], [t])
            break
        end
    end
    error("Time $r out of interval [$(minimum(T)),$(maximum(T))]")
end


text:
export ScaledInterpolation, eachvalue

import Base: done, next, start

immutable ScaledInterpolation{T,N,ITPT,IT,GT,RT} <: AbstractInterpolationWrapper{T,N,ITPT,IT,GT}
    itp::ITPT
    ranges::RT
end
@generated function ScaledInterpolation{ITPT,RT}(itp::ITPT, ranges::RT)
    T = eltype(itp)
    N = ndims(itp)
    IT = itptype(itp)
    GT = gridtype(itp)
    :(ScaledInterpolation{$T,$N,$ITPT,$IT,$GT,$RT}(itp, ranges))
end

"""
`scale(itp, xs, ys, ...)` scales an existing interpolation object to allow for indexing using other coordinate axes than unit ranges, by wrapping the interpolation object and transforming the indices from the provided axes onto unit ranges upon indexing.

The parameters `xs` etc must be either ranges or linspaces, and there must be one coordinate range/linspace for each dimension of the interpolation object.

For every `NoInterp` dimension of the interpolation object, the range must be exactly `1:size(itp, d)`.
"""
function scale{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, ranges::Range...)
    length(ranges) == N || throw(ArgumentError("Must scale $N-dimensional interpolation object with exactly $N ranges (you used $(length(ranges)))"))
    for d in 1:N
        if iextract(IT,d) != NoInterp
            length(ranges[d]) == size(itp,d) || throw(ArgumentError("The length of the range in dimension $d ($(length(ranges[d]))) did not equal the size of the interpolation object in that direction ($(size(itp,d)))"))
        elseif ranges[d] != 1:size(itp,d)
            throw(ArgumentError("NoInterp dimension $d must be scaled with unit range 1:$(size(itp,d))"))
        end
    end

    ScaledInterpolation(itp, ranges)
end

@generated function getindex{T,N,ITPT,IT<:DimSpec}(sitp::ScaledInterpolation{T,N,ITPT,IT}, xs::Number...)
    length(xs) == N || throw(ArgumentError("Must index into $N-dimensional scaled interpolation object with exactly $N indices (you used $(length(xs)))"))
    interp_indices = map(i -> iextract(IT, i) != NoInterp ? :(coordlookup(sitp.ranges[$i], xs[$i])) : :(xs[$i]), 1:N)
    return :($(Expr(:meta,:inline)); getindex(sitp.itp, $(interp_indices...)))
end

getindex{T}(sitp::ScaledInterpolation{T,1}, x::Number, y::Int) = y == 1 ? sitp[x] : throw(BoundsError())

size(sitp::ScaledInterpolation, d) = size(sitp.itp, d)
lbound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnGrid}, d) = 1 <= d <= N ? sitp.ranges[d][1] : throw(BoundsError())
lbound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnCell}, d) = 1 <= d <= N ? sitp.ranges[d][1] - boundstep(sitp.ranges[d]) : throw(BoundsError())
ubound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnGrid}, d) = 1 <= d <= N ? sitp.ranges[d][end] : throw(BoundsError())
ubound{T,N,ITPT,IT}(sitp::ScaledInterpolation{T,N,ITPT,IT,OnCell}, d) = 1 <= d <= N ? sitp.ranges[d][end] + boundstep(sitp.ranges[d]) : throw(BoundsError())

boundstep(r::LinSpace) = ((r.stop - r.start) / r.divisor) / 2
boundstep(r::FloatRange) = r.step / 2
boundstep(r::StepRange) = r.step / 2
boundstep(r::UnitRange) = 1//2

"""
Returns *half* the width of one step of the range.

This function is used to calculate the upper and lower bounds of `OnCell` interpolation objects.
""" boundstep

coordlookup(r::LinSpace, x) = (r.divisor * x + r.stop - r.len * r.start) / (r.stop - r.start)
coordlookup(r::FloatRange, x) = (r.divisor * x - r.start) / r.step + one(eltype(r))
coordlookup(r::StepRange, x) = (x - r.start) / r.step + one(eltype(r))
coordlookup(r::UnitRange, x) = x - r.start + one(eltype(r))
coordlookup(i::Bool, r::Range, x) = i ? coordlookup(r, x) : convert(typeof(coordlookup(r,x)), x)

basetype{T,N,ITPT,IT,GT,RT}(::Type{ScaledInterpolation{T,N,ITPT,IT,GT,RT}}) = ITPT
basetype(sitp::ScaledInterpolation) = basetype(typeof(sitp))

# @eval uglyness required for disambiguation with method in b-splies/indexing.jl
# also, GT is only specified to avoid disambiguation warnings on julia 0.4
gradient{T,N,ITPT,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType}}(sitp::ScaledInterpolation{T,N,ITPT,IT,GT}, xs::Real...) =
        gradient!(Array(T,count_interp_dims(IT,N)), sitp, xs...)
gradient{T,N,ITPT,IT<:DimSpec{InterpolationType},GT<:DimSpec{GridType}}(sitp::ScaledInterpolation{T,N,ITPT,IT,GT}, xs...) =
        gradient!(Array(T,count_interp_dims(IT,N)), sitp, xs...)
@generated function gradient!{T,N,ITPT,IT}(g, sitp::ScaledInterpolation{T,N,ITPT,IT}, xs::Number...)
    ndims(g) == 1 || throw(DimensionMismatch("g must be a vector (but had $(ndims(g)) dimensions)"))
    length(xs) == N || throw(DimensionMismatch("Must index into $N-dimensional scaled interpolation object with exactly $N indices (you used $(length(xs)))"))

    interp_types = length(IT.parameters) == N ? IT.parameters : tuple([IT.parameters[1] for _ in 1:N]...)
    interp_dimens = map(it -> interp_types[it] != NoInterp, 1:N)
    interp_indices = map(i -> interp_dimens[i] ? :(coordlookup(sitp.ranges[$i], xs[$i])) : :(xs[$i]), 1:N)

    quote
        length(g) == $(count_interp_dims(IT, N)) || throw(ArgumentError(string("The length of the provided gradient vector (", length(g), ") did not match the number of interpolating dimensions (", $(count_interp_dims(IT, N)), ")")))
        gradient!(g, sitp.itp, $(interp_indices...))
        for i in eachindex(g)
            g[i] = rescale_gradient(sitp.ranges[i], g[i])
        end
        g
    end
end

rescale_gradient(r::LinSpace, g) = g * r.divisor / (r.stop - r.start)
rescale_gradient(r::FloatRange, g) = g * r.divisor / r.step
rescale_gradient(r::StepRange, g) = g / r.step
rescale_gradient(r::UnitRange, g) = g

"""
`rescale_gradient(r::Range)`

Implements the chain rule dy/dx = dy/du * du/dx for use when calculating gradients with scaled interpolation objects.
""" rescale_gradient

### Iteration
type ScaledIterator{CR<:CartesianRange,SITPT,X1,Deg,T}
    rng::CR
    sitp::SITPT
    dx_1::X1
    nremaining::Int
    fx_1::X1
    itp_tail::NTuple{Deg,T}
end

"""
`eachvalue(sitp)` constructs an iterator for efficiently visiting each
grid point of a ScaledInterpolation object in which a small grid is
being "scaled up" to a larger one.  For example, suppose you have a
core `BSpline` object defined on a 5x7x4 grid, and you are scaling it
to a 100x120x20 grid (via `linspace(1,5,100), linspace(1,7,120),
linspace(1,4,20)`).  You can perform interpolation at each of these
grid points via

```
    function foo!(dest, sitp)
        i = 0
        for s in eachvalue(sitp)
            dest[i+=1] = s
        end
        dest
    end
```

which should be more efficient than

```
    function bar!(dest, sitp)
        for I in CartesianRange(size(dest))
            dest[I] = sitp[I]
        end
        dest
    end
```
"""
@generated function eachvalue{T,N}(sitp::ScaledInterpolation{T,N})
    ITPT = basetype(sitp)
    IT = itptype(ITPT)
    itp_tail = ntuple(i->zero(getindex_return_type(ITPT, ntuple(i->Int, N-1))), nelements(bsplinetype(iextract(IT, 1))))
    quote
        dx_1 = coordlookup(sitp.ranges[1], 2) - coordlookup(sitp.ranges[1], 1)
        ScaledIterator(CartesianRange(ssize(sitp)), sitp, dx_1, 0, zero(dx_1), $itp_tail)
    end
end

start(iter::ScaledIterator) = start(iter.rng)
done(iter::ScaledIterator, state) = done(iter.rng, state)

@generated function next{CR,ITPT,N}(iter::ScaledIterator{CR,ITPT}, state::CartesianIndex{N})
    value_expr = next_gen(iter)
    quote
        $value_expr
        (value, next(iter.rng, state)[2])
    end
end

ssize{T,N}(sitp::ScaledInterpolation{T,N}) = map(r->round(Int, last(r)-first(r)+1), sitp.ranges)::NTuple{N,Int}

nelements(::Union{Type{NoInterp},Type{Constant}}) = 1
nelements(::Type{Linear}) = 2
nelements{Q<:Quadratic}(::Type{Q}) = 3

function next_gen{CR,SITPT,X1,Deg,T}(::Type{ScaledIterator{CR,SITPT,X1,Deg,T}})
    N = ndims(CR)
    ITPT = basetype(SITPT)
    IT = itptype(ITPT)
    BS1 = iextract(IT, 1)
    BS1 == NoInterp && error("eachvalue is not implemented (and does not make sense) for NoInterp along the first dimension")
    pad = padding(ITPT)
    x_syms = [Symbol("x_", i) for i = 1:N]
    interp_index(IT, i) = iextract(IT, i) != NoInterp ?
        :($(x_syms[i]) = coordlookup(sitp.ranges[$i], state[$i])) :
        :($(x_syms[i]) = state[$i])
    # Calculations for the first dimension
    interp_index1 = interp_index(IT, 1)
    indices1 = define_indices_d(BS1, 1, padextract(pad, 1))
    coefexprs1 = coefficients(BS1, N, 1)
    nremaining_expr = nremaining_gen(BS1)
    # Calculations for the rest of the dimensions
    interp_indices_tail = map(i -> interp_index(IT, i), 2:N)
    indices_tail = [define_indices_d(iextract(IT, i), i, padextract(pad, i)) for i = 2:N]
    coefexprs_tail = [coefficients(iextract(IT, i), N, i) for i = 2:N]
    value_exprs_tail = index_gen_tail(BS1, IT, N)
    quote
        sitp = iter.sitp
        itp = sitp.itp
        if iter.nremaining > 0
            iter.nremaining -= 1
            iter.fx_1 += iter.dx_1
        else
            range1 = sitp.ranges[1]
            $interp_index1
            $indices1
            iter.nremaining = $nremaining_expr
            iter.fx_1 = fx_1
            $(interp_indices_tail...)
            $(indices_tail...)
            $(coefexprs_tail...)
            @inbounds iter.itp_tail = ($(value_exprs_tail...),)
        end
        fx_1 = iter.fx_1
        $coefexprs1
        $(index_gen1(BS1))
    end
end

function index_gen1(::Union{Type{NoInterp}, Type{BSpline{Constant}}})
    quote
        value = iter.itp_tail[1]
    end
end

function index_gen1(::Type{BSpline{Linear}})
    quote
        p = iter.itp_tail
        value = c_1*p[1] + cp_1*p[2]
    end
end

function index_gen1{Q<:Quadratic}(::Type{BSpline{Q}})
    quote
        p = iter.itp_tail
        value = cm_1*p[1] + c_1*p[2] + cp_1*p[3]
    end
end

function index_gen_tail{IT}(B::Union{Type{NoInterp}, Type{BSpline{Constant}}}, ::Type{IT}, N)
    [index_gen(B, IT, N, 0)]
end

function index_gen_tail{IT}(::Type{BSpline{Linear}}, ::Type{IT}, N)
    [index_gen(BS1, IT, N, i) for i = 0:1]
end

function index_gen_tail{IT,Q<:Quadratic}(::Type{BSpline{Q}}, ::Type{IT}, N)
    [index_gen(BSpline{Q}, IT, N, i) for i = -1:1]
end

function nremaining_gen{Q<:Quadratic}(::Union{Type{BSpline{Constant}}, Type{BSpline{Q}}})
    quote
        EPS = 0.001*iter.dx_1
        floor(Int, iter.dx_1 >= 0 ?
              (min(length(range1)+EPS, round(Int,x_1) + 0.5) - x_1)/iter.dx_1 :
              (max(1-EPS, round(Int,x_1) - 0.5) - x_1)/iter.dx_1)
    end
end

function nremaining_gen(::Type{BSpline{Linear}})
    quote
        EPS = 0.001*iter.dx_1
        floor(Int, iter.dx_1 >= 0 ?
              (min(length(range1)+EPS, floor(Int,x_1) + 1) - x_1)/iter.dx_1 :
              (max(1-EPS, floor(Int,x_1)) - x_1)/iter.dx_1)
    end
end


text:
# rand_inertia() = (A=randn(3,3);return A*A')

# function create_random_mechanism(nbodies; require_fixed_joint=false)
#     origin = Origin()
#     bodies = [Body(rand(),rand_inertia()) for i=1:nbodies]
#     parent_list = [origin;bodies]
#     joints = JointConstraint{Float64}[]
#     fixed_joint_flag = false

#     while !fixed_joint_flag
#         joints = JointConstraint{Float64}[]
#         for i=1:nbodies
#             parent = parent_list[rand(1:i)]
#             push!(joints, JointConstraint(
#                 (Translational{Float64,rand(0:3)}(parent, bodies[i]; parent_vertex=randn(3), child_vertex=randn(3), axis=randn(3)),
#                 Rotational{Float64,rand(0:3)}(parent, bodies[i]; orientation_offset=Quaternion(LinearAlgebra.normalize(randn(4))...)))))
#         end

#         fixed_joint_flag = !require_fixed_joint || any(typeof.(joints) .<: JointConstraint{T,6} where T)
#     end

#     return Mechanism(origin,bodies,joints)
# end

# function create_random_mechanism(nbodies; require_fixed_joint=false)
#     origin = Origin()
#     bodies = [Body(1.0,I(3)) for i=1:nbodies]
#     parent_list = [origin;bodies]
#     joints = JointConstraint{Float64}[]
#     fixed_joint_flag = false

#     while !fixed_joint_flag
#         joints = JointConstraint{Float64}[]
#         for i=1:nbodies
#             parent = parent_list[rand(1:i)]
#             push!(joints, JointConstraint(
#                 (Translational{Float64,rand(0:3)}(parent, bodies[i]; parent_vertex=zeros(3), child_vertex=zeros(3), axis=[1;0;0]),
#                 Rotational{Float64,rand(0:3)}(parent, bodies[i]; ))))
#         end

#         fixed_joint_flag = !require_fixed_joint || any(typeof.(joints) .<: JointConstraint{T,6} where T)
#     end

#     return Mechanism(origin,bodies,joints)
# end

# @testset "Merge Fixed Joints" begin
#     mech_full = create_random_mechanism(2, require_fixed_joint=true)
#     origin_reduced, bodies_reduced, joints_reduced = Dojo.reduce_fixed_joints(deepcopy(mech_full.origin),deepcopy.(mech_full.bodies),deepcopy.(mech_full.joints))
#     mech_reduced = Mechanism(origin_reduced,bodies_reduced,joints_reduced)
#     storage_full = simulate!(mech_full,5;record=true)
#     storage_reduced = simulate!(mech_reduced,5;record=true)
#     iddict = Dict{Int64,Int64}()
#     for (i_r,body_reduced) in enumerate(mech_reduced.bodies)
#         for (i_f,body_full) in enumerate(mech_full.bodies)
#             if string(body_reduced.name)[1:9] == string(body_full.name)
#                 push!(iddict,i_r => i_f)
#                 break
#             end
#         end
#     end
#     for pair in pairs(iddict)
#         # display(norm(norm.(storage_reduced.x[pair[1]] .- storage_full.x[pair[2]])))
#         # display(norm(norm.(storage_reduced.q[pair[1]] .- storage_full.q[pair[2]])))
#         @test all(norm.(storage_reduced.x[pair[1]] .- storage_full.x[pair[2]]) .< 1e-8)
#         @test all(norm.(storage_reduced.q[pair[1]] .- storage_full.q[pair[2]]) .< 1e-8)
#     end
# end

@testset "Pendulum" begin
    # get pendulum environment and simulate
    timestep=0.1
    env = DojoEnvironments.get_environment("pendulum";
        timestep,
        gravity=-10.0);
    reset(env);
    DojoEnvironments.initialize_pendulum!(env.mechanism,
        angle=0.25 * π)
    u1 = rand(Dojo.input_dimension(env.mechanism))
    z1 = Dojo.get_maximal_state(env.mechanism)

    storage = Dojo.simulate!(env.mechanism, 1.0,
        record=true)
    @test norm(z1 - Dojo.get_maximal_state(storage, 1)) < 1.0e-6
    zTs = Dojo.get_maximal_state(env.mechanism)
    @test norm(zTs - Dojo.get_maximal_state(storage, Dojo.length(storage))) < 1.0e-6
    zT = Dojo.get_next_state(env.mechanism)
    DojoEnvironments.initialize_pendulum!(env.mechanism,
        angle=0.25 * π)
    storage = Dojo.simulate!(env.mechanism, 1.0 + timestep,
        record=true)
    @test norm(zT - Dojo.get_maximal_state(storage, Dojo.length(storage))) < 1.0e-6

    # maximal gradients
    Dojo.step!(env.mechanism, z1, u1)
    Fx1, Fu1 = Dojo.get_maximal_gradients(env.mechanism)
    Fx2, Fu2 = Dojo.get_maximal_gradients!(env.mechanism, z1, u1)
    @test norm(Fx1 - Fx2, Inf) < 1.0e-6
    @test norm(Fu1 - Fu2, Inf) < 1.0e-6

    timestep=0.1
    env = DojoEnvironments.get_environment("halfcheetah";
        timestep,
        gravity=-10.0);
    reset(env);

    # get body
    @test Dojo.get_body(env.mechanism, :origin).name == :origin

    # get contact
    contact_name = env.mechanism.contacts[1].name
    @test Dojo.get_contact(env.mechanism, contact_name).name == contact_name

    # velocity indices
    @test Dojo.velocity_index(env.mechanism) == [4, 5, 6, 8, 10, 12, 14, 16, 18]

    # set state
    z = Dojo.get_maximal_state(env.mechanism)
    env2 = DojoEnvironments.get_environment(:halfcheetah;
        timestep,
        gravity=-10.0);
    reset(env);
    Dojo.set_maximal_state!(env2.mechanism, zeros(Dojo.maximal_dimension(env2.mechanism)))
    z2 = Dojo.get_maximal_state(env2.mechanism)
    @test norm(z2) < 1.0e-8
    Dojo.set_maximal_state!(env2.mechanism, z)
    @test norm(z - Dojo.get_maximal_state(env2.mechanism)) < 1.0e-8

    # change atlas floating base
    mech = get_mechanism(:atlas; parse_damper=false)
    j = get_joint(mech, :floating_base)
    base_origin = get_body(mech, j.child_id).name

    mech = set_floating_base(mech, :head)
    j = get_joint(mech, :floating_base)
    base_new = get_body(mech, j.child_id).name

    @test base_origin != base_new
    @test base_new == :head
end

#TODO: get and set methods


text:
function create_initial_condition(kx, kz, k2cut, N, parameters::GeneralParams)
    datah, datap, dataq, U = initial_condition(parameters);

    datas1s2r1r2 = zeros(length(datah))'

    # remove high frequency parts by aliasing filter
    datah = datah .* k2cut;
    datap = datap .* k2cut;
    dataq = dataq .* k2cut;

    # reshape data for U vector
    datahSV = reshape(datah', 1, :)
    datapSV = reshape(datap', 1, :)
    dataqSV = reshape(dataq', 1, :)

    U = vcat(datahSV[1:N ÷ 2], dataqSV[1:N ÷ 2], datapSV[1:N ÷ 2])
    U = reshape(U, 1, :)
    return U
end

function initial_condition(parameters::GeneralParams)
    @unpack NN1, NN2 = parameters
    @unpack stSine, spSine, stCosine, spCosine, stNwaves, spNwaves, stNoise, spNoise = parameters
    @unpack k0x, k0z, h0, Lx, Lz = parameters

    Ax = stSine
    Az = spSine
    Bx = stCosine
    Bz = spCosine
    nwx = stNwaves
    nwz = spNwaves

    if (NN1 == 0) k0x = 0; Lx = 0; Ax = 0; Bx = 0 end
    if (NN2 == 0) k0z = 0; Lz = 0; Az = 0; Bz = 0 end

    N1 = 2^NN1; M1 = 2 * N1; N2 = 2^NN2; M2 = 2 * N2; N = 4 * N2 * (N1 + 1);

    datah = zeros((M1, M2))
    dataq = zeros((M1, M2))
    datap = zeros((M1, M2))

    randz = zeros(M2)
    for n = 1:M2
        randz[n] = spNoise * rand()
    end

    for m = 1:M1
        x = (m - 1) * Lx / M1
        for n = 1:M2
            z = (n - 1) * Lz / M2
            datah[m, n] = h0 * (1 + Ax * sin(nwx * k0x * x) + Bx * cos(nwx * k0x * x) + Az * sin(nwz * k0z * z) + Bz * cos(nwz * k0z * z) + stNoise * rand()) + randz[n]
            datap[m, n] = 0
            dataq[m, n] = datah[m, n].^3 / 3
        end
    end

    # calculate mass
    mass = sum(datah)
    mass = mass / (M1 * M2)

    # set mass to h0 for random function
    if spNoise != 0
        datah = datah - mass + h0
    end

    # Fourier transform of the initial condition
    datah = fft(datah)
    datap = fft(datap)
    dataq = fft(dataq)

    datah = conj(datah)
    datap = conj(datap)
    dataq = conj(dataq)

    datahSV = vec(datah')
    datapSV = vec(datap')
    dataqSV = vec(dataq')

    M1M2 = M1 * M2;

    # note ÷ for integer division
    U = vcat(datahSV[1:N ÷ 2], dataqSV[1:N ÷ 2], datapSV[1:N ÷ 2]) ./ M1M2

    datah = datah ./ M1M2;
    dataq = dataq ./ M1M2;
    datap = datap ./ M1M2;

    return datah, datap, dataq, U

end

function aliasing_filter(radius::Float64, parameters::GeneralParams)
    @unpack k0x, k0z, N1, N2, M1, M2 = parameters;

    m = round(radius * N1)
    n = round(radius * N2)
    k2cutCriterion = (m * k0x)^2 + ((-N2 + mod(n + N2, M2)) * k0z)^2
    kx, kz, k2, k2cut = [zeros((N1 + 1, M2)) for i = 1:4]  # initialise
    for m = 1:N1 + 1
        for n = 1:M2
            kx[m, n] = (m - 1) * k0x
            kz[m, n] = (-N2 + mod((n - 1) + N2, M2)) * k0z
            k2[m, n] = kx[m, n] * kx[m, n] + kz[m, n] * kz[m, n]
            if (k2[m, n] > k2cutCriterion)
                k2cut[m, n] = 0
            else
                k2cut[m, n] = 1
            end
        end
    end
    k2cut = [k2cut; [conj(reverse(k2cut[2:N1,1], dims=1)) conj(reverse(reverse(k2cut[2:N1,2:M2], dims=1), dims=2))]]
    return kx, kz, k2cut
end


text:
#=
# Examples
=#

import Pastebin

# Fetch a paste by key.

response = Pastebin.getraw("Lt9PtFc4")
Pastebin.content(response)

# Create a Pastebin client using only a `devKey`.

client = Pastebin.Client(ENV["PASTEBIN_DEVKEY"])

# Paste something (to expire in 1 hour)

response = Pastebin.paste(client, "A test", "Hola Pastebin"; expire="1H")

# Extract the key of the created paste

key = Pastebin.parsekey(response)

# We can now access the created paste using this key

Pastebin.content(Pastebin.getraw(key))


text:
using LinearAlgebra

# k_local(stiffness coefficient,number of nodes,row/column one, row/column two, ==1 adjusts for fixed node)

#creates local stiffness matrix and maps to global
function k_local(k,i,n1,n2,fixed=0)
      k_global = zeros(i,i)
      k_global[n1,n1] = k
      k_global[n1,n2] = -k
      k_global[n2,n1] = -k
      k_global[n2,n2] = k
      if fixed == 1
        adj = zeros(1,i)
        adj[n1] = 1
        k_global[n1,:] = adj
        k_global[:,n1] = adj
      end
      display(k_global)
      return k_global
  end

# constants from problem statement
k1 = 5
k2 = 8
k3 = 5
k4 = 20
k5 = 10
k6 = 20
# combining all globally mapped matrices
kG = k_local(k1,5,1,2,1)+
     k_local(k2,5,2,3)+
     k_local(k3,5,2,3)+
     k_local(k4,5,2,4)+
     k_local(k5,5,3,4)+
     k_local(k6,5,5,4,1)

# Force Vector
f = Float64[0 10 0 10 0]
F = transpose(f)

# Displacement vector
u = kG\F
println("Nodal displacement vector (inches):")
display(u)

# Strain Function
ϵ(x2,x1,l) = (x2 - x1)/l

#Reaction Forces vector
Fapp(k,x2,x) = k*(x2-x)
FappM = [Fapp(k1,u[2],u[1])
         Fapp(k2,u[3],u[2])
         Fapp(k3,u[3],u[2])
         Fapp(k4,u[4],u[2])
         Fapp(k5,u[4],u[3])
         Fapp(k6,u[5],u[4])]
println("Reaction Forces vector (lbs):")
display(FappM)


text:
using Trixi, OrdinaryDiffEq

dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = SBP(),
             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
             volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))

equations = CompressibleEulerEquations2D(1.4)
initial_condition = initial_condition_convergence_test
source_terms = source_terms_convergence_test

# example where we tag two separate boundary segments of the mesh
top_boundary(x, tol=50*eps()) = abs(x[2]-1)<tol
rest_of_boundary(x, tol=50*eps()) = !top_boundary(x, tol)
is_on_boundary = Dict(:top => top_boundary, :rest => rest_of_boundary)

function mapping(xi, eta)
  x = xi  + 0.1 * sin(pi * xi) * sin(pi * eta)
  y = eta + 0.1 * sin(pi * xi) * sin(pi * eta)
  return SVector(x, y)
end
cells_per_dimension = (16, 16)
vertex_coordinates, EToV = StartUpDG.uniform_mesh(dg.basis.elementType, cells_per_dimension...)
for i in eachindex(vertex_coordinates[1])
  vx, vy = getindex.(vertex_coordinates, i)
  setindex!.(vertex_coordinates, mapping(vx, vy), i)
end

mesh = DGMultiMesh(vertex_coordinates, EToV, dg, is_on_boundary=is_on_boundary)

boundary_condition_convergence_test = BoundaryConditionDirichlet(initial_condition)
boundary_conditions = (; :top => boundary_condition_convergence_test,
                        :rest => boundary_condition_convergence_test)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
                                    source_terms = source_terms,
                                    boundary_conditions = boundary_conditions)

tspan = (0.0, 0.4)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()
alive_callback = AliveCallback(alive_interval=10)
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)

###############################################################################
# run the simulation

sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-6, reltol=1.0e-6,
            save_everystep=false, callback=callbacks);

summary_callback() # print the timer summary