text:
Europe = [
    # "Country"                   "GADM"
      "Albania"                   GADM("Albania")
      "Armenia"                   GADM("Armenia")
      "Austria"                   GADM("Austria")
      "Azerbaijan"                GADM("Azerbaijan")
      "Belarus"                   GADM("Belarus")
      "Belgium"                   GADM("Belgium")
      "Bosnia and Herzegovina"    GADM("Bosnia and Herzegovina")
      "Bulgaria"                  GADM("Bulgaria")
      "Croatia"                   GADM("Croatia")
      "Cyprus"                    GADM("Cyprus")
      "Czech Republic"            GADM("Czech Republic")
      "Denmark"                   GADM("Denmark")
      "Estonia"                   GADM("Estonia")
      "Finland"                   GADM("Finland")
      "France"                    GADM("France")
      "Germany"                   GADM("Germany")
      "Greece"                    GADM("Greece")
      "Hungary"                   GADM("Hungary")
      "Iceland"                   GADM("Iceland")
      "Ireland"                   GADM("Ireland")
      "Italy"                     GADM("Italy")
      "Kosovo"                    GADM("Kosovo")
      "Latvia"                    GADM("Latvia")
      "Lithuania"                 GADM("Lithuania")
      "Luxembourg"                GADM("Luxembourg")
      "Malta"                     GADM("Malta")
      "Moldova"                   GADM("Moldova")
      "Montenegro"                GADM("Montenegro")
      "Netherlands"               GADM("Netherlands")
      "North Macedonia"           GADM("Macedonia")
      "Norway"                    GADM("Norway")
      "Poland"                    GADM("Poland")
      "Portugal"                  GADM("Portugal")
      "Romania"                   GADM("Romania")
      "Russia"                    GADM("Russia")
      "Serbia"                    GADM("Serbia")
      "Slovakia"                  GADM("Slovakia")
      "Slovenia"                  GADM("Slovenia")
      "Spain"                     GADM("Spain")
      "Sweden"                    GADM("Sweden")
      "Switzerland"               GADM("Switzerland")
      "Turkey"                    GADM("Turkey")
      "Ukraine"                   GADM("Ukraine")
      "United Kingdom"            GADM("United Kingdom")
  ]

Regions_dict = Dict( "Nordics" => (["Finland", "Sweden", "Norway", "Denmark"]),
  "Central" => (["Germany", "Austria", "Switzerland", "Czech Republic"]),
  "Western" => (["France", "United Kingdom", "Ireland", "Netherlands", "Belgium", "Luxembourg"]),
  "Mediterranian" => (["Spain", "Portugal", "Italy", "Greece", "Croatia", "Malta", "Albania", "Bosnia and Herzegovina"]),
  "Eastern" => (["Poland", "Slovakia", "Hungary", "Lithuania", "Latvia", "Estonia"]))

"""Takes the list of countries specified in run_data_generation.jl and turns them into a format that GlobalEnergyGIS can work with.
# Arguments
- `Regions`: List of regions used for the model.
"""
# function if one wants specific countries
# Check if the regions specified in run_data_generation.jl are in Europe or in Regions_dict
# Get corresponding values
function get_countries(Regions)
    if sum(occursin.(Regions[1,1], Europe[:,1]))>0
        n = length(Europe[:,1])
        m = length(Regions)
        Dataset_Countries = []
        GADM_Name = []
        for j in 1:m
            for i in 1:n
                    if occursin.(Regions[j,1], Europe[i,1])
                        Dataset_Countries = [Dataset_Countries; permutedims(Regions[j,:])]
                        GADM_Name = [GADM_Name; GADM(Europe[i])]
                    end
            end
        end
        Countries = hcat(Dataset_Countries, GADM_Name)
    elseif sum(occursin.(Regions[1,1], collect(keys(Regions_dict))))>0
        Keys = []
        Values = []
        m = length(Regions)
        R = getindex.(Ref(Regions_dict),(Regions))
        for i in 1:m
            RVal = getindex.(Ref(Regions_dict),(Regions))[i]
            Keys = [Keys; Regions[i]]
            Values = [Values; GADM(RVal...)]
        end
        GADM_List = hcat(Keys, Values)
    end
end

# function if one wants to create whole regions
# if you want all countries in Europe
function get_countries()
   Europe
end

function replace_nans!(array::Array{Float64, N}) where N
    for i in eachindex(array)
        if isnan(array[i])
            array[i] = zero(1)
        end
    end
end

"""Uses wind, solar and hydro data produced by the GlobalEnergyGIS package to generate CSV files with the data used in the model.
# Arguments
- `inputdata`: Path to the output folder of GlobalEnergyGIS.
- `sspscenario_input`: Specified ssp scenario used to generate data.
- `sspyear_input`: Target year of the model.
- `era_year_input`: Year of the base data.
- `Dataset`: Countries generated by the function `get_countries`.
- `folder`: Sets the folder in which the generated data is put.
- `subfolder`: Sets the subfolder in which the generated data is put.
- `instance`: Sets the name of the instance for which the data is being generated.
- `T`: Time period of the data set in hours.
- `t`: number of technologies used in the model
"""

function create_data_sets(inputdata, sspscenario_input, sspyear_input, era_year_input, Dataset, folder, subfolder, instance, T, t, Fmin)

    # create path for the instance
    structure_path = joinpath(folder,subfolder)
    instance_path = mkpath(joinpath(structure_path, instance))
    nodes_path =  mkpath(joinpath(instance_path, "nodes"))

# Generate datasets for regions using GlobalEnergyGIS
    # define regions and countries and name the instance

   saveregions(instance, Dataset)
   makedistances(instance)
   createmaps(instance)

    # generate VRE data and demand
   GISsolar(gisregion = instance)
   GISwind(gisregion = instance)
   GIShydro(gisregion = instance)
   predictdemand(gisregion = instance, sspscenario=sspscenario_input, sspyear=sspyear_input, era_year=era_year_input)

# function get data
    # Name of the instance you defined in Inputdata.jl
    SolarData = matread(joinpath(inputdata, "GISdata_solar$(era_year_input)_$instance.mat"))
    WindData = matread(joinpath(inputdata, "GISdata_wind$(era_year_input)_$instance.mat"))
    HydroData = matread(joinpath(inputdata, "GISdata_hydro_$instance.mat"))
    TransmissionData = matread(joinpath(inputdata, "distances_$instance.mat"))

    Regionlist = typeof(TransmissionData["regionlist"]) == Float64 ? [TransmissionData["regionlist"]] : TransmissionData["regionlist"]
    n = length(Regionlist)

    # Solar

    PVCapA = typeof(SolarData["capacity_pvplantA"]) == Float64 ? [SolarData["capacity_pvplantA"]] : SolarData["capacity_pvplantA"]
    PVCapB = typeof(SolarData["capacity_pvplantB"]) == Float64 ? [SolarData["capacity_pvplantB"]] : SolarData["capacity_pvplantB"]
    CSPCapA = typeof(SolarData["capacity_cspplantA"]) == Float64 ? [SolarData["capacity_cspplantA"]] : SolarData["capacity_cspplantA"]
    CSPCapB = typeof(SolarData["capacity_cspplantB"]) == Float64 ? [SolarData["capacity_cspplantB"]] : SolarData["capacity_cspplantB"]
    RooftopCap = typeof(SolarData["capacity_pvrooftop"]) == Float64 ? [SolarData["capacity_pvrooftop"]] : SolarData["capacity_pvrooftop"]
    SolarCapTotal = PVCapA .+ PVCapB .+ CSPCapA .+ CSPCapB .+ RooftopCap
    replace_nans!(PVCapA)
    replace_nans!(PVCapB)
    replace_nans!(CSPCapA)
    replace_nans!(CSPCapB)
    replace_nans!(RooftopCap)
    replace_nans!(SolarCapTotal)

    CFtime_PVA = typeof(SolarData["CFtime_pvplantA"]) == Float64 ? [SolarData["CFtime_pvplantA"]] : SolarData["CFtime_pvplantA"]
    CFtime_PVB = typeof(SolarData["CFtime_pvplantB"]) == Float64 ? [SolarData["CFtime_pvplantB"]] : SolarData["CFtime_pvplantB"]
    CFtime_CSPA = typeof(SolarData["CFtime_cspplantA"]) == Float64 ? [SolarData["CFtime_cspplantA"]] : SolarData["CFtime_cspplantA"]
    CFtime_CSPB = typeof(SolarData["CFtime_cspplantB"]) == Float64 ? [SolarData["CFtime_cspplantB"]] : SolarData["CFtime_cspplantB"]
    CFtime_PVrooftop = typeof(SolarData["CFtime_pvrooftop"]) == Float64 ? [SolarData["CFtime_pvrooftop"]] : SolarData["CFtime_pvrooftop"]
    CFtime_solar = CFtime_PVA .+ CFtime_PVB .+ CFtime_CSPA .+ CFtime_CSPB .+ CFtime_PVrooftop
    replace_nans!(CFtime_PVA)
    replace_nans!(CFtime_PVB)
    replace_nans!(CFtime_CSPA)
    replace_nans!(CFtime_CSPB)
    replace_nans!(CFtime_PVrooftop)
    replace_nans!(CFtime_solar)

    #Preallocate arrays for results
    avail_sol_cspA = zeros(Float64, T, n)
    avail_sol_cspB = zeros(Float64, T, n)
    avail_sol_pvA = zeros(Float64, T, n)
    avail_sol_pvB = zeros(Float64, T, n)
    avail_sol_pvrooftop = zeros(Float64, T, n)

    #Calculate absolute availability values and sum up the different classes for each solar power type
    for i in 1:T
        avail_sol_cspA[i,:,:] = permutedims(sum(CSPCapA .* CFtime_CSPA[i,:,:], dims=2))
        avail_sol_cspB[i,:,:] = permutedims(sum(CSPCapB .* CFtime_CSPB[i,:,:], dims=2))
        avail_sol_pvA[i,:,:] = permutedims(sum(PVCapA .* CFtime_PVA[i,:,:], dims=2))
        avail_sol_pvB[i,:,:] = permutedims(sum(PVCapB .* CFtime_PVB[i,:,:], dims=2))
        avail_sol_pvrooftop[i,:,:] = permutedims(sum(RooftopCap .* CFtime_PVrooftop[i,:,:], dims=2))
    end
    #Calculate total solar capacity
    capacity_solar = permutedims(sum(CSPCapA, dims=2)) + permutedims(sum(CSPCapB, dims=2)) + permutedims(sum(PVCapA, dims=2)) +
                    permutedims(sum(PVCapB, dims=2))  + permutedims(sum(RooftopCap, dims=2))
    #Sum the different solar types and divide with total capacity to get relative availability values
    avail_sol = (avail_sol_cspA + avail_sol_cspB + avail_sol_pvA + avail_sol_pvB + avail_sol_pvrooftop) ./ capacity_solar

    # Wind
    WindCapA = typeof(WindData["capacity_onshoreA"]) == Float64 ? [WindData["capacity_onshoreA"]] : WindData["capacity_onshoreA"]
    WindCapB = typeof(WindData["capacity_onshoreB"]) == Float64 ? [WindData["capacity_onshoreB"]] : WindData["capacity_onshoreB"]
    WindCapTotal = WindCapA .+ WindCapB
    replace_nans!(WindCapA)
    replace_nans!(WindCapB)
    replace_nans!(WindCapTotal)

    WindCapOff = typeof(WindData["capacity_offshore"]) == Float64 ? [WindData["capacity_offshore"]] : WindData["capacity_offshore"]
    replace_nans!(WindCapOff)

    ## Wind Onshore
    CFtime_windonshoreA = typeof(WindData["CFtime_windonshoreA"]) == Float64 ? [WindData["CFtime_windonshoreA"]] : WindData["CFtime_windonshoreA"]
    CFtime_windonshoreB = typeof(WindData["CFtime_windonshoreB"]) == Float64 ? [WindData["CFtime_windonshoreB"]] : WindData["CFtime_windonshoreB"]
    CFtime_windonshore = CFtime_windonshoreA .+ CFtime_windonshoreB
    replace_nans!(CFtime_windonshoreA)
    replace_nans!(CFtime_windonshoreB)
    replace_nans!(CFtime_windonshore)

    ## Wind Offshore
    CFtime_windoffshore = typeof(WindData["CFtime_windoffshore"]) == Float64 ? [WindData["CFtime_windoffshore"]] : WindData["CFtime_windoffshore"]
    replace_nans!(CFtime_windoffshore)

    #Preallocate arrays for results
    avail_wind_offshore = zeros(Float64, T, n)
    avail_wind_onshoreA = zeros(Float64, T, n)
    avail_wind_onshoreB = zeros(Float64, T, n)

    # Calculate absolute availability values and sum up the different classes for each wind power type
    for i in 1:T
        avail_wind_offshore[i,:,:] = permutedims(sum(WindCapOff .* CFtime_windoffshore[i,:,:], dims=2))
        avail_wind_onshoreA[i,:,:] = permutedims(sum(WindCapA .* CFtime_windonshoreA[i,:,:], dims=2))
        avail_wind_onshoreB[i,:,:] = permutedims(sum(WindCapB .* CFtime_windonshoreB[i,:,:], dims=2))
    end
    capacity_wind_on = permutedims(sum(WindCapA, dims=2)) + permutedims(sum(WindCapB, dims=2))
    capacity_wind_off = permutedims(sum(WindCapOff, dims=2))
    avail_wind_on = (avail_wind_onshoreA + avail_wind_onshoreB) ./ capacity_wind_on
    avail_wind_off = avail_wind_offshore ./ capacity_wind_off
    replace_nans!(avail_wind_off)

    # demand
    gisdemand = JLD.load(joinpath(inputdata, "SyntheticDemand_$(instance)_$(sspscenario_input)-$(sspyear_input)_$(era_year_input).jld"), "demand")

    demand = zeros(Float64, T, n)
    for i in 1:T
        demand[i,:] = gisdemand[i,:]
    end

    # Transmission distances and connected nodes
    Distances = typeof(TransmissionData["distances"]) == Float64 ? [TransmissionData["distances"]] : TransmissionData["distances"]
    connected = typeof(TransmissionData["connected"]) == Float64 ? [TransmissionData["connected"]] : TransmissionData["connected"]

    # Hydro
    HydroCap = typeof(HydroData["existingcapac"]) == Float64 ? [HydroData["existingcapac"]] : HydroData["existingcapac"]
    replace_nans!(HydroCap)

    existingcapac = permutedims(HydroData["existingcapac"]) .* 1000
    potentialcapacities = HydroData["potentialcapac"]
    potentialcapac = []
    for i in 1:n
        potentialcapac = vcat(potentialcapac, sum(potentialcapacities[i,:,:])*1000 + existingcapac[1,i])
    end
    #Monthly inflow
    existinginflow = permutedims(HydroData["existinginflowcf"])
    # Replace NaNs from monthly data with 0's
    replace_nans!(existinginflow)
    # Preallocate array for intermediate result
    avail_inflow = zeros(Float64, T, n)           ## Make it variable

    #Turn monthly inflow data into hourly data
    days = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
    lasthours = 24 * cumsum(days)
    firsthours = [1; 1 .+ lasthours[1:end-1]]

    for m in 1:12
        if lasthours[m] <= T
            for i = firsthours[m]:lasthours[m]
                avail_inflow[i,:] = existinginflow[m,:]
            end
        elseif lasthours[m] >= T
            for i = firsthours[m]:T
                avail_inflow[i,:] = existinginflow[m,:]
            end
        end
    end

    # Get the installed capacity of Run of river and Hydro reservoir plants from ENTSO-E datasets: https://transparency.entsoe.eu/generation/r2/installedGenerationCapacityAggregation/show
    # Create a table (RoR_Res) with values from ENTSOE for Pumped Storage, Run of River and Reservoir for each values
    data_path = joinpath("ENTSO-E_data", "Capacities")
    if sum(occursin.(Regions[1,1], Europe[:,1]))>0
        Countries = Regions[:,1]
        Hydro_Matrix = Array{Any}(undef, 0, 3)
        # take the values from the ENTSO-E CSV files, put them into a matrix
        for i in 1:length(Countries)
            Country = Regions[:,1][i]
            Country_ENTSOE = joinpath(data_path, "$(Country).csv") |> CSV.File |> DataFrame
            RoR_Res_data = Country_ENTSOE[in(["Hydro Pumped Storage", "Hydro Run-of-river and poundage", "Hydro Water Reservoir"]).(Country_ENTSOE."Production Type"), 2]
            replace!(RoR_Res_data, "n/e"=>"0")
            RoR_Res_data = parse.(Float64, string.(RoR_Res_data))
            Water = RoR_Res_data
            Hydro_Matrix = vcat(Hydro_Matrix, permutedims(Water))
        end
        Hydro = hcat(Countries, Hydro_Matrix)
        Header = ["Country" "Pumped Storage" "Run of River" "Reservoir"]
        RoR_Res = vcat(Header, Hydro)
        replace!(RoR_Res, "N/A"=>0)
        # Reservoir percentage is the reservoir value devided by the sum of all three hydro power generation capacities
        reservoir = RoR_Res[2:end,4]./(RoR_Res[2:end,2] .+ RoR_Res[2:end,3] .+ RoR_Res[2:end,4])
        replace_nans!(reservoir)
        reservoir_table = hcat(Regions, reservoir)
        ## Percentage between PHS and reservoir for later calculation of Wmin and Wmax
        reservoir_PHS = RoR_Res[2:end,4]./(RoR_Res[2:end,2] .+ RoR_Res[2:end,4])
        replace_nans!(reservoir_PHS)
    elseif sum(occursin.(Regions[1,1], collect(keys(Regions_dict))))>0
        Hydro_Matrix = Array{Any}(undef, 0, 3)
        ## get the countries, that are in the regions
        Countries = []
        for i in 1:length(Regions)
            R = getindex.(Ref(Regions_dict),(Regions))[i]
            Countries = [Countries; R]
        end
        Hydro_Matrix = Array{Any}(undef, 0, 3)
        ## get ENTSO-E values for the countries in the regions
        for i in 1:length(Countries)
            Country = Countries[:,1][i]
            Country_ENTSOE = joinpath(data_path, "$(Country).csv") |> CSV.File |> DataFrame
            RoR_Res_data = Country_ENTSOE[in(["Hydro Pumped Storage", "Hydro Run-of-river and poundage", "Hydro Water Reservoir"]).(Country_ENTSOE."Production Type"), 2]
            replace!(RoR_Res_data, "n/e"=>"0")
            replace!(RoR_Res_data, "N/A"=>"0")
            RoR_Res_data = parse.(Float64, string.(RoR_Res_data))
            Water = RoR_Res_data
            Hydro_Matrix = vcat(Hydro_Matrix, permutedims(Water))
        end
        Hydro = hcat(Countries, Hydro_Matrix)
        Header = ["Country" "Pumped Storage" "Run of River" "Reservoir"]
        RoR_Res = vcat(Header, Hydro)
        replace!(RoR_Res, "n/e"=>0)
        replace!(RoR_Res, "N/A"=>0)
        RoR_Res[2:end,2:end] = parse.(Float64, string.(RoR_Res[2:end,2:end]))
        res_table = []
        ## get the region of each country
        for i in 1:length(Regions)
            R = getindex.(Ref(Regions_dict),(Regions))[i]
            Key = [k for (k,v) in Regions_dict if v==R]
            Key_List = repeat(Key, inner=length(R))
            res_table = [res_table; Key_List]
        end
        # List all countries with their region
        List_full = DataFrame(hcat(res_table, RoR_Res[2:end,2:end]), :auto)
        # sum the different regions, so you get one value per region
        Hydro_summed = combine(groupby(List_full, :x1), names(List_full, Not(:x1)) .=> sum, renamecols=false)
        reservoir = Hydro_summed[!,:x4]./(Hydro_summed[!,:x2] .+ Hydro_summed[!,:x3] .+ Hydro_summed[!,:x4])
        reservoir_table = hcat(Regions, reservoir)
        ## Percentage between PHS and reservoir for later calculation of Wmin and Wmax
        reservoir_PHS = Hydro_summed[!,:x4]./(Hydro_summed[!,:x2] .+ Hydro_summed[!,:x4])
        replace_nans!(reservoir_PHS)
    end

    reservoirp =  permutedims(reservoir_table[:,2])

    #Hydro capacities
    hydroCap = existingcapac .* reservoirp
    hydroRoRCap = existingcapac - hydroCap
    hydrocapacity = [hydroCap; hydroRoRCap]

    #Reservoir and RoR inflow
    hyd_in = avail_inflow .* hydroCap           ## availability is not what is installed, but the amount of water available
    hydRoR_in = avail_inflow .* hydroRoRCap

