text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    DecisionVariable <: AbstractVariableRef

Identifier for a decision in a stochastic program. Holds a reference to the stochastic program, the stage the decision is taken in, and its corresponding MOI.VariableIndex.
"""
struct DecisionVariable <: JuMP.AbstractVariableRef
    stochasticprogram::StochasticProgram
    stage::Int
    index::MOI.VariableIndex
end

# Getters #
# ========================== #
"""
    decision(dvar::DecisionVariable)

Return the internal `Decision` associated with the first-stage `dvar`.
"""
function decision(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `decision(dvar, scenario_index)`.")
    # Dispatch to structure
    return decision_dispatch(decision,
                             structure(owner_model(dvar)),
                             index(dvar),
                             stage(dvar))
end
"""
    decision(dvar::DecisionVariable, scenario_index::Integer)

Return the scenario-dependent internal `Decision` associated with `dvar` at `scenario_index`.
"""
function decision(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `decision(dvar)`.")
    # Dispatch to structure
    return scenario_decision_dispatch(decision,
                                      structure(owner_model(dvar)),
                                      index(dvar),
                                      stage(dvar),
                                      scenario_index)
end
"""
    stage(dvar::DecisionVariable)

Return the stage of `dvar`.
"""
function stage(dvar::DecisionVariable)
    return dvar.stage
end
"""
    state(dvar::DecisionVariable)

Return the `DecisionState` of the first-stage `dvar`.
"""
function state(dvar::DecisionVariable)
    return decision(dvar).state
end
"""
    state(dvar::DecisionVariable, scenario_index::Integer)

Return the scenario-dependent `DecisionState` of `dvar` at `scenario_index`.
"""
function state(dvar::DecisionVariable, scenario_index::Integer)
    return decision(dvar,scenario_index).state
end

# Setters #
# ========================== #
function take_decisions!(stochasticprogram::StochasticProgram, dvars::Vector{DecisionVariable}, vals::AbstractVector)
    if isempty(dvars)
        @warn "No decision variables specified. Nothing to do."
        return nothing
    end
    # Check that all given decisions are in model
    map(dvar -> check_belongs_to_model(dvar, stochasticprogram), dvars)
    # Check that all decisions belong to the same stage
    all(dvar -> stage(dvar) == stage(dvars[1]), dvars) || error("All decisions should be taken in the same stage.")
    # Check that decisions are first-stage
    stage(dvars[1]) == 1 || error("Decisions are scenario-dependent, consider `take_decisions!(sp, dvars, vals, scenario_index)`")
    # Check decision length
    length(dvars) == length(vals) || error("Given decision of length $(length(vals)) not compatible with number of decision variables $(length(dvars)).")
    # Update decisions
    for (dvar, val) in zip(dvars, vals)
        fix(dvar, val)
    end
    return nothing
end
function take_decisions!(stochasticprogram::StochasticProgram, dvars::Vector{DecisionVariable}, vals::AbstractVector, scenario_index::Integer)
    if isempty(dvars)
        @warn "No decision variables specified. Nothing to do."
        return nothing
    end
    # Check that all given decisions are in model
    map(dvar -> check_belongs_to_model(dvar, stochasticprogram), dvars)
    # Check that all decisions belong to the same stage
    all(dvar -> stage(dvar) == stage(dvars[1]), dvars) || error("All decisions should be taken in the same stage.")
    # Check that decisions are not first-stage
    stage(dvars[1]) > 1 || error("Decisions are not scenario-dependent, consider `take_decisions!(sp, dvars, vals)`")
    # Check decision length
    length(dvars) == length(vals) || error("Given decision of length $(length(vals)) not compatible with number of decision variables $(length(dvars)).")
    # Update decisions
    for (dvar, val) in zip(dvars, vals)
        fix(dvar, scenario_index, val)
    end
    return nothing
end

function untake_decisions!(stochasticprogram::StochasticProgram, dvars::Vector{DecisionVariable})
    if isempty(dvars)
        @warn "No decision variables specified. Nothing to do."
        return nothing
    end
    # Check that all given decisions are in model
    map(dvar -> check_belongs_to_model(dvar, model), dvars)
    # Check that all decisions belong to the same stage
    all(dvar -> stage(dvar) == stage(dvars[1]), dvars) || error("All decisions should be taken in the same stage.")
    # Check that decisions are first-stage
    stage(dvars[1]) == 1 || error("Decisions are scenario-dependent, consider `unttake_decisions!(sp, dvars, scenario_index)`")
    # Update decisions
    need_update = false
    for dvar in dvars
        unfix(dvar)
    end
    return nothing
end
function untake_decisions!(stochasticprogram::StochasticProgram, dvars::Vector{DecisionVariable}, scenario_index::Integer)
    if isempty(dvars)
        @warn "No decision variables specified. Nothing to do."
        return nothing
    end
    # Check that all given decisions are in model
    map(dvar -> check_belongs_to_model(dvar, model), dvars)
    # Check that all decisions belong to the same stage
    all(dvar -> stage(dvar) == stage(dvars[1]), dvars) || error("All decisions should be taken in the same stage.")
    # Check that decisions are first-stage
    stage(dvars[1]) > 1 || error("Decisions are not scenario-dependent, consider `take_decisions!(sp, dvars, vals)`")
    # Update decisions
    for dvar in dvars
        unfix(dvar, scenario_index)
    end
    return nothing
end

# MOI #
# ========================== #
function MOI.get(stochasticprogram::StochasticProgram, attr::MOI.AbstractVariableAttribute,
                 dvar::DecisionVariable)
    check_belongs_to_model(dvar, stochasticprogram)
    if MOI.is_set_by_optimize(attr)
        # Check if there is a cached solution
        cache = solutioncache(stochasticprogram)
        if haskey(cache, :solution)
            # Returned cached solution if possible
            try
                return MOI.get(cache[:solution], attr, index(dvar))
            catch
            end
        end
        if haskey(cache, :node_solution_1)
            # Value was possibly only cached in first-stage solution
            try
                return MOI.get(cache[:node_solution_1], attr, index(dvar))
            catch
            end
        end
        check_provided_optimizer(stochasticprogram.optimizer)
        if MOI.get(stochasticprogram, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
            throw(OptimizeNotCalled())
        end
        return MOI.get(optimizer(stochasticprogram), attr, index(dvar))
    end
    return MOI.get(backend(proxy(stochasticprogram, stage(dvar))), attr, index(dvar))
end
function MOI.get(stochasticprogram::StochasticProgram, attr::ScenarioDependentVariableAttribute,
                 dvar::DecisionVariable)
    check_belongs_to_model(dvar, stochasticprogram)
    if MOI.is_set_by_optimize(attr)
        # Check if there is a cached solution
        cache = solutioncache(stochasticprogram)
        key = Symbol(:node_solution_, attr.stage, :_, attr.scenario_index)
        if haskey(cache, key)
            # Returned cached solution if possible
            try
                return MOI.get(cache[key], attr.attr, index(dvar))
            catch
            end
        end
        check_provided_optimizer(stochasticprogram.optimizer)
        if MOI.get(stochasticprogram, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
            throw(OptimizeNotCalled())
        end
        try
            # Try to get scenario-dependent value directly
            return MOI.get(optimizer(stochasticprogram), attr, index(dvar))
        catch
            # Fallback to resolving scenario-dependence in structure if
            # not supported natively by optimizer
            return scenario_decision_dispatch(
                structure(stochasticprogram),
                index(dvar),
                stage(dvar),
                attr.scenario_index,
                attr.attr,
            ) do dref, attr
                return MOI.get(owner_model(dref), attr, dref)
            end
        end
    end
    # Get value from structure if not set by optimizer
    return decision_dispatch(
        structure(stochasticprogram),
        index(dvar),
        stage(dvar),
        attr.scenario_index,
        attr.attr,
    ) do dref, attr
        return MOI.get(owner_model(dref), attr, dref)
    end
end

function MOI.set(stochasticprogram::StochasticProgram, attr::MOI.AbstractVariableAttribute,
                 dvar::DecisionVariable, value)
    check_belongs_to_model(dvar, stochasticprogram)
    # Dispatch to structure
    return decision_dispatch!(
        structure(owner_model(dvar)),
        index(dvar),
        stage(dvar),
        attr,
        value,
    ) do dref, attr, value
        MOI.set(owner_model(dref), attr, dref, value)
    end
    return nothing
end
function MOI.set(stochasticprogram::StochasticProgram, attr::ScenarioDependentVariableAttribute,
                 dvar::DecisionVariable, value)
    check_belongs_to_model(dvar, stochasticprogram)
    # Dispatch to structure
    return scenario_decision_dispatch!(
        structure(owner_model(dvar)),
        index(dvar),
        stage(dvar),
        scenario_index,
        attr,
        value,
    ) do dref, attr, value
        MOI.set(owner_model(dref), attr, dref, value)
    end
end

# JuMP variable interface #
# ========================== #
"""
    name(dvar::DecisionVariable, scenario_index::Integer)::String

Get the name of the decision variable `dvar`.
"""
function JuMP.name(dvar::DecisionVariable)::String
    return MOI.get(backend(proxy(owner_model(dvar), stage(dvar))), MOI.VariableName(), index(dvar))::String
end
"""
    name(dvar::DecisionVariable, scenario_index::Integer)::String

