text:
# API (Two-stage) #
# ========================== #
"""
    instantiate(stochasticmodel::StochasticModel{2},
                scenarios::Vector{<:AbstractScenario};
                solver = JuMP.UnsetSolver(),
                procs = workers(),
                kw...)

Instantiate a new two-stage stochastic program using the model definition stored in the two-stage `stochasticmodel`, and the given collection of `scenarios`.
"""
function instantiate(sm::StochasticModel{2}, scenarios::Vector{<:AbstractScenario}; solver = JuMP.UnsetSolver(), procs = workers(), kw...)
    sp = StochasticProgram(parameters(sm.parameters[1]; kw...),
                           parameters(sm.parameters[2]; kw...),
                           scenarios,
                           solver,
                           procs)
    sm.generator(sp)
    return sp
end
"""
    instantiate(stochasticmodel::StochasticModel{2};
                scenariotype::Type{S} = Scenario,
                solver = JuMP.UnsetSolver(),
                procs = workers(),
                kw...) where S <: AbstractScenario

Instantiate a deferred two-stage stochastic program using the model definition stored in the two-stage `stochasticmodel` over the scenario type `S`.
"""
function instantiate(sm::StochasticModel{2}; scenariotype::Type{S} = Scenario, solver = JuMP.UnsetSolver(), procs = workers(), kw...) where S <: AbstractScenario
    sp = StochasticProgram(parameters(sm.parameters[1]; kw...),
                           parameters(sm.parameters[2]; kw...),
                           scenariotype,
                           solver,
                           procs)
    sm.generator(sp)
    return sp
end
"""
    instantiate(stochasticmodel::StochasticModel,
                scenarios::Vector{<:AbstractScenario};
                solver = JuMP.UnsetSolver(),
                procs = workers(),
                kw...)

Instantiate a new stochastic program using the model definition stored in `stochasticmodel`, and the given collection of `scenarios`.
"""
function instantiate(sm::StochasticModel{N},
                     scenarios::NTuple{M,Vector{<:AbstractScenario}};
                     solver = JuMP.UnsetSolver(),
                     procs = workers(),
                     kw...) where {N,M}
    M == N - 1 || error("Inconsistent number of stages $N and number of scenario types $M")
    params = ntuple(Val(N)) do i
        parameters(sm.parameters[i]; kw...)
    end
    sp = StochasticProgram(params,
                           scenarios,
                           solver,
                           procs)
    sm.generator(sp)
    return sp
end
"""
    optimize!(stochasticprogram::StochasticProgram; solver::SPSolverType = JuMP.UnsetSolver())

Optimize the `stochasticprogram` in expectation using `solver`.

## Examples

The following solves the stochastic program `sp` using the L-shaped algorithm.

```julia
using LShapedSolvers
using GLPKMathProgInterface

optimize!(sp, solver = LShapedSolver(:ls, GLPKSolverLP()));

# output

L-Shaped Gap  Time: 0:00:01 (4 iterations)
  Objective:       -855.8333333333339
  Gap:             0.0
  Number of cuts:  7
:Optimal
```

The following solves the stochastic program `sp` using GLPK on the extended form.

```julia
using GLPKMathProgInterface

optimize!(sp, solver = GLPKSolverLP())

:Optimal
```

See also: [`VRP`](@ref)
"""
function optimize!(stochasticprogram::StochasticProgram{2}; solver::SPSolverType = JuMP.UnsetSolver(), kwargs...)
    if deferred(stochasticprogram)
        generate!(stochasticprogram)
    end
    # Prefer cached solver if available
    supplied_solver = pick_solver(stochasticprogram, solver)
    # Abort if no solver was given
    if isa(supplied_solver, JuMP.UnsetSolver)
        error("Cannot optimize without a solver.")
    end
    # Switch on solver type
    return _optimize!(stochasticprogram, supplied_solver; kwargs...)
end
function _optimize!(stochasticprogram::StochasticProgram{2}, solver::MathProgBase.AbstractMathProgSolver; kwargs...)
    # Standard mathprogbase solver. Fallback to solving DEP, relying on JuMP.
    dep = DEP(stochasticprogram; solver = solver)
    status = solve(dep; kwargs...)
    stochasticprogram.spsolver.internal_model = dep.internalModel
    fill_solution!(stochasticprogram)
    return status
end
function _optimize!(stochasticprogram::StochasticProgram{2}, solver::AbstractStructuredSolver; kwargs...)
    # Use structured solver
    structuredmodel = StructuredModel(stochasticprogram, solver)
    stochasticprogram.spsolver.internal_model = structuredmodel
    status = optimize_structured!(structuredmodel)
    fill_solution!(stochasticprogram, structuredmodel)
    # Now safe to generate the objective value of the stochastic program
    calculate_objective_value!(stochasticprogram)
    return status
end
"""
    optimize(stochasticmodel::StochasticModel, sampler::AbstractSampler; solver::SPSolverType = JuMP.UnsetSolver(), confidence::AbstractFloat = 0.95, kwargs...)

Approximately optimize the `stochasticmodel` using `solver` when the underlying scenario distribution is inferred by `sampler` and return a `StochasticSolution` with the given `confidence` level.

See also: [`StochasticSolution`](@ref)
"""
function optimize(stochasticmodel::StochasticModel, sampler::AbstractSampler; solver::SPSolverType = JuMP.UnsetSolver(), confidence::AbstractFloat = 0.95, kwargs...)
    # Abort if no solver was given
    if isa(solver, JuMP.UnsetSolver)
        error("Cannot optimize without a solver.")
    end
    # Switch on solver type
    return _optimize(stochasticmodel, sampler, solver, confidence; kwargs...)
end
function _optimize(stochasticmodel::StochasticModel, sampler::AbstractSampler, solver::MathProgBase.AbstractMathProgSolver, confidence::AbstractFloat; kwargs...)
    return _optimize(stochasticmodel, sampler, SAASolver(solver), confidence; kwargs...)
end
function _optimize(stochasticmodel::StochasticModel, sampler::AbstractSampler, solver::AbstractStructuredSolver, confidence::AbstractFloat; kwargs...)
    return _optimize(stochasticmodel, sampler, SAASolver(solver), confidence; kwargs...)
end
function _optimize(stochasticmodel::StochasticModel, sampler::AbstractSampler, solver::AbstractSampledSolver, confidence::AbstractFloat; kwargs...)
    sampledmodel = SampledModel(stochasticmodel, solver)
    status = optimize_sampled!(sampledmodel, sampler, confidence; kwargs...)
    if status != :Optimal
        @warn "Optimal solution not found. Returned status $status"
    end
    return stochastic_solution(sampledmodel)
end
"""
    optimal_decision(stochasticprogram::StochasticProgram)

Return the optimal first stage decision of `stochasticprogram`, after a call to `optimize!(stochasticprogram)`.
"""
function optimal_decision(stochasticprogram::StochasticProgram)
    decision = get_stage_one(stochasticprogram).colVal
    if any(isnan.(decision))
        @warn "Optimal decision not defined. Check that the model was properly solved."
    end
    return decision
end
"""
    optimal_decision(stochasticprogram::StochasticProgram, var::Symbol)

Return the optimal first stage variable `var` of `stochasticprogram`, after a call to `optimize!(stochasticprogram)`.
"""
function optimal_decision(stochasticprogram::StochasticProgram, var::Symbol)
    return getvalue(get_stage_one(stochasticprogram).objDict[var])
end
"""
    optimal_decision(stochasticprogram::StochasticProgram, i::Integer)

Return the optimal second stage decision of `stochasticprogram` in the `i`th scenario, after a call to `optimize!(stochasticprogram)`.
"""
function optimal_decision(stochasticprogram::StochasticProgram, i::Integer)
    submodel = subproblem(stochasticprogram, i)
    decision = submodel.colVal
    if any(isnan.(decision))
        @warn "Optimal decision not defined in subproblem $i. Check that the model was properly solved."
    end
    return decision
end
"""
    optimal_decision(stochasticprogram::StochasticProgram, var::Symbol)

Return the optimal second stage variable `var` of `stochasticprogram` in the `i`th scenario, after a call to `optimize!(stochasticprogram)`.
"""
function optimal_decision(stochasticprogram::StochasticProgram, i::Integer, var::Symbol)
    submodel = subproblem(stochasticprogram, i)
    return getvalue(submodel.objDict[var])
end
"""
    optimal_value(stochasticprogram::StochasticProgram)

Return the optimal value of `stochasticprogram`, after a call to `optimize!(stochasticprogram)`.
"""
function optimal_value(stochasticprogram::StochasticProgram)
    return get_stage_one(stochasticprogram).objVal
end
"""
    optimal_value(stochasticprogram::StochasticProgram, i::Integer)

Return the optimal value of the `i`th subproblem in `stochasticprogram`, after a call to `optimize!(stochasticprogram)`.
"""
function optimal_value(stochasticprogram::StochasticProgram, i::Integer)
    submodel = subproblem(stochasticprogram, i)
    return submodel.objVal
end
"""
    stage_parameters(stochasticprogram::StochasticProgram, s::Integer)

Return the parameters at stage `s` in `stochasticprogram`.
"""
function stage_parameters(stochasticprogram::StochasticProgram{N}, s::Integer) where N
    1 <= s <= N || error("Stage $s not in range 1 to $N.")
    return stochasticprogram.stages[s].parameters
end
"""
    scenarioproblems(stochasticprogram::StochasticProgram, s::Integer)

Return the scenario problems at stage `s` in `stochasticprogram`.
"""
function scenarioproblems(stochasticprogram::StochasticProgram{N}, s::Integer) where N
    1 < s <= N || error("Stage $s not in range 2 to $N.")
    return stochasticprogram.scenarioproblems[s-1]
end
"""
    scenarioproblems(stochasticprogram::TwoStageStochasticProgram)

Return the scenario problems in `stochasticprogram`.
"""
function scenarioproblems(stochasticprogram::StochasticProgram{2}, s::Integer = 2)
    s == 1 && error("Stage 1 does not have scenario problems.")
    s == 2 || error("Stage $s not available in two-stage model.")
    return stochasticprogram.scenarioproblems
end
"""
    decision_length(stochasticprogram::StochasticProgram, s::Integer)

Return the length of the decision at stage `s` in the `stochasticprogram`.
"""
function decision_length(stochasticprogram::StochasticProgram{N}, s::Integer) where N
    if s == 1
        haskey(stochasticprogram.problemcache, :stage_1) || return 0
        return stochasticprogram.problemcache[:stage_1].numCols
    end
    nsubproblems(stochasticprogram, s) == 0 && return 0
    return recourse_length(scenarioproblems(stochasticprogram, s))
end
"""
    decision_length(stochasticprogram::TwoStageStochasticProgram)

Return the length of the first-stage decision of the two-stage `stochasticprogram`.
"""
function decision_length(stochasticprogram::StochasticProgram{2})
    return decision_length(stochasticprogram, 1)
end
"""
    first_stage_nconstraints(stochasticprogram::StochasticProgram)

Return the number of constraints in the the first stage of `stochasticprogram`.
"""
function first_stage_nconstraints(stochasticprogram::StochasticProgram)
    !haskey(stochasticprogram.problemcache, :stage_1) && return 0
    first_stage = get_stage_one(stochasticprogram)
    return length(first_stage.linconstr)
end
"""
    first_stage_dims(stochasticprogram::StochasticProgram)

Return a the number of variables and the number of constraints in the the first stage of `stochasticprogram` as a tuple.
"""
function first_stage_dims(stochasticprogram::StochasticProgram)
    !haskey(stochasticprogram.problemcache, :stage_1) && return 0, 0
    first_stage = get_stage_one(stochasticprogram)
    return length(first_stage.linconstr), first_stage.numCols
end
"""
    recourse_length(stochasticprogram::TwoStageStochasticProgram)

Return the length of the second-stage decision in the two-stage `stochasticprogram`.
"""
function recourse_length(stochasticprogram::StochasticProgram{2})
    has_generator(stochasticprogram, :stage_2) || error("Second-stage problem not defined in stochastic program. Consider @second_stage.")
    nsubproblems(stochasticprogram) == 0 && return 0
    return recourse_length(scenarioproblems(stochasticprogram))
end
"""
    scenario(stochasticprogram::StochasticProgram, i::Integer, s::Integer = 2)

Return the `i`th scenario of stochasticprogram` at stage `s`. Defaults to the second stage.
"""
function scenario(stochasticprogram::StochasticProgram, i::Integer, s::Integer = 2)
    return scenario(scenarioproblems(stochasticprogram, s), i)
end
"""
    scenarios(stochasticprogram::StochasticProgram, s::Integer = 2)

Return an array of all scenarios of the `stochasticprogram` at stage `s`. Defaults to the second stage.
"""
function scenarios(stochasticprogram::StochasticProgram, s::Integer = 2)
    return scenarios(scenarioproblems(stochasticprogram, s))
end
"""
    expected(stochasticprogram::StochasticProgram, s::Integer = 2)

Return the exected scenario of all scenarios of the `stochasticprogram` at stage `s`. Defaults to the second stage.
"""
function expected(stochasticprogram::StochasticProgram, s::Integer = 2)
    return expected(scenarioproblems(stochasticprogram, s)).scenario
end
"""
    scenariotype(stochasticprogram::StochasticProgram, s::Integer = 2)

Return the type of the scenario structure associated with `stochasticprogram` at stage `s`. Defaults to the second stage.
"""
function scenariotype(stochasticprogram::StochasticProgram, s::Integer = 2)
    return scenariotype(scenarioproblems(stochasticprogram, s))
