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

# Extend _add_data_object
function _add_data_object(model::InfiniteModel,
                          object::MeasureData)::MeasureIndex
    return MOIUC.add_item(model.measures, object)
end

# Extend _data_dictionary (type based)
function _data_dictionary(model::InfiniteModel, ::Type{Measure}
    )::MOIUC.CleverDict{MeasureIndex, MeasureData}
    return model.measures
end

# Extend _data_dictionary (reference based)
function _data_dictionary(mref::MeasureRef
    )::MOIUC.CleverDict{MeasureIndex, MeasureData}
    return JuMP.owner_model(mref).measures
end

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

# Extend _core_variable_object
function _core_variable_object(mref::MeasureRef)::Measure
    return _data_object(mref).measure
end

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

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

# Extend _set_core_variable_object
function _set_core_variable_object(mref::MeasureRef, object::Measure)::Nothing
    _data_object(mref).measure = object
    return
end

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

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

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

################################################################################
#                              MEASURE DATA METHODS
################################################################################
# Extend Base.copy for DiscreteMeasureData
function Base.copy(data::DiscreteMeasureData)::DiscreteMeasureData
    return DiscreteMeasureData(copy(data.parameter_refs), copy(data.coefficients),
                               copy(data.supports), data.label,
                               data.weight_function, copy(data.lower_bounds),
                               copy(data.upper_bounds), data.is_expect)
end

# Extend Base.copy for FunctionalDiscreteMeasureData
function Base.copy(data::FunctionalDiscreteMeasureData)::FunctionalDiscreteMeasureData
    return FunctionalDiscreteMeasureData(copy(data.parameter_refs),
                                         data.coeff_function, data.min_num_supports,
                                         data.label, copy(data.generative_supp_info),
                                         data.weight_function,
                                         copy(data.lower_bounds),
                                         copy(data.upper_bounds), data.is_expect)
end

"""
    default_weight(t) = 1

Default weight function for [`DiscreteMeasureData`](@ref) and
[`FunctionalDiscreteMeasureData`](@ref). Returns 1 regardless of the input value.
"""
default_weight(t) = 1

# Ensure that a pref is valid for use in measure data
function _check_params(pref::GeneralVariableRef)::Nothing
    if !(_index_type(pref) <: InfiniteParameterIndex)
        error("$pref is not an infinite parameter.")
    end
    return
end

"""
    DiscreteMeasureData(pref::GeneralVariableRef,
        coefficients::Vector{<:Real},
        supports::Vector{<:Real};
        [label::Type{<:AbstractSupportLabel} = generate_unique_label(),
        weight_function::Function = [`default_weight`](@ref),
        lower_bound::Real = NaN,
        upper_bound::Real = NaN,
        is_expect::Bool = false]
        )::DiscreteMeasureData

Returns a 1-dimensional `DiscreteMeasureData` object that can be utilized
to define measures using [`measure`](@ref). This accepts input for a scalar (single)
infinite parameter. A description of the other arguments is provided in the
documentation for [`DiscreteMeasureData`](@ref). Errors if supports are out
bounds or an unequal number of supports and coefficients are given. Note that by
default a unique `label` is generated via `generate_unique_label` to ensure the supports can
be located in the infinite parameter support storage. Advanced implementations,
may choose a different behavior but should do so with caution.

**Example**
```julia-repl
julia> data = DiscreteMeasureData(pref, [0.5, 0.5], [1, 2])
DiscreteMeasureData{GeneralVariableRef,1,Float64}(pref, [0.5, 0.5], [1.0, 2.0], UniqueMeasure{Symbol("##373")}, default_weight, NaN, NaN, false)
```
"""
function DiscreteMeasureData(pref::GeneralVariableRef,
    coefficients::Vector{<:Real},
    supports::Vector{<:Real};
    label::Type{<:AbstractSupportLabel} = generate_unique_label(),
    weight_function::Function = default_weight,
    lower_bound::Real = NaN,
    upper_bound::Real = NaN,
    is_expect::Bool = false
    )::DiscreteMeasureData{GeneralVariableRef, 1, Float64}
    _check_params(pref)
    if length(coefficients) != length(supports)
        error("The amount of coefficients must match the amount of " *
              "support points.")
    end
    set = infinite_set(pref)
    if !supports_in_set(supports, set)
        error("Support points violate parameter domain.")
    end
    return DiscreteMeasureData(pref, coefficients, supports, label,
                               weight_function, lower_bound, upper_bound,
                               is_expect)
end

# Ensure that an array of prefs is valid for use in multi-dimensional data
function _check_params(prefs::AbstractArray{GeneralVariableRef})::Nothing
    only_dep = all(_index_type(pref) == DependentParameterIndex for pref in prefs)
    only_indep = !only_dep && all(_index_type(pref) == IndependentParameterIndex
                                  for pref in prefs)
    if !only_dep && !only_indep
        error("Cannot specify a mixture of infinite parameter types.")
    elseif only_dep && any(_raw_index(pref) != _raw_index(first(prefs)) for pref in prefs)
        error("Cannot specify multiple dependent parameter groups for one measure.")
    elseif only_dep && length(prefs) != _num_parameters(first(dispatch_variable_ref.(prefs)))
        error("Cannot specify a subset of dependent parameters, consider using " *
              "nested one-dimensional measures instead.")
    end
    return
end

## Define function to intelligently format prefs and supports properly
# Vector
function _convert_param_refs_and_supports(
    prefs::Vector{GeneralVariableRef},
    supports::Vector{Vector{T}}
    )::Tuple{Vector{GeneralVariableRef}, Matrix{T}} where {T <: Real}
    return prefs, reduce(hcat, supports)
end

# DenseAxisArray and Arrays
function _convert_param_refs_and_supports(
    prefs::Union{Array{GeneralVariableRef},
                 JuMPC.DenseAxisArray{GeneralVariableRef}},
#    supports::Vector{<:JuMPC.DenseAxisArray{T}}
    supports::Union{Array{<:Array{T}},Vector{<:JuMPC.DenseAxisArray{T}}}
    )::Tuple{Vector{GeneralVariableRef}, Matrix{T}} where {T <: Real}
    return reduce(vcat, prefs), [supp[i] for i in eachindex(first(supports)), supp in supports]
end

# SparseAxisArray
function _convert_param_refs_and_supports(
    prefs::JuMPC.SparseAxisArray{GeneralVariableRef},
    supports::Vector{<:JuMPC.SparseAxisArray{T}}
    )::Tuple{Vector{GeneralVariableRef}, Matrix{T}} where {T <: Real}
    indices = Collections._get_indices(prefs)
    ordered_prefs = Collections._make_ordered(prefs, indices)
    supps = [Collections._make_ordered(supp, indices) for supp in supports]
    return _convert_param_refs_and_supports(ordered_prefs, supps)
end

## Check if a matrix of supports repsect infinite set(s)
# IndependentParameterRefs
function _check_supports_in_bounds(prefs::Vector{IndependentParameterRef},
                                   supports::Matrix{<:Real})::Nothing
    for i in eachindex(prefs)
        set = infinite_set(prefs[i])
        if !supports_in_set(supports[i, :], set)
            error("Support points violate parameter domain.")
        end
    end
    return
end

# DependentParameterRefs
function _check_supports_in_bounds(prefs::Vector{DependentParameterRef},
                                   supports::Matrix{<:Real})::Nothing
    set = _parameter_set(first(prefs))
    if !supports_in_set(supports, set)
        error("Support points violate parameter domain.")
    end
    return
end

"""
    DiscreteMeasureData(prefs::AbstractArray{GeneralVariableRef},
        coefficients::Vector{<:Real},
        supports::Vector{<:AbstractArray{<:Real}};
        label::Type{<:AbstractSupportLabel} = generate_unique_label(),
        weight_function::Function = [`default_weight`](@ref),
        lower_bounds::AbstractArray{<:Real} = [NaN...],
        upper_bounds::AbstractArray{<:Real} = [NaN...],
        is_expect::Bool = false
        )::DiscreteMeasureData

Returns a `DiscreteMeasureData` object that can be utilized to
define measures using [`measure`](@ref). This accepts input for an array (multi)
parameter. The inner arrays in the supports vector need to match the formatting
of the array used for `parameter_refs`. A description of the other arguments is
provided in the
documentation for [`DiscreteMeasureData`](@ref). Errors if supports are out
bounds, an unequal number of supports and coefficients are given, the array
formats do not match, or if mixed infinite parameter types are given. Note that by
default a unique `label` is generated via `generate_unique_label` to ensure the supports can
be located in the infinite parameter support storage. Advanced implementations,
may choose a different behavior but should do so with caution.

**Example**
```julia-repl
julia> data = DiscreteMeasureData(prefs, [0.5, 0.5], [[1, 1], [2, 2]]);

julia> typeof(data)
DiscreteMeasureData{Array{GeneralVariableRef,1},2,Array{Float64,1}}
```
"""
function DiscreteMeasureData(prefs::AbstractArray{GeneralVariableRef},
    coefficients::Vector{<:Real},
    supports::Vector{<:AbstractArray{<:Real}};
    label::Type{<:AbstractSupportLabel} = generate_unique_label(),
    weight_function::Function = default_weight,
    lower_bounds::AbstractArray{<:Real} = map(e -> NaN, prefs),
    upper_bounds::AbstractArray{<:Real} = map(e -> NaN, prefs),
    is_expect::Bool = false
    )::DiscreteMeasureData{Vector{GeneralVariableRef}, 2, Vector{Float64}}
    _check_params(prefs) # ensures that prefs are valid
    if _keys(prefs) != _keys(first(supports))
        error("Parameter references and supports must use same container type.")
    elseif _keys(prefs) != _keys(lower_bounds)
        error("Parameter references and bounds must use same container type.")
    end
    if length(coefficients) != length(supports)
        error("The amount of coefficients must match the amount of " *
              "support points.")
    end
    vector_prefs, supps = _convert_param_refs_and_supports(prefs, supports)
    _check_supports_in_bounds(dispatch_variable_ref.(vector_prefs), supps)
    lb = _make_ordered_vector(lower_bounds)
    ub = _make_ordered_vector(upper_bounds)
    return DiscreteMeasureData(vector_prefs, coefficients, supps, label,
                               weight_function, lb, ub, is_expect)
end

"""
    FunctionalDiscreteMeasureData(pref::GeneralVariableRef,
        coeff_func::Function,
        min_num_supports::Int,
        label::Type{<:AbstractSupportLabel};
        [weight_function::Function = [`default_weight`](@ref),
        lower_bound::Real = NaN,
        upper_bound::Real = NaN,
        is_expect::Bool = false,
        generative_support_info::AbstractGenerativeInfo = NoGenerativeSupports()]
        )::FunctionalDiscreteMeasureData

Returns a 1-dimensional `FunctionalDiscreteMeasureData` object that can be utilized
to define measures using [`measure`](@ref). This accepts input for a scalar (single)
infinite parameter. A description of the other arguments is provided in the
documentation for [`FunctionalDiscreteMeasureData`](@ref). Errors if `pref` is
not an infinite parameter. Built-in choices for `label` include:
- `All`: Use all of the supports stored in `pref`
- `MCSample`: Use Monte Carlo samples associated with `pref`
- `WeightedSample`: Use weighted Monte Carlo samples associated with `pref`
- `UniformGrid`: Use uniform grid points associated with `pref`.

**Example**
```julia-repl
julia> data = FunctionalDiscreteMeasureData(pref, my_func, 20, UniformGrid)
FunctionalDiscreteMeasureData{GeneralVariableRef,Float64,NoGenerativeSupports}(pref, my_func, 20, UniformGrid, NoGenerativeSupports(), default_weight, NaN, NaN, false)
```
"""
function FunctionalDiscreteMeasureData(
    pref::GeneralVariableRef,
    coeff_func::Function,
    min_num_supports::Int,
    label::Type{<:AbstractSupportLabel};
    weight_function::Function = default_weight,
    lower_bound::Real = NaN,
    upper_bound::Real = NaN,
    is_expect::Bool = false,
    generative_support_info::AbstractGenerativeInfo = NoGenerativeSupports()
    )
    _check_params(pref)
    if _index_type(pref) == DependentParameterIndex && min_num_supports != 0
        error("`min_num_supports` must be 0 for individual dependent parameters.")
    end
    min_num_supports >= 0 || error("Number of supports must be nonnegative.")
    if !isnan(lower_bound) && !isnan(upper_bound)
        if !supports_in_set([lower_bound, upper_bound], infinite_set(pref))
            error("Bounds violate infinite set bounds.")
        end
    end
    return FunctionalDiscreteMeasureData(pref, coeff_func, min_num_supports,
                                         label, generative_support_info,
                                         weight_function, lower_bound,
                                         upper_bound, is_expect)
end

## Check if the integral bounds satisfy the parameter set(s)
# IndependentParameterRefs
function _check_bounds_in_set(prefs::Vector{IndependentParameterRef}, lbs,
                              ubs)::Nothing
    for i in eachindex(prefs)
        if !supports_in_set([lbs[i], ubs[i]], infinite_set(prefs[i]))
            error("Bounds violate the infinite domain.")
        end
    end
    return
end

# DependentParameterRefs
function _check_bounds_in_set(prefs::Vector{DependentParameterRef}, lbs,
                              ubs)::Nothing
    if !supports_in_set(hcat(lbs, ubs), infinite_set(prefs))
        error("Bounds violate the infinite domain.")
    end
    return
end

"""
    FunctionalDiscreteMeasureData(prefs::AbstractArray{GeneralVariableRef},
        coeff_func::Function,
        min_num_supports::Int,
        label::Type{<:AbstractSupportLabel};
        [weight_function::Function = [`default_weight`](@ref),
        lower_bounds::AbstractArray{<:Real} = [NaN...],
        upper_bounds::AbstractArray{<:Real} = [NaN...],
        is_expect::Bool = false]
        )::FunctionalDiscreteMeasureData

Returns a multi-dimensional `FunctionalDiscreteMeasureData` object that can be utilized
to define measures using [`measure`](@ref). This accepts input for an array of
infinite parameters. A description of the other arguments is provided in the
documentation for [`FunctionalDiscreteMeasureData`](@ref). Errors if `prefs` are
not infinite parameters or if the mixed parameter types are provided.
Built-in choices for `label` include:
- `All`: Use all of the supports stored in `prefs`
- `MCSample`: Use Monte Carlo samples associated with `prefs`
- `WeightedSample`: Use weighted Monte Carlo samples associated with `prefs`
- `UniformGrid`: Use uniform grid points associated with `prefs`.

**Example**
```julia-repl
julia> data = FunctionalDiscreteMeasureData(prefs, my_func, 20, MCSample);
```
"""
function FunctionalDiscreteMeasureData(
    prefs::AbstractArray{GeneralVariableRef},
    coeff_func::Function,
    min_num_supports::Int,
    label::Type{<:AbstractSupportLabel};
    weight_function::Function = default_weight,
    lower_bounds::AbstractArray{<:Real} = map(e -> NaN, prefs),
    upper_bounds::AbstractArray{<:Real} = map(e -> NaN, prefs),
    is_expect::Bool = false
    )
    _check_params(prefs)
    if _keys(prefs) != _keys(lower_bounds)
        error("Parameter references and bounds must use same container type.")
    end
    min_num_supports >= 0 || error("Number of supports must be nonnegative.")
    vector_prefs = _make_ordered_vector(prefs)
    vector_lbs = _make_ordered_vector(lower_bounds)
    vector_ubs = _make_ordered_vector(upper_bounds)
    if !isnan(first(vector_lbs)) && !isnan(first(vector_ubs))
        dprefs = map(p -> dispatch_variable_ref(p), vector_prefs)
        _check_bounds_in_set(dprefs, vector_lbs, vector_ubs)
    end
    return FunctionalDiscreteMeasureData(vector_prefs, coeff_func, min_num_supports,
                                         label, weight_function, lower_bounds,
                                         upper_bounds, is_expect)
end

"""
    parameter_refs(data::AbstractMeasureData)::Union{GeneralVariableRef,
                                                     AbstractArray{GeneralVariableRef}}

Return the infinite parameter reference(s) in `data`. This is intended as an
internal function to be used with measure addition. User-defined measure data types
will need to extend this function otherwise an error is thrown.
"""
function parameter_refs(data::AbstractMeasureData)
    error("Function `parameter_refs` not extended for measure data of type $(typeof(data)).")
end

# DiscreteMeasureData
function parameter_refs(data::DiscreteMeasureData{T})::T where {T}
    return data.parameter_refs
end

# FunctionalDiscreteMeasureData
function parameter_refs(data::FunctionalDiscreteMeasureData{T})::T where {T}
    return data.parameter_refs
end

"""
    support_label(data::AbstractMeasureData)::Type{<:AbstractSupportLabel}

Return the label stored in `data` associated with its supports.
This is intended as en internal method for measure creation and ensures any
new supports are added to parameters with such a label.
User-defined measure data types should extend this functionif supports are used,
otherwise an error is thrown.
"""
function support_label(data::AbstractMeasureData)
    error("Function `support_label` not defined for measure data of type " *
          "$(typeof(data)).")
end

# DiscreteMeasureData and FunctionalDiscreteMeasureData
function support_label(data::Union{FunctionalDiscreteMeasureData,
                       DiscreteMeasureData})::Type{<:AbstractSupportLabel}
    return data.label
end

"""
    JuMP.lower_bound(data::AbstractMeasureData)::Union{Float64, Vector{Float64}}

Return the lower bound associated with `data` that defines its domain. This is
intended as an internal method, but may be useful for extensions. User-defined
measure data types should extend this function if desired, otherwise `NaN` is
returned
"""
function JuMP.lower_bound(data::AbstractMeasureData)::Float64
    return NaN
end

# DiscreteMeasureData and FunctionalDiscreteMeasureData
function JuMP.lower_bound(data::Union{FunctionalDiscreteMeasureData,
                                      DiscreteMeasureData})
    return data.lower_bounds
end

"""
    JuMP.upper_bound(data::AbstractMeasureData)::Union{Float64, Vector{Float64}}

Return the lower bound associated with `data` that defines its domain. This is
intended as an internal method, but may be useful for extensions. User-defined
measure data types should extend this function if desired, otherwise `NaN` is
returned.
"""
function JuMP.upper_bound(data::AbstractMeasureData)::Float64
    return NaN
end

# DiscreteMeasureData and FunctionalDiscreteMeasureData
function JuMP.upper_bound(data::Union{FunctionalDiscreteMeasureData,
                                      DiscreteMeasureData})
    return data.upper_bounds
end

## Indicate is the data is from an expectation call
# Fallback
function _is_expect(data::AbstractMeasureData)::Bool
    return false
end

# DiscreteMeasureData and FunctionalDiscreteMeasureData
function _is_expect(data::Union{FunctionalDiscreteMeasureData,
                                DiscreteMeasureData})::Bool
    return data.is_expect
end

"""
    generative_support_info(data::AbstractMeasureData)::AbstractGenerativeInfo

Return the generative support creation info that corresponds to `data`. This is
intended as an internal method and only needs to be extended for user-defined
measure data types that use generative supports.
"""
function generative_support_info(data::AbstractMeasureData)::NoGenerativeSupports
    return NoGenerativeSupports()
end

# FunctionalDiscreteMeasureData
function generative_support_info(data::FunctionalDiscreteMeasureData)
    return data.generative_supp_info
end

"""
    supports(data::AbstractMeasureData)::Array{Float64}

Return the supports associated with `data` and its infinite parameters.
This is intended as en internal method for measure creation and ensures any
new supports are added to parameters. User-defined measure data types should
extend this function if appropriate, otherwise an empty vector is returned.
"""
function supports(data::AbstractMeasureData)::Vector{Float64}
    return Float64[]
end

# DiscreteMeasureData
function supports(data::DiscreteMeasureData{T, N})::Array{Float64, N} where {T, N}
    return data.supports
end

## Define helper functions for generative functional data
# Include generative and uses All
function _get_supports(pref, gen::Val{true}, info, label::Type{All})::Vector{Float64}
    return supports(pref, label = label)
end

# Include generative and does not use All
function _get_supports(pref, gen::Val{true}, info, label)::Vector{Float64}
    return supports(pref, label = Union{support_label(info), label})
end

# Exclude generative and uses All
function _get_supports(pref, gen::Val{false}, info, label::Type{All})::Vector{Float64}
    gen_label = support_label(info)
    return findall(x -> any(v -> v != gen_label, x), _parameter_supports(pref))
end

# Exclude generative and does not use All
function _get_supports(pref, gen::Val{false}, info, label)::Vector{Float64}
    gen_label = support_label(info)
    return supps = findall(x -> any(v -> v <: label && v != gen_label, x),
                           _parameter_supports(pref))
end

# 1D FunctionalDiscreteMeasureData (not NoGenerativeSupports)
function supports(
    data::FunctionalDiscreteMeasureData{GeneralVariableRef};
    include_generative::Bool = true
    )::Vector{Float64}
    info = generative_support_info(data)
    label = support_label(data)
    pref = dispatch_variable_ref(parameter_refs(data))
    supps = _get_supports(pref, Val(include_generative), info, label)
    lb = JuMP.lower_bound(data)
    ub = JuMP.upper_bound(data)
    if isnan(lb) && isnan(ub)
        return supps
    else
        return filter!(i -> lb <= i <= ub, supps)
    end
end

# 1D FunctionalDiscreteMeasureData (NoGenerativeSupports)
function supports(
    data::FunctionalDiscreteMeasureData{GeneralVariableRef, Float64, NoGenerativeSupports};
    include_generative::Bool = true
    )::Vector{Float64}
    supps = supports(parameter_refs(data), label = support_label(data))
    lb = JuMP.lower_bound(data)
    ub = JuMP.upper_bound(data)
    if isnan(lb) && isnan(ub)
        return supps
    else
        return filter!(i -> lb <= i <= ub, supps)
    end
end

# Multi FunctionalDiscreteMeasureData
function supports(
    data::FunctionalDiscreteMeasureData{Vector{GeneralVariableRef}};
    include_generative::Bool = true
    )::Matrix{Float64}
    supps = supports(parameter_refs(data), label = support_label(data))
    lb = JuMP.lower_bound(data)
    ub = JuMP.upper_bound(data)
    if isnan(first(lb)) && isnan(first(ub))
        return supps
    else
        inds = [all(lb .<= @view(supps[:, i]) .<= ub) for i in 1:size(supps, 2)]
        return supps[:, inds]
    end
end

"""
    num_supports(data::AbstractMeasureData)::Int

Return the number supports associated with `data` and its infinite parameters.
This is intended as an internal method for measure creation. User-defined
measure data types should extend this function if appropriate, otherwise
0 is returned.
"""
function num_supports(data::AbstractMeasureData)::Int
    return 0
end

# DiscreteMeasureData/FunctionalDiscreteMeasureData
function num_supports(data::Union{FunctionalDiscreteMeasureData,
                                  DiscreteMeasureData})::Int
    return size(supports(data))[end]
end

"""
    min_num_supports(data::AbstractMeasureData)::Int

Return the minimum number of supports associated with `data`. By fallback, this
will just return `num_supports(data)`. This is primarily intended for internal
queries of `FunctionalDiscreteMeasureData`, but can be extended for other
measure data types if needed.
"""
function min_num_supports(data::AbstractMeasureData)::Int
    return num_supports(data)
end

# FunctionalDiscreteMeasureData
function min_num_supports(data::FunctionalDiscreteMeasureData)::Int
    return data.min_num_supports
end

"""
    coefficient_function(data::AbstractMeasureData)::Function

Return the coefficient function stored in `data` associated with its
expansion abstraction is there is such a function. This is intended as an
internal method for measure creation. User-defined measure
data types should extend this function if appropriate, otherwise an error is
thrown for unsupported types.
"""
function coefficient_function(data::AbstractMeasureData)::Function
    error("Function `coefficient_function` not defined for measure data of type " *
          "$(typeof(data)).")
end

# FunctionalDiscreteMeasureData
function coefficient_function(data::FunctionalDiscreteMeasureData)::Function
    return data.coeff_function
end

"""
    coefficients(data::AbstractMeasureData)::Vector{<:Real}

Return the coefficients associated with `data` associated with its expansion abstraction.
This is intended as en internal method for measure creation. User-defined measure
data types should extend this function if appropriate, otherwise an empty vector
is returned.
"""
function coefficients(data::AbstractMeasureData)::Vector{Float64}
    return Float64[]
end

# DiscreteMeasureData
function coefficients(data::DiscreteMeasureData)::Vector{Float64}
    return data.coefficients
end

# FunctionalDiscreteMeasureData
function coefficients(data::FunctionalDiscreteMeasureData)::Vector{Float64}
    return data.coeff_function(supports(data, include_generative = false))
end

"""
    weight_function(data::AbstractMeasureData)::Function

Return the weight function stored in `data` associated with its expansion abstraction.
This is intended as en internal method for measure creation. User-defined measure
data types should extend this function if appropriate, otherwise an error is thrown.
"""
function weight_function(data::AbstractMeasureData)::Function
    error("Function `weight_function` not defined for measure data of type " *
          "$(typeof(data)).")
end

# DiscreteMeasureData and FunctionalDiscreteMeasureData
function weight_function(data::Union{DiscreteMeasureData,
                                     FunctionalDiscreteMeasureData})::Function
    return data.weight_function
end

################################################################################
#                            MEASURE CONSTRUCTION METHODS
################################################################################
"""
    measure_data_in_hold_bounds(data::AbstractMeasureData,
                                bounds::ParameterBounds)::Bool

Return a `Bool` whether the domain of `data` is valid in accordance with
`bounds`. This is intended as an internal method and is used to check hold
variables used in measures. User-defined measure data types will need to
extend this function to enable this error checking, otherwise it is skipped and
a warning is given.
"""
function measure_data_in_hold_bounds(data::AbstractMeasureData,
                                     bounds::ParameterBounds)::Bool
    @warn "Unable to check if hold variables bounds are valid in measure " *
           "with measure data type `$(typeof(data))`. This can be resolved by " *
           "extending `measure_data_in_hold_bounds`."
    return true
end

# Scalar DiscreteMeasureData and FunctionalDiscreteMeasureData
function measure_data_in_hold_bounds(
    data::Union{DiscreteMeasureData{P, 1}, FunctionalDiscreteMeasureData{P}},
    bounds::ParameterBounds{GeneralVariableRef}
    )::Bool where {P <: GeneralVariableRef}
    pref = parameter_refs(data)
    supps = supports(data)
    if haskey(bounds, pref) && length(supps) != 0
        return supports_in_set(supps, bounds[pref])
    end
    return true
end

# Multi-dimensional DiscreteMeasureData and FunctionalDiscreteMeasureData
function measure_data_in_hold_bounds(
    data::Union{DiscreteMeasureData{P, 2}, FunctionalDiscreteMeasureData{P}},
    bounds::ParameterBounds{GeneralVariableRef}
    )::Bool where {P <: Vector{GeneralVariableRef}}
    prefs = parameter_refs(data)
    supps = supports(data)
    if length(supps) != 0
        for i in eachindex(prefs)
            if haskey(bounds, prefs[i])
                if !supports_in_set(supps[i, :], bounds[prefs[i]])
                    return false
                end
            end
        end
    end
    return true
end

# Check that variables don't violate the parameter bounds
function _check_var_bounds(vref::GeneralVariableRef, data::AbstractMeasureData)
    if _index_type(vref) == HoldVariableIndex
        bounds = parameter_bounds(vref)
        if !measure_data_in_hold_bounds(data, bounds)
            error("Measure bounds violate hold variable bounds.")
        end
    elseif _index_type(vref) == MeasureIndex
        vrefs = _all_function_variables(measure_function(vref))
        for vref in vrefs
            _check_var_bounds(vref, data)
        end
    end
    return
end

"""
    build_measure(expr::JuMP.AbstractJuMPScalar,
                  data::AbstractMeasureData)::Measure

Build and return a [`Measure`](@ref) given the expression to be measured `expr`
using measure data `data`. This principally serves as an internal method for
measure definition. Errors if the supports associated with `data` violate
an hold variable parameter bounds of hold variables that are included in the
measure.
"""
function build_measure(expr::T, data::D;
    )::Measure{T, D} where {T <: JuMP.AbstractJuMPScalar, D <: AbstractMeasureData}
    vrefs = _all_function_variables(expr)
    model = _model_from_expr(vrefs)
    if model.has_hold_bounds
        for vref in vrefs
            _check_var_bounds(vref, data)
        end
    end
    expr_obj_nums = _object_numbers(expr)
    expr_param_nums = _parameter_numbers(expr)
    prefs = parameter_refs(data)
    data_obj_nums = _object_numbers(prefs)
    data_param_nums = [_parameter_number(pref) for pref in prefs]
    # NOTE setdiff! cannot be used here since it modifies object_nums of expr if expr is a single infinite variable
    obj_nums = sort(setdiff(expr_obj_nums, data_obj_nums))
    param_nums = sort(setdiff(expr_param_nums, data_param_nums))
    # check if analytic method should be applied
    lb_nan = isnan(first(JuMP.lower_bound(data)))
    ub_nan = isnan(first(JuMP.upper_bound(data)))
    # NOTE intersect! cannot be used here since it modifies parameter_nums of expr if expr is a single infinite variable
    constant_func = isempty(intersect(expr_param_nums, data_param_nums)) &&
                    ((!lb_nan && !ub_nan) || _is_expect(data))
    return Measure(expr, data, obj_nums, param_nums, constant_func)
end

################################################################################
#                               DEFINITION METHODS
################################################################################
## Define methods to check a measure's generative info and add it if needed
# DependentParameterRef and NoGenerativeSupports
function _check_and_set_generative_info(pref::DependentParameterRef,
    info::NoGenerativeSupports)::Nothing
    return
end

# DependentParameterRef and not NoGenerativeSupports
function _check_and_set_generative_info(pref::DependentParameterRef, info)::Nothing
    error("Cannot use a measure that requires generative supports with respect to " *
          "a dependent infinite parameter.")
end

# IndependentParameterRef and NoGenerativeSupports
function _check_and_set_generative_info(pref::IndependentParameterRef,
    info::NoGenerativeSupports)::Nothing
    return
end

# IndependentParameterRef w/ NoGenerativeSupports
function _check_and_set_generative_info(pref, info1::NoGenerativeSupports, info2)::Nothing
    _set_generative_support_info(pref, info2)
    return
end

# IndependentParameterRef w/ info already
function _check_and_set_generative_info(pref, info1, info2)::Nothing
    info1 != info2 && error("The generative support structure associated with " *
                            "the desired measure conflicts with the current scheme. " *
                            "This is likely due to the derivative method and/or other measures.")
    return
end

# IndependentParameterRef and not NoGenerativeSupports
function _check_and_set_generative_info(pref::IndependentParameterRef,
    info::AbstractGenerativeInfo)::Nothing
    current_info = generative_support_info(pref)
    _check_and_set_generative_info(pref, current_info, info)
    return
end

## Helper methods for adding supports
function _add_supports_to_multiple_parameters(
    prefs::Vector{DependentParameterRef},
    supps::Array{Float64, 2},
    label::Type{<:AbstractSupportLabel}
    )::Nothing
    add_supports(prefs, supps, label = label, check = false)
    return
end

function _add_supports_to_multiple_parameters(
    prefs::Vector{IndependentParameterRef},
    supps::Array{Float64, 2},
    label::Type{<:AbstractSupportLabel}
    )::Nothing
    for i in eachindex(prefs)
        add_supports(prefs[i], supps[i, :], label = label, check = false)
    end
    return
end

"""
    add_supports_to_parameters(data::AbstractMeasureData)::Nothing

Add supports as appropriate with `data` to the underlying infinite parameters.
This is an internal method with by [`add_measure`](@ref) and should be defined
for user-defined measure data types.
"""
function add_supports_to_parameters(data::AbstractMeasureData)
    error("`add_supports_to_parameters` not defined for measures with " *
          "measure data type $(typeof(data)).")
end

## Internal functions for adding measure data supports to the parameter supports
# scalar DiscreteMeasureData
function add_supports_to_parameters(
    data::DiscreteMeasureData{GeneralVariableRef, 1}
    )::Nothing
    pref = parameter_refs(data)
    supps = supports(data)
    label = support_label(data)
    add_supports(pref, supps, label = label, check = false)
    return
