text:
module TestOverlap
  using Test
  using EcologicalNetworks

  # Generate some data

  A = BipartiteNetwork([true false; true true; false true], [:a, :b, :c], [:d, :e])
  @test length(AJS(A)) == 2
  @test length(EAJS(A)) == 2
  @test AJS(A)[1].second == 0.5
  @test EAJS(A)[1].second == 0.5
  @test length(overlap(A)) == 2
  @test overlap(A)[1].second == 1

  B = BipartiteNetwork([true true true; true true false; true false false], [:a, :b, :c], [:d, :e, :f])
  @test AJS(B) == EAJS(B)
  @test length(AJS(B)) == 3
  @test length(overlap(B)) == 3
  @test filter(x -> x.first == Set((:a, :b)), overlap(B))[1].second == 2

end


text:
## This file is part of Dieter.jl : Model definition

# const hoursInYear = 8760     # Units: hours/yr
# const H2_energy_density = 39405.6 # Units : MWh / kg
# const infeas_cost = 1000     # Currency (cost)

"""
Build the JuMP model describing the optimization problem, specifying the `solver` to use.
The `Timestep` parameter should be 1 for half-hourly steps, 2 for hourly steps,
and 4 for 2-hourly steps. The data must match the time-steps of the `Timestep` parameter.
"""
function build_model!(dtr::DieterModel)

    dtr.model = Model(dtr.settings[:solver]; bridge_constraints = false)

    m = dtr.model

    ## Sync. the following with the file "load_names.jl"

    Timestep = dtr.settings[:timestep]

    periods = round(Int,hoursInYear*(2/Timestep))
    Hours = Base.OneTo(periods)
    Hours2 = Hours[2:end]
    # time_ratio relates generation levels in one time-step (e.g 1/2-hourly) to energy in MWh on an hourly basis
    # Used to compare energy to capacity in MW. A 1/2-hourly resolution means time_ratio = 1//2
    time_ratio = hoursInYear//length(Hours)
    dtr.settings[:time_ratio] = time_ratio
    # corr_factor = 1//time_ratio  # length(Hours)//hoursInYear

    dtr.sets[:Hours] = Hours

    cost_scaling = dtr.settings[:cost_scaling]
    loss_factor_tx = dtr.settings[:loss_factor_tx] # e.g 1.06 means a 6% increase in given Demand to account for transmission losses.

    # Set definitions

    # # Available sets, but not used directly in model:
        # Technologies = dtr.sets[:Technologies]  # Generation technologies
        # Storages = dtr.sets[:Storages]          # Storage technologies
        # Renewables = dtr.sets[:Renewables]           # Renewable generation technologies
        # Conventional = dtr.sets[:Conventional]       # Conventional generation technologies
        # Dispatchable = dtr.sets[:Dispatchable]       # Dispatchable generation technologies
        # NonDispatchable = dtr.sets[:NonDispatchable] # Non-dispatchable generation technologies

        # Nodes = dtr.sets[:Nodes]

    TxZones = dtr.sets[:TxZones]
    REZones = dtr.sets[:REZones]

    DemandZones = dtr.sets[:DemandZones]
    DemandRegions = dtr.sets[:DemandRegions]
    node2DemReg = dtr.parameters[:node_demreg_dict]

    # Mapping set definitions
    Arcs = dtr.sets[:Arcs]
    Arcs_From = dtr.sets[:Arcs_From]
    Arcs_REZones = dtr.sets[:Arcs_REZones]  # Connections between REZones and TxZones

    # ArcREZpairing = dtr.parameters[:ArcREZpairing]
    # Arcs_REZ_Pairs = intersect(Arcs,keys(ArcREZpairing))  # Filter pairings by allowed arcs

    Nodes_Techs = dtr.sets[:Nodes_Techs]
    Nodes_Storages = dtr.sets[:Nodes_Storages]

    # Nodes_Techs_Exi = dtr.sets[:Nodes_Techs_Exi]
    Nodes_Techs_New = dtr.sets[:Nodes_Techs_New]
    Nodes_Storages_New = dtr.sets[:Nodes_Storages_New]

    Nodes_Renew = dtr.sets[:Nodes_Renew]
    # Nodes_Conven = dtr.sets[:Nodes_Conven]
    Nodes_Dispatch = dtr.sets[:Nodes_Dispatch]
    Nodes_NonDispatch = dtr.sets[:Nodes_NonDispatch]

    # Subset of Nodes_Techs with availability traces
    Nodes_Avail_Techs = dtr.sets[:Nodes_Avail_Techs]

    # Test (TODO - relocate elsewhere):
    if !isempty(intersect(Nodes_Dispatch,Nodes_Avail_Techs))
        @warn "The sets Nodes_Dispatch and Nodes_Avail_Techs are not disjoint, and are required to be for correctness."
    end

    Nodes_Types = dtr.sets[:Nodes_Types]
    Nodes_Promotes = dtr.sets[:Nodes_Promotes]

    Nodes_Demand = dtr.sets[:Nodes_Demand]

    Nodes_RampingTechs = dtr.sets[:Nodes_RampingTechs]

    # Synchronous condensers
    Nodes_SynCons = dtr.sets[:Nodes_SynCons]

    ## Electric Vehicles
    EV = dtr.sets[:ElectricVehicles]

    ## Building heating
    BU = dtr.sets[:BuildingType]  # Building archtypes
    HP = dtr.sets[:HeatingType]   # Heating combination type

    ## Power-to-gas with Hydrogen
        # P2G = dtr.sets[:P2G]
        # G2P = dtr.sets[:G2P]
        # GasStorages = dtr.sets[:GasStorages]
    Nodes_P2G = dtr.sets[:Nodes_P2G]
    Nodes_G2P = dtr.sets[:Nodes_G2P]
    Nodes_GasStorages = dtr.sets[:Nodes_GasStorages]

    # Parameter definitions

    MarginalCost = dtr.parameters[:MarginalCost] # Units: currency/MWh; Marginal (variable) cost per unit of generated energy
    InvestmentCost = dtr.parameters[:InvestmentCost] # Units: currency/MW; Investment cost per unit of generation power capacity
    InvestmentCostPower = dtr.parameters[:InvestmentCostPower] # Units: currency/MW; Annualised investment cost per unit of storage power capacity
    InvestmentCostEnergy = dtr.parameters[:InvestmentCostEnergy] # Units: currency/MWh; Annualised investment cost per unit of storage energy capacity
    FixedCost = dtr.parameters[:FixedCost] # Units: currency/MW; Fixed cost per unit of power capacity
    Availability = dtr.parameters[:Availability] # Units: none; Available energy from renewables as fraction in [0,1] of installed capacity within a time-interval
    CapacityDerating = dtr.parameters[:CapacityDerating] # Units: [0,1]; Reduction in nominal capacity as fraction in [0,1] of installed capacity within a time-interval
    PeakContribution = dtr.parameters[:PeakContribution] # Units: [0,1]; Fraction of nominal capacity available from renewable tech. at peak demand.
    MaxCapacity = dtr.parameters[:MaxCapacity] # Units: MW; Maximum installable capacity - power
    MaxEnergy = dtr.parameters[:MaxEnergy] # Units: MWh; Maximum installable storage capacity energy
    MaxEtoP_ratio = dtr.parameters[:MaxEnergyToPowerRatio] # Units: hours; Maximum ratio of stored energy to power delivery
    Efficiency = dtr.parameters[:Efficiency] # Units: [0,1]; Combustion/Storage roundtrip efficiency
    StartLevel = dtr.parameters[:StartLevel] # Units: [0,1]; Initial storage level as fraction of storage energy installed

    CarbonContent = dtr.parameters[:CarbonContent] # Units: t-CO2/MWh-thermal; CO2 equivalent content per unit fuel used by tech.
    CarbonBudget = dtr.settings[:carbon_budget] # Units: t-CO2

    InvestmentCostSynCon = dtr.parameters[:InvestmentCostSynCon]

    OperatingReserve = dtr.parameters[:OperatingReserve] # Units: MW; Regional operating reserve requirement (e.g. within state)
    MinStableGen = dtr.parameters[:MinStableGen] # Units: [0,1]; Minimum stable operational level for generation as fraction of Capacity in MW.

    # WindLimit = dtr.parameters[:WindLimit]
    # SolarLimit = dtr.parameters[:SolarLimit]
    TotalBuildCap = dtr.parameters[:TotalBuildCap]

    ExpansionLimit_REZ = dtr.parameters[:ExpansionLimit]
    ExpansionLimit_Tx = dtr.parameters[:ExpansionLimit_Tx]
    TransferCapacity = dtr.parameters[:TransferCapacity] # Units: MW; Interconnector power transfer capability

    InvestmentCostTransExp = filter(x -> !ismissing(x[2]) && x[1] in Arcs, dtr.parameters[:InvestmentCostTransExp]) # Units:  currency/MW; ; REZ capacity expansion costs, amortised to year
    InvestmentCostREZ_Exp  = filter(x -> !ismissing(x[2]) && x[1] in Arcs_REZones, dtr.parameters[:InvestmentCostTransExp]) # Units:  currency/MW; ; Transmission interconnector capacity expansion costs, amortised to year

    ArcWeight = dtr.parameters[:ArcWeight]

    Load = dtr.parameters[:Load] # Units: MWh per time-interval; wholesale energy demand within a time-interval (e.g. hourly or 1/2-hourly)
    NegOpDemand = dtr.parameters[:NegOpDemand] # Units: MWh per time-interval; wholesale energy produced from behind-the-meter generation (e.g. PV or EV)

    CurtailmentCost = dtr.settings[:cu_cost] # Units: currency/MWh; Cost per unit of generated energy that is curtailed

    LoadIncreaseCost = dtr.parameters[:LoadIncreaseCost] # Units: $/MW; Load change costs for changing generation upward
    LoadDecreaseCost = dtr.parameters[:LoadDecreaseCost] # Units: $/MW; Load change costs for changing generation downward

    MaxRampUpPerHour = dtr.parameters[:MaxRampUpPerHour]      # Units MW/h; Maximum ramping up of technology
    MaxRampDownPerHour = dtr.parameters[:MaxRampDownPerHour]  # Units MW/h; Maximum ramping down of technology

    MinimumRenewShare = dtr.settings[:min_res_system]

    CapAdd = dtr.parameters[:CapAdd]

    ## End of sync. with "load_names.jl"

    # Transformations on the minimum renewable share settings:
    min_res_dict = Dict{String,Float64}()
    if isa(dtr.settings[:min_res],Number)
        for n in DemandRegions
            min_res_dict[n] = dtr.settings[:min_res]
        end
    elseif isa(dtr.settings[:min_res],Dict)
        min_res_dict = dtr.settings[:min_res]
        dtr.settings[:min_res] = minimum(collect(values(min_res_dict)))  # overwrite as a number, the smallest share.
        println("The `min_res` setting is now equal to $(dtr.settings[:min_res])")
    else
        error("The parameter `min_res_dict` is not of type `Number` or `Dict`.")
    end
    dtr.settings[:min_res_dict] = min_res_dict
    @info "The `min_res_dict` setting is currently equal to "
    @info min_res_dict
    for (k,v) in min_res_dict
        if !(0 <= v <= 100)
            error("The setting in `min_res_dict` of value $v for Demand Region $k is not in the interval [0,100].")
        elseif (0 < v <= 1)
            @warn "Note: The setting in `min_res_dict` of value $v for Demand Region $k is a percentage, not necessarily a fraction in the interval (0,1]."
        end
    end

    @info "Start of model building:"

    prog = Progress(4, dt=0.01, desc="Building Model...         \n", barlen=30)

# %% * ----------------------------------------------------------------------- *
#    ***** Variable definitions *****
#    * ----------------------------------------------------------------------- *

    shorter =  Dict("Total_cost_objective" => "Z",
                    "Generation_level" => "G",
                    "Generation_upshift" => "G_UP",
                    "Generation_downshift" => "G_DO",
                    "Generation_renewable" => "G_REZ",
                    "Generation_transmission_zones" => "G_TxZ",
                    "Storage_inflow" => "STO_IN",
                    "Storage_outflow" => "STO_OUT",
                    "Storage_level" => "STO_L",
                    "Technology_capacity" => "N_TECH",
                    "Renewable_capacity_expand" => "N_REZ_EXP",
                    "REZ_expand_by_transmission" => "N_REZ_EXP_TX",
                    "Storage_build_energy" => "N_STO_E",
                    "Storage_capacity" => "N_STO_P",
                    "SynCon_capacity" => "N_SYNC",
                    "Inertia_level" => "IN_TxZ",
                    "Internodal_flow" => "FLOW",
                    "Internodal_flow_expansion" => "N_IC_EXP",
                    "H2_power_to_gas" => "H2_P2G",
                    "H2_gas_to_power" => "H2_G2P",
                    "H2_storage_level" => "H2_GS_L",
                    "H2_storage_inflow" => "H2_GS_IN",
                    "H2_storage_outflow" => "H2_GS_OUT",
                    "H2_P2G_capacity" => "N_P2G",
                    "H2_G2P_capacity" => "N_G2P",
                    "H2_storage_capacity" => "N_GS",
                    "EV_charging" => "EV_CHARGE",
                    "EV_discharging" => "EV_DISCHARGE",
                    "EV_charge_level" => "EV_L",
                    "EV_PHEV_fuel_use" => "EV_PHEVFUEL",
                    "EV_infeasible" => "EV_INF",
                    "Heat_storage_level" => "HEAT_STO_L",
                    "Heat_heat_pump_" => "HEAT_HP_IN",
                    "Heat_infeasible" => "HEAT_INF")

    @info "Variable definitions."
    @variables(m, begin
            Z, (base_name=shorter["Total_cost_objective"], lower_bound=0)
            G[Nodes_Techs, Hours], (base_name=shorter["Generation_level"], lower_bound=0) # Units: MWh per time-interval; Generation level - all generation tech.
            G_UP[Nodes_Dispatch, Hours] , (base_name=shorter["Generation_upshift"], lower_bound=0)  # Units: MWh per time-interval; Generation level change up
            G_DO[Nodes_Dispatch, Hours], (base_name=shorter["Generation_downshift"], lower_bound=0) # Units: MWh per time-interval; Generation level change down
            # G_INF[Nodes, Hours], (base_name=shorter["Generation_infeasible"], lower_bound=0) # Units: MWh per time-interval; Infeasibility term for Energy Balance
            # CU[Nodes_NonDispatch, Hours], (base_name=shorter["Curtailment_renewables"], lower_bound=0) # Units: MWh per time-interval; Non-dispatchable curtailment
            G_REZ[REZones,Hours], (base_name=shorter["Generation_renewable"], lower_bound=0) # Units: MWh per time-interval; Generation level - renewable energy zone tech. & storage
            G_TxZ[TxZones,Hours], (base_name=shorter["Generation_transmission_zones"], lower_bound=0) # Units: MWh per time-interval; Generation level - transmission zone tech. & storage
            STO_IN[Nodes_Storages, Hours], (base_name=shorter["Storage_inflow"], lower_bound=0) # Units: MWh per time-interval; Storage energy inflow
            STO_OUT[Nodes_Storages, Hours], (base_name=shorter["Storage_outflow"], lower_bound=0) # Units: MWh per time-interval; Storage energy outflow
            STO_L[Nodes_Storages, Hours], (base_name=shorter["Storage_level"], lower_bound=0) # Units: MWh at a given time-interval; Storage energy level
            N_TECH[Nodes_Techs], (base_name=shorter["Technology_capacity"], lower_bound=0) # Units: MW; Technology capacity built
            N_REZ_EXP[REZones], (base_name=shorter["Renewable_capacity_expand"], lower_bound=0) # Units: MW; Renewable technology transmission capacity built
            N_REZ_EXP_TX[REZones], (base_name=shorter["REZ_expand_by_transmission"], lower_bound=0) # Units: MW; Free REZ transmission capacity built when tranmission network is expanded
            N_STO_E[Nodes_Storages], (base_name=shorter["Storage_build_energy"], lower_bound=0) # Units: MWh; Storage energy technology built
            N_STO_P[Nodes_Storages], (base_name=shorter["Storage_capacity"], lower_bound=0) # Units: MW; Storage loading and discharging power capacity built
            N_SYNC[Nodes_SynCons], (base_name=shorter["SynCon_capacity"], lower_bound=0) # Units: MW; synchronous condenser capacity
            FLOW[Arcs,Hours], (base_name=shorter["Internodal_flow"]) # Units: MWh; Power flow between nodes in topology
            N_IC_EXP[Arcs], (base_name=shorter["Internodal_flow_expansion"], lower_bound=0) # Units: MW; Power flow expansion betweeen TxZones, and REZones to TxZones
            #
            EV_CHARGE[EV, Hours], (base_name=shorter["EV_charging"], lower_bound=0) # Units: MWh per time-interval; Electric vehicle charge for vehicle profile in set EV
            EV_DISCHARGE[EV, Hours], (base_name=shorter["EV_discharging"], lower_bound=0) # Units: MWh per time-interval; Electric vehicle dischargw for vehicle profile in set EV
            EV_L[EV, Hours], (base_name=shorter["EV_charge_level"], lower_bound=0) # Units: MWh at a given time-interval; Electric vehicle charging level for vehicle profile in set EV
            EV_PHEVFUEL[EV, Hours], (base_name=shorter["EV_PHEV_fuel_use"], lower_bound=0) #  Plug-in hybrid electric vehicle conventional fuel use
            EV_INF[EV, Hours], (base_name=shorter["EV_infeasible"], lower_bound=0) # Units: MWh per time-interval; Infeasibility term for Electric vehicle energy balance
            #
            H2_P2G[Nodes_P2G, Hours], (base_name=shorter["H2_power_to_gas"], lower_bound=0) # Units: MWh per time-interval; Power-to-gas energy conversion
            H2_G2P[Nodes_G2P, Hours], (base_name=shorter["H2_gas_to_power"], lower_bound=0) # Units: MWh per time-interval; Gas-to-power energy conversion
            H2_GS_L[Nodes_GasStorages, Hours], (base_name=shorter["H2_storage_level"], lower_bound=0) # Units: tonne-H2 at a given time-interval; Current gas storage level
            H2_GS_IN[Nodes_GasStorages, Hours], (base_name=shorter["H2_storage_inflow"], lower_bound=0) # Units: tonne-H2 at a given time-interval; Current gas storage input
            H2_GS_OUT[Nodes_GasStorages, Hours], (base_name=shorter["H2_storage_outflow"], lower_bound=0) # Units: tonne-H2 at a given time-interval; Current gas storage output
            N_P2G[Nodes_P2G], (base_name=shorter["H2_P2G_capacity"], lower_bound=0) # Units: MW; Power-to-gas capacity
            N_G2P[Nodes_G2P], (base_name=shorter["H2_G2P_capacity"], lower_bound=0) # Units: MW; Gas-to-power capacity
            N_GS[Nodes_GasStorages], (base_name=shorter["H2_storage_capacity"], lower_bound=0) # Units; tonne-H2 ; Gas storage (energy) capacity
            #
            HEAT_STO_L[BU,HP,Hours], (base_name=shorter["Heat_storage_level"], lower_bound=0) # Units: MWh at a given time-interval; Heating: storage level
            HEAT_HP_IN[BU,HP, Hours], (base_name=shorter["Heat_heat_pump_"], lower_bound=0)   # Units: MWh per time-interval; Heating: electricity demand from heat pump
            HEAT_INF[BU,HP, Hours], (base_name=shorter["Heat_infeasible"], lower_bound=0)  # Units: MWh per time-interval; Heating: Infeasibility term for Electric vehicle energy balance
    end)

# %% * ----------------------------------------------------------------------- *
#    ***** Objective function *****
#    * ----------------------------------------------------------------------- *

    @objective(m, Min, Z)

    @info "Objective function."
    @constraint(m, ObjectiveFunction,
cost_scaling*(sum(MarginalCost[n,t] * G[(n,t),h] for (n,t) in Nodes_Dispatch, h in Hours)

            + LoadIncreaseCost[n,t] * sum(G_UP[(n,t),h] for (n,t) in Nodes_Dispatch, h in Hours2)
            + LoadDecreaseCost[n,t] * sum(G_DO[(n,t),h] for (n,t) in Nodes_Dispatch, h in Hours)

            # + sum(CurtailmentCost * CU[(n,t),h] for (n,t) in Nodes_NonDispatch, h in Hours)

            # + sum(infeas_cost * G_INF[n,h] for n in Nodes, h in Hours)

            + sum(MarginalCost[n,sto] * (STO_OUT[(n,sto),h] + STO_IN[(n,sto),h])
                for (n,sto) in Nodes_Storages, h in Hours)

            + sum(MarginalCost[ev] * EV_DISCHARGE[ev, h]
                  + EvFuel[ev] * EV_PHEVFUEL[ev, h]
                  + infeas_cost * EV_INF[ev, h] for ev in EV, h in Hours)

            + sum(MarginalCost[n,p2g] * H2_P2G[(n,p2g),h] for (n,p2g) in Nodes_P2G, h in Hours)
            + sum(MarginalCost[n,g2p] * H2_G2P[(n,g2p),h] for (n,g2p) in Nodes_G2P, h in Hours)

            + sum(infeas_cost * HEAT_INF[bu,hp,h] for bu in BU, hp in HP, h in Hours)
            )

        +
cost_scaling*(sum(InvestmentCost[n,t] * N_TECH[(n,t)] for (n,t) in Nodes_Techs_New)
            + sum(InvestmentCostPower[n,sto] * N_STO_P[(n,sto)] for (n,sto) in Nodes_Storages_New)
            + sum(InvestmentCostEnergy[n,sto] * N_STO_E[(n,sto)] for (n,sto) in Nodes_Storages_New)

            + sum(InvestmentCostSynCon[n,syn] * N_SYNC[(n,syn)] for (n,syn) in Nodes_SynCons)

            + sum(FixedCost[n,t] * N_TECH[(n,t)] for (n,t) in Nodes_Techs)
            + sum(FixedCost[n,sto] * N_STO_P[(n,sto)] for (n,sto) in Nodes_Storages)

            + sum(InvestmentCostREZ_Exp[rez,txz] * N_REZ_EXP[rez] for (rez,txz) in keys(InvestmentCostREZ_Exp) )
            + sum(ArcWeight[from,to]*InvestmentCostTransExp[from,to] * N_IC_EXP[(from,to)] for (from,to) in Arcs)

            + sum(InvestmentCost[n,p2g] * N_P2G[(n,p2g)] for (n,p2g) in Nodes_P2G)
            + sum(InvestmentCost[n,g2p] * N_G2P[(n,g2p)] for (n,g2p) in Nodes_G2P)
            + sum(InvestmentCost[n,gs] * N_GS[(n,gs)] for (n,gs) in Nodes_GasStorages)

            + sum(FixedCost[n,p2g] * N_P2G[(n,p2g)] for (n,p2g) in Nodes_P2G)
            + sum(FixedCost[n,g2p] * N_G2P[(n,g2p)] for (n,g2p) in Nodes_G2P)
            + sum(FixedCost[n,gs] * N_GS[(n,gs)] for (n,gs) in Nodes_GasStorages)
           )
        == Z
    );

    next!(prog)
    println("\n")

# %% * ----------------------------------------------------------------------- *
#    ***** Energy balance and load levels *****
#    * ----------------------------------------------------------------------- *

## Toy example to understand index sets made from tuples that are used to index constraints:
# D = ["a","b"]; T = [3,4]; NR = [("a",1),("a",2),("b",1)]; mod = JuMP.Model();
# @variable(mod, x[NR,T])  # Note NR is a set of tuples.
# This constraint form is correct for using the constraint index 'for n in D'
# and the first tuple entry is then also `n`:
# @constraint(mod, [n=D], sum(x[(i,j),t] for (i,j) in NR if i == n for t in T) <= 1)
# but NOT this one (a new scope is created inside the sum for `n`):
# @constraint(mod, Test[n=D], sum(x[(n,j)] for (n,j) in NR) <= 1)
# and not this one (do not use commas in comprehension statements, use `for` to concatenate):
# @constraint(mod, [n=D], sum(x[(i,j),t] for (i,j) in NR if i == n, t in T) <= 1)

    # The implicit design assumption of this constraint is that DemandRegions are directly related to the Dispatch region level.
    # We may overcome this with nested sums from lower levels.
    # In particular, G_REZ is an aggregating variable for generation from Renewables and used directly in other constraints
    # while G_TxZ aggregates Dispatchable generation from Transmission regions.

    @info "Definition of REZone generation book-keeping variables"
    @constraint(m, REZoneGen[rez=REZones,h=Hours],
        G_REZ[rez,h] ==
            sum(G[(z,t),h] for (z,t) in Nodes_Techs if z == rez)
         +  sum(STO_OUT[(q,sto),h] - STO_IN[(q,sto),h] for (q,sto) in Nodes_Storages if q == rez)
    );

    @info "Definition of TxZone generation book-keeping variables"
    @constraint(m, TxZoneGen[zone=TxZones,h=Hours],
        G_TxZ[zone,h] ==
            sum(G[(z,t),h] for (z,t) in Nodes_Techs if z == zone)
         +  sum(G_REZ[rez,h] for (rez, z) in Nodes_Promotes if z == zone)
         +  sum(STO_OUT[(z,sto),h] - STO_IN[(z,sto),h] for (z,sto) in Nodes_Storages if z == zone)
         -  sum(FLOW[(from,to),h] for (from,to) in Arcs if from == zone)
    );

# # TODO : An explicit assumption in EnergyBalance is that H2_G2P and H2_P2G are
# indexed over DemandZones; we should add additional checks to verify this holds.

    # Energy balance at each demand node:
    @info "Energy balance at each demand node: supply equals or exceeds demand."
    @constraint(m, EnergyBalance[n=DemandZones,h=Hours],
      sum(G_TxZ[zone,h] for (zone, d) in Nodes_Demand if d == n)
        # + sum(EV_DISCHARGE[ev,h] for ev in EV)
        + sum(H2_G2P[(zone,g2p),h] for (zone,g2p) in Nodes_G2P if zone == n)
        - sum(NegOpDemand[zone,hour] for (zone,hour) in keys(NegOpDemand) if (zone == n && hour == h))
        >=
      loss_factor_tx*(Load[n,h]
                    + sum(H2_P2G[(zone,p2g),h] for (zone,p2g) in Nodes_P2G if zone == n))
      );
      # + sum(HEAT_HP_IN[bu,hp,h] for bu in BU for hp in HP)
      # + sum(EV_CHARGE[ev,h] for ev in EV)

    # Energy flow reflexive constraint:
    @info "Energy flow reflexive constraint."
    @constraint(m, FlowEnergyReflex[(from,to)=Arcs,h=Hours; from in Arcs_From],
        FLOW[(from,to),h] + FLOW[(to,from),h] == 0
    );

    # Energy flow bounds:
    @info "Energy flow bounds."
    @constraint(m, FlowEnergyUpperBound[(from,to)=Arcs,h=Hours],
        FLOW[(from,to),h] <= time_ratio * ( TransferCapacity[(from,to)] + N_IC_EXP[(from,to)] + CapAdd[:N_IC_EXP][(from,to)])
    );

    # @info "Transmission expansion upper bounds."
    # @constraint(m, FlowExpandUpperBound[(from,to)=Arcs],
    #     N_IC_EXP[(from,to)] + trunc(CapAdd[:N_IC_EXP][(from,to)]) <= ExpansionLimit_Tx[(from,to)]
    # );

    @info "Energy flow expansion symmetry."
    @constraint(m, FlowEnergySymmetry[(from,to)=Arcs; from in Arcs_From],
        N_IC_EXP[(from,to)] == N_IC_EXP[(to,from)]
    );

    # Generation level start - con2b_loadlevelstart
    @info "Generation level start."
    @constraint(m, GenLevelStart[(n,t)=Nodes_Dispatch],
               G[(n,t),Hours[1]] == G_UP[(n,t),Hours[1]]
    );

    # Generation level - con2a_loadlevel
    @info "Generation level."
    @constraint(m, GenLevelUpdate[(n,t)=Nodes_Dispatch,h=Hours2],
        G[(n,t),h] == G[(n,t),h-1] + G_UP[(n,t),h] - G_DO[(n,t),h]
    );

    next!(prog)
    println("\n")

    # Variable upper bound on dispatchable generation by capacity
    @info "Variable upper bound on dispatchable generation by capacity."
    @constraint(m, MaxGenerationDisp[(n,t)=Nodes_Dispatch,h=Hours],
        G[(n,t),h] <= CapacityDerating[n,t,h] * time_ratio * (N_TECH[(n,t)] + CapAdd[:N_TECH][(n,t)])
    );

    @info "Variable upper bound on non-dispatchable generation by capacity."
    @constraint(m, MaxGenerationNonDisp[(n,t)=Nodes_Avail_Techs,h=Hours],
        # G[(n,t),h] + CU[(n,t),h] == Availability[n,t,h] * time_ratio * N_TECH[(n,t)]
        G[(n,t),h] <= Availability[n,t,h] * time_ratio * (N_TECH[(n,t)] + CapAdd[:N_TECH][(n,t)])
    );

    # Maximum capacity allowed
    if !(isDictAllMissing(MaxCapacity))
        @info "Maximum capacity allowed."
        @constraint(m, MaxCapacityBound[(n,t)=Nodes_Techs; !(MaxCapacity[n,t] |> ismissing)],
            N_TECH[(n,t)] + CapAdd[:N_TECH][(n,t)] <= MaxCapacity[n,t]
        );
    end

    # Maximum generated energy allowed.
    if !(isDictAllMissing(MaxEnergy))
        @info "Maximum generated energy allowed."
        @constraint(m, MaxEnergyGenerated[(n,t)=Nodes_Techs; !(MaxEnergy[n,t] |> ismissing)],
            sum(G[(n,t),h] for h in Hours) <= MaxEnergy[n,t]
        );
    end

    # Renewable energy zone build limits.
    @info "Renewable energy zone build limits."
    @constraint(m, REZBuildLimits[rez=REZones,h=Hours],
        # sum(G[(z,t),h] for (z,t) in Nodes_Avail_Techs if z == rez)
        sum(G[(z,t),h] for (z,t) in Nodes_Techs if (z == rez && !occursin(r"Hydro_",t))) ## TODO: remove this hard-coding!
            <= time_ratio*(TotalBuildCap[rez] +  N_REZ_EXP[rez] + CapAdd[:N_REZ_EXP][rez]) # + N_REZ_EXP_TX[rez]
    );

    @info "Renewable energy zone expansion limits."
    @constraint(m, REZExpansionBound[rez=REZones],
        N_REZ_EXP[rez] + trunc(CapAdd[:N_REZ_EXP][rez]) <= ExpansionLimit_REZ[rez]
    );

    @info "Renewable energy zone expansion link to transmission expansion."
    @constraint(m, REZpairExpansion[rez=REZones],
        N_REZ_EXP_TX[rez] == 0  # NOTE: This constraint is now inactive.
            # sum(N_IC_EXP[(from,to)]/0.75 for (from,to) in Arcs_REZ_Pairs if ArcREZpairing[(from,to)] == rez)
    );

    next!(prog)
    println("\n")

# %% * ----------------------------------------------------------------------- *
#    ***** Operational conditions and constraints *****
#    * ----------------------------------------------------------------------- *

    # Ramping Limits

    if !(isDictAllMissing(MaxRampUpPerHour))
        @info "Maximum ramp up rates"
        @constraint(m, RampingUpLimits[(n,t)=Nodes_RampingTechs, h=Hours2],
            G_UP[(n,t),h] <= time_ratio * MaxRampUpPerHour[t]
        );
    end

    if !(isDictAllMissing(MaxRampDownPerHour))
        @info "Maximum ramp down rates"
        @constraint(m, RampingDownLimits[(n,t)=Nodes_RampingTechs, h=Hours2],
             G_DO[(n,t),h] <= time_ratio * MaxRampDownPerHour[t]
        );
    end

    # Operating reserves
    @info "Operating reserve margin"
    @constraint(m, OperatingReserve[dr=DemandRegions, h=Hours],
        sum( CapacityDerating[n,t,h] * time_ratio * (N_TECH[(n,t)] + CapAdd[:N_TECH][(n,t)]) - G[(n,t),h]
            for (n,t) in Nodes_Dispatch if node2DemReg[n] == dr)
      + sum( PeakContribution[n,t] * time_ratio * (N_TECH[(n,t)] + CapAdd[:N_TECH][(n,t)])
            for (n,t) in Nodes_Avail_Techs if node2DemReg[n] == dr)  # or replace `Nodes_Avail_Techs` with `setdiff(Nodes_Avail_Techs,Nodes_Dispatch)`
      + sum( sqrt(Efficiency[n,sto])*STO_L[(n,sto), h]
            for (n,sto) in Nodes_Storages if node2DemReg[n] == dr)
      >= time_ratio * OperatingReserve[dr]
    );

    # Minimum stable generation levels (for inflexible technologies)

    if !(isempty(keys(MinStableGen)))
        @info "Minimum stable generation levels"
        @constraint(m, MinStableGeneration[(n,t)=keys(MinStableGen), h=Hours; !(MinStableGen[n,t] |> ismissing)],
            G[(n,t),h] >= MinStableGen[n,t] * CapacityDerating[n,t,h] * time_ratio * (N_TECH[(n,t)] + CapAdd[:N_TECH][(n,t)])
        );
        ## Average generation version:
        # @constraint(m, MinStableGeneration[(n,t)=Nodes_Techs; !(MinStableGen[n,t] |> ismissing)],
            # sum(G[(n,t),h] for h in Hours) >= MinStableGen[n,t]*length(Hours) )
    end

#=
    @info "Carbon budget upper limit (annual)"
    @constraint(m, CarbonBudgetLimit,
        sum((CarbonContent[n,t]/Efficiency[n,t]) * G[(n,t),h]
                for (n,t) in Nodes_Techs, h in Hours)
        <= CarbonBudget
    );

    next!(prog)
    println("\n")
=#
# %% * ----------------------------------------------------------------------- *
#    ***** Quotas for renewable technologies *****
#    * ----------------------------------------------------------------------- *

    #  Minimum yearly renewables requirement (con5a_minRES)
    @info "Minimum yearly renewables requirement."
    @constraint(m, MinRES[n=DemandRegions],
        sum(
           sum(G[(z,t),h] for (z,t) in Nodes_Renew if z == zone)  # Any renewable TxZone-level tech.
         + sum(G_REZ[rez,h] for (rez, z) in Nodes_Promotes if z == zone)
        for (zone, d) in Nodes_Demand if node2DemReg[d] == n
        for h in Hours
        )
        >=
        (min_res_dict[n]/100)*(
            sum(
               sum(G[(z,t),h] for (z,t) in Nodes_Techs if z == zone)
             + sum(G_REZ[rez,h] for (rez, z) in Nodes_Promotes if z == zone)
            for (zone, d) in Nodes_Demand if node2DemReg[d] == n
            for h in Hours
         )
       )
    );
    # Note: if there is NO renewable energy associated to the DemandRegion, then this constraint will
    # act to zero out ALL generation associated to the DemandRegion.

            #  + sum(STO_OUT[(z,sto),h] - STO_IN[(z,sto),h] for (z,sto) in Nodes_Storages if z == zone)

            #  + sum(H2_G2P[(z,g2p),h] for (z,g2p) in Nodes_G2P if z == zone)
            #  - sum(H2_P2G[(z,p2g),h] for (z,p2g) in Nodes_P2G if z == zone)

    #  Minimum yearly renewables requirement for whole of system.
    @info "Minimum yearly renewables requirement for whole of system."
    @constraint(m, MinRESsystem,
        sum(
            sum(
                sum(G[(z,t),h] for (z,t) in Nodes_Renew if z == zone)  # Any renewable TxZone-level tech.
                + sum(G_REZ[rez,h] for (rez, z) in Nodes_Promotes if z == zone)
                for (zone, d) in Nodes_Demand if node2DemReg[d] == n
            for h in Hours)
        for n in DemandRegions)
            >=
            (MinimumRenewShare/100)*sum(
                sum(
                    sum(G[(z,t),h] for (z,t) in Nodes_Techs if z == zone)
                    + sum(G_REZ[rez,h] for (rez, z) in Nodes_Promotes if z == zone)
                    for (zone, d) in Nodes_Demand if node2DemReg[d] == n
                for h in Hours)
            for n in DemandRegions)
    );

# %% * ----------------------------------------------------------------------- *
#    ***** Storage constraints *****
#    * ----------------------------------------------------------------------- *

    # Storage Level Dynamics Initial Condition (con4a_stolev_start)
    @info "Storage level dynamics: initial condition."
    @constraint(m, StorageLevelStart[(n,sto)=Nodes_Storages],
        STO_L[(n,sto),Hours[1]]
         ==
        StartLevel[n,sto] * (N_STO_E[(n,sto)] + CapAdd[:N_STO_E][(n,sto)])
        +   sqrt(Efficiency[n,sto])*STO_IN[(n,sto), Hours[1]]
        - 1/sqrt(Efficiency[n,sto])*STO_OUT[(n,sto), Hours[1]]
    );

    # End level equal to initial level (con4j_ending)
    @info "Storage end level equal to initial level."
    @constraint(m, StorageLevelEnd[(n,sto)=Nodes_Storages],
        STO_L[(n,sto),Hours[end]]
         ==
        StartLevel[n,sto] * (N_STO_E[(n,sto)] + CapAdd[:N_STO_E][(n,sto)])
    );

    # @constraint(m, StorageBalanceFirstHour[(n,sto) in Nodes_Storages],
    #     STO_L[(n,sto),Hours[1]]
    #     ==
    #     STO_L[(n,sto),Hours[end]]
    #     +   sqrt(Efficiency[n,sto])*STO_IN[(n,sto), Hours[1]]
    #     - 1/sqrt(Efficiency[n,sto])*STO_OUT[(n,sto), Hours[1]]
    #     );

    # "Storage Level Dynamics
    @info "Storage level dynamics."
    @constraint(m, StorageBalance[(n,sto)=Nodes_Storages,h=Hours2],
        STO_L[(n,sto), h]
        ==
        STO_L[(n,sto), h-1]
        +   sqrt(Efficiency[n,sto])*STO_IN[(n,sto), h]
        - (1/sqrt(Efficiency[n,sto]))*STO_OUT[(n,sto), h]
    );

    # Storage Power Capacity (con4c_stolev_max)
    @info "Storage power capacity."
    @constraint(m, MaxLevelStorage[(n,sto)=Nodes_Storages,h=Hours],
        STO_L[(n,sto),h] <= (N_STO_E[(n,sto)] + CapAdd[:N_STO_E][(n,sto)])
    );

    # Storage maximum inflow (con4d_maxin_sto)
    @info "Storage maximum inflow."
    @constraint(m, MaxWithdrawStorage[(n,sto)=Nodes_Storages,h=Hours],
        STO_IN[(n,sto),h] <= time_ratio * (N_STO_P[(n,sto)] + CapAdd[:N_STO_P][(n,sto)])
    );

    # Storage maximum outflow by capacity (con4e_maxout_sto)
    @info "Storage generation outflow by capacity."
    @constraint(m, MaxGenerationStorage[(n,sto)=Nodes_Storages,h=Hours],
        STO_OUT[(n,sto),h] <= time_ratio * (N_STO_P[(n,sto)] + CapAdd[:N_STO_P][(n,sto)])
    );

    if !(isDictAllMissing(dtr.parameters[:MaxEnergy]))
        @info "Storage: maximum energy allowed."
        @constraint(m, MaxEnergyStorage[(n,sto)=Nodes_Storages; !(MaxEnergy[n,sto] |> ismissing)],
            N_STO_E[(n,sto)] <= MaxEnergy[n,sto]
        );
    end

    @info "Storage: maximum power allowed."
    @constraint(m, MaxPowerStorage[(n,sto)=Nodes_Storages; !(MaxCapacity[n,sto] |> ismissing)],
        N_STO_P[(n,sto)] + CapAdd[:N_STO_P][(n,sto)] <= MaxCapacity[n,sto]
    );

    # Maximum Energy to Power ratio for certain storage technologies (con4k_PHS_EtoP)
    @info "Storage: maximum energy-to-power ratio (use time)"
    @constraint(m, EnergyToPowerRatio[(n,sto)=Nodes_Storages; !(MaxEtoP_ratio[n,sto] |> ismissing)],
        N_STO_E[(n,sto)] + CapAdd[:N_STO_E][(n,sto)] <= MaxEtoP_ratio[n,sto] * (N_STO_P[(n,sto)] + CapAdd[:N_STO_P][(n,sto)])
    );

    # Maximum storage outflow - no more than level of last period (con4h_maxout_lev)
    @info "Storage: maximum outflow - no more than level of last period."
    @constraint(m, MaxOutflowStorage[(n,sto)=Nodes_Storages,h=Hours2],
         (1/sqrt(Efficiency[n,sto]))*STO_OUT[(n,sto),h] <= STO_L[(n,sto),h-1]
    );

    # Maximum storage inflow - no more than energy capacity minus level of last period (con4i_maxin_lev)
    @info "Storage: maximum inflow - no more than energy capacity minus level of last period"
    @constraint(m, MaxInflowStorage[(n,sto)=Nodes_Storages,h=Hours2],
         sqrt(Efficiency[n,sto])*STO_IN[(n,sto),h] <= N_STO_E[(n,sto)] + CapAdd[:N_STO_E][(n,sto)] - STO_L[(n,sto),h-1]
    );

    next!(prog)
    println("\n")

    return dtr
end

