text:
"""
    InverseGaussianExpert(μ, λ)

PDF:

```math
f(x; \\mu, \\lambda) = \\sqrt{\\frac{\\lambda}{2\\pi x^3}}
\\exp\\left(\\frac{-\\lambda(x-\\mu)^2}{2\\mu^2x}\\right),
\\quad x > 0
```

See also: [Inverse Gaussian Distribution](https://en.wikipedia.org/wiki/Inverse_Gaussian_distribution) (Wikipedia)

"""
struct InverseGaussianExpert{T<:Real} <: NonZIContinuousExpert
    μ::T
    λ::T
    InverseGaussianExpert{T}(µ::T, λ::T) where {T<:Real} = new{T}(µ, λ)
end

function InverseGaussianExpert(μ::T, λ::T; check_args=true) where {T <: Real}
    check_args && @check_args(InverseGaussianExpert, μ >= zero(μ) && λ > zero(λ))
    return InverseGaussianExpert{T}(μ, λ)
end

## Outer constructors
InverseGaussianExpert(μ::Real, λ::Real) = InverseGaussianExpert(promote(μ, λ)...)
InverseGaussianExpert(μ::Integer, λ::Integer) = InverseGaussianExpert(float(μ), float(λ))
InverseGaussianExpert() = InverseGaussianExpert(1.0, 1.0)

## Conversion
function convert(::Type{InverseGaussianExpert{T}}, μ::S, λ::S) where {T <: Real, S <: Real}
    InverseGaussianExpert(T(μ), T(λ))
end
function convert(::Type{InverseGaussianExpert{T}}, d::InverseGaussianExpert{S}) where {T <: Real, S <: Real}
    InverseGaussianExpert(T(d.μ), T(d.λ), check_args=false)
end
copy(d::InverseGaussianExpert) = InverseGaussianExpert(d.μ, d.λ, check_args=false)

## Loglikelihood of Expert
logpdf(d::InverseGaussianExpert, x...) = isinf(x...) ? -Inf : Distributions.logpdf.(Distributions.InverseGaussian(d.μ, d.λ), x...)
pdf(d::InverseGaussianExpert, x...) = isinf(x...) ? 0.0 : Distributions.pdf.(Distributions.InverseGaussian(d.μ, d.λ), x...)
logcdf(d::InverseGaussianExpert, x...) = isinf(x...) ? 0.0 : Distributions.logcdf.(Distributions.InverseGaussian(d.μ, d.λ), x...)
cdf(d::InverseGaussianExpert, x...) = isinf(x...) ? 1.0 : Distributions.cdf.(Distributions.InverseGaussian(d.μ, d.λ), x...)

## expert_ll, etc
expert_ll_exact(d::InverseGaussianExpert, x::Real) = LRMoE.logpdf(d, x)
function expert_ll(d::InverseGaussianExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_ll = (yl == yu) ? logpdf.(d, yl) : logcdf.(d, yu) + log1mexp.(logcdf.(d, yl) - logcdf.(d, yu))
    expert_ll = (tu == 0.) ? -Inf : expert_ll
    return expert_ll
end
function expert_tn(d::InverseGaussianExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn = (tl == tu) ? logpdf.(d, tl) : logcdf.(d, tu) + log1mexp.(logcdf.(d, tl) - logcdf.(d, tu))
    expert_tn = (tu == 0.) ? -Inf : expert_tn
    return expert_tn
end
function expert_tn_bar(d::InverseGaussianExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn_bar = (tl == tu) ? 0.0 : log1mexp.(logcdf.(d, tu) + log1mexp.(logcdf.(d, tl) - logcdf.(d, tu)))
    return expert_tn_bar
end

exposurize_expert(d::InverseGaussianExpert; exposure = 1) = d

## Parameters
params(d::InverseGaussianExpert) = (d.μ, d.λ)
function params_init(y, d::InverseGaussianExpert)
    pos_idx = (y .> 0.0)
    μ, σ2 = mean(y[pos_idx]), var(y[pos_idx])
    μ_init = μ
    λ_init = μ^3 / σ2
    if isnan(μ_init) || isnan(λ_init)
        return InverseGaussianExpert()
    else
        return InverseGaussianExpert(μ_init, λ_init)
    end
end

## KS stats for parameter initialization
function ks_distance(y, d::InverseGaussianExpert)
    p_zero = sum(y .== 0.0) / sum(y .>= 0.0)
    return max(abs(p_zero-0.0), (1-0.0)*HypothesisTests.ksstats(y[y .> 0.0], Distributions.InverseGaussian(d.μ, d.λ))[2])
end

## Simululation
sim_expert(d::InverseGaussianExpert) = Distributions.rand(Distributions.InverseGaussian(d.μ, d.λ), 1)[1]

## penalty
penalty_init(d::InverseGaussianExpert) = [1.0 Inf 1.0 Inf]
no_penalty_init(d::InverseGaussianExpert) = [1.0 Inf 1.0 Inf]
penalize(d::InverseGaussianExpert, p) = (p[1]-1)*log(d.μ) - d.μ/p[2] + (p[3]-1)*log(d.λ) - d.λ/p[4]

## statistics
mean(d::InverseGaussianExpert) = mean(Distributions.InverseGaussian(d.μ, d.λ))
var(d::InverseGaussianExpert) = var(Distributions.InverseGaussian(d.μ, d.λ))
quantile(d::InverseGaussianExpert, p) = quantile(Distributions.InverseGaussian(d.μ, d.λ), p)
function lev(d::InverseGaussianExpert, u)
    z = (u-d.μ)/d.μ
    y = (u+d.μ)/d.μ
    return u - d.μ*z*cdf.(Normal(0, 1), z*(d.λ/u)^0.5) - d.μ*y*exp(2*d.λ/d.μ)*cdf.(Normal(0, 1), -y*(d.λ/u)^0.5)
end
excess(d::InverseGaussianExpert, u) = mean(d) - lev(d, u)

## Misc functions for E-Step
function _int_y_func(d::InverseGaussianExpert, x)
    return (iszero(x) || isinf(x)) ? 0.0 : x * pdf.(d, x)
end

function _int_obs_Y_raw(d::InverseGaussianExpert, yl, yu)
    if yl == yu
        return yl * pdf.(d, yl)
    else
        return quadgk.(x -> _int_y_func(d, x), yl, yu, rtol = 1e-8)[1]
    end
end

function _int_lat_Y_raw(d::InverseGaussianExpert, tl, tu)
    return (tl==0 ? 0.0 : quadgk.(x -> _int_y_func(d, x), 0.0, tl, rtol = 1e-8)[1]) + (isinf(tu) ? 0.0 : quadgk.(x -> _int_y_func(d, x), tu, Inf, rtol = 1e-8)[1])
end

function _int_logy_func(d::InverseGaussianExpert, x)
    return (iszero(x) || isinf(x)) ? 0.0 : log(x) * pdf.(d, x)
end

function _int_obs_logY_raw(d::InverseGaussianExpert, yl, yu)
    if yl == yu
        return log(yl) * pdf.(d, yl)
    else
        return quadgk.(x -> _int_logy_func(d, x), yl, yu, rtol = 1e-8)[1]
    end
end

function _int_lat_logY_raw(d::InverseGaussianExpert, tl, tu)
    return (tl==0 ? 0.0 : quadgk.(x -> _int_logy_func(d, x), 0.0, tl, rtol = 1e-8)[1]) + (isinf(tu) ? 0.0 : quadgk.(x -> _int_logy_func(d, x), tu, Inf, rtol = 1e-8)[1])
end

function _int_invy_func(d::InverseGaussianExpert, x)
    return (iszero(x) || isinf(x)) ? 0.0 : 1/(x) * pdf.(d, x)
end

function _int_obs_invY_raw(d::InverseGaussianExpert, yl, yu)
    if yl == yu
        return 1/yl * pdf.(d, yl)
    else
        return quadgk.(x -> _int_invy_func(d, x), yl, yu, rtol = 1e-8)[1]
    end
end

function _int_lat_invY_raw(d::InverseGaussianExpert, tl, tu)
    return (tl==0 ? 0.0 : quadgk.(x -> _int_invy_func(d, x), 0.0, tl, rtol = 1e-8)[1]) + (isinf(tu) ? 0.0 : quadgk.(x -> _int_invy_func(d, x), tu, Inf, rtol = 1e-8)[1])
end

## EM: M-Step
function EM_M_expert(d::InverseGaussianExpert,
                     tl, yl, yu, tu,
                     exposure,
                     z_e_obs, z_e_lat, k_e;
                     penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf])

    expert_ll_pos = expert_ll.(d, tl, yl, yu, tu)
    expert_tn_bar_pos = expert_tn_bar.(d, tl, yl, yu, tu)

    # Further E-Step
    yl_yu_unique = unique_bounds(yl, yu)

    int_obs_Y_tmp = _int_obs_Y_raw.(d, yl_yu_unique[:,1], yl_yu_unique[:,2])
    # int_obs_logY_tmp = _int_obs_logY_raw.(d, yl_yu_unique[:,1], yl_yu_unique[:,2])
    int_obs_invY_tmp = _int_obs_invY_raw.(d, yl_yu_unique[:,1], yl_yu_unique[:,2])

    Y_e_obs = exp.(-expert_ll_pos) .* int_obs_Y_tmp[match_unique_bounds(hcat(vec(yl), vec(yu)), yl_yu_unique)]
    # logY_e_obs = exp.(-expert_ll_pos) .* int_obs_logY_tmp[match_unique_bounds(hcat(vec(yl), vec(yu)), yl_yu_unique)]
    invY_e_obs = exp.(-expert_ll_pos) .* int_obs_invY_tmp[match_unique_bounds(hcat(vec(yl), vec(yu)), yl_yu_unique)]

    nan2num(Y_e_obs, 0.0) # get rid of NaN
    # nan2num(logY_e_obs, 0.0) # get rid of NaN
    nan2num(invY_e_obs, 0.0) # get rid of NaN

    tl_tu_unique = unique_bounds(tl, tu)

    int_lat_Y_tmp = _int_lat_Y_raw.(d, tl_tu_unique[:,1], tl_tu_unique[:,2])
    # int_lat_logY_tmp = _int_lat_logY_raw.(d, tl_tu_unique[:,1], tl_tu_unique[:,2])
    int_lat_invY_tmp = _int_lat_invY_raw.(d, tl_tu_unique[:,1], tl_tu_unique[:,2])

    Y_e_lat = exp.(-expert_tn_bar_pos) .* int_lat_Y_tmp[match_unique_bounds(hcat(vec(tl), vec(tu)), tl_tu_unique)]
    # logY_e_lat = exp.(-expert_tn_bar_pos) .* int_lat_logY_tmp[match_unique_bounds(hcat(vec(tl), vec(tu)), tl_tu_unique)]
    invY_e_lat = exp.(-expert_tn_bar_pos) .* int_lat_invY_tmp[match_unique_bounds(hcat(vec(tl), vec(tu)), tl_tu_unique)]

    nan2num(Y_e_lat, 0.0) # get rid of NaN
    # nan2num(logY_e_lat, 0.0) # get rid of NaN
    nan2num(invY_e_lat, 0.0) # get rid of NaN

    # Update parameters
    pos_idx = (yu .!= 0.0)
    term_zkz = z_e_obs[pos_idx] .+ (z_e_lat[pos_idx] .* k_e[pos_idx])
    term_zkz_Y = (z_e_obs[pos_idx] .* Y_e_obs[pos_idx]) .+ (z_e_lat[pos_idx] .* k_e[pos_idx] .* Y_e_lat[pos_idx])
    # term_zkz_logY = (z_e_obs[pos_idx] .* logY_e_obs[pos_idx]) .+ (z_e_lat[pos_idx] .* k_e[pos_idx] .* logY_e_lat[pos_idx])
    term_zkz_invY = (z_e_obs[pos_idx] .* invY_e_obs[pos_idx]) .+ (z_e_lat[pos_idx] .* k_e[pos_idx] .* invY_e_lat[pos_idx])

    μ_new = sum(term_zkz_Y)[1] / sum(term_zkz)[1]
    λ_new = sum(term_zkz)[1] / ( sum(term_zkz_Y)[1]/(μ_new)^2 - 2*sum(term_zkz)[1]/μ_new + sum(term_zkz_invY)[1] )

    return InverseGaussianExpert(μ_new, λ_new)
end

## EM: M-Step, exact observations
function EM_M_expert_exact(d::InverseGaussianExpert,
                            ye, exposure,
                            z_e_obs;
                            penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf])

    # Further E-Step
    Y_e_obs = ye
    invY_e_obs = 1.0 ./ ye

    nan2num(Y_e_obs, 0.0) # get rid of NaN

    # Update parameters
    pos_idx = (ye .!= 0.0)
    term_zkz = z_e_obs[pos_idx]
    term_zkz_Y = z_e_obs[pos_idx] .* Y_e_obs[pos_idx]
    term_zkz_invY = z_e_obs[pos_idx] .* invY_e_obs[pos_idx]

    μ_new = sum(term_zkz_Y)[1] / sum(term_zkz)[1]
    λ_new = sum(term_zkz)[1] / ( sum(term_zkz_Y)[1]/(μ_new)^2 - 2*sum(term_zkz)[1]/μ_new + sum(term_zkz_invY)[1] )

    return InverseGaussianExpert(μ_new, λ_new)
end


text:
# TODO: better names for the functions

function updatecoeffs!{T<:Real}(coeff0, coeff1, coeff2, rhs, model, v, t::T, x,
                                a::T, Δτ::T, Δx::T)
    # TODO: should we just remove this dispatch and pass in ones(n-2)*a instead?
    # Updates
    # coeffn = values in linear system
    # rhs    = value function at previous timestep + f at current timestep

    # Input
    # model  = HJBOneDim object
    # v      = value function at previous timestep
    # t      = value of (forward) time
    # x      = vector of x-values
    # a      = policy-values on interior
    # Δτ     = time-step size
    # Δx     = spacial step length
    taux = Δτ/Δx
    htaux2 = 0.5*Δτ/Δx^2

    n = length(coeff0)
    for j = 2:n-1
        @inbounds begin
            bval = model.b(t,x[j],a)
            sval2 = model.σ(t,x[j],a)^2
            coeff1[j] = -(sval2*htaux2 + max(bval,0.)*taux)
            coeff2[j-1] = -(sval2*htaux2 - min(bval,0.)*taux)
            coeff0[j] = 1.-coeff1[j]-coeff2[j-1]
            rhs[j] = v[j] + Δτ*model.f(t,x[j],a)
        end
    end

end

function policytimestep(model::HJBOneDim,
                        v, avals, x, Δx, Δτ, ti::Int)
    taux = Δτ/Δx
    htaux2 = 0.5*Δτ/Δx^2
    t = (ti-1)*Δτ
    n = length(v)

    # TODO: redo this thing
    newind = ones(Int, n)
    @inbounds vnew = fill(-maxintfloat(typeof(v[1])), n)

    ind12 = zeros(Bool, length(vnew))

    coeff0 = ones(v)   # v_i
    coeff1 = zeros(n-1) # v_{i+1} # TODO: type stability
    coeff2 = zeros(n-1) # v_{i-1} # TODO: type stability
    rhs = zeros(v)
    # Dirichlet conditions
    if model.bcond[1] == true
        @inbounds rhs[1] = model.Dfun[1](t)
    end
    if model.bcond[2] == true
        @inbounds rhs[end] = model.Dfun[2](t)
    end

    for i = 1:length(avals)
        @inbounds begin
            a = avals[i]
            vold = vnew

            updatecoeffs!(coeff0, coeff1, coeff2, rhs, model, v, t, x, a, Δτ, Δx)
            Mat = spdiagm((coeff2, coeff0, coeff1), -1:1, n, n)
            # Move to sparse(I,J,K) instead?
            if model.bcond[1] == false
                bval = taux*model.b(t,x[1],a)
                sval2 = htaux2*model.σ(t,x[1],a)^2
                Mat[1,1:3] = [1.-(sval2-bval), -(bval-2*sval2), -sval2]
                rhs[1] = v[1] + Δτ*model.f(t,x[1],a)
            end

            if model.bcond[2] == false
                bval = taux*model.b(t,x[end],a)
                sval2 = htaux2*model.σ(t,x[end],a)^2
                Mat[end,end-2:end] = [-sval2, bval+2*sval2, 1.-(sval2+bval)]
                rhs[end] = v[end] + Δτ*model.f(t,x[end],a)
            end

            # TODO: Use Krylov solver for high-dimensional PDEs
            vnew = Mat\rhs

            ind12[:] = vold .> vnew
            vnew[ind12] = vold[ind12]
            newind[!ind12] = i
        end
    end

    @inbounds pol = avals[newind]

    return vnew, pol
end

function timeloopconstant(model::HJBOneDim, N::Int,
                          Δτ, vinit, avals, x, Δx)
    # Pass v and pol by reference?
    K = length(x)
    v = zeros(K, N+1)
    pol = zeros(K, N)

    @inbounds v[:,end] = vinit # We use forward time t instead of backward time τ

    for j = N:-1:1
        @inbounds begin
            # t = (j-1)*Δτ
            # TODO: pass v-column, pol-column by reference?
            v[:,j], pol[:,j] = policytimestep(model, v[:,j+1],
                                              avals, x, Δx, Δτ, j)
        end
    end

    return v, pol
end

ext: jl

max_stars_repo_name: anriseth/HJBSolver

max_issues_repo_name: anriseth/HJBSolver

max_forks_repo_name: anriseth/HJBSolver


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

doc"""
    AdaGrad(; kwargs...)

Scale learning rates by dividing with the square root of accumulated
squared gradients. See [1] for further description.

### Arguments
* `η`: default `0.1`, learning rate.
* `ϵ`: default `1e-6`, small value added for numerical stability.
* `clip`: default `0`, gradient clipping.
  If positive, will clip the gradient into the range `[-clip, clip]`.
* `scale`: default `0`, gradient rescaling.
  If != 0, multiply the gradient with `scale` before updating.
  Often choose to be `1.0 / batch_size`.
  If leave it default, high-level API like `fit!` will set it to
  `1.0 / batch_size`, since `fit!` knows the `batch_size`.
* `λ`: default `0.00001`, weight decay is equivalent
  to adding a global l2 regularizer for all the parameters.

### Notes
Using step size `η` AdaGrad calculates the learning rate for feature `i` at
time step t as:

```math
η_{t,i} = \frac{lr}{\sqrt{\sum^t_{t^\prime} g^2_{t^\prime,i} + ϵ}} g_{t,i}
```

as such the learning rate is monotonically decreasing.
Epsilon is not included in the typical formula, see [2].

### References
1. Duchi, J., Hazan, E., & Singer, Y. (2011):
   Adaptive subgradient methods for online learning and
   stochastic optimization. JMLR, 12:2121-2159.
2. Chris Dyer: Notes on AdaGrad.
   [http://www.ark.cs.cmu.edu/cdyer/adagrad.pdf]
   (http://www.ark.cs.cmu.edu/cdyer/adagrad.pdf)
"""
AdaGrad

@defstruct AdaGrad <: AbstractOptimizer (
  (η      :: Real = 0.1,  η > 0),
  (ϵ      :: Real = 1e-6, ϵ > 0),
  (clip   :: Real = 0,    clip >= 0),
   scale  :: Real = 0,
  (λ      :: Real = 1e-5, λ >= 0),
  η_sched :: Any  = initlrsched(η)
)

create_state(::AdaGrad, ::Int, W::NDArray) = zeros(size(W), context(W))

function update!(ada::AdaGrad, ::Int, W::NDArray, ∇::NDArray, x::NDArray)
  η = get(ada.η_sched)
  ϵ = ada.ϵ

  normgrad!(ada, W, ∇)

  @inplace x .+= ∇.^2  # update state
  @inplace W .+= -η .* ∇ ./ sqrt(x .+ ϵ)  # FIXME: sqrt dot-call
end

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet


text:
using CSV

@info("Loading mnist data from csv file")
train_df = CSV.file("Data/train.csv") |> DataFrame!
train_label = train_df["label"]


text:
@rule AR(:θ, Marginalisation) (q_y_x::MultivariateNormalDistributionsFamily, q_γ::GammaShapeRate, meta::ARMeta) = begin
    order = getorder(meta)
    F     = getvform(meta)

    myx, Vyx = mean_cov(q_y_x)
    my, Vy   = ar_slice(F, myx, 1:order), ar_slice(F, Vyx, 1:order, 1:order)
    mx, Vx   = ar_slice(F, myx, (order + 1):2order), ar_slice(F, Vyx, (order + 1):2order, (order + 1):2order)
    Vyx      = ar_slice(F, Vyx, (order + 1):2order, 1:order)

    mγ = mean(q_γ)

    # Equivalent to W = mγ * (Vx + mx * mx')
    W = rank1update(Vx, mx)
    W = mul_inplace!(mγ, W)

    c = ar_unit(getvform(meta), order)

    # Equivalent to ξ = (Vyx + mx * my') * c * mγ
    ξ = rank1update(Vyx, mx, my) * c
    ξ = mul_inplace!(mγ, ξ)

    return convert(promote_variate_type(F, NormalWeightedMeanPrecision), ξ, W)
end


text:
mutable struct JumpProblem{iip,P,A,C,J<:Union{Nothing,AbstractJumpAggregator},J2,J3,J4} <: DiffEqBase.AbstractJumpProblem{P,J}
  prob::P
  aggregator::A
  discrete_jump_aggregation::J
  jump_callback::C
  variable_jumps::J2
  regular_jump::J3
  massaction_jump::J4
end

DiffEqBase.isinplace(::JumpProblem{iip}) where {iip} = iip
JumpProblem(prob::JumpProblem) = prob

JumpProblem(prob,jumps::ConstantRateJump;kwargs...) = JumpProblem(prob,JumpSet(jumps);kwargs...)
JumpProblem(prob,jumps::VariableRateJump;kwargs...) = JumpProblem(prob,JumpSet(jumps);kwargs...)
JumpProblem(prob,jumps::RegularJump;kwargs...) = JumpProblem(prob,JumpSet(jumps);kwargs...)
JumpProblem(prob,jumps::MassActionJump;kwargs...) = JumpProblem(prob,JumpSet(jumps);kwargs...)
JumpProblem(prob,jumps::AbstractJump...;kwargs...) = JumpProblem(prob,JumpSet(jumps...);kwargs...)

JumpProblem(prob,aggregator::AbstractAggregatorAlgorithm,jumps::ConstantRateJump;kwargs...) = JumpProblem(prob,aggregator,JumpSet(jumps);kwargs...)
JumpProblem(prob,aggregator::AbstractAggregatorAlgorithm,jumps::VariableRateJump;kwargs...) = JumpProblem(prob,aggregator,JumpSet(jumps);kwargs...)
JumpProblem(prob,aggregator::AbstractAggregatorAlgorithm,jumps::RegularJump;kwargs...) = JumpProblem(prob,aggregator,JumpSet(jumps);kwargs...)
JumpProblem(prob,aggregator::AbstractAggregatorAlgorithm,jumps::MassActionJump;kwargs...) = JumpProblem(prob,aggregator,JumpSet(jumps);kwargs...)
JumpProblem(prob,aggregator::AbstractAggregatorAlgorithm,jumps::AbstractJump...;kwargs...) = JumpProblem(prob,aggregator,JumpSet(jumps...);kwargs...)
JumpProblem(prob,jumps::JumpSet;kwargs...) = JumpProblem(prob,NullAggregator(),jumps;kwargs...)

function JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, jumps::JumpSet;
                     save_positions = typeof(prob) <: DiffEqBase.AbstractDiscreteProblem ? (false,true) : (true,true),
                     rng = Xorshifts.Xoroshiro128Star(rand(UInt64)), kwargs...)

  ## Constant Rate Handling
  t,end_time,u = prob.tspan[1],prob.tspan[2],prob.u0
  if (typeof(jumps.constant_jumps) <: Tuple{}) && (jumps.massaction_jump === nothing)
    disc = nothing
    constant_jump_callback = CallbackSet()
  else
    disc = aggregate(aggregator,u,prob.p,t,end_time,jumps.constant_jumps,jumps.massaction_jump,save_positions,rng;kwargs...)
    constant_jump_callback = DiscreteCallback(disc)
  end

  iip = if prob isa DiscreteProblem && prob.f === DiffEqBase.DISCRETE_INPLACE_DEFAULT && jumps.regular_jump !== nothing
    # Just a default discrete problem f, so don't use it for iip
    DiffEqBase.isinplace(jumps.regular_jump)
  else
    DiffEqBase.isinplace(prob)
  end

  ## Variable Rate Handling
  if typeof(jumps.variable_jumps) <: Tuple{}
    new_prob = prob
    variable_jump_callback = CallbackSet()
  else
    new_prob = extend_problem(prob,jumps)
    variable_jump_callback = build_variable_callback(CallbackSet(),0,jumps.variable_jumps...)
  end
  callbacks = CallbackSet(constant_jump_callback,variable_jump_callback)
  JumpProblem{iip,typeof(new_prob),typeof(aggregator),typeof(callbacks),
              typeof(disc),typeof(jumps.variable_jumps),
              typeof(jumps.regular_jump),typeof(jumps.massaction_jump)}(
                        new_prob,aggregator,disc,
                        callbacks,
                        jumps.variable_jumps,
                        jumps.regular_jump, jumps.massaction_jump)