end
"""
    probability(stochasticprogram::StochasticProgram, i::Integer, s::Integer = 2)

Return the probability of scenario `i`th scenario in the `stochasticprogram` at stage `s` occuring. Defaults to the second stage.
"""
function probability(stochasticprogram::StochasticProgram, i::Integer, s::Integer = 2)
    return probability(scenario(stochasticprogram, s, i))
end
"""
    stage_probability(stochasticprogram::StochasticProgram, s::Integer = 2)

Return the probability of any scenario in the `stochasticprogram` at stage `s` occuring. A well defined model should return 1. Defaults to the second stage.
"""
function stage_probability(stochasticprogram::StochasticProgram, s::Integer = 2)
    return probability(scenarioproblems(stochasticprogram, s))
end
"""
    has_generator(stochasticprogram::StochasticProgram, key::Symbol)

Return true if a problem generator with `key` exists in `stochasticprogram`.
"""
function has_generator(stochasticprogram::StochasticProgram, key::Symbol)
    return haskey(stochasticprogram.generator, key)
end
"""
    generator(stochasticprogram::StochasticProgram, key::Symbol)

Return the problem generator associated with `key` in `stochasticprogram`.
"""
function generator(stochasticprogram::StochasticProgram, key::Symbol)
    return stochasticprogram.generator[key]
end
"""
    subproblem(stochasticprogram::StochasticProgram, i::Integer, s::Integer = 2)

Return the `i`th subproblem of the `stochasticprogram` at stage `s`. Defaults to the second stage.
"""
function subproblem(stochasticprogram::StochasticProgram, i::Integer, s::Integer = 2)
    return subproblem(scenarioproblems(stochasticprogram, s), i)
end
"""
    subproblems(stochasticprogram::StochasticProgram, s::Integer = 2)

Return an array of all subproblems of the `stochasticprogram` at stage `s`. Defaults to the second stage.
"""
function subproblems(stochasticprogram::StochasticProgram, s::Integer = 2)
    return subproblems(scenarioproblems(stochasticprogram, s))
end
"""
    nsubproblems(stochasticprogram::StochasticProgram, s::Integer = 2)

Return the number of subproblems in the `stochasticprogram` at stage `s`. Defaults to the second stage.
"""
function nsubproblems(stochasticprogram::StochasticProgram, s::Integer = 2)
    s == 1 && return 0
    return nsubproblems(scenarioproblems(stochasticprogram, s))
end
"""
    masterterms(stochasticprogram::StochasticProgram, i::Integer, s::Integer = 2)

Return the first stage terms appearing in scenario `i` in the `stochasticprogram` at stage `s`. Defaults to the second stage.

The master terms are given in sparse format as an array of tuples `(row,col,coeff)` which specify the occurance of master problem variables in the second stage constraints.
"""
function masterterms(stochasticprogram::StochasticProgram, i::Integer, s::Integer = 2)
    return masterterms(scenarioproblems(stochasticprogram, s), i)
end
"""
    nscenarios(stochasticprogram::StochasticProgram, s::Integer = 2)

Return the number of scenarios in the `stochasticprogram` at stage `s`. Defaults to the second stage.
"""
function nscenarios(stochasticprogram::StochasticProgram, s::Integer = 2)
    s == 1 && return 0
    return nscenarios(scenarioproblems(stochasticprogram, s))
end
"""
    nstages(stochasticprogram::StochasticProgram)

Return the number of stages in `stochasticprogram`.
"""
nstages(::StochasticProgram{N}) where N = N
"""
    distributed(stochasticprogram::StochasticProgram, s::Integer = 2)

Return true if the `stochasticprogram` is memory distributed at stage `s`. Defaults to the second stage.
"""
distributed(stochasticprogram::StochasticProgram, s::Integer = 2) = distributed(scenarioproblems(stochasticprogram, s))
"""
    deferred(stochasticprogram::StochasticProgram)

Return true if `stochasticprogram` is not fully generated.
"""
deferred(stochasticprogram::StochasticProgram{N}) where N = deferred(stochasticprogram, Val(N))
deferred(stochasticprogram::StochasticProgram, ::Val{1}) = deferred_first_stage(stochasticprogram)
function deferred(stochasticprogram::StochasticProgram, ::Val{N}) where N

    return deferred_stage(stochasticprogram, N) || deferred(stochasticprogram, Val(N-1))
end
deferred_first_stage(stochasticprogram::StochasticProgram) = has_generator(stochasticprogram, :stage_1) && !haskey(stochasticprogram.problemcache, :stage_1)
function deferred_stage(stochasticprogram::StochasticProgram{N}, s::Integer) where N
    1 <= s <= N || error("Stage $s not in range 1 to $N.")
    s == 1 && return deferred_first_stage(stochasticprogram)
    nsubproblems(stochasticprogram, s) < nscenarios(stochasticprogram, s)
end
"""
    spsolver(stochasticprogram::StochasticProgram)

Return the stochastic program solver `spsolver` of the `stochasticprogram`.
"""
function spsolver(stochasticprogram::StochasticProgram)
    return stochasticprogram.spsolver.solver
end
"""
    internal_model(stochasticprogram::StochasticProgram)

Return the internal solve model object of the `stochasticprogram`, after a call to `optimize!(stochasticprogram)`.
"""
function spsolver_model(stochasticprogram::StochasticProgram)
    return stochasticprogram.spsolver.internal_model
end
# ========================== #

# Setters
# ========================== #
"""
    set_spsolver(stochasticprogram::StochasticProgram, spsolver::Union{MathProgBase.AbstractMathProgSolver,AbstractStructuredSolver})

Store the stochastic program solver `spsolver` of the `stochasticprogram`.
"""
function set_spsolver(stochasticprogram::StochasticProgram, spsolver::Union{MathProgBase.AbstractMathProgSolver,AbstractStructuredSolver})
    stochasticprogram.spsolver.solver = spsolver
    nothing
end
"""
    add_scenario!(stochasticprogram::StochasticProgram, scenario::AbstractScenario, stage::Integer = 2; defer::Bool = false)

Store the second stage `scenario` in the `stochasticprogram` at `stage`. Defaults to the second stage.

If `defer` is true, then model creation is deferred until `generate!(stochasticprogram)` is called. If the `stochasticprogram` is distributed, the scenario will be defined on the node that currently has the fewest scenarios.
"""
function add_scenario!(stochasticprogram::StochasticProgram, scenario::AbstractScenario, stage::Integer = 2; defer::Bool = false)
    add_scenario!(scenarioproblems(stochasticprogram, stage), scenario)
    invalidate_cache!(stochasticprogram)
    if !defer
        generate!(stochasticprogram)
    end
    return stochasticprogram
end
"""
    add_worker_scenario!(stochasticprogram::StochasticProgram, scenario::AbstractScenario, w::Integer, stage::Integer = 2; defer::Bool = false)

Store the second stage `scenario` in worker node `w` of the `stochasticprogram` at `stage`. Defaults to the second stage.

If `defer` is true, then model creation is deferred until `generate!(stochasticprogram)` is called.
"""
function add_worker_scenario!(stochasticprogram::StochasticProgram, scenario::AbstractScenario, w::Integer, stage::Integer = 2; defer::Bool = false)
    add_scenario!(scenarioproblems(stochasticprogram, stage), scenario, w)
    invalidate_cache!(stochasticprogram)
    if !defer
        generate!(stochasticprogram)
    end
    return stochasticprogram
end
"""
    add_scenario!(scenariogenerator::Function, stochasticprogram::StochasticProgram, stage::Integer = 2; defer::Bool = false)

Store the second stage scenario returned by `scenariogenerator` in the second stage of the `stochasticprogram`. Defaults to the second stage.

If `defer` is true, then model creation is deferred until `generate!(stochasticprogram)` is called. If the `stochasticprogram` is distributed, the scenario will be defined on the node that currently has the fewest scenarios.
"""
function add_scenario!(scenariogenerator::Function, stochasticprogram::StochasticProgram, stage::Integer = 2; defer::Bool = false)
    add_scenario!(scenariogenerator, scenarioproblems(stochasticprogram, stage))
    invalidate_cache!(stochasticprogram)
    if !defer
        generate!(stochasticprogram)
    end
    return stochasticprogram
end
"""
    add_worker_scenario!(scenariogenerator::Function, stochasticprogram::StochasticProgram, w::Integer, stage::Integer = 2; defer::Bool = false)

Store the second stage scenario returned by `scenariogenerator` in worker node `w` of the `stochasticprogram` at `stage`. Defaults to the second stage.

If `defer` is true, then model creation is deferred until `generate!(stochasticprogram)` is called.
"""
function add_worker_scenario!(scenariogenerator::Function, stochasticprogram::StochasticProgram, w::Integer, stage::Integer = 2; defer::Bool = false)
    add_scenario!(scenariogenerator, scenarioproblems(stochasticprogram, stage), w)
    invalidate_cache!(stochasticprogram)
    if !defer
        generate!(stochasticprogram)
    end
    return stochasticprogram
end
"""
    add_scenarios!(stochasticprogram::StochasticProgram, scenarios::Vector{<:AbstractScenario}, stage::Integer = 2; defer::Bool = false)

Store the collection of second stage `scenarios` in the `stochasticprogram` at `stage`. Defaults to the second stage.

If `defer` is true, then model creation is deferred until `generate!(stochasticprogram)` is called. If the `stochasticprogram` is distributed, scenarios will be distributed evenly across workers.
"""
function add_scenarios!(stochasticprogram::StochasticProgram, scenarios::Vector{<:AbstractScenario}, stage::Integer = 2; defer::Bool = false)
    add_scenarios!(scenarioproblems(stochasticprogram, stage), scenarios)
    invalidate_cache!(stochasticprogram)
    if !defer
        generate!(stochasticprogram)
    end
    return stochasticprogram
end
"""
    add_worker_scenarios!(stochasticprogram::StochasticProgram, scenarios::Vector{<:AbstractScenario}, w::Integer, stage::Integer = 2; defer::Bool = false)

Store the collection of second stage `scenarios` in in worker node `w` of the `stochasticprogram` at `stage`. Defaults to the second stage.

If `defer` is true, then model creation is deferred until `generate!(stochasticprogram)` is called.
"""
function add_worker_scenarios!(stochasticprogram::StochasticProgram, scenarios::Vector{<:AbstractScenario}, w::Integer, stage::Integer = 2; defer::Bool = false)
    add_scenarios!(scenarioproblems(stochasticprogram, stage), scenarios, w)
    invalidate_cache!(stochasticprogram)
    if !defer
        generate!(stochasticprogram)
    end
    return stochasticprogram
end
"""
    add_scenarios!(scenariogenerator::Function, stochasticprogram::StochasticProgram, n::Integer, stage::Integer = 2; defer::Bool = false)

Generate `n` second-stage scenarios using `scenariogenerator`and store in the `stochasticprogram` at `stage`. Defaults to the second stage.

If `defer` is true, then model creation is deferred until `generate!(stochasticprogram)` is called. If the `stochasticprogram` is distributed, scenarios will be distributed evenly across workers.
"""
function add_scenarios!(scenariogenerator::Function, stochasticprogram::StochasticProgram, n::Integer, stage::Integer = 2; defer::Bool = false)
    add_scenarios!(scenariogenerator, scenarioproblems(stochasticprogram, stage), n)
    invalidate_cache!(stochasticprogram)
    if !defer
        generate!(stochasticprogram)
    end
    return stochasticprogram
end
"""
    add_worker_scenarios!(scenariogenerator::Function, stochasticprogram::StochasticProgram, n::Integer, w::Integer, stage::Integer = 2; defer::Bool = false)

Generate `n` second-stage scenarios using `scenariogenerator`and store them in worker node `w` of the `stochasticprogram` at `stage`. Defaults to the second stage.

If `defer` is true, then model creation is deferred until `generate!(stochasticprogram)` is called.
"""
function add_worker_scenarios!(scenariogenerator::Function, stochasticprogram::StochasticProgram, n::Integer, w::Integer, stage::Integer = 2; defer::Bool = false)
    add_scenarios!(scenariogenerator, scenarioproblems(stochasticprogram, stage), n, w)
    invalidate_cache!(stochasticprogram)
    if !defer
        generate!(stochasticprogram)
    end
    return stochasticprogram
end
"""
    sample!(stochasticprogram::StochasticProgram, sampler::AbstractSampler, n::Integer, stage::Integer = 2; defer::Bool = false)

Sample `n` scenarios using `sampler` and add to the `stochasticprogram` at `stage`. Defaults to the second stage.

If `defer` is true, then model creation is deferred until `generate!(stochasticprogram)` is called. If the `stochasticprogram` is distributed, scenarios will be distributed evenly across workers.
"""
function sample!(stochasticprogram::StochasticProgram, sampler::AbstractSampler, n::Integer, stage::Integer = 2; defer::Bool = false)
    sample!(scenarioproblems(stochasticprogram, stage), sampler, n)
    if !defer
        generate!(stochasticprogram)
    end
    return stochasticprogram
end
# ========================== #


text:
module DeepcopyModules

include("deepcopy_functions.jl")
include("deepcopy.jl")

export deepcopy_module

end # module

ext: jl


text:
using DelayDiffEq, DiffEqProblemLibrary.DDEProblemLibrary
using Test

DDEProblemLibrary.importddeproblems()

