text:
"""
Function to create a unique index of time series names for each device model. For example,
if two parameters each reference the same time series name, this function will return a
different value for each parameter entry
"""
function create_time_series_multiplier_index(
    model,
    ::Type{T},
) where {T <: TimeSeriesParameter}
    ts_names = get_time_series_names(model)
    if length(ts_names) > 1
        ts_name = ts_names[T]
        ts_id = findfirst(x -> x == T, [k for (k, v) in ts_names if v == ts_name])
    else
        ts_id = 1
    end
    return ts_id
end

function add_parameters!(
    container::OptimizationContainer,
    ::Type{T},
    devices::U,
    model::DeviceModel{D, W},
) where {
    T <: TimeSeriesParameter,
    U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
    W <: AbstractDeviceFormulation,
} where {D <: PSY.Component}
    add_parameters!(container, T(), devices, model)
    return
end

function add_parameters!(
    container::OptimizationContainer,
    ::Type{T},
    devices::U,
    model::DeviceModel{D, W},
) where {
    T <: TimeSeriesParameter,
    U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
    W <: AbstractDeviceFormulation,
} where {D <: PSY.HybridSystem}
    _devices = [d for d in devices if PSY.get_renewable_unit(d) !== nothing]
    add_parameters!(container, T(), _devices, model)
    return
end

function add_parameters!(
    container::OptimizationContainer,
    ::Type{T},
    service::U,
    model::ServiceModel{U, V},
) where {T <: TimeSeriesParameter, U <: PSY.Service, V <: AbstractReservesFormulation}
    add_parameters!(container, T(), service, model)
    return
end

function add_parameters!(
    container::OptimizationContainer,
    ::Type{T},
    key::VariableKey{U, D},
    model::DeviceModel{D, W},
    devices::V,
) where {
    T <: VariableValueParameter,
    U <: VariableType,
    V <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
    W <: AbstractDeviceFormulation,
} where {D <: PSY.Component}
    add_parameters!(container, T(), key, model, devices)
    return
end

function add_parameters!(
    container::OptimizationContainer,
    ::T,
    devices::U,
    model::DeviceModel{D, W},
) where {
    T <: TimeSeriesParameter,
    U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
    W <: AbstractDeviceFormulation,
} where {D <: PSY.Component}
    ts_type = get_default_time_series_type(container)
    if !(ts_type <: Union{PSY.AbstractDeterministic, PSY.StaticTimeSeries})
        error("add_parameters! for TimeSeriesParameter is not compatible with $ts_type")
    end
    time_steps = get_time_steps(container)
    names = [PSY.get_name(d) for d in devices]
    ts_name = get_time_series_names(model)[T]
    time_series_mult_id = create_time_series_multiplier_index(model, T)
    @debug "adding" T ts_name ts_type time_series_mult_id _group =
        LOG_GROUP_OPTIMIZATION_CONTAINER
    parameter_container =
        add_param_container!(container, T(), D, ts_type, ts_name, names, time_steps)
    set_time_series_multiplier_id!(get_attributes(parameter_container), time_series_mult_id)
    jump_model = get_jump_model(container)
    for d in devices
        name = PSY.get_name(d)
        ts_vector = get_time_series(container, d, T())
        multiplier = get_multiplier_value(T(), d, W())
        for t in time_steps
            set_parameter!(
                parameter_container,
                jump_model,
                ts_vector[t],
                multiplier,
                name,
                t,
            )
        end
    end
    return
end

function add_parameters!(
    container::OptimizationContainer,
    ::T,
    devices::U,
    model::DeviceModel{D, W},
) where {
    T <: ActivePowerTimeSeriesParameter,
    U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
    W <: AbstractDeviceFormulation,
} where {D <: PSY.HybridSystem}
    ts_type = get_default_time_series_type(container)
    if !(ts_type <: Union{PSY.AbstractDeterministic, PSY.StaticTimeSeries})
        error("add_parameters! for TimeSeriesParameter is not compatible with $ts_type")
    end
    time_steps = get_time_steps(container)
    names = [PSY.get_name(d) for d in devices]
    ts_name = get_time_series_names(model)[T]
    time_series_mult_id = create_time_series_multiplier_index(model, T)
    @debug "adding" T ts_name ts_type time_series_mult_id _group =
        LOG_GROUP_OPTIMIZATION_CONTAINER
    sub_comp_type = [PSY.RenewableGen, PSY.ElectricLoad]
    parameter_container = add_param_container!(
        container,
        T(),
        D,
        ts_type,
        ts_name,
        names,
        string.(sub_comp_type),
        time_steps;
        sparse=true,
    )
    set_time_series_multiplier_id!(get_attributes(parameter_container), time_series_mult_id)
    jump_model = get_jump_model(container)
    for d in devices, comp_type in sub_comp_type
        name = PSY.get_name(d)
        if does_subcomponent_exist(d, comp_type)
            ts_vector = get_time_series(container, d, comp_type, T())
            multiplier = get_multiplier_value(T(), d, comp_type, W())
        else
            ts_vector = zeros(time_steps[end])
            multiplier = 0.0
        end
        for t in time_steps
            set_parameter!(
                parameter_container,
                jump_model,
                ts_vector[t],
                multiplier,
                name,
                string(comp_type),
                t,
            )
        end
    end
    return
end

function add_parameters!(
    container::OptimizationContainer,
    ::T,
    devices::U,
    model::DeviceModel{D, W},
) where {
    T <: ObjectiveFunctionParameter,
    U <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
    W <: AbstractDeviceFormulation,
} where {D <: PSY.Component}
    ts_type = get_default_time_series_type(container)
    if !(ts_type <: Union{PSY.AbstractDeterministic, PSY.StaticTimeSeries})
        error(
            "add_parameters! for ObjectiveFunctionParameter is not compatible with $ts_type",
        )
    end
    time_steps = get_time_steps(container)
    names = [PSY.get_name(d) for d in devices]
    ts_name = get_time_series_names(model)[T]
    time_series_mult_id = create_time_series_multiplier_index(model, T)
    @debug "adding" T ts_name ts_type time_series_mult_id _group =
        LOG_GROUP_OPTIMIZATION_CONTAINER
    parameter_container =
        add_param_container!(container, T(), D, ts_type, ts_name, names, time_steps)
    set_time_series_multiplier_id!(get_attributes(parameter_container), time_series_mult_id)
    jump_model = get_jump_model(container)
    for d in devices
        name = PSY.get_name(d)
        ts_vector = get_time_series(container, d, T())
        multiplier = get_multiplier_value(T(), d, W())
        for t in time_steps
            set_parameter!(
                parameter_container,
                jump_model,
                ts_vector[t],
                multiplier,
                name,
                t,
            )
        end
    end
    return
end

function add_parameters!(
    container::OptimizationContainer,
    ::Type{T},
    key::VariableKey{U, S},
    model::ServiceModel{S, W},
    devices::V,
) where {
    S <: PSY.AbstractReserve,
    T <: VariableValueParameter,
    U <: VariableType,
    V <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
    W <: AbstractReservesFormulation,
} where {D <: PSY.Component}
    add_parameters!(container, T(), key, model, devices)
    return
end

function add_parameters!(
    container::OptimizationContainer,
    ::T,
    service::U,
    model::ServiceModel{U, V},
) where {T <: TimeSeriesParameter, U <: PSY.Service, V <: AbstractReservesFormulation}
    ts_type = get_default_time_series_type(container)
    if !(ts_type <: Union{PSY.AbstractDeterministic, PSY.StaticTimeSeries})
        error("add_parameters! for TimeSeriesParameter is not compatible with $ts_type")
    end
    ts_name = get_time_series_names(model)[T]
    time_series_mult_id = create_time_series_multiplier_index(model, T)
    time_steps = get_time_steps(container)
    name = PSY.get_name(service)
    @debug "adding" parameter_type U _group = LOG_GROUP_OPTIMIZATION_CONTAINER
    parameter_container = add_param_container!(
        container,
        T(),
        U,
        ts_type,
        ts_name,
        [name],
        time_steps;
        meta=name,
    )
    set_time_series_multiplier_id!(get_attributes(parameter_container), time_series_mult_id)
    jump_model = get_jump_model(container)
    ts_vector = get_time_series(container, service, T(), name)
    multiplier = get_multiplier_value(T(), service, V())
    for t in time_steps
        set_parameter!(parameter_container, jump_model, ts_vector[t], multiplier, name, t)
    end

    return
end

function add_parameters!(
    container::OptimizationContainer,
    ::T,
    key::VariableKey{U, D},
    model::DeviceModel{D, W},
    devices::V,
) where {
    T <: VariableValueParameter,
    U <: VariableType,
    V <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
    W <: AbstractDeviceFormulation,
} where {D <: PSY.Component}
    @debug "adding" T D U _group = LOG_GROUP_OPTIMIZATION_CONTAINER
    names = [PSY.get_name(device) for device in devices]
    time_steps = get_time_steps(container)
    parameter_container = add_param_container!(container, T(), D, key, names, time_steps)
    jump_model = get_jump_model(container)

    for d in devices
        name = PSY.get_name(d)
        for t in time_steps
            set_parameter!(
                parameter_container,
                jump_model,
                get_initial_parameter_value(T(), d, W()),
                get_parameter_multiplier(T(), d, W()),
                name,
                t,
            )
        end
    end
    return
end

function add_parameters!(
    container::OptimizationContainer,
    ::T,
    key::AuxVarKey{U, D},
    model::DeviceModel{D, W},
    devices::V,
) where {
    T <: VariableValueParameter,
    U <: AuxVariableType,
    V <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
    W <: AbstractDeviceFormulation,
} where {D <: PSY.Component}
    @debug "adding" T D U _group = LOG_GROUP_OPTIMIZATION_CONTAINER
    names = [PSY.get_name(device) for device in devices]
    time_steps = get_time_steps(container)
    parameter_container = add_param_container!(container, T(), D, key, names, time_steps)
    jump_model = get_jump_model(container)

    for d in devices
        name = PSY.get_name(d)
        for t in time_steps
            set_parameter!(
                parameter_container,
                jump_model,
                get_initial_parameter_value(T(), d, W()),
                get_parameter_multiplier(T(), d, W()),
                name,
                t,
            )
        end
    end
    return
end

function add_parameters!(
    container::OptimizationContainer,
    ::T,
    devices::V,
    model::DeviceModel{D, W},
) where {
    T <: OnStatusParameter,
    V <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
    W <: AbstractDeviceFormulation,
} where {D <: PSY.Component}
    @debug "adding" T D V _group = LOG_GROUP_OPTIMIZATION_CONTAINER

    # We do this to handle cases where the same parameter is also added as a Feedforward.
    # When the OnStatusParameter is added without a feedforward it takes a Float value.
    # This is used to handle the special case of compact formulations.
    !isempty(get_feedforwards(model)) && return
    names = [PSY.get_name(device) for device in devices]
    time_steps = get_time_steps(container)
    parameter_container = add_param_container!(
        container,
        T(),
        D,
        VariableKey(OnVariable, D),
        names,
        time_steps,
    )
    jump_model = get_jump_model(container)

    for d in devices
        name = PSY.get_name(d)
        for t in time_steps
            set_parameter!(
                parameter_container,
                jump_model,
                get_initial_parameter_value(T(), d, W()),
                get_parameter_multiplier(T(), d, W()),
                name,
                t,
            )
        end
    end
    return
end

function add_parameters!(
    container::OptimizationContainer,
    ::T,
    key::VariableKey{U, S},
    model::ServiceModel{S, W},
    devices::V,
) where {
    S <: PSY.AbstractReserve,
    T <: VariableValueParameter,
    U <: VariableType,
    V <: Union{Vector{D}, IS.FlattenIteratorWrapper{D}},
    W <: AbstractReservesFormulation,
} where {D <: PSY.Component}
    @debug "adding" T D U _group = LOG_GROUP_OPTIMIZATION_CONTAINER
    names = [PSY.get_name(device) for device in devices]
    time_steps = get_time_steps(container)
    parameter_container = add_param_container!(container, T(), D, key, names, time_steps)
    jump_model = get_jump_model(container)

    for d in devices
        name = PSY.get_name(d)
        for t in time_steps
            set_parameter!(
                parameter_container,
                jump_model,
                get_initial_parameter_value(T(), S, W()),
                get_parameter_multiplier(T(), S, W()),
                name,
                t,
            )
        end
    end
    return
end


text:
@info("Precompiling headers...")
using Flux, Zygote, CUDA, Random
using Flux:@epochs
using Flux.Data:DataLoader
using Flux:onehotbatch
using Printf, BSON, Dates, DelimitedFiles, LaTeXStrings

include("data.jl")
include("train.jl")
include("eval.jl")
include("plot.jl")

# Parameters for data I/O
TOTAL_DATA = 473_800_775
TOTAL_NUM_DB = 474
DIM = 4
SAMPLE_SIZE = 1_000_000
DB_PATH = "/media/share/Dev/CalabiYau/data/polytopes_db_4d"
PATH = "/home/oppenheimer/Dev/calabiyau/trained"
NUMBERS_ONLY = true

# Parameters for model building

NUMBERS_ONLY ? MAX_SIZE = 1 : MAX_SIZE = 34;

ONTOLOGY = "euler"
if ONTOLOGY == "h11"
    ontology_index = 2
    CLASS_NUM = 480
elseif ONTOLOGY == "h21"
    ontology_index = 3
    CLASS_NUM = 480
elseif ONTOLOGY == "euler"
    ontology_index = 4
    CLASS_NUM = 961
end
if CLASS_NUM > 1
    CLASSES = collect(1:CLASS_NUM)
else
    # CLASSES = collect(-CLASS_NUM:2:960)
end
BATCH_SIZE = 1000
DATA_ID = "vector_data_with_features"

function __main__()
    ############################
    # For startng a new data set
    ############################
    # new_sample(export_path="$PATH/$DATA_ID", split=true)

    @info("Importing Data...")
    train_set, test_set = import_split("$PATH/$DATA_ID/data");

    # @info("Augmenting Data...")
    # augment!(train_set)
    # augment!(test_set)

    @info("Unpacking Data...")
    train_x, train_y = unpackage(train_set)
    test_x, test_y =   unpackage(test_set)

    ######################################
    # Debugging and verification
    sample = rand(1:size(train_set,1))
    display(train_set[sample][1])
    display(train_x[:,sample])
    ######################################

    if CLASS_NUM > 1 && ONTOLOGY == "euler"
        @info("Mapping ontology to classification...")
        train_y = euler_to_index(train_y)
        test_y = euler_to_index(test_y)
    end

    @info("Data loaded.")
    train_data = DataLoader(train_x, train_y, batchsize=BATCH_SIZE)
    global test_data = (test_x, test_y)

    @info("Building Model...")
    layers = 3
    neuron_density = 500
    hiddenLayer = Dense(neuron_density, neuron_density, relu)
    hidden(X) = Chain(hiddenLayer, X)
    hidden_layers = hidden(hiddenLayer)
    if l > 2
        for l ∈ 1:(layers - 2)
            hidden_layers = hidden(hidden_layers)
        end
    end
    model_label="$layers Hidden Layers w/ $neuron_density Neurons Per Layer"
    CLASS_NUM > 1 ? Output = softmax : Output = identity
    global model = Chain(
                   Dense(MAX_SIZE*DIM, neuron_density,relu),
                   hidden_layers,
                   Dense(neuron_density, CLASS_NUM, relu),
                   Output);

    @info("Passing to GPU...")
    model |> gpu
    train_data |> gpu
    test_data |> gpuds

    # sqnorm(x) = sum(abs2, x)
    opt = ADAM()
    acc(x,y) = sum(map(argmax, eachcol(model(x))) .== y)/size(y,1)
    if CLASS_NUM > 1
        loss(x,y) = Flux.crossentropy(model(x), Flux.onehotbatch(y, CLASSES)) #+ sum(sqnorm, Flux.params(model))
    else
        loss(x, y) = Flux.mse(model(x), Flux.onehotbatch(y,CLASSES))
    end

    Train(model, "$PATH/$DATA_ID", train_data, test_data, opt, loss, acc)

    model_dir = "$PATH/$DATA_ID/models/$ONTOLOGY/h11-2021-04-25T14:50:28.434_acc_37.141.bson"
    eval_data = [test_x[:,i] for i ∈ rand(1:size(test_x,2), 10)]
    R = LRP(model_dir, eval_data)
    R_max = map(maximum, eachcol(abs.(R)))
    R_norm = R
    for i = 1:10
        R_norm[:,i] = abs.(R_norm[:,i] / R_max[i])
    end

    xs = string.([test_y[i] - argmax(model(eval_data[i])) for i=1:length(eval_data)])
    ys = [1,2,3,4]
    heatmap(xs,ys,R_norm, c=cgrad(["navyblue", "white", "red"]), size=(340,125), dpi=600, xlabel=L"$h^{1,1}$", ylabel="Number of", colorbarlabel="Relevance", fontfamily="Computer Modern")
end

@info("Headers compiled.")

# __main__()

ext: jl


text:
### Power System
function psheader(system)
    basePower = [["Base Power" ""]; ["Power [MVA]" ""]]

    bus = [["Bus" "Type" "Demand" "Demand" "Shunt Conductance" "Shunt Susceptance" "Area" "Voltage" "Voltage" "Base Voltage" "Loss Zone" "Maximum Voltage" "Minimum Voltage" "Lagrange multiplier" "Lagrange multiplier" "Kuhn-Tucker multiplier" "Kuhn-Tucker multiplier"]
           ["Integer" "PQ(1), PV(2), Slack(3)" "Active Power [MW]" "Reactive Power [MVAr]" "Active Power [MW]" "Reactive Power [MVAr]" "Integer" "Magnitude [p.u.]" "Angle [deg]" "Magnitude [kV]" "Integer" "Magnitude [p.u.]" "Magnitude [p.u.]" "Active Power Mismatch [u/MW]" "Reactive Power Mismatch [u/MVAr]" "Upper Voltage Magnitude Limit [u/p.u.]" "Lower Voltage Magnitude Limit [u/p.u.]"]]

    branch = [["Branch" "From Bus" "To Bus" "Series Parameter" "Series Parameter" "Charging Parameter" "Long Term Rate" "Short Term Rate" "Emergency Rate" "Transformer" "Transformer" "Status" "Minimum Voltage Difference" "Maximum Voltage Difference" "Injected From Bus" "Injected From Bus" "Injected To Bus" "Injected To Bus" "Kuhn-Tucker Multiplier" "Kuhn-Tucker Multiplier" "Kuhn-Tucker Multiplier" "Kuhn-Tucker Multiplier"]
              ["Integer" "Integer" "Integer" "Resistance [p.u.]" "Reactance [p.u.]" "Susceptance [p.u.]" "Power [MVA]" "Power [MVA]" "Power [MVA]" "Turns Ratio" "Shift Angle [deg]" "Integer" "Angle [deg]" "Angle [deg]" "Active Power [MW]" "Reactive Power [MVAr]" "Active Power [MW]" "Reactive Power [MVAr]" "Power From Bus Limit [u/MVA]" "Power To Bus Limit [u/MVA]" "Lower Angle Difference Limit [u/deg]" "Upper Angle Difference Limit [u/deg]"]]

    generator = [["Bus" "Generation" "Generation" "Maximum Generation" "Minimum Generation" "Voltage" "Base" "Status" "Maximum Generation" "Minimum Generation" "Lower of PQ Curve" "Upper of PQ Curve" "Minimum at PC1" "Maximum at PC1" "Minimum at PC2" "Maximum at PC2" "Ramp Rate ACG" "Ramp Rate 10" "Ramp Rate 30" "Ramp Rate Q" "Area Factor" "Kuhn-Tucker Multiplier" "Kuhn-Tucker Multiplier" "Kuhn-Tucker Multiplier" "Kuhn-Tucker Multiplier"]
                  ["Integer" "Active Power [MW]" "Reactive Power [MVAr]" "Reactive Power [MVAr]" "Reactive Power [MVAr]" "Magnitude [p.u.]" "Power [MVA]" "Integer" "Active Power [MW]" "Active Power [MW]" "Active Power [MW]" "Active Power [MW]" "Reactive Power [MVAr]" "Reactive Power [MVAr]" "Reactive Power [MVAr]" "Reactive Power [MVAr]" "Active Power per Minut [MW/min]" "Active Power [MW]" "Active Power [MW]" "Reactive Power per Minut [MVAr/min]" "Integer" "Upper Genration Active Power Limit [u/MW]" "Lower Genration Active Power Limit [u/MW]" "Upper Genration Reactive Power Limit [u/MVAr]" "Lower Genration Reactive Power Limit [u/MVAr]"]]

    N = size(system.generatorcost, 2) - 4
    generatorcost = [["Cost Model" "Cost" "Cost" "Cost Model" hcat(["Cost Model" for n = 1:N]...)]
                     ["Piecewise(1), Polynomial(2)" "Startup [currency]" "Shutdown [currency]" "Number of Data Points" hcat(["Coefficient c$(trunc(Int, (n-1)))" for n = N:-1:1]...)]]

    return (bus = bus, branch = branch, generator = generator, generatorcost = generatorcost, basePower = basePower)
end

### DC Power Flow
function dcpfheader()
    main = [["Bus" "Voltage" "Injection" "Generation" "Demand" "Shunt"]
            ["Integer" "Angle [deg]" "Active Power [MW]" "Active Power [MW]" "Active Power [MW]" "Active Power [MW]"]]

    flow = [["Branch" "From Bus" "To Bus" "From Bus Flow" "To Bus Flow"]
            ["Integer" "Integer" "Integer" "Active Power [MW]" "Active Power [MW]"]]

    generation = [["Bus" "Generation"]; ["Integer" "Active Power [MW]"]]

    return (main = main, flow = flow, generation = generation)
end

### AC Power Flow
function acpfheader()
    main = [["Bus" "Voltage" "Voltage" "Injection" "Injection" "Generation" "Generation" "Demand" "Demand" "Shunt" "Shunt"]
            ["Integer" "Magnitude [p.u.]" "Angle [deg]" "Active Power [MW]" "Reactive Power [MVAr]" "Active Power [MW]" "Reactive Power [MVAr]" "Active Power [MW]" "Reactive Power [MVAr]" "Active Power [MW]" "Reactive Power [MVAr]"]]

    flow = [["Branch" "From Bus" "To Bus" "From Bus Flow" "From Bus Flow" "To Bus Flow" "To Bus Flow" "Branch Injection" "Loss" "Loss" "From Bus Current" "From Bus Current" "To Bus Current" "To Bus Current"]
            ["Integer" "Integer" "Integer" "Active Power [MW]" "Reactive Power [MVAr]" "Active Power [MW]" "Reactive Power [MVAr]" "Reactive Power [MVAr]" "Active Power [MW]" "Reactive Power [MVAr]" "Magnitude [p.u.]" "Angle [deg]" "Magnitude [p.u.]" "Angle [deg]"]]

    generation = [["Bus" "Generation" "Generation"]; ["Integer" "Active Power [MW]" "Reactive Power [MVAr]"]]

    return (main = main, flow = flow, generation = generation)
end

### Measurements
function measureheader()
    pmuVoltage = [["Bus" "Measurement" "Noise" "Status" "Measurement" "Noise" "Status" "Exact" "Exact"]
                  ["Integer" "Magnitude [p.u.]" "Variance [p.u.]" "Integer" "Angle [rad]" "Variance [p.u.]" "Integer" "Magnitude [p.u.]" "Angle [rad]"]]

    pmuCurrent = [["Branch" "From Bus" "To Bus" "Measurement" "Noise" "Status" "Measurement" "Noise" "Status" "Exact" "Exact"]
                  ["Integer" "Integer" "Integer" "Magnitude [p.u.]" "Variance [p.u.]" "Integer" "Angle [rad]" "Variance [p.u.]" "Integer" "Magnitude [p.u.]" "Angle [rad]"]]

    legacyFlow = [["Branch" "From Bus" "To Bus" "Measurement" "Noise" "Status" "Measurement" "Noise" "Status" "Exact" "Exact"]
                  ["Integer" "Integer" "Integer" "Active Power [p.u.]" "Variance [p.u.]" "Integer" "Reactive Power [p.u.]" "Variance [p.u.]" "Integer" "Active Power [p.u.]" "Reactive Power [p.u.]"]]

    legacyCurrent = [["Branch" "From Bus" "To Bus" "Measurement" "Noise" "Status" "Exact"]
                     ["Integer" "Integer" "Integer" "Magnitude [p.u.]" "Variance [p.u.]" "Integer" "Magnitude [p.u.]"]]

    legacyInjection = [["Bus" "Measurement" "Noise" "Status" "Measurement" "Noise" "Status" "Exact" "Exact"]
                       ["Integer" "Active Power [p.u.]" "Variance [p.u.]" "Integer" "Reactive Power [p.u.]" "Variance [p.u.]" "Integer" "Active Power [p.u.]" "Reactive Power [p.u.]"]]

    legacyVoltage = [["Bus" "Measurement" "Noise" "Status" "Exact"]
                     ["Integer" "Magnitude [p.u.]" "Variance [p.u.]" "Integer" "Magnitude [p.u.]"]]

    return (pmuVoltage = pmuVoltage, pmuCurrent = pmuCurrent, legacyFlow = legacyFlow, legacyCurrent = legacyCurrent, legacyInjection = legacyInjection, legacyVoltage = legacyVoltage)
end

function dcseheader()
    main = [["Bus" "Voltage" "Injection"]
            ["Integer" "Angle [deg]" "Active Power [MW]"]]

    flow = [["Branch" "From Bus" "To Bus" "From Bus Flow" "To Bus Flow"]
           ["Integer" "Integer" "Integer" "Active Power [MW]" "Active Power [MW]"]]

    estimate = [["Number" "Device" "Device" "Device" "Device" "Device" "Device" "Algorithm" "Residual" "User" "Residual" ]
                ["Integer" "In-service(1), Bad Data(2), Pseudo(3)" "Legacy(1), PMU(2)" "Pij(1), Pi(4), Ti(8)" "Local Index" "Measure" "Variance" "Estimate" "Estimate to Measure" "Exact" "Estimate to Exact"]]

    estimatedisplay = [["Device" "Device" "Device" "Device" "Device" "Device" "Algorithm" "Residual" "User" "Residual" ]
                       ["Status" "Class" "Type" "Local Index" "Measure" "Variance" "Estimate" "Estimate to Measure" "Exact" "Estimate to Exact"]]

    error = [["MAE" "RMS" "WRSSE" "MAE" "RMSE" "WRSSE"]
            ["estimates to measurements" "estimates to measurements" "estimates to measurements" "estimates to exacts"  "estimates to exacts"  "estimates to exacts"]]

    errordisplay = [["Estimate and corresponding measurement values in per-unit system" ""]
                    ["Mean absolute error" "."^10]
                    ["Root mean square error" "."^10]
                    ["Weighted residual sum of squares error" "."^10]
                    ["" ""]
                    ["Estimate and corresponding exact values in per-unit system" ""]
                    ["Mean absolute error" "."^10]
                    ["Root mean square error" "."^10]
                    ["Weighted residual sum of squares error" "."^10]]

    bad = [["Algorithm" "Device" "Device" "Device"  "Algorithm" "Device"]
           ["Pass" "Class" "Suspected Bad Data" "Local Index" "Normalized Residual" "Status"]]

    observe  = [["Island" "Bus in Island"  "Device" "Device" "Device" "Device" "Device"]
                ["Integer" "Integer" "Class" "Pseudo-measurement" "Local Index" "Measure" "Variance"]]

    return (main = main, flow = flow, estimatedisplay = estimatedisplay, estimate = estimate, error = error, errordisplay = errordisplay, bad = bad, observe = observe)
end

### AC Power Flow
function pmuseheader()
    main = [["Bus" "Voltage" "Voltage" "Injection" "Injection" "Shunt" "Shunt"]
            ["Integer" "Magnitude [p.u.]" "Angle [deg]" "Active Power [MW]" "Reactive Power [MVAr]" "Active Power [MW]" "Reactive Power [MVAr]"]]

    flow = [["Branch" "From Bus" "To Bus" "From Bus Flow" "From Bus Flow" "To Bus Flow" "To Bus Flow" "Branch Injection" "Loss" "Loss" "From Bus Current" "From Bus Current" "To Bus Current" "To Bus Current"]
            ["Integer" "Integer" "Integer" "Active Power [MW]" "Reactive Power [MVAr]" "Active Power [MW]" "Reactive Power [MVAr]" "Reactive Power [MVAr]" "Active Power [MW]" "Reactive Power [MVAr]" "Magnitude [p.u.]" "Angle [deg]" "Magnitude [p.u.]" "Angle [deg]"]]

    estimate = [["Number" "Device" "Device" "Device" "Device" "Device" "Device" "Algorithm" "Residual" "User" "Residual" ]
                ["Integer" "In-service(1), Bad Data(2), Pseudo(3)" "Legacy(1), PMU(2)" "Pij(1), Pi(4), Ti(8)" "Local Index" "Measure" "Variance" "Estimate" "Estimate to Measure" "Exact" "Estimate to Exact"]]

    estimatedisplay = [["Device" "Device" "Device" "Device" "Device" "Device" "Algorithm" "Residual" "User" "Residual" ]
                       ["Status" "Class" "Type" "Local Index" "Measure" "Variance" "Estimate" "Estimate to Measure" "Exact" "Estimate to Exact"]]

    error = [["MAE" "RMS" "WRSSE" "MAE" "RMSE" "WRSSE"]
            ["estimates to measurements" "estimates to measurements" "estimates to measurements" "estimates to exacts"  "estimates to exacts"  "estimates to exacts"]]

    errordisplay = [["Estimate and corresponding measurement values in per-unit system" ""]
                    ["Mean absolute error" "."^10]
                    ["Root mean square error" "."^10]
                    ["Weighted residual sum of squares error" "."^10]
                    ["" ""]
                    ["Estimate and corresponding exact values in per-unit system" ""]
                    ["Mean absolute error" "."^10]
                    ["Root mean square error" "."^10]
                    ["Weighted residual sum of squares error" "."^10]]

    bad = [["Algorithm" "Device" "Device" "Device"  "Algorithm" "Device"]
            ["Pass" "Class" "Suspected Bad Data" "Local Index" "Normalized Residual" "Status"]]

    observe  = [["Flow Island" "Bus in Island"  "Device" "Device" "Device" "Device" "Device"]
                ["Integer" "Integer" "Class" "Pseudo-measurement" "Local Index" "Measure" "Variance"]]

    return (main = main, flow = flow, estimatedisplay = estimatedisplay, estimate = estimate, error = error, errordisplay = errordisplay, bad = bad, observe = observe)
end

### AC Power Flow
function acseheader()
    main = [["Bus" "Voltage" "Voltage" "Injection" "Injection" "Shunt" "Shunt"]
            ["Integer" "Magnitude [p.u.]" "Angle [deg]" "Active Power [MW]" "Reactive Power [MVAr]" "Active Power [MW]" "Reactive Power [MVAr]"]]

    flow = [["Branch" "From Bus" "To Bus" "From Bus Flow" "From Bus Flow" "To Bus Flow" "To Bus Flow" "Branch Injection" "Loss" "Loss" "From Bus Current" "From Bus Current" "To Bus Current" "To Bus Current"]
            ["Integer" "Integer" "Integer" "Active Power [MW]" "Reactive Power [MVAr]" "Active Power [MW]" "Reactive Power [MVAr]" "Reactive Power [MVAr]" "Active Power [MW]" "Reactive Power [MVAr]" "Magnitude [p.u.]" "Angle [deg]" "Magnitude [p.u.]" "Angle [deg]"]]

    estimate = [["Number" "Device" "Device" "Device" "Device" "Device" "Device" "Algorithm" "Residual" "User" "Residual" ]
                ["Integer" "In-service(1), Bad Data(2), Pseudo(3)" "Legacy(1), PMU(2)" "Pij(1), Pi(4), Ti(8)" "Local Index" "Measure" "Variance" "Estimate" "Estimate to Measure" "Exact" "Estimate to Exact"]]

    estimatedisplay = [["Device" "Device" "Device" "Device" "Device" "Device" "Algorithm" "Residual" "User" "Residual" ]
                        ["Status" "Class" "Type" "Local Index" "Measure" "Variance" "Estimate" "Estimate to Measure" "Exact" "Estimate to Exact"]]

    error = [["MAE" "RMS" "WRSSE" "MAE" "RMSE" "WRSSE"]
            ["estimates to measurements" "estimates to measurements" "estimates to measurements" "estimates to exacts"  "estimates to exacts"  "estimates to exacts"]]

    errordisplay = [["Estimate and corresponding measurement values in per-unit system" ""]
                    ["Mean absolute error" "."^10]
                    ["Root mean square error" "."^10]
                    ["Weighted residual sum of squares error" "."^10]
                    ["" ""]
                    ["Estimate and corresponding exact values in per-unit system" ""]
                    ["Mean absolute error" "."^10]
                    ["Root mean square error" "."^10]
                    ["Weighted residual sum of squares error" "."^10]]

    bad = [["Algorithm" "Device" "Device" "Device"  "Algorithm" "Device"]
            ["Pass" "Class" "Suspected Bad Data" "Local Index" "Normalized Residual" "Status"]]

    observe  = [["Island" "Bus in Island"  "Device" "Device" "Device" "Device" "Device"]
                ["Integer" "Integer" "Class" "Pseudo-measurement" "Local Index" "Measure" "Variance"]]

    return (main = main, flow = flow, estimatedisplay = estimatedisplay, estimate = estimate, error = error, errordisplay = errordisplay, bad = bad, observe = observe)
end

max_stars_repo_name: mcosovic/JuliaGrid

max_issues_repo_name: mcosovic/JuliaGrid

max_forks_repo_name: mcosovic/JuliaGrid


text:
using SIIPExamples
pkgpath = dirname(dirname(pathof(SIIPExamples)))
include(joinpath(pkgpath, "test", "3_PowerSimulations_examples", "01_operations_problems.jl"));

using Ipopt
solver = optimizer_with_attributes(Ipopt.Optimizer)

ed_template = template_economic_dispatch(network = ACPPowerModel)

ed_template.devices[:Hydro] = DeviceModel(HydroEnergyReservoir, HydroDispatchRunOfRiver)

problem = OperationsProblem(
    EconomicDispatchProblem,
    ed_template,
    sys,
    horizon = 4,
    optimizer = solver,
    balance_slack_variables = true,
)

solve!(problem)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
## much of this is modified from legacy Amanda Mason code and appears to have a lot of
## unnecessary and/or redundant code, JJS, 12/5/19

"""
This splits the wntr dictionary (links, nodes) into its sub compononents
"""
function wntr_dict(inp_file::String)
    wn = model.WaterNetworkModel(inp_file)
    wns =  sim.EpanetSimulator(wn)
    results = wns.run_sim()  # why run a wntr simulation? JJS 12/5/19
    # dictionary of head, pressure, demand, quality
    node_results = results.node
    # dictionary of status, flowrate, velocity, {headloss, setting, friciton factor,
    # reaction rate, link quality} = epanet simulator only
    link_results = results.link
    wntr_dict = wn.todict()
    node_dicts(wntr_dict)
    link_dicts(wntr_dict)
    curves_dict(wntr_dict)
    # why have these results? are they useful for constructing the julia-version of a
    # WaterSystem? JJS 12/5/19
    wntr_dict["node_results"] = node_results
    wntr_dict["link_results"] = link_results
    return wntr_dict
end

# these functions should all have '!' appended to their names, JJS 12/9/19
function link_dicts(wntr_dict::Dict{Any, Any})
    links = wntr_dict["links"]
    # why a dict for pipes, but a vector for pumps? JJS 12/11/19
    links_vec = Vector{Any}()
    pipe_dict = Dict{Int64,Any}()
    valves = Vector{Any}()
    pumps = Vector{Any}()
    control_pipe = Array{Union{Nothing,String}}(nothing,0)
    for (i, (name,control)) in enumerate(wntr_dict["controls"])
        target = control._then_actions[1]._target_obj
        target.link_type == "Pipe" ? push!(control_pipe, target.name) : nothing
    end
    l = 0
    for link in links
        push!(links_vec, link)
        if link["link_type"] == "Pipe"
            l = l + 1
            pipe_dict[l] = link
            pipe_dict[l]["control_pipe"] = link["name"] in control_pipe
        elseif link["link_type"] == "Pump"
            push!(pumps, link)
        else
            push!(valves,link)
        end
    end
    wntr_dict["pipes"] = pipe_dict
    wntr_dict["valves"] = valves
    wntr_dict["pumps"] = pumps
    wntr_dict["links_vec"] = links_vec # container for all links, JJS 12/11/19
end

function node_dicts(wntr_dict::Dict{Any,Any})
    nodes = wntr_dict["nodes"]
    junctions= Vector{Any}()
    tanks =  Vector{Any}() # TODO: Does this need to be a dict?
    reservoirs= Vector{Any}() # TODO: Does this need to be a dict?
    for node in nodes
        if node["node_type"] == "Junction"
            push!(junctions, node)
        elseif node["node_type"] == "Tank"
            push!(tanks, node)
        else
            push!(reservoirs, node)
        end
    end
    wntr_dict["junctions"] = junctions
    wntr_dict["tanks"] = tanks
    wntr_dict["reservoirs"] = reservoirs
end

function curves_dict(wntr_dict::Dict{Any,Any})
    curves = wntr_dict["curves"]
    curves_dict = Dict{String,Any}()
    for curve in curves
        curves_dict[curve["name"]] = curve
    end
    delete!(wntr_dict, "curves")
    wntr_dict["curves"] = curves_dict
end