end

function extend_problem(prob::DiffEqBase.AbstractODEProblem,jumps)
  function jump_f(du::ExtendedJumpArray,u::ExtendedJumpArray,p,t)
    prob.f(du.u,u.u,p,t)
    update_jumps!(du,u,p,t,length(u.u),jumps.variable_jumps...)
  end
  u0 = ExtendedJumpArray(prob.u0,[-randexp() for i in 1:length(jumps.variable_jumps)])
  remake(prob,f=ODEFunction{true}(jump_f),u0=u0)
end

function extend_problem(prob::DiffEqBase.AbstractSDEProblem,jumps)
  function jump_f(du,u,p,t)
    prob.f(du.u,u.u,p,t)
    update_jumps!(du,u,p,t,length(u.u),jumps.variable_jumps...)
  end

  if prob.noise_rate_prototype === nothing
    jump_g = function (du,u,p,t)
      prob.g(du.u,u.u,p,t)
    end
  else
    jump_g = function (du,u,p,t)
      prob.g(du,u.u,p,t)
    end
  end

  u0 = ExtendedJumpArray(prob.u0,[-randexp() for i in 1:length(jumps.variable_jumps)])
  remake(prob,f=SDEFunction{true}(jump_f,jump_g),g=jump_g,u0=u0)
end

function extend_problem(prob::DiffEqBase.AbstractDDEProblem,jumps)
  jump_f = function (du,u,h,p,t)
    prob.f(du.u,u.u,h,p,t)
    update_jumps!(du,u,p,t,length(u.u),jumps.variable_jumps...)
  end
  u0 = ExtendedJumpArray(prob.u0,[-randexp() for i in 1:length(jumps.variable_jumps)])
  ramake(prob,f=DDEFunction{true}(jump_f),u0=u0)
end

# Not sure if the DAE one is correct: Should be a residual of sorts
function extend_problem(prob::DiffEqBase.AbstractDAEProblem,jumps)
  jump_f = function (out,du,u,p,t)
    prob.f(out.u,du.u,u.u,t)
    update_jumps!(du,u,t,length(u.u),jumps.variable_jumps...)
  end
  u0 = ExtendedJumpArray(prob.u0,[-randexp() for i in 1:length(jumps.variable_jumps)])
  remake(prob,f=DAEFunction{true}(jump_f),u0=u0)
end

function build_variable_callback(cb,idx,jump,jumps...)
  idx += 1
  condition = function (u,t,integrator)
    u.jump_u[idx]
  end
  affect! = function (integrator)
    jump.affect!(integrator)
    integrator.u.jump_u[idx] = -randexp()
  end
  new_cb = ContinuousCallback(condition,affect!;
                      idxs = jump.idxs,
                      rootfind = jump.rootfind,
                      interp_points = jump.interp_points,
                      save_positions = jump.save_positions,
                      abstol = jump.abstol,
                      reltol = jump.reltol)
  build_variable_callback(CallbackSet(cb,new_cb),idx,jumps...)
end

function build_variable_callback(cb,idx,jump)
  idx += 1
  condition = function (u,t,integrator)
    u.jump_u[idx]
  end
  affect! = function (integrator)
    jump.affect!(integrator)
    integrator.u.jump_u[idx] = -randexp()
  end
  new_cb = ContinuousCallback(condition,affect!;
                      idxs = jump.idxs,
                      rootfind = jump.rootfind,
                      interp_points = jump.interp_points,
                      save_positions = jump.save_positions,
                      abstol = jump.abstol,
                      reltol = jump.reltol)
  CallbackSet(cb,new_cb)
end

aggregator(jp::JumpProblem{P,A,C,J,J2}) where {P,A,C,J,J2} = A

@inline function extend_tstops!(tstops,jp::JumpProblem{P,A,C,J,J2}) where {P,A,C,J,J2}
  !(typeof(jp.jump_callback.discrete_callbacks) <: Tuple{}) && push!(tstops,jp.jump_callback.discrete_callbacks[1].condition.next_jump_time)
end

@inline function update_jumps!(du,u,p,t,idx,jump)
  idx += 1
  du[idx] = jump.rate(u,p,t)
end

@inline function update_jumps!(du,u,p,t,idx,jump,jumps...)
  idx += 1
  du[idx] = jump.rate(u,p,t)
  update_jumps!(du,u,p,t,idx,jumps...)
end

### Displays

Base.summary(prob::JumpProblem) = string(DiffEqBase.parameterless_type(prob)," with problem ",DiffEqBase.parameterless_type(prob.prob)," and aggregator ",typeof(prob.aggregator))
function Base.show(io::IO, A::JumpProblem)
  println(io,summary(A))
  println(io,"Number of constant rate jumps: ",A.discrete_jump_aggregation === nothing ? 0 : length(A.discrete_jump_aggregation.rates))
  println(io,"Number of variable rate jumps: ",length(A.variable_jumps))
  if A.regular_jump !== nothing
    println(io,"Have a regular jump")
  end
  if (A.massaction_jump !== nothing) && (get_num_majumps(A.massaction_jump) > 0)
    println(io,"Have a mass action jump")
  end