end

# multi-dimensional DiscreteMeasureData
function add_supports_to_parameters(
    data::DiscreteMeasureData{Vector{GeneralVariableRef}, 2}
    )::Nothing
    prefs = map(p -> dispatch_variable_ref(p), parameter_refs(data))
    supps = supports(data)
    label = support_label(data)
    _add_supports_to_multiple_parameters(prefs, supps, label)
    return
end

# scalar FunctionalDiscreteMeasureData
function add_supports_to_parameters(
    data::FunctionalDiscreteMeasureData{GeneralVariableRef}
    )::Nothing
    # determine if we need to add more supports
    num_supps = min_num_supports(data)
    curr_num_supps = num_supports(data)
    if curr_num_supps < num_supps
        # prepare the parameter reference
        pref = dispatch_variable_ref(parameter_refs(data))
        if pref isa DependentParameterRef # This is just a last line of defense
            error("min_num_supports must be 0 for individual dependent parameters.")
        end
        # prepare the generation set
        lb = JuMP.lower_bound(data)
        ub = JuMP.upper_bound(data)
        if isnan(lb) || isnan(ub)
            set = infinite_set(pref)
        else
            set = IntervalSet(lb, ub) # assumes lb and ub are in the set
        end
        # generate the needed supports
        label = support_label(data)
        generate_and_add_supports!(pref, set, label,
                                   num_supports = num_supps - curr_num_supps)
    end
    return
end

## generate more supports as needed in accordance with a label
# DependentParameterRefs
function _generate_multiple_functional_supports(
    prefs::Vector{DependentParameterRef},
    num_supps::Int, label::Type{<:AbstractSupportLabel},
    lower_bounds::Vector{<:Number},
    upper_bounds::Vector{<:Number}
    )::Nothing
    # prepare the generation set
    if isnan(first(lower_bounds)) || isnan(first(upper_bounds))
        set = infinite_set(prefs)
    else
        # assumes we have valid bounds
        set = CollectionSet([IntervalSet(lower_bounds[i], upper_bounds[i])
                             for i in eachindex(lower_bounds)])
    end
    # generate the supports
    generate_and_add_supports!(prefs, set, label, num_supports = num_supps)
    return
end

# IndependentParameterRefs
function _generate_multiple_functional_supports(
    prefs::Vector{IndependentParameterRef},
    num_supps::Int, label::Type{<:AbstractSupportLabel},
    lower_bounds::Vector{<:Number},
    upper_bounds::Vector{<:Number}
    )::Nothing
    # we are gauranteed that each have the same number of supports
    for i in eachindex(prefs)
        # prepare the generation set
        if isnan(lower_bounds[i]) || isnan(upper_bounds[i])
            set = infinite_set(prefs[i])
        else
            # assumes lb and ub are in the set
            set = IntervalSet(lower_bounds[i], upper_bounds[i])
        end
        # generate the supports
        generate_and_add_supports!(prefs[i], set, label, num_supports = num_supps)
    end
    return
end

# multi-dimensional FunationalDiscreteMeasureData
function add_supports_to_parameters(
    data::FunctionalDiscreteMeasureData{Vector{GeneralVariableRef}}
    )::Nothing
    min_num_supps = min_num_supports(data)
    curr_num_supps = num_supports(data) # this will error check the support dims
    needed_supps = min_num_supps - curr_num_supps
    if needed_supps > 0
        prefs = map(p -> dispatch_variable_ref(p), parameter_refs(data))
        label = support_label(data)
        lbs = JuMP.lower_bound(data)
        ubs = JuMP.upper_bound(data)
        _generate_multiple_functional_supports(prefs, needed_supps, label, lbs,
                                               ubs)
    end
    return
end

## Helper methods for updating the generative measure dependencies
# NoGenerativeSupports
function _update_generative_measures(prefs, info::NoGenerativeSupports,
                                     mindex)::Nothing
    return
end

# not NoGenerativeSupports and must be an IndependentParameterRef
function _update_generative_measures(pref, info, mindex)::Nothing
    push!(_generative_measures(pref), mindex)
    return
end

"""
    add_measure(model::InfiniteModel, meas::Measure,
                name::String = "measure")::GeneralVariableRef

Add a measure to `model` and return the corresponding measure reference. This
operates in a manner similar to `JuMP.add_variable`. Note this intended
as an internal method.
"""
function add_measure(model::InfiniteModel, meas::Measure,
                     name::String = "measure")::GeneralVariableRef
    # get the expression variables and check validity
    vrefs = _all_function_variables(meas.func)
    for vref in vrefs
        JuMP.check_belongs_to_model(vref, model)
    end
    # get the measure data info and check validity
    data = meas.data
    prefs = parameter_refs(data)
    for pref in prefs
        JuMP.check_belongs_to_model(pref, model)
    end
    # check generative info is compatible and update model if needed
    info = generative_support_info(data)
    _check_and_set_generative_info(dispatch_variable_ref(first(prefs)), info)
    # add supports to the model as needed
    if !meas.constant_func
        add_supports_to_parameters(data)
    end
    # add the measure to the model
    object = MeasureData(meas, name)
    mindex = _add_data_object(model, object)
    mref = _make_variable_ref(model, mindex)
    # update mappings
    for vref in union!(vrefs, prefs)
        push!(_measure_dependencies(vref), mindex)
    end
    _update_generative_measures(prefs, info, mindex)
    return mref
end

"""
    measure_function(mref::MeasureRef)::JuMP.AbstractJuMPScalar

Return the function associated with `mref`.

**Example**
```julia-repl
julia> measure_function(meas)
y(x, t) + 2
```
"""
function measure_function(mref::MeasureRef)::JuMP.AbstractJuMPScalar
    return _core_variable_object(mref).func
end

"""
    measure_data(mref::MeasureRef)::AbstractMeasureData

Return the measure data associated with `mref`.

**Example**
```julia-repl
julia> data = measure_data(meas);

julia> typeof(data)
FunctionalDiscreteMeasureData{Vector{GeneralVariableRef},Vector{Float64}}
```
"""
function measure_data(mref::MeasureRef)::AbstractMeasureData
    return _core_variable_object(mref).data
end

"""
    is_analytic(mref::MeasureRef)::Bool

Return if `mref` is evaluated analytically.

**Example**
```julia-repl
julia> is_analytic(meas)
false
```
"""
function is_analytic(mref::MeasureRef)::Bool
    return _core_variable_object(mref).constant_func
end

## Return an element of a parameter reference tuple given the model, index, and parameter numbers
# IndependentParameterIndex
function _make_param_tuple_element(model::InfiniteModel,
    idx::IndependentParameterIndex,
    param_nums::Vector{Int}
    )::GeneralVariableRef
    return _make_parameter_ref(model, idx)
end

# DependentParametersIndex
function _make_param_tuple_element(model::InfiniteModel,
    idx::DependentParametersIndex,
    param_nums::Vector{Int}
    )::Union{GeneralVariableRef, Vector{GeneralVariableRef}}
    dpref = DependentParameterRef(model, DependentParameterIndex(idx, 1))
    el_param_nums = _data_object(dpref).parameter_nums
    prefs = [GeneralVariableRef(model, idx.value, DependentParameterIndex, i)
             for i in eachindex(el_param_nums) if el_param_nums[i] in param_nums]
    return length(prefs) > 1 ? prefs : first(prefs)
end

"""
    parameter_refs(mref::MeasureRef)::Tuple

Return the tuple of infinite parameters that the measured expression associated
`mref` depends on once the measure has been evaluated. Note that this will
correspond to the parameter dependencies of the measure function excluding those
included in the measure data.

**Example**
```julia-repl
julia> parameter_refs(meas)
(t,)
```
"""
function parameter_refs(mref::MeasureRef)::Tuple
    model = JuMP.owner_model(mref)
    obj_indices = _param_object_indices(model)[_object_numbers(mref)]
    param_nums = _parameter_numbers(mref)
    return Tuple(_make_param_tuple_element(model, idx, param_nums)
                 for idx in obj_indices)
end

# Extend raw_parameter_refs (this is helpful for defining derivatives)
function raw_parameter_refs(mref::MeasureRef)::VectorTuple{GeneralVariableRef}
    return VectorTuple(parameter_refs(mref))
end

# Extend parameter_list (this is helpful for defining derivatives)
function parameter_list(mref::MeasureRef)::Vector{GeneralVariableRef}
    return raw_parameter_refs(mref).values
end

"""
    measure(expr::JuMP.AbstractJuMPScalar,
            data::AbstractMeasureData;
            [name::String = "measure"])::GeneralVariableRef

Return a measure reference that evaluates `expr` using according to `data`.
The measure data `data` determines how the measure is to be evaluated.
Typically, the `DiscreteMeasureData` and the `FunctionalDiscreteMeasureData`
constructors can be used to for `data`. The variable expression `expr` can contain
`InfiniteOpt` variables, infinite parameters, other measure references (meaning
measures can be nested), and constants. Typically, this is called inside of
[`JuMP.@expression`](@ref), [`JuMP.@objective`](@ref), and
[`JuMP.@constraint`](@ref) in a manner similar to `sum`. Note measures are not
explicitly evaluated until [`build_optimizer_model!`](@ref) is called or unless
they are expanded via [`expand`](@ref) or [`expand_all_measures!`](@ref).

**Example**
```julia-repl
julia> tdata = DiscreteMeasureData(t, [0.5, 0.5], [1, 2]);

julia> xdata = DiscreteMeasureData(xs, [0.5, 0.5], [[-1, -1], [1, 1]]);

julia> constr_RHS = @expression(model, measure(g - s + 2, tdata) + s^2)
measure{t}[g(t) - s + 2] + s²

julia> @objective(model, Min, measure(g - 1  + measure(T, xdata), tdata))
measure{xs}[g(t) - 1 + measure{xs}[T(t, x)]]
```
"""
function measure(expr::JuMP.AbstractJuMPScalar,
                 data::AbstractMeasureData;
                 name::String = "measure")::GeneralVariableRef
    model = _model_from_expr(expr)
    if model === nothing
        error("Expression contains no variables or parameters.")
    end
    meas = build_measure(expr, data)
    return add_measure(model, meas, name)
end

"""
    @measure(expr::JuMP.AbstractJuMPScalar,
             data::AbstractMeasureData;
             [name::String = "measure"])::GeneralVariableRef

An efficient wrapper for [`measure`](@ref), please see its doc string for more
information.
"""
macro measure(expr, data, args...)
    _error(str...) = _macro_error(:measure, (expr, args...), str...)
    extra, kw_args, requestedcontainer = _extract_kw_args(args)
    if length(extra) != 0
        _error("Incorrect number of arguments. Must be of form " *
               "@measure(expr, data, name = ...).")
    end
    if !isempty(filter(kw -> kw.args[1] != :name, kw_args))
        _error("Invalid keyword arguments. Must be of form " *
               "@measure(expr, data, name = ...).")
    end
    expression = :( JuMP.@expression(InfiniteOpt._Model, $expr) )
    mref = :( measure($expression, $data; ($(kw_args...))) )
    return esc(mref)
end

################################################################################
#                               NAMING METHODS
################################################################################
"""
    JuMP.name(mref::MeasureRef)::String

Extend `JuMP.name` to return the name associated with a measure
reference.
"""
function JuMP.name(mref::MeasureRef)::String
    object = get(_data_dictionary(mref), JuMP.index(mref), nothing)
    return object === nothing ? "" : object.name
end

"""
    JuMP.set_name(mref::MeasureRef, name::String)::Nothing

Extend `JuMP.set_name` to specify the name of a measure reference.
"""
function JuMP.set_name(mref::MeasureRef, name::String)::Nothing
    _data_object(mref).name = name
    return
end

################################################################################
#                              DEPENDENCY METHODS
################################################################################
"""
    used_by_measure(mref::MeasureRef)::Bool

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

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

"""
    used_by_constraint(mref::MeasureRef)::Bool

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

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

"""
    used_by_objective(mref::MeasureRef)::Bool

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

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

"""
    used_by_derivative(mref::MeasureRef)::Bool

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

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

"""
    is_used(mref::MeasureRef)::Bool

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

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

################################################################################
#                                GENERAL QUERIES
################################################################################
"""
    num_measures(model::InfiniteModel)::Int

Return the number of measures defined in `model`.

**Example**
```julia-repl
julia> num_measures(model)
2
```
"""
function num_measures(model::InfiniteModel)::Int
    return length(_data_dictionary(model, Measure))
end

"""
    all_measures(model::InfiniteModel)::Vector{GeneralVariableRef}

Return the list of all measures added to `model`.

**Examples**
```julia-repl
julia> all_measures(model)
2-element Array{GeneralVariableRef,1}:
 ∫{t ∈ [0, 6]}[w(t, x)]
 𝔼{x}[w(t, x)]
```
"""
function all_measures(model::InfiniteModel)::Vector{GeneralVariableRef}
    vrefs_list = Vector{GeneralVariableRef}(undef, num_measures(model))
    for (i, (index, _)) in enumerate(_data_dictionary(model, Measure))
        vrefs_list[i] = _make_variable_ref(model, index)
    end
    return vrefs_list
end

################################################################################
#                                   DELETION
################################################################################
"""
    JuMP.delete(model::InfiniteModel, mref::MeasureRef)::Nothing

Extend [`JuMP.delete`](@ref) to delete measures. Errors if measure is invalid,
meaning it does not belong to the model or it has already been deleted.

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

julia> delete(model, meas)

julia> print(model)
Min z
Subject to
 z ≥ 0.0
 0 = 0
 g(t) + z ≥ 42.0, ∀ t ∈ [0, 6]
 g(0.5) = 0
