text:
"Holds the tableau of an Specialised Partitioned Additive Runge-Kutta method."
const TableauSPARK = AbstractTableauSPARK{:spark}
const ParametersSPARK = AbstractParametersSPARK{:spark}

function RungeKutta.check_symplecticity(tab::TableauSPARK{T}; atol=16*eps(T), rtol=16*eps(T)) where {T}
    s_αa_qp̃ = [isapprox(tab.q.b[i] * tab.p.α[i,j] + tab.p.β[j] * tab.q̃.a[j,i], tab.q.b[i] * tab.p.β[j]; atol=atol, rtol=rtol) for i in 1:tab.s, j in 1:tab.r]
    s_α_q̃p̃  = [isapprox(tab.p.β[i] * tab.q̃.α[i,j] + tab.q.β[j] * tab.p̃.α[j,i], tab.p.β[i] * tab.q.β[j]; atol=atol, rtol=rtol) for i in 1:tab.r, j in 1:tab.r]
    s_b_qp  = isapprox.(tab.q.b, tab.p.b; atol=atol, rtol=rtol)
    s_ω     = [isapprox(tab.ω[i,j], (i == j ? 1 : 0); atol=atol, rtol=rtol) for i in 1:tab.s, j in 1:tab.s+1]

    return (s_αa_qp̃, s_α_q̃p̃, s_b_qp, s_ω)
end

function Integrators.symplecticity_conditions(::TableauSPARK)
    (
        """`` b^{1}_{i} b^{3}_{j} = b^{3}_{j} \\tilde{a}^{1}_{ji} + b^{1}_{i} a^{3}_{ij} ``""",
        """`` b^{2}_{i} b^{3}_{j} = b^{3}_{j} \\tilde{a}^{2}_{ji} + b^{2}_{i} \\tilde{a}^{3}_{ij} ``""",
        """`` b^{3}_{i} = b^{2}_{i} ``""",
        """`` \\omega_{ij} = \\begin{cases}  1 & i = j , \\\\  0 & i \\ne j . \\\\  \\end{cases} `` """,
    )
end

@doc raw"""
Specialised Partitioned Additive Runge-Kutta integrator for index-two DAE systems.

This integrator solves the following system of equations for the internal stages,
```math
\begin{aligned}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} U_{n,j} , & i &= 1, ..., s , \\
P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} G_{n,j} , & i &= 1, ..., s , \\
\tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} U_{n,j} , & i &= 1, ..., r , \\
\tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} G_{n,j} , & i &= 1, ..., r , \\
0 &= \sum \limits_{j=1}^{r} \omega_{ij} \tilde{\Phi}_{n,j} , & i &= 1, ..., r-1 ,
\end{aligned}
```
with definitions
```math
\begin{aligned}
P_{n,i} &= \frac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
F_{n,i} &= \frac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
U_{n,i} &= \hphantom{-} \frac{\partial \phi}{\partial p} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
G_{n,i} &=           -  \frac{\partial \phi}{\partial q} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
\tilde{\Phi}_{n,i} &= \phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) , & i &= 1, ..., r ,
\end{aligned}
```
and update rule
```math
\begin{aligned}
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} V_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} U_{n,i} , \\
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} F_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} G_{n,i} , \\
0 &= \phi (q_{n+1}, p_{n+1}) .
\end{aligned}
```
"""
struct IntegratorSPARK{DT, TT, D, S, R, PT <: ParametersSPARK{DT,TT,D,S,R},
                                        ST <: NonlinearSolver{DT},
                                        IT <: InitialGuessIODE{TT}} <: AbstractIntegratorSPARK{DT,TT,D,S,R}
    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    function IntegratorSPARK(params::ParametersSPARK{DT,TT,D,S,R}, solver::ST, iguess::IT, caches) where {DT,TT,D,S,R,ST,IT}
        new{DT, TT, D, S, R, typeof(params), ST, IT}(params, solver, iguess, caches)
    end

    function IntegratorSPARK{DT,D}(equations::NamedTuple, tableau::Union{TableauSPARK{TT},TableauSPARK{TT}}, Δt::TT) where {DT,TT,D}
        # @assert tableau.ρ == tableau.r-1

        # get number of stages
        S = tableau.s
        R = tableau.r
        P = tableau.ρ

        N = 2*D*S + 3*D*R

        if isdefined(tableau, :d) && length(tableau.d) > 0
            N += D
        end

        # create params
        params = ParametersSPARK{DT,D}(equations, tableau, Δt)

        # create cache dict
        caches = CacheDict(params)

        # create solver
        solver = create_nonlinear_solver(DT, N, params, caches)

        # create initial guess
        iguess = InitialGuessIODE(get_config(:ig_extrapolation), equations[:v̄], equations[:f̄], Δt)

        # create integrator
        IntegratorSPARK(params, solver, iguess, caches)
    end

    function IntegratorSPARK(equation::Union{IDAE{DT}, LDAE{DT}}, tableau::Union{TableauSPARK,TableauSPARK}, Δt; kwargs...) where {DT}
        IntegratorSPARK{DT, ndims(equation)}(get_functions(equation), tableau, Δt; kwargs...)
    end
end

GeometricBase.nconstraints(::IntegratorSPARK{DT,TT,D}) where {DT,TT,D} = D

function Integrators.initialize!(int::IntegratorSPARK, sol::AtomicSolutionPDAE)
    sol.t̄ = sol.t - timestep(int)

    equation(int, :v̄)(sol.t, sol.q, sol.v)
    equation(int, :f̄)(sol.t, sol.q, sol.v, sol.f)

    initialize!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f,
                            sol.t̄, sol.q̄, sol.p̄, sol.v̄, sol.f̄)
end

function initial_guess!(int::IntegratorSPARK{DT}, sol::AtomicSolutionPDAE{DT},
                        cache::IntegratorCacheSPARK{DT}=int.caches[DT]) where {DT}
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q̄, sol.p̄, sol.v̄, sol.f̄,
                              sol.q, sol.p, sol.v, sol.f,
                              cache.q̃, cache.p̃, cache.ṽ, cache.f̃,
                              tableau(int).q.c[i], tableau(int).p.c[i])

        for k in eachdim(int)
            int.solver.x[2*(ndims(int)*(i-1)+k-1)+1] = (cache.q̃[k] - sol.q[k])/timestep(int)
            int.solver.x[2*(ndims(int)*(i-1)+k-1)+2] = (cache.p̃[k] - sol.p[k])/timestep(int)
        end
    end

    for i in 1:pstages(int)
        evaluate!(int.iguess, sol.q̄, sol.p̄, sol.v̄, sol.f̄,
                              sol.q, sol.p, sol.v, sol.f,
                              cache.q̃, cache.p̃, cache.ṽ, cache.f̃,
                              tableau(int).q̃.c[i], tableau(int).p̃.c[i])

        for k in eachdim(int)
            int.solver.x[2*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+1] = (cache.q̃[k] - sol.q[k])/timestep(int)
            int.solver.x[2*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+2] = (cache.p̃[k] - sol.p[k])/timestep(int)
            int.solver.x[2*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+3] = cache.ṽ[k]
        end
    end

    # if isdefined(tableau(int), :λ) && tableau(int).λ.c[1] == 0
    #     for k in eachdim(int)
    #         int.solver.x[2*ndims(int)*nstages(int)+3*(k-1)+1] = cache.λ[k]
    #     end
    # end

    if isdefined(tableau(int), :d) && length(tableau(int).d) > 0
        for k in eachdim(int)
            int.solver.x[2*ndims(int)*nstages(int)+3*ndims(int)*pstages(int)+k] = 0
        end
    end
end

function compute_stages!(x::Vector{ST}, cache::IntegratorCacheSPARK{ST,D,S,R},
                                        params::ParametersSPARK{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    local tpᵢ::TT
    local tλᵢ::TT

    for i in 1:S
        for k in 1:D
            # copy x to Y, Z
            cache.Yi[i][k] = x[2*(D*(i-1)+k-1)+1]
            cache.Zi[i][k] = x[2*(D*(i-1)+k-1)+2]

            # compute Q and P
            cache.Qi[i][k] = params.q[k] + params.Δt * cache.Yi[i][k]
            cache.Pi[i][k] = params.p[k] + params.Δt * cache.Zi[i][k]
        end

        # compute f(X)
        tpᵢ = params.t + params.Δt * params.tab.p.c[i]
        params.equs[:f](tpᵢ, cache.Qi[i], cache.Vi[i], cache.Fi[i])
        params.equs[:ϑ](tpᵢ, cache.Qi[i], cache.Vi[i], cache.Φi[i])
        cache.Φi[i] .-= cache.Pi[i]
    end

    for i in 1:R
        for k in 1:D
            # copy y to Y, Z and Λ
            cache.Yp[i][k] = x[2*D*S+3*(D*(i-1)+k-1)+1]
            cache.Zp[i][k] = x[2*D*S+3*(D*(i-1)+k-1)+2]
            cache.Vp[i][k] = x[2*D*S+3*(D*(i-1)+k-1)+3]
            cache.Λp[i][k] = cache.Vp[i][k]

            # compute Q and V
            cache.Qp[i][k] = params.q[k] + params.Δt * cache.Yp[i][k]
            cache.Pp[i][k] = params.p[k] + params.Δt * cache.Zp[i][k]
        end

        # compute f(X)
        tλᵢ = params.t + params.Δt * params.tab.λ.c[i]
        params.equs[:u](tλᵢ, cache.Qp[i], cache.Pp[i], cache.Vp[i], cache.Up[i])
        params.equs[:g](tλᵢ, cache.Qp[i], cache.Pp[i], cache.Vp[i], cache.Gp[i])
        params.equs[:ϕ](tλᵢ, cache.Qp[i], cache.Pp[i], cache.Φp[i])
    end

    if isdefined(params.tab, :d) && length(params.tab.d) > 0
        for k in 1:D
            cache.μ[k] = x[2*D*S+3*D*R+k]
        end
    end

    # compute q and p
    cache.q̃ .= params.q
    cache.p̃ .= params.p
    for i in 1:S
        cache.p̃ .+= params.Δt .* params.tab.p.b[i] .* cache.Fi[i]
    end
    for i in 1:R
        cache.q̃ .+= params.Δt .* params.tab.q.β[i] .* cache.Up[i]
        cache.p̃ .+= params.Δt .* params.tab.p.β[i] .* cache.Gp[i]
    end

    # compute ϕ(q,p)
    tλᵢ = params.t + params.Δt
    params.equs[:ϕ](tλᵢ, cache.q̃, cache.p̃, cache.ϕ̃)
end

"Compute stages of specialised partitioned additive Runge-Kutta methods for variational systems."
function Integrators.function_stages!(y::Vector{ST}, b::Vector{ST}, params::ParametersSPARK{DT,TT,D,S,R,P},
                                      caches::CacheDict) where {ST,DT,TT,D,S,R,P}

    # get cache for internal stages
    cache = caches[ST]

    compute_stages!(y, cache, params)

    # compute b = - [(Y-AV-AU), (Z-AF-AG), Φ]
    for i in 1:S
        for k in 1:D
            b[2*(D*(i-1)+k-1)+1] = - cache.Yi[i][k]
            b[2*(D*(i-1)+k-1)+2] = - cache.Zi[i][k]
            for j in 1:S
                b[2*(D*(i-1)+k-1)+2] += params.tab.p.a[i,j] * cache.Fi[j][k]
            end
            for j in 1:R
                b[2*(D*(i-1)+k-1)+1] += params.tab.q.α[i,j] * cache.Up[j][k]
                b[2*(D*(i-1)+k-1)+2] += params.tab.p.α[i,j] * cache.Gp[j][k]
            end
        end
    end

    # compute b = - [(Y-AV-AU), (Z-AF-AG)]
    for i in 1:R
        for k in 1:D
            b[2*D*S+3*(D*(i-1)+k-1)+1] = - cache.Yp[i][k]
            b[2*D*S+3*(D*(i-1)+k-1)+2] = - cache.Zp[i][k]
            b[2*D*S+3*(D*(i-1)+k-1)+3] = 0
            for j in 1:S
                b[2*D*S+3*(D*(i-1)+k-1)+2] += params.tab.p̃.a[i,j] * cache.Fi[j][k]
            end
            for j in 1:R
                b[2*D*S+3*(D*(i-1)+k-1)+1] += params.tab.q̃.α[i,j] * cache.Up[j][k]
                b[2*D*S+3*(D*(i-1)+k-1)+2] += params.tab.p̃.α[i,j] * cache.Gp[j][k]
            end
        end
    end

    # compute b = - ωΦ
    for i in 1:R-P
        for k in 1:D
            for j in 1:R
                b[2*D*S+3*(D*(i-1)+k-1)+3] -= params.tab.ω[i,j] * cache.Φp[j][k]
            end
            b[2*D*S+3*(D*(i-1)+k-1)+3] -= params.tab.ω[i,R+1] * cache.ϕ̃[k]
        end
    end

    # compute b = d_λ ⋅ V
    for i in R-P+1:R
        for k in 1:D
            for j in 1:R
                b[2*D*S+3*(D*(i-1)+k-1)+3] -= params.tab.δ[j] * cache.Vp[j][k]
            end
        end
    end

    if isdefined(params.tab, :d) && length(params.tab.d) > 0
        for i in 1:R
            for k in 1:D
                b[2*(D*(i-1)+k-1)+3] -= cache.μ[k] * params.tab.d[i] / params.tab.p.b[i]
            end
        end

        for k in 1:D
            b[2*D*S+3*D*R+k] = 0
            for i in 1:R
                b[2*D*S+3*D*R+k] -= cache.Vp[i][k] * params.tab.d[i]
            end
        end
    end
end


text:
@recipe(FlexpartPlot) do scene
    Attributes(
        time = Ti(1),
        height = Dim{:height}(1),
        nageclass = Dim{:nageclass}(1),
        pointspec = Dim{:pointspec}(1),
        lons = X(:),
        lats = Y(:),
        scale = :lin,
        title = ""
    )
end

@recipe(SliderPlot) do scene
    Attributes(
        lons = X(:),
        lats = Y(:),
        scale = :lin,
        title = ""
    )
end

Makie.convert_arguments(P::Union{Type{<:Heatmap}, Type{<:Contourf}, Type{<:Contour}}, r::AbstractRaster) = convert_arguments(P, dims(r, :X) |> collect, dims(r, :Y) |> collect, Matrix(r))
# Makie.convert_single_argument(r::AbstractRaster) = dims(r, :X) |> collect, dims(r, :Y) |> collect, Matrix(r)
# Makie.convert_arguments(P::Type{<:Contourf}, r::AbstractRaster) = convert_arguments(P, dims(r, :X) |> collect, dims(r, :Y) |> collect, Matrix(r))

# _convert(r::AbstractRaster) = Vector(dims(r, :X) |> collect), Vector(dims(r, :Y) |> collect), Matrix(r)

# function Makie.plot!(ensplot::EnsemblePlot{<:Tuple{AbstractVector{<:EnsembleOutput}}})
function Makie.plot!(fpplot::FlexpartPlot{<:Tuple{<:MixingRatio}})
    # ax = Axis(ensplot)
    # global access = ensplot
    mr = fpplot[1].val
    raster = mr.raster
    # filename = convert(String, fpoutput)
    # raster = Raster(filename, name= :spec001_mr)
    scale =  fpplot.scale.val
    filtrast = view(raster,
        fpplot.time.val,
        fpplot.height.val,
        fpplot.nageclass.val,
        fpplot.pointspec.val,
        fpplot.lons.val,
        fpplot.lats.val,
        )
    lons = dims(filtrast, :X) |> collect
    lats = dims(filtrast, :Y) |> collect
    cmap = cgrad(:viridis, scale = scale)

    matrix = replace!(x -> isapprox(x, 0) ? NaN : x, Matrix(filtrast))
    scale == :log10 && replace!(x -> log10(x), matrix)
    heatmap!(fpplot, lons, lats,
        matrix,
        colormap = cmap,
        axis = (aspect = DataAspect(), title = fpplot.title, titlegap = 0),
        # cscale = log10
        )
    # Colorbar(ax[1, 2], heat)
end

# function Makie.plot!(sp::SliderPlot{<:Tuple{<:AbstractRaster}})
#     # ax = Axis(ensplot)
#     # global access = ensplot
#     rasterobs = sp[1]
#     raster = rasterobs[]
#     times = dims(raster, Ti) |> collect

#     mr = fpplot[1].val
#     raster = mr.raster
#     # filename = convert(String, fpoutput)
#     # raster = Raster(filename, name= :spec001_mr)
#     scale =  fpplot.scale.val
#     filtrast = view(raster,
#         fpplot.time.val,
#         fpplot.height.val,
#         fpplot.nageclass.val,
#         fpplot.pointspec.val,
#         fpplot.lons.val,
#         fpplot.lats.val,
#         )
#     lons = dims(filtrast, :X) |> collect
#     lats = dims(filtrast, :Y) |> collect
#     cmap = cgrad(:viridis, scale = scale)

#     matrix = replace!(x -> isapprox(x, 0) ? NaN : x, Matrix(filtrast))
#     scale == :log10 && replace!(x -> log10(x), matrix)
#     heatmap!(fpplot, lons, lats,
#         matrix,
#         colormap = cmap,
#         axis = (aspect = DataAspect(), title = fpplot.title, titlegap = 0),
#         # cscale = log10
#         )
#     # Colorbar(ax[1, 2], heat)
# end

function sliderplot(raster::AbstractRaster)
    fig = Figure()
    # ax = Axis(fig[1, 1])
    times = dims(raster, Ti) |> collect
    heights = dims(raster, Dim{:height}) |> collect
    lsgrid = labelslidergrid!(
        fig,
        ["Time", "Height"],
        [1:length(times), 1:length(heights)];
        # formats = "{:.1f}" .* ["V", "A", "Ω"],
        formats = [i -> Dates.format(times[i], "yyyy-mm-ddTHH:MM:SS"), i -> "$(heights[i])m"],
        # formats = [x->"$x"],
        width = 600,
        height = 100)

    fig[1, 1] = lsgrid.layout

    sliderobs = [s.value for s in lsgrid.sliders]
    viewed = map(sliderobs...) do t, h
        view(raster,
            Ti(t),
            Dim{:height}(h),
            )
    end

    ax, hm = heatmap(fig[2, 1], viewed, colormap = (:Spectral, 0.3))
    Colorbar(fig[2, 2], hm)
    fig
end

function sliderplot(raster::AbstractRaster, tiff::GeoTiff)
    fig = Figure()
    # ax = Axis(fig[1, 1])
    times = dims(raster, Ti) |> collect
    heights = dims(raster, Dim{:height}) |> collect
    lsgrid = labelslidergrid!(
        fig,
        ["Time", "Height"],
        [1:length(times), 1:length(heights)];
        # formats = "{:.1f}" .* ["V", "A", "Ω"],
        formats = [i -> Dates.format(times[i], "yyyy-mm-ddTHH:MM:SS"), i -> "$(heights[i])m"],
        # formats = [x->"$x"],
        width = 600,
        height = 100)

    fig[1, 1] = lsgrid.layout

    sliderobs = [s.value for s in lsgrid.sliders]
    viewed = map(sliderobs...) do t, h
        view(raster,
            Ti(t),
            Dim{:height}(h),
            )
    end

    ax = Axis(fig[2, 1])
    hmtiff = heatmap!(ax, tiff.raster)
    hm = heatmap!(ax, viewed, colormap = (:viridis, 0.8))
    Colorbar(fig[2, 2], hm)
    # translate!(hm, 0, 0, -1)
    fig
end

struct LogMinorTicks end

function MakieLayout.get_minor_tickvalues(
        ::LogMinorTicks, scale, tickvalues, vmin, vmax
)
    vals = Float64[]
    extended_tickvalues = [
        tickvalues[1] - (tickvalues[2] - tickvalues[1]);
        tickvalues;
        tickvalues[end] + (tickvalues[end] - tickvalues[end-1]);
    ]
    # @show extended_tickvalues
    # @show tickvalues
    for (lo, hi) in zip(
            @view(extended_tickvalues[1:end-1]),
            @view(extended_tickvalues[2:end])
        )
        interval = hi-lo
        steps = log10.(LinRange(10^lo, 10^hi, 11))
        append!(vals, steps[2:end-1])
    end
    return filter(x -> vmin < x < vmax, vals)
end

custom_formatter(values) = map(
    v -> "10" * Makie.UnicodeFun.to_superscript(round(Int64, v)),
    values
)

ext: jl


text:
function get_nullable{K, V}(dict :: Dict{K, V}, key :: K)
  if haskey(dict, key)
    Nullable{V}(dict[key])
  else
    Nullable{V}()
  end
end


text:
#
# Functions for handling unassigned shifts
#

# Use this when there are not enough employees (or no employees)
# so we don't have 80 loops of unassigned shifts
function meet_base_coverage(env, employees)
    employees = deepcopy(employees)
    days_per_week = size(env["coverage"], 1)
    sum_coverage = sum(sum(env["coverage"]))
    shift_count = 0

    # Sum the hours of employees
    employee_hour_sum = 0
    for e in keys(employees)
        employee_hour_sum += employees[e]["hours_max"]
    end

    # min_lift actually raises the minimum number of hours needed.
    # It helps feasibility, but could hurt the one case of a perfect, single-
    # solution optimal schedule.
    while employee_hour_sum < (sum_coverage*MIN_LYFT)
        name, value = generate_unassigned_shift(env)
        employees[name] = value
        employee_hour_sum += value["hours_max"]
        shift_count += 1
    end

    Logging.info( "Generated $shift_count shifts during pre-processing")
    return employees
end

# Generate lots of unassigned shifts - more than we need -
# so it is fast
function meet_unassigned_base_coverage(env, day)
    sum_coverage = sum(env["coverage"][day])

    average_unassigned_shift_length = env["shift_time_min"] + (env["shift_time_max"] - env["shift_time_min"]) * UNASSIGNED_RATIO
    unassigned_shift_count = ceil(sum_coverage / average_unassigned_shift_length)

    # Meet peak staffing level
    peak = 0
    for t in 1:length(env["coverage"][day])
        if env["coverage"][day][t] > peak
            peak = env["coverage"][day][t]
        end
    end
    if peak > unassigned_shift_count
        Logging.info("Peak coverage ($peak shifts) chosen instead for unassigned shift count ($unassigned_shift_count shifts) as starting point")
        unassigned_shift_count = peak
    end

    employees = Dict()
    for i in 1:unassigned_shift_count
        name, value = generate_unassigned_shift(env)
        employees[name] = value
    end

    Logging.info( "Generated $unassigned_shift_count unassigned shifts during pre-processing for day $day")
    return employees
end

function generate_unassigned_shift(env)
    shift_time_min = env["shift_time_min"]
    shift_time_max = env["shift_time_max"]
    shift_count = 1 # Each unassigned shift is a single shift

    # need unique names
    name = string(UNASSIGNED_PREFIX, randstring(8))

    # Available all the time
    availability = Array[]
    days_per_week = size(env["coverage"], 1)
    for day in 1:days_per_week
        push!(availability, [ones(Int, length(env["coverage"][day]))])
    end

    # return tuple so you can do employees[name] = values
    return name, {
        "hours_min" => shift_time_min,
        "shift_time_min" => shift_time_min,
        "hours_max" => shift_time_max,
        "shift_time_max" => shift_time_max,
        "shift_count_max" => shift_count,
        "shift_count_min" => shift_count,
        "availability" => availability,
    }
end

function is_unassigned_shift(name)
    # Julia doesn't wrap bounds :-(
    if length(UNASSIGNED_PREFIX) > length(name)
        return false
    end
    return name[1:length(UNASSIGNED_PREFIX)] == UNASSIGNED_PREFIX
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
using Flux
using Flux: crossentropy, throttle
using Flux.Data: Tree, children, isleaf

include("data.jl")

N = 300

embedding = param(randn(Float32, N, length(alphabet)))

W = Dense(2N, N, tanh)
combine(a, b) = W([a; b])

sentiment = Chain(Dense(N, 5), softmax)

function forward(tree)
  if isleaf(tree)
    token, sent = tree.value
    phrase = embedding * token
    phrase, crossentropy(sentiment(phrase), sent)
  else
    _, sent = tree.value
    c1, l1 = forward(tree[1])
    c2, l2 = forward(tree[2])
    phrase = combine(c1, c2)
    phrase, l1 + l2 + crossentropy(sentiment(phrase), sent)
  end
end

loss(tree) = forward(tree)[2]

opt = ADAM()
ps = params(embedding, W, sentiment)
evalcb = () -> @show loss(train[1])

Flux.train!(loss, ps, zip(train), opt,
           cb = throttle(evalcb, 10))

max_stars_repo_name: jstrube/model-zoo

max_issues_repo_name: jstrube/model-zoo

max_forks_repo_name: jstrube/model-zoo


text:
module Alert
using Base64
using Dates
using Logging
using Printf
using REPL
using MacroTools

export alert_REPL!, alert, @alert

function __init__()
    init_alert_backends()
    init_alert_REPL!()
end

include("core.jl")
include("at_alert.jl")
include("alert_repl.jl")
include("apple.jl")

end


text:
VcovFormula(::Type{Val{:simple}}) = VcovSimpleFormula()

struct VcovSimpleFormula <: AbstractVcovFormula end

struct VcovSimpleMethod <: AbstractVcovMethod end
VcovMethod(::AbstractDataFrame, ::VcovSimpleFormula) = VcovSimpleMethod()

function vcov!(::VcovSimpleMethod, x::VcovData)
    invcrossmatrix = Matrix(inv(x.crossmatrix))
    rmul!(invcrossmatrix, sum(abs2, x.residuals) /  x.dof_residual)
    return Symmetric(invcrossmatrix)
end
shat!(::VcovSimpleMethod, x::VcovData) = scale(x.crossmatrix, sumabs2(x.residuals))


text:
"""
    calibrate(d::Dict, set::Dict; save_build=true, overwrite=false)
    calibrate(year::Int, d::Dict, set::Dict)

# Arguments
- `d::Dict` of DataFrames containing the model data.
- `set::Dict` of Arrays describing region, sector, final demand, etc.
- `year::Int`: year for which to perform calibration

# Keywords
- `save_build = true`
- `overwrite = false`
See [`SLiDE.build`](@ref) for keyword argument descriptions.

# Returns
- `d::Dict` of DataFrames containing the model data at the calibration step.
"""
function calibrate(
    dataset::String,
    d::Dict,
    set::Dict;
    save_build::Bool=DEFAULT_SAVE_BUILD,
    overwrite::Bool=DEFAULT_OVERWRITE,
    penalty_nokey::AbstractFloat=DEFAULT_PENALTY_NOKEY,
)
    CURR_STEP = "calibrate"

    # If there is already calibration data, read it and return.
    d_read = read_build(dataset, CURR_STEP; overwrite = overwrite)
    !(isempty(d_read)) && (return d_read)

    # Copy the relevant input DataFrames before making any changes.
    _calibration_set!(set)
    d = Dict(k => copy(d[k]) for k in set[:cal])

    # Initialize a DataFrame to contain results and do the calibration iteratively.
    cal = Dict(k => DataFrame() for k in set[:cal])

    for year in set[:yr]
        cal_yr = calibrate(year, d, set; penalty_nokey = penalty_nokey)
        [cal[k] = [cal[k]; cal_yr[k]] for k in set[:cal]]
    end

    write_build!(dataset, CURR_STEP, cal; save_build = save_build)
    return cal
end

function calibrate(
    year::Int,
    io::Dict,
    set::Dict;
    penalty_nokey::AbstractFloat=DEFAULT_PENALTY_NOKEY,
)
    @info("Calibrating $year data")

    # Prepare the data and initialize the model.
    _calibration_set!(set)
    (cal, idx) = _calibration_input(year, io, set);
    calib = Model(optimizer_with_attributes(Ipopt.Optimizer, "max_cpu_time" => 60.0))

    @variable(calib, ys0_est[j in set[:j], i in set[:i]]   >= 0, start = 0);
    @variable(calib, fs0_est[i in set[:i]]                 >= 0, start = 0);
    @variable(calib, ms0_est[i in set[:i], m in set[:m]]   >= 0, start = 0);
    @variable(calib, y0_est[i in set[:i]]                  >= 0, start = 0);
    @variable(calib, id0_est[i in set[:i],j in set[:j]]    >= 0, start = 0);
    @variable(calib, fd0_est[i in set[:i], fd in set[:fd]] >= 0, start = 0);
    @variable(calib, va0_est[va in set[:va], j in set[:j]] >= 0, start = 0);
    @variable(calib, a0_est[i in set[:i]]                  >= 0, start = 0);
    @variable(calib, x0_est[i in set[:i]]                  >= 0, start = 0);
    @variable(calib, m0_est[i in set[:i]]                  >= 0, start = 0);
    @variable(calib, md0_est[m in set[:m], i in set[:i]]   >= 0, start = 0);

    # --- DEFINE CONSTRAINTS ---------------------------------------------------------------
    @constraint(calib,mkt_py[i in set[:i]],
        sum(ys0_est[j,i] for j in set[:j]) + fs0_est[i] ==
        sum(ms0_est[i,m] for m in set[:m]) + y0_est[i]
    );

    @constraint(calib,mkt_pa[i in set[:i]],
        a0_est[i] ==
        sum(id0_est[i,j] for j in set[:j]) + sum(fd0_est[i,fd] for fd in set[:fd])
    );

    @constraint(calib,mkt_pm[m in set[:m]],
        sum(ms0_est[i,m] for i in set[:i]) ==
        sum(md0_est[m,i] for i in set[:i])
    );

    @constraint(calib,prf_y[j in set[:j]],
        sum(ys0_est[j,i] for i in set[:i]) ==
        sum(id0_est[i,j] for i in set[:i]) + sum(va0_est[va,j] for va in set[:va])
    );

    @constraint(calib,prf_a[i in set[:i]],
        a0_est[i] * (1 - cal[:ta0][i]) + x0_est[i] ==
        y0_est[i] + m0_est[i] * (1 + cal[:tm0][i]) + sum(md0_est[m,i] for m in set[:m])
    );

    # --- DEFINE OBJECTIVE -----------------------------------------------------------------
    @objective(calib,Min,
        + sum(abs(cal[:ys0][j,i]) * (ys0_est[j,i] / cal[:ys0][j,i]   - 1)^2 for (j,i)  in set[:j,:i]  if cal[:ys0][j,i]  != 0)
        + sum(abs(cal[:id0][i,j]) * (id0_est[i,j] / cal[:id0][i,j]   - 1)^2 for (i,j)  in set[:i,:j]  if cal[:id0][i,j]  != 0)
        + sum(abs(cal[:fs0][i])   * (fs0_est[i]/ cal[:fs0][i]        - 1)^2 for  i     in set[:i]     if cal[:fs0][i]    != 0)
        + sum(abs(cal[:ms0][i,m]) * (ms0_est[i,m]/ cal[:ms0][i,m]    - 1)^2 for (i,m)  in set[:i,:m]  if cal[:ms0][i,m]  != 0)
        + sum(abs(cal[:y0][i])    * (y0_est[i]/ cal[:y0][i]          - 1)^2 for  i     in set[:i]     if cal[:y0][i]     != 0)
        + sum(abs(cal[:fd0][i,fd])* (fd0_est[i,fd] / cal[:fd0][i,fd] - 1)^2 for (i,fd) in set[:i,:fd] if cal[:fd0][i,fd] != 0)
        + sum(abs(cal[:va0][va,j])* (va0_est[va,j] / cal[:va0][va,j] - 1)^2 for (va,j) in set[:va,:j] if cal[:va0][va,j] != 0)
        + sum(abs(cal[:a0][i])    * (a0_est[i] / cal[:a0][i]         - 1)^2 for  i     in set[:i]     if cal[:a0][i]     != 0)
        + sum(abs(cal[:x0][i])    * (x0_est[i] / cal[:x0][i]         - 1)^2 for  i     in set[:i]     if cal[:x0][i]     != 0)
        + sum(abs(cal[:m0][i])    * (m0_est[i] / cal[:m0][i]         - 1)^2 for  i     in set[:i]     if cal[:m0][i]     != 0)
        + sum(abs(cal[:md0][m,i]) * (md0_est[m,i] / cal[:md0][m,i]   - 1)^2 for (m,i)  in set[:m,:i]  if cal[:md0][m,i]  != 0)

    + penalty_nokey * (
        + sum(ys0_est[j,i]  for (j,i)  in set[:j,:i]  if cal[:ys0][j,i]  == 0)
        + sum(id0_est[i,j]  for (i,j)  in set[:i,:j]  if cal[:id0][i,j]  == 0)
        + sum(fs0_est[i]    for  i     in set[:i]     if cal[:fs0][i]    == 0)
        + sum(ms0_est[i,m]  for (i,m)  in set[:i,:m]  if cal[:ms0][i,m]  == 0)
        + sum(y0_est[i]     for  i     in set[:i]     if cal[:y0][i]     == 0)
        + sum(fd0_est[i,fd] for (i,fd) in set[:i,:fd] if cal[:fd0][i,fd] == 0)
        + sum(va0_est[va,j] for (va,j) in set[:va,:j] if cal[:va0][va,j] == 0)
        + sum(a0_est[i]     for  i     in set[:i]     if cal[:a0][i]     == 0)
        + sum(x0_est[i]     for  i     in set[:i]     if cal[:x0][i]     == 0)
        + sum(m0_est[i]     for  i     in set[:i]     if cal[:m0][i]     == 0)
        + sum(md0_est[m,i]  for (m,i)  in set[:m,:i]  if cal[:md0][m,i]  == 0)
        )
    );

    # --- SET START VALUE ------------------------------------------------------------------
    [set_start_value(ys0_est[j,i], cal[:ys0][j,i])  for (j,i)  in set[:j,:i]  ]
    [set_start_value(id0_est[i,j], cal[:id0][i,j])  for (i,j)  in set[:i,:j]  ]
    [set_start_value(fs0_est[i],   cal[:fs0][i])    for  i     in set[:i]     ]
    [set_start_value(ms0_est[i,m], cal[:ms0][i,m])  for (i,m)  in set[:i,:m]  ]
    [set_start_value(y0_est[i],    cal[:y0][i])     for  i     in set[:i]     ]
    [set_start_value(fd0_est[i,fd],cal[:fd0][i,fd]) for (i,fd) in set[:i,:fd] ]
    [set_start_value(va0_est[va,j],cal[:va0][va,j]) for (va,j) in set[:va,:j] ]
    [set_start_value(a0_est[i],    cal[:a0][i])     for  i     in set[:i]     ]
    [set_start_value(x0_est[i],    cal[:x0][i])     for  i     in set[:i]     ]
    [set_start_value(m0_est[i],    cal[:m0][i])     for  i     in set[:i]     ]
    [set_start_value(md0_est[m,i], cal[:md0][m,i])  for (m,i)  in set[:m,:i]  ]

    # --- SET BOUNDS -----------------------------------------------------------------------
    # multipliers for lower and upper bound relative
    # to each respective variables reference parameter
    lb = DEFAULT_CALIBRATE_LOWER_BOUND
    ub = DEFAULT_CALIBRATE_UPPER_BOUND

    [set_lower_bound(ys0_est[j,i], max(0, lb * cal[:ys0][j,i]))  for (j,i)  in set[:j,:i]  ]
    [set_lower_bound(id0_est[i,j], max(0, lb * cal[:id0][i,j]))  for (i,j)  in set[:i,:j]  ]
    [set_lower_bound(fs0_est[i],   max(0, lb * cal[:fs0][i]))    for  i     in set[:i]     ]
    [set_lower_bound(ms0_est[i,m], max(0, lb * cal[:ms0][i,m]))  for (i,m)  in set[:i,:m]  ]
    [set_lower_bound(y0_est[i],    max(0, lb * cal[:y0][i]))     for  i     in set[:i]     ]
    [set_lower_bound(fd0_est[i,fd],max(0, lb * cal[:fd0][i,fd])) for (i,fd) in set[:i,:fd] ]
    [set_lower_bound(va0_est[va,j],max(0, lb * cal[:va0][va,j])) for (va,j) in set[:va,:j] ]
    [set_lower_bound(a0_est[i],    max(0, lb * cal[:a0][i]))     for  i     in set[:i]     ]
    [set_lower_bound(x0_est[i],    max(0, lb * cal[:x0][i]))     for  i     in set[:i]     ]
    [set_lower_bound(m0_est[i],    max(0, lb * cal[:m0][i]))     for  i     in set[:i]     ]
    [set_lower_bound(md0_est[m,i], max(0, lb * cal[:md0][m,i]))  for (m,i)  in set[:m,:i]  ]

    [set_upper_bound(ys0_est[j,i], abs(ub * cal[:ys0][j,i]))  for (i,j)  in set[:i,:j]  ]
    [set_upper_bound(id0_est[i,j], abs(ub * cal[:id0][i,j]))  for (i,j)  in set[:i,:j]  ]
    [set_upper_bound(fs0_est[i],   abs(ub * cal[:fs0][i]))    for  i     in set[:i]     ]
    [set_upper_bound(ms0_est[i,m], abs(ub * cal[:ms0][i,m]))  for (i,m)  in set[:i,:m]  ]
    [set_upper_bound(y0_est[i],    abs(ub * cal[:y0][i]))     for  i     in set[:i]     ]
    [set_upper_bound(fd0_est[i,fd],abs(ub * cal[:fd0][i,fd])) for (i,fd) in set[:i,:fd] ]
    [set_upper_bound(va0_est[va,j],abs(ub * cal[:va0][va,j])) for (va,j) in set[:va,:j] ]
    [set_upper_bound(a0_est[i],    abs(ub * cal[:a0][i]))     for  i     in set[:i]     ]
    [set_upper_bound(x0_est[i],    abs(ub * cal[:x0][i]))     for  i     in set[:i]     ]
    [set_upper_bound(m0_est[i],    abs(ub * cal[:m0][i]))     for  i     in set[:i]     ]
    [set_upper_bound(md0_est[m,i], abs(ub * cal[:md0][m,i]))  for (m,i)  in set[:m,:i]  ]

    # Fix "certain parameters" to their original values: fs0, va0, m0.
    [fix(fs0_est[i],    cal[:fs0][i],    force=true) for  i     in set[:i]     ]
    [fix(va0_est[va,j], cal[:va0][va,j], force=true) for (va,j) in set[:va,:j] ]
    [fix(m0_est[i],     cal[:m0][i],     force=true) for  i     in set[:i]     ]

    # Fix other/use sector output to zero.
    [fix(ys0_est[j,i], 0, force = true) for j in set[:oth,:use] for i in set[:i]]

    # --- OPTIMIZE AND SAVE RESULTS --------------------------------------------------------
    JuMP.optimize!(calib)

    # Populate resultant Dictionary.
    cal = Dict(k => filter_with(io[k], (yr = year,); drop = true) for k in [:ta0,:tm0])
    cal[:ys0] = convert_type(DataFrame, ys0_est; cols=idx[:ys0])
    cal[:fs0] = convert_type(DataFrame, fs0_est; cols=idx[:fs0])
    cal[:ms0] = convert_type(DataFrame, ms0_est; cols=idx[:ms0])
    cal[:y0]  = convert_type(DataFrame, y0_est;  cols=idx[:y0])
    cal[:id0] = convert_type(DataFrame, id0_est; cols=idx[:id0])
    cal[:fd0] = convert_type(DataFrame, fd0_est; cols=idx[:fd0])
    cal[:va0] = convert_type(DataFrame, va0_est; cols=idx[:va0])
    cal[:a0]  = convert_type(DataFrame, a0_est;  cols=idx[:a0])
    cal[:x0]  = convert_type(DataFrame, x0_est;  cols=idx[:x0])
    cal[:m0]  = convert_type(DataFrame, m0_est;  cols=idx[:m0])
    cal[:md0] = convert_type(DataFrame, md0_est; cols=idx[:md0])

    # Add the year back to the DataFrame and return.
    x = Add(:yr, year)
    [cal[k] = edit_with(df, x)[:, [:yr; idx[k]; :value]] for (k, df) in cal]
    return cal
end

"""
    _calibration_input(year::Int, d::Dict, set::Dict)
This function prepares the input for the calibration routine:
    1. Select parameters relevant to the calibration routine.
    2. For all parameters except taxes (ta0, tm0), set negative values to zero.
        In the case of final demand, only set negative values to zero for
        `fd = pce`.
    3. Fill all "missing" values with zeros to generate a complete dataset. This is relevant
        to how the penalty for missing keys is applied in the objective function.

# Arguments
- `yr::Int`
- `d::Dict{Symbol,DataFrame}` of input DataFrames.

# Returns
- `cal::Dict{Symbol,Dict}`
"""
function _calibration_input(year::Int, d::Dict{Symbol,DataFrame}, set::Dict)
    param = Dict()

    param[:cal] = set[:cal]
    param[:tax] = [:ta0, :tm0]
    param[:var] = setdiff(param[:cal], param[:tax])

    # Isolate the current year.
    d = Dict(k => fill_zero(filter_with(d[k], (yr=year,); drop=true); with=set)
        for k in param[:cal])

    # Save the DataFrame indices in correct order. This will be used to convert the
    # calibration output back into DataFrames.
    idx = Dict(k => findindex(df) for (k,df) in d)

    # Set negative values to zero. In the case of final demand,
    # only set negative values to zero for fd = pce.
     d[:fd0][.&(d[:fd0][:,:fd] .== "pce", d[:fd0][:,:value] .< 0),:value] .= 0.0
     d[:fd0][.&(d[:fd0][:,:fd] .== "C", d[:fd0][:,:value] .< 0),:value] .= 0.0
    [d[k][d[k][:,:value] .< 0, :value] .= 0.0 for k in param[:var] if k !== :fd0]

    # Convert the calibration input into DataFrames. Fill zeros for taxes; drop otherwise.
    cal = Dict(k => convert_type(Dict, d[k]) for k in param[:cal])
    return (cal, idx)
end

"""
"""
function _calibration_set!(set)
    # !!!! replace in usage.
    !(:i in keys(set)) && (set[:i] = set[:g])
    !(:j in keys(set)) && (set[:j] = set[:s])

    !((:i,:j)  in keys(set)) && add_permutation!(set, (:i,:j))
    !((:j,:i)  in keys(set)) && add_permutation!(set, (:j,:i))
    !((:i,:m)  in keys(set)) && add_permutation!(set, (:i,:m))
    !((:m,:i)  in keys(set)) && add_permutation!(set, (:m,:i))
    !((:i,:fd) in keys(set)) && add_permutation!(set, (:i,:fd))
    !((:va,:j) in keys(set)) && add_permutation!(set, (:va,:j))

    if !(:cal in keys(set))
        set[:cal] = Symbol.(read_file([SLIDE_DIR,"src","build","parameters"],
            SetInput("list_calibrate.csv", :cal)))
    end
    return set
end

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
#!/usr/bin/env julia6
using ReadGlobal, MyStats

function probdf(file::String)
    writedlm(file*".pdf.txt",hist(readfield(file)))
    return 0
end

function probdf(file::String,condfile::String)
    N = parse(Int,split(condfile,"_")[end][2:end])
    writedlm(file*".cond."*condfile*".pdf.txt",hist(readfield(file),read(condfile,Int64,(N,))))
    return 0
end

probdf(ARGS...)


text:
using Test

@testset "ITensors.ContractionSequenceOptimization" begin
  @testset "$filename" for filename in ["itensor_contract.jl"]
    println("Running $filename")
    include(filename)
  end
end

nothing


text:
module TestFunctions

using Test

using MultivariatePolynomials
const MP = MultivariatePolynomials

using JuMP
using PolyJuMP

function test_functions(var)
    x = MP.similarvariable(var, Val{:x})
    α = MOI.VariableIndex(1)
    @test α * x isa AbstractTerm{MOI.VariableIndex}
    @test x * α isa AbstractTerm{MOI.VariableIndex}
    @test α * x^2 isa AbstractTerm{MOI.VariableIndex}
    @test x^2 * α isa AbstractTerm{MOI.VariableIndex}
    @test (1α) * x^2 isa AbstractTerm{MOI.ScalarAffineFunction{Int}}
    @test (1x)^2 * α isa AbstractTerm{MOI.ScalarAffineFunction{Int}}
    @test α * (1x) isa AbstractTerm{MOI.ScalarAffineFunction{Int}}
    @test (1x) * α isa AbstractTerm{MOI.ScalarAffineFunction{Int}}
    @test α * (1.0x) isa AbstractTerm{MOI.ScalarAffineFunction{Float64}}
    @test (1.0x) * α isa AbstractTerm{MOI.ScalarAffineFunction{Float64}}
    @test α * (1x + 2) isa AbstractPolynomial{MOI.ScalarAffineFunction{Int}}
    @test (1x + 2) * α isa AbstractPolynomial{MOI.ScalarAffineFunction{Int}}
    @test α * (1.0x + 2.0) isa AbstractPolynomial{MOI.ScalarAffineFunction{Float64}}
    @test (1.0x + 2.0) * α isa AbstractPolynomial{MOI.ScalarAffineFunction{Float64}}
end

function runtests(var)
    for name in names(@__MODULE__; all = true)
        if startswith("$name", "test_")
            @testset "$(name) $(typeof(var))" begin
                getfield(@__MODULE__, name)(var)
            end
        end
    end
end

end

import DynamicPolynomials
DynamicPolynomials.@polyvar(x)
TestFunctions.runtests(x)
import TypedPolynomials
TypedPolynomials.@polyvar(y)
TestFunctions.runtests(y)


text:
"""
*Types and constructors that represent Bell scenarios, their statistics, and their bounds.*

The BellScenario.jl module is the base library for Bell scenario analysis.

### Features:
* [`Scenario`](@ref): abstract type describing Bell scenarios.
* [`AbstractStrategy`](@ref): abstract type for Bell scenario statistics.
* [`AbstractGame`](@ref): abstract type for Bell scenario bounds.

A game theoretic framework is used to evaluate the performance of black-box systems.
A cost function testing the black-box system is regarded as a *game* while the statistics
generated by the black-box system are regarded as a *strategy*.
A strategy is played against a game to achieve a score, hence, this framework
provides a quantitative metric of the performance of a black-box system
in regards to a particular task.
"""
module BellScenario

using QBase, LinearAlgebra, Combinatorics

import Base: *, convert

# exported modules
export LocalPolytope, Nonlocality

# Legacy modules
export ConvexPolytope, Degeneracy, Behavior, QuantumBehavior, QuantumOpt

# Scenarios
include("./scenarios.jl")

# combinatorics math utilities
include("./combinatorics.jl")

# Strategies
include("./strategies.jl")
include("./deterministic_strategies.jl")
include("./strategy_conversions.jl")
include("./quantum_strategies.jl")

# Games
include("./games.jl")
include("./game_conversions.jl")

# read/write and printing
include("./file_io.jl")

include("./LocalPolytope/LocalPolytope.jl")
using .LocalPolytope

include("./Nonlocality/Nonlocality.jl")
using .Nonlocality

# include legacy modules
include("./Legacy/ConvexPolytope.jl")
using .ConvexPolytope

include("./Legacy/Behavior.jl")
using .Behavior

include("./Legacy/Degeneracy.jl")
using .Degeneracy

include("./Legacy/DichotomicLocalPolytope.jl")
using .DichotomicLocalPolytope

include("./Legacy/QuantumBehavior.jl")
include("./Legacy/QuantumOpt.jl")
include("./Legacy/LocalSignaling.jl")

end


text:
const A9 = 1
const B9 = 2
const C9 = 3
const D9 = 4
const E9 = 5
const F9 = 6
const G9 = 7
const H9 = 8
const I9 = 9
const A8 = 10
const B8 = 11
const C8 = 12
const D8 = 13
const E8 = 14
const F8 = 15
const G8 = 16
const H8 = 17
const I8 = 18
const A7 = 19
const B7 = 20
const C7 = 21
const D7 = 22
const E7 = 23
const F7 = 24
const G7 = 25
const H7 = 26
const I7 = 27
const A6 = 28
const B6 = 29
const C6 = 30
const D6 = 31
const E6 = 32
const F6 = 33
const G6 = 34
const H6 = 35
const I6 = 36
const A5 = 37
const B5 = 38
const C5 = 39
const D5 = 40
const E5 = 41
const F5 = 42
const G5 = 43
const H5 = 44
const I5 = 45
const A4 = 46
const B4 = 47
const C4 = 48
const D4 = 49
const E4 = 50
const F4 = 51
const G4 = 52
const H4 = 53
const I4 = 54
const A3 = 55
const B3 = 56
const C3 = 57
const D3 = 58
const E3 = 59
const F3 = 60
const G3 = 61
const H3 = 62
const I3 = 63
const A2 = 64
const B2 = 65
const C2 = 66
const D2 = 67
const E2 = 68
const F2 = 69
const G2 = 70
const H2 = 71
const I2 = 72
const A1 = 73
const B1 = 74
const C1 = 75
const D1 = 76
const E1 = 77
const F1 = 78
const G1 = 79
const H1 = 80
const I1 = 81

max_stars_repo_name: kimrin/JapaneseChess

max_issues_repo_name: kimrin/MechajyoForWindows

max_forks_repo_name: kimrin/MechajyoForWindows


text:
# This file was formerly a part of Julia. License is MIT: https://julialang.org/license

x="println(\"# This file was formerly a part of Julia. License is MIT: https://julialang.org/license\\n\\nx=\$(repr(x))\\n\$x\")"
println("# This file was formerly a part of Julia. License is MIT: https://julialang.org/license\n\nx=$(repr(x))\n$x")

max_stars_repo_name: vtjnash/Examples

max_issues_repo_name: JuliaArchive/Examples

max_forks_repo_name: JuliaArchive/Examples


text:
module SplineMakie

using GeometryBasics
using GeometryBasics: VecTypes
using Makie
using LinearAlgebra: Tridiagonal

export BezierPath, MoveTo, LineTo, CurveTo, EllipticalArc, ClosePath, discretize

include("bezier.jl")
include("bezier_extension.jl")

Makie.convert_arguments(P::PointBased, x::BezierPath) = (discretize(x),)

@recipe(Spline, positions) do scene
    Attributes(default_theme(scene, Lines)...)
end

Makie.conversion_trait(::Type{<:Spline}) = PointBased()

function Makie.plot!(p::Combined{spline, S}) where {S<:Tuple{AbstractVector{<:AbstractPoint}}}
    path = @lift splinepath($(p.positions))
    lines!(p, path; p.attributes...)
    plot
end

function splinepath(P::AbstractVector{<:AbstractPoint{Dim, T}}) where {Dim, T}
    if length(P) == 2
        return BezierPath([MoveTo(P[1]), LineTo(P[2])])
    end
    N = length(P)
    pxy = Tuple(cubic_spline(map(p -> p[i], P)) for i in 1:Dim)
    WP = Point{Dim, T}.(zip(pxy...))

    commands = Vector{PathCommand}(undef, N)
    commands[1] = MoveTo(P[1])
    for i in 2:(N-1)
        commands[i] = CurveTo(WP[i-1],
                              2*P[i] - WP[i],
                              P[i])
    end
    commands[N] = CurveTo(WP[N-1],
                          (P[N] + WP[N-1])/2,
                          P[N])
    BezierPath(commands)
end

function cubic_spline(p::Vector{<:Number})
    N = length(p) - 1

    dl = [if i==N-1; 2; else 1 end for i in 1:(N-1)]
    d  = [if i==1; 2; elseif i==N; 7; else 4 end for i in 1:N]
    du = [1 for i in 1:(N-1)]
    M = Tridiagonal(dl, d, du)

    b = [if i == 1
             p[i] + 2p[i+1]
         elseif i == N
             8p[i] + p[i+1]
         else
             4p[i] + 2p[i+1]
         end for i in 1:N]
    return M \ b
end

# TODO: add other inter types like hobby splines

end


text:
# Observer tests

# Test constructors
@test size(LinearObserver(ones(2,2)).H) == (2,2)
@test size(LinearObserver(1.0).H) == (1,1)
@test size(LinearObserver(1.0, 1.0).R) == (1,1)
@test size(LinearObserver(eye(3), eye(3)).H) == (3,3)
@test size(LinearObserver(ones(3)).R) == (3,3)
@test size(LinearObserver(ones(3), eye(3)).H) == (3,1)
@test size(LinearObserver(ones(3)').R) == (1,1)
@test size(LinearObserver(ones(3)', 1.0).R) == (1,1)
@test size(LinearObserver(ones(3)', ones(1,1)').R) == (1,1)
@test size(LinearObserver(ones(3)', 0.1).H) == (1,3)

# Test compatibility methods
lin_obs = LinearObserver(ones(2,2), eye(2))
StateEstimation.assert_compatibility(lin_obs, DiscreteState(ones(2)))
StateEstimation.assert_compatibility(DiscreteState(ones(2)), lin_obs)
@test_throws DimensionMismatch StateEstimation.assert_compatibility(lin_obs,
    DiscreteState(ones(3)))
lin_obs = LinearObserver(ones(2)', 2.0)
StateEstimation.assert_compatibility(DiscreteState([1.0]), lin_obs)
@test_throws DimensionMismatch StateEstimation.assert_compatibility(
    DiscreteState(ones(2)), lin_obs)

# Test discrete prediction methods
lin_obs = LinearObserver(ones(2,2), eye(2))
@test predict(DiscreteState(ones(2)), lin_obs) == DiscreteState(2*ones(2), 0)
@test predict(UncertainDiscreteState(ones(2), eye(2)), lin_obs) ==
    UncertainDiscreteState(2*ones(2), lin_obs.H*lin_obs.H'+lin_obs.R, 0)

# Test continuous prediction methods
lin_obs = LinearObserver(ones(2)', 1.0)
@test predict(ContinuousState(ones(2)), lin_obs) == ContinuousState(2.0, 0.0)
@test predict(UncertainContinuousState(ones(2), eye(2)), lin_obs) ==
    UncertainContinuousState(2.0*ones(1), lin_obs.H*lin_obs.H'+lin_obs.R, 0.0)

# Test observability methods
linear_sys = LinearSystem(Float64[[1, -2] [-3, -4]]', 0.1*eye(2))
linear_obs = LinearObserver(reshape(Float64[1, 2], 1, 2), reshape([0.1], 1, 1))
@test observable(linear_sys, linear_obs) == false
@test observable(eye(2), eye(2))

# Test simulation methods
lin_sys = LinearSystem([[0.0, 1.0] [-1.0, 1.0]], 2*eye(2))
lin_obs = LinearObserver(ones(2)', 1.0)
continuous_state = ContinuousState(ones(2))
uncertain_discrete_state = UncertainDiscreteState(ones(2), eye(2))
@test typeof(simulate(continuous_state, lin_obs)) <: ContinuousState
@test typeof(simulate(uncertain_discrete_state, lin_obs)) <: DiscreteState
sim_state, sim_meas = simulate(continuous_state, lin_sys, lin_obs, 1.0)
@test typeof(sim_state) <: ContinuousState
@test typeof(sim_meas) <: ContinuousState

# Test nonlinear constructors
discrete_nl_fcn(t, x::Vector) = x.^2;
discrete_nl_jac(t, x::Vector) = diagm(2*x)
NonlinearObserver(discrete_nl_fcn, discrete_nl_jac, eye(2))
NonlinearObserver(discrete_nl_fcn, discrete_nl_jac, 1.0)
NonlinearObserver(LinearObserver(ones(2,2), eye(2)))

# Test nonlinear absolute state prediction methods
nonlin_obs = NonlinearObserver(discrete_nl_fcn, discrete_nl_jac, eye(3))
@test predict(DiscreteState([0.0, 1.0, 2.0]), nonlin_obs) ==
    DiscreteState([0.0, 1.0, 4.0], 0)
@test predict(ContinuousState([0.0, 2.0, 4.0], 1.0), nonlin_obs) ==
    ContinuousState([0.0, 4.0, 16.0], 1.0)

# Test nonlinear uncertain state prediction methods
@test predict(UncertainDiscreteState([0.0, 1.0, 2.0], eye(3)), nonlin_obs) ==
    UncertainDiscreteState([0.0, 1.0, 4.0], diagm([1.0, 5.0, 17.0]), 0)
@test predict(UncertainContinuousState([0.0, 1.0, 2.0], eye(3)), nonlin_obs) ==
    UncertainContinuousState([0.0, 1.0, 4.0], diagm([1.0, 5.0, 17.0]), 0.0)


text:
"""
    arc(origin, radius, start_angle, stop_angle; kwargs...)

This function plots a circular arc, centered at `origin` with radius `radius`,
from `start_angle` to `stop_angle`.
`origin` must be a coordinate in 2 dimensions (i.e., a `Point2`); the rest of the arguments must be
`<: Number`.

Examples:

`arc(Point2f0(0), 1, 0.0, π)`
`arc(Point2f0(1, 2), 0.3. π, -π)`

## Attributes
$(ATTRIBUTES)
"""
@recipe(Arc, origin, radius, start_angle, stop_angle) do scene
    Attributes(;
        default_theme(scene, Lines)...,
        resolution = 361,
    )
end

function plot!(p::Arc)
    args = getindex.(p, (:origin, :radius, :start_angle, :stop_angle, :resolution))
    positions = lift(args...) do origin, radius, start_angle, stop_angle, resolution
        map(range(start_angle, stop=stop_angle, length=resolution)) do angle
            origin .+ Point2f0((cos(angle), sin(angle)) .* radius)
        end
    end
    lines!(p, Attributes(p), positions)
end


text:
# scratchwork!
function initialize_interaction_mat( h::Array{Int64, 2} )
    # return a matrix of interactions indicating who interacts with whom
    N, M = size(h)
    interaction_mat = zeros(Int64, N, N, M) # initialize
    for k in 1:M
        for i in 1:N
            if h[i,k] in [-1,1] # check if i belongs to k
                for j in i+1:N
                    if h[j,k] in [-1,1] # check if j both belong to set k
                        interaction_mat[i,j,k] = 1
                        interaction_mat[j,i,k] = 1
                    end
                end
            end
        end
    end
    return interaction_mat
end

function update_interaction_mat!( pop::Population )
    # function to upddate interaction_mat after each round
    if pop.sets_changed
        pop.sets.interaction_mat .= 0
        for k = 1:M
            for i = 1:N
                if h[i,k] in [-1,1] # check if i belongs to k
                    for j = i+1:N
                        if h[j,k] in [-1,1] # check if j both belong to set k
                            pop.sets.interaction_mat[i,j,k] = 1
                            pop.sets.interaction_mat[j,i,k] = 1
                        end
                    end
                end
            end
        end
    end
end

function update_interaction_mat2!( pop::Population )
    # function to upddate interaction_mat after each round
    if pop.sets_changed
        # only update interaction matrix if a set change has occurred
        set_members = [Int64[] for k in 1:M]
        for indices in findall(x->x in [1,-1], pop.sets.h)
            push!(set_members[indices[2]],indices[1])
        end

        # update interaction matrix
        pop.sets.interaction_mat .= 0
        for k = 1:M
            for i in set_members[k], j in set_members[k]
                if i < j
                    pop.sets.interaction_mat[i,j,k] = 1
                    pop.sets.interaction_mat[j,i,k] = 1
                end
            end
        end
    end
end

for k in 1:pop.sets.M  # for each set
    for (i, j) in pop.sets.set_pairs[k]   # for each pair within each set
        compute_payoffs!( pop, i, j, k )  # compute payoffs for i & j when intereacting in set k
                                          # this also updates pop.prev_actions accordingly
    end
end

for k in 1:pop.sets.M
    for i in 1:pop.sets.N
        for indices in findall(x->x==1, pop.sets.interaction_mat[i,:,k])
            # fill in here
        end
    end
end

function update_set_members_and_pairs!( sets::Sets )
    # return a list of list of set members given a boolean array of memberships
    # tried this, but it's not faster than set_members_and_pairs

    sets.set_members = [Int64[] for k in 1:pop.sets.M]
    for indices in findall(x->x in [1,-1], pop.sets.h)
        push!(sets.set_members[indices[2]],indices[1])
    end

    sets.set_pairs = [Tuple{Int64, Int64}[] for k in 1:pop.sets.M]
    for k in 1:M
        for i in sets.set_members[k]
            for j in sets.set_members[k]
                if i < j push!(sets.set_pairs[k], (i, j)) end
            end
        end
    end
end


text:
@testset "NumFieldOrd" begin
  x = PolynomialRing(FlintQQ, "x", cached = false)[2]
  K, a = number_field(x^2+1, check = false, cached = false)
  Kns, gKns = number_field([x^2+1], check = false, cached = false)
  Kt, t = PolynomialRing(K, "t", cached = false)
  L, b = number_field(t^2+3, "b", cached = false, check = false)
  Lns, gLns = number_field([t^2+3], cached = false, check = false)

  OK = maximal_order(K)
  @inferred nf(OK)
  @test issimple(OK)
  @test iscommutative(OK)
  @test @inferred isequation_order(OK)
  @test @inferred ismaximal(OK)
  @test @inferred degree(OK) == 2
  @test @inferred absolute_degree(OK) == 2
  BOK = absolute_basis(OK)
  @test @inferred discriminant(BOK) == @inferred discriminant(OK)
  @test @inferred signature(OK) == (0, 1)
  @test @inferred discriminant(OK, FlintQQ) == discriminant(OK)

  OKns = maximal_order(Kns)
  @inferred nf(OKns)
  @test !issimple(OKns)
  @test iscommutative(OKns)
  @test @inferred !isequation_order(OKns)
  @test @inferred ismaximal(OKns)
  @test @inferred degree(OKns) == 2
  @test @inferred absolute_degree(OKns) == 2
  BOKns = absolute_basis(OKns)
  @test @inferred discriminant(BOKns) == @inferred discriminant(OKns)
  @test @inferred signature(OKns) == (0, 1)
  @test @inferred discriminant(OKns, FlintQQ) == discriminant(OKns)

  OL = maximal_order(L)
  @inferred nf(OL)
  @test issimple(OL)
  @test iscommutative(OL)
  @test @inferred !isequation_order(OL)
  @test @inferred ismaximal(OL)
  @test @inferred degree(OL) == 2
  @test @inferred absolute_degree(OL) == 4
  BOL = absolute_basis(OL)
  @test numerator(det(trace_matrix(map(elem_in_nf, BOL)))) == @inferred absolute_discriminant(OL)
  @test @inferred signature(OL) == (0, 2)
  @test @inferred discriminant(OL, FlintQQ) == absolute_discriminant(OL)

  OLns = maximal_order(Lns)
  @inferred nf(OLns)
  @test !issimple(OLns)
  @test iscommutative(OLns)
  @test @inferred !isequation_order(OLns)
  @test @inferred ismaximal(OLns)
  @test @inferred degree(OLns) == 2
  @test @inferred absolute_degree(OLns) == 4
  BOLns = absolute_basis(OLns)
  @test numerator(det(trace_matrix(map(elem_in_nf, BOLns)))) == @inferred absolute_discriminant(OLns)
  @test @inferred signature(OLns) == (0, 2)
  @test @inferred discriminant(OLns, FlintQQ) == absolute_discriminant(OLns)

end


text:
using Test

@test true == true


text:
#included within module Time

export easterMonday, isWeekend, isWorkingDay

abstract BusinessCalendar

abstract WesternCalendar <: BusinessCalendar
abstract OrthodoxCalendar <: BusinessCalendar

const w_easter_monday = [
                     98,  90, 103,  95, 114, 106,  91, 111, 102,   # 1901-1909
                87, 107,  99,  83, 103,  95, 115,  99,  91, 111,   # 1910-1919
                96,  87, 107,  92, 112, 103,  95, 108, 100,  91,   # 1920-1929
               111,  96,  88, 107,  92, 112, 104,  88, 108, 100,   # 1930-1939
                85, 104,  96, 116, 101,  92, 112,  97,  89, 108,   # 1940-1949
               100,  85, 105,  96, 109, 101,  93, 112,  97,  89,   # 1950-1959
               109,  93, 113, 105,  90, 109, 101,  86, 106,  97,   # 1960-1969
                89, 102,  94, 113, 105,  90, 110, 101,  86, 106,   # 1970-1979
                98, 110, 102,  94, 114,  98,  90, 110,  95,  86,   # 1980-1989
               106,  91, 111, 102,  94, 107,  99,  90, 103,  95,   # 1990-1999
               115, 106,  91, 111, 103,  87, 107,  99,  84, 103,   # 2000-2009
                95, 115, 100,  91, 111,  96,  88, 107,  92, 112,   # 2010-2019
               104,  95, 108, 100,  92, 111,  96,  88, 108,  92,   # 2020-2029
               112, 104,  89, 108, 100,  85, 105,  96, 116, 101,   # 2030-2039
                93, 112,  97,  89, 109, 100,  85, 105,  97, 109,   # 2040-2049
               101,  93, 113,  97,  89, 109,  94, 113, 105,  90,   # 2050-2059
               110, 101,  86, 106,  98,  89, 102,  94, 114, 105,   # 2060-2069
                90, 110, 102,  86, 106,  98, 111, 102,  94, 114,   # 2070-2079
                99,  90, 110,  95,  87, 106,  91, 111, 103,  94,   # 2080-2089
               107,  99,  91, 103,  95, 115, 107,  91, 111, 103,   # 2090-2099
                88, 108, 100,  85, 105,  96, 109, 101,  93, 112,   # 2100-2109
                97,  89, 109,  93, 113, 105,  90, 109, 101,  86,   # 2110-2119
               106,  97,  89, 102,  94, 113, 105,  90, 110, 101,   # 2120-2129
                86, 106,  98, 110, 102,  94, 114,  98,  90, 110,   # 2130-2139
                95,  86, 106,  91, 111, 102,  94, 107,  99,  90,   # 2140-2149
               103,  95, 115, 106,  91, 111, 103,  87, 107,  99,   # 2150-2159
                84, 103,  95, 115, 100,  91, 111,  96,  88, 107,   # 2160-2169
                92, 112, 104,  95, 108, 100,  92, 111,  96,  88,   # 2170-2179
               108,  92, 112, 104,  89, 108, 100,  85, 105,  96,   # 2180-2189
               116, 101,  93, 112,  97,  89, 109, 100,  85, 105    # 2190-2199
        ]
const o_easter_monday = [
                     105, 118, 110, 102, 121, 106, 126, 118, 102,   # 1901-1909
                122, 114,  99, 118, 110,  95, 115, 106, 126, 111,   # 1910-1919
                103, 122, 107,  99, 119, 110, 123, 115, 107, 126,   # 1920-1929
                111, 103, 123, 107,  99, 119, 104, 123, 115, 100,   # 1930-1939
                120, 111,  96, 116, 108, 127, 112, 104, 124, 115,   # 1940-1949
                100, 120, 112,  96, 116, 108, 128, 112, 104, 124,   # 1950-1959
                109, 100, 120, 105, 125, 116, 101, 121, 113, 104,   # 1960-1969
                117, 109, 101, 120, 105, 125, 117, 101, 121, 113,   # 1970-1979
                 98, 117, 109, 129, 114, 105, 125, 110, 102, 121,   # 1980-1989
                106,  98, 118, 109, 122, 114, 106, 118, 110, 102,   # 1990-1999
                122, 106, 126, 118, 103, 122, 114,  99, 119, 110,   # 2000-2009
                 95, 115, 107, 126, 111, 103, 123, 107,  99, 119,   # 2010-2019
                111, 123, 115, 107, 127, 111, 103, 123, 108,  99,   # 2020-2029
                119, 104, 124, 115, 100, 120, 112,  96, 116, 108,   # 2030-2039
                128, 112, 104, 124, 116, 100, 120, 112,  97, 116,   # 2040-2049
                108, 128, 113, 104, 124, 109, 101, 120, 105, 125,   # 2050-2059
                117, 101, 121, 113, 105, 117, 109, 101, 121, 105,   # 2060-2069
                125, 110, 102, 121, 113,  98, 118, 109, 129, 114,   # 2070-2079
                106, 125, 110, 102, 122, 106,  98, 118, 110, 122,   # 2080-2089
                114,  99, 119, 110, 102, 115, 107, 126, 118, 103,   # 2090-2099
                123, 115, 100, 120, 112,  96, 116, 108, 128, 112,   # 2100-2109
                104, 124, 109, 100, 120, 105, 125, 116, 108, 121,   # 2110-2119
                113, 104, 124, 109, 101, 120, 105, 125, 117, 101,   # 2120-2129
                121, 113,  98, 117, 109, 129, 114, 105, 125, 110,   # 2130-2139
                102, 121, 113,  98, 118, 109, 129, 114, 106, 125,   # 2140-2149
                110, 102, 122, 106, 126, 118, 103, 122, 114,  99,   # 2150-2159
                119, 110, 102, 115, 107, 126, 111, 103, 123, 114,   # 2160-2169
                 99, 119, 111, 130, 115, 107, 127, 111, 103, 123,   # 2170-2179
                108,  99, 119, 104, 124, 115, 100, 120, 112, 103,   # 2180-2189
                116, 108, 128, 119, 104, 124, 116, 100, 120, 112    # 2190-2199
        ]

easterMonday(c::WesternCalendar, y::Int) = w_easter_monday[y-1900]
easterMonday(c::OrthodoxCalendar, y::Int) = o_easter_monday[y-1900]

isWeekend(c::Union(WesternCalendar, OrthodoxCalendar), dow::Int) = (dow == 7 || dow == 1)

businessDaysBetween(c::BusinessCalendar, from::CalendarTime, to::CalendarTime) = businessDaysBetween(c, from, to, true, false)

function businessDaysBetween(c::BusinessCalendar, from::CalendarTime, to::CalendarTime, includeFirst::Bool,  includeLast::Bool)
  wd = 0
  if (from != to)
      if (from < to)
         for  d in from:days(1):to
              if (isBusinessDay(c, d))
                  wd += 1
              end
          end
      elseif (from > to)
          for d in to:days(1):from
              if (isBusinessDay(c, d))
                  wd += 1
              end
          end

      end

      if (isBusinessDay(c, from) && !includeFirst)
          wd -= 1
      end
      if (isBusinessDay(c, to) && !includeLast)
          wd -= 1
      end

      if (from > to)
          wd = -wd
      end
  end

  return wd
end

include("calendars/UnitedKingdom.jl")
include("calendars/UnitedStates.jl")
include("calendars/India.jl")


text:
##### temperature function
# @inline function tempFunc(temp, p)
#     TFunc = max(1.0e-30, exp(p.TAe * (1.0/(temp+273.15) - (1.0/p.TRef)))) * p.TCoeff
#     return TFunc
# end
@inline function tempFunc(temp, p)
    k = exp(-p.Ea/(8.3145*(temp+273.15)))*(1.0-exp(temp+273.15 - p.T⁺))
    k = max(0.0, k)
    OGT_rate = exp(-p.Ea/(8.3145*(p.T⁺-2)))
    return k/OGT_rate
end

##### calculate photosynthesis rate (mmolC/individual/second)
@inline function calc_PS(par, temp, Chl, Bm, Sz, p)
    αI  = par * p.α * p.Φ
    PCm = p.PCmax * Sz^p.PC_b * tempFunc(temp, p)
    PS  = PCm * (1.0 - exp(-αI / max(1.0e-30, PCm) * Chl / max(1.0e-30, Bm))) * Bm
    return PS
end

##### calculate nutrient uptake rate (mmolN/individual/second)
@inline function calc_NH4_uptake(NH4, temp, Cq, Nq, Bm, Sz, p)
    regQ = max(0.0, min(1.0, (p.Nqmax - (Nq + Bm * p.R_NC) / max(1.0e-30, Bm + Cq)) / (p.Nqmax - p.Nqmin)))
    VN = p.VNH4max * Sz^p.VN_b * regQ * NH4/max(1.0e-30, NH4+p.KsatNH4) * tempFunc(temp, p) * Bm
    return VN
end
@inline function calc_NO3_uptake(NO3, temp, Cq, Nq, Bm, Sz, p)
    regQ = max(0.0, min(1.0, (p.Nqmax - (Nq + Bm * p.R_NC) / max(1.0e-30, Bm + Cq)) / (p.Nqmax - p.Nqmin)))
    VN = p.VNO3max * Sz^p.VN_b * regQ * NO3/max(1.0e-30, NO3+p.KsatNO3) * tempFunc(temp, p) * Bm
    return VN
end
@inline function calc_PO4_uptake(PO4, temp, Cq, Pq, Bm, Sz, p)
    regQ = max(0.0, min(1.0, (p.Pqmax - (Pq + Bm * p.R_PC) / max(1.0e-30, Bm + Cq)) / (p.Pqmax - p.Pqmin)))
    VN = p.VPO4max * Sz^p.VP_b * regQ * PO4/max(1.0e-30, PO4+p.KsatPO4) * tempFunc(temp, p) * Bm
    return VN
end
@inline function calc_DOC_uptake(DOC, temp, Cq, Bm, Sz, p)
    regQ = max(0.0, min(1.0, (p.Cqmax - Cq / max(1.0e-30, Bm + Cq)) / (p.Cqmax - p.Cqmin)))
    VN = p.VDOCmax * Sz^p.VDOC_b * regQ * DOC/max(1.0e-30, DOC+p.KsatDOC) * tempFunc(temp, p) * Bm
    return VN
end

@kernel function calc_inorganic_uptake_kernel!(plank, nuts, p)
    i = @index(Global)
    @inbounds plank.PS[i] = calc_PS(nuts.par[i], nuts.T[i], plank.Chl[i], plank.Bm[i], plank.Sz[i], p) * plank.ac[i]

    @inbounds plank.VNH4[i] = calc_NH4_uptake(nuts.NH4[i], nuts.T[i], plank.Cq[i], plank.Nq[i], plank.Bm[i],
                                             plank.Sz[i], p) * plank.ac[i]

    @inbounds plank.VNO3[i] = calc_NO3_uptake(nuts.NO3[i], nuts.T[i], plank.Cq[i], plank.Nq[i], plank.Bm[i],
                                             plank.Sz[i], p) * plank.ac[i]

    @inbounds plank.VPO4[i] = calc_PO4_uptake(nuts.PO4[i], nuts.T[i], plank.Cq[i], plank.Pq[i], plank.Bm[i],
                                             plank.Sz[i], p) * plank.ac[i]
end
function calc_inorganic_uptake!(plank, nuts, p, arch::Architecture)
    kernel! = calc_inorganic_uptake_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, nuts, p)
    wait(device(arch), event)
    return nothing
end

##### update C, N, P quotas for the first time of each time step
@kernel function update_quotas_1_kernel!(plank, ΔT)
    i = @index(Global)
    @inbounds plank.Cq[i] += plank.PS[i]   * ΔT
    @inbounds plank.Pq[i] += plank.VPO4[i] * ΔT
    @inbounds plank.Nq[i] +=(plank.VNH4[i] + plank.VNO3[i]) * ΔT
end
function update_quotas_1!(plank, ΔT, arch)
    kernel! = update_quotas_1_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, ΔT)
    wait(device(arch), event)
    return nothing
end

##### calculate DOC uptake rate (mmolC/individual/second)
##### DOC uptake needs support of photosynthesis for at least 1% of total C acquisition.
@kernel function calc_organic_uptake_kernel!(plank, nuts, p)
    i = @index(Global)
    @inbounds plank.VDOC[i] = calc_DOC_uptake(nuts.DOC[i], nuts.T[i], plank.Cq[i], plank.Bm[i], plank.Sz[i], p) * plank.ac[i]

    @inbounds plank.VDOC[i] = plank.VDOC[i] * isless(0.01, plank.PS[i]/(plank.VDOC[i]+plank.PS[i]))
end
function calc_organic_uptake!(plank, nuts, p, arch::Architecture)
    kernel! = calc_organic_uptake_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, nuts, p)
    wait(device(arch), event)
    return nothing
end

##### calculate ρChl
@kernel function calc_ρChl_kernel!(plank, par, p)
    i = @index(Global)
    @inbounds plank.ρChl[i] = plank.PS[i]/max(1.0e-30, plank.Bm[i]) * p.Chl2N /
                             max(1.0e-30, par[i] * p.α * p.Φ * plank.Chl[i]/plank.Bm[i]) *
                             isless(1.0e-1, par[i]) * plank.ac[i]
end
function calc_ρChl!(plank, par, p, arch)
    kernel! = calc_ρChl_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, par, p)
    wait(device(arch), event)
    return nothing
end

##### calculate respiration (mmolC/individual/second)
@kernel function calc_respir_kernel!(plank, T, p)
    i = @index(Global)
    @inbounds plank.resp[i] = p.respir_a * plank.Sz[i]^p.respir_b * plank.Bm[i] * tempFunc(T[i], p) * plank.ac[i]
end
function calc_respir!(plank, T, p, arch)
    kernel! = calc_respir_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, T, p)
    wait(device(arch), event)
    return nothing
end

##### update C, N, P quotas for the second time of each time step
@kernel function update_quotas_2_kernel!(plank, ΔT, p)
    i = @index(Global)
    @inbounds plank.Cq[i] = plank.Cq[i] + (plank.VDOC[i] - plank.resp[i]) * ΔT
    @inbounds plank.Cq[i] = plank.Cq[i] + max(0.0, (0.0 - plank.Cq[i]))
    @inbounds plank.Bm[i] = plank.Bm[i] - max(0.0, (0.0 - plank.Cq[i]))
    @inbounds plank.Nq[i] = plank.Nq[i] + max(0.0, (0.0 - plank.Cq[i])) * p.R_NC
    @inbounds plank.Pq[i] = plank.Pq[i] + max(0.0, (0.0 - plank.Cq[i])) * p.R_PC
end
function update_quotas_2!(plank, ΔT, p, arch)
    kernel! = update_quotas_2_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, ΔT, p)
    wait(device(arch), event)
    return nothing
end

##### calculate biosynthesis and exudation (mmolC/individual/second)
@kernel function calc_BS_kernel!(plank, p)
    i = @index(Global)
    @inbounds plank.BS[i] = min(plank.Cq[i], plank.Nq[i]/p.R_NC, plank.Pq[i]/p.R_PC) *
                           p.k_mtb * plank.Sz[i]^p.k_mtb_b * plank.ac[i]
    @inbounds plank.exu[i]= max(0.0, plank.Cq[i] - min(plank.Cq[i], plank.Nq[i]/p.R_NC, plank.Pq[i]/p.R_PC)) *
                           p.k_mtb * plank.Sz[i]^p.k_mtb_b * plank.ac[i]
end
function calc_BS!(plank, p, arch)
    kernel! = calc_BS_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, p)
    wait(device(arch), event)
    return nothing
end

##### update C, N, P quotas, biomass, Chla
@kernel function update_biomass_kernel!(plank, p, ΔT)
    i = @index(Global)
    @inbounds plank.Bm[i]  += ΔT * plank.BS[i]
    @inbounds plank.Cq[i]  -= ΔT *(plank.BS[i] + plank.exu[i])
    @inbounds plank.Nq[i]  -= ΔT * plank.BS[i] * p.R_NC
    @inbounds plank.Pq[i]  -= ΔT * plank.BS[i] * p.R_PC
    @inbounds plank.Chl[i] += ΔT * plank.BS[i] * p.R_NC * plank.ρChl[i]
    @inbounds plank.age[i] += ΔT / 3600.0 * plank.ac[i]
end
function update_biomass!(plank, p, ΔT, arch)
    kernel! = update_biomass_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, p, ΔT)
    wait(device(arch), event)
    return nothing
end

##### update cell size
@kernel function update_cellsize_kernel!(plank, p)
    i = @index(Global)
    @inbounds plank.Sz[i] = (plank.Bm[i] + plank.Cq[i]) / (p.Cquota * p.Nsuper)
end
function update_cellsize!(plank, p, arch)
    kernel! = update_cellsize_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, p)
    wait(device(arch), event)
    return nothing
end


text:
"""
    arrows(points, directions; kwargs...)
    arrows(x, y, u, v)
    arrows(x::AbstractVector, y::AbstractVector, u::AbstractMatrix, v::AbstractMatrix)
    arrows(x, y, z, u, v, w)

Plots arrows at the specified points with the specified components.
`u` and `v` are interpreted as vector components (`u` being the x
and `v` being the y), and the vectors are plotted with the tails at
`x`, `y`.

If `x, y, u, v` are `<: AbstractVector`, then each 'row' is plotted
as a single vector.

If `u, v` are `<: AbstractMatrix`, then `x` and `y` are interpreted as
specifications for a grid, and `u, v` are plotted as arrows along the
grid.

`arrows` can also work in three dimensions.

## Attributes
$(ATTRIBUTES)
"""
@recipe(Arrows, points, directions) do scene
    attr = merge!(
        default_theme(scene),
        Attributes(
            arrowhead = automatic,
            arrowtail = automatic,
            color = :black,
            linecolor = automatic,
            arrowsize = automatic,
            linestyle = nothing,
            align = :origin,
            normalize = false,
            lengthscale = 1f0,
            colormap = theme(scene, :colormap),
            quality = 32,
            inspectable = theme(scene, :inspectable),
            markerspace = Pixel,
        )
    )
    attr[:fxaa] = automatic
    attr[:linewidth] = automatic
    # connect arrow + linecolor by default
    get!(attr, :arrowcolor, attr[:linecolor])
    attr
end

# For the matlab/matplotlib users
const quiver = arrows
const quiver! = arrows!
export quiver, quiver!

arrow_head(N, marker, quality) = marker
function arrow_head(N, marker::Automatic, quality)
    if N == 2
        return '▲'
    else
        merge([
           _circle(Point3f0(0), 0.5f0, Vec3f0(0,0,-1), quality),
           _mantle(Point3f0(0), Point3f0(0,0,1), 0.5f0, 0f0, quality)
        ])
    end
end

arrow_tail(N, marker, quality) = marker
function arrow_tail(N, marker::Automatic, quality)
    if N == 2
        nothing
    else
        merge([
            _circle(Point3f0(0,0,-1), 0.5f0, Vec3f0(0,0,-1), quality),
            _mantle(Point3f0(0,0,-1), Point3f0(0), 0.5f0, 0.5f0, quality)
        ])
    end
end

function _mantle(origin, extremity, r1, r2, N)
    dphi = 2pi / N

    # Equivalent to
    # xy = cos(atan(temp))
    # z  = sin(atan(temp))
    temp = -(r2-r1) / norm(extremity .- origin)
    xy = 1.0 / sqrt(temp^2 + 1)
    z = temp / sqrt(temp^2 + 1)

    coords = Vector{Point3f0}(undef, 2N)
    normals = Vector{Vec3f0}(undef, 2N)
    faces = Vector{GLTriangleFace}(undef, 2N)

    for (i, phi) in enumerate(0:dphi:2pi-0.5dphi)
        coords[2i - 1] = origin .+ r1 * Vec3f0(cos(phi), sin(phi), 0)
        coords[2i] = extremity .+ r2 * Vec3f0(cos(phi), sin(phi), 0)
        normals[2i - 1] = Vec3f0(xy*cos(phi), xy*sin(phi), z)
        normals[2i] = Vec3f0(xy*cos(phi), xy*sin(phi), z)
        faces[2i - 1] = GLTriangleFace(2i-1, mod1(2i+1, 2N), 2i)
        faces[2i] = GLTriangleFace(mod1(2i+1, 2N), mod1(2i+2, 2N), 2i)
    end

    GeometryBasics.Mesh(meta(coords; normals=normals), faces)
end

# GeometryBasics.Circle doesn't work with Point3f0...
function _circle(origin, r, normal, N)
    dphi = 2pi / N

    coords = Vector{Point3f0}(undef, N+1)
    normals = fill(normal, N+1)
    faces = Vector{GLTriangleFace}(undef, N)

    for (i, phi) in enumerate(0:dphi:2pi-0.5dphi)
        coords[i] = origin .+ r * Vec3f0(cos(phi), sin(phi), 0)
        faces[i] = GLTriangleFace(N+1, mod1(i+1, N), i)
    end
    coords[N+1] = origin

    GeometryBasics.Mesh(meta(coords; normals=normals), faces)
end

convert_arguments(::Type{<: Arrows}, x, y, u, v) = (Point2f0.(x, y), Vec2f0.(u, v))
function convert_arguments(::Type{<: Arrows}, x::AbstractVector, y::AbstractVector, u::AbstractMatrix, v::AbstractMatrix)
    (vec(Point2f0.(x, y')), vec(Vec2f0.(u, v)))
end
convert_arguments(::Type{<: Arrows}, x, y, z, u, v, w) = (Point3f0.(x, y, z), Vec3f0.(u, v, w))

function plot!(arrowplot::Arrows{<: Tuple{AbstractVector{<: Point{N, T}}, V}}) where {N, T, V}
    @extract arrowplot (
        points, directions, colormap, normalize, align,
        arrowtail, color, linecolor, linestyle, linewidth, lengthscale,
        arrowhead, arrowsize, arrowcolor, quality,
        # passthrough
        lightposition, ambient, diffuse, specular, shininess,
        fxaa, ssao, transparency, visible, inspectable
    )

    arrow_c = map((a, c)-> a === automatic ? c : a , arrowcolor, color)
    line_c = map((a, c)-> a === automatic ? c : a , linecolor, color)

    if N == 2
        fxaa_bool = @lift($fxaa == automatic ? false : $fxaa)
        headstart = lift(points, directions, normalize, align, lengthscale) do points, dirs, n, align, s
            map(points, dirs) do p1, dir
                dir = n ? StaticArrays.normalize(dir) : dir
                if align in (:head, :lineend, :tailend, :headstart, :center)
                    shift = s .* dir
                else
                    shift = Vec2f0(0)
                end
                Point2f0(p1 .- shift) => Point2f0(p1 .- shift .+ (dir .* s))
            end
        end

        scene = parent_scene(arrowplot)
        rotations = directions

        # for 2D arrows, compute the correct marker rotation given the projection / scene size
        # for the screen-space marker
        if arrowplot.markerspace[] == Pixel
            rotations = lift(scene.camera.projectionview, scene.px_area, headstart) do pv, pxa, hs
                angles = map(hs) do (start, stop)
                    pstart = project(scene, start)
                    pstop = project(scene, stop)
                    diff = pstop - pstart
                    n = norm(diff)
                    if n == 0
                        zero(n)
                    else
                        angle = acos(diff[2] / norm(diff))
                        angle = ifelse(diff[1] > 0, 2pi - angle, angle)
                    end
                end
                Billboard(angles)
            end
        end

        linesegments!(
            arrowplot, headstart,
            color = line_c, colormap = colormap, linestyle = linestyle,
            linewidth = @lift($linewidth === automatic ? 1f0 : $linewidth),
            fxaa = fxaa_bool, inspectable = inspectable,
            transparency = transparency, visible = visible,
        )
        scatter!(
            arrowplot,
            lift(x-> last.(x), headstart),
            marker = @lift(arrow_head(2, $arrowhead, $quality)),
            markersize = @lift($arrowsize === automatic ? theme(scene, :markersize)[] : $arrowsize),
            color = arrow_c, rotations = rotations, strokewidth = 0.0,
            colormap = colormap, markerspace = arrowplot.markerspace,
            fxaa = fxaa_bool, inspectable = inspectable,
            transparency = transparency, visible = visible
        )
    else
        fxaa_bool = @lift($fxaa == automatic ? true : $fxaa)
        start = lift(points, directions, align, lengthscale) do points, dirs, align, s
            map(points, dirs) do p, dir
                if align in (:head, :lineend, :tailend, :headstart, :center)
                    shift = Vec3f0(0)
                else
                    shift = -s .* dir
                end
                Point3f0(p .- shift)
            end
        end
        meshscatter!(
            arrowplot,
            start, rotations = directions,
            marker = @lift(arrow_tail(3, $arrowtail, $quality)),
            markersize = lift(directions, normalize, linewidth, lengthscale) do dirs, n, linewidth, ls
                lw = linewidth === automatic ? 0.05f0 : linewidth
                if n
                    Vec3f0(lw, lw, ls)
                else
                    map(dir -> Vec3f0(lw, lw, norm(dir) * ls), dirs)
                end
            end,
            color = line_c, colormap = colormap,
            fxaa = fxaa_bool, ssao = ssao,
            lightposition = lightposition, ambient = ambient, diffuse = diffuse,
            specular = specular, shininess = shininess, inspectable = inspectable,
            transparency = transparency, visible = visible
        )
        meshscatter!(
            arrowplot,
            start, rotations = directions,
            marker = @lift(arrow_head(3, $arrowhead, $quality)),
            markersize = lift(Any, arrowsize) do as
                as === automatic ? Vec3f0(0.2, 0.2, 0.3) : as
            end,
            color = arrow_c, colormap = colormap,
            fxaa = fxaa_bool, ssao = ssao,
            lightposition = lightposition, ambient = ambient, diffuse = diffuse,
            specular = specular, shininess = shininess, inspectable = inspectable,
            transparency = transparency, visible = visible
        )
    end

end


text:
module TestConstraints

using Test
import MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

function test_normalize_and_add_constrant_VariableIndex()
    model = MOI.Utilities.Model{Float64}()
    x = MOI.add_variable(model)
    ci = MOI.Utilities.normalize_and_add_constraint(
        model,
        x,
        MOI.EqualTo(1.0),
        allow_modify_function = false,
    )
    @test MOI.get(model, MOI.ConstraintFunction(), ci) == x
    @test MOI.get(model, MOI.ConstraintSet(), ci) == MOI.EqualTo(1.0)
    return
end

function test_normalize_and_add_constrant_ScalarAffineFunction()
    model = MOI.Utilities.Model{Float64}()
    x = MOI.add_variable(model)
    f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 2.0)
    g = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0)
    ci = MOI.Utilities.normalize_and_add_constraint(model, f, MOI.EqualTo(3.0))
    @test f.constant == 2.0
    @test MOI.get(model, MOI.ConstraintFunction(), ci) ≈ g
    @test MOI.get(model, MOI.ConstraintSet(), ci) == MOI.EqualTo(1.0)
    ci = MOI.Utilities.normalize_and_add_constraint(
        model,
        f,
        MOI.Interval(-1.0, 1.0),
        allow_modify_function = true,
    )
    @test f.constant == 0.0
    @test MOI.get(model, MOI.ConstraintFunction(), ci) ≈ g
    @test MOI.get(model, MOI.ConstraintSet(), ci) == MOI.Interval(-3.0, -1.0)
    return
end

end  # module

TestConstraints.runtests()