function build_ev_constraints(dtr::DieterModel)
# %% * ----------------------------------------------------------------------- *
#    ***** Electric vehicle constraints *****
#    * ----------------------------------------------------------------------- *

    # Sets:
    Hours = dtr.sets[:Hours]
    Hours2 = Hours[2:end]

    EV = dtr.sets[:ElectricVehicles]

    # Parameters:
    EvDemand = dtr.parameters[:AbsoluteEvDemand]
    EvPower = dtr.parameters[:AbsoluteEvPower]
    EvCapacity = dtr.parameters[:AbsoluteEvCapacity]
    EvFuel = dtr.parameters[:EvFuel]
    EvType = dtr.parameters[:EvType]

    # Variables:
    EV_CHARGE = dtr.model.obj_dict[:EV_CHARGE]
    EV_DISCHARGE = dtr.model.obj_dict[:EV_DISCHARGE]
    EV_L = dtr.model.obj_dict[:EV_L]
    EV_PHEVFUEL = dtr.model.obj_dict[:EV_PHEVFUEL]
    EV_INF = dtr.model.obj_dict[:EV_INF]

    @info "Electric Vehicles: electric charge allowable."
    @constraint(m, MaxWithdrawEV[ev=EV,h=Hours],
        EV_CHARGE[ev,h] <= EvPower[ev][h]
    );

    @info "Electric Vehicles: electric discharge allowable."
    @constraint(m, MaxGenerationEV[ev=EV,h=Hours],
        EV_DISCHARGE[ev,h] <= EvPower[ev][h]
    );

    @info "Electric Vehicles: charge level limited by capacity."
    @constraint(m, MaxLevelEV[ev=EV,h=Hours],
        EV_L[ev,h] <= EvCapacity[ev]
    );

    @info "Electric Vehicles: energy balance."
    @constraint(m, EVBalance[ev=EV,h=Hours2],
        EV_L[ev, h]
        ==
        EV_L[ev, h-1]
        + EV_CHARGE[ev, h]*sqrt(Efficiency[ev])
        - EV_DISCHARGE[ev, h]/sqrt(Efficiency[ev])
        - EvDemand[ev][h]
        + EvType[ev]*EV_PHEVFUEL[ev, h]
        + EV_INF[ev, h]
    );

    @info "Electric Vehicles: energy balance at first and last time-steps."
    @constraint(m, EVBalanceFirstHour[ev=EV],
        EV_L[ev,Hours[end]]    # or EV_L[ev,Hours[1]]
        ==
        EV_L[ev,Hours[1]]  # or EV_L[ev,Hours[end]]
        + EV_CHARGE[ev,Hours[end]]*sqrt(Efficiency[ev])
        - EV_DISCHARGE[ev,Hours[end]]/sqrt(Efficiency[ev])
        - EvDemand[ev][Hours[end]]
        + EvType[ev]*EV_PHEVFUEL[ev, Hours[end]]
        + EV_INF[ev, Hours[end]]
    );

    # next!(prog)
    # println("\n")

    return dtr
end

function build_h2_constraints(dtr::DieterModel)
# %% * ----------------------------------------------------------------------- *
#    ***** Hydrogen constraints *****
#    * ----------------------------------------------------------------------- *

    Timestep = dtr.settings[:timestep]
    periods = round(Int,hoursInYear*(2/Timestep))

    # Sets:
    Hours = dtr.sets[:Hours]
    Hours2 = Hours[2:end]

    HoursToMonths = dtr.parameters[:HoursToMonths]

    Nodes_P2G = dtr.sets[:Nodes_P2G]
    Nodes_P2G_PEM = dtr.sets[:Nodes_P2G_PEM]
    Nodes_P2G_AE = dtr.sets[:Nodes_P2G_AE]

    Nodes_G2P = dtr.sets[:Nodes_G2P]
    Nodes_GasStorages = dtr.sets[:Nodes_GasStorages]

    # Parameters:
    time_ratio = dtr.settings[:time_ratio]

    H2Conversion = dtr.parameters[:H2Conversion]  # Units: MWh / tonne-H2 for a given P2G tech.
    # H2Demand = coalesce((dtr.settings[:h2]*1e6)/hoursInYear,0)
    H2Demand = dtr.parameters[:H2Demand] # Units: MWh / year for a given Node and P2G tech.

    CapAdd = dtr.parameters[:CapAdd]

    Capacity_Factor_AE = dtr.settings[:capacity_factor_ae] # e.g 0.8 means must produce at least 80% of hourly quota each hour.

    # if dtr.settings[:H2_G2P_Allowed_flag] == true
    RecipH2_CF_UB = dtr.settings[:capacity_factor_recipH2_ub]
    # end

    # Variables:
    H2_P2G = dtr.model.obj_dict[:H2_P2G]
    H2_G2P = dtr.model.obj_dict[:H2_G2P]
    H2_GS_L = dtr.model.obj_dict[:H2_GS_L]
    H2_GS_IN = dtr.model.obj_dict[:H2_GS_IN]
    H2_GS_OUT = dtr.model.obj_dict[:H2_GS_OUT]
    N_P2G = dtr.model.obj_dict[:N_P2G]
    N_G2P = dtr.model.obj_dict[:N_G2P]
    N_GS = dtr.model.obj_dict[:N_GS]

    @info "Hydrogen: Minimum yearly lower bound on power-to-gas."
    @constraint(dtr.model, MinYearlyP2G[(n,p2g)=Nodes_P2G],
        sum(H2_P2G[(n,p2g),h] for h in Hours) >= H2Demand[n,p2g]
    );

    # @info "Hydrogen: Minimum monthly lower bound on power-to-gas for PEM tech."
    # @constraint(dtr.model, MinMonthlyP2G_PEM[(n,p2g)=Nodes_P2G_PEM, mth=1:12],
    #     sum(H2_P2G[(n,p2g),h] for h in Hours if HoursToMonths[h] == mth)
    #         >= (H2Demand[n,p2g]/periods)*sum(1.0 for h in Hours if HoursToMonths[h] == mth)
    # );

    @info "Hydrogen: Constant hourly lower bound on power-to-gas for AE tech."
    @constraint(dtr.model, MinHourlyP2G_AE[(n,p2g)=Nodes_P2G_AE, h=Hours],
        H2_P2G[(n,p2g),h] >= Capacity_Factor_AE*(H2Demand[n,p2g]/periods)
    );

    @info "Hydrogen: Variable upper bound on power-to-gas."
    @constraint(dtr.model, MaxP2G[(n,p2g)=Nodes_P2G,h=Hours],
        H2_P2G[(n,p2g),h] <= time_ratio * (N_P2G[(n,p2g)] + CapAdd[:N_P2G][(n,p2g)])
    );

    @info "Hydrogen: Variable upper bound on gas-to-power."
    @constraint(dtr.model, MaxG2P[(n,g2p)=Nodes_G2P,h=Hours],
        H2_G2P[(n,g2p),h] <= time_ratio * (N_G2P[(n,g2p)] + CapAdd[:N_G2P][(n,g2p)])
    );

    @constraint(dtr.model, MaxCapFactorH2_G2P[(n,g2p)=Nodes_G2P],
        sum(H2_G2P[(n,g2p),h] for h in Hours) <= RecipH2_CF_UB*length(Hours)*(N_G2P[(n,g2p)] + CapAdd[:N_G2P][(n,g2p)])
    );

    @info "Hydrogen: Variable upper bound on gas storage."
    @constraint(dtr.model, MaxLevelGasStorage[(n,gs)=Nodes_GasStorages,h=Hours],
        H2_GS_L[(n,gs),h] <= N_GS[(n,gs)]  # + CapAdd[:N_GS][(n,p2g)] )
    );

    @constraint(dtr.model, GasStorageIn[(n,gs)=Nodes_GasStorages,h=Hours],
        H2_GS_IN[(n,gs),h] == sum(H2_P2G[(n,p2g),h]/H2Conversion[n,p2g] for (n,p2g) in Nodes_P2G)
    );

    @constraint(dtr.model, GasStorageOut[(n,gs)=Nodes_GasStorages,h=Hours],
        H2_GS_OUT[(n,gs),h] == sum(H2_G2P[(n,g2p),h]/H2Conversion[n,g2p] for (n,g2p) in Nodes_G2P)
    );

#=
    @info "Hydrogen: energy balance."
    @constraint(dtr.model, H2Balance[h=Hours],
        sum(Efficiency[p2g]*H2_P2G[p2g,h] for p2g in Nodes_P2G)  # No sqrt on Efficiency, not a roundtrip
        + sum(H2_GS_IN[gs,h] for gs in Nodes_GasStorages)   ## Changed: H2_GS_OUT -> H2_GS_IN
        # + sum(H2_GS_OUT[gs,h] for gs in GasStorages)   ## Should this be changed: H2_GS_OUT -> H2_GS_IN ?
        ==
        sum((1/Efficiency[g2p])*H2_G2P[g2p,h] for g2p in Nodes_G2P)   # No sqrt on Efficiency, not a roundtrip
        + sum(H2_GS_OUT[gs,h] for gs in Nodes_GasStorages)    ## Changed: H2_GS_IN -> H2_GS_OUT
        # + sum(H2_GS_IN[gs,h] for gs in GasStorages)    ## Should this be changed: H2_GS_IN -> H2_GS_OUT ?
        + H2Demand
    );
=#
    @info "Hydrogen: gas storage mass balance."
    @constraint(dtr.model, GasStorageBalance[(n,gs)=Nodes_GasStorages,h=Hours2],
        H2_GS_L[(n,gs), h] == H2_GS_L[(n,gs), h-1] + H2_GS_IN[(n,gs), h] - H2_GS_OUT[(n,gs), h]
    );

    @info "Hydrogen: gas storage balance at first time-steps."
    @constraint(dtr.model, GasStorageBalanceFirstHours[(n,gs)=Nodes_GasStorages],
        H2_GS_L[(n,gs), Hours[1]] == StartLevel[n,gs] * N_GS[(n,gs)] + H2_GS_IN[(n,gs),Hours[1]] - H2_GS_OUT[(n,gs),Hours[1]]  # + CapAdd[:N_GS][(n,p2g)] )
    );

    # End level equal to initial level
    @info "Hydrogen: gas storage end level equal to initial level."
    @constraint(dtr.model, GasStorageLevelEnd[(n,gs)=Nodes_GasStorages],
        H2_GS_L[(n,gs), Hours[end]] == StartLevel[n,gs] * N_GS[(n,gs)]  # + CapAdd[:N_GS][(n,p2g)] )
    );

    # next!(prog)
    # println("\n")

    return dtr
end

function build_heat_load_constraints(dtr::DieterModel)
# %% * ----------------------------------------------------------------------- *
#    ***** Heat consumption constraints *****
#    * ----------------------------------------------------------------------- *
    # Sets:
    Hours = dtr.sets[:Hours]
    Hours2 = Hours[2:end]

    BU = dtr.sets[:BuildingType]  # Building archtypes
    HP = dtr.sets[:HeatingType]   # Heating combination type

    # Parameters:
    time_ratio = dtr.settings[:time_ratio]

    HeatConsumption = dtr.parameters[:HeatConsumption]  # Units: MWh per m^2 per hour
    HeatShare = dtr.parameters[:HeatShare]
    MaxLevel = dtr.parameters[:MaxLevel]
    HeatMaxPower = dtr.parameters[:HeatMaxPower]
    StaticEfficiency = dtr.parameters[:StaticEfficiency]
    CoP = dtr.parameters[:CoP]

    # Variables:
    HEAT_STO_L = dtr.model.obj_dict[:HEAT_STO_L]
    HEAT_HP_IN = dtr.model.obj_dict[:HEAT_HP_IN]
    HEAT_INF = dtr.model.obj_dict[:HEAT_INF]

    @constraint(m, HeatBalance[bu=BU, hp=HP, h=Hours2],
        HEAT_STO_L[bu,hp,h]
        ==
        StaticEfficiency[hp] * HEAT_STO_L[bu,hp,h-1]  # h-timestep replaced with h-1
        + CoP[bu,hp][h] * HEAT_HP_IN[bu,hp,h]
        - HeatConsumption[bu,hp][h]
        + HEAT_INF[bu,hp,h]
    );

    @constraint(m, HeatBalanceFirstHour[bu=BU, hp=HP],
        HEAT_STO_L[bu,hp,Hours[1]]
        ==
        StaticEfficiency[hp] * HEAT_STO_L[bu,hp,Hours[end]]
        + CoP[bu,hp][Hours[1]] * HEAT_HP_IN[bu,hp,Hours[1]]
        - HeatConsumption[bu,hp][Hours[1]]
        + HEAT_INF[bu,hp,Hours[1]]
    );

    @constraint(m, MaxHeatLevel[bu=BU, hp=HP, h=Hours],
        HEAT_STO_L[bu,hp,h] <= MaxLevel[bu]
    );

    @constraint(m, MaxHeatPower[bu=BU, hp=HP, h=Hours],
        HEAT_HP_IN[bu,hp,h] <= time_ratio * HeatMaxPower[bu,hp]
    );

    # next!(prog)

    return dtr
end

function solve_model!(dtr::DieterModel,solver)

    @info "Starting optimization..."
    JuMP.optimize!(dtr.model)

    return dtr
end

function generate_results!(dtr::DieterModel,varNamesDict::Dict{Symbol,Array{Symbol,1}})
    @info "Storing results"

    mo = dtr.model.obj_dict

    dtr.results = Dict(v[1] => convert_jump_container_to_df(mo[v[1]], dim_names=v[2]) for v in varNamesDict)
            # value_col=v[1])
    # dtr.results = [convert_jump_container_to_df(value.(model_dict[v[1]]), dim_names=v[2]) for v in varNamesDict]

    # if abs(sum(dtr.results[:G_INF][!, :Value]))  > (1e-5) ||
    #    ( !ismissing(dtr.settings[:ev]) &&
    #    abs(sum(dtr.results[:EV_INF][!, :Value])) > (1e-5) )
    #       @warn "Problem might be infeasable"
    # end

    return nothing
end

ext: jl


text:
module GdOpt

import LsoBase
import Obj
import Obj.Objective
import Ls
import Ls.LineSearch

"""
    Arbitrary state of optimizer
"""
abstract State

"""
    State for stateless optimizers
"""
type EmptyState <: State
end

"""
    GdOptimizer(name, init, update)

      name: The name of the optimization method
      init: The optimizers initial state.
    update: (obj, k, w, b, inf, state) -> (f(w), ∇ f(w), s, state)
"""
type GdOptimizer
      name::ASCIIString
      init::State
    update::Function
end

"""
    gd(obj, w [, ls; ϵ, maxiter, storeiter, maxtime])

    Performs Steepest Descent on objective function with initial w and the
    given Line Search function. Returns info DataFrame.
"""
@fastmath function opt(optimizer::GdOptimizer, ls::LineSearch, obj::Objective, w_0::Array{Float64,1};
                         ϵ::Float64=1e-6, maxiter::Int=typemax(Int), maxtime::Float64=60.0, batchsize::Int=-1)
    inf = LsoBase.new_inf()

    # print info header
    headline = @sprintf "\n%6s | %6s | %3s | %9s | %9s"  "k" "sec" "i" "f" "infeas"
    println(headline, "\n", repeat("-", length(headline)))

    # optimization
    lsiter      = 0
    start       = Base.time()
    storagetime = 0.0
    printtime   = 0.0
    minopt      = Inf
    maxf        = -Inf
    state       = optimizer.init
    try

        for k = 1:maxiter

            b = Int[]
            if batchsize > 0
                b = Obj.randbatch(obj, batchsize)   # random sbt index batch
            end
            fw, gw, s, state = optimizer.update(obj, k, w_0, b, inf, state)

            # obtain opt, push info to array
            opt = vecnorm(gw, Inf)

            # store and print info
            time   = Base.time() - start
            maxf   = max(fw,  maxf)
            minopt = min(opt, minopt)
            if time - storagetime > .1 || k == 1 || time > maxtime || opt < ϵ
                LsoBase.push_inf!(inf, w_0, maxf, minopt, k-1, lsiter, time)
                storagetime = floor(time*10) / 10
                if time - printtime > 1.0 || k == 1 || time > maxtime
                    println(@sprintf "%6d | %6.3f | %3d | %9.3e | %9.3e"  k-1 time lsiter maxf minopt)
                    printtime = floor(time)
                end
                minopt = Inf
                maxf   = -Inf
            end

            # take step or stop
            if time > maxtime || opt < ϵ    # stopping criterion satisfied?
                break
            else
                α, lsiter = Ls.ls(ls, w_0, s, obj, b, fw, gw)
                w_0 += α*s
            end

        end # end of optimization

    catch e

        if isa(e, InterruptException)
            println("Optimization aborted due to InterruptException")
        else
            throw(e)
        end

    end

    return inf
end

# include ./opt/*
LsoBase.includedir(dirname(@__FILE__)*"/gdopt")

end

max_stars_repo_name: mirkobunse/lso

max_issues_repo_name: mirkobunse/lso

max_forks_repo_name: mirkobunse/lso


text:
function SerializeArray(array::Array{Float64,2})
    out=zeros(size(array)[1]*3)
    for i in 0:size(array)[1]-1
        out[i*3+1]=array[i+1,1]
        out[i*3+2]=array[i+1,2]
        out[i*3+3]=array[i+1,3]
    end
    return out
end

function SerializeArray(array::Array{Array{Float64,1}})
    out=zeros(length(array)*3)
    for i in 0:length(array)-1
        out[i*3+1]=array[i+1][1]
        out[i*3+2]=array[i+1][2]
        out[i*3+3]=array[i+1][3]
    end
    return out
end

function SerializeArray(array::Array{Array{Int64,1}})
    out=zeros(Int64,length(array)*3)
    for i in 0:length(array)-1
        out[i*3+1]=array[i+1][1]
        out[i*3+2]=array[i+1][2]
        out[i*3+3]=array[i+1][3]
    end
    return out
end

function SerializeArray(array::Array{Int64,2})
    out=zeros(Int64,size(array)[1]*3)
    for i in 0:size(array)[1]-1
        out[i*3+1]=array[i+1,1]
        out[i*3+2]=array[i+1,2]
        out[i*3+3]=array[i+1,3]
    end
    return out
end

function SerializeArray(array::Array{Int64,1})
    out=zeros(Int64,size(array)[1])
    for i in 0:size(array)[1]-1
        out[i+1]=array[i+1]
    end
    return out
end

function SerializeArray(array::Array{Float64,1})
    out=zeros(size(array)[1])
    for i in 0:size(array)[1]-1
        out[i+1]=array[i+1]
    end
    return out
end

function DeSerializeArray(array::Nothing)
end

function DeSerializeArray(array::Array{Float64,1})
    out=zeros(div(size(array)[1],3),3)
    for i in 0:div(size(array)[1],3)-1
        out[i+1,1]=array[i*3+1]
        out[i+1,2]=array[i*3+2]
        out[i+1,3]=array[i*3+3]
    end
    return out
end

function DeSerializeArray(array::Array{Int64,1})
    out=zeros(Int64,div(size(array)[1],3),3)
    for i in 0:div(size(array)[1],3)-1
        out[i+1,1]=array[i*3+1]
        out[i+1,2]=array[i*3+2]
        out[i+1,3]=array[i*3+3]
    end
    return out
end

function DeSerializeBondArray(array::Array{Int64,1})
    out=zeros(Int64,div(size(array)[1],3),3)
    for i in 0:div(size(array)[1],3)-1
        out[i+1,1]=array[i*3+1]
        out[i+1,2]=array[i*3+2]
        out[i+1,3]=array[i*3+3]
    end
    return out
end

function first_fit(items::Array{Int64,1},bincount::Int64)
    sort_idex = sortperm(items,rev=true)
    items = items[sort_idex] # New - improves first fit.
    bins     = [Int64[] for c =1:bincount]
    index_bins     = [Int64[] for c =1:bincount]
    binsizes = [0 for c =1:bincount]
    for (index,item) in enumerate(items)
        minbinindex = findall(x->x==minimum(binsizes),binsizes)[1]
        push!(bins[minbinindex],item)
        push!(index_bins[minbinindex],sort_idex[index])
        binsizes[minbinindex] += item
    end
    average = sum(binsizes) / bincount
    maxdeviation = maximum([abs(average - bs) for bs in binsizes])
    return bins, index_bins,binsizes, average, maxdeviation
end

function greedy_swap(columns, colsize, average, margin)
    # See if you can do a swap to smooth the heights
    a=0
end

# Target: distribute atoms on tasks as even as possible
struct BalancingMPI
    comms::Array{Array{Int64,1},1}
    mole_index_comms::Array{Array{Int64,1},1}
    commsizes::Array{Int64,1}
    average::Float64
    maxdeviation::Float64
    bond_commsizes::Array{Int64,1}
    moles_length::Array{Int64,1}
    bond_moles_length::Array{Int64,1}
    function BalancingMPI(config::Configuration,comm_size::Int64)
        comms, mole_index_comms,commsizes, average, maxdeviation=first_fit(config.moles_length,comm_size)
        bond_moles_length=zeros(Int64,length(config.moles))
        bond_commsizes=zeros(Int64,comm_size)
        if config.moles_length!=zeros(Int64,0)
            @info "CPU Loading Information: CPU_index Num_Atoms"
            for (index,tasksize_i) in enumerate(commsizes)
                @info "CPU Loading Information: $(index) $(tasksize_i)"
            end
            @info "CPU Loading Average: $(average) Atoms with Deviation: $(maxdeviation)"
            for (index,moles_comm) in enumerate(mole_index_comms)
                n_bonds=0
                for mole_i in moles_comm
                    n_bonds+=size(config.moles[mole_i].bonds)[1]
                    bond_moles_length[mole_i]=size(config.moles[mole_i].bonds)[1]
                end
                bond_commsizes[index]=n_bonds
            end
        end

        new(
            comms::Array{Array{Int64,1},1},
            mole_index_comms::Array{Array{Int64,1},1},
            commsizes::Array{Int64,1},
            average::Float64,
            maxdeviation::Float64,
            bond_commsizes::Array{Int64,1},
            config.moles_length::Array{Int64,1},
            bond_moles_length::Array{Int64,1}
        )
    end
end


text:
using HyperbolicPlane, Plots

"""
`concentric()` draws a picture of several concentric circles.
"""
function concentric()
    Plots.plot()
    P = HPoint(1,pi/6)
    for r in 0.25:0.25:3
        C = HCircle(P,r)
        set_thickness(C, 1.5)
        draw(C)
    end
    draw(HPlane())
    finish()
end

concentric()

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
module SpatioTemporalTraits

using ArrayInterface
using ArrayInterface: dimnames, has_dimnames, to_dims, parent_type, to_parent_dims
using Dates
using Metadata
using Metadata: MetaArray, NoMetadata
using Static
using Base: tail

using LinearAlgebra

export
    assert_timedim_last,
    duration,
    has_timedim,
    ntimes,
    origin,
    pixel_spacing,
    sampling_rate,
    sdims,
    spatialdims,
    spatial_order,
    spatial_size,
    spatial_indices,
    spatial_first,
    spatial_last,
    spatial_directions,
    timedim,
    times,
    time_first,
    time_last,
    time_step,
    width,
    depth,
    height

include("spatial.jl")
include("temporal.jl")

end # module


text:
# As with the first half of this assignment, `learning_julia1_exercises.jl`,
# fill in the missing implementation and pass the `@test`s.
#
# If you're writing your code within the `@testset`, the functions will not be available outside it.
# (A `@testset` defines an enclosed scope that doesn't "pollute" the outside.)
# You may find this makes it harder to develop and debug your code.
# Remember that in VS Code you can select a range of lines and hit Ctrl-Enter to evaluate it,
# and that may make it easier to experiment in VS Code's Julia REPL.
#
# We're also extending your usage of Test by introducing `@test_throws`.
# Read its help to see what it does. When you pass argument(s) to the exception type, you check that, e.g.,
# a particular message will be passed to the user.
# Starting in Julia 1.8 it will be possible to check the error message without being concerned about
# the exception type. (https://github.com/JuliaLang/julia/pull/41888)

using Test

@testset "Learning Julia 2" begin
    # Problem 1 [Adapted from ThinkJulia, exercise 3-2]:
    # Create a function `rightjustify` that prints a string with the right number of leading spaces
    # so that the final character aligns with a chosen column.
    # I've given you the signature of one method; also create another method that supplies `stdout`
    # as the default `io` choice. See the style guide for more info about recommended best
    # practices for ordering arguments for Julia functions:
    # https://docs.julialang.org/en/v1/manual/style-guide/#Write-functions-with-argument-ordering-similar-to-Julia-Base
    function rightjustify(io::IO, str::AbstractString, col::Integer = 70)
        # fill me in!
    end
    # also define a method that doesn't have `io` supplied

    # here are my tests (don't modify this code)
    checkstr(ostr, refstr, len) = @test length(ostr) == len && endswith(ostr, refstr) && strip(ostr) == refstr
    for str in ("short", "kind of medium", "getting longer, but not too huge")
        checkstr(sprint(rightjustify, str),     str, 70)
        checkstr(sprint(rightjustify, str, 80), str, 80)
    end
    mktemp() do fname, io
        str = "send me to stdout!"
        redirect_stdout(io) do
            rightjustify(str)
        end
        seek(io, 0)
        checkstr(chomp(read(io, String)), str, 70)
    end

    # Problem 2: adapted from https://exercism.org/tracks/julia/exercises/robot-simulator
    # (Exercism is a site offering programming exercises in many different languages, with a nice
    #  framework for in-browser coding, automated testing, progress-monitoring, and even mentoring.)
    # The version of this problem on Exercism encourages you to develop a custom type, which is great,
    # but here we're going to treat it as an exercise in control-flow.
    # If you wish, you can put everything in a single function.
    #
    # A robot starts at a location `(x, y)`, where both `x` and `y` are integers, facing in a specified direction
    # ('N'=North, 'S'=South, 'E'=East, 'W'=West). `x` is the East-West axis and `y` the North-South axis.
    # You also supply a string, where each character of the string is one of the following:
    # - 'A': advance one step in the direction it's facing
    # - 'R': turn right
    # - 'L': turn left
    # When the robot turns, it pivots without changing its location. For example, "RAALAL" means:
    # 1. turn right
    # 2. advance two spaces
    # 3. turn left
    # 4. advance one space
    # 5. turn left
    #
    # Your code should return both the final position and the direction the robot is facing.
    # See if you can come up with a more efficient approach than a huge nest of `if...else...end`
    # statements. My solution is less than 20 lines of code.
    # Hint: some of the better solutions might use `%` or `mod1` (read their help strings).

    function robot_command(cmd, (x, y), dir)  # argument destructuring syntax: https://docs.julialang.org/en/v1/manual/functions/#Argument-destructuring
        # you write the body
    end

    # Single steps (these are deliberately written out in long form to avoid giving too much of the solution away)
    @test robot_command("R", (0, 0), 'N') == ((0, 0), 'E')
    @test robot_command("R", (0, 0), 'E') == ((0, 0), 'S')
    @test robot_command("R", (0, 0), 'S') == ((0, 0), 'W')
    @test robot_command("R", (0, 0), 'W') == ((0, 0), 'N')
    @test robot_command("L", (0, 0), 'N') == ((0, 0), 'W')
    @test robot_command("L", (0, 0), 'E') == ((0, 0), 'N')
    @test robot_command("L", (0, 0), 'S') == ((0, 0), 'E')
    @test robot_command("L", (0, 0), 'W') == ((0, 0), 'S')
    @test robot_command("A", (0, 0), 'N') == ((0, 1), 'N')
    @test robot_command("A", (0, 0), 'E') == ((1, 0), 'E')
    @test robot_command("A", (0, 0), 'S') == ((0, -1), 'S')
    @test robot_command("A", (0, 0), 'W') == ((-1, 0), 'W')
    @test_throws ErrorException("robot command K not recognized") robot_command("K", (0, 0), 'W')
    # Sequences
    for x in (7, 22), y in (3, -7)
        @test robot_command("RAALAL", (x, y), 'N') == ((x+2, y+1), 'W')
        @test robot_command("ARAAALA", (x, y), 'W') == ((x-2, y+3), 'W')
        for dir in ('N', 'E', 'S', 'W')
            @test robot_command("ALALALAL", (x, y), dir) == ((x, y), dir)
        end
    end

    # Problem 3: introduction to structures
    # Note that when you redefine structures, you may have to shut down your Julia session and restart it.
    # (Julia won't let you change an existing `struct` definition.)
    # a) create a new subtype of `Integer` called `MyInt`, one that stores an `Int` internally

    @test MyInt(3) isa Integer
    @test !ismutable(MyInt(3))

    # b) extend `Base`'s `==` so that you can compare `MyInt` against any other `Integer`s
    #    You can either `import Base: ==` and write `==(a, b)` methods, or write your new
    #    methods using module qualification, `Base.:(==)(a, b)`. The `:` is needed only for operators;
    #    when extending other functions, you don't need `:`.

    @test MyInt(3) == 3
    @test 2 == MyInt(2)
    @test MyInt(1) != MyInt(2)   # this will follow automatically if you define `==`

    # c) extend `+` so that `MyInt` "wraps around" at 5: `MyInt(5) + 1` return `MyInt(1)` and then it starts counting up again.

    @test MyInt(1) + 1  === MyInt(2)
    @test MyInt(2) + 1  === MyInt(3)
    @test MyInt(3) + 1  === MyInt(4)
    @test MyInt(4) + 1  === MyInt(5)
    @test MyInt(5) + 1  === MyInt(1)
    @test MyInt(4) + 2  === MyInt(1)
    @test MyInt(3) + 10 === MyInt(3)

    # d) also extend the `typemax` "trait" for `MyInt`. Read the docs for `typemax`, and then do `@edit typemax(UInt8)`
    #    to see how such methods are implemented. Mimic that definition.

    @test typemax(MyInt)    === MyInt(5)
    @test typemax(MyInt(2)) === MyInt(5)   # note you didn't have to write a method specifically for this test
    @test typemax(UInt8)    === 0xff       # just to make sure it's not broken

    # Problem 4: encapsulating data in objects & functions
    # a) Create a parametric `struct` representing a unit Gaussian in 1d: if `g isa Gaussian`, then
    #    - `g.μ` should return the mean (center) of the Gaussian, a scalar of type `T`
    #    - `g.σ` should return the standard deviation, a scalar also of type `T`
    #    Also, `T` should only be allowed to be a subtype of `Real`

    g = Gaussian(1.0, 0.5)
    @test g isa Gaussian{Float64}
    @test g.μ === 1.0
    @test g.σ === 0.5
    g = Gaussian(2, 5)
    @test g isa Gaussian{Int}
    @test g.μ === 2
    @test g.σ === 5
    @test_throws TypeError Gaussian{String}

    # b) make objects of type Gaussian callable (review https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects
    #    as needed), so that `g(x)` returns `exp(-(x - μ)^2/(2*σ^2))`.

    g = Gaussian(1.0, 0.5)
    @test g(2) ≈ 0.1353352832366127
    @test g(2.0) == g(2)
    @test g(1.25) ≈ 0.8824969025845955
    g = Gaussian(2, 5)
    @test g(-3) ≈ 0.6065306597126334

    # c) now do the same thing with a *closure*, a function that wraps data.  Here's a demo:
    """
        f = adder(x)

    Return a function `f` that adds `x` to any number.
    """
    function adder(x)
        return y::Number -> x+y
    end
    f = adder(3)
    @test f isa Function
    @test f(2) == 5
    @test_throws MethodError f("hello")   # f isn't defined on strings because of the `y::Number`

    # Another way you could have written the above is
    # function adder(x)
    #     f = function(y::Number)
    #         return x + y
    #     end
    #     return f
    # end
    # If that's clearer to you, feel free to use that form.

    # Now you write one: `f = gaussian_closure(μ, σ)` should return a function `f` such that
    # `f(x) == exp(-(x - μ)^2/(2*σ^2))`. Enforce the fact that both `μ` and `σ` must be `Real`,
    # as must `x`, but don't require them to be the same type.

    f = gaussian_closure(2, 5)
    @test f isa Function
    @test f(-3) ≈ 0.6065306597126334
    @test_throws MethodError f("hello")
    @test_throws MethodError gaussian_closure("hello", 5)
    @test_throws MethodError gaussian_closure(2, [5])
end

max_stars_repo_name: HsupoLeng/AdvancedScientificComputing

max_issues_repo_name: HsupoLeng/AdvancedScientificComputing

max_forks_repo_name: HsupoLeng/AdvancedScientificComputing


text:
using CautiousEngine
using BSON

@testset "serialization and deserialization test" begin

	region_dict = Dict()
	region_post_dict = Dict()
	N_regions = 10000

	for i=1:N_regions-1
		region_dict[i] = Dict("x1" => [-0.1, 0.], "x2" => [-0.1, 0.])
		region_post_dict[i] = (Dict("x1" => [0.1, 0.2], "x2" => [0.3, 0.4]), [0.5, 0.6])
	end
	region_dict[N_regions] = Dict("x1" => [-0.1, 0.], "x2" => [-0.1, 0.])

	region_array = CautiousEngine.serialize_region_data(region_dict, region_post_dict)
	test_filename = "test.bson"
	bson(test_filename, Dict(:res => region_array))

	r_region_array = BSON.load(test_filename)
	@test r_region_array[:res] == region_array

	etime_new = @elapsed CautiousEngine.save_region_data(region_dict, region_post_dict, test_filename)
	etime_old = @elapsed bson("foo.bson", Dict(:regions => region_dict, :region_posts => region_post_dict))

	println("new save time: ", etime_new)
	println("old save time: ", etime_old)

	r_region_dict, r_region_post_dict, _ = CautiousEngine.deserialize_region_data(test_filename)

	@test r_region_dict == region_dict
	@test r_region_post_dict == region_post_dict
end


text:
# ---
# title: 177. Nth Highest Salary
# id: problem177
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories:
# link: <https://leetcode.com/problems/nth-highest-salary/description/>
# hidden: true
# ---
#
# Write a SQL query to get the _n_ th highest salary from the `Employee` table.
#
#
#
#     +----+--------+
#     | Id | Salary |
#     +----+--------+
#     | 1  | 100    |
#     | 2  | 200    |
#     | 3  | 300    |
#     +----+--------+
#
#
# For example, given the above Employee table, the _n_ th highest salary where
# _n_ = 2 is `200`. If there is no _n_ th highest salary, then the query should
# return `null`.
#
#
#
#     +------------------------+
#     | getNthHighestSalary(2) |
#     +------------------------+
#     | 200                    |
#     +------------------------+
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Distributions
using Optim
using HTTP
using GLM
using LinearAlgebra
using Random
using Statistics
using DataFrames
using CSV
using FreqTables

 #:::::::::::::::::::::::::::::::::::::::::::::::::::
 # question 1
 #:::::::::::::::::::::::::::::::::::::::::::::::::::

url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2021/master/ProblemSets/PS4-mixture/nlsw88t.csv"
df = CSV.read(HTTP.get(url).body,DataFrame)
X = [df.age df.white df.collgrad]
Z = hcat(df.elnwage1, df.elnwage2, df.elnwage3, df.elnwage4,
df.elnwage5, df.elnwage6, df.elnwage7, df.elnwage8)
y = df.occ_code

function mlogit_with_Z(theta, X, Z, y)

    alpha = theta[1:end-1]
    gamma = theta[end]
    K = size(X,2)
    J = length(unique(y))
    N = length(y)
    bigY = zeros(N,J)
    for j=1:J
        bigY[:,j] = y.==j
    end
    bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

    T = promote_type(eltype(X),eltype(theta))
    num   = zeros(T,N,J)
    dem   = zeros(T,N)
    for j=1:J
        num[:,j] = exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)
        dem .+= num[:,j]
    end

    P = num./repeat(dem,1,J)

    loglike = -sum( bigY.*log.(P) )

    return loglike
end
startvals = [2*rand(7*size(X,2)).-1; .1]
td = TwiceDifferentiable(theta -> mlogit_with_Z(theta, X, Z, y), startvals; autodiff = :forward)
# run the optimizer
theta_hat_optim_ad = optimize(td, startvals, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
theta_hat_mle_ad = theta_hat_optim_ad.minimizer
# evaluate the Hessian at the estimates
H  = Optim.hessian!(td, theta_hat_mle_ad)
theta_hat_mle_ad_se = sqrt.(diag(inv(H)))
println([theta_hat_mle_ad theta_hat_mle_ad_se]) # these standard errors match Stata

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 2
#:::::::::::::::::::::::::::::::::::::::::::::::::::

#Yes. Estimated gamma in PS3 is -0.09419383447879262. Yet gamma/100 is the change in utility with a 1% increase in expected wage.
#It should not be negative. Here I got a postitive number(1.3074804109794527) so it is more reasonable.

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 3
#:::::::::::::::::::::::::::::::::::::::::::::::::::

#(a)
include("lgwt.jl") # make sure the function gets read in
# define distribution
d = Normal(0,1) # mean=0, standard deviation=1
# get quadrature nodes and weights for 7 grid points
nodes, weights = lgwt(7,-4,4)
# now compute the integral over the density and verify it
println(sum(weights.*pdf.(d,nodes)))
# now compute the expectation and verify it
println(sum(weights.*nodes.*pdf.(d,nodes)))

#(b)
#(1)
d = Normal(0,2)
nodes, weights = lgwt(7,-5,5)
println(sum(weights.*(nodes.^2).*pdf.(d,nodes)))
#(2)
d = Normal(0,2)
nodes, weights = lgwt(10,-5,5)
println(sum(weights.*(nodes.^2).*pdf.(d,nodes)))
#(3)
#As the number of quadrature points increase, the quadrature approximates closer to the true value.

#(c)
#(1)
d = Normal(0,2)
a = -5*2
b = 5*2
z = Uniform(a,b)
x = rand(z,1000000)
println((b-a)*mean((x.^2).*pdf.(d,x)))
#(2)
println((b-a)*mean(x.*pdf.(d,x)))
#(3)
println((b-a)*mean(pdf.(d,x)))
#(4)
x = rand(z,1000)
println((b-a)*mean((x.^2).*pdf.(d,x)))
println((b-a)*mean(x.*pdf.(d,x)))
println((b-a)*mean(pdf.(d,x)))
#when D = 1000000, the simulated integral approximates the true value better than that when D = 1000

function mixed_logit_normal(theta, X, Z, y)

    alpha = theta[1:end-2]
    mu_gamma = theta[end-1]
    sigma_gamma = theta[end]
    K = size(X,2)
    J = length(unique(y))
    N = length(y)
    bigY = zeros(N,J)
    for j=1:J
        bigY[:,j] = y.==j
    end
    bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

    T = promote_type(eltype(X),eltype(theta))
    num   = zeros(T,N,J)
    dem   = zeros(T,N)
    d = Normal(mu_gamma,sigma_gamma)
    gammas, weights = lgwt(10,-4,4)
    M = size(gammas,1)
    F = ones(N,M)
    for i =1:M
        for j=1:J
            num[:,j] = exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J]).*gammas[i])
            dem .+= num[:,j]
        end
        P = (num./repeat(dem,1,J))

        for j=1:J
            P[:,j]=P[:,j].^bigY[:,j]
            F[:,i].*=P[:,j]
        end
    end

    loglike = -sum([log(sum(weights.*F[i,:].*pdf.(d,gammas))) for i = 1:N])

    return loglike
end

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 4
#:::::::::::::::::::::::::::::::::::::::::::::::::::

startvals = [2*rand(7*size(X,2)).-1; 0; 1]
obj = optimize(alpha->mixed_logit_normal(alpha, X, Z, y), startvals, LBFGS(), Optim.Options(g_tol=1e-5, iterations=100_000, show_trace=true, show_every=50))
beta = obj.minimizer[1:end-2]
mu_gamma = obj.minimizer[end-1]
sigma_gamma = obj.minimizer[end]
println('beta is ',beta)
println('mu_gamma is ',mu_gamma)
println('sigma_gamma is ', sigma_gamma)

#I also tried this code yet MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{var"#7#8", Float64}, Float64, 12})
# startvals = [2*rand(7*size(X,2)).-1; 0; 1]
# td = TwiceDifferentiable(theta -> mixed_logit_normal(theta, X, Z, y), startvals; autodiff = :forward)
# theta_hat_optim_ad = optimize(td, startvals, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100000, show_trace=true, show_every=50))
# theta_hat_mle_ad = theta_hat_optim_ad.minimizer
# H  = Optim.hessian!(td, theta_hat_mle_ad)
# theta_hat_mle_ad_se = sqrt.(diag(inv(H)))
# println([theta_hat_mle_ad theta_hat_mle_ad_se]) # these standard errors match Stata

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 5
#:::::::::::::::::::::::::::::::::::::::::::::::::::

function mixed_logit_normal2(theta, X, Z, y)

    alpha = theta[1:end-2]
    mu_gamma = theta[end-1]
    sigma_gamma = theta[end]
    K = size(X,2)
    J = length(unique(y))
    N = length(y)
    bigY = zeros(N,J)
    for j=1:J
        bigY[:,j] = y.==j
    end
    bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

    T = promote_type(eltype(X),eltype(theta))
    num   = zeros(T,N,J)
    dem   = zeros(T,N)
    d = Normal(mu_gamma,sigma_gamma)
    gammas, weights = lgwt(10,-4,4)
    M = size(gammas,1)
    F = ones(N,M)
    for i =1:M
        for j=1:J
            num[:,j] = exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J]).*gammas[i])
            dem .+= num[:,j]
        end
        P = (num./repeat(dem,1,J))

        for j=1:J
            P[:,j]=P[:,j].^bigY[:,j]
            F[:,i].*=P[:,j]
        end
    end

    a = -5*sigma_gamma
    b = 5*sigma_gamma
    z = Uniform(a,b)
    x = rand(z,1000000)
    loglike = -sum([log((b-a)*mean(F[i,:].*pdf.(d,gammas))) for i = 1:N])

    return loglike
end

startvals = [2*rand(7*size(X,2)).-1; 0; 1]
obj = optimize(alpha->mixed_logit_normal2(alpha, X, Z, y), startvals, LBFGS(), Optim.Options(g_tol=1e-5, iterations=100_000, show_trace=true, show_every=50))
beta = obj.minimizer[1:end-2]
mu_gamma = obj.minimizer[end-1]
sigma_gamma = obj.minimizer[end]
println('beta is ',beta)
println('mu_gamma is ',mu_gamma)
println('sigma_gamma is ', sigma_gamma)


text:
using CSV, JLD2, FileIO, DataFrames, Dates

println("Iniciando a redução das matrizes de 3 colunas para um CSV único.")

root_dir = dirname(pwd())
unzipped_files_path = *(root_dir, "/Datasets/Dataset_Bovespa/Unzipped/")
matrix_jld2_file_name = (unzipped_files_path * "COTAHIST" * ".jld2")

# caminho para matriz reduzida
reduced_matrix_file_name = (unzipped_files_path * "COTAHIST_REDUCED" * ".jld2")

csv = open((unzipped_files_path*"COTAHIST_PIVOTED"*".csv"),"w")

file_begin_counter = 1995
file_end_counter = 2019
file_counter = file_begin_counter

dataframe_total = DataFrame(data = [], codneg = [], preult = [])

time_inicio = Dates.now()

while (file_counter <= file_end_counter)
    local file_counter_string = string(file_counter)

    # matriz com dados completos
    local matrix = load(reduced_matrix_file_name, "COTAHIST_REDUCED/" * file_counter_string)

    local dataframe_reduced = convert(DataFrame, matrix)

    # remapea os dados
    dataframe_reduced[!,1] = map(x -> convert(Date,x), dataframe_reduced[!,1])
    dataframe_reduced[!,2] = map(x -> convert(String,x), dataframe_reduced[!,2])
    dataframe_reduced[!,3] = map(x -> convert(Float32,x), dataframe_reduced[!,3])

    # renomea as colunas para bater com o outro dataframe
    names!(dataframe_reduced, [:data, :codneg, :preult])
    append!(dataframe_total, dataframe_reduced)

    global file_counter += 1
    println(file_counter_string)
end

# pivota a tabela
dataframe_pivotado = unstack(dataframe_total, :data, :codneg, :preult) # pivota a tabela para as acoes
dataframe_pivotado = coalesce.(dataframe_pivotado, 0)  # substitui missing por 0
#println(names(dataframe_pivotado))
CSV.write(csv, dataframe_pivotado)

time_fim = Dates.now()
println("TI: " * Dates.format(time_inicio, "S:s"))
println("TF: " * Dates.format(time_fim, "S:s"))
println(time_fim - time_inicio)

println("Finalizando a redução das matrizes de 3 colunas para um CSV único.")


text:
@nospecialize
"""the view point of the type tag for each term
"""
function tag_extract(points_of_view::Dict{Function,Int})
    viewpoint = points_of_view[tag_extract]

    function and(many)
        @assert !isempty(many)
        ts = getindex.(many, viewpoint)
        t = reduce(typeintersect, ts)
        if t === Base.Bottom
            core_msg = "and patterns require an intersection of $(ts), which seems empty!"
            error(core_msg)
        end
        t
    end

    function or(many)
        ts = getindex.(many, viewpoint)
        Union{ts...}
    end

    function literal(val)
        typeof(val)
    end

    wildcard = Any

    function decons(comp::PComp, _, ns)
        targs = getindex.(ns, viewpoint)
        try
            comp.tcons(targs...)
        catch e
            join(map(repr, targs), ",")
            if e isa MethodError && e.f === comp.tcons
                argstr = join(repeat(String["_"], length(targs)), ", ")
                error("invalid deconstructor $(comp.repr)($(argstr))")
            end
            rethrow()
        end
    end

    guard(_) = Any
    effect(_) = Any
    (
        and = and,
        or = or,
        literal = literal,
        wildcard = wildcard,
        decons = decons,
        guard = guard,
        effect = effect
    )
end
@specialize


text:
function doji{T,N}(ohlc::TimeArray{T,N}; width=.01, op="Open", hi="High", lo="Low", cl="Close")
  res = abs((ohlc[op] .- ohlc[cl])) ./ (ohlc[hi] .- ohlc[lo]) .< width
  TimeArray(res.timestamp, res.values, ["doji"], ohlc.meta)
end

function hammer(x)
  #code here
end

function inverted_hammer(x)
  #code here
end

function hanging_man(x)
  #code here
end

function mirabosu(x)
  #code here
end

function shooting_star(x)
  #code here
end

function spinning_top(x)
  #code here
end

function three_white_soldiers(x)
  #code here
end

function three_black_crows(x)
  #code here
end

function morning_star(x)
  #code here
end


text:
println("##### methods #####")

using Test
using TreeTools

## Testing equality operator
root_1 = TreeTools.read_newick("$(dirname(pathof(TreeTools)))/../test/methods/tree1.nwk")
root_2 = TreeTools.read_newick("$(dirname(pathof(TreeTools)))/../test/methods/tree1_reordered.nwk")
@testset "Equality `==`" begin
	@test root_1 == root_2
end

@testset "node2tree" begin
	@test typeof(node2tree(root_1)) <: Tree
	@test typeof(node2tree(root_2)) <: Tree
end

# Testing ancestors
@testset "Ancestors" begin
    root_1 = TreeTools.read_newick("$(dirname(pathof(TreeTools)))/../test/methods/tree1.nwk")
    @test TreeTools.isancestor(root_1, root_1.child[1])
    @test TreeTools.isancestor(root_1, root_1.child[1].child[1])
    @test !TreeTools.isancestor(root_1.child[1],root_1.child[2])
    root_2 = TreeTools.read_newick("$(dirname(pathof(TreeTools)))/../test/methods/tree2.nwk")
    @test lca((root_2.child[1].child[1], root_2.child[1].child[2])).label == "ABC"
    @test lca((root_2.child[1].child[1].child[1], root_2.child[1].child[1].child[2], root_2.child[1].child[2])).label == "ABC"
    @test lca((root_2.child[1].child[1], root_2.child[2])).label == "ABCD"
end

@testset "Count" begin
	t1 = node2tree(root_1)
	@test count(isleaf, t1) == 4
	@test count(n -> n.label[1] == 'A', t1) == 3
	@test count(isleaf, t1.lnodes["AB"]) == 2
	@test count(n -> n.label[1] == 'A', t1.lnodes["AB"]) == 2
end

@testset "Copy type" begin
	t1 = Tree(TreeNode(MiscData(Dict(1=>2))))
	tc = copy(t1)
	@test tc.root.data[1] == 2
	t1.root.data[1] = 3
	@test tc.root.data[1] == 2
end

@testset "Copy" begin
	t1 = node2tree(root_1)
	t2 = copy(t1)
	@test typeof(t1) == typeof(t2)
	prunesubtree!(t2, ["A"])
	@test haskey(t1.lnodes, "A")
	@test !haskey(t2.lnodes, "A")
end

@testset "Convert" begin
	t1 = Tree(TreeNode(MiscData(Dict(1=>2))))
	# No op
	@test convert(Tree{MiscData}, t1) === t1
	@test convert(Tree{MiscData}, t1).root.data === t1.root.data

	# Converting to EmptyData and back
	t2 = convert(Tree{TreeTools.EmptyData}, t1)
	t3 = convert(Tree{MiscData}, t2)
	@test typeof(t2) == Tree{TreeTools.EmptyData}
	@test t2.root.data == TreeTools.EmptyData()

	@test typeof(t3) == Tree{TreeTools.MiscData}
	@test !haskey(t3.root.data, 1)

	###

	t1 = Tree(TreeNode(TreeTools.EmptyData()))
	# No op
	@test convert(Tree{TreeTools.EmptyData}, t1) === t1

	# Converting to MiscData and back
	t2 = convert(Tree{TreeTools.MiscData}, t1)
	@test typeof(t2) == Tree{TreeTools.MiscData}
	@test isempty(t2.root.data)

	@test typeof(convert(Tree{TreeTools.EmptyData}, t2)) == Tree{TreeTools.EmptyData}
end

nwk = "(A:3,(B:1,C:1):2)"
@testset "Distance" begin
	t = parse_newick_string(nwk)
	# Branch length
	@test distance(t, "A", "B") == 6
	@test distance(t.lnodes["A"], t.lnodes["B"]) == 6
	@test distance(t, "A", "B") == distance(t, "A", "C")
	@test distance(t.root, t.lnodes["A"]) == 3
	# Topological
	@test distance(t.root, t.lnodes["A"]; topological=true) == 1
	@test distance(t, "A", "B"; topological=true) == distance(t, "A", "C"; topological=true)
	@test distance(t, "A", "B"; topological=true) == 3
	for n in nodes(t)
		@test distance(t[n.label], t[n.label]; topological=true) == 0
		@test distance(t[n.label], t[n.label]; topological=false) == 0
	end
	# tests below can be removed when `divtime` is removed
	@test divtime(t.lnodes["A"], t.lnodes["B"]) == 6
	@test divtime(t.root, t.lnodes["A"]) == 3
end


text:
struct ProjectionOptSet{T, RepT<:Rep{T}, I} <: MOI.AbstractVectorSet
    p::Projection{RepT, I}
end
MOI.dimension(set::ProjectionOptSet) = fulldim(set.p)
Base.copy(set::ProjectionOptSet) = set

struct ProjectionBridge{T, F, RepT, I} <: MOI.Bridges.Constraint.AbstractBridge
    variables::Vector{MOI.VariableIndex}
    constraint::MOI.ConstraintIndex{F, PolyhedraOptSet{T, RepT}}
    dimensions::I
end
function MOI.Bridges.Constraint.bridge_constraint(
    ::Type{ProjectionBridge{T, F, RepT, I}}, model::MOI.ModelLike,
    f::MOI.AbstractVectorFunction, p::ProjectionOptSet) where {T, F, RepT, I}
    vf = MOIU.eachscalar(f)
    func = Vector{eltype(vf)}(undef, fulldim(p.p.set))
    for (i, j) in enumerate(p.p.dimensions)
        func[j] = vf[i]
    end
    N = fulldim(p.p.set)
    variables = MOI.add_variables(model, N - fulldim(p.p))
    for (i, j) in enumerate(setdiff(1:N, p.p.dimensions))
        func[j] = variables[i]
    end
    constraint = MOI.add_constraint(model, MOI.Utilities.vectorize(func), PolyhedraOptSet(p.p.set))
    return ProjectionBridge{T, F, RepT, I}(variables, constraint, p.p.dimensions)
end

MOI.supports_constraint(::Type{ProjectionBridge{T}},
                        ::Type{<:MOI.AbstractVectorFunction},
                        ::Type{<:ProjectionOptSet{T}}) where {T} = true
function MOI.Bridges.added_constrained_variable_types(::Type{<:ProjectionBridge})
    return Tuple{DataType}[]
end
function MOI.Bridges.added_constraint_types(::Type{ProjectionBridge{T, F, RepT, I}}) where {T, F, RepT, I}
    return [(F, PolyhedraOptSet{T, RepT})]
end
function MOI.Bridges.Constraint.concrete_bridge_type(
    ::Type{<:ProjectionBridge{T}},
    F::Type{<:MOI.AbstractVectorFunction},
    ::Type{ProjectionOptSet{T, RepT, I}}) where {T, RepT, I}
    return ProjectionBridge{T, F, RepT, I}
end

# Attributes, Bridge acting as an model
function MOI.get(b::ProjectionBridge{T, F, RepT, I}, ::MOI.NumberOfConstraints{F, ProjectionOptSet{RepT, I}}) where {T, F, RepT, I}
    return 1
end
function MOI.get(b::ProjectionBridge{T, F, RepT, I}, ::MOI.ListOfConstraintIndices{F, ProjectionOptSet{RepT, I}}) where {T, F, RepT, I}
    return [b.constraint]
end

# Indices
function MOI.delete(model::MOI.ModelLike, b::ProjectionBridge)
    MOI.delete(model, b.constraint)
    for vi in b.variables
        MOI.delete(model, vi)
    end
end
_moi_set(set::Projection) = ProjectionOptSet(set)
function JuMP.build_constraint(error_fun::Function, func, set::Projection)
    return JuMP.BridgeableConstraint(JuMP.BridgeableConstraint(
        JuMP.build_constraint(error_fun, func, _moi_set(set)),
        ProjectionBridge), PolyhedraToLPBridge)
end

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
import Oceananigans: short_show

location_str(::Type{Face})    = "Face"
location_str(::Type{Center})  = "Center"
location_str(::Type{Nothing}) = "⋅"

show_location(X, Y, Z) = "($(location_str(X)), $(location_str(Y)), $(location_str(Z)))"

show_location(field::AbstractField{X, Y, Z}) where {X, Y, Z} = show_location(X, Y, Z)

short_show(m::Missing) = "$m"

short_show(field::AbstractField) = string(typeof(field).name.wrapper, " located at ", show_location(field))
short_show(field::AveragedField) = string("AveragedField over dims=$(field.dims) located at ", show_location(field), " of ", short_show(field.operand))
short_show(field::ComputedField) = string("ComputedField located at ", show_location(field), " of ", short_show(field.operand))

Base.show(io::IO, field::AbstractField{X, Y, Z, A, G}) where {X, Y, Z, A, G} =
    print(io, "$(short_show(field))\n",
          "├── architecture: $A",
          "└── grid: $G")

Base.show(io::IO, field::Field) =
    print(io, "$(short_show(field))\n",
          "├── data: $(typeof(field.data)), size: $(size(field.data))\n",
          "├── grid: $(short_show(field.grid))\n",
          "└── boundary conditions: $(short_show(field.boundary_conditions))")

show_status(::Nothing) = "nothing"
show_status(status) = "time=$(status.time)"

Base.show(io::IO, field::AveragedField) =
    print(io, "$(short_show(field))\n",
          "├── data: $(typeof(field.data)), size: $(size(field.data))\n",
          "├── grid: $(short_show(field.grid))", '\n',
          "├── dims: $(field.dims)", '\n',
          "├── operand: $(short_show(field.operand))", '\n',
          "└── status: ", show_status(field.status), '\n')

Base.show(io::IO, field::ComputedField) =
    print(io, "$(short_show(field))\n",
          "├── data: $(typeof(field.data)), size: $(size(field.data))\n",
          "├── grid: $(short_show(field.grid))", '\n',
          "├── operand: $(short_show(field.operand))", '\n',
          "└── status: ", show_status(field.status), '\n')

Base.show(io::IO, field::KernelComputedField) =
    print(io, "$(short_show(field))\n",
          "├── data: $(typeof(field.data)), size: $(size(field.data))\n",
          "├── grid: $(short_show(field.grid))", '\n',
          "├── computed_dependencies: $(Tuple(short_show(d) for d in field.computed_dependencies))", '\n',
          "├── kernel: $(short_show(field.kernel))", '\n',
          "└── status: ", show_status(field.status), '\n')

short_show(array::OffsetArray{T, D, A}) where {T, D, A} = string("OffsetArray{$T, $D, $A}")

Base.show(io::IO, ::MIME"text/plain", f::AbstractField) = show(io, f)


text:
@adddevice(StorageUnit, OnePort, LinearExpansionForm{LinearDispatchForm},
           :storageunits, (:S, :T=>:snapshots), joinpath(@__DIR__, "attrs", "storageunits.csv"))

## StorageUnit
cost(d::StorageUnit) = sum(d[:marginal_cost] .* d[:p_dispatch]) + sum(d[:capital_cost] .* (d[:p_nom] - getparam(d, :p_nom)))

function addtoexpr!(injection::AxisArray, ::NodalActivePower,
                    m::EnergyModel, d::StorageUnit)
    B, T = AxisArrays.axes(injection)
    D = axis(m, d)

    buses = get(d, :bus, D)
    p_dispatch = get(d, :p_dispatch, D, T)
    p_store = get(d, :p_store, D, T)

    for (i, idx) in enumerate(indexin(buses, B)), t in T
        if !isnothing(idx)
            add_to_expression!(injection[idx,t], p_dispatch[i,t])
            add_to_expression!(injection[idx,t], -1, p_store[i,t])
        end
    end
end

function addto!(jm::ModelView, m::EnergyModel, d::StorageUnit{DF}) where
      {DDF, DF <: LinearExpansionForm{DDF}}

    addto!(jm, m, with_formulation(d, LinearExpansionInvestmentForm))
    addto!(jm, m, with_formulation(d, LinearExpansionDispatchForm{DDF}))
end

function addto!(jm::ModelView, m::EnergyModel, d::StorageUnit{DF}) where
      {DF <: LinearExpansionInvestmentForm}

    S = axis(m, d)

    p_nom_min = get(d, :p_nom_min, S)
    p_nom_max = get(d, :p_nom_max, S)

    @variable(jm, p_nom_min[s] <= p_nom[s=S] <= p_nom_max[s])
end

function addto!(jm::ModelView, m::EnergyModel, d::StorageUnit{DF}) where
      {DDF, DF <: LinearExpansionDispatchForm{DDF}}

    S = axis(m, d)
    T = axis(m, :snapshots)

    p_nom = get(d, :p_nom, S)
    p_min_pu = get(d, :p_min_pu, S, T)
    p_max_pu = get(d, :p_max_pu, S, T)

    @variables jm begin
        p_store[s=S,t=T] >= 0
        p_dispatch[s=S,t=T] >= 0
    end

    @constraints jm begin
        p_lower[s=S,t=T], p_store[s,t] <= - p_min_pu[s,t] * p_nom[s]
        p_upper[s=S,t=T], p_dispatch[s,t] <= p_max_pu[s,t] * p_nom[s]
    end

    addto!_state_of_charge(jm, m, d)
end

function addto!(jm::ModelView, m::EnergyModel, d::StorageUnit{DF}) where
      {DF <: DispatchForm}

    S = axis(m, d)
    T = axis(m, :snapshots)

    p_nom = get(d, :p_nom, S)
    p_min_pu = get(d, :p_min_pu, S, T)
    p_max_pu = get(d, :p_max_pu, S, T)

    @variables jm begin
        0 <= p_store[s=S,t=T] <= - p_min_pu[s,t] * p_nom[s]
        0 <= p_dispatch[s=S,t=T] <= p_max_pu[s,t] * p_nom[s]
    end

    addto!_state_of_charge(jm, m, d)
end

function addto!_state_of_charge(jm::ModelView, m::EnergyModel, d::StorageUnit)
    S = axis(m, d)
    T = axis(m, :snapshots)

    p_nom = get(d, :p_nom, S)
    p_store = get(d, :p_store)
    p_dispatch = get(d, :p_dispatch)

    inflow = get(d, :inflow, S, T)
    max_hours = get(d, :max_hours, S, T)
    efficiency_store = get(d, :efficiency_store, S, T)
    efficiency_dispatch = get(d, :efficiency_dispatch, S, T)

    @variables jm begin
        0 <= spill[s=S,t=T] <= inflow[s,t]
        0 <= state_of_charge[s=S,t=T]
    end

    @constraint(jm, state_of_charge_upper[s=S,t=T], state_of_charge[s,t] <= max_hours[s,t] * p_nom[s])

    if d[:cyclic_state_of_charge]
        soc_prev = circshift(state_of_charge, :snapshots=>1)
    else
        soc_prev = similar(state_of_charge, Union{Float64,eltype(soc)})
        soc_prev[:,1] .= get(d, :state_of_charge_initial, S)
        soc_prev[:,2:end] .= state_of_charge[:,1:end-1]
    end

    @constraint(jm, state_of_charge_eq[s=S, t=T],
                state_of_charge[s,t] - soc_prev[s,t]
                == p_store[s,t] * efficiency_store[s,t]
                - p_dispatch[s,t] / efficiency_dispatch[s,t]
                + inflow[s,t] - spill[s,t])
end

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
export hsvtorgb

"""
   hsvtorgb(color)

Convert a `color` from HSV space to RGB.
"""
hsvtorgb(color) = begin
    H, S, V = color
    C = V * S
    X = C * (1 - Base.abs((H / 60) % 2 - 1))
    m = V - C
    if 0 ≤ H < 60
        R′, G′, B′ = C, X, 0
    elseif 60 ≤ H < 120
        R′, G′, B′ = X, C, 0
    elseif 120 ≤ H < 180
        R′, G′, B′ = 0, C, X
    elseif 180 ≤ H < 240
        R′, G′, B′ = 0, X, C
    elseif 240 ≤ H < 300
        R′, G′, B′ = X, 0, C
    elseif 300 ≤ H < 360
        R′, G′, B′ = C, 0, X
    else
        R′, G′, B′ = rand(3)
    end
    R, G, B = R′ + m, G′ + m, B′ + m
    [R; G; B]
end

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
# # Providing uncertainty with Measurements.jl
#
# [Measurements.jl](https://github.com/JuliaPhysics/Measurements.jl) provides automatic
# error propagation, and integrates seamlessly with much of the Julia ecosystem.
#
# Here, we'll slightly modify the [Daisyworld](@ref) example, to simulate some measurement
# uncertainty in our world's parameters.
#
# ## Setup
#
# First we'll construct our agents.

using Agents
using Measurements

mutable struct Daisy <: AbstractAgent
    id::Int
    pos::Dims{2}
    breed::Symbol
    age::Int
    albedo::AbstractFloat # Allow Measurements
end

mutable struct Land <: AbstractAgent
    id::Int
    pos::Dims{2}
    temperature::AbstractFloat # Allow Measurements
end

# Notice that there is only one small difference between this version and the original
# example model: the use of `AbstractFloat` instead of `Float64` for the `albedo` and
# `temperature` parameters. Behaviour between these two types is practically equivalent
# from our perspective, but it allows us to use an uncertain value for our two parameters.
# `1.0 ± 0.1` rather than `1.0` for example. We could also be specific here and bind the
# parameters with type `Measurement{Float64}` as well.
#
# Next, we'll implement all the important functions for DaisyWorld. If you want to know what
# each of these functions do, see the [Daisyworld](@ref) example, as they are copied directly
# from there.

using Plots
using Statistics: mean
import DrWatson: @dict
import StatsBase
using Random # hide

const DaisyWorld = ABM{<:GridSpace,Union{Daisy,Land}}

function update_surface_temperature!(pos::Dims{2}, model::DaisyWorld)
    ids = ids_in_position(pos, model)
    absorbed_luminosity = if length(ids) == 1
        (1 - model.surface_albedo) * model.solar_luminosity
    else
        (1 - model[ids[2]].albedo) * model.solar_luminosity
    end
    local_heating = absorbed_luminosity > 0 ? 72 * log(absorbed_luminosity) + 80 : 80
    T0 = model[ids[1]].temperature
    model[ids[1]].temperature = (T0 + local_heating) / 2
end

function diffuse_temperature!(pos::Dims{2}, model::DaisyWorld)
    ratio = get(model.properties, :ratio, 0.5)
    ids = nearby_ids(pos, model)
    meantemp = sum(model[i].temperature for i in ids if model[i] isa Land) / 8
    land = model[ids_in_position(pos, model)[1]]
    land.temperature = (1 - ratio) * land.temperature + ratio * meantemp
end

function propagate!(pos::Dims{2}, model::DaisyWorld)
    ids = ids_in_position(pos, model)
    if length(ids) > 1
        daisy = model[ids[2]]
        temperature = model[ids[1]].temperature
        seed_threshold = (0.1457 * temperature - 0.0032 * temperature^2) - 0.6443
        if rand() < seed_threshold
            empty_neighbors = Tuple{Int,Int}[]
            neighbors = nearby_positions(pos, model)
            for n in neighbors
                if length(ids_in_position(n, model)) == 1
                    push!(empty_neighbors, n)
                end
            end
            if !isempty(empty_neighbors)
                seeding_place = rand(empty_neighbors)
                a = Daisy(nextid(model), seeding_place, daisy.breed, 0, daisy.albedo)
                add_agent_pos!(a, model)
            end
        end
    end
end

function agent_step!(agent::Daisy, model::DaisyWorld)
    agent.age += 1
    agent.age >= model.max_age && kill_agent!(agent, model)
end

agent_step!(agent::Land, model::DaisyWorld) = nothing

function model_step!(model)
    for p in positions(model)
        update_surface_temperature!(p, model)
        diffuse_temperature!(p, model)
        propagate!(p, model)
    end
    model.tick += 1
    solar_activity!(model)
end

function solar_activity!(model::DaisyWorld)
    if model.scenario == :ramp
        if model.tick > 200 && model.tick <= 400
            model.solar_luminosity += model.solar_change
        end
        if model.tick > 500 && model.tick <= 750
            model.solar_luminosity -= model.solar_change / 2
        end
    elseif model.scenario == :change
        model.solar_luminosity += model.solar_change
    end
end

# ## Adding Uncertainty
#
# Now, we can write a constructor function, and use uncertainly values which will propagate
# automatically through our model.

function daisyworld(;
    griddims = (30, 30),
    max_age = 25,
    init_white = 0.2,
    init_black = 0.2,
    albedo_white = 0.75,
    albedo_black = 0.25,
    ## Surface albedo measurements are complicated for our satellites perhaps
    surface_albedo = 0.4 ± 0.15,
    ## Measurements from the sun are generally stable, but fluctuate around 10%
    solar_change = 0.005 ± 0.002,
    solar_luminosity = 1.0 ± 0.1,
    scenario = :default,
)

    space = GridSpace(griddims)
    properties = @dict max_age surface_albedo solar_luminosity solar_change scenario
    properties[:tick] = 0
    daisysched(model) = [a.id for a in allagents(model) if a isa Daisy]
    model = ABM(
        Union{Daisy,Land},
        space;
        scheduler = daisysched,
        properties = properties,
        warn = false,
    )

    ## An uncertain initial temperature, solely for type stability
    fill_space!(Land, model, 0.0 ± 0.0)
    grid = collect(positions(model))
    num_positions = prod(griddims)
    white_positions =
        StatsBase.sample(grid, Int(init_white * num_positions); replace = false)
    for wp in white_positions
        wd = Daisy(nextid(model), wp, :white, rand(0:max_age), albedo_white)
        add_agent_pos!(wd, model)
    end
    allowed = setdiff(grid, white_positions)
    black_positions =
        StatsBase.sample(allowed, Int(init_black * num_positions); replace = false)
    for bp in black_positions
        wd = Daisy(nextid(model), bp, :black, rand(0:max_age), albedo_black)
        add_agent_pos!(wd, model)
    end

    return model
end

# You see we've included uncertainty in four places: surface albedo and initial temperature, and
# the two solar luminosity values. We do not require changes to any model code, nor handle
# these parameters in any special way; for example `2.0 * surface_albedo` is a regular operation.
# Errors will be propagated under the hood automatically.

# ## Visualizing the Result
#
# Similar to the [Daisyworld](@ref) example, we will now check out how the surface temperature
# and daisy count fares when solar luminosity ramps up.
#
# First, some helper functions

black(a) = a.breed == :black
white(a) = a.breed == :white
daisies(a) = a isa Daisy

land(a) = a isa Land
adata = [(black, count, daisies), (white, count, daisies), (:temperature, mean, land)]

mdata = [:solar_luminosity]

# And now the simulation

Random.seed!(19) # hide
model = daisyworld(scenario = :ramp)
agent_df, model_df =
    run!(model, agent_step!, model_step!, 1000; adata = adata, mdata = mdata)

p1 = plot(
    agent_df.step,
    agent_df.count_black_daisies,
    label = "black",
    ylabel = "Daisy count",
)
plot!(p1, agent_df.step, agent_df.count_white_daisies, label = "white")

p2 = plot(
    agent_df.step,
    Measurements.value.(agent_df[!, aggname(adata[3])]),
    ribbon = Measurements.uncertainty.(agent_df[!, aggname(adata[3])]),
    ylabel = "Temperature",
    legend = false,
)

p3 = plot(
    model_df.step,
    Measurements.value.(model_df.solar_luminosity),
    ribbon = Measurements.uncertainty.(model_df.solar_luminosity),
    ylabel = "Luminosity",
    xlabel = "Years",
    legend = false,
)

plot(p1, p2, p3, layout = (3, 1), size = (600, 700))

# [Plots.jl](https://github.com/JuliaPlots/Plots.jl) automatically adds error bars onto series that
# have `Measurement{Float64}` types. However, we use `ribbon`s above to display the uncertainty since
# we have so many points. Notice that the error has automatically propagated throughout the model.


text:
@slackmethod(TestMethodMacro, "method.macro",
    begin
        a::Int64
    end,

    begin
        b::UTF8String
        c::Int64
    end)

facts("Defining requests and responses with a macro") do
    req = TestMethodMacro(Token(utf8("foo")), Int64(42))
    @fact req.token --> utf8("foo")
    @fact req.a --> 42

    resp = TestMethodMacroResponse(utf8("bar"), 17)
    @fact resp.b --> utf8("bar")
    @fact resp.c --> 17

    @fact getresponsetype(TestMethodMacro) --> TestMethodMacroResponse
    @fact method_name(TestMethodMacro) --> "method.macro"
end