text:
mutable struct PathComboBox <: GtkComboBoxText
    handle::Ptr{Gtk.GObject}
    entry::GtkEntry
    time_last_keypress::AbstractFloat
    main_window::MainWindow

    function PathComboBox(main_window::MainWindow)
        cbox = GtkComboBoxText(true)
        entry = cbox[1]

        p = new(cbox.handle, entry, time(), main_window)
        Gtk.gobject_move_ref(p, cbox)
    end
end

update_pathEntry(pathCBox::PathComboBox, path) = set_gtk_property!(pathCBox.entry, :text, path)
update_pathEntry(main_window::MainWindow, path) = update_pathEntry(main_window.pathCBox, path)

function pathEntry_key_press_cb(widgetptr::Ptr, eventptr::Ptr, user_data)
    widget = convert(GtkEntry, widgetptr)
    event = convert(Gtk.GdkEvent, eventptr)

    pathCBox = user_data
    pathCBox.time_last_keypress = time()

    if event.keyval == GdkKeySyms.Return
        path = get_gtk_property(widget, :text, AbstractString)
        try
            remotecall_fetch(cd, worker(current_console(pathCBox.main_window)), path)
        catch err
            println(string(err))
        end
        on_path_change(pathCBox.main_window, true)
    end

    return convert(Cint, false)
end

# I think I don't need this anymore but not 100% sure
@guarded (nothing) function pathDbox_changed_cb(ptr::Ptr, user_data)

    pathCBox = convert(GtkComboBoxText, ptr)

    #@show time() - pathCBox.time_last_keypress

    if time() - pathCBox.time_last_keypress < 0.1
        return nothing
    end

    pth = unsafe_string(Gtk.G_.active_text(pathCBox))
    try
        #cd(pth)
    catch err
        println(string(err))
    end

    #on_path_change(pathCBox.main_window)
    nothing
end

function init!(pathCBox::PathComboBox)

    signal_connect(pathEntry_key_press_cb, pathCBox.entry, "key-press-event",
    Cint, (Ptr{Gtk.GdkEvent}, ), false, pathCBox)

    signal_connect(pathDbox_changed_cb, pathCBox, "changed", Nothing, (), false)

    set_gtk_property!(pathCBox.entry, :width_request, 400)
    update_pathEntry(pathCBox, pwd())
    set_gtk_property!(pathCBox.entry, :hexpand, true)

    sc = Gtk.G_.style_context(pathCBox.entry)
    push!(sc, style_provider(pathCBox.main_window), 600)
end


text:
# SparseArrays.dropstored!
size(env::NatureEnv) = env.world_size[1:2]

function NatureEnv12(;
        num_starting_players=2,
        world_size=(32, 32),
        window=3,
        episode_len=100,
        num_frames=4,
        player_starting_food=5f0,
        food_generators=[
            FoodGen([5,5],[1,1]),
            FoodGen([25,25],[1,1]),
           ],
        vocab_size=0) #width, heigh, channels, batch

    num_food_types = length(food_generators)
    # self, other, x pos, y pos, food_poses, food_counts
    num_channels = num_frames * (2+2+num_food_types+num_food_types+vocab_size)

    observation_space = Space(
        ClosedInterval.(
            fill(typemin(Int), 2window+1, 2window+1, num_channels),
            fill(typemax(Int), 2window+1, 2window+1, num_channels)
        )
    )
    place_records = bigger_dd_builder()

    NatureEnv12{num_food_types}(
        1,
        episode_len,
        num_starting_players,
        [],
        num_food_types,
        Vector{Matrix{Float32}}(),
        food_generators,
        player_starting_food,
        world_size,
        window,
        (2*window+1, 2*window+1, num_channels),
        observation_space,
        place_records,
        zeros(Float32, num_food_types),
        vocab_size,
        [[] for _ in 1:episode_len],
        num_frames,
        [[] for _ in 1:num_starting_players]
       )
end

RLBase.state_space(env::NatureEnv,::Observation{Any}, players::Dict{Int,Player}) = Dict(p=>state_space(env, p) for p in keys(players))
RLBase.state(env::NatureEnv,::Observation{Any}, players::Vector{Int}) = Dict(p=>state(env, p) for p in players)
RLBase.action_space(env::NatureEnv, players::Dict{Int,Player}) = Dict(p=>action_space(env, p) for p in keys(players))
RLBase.reward(env::NatureEnv, players::Dict{Int,Player}) = Dict(p=>reward(env, p) for p in keys(players))
RLBase.reward(env::NatureEnv, players::Vector{Int}) = Dict(p=>reward(env, p) for p in players)

RLBase.state_space(env::NatureEnv,::Observation{Any}, player::Int) = env.observation_space
RLBase.action_space(env::NatureEnv, player::Int) = Base.OneTo(4 + 2env.food_types + env.vocab_size)

function RLBase.reward(env::NatureEnv, player::Int)
    if all(env.players[player].food_counts .> 0f0) && !env.players[player].dead
        return 1f0
    else
        return 0f0
    end
end

RLBase.is_terminated(env::NatureEnv, player::Int) = env.players[player].dead
RLBase.is_terminated(env::NatureEnv, players::Vector{Int}) = Dict(p=>is_terminated(env, p) for p in players)
function RLBase.is_terminated(env::NatureEnv)
    if all(is_terminated(env, p) for p in keys(env.players)) ||
        env.step >= env.episode_len
        return true
    else
        return false
    end
end
function RLBase.reset!(env::NatureEnv)
    # Add Food according to each generator
    env.step=1
    env.food_frames = []
    for type in 1:env.food_types
        push!(env.food_frames, zeros(Float32, size(env)...))
        for nf::Tuple in new_food(env.food_generators[type], 10000)
            env.food_frames[type][nf...] += 0.1f0
        end
    end

    env.players = [Player(rand_pos(env)...,env.food_types, env.player_starting_food) for _ in 1:env.num_starting_players]
    env.comms = [[] for _ in 1:env.episode_len]
    frame_channels = Int(env.obs_size[3] / env.num_frames)
    env.past_states = [
        [[zeros(Float32, env.obs_size[1:2]..., frame_channels) for _ in 1:env.num_frames-1]...] for _ in 1:env.num_starting_players]

    empty!(env.place_record)
    env.exchanges = zeros(Float32, env.food_types)
    nothing
end

function RLBase.state(env::NatureEnv, player::Int)
    if env.players[player].dead
        return zeros(Float32, env.obs_size)
    end
    # Don't create a new frame if episode is over
    if env.step > env.episode_len
        return cat(env.past_states[player][end-env.num_frames+1:end]..., dims=3)
    end
    w = env.window
    px, py = env.players[player].pos
    self_frame = make_frame(size(env)..., w)
    other_frame = make_frame(size(env)..., w)
    if env.players[player].dead
        return zeros(Float32, env.obs_size...)
    end
    # Player frames
    for i in 1:length(env.players)
        if i == player
            self_frame[(w .+ env.players[i].pos)...] = 1
        else
            other_frame[(w .+ env.players[i].pos)...] = 1
        end
    end

    # positions
    xpos = hcat([i*ones(Float32, env.world_size[1]+(2*w)) for i in 1-w:env.world_size[1]+w]...) ./ (env.world_size[1]+w)
    ypos = vcat([i*ones(Float32, 1, env.world_size[2]+(2*w)) for i in 1-w:env.world_size[2]+w]...) ./ (env.world_size[2]+w)

    # Food frames and food counts
    food_frames = []
    food_counts = []
    for (f, env_food_frame) in enumerate(env.food_frames)
        ff = make_frame(size(env)..., w)
        fc = make_frame(size(env)..., w) .* 0f0
        ff[w+1:end-w,w+1:end-w] = env_food_frame ./ 10f0
        for p in env.players
            p.dead && continue
            fc[w+p.pos[1],w+p.pos[2]] = p.food_counts[f] / 10f0
        end
        push!(food_frames, ff)
        push!(food_counts, fc)
    end
    comm_frames = [make_frame(size(env)...,w) for _ in 1:env.vocab_size]
    for (x, y, _, symbol) in env.comms[env.step]
        comm_frames[symbol][w+x,w+y] = 1
    end

    frames = [self_frame, other_frame, xpos, ypos, food_frames..., food_counts..., comm_frames...]
    frames = cat(frames..., dims=3)
    frames = frames[px:px+(2*w), py:py+(2*w), :]
    if length(env.past_states[player]) == env.num_frames + env.step - 2
        push!(env.past_states[player], frames)
    elseif length(env.past_states[player]) == env.num_frames + env.step - 1
        env.past_states[player][end] = frames
    end
    frame_history = cat(env.past_states[player][end-env.num_frames+1:end]..., dims=3)
    frame_history
end

function (env::NatureEnv)(actions::Dict)
    if is_terminated(env)
        return nothing
    end
    env.step > env.episode_len && begin
        println("exceeded")
        return nothing
    end
    env.step += 1
    Dict(p=>env(a.action[1], p) for (p, a) in actions)
end

function (env::NatureEnv)(action::Int, player::Int)

    if is_terminated(env, player)
        return nothing
    end

    MOVE_RANGE = 1:4
    FOOD_RANGE = (MOVE_RANGE.stop+1)$(2*env.food_types)
    VOCAB_RANGE = (FOOD_RANGE.stop+1)$(env.vocab_size+1)

    # Players lose 0.1 food per tick, floored at 0
    if action in MOVE_RANGE
        env.players[player].food_counts = env.players[player].food_counts .- 0.15f0 .|> x->max(x, 0f0)
    else
        env.players[player].food_counts = env.players[player].food_counts .- 0.1f0 .|> x->max(x, 0f0)
    end

    (action in MOVE_RANGE) && move(env, player, action)
    (action in FOOD_RANGE) && food(env, player, action - FOOD_RANGE.start + 1)
    (action in VOCAB_RANGE) && comm(env, player, action - VOCAB_RANGE.start + 1)
    if min(env.players[player].food_counts...) <= 0f0
        env.players[player].dead = true
    end
    nothing
end

RLBase.legal_action_space_mask(::NatureEnv) = nothing
RLBase.current_player(env::NatureEnv) =  Dict{Int, Player}(i=>p for (i,p) in enumerate(env.players) if !p.dead)
# https://juliareinforcementlearning.org/docs/How_to_write_a_customized_environment/#More-Complicated-Environments
RLBase.NumAgentStyle(env::NatureEnv) = MultiAgent(length(env.players))
RLBase.DynamicStyle(::NatureEnv) = SIMULTANEOUS
RLBase.ActionStyle(::NatureEnv) = FULL_ACTION_SET
RLBase.InformationStyle(::NatureEnv) = IMPERFECT_INFORMATION
# RLBase.StateStyle(::TicTacToeEnv) =
#     (Observation{Any}(), Observation{Int}(), Observation{BitArray{3}}())
# RLBase.RewardStyle(::TicTacToeEnv) = TERMINAL_REWARD
# RLBase.UtilityStyle(::TicTacToeEnv) = GENERAL_SUM
RLBase.ChanceStyle(::NatureEnv) = DETERMINISTIC

ext: jl


text:
using Oceananigans.Solvers

using Oceananigans.ImmersedBoundaries: mask_immersed_velocities!, mask_immersed_field!

import Oceananigans.TimeSteppers: calculate_pressure_correction!, pressure_correct_velocities!

"""
    calculate_pressure_correction!(model::IncompressibleModel, Δt)

Calculate the (nonhydrostatic) pressure correction associated `tendencies`, `velocities`, and step size `Δt`.
"""
function calculate_pressure_correction!(model::IncompressibleModel, Δt)

    # Mask immersed velocities
    velocity_masking_events = mask_immersed_velocities!(model.velocities, model.architecture, model.grid)
    wait(device(model.architecture), MultiEvent(velocity_masking_events))

    fill_halo_regions!(model.velocities, model.architecture, model.clock, fields(model))

    solve_for_pressure!(model.pressures.pNHS, model.pressure_solver, model.architecture, model.grid, Δt, model.velocities)

    fill_halo_regions!(model.pressures.pNHS, model.architecture)

    return nothing
end

#####
##### Fractional and time stepping
#####

"""
Update the predictor velocities u, v, and w with the non-hydrostatic pressure via

    `u^{n+1} = u^n - δₓp_{NH} / Δx * Δt`
"""
@kernel function _pressure_correct_velocities!(U, grid, Δt, pNHS)
    i, j, k = @index(Global, NTuple)

    @inbounds U.u[i, j, k] -= ∂xᶠᶜᵃ(i, j, k, grid, pNHS) * Δt
    @inbounds U.v[i, j, k] -= ∂yᶜᶠᵃ(i, j, k, grid, pNHS) * Δt
    @inbounds U.w[i, j, k] -= ∂zᵃᵃᶠ(i, j, k, grid, pNHS) * Δt
end

"Update the solution variables (velocities and tracers)."
function pressure_correct_velocities!(model::IncompressibleModel, Δt)

    event = launch!(model.architecture, model.grid, :xyz,
                    _pressure_correct_velocities!,
                    model.velocities,
                    model.grid,
                    Δt,
                    model.pressures.pNHS,
                    dependencies=Event(device(model.architecture)))

    wait(device(model.architecture), event)

    return nothing
end


text:
"This module is created to build the SCR-NLP version of the MPPDC with quanity bids, with dual objective function, for all UL agents,
Selection of social planner's objective is based on the string SP_obj"
function build_EPEC_quantity_mod(t_steps, non_str_gens, non_str_res, str_agents, optimizer, SP_obj)
    ηch = 0.95 #Needs to be passed in
    ηdch = 0.95 #Needs to be passed in
    soc_init = dat_str.ESS_STR.E_max[1] * 0.5
    non_str_agents = vcat(generators, renewables)
    ####INIT ABSTRACT MODEL######
    println("...Initializing EPEC quantity abstract model...")
    function init_EPEC_quantity()
        ############INIT THE MODEL INSTANCE############
        EPEC =  Model()
        EPEC.ext[:param] = Dict()
        param = EPEC.ext[:param]
        EPEC.ext[:eq] = Dict()
        eq = EPEC.ext[:eq]
        #PRICE-TAKING AGENTS QUANTITY BIDS
        param[:G] = @NLparameter(EPEC, G[i in non_str_gens, t in t_steps] == 0) #Generation capacity of non-str generator
        param[:W] = @NLparameter(EPEC, W[k in non_str_res, t in t_steps]  == 0) #Generation capacity of non-str RES
        param[:D] = @NLparameter(EPEC, D[t in t_steps]  == 0) #Max load
        #########MAX POWER OF THE UL AGENTS###########
        param[:P_G_ST] = @NLparameter(EPEC, P_G_ST[t in t_steps] == 0)
        param[:P_W_ST] = @NLparameter(EPEC, P_W_ST[t in t_steps] == 0)
        param[:P_DCH_ST] = @NLparameter(EPEC, P_DCH_ST[t in t_steps] == 0)
        param[:P_CH_ST] = @NLparameter(EPEC, P_CH_ST[t in t_steps] == 0)
        param[:E_max] = @NLparameter(EPEC, E_max[t in t_steps] == 0) #max capacity of the ESS
        #########PRICE BIDS#############
        #NOTE:ESS HAS BOTH PRCE AND QUANTITY AS DEC VARS!
        param[:c_g_st] = @NLparameter(EPEC, c_g_st[t in t_steps] == 0) #Generation price bid str
        param[:c_w_st] = @NLparameter(EPEC, c_w_st[t in t_steps]  == 0)#Wind price bid str
        param[:c_g] = @NLparameter(EPEC, c_g[i in non_str_gens, t in t_steps] == 0) #Generation price bid nonstr
        param[:c_w] = @NLparameter(EPEC, c_w[k in non_str_res, t in t_steps] == 0) #Wind price bid nonstr
        param[:c_d] = @NLparameter(EPEC, c_d[t in t_steps] == 0) #WTP of the load
        #############UPDATING PARAMS#############
        for t in t_steps
            #MAXIMUM ENERGY CONTENT OF ESS
            JuMP.set_value(E_max[t], dat_str.ESS_STR.E_max[t]) #By this the upper bound var. on SoC is set
            #Updating max powers of UL agents
            JuMP.set_value(P_G_ST[t], dat_str.GEN_STR.P_max[t])
            JuMP.set_value(P_W_ST[t], dat_str.WIND_STR.P_max[t])
            JuMP.set_value(P_DCH_ST[t], dat_str.ESS_STR.P_max[t])
            JuMP.set_value(P_CH_ST[t], dat_str.ESS_STR.P_max[t])

            for i in non_str_gens
                JuMP.set_value(c_g[i,t], dat_nonstr[i].c_g[t]) #MC OF GEN
                JuMP.set_value(G[i,t], dat_nonstr[i].G[t]) #MAX CAPACITY OF GEN
            end
            for k in non_str_res
                JuMP.set_value(c_w[k,t], dat_nonstr[k].c_g[t]) #COST OF WIND
                JuMP.set_value(W[k,t], dat_nonstr[k].G[t]) #MAX WIND CAPACITY
            end
            JuMP.set_value(c_d[t], p_cap) #WTP SET TO PRICE CAP
            JuMP.set_value(D[t], load[t]) #MAX LOAD
            #W.R.T OTHER STRATEGIC AGENTS
            #PRICES ARE EQUAL TO THE MARGINAL PRICES
            JuMP.set_value(c_g_st[t], dat_str.GEN_STR.MC_g[t])
            JuMP.set_value(c_w_st[t], dat_str.WIND_STR.MC_g[t])
        end
        #########ADDING VARS###########
        @variables(EPEC, begin
            #REGULARIZATION A PARAMETER
            ϵ1 == 0.0001
            ϵ2 == 0.0001
            #SLACKS FOR REGULARIZATION
            0 <= r_dch_st[t in t_steps]
            0 <= r_ch_st[t in t_steps]
            0 <= r_g_st[t in t_steps]
            0 <= r_w_st[t in t_steps]
            #THE ASSOCIATED DUALS
            0 <= γ⁺ch_st[j in str_agents, t in t_steps] #Dual of the upper bound of the str charging
            0 <= γ⁻ch_st[j in str_agents, t in t_steps] #Dual of the lower bound of the str charging
            0 <= γ⁺dch_st[j in str_agents, t in t_steps] #Dual of the upper bound of the str discharging
            0 <= γ⁻dch_st[j in str_agents, t in t_steps] #Dual of the lower bound of the str discharging
            0 <= γ⁺g_st[j in str_agents, t in t_steps]
            0 <= γ⁻g_st[j in str_agents, t in t_steps]
            0 <= γ⁺w_st[j in str_agents, t in t_steps]
            0 <= γ⁻w_st[j in str_agents, t in t_steps]
            ############UL DECISION VARS################
            ############PRICE BIDS###########
            c_dch_st[t in t_steps] #ESS price bid for discharging
            c_ch_st[t in t_steps] #ESS price bid for charging
            #########QUANTITY BIDS############
            0 <= G_ST[t in t_steps] #Generator str quantity bid
            0 <= W_ST[t in t_steps] #Wind str quantity bid
            0 <= DCH_ST[t in t_steps] #ESS discharge quantity bid
            0 <= CH_ST[t in t_steps] #ESS charge quantity bid
            0 <= SOC[t in t_steps] #ESS State of Charge var
            #############LL DECISION VARS###########
            ######PRIMALS#######
            0 <= g_st[t in t_steps] #Dispatched generation of the strategic agent
            0 <= w_st[t in t_steps] #Dispatched wind of the strategic agent
            0 <= g[i in non_str_gens, t in t_steps] #Dispatched generation
            0 <= w[k in non_str_res, t in t_steps] #Dispatched wind
            0 <= dch_st[t in t_steps] #Dispatched strategic discharge
            0 <= ch_st[t in t_steps] #Dispatched strategic charge
            0 <= d[t in t_steps] #Dispatched demand
            #######DUALS#########
            0 <= β⁺g[i in non_str_gens, t in t_steps] #Dual of the upper bound of the generation
            0 <= β⁻g[i in non_str_gens, t in t_steps] #Dual of the lower bound of the generation
            0 <= β⁺w[k in non_str_res, t in t_steps] #Dual of the upper bound of the RE generation
            0 <= β⁻w[k in non_str_res, t in t_steps] #Dual of the lower bound of the RE generation
            0 <= β⁺d[t in t_steps] #Dual of the upper bound of the demand
            0 <= β⁻d[t in t_steps] #Dual of the lower bound of the demand
            0 <= β⁺g_st[t in t_steps] #Dual of the upper bound of the str generation
            0 <= β⁻g_st[t in t_steps] #Dual of the lower bound of the str generation
            0 <= β⁺w_st[t in t_steps] #Dual of the upper bound of the str RE generation
            0 <= β⁻w_st[t in t_steps] #Dual of the lower bound of the str RE generation
            0 <= β⁺ch_st[t in t_steps] #Dual of the upper bound of the str charging
            0 <= β⁻ch_st[t in t_steps] #Dual of the lower bound of the str charging
            0 <= β⁺dch_st[t in t_steps] #Dual of the upper bound of the str discharging
            0 <= β⁻dch_st[t in t_steps] #Dual of the lower bound of the str discharging
            λ[t in t_steps] #Market price, dual of the energy balance cstr of LL
            ########SECOND STAGE DUALS#############
            0 <= ν⁺g[j in str_agents, i in non_str_gens, t in t_steps]
            0 <= ν⁻g[j in str_agents, i in non_str_gens, t in t_steps]
            0 <= ν⁺w[j in str_agents, k in non_str_res, t in t_steps]
            0 <= ν⁻w[j in str_agents, k in non_str_res, t in t_steps]
            0 <= ν⁺d[j in str_agents, t in t_steps]
            0 <= ν⁻d[j in str_agents, t in t_steps]
            0 <= ν⁺g_st[j in str_agents, t in t_steps]
            0 <= ν⁻g_st[j in str_agents, t in t_steps]
            0 <= ν⁺w_st[j in str_agents, t in t_steps]
            0 <= ν⁻w_st[j in str_agents, t in t_steps]
            0 <= ν⁺ch_st[j in str_agents, t in t_steps]
            0 <= ν⁻ch_st[j in str_agents, t in t_steps]
            0 <= ν⁺dch_st[j in str_agents, t in t_steps]
            0 <= ν⁻dch_st[j in str_agents, t in t_steps]
            #DUAL OF THE STRONG DUALITY EQ DUAL
            0 <= σ_PD[j in str_agents]
            #DUAL FOR THE MARKET PRICE (2nd stage)
            π[j in str_agents, t in t_steps]
            #DUAL OF THE SoC rules
            μsoc[t in t_steps]
            μ_0 #Dual assiciated with the SOC init
            #Duals of the LL equalities for KKTs
            ψg[j in str_agents, i in non_str_gens, t in t_steps]
            ψw[j in str_agents, k in non_str_res, t in t_steps]
            ψd[j in str_agents, t in t_steps]
            #For strategic agents:
            ψch_st[j in str_agents, t in t_steps]
            ψdch_st[j in str_agents, t in t_steps]
            ψg_st[j in str_agents, t in t_steps]
            ψw_st[j in str_agents, t in t_steps]
            0 <= δ⁺g[j in str_agents, i in non_str_gens, t in t_steps]
            0 <= δ⁻g[j in str_agents, i in non_str_gens, t in t_steps]
            0 <= δ⁺w[j in str_agents, k in non_str_res, t in t_steps]
            0 <= δ⁻w[j in str_agents, k in non_str_res, t in t_steps]
            0 <= δ⁺d[j in str_agents, t in t_steps]
            0 <= δ⁻d[j in str_agents, t in t_steps]
            0 <= δ⁺g_st[j in str_agents, t in t_steps]
            0 <= δ⁻g_st[j in str_agents, t in t_steps]
            0 <= δ⁺w_st[j in str_agents, t in t_steps]
            0 <= δ⁻w_st[j in str_agents, t in t_steps]
            0 <= δ⁺ch_st[j in str_agents, t in t_steps]
            0 <= δ⁻ch_st[j in str_agents, t in t_steps]
            0 <= δ⁺dch_st[j in str_agents, t in t_steps]
            0 <= δ⁻dch_st[j in str_agents, t in t_steps]
            #Duals of the UL bounds, specific to each agent
            #FOR THE ESS
            0 <= κ⁺soc[t in t_steps] #Dual of the upper bound of the SoC cstr
            0 <= κ⁻soc[t in t_steps] #Dual of the lower bound of the SoC cstr
            0 <= κ⁺dch[t in t_steps] #Dual of the upper bound of the dch limit cstr
            0 <= κ⁻dch[t in t_steps] #Dual of the lower bound of the dch limit cstr
            0 <= κ⁺ch[t in t_steps] #Dual of the upper bound of the ch limit cstr
            0 <= κ⁻ch[t in t_steps] #Dual of the lower bound of the ch limit cstr
            #FOR THE GEN
            0 <= κ⁺g[t in t_steps] #Dual of the upper bound of the gen limit cstr
            0 <= κ⁻g[t in t_steps] #Dual of the lower bound of the gen limit cstr
            #FOR THE GEN
            0 <= κ⁺w[t in t_steps] #Dual of the upper bound of the gen limit cstr
            0 <= κ⁻w[t in t_steps] #Dual of the lower bound of the gen limit cstr
            # #FOR THE PRICE LIMITS
            0 <= ϵ⁺ch[t in t_steps]
            0 <= ϵ⁻ch[t in t_steps]
            0 <= ϵ⁺dch[t in t_steps]
            0 <= ϵ⁻dch[t in t_steps]
            #############SLACKS#############
            #LL Slacks
            0 <= s_g[i in non_str_gens, t in t_steps]
            0 <= s_w[k in non_str_res, t in t_steps]
            0 <= s_d[t in t_steps]
            0 <= s_ch_st[t in t_steps]
            0 <= s_dch_st[t in t_steps]
            0 <= s_g_st[t in t_steps]
            0 <= s_w_st[t in t_steps]
            #UL slacks
            #Quantity bounds
            0 <= t_g_bid[t in t_steps]
            0 <= t_w_bid[t in t_steps]
            0 <= t_ch_bid[t in t_steps]
            0 <= t_dch_bid[t in t_steps]
            # #Price bounds
            0 <= t⁺ch[t in t_steps]
            0 <= t⁻ch[t in t_steps]
            0 <= t⁺dch[t in t_steps]
            0 <= t⁻dch[t in t_steps]
            #SoC
            0 <= t_soc[t in t_steps]
            #THe param for regularization
        end)

        ##########ENERGY BALANCE##########
        eq[:en_balance] = balance_eq = @expression(EPEC, [t in t_steps],
        - g_st[t] - w_st[t] -dch_st[t] + ch_st[t]
        - sum(g[i,t] for i in non_str_gens) - sum(w[k,t] for k in non_str_res) + d[t])

        #######STRONG DUALITY EQ##########
        #PRIMAL OBJ - DUAL OBJ
        eq[:SDE_eq] = SDE_eq = @NLexpression(EPEC,
                            sum(
                            g_st[t] * c_g_st[t] +
                            w_st[t] * c_w_st[t] +
                            dch_st[t] * c_dch_st[t] +
                            - ch_st[t] * c_ch_st[t] +
                            sum(w[k,t] * c_w[k,t] for k in non_str_res) +
                            sum(g[i,t] * c_g[i,t] for i in non_str_gens) -
                            d[t] * c_d[t] +
                            sum(β⁺g[i,t] * G[i,t] for i in non_str_gens) +
                            sum(β⁺w[k,t] * W[k,t] for k in non_str_res) +
                            r_g_st[t] +
                            r_w_st[t] +
                            r_ch_st[t] +
                            r_dch_st[t] +
                            β⁺d[t] * D[t]
                            for t in t_steps)
        )

        ###########COMMON CONSTRAINTS#########
        #NOTE: MOVE HERE THE SHARED CSTRS
        ############PRIMAL FEASIBILITY OF THE MPPDC#####################
        @constraints(EPEC, begin
            # ENERGY BALANCE CSTR
            en_balance[t in t_steps], balance_eq[t] == 0
            #SOC RULES
            #NOTE: SOC RULES ARE NOT COMMON TO ALL AGENT BUT THERE IS NO DIFFERENTIATION IN J VARS
            soc_rule_0[t = 1], SOC[end] - SOC[t] + ch_st[t] * ηch - (dch_st[t]/ηdch) == 0
            soc_rule[t = 2:t_steps[end]], SOC[t-1] - SOC[t] + ch_st[t] * ηch - (dch_st[t]/ηdch) == 0
            soc_init, SOC[end] == soc_init
        end)

        @NLconstraints(EPEC, begin
            #ENFORCING STRONG DUALITY
            SDE_eq == 0
            #LIMITS OF THE DISPATCHED QUANTITIES OF THE STRATEGIC AGENTS, THEY CANT BE SET FROM OUTSIDE THE MODEL
            [t in t_steps], g_st[t] <= G_ST[t]
            [t in t_steps], w_st[t] <= W_ST[t]
            [t in t_steps], dch_st[t] <= DCH_ST[t]
            [t in t_steps], ch_st[t] <= CH_ST[t]
            [t in t_steps], SOC[t] <= E_max[t]
            #NOTE: THIS IS STILL NEEDED TO MAKE SURE CS WORKS AS A MULTIPLICATION
            [t in t_steps], CH_ST[t] <= P_CH_ST[t]
            [t in t_steps], DCH_ST[t] <= P_DCH_ST[t]
            [t in t_steps], G_ST[t] <= P_G_ST[t]
            [t in t_steps], W_ST[t] <= P_W_ST[t]
            #FOR THE NON STR AGENTS
            [i in non_str_gens, t in t_steps], g[i,t] <= G[i,t]
            [k in non_str_res, t in t_steps], w[k,t] <= W[k,t]

            [t in t_steps], d[t] <= D[t]
            #LAGRANGIAN STATIONARITY
            #w.r.t g_{i,t}:
            [i in non_str_gens, t in t_steps], c_g[i,t] - λ[t] - β⁻g[i,t] + β⁺g[i,t] == 0
            #w.r.t g_st_{i,t}:
            [t in t_steps], c_g_st[t] - λ[t] - β⁻g_st[t] + β⁺g_st[t] == 0
            #w.r.t w_{k,t}:
            [k in non_str_res, t in t_steps], c_w[k,t] - λ[t] - β⁻w[k,t] + β⁺w[k,t] == 0
            #w.r.t w_st_{t}:
            [t in t_steps], c_w_st[t] - λ[t] - β⁻w_st[t] + β⁺w_st[t] == 0
            #w.r.t dch_st_{t}:
            [t in t_steps], c_dch_st[t] - λ[t] - β⁻dch_st[t] + β⁺dch_st[t] == 0
            #w.r.t ch_st_{t}:
            [t in t_steps], -c_ch_st[t] + λ[t] - β⁻ch_st[t] + β⁺ch_st[t] == 0
            #w.r.t d_{t}:
            [t in t_steps], -c_d[t] + λ[t] - β⁻d[t] + β⁺d[t] == 0

            #NOTE: ESS CAN RESPECT IT BY USING ECONOMIC SENSE, FOR OTHERS NEEDED
            #EXCLUSIVE CH/DCH
            # excl[t in t_steps], dch_st[t] * ch_st[t] == 0
            #########ADDING THE SLACKS TO THE CS CONSTRAINTS#######
            #INHERITING FROM LL
            [i in non_str_gens, t in t_steps], G[i,t] -  g[i,t] - s_g[i,t] == 0
            [t in t_steps], G_ST[t] -  g_st[t] - s_g_st[t] == 0
            [k in non_str_res, t in t_steps], W[k,t] -  w[k,t] - s_w[k,t] == 0
            [t in t_steps], W_ST[t] -  w_st[t] - s_w_st[t] == 0
            [t in t_steps], DCH_ST[t] -  dch_st[t] - s_dch_st[t] == 0
            [t in t_steps], CH_ST[t] -  ch_st[t] - s_ch_st[t] == 0
            [t in t_steps], D[t] -  d[t] - s_d[t] == 0
        end)

        #################FOR THE ESS##########
            j = "ESS_STR"

            ##############CONSTRAINTS INHERITING FROM SECOND-STAGE#########
            @NLconstraints(EPEC, begin
                ############LAGRANGIAN STATIONARITY OF THE SECOND STAGE KKTs#########
                #W.R.T. PRIMALS
                #w.r.t c_dch_st{t}:
                [t in t_steps], - ϵ⁻dch[t] + ϵ⁺dch[t] + σ_PD[j] * dch_st[t] + ψdch_st[j,t] == 0
                #w.r.t c_ch_st{t}:
                [t in t_steps], - ϵ⁻ch[t] + ϵ⁺ch[t] - σ_PD[j] * ch_st[t] - ψch_st[j,t] == 0
                #w.r.t DCH_st{t}:
                [t in t_steps],  - κ⁻dch[t] + κ⁺dch[t] - ν⁺dch_st[j,t] - γ⁺dch_st[j,t] * β⁺dch_st[t] + γ⁻dch_st[j,t] * β⁺dch_st[t]  == 0
                #w.r.t CH_st{t}:
                [t in t_steps],  - κ⁻ch[t] + κ⁺ch[t] - ν⁺ch_st[j,t] - γ⁺ch_st[j,t] * β⁺ch_st[t] + γ⁻ch_st[j,t] * β⁺ch_st[t] == 0
                #w.r.t dch_st{t}:
                [t in t_steps],  - π[j,t] - ν⁻dch_st[j,t] + ν⁺dch_st[j,t] + σ_PD[j] * c_dch_st[t] - (μsoc[t]/ηdch)  == 0
                #w.r.t ch_st{t}:
                [t in t_steps],  + π[j,t] - ν⁻ch_st[j,t] + ν⁺ch_st[j,t] - σ_PD[j] * c_ch_st[t] + (μsoc[t]*ηch) == 0
                #w.r.t g_st{t}:
                [t in t_steps], + c_g_st[t] - π[j,t] - ν⁻g_st[j,t] + ν⁺g_st[j,t] + σ_PD[j] * c_g_st[t] == 0
                #w.r.t g{t}:
                [i in non_str_gens, t in t_steps], + c_g[i,t] - π[j,t] - ν⁻g[j,i,t] + ν⁺g[j,i,t] + σ_PD[j] * c_g[i,t] == 0
                #w.r.t w_st{t}:
                [t in t_steps], + c_w_st[t] - π[j,t] - ν⁻w_st[j,t] + ν⁺w_st[j,t] + σ_PD[j] * c_w_st[t] == 0
                #w.r.t w{t}:
                [k in non_str_res, t in t_steps], + c_w[k,t] - π[j,t] - ν⁻w[j,k,t] + ν⁺w[j,k,t] + σ_PD[j] * c_w[k,t] == 0
                #w.r.t d{t}:
                [t in t_steps], - c_d[t] + π[j,t] - ν⁻d[j,t] + ν⁺d[j,t] - σ_PD[j] * c_d[t] == 0
                #w.r.t SOC{t}:
                [t = t_steps[end]], + μsoc[1] - μsoc[t] - κ⁻soc[t] + κ⁺soc[t] + μ_0 == 0
                [t in 1:t_steps[end-1]], + μsoc[t+1] - μsoc[t] - κ⁻soc[t] + κ⁺soc[t] == 0
                #W.R.T. DUALS
                #w.r.t λ{t}:
                [t in t_steps],  + ψd[j,t] - ψg_st[j,t]  - ψw_st[j,t] - sum(ψg[j,i,t] for i in non_str_gens) - sum(ψw[j,k,t] for k in non_str_res) + ψch_st[j,t] - ψdch_st[j,t] == 0
                #w.r.t β⁺g_st{t}:
                [t in t_steps], - γ⁺g_st[j,t] * G_ST[t] + γ⁻g_st[j,t] * G_ST[t] + ψg_st[j,t] - δ⁺g_st[j,t] == 0
                #w.r.t β⁻g_st{t}:
                [t in t_steps], - ψg_st[j,t] - δ⁻g_st[j,t] == 0
                #w.r.t β⁺g{i,t}:
                [i in non_str_gens, t in t_steps], + G[i,t] + σ_PD[j] * G[i,t] + ψg[j,i,t] - δ⁺g[j,i,t] == 0
                #w.r.t β⁻g{i,t}:
                [i in non_str_gens, t in t_steps], - ψg[j,i,t] - δ⁻g[j,i,t] == 0
                #w.r.t β⁺w_st{t}:
                [t in t_steps], - γ⁺w_st[j,t] * W_ST[t] + γ⁻w_st[j,t] * W_ST[t] + ψw_st[j,t] - δ⁺w_st[j,t] == 0
                #w.r.t β⁻w_st{t}:
                [t in t_steps], - ψw_st[j,t] - δ⁻w_st[j,t] == 0
                #w.r.t β⁺w{k,t}:
                [k in non_str_res, t in t_steps], + W[k,t] + σ_PD[j] * W[k,t] + ψw[j,k,t] - δ⁺w[j,k,t] == 0
                #w.r.t β⁻w{k,t}:
                [k in non_str_res, t in t_steps], - ψw[j,k,t] - δ⁻w[j,k,t] == 0
                #w.r.t β⁺dch_st{t}:
                [t in t_steps], - γ⁺dch_st[j,t] * DCH_ST[t] + γ⁻dch_st[j,t] * DCH_ST[t] + ψdch_st[j,t] - δ⁺dch_st[j,t] == 0
                #w.r.t β⁻dch_st{t}:
                [t in t_steps], - ψdch_st[j,t] - δ⁻dch_st[j,t] == 0
                #w.r.t β⁺ch_st{t}:
                [t in t_steps], - γ⁺ch_st[j,t] * CH_ST[t] + γ⁻ch_st[j,t] * CH_ST[t] + ψch_st[j,t] - δ⁺ch_st[j,t] == 0
                #w.r.t β⁻ch_st{t}:
                [t in t_steps], - ψch_st[j,t] - δ⁻ch_st[j,t] == 0
                #w.r.t β⁺d{t}:
                [t in t_steps], + D[t] + σ_PD[j] * D[t] + ψd[j,t] - δ⁺d[j,t] == 0
                #w.r.t β⁻d{t}:
                [t in t_steps], - ψd[j,t] - δ⁻d[j,t] == 0

                #w.r.t. s_dch:
                [t in t_steps],  + γ⁺dch_st[j,t] - γ⁻dch_st[j,t] + σ_PD[j] == 0
                #w.r.t. s_ch:
                [t in t_steps],  + γ⁺ch_st[j,t] - γ⁻ch_st[j,t] + σ_PD[j] == 0
                #w.r.t. s_g:
                [t in t_steps],  + γ⁺g_st[j,t] - γ⁻g_st[j,t] + σ_PD[j] == 0
                #w.r.t. s_w:
                [t in t_steps],  + γ⁺w_st[j,t] - γ⁻w_st[j,t] + σ_PD[j] == 0

                #########ADDING THE SLACKS TO THE CS CONSTRAINTS#######
                #INHERITING FROM UL
                #FROM QUANTITY BOUNDS
                [t in t_steps], P_CH_ST[t] - CH_ST[t] - t_ch_bid[t] == 0
                [t in t_steps], P_DCH_ST[t] - DCH_ST[t] - t_dch_bid[t] == 0
                [t in t_steps], E_max[t] - SOC[t] - t_soc[t] == 0
                #FROM PRICE LIMITS
                [t in t_steps], c_ch_st[t] - p_floor - t⁻ch[t] == 0
                [t in t_steps], p_cap -  c_ch_st[t] - t⁺ch[t] == 0
                [t in t_steps], c_dch_st[t] - p_floor - t⁻dch[t] == 0
                [t in t_steps], p_cap -  c_dch_st[t] - t⁺dch[t] == 0

            end)
            ############COMPLEMENTARITY SLACKNESS##########
            @NLconstraints(EPEC, begin
                #RESULTING FROM THE RELAXATION
                [t in t_steps], β⁺ch_st[t] * CH_ST[t] - r_ch_st[t] <= ϵ1
                [t in t_steps], -β⁺ch_st[t] * CH_ST[t] + r_ch_st[t] <= ϵ1
                [t in t_steps], β⁺dch_st[t] * DCH_ST[t] - r_dch_st[t]  <= ϵ1
                [t in t_steps], -β⁺dch_st[t] * DCH_ST[t] + r_dch_st[t]  <= ϵ1
                [t in t_steps], β⁺g_st[t] * G_ST[t] - r_g_st[t] <= ϵ1
                [t in t_steps], -β⁺g_st[t] * G_ST[t] + r_g_st[t] <= ϵ1
                [t in t_steps], β⁺w_st[t] * W_ST[t] - r_w_st[t] <= ϵ1
                [t in t_steps], -β⁺w_st[t] * W_ST[t] + r_w_st[t] <= ϵ1
                #CS PART
                [t in t_steps], (ϵ1 - (β⁺ch_st[t] * CH_ST[t] - r_ch_st[t]) * γ⁻ch_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (-β⁺ch_st[t] * CH_ST[t] + r_ch_st[t]) * γ⁺ch_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (β⁺dch_st[t] * DCH_ST[t] - r_dch_st[t]) * γ⁻dch_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (-β⁺dch_st[t] * DCH_ST[t] + r_dch_st[t]) * γ⁺dch_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (β⁺g_st[t] * G_ST[t] - r_g_st[t]) * γ⁻g_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (-β⁺g_st[t] * G_ST[t] + r_g_st[t]) * γ⁺g_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (β⁺w_st[t] * W_ST[t] - r_w_st[t]) * γ⁻w_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (-β⁺w_st[t] * W_ST[t] + r_w_st[t]) * γ⁺w_st[j,t]) <= ϵ2

                #LL related
                [i in non_str_gens, t in t_steps], g[i,t] * ν⁻g[j,i,t]  <= ϵ2
                [i in non_str_gens, t in t_steps], s_g[i,t] * ν⁺g[j,i,t]  <= ϵ2
                [i in non_str_gens, t in t_steps], β⁻g[i,t] *  δ⁻g[j,i,t]  <= ϵ2
                [i in non_str_gens, t in t_steps], β⁺g[i,t] *  δ⁺g[j,i,t] <= ϵ2
                [k in non_str_res, t in t_steps], w[k,t] * ν⁻w[j,k,t]  <= ϵ2
                [k in non_str_res, t in t_steps], s_w[k,t] * ν⁺w[j,k,t]  <= ϵ2
                [k in non_str_res, t in t_steps], β⁻w[k,t] *  δ⁻w[j,k,t]  <= ϵ2
                [k in non_str_res, t in t_steps], β⁺w[k,t] *  δ⁺w[j,k,t] <= ϵ2
                [t in t_steps], g_st[t] * ν⁻g_st[j,t]  <= ϵ2
                [t in t_steps], s_g_st[t] * ν⁺g_st[j,t]  <= ϵ2
                [t in t_steps], β⁻g_st[t] *  δ⁻g_st[j,t]  <= ϵ2
                [t in t_steps], β⁺g_st[t] *  δ⁺g_st[j,t] <= ϵ2
                [t in t_steps], w_st[t] * ν⁻w_st[j,t]  <= ϵ2
                [t in t_steps], s_w_st[t] * ν⁺w_st[j,t]  <= ϵ2
                [t in t_steps], β⁻w_st[t] *  δ⁻w_st[j,t]  <= ϵ2
                [t in t_steps], β⁺w_st[t] *  δ⁺w_st[j,t] <= ϵ2
                [t in t_steps], dch_st[t] * ν⁻dch_st[j,t]  <= ϵ2
                [t in t_steps], s_dch_st[t] * ν⁺dch_st[j,t]  <= ϵ2
                [t in t_steps], β⁻dch_st[t] *  δ⁻dch_st[j,t]  <= ϵ2
                [t in t_steps], β⁺dch_st[t] *  δ⁺dch_st[j,t] <= ϵ2
                [t in t_steps], ch_st[t] * ν⁻ch_st[j,t]  <= ϵ2
                [t in t_steps], s_ch_st[t] * ν⁺ch_st[j,t]  <= ϵ2
                [t in t_steps], β⁻ch_st[t] *  δ⁻ch_st[j,t]  <= ϵ2
                [t in t_steps], β⁺ch_st[t] *  δ⁺ch_st[j,t] <= ϵ2
                [t in t_steps], d[t] * ν⁻d[j,t]  <= ϵ2
                [t in t_steps], s_d[t] * ν⁺d[j,t]  <= ϵ2
                [t in t_steps], β⁻d[t] *  δ⁻d[j,t]  <= ϵ2
                [t in t_steps], β⁺d[t] *  δ⁺d[j,t] <= ϵ2
                #UL related
                #Quantity related
                [t in t_steps], t_ch_bid[t] * κ⁺ch[t]  <= ϵ2
                [t in t_steps], DCH_ST[t] * κ⁻ch[t]  <= ϵ2
                [t in t_steps], t_dch_bid[t] * κ⁺dch[t]  <= ϵ2
                [t in t_steps], CH_ST[t] * κ⁻dch[t]  <= ϵ2
                [t in t_steps], t_soc[t] * κ⁺soc[t]  <= ϵ2
                [t in t_steps], SOC[t] * κ⁻soc[t] <= ϵ2
                #Price related
                [t in t_steps], t⁻ch[t] * ϵ⁻ch[t]  <= ϵ2
                [t in t_steps], t⁺ch[t] * ϵ⁺ch[t]  <= ϵ2
                [t in t_steps], t⁻dch[t] * ϵ⁻dch[t]  <= ϵ2
                [t in t_steps], t⁺dch[t] * ϵ⁺dch[t]  <= ϵ2
            end)

        #################FOR THE GENERATOR##########
            j = "GEN_STR"

            ##############CONSTRAINTS INHERITING FROM SECOND-STAGE#########
            @NLconstraints(EPEC, begin
                ############LAGRANGIAN STATIONARITY OF THE SECOND STAGE KKTs#########
                #W.R.T. PRIMALS
                #w.r.t G_st{t}:
                [t in t_steps],  - κ⁻g[t] + κ⁺g[t] - ν⁺g_st[j,t] - γ⁺g_st[j,t] * β⁺g_st[t] + γ⁻g_st[j,t] * β⁺g_st[t] == 0
                #w.r.t dch_st{t}:
                [t in t_steps], + c_dch_st[t] - π[j,t] - ν⁻dch_st[j,t] + ν⁺dch_st[j,t] + σ_PD[j] * c_dch_st[t]  == 0
                #w.r.t ch_st{t}:
                [t in t_steps], - c_ch_st[t]  + π[j,t] - ν⁻ch_st[j,t] + ν⁺ch_st[j,t] - σ_PD[j] * c_ch_st[t] == 0
                #w.r.t g_st{t}:
                [t in t_steps], + dat_str.GEN_STR.MC_g[t] - π[j,t] - ν⁻g_st[j,t] + ν⁺g_st[j,t] + σ_PD[j] * c_g_st[t] == 0
                #w.r.t g{t}:
                [i in non_str_gens, t in t_steps], + c_g[i,t] - π[j,t] - ν⁻g[j,i,t] + ν⁺g[j,i,t] + σ_PD[j] * c_g[i,t] == 0
                #w.r.t w_st{t}:
                [t in t_steps], + c_w_st[t] - π[j,t] - ν⁻w_st[j,t] + ν⁺w_st[j,t] + σ_PD[j] * c_w_st[t] == 0
                #w.r.t w{t}:
                [k in non_str_res, t in t_steps], + c_w[k,t] - π[j,t] - ν⁻w[j,k,t] + ν⁺w[j,k,t] + σ_PD[j] * c_w[k,t] == 0
                #w.r.t d{t}:
                [t in t_steps], - c_d[t] + π[j,t] - ν⁻d[j,t] + ν⁺d[j,t] - σ_PD[j] * c_d[t] == 0
                #W.R.T. DUALS
                #w.r.t λ{t}:
                [t in t_steps],  + ψd[j,t] - ψg_st[j,t]  - ψw_st[j,t] - sum(ψg[j,i,t] for i in non_str_gens) - sum(ψw[j,k,t] for k in non_str_res) + ψch_st[j,t] - ψdch_st[j,t] == 0
                #w.r.t β⁺g_st{t}:
                [t in t_steps], - γ⁺g_st[j,t] * G_ST[t] + γ⁻g_st[j,t] * G_ST[t] + ψg_st[j,t] - δ⁺g_st[j,t] == 0
                #w.r.t β⁻g_st{t}:
                [t in t_steps], - ψg_st[j,t] - δ⁻g_st[j,t] == 0
                #w.r.t β⁺g{i,t}:
                [i in non_str_gens, t in t_steps], + G[i,t] + σ_PD[j] * G[i,t] + ψg[j,i,t] - δ⁺g[j,i,t] == 0
                #w.r.t β⁻g{i,t}:
                [i in non_str_gens, t in t_steps], - ψg[j,i,t] - δ⁻g[j,i,t] == 0
                #w.r.t β⁺w_st{t}:
                [t in t_steps], - γ⁺w_st[j,t] * W_ST[t] + γ⁻w_st[j,t] * W_ST[t] + ψw_st[j,t] - δ⁺w_st[j,t] == 0
                #w.r.t β⁻w_st{t}:
                [t in t_steps], - ψw_st[j,t] - δ⁻w_st[j,t] == 0
                #w.r.t β⁺w{k,t}:
                [k in non_str_res, t in t_steps], + W[k,t] + σ_PD[j] * W[k,t] + ψw[j,k,t] - δ⁺w[j,k,t] == 0
                #w.r.t β⁻w{k,t}:
                [k in non_str_res, t in t_steps], - ψw[j,k,t] - δ⁻w[j,k,t] == 0
                #w.r.t β⁺dch_st{t}:
                [t in t_steps], - γ⁺dch_st[j,t] * DCH_ST[t] + γ⁻dch_st[j,t] * DCH_ST[t] + ψdch_st[j,t] - δ⁺dch_st[j,t] == 0
                #w.r.t β⁻dch_st{t}:
                [t in t_steps], - ψdch_st[j,t] - δ⁻dch_st[j,t] == 0
                #w.r.t β⁺ch_st{t}:
                [t in t_steps], - γ⁺ch_st[j,t] * CH_ST[t] + γ⁻ch_st[j,t] * CH_ST[t] + ψch_st[j,t] - δ⁺ch_st[j,t] == 0
                #w.r.t β⁻ch_st{t}:
                [t in t_steps], - ψch_st[j,t] - δ⁻ch_st[j,t] == 0
                #w.r.t β⁺d{t}:
                [t in t_steps], + D[t] + σ_PD[j] * D[t] + ψd[j,t] - δ⁺d[j,t] == 0
                #w.r.t β⁻d{t}:
                [t in t_steps], - ψd[j,t] - δ⁻d[j,t] == 0

                #w.r.t. s_dch:
                [t in t_steps],  + γ⁺dch_st[j,t] - γ⁻dch_st[j,t] + σ_PD[j] == 0
                #w.r.t. s_ch:
                [t in t_steps],  + γ⁺ch_st[j,t] - γ⁻ch_st[j,t] + σ_PD[j] == 0
                #w.r.t. s_g:
                [t in t_steps],  + γ⁺g_st[j,t] - γ⁻g_st[j,t] + σ_PD[j] == 0
                #w.r.t. s_w:
                [t in t_steps],  + γ⁺w_st[j,t] - γ⁻w_st[j,t] + σ_PD[j] == 0
                #########ADDING THE SLACKS TO THE CS CONSTRAINTS#######
                #INHERITING FROM UL
                #FROM QUANTITY BOUNDS
                [t in t_steps], P_G_ST[t] - G_ST[t] - t_g_bid[t] == 0

            end)
            ############COMPLEMENTARITY SLACKNESS##########
            @NLconstraints(EPEC, begin
                #RESULTING FROM THE RELAXATION
                [t in t_steps], β⁺ch_st[t] * CH_ST[t] - r_ch_st[t] <= ϵ1
                [t in t_steps], -β⁺ch_st[t] * CH_ST[t] + r_ch_st[t] <= ϵ1
                [t in t_steps], β⁺dch_st[t] * DCH_ST[t] - r_dch_st[t]  <= ϵ1
                [t in t_steps], -β⁺dch_st[t] * DCH_ST[t] + r_dch_st[t]  <= ϵ1
                [t in t_steps], β⁺g_st[t] * G_ST[t] - r_g_st[t] <= ϵ1
                [t in t_steps], -β⁺g_st[t] * G_ST[t] + r_g_st[t] <= ϵ1
                [t in t_steps], β⁺w_st[t] * W_ST[t] - r_w_st[t] <= ϵ1
                [t in t_steps], -β⁺w_st[t] * W_ST[t] + r_w_st[t] <= ϵ1
                #CS PART
                [t in t_steps], (ϵ1 - (β⁺ch_st[t] * CH_ST[t] - r_ch_st[t]) * γ⁻ch_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (-β⁺ch_st[t] * CH_ST[t] + r_ch_st[t]) * γ⁺ch_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (β⁺dch_st[t] * DCH_ST[t] - r_dch_st[t]) * γ⁻dch_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (-β⁺dch_st[t] * DCH_ST[t] + r_dch_st[t]) * γ⁺dch_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (β⁺g_st[t] * G_ST[t] - r_g_st[t]) * γ⁻g_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (-β⁺g_st[t] * G_ST[t] + r_g_st[t]) * γ⁺g_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (β⁺w_st[t] * W_ST[t] - r_w_st[t]) * γ⁻w_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (-β⁺w_st[t] * W_ST[t] + r_w_st[t]) * γ⁺w_st[j,t]) <= ϵ2

                #LL related
                [i in non_str_gens, t in t_steps], g[i,t] * ν⁻g[j,i,t]  <= ϵ2
                [i in non_str_gens, t in t_steps], s_g[i,t] * ν⁺g[j,i,t]  <= ϵ2
                [i in non_str_gens, t in t_steps], β⁻g[i,t] *  δ⁻g[j,i,t]  <= ϵ2
                [i in non_str_gens, t in t_steps], β⁺g[i,t] *  δ⁺g[j,i,t] <= ϵ2
                [k in non_str_res, t in t_steps], w[k,t] * ν⁻w[j,k,t]  <= ϵ2
                [k in non_str_res, t in t_steps], s_w[k,t] * ν⁺w[j,k,t]  <= ϵ2
                [k in non_str_res, t in t_steps], β⁻w[k,t] *  δ⁻w[j,k,t]  <= ϵ2
                [k in non_str_res, t in t_steps], β⁺w[k,t] *  δ⁺w[j,k,t] <= ϵ2
                [t in t_steps], g_st[t] * ν⁻g_st[j,t]  <= ϵ2
                [t in t_steps], s_g_st[t] * ν⁺g_st[j,t]  <= ϵ2
                [t in t_steps], β⁻g_st[t] *  δ⁻g_st[j,t]  <= ϵ2
                [t in t_steps], β⁺g_st[t] *  δ⁺g_st[j,t] <= ϵ2
                [t in t_steps], w_st[t] * ν⁻w_st[j,t]  <= ϵ2
                [t in t_steps], s_w_st[t] * ν⁺w_st[j,t]  <= ϵ2
                [t in t_steps], β⁻w_st[t] *  δ⁻w_st[j,t]  <= ϵ2
                [t in t_steps], β⁺w_st[t] *  δ⁺w_st[j,t] <= ϵ2
                [t in t_steps], dch_st[t] * ν⁻dch_st[j,t]  <= ϵ2
                [t in t_steps], s_dch_st[t] * ν⁺dch_st[j,t]  <= ϵ2
                [t in t_steps], β⁻dch_st[t] *  δ⁻dch_st[j,t]  <= ϵ2
                [t in t_steps], β⁺dch_st[t] *  δ⁺dch_st[j,t] <= ϵ2
                [t in t_steps], ch_st[t] * ν⁻ch_st[j,t]  <= ϵ2
                [t in t_steps], s_ch_st[t] * ν⁺ch_st[j,t]  <= ϵ2
                [t in t_steps], β⁻ch_st[t] *  δ⁻ch_st[j,t]  <= ϵ2
                [t in t_steps], β⁺ch_st[t] *  δ⁺ch_st[j,t] <= ϵ2
                [t in t_steps], d[t] * ν⁻d[j,t]  <= ϵ2
                [t in t_steps], s_d[t] * ν⁺d[j,t]  <= ϵ2
                [t in t_steps], β⁻d[t] *  δ⁻d[j,t]  <= ϵ2
                [t in t_steps], β⁺d[t] *  δ⁺d[j,t] <= ϵ2
                #UL related
                #Quantity related
                [t in t_steps], t_g_bid[t] * κ⁺g[t]  <= ϵ2
                [t in t_steps], G_ST[t] * κ⁻g[t]  <= ϵ2
            end)

        #################FOR THE WIND##########
            j = "WIND_STR"
            ##############CONSTRAINTS INHERITING FROM SECOND-STAGE#########
            @NLconstraints(EPEC, begin
                ############LAGRANGIAN STATIONARITY OF THE SECOND STAGE KKTs#########
                #W.R.T. PRIMALS
                #w.r.t W_st{t}:
                [t in t_steps],  - κ⁻w[t] + κ⁺w[t] - ν⁺w_st[j,t] - γ⁺w_st[j,t] * β⁺w_st[t] + γ⁻w_st[j,t] * β⁺w_st[t] == 0
                #w.r.t dch_st{t}:
                [t in t_steps], + c_dch_st[t]  - π[j,t] - ν⁻dch_st[j,t] + ν⁺dch_st[j,t] + σ_PD[j] * c_dch_st[t]  == 0
                #w.r.t ch_st{t}:
                [t in t_steps], - c_ch_st[t]  + π[j,t] - ν⁻ch_st[j,t] + ν⁺ch_st[j,t] - σ_PD[j] * c_ch_st[t] == 0
                #w.r.t g_st{t}:
                [t in t_steps], + c_g_st[t]  - π[j,t] - ν⁻g_st[j,t] + ν⁺g_st[j,t] + σ_PD[j] * c_g_st[t] == 0
                #w.r.t g{t}:
                [i in non_str_gens, t in t_steps], + c_g[i,t] - π[j,t] - ν⁻g[j,i,t] + ν⁺g[j,i,t] + σ_PD[j] * c_g[i,t] == 0
                #w.r.t w_st{t}:
                [t in t_steps], + dat_str.WIND_STR.MC_g[t] - π[j,t] - ν⁻w_st[j,t] + ν⁺w_st[j,t] + σ_PD[j] * c_w_st[t] == 0
                #w.r.t w{t}:
                [k in non_str_res, t in t_steps], + c_w[k,t] - π[j,t] - ν⁻w[j,k,t] + ν⁺w[j,k,t] + σ_PD[j] * c_w[k,t] == 0
                #w.r.t d{t}:
                [t in t_steps], - c_d[t] + π[j,t] - ν⁻d[j,t] + ν⁺d[j,t] - σ_PD[j] * c_d[t] == 0
                #W.R.T. DUALS
                #w.r.t λ{t}:
                [t in t_steps], + ψd[j,t] - ψg_st[j,t]  - ψw_st[j,t] - sum(ψg[j,i,t] for i in non_str_gens) - sum(ψw[j,k,t] for k in non_str_res) + ψch_st[j,t] - ψdch_st[j,t] == 0
                #w.r.t β⁺g_st{t}:
                [t in t_steps], - γ⁺g_st[j,t] * G_ST[t] + γ⁻g_st[j,t] * G_ST[t] + ψg_st[j,t] - δ⁺g_st[j,t] == 0
                #w.r.t β⁻g_st{t}:
                [t in t_steps], - ψg_st[j,t] - δ⁻g_st[j,t] == 0
                #w.r.t β⁺g{i,t}:
                [i in non_str_gens, t in t_steps], + G[i,t] + σ_PD[j] * G[i,t] + ψg[j,i,t] - δ⁺g[j,i,t] == 0
                #w.r.t β⁻g{i,t}:
                [i in non_str_gens, t in t_steps], - ψg[j,i,t] - δ⁻g[j,i,t] == 0
                #w.r.t β⁺w_st{t}:
                [t in t_steps], - γ⁺w_st[j,t] * W_ST[t] + γ⁻w_st[j,t] * W_ST[t] + ψw_st[j,t] - δ⁺w_st[j,t] == 0
                #w.r.t β⁻w_st{t}:
                [t in t_steps], - ψw_st[j,t] - δ⁻w_st[j,t] == 0
                #w.r.t β⁺w{k,t}:
                [k in non_str_res, t in t_steps], + W[k,t] + σ_PD[j] * W[k,t] + ψw[j,k,t] - δ⁺w[j,k,t] == 0
                #w.r.t β⁻w{k,t}:
                [k in non_str_res, t in t_steps], - ψw[j,k,t] - δ⁻w[j,k,t] == 0
                #w.r.t β⁺dch_st{t}:
                [t in t_steps], - γ⁺dch_st[j,t] * DCH_ST[t] + γ⁻dch_st[j,t] * DCH_ST[t] + ψdch_st[j,t] - δ⁺dch_st[j,t] == 0
                #w.r.t β⁻dch_st{t}:
                [t in t_steps], - ψdch_st[j,t] - δ⁻dch_st[j,t] == 0
                #w.r.t β⁺ch_st{t}:
                [t in t_steps], - γ⁺ch_st[j,t] * CH_ST[t] + γ⁻ch_st[j,t] * CH_ST[t] + ψch_st[j,t] - δ⁺ch_st[j,t] == 0
                #w.r.t β⁻ch_st{t}:
                [t in t_steps], - ψch_st[j,t] - δ⁻ch_st[j,t] == 0
                #w.r.t β⁺d{t}:
                [t in t_steps], + D[t] + σ_PD[j] * D[t] + ψd[j,t] - δ⁺d[j,t] == 0
                #w.r.t β⁻d{t}:
                [t in t_steps], - ψd[j,t] - δ⁻d[j,t] == 0

                #w.r.t. s_dch:
                [t in t_steps],  + γ⁺dch_st[j,t] - γ⁻dch_st[j,t] + σ_PD[j] == 0
                #w.r.t. s_ch:
                [t in t_steps],  + γ⁺ch_st[j,t] - γ⁻ch_st[j,t] + σ_PD[j] == 0
                #w.r.t. s_g:
                [t in t_steps],  + γ⁺g_st[j,t] - γ⁻g_st[j,t] + σ_PD[j] == 0
                #w.r.t. s_w:
                [t in t_steps],  + γ⁺w_st[j,t] - γ⁻w_st[j,t] + σ_PD[j] == 0

                #########ADDING THE SLACKS TO THE CS CONSTRAINTS#######
                #INHERITING FROM UL
                #FROM QUANTITY BOUNDS
                [t in t_steps], P_W_ST[t] - W_ST[t] - t_w_bid[t] == 0
                #FROM PRICE LIMITS

            end)
            ############COMPLEMENTARITY SLACKNESS##########
            @NLconstraints(EPEC, begin
                #RESULTING FROM THE RELAXATION
                [t in t_steps], β⁺ch_st[t] * CH_ST[t] - r_ch_st[t] <= ϵ1
                [t in t_steps], -β⁺ch_st[t] * CH_ST[t] + r_ch_st[t] <= ϵ1
                [t in t_steps], β⁺dch_st[t] * DCH_ST[t] - r_dch_st[t]  <= ϵ1
                [t in t_steps], -β⁺dch_st[t] * DCH_ST[t] + r_dch_st[t]  <= ϵ1
                [t in t_steps], β⁺g_st[t] * G_ST[t] - r_g_st[t] <= ϵ1
                [t in t_steps], -β⁺g_st[t] * G_ST[t] + r_g_st[t] <= ϵ1
                [t in t_steps], β⁺w_st[t] * W_ST[t] - r_w_st[t] <= ϵ1
                [t in t_steps], -β⁺w_st[t] * W_ST[t] + r_w_st[t] <= ϵ1
                #CS PART
                [t in t_steps], (ϵ1 - (β⁺ch_st[t] * CH_ST[t] - r_ch_st[t]) * γ⁻ch_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (-β⁺ch_st[t] * CH_ST[t] + r_ch_st[t]) * γ⁺ch_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (β⁺dch_st[t] * DCH_ST[t] - r_dch_st[t]) * γ⁻dch_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (-β⁺dch_st[t] * DCH_ST[t] + r_dch_st[t]) * γ⁺dch_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (β⁺g_st[t] * G_ST[t] - r_g_st[t]) * γ⁻g_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (-β⁺g_st[t] * G_ST[t] + r_g_st[t]) * γ⁺g_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (β⁺w_st[t] * W_ST[t] - r_w_st[t]) * γ⁻w_st[j,t]) <= ϵ2
                [t in t_steps], (ϵ1 - (-β⁺w_st[t] * W_ST[t] + r_w_st[t]) * γ⁺w_st[j,t]) <= ϵ2

                #LL related
                [i in non_str_gens, t in t_steps], g[i,t] * ν⁻g[j,i,t]  <= ϵ2
                [i in non_str_gens, t in t_steps], s_g[i,t] * ν⁺g[j,i,t]  <= ϵ2
                [i in non_str_gens, t in t_steps], β⁻g[i,t] *  δ⁻g[j,i,t]  <= ϵ2
                [i in non_str_gens, t in t_steps], β⁺g[i,t] *  δ⁺g[j,i,t] <= ϵ2
                [k in non_str_res, t in t_steps], w[k,t] * ν⁻w[j,k,t]  <= ϵ2
                [k in non_str_res, t in t_steps], s_w[k,t] * ν⁺w[j,k,t]  <= ϵ2
                [k in non_str_res, t in t_steps], β⁻w[k,t] *  δ⁻w[j,k,t]  <= ϵ2
                [k in non_str_res, t in t_steps], β⁺w[k,t] *  δ⁺w[j,k,t] <= ϵ2
                [t in t_steps], g_st[t] * ν⁻g_st[j,t]  <= ϵ2
                [t in t_steps], s_g_st[t] * ν⁺g_st[j,t]  <= ϵ2
                [t in t_steps], β⁻g_st[t] *  δ⁻g_st[j,t]  <= ϵ2
                [t in t_steps], β⁺g_st[t] *  δ⁺g_st[j,t] <= ϵ2
                [t in t_steps], w_st[t] * ν⁻w_st[j,t]  <= ϵ2
                [t in t_steps], s_w_st[t] * ν⁺w_st[j,t]  <= ϵ2
                [t in t_steps], β⁻w_st[t] *  δ⁻w_st[j,t]  <= ϵ2
                [t in t_steps], β⁺w_st[t] *  δ⁺w_st[j,t] <= ϵ2
                [t in t_steps], dch_st[t] * ν⁻dch_st[j,t]  <= ϵ2
                [t in t_steps], s_dch_st[t] * ν⁺dch_st[j,t]  <= ϵ2
                [t in t_steps], β⁻dch_st[t] *  δ⁻dch_st[j,t]  <= ϵ2
                [t in t_steps], β⁺dch_st[t] *  δ⁺dch_st[j,t] <= ϵ2
                [t in t_steps], ch_st[t] * ν⁻ch_st[j,t]  <= ϵ2
                [t in t_steps], s_ch_st[t] * ν⁺ch_st[j,t]  <= ϵ2
                [t in t_steps], β⁻ch_st[t] *  δ⁻ch_st[j,t]  <= ϵ2
                [t in t_steps], β⁺ch_st[t] *  δ⁺ch_st[j,t] <= ϵ2
                [t in t_steps], d[t] * ν⁻d[j,t]  <= ϵ2
                [t in t_steps], s_d[t] * ν⁺d[j,t]  <= ϵ2
                [t in t_steps], β⁻d[t] *  δ⁻d[j,t]  <= ϵ2
                [t in t_steps], β⁺d[t] *  δ⁺d[j,t] <= ϵ2
                #UL related
                #Quantity related
                [t in t_steps], t_w_bid[t] * κ⁺w[t]  <= ϵ2
                [t in t_steps], W_ST[t] * κ⁻w[t]  <= ϵ2
            end)

        return EPEC
    end

    EPEC = init_EPEC_quantity()
    ##########GETTING THE SP OBJECTIVE

    function get_EPEC_obj()
        if SP_obj == "empty"
            obj = 0
        elseif SP_obj == "competitive"
            #COMPETITIVE OBJECTIVE CORRESPONDS TO THE TRUE SW MAXIMUM
            #CALCULATING SOCIAL WELFARE
            obj =(
                    #Maximize consumers surplus
                    - sum(EPEC[:d][t] .* p_cap for t in t_steps)
                    #Minimize generation cost
                    + sum(EPEC[:g_st][t] * dat_str.GEN_STR.MC_g[t] for t in t_steps) #GEN_STR
                    + sum(EPEC[:w_st][t] * dat_str.WIND_STR.MC_g[t] for t in t_steps) #WIND_STR
                    + sum(sum(EPEC[:g][i,t] * dat_nonstr[i].c_g[t] for t in t_steps) for i in non_str_gens) #GENs
                    + sum(sum(EPEC[:w][k,t] * dat_nonstr[k].c_g[t] for t in t_steps) for k in non_str_res) #GENs
                )

        elseif SP_obj == "fav_demand"
            obj =(
                    #Minimize payment
                    #NOTE: As demand is always fully satisfied this can be linearized by using the maximum demand
                    sum(EPEC[:λ][t] .* load[t] for t in t_steps)
                )

        elseif SP_obj == "collusive"
            #COLLUSIVE EQUILIBRIUM BY MAXIMIZING PRODUCERS PROFITS
            obj = (
            sum(- EPEC[:dch_st][t] * EPEC[:λ][t]
            + EPEC[:ch_st][t] * EPEC[:λ][t] for t in t_steps)
            +
            sum(dat_str.GEN_STR.MC_g[t] * EPEC[:g_st][t] -
            EPEC[:g_st][t] * EPEC[:λ][t] for t in t_steps)
            +
            sum(dat_str.GEN_STR.MC_g[t] * EPEC[:w_st][t] -
            EPEC[:w_st][t] * EPEC[:λ][t] for t in t_steps)
            )

            # obj = (
            #     #ESS
            #     - (- sum(EPEC[:g_st][t] * dat_str.GEN_STR.MC_g[t] for t in t_steps) #GEN_STR
            #     - sum(EPEC[:w_st][t] * dat_str.WIND_STR.MC_g[t] for t in t_steps) #WIND_STR
            #     - sum(sum(EPEC[:g][i,t] * dat_nonstr[i].c_g[t] for t in t_steps) for i in non_str_gens) #GENs
            #     - sum(sum(EPEC[:w][k,t] * dat_nonstr[k].c_g[t] for t in t_steps) for k in non_str_res) #GENs
            #     + sum(EPEC[:d][t] .* p_cap for t in t_steps)
            #     - sum(EPEC[:β⁺g_st][t] * dat_str.GEN_STR.P_max[t] for t in t_steps) #GEN_STR
            #     - sum(EPEC[:β⁺w_st][t] * dat_str.WIND_STR.P_max[t] for t in t_steps) #WIND_STR
            #     - sum(sum(EPEC[:β⁺g][i,t] * dat_nonstr[i].G[t] for t in t_steps) for i in non_str_gens) #GENs
            #     - sum(sum(EPEC[:β⁺w][k,t] * dat_nonstr[k].G[t] for t in t_steps) for k in non_str_res) #GENs
            #     - sum(EPEC[:β⁺d][t] * load[t] for t in t_steps)
            #     )
            #     #GEN
            #     - (
            #     # - sum(EPEC[:g_st][t] * dat_str.GEN_STR.MC_g[t] for t in t_steps) #GEN_STR
            #     - sum(EPEC[:dch_st][t] * dat_str.ESS_STR.c_dch[t] for t in t_steps) #DCH_STR
            #     + sum(EPEC[:ch_st][t] * dat_str.ESS_STR.c_ch[t] for t in t_steps) #DCH_STR
            #     - sum(EPEC[:w_st][t] * dat_str.WIND_STR.MC_g[t] for t in t_steps) #WIND_STR
            #     - sum(sum(EPEC[:g][i,t] * dat_nonstr[i].c_g[t] for t in t_steps) for i in non_str_gens) #GENs
            #     - sum(sum(EPEC[:w][k,t] * dat_nonstr[k].c_g[t] for t in t_steps) for k in non_str_res) #GENs
            #     + sum(EPEC[:d][t] .* p_cap for t in t_steps)
            #     # - sum(EPEC[:β⁺g_st][t] * dat_str.GEN_STR.P_max[t] for t in t_steps) #GEN_STR
            #     - sum(EPEC[:β⁺dch_st][t] * dat_str.ESS_STR.P_max[t] for t in t_steps) #WIND_STR
            #     - sum(EPEC[:β⁺ch_st][t] * dat_str.ESS_STR.P_max[t] for t in t_steps) #WIND_STR
            #     - sum(EPEC[:β⁺w_st][t] * dat_str.WIND_STR.P_max[t] for t in t_steps) #WIND_STR
            #     - sum(sum(EPEC[:β⁺g][i,t] * dat_nonstr[i].G[t] for t in t_steps) for i in non_str_gens) #GENs
            #     - sum(sum(EPEC[:β⁺w][k,t] * dat_nonstr[k].G[t] for t in t_steps) for k in non_str_res) #GENs
            #     - sum(EPEC[:β⁺d][t] * load[t] for t in t_steps)
            #     )
            #     #WIND
            #     - (
            #     - sum(EPEC[:g_st][t] * dat_str.GEN_STR.MC_g[t] for t in t_steps) #GEN_STR
            #     - sum(EPEC[:dch_st][t] * dat_str.ESS_STR.c_dch[t] for t in t_steps) #DCH_STR
            #     + sum(EPEC[:ch_st][t] * dat_str.ESS_STR.c_ch[t] for t in t_steps) #DCH_STR
            #     # - sum(EPEC[:w_st][t] * dat_str.WIND_STR.MC_g[t] for t in t_steps) #WIND_STR
            #     - sum(sum(EPEC[:g][i,t] * dat_nonstr[i].c_g[t] for t in t_steps) for i in non_str_gens) #GENs
            #     - sum(sum(EPEC[:w][k,t] * dat_nonstr[k].c_g[t] for t in t_steps) for k in non_str_res) #GENs
            #     + sum(EPEC[:d][t] .* p_cap for t in t_steps)
            #     - sum(EPEC[:β⁺g_st][t] * dat_str.GEN_STR.P_max[t] for t in t_steps) #GEN_STR
            #     - sum(EPEC[:β⁺dch_st][t] * dat_str.ESS_STR.P_max[t] for t in t_steps) #WIND_STR
            #     - sum(EPEC[:β⁺ch_st][t] * dat_str.ESS_STR.P_max[t] for t in t_steps) #WIND_STR
            #     # - sum(EPEC[:β⁺w_st][t] * dat_str.WIND_STR.P_max[t] for t in t_steps) #WIND_STR
            #     - sum(sum(EPEC[:β⁺g][i,t] * dat_nonstr[i].G[t] for t in t_steps) for i in non_str_gens) #GENs
            #     - sum(sum(EPEC[:β⁺w][k,t] * dat_nonstr[k].G[t] for t in t_steps) for k in non_str_res) #GENs
            #     - sum(EPEC[:β⁺d][t] * load[t] for t in t_steps)
            #     )
            # )

        elseif SP_obj == "ESS_fav"
            #THE OBJECTIVE FUNCTION THAT MAXIMIZES THE ESS PROFIT
            obj = sum(- EPEC[:dch_st][t] * EPEC[:λ][t]
            + EPEC[:ch_st][t] * EPEC[:λ][t] for t in t_steps)
            #NOTE: THIS IS THE LINEARIZED VERSION
            # obj =
            #     - (- sum(EPEC[:g_st][t] * dat_str.GEN_STR.MC_g[t] for t in t_steps) #GEN_STR
            #     - sum(EPEC[:w_st][t] * dat_str.WIND_STR.MC_g[t] for t in t_steps) #WIND_STR
            #     - sum(sum(EPEC[:g][i,t] * dat_nonstr[i].c_g[t] for t in t_steps) for i in non_str_gens) #GENs
            #     - sum(sum(EPEC[:w][k,t] * dat_nonstr[k].c_g[t] for t in t_steps) for k in non_str_res) #GENs
            #     + sum(EPEC[:d][t] .* p_cap for t in t_steps)
            #     - sum(EPEC[:β⁺g_st][t] * dat_str.GEN_STR.P_max[t] for t in t_steps) #GEN_STR
            #     - sum(EPEC[:β⁺w_st][t] * dat_str.WIND_STR.P_max[t] for t in t_steps) #WIND_STR
            #     - sum(sum(EPEC[:β⁺g][i,t] * dat_nonstr[i].G[t] for t in t_steps) for i in non_str_gens) #GENs
            #     - sum(sum(EPEC[:β⁺w][k,t] * dat_nonstr[k].G[t] for t in t_steps) for k in non_str_res) #GENs
            #     - sum(EPEC[:β⁺d][t] * load[t] for t in t_steps)
            # )

        elseif SP_obj == "GEN_fav"
            #FAVORING THE CONVENTIONAL GENERATOR
            obj = (
                #GEN
                - (
                    - sum(EPEC[:w_st][t] * dat_str.WIND_STR.MC_g[t] for t in t_steps) #WIND_STR
                    - sum(sum(EPEC[:g][i,t] * dat_nonstr[i].c_g[t] for t in t_steps) for i in non_str_gens) #GENs
                    - sum(sum(EPEC[:w][k,t] * dat_nonstr[k].c_g[t] for t in t_steps) for k in non_str_res) #GENs
                    + sum(EPEC[:d][t] .* p_cap for t in t_steps)
                    - sum(EPEC[:β⁺dch_st][t] * dat_str.ESS_STR.P_max[t] for t in t_steps) #WIND_STR
                    - sum(EPEC[:β⁺ch_st][t] * dat_str.ESS_STR.P_max[t] for t in t_steps) #WIND_STR
                    - sum(EPEC[:β⁺w_st][t] * dat_str.WIND_STR.P_max[t] for t in t_steps) #WIND_STR
                    - sum(sum(EPEC[:β⁺g][i,t] * dat_nonstr[i].G[t] for t in t_steps) for i in non_str_gens) #GENs
                    - sum(sum(EPEC[:β⁺w][k,t] * dat_nonstr[k].G[t] for t in t_steps) for k in non_str_res) #GENs
                    - sum(EPEC[:β⁺d][t] * load[t] for t in t_steps)
                    - sum(EPEC[:κ⁺soc][t] * dat_str.ESS_STR.E_max[t] for t in t_steps)
                )
            )
        else
            println("ERROR:The chosen Social Planner objective is not part of the list!!!")
        end

        return obj
    end

    println("...Setting social planner's objective: $SP_obj...")
    obj_EPEC = get_EPEC_obj()
    JuMP.set_objective_function(EPEC, obj_EPEC)
    JuMP.set_objective_sense(EPEC, MOI.MIN_SENSE)

    ##########SETTING THE OPTIMIZER#########
    set_optimizer(EPEC, optimizer)

    return EPEC
end


text:
# SKCE

# predicted Laplace distributions with exponential kernel for the targets
function CalibrationErrors.unsafe_skce_eval_targets(
    kernel::ExponentialKernel{Euclidean}, p::Laplace, y::Real, p̃::Laplace, ỹ::Real
)
    # extract the parameters
    μ = p.μ
    β = p.θ
    μ̃ = p̃.μ
    β̃ = p̃.θ

    res =
        kernel(y, ỹ) - laplace_laplacian_kernel(β, abs(μ - ỹ)) -
        laplace_laplacian_kernel(β̃, abs(μ̃ - y)) +
        laplace_laplacian_kernel(β, β̃, abs(μ - μ̃))

    return res
end

# z = abs(μ - y)
function laplace_laplacian_kernel(β, z)
    if isone(β)
        (1 + z) * exp(-z) / 2
    else
        (β * exp(-z / β) - exp(-z)) / (β^2 - 1)
    end
end

# z = abs(μ - μ̃)
function laplace_laplacian_kernel(β, β̃, z)
    if isone(β)
        if isone(β̃)
            return (3 + 3 * z + z^2) * exp(-z) / 8
        else
            c = β̃^2 - 1
            csq = c^2
            return β̃^3 * exp(-z / β̃) / csq - ((1 + z) / (2 * c) + β̃^2 / csq) * exp(-z)
        end
    end

    if isone(β̃)
        c = β^2 - 1
        csq = c^2
        return β^3 * exp(-z / β) / csq - ((1 + z) / (2 * c) + β^2 / csq) * exp(-z)
    elseif β̃ == β
        c = β^2 - 1
        csq = c^2
        return exp(-z) / csq + ((β + z) / (2 * c) - β / csq) * exp(-z / β)
    else
        c1 = β^2 - 1
        c2 = β̃^2 - 1
        c3 = β^2 - β̃^2
        return β^3 * exp(-z / β) / (c1 * c3) - β̃^3 * exp(-z / β̃) / (c2 * c3) +
               exp(-z) / (c1 * c2)
    end