# Hydro_capacity.csv file

    ## get the reservoir levels from ENTSO-E, base year 2018: https://transparency.entsoe.eu/generation/r2/waterReservoirsAndHydroStoragePlants/show
    ## Wmax is the maximum value of year 2018 for each country
    ## Wmin is the minimum value of year 2018 for each country
    data_path_reslvl = joinpath("ENTSO-E_data", "Reservoir_level")
    if sum(occursin.(Regions[1,1], Europe[:,1]))>0
        Countries = Regions[:,1]
        Hydro_Matrix = Array{Any}(undef, 0, 3)
        W_Min = Array{Any}(undef, 0, 1)
        W_Max = Array{Any}(undef, 0, 1)
        # take the values from the ENTSO-E CSV files, put them into a matrix
        for i in 1:length(Countries)
            Country = Regions[:,1][i]
            Country_ENTSOE = joinpath(data_path_reslvl, "$(Country).csv") |> CSV.File |> DataFrame
            W_data_total = Matrix(Country_ENTSOE[!, Not("Week")])
            replace!(W_data_total, "n/e"=>"0")
            replace!(W_data_total, "N/A"=>"0")
            W_data_total = parse.(Float64, string.(W_data_total))
            W_Min = vcat(W_Min, minimum(W_data_total, dims=1))
            W_Max = vcat(W_Max, maximum(W_data_total, dims=1))
        end
        Wmin_summed = W_min_total = hcat(Countries, W_Min)
        Wmax_summed = W_max_total = hcat(Countries, W_Max)
    elseif sum(occursin.(Regions[1,1], collect(keys(Regions_dict))))>0
        ## get the countries, that are in the regions
        Countries = []
        for i in 1:length(Regions)
            R = getindex.(Ref(Regions_dict),(Regions))[i]
            Countries = [Countries; R]
        end
        W_Min = Array{Any}(undef, 0, 1)
        W_Max = Array{Any}(undef, 0, 1)
        for i in 1:length(Countries)
            Country = Countries[:,1][i]
            Country_ENTSOE = joinpath(data_path_reslvl, "$(Country).csv") |> CSV.File |> DataFrame
            W_data_total = Matrix(Country_ENTSOE[!, Not("Week")])
            replace!(W_data_total, "n/e"=>"0")
            replace!(W_data_total, "N/A"=>"0")
            W_data_total = parse.(Float64, string.(W_data_total))
            W_Min = vcat(W_Min, minimum(W_data_total, dims=1))
            W_Max = vcat(W_Max, maximum(W_data_total, dims=1))
        end
        W_min_total = hcat(Countries, W_Min)
        W_max_total = hcat(Countries, W_Max)
        W_table = []
        for i in 1:length(Regions)
            R = getindex.(Ref(Regions_dict),(Regions))[i]
            Key = [k for (k,v) in Regions_dict if v==R]
            Key_List = repeat(Key, inner=length(R))
            W_table = [W_table; Key_List]
        end
        Wmin_full = DataFrame(hcat(W_table, W_min_total[:,2]), :auto)
        Wmin_summed = combine(groupby(Wmin_full, :x1), names(Wmin_full, Not(:x1)) .=> sum, renamecols=false)
        Wmax_full = DataFrame(hcat(W_table, W_max_total[:,2]), :auto)
        Wmax_summed = combine(groupby(Wmax_full, :x1), names(Wmax_full, Not(:x1)) .=> sum, renamecols=false)
    end

    H_max = potentialcapac .* permutedims(reservoirp) # reservoir production capacity per node [MWh]
    NodeNr = [1:n;]
    H_min = permutedims(hydroCap)      # existing reservoir capacity per node (percentage)
    HR_max = permutedims(hydroRoRCap)        # RoR counterpart of existing capacity
    ## ENTSO-E data includes both PHS and reservoir, so we need the reservoir percentage to calculate acutal Wmin and Wmax for reservoir only
    W_max = Wmax_summed[:,2].*reservoir_PHS
    W_min = Wmin_summed[:,2].*reservoir_PHS
    hydro_tech = repeat(1:1, inner=n)
    Hydro_cap = DataFrame(hcat(NodeNr, hydro_tech, H_min, H_max, W_min, W_max), :auto)
    rename!(Hydro_cap, ["node", "hydro_tech", "hcap_min", "hcap_max", "wcap_min" ,"wcap_max"])
    CSV.write(joinpath(instance_path, "hydro_capacity.csv"), Hydro_cap)

# hydro.csv file
    node = [1:n;]
    HydroRoR = permutedims(hydroRoRCap)
    Fmin_mat = fill(Fmin, n)
    hyd_in_avg = []
    for i in 1:n
        hyd_in_avg = vcat(hyd_in_avg, mean(hyd_in[:,i]))
    end
    hydror_avg = []
    for i in 1:n
        hydror_avg = vcat(hydror_avg, mean(hydRoR_in[:,i]))
    end
    hyd_flow_min = []
    for i in 1:n
       hyd_flow_min = vcat(hyd_flow_min, argmin((avail_inflow .* existingcapac)[:,i]))
    end
    Hydro = DataFrame(hcat(node, HydroRoR, Fmin_mat, hyd_in_avg, hydror_avg, hyd_flow_min), :auto)
    rename!(Hydro, ["node", "HydroRoR", "Fmin", "hyd_in_avg", "hydror_avg" ,"hyd_flow_min"])
    CSV.write(joinpath(instance_path, "hydro.csv"), Hydro)

# function create_node_data()
    # Create a Vector for every technology spanning all nodes (at 5 nodes: 43.800 columns)
    Demand = []
    WindOff = []
    WindOn = []
    Solar = []
    hydro_in = []
    hydroRoR_in = []

    for i in 1:n
        Demand$i = sum(demand[:,i], dims = 2)
        append!(Demand, Demand$i)
    end
    for i in 1:n
        WindOff$i = sum(avail_wind_off[:,i], dims = 2)
        append!(WindOff, WindOff$i)
    end
    for i in 1:n
        WindOn$i = sum(avail_wind_on[:,i], dims = 2)
        append!(WindOn, WindOn$i)
    end
    for i in 1:n
        Solar$i = sum(avail_sol[:,i], dims = 2)
        append!(Solar, Solar$i)
    end
    for i in 1:n
        hyd_in$i = sum(hyd_in[:,i], dims = 2)
        append!(hydro_in, hyd_in$i)
    end
    for i in 1:n
        hydRoR_in$i = sum(hydRoR_in[:,i], dims = 2)
        append!(hydroRoR_in, hydRoR_in$i)
    end
    All = hcat(Demand, Solar, WindOn, WindOff, hydro_in, hydroRoR_in)

    # Generate CSV files for every single node
    x = 1
    for i = 1:n
        NodeData = All[x:T*i,:]
        x = x+T
        N_m = DataFrame(NodeData, :auto)
        rename!(N_m, ["Demand", "Avail_Sol", "Avail_Wind_On" ,"Avail_Wind_Off", "Hyd_In" ,"HydRoR_In"])
        CSV.write(joinpath(nodes_path, "$i.csv"), N_m)
    end

    x = 1
    # create file for node_specs
    NodeDemand = []
    for i = 1:n
        Max_Node_Demand = maximum(Demand[x:T*i], dims=1)
        x = x+T
        NodeDemand = append!(NodeDemand, Max_Node_Demand)
    end

    NodeNr = [1:n;]
    Max_Demand = NodeDemand
    nodes_specs = DataFrame(hcat(NodeNr, Regionlist, Max_Demand), :auto)
    rename!(nodes_specs, ["Node", "Name", "Max_Demand" ])
    CSV.write(joinpath(instance_path, "nodes_specs.csv"), nodes_specs)

