text:
struct ScenarioProblems{S <: AbstractScenario} <: AbstractScenarioProblems{S}
    scenarios::Vector{S}
    problems::Vector{JuMP.Model}

    function ScenarioProblems(scenarios::Vector{S}) where S <: AbstractScenario
        # ScenarioProblems are initialized without any subproblems.
        # These are added during generation.
        return new{S}(scenarios, Vector{JuMP.Model}())
    end
end
ScenarioProblemChannel{S} = RemoteChannel{Channel{ScenarioProblems{S}}}
DecisionChannel = RemoteChannel{Channel{Decisions}}
struct DistributedScenarioProblems{S <: AbstractScenario} <: AbstractScenarioProblems{S}
    scenario_distribution::Vector{Int}
    scenarioproblems::Vector{ScenarioProblemChannel{S}}
    decisions::Vector{DecisionChannel}

    function DistributedScenarioProblems(scenario_distribution::Vector{Int},
                                         scenarioproblems::Vector{ScenarioProblemChannel{S}},
                                         decisions::Vector{DecisionChannel}) where S <: AbstractScenario
        return new{S}(scenario_distribution, scenarioproblems, decisions)
    end
end

function DistributedScenarioProblems(_scenarios::Vector{S}) where S <: AbstractScenario
    scenarioproblems = Vector{ScenarioProblemChannel{S}}(undef, nworkers())
    decisions = Vector{DecisionChannel}(undef, nworkers())
    (nscen, extra) = divrem(length(_scenarios), nworkers())
    start = 1
    stop = nscen + (extra > 0)
    scenario_distribution = zeros(Int, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            n = nscen + (extra > 0)
            scenarioproblems[i] = RemoteChannel(() -> Channel{ScenarioProblems{S}}(1), w)
            decisions[i] = RemoteChannel(() -> Channel{Decisions}(1), w)
            scenario_range = start:stop
            @async remotecall_fetch(
                w,
                scenarioproblems[i],
                _scenarios[scenario_range]) do sp, scenarios
                    put!(sp, ScenarioProblems(scenarios))
                end
            @async remotecall_fetch(
                w,
                decisions[i]) do channel
                    put!(channel, Decisions())
                end
            scenario_distribution[i] = n
            start = stop + 1
            stop += n
            stop = min(stop, length(_scenarios))
            extra -= 1
        end
    end
    return DistributedScenarioProblems(scenario_distribution, scenarioproblems, decisions)
end

ScenarioProblems(::Type{S}, instantiation) where S <: AbstractScenario = ScenarioProblems(Vector{S}(), instantiation)

function ScenarioProblems(scenarios::Vector{S}, ::Union{Vertical, Horizontal}) where S <: AbstractScenario
    ScenarioProblems(scenarios)
end

function ScenarioProblems(scenarios::Vector{S}, ::Union{DistributedVertical, DistributedHorizontal}) where S <: AbstractScenario
    DistributedScenarioProblems(scenarios)
end

# Base overloads #
# ========================== #
Base.getindex(sp::DistributedScenarioProblems, i::Integer) = sp.scenarioproblems[i]
# ========================== #

# MOI #
# ========================== #
function MOI.get(scenarioproblems::ScenarioProblems, attr::ScenarioDependentModelAttribute)
    MOI.get(backend(subproblem(scenarioproblems, attr.scenario_index)), attr.attr)
end
function MOI.get(scenarioproblems::DistributedScenarioProblems, attr::ScenarioDependentModelAttribute)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if attr.scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], attr.scenario_index - j, attr.attr) do sp, i, attr
                    MOI.get(backend(fetch(sp).problems[i]), attr)
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end
function MOI.get(scenarioproblems::ScenarioProblems, attr::ScenarioDependentVariableAttribute, index::MOI.VariableIndex)
    MOI.get(backend(subproblem(scenarioproblems, attr.scenario_index)), attr.attr, index)
end
function MOI.get(scenarioproblems::DistributedScenarioProblems, attr::ScenarioDependentVariableAttribute, index::MOI.VariableIndex)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if attr.scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], attr.scenario_index - j, attr.attr, index) do sp, i, attr, index
                    MOI.get(backend(fetch(sp).problems[i]), attr, index)
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end
function MOI.get(scenarioproblems::ScenarioProblems, attr::ScenarioDependentConstraintAttribute, ci::CI)
    MOI.get(backend(subproblem(scenarioproblems, attr.scenario_index)), attr.attr, ci)
end
function MOI.get(scenarioproblems::DistributedScenarioProblems, attr::ScenarioDependentConstraintAttribute, ci::CI)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if attr.scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], attr.scenario_index - j, attr.attr, ci) do sp, i, attr, ci
                    MOI.get(backend(fetch(sp).problems[i]), attr, ci)
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end

function MOI.set(scenarioproblems::ScenarioProblems, attr::MOI.AbstractModelAttribute, value)
    for problem in subproblems(scenarioproblems)
        MOI.set(backend(problem), attr, value)
    end
    return nothing
end
function MOI.set(scenarioproblems::DistributedScenarioProblems, attr::MOI.AbstractModelAttribute, value)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    @sync begin
        for (i,w) in enumerate(workers())
            @async remotecall_fetch(
                w, scenarioproblems[w-1], attr, value) do sp, attr, value
                    MOI.set(fetch(sp), attr, value)
                end
        end
    end
    return nothing
end
function MOI.set(scenarioproblems::ScenarioProblems, attr::ScenarioDependentModelAttribute, value)
    MOI.set(backend(subproblem(scenarioproblems, attr.scenario_index)), attr.attr, value)
    return nothing
end
function MOI.set(scenarioproblems::DistributedScenarioProblems, attr::ScenarioDependentModelAttribute, value)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j, attr.attr, value) do sp, i, attr, value
                    MOI.set(backend(fetch(sp).problems[i]), attr, value)
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end
function MOI.set(scenarioproblems::ScenarioProblems, attr::MOI.AbstractOptimizerAttribute, value)
    for problem in subproblems(scenarioproblems)
        MOI.set(backend(problem), attr, value)
    end
    return nothing
end
function MOI.set(scenarioproblems::DistributedScenarioProblems, attr::MOI.AbstractOptimizerAttribute, value)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    @sync begin
        for (i,w) in enumerate(workers())
            @async remotecall_fetch(
                w, scenarioproblems[w-1], attr, value) do sp, attr, value
                    MOI.set(fetch(sp), attr, value)
                end
        end
    end
    return nothing
end
function MOI.set(scenarioproblems::ScenarioProblems, attr::MOI.AbstractVariableAttribute,
                 index::MOI.VariableIndex, value)
    for problem in subproblems(scenarioproblems)
        MOI.set(backend(problem), attr, index, value)
    end
    return nothing
end
function MOI.set(scenarioproblems::DistributedScenarioProblems, attr::MOI.AbstractVariableAttribute,
                 index::MOI.VariableIndex, value)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    @sync begin
        for (i,w) in enumerate(workers())
            @async remotecall_fetch(
                w, scenarioproblems[w-1], attr, index, value) do sp, attr, index, value
                    MOI.set(fetch(sp), attr, index, value)
                end
        end
    end
    return nothing
end
function MOI.set(scenarioproblems::ScenarioProblems, attr::ScenarioDependentVariableAttribute,
                 index::MOI.VariableIndex, value)
    MOI.set(backend(subproblem(scenarioproblems, attr.scenario_index)), attr.attr, index, value)
    return nothing
end
function MOI.set(scenarioproblems::DistributedScenarioProblems, attr::ScenarioDependentVariableAttribute,
                 index::MOI.VariableIndex, value)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if attr.scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], attr.scenario_index - j, attr.attr, index, value) do sp, i, attr, index, value
                    MOI.set(backend(fetch(sp).problems[i]), attr, index, value)
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end
function MOI.set(scenarioproblems::ScenarioProblems, attr::MOI.AbstractConstraintAttribute, ci::MOI.ConstraintIndex, value)
    for problem in subproblems(scenarioproblems)
        MOI.set(backend(problem), attr, ci, value)
    end
    return nothing
end
function MOI.set(scenarioproblems::DistributedScenarioProblems, attr::MOI.AbstractConstraintAttribute,
                 ci::MOI.ConstraintIndex, value)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    @sync begin
        for (i,w) in enumerate(workers())
            @async remotecall_fetch(
                w, scenarioproblems[w-1], attr, ci, value) do sp, attr, ci, value
                    MOI.set(fetch(sp), attr, ci, value)
                end
        end
    end
    return nothing
end
function MOI.set(scenarioproblems::ScenarioProblems, attr::ScenarioDependentConstraintAttribute,
                 ci::MOI.ConstraintIndex, value)
    MOI.set(backend(subproblem(scenarioproblems, attr.scenario_index)), attr.attr, ci, value)
    return nothing
end
function MOI.set(scenarioproblems::DistributedScenarioProblems, attr::ScenarioDependentConstraintAttribute,
                 ci::MOI.ConstraintIndex, value)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if attr.scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], attr.scenario_index - j, attr.attr, ci, value) do sp, i, attr, ci, value
                    MOI.set(backend(fetch(sp).problems[i]), attr, ci, value)
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end

function MOI.is_valid(scenarioproblems::ScenarioProblems, index::MOI.VariableIndex, scenario_index::Integer)
    return MOI.is_valid(backend(subproblem(scenarioproblems, scenario_index)), index)
end
function MOI.is_valid(scenarioproblems::DistributedScenarioProblems, index::MOI.VariableIndex, scenario_index::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j, index) do sp, i, index
                    MOI.is_valid(backend(fetch(sp).problems[i]), index)
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end

function MOI.is_valid(scenarioproblems::ScenarioProblems, ci::MOI.ConstraintIndex, scenario_index::Integer)
    return MOI.is_valid(backend(subproblem(scenarioproblems, scenario_index)), ci)
end
function MOI.is_valid(scenarioproblems::DistributedScenarioProblems, ci::MOI.ConstraintIndex, scenario_index::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j, ci) do sp, i, ci
                    MOI.is_valid(backend(fetch(sp).problems[i]), ci)
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end

function MOI.add_constraint(scenarioproblems::ScenarioProblems, f::SingleDecision, s::MOI.AbstractSet, scenario_index::Integer)
    return MOI.add_constraint(backend(subproblem(scenarioproblems, scenario_index)), f, s)
end
function MOI.add_constraint(scenarioproblems::DistributedScenarioProblems, f::SingleDecision, s::MOI.AbstractSet, scenario_index::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j, f, s) do sp, i, f, s
                    MOI.add_constraint(backend(fetch(sp).problems[i]), f, s)
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end

function MOI.delete(scenarioproblems::ScenarioProblems, index::MOI.VariableIndex, scenario_index::Integer)
    subprob = subproblem(scenarioproblems, scenario_index)
    JuMP.delete(subprob, DecisionRef(subprob, index))
    return nothing
end
function MOI.delete(scenarioproblems::DistributedScenarioProblems, index::MOI.VariableIndex, scenario_index::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j, index) do sp, i, index
                    subprob = fetch(sp).problems[i]
                    JuMP.delete(subprob, DecisionRef(subprob, index))
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end
function MOI.delete(scenarioproblems::ScenarioProblems, indices::Vector{MOI.VariableIndex}, scenario_index::Integer)
    subprob = subproblem(scenarioproblems, scenario_index)
    JuMP.delete(subprob, DecisionRef.(subprob, indices))
    return nothing
end
function MOI.delete(scenarioproblems::DistributedScenarioProblems, indices::Vector{MOI.VariableIndex}, scenario_index::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j, indices) do sp, i, indices
                    subprob = fetch(sp).problems[i]
                    JuMP.delete(subprob, DecisionRef.(subprob, indices))
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end
function MOI.delete(scenarioproblems::ScenarioProblems, ci::MOI.ConstraintIndex, scenario_index::Integer)
    MOI.delete(backend(subproblem(scenarioproblems, scenario_index)), ci)
    return nothing
end
function MOI.delete(scenarioproblems::DistributedScenarioProblems, ci::MOI.ConstraintIndex, scenario_index::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j, ci) do sp, i, ci
                    MOI.delete(backend(fetch(sp).problems[i]), ci)
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end
function MOI.delete(scenarioproblems::ScenarioProblems, cis::Vector{<:MOI.ConstraintIndex}, scenario_index::Integer)
    MOI.delete(backend(subproblem(scenarioproblems, scenario_index)), cis)
    return nothing
end
function MOI.delete(scenarioproblems::DistributedScenarioProblems, cis::Vector{<:MOI.ConstraintIndex}, scenario_index::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j, cis) do sp, i, cis
                    MOI.delete(backend(fetch(sp).problems[i]), cis)
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end

# JuMP #
# ========================== #
function JuMP.fix(scenarioproblems::ScenarioProblems, index::MOI.VariableIndex, scenario_index::Integer, val::Number)
    subprob = scenarioproblems.problems[scenario_index]
    dref = DecisionRef(subprob, index)
    fix(dref, val)
    return nothing
end
function JuMP.fix(scenarioproblems::DistributedScenarioProblems, index::MOI.VariableIndex, scenario_index::Integer, val::Number)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j, index, val) do sp, i, index, val
                    subprob = fetch(sp).problems[i]
                    dref = DecisionRef(subprob, index)
                    fix(dref, val)
                    return nothing
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end
function JuMP.unfix(scenarioproblems::ScenarioProblems, index::MOI.VariableIndex, scenario_index::Integer)
    subprob = scenarioproblems.problems[scenario_index]
    dref = DecisionRef(subprob, index)
    unfix(dref)
    return nothing
end
function JuMP.unfix(scenarioproblems::DistributedScenarioProblems, index::MOI.VariableIndex, scenario_index::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j, index) do sp, i, index
                    subprob = fetch(sp).problems[i]
                    dref = DecisionRef(subprob, index)
                    unfix(dref)
                    return nothing
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end

function JuMP.objective_function_type(scenarioproblems::ScenarioProblems, scenario_index::Integer)
    subprob = subproblem(scenarioproblems, scenario_index)
    return jump_function_type(subprob, MOI.get(backend(subprob), MOI.ObjectiveFunctionType()))
end
function JuMP.objective_function_type(scenarioproblems::DistributedScenarioProblems, scenario_index::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j) do sp, i
                    s = fetch(sp).problems[i]
                    return jump_function_type(s, MOI.get(backend(s), MOI.ObjectiveFunctionType()))
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end

function JuMP.objective_function(scenarioproblems::ScenarioProblems,
                                 proxy::JuMP.Model,
                                 scenario_index::Integer,
                                 FunType::Type{<:AbstractJuMPScalar})
    MOIFunType = moi_function_type(FunType)
    subprob = subproblem(scenarioproblems, scenario_index)
    func = MOI.get(subprob, MOI.ObjectiveFunction{MOIFunType}())::MOIFunType
    return jump_function(proxy, func)
end
function JuMP.objective_function(scenarioproblems::DistributedScenarioProblems,
                                 proxy::JuMP.Model,
                                 scenario_index::Integer,
                                 FunType::Type{<:AbstractJuMPScalar})
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            f = remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j, FunType) do sp, i, FunType
                    MOIFunType = moi_function_type(FunType)
                    subprob = fetch(sp).problems[i]
                    func = MOI.get(subprob, MOI.ObjectiveFunction{MOIFunType}())::MOIFunType
                    return func
                end
            return jump_function(proxy, f)
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end

function JuMP._moi_optimizer_index(scenarioproblems::ScenarioProblems, index::MOI.VariableIndex, scenario_index::Integer)
    return decision_index(backend(subproblem(scenarioproblems, scenario_index)), index)
end

function JuMP._moi_optimizer_index(scenarioproblems::ScenarioProblems, ci::CI, scenario_index::Integer)
    return decision_index(backend(subproblem(scenarioproblems, scenario_index)), ci)
end

function JuMP.set_objective_coefficient(scenarioproblems::ScenarioProblems, index::VI, scenario_index::Integer, coeff::Real)
    subprob = subproblem(scenarioproblems, scenario_index)
    dref = DecisionRef(subprob, index)
    set_objective_coefficient(subprob, dref, coeff)
    return nothing
end
function JuMP.set_objective_coefficient(scenarioproblems::DistributedScenarioProblems, index::VI, scenario_index::Integer, coeff::Real)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j, index, coeff) do sp, i, index, coeff
                    subprob = fetch(sp).problems[i]
                    dref = DecisionRef(subprob, index)
                    set_objective_coefficient(subprob, dref, coeff)
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end

function JuMP.set_normalized_coefficient(scenarioproblems::ScenarioProblems,
                                         ci::CI{F,S},
                                         index::VI,
                                         scenario_index::Integer,
                                         value) where {T, F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}}, S}
    MOI.modify(backend(subproblem(scenarioproblems, scenario_index)), ci,
               DecisionCoefficientChange(index, convert(T, value)))
    return nothing
end
function JuMP.set_normalized_coefficient(scenarioproblems::DistributedScenarioProblems,
                                         ci::CI{F,S},
                                         index::VI,
                                         scenario_index::Integer,
                                         value) where {T, F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}}, S}
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j, ci, index, value) do sp, i, ci, index, value
                    MOI.modify(backend(fetch(sp).problems[i]), ci,
                               DecisionCoefficientChange(index, convert(T, value)))
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
    return nothing
end

function JuMP.normalized_coefficient(scenarioproblems::ScenarioProblems,
                                     ci::CI{F,S},
                                     index::VI,
                                     scenario_index::Integer) where {T, F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}}, S}
    subprob = subproblem(scenarioproblems, scenario_index)
    f = MOI.get(backend(subprob), MOI.ConstraintFunction(), ci)::F
    dref = DecisionRef(subprob, index)
    return JuMP._affine_coefficient(jump_function(subprob, f), dref)
end
function JuMP.normalized_coefficient(scenarioproblems::DistributedScenarioProblems,
                                     ci::CI{F,S},
                                     index::VI,
                                     scenario_index::Integer) where {T, F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}}, S}
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j, ci, index) do sp, i, ci, index
                    subprob = fetch(sp).problems[i]
                    f = MOI.get(backend(subprob), MOI.ConstraintFunction(), ci)::F
                    dref = DecisionRef(subprob, index)
                    return JuMP._affine_coefficient(jump_function(subprob, f), dref)
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end

function JuMP.set_normalized_rhs(scenarioproblems::ScenarioProblems,
                                 ci::CI{F,S},
                                 scenario_index::Integer,
                                 value) where {T,
                                               F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}},
                                               S <: MOIU.ScalarLinearSet{T}}
    MOI.set(backend(subproblem(scenarioproblems, scenario_index)), MOI.ConstraintSet(), ci,
            S(convert(T, value)))
    return nothing
end
function JuMP.set_normalized_rhs(scenarioproblems::DistributedScenarioProblems,
                                 ci::CI{F,S},
                                 scenario_index::Integer,
                                 value) where {T,
                                               F <: Union{AffineDecisionFunction{T}, QuadraticDecisionFunction{T}},
                                               S <: MOIU.ScalarLinearSet{T}}
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j, ci, value) do sp, i, ci, value
                    MOI.set(backend(fetch(sp).problems[i]), MOI.ConstraintSet(), ci,
                            S(convert(T, value)))
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end

function JuMP.delete(scenarioproblems::ScenarioProblems, index::MOI.VariableIndex)
    for subprob in subproblems(scenarioproblems)
        dref = DecisionRef(subprob, index)
        delete(subprob, dref)
    end
    return nothing
end
function JuMP.delete(scenarioproblems::DistributedScenarioProblems, index::MOI.VariableIndex)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    @sync begin
        for w in workers()
            @async remotecall_fetch(
                w, scenarioproblems[w-1], index) do sp, index
                    for subprob in fetch(sp).problems
                        dref = DecisionRef(subprob, index)
                        delete(subprob, dref)
                    end
                end
        end
    end
    return nothing
end
function JuMP.delete(scenarioproblems::ScenarioProblems, indices::Vector{MOI.VariableIndex})
    for subprob in subproblems(scenarioproblems)
        drefs = DecisionRef.(subprob, indices)
        delete(subprob, drefs)
    end
    return nothing
end
function JuMP.delete(scenarioproblems::DistributedScenarioProblems, indices::Vector{MOI.VariableIndex})
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    @sync begin
        for w in workers()
            @async remotecall_fetch(
                w, scenarioproblems[w-1], indices) do sp, indices
                    for subprob in fetch(sp).problems
                        drefs = DecisionRef.(subprob, indices)
                        delete(subprob, drefs)
                    end
                end
        end
    end
    return nothing
end

# Getters #
# ========================== #
function decision(scenarioproblems::ScenarioProblems, index::MOI.VariableIndex, scenario_index::Integer) where N
    subprob = subproblem(scenarioproblems, scenario_index)
    return decision(DecisionRef(subprob, index))
end
function decision(scenarioproblems::DistributedScenarioProblems, index::MOI.VariableIndex, scenario_index::Integer) where N
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j, index) do sp, i, index
                    subprob = fetch(sp).problems[i]
                    return decision(DecisionRef(subprob, index))
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end
function scenario(scenarioproblems::ScenarioProblems, scenario_index::Integer)
    return scenarioproblems.scenarios[scenario_index]
end
function scenario(scenarioproblems::DistributedScenarioProblems, scenario_index::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j) do sp, i
                    fetch(sp).scenarios[i]
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end
function scenarios(scenarioproblems::ScenarioProblems)
    return scenarioproblems.scenarios
end
function scenarios(scenarioproblems::DistributedScenarioProblems{S}) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    partial_scenarios = Vector{Vector{S}}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_scenarios[i] = remotecall_fetch(
                w, scenarioproblems[w-1]) do sp
                    fetch(sp).scenarios
                end
        end
    end
    return reduce(vcat, partial_scenarios)
end
function expected(scenarioproblems::ScenarioProblems)
    return expected(scenarioproblems.scenarios)
end
function expected(scenarioproblems::DistributedScenarioProblems{S}) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    partial_expecations = Vector{S}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_expecations[i] = remotecall_fetch(
                w, scenarioproblems[w-1]) do sp
                    expected(fetch(sp)).scenario
                end
        end
    end
    return expected(partial_expecations)
end
function scenario_type(scenarioproblems::ScenarioProblems{S}) where S <: AbstractScenario
    return S
end
function scenario_type(scenarioproblems::DistributedScenarioProblems{S}) where S <: AbstractScenario
    return S
end
function subproblem(scenarioproblems::ScenarioProblems, scenario_index::Integer)
    return scenarioproblems.problems[scenario_index]
end
function subproblem(scenarioproblems::DistributedScenarioProblems, scenario_index::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j) do sp, i
                    fetch(sp).problems[i]
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end
function subproblems(scenarioproblems::ScenarioProblems)
    return scenarioproblems.problems
end
function subproblems(scenarioproblems::DistributedScenarioProblems)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    partial_subproblems = Vector{Vector{JuMP.Model}}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_subproblems[i] = remotecall_fetch(
                w,scenarioproblems[w-1]) do sp
                    fetch(sp).problems
                end
        end
    end
    return reduce(vcat, partial_subproblems)
end
function num_subproblems(scenarioproblems::ScenarioProblems)
    return length(scenarioproblems.problems)
end
function num_subproblems(scenarioproblems::DistributedScenarioProblems)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    partial_lengths = Vector{Int}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_lengths[i] = remotecall_fetch(
                w,scenarioproblems[w-1]) do sp
                    num_subproblems(fetch(sp))
                end
        end
    end
    return sum(partial_lengths)
end
function decision_variables(scenarioproblems::ScenarioProblems)
    return scenarioproblems.decision_variables
end
function probability(scenarioproblems::ScenarioProblems, scenario_index::Integer)
    return probability(scenario(scenarioproblems, scenario_index))
end
function probability(scenarioproblems::DistributedScenarioProblems, scenario_index::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j) do sp, i
                    probability(fetch(sp).scenarios[i])
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end
function probability(scenarioproblems::ScenarioProblems)
    return probability(scenarioproblems.scenarios)
end
function probability(scenarioproblems::DistributedScenarioProblems)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    partial_probabilities = Vector{Float64}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_probabilities[i] = remotecall_fetch(
                w, scenarioproblems[w-1]) do sp
                    probability(fetch(sp))
                end
        end
    end
    return sum(partial_probabilities)
end
function num_scenarios(scenarioproblems::ScenarioProblems)
    return length(scenarioproblems.scenarios)
end
function num_scenarios(scenarioproblems::DistributedScenarioProblems)
    return sum(scenarioproblems.scenario_distribution)
end
distributed(scenarioproblems::ScenarioProblems) = false
distributed(scenarioproblems::DistributedScenarioProblems) = true
# ========================== #

# Setters
# ========================== #
function update_decision_state!(scenarioproblems::ScenarioProblems, index::MOI.VariableIndex, state::DecisionState)
    map(subproblems(scenarioproblems)) do subprob
        dref = DecisionRef(subprob, index)
        update_decision_state!(dref, state)
    end
    return nothing
end
function update_decision_state!(scenarioproblems::DistributedScenarioProblems, index::MOI.VariableIndex, state::DecisionState)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    @sync begin
        for (i,w) in enumerate(workers())
            @async remotecall_fetch(
                w, scenarioproblems[w-1], index, state) do sp, index, state
                    update_decision_state!(fetch(sp), index, state)
                end
        end
    end
    return nothing
end
function update_known_decisions!(scenarioproblems::ScenarioProblems, scenario_index::Integer)
    update_known_decisions!(subproblem(scenarioproblems, scenario_index))
    return nothing
end
function update_known_decisions!(scenarioproblems::DistributedScenarioProblems, scenario_index::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j) do sp, i
                    update_known_decisions!(fetch(sp).problems[i])
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end

function set_optimizer!(scenarioproblems::ScenarioProblems, optimizer)
    map(subproblems(scenarioproblems)) do subprob
        set_optimizer(subprob, optimizer)
    end
    return nothing
end
function set_optimizer!(scenarioproblems::DistributedScenarioProblems, optimizer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    @sync begin
        for (i,w) in enumerate(workers())
            @async remotecall_fetch(
                w, scenarioproblems[w-1], optimizer) do sp, opt
                    set_optimizer!(fetch(sp), opt)
                end
        end
    end
    return nothing
end
function add_scenario!(scenarioproblems::ScenarioProblems{S}, scenario::S) where S <: AbstractScenario
    push!(scenarioproblems.scenarios, scenario)
    return nothing
end
function add_scenario!(scenarioproblems::DistributedScenarioProblems{S}, scenario::S) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    _, w = findmin(scenarioproblems.scenario_distribution)
    add_scenario!(scenarioproblems, scenario, w+1)
    return nothing
end
function add_scenario!(scenarioproblems::DistributedScenarioProblems{S}, scenario::S, w::Integer) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    remotecall_fetch(
        w,
        scenarioproblems[w-1],
        scenario) do sp, scenario
            add_scenario!(fetch(sp), scenario)
        end
    scenarioproblems.scenario_distribution[w-1] += 1
    return nothing
end
function add_scenario!(scenariogenerator::Function, scenarioproblems::ScenarioProblems)
    add_scenario!(scenarioproblems, scenariogenerator())
    return nothing
end
function add_scenario!(scenariogenerator::Function, scenarioproblems::DistributedScenarioProblems)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    _, w = findmin(scenarioproblems.scenario_distribution)
    add_scenario!(scenariogenerator, scenarioproblems, w + 1)
    return nothing
end
function add_scenario!(scenariogenerator::Function, scenarioproblems::DistributedScenarioProblems, w::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    remotecall_fetch(
        w,
        scenarioproblems[w-1],
        scenariogenerator) do sp, generator
            add_scenario!(fetch(sp), generator())
        end
    scenarioproblems.scenario_distribution[w] += 1
    return nothing
end
function add_scenarios!(scenarioproblems::ScenarioProblems{S}, scenarios::Vector{S}) where S <: AbstractScenario
    append!(scenarioproblems.scenarios, scenarios)
    return nothing
end
function add_scenarios!(scenariogenerator::Function, scenarioproblems::ScenarioProblems{S}, n::Integer) where S <: AbstractScenario
    for i = 1:n
        add_scenario!(scenarioproblems) do
            return scenariogenerator()
        end
    end
    return nothing
end
function add_scenarios!(scenarioproblems::DistributedScenarioProblems{S}, scenarios::Vector{S}) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    (nscen, extra) = divrem(length(scenarios), nworkers())
    start = 1
    stop = 0
    @sync begin
        for w in workers()
            n = nscen + (extra > 0)
            stop += n
            stop = min(stop, length(scenarios))
            scenario_range = start:stop
            @async remotecall_fetch(
                w,
                scenarioproblems[w-1],
                scenarios[scenario_range]) do sp, scenarios
                    add_scenarios!(fetch(sp), scenarios)
                end
            scenarioproblems.scenario_distribution[w-1] += n
            start = stop + 1
            extra -= 1
        end
    end
    return nothing
end
function add_scenarios!(scenarioproblems::DistributedScenarioProblems{S}, scenarios::Vector{S}, w::Integer) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    remotecall_fetch(
        w,
        scenarioproblems[w-1],
        scenarios) do sp, scenarios
            add_scenarios!(fetch(sp), scenarios)
        end
    scenarioproblems.scenario_distribution[w-1] += length(scenarios)
    return nothing
end
function add_scenarios!(scenariogenerator::Function, scenarioproblems::DistributedScenarioProblems, n::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    (nscen, extra) = divrem(n, nworkers())
    @sync begin
        for w in workers()
            m = nscen + (extra > 0)
            @async remotecall_fetch(
                w,
                scenarioproblems[w-1],
                scenariogenerator,
                m) do sp, gen, n
                    add_scenarios!(gen, fetch(sp), n)
                end
            scenarioproblems.scenario_distribution[w-1] += m
            extra -= 1
        end
    end
    return nothing
end
function add_scenarios!(scenariogenerator::Function, scenarioproblems::DistributedScenarioProblems, n::Integer, w::Integer)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    remotecall_fetch(
        w,
        scenarioproblems[w-1],
        scenariogenerator,
        n) do sp, gen
            add_scenarios!(gen, fetch(sp), n)
        end
    scenarioproblems.scenario_distribution[w-1] += n
    return nothing
end
function clear_scenarios!(scenarioproblems::ScenarioProblems)
    empty!(scenarioproblems.scenarios)
    return nothing
end
function clear_scenarios!(scenarioproblems::DistributedScenarioProblems)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    @sync begin
        for w in workers()
            @async remotecall_fetch(
                w, scenarioproblems[w-1]) do sp
                    remove_scenarios!(fetch(sp))
                end
            scenarioproblems.scenario_distribution[w-1] = 0
        end
    end
    return nothing
end
function clear!(scenarioproblems::ScenarioProblems)
    map(scenarioproblems.problems) do subprob
        # Clear decisions
        if haskey(subprob.ext, :decisions)
            map(clear!, subprob.ext[:decisions])
        end
        # Clear model
        empty!(subprob)
    end
    empty!(scenarioproblems.problems)
    return nothing
end
function clear!(scenarioproblems::DistributedScenarioProblems)
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    @sync begin
        for w in workers()
            @async remotecall_fetch(
                w, scenarioproblems[w-1]) do sp
                    clear!(fetch(sp))
                end
        end
    end
    return nothing
end

function cache_solution!(cache::Dict{Symbol,SolutionCache}, scenarioproblems::ScenarioProblems, optimizer::MOI.AbstractOptimizer, stage::Integer)
    for i in 1:num_scenarios(scenarioproblems)
        key = Symbol(:node_solution_, stage, :_, i)
        subprob = subproblem(scenarioproblems, i)
        cache[key] = SolutionCache(backend(subprob))
        cache_model_attributes!(cache[key], optimizer, stage, i)
        variables = MOI.get(backend(subprob), MOI.ListOfVariableIndices())
        cache_variable_attributes!(cache[key], optimizer, variables, stage, i)
        ctypes = filter(t -> is_decision_type(t[1]), MOI.get(backend(subprob), MOI.ListOfConstraints()))
        constraints = mapreduce(vcat, ctypes) do (F, S)
            return MOI.get(backend(subprob), MOI.ListOfConstraintIndices{F,S}())
        end
        cache_constraint_attributes!(cache[key], optimizer, constraints, stage, i)
    end
end
function cache_solution!(cache::Dict{Symbol,SolutionCache}, scenarioproblems::DistributedScenarioProblems, optimizer::MOI.AbstractOptimizer, stage::Integer)
    for scenario_index in 1:num_scenarios(scenarioproblems)
        key = Symbol(:node_solution_, stage, :_, scenario_index)
        cache[key], variables, constraints = _prepare_subproblem_cache(scenarioproblems, scenario_index)
        cache_model_attributes!(cache[key], optimizer, stage, scenario_index)
        cache_variable_attributes!(cache[key], optimizer, variables, stage, scenario_index)
        cache_constraint_attributes!(cache[key], optimizer, constraints, stage, scenario_index)
    end
end
function _prepare_subproblem_cache(scenarioproblems::DistributedScenarioProblems, scenario_index::Integer)
    j = 0
    for w in workers()
        n = scenarioproblems.scenario_distribution[w-1]
        if scenario_index <= n + j
            return remotecall_fetch(
                w, scenarioproblems[w-1], scenario_index - j) do sp, i
                    subprob = fetch(sp).problems[i]
                    subcache = SolutionCache(backend(subprob))
                    variables = MOI.get(backend(subprob), MOI.ListOfVariableIndices())
                    ctypes = filter(t -> is_decision_type(t[1]), MOI.get(backend(subprob), MOI.ListOfConstraints()))
                    constraints = mapreduce(vcat, ctypes) do (F, S)
                        return MOI.get(backend(subprob), MOI.ListOfConstraintIndices{F,S}())
                    end
                    return subcache, variables, constraints
                end
        end
        j += n
    end
    throw(BoundsError(scenarioproblems, scenario_index))
end
# ========================== #

# Sampling #
# ========================== #
function sample!(scenarioproblems::ScenarioProblems{S}, sampler::AbstractSampler{S}, n::Integer) where S <: AbstractScenario
    _sample!(scenarioproblems, sampler, n, num_scenarios(scenarioproblems), 1/n)
end
function sample!(scenarioproblems::ScenarioProblems{S}, sampler::AbstractSampler{Scenario}, n::Integer) where S <: AbstractScenario
    _sample!(scenarioproblems, sampler, n, num_scenarios(scenarioproblems), 1/n)
end
function sample!(scenarioproblems::DistributedScenarioProblems{S}, sampler::AbstractSampler{S}, n::Integer) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    m = nscenarios(scenarioproblems)
    (nscen, extra) = divrem(n, nworkers())
    @sync begin
        for w in workers()
            d = nscen + (extra > 0)
            @async remotecall_fetch(
                w,
                scenarioproblems[w-1],
                sampler,
                d,
                m,
                1/n) do sp, sampler, n, m, π
                    _sample!(fetch(sp), sampler, n, m, π)
                end
            scenarioproblems.scenario_distribution[w-1] += d
            extra -= 1
        end
    end
    return nothing
end
function sample!(scenarioproblems::DistributedScenarioProblems{S}, sampler::AbstractSampler{Scenario}, n::Integer) where S <: AbstractScenario
    isempty(scenarioproblems.scenarioproblems) && error("No remote scenario problems.")
    m = nscenarios(scenarioproblems)
    (nscen, extra) = divrem(n, nworkers())
    @sync begin
        for w in workers()
            d = nscen + (extra > 0)
            @async remotecall_fetch(
                w,
                scenarioproblems[w-1],
                sampler,
                d,
                m,
                1/n) do sp, sampler, n, m, π
                    _sample!(fetch(sp), sampler, n, m, π)
                end
            scenarioproblems.scenario_distribution[w-1] += d
            extra -= 1
        end
    end
    return nothing
end
function _sample!(scenarioproblems::ScenarioProblems, sampler::AbstractSampler, n::Integer, m::Integer, π::AbstractFloat)
    if m > 0
        # Rescale probabilities of existing scenarios
        for scenario in scenarioproblems.scenarios
            p = probability(scenario) * m / (m+n)
            set_probability!(scenario, p)
        end
        π *= n/(m+n)
    end
    for i = 1:n
        add_scenario!(scenarioproblems) do
            return sample(sampler, π)
        end
    end
    return nothing
end
# ========================== #


text:
function cost_equation(cur_reactor::AbstractReactor, cur_value::Number)
  cur_G_T = cur_value

  cur_G_T *= K_F(cur_reactor)

  cur_G_T *= cur_reactor.sigma_v

  cur_G_T /= K_WM(cur_reactor)

  cur_G_T ^= (1/2)

  cur_G_T *= K_n(cur_reactor)

  cur_equation = Equation(2, 1, -1, cur_G_T)

  cur_equation
end

ext: jl


text:
module TestAvgFilter2D

using TapirBenchmarks
using Test

@testset for n in [3, 5, 100, 2^8]
    ys0, xs = avgfilter2d_setup(n)
    ys1 = avgfilter2d_seq!(copy(ys0), xs)
    @test avgfilter2d_threads!(copy(ys0), xs) == ys1
    @test avgfilter2d_tapir_dac!(copy(ys0), xs) == ys1
    @test avgfilter2d_tapir_seq!(copy(ys0), xs) == ys1
end

end


text:
immutable Attribute{T, Read, Write}
    name::Symbol
    stream::IOStream
end

function call{T, Read, Write}(::Type{Attribute{T, Read, Write}}, brick::Brick, relative_path::AbstractString, name::Symbol)
    path = joinpath(brick.root_path, relative_path, string(name))
    isfile(path) || error("file not found: $path")
    Read && (isreadable(path) || error("read access for attribute $(path) was requested, but file is not readable"))
    Write && (iswritable(path) || error("write access for attribute $(path) was requested, but file is not writable"))
    stream = open(path, Read, Write, false, Write, false)
    Attribute{T, Read, Write}(name, stream)
end

typealias ReadOnly{T} Attribute{T, true, false}
typealias WriteOnly{T} Attribute{T, false, true}
typealias ReadWrite{T} Attribute{T, true, true}

parse(::Type{ASCIIString}, s::ASCIIString) = s
parse(::Type{Vector{ASCIIString}}, s::ASCIIString) = split(s)

call{T, X}(attr::Attribute{T, true, X}) = parse(T, chomp(strip(readline(seekstart(attr.stream)), '\0')))
function call{T, X}(attr::Attribute{T, X, true}, value::T)
    write(seekstart(attr.stream), string(value))
    flush(attr.stream)
end


text:
module RadiativeEnergyFlux

using ...VariableTemplates
using DocStringExtensions

export AbstractNetSwFluxModel,
    PrescribedSwFluxAndAlbedo,
    PrescribedNetSwFlux,
    compute_net_sw_flux,
    compute_net_radiative_energy_flux

"""
   AbstractNetSwFluxModel{FT <: AbstractFloat}
"""
abstract type AbstractNetSwFluxModel{FT <: AbstractFloat} end

"""
    PrescribedSwFluxAndAlbedo{FT, FN1, FN2} <: AbstractNetSwFluxModel{FT}

Structure which contains functions for shortwave albedo and
shortwave fluxes. They are user-defined, constant across the domain
and can be a function of time.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct PrescribedSwFluxAndAlbedo{FT, FN1, FN2} <: AbstractNetSwFluxModel{FT}
    "Shortwave albedo in grid. No units."
    α::FN1
    "Shortwave flux in grid. Units of J m-2 s-1"
    swf::FN2
    function PrescribedSwFluxAndAlbedo(
        ::Type{FT};
        α::FN1 = (t) -> FT(NaN),
        swf::FN2 = (t) -> FT(NaN),
    ) where {FT, FN1, FN2}
        args = (α, swf)
        return new{FT, FN1, FN2}(args...)
    end
end

"""
    PrescribedNetSwFlux{FT,FN3} <: AbstractNetSwFluxModel{FT}

Structure which contains net shortwave flux values. The function is
user-defined, constant across the domain and can be a function of time.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct PrescribedNetSwFlux{FT, FN3} <: AbstractNetSwFluxModel{FT}
    "Net shortwave flux. Units of J m-2 s-1"
    nswf::FN3
    function PrescribedNetSwFlux(
        ::Type{FT};
        nswf::FN3 = (t) -> FT(NaN),
    ) where {FT, FN3}
        return new{FT, FN3}(nswf)
    end
end

"""
    compute_net_sw_flux(
        nswf_model::PrescribedSwFluxAndAlbedo{FT},
        t::Real,
    ) where {FT}

Computes the net shortwave flux as a function of time.
"""
function compute_net_sw_flux(
    nswf_model::PrescribedSwFluxAndAlbedo{FT},
    t::Real,
) where {FT}
    net_sw_flux = (FT(1) - nswf_model.α(t)) * nswf_model.swf(t)
    return net_sw_flux
end

"""
    compute_net_sw_flux(
        nswf_model:: PrescribedNetSwFlux{FT},
        t::Real
    ) where {FT}

Computes the net shortwave flux as a function of time.
"""
function compute_net_sw_flux(
    nswf_model::PrescribedNetSwFlux{FT},
    t::Real,
) where {FT}
    net_sw_flux = nswf_model.nswf(t)
    return net_sw_flux
end

"""
    compute_net_radiative_energy_flux(
        nswf_model::AbstractNetSwFluxModel{FT},
        t::Real
    )

Returns the net radiative energy flux as a function of time. Will
include long wave fluxes in future version.
"""
function compute_net_radiative_energy_flux(
    nswf_model::AbstractNetSwFluxModel{FT},
    t::Real,
) where {FT}
    net_radiative_flux = compute_net_sw_flux(nswf_model, t)
    return net_radiative_flux
end

end


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#
using DiscreteEvents, Random

mutable struct Server
  id::Int
  name::String
  input::Channel
  output::Channel
  op     # operation to take
  token  # current token

  Server(id, name, input, output, op) = new(id, name, input, output, op, nothing)
end

A = []

function take(S::Server)
    if isready(S.input)
        S.token = take!(S.input)
        push!(A, (tau(), S.name, S.id, S.token))
        event!(fun(put, S), after, rand())         # call put after some time
    else
        event!(fun(take, S), fun(isready, S.input)) # call again if input is ready
    end
end

function put(S::Server)
    put!(S.output, S.op(S.id, S.token))
    S.token = nothing
    take(S)
end

resetClock!(𝐶)
Random.seed!(123)

ch1 = Channel(32)  # create two channels
ch2 = Channel(32)

si = shuffle(1:8)
for i in 1:2:8
    take(Server(si[i], "foo", ch1, ch2, +))
    take(Server(si[i+1], "bar", ch2, ch1, *))
end

put!(ch1, 1) # put first token into channel 1
run!(𝐶, 10)

@test length(A) > 20
p = [i[3] for i in A]
for i in 1:8
    @test i ∈ p  # all processes did something
end


text:
# Keeping these dependencies separate as the are slow to load
using Plots
using VisualRegressionTests

@testset "plotting" begin
    @testset "Interval{Float64,$L,$R}" for (L, R) in BOUND_PERMUTATIONS
        intervals = [Interval{L,R}(float(x), float(x + 0.5)) for x in 1:11]
        plot_file = "interval_$(lowercase("$(L)_$(R)")).png"
        @plottest plot(intervals, 1:11) "references/$plot_file" false 0.01

        @testset "scatter" begin
            # Earlier versions of this functionality showed only end-points if plotted in
            # scatter, but for intervals the connect-line is part of the "marker"
            @plottest scatter(intervals, 1:11) "references/$plot_file" false 0.01
        end
    end

    @testset "DateTime intervals" begin
        start_dt = DateTime(2017, 1 ,1, 0, 0)
        end_dt = DateTime(2017, 1, 1, 10, 30)
        datetimes = start_dt:Hour(1):end_dt

        @testset "Interval{DateTime}" begin
            date_intervals = [dt .. (dt + Hour(1)) for dt in datetimes]
            @plottest plot(date_intervals, 1:11) "references/interval_datetime.png" false 0.01
        end

        @testset "AnchoredInterval" begin
            @plottest plot(HE.(datetimes), 1:11) "references/HE.png" false 0.01
            @plottest plot(HB.(datetimes), 1:11) "references/HB.png" false 0.01
        end
    end
end


text:
using DFTK
using WriteVTK

function build_magnetic_moments(atoms::AbstractArray; magmoms...)
    magmoms = Dict(magmoms)
    map(atoms) do (element, positions)
        element => fill(magmoms[element.symbol], length(positions))
    end
end

function attach_pseudos(atoms::AbstractArray; pseudomap...)
    pseudomap = Dict(pseudomap)

    map(atoms) do (element, position)
        pspfile = get(pseudomap, element.symbol, nothing)
        ElementPsp(element.symbol, psp=load_psp(pspfile)) => position
    end
end

lattice = load_lattice(joinpath(@__DIR__, "./Mn3Si.in"))
atoms   = load_atoms(  joinpath(@__DIR__, "./Mn3Si.in"))
atoms   = attach_pseudos(atoms, Mn="hgh/pbe/mn-q15.hgh",
                                Si="hgh/pbe/si-q4.hgh")
magnetic_moments = build_magnetic_moments(atoms; Si=0.0, Mn=5.0)

model = model_PBE(lattice, atoms; magnetic_moments=magnetic_moments,
                  smearing=Smearing.Gaussian(), temperature=0.01)
basis = PlaneWaveBasis(model; Ecut=45, kgrid=[13, 13, 13])

determine_diagtol = DFTK.ScfDiagtol(; ratio_ρdiff=0.01, diagtol_max=1e-3, diagtol_min=1e-12)
ρ = guess_density(basis, magnetic_moments)
scfres = DFTK.scf_potential_mixing(basis; tol=1e-11, determine_diagtol=determine_diagtol,
                                   mixing=KerkerMixing(), damping=DFTK.FixedDamping(0.7),
                                   diag_miniter=2, ρ=ρ)

save_scfres("Mn3Si.vts", scfres; save_ψ=false);

max_stars_repo_name: mfherbst/supporting-adaptive-damping

max_issues_repo_name: mfherbst/supporting-adaptive-damping

max_forks_repo_name: mfherbst/supporting-adaptive-damping


text:
module TypeEmulator

export
    emulate, MainModule, CoreModule, NewEmpty,
    NewModule, NewDataType, NewAny, NewBottom, NewFunction,
    NewType, NewTypeName, NewTypeVar, NewUnionAll, NewUnion,
    NewFunction, NewMethod, NewMethodList, NewMethodTable, NewTypeMapEntry,
    supertypes, isnewsubtype, linearize, register!

import Base:    show, subtypes, supertype

include("metadata.jl")
include("types.jl")
include("typeconversions.jl")
include("linear.jl")
include("relations.jl")
include("variables.jl")

end # module


text:
import AoC2019
using Test
using LightGraphs, SimpleWeightedGraphs
#using Plots
input_data = read(open("input20.txt"), String)
test_input = """
         A
         A
  #######.#########
  #######.........#
  #######.#######.#
  #######.#######.#
  #######.#######.#
  #####  B    ###.#
BC...##  C    ###.#
  ##.##       ###.#
  ##...DE  F  ###.#
  #####    G  ###.#
  #########.#####.#
DE..#######...###.#
  #.#########.###.#
FG..#########.....#
  ###########.#####
             Z
             Z
"""
test_input2 = """
                   A
                   A
  #################.#############
  #.#...#...................#.#.#
  #.#.#.###.###.###.#########.#.#
  #.#.#.......#...#.....#.#.#...#
  #.#########.###.#####.#.#.###.#
  #.............#.#.....#.......#
  ###.###########.###.#####.#.#.#
  #.....#        A   C    #.#.#.#
  #######        S   P    #####.#
  #.#...#                 #......VT
  #.#.#.#                 #.#####
  #...#.#               YN....#.#
  #.###.#                 #####.#
DI....#.#                 #.....#
  #####.#                 #.###.#
ZZ......#               QG....#..AS
  ###.###                 #######
JO..#.#.#                 #.....#
  #.#.#.#                 ###.#.#
  #...#..DI             BU....#..LF
  #####.#                 #.#####
YN......#               VT..#....QG
  #.###.#                 #.###.#
  #.#...#                 #.....#
  ###.###    J L     J    #.#.###
  #.....#    O F     P    #.#...#
  #.###.#####.#.#####.#####.###.#
  #...#.#.#...#.....#.....#.#...#
  #.#####.###.###.#.#.#########.#
  #...#.#.....#...#.#.#.#.....#.#
  #.###.#####.###.###.#.#.#######
  #.#.........#...#.............#
  #########.###.###.#############
           B   J   C
           U   P   P               """

test_input3 = """
             Z L X W       C
             Z P Q B       K
  ###########.#.#.#.#######.###############
  #...#.......#.#.......#.#.......#.#.#...#
  ###.#.#.#.#.#.#.#.###.#.#.#######.#.#.###
  #.#...#.#.#...#.#.#...#...#...#.#.......#
  #.###.#######.###.###.#.###.###.#.#######
  #...#.......#.#...#...#.............#...#
  #.#########.#######.#.#######.#######.###
  #...#.#    F       R I       Z    #.#.#.#
  #.###.#    D       E C       H    #.#.#.#
  #.#...#                           #...#.#
  #.###.#                           #.###.#
  #.#....OA                       WB..#.#..ZH
  #.###.#                           #.#.#.#
CJ......#                           #.....#
  #######                           #######
  #.#....CK                         #......IC
  #.###.#                           #.###.#
  #.....#                           #...#.#
  ###.###                           #.#.#.#
XF....#.#                         RF..#.#.#
  #####.#                           #######
  #......CJ                       NM..#...#
  ###.#.#                           #.###.#
RE....#.#                           #......RF
  ###.###        X   X       L      #.#.#.#
  #.....#        F   Q       P      #.#.#.#
  ###.###########.###.#######.#########.###
  #.....#...#.....#.......#...#.....#.#...#
  #####.#.###.#######.#######.###.###.#.#.#
  #.......#.......#.#.#.#.#...#...#...#.#.#
  #####.###.#####.#.#.#.#.###.###.#.###.###
  #.......#.....#.#...#...............#...#
  #############.#.#.###.###################
               A O F   N
               A A D   M                     """

function find_token(arr, position)
  dims = size(arr)
  pos_vectors = [CartesianIndex(0,1), CartesianIndex(1,0)]
  neg_vectors = [CartesianIndex(0,-1), CartesianIndex(-1,0)]
  token = ""
  for vec in pos_vectors
    if isuppercase(arr[(position + vec)])
      token = arr[(position + vec)] * arr[(position + 2 * vec)]
    end
  end
  for vec in neg_vectors
    if isuppercase(arr[(position + vec)])
      token = arr[(position + 2 * vec)] * arr[(position + vec)]
    end
  end
  if token != ""
    if position[1] == 3 || position[1] == dims[1] - 2 || position[2] == 3 || position[2] == dims[2] - 2
      return (token, false)
    else
      return (token, true)
    end
  end
  return ("", false)
end

function convert_labyrinth(arr)
  found_symbols_offset = 2
  lookup_table = Dict("" => 1)
  dims = size(arr)
  labyrinth_int = zeros(Int16, dims...)

  for x = 3:dims[1]-2
    for y = 3:dims[2]-2
      if arr[x,y] == '.'
        (token, inside) = find_token(arr, CartesianIndex(x,y))
        if token ∉ keys(lookup_table)
          lookup_table[token] = found_symbols_offset
          found_symbols_offset += 1
        end
        if ! inside
          labyrinth_int[x,y] = lookup_table[token]
        else
          labyrinth_int[x,y] = lookup_table[token] + 100
        end

      end
    end
  end
  return (labyrinth=labyrinth_int, tokens=lookup_table)
end

function create_edges(input_data, tok)
  reverse_tok= Dict(value => key for (key, value) in tok)
  mapping = Dict()
  tokens = unique(input_data)
  println(tokens)

  for i in tokens
    if i ∉ [0,1]
      start_pos = findfirst(x -> x == i, input_data)
      distances = AoC2019.bfs_depth(input_data, start_pos, x -> x != 0)
      portal_map = map(x -> x ∉ [0,1,i], input_data)
      distances .*= portal_map
      target_positions = findall(x -> x != 0, distances)
      for target_pos in target_positions
        mapping[(i, input_data[target_pos])] = distances[target_pos]
      end
      #println(reverse_tok[rem(i,100)])
      if reverse_tok[rem(i, 100)] ∉ ["AA", "ZZ"]
        other_token = (i < 100) ? i + 100 : i - 100
        #println(other_token)
        mapping[(i, other_token)] = 1
        mapping[(other_token, i)] = 1
      end
      # find all paths
    end
  end
  return mapping
end

function debug_labyrinth(arr, status, mapping)
  for current_counter = 1:maximum(status)
    pos = findfirst(x-> x == current_counter, status)
    if arr[pos] ∉ [0,1]
      suffix = (arr[pos] < 100) ? "_i" : "_o"
      println("Portal: $(mapping[rem(arr[pos], 100)])$(suffix)" )
    end
  end
end
function workflow(text)
  labyrinth_arr = AoC2019.text2arr(text)
  (lab, tok) = convert_labyrinth(labyrinth_arr)
  println(lab)
  println(tok)
  println(unique(lab))
  status = bfs_depth_warping(lab, findfirst(x -> x == tok["AA"], lab), x-> x == (tok["ZZ"]))
  reverse_mapping = Dict(value => key for (key, value) in tok)
  debug_labyrinth(lab, status, reverse_mapping)
  return maximum(status)
end

function workflow_dijkstra(text)
  labyrinth_arr = AoC2019.text2arr(text)
  (lab, tok) = convert_labyrinth(labyrinth_arr)
  distances = create_edges(lab, tok)
  println(distances)
  shortest_distances = AoC2019.dijkstra(distances, tok["AA"])
  return shortest_distances[tok["ZZ"]]
end

function workflow_recursive_dijkstra(text, do_recursive=true)
  labyrinth_arr = AoC2019.text2arr(text)
  (lab, tok) = convert_labyrinth(labyrinth_arr)
  println(tok)
  reverse_tok = Dict(value => key for (key, value) in tok)

  distances = create_edges(lab, tok)
  mg = AoC2019.generate_distance_graph(distances)

  if do_recursive
    lab_outer = deepcopy(lab)
    lab_inner = deepcopy(lab)

    lab_outer[findall(x -> x ∉ [tok["AA"], tok["ZZ"], 0, 1] && x < 100, lab_outer)] .= 0
    lab_inner[findall(x -> x in [tok["AA"], tok["ZZ"]], lab_inner)] .= 0

    distances_outer = create_edges(lab_outer, tok)
    #println(distances_outer)
    distances_inner = create_edges(lab_inner, tok)
    println("creating outer edges")
    mg_outer = AoC2019.generate_distance_graph(distances_outer)
    println(mg_outer)
    println(mg_outer.weights[113,13])
    println("creating inner edges")
    mg_inner = AoC2019.generate_distance_graph(distances_inner)
    println(mg_inner)
  end

  ind_beg = tok["AA"]
  ind_end = tok["ZZ"]

  heap = []
  push!(heap, (ind_beg, 0, 0, (ind_beg, 0, 0)))
  predecessors = Dict()

  while length(heap) > 0
    (index, num_steps, depth, pred) = popfirst!(heap)
    if (index, depth) in keys(predecessors)
      continue
    end
    predecessors[(index, depth)] = pred
    #println("$(reverse_tok[rem(index, 100)]), $num_steps, $depth")
    if index == ind_end
      return (steps=num_steps, pred=predecessors, tok=tok, mg_outer=mg_outer, mg_inner=mg_inner)
    end
    if do_recursive
      if depth > 0
        mg = mg_inner
      else
        mg = mg_outer
      end
    end
    for neighbor in outneighbors(mg, index)
      new_num_steps = num_steps + mg.weights[index, neighbor]
      if mg.weights[index, neighbor] == 0
        #error("$index, $neighbor has weight 0")
      end
      #if neighbor == 13
      #  println(mg.weights[index, neighbor])
      #  println(depth)
      #  break
      #end
      dist = neighbor - index
      if dist == 100
        new_depth = depth - 1
      elseif dist == -100
        new_depth = depth + 1
      else
        new_depth = depth
      end
      if do_recursive
        push!(heap, (neighbor, new_num_steps, new_depth, (index, depth, Int(mg.weights[index, neighbor]))))
      else
        push!(heap, (neighbor, new_num_steps, depth, (index, depth, Int(mg.weights[index, neighbor]))))
      end
    end
    sort!(heap, by=x->x[2])

  end
  error("Path not found")
end

function print_predecessors(predecessors, start, tok)
  #println(predecessors)
  reverse_tok = Dict(value => key for (key, value) in tok)
  current_node = (start, 0)
  while true
    #println(current_node)
    if reverse_tok[rem(current_node[1], 100)] == "AA"
      break
    end
    current_node = predecessors[current_node[1], current_node[2]]
    nodename = reverse_tok[rem(current_node[1], 100)]
    suffix = (current_node[1] < 100) ? "_o" : "_i"
    println("node: $(nodename)$(suffix), depth=$(current_node[2]), steps=$(current_node[3])")
  end
end

#@test workflow_dijkstra(test_input) == 23
#@test workflow_recursive_dijkstra(test_input, false).steps == 23
##@test workflow(test_input) == 23
#@test workflow_recursive_dijkstra(test_input2, false).steps == 58
#@test workflow_recursive_dijkstra(test_input, true).steps == 26
#@test workflow_recursive_dijkstra(test_input3, true).steps == 396
results = workflow_recursive_dijkstra(test_input3, true)
#print_predecessors(results.pred, results.tok["ZZ"], results.tok)

#results = workflow_recursive_dijkstra(input_data, true)
print_predecessors(results.pred, results.tok["ZZ"], results.tok)
#status = workflow(input_data)

max_stars_repo_name: buotex/aoc

max_issues_repo_name: buotex/aoc

max_forks_repo_name: buotex/aoc


text:
#-----------------------------------------------------------------------# Counter
"""
    Counter(T=Number)

Count the number of items in a data stream with elements of type `T`.

# Example

    fit!(Counter(Int), 1:100)
"""
mutable struct Counter{T} <: OnlineStat{T}
    n::Int
    Counter{T}() where {T} = new{T}(0)
end
Counter(T = Number) = Counter{T}()
_fit!(o::Counter{T}, y) where {T} = (o.n += 1)
_merge!(a::Counter, b::Counter) = (a.n += b.n)

#-----------------------------------------------------------------------# CountMap
"""
    CountMap(T::Type)
    CountMap(dict::AbstractDict{T, Int})

Track a dictionary that maps unique values to its number of occurrences.  Similar to
`StatsBase.countmap`.

# Example

    o = fit!(CountMap(Int), rand(1:10, 1000))
    value(o)
    probs(o)
    OnlineStats.pdf(o, 1)
    collect(keys(o))
"""
mutable struct CountMap{T, A <: AbstractDict{T, Int}} <: OnlineStat{T}
    value::A  # OrderedDict by default
    n::Int
end
CountMap{T}() where {T} = CountMap{T, OrderedDict{T,Int}}(OrderedDict{T,Int}(), 0)
CountMap(T::Type = Any) = CountMap{T, OrderedDict{T,Int}}(OrderedDict{T, Int}(), 0)
CountMap(d::D) where {T,D<:AbstractDict{T, Int}} = CountMap{T, D}(d, 0)
function _fit!(o::CountMap, x)
    o.n += 1
    o.value[x] = get!(o.value, x, 0) + 1
end
_merge!(o::CountMap, o2::CountMap) = (merge!(+, o.value, o2.value); o.n += o2.n)
function probs(o::CountMap, kys = keys(o.value))
    out = zeros(Int, length(kys))
    valkeys = keys(o.value)
    for (i, k) in enumerate(kys)
        out[i] = k in valkeys ? o.value[k] : 0
    end
    sum(out) == 0 ? Float64.(out) : out ./ sum(out)
end
pdf(o::CountMap, y) = y in keys(o.value) ? o.value[y] / nobs(o) : 0.0
Base.keys(o::CountMap) = keys(o.value)
nkeys(o::CountMap) = length(o.value)
Base.values(o::CountMap) = values(o.value)
Base.getindex(o::CountMap, i) = o.value[i]

#-----------------------------------------------------------------------------# CountMissing
"""
    CountMissing(stat)

Calculate a `stat` along with the count of `missing` values.

# Example

    o = CountMissing(Mean())
    fit!(o, [1, missing, 3])
"""
mutable struct CountMissing{T, O<:OnlineStat{T}} <: OnlineStat{Union{Missing,T}}
    stat::O
    nmissing::Int
end
CountMissing(stat::OnlineStat) = CountMissing(stat, 0)
value(o::CountMissing) = (nmissing=o.nmissing, stat=o.stat)
nobs(o::CountMissing) = nobs(o.stat) + o.nmissing

_fit!(o::CountMissing, x) = _fit!(o.stat, x)
_fit!(o::CountMissing, ::Missing) = (o.nmissing += 1)

_merge!(a::CountMissing, b::CountMissing) = (merge!(a.stat, b.stat); a.nmissing += b.nmissing)

#-----------------------------------------------------------------------# CovMatrix
"""
    CovMatrix(p=0; weight=EqualWeight())
    CovMatrix(::Type{T}, p=0; weight=EqualWeight())

Calculate a covariance/correlation matrix of `p` variables.  If the number of variables is
unknown, leave the default `p=0`.

# Example

    o = fit!(CovMatrix(), randn(100, 4) |> eachrow)
    cor(o)
    cov(o)
    mean(o)
    var(o)
"""
mutable struct CovMatrix{T,W} <: OnlineStat{Union{Tuple, NamedTuple, AbstractVector}} where T<:Number
    value::Matrix{T}
    A::Matrix{T}  # x'x/n
    b::Vector{T}  # 1'x/n
    weight::W
    n::Int
end
function CovMatrix(::Type{T}, p::Int=0; weight = EqualWeight()) where T<:Number
    CovMatrix(zeros(T,p,p), zeros(T,p,p), zeros(T,p), weight, 0)
end
CovMatrix(p::Int=0; weight = EqualWeight()) = CovMatrix(zeros(p,p), zeros(p,p), zeros(p), weight, 0)
function _fit!(o::CovMatrix{T}, x) where {T}
    γ = o.weight(o.n += 1)
    if isempty(o.A)
        p = length(x)
        o.b = zeros(T, p)
        o.A = zeros(T, p, p)
        o.value = zeros(T, p, p)
    end
    smooth!(o.b, x, γ)
    smooth_syr!(o.A, x, γ)
end
nvars(o::CovMatrix) = size(o.A, 1)
function value(o::CovMatrix; corrected::Bool = true)
    o.value[:] = Matrix(Hermitian((o.A - o.b * o.b')))
    corrected && rmul!(o.value, bessel(o))
    o.value
end
function _merge!(o::CovMatrix, o2::CovMatrix)
    γ = o2.n / (o.n += o2.n)
    smooth!(o.A, o2.A, γ)
    smooth!(o.b, o2.b, γ)
end
Statistics.cov(o::CovMatrix; corrected::Bool = true) = value(o; corrected=corrected)
Statistics.mean(o::CovMatrix) = o.b
Statistics.var(o::CovMatrix; kw...) = diag(value(o; kw...))
function Statistics.cor(o::CovMatrix; kw...)
    value(o; kw...)
    v = 1.0 ./ sqrt.(diag(o.value))
    rmul!(o.value, Diagonal(v))
    lmul!(Diagonal(v), o.value)
    o.value
end

#-----------------------------------------------------------------------# Extrema
"""
    Extrema(T::Type = Float64)

Maximum and minimum (and number of occurrences for each) for a data stream of type `T`.

# Example

    o = fit!(Extrema(), rand(10^5))
    extrema(o)
    maximum(o)
    minimum(o)
"""
# T is type to store data, S is type of single observation.
# E.g. you may want to accept any Number even if you are storing values as Float64
mutable struct Extrema{T,S} <: OnlineStat{S}
    min::T
    max::T
    nmin::Int
    nmax::Int
    n::Int
end
function Extrema(T::Type = Float64)
    a, b, S = extrema_init(T)
    Extrema{T,S}(a, b, 0, 0, 0)
end
extrema_init(T::Type{<:Number}) = typemax(T), typemin(T), Number
extrema_init(T::Type{String}) = "", "", String
extrema_init(T::Type{Date}) = typemax(Date), typemin(Date), Union{Date, Dates.AbstractDateTime}
extrema_init(T::Type) = rand(T), rand(T), T
function _fit!(o::Extrema, y)
    (o.n += 1) == 1 && (o.min = o.max = y)
    if y < o.min
        o.min = y
        o.nmin = 0
    elseif y > o.max
        o.max = y
        o.nmax = 0
    end
    y == o.min && (o.nmin += 1)
    y == o.max && (o.nmax += 1)
end
function _merge!(a::Extrema, b::Extrema)
    if a.min == b.min
        a.nmin += b.nmin
    elseif b.min < a.min
        a.min = b.min
        a.nmin = b.nmin
    end
    if a.max == b.max
        a.nmax += b.nmax
    elseif b.max > a.max
        a.max = b.max
        a.nmax = b.nmax
    end
    a.n += b.n
end
value(o::Extrema) = (min=o.min, max=o.max, nmin=o.nmin, nmax=o.nmax)
Base.extrema(o::Extrema) = (o.min, o.max)
Base.maximum(o::Extrema) = o.max
Base.minimum(o::Extrema) = o.min

#-----------------------------------------------------------------------# Group
"""
    Group(stats::OnlineStat...)
    Group(; stats...)
    Group(collection)

Create a vector-input stat from several scalar-input stats.  For a new
observation `y`, `y[i]` is sent to `stats[i]`.

# Examples

    x = randn(100, 2)

    fit!(Group(Mean(), Mean()), eachrow(x))
    fit!(Group(Mean(), Variance()), eachrow(x))

    o = fit!(Group(m1 = Mean(), m2 = Mean()), eachrow(x))
    o.stats.m1
    o.stats.m2
"""
struct Group{T, S} <: StatCollection{S}
    stats::T
    function Group(stats::T) where {T}
        inputs = map(input, stats)
        tup = Tuple{inputs...}
        S = Union{tup, NamedTuple{names, R} where R<:tup, AbstractVector{<: promote_type(inputs...)}} where names
        new{T,S}(stats)
    end
end
Group(o::OnlineStat...) = Group(o)
Group(;o...) = Group(o.data)
nobs(o::Group) = nobs(first(o.stats))
Base.:(==)(a::Group, b::Group) = all(a.stats .== b.stats)

Base.getindex(o::Group, i) = o.stats[i]
Base.first(o::Group) = first(o.stats)
Base.last(o::Group) = last(o.stats)
Base.lastindex(o::Group) = length(o)
Base.length(o::Group) = length(o.stats)
Base.values(o::Group) = map(value, o.stats)

Base.iterate(o::Group) = (o.stats[1], 2)
Base.iterate(o::Group, i) = i > length(o) ? nothing : (o.stats[i], i + 1)

@generated function _fit!(o::Group{T}, y) where {T}
    N = fieldcount(T)
    :(Base.Cartesian.@nexprs $N i -> @inbounds(_fit!(o.stats[i], y[i])))
end
function _fit!(o::Group{T}, y) where {T<:AbstractVector}
    for (i,yi) in enumerate(y)
        _fit!(o.stats[i], yi)
    end
end

_merge!(o::Group, o2::Group) = map(merge!, o.stats, o2.stats)

Base.:*(n::Integer, o::OnlineStat) = Group([copy(o) for i in 1:n]...)

#-----------------------------------------------------------------------# GroupBy
"""
    GroupBy(T, stat)

Update `stat` for each group (of type `T`).  A single observation is either a (named)
tuple with two elements or a Pair.

# Example

    x = rand(Bool, 10^5)
    y = x .+ randn(10^5)
    fit!(GroupBy(Bool, Series(Mean(), Extrema())), zip(x,y))
"""
mutable struct GroupBy{T, S, O <: OnlineStat{S}} <: StatCollection{TwoThings{T,S}}
    value::OrderedDict{T, O}
    init::O
    n::Int
    function GroupBy(value::OrderedDict{T,O}, init::O, n::Int) where {T,S,O<:OnlineStat{S}}
        new{T,S,O}(value, init, n)
    end
end
GroupBy(T::Type, stat::O) where {O<:OnlineStat} = GroupBy(OrderedDict{T, O}(), stat, 0)
function _fit!(o::GroupBy, xy)
    o.n += 1
    x, y = xy
    x in keys(o.value) ? fit!(o.value[x], y) : (o.value[x] = fit!(copy(o.init), y))
end
Base.getindex(o::GroupBy{T}, i::T) where {T} = o.value[i]
AbstractTrees.children(o::GroupBy) = collect(o.value)
AbstractTrees.printnode(io::IO, o::GroupBy{T,S,O}) where {T,S,O} = print(io, "GroupBy: $T => $(name(O,false,false))")
Base.sort!(o::GroupBy) = (sort!(o.value); o)
function _merge!(a::GroupBy{T,O}, b::GroupBy{T,O}) where {T,O}
    a.init == b.init || error("Cannot merge GroupBy objects with different inits")
    a.n += b.n
    merge!((o1, o2) -> merge!(o1, o2), a.value, b.value)
end

#-----------------------------------------------------------------------# Mean
"""
    Mean(T = Float64; weight=EqualWeight())

Track a univariate mean, stored as type `T`.

# Example

    @time fit!(Mean(), randn(10^6))
"""
mutable struct Mean{T,W} <: OnlineStat{Number}
    μ::T
    weight::W
    n::Int
end
Mean(T::Type{<:Number} = Float64; weight = EqualWeight()) = Mean(zero(T), weight, 0)
_fit!(o::Mean{T}, x) where {T} = (o.μ = smooth(o.μ, x, T(o.weight(o.n += 1))))
function _merge!(o::Mean, o2::Mean)
    o.n += o2.n
    o.μ = smooth(o.μ, o2.μ, o2.n / o.n)
end
Statistics.mean(o::Mean) = o.μ
Base.copy(o::Mean) = Mean(o.μ, o.weight, o.n)

#-----------------------------------------------------------------------# Moments
"""
    Moments(; weight=EqualWeight())

First four non-central moments.

# Example

    o = fit!(Moments(), randn(1000))
    mean(o)
    var(o)
    std(o)
    skewness(o)
    kurtosis(o)
"""
mutable struct Moments{W} <: OnlineStat{Number}
    m::Vector{Float64}
    weight::W
    n::Int
end
Moments(;weight = EqualWeight()) = Moments(zeros(4), weight, 0)
function _fit!(o::Moments, y::Real)
    γ = o.weight(o.n += 1)
    y2 = y * y
    @inbounds o.m[1] = smooth(o.m[1], y, γ)
    @inbounds o.m[2] = smooth(o.m[2], y2, γ)
    @inbounds o.m[3] = smooth(o.m[3], y * y2, γ)
    @inbounds o.m[4] = smooth(o.m[4], y2 * y2, γ)
end
Statistics.mean(o::Moments) = o.m[1]
function Statistics.var(o::Moments; corrected=true)
    out = (o.m[2] - o.m[1] ^ 2)
    corrected ? bessel(o) * out : out
end
function StatsBase.skewness(o::Moments)
    v = value(o)
    vr = o.m[2] - o.m[1]^2
    (v[3] - 3.0 * v[1] * vr - v[1] ^ 3) / vr ^ 1.5
end
function StatsBase.kurtosis(o::Moments)
    # v = value(o)
    # (v[4] - 4.0 * v[1] * v[3] + 6.0 * v[1] ^ 2 * v[2] - 3.0 * v[1] ^ 4) / var(o) ^ 2 - 3.0
    m1, m2, m3, m4 = value(o)
    (m4 - 4.0 * m1 * m3 + 6.0 * m1^2 * m2 - 3.0 * m1 ^ 4) / var(o; corrected=false) ^ 2 - 3.0
end
function _merge!(o::Moments, o2::Moments)
    γ = o2.n / (o.n += o2.n)
    smooth!(o.m, o2.m, γ)
end

#-----------------------------------------------------------------------# Sum
"""
    Sum(T::Type = Float64)

Track the overall sum.

# Example

    fit!(Sum(Int), fill(1, 100))
"""
mutable struct Sum{T} <: OnlineStat{Number}
    sum::T
    n::Int
end
Sum(T::Type = Float64) = Sum(T(0), 0)
Base.sum(o::Sum) = o.sum
_fit!(o::Sum{T}, x::Real) where {T<:AbstractFloat} = (o.sum += convert(T, x); o.n += 1)
_fit!(o::Sum{T}, x::Real) where {T<:Integer} =       (o.sum += round(T, x); o.n += 1)
_merge!(o::T, o2::T) where {T <: Sum} = (o.sum += o2.sum; o.n += o2.n; o)

#-----------------------------------------------------------------------# Variance
"""
    Variance(T = Float64; weight=EqualWeight())

Univariate variance, tracked as type `T`.

# Example

    o = fit!(Variance(), randn(10^6))
    mean(o)
    var(o)
    std(o)
"""
mutable struct Variance{T, W} <: OnlineStat{Number}
    σ2::T
    μ::T
    weight::W
    n::Int
end
function Variance(T::Type{<:Number} = Float64; weight = EqualWeight())
    Variance(zero(T), zero(T), weight, 0)
end
Base.copy(o::Variance) = Variance(o.σ2, o.μ, o.weight, o.n)
function _fit!(o::Variance{T}, x) where {T}
    μ = o.μ
    γ = T(o.weight(o.n += 1))
    o.μ = smooth(o.μ, T(x), γ)
    o.σ2 = smooth(o.σ2, (T(x) - o.μ) * (T(x) - μ), γ)
end
function _merge!(o::Variance, o2::Variance)
    γ = o2.n / (o.n += o2.n)
    δ = o2.μ - o.μ
    o.σ2 = smooth(o.σ2, o2.σ2, γ) + δ ^ 2 * γ * (1.0 - γ)
    o.μ = smooth(o.μ, o2.μ, γ)
end
value(o::Variance) = o.n > 1 ? o.σ2 * bessel(o) : 1.0
Statistics.var(o::Variance) = value(o)
Statistics.mean(o::Variance) = o.μ

#-----------------------------------------------------------------------# Series
"""
    Series(stats)
    Series(stats...)
    Series(; stats...)

Track a collection stats for one data stream.

# Example

    s = Series(Mean(), Variance())
    fit!(s, randn(1000))
"""
struct Series{IN, T} <: StatCollection{IN}
    stats::T
    Series(stats::T) where {T} = new{Union{map(input, stats)...}, T}(stats)
end
Series(t::OnlineStat...) = Series(t)
Series(; t...) = Series(t.data)

value(o::Series) = map(value, o.stats)
nobs(o::Series) = nobs(o.stats[1])
@generated function _fit!(o::Series{IN, T}, y) where {IN, T}
    n = length(fieldnames(T))
    :(Base.Cartesian.@nexprs $n i -> _fit!(o.stats[i], y))
end
_merge!(o::Series, o2::Series) = map(_merge!, o.stats, o2.stats)

#-----------------------------------------------------------------------# FTSeries
"""
    FTSeries(stats...; filter=x->true, transform=identity)

Track multiple stats for one data stream that is filtered and transformed before being
fitted.

    FTSeries(T, stats...; filter, transform)

Create an FTSeries and specify the type `T` of the pre-transformed values.

# Example

    o = FTSeries(Mean(), Variance(); transform=abs)
    fit!(o, -rand(1000))

    # Remove missing values represented as DataValues
    using DataValues
    y = DataValueArray(randn(100), rand(Bool, 100))
    o = FTSeries(DataValue, Mean(); transform=get, filter=!isna)
    fit!(o, y)

    # Remove missing values represented as Missing
    y = [rand(Bool) ? rand() : missing for i in 1:100]
    o = FTSeries(Mean(); filter=!ismissing)
    fit!(o, y)

    # Alternatively for Missing:
    fit!(Mean(), skipmissing(y))
"""
mutable struct FTSeries{IN, OS, F, T} <: StatCollection{Union{IN,Missing}}
    stats::OS
    filter::F
    transform::T
    nfiltered::Int
end
function FTSeries(stats::OnlineStat...; filter=x->true, transform=identity)
    IN, OS = Union{map(input, stats)...}, typeof(stats)
    FTSeries{IN, OS, typeof(filter), typeof(transform)}(stats, filter, transform, 0)
end
function FTSeries(T::Type, stats::OnlineStat...; filter=x->true, transform=identity)
    FTSeries{T, typeof(stats), typeof(filter), typeof(transform)}(stats, filter, transform, 0)
end
value(o::FTSeries) = value.(o.stats)
nobs(o::FTSeries) = nobs(o.stats[1])
@generated function _fit!(o::FTSeries{N, OS}, y) where {N, OS}
    n = length(fieldnames(OS))
    quote
        if o.filter(y)
            yt = o.transform(y)
            Base.Cartesian.@nexprs $n i -> @inbounds begin
                _fit!(o.stats[i], yt)
            end
        else
            o.nfiltered += 1
        end
    end
end
function _merge!(o::FTSeries, o2::FTSeries)
    o.nfiltered += o2.nfiltered
    _merge!.(o.stats, o2.stats)
end


text:
#contains functions related to building, loading, and saving metadata
#currently metadata is a vector of MetaDictionary, but it could eventually
#become its own object

#reads in the metadata and creates type tuples
function readMetaIndex(;indexPath::String = DOCUMENTATION_PATH,
  fileName::String = META_INDEX_NAME)::DataFrame

  iStream::IOStream = open("$indexPath\\$fileName.csv")
  df::DataFrame = CSV.read(iStream, strings=:raw, categorical=false)
  close(iStream)

  return df

end

#froms the vector of Metadata Dictionaries
function formMetadata(indexPath::String = DOCUMENTATION_PATH,
  indexName::String = META_INDEX_NAME)::Vector{MetaDictionary}

  #first get a list of the paths
  indexDF::DataFrame = readMetaIndex()
  fileNames::Vector{String} = indexDF[:name]
  NFiles::Int = length(fileNames)

  metadata::Vector{MetaDictionary} = (MetaDictionary).(fileNames)

  return metadata
end

#checks for and removes metadata fields not actually in the file
function validateMetaDictionary!(meta::MetaDictionary, keeperList::Vector{Symbol};
  nccsPath::String = NCCS_PATH,  streamer::Function = I_CSV_STREAMER)::MetaDictionary

  #these are the column names that we expect
  expectedCols::Vector{Symbol} = collect(keys(meta.dictionary))

  path::String = "$nccsPath\\$(meta.fileName).csv.gz"
  iStream = streamer(path)
  df::DataFrame = CSV.read(iStream, strings=:raw, categorical=false, rows=2)
  close(iStream)

  #check first that the column names are correct
  actualCols::Vector{Symbol} = (s::Symbol->Symbol(lowercase(string(s)))).(names(df))

  #get the list of files we don't want
  expectedNotFound::Vector{Symbol} = setdiff(expectedCols, actualCols)
  expectedNotWanted::Vector{Symbol} = setdiff(expectedCols, keeperList)
  fieldsToRemove::Vector{Symbol} = unique([expectedNotFound; expectedNotWanted])

  #remove the offending fields
  if length(fieldsToRemove) > 0

    removeFields!(meta, fieldsToRemove)
    println("Validated $(meta.fileName). Removed $(length(fieldsToRemove)) fields from metadata.")
  else
    println("Validated $(meta.fileName).")
  end

  return meta

end

function validateMetadata(metadata::Vector{MetaDictionary};
    indexPath::String = DOCUMENTATION_PATH, keeperName::String=KEEPER_NAME)::Vector{MetaDictionary}

  metadata = deepcopy(metadata) #its hard to do this without some modifications AND copies
                                #hence err on the side of caution here

  #metadata = removeDictionary(metadata, "nccs_core_2013_co_legacy")
  NMeta::Int = length(metadata)

  #get list of kept columns
  iStream::IOStream = open("$indexPath\\$keeperName.csv")
  keeperDF::DataFrame = CSV.read(iStream, strings=:raw, categorical=false)
  close(iStream)
  keeperList::Vector{Symbol} = (Symbol).(keeperDF[:keeper])

  #this can take a while so multi-thread it
  println("Beginning file validation")

  futures::Vector{Future} = Vector{Future}()

  #use this to check orders
  metaNames::Vector{String} = (m::MetaDictionary->m.fileName).(metadata)

  #assign and retrieve the work
  for i::Int ∈ 1:NMeta
    push!(futures, (@spawn validateMetaDictionary!(metadata[i], keeperList)))
  end
  metadata .= (fetch).(futures)

  return metadata
end

#this function either forms or loads the metadata
function getMetadata(;refreshMetadata::Bool = true, metadataPath::String = DOCUMENTATION_PATH,
  metadataName::String = METADATA_NAME)

  metadataFilePath::String = "$metadataPath\\$metadataName.jls"

  #load the metadata if we don't need to redo the validation
  if refreshMetadata || !(isfile(metadataFilePath))
    metadata::Vector{MetaDictionary}= formMetadata()
    metadata = validateMetadata(metadata)

    oStream::IOStream = open(metadataFilePath, "w")
    serialize(oStream, metadata)
    close(oStream)
  else
    iStream::IOStream = open(metadataFilePath)
    metadata = deserialize(iStream)
    close(iStream)
  end

  return metadata
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA