text:
const DATA_DIR = joinpath(BASE_DIR, "test", "test_data")
const TEST_OUTPUTS = joinpath(BASE_DIR, "test", "test_results")
!isdir(TEST_OUTPUTS) && mkdir(TEST_OUTPUTS)
const TEST_RESULT_DIR = joinpath(TEST_OUTPUTS, "results")
!isdir(TEST_RESULT_DIR) && mkdir(TEST_RESULT_DIR)

include(joinpath(DATA_DIR, "data_5bus_pu.jl"))
include(joinpath(DATA_DIR, "data_14bus_pu.jl"))
include(joinpath(DATA_DIR, "results_data.jl"))

const TS = TimeSeries
# Test Systems

# The code below provides a mechanism to optimally construct test systems. The first time a
# test builds a particular system name, the code will construct the system from raw files
# and then serialize it to storage.
# When future tests ask for the same system the code will deserialize it from storage.
#
# If you add a new system then you need to add an entry to TEST_SYSTEMS.
# The build function should accept `kwargs...` instead of specific named keyword arguments.
# This will allow easy addition of new parameters in the future.

struct TestSystemLabel
    name::String
    add_forecasts::Bool
    add_reserves::Bool
end

mutable struct SystemBuildStats
    count::Int
    initial_construct_time::Float64
    serialize_time::Float64
    min_deserialize_time::Float64
    max_deserialize_time::Float64
    total_deserialize_time::Float64
end

function SystemBuildStats(initial_construct_time::Float64, serialize_time::Float64)
    return SystemBuildStats(1, initial_construct_time, serialize_time, 0.0, 0.0, 0.0)
end

function update_stats!(stats::SystemBuildStats, deserialize_time::Float64)
    stats.count += 1
    if stats.min_deserialize_time == 0 || deserialize_time < stats.min_deserialize_time
        stats.min_deserialize_time = deserialize_time
    end
    if deserialize_time > stats.max_deserialize_time
        stats.max_deserialize_time = deserialize_time
    end
    stats.total_deserialize_time += deserialize_time
end

avg_deserialize_time(stats::SystemBuildStats) = stats.total_deserialize_time / stats.count

g_system_serialized_files = Dict{TestSystemLabel, String}()
g_system_build_stats = Dict{TestSystemLabel, SystemBuildStats}()

function initialize_system_serialized_files()
    empty!(g_system_serialized_files)
    empty!(g_system_build_stats)
end

function summarize_system_build_stats()
    @info "System Build Stats"
    labels = sort!(collect(keys(g_system_build_stats)), by = x -> x.name)
    for label in labels
        x = g_system_build_stats[label]
        system = "$(label.name) add_forecasts=$(label.add_forecasts) add_reserves=$(label.add_reserves)"
        @info system x.count x.initial_construct_time x.serialize_time x.min_deserialize_time x.max_deserialize_time avg_deserialize_time(
            x,
        )
    end
end

function build_system(name::String; add_forecasts = true, add_reserves = false)
    !haskey(TEST_SYSTEMS, name) && error("invalid system name: $name")
    label = TestSystemLabel(name, add_forecasts, add_reserves)
    sys_params = TEST_SYSTEMS[name]
    if !haskey(g_system_serialized_files, label)
        @debug "Build new system" label sys_params.description
        build_func = sys_params.build
        start = time()
        sys = build_func(;
            add_forecasts = add_forecasts,
            add_reserves = add_reserves,
            time_series_in_memory = sys_params.time_series_in_memory,
        )
        construct_time = time() - start
        serialized_file = joinpath(mktempdir(), "sys.json")
        start = time()
        PSY.to_json(sys, serialized_file)
        serialize_time = time() - start
        g_system_build_stats[label] = SystemBuildStats(construct_time, serialize_time)
        g_system_serialized_files[label] = serialized_file
    else
        @debug "Deserialize system from file" label
        start = time()
        sys = System(
            g_system_serialized_files[label];
            time_series_in_memory = sys_params.time_series_in_memory,
        )
        update_stats!(g_system_build_stats[label], time() - start)
    end

    return sys
end

function build_c_sys5(; kwargs...)
    nodes = nodes5()
    c_sys5 =
        System(100.0, nodes, thermal_generators5(nodes), loads5(nodes), branches5(nodes))

    if get(kwargs, :add_forecasts, true)
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys5))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(load_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = load_timeseries_DA[t][ix]
            end
            add_time_series!(c_sys5, l, Deterministic("max_active_power", forecast_data))
        end
    end

    return c_sys5
end

function build_c_sys5_ml(; kwargs...)
    nodes = nodes5()
    c_sys5_ml = System(
        100.0,
        nodes,
        thermal_generators5(nodes),
        loads5(nodes),
        branches5(nodes);
        time_series_in_memory = get(kwargs, :time_series_in_memory, true),
    )

    if get(kwargs, :add_forecasts, true)
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_ml))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(load_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = load_timeseries_DA[t][ix]
            end
            add_time_series!(c_sys5_ml, l, Deterministic("max_active_power", forecast_data))
        end
    end

    return c_sys5_ml
end

function build_c_sys14(; kwargs...)
    nodes = nodes14()
    c_sys14 = System(
        100.0,
        nodes,
        thermal_generators14(nodes),
        loads14(nodes),
        branches14(nodes);
        time_series_in_memory = get(kwargs, :time_series_in_memory, true),
    )

    if get(kwargs, :add_forecasts, true)
        forecast_data = SortedDict{Dates.DateTime, TimeArray}()
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys14))
            ini_time = TS.timestamp(timeseries_DA14[ix])[1]
            forecast_data[ini_time] = timeseries_DA14[ix]
            add_time_series!(c_sys14, l, Deterministic("max_active_power", forecast_data))
        end
    end

    return c_sys14
end

function build_c_sys5_re(; kwargs...)
    nodes = nodes5()
    c_sys5_re = System(
        100.0,
        nodes,
        thermal_generators5(nodes),
        renewable_generators5(nodes),
        loads5(nodes),
        branches5(nodes);
        time_series_in_memory = get(kwargs, :time_series_in_memory, true),
    )

    if get(kwargs, :add_forecasts, true)
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_re))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(load_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = load_timeseries_DA[t][ix]
            end
            add_time_series!(c_sys5_re, l, Deterministic("max_active_power", forecast_data))
        end
        for (ix, r) in enumerate(get_components(RenewableGen, c_sys5_re))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(ren_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = ren_timeseries_DA[t][ix]
            end
            add_time_series!(c_sys5_re, r, Deterministic("max_active_power", forecast_data))
        end
    end

    if get(kwargs, :add_reserves, false)
        reserve_re = reserve5_re(get_components(RenewableDispatch, c_sys5_re))
        add_service!(c_sys5_re, reserve_re[1], get_components(RenewableDispatch, c_sys5_re))
        add_service!(
            c_sys5_re,
            reserve_re[2],
            [collect(get_components(RenewableDispatch, c_sys5_re))[end]],
        )
        # ORDC
        add_service!(c_sys5_re, reserve_re[3], get_components(RenewableDispatch, c_sys5_re))
        for (ix, serv) in enumerate(get_components(VariableReserve, c_sys5_re))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(Reserve_ts[t])[1]
                forecast_data[ini_time] = Reserve_ts[t]
            end
            add_time_series!(c_sys5_re, serv, Deterministic("requirement", forecast_data))
        end
        for (ix, serv) in enumerate(get_components(ReserveDemandCurve, c_sys5_re))
            forecast_data = SortedDict{Dates.DateTime, Vector{IS.PWL}}()
            for t in 1:2
                ini_time = TS.timestamp(ORDC_cost_ts[t])[1]
                forecast_data[ini_time] = TimeSeries.values(ORDC_cost_ts[t])
            end
            resolution = TS.timestamp(ORDC_cost_ts[1])[2] - TS.timestamp(ORDC_cost_ts[1])[1]
            set_variable_cost!(
                c_sys5_re,
                serv,
                Deterministic("variable_cost", forecast_data, resolution),
            )
        end
    end

    return c_sys5_re
end

function build_c_sys5_re_only(; kwargs...)
    nodes = nodes5()
    c_sys5_re_only = System(
        100.0,
        nodes,
        renewable_generators5(nodes),
        loads5(nodes),
        branches5(nodes);
        time_series_in_memory = get(kwargs, :time_series_in_memory, true),
    )

    if get(kwargs, :add_forecasts, true)
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_re_only))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(load_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = load_timeseries_DA[t][ix]
            end
            add_time_series!(
                c_sys5_re_only,
                l,
                Deterministic("max_active_power", forecast_data),
            )
        end
        for (ix, r) in enumerate(get_components(RenewableGen, c_sys5_re_only))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(ren_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = ren_timeseries_DA[t][ix]
            end
            add_time_series!(
                c_sys5_re_only,
                r,
                Deterministic("max_active_power", forecast_data),
            )
        end
    end

    return c_sys5_re_only
end

function build_c_sys5_hy(; kwargs...)
    nodes = nodes5()
    c_sys5_hy = System(
        100.0,
        nodes,
        thermal_generators5(nodes),
        [hydro_generators5(nodes)[1]],
        loads5(nodes),
        branches5(nodes);
        time_series_in_memory = get(kwargs, :time_series_in_memory, true),
    )

    if get(kwargs, :add_forecasts, true)
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_hy))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(load_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = load_timeseries_DA[t][ix]
            end
            add_time_series!(c_sys5_hy, l, Deterministic("max_active_power", forecast_data))
        end
        for (ix, r) in enumerate(get_components(HydroGen, c_sys5_hy))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(hydro_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = hydro_timeseries_DA[t][ix]
            end
            add_time_series!(c_sys5_hy, r, Deterministic("max_active_power", forecast_data))
        end
    end

    return c_sys5_hy
end

function build_c_sys5_hyd(; kwargs...)
    nodes = nodes5()
    c_sys5_hyd = System(
        100.0,
        nodes,
        thermal_generators5(nodes),
        [hydro_generators5(nodes)[2]],
        loads5(nodes),
        branches5(nodes);
        time_series_in_memory = get(kwargs, :time_series_in_memory, true),
    )

    if get(kwargs, :add_forecasts, true)
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_hyd))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(load_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = load_timeseries_DA[t][ix]
            end
            add_time_series!(
                c_sys5_hyd,
                l,
                Deterministic("max_active_power", forecast_data),
            )
        end
        for (ix, h) in enumerate(get_components(HydroGen, c_sys5_hyd))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(hydro_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = hydro_timeseries_DA[t][ix]
            end
            add_time_series!(
                c_sys5_hyd,
                h,
                Deterministic("max_active_power", forecast_data),
            )
        end
        for (ix, h) in enumerate(get_components(HydroEnergyReservoir, c_sys5_hyd))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(hydro_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = hydro_timeseries_DA[t][ix]
            end
            add_time_series!(c_sys5_hyd, h, Deterministic("hydro_budget", forecast_data))
        end
        for (ix, h) in enumerate(get_components(HydroEnergyReservoir, c_sys5_hyd))
            forecast_data_inflow = SortedDict{Dates.DateTime, TimeArray}()
            forecast_data_target = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(hydro_timeseries_DA[t][ix])[1]
                forecast_data_inflow[ini_time] = hydro_timeseries_DA[t][ix] .* 0.8
                forecast_data_target[ini_time] = hydro_timeseries_DA[t][ix] .* 0.5
            end
            add_time_series!(c_sys5_hyd, h, Deterministic("inflow", forecast_data_inflow))
            add_time_series!(
                c_sys5_hyd,
                h,
                Deterministic("storage_target", forecast_data_target),
            )
        end
    end

    if get(kwargs, :add_reserves, false)
        reserve_hy = reserve5_hy(get_components(HydroEnergyReservoir, c_sys5_hyd))
        add_service!(
            c_sys5_hyd,
            reserve_hy[1],
            get_components(HydroEnergyReservoir, c_sys5_hyd),
        )
        add_service!(
            c_sys5_hyd,
            reserve_hy[2],
            [collect(get_components(HydroEnergyReservoir, c_sys5_hyd))[end]],
        )
        # ORDC curve
        add_service!(
            c_sys5_hyd,
            reserve_hy[3],
            get_components(HydroEnergyReservoir, c_sys5_hyd),
        )
        for (ix, serv) in enumerate(get_components(VariableReserve, c_sys5_hyd))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(Reserve_ts[t])[1]
                forecast_data[ini_time] = Reserve_ts[t]
            end
            add_time_series!(c_sys5_hyd, serv, Deterministic("requirement", forecast_data))
        end
        for (ix, serv) in enumerate(get_components(ReserveDemandCurve, c_sys5_hyd))
            forecast_data = SortedDict{Dates.DateTime, Vector{IS.PWL}}()
            for t in 1:2
                ini_time = TS.timestamp(ORDC_cost_ts[t])[1]
                forecast_data[ini_time] = TimeSeries.values(ORDC_cost_ts[t])
            end
            resolution = TS.timestamp(ORDC_cost_ts[1])[2] - TS.timestamp(ORDC_cost_ts[1])[1]
            set_variable_cost!(
                c_sys5_hyd,
                serv,
                Deterministic("variable_cost", forecast_data, resolution),
            )
        end
    end

    return c_sys5_hyd
end

function build_c_sys5_bat(; kwargs...)
    time_series_in_memory = get(kwargs, :time_series_in_memory, true)
    nodes = nodes5()
    c_sys5_bat = System(
        100.0,
        nodes,
        thermal_generators5(nodes),
        renewable_generators5(nodes),
        loads5(nodes),
        branches5(nodes),
        battery5(nodes);
        time_series_in_memory = time_series_in_memory,
    )

    if get(kwargs, :add_forecasts, true)
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_bat))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(load_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = load_timeseries_DA[t][ix]
            end
            add_time_series!(
                c_sys5_bat,
                l,
                Deterministic("max_active_power", forecast_data),
            )
        end
        for (ix, r) in enumerate(get_components(RenewableGen, c_sys5_bat))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(ren_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = ren_timeseries_DA[t][ix]
            end
            add_time_series!(
                c_sys5_bat,
                r,
                Deterministic("max_active_power", forecast_data),
            )
        end
    end

    if get(kwargs, :add_reserves, false)
        reserve_bat = reserve5_re(get_components(RenewableDispatch, c_sys5_bat))
        add_service!(c_sys5_bat, reserve_bat[1], get_components(GenericBattery, c_sys5_bat))
        add_service!(c_sys5_bat, reserve_bat[2], get_components(GenericBattery, c_sys5_bat))
        # ORDC
        add_service!(c_sys5_bat, reserve_bat[3], get_components(GenericBattery, c_sys5_bat))
        for (ix, serv) in enumerate(get_components(VariableReserve, c_sys5_bat))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(Reserve_ts[t])[1]
                forecast_data[ini_time] = Reserve_ts[t]
            end
            add_time_series!(c_sys5_bat, serv, Deterministic("requirement", forecast_data))
        end
        for (ix, serv) in enumerate(get_components(ReserveDemandCurve, c_sys5_bat))
            forecast_data = SortedDict{Dates.DateTime, Vector{IS.PWL}}()
            for t in 1:2
                ini_time = TS.timestamp(ORDC_cost_ts[t])[1]
                forecast_data[ini_time] = TimeSeries.values(ORDC_cost_ts[t])
            end
            resolution = TS.timestamp(ORDC_cost_ts[1])[2] - TS.timestamp(ORDC_cost_ts[1])[1]
            set_variable_cost!(
                c_sys5_bat,
                serv,
                Deterministic("variable_cost", forecast_data, resolution),
            )
        end
    end

    return c_sys5_bat
end

function build_c_sys5_il(; kwargs...)
    nodes = nodes5()
    c_sys5_il = System(
        100.0,
        nodes,
        thermal_generators5(nodes),
        loads5(nodes),
        interruptible(nodes),
        branches5(nodes);
        time_series_in_memory = get(kwargs, :time_series_in_memory, true),
    )

    if get(kwargs, :add_forecasts, true)
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_il))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(load_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = load_timeseries_DA[t][ix]
            end
            add_time_series!(c_sys5_il, l, Deterministic("max_active_power", forecast_data))
        end
        for (ix, i) in enumerate(get_components(InterruptibleLoad, c_sys5_il))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(Iload_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = Iload_timeseries_DA[t][ix]
            end
            add_time_series!(c_sys5_il, i, Deterministic("max_active_power", forecast_data))
        end
    end

    if get(kwargs, :add_reserves, false)
        reserve_il = reserve5_il(get_components(InterruptibleLoad, c_sys5_il))
        add_service!(c_sys5_il, reserve_il[1], get_components(InterruptibleLoad, c_sys5_il))
        add_service!(
            c_sys5_il,
            reserve_il[2],
            [collect(get_components(InterruptibleLoad, c_sys5_il))[end]],
        )
        # ORDC
        add_service!(c_sys5_il, reserve_il[3], get_components(InterruptibleLoad, c_sys5_il))
        for (ix, serv) in enumerate(get_components(VariableReserve, c_sys5_il))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(Reserve_ts[ix])[1]
                forecast_data[ini_time] = Reserve_ts[t]
            end
            add_time_series!(c_sys5_il, serv, Deterministic("requirement", forecast_data))
        end
        for (ix, serv) in enumerate(get_components(ReserveDemandCurve, c_sys5_il))
            forecast_data = SortedDict{Dates.DateTime, Vector{IS.PWL}}()
            for t in 1:2
                ini_time = TS.timestamp(ORDC_cost_ts[t])[1]
                forecast_data[ini_time] = TimeSeries.values(ORDC_cost_ts[t])
            end
            resolution = TS.timestamp(ORDC_cost_ts[1])[2] - TS.timestamp(ORDC_cost_ts[1])[1]
            set_variable_cost!(
                c_sys5_il,
                serv,
                Deterministic("variable_cost", forecast_data, resolution),
            )
        end
    end

    return c_sys5_il
end

function build_c_sys5_dc(; kwargs...)
    nodes = nodes5()
    c_sys5_dc = System(
        100.0,
        nodes,
        thermal_generators5(nodes),
        renewable_generators5(nodes),
        loads5(nodes),
        branches5_dc(nodes);
        time_series_in_memory = get(kwargs, :time_series_in_memory, true),
    )

    if get(kwargs, :add_forecasts, true)
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_dc))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(load_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = load_timeseries_DA[t][ix]
            end
            add_time_series!(c_sys5_dc, l, Deterministic("max_active_power", forecast_data))
        end
        for (ix, r) in enumerate(get_components(RenewableGen, c_sys5_dc))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(ren_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = ren_timeseries_DA[t][ix]
            end
            add_time_series!(c_sys5_dc, r, Deterministic("max_active_power", forecast_data))
        end
    end

    return c_sys5_dc
end

function build_c_sys14_dc(; kwargs...)
    nodes = nodes14()
    c_sys14_dc = System(
        100.0,
        nodes,
        thermal_generators14(nodes),
        loads14(nodes),
        branches14_dc(nodes);
        time_series_in_memory = get(kwargs, :time_series_in_memory, true),
    )

    if get(kwargs, :add_forecasts, true)
        forecast_data = SortedDict{Dates.DateTime, TimeArray}()
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys14_dc))
            ini_time = TS.timestamp(timeseries_DA14[ix])[1]
            forecast_data[ini_time] = timeseries_DA14[ix]
            add_time_series!(
                c_sys14_dc,
                l,
                Deterministic("max_active_power", forecast_data),
            )
        end
    end

    return c_sys14_dc
end

function build_c_sys5_reg(; kwargs...)
    nodes = nodes5()

    c_sys5_reg =
        System(100.0, nodes, thermal_generators5(nodes), loads5(nodes), branches5(nodes))

    area = Area("1")
    add_component!(c_sys5_reg, area)
    [set_area!(b, area) for b in get_components(Bus, c_sys5_reg)]
    AGC_service = PSY.AGC(
        name = "AGC_Area1",
        available = true,
        bias = 739.0,
        K_p = 2.5,
        K_i = 0.1,
        K_d = 0.0,
        delta_t = 4,
        area = first(get_components(Area, c_sys5_reg)),
    )
    #add_component!(c_sys5_reg, AGC_service)
    if get(kwargs, :add_forecasts, true)
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_reg))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(load_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = load_timeseries_DA[t][ix]
            end
            add_time_series!(
                c_sys5_reg,
                l,
                Deterministic("max_active_power", forecast_data),
            )
        end
        for (_, l) in enumerate(get_components(ThermalStandard, c_sys5_reg))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(load_timeseries_DA[t][1])[1]
                forecast_data[ini_time] = load_timeseries_DA[t][1]
            end
            add_time_series!(
                c_sys5_reg,
                l,
                Deterministic("max_active_power", forecast_data),
            )
        end
    end

    contributing_devices = Vector()
    for g in get_components(Generator, c_sys5_reg)
        droop =
            isa(g, ThermalStandard) ? 0.04 * get_base_power(g) : 0.05 * get_base_power(g)
        p_factor = (up = 1.0, dn = 1.0)
        t = RegulationDevice(g, participation_factor = p_factor, droop = droop)
        add_component!(c_sys5_reg, t)
        push!(contributing_devices, t)
    end
    add_service!(c_sys5_reg, AGC_service, contributing_devices)
    return c_sys5_reg
end

# System to test UC Forms
#Park City and Sundance Have non-binding Ramp Limitst at an Hourly Resolution
# Solitude, Sundance and Brighton have binding time_up constraints.
# Solitude and Brighton have binding time_dn constraints.
# Sundance has non-binding Time Down constraint at an Hourly Resolution
# Alta, Park City and Brighton start at 0.

thermal_generators5_uc_testing(nodes) = [
    ThermalStandard(
        name = "Alta",
        available = true,
        status = false,
        bus = nodes[1],
        active_power = 0.0,
        reactive_power = 0.0,
        rating = 0.5,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 0.2, max = 0.40),
        reactive_power_limits = (min = -0.30, max = 0.30),
        ramp_limits = nothing,
        time_limits = nothing,
        operation_cost = ThreePartCost((0.0, 14.0), 0.0, 4.0, 2.0),
        base_power = 100.0,
    ),
    ThermalStandard(
        name = "Park City",
        available = true,
        status = false,
        bus = nodes[1],
        active_power = 0.0,
        reactive_power = 0.0,
        rating = 2.2125,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 0.65, max = 1.70),
        reactive_power_limits = (min = -1.275, max = 1.275),
        ramp_limits = (up = 0.02 * 2.2125, down = 0.02 * 2.2125),
        time_limits = nothing,
        operation_cost = ThreePartCost((0.0, 15.0), 0.0, 1.5, 0.75),
        base_power = 100.0,
    ),
    ThermalStandard(
        name = "Solitude",
        available = true,
        status = true,
        bus = nodes[3],
        active_power = 2.7,
        reactive_power = 0.00,
        rating = 5.20,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 1.0, max = 5.20),
        reactive_power_limits = (min = -3.90, max = 3.90),
        ramp_limits = (up = 0.0012 * 5.2, down = 0.0012 * 5.2),
        time_limits = (up = 5.0, down = 3.0),
        operation_cost = ThreePartCost((0.0, 30.0), 0.0, 3.0, 1.5),
        base_power = 100.0,
    ),
    ThermalStandard(
        name = "Sundance",
        available = true,
        status = false,
        bus = nodes[4],
        active_power = 0.0,
        reactive_power = 0.00,
        rating = 2.5,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 1.0, max = 2.0),
        reactive_power_limits = (min = -1.5, max = 1.5),
        ramp_limits = (up = 0.015 * 2.5, down = 0.015 * 2.5),
        time_limits = (up = 2.0, down = 1.0),
        operation_cost = ThreePartCost((0.0, 40.0), 0.0, 4.0, 2.0),
        base_power = 100.0,
    ),
    ThermalStandard(
        name = "Brighton",
        available = true,
        status = true,
        bus = nodes[5],
        active_power = 6.0,
        reactive_power = 0.0,
        rating = 7.5,
        prime_mover = PrimeMovers.ST,
        fuel = ThermalFuels.COAL,
        active_power_limits = (min = 3.0, max = 6.0),
        reactive_power_limits = (min = -4.50, max = 4.50),
        ramp_limits = (up = 0.0015 * 7.5, down = 0.0015 * 7.5),
        time_limits = (up = 5.0, down = 3.0),
        operation_cost = ThreePartCost((0.0, 10.0), 0.0, 1.5, 0.75),
        base_power = 100.0,
    ),
];