@testset "Constant delays" begin
    prob = DDEProblemLibrary.prob_dde_constant_2delays_ip
    prob_scalar = DDEProblemLibrary.prob_dde_constant_2delays_scalar
    # disable reuse
    nlsolve = nlsolve=NLNewton(fast_convergence_cutoff=0)

    algdict = Dict(BS3() => 2.4e-6,
                   Tsit5() => 4.5e-3,
                   RK4() => 1.1e-4,
                   Vern6() => 1.0e-3,
                   SDIRK2(nlsolve=nlsolve) => 2.3e-1,
                   TRBDF2(nlsolve=nlsolve) => 6.2e-2,
                   KenCarp4(nlsolve=nlsolve) => 5.6e-2,
                   Rosenbrock23() => 6.5e-4,
                   Rodas4() => 5.4e-4
                  )

    for (alg, error) in algdict
        ddealg = MethodOfSteps(alg)

        sol = solve(prob, ddealg)
        @test sol.errors[:l∞] < error

        sol_scalar = solve(prob_scalar, ddealg)
        @test sol.t ≈ sol_scalar.t
        @test sol[1, :] ≈ sol_scalar.u
    end
end


text:
# module for calculating std of sensor observation
# observe a landmark by a sensor
# iterate for 1000 times

module SensorTestStatic
  using DataFrames, Statistics

  include(joinpath(split(@__FILE__, "src")[1], "src/model/world/world.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/robot/differential_wheeled_robot/differential_wheeled_robot.jl"))

  function main()
    m = Map()
    add_object(m, Object(1.0, 0.0))

    distances = []
    directions = []

    for i in 1:1000
      s = Sensor(m, dist_noise_rate=0.1,
                 dir_noise=pi/90, dist_bias_rate_stddev=0.1,
                 dir_bias_stddev=pi/90)
      data(s, [0.0, 0.0, 0.0])
      if length(s.last_data) > 0
        push!(distances, s.last_data[1][1][1])
        push!(directions, s.last_data[1][1][2])
      end
    end

    df = DataFrame(distance=distances, direction=directions)
    println("Observation data: ", df)
    println("standard deviation of distance: ", std(df.distance, corrected=false))
    println("standard deviation of direction: ", std(df.direction, corrected=false))
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
using ADCME

save_tensor = load_op_and_grad("./build/libSaveTensor","save_tensor")
get_tensor = load_op_and_grad("./build/libGetTensor","get_tensor")
delete_tensor = load_op_and_grad("./build/libDeleteTensor","delete_tensor")

val = constant(rand(10))
t1 = tf.constant("tensor1")
t2 = tf.constant("tensor2")
t3 = tf.constant("tensor3")
u1 = save_tensor(t1,val)
u2 = save_tensor(t2,2*val)
u3 = save_tensor(t3,3*val)

z1 = get_tensor(t1);
z2 = get_tensor(t2);
z3 = get_tensor(t3);

d1 = delete_tensor(t1);
d2 = delete_tensor(t2);
d3 = delete_tensor(t3);
sess = Session();
run(sess, [u1,u2,u3])

run(sess, z1)
run(sess, z2)
run(sess, z3)
run(sess, d2)


text:
using Test
using NeXLCore

@testset "Electron" begin

    @test isapprox(λₑ(1.0), 1.227e-7, atol = 0.001e-7)
    @test isapprox(λₑ(100.0), 1.227e-8, atol = 0.001e-8)
    @test isapprox(λₑ(1.0e4), 1.227e-9, atol = 0.001e-9)
    @test isapprox(λₑ(5.0e3), 0.0173 * 1.0e-7, atol = 0.0001 * 1.0e-7)
    @test isapprox(λₑ(10.0e3), 0.0122 * 1.0e-7, atol = 0.0001 * 1.0e-7)
    @test isapprox(λₑ(200.0e3), 0.0027 * 1.0e-7, atol = 0.0001 * 1.0e-7)

    @test isapprox(
        kₑ(1.0),
        0.01 * sqrt(2.0 * 9.109e-31 * 1.602e-19 * 1.0) / 1.05457e-34,
        rtol = 1.0e-3,
    )
    @test isapprox(
        kₑ(10.0),
        0.01 * sqrt(2.0 * 9.109e-31 * 1.602e-19 * 10.0) / 1.05457e-34,
        rtol = 1.0e-3,
    )
    @test isapprox(
        kₑ(20.0),
        0.01 * sqrt(2.0 * 9.109e-31 * 1.602e-19 * 20.0) / 1.05457e-34,
        rtol = 1.0e-3,
    )
    @test isapprox(
        kₑ(200.0),
        0.01 * sqrt(2.0 * 9.109e-31 * 1.602e-19 * 200.0) / 1.05457e-34,
        rtol = 1.0e-3,
    )
    @test isapprox(
        kₑ(2000.0),
        0.01 * sqrt(2.0 * 9.109e-31 * 1.602e-19 * 2000.0) / 1.05457e-34,
        rtol = 1.0e-3,
    )

    @test isapprox(mₑ, 0.511e6, atol = 0.001e6)  # Must be in eV!!!!
end


text:
#!/usr/bin/env julia

"""
# Problem 7: 10001st prime

By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that
the 6th prime is 13.

What is the `n`th prime number?
"""
module Problem007

include("Primes.jl")

"""
    solve(; n::Integer = 10001)

Returns the solution for [`Problem007`](@ref) with the given parameters.
"""
function solve(; n::Integer = 10001)
    return Primes.prime(n)
end

end

max_stars_repo_name: curtislb/ProjectEuler

max_issues_repo_name: curtislb/ProjectEuler

max_forks_repo_name: curtislb/ProjectEuler


text:
using open62541
using Dates
using Printf

client = UA_Client_new()
UA_ClientConfig_setDefault(UA_Client_getConfig(client))
retval = UA_Client_connect(client, "opc.tcp://127.0.0.1:4840")

nodeid = UA_NODEID_NUMERIC(0, UA_NS0ID_SERVER_SERVERSTATUS_CURRENTTIME)

value = UA_Client_readValueAttribute(client, nodeid)
UA_Variant_hasScalarType(value, UA_TYPES_PTRS[UA_TYPES_DATETIME])
raw_date = unsafe_wrap(value)

dts = UA_DateTime_toStruct(raw_date)
Printf.@printf("date is: %u-%u-%u %u:%u:%u.%03u\n", dts.day, dts.month, dts.year, dts.hour, dts.min, dts.sec, dts.milliSec)


text:
# Generate a functional interface from the Validator and Imputor types.
"""
    _splitkwargs(::Type{T}, kwargs...) where T -> (imp, rem)

Takes a type with kwargs and returns the constructed type and the
unused kwargs which should be passed to the `impute!` call.

NOTE: This is only intended to be used internally
"""
function _splitkwargs(::Type{T}, kwargs...) where T
    rem = Dict(kwargs...)
    kwdef = empty(rem)

    for f in fieldnames(T)
        if haskey(rem, f)
            kwdef[f] = pop!(rem, f)
        end
    end

    return (T(; kwdef...), rem)
end

# Specialcase kwargs constructor for substitute.
# TODO: Add an imputor method that types should overwrite when necessary or have it fallback to `fieldnames`
function _splitkwargs(::Type{Substitute}, kwargs...)
    rem = Dict(kwargs...)
    kwdef = empty(rem)

    for f in (:statistic, :weights)
        if haskey(rem, f)
            kwdef[f] = pop!(rem, f)
        end
    end

    return (Substitute(; kwdef...), rem)
end

const global validation_methods = (
    threshold = Threshold,
    wthreshold = WeightedThreshold,
)

const global imputation_methods = (
    dropobs = DropObs,
    dropvars = DropVars,
    interp = Interpolate,
    interpolate = Interpolate,
    fill = Fill,
    locf = LOCF,
    nocb = NOCB,
    replace = Replace,
    srs = SRS,
    substitute = Substitute,
    wsubstitute = WeightedSubstitute,
    svd = SVD,
    knn = KNN,
)

for (func, type) in pairs(validation_methods)
    typename = nameof(type)
    @eval begin
        function $func(data; kwargs...)
            a, rem = _splitkwargs($typename, kwargs...)
            return validate(data, a; rem...)
        end
    end
end

for (func, type) in pairs(imputation_methods)
    typename = nameof(type)
    func! = Symbol(func, :!)

    @eval begin
        function $func(data; kwargs...)
            imp, rem = _splitkwargs($typename, kwargs...)
            return impute(data, imp; rem...)
        end
        function $func!(data; kwargs...)
            imp, rem = _splitkwargs($typename, kwargs...)
            return impute!(data, imp; rem...)
        end
        @deprecate $func(; kwargs...) data -> $func(data; kwargs...) false
        @deprecate $func!(; kwargs...) data -> $func!(data; kwargs...) false
    end
end

declaremissings(data; kwargs...) = apply(data, DeclareMissings(; kwargs...))
declaremissings!(data; kwargs...) = apply!(data, DeclareMissings(; kwargs...))

# Provide a specific functional API for Impute.Filter.
filter(data; kwargs...) = apply(data, Filter(); kwargs...)
filter!(data; kwargs...) = apply!(data, Filter(); kwargs...)
filter(f::Function, data; kwargs...) = apply(data, Filter(f); kwargs...)
filter!(f::Function, data; kwargs...) = apply!(data, Filter(f); kwargs...)

@doc """
    Impute.threshold(data; limit=0.1, kwargs...)

Assert that proportion of missing values in the `data` do not exceed the `limit`.

# Examples
```julia-repl
julia> using DataFrames, Impute

julia> df = DataFrame(:a => [1.0, 2.0, missing, missing, 5.0], :b => [1.1, 2.2, 3.3, missing, 5.5])
5×2 DataFrames.DataFrame
│ Row │ a        │ b        │
│     │ Float64  │ Float64  │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ missing  │ 3.3      │
│ 4   │ missing  │ missing  │
│ 5   │ 5.0      │ 5.5      │

julia> Impute.threshold(df)
ERROR: ThresholdError: Missing data limit exceeded 0.1 (0.4)
Stacktrace:
...

julia> Impute.threshold(df; limit=0.8)
5×2 DataFrames.DataFrame
│ Row │ a        │ b        │
│     │ Float64  │ Float64  │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ missing  │ 3.3      │
│ 4   │ missing  │ missing  │
│ 5   │ 5.0      │ 5.5      │
```
"""
threshold

@doc """
    Impute.wthreshold(data; ratio, weights, kwargs...)

Assert that the weighted proportion of missing values in the `data` do not exceed the `limit`.

# Examples
```julia-repl
julia> using DataFrames, Impute

julia> df = DataFrame(:a => [1.0, 2.0, missing, missing, 5.0], :b => [1.1, 2.2, 3.3, missing, 5.5])
5×2 DataFrames.DataFrame
│ Row │ a        │ b        │
│     │ Float64  │ Float64  │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ missing  │ 3.3      │
│ 4   │ missing  │ missing  │
│ 5   │ 5.0      │ 5.5      │

julia> Impute.wthreshold(df; limit=0.4, weights=0.1:0.1:0.5)
ERROR: ThresholdError: Missing data limit exceeded 0.4 (0.4666666666666666)
Stacktrace:
...

julia> Impute.wthreshold(df; limit=0.4, weights=0.5:-0.1:0.1)
5×2 DataFrames.DataFrame
│ Row │ a        │ b        │
│     │ Float64  │ Float64  │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ missing  │ 3.3      │
│ 4   │ missing  │ missing  │
│ 5   │ 5.0      │ 5.5      │
```
"""
wthreshold

@doc """
    Impute.dropobs(data; dims=1)

Removes missing observations from the `AbstractArray` or `Tables.table` provided.
See [DropObs](@ref) for details.

!!! Use `Impute.filter(data; dims=1)` instead.

# Example
```julia-repl
julia> using DataFrames; using Impute: Impute

julia> df = DataFrame(:a => [1.0, 2.0, missing, missing, 5.0], :b => [1.1, 2.2, 3.3, missing, 5.5])
5×2 DataFrames.DataFrame
│ Row │ a        │ b        │
│     │ Float64  │ Float64  │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ missing  │ 3.3      │
│ 4   │ missing  │ missing  │
│ 5   │ 5.0      │ 5.5      │

julia> Impute.dropobs(df; dims=2)
3×2 DataFrames.DataFrame
│ Row │ a       │ b       │
│     │ Float64 │ Float64 │
├─────┼─────────┼─────────┤
│ 1   │ 1.0     │ 1.1     │
│ 2   │ 2.0     │ 2.2     │
│ 3   │ 5.0     │ 5.5     │
```
""" dropobs

@doc """
    Impute.dropvars(data; dims=1)

Finds variables with missing values in a `AbstractMatrix` or `Tables.table` and
removes them from the input data. See [DropVars](@ref) for details.

!!! Use `Impute.filter(data; dims=2)` instead.

# Example
```julia-repl
julia> using DataFrames; using Impute: Impute

julia> df = DataFrame(:a => [1.0, 2.0, missing, missing, 5.0], :b => [1.1, 2.2, 3.3, missing, 5.5])
5×2 DataFrame
│ Row │ a        │ b        │
│     │ Float64? │ Float64? │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ missing  │ 3.3      │
│ 4   │ missing  │ missing  │
│ 5   │ 5.0      │ 5.5      │

julia> Impute.dropvars(df)
0×0 DataFrame
```
""" dropvars

@doc """
    Impute.filter([f,] data; dims)

Filters values, rows, columns or slices of data that should be removed.
The default function `f` will removing `missing`s, or any rows, columns or slices
containing `missing`s.

# Examples
```jldoctest
julia> using DataFrames; using Impute: Impute

julia> df = DataFrame(:a => [1.0, 2.0, missing, missing, 5.0], :b => [1.1, 2.2, 3.3, missing, 5.5])
5×2 DataFrame
 Row │ a          b
     │ Float64?   Float64?
─────┼──────────────────────
   1 │       1.0        1.1
   2 │       2.0        2.2
   3 │ missing          3.3
   4 │ missing    missing
   5 │       5.0        5.5

julia> Impute.filter(df; dims=:cols)
0×0 DataFrame

julia> Impute.filter(df; dims=:rows)
3×2 DataFrame
 Row │ a        b
     │ Float64  Float64
─────┼──────────────────
   1 │     1.0      1.1
   2 │     2.0      2.2
   3 │     5.0      5.5
```
""" filter

@doc """
    Impute.interp(data; dims=1)

Performs linear interpolation between the nearest values in an vector.
See [`Impute.Interpolate`](@ref) for details.

# Example
```jldoctest
julia> using DataFrames; using Impute: Impute

julia> df = DataFrame(:a => [1.0, 2.0, missing, missing, 5.0], :b => [1.1, 2.2, 3.3, missing, 5.5])
5×2 DataFrame
 Row │ a          b
     │ Float64?   Float64?
─────┼──────────────────────
   1 │       1.0        1.1
   2 │       2.0        2.2
   3 │ missing          3.3
   4 │ missing    missing
   5 │       5.0        5.5

julia> Impute.interp(df)
5×2 DataFrame
 Row │ a         b
     │ Float64?  Float64?
─────┼────────────────────
   1 │      1.0       1.1
   2 │      2.0       2.2
   3 │      3.0       3.3
   4 │      4.0       4.4
   5 │      5.0       5.5
```
""" interp

@doc """
    Impute.fill(data; value=mean, dims=1)

Fills in the missing data with a specific value. See [Fill](@ref) for details.

!!! Use `Impute.replace` for constants or `Impute.substitue` for median/mode substitution.

# Example
```julia-repl
julia> using DataFrames; using Impute: Impute

julia> df = DataFrame(:a => [1.0, 2.0, missing, missing, 5.0], :b => [1.1, 2.2, 3.3, missing, 5.5])
5×2 DataFrame
│ Row │ a        │ b        │
│     │ Float64? │ Float64? │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ missing  │ 3.3      │
│ 4   │ missing  │ missing  │
│ 5   │ 5.0      │ 5.5      │

julia> Impute.fill(df; value=-1.0)
5×2 DataFrame
│ Row │ a        │ b        │
│     │ Float64? │ Float64? │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ -1.0     │ 3.3      │
│ 4   │ -1.0     │ -1.0     │
│ 5   │ 5.0      │ 5.5      │
```
""" fill

@doc """
    Impute.locf(data; dims=1)

Iterates forwards through the `data` and fills missing data with the last existing
observation. See [`Impute.LOCF`](@ref) for details.

# Example
```jldoctest
julia> using DataFrames; using Impute: Impute

julia> df = DataFrame(:a => [1.0, 2.0, missing, missing, 5.0], :b => [1.1, 2.2, 3.3, missing, 5.5])
5×2 DataFrame
 Row │ a          b
     │ Float64?   Float64?
─────┼──────────────────────
   1 │       1.0        1.1
   2 │       2.0        2.2
   3 │ missing          3.3
   4 │ missing    missing
   5 │       5.0        5.5

julia> Impute.locf(df)
5×2 DataFrame
 Row │ a         b
     │ Float64?  Float64?
─────┼────────────────────
   1 │      1.0       1.1
   2 │      2.0       2.2
   3 │      2.0       3.3
   4 │      2.0       3.3
   5 │      5.0       5.5
```
""" locf

@doc """
    Impute.nocb(data; dims=1)

Iterates backwards through the `data` and fills missing data with the next existing
observation. See [`Impute.NOCB`](@ref) for details.

# Example
```jldoctest
julia> using DataFrames; using Impute: Impute

julia> df = DataFrame(:a => [1.0, 2.0, missing, missing, 5.0], :b => [1.1, 2.2, 3.3, missing, 5.5])
5×2 DataFrame
 Row │ a          b
     │ Float64?   Float64?
─────┼──────────────────────
   1 │       1.0        1.1
   2 │       2.0        2.2
   3 │ missing          3.3
   4 │ missing    missing
   5 │       5.0        5.5

julia> Impute.nocb(df)
5×2 DataFrame
 Row │ a         b
     │ Float64?  Float64?
─────┼────────────────────
   1 │      1.0       1.1
   2 │      2.0       2.2
   3 │      5.0       3.3
   4 │      5.0       5.5
   5 │      5.0       5.5
```
""" nocb

@doc """
    Impute.srs(data; rng=Random.GLOBAL_RNG)

Simple Random Sampling (SRS) imputation is a method for imputing both continuous and
categorical variables. Furthermore, it completes imputation while preserving the
distributional properties of the variables (e.g., mean, standard deviation).

# Example
```julia-repl
julia> using DataFrames; using Random; using Impute: Impute

julia> df = DataFrame(:a => [1.0, 2.0, missing, missing, 5.0], :b => [1.1, 2.2, 3.3, missing, 5.5])
5×2 DataFrame
│ Row │ a        │ b        │
│     │ Float64? │ Float64? │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ missing  │ 3.3      │
│ 4   │ missing  │ missing  │
│ 5   │ 5.0      │ 5.5      │

julia> Impute.srs(df; rng=MersenneTwister(1234))
5×2 DataFrame
│ Row │ a        │ b        │
│     │ Float64? │ Float64? │
├─────┼──────────┼──────────┤
│ 1   │ 1.0      │ 1.1      │
│ 2   │ 2.0      │ 2.2      │
│ 3   │ 1.0      │ 3.3      │
│ 4   │ 2.0      │ 3.3      │
│ 5   │ 5.0      │ 5.5      │
```
""" srs

@doc """
    Impute.declaremissings(data; values)

DeclareMissings (or replace) various missing data representations with `missing`.

# Keyword Arguments
* `value::Tuple`: A tuple of values that should be considered `missing`

# Example
```jldoctest
julia> using DataFrames, Impute

julia> df = DataFrame(
           :a => [1.1, 2.2, NaN, NaN, 5.5],
           :b => [1, 2, 3, -9999, 5],
           :c => ["v", "w", "x", "y", "NULL"],
       )
5×3 DataFrame
 Row │ a        b      c
     │ Float64  Int64  String
─────┼────────────────────────
   1 │     1.1      1  v
   2 │     2.2      2  w
   3 │   NaN        3  x
   4 │   NaN    -9999  y
   5 │     5.5      5  NULL

julia> Impute.declaremissings(df; values=(NaN, -9999, "NULL"))
5×3 DataFrame
 Row │ a          b        c
     │ Float64?   Int64?   String?
─────┼─────────────────────────────
   1 │       1.1        1  v
   2 │       2.2        2  w
   3 │ missing          3  x
   4 │ missing    missing  y
   5 │       5.5        5  missing
```
""" declaremissings

@doc """
    Impute.replace(data; values)

Replace `missing`s with one of the specified constant values, depending on the input type.
If multiple values of the same type are provided then the first one will be used.
If the input data is of a different type then the no replacement will be performed.

# Keyword Arguments
* `values::Tuple`: A scalar or tuple of different values that should be used to replace
  missings. Typically, one value per type you're considering imputing for.

# Example
```jldoctest
julia> using DataFrames, Impute

julia> df = DataFrame(
           :a => [1.1, 2.2, missing, missing, 5.5],
           :b => [1, 2, 3, missing, 5],
           :c => ["v", "w", "x", "y", missing],
       )
5×3 DataFrame
 Row │ a          b        c
     │ Float64?   Int64?   String?
─────┼─────────────────────────────
   1 │       1.1        1  v
   2 │       2.2        2  w
   3 │ missing          3  x
   4 │ missing    missing  y
   5 │       5.5        5  missing

julia> Impute.replace(df; values=(NaN, -9999, "NULL"))
5×3 DataFrame
 Row │ a         b       c
     │ Float64?  Int64?  String?
─────┼───────────────────────────
   1 │      1.1       1  v
   2 │      2.2       2  w
   3 │    NaN         3  x
   4 │    NaN     -9999  y
   5 │      5.5       5  NULL
```
""" replace

@doc """
    Impute.substitute(data; statistic=nothing)
    Impute.substitute(data; weights=nothing)

Substitute missing values with a summary statistic over the non-missing values.

# Keyword Arguments
* `statistic`: A summary statistic function to be applied to the non-missing values.
  This function should return a value of the same type as the input data `eltype`.
  If this function isn't passed in then the `defaultstats` function is used to make a
  best guess.
* `weights`: A set of statistical weights to apply to the `mean` or `median` in `defaultstats`.

See [Substitute](@ref) for details on substitution rules defined in `defaultstats`.

# Example
```jldoctest
julia> using DataFrames, Impute

julia> df = DataFrame(
                  :a => [8.9, 2.2, missing, missing, 1.3, 6.2, 3.7, 4.8],
                  :b => [2, 6, 3, missing, 7, 1, 9, missing],
                  :c => [true, false, true, true, false, missing, false, true],
              )
8×3 DataFrame
 Row │ a          b        c
     │ Float64?   Int64?   Bool?
─────┼─────────────────────────────
   1 │       8.9        2     true
   2 │       2.2        6    false
   3 │ missing          3     true
   4 │ missing    missing     true
   5 │       1.3        7    false
   6 │       6.2        1  missing
   7 │       3.7        9    false
   8 │       4.8  missing     true

julia> Impute.substitute(df)
8×3 DataFrame
 Row │ a         b       c
     │ Float64?  Int64?  Bool?
─────┼─────────────────────────
   1 │     8.9        2   true
   2 │     2.2        6  false
   3 │     4.25       3   true
   4 │     4.25       4   true
   5 │     1.3        7  false
   6 │     6.2        1   true
   7 │     3.7        9  false
   8 │     4.8        4   true
```
""" substitute

@doc """
    Impute.knn(; k=1, threshold=0.5, dist=Euclidean())

Imputation using k-Nearest Neighbor algorithm.

# Keyword Arguments
* `k::Int`: number of nearest neighbors
* `dist::MinkowskiMetric`: distance metric suppports by `NearestNeighbors.jl` (Euclidean, Chebyshev, Minkowski and Cityblock)
* `threshold::AbsstractFloat`: thershold for missing neighbors

# Reference
* Troyanskaya, Olga, et al. "Missing value estimation methods for DNA microarrays." Bioinformatics 17.6 (2001): 520-525.

# Example
```jldoctest
julia> using Impute, Missings

julia> data = allowmissing(reshape(sin.(1:20), 5, 4)); data[[2, 3, 7, 9, 13, 19]] .= missing; data
5×4 Matrix{Union{Missing, Float64}}:
  0.841471  -0.279415  -0.99999   -0.287903
   missing    missing  -0.536573  -0.961397
   missing   0.989358    missing  -0.750987
 -0.756802    missing   0.990607    missing
 -0.958924  -0.544021   0.650288   0.912945

julia> result = Impute.knn(data; dims=:cols)
5×4 Matrix{Union{Missing, Float64}}:
  0.841471  -0.279415  -0.99999   -0.287903
 -0.756802   0.989358  -0.536573  -0.961397
 -0.756802   0.989358  -0.536573  -0.750987
 -0.756802  -0.544021   0.990607   0.912945
 -0.958924  -0.544021   0.650288   0.912945
```
""" knn

@doc """
    Impute.svd(; kwargs...)

Imputes the missing values in a matrix using an expectation maximization (EM) algorithm
over low-rank SVD approximations.

# Keyword Arguments
* `init::Imputor`: initialization method for missing values (default: Substitute())
* `rank::Union{Int, Nothing}`: rank of the SVD approximation (default: nothing meaning start and 0 and increase)
* `tol::Float64`: convergence tolerance (default: 1e-10)
* `maxiter::Int`: Maximum number of iterations if convergence is not achieved (default: 100)
* `limits::Unoin{Tuple{Float64, Float64}, Nothing}`: Bound the possible approximation values (default: nothing)
* `verbose::Bool`: Whether to display convergence progress (default: true)

# References
* Troyanskaya, Olga, et al. "Missing value estimation methods for DNA microarrays." Bioinformatics 17.6 (2001): 520-525.

# Example
```jldoctest
julia> using Impute, Missings

julia> data = allowmissing(reshape(sin.(1:20), 5, 4)); data[[2, 3, 7, 9, 13, 19]] .= missing; data
5×4 Matrix{Union{Missing, Float64}}:
  0.841471  -0.279415  -0.99999   -0.287903
   missing    missing  -0.536573  -0.961397
   missing   0.989358    missing  -0.750987
 -0.756802    missing   0.990607    missing
 -0.958924  -0.544021   0.650288   0.912945

julia> result = Impute.svd(data; dims=:cols)
5×4 Matrix{Union{Missing, Float64}}:
  0.841471  -0.279415  -0.99999   -0.287903
  0.220258   0.555829  -0.536573  -0.961397
 -0.372745   0.989358   0.533193  -0.750987
 -0.756802   0.253309   0.990607   0.32315
 -0.958924  -0.544021   0.650288   0.912945
```
""" svd


text:
"""
    StaticOrdinalRange

"""
abstract type StaticOrdinalRange{T,B,S,E,L} <: StaticStartRange{T,B,E,L} end

