text:
# Define symbol inputs for different variable types
const Parameter = :Parameter

# Extend Base.copy for infinite parameters
function Base.copy(p::ParameterRef, new_model::InfiniteModel)::ParameterRef
    return ParameterRef(new_model, p.index)
end

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

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

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

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

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

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

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

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

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

"""
    build_parameter(_error::Function, set::AbstractInfiniteSet,
                    [num_params::Int = 1; num_supports::Int = 0,
                    supports::Union{Number, Vector{<:Number}} = Number[],
                    independent::Bool = false,
                    sig_fig::Int = 5])::InfOptParameter

Returns a [`InfOptParameter`](@ref) given the appropriate information. This is
analagous to `JuMP.build_variable`. Errors if supports violate the bounds
associated `set`. Also errors if `set` contains a multivariate distribution with
a different dimension than `num_params`. This is meant to primarily serve as a
helper method for [`@infinite_parameter`](@ref).

**Example**
```jldoctest; setup = :(using InfiniteOpt)
julia> build_parameter(error, IntervalSet(0, 3), supports = Vector(0:3))
InfOptParameter{IntervalSet}([0, 3], [0, 1, 2, 3], false)
```
"""
function build_parameter(_error::Function, set::AbstractInfiniteSet,
                         num_params::Int = 1;
                         num_supports::Int = 0,
                         supports::Union{Number, Vector{<:Number}} = Number[],
                         independent::Bool = false,
                         sig_fig::Int = 5,
                         extra_kw_args...)::InfOptParameter
    for (kwarg, _) in extra_kw_args
        _error("Unrecognized keyword argument $kwarg")
    end
    length_supports = length(supports)
    if num_supports == 0 && length_supports != 0
        _check_supports_in_bounds(_error, supports, set)
    elseif num_supports != 0 && length_supports != 0
        @warn("Ignoring num_supports since supports is not empty.")
        _check_supports_in_bounds(_error, supports, set)
    elseif num_supports != 0 && length_supports == 0
        if isa(set, DistributionSet{<:Distributions.MultivariateDistribution})
            _error("Support generation is not available for multivariate " *
                   "distributions.")
        end
        supports = generate_support_values(set, num_supports = num_supports,
                                           sig_fig = sig_fig)
    end
    if isa(set, DistributionSet{<:Distributions.MultivariateDistribution})
        if num_params != length(set.distribution)
            _error("Multivariate distribution dimension must match dimension " *
                   "of parameter.")
        end
    end
    if num_params == 1 || independent # double check num_params for other types of containers
        unique_supports = unique(supports)
        if length(unique_supports) != length(supports)
            @warn("Support points are not unique, eliminating redundant points.")
        end
        return InfOptParameter(set, unique_supports, independent)
    end
    if isa(supports, Number)
        supports = [supports]
    end
    return InfOptParameter(set, supports, independent)
end

# Check the number of supports of one dimension matches the other dimension
# within the same multi-dimensional parameter
function _check_supports_dimensions(model::InfiniteModel, p::InfOptParameter,
    index::Int)
    if haskey(model.params, index - 1) &&
       model.param_to_group_id[index - 1] == model.next_param_id
        if length(p.supports) != length(model.params[index - 1].supports)
            error("Support dimension mismatch. Make sure support dimension of " *
                  "each parameter element match.")
        end
    end
    return
end

"""
    add_parameter(model::InfiniteModel, p::InfOptParameter,
                  [name::String = ""])::ParameterRef

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

**Example**
```jldoctest; setup = :(using InfiniteOpt; model = InfiniteModel())
julia> p = build_parameter(error, IntervalSet(0, 3), supports = Vector(0:3))
InfOptParameter{IntervalSet}([0, 3], [0, 1, 2, 3], false)

julia> param_ref = add_parameter(model, p, "name")
name
```
"""
function add_parameter(model::InfiniteModel, p::InfOptParameter,
                       name::String=""; multi_dim = false,
                       macro_call = false)::ParameterRef
    index = model.next_param_index += 1
    pref = ParameterRef(model, index)
    if !macro_call
        model.next_param_id += 1
    elseif !(p.independent) && multi_dim
        _check_supports_dimensions(model, p, index)
    end
    model.params[JuMP.index(pref)] = p
    model.param_to_group_id[JuMP.index(pref)] = model.next_param_id
    JuMP.set_name(pref, name)
    return pref
end

"""
    used_by_constraint(pref::ParameterRef)::Bool

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

**Example**
```julia-repl
julia> used_by_constraint(t)
true
```
"""
function used_by_constraint(pref::ParameterRef)::Bool
    return haskey(JuMP.owner_model(pref).param_to_constrs, JuMP.index(pref))
end

"""
    used_by_measure(pref::ParameterRef)::Bool

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

**Example**
```julia-repl
julia> used_by_measure(t)
false
```
"""
function used_by_measure(pref::ParameterRef)::Bool
    return haskey(JuMP.owner_model(pref).param_to_meas, JuMP.index(pref))
end

"""
    used_by_variable(pref::ParameterRef)::Bool

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

**Example**
```julia-repl
julia> used_by_variable(t)
true
```
"""
function used_by_variable(pref::ParameterRef)::Bool
    return haskey(JuMP.owner_model(pref).param_to_vars, JuMP.index(pref))
end

"""
    is_used(pref::ParameterRef)::Bool

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

**Example**
```julia-repl
julia> is_used(t)
true
```
"""
function is_used(pref::ParameterRef)::Bool
    return used_by_measure(pref) || used_by_constraint(pref) || used_by_variable(pref)
end

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

## Determine if a tuple element contains a particular parameter
# ParameterRef
function _contains_pref(search_pref::ParameterRef, pref::ParameterRef)::Bool
    return search_pref == pref
end

# SparseAxisArray
function _contains_pref(arr::JuMPC.SparseAxisArray{<:ParameterRef},
                        pref::ParameterRef)::Bool
    return pref in collect(values(arr.data))
end

# Return parameter tuple without a particular parameter and return location of
# where it was
function _remove_parameter(prefs::Tuple, delete_pref::ParameterRef)::Tuple
    for i in eachindex(prefs)
        if _contains_pref(prefs[i], delete_pref)
            if isa(prefs[i], ParameterRef)
                return Tuple(prefs[j] for j = 1:length(prefs) if j != i), (i, )
            else
                for (k, v) in prefs[i].data
                    if v == delete_pref
                        new_dict = filter(x -> x.second != delete_pref,
                                         prefs[i].data)
                        if length(new_dict) != 0
                            pref_list = [prefs...]
                            pref_list[i] = JuMPC.SparseAxisArray(new_dict)
                            return Tuple(pref_list[j] for j = 1:length(prefs)), (i, k)
                        else
                            return Tuple(prefs[j] for j = 1:length(prefs) if j != i), (i, )
                        end
                    end
                end
            end
        end
    end
end

# Used to update infinite variable when one of its parameters is deleted
function _update_infinite_variable(vref::InfiniteVariableRef, new_prefs::Tuple)
    _update_variable_param_refs(vref, new_prefs)
    JuMP.set_name(vref, _root_name(vref))
    if used_by_measure(vref)
        for mindex in JuMP.owner_model(vref).var_to_meas[JuMP.index(vref)]
            JuMP.set_name(MeasureRef(JuMP.owner_model(vref), mindex),
                       _make_meas_name(JuMP.owner_model(vref).measures[mindex]))
        end
    end
    return
end

# Return a parameter value tuple without element at a particular location
function _remove_parameter_values(pref_vals::Tuple, location::Tuple)::Tuple
    # removed parameter was a scalar parameter
    if length(location) == 1
        return Tuple(pref_vals[i] for i in eachindex(pref_vals) if i != location[1])
    # removed parameter was part of an array
    else
        new_dict = filter(x -> x.first != location[2], pref_vals[location[1]].data)
        val_list = [pref_vals...]
        val_list[location[1]] = JuMPC.SparseAxisArray(new_dict)
        return Tuple(val_list[i] for i in eachindex(pref_vals))
    end
end

# Update point variable for which a parameter is deleted
function _update_point_variable(pvref::PointVariableRef, pref_vals::Tuple)
    _update_variable_param_values(pvref, pref_vals)
    # update name if no alias was provided
    if !isa(findfirst(isequal('('), JuMP.name(pvref)), Nothing)
        JuMP.set_name(pvref, "")
    end
    if used_by_measure(pvref)
        for mindex in JuMP.owner_model(pvref).var_to_meas[JuMP.index(pvref)]
            JuMP.set_name(MeasureRef(JuMP.owner_model(pvref), mindex),
                      _make_meas_name(JuMP.owner_model(pvref).measures[mindex]))
        end
    end
    return
end

# Update a reduced variable associated with an infinite variable whose parameter
# was removed
function _update_reduced_variable(vref::ReducedInfiniteVariableRef,
                                  location::Tuple)
    eval_supps = eval_supports(vref)
    # removed parameter was a scalar
    if length(location) == 1
        new_supports = Dict{Int, Union{Number, JuMPC.SparseAxisArray}}()
        for (index, support) in eval_supps
            if index < location[1]
                new_supports[index] = support
            elseif index > location[1]
                new_supports[index - 1] = support
            end
        end
        JuMP.owner_model(vref).reduced_info[JuMP.index(vref)] = ReducedInfiniteInfo(infinite_variable_ref(vref), new_supports)
    # removed parameter was part of an array and was reduced previously
    elseif haskey(eval_supps, location[1])
        new_dict = filter(x -> x.first != location[2], eval_supps[location[1]].data)
        new_supports = copy(eval_supps)
        new_supports[location[1]] = JuMPC.SparseAxisArray(new_dict)
        JuMP.owner_model(vref).reduced_info[JuMP.index(vref)] = ReducedInfiniteInfo(infinite_variable_ref(vref), new_supports)
    end
    if used_by_measure(vref)
        for mindex in JuMP.owner_model(vref).reduced_to_meas[JuMP.index(vref)]
            JuMP.set_name(MeasureRef(JuMP.owner_model(vref), mindex),
                       _make_meas_name(JuMP.owner_model(vref).measures[mindex]))
        end
    end
    return
end

"""
    JuMP.delete(model::InfiniteModel, pref::ParameterRef)

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

**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
 g(0.5, x) == 0
 t in [0, 6]
 x in [0, 1]

julia> delete(model, x)

julia> print(model)
Min measure(g(t)*t) + z
Subject to
 g(t) + z >= 42.0
 g(0.5) == 0
 z >= 0.0
 t in [0, 6]
```
"""
function JuMP.delete(model::InfiniteModel, pref::ParameterRef)
    @assert JuMP.is_valid(model, pref) "Parameter reference is invalid."
    # update optimizer model status
    if is_used(pref)
        set_optimizer_model_ready(model, false)
    end
    # update measures
    if used_by_measure(pref)
        # ensure deletion is okay (pref isn't used by measure data)
        for mindex in model.param_to_meas[JuMP.index(pref)]
            _check_param_in_data(pref, model.measures[mindex].data)
        end
        # delete dependence of measures on pref
        for mindex in model.param_to_meas[JuMP.index(pref)]
            if isa(model.measures[mindex].func, ParameterRef)
                model.measures[mindex] = Measure(zero(JuMP.AffExpr),
                                                 model.measures[mindex].data)
            else
                _remove_variable(model.measures[mindex].func, pref)
            end
            JuMP.set_name(MeasureRef(model, mindex),
                          _make_meas_name(model.measures[mindex]))
        end
        # delete mapping
        delete!(model.param_to_meas, JuMP.index(pref))
    end
    # update variables
    if used_by_variable(pref)
        # update infinite variables that depend on pref
        for vindex in model.param_to_vars[JuMP.index(pref)]
            # find location of parameter in storage tuple
            prefs, location = _remove_parameter(model.vars[vindex].parameter_refs,
                                                pref)
            vref = InfiniteVariableRef(model, vindex)
            # remove the parameter dependence
            _update_infinite_variable(vref, prefs)
            # update any point variables that depend on vref accordingly
            if used_by_point_variable(vref)
                for pindex in model.infinite_to_points[vindex]
                    pvref = PointVariableRef(model, pindex)
                    pref_vals = _remove_parameter_values(parameter_values(pvref),
                                                         location)
                    _update_point_variable(pvref, pref_vals)
                end
            end
            # update any reduced variables that depend on vref accordingly
            if used_by_reduced_variable(vref)
                for rindex in model.infinite_to_reduced[vindex]
                    _update_reduced_variable(ReducedInfiniteVariableRef(model,
                                                              rindex), location)
                end
            end
        end
        # delete mapping
        delete!(model.param_to_vars, JuMP.index(pref))
    end
    # update constraints
    if used_by_constraint(pref)
        # update constraints in mapping to remove the parameter
        for cindex in model.param_to_constrs[JuMP.index(pref)]
            if isa(model.constrs[cindex].func, ParameterRef)
                model.constrs[cindex] = JuMP.ScalarConstraint(zero(JuMP.AffExpr),
                                                      model.constrs[cindex].set)
            else
                _remove_variable(model.constrs[cindex].func, pref)
            end
        end
        # delete mapping
        delete!(model.param_to_constrs, JuMP.index(pref))
    end
    # delete parameter information stored in model
    delete!(model.params, JuMP.index(pref))
    delete!(model.param_to_name, JuMP.index(pref))
    delete!(model.param_to_group_id, JuMP.index(pref))
    return
end

"""
    JuMP.is_valid(model::InfiniteModel, pref::ParameterRef)::Bool

Extend the [`JuMP.is_valid`](@ref JuMP.is_valid(::JuMP.Model, ::JuMP.VariableRef))
function to accomodate infinite parameters.
Returns true if the `InfiniteModel` stored in `pref` matches `model` and if
the parameter index is used by `model`. It returns false otherwise.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1]))
julia> is_valid(model, t)
true
```
"""
function JuMP.is_valid(model::InfiniteModel, pref::ParameterRef)::Bool
    check1 = model === JuMP.owner_model(pref)
    check2 = JuMP.index(pref) in keys(model.params)
    return check1 && check2
end

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

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

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1]))
julia> name(t)
"t"
```
"""
function JuMP.name(pref::ParameterRef)::String
    return JuMP.owner_model(pref).param_to_name[JuMP.index(pref)]
end

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

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

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

julia> name(t)
"time"
```
"""
function JuMP.set_name(pref::ParameterRef, name::String)
    JuMP.owner_model(pref).param_to_name[JuMP.index(pref)] = name
    JuMP.owner_model(pref).name_to_param = nothing
    return
end

"""
    num_parameters(model::InfiniteModel)::Int

Return the number of infinite parameters currently present in `model`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1]))
julia> num_parameters(model)
1
```
"""
function num_parameters(model::InfiniteModel)::Int
    return length(model.params)
end

# Internal functions
_parameter_set(pref::ParameterRef) = JuMP.owner_model(pref).params[JuMP.index(pref)].set
_parameter_supports(pref::ParameterRef) = JuMP.owner_model(pref).params[JuMP.index(pref)].supports
function _update_parameter_set(pref::ParameterRef, set::AbstractInfiniteSet)
    supports = JuMP.owner_model(pref).params[JuMP.index(pref)].supports
    independent = JuMP.owner_model(pref).params[JuMP.index(pref)].independent
    JuMP.owner_model(pref).params[JuMP.index(pref)] = InfOptParameter(set, supports, independent)
    if is_used(pref)
        set_optimizer_model_ready(JuMP.owner_model(pref), false)
    end
    return
end
function _update_parameter_supports(pref::ParameterRef, supports::Vector{<:Number})
    set = JuMP.owner_model(pref).params[JuMP.index(pref)].set
    independent = JuMP.owner_model(pref).params[JuMP.index(pref)].independent
    JuMP.owner_model(pref).params[JuMP.index(pref)] = InfOptParameter(set, supports, independent)
    if is_used(pref)
        set_optimizer_model_ready(JuMP.owner_model(pref), false)
    end
    return
end

"""
    infinite_set(pref::ParameterRef)::AbstractInfiniteSet

Return the infinite set associated with `pref`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1]))
julia> infinite_set(t)
[0, 1]
```
"""
function infinite_set(pref::ParameterRef)::AbstractInfiniteSet
    return _parameter_set(pref)
end

"""
    set_infinite_set(pref::ParameterRef, set::AbstractInfiniteSet)

Specify the infinite set of `pref`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1]))
julia> set_infinite_set(t, IntervalSet(0, 2))

julia> infinite_set(t)
[0, 2]
```
"""
function set_infinite_set(pref::ParameterRef, set::AbstractInfiniteSet)
    _update_parameter_set(pref, set)
    return
end

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

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

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1]))
julia> has_lower_bound(t)
true
```
"""
function JuMP.has_lower_bound(pref::ParameterRef)::Bool
    set = _parameter_set(pref)
    return JuMP.has_lower_bound(set)
end

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

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

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1]))
julia> lower_bound(t)
0.0
```
"""
function JuMP.lower_bound(pref::ParameterRef)::Number
    set = _parameter_set(pref)
    if !JuMP.has_lower_bound(pref)
        error("Parameter $(pref) does not have a lower bound.")
    end
    return JuMP.lower_bound(set)
end

"""
    JuMP.set_lower_bound(pref::ParameterRef, lower::Number)

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

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1]))
julia> set_lower_bound(t, -1)

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

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

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

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1]))
julia> has_upper_bound(t)
true
```
"""
function JuMP.has_upper_bound(pref::ParameterRef)::Bool
    set = _parameter_set(pref)
    return JuMP.has_upper_bound(set)
end

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

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

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1]))
julia> upper_bound(t)
1.0
```
"""
function JuMP.upper_bound(pref::ParameterRef)::Number
    set = _parameter_set(pref)
    if !JuMP.has_upper_bound(pref)
        error("Parameter $(pref) does not have a upper bound.")
    end
    return JuMP.upper_bound(set)
end

"""
    JuMP.set_upper_bound(pref::ParameterRef, lower::Number)

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

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1]))
julia> set_upper_bound(t, 2)

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

"""
    num_supports(pref::ParameterRef)::Int

Return the number of support points associated with `pref`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1], supports = [0, 1]))
julia> num_supports(t)
2
```
"""
function num_supports(pref::ParameterRef)::Int
    return length(_parameter_supports(pref))
end

"""
    has_supports(pref::ParameterRef)::Bool

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

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1], supports = [0, 1]))
julia> has_supports(t)
true
```
"""
has_supports(pref::ParameterRef)::Bool = num_supports(pref) > 0

"""
    supports(pref::ParameterRef)::Vector

Return the support points associated with `pref`. Errors if there are no
supports.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1], supports = [0, 1]))
julia> supports(t)
2-element Array{Int64,1}:
 0
 1
```
"""
function supports(pref::ParameterRef)::Vector
    !has_supports(pref) && error("Parameter $pref does not have supports.")
    return _parameter_supports(pref)
end

"""
    supports(prefs::AbstractArray{<:ParameterRef})::Vector

Return the support points associated with an array of `prefs` formatted as a
vector of SparseAxisArrays following the format of the input array. If the
parameters are not independent then the supports of each parameter are simply
spliced together. Alternatively can call `supports.` to more efficiently obtain
an array of the same input format whose parameter references have been replaced
with their supports. Errors if all the parameter references do not have the same
group ID number (were intialized together as an array) or if the nonindependent
parameters have support vectors of different lengths. If the parameters are
independent then all the unique combinations are identified and returned as
supports. Warning this operation is computationally expensive if there exist a
large number of combinations.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel())
julia> x = @infinite_parameter(model, [i = 1:2], set = IntervalSet(-1, 1),
                               base_name = "x", independent = true)
2-element Array{ParameterRef,1}:
 x[1]
 x[2]

julia> for i = 1:length(x)
           set_supports(x[i], [-1, 1])
       end

julia> supports(x)
4-element Array{JuMP.Containers.SparseAxisArray,1}:
   [2]  =  -1
  [1]  =  -1
   [2]  =  1
  [1]  =  -1
   [2]  =  -1
  [1]  =  1
   [2]  =  1
  [1]  =  1
```
"""
function supports(prefs::AbstractArray{<:ParameterRef})::Vector
    prefs = convert(JuMPC.SparseAxisArray, prefs)
    !_only_one_group(prefs) && error("Array contains parameters from multiple" *
                                     " groups.")
    for (k, pref) in prefs.data
        !has_supports(pref) && error("Parameter $pref does not have supports.")
    end
    lengths = [num_supports(pref) for (k, pref) in prefs.data]
    if !is_independent(collect(values(prefs.data))[1])
        length(unique(lengths)) != 1 && error("Each nonindependent parameter " *
                                              "must have the same number of " *
                                              "support points.")
        support_list = Vector{JuMPC.SparseAxisArray}(undef, lengths[1])
        for i in eachindex(support_list)
            support_list[i] = JuMPC.SparseAxisArray(Dict(k => supports(pref)[i] for (k, pref) in prefs.data))
        end
        # unique-nize support_list using the dictionary unique method
        # TODO: make it work with JuMPC.SparseAxisArray
        dict_list = unique([arr.data for arr in support_list])
        support_list = [JuMPC.SparseAxisArray(dict) for dict in dict_list]
    else
        all_keys = collect(keys(prefs))
        all_supports = [supports(pref) for (k, pref) in prefs.data]
        support_list = Vector{JuMPC.SparseAxisArray}(undef, prod(lengths))
        counter = 1
        for combo in Iterators.product(all_supports...)
            support_list[counter] = JuMPC.SparseAxisArray(Dict(all_keys[i] => combo[i] for i in eachindex(combo)))
            counter += 1
        end
    end
    return support_list
end

"""
    set_supports(pref::ParameterRef, supports::Vector{<:Number}; [force = false])

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

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

julia> supports(t)
2-element Array{Int64,1}:
 0
 1
```
"""
function set_supports(pref::ParameterRef, supports::Vector{<:Number};
                      force = false)
    set = _parameter_set(pref)
    _check_supports_in_bounds(error, supports, set)
    if has_supports(pref) && !force
        error("Unable set supports for $pref since it already has supports." *
              " Consider using `add_supports` or use set `force = true` to " *
              "overwrite the existing supports.")
    end
    if !(is_independent(pref)) &&
       sum(values(pref.model.param_to_group_id) .== group_id(pref)) > 1
        _update_parameter_supports(pref, supports)
    else
        unique_supports = unique(supports)
        if length(unique_supports) != length(supports)
            @warn("Support points are not unique, eliminating redundant points.")
        end
        _update_parameter_supports(pref, unique_supports)
    end
    return
end

"""
    add_supports(pref::ParameterRef, supports::Union{Number, Vector{<:Number}})

Add additional support points for `pref`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1], supports = [0, 1]))
julia> add_supports(t, 0.5)

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

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

julia> supports(t)
4-element Array{Float64,1}:
 0.0
 1.0
 0.5
 0.25
```
"""
function add_supports(pref::ParameterRef, supports::Union{Number,
                                                          Vector{<:Number}})
    set = _parameter_set(pref)
    _check_supports_in_bounds(error, supports, set)
    current_supports = _parameter_supports(pref)
    if !(is_independent(pref)) &&
       sum(values(pref.model.param_to_group_id) .== group_id(pref)) > 1
        new_supports = [current_supports; supports]
    else
        new_supports = unique([current_supports; supports])
    end
    _update_parameter_supports(pref, new_supports)
    return
end

"""
    delete_supports(pref::ParameterRef)

Delete the support points for `pref`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1], supports = [0, 1]))
julia> delete_supports(t)

julia> supports(t)
ERROR: Parameter t does not have supports.
```
"""
function delete_supports(pref::ParameterRef)
    _update_parameter_supports(pref, Int[])
    return
end

"""
    is_finite_parameter(pref::ParameterRef)::Bool