function build_sys_ramp_testing(; kwargs...)
    node = Bus(1, "nodeA", "REF", 0, 1.0, (min = 0.9, max = 1.05), 230, nothing, nothing)
    load = PowerLoad("Bus1", true, node, nothing, 0.4, 0.9861, 100.0, 1.0, 2.0)
    gen_ramp = [
        ThermalStandard(
            name = "Alta",
            available = true,
            status = true,
            bus = node,
            active_power = 0.20, # Active power
            reactive_power = 0.010,
            rating = 0.5,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.0, max = 0.40),
            reactive_power_limits = nothing,
            ramp_limits = nothing,
            time_limits = nothing,
            operation_cost = ThreePartCost((0.0, 14.0), 0.0, 4.0, 2.0),
            base_power = 100.0,
        ),
        ThermalStandard(
            name = "Park City",
            available = true,
            status = true,
            bus = node,
            active_power = 0.70, # Active Power
            reactive_power = 0.20,
            rating = 2.0,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.7, max = 2.20),
            reactive_power_limits = nothing,
            ramp_limits = (up = 0.010625 * 2.0, down = 0.010625 * 2.0),
            time_limits = nothing,
            operation_cost = ThreePartCost((0.0, 15.0), 0.0, 1.5, 0.75),
            base_power = 100.0,
        ),
    ]
    DA_ramp = collect(
        DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):Hour(1):DateTime(
            "1/1/2024  4:00:00",
            "d/m/y  H:M:S",
        ),
    )
    ramp_load = [0.9, 1.1, 2.485, 2.175, 0.9]
    ts_dict = SortedDict(DA_ramp[1] => ramp_load)
    load_forecast_ramp = Deterministic("max_active_power", ts_dict, Hour(1))
    ramp_test_sys = System(100.0)
    add_component!(ramp_test_sys, node)
    add_component!(ramp_test_sys, load)
    add_component!(ramp_test_sys, gen_ramp[1])
    add_component!(ramp_test_sys, gen_ramp[2])
    add_time_series!(ramp_test_sys, load, load_forecast_ramp)
    return ramp_test_sys
end

function build_c_sys5_uc(; kwargs...)
    nodes = nodes5()
    c_sys5_uc = System(
        100.0,
        nodes,
        thermal_generators5_uc_testing(nodes),
        loads5(nodes),
        branches5(nodes);
        time_series_in_memory = get(kwargs, :time_series_in_memory, true),
    )

    if get(kwargs, :add_forecasts, true)
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_uc))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(load_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = load_timeseries_DA[t][ix]
            end
            add_time_series!(c_sys5_uc, l, Deterministic("max_active_power", forecast_data))
        end
        for (ix, r) in enumerate(get_components(RenewableGen, c_sys5_uc))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(ren_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = ren_timeseries_DA[t][ix]
            end
            add_time_series!(c_sys5_uc, r, Deterministic("max_active_power", forecast_data))
        end
        for (ix, i) in enumerate(get_components(InterruptibleLoad, c_sys5_uc))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(Iload_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = Iload_timeseries_DA[t][ix]
            end
            add_time_series!(c_sys5_uc, r, Deterministic("max_active_power", forecast_data))
        end
    end

    if get(kwargs, :add_reserves, false)
        reserve_uc = reserve5(get_components(ThermalStandard, c_sys5_uc))
        add_service!(c_sys5_uc, reserve_uc[1], get_components(ThermalStandard, c_sys5_uc))
        add_service!(
            c_sys5_uc,
            reserve_uc[2],
            [collect(get_components(ThermalStandard, c_sys5_uc))[end]],
        )
        add_service!(c_sys5_uc, reserve_uc[3], get_components(ThermalStandard, c_sys5_uc))
        # ORDC Curve
        add_service!(c_sys5_uc, reserve_uc[4], get_components(ThermalStandard, c_sys5_uc))
        for serv in get_components(VariableReserve, c_sys5_uc)
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(Reserve_ts[t])[1]
                forecast_data[ini_time] = Reserve_ts[t]
            end
            add_time_series!(c_sys5_uc, serv, Deterministic("requirement", forecast_data))
        end
        for (ix, serv) in enumerate(get_components(ReserveDemandCurve, c_sys5_uc))
            forecast_data = SortedDict{Dates.DateTime, Vector{IS.PWL}}()
            for t in 1:2
                ini_time = TS.timestamp(ORDC_cost_ts[t])[1]
                forecast_data[ini_time] = TimeSeries.values(ORDC_cost_ts[t])
            end
            resolution = TS.timestamp(ORDC_cost_ts[1])[2] - TS.timestamp(ORDC_cost_ts[1])[1]
            set_variable_cost!(
                c_sys5_uc,
                serv,
                Deterministic("variable_cost", forecast_data, resolution),
            )
        end
    end

    return c_sys5_uc
end

function build_c_sys5_pwl_uc(; kwargs...)
    c_sys5_uc = build_c_sys5_uc(; kwargs...)
    thermal = thermal_generators5_pwl(collect(get_components(Bus, c_sys5_uc)))
    for d in thermal
        PSY.add_component!(c_sys5_uc, d)
    end
    return c_sys5_uc
end

function build_c_sys5_ed(; kwargs...)
    nodes = nodes5()
    c_sys5_ed = System(
        100.0,
        nodes,
        thermal_generators5_uc_testing(nodes),
        renewable_generators5(nodes),
        loads5(nodes),
        interruptible(nodes),
        branches5(nodes);
        time_series_in_memory = get(kwargs, :time_series_in_memory, true),
    )

    if get(kwargs, :add_forecasts, true)
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_ed))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2 # loop over days
                ta = load_timeseries_DA[t][ix]
                for i in 1:length(ta) # loop over hours
                    ini_time = TS.timestamp(ta[i]) #get the hour
                    data = when(load_timeseries_RT[t][ix], hour, hour(ini_time[1])) # get the subset ts for that hour
                    forecast_data[ini_time[1]] = data
                end
            end
            add_time_series!(c_sys5_ed, l, Deterministic("max_active_power", forecast_data))
        end
        for (ix, l) in enumerate(get_components(RenewableGen, c_sys5_ed))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2 # loop over days
                ta = load_timeseries_DA[t][ix]
                for i in 1:length(ta) # loop over hours
                    ini_time = TS.timestamp(ta[i]) #get the hour
                    data = when(load_timeseries_RT[t][ix], hour, hour(ini_time[1])) # get the subset ts for that hour
                    forecast_data[ini_time[1]] = data
                end
            end
            add_time_series!(c_sys5_ed, l, Deterministic("max_active_power", forecast_data))
        end
        for (ix, l) in enumerate(get_components(InterruptibleLoad, c_sys5_ed))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2 # loop over days
                ta = load_timeseries_DA[t][ix]
                for i in 1:length(ta) # loop over hours
                    ini_time = TS.timestamp(ta[i]) #get the hour
                    data = when(load_timeseries_RT[t][ix], hour, hour(ini_time[1])) # get the subset ts for that hour
                    forecast_data[ini_time[1]] = data
                end
            end
            add_time_series!(c_sys5_ed, l, Deterministic("max_active_power", forecast_data))
        end
    end

    return c_sys5_ed
end

function build_c_sys5_pwl_ed(; kwargs...)
    c_sys5_ed = build_c_sys5_ed(; kwargs...)
    thermal = thermal_generators5_pwl(collect(get_components(Bus, c_sys5_ed)))
    for d in thermal
        PSY.add_component!(c_sys5_ed, d)
    end
    return c_sys5_ed
end

function build_c_sys5_pwl_ed_nonconvex(; kwargs...)
    c_sys5_ed = build_c_sys5_ed(; kwargs...)
    thermal = thermal_generators5_pwl_nonconvex(collect(get_components(Bus, c_sys5_ed)))
    for d in thermal
        PSY.add_component!(c_sys5_ed, d)
    end
    return c_sys5_ed
end

function build_init(gens, data)
    init = Vector{InitialCondition}(undef, length(collect(gens)))
    for (ix, g) in enumerate(gens)
        init[ix] = InitialCondition(
            g,
            PSI.UpdateRef{JuMP.VariableRef}(PSI.ACTIVE_POWER),
            data[ix],
            TimeStatusChange,
        )
    end
    return init
end

function build_c_sys5_hy_uc(; kwargs...)
    nodes = nodes5()
    c_sys5_hy_uc = System(
        100.0,
        nodes,
        thermal_generators5_uc_testing(nodes),
        hydro_generators5(nodes),
        renewable_generators5(nodes),
        loads5(nodes),
        branches5(nodes);
        time_series_in_memory = get(kwargs, :time_series_in_memory, true),
    )

    if get(kwargs, :add_forecasts, true)
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_hy_uc))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(load_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = load_timeseries_DA[t][ix]
            end
            add_time_series!(
                c_sys5_hy_uc,
                l,
                Deterministic("max_active_power", forecast_data),
            )
        end
        for (ix, h) in enumerate(get_components(HydroEnergyReservoir, c_sys5_hy_uc))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(hydro_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = hydro_timeseries_DA[t][ix]
            end
            add_time_series!(
                c_sys5_hy_uc,
                h,
                Deterministic("max_active_power", forecast_data),
            )
        end
        for (ix, h) in enumerate(get_components(HydroEnergyReservoir, c_sys5_hy_uc))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(hydro_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = hydro_timeseries_DA[t][ix]
            end
            add_time_series!(
                c_sys5_hy_uc,
                h,
                Deterministic("storage_capacity", forecast_data),
            )
            add_time_series!(
                c_sys5_hy_uc,
                h,
                Deterministic("storage_target", forecast_data),
            )
        end
        for (ix, h) in enumerate(get_components(HydroEnergyReservoir, c_sys5_hy_uc))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(hydro_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = hydro_timeseries_DA[t][ix] .* 0.8
            end
            add_time_series!(c_sys5_hy_uc, h, Deterministic("inflow", forecast_data))
        end
        for (ix, h) in enumerate(get_components(HydroDispatch, c_sys5_hy_uc))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(hydro_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = hydro_timeseries_DA[t][ix]
            end
            add_time_series!(
                c_sys5_hy_uc,
                h,
                Deterministic("max_active_power", forecast_data),
            )
        end
        for (ix, r) in enumerate(get_components(RenewableGen, c_sys5_hy_uc))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(ren_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = ren_timeseries_DA[t][ix]
            end
            add_time_series!(
                c_sys5_hy_uc,
                r,
                Deterministic("max_active_power", forecast_data),
            )
        end
        for (ix, i) in enumerate(get_components(InterruptibleLoad, c_sys5_hy_uc))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(Iload_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = Iload_timeseries_DA[t][ix]
            end
            add_time_series!(
                c_sys5_hy_uc,
                i,
                Deterministic("max_active_power", forecast_data),
            )
        end
    end

    if get(kwargs, :add_reserves, false)
        reserve_uc = reserve5(get_components(ThermalStandard, c_sys5_hy_uc))
        add_service!(
            c_sys5_hy_uc,
            reserve_uc[1],
            get_components(ThermalStandard, c_sys5_hy_uc),
        )
        add_service!(
            c_sys5_hy_uc,
            reserve_uc[2],
            [collect(get_components(ThermalStandard, c_sys5_hy_uc))[end]],
        )
        add_service!(
            c_sys5_hy_uc,
            reserve_uc[3],
            get_components(ThermalStandard, c_sys5_hy_uc),
        )
        # ORDC Curve
        add_service!(
            c_sys5_hy_uc,
            reserve_uc[4],
            get_components(ThermalStandard, c_sys5_hy_uc),
        )
        for serv in get_components(VariableReserve, c_sys5_hy_uc)
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(Reserve_ts[t])[1]
                forecast_data[ini_time] = Reserve_ts[t]
            end
            add_time_series!(
                c_sys5_hy_uc,
                serv,
                Deterministic("requirement", forecast_data),
            )
        end
        for (ix, serv) in enumerate(get_components(ReserveDemandCurve, c_sys5_hy_uc))
            forecast_data = SortedDict{Dates.DateTime, Vector{IS.PWL}}()
            for t in 1:2
                ini_time = TS.timestamp(ORDC_cost_ts[t])[1]
                forecast_data[ini_time] = TimeSeries.values(ORDC_cost_ts[t])
            end
            resolution = TS.timestamp(ORDC_cost_ts[1])[2] - TS.timestamp(ORDC_cost_ts[1])[1]
            set_variable_cost!(
                c_sys5_hy_uc,
                serv,
                Deterministic("variable_cost", forecast_data, resolution),
            )
        end
    end

    return c_sys5_hy_uc
end

function build_c_sys5_hy_ed(; kwargs...)
    nodes = nodes5()
    c_sys5_hy_ed = System(
        100.0,
        nodes,
        thermal_generators5_uc_testing(nodes),
        hydro_generators5(nodes),
        renewable_generators5(nodes),
        loads5(nodes),
        interruptible(nodes),
        branches5(nodes);
        time_series_in_memory = get(kwargs, :time_series_in_memory, true),
    )

    if get(kwargs, :add_forecasts, true)
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_hy_ed))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2 # loop over days
                ta = load_timeseries_DA[t][ix]
                for i in 1:length(ta) # loop over hours
                    ini_time = TS.timestamp(ta[i]) #get the hour
                    data = when(load_timeseries_RT[t][ix], hour, hour(ini_time[1])) # get the subset ts for that hour
                    forecast_data[ini_time[1]] = data
                end
            end
            add_time_series!(
                c_sys5_hy_ed,
                l,
                Deterministic("max_active_power", forecast_data),
            )
        end
        for (ix, l) in enumerate(get_components(HydroEnergyReservoir, c_sys5_hy_ed))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ta = hydro_timeseries_DA[t][ix]
                for i in 1:length(ta)
                    ini_time = TS.timestamp(ta[i])
                    data = when(hydro_timeseries_RT[t][ix], hour, hour(ini_time[1]))
                    forecast_data[ini_time[1]] = data
                end
            end
            add_time_series!(
                c_sys5_hy_ed,
                l,
                Deterministic("max_active_power", forecast_data),
            )
        end
        for (ix, l) in enumerate(get_components(RenewableGen, c_sys5_hy_ed))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ta = load_timeseries_DA[t][ix]
                for i in 1:length(ta)
                    ini_time = TS.timestamp(ta[i])
                    data = when(load_timeseries_RT[t][ix], hour, hour(ini_time[1]))
                    forecast_data[ini_time[1]] = data
                end
            end
            add_time_series!(
                c_sys5_hy_ed,
                l,
                Deterministic("max_active_power", forecast_data),
            )
        end
        for (ix, l) in enumerate(get_components(HydroEnergyReservoir, c_sys5_hy_ed))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ta = hydro_timeseries_DA[t][ix]
                for i in 1:length(ta)
                    ini_time = TS.timestamp(ta[i])
                    data = when(hydro_timeseries_RT[t][ix], hour, hour(ini_time[1]))
                    forecast_data[ini_time[1]] = data
                end
            end
            add_time_series!(
                c_sys5_hy_ed,
                l,
                Deterministic("storage_capacity", forecast_data),
            )
            add_time_series!(
                c_sys5_hy_ed,
                l,
                Deterministic("storage_target", forecast_data),
            )
        end
        for (ix, l) in enumerate(get_components(HydroEnergyReservoir, c_sys5_hy_ed))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ta = hydro_timeseries_DA[t][ix]
                for i in 1:length(ta)
                    ini_time = TS.timestamp(ta[i])
                    data = when(hydro_timeseries_RT[t][ix] .* 0.8, hour, hour(ini_time[1]))
                    forecast_data[ini_time[1]] = data
                end
            end
            add_time_series!(c_sys5_hy_ed, l, Deterministic("inflow", forecast_data))
        end
        for (ix, l) in enumerate(get_components(InterruptibleLoad, c_sys5_hy_ed))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ta = load_timeseries_DA[t][ix]
                for i in 1:length(ta)
                    ini_time = TS.timestamp(ta[i])
                    data = when(load_timeseries_RT[t][ix], hour, hour(ini_time[1]))
                    forecast_data[ini_time[1]] = data
                end
            end
            add_time_series!(
                c_sys5_hy_ed,
                l,
                Deterministic("max_active_power", forecast_data),
            )
        end
        for (ix, l) in enumerate(get_components(HydroDispatch, c_sys5_hy_ed))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ta = hydro_timeseries_DA[t][ix]
                for i in 1:length(ta)
                    ini_time = TS.timestamp(ta[i])
                    data = when(hydro_timeseries_RT[t][ix], hour, hour(ini_time[1]))
                    forecast_data[ini_time[1]] = data
                end
            end
            add_time_series!(
                c_sys5_hy_ed,
                l,
                Deterministic("max_active_power", forecast_data),
            )
        end
    end

    return c_sys5_hy_ed
end

function build_c_sys5_phes_ed(; kwargs...)
    nodes = nodes5()
    c_sys5_phes_ed = System(
        100.0,
        nodes,
        thermal_generators5_uc_testing(nodes),
        phes5(nodes),
        renewable_generators5(nodes),
        loads5(nodes),
        interruptible(nodes),
        branches5(nodes);
        time_series_in_memory = get(kwargs, :time_series_in_memory, true),
    )

    if get(kwargs, :add_forecasts, true)
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_phes_ed))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2 # loop over days
                ta = load_timeseries_DA[t][ix]
                for i in 1:length(ta) # loop over hours
                    ini_time = TS.timestamp(ta[i]) #get the hour
                    data = when(load_timeseries_RT[t][ix], hour, hour(ini_time[1])) # get the subset ts for that hour
                    forecast_data[ini_time[1]] = data
                end
            end
            add_time_series!(
                c_sys5_phes_ed,
                l,
                Deterministic("max_active_power", forecast_data),
            )
        end
        for (ix, l) in enumerate(get_components(HydroGen, c_sys5_phes_ed))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ta = hydro_timeseries_DA[t][ix]
                for i in 1:length(ta)
                    ini_time = TS.timestamp(ta[i])
                    data = when(hydro_timeseries_RT[t][ix], hour, hour(ini_time[1]))
                    forecast_data[ini_time[1]] = data
                end
            end
            add_time_series!(
                c_sys5_phes_ed,
                l,
                Deterministic("max_active_power", forecast_data),
            )
        end
        for (ix, l) in enumerate(get_components(RenewableGen, c_sys5_phes_ed))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ta = load_timeseries_DA[t][ix]
                for i in 1:length(ta)
                    ini_time = TS.timestamp(ta[i])
                    data = when(load_timeseries_RT[t][ix], hour, hour(ini_time[1]))
                    forecast_data[ini_time[1]] = data
                end
            end
            add_time_series!(
                c_sys5_phes_ed,
                l,
                Deterministic("max_active_power", forecast_data),
            )
        end
        for (ix, l) in enumerate(get_components(HydroPumpedStorage, c_sys5_phes_ed))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ta = hydro_timeseries_DA[t][ix]
                for i in 1:length(ta)
                    ini_time = TS.timestamp(ta[i])
                    data = when(hydro_timeseries_RT[t][ix], hour, hour(ini_time[1]))
                    forecast_data[ini_time[1]] = data
                end
            end
            add_time_series!(
                c_sys5_phes_ed,
                l,
                Deterministic("storage_capacity", forecast_data),
            )
        end
        for (ix, l) in enumerate(get_components(HydroPumpedStorage, c_sys5_phes_ed))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ta = hydro_timeseries_DA[t][ix]
                for i in 1:length(ta)
                    ini_time = TS.timestamp(ta[i])
                    data = when(hydro_timeseries_RT[t][ix] .* 0.8, hour, hour(ini_time[1]))
                    forecast_data[ini_time[1]] = data
                end
            end
            add_time_series!(c_sys5_phes_ed, l, Deterministic("inflow", forecast_data))
            add_time_series!(c_sys5_phes_ed, l, Deterministic("outflow", forecast_data))
        end
        for (ix, l) in enumerate(get_components(InterruptibleLoad, c_sys5_phes_ed))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ta = load_timeseries_DA[t][ix]
                for i in 1:length(ta)
                    ini_time = TS.timestamp(ta[i])
                    data = when(load_timeseries_RT[t][ix], hour, hour(ini_time[1]))
                    forecast_data[ini_time[1]] = data
                end
            end
            add_time_series!(
                c_sys5_phes_ed,
                l,
                Deterministic("max_active_power", forecast_data),
            )
        end
    end

    return c_sys5_phes_ed
end

function build_c_sys5_pglib(; kwargs...)
    nodes = nodes5()
    c_sys5_uc = System(
        100.0,
        nodes,
        thermal_generators5_uc_testing(nodes),
        thermal_pglib_generators5(nodes),
        loads5(nodes),
        branches5(nodes);
        time_series_in_memory = get(kwargs, :time_series_in_memory, true),
    )

    if get(kwargs, :add_forecasts, true)
        forecast_data = SortedDict{Dates.DateTime, TimeArray}()
        for (ix, l) in enumerate(get_components(PowerLoad, c_sys5_uc))
            for t in 1:2
                ini_time = TS.timestamp(load_timeseries_DA[t][ix])[1]
                forecast_data[ini_time] = load_timeseries_DA[t][ix]
            end
            add_time_series!(c_sys5_uc, l, Deterministic("max_active_power", forecast_data))
        end
    end

    if get(kwargs, :add_reserves, false)
        reserve_uc = reserve5(get_components(ThermalStandard, c_sys5_uc))
        add_service!(c_sys5_uc, reserve_uc[1], get_components(ThermalStandard, c_sys5_uc))
        add_service!(
            c_sys5_uc,
            reserve_uc[2],
            [collect(get_components(ThermalStandard, c_sys5_uc))[end]],
        )
        add_service!(c_sys5_uc, reserve_uc[3], get_components(ThermalStandard, c_sys5_uc))
        for (ix, serv) in enumerate(get_components(VariableReserve, c_sys5_uc))
            forecast_data = SortedDict{Dates.DateTime, TimeArray}()
            for t in 1:2
                ini_time = TS.timestamp(Reserve_ts[ix])[1]
                forecast_data[ini_time] = Reserve_ts[t]
            end
            add_time_series!(c_sys5_uc, serv, Deterministic("requirement", forecast_data))
        end
    end

    return c_sys5_uc
end

function build_sos_test_sys(; kwargs...)
    node = Bus(1, "nodeA", "PV", 0, 1.0, (min = 0.9, max = 1.05), 230, nothing, nothing)
    load = PowerLoad("Bus1", true, node, nothing, 0.4, 0.9861, 100.0, 1.0, 2.0)
    gens_cost_sos = [
        ThermalStandard(
            name = "Alta",
            available = true,
            status = true,
            bus = node,
            active_power = 0.52,
            reactive_power = 0.010,
            rating = 0.5,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.22, max = 0.55),
            reactive_power_limits = nothing,
            time_limits = nothing,
            ramp_limits = nothing,
            operation_cost = ThreePartCost(
                [(1122.43, 22.0), (1617.43, 33.0), (1742.48, 44.0), (2075.88, 55.0)],
                0.0,
                5665.23,
                0.0,
            ),
            base_power = 100.0,
        ),
        ThermalStandard(
            name = "Park City",
            available = true,
            status = true,
            bus = node,
            active_power = 0.62,
            reactive_power = 0.20,
            rating = 2.2125,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.62, max = 1.55),
            reactive_power_limits = nothing,
            time_limits = nothing,
            ramp_limits = nothing,
            operation_cost = ThreePartCost(
                [(1500.19, 62.0), (2132.59, 92.9), (2829.875, 124.0), (2831.444, 155.0)],
                0.0,
                5665.23,
                0.0,
            ),
            base_power = 100.0,
        ),
    ]

    #Checks the data remains non-convex
    for g in gens_cost_sos
        @assert PSI.pwlparamcheck(PSY.get_operation_cost(g).variable) == false
    end

    DA_load_forecast = SortedDict{Dates.DateTime, TimeArray}()
    ini_time = DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S")
    cost_sos_load = [[1.3, 2.1], [1.3, 2.1]]
    for (ix, date) in enumerate(range(ini_time; length = 2, step = Hour(1)))
        DA_load_forecast[date] = TimeArray([date, date + Hour(1)], cost_sos_load[ix])
    end
    load_forecast_cost_sos = Deterministic("max_active_power", DA_load_forecast)
    cost_test_sos_sys =
        System(100.0; time_series_in_memory = get(kwargs, :time_series_in_memory, true))
    add_component!(cost_test_sos_sys, node)
    add_component!(cost_test_sos_sys, load)
    add_component!(cost_test_sos_sys, gens_cost_sos[1])
    add_component!(cost_test_sos_sys, gens_cost_sos[2])
    add_time_series!(cost_test_sos_sys, load, load_forecast_cost_sos)

    return cost_test_sos_sys
end

function build_pwl_test_sys(; kwargs...)
    node = Bus(1, "nodeA", "PV", 0, 1.0, (min = 0.9, max = 1.05), 230, nothing, nothing)
    load = PowerLoad("Bus1", true, node, nothing, 0.4, 0.9861, 100.0, 1.0, 2.0)
    gens_cost = [
        ThermalStandard(
            name = "Alta",
            available = true,
            status = true,
            bus = node,
            active_power = 0.52,
            reactive_power = 0.010,
            rating = 0.5,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.22, max = 0.55),
            reactive_power_limits = nothing,
            time_limits = nothing,
            ramp_limits = nothing,
            operation_cost = ThreePartCost(
                [(589.99, 22.0), (884.99, 33.0), (1210.04, 44.0), (1543.44, 55.0)],
                532.44,
                5665.23,
                0.0,
            ),
            base_power = 100.0,
        ),
        ThermalStandard(
            name = "Park City",
            available = true,
            status = true,
            bus = node,
            active_power = 0.62,
            reactive_power = 0.20,
            rating = 221.25,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.62, max = 1.55),
            reactive_power_limits = nothing,
            time_limits = nothing,
            ramp_limits = nothing,
            operation_cost = ThreePartCost(
                [(1264.80, 62.0), (1897.20, 93.0), (2594.4787, 124.0), (3433.04, 155.0)],
                235.397,
                5665.23,
                0.0,
            ),
            base_power = 100.0,
        ),
    ]

    DA_load_forecast = SortedDict{Dates.DateTime, TimeArray}()
    ini_time = DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S")
    cost_sos_load = [[1.3, 2.1], [1.3, 2.1]]
    for (ix, date) in enumerate(range(ini_time; length = 2, step = Hour(1)))
        DA_load_forecast[date] = TimeArray([date, date + Hour(1)], cost_sos_load[ix])
    end
    load_forecast_cost_sos = Deterministic("max_active_power", DA_load_forecast)
    cost_test_sys =
        System(100.0; time_series_in_memory = get(kwargs, :time_series_in_memory, true))
    add_component!(cost_test_sys, node)
    add_component!(cost_test_sys, load)
    add_component!(cost_test_sys, gens_cost[1])
    add_component!(cost_test_sys, gens_cost[2])
    add_time_series!(cost_test_sys, load, load_forecast_cost_sos)
    return cost_test_sys
end

function build_duration_test_sys(; kwargs...)
    node = Bus(1, "nodeA", "PV", 0, 1.0, (min = 0.9, max = 1.05), 230, nothing, nothing)
    load = PowerLoad("Bus1", true, node, nothing, 0.4, 0.9861, 100.0, 1.0, 2.0)
    DA_dur = collect(
        DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):Hour(1):DateTime(
            "1/1/2024  6:00:00",
            "d/m/y  H:M:S",
        ),
    )
    gens_dur = [
        ThermalStandard(
            name = "Alta",
            available = true,
            status = true,
            bus = node,
            active_power = 0.40,
            reactive_power = 0.010,
            rating = 0.5,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.3, max = 0.9),
            reactive_power_limits = nothing,
            ramp_limits = nothing,
            time_limits = (up = 4, down = 2),
            operation_cost = ThreePartCost((0.0, 14.0), 0.0, 4.0, 2.0),
            base_power = 100.0,
            time_at_status = 2.0,
        ),
        ThermalStandard(
            name = "Park City",
            available = true,
            status = false,
            bus = node,
            active_power = 1.70,
            reactive_power = 0.20,
            rating = 2.2125,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.7, max = 2.2),
            reactive_power_limits = nothing,
            ramp_limits = nothing,
            time_limits = (up = 6, down = 4),
            operation_cost = ThreePartCost((0.0, 15.0), 0.0, 1.5, 0.75),
            base_power = 100.0,
            time_at_status = 3.0,
        ),
    ]

    duration_load = [0.3, 0.6, 0.8, 0.7, 1.7, 0.9, 0.7]
    load_data = SortedDict(DA_dur[1] => TimeArray(DA_dur, duration_load))
    load_forecast_dur = Deterministic("max_active_power", load_data)
    duration_test_sys =
        System(100.0; time_series_in_memory = get(kwargs, :time_series_in_memory, true))
    add_component!(duration_test_sys, node)
    add_component!(duration_test_sys, load)
    add_component!(duration_test_sys, gens_dur[1])
    add_component!(duration_test_sys, gens_dur[2])
    add_time_series!(duration_test_sys, load, load_forecast_dur)

    return duration_test_sys
end

function build_pwl_marketbid_sys(; kwargs...)
    node = Bus(1, "nodeA", "PV", 0, 1.0, (min = 0.9, max = 1.05), 230, nothing, nothing)
    load = PowerLoad("Bus1", true, node, nothing, 0.4, 0.9861, 100.0, 1.0, 2.0)
    gens_cost = [
        ThermalStandard(
            name = "Alta",
            available = true,
            status = true,
            bus = node,
            active_power = 0.52,
            reactive_power = 0.010,
            rating = 0.5,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.22, max = 0.55),
            reactive_power_limits = nothing,
            time_limits = nothing,
            ramp_limits = nothing,
            operation_cost = MarketBidCost(
                no_load = 0.0,
                start_up = (hot = 0.0, warm = 0.0, cold = 0.0),
                shut_down = 0.0,
            ),
            base_power = 100.0,
        ),
        ThermalMultiStart(
            name = "115_STEAM_1",
            available = true,
            status = true,
            bus = node,
            active_power = 0.05,
            reactive_power = 0.010,
            rating = 0.12,
            prime_mover = PrimeMovers.ST,
            fuel = ThermalFuels.COAL,
            active_power_limits = (min = 0.05, max = 0.12),
            reactive_power_limits = (min = -0.30, max = 0.30),
            ramp_limits = (up = 0.2 * 0.12, down = 0.2 * 0.12),
            power_trajectory = (startup = 0.05, shutdown = 0.05),
            time_limits = (up = 4.0, down = 2.0),
            start_time_limits = (hot = 2.0, warm = 4.0, cold = 12.0),
            start_types = 3,
            operation_cost = MarketBidCost(
                no_load = 0.0,
                start_up = (hot = 393.28, warm = 455.37, cold = 703.76),
                shut_down = 0.0,
            ),
            base_power = 100.0,
        ),
    ]
    ini_time = DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S")
    DA_load_forecast = Dict{Dates.DateTime, TimeArray}()
    market_bid_gen1_data = Dict(
        ini_time => [
            [(589.99, 22.0), (884.99, 33.0), (1210.04, 44.0), (1543.44, 55.0)],
            [(589.99, 22.0), (884.99, 33.0), (1210.04, 44.0), (1543.44, 55.0)],
        ],
        ini_time + Hour(1) => [
            [(589.99, 22.0), (884.99, 33.0), (1210.04, 44.0), (1543.44, 55.0)],
            [(589.99, 22.0), (884.99, 33.0), (1210.04, 44.0), (1543.44, 55.0)],
        ],
    )
    market_bid_gen1 = Deterministic(
        name = "variable_cost",
        data = market_bid_gen1_data,
        resolution = Hour(1),
    )
    market_bid_gen2_data = Dict(
        ini_time => [
            [(0.0, 5.0), (290.1, 7.33), (582.72, 9.67), (894.1, 12.0)],
            [(0.0, 5.0), (300.1, 7.33), (600.72, 9.67), (900.1, 12.0)],
        ],
        ini_time + Hour(1) => [
            [(0.0, 5.0), (290.1, 7.33), (582.72, 9.67), (894.1, 12.0)],
            [(0.0, 5.0), (300.1, 7.33), (600.72, 9.67), (900.1, 12.0)],
        ],
    )
    market_bid_gen2 = Deterministic(
        name = "variable_cost",
        data = market_bid_gen2_data,
        resolution = Hour(1),
    )
    market_bid_load = [[1.3, 2.1], [1.3, 2.1]]
    for (ix, date) in enumerate(range(ini_time; length = 2, step = Hour(1)))
        DA_load_forecast[date] = TimeArray([date, date + Hour(1)], market_bid_load[ix])
    end
    load_forecast_cost_market_bid = Deterministic("max_active_power", DA_load_forecast)
    cost_test_sys =
        System(100.0; time_series_in_memory = get(kwargs, :time_series_in_memory, true))
    add_component!(cost_test_sys, node)
    add_component!(cost_test_sys, load)
    add_component!(cost_test_sys, gens_cost[1])
    add_component!(cost_test_sys, gens_cost[2])
    add_time_series!(cost_test_sys, load, load_forecast_cost_market_bid)
    set_variable_cost!(cost_test_sys, gens_cost[1], market_bid_gen1)
    set_variable_cost!(cost_test_sys, gens_cost[2], market_bid_gen2)
    return cost_test_sys
end

TEST_SYSTEMS = Dict(
    "c_sys14" => (
        description = "14-bus system",
        build = build_c_sys14,
        time_series_in_memory = true,
    ),
    "c_sys14_dc" => (
        description = "14-bus system with DC line",
        build = build_c_sys14_dc,
        time_series_in_memory = true,
    ),
    "c_sys5" => (
        description = "5-bus system",
        build = build_c_sys5,
        time_series_in_memory = true,
    ),
    "c_sys5_bat" => (
        description = "5-bus system with Storage Device",
        build = build_c_sys5_bat,
        time_series_in_memory = true,
    ),
    "c_sys5_dc" => (
        description = "Systems with HVDC data in the branches",
        build = build_c_sys5_dc,
        time_series_in_memory = true,
    ),
    "c_sys5_ed" =>
        (description = "", build = build_c_sys5_ed, time_series_in_memory = true),
    "c_sys5_hy" => (
        description = "5-bus system with HydroPower Energy",
        build = build_c_sys5_hy,
        time_series_in_memory = true,
    ),
    "c_sys5_hy_ed" =>
        (description = "", build = build_c_sys5_hy_ed, time_series_in_memory = true),
    "c_sys5_phes_ed" =>
        (description = "", build = build_c_sys5_phes_ed, time_series_in_memory = true),
    "c_sys5_hy_uc" =>
        (description = "", build = build_c_sys5_hy_uc, time_series_in_memory = true),
    "c_sys5_hyd" =>
        (description = "", build = build_c_sys5_hyd, time_series_in_memory = true),
    "c_sys5_il" => (
        description = "System with Interruptible Load",
        build = build_c_sys5_il,
        time_series_in_memory = true,
    ),
    "c_sys5_ml" =>
        (description = "", build = build_c_sys5_ml, time_series_in_memory = true),
    "c_sys5_re" => (
        description = "5-bus system with Renewable Energy",
        build = build_c_sys5_re,
        time_series_in_memory = true,
    ),
    "c_sys5_re_only" =>
        (description = "", build = build_c_sys5_re_only, time_series_in_memory = true),
    "c_sys5_uc" =>
        (description = "", build = build_c_sys5_uc, time_series_in_memory = true),
    "c_sys5_pglib" => (
        description = "5-bus with ThermalMultiStart",
        build = build_c_sys5_pglib,
        time_series_in_memory = true,
    ),
    "c_sys5_pwl_uc" => (
        description = "5-bus with SOS cost function",
        build = build_c_sys5_pwl_uc,
        time_series_in_memory = true,
    ),
    "c_sys5_pwl_ed" => (
        description = "5-bus with pwl cost function",
        build = build_c_sys5_pwl_ed,
        time_series_in_memory = true,
    ),
    "c_sys5_pwl_ed_nonconvex" => (
        description = "5-bus with SOS cost function",
        build = build_c_sys5_pwl_ed_nonconvex,
        time_series_in_memory = true,
    ),
    "c_sys5_reg" => (
        description = "5-bus with regulation devices and AGC",
        build = build_c_sys5_reg,
        time_series_in_memory = true,
    ),
    "c_ramp_test" => (
        description = "1-bus for ramp testing",
        build = build_sys_ramp_testing,
        time_series_in_memory = true,
    ),
    "c_duration_test" => (
        description = "1 Bus for durantion testing",
        build = build_duration_test_sys,
        time_series_in_memory = true,
    ),
    "c_linear_pwl_test" => (
        description = "1 Bus lineal PWL linear testing",
        build = build_pwl_test_sys,
        time_series_in_memory = true,
    ),
    "c_sos_pwl_test" => (
        description = "1 Bus lineal PWL sos testing",
        build = build_sos_test_sys,
        time_series_in_memory = true,
    ),
    "c_market_bid_cost" => (
        description = "1 bus system with MarketBidCost Model",
        build = build_pwl_marketbid_sys,
        time_series_in_memory = true,
    ),
)

build_PTDF5() = PTDF(build_system("c_sys5"))
build_PTDF14() = PTDF(build_system("c_sys14"))
build_PTDF5_dc() = PTDF(build_system("c_sys5_dc"))
build_PTDF14_dc() = PTDF(build_system("c_sys14_dc"))


text:
# These define majority voting functions;
# Majority voting function theory

cumulative_sum(x::DBool, y::Array{DBool}) =
    length(y) > 2 ? x.bit + cumulative_sum(y) : x.bit + y[1].bit
cumulative_sum(x::Array{DBool}) = cumulative_sum(x[1], x[2:end])

majorityN(x::Array{DBool}) = let y = floor(0.5 + (cumulative_sum(x) - 0.5)/length(x));
    @assert  y in [0.0, 1.0]
    return DBool(y == 1.0)
end

#=
And Voting procedure
=#

andN(x::Array{DBool}) = DBool(all([i.bit for i in x]))

#=
Or Voting procedure
=#

orN(x::Array{DBool}) = DBool(any([i.bit for i in x]))

#=
Dictator voting procedure
=#

dictatorN(x::Array{DBool}; idx=1) = x[idx]

#=
K-junta voting procedure
=#
# TODO Union statement is not typechecking

kjuntaN(x::Array{DBool}; k =10, voting=Union{andN, orN, majorityN}, kidxs=nothing) = begin(x)
    if kidxs == nothing
        kidxs = CartesianIndices(1:k)
    end
    voting(x[kidxs])
end

# TODO small utility function needs to be moved to utility.jl
*(a::DBool, b::Number) = a.bit * true

#=
Weighted Majority or linear thresholding function
=#

weightedMajorityN(x::Array{DBool}, weights; bias=0.0, threshold=0.5) = bias + sum( weights.*[i.bit for i in x]) > threshold

ext: jl


text:
# This file only forward to Pkg's command but under --project by default

function detect_user_julia_binary()::String
    # 1. use user specified julia binary
    if haskey(ENV, "JULIA_EXECUTABLE_PATH")
        return ENV["JULIA_EXECUTABLE_PATH"]
    end

    # 2. try `julia` command
    try
        user_julia = readchomp(Cmd(`julia -E "joinpath(Sys.BINDIR, Base.julia_exename())"`))
        return Meta.parse(user_julia)
    catch
        error(
            "cannot detect Julia compiler binary, " *
            "please specify Julia compiler binary " *
            "via environment variable JULIA_EXECUTABLE_PATH"
        )
    end
end

function withproject(command, glob, action_msg, compile_min=true)
    script = "using Pkg;"
    if !glob
        msg = "cannot $action_msg in global environment, use -g, --glob to $action_msg to global environment"
        script *= "(dirname(dirname(Pkg.project().path)) == joinpath(DEPOT_PATH[1], \"environments\")) && error(\"$msg\");"
    end

    script *= command
    exename = detect_user_julia_binary()

    options = []
    if compile_min
        push!(options, "--compile=min")
    end
    push!(options, "-g1", "--color=yes", "--startup-file=no", "-e", script)

    if glob
        GLOB_ENV = filter(x->x.first!="JULIA_PROJECT", ENV)
        run(setenv(Cmd([exename, options...]), GLOB_ENV))
    else
        withenv("JULIA_PROJECT"=>"@.") do
            run(Cmd([exename, options...]))
        end
    end
    return
end

"""
add package/project to the closest project.

# Arguments

- `urls`: package names or urls to add.

# Flags

- `-g, --glob`: add package to global shared environment.

"""
@cast function add(urls...; glob::Bool=false)
    isempty(urls) && error("expect a package name or url")
    packages = join(urls, " ")

    withproject(
        "pkg\"add $packages\"",
        glob,
        "install a package",
    )
end

"""
Make a package available for development. If pkg is an existing local path, that path will be recorded in the manifest and used.
Otherwise, a full git clone of pkg is made. Unlike the `dev/develop` command in Julia REPL pkg mode, `ion` will clone the package
to the dev folder of the current project. You can specify `--shared` flag to use shared `dev` folder under `~/.julia/dev`
(specified by `JULIA_PKG_DEVDIR`).

# Arguments

- `url`: URL or local path to the package.

# Flags

- `-s, --shared`: controls whether to use the shared develop folder.

"""
@cast function dev(url; shared::Bool=false, glob::Bool=false)
    shared_flag = shared ? "--shared" : "--local"
    cmd = Cmd(["-e", "using Pkg; pkg\"dev $shared_flag $url\" "])
    withenv("JULIA_PROJECT"=>"@.") do
        run(`$(Base.julia_cmd()) $cmd`)
    end
end

"""
Update a package. If no posistional argument is given, update all packages in current project.

# Arguments

- `pkg`: package name.
"""
@cast function update(pkg=""; glob::Bool=false)
    if isempty(pkg)
        cmd = "pkg\"up\""
    else
        cmd = "pkg\"up $pkg\""
    end

    withproject(cmd, glob, "update dependencies")
end

"""
build package/project/environment

# Arguments

- `pkg`: package name.

# Flags

- `-v, --verbose`: print verbose log
- `-g, --glob`: enable to build in global shared environment
"""
@cast function build(pkg=""; verbose::Bool=false, glob::Bool=false)
    if isempty(pkg)
        cmd = "pkg\"build\""
    else
        cmd = "pkg\"build $pkg\""
    end

    withproject(cmd, glob, "build", false)
end

"""
test package/project

# Arguments

- `pkg`: package name.

# Flags

- `-g, --glob`: enable to test in global shared environment
"""
@cast function test(pkg=""; glob::Bool=false)
    if isempty(pkg)
        cmd = "Pkg.test()"
    else
        cmd = "Pkg.test(\"$pkg\")"
    end

    withproject(cmd, glob, "test", false)
end

"""
show current environment status

# Arguments

- `pkg`: package name.

# Flags

- `-d, --diff`: show diff to last git commit
- `-g, --glob`: enable to show status in global shared environment
"""
@cast function status(pkg=""; diff::Bool=false, glob::Bool=false)
    if isempty(pkg)
        cmd = "Pkg.status(;diff=$diff)"
    else
        cmd = "Pkg.status(\"$pkg\"; diff=$diff)"
    end

    withproject(cmd, glob, "show status")
end

"""
Update the current manifest. It will update manifest
with potential changes to the dependency graph from
packages that are tracking a path.

# Flags
- `-g, --glob`: enable to resolve in global shared environment
"""
@cast function resolve(; glob::Bool=false)
    withproject("Pkg.resolve()", glob, "resolve dependencies")
end

"""
Remove a package from the current project. If the mode of pkg is PKGMODE_MANIFEST also
remove it from the manifest including all recursive dependencies of pkg.

# Arguments

- `pkg`: package name to remove.

# Flags

- `-g, --glob`: enable to remove package in global shared environment
"""
@cast function rm(pkg, pkgs...; glob::Bool=false)
    withproject("pkg\"rm $(join([pkg, pkgs...], " "))\"", glob, "rm package")
end

"""
Undoes the latest change to the active project. Only states in the current session are
stored, up to a maximum of 50 states.

# Flags

- `-g, --glob`: enable to execute in global shared environment
"""
@cast undo(; glob::Bool=false) = withproject("Pkg.undo()", glob, "undo")

"""
Redoes the changes from the latest undo.

# Flags

- `-g, --glob`: enable to execute in global shared environment
"""
@cast redo(; glob::Bool=false) = withproject("Pkg.redo()", glob, "redo")

"""
If a Manifest.toml file exists in the active project, download all the packages
declared in that manifest. Otherwise, resolve a set of feasible packages from the
Project.toml files and install them. If no Project.toml exist
in the current active project, create one with all the dependencies in the manifest
and instantiate the resulting project.

# Flags

- `-v, --verbose`: prints the build output to stdout/stderr instead of redirecting to the build.log file.
"""
@cast function instantiate(;verbose::Bool=false)
    withproject("Pkg.instantiate(;verbose=$verbose)", false, "instantiate")
end

"""
If pkg is pinned, remove the pin. If pkg is tracking a path, e.g. after Pkg.develop,
go back to tracking registered versions.

# Flags

- `-g, --glob`: enable to execute in global shared environment

"""
@cast function free(pkg; glob::Bool=false)
    withproject("Pkg.free(\"$pkg\")", glob, "free $pkg")
end

# @cast function register()
# end


text:
module LocalMethods

using SpecialFunctions: besselk, gamma

#M_vρσ is the matern covariance and G_σρ is Gaussian autocovariance model
export ℳ_ν, M_νρσ, G_σρ

tν𝒦t(t,ν) = t^ν * besselk(ν, t)

"""
    ℳ_ν(t, ν)

Compute (√(2ν)*t) ^ ν * 𝒦ν(√(2ν)*t) where 𝒦ν is the modified Bessel function of the second kind of order ν.

# Example
```julia-repl
julia> ℳ_ν(0.0, 0.5)
1.0
```
"""
function ℳ_ν(t, ν)
	pt, pν, p0, p1 = promote(t, ν, 0, 1)
	return (pt==p0) ? p1 : tν𝒦t(√(2pν)*pt,pν) * 2^(1-pν) / gamma(pν)
end

function M_νρσ(t; ν, ρ, σ)
	return σ * σ * ℳ_ν(t/ρ, ν)
end

function G_σρ(t; σ,ρ)
	return σ * σ *exp(-t^2/ρ^2)

end
end # end module


text:
## Foil type
#==========================================================================================#

abstract type AbstractFoil end

"""
    Foil(x, y, name = "")
    Foil(coordinates, name = "")

Structure consisting of foil coordinates in 2 dimensions with an optional name.

The coordinates should be provided in counter-clockwise format, viz. from the trailing edge of the upper surface to the trailing edge of the lower surface.
"""
struct Foil{T <: Real} <: AbstractFoil
    x      :: Vector{T}
    y      :: Vector{T}
    name   :: String
end

function Foil(xs, ys, name = "")
    @assert length(xs) == length(ys) "Lengths of xs and ys must match!"

    T = promote_type(eltype(xs), eltype(ys))
    Foil{T}(xs, ys, name)
end

Foil(coords :: Vector{<: FieldVector{2,<: Real}}, name = "") = Foil(getindex.(coords, 1), getindex.(coords, 2), name)
Foil((name, coords)) = Foil(coords, name)

function Foil(coords :: AbstractMatrix{<: Real}, name = "")
    @assert size(coords, 2) == 2 "The array must have only two columns for coordinates!"

    @views Foil(coords[:,1], coords[:,2], name)
end

function Base.show(io :: IO, foil :: Foil)
    print(io, foil.name, " ", typeof(foil), " with ", length(foil.x), " points.")
end

## Foil processing
#==========================================================================================#

"""
    coordinates(foil :: Foil)

Generate the array of `Foil` coordinates.
"""
coordinates(foil :: Foil) = [ foil.x foil.y ]

"""
    arc_length(foil :: Foil)

Compute the arc-length of a `Foil`.
"""
arc_length(foil :: Foil) = let c = coordinates(foil); @views norm(c[2:end,:] .- c[1:end-1,:]) end
chord_length(foil :: Foil) = maximum(foil.x)

"""
    scale(foil :: Foil, scale)

Scale the coordinates of a `Foil` to a scaling value.
"""
scale(foil :: Foil, scale) = Foil(scale .* coordinates(foil), foil.name)

translate(foil :: Foil; vector) = Foil(foil.x .+ vector[1], foil.y .+ vector[2])

function rotate(foil :: Foil; angle, center = zeros(2))
    trans  = @views [ foil.x .- center[1] foil.y .- center[2] ]     # Translate
    rotate = trans * RotMatrix{2,eltype(angle)}(deg2rad(angle))'    # Rotate
    @views Foil(rotate[:,1] .+ center[1], rotate[:,2] .+ center[2], foil.name) # Inverse translate
end

function interpolate(foil :: Foil, xs)
    upper, lower = split_surface(foil)

    y_u = @views LinearInterpolation(upper[:,1], upper[:,2]).(xs)
    y_l = @views LinearInterpolation(lower[:,1], lower[:,2]).(xs)

    @views Foil([ xs[end:-1:2]; xs ], [ y_u[end:-1:2]; y_l ], foil.name)
end

affine(foil :: Foil; angle, vector) = translate(rotate(foil; angle = angle); vector = vector)

"""
    camber_thickness(foil :: Foil, num :: Integer)

Compute the camber-thickness distribution of a `Foil` with cosine spacing.
"""
camber_thickness(foil :: Foil, num = 40) = coordinates_to_camber_thickness(foil, num + 1)

leading_edge_index(foil :: Foil) = argmin(coordinates(foil)[:,1])

"""
    upper_surface(foil :: Foil)

Get the upper surface coordinates of a `Foil` from leading to trailing edge.
"""
upper_surface(foil :: Foil) = @views reverse(coordinates(foil)[1:leading_edge_index(foil),:], dims = 1)

"""
    lower_surface(foil :: Foil)

Get the lower surface coordinates of a `Foil` from leading to trailing edge.
"""
lower_surface(foil :: Foil) = @view coordinates(foil)[leading_edge_index(foil):end,:]

"""
    split_surface(foil :: Foil)

Split the `Foil` coordinates into upper and lower surfaces.
"""
split_surface(foil :: Foil) = upper_surface(foil), lower_surface(foil)

"""
    cosine_interpolation(foil :: Foil, n :: Integer = 40)

Interpolate a `Foil` profile's coordinates to a cosine by projecting the x-coordinates of a circle onto the geometry with ``2n`` points.
"""
function cosine_interpolation(foil :: Foil, n :: Integer = 40)
    x_min, x_max = extrema(foil.x)
    x_circ = cosine_spacing((x_min + x_max) / 2, x_max - x_min, n)

    interpolate(foil, x_circ)
end

function camber_line(foil :: Foil, n = 60)
    upper, lower = split_surface(foil)
    xs  = LinRange(minimum(foil.x), maximum(foil.x), n + 1)
    y_u = @views LinearInterpolation(upper[:,1], upper[:,2]).(xs)
    y_l = @views LinearInterpolation(lower[:,1], lower[:,2]).(xs)

    [ xs (y_u + y_l) / 2 ]
end

"""
    make_panels(foil :: Foil)
    make_panels(foil :: Foil, n :: Integer)

Generate a vector of `Panel2D`s from a `Foil`, additionally with cosine interpolation using (approximately) ``n`` points if provided.
"""
function make_panels(foil :: Foil)
    vecs = SVector.(foil.x, foil.y)
    @views Panel2D.(vecs[2:end,:], vecs[1:end-1,:])[end:-1:1]
end

make_panels(foil :: Foil, n :: Integer) = make_panels(cosine_interpolation(foil, n ÷ 2))

function camber(foil :: Foil, x_by_c)
    upper, lower = split_surface(foil)
    y_u = @views LinearInterpolation(upper[:,1], upper[:,2])(x_by_c * chord_length(foil))
    y_l = @views LinearInterpolation(lower[:,1], lower[:,2])(x_by_c * chord_length(foil))

    (y_u + y_l) / 2
end

function control_surface(foil :: Foil, δ, xc_hinge)
    y_hinge  = camber(foil, xc_hinge)
    rot_foil = rotate(foil; angle = -δ, center = [xc_hinge, y_hinge])
    coords   = coordinates(foil)
    @views coords[foil.x .>= xc_hinge,1] = rot_foil.x[foil.x .>= xc_hinge]
    @views coords[foil.x .>= xc_hinge,2] = rot_foil.y[foil.x .>= xc_hinge]
    Foil(coords, foil.name)
end

"""
    control_surface(foil :: Foil; angle, hinge)

Modify a `Foil` to mimic a control surface by specifying a deflection angle (in degrees, clockwise-positive convention) and a normalized hinge ``x``-coordinate ``∈ [0,1]`` in terms of the chord length.
"""
control_surface(foil :: Foil; angle, hinge) = control_surface(foil, angle, hinge)

## Camber-thickness representation
#==========================================================================================#

"""
    coordinates_to_camber_thickness(coords, n = 40)

Convert 2-dimensional coordinates to its camber-thickness representation after cosine interpolation with ``2n`` points.
"""
function coordinates_to_camber_thickness(foil, n = 40)
    # Cosine interpolation and splitting
    upper, lower = split_surface(cosine_interpolation(foil, n))

    camber       = @views (upper[:,2] + lower[:,2]) / 2
    thickness    = @views  upper[:,2] - lower[:,2]

    @views [ upper[:,1] camber thickness ]
end

"""
    camber_thickness_to_coordinates(xs, camber, thickness)

Convert the camber-thickness representation to 2-dimensional coordinates given the ``x``-locations and their corresponding camber and thickness values.
"""
camber_thickness_to_coordinates(xs, camber, thickness) =
    @views [ [xs camber + thickness / 2][end:-1:1,:];
              xs camber - thickness / 2             ]

camber_thickness_to_coordinates(coords) = @views camber_thickness_to_coordinates(coords[:,1], coords[:,2], coords[:,3])

"""
    camber_coordinates(coords :: Array{2, <: Real})

Generate the camber coordinates on the ``x``-``z`` plane at ``y = 0``.
"""
camber_coordinates(coords) = @views [ coords[:,1] zero(coords[:,1]) coords[:,2] ]

"""
    thickness_coordinates(coords :: Array{2, <: Real})

Generate the thickness coordinates on the ``x``-``z`` plane at ``y = 0``.
"""
thickness_coordinates(coords) = @views [ coords[:,1] zero(coords[:,1]) coords[:,3] ]

function maximum_thickness_to_chord(coords)
    xs, thiccs    = @views coords[:,1], coords[:,3]
    max_thick_arg = argmax(thiccs)
    chord         = @views maximum(coords[:,1])
    @views xs[max_thick_arg] / chord, thiccs[max_thick_arg] / chord
end

"""
    read_foil(path :: String; header = true; name = "")

Generate a `Foil` from a file consisting of 2D coordinates with named arguments to skip the header (first line of the file) or assign a name.

By default, the header is assumed to exist and should contain the airfoil name, which is assigned to the name of the `Foil`.
"""
function read_foil(path :: String; name = "")
    coords, foil_name = readdlm(path, header = true)
    if name != ""
        return Foil(coords, name)
    else
        return Foil(coords, strip(join(foil_name, " ")))
    end
end


text:
"""
`boxplot(data; nargs...)` → `Plot`

Description
============

Draws a box-and-whisker plot.
The first argument specifies the data to plot. This is a vector of vectors, with each
inner vector representing a data series. We use a vector of vectors over a matrix
to allow series of different lengths.
Optionally, a list of labels may be provided, with length equal to the number of series.
Alternatively, one can specify a boxplot using a dictionary.
In that case, the values, which have to be numeric, will be used as the data series,
and the keys, which have to be strings, will be used as the labels.

Usage
======

    boxplot(data; labels = [" " for _ in 1:size(data, 1)], border = :solid, title = "",
            margin = 3, padding = 1, color = :blue, width = 40,
            left=minimum(map(minimum, data)) - 1, right=maximum(map(maximum, data)) + 1)

    boxplot(dictionary; nargs...)

Arguments
==========

- **`data`** : The data the box plot is based on. A vector of vectors, with each
  inner vector representing a data series. Choose a vector of vectors over a matrix
  to allow series of different lengths.

- **`labels`** : A list of labels for the data series. Must be the same length as the number of series.

- **`dictionary`** : A dictonary in which the keys will be used as `labels`
  and the values will be utilized as `data`.

- **`border`** : The style of the bounding box of the plot.
  Supports `:solid`, `:bold`, `:dashed`, `:dotted`, `:ascii`, and `:none`.

- **`title`** : Text to display on the top of the plot.

- **`margin`** : Number of empty characters to the left of the whole plot.

- **`padding`** : Space of the left and right of the plot between the labels and the canvas.

- **`color`** : Colour of the drawing. Can be any of `:black`, `:blue`, `:cyan`,
  `:green`, `:magenta`, `:red`, `:yellow`, `:white`, or a light version of the above (`:light_colour`).
  By default no colouring is applied.

- **`width`** : Number of characters per row that should be used for plotting.

- **`left`** : The value of the left-hand edge of the plot.

- **`right`** : The value of the right-hand edge of the plot.

Returns
========

A plot object of type `Plot{BoxplotGraphics}`

Author(s)
==========

- Matthew Lake (Github: https://github.com/mgtlake)

Examples
=========

    julia> boxplot([1, 2, 3, 7], title = "Test")

                       Test
     ┌────────────────────────────────────────┐
     │    │   ┌──┬───────┐              │     │
     │    ├───┤  │       ├──────────────┤     │
     │    │   └──┴───────┘              │     │
     └────────────────────────────────────────┘
     0                  4.0                   8

See also
=========

`Plot`, `histogram`, `BoxplotGraphics`
"""
function boxplot(
        data::AbstractVector;
        labels::AbstractVector{<:AbstractString} = [" " for _ in 1:size(data, 1)],
        border = :solid,
        title::AbstractString = "",
        margin::Int = 3,
        padding::Int = 1,
        color::Symbol = :normal,
        width::Int = 40,
        left::Number = minimum(map(minimum, data)) - 1,
        right::Number = maximum(map(maximum, data)) + 1)
    margin >= 0 || throw(ArgumentError("Margin must be greater than or equal to 0"))
    length(labels) == length(data) || throw(DimensionMismatch("Wrong number of labels"))

    width = max(width, 10)

    if left > minimum(map(minimum, data)) || right < maximum(map(maximum, data))
        throw(ArgumentError("Plot range ($left, $right) too restrictive for data range ($(min(data...)), $(max(data...)))"))
    end

    area = BoxplotGraphics(data[1], width, color = color,
                           left = left, right = right, labels = labels)
    for i in 2:length(data)
        addseries!(area, data[i])
    end

    new_plot = Plot(area, title = title, margin = margin,
                   padding = padding, border = border)

    annotate!(new_plot, :bl, string(left))
    annotate!(new_plot, :b, string((left + right) / 2))
    annotate!(new_plot, :br, string(right))

    for (i, label) in enumerate(labels)
        # Find end of last 3-line region, then add 2 for centre of current
        annotate!(new_plot, :l, (i-1)*3+2, label)
    end

    new_plot
end

"""
`boxplot!(plot, data; nargs)` → `Plot`

Mutating variant of `boxplot`, in which the first parameter (`plot`) specifies
the existing plot to draw on.

See `boxplot` for more information.
"""
function boxplot!(
        plot::Plot{<:BoxplotGraphics},
        data::AbstractVector{<:Number};
        label = " ",
        kw...)
    !isempty(data)|| throw(ArgumentError("Can't append empty array to boxplot"))

    if plot.graphics.left > minimum(data)
        plot.graphics.left = minimum(data) - 1
    end
    if plot.graphics.right < maximum(data)
        plot.graphics.right = maximum(data) + 1
    end

    addseries!(plot.graphics, data)

    # Find end of last 3-line region, then add 2 for centre of current
    annotate!(plot, :l, (length(plot.graphics.data)-1)*3+2, label)

    annotate!(plot, :bl, string(plot.graphics.left))
    annotate!(plot, :b, string((plot.graphics.left + plot.graphics.right) / 2))
    annotate!(plot, :br, string(plot.graphics.right))

    plot
end

function boxplot(data::AbstractVector{<:Number}; kw...) where {T, N <: Number}
    boxplot([data]; kw...)
end

function boxplot(dict::Dict; kw...) where {T, N <: Number}
    boxplot(collect(values(dict)); labels=collect(keys(dict)), kw...)
end


text:
using GeometricEquations: function_dummy_v, initial_multiplier, symplectic_matrix,
                          promote_tspan, promote_tspan_and_tstep, parameter_types
using Test

include("initial_conditions.jl")

@test function_dummy_v(t₀, q₀, λ₀) === nothing
@test function_dummy_v(t₀, q₀, p₀, λ₀) === nothing

@test promote_tspan(0) == (0,0)
@test promote_tspan(1) == (0,1)
@test promote_tspan(0.0) == (0.0,0.0)
@test promote_tspan(1.0) == (0.0,1.0)
@test promote_tspan((0,   1.0)) == (0.0,1.0)
@test promote_tspan((0.0, 1  )) == (0.0,1.0)
@test promote_tspan([0,   1.0]) == (0.0,1.0)
@test promote_tspan([0.0, 1  ]) == (0.0,1.0)
@test promote_tspan([0,   1  ]) == (0,  1  )
@test promote_tspan([0.0, 1.0]) == (0.0,1.0)

@test_throws ErrorException promote_tspan([0.0])
@test_throws ErrorException promote_tspan([0.0, 0.5, 1.0])

@test promote_tspan_and_tstep((0,   1.0), 0  ) == ((0.0,1.0), 0.0)
@test promote_tspan_and_tstep((0,   1.0), 0.0) == ((0.0,1.0), 0.0)
@test promote_tspan_and_tstep((0.0, 1  ), 0  ) == ((0.0,1.0), 0.0)
@test promote_tspan_and_tstep((0.0, 1  ), 0.0) == ((0.0,1.0), 0.0)
@test promote_tspan_and_tstep((0,   1  ), 0  ) == ((0,  1  ), 0  )
@test promote_tspan_and_tstep((0,   1  ), 0.0) == ((0.0,1.0), 0.0)
@test promote_tspan_and_tstep((0.0, 1.0), 0  ) == ((0.0,1.0), 0.0)
@test promote_tspan_and_tstep((0.0, 1.0), 0.0) == ((0.0,1.0), 0.0)

params = (a = 1.0, b = 0, c = :c, d = [0,1])
types  = (a = Float64, b = Int, c = Symbol, d = Vector{Int})
@test parameter_types(params) == types

zero_vec = [zeros(3) for i in 1:3]

@test initial_multiplier(zeros(3), zeros(3)) == zeros(3)
@test initial_multiplier(zeros(3), zero_vec) == zeros(3)
@test initial_multiplier(zero_vec, zeros(3)) == zero_vec
@test initial_multiplier(zero_vec, zero_vec) == zero_vec

Ω = symplectic_matrix(Float64, 4)

@test Ω == -Ω'
@test Ω == -inv(Ω)

@test symplectic_matrix(eltype(q₀), length(q₀)) == symplectic_matrix(t₀, q₀, p₀)


text:
#### Operators for a system with a body and stationary points ####

export temperature_bc_constraint_rhs!, heatflux_bc_constraint_rhs!

function temperature_bc_constraint_rhs!(Ts::ScalarData{N},sys::HeatConduction{NX,NY,N},t::Real) where {NX,NY,N}
    surface_temperature!(sys.ΔTs,sys,t)
    Ts .= sys.ΔTs
    return Ts
end

function heatflux_bc_constraint_rhs!(qs::ScalarData{N},sys::HeatConduction{NX,NY,N},t::Real) where {NX,NY,N}
    fill!(qs,0.0)
    # for now, it's adiabatic only
    return qs
end

# Constraint operators, using stored regularization and interpolation operators
# B₁ᵀ = Eᵀ
function temperature_constraint_force!(out::Nodes{Primal,NX,NY},σ::ScalarData{N},sys::HeatConduction{NX,NY,N}) where {NX,NY,N}
    out .= sys.Rc*σ
    return out
end

function temperature_constraint_force!(out::Nodes{Primal,NX,NY},σ,sys::HeatConduction{NX,NY,0}) where {NX,NY}
    return out
end

function heatflux_constraint_force!(out::Nodes{Primal,NX,NY},τ::ScalarData{N},sys::HeatConduction{NX,NY,N}) where {NX,NY,N}
    @unpack params, normals, Rf, Vb, Vf = sys
    @unpack d, α = params
    fact = α/cellsize(sys)
    #fact = isnothing(d) ? fact : fact/d

    fill!(out,0.0)
    product!(Vb,normals,τ)
    Vf .= Rf*Vb
    divergence!(out,Vf)
    out .*= fact
    return out
end

function heatflux_constraint_force!(out::Nodes{Primal,NX,NY},σ,sys::HeatConduction{NX,NY,0}) where {NX,NY}
    return out
end

#=
function _vel_ns_op_constraint_force!(u::Edges{Primal,NX,NY},τ,sys::HeatConduction{NX,NY,0}) where {NX,NY}
    return u
end

function _vel_ns_op_constraint_force!(u::Edges{Primal,NX,NY},τ::VectorData{N},sys::HeatConduction{NX,NY,N}) where {NX,NY,N}
    u .= sys.Rf*τ
    return u
end
=#

# B₂ = E
function temperature_constraint_op!(out::ScalarData{N},T::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY,N}) where {NX,NY,N}
    out .= sys.Ec*T
    return out
end

function heatflux_constraint_op!(out::ScalarData{N},T::Nodes{Primal,NX,NY},sys::HeatConduction{NX,NY,N}) where {NX,NY,N}
    @unpack params, normals, Ef, gradTb, gradTf = sys
    @unpack k = params
    fill!(gradTf,0.0)

    grad!(gradTf,T)
    gradTb .= Ef*gradTf
    pointwise_dot!(out,gradTb,normals)
    out .*= -k/cellsize(sys)
    return out
end