@pure step(::StaticOrdinalRange{T,B,S}) where {T,B,S<:SVal} = S()
@pure step(::Type{<:StaticOrdinalRange{T,B,S}}) where {T,B,S<:SVal} = S()

@pure step(::StaticOrdinalRange{T,B,TPVal{H,L}}) where {T,B,H,L} = H() + L()
@pure step(::Type{<:StaticOrdinalRange{T,B,TPVal{H,L}}}) where {T,B,H,L} = H() + L()

@inline step(r::StaticStartRange{T,B,S}) where {T,B,S<:BaseNumber} = getfield(r, :step)::S
@inline step(r::StaticStartRange{T,B,S}) where {T,B,S<:TwicePrecision} = T(getfield(r, :step))

function length(r::StaticOrdinalRange{T,B,S,E,Dynamic}) where {T<:Union{Int,UInt,Int64,UInt64,Int128,UInt128},B,S,E}
    isempty(r) && return zero(T)
    if step(r) > 1
        return convert(T, div(unsigned(last(r) - first(r)), step(r))) + one(T)
    elseif step(r) < -1
        return convert(T, div(unsigned(first(r) - last(r)), -step(r))) + one(T)
    elseif step(r) > 0
        return div((last(r) - first(r)), step(r)) + one(T)
    else
        return div(first(r) - last(r), -step(r)) + one(T)
    end
end

function length(r::StaticOrdinalRange{T,B,S,E,Dynamic}) where {T,B,S,E}
    n = Integer(div((last(r) - first(r)) + step(r), step(r)))
    isempty(r) ? zero(n) : n
end

showrange(io::IO, r::StaticOrdinalRange) = print(io, "$(first(r)):$(step(r)):$(last(r))")


text:
@polyvar x y
@variable(model, s)
@variable(model, t)
p = s*x^6+t*y^6-x^4*y^2-x^2*y^4-x^4+3*x^2*y^2-y^4-x^2-y^2+1
@constraint(model, p >= 0)
@objective(model, Min, s+t)
optimize!(model)
println("Solution: [ $(value(s)), $(value(t)) ]")
# returns s ~ 1.203, t ~ 1.203

max_stars_repo_name: sums-of-squares/tools

max_issues_repo_name: sums-of-squares/tools

max_forks_repo_name: sums-of-squares/tools


text:
module GalacticBBO

using BlackBoxOptim, GalacticOptim, GalacticOptim.SciMLBase

abstract type BBO end

for j = string.(BlackBoxOptim.SingleObjectiveMethodNames)
    eval(Meta.parse("Base.@kwdef struct BBO_" * j * " <: BBO method=:" * j * " end"))
    eval(Meta.parse("export BBO_" * j))
end

decompose_trace(opt::BlackBoxOptim.OptRunController) = BlackBoxOptim.best_candidate(opt)

function __map_optimizer_args(prob::SciMLBase.OptimizationProblem, opt::BBO;
    callback=nothing,
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    kwargs...)

    if !isnothing(reltol)
        @warn "common reltol is currently not used by $(opt)"
    end

    mapped_args = (; Method=opt.method,
        SearchRange=[(prob.lb[i], prob.ub[i]) for i in 1:length(prob.lb)])

    if !isnothing(callback)
        mapped_args = (; mapped_args..., CallbackFunction=callback, CallbackInterval=0.0)
    end

    mapped_args = (; mapped_args..., kwargs...)

    if !isnothing(maxiters)
        mapped_args = (; mapped_args..., MaxSteps=maxiters)
    end

    if !isnothing(maxtime)
        mapped_args = (; mapped_args..., MaxTime=maxtime)
    end

    if !isnothing(abstol)
        mapped_args = (; mapped_args..., MinDeltaFitnessTolerance=abstol)
    end

    return mapped_args
end

function SciMLBase.__solve(prob::SciMLBase.OptimizationProblem, opt::BBO, data=GalacticOptim.DEFAULT_DATA;
    callback=(args...) -> (false),
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    progress=false, kwargs...)

    local x, cur, state

    if data != GalacticOptim.DEFAULT_DATA
        maxiters = length(data)
    end

    cur, state = iterate(data)

    function _cb(trace)
        cb_call = callback(decompose_trace(trace), x...)
        if !(typeof(cb_call) <: Bool)
            error("The callback should return a boolean `halt` for whether to stop the optimization process.")
        end
        if cb_call == true
            BlackBoxOptim.shutdown_optimizer!(trace) #doesn't work
        end
        cur, state = iterate(data, state)
        cb_call
    end

    maxiters = GalacticOptim._check_and_convert_maxiters(maxiters)
    maxtime = GalacticOptim._check_and_convert_maxtime(maxtime)

    _loss = function (θ)
        x = prob.f(θ, prob.p, cur...)
        return first(x)
    end

    opt_args = __map_optimizer_args(prob, opt, callback=_cb, maxiters=maxiters, maxtime=maxtime, abstol=abstol, reltol=reltol; kwargs...)

    opt_setup = BlackBoxOptim.bbsetup(_loss; opt_args...)

    t0 = time()
    opt_res = BlackBoxOptim.bboptimize(opt_setup)
    t1 = time()

    opt_ret = Symbol(opt_res.stop_reason)

    SciMLBase.build_solution(prob, opt, BlackBoxOptim.best_candidate(opt_res),
        BlackBoxOptim.best_fitness(opt_res); original=opt_res, retcode=opt_ret)
end

end


text:
#=
R ploting function
=#

function create_plot_ClaDS_ape()
    reval("""

    list.of.packages <- c("ape", "fields", "RColorBrewer")
    new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
    if(length(new.packages)){
		print(paste0("package ", list.of.packages, " should be installed in R."))
		print(paste0("Please install these packages within a R session before continuing."))
		}
	require(fields)
	require(ape)
	require(RColorBrewer)

    plot_ClaDS=function(phylo,rate1,rate2=NULL,same.scale=T,main=NULL,lwd=3,log=T, show_labels=F, minr = Inf, maxr = 0, show_legend = T,...){
        Colors = colorRampPalette(rev(c('darkred',brewer.pal(n = 8, name = "Spectral"),'darkblue')))(100)

     if(nrow(phylo[[1]]) <=1){
        plot(1000,xlim = c(0,1), ylim = c(0,2), axes = F, xlab = "", ylab = "")
        lines(c(0,1), c(1,1), lwd=lwd, col="steelblue2")
      }else{
        if(is.null(rate2)){
          if(log) {
              rate1=log(rate1)
              minr = log(minr)
              maxr = log(maxr)
              minr = min(c(minr, log(0.95)+max(rate1),min(rate1)))
              maxr = max(c(maxr, log(1.05)+min(rate1), max(rate1)))
              }else{
                minr = min(c(minr, 0.99*max(rate1),min(rate1)))
                maxr = max(c(maxr, 1.01*min(rate1), max(rate1)))
              }

            col = round( (rate1 - minr) / (maxr - minr)*99   )+1
            plot(phylo, edge.color = Colors[col], show.tip.label = show_labels,main=main,edge.width =lwd,...)

            if(log){
              m10=floor(minr/log(10))
              M10=ceiling(maxr/log(10))
              if((M10-m10)<4){
                ticks=c(1,2,5)
              }else{
                ticks=1
              }

              ticks=as.vector(sapply(m10:M10,function(k){return(ticks*10^k)}))
              lt=length(ticks[ticks>exp(minr) & ticks<exp(maxr)])
              if(lt<4){
                  ticks=c(round(exp(minr),digit=2),max(0,-1*m10+(lt<2)),ticks,round(exp(maxr),digit = 2))
                  }
              if(show_legend) image.plot(z = c(minr,maxr),col = Colors, horizontal=T,legend.only = T,axis.args=list( at=log(ticks), labels=ticks))
            }else{
              if(show_legend) image.plot(z = as.matrix(rate1),col = Colors, horizontal=T,legend.only = T)
            }

        }else{
          if(log){
            rate1=log(rate1)
            rate2=log(rate2)
          }
          if(same.scale){
            min=min(min(rate1),min(rate2))
            max=max(max(rate1),max(rate2))
            par(mfrow=c(1,2))
            col = round(( (rate1 - min) / (max-min))*99   )+1
            plot(phylo, edge.color = Colors[col], show.tip.label = show_labels,edge.width =lwd,...)
            col = round(( (rate2 - min) / (max-min))*99   )+1
            plot(phylo, edge.color = Colors[col], show.tip.label = show_labels,edge.width =lwd,...)
            par(mfrow=c(1,1))
            if(log){
              m10=floor(min/log(10))
              M10=ceiling(max/log(10))
              if((M10-m10)<4){
                ticks=c(1,2,5)
              }else{
                ticks=1
              }
              ticks=as.vector(sapply(m10:M10,function(k){return(ticks*10^k)}))
              lt=length(ticks[ticks>exp(min) & ticks<exp(max)])
              if(lt<4){

                  ticks=c(round(exp(min),max(0,-1*m10+(lt<2))),ticks,round(exp(max),max(0,-1*M10+1+(lt<2))))
                  }
              # ticks=seq(min,max,length.out = 5)
              if(show_legend) image.plot(z = c(min,max),col = Colors, horizontal=T,legend.only = T,axis.args=list( at=log(ticks), labels=ticks))
            }else{
              if(show_legend) image.plot(z = c(min,max),col = Colors, horizontal=T,legend.only = T)
            }
          }else{
            par(mfrow=c(1,2))
            if(isTRUE(all.equal(rep(rate1[1],length(rate1)),rate1))){
              col=rep(1,length(rate1))
              plot(phylo, edge.color = Colors[col], show.tip.label = show_labels,edge.width =lwd,...)
              if(log){

                if(show_legend) image.plot(z = c(exp(rate1[1]),2*exp(rate1[1])),col = Colors, horizontal=T,legend.only = T)
              }else{
                if(show_legend) image.plot(z = c(rate1[1],2*rate1[1]),col = Colors, horizontal=T,legend.only = T)
              }
            }else{
              col = round(( (rate1 - min(rate1)) / (max(rate1)-min(rate1)))*99   )+1
              plot(phylo, edge.color = Colors[col], show.tip.label = show_labels,edge.width =lwd,...)
              if(log){
                min=min(rate1)
                max=max(rate1)
                m10=floor(min/log(10))
                M10=ceiling(max/log(10))
                if((M10-m10)<4){
                  ticks=c(1,2,5)
                }else{
                  ticks=1
                }
                ticks=as.vector(sapply(m10:M10,function(k){return(ticks*10^k)}))
                lt=length(ticks[ticks>exp(min) & ticks<exp(max)])
                if(lt<4){ticks=c(round(exp(min),max(0,-1*m10+(lt<2))),ticks,round(exp(max),max(0,-1*M10+1+(lt<2))))}
                if(show_legend) image.plot(z = c(min,max),col = Colors, horizontal=T,legend.only = T,axis.args=list( at=log(ticks), labels=ticks))
              }else{
                if(show_legend) image.plot(z = as.matrix(rate1),col = Colors, horizontal=T,legend.only = T)
              }
            }
            if(isTRUE(all.equal(rep(rate2[1],length(rate2)),rate2))){
              col=rep(1,length(rate2))
              plot(phylo, edge.color = Colors[col], show.tip.label = show_labels,edge.width =lwd,...)
              if(log){
                if(show_legend) image.plot(z = c(exp(rate2[1]),2*exp(rate2[1])),col = Colors, horizontal=T,legend.only = T)
              }else{
                if(show_legend) image.plot(z = c(rate2[1],2*rate2[1]),col = Colors, horizontal=T,legend.only = T)
              }
            }else{
              col = round(( (rate2 - min(rate2)) / (max(rate2)-min(rate2)))*99   )+1
              plot(phylo, edge.color = Colors[col], show.tip.label = show_labels,edge.width =lwd,...)
              if(log){
                minr=min(rate2)
                maxr=max(rate2)
                m10=floor(minr/log(10))
                M10=ceiling(maxr/log(10))
                if((M10-m10)<4){
                  ticks=c(1,2,5)
                }else{
                  ticks=1
                }
                ticks=as.vector(sapply(m10:M10,function(k){return(ticks*10^k)}))
                lt=length(ticks[ticks>exp(min) & ticks<exp(max)])
                if(lt<4){ticks=c(round(exp(min),max(0,-1*m10+(lt<2))),ticks,round(exp(max),max(0,-1*M10+1+(lt<2))))}
                if(show_legend) image.plot(z = c(min,max),col = Colors, horizontal=T,legend.only = T,axis.args=list( at=log(ticks), labels=ticks))
              }else{
                if(show_legend) image.plot(z = as.matrix(rate2),col = Colors, horizontal=T,legend.only = T)
              }
            }
          }
          par(mfrow=c(1,1))
        }
      }
      if(!show_legend) return(list(Colors = Colors, col = col, minr = minr, maxr = maxr))
    }
    """)
end

#=
Plot a tree with branch rates from Julia
3 methods
=#

# specifying only the tree, .attribute[id] used as rates
"""
    plot_ClaDS(tree::Tree ; ln=true, show_labels = false, options="")

Plots a tree with branches colored according to their rates

# Arguments
- `tree::Tree`: the phylogeny to be plotted.

# Keyword arguments
- `ln::Bool`: Should rates be plotted on a log scale? Default to true.
- `show_labels::Bool`: Should tip labels be printed? Default to false.
- `options::String`: Additonal options for the ploting function.

	plot_ClaDS(tree::Tree, rates::Array{Number,1} ; ln=true, show_labels = false, options="")

# Additional arguments
- `rates::Array{Number,1} `: the speciation rates.
"""
function plot_ClaDS(tree::Tree ; id = 1, ln=true, show_labels = false, options="")
    plot_tree = Tree(tree.offsprings, 0., tree.attributes, tree.n_nodes)

    opt = options
    if length(options) > 0
        if options[1] != ","[1]
            opt = ", $options"
        end
    end

    edges, branch_lengths, rates, tip_labels = Tree2ape(plot_tree, id = 1)
    ntip = (tree.n_nodes + 1)/2
    @rput edges
    @rput branch_lengths
    @rput rates
    @rput ntip
    @rput ln
    @rput show_labels
    @rput tip_labels

    create_plot_ClaDS_ape()
    reval("""
        tree = list(edge = edges, Nnode = ntip - 1, edge.lengths = branch_lengths, tip.labels = tip_labels)
        class(tree) = "phylo"
        colors = plot_ClaDS(tree, rates, log=ln, show_labels = show_labels $opt);
    """)
end

# specifying the tree and a vector of rates

function plot_ClaDS(tree::Tree, rates::Array{Float64,1} ; id = 1, ln=true, round = false, options = "", show_labels=false)
    plot_tree = Tree(tree.offsprings, 0., tree.attributes, tree.n_nodes)

    opt = options
    if length(options) > 0
        if options[1] != ","[1]
            opt = ", $options"
        end
    end

    edges, branch_lengths, tree_rates, tip_labels = Tree2ape(plot_tree, id = 1)
    ntip = (tree.n_nodes + 1)/2
    @rput edges
    @rput branch_lengths
    @rput rates
    @rput ntip
    @rput ln
    @rput show_labels
    @rput tip_labels

    create_plot_ClaDS_ape()
    reval("""
        tree = list(edge = edges, Nnode = ntip - 1, edge.lengths = branch_lengths, tip.labels = tip_labels)
        class(tree) = "phylo"
        plot_ClaDS(tree, rates, log=ln, show_labels = show_labels $opt)
    """)
end

# diversity through time
function plot_DTT(co::CladsOutput; nyl = 1000, n_ltt = 100, alpha_col = 0.05, burn = 0.25,
	thin = 1, axes_cex = 1., lwd = 1., lab_cex = 1., force_ylim = [2.]) where {T <: Number}

	extant_tree = co.tree
	chains = co.chains
    npar = extant_tree.n_nodes + 3 + (extant_tree.n_nodes+1)/2
    ltt_times = co.time_points
    maps = co.DTT_mean

    new_tree = Tree(extant_tree.offsprings, 0., extant_tree.attributes)
    live_ltt = LTT(new_tree, ltt_times)[2]

    @rput chains
    @rput maps
    @rput npar
    @rput ltt_times
    @rput n_ltt
    @rput alpha_col
    @rput burn
    @rput thin
    @rput live_ltt
	@rput axes_cex
	@rput lab_cex
    @rput nyl
	@rput force_ylim
    reval("""
        require(coda)
        id_ltt = (npar+3):length(chains[[1]])
        n_row = length(chains[[1]][[1]])
        ini = floor(burn * n_row)
        if(ini == 0) ini = 1
        it = seq(ini, n_row, thin)
        npar2 = length(chains[[1]])
        map_chains = mcmc.list(lapply(1:3, function(i){
            mcmc(sapply(id_ltt, function(j){
                if(j<4 || (j>npar )){#}& j<(npar +3))) {
                    chains[[i]][[j]][it]
                }else{
                    log(chains[[i]][[j]][it])
                }
        }))}))
        unlist_chains = (sapply(1:length(id_ltt), function(k){
            c(map_chains[[1]][,k], map_chains[[2]][,k], map_chains[[3]][,k])
            }))
		y_max = max(unlist_chains)
        #maps_log=sapply(1:npar2, function(i){D=density(unlist_chains[,i]);
        #    return(D[[1]][which.max(D[[2]])])})

        means=sapply(1:length(id_ltt), function(i){if(T){log(mean(unlist_chains[,i]))}else{mean(unlist_chains[,i])}})
    """)

    reval("""
        library(scales)
        plot(100000, axes = F, xlim = range(ltt_times), ylim = range(log(c(2,force_ylim,y_max))), xlab = "time",
        ylab = "nb lineages", cex.lab = lab_cex)
        y_lab = c(2,5,10,50,100,200,500,1000,2000,5000,10000,20000,50000)
        y_lab = y_lab[y_lab<=y_max]
        if(nyl < length(y_lab)){
            y_lab = y_lab[unique(floor(seq(1,length(y_lab), length.out = nyl)))]
            }
        lines(ltt_times, log(live_ltt), col = "black", lwd = 6, lty = 1)

        id_plot = unique(floor(seq(ini, n_row, length.out = n_ltt)))#unique(floor(seq(2,length(chains[[1]][[1]]), length.out = n_ltt)))

        Ys = c()
        for (i in 1:3){
            for (j in id_plot){
                y = sapply(id_ltt, function(k){log(chains[[i]][[k]][j])})
                Ys = rbind(Ys,y)
                lines(ltt_times, y, col = alpha("deepskyblue2", alpha = alpha_col), lwd = 2)
                }
            }
		quant = sapply(1:ncol(Ys), function(i){quantile(log(unlist_chains[,i]), probs = c(0.05,0.95))})

        lines(ltt_times, quant[1,], col = "deepskyblue3", lwd = 2)
        lines(ltt_times, quant[2,], col = "deepskyblue3", lwd = 2)
        lines(ltt_times, log(maps[id_ltt]), col = "cadetblue1", lwd = 5, lty = 1)
        lines(ltt_times, log(maps[id_ltt]), col = "cadetblue4", lwd = 3, lty = 6)

        lines(ltt_times, means[1:length(id_ltt)], col = "darkseagreen1", lwd = 5, lty = 1)
        lines(ltt_times, means[1:length(id_ltt)], col = "darkseagreen4", lwd = 3, lty = 6)

        axis(1, cex.axis = axes_cex, lwd = 2, lwd.ticks = 2)
        axis(2, at = log(y_lab), lab = y_lab, cex.axis = axes_cex, lwd = 2, lwd.ticks = 2)

    """)
end

function plot_RTT(co::CladsOutput ;nplot = 50, miny = -1, maxy = 1, alpha_col = 0.05, burn = 0., axes_cex = 1., lwd = 1., lab_cex = 1., lay = 5,
	force_ylim = Array{Float64,1}(undef,0))

	times = co.time_points[2:end]
	RTT_map = co.RTT_map
    N = length(co.rtt_chains[1])
    t = times#[1:length(tr)]
	ini_plot = Int64(2+ floor(burn * N))
    plot_id = Array{Int64,1}(floor.(range(ini_plot, length=nplot, stop=N)))
    mean_mr = zeros(length(times))
    a = Array{Float64,1}(undef,0)
    map_mr = [deepcopy(a) for i in 1:length(times)]
	range_RTT = [minimum(minimum(minimum(co.rtt_chains)[2:end]));maximum(maximum(maximum(co.rtt_chains)[2:end]))]
    @rput t
    @rput miny
    @rput maxy
    @rput alpha_col
	@rput lab_cex
	@rput axes_cex
	@rput lay
	@rput RTT_map
	@rput range_RTT
	@rput force_ylim

    reval("""
		ylim = log(range(c(range_RTT, force_ylim)))#+c(miny,maxy)
        library(scales)
        plot(100000, axes = F, xlim = range(t), ylim = ylim,
        xlab = "time", ylab = "mean rate", cex.lab = lab_cex)
    """)
    n = 0
    colors = ["deepskyblue2" ; "deepskyblue2" ; "deepskyblue2"]

    #colors = ["deepskyblue2" ; "orange2" ; "orangered3"]
    for i in 1:3
        color = colors[i]
        @rput color
        for k in ini_plot:N
            y = co.rtt_chains[i][k]
            for j in 1:length(y)
                push!(map_mr[j],log(y[j]))
            end
			if k ∈ plot_id
            	@rput y
            	reval("""lines(t, log(y), col = alpha(color, alpha = alpha_col), lwd = 2)""")
			end
            mean_mr .+= y
            n += 1
        end
    end

    mean_mr ./= n
    @rput mean_mr
    @rput map_mr
    reval("""
        maps = log(RTT_map) #sapply(map_mr, function(x){D=density(x); return(D[[1]][which.max(D[[2]])])})
		#quant = sapply(map_mr, function(x){quantile(x, probs = c(0.05,0.95))})
		quant = sapply(map_mr, function(x){quantile(x, probs = c(0.05,0.95))})

    """)
    reval("""
		lwd = 2
        lines(t, quant[1,], col = "deepskyblue3", lwd = 2)
        lines(t, quant[2,], col = "deepskyblue3", lwd = 2)
        lines(t, maps, col = "darkseagreen1", lwd = 5, lty = 1)
        lines(t, maps, col = "darkseagreen4", lwd = 3, lty = 6)
        axis(1, cex.axis = axes_cex, lwd = lwd, lwd.ticks = lwd)
		labels_a2 = unique(signif(exp(seq(ylim[1], ylim[2], length.out = lay)), digits = 1))
		axis(2, at = log(labels_a2), labels = labels_a2, cex.axis = axes_cex, lwd = lwd, lwd.ticks = lwd)
    """)

end

function plot_density(co::CladsOutput, id_par::Int64; burn = 0.25)
	n_par = length(co.chains[1])
	n_tips = n_extant_tips(co.tree)
	n_edges = (n_tips - 1)*2
	nit = length(co.chains[1][1])
	init = Int64(ceil(burn * nit))
	color = "orange"

	if id_par >= 5 + n_edges + n_tips
		id_par += 2
	end

	if id_par <= n_par
		if id_par == 1				#σ
			name = "sigma"
			map = co.σ_map
		elseif id_par == 2			#α
			name = "alpha"
			map = co.α_map
		elseif id_par == 3			#ε
			name = "epsilon"
			map = co.ε_map
		elseif id_par == 4			#λ_0
			name = "lambda_0"
			map = co.λ0_map
		elseif id_par < 5 + n_edges #λ_i
			i = id_par - 4
			name = "lambda_$i"
			map = co.λi_map[i]
		elseif id_par < 5 + n_edges + n_tips						#λ_tip_i
			i = id_par - 4 - n_edges
			tip_name = tip_labels(co.tree)[i]
			name = "lambda_tip $i"
			if length(tip_name)>0
				name = "lambda_tip $tip_name"
			end
			map = co.λtip_map[i]
		else
			color = "#3498db"
			i = id_par - (6 + n_edges + n_tips)
			name = "number of lineages (time $(co.time_points[i]))"
			map = co.DTT_mean[i]
		end

		chain = [co.chains[i][id_par][j] for i in 1:length(co.chains) for j in init:length(co.chains[1][1])]
	else

		i = id_par - n_par
		name = "mean rate (time $(co.time_points[i]))"
		map = co.RTT_map[i]

		chain = [co.rtt_chains[i][j][id_par - n_par] for i in 1:length(co.rtt_chains) for j in init:length(co.rtt_chains[1])]
	end
	@rput chain
	@rput name
	@rput map
	@rput id_par
	@rput n_par
	@rput color
	reval("""
		if (id_par < 5 ){
			d = density(chain)
			plot(d\$x, d\$y, type = 'l', lwd = 3, xlab = name, ylab = "density")
		}else{
			d = density(log(chain))
			plot(exp(d\$x), d\$y, type = 'l', lwd = 3, xlab = name, ylab = "density", log = 'x')
		}
		abline(v = map, lwd = 3, col = color)

		color_q = adjustcolor("coral3", alpha.f = 0.1)
		quant = quantile(chain, c(0.025,0.975))
		polygon(quant[c(1,2,2,1)], c(-1,-1,100000,100000), col = color_q, border = NA)

	""")
end

function plot_chain(co::CladsOutput, id_par::Int64; burn = 0.25)
	n_par = length(co.chains[1])
	n_tips = n_extant_tips(co.tree)
	n_edges = (n_tips - 1)*2
	nit = length(co.chains[1][1])
	init = Int64(ceil(burn * nit))
	if id_par >= 5 + n_edges + n_tips
		id_par += 2
	end

	if id_par <= n_par
		if id_par == 1				#σ
			name = "sigma"
			map = co.σ_map
		elseif id_par == 2			#α
			name = "alpha"
			map = co.α_map
		elseif id_par == 3			#ε
			name = "epsilon"
			map = co.ε_map
		elseif id_par == 4			#λ_0
			name = "lambda_0"
			map = co.λ0_map
		elseif id_par < 5 + n_edges #λ_i
			i = id_par - 4
			name = "lambda_$i"
			map = co.λi_map[i]
		elseif id_par < 5 + n_edges + n_tips						#λ_tip_i
			i = id_par - 4 - n_edges
			tip_name = tip_labels(co.tree)[i]
			name = "lambda_tip $i"
			if length(tip_name)>0
				name = "lambda_tip $tip_name"
			end
			map = co.λtip_map[i]
		else
			i = id_par - (6 + n_edges + n_tips)
			name = "number of lineages (time $(co.time_points[i]))"
			map = co.DTT_mean[i]
		end

		chain = [co.chains[i][id_par][init:length(co.chains[1][1])] for i in 1:length(co.chains)]
	else

		i = id_par - n_par
		name = "mean rate (time $(co.time_points[i]))"
		map = co.RTT_map[i]

		chain = Array{Array{Float64,1},1}(undef,0)
		for i in 1:length(co.rtt_chains)
			push!(chain, [co.rtt_chains[i][j][id_par - n_par] for j in init:length(co.rtt_chains[1])])
		end
	end
	@rput chain
	@rput name
	@rput map
	@rput id_par
	@rput n_par

	reval("""
		require("coda")
		traceplot(mcmc.list(lapply(chain, mcmc)), ylab = name)

	""")
end

function convert_id(co::CladsOutput, id::String)
	n_par = length(co.chains[1])
	n_tips = n_extant_tips(co.tree)
	n_edges = (n_tips - 1)*2

	if occursin("σ",id) | occursin("sig",id)
		par_id = 1
	elseif occursin("α",id) | occursin("al",id)
		par_id = 2
	elseif occursin("ε",id) | occursin("ϵ",id) | occursin("eps",id) | occursin("turn",id)
		par_id = 3
	elseif occursin("λ0",id) | occursin("lambda0",id) | occursin("lambda_0",id) | occursin("lamb_0",id)
		par_id = 4
	elseif occursin(r"λ*tip",id) | occursin(r"lamb*tip",id)
		number = 0
		tip_id=split(id,"tip_")[end]
		re = r"^[+-]?([0-9]+([.][0-9]*)?|[.][0-9]+)$";
		if occursin(re, tip_id)
			number=parse(Int, tip_id)
		else
			dist = Inf
			tip_labs = tip_labels(co.tree)
			for itl in 1:length(tip_labs)
				d = Levenshtein()(tip_id, tip_labs[itl])
				if d < dist
					dist = d
					number = itl
					if d == 0
						break
					end
				end
			end
		end
		par_id = 4 + n_edges + number
	elseif occursin("λ",id) | occursin("lamb",id)
		number_as_string=split(id,"_")[end]
		number=parse(Int, number_as_string)
		par_id = 4 + number
	elseif occursin("lin",id) | occursin("div",id)
		number_as_string=split(id,"_")[end]
		number=parse(Int, number_as_string)
		par_id = 4 + n_edges + n_tips + number
	elseif occursin("rat",id)
		number_as_string=split(id,"_")[end]
		number=parse(Int, number_as_string)
		par_id = n_par - 2 + number
	else
		par_id = 1
	end

	return par_id
end

function plot_density(co::CladsOutput, id_par::String; burn = 0.25)
	id = convert_id(co, id_par)
	plot_density(co, id, burn = burn)
end

function plot_chain(co::CladsOutput, id_par::String; burn = 0.25)
	id = convert_id(co, id_par)
	plot_chain(co, id, burn = burn)
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
# PROBABILISTIC CIRCUITS LIBRARY ROOT

module ProbabilisticCircuits

# USE EXTERNAL MODULES

using Reexport

using LogicCircuits
# only reexport selectively from LogicCircuits
export pos_literals, neg_literals

include("Utils/Utils.jl")
@reexport using .Utils

include("abstract_prob_nodes.jl")
include("shared_prob_nodes.jl")
include("plain_prob_nodes.jl")
include("structured_prob_nodes.jl")
include("logistic_nodes.jl")
include("param_bit_circuit.jl")
include("param_bit_circuit_pair.jl")
include("parameters.jl")

include("queries/likelihood.jl")
include("queries/marginal_flow.jl")
include("queries/map.jl")
include("queries/sample.jl")
include("queries/pr_constraint.jl")
include("queries/information.jl")
include("queries/expectation_rec.jl")
include("queries/expectation_graph.jl")
include("queries/expectation_bit.jl")

include("Logistic/Logistic.jl")
@reexport using .Logistic

include("mixtures/em.jl")

include("structurelearner/chow_liu_tree.jl")
include("structurelearner/init.jl")
include("structurelearner/heuristics.jl")
include("structurelearner/learner.jl")
include("structurelearner/vtree_learner.jl")

include("LoadSave/LoadSave.jl")
@reexport using .LoadSave

end


text:
# Piracy that used to be requires, allowing Tracker.jl to be specialized for SciML

function RecursiveArrayTools.recursivecopy!(b::AbstractArray{T,N}, a::AbstractArray{T2,N}) where {T<:Tracker.TrackedArray,T2<:Tracker.TrackedArray,N}
    @inbounds for i in eachindex(a)
        b[i] = copy(a[i])
    end
end

DiffEqBase.value(x::Type{Tracker.TrackedReal{T}}) where {T} = T
DiffEqBase.value(x::Type{Tracker.TrackedArray{T,N,A}}) where {T,N,A} = Array{T,N}
DiffEqBase.value(x::Tracker.TrackedReal) = x.data
DiffEqBase.value(x::Tracker.TrackedArray) = x.data

DiffEqBase.promote_u0(u0::Tracker.TrackedArray, p::Tracker.TrackedArray, t0) = u0
DiffEqBase.promote_u0(u0::AbstractArray{<:Tracker.TrackedReal}, p::Tracker.TrackedArray, t0) = u0
DiffEqBase.promote_u0(u0::Tracker.TrackedArray, p::AbstractArray{<:Tracker.TrackedReal}, t0) = u0
DiffEqBase.promote_u0(u0::AbstractArray{<:Tracker.TrackedReal}, p::AbstractArray{<:Tracker.TrackedReal}, t0) = u0
DiffEqBase.promote_u0(u0, p::Tracker.TrackedArray, t0) = Tracker.track(u0)
DiffEqBase.promote_u0(u0, p::AbstractArray{<:Tracker.TrackedReal}, t0) = eltype(p).(u0)

@inline DiffEqBase.fastpow(x::Tracker.TrackedReal, y::Tracker.TrackedReal) = x^y
@inline Base.any(f::Function, x::Tracker.TrackedArray) = any(f, Tracker.data(x))

# Support adaptive with non-tracked time
@inline function DiffEqBase.ODE_DEFAULT_NORM(u::Tracker.TrackedArray, t) where {N}
    sqrt(sum(abs2, DiffEqBase.value(u)) / length(u))
end
@inline function DiffEqBase.ODE_DEFAULT_NORM(u::AbstractArray{<:Tracker.TrackedReal,N}, t) where {N}
    sqrt(sum(x -> DiffEqBase.ODE_DEFAULT_NORM(x[1], x[2]), zip((DiffEqBase.value(x) for x in u), Iterators.repeated(t))) / length(u))
end
@inline function DiffEqBase.ODE_DEFAULT_NORM(u::Array{<:Tracker.TrackedReal,N}, t) where {N}
    sqrt(sum(x -> DiffEqBase.ODE_DEFAULT_NORM(x[1], x[2]), zip((DiffEqBase.value(x) for x in u), Iterators.repeated(t))) / length(u))
end
@inline DiffEqBase.ODE_DEFAULT_NORM(u::Tracker.TrackedReal, t) = abs(DiffEqBase.value(u))

# Support TrackedReal time, don't drop tracking on the adaptivity there
@inline function DiffEqBase.ODE_DEFAULT_NORM(u::Tracker.TrackedArray, t::Tracker.TrackedReal) where {N}
    sqrt(sum(abs2, u) / length(u))
end
@inline function DiffEqBase.ODE_DEFAULT_NORM(u::AbstractArray{<:Tracker.TrackedReal,N}, t::Tracker.TrackedReal) where {N}
    sqrt(sum(x -> DiffEqBase.ODE_DEFAULT_NORM(x[1], x[2]), zip(u, Iterators.repeated(t))) / length(u))
end
@inline function DiffEqBase.ODE_DEFAULT_NORM(u::Array{<:Tracker.TrackedReal,N}, t::Tracker.TrackedReal) where {N}
    sqrt(sum(x -> DiffEqBase.ODE_DEFAULT_NORM(x[1], x[2]), zip(u, Iterators.repeated(t))) / length(u))
end
@inline DiffEqBase.ODE_DEFAULT_NORM(u::Tracker.TrackedReal, t::Tracker.TrackedReal) = abs(u)

function DiffEqBase.solve_up(prob::DiffEqBase.DEProblem, sensealg::Union{DiffEqBase.AbstractSensitivityAlgorithm,Nothing}, u0::Tracker.TrackedArray, p::Tracker.TrackedArray, args...; kwargs...)
    Tracker.track(solve_up, prob, sensealg, u0, p, args...; kwargs...)
end

function DiffEqBase.solve_up(prob::DiffEqBase.DEProblem, sensealg::Union{DiffEqBase.AbstractSensitivityAlgorithm,Nothing}, u0::Tracker.TrackedArray, p, args...; kwargs...)
    Tracker.track(solve_up, prob, sensealg, u0, p, args...; kwargs...)
end

function DiffEqBase.solve_up(prob::DiffEqBase.DEProblem, sensealg::Union{DiffEqBase.AbstractSensitivityAlgorithm,Nothing}, u0, p::Tracker.TrackedArray, args...; kwargs...)
    Tracker.track(solve_up, prob, sensealg, u0, p, args...; kwargs...)
end

Tracker.@grad function DiffEqBase.solve_up(prob, sensealg::Union{Nothing,DiffEqBase.AbstractSensitivityAlgorithm},
    u0, p, args...;
    kwargs...)
    _solve_adjoint(prob, sensealg, Tracker.data(u0), Tracker.data(p), args...; kwargs...)
end


text:
# Deterministic Value Function Iteration
# Greg Kaplan 2017
# Translated by Tom Sweeney Dec 2020

using Random, Interpolations, Plots

# PARAMETERS

## preferences
risk_aver = 2
beta = 0.95

## returns
r = 0.03
R = 1+r

## income
y = 1

## asset grids
na = 1000
amax = 20
borrow_lim = 0
agrid_par = 1 # 1 for linear, 0 for L-shaped

## computation
max_iter = 1000
tol_iter = 1.0e-6
Nsim = 100
Tsim = 500

# OPTIONS
Display = 1
DoSimulate = 1
MakePlots = 1

# DRAW RANDOM NUMBERS
Random.seed!(2020)
arand = rand(Nsim)

# SET UP GRIDS

## assets
agrid = range(0,1,length=na)
agrid = agrid.^(1 ./ agrid_par)
agrid = borrow_lim .+ (amax.-borrow_lim).*agrid

# UTILITY FUNCTION

if risk_aver==1
    u(c) = log.(c)
else
    u(c) = (c.^(1-risk_aver).-1)./(1-risk_aver)
end

u1(c) = c.^(-risk_aver)

# INITIALIZE VALUE FUNCTION

Vguess = u(r.*agrid.+y)./(1-beta)

# ITERATE ON VALUE FUNCTION

V = copy(Vguess)

Vdiff = 1
iter = 0

while iter <= max_iter && Vdiff > tol_iter
    iter = iter + 1
    Vlast = copy(V)
    V = zeros(na)
    global sav = zeros(na)
    global savind = zeros(Int,na)
    global con = zeros(na)

    ## loop over assets
    for ia = 1:na

        cash = R.*agrid[ia] + y
        Vchoice = u(max.(cash.-agrid,1.0e-10)) + beta.*Vlast
        V[ia] = maximum(Vchoice)
        savind[ia] = argmax(Vchoice)[1]
        sav[ia] = agrid[savind[ia]]
        con[ia] = cash .- sav[ia]
    end

    Vdiff = maximum(abs.(V-Vlast))
    if Display>=1
        println("Iteration no. " * string(iter), " max val fn diff is " * string(Vdiff))
    end
end

# SIMULATE
if DoSimulate==1
    yindsim = zeros(Int,Nsim,Tsim)
    aindsim = zeros(Int,Nsim,Tsim)

    ## initial assets: uniform on [borrow_lim, amax]
    ainitial = borrow_lim .+ arand.*(amax-borrow_lim)

    ## allocate to nearest point on agrid;
    aindsim[:,1] = interpolate((agrid,), 1:na, Gridded(Constant())).(ainitial)

    ## loop over time periods
    for it = 1:Tsim
        if Display >= 1 && mod(it,100)==0
            println(" Simulating, time period " * string(it))
        end
        ## asset choice
        if it<Tsim
            aindsim[:,it+1] = savind[aindsim[:,it]]
        end
    end

    ## assign actual asset and income values;
    asim = agrid[aindsim]
    csim = R.*asim[:,1:Tsim-1] .+ y .- asim[:,2:Tsim]
end