end

# UCME

function CalibrationErrors.unsafe_ucme_eval_targets(
    kernel::ExponentialKernel{Euclidean}, p::Laplace, y::Real, ::Laplace, testy::Real
)
    return kernel(y, testy) - laplace_laplacian_kernel(p.θ, abs(p.μ - testy))
end

# kernels with input transformations
# TODO: scale upfront?
function CalibrationErrors.unsafe_skce_eval_targets(
    kernel::TransformedKernel{
        ExponentialKernel{Euclidean},<:Union{ScaleTransform,ARDTransform}
    },
    p::Laplace,
    y::Real,
    p̃::Laplace,
    ỹ::Real,
)
    # obtain the transform
    t = kernel.transform

    return CalibrationErrors.unsafe_skce_eval_targets(
        ExponentialKernel(), apply(t, p), t(y), apply(t, p̃), t(ỹ)
    )
end

function CalibrationErrors.unsafe_ucme_eval_targets(
    kernel::TransformedKernel{
        ExponentialKernel{Euclidean},<:Union{ScaleTransform,ARDTransform}
    },
    p::Laplace,
    y::Real,
    testp::Laplace,
    testy::Real,
)
    # obtain the transform
    t = kernel.transform

    # `testp` is irrelevant for the evaluation and therefore not transformed
    return CalibrationErrors.unsafe_ucme_eval_targets(
        ExponentialKernel(), apply(t, p), t(y), testp, t(testy)
    )
end

# utilities

# internal `apply` avoids type piracy of transforms
apply(t::Union{ScaleTransform,ARDTransform}, d::Laplace) = Laplace(t(d.μ), t(d.θ))


text:
module Day5

export day5

function findposition(str, upper)
    reduce(str, init = 0) do acc c
        c == upper ? acc+1 : acc
    end
end

function day5(::Val{:a}, filename)
    highest = 0
    for line in readlines(filename)
        line = strip(line)
        row_position = line[1:7]
        seat_position = line[8:end]
        row = findposition(row_position, 'B')
        seat = findposition(seat_position, 'R')
        seatid = row * 8 + seat
        if seatid > highest
            highest = seatid
        end
    end
    highest
end

function day5(::Val{:b}, filename)
    seats = Set{Int}()
    for line in readlines(filename)
        line = strip(line)
        row_position = line[1:7]
        seat_position = line[8:end]
        row = findposition(row_position, 'B')
        col = findposition(seat_position, 'R')
        push!(seats, row * 8 + col)
    end
    miss = setdiff(Set{Int}(0:127*8-1), seats)
    for seat in miss
        if in(seat-1, seats) && in(seat+1, seats)
            println(seat)
        end
    end
end

end


text:
using Distributions
using LinearAlgebra
using LinearManifoldModels
using Statistics
import StatsBase
import KernelDensity
import LMCLUS
import Random
using Test
# using Logging; global_logger(ConsoleLogger(stderr, LMCLUS.Base.CoreLogging.Debug))