end

TreeViews.hastreeview(x::JumpProblem) = true


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "MiniFB"
version = v"0.1.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/emoon/minifb.git", "dda1275bd741752d1b7cdc1ee4d9941887a1891a") #master as of 10Sep2020
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd minifb/
sed -i 's/add_library(minifb STATIC/add_library(minifb SHARED/' CMakeLists.txt
mkdir build
cd build/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DMINIFB_BUILD_EXAMPLES=OFF ..
make -j${nproc}
mv libminifb* ${libdir}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64),
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libminifb", :libminifb)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Xorg_libX11_jll", uuid="4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"))
    BuildDependency(PackageSpec(name="Xorg_xorgproto_jll", uuid="c4d99508-4286-5418-9131-c86396af500b"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"5.2.0")

max_stars_repo_name: scottmishra/Yggdrasil

max_issues_repo_name: scottmishra/Yggdrasil

max_forks_repo_name: scottmishra/Yggdrasil


text:
# This file was generated, do not modify it. # hide
y_hat = predict(m_linReg, rows=test);


text:
println("##### reading #####")

using Test
using TreeTools

@testset "Basic" begin
	@test typeof(read_tree("$(dirname(pathof(TreeTools)))/../test/reading/trees.nwk")) <: Vector{<:Tree}
	@test typeof(read_tree("$(dirname(pathof(TreeTools)))/../test/reading/tree.nwk")) <: Tree
end


text:
"""
Fisher's classic iris dataset.
Measurements from 3 different species of iris: setosa, versicolor and
virginica. There are 50 examples of each species.
There are 4 measurements for each example: sepal length, sepal width,
petal length and petal width. The measurements are in centimeters.
The module retrieves the data from the [UCI Machine Learning Repository](https://archive.ics.uci.edu/ml/datasets/iris).
"""
module Iris

using DelimitedFiles
using ..Data: deps, download_and_verify, deprecation_message

# Uncomment if the iris.data file is cached to cache.julialang.org.
const cache_prefix = "https://cache.julialang.org/"

function load()
    isfile(deps("iris.data")) && return

    @info "Downloading iris dataset."
    download_and_verify("$(cache_prefix)https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data",
                        deps("iris.data"),
                        "6f608b71a7317216319b4d27b4d9bc84e6abd734eda7872b71a458569e2656c0")
end

"""
    labels()
Get the labels of the iris dataset, a 150 element array of strings listing the
species of each example.
```julia
julia> labels = Flux.Data.Iris.labels();
julia> summary(labels)
"150-element Array{String,1}"
julia> labels[1]
"Iris-setosa"
```
"""
function labels()
    deprecation_message()
    load()
    iris = readdlm(deps("iris.data"), ',')
    Vector{String}(iris[1:end, end])
end

"""
    features()
Get the features of the iris dataset. This is a 4x150 matrix of Float64
elements. It has a row for each feature (sepal length, sepal width,
petal length, petal width) and a column for each example.
```julia
julia> features = Flux.Data.Iris.features();
julia> summary(features)
"4×150 Array{Float64,2}"
julia> features[:, 1]
4-element Array{Float64,1}:
 5.1
 3.5
 1.4
 0.2
```
"""
function features()
    deprecation_message()
    load()
    iris = readdlm(deps("iris.data"), ',')
    Matrix{Float64}(iris[1:end, 1:4]')
end

end


text:
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    ::DeviceModel{T, FixedOutput},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(T, sys)
    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)
    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerTimeSeriesParameter,
        devices,
        model,
        S,
    )
    return
end

function construct_device!(
    ::OptimizationContainer,
    ::PSY.System,
    ::ModelConstructStage,
    ::DeviceModel{T, FixedOutput},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractPowerModel}
    # FixedOutput doesn't add any constraints to the model. This function covers
    # AbstractPowerModel and AbtractActivePowerModel
    return
end

"""
This function creates the arguments for the model for a full thermal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, D},
    ::Type{S},
) where {
    T <: PSY.ThermalGen,
    D <: AbstractStandardUnitCommitment,
    S <: PM.AbstractPowerModel,
}
    devices = get_available_components(T, sys)

    add_variables!(container, ActivePowerVariable, devices, D())
    add_variables!(container, ReactivePowerVariable, devices, D())
    add_variables!(container, OnVariable, devices, D())
    add_variables!(container, StartVariable, devices, D())
    add_variables!(container, StopVariable, devices, D())

    add_variables!(container, TimeDurationOn, devices, D())
    add_variables!(container, TimeDurationOff, devices, D())

    initial_conditions!(container, devices, D())

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ActivePowerRangeExpressionLB,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerRangeExpressionUB,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_feedforward_arguments!(container, model, devices)
    return
end

"""
This function creates the constraints for the model for a full thermal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, D},
    ::Type{S},
) where {
    T <: PSY.ThermalGen,
    D <: AbstractStandardUnitCommitment,
    S <: PM.AbstractPowerModel,
}
    devices = get_available_components(T, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ReactivePowerVariableLimitsConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(container, CommitmentConstraint, devices, model, S)
    add_constraints!(container, RampConstraint, devices, model, S)
    add_constraints!(container, DurationConstraint, devices, model, S)

    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

"""
This function creates the arguments model for a full thermal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, D},
    ::Type{S},
) where {
    T <: PSY.ThermalGen,
    D <: AbstractStandardUnitCommitment,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(T, sys)

    add_variables!(container, ActivePowerVariable, devices, D())
    add_variables!(container, OnVariable, devices, D())
    add_variables!(container, StartVariable, devices, D())
    add_variables!(container, StopVariable, devices, D())

    add_variables!(container, TimeDurationOn, devices, D())
    add_variables!(container, TimeDurationOff, devices, D())

    initial_conditions!(container, devices, D())

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ActivePowerRangeExpressionLB,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerRangeExpressionUB,
        ActivePowerVariable,
        devices,
        model,
        S,
    )

    add_feedforward_arguments!(container, model, devices)
    return
end

"""
This function creates the constraints for the model for a full thermal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, D},
    ::Type{S},
) where {
    T <: PSY.ThermalGen,
    D <: AbstractStandardUnitCommitment,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(T, sys)
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )

    add_constraints!(container, CommitmentConstraint, devices, model, S)
    add_constraints!(container, RampConstraint, devices, model, S)
    add_constraints!(container, DurationConstraint, devices, model, S)

    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)

    add_constraint_dual!(container, sys, model)
    return
end

"""
This function creates the model for a full thermal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, ThermalBasicUnitCommitment},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractPowerModel}
    devices = get_available_components(T, sys)

    add_variables!(container, ActivePowerVariable, devices, ThermalBasicUnitCommitment())
    add_variables!(container, ReactivePowerVariable, devices, ThermalBasicUnitCommitment())
    add_variables!(container, OnVariable, devices, ThermalBasicUnitCommitment())
    add_variables!(container, StartVariable, devices, ThermalBasicUnitCommitment())
    add_variables!(container, StopVariable, devices, ThermalBasicUnitCommitment())

    initial_conditions!(container, devices, ThermalBasicUnitCommitment())

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ActivePowerRangeExpressionLB,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerRangeExpressionUB,
        ActivePowerVariable,
        devices,
        model,
        S,
    )

    add_feedforward_arguments!(container, model, devices)
    return
end

"""
This function creates the model for a full thermal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, ThermalBasicUnitCommitment},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractPowerModel}
    devices = get_available_components(T, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )

    add_constraints!(
        container,
        ReactivePowerVariableLimitsConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(container, CommitmentConstraint, devices, model, S)

    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

"""
This function creates the arguments for the model for a full thermal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, ThermalBasicUnitCommitment},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(T, sys)

    add_variables!(container, ActivePowerVariable, devices, ThermalBasicUnitCommitment())
    add_variables!(container, OnVariable, devices, ThermalBasicUnitCommitment())
    add_variables!(container, StartVariable, devices, ThermalBasicUnitCommitment())
    add_variables!(container, StopVariable, devices, ThermalBasicUnitCommitment())

    initial_conditions!(container, devices, ThermalBasicUnitCommitment())

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ActivePowerRangeExpressionLB,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerRangeExpressionUB,
        ActivePowerVariable,
        devices,
        model,
        S,
    )

    add_feedforward_arguments!(container, model, devices)
    return
end

"""
This function creates the constraints for the model for a full thermal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, ThermalBasicUnitCommitment},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(T, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )

    add_constraints!(container, CommitmentConstraint, devices, model, S)

    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

"""
This function creates the model for a full thermal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, ThermalStandardDispatch},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractPowerModel}
    devices = get_available_components(T, sys)

    add_variables!(container, ActivePowerVariable, devices, ThermalStandardDispatch())
    add_variables!(container, ReactivePowerVariable, devices, ThermalStandardDispatch())

    initial_conditions!(container, devices, ThermalStandardDispatch())

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ActivePowerRangeExpressionLB,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerRangeExpressionUB,
        ActivePowerVariable,
        devices,
        model,
        S,
    )

    add_feedforward_arguments!(container, model, devices)
    return
end

"""
This function creates the model for a full thermal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, ThermalStandardDispatch},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractPowerModel}
    devices = get_available_components(T, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )

    add_constraints!(
        container,
        ReactivePowerVariableLimitsConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(container, RampConstraint, devices, model, S)

    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

"""
This function creates the arguments for the model for a full thermal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, ThermalStandardDispatch},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(T, sys)

    add_variables!(container, ActivePowerVariable, devices, ThermalStandardDispatch())

    initial_conditions!(container, devices, ThermalStandardDispatch())

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ActivePowerRangeExpressionLB,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerRangeExpressionUB,
        ActivePowerVariable,
        devices,
        model,
        S,
    )

    add_feedforward_arguments!(container, model, devices)
    return
end

"""
This function creates the constraints for the model for a full thermal dispatch formulation depending on combination of devices, device_formulation and system_formulation
"""
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, ThermalStandardDispatch},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(T, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )

    add_constraints!(container, RampConstraint, devices, model, S)

    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, D},
    ::Type{S},
) where {
    T <: PSY.ThermalGen,
    D <: AbstractThermalDispatchFormulation,
    S <: PM.AbstractPowerModel,
}
    devices = get_available_components(T, sys)

    add_variables!(container, ActivePowerVariable, devices, D())
    add_variables!(container, ReactivePowerVariable, devices, D())

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ActivePowerRangeExpressionLB,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerRangeExpressionUB,
        ActivePowerVariable,
        devices,
        model,
        S,
    )

    add_feedforward_arguments!(container, model, devices)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, D},
    ::Type{S},
) where {
    T <: PSY.ThermalGen,
    D <: AbstractThermalDispatchFormulation,
    S <: PM.AbstractPowerModel,
}
    devices = get_available_components(T, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )

    add_constraints!(
        container,
        ReactivePowerVariableLimitsConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )

    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, D},
    ::Type{S},
) where {
    T <: PSY.ThermalGen,
    D <: AbstractThermalDispatchFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(T, sys)

    add_variables!(container, ActivePowerVariable, devices, D())

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ActivePowerRangeExpressionLB,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerRangeExpressionUB,
        ActivePowerVariable,
        devices,
        model,
        S,
    )

    add_feedforward_arguments!(container, model, devices)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, D},
    ::Type{S},
) where {
    T <: PSY.ThermalGen,
    D <: AbstractThermalDispatchFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(T, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )

    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{PSY.ThermalMultiStart, ThermalMultiStartUnitCommitment},
    ::Type{S};
    kwargs...,
) where {S <: PM.AbstractPowerModel}
    devices = get_available_components(PSY.ThermalMultiStart, sys)

    add_variables!(
        container,
        PowerAboveMinimumVariable,
        devices,
        ThermalMultiStartUnitCommitment(),
    )
    add_variables!(
        container,
        ReactivePowerVariable,
        devices,
        ThermalMultiStartUnitCommitment(),
    )
    add_variables!(container, OnVariable, devices, ThermalMultiStartUnitCommitment())
    add_variables!(container, StopVariable, devices, ThermalMultiStartUnitCommitment())
    add_variables!(container, StartVariable, devices, ThermalMultiStartUnitCommitment())
    add_variables!(container, ColdStartVariable, devices, ThermalMultiStartUnitCommitment())
    add_variables!(container, WarmStartVariable, devices, ThermalMultiStartUnitCommitment())
    add_variables!(container, HotStartVariable, devices, ThermalMultiStartUnitCommitment())

    add_variables!(container, TimeDurationOn, devices, ThermalMultiStartUnitCommitment())
    add_variables!(container, TimeDurationOff, devices, ThermalMultiStartUnitCommitment())
    add_variables!(container, PowerOutput, devices, ThermalMultiStartUnitCommitment())

    initial_conditions!(container, devices, ThermalMultiStartUnitCommitment())

    add_to_expression!(
        container,
        ActivePowerBalance,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(container, ActivePowerBalance, OnVariable, devices, model, S)
    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ActivePowerRangeExpressionLB,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerRangeExpressionUB,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )

    add_feedforward_arguments!(container, model, devices)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{PSY.ThermalMultiStart, ThermalMultiStartUnitCommitment},
    ::Type{S};
    kwargs...,
) where {S <: PM.AbstractPowerModel}
    devices = get_available_components(PSY.ThermalMultiStart, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )

    add_constraints!(
        container,
        ReactivePowerVariableLimitsConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(container, CommitmentConstraint, devices, model, S)
    add_constraints!(container, RampConstraint, devices, model, S)
    add_constraints!(container, DurationConstraint, devices, model, S)
    add_constraints!(container, StartupTimeLimitTemperatureConstraint, devices, model, S)
    add_constraints!(container, StartTypeConstraint, devices, model, S)
    add_constraints!(container, StartupInitialConditionConstraint, devices, model, S)
    add_constraints!(container, MustRunConstraint, devices, model, S)
    add_constraints!(container, ActiveRangeICConstraint, devices, model, S)

    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{PSY.ThermalMultiStart, ThermalMultiStartUnitCommitment},
    ::Type{S},
) where {S <: PM.AbstractActivePowerModel}
    devices = get_available_components(PSY.ThermalMultiStart, sys)

    add_variables!(
        container,
        PowerAboveMinimumVariable,
        devices,
        ThermalMultiStartUnitCommitment(),
    )
    add_variables!(container, OnVariable, devices, ThermalMultiStartUnitCommitment())
    add_variables!(container, StopVariable, devices, ThermalMultiStartUnitCommitment())
    add_variables!(container, StartVariable, devices, ThermalMultiStartUnitCommitment())
    add_variables!(container, ColdStartVariable, devices, ThermalMultiStartUnitCommitment())
    add_variables!(container, WarmStartVariable, devices, ThermalMultiStartUnitCommitment())
    add_variables!(container, HotStartVariable, devices, ThermalMultiStartUnitCommitment())

    add_variables!(container, TimeDurationOn, devices, ThermalMultiStartUnitCommitment())
    add_variables!(container, TimeDurationOff, devices, ThermalMultiStartUnitCommitment())
    add_variables!(container, PowerOutput, devices, ThermalMultiStartUnitCommitment())

    add_to_expression!(
        container,
        ActivePowerBalance,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(container, ActivePowerBalance, OnVariable, devices, model, S)

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ActivePowerRangeExpressionLB,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerRangeExpressionUB,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )

    add_feedforward_arguments!(container, model, devices)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{PSY.ThermalMultiStart, ThermalMultiStartUnitCommitment},
    ::Type{S},
) where {S <: PM.AbstractActivePowerModel}
    devices = get_available_components(PSY.ThermalMultiStart, sys)

    initial_conditions!(container, devices, ThermalMultiStartUnitCommitment())

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )

    add_constraints!(container, CommitmentConstraint, devices, model, S)
    add_constraints!(container, RampConstraint, devices, model, S)
    add_constraints!(container, DurationConstraint, devices, model, S)
    add_constraints!(container, StartupTimeLimitTemperatureConstraint, devices, model, S)
    add_constraints!(container, StartTypeConstraint, devices, model, S)
    add_constraints!(container, StartupInitialConditionConstraint, devices, model, S)
    add_constraints!(container, MustRunConstraint, devices, model, S)
    add_constraints!(container, ActiveRangeICConstraint, devices, model, S)

    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, ThermalCompactUnitCommitment},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractPowerModel}
    devices = get_available_components(T, sys)

    add_variables!(
        container,
        PowerAboveMinimumVariable,
        devices,
        ThermalCompactUnitCommitment(),
    )
    add_variables!(
        container,
        ReactivePowerVariable,
        devices,
        ThermalCompactUnitCommitment(),
    )
    add_variables!(container, OnVariable, devices, ThermalCompactUnitCommitment())
    add_variables!(container, StartVariable, devices, ThermalCompactUnitCommitment())
    add_variables!(container, StopVariable, devices, ThermalCompactUnitCommitment())

    add_variables!(container, TimeDurationOn, devices, ThermalCompactUnitCommitment())
    add_variables!(container, TimeDurationOff, devices, ThermalCompactUnitCommitment())
    add_variables!(container, PowerOutput, devices, ThermalCompactUnitCommitment())

    initial_conditions!(container, devices, ThermalCompactUnitCommitment())

    add_to_expression!(
        container,
        ActivePowerBalance,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(container, ActivePowerBalance, OnVariable, devices, model, S)

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ActivePowerRangeExpressionLB,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerRangeExpressionUB,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )

    add_feedforward_arguments!(container, model, devices)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, ThermalCompactUnitCommitment},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractPowerModel}
    devices = get_available_components(T, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )

    add_constraints!(
        container,
        ReactivePowerVariableLimitsConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(container, CommitmentConstraint, devices, model, S)
    add_constraints!(container, RampConstraint, devices, model, S)
    add_constraints!(container, DurationConstraint, devices, model, S)

    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, ThermalCompactUnitCommitment},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(T, sys)

    add_variables!(
        container,
        PowerAboveMinimumVariable,
        devices,
        ThermalCompactUnitCommitment(),
    )
    add_variables!(container, OnVariable, devices, ThermalCompactUnitCommitment())
    add_variables!(container, StartVariable, devices, ThermalCompactUnitCommitment())
    add_variables!(container, StopVariable, devices, ThermalCompactUnitCommitment())

    add_variables!(container, TimeDurationOn, devices, ThermalCompactUnitCommitment())
    add_variables!(container, TimeDurationOff, devices, ThermalCompactUnitCommitment())
    add_variables!(container, PowerOutput, devices, ThermalCompactUnitCommitment())

    initial_conditions!(container, devices, ThermalCompactUnitCommitment())

    add_to_expression!(
        container,
        ActivePowerBalance,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(container, ActivePowerBalance, OnVariable, devices, model, S)

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ActivePowerRangeExpressionLB,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerRangeExpressionUB,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    add_feedforward_arguments!(container, model, devices)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, ThermalCompactUnitCommitment},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(T, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )

    add_constraints!(container, CommitmentConstraint, devices, model, S)
    add_constraints!(container, RampConstraint, devices, model, S)
    add_constraints!(container, DurationConstraint, devices, model, S)

    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, ThermalBasicCompactUnitCommitment},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractPowerModel}
    devices = get_available_components(T, sys)

    add_variables!(
        container,
        PowerAboveMinimumVariable,
        devices,
        ThermalBasicCompactUnitCommitment(),
    )
    add_variables!(
        container,
        ReactivePowerVariable,
        devices,
        ThermalBasicCompactUnitCommitment(),
    )
    add_variables!(container, OnVariable, devices, ThermalBasicCompactUnitCommitment())
    add_variables!(container, StartVariable, devices, ThermalBasicCompactUnitCommitment())
    add_variables!(container, StopVariable, devices, ThermalBasicCompactUnitCommitment())

    add_variables!(container, PowerOutput, devices, ThermalBasicCompactUnitCommitment())

    initial_conditions!(container, devices, ThermalBasicCompactUnitCommitment())

    add_to_expression!(
        container,
        ActivePowerBalance,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(container, ActivePowerBalance, OnVariable, devices, model, S)

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ActivePowerRangeExpressionLB,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerRangeExpressionUB,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )

    add_feedforward_arguments!(container, model, devices)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, ThermalBasicCompactUnitCommitment},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractPowerModel}
    devices = get_available_components(T, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )

    add_constraints!(
        container,
        ReactivePowerVariableLimitsConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(container, CommitmentConstraint, devices, model, S)

    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, ThermalBasicCompactUnitCommitment},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(T, sys)

    add_variables!(
        container,
        PowerAboveMinimumVariable,
        devices,
        ThermalBasicCompactUnitCommitment(),
    )
    add_variables!(container, OnVariable, devices, ThermalBasicCompactUnitCommitment())
    add_variables!(container, StartVariable, devices, ThermalBasicCompactUnitCommitment())
    add_variables!(container, StopVariable, devices, ThermalBasicCompactUnitCommitment())

    add_variables!(container, PowerOutput, devices, ThermalBasicCompactUnitCommitment())

    initial_conditions!(container, devices, ThermalBasicCompactUnitCommitment())

    add_to_expression!(
        container,
        ActivePowerBalance,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(container, ActivePowerBalance, OnVariable, devices, model, S)

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ActivePowerRangeExpressionLB,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerRangeExpressionUB,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    add_feedforward_arguments!(container, model, devices)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, ThermalBasicCompactUnitCommitment},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(T, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )

    add_constraints!(container, CommitmentConstraint, devices, model, S)

    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, ThermalCompactDispatch},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractPowerModel}
    devices = get_available_components(T, sys)

    add_variables!(container, PowerAboveMinimumVariable, devices, ThermalCompactDispatch())
    add_variables!(container, ReactivePowerVariable, devices, ThermalCompactDispatch())

    add_variables!(container, PowerOutput, devices, ThermalCompactDispatch())

    add_parameters!(container, OnStatusParameter(), devices, model)

    add_feedforward_arguments!(container, model, devices)

    initial_conditions!(container, devices, ThermalCompactDispatch())

    add_to_expression!(
        container,
        ActivePowerBalance,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(container, ActivePowerBalance, OnStatusParameter, devices, model, S)
    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(container, ActivePowerBalance, OnStatusParameter, devices, model, S)
    add_to_expression!(
        container,
        ActivePowerRangeExpressionLB,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerRangeExpressionUB,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, ThermalCompactDispatch},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractPowerModel}
    devices = get_available_components(T, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )

    add_constraints!(
        container,
        ReactivePowerVariableLimitsConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(container, RampConstraint, devices, model, S)

    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{T, ThermalCompactDispatch},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(T, sys)

    add_variables!(container, PowerAboveMinimumVariable, devices, ThermalCompactDispatch())

    add_variables!(container, PowerOutput, devices, ThermalCompactDispatch())

    add_parameters!(container, OnStatusParameter(), devices, model)

    add_feedforward_arguments!(container, model, devices)

    add_to_expression!(
        container,
        ActivePowerBalance,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )

    add_to_expression!(container, ActivePowerBalance, OnStatusParameter, devices, model, S)

    initial_conditions!(container, devices, ThermalCompactDispatch())

    add_expressions!(container, ProductionCostExpression, devices, model)

    add_to_expression!(
        container,
        ActivePowerRangeExpressionLB,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ActivePowerRangeExpressionUB,
        PowerAboveMinimumVariable,
        devices,
        model,
        S,
    )
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{T, ThermalCompactDispatch},
    ::Type{S},
) where {T <: PSY.ThermalGen, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(T, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionLB,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerRangeExpressionUB,
        devices,
        model,
        S,
    )

    add_constraints!(container, RampConstraint, devices, model, S)

    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)
    add_constraint_dual!(container, sys, model)
    return
end


text:
using Bokeh
using Compat
if VERSION < v"0.4-"
    using Dates
else
    using Base.Dates
end

autoopen(true)

plotfile("dates1.html")

# this simple example uses unix time stamps for datetimes
epoch = DateTime(1970, 1, 1)
@compat n = Int(now() - epoch) #  ms
# 1 month range
day = 1000 * 3600 * 24
x = linspace(n, n + day * 30, 1000)
y = 0.5sin(x * pi / day) + cos(x * pi / (day * 7))
# autoopen is set to true so `plot` will try to open a browser to show the plot
plot(x, y, x_axis_type=:datetime)


text:
#Function for getting atom specific parameters for Rydberg code. Based on Matlab code written in the summer of 2015 migrated to Julia on 23/07/2017
#The intention is not that this be a complete database but rather that more parameters be added as they are needed

function GetAtomParams(atom,nn,ll,jj)

  #Parameters for Rubidium-8
  if atom=="87Rb"
    ZZ, spin, eneigval, params, delta_nlj = Rb87numbers(nn,ll,jj)
  end

  return ZZ, spin, eneigval, params, delta_nlj
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
using Graphs
using GraphViz

g = simple_graph(3)  #Directed Graph (3 vertices, 0 edges)
add_edge!(g, 1, 2)   #edge [1]: 1 -- 2
add_edge!(g, 3, 2)   #edge [2]: 3 -- 2
add_edge!(g, 3, 1)   #edge [3]: 3 -- 1

GraphViz.Graph(Graphs.to_dot(g))


text:
using IndividualDisplacements, MeshArrays, MITgcmTools, OrdinaryDiffEq, Statistics

"""
    example1()

Pre-computed global ocean case. Here we just re-read data from a file produced
earlier, rather than computing trajectories as in the other examples.

```
df=example1()

p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/recipes_pyplot.jl"))
PyPlot.figure(); PlotMapProj(df,300); gcf()
```
"""
function example1()
   p=dirname(pathof(IndividualDisplacements))
   dirIn=joinpath(p,"../examples/run_offflt/")
   prec=Float32
   df=read_flt(dirIn,prec)
end

"""
    example2()

Reproducing `MITgcm/verification/flt_example/` case. This is based on an
extended and modified configuration of the standard MITgcm test case.

```
(𝐼,df,ref)=example2();

p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/recipes_plots.jl"))
PlotBasic(df,300,100000.0)

using Plots
Plots.plot(𝐼.🔴.x,𝐼.🔴.y,linewidth=5,lc=:black, title="One Trajectory Example",
xaxis="x",yaxis="y",label="Julia Solution") # legend=false
pl=Plots.plot!(ref[1,:],ref[2,:],lw=3,ls=:dash,lc=:red,label="MITgcm Solution")
```
"""
function example2()
   𝑃,Γ=example2_setup()
   (xy,df,ref,nSteps)=example2_xy(𝑃)

   𝑃.𝑇[:] = [0.0,nSteps*3600.0]
   solv(prob) = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)

   tr = DataFrame(ID=Int[], x=Float64[], y=Float64[], t=Float64[])

   𝐼 = Individuals{Float64}(📌=xy[:,:], 🔴=tr, 🚄 = dxdt!, ∫ = solv, 🔧 = postprocess_xy, 𝑃=𝑃)
   𝑇=(0.0,𝐼.𝑃.𝑇[2])
   ∫!(𝐼,𝑇)

   return 𝐼, df,ref
end

"""
example2_xy()

Read MITgcm/pkg/flt output
"""
function example2_xy(𝑃)
dirIn=IndividualDisplacements.flt_example
prec=Float32
df=read_flt(dirIn,prec)
#
tmp=df[df.ID .== 200, :]
nSteps=Int32(tmp[end,:time]/3600)-2
ref=transpose([tmp[1:nSteps,:lon] tmp[1:nSteps,:lat]])
dx=5000.0
maxLon=80*dx
maxLat=42*dx
for i=1:nSteps-1
    ref[1,i+1]-ref[1,i]>maxLon/2 ? ref[1,i+1:end]-=fill(maxLon,(nSteps-i)) : nothing
    ref[1,i+1]-ref[1,i]<-maxLon/2 ? ref[1,i+1:end]+=fill(maxLon,(nSteps-i)) : nothing
    ref[2,i+1]-ref[2,i]>maxLat/2 ? ref[2,i+1:end]-=fill(maxLat,(nSteps-i)) : nothing
    ref[2,i+1]-ref[2,i]<-maxLat/2 ? ref[2,i+1:end]+=fill(maxLat,(nSteps-i)) : nothing
end
ref=ref./dx
xy=[tmp[1,:lon];tmp[1,:lat]]./dx
return xy,df,ref,nSteps
end

"""
example2_setup()

Read gridded variables from file using MeshArrays and
return result in uvetc Dictionary.
"""
function example2_setup()

   ###### 1) Get gridded variables via MeshArrays.jl

   dirIn=IndividualDisplacements.flt_example
   γ=gcmgrid(dirIn,"PeriodicChannel",1,[(80,42)], [80 42], Float32, read, write)
   nr=8

   ## Put grid variables in a dictionary:

   Γ=Dict("XC" => read_mds(γ.path*"XC",MeshArray(γ,Float32)),
   "YC" => read_mds(γ.path*"YC",MeshArray(γ,Float32)),
   "XG" => read_mds(γ.path*"XG",MeshArray(γ,Float32)),
   "YG" => read_mds(γ.path*"YG",MeshArray(γ,Float32)),
   "dx" => 5000.0)

   ## Put velocity fields in a dictionary:

   t0=0.0 #approximation / simplification
   t1=18001.0*3600.0

   u0=read_mds(γ.path*"U.0000000001",MeshArray(γ,Float32,nr))
   u1=read_mds(γ.path*"U.0000018001",MeshArray(γ,Float32,nr))
   v0=read_mds(γ.path*"V.0000000001",MeshArray(γ,Float32,nr))
   v1=read_mds(γ.path*"V.0000018001",MeshArray(γ,Float32,nr))

   kk=3 #3 to match -1406.25 in pkg/flt output
   u0=u0[:,kk]; u1=u1[:,kk];
   v0=v0[:,kk]; v1=v1[:,kk];

   u0=u0./Γ["dx"]
   u1=u1./Γ["dx"]
   v0=v0./Γ["dx"]
   v1=v1./Γ["dx"]

   ## Visualize velocity fields

   mskW=read_mds(γ.path*"hFacW",MeshArray(γ,Float32,nr))
   mskW=1.0 .+ 0.0 * mask(mskW[:,kk],NaN,0.0)
   mskS=read_mds(γ.path*"hFacS",MeshArray(γ,Float32,nr))
   mskS=1.0 .+ 0.0 * mask(mskS[:,kk],NaN,0.0)
   Γ=merge(Γ,Dict("mskW" => mskW, "mskS" => mskS))

   𝑃=FlowFields(u0[1], u1[1], v0[1], v1[1], [t0,t1])
   return 𝑃,Γ
end


text:
export skip_error

import Base: show

"""
    skip_error()

Creates a `skip_error` operator, which filters out `error` event by the source Observable by emitting only
`next` and `complete` messages.

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = faulted("error")
subscribe!(source |> skip_error(), logger())
;

# output

```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref), [`skip_error`](@ref), [`skip_complete`](@ref), [`logger`](@ref)
"""
skip_error() = SkipErrorOperator()

struct SkipErrorOperator <: InferableOperator end

function on_call!(::Type{L}, ::Type{L}, operator::SkipErrorOperator, source) where L
    return proxy(L, source, SkipErrorProxy())
end

operator_right(operator::SkipErrorOperator, ::Type{L}) where L = L

struct SkipErrorProxy <: ActorProxy end

actor_proxy!(::Type{L}, proxy::SkipErrorProxy, actor::A) where { L, A } = SkipErrorActor{L, A}(actor)

struct SkipErrorActor{L, A} <: Actor{L}
    actor :: A
end

on_next!(actor::SkipErrorActor{L}, data::L) where L = next!(actor.actor, data)
on_error!(actor::SkipErrorActor, err)               = begin end
on_complete!(actor::SkipErrorActor)                 = complete!(actor.actor)

Base.show(io::IO, ::SkipErrorOperator)         = print(io, "SkipErrorOperator()")
Base.show(io::IO, ::SkipErrorProxy)            = print(io, "SkipErrorProxy()")
Base.show(io::IO, ::SkipErrorActor{L}) where L = print(io, "SkipErrorActor($L)")


text:
# # One-dimensional cellular automata
# Agents.jl provides a module (CA1D) to create and plot 1D cellular automata.

using Agents
using Agents.CA1D
using AgentsPlots

# ## 1. Define the rule
# Here is Wolfram's rule 22

rules = Dict(
    "111" => "0",
    "110" => "0",
    "101" => "0",
    "100" => "1",
    "011" => "0",
    "010" => "1",
    "001" => "1",
    "000" => "0",
);

# ## 2. Build the model
# All the cells are by default initially "off"
model = CA1D.build_model(rules = rules, ncols = 21)

# This turns on the middle cell
model.agents[11].status = "1";

# ## 3. Run the model and collect data.
runs = 20
data, _ = CA1D.ca_run(model, runs);

# ## 4. Visualize the data
p = plot_CA1D(data, nodesize = 2)


text:
include("day6_solution.jl")

entry = "abcx
abcy
abcz"
@test get_unique_yes_questions_in_group(entry) == Set(['a','b','c','x','y','z'])

shared_entry = "ymw
w
wm
vsw
wm"
@test get_shared_yes_questions_in_group(shared_entry) == Set(['w'])


text:
# This file contains code that was a part of Julia. License is MIT: https://julialang.org/license
# It fixes problems with dealing with AbstractString types that don't use the Char type

import Base: parseint_iterate, parseint_preamble, tryparse_internal
import Base.Checked: add_with_overflow, mul_with_overflow

function parseint_iterate(s::T, startpos::Int, endpos::Int) where {T<:Str}
    C = eltype(T)
    (0 < startpos <= endpos) || (return C(0), 0, 0)
    j = startpos
    c, startpos = iterate(s,startpos)::Tuple{C, Int}
    c, startpos, j
end

function parseint_preamble(signed::Bool, base::Int, s::T, startpos::Int, endpos::Int
                           ) where {T<:Str}
    c, i, j = parseint_iterate(s, startpos, endpos)

    while isspace(c)
        c, i, j = parseint_iterate(s,i,endpos)
    end
    (j == 0) && (return 0, 0, 0)

    sgn = 1
    if signed
        if c == '-' || c == '+'
            (c == '-') && (sgn = -1)
            c, i, j = parseint_iterate(s,i,endpos)
        end
    end

    while isspace(c)
        c, i, j = parseint_iterate(s,i,endpos)
    end
    (j == 0) && (return 0, 0, 0)

    if base == 0
        if c == '0' && i <= ncodeunits(s)
            c, i = iterate(s,i)::Tuple{eltype(s), Int}
            base = c=='b' ? 2 : c=='o' ? 8 : c=='x' ? 16 : 10
            if base != 10
                c, i, j = parseint_iterate(s,i,endpos)
            end
        else
            base = 10
        end
    end
    return sgn, base, j
end

endparse(str::String, raise::Bool, err) = raise ? throw(err(str)) : nothing
endparse(str::String, raise::Bool) = endparse(str,raise, ArgumentError)

@inline _rs(str, startpos, endpos) = repr(SubString(str, startpos, endpos))

function tryparse_internal(::Type{BigInt}, s::S, startpos::Int, endpos::Int,
                           base_::Integer, raise::Bool) where {S<:Str}
    str = convert(String, SubString(s, startpos, endpos))
    Base.tryparse_internal(BigInt, str, 1, lastindex(str), base_, raise)
end

function tryparse_internal(::Type{T}, s::S, startpos::Int, endpos::Int,
                           base_::Integer, raise::Bool) where {T<:Integer, S<:Str}
    C = eltype(S)
    sgn, base, i = parseint_preamble(T<:Signed, Int(base_), s, startpos, endpos)
    sgn == 0 && base == 0 && i == 0 &&
        return endparse("input string is empty or only contains whitespace", raise)
    (2 <= base <= 62) ||
        return endparse("invalid base: base must be 2 ≤ base ≤ 62, got $base", raise)
    i == 0 &&
        return endparse("premature end of integer: " * _rs(s, startpos, endpos), raise)
    c, i = parseint_iterate(s,i,endpos)
    i == 0 &&
        return endparse("premature end of integer: " * _rs(s, startpos, endpos), raise)

    base = convert(T,base)
    m::T = div(typemax(T)-base+1,base)
    n::T = 0
    a::Int = base <= 36 ? 10 : 36
    while n <= m
        d::T = '0' <= c <= '9' ? c-'0'    :
               'A' <= c <= 'Z' ? c-'A'+10 :
               'a' <= c <= 'z' ? c-'a'+a  : base
        d >= base &&
            return endparse("invalid base $base digit $(repr(c)) in " * _rs(s,startpos,endpos),
                            raise)
        n *= base
        n += d
        if i > endpos
            n *= sgn
            return n
        end
        c, i = iterate(s,i)::Tuple{C, Int}
        isspace(c) && break
    end
    (T <: Signed) && (n *= sgn)
    while !isspace(c)
        d::T = '0' <= c <= '9' ? c-'0'    :
        'A' <= c <= 'Z' ? c-'A'+10 :
            'a' <= c <= 'z' ? c-'a'+a  : base
        d >= base &&
            return endparse("invalid base $base digit $(repr(c)) in " * _rs(s, startpos, endpos),
                            raise)
        (T <: Signed) && (d *= sgn)

        n, ov_mul = mul_with_overflow(n, base)
        n, ov_add = add_with_overflow(n, d)
        (ov_mul | ov_add) &&
            return endparse("overflow parsing " * _rs(s,startpos,endpos), raise, OverflowError)
        (i > endpos) && return n
        c, i = iterate(s,i)::Tuple{C, Int}
    end
    while i <= endpos
        c, i = iterate(s,i)::Tuple{C, Int}
        isspace(c) ||
            return endparse("extra characters after whitespace in " * _rs(s,startpos,endpos),
                            raise)
    end
    return n
end


text:
"""
    ECA <: AbstractAutomaton{Int,1}

An elementary cellular automaton.

ECA's consist of an 8-bit Wolfram code, and a boundary condition. If the boundary condition
`isnull` then a circular boundary is used, otherwise the tuple represents the (left,right)
boundary condition.

# Examples
```julia
julia> ECA(30)
Annette.ECA(0x1e, Nullable{Tuple{Int64,Int64}}())

julia> ECA(30, 1, 0)
Annette.ECA(0x1e, Nullable{Tuple{Int64,Int64}}((1, 0)))

julia> ECA(30, (1,0))
Annette.ECA(0x1e, Nullable{Tuple{Int64,Int64}}((1, 0)))
```
"""
type ECA <: AbstractAutomaton{Int,1}
    "The 8-bit Wolfram ECA code"
    code :: UInt8
    "The boundary conditions for the ECA"
    boundary :: Nullable{Tuple{Int,Int}}

    ECA(code) = new(code, Nullable{Tuple{Int,Int}}())
    ECA(code, left, right) = new(code, Nullable((1 & left, 1 & right)))
end

ECA(code, boundary::Tuple{Int,Int}) = ECA(code, boundary...)

"""
    iscircular(ca::ECA)

Does the ECA have circular or open boundary conditions?

# Examples
```julia
julia> iscircular(ECA(30))
true

julia> iscircular(ECA(30, (0,0)))
false
```
"""
iscircular(ca::ECA) = isnull(ca.boundary)

function fire!(ca::ECA, state::AbstractVector{Int}, dest::AbstractVector{Int})
    if length(dest) != length(state)
        throw(BoundsError())
    end

    left, right = iscircular(ca) ? (state[end], state[1]) : get(ca.boundary)

    shift = 2 * left + state[1]
    for i in 1:length(state)-1
        shift = 7 & (2 * shift + state[i+1])
        dest[i] = 1 & (ca.code >> shift)
    end
    shift = 7 & (2 * shift + right)
    dest[end] = 1 & (ca.code >> shift)

    dest
end


text:
# push!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))
using JLD2
jlfile = jldopen("data/local.jld2", "a+")
# jlfile = jldopen("data/nonlocal.jld2", "a+")
# jlfile = jldopen("data/nonlocal_more_velocity.jld2", "a+")

κ¹¹ = jlfile["averagelocaldiffusivity"]["κ11"]
κ¹² = jlfile["averagelocaldiffusivity"]["κ12"]
κ²¹ = jlfile["averagelocaldiffusivity"]["κ21"]
κ²² = jlfile["averagelocaldiffusivity"]["κ22"]
z = jlfile["grid"]["z"][:]
x = jlfile["grid"]["x"][:]
##
# Contruct local diffusivity estimate from flow field and γ
γ = jlfile["parameters"]["γ"]
ω = jlfile["parameters"]["ω"]
u¹ = jlfile["velocities"]["u¹"]
u² = jlfile["velocities"]["u²"]
v¹ = jlfile["velocities"]["v¹"]
v² = jlfile["velocities"]["v²"]

# integrating out the first dimension amounts to assuming that the dominant mode
# is the zero'th.
N = length(x)
scale = (2 * (ω^2 + γ^2))
analytic_κ¹¹ = sum((γ * u¹ .* u¹ + γ * u² .* u²) ./ scale, dims = 1)[:] ./ N
analytic_κ¹² = sum((u¹ .* (γ * v¹ + ω * v²) + u² .* (γ * v² - ω * v¹)) ./ scale, dims = 1)[:] ./ N
analytic_κ²¹ = sum((v¹ .* (γ * u¹ + ω * u²) + v² .* (γ * u² - ω * u¹)) ./ scale, dims = 1)[:] ./ N
analytic_κ²² = sum((γ * v¹ .* v¹ + γ * v² .* v²) ./ scale, dims = 1)[:] ./ N

##
using GLMakie

options = (; ylabel = "z", ylabelsize = 32,
    xlabelsize = 32, xgridstyle = :dash, ygridstyle = :dash, xtickalign = 1,
    xticksize = 30, ytickalign = 1, yticksize = 30,
    xticklabelsize = 30, yticklabelsize = 30)

fig = Figure(resolution = (1800, 1300), title = "Local Operators")
titlestring = "Kˣˣ"
ax1 = Axis(fig[1, 1]; options..., title = titlestring, titlesize = 30)
titlestring = "Kˣᶻ"
ax2 = Axis(fig[1, 2]; options..., title = titlestring, titlesize = 30)
titlestring = "Kᶻˣ"
ax3 = Axis(fig[2, 1]; options..., title = titlestring, titlesize = 30)
titlestring = "Kᶻᶻ"
ax4 = Axis(fig[2, 2]; options..., title = titlestring, titlesize = 30)

plot_string_1 = "analytic"
plot_string_2 = "numerical"

ln1 = lines!(ax1, analytic_κ¹¹, z, color = :red)
sc1 = scatter!(ax1, κ¹¹, z, color = :blue)
axislegend(ax1, [ln1, sc1], [plot_string_1, plot_string_2], position = :rc, labelsize = 30)

ln2 = lines!(ax2, analytic_κ¹², z, color = :red)
sc2 = scatter!(ax2, κ¹², z, color = :blue)
# axislegend(ax2, [ln2, sc2], [plot_string_1, plot_string_2], position = :rt, labelsize = 30)
hideydecorations!(ax2)

ln3 = lines!(ax3, analytic_κ²¹, z, color = :red)
sc3 = scatter!(ax3, κ²¹, z, color = :blue)
# axislegend(ax3, [ln3, sc3], [plot_string_1, plot_string_2], position = :rb, labelsize = 30)

ln4 = lines!(ax4, analytic_κ²², z, color = :red)
sc4 = scatter!(ax4, κ²², z, color = :blue)
hideydecorations!(ax4)
# axislegend(ax4, [ln4, sc4], [plot_string_1, plot_string_2], position = :rt, labelsize = 30)
# colgap!(ax1, 10)
# rowgap!(ax1, 10)

for ax in [ax1, ax2, ax3, ax4]
    # xlims!(ax, (-0.6, 0.6))
end

display(fig)

##
# Check x-independence of local diffusivity estimate
κ¹¹ = jlfile["localdiffusivity"]["κ11"]
κ¹² = jlfile["localdiffusivity"]["κ12"]
κ²¹ = jlfile["localdiffusivity"]["κ21"]
κ²² = jlfile["localdiffusivity"]["κ22"]
z = jlfile["grid"]["z"][:]
x = jlfile["grid"]["x"][:]
fig = Figure()
ax1 = Axis(fig[1, 1])
ax2 = Axis(fig[2, 1])
ax3 = Axis(fig[1, 2])
ax4 = Axis(fig[2, 2])
cmap = :balance
heatmap!(ax1, x, z, κ¹¹, colormap = cmap, interpolate = true)
heatmap!(ax2, x, z, κ²¹, colormap = cmap, interpolate = true)
heatmap!(ax3, x, z, κ¹², colormap = cmap, interpolate = true)
heatmap!(ax4, x, z, κ²², colormap = cmap, interpolate = true)
display(fig)

max_stars_repo_name: sandreza/StatisticalNonlocality

max_issues_repo_name: sandreza/StatisticalNonlocality

max_forks_repo_name: sandreza/StatisticalNonlocality


text:
using ReuseDistance

using Random
using Test

import ProgressMeter

include("redblack.jl")

@testset "Testing Node" begin
    Node = ReuseDistance.Node
    lchild = ReuseDistance.lchild
    rchild = ReuseDistance.rchild

    #####
    ##### rotate_right!
    #####

    # only left child
    node = Node(10)
    left = Node(5)
    ReuseDistance.setleft!(node, left)

    @test lchild(node) === left
    @test rchild(node) === nothing

    newnode = ReuseDistance.rotate_right!(node)
    @test newnode === left
    @test rchild(left) === node
    @test lchild(left) === nothing
    @test rchild(node) === nothing
    @test lchild(node) === nothing

    # two children
    node = Node(10)
    left = Node(5)
    right = Node(15)
    ReuseDistance.setleft!(node, left)
    ReuseDistance.setright!(node, right)

    @test lchild(node) === left
    @test rchild(node) === right

    newnode = ReuseDistance.rotate_right!(node)
    @test newnode === left
    @test rchild(left) === node
    @test lchild(left) === nothing
    @test rchild(node) === right
    @test lchild(node) === nothing

    # hierarchy
    node = Node(10)
    left = Node(5)
    leftleft = Node(1)
    leftright = Node(7)

    right = Node(20)

    ReuseDistance.setleft!(node, left)
    ReuseDistance.setleft!(left, leftleft)
    ReuseDistance.setright!(left, leftright)

    ReuseDistance.setright!(node, right)

    newnode = ReuseDistance.rotate_right!(node)
    @test newnode === left
    @test rchild(left) === node
    @test lchild(left) === leftleft

    @test rchild(node) === right
    @test lchild(node) === leftright

    #####
    ##### rotate_left!
    #####

    # only right child
    node = Node(10)
    right = Node(15)
    ReuseDistance.setright!(node, right)

    @test lchild(node) === nothing
    @test rchild(node) === right

    newnode = ReuseDistance.rotate_left!(node)
    @test newnode === right
    @test rchild(right) === nothing
    @test lchild(right) === node
    @test rchild(node) === nothing
    @test lchild(node) === nothing

    # two children
    node = Node(10)
    left = Node(5)
    right = Node(15)
    ReuseDistance.setleft!(node, left)
    ReuseDistance.setright!(node, right)

    @test lchild(node) === left
    @test rchild(node) === right

    newnode = ReuseDistance.rotate_left!(node)
    @test newnode === right
    @test rchild(right) === nothing
    @test lchild(right) === node
    @test rchild(node) === nothing
    @test lchild(node) === left

    # hierarchy
    node = Node(10)
    left = Node(5)

    right = Node(20)
    rightleft = Node(15)
    rightright = Node(25)

    ReuseDistance.setleft!(node, left)

    ReuseDistance.setright!(node, right)
    ReuseDistance.setleft!(right, rightleft)
    ReuseDistance.setright!(right, rightright)

    newnode = ReuseDistance.rotate_left!(node)
    @test newnode === right
    @test rchild(right) === rightright
    @test lchild(right) === node

    @test rchild(node) === rightleft
    @test lchild(node) === left
end

@testset "Testing Insert, Search, and Delete" begin
    treap = ReuseDistance.Treap{Int}()
    for i in shuffle(1:100)
        push!(treap, i)
    end

    @test all(in(treap), 1:100)
    @test !any(in(treap), 101:200)
    @test !any(in(treap), -100:0)
    @test ReuseDistance.isheap(treap)
    @test ReuseDistance.istree(treap)

    alldeleted = true
    for i in shuffle(2:2:100)
        alldeleted &= delete!(treap, i)
    end
    @test alldeleted
    @test length(treap) == 50
    @test all(in(treap), 1:2:99)
    @test !any(in(treap), 2:2:100)
    @test ReuseDistance.isheap(treap)
    @test ReuseDistance.istree(treap)

    anydeleted = false
    for i in shuffle(2:2:100)
        anydeleted |= delete!(treap, i)
    end
    @test anydeleted == false

    alldeleted = true
    for i in shuffle(1:2:99)
        alldeleted &= delete!(treap, i)
    end
    @test alldeleted
    @test length(treap) == 0
    @test !any(in(treap), 1:100)
    @test ReuseDistance.root(treap) === nothing
    @test ReuseDistance.isheap(treap)
    @test ReuseDistance.istree(treap)
end


text:
using VLMinorityGameKit
using BSON

# setup call -
memory = 10
numberOfStrategiesPerAgent = 5
numberOfAgents = 1001

# setup -
# path_to_game_world = "/Users/jeffreyvarner/Desktop/julia_work/VLMinorityGameKit.jl/test/data/GW-T-M10-S5-NA1001.bson"
path_to_game_world = "/Users/jeffreyvarner/Desktop/proposals/VLMinorityGameKit.jl/test/data/GW-GC-M10-S5-NA1001.bson"

# build the game world -
# args: numberOfAgents::Int64, agentMemorySize::Int64, numberOfStrategiesPerAgent::Int64; actions::Array{Int64,1} = [-1,0,1]
game_world = build_minority_game_world(VLGCMinorityGameWorld; numberOfAgents=numberOfAgents, memory=memory, numberOfStrategiesPerAgent=numberOfStrategiesPerAgent)

# write the gameworld to disk -
bson(path_to_game_world, gameworld=game_world)