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

# Extend _add_data_object
function _add_data_object(model::InfiniteModel,
                          object::MultiParameterData
                          )::DependentParametersIndex
    index = MOIUC.add_item(model.dependent_params, object)
    push!(model.param_object_indices, index)
    return index
end

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

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

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

# Extend _core_variable_object
function _core_variable_object(pref::DependentParameterRef)::DependentParameters
    return _data_object(pref).parameters
end

# Return the number of dependent parameters involved
function _num_parameters(pref::DependentParameterRef)::Int
    return length(_data_object(pref).names)
end

# Extend _delete_data_object
function _delete_data_object(vref::DependentParameterRef)::Nothing
    delete!(_data_dictionary(vref), JuMP.index(vref).object_index)
    return
end

################################################################################
#                             PARAMETER DEFINITION
################################################################################
# Store partially processed individual dependent parameters
struct _DependentParameter{S <: AbstractInfiniteSet, M <: AbstractDerivativeMethod}
    set::S
    supports::Vector{Float64}
    name::String
    deriv_method::M
    function _DependentParameter(set::S,
        supports::Union{Vector{<:Real}, Real},
        name::String,
        method::M
        )::_DependentParameter{S, M} where {S <: AbstractInfiniteSet, M <: AbstractDerivativeMethod}
        if supports isa Real
            return new{S, M}(set, [supports], name, method)
        else
            return new{S, M}(set, supports, name, method)
        end
    end
end

## Use type dispatch to efficiently check the set(s)
# All InfiniteScalarSets
function _check_param_sets(_error::Function,
    params::AbstractArray{<:_DependentParameter{<:InfiniteScalarSet}}
    )::Nothing
    return
end

# All MultiDistributionSet with non SparseAxisArray
function _check_param_sets(_error::Function,
    params::AbstractArray{<:_DependentParameter{<:MultiDistributionSet}}
    )::Nothing
    dist = first(params).set.distribution
    set_type = typeof(first(params).set)
    if size(dist) != size(params)
        _error("The dimensions of the parameters and the multi-dimensional " *
               "distribution $dist.")
    end
    return
end

# All MultiDistributionSet with SparseAxisArray
function _check_param_sets(_error::Function,
    params::JuMPC.SparseAxisArray{<:_DependentParameter{<:MultiDistributionSet}}
    )
    _error("Cannot specify multiple-dimensional distribution set with a " *
           "`SparseAxisArray` of dependent infinite parameters.")
end

# All CollectionSet
function _check_param_sets(_error::Function,
    params::AbstractArray{<:_DependentParameter{<:CollectionSet}}
    )::Nothing
    sets = collection_sets(first(params).set)
    set_type = typeof(first(params).set)
    if length(sets) != length(params)
        _error("The dimensions of the parameters and the specified CollectionSet " *
               "do not match.")
    elseif params isa JuMPC.SparseAxisArray
        @warn("CollectionSet order may not match the given `SparseAxisArray` " *
              "of specified dependent infinite parameters, consider instead " *
              "specifying the `InfiniteScalarSet` for each parameter using " *
              "the `set` keyword and the appropriate indices.")
    end
    return
end

# All some InfiniteArraySet (for extensions)
function _check_param_sets(_error::Function,
    params::AbstractArray{<:_DependentParameter{<:InfiniteArraySet}}
    )::Nothing
    return
end

# Mixed sets
function _check_param_sets(_error::Function,
    params::AbstractArray{<:_DependentParameter}
    )::Nothing
    set_type = typeof(first(params).set)
    if !all(param.set isa InfiniteScalarSet for param in params)
        _error("Cannot specify multiple `InfiniteScalarSets` for one container " *
               "of infinite dependent parameters.")
    end
    return
end

# Fallback
function _check_param_sets(_error::Function, params)
    _error("Unrecognized input for infinite set.")
end

## Define methods for checking the the derivative methods
# Expected format
function _check_derivative_methods(_error::Function,
    params::AbstractArray{<:_DependentParameter})::Nothing
    if !all(p.deriv_method isa NonGenerativeDerivativeMethod for p in params)
        _error("Cannot use generative derivative evaluation methods with dependent " *
               "infinite parameters. Only subtypes of `NonGenerativeDerivativeMethod` " *
               "can be used.")
    end
    return
end

# Fallback
function _check_derivative_methods(_error::Function, params)
    _error("Unrecognized input for derivative method.")
end

## Use set type dispatch to make the proper InfiniteArraySet
# InfiniteArraySet
function _make_array_set(params::Vector{<:_DependentParameter{T}}
                         )::T where {T <: InfiniteArraySet}
    return first(params).set
end

# InfiniteScalarSets
function _make_array_set(params::Vector{<:_DependentParameter{T}}
                         )::CollectionSet{T} where {T <: InfiniteScalarSet}
    return CollectionSet([p.set for p in params])
end

# Build a DependentParameters object given an array of _DependentParameters
function _build_parameters(_error::Function,
                           params::AbstractArray{<:_DependentParameter};
                           num_supports::Int = 0,
                           sig_digits::Int = DefaultSigDigits,
                           extra_kw_args...)
    # error with extra keywords
    for (kwarg, _) in extra_kw_args
       _error("Unrecognized keyword argument $kwarg")
    end
    # check the formatting
    _check_param_sets(_error, params)
    _check_derivative_methods(_error, params)
    # vectorize the parameter array
    indices = Collections._get_indices(params)
    ordered_params = Collections._make_ordered(params, indices)
    vector_params = _make_vector(ordered_params)
    # make the set
    set = _make_array_set(vector_params)
    # make the supports and labels
    lens = map(p -> length(p.supports), vector_params)
    _allequal(lens) || _error("Inconsistent support dimensions.")
    # we have supports
    if first(lens) != 0
        # build the support array transpose to fill in column order (leverage locality)
        trans_supps = Array{Float64}(undef, first(lens), length(vector_params))
        for i in 1:size(trans_supps, 2)
            trans_supps[:, i] = vector_params[i].supports
        end
        supps = permutedims(trans_supps)
        supports_in_set(supps, set) || _error("Supports violate infinite set domain.")
        supps = round.(supps, sigdigits = sig_digits)
        label = UserDefined
        supp_dict = Dict{Vector{Float64}, Set{DataType}}(@views supps[:, i] =>
                                         Set([label]) for i in 1:size(supps, 2))
    # we want to generate supports
    elseif num_supports != 0
        supps, label = generate_support_values(set, num_supports = num_supports,
                                               sig_digits = sig_digits)
        supp_dict = Dict{Vector{Float64}, Set{DataType}}(@views supps[:, i] =>
                                         Set([label]) for i in 1:size(supps, 2))
    # no supports are specified
    else
        supp_dict = Dict{Vector{Float64}, Set{DataType}}()
    end
    # make the parameter object
    names = map(p -> p.name, vector_params)
    methods = map(p -> p.deriv_method, vector_params)
    return DependentParameters(set, supp_dict, sig_digits, methods), names, indices
end

"""
    add_parameters(model::InfiniteModel,
                   params::DependentParameters,
                   names::Vector{String} = ["noname", "noname", ...],
                   indices = nothing
                   )::AbstractArray{<:GeneralVariableRef}

Add `params` to `model` and return an appropriate container of the dependent
infinite parameter references. This is intended as an internal method for use
with [`@dependent_parameters`](@ref). However, if desired users can use this
to create a container of infinite dependent parameter without the use of a
macro. `names` denote the name of each parameter and `indices` denote the
indices of the expected container as used by `Containers._make_array`
(implemented by `VectorTuple`s), by default a `Vector` is returned.

**Example**
```julia-repl
julia> using Distributions

julia> dist = MvNormal(ones(3)); # 3 dimensional

julia> set = MultiDistributionSet(dist); # 3 dimensional

julia> params = DependentParameters(set, Dict{Vector{Float64}, Set{DatatType}}(), 10);

julia> prefs = add_parameters(model, params, ["par1", "par2", "par3"])
3-element Array{GeneralVariableRef,1}:
 par1
 par2
 par3
```
"""
function add_parameters(model::InfiniteModel,
                        params::DependentParameters,
                        names::Vector{String} = String[],
                        indices = nothing
                        )::AbstractArray{<:GeneralVariableRef}
    # get the number of parameters
    num_params = length(params.set)
    # process the names
    if isempty(names)
        names = ["noname" for i in 1:num_params]
    end
    # process the indices
    if indices === nothing
        indices = CartesianIndices(1:num_params)
    end
    # make the parameter model object
    obj_num = length(_param_object_indices(model)) + 1
    first_param_num = model.last_param_num + 1
    last_param_num = model.last_param_num += num_params
    param_nums = first_param_num:last_param_num
    data_object = MultiParameterData(params, obj_num, param_nums, names)
    # add the data object to the model and make the references
    obj_index = _add_data_object(model, data_object)
    prefs = [GeneralVariableRef(model, obj_index.value, DependentParameterIndex, i)
             for i in 1:num_params]
    return Collections._make_array(prefs, indices)
end

# Construct an expression to build an infinite set(s) (use with @dependent_parameters)
function _construct_array_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 real number."))
    elseif info.has_dist
        check = :(isa($(info.distribution), Distributions.UnivariateDistribution))
        return :($(check) ? UniDistributionSet($(info.distribution)) : MultiDistributionSet($(info.distribution)))
    elseif info.has_set
        check1 = :(isa($(info.set), AbstractInfiniteSet))
        check2 = :(isa($(info.set), Distributions.UnivariateDistribution))
        return :($(check1) ? $(info.set) : ($(check2) ? UniDistributionSet($(info.set)) : MultiDistributionSet($(info.set))))
    else
        _error("Must specify upper/lower bounds, a distribution, or a set")
    end
end

################################################################################
#                                  NAMING
################################################################################
# Get the parameter index in the DependentParameters object
_param_index(pref::DependentParameterRef)::Int = JuMP.index(pref).param_index

"""
    JuMP.name(pref::DependentParameterRef)::String

Extend [`JuMP.name`](@ref JuMP.name(::JuMP.VariableRef)) to return the names of
infinite dependent parameters.

**Example**
```julia-repl
julia> name(pref)
"par_name"
```
"""
function JuMP.name(pref::DependentParameterRef)::String
    object = get(_data_dictionary(pref), JuMP.index(pref).object_index, nothing)
    return object === nothing ? "" : object.names[_param_index(pref)]
end

"""
    JuMP.set_name(pref::DependentParameterRef, name::String)::Nothing

Extend [`JuMP.set_name`](@ref JuMP.set_name(::JuMP.VariableRef, ::String)) to set
names of dependent infinite parameters.

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

julia> name(vref)
"para_name"
```
"""
function JuMP.set_name(pref::DependentParameterRef, name::String)::Nothing
    _data_object(pref).names[_param_index(pref)] = name
    JuMP.owner_model(pref).name_to_param = nothing
    return