const SEED = 986163731

# generate dataset
Random.seed!(SEED)
X1 = rand(Normal(0., 0.7), 2, 100)
X2 = rand(MvNormal([2., 2.], diagm(0=>[0.02, 0.6])), 100)
X = [X1 X2]
Xn = X .- mean(X, dims=2)

# generate clusters
p = LMCLUS.Parameters(1)
p.random_seed = SEED
p.sampling_heuristic = 1
p.sampling_factor = 1.5
p.max_bin_portion = 0.2
p.best_bound = 1.0
p.basis_alignment = true
p.dim_adjustment = true
p.dim_adjustment_ratio = 0.1
res = LMCLUS.lmclus(Xn, p)

# test distribution properties
@testset for E in [StatsBase.Histogram, KernelDensity.UnivariateKDE]
    for C in LMCLUS.manifolds(res)
        L = LMCLUS.points(C)
        lmcm = estimate(NonparametricEstimator(EmpiricalLinearManifold, E), C, Xn[:,L])
        @test length(lmcm) == 2
        @test params(lmcm) == (mean(C), LMCLUS.projection(C))
        @test partype(lmcm) == typeof(Xn[1])
        @test all(pdf(lmcm, Xn[:,L]) .> 0.)
        @test all(.!isinf.(logpdf(lmcm, Xn[:,L])))
        @test entropy(lmcm) > 0.
    end
end


text:
module Atp

export atp_ranking_get, atp_tourn_get

using Dates
import HTTP
using Gumbo
using Cascadia
using DataFrames

atprankingstart = Date(1973, 8, 23)

function http_get(url)
    rtn = nothing
    try
        rtn = HTTP.get(url)
    catch
        rtn = HTTP.get(url; require_ssl_verification = false)
    end
    return rtn
end

function atp_ranking_get(date_wanted::Date; race=false)
    if date_wanted < atprankingstart
        @warn "$date_wanted < $atprankingstart : " *
          "ranking for atprankingstart will be retrieved"
    end
    url_atp = ""
    sel_player = Selector(".player-cell")
    sel_points = Selector(".points-cell")
    #ranking date
    dropdown_class = "dropdown-holder-wrapper"
    select_dropdown = Selector(".$dropdown_class")
    if race
        url_atp = "https://www.atptour.com/en/rankings/singles-race-to-turin?"  # FIXME to-turin might depends on date
    else
        url_atp = "https://www.atptour.com/en/rankings/singles?"
    end
    http_atp = http_get(url_atp);
    html_str_atp = String(http_atp.body);
    html_atp = parsehtml(html_str_atp);
    dropdown_lists = eachmatch(select_dropdown, html_atp.root)
    dates_dropdown = dropdown_lists[3]
    n_dates = size(Gumbo.children(dates_dropdown[1][2]), 1)
    if size(eachmatch(sel_player, html_atp.root), 1) == 0
        date_wanted -= Day(7)
    end
    ranking_dates = Array{Date}(undef, n_dates)
    for date_i in 1:n_dates
        date_str = replace(text(dates_dropdown[1][2][date_i][1]),
                            r"\r|\t|\s"=>"")
        ranking_dates[date_i] = Date(date_str, dateformat"y.m.d")
    end
    ranking_date_idx = findfirst(date_wanted .- ranking_dates .>= Day(0))
    if isnothing(ranking_date_idx)
        ranking_date = ranking_dates[end]
    else
        ranking_date = ranking_dates[ranking_date_idx]
    end
    url_atp *= "rankDate=" * string(ranking_date) * "&rankRange=1-500"
    http_atp = http_get(url_atp);
    html_str_atp = String(http_atp.body);
    table_start = findfirst("<div id=\"singlesRanking", html_str_atp)[1];
    table_end = findnext("
", html_str_atp, table_start)[end];
    table_str = html_str_atp[table_start:table_end];
    table_html = parsehtml(table_str);
    player_table_html = eachmatch(sel_player, table_html.root)
    points_table_html = eachmatch(sel_points, table_html.root)
    rank_df = DataFrame(player=String[], points=Int[])
    for player_i in 1:size(player_table_html, 1)
        if player_i == 1 &&
          occursin(r"[A-Za-z]+", text(points_table_html[player_i]))
            continue
        end
        player = text(player_table_html[player_i][1][1])
        player = replace(player,  r"^[\n\s]*"=>"")
        points_str = text(points_table_html[player_i][1][1])
        points_str = replace(points_str, "," => "")
        points = parse(Int, points_str)
        push!(rank_df, (player, points))
    end
    return rank_df
end

function atp_tourn_get()
    url_atp = "https://www.atptour.com/en/tournaments"
    url_atp = url_atp
    http_atp = http_get(url_atp);
    html_str_atp = String(http_atp.body);
    html_atp = parsehtml(html_str_atp);
    select_tourn = Selector(".tourney-result")
    tourn_list = eachmatch(select_tourn, html_atp.root)
    tourn_GS_img = "/assets/atpwt/images/tournament/badges/categorystamps_grandslam.png"
    tourn_M1000_img = "/assets/atpwt/images/tournament/badges/categorystamps_1000.png"
    tourn_ATPF_img = "/assets/atpwt/images/tournament/badges/categorystamps_finals.svg"
    tourn_df = DataFrame(name=String[], cat=String[], date0=Date[],
                            date1=Date[], url=String[])
    for i_tourn in 1:size(tourn_list, 1)
        select_cat = Selector(".tourney-badge-wrapper")
        select_name = Selector(".tourney-title")
        select_loc = Selector(".tourney-location")
        select_dates = Selector(".tourney-dates")
        tourn_name = eachmatch(select_name, tourn_list[i_tourn])
        tourn_loc = eachmatch(select_loc, tourn_list[i_tourn])
        tourn_loc = text(tourn_loc[1][1])
        tourn_loc = replace(tourn_loc, r"\r|\t|\n"=>"")
        idx0 = findfirst(!isequal(' '), tourn_loc)
        idx1 = findlast(!isequal(' '), tourn_loc)
        tourn_loc = tourn_loc[idx0:idx1]
        tourn_cat = eachmatch(select_cat, tourn_list[i_tourn])
        tourn_dates = eachmatch(select_dates, tourn_list[i_tourn])
        name = text(tourn_name[1][1]) * " (" * tourn_loc * ")"
        name = replace(replace(name, r"\n[\s]*"=>" "), r"^[\s]*"=>"")
        if occursin("Olympic", name)
            tourn_cat_img = "OG"
        else
            tourn_cat_img = Gumbo.attrs(tourn_cat[1][1])["src"]
        end
        if !haskey(Gumbo.attrs(tourn_name[1]), "href")
            continue
        end
        url = Gumbo.attrs(tourn_name[1])["href"]
        url = url
        dates = text(tourn_dates[1][1])
        dates = replace(dates, r"\r|\t|\s"=>"")
        sep_idx = findfirst("-", dates)[1]
        date0 = Date(dates[1:sep_idx-1], dateformat"y.m.d")
        date1 = Date(dates[sep_idx+1:end], dateformat"y.m.d")
        if occursin(tourn_GS_img, tourn_cat_img)
            cat = "GS"
        elseif occursin(tourn_M1000_img, tourn_cat_img)
            cat = "M1000"
        elseif occursin("OG", tourn_cat_img)
            cat = "OG"
        elseif occursin(tourn_ATPF_img, tourn_cat_img)
            cat = "ATPF"
        else
            cat = "other"
        end
        cv_tourn = ["Indian Wells", "Miami", "Monte Carlo", "Madrid"]
        corona_virus = (year(date0) == 2020) && any(occursin.(cv_tourn, name))
        suspended = occursin("Suspended", name)
        postponed = occursin("Postponed", name)
        if cat != "other" && !corona_virus && !suspended && !postponed
            push!(tourn_df, Dict(:name=>name, :cat=>cat, :date0=>date0,
                                    :date1=>date1, :url=>url))
        end
    end
    return tourn_df
end

end


text:
module test_octo_functions

using Test
using Octo.Adapters.SQL

struct T
end

t = from(T)
@test to_sql([AVG(t.salary)]) == "AVG(salary)"

s = as(AVG(t.salary), :s)
@test to_sql([s]) == "AVG(salary) AS s"

using Octo: @sql_functions

@sql_functions A

@test to_sql([A(2,3)]) == "A(2, 3)"

@test to_sql([NOT(t.salary > 100)]) == "NOT(salary > 100)"
@test to_sql([SOME(t.salary > 100)]) == "SOME(salary > 100)"
@test to_sql([EVERY(t.salary > 100)]) == "EVERY(salary > 100)"

end # module test_octo_functions


text:
using Documenter, Luxor

makedocs(
  modules = [Luxor],
  format = Documenter.Formats.HTML,
  sitename = "Luxor",
  pages    = Any[
    "Introduction to Luxor"           => "index.md",
    "A few examples"                  => "examples.md",
    "Basic graphics"                  => "basics.md",
    "Transforms and matrices"         => "transforms.md",
    "Clipping"                        => "clipping.md",
    "Images"                          => "images.md",
    "Turtle graphics"                 => "turtle.md",
    "More examples"                   => "moreexamples.md",
    "Index"                           => "functionindex.md"
    ]
  )

#=
after this has built the source files, some of them
need to be moved up from build/ into docs/ so
that they're served at http://cormullion.github.io/Luxor.jl/:

docs
  build
    assets                ->
      documenter.css      -> not moved, we've customized it
      documenter.js       -> moved to docs/assets/documenter
      search.js           -> moved to docs/assets/search.js
      style.css           -> not moved, we've customized it
    index.html            -> moved to docs/index.html
    search_index.js       -> moved to docs/index.html
    search.html           -> moved to docs/index.html

=#

info("moving files")

#mv("build/assets/documenter.js", "assets/documenter.js", remove_destination=true)
#mv("build/assets/search.js", "assets/search.js", remove_destination=true)
#mv("build/assets/logo.png", "assets/logo.png", remove_destination=true)

mv("build/search_index.js", "search_index.js", remove_destination=true)

for f in filter(f -> endswith(f, ".html"), readdir("build/"))
  info("moving $f")
  mv(string("build/", f), f, remove_destination=true)
end


text:
module DelayDiffEq

using Reexport
@reexport using OrdinaryDiffEq

using DataStructures
using LinearAlgebra
using Logging
using Printf
using RecursiveArrayTools
using UnPack

import NonlinearSolve

using DiffEqBase: AbstractDDEAlgorithm, AbstractDDEIntegrator, AbstractODEIntegrator, DEIntegrator, AbstractDDEProblem

using DiffEqBase: @..

using DiffEqBase: FastConvergence, Convergence, SlowConvergence, VerySlowConvergence, Divergence

using OrdinaryDiffEq: RosenbrockMutableCache

export Discontinuity, MethodOfSteps

include("discontinuity_type.jl")
include("functionwrapper.jl")
include("ddestats.jl")

include("integrators/type.jl")
include("integrators/utils.jl")
include("integrators/interface.jl")

include("fpsolve/type.jl")
include("fpsolve/fpsolve.jl")
include("fpsolve/utils.jl")
include("fpsolve/functional.jl")
include("cache_utils.jl")
include("interpolants.jl")
include("history_function.jl")
include("algorithms.jl")
include("track.jl")
include("alg_utils.jl")
include("solve.jl")
include("utils.jl")

end # module


text:
module h20

  using FinEtools
  using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
  using FinEtoolsDeforLinear: FEMMDeforLinear
  using FinEtoolsDeforNonlinear
  using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
  using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
  using FinEtoolsDeforNonlinear.AlgoDeforNonlinearModule: nonlinearstatics
  using LinearAlgebra: norm
  using Statistics: mean
  using SparseArrays
  using DelimitedFiles
  using Interpolations
  using UnicodePlots

  function neohookean_h20()
      mr = DeforModelRed3D
      E, nu = 7.0*phun("MPa"), 0.3
      m = MatDeforNeohookean(mr, E, nu)
      L= 6/2*phun("mm");
      H = 2/2*phun("mm");
      W = 2/2*phun("mm");
      tmag = 0.35*phun("MPa");# Magnitude of the traction
      nincr = 8
      utol = 10e-7;
      graphics = ~true;
      maxdu_tol = W/1e7;
      maxiter = 9
      tolerance = W / 1000
      traction_vector = [0.0, 0.0, -tmag]

      fens, fes = H20block(L, W, H, 60, 3, 3)

      l1  = selectnode(fens; box = [0,0,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
      e1 = FDataDict("node_list"=>l1, "component"=>1)
      e2 = FDataDict("node_list"=>l1, "component"=>2)
      e3 = FDataDict("node_list"=>l1, "component"=>3)

      bfes = meshboundary(fes)
      el1 = selectelem(fens, bfes, box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
      setvector!(v, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0) = begin
          v .= time .* traction_vector
          return v
      end
      fi = ForceIntensity(FFlt, length(traction_vector), setvector!, 0.0);
      t1 = FDataDict("femm"=>FEMMBase(IntegDomain(subset(bfes, el1), GaussRule(2, 2))), "traction_vector"=>fi)

      movel1  = selectnode(fens; box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)

      femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), m)

      region1 = FDataDict("femm"=>femm)
      modeldata =  FDataDict("fens"=> fens, "regions"=>  [region1], "traction_bcs"=>[t1], "essential_bcs"=>[e1, e2, e3])

      modeldata["load_multipliers"] = (1:nincr)./nincr*1.0;
      modeldata["maxdu_tol"] = maxdu_tol;
      modeldata["maxiter"] = maxiter;
      modeldata["line_search"]  = true;
      modeldata["iteration_observer"] = (lambda, iter, du, modeldata) -> begin
          @show lambda, iter, modeldata["maxdu"], modeldata["maxbal"]
      end
      Ux = FFlt[]; lambdas = FFlt[]
      modeldata["increment_observer"] = (lambda, modeldata) -> begin
          push!(Ux, mean(modeldata["un1"].values[movel1,3]));
          push!(lambdas, lambda);
      end

      modeldata = nonlinearstatics(modeldata);
      @show Ux / phun("mm"), lambdas

      pl = lineplot(Ux / phun("mm"), lambdas)
      display(pl)

      vtkexportmesh("neohookean_h20.vtk", fens, fes; vectors = [("u", modeldata["un1"].values)])
  end # function neohookean_h8

end #


text:
import TransportBasedInference: optimize

@testset "Verify clearcoeff!" begin
  m = 5
  Nx = 10
  Nψ = 5

  coeff = randn(Nψ)

  C = HermiteMapComponent[]

  for i=1:Nx
      idxi = rand(0:m, Nψ, i)
      coeffi = randn(Nψ)
      push!(C, HermiteMapComponent(m, i, deepcopy(idxi), deepcopy(coeffi)))
  end

  M = HermiteMap(m, Nx, LinearTransform(Nx), C)

  clearcoeff!(M)

  for i=1:Nx
      @test norm(getcoeff(M.C[i]) - zeros(Nψ))<1e-12
  end
end

@testset "Test evaluation of HermiteMap" begin

    Nx = 3
    Ne = 50
    m = 10
    X = randn(Nx, Ne) .* randn(Nx, Ne) .+ rand(Nx);
    X0 = deepcopy(X)

    B = CstProHermiteBasis(m-2)

    idx1   = reshape([0; 1; 2], (3,1))
    coeff1 = randn(3)
    MB1 = MultiBasis(B, 1)
    H1 = HermiteMapComponent(IntegratedFunction(ExpandedFunction(MB1, idx1, coeff1)))

    idx2   = reshape([0 0 ;0 1; 1 0; 2 0], (4,2))
    coeff2 = randn(4)
    MB2 = MultiBasis(B, 2)
    H2 = HermiteMapComponent(IntegratedFunction(ExpandedFunction(MB2, idx2, coeff2)))

    idx3   = reshape([0 0 0; 0 0 1; 0 1 0; 1 0 1], (4,3))
    coeff3 = randn(4)
    MB3 = MultiBasis(B, 3)
    H3 = HermiteMapComponent(IntegratedFunction(ExpandedFunction(MB3, idx3, coeff3)))

    L = LinearTransform(X; diag = true)
    M = HermiteMap(L, HermiteMapComponent[H1; H2; H3])

    out = zero(X)
    # Without rescaling of the variables

    evaluate!(out, M, X; apply_rescaling = false)

    @test norm(out[1,:] - evaluate(H1, X[1:1,:]))<1e-10
    @test norm(out[2,:] - evaluate(H2, X[1:2,:]))<1e-10
    @test norm(out[3,:] - evaluate(H3, X[1:3,:]))<1e-10

    # With rescaling of the variables

    evaluate!(out, M, X; apply_rescaling = true)

    transform!(L, X)
    @test norm(out[1,:] - evaluate(H1, X[1:1,:]))<1e-10
    @test norm(out[2,:] - evaluate(H2, X[1:2,:]))<1e-10
    @test norm(out[3,:] - evaluate(H3, X[1:3,:]))<1e-10

    itransform!(L, X)

    @test norm(X0 - X)<1e-12
end

@testset "Validate evalution of HermiteMap with multi-threading" begin

    Nx = 3
    Ne = 50
    m = 10
    X = randn(Nx, Ne) .* randn(Nx, Ne) .+ rand(Nx);
    X0 = deepcopy(X)

    B = CstProHermiteBasis(m-2)

    idx1   = reshape([0; 1; 2], (3,1))
    coeff1 = randn(3)
    MB1 = MultiBasis(B, 1)
    H1 = HermiteMapComponent(IntegratedFunction(ExpandedFunction(MB1, idx1, coeff1)))

    idx2   = reshape([0 0 ;0 1; 1 0; 2 0], (4,2))
    coeff2 = randn(4)
    MB2 = MultiBasis(B, 2)
    H2 = HermiteMapComponent(IntegratedFunction(ExpandedFunction(MB2, idx2, coeff2)))

    idx3   = reshape([0 0 0; 0 0 1; 0 1 0; 1 0 1], (4,3))
    coeff3 = randn(4)
    MB3 = MultiBasis(B, 3)
    H3 = HermiteMapComponent(IntegratedFunction(ExpandedFunction(MB3, idx3, coeff3)))

    L = LinearTransform(X; diag = true)
    M = HermiteMap(L, HermiteMapComponent[H1; H2; H3])

    out = zero(X)
    out_thread = zero(X)
    # Without rescaling of the variables

    evaluate!(out, M, X; apply_rescaling = false, P = serial)
    evaluate!(out_thread, M, X; apply_rescaling = false, P = thread)

    @test norm(out[1,:] - evaluate(H1, X[1:1,:]))<1e-10
    @test norm(out[2,:] - evaluate(H2, X[1:2,:]))<1e-10
    @test norm(out[3,:] - evaluate(H3, X[1:3,:]))<1e-10

    @test norm(out_thread[1,:] - evaluate(H1, X[1:1,:]))<1e-10
    @test norm(out_thread[2,:] - evaluate(H2, X[1:2,:]))<1e-10
    @test norm(out_thread[3,:] - evaluate(H3, X[1:3,:]))<1e-10

    # With rescaling of the variables

    evaluate!(out, M, X; apply_rescaling = true, P = serial)
    evaluate!(out_thread, M, X; apply_rescaling = true, P = thread)

    transform!(L, X)
    @test norm(out[1,:] - evaluate(H1, X[1:1,:]))<1e-10
    @test norm(out[2,:] - evaluate(H2, X[1:2,:]))<1e-10
    @test norm(out[3,:] - evaluate(H3, X[1:3,:]))<1e-10

    @test norm(out_thread[1,:] - evaluate(H1, X[1:1,:]))<1e-10
    @test norm(out_thread[2,:] - evaluate(H2, X[1:2,:]))<1e-10
    @test norm(out_thread[3,:] - evaluate(H3, X[1:3,:]))<1e-10

    itransform!(L, X)

    @test norm(X0 - X)<1e-12
end

@testset "Verify that Serial and Multi-threading optimization are working properly without/with QR" begin
    Nx = 20
    Ne = 300
    m = 10

    X = randn(Nx, Ne) .* randn(Nx, Ne) .+ randn(Nx)
    M_noqr = HermiteMap(m, X)
    Mthread_noqr = HermiteMap(m,X)

    M_qr = HermiteMap(m, X)
    Mthread_qr = HermiteMap(m,X)

    optimize(M_noqr, X, 10; withqr = false, P = serial)
    optimize(Mthread_noqr, X, 10; withqr = false, P = thread)

    optimize(M_qr, X, 10; withqr = true, P = serial)
    optimize(Mthread_qr, X, 10; withqr = true, P = thread)

    for i=1:Nx
        @test norm(getcoeff(M_noqr.C[i]) - getcoeff(Mthread_noqr.C[i]))<1e-4
        @test norm(getcoeff(M_noqr.C[i]) - getcoeff(M_qr.C[i]))<1e-4
        @test norm(getcoeff(M_noqr.C[i]) - getcoeff(Mthread_qr.C[i]))<1e-4
    end

    @test norm(evaluate(M_noqr, X; P = serial) - evaluate(Mthread_noqr, X; P = thread))<1e-4
    @test norm(evaluate(M_noqr, X; P = serial) - evaluate(M_qr, X; P = thread))<1e-4
    @test norm(evaluate(M_noqr, X; P = serial) - evaluate(Mthread_qr, X; P = thread))<1e-4

end

@testset "Verify log_pdf function" begin
    # For diagonal rescaling
    Nx = 2
    m = 10
    X  =  Matrix([0.267333   1.43021;
              0.364979   0.607224;
             -1.23693    0.249277;
             -2.0526     0.915629;
             -0.182465   0.415874;
              0.412907   1.01672;
              1.41332   -0.918205;
              0.766647  -1.00445]');
    X0 = deepcopy(X)

    B = MultiBasis(CstProHermiteBasis(m), Nx);

    M = HermiteMap(m, X);

    idx1 = reshape([0; 1; 2;3],(4,1))
    coeff1 =[-7.37835916713439;
             16.897974063168977;
             -0.20935914261414929;
              8.116032140446544];

    M[1] = HermiteMapComponent(m, 1,  idx1, coeff1)

    idx2 = [ 0  0; 0 1; 1 0; 2 0];
    coeff2 = [2.995732886294909; -2.2624558703623903; -3.3791974895855486; 1.3808989978516617]

    M[2] = HermiteMapComponent(m, 2,  idx2, coeff2);

    @test norm(log_pdf(M, X) -   [-1.849644462509894;
                                  -1.239419791201162;
                                  -3.307440928766559;
                                  -2.482615159224171;
                                  -2.392017816948268;
                                  -1.331903266052277;
                                  -3.638709417555485;
                                  -2.242238529589172])<1e-6

    @test norm(X - X0)<1e-10

    # For Cholesky refactorization

    Nx = 2
    m = 10
    X  =  Matrix([0.267333   1.43021;
              0.364979   0.607224;
             -1.23693    0.249277;
             -2.0526     0.915629;
             -0.182465   0.415874;
              0.412907   1.01672;
              1.41332   -0.918205;
              0.766647  -1.00445]');
    X0 = deepcopy(X)

    B = MultiBasis(CstProHermiteBasis(m), Nx);

    M = HermiteMap(m, X; diag = false);

    idx1 = reshape([0; 1; 2;3],(4,1))
    coeff1 =[-7.37835916713439;
             16.897974063168977;
             -0.20935914261414929;
              8.116032140446544];

    M[1] = HermiteMapComponent(m, 1,  idx1, coeff1)

    idx2 = [ 0  0; 0 1; 1 0; 2 0];
    coeff2 = [2.995732886294909; -2.2624558703623903; -3.3791974895855486; 1.3808989978516617]

    M[2] = HermiteMapComponent(m, 2,  idx2, coeff2);

    @test norm(log_pdf(M, X) - [  -2.284420690193831;
                                  -1.113693643793066;
                                  -3.927922581187332;
                                  -2.461816757918692;
                                  -2.328286973844927;
                                  -1.529996011165051;
                                  -3.444940697346834;
                                  -1.999419335266747])<1e-6

    @test norm(X - X0)<1e-10

end

@testset "Verify the different optimization max_terms ∈ Int64" begin

    Nx = 3
    Ne = 8

    X = Matrix([  0.112122    0.371145  -1.09587;
                 -0.124971    0.563166   1.21132;
                 -0.0781768   0.974237   1.16974;
                  0.7536      3.67691    1.31592;
                  0.329414   -1.57155    1.02031;
                  0.0263465   0.445467   1.33992;
                  0.457334    0.20462    1.91469;
                  0.190431   -0.22547    1.25082]');

    m = 30

    M = HermiteMap(m, X)

    optimize(M, X, 5; withconstant = false)

    @test norm(log_pdf(M, X) - [  -1.648733127003676;
                                  -0.524080373136291;
                                   0.466815189258083;
                                  -1.265418705621051;
                                  -4.178674562888346;
                                   0.977833056870988;
                                  -3.449875825411972;
                                   0.460274224407305])<1e-3

    @test getidx(M[1]) == reshape([1; 2; 3; 4; 5],(5,1))
    @test norm(getcoeff(M[1]) - [-1.485021686942572;
                                  1.741440183812295;
                                 -3.814293233298080;
                                  2.150582375812929;
                                 -1.418659400598183])<1e-3

    @test getidx(M[2]) == [1 0; 0 1; 0 2; 0 3; 2 0]
    @test norm(getcoeff(M[2]) - [18.385719748659056 ;
                                -20.079305399569837 ;
                                 5.208177303845926  ;
                                -7.543285947028642  ;
                                 1.043264254165708])<1e-3

    @test getidx(M[3]) == [0 0 1; 0 0 2; 0 0 3; 0 0 4; 0 1 0]
    @test norm(getcoeff(M[3]) - [  9.801288179417845;
                                 -18.297926511077648;
                                  14.512139032442953;
                                 -16.509726071255606;
                                  -0.962812577491049])<1e-2

end

@testset "Verify the different optimization max_terms = Split" begin

    Nx = 3
    Ne = 8

    X = Matrix([  0.112122    0.371145  -1.09587;
                 -0.124971    0.563166   1.21132;
                 -0.0781768   0.974237   1.16974;
                  0.7536      3.67691    1.31592;
                  0.329414   -1.57155    1.02031;
                  0.0263465   0.445467   1.33992;
                  0.457334    0.20462    1.91469;
                  0.190431   -0.22547    1.25082]');

    m = 30

    M = HermiteMap(m, X)

    optimize(M, X, "Split"; withconstant = false)

    @test norm(log_pdf(M, X) - [ -21.745907202610731;
                                  -0.464429325572978;
                                  -0.155330841509541;
                                  -4.329981710326985;
                                  -2.539350574125156;
                                   0.206049495114118;
                                  -3.760956579928015;
                                  -0.463356887488941])<1e-4

    @test getidx(M[1]) == reshape([1; 2; 3],(3,1))
    @test norm(getcoeff(M[1]) - [-0.907122832552077;
                                 -0.471583533261757;
                                 -1.653849103515142])<1e-4

    @test getidx(M[2]) == [1 0; 0 1; 0 2]
    @test norm(getcoeff(M[2]) - [ 6.254057966900727;
                                 -5.202600223050682;
                                  1.829358246492261])<1e-4

    @test getidx(M[3]) == [0 0 1; 0 0 2; 0 1 0]
    @test norm(getcoeff(M[3]) - [ -0.348002361230585;
                                   6.995417803422829;
                                  -2.169327382026170])<1e-4

end

@testset "Verify the different optimization max_terms = kfold" begin

    Nx = 3
    Ne = 16

    X = Matrix([  0.112122    0.371145  -1.09587  ;
                 -0.124971    0.563166   1.21132  ;
                 -0.0781768   0.974237   1.16974  ;
                  0.7536      3.67691    1.31592  ;
                  0.329414   -1.57155    1.02031  ;
                  0.0263465   0.445467   1.33992  ;
                  0.457334    0.20462    1.91469  ;
                  0.190431   -0.22547    1.25082  ;
                 -0.383296    0.838512  -0.687386 ;
                 -1.57094    -0.808855   0.0632401;
                 -0.684854   -0.12027    0.187772 ;
                  0.767601    2.35813   -0.658357 ;
                  0.338965   -0.168985  -0.00250368;
                 -0.27756    -0.790416  -0.133158 ;
                 -0.824008    0.459429  -0.235808 ;
                  0.0145145   0.123892   0.0831142]')

    m = 30

    M = HermiteMap(m, X)

    M = HermiteMap(m, X)

    optimize(M, X, "kfold"; withconstant = false)

    @test norm(log_pdf(M, X) - [  -3.880921471282678;
                                  -2.753134673606892;
                                  -2.839021221004361;
                                  -7.308105457845060;
                                  -3.838700341109559;
                                  -2.845755023406152;
                                  -3.982132693868225;
                                  -2.869826713565865;
                                  -3.474164717698875;
                                  -6.086150330191086;
                                  -3.064117775500836;
                                  -5.349287829523062;
                                  -2.704253030752519;
                                  -3.037988603791519;
                                  -3.557543327559524;
                                  -2.443326419210569])<1e-4

    @test getidx(M[1]) == reshape([1],(1,1))
    @test norm(getcoeff(M[1]) - [-0.118440805265421])<1e-4

    @test getidx(M[2]) == [0 1]
    @test norm(getcoeff(M[2]) - [0.146373981523974])<1e-4

    @test getidx(M[3]) == [1 0 0]
    @test norm(getcoeff(M[3]) - [-0.067872411753720])<1e-4
end

@testset "Test inversion of the Hermite Map I" begin
    Nx = 100
    m = 20
    Ne = 500
    Xprior = randn(Nx, Ne).*randn(Nx, Ne)
    Xpost = deepcopy(Xprior) .+ randn(Nx, Ne)

    M = HermiteMap(m, Xprior; diag = true)
    optimize(M, Xprior, 5; withconstant = false)

    F = evaluate(M, Xprior)
    inverse!(Xpost, F, M)

    @test norm(Xprior - Xpost)/norm(Xpost)<1e-6
    @test norm(Xprior - Xpost)<1e-6
    @test norm(evaluate(M, Xprior) - evaluate(M, Xpost))<1e-6
end

@testset "Test inversion of the Hermite Map II" begin

    Nx = 100
    Ny = 50
    m = 20
    Ne = 500
    Xprior = randn(Nx, Ne).*randn(Nx, Ne)
    Xpost = deepcopy(Xprior) .+ randn(Nx, Ne)

    M = HermiteMap(m, Xprior; diag = true)
    Ystar = deepcopy(Xpost[1:Ny,:])# + 0.05*randn(Ny,Ne);

    M = HermiteMap(m, Xprior)
    optimize(M, Xprior, 5; withconstant = false, start = Ny+1)

    F = evaluate(M, Xpost; start = Ny+1)
    inverse!(Xprior, F, M, Ystar; start = Ny+1, P = serial)

    @test norm(evaluate(M, Xprior; start = Ny+1)-evaluate(M, Xpost; start = Ny+1))/norm(evaluate(M, Xpost; start = Ny+1))<1e-6

    @test norm(Xprior[Ny+1:end,:] - Xpost[Ny+1:end,:])/norm(Xpost[Ny+1:end,:])<1e-6
end

@testset "Test hybrid inversion of the Hermite Map I" begin
    Nx = 100
    m = 20
    Ne = 500
    Xprior = randn(Nx, Ne).*randn(Nx, Ne)
    Xpost = deepcopy(Xprior) .+ randn(Nx, Ne)

    M = HermiteMap(m, Xprior; diag = true)
    optimize(M, Xprior, 5; withconstant = false)

    F = evaluate(M, Xprior)
    hybridinverse!(Xpost, F, M)

    @test norm(Xprior - Xpost)/norm(Xpost)<1e-6
    @test norm(Xprior - Xpost)<1e-6
    @test norm(evaluate(M, Xprior) - evaluate(M, Xpost))<1e-6
end

@testset "Test hybrid inversion of the Hermite Map II" begin

    Nx = 100
    Ny = 50
    m = 20
    Ne = 500
    Xprior = randn(Nx, Ne).*randn(Nx, Ne)
    Xpost = deepcopy(Xprior) .+ randn(Nx, Ne)

    M = HermiteMap(m, Xprior; diag = true)
    Ystar = deepcopy(Xpost[1:Ny,:])# + 0.05*randn(Ny,Ne);

    M = HermiteMap(m, Xprior)
    optimize(M, Xprior, 5; withconstant = false, start = Ny+1)

    F = evaluate(M, Xpost; start = Ny+1)
    inverse!(Xprior, F, M, Ystar; start = Ny+1, P = serial)

    @test norm(evaluate(M, Xprior; start = Ny+1)-evaluate(M, Xpost; start = Ny+1))/norm(evaluate(M, Xpost; start = Ny+1))<1e-6

    @test norm(Xprior[Ny+1:end,:] - Xpost[Ny+1:end,:])/norm(Xpost[Ny+1:end,:])<1e-6
end


text:
include("../set_loadpath.jl")
using LinAlgUtils
using JacobianUtils
using JSON
using Rotations
using LinearAlgebra

###########
# OPTIONS #
###########
# bal_problem = "ladybug/problem-372-47423-post.txt"
# output_json = "vis/sample_scene.json"

bal_problem = "venice/problem-1490-935273-post.txt"
output_json = "vis/big_sample_scene.json"

# bal_problem = "trafalgar/problem-257-65132-post.txt"
# bal_problem = "dubrovnik/problem-253-163691-post.txt"
# output_json = "vis/other_sample_scene.json"
num_vecs = 20

##########
# SCRIPT #
##########

# extract the jacobian from a problem instance
problem_file = normpath(joinpath(@__DIR__, "..", "dataset", "bal", bal_problem))
compute_jacobian = normpath(joinpath(@__DIR__, "..", "build", "bin", "compute_jacobian"))
jacobian_file = tempname() * ".h5" # this is string concat in julia??? !!!
cmd = `$compute_jacobian "$problem_file" "$jacobian_file"`
run(cmd)

# get the eigenvectors from this jacobian
problem = readproblem(jacobian_file)
rm(jacobian_file)
λ, X = compute_covariance_evecs(problem, num_vecs, method=:dropping)

# TEMP: camera center calculations so we don't have to figure out how to do
# this in Javascript
camera_centers = zeros(Float64, 3*problem.num_cameras)
for k = 0:problem.num_cameras-1
    rot = problem.camera_parameters[9*k+1:9*k+3]
    θ = norm(rot)
    rot_aa = AngleAxis(θ, rot[1], rot[2], rot[3])
    t = problem.camera_parameters[9*k+4:9*k+6]
    camera_centers[3*k+1 : 3*k+3] = -inv(rot_aa) * t
end

# stuff this into a json file
json_data = Dict()
json_data["num_cameras"] = problem.num_cameras
json_data["num_points"] = problem.num_points
json_data["camera_parameters"] = problem.camera_parameters
json_data["point_parameters"] = problem.point_parameters
json_data["camera_centers"] = camera_centers

evecs = []
for i=1:size(X,2)
    push!(evecs, X[:,i])
end
json_data["eigenvectors"] = evecs
json_data["eigenvalues"] = λ

# write the json structure to a file
open(output_json, "w") do f
    write(f, "sample_scene = ")
    JSON.print(f, json_data)
end

max_stars_repo_name: wilsonkl/sfmflex-release

max_issues_repo_name: wilsonkl/sfmflex-release

max_forks_repo_name: wilsonkl/sfmflex-release


text:
"""This file initialize module robox."""
# Author: Henri Gérard <hgerard.proy@gmail.com>
# License: MIT

__precompile__()

module RobustOptimization

using JuMP, StatsBase
using Distributed
using DelimitedFiles
using LinearAlgebra

export DRO_opt, normal_opt, KL_opt, run_algo,
        error_loss, error_prime, equivalent,
        fconstraint, subgradient, typeof_constraint,
        x_in_inter, aux_proj, QR, proj_pi,
        GeneralConstraint, KLConstraint, DROConstraint, PositiveConstraint, EntropicConstraint,
        LinearRegression, LogisticRegression,
        Parallel, Sequential,
        OptParams, ProjParams,
        RobustModel,
        DRO_opt, normal_opt, KL_opt,
        read_data_libsvm, getsolution, pred, initialize, init_proj, create_data,
        dispatch_index, proj_pi, algo_proj_bis, algo_proj, run_algo_bis, positive_rate

include("objects.jl")
include("projection.jl")
include("ipsolve.jl")
include("preprocess.jl")
include("postprocess.jl")
include("constraint_functions.jl")
include("error_functions.jl")
include("utils.jl")

end

max_stars_repo_name: Henri-Gerard/robox

max_issues_repo_name: Henri-Gerard/robox

max_forks_repo_name: Henri-Gerard/robox


text:
using SymbolicUtils
import SymbolicUtils: Chain, Postwalk
import SymbolicUtils: Sym, Term, istree, operation, arguments, to_symbolic, Fixpoint

for unary_operator in unary_operators
    b_name, b_symbol = Meta.parse.(unary_operator)
    @eval SymbolicUtils.istree(a::$b_name) = true
    @eval SymbolicUtils.arguments(a::$b_name) = [a.term]
    @eval SymbolicUtils.operation(a::$b_name) = $b_symbol
    @eval SymbolicUtils.symtype(a::$b_name) = Number
end

# Define Struct and Symbol Overload for Binary Operators
for binary_operator in binary_operators
    b_name, b_symbol = Meta.parse.(binary_operator)
    @eval SymbolicUtils.istree(a::$b_name) = true
    @eval SymbolicUtils.arguments(a::$b_name) = [a.term1, a.term2]
    @eval SymbolicUtils.operation(a::$b_name) = $b_symbol
    @eval SymbolicUtils.symtype(a::$b_name) = Number
end

# Define Struct and Symbol Overload for n-ary Operators
for nary_operator in nary_operators
    b_name, b_symbol = Meta.parse.(nary_operator)
    @eval SymbolicUtils.istree(a::$b_name) = true
    @eval SymbolicUtils.arguments(a::$b_name) = [a.terms...]
    @eval SymbolicUtils.operation(a::$b_name) = $b_symbol
    @eval SymbolicUtils.symtype(a::$b_name) = Number
end