# MAKE PLOTS
if MakePlots==1

    ## consumption policy function
    p1 = plot(agrid, con, xlims=(0,amax), title="Consumption", color=:blue, legend=false)
    display(p1)

    ## savings policy function
    p2 = plot(agrid, sav.-agrid, xlims=(0,amax), title="Savings", color=:blue, legend=false)
    plot!(agrid, zeros(na,1), color=:black, lw=0.5)
    display(p2)

    ## nice zoom
    xlimits = (0,1)
    xlimind = trues(na)
    if minimum(agrid) < xlimits[1]
        xlimind = xlimind .& (agrid.>=maximum(agrid[agrid<xlimits[1]]))
    elseif minimum(agrid) > xlimits[2]
        xlimind .= 0
    end
    if maximum(agrid) > xlimits[2]
        xlimind = xlimind .& (agrid.<=minimum(agrid[agrid.>xlimits[2]]))
    elseif maximum(agrid) < xlimits[1]
        xlimind .= 0
    end

    ## consumption policy function: zoomed in
    p3 = plot(agrid[xlimind], con[xlimind], xlims=xlimits, title="Consumption: Zoomed", marker=:circle, color=:blue, linewidth=2, legend=false)
    plot!(show=true)
    display(p3)

    ## savings policy function: zoomed in
    p4 = plot(agrid[xlimind], sav[xlimind].-agrid[xlimind], xlims=xlimits, title="Savings: Zoomed (a'-a)", marker=:circle, color=:blue, linewidth=2, legend=false)
    plot!(agrid, zeros(na), color=:black, lw=0.5)
    display(p4)

    ## asset dynamics distribution
    p5 = plot(repeat((1:Tsim),1,Nsim), asim', title="Asset Dynamics", legend=false)
    display(p5)

    ## consumption dynamics distribution
    p6 = plot(repeat((1:Tsim-1),1,Nsim), csim', title="Consumption Dynamics", legend=false)
    display(p6)

end


text:
function solve!(model::MoominModel, optimizer; enumerate=1, stoichiometry=true, printLevel=1, timeLimit=1000)
  printLevel == 0 || @info "Creating MILP..."
  MILP = createMILP(model, optimizer, stoichiometry=stoichiometry, timeLimit=timeLimit)
  printLevel < 2 && set_silent(MILP)
  cont = true
  counter = 1
  outputColours = []
  while cont
    printLevel == 0 || @info "Solving MILP #$counter..."
    optimize!(MILP)
    solFound = termination_status(MILP) == MOI.OPTIMAL
    timeLimitReached = termination_status(MILP) == MOI.TIME_LIMIT
    (solFound & (printLevel > 0)) && @info "Found solution!"
    (!solFound & (printLevel > 0) & (counter > 1) ) && @info "No optimum found."
    (!solFound & (counter == 1) ) && @warn "Couldn't solve MILP #1."
    timeLimitReached && @warn "Solver time limit reached."
    if solFound
      if isempty(outputColours)
        outputColours = interpretSolution(value.(MILP[:xPlus]),
                                          value.(MILP[:xMinus]),
                                          model.reactions.inputColours,
                                          model.reactions.reversible)
      else
        outputColours = [outputColours interpretSolution(value.(MILP[:xPlus]),
                                          value.(MILP[:xMinus]),
                                          model.reactions.inputColours,
                                          model.reactions.reversible)]
      end
      updateMILP!(MILP, outputColours)
    end
    cont = solFound & (counter<enumerate)
    counter = counter+1
  end

  fillOutputs!(model, outputColours)
end

function interpretSolution(xPlus, xMinus, inputColours, reversible)
  all(.!((xPlus .> 0.5) .& (xMinus .> 0.5))) || error("Something is wrong with a solution: `x+` and `x-` should be mutually exclusive.")

  outputColours = zeros(length(xPlus))
  outputColours[xPlus .> 0.5] .= 1
  outputColours[xMinus .> 0.5] .= -1
  outputColours[(xPlus .> 0.5) .& (inputColours .== -1)] .= -2
  outputColours[(xMinus .> 0.5) .& (inputColours .== 1)] .= 2
  outputColours[((xPlus .> 0.5) .| (xMinus .> 0.5)) .& (inputColours .==0) .& reversible] .= 6

  return outputColours
end

function fillOutputs!(model::MoominModel, outputColours)
  if isempty(outputColours)
    return
  end
  model.reactions.outputColours = convert.(Int64, outputColours)
  model.reactions.outputFrequency = round.(sum(outputColours .!= 0, dims=2) ./ size(outputColours, 2), digits=3)[:]
  combined = zeros(size(outputColours, 1))
  for row in 1:size(outputColours, 1)
    nonZ = findall(outputColours[row, :] .!= 0)
    if !isempty(nonZ)
      if all(outputColours[row, nonZ[1]] .== outputColours[row, nonZ])
        combined[row] = outputColours[row, nonZ[1]]
      else
        combined[row] = 6
      end
    end
  end
  model.reactions.combinedOutput = convert.(Int64, combined)
end

function createMILP(model::MoominModel, optimizer; stoichiometry=true, timeLimit=1000)
  MILP = JuMP.Model(optimizer)
  if solver_name(MILP) == "CPLEX"
    set_optimizer_attribute(MILP, "CPXPARAM_TimeLimit", timeLimit)
  elseif solver_name(MILP) == "Gurobi"
    set_optimizer_attribute(MILP, "TimeLimit", timeLimit)
  end

  if stoichiometry
    MILPstoich!(MILP, model)
  else
    MILPtopo!(MILP, model)
  end

  return MILP
end

function MILPstoich!(MILP, model::MoominModel)
  epsilon = 1
  (m, n) = size(model.reactions.S)

  # impose a priori colours
  lb = fill(-100., n)
  lb[.!model.reactions.reversible .& (model.reactions.inputColours.==1)] .= 0
  ub = fill(100., n)
  ub[.!model.reactions.reversible .& (model.reactions.inputColours.==-1)] .= 0

  @variable(MILP, v[i=1:n], lower_bound=lb[i], upper_bound=ub[i])
  @variable(MILP, xPlus[1:n], Bin)
  @variable(MILP, xMinus[1:n], Bin)

  # stoichiometry
  @constraint(MILP, model.reactions.S*v .== zeros(m))
  # x+=1 -> v>=epsilon
  @constraint(MILP, v .+ xPlus.*(lb .- epsilon) .>= lb)
  # x+=0 -> v<=0
  @constraint(MILP, v .- xPlus.*ub .<= zeros(n))
  # x-=1 -> v<=-epsilon
  @constraint(MILP, v .+ xMinus.*(ub .+ epsilon) .<= ub)
  # x-=0 -> v>=0
  @constraint(MILP, v .+ xMinus.*(-lb) .>= zeros(n))

  @objective(MILP, Max, model.reactions.weights' * xPlus + model.reactions.weights' * xMinus)
end

function MILPtopo!(MILP, model::MoominModel)
  (m, n) = size(model.reactions.S)

  @variable(MILP, y[1:m], Bin)
  @variable(MILP, xPlus[1:n], Bin)
  @variable(MILP, xMinus[1:n], Bin)

  # impose a priori colours
  notRed = (.!model.reactions.reversible) .& (model.reactions.inputColours.==-1)
  @constraint(MILP, xPlus[notRed] .== 0)
  notBlue = (.!model.reactions.reversible) .& (model.reactions.inputColours.==1)
  @constraint(MILP, xMinus[notBlue] .== 0)

  # x+ and x- cannot be 1 at the same time
  @constraint(MILP, xPlus .+ xMinus .<= ones(n))
  # if a connected arc is included, a node is included
  @constraint(MILP, (model.reactions.S .!= 0) * (xPlus.+xMinus) .<= sum(model.reactions.S .!= 0, dims=2).*y)
  # if a node is included, it has to have an outgoing arc
  @constraint(MILP, (model.reactions.S .< 0) * xPlus .+ (model.reactions.S .> 0) * xMinus .>= y)
  # if a node is included, it has to have an incoming arc
  @constraint(MILP, (model.reactions.S .> 0) * xPlus .+ (model.reactions.S .< 0) * xMinus .>= y)

  @objective(MILP, Max, model.reactions.weights' * xPlus + model.reactions.weights' * xMinus)
end

function updateMILP!(MILP, outputColours)
  if isempty(outputColours)
    return
  end
  if isnothing(constraint_by_name(MILP, "optimality"))
    @constraint(MILP, optimality, objective_function(MILP) == objective_value(MILP))
  end
  prevSol = outputColours[:, end] .!= 0
  @constraint(MILP, 2*(MILP[:xPlus]'*prevSol + MILP[:xMinus]'*prevSol)
                      - sum(MILP[:xPlus] .+ MILP[:xMinus])
                      <= sum(prevSol) - 1)
end

max_stars_repo_name: htpusa/moominJL

max_issues_repo_name: htpusa/moominJL

max_forks_repo_name: htpusa/moominJL


text:
module Dimensions
using ...UnitfulHartree
using ...Dispatch: Dimensioned, DFTAxisArray, DFTArray, Dimensioned
using AxisArrays, Unitful

const UH = UnitfulHartree

macro lintpragma(s) end
@lintpragma("Ignore use of undeclared variable T")
@lintpragma("Ignore use of undeclared variable U")
@lintpragma("Ignore use of undeclared variable N")
@lintpragma("Ignore use of undeclared variable A")
@lintpragma("Ignore use of undeclared variable Eₕ")

for (name, abbr) in Dimensioned
    DAA = Symbol("$(name)AxisArray")
    DDA = Symbol("$(name)DenseArray")
    DA = Symbol("$(name)Array")
    abbr = ifelse(abbr == :σ, :σₑ, abbr)

    @eval begin
        const $name = Quantity{T, typeof(dimension(UH.$abbr)), U} where U where T
        const $DAA = DFTAxisArray{$name{T, U}, N, A} where U where A where N where T
        const $DDA = DenseArray{$name{T, U}, N} where U where N where T
        const $DA = DFTArray{$name{T, U}, N} where U where N where T
    end
end

""" Dimenional dispatch types for scalars """
module Scalars
    using ...Dispatch: Dimensioned, Dimensions
    for (name, abbr) in Dimensioned
        @eval begin
            const $name = Dimensions.$name
            const $abbr = $name
        end
    end
    const All = @eval Union{$((Expr(:curly, n, :T) for (n, a) in Dimensioned)...)} where T
    const ϵ = Eₕ
end

""" Dimensional dispatch types for Arrays """
module Arrays
    using ...Dispatch: Dimensioned, Dimensions
    for (name, abbr) in Dimensioned
        @eval begin
            const $name = Dimensions.$(Symbol("$(name)Array"))
            const $abbr = $name
        end
    end
    const All = @eval Union{
            $((Expr(:curly, n, :T, :N) for (n, a) in Dimensioned)...)} where {T, N}
    const ϵ = Eₕ
end

""" Dimenional dispatch types for AxisArrays """
module AxisArrays
    using ...Dispatch: Dimensioned, Dimensions
    for (name, abbr) in Dimensioned
        @eval begin
            const $name = Dimensions.$(Symbol("$(name)AxisArray"))
            const $abbr = $name
        end
    end
    const All = @eval Union{
            $((Expr(:curly, n, :T, :N, :A) for (n, a) in Dimensioned)...)} where {T, N, A}
    const ϵ = Eₕ
end

""" Dimenional dispatch types for DenseArrays """
module DenseArrays
    using ...Dispatch: Dimensioned, Dimensions
    for (name, abbr) in Dimensioned
        @eval begin
            const $name = Dimensions.$(Symbol("$(name)DenseArray"))
            const $abbr = $name
        end
    end
    const All = @eval Union{
                    $((Expr(:curly, n, :T, :N) for (n, a) in Dimensioned)...)} where {T, N}
    const ϵ = Eₕ
end

end

max_stars_repo_name: mdavezac/DFTShims

max_issues_repo_name: mdavezac/DFTShims

max_forks_repo_name: mdavezac/DFTShims


text:
result = [
    2.32551,
    -0.155293,
    0.469214,
    -0.0870457,
    0.271364,
    -0.222398,
    1.01423,
    -0.179009,
    1.01724,
    -0.152972,
    0.216039,
    -0.251637,
    1.05034,
    -0.231289,
    0.245388,
    -0.231289,
    1.03371,
    -0.258872,
    1.03256,
    -0.262519,
    1.04748,
    -0.259143,
    1.0535,
    -0.266484,
    1.04711,
    -0.267177,
    1.02131,
    -0.280381,
]

include(joinpath(BASE_DIR, "test", "data_5bus_pu.jl"))
include(joinpath(BASE_DIR, "test", "data_14bus_pu.jl"))

nodes_5 = nodes5()
nodes_14 = nodes14()
c_sys14() = System(
    nodes_14,
    thermal_generators14(nodes_14),
    loads14(nodes_14),
    branches14(nodes_14),
    nothing,
    100.0,
    nothing,
    nothing,
);

c_sys5_re() = System(
    nodes_5,
    vcat(thermal_generators5(nodes_5), renewable_generators5(nodes_5)),
    loads5(nodes_5),
    nothing,
    nothing,
    100.0,
    nothing,
    nothing,
)

@testset begin
    using NLsolve
    # This is a negative test. The data passed for sys5_re is known to be infeasible.
    @test_logs(
        (:error, "The powerflow solver returned convergence = false"),
        match_mode = :any,
        @test !solve_powerflow!(c_sys5_re(), nlsolve)
    )
    @test solve_powerflow!(c_sys14(), nlsolve, method = :newton)
end


text:
const N = 40
const F = 8.0

const Tmax = 10.0
const dT = 0.05
const nT = convert(Int32, Tmax/dT)

function lorenz96(x::Vector{Float64}, _, _)
    (circshift(x, -1) - circshift(x, 2)) .* circshift(x, 1) - x .+ F
end

max_stars_repo_name: Typas/Data-Assimilation-Project

max_issues_repo_name: Typas/Data-Assimilation-Project

max_forks_repo_name: Typas/Data-Assimilation-Project


text:
let
    pomdp = TigerPOMDP()

    mdp = UnderlyingMDP(pomdp)

    @test n_states(mdp) == n_states(pomdp)
    @test states(mdp) == states(pomdp)
    s_mdp = rand(MersenneTwister(1), initialstate_distribution(mdp))
    s_pomdp = rand(MersenneTwister(1), initialstate_distribution(pomdp))

    @test s_mdp == s_pomdp

    solver = ValueIterationSolver(max_iterations = 100)
    mdp_policy = solve(solver, mdp)
    pomdp_policy = solve(solver, pomdp)
    @test mdp_policy.util == pomdp_policy.util

    actionindex(mdp, 1)
end