```
"""
function JuMP.delete(model::InfiniteModel, mref::MeasureRef)::Nothing
    @assert JuMP.is_valid(model, mref) "Invalid measure reference."
    # Reset the transcription status
    if is_used(mref)
        set_optimizer_model_ready(model, false)
    end
    gvref = _make_variable_ref(JuMP.owner_model(mref), JuMP.index(mref))
    # Remove from dependent measures if there are any
    for mindex in _measure_dependencies(mref)
        meas_ref = dispatch_variable_ref(model, mindex)
        func = measure_function(meas_ref)
        data = measure_data(meas_ref)
        if func isa GeneralVariableRef
            new_func = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
            new_meas = Measure(new_func, data, Int[], Int[], false)
        else
            _remove_variable(func, gvref)
            new_meas = build_measure(func, data)
        end
        _set_core_variable_object(meas_ref, new_meas)
    end
    # Remove from dependent constraints if there are any
    for cindex in _constraint_dependencies(mref)
        cref = _temp_constraint_ref(model, cindex)
        func = JuMP.jump_function(JuMP.constraint_object(cref))
        if func isa GeneralVariableRef
            set = JuMP.moi_set(JuMP.constraint_object(cref))
            new_func = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
            new_constr = JuMP.ScalarConstraint(new_func, set)
            _set_core_constraint_object(cref, new_constr)
            empty!(_object_numbers(cref))
            empty!(_measure_dependencies(cref))
        else
            _remove_variable(func, gvref)
            _data_object(cref).object_nums = sort(_object_numbers(func))
            filter!(e -> e != JuMP.index(mref), _measure_dependencies(cref))
        end
    end
    # Remove from objective if used there
    if used_by_objective(mref)
        if JuMP.objective_function(model) isa GeneralVariableRef
            new_func = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
            JuMP.set_objective_function(model, new_func)
            JuMP.set_objective_sense(model, MOI.FEASIBILITY_SENSE)
        else
            _remove_variable(JuMP.objective_function(model), gvref)
        end
    end
    # delete associated derivative variables and mapping
    for index in _derivative_dependencies(mref)
        JuMP.delete(model, dispatch_variable_ref(model, index))
    end
    # Update that the variables used by it are no longer used by it
    prefs = parameter_refs(measure_data(mref))
    vrefs = _all_function_variables(measure_function(mref))
    union!(vrefs, prefs)
    for vref in vrefs
        filter!(e -> e != JuMP.index(mref), _measure_dependencies(vref))
    end
    # Remove any unique supports associated with this measure
    if num_supports(measure_data(mref)) > 0
        label = support_label(measure_data(mref))
        if label <: UniqueMeasure
            delete_supports(prefs, label = label)
        end
    end
    # update the generative info if needed
    info = generative_support_info(measure_data(mref))
    if prefs isa GeneralVariableRef && !(info isa NoGenerativeSupports)
        filter!(e -> e != JuMP.index(mref), _generative_measures(prefs))
        if isempty(_generative_measures(prefs)) && derivative_method(prefs) isa NonGenerativeDerivativeMethod
            _set_generative_support_info(dispatch_variable_ref(prefs),
                                         NoGenerativeSupports())
        end
    end
    # delete remaining measure information
    _delete_data_object(mref)
    return
end


text:
using TensND
using Documenter
using SymPy

DocMeta.setdocmeta!(
    TensND,
    :DocTestSetup,
    :(using TensND, LinearAlgebra, SymPy, Tensors, OMEinsum, Rotations);
    recursive = true,
)

use_KaTeX = false

# Use MathJax syntax to define macros:
# - macro names should not be preceded by "\\"
macros = Dict(
    # "cellindices" => "\\mathcal{P}",
    # "conj" => "\\operatorname{conj}",
    # "D" => "\\mathrm d",
    # "dbldot" => "\\mathbin{\\mathord{:}}",
    # "dft" => "\\operatorname{DFT}",
    # "element" => "\\mathrm{e}",
    # "fftfreq" => "\\operatorname{Z}",
    # "I" => "\\mathrm i",
    # "integers" => "\\mathbb Z",
    # "naturals" => "\\mathbb N",
    # "PI" => "\\mathrm{\\pi}",
    # "reals" => "\\mathbb R",
    # "sinc" => "\\operatorname{sinc}",
    # "symotimes" => "\\stackrel{\\mathrm{s}}{\\otimes}",
    # "strain" => "\\boldsymbol{\\varepsilon}",
    # "strains" => "\\mathcal E",
    # "stress" => "\\boldsymbol{\\sigma}",
    # "stresses" => "\\mathcal S",
    # #"tens" => "\\@ifstar{\\boldsymbol}{\\mathbf}",
    # "tens" => "\\boldsymbol",
    # "tensors" => "\\mathcal T",
    # "tr" => "\\operatorname{tr}",
    # "tuple" => "\\mathsf",
    # "vec" => "\\boldsymbol",
)

if use_KaTeX
    macros2 = Dict()
    for (key, val) ∈ pairs(macros)
        macros2["\\"*key] = val
    end
    macros = macros2
end

mathengine =
    use_KaTeX ?
    KaTeX(
        Dict{Symbol,Any}(
            :delimiters => Dict{Any,Any}[
                Dict(:left => "\$", :right => "\$", display => false),
                Dict(:left => "\$\$", :right => "\$\$", display => true),
                Dict(:left => "\\[", :right => "\\]", display => true),
            ],
            :macros => macros,
        ),
    ) :
    MathJax3(
        Dict{Symbol,Any}(
            :tex => Dict{String,Any}(
                "macros" => macros,
                #"packages" => ["base", "ams", "autoload", "bm"],
                "inlineMath" => [["\$", "\$"], ["\\(", "\\)"]],
                "tags" => "ams",
            ),
            :options => Dict(
                "ignoreHtmlClass" => "tex2jax_ignore",
                "processHtmlClass" => "tex2jax_process",
            ),
        ),
        true,
    )

format = Documenter.HTML(;
    prettyurls = get(ENV, "CI", "false") == "true",
    canonical = "https://jfbarthelemy.github.io/TensND.jl",
    assets = String[],
    mathengine = mathengine,
)

# format = Documenter.LaTeX(platform = "none")

# open(joinpath(@__DIR__, "src", "assets", "custom.sty"), "w") do io
#     write(io, "% WARNING: Automatically generated -- Do NOT edit by hand!\n")
#     write(io, "% Script: $(@__FILE__)\n")
#     write(io, "% Date: $(now())\n\n")
#     for (name, definition) in
#         macros write(io, "\\providecommand{\\$name}{$definition}\n")
#     end
# end

makedocs(;
    modules = [TensND],
    authors = "Jean-François Barthélémy <jfbarthelemy@users.noreply.github.com> and contributors",
    repo = "https://github.com/jfbarthelemy/TensND.jl/blob/{commit}{path}#{line}",
    sitename = "TensND.jl",
    format = format,
    pages = Any[
        "Home" => "index.md",
        "Manual" => [
            "man/getting_started.md",
            "man/bases.md",
            "man/tensors.md",
            "man/coorsystems.md",
            ],
        "Tutorials" => [
            "tuto/nlayersphere.md",
            ],
        "API" => "api.md",
    ],
)

deploydocs(; repo = "github.com/jfbarthelemy/TensND.jl", devbranch = "main")

ext: jl


text:
#!/usr/bin/env julia

const N = 440_000_000

get_cache() = (0, ntuple(i -> i^i, 9)...)

function is_munchausen(number, cache)
    n = number
    total = 0

    while n > 0
        (n, digit) = divrem(n, 10)
        @inbounds total += cache[digit+1]
        if total > number
            return false
        end
    end

    return total == number
end

function main()
    cache = get_cache()

    for n in 0:N
        if is_munchausen(n, cache)
            println(n)
        end
    end
end

main()

max_stars_repo_name: uninhm/SpeedTests

max_issues_repo_name: uninhm/SpeedTests

max_forks_repo_name: uninhm/SpeedTests


text:
# first pass: read whole file, save in binary format, count IDs for
# collation, store AM spell strings in the order of appearance

using AMDB: data_path, all_data_files, data_file
using ByteParsers:
    Line, DateYYYYMMDD, PosInteger, ViewBytes, parsedtype, ByteVector, parsenext,
    isparsed, Skip
using DocStringExtensions: SIGNATURES
using CodecZlib: GzipDecompressorStream
using WallTimeProgress: WallTimeTracker, increment!

"""
A parser for just the columns we need for labor market transitions.
"""
const lineparser = Line(PosInteger(Int64), # PENR
                        Skip(),            # start, skip
                        Skip(),            # end, skip
                        Skip(),            # firm, skip
                        DateYYYYMMDD(),    # STICHTAG
                        ViewBytes())       # AM

"""
   $SIGNATURES

Read stream from `io` line by line, extract

1. person id,

2. sample date,

3. labor market status

then reclassify the last one with `aggregator`, then return the assembled
histories as a dictionary of person id to 12-element vectors.

The latter contain the histories of individuals. The value `missing` is used
when no observation is found in the sample.
"""
function read_aggregated_AMs(io::IO, aggregator::Associative{ByteVector, T}, missing::T,
                             lineparser = lineparser,
                             tracker = WallTimeTracker(10_000_000;
                                                       item_name = "line")) where T
    individuals = Dict{Int32, Vector{T}}()
    while !eof(io)
        line_content = readuntil(io, 0x0a)
        record = parsenext(lineparser, line_content, 1, UInt8(';'))
        if isparsed(record)
            (penr, stich, am) = unsafe_get(record)
            mon = Dates.month(stich) # just the month
            state = aggregator[am]   # look up from dictionary
            if !haskey(individuals, penr)
                record = individuals[penr] = fill(missing, 12)
            else
                record = individuals[penr]
            end
            record[mon] = state
            increment!(tracker) # display progress
        end
    end
    individuals
end

function read_aggregated_AMs(filename::String, args...)
    io_gz = open(filename, "r")
    io = GzipDecompressorStream(io_gz)
    individuals = read_aggregated_AMs(io, args...)
    close(io)
    close(io_gz)
    individuals
end

"""
    $SIGNATURES

Expand the shell mapping from a compact format of a dictionary of

```julia
classification => [keys...]
```

to a dictionary of

```julia
key => classification
```
"""
function expand_spell_mapping(spell_mapping::Dict{T, Vector{String}}) where T
    aggregator = Dict{ByteVector, T}()
    for (k, vs) in spell_mapping
        for v in vs
            aggregator[convert(ByteVector, v)] = k
        end
    end
    aggregator
end

# runtime code here

spell_mapping = Dict(# unemployed
                     Int8(1) => String["D2", "SR", "VM", "AS", "LS", "SC", "AL",
                                       "AO","66"],
                     # NA
                     Int8(0) => String["KD", "SV", "TO",
                                       # FIXME PLEASE PUT THESE SOMETHERE
                                       "AO", "AL", "SC", "VM", "66", "LS"],
                     # employed
                     Int8(2) => String["FBENG", "FBUN1", "FBKOM", "LE", "FBEB2",
                                       "FBBS1", "FBEB1", "FBSOL", "LFIBA", "FBKUA",
                                       "LFNRM", "FBEPU", "LFP30", "FBSOB", "SO",
                                       "FBES", "LW", "LFUBA", "SBSVA", "FBBS3",
                                       "FBUN3", "LFVOL", "FBUN2", "S1S2", "G1", "FU",
                                       "FBBP", "FBEB", "FBARB", "BE", "AA", "LFLST",
                                       "LFTEL", "LFVRL", "FD", "LFJAS", "FBGEB",
                                       "FBBEB", "FBES1", "FBBS2", "FBGBP"],
                     # other (inactive)
                     Int8(1) => String["SF", "KG", "AU", "W1", "TA", "MS", "BA",
                                       "ED", "LL", "AG", "AF", "EO", "KO", "RE",
                                       "MK", "PZ", "AM", "LF", "MP", "SG", "W2"])

aggregator = expand_spell_mapping(spell_mapping)
test_file = data_file(2000)
histories = read_aggregated_AMs(test_file, aggregator, Int8(-1))


text:
approximator_GNN = GeometricFlux.GraphConv(64 => 64, Flux.leakyrelu)
target_approximator_GNN = GeometricFlux.GraphConv(64 => 64, Flux.leakyrelu)
gnnlayers = 1
approximator_model = SeaPearl.CPNN(
    graphChain = Flux.Chain(
        GeometricFlux.GraphConv(3 => 64, Flux.leakyrelu),
        [approximator_GNN for i = 1:gnnlayers]...
    ),
    nodeChain = Flux.Chain(
        Flux.Dense(64, 32, Flux.leakyrelu),
        Flux.Dense(32, 32, Flux.leakyrelu),
        Flux.Dense(32, 16, Flux.leakyrelu),
    ),
    outputChain = Flux.Dense(16, 4),
) |> cpu
target_approximator_model = SeaPearl.CPNN(
    graphChain = Flux.Chain(
        GeometricFlux.GraphConv(3 => 64, Flux.leakyrelu),
        [target_approximator_GNN for i = 1:gnnlayers]...
    ),
    nodeChain = Flux.Chain(
        Flux.Dense(64, 32, Flux.leakyrelu),
        Flux.Dense(32, 32, Flux.leakyrelu),
        Flux.Dense(32, 16, Flux.leakyrelu),
    ),
    outputChain = Flux.Dense(16, 4),
) |> cpu

agent = RL.Agent(
    policy = RL.QBasedPolicy(
        learner = RL.DQNLearner(
            approximator = RL.NeuralNetworkApproximator(
                model = approximator_model,
                optimizer = ADAM(0.001f0)
            ),
            target_approximator = RL.NeuralNetworkApproximator(
                model = target_approximator_model,
                optimizer = ADAM(0.001f0)
            ),
            loss_func = Flux.Losses.huber_loss,
            stack_size = nothing,
            γ = 0.99f0,
            batch_size = 32,
            update_horizon = 1,
            min_replay_history = 1,
            update_freq = 1,
            target_update_freq = 100,
        ),
        explorer = RL.EpsilonGreedyExplorer(
            ϵ_stable = 0.01,
            kind = :exp,
            ϵ_init = 1.0,
            warmup_steps = 0,
            decay_steps = 500,
            step = 1,
            is_break_tie = false,
            #is_training = true,
        )
    ),
    trajectory = RL.CircularArraySLARTTrajectory(
        capacity = 500,
        state = SeaPearl.DefaultTrajectoryState[] => (),
        action = Int => (),
        legal_actions_mask = Vector{Bool} => (4, ),
    )
)

@testset "learning.jl" begin
    include("environment.jl")
    include("rewards/rewards.jl")
    include("utils.jl")
    include("learnedheuristic.jl")

    LearnedHeristicBasicConstructor = SeaPearl.LearnedHeuristic(agent)
    @test isa(LearnedHeristicBasicConstructor,SeaPearl.LearnedHeuristic{SeaPearl.DefaultStateRepresentation{SeaPearl.DefaultFeaturization, SeaPearl.DefaultTrajectoryState}, SeaPearl.DefaultReward, SeaPearl.FixedOutput})
    @test isnothing(LearnedHeristicBasicConstructor.fitted_problem)
    @test isnothing(LearnedHeristicBasicConstructor.fitted_strategy)
    @test isnothing(LearnedHeristicBasicConstructor.action_space)
    @test isnothing(LearnedHeristicBasicConstructor.current_state)
    @test isnothing(LearnedHeristicBasicConstructor.reward)
    @test isnothing(LearnedHeristicBasicConstructor.search_metrics)

end


text:
# just-plot-it.jl

using Plots

# Choose a friendly plotting package. PyPlot behaves nicely.
# pyplot() # Stopped working after a recent update.
pyplot()

# Make data shaped like the data that had the problem.
t      = collect(0:3600) / 60.
q_BI   = randn(4, 3601)
ω_BI_B = randn(3, 3601)

# Define each plot that we'll need.
p1 = plot(t, q_BI',
            label  = ["q1" "q2" "q3" "q4"],
            ylabel = "Attitude Quaternion")
p2 = plot(t, 180/π * ω_BI_B',
            label  = ["ω1" "ω2" "ω3"],
            ylabel = "Rotation Rate (deg/s)")

# Display those plots as subplots using the default "2" layout.
display(plot(p1, p2,
             layout = 2,
             xlabel = "Time (s)"))

# If we comment out the above and only display one plot, it works the right way.
# display(p1)


text:
using ExecutableSpecifications.Gherkin: issuccessful, parsescenario, Given, When, Then, ByLineParser, ScenarioStep

@testset "Scenario        " begin
    @testset "Scenario has a Given step; the parsed scenario has a Given struct" begin
        text = """
        Scenario: Some description
            Given a precondition
        """

        byline = ByLineParser(text)
        result = parsescenario(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps == ScenarioStep[Given("a precondition")]
    end

    @testset "Scenario has a When step; the parsed scenario has a When struct" begin
        text = """
        Scenario: Some description
            When some action
        """

        byline = ByLineParser(text)
        result = parsescenario(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps == ScenarioStep[When("some action")]
    end

    @testset "Scenario has a Then step; the parsed scenario has a Then struct" begin
        text = """
        Scenario: Some description
            Then a postcondition
        """

        byline = ByLineParser(text)
        result = parsescenario(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps == ScenarioStep[Then("a postcondition")]
    end

    @testset "Scenario has an And following a Given; the And step becomes a Given" begin
        text = """
        Scenario: Some description
            Given a precondition
              And another precondition
        """

        byline = ByLineParser(text)
        result = parsescenario(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps == ScenarioStep[Given("a precondition"),
                                                 Given("another precondition")]
    end

    @testset "Scenario has an And following a When; the And step becomes a When" begin
        text = """
        Scenario: Some description
            When some action
             And another action
        """

        byline = ByLineParser(text)
        result = parsescenario(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps == ScenarioStep[When("some action"),
                                                 When("another action")]
    end

    @testset "Scenario has an And following a Then; the And step becomes a Then" begin
        text = """
        Scenario: Some description
            Then some postcondition
             And another postcondition
        """

        byline = ByLineParser(text)
        result = parsescenario(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps == ScenarioStep[Then("some postcondition"),
                                                 Then("another postcondition")]
    end

    @testset "Scenario is not terminated by newline; EOF is also an OK termination" begin
        text = """
        Scenario: Some description
            Then some postcondition
            And another postcondition"""

        byline = ByLineParser(text)
        result = parsescenario(byline)

        @test issuccessful(result)
    end

    @testset "Malformed scenarios" begin
        @testset "And as a first step; Expected Given, When, or Then before that" begin
            text = """
            Scenario: Some description
                And another postcondition
            """

            byline = ByLineParser(text)
            result = parsescenario(byline)

            @test !issuccessful(result)
            @test result.reason == :leading_and
            @test result.expected == :specific_step
            @test result.actual == :and_step
        end

        @testset "Given after a When; Expected When or Then" begin
            text = """
            Scenario: Some description
                 When some action
                Given some precondition
            """

            byline = ByLineParser(text)
            result = parsescenario(byline)

            @test !issuccessful(result)
            @test result.reason == :bad_step_order
            @test result.expected == :NotGiven
            @test result.actual == :Given
        end

        @testset "Given after Then; Expected Then" begin
            text = """
            Scenario: Some description
                Then some postcondition
                Given some precondition
            """

            byline = ByLineParser(text)
            result = parsescenario(byline)

            @test !issuccessful(result)
            @test result.reason == :bad_step_order
            @test result.expected == :NotGiven
            @test result.actual == :Given
        end

        @testset "When after Then; Expected Then" begin
            text = """
            Scenario: Some description
                Then some postcondition
                When some action
            """

            byline = ByLineParser(text)
            result = parsescenario(byline)

            @test !issuccessful(result)
            @test result.reason == :bad_step_order
            @test result.expected == :NotWhen
            @test result.actual == :When
        end

        @testset "Invalid step definition NotAStep; Expected a valid step definition" begin
            text = """
            Scenario: Some description
                NotAStep some more text
            """

            byline = ByLineParser(text)
            result = parsescenario(byline)

            @test !issuccessful(result)
            @test result.reason == :invalid_step
            @test result.expected == :step_definition
            @test result.actual == :invalid_step_definition
        end

        @testset "A step definition without text; Expected a valid step definition" begin
            text = """
            Scenario: Some description
                Given
            """

            byline = ByLineParser(text)
            result = parsescenario(byline)

            @test !issuccessful(result)
            @test result.reason == :invalid_step
            @test result.expected == :step_definition
            @test result.actual == :invalid_step_definition
        end
    end

    @testset "Block text" begin
        @testset "Block text in a Given; Block text is present in step" begin
            text = """
            Scenario: Some description
                Given some precondition
                \"\"\"
                This is block text.
                There are two lines.
                \"\"\"
            """

            byline = ByLineParser(text)
            result = parsescenario(byline)

            @test issuccessful(result)
            scenario = result.value

            @test scenario.steps[1].block_text == """
            This is block text.
            There are two lines."""
        end

        @testset "Another block text in a Given; Block text is present in step" begin
            text = """
            Scenario: Some description
                Given some precondition
                \"\"\"
                This is another block text.
                There are three lines.
                This is the last line.
                \"\"\"
            """

            byline = ByLineParser(text)
            result = parsescenario(byline)

            @test issuccessful(result)
            scenario = result.value

            @test scenario.steps[1].block_text == """
            This is another block text.
            There are three lines.
            This is the last line."""
        end

        @testset "Block text in a When step; Block text is present in step" begin
            text = """
            Scenario: Some description
                When some action
                \"\"\"
                This is block text.
                \"\"\"
            """

            byline = ByLineParser(text)
            result = parsescenario(byline)

            @test issuccessful(result)
            scenario = result.value

            @test scenario.steps[1] == When("some action"; block_text="""This is block text.""")
        end

        @testset "Block text in a Then step; Block text is present in step" begin
            text = """
            Scenario: Some description
                Then some postcondition
                \"\"\"
                This is block text.
                \"\"\"
            """

            byline = ByLineParser(text)
            result = parsescenario(byline)

            @test issuccessful(result)
            scenario = result.value

            @test scenario.steps[1] == Then("some postcondition"; block_text="""This is block text.""")
        end
    end
end


text:
using IntroLinearAlgebra
using Test

# write your own tests here
@test rowswitch([1 2 3 4; 5 6 7 8],1,2) == [5 6 7 8; 1 2 3 4]

@test rowscale([1 2 3; 4 5 6],2,-3//2) == [1 2 3; -6 -15//2 -9]

@test rowadd([0.3 0.6; π 2π],2,1,-π/0.3) == [0.3 0.6; 0.0 0.0]

@test rref([1 2 3 4 5; 6 7 8 9 10]) == [1 0 -1 -2 -3; 0 1 2 3 4]

@test length(rref(randn(8,8))) == 64

A = [1//1 2 3 4;
        5 6 7 8]

zerooutbelow!(A,1,1)

@test A == [1 2 3 4; 0 -4 -8 -12]


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using DelimitedFiles
using Random
Random.seed!(233)

function halo_exchange_neighbor_two(u,fill_value,m,n,tag,w)
    halo_exchange_neighbor_two_ = load_op_and_grad("./build/libHaloExchangeNeighborTwo","halo_exchange_neighbor_two")
    u,fill_value,m,n,tag,w = convert_to_tensor(Any[u,fill_value,m,n,tag,w], [Float64,Float64,Int64,Int64,Int64,Float64])
    halo_exchange_neighbor_two_(u,fill_value,m,n,tag,w)
end

mpi_init()
u = reshape(Array(1:36), 6, 6)|>Array

if mpi_rank()==0
    u = u[1:3,1:3]
end
if mpi_rank()==1
    u = u[1:3,4:end]
end

if mpi_rank()==2
    u = u[4:end,1:3]
end
if mpi_rank()==3
    u = u[4:end,4:end]
end

fill_value = 10.0
m = 2
n = 2
tag = 1
w = 1.0
# TODO: specify your input parameters
uext = halo_exchange_neighbor_two(u,fill_value,m,n,tag,w)
sess = Session(); init(sess)
uval =  run(sess, uext)

sleep(mpi_rank())
@info "rank = $(mpi_rank())===================================================="
writedlm(stdout, round.(uval, digits=3))
# uncomment it for testing gradients

mpi_finalize()


text:
function velocity_normalization_constant(data::AbstractVector)
    T = length(data)
    c = 0.0
    for δ = 1:T-1
        c += velocity(data, δ)^2
    end
    c /= (T - 1)
    sqrt(c)
end

function velocity(data::AbstractVector, n::Int)
    data[n+1] - data[n]
end

# n（タイムラグ）
function _dynamical_functional(data::AbstractVector, n::Int, ω)
    @argcheck (n >= 1) && (ω >= 0)
    T = length(data)

    e1 = 0.0
    @inbounds for k = 1:(T-n) # kは開始地点
        @views e1 += exp(ω * im * (data[k+n] - data[k]))
    end
    e1 /= (T - n)

    e2 = 0.0
    @inbounds for k = 1:T
        @views e2 += exp(ω * im * (data[k] - data[1]))
    end
    e2 = abs2(e2) / (T * (T - 1))
    e1 - e2 + (1 / (T + 1))
end

function dynamical_functional(
    df::DataFrame,
    id::Symbol = :TrackID,
    x::Symbol = :POSITION_X;
    # frame::Symbol = :FRAME2;
    ω = 1,
)
    N = maximum(df[!, id])
    # T = Int(maximum(df[!, frame])) - 1
    T = 1000-1
    E = Matrix{Union{ComplexF64,Nothing}}(nothing, T, N)
    @inbounds for n = 1:N
        data = spt.extract(df, n, id, x)
        _T = size(data, 1)-1
        @simd for t = 1:_T
            E[t, n] = spt._dynamical_functional(data, t, ω)
        end
    end
    E
end

function ergodicity_estimator(E)
    T, N = size(E)
    F = Matrix{Union{ComplexF64,Nothing}}(nothing, T, N)
    @inbounds for n = 1:N
        @views @simd for t = 1:size(filter(!isnothing, E[:, n]), 1)
            @views F[t, n] = sum(E[1:t, n])/t
        end
    end
    F
end

max_stars_repo_name: SosUts/spt

max_issues_repo_name: SosUts/spt

max_forks_repo_name: SosUts/spt


text:
using Revise
using IsopycnalSurfaces, Test

@testset "IsopycnalSurfaces.jl" begin
    #############################
    # Test dedup! function
    @testset "deduplication" begin
        na = 10
        for ndupes = 1:na-1
            println(ndupes)
            a = sort(randn(na))
            b = randn(na)
            # make duplicates at end
            counter = 0
            while counter < ndupes
                a[end-counter-1] = a[end-counter]
                counter += 1
            end
            dedup!(a,b)
            println(a)
            println(b)
            println(@test issorted(a))
        end
    end
end


text:
# supported op: +, -, *, /, max, min, &, |, mean

function scatter_kernel!(op, dst, src, idx)
    index = threadIdx().x + (blockIdx().x - 1) * blockDim().x

    @inbounds if index <= length(idx)
        CUDA.@atomic dst[idx[index]...] = op(dst[idx[index]...], src[index])
    end
    return nothing
end

function scatter_kernel!(op, dst, src, idx::CUDA.CuDeviceArray{<:CartesianIndex})
    index = threadIdx().x + (blockIdx().x - 1) * blockDim().x

    @inbounds if index <= length(idx)
        li = Base._to_linear_index(dst, Tuple(idx[index])...)
        CUDA.@atomic dst[li] = op(dst[li], src[index])
    end
    return nothing
end

function scatter_kernel!(op, dst, src, idx, max_idx, max_dims_idx, dims_size)
    index = threadIdx().x + (blockIdx().x - 1) * blockDim().x

    @inbounds if index <= max_idx
        j, k = divrem(index-1, max_dims_idx)
        dims_i = CartesianIndices(dims_size)[k+1]
        CUDA.@atomic dst[Tuple(dims_i)..., idx[j+1]...] = op(dst[Tuple(dims_i)..., idx[j+1]...], src[index])
    end
    return nothing
end

function scatter_kernel!(op, dst, src, idx::CUDA.CuDeviceArray{<:CartesianIndex}, max_idx, max_dims_idx, dims_size)
    index = threadIdx().x + (blockIdx().x - 1) * blockDim().x

    @inbounds if index <= max_idx
        j, k = divrem(index-1, max_dims_idx)
        dims_i = CartesianIndices(dims_size)[k+1]
        li = Base._to_linear_index(dst, Tuple(dims_i)..., Tuple(idx[j+1])...)
        CUDA.@atomic dst[li] = op(dst[li], src[index])
    end
    return nothing
end

function NNlib.scatter!(op, dst::AnyCuArray, src::AnyCuArray, idx::AnyCuArray)
    dims = NNlib.scatter_dims(dst, src, idx)
    args = if dims == 0
        max_idx = length(idx)
        op, dst, src, idx
    else
        dims_size = size(dst)[1:dims]
        max_dims_idx = prod(dims_size)
        max_idx = max_dims_idx * length(idx)
        op, dst, src, idx, max_idx, max_dims_idx, dims_size
    end

    kernel = @cuda launch=false scatter_kernel!(args...)
    config = launch_configuration(kernel.fun; max_threads=256)
    threads = min(max_idx, config.threads)
    blocks = cld(max_idx, threads)
    kernel(args...; threads=threads, blocks=blocks)
    return dst
end

function NNlib.scatter!(op::typeof(mean), dst::AnyCuArray, src::AnyCuArray, idx::AnyCuArray)
    Ns = NNlib.scatter!(+, zero(dst), one.(src), idx)
    dst_ = NNlib.scatter!(+, zero(dst), src, idx)
    dst .+= NNlib.safe_div.(dst_, Ns)
    return dst
end

## Gradients

function ∇scatter_src_kernel!(op, Δsrc, src, idx, rev_idx, max_idx, T)
    index = threadIdx().x + (blockIdx().x - 1) * blockDim().x

    @inbounds if index <= max_idx
        cart_j = CartesianIndices(idx)[index]
        # get aggregating indeices, which is to be aggregated together, and itself index
        inds = rev_idx[idx[cart_j]...]
        # multiply all values to be aggregated but not itself
        x = one(T)
        for k in inds
            x *= src[k]
        end
        x /= src[cart_j]
        # apply `op` on `Δsrc[i, k]` and `x`
        Δsrc[cart_j] = op(Δsrc[cart_j], x)
    end
    return nothing
end

function ∇scatter_src_kernel!(op, Δsrc, src, idx::CUDA.CuDeviceArray{<:CartesianIndex}, rev_idx, max_idx, T)
    index = threadIdx().x + (blockIdx().x - 1) * blockDim().x

    @inbounds if index <= max_idx
        cart_j = CartesianIndices(idx)[index]
        # get aggregating indeices, which is to be aggregated together, and itself index
        inds = rev_idx[Tuple(idx[cart_j])...]
        # multiply all values to be aggregated but not itself
        x = one(T)
        for k in inds
            x *= src[k]
        end
        x /= src[cart_j]
        # apply `op` on `Δsrc[i, k]` and `x`
        Δsrc[cart_j] = op(Δsrc[cart_j], x)
    end
    return nothing
end

function ∇scatter_src_kernel!(op, Δsrc, src, idx, rev_idx, pre_cart_idx, max_dims_idx, max_idx, T)
    index = threadIdx().x + (blockIdx().x - 1) * blockDim().x

    @inbounds if index <= max_idx
        i, j = fldmod1(index, max_dims_idx)
        cart_i = CartesianIndices(idx)[i]
        cart_j = pre_cart_idx[j]
        # get aggregating indeices, which is to be aggregated together, and itself index
        inds = rev_idx[idx[cart_i]...]
        # multiply all values to be aggregated but not itself
        x = one(T)
        for k in inds
            jk = Base._to_linear_index(src, Tuple(cart_j)..., Tuple(k)...)
            x *= src[jk]
        end
        x /= src[index]
        # apply `op` on `Δsrc[i, k]` and `x`
        Δsrc[index] = op(Δsrc[index], x)
    end
    return nothing
end

function ∇scatter_src_kernel!(op, Δsrc, src, idx::CUDA.CuDeviceArray{<:CartesianIndex}, rev_idx, pre_cart_idx, max_dims_idx, max_idx, T)
    index = threadIdx().x + (blockIdx().x - 1) * blockDim().x

    @inbounds if index <= max_idx
        i, j = fldmod1(index, max_dims_idx)
        cart_i = CartesianIndices(idx)[i]
        cart_j = pre_cart_idx[j]
        # get aggregating indeices, which is to be aggregated together, and itself index
        inds = rev_idx[Tuple(idx[cart_i])...]
        # multiply all values to be aggregated but not itself
        x = one(T)
        for k in inds
            jk = Base._to_linear_index(src, Tuple(cart_j)..., Tuple(k)...)
            x *= src[jk]
        end
        x /= src[index]
        # apply `op` on `Δsrc[i, k]` and `x`
        Δsrc[index] = op(Δsrc[index], x)
    end
    return nothing
end

function NNlib.∇scatter_src(op::Union{typeof(*),typeof(/)}, Δ, dst,
                            src::AnyCuArray{Tsrc,Nsrc},
                            idx::AnyCuArray{Tidx,Nidx}) where {Tsrc,Tidx,Nsrc,Nidx}
    dims = Nsrc - Nidx
    Δsrc = NNlib.modify_src(op, NNlib.gather(Δ, idx), src)
    rev_idx = NNlib.reverse_indices(idx)
    rev_idx = CuArray(map(CUDA.cudaconvert, rev_idx))

    if dims == 0
        max_idx = length(idx)
        args = op, Δsrc, src, idx, rev_idx, max_idx, Tsrc
    else
        pre_cart_idx = CartesianIndices(axes(src)[1:dims])
        max_dims_idx = length(pre_cart_idx)
        max_idx = max_dims_idx * length(idx)
        args = op, Δsrc, src, idx, rev_idx, pre_cart_idx, max_dims_idx, max_idx, Tsrc
    end

    kernel = @cuda launch=false ∇scatter_src_kernel!(args...)
    config = launch_configuration(kernel.fun; max_threads=256)
    threads = min(max_idx, config.threads)
    blocks = cld(max_idx, threads)
    kernel(args...; threads=threads, blocks=blocks)

    CUDA.unsafe_free!(rev_idx)
    return Δsrc
end


text:
@testset "Continuum mechanics" begin
    @testset "mean/vol/dev" begin
        for T in (Float32, Float64)
            x = rand(SecondOrderTensor{3, T})
            y = rand(SymmetricSecondOrderTensor{3, T})
            # mean
            @test (@inferred mean(x))::T ≈ tr(Array(x)) / 3
            @test (@inferred mean(y))::T ≈ tr(Array(y)) / 3
            # vol/dev for 2nd-order tensors
            @test (@inferred vol(x))::typeof(x) + (@inferred dev(x))::typeof(x) ≈ x
            @test (@inferred vol(y))::typeof(y) + (@inferred dev(y))::typeof(y) ≈ y
            @test (@inferred vol(x))::typeof(x) ⊡ (@inferred dev(x))::typeof(x) ≈ zero(T)  atol = sqrt(eps(T))
            @test (@inferred vol(y))::typeof(y) ⊡ (@inferred dev(y))::typeof(y) ≈ zero(T)  atol = sqrt(eps(T))

            # vol/dev for vectors
            v = rand(Vec{3, T})
            @test (@inferred vol(v)) + (@inferred dev(v)) ≈ v
            @test sum(vol(v)) ≈ sum(v)
            @test dev(vol(v)) ≈ zero(v) atol = sqrt(eps(T))
            @test vol(dev(v)) ≈ zero(v) atol = sqrt(eps(T))

            # vol/dev for 4th-order tensors
            δ = one(SymmetricSecondOrderTensor{3, T})
            I_vol = @einsum (i,j,k,l) -> δ[i,j]*δ[k,l]/3
            I_dev = @einsum (i,j,k,l) -> δ[i,k]*δ[j,l] - δ[i,j]*δ[k,l]/3
            Isym_dev = @einsum (i,j,k,l) -> (δ[i,k]*δ[j,l]+δ[i,l]*δ[j,k])/2 - δ[i,j]*δ[k,l]/3
            @test (@inferred vol(FourthOrderTensor{3}))::FourthOrderTensor{3, Float64}                   ≈ I_vol
            @test (@inferred dev(FourthOrderTensor{3}))::FourthOrderTensor{3, Float64}                   ≈ I_dev
            @test (@inferred vol(FourthOrderTensor{3, T}))::FourthOrderTensor{3, T}                      ≈ I_vol
            @test (@inferred dev(FourthOrderTensor{3, T}))::FourthOrderTensor{3, T}                      ≈ I_dev
            @test (@inferred vol(SymmetricFourthOrderTensor{3}))::SymmetricFourthOrderTensor{3, Float64} ≈ I_vol
            @test (@inferred dev(SymmetricFourthOrderTensor{3}))::SymmetricFourthOrderTensor{3, Float64} ≈ Isym_dev
            @test (@inferred vol(SymmetricFourthOrderTensor{3, T}))::SymmetricFourthOrderTensor{3, T}    ≈ I_vol
            @test (@inferred dev(SymmetricFourthOrderTensor{3, T}))::SymmetricFourthOrderTensor{3, T}    ≈ Isym_dev
            @test (vol(FourthOrderTensor{3}) + dev(FourthOrderTensor{3}))::FourthOrderTensor{3} ≈ one(FourthOrderTensor{3})
            @test (vol(SymmetricFourthOrderTensor{3}) + dev(SymmetricFourthOrderTensor{3}))::SymmetricFourthOrderTensor{3} ≈ one(SymmetricFourthOrderTensor{3})
        end
    end
    @testset "stress invariants" begin
        for T in (Float32, Float64), x in (rand(SecondOrderTensor{3, T}),
                                           rand(SymmetricSecondOrderTensor{3, T}))
            # stress invariants
            I₁, I₂, I₃ = (@inferred stress_invariants(x))::NamedTuple{(:I1,:I2,:I3), NTuple{3, T}}
            @test x^3 - I₁*x^2 + I₂*x - I₃*I ≈ zero(x)  atol = sqrt(eps(T))
            # deviatoric stress invariants
            J₁, J₂, J₃ = (@inferred deviatoric_stress_invariants(x))::NamedTuple{(:J1,:J2,:J3), NTuple{3, T}}
            @test dev(x)^3 - J₁*dev(x)^2 - J₂*dev(x) - J₃*I ≈ zero(x)  atol = sqrt(eps(T))
        end
    end
end


text:
# Autogenerated wrapper script for capnproto_jll for i686-linux-musl
export capnp, libcapnp

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"

# Relative path to `capnp`
const capnp_splitpath = ["bin", "capnp"]

# This will be filled out by __init__() for all products, as it must be done at runtime
capnp_path = ""

# capnp-specific global declaration
function capnp(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        if !isempty(get(ENV, LIBPATH_env, ""))
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', ENV[LIBPATH_env])
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(capnp_path)
    end
end

# Relative path to `libcapnp`
const libcapnp_splitpath = ["lib", "libcapnp.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libcapnp_path = ""

# libcapnp-specific global declaration
# This will be filled out by __init__()
libcapnp_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libcapnp = "libcapnp-0.7.0.so"

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"capnproto")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    global capnp_path = normpath(joinpath(artifact_dir, capnp_splitpath...))

    push!(PATH_list, dirname(capnp_path))
    global libcapnp_path = normpath(joinpath(artifact_dir, libcapnp_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libcapnp_handle = dlopen(libcapnp_path)
    push!(LIBPATH_list, dirname(libcapnp_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
using Yao
using Yao.Intrinsics
################## Benchmarks #################
using BenchmarkTools
v = collect(0:1<<16-1)
orders = randperm(16)
#bres = @benchmark  reorder($v, $orders)
#@code_warntype  collect((reordered_basis(16, orders)))
#bres = @benchmark  collect((reordered_basis(16, orders)))
#bres = @benchmark  collect((reordered_basis(16, orders)))
Pm = pmrand(1<<16)
#@benchmark $Pm |> invorder
Dv = Diagonal(v)
@benchmark $Dv |> invorder