Return a `Bool` indicating if `pref` is a finite parameter.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @finite_parameter(model, cost, 42))
julia> is_finite_parameter(cost)
true
```
"""
function is_finite_parameter(pref::ParameterRef)::Bool
    set = infinite_set(pref)
    if isa(set, IntervalSet) && set.lower_bound == set.upper_bound
        return true
    end
    return false
end

"""
    JuMP.value(pref::ParameterRef)::Number

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

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @finite_parameter(model, cost, 42))
julia> value(cost)
42
```
"""
function JuMP.value(pref::ParameterRef)::Number
    !is_finite_parameter(pref) && error("$pref is an infinite parameter.")
    return supports(pref)[1]
end

"""
    JuMP.set_value(pref::ParameterRef, value::Number)

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

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @finite_parameter(model, cost, 42))
julia> set_value(cost, 27)

julia> value(cost)
27
```
"""
function JuMP.set_value(pref::ParameterRef, value::Number)
    !is_finite_parameter(pref) && error("$pref is an infinite parameter.")
    set_infinite_set(pref, IntervalSet(value, value))
    set_supports(pref, [value], force = true)
    return
end

"""
    fill_in_supports!(model::InfiniteModel; [num_supports::Int = 10,
                      sig_fig::Int = 5])

Automatically generate support points for all infinite parameters in model. User
can specify the number of significant figures kept after decimal point for the
auto-generated supports wtih `sig_fig`. This calls
[`fill_in_supports!`](@ref fill_in_supports!(::ParameterRef)) for each parameter
in the model. See [`fill_in_supports!`](@ref fill_in_supports!(::ParameterRef))
for more information. Errors if one of the infinite set types is unrecognized.

**Example**
```jldoctest; setup = :(using InfiniteOpt; model = InfiniteModel(); @infinite_parameter(model, 0 <= x <= 1);)
julia> fill_in_supports!(model, num_supports = 4, sig_fig = 3)

julia> supports(x)
4-element Array{Number,1}:
 0.0
 0.333
 0.667
 1.0
```
"""
function fill_in_supports!(model::InfiniteModel; num_supports::Int = 10,
                           sig_fig::Int = 5)
    for key in keys(model.params)
        pref = ParameterRef(model, key)
        fill_in_supports!(pref, num_supports = num_supports, sig_fig = sig_fig)
    end
    return
end

"""
    fill_in_supports!(pref::ParameterRef; [num_supports::Int = 10,
                                           sig_fig::Int = 5])

Automatically generate support points for a particular infinite parameter `pref`.
Generating `num_supports` for the parameter. The supports are generated uniformly
if the underlying infinite set is an `IntervalSet` or they are generating randomly
accordingly to the distribution if the set is a `DistributionSet`.
User can specify the number of digits kept after decimal point for the
auto-generated supports wtih `sig_fig`. 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**
```jldoctest; setup = :(using InfiniteOpt; model = InfiniteModel(); @infinite_parameter(model, 0 <= x <= 1);)
julia> fill_in_supports!(x, num_supports = 4, sig_fig = 3)

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

```
"""
function fill_in_supports!(pref::ParameterRef; num_supports::Int = 10,
                           sig_fig::Int = 5)
    p = JuMP.owner_model(pref).params[JuMP.index(pref)]
    if length(p.supports) == 0
        generate_and_add_supports!(pref, p.set, num_supports = num_supports,
                                   sig_fig = sig_fig)
    end
    return
end

"""
    generate_and_add_supports!(pref::ParameterRef, set::AbstractInfiniteSet;
                               [num_supports::Int = 10, sig_fig::Int = 5])

Generate supports for `pref` via [`generate_support_values`](@ref) and add them
to `pref`. This is intended as an extendable internal method for
[`fill_in_supports!`](@ref fill_in_supports!(::ParameterRef)). Note that if
`pref` is part of a `DistributionSet` that features a multivariate distribution,
all the associated parameters with `pref` will also have supports added to them.
Most extensions that empoy 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 single infinite parameter (e.g., how we enable multivariate
distribution sets). Errors if the infinite set type is not recognized.
"""
function generate_and_add_supports!(pref::ParameterRef, set::AbstractInfiniteSet;
                                    num_supports::Int = 10, sig_fig::Int = 5)
    add_supports(pref, generate_support_values(set, num_supports = num_supports,
                                               sig_fig = sig_fig))
    return
end

# Multivariate distribution sets
function generate_and_add_supports!(pref::ParameterRef,
                                    set::DistributionSet{<:Distributions.MultivariateDistribution};
                                    num_supports::Int = 10, sig_fig::Int = 5)
    pref_group_id = group_id(pref)
    model = JuMP.owner_model(pref)
    associated_p_index = sort([i for i in 1:length(model.params)
                               if model.param_to_group_id[i] == pref_group_id])
    new_supports = generate_support_values(set, num_supports = num_supports, sig_fig = sig_fig)

    for i in 1:length(associated_p_index)
        pref_i = ParameterRef(model, associated_p_index[i])
        add_supports(pref_i, new_supports[i, :])
    end
    return
end

"""
    group_id(pref::ParameterRef)::Int

Return the group ID number for `pref`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1], supports = [0, 1]))
julia> group_id(t)
1
```
"""
function group_id(pref::ParameterRef)::Int
    return JuMP.owner_model(pref).param_to_group_id[JuMP.index(pref)]
end

"""
    group_id(prefs::AbstractArray{<:ParameterRef})::Int

Return the group ID number for a group of `prefs`. Error if contains multiple
groups.

**Example**
```julia-repl
julia> group_id([x[1], x[2]])
2
```
"""
function group_id(prefs::AbstractArray{<:ParameterRef})::Int
    groups = group_id.(prefs)
    length(unique(groups)) != 1 && error("Array contains parameters from " *
                                         "multiple groups.")
    return first(groups)
end

"""
    is_independent(pref::ParameterRef)::Bool

Returns true for `pref` if it is independent or false otherwise.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1], supports = [0, 1]))
julia> is_independent(t)
false
```
"""
function is_independent(pref::ParameterRef)::Bool
    return JuMP.owner_model(pref).params[JuMP.index(pref)].independent
end

"""
    set_independent(pref::ParameterRef)

Specify that `pref` be independent.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1], supports = [0, 1]))
julia> set_independent(t)

julia> is_independent(t)
true
```
"""
function set_independent(pref::ParameterRef)
    old_param = JuMP.owner_model(pref).params[JuMP.index(pref)]
    new_param = InfOptParameter(old_param.set, old_param.supports, true)
    JuMP.owner_model(pref).params[JuMP.index(pref)] = new_param
    return
end

"""
    unset_independent(pref::ParameterRef)

Specify that `pref` be not independent.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1], supports = [0, 1]))
julia> unset_independent(t)

julia> is_independent(t)
false
```
"""
function unset_independent(pref::ParameterRef)
    old_param = JuMP.owner_model(pref).params[JuMP.index(pref)]
    new_param = InfOptParameter(old_param.set, old_param.supports, false)
    JuMP.owner_model(pref).params[JuMP.index(pref)] = new_param
    return
end

"""
    parameter_by_name(model::InfiniteModel, name::String)::Union{ParameterRef,
                                                                 Nothing}

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

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP; model = InfiniteModel(); @infinite_parameter(model, t in [0, 1], supports = [0, 1]))
julia> parameter_by_name(model, "t")
t
```
"""
function parameter_by_name(model::InfiniteModel,
                           name::String)::Union{ParameterRef, Nothing}
    if model.name_to_param === nothing
        # Inspired from MOI/src/Utilities/model.jl
        model.name_to_param = Dict{String, Int}()
        for (param, param_name) in model.param_to_name
            if haskey(model.name_to_param, param_name)
                # -1 is a special value that means this string does not map to
                # a unique variable name.
                model.name_to_param[param_name] = -1
            else
                model.name_to_param[param_name] = param
            end
        end
    end
    index = get(model.name_to_param, name, nothing)
    if index isa Nothing
        return
    elseif index == -1
        error("Multiple parameters have the name $name.")
    else
        return ParameterRef(model, index)
    end
end

"""
    all_parameters(model::InfiniteModel)::Vector{ParameterRef}

Return all of the infinite parameter references currently in `model`.

**Example**
```julia-repl
julia> all_parameters(model)
3-element Array{ParameterRef,1}:
 t
 x[1]
 x[2]
```
"""
function all_parameters(model::InfiniteModel)::Vector{ParameterRef}
    pref_list = Vector{ParameterRef}(undef, num_parameters(model))
    indexes = sort([index for index in keys(model.params)])
    counter = 1
    for index in indexes
        pref_list[counter] = ParameterRef(model, index)
        counter += 1
    end
    return pref_list
end

## Define functions to extract the names of parameters
# Extract the root name of a parameter reference
function _root_name(pref::ParameterRef)::String
    name = JuMP.name(pref)
    first_bracket = findfirst(isequal('['), name)
    if first_bracket == nothing
        return name
    else
        # Hacky fix to handle invalid Unicode
        try
            return name[1:first_bracket-1]
        catch
            return name[1:first_bracket-2]
        end
    end
end

# Return the root names of a tuple parameter of references
function _root_names(prefs::Tuple)::Tuple
    return _root_name.(first.(prefs))
end

## Internal functions for group checking
# Return group id of ParameterRef
function _group(pref::ParameterRef)::Int
    return group_id(pref)
end

# Return the group if of the first element in an array (assuming all same)
function _group(arr::AbstractArray{<:ParameterRef})::Int
    return group_id(first(arr))
end

# Return true if SparseAxisArray only has one group
function _only_one_group(arr::JuMPC.SparseAxisArray{<:ParameterRef})::Bool
    return length(unique(group_id.(arr))) == 1
end

# Return true to have one group ID since is singular
_only_one_group(pref::ParameterRef)::Bool = true

## Internal function for extracting parameter references
# Return a vector of parameter references from a tuple of references
function _list_parameter_refs(prefs::Tuple)
    list = ParameterRef[]
    for pref in prefs
        if isa(pref, ParameterRef)
            push!(list, pref)
        else
            for k in keys(pref.data)
                push!(list, pref[k])
            end
        end
    end
    return list
end


text:
function [CG]=DerivativeMatrix(CG)
OrdPolyX=CG.OrdPolyX;
CG.DSX=zeros(OrdPolyX+1,OrdPolyX+1);
for i=1:OrdPolyX+1
  for j=1:OrdPolyX+1
    CG.DSX(i,j)=DLagrange(CG.xwX(i),CG.xwX,j);
  end
end
CG.DWX=-inv(diag(CG.wX))*CG.DSX'*diag(CG.wX);

OrdPolyY=CG.OrdPolyY;
CG.DSY=zeros(OrdPolyY+1,OrdPolyY+1);
for i=1:OrdPolyY+1
  for j=1:OrdPolyY+1
    CG.DSY(i,j)=DLagrange(CG.xwY(i),CG.xwY,j);
  end
end
CG.DWY=-inv(diag(CG.wY))*CG.DSY'*diag(CG.wY);

if strcmp(CG.FeX,'LG')
  OrdPolyX=CG.OrdPolyX;
  CG.DSXLGL=zeros(OrdPolyX+2,OrdPolyX+1);
  CG.IntXLGL=zeros(OrdPolyX+2,OrdPolyX+1);
  CG.IntXLG=zeros(OrdPolyX+1,OrdPolyX+2);
  [CG.wXLGL,CG.xwXLGL]=GaussLobattoQuad(CG.OrdPolyX+1);
  for i=1:OrdPolyX+2
    for j=1:OrdPolyX+1
      CG.DSXLGL(i,j)=DLagrange(CG.xwXLGL(i),CG.xwX,j);
      CG.IntXLGL(i,j)=Lagrange(CG.xwXLGL(i),CG.xwX,j);
      CG.IntXLG(j,i)=Lagrange(CG.xwX(j),CG.xwXLGL,i);
    end
  end
  %CG.IntXLG=inv(diag(CG.wX))*CG.IntXLG*diag(CG.wXLGL);
  CG.DWXLGL=-inv(diag(CG.wX))*CG.DSXLGL'*diag(CG.wXLGL);
end
if strcmp(CG.FeX,'LG') && strcmp(CG.FeY,'LGL')
  DSXLGL=zeros(OrdPolyX+2,OrdPolyX+1);
  for i=1:OrdPolyX+2
    for j=1:OrdPolyX+1
      DSXLGL(i,j)=DLagrange(CG.xwXLGL(i),CG.xwX,j);
    end
  end
  CG.DWXLG=-inv(diag(CG.wX))*DSXLGL'*diag(CG.wXLGL);
end
if strcmp(CG.FeX,'LGL')
  [CG.wXLG,CG.xwXLG]=GaussLegendreQuad(CG.OrdPolyX-1);
  CG.IntXLG=zeros(CG.OrdPolyX,CG.OrdPolyX+1);
  for i=1:CG.OrdPolyX
    for j=1:CG.OrdPolyX+1
      CG.IntXLG(i,j)=Lagrange(CG.xwXLG(i),CG.xwX,j);
      CG.IntXLGL(j,i)=Lagrange(CG.xwX(j),CG.xwXLG,i);
    end
  end
  CG.IntXLGL=eye(CG.OrdPolyX+1);
  CG.IntXLG=eye(CG.OrdPolyX+1);
end

if strcmp(CG.FeY,'LGL')
  [CG.wYLG,CG.xwYLG]=GaussLegendreQuad(CG.OrdPolyY-1);
  CG.IntYLG=zeros(CG.OrdPolyY,CG.OrdPolyY+1);
  for i=1:CG.OrdPolyY
    for j=1:CG.OrdPolyY+1
      CG.IntYLG(i,j)=Lagrange(CG.xwYLG(i),CG.xwY,j);
      CG.IntYLGL(j,i)=Lagrange(CG.xwY(j),CG.xwYLG,i);
    end
  end
  CG.IntYLGL=eye(CG.OrdPolyY+1);
  CG.IntYLG=eye(CG.OrdPolyY+1);
end

if strcmp(CG.FeY,'LG')
  OrdPolyY=CG.OrdPolyY;
  CG.DSYLGL=zeros(OrdPolyY+2,OrdPolyY+1);
  CG.IntYLGL=zeros(OrdPolyY+2,OrdPolyY+1);
  CG.IntYLG=zeros(OrdPolyY+1,OrdPolyY+2);
  [CG.wYLGL,CG.xwYLGL]=GaussLobattoQuad(CG.OrdPolyY+1);
  for i=1:OrdPolyY+2
    for j=1:OrdPolyY+1
      CG.DSYLGL(i,j)=DLagrange(CG.xwYLGL(i),CG.xwY,j);
      CG.IntYLGL(i,j)=Lagrange(CG.xwYLGL(i),CG.xwY,j);
      CG.IntYLG(j,i)=Lagrange(CG.xwY(j),CG.xwYLGL,i);
    end
  end
  %CG.IntYLG=inv(diag(CG.wY))*CG.IntYLG*diag(CG.wYLGL);
  CG.DWYLGL=-inv(diag(CG.wY))*CG.DSYLGL'*diag(CG.wYLGL);
end
end

ext: jl


text:
# Test statistics related functions

using NumericExtensions
using Base.Test

# data

x = randn(3, 4)
y = randn(3, 4)

x1 = randn(6)
y1 = randn(6)

x2 = randn(5, 6)
y2 = randn(5, 6)

x3 = randn(3, 4, 5)
y3 = randn(3, 4, 5)

x4 = randn(3, 4, 5, 3)
y4 = randn(3, 4, 5, 3)

# varm

safe_varm(x, u) = sum(abs2(x - u)) / (length(x) - 1)
safe_varm(x, u, d::Int) = sum(abs2(x .- u), d) / (size(x, d) - 1)

@test_approx_eq varm(x, mean(y)) safe_varm(x, mean(y))
@test_approx_eq varm(x1, mean(y1,1), 1) safe_varm(x1, mean(y1,1), 1)

@test_approx_eq varm(x2, mean(y2,1), 1) safe_varm(x2, mean(y2,1), 1)
@test_approx_eq varm(x2, mean(y2,2), 2) safe_varm(x2, mean(y2,2), 2)

@test_approx_eq varm(x3, mean(y3,1), 1) safe_varm(x3, mean(y3,1), 1)
@test_approx_eq varm(x3, mean(y3,2), 2) safe_varm(x3, mean(y3,2), 2)
@test_approx_eq varm(x3, mean(y3,3), 3) safe_varm(x3, mean(y3,3), 3)

@test_approx_eq varm(x4, mean(y4,1), 1) safe_varm(x4, mean(y4,1), 1)
@test_approx_eq varm(x4, mean(y4,2), 2) safe_varm(x4, mean(y4,2), 2)
@test_approx_eq varm(x4, mean(y4,3), 3) safe_varm(x4, mean(y4,3), 3)
@test_approx_eq varm(x4, mean(y4,4), 4) safe_varm(x4, mean(y4,4), 4)

# stdm

safe_stdm(x, u) = sqrt(safe_varm(x, u))
safe_stdm(x, u, d::Int) = sqrt(safe_varm(x, u, d))

@test_approx_eq stdm(x, mean(y)) safe_stdm(x, mean(y))
@test_approx_eq stdm(x1, mean(y1,1), 1) safe_stdm(x1, mean(y1,1), 1)

@test_approx_eq stdm(x2, mean(y2,1), 1) safe_stdm(x2, mean(y2,1), 1)
@test_approx_eq stdm(x2, mean(y2,2), 2) safe_stdm(x2, mean(y2,2), 2)

@test_approx_eq stdm(x3, mean(y3,1), 1) safe_stdm(x3, mean(y3,1), 1)
@test_approx_eq stdm(x3, mean(y3,2), 2) safe_stdm(x3, mean(y3,2), 2)
@test_approx_eq stdm(x3, mean(y3,3), 3) safe_stdm(x3, mean(y3,3), 3)

@test_approx_eq stdm(x4, mean(y4,1), 1) safe_stdm(x4, mean(y4,1), 1)
@test_approx_eq stdm(x4, mean(y4,2), 2) safe_stdm(x4, mean(y4,2), 2)
@test_approx_eq stdm(x4, mean(y4,3), 3) safe_stdm(x4, mean(y4,3), 3)
@test_approx_eq stdm(x4, mean(y4,4), 4) safe_stdm(x4, mean(y4,4), 4)

# var

safe_var(x) = safe_varm(x, mean(x))
safe_var(x, d::Int) = safe_varm(x, mean(x, d), d)

@test_approx_eq var(x) safe_var(x)
@test_approx_eq var(x1, 1) safe_var(x1, 1)

@test_approx_eq var(x2, 1) safe_var(x2, 1)
@test_approx_eq var(x2, 2) safe_var(x2, 2)

@test_approx_eq var(x3, 1) safe_var(x3, 1)
@test_approx_eq var(x3, 2) safe_var(x3, 2)
@test_approx_eq var(x3, 3) safe_var(x3, 3)

@test_approx_eq var(x4, 1) safe_var(x4, 1)
@test_approx_eq var(x4, 2) safe_var(x4, 2)
@test_approx_eq var(x4, 3) safe_var(x4, 3)
@test_approx_eq var(x4, 4) safe_var(x4, 4)

r = zeros(size(x2, 2)); var!(r, x2, 1)
@test_approx_eq r vec(var(x2, 1))

# std

safe_std(x) = sqrt(var(x))
safe_std(x, d) = sqrt(var(x, d))

@test_approx_eq std(x) safe_std(x)
@test_approx_eq std(x1, 1) safe_std(x1, 1)

@test_approx_eq std(x2, 1) safe_std(x2, 1)
@test_approx_eq std(x2, 2) safe_std(x2, 2)

@test_approx_eq std(x3, 1) safe_std(x3, 1)
@test_approx_eq std(x3, 2) safe_std(x3, 2)
@test_approx_eq std(x3, 3) safe_std(x3, 3)

@test_approx_eq std(x4, 1) safe_std(x4, 1)
@test_approx_eq std(x4, 2) safe_std(x4, 2)
@test_approx_eq std(x4, 3) safe_std(x4, 3)
@test_approx_eq std(x4, 4) safe_std(x4, 4)

r = zeros(size(x2, 2)); std!(r, x2, 1)
@test_approx_eq r vec(std(x2, 1))

# logsumexp

safe_logsumexp(x) = log(sum(exp(x)))
safe_logsumexp(x, d::Int) = log(sum(exp(x), d))

@test_approx_eq logsumexp(x) safe_logsumexp(x)
@test_approx_eq logsumexp(x1, 1) safe_logsumexp(x1, 1)
@test_approx_eq logsumexp(x2, 1) safe_logsumexp(x2, 1)
@test_approx_eq logsumexp(x2, 2) safe_logsumexp(x2, 2)
@test_approx_eq logsumexp(x3, 1) safe_logsumexp(x3, 1)
@test_approx_eq logsumexp(x3, 2) safe_logsumexp(x3, 2)
@test_approx_eq logsumexp(x3, 3) safe_logsumexp(x3, 3)

r = zeros(size(x2, 2)); logsumexp!(r, x2, 1)
@test_approx_eq r vec(logsumexp(x2, 1))

@test_approx_eq logsumexp(x + 1000.) logsumexp(x) + 1000.
@test_approx_eq logsumexp(x2 + 1000., 1) logsumexp(x2, 1) + 1000.
@test_approx_eq logsumexp(x2 + 1000., 2) logsumexp(x2, 2) + 1000.
@test_approx_eq logsumexp(x3 + 1000., 2) logsumexp(x3, 2) + 1000.

# softmax

safe_softmax(x) = exp(x) / sum(exp(x))
safe_softmax(x, d::Int) = exp(x) ./ sum(exp(x), d)

@test_approx_eq softmax(x) safe_softmax(x)
@test_approx_eq softmax(x1, 1) safe_softmax(x1, 1)
@test_approx_eq softmax(x2, 1) safe_softmax(x2, 1)
@test_approx_eq softmax(x2, 2) safe_softmax(x2, 2)
@test_approx_eq softmax(x3, 1) safe_softmax(x3, 1)
@test_approx_eq softmax(x3, 2) safe_softmax(x3, 2)
@test_approx_eq softmax(x3, 3) safe_softmax(x3, 3)

r = zeros(size(x2)); softmax!(r, x2, 1)
@test_approx_eq r softmax(x2, 1)

@test_approx_eq softmax(x + 1000.) softmax(x)
@test_approx_eq softmax(x2 + 1000., 1) softmax(x2, 1)
@test_approx_eq softmax(x2 + 1000., 2) softmax(x2, 2)
@test_approx_eq softmax(x3 + 1000., 2) softmax(x3, 2)


text:
module Magrad

export Tensor, backward!

include("tensor.jl")

end

max_stars_repo_name: leonvol/magrad

max_issues_repo_name: leonvol/magrad

max_forks_repo_name: leonvol/magrad


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

@testset "Environment-Bond-Basis" begin

@info("Testset Environment-Bond-Basis")

##

using StaticArrays, Test, Random
using LinearAlgebra
using JuLIP.Testing: print_tf
using JuLIP: alloc_temp
using JuLIP.FIO: read_dict, write_dict, load_dict, save_dict

using ACE
using JuLIP: evaluate!, evaluate, evaluate_d!
using ACE: alloc_B, alloc_dB
using ACE.Bonds: envpairbasis

##

function randR(rin, rout)
   R = @SVector randn(3)
   R /= norm(R) * (rin +  (rout-rin) * rand())
end

function rand_env(Nneigs, rnn)
   R0 = randR(0.8*rnn, 1.2*rnn)
   Renv = [ ((0.66*R0) +  randR(0.6*rnn, 2.0 * rnn)) for _=1:Nneigs ]
   return  R0, Renv
end

"""
This applies several random symmetry operations to the bond configuration
(R0, Renv) :
  - rotate it about the R0 axis
  - reflect along the z axis
  - rotate the entire configuration (R0, Renv) -> (Q R0, Q Renv)
  - random permutatio of Renv
"""
function randiso(R0, Renv)
   r̂ = R0 / norm(R0)
   o = R0 / 2
   t = rand() * 2*pi
   Q =  [ cos(t) sin(t) 0; -sin(t) cos(t) 0; 0 0 1 ]
   v = (r̂ - [0,0,1])
   v /= norm(v)
   H = I - 2 * v * v'
   Rot = H * Q * H
   iso = R -> Rot * (R - o) + o
   isot = R -> Rot' * (R - o) + o
   @assert Rot' * Rot ≈ I
   @assert iso(R0) ≈ R0
   @assert isot.(iso.(Renv)) ≈ Renv
   return randglobalrot(randiso_z(R0, shuffle(iso.(Renv)))...)
end

function randglobalrot(R0, Renv)
   Q = qr(randn(3,3)).Q
   iso = R -> SVector((Q * R)...)
   return iso(R0), iso.(Renv)
end

function randiso_z(R0, Renv)
   σ = rand([0, 1])
   r̂ = R0 / norm(R0)
   iso = R -> SVector( R - σ * 2 * dot(r̂, R-R0/2) * r̂  )
   return R0, shuffle(iso.(Renv))
end

##

@info("Basic setup and evaluation test")
r0 = 1.0
Benv = envpairbasis(:X, 3; rnn = r0, rcut0 = 2.0, degree = 5, wenv = 1
   )
@show length(Benv)
tmp = alloc_temp(Benv)
R0, Renv = rand_env(10, r0)
B = alloc_B(Benv)
B1 = evaluate!(B, tmp, Benv, R0, Renv)
B2 = evaluate(Benv, R0, Renv)
println(@test B1 ≈ B2)

##

@info("Symmetry-invariance test")
for ntest = 1:50
   R0_, Renv_ = randiso(R0, Renv)
   print_tf(@test evaluate(Benv, R0_, Renv_) ≈ B1)
end

##

@info("Generate and test an EnvPairPot")
c = rand(ComplexF64, length(Benv))
V = JuLIP.MLIPs.combine(Benv, c)
tmp = alloc_temp(V)

R0, Renv = rand_env(10, r0)
evaluate!(tmp, V, R0, Renv)

for ntest = 1:50
   R0, Renv = rand_env(10, r0)
   B = evaluate(Benv, R0, Renv)
   print_tf(@test evaluate!(tmp, V, R0, Renv) ≈ real(sum(c .* B)))
end

##

@info("Test (de-)serialisation")
DV = write_dict(V)
V1 = read_dict(DV)
println(@test V1 == V)
tmpname = tempname()
save_dict(tmpname, DV)
V2 = read_dict(load_dict(tmpname))
println(@test V == V2)

##
end


text:
function rectangle_assets(Nx, Nz, q, a, R0, k, n0, T0, B0, dt)

    # define flux surfaces
    psi(x, y) = -y
    rmap_in(x, y) = [x, -y]
    rmap_out(x, y) = [x, 2 - y]
    inner_flux = Barrier(psi, 0, rmap_in)
    outer_flux = Barrier(psi, -1, rmap_out)

    # define field-aligned sheaths
    vert = zeros(2, 5)
    vert[:,1] = Float64[0, -10]
    vert[:,2] = Float64[0, 10]
    vert[:,3] = Float64[1, 10]
    vert[:,4] = Float64[1, -10]
    vert[:,5] = Float64[0, -10]
    sheaths = PolyBarrier(vert)

    # define grid
    bigdeltas = Float64[-0.3, 0.3, -0.3]
    bars = Barrier[inner_flux, outer_flux, sheaths]
    r0 = Float64[-0.5, -0.5]
    r1 = Float64[1.5, 1.5]
    grd = Grid((x,y) -> check_if_inside(x, y, bigdeltas, bars), r0, r1, Nx, Nx, Nz)

    @info "Computed Grid" Nk=grd.Nk

    # make 5 points stencils
    m = 5
    Minv = stencil1d(m)
    MinvT = Matrix(stencil2d(m, m))

    # create derivative matrices
    Dx = derivative_matrix(1, 0, Minv, Minv, grd)
    Dy = derivative_matrix(0, 1, Minv, Minv, grd)
    Dxy = derivative_matrix(1, 1, Minv, Minv, grd)
    Dxx = derivative_matrix(2, 0, Minv, Minv, grd)
    Dyy = derivative_matrix(0, 2, Minv, Minv, grd)
    Dxxx = derivative_matrix(3, 0, Minv, Minv, grd)
    Dyyy = derivative_matrix(0, 3, Minv, Minv, grd)
    Dxxy = derivative_matrix(2, 1, Minv, Minv, grd)
    Dxyy = derivative_matrix(1, 2, Minv, Minv, grd)

    @info "Generated x and y derivatives"

    # create field-line derivative matrices
    bx(x, y) = 1 / sqrt(1 + (R0*q/a)^2)
    by(x, y) = 0.0
    bz(x, y) = R0*q*bx(x, y) / a
    Ds, Dss = fieldline_derivatives(bx, by, bz, 0.001, m, m, MinvT, Nz, grd)

    @info "Generated field-line derivatives"

    # generate boundary operators
    smallerdeltas = Float64[-0.1, 0.1, -0.1]
    qs = Float64[0.01, 0.01, 0.01]
    PENs, REFs, DCHLTs, NMANNs = boundary_operators(m, m, MinvT, smallerdeltas, bars, qs, grd)
    P0, P1, P2, P3 = PENs
    R1, R2, R3 = REFs
    DCHLT1, DCHLT2, DCHLT3 = DCHLTs
    NMANN1, NMANN2, NMANN3 = NMANNs

    @info "Generated boundary operators"

    # compute LAM matrix
    LAM = Diagonal(1 ./ (1 .+ dt*diag(P1 + P2 + P3)))

    # setup diffusion problems with diffusion constant k
    L = I - k*(Dxx + Dyy)
    DIFF_lnn = LinearLeftHandSide(P0*L + NMANN1 + NMANN2 + NMANN3, 2/3)
    DIFF_lnTe = LinearLeftHandSide(P0*L + NMANN1 + NMANN2 + NMANN3, 2/3)
    DIFF_lnTi = LinearLeftHandSide(P0*L + NMANN1 + NMANN2 + NMANN3, 2/3)
    DIFF_u = LinearLeftHandSide(P0*L + NMANN1 + NMANN2 + DCHLT3, 2/3)
    DIFF_w = LinearLeftHandSide(P0*L + DCHLT1 + DCHLT2 + DCHLT3, 2/3)
    DIFF_A = LinearLeftHandSide(P0*L + DCHLT1 + DCHLT2 + DCHLT3, 2/3)

    # compute parameters
    params = dimensionless_parameters(a, R0, n0, T0, B0)

    # setup helmoltz problem
    de2 = params[8]
    Lhelm = I - de2*(Dxx + Dyy)
    HHOLTZ = LinearLeftHandSide(P0*Lhelm + DCHLT1 + DCHLT2 + DCHLT3, 2/3)

    @info "Setup linear problems"

    # compute FLXAVG, from middle flux surface
    flxpoints = hcat([Float64[x, 0.5] for x in LinRange(0, 1, 15)]...)
    rows = SparseMatrixCSC[]
    for i=1:15
        dat, js = interpolation_row(flxpoints[:,i]..., m, m, MinvT, grd)
        is = ones(Int64, length(js))
        row = sparse(is, js, dat, 1, grd._Nx*grd._Ny) * transpose(grd.Proj)
        append!(rows, [row])
    end
    FLXAVG = +(rows...) / 15 # XXX

    # define TRGT
    TRGT = f_to_grid((x, y) -> sign(x - 0.5), grd)

    # define source terms
    Sn = f_to_grid((x,y) -> 1/(100*dt), grd)
    STe = f_to_grid((x,y) -> 1/(100*dt), grd)
    STi = f_to_grid((x,y) -> 1/(100*dt), grd)

    # subcycle 20 times per timestep
    N_subcycle = 20

    # placeholder seed for now
    # seed = Float64[0, 0, 0]

    return Assets(grd, Dx, Dy, Dxy, Dxx, Dyy, Dxxx, Dyyy, Dxxy, Dxyy, Ds, Dss,
                     DIFF_lnn, DIFF_lnTe, DIFF_lnTi, DIFF_u, DIFF_w, DIFF_A,
                     HHOLTZ, P0, P1, P2, P3, R1, R2, R3, LAM, DCHLT1, DCHLT2,
                     DCHLT3, NMANN1, NMANN2, NMANN3, FLXAVG, TRGT, Sn, STe,
                     STi, params, dt, N_subcycle)#, seed)
end


text:
using Query
using DataFrames
using NamedTuples

# We want a DataFrame without Nullable's, so
# we need to construct a DataFrame with columns of Array
columns=[]
push!(columns, ["John", "Sally", "Kirk"])
push!(columns, [23., 42., 59.])
push!(columns, [3,5,2])
df = DataFrame(columns, [:name, :age, :children])

x = @from i in df begin
    @where i.age>30. && i.children > 2
    @select {Name=lowercase(i.name)}
    @collect DataFrame
end

println(x)


text:
# site_counting.jl
# ================
#
# Extension to the site counting framework for biological sequences.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/GeneticVariation.jl/blob/master/LICENSE.md

abstract type Mutation <: Position end

"""
A `Conserved` site describes a site where two aligned nucleotides are definately
conserved. By definately conserved this means that the symbols of the site are
non-ambiguity symbols, and they are the same symbol.
"""
struct Conserved <: Mutation end

"""
A `Mutated` site describes a site where two aligned nucleotides are definately
mutated. By definately mutated this means that the symbols of the site are
non-ambiguity symbols, and they are not the same symbol.
"""
struct Mutated <: Mutation end

"""
A `Transition` site describes a site where two aligned nucleotides are definately
mutated, and the type of mutation is a transition mutation.
In other words, the symbols must not be ambiguity symbols, and they must
be different such that they constitute a transition mutation: i.e. A<->G, or C<->T.
"""
struct Transition <: Mutation end

"""
A `Transversion` site describes a site where two aligned nucleotides are
definately mutated, and the type of mutation is a transversion mutation.
In other words, the symbols must not be ambiguity symbols, and they must
be different such that they constitute a transversion mutation: i.e. A<->C,
A<->T, G<->T, G<->C.
"""
struct Transversion <: Mutation end

# Masking functions
# -----------------

@inline function nibble_mask(::Type{Certain}, x::UInt64)
    return nibble_mask(enumerate_nibbles(x), 0x1111111111111111)
end

@inline function nibble_mask(::Type{Certain}, a::UInt64, b::UInt64)
    return nibble_mask(Certain, a) & nibble_mask(Certain, b)
end

# BioSequences.count_sites_bitpar extension
# -----------------------------------------

for A in (DNAAlphabet, RNAAlphabet)
    @eval begin

        # Counter types
        @inline BioSequences.bp_counter_type(::Type{M}, ::Type{$A{4}}) where M <: Mutation = Tuple{Int, Int}
        @inline BioSequences.bp_start_counter(::Type{M}, ::Type{$A{4}}) where M <: Mutation = Int(0), Int(0)

        # Conserved
        @inline bp_chunk_count(::Type{Conserved}, ::Type{$A{2}}, a::UInt64, b::UInt64) = bp_chunk_count(Match, $A{2}, a, b)
        @inline BioSequences.bp_correct_emptyspace(::Type{Conserved}, ::Type{$A{2}}) = true

        @inline function bp_chunk_count(::Type{Conserved}, ::Type{$A{4}}, a::UInt64, b::UInt64)
            k, c = bp_chunk_count(Mutated, $A{4}, a, b)
            return c - k, c
        end

        # Mutated
        @inline bp_chunk_count(::Type{Mutated}, ::Type{$A{2}}, a::UInt64, b::UInt64) = bp_chunk_count(Mismatch, $A{2}, a, b)

        @inline function bp_chunk_count(::Type{Mutated}, ::Type{$A{4}}, a::UInt64, b::UInt64)
            m = nibble_mask(Certain, a, b)
            d = (a ⊻ b) & m
            return count_nonzero_nibbles(d), count_one_nibbles(m)
        end

        # TODO Finish transitions and transversions
        # Transition
        @inline function bp_chunk_count(::Type{Transition}, ::Type{$A{4}}, a::UInt64, b::UInt64)
            m = nibble_mask(Certain, a, b)
            d = (a ⊻ b) & m
            remcount = count_zero_nibbles(d) # Count of the identical nucleotides.
            tve = count_one_nibbles(nibble_mask(0x9999999999999999, d)) # Count the 1001 transversion edge case.
            d &= (d >> 1)
            d &= 0x7777777777777777
            tvc = count_ones(d)
            tve += (16 - tvc - remcount)
        end

        # Transversion
        # TODO

    end
end

@inline BioSequences.bp_update_counter(acc::Tuple{Int,Int}, up::Tuple{Int,Int}) = acc[1] + up[1], acc[2] + up[2]
@inline BioSequences.diag_val(::Type{Tuple{Int,Int}}) = Int(0), Int(0)


text:
##################################################################################
# This file is part of StateSpaceEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020, Bank of Canada
# All rights reserved.
##################################################################################

@inline ModelBaseEcon.geteqn(i::Integer, m::Model) = geteqn(i, m.sstate)

"""
    SolverData

A data structure used during the solution of the steady state problem.
It contains some current state information and some buffers.

!!! warning
    This is an internal data structure. Do not use directly.

### Why is this necessary?
The solution of the steady state problem is done in steps:
  1. The user sets the initial guess. This can be done with either `initial_sstate!` or `clear_sstate!`
  2. Pre-solve step. In this step we look for equations that have only one unknown and attempt to solve it.
  If successful, the unknown is marked as "solved" and is not an unknown anymore, also the equation is marked as
  "solved" and is not considered anymore. The process repeats for as long as it keeps solving.
  3. Solve step. This is where we solve the non-linear system composed of the remaining equations for the
  remaining variables.

This data structure keeps track of which unknowns and equations are solved and which remain to be solved.
It also holds an indexing map that translates the indexes of the variables and equations we solve for
to their original indexes in the full steady state system.

### Fields
  * `point` - a buffer for the current solution values. Some of these may be presolved, which are kept
  fixed while solving the system, while the rest are updated during solver iterations. The length equals
  the total number of steady state variables.
  * `resid` - a buffer for the current residual vector. The entries corresponding to presolved equations
  would normally be all zeros, while the ones corresponding to "active" equations would be updated during
  solver iterations. Lentgh equals the total number of steady state equations.
  * `solve_var` - a Boolean vector, same size as `point`. Value of `true` means that the unknown is "active", while `false` indicates
  that it has been pre-solved.
  * `solve_eqn` - a Boolean vector, same size as `resid`. Value of `true` means that the equation is "active", while `false` indicates
  that it has been pre-solved.
  * `vars_index` - an Integer vector, same length as `point`. Entries corresponding to pre-solved variables hold zeros.
  Entries for active vars are numbered sequentially from 1 to `nvars`
  * `eqns_index` - an Integer vector, same length as `resid`. Presolved equations have a zero here, while active equations
  are numbered sequentially from 1 to `neqns`.
"""
struct SolverData
    "Buffer holding the current solution."
    point::Vector{Float64}
    resid::Vector{Float64}
    "`true` for active vars and `false` for pre-solved vars."
    solve_var::Vector{Bool}
    "`true` for active equations and `false` for pre-solved equations."
    solve_eqn::Vector{Bool}
    "Sequential indexes of active vars. Pre-solved vars have index 0."
    vars_index::Vector{Int64}
    "Sequential indexes of active equations. Pre-solved equations have index 0."
    eqns_index::Vector{Int64}
    "The model equations to solve"
    alleqns::Vector{SteadyStateEquation}
end

function Base.getproperty(sd::SolverData, pname::Symbol)
    if pname == :nvars
        return sum(getfield(sd, :solve_var))
    elseif pname == :neqns
        return sum(getfield(sd, :solve_eqn))
    else
        return getfield(sd, pname)
    end
end

Base.propertynames(::SolverData) = (:nvars, :neqns, fieldnames(SolverData)...)

"""
    SolverData(model, presolve=Val(false); <options>)

Construct a SolverData instance from all variables and equations in the model,
ignoring anything pre-solved.

### Options
  `verbose::Bool` - if not specified it's taken from the model options.
  `tol::Float64` - desired tolerance when checking the residual of presolved equation.
  `presolve::Bool` - if `false`, any pre-solved information is ignored and
  the solver data is set up to solve all equations for all variables.
"""
function SolverData(model::Model; presolve::Bool=true,
            verbose::Bool=model.options.verbose,
            tol::Float64=model.options.tol
)
    # constructor where we're solving all equations for all variables
    local sstate = model.sstate
    local alleqns = ModelBaseEcon.alleqns(sstate)
    local neqns = length(alleqns)
    local nvars = length(sstate.values)
    sd = SolverData(copy(model.sstate.values), Vector{Float64}(undef, neqns),   # point and resid
                    Vector{Bool}(undef, nvars), Vector{Bool}(undef, neqns),     # solve_var and solve_eqn
                    Vector{Int64}(undef, nvars), Vector{Bool}(undef, neqns),    # vars_index and eqns_index
                    alleqns
    )
    if presolve
        # active vars are the ones with `sstate.mask` equal to `false`
        sd.solve_var .= .! sstate.mask
        # active equations are the ones that have any active vars
        for (i, eqn) in enumerate(sd.alleqns)
            sd.solve_eqn[i] = any(sd.solve_var[eqn.vinds])
        end
        # compute the residual at the initial point.
        global_SS_R!(sd.resid, sd.point, sd.alleqns)
        # collect the indexes of pre-solved equation with non-zero residuals.
        bad_eqn_inds = findall(@. (!sd.solve_eqn) & (abs(sd.resid) > tol))
        if !isempty(bad_eqn_inds)
            # Sort so that largest residuals are at the top
            sort!(bad_eqn_inds, lt=(l, r) -> abs(sd.resid[l]) > abs(sd.resid[r]))
            if verbose
                # print the list of bad equations
                sep = "\n    "
                bad_eqn_str = join(("E$i  res=$(sd.resid[i])  $(geteqn(i, sstate))" for i in bad_eqn_inds), sep)
                @warn "The following presolved equations are not satisfied.$(sep)$(bad_eqn_str)"
            end

            # Mark all bad equations and *all* their variables as active
            sd.solve_eqn[bad_eqn_inds] .= true
            for eqn in sd.alleqns[bad_eqn_inds]
                sd.solve_var[eqn.vinds] .= true
            end

            # Equations that are still pre-solved must have only pre-solved variables
            for eqn in sd.alleqns[.!sd.solve_eqn]
                sd.solve_var[eqn.vinds] .= false
            end

            # All masked variables (where sstate.mask == true) are still inactive as well
            sd.solve_var[sstate.mask] .= false

            # # if ssZeroSlope, make sure any slopes that were marked as active are restored.
            # if model.flags.ssZeroSlope
            #     sd.solve_var[2:2:end] .= false
            # end
        end
    else
        if model.flags.ssZeroSlope
            # ssZeroSlope means that all slopes are 0, so only the levels are active.
            sd.solve_var[1:2:end] .= true
            sd.solve_var[2:2:end] .= false
        else
            sd.solve_var .= true
        end
        # all shocks and exogenous are pre-solved
        for (i,v) in enumerate(model.allvars)
            if isshock(v) || isexog(v)
                sd.solve_var[2i .+ (-1:0)] .= false
            elseif issteady(v)
                # steady have presolved slopes
                sd.solve_var[2i] = false
            end
        end
        # all equations are active regardless.
        sd.solve_eqn .= true
        # compute the residual at the initial point.
        global_SS_R!(sd.resid, sd.point, sd.alleqns)
    end
    sd.vars_index .= 0
    sd.vars_index[sd.solve_var] .= 1:sum(sd.solve_var)
    sd.eqns_index .= 0
    sd.eqns_index[sd.solve_eqn] .= 1:sum(sd.solve_eqn)
    return sd
end
@assert precompile(SolverData, (Model,))

"""
    global_SS_RJ(point, sd::SolverData)

When applied to a solver data, computes the residual and Jacobian of the active set of equations
with respect to the active set of variables.
"""
function global_SS_RJ(point::AbstractVector{Float64}, sd::SolverData)
    sd.point[sd.solve_var] .= point
    R = zeros(sd.neqns)
    J = zeros(sd.neqns, sd.nvars)
    if sd.neqns == 0
        if sd.nvars > 0
            # Show an error message, but no exception is thrown
            @error "System is underdetermined"
        end
        return R, J
    end
    for (i, (solve, ind, eqn)) in enumerate(zip(sd.solve_eqn, sd.eqns_index, sd.alleqns))
        solve || continue
        rr, jj = try
            eqn.eval_RJ(sd.point[eqn.vinds])
        catch
            (NaN64, fill(NaN64, size(eqn.vinds)))
        end
        (isnan(rr) || isinf(rr)) && inadmissible_error(i, eqn, sd.point, rr)
        any(@. isnan(jj) | isinf(jj)) && inadmissible_error(i, eqn, sd.point, jj)
        R[ind] = rr
        # assign jj to the `ind`-th row in J.
        # jj contains the entire gradient. we need only the partials w.r.t. the active variables
        l_active = sd.solve_var[eqn.vinds]  # local mask the active variables
        l_index = sd.vars_index[eqn.vinds]  # local set of indexes (only active are valid, pre-solved are zero)
        J[ind, l_index[l_active]] .= jj[l_active]
    end
    return R, J
end
@assert precompile(global_SS_RJ, (Vector{Float64}, SolverData))

"""
    global_SS_RJ(point, sd::SolverData)

When a solver data is given, we compute the residual of the active equations only.
"""
function global_SS_R!(resid::AbstractVector{Float64}, point::AbstractVector{Float64},  sd::SolverData)
    sd.point[sd.solve_var] .= point
    for (i, (solve, ind, eqn)) in enumerate(zip(sd.solve_eqn, sd.eqns_index, sd.alleqns))
        solve || continue
        rr = try
            eqn.eval_resid(sd.point[eqn.vinds])
        catch
            NaN64
        end
        (isnan(rr) || isinf(rr)) && inadmissible_error(i, eqn, sd.point, rr)
        resid[ind] = rr
    end
    return nothing
end
@assert precompile(global_SS_R!, (Vector{Float64}, Vector{Float64}, SolverData))


text:
# Day 12: Passage Pathing
load(::Day{12}, path) = split.(readlines(path), "-")

# Part 1 - count number of paths through caves
function connected_caves(str, data)
    tuples = filter(t -> (str in t), data)
    return delete!(Set(reduce(vcat, filter.((!=)(str), tuples))), "start")
end

function solve1(::Day{12}, data)
    # Pre-compute look-up table of connected caves for each cave
    caves = unique(vcat(data...))
    connections = Dict([(cave, connected_caves(cave, data)) for cave in caves])

    # We can use recursion as there are no cyclic paths
    function count_leaves(visited::Vector{<:AbstractString})
        last_cave = visited[end]
        last_cave == "end" && (return 1) # is leaf itself

        visited_small_caves = filter(c -> islowercase(c[1]), visited)
        next_caves = filter(!in(visited_small_caves), connections[last_cave])
        length(next_caves) == 0 && (return 0) # is dead end

        return sum(count_leaves(vcat(visited, c)) for c in next_caves)
    end

    # Call recursion on start node
    return count_leaves(["start"])
end

# Part 2 – we can visit _a single_ small cave twice
function solve2(::Day{12}, data)
    caves = unique(vcat(data...))
    connections = Dict([(cave, connected_caves(cave, data)) for cave in caves])

    function count_leaves(visited::Vector{<:AbstractString})
        last_cave = visited[end]
        last_cave == "end" && (return 1) # is leaf itself

        visited_small_caves = filter(c -> islowercase(c[1]), visited)

        # The only change here is that we can visit a single small cave twice
        if any(count(i -> i == c, visited_small_caves) >= 2 for c in visited_small_caves)
            next_caves = filter(!in(visited_small_caves), connections[last_cave])
        else
            next_caves = connections[last_cave]
        end

        length(next_caves) == 0 && (return 0) # is dead end
        return sum(count_leaves(vcat(visited, c)) for c in next_caves)
    end

    # Call recursion on start node
    return count_leaves(["start"])
end

testresult1(::Day{12}) = 10
testresult2(::Day{12}) = 36


text:
# just show body, not standalone
function Plots._show(io::IO, ::MIME"text/html", plt::Plots.Plot{Plots.PlotlyBackend})
    write(io, Plots.html_body(plt))
end

## -----

export plotif, trimplot, signchart,
       plot_polar, plot_polar!,
       plot_parametric,   plot_parametric!,
       vectorfieldplot,   vectorfieldplot!,
       vectorfieldplot3d, vectorfieldplot3d,
       arrow, arrow!

export newton_vis

##
## --------------------------------------------------
##

"""
   trimplot(f, a, b, c=20; kwargs...)

Plot f over [a,b] but break graph if it exceeds c in absolute value.
"""
function trimplot(f, a, b, c=20; color=:black, legend=false, kwargs...)
    F = (a,b) -> begin
        fa, fb = f(a), f(b)
        M = max(fa, fb)
        m = min(fa, fb)
        m < -c && return false
        M > c && return false
        true
    end
    xs = range(a, b, length=251)
    cols = find_colors(F, xs, (color, :transparent, :red))
    Plots.plot(xs, f.(xs), colors=cols, legend=legend, kwargs...)
end

"""
    plotif(f, g, a, b)

Plot f colored depending on g >= 0 or not.
"""
function plotif(f, g, a, b)
    xs = range(a, b, length=251)
    cols = identify_colors(g, xs)
    Plots.plot(xs, f.(xs), color=cols, legend=false)
end

"""
   signchart(f, a, b)

Plot f over a,b with different color when negative.
"""
function signchart(f, a, b)
    p = plotif(f, f, a, b)
    Plots.plot!(p, zero)
    p
end

## --- parametric plots

adapted_grid = CalculusWithJulia.PlotUtils.adapted_grid
function xyzs(ab::ClosedInterval, r)
    a, b = abs = extrema(ab)
    xs = sort(union([adapted_grid(t -> r(t)[i], abs)[1] for i in eachindex(r(a))]...))
    unzip(r.(xs))
end

"""
    plot_parametric(ab, r; kwargs...)
    plot_parametric!(ab, r; kwargs...)
    plot_parametric(u, v, F; kwargs...)
    plot_parametric!(u, v, F; kwargs...)

Make a parametric plot of a space curve or parametrized surface

The intervals to plot over are specifed using `a..b` notation, from IntervalSets
"""
plot_parametric(ab::ClosedInterval, r; kwargs...) = Plots.plot(xyzs(ab, r)...; kwargs...)
plot_parametric!(ab::ClosedInterval, r; kwargs...) = Plots.plot!(xyzs(ab, r)...; kwargs...)

plot_polar(ab::ClosedInterval, r; kwargs...) =
    plot_parametric(ab, t->[r(t)*cos(t), r(t)*sin(t)]; kwargs...)
plot_polar!(ab::ClosedInterval, r; kwargs...) =
    plot_parametric!(ab, t->[r(t)*cos(t), r(t)*sin(t)]; kwargs...)

## ----

function XYZs(u::ClosedInterval, v::ClosedInterval, r; n=51)
    us = range(extrema(u)..., length=n)
    vs = range(extrema(v)..., length=n)
    unzip(r.(us', vs))  # Plots.jl style; flip for Makie
end

# needs plotly(); not gr()
plot_parametric(u::ClosedInterval, v::ClosedInterval, F;  kwargs...) =
    Plots.surface(XYZs(u,v,F)...;  kwargs..., )
plot_parametric!(u::ClosedInterval, v::ClosedInterval, F; kwargs...) =
    Plots.surface!(XYZs(u,v,F)...;  kwargs...)

import Contour

"""
    Visualize `F(x,y,z) = c` by plotting assorted contour lines

This graphic makes slices in the `x`, `y`, and/or `z` direction of the 3-D level surface and plots them accordingly. Which slices (and their colors) are specified through a dictionary.

Examples:
```
F(x,y,z) = x^2 + y^2 + x^2
plot_implicit_surface(F, 20)  # 20 slices in z direction
plot_implicit_surface(F, 20, slices=Dict(:x=>:blue, :y=>:red, :z=>:green), nlevels=6) # all 3 shown

# A heart
a,b = 1,3
F(x,y,z) = (x^2+((1+b)*y)^2+z^2-1)^3-x^2*z^3-a*y^2*z^3
plot_implicit_surface(F, xlims=-2..2,ylims=-1..1,zlims=-1..2)
```