# function gen_capacity()
    # capacity of the remaining non-renewable energy sources fixed at 1000 GW (biomass, nuclear, coal, gas_cc, gas_oc)
    Capacity_else = 1000*1000
    Capacity = zeros(5,1)
    gcap_max = []

    # create a vector with the capacity of each energy source per node
    for i in 1:n
        Capacity$i = [sum(WindCapTotal[i,:]), sum(WindCapOff[i,:]), sum(SolarCapTotal[i,:]),  Capacity_else, Capacity_else, Capacity_else, Capacity_else, Capacity_else].*1000
        append!(gcap_max, Capacity$i)
    end

    # Zero matrix with 4 columns
    gen_capacity = zeros(n*t, 4)

    # Column with node numbers from 1 to node n, which repeats itself as often as we have technologies
    node = repeat(1:n, inner=t)
    # repeating the technologies 1 - 8 for every node
    gen_tech = repeat([1:t;], n)
    # fill the whole column of gcap_min with 0
    gcap_min = fill(0, (n*t,1))
    Gen_capac = hcat(node, gen_tech, gcap_min, gcap_max)
    gen_capacity = DataFrame(Gen_capac, :auto)
    rename!(gen_capacity, ["node", "gen_tech", "gcap_min" ,"gcap_max"])
    CSV.write(joinpath(instance_path, "gen_capacity.csv"), gen_capacity)

# sto_capacity()
    st = 1
    Max_capacity = 1000000000
    s = repeat([1:st;], n)
    node = repeat(1:n, inner=st)
    scap_min = zeros(n)
    scap_max = fill(Max_capacity, (n*st, 1))
    sto_cap = DataFrame(hcat(s, node, scap_min, scap_max), :auto)
    rename!(sto_cap, ["s", "node", "scap_min", "scap_max"])
    CSV.write(joinpath(instance_path, "sto_capacity.csv"), sto_cap)

# transmission()
    Dist = Distances[tril!(trues(size(Distances)), -1)]
    # Connected nodes
    Con = connected[tril!(trues(size(connected)), -1)]
    # Find distances for connected nodes only
    Dist_Con = hcat(Dist, Con)
    Dist_Con_df = DataFrame(Dist_Con, :auto)
    Dist_con = Dist_Con_df[Dist_Con_df[!,:x2].!=0,:]
    Distance_Connected = Matrix(Dist_con)[:,1]
    # Ged rid of double nodes by only using lower triangular matrix
    LinesLT = LowerTriangular(connected)
    # Find connections by getting the CartesianIndex of the connected nodes
    LinesCI = findall(x->x==1, LinesLT)
    Lines = hcat(getindex.(LinesCI, 1), getindex.(LinesCI, 2))
    LinesDoubledf = DataFrame(Lines, :auto)
    # Delete entries of connections between only one node
    Linesdf = LinesDoubledf[LinesDoubledf.x1 .!= LinesDoubledf.x2, :]
    LinesMatrix = Matrix(Linesdf)
    # Reshape into a 1d array
    LinesSingle = collect(eachrow(LinesMatrix))
    L_ind = [1:size(LinesMatrix)[1];]

    # Create Matrix for remaining parameters of transmission
    TransData = [
       # :cost       :converter_cost         :M      :C      :B      :efficiency         :lifetime       :tcap_min       :tcap_max
        400          150000                  0.02    0       1       0.95                40              0               1000000000
    ]

    # Fill the matrix according to the number of pairs in the system
    TransDatafull = repeat(TransData; outer=[size(LinesMatrix)[1]])
    # combine the tables, add header and create CSV file
    TransmissionData = DataFrame(hcat(L_ind, LinesSingle, Distance_Connected, TransDatafull), :auto)
    rename!(TransmissionData, ["L_ind" ,"line" ,"dist" , "cost", "converter_cost", "M" ,"C" , "B", "efficiency", "lifetime" ,"tcap_min" ,"tcap_max"])
    CSV.write(joinpath(instance_path, "transmission.csv"), TransmissionData)

# function generate_tech_tables()
    techtable = [
        :name          :g       :investment_cost    :fixedOM    :varOM      :fuel_cost      :efficiency     :emissions      :lifetime       :r_minus    :r_plus
        :wind_on        1       1473                35          0           0               1               0               25              1           1
        :wind_off       2       3800                80          0           0               1               0               25              1           1
        :solar          3       995                 26          0           0               1               0               25              1           1
        :biomass        4       2141                100         0           7               0.448           0.39            30              1           1
        :nuclear        5       5000                150         3           3               0.34            0               60              0.05        0.05
        :coal           6       1800                25          6           8               0.466           0.34            40              0.15        0.15
        :gas_cc         7       800                 20          4           22              0.615           0.20            30              0.3         0.3
        :gas_oc         8       400                 15          3           22              0.395           0.20            30              1           1
    ]
    CSV.write(joinpath(instance_path, "gen_technology.csv"), Tables.table(techtable), writeheader=false)

    ## sources investment cost
    # Wind onshore, Wind offshore, Solar, Biomass: https://www.irena.org/publications/2020/Jun/Renewable-Power-Costs-in-2019
    # Nuclear, coal, gas_cc: https://www.econstor.eu/bitstream/10419/80348/1/757528015.pdf (rough average of the different estimations)
    #

    # Hydro
    hydro_tech = [
        :name          :hydro_tech       :investment_cost    :fixedOM    :varOM      :fuel_cost      :efficiency     :emissions      :lifetime       :r_minus    :r_plus
        :hydropower     1                 1704                200          0           0               1               0               60              1           1
    ]
    CSV.write(joinpath(instance_path, "hydro_technology.csv"), Tables.table(hydro_tech), writeheader=false)
    ## sources investment cost
    # https://www.irena.org/publications/2020/Jun/Renewable-Power-Costs-in-2019

    # Storage
    storage = [
        :s          :xi      :C    :cost    :lifetime
        1           0.86     15     185     20
    ]
    CSV.write(joinpath(instance_path, "storage.csv"), Tables.table(storage), writeheader=false)

    # create weights.csv if hours = 8760
    if T == 8760
        weights = repeat([1:1;], 8760)
        df = DataFrame(col1 = weights)
        rename!(df, ["Weights"])
        CSV.write(joinpath(instance_path, "weights.csv"), df)
    end

    # indices.json
    G = [1:t;]
    G_r = [1:4;]
    N = [1:n;]
    L = LinesSingle
    L_ind = L_ind
    S = [1]
    H = [1]
    T = T
    indices = Dict("G" => G, "G_r" => G_r, "N" => N, "L" => L, "L_ind" => L_ind, "S" => S, "H" => H, "T" => T)
    open(joinpath(instance_path, "indices.json"), "w") do f
        JSON.print(f, indices)
    end

end


text:
# scan into a sorted list of free buffers, splitting buffers along the way

using DataStructures

## tunables

# pool boundaries
const SMALL_CUTOFF   = 2^20    # 1 MiB
const LARGE_CUTOFF   = 2^25    # 32 MiB

# memory size rounding
const SMALL_ROUNDOFF = 2^9     # 512 bytes
const LARGE_ROUNDOFF = 2^17    # 128 KiB
const HUGE_ROUNDOFF  = 2^20    # 1 MiB

# maximum overhead (unused space) when selecting a buffer
# small and large buffers with large overhead will be split, but huge buffers never are
const SMALL_OVERHEAD = typemax(Int)
const LARGE_OVERHEAD = typemax(Int)
const HUGE_OVERHEAD  = 0

## block utilities

# split a block at size `sz`, returning the newly created block
function split!(block, sz)
    @assert sz < block.sz "Cannot split a $block at too-large offset $sz"
    split = Block(block.buf, sizeof(block) - sz; off = block.off + sz)
    block.sz = sz

    # update links
    split.prev = block
    split.next = block.next
    if block.next !== nothing
        block.next.prev = split
    end
    block.next = split

    return split
end

# merge a sequence of blocks that starts at `head` and ends with `tail` (inclusive)
function merge!(head, tail)
    head.next = tail.next
    if tail.next !== nothing
        tail.next.prev = head
    end
    head.sz = tail.sz + (tail.off - head.off)

    return head
end

## pooling

const pool_lock = ReentrantLock()

function pool_scan(dev, pool, sz, max_overhead=typemax(Int))
    max_sz = Base.max(sz + max_overhead, max_overhead)   # protect against overflow
    @lock pool_lock begin
        # get the first entry that is sufficiently large
        i = searchsortedfirst(pool, Block(Mem.DeviceBuffer(CU_NULL, 0), sz))
        if i != pastendsemitoken(pool)
            block = deref((pool,i))
            @assert sizeof(block) >= sz
            if sz <= max_sz
                delete!((pool,i))   # FIXME: this allocates
                return block
            end
        end

        return nothing
    end
end

# compact sequences of blocks into larger ones.
# looks up possible sequences based on each block in the input set.
# destroys the input set.
# returns the net difference in amount of blocks.
function pool_compact(dev, blocks)
    compacted = 0
    @lock pool_lock begin
        while !isempty(blocks)
            block = pop!(blocks)
            @assert block.state == AVAILABLE
            pool = get_pool(dev, sizeof(block))
            @assert in(block, pool)

            # find the head of a sequence
            head = block
            while head.prev !== nothing && head.prev.state == AVAILABLE
                head = head.prev
                @assert in(head, pool)
            end
            szclass = size_class(sizeof(head))

            if head.next !== nothing && head.next.state == AVAILABLE
                delete!(pool, head)

                # find the tail (from the head, removing blocks as we go)
                tail = head.next
                while true
                    @assert szclass === size_class(sizeof(tail)) "block $tail should not have been split to a different pool than $head"
                    delete!(pool, tail)    # FIXME: allocates
                    delete!(blocks, tail)
                    tail.state = INVALID
                    compacted += 1
                    if tail.next !== nothing && tail.next.state == AVAILABLE
                        tail = tail.next
                    else
                        break
                    end
                end

                # compact
                head = merge!(head, tail)
                @assert !in(head, pool) "$head should not be in the pool"
                @assert head.state == AVAILABLE
                @assert szclass === size_class(sizeof(head)) "compacted $head should not end up in a different pool"
                push!(pool, head)
            end
        end
    end
    return compacted
end

function pool_reclaim_single(dev, pool, sz=typemax(Int))
    freed = 0

    @lock pool_lock begin
        # mark non-split blocks
        candidates = Block[]
        for block in pool
            if iswhole(block)
                push!(candidates, block)
            end
        end

        # free them
        for block in candidates
            delete!(pool, block)
            freed += sizeof(block)
            actual_free(dev, block)
            freed >= sz && break
        end
    end

    return freed
end

# repopulate the pools from the list of freed blocks
function pool_repopulate(dev)
    blocks = @safe_lock freed_lock begin
        isempty(freed[dev]) && return
        blocks = Set(freed[dev])
        empty!(freed[dev])
        blocks
    end

    @lock pool_lock begin
        for block in blocks
            pool = get_pool(dev, sizeof(block))
            @assert !in(block, pool) "$block should not be in the pool"
            @assert block.state == FREED "$block should have been marked freed"
            block.state = AVAILABLE
            push!(pool, block) # FIXME: allocates
        end

        pool_compact(dev, blocks)
    end

    return
end

const SMALL = 1
const LARGE = 2
const HUGE  = 3

# sorted containers need unique keys, which the size of a block isn't.
# mix in the block address to keep the key sortable, but unique.
unique_sizeof(block::Block) = (UInt128(sizeof(block))<<64) | UInt64(pointer(block))
const UniqueIncreasingSize = Base.By(unique_sizeof)

const pool_small = PerDevice{SortedSet{Block}}((dev)->SortedSet{Block}(UniqueIncreasingSize))
const pool_large = PerDevice{SortedSet{Block}}((dev)->SortedSet{Block}(UniqueIncreasingSize))
const pool_huge  = PerDevice{SortedSet{Block}}((dev)->SortedSet{Block}(UniqueIncreasingSize))

const freed = PerDevice{Vector{Block}}((dev)->Vector{Block}())
const freed_lock = NonReentrantLock()

function size_class(sz)
    if sz <= SMALL_CUTOFF
        SMALL
    elseif SMALL_CUTOFF < sz <= LARGE_CUTOFF
        LARGE
    else
        HUGE
    end
end

@inline function get_pool(dev, sz)
    szclass = size_class(sz)
    if szclass == SMALL
        return pool_small[dev]
    elseif szclass == LARGE
        return pool_large[dev]
    elseif szclass == HUGE
        return pool_huge[dev]
    end
end

function pool_alloc(dev, sz)
    szclass = size_class(sz)

    # round off the block size
    req_sz = sz
    roundoff = if szclass == SMALL
        SMALL_ROUNDOFF
    elseif szclass == LARGE
        LARGE_ROUNDOFF
    elseif szclass == HUGE
        HUGE_ROUNDOFF
    end
    sz = cld(sz, roundoff) * roundoff
    szclass = size_class(sz)

    # select a pool
    pool = get_pool(dev, sz)

    # determine the maximum scan overhead
    max_overhead = if szclass == SMALL
        SMALL_OVERHEAD
    elseif szclass == LARGE
        LARGE_OVERHEAD
    elseif szclass == HUGE
        HUGE_OVERHEAD
    end

    block = nothing
    for phase in 1:3
        if phase == 2
            @pool_timeit "$phase.0 gc (incremental)" GC.gc(false)
        elseif phase == 3
            @pool_timeit "$phase.0 gc (full)" GC.gc(true)
        end

        @pool_timeit "$phase.1 repopulate" pool_repopulate(dev)

        @pool_timeit "$phase.2 scan" begin
            block = pool_scan(dev, pool, sz, max_overhead)
        end
        block === nothing || break

        @pool_timeit "$phase.3 alloc" begin
            block = actual_alloc(dev, sz)
        end
        block === nothing || break

        # we're out of memory, try freeing up some memory. this is a fairly expensive
        # operation, so start with the largest pool that is likely to free up much memory
        # without requiring many calls to free.
        for pool in (pool_huge[dev], pool_large[dev], pool_small[dev])
            @pool_timeit "$phase.4a reclaim" pool_reclaim_single(dev, pool, sz)
            @pool_timeit "$phase.4b alloc" block = actual_alloc(dev, sz)
            block === nothing || break
        end
        block === nothing || break

        # last-ditch effort, reclaim everything
        @pool_timeit "$phase.5a reclaim" begin
            pool_reclaim_single(dev, pool_huge[dev])
            pool_reclaim_single(dev, pool_large[dev])
            pool_reclaim_single(dev, pool_small[dev])
        end
        @pool_timeit "$phase.5b alloc" block = actual_alloc(dev, sz)
    end

    if block !== nothing
        @assert size_class(sizeof(block)) == szclass "Cannot satisfy an allocation in pool $szclass with a buffer in pool $(size_class(sizeof(block)))"

        # split the block if that would yield one from the same pool
        # (i.e. don't split off chunks that would be way smaller)
        # TODO: avoid creating unaligned blocks here (doesn't currently seem to happen
        #       because of the roundoff, but we should take care anyway)
        remainder = sizeof(block) - sz
        if szclass != HUGE && remainder > 0 && size_class(remainder) == szclass
            split = split!(block, sz)
            split.state = AVAILABLE
            @lock pool_lock begin
                push!(pool, split)
            end
        end
    end

    if block !== nothing
        block.state = ALLOCATED
    end
    return block
end

function pool_free(dev, block)
    # we don't do any work here to reduce pressure on the GC (spending time in finalizers)
    # and to simplify locking (preventing concurrent access during GC interventions)
    block.state == ALLOCATED || error("Cannot free a $(block.state) block")
    block.state = FREED
    @safe_lock_spin freed_lock begin
        push!(freed[dev], block)
    end
end

function pool_init()
    initialize!(freed, ndevices())

    initialize!(pool_small, ndevices())
    initialize!(pool_large, ndevices())
    initialize!(pool_huge, ndevices())
end

function pool_reclaim(dev, sz::Int=typemax(Int))
    pool_repopulate(dev)

    freed_sz = 0
    for pool in (pool_huge[dev], pool_large[dev], pool_small[dev])
        freed_sz >= sz && break
        freed_sz += pool_reclaim_single(dev, pool, sz-freed_sz)
    end
    return freed_sz
end

function cached_memory(dev=device())
    sz = @safe_lock freed_lock mapreduce(sizeof, +, freed[dev]; init=0)
    @lock pool_lock for pool in (pool_small[dev], pool_large[dev], pool_huge[dev])
        sz += mapreduce(sizeof, +, pool; init=0)
    end
    return sz
end

ext: jl


text:
using SafeTestsets

@safetestset "Test Function Bounds" begin
    include("test_bounds.jl")
end

@safetestset "Test Utilities" begin
    include("test_utilities.jl")
end

@safetestset "Test Quadrature Structs" begin
    include("test_quadrature_structs.jl")
end

@safetestset "Test One-Dimensional Quadrature" begin
    include("test_one_dimensional_quadrature.jl")
end

@safetestset "Test Area Quadrature" begin
    include("test_area_quadrature.jl")
end

@safetestset "Test Surface Quadrature" begin
    include("test_surface_quadrature.jl")
end

@safetestset "Test Subdivision Algorithm" begin
    include("test_subdivision.jl")
end

max_stars_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_issues_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_forks_repo_name: ArjunNarayanan/ImplicitDomainQuadrature


text:
module CEDICT

export DictionaryEntry, ChineseDictionary,

search_headwords, search_senses, search_pinyin,

idioms

include("dictionary.jl")
include("searching.jl")

"""
    idioms([dict])

Retrieves the set of idioms in the provided dictionary (by looking for a label of "(idiom)" in any
of the senses) or in the default dictionary if none provided.
"""
function idioms end
idioms(dict) = search_senses(dict, "(idiom)")
idioms() = idioms(ChineseDictionary())

end


text:
using Erdos
using ErdosExtras
using Test
using SparseArrays: sprand
using LinearAlgebra: Diagonal
using Random

testdir = dirname(@__FILE__)

tests = [
    "tsp/TSP",
    "matching/matching",
    "matching/b-matching",
    ]

GLIST =    [
            (Graph{Int64}, DiGraph{Int64}),
            (Graph{UInt32}, DiGraph{UInt32}),
            (Network, DiNetwork)
            ]

for GDG in GLIST
    global G = GDG[1]
    global DG = GDG[2]
    global E = edgetype(G)
    global V = vertextype(G)

    for t in tests
        include(joinpath(@__DIR__,"$(t).jl"))
    end
end


text:
println("LATHE.MODELS TESTS")
println("......................")
using Lathe.models: LinearLeastSquare, LinearRegression
x = [5, 10, 15, 20]
y = [5, 10, 15, 20]
@testset "Regressors" begin
m = LinearRegression(x, y)
testy = [5, 10, 15]
yhat = m.predict(testy)
@test yhat == testy
m = LinearLeastSquare(x, y)
yhat = m.predict(testy)
@test yhat == testy
end


text:
using MovingWeightedLeastSquares

ins = [1. 1.; 1. -1.0; -1.0 1.0; -1.0 -1.0; 0. 0.; 1. 0.; -1. 0.; 0. 1.; 0. -1]
outs = [1., -0.5, 1., 1., -1., 0., 0., 0., 0.]

m = mwls(ins, outs, 1., (d, e) -> exp(-d^2))

println(MovingWeightedLeastSquares.diff(m, [0., 0.], (1, 1); dist = 10.))


text:
export ATARI_GAMES
export atari_init

const atari_checksum = "e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855"

const ATARI_GAMES = [
    "air-raid", "alien", "amidar", "assault", "asterix",
    "asteroids", "atlantis", "bank-heist", "battle-zone", "beam-rider",
    "berzerk", "bowling", "boxing", "breakout", "carnival", "centipede",
    "chopper-command", "crazy-climber", "demon-attack",
    "double-dunk", "elevator-action", "enduro", "fishing-derby", "freeway",
    "frostbite", "gopher", "gravitar", "hero", "ice-hockey", "jamesbond",
    "journey-escape", "kangaroo", "krull", "kung-fu-master",
    "montezuma-revenge", "ms-pacman", "name-this-game", "phoenix",
    "pitfall", "pong", "pooyan", "private-eye", "qbert", "riverraid",
    "road-runner", "robotank", "seaquest", "skiing", "solaris",
    "space-invaders", "star-gunner", "tennis", "time-pilot", "tutankham",
    "up-n-down", "venture", "video-pinball", "wizard-of-wor",
    "yars-revenge", "zaxxon"
]

function fetch_atari_ds(src, dest)
    try run(`which gsutil`) catch x throw("gsutil not found, install gsutil to proceed further") end

    run(`gsutil -m cp -r $src $dest`)
    return dest
end

game_name(game) = join(titlecase.(split(game, "-")))
# use function to initialise atari
function atari_init()
    for game in ATARI_GAMES
        for index in 1:5
            register(
                DataDep(
                    "atari-replay-datasets-$game-$index",
                    """
                    Dataset: The DQN Replay Dataset from Google Research
                    Authors: Rishabh Agarwal, Dale Schuurmans, Mohammad Norouzi
                    Publication Year: 2020
                    Related Publications: https://research.google/pubs/pub49020/
                    Credits: https://arxiv.org/abs/1907.04543
                    Journal: International Conference on Machine Learning

                    The DQN Replay Dataset is generated using DQN agents trained on 60 Atari 2600
                    games for 200 million frames each, while using sticky actions (with 25% probability
                    that the agent’s previous action is executed instead of the current action) to make
                    the problem more challenging. For each of the 60 games, we train 5 DQN agents with different
                    random initializations, and store all of the (state, action, reward, next state) tuples
                    encountered during training into 5 replay datasets per game, resulting in a total of 300 datasets.
                    """,
                    "gs://atari-replay-datasets/dqn/$(game_name(game))/$index/replay_logs/";
                    fetch_method = fetch_atari_ds
                )
            )
        end
    end
end


text:
using HOODESolver
using Test

include("henon_heiles.jl")

function testprepare()
    n_tau=32
    A = [0 0 1 0; 0 0 0 0;-1 0 0 0; 0 0 0 0]
    u =rand(BigFloat, 4)
    for i=1:5
        prec = 256*2^(i-1)
        setprecision(256*2^(i-1))
        for i_eps = 1:40
            epsilon = big"1.25"^(-i_eps)
            println("prec=$prec epsilon=$epsilon")
            for order=3:10
                @time @testset "test PrepareU0 prec=$prec order=$order epsilon=$epsilon" begin
                    newprec = prec + 32 + div(-exponent(epsilon)*order^2, 3)
                    precref = newprec + 10
                    parphi = PreparePhi(epsilon, n_tau, A, henon_heiles)
                    pref=PrepareU0(parphi, order, u, precref)
                    pnew=PrepareU0(parphi, order, u, newprec)
                    @test pref.ut0 == pnew.ut0
                end
            end
        end
    end
end
testprepare()


text:
#=
dynamic_nelson_siegel_model.jl

    Dynamic Nelson Siegel model specification in state space formulation

@author: Quint Wiersma <q.wiersma@vu.nl>

@date: 2022/03/28
=#

"""
    DynamicNelsonSiegelModel <: StateSpaceModel

Constructor for a dynamic Nelson Siegel model instance of the state space model type
with hyperparameters `λ`, `ϕ`, and error specifications
"""
mutable struct DynamicNelsonSiegelModel{Ty, Tτ, Tλ, Tϕ, Teobs, Tefac} <: StateSpaceModel
    y::Ty               # data
    τ::Tτ               # maturities
    λ::Tλ               # decay
    ϕ::Tϕ               # autoregressive parameters
    error_obs::Teobs    # obs. eq. error specification
    error_factor::Tefac # factor eq. error specification
end
# Constructors
function DynamicNelsonSiegelModel(y::AbstractMatrix, τ::AbstractVector,
                                    error_obs::AbstractErrorModel, error_factor::AbstractErrorModel)
    # hyper paremeters
    λ= .0609
    ϕ= similar(y, 3, 3)

    return DynamicNelsonSiegelModel(y, τ, λ, ϕ, error_obs, error_factor)
end

# Methods
cov(model::DynamicNelsonSiegelModel)= cov(model.error_obs)
prec(model::DynamicNelsonSiegelModel)= prec(model.error_obs)
resid(model::DynamicNelsonSiegelModel)= resid(model.error_obs)
maturities(model::DynamicNelsonSiegelModel)= model.τ
decay(model::DynamicNelsonSiegelModel)= model.λ

function loadings!(Λ::AbstractMatrix, model::DynamicNelsonSiegelModel)
    λ= decay(model)
    τ= maturities(model)
    n= length(τ)

    Λ.= one(eltype(Λ))
    @inbounds @fastmath for i in 1:n
        Λ[i,2]= (1 - exp(-λ * τ[i])) * inv(λ * τ[i])
        Λ[i,3]= Λ[i,2] - exp(-λ * τ[i])
    end

    return nothing
end
function loadings(model::DynamicNelsonSiegelModel)
    # get dims
    n= length(maturities(model))

    Λ= Matrix{typeof(decay(model))}(undef, n, 3)
    loadings!(Λ, model)

    return Λ
end

# State space system and hyperparameters
nparams(model::DynamicNelsonSiegelModel)= 10 + nparams(model.error_obs) + nparams(model.error_factor)

function get_params!(ψ::AbstractVector, model::DynamicNelsonSiegelModel)
    idx= 0  # index counter

    # Store values
    # λ
    ψ[1]= model.λ
    idx+= 1

    # ϕ
    ψ[idx+1:idx+9].= vec(model.ϕ)
    idx+= 9

    # observation equation error model
    get_params!(view(ψ,idx+1:idx+nparams(model.error_obs)), model.error_obs)
    idx+= nparams(model.error_obs)

    # factor equation error model
    get_params!(view(ψ,idx+1:idx+nparams(model.error_fac)), model.error_factor)

    return nothing
end

function get_system!(sys::LinearTimeInvariant, model::DynamicNelsonSiegelModel, method::Symbol)
    Λ= loadings(model)  # loadings

    # Store values
    sys.T.= model.ϕ
    sys.Q.= cov(model.error_factor)
    if method === :univariate
        # Cholesky decomposition of H
        C= cholesky(cov(model))

        # Store values
        ldiv!(sys.y, C.L, model.y)
        ldiv!(sys.Z, C.L, Λ)
    elseif method === :collapsed
        # Collapsing transformation
        # Z'×H⁻¹×Z
        tmp_pn= transpose(Λ) * prec(model)
        tmp= tmp_pn * Λ
        # (Z'×H⁻¹×Z)⁻¹
        # perform pseudo inverse as Λ can contain zero columns
        pseudo= pinv(tmp)
        # pivoted Cholesky decomposition
        C= cholesky!(Hermitian(pseudo), Val(true), check=false)
        # transformation
        ip= invperm(C.p)
        U= C.U[:,ip]
        A= U * tmp_pn

        # Store values
        mul!(sys.y, A, model.y)
        sys.Z.= pinv(transpose(U))
    else
        # Store values
        sys.Z.= Λ
        sys.H.= cov(model)
    end

    return nothing
end

function store_params!(model::DynamicNelsonSiegelModel, ψ::AbstractVector)
    idx= 0  # index counter

    # Store values
    # λ
    model.λ= ψ[1]
    idx+= 1

    # ϕ
    vec(model.ϕ).= view(ψ,idx+1:idx+9)
    idx+= 9

    # observation equation error model
    store_params!(model.error_obs, view(ψ,idx+1:idx+nparams(model.error_obs)))
    idx+= nparams(model.error_obs)

    # factor equation error model
    store_params!(model.error_factor, view(ψ,idx+1:idx+nparams(model.error_factor)))

    return nothing
end

# Initialization
function init!(model::DynamicNelsonSiegelModel, init::NamedTuple, method::Symbol)
    # Model
    init_model!(model, init)

    # State space system
    create_system(model, method)
    fix_system!(sys, model, method)
    init_system!(sys, model)

    return sys
end

function init_factors(y::AbstractMatrix, Λ::AbstractMatrix)
    # Get dims
    T= size(y,2)

    # Result container
    β= similar(y, 3, T)

    # OLS
    A= transpose(Λ) * Λ
	B= A \ transpose(Λ)
    mul!(β, B, y)

    return β
end

"""
    init_ϕ!(ϕ, β)

Initialize autoregressive parameters `ϕ` using OLS with factors, storing the
result in `ϕ`.

#### Arguments
  - `β::AbstractMatrix` : factors

#### Returns
  - `ϕ::AbstractMatrix` : autoregressive parameters
"""
function init_ϕ!(ϕ::AbstractMatrix, β::AbstractMatrix)
    # Get dims
    T= size(β,2)

    # OLS
    X= view(β,:,1:T-1)
    y= view(β,:,2:T)
    ϕ.= y / X

    return nothing
end

function init_model!(model::DynamicNelsonSiegelModel, init::NamedTuple)
    # get dims
    T= size(model.y,2)

    # λ
    haskey(init, :λ) ? model.λ= init.λ : nothing
    # factors
    Λ= loadings(model)
    β= init_factors(model.y, Λ)

    # observation equation residuals
    resid(model).= model.y
    mul!(resid(model), Λ, β, -1., 1.)

    # ϕ
    haskey(init, :ϕ) ? model.ϕ.= init.ϕ : init_ϕ!(model.ϕ, β)

    # factor equation residuals
    resid(model.error_factor).= view(β,:,2:T)
    mul!(resid(model.error_factor), model.ϕ, view(β,:,1:T-1), -1., 1.)

    # observation equation error model
    init_error!(model.error_obs, init)

    # factor equation error model
    init_error!(model.error_factor, init)

    return nothing
end

function create_system(model::DynamicNelsonSiegelModel, method::Symbol)
    # get dims
    (n,T_len)= method === :collapsed ? (3, size(model.y,2)) : size(model.y)

    # types
    Te= eltype(model.y)
    Tv= Vector{Te}
    Tm= Matrix{Te}
    Td= Diagonal{Te}
    TQ= cov(model.error_factor) isa Symmetric ? Symmetric{Te} : Td
    if method === :univariate || method === :collapsed
        TH= Td
    else
        TH= cov(model) isa Symmetric ? Symmetric{Te} : Td
    end

    # State space system
    sys= LinearTimeInvariant{Tm, Tm, Tm, Tv, Tv, TH, TQ, Tv, Tm}(n, 3, T_len)

    return sys
end

function fix_system!(sys::LinearTimeInvariant, model::DynamicNelsonSiegelModel, method::Symbol)
    # infer type
    T= eltype(model.y)

    if method === :multivariate || method === :woodbury
        sys.y.= model.y
    end
    if method === :univariate || method === :collapsed
        sys.H.diag.= one(T)
    end
    sys.d.= zero(T)
    sys.c.= zero(T)
end

function init_system!(sys::StateSpaceSystem, model::DynamicNelsonSiegelModel)
    T= eltype(model.y)
    # a
    sys.a1.= zero(T)
    # P
    sys.P1.= zero(T)
    @inbounds @fastmath for i in 1:3
        sys.P1[i,i]= one(T)
    end

    return nothing
end

# Estimation
# Log-likelihood
function loglik(filter::KalmanFilter, model::DynamicNelsonSiegelModel, method::Symbol)
    # Compute log-likelihood
    ll= loglik(filter)

    # Account for transformation
    if method === :univariate
        # get dims
        T_len= size(model.y,2)

        # Cholesky decomposition of H
        C= cholesky(cov(model))

        # Add Jacobian determinant term
        d,s= logabsdet(C.L)
        ll-= T_len * (d + log(s))
    elseif method === :collapsed
        # get dims
        (n,T_len)= size(model.y)

        Λ= loadings(model)  # loadings

        # Cholesky decomposition of H
        C= cholesky(cov(model))

        # Add Jacobian determinant term
        ll-= .5 * T_len * logdet(C)

        # Z'×H⁻¹×Z
        tmp_pn= transpose(Λ) * prec(model)
        tmp= tmp_pn * Λ
        # (Z'×H⁻¹×Z)⁻¹
        # perform pseudo inverse as Λ can contain zero columns
        pseudo= pinv(tmp)
        # Projection
        M= I - Λ * pseudo * tmp_pn

        # Add projected out term
        ll-= .5 * (n - model.r) * T_len * log(2*π)
        e= similar(model.y, n)
        @inbounds @fastmath for t in 1:T_len
            mul!(e, M, view(model.y,:,t))
            ll-= .5 * dot(e, prec(model), e)
        end
    end

    return ll
end

# Forecast
function reinstantiate(model::DynamicNelsonSiegelModel, h::Integer)
    # number of time series
    n= size(model.y,1)

    # expand
    y_f= hcat(model.y, fill(NaN, n, h))

    return DynamicNelsonSiegelModel(y_f, model.τ, model.λ, model.ϕ, model.error_obs, model.error_factor)
end

function forecast(model::DynamicNelsonSiegelModel, h::Integer)
    # reinstantiate
    forecast_model= reinstantiate(model, h)

    # State space system
    sys= create_system(forecast_model, :multivariate)
    fix_system!(sys, forecast_model, :multivariate)
    init_system!(sys, forecast_model)
    get_system!(sys, forecast_model, :multivariate)

    # forecasts
    f= forecast(sys, h)

    return f
end


text:
ei(::Type{T}, i::Integer, n::Integer) where T = (v = zeros(T, n); v[i] = one(T); v)
ei(i::Integer, n::Integer) = ei(Float64, i, n)

function pivot!(M::Array{T,2}, i::Int, j::Int) where T
    m,n = size(M)
    @assert M[i, j] != 0 "pivot  illegal"
    M[i, :] = M[i, :] / M[i, j]
    for k in setdiff(1:m, i)
        M[k, :] -= M[k, j] * M[i, :]
    end
    return M
end


text:
export UpdateRevert

"""
    UpdateRevert

 Revert a set of ratings to a baseline rating using weighted average
    -- repeated used could be used to do EWMA smoother

## Parameters
* `r_base::Float64`: baseline rating to revert towards
* `α::Float64`:      weight to give to existing rating, (1-α) is weight towards baseline r_base

```
"""
struct UpdateRevert <: UpdateRule
    r_base::Float64
    α::Float64
    function UpdateRevert(r_base, α)
        @check_args(UpdateElo, α >= zero(α))
        @check_args(UpdateElo, α <= one(α))
        new(r_base, α)
    end
end
UpdateRevert(;r_base=1500.0, α=0.75) = UpdateRevert(r_base, α)

function update_info( rule::UpdateRevert )
    info = Dict(
                :name => "Revert",
                :reference => "Patterned after 538's approach to dealing with season changes: see https://fivethirtyeight.com/features/how-we-calculate-nba-elo-ratings/",
                :computation => "simultaneous",# could be argued this is sequential
                :state_model => "recursive",
                :input => "none",
                :output => "either",
                :model => "single",
                :ties => true,
                :factors => false,
                :parameters => ["r_base", "α"],
                :record => true
                )
    return info
end
update_info( ::Type{UpdateRevert} ) =  update_info( UpdateRevert() )

function update_ratings( rule::UpdateRevert,
                         input_ratings::RatingsList,
                         input_competitions::DataFrame)
    m = length( input_ratings.players )
    I = player_indexes( input_ratings.players )

    r = copy( input_ratings.ratings )
    for p in input_ratings.players
        if !haskey( r, p )# set initial ratings if they aren't already
            r[p] = rule.r_base
        else
            r[p] = (1-rule.α) * rule.r_base  +   rule.α * r[p]
        end
    end

    # output ratings list
    output_ratings = RatingsList(input_ratings.players, r )
    return output_ratings
end

function update_ratings( rule::UpdateRevert,
                         input_ratings::RatingsList )
    output_ratings = update_ratings( rule, input_ratings, DataFrame() )
end


text:
mutable struct ResidueRingPolyMap{D, C, T} <: Map{D, C, HeckeMap, ResidueRingPolyMap}
  header::MapHeader{D, C}
  gen_image::Generic.Res{T}
  coeff_map::Map # can be missing if domain and codomain have the same
                 # base_ring(base_ring())
  function ResidueRingPolyMap{D, C, T}(domain::D, codomain::C, gen_image::Generic.Res{T}, coeff_map::Map) where {D, C, T}
    z = new{D, C, T}()
    z.gen_image = gen_image
    z.coeff_map = coeff_map

    image = function(a::Generic.Res)
      #a should be in the domain of M...
      I = codomain(0)
      for i in degree(a.data):-1:0
        I = I*z.gen_image + z.coeff_map(coeff(a.data, i))
      end
      return I::elem_type(C)
    end

    # I need to call preimage in _preimage
    _preimage = function(a::Generic.Res)
      R = codomain
      parent(a) == R || throw("mixed rings in preimage")
      g = gens(domain)
      im_gen = map(z, g) # apply x -> z(x) to the generatotrs ## possibly should be cached and stored
          ## need to write the elements in a matrix, solve the eqn for a
      Mt = zero_matrix(base_ring(base_ring(R)), degree(R.modulus), degree(R.modulus))

      for i=1:degree(R.modulus)
        elem_to_mat_row!(Mt, i, im_gen[i])
      end
      b = zero_matrix(base_ring(base_ring(R)), 1, degree(R.modulus))
      elem_to_mat_row!(b, 1, a)
      s = solve_rational(Mt', b') # why, oh why is solve operating on columns????
      # This is the worst
      if isa(s, Tuple) ## again, why, oh why is solve doing things differently
                   ## over rings than fields?
        s = s[1] * inv(s[2]) # all rings here (type) are actually fields (math)
      end

      s = sum([preimage(z.coeff_map, s[i,1])*g[i] for i=1:length(g)])
      @assert parent(s) == domain
      return s::elem_type(domain)
    end

    z.header = MapHeader(domain, codomain, image, _preimage)
    return z
  end

  function ResidueRingPolyMap{D, C, T}(domain::D, codomain::C, gen_image::Generic.Res{T}) where {D, C, T}
    z = new{D, C, T}()
    z.gen_image = gen_image

    image = function(a::Generic.Res)
      I = z.codomain(0)
      for i in degree(a.data):-1:0
        I = I*z.gen_image + coeff(a.data, i)
      end
      return I::elem_type(C)
    end

    preimage = function(a::Generic.Res)
      R = z.codomain
      parent(a) == R || throw("mixed rings in preimage")
      g = gens(domain)
      im_gen = map(x -> z(x), g) # apply x -> z(x) to the generatotrs ## possibly should be cached and stored
          ## need to write the elements in a matrix, solve the eqn for a
      Mt = zero_matrix(base_ring(base_ring(R)), degree(R.modulus), degree(R.modulus))

      for i=1:degree(R.modulus)
        elem_to_mat_row!(Mt, i, im_gen[i])
      end
      b = zero_matrix(base_ring(base_ring(R)), 1, degree(R.modulus))
      elem_to_mat_row!(b, 1, a)
      s = solve_rational(Mt', b') # why, oh why is solve operating on columns????
      # This is the worst
      if isa(s, Tuple) ## again, why, oh why is solve doing things differently
                   ## over rings than fields?
        s = s[1] * inv(s[2]) # all rings here (type) are actually fields (math)
      end

      s = sum([s[i,1]*g[i] for i=1:length(g)])

      @assert parent(s) == domain
      return s::elem_type(domain)
    end

    z.header = MapHeader(domain, codomain, image, preimage)
    return z
  end
end

function ResidueRingPolyMap(domain::D, codomain::C, i::Generic.Res{T}) where {D, C, T}
  return ResidueRingPolyMap{D, C, T}(domain, codomain, i)
end

function ResidueRingPolyMap(domain::D, codomain::C, i::Generic.Res{T}, coeff_map::Map) where {D, C, T}
  return ResidueRingPolyMap{D, C, T}(domain, codomain, i, coeff_map)
end

mutable struct CoerceMap{D, C} <: Map{D, C, HeckeMap, CoerceMap}
  header::MapHeader{D, C}

  function CoerceMap{D, C}(domain::D, codomain::C) where {D, C}
    z = new{D, C}()
    z.header = MapHeader(domain, codomain)
    return z
  end

  function CoerceMap(domain::Generic.ResRing{S}, codomain::Generic.ResRing{T}) where {S, T <: PolyElem}
    z = new{Generic.ResRing{S}, Generic.ResRing{T}}()

    image = function(a::Generic.Res)
      return codomain(a)::elem_type(codomain)
    end

    preimage = function(a::Generic.Res)
      while parent(a) != domain
        degree(a.data)>0 && throw("Element not in subfield")
        a = coeff(a.data, 0)
      end
      return a::elem_type(domain)
    end

    z.header = MapHeader(domain, codomain, image, preimage)
    return z
  end
end

function CoerceMap(domain::Generic.ResRing{fmpz}, codomain::FqNmodFiniteField)
  z = CoerceMap{Generic.ResRing{fmpz}, FqNmodFiniteField}()

  image = function(a::Generic.Res{fmpz})
      parent(a) != domain && error("Element not in domain")
      return codomain(fmpz(a))
  end

  preimage = function(a::fq_nmod)
    parent(a) != codomain && error("Element not in codomain")
    a.length > 1 && throw("Element not in image")
    return domain(coeff(a, 0))::Generic.Res{fmpz}
  end

  z.header = MapHeader(domain, codomain, image, preimage)
  return z
end

function CoerceMap(domain::FqNmodFiniteField, codomain::Generic.ResRing{fq_nmod_poly})
  z = CoerceMap{FqNmodFiniteField, Generic.ResRing{fq_nmod_poly}}()

  image = function(a::fq_nmod)
    parent(a) != domain && error("Element not in domain")
    return codomain(a)::Generic.Res{fq_nmod_poly}
  end

  preimage = function(a::Generic.Res{fq_nmod_poly})
    degree(a.data) > 0 && throw("Element not in subfield")
    return domain(coeff(a.data, 0))::fq_nmod
  end

  z.header = MapHeader(domain, codomain, image, preimage)
  return z
end

function CoerceMap(domain::D, codomain::C) where {D, C}
  return CoerceMap{D, C}(domain, codomain)
end


text:
using Statistics
using LinearAlgebra
using Random

using Distributions

include("etkf.jl")
import .ETKF

include("ens_forecast.jl")
import .ens_forecast

include("models.jl")
import .Models

include("integrators.jl")
import .Integrators

Random.seed!(1)

models = [Models.lorenz63_err.func, Models.lorenz63_err3.func]
model_true = Models.lorenz63_true.func
n_models = length(models)
D = 3
obs_ops = [I(D), I(D)]
H = I(D)
ens_sizes = [40, 40]
model_sizes = [D, D]
integrator = Integrators.rk4
x0 = rand(D)
t0 = 0.0
Δt = 0.05
outfreq = 1
window = 1
transient = 2000
x0 = integrator(models[1], x0, t0, transient*outfreq*Δt, Δt, inplace=false)
R = Symmetric(diagm(0=>(0.1*std(x0, dims=1)[:]).^2))
x0 = x0[end, :]
#ensembles = [ens_forecast.init_ens(model=models[model], integrator=integrator,
#                                   x0=x0, t0=t0, outfreq=outfreq, Δt=Δt,
#                                   ens_size=ens_sizes[model]) for model=1:n_models]
#x0 = ensembles[1][:, end]
n_cycles = 1000
ρ = 1e-3

#model_errs = [ens_forecast.model_err(model_true=model_true, model_err=models[model],
#                                     integrator=integrator, x0=x0, t0=t0,
#                                     outfreq=outfreq, Δt=Δt, window=window,
#                                     n_samples=100)[1] for model=1:n_models]

infos = Vector(undef, n_models)
for model=1:length(models)
    model_errs = [0.1*diagm(0=>ones(D))]#Vector{Matrix{Float64}}(undef, 1)]
    biases = [zeros(D)]

    ensembles = [x0 .+ rand(MvNormal(R), ens_sizes[model])]

    info, _, _ = ens_forecast.mmda(x0=x0, ensembles=ensembles,
                         models=[models[model]],
                         model_true=model_true,
                         obs_ops=[obs_ops[model]], H=H,
                         model_errs=model_errs, biases=biases, integrator=integrator,
                         ens_sizes=[ens_sizes[model]], Δt=Δt, window=window,
                         n_cycles=n_cycles, outfreq=outfreq,
                         model_sizes=model_sizes, R=R, ρ=ρ)
    infos[model] = info
end

#model_errs = [0.1*diagm(0=>ones(D)), 0.1*diagm(0=>ones(D))]
# #biases = [nothing, nothing]
biases = [zeros(3), zeros(3)]

# #incs1 = hcat(info1.increments...)
# #incs2 = hcat(info2.increments...)
# #model_errs = [cov(incs1'), cov(incs2')]
# #biases = [mean(incs1, dims=2)[:], mean(incs2, dims=2)[:]]
# #biases = [zeros(3), zeros(3)]

# #ensembles = [x0 .+ rand(MvNormal(R), ens_sizes[model]) for model=1:n_models]
ensembles = [x0 .+ rand(MvNormal(R), ens_sizes[model]) for model=1:n_models]

model_errs = [mean(infos[1].Q_hist[500:end]), mean(infos[2].Q_hist[500:end])]

info_mm, _, _ = ens_forecast.mmda(x0=x0, ensembles=ensembles, models=models,
                         model_true=model_true, obs_ops=obs_ops, H=H,
                         model_errs=model_errs, biases=biases, integrator=integrator,
                         ens_sizes=ens_sizes, Δt=Δt, window=window,
                         n_cycles=n_cycles, outfreq=outfreq,
                         model_sizes=model_sizes, R=R, ρ=ρ, fixed=false)

max_stars_repo_name: eviatarbach/mmda

max_issues_repo_name: eviatarbach/mmda

max_forks_repo_name: eviatarbach/mmda


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct GenericDER <: DynamicInjection
        name::String
        Qref_Flag::Int
        PQ_Flag::Int
        Gen_Flag::Int
        PerOp_Flag::Int
        Recon_Flag::Int
        Trv::Float64
        VV_pnts::NamedTuple{(:V1, :V2, :V3, :V4), Tuple{Float64, Float64, Float64, Float64}}
        Q_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Tp::Float64
        e_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Kpq::Float64
        Kiq::Float64
        Iqr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        I_max::Float64
        Tg::Float64
        kWh_Cap::Float64
        SOC_ini::Float64
        SOC_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Trf::Float64
        fdbd_pnts::NamedTuple{(:fdbd1, :fdbd2), Tuple{Float64, Float64}}
        D_dn::Float64
        D_up::Float64
        fe_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Kpp::Float64
        Kip::Float64
        P_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        dP_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        T_pord::Float64
        rrpwr::Float64
        VRT_pnts::NamedTuple{(:vrt1, :vrt2, :vrt3, :vrt4, :vrt5), Tuple{Float64, Float64, Float64, Float64, Float64}}
        TVRT_pnts::NamedTuple{(:tvrt1, :tvrt2, :tvrt3), Tuple{Float64, Float64, Float64}}
        tV_delay::Float64
        VES_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        FRT_pnts::NamedTuple{(:frt1, :frt2, :frt3, :frt4), Tuple{Float64, Float64, Float64, Float64}}
        TFRT_pnts::NamedTuple{(:tfrt1, :tfrt2), Tuple{Float64, Float64}}
        tF_delay::Float64
        FES_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Pfa_ref::Float64
        Q_ref::Float64
        P_ref::Float64
        base_power::Float64
        states::Vector{Symbol}
        n_states::Int
        ext::Dict{String, Any}
        internal::InfrastructureSystemsInternal
    end

Parameters of a Generic Distributed Energy Resource Model. Based on https://scholarspace.manoa.hawaii.edu/bitstream/10125/70994/0304.pdf

# Arguments
- `name::String`
- `Qref_Flag::Int`: Reactive Power Control Mode. 1 VoltVar Control, 2 Constant Q Control, 3 Constant PF Control, validation range: `(1, 3)`
- `PQ_Flag::Int`: Active and reactive power priority mode. 0 for Q priority, 1 for P priority, validation range: `(0, 1)`
- `Gen_Flag::Int`: Define generator or storage system. 0 unit is a storage device, 1 unit is a generator, validation range: `(0, 1)`
- `PerOp_Flag::Int`: Defines operation of permisible region in VRT characteristic. 0 for cease, 1 for continuous operation, validation range: `(0, 1)`
- `Recon_Flag::Int`: Defines if DER can reconnect after voltage ride-through disconnection, validation range: `(0, 1)`
- `Trv::Float64`: Voltage measurement transducer's time constant, in s, validation range: `(0, nothing)`
- `VV_pnts::NamedTuple{(:V1, :V2, :V3, :V4), Tuple{Float64, Float64, Float64, Float64}}`: Y-axis Volt-var curve points (V1,V2,V3,V4)
- `Q_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Reactive power limits in pu (Q_min, Q_max)
- `Tp::Float64`: Power measurement transducer's time constant, in s, validation range: `(0, nothing)`
- `e_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Error limit in PI controller for q control (e_min, e_max)
- `Kpq::Float64`: PI controller proportional gain for q control, validation range: `(0, nothing)`
- `Kiq::Float64`: PI controller integral gain for q control, validation range: `(0, nothing)`
- `Iqr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Limit on rate of change for reactive current (pu/s) (Iqr_min, Iqr_max)
- `I_max::Float64`: Max. inverter's current, validation range: `(0, nothing)`
- `Tg::Float64`: Current control's time constant, in s, validation range: `(0, nothing)`
- `kWh_Cap::Float64`: BESS capacity in kWh, validation range: `(0, nothing)`
- `SOC_ini::Float64`: Initial state of charge (SOC) in pu, validation range: `(0, 1)`
- `SOC_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Battery's SOC limits (SOC_min, SOC_max)
- `Trf::Float64`: Time constant to estimate system frequency, in s, validation range: `(0, nothing)`
- `fdbd_pnts::NamedTuple{(:fdbd1, :fdbd2), Tuple{Float64, Float64}}`: Frequency error dead band thresholds `(fdbd1, fdbd2)`
- `D_dn::Float64`: reciprocal of droop for over-frequency conditions, in pu, validation range: `(0, nothing)`
- `D_up::Float64`: reciprocal of droop for under-frequency conditions, in pu, validation range: `(0, nothing)`
- `fe_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Frequency error limits in pu (fe_min, fe_max)
- `Kpp::Float64`: PI controller proportional gain for p control, validation range: `(0, nothing)`
- `Kip::Float64`: PI controller integral gain for p control, validation range: `(0, nothing)`
- `P_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Active power limits in pu (P_min, P_max)
- `dP_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Ramp rate limits for active power in pu/s (dP_min, dP_max)
- `T_pord::Float64`: Power filter time constant in s, validation range: `(0, nothing)`
- `rrpwr::Float64`: Ramp rate for real power increase following a fault, in pu/s, validation range: `(0, nothing)`
- `VRT_pnts::NamedTuple{(:vrt1, :vrt2, :vrt3, :vrt4, :vrt5), Tuple{Float64, Float64, Float64, Float64, Float64}}`: Voltage ride through v points (vrt1,vrt2,vrt3,vrt4,vrt5)
- `TVRT_pnts::NamedTuple{(:tvrt1, :tvrt2, :tvrt3), Tuple{Float64, Float64, Float64}}`: Voltage ride through time points (tvrt1,tvrt2,tvrt3)
- `tV_delay::Float64`: Time delay for reconnection after voltage ride-through disconnection, validation range: `(0, nothing)`
- `VES_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Min and max voltage for entering service (VES_min,VES_max)
- `FRT_pnts::NamedTuple{(:frt1, :frt2, :frt3, :frt4), Tuple{Float64, Float64, Float64, Float64}}`: Frequency ride through v points (frt1,frt2,frt3,frt4)
- `TFRT_pnts::NamedTuple{(:tfrt1, :tfrt2), Tuple{Float64, Float64}}`: Frequency ride through time points (tfrt1,tfrt2)
- `tF_delay::Float64`: Time delay for reconnection after frequency ride-through disconnection, validation range: `(0, nothing)`
- `FES_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Min and max frequency for entering service (FES_min,FES_max)
- `Pfa_ref::Float64`: Reference power factor, validation range: `(0, nothing)`
- `Q_ref::Float64`: Reference reactive power, in pu, validation range: `(0, nothing)`
- `P_ref::Float64`: Reference active power, in pu, validation range: `(0, nothing)`
- `base_power::Float64`: Base power
- `states::Vector{Symbol}`: The states of GenericDER depend on the Flags
- `n_states::Int`: The states of GenericDER depend on the Flags
- `ext::Dict{String, Any}`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct GenericDER <: DynamicInjection
    name::String
    "Reactive Power Control Mode. 1 VoltVar Control, 2 Constant Q Control, 3 Constant PF Control"
    Qref_Flag::Int
    "Active and reactive power priority mode. 0 for Q priority, 1 for P priority"
    PQ_Flag::Int
    "Define generator or storage system. 0 unit is a storage device, 1 unit is a generator"
    Gen_Flag::Int
    "Defines operation of permisible region in VRT characteristic. 0 for cease, 1 for continuous operation"
    PerOp_Flag::Int
    "Defines if DER can reconnect after voltage ride-through disconnection"
    Recon_Flag::Int
    "Voltage measurement transducer's time constant, in s"
    Trv::Float64
    "Y-axis Volt-var curve points (V1,V2,V3,V4)"
    VV_pnts::NamedTuple{(:V1, :V2, :V3, :V4), Tuple{Float64, Float64, Float64, Float64}}
    "Reactive power limits in pu (Q_min, Q_max)"
    Q_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Power measurement transducer's time constant, in s"
    Tp::Float64
    "Error limit in PI controller for q control (e_min, e_max)"
    e_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "PI controller proportional gain for q control"
    Kpq::Float64
    "PI controller integral gain for q control"
    Kiq::Float64
    "Limit on rate of change for reactive current (pu/s) (Iqr_min, Iqr_max)"
    Iqr_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Max. inverter's current"
    I_max::Float64
    "Current control's time constant, in s"
    Tg::Float64
    "BESS capacity in kWh"
    kWh_Cap::Float64
    "Initial state of charge (SOC) in pu"
    SOC_ini::Float64
    "Battery's SOC limits (SOC_min, SOC_max)"
    SOC_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Time constant to estimate system frequency, in s"
    Trf::Float64
    "Frequency error dead band thresholds `(fdbd1, fdbd2)`"
    fdbd_pnts::NamedTuple{(:fdbd1, :fdbd2), Tuple{Float64, Float64}}
    "reciprocal of droop for over-frequency conditions, in pu"
    D_dn::Float64
    "reciprocal of droop for under-frequency conditions, in pu"
    D_up::Float64
    "Frequency error limits in pu (fe_min, fe_max)"
    fe_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "PI controller proportional gain for p control"
    Kpp::Float64
    "PI controller integral gain for p control"
    Kip::Float64
    "Active power limits in pu (P_min, P_max)"
    P_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Ramp rate limits for active power in pu/s (dP_min, dP_max)"
    dP_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Power filter time constant in s"
    T_pord::Float64
    "Ramp rate for real power increase following a fault, in pu/s"
    rrpwr::Float64
    "Voltage ride through v points (vrt1,vrt2,vrt3,vrt4,vrt5)"
    VRT_pnts::NamedTuple{(:vrt1, :vrt2, :vrt3, :vrt4, :vrt5), Tuple{Float64, Float64, Float64, Float64, Float64}}
    "Voltage ride through time points (tvrt1,tvrt2,tvrt3)"
    TVRT_pnts::NamedTuple{(:tvrt1, :tvrt2, :tvrt3), Tuple{Float64, Float64, Float64}}
    "Time delay for reconnection after voltage ride-through disconnection"
    tV_delay::Float64
    "Min and max voltage for entering service (VES_min,VES_max)"
    VES_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Frequency ride through v points (frt1,frt2,frt3,frt4)"
    FRT_pnts::NamedTuple{(:frt1, :frt2, :frt3, :frt4), Tuple{Float64, Float64, Float64, Float64}}
    "Frequency ride through time points (tfrt1,tfrt2)"
    TFRT_pnts::NamedTuple{(:tfrt1, :tfrt2), Tuple{Float64, Float64}}
    "Time delay for reconnection after frequency ride-through disconnection"
    tF_delay::Float64
    "Min and max frequency for entering service (FES_min,FES_max)"
    FES_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Reference power factor"
    Pfa_ref::Float64
    "Reference reactive power, in pu"
    Q_ref::Float64
    "Reference active power, in pu"
    P_ref::Float64
    "Base power"
    base_power::Float64
    "The states of GenericDER depend on the Flags"
    states::Vector{Symbol}
    "The states of GenericDER depend on the Flags"
    n_states::Int
    ext::Dict{String, Any}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function GenericDER(name, Qref_Flag, PQ_Flag, Gen_Flag, PerOp_Flag, Recon_Flag, Trv, VV_pnts, Q_lim, Tp, e_lim, Kpq, Kiq, Iqr_lim, I_max, Tg, kWh_Cap, SOC_ini, SOC_lim, Trf, fdbd_pnts, D_dn, D_up, fe_lim, Kpp, Kip, P_lim, dP_lim, T_pord, rrpwr, VRT_pnts, TVRT_pnts, tV_delay, VES_lim, FRT_pnts, TFRT_pnts, tF_delay, FES_lim, Pfa_ref=0.0, Q_ref=0.0, P_ref=1.0, base_power=100.0, ext=Dict{String, Any}(), )
    GenericDER(name, Qref_Flag, PQ_Flag, Gen_Flag, PerOp_Flag, Recon_Flag, Trv, VV_pnts, Q_lim, Tp, e_lim, Kpq, Kiq, Iqr_lim, I_max, Tg, kWh_Cap, SOC_ini, SOC_lim, Trf, fdbd_pnts, D_dn, D_up, fe_lim, Kpp, Kip, P_lim, dP_lim, T_pord, rrpwr, VRT_pnts, TVRT_pnts, tV_delay, VES_lim, FRT_pnts, TFRT_pnts, tF_delay, FES_lim, Pfa_ref, Q_ref, P_ref, base_power, ext, PowerSystems.get_GenericDER_states(Qref_Flag)[1], PowerSystems.get_GenericDER_states(Qref_Flag)[2], InfrastructureSystemsInternal(), )
end

function GenericDER(; name, Qref_Flag, PQ_Flag, Gen_Flag, PerOp_Flag, Recon_Flag, Trv, VV_pnts, Q_lim, Tp, e_lim, Kpq, Kiq, Iqr_lim, I_max, Tg, kWh_Cap, SOC_ini, SOC_lim, Trf, fdbd_pnts, D_dn, D_up, fe_lim, Kpp, Kip, P_lim, dP_lim, T_pord, rrpwr, VRT_pnts, TVRT_pnts, tV_delay, VES_lim, FRT_pnts, TFRT_pnts, tF_delay, FES_lim, Pfa_ref=0.0, Q_ref=0.0, P_ref=1.0, base_power=100.0, states=PowerSystems.get_GenericDER_states(Qref_Flag)[1], n_states=PowerSystems.get_GenericDER_states(Qref_Flag)[2], ext=Dict{String, Any}(), internal=InfrastructureSystemsInternal(), )
    GenericDER(name, Qref_Flag, PQ_Flag, Gen_Flag, PerOp_Flag, Recon_Flag, Trv, VV_pnts, Q_lim, Tp, e_lim, Kpq, Kiq, Iqr_lim, I_max, Tg, kWh_Cap, SOC_ini, SOC_lim, Trf, fdbd_pnts, D_dn, D_up, fe_lim, Kpp, Kip, P_lim, dP_lim, T_pord, rrpwr, VRT_pnts, TVRT_pnts, tV_delay, VES_lim, FRT_pnts, TFRT_pnts, tF_delay, FES_lim, Pfa_ref, Q_ref, P_ref, base_power, states, n_states, ext, internal, )
end

# Constructor for demo purposes; non-functional.
function GenericDER(::Nothing)
    GenericDER(;
        name="init",
        Qref_Flag=1,
        PQ_Flag=0,
        Gen_Flag=0,
        PerOp_Flag=0,
        Recon_Flag=0,
        Trv=0,
        VV_pnts=(V1=0.0, V2=0.0, V3=0.0, V4=0.0),
        Q_lim=(min=0.0, max=0.0),
        Tp=0,
        e_lim=(min=0.0, max=0.0),
        Kpq=0,
        Kiq=0,
        Iqr_lim=(min=0.0, max=0.0),
        I_max=0,
        Tg=0,
        kWh_Cap=0,
        SOC_ini=0,
        SOC_lim=(min=0.0, max=0.0),
        Trf=0,
        fdbd_pnts=(fdbd1=0.0, fdbd2=0.0),
        D_dn=0,
        D_up=0,
        fe_lim=(min=0.0, max=0.0),
        Kpp=0,
        Kip=0,
        P_lim=(min=0.0, max=0.0),
        dP_lim=(min=0.0, max=0.0),
        T_pord=0,
        rrpwr=0,
        VRT_pnts=(vrt1=0.0, vrt2=0.0, vrt3=0.0, vrt4=0.0, vrt5=0.0),
        TVRT_pnts=(tvrt1=0.0, tvrt2=0.0, tvrt3=0.0),
        tV_delay=0,
        VES_lim=(min=0.0, max=0.0),
        FRT_pnts=(frt1=0.0, frt2=0.0, frt3=0.0, frt4=0.0),
        TFRT_pnts=(tfrt1=0.0, tfrt2=0.0),
        tF_delay=0,
        FES_lim=(min=0.0, max=0.0),
        Pfa_ref=0,
        Q_ref=0,
        P_ref=0,
        base_power=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`GenericDER`](@ref) `name`."""
get_name(value::GenericDER) = value.name
"""Get [`GenericDER`](@ref) `Qref_Flag`."""
get_Qref_Flag(value::GenericDER) = value.Qref_Flag
"""Get [`GenericDER`](@ref) `PQ_Flag`."""
get_PQ_Flag(value::GenericDER) = value.PQ_Flag
"""Get [`GenericDER`](@ref) `Gen_Flag`."""
get_Gen_Flag(value::GenericDER) = value.Gen_Flag
"""Get [`GenericDER`](@ref) `PerOp_Flag`."""
get_PerOp_Flag(value::GenericDER) = value.PerOp_Flag
"""Get [`GenericDER`](@ref) `Recon_Flag`."""
get_Recon_Flag(value::GenericDER) = value.Recon_Flag
"""Get [`GenericDER`](@ref) `Trv`."""
get_Trv(value::GenericDER) = value.Trv
"""Get [`GenericDER`](@ref) `VV_pnts`."""
get_VV_pnts(value::GenericDER) = value.VV_pnts
"""Get [`GenericDER`](@ref) `Q_lim`."""
get_Q_lim(value::GenericDER) = value.Q_lim
"""Get [`GenericDER`](@ref) `Tp`."""
get_Tp(value::GenericDER) = value.Tp
"""Get [`GenericDER`](@ref) `e_lim`."""
get_e_lim(value::GenericDER) = value.e_lim
"""Get [`GenericDER`](@ref) `Kpq`."""
get_Kpq(value::GenericDER) = value.Kpq
"""Get [`GenericDER`](@ref) `Kiq`."""
get_Kiq(value::GenericDER) = value.Kiq
"""Get [`GenericDER`](@ref) `Iqr_lim`."""
get_Iqr_lim(value::GenericDER) = value.Iqr_lim
"""Get [`GenericDER`](@ref) `I_max`."""
get_I_max(value::GenericDER) = value.I_max
"""Get [`GenericDER`](@ref) `Tg`."""
get_Tg(value::GenericDER) = value.Tg
"""Get [`GenericDER`](@ref) `kWh_Cap`."""
get_kWh_Cap(value::GenericDER) = value.kWh_Cap
"""Get [`GenericDER`](@ref) `SOC_ini`."""
get_SOC_ini(value::GenericDER) = value.SOC_ini
"""Get [`GenericDER`](@ref) `SOC_lim`."""
get_SOC_lim(value::GenericDER) = value.SOC_lim
"""Get [`GenericDER`](@ref) `Trf`."""
get_Trf(value::GenericDER) = value.Trf
"""Get [`GenericDER`](@ref) `fdbd_pnts`."""
get_fdbd_pnts(value::GenericDER) = value.fdbd_pnts
"""Get [`GenericDER`](@ref) `D_dn`."""
get_D_dn(value::GenericDER) = value.D_dn
"""Get [`GenericDER`](@ref) `D_up`."""
get_D_up(value::GenericDER) = value.D_up
"""Get [`GenericDER`](@ref) `fe_lim`."""
get_fe_lim(value::GenericDER) = value.fe_lim
"""Get [`GenericDER`](@ref) `Kpp`."""
get_Kpp(value::GenericDER) = value.Kpp
"""Get [`GenericDER`](@ref) `Kip`."""
get_Kip(value::GenericDER) = value.Kip
"""Get [`GenericDER`](@ref) `P_lim`."""
get_P_lim(value::GenericDER) = value.P_lim
"""Get [`GenericDER`](@ref) `dP_lim`."""
get_dP_lim(value::GenericDER) = value.dP_lim
"""Get [`GenericDER`](@ref) `T_pord`."""
get_T_pord(value::GenericDER) = value.T_pord
"""Get [`GenericDER`](@ref) `rrpwr`."""
get_rrpwr(value::GenericDER) = value.rrpwr
"""Get [`GenericDER`](@ref) `VRT_pnts`."""
get_VRT_pnts(value::GenericDER) = value.VRT_pnts
"""Get [`GenericDER`](@ref) `TVRT_pnts`."""
get_TVRT_pnts(value::GenericDER) = value.TVRT_pnts
"""Get [`GenericDER`](@ref) `tV_delay`."""
get_tV_delay(value::GenericDER) = value.tV_delay
"""Get [`GenericDER`](@ref) `VES_lim`."""
get_VES_lim(value::GenericDER) = value.VES_lim
"""Get [`GenericDER`](@ref) `FRT_pnts`."""
get_FRT_pnts(value::GenericDER) = value.FRT_pnts
"""Get [`GenericDER`](@ref) `TFRT_pnts`."""
get_TFRT_pnts(value::GenericDER) = value.TFRT_pnts
"""Get [`GenericDER`](@ref) `tF_delay`."""
get_tF_delay(value::GenericDER) = value.tF_delay
"""Get [`GenericDER`](@ref) `FES_lim`."""
get_FES_lim(value::GenericDER) = value.FES_lim
"""Get [`GenericDER`](@ref) `Pfa_ref`."""
get_Pfa_ref(value::GenericDER) = value.Pfa_ref
"""Get [`GenericDER`](@ref) `Q_ref`."""
get_Q_ref(value::GenericDER) = value.Q_ref
"""Get [`GenericDER`](@ref) `P_ref`."""
get_P_ref(value::GenericDER) = value.P_ref
"""Get [`GenericDER`](@ref) `base_power`."""
get_base_power(value::GenericDER) = value.base_power
"""Get [`GenericDER`](@ref) `states`."""
get_states(value::GenericDER) = value.states
"""Get [`GenericDER`](@ref) `n_states`."""
get_n_states(value::GenericDER) = value.n_states
"""Get [`GenericDER`](@ref) `ext`."""
get_ext(value::GenericDER) = value.ext
"""Get [`GenericDER`](@ref) `internal`."""
get_internal(value::GenericDER) = value.internal

"""Set [`GenericDER`](@ref) `Qref_Flag`."""
set_Qref_Flag!(value::GenericDER, val) = value.Qref_Flag = val
"""Set [`GenericDER`](@ref) `PQ_Flag`."""
set_PQ_Flag!(value::GenericDER, val) = value.PQ_Flag = val
"""Set [`GenericDER`](@ref) `Gen_Flag`."""
set_Gen_Flag!(value::GenericDER, val) = value.Gen_Flag = val
"""Set [`GenericDER`](@ref) `PerOp_Flag`."""
set_PerOp_Flag!(value::GenericDER, val) = value.PerOp_Flag = val
"""Set [`GenericDER`](@ref) `Recon_Flag`."""
set_Recon_Flag!(value::GenericDER, val) = value.Recon_Flag = val
"""Set [`GenericDER`](@ref) `Trv`."""
set_Trv!(value::GenericDER, val) = value.Trv = val
"""Set [`GenericDER`](@ref) `VV_pnts`."""
set_VV_pnts!(value::GenericDER, val) = value.VV_pnts = val
"""Set [`GenericDER`](@ref) `Q_lim`."""
set_Q_lim!(value::GenericDER, val) = value.Q_lim = val
"""Set [`GenericDER`](@ref) `Tp`."""
set_Tp!(value::GenericDER, val) = value.Tp = val
"""Set [`GenericDER`](@ref) `e_lim`."""
set_e_lim!(value::GenericDER, val) = value.e_lim = val
"""Set [`GenericDER`](@ref) `Kpq`."""
set_Kpq!(value::GenericDER, val) = value.Kpq = val
"""Set [`GenericDER`](@ref) `Kiq`."""
set_Kiq!(value::GenericDER, val) = value.Kiq = val
"""Set [`GenericDER`](@ref) `Iqr_lim`."""
set_Iqr_lim!(value::GenericDER, val) = value.Iqr_lim = val
"""Set [`GenericDER`](@ref) `I_max`."""
set_I_max!(value::GenericDER, val) = value.I_max = val
"""Set [`GenericDER`](@ref) `Tg`."""
set_Tg!(value::GenericDER, val) = value.Tg = val
"""Set [`GenericDER`](@ref) `kWh_Cap`."""
set_kWh_Cap!(value::GenericDER, val) = value.kWh_Cap = val
"""Set [`GenericDER`](@ref) `SOC_ini`."""
set_SOC_ini!(value::GenericDER, val) = value.SOC_ini = val
"""Set [`GenericDER`](@ref) `SOC_lim`."""
set_SOC_lim!(value::GenericDER, val) = value.SOC_lim = val
"""Set [`GenericDER`](@ref) `Trf`."""
set_Trf!(value::GenericDER, val) = value.Trf = val
"""Set [`GenericDER`](@ref) `fdbd_pnts`."""
set_fdbd_pnts!(value::GenericDER, val) = value.fdbd_pnts = val
"""Set [`GenericDER`](@ref) `D_dn`."""
set_D_dn!(value::GenericDER, val) = value.D_dn = val
"""Set [`GenericDER`](@ref) `D_up`."""
set_D_up!(value::GenericDER, val) = value.D_up = val
"""Set [`GenericDER`](@ref) `fe_lim`."""
set_fe_lim!(value::GenericDER, val) = value.fe_lim = val
"""Set [`GenericDER`](@ref) `Kpp`."""
set_Kpp!(value::GenericDER, val) = value.Kpp = val
"""Set [`GenericDER`](@ref) `Kip`."""
set_Kip!(value::GenericDER, val) = value.Kip = val
"""Set [`GenericDER`](@ref) `P_lim`."""
set_P_lim!(value::GenericDER, val) = value.P_lim = val
"""Set [`GenericDER`](@ref) `dP_lim`."""
set_dP_lim!(value::GenericDER, val) = value.dP_lim = val
"""Set [`GenericDER`](@ref) `T_pord`."""
set_T_pord!(value::GenericDER, val) = value.T_pord = val
"""Set [`GenericDER`](@ref) `rrpwr`."""
set_rrpwr!(value::GenericDER, val) = value.rrpwr = val
"""Set [`GenericDER`](@ref) `VRT_pnts`."""
set_VRT_pnts!(value::GenericDER, val) = value.VRT_pnts = val
"""Set [`GenericDER`](@ref) `TVRT_pnts`."""
set_TVRT_pnts!(value::GenericDER, val) = value.TVRT_pnts = val
"""Set [`GenericDER`](@ref) `tV_delay`."""
set_tV_delay!(value::GenericDER, val) = value.tV_delay = val
"""Set [`GenericDER`](@ref) `VES_lim`."""
set_VES_lim!(value::GenericDER, val) = value.VES_lim = val
"""Set [`GenericDER`](@ref) `FRT_pnts`."""
set_FRT_pnts!(value::GenericDER, val) = value.FRT_pnts = val
"""Set [`GenericDER`](@ref) `TFRT_pnts`."""
set_TFRT_pnts!(value::GenericDER, val) = value.TFRT_pnts = val
"""Set [`GenericDER`](@ref) `tF_delay`."""
set_tF_delay!(value::GenericDER, val) = value.tF_delay = val
"""Set [`GenericDER`](@ref) `FES_lim`."""
set_FES_lim!(value::GenericDER, val) = value.FES_lim = val
"""Set [`GenericDER`](@ref) `Pfa_ref`."""
set_Pfa_ref!(value::GenericDER, val) = value.Pfa_ref = val
"""Set [`GenericDER`](@ref) `Q_ref`."""
set_Q_ref!(value::GenericDER, val) = value.Q_ref = val
"""Set [`GenericDER`](@ref) `P_ref`."""
set_P_ref!(value::GenericDER, val) = value.P_ref = val
"""Set [`GenericDER`](@ref) `base_power`."""
set_base_power!(value::GenericDER, val) = value.base_power = val
"""Set [`GenericDER`](@ref) `ext`."""
set_ext!(value::GenericDER, val) = value.ext = val


text:
### A Pluto.jl notebook ###
# v0.17.2

using Markdown
using InteractiveUtils

# ╔═╡ d82325be-ee76-11eb-255e-21790b024173
# Activate the shared project environment
begin
    import Pkg
    Pkg.activate(Base.current_project())

	using IntrinsicDecoherence
end

# ╔═╡ d78603ff-5546-47e8-8f11-3f37bbd78e3b
# Calling plotting suite
begin
	using PyPlot
	using LaTeXStrings
	ion()
	rc("text", usetex = true)
end

# ╔═╡ e5378250-e42e-47a3-85bf-2c214abca6e9
md"## Intrinsic decoherence dynamics in the three-coupled harmonic oscillators interaction"

# ╔═╡ a6d194a6-5336-4495-9555-bd1b084b4203
# Set the constant parameters
begin
	const α = 4
	const ω = 4
	const λ = 0.5
end;

# ╔═╡ a5f8d238-8f9e-48a8-9e8c-292e0a884bd0
# For γ = 10
begin
	t_list = range(0, stop = 70, length = 1000)

	texv1_a = [exv_a(ω, λ, 0.1, 10, α, t_list[τ]) for τ in 1:length(t_list)]
	texv1_a1 = [exv_a₁(ω, λ, 0.1, 10, α, t_list[τ]) for τ in 1:length(t_list)]
	texv1_a2 = [exv_a₂(ω, λ, 0.1, 10, α, t_list[τ]) for τ in 1:length(t_list)]

	texv2_a = [exv_a(ω, λ, 0.5, 10, α, t_list[τ]) for τ in 1:length(t_list)]
	texv2_a1 = [exv_a₁(ω, λ, 0.5, 10, α, t_list[τ]) for τ in 1:length(t_list)]
	texv2_a2 = [exv_a₂(ω, λ, 0.5, 10, α, t_list[τ]) for τ in 1:length(t_list)]

	texv3_a = [exv_a(ω, λ, 1, 10, α, t_list[τ]) for τ in 1:length(t_list)]
	texv3_a1 = [exv_a₁(ω, λ, 1, 10, α, t_list[τ]) for τ in 1:length(t_list)]
	texv3_a2 = [exv_a₂(ω, λ, 1, 10, α, t_list[τ]) for τ in 1:length(t_list)]
end;

# ╔═╡ a31563ba-6c02-4092-bad0-765f820f7a38
begin
	figs, (axs1, axs2, axs3) = subplots(3,1, figsize = (6,10), sharex = true)

	axs1.plot(t_list, real(texv1_a1), "r", lw = 3, label = L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle")
	axs1.plot(t_list, real(texv1_a2), "b", lw = 3, label = L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle")
	axs1.plot(t_list, real(texv1_a), "g", lw = 3, label = L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle")

	axs2.plot(t_list, real(texv2_a1), "r", lw = 3, label = L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle")
	axs2.plot(t_list, real(texv2_a2), "b", ms = 4, lw = 3, label = L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle")
	axs2.plot(t_list, real(texv2_a), "g", linestyle = (0, (5, 5)), lw = 3, label = L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle")

	axs3.plot(t_list, real(texv3_a1), "r", lw = 3, label = L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle")
	axs3.plot(t_list, real(texv3_a2), "b", lw = 3, label = L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle")
	axs3.plot(t_list, real(texv3_a), "g", lw = 3, label = L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle")

	axs1.text(69, 0, L"(a)", fontsize=16)
	axs2.text(69, 0, L"(b)", fontsize=16)
	axs3.text(69, 0, L"(c)", fontsize=16)

	axs1.set_yticks([0, 4, 8, 12, 16])
	axs2.set_yticks([0, 4, 8, 12, 16])
	axs3.set_yticks([0, 4, 8, 12, 16])

	axins = axs3.inset_axes([0.4, 0.6, 0.25, 0.25])
	axins.plot(t_list, real(texv3_a1), "r")
	axins.plot(t_list, real(texv3_a2), "b")
	axins.plot(t_list, real(texv3_a), "g")
	axins.set_xlim([0,8])
	axins.set_title(L"Close view $0 \leq t \leq 8$", fontsize = 12)
	setp(axins.get_xticklabels(), fontsize = 12)
	setp(axins.get_yticklabels(), fontsize = 12)
	axins.set_yticks([0,8,16])
	axins.set_xticks([0,4,8])

	axs1.set_ylabel(L"\textrm{Number of photons}", fontsize = 20)
	setp(axs1.get_xticklabels(), fontsize = 18)
	setp(axs1.get_yticklabels(), fontsize = 18)
	axs1.legend(fontsize = 13)

	axs2.set_ylabel(L"\textrm{Number of photons}", fontsize = 20)
	setp(axs2.get_xticklabels(), fontsize = 18)
	setp(axs2.get_yticklabels(), fontsize = 18)
	axs2.legend(fontsize = 13)

	axs3.set_xlabel(L"t\; [s]", fontsize = 18)
	axs3.set_ylabel(L"\textrm{Number of photons}", fontsize = 20)
	setp(axs3.get_xticklabels(), fontsize = 18)
	setp(axs3.get_yticklabels(), fontsize = 18)
	axs3.legend(fontsize = 13)

	tight_layout()

	savefig("γ10_sp.pdf", dpi = 300)
	savefig("γ10_sp.png", dpi = 300)
end

# ╔═╡ e5d45e47-c94a-41f6-80fb-6874e3b84fbf
# For γ = 100
begin
	t_lisT = range(0, stop = 7, length = 500)

	texV1_a = [exv_a(ω, λ, 0.1, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]
	texV1_a1 = [exv_a₁(ω, λ, 0.1, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]
	texV1_a2 = [exv_a₂(ω, λ, 0.1, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]

	texV2_a = [exv_a(ω, λ, 0.5, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]
	texV2_a1 = [exv_a₁(ω, λ, 0.5, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]
	texV2_a2 = [exv_a₂(ω, λ, 0.5, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]

	texV3_a = [exv_a(ω, λ, 1, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]
	texV3_a1 = [exv_a₁(ω, λ, 1, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]
	texV3_a2 = [exv_a₂(ω, λ, 1, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]
end;

# ╔═╡ 07961c9f-34fd-401b-b7a9-f4ab25e6c94b
begin
	figS, (axS1, axS2, axS3) = subplots(3,1, figsize = (6,10), sharex = true)

	axS1.plot(t_lisT, real(texV1_a1), "r", lw = 3, label = L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle")
	axS1.plot(t_lisT, real(texV1_a2), "b", lw = 3, label = L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle")
	axS1.plot(t_lisT, real(texV1_a), "g", lw = 3, label = L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle")

	axS2.plot(t_lisT, real(texV2_a1), "r", lw = 3, label = L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle")
	axS2.plot(t_lisT, real(texV2_a2), "b", ms = 3, lw = 3, label = L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle")
	axS2.plot(t_lisT, real(texV2_a), "g", linestyle = (0, (5, 5)), lw = 3, label = L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle")

	axS3.plot(t_lisT, real(texV3_a1), "r", lw = 3, label = L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle")
	axS3.plot(t_lisT, real(texV3_a2), "b", lw = 3, label = L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle")
	axS3.plot(t_lisT, real(texV3_a), "g", lw = 3, label = L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle")

	axS1.text(6.9, 0, L"(a)", fontsize=16)
	axS2.text(6.9, 0, L"(b)", fontsize=16)
	axS3.text(6.9, 0, L"(c)", fontsize=16)

	axS1.set_yticks([0, 4, 8, 12, 16])
	axS2.set_yticks([0, 4, 8, 12, 16])
	axS3.set_yticks([0, 4, 8, 12, 16])

	axS1.set_ylabel(L"\textrm{Number of photons}", fontsize = 20)
	setp(axS1.get_xticklabels(), fontsize = 18)
	setp(axS1.get_yticklabels(), fontsize = 18)
	axS1.legend(fontsize = 13)

	axS2.set_ylabel(L"\textrm{Number of photons}", fontsize = 20)
	setp(axS2.get_xticklabels(), fontsize = 18)
	setp(axS2.get_yticklabels(), fontsize = 18)
	axS2.legend(fontsize = 13)

	axS3.set_xlabel(L"t\; [s]", fontsize = 18)
	axS3.set_ylabel(L"\textrm{Number of photons}", fontsize = 20)
	setp(axS3.get_xticklabels(), fontsize = 18)
	setp(axS3.get_yticklabels(), fontsize = 18)
	axS3.legend(fontsize = 13)

	tight_layout()

	savefig("γ100_sp.pdf", dpi = 300)
	savefig("γ100_sp.png", dpi = 300)
end

# ╔═╡ Cell order:
# ╟─e5378250-e42e-47a3-85bf-2c214abca6e9
# ╠═d82325be-ee76-11eb-255e-21790b024173
# ╠═d78603ff-5546-47e8-8f11-3f37bbd78e3b
# ╠═a6d194a6-5336-4495-9555-bd1b084b4203
# ╠═a5f8d238-8f9e-48a8-9e8c-292e0a884bd0
# ╠═a31563ba-6c02-4092-bad0-765f820f7a38
# ╠═e5d45e47-c94a-41f6-80fb-6874e3b84fbf
# ╠═07961c9f-34fd-401b-b7a9-f4ab25e6c94b

max_stars_repo_name: rurz/IntrinsicDecoherence

max_issues_repo_name: rurz/IntrinsicDecoherence

max_forks_repo_name: rurz/IntrinsicDecoherence


text:
using ControlSystems
using MAT

@testset "mapABCD" begin
    ntf = synthesizeNTF(5, 42, 1)
    forms = ["CRFB", "CRFF", "CRFBD", "CRFFD", "Stratos"]
    for i = 1:length(forms)
        form = forms[i]
        a1, g1, b1, c1 = realizeNTF(ntf, form)
        ABCD = stuffABCD(a1, g1, b1, c1, form)
        a2, g2, b2, c2 = mapABCD(ABCD, form)
        # test that mapABCD undoes stuffABCD, allow error in last digit
        @test isapprox(a1, a2, rtol=eps(Float64))
        @test isapprox(g1, g2, rtol=eps(Float64))
        @test isapprox(b1, b2, rtol=eps(Float64))
        @test isapprox(c1, c2, rtol=eps(Float64))
    end

    # "CIFB" and "CIFF" expect the real part of zeros to be 1
    ntf = zpk([1, 1+2im, 1-2im, 1+4im, 1-4im], [2, 2+1im, 2-1im, 2+3im, 2-3im], 1, 1)
    forms = ["CIFB", "CIFF"]
    for i = 1:length(forms)
        form = forms[i]
        a1, g1, b1, c1 = realizeNTF(ntf, form)
        ABCD = stuffABCD(a1, g1, b1, c1, form)
        a2, g2, b2, c2 = mapABCD(ABCD, form)
        # test that mapABCD undoes stuffABCD, allow error in last digit
        @test isapprox(a1, a2, rtol=eps(Float64))
        @test isapprox(g1, g2, rtol=eps(Float64))
        @test isapprox(b1, b2, rtol=eps(Float64))
        @test isapprox(c1, c2, rtol=eps(Float64))
    end

    # "PFF" and "DSFB" aren't supported by MATLAB toolbox yet
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma