text:
using Test

struct Entry
    Ingredients::Set{String}
    Allergens::Set{String}
end

function find_safe_ingredients(a, b)
    out = Set()
    # If they *dont* have overlapping allergens
    if length(intersect(a.Allergens, b.Allergens)) > 0
        return out
    end
    # the overlapping ingredients cannot be the cause of allergens
    return intersect(a.Ingredients, b.Ingredients)
end

# TODO
# INSTEAD, try this algo
# 1. find everything that has the allergen
# 2. set if they have exactly one common ingredient

function run(fname)
    f = open(fname, "r")
    lines = readlines(f)
    close(f)

    entries = []
    for l in lines
        l = l[1:end - 1]
        parts = split(l, " (contains ")
        ingredients = split(parts[1])
        allergens = split(parts[2], ", ")
        push!(entries, Entry(Set(ingredients), Set(allergens)))
    end
    entries_bak = copy(entries)

    for e in entries
        @show e
    end
    all_allergens = reduce(union, map((e) -> e.Allergens, entries))
    @show all_allergens
    ingredient_to_allergen = Dict()
    while length(all_allergens) > 0
        for a in all_allergens
            @show a
            entries_w_allergen = filter((e) -> in(a, e.Allergens), entries)
            @show entries_w_allergen
            ingredients_that_may_cause = intersect(map((e) -> e.Ingredients, entries_w_allergen)...)
            @show ingredients_that_may_cause
            if length(ingredients_that_may_cause) == 1
                ingr = first(ingredients_that_may_cause)
                # we found it!
                ingredient_to_allergen[ingr] = a
                delete!(all_allergens, a)
                # remove it from entries as well
                for e in entries
                    delete!(e.Ingredients, ingr)
                    delete!(e.Allergens, a)
                end
                # break, since we modified all_allergens and are looping within it
                break
            end
        end
    end

    # Part 2: Sort alphabetically *by their allergen*
    dangerous_ingredients = map(x -> x[1], sort(collect(ingredient_to_allergen), by=x -> x[2]))
    result = 0
    for e in entries_bak
        result += length(setdiff(e.Ingredients, dangerous_ingredients))
    end

    # remove inert ingredients
    result2 = join(dangerous_ingredients, ",")
    println("[part 1] File = ", fname, " got result = ", result)
    println("[part 2] File = ", fname, " got result = ", result2)
    return result, result2
end

# mxmxvkd kfcds sqjhc nhms (contains dairy, fish)
# trh fvjkl sbzzf mxmxvkd (contains dairy)
# sqjhc fvjkl (contains soy)
# sqjhc mxmxvkd sbzzf (contains fish)

# manually, seems like
# mxmxvkd => dairy
# sqjhc => fish
# fvjkl => soy

# Each allergen is found in exactly one ingredient.
# Each ingredient contains zero or one allergen.
# **Allergens aren't always marked**, but if they are marked, the corresponding ingredien t is definitely present
@test run("21ex.txt") == (5, "mxmxvkd,sqjhc,fvjkl")

run("21.txt")