Get the name of the scenario-dependent decision variable `dvar` in scenario `scenario_index`.
"""
function JuMP.name(dvar::DecisionVariable, scenario_index::Integer)::String
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `name(dvar)`.")
    # Dispatch to structure
    scenario_decision_dispatch(JuMP.name,
                               structure(owner_model(dvar)),
                               index(dvar),
                               stage(dvar),
                               scenario_index)::String
end
"""
    set_name(dvar::DecisionVariable, scenario_index::Integer, name::String)

Set the name of the decision variable `dvar` to `name`.
"""
function JuMP.set_name(dvar::DecisionVariable, name::String)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `set_name(dvar, scenario_index, name)`.")
    # Dispatch to structure
    decision_dispatch!(JuMP.set_name,
                       structure(owner_model(dvar)),
                       index(dvar),
                       stage(dvar),
                       name)
    return nothing
end
"""
    set_name(dvar::DecisionVariable, scenario_index::Integer, name::String)

Set the name of the scenario-dependent decision variable `dvar` in scenario `scenario_index` to `name`.
"""
function JuMP.set_name(dvar::DecisionVariable, scenario_index::Integer, name::String)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `set_name(dvar, name)`.")
    # Dispatch to structure
    scenario_decision_dispatch!(JuMP.set_name,
                                structure(owner_model(dvar)),
                                index(dvar),
                                stage(dvar),
                                scenario_index,
                                name)
    return nothing
end
"""
    decision_by_name(stochasticprogram::Stochasticprogram,
                     stage::Integer,
                     name::String)::Union{AbstractVariableRef, Nothing}

Returns the reference of the variable with name attribute `name` at `stage` of `stochasticprogram` or `Nothing` if
no variable has this name attribute. Throws an error if several variables have
`name` as their name attribute at stage `s`.
"""
function decision_by_name(stochasticprogram::StochasticProgram{N}, stage::Integer, name::String) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    index = MOI.get(backend(proxy(stochasticprogram, stage)), MOI.VariableIndex, name)
    if index isa Nothing
        return nothing
    else
        return DecisionVariable(stochasticprogram, stage, index)
    end
end
"""
    index(dvar::DecisionVariable)::MOI.VariableIndex

Return the index of the decision variable that corresponds to `dvar` in the MOI backend.
"""
function JuMP.index(dvar::DecisionVariable)
    return dvar.index
end
"""
    optimizer_index(dvar::DecisionVariable)::MOI.VariableIndex

Return the index of the variable that corresponds to `dvar` in the optimizer model.
"""
function JuMP.optimizer_index(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `optimizer_index(dvar, scenario_index)`.")
    return JuMP._moi_optimizer_index(structure(owner_model(dvar)), index(dvar))
end
"""
    optimizer_index(dvar::DecisionVariable, scenario_index)::MOI.VariableIndex

Return the index of the variable that corresponds to the scenario-dependent `dvar` in the optimizer model at `scenario_index`.
"""
function JuMP.optimizer_index(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `optimizer_index(dvar)`.")
    return JuMP._moi_optimizer_index(structure(owner_model(dvar)), index(dvar), scenario_index)
end
"""
    has_duals(stochasticprogram::StochasticProgram; result::Int = 1)

Return `true` if the solver has a primal solution in the first-stage of
`stochasticprogram` in result index `result` available to query,
otherwise return `false`.
"""
function JuMP.has_values(stochasticprogram::StochasticProgram; result::Int = 1)
    return primal_status(stochasticprogram; result = result) != MOI.NO_SOLUTION
end
"""
    has_values(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer; result::Int = 1)

Return `true` if the solver has a primal solution in the node at stage
`stage` and scenario `scenario_index` in result index `result`
available to query, otherwise return `false`.
"""
function JuMP.has_values(stochasticprogram::StochasticProgram, stage::Integer, scenario_index::Integer; result::Int = 1)
    return primal_status(stochasticprogram, stage, scenario_index; result = result) != MOI.NO_SOLUTION
end
"""
    has_values(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer; result::Int = 1)

Return `true` if the solver has a primal solution in scenario `scenario_index`
in result index `result` available to query, otherwise return `false`.
"""
function JuMP.has_values(stochasticprogram::TwoStageStochasticProgram, scenario_index::Integer; result::Int = 1)
    return has_values(stochasticprogram, 2, scenario_index; result = result)
end
"""
    value(dvar::DecisionVariable; result = 1)

Return the value of the first-stage decision variable `dvar`
associated with result index `result` of the most-recent
returned by the solver.
"""
function JuMP.value(dvar::DecisionVariable; result::Int = 1)::Float64
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `value(dvar, scenario_index)`.")
    d = decision(dvar)
    if d.state == Taken
        # If decision has been fixed the value can be fetched
        # directly
        return d.value
    end
    return MOI.get(owner_model(dvar), MOI.VariablePrimal(result), dvar)
end
"""
    value(dvar::DecisionVariable, scenario_index::Integer; result = 1)

Return the scenario-dependent value of the decision variable `dvar`
associated with result index `result` at `scenario_index` of the
most-recent returned by the solver.
"""
function JuMP.value(dvar::DecisionVariable, scenario_index::Integer; result::Int = 1)::Float64
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `value(dvar)`.")
    d = decision(dvar, scenario_index)
    if d.state == Taken
        # If decision has been fixed the value can be fetched
        # directly
        return d.value
    end
    attr = ScenarioDependentVariableAttribute(stage(dvar), scenario_index, MOI.VariablePrimal(result))
    return MOI.get(owner_model(dvar), attr, dvar)
end
"""
    value(dvar_expr::Union{GenericAffExpr{T,DecisionVariable}, GenericQuadExpr{T,DecisionVariable}}, stage_to_scenario::Dict{Int,Int}) where T

Evaluate `dvar_expr` where the value of a given `dvar` is found in the scenario returned by the provided `stage_to_scenario` map.
"""
function JuMP.value(dvar_expr::Union{GenericAffExpr{T,DecisionVariable}, GenericQuadExpr{T,DecisionVariable}}, stage_to_scenario::Dict{Int,Int}; result::Int = 1)::Float64 where T
    var_value = (dvar) -> begin
        if stage(dvar) == 1
            return value(dvar; result = result)
        end
        return value(dvar, stage_to_scenario[stage(dvar)]; result = result)
    end
    return value(dvar_expr, var_value)
end
"""
    reduced_cost(dvar::DecisionVariable)::Float64

Return the reduced cost associated with the decision variable `dvar`.
"""
function JuMP.reduced_cost(dvar::DecisionVariable)::Float64
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `reduced_cost(dvar, scenario_index)`.")
    sp = owner_model(dvar)
    if !has_duals(sp)
        error("Unable to query reduced cost of variable because stochastic program does" *
              " not have duals available in the first stage.")
    end
    sign = objective_sense(sp) == MOI.MIN_SENSE ? 1.0 : -1.0
    if is_fixed(dvar)
        return sign * dual(FixRef(dvar))
    end
    rc = 0.0
    if has_upper_bound(dvar)
        rc += dual(UpperBoundRef(dvar))
    end
    if has_lower_bound(dvar)
        rc += dual(LowerBoundRef(dvar))
    end
    return sign * rc
end
"""
    reduced_cost(dvar::DecisionVariable)::Float64

Return the reduced cost associated with the scenario-dependent decision variable `dvar` at `scenario_index`.
"""
function JuMP.reduced_cost(dvar::DecisionVariable, scenario_index::Integer)::Float64
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `reduced_cost(dvar)`.")
    sp = owner_model(dvar)
    if !has_duals(sp, scenario_index)
        error("Unable to query reduced cost of variable because stochastic program does" *
              " not have duals available in scenario `scenario_index`.")
    end
    sign = objective_sense(sp) == MOI.MIN_SENSE ? 1.0 : -1.0
    if is_fixed(dvar, scenario_index)
        return sign * dual(FixRef(dvar), scenario_index)
    end
    rc = 0.0
    if has_upper_bound(dvar, scenario_index)
        rc += dual(UpperBoundRef(dvar), scenario_index)
    end
    if has_lower_bound(dvar, scenario_index)
        rc += dual(LowerBoundRef(dvar), scenario_index)
    end
    return sign * rc
end
"""
    is_fixed(dvar::DecisionVariable)

Return `true` if `dvar` is a fixed first-stage decision variable.
"""
function JuMP.is_fixed(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `is_fixed(dvar, scenario_index)`.")
    if state(dvar) == Taken
        return true
    end
    return false
end
"""
    is_fixed(dvar::DecisionVariable, scenario_index)

Return `true` if `dvar` is a fixed decision variable in `scenario_index`.
"""
function JuMP.is_fixed(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `is_fixed(dvar)`.")
    if state(dvar, scenario_index) == Taken
        return true
    end
    return false
end
"""
    fix(dvar::DecisionVariable, val::Number)

Fix the first-stage decision associated with `dvar` to `val`. In contexts
where `dvar` is a variable, the variable is fixed to the value. In
contexts where `dvar` is a known parameter value, the value is updated.
"""
function JuMP.fix(dvar::DecisionVariable, val::Number)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `fix(dvar, scenario_index, val)`.")
    fix(structure(owner_model(dvar)), index(dvar), stage(dvar), val)
    return nothing
end
"""
    fix(dvar::DecisionVariable, scenario_index::Integer, val::Number)

Fix the scenario-dependent decision associated with `dvar` at `scenario_index`
to `val`. In contexts where `dvar` is a variable, the variable is fixed to the
value. In contexts where `dvar` is a known parameter value, the value is updated.
"""
function JuMP.fix(dvar::DecisionVariable, scenario_index::Integer, val::Number)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `fix(dvar, val)`.")
    fix(structure(owner_model(dvar)), index(dvar), stage(dvar), scenario_index, val)
    return nothing
end
"""
    unfix(dvar::DecisionVariable)

Unfix the first-stage decision associated with `dvar`. If the decision
is already in a `NotTaken` state, this does nothing.
"""
function JuMP.unfix(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `unfix(dvar, scenario_index)`.")
    if state(dvar) == NotTaken
        # Nothing to do, just return
        return nothing
    end
    unfix(structure(owner_model(dvar)), index(dvar), stage(dvar))
    return nothing
end
"""
    unfix(dvar::DecisionVariable, scenario_index::Integer)

Unfix the scenario-dependent decision associated with `dvar`
at `scenario_index`. If the decision is already in a `NotTaken`
state, this does nothing.
"""
function JuMP.unfix(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `unfix(dvar)`.")
    if state(dvar, scenario_index) == NotTaken
        # Nothing to do, just return
        return nothing
    end
    unfix(structure(owner_model(dvar)), index(dvar), stage(dvar), scenario_index)
    return nothing
end

"""
    FixRef(dvar::DecisionVariable)

Return a constraint reference to the constraint fixing the value of the decision `dvar`. Errors
if one does not exist.
"""
function JuMP.FixRef(dvar::DecisionVariable)
    moi_fix =  MOI.ConstraintIndex{SingleDecision, MOI.EqualTo{Float64}}
    return SPConstraintRef{moi_fix, ScalarShape}(owner_model(dvar),
                                                 stage(dvar),
                                                 moi_fix(index(dvar).value),
                                                 ScalarShape())
end

function JuMP.fix_value(dvar::DecisionVariable)
    cset = MOI.get(owner_model(v), MOI.ConstraintSet(), FixRef(dvar))::MOI.EqualTo{Float64}
    return cset.value
end

function JuMP.owner_model(dvar::DecisionVariable)
    return dvar.stochasticprogram
end

struct DecisionVariableNotOwned <: Exception
    dvar::DecisionVariable
end

function JuMP.check_belongs_to_model(dvar::DecisionVariable, stochasticprogram::StochasticProgram)
    if owner_model(dvar) !== stochasticprogram
        throw(DecisionNotOwned(dvar))
    end
end
"""
    is_valid(stochasticprogram::StochasticProgram, dvar::DecisionVariable)

Return `true` if `dvar` refers to a valid first-stage decision variable in `stochasticprogram`.
"""
function JuMP.is_valid(stochasticprogram::StochasticProgram, dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `is_valid(dvar, scenario_index)`.")
    return stochasticprogram === owner_model(dvar) &&
        MOI.is_valid(structure(stochasticprogram), index(dvar), stage(dvar))
end
"""
    is_valid(stochasticprogram::StochasticProgram, dvar::DecisionVariable, scenario_index::Integer)

Return `true` if the scenario-dependent `dvar` refers to a valid decision variable in `stochasticprogram` at `scenario_index`.
"""
function JuMP.is_valid(stochasticprogram::StochasticProgram, dvar::DecisionVariable, scenario_index::Integer)::Bool
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `is_valid(dvar)`.")
    return stochasticprogram === owner_model(dvar) &&
        MOI.is_valid(structure(stochasticprogram), index(dvar), stage(dvar), scenario_index)
end
"""
    delete(stochasticprogram::StochasticProgram, dvar::DecisionVariable)

Delete the first-stage decision variable associated with `dvar` from the `stochasticprogram`.
"""
function JuMP.delete(stochasticprogram::StochasticProgram, dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `delete(stochasticprogram, dvar, scenario_index)`.")
    if stochasticprogram !== owner_model(dvar)
        error("The decision variable you are trying to delete does not " *
              "belong to the stochastic program.")
    end
    # Dispatch to structure
    decision_dispatch!(structure(owner_model(dvar)),
                       index(dvar),
                       stage(dvar)) do dref
                           JuMP.delete(owner_model(dref), dref)
                           return nothing
                       end
    return nothing
end
"""
    delete(stochasticprogram::StochasticProgram, dvar::DecisionVariable, scenario_index::Integer)

Delete the scenario-dependent decision variable associated with `dvar` from the `stochasticprogram` at `scenario_index`.
"""
function JuMP.delete(stochasticprogram::StochasticProgram, dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `delete(stochasticprogram, dvar)`.")
    if stochasticprogram !== owner_model(dvar)
        error("The decision variable you are trying to delete does not " *
              "belong to the stochastic program.")
    end
    # Dispatch to structure
    scenario_decision_dispatch!(structure(owner_model(dvar)),
                                index(dvar),
                                stage(dvar),
                                scenario_index) do dref
                                    JuMP.delete(owner_model(dref), dref)
                                end
    return nothing
end
"""
    delete(stochasticprogram::StochasticProgram, dvars::Vector{DecisionVariable})

Delete the decisions associated with `dvars` from the `stochasticprogram`.
"""
function JuMP.delete(stochasticprogram::StochasticProgram, dvars::Vector{DecisionVariable})
    all(stage.(dvars) .== stage(dvars[1])) || error("$dvars are not all from the same stage")
    stage(dvars[1]) > 1 && error("Some of $dvars are scenario dependent, consider `delete(stochasticprogram, dvars, scenario_index)`.")
    if any(stochasticprogram !== owner_model(dvar) for dvar in dvars)
        error("A decision variable you are trying to delete does not " *
              "belong to the stochastic program.")
    end
    MOI.delete(structure(stochasticprogram), index.(dvars), stage(dvars[1]))
    return nothing
end
"""
    delete(stochasticprogram::StochasticProgram, dvars::Vector{DecisionVariable}, scenario_index::Integer)

Delete the scenario-dependent decisions associated with `dvars` from the `stochasticprogram` at `scenario_index`.
"""
function JuMP.delete(stochasticprogram::StochasticProgram, dvars::Vector{DecisionVariable}, scenario_index::Integer)
    all(stage.(dvars) .== stage(dvars[1])) || error("$dvars are not all from the same stage")
    stage(dvars[1]) == 1 && error("$dvars are not scenario dependent, consider `delete(stochasticprogram, dvars)`.")
    if any(stochasticprogram !== owner_model(dvar) for dvar in dvars)
        error("A decision variable you are trying to delete does not " *
              "belong to the stochastic program.")
    end
    MOI.delete(structure(stochasticprogram), index.(dvars), stage(dvars[1]), scenario_index)
    return nothing
end
"""
    has_lower_bound(dvar::DecisionVariable)

Return `true` if the first-stage decision variable `dvar` has a lower bound.
"""
function JuMP.has_lower_bound(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `has_lower_bound(dvar, scenario_index)`.")
    # Dispatch to structure
    return decision_dispatch(JuMP.has_lower_bound,
                             structure(owner_model(dvar)),
                             index(dvar),
                             stage(dvar))::Bool
end
"""
    has_lower_bound(dvar::DecisionVariable, scenario_index::Integer)

Return `true` if the scenario-dependent decision variable `dvar` has a lower bound at `scenario_index`.
"""
function JuMP.has_lower_bound(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `has_lower_bound(dvar)`.")
    # Dispatch to structure
    return scenario_decision_dispatch(JuMP.has_lower_bound,
                                      structure(owner_model(dvar)),
                                      index(dvar),
                                      stage(dvar),
                                      scenario_index)::Bool
end
"""
    LowerBoundRef(dvar::DecisionVariable)

Return a constraint reference to the lower bound constraint of the decision variable `dvar`.
Errors if one does not exist.
"""
function JuMP.LowerBoundRef(dvar::DecisionVariable)
    moi_lb =  MOI.ConstraintIndex{SingleDecision, MOI.GreaterThan{Float64}}
    return SPConstraintRef{moi_lb, ScalarShape}(owner_model(dvar),
                                                stage(dvar),
                                                moi_lb(index(dvar).value),
                                                ScalarShape())
end
"""
    set_lower_bound(dvar::DecisionVariable)

Set the lower bound of the first-stage decision variable `dvar` to `lower`. If one does not exist, create a new lower bound constraint.
"""
function JuMP.set_lower_bound(dvar::DecisionVariable, lower::Number)
    new_set = MOI.GreaterThan(convert(Float64, lower))
    # Update proxy
    proxy_ = proxy(owner_model(dvar), stage(dvar))
    dref = DecisionRef(proxy_, index(dvar))
    set_lower_bound(dref, lower)
    # Update structure
    if stage(dvar) == 1
        # Dispatch to structure
        decision_dispatch!(JuMP.set_lower_bound,
                           structure(owner_model(dvar)),
                           index(dvar),
                           stage(dvar),
                           lower)
    else
        # Update all constraints at stage
        for scenario_index in 1:num_scenarios(owner_model(dvar), stage(dvar))
            set_lower_bound(dvar, scenario_index, lower)
        end
    end
    return nothing
end
"""
    set_lower_bound(dvar::DecisionVariable, scenario_index::Integer, lower::Number)

Set the lower bound of the scenario-dependent decision variable `dvar` at `scenario_index` to `lower`. If one does not exist, create a new lower bound constraint.
"""
function JuMP.set_lower_bound(dvar::DecisionVariable, scenario_index::Integer, lower::Number)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `set_lower_bound(dvar, lower)`.")
    new_set = MOI.GreaterThan(convert(Float64, lower))
    # Dispatch to structure
    scenario_decision_dispatch!(JuMP.set_lower_bound,
                                structure(owner_model(dvar)),
                                index(dvar),
                                stage(dvar),
                                scenario_index,
                                lower)
    return nothing
end
"""
    delete_lower_bound(dvar::DecisionVariable)

Delete the lower bound constraint of the first-stage decision variable `dvar`.
"""
function JuMP.delete_lower_bound(dvar::DecisionVariable)
    # Update proxy
    proxy_ = proxy(owner_model(dvar), stage(dvar))
    dref = DecisionRef(proxy_, index(dvar))
    delete_lower_bound(dref)
    # Update structure
    if stage(dvar) == 1
        # Dispatch to structure
        decision_dispatch!(JuMP.delete_lower_bound,
                           structure(owner_model(dvar)),
                           index(dvar),
                           stage(dvar))
    else
        # Update all constraints at stage
        for scenario_index in 1:num_scenarios(owner_model(dvar), stage(dvar))
            delete_lower_bound(dvar, scenario_index, lower)
        end
    end
    return nothing
end
"""
    delete_lower_bound(dvar::DecisionVariable, scenario_index::Integer)

Delete the lower bound constraint of the scenario-dependent decision variable `dvar` at `scenario_index`.
"""
function JuMP.delete_lower_bound(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `delete_lower_bound(dvar)`.")
    # Dispatch to structure
    scenario_decision_dispatch!(JuMP.delete_lower_bound,
                                structure(owner_model(dvar)),
                                index(dvar),
                                stage(dvar),
                                scenario_index)
    return nothing
end
"""
    lower_bound(dvar::DecisionVariable)

Return the lower bound of the first-stage decision variable `dvar`. Error if one does not exist.
"""
function JuMP.lower_bound(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `lower_bound(dvar, scenario_index)`.")
    if !has_lower_bound(dvar)
        error("Decision variable $(dvar) does not have a lower bound.")
    end
    # Dispatch to structure
    return decision_dispatch(JuMP.lower_bound,
                             structure(owner_model(dvar)),
                             index(dvar),
                             stage(dvar))::Float64
end
"""
    lower_bound(dvar::DecisionVariable, scenario_index::Integer)

Return the lower bound of the scenario-dependent decision variable `dvar` at `scenario_index`. Error if one does not exist.
"""
function JuMP.lower_bound(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `lower_bound(dvar)`.")
    if !has_lower_bound(dvar, scenario_index)
        error("Decision variable $(dvar) at $scenario_index does not have a lower bound.")
    end
    # Dispatch to structure
    return scenario_decision_dispatch(JuMP.lower_bound,
                                      structure(owner_model(dvar)),
                                      index(dvar),
                                      stage(dvar),
                                      scenario_index)::Float64
end
"""
    has_upper_bound(dvar::DecisionVariable)

Return `true` if the first-stage decision variable `dvar` has a upper bound.
"""
function JuMP.has_upper_bound(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `upper_bound(dvar, scenario_index)`.")
    # Dispatch to structure
    return decision_dispatch(JuMP.has_upper_bound,
                             structure(owner_model(dvar)),
                             index(dvar),
                             stage(dvar))::Bool
end
"""
    has_upper_bound(dvar::DecisionVariable, scenario_index::Integer)

Return `true` if the scenario-dependent decision variable `dvar` has a upper bound at `scenario_index`.
"""
function JuMP.has_upper_bound(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `upper_bound(dvar)`.")
    # Dispatch to structure
    return scenario_decision_dispatch(JuMP.has_upper_bound,
                                      structure(owner_model(dvar)),
                                      index(dvar),
                                      stage(dvar),
                                      scenario_index)::Bool
end
"""
    LowerBoundRef(dvar::DecisionVariable)

Return a constraint reference to the upper bound constraint of the decision variable `dvar`.
Errors if one does not exist.
"""
function JuMP.UpperBoundRef(dvar::DecisionVariable)
    moi_ub =  MOI.ConstraintIndex{SingleDecision, MOI.LessThan{Float64}}
    return SPConstraintRef{moi_ub, ScalarShape}(owner_model(dvar),
                                                stage(dvar),
                                                moi_ub(index(dvar).value),
                                                ScalarShape())
end
"""
    set_upper_bound(dvar::DecisionVariable, upper::Number)

Set the upper bound of the first-stage decision variable `dvar` to `upper`. If one does not exist, create a new upper bound constraint.
"""
function JuMP.set_upper_bound(dvar::DecisionVariable, upper::Number)
    new_set = MOI.LessThan(convert(Float64, upper))
    # Update proxy
    proxy_ = proxy(owner_model(dvar), stage(dvar))
    dref = DecisionRef(proxy_, index(dvar))
    set_upper_bound(dref, upper)
    # Update structure
    if stage(dvar) == 1
        # Dispatch to structure
        decision_dispatch!(JuMP.set_upper_bound,
                           structure(owner_model(dvar)),
                           index(dvar),
                           stage(dvar),
                           upper)
    else
        # Update all constraints at stage
        for scenario_index in 1:num_scenarios(owner_model(dvar), stage(dvar))
            set_upper_bound(dvar, scenario_index, upper)
        end
    end
    return nothing
end
"""
    set_upper_bound(dvar::DecisionVariable, scenario_index::Integer, upper::Number)

Set the upper bound of the scenario-dependent decision variable `dvar` at `scenario_index` to `upper`. If one does not exist, create a new upper bound constraint.
"""
function JuMP.set_upper_bound(dvar::DecisionVariable, scenario_index::Integer, upper::Number)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `set_upper_bound(dvar, upper)`.")
    new_set = MOI.LessThan(convert(Float64, upper))
    # Dispatch to structure
    scenario_decision_dispatch!(JuMP.set_upper_bound,
                                structure(owner_model(dvar)),
                                index(dvar),
                                stage(dvar),
                                scenario_index,
                                upper)
    return nothing
end
"""
    delete_upper_bound(dvar::DecisionVariable)

Delete the upper bound constraint of the first-stage decision variable `dvar`.
"""
function JuMP.delete_upper_bound(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `delete_upper_bound(dvar, scenario_index)`.")
    # Update proxy
    proxy_ = proxy(owner_model(dvar), stage(dvar))
    dref = DecisionRef(proxy_, index(dvar))
    delete_upper_bound(dref)
    # Update structure
    if stage(dvar) == 1
        # Dispatch to structure
        decision_dispatch!(JuMP.delete_upper_bound,
                           structure(owner_model(dvar)),
                           index(dvar),
                           stage(dvar))
    else
        # Update all constraints at stage
        for scenario_index in 1:num_scenarios(owner_model(dvar), stage(dvar))
            delete_upper_bound(dvar, scenario_index)
        end
    end
end
"""
    delete_upper_bound(dvar::DecisionVariable, scenario_index::Integer)

Delete the upper bound constraint of the scenario-dependent decision variable `dvar` at `scenario_index`.
"""
function JuMP.delete_upper_bound(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `delete_upper_bound(dvar)`.")
    # Dispatch to structure
    scenario_decision_dispatch!(JuMP.delete_upper_bound,
                                structure(owner_model(dvar)),
                                index(dvar),
                                stage(dvar),
                                scenario_index)
end
"""
    upper_bound(dvar::DecisionVariable)

Return the upper bound of the first-stage decision variable `dvar`. Error if one does not exist.
"""
function JuMP.upper_bound(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `upper_bound(dvar, scenario_index)`.")
    if !has_upper_bound(dvar)
        error("Decision $(dvar) does not have a upper bound.")
    end
    # Dispatch to structure
    return decision_dispatch(JuMP.upper_bound,
                             structure(owner_model(dvar)),
                             index(dvar),
                             stage(dvar))::Float64
end
"""
    upper_bound(dvar::DecisionVariable, scenario_index::Integer)

Return the upper bound of the scenario-dependent decision variable `dvar` at `scenario_index`. Error if one does not exist.
"""
function JuMP.upper_bound(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `upper_bound(dvar)`.")
    if !has_upper_bound(dvar, scenario_index)
        error("Decision $(dvar) at `scenario_index` does not have a upper bound.")
    end
    # Dispatch to structure
    return scenario_decision_dispatch(JuMP.upper_bound,
                                      structure(owner_model(dvar)),
                                      index(dvar),
                                      stage(dvar),
                                      scenario_index)::Float64
end
"""
    is_integer(dvar::DecisionVariable)

Return `true` if the first-stage decision variable `dvar` is constrained to be integer.
"""
function JuMP.is_integer(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `is_integer(dvar, scenario_index)`.")
    # Dispatch to structure
    return decision_dispatch(JuMP.is_integer,
                             structure(owner_model(dvar)),
                             index(dvar),
                             stage(dvar))::Bool
end
"""
    is_integer(dvar::DecisionVariable, scenario_index::Integer)

Return `true` if the scenario-dependent decision variable `dvar` is constrained to be integer at `scenario_index`.
"""
function JuMP.is_integer(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `is_integer(dvar)`.")
    # Dispatch to structure
    return scenario_decision_dispatch(JuMP.is_integer,
                                      structure(owner_model(dvar)),
                                      index(dvar),
                                      stage(dvar),
                                      scenario_index)::Bool
end
"""
    IntegerRef(dvar::DecisionVariable)

Return a constraint reference to the integrality constraint of the decision variable `dvar`.
Errors if one does not exist.
"""
function JuMP.IntegerRef(dvar::DecisionVariable)
    moi_int =  MOI.ConstraintIndex{SingleDecision, MOI.Integer}
    return SPConstraintRef{moi_int, ScalarShape}(owner_model(dvar),
                                                 stage(dvar),
                                                 moi_int(index(dvar).value),
                                                 ScalarShape())
end
"""
    set_integer(dvar::DecisionVariable)

Add an integrality constraint on the first-stage decision variable `dvar`.
"""
function JuMP.set_integer(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `set_integer(dvar, scenario_index)`.")
    # Update proxy
    proxy_ = proxy(owner_model(dvar), stage(dvar))
    dref = DecisionRef(proxy_, index(dvar))
    set_integer(dref)
    # Update structure
    if stage(dvar) == 1
        # Dispatch to structure
        decision_dispatch!(JuMP.set_integer,
                           structure(owner_model(dvar)),
                           index(dvar),
                           stage(dvar))
    else
        # Update all constraints at stage
        for scenario_index in 1:num_scenarios(owner_model(dvar), stage(dvar))
            set_integer(dvar, scenario_index)
        end
    end
    return nothing
end
"""
    set_integer(dvar::DecisionVariable, scenario_index::Integer)

Add an integrality constraint on the scenario-dependent decision variable `dvar` at `scenario_index`.
"""
function JuMP.set_integer(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `set_integer(dvar)`.")
    # Dispatch to structure
    scenario_decision_dispatch!(JuMP.set_integer,
                                structure(owner_model(dvar)),
                                index(dvar),
                                stage(dvar),
                                scenario_index)
    return nothing
end
"""
    unset_integer(dvar::DecisionVariable)

Delete the integrality constraint of the first-stage decision variable `dvar`.
"""
function JuMP.unset_integer(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `unset_integer(dvar, scenario_index)`.")
    # Update proxy
    proxy_ = proxy(owner_model(dvar), stage(dvar))
    dref = DecisionRef(proxy_, index(dvar))
    unset_integer(dref)
    # Update structure
    if stage(dvar) == 1
        # Dispatch to structure
        decision_dispatch!(JuMP.unset_integer,
                           structure(owner_model(dvar)),
                           index(dvar),
                           stage(dvar))
    else
        # Update all constraints at stage
        for scenario_index in 1:num_scenarios(owner_model(dvar), stage(dvar))
            unset_integer(dvar, scenario_index)
        end
    end
    return nothing
end
"""
    unset_integer(dvar::DecisionVariable, scenario_index::Integer)

Delete the integrality constraint of the scenario-dependent decision variable `dvar` at `scenario_index`.
"""
function JuMP.unset_integer(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `unset_integer(dvar)`.")
    # Dispatch to structure
    scenario_decision_dispatch!(JuMP.unset_integer,
                                structure(owner_model(dvar)),
                                index(dvar),
                                stage(dvar),
                                scenario_index)
    return nothing
end
"""
    is_binary(dvar::DecisionVariable)

Return `true` if the first-stage decision variable `dvar` is constrained to be binary.
"""
function JuMP.is_binary(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `is_binary(dvar, scenario_index)`.")
    # Dispatch to structure
    return decision_dispatch(JuMP.is_binary,
                             structure(owner_model(dvar)),
                             index(dvar),
                             stage(dvar))::Bool
end
"""
    is_binary(dvar::DecisionVariable, scenario_index::Integer)

Return `true` if the scenario-dependent decision variable `dvar` is constrained to be binary at `scenario_index`.
"""
function JuMP.is_binary(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `is_binary(dvar)`.")
    # Dispatch to structure
    return scenario_decision_dispatch(JuMP.is_binary,
                                      structure(owner_model(dvar)),
                                      index(dvar),
                                      stage(dvar),
                                      scenario_index)::Bool
end
"""
    BinaryRef(dvar::DecisionVariable)

Return a constraint reference to the binary constraint of the decision variable `dvar`.
Errors if one does not exist.
"""
function JuMP.BinaryRef(dvar::DecisionVariable)
    moi_bin =  MOI.ConstraintIndex{SingleDecision, MOI.ZeroOne}
    return SPConstraintRef{moi_bin, ScalarShape}(owner_model(dvar),
                                                 stage(dvar),
                                                 moi_bin(index(dvar).value),
                                                 ScalarShape())
end
"""
    set_binary(dvar::DecisionVariable)

Constrain the first-stage decision variable `dvar` to the set ``\\{0,1\\}``.
"""
function JuMP.set_binary(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `set_binary(dvar, scenario_index)`.")
    # Update proxy
    proxy_ = proxy(owner_model(dvar), stage(dvar))
    dref = DecisionRef(proxy_, index(dvar))
    set_binary(dref)
    # Update structure
    if stage(dvar) == 1
        # Dispatch to structure
        decision_dispatch!(JuMP.set_binary,
                           structure(owner_model(dvar)),
                           index(dvar),
                           stage(dvar))
    else
        # Update all constraints at stage
        for scenario_index in 1:num_scenarios(owner_model(dvar), stage(dvar))
            set_binary(dvar, scenario_index)
        end
    end
    return nothing
end
"""
    set_binary(dvar::DecisionVariable, scenario_index::Integer)

Constrain the scenario-dependent decision variable `dvar` to the set ``\\{0,1\\}`` at `scenario_index`.
"""
function JuMP.set_binary(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `set_binary(dvar)`.")
    # Dispatch to structure
    scenario_decision_dispatch!(JuMP.set_binary,
                                structure(owner_model(dvar)),
                                index(dvar),
                                stage(dvar),
                                scenario_index)
end
"""
    unset_binary(dvar::DecisionVariable)

Delete the binary constraint of the first-stage decision variable `dvar`.
"""
function JuMP.unset_binary(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `unset_binary(dvar, scenario_index)`.")
    # Update proxy
    proxy_ = proxy(owner_model(dvar), stage(dvar))
    dref = DecisionRef(proxy_, index(dvar))
    unset_binary(dref)
    # Update structure
    if stage(dvar) == 1
        # Dispatch to structure
        decision_dispatch!(JuMP.unset_binary,
                           structure(owner_model(dvar)),
                           index(dvar),
                           stage(dvar))
    else
        # Update all constraints at stage
        for scenario_index in 1:num_scenarios(owner_model(dvar), stage(dvar))
            unset_binary(dvar, scenario_index)
        end
    end
    return nothing
end
"""
    unset_binary(dvar::DecisionVariable, scenario_index::Integer)

Delete the binary constraint of the scenario-dependent decision variable `dvar` at `scenario_index`.
"""
function JuMP.unset_binary(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `unset_binary(dvar)`.")
    # Dispatch to structure
    scenario_decision_dispatch!(JuMP.unset_binary,
                                structure(owner_model(dvar)),
                                index(dvar),
                                stage(dvar),
                                scenario_index)
    return nothing
end
"""
    start_value(dvar::DecisionVariable)

Return the start value of the first-stage decision variable `dvar`.
"""
function JuMP.start_value(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `start_value(dvar, scenario_index)`.")
    # Dispatch to structure
    decision_dispatch(JuMP.start_value,
                      structure(owner_model(dvar)),
                      index(dvar),
                      stage(dvar))
end
"""
    start_value(dvar::DecisionVariable, scenario_index::Integer)

Return the start value of the scenario-dependent decision variable `dvar`
at `scenario_index`.
"""
function JuMP.start_value(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `start_value(dvar)`.")
    # Dispatch to structure
    scenario_decision_dispatch(JuMP.start_value,
                               structure(owner_model(dvar)),
                               index(dvar),
                               stage(dvar),
                               scenario_index)
end
"""
    set_start_value(dvar::DecisionVariable)

Set the start value of the first-stage decision variable `dvar` to `value`.
"""
function JuMP.set_start_value(dvar::DecisionVariable, value::Number)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `set_start_value(dvar, scenario_index, value)`.")
    # Update proxy
    proxy_ = proxy(owner_model(dvar), stage(dvar))
    dref = DecisionRef(proxy_, index(dvar))
    set_start_value(dref, value)
    # Dispatch to structure
    decision_dispatch!(JuMP.set_start_value,
                       structure(owner_model(dvar)),
                       index(dvar),
                       stage(dvar),
                       value)
    return nothing
end
"""
    set_start_value(dvar::DecisionVariable, scenario_index::Integer, value::Number)

Set the start value of the scenario-dependent decision variable `dvar`
at `scenario_index` to `value`.
"""
function JuMP.set_start_value(dvar::DecisionVariable, scenario_index::Integer, value::Number)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `set_start_value(dvar, value)`.")
    # Dispatch to structure
    scenario_decision_dispatch!(JuMP.set_start_value,
                                structure(owner_model(dvar)),
                                index(dvar),
                                stage(dvar),
                                scenario_index,
                                value)
end

function Base.hash(dvar::DecisionVariable, h::UInt)
    return hash(objectid(owner_model(dvar)), hash(dvar.index, h))
end

JuMP.isequal_canonical(d::DecisionVariable, other::DecisionVariable) = isequal(d, other)
function Base.isequal(dvar::DecisionVariable, other::DecisionVariable)
    return owner_model(dvar) === owner_model(other) && dvar.index == other.index
end

Base.iszero(::DecisionVariable) = false
Base.copy(dvar::DecisionVariable) = DecisionVariable(dvar.stochasticprogram, stage(dvar), dvar.index)
Base.broadcastable(dvar::DecisionVariable) = Ref(dvar)

function DecisionRef(dvar::DecisionVariable)
    stage(dvar) > 1 && error("$dvar is scenario dependent, consider `DecisionRef(dvar, scenario_index)`.")
    sp = owner_model(dvar)
    return DecisionRef(structure(sp), index(dvar))
end
function DecisionRef(dvar::DecisionVariable, scenario_index::Integer)
    stage(dvar) == 1 && error("$dvar is not scenario dependent, consider `DecisionRef(dvar)`.")
    sp = owner_model(dvar)
    return DecisionRef(structure(sp), index(dvar), stage(dvar), scenario_index)
end
function DecisionRef(dvar::DecisionVariable, at_stage::Integer, scenario_index::Integer)
    at_stage > stage(dvar) || error("$dvar can only be known after stage $(stage(dvar)).")
    sp = owner_model(dvar)
    return DecisionRef(structure(sp), index(dvar), at_stage, stage(dvar), scenario_index)
end

is_decision_type(::Type{DecisionVariable}) = true

function JuMP.moi_function_type(::Type{DecisionVariable})
    return SingleDecision
end
function JuMP.moi_function_type(::Type{<:Vector{<:DecisionVariable}})
    return VectorOfDecisions
end

"""
    relax_integrality(stochasticprogram::StochasticProgram)

Modifies `stochasticprogram` to "relax" all binary and integrality constraints on decisions and auxiliary variables.

Returns a function that can be called without any arguments to restore the
original stochasticprogram. The behavior of this function is undefined if additional
changes are made to the affected decisions and variables in the meantime.
"""
function JuMP.relax_integrality(stochasticprogram::StochasticProgram)
    unrelax = relax_integrality(structure(stochasticprogram))
    return unrelax
end


text:
using ProjAssistant
@quickactivate

## -------------------------------------------------------------------
@time begin

    import SparseArrays
    using ArgParse

    import Chemostat_Nanchen2006
    const ChN = Chemostat_Nanchen2006

    const iJR = ChN.iJR904
    const Nd = ChN.NanchenData # experimental data
    const Bd = ChN.BegData    # cost data

    import Chemostat
    const Ch = Chemostat
    const ChU = Ch.Utils
    const ChSS = Ch.SteadyState
    const ChLP = Ch.LP
end

## ----------------------------------------------------------------------------
# arg settings
ARGSET = ArgParseSettings()
@ArgParse.add_arg_table! ARGSET begin
    "--ignore-cached"
        help = "Ingnore on disk version of data"
        action = :store_true
end

ARGS_DICT = ArgParse.parse_args(ARGSET)
if isinteractive()
    # dev values
    ignore_cached = true
else
    ignore_cached = ARGS_DICT["ignore-cached"]
end
@info("ARGS", ignore_cached)

## -------------------------------------------------------------------
# Helper function
include("1.0.1_functions.jl")

## -------------------------------------------------------------------
MODELS_FILE = procdir(iJR, "base_models.bson")
MODELS = ldat(() -> Dict{Any, Any}(), iJR, MODELS_FILE);

## -------------------------------------------------------------------
# file globals
use_cached = false

## -------------------------------------------------------------------
let
    compress(model) = ChU.compressed_model(model)

    ## --------------------------------------
    @info "COST_CONSTRAINED MODELS";

    base_model = load_raw_model()
    exchs = ChU.exchanges(base_model)

    rescale_bounds!(base_model)
    close_exchanges!(base_model, exchs)
    base_model = add_cost(base_model)
    base_model = reset_exchanges(base_model, exchs)

    MODELS["base_model"] = compress(base_model)

    ## --------------------------------------
    scale_factor = 1000.0
    # fva_preprocessed models
    for (exp, D) in Nd.val(:D) |> enumerate

        DAT = get!(MODELS, "fva_models", Dict())
        !ignore_cached && haskey(DAT, exp) && continue # cached

        fva_model = make_fva_model(base_model, exp, D;
            scale_factor
        )

        # storing
        DAT[exp] = compress(fva_model)

        # caching
        sdat(iJR, MODELS, MODELS_FILE);
        GC.gc()
    end

    ## --------------------------------------
    # MAX MODEL
    max_model = make_max_model(base_model; scale_factor)

    # saving
    MODELS["max_model"] = compress(max_model)
    sdat(iJR, MODELS, MODELS_FILE)

    println()
end

ext: jl


text:
###################################################
# DataCubes.jl
# A general array library for Julia.
# https://github.com/c-s/DataCubes.jl.git
# MIT Licensed
###################################################

module TestDataCubes

using FactCheck

subtests = ["test_ldict.jl", "test_dict_array.jl", "test_labeled_array.jl",
            "test_na.jl", "test_naarray_operators.jl",
            "test_array_util.jl", "test_dataframe_interface.jl",
            "test_join.jl", "test_select.jl", "test_sort.jl", "test_ungroup.jl",
            "test_setops.jl", "test_array_helper_functions.jl",
            "test_enumeration_array.jl"]

FactCheck.onlystats(true)

for subtest in subtests
  include(subtest)
end

exitstatus()

end


text:
module matrixcore

export StoM, MtoS, CascadeScattering, Repeating

using LinearAlgebra

function StoM(S::Matrix{T})::Matrix{<:Number} where T <: Number
    # Converts a 2x2 scattering matrix to a transfer matrix
    # See Saleh & Teich 3.ed eq.7.1-6
    # Verified manually
    @assert(size(S) == (2, 2))
    @inbounds begin
        t12 = S[1, 1]
        r21 = S[1, 2]
        r12 = S[2, 1]
        t21 = S[2, 2]
    end
    M = [
        t12*t21 - r12*r21 r21;
        -r12              1
    ]
    (1/t21) .* M
end

function MtoS(M::Matrix{T})::Matrix{<:Number} where T <: Number
    # Converts a 2x2 transfer matrix to a scattering matrix
    # See Saleh & Teich 3.ed eq.7.1-5
    # Verified by tests as inverfse of StoM
    StoM(M)
end

function CascadeScattering(layers::Vector)::Matrix{Number}
    # Generates the total S matrix for a layered system
    # Matrixes are given in the order the elements appear, the array is inversed inside the function
    accumulated = I

    for S in view(layers, length(layers):-1:1)
        M = StoM(S)
        accumulated *= M
    end

    MtoS(accumulated)
end

function Repeating(cell::Matrix{T}, repetition::Number)::Matrix{T} where T <: Number
    # Simple helpertool to make sure gratings are handleded correctly
    # This helper helps achieve the goal that transfer matrices are only handled in this module
    M = StoM(cell)
    MtoS(M ^ repetition)
end

end # matrixcore


text:
@static if VERSION < v"0.7-"
    const Nothing = Void
    const Cvoid   = Void
    using Missings

    firstindex(s::AbstractString) = start(s)
    lastindex(s::AbstractString) = endof(s)

else
    using Dates
end

@inline function undef_vector(::Type{T}, len::Integer) where {T}
    @static if VERSION < v"0.7-"
        Vector{T}(len)
    else
        Vector{T}(undef, len)
    end
end

#
# Code snippets from https://github.com/JuliaLang/Compat.jl
#

function _compat(ex::Expr)
    if VERSION < v"0.7.0-DEV.2562"
        if ex.head == :call && ex.args[1] == :finalizer
            ex.args[2], ex.args[3] = ex.args[3], ex.args[2]
        end
    end
    return Expr(ex.head, map(_compat, ex.args)...)
end

_compat(ex) = ex

macro compat(ex)
    esc(_compat(ex))
end


text:
"""
Interpolations.jl requires that one(::Type{T}) be defined for any data
type we want to interpolate. Rather than defining one(::Type{Vector}) here,
which might have unforeseen consequences in other packages, we'll create
a very simple wrapper type that just knows one() and *
"""
struct InterpolatableArray{A <: AbstractArray}
    data::A
end

one(::Type{InterpolatableArray{A}}) where {A} = 1
*(n::Number, a::InterpolatableArray) = n * a.data

normalize_configuration!(joint_type::JointType, q) = nothing
function normalize_configuration!(joint_type::QuaternionFloating, q)
    n = norm(q[1:4])
    for i = 1:4
        q[i] /= n
    end
end

"""
    animate(vis::Visualizer, mechanism::Mechanism{Float64},
            times::Vector{Float64},
            configurations::Vector{Vector{Float64}};
            fps::Float64=60, realtimerate::Float64=1.)

Animate the given mechanism passing through a time-coded series of
configurations by linearly interpolating the configuration vectors.
"""
function animate(vis::Visualizer, mechanism::Mechanism{Float64},
                 times::Vector{Float64},
                 configurations::AbstractVector{<:AbstractVector{Float64}};
                 fps::Float64 = 60., realtimerate::Float64 = 1.)
    @assert fps > 0
    @assert 0 < realtimerate < Inf

    state = MechanismState(mechanism)
    interp_values = [InterpolatableArray(c) for c in configurations]
    interpolated_configurations = interpolate((times,), interp_values, Gridded(Linear()))
    t0, tf = first(times), last(times)
    framenum = 0
    walltime0 = time()
    @throttle framenum while true
        t = min(tf, t0 + (time() - walltime0) * realtimerate)
        q = interpolated_configurations(t)
        set_configuration!(state, q)
        RigidBodyDynamics.normalize_configuration!(state)
        settransform!(vis, state)
        framenum += 1
        t == tf && break
    end max_rate = fps
end

animate(mechanism::Mechanism, times::Vector{Float64},
        configurations::Vector{Vector{Float64}}) =
    animate(Visualizer(mechanism), mechanism, times, configurations)


text:
using CNNVisualize, Metalhead

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

img = load_image("../images/original/cat_dog.png")
model = VGG19()

@testset "BackPropagation" begin
  @test_nowarn begin
    vbprop = VanillaBackprop(model.layers)
    grads = vbprop(img)
    pgrads, ngrads = positive_negative_saliency(grads[1][1])
  end

  @test_nowarn begin
    gbprop = GuidedBackprop(model.layers)
    grads = gbprop(img)
  end

  @test_nowarn begin
    deconv = Deconvolution(model.layers)
    grads = deconv(img, 2) # Test in case multiple visualizations are to be generated
  end

  save_gradient_images(grads[1][1], "../images/results/save1.png")
  save_grayscale_gradient(grads[2][1], "../images/results/save1.png")
end

@testset "Grad CAM" begin
  @test_nowarn begin
    gcam = GradCAM(model.layers, 12)
    grads = gcam(img)
    save_gradcam(grads[1][1], "../images/original/cat_dog.png", "../images/results/save2.png", "../images/results/save3.png", "../images/results/save4.png")
  end

  @test_warn "target should be less than the model length" GradCAM(model.layers, length(model.layer.layers) + 5)

  @test_nowarn begin
    ggcam = GuidedGradCAM(model.layers, 12)
    grads = ggcam(img)
  end
end


text:
#Make sure InspectDR does not interpolate in points not in data.
using InspectDR
using Colors

#==Input
===============================================================================#
#Constants
#-------------------------------------------------------------------------------
black = RGB24(0, 0, 0)
white = RGB24(1, 1, 1)
red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)

#Input data
#-------------------------------------------------------------------------------
npts = 335
tmax = 1
Δ = tmax/(npts-1)
t = collect(0:Δ:tmax)-pi*Δ
y = zeros(t)
y[div(length(t), 2)] = 1

#==Generate plot
===============================================================================#
plot = InspectDR.Plot2D(title="Sample Plot (λ)")
strip = plot.strips[1]
plot.xext_full = InspectDR.PExtents1D(0, NaN)
strip.yext_full = InspectDR.PExtents1D(-.2, 1.2)
plot.xres = 100 #coarse resolution
#plot.layout.grid = grid(vmajor=true, vminor=true, hmajor=true)

style = :dashdot #solid/dashdot/...
wfrm = add(plot, t, y, id="y")
	wfrm.line = line(style=:none)
	wfrm.line = line(style=:solid)
	wfrm.glyph = glyph(shape=:o, size=10, color=blue)

a = plot.annotation
a.xlabel = "Time (s)"
a.ylabels = ["Signal Voltage (V)"]

gplot = display(InspectDR.GtkDisplay(), plot)

:Test_Complete


text:
"""
CodeTracking can be thought of as an extension of InteractiveUtils, and pairs well with Revise.jl.

- `code_string`, `@code_string`: fetch the source code (as a string) for a method definition
- `code_expr`, `@code_expr`: fetch the expression for a method definition
- `definition`: a lower-level variant of the above
- `pkgfiles`: return information about the source files that define a package
- `whereis`: Return location information about methods (with Revise, it updates as you edit files)
- `signatures_at`: return the signatures of all methods whose definition spans the specified location
"""
module CodeTracking

using Base: PkgId
using Core: LineInfoNode
using Base.Meta: isexpr
using UUIDs
using InteractiveUtils

export code_expr, @code_expr, code_string, @code_string, whereis, definition, pkgfiles, signatures_at

# More recent Julia versions assign the line number to the line with the function declaration,
# not the first non-comment line of the body.
const line_is_decl = VERSION >= v"1.5.0-DEV.567"

include("pkgfiles.jl")
include("utils.jl")

### Global storage

# These values get populated by Revise

# `method_info[sig]` is either:
#   - `missing`, to indicate that the method cannot be located
#   - a list of `(lnn,ex)` pairs. In almost all cases there will be just one of these,
#     but "mistakes" in moving methods from one file to another can result in more than
#     definition. The last pair in the list is the currently-active definition.
const method_info = IdDict{Type,Union{Missing,Vector{Tuple{LineNumberNode,Expr}}}}()

const _pkgfiles = Dict{PkgId,PkgFiles}()

# Callback for method-lookup. `lookupfunc = method_lookup_callback[]` must have the form
#     ret = lookupfunc(method)
# where `ret` is either `nothing` or `(lnn, def)`. `lnn` is a LineNumberNode (or any valid
# input to `CodeTracking.fileline`) and `def` is the expression defining the method.
const method_lookup_callback = Ref{Any}(nothing)

# Callback for `signatures_at` (lookup by file/lineno). `lookupfunc = expressions_callback[]`
# must have the form
#    mod, exsigs = lookupfunc(id, relpath)
# where
#    id is the PkgId of the corresponding package
#    relpath is the path of the file from the basedir of `id`
#    mod is the "active" module at that point in the source
#    exsigs is a ex=>sigs dictionary, where `ex` is the source expression and `sigs`
#        a list of method-signatures defined by that expression.
const expressions_callback = Ref{Any}(nothing)

const juliabase = joinpath("julia", "base")
const juliastdlib = joinpath("julia", "stdlib", "v$(VERSION.major).$(VERSION.minor)")

### Public API

"""
    filepath, line = whereis(method::Method)

Return the file and line of the definition of `method`. The meaning of `line`
depends on the Julia version: on Julia 1.5 and higher it is the line number of
the method declaration, otherwise it is the first line of the method's body.
"""
function whereis(method::Method)
    file, line = String(method.file), method.line
    startswith(file, "REPL[") && return file, line
    lin = get(method_info, method.sig, nothing)
    if lin === nothing
        f = method_lookup_callback[]
        if f !== nothing
            try
                Base.invokelatest(f, method)
                lin = get(method_info, method.sig, nothing)
            catch
            end
        end
    end
    if lin === nothing || ismissing(lin)
    else
        file, line = fileline(lin[end][1])
    end
    file = maybe_fix_path(file)
    return file, line
end

"""
    loc = whereis(sf::StackFrame)

Return location information for a single frame of a stack trace.
If `sf` corresponds to a frame that was inlined, `loc` will be `nothing`.
Otherwise `loc` will be `(filepath, line)`.
"""
function whereis(sf::StackTraces.StackFrame)
    sf.linfo === nothing && return nothing
    return whereis(sf, sf.linfo.def)
end

"""
    filepath, line = whereis(lineinfo, method::Method)

Return the file and line number associated with a specific statement in `method`.
`lineinfo.line` should contain the line number of the statement at the time `method`
was compiled. The current location is returned.
"""
function whereis(lineinfo, method::Method)
    file, line1 = whereis(method)
    # We could be in an expanded macro. Apply the correction only if the filename checks out.
    # (We're not super-fastidious here because of symlinks and other path ambiguities)
    samefile = basename(file) == basename(String(lineinfo.file))
    if !samefile
        return maybe_fix_path(String(lineinfo.file)), lineinfo.line
    end
    return file, lineinfo.line - method.line + line1
end
function whereis(lineinfo::Core.LineInfoNode, method::Method)
    # With LineInfoNode we have certainty about whether we're in a macro expansion
    meth = lineinfo.method
    if isa(meth, WeakRef)
        meth = meth.value
    end
    if meth === Symbol("macro expansion")
        return maybe_fix_path(String(lineinfo.file)), lineinfo.line
    end
    file, line1 = whereis(method)
    return file, lineinfo.line - method.line + line1
end

"""
    sigs = signatures_at(filename, line)

Return the signatures of all methods whose definition spans the specified location.
Prior to Julia 1.5, `line` must correspond to a line in the method body
(not the signature or final `end`).

Returns `nothing` if there are no methods at that location.
"""
function signatures_at(filename::AbstractString, line::Integer)
    if !startswith(filename, "REPL[")
        filename = abspath(filename)
    end
    if occursin(juliabase, filename)
        rpath = postpath(filename, juliabase)
        id = PkgId(Base)
        return signatures_at(id, rpath, line)
    elseif occursin(juliastdlib, filename)
        rpath = postpath(filename, juliastdlib)
        spath = splitpath(rpath)
        libname = spath[1]
        project = Base.active_project()
        id = getpkgid(project, libname)
        return signatures_at(id, joinpath(spath[2:end]...), line)
    end
    if startswith(filename, "REPL[")
        id = PkgId("@REPL")
        return signatures_at(id, filename, line)
    end
    for (id, pkgfls) in _pkgfiles
        if startswith(filename, basedir(pkgfls)) || id.name == "Main"
            bdir = basedir(pkgfls)
            rpath = isempty(bdir) ? filename : relpath(filename, bdir)
            if rpath ∈ pkgfls.files
                return signatures_at(id, rpath, line)
            end
        end
    end
    throw(ArgumentError("$filename not found in internal data, perhaps the package is not loaded (or not loaded with `includet`)"))
end

"""
    sigs = signatures_at(mod::Module, relativepath, line)

For a package that defines module `mod`, return the signatures of all methods whose definition
spans the specified location. `relativepath` indicates the path of the file relative to
the packages top-level directory, e.g., `"src/utils.jl"`.
`line` must correspond to a line in the method body (not the signature or final `end`).