end

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

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

# Extend _measure_dependencies
function _measure_dependencies(pref::DependentParameterRef
    )::Vector{MeasureIndex}
    return _data_object(pref).measure_indices[_param_index(pref)]
end

# Extend _constraint_dependencies
function _constraint_dependencies(pref::DependentParameterRef
    )::Vector{ConstraintIndex}
    return _data_object(pref).constraint_indices[_param_index(pref)]
end

# Extend _derivative_dependencies
function _derivative_dependencies(pref::DependentParameterRef
    )::Vector{DerivativeIndex}
    return _data_object(pref).derivative_indices[_param_index(pref)]
end

"""
    used_by_infinite_variable(pref::DependentParameterRef)::Bool

Return a `Bool` indicating if the dependent infinite parameter `pref` is used by
an infinite variable.

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

"""
    used_by_parameter_function(pref::DependentParameterRef)::Bool

Return a `Bool` indicating if the dependent infinite parameter `pref` is used by
an infinite parameter function.

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

"""
    used_by_measure(pref::DependentParameterRef)::Bool

Return a `Bool` indicating if the dependent infinite parameter `pref` is used
by a measure.

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

"""
    used_by_constraint(pref::DependentParameterRef)::Bool

Return a `Bool` indicating if the dependent infinite parameter `pref` is used by
a constraint.

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

"""
    used_by_derivative(pref::DependentParameterRef)::Bool

Return a `Bool` indicating if the dependent infinite parameter `pref` is used by
a derivative.

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

# Extend used by objective
used_by_objective(pref::DependentParameterRef)::Bool = false

"""
    is_used(pref::DependentParameterRef)::Bool

Return a `Bool` indicating if the dependent infinite parameter `pref` is used in
the model.

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

################################################################################
#                          PARAMETER OBJECT METHODS
################################################################################
# Extend _parameter_number
function _parameter_number(pref::DependentParameterRef)::Int
    return _data_object(pref).parameter_nums[_param_index(pref)]
end

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

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

# Extend _object_numbers
function _object_numbers(pref::DependentParameterRef)::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::DependentParameterRef, params::P,
    data::MultiParameterData{P})::Nothing where {P <: DependentParameters}
    data.parameters = params
    return
end

# Reconstruction is necessary
function _adaptive_data_update(pref::DependentParameterRef, params::P1,
    data::MultiParameterData{P2})::Nothing  where {P1, P2}
    new_data = MultiParameterData(params, data.object_num, data.parameter_nums,
                                  data.names, data.parameter_func_indices,
                                  data.infinite_var_indices,
                                  data.derivative_indices, data.measure_indices,
                                  data.constraint_indices,
                                  data.has_internal_supports,
                                  data.has_deriv_constrs)
    _data_dictionary(pref)[JuMP.index(pref).object_index] = new_data
    return
end

# Extend _set_core_variable_object
function _set_core_variable_object(pref::DependentParameterRef,
                                   params::DependentParameters)::Nothing
    _adaptive_data_update(pref, params, _data_object(pref))
    return
end

################################################################################
#                        DERIVATIVE METHOD FUNCTIONS
################################################################################
# Extend fallback for dependent parameters
function has_derivative_supports(pref::DependentParameterRef)::Bool
    return false
end

# Extend fallback for dependent parameters
function _set_has_derivative_supports(pref::DependentParameterRef,
                                      status::Bool)::Nothing
    return
end

# Extend has derivative constraints
function has_derivative_constraints(pref::DependentParameterRef)::Bool
    return _data_object(pref).has_deriv_constrs[_param_index(pref)]
end

# Extend setting if has derivative constraints
function _set_has_derivative_constraints(pref::DependentParameterRef,
                                         status::Bool)::Nothing
    _data_object(pref).has_deriv_constrs[_param_index(pref)] = status
    return
end

# Get the raw derivative method vector
function _derivative_methods(pref::DependentParameterRef)
    return _core_variable_object(pref).derivative_methods
end

"""
    derivative_method(pref::DependentParameterRef)::NonGenerativeDerivativeMethod

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
```
"""
function derivative_method(pref::DependentParameterRef)::NonGenerativeDerivativeMethod
    return _derivative_methods(pref)[_param_index(pref)]
end

## Define helper methods for setting the derivative method efficiently
# Compatible with vector type
function _adaptive_method_update(pref,
    p::DependentParameters{S, M1},
    method::M2
    )::Nothing where {S, M1 <: NonGenerativeDerivativeMethod, M2 <: M1}
    p.derivative_methods[_param_index(pref)] = method
    return
end

# Not compatible
function _adaptive_method_update(pref,
    p::DependentParameters{S, M1},
    method::M2
    )::Nothing where {S, M1, M2}
    methods = p.derivative_methods
    new_methods = [i == _param_index(pref) ? method : m
                   for (i, m) in enumerate(methods)]
    new_params = DependentParameters(p.set, p.supports, p.sig_digits, new_methods)
    _set_core_variable_object(pref, new_params)
    return
end

"""
    set_derivative_method(pref::DependentParameterRef,
                          method::NonGenerativeDerivativeMethod)::Nothing

Specfies the desired derivative evaluation method `method` for derivatives that are
taken with respect to `pref`. Errors if `method` is generative (i.e., it requires
the definition of additional supports)

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

```
"""
function set_derivative_method(pref::DependentParameterRef,
    method::AbstractDerivativeMethod
    )::Nothing
    if !(method isa NonGenerativeDerivativeMethod)
        error("Must specify a subtype of `NonGenerativeDerivativeMethod` for " *
              "for a dependent parameter.")
    end
    _adaptive_method_update(pref, _core_variable_object(pref), method)
    _reset_derivative_evaluations(pref)
    if is_used(pref)
        set_optimizer_model_ready(JuMP.owner_model(pref), false)
    end
    return
end

"""
    set_all_derivative_methods(model::InfiniteModel,
                               method::AbstractDerivativeMethod)::Nothing

Sets the desired evaluation method `method` for all the derivatives currently added
to `model`. Note that this is done with respect to the infinite parameters. Errors
if a generative method is specified and the model contains dependent parameters.

**Example**
```julia-repl
julia> set_all_derivative_methods(model, OrthogonalCollocation(2))

```
"""
function set_all_derivative_methods(model::InfiniteModel,
    method::AbstractDerivativeMethod
    )::Nothing
    for pref in all_parameters(model, InfiniteParameter)
        set_derivative_method(pref, method)
    end
    return
end

################################################################################
#                             INFINITE SET METHODS
################################################################################
## Get the individual infinite set if possible
# raw_set
function _parameter_set(pref::DependentParameterRef)::InfiniteArraySet
    return _core_variable_object(pref).set
end

# CollectionSet
function _parameter_set(set::CollectionSet{S},
                        pref::DependentParameterRef
                        )::S where {S <: InfiniteScalarSet}
    return collection_sets(set)[_param_index(pref)]
end

# InfiniteArraySet (Fallback)
function _parameter_set(set::InfiniteArraySet, pref::DependentParameterRef)
    error("An individual infinite set is not well-defined for $pref which " *
          "is part of a group of dependent infinite parameters that correspond " *
          "to an multi-dimensional infinite set of type `$(typeof(set))`.")
end

"""
    infinite_set(pref::DependentParameterRef)::InfiniteScalarSet

Return the infinite set associated with the particular infinite dependent
parameter `pref` if valid. Errors if the underlying [`DependentParameters`](@ref)
object does not use a [`CollectionSet`](@ref).

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

# Check that prefs are complete
function _check_complete_param_array(
    prefs::AbstractArray{<:DependentParameterRef}
    )::Nothing
    if length(prefs) != _num_parameters(first(prefs))
        error("Dimensions of parameter container and the infinite set do not " *
              "match, ensure all related dependent parameters are included.")
    end
    return
end

"""
    infinite_set(prefs::AbstractArray{<:DependentParameterRef})::InfiniteArraySet

Return the infinite set associated with the container of infinite dependent
parameters `prefs`. Errors if the container `prefs` is incomplete.

**Example**
```julia-repl
julia> infinite_set(x)
ZeroMeanDiagNormal(
dim: 2
μ: [0.0, 0.0]
Σ: [1.0 0.0; 0.0 1.0]
)
```
"""
function infinite_set(prefs::AbstractArray{<:DependentParameterRef}
                      )::InfiniteArraySet
    _check_complete_param_array(prefs)
    return _parameter_set(first(prefs))
end

# Update the underlying set and delete the supports
function _update_parameter_set(pref::DependentParameterRef,
                               new_set::InfiniteArraySet)::Nothing
    old_params = _core_variable_object(pref)
    new_supports = Dict{Vector{Float64}, Set{DataType}}()
    sig_figs = significant_digits(pref)
    methods = _derivative_methods(pref)
    new_params = DependentParameters(new_set, new_supports, sig_figs, methods)
    _set_core_variable_object(pref, new_params)
    for i in 1:length(new_set)
        idx = DependentParameterIndex(JuMP.index(pref).object_index, i)
        p = DependentParameterRef(JuMP.owner_model(pref), idx)
        _reset_derivative_evaluations(p)
    end
    _set_has_internal_supports(pref, false)
    if is_used(pref)
        set_optimizer_model_ready(JuMP.owner_model(pref), false)
    end
    return
end

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

Specify the scalar infinite set of the dependent infinite parameter `pref` to
`set` if `pref` is part of a [`CollectionSet`](@ref), otherwise an error is
thrown. Note this will reset/delete all the supports contained in the
underlying [`DependentParameters`](@ref) object. Also, errors if `pref` is used
by a measure.

**Example**
```julia-repl
julia> set_infinite_set(x[1], IntervalSet(0, 2))

julia> infinite_set(x[1])
[0, 2]
```
"""
function set_infinite_set(pref::DependentParameterRef,
                          set::InfiniteScalarSet)::Nothing
    old_set = _parameter_set(pref)
    if !(old_set isa CollectionSet)
        error("Cannot set the individual infinite set of $pref if the " *
              "underlying set is not a CollectionSet.")
    elseif used_by_measure(pref)
        error("Cannot override the infinite set of $pref since it is used by " *
              "a measure.")
    end
    param_idx = _param_index(pref)
    new_set = CollectionSet([i != param_idx ? collection_sets(old_set)[i] : set
                             for i in eachindex(collection_sets(old_set))])
    _update_parameter_set(pref, new_set)
    return
end

"""
    set_infinite_set(prefs::AbstractArray{<:DependentParameterRef},
                     set::InfiniteArraySet)::Nothing

Specify the multi-dimensional infinite set of the dependent infinite parameters
`prefs` to `set`. Note this will reset/delete all the supports contained in the
underlying [`DependentParameters`](@ref) object. This will error if the not all
of the dependent infinite parameters are included, if any of them are used by
measures.

**Example**
```julia-repl
julia> set_infinite_set(x, CollectionSet([IntervalSet(0, 1), IntervalSet(0, 2)]))
```
"""
function set_infinite_set(prefs::AbstractArray{<:DependentParameterRef},
                          set::InfiniteArraySet)::Nothing
    if any(used_by_measure(pref) for pref in prefs)
        error("Cannot override the infinite set of $prefs since it is used by " *
              "a measure.")
    end
    _check_complete_param_array(prefs)
    _update_parameter_set(first(prefs), set)
    return
end

"""
    JuMP.has_lower_bound(pref::DependentParameterRef)::Bool

Extend the `JuMP.has_lower_bound` function to accomodate a single dependent
infinite parameter.
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 scalar infinite set types
should extend `JuMP.has_lower_bound(set::NewType)`.

**Example**
```julia-repl
julia> has_lower_bound(x[1])
true
```
"""
function JuMP.has_lower_bound(pref::DependentParameterRef)::Bool
    set = _parameter_set(pref)
    if set isa CollectionSet
        return JuMP.has_lower_bound(collection_sets(set)[_param_index(pref)])
    else
        return false
    end
end

"""
    JuMP.lower_bound(pref::DependentParameterRef)::Number

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

**Example**
```julia-repl
julia> lower_bound(x[1])
0.0
```
"""
function JuMP.lower_bound(pref::DependentParameterRef)::Number
    if !JuMP.has_lower_bound(pref)
        error("Parameter $(pref) does not have a lower bound.")
    end
    return JuMP.lower_bound(infinite_set(pref))
end

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

Extend the `JuMP.set_lower_bound` function to accomodate a single dependent
infinite parameter. Updates the infinite set lower bound if such an operation
is supported. Infinite scalar set extensions that seek to employ this should extend
`JuMP.set_lower_bound(set::NewType, lower::Number)`. This will call
[`set_infinite_set`](@ref) and will error if this is not well-defined. Note
that existing supports will be deleted.

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

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

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

Extend the `JuMP.has_upper_bound` function to accomodate a single dependent
infinite parameter.
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 scalar infinite set types
should extend `JuMP.has_upper_bound(set::NewType)`.

**Example**
```julia-repl
julia> has_upper_bound(x[1])
true
```
"""
function JuMP.has_upper_bound(pref::DependentParameterRef)::Bool
    set = _core_variable_object(pref).set
    if set isa CollectionSet
        return JuMP.has_upper_bound(collection_sets(set)[_param_index(pref)])
    else
        return false
    end
end

"""
    JuMP.upper_bound(pref::DependentParameterRef)::Number

Extend the `JuMP.upper_bound` function to accomodate a single dependent infinite
parameter. Returns the upper bound associated with the infinite set. Errors if
such a bound is not well-defined.

**Example**
```julia-repl
julia> upper_bound(x[1])
0.0
```
"""
function JuMP.upper_bound(pref::DependentParameterRef)::Number
    if !JuMP.has_upper_bound(pref)
        error("Parameter $(pref) does not have a upper bound.")
    end
    return JuMP.upper_bound(infinite_set(pref))
end

"""
    JuMP.set_upper_bound(pref::DependentParameterRef, upper::Real)::Nothing

Extend the `JuMP.set_upper_bound` function to accomodate a single dependent
infinite parameter. Updates the infinite set upper bound if such an operation
is supported. Infinite scalar set extensions that seek to employ this should extend
`JuMP.set_upper_bound(set::NewType, upper::Number)`. This will call
[`set_infinite_set`](@ref) and will error if this is not well-defined. Note
that existing supports will be deleted.

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

julia> upper_bound(t)
-1.0
```
"""
function JuMP.set_upper_bound(pref::DependentParameterRef, upper::Real)::Nothing
    set = infinite_set(pref)
    new_set = JuMP.set_upper_bound(set, upper)
    set_infinite_set(pref, new_set)
    return
end

################################################################################
#                              SUPPORT METHODS
################################################################################
# Get the raw supports
function _parameter_supports(pref::DependentParameterRef
                             )::Dict{Vector{Float64}, Set{DataType}}
    return _core_variable_object(pref).supports
end

"""
    significant_digits(pref::DependentParameterRef)::Int

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

**Example**
```julia-repl
julia> significant_digits(x[1])
12
```
"""
function significant_digits(pref::DependentParameterRef)::Int
    return _core_variable_object(pref).sig_digits
end

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

Return the number of support points associated with a single dependent infinite
parameter `pref`. Specify a subset of supports via `label` to only count the
supports with `label`. By default only the amount of public supports are given, but
the full amount is obtained via `label == All`.

**Example**
```julia-repl
julia> num_supports(x[1])
2

julia> num_supports(x[1], label = MCSample)
0
```
"""
function num_supports(pref::DependentParameterRef;
                      label::Type{<:AbstractSupportLabel} = PublicLabel)::Int
    supp_dict = _parameter_supports(pref)
    if label == All || (!has_internal_supports(pref) && label == PublicLabel)
        return length(supp_dict)
    else
        return count(p -> any(v -> v <: label, p[2]), supp_dict)
    end
end

"""
    num_supports(prefs::AbstractArray{<:DependentParameterRef};
                 [label::Type{<:AbstractSupportLabel} = PublicLabel])::Int

Return the number of support points associated with dependent infinite
parameters `prefs`. Errors if not all from the same underlying object.
Specify a subset of supports via `label` to only count the supports with `label`.
By default only the amount of public supports are given, but the full amount is
obtained via `label == All`.

**Example**
```julia-repl
julia> num_supports(x)
2
```
"""
function num_supports(prefs::AbstractArray{<:DependentParameterRef};
    label::Type{<:AbstractSupportLabel} = PublicLabel
    )::Int
    _check_complete_param_array(prefs)
    return num_supports(first(prefs), label = label)
end

"""
    has_supports(pref::DependentParameterRef)::Bool

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

**Example**
```julia-repl
julia> has_supports(x[1])
true
```
"""
has_supports(pref::DependentParameterRef)::Bool = !isempty(_parameter_supports(pref))

"""
    has_supports(prefs::AbstractArray{<:DependentParameterRef})::Bool

Return true if `prefs` have supports or false otherwise. Errors if not all of the
infinite dependent parameters are from the same object.

**Example**
```julia-repl
julia> has_supports(x)
true
```
"""
function has_supports(prefs::AbstractArray{<:DependentParameterRef})::Bool
    _check_complete_param_array(prefs)
    return has_supports(first(prefs))
end

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

Return the support points associated with `pref`. A subset of supports can be
returned via `label` to return just the supports associated with `label`. By
default only the public supports are given, but the full set is
obtained via `label == All`.

**Example**
```julia-repl
julia> supports(x[1])
2-element Array{Float64,1}:
 0.0
 1.0
```
"""
function supports(pref::DependentParameterRef;
                  label::Type{<:AbstractSupportLabel} = PublicLabel)::Vector{Float64}
    supp_dict = _parameter_supports(pref)
    pindex = _param_index(pref)
    if label == All || (!has_internal_supports(pref) && label == PublicLabel)
        return Float64[supp[pindex] for supp in keys(supp_dict)]
    else
        reduced_supps = findall(e -> any(v -> v <: label, e), supp_dict)
        return Float64[supp[pindex] for supp in reduced_supps]
    end
end

"""
    supports(prefs::AbstractArray{<:DependentParameterRef};
             [label::Type{<:AbstractSupportLabel} = PublicLabel]
             )::Union{AbstractArray{<:Vector{<:Float64}}, Array{Float64, 2}}

Return the support points associated with `prefs`. Errors if not all of the
infinite dependent parameters are from the same object. This will return a
matrix if `prefs` is `Vector`, otherwise an array of vectors is returned by
calling `supports.(prefs)`. A subset of supports can be
returned via `label` to return just the supports associated with `label`. By
default only the public supports are given, but the full set is  obtained via
`label == All`.

**Example**
```julia-repl
julia> supports(x) # columns are supports
2×2 Array{Float64,2}:
 0.0  1.0
 0.0  1.0
```
"""
function supports(prefs::AbstractArray{<:DependentParameterRef};
                  label::Type{<:AbstractSupportLabel} = PublicLabel
                  )::AbstractArray{<:Vector{<:Float64}}
    _check_complete_param_array(prefs)
    return supports.(prefs, label = label) # TODO make more efficient
end

# More efficient dispatch for Vectors
function supports(prefs::Vector{DependentParameterRef};
                  label::Type{<:AbstractSupportLabel} = PublicLabel
                  )::Array{Float64, 2}
    if !has_supports(prefs)
        return zeros(Float64, _num_parameters(first(prefs)), 0)
    elseif label == All || (!has_internal_supports(first(prefs)) && label == PublicLabel)
        raw_supps = keys(_parameter_supports(first(prefs)))
        if length(raw_supps) == 1
            return reduce(hcat, collect(raw_supps))
        else
            return reduce(hcat, raw_supps)
        end
    else
        raw_supps = findall(e -> any(v -> v <: label, e),
                            _parameter_supports(first(prefs)))
        if isempty(raw_supps)
            return zeros(Float64, _num_parameters(first(prefs)), 0)
        else
            return reduce(hcat, raw_supps)
        end
    end
end

# Define method for overriding the current supports
function _update_parameter_supports(prefs::AbstractArray{<:DependentParameterRef},
                                    supports::Array{<:Real, 2},
                                    label::Type{<:AbstractSupportLabel})::Nothing
    set = _parameter_set(first(prefs))
    new_supps = Dict{Vector{Float64}, Set{DataType}}(@views supports[:, i] =>
                                      Set([label]) for i in 1:size(supports, 2))
    sig_figs = significant_digits(first(prefs))
    methods = _derivative_methods(first(prefs))
    new_params = DependentParameters(set, new_supps, sig_figs, methods)
    _set_core_variable_object(first(prefs), new_params)
    _set_has_internal_supports(first(prefs), label <: InternalLabel)
    for pref in prefs
        _reset_derivative_evaluations(pref)
    end
    if any(is_used(pref) for pref in prefs)
        set_optimizer_model_ready(JuMP.owner_model(first(prefs)), false)
    end
    return
end

# Process an array of vectors into a support matrix
function _make_support_matrix(prefs::AbstractArray{<:DependentParameterRef},
                              supports::AbstractArray{<:Vector{<:Real}}
                              )::Array{<:Real, 2}
    _keys(supports) == _keys(prefs) || error("Inconsistent support indices")
    lens = [length(supp) for supp in supports]
    _allequal(lens) || error("Inconsistent support dimensions.")
    trans_supps = Array{Float64}(undef, first(lens), length(prefs))
    for k in eachindex(prefs)
        trans_supps[:, _param_index(prefs[k])] = supports[k]
    end
    return permutedims(trans_supps)
end

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

Specify the support points for `prefs`. Errors if the supports violate the domain
of the infinite set, if the dimensions don't match up properly,
if `prefs` and `supports` have different indices, not all of the `prefs` are
from the same dependent infinite parameter container, there are existing
supports and `force = false`. Note that it is strongly preferred to use
`add_supports` if possible to avoid destroying measure dependencies.

```julia
    set_supports(prefs::Vector{DependentParameterRef},
                 supports::Array{<:Real, 2};
                 [force::Bool = false,
                 label::Type{<:AbstractSupportLabel} = UserDefined])::Nothing
```
Specify the supports for a vector `prefs` of dependent infinite parameters.
Here rows of `supports` correspond to `prefs` and the columns correspond to the
supports. This is more efficient than the above method and will error for the
same reasons.

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

julia> set_supports(x, [0 1; 0 1])

julia> supports(x)
2×2 Array{Float64,2}:
 0.0  1.0
 0.0  1.0
```
"""
function set_supports(prefs::AbstractArray{<:DependentParameterRef},
                      supports::AbstractArray{<:Vector{<:Real}};
                      force::Bool = false,
                      label::Type{<:AbstractSupportLabel} = UserDefined
                      )::Nothing
    supps = _make_support_matrix(prefs, supports)
    set_supports(_make_vector(prefs), supps, force = force, label = label)
    return
end

# Efficient method for vector prefs and matrix of supports
function set_supports(prefs::Vector{DependentParameterRef},
                      supports::Array{<:Real, 2};
                      force::Bool = false,
                      label::Type{<:AbstractSupportLabel} = UserDefined
                      )::Nothing
    set = infinite_set(prefs) # this does a check on prefs
    if has_supports(prefs) && !force
        error("Unable set supports for $prefs since they already have supports." *
              " Consider using `add_supports` or use set `force = true` to " *
              "overwrite the existing supports.")
    elseif !supports_in_set(supports, set)
        error("Supports violate the domain of the infinite set.")
    end
    supports = round.(supports, sigdigits = significant_digits(first(prefs)))
    _update_parameter_supports(prefs, supports, label)
    return
end

# Error for single dependent parameters
function set_supports(pref::DependentParameterRef, supports; kwargs...)
    error("Cannot modify the supports of a single dependent infinite parameter.")
end

"""
    add_supports(prefs::AbstractArray{<:DependentParameterRef},
                 supports::AbstractArray{<:Vector{<:Real}};
                 [label::Type{<:AbstractSupportLabel} = UserDefined])::Nothing

Add additional support points for `prefs`. Errors if the supports violate the domain
of the infinite set, if the dimensions don't match up properly,
if `prefs` and `supports` have different indices, or not all of the `prefs` are
from the same dependent infinite parameter container.

```julia
    add_supports(prefs::Vector{DependentParameterRef},
                 supports::Array{<:Real, 2};
                 [label::Type{<:AbstractSupportLabel} = UserDefined])::Nothing
```
Specify the supports for a vector `prefs` of dependent infinite parameters.
Here rows of `supports` correspond to `prefs` and the columns correspond to the
supports. This is more efficient than the above method and will error for the
same reasons.

**Example**
```julia-repl
julia> add_supports(x, [[1], [1]])

julia> supports(x)
2×2 Array{Float64,2}:
 0.0  1.0
 0.0  1.0

julia> add_supports(x, ones(2, 1) * 0.5)

julia> supports(t)
2×3 Array{Float64,2}:
 0.0  1.0  0.5
 0.0  1.0  0.5
```
"""
function add_supports(prefs::AbstractArray{<:DependentParameterRef},
                      supports::AbstractArray{<:Vector{<:Real}};
                      label::Type{<:AbstractSupportLabel} = UserDefined, # interal keyword args
                      check::Bool = true)::Nothing
    supps = _make_support_matrix(prefs, supports)
    add_supports(_make_vector(prefs), supps, label = label, check = check)
    return
end

# More efficient version for supports in the correct format
function add_supports(prefs::Vector{DependentParameterRef},
                      supports::Array{<:Real, 2};
                      label::Type{<:AbstractSupportLabel} = UserDefined, # internal keyword args
                      check::Bool = true)::Nothing
    set = infinite_set(prefs) # this does a check on prefs
    if check && !supports_in_set(supports, set)
        error("Supports violate the domain of the infinite set.")
    end
    supports = round.(supports, sigdigits = significant_digits(first(prefs)))
    current_supports = _parameter_supports(first(prefs))
    for i in 1:size(supports, 2)
        s = @view(supports[:, i])
        if haskey(current_supports, s)
            push!(current_supports[s], label)
        else
            current_supports[s] = Set([label])
        end
    end
    if label <: InternalLabel
        _set_has_internal_supports(first(prefs), true)
    end
    for pref in prefs
        _reset_derivative_evaluations(pref)
    end
    if any(is_used(pref) for pref in prefs)
        set_optimizer_model_ready(JuMP.owner_model(first(prefs)), false)
    end
    return
end

# Error for single dependent parameters
function add_supports(pref::DependentParameterRef, supports; kwargs...)
    error("Cannot modify the supports of a single dependent infinite parameter.")
end

"""
    delete_supports(prefs::AbstractArray{<:DependentParameterRef};
                    [label::Type{<:AbstractSupportLabel} = All])::Nothing

Delete the support points for `prefs`. Errors if any of the parameters are
used by a measure or if not all belong to the same set of dependent parameters.
If `label != All` then that label is removed along with any supports that solely
contain that label.

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

```
"""
function delete_supports(prefs::AbstractArray{<:DependentParameterRef};
                         label::Type{<:AbstractSupportLabel} = All)::Nothing
    _check_complete_param_array(prefs)
    supp_dict = _parameter_supports(first(prefs))
    for pref in prefs
        _reset_derivative_evaluations(pref)
    end
    if label == All
        if any(used_by_measure(pref) for pref in prefs)
            error("Cannot delete supports with measure dependencies.")
        end
        empty!(supp_dict)
        _set_has_internal_supports(first(prefs), false)
    else
        filter!(p -> !all(v -> v <: label, p[2]), supp_dict)
        for (k, v) in supp_dict
            filter!(l -> !(l <: label), v)
        end
        pref1 = first(prefs)
        if has_internal_supports(pref1) && num_supports(pref1, label = InternalLabel) == 0
            _set_has_internal_supports(pref1, false)
        end
    end
    if any(is_used(pref) for pref in prefs)
        set_optimizer_model_ready(JuMP.owner_model(first(prefs)), false)
    end
    return
end

# Error for single dependent parameters
function delete_supports(pref::DependentParameterRef)
    error("Cannot delete the supports of a single dependent infinite parameter.")
end

# TODO resolve case that there are existing UniformGrid supports
"""
    generate_and_add_supports!(prefs::AbstractArray{<:DependentParameterRef},
                               set::InfiniteArraySet,
                               [method::Type{<:AbstractSupportLabel}];
                               [num_supports::Int = DefaultNumSupports])::Nothing

Generate supports for `prefs` 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!(::AbstractArray{<:DependentParameterRef})).
Most extensions that employ user-defined infinite sets can typically enable this
by extending [`generate_support_values`](@ref). However, in some cases it may be
necessary to extend this when more complex operations need to take place then just
adding supports to a set of infinite parameters. Errors if the
infinite set type is not recognized.
"""
function generate_and_add_supports!(prefs::AbstractArray{<:DependentParameterRef},
                                    set::InfiniteArraySet;
                                    num_supports::Int = DefaultNumSupports)::Nothing
    new_supps, label = generate_supports(set,
                                         num_supports = num_supports,
                                    sig_digits = significant_digits(first(prefs)))
    add_supports(_make_vector(prefs), new_supps, check = false, label = label)
    return
end

# Method dispatch
function generate_and_add_supports!(prefs::AbstractArray{<:DependentParameterRef},
                                    set::InfiniteArraySet,
                                    method::Type{<:AbstractSupportLabel};
                                    num_supports::Int = DefaultNumSupports)::Nothing
    new_supps, label = generate_supports(set, method,
                                         num_supports = num_supports,
                                    sig_digits = significant_digits(first(prefs)))
    add_supports(_make_vector(prefs), new_supps, check = false, label = label)
    return
end

"""
    fill_in_supports!(prefs::AbstractArray{<:DependentParameterRef};
                      [num_supports::Int = DefaultNumSupports,
                       modify::Bool = true])::Nothing

Automatically generate support points for a container of dependent infinite
parameters `prefs`. Generating up to `num_supports` for the parameters in accordance
with `generate_and_add_supports!`. 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)
2×4 Array{Float64,2}:
 0.0  0.333  0.667  1.0
 0.0  0.333  0.667  1.0
```
"""
function fill_in_supports!(prefs::AbstractArray{<:DependentParameterRef};
                           num_supports::Int = DefaultNumSupports,
                           modify::Bool = true)::Nothing
    set = infinite_set(prefs) # does check for bad container
    current_amount = InfiniteOpt.num_supports(first(prefs))
    if (modify || current_amount == 0) && (current_amount < num_supports)
        generate_and_add_supports!(prefs, set,
                                   num_supports = num_supports - current_amount)
    end
    return
end

# Error for single dependent parameters
function fill_in_supports!(pref::DependentParameterRef; kwargs...)
    error("Cannot modify the supports of a single dependent infinite parameter.")
end

"""
    fill_in_supports!(model::InfiniteModel; [num_supports::Int = DefaultNumSupports,
                      modify::Bool = true])::Nothing

Automatically generate support points for all infinite parameters in model.
This calls `fill_in_supports!` for each parameter in the model.
See [`fill_in_supports!`](@ref)
for more information. Errors if one of the infinite set types is unrecognized.
Note that no supports will be added to a particular parameter if it already has
some and `modify = false`.

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

julia> supports(t)
4-element Array{Float64,1}:
 0.0
 0.333
 0.667
 1.0
```
"""
function fill_in_supports!(model::InfiniteModel; num_supports::Int = DefaultNumSupports,
                           modify::Bool = true)::Nothing
    # fill in the the supports of each independent parameter
    for (key, data_object) in model.independent_params
        pref = dispatch_variable_ref(model, key)
        fill_in_supports!(pref, num_supports = num_supports, modify = modify)
    end
    # fill in the supports of each dependent parameter set
    for (key, data_object) in model.dependent_params
        prefs = [dispatch_variable_ref(model, DependentParameterIndex(key, i))
                 for i in 1:length(data_object.names)]
        fill_in_supports!(prefs, num_supports = num_supports, modify = modify)
    end
    return
end

################################################################################
#                          MODEL PARAMETER QUERIES
################################################################################
"""
    num_parameters(model::InfiniteModel,
                   [type::Type{InfOptParameter} = InfOptParameter])::Int

Return the number of `InfiniteOpt` parameters assigned to `model`. By default,
the total number of infinite and finite parameters is returned. The amount
of a particular type is obtained by specifying the concrete parameter type
of [`InfOptParameter`](@ref) via `type`. Type options include:
 - `InfOptParameter`: all parameters
 - `ScalarParameter`: all scalar parameters
 - `InfiniteParameter`: all infinite parameters
 - `FiniteParameter`: all finite parameters
 - `IndependentParameter`: all independent infinite parameters
 - `DependentParameters`: all dependent infinite parameters

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

julia> num_parameters(model, IndependentParameter)
2
```
"""
function num_parameters(model::InfiniteModel,
                        type::Type{InfOptParameter} = InfOptParameter
                        )::Int
    num_pars = num_parameters(model, IndependentParameter)
    num_pars += num_parameters(model, FiniteParameter)
    num_pars += num_parameters(model, DependentParameters)
    return num_pars
end

# Particular scalar parameter types
function num_parameters(model::InfiniteModel,
                        type::Type{C})::Int where {C <: ScalarParameter}
    return length(_data_dictionary(model, type))
end

# ScalarParameter
function num_parameters(model::InfiniteModel,
                        type::Type{ScalarParameter})::Int
    num_pars = num_parameters(model, FiniteParameter)
    num_pars += num_parameters(model, IndependentParameter)
    return num_pars
end

# DependentParameters
function num_parameters(model::InfiniteModel,
                        type::Type{DependentParameters})::Int
    num_pars = 0
    for (_, object) in _data_dictionary(model, type)
        num_pars += length(object.names)
    end
    return num_pars
end

# InfiniteParameter
function num_parameters(model::InfiniteModel,
                        type::Type{InfiniteParameter})::Int
    num_pars = num_parameters(model, IndependentParameter)
    num_pars += num_parameters(model, DependentParameters)
    return num_pars
end

"""
    all_parameters(model::InfiniteModel,
                   type::Type{InfOptParameter} = InfOptParameter
                   )::Vector{GeneralVariableRef}

Return a list of all the `InfiniteOpt` parameters assigned to `model`. By default,
all of the infinite and finite parameters is returned. The search is reduced to
a particular type is obtained by specifying the concrete parameter type
of [`InfOptParameter`](@ref) via `type`. Type options include:
- `InfOptParameter`: all parameters
- `ScalarParameter`: all scalar parameters
- `InfiniteParameter`: all infinite parameters
- `FiniteParameter`: all finite parameters
- `IndependentParameter`: all independent infinite parameters
- `DependentParameters`: all dependent infinite parameters

**Examples**
```julia-repl
julia> all_parameters(model)
4-element Array{GeneralVariableRef,1}:
 t
 x[1]
 x[2]
 alpha

julia> all_parameters(model, FiniteParameter)
1-element Array{GeneralVariableRef,1}:
 alpha
```
"""
function all_parameters(model::InfiniteModel,
                        type::Type{InfOptParameter} = InfOptParameter
                        )::Vector{GeneralVariableRef}
    prefs_list = all_parameters(model, IndependentParameter)
    append!(prefs_list, all_parameters(model, DependentParameters))
    append!(prefs_list, all_parameters(model, FiniteParameter))
    return prefs_list
end

# Particular scalar parameter types
function all_parameters(model::InfiniteModel,
                        type::Type{C}
                        )::Vector{GeneralVariableRef} where {C <: InfOptParameter}
    prefs_list = Vector{GeneralVariableRef}(undef, num_parameters(model, type))
    for (i, (index, _)) in enumerate(_data_dictionary(model, type))
        prefs_list[i] = _make_parameter_ref(model, index)
    end
    return prefs_list
end

# ScalarParameter
function all_parameters(model::InfiniteModel,
                        type::Type{ScalarParameter})::Vector{GeneralVariableRef}
    prefs_list = all_parameters(model, IndependentParameter)
    append!(prefs_list, all_parameters(model, FiniteParameter))
    return prefs_list
end

# DependentParameters
function all_parameters(model::InfiniteModel,
                        type::Type{DependentParameters}
                        )::Vector{GeneralVariableRef}
    prefs_list = Vector{GeneralVariableRef}(undef, num_parameters(model, type))
    counter = 1
    for (index, object) in _data_dictionary(model, type)
        for i in eachindex(object.names)
            dep_idx = DependentParameterIndex(index, i)
            prefs_list[counter] = _make_parameter_ref(model, dep_idx)
            counter += 1
        end
    end
    return prefs_list
end

# InfiniteParameter
function all_parameters(model::InfiniteModel,
                        type::Type{InfiniteParameter}
                        )::Vector{GeneralVariableRef}
    prefs_list = all_parameters(model, DependentParameters)
    append!(prefs_list, all_parameters(model, IndependentParameter))
    return prefs_list
end

################################################################################
#                                 DELETION
################################################################################
"""
    JuMP.delete(model::InfiniteModel,
                prefs::AbstractArray{<:DependentParameterRef})::Nothing

Extend `JuMP.delete` to delete
dependent infinite parameters and their dependencies. All variables, constraints, and
measure functions that depend on `prefs` are updated to exclude them. Errors if the
parameters are 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 to the default 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, 1], x[2] ∈ [-1, 1]
 g(0.5, x) = 0, x[1] ∈ [-1, 1], x[2] ∈ [-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,
                     prefs::AbstractArray{<:DependentParameterRef})::Nothing
    @assert JuMP.is_valid(model, first(prefs)) "Parameter references are invalid."
    _check_complete_param_array(prefs)
    gvrefs = [_make_parameter_ref(model, JuMP.index(pref)) for pref in prefs]
    # ensure deletion is okay (prefs are not used by measure data)
    for pref in gvrefs
        for mindex in _measure_dependencies(pref)
            data = measure_data(dispatch_variable_ref(model, mindex))
            _check_param_in_data(pref, data)
        end
    end
    # make sure it isn't used by parameter function
    if used_by_parameter_function(first(prefs))
        error("Cannot delete `$prefs` since they are used by an infinite " *
              "parameter function.")
    end
    # update optimizer model status
    if any(is_used(pref) for pref in prefs)
        set_optimizer_model_ready(model, false)
    end
    # delete dependence of measures and constraints on prefs
    for pref in gvrefs
        _update_measures(model, pref)
        _update_constraints(model, pref)
    end
    # get the object and parameter numbers
    obj_num = _object_number(first(prefs))
    param_nums = collect(_data_object(first(prefs)).parameter_nums)
    # update infinite variables that depend on pref
    for vindex in _infinite_variable_dependencies(first(prefs))
        # remove the parameter dependences
        vref = InfiniteVariableRef(model, vindex)
        vprefs = raw_parameter_refs(vref)
        tup_index = findfirst(isequal(obj_num), _object_number.(vprefs[:, 1]))
        delete_indices = vprefs.ranges[tup_index]
        deleteat!(vprefs, tup_index, tuple_index = true)
        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_indices)
        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_indices)
        end
    end
    # delete derivatives that depend on any of these parameters
    for pref in gvrefs
        for index in _derivative_dependencies(pref)
            JuMP.delete(model, dispatch_variable_ref(model, index))
        end
    end
    # delete parameter information stored in model
    _delete_data_object(first(prefs))
    # update the object numbers and parameter numbers
    _update_model_numbers(model, obj_num, param_nums)
    return
end


text:
module TestSpacePrint

using Test
using TreeParzen

# Print out a space for visual inspection

space = Dict(
    :args => HP.PChoice(:probchoice, [
        Prob(
            0.1,
            Dict(
                :symbol => :test,
                :string => "test",
                :int => 1,
                :float => 1.0,
                :nothing => nothing,
            )
        ),
        Prob(
            0.2,
            Dict(
                :dists => HP.Choice(:dists, [
                    HP.Uniform(:uniform, -10.0, 10.0) + HP.Normal(:normal, -10.0, 10.0),
                    float(HP.LogNormal(:lognormal, -10.0, 10.0)),
                ]),
            )
        ),
        Prob(
            0.7,
            Dict(
                :dictmiddle => Dict(
                    :tuple => (0.1, 0.3, 0.6),
                    :set => Set([4, 5, 6]),
                    :vector => [7, 8, 9],
                ),
                :dictend => Dict(:a => 1, :b => 4),
            )
        ),
    ]),
)
spaceprint(space)

end # module TestSpacePrint
true

ext: jl


text:
using Documenter, IntervalMatrices

DocMeta.setdocmeta!(IntervalMatrices, :DocTestSetup, :(using IntervalMatrices);
                    recursive=true)

makedocs(
    sitename = "IntervalMatrices.jl",
    modules = [IntervalMatrices],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets = ["assets/juliareach.css"]),
    pages = [
        "Home" => "index.md",
        "Library" => Any[
        "Types" => "lib/types.md",
        "Methods" => "lib/methods.md"],
        "About" => "about.md"
    ],
    strict = true
)

deploydocs(
    repo = "github.com/JuliaReach/IntervalMatrices.jl.git"
)


text:
repo = "bsuite"
export bsuite_params

const BSUITE_DATASETS = [
    "cartpole",
    "catch",
    "mountain_car"
]

types = [
    "full",
    "full_train",
    "full_valid"
]

function bsuite_params()
    game = BSUITE_DATASETS
    shards = 0:4
    type = types

    @info game shards type
end

function bsuite_init()
    for env in BSUITE_DATASETS
        for stochasticity in 0.0:0.1:0.5
            for index in 0:4
                for type in types
                    register(
                        DataDep(
                            "rl-unplugged-$repo-$env-$stochasticity-$index-$type",
                            """
                            Dataset: RL Unplugged bsuite
                            Credits: https://arxiv.org/abs/2006.13888
                            Url: https://github.com/deepmind/deepmind-research/tree/master/rl_unplugged
                            Authors: Caglar Gulcehre, Ziyu Wang, Alexander Novikov, Tom Le Paine,
                            Sergio Gómez Colmenarejo, Konrad Zolna, Rishabh Agarwal,
                            Josh Merel, Daniel Mankowitz, Cosmin Paduraru, Gabriel
                            Dulac-Arnold, Jerry Li, Mohammad Norouzi, Matt Hoffman,
                            Ofir Nachum, George Tucker, Nicolas Heess, Nando deFreitas
                            Year: 2020

                            Data accompanying [RL Unplugged: Benchmarks for Offline Reinforcement Learning].
                            bsuite datasets are fairly light-weight and running experiments doesn't require
                            too much compute. We recommend you to try bsuite, if you are interested in
                            small-scale and easy to run offline RL datasets generated by stochastic environments
                            where the stochasticity of the environment is easy to control.
                            """,
                            "gs://rl_unplugged/$repo/$env/0_$stochasticity/$(index)_$type-00000-of-00001",
                            fetch_method = fetch_gc_file
                        )
                    )
                end
            end
        end
    end
end


text:
function plot_lines!(vis::Visualizer, model::Hopper3D, q::AbstractVector;
		r_foot=0.04, offset=0.04, size=10, name::Symbol=:hopper_3D, col::Bool=true)
	p_shift = [0.0, 0.0, r_foot]
	p_offset = [0.0, -offset, 0.0]

	# Point Traj
	top_point = Vector{Point{3,Float64}}()
	bot_point = Vector{Point{3,Float64}}()

	for qi in q
		push!(top_point, Point(qi[1:3]...) + p_shift + p_offset)
		push!(bot_point, Point(kinematics(model, qi)...) + p_shift + p_offset)
	end

	# Set lines
	orange_mat, blue_mat, black_mat = get_line_material(size)
	if col
		setobject!(vis[name]["/lines/top"], MeshCat.Line(top_point, orange_mat))
		setobject!(vis[name]["/lines/bot"], MeshCat.Line(bot_point, blue_mat))
	else
		setobject!(vis[name]["/lines/top"], MeshCat.Line(top_point, black_mat))
		setobject!(vis[name]["/lines/bot"], MeshCat.Line(bot_point, black_mat))
	end
	return nothing
end

function build_robot!(vis::Visualizer, model::Hopper3D; name::Symbol=:Hopper3D, r=0.04, α=1.0)
	n_leg = 100
	r = convert(Float32, r)
	r_leg = convert(Float32, 0.5 * r)
	body_mat = MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, α))
	contact_mat = MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0.0, α))

	default_background!(vis)

    setobject!(vis[name][:robot]["body"], Sphere(Point3f0(0),
        convert(Float32, 0.1)), body_mat)
    setobject!(vis[name][:robot]["foot"], Sphere(Point3f0(0), r), contact_mat)
    for i = 1:n_leg
        setobject!(vis[name][:robot]["leg$i"], Sphere(Point3f0(0), r_leg), body_mat)
    end
	return nothing
end

function set_robot!(vis::Visualizer, model::Hopper3D, q::AbstractVector;
		name::Symbol=:Hopper3D, r=0.04)
	n_leg = 100
	r = convert(Float32, r)
	r_leg = convert(Float32, 0.5 * r)
	p_shift = [0.0; 0.0; r]

    p_body = q[1:3]
    p_foot = kinematics(model, q)

    settransform!(vis[name][:robot]["body"], Translation(p_body + p_shift))
    settransform!(vis[name][:robot]["foot"], Translation(p_foot + p_shift))
	r_range = range(0, stop = q[7], length = n_leg)
    for i = 1:n_leg
		q_tmp = Array(copy(q))
		q_tmp[7] = r_range[i]
		p_leg = kinematics(model, q_tmp)
        settransform!(vis[name][:robot]["leg$i"], Translation(p_leg + p_shift))
    end
	return nothing
end

function contact_point(model::Hopper3D, q::AbstractVector)
	p_foot = kinematics(model, q)
	pc = [p_foot]
	return pc
end


text:
##########
##########
# Functions to perform perturbations on the food webs
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2021-11-20
##########
##########

# set up =======================================================================
using DrWatson
@quickactivate "trait-based-rewiring"
# logger set up
using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())

include(srcdir("01_eqn_component_functions.jl"))
include(srcdir("00_species_list_creation.jl"))
include(srcdir("04_perturbation_sims_functions.jl"))

# set all positional values
organism_id_pos = 1
organism_type_pos = 2
body_size_pos = 3
lower_prey_limit_pos = 4
upper_prey_limit_pos = 5
habitat_midpoint_pos = 6
habitat_lower_limit_pos = 7
habitat_upper_limit_pos = 8
nutritional_value_pos = 9
biomass_pos = 10
network_id_pos = 11
alive_pos = 12

# periodic pulse perturbations =================================================

# little line to set the network id (this happens from the shell line)
network_id = parse(Int64, ARGS[1])

species_list = readdlm(datadir("./intermediate-files/species-lists-pre-perturbation/species_list_0001.csv"), ',', Any)
network = create_interaction_matrix(species_list)

# set vector callback
extinction_callback =
    VectorContinuousCallback(extinction_threshold,
                            extinction_affect!,
                            size(species_list, 1),
                            save_positions = ntuple(x->true,
                                                    size(species_list, 1)))
steady_callback = TerminateSteadyState(1e-7,
                                        1e-5,
                                        DiffEqCallbacks.allDerivPass,
                                        min_t = 1205)

# set callback set
save_positions = (true, true)
cb_inc = DiscreteCallback(inc_mort_cond, affect_inc!, save_positions=save_positions)
save_positions = (false, true)
cb_dec = DiscreteCallback(dec_mort_cond, affect_dec!, save_positions=save_positions)
callback_set = CallbackSet(extinction_callback,
                            #steady_callback, #keep this out in the period ones
                            cb_inc,
                            cb_dec)

#  allocate x0
u0 = Float64[]
for i in 1:size(species_list, 1)
    append!(u0, species_list[i, biomass_pos])
end

# set perturbation tspan
tspan = (1000.0, 5000.0)

# set constant stop values
const stop1 = [1000.0, 2000.0, 3000.0, 4000.0]
const stop2 = [1100.0, 2100.0, 3100.0, 4100.0]
const tstop = [1000.;1100.;2000.;2100;3000.;3100.;4000.;4100.]

c = 210
x = 100
f =
z = c/(x*(1/f))

# loop through random selections of species to perturb
#@info{"starting loop..."}

for j in 1:10

    # pick set of random species to perturb
    all = Array(1:size(species_list, 1))
    perturb_sp = rand(Array(1:size(species_list, 1)),
                        Int(round(size(species_list, 1)*(1/3))))
    non_perturb = setdiff(Set(all), Set(perturb_sp))

    # set times to save from the simulation
    save_times = [collect(range(1000.0, 1250.0, step = 0.1));
                collect(range(1250, 1999, step = sample(200:400, 1)[1]));
                collect(range(2000.0, 2250.0, step = 0.1));
                collect(range(2250, 2999, step = sample(200:400, 1)[1]));
                collect(range(3000.0, 3250.0, step = 0.1));
                collect(range(3250, 3999, step = sample(200:400, 1)[1]));
                collect(range(4000.0, 4250.0, step = 0.1));
                collect(range(4250, 5000, step = sample(200:400, 1)[1]));]

    # set p
    p = [species_list, network, 0.3, z, perturb_sp, non_perturb]

    prob_period = ODEProblem(biomass_perturb_random,
                                        u0,
                                        tspan,
                                        p,
                                        tstops=tstop)
    @time sol_prob_period = solve(prob_period,
                                        Tsit5(),
                                        callback = callback_set,
                                        # save every part for the perturbation and
                                        # the first 150 time steps after the
                                        # perturbation, but then just grab a few
                                        # from there to the end
                                        saveat = save_times)
    # save output of the solver
    output = hcat(sol_prob_period.t, sol_prob_period.u)

    # write out the biomass outputs
    name = string("periodic-", network_id, "-", j, "-", "sims.csv")
    path = "./data/simulation-outputs/period-sims/"
    writedlm(string(path,name), output, ",")
end

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
# eventually these should go inside the module or their own modules and I should choose what to export. not yet
include("LJH.jl")
using LJH
using Dierckx, PyCall, HDF5, ArrayViews, Distributions, KernelDensity, PyPlot, JLD
import Graphs

include("summarize.jl")
include("mockpulses.jl")
include("runningvectors.jl")
include("histogram.jl")
include("steps.jl")
include("steps_macros.jl")
include("calibration.jl")
include("corrections.jl")
include("ljhutil.jl")
include("tes_optimal_filtering.jl")
include("glue.jl")
include("matterwatcher.jl")
module Mass2

# package code goes here

end # module


text:
using Test

@testset "generated_weights" begin
    include("generated_weights/mfc+mfc+softmax.jl")
    include("generated_weights/conv+fc+softmax.jl")
end


text:
"""
    AbstractDiffMethod

Abstract type representing different differentiation methods
"""
abstract type AbstractDiffMethod end

# --- differentiation method --- #

"""
    AbstractAD

Abstract type representing different automatic differentation methods
"""
abstract type AbstractAD <: AbstractDiffMethod end

"""
    ForwardAD{T} <: AbstractAD

Forward automatic differentiation using the ForwardDiff package.

# Fields:
 - `chunk::T`: Chunk size
"""
struct ForwardAD{T} <: AbstractAD
    chunk::T
end

# use default chunk size if chunk is not specified
ForwardAD() = ForwardAD(nothing)

"""
    ReverseAD <: AbstractAD

Reverse automatic differentiation using the ReverseDiff package.

"""
struct ReverseAD <: AbstractAD end

"""
    AbstractFD

Abstract type representing different finite difference methods
"""
abstract type AbstractFD <: AbstractDiffMethod end

"""
    ForwardFD <: AbstractFD

Forward finite differencing as implemented by the FiniteDiff package
"""
struct ForwardFD <: AbstractFD end

"""
    CentralFD <: AbstractFD

Central finite differencing as implemented by the FiniteDiff package
"""
struct CentralFD <: AbstractFD end

"""
    ComplexFD <: AbstractFD

Complex step finite differencing as implemented by the FiniteDiff package
"""
struct ComplexFD <: AbstractFD end

"""
    AbstractChainRuleMode

Type representing direction in which to to apply the chain rule when computing
the derivatives of an explicit system.
"""
abstract type AbstractChainRuleMode end

"""
    Forward <: AbstractChainRuleMode

Type indicating that the chain rule should be applied in forward mode
"""
struct Forward <: AbstractChainRuleMode end

"""
    Reverse <: AbstractChainRuleMode

Type indicating that the chain rule should be applied in reverse mode
"""
struct Reverse <: AbstractChainRuleMode end

"""
    AbstractSensitivityMode

Abstract type representing mode by which to apply the analytic sensitivity
equations when computing derivatives of an explicit component created by
combining an implicit component, output component, and solver.
"""
abstract type AbstractSensitivityMode end

"""
    Direct <: AbstractSensitivityMode

Type indicating that the analytic sensitivity equations should be applied using
the direct method.
"""
struct Direct <: AbstractSensitivityMode end

"""
    Adjoint <: AbstractSensitivityMode

Type indicating that the analytic sensitivity equations should be applied using
the adjoint method.
"""
struct Adjoint <: AbstractSensitivityMode end

"""
    finitediff_type(mode::AbstractFD)

Converts the types used by this package to the inputs expected by the FiniteDiff
package.
"""
finitediff_type
finitediff_type(::ForwardFD) = Val(:forward)
finitediff_type(::CentralFD) = Val(:central)
finitediff_type(::ComplexFD) = Val(:complex)

# --- sparsity patterns --- #

"""
    AbstractSparsityPattern

Abstract type representing the sparsity pattern of a jacobian matrix
"""
abstract type AbstractSparsityPattern end

"""
    DensePattern <: AbstractSparsityPattern

Type which indicates that a jacobian is dense
"""
struct DensePattern <: AbstractSparsityPattern end

"""
    SparsePattern{TI} <: AbstractSparsityPattern

Type which indicates that a jacobian is sparse and indicates the non-zero rows
and columns

# Fields:
 - `rows::Vector{TI}`: Row of each non-zero matrix element
 - `cols::Vector{TI}`: Column of each non-zero matrix element
"""
struct SparsePattern{TI} <: AbstractSparsityPattern
    rows::Vector{TI}
    cols::Vector{TI}
end

# --- jacobian function creation --- #

"""
    create_output_jacobian_functions(f!, x0, y0, dtype, sp)

Creates functions to compute the jacobian using the specified derivative
type and sparsity pattern.

# Arguments
- `f!`: In-place version of the output function: `y = f(x)`
- `x0::AbstractVector`: Initial values for inputs (used for size and type information)
- `y0::AbstractVector`: Initial values for outputs (used for size and type information)
- `dtype::AbstractDiffMethod`: Method used to calculate derivatives
- `sparsity::AbstractSparsityPattern`: Jacobian sparsity structure
"""
create_output_jacobian_functions

function create_output_jacobian_functions(f!, x0, y0, dtype::ForwardAD, sp::DensePattern)

    # set chunk size
    chunk = isnothing(dtype.chunk) ? ForwardDiff.Chunk(x0) : dtype.chunk

    # configure ForwardDiff
    cfg = ForwardDiff.JacobianConfig(f!, y0, x0, chunk)

    # check tag now since we disable it when running our function
    ForwardDiff.checktag(cfg, f!, x0)

    # define jacobian function
    ycache = copy(y0)
    df! = function(dydx, x)
        ForwardDiff.jacobian!(dydx, f!, ycache, x, cfg, Val{false}())
        return dydx
    end

    # define combined output and jacobian function
    fdf! = function(y, dydx, x)
        ForwardDiff.jacobian!(dydx, f!, y, x, cfg, Val{false}())
        return y, dydx
    end

    return df!, fdf!
end

function create_output_jacobian_functions(f!, x0, y0, dtype::ReverseAD, sp::DensePattern)

    f_tape = ReverseDiff.JacobianTape(f!, y0, x0)
    cache = ReverseDiff.compile(f_tape)
    result = DiffResults.JacobianResult(y0, x0)

    # define jacobian function
    df! = function(dydx, x)
        ReverseDiff.jacobian!(result, cache, x)
        copyto!(dydx, DiffResults.jacobian(result))
        return dydx
    end

    # define combined output and jacobian function
    fdf! = function(y, dydx, x)
        ReverseDiff.jacobian!(result, cache, x)
        copyto!(y, DiffResults.value(result))
        copyto!(dydx, DiffResults.jacobian(result))
        return y, dydx
    end

    return df!, fdf!
end

function create_output_jacobian_functions(f!, x0, y0, dtype::AbstractFD, sp::DensePattern)

    # Construct FiniteDiff Cache
    xcache = copy(x0)
    ycache = copy(y0) # this needs to be updated before calling FiniteDiff
    ycache1 = copy(y0)
    fdtype = finitediff_type(dtype)
    cache = FiniteDiff.JacobianCache(xcache, ycache, ycache1, fdtype)

    # define jacobian function
    df! = function(dydx, x)
        TF = promote_type(eltype(x), eltype(ycache))
        if TF <: eltype(ycache)
            # don't allocate because we can store the result
            f!(ycache, x)
            copyto!(xcache, x)
            FiniteDiff.finite_difference_jacobian!(dydx, f!, xcache, cache, ycache)
        else
            # allocate because we can't store the result
            FiniteDiff.finite_difference_jacobian!(dydx, f!, x)
        end

        return dydx
    end

    # define combined output and jacobian function
    fdf! = function(y, dydx, x)
        f!(y, x)
        TF = promote_type(eltype(y), eltype(ycache))
        if TF <: eltype(ycache)
            # don't allocate because we can store the result
            copyto!(xcache, x)
            copyto!(ycache, y)
            FiniteDiff.finite_difference_jacobian!(dydx, f!, xcache, cache)
        else
            # allocate because we can't store the result
            FiniteDiff.finite_difference_jacobian!(dydx, f!, x)
        end
        return y, dydx
    end

    return df!, fdf!
end

function create_output_jacobian_functions(f!, x0, y0, dtype::ForwardAD, sp::SparsePattern)

    # get chunk size for ForwardDiff
    chunk = isnothing(dtype.chunk) ? ForwardDiff.Chunk(x0) : dtype.chunk

    # construct cache
    sparsity = sparse(sp.rows, sp.cols, ones(length(sp.rows)))
    colorvec = SparseDiffTools.matrix_colors(sparsity)
    cache = SparseDiffTools.ForwardColorJacCache(f!, x0, chunk; dx=y0,
        colorvec=colorvec, sparsity=sparsity)

    # define jacobian function
    df! = function(dydx, x)
        SparseDiffTools.forwarddiff_color_jacobian!(dydx, f!, x, cache)
        return dydx
    end

    # define combined output and jacobian function
    fdf! = function(y, dydx, x)
        f!(y, x)
        SparseDiffTools.forwarddiff_color_jacobian!(dydx, f!, x, cache)
        return y, dydx
    end

    return df!, fdf!
end

function create_output_jacobian_functions(f!, x0, y0, dtype::AbstractFD, sp::SparsePattern)

    # Construct FiniteDiff Cache
    xcache = copy(x0)
    ycache = copy(y0)
    ycache1 = copy(y0)
    fdtype = finitediff_type(dtype)
    sparsity = sparse(sp.rows, sp.cols, ones(length(sp.rows)))
    colorvec = SparseDiffTools.matrix_colors(sparsity)
    cache = FiniteDiff.JacobianCache(xcache, ycache, ycache1, fdtype, sparsity,
        colorvec)

    # define jacobian function
    df! = function(dydx, x)
        f!(ycache, x)
        FiniteDiff.finite_difference_jacobian!(dydx, f!, x, cache)
        return dydx
    end

    # define combined output and jacobian function
    fdf! = function(y, dydx, x)
        f!(y, x)
        copyto!(ycache, y)
        FiniteDiff.finite_difference_jacobian!(dydx, f!, x, cache)
        return y, dydx
    end

    return df!, fdf!
end

"""
    create_residual_input_jacobian_functions(f!, x0, y0, r0, dtype, sp)

Creates functions to compute the jacobian of the residual with respect to the
inputs using the specified derivative type and sparsity pattern.

# Arguments
- `f!`: In-place residual function `f(r, x, y)`.
- `x0`: Initial values for inputs (used for size and type information)
- `y0`: Initial values for outputs (used for size and type information)
- `r0`: Initial values for residuals (used for size and type information)
- `dtype`: Method to use to calculate derivatives
- `sparsity`: Jacobian sparsity structure
"""
create_residual_input_jacobian_functions

function create_residual_input_jacobian_functions(f!, x0, y0, r0, dtype::ForwardAD, sp::DensePattern)

    # construct function to take derivative of
    ycache = copy(y0)
    fx! = (r, x) -> f!(r, x, ycache)

    # set chunk size
    chunk = isnothing(dtype.chunk) ? ForwardDiff.Chunk(x0) : dtype.chunk

    # configure ForwardDiff
    cfg = ForwardDiff.JacobianConfig(fx!, r0, x0, chunk)

    # check tag now since we disable it when running our function
    ForwardDiff.checktag(cfg, fx!, x0)

    # define jacobian function
    rcache = copy(r0)
    dfdx! = function(drdx, x, y)
        copyto!(ycache, y)
        ForwardDiff.jacobian!(drdx, fx!, rcache, x, cfg, Val{false}())
        return drdx
    end

    # define combined output and jacobian function
    fdfdx! = function(r, drdx, x, y)
        copyto!(ycache, y)
        ForwardDiff.jacobian!(drdx, fx!, r, x, cfg, Val{false}())
        return r, drdx
    end

    return dfdx!, fdfdx!
end

function create_residual_input_jacobian_functions(f!, x0, y0, r0, dtype::ReverseAD, sp::DensePattern)

    f_tape = ReverseDiff.JacobianTape(f!, r0, (x0, y0))
    cache = ReverseDiff.compile(f_tape)
    xres = DiffResults.JacobianResult(r0, x0)
    yres = DiffResults.JacobianResult(r0, y0)

    # define jacobian function
    dfdx! = function(drdx, x, y)
        ReverseDiff.jacobian!((xres, yres), cache, (x, y))
        copyto!(drdx, DiffResults.jacobian(xres))
        return drdx
    end

    # define combined output and jacobian function
    fdfdx! = function(r, drdx, x, y)
        ReverseDiff.jacobian!((xres, yres), cache, (x, y))
        copyto!(r, DiffResults.value(xres))
        copyto!(drdx, DiffResults.jacobian(xres))
        return r, drdx
    end

    return dfdx!, fdfdx!
end

function create_residual_input_jacobian_functions(f!, x0, y0, r0, dtype::AbstractFD, sp::DensePattern)

    # construct function to take derivative of
    ycache = copy(y0)
    fx! = (r, x) -> f!(r, x, ycache)

    # Construct FiniteDiff Cache
    xcache = copy(x0)
    rcache = copy(r0) # this needs to be updated
    rcache1 = copy(r0)
    fdtype = finitediff_type(dtype)
    cache = FiniteDiff.JacobianCache(xcache, rcache, rcache1, fdtype)

    # define jacobian function
    dfdx! = function(drdx, x, y)
        TF = promote_type(eltype(y), eltype(ycache))
        if TF <: eltype(ycache)
            # don't allocate because we can store the result
            f!(rcache, x, y)
            copyto!(xcache, x)
            copyto!(ycache, y)
            FiniteDiff.finite_difference_jacobian!(drdx, fx!, xcache, cache)
        else
            # allocate because we can't store the result
            fx_new! = (r, x) -> f!(r, x, y)
            FiniteDiff.finite_difference_jacobian!(drdx, fx_new!, x, fdtype)
        end
        return drdx
    end

    # define combined output and jacobian function
    fdfdx! = function(r, drdx, x, y)
        TF = promote_type(eltype(y), eltype(ycache))
        if TF <: eltype(ycache)
            # don't allocate because we can store the result
            f!(r, x, y)
            copyto!(xcache, x)
            copyto!(ycache, y)
            copyto!(rcache, r)
            FiniteDiff.finite_difference_jacobian!(drdx, fx!, xcache, cache)
        else
            # allocate because we can't store the result
            fx_new! = (r, x) -> f!(r, x, y)
            f!(r, x, y)
            FiniteDiff.finite_difference_jacobian!(drdx, fx_new!, x, fdtype)
        end
        return r, drdx
    end

    return dfdx!, fdfdx!
end

function create_residual_input_jacobian_functions(f!, x0, y0, r0, dtype::ForwardAD, sp::SparsePattern)

    # construct function to take derivative of
    ycache = copy(y0)
    fx! = (r, x) -> f!(r, x, ycache)

    # get chunk size for ForwardDiff
    chunk = isnothing(dtype.chunk) ? ForwardDiff.Chunk(x0) : dtype.chunk

    # construct cache
    sparsity = sparse(sp.rows, sp.cols, ones(length(sp.rows)))
    colorvec = SparseDiffTools.matrix_colors(sparsity)
    cache = SparseDiffTools.ForwardColorJacCache(fx!, x0, chunk; dx=r0,
        colorvec=colorvec, sparsity=sparsity)

    # define jacobian function
    dfdx! = function(drdx, x, y)
        copyto!(ycache, y)
        SparseDiffTools.forwarddiff_color_jacobian!(drdx, fx!, x, cache)
        return drdx
    end

    # define combined output and jacobian function
    fdfdx! = function(r, drdx, x, y)
        copyto!(ycache, y)
        f!(r, x, y)
        SparseDiffTools.forwarddiff_color_jacobian!(drdx, fx!, x, cache)
        return r, drdx
    end

    return dfdx!, fdfdx!
end

function create_residual_input_jacobian_functions(f!, x0, y0, r0, dtype::AbstractFD, sp::SparsePattern)

    # construct function to take derivative of
    ycache = copy(y0)
    fx! = (r, x) -> f!(r, x, ycache)

    # Construct FiniteDiff Cache
    xcache = copy(x0)
    rcache = copy(r0) # this needs to be updated before calling FiniteDiff
    rcache1 = copy(r0)
    fdtype = finitediff_type(dtype)
    sparsity = sparse(sp.rows, sp.cols, ones(length(sp.rows)))
    colorvec = SparseDiffTools.matrix_colors(sparsity)
    cache = FiniteDiff.JacobianCache(xcache, rcache, rcache1, fdtype, sparsity,
        colorvec)

    # define jacobian function
    dfdx! = function(drdx, x, y)
        copyto!(ycache, y)
        f!(rcache, x, y)
        FiniteDiff.finite_difference_jacobian!(drdx, fx!, x, cache)
        return drdx
    end

    # define combined output and jacobian function
    fdfdx! = function(r, drdx, x, y)
        f!(r, x, y)
        copyto!(ycache, y)
        copyto!(rcache, r)
        FiniteDiff.finite_difference_jacobian!(drdx, fx!, x, cache)
        return r, drdx
    end

    return dfdx!, fdfdx!
end

"""
    create_residual_output_jacobian_functions(f!, x0, y0, r0, dtype, sp)

Creates functions to compute the jacobian of the residual with respect to the
inputs using the specified derivative type and sparsity pattern.

# Arguments
- `f!`: In-place residual function `f(r, x, y)`.
- `x0`: Initial values for inputs (used for size and type information)
- `y0`: Initial values for outputs (used for size and type information)
- `r0`: Initial values for residuals (used for size and type information)
- `dtype`: Method to use to calculate derivatives
- `sparsity`: Jacobian sparsity structure
"""
create_residual_output_jacobian_functions

function create_residual_output_jacobian_functions(f!, x0, y0, r0, dtype::ForwardAD, sp::DensePattern)

    # construct function to take derivative of
    xcache = copy(x0)
    fy! = (r, y) -> f!(r, xcache, y)

    # set chunk size
    chunk = isnothing(dtype.chunk) ? ForwardDiff.Chunk(y0) : dtype.chunk

    # configure ForwardDiff
    cfg = ForwardDiff.JacobianConfig(fy!, r0, y0, chunk)

    # check tag now since we disable it when running our function
    ForwardDiff.checktag(cfg, fy!, y0)

    # define jacobian function
    rcache = copy(r0)
    dfdy! = function(drdy, x, y)
        copyto!(xcache, x)
        ForwardDiff.jacobian!(drdy, fy!, rcache, y, cfg, Val{false}())
        return drdy
    end

    # define combined output and jacobian function
    fdfdy! = function(r, drdy, x, y)
        copyto!(xcache, x)
        ForwardDiff.jacobian!(drdy, fy!, r, y, cfg, Val{false}())
        return r, drdy
    end

    return dfdy!, fdfdy!
end

function create_residual_output_jacobian_functions(f!, x0, y0, r0, dtype::ReverseAD, sp::DensePattern)

    f_tape = ReverseDiff.JacobianTape(f!, r0, (x0, y0))
    cache = ReverseDiff.compile(f_tape)
    xres = DiffResults.JacobianResult(r0, x0)
    yres = DiffResults.JacobianResult(r0, y0)

    # define jacobian function
    dfdy! = function(drdy, x, y)
        ReverseDiff.jacobian!((xres, yres), cache, (x, y))
        copyto!(drdy, DiffResults.jacobian(yres))
        return drdy
    end

    # define combined output and jacobian function
    fdfdy! = function(r, drdy, x, y)
        ReverseDiff.jacobian!((xres, yres), cache, (x, y))
        copyto!(r, DiffResults.value(yres))
        copyto!(drdy, DiffResults.jacobian(yres))
        return r, drdy
    end

    return dfdy!, fdfdy!
end

function create_residual_output_jacobian_functions(f!, x0, y0, r0, dtype::AbstractFD, sp::DensePattern)

    # construct function to take derivative of
    xcache = copy(x0)
    fy! = (r, y) -> f!(r, xcache, y)

    # Construct FiniteDiff Cache
    ycache = copy(y0)
    rcache = copy(r0) # this needs to be updated
    rcache1 = copy(r0)
    fdtype = finitediff_type(dtype)
    cache = FiniteDiff.JacobianCache(ycache, rcache, rcache1, fdtype)

    # define jacobian function
    dfdy! = function(drdy, x, y)
        TF = promote_type(eltype(y), eltype(ycache))
        if TF <: eltype(xcache)
            # don't allocate because we can store the result
            f!(rcache, x, y)
            copyto!(xcache, x)
            copyto!(ycache, y)
            FiniteDiff.finite_difference_jacobian!(drdy, fy!, ycache, cache)
        else
            # allocate because we can't store the result
            fy_new! = (r, y) -> f!(r, x, y)
            FiniteDiff.finite_difference_jacobian!(drdy, fy!, y, fdtype)
        end
        return drdy
    end

    # define combined output and jacobian function
    fdfdy! = function(r, drdy, x, y)
        TF = promote_type(eltype(x), eltype(xcache))
        if TF <: eltype(xcache)
            # don't allocate because we can store the result
            f!(r, x, y)
            copyto!(xcache, x)
            copyto!(ycache, y)
            copyto!(rcache, r)
            FiniteDiff.finite_difference_jacobian!(drdy, fy!, ycache, cache)
        else
            # allocate because we can't store the result
            fy_new! = (r, y) -> f!(r, x, y)
            f!(r, x, y)
            FiniteDiff.finite_difference_jacobian!(drdy, fy!, y, fdtype)
        end
        return r, drdy
    end

    return dfdy!, fdfdy!
end

function create_residual_output_jacobian_functions(f!, x0, y0, r0, dtype::ForwardAD, sp::SparsePattern)

    # construct function to take derivative of
    xcache = copy(x0)
    fy! = (r, y) -> f!(r, xcache, y)

    # get chunk size for ForwardDiff
    chunk = isnothing(dtype.chunk) ? ForwardDiff.Chunk(y0) : dtype.chunk

    # construct cache
    sparsity = sparse(sp.rows, sp.cols, ones(length(sp.rows)))
    colorvec = SparseDiffTools.matrix_colors(sparsity)
    cache = SparseDiffTools.ForwardColorJacCache(fy!, y0, chunk; dx=r0,
        colorvec=colorvec, sparsity=sparsity)

    # define jacobian function
    dfdy! = function(drdy, x, y)
        copyto!(xcache, x)
        SparseDiffTools.forwarddiff_color_jacobian!(drdy, fy!, y, cache)
        return drdy
    end

    # define combined output and jacobian function
    fdfdy! = function(r, drdy, x, y)
        copyto!(xcache, x)
        f!(r, x, y)
        SparseDiffTools.forwarddiff_color_jacobian!(drdy, fy!, y, cache)
        return r, drdy
    end

    return dfdy!, fdfdy!
end

function create_residual_output_jacobian_functions(f!, x0, y0, r0, dtype::AbstractFD, sp::SparsePattern)

    # construct function to take derivative of
    xcache = copy(x0)
    fy! = (r, y) -> f!(r, xcache, y)

    # Construct FiniteDiff Cache
    ycache = copy(y0)
    rcache = copy(r0) # this needs to be updated before calling FiniteDiff
    rcache1 = copy(r0)
    fdtype = finitediff_type(dtype)
    sparsity = sparse(sp.rows, sp.cols, ones(length(sp.rows)))
    colorvec = SparseDiffTools.matrix_colors(sparsity)
    cache = FiniteDiff.JacobianCache(ycache, rcache, rcache1, fdtype, sparsity,
        colorvec)

    # define jacobian function
    dfdy! = function(drdy, x, y)
        copyto!(xcache, x)
        f!(rcache, x, y)
        FiniteDiff.finite_difference_jacobian!(drdy, fy!, y, cache)
        return drdy
    end

    # define combined output and jacobian function
    fdfdy! = function(r, drdy, x, y)
        f!(r, x, y)
        copyto!(xcache, x)
        copyto!(rcache, r)
        FiniteDiff.finite_difference_jacobian!(drdy, fy!, y, cache)
        return r, drdy
    end

    return dfdy!, fdfdy!
end


text:
import FEM.write_vtk_file

facts("FEM.vtkexport") do

context("FEM.vtkexport") do

mesh = GeoMesh()

node_1 = GeoNode2([0, 0], 1)
node_2 = GeoNode2([1, 1], 2)
node_3 = GeoNode2([1, 2], 3)
node_4 = GeoNode2([0, 1], 4)
addnode!(mesh, node_1)
addnode!(mesh, node_2)
addnode!(mesh, node_3)
addnode!(mesh, node_4)

# Elements
element_1 = GeoTrig([1, 2, 3], 1)
element_2 = GeoTrig([1, 2, 4], 2)
addelem!(mesh, element_1)
addelem!(mesh, element_2)

write_vtk_file(mesh, "test.vtp", true)

end # context

end # facts


text:
module Cameras

using ResourcePools

import Base: isopen, take!

import ResourcePools:
    ref_count,
    release!,
    resource,
    retain!

export Camera,
    open!,
    close!,
    isopen,
    isrunning,
    start!,
    stop!,
    take!,
    trigger!
include("camera.jl")

export AbstractAcquiredImage,
    id,
    timestamp,
    ref_count,
    release!,
    resource,
    retain!
include("acquired_image.jl")

import Base: iterate, IteratorSize
export iterate,
    IteratorSize
include("iteration.jl")

export SimulatedCamera,
    AcquiredImage
include("simulated_camera.jl")

end # module