Note: Idea [from](https://stackoverflow.com/questions/4680525/plotting-implicit-equations-in-3d).

Not exported.
"""
function plot_implicit_surface(F, c=0;
                       xlim=(-5,5), ylim=xlim, zlim=xlim,
                       nlevels=25,         # number of levels in a direction
                       slices=Dict(:z => :blue), # Dict(:x => :color, :y=>:color, :z=>:color)
                       kwargs...          # passed to initial `plot` call
                       )

    _linspace(rng, n=150) = range(extrema(rng)[1], stop=extrema(rng)[2], length=n)

    X1, Y1, Z1 = _linspace(xlim), _linspace(ylim), _linspace(zlim)

    p = Plots.plot(;legend=false,kwargs...)

    if :x ∈ keys(slices)
        for x in _linspace(xlim, nlevels)
            local X1 = [F(x,y,z) for y in Y1, z in Z1]
            cnt = Contour.contours(Y1,Z1,X1, [c])
            for line in Contour.lines(Contour.levels(cnt)[1])
                ys, zs = Contour.coordinates(line) # coordinates of this line segment
                Plots.plot!(p, x .+ 0 * ys, ys, zs, color=slices[:x])
          end
        end
    end

    if :y ∈ keys(slices)
        for y in _linspace(ylim, nlevels)
            local Y1 = [F(x,y,z) for x in X1, z in Z1]
            cnt = Contour.contours(Z1,X1,Y1, [c])
            for line in Contour.lines(Contour.levels(cnt)[1])
                xs, zs = Contour.coordinates(line) # coordinates of this line segment
                Plots.plot!(p, xs, y .+ 0 * xs, zs, color=slices[:y])
            end
        end
    end

    if :z ∈ keys(slices)
        for z in _linspace(zlim, nlevels)
            local Z1 = [F(x, y, z) for x in X1, y in Y1]
            cnt = Contour.contours(X1, Y1, Z1, [c])
            for line in Contour.lines(Contour.levels(cnt)[1])
                xs, ys = Contour.coordinates(line) # coordinates of this line segment
                Plots.plot!(p, xs, ys, z .+ 0 * xs, color=slices[:z])
            end
        end
    end

    p
end

## ----

"""
   `arrow!(p, v)`

Add the vector `v` to the plot anchored at `p`.

This would just be a call to `quiver`, but there is no 3-D version of that. As well, the syntax for quiver is a bit awkward for plotting just a single arrow. (Though efficient if plotting many).

```
using Plots
r(t) = [sin(t), cos(t), t]
rp(t) = [cos(t), -sin(t), 1]
plot(unzip(r, 0, 2pi)...)
t0 = 1
arrow!(r(t0), r'(t0))
```
"""
function arrow!(plt, p, v; kwargs...)
  if length(p) == 2
      Plots.quiver!(plt, unzip([p])..., quiver=Tuple(unzip([v])); kwargs...)
  elseif length(p) == 3
    # 3d quiver needs support
    # https://github.com/JuliaPlots/Plots.jl/issues/319#issue-159652535
    # headless arrow instead
      Plots.plot!(plt, unzip(p, p+v)...; kwargs...)
	end
end
arrow!(p,v;kwargs...) = arrow!(Plots.current(), p, v; kwargs...)

"""

    vectorfieldplot(V; xlim=(-5,5), ylim=(-5,5), n=10; kwargs...)

V is a function that takes a point and returns a vector (2D dimensions), such as `V(x) = x[1]^2 + x[2]^2`.

The grid `xlim × ylim` is paritioned into (n+1) × (n+1) points. At each point, `pt`, a vector proportional to `V(pt)` is drawn.

This is written to add to an existing plot.

```
plot()  # make a plot
V(x,y) = [x, y-x]
vectorfield_plot!(p, V)
p
```
"""
function vectorfieldplot!(plt, V; xlim=(-5,5), ylim=(-5,5), n=10, kwargs...)

    dx, dy = (xlim[2]-xlim[1])/n, (ylim[2]-ylim[1])/n
    xs, ys = xlim[1]:dx:xlim[2], ylim[1]:dy:ylim[2]

    ps = [[x,y] for x in xs for y in ys]
    vs = V.(ps)
    λ = 0.9 * min(dx, dy) /maximum(norm.(vs))

    Plots.quiver!(plt, unzip(ps)..., quiver=unzip(λ * vs))

end
vectorfieldplot!(V; kwargs...) = vectorfieldplot!(Plots.current(), V; kwargs...)

##
# visualize newtons method
function newton_vis(f, x0, a=Inf,b=-Inf; steps=5, kwargs...)
    xs = Float64[x0]
    for i in 1:steps
        push!(xs, xs[end] - f(xs[end]) / f'(xs[end]))
    end

    m,M = extrema(xs)
    m = min(m, a)
    M = max(M, b)

    p = Plots.plot(f, m, M; linewidth=3, legend=false, kwargs...)
    Plots.plot!(p, zero)
    for i in 1:steps
        Plots.plot!(p, [xs[i],xs[i],xs[i+1]], [0,f(xs[i]), 0])
        Plots.scatter!(p, xs[i:i],[0])
    end
    Plots.scatter!(p, [xs[steps+1]], [0])
    p
end


text:
using WeaklyHard
"""
Example for generating automata for (sets of) weakly-hard constraints
"""

# Constraints
lambda1 = AnyHitConstraint(12, 15)
lambda2 = RowHitConstraint(15, 47)
lambda3 = RowMissConstraint(3)
lambda4 = AnyMissConstraint(2, 3)
lambda5 = RowHitConstraint(2, 6)

# Generate automata
G1 = build_automaton(lambda1)
G2 = build_automaton(lambda2)
G3 = build_automaton(lambda3)

# set construction
Lambda_star = dominant_set(Set([lambda4, lambda5]))
G = build_automaton(Lambda_star)
minimize_automaton!(G)


text:
experiments = Any[]
experiments_dir = "experiments"
experiments_doc_dir = joinpath("generated", "experiments")

for subdir in readdir(joinpath(@__DIR__, "..", experiments_dir))
    experiments_sub = Any[]
    mkpath(joinpath(@__DIR__, "src", experiments_doc_dir, subdir))
    for experiment in readdir(joinpath(@__DIR__, "..", experiments_dir, subdir))
        endswith(experiment, ".jl") || continue

        experiment_file = joinpath(experiments_dir, subdir, experiment)
        experiment_doc_file = joinpath(
            experiments_doc_dir,
            subdir,
            replace(experiment, ".jl" => ".md"),
        )

        helpstring = cd(joinpath(@__DIR__, "..")) do # run from top-level so we get nice paths
            read(`$(Base.julia_cmd()) --project=$(Base.active_project()) $experiment_file --help`)
        end

        open(joinpath(@__DIR__, "src", experiment_doc_file), write = true) do f
            println(f, "# `$experiment`")
            println(f, "")
            println(f, "```")
            write(f, helpstring)
            println(f, "```")
        end
        push!(experiments_sub, experiment => experiment_doc_file)
    end
    push!(experiments, subdir => experiments_sub)
end


text:
function time_series_intraday_extended(symbol::String, interval::String="60min", slice::String="year1month1"; client = GLOBAL[], parser = "default")
    @argcheck in(interval, ["1min", "5min", "15min", "30min", "60min"])
    sliceMatch = match(r"year(?<year>\d+)month(?<month>\d+)", slice)
    @argcheck !Compat.isnothing(sliceMatch)
    @argcheck parse(Int, sliceMatch["year"]) > 0
    @argcheck parse(Int, sliceMatch["year"]) < 3
    @argcheck parse(Int, sliceMatch["month"]) > 0
    @argcheck parse(Int, sliceMatch["month"]) < 13
    uri = _form_uri_head(client, "TIME_SERIES_INTRADAY_EXTENDED") * "&symbol=$symbol&interval=$interval&slice=$slice" * _form_uri_tail(client, nothing, nothing)
    data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
    p = _parser(parser, "csv")
    return p(data)
end

function time_series_intraday(symbol::String, interval::String="1min"; client = GLOBAL[], outputsize::String="compact", datatype::String="json", parser = "default")
    @argcheck in(interval, ["1min", "5min", "15min", "30min", "60min"])
    @argcheck in(outputsize, ["compact", "full"])
    @argcheck in(datatype, ["json", "csv"])
    uri = _form_uri_head(client, "TIME_SERIES_INTRADAY") * "&symbol=$symbol&interval=$interval" * _form_uri_tail(client, outputsize, datatype)
    data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
    p = _parser(parser, datatype)
    return p(data)
end

for func in (:daily, :daily_adjusted, :weekly, :weekly_adjusted, :monthly, :monthly_adjusted)
    x = "time_series_$(func)"
    fname = Symbol(x)
    @eval begin
        function ($fname)(symbol::String; client = GLOBAL[], outputsize::String="compact", datatype::String="json", parser = "default")
            @argcheck in(outputsize, ["compact", "full"])
            @argcheck in(datatype, ["json", "csv"])
            uri = _form_uri_head(client, uppercase($x)) * "&symbol=$symbol" * _form_uri_tail(client, outputsize, datatype)
            data = retry(_get_request, delays=Base.ExponentialBackOff(n=3, first_delay=5, max_delay=1000))(uri)
            p = _parser(parser, datatype)
            return p(data)
        end

        export $fname
    end
end


text:
# polar plots

#JL Plots
using Plots

r(θ) = 1 + cos(θ) * sin(θ)^2
Plots.plot(r, 0, 2π, proj=:polar, lims=(0,1.5))

r(θ) = 1 + cos(θ) * sin(θ)^2
θ = range(0, stop=2π, length=50)
Plots.plot(θ, r.(θ), proj=:polar, lims=(0,1.5), c=:steelblue,  linewidth=2)

# Polar with markers

Θ = range(0, stop = 2π, length = 100)
pr = abs.(0.1 * randn(100) + sin.(3Θ))
Plots.plot(Θ, pr, proj = :polar, m = 2, c=:steelblue,  linewidth=2)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
export CsvChunkIterator
using CSV
using DataFrames: DataFrame, names

import Base: iterate, length, IteratorSize
using Base.Iterators

"""
    CsvChunkIterator("path/to/file.csv")

Define a Chunking iterator on CSV file
"""
mutable struct CsvChunkIterator
    file::IOStream
    step::Int
    column_headers::Union{Vector{String}, Vector{Symbol}}
    csv_rows_params

    function CsvChunkIterator(path::String, chunk_byte_size = 2^30; csv_rows_params...)
         new(open(path, "r"), chunk_byte_size, Symbol[], csv_rows_params)
    end
end

Base.iterate(chunk_iterator::CsvChunkIterator) = begin
    first_read = position(chunk_iterator.file) == 0
    bytes_read = read(chunk_iterator.file, chunk_iterator.step)

    # try to find the newline character
    last_newline_pos = findlast(x->x==UInt8('\n'), bytes_read)

    # no more to be read
    # if it's nothing then

    if isnothing(last_newline_pos) & (length(bytes_read) == 0)
        close(chunk_iterator.file)
        return nothing
    # have not found a new line, so continue
    elseif isnothing(last_newline_pos)
        if eof(chunk_iterator.file)
            close(chunk_iterator.file)
            # do nothing and go to CSV.read
            last_newline_pos = length(bytes_read)
        else
            # increase step size but doubling
            chunk_iterator.step = 2chunk_iterator.step
            # go back to the beginning
            seek(chunk_iterator.file, position(chunk_iterator.file) - length(bytes_read))
            return iterate(chunk_iterator)
        end
    end

    if first_read
        df =
            CSV.read(
                IOBuffer(
                    @view bytes_read[1:last_newline_pos]
                ), DataFrame;
                chunk_iterator.csv_rows_params...
            )

        chunk_iterator.column_headers = names(df)

        # removes header options from table
        c = chunk_iterator.csv_rows_params
        #d = Dict(zip(keys(c), values(c))...)
        d = Dict(c)
        delete!(d, :head)
        chunk_iterator.csv_rows_params = (;d...)
    else
        df =
            CSV.read(
                IOBuffer(
                    @view bytes_read[1:last_newline_pos]
                ), DataFrame;
                header=chunk_iterator.column_headers,
                chunk_iterator.csv_rows_params...
            )
    end

    new_pos = position(chunk_iterator.file) - (length(bytes_read) - last_newline_pos)
    seek(chunk_iterator.file, new_pos)
    return df, nothing
end

Base.iterate(chunk_iterator::CsvChunkIterator, _) = Base.iterate(chunk_iterator)

# this is needed for `[a for a in chunk_iterator]` to work properly
Base.IteratorSize(chunk_iterator::CsvChunkIterator) = Base.SizeUnknown()

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
### A Pluto.jl notebook ###
# v0.12.6

using Markdown
using InteractiveUtils

# ╔═╡ ce2dfd4a-1deb-11eb-31f5-4d237b2d0610
include("dg_operators.jl");

# ╔═╡ 8e8eb794-1df0-11eb-3f3b-9f2c61ce3c3a
import Printf: @sprintf

# ╔═╡ 8e8f981c-1df0-11eb-2a06-c5686fc02af0
import LinearAlgebra: eigvals

# ╔═╡ d1a69f42-1df0-11eb-09b3-215314f63339
import Plots

# ╔═╡ 8e9477c4-1df0-11eb-188b-2ff136c41a1a
# First we define the computational parameters
begin
  # Floating Point Type
  FT = Float64

  # Polynomial Order
  N = 1

  # number of elements
  num_elems = 11

  # Final time for PDE simulation
  t_final = FT(1 // 10)

  # create the dt data for the simulation
  dg = dg_data(N, num_elems, FT(0), 2FT(π))
end;

# ╔═╡ 8ea3c5d0-1df0-11eb-0d3b-8bfedcbaa3ad
begin
  # First we test the function
  ρ = sin.(dg.x)

  # which has the property that ∂^4 ρ / ∂_x^4 = ρ
  D4ρ = reshape(dg.D4 * ρ[:], size(ρ))
end;

# ╔═╡ e083d63c-1df1-11eb-21c0-f5312717ad4d
begin
  # Plot the solution and the approximation... not so good!
  Plots.plot(dg.x, ρ, color = :red, legend = :none)
  Plots.plot!(dg.x, D4ρ, color = :blue, legend = :none)
end

# ╔═╡ c3740ae6-1df4-11eb-333f-a7ae14f8982e
begin
  # Looking at the error we see the obvious badnees :(
  Plots.plot(dg.x, ρ - D4ρ, color = :red, legend = :none)
end

# ╔═╡ 2ae63aa8-1df4-11eb-00d7-c328b731c322
begin
  # Let's instead evolve the PDE ρ_t = ρ_xxxx some small time
  # with ρ₀(x) = sin(x) the exact solution is ρ(x, t) = exp(-t) sin(x)
  ρ2 = reshape(exp(-t_final * dg.D4) * ρ[:], size(ρ))
  D4ρ2 = reshape(dg.D4 * ρ2[:], size(ρ2))
  D4ρ2_exact = exp(-t_final) * ρ
end;

# ╔═╡ f5a5924e-1df1-11eb-30e4-3908a1ac07cc
begin
  # Plot the solution and the approximation... looking better!
  Plots.plot(dg.x, ρ2, color = :red, legend = :none)
  Plots.plot!(dg.x, D4ρ2, color = :blue, legend = :none)
end

# ╔═╡ d278db02-1df4-11eb-0a90-d90e75af7ea3
begin
  # Looking at the error things are converging
  Plots.plot(dg.x, ρ2 - D4ρ2, color = :red, legend = :none)
end

# ╔═╡ fec8c72a-1e29-11eb-3ce4-af2701a644a7
# Now let's do a refinement study
let
  # These are the levels we will compute
  Ks = 2 .^ (1:8)

  # Loop over refinement levels and compute error
  err = zeros(FT, length(Ks))
  for (iter, K) in enumerate(Ks)
    dg = dg_data(N, K, FT(0), 2FT(π))
    ρ = sin.(dg.x)
    ρ2 = reshape(exp(-t_final * dg.D4) * ρ[:], size(ρ))
    err[iter] = sqrt(sum(dg.M .* (ρ2 - ρ * exp(-t_final)).^2))
  end

  # Compute the rates
  rate = (log2.(err[1:end-1]) - log2.(err[2:end])) ./
  (log2.(Ks[2:end])    - log2.(Ks[1:end-1]))

  # Make a nice ouput string
  str = @sprintf "```polynomial order = %d```\n" N
  str *= @sprintf "
  ```num elements = %03d :: rate (error) = XXXX (%.3e)```\n
  " Ks[1] err[1]
  for (K, e, r) in zip(Ks[2:end], err[2:end], rate)
    str *= @sprintf "```num elements = %03d :: rate (error) = %7.2f (%.3e)```\n\n" K r e
  end
  Markdown.parse(str)
end

# ╔═╡ Cell order:
# ╠═ce2dfd4a-1deb-11eb-31f5-4d237b2d0610
# ╠═8e8eb794-1df0-11eb-3f3b-9f2c61ce3c3a
# ╠═8e8f981c-1df0-11eb-2a06-c5686fc02af0
# ╠═d1a69f42-1df0-11eb-09b3-215314f63339
# ╠═8e9477c4-1df0-11eb-188b-2ff136c41a1a
# ╠═8ea3c5d0-1df0-11eb-0d3b-8bfedcbaa3ad
# ╠═e083d63c-1df1-11eb-21c0-f5312717ad4d
# ╠═c3740ae6-1df4-11eb-333f-a7ae14f8982e
# ╠═2ae63aa8-1df4-11eb-00d7-c328b731c322
# ╠═f5a5924e-1df1-11eb-30e4-3908a1ac07cc
# ╠═d278db02-1df4-11eb-0a90-d90e75af7ea3
# ╠═fec8c72a-1e29-11eb-3ce4-af2701a644a7


text:
rwalk = function(N, M)
    result = zeros(Int, N)
    for i in 1:N
        steps = 2
        rand() >= 0.5 ?  a = 1 : a = -1
        rand() >= 0.5 ? a += 1 : a -= 1
        while a != 0
            steps += 2
            rand() >= 0.5 ? a += 1 : a -= 1
            rand() >= 0.5 ? a += 1 : a -= 1
            if steps >= M
                break
            end
        end
        result[i] = steps
    end
    return result
end

max_stars_repo_name: hwborchers/Presentations

max_issues_repo_name: hwborchers/Presentations

max_forks_repo_name: hwborchers/Presentations