Returns `nothing` if there are no methods at that location.
"""
function signatures_at(mod::Module, relpath::AbstractString, line::Integer)
    id = PkgId(mod)
    return signatures_at(id, relpath, line)
end

function signatures_at(id::PkgId, relpath::AbstractString, line::Integer)
    expressions = expressions_callback[]
    expressions === nothing && error("cannot look up methods by line number, try `using Revise` before loading other packages")
    try
        for (mod, exsigs) in Base.invokelatest(expressions, id, relpath)
            for (ex, sigs) in exsigs
                lr = linerange(ex)
                lr === nothing && continue
                line ∈ lr && return sigs
            end
        end
    catch
    end
    return nothing
end

"""
    src, line1 = definition(String, method::Method)

Return a string with the code that defines `method`. Also return the first line of the
definition, including the signature (which may not be the same line number returned
by `whereis`). If the method can't be located (line number 0), then `definition`
instead returns `nothing.`

Note this may not be terribly useful for methods that are defined inside `@eval` statements;
see [`definition(Expr, method::Method)`](@ref) instead.

See also [`code_string`](@ref).
"""
function definition(::Type{String}, method::Method)
    file, line = whereis(method)
    line == 0 && return nothing
    src = src_from_file_or_REPL(file)
    src === nothing && return nothing
    src = replace(src, "\r"=>"")
    eol = isequal('\n')
    linestarts = Int[]
    istart = 1
    for i = 1:line-1
        push!(linestarts, istart)
        istart = findnext(eol, src, istart) + 1
    end
    ex, iend = Meta.parse(src, istart; raise=false)
    iend = prevind(src, iend)
    if isfuncexpr(ex, method.name)
        iend = min(iend, lastindex(src))
        return strip(src[istart:iend], '\n'), line
    end
    # The function declaration was presumably on a previous line
    lineindex = lastindex(linestarts)
    linestop = max(0, lineindex - 20)
    while !isfuncexpr(ex, method.name) && lineindex > linestop
        istart = linestarts[lineindex]
        try
            ex, iend = Meta.parse(src, istart)
        catch
        end
        lineindex -= 1
        line -= 1
    end
    lineindex <= linestop && return nothing
    return chomp(src[istart:iend-1]), line
end

"""
    ex = definition(Expr, method::Method)
    ex = definition(method::Method)

Return an expression that defines `method`. If the definition can't be found,
returns `nothing`.

See also [`code_expr`](@ref).
"""
function definition(::Type{Expr}, method::Method)
    file = String(method.file)
    def = startswith(file, "REPL[") ? nothing : get(method_info, method.sig, nothing)
    if def === nothing
        f = method_lookup_callback[]
        if f !== nothing
            try
                Base.invokelatest(f, method)
                def = get(method_info, method.sig, nothing)
            catch
            end
        end
    end
    return def === nothing || ismissing(def) ? nothing : copy(def[end][2])
end

definition(method::Method) = definition(Expr, method)

"""
    code_expr(f, types)

Returns the expression for the method definition for `f` with the specified types.

May return `nothing` if Revise isn't loaded. In such cases, calling
`Meta.parse(code_string(f, types))` can sometimes be an alternative.
"""
code_expr(f, t) = definition(Expr, which(f, t))
macro code_expr(ex0...)
    InteractiveUtils.gen_call_with_extracted_types_and_kwargs(__module__, :code_expr, ex0)
end

"""
    code_string(f, types)

Returns the code-string for the method definition for `f` with the specified types.
"""
code_string(f, t) = definition(String, which(f, t))[1]
macro code_string(ex0...)
    InteractiveUtils.gen_call_with_extracted_types_and_kwargs(__module__, :code_string, ex0)
end

"""
    info = pkgfiles(name::AbstractString)
    info = pkgfiles(name::AbstractString, uuid::UUID)

Return a [`CodeTracking.PkgFiles`](@ref) structure with information about the files that
define the package specified by `name` and `uuid`.
Returns `nothing` if this package has not been loaded.
"""
pkgfiles(name::AbstractString, uuid::UUID) = pkgfiles(PkgId(uuid, name))
function pkgfiles(name::AbstractString)
    project = Base.active_project()
    # The value returned by Base.project_deps_get depends on the Julia version
    id = isdefined(Base, :TOMLCache) && Base.VERSION < v"1.6.0-DEV.1180" ? Base.project_deps_get(project, name, Base.TOMLCache()) :
                                                                           Base.project_deps_get(project, name)
    (id == false || id === nothing) && error("no package ", name, " recognized")
    return isa(id, PkgId) ? pkgfiles(id) : pkgfiles(name, id)
end
pkgfiles(id::PkgId) = get(_pkgfiles, id, nothing)

"""
    info = pkgfiles(mod::Module)

Return a [`CodeTracking.PkgFiles`](@ref) structure with information about the files that
were loaded to define the package that defined `mod`.
"""
pkgfiles(mod::Module) = pkgfiles(PkgId(mod))

if ccall(:jl_generating_output, Cint, ()) == 1
    precompile(Tuple{typeof(setindex!), Dict{PkgId,PkgFiles}, PkgFiles, PkgId})
end

end # module


text:
"""
    tnmean(a, b)

Mean of the truncated standard normal distribution in the interval [a, b].
"""
function tnmean(a::Real, b::Real)
    if !(a ≤ b)
        return oftype(middle(a, b), NaN)
    elseif a == b
        return middle(a, b)
    elseif abs(a) > abs(b)
        return -tnmean(-b, -a)
    elseif isinf(a) && isinf(b)
        return zero(middle(a, b))
    end

    @assert a < b && abs(a) ≤ abs(b)
    @assert a ≤ 0 ≤ b || 0 < a < b

    Δ = (b - a) * middle(a, b)
    #Δ = one(Δm1) + Δm1

    if a ≤ 0 ≤ b
        m = √(2/π) * expm1(-Δ) * exp(-a^2 / 2) / erf(b/√2, a/√2)
    elseif 0 < a < b
        z = exp(-Δ) * erfcx(b/√2) - erfcx(a/√2)
        iszero(z) && return middle(a, b)
        m = √(2/π) * expm1(-Δ) / z
    end
    return clamp(m, a, b)
end

"""
    tnmean(a, b, μ, σ)

Mean of the truncated normal distribution, where μ, σ are the mean and standard
deviation of the untruncated distribution.
"""
function tnmean(a, b, μ, σ)
    α = (a - μ) / σ
    β = (b - μ) / σ
    return μ + tnmean(α, β) * σ
end

"""
    tnmean(a, b)

Mean of the truncated standard normal distribution.

    tnmean(a, b, μ, σ)

Mean of the truncated normal distribution, where μ, σ are the mean and standard
deviation of the untruncated distribution.
"""
tnmean(a, b) = tnmean(a, b)

"""
    tnmom1i(a, b)

Mean of the normal distribution with variance -1 and mean 0, truncated to [a,b].
"""
function tnmom1i(a::Real, b::Real)
    if !(-Inf < a ≤ b < Inf)
        return oftype(middle(a, b), NaN)
    elseif a == b
        return middle(a, b)
    elseif abs(a) > abs(b)
        return -tnmom1i(-b, -a)
    end

    @assert -Inf < a < b < Inf && abs(a) ≤ abs(b)
    @assert a ≤ 0 ≤ b || 0 < a < b

    Δ = (b - a) * middle(a, b)
    m = (1/√2) * expm1(-Δ) / (dawson(a/√2) * exp(-Δ) - dawson(b/√2))
    return clamp(m, a, b)
end

"""
    tnmom1i(a, b, μ, σ)

Mean of the normal distribution with variance -σ^2 and mean μ,
truncated to [a,b].
"""
tnmom1i(a, b, μ, σ) = μ + σ * tnmom1i((a - μ)/σ, (b - μ)/σ)


text:
#generate community
p = random_micrm_params(2,2,0.5)

#convert to ODESystem
@named sys = micrm_system(p)

@testset "MTK system" begin
    @test length(states(sys)) == 4
    @test length(parameters(sys)) == 14
end

#convert to problem
#define starting mass
u0 = fill(0.1, 4)
u0 = [states(sys)[i] => u0[i] for i = eachindex(u0)]

tspan = (0.0, 10.0) #define tspan

prob = ODEProblem(sys,u0,(0.0,100.0),[], jac = true)
sol = solve(prob)


text:
# TODO:
# * Force seeds.
# * Add MC progress fcn.
# * Connect to needs_init.
# * Make log name.
# * Show logged results.

# Monte Carlo runs
function mc(scenario::Scenario, n::Int)

    # Sets up each individual scenario (seed and log file) and runs it.
    function mchelper(scenario::Scenario, k::Int)
        println("Running MC variation ", k, " on worker ", myid())
        scenario.sim.seed = k
        scenario.sim.log_file = "out/mc" * string(k) * ".h5"
        simulate(scenario)
    end

    # Run the scenarios in parallel.
    pmap(k -> mchelper(scenario, k), 1:n)

end


text:
#

function jet(N::Integer)
	cm = RGB{Float64}[
	    RGB(
	        clamp(min(4x - 1.5, -4x + 4.5) ,0.0,1.0),
	        clamp(min(4x - 0.5, -4x + 3.5) ,0.0,1.0),
	        clamp(min(4x + 0.5, -4x + 2.5) ,0.0,1.0))
	    for x in linspace(0.0,1.0,N)]
	end

jet() = jet(64);

access_dict["jet"] = jet

max_stars_repo_name: lobingera/ColorMaps