text:
function AUSM_plus(Qbase, Qcon, cellxmax, cellymax, vecAx, vecAy, specific_heat_ratio, volume, nval)
    E_adv_hat = zeros(cellxmax+1,   cellymax, nval)
    F_adv_hat = zeros(  cellxmax, cellymax+1, nval)
    g         = specific_heat_ratio

    for i in 2:cellxmax+1 -1
        for j in 2:cellymax -1
            # i-1セル
            cell_vecAx = zeros(2)
            cell_vecAx[1] = 0.5*(vecAx[i-1,j,1]+vecAx[i,j,1])
            cell_vecAx[2] = 0.5*(vecAx[i-1,j,2]+vecAx[i,j,2])

            rhoL = Qbase[i-1,j,1]
            UL = Qcon[i-1,j,2] / Qcon[i-1,j,1]*cell_vecAx[1] + Qcon[i-1,j,3] / Qcon[i-1,j,1]*cell_vecAx[2]
            pL = Qbase[i-1,j,4]

            # iセル
            cell_vecAx = zeros(2)
            cell_vecAx[1] = 0.5*(vecAx[i,j,1]+vecAx[i+1,j,1])
            cell_vecAx[2] = 0.5*(vecAx[i,j,2]+vecAx[i+1,j,2])

            rhoR = Qbase[i,j,1]
            #UR = Qbase[i,j,2]*cell_vecAx[1] + Qbase[i,j,3]*cell_vecAx[2]
            UR = Qcon[i,j,2] / Qcon[i,j,1]*cell_vecAx[1] + Qcon[i,j,3] / Qcon[i,j,1]*cell_vecAx[2]
            pR = Qbase[i,j,4]

            # AUSM+
            mdot, ph = AUSM_plus_half(rhoL, rhoR, UL, UR, pL, pR, g)

            # flux half
            temp_vecX = zeros(nval)
            temp_vecX[2] = vecAx[i,j,1]
            temp_vecX[3] = vecAx[i,j,2]

            Lpsi = zeros(nval)
            Rpsi = zeros(nval)

            for l in 1:nval
                Lpsi[l] = Qcon[i-1,j,l] / Qcon[i-1,j,1]
                Rpsi[l] = Qcon[i,j,l]   / Qcon[i,j,1]
            end
            Lpsi[4] = (Qcon[i-1,j,4] + Qbase[i-1,j,4]) / Qcon[i-1,j,1]
            Rpsi[4] = (Qcon[i,j,4]   + Qbase[i,j,4])   / Qcon[i,j,1]

            if mdot > 0
                for l in 1:nval
                    E_adv_hat[i,j,l] = mdot * Lpsi[l] + ph * temp_vecX[l]
                end
            else
                for l in 1:nval
                    E_adv_hat[i,j,l] = mdot * Rpsi[l] + ph * temp_vecX[l]
                end
            end
        end
    end

    for i in 2:cellxmax -1
        for j in 2:cellymax+1 -1

            # j-1セル
            cell_vecAy = zeros(2)
            cell_vecAy[1] = 0.5*(vecAy[i,j-1,1]+vecAy[i,j,1])
            cell_vecAy[2] = 0.5*(vecAy[i,j-1,2]+vecAy[i,j,2])

            rhoL = Qbase[i,j-1,1]
            VL = Qbase[i,j-1,2]*cell_vecAy[1] + Qbase[i,j-1,3]*cell_vecAy[2]
            pL = Qbase[i,j-1,4]

            # jセル
            cell_vecAy = zeros(2)
            cell_vecAy[1] = 0.5*(vecAy[i,j,1]+vecAy[i,j+1,1])
            cell_vecAy[2] = 0.5*(vecAy[i,j,2]+vecAy[i,j+1,2])

            rhoR = Qbase[i,j,1]
            VR = Qbase[i,j,2]*cell_vecAy[1] + Qbase[i,j,3]*cell_vecAy[2]
            pR = Qbase[i,j,4]

            # AUSM+up
            mdot, ph = AUSM_plus_half(rhoL, rhoR, VL, VR, pL, pR, g)

            # flux half
            temp_vecY = zeros(nval)
            temp_vecY[2] = vecAy[i,j,1]
            temp_vecY[3] = vecAy[i,j,2]

            Lpsi = zeros(nval)
            Rpsi = zeros(nval)

            for l in 1:nval
                Lpsi[l] = Qcon[i,j-1,l] / Qcon[i,j-1,1]
                Rpsi[l] = Qcon[i,j,l]   / Qcon[i,j,1]
            end
            Lpsi[4] = (Qcon[i,j-1,4] + Qbase[i,j-1,4]) / Qcon[i,j-1,1]
            Rpsi[4] = (Qcon[i,j,4]   + Qbase[i,j,4])   / Qcon[i,j,1]

            if mdot > 0
                for l in 1:nval
                    F_adv_hat[i,j,l] = mdot * Lpsi[l] + ph * temp_vecY[l]
                end
            else
                for l in 1:nval
                    F_adv_hat[i,j,l] = mdot * Rpsi[l] + ph * temp_vecY[l]
                end
            end
        end
    end

    return E_adv_hat, F_adv_hat
end

function AUSM_plus_half(rhoL, rhoR, UL, UR, pL, pR, g)
    # param
    beta  = 1/8
    alpha = 3/16

    # L, R
    aL = (g * pL / rhoL)^0.5
    aR = (g * pR / rhoR)^0.5

    # half
    ah   = 0.5*( aL + aR )
    rhoh = 0.5*( rhoL + rhoR )

    ML = UL/ah
    MR = UR/ah

    M_p4 = 0
    M_m4 = 0
    p_p5 = 0
    p_m5 = 0
    if abs(ML) >= 1
        M_p4 = 0.5*(ML + abs(ML))
        p_p5 = 0.5*(ML + abs(ML)) / ML
    else
        M_p4 = 0.5*(ML+1)^2 + beta*(ML^2-1)^2
        p_p5 = 0.25*(ML+1)^2 * (2-ML) + alpha*ML*(ML^2-1)^2
    end

    if abs(MR) >= 1
        M_m4 = 0.5*(MR - abs(MR))
        p_m5 = 0.5*(MR - abs(MR)) / MR
    else
        M_m4 = -0.5*(MR-1)^2 - beta*(MR^2-1)^2
        p_m5 = 0.25*(MR-1)^2 * (2+MR) - alpha*MR*(MR^2-1)^2
    end

    # Mh = M_p4 + M_m4 + Mp/fa
    Mh = M_p4 + M_m4

    # mdot half
    mdot = ah * Mh
    if Mh > 0
        mdot = mdot * rhoL
    else
        mdot = mdot * rhoR
    end

    # p half
    ph = p_p5*pL + p_m5*pR

    return mdot, ph
end

ext: jl

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
"""
Sub-module of the module `DynamicalSystems`, which contains pre-defined
famous systems.
"""
module Systems
using DynamicalSystems, StaticArrays
#######################################################################################
#                                    Continuous                                       #
#######################################################################################
"""
    lorenz(u0=[0.0, 10.0, 0.0]; σ = 10.0, ρ = 28.0, β = 8/3)
The famous three dimensional system due to Lorenz [1], shown to exhibit
so-called "deterministic nonperiodic flow". It was originally invented to study a
simplified form of atmospheric convection.

Currently, it is most famous for its strange attractor (occuring at the default
parameters), which resembles a butterfly. For the same reason it is
also associated with the term "butterfly effect" (a term which Lorenz himself disliked)
even though the effect applies generally to dynamical systems.

Default values are the ones used in the original paper.

[1] E. N. Lorenz, J. atmos. Sci. **20**, pp 130 (1963)
"""
function lorenz(u0=[0.0, 10.0, 0.0]; σ = 10.0, ρ = 28.0, β = 8/3)
  @inline @inbounds eom_lorenz(u) =
  SVector{3}(σ*(u[2]-u[1]), u[1]*(ρ-u[3]) - u[2], u[1]*u[2] - β*u[3])
  @inline @inbounds function jacob_lorenz(u)
    i = one(eltype(u))
    o = zero(eltype(u))
    @SMatrix [-σ*i           σ*i    zero(i);
              (ρ*i - u[3])   (-i)   (-u[1]);
              u[2]           u[1]   (-β*i) ]
  end# should give exponents [0.9056, 0, -14.5723]
  return ContinuousDS(u0, eom_lorenz, jacob_lorenz)
end

"""
```julia
roessler(u0=rand(3); a = 0.2, b = 0.2, c = 5.7)
```
This three-dimensional continuous system is due to Rössler [1].
It is a system that by-design behaves similarly
to the `lorenz` system and displays a (fractal)
strange attractor. However, it is easier to analyze qualitatively, as for example
the attractor is composed of a single manifold.

Default values are the same as the original paper.

[1] O. E. Rössler, Phys. Lett. **57A**, pp 397 (1976)
"""
function roessler(u0=rand(3); a = 0.2, b = 0.2, c = 5.7)
  @inline @inbounds eom_roessler(u) =
  SVector{3}(-u[2]-u[3], u[1] + a*u[2], b + u[3]*(u[1] - c))
  @inline @inbounds function jacob_roessler(u)
    i = one(eltype(u))
    o = zero(eltype(u))
    @SMatrix [o     -i      -i;
              i      a       o;
              u[3]   o       u[1] - c]
  end
  return ContinuousDS(u0, eom_roessler, jacob_roessler)
end

#######################################################################################
#                                     Discrete                                        #
#######################################################################################
"""
    towel(u0 = [0.085, -0.121, 0.075])

The folded-towel map is a hyperchaotic mapping due to Rössler [1]. It is famous
for being a mapping that has the smallest possible dimensions necessary for hyperchaos,
having two positive and one negative lyapunov exponent.

The name comes from the fact that when plotted looks like a folded towel, in every
projection.

Default values are the ones used in the original paper.

[1] : O. E. Rössler, Phys. Lett. **71A**, pp 155 (1979)
"""
function towel(u0=[0.085, -0.121, 0.075])
  @inline @inbounds function eom_towel(x)
    x1, x2, x3 = x[1], x[2], x[3]
    SVector( 3.8*x1*(1-x1) - 0.05*(x2+0.35)*(1-2*x3),
    0.1*( (x2+0.35)*(1-2*x3) - 1 )*(1 - 1.9*x1),
    3.78*x3*(1-x3)+0.2*x2 )
  end

  @inline @inbounds function jacob_towel(x)
    @SMatrix [3.8*(1 - 2x[1]) -0.05*(1-2x[3]) 0.1*(x[2] + 0.35);
    -0.19((x[2] + 0.35)*(1-2x[3]) - 1)  0.1*(1-2x[3])*(1-1.9x[1])  -0.2*(x[2] + 0.35)*(1-1.9x[1]);
    0.0  0.2  3.78(1-2x[3]) ]
  end
  return DiscreteDS(u0, eom_towel, jacob_towel)
end# should result in lyapunovs: [0.432207,0.378834,-3.74638]

"""
    standardmap(u0=0.001rand(2); k = 0.971635)

The standard map (also known as Chirikov standard map) is a two dimensional,
area-preserving chaotic mapping due to Chirikov [1]. It is one of the most studied
chaotic systems and by far the most studied Hamiltonian (area-preserving) mapping.

The map corresponds to the  Poincaré's surface of section of the kicked rotor system.
Changing the non-linearity parameter `k` transitions the system from completely periodic
motion, to quasi-periodic, to local chaos (mixed phase-space) and finally to global
chaos.

The default parameter `k` is the critical parameter where the golden-ratio torous is
destroyed, as was calculated by Greene [2]. The e.o.m. considers the angle variable
`θ` to be the first, and the angular momentum `p` to be the second variable, while
both variables
are always taken modulo 2π (the mapping is on the [0,2π)² torus).

[1] : B. V. Chirikov, Preprint N. **267**, Institute of Nuclear Physics, Novosibirsk (1969)

[2] : J. M. Greene, J. Math. Phys. **20**, pp 1183 (1979)
"""
function standardmap(u0=0.001rand(2); k = 0.971635)
  @inline @inbounds eom_standard(x) =
  SVector{2}(mod2pi(x[1] + x[2] + k*sin(x[1])), mod2pi(x[2] + k*sin(x[1])))
  @inline @inbounds jacob_standard(x) =
  @SMatrix [1 + k*cos(x[1])    1;
            k*cos(x[1])        1]
  return DiscreteDS(u0, eom_standard, jacob_standard)
end

"""
    henon(u0=zeros(2); a = 1.4, b = 0.3)
The Hénon map is a two-dimensional mapping due to Hénon [1] that can display a strange
attractor (at the default parameters). In addition, it also displays many other aspects
of chaos, like period doubling or intermittency, for other parameters.

According to the author, it is a system displaying all the properties of the
Lorentz system (1963) while being as simple as possible.

Default values are the ones used in the original paper.

[1] : M. Hénon, Commun.Math. Phys. **50**, pp 69 (1976)
"""
function henon(u0=zeros(2); a = 1.4, b = 0.3)
  @inline @inbounds eom_henon(x) = SVector{2}(1.0 - a*x[1]^2 + x[2], b*x[1])
  @inline @inbounds jacob_henon(x) = @SMatrix [-2*a*x[1] 1.0; b 0.0]
  # should give exponents 0.4189 -1.6229
  return DiscreteDS(u0, eom_henon, jacob_henon)
end

"""
    logistic(x0 = rand(); r = 4.0)

The logistic map is an one dimensional unimodal mapping due to May [1] and is used by
many as the archetypal example of how chaos can arise from very simple equations.

Originally intentend to be a discretized model of polulation dynamics, it is now famous
for its bifurcation diagram, an immensly complex graph that that was shown
be universal by Feigenbaum [2].

[1] : R. M. May, Nature **261**, pp 459 (1976)

[2] : M. J. Feigenbaum, J. Stat. Phys. **19**, pp 25 (1978)
"""
function logistic(x0=rand(); r = 4.0)
  @inline eom_logistic(x) = r*x*(1-x)
  @inline deriv_logistic(x) = r*(1-2x)
  return DiscreteDS1D(x0, eom_logistic, deriv_logistic)
end

end# Systems module


text:
import Pkg

BACKEND = get(ENV, "KERNELABSTRACTIONS_TEST_BACKEND", "all")
CI = parse(Bool, get(ENV, "CI", "false"))

pkgs = [
    "KernelAbstractions",
]
if !CI || BACKEND == "ROCM"
    push!(pkgs, "ROCKernels")
end
if !CI || BACKEND == "CUDA"
    push!(pkgs, "CUDAKernels")
end
if !CI || haskey(ENV, "TEST_KERNELGRADIENTS")
    if VERSION < v"1.8"
        push!(pkgs, "KernelGradients")
    end
end

Pkg.test(pkgs; coverage = true)


text:
# the fundamental struct for plotting
struct PGData
    data::Dict{Symbol, DataFrames.DataFrame}
    time::Union{StepRange{Dates.DateTime}, Vector{Dates.DateTime}}
end

function PGData(
    data::Dict{PSI.OptimizationContainerKey, DataFrames.DataFrame},
    time::Union{StepRange{Dates.DateTime}, Vector{Dates.DateTime}},
)
    d = Dict(zip(Symbol.(PSI.encode_keys_as_strings(keys(data))), values(data)))

    rename_load!(d)
    return PGData(d, time)
end

function PGData(
    data::Dict{String, DataFrames.DataFrame},
    time::Union{StepRange{Dates.DateTime}, Vector{Dates.DateTime}},
)
    d = Dict(zip(Symbol.(keys(data))), values(data))

    rename_load!(d)
    return PGData(d, time)
end

function PGData(data::Dict{String, DataFrames.DataFrame})
    d = Dict(zip(Symbol.(keys(data)), no_datetime.(values(data))))
    return PGData(d, first(values(data)).DateTime)
end

# Rename Load variables: TODO: find a better way to do this
function rename_load!(load_values::Dict)
    for (k, v) in load_values
        if haskey(LOAD_RENAMING, k)
            @debug "renaming" k => LOAD_RENAMING[k]
            load_values[LOAD_RENAMING[k]] = v
            pop!(load_values, k)
        end
    end
end

#### Generation Names ####
function get_generation_variable_keys(
    results::IS.Results;
    variable_keys::Vector{T} = PSI.list_variable_keys(results),
) where {T <: PSI.OptimizationContainerKey}
    # TODO: add slacks
    filter_keys = Vector{PSI.OptimizationContainerKey}()
    for k in variable_keys
        if (
            PSI.get_component_type(k) <: PSY.Generator &&
            PSI.get_entry_type(k) == PSI.ActivePowerVariable
        ) || PSI.get_entry_type(k) ∈ keys(BALANCE_SLACKVARS)
            push!(filter_keys, k)
        end
    end

    return filter_keys
end

function get_storage_variable_keys(
    results::IS.Results;
    variable_keys::Vector{T} = PSI.list_variable_keys(results),
) where {T <: PSI.OptimizationContainerKey}
    filter_keys = Vector{PSI.OptimizationContainerKey}()
    for k in variable_keys
        if PSI.get_component_type(k) <: PSY.Storage &&
           PSI.get_entry_type(k) ∈ SUPPORTED_STORAGE_VARIABLES
            push!(filter_keys, k)
        end
    end
    return filter_keys
end

function get_generation_parameter_keys(
    results::IS.Results;
    parameter_keys::Vector{T} = PSI.list_parameter_keys(results),
) where {T <: PSI.OptimizationContainerKey}
    filter_keys = Vector{PSI.OptimizationContainerKey}()
    for k in parameter_keys
        if PSI.get_component_type(k) <: PSY.Generator &&
           PSI.get_entry_type(k) == PSI.ActivePowerTimeSeriesParameter
            push!(filter_keys, k)
        end
    end

    return filter_keys
end

#### Load Names ####
function get_load_variable_keys(
    results::IS.Results;
    variable_keys::Vector{T} = PSI.list_variable_keys(results),
) where {T <: PSI.OptimizationContainerKey}
    filter_keys = Vector{PSI.OptimizationContainerKey}()
    for k in variable_keys
        if PSI.get_component_type(k) <: PSY.ElectricLoad &&
           PSI.get_entry_type(k) ∈ SUPPORTED_LOAD_VARIABLES
            push!(filter_keys, k)
        end
    end
    return filter_keys
end

function get_load_parameter_keys(
    results::IS.Results;
    parameter_keys::Vector{T} = PSI.list_parameter_keys(results),
) where {T <: PSI.OptimizationContainerKey}
    filter_keys = Vector{PSI.OptimizationContainerKey}()
    for k in parameter_keys
        if PSI.get_component_type(k) <: PSY.ElectricLoad &&
           PSI.get_entry_type(k) == PSI.ActivePowerTimeSeriesParameter
            push!(filter_keys, k)
        end
    end

    return filter_keys
end

#### Service Names ####
function get_service_variable_keys(
    results::IS.Results;
    variable_keys::Vector{T} = PSI.list_variable_keys(results),
) where {T <: PSI.OptimizationContainerKey}
    filter_keys = Vector{PSI.OptimizationContainerKey}()
    for k in variable_keys
        if PSI.get_component_type(k) <: PSY.Service &&
           PSI.get_entry_type(k) ∈ SUPPORTED_SERVICE_VARIABLES
            push!(filter_keys, k)
        end
    end
    return filter_keys
end

function get_service_parameter_names(
    results::IS.Results;
    parameter_keys::Vector{T} = PSI.list_parameter_keys(results),
) where {T <: PSI.OptimizationContainerKey}
    filter_keys = Vector{PSI.OptimizationContainerKey}()
    for k in parameter_keys
        if PSI.get_component_type(k) <: PSY.ElectricLoad &&
           PSI.get_entry_type(k) == PSI.RequirementTimeSeriesParameter
            push!(filter_keys, k)
        end
    end

    return filter_keys
end

no_datetime(df::DataFrames.DataFrame) = df[:, propertynames(df) .!== :DateTime]

function add_fixed_parameters!(
    variables::Dict{V, DataFrames.DataFrame},
    parameters::Dict{P, DataFrames.DataFrame},
) where {V <: PSI.OptimizationContainerKey, P <: PSI.OptimizationContainerKey}
    # fixed output should be added to plots when there exists a parameter of the form
    # :P__max_active_power__* but there is no corresponding :P__* variable
    for (param_key, param) in parameters
        PSI.get_component_type(param_key) ∈ PSI.get_component_type.(keys(variables)) &&
            continue
        if !haskey(variables, param_key)
            mult = PSI.get_component_type(param_key) ∈ NEGATIVE_PARAMETERS ? -1.0 : 1.0
            variables[param_key] = param
            variables[param_key][:, propertynames(param) .!== :DateTime] .*= mult
        end
    end
end

# finds the parameters corresponding to variables for curtailment calculations
function _curtailment_parameters(
    parameter_keys::Vector{PSI.OptimizationContainerKey},
    variable_keys::Vector{PSI.OptimizationContainerKey},
)
    curtailable_parameters = parameter_keys[findall(
        in(SUPPORTED_CURTAILMENT_PARAMETERS),
        PSI.get_entry_type.(parameter_keys),
    )]
    curtailable_variables = variable_keys[findall(
        in(SUPPORTED_CURTAILMENT_VARIABLES),
        PSI.get_entry_type.(variable_keys),
    )]

    curtailment_parameters = Vector{
        NamedTuple{
            (:parameter, :variable),
            Tuple{PSI.OptimizationContainerKey, PSI.OptimizationContainerKey},
        },
    }()
    for pk in curtailable_parameters
        for cv in curtailable_variables[PSI.get_component_type.(
            curtailable_variables,
        ) .== PSI.get_component_type(pk)]
            push!(curtailment_parameters, (parameter = pk, variable = cv))
        end
    end
    return unique(curtailment_parameters)
end

function _filter_curtailment!(
    variable_values::Dict,
    parameter_values::Dict,
    curtailment_parameters::Vector{
        NamedTuple{
            (:parameter, :variable),
            Tuple{PSI.OptimizationContainerKey, PSI.OptimizationContainerKey},
        },
    },
)
    for curtailment in curtailment_parameters
        curtailment_var_key = PSI.VariableKey(
            PSI.get_entry_type(curtailment.variable),
            PSI.get_component_type(curtailment.variable),
            "Curtailment",
        )

        curt =
            parameter_values[curtailment.parameter] .- variable_values[curtailment.variable]
        if haskey(variable_values, curtailment_var_key)
            variable_values[curtailment_var_key] =
                hcat(variable_values[curtailment_var_key], no_datetime(curt))
        else
            variable_values[curtailment_var_key] = curt
        end
    end
end

function get_generation_data(results::R; kwargs...) where {R <: IS.Results}
    initial_time = get(kwargs, :initial_time, get(kwargs, :start_time, nothing))
    len = get(kwargs, :horizon, get(kwargs, :len, nothing))
    variable_keys = get(kwargs, :variable_keys, PSI.list_variable_keys(results))
    parameter_keys = get(kwargs, :parameter_keys, PSI.list_parameter_keys(results))
    curtailment = get(kwargs, :curtailment, true)
    storage = get(kwargs, :storage, true)

    if curtailment && (haskey(kwargs, :variable_keys) || haskey(kwargs, :parameter_keys))
        @warn "Cannot guarantee curtailment calculations with specified keys"
    end

    injection_keys = get_generation_variable_keys(results; variable_keys = variable_keys)
    if storage
        injection_keys = vcat(
            injection_keys,
            get_storage_variable_keys(results; variable_keys = variable_keys),
        )
    end

    parameter_keys = get_generation_parameter_keys(results; parameter_keys = parameter_keys)

    variables = PSI.read_variables_with_keys(
        results,
        injection_keys;
        start_time = initial_time,
        len = len,
    )
    parameters = PSI.read_parameters_with_keys(
        results,
        parameter_keys;
        start_time = initial_time,
        len = len,
    )

    add_fixed_parameters!(variables, parameters)

    if curtailment
        curtailment_parameters = _curtailment_parameters(parameter_keys, injection_keys)
        _filter_curtailment!(variables, parameters, curtailment_parameters)
    end

    timestamps = PSI.get_realized_timestamps(results; start_time = initial_time, len = len)
    return PGData(variables, timestamps)
end

function get_load_data(results::R; kwargs...) where {R <: IS.Results}
    initial_time = get(kwargs, :initial_time, get(kwargs, :start_time, nothing))
    len = get(kwargs, :horizon, get(kwargs, :len, nothing))
    variable_keys = get(kwargs, :variable_keys, PSI.list_variable_keys(results))
    parameter_keys = get(kwargs, :parameter_keys, PSI.list_parameter_keys(results))

    variable_keys = get_load_variable_keys(results; variable_keys = variable_keys)
    parameter_keys = get_load_parameter_keys(results; parameter_keys = parameter_keys)

    variables = PSI.read_variables_with_keys(
        results,
        variable_keys;
        start_time = initial_time,
        len = len,
    )
    parameters = PSI.read_parameters_with_keys(
        results,
        parameter_keys;
        start_time = initial_time,
        len = len,
    )

    add_fixed_parameters!(variables, parameters)

    timestamps = PSI.get_realized_timestamps(results; start_time = initial_time, len = len)
    return PGData(variables, timestamps)
end

################################### INPUT DEMAND #################################

function _get_loads(system::PSY.System, bus::PSY.Bus)
    return [
        load for load in PSY.get_components(PSY.PowerLoad, system, PSY.get_available) if
        PSY.get_bus(load) == bus
    ]
end
function _get_loads(system::PSY.System, agg::T) where {T <: PSY.AggregationTopology}
    return PSY.get_components_in_aggregation_topology(PSY.PowerLoad, system, agg)
end
function _get_loads(system::PSY.System, load::PSY.PowerLoad)
    return [load]
end
function _get_loads(system::PSY.System, sys::PSY.System)
    return PSY.get_components(PSY.PowerLoad, system, PSY.get_available)
end

get_base_power(system::PSY.System) = PSY.get_base_power(system)
get_base_power(results::PSI.SimulationProblemResults) = IS.get_base_power(results)
get_base_power(results::PSI.ProblemResults) = results.base_power

function get_load_data(
    system::PSY.System;
    aggregation::Union{
        Type{PSY.PowerLoad},
        Type{PSY.Bus},
        Type{PSY.System},
        Type{<:PSY.AggregationTopology},
    } = PSY.PowerLoad,
    kwargs...,
)
    aggregation_components =
        aggregation == PSY.System ? [system] : PSY.get_components(aggregation, system)
    if isempty(aggregation_components)
        throw(ArgumentError("System does not have type $aggregation."))
    end
    horizon = get(kwargs, :horizon, PSY.get_forecast_horizon(system))
    initial_time = get(kwargs, :initial_time, PSY.get_forecast_initial_timestamp(system))
    parameters = Dict{Symbol, DataFrames.DataFrame}()
    PSY.set_units_base_system!(system, "SYSTEM_BASE")
    for agg in aggregation_components
        loads = _get_loads(system, agg)
        length(loads) == 0 && continue
        colname = aggregation == PSY.System ? "System" : PSY.get_name(agg)
        load_values = DataFrames.DataFrame()
        for load in loads
            f = PSY.get_time_series_values( # TODO: this isn't applying the scaling factors
                PSY.Deterministic,
                load,
                "max_active_power",
                start_time = initial_time,
                len = horizon,
            )
            load_values[:, PSY.get_name(load)] = f
        end
        parameters[Symbol(colname)] = load_values
    end
    time_range =
        range(initial_time, step = PSY.get_time_series_resolution(system), length = horizon)

    return PGData(parameters, time_range)
end

function get_service_data(results::R; kwargs...) where {R <: IS.Results}
    initial_time = get(kwargs, :initial_time, get(kwargs, :start_time, nothing))
    len = get(kwargs, :horizon, get(kwargs, :len, nothing))
    variable_keys = get(kwargs, :variable_keys, PSI.list_variable_keys(results))
    #parameter_keys = get(kwargs, :parameter_keys, PSI.list_parameter_keys(results))

    variable_keys = get_service_variable_keys(results; variable_keys = variable_keys)

    variables = PSI.read_variables_with_keys(
        results,
        variable_keys;
        start_time = initial_time,
        len = len,
    )

    timestamps = PSI.get_realized_timestamps(results; start_time = initial_time, len = len)

    return PGData(variables, timestamps)
end

#### result combination and aggregation ####

"""
aggregates and combines data into single DataFrame

# Example

```julia
PG.combine_categories(gen_uc.data)
```

"""
function combine_categories(
    data::Union{Dict{Symbol, DataFrames.DataFrame}, Dict{String, DataFrames.DataFrame}};
    names::Union{Vector{String}, Vector{Symbol}, Nothing} = nothing,
    aggregate::Union{Function, Nothing} = nothing,
)
    aggregate = isnothing(aggregate) ? x -> sum(x, dims = 2) : aggregate
    names = isnothing(names) ? keys(data) : names
    values = []
    keep_names = []
    for k in names
        if !isempty(data[k])
            push!(values, aggregate(Matrix(no_datetime(data[k]))))
            push!(keep_names, k)
        end
    end
    data = hcat(values...)
    keep_names = string.(keep_names)
    isempty(data) && return DataFrames.DataFrame()
    return DataFrames.DataFrame(data, keep_names)
end

"""
Re-categorizes data according to an aggregation dictionary
* makes no guarantee of complete data collection *

# Example

```julia
aggregation = PG.make_fuel_dictionary(results_uc.system)
categorize_data(gen_uc.data, aggregation)
```

"""
function categorize_data(
    data::Dict{Symbol, DataFrames.DataFrame},
    aggregation::Dict;
    curtailment = true,
    slacks = true,
)
    category_dataframes = Dict{String, DataFrames.DataFrame}()
    var_types = Dict(zip(last.(split.(string.(keys(data)), "_")), keys(data)))
    for (category, list) in aggregation
        category_df = DataFrames.DataFrame()
        for (component_type, variable) in list
            if haskey(var_types, component_type)
                category_data = data[var_types[component_type]]
                colname = typeof(names(category_data)[1]) == String ? "$variable" : variable
                DataFrames.insertcols!(
                    category_df,
                    (colname => category_data[:, colname]),
                    makeunique = true,
                )
            end
        end
        if !isempty(category_df)
            category_dataframes[string(category)] = category_df
        end
    end
    if curtailment
        dfs = []
        for (key, val) in data
            if endswith(string(key), "Curtailment")
                push!(dfs, no_datetime(val))
            end
        end
        if !isempty(dfs)
            category_dataframes["Curtailment"] = hcat(dfs...)
        end
    end
    if slacks
        for (slack, slack_name) in BALANCE_SLACKVARS
            for id in findall(x -> occursin(string(slack), x), string.(keys(data)))
                slack_key = collect(keys(data))[id]
                category_dataframes[slack_name] = data[slack_key]
            end
        end
    end

    return category_dataframes
end

max_stars_repo_name: NREL-SIIP/InfrastructureSystemsPlots

max_issues_repo_name: NREL-SIIP/InfrastructureSystemsPlots

max_forks_repo_name: NREL-SIIP/InfrastructureSystemsPlots


text:
"""
Start a notebook server using:

```julia
julia> Pluto.run()
```

Have a look at the FAQ:
https://github.com/fonsp/Pluto.jl/wiki
"""
module Pluto

import RelocatableFolders: @path
const ROOT_DIR = normpath(joinpath(@__DIR__, ".."))
const FRONTEND_DIR = @path(joinpath(ROOT_DIR, "frontend"))
const FRONTEND_DIST_DIR = let dir = joinpath(ROOT_DIR, "frontend-dist")
    isdir(dir) ? @path(dir) : FRONTEND_DIR
end
const frontend_dist_exists = FRONTEND_DIR !== FRONTEND_DIST_DIR
const SAMPLE_DIR = @path(joinpath(ROOT_DIR, "sample"))
const RUNNER_DIR = @path(joinpath(ROOT_DIR, "src", "runner"))
function project_relative_path(root, xs...)
    root == joinpath("src", "runner") ? joinpath(RUNNER_DIR, xs...) :
    root == "frontend-dist" && frontend_dist_exists ? joinpath(FRONTEND_DIST_DIR, xs...) :
    root == "frontend" ? joinpath(FRONTEND_DIR, xs...) :
    root == "sample" ? joinpath(SAMPLE_DIR, xs...) :
        normpath(joinpath(pkgdir(Pluto), root, xs...))
end

import Pkg

include_dependency("../Project.toml")
const PLUTO_VERSION = VersionNumber(Pkg.TOML.parsefile(joinpath(ROOT_DIR, "Project.toml"))["version"])
const PLUTO_VERSION_STR = 'v' * string(PLUTO_VERSION)
const JULIA_VERSION_STR = 'v' * string(VERSION)

include("./notebook/PathHelpers.jl")
include("./notebook/Export.jl")
include("./Configuration.jl")

include("./evaluation/Tokens.jl")
include("./evaluation/Throttled.jl")
include("./runner/PlutoRunner.jl")
include("./analysis/ExpressionExplorer.jl")
include("./analysis/FunctionDependencies.jl")
include("./analysis/ReactiveNode.jl")
include("./packages/PkgCompat.jl")

include("./notebook/Cell.jl")
include("./analysis/data structures.jl")
include("./analysis/Topology.jl")
include("./analysis/Errors.jl")
include("./analysis/TopologicalOrder.jl")
include("./notebook/Notebook.jl")
include("./notebook/frontmatter.jl")
include("./notebook/Events.jl")
include("./webserver/Session.jl")
include("./webserver/PutUpdates.jl")

include("./analysis/Parse.jl")
include("./analysis/topological_order.jl")
include("./analysis/is_just_text.jl")
include("./analysis/TopologyUpdate.jl")
include("./analysis/DependencyCache.jl")
include("./analysis/MoreAnalysis.jl")

include("./evaluation/WorkspaceManager.jl")
include("./packages/Packages.jl")
include("./packages/PkgUtils.jl")
include("./evaluation/Run.jl")
include("./evaluation/RunBonds.jl")

module DownloadCool include("./webserver/data_url.jl") end
include("./webserver/MsgPack.jl")
include("./webserver/SessionActions.jl")
include("./webserver/Static.jl")
include("./webserver/Dynamic.jl")
include("./webserver/REPLTools.jl")
include("./webserver/WebServer.jl")

const reset_notebook_environment = PkgUtils.reset_notebook_environment
const update_notebook_environment = PkgUtils.update_notebook_environment
const activate_notebook_environment = PkgUtils.activate_notebook_environment
export reset_notebook_environment
export update_notebook_environment
export activate_notebook_environment

if get(ENV, "JULIA_PLUTO_SHOW_BANNER", "1") != "0" && get(ENV, "CI", "🍄") != "true"
@info """\n
    Welcome to Pluto $(PLUTO_VERSION_STR) 🎈
    Start a notebook server using:

  julia> Pluto.run()

    Have a look at the FAQ:
    https://github.com/fonsp/Pluto.jl/wiki
\n"""
end

end