text:
include("./downloader.jl")
include("./utils.jl")
include("./menu.jl")

using Base.Threads
using Dates
using DataFrames
using CSV
using YAML
using ProgressMeter

function main()
    fill_value_regional, fill_value_point_daily, fill_value_point_hourly  = nothing, nothing, nothing
    parsed_args = Menu.main_menu()
    locations = YAML.load_file("config/locations.yml")

    # Prune the locations list
    for (location_name_A, point_A) in locations["target_locations"]
        for (location_name_B, point_B) in locations["target_locations"]
            if location_name_A != location_name_B
                region_A = Utils.get_area(point_A["location"], parsed_args["width"], parsed_args["height"])
                region_B = Utils.get_area(point_B["location"], parsed_args["width"], parsed_args["height"])
                iou = Utils.intersection_over_union(region_A, region_B)

                # Choose one of two locations to keep by it's power or area if defined.
                # The last option is random selection.
                if iou > 0.1
                    println("Point A: $(location_name_A)")
                    println("Point B: $(location_name_B)")
                    println("IoU: $(iou)")
                    println()

                    if (haskey(point_A, "permanent") && point_A["permanent"] == true)
                        if (haskey(point_B, "permanent") == false || point_B["permanent"] == false)
                            delete!(locations["target_locations"], location_name_B)
                        end
                    elseif (haskey(point_B, "permanent") && point_B["permanent"] == true)
                        delete!(locations["target_locations"], location_name_A)
                    else
                        if (haskey(point_A, "power") && haskey(point_B, "power"))
                            if point_A["power"] > point_B["power"]
                                delete!(locations["target_locations"], location_name_B)
                            else
                                delete!(locations["target_locations"], location_name_A)
                            end
                        elseif (haskey(point_A, "area") && haskey(point_B, "area"))
                            if point_A["area"] > point_B["area"]
                                delete!(locations["target_locations"], location_name_B)
                            else
                                delete!(locations["target_locations"], location_name_A)
                            end
                        else
                            if rand() > 0.5
                                delete!(locations["target_locations"], location_name_B)
                            else
                                delete!(locations["target_locations"], location_name_A)
                            end
                        end
                    end
                end
           end
        end
    end

    # Save pruned version of the locations list
    YAML.write_file("config/locations_pruned.yml", locations)

    # Init progress bar
    progress_bar = Progress(((parsed_args["end"] - parsed_args["start"] + 1) * length(locations["target_locations"])), 1, "Downloading:")

    # Temporal dataset per thread
    df_regional_daily = [ DataFrame(
        DateTime = Date[],
        Name = String[],
        Latitude = Float32[],
        Longitude = Float32[],
        Day = Int32[],
    ) for _ in 1:nthreads()]
    for i in 1:nthreads()
        for j::Int in 1:parsed_args["width"] * 2 * parsed_args["height"] * 2
            df_regional_daily[i][!, "Irradiance$(j)"] = Float32[]
            df_regional_daily[i][!, "Temp$(j)"] = Float32[]
            df_regional_daily[i][!, "TempMin$(j)"] = Float32[]
            df_regional_daily[i][!, "TempMax$(j)"] = Float32[]
            df_regional_daily[i][!, "Humidity$(j)"] = Float32[]
            df_regional_daily[i][!, "WindSpeed$(j)"] = Float32[]
            df_regional_daily[i][!, "WindSpeedMin$(j)"] = Float32[]
            df_regional_daily[i][!, "WindSpeedMax$(j)"] = Float32[]
            df_regional_daily[i][!, "WindDirection$(j)"] = Float32[]
            df_regional_daily[i][!, "WindX$(j)"] = Float32[]
            df_regional_daily[i][!, "WindY$(j)"] = Float32[]
            df_regional_daily[i][!, "WindXMin$(j)"] = Float32[]
            df_regional_daily[i][!, "WindYMin$(j)"] = Float32[]
            df_regional_daily[i][!, "WindXMax$(j)"] = Float32[]
            df_regional_daily[i][!, "WindYMax$(j)"] = Float32[]
            df_regional_daily[i][!, "Pressure$(j)"] = Float32[]
        end
    end
    df_point_daily = [ DataFrame(
        DateTime = Date[],
        Name = String[],
        Latitude = Float32[],
        Longitude = Float32[],
        Day = Int32[],
        Irradiance = Float32[]
    ) for _ in 1:nthreads()]
    df_point_hourly = [ DataFrame(
        DateTime = DateTime[],
        Name = String[],
        Latitude = Float32[],
        Longitude = Float32[],
        Day = Int32[],
        Hour = Int32[],
        Irradiance = Float32[]
    ) for _ in 1:nthreads()]

    # Information about the downloaded dataset
    println("\u001b[33;1m----------------------------------------------------------\u001b[0m")
    println("\u001b[34;1mNASA \u001b[31;1mPower \u001b[32;1mBot\u001b[0m ⛅ 🌞 ⚡ 🛰")
    println("Years range: $(parsed_args["start"]) - $(parsed_args["end"])")
    println("Locations:")
    for loc in keys(locations["target_locations"])
        println("          * $(loc)")
    end
    println("\u001b[33;1m----------------------------------------------------------\u001b[0m\n")

    # Downloading data
    @threads :dynamic for year in parsed_args["start"]:parsed_args["end"]
        # Region - daily
        for (location_name, point) in locations["target_locations"]
            data_regional = NASAPowerDownloader.download_regional(year, Utils.get_area(point["location"], parsed_args["width"], parsed_args["height"]), "daily", parsed_args["timeout"])
            long_name = data_regional["parameters"]["ALLSKY_SFC_SW_DWN"]["longname"]
            units = data_regional["parameters"]["ALLSKY_SFC_SW_DWN"]["units"]
            fill_value_regional = data_regional["header"]["fill_value"]

            X = Dict{String, Array{Float32}}()     # temporary
            features = data_regional["features"]
            for f in keys(features)
                irradiance = features[f]["properties"]["parameter"]["ALLSKY_SFC_SW_DWN"]
                temp = features[f]["properties"]["parameter"]["T2M"]
                temp_min = features[f]["properties"]["parameter"]["T2M_MIN"]
                temp_max = features[f]["properties"]["parameter"]["T2M_MAX"]
                humidity = features[f]["properties"]["parameter"]["RH2M"]
                wind_speed = features[f]["properties"]["parameter"]["WS10M"]
                wind_speed_min = features[f]["properties"]["parameter"]["WS10M_MIN"]
                wind_speed_max = features[f]["properties"]["parameter"]["WS10M_MAX"]
                wind_direction = features[f]["properties"]["parameter"]["WD10M"]
                pressure = features[f]["properties"]["parameter"]["PS"]
                for (i, t, t_min, t_max, h, ws, ws_min, ws_max, wd, p) in zip(irradiance, temp, temp_min, temp_max, humidity, wind_speed, wind_speed_min, wind_speed_max, wind_direction, pressure)
                    time = i[1]
                    if haskey(X, time)
                        push!(X[time], i[2])        # copying
                        push!(X[time], t[2])        # copying
                        push!(X[time], t_min[2])        # copying
                        push!(X[time], t_max[2])        # copying
                        push!(X[time], h[2])        # copying
                        push!(X[time], ws[2])       # copying
                        push!(X[time], ws_min[2])   # copying
                        push!(X[time], ws_max[2])   # copying
                        push!(X[time], wd[2])       # copying
                        push!(X[time], ws[2]*cosd(wd[2]))       # copying
                        push!(X[time], ws[2]*sind(wd[2]))       # copying
                        push!(X[time], ws_min[2]*cosd(wd[2]))   # copying
                        push!(X[time], ws_min[2]*sind(wd[2]))   # copying
                        push!(X[time], ws_max[2]*cosd(wd[2]))   # copying
                        push!(X[time], ws_max[2]*sind(wd[2]))   # copying

                        push!(X[time], p[2])        # copying
                    else
                        X[time] = [
                            i[2],
                            t[2],
                            t_min[2],
                            t_max[2],
                            h[2],
                            ws[2],
                            ws_min[2],
                            ws_max[2],
                            wd[2],
                            ws[2]*cosd(wd[2]),
                            ws[2]*sind(wd[2]),
                            ws_min[2]*cosd(wd[2]),
                            ws_min[2]*sind(wd[2]),
                            ws_max[2]*cosd(wd[2]),
                            ws_max[2]*sind(wd[2]),
                            p[2]
                        ]            # creating
                    end
                end
            end
            #println(X)

            for (t, value) in X
                t = Date(t, "yyyymmdd")
                push!(df_regional_daily[threadid()], [
                    t,
                    location_name,
                    point["location"][1],
                    point["location"][2],
                    dayofyear(t),
                    value...
                ])
            end

            # Point - daily
            data_point_daily = NASAPowerDownloader.download_point(year, Utils.Point(point["location"][1], point["location"][2]), "daily", parsed_args["timeout"])
            long_name = data_point_daily["parameters"]["ALLSKY_SFC_SW_DWN"]["longname"]
            units = data_point_daily["parameters"]["ALLSKY_SFC_SW_DWN"]["units"]
            fill_value_point_daily = data_point_daily["header"]["fill_value"]
            irradiance = data_point_daily["properties"]["parameter"]["ALLSKY_SFC_SW_DWN"]
            for (t, value) in irradiance
                t = Date(t, "yyyymmdd")
                push!(df_point_daily[threadid()], [
                    t,
                    location_name,
                    point["location"][1],
                    point["location"][2],
                    dayofyear(t),
                    value
                ])
            end

            # Point - hourly
            if year >= 2001
                data_point_hourly = NASAPowerDownloader.download_point(year, Utils.Point(point["location"][1], point["location"][2]), "hourly", parsed_args["timeout"])
                long_name = data_point_hourly["parameters"]["ALLSKY_SFC_SW_DWN"]["longname"]
                units = data_point_hourly["parameters"]["ALLSKY_SFC_SW_DWN"]["units"]
                fill_value_point_hourly = data_point_hourly["header"]["fill_value"]
                irradiance = data_point_hourly["properties"]["parameter"]["ALLSKY_SFC_SW_DWN"]
                for (t, value) in irradiance
                    t = DateTime(t, "yyyymmddHH")
                    push!(df_point_hourly[threadid()], [
                        t,
                        location_name,
                        point["location"][1],
                        point["location"][2],
                        dayofyear(t),
                        hour(t),
                        value
                    ])
                end
            end

            next!(progress_bar)
        end
    end

    # Combine dataframes
    X_all_daily = vcat(df_regional_daily...)
    y_all_daily = vcat(df_point_daily...)
    y_all_hourly = vcat(df_point_hourly...)

    # remove bad data
    if !isnothing(fill_value_regional)
        for j::Int in 1:parsed_args["width"] * 2 * parsed_args["height"] * 2
            X_all_daily = filter("Irradiance$(j)" => v -> v != fill_value_regional, X_all_daily)
        end
    end
    if !isnothing(fill_value_point_daily)
        y_all_daily = filter(:Irradiance => v -> v != fill_value_point_daily, y_all_daily)
    end
    if !isnothing(fill_value_point_hourly)
        y_all_hourly = filter(:Irradiance => v -> v != fill_value_point_hourly, y_all_hourly)
    end

    # Convert to kW/m^2
    y_all_hourly.Irradiance = y_all_hourly.Irradiance / 1000

    sort!(X_all_daily, [:DateTime, :Name])    # (timestep, patch, features)
    sort!(y_all_daily, [:DateTime, :Name])    # (timestep, patch, features)
    sort!(y_all_hourly, [:DateTime, :Name])   # (timestep, patch, features)

    println("\nDataset downloaded 🙂🙂🙂\n")

    # Show summary statistics
    println(describe(X_all_daily, :mean, :std, :median, :min, :max))
    println(describe(y_all_daily, :mean, :std, :median, :min, :max))
    println(describe(y_all_hourly, :mean, :std, :median, :min, :max))

    # write DataFrame out to CSV file
    CSV.write("dataset/X_all_daily.csv", X_all_daily)
    CSV.write("dataset/y_all_daily.csv", y_all_daily)
    CSV.write("dataset/y_all_hourly.csv", y_all_hourly)

    # Filter
    X_filtered = filter(:Name => n -> n == "Adani Green Energy Tamilnadu Limited", filter(:DateTime => d -> Dates.month(d) == 11, X_all_daily))
    Y_filtered = filter(:Name => n -> n == "Indira Paryavaran Bhawan", filter(:DateTime => d -> Dates.month(d) == 11, y_all_daily))
    # center
    # Target 1

    # north - south
    # Indira Paryavaran Bhawan
    # Adani Green Energy Tamilnadu Limited

    # west - east
    # Bitta Solar Power Plant
    # Rewa Ultra Mega Solar

    Utils.create_heatmap(X_filtered, Y_filtered, round(Int, parsed_args["width"] * 2), round(Int, parsed_args["height"] * 2), "Irradiance")
    Utils.create_heatmap(X_filtered, Y_filtered, round(Int, parsed_args["width"] * 2), round(Int, parsed_args["height"] * 2), "Temp")
    Utils.create_heatmap(X_filtered, Y_filtered, round(Int, parsed_args["width"] * 2), round(Int, parsed_args["height"] * 2), "TempMin")
    Utils.create_heatmap(X_filtered, Y_filtered, round(Int, parsed_args["width"] * 2), round(Int, parsed_args["height"] * 2), "TempMax")
    Utils.create_heatmap(X_filtered, Y_filtered, round(Int, parsed_args["width"] * 2), round(Int, parsed_args["height"] * 2), "Humidity")
    Utils.create_heatmap(X_filtered, Y_filtered, round(Int, parsed_args["width"] * 2), round(Int, parsed_args["height"] * 2), "WindSpeed", windDir = true)
    Utils.create_heatmap(X_filtered, Y_filtered, round(Int, parsed_args["width"] * 2), round(Int, parsed_args["height"] * 2), "WindSpeedMin", windDir = true)
    Utils.create_heatmap(X_filtered, Y_filtered, round(Int, parsed_args["width"] * 2), round(Int, parsed_args["height"] * 2), "WindSpeedMax", windDir = true)
    Utils.create_heatmap(X_filtered, Y_filtered, round(Int, parsed_args["width"] * 2), round(Int, parsed_args["height"] * 2), "WindDirection", windDir = true)
    Utils.create_heatmap(X_filtered, Y_filtered, round(Int, parsed_args["width"] * 2), round(Int, parsed_args["height"] * 2), "Pressure")
end

# check the num. of threads
if Threads.nthreads() == 1
    println("Warning: The number of threads is only 1. It is recommended to use at least 2 threads.")
end

@time main()


text:
using GPR

function experimentVarInt(config, id, eid, mechanism, varargs...)
    testdfs = config["datasets"][2]  # Each thread operates on its own dataset -> no races
    testdf = testdfs.df[id][shuffle(1:nrow(testdfs.df[id]))[1:config["testsamples"]], :]
    xtest_future_true = [CState(x) for x in testdf.sfuture]
    # Add noise to the dataset
    applynoise!(testdf, config["Σ"], eid, config["Δtsim"], varargs...)
    xtest_old = [CState(x) for x in testdf.sold]
    predictedstates = Vector{CState{Float64, length(mechanism.bodies)}}()
    for i in 1:length(xtest_old)
        predictedstate = predictdynamics(mechanism, xtest_old[i], config["simsteps"])
        push!(predictedstates, predictedstate)
    end
    return predictedstates, xtest_future_true
end

function experimentVarIntP1(config, id)
    mechanism = simplependulum2D(1, Δt=0.01, threadlock = config["mechanismlock"])[2]
    return experimentVarInt(config, id, "P1", mechanism, mechanism.bodies[1].shape.rh[2])
end

function experimentVarIntP2(config, id)
    mechanism = doublependulum2D(1, Δt=0.01, threadlock = config["mechanismlock"])[2]
    return experimentVarInt(config, id, "P2", mechanism, mechanism.bodies[1].shape.xyz[3], mechanism.bodies[2].shape.xyz[3])
end

function experimentVarIntCP(config, id)
    mechanism = cartpole(1, Δt=0.01, threadlock = config["mechanismlock"])[2]
    return experimentVarInt(config, id, "CP", mechanism, mechanism.bodies[2].shape.rh[2])
end

function experimentVarIntFB(config, id)
    mechanism = fourbar(1, Δt=0.01, threadlock = config["mechanismlock"])[2]
    return experimentVarInt(config, id, "FB", mechanism, mechanism.bodies[1].shape.xyz[3])
end

ext: jl


text:
using JuMP
using StructJuMP
using ECOS
using Cbc
using Base.Test

misocp_solver = CbcSolver()
socp_solver = ECOS.ECOSSolver(verbose=false)

@testset "[Benders] conicconstraintdata with more variables in parent" begin
    m = StructuredModel(num_scenarios=1)
    @variable(m, x[1:2])
    @objective(m, :Min, sum(x))

    bl = StructuredModel(parent=m, id=1)
    @variable(bl, y)
    @constraint(bl, 4y + 5x[1] + 6x[2] >= 2)
    @objective(bl, :Max, 3y)

    c, A, B, b, var_cones, con_cones, v = StructJuMP.conicconstraintdata(bl)
    @test c == [-3]
    @test A == [5 6]
    @test B == reshape([4], 1, 1)
    @test b == [2]
    @test length(var_cones) == 1
    @test var_cones[1][1] == :Free
    @test collect(var_cones[1][2]) == [1]
    @test con_cones[1][1] == :NonPos
    @test collect(con_cones[1][2]) == [1]
    @test v == [:Cont]
end

@testset "[Benders] Empty scenario test" begin

    m = StructuredModel(num_scenarios=0)
    @variable(m, x, Int)
    @constraint(m, x <= 4)
    @objective(m, :Min, -5*x)

    output = BendersBridge(m, misocp_solver, socp_solver)

    @test output[1] == :Optimal
    @test isapprox(output[2], -20.0)

end

@testset "[Benders] Infeasible problem test" begin

    numScen = 1
    m = StructuredModel(num_scenarios=numScen)

    @variable(m, x, Int)

    @constraint(m, x <= 1)
    @objective(m, :Min, -5*x)

    bl = StructuredModel(parent=m, id=1)
    @variable(bl, y1 >= 2)
    @variable(bl, y2 <= 2)
    @constraint(bl, x >= y1)
    @constraint(bl, norm(y1) <= y2)
    @objective(bl, :Min, 2*y1 + y2)

    output = BendersBridge(m, misocp_solver, socp_solver)

    @test output[1] == :Infeasible

end

@testset "[Benders] Infeasibility cut execution test #1" begin

    numScen = 1
    m = StructuredModel(num_scenarios=numScen)

    @variable(m, x, Int)

    @constraint(m, x <= 4)
    @objective(m, :Min, -5*x)

    bl = StructuredModel(parent=m, id=1)
    @variable(bl, y1 >= 0)
    @variable(bl, y2 <= 2)
    @constraint(bl, x <= y1)
    @constraint(bl, norm(y1) <= y2)
    @objective(bl, :Min, 2*y1 + y2)

    output = BendersBridge(m, misocp_solver, socp_solver)

    @test output[1] == :Optimal
    @test isapprox(output[2], -4.0)

end

@testset "[Benders] Optimality cut execution test #1" begin

    numScen = 1
    m = StructuredModel(num_scenarios=numScen)

    @variable(m, x, Int)

    @constraint(m, x <= 4)
    @objective(m, :Min, -5*x)

    bl = StructuredModel(parent=m, id=1)
    @variable(bl, y1 >= 2)
    @variable(bl, y2 <= 4)
    @constraint(bl, x <= y1)
    @constraint(bl, norm(y1) <= y2)
    @objective(bl, :Min, 2*y1 + y2)

    output = BendersBridge(m, misocp_solver, socp_solver)

    @test output[1] == :Optimal
    @test isapprox(output[2], -8.0)

end


text:
using JuMP,CPLEX,Statistics,MAT,Redis,Printf

#connect to redi
conn = RedisConnection()

function min_(X,model)
    if(length(X)!=2)
        error("supperted only two sided minimum")
    end
    d=@variable(model,binary = true)
    dmin=@variable(model)
    @constraints(model, begin
           dmin<=X[1]
           dmin<=X[2]
           dmin>=X[1]-1000*(d)
           dmin>=X[2]-1000*(1-d)
    end)
    return dmin
end

function rate(F,T,h)
    return [(1-F)*T[1,h],
        (1-F)*T[3,h],
        (F)*T[2,h],
        ];
end

jump=[[-1,1,0],
      [1,0,-1],
      [0,-1,1]];

isMIP=true
delta=10^2
dt=10^-2  #time step lento
dt2=0.5*dt #time step veloce
iH=1
oH=parse(Int64, ARGS[3])
H=oH*(1+iH)
#MU=round.([rand()*100,rand()*100],digits=2)
MU=[parse(Float64, ARGS[1]),parse(Float64, ARGS[2])]
#X0=[rand(1:60),rand(1:60),0]
X0=parse.(Int64, mget(conn, "driver","acquire","exec"))
tgt=round((sum(X0)*((1/MU[1])/(1/MU[1]+1/MU[2]))),digits=2)
tgt=tgt
N=3
Nr=3
isStarted=false

model = Model(CPLEX.Optimizer)
set_optimizer_attribute(model, "CPX_PARAM_SCRIND", 0)
set_optimizer_attribute(model, "CPX_PARAM_THREADS", 4)

@variable(model,X[i=1:N,j=0:H]>=0)
@variable(model,T[i=1:Nr,j=0:H]>=0)
@variable(model,E_abs[i=1:H]>=0)
@variable(model,S[i=1:2]>=1)

#@constraint(model,[i=1:N,j=0:H],X[i,j]<=sum(X0))
#@constraint(model,S[1]<=32)
#@constraint(model,S[2]<=32)

@variable(model,D[i=1:2,j=0:H]>=0)
@constraint(model,[h=0:H],D[1,h]<=X[2,h])
@constraint(model,[h=0:H],D[1,h]<=S[1]-X[3,h])
@constraint(model,[h=0:H],D[2,h]<=X[3,h])
@constraint(model,[h=0:H],D[2,h]<=S[2])

#INIT STATE constraint
@constraint(model,X0c[i=1:N],X[i,0]==X0[i])
#Constraints per lo sforzo di controllo
@constraint(model,SC_u[i=1:2],S[i]<=1000);
@constraint(model,SC_l[i=1:2],S[i]>=1);

for h=0:H
    @constraint(model,T[1,h]==MU[1]*X[1,h])
    if(h<=3)
        @constraint(model,T[2,h]==delta*min_([X[2,h],S[1]-X[3,h]],model))
        @constraint(model,T[3,h]==MU[2]*min_([X[3,h],S[2]],model))
    else
        @constraint(model,T[2,h]==delta*D[1,h])
        @constraint(model,T[3,h]==MU[2]*D[2,h])
    end
end

#integration of the fast part
h=0
for oIdx=1:(oH)
    for iIdx=1:iH
        local dy=(jump'*rate(1,T,h))'
        @constraint(model,[i=1:N],X[i,h+1]==dy[i]*dt2+X[i,h])
        global h=h+1;
    end
    #integrate the slow system for one step
    local dy2=(jump'*rate(0,T,h))'
    @constraint(model,[i=1:N],X[i,h+1]==dy2[i]*dt+X[i,h])
    global h=h+1;
end

#add absolute value constraint
for i=1:H
    @constraints(model, begin
           E_abs[i]>=(X[1,i]-tgt)
           E_abs[i]>=-(X[1,i]-tgt)
           #E_abs[i,2]==((X[2,i]+X[3,i])-(sum(X0)-tgt))
    end)
end

for i=0:H
    @constraints(model, begin
           X[3,i]<=S[1]
    end)
end

@objective(model,Min,sum(E_abs)+0.1*sum(S))

X0=parse.(Int64, mget(conn, "driver","acquire","exec"))
prevS=[0,0];
while(true)
    if(isStarted==true)
        set_normalized_rhs(SC_l[1],max(1,prevS[1]-1))
        set_normalized_rhs(SC_l[2],max(1,prevS[2]-1))
        set_normalized_rhs(SC_u[1],min(100,prevS[1]+5))
        set_normalized_rhs(SC_u[2],min(100,prevS[2]+1))
    end

    if(isStarted==false)
        global isStarted=true
        set(conn, "ctrl","1")
    end
    for j=1:N
        set_normalized_rhs(X0c[j],X0[j])
    end
    JuMP.optimize!(model)
    local soltime=MOI.get(model, MOI.SolveTime())

    println(round.(value.(S)),X0,soltime)

    Sopt=Int.(round.(value.(S)))
    global prevS=Sopt

    set(conn, "swthread",@sprintf("%d",Sopt[1]))
    set(conn, "hwthread",@sprintf("%d",Sopt[2]))
    global X0=parse.(Int64, mget(conn, "driver","acquire","exec"))
end

#Controlled layered three-tier system by varying #threads HW/SW at each tier. Simulation (i.e., 100 runs) Vs Target reguirement.
#println(mean(mean(XS[:,:,:],dims=3)[end-100:end,1])," ",abs(mean(mean(XS[:,:,:],dims=3)[end-100:end,1])-tgt)*100/tgt," ",mean(stimes)," ",maximum(stimes))
#println(mean(mean(Shw,dims=3),dims=1))
# plot(t,mean(XS[:,:,:],dims=3)[:,:],label = "Controlled-simulation",legend=false)
# hline!([tgt],label = "Target",reuse=true,legend=false)
# xlabel!("Time(s)")
# ylabel!("QueueLength")

#p=plot(t,mean(Shw,dims=3)[:,:],reuse = false,title="SW/HW Core")

# matwrite("LP_LQNMPC.mat", Dict(
# 	"XS" => XS,
# 	"NT" => Tsw,
#     "NC" => Shw,
#     "Time"=>stimes
# ); compress = false)

max_stars_repo_name: giulio-garbi/fitMS

max_issues_repo_name: giulio-garbi/fitMS

max_forks_repo_name: giulio-garbi/fitMS


text:
tickPrice= function(tickerId::Int, field::String, price::Float64, size::Float64, attrib::TickAttrib)

tickSize= function(tickerId::Int, field::String, size::Float64)

tickOptionComputation= function(tickerId::Int, tickType::String, tickAttrib::Union{Int,Nothing}, impliedVol::Union{Float64,Nothing}, delta::Union{Float64,Nothing}, optPrice::Union{Float64,Nothing}, pvDividend::Union{Float64,Nothing}, gamma::Union{Float64,Nothing}, vega::Union{Float64,Nothing}, theta::Union{Float64,Nothing}, undPrice::Union{Float64,Nothing})

tickGeneric= function(tickerId::Int, tickType::String, value::Float64)

tickString= function(tickerId::Int, tickType::String, value::String)

tickEFP= function(tickerId::Int, tickType::String, basisPoints::Float64, formattedBasisPoints::String, totalDividends::Float64, holdDays::Int, futureLastTradeDate::String, dividendImpact::Float64, dividendsToLastTradeDate::Float64)

orderStatus= function(orderId::Int, status::String, filled::Float64, remaining::Float64, avgFillPrice::Float64, permId::Int, parentId::Int, lastFillPrice::Float64, clientId::Int, whyHeld::String, mktCapPrice::Float64)

openOrder= function(orderId::Int, contract::Contract, order::Order, orderstate::OrderState)

openOrderEnd= function()

updateAccountValue= function(key::String, val::String, currency::String, accountName::String)

updatePortfolio= function(contract::Contract, position::Float64, marketPrice::Float64, marketValue::Float64, averageCost::Float64, unrealizedPNL::Float64, realizedPNL::Float64, accountName::String)

updateAccountTime= function(timeStamp::String)

accountDownloadEnd= function(accountName::String)

nextValidId= function(orderId::Int)

contractDetails= function(reqId::Int, contractDetails::ContractDetails)

bondContractDetails= function(reqId::Int, contractDetails::ContractDetails)

contractDetailsEnd= function(reqId::Int)

execDetails= function(reqId::Int, contract::Contract, execution::Execution)

execDetailsEnd= function(reqId::Int)

error= function(id::Union{Int,Nothing}, errorCode::Union{Int,Nothing}, errorString::String)

updateMktDepth= function(id::Int, position::Int, operation::Int, side::Int, price::Float64, size::Float64)

updateMktDepthL2= function(id::Int, position::Int, marketMaker::String, operation::Int, side::Int, price::Float64, size::Float64, isSmartDepth::Bool)

updateNewsBulletin= function(msgId::Int, msgType::Int, newsMessage::String, originExch::String)

managedAccounts= function(accountsList::String)

receiveFA= function(faDataType::FaDataType, xml::String)

historicalData= function(reqId::Int, bar::DataFrame)

scannerParameters= function(xml::String)

scannerData= function(reqId::Int, rank::Vector{Int}, contractDetails::Vector{ContractDetails}, distance::Vector{String}, benchmark::Vector{String}, projection::Vector{String}, legsStr::Vector{String})

realtimeBar= function(reqId::Int, time::Int, open::Float64, high::Float64, low::Float64, close::Float64, volume::Float64, wap::Float64, count::Int)

currentTime= function(time::Int)

fundamentalData= function(reqId::Int, data::String)

deltaNeutralValidation= function(reqId::Int, deltaNeutralContract::DeltaNeutralContract)

tickSnapshotEnd= function(reqId::Int)

marketDataType= function(reqId::Int, marketDataType::MarketDataType)

commissionReport= function(commissionReport::CommissionReport)

position= function(account::String, contract::Contract, position::Float64, avgCost::Float64)

positionEnd= function()

accountSummary= function(reqId::Int, account::String, tag::String, value::String, currency::String)

accountSummaryEnd= function(reqId::Int)

verifyMessageAPI= function(apiData::String)

verifyCompleted= function(isSuccessful::Bool, errorText::String)

displayGroupList= function(reqId::Int, groups::String)

displayGroupUpdated= function(reqId::Int, contractInfo::String)

verifyAndAuthMessageAPI= function(apiData::String, xyzChallange::String)

verifyAndAuthCompleted= function(isSuccessful::Bool, errorText::String)

positionMulti= function(reqId::Int, account::String, modelCode::String, contract::Contract, position::Float64, avgCost::Float64)

positionMultiEnd= function(reqId::Int)

accountUpdateMulti= function(reqId::Int, account::String, modelCode::String, key::String, value::String, currency::String)

accountUpdateMultiEnd= function(reqId::Int)

securityDefinitionOptionalParameter= function(reqId::Int, exchange::String, underlyingConId::Int, tradingClass::String, multiplier::String, expirations::Vector{String}, strikes::Vector{Float64})

securityDefinitionOptionalParameterEnd= function(reqId::Int)

softDollarTiers= function(reqId::Int, tiers::Vector{SoftDollarTier})

familyCodes= function(familyCodes::Vector{FamilyCode})

symbolSamples= function(reqId::Int, contractDescriptions::Vector{ContractDescription})

mktDepthExchanges= function(depthMktDataDescriptions::DataFrame)

tickNews= function(tickerId::Int, timeStamp::Int, providerCode::String, articleId::String, headline::String, extraData::String)

smartComponents= function(reqId::Int, theMap::DataFrame)

tickReqParams= function(tickerId::Int, minTick::Union{Float64,Nothing}, bboExchange::String, snapshotPermissions::Int)

newsProviders= function(newsProviders::DataFrame)

newsArticle= function(requestId::Int, articleType::Int, articleText::String)

historicalNews= function(requestId::Int, time::String, providerCode::String, articleId::String, headline::String)

historicalNewsEnd= function(requestId::Int, hasMore::Bool)

headTimestamp= function(reqId::Int, headTimestamp::String)

histogramData= function(reqId::Int, data::DataFrame)

historicalDataUpdate= function(reqId::Int, bar::NamedTuple)

rerouteMktDataReq= function(reqId::Int, conid::Int, exchange::String)

rerouteMktDepthReq= function(reqId::Int, conid::Int, exchange::String)

marketRule= function(marketRuleId::Int, priceIncrements::DataFrame)

pnl= function(reqId::Int, dailyPnL::Float64, unrealizedPnL::Float64, realizedPnL::Float64)

pnlSingle= function(reqId::Int, pos::Int, dailyPnL::Float64, unrealizedPnL::Union{Float64,Nothing}, realizedPnL::Union{Float64,Nothing}, value::Float64)

historicalTicks= function(reqId::Int, ticks::DataFrame, done::Bool)

historicalTicksBidAsk= function(reqId::Int, ticks::DataFrame, done::Bool)

historicalTicksLast= function(reqId::Int, ticks::DataFrame, done::Bool)

tickByTickAllLast= function(reqId::Int, tickType::Int, time::Int, price::Float64, size::Float64, attribs::TickAttribLast, exchange::String, specialConditions::String)

tickByTickBidAsk= function(reqId::Int, time::Int, bidPrice::Float64, askPrice::Float64, bidSize::Float64, askSize::Float64, attribs::TickAttribBidAsk)

tickByTickMidPoint= function(reqId::Int, time::Int, midPoint::Float64)

orderBound= function(orderId::Int, apiClientId::Int, apiOrderId::Int)

completedOrder= function(contract::Contract, order::Order, orderState::OrderState)

completedOrdersEnd= function()

replaceFAEnd= function(reqId::Int, text::String)


text:
using DelimitedFiles
using Test
using Base.Iterators

# Those corresponding verifiable data are obtained from orginal fortran subroutines at
# http://www.bosai.go.jp/study/application/dc3d/DC3Dhtml_E.html
@testset "Test Okada's displacement" begin
    @testset "Normal computations" begin
        u_truth = readdlm(joinpath(@__DIR__, "data/test_dc3d.dat"), ',', Float64)
        pos = product(-40.: -30., 20.: 30., 10.: 15.)
        fu = (x) -> dc3d(reverse(x)..., 2. /3, 50., 70., [-80., 120.], [-30., 25.], [200., -150., 100.])
        u_cal = map(fu, pos) |> vec
        cache = GeoGreensFunctions.dc3d_cache(Float64)
        for (i, x) in enumerate(pos)
            dc3d(reverse(x)..., 2. /3, 50., 70., [-80., 120.], [-30., 25.], [200., -150., 100.], cache)
            @test u_truth[i,:] ≈ u_cal[i] ≈ cache[1]
        end
    end

    @testset "Negative depth" begin
        u = dc3d(10., 20., 30., 2. /3, 50., 70., [-80., 120.], [-30., 25.], [200., -150., 100.])
        @test isapprox(zeros(Float64, 12), u, rtol=1e-8)
    end

    @testset "Singularity" begin
        u_singular = zeros(Float64, 12)

        u = dc3d(-80., 0., -50., 2. /3, 50., 70., [-80., 120.], [-30., 25.], [200., -150., 100.])
        @test isapprox(u_singular, u, atol=1e-8)

        u = dc3d(-70., 0., -10., 2. /3, 10., 90., [-80., 120.], [-20., 20.], [200., -150., 100.])
        @test isapprox(u_singular, u, atol=1e-8)

        u = dc3d(-80., 20., -50., 2. /3, 50., 70., [-80., 120.], [-30., 25.], [200., -150., 100.])
        u_truth = [
            -32.60576975, 59.39436703, 1.61816100, -0.43682428, 1.66503564, 0.27751175,
            0.64126234, -1.05836342, 0.03630938, -0.16974972, 0.61800758, 0.22842501]
        @test isapprox(u_truth, u, rtol=1e-8)

        u = dc3d(-10., 20., -50., 2. /3, 50., 90., [-80., 120.], [-30., 25.], [200., -150., 100.])
        u_truth = [
            -56.04511679, 52.38458483, 40.87786501, -0.10140811, -0.09099896, 0.02371607,
            1.60327214, -0.63369117, -1.13609685, 0.05783881, 0.88665943, 0.21946899,
        ]
        @test isapprox(u_truth, u, rtol=1e-8)

        u = dc3d(-10., 20., 0., 2. /3, 0., 0., [-80., 120.], [-30., 25.], [200., -150., 100.])
        # take a look: https://github.com/JuliaLang/julia/issues/23376
        @test isapprox(u_singular, u, atol=1e-8)

        u = dc3d(-100., 20., 0., 2. /3, 0., 0., [-80., -90.], [20., 25.], [200., -150., 100.])
        @test isapprox(u_singular, u, atol=1e-8)

        u = dc3d(-80., 20., 0., 2. /3, 0., 0., [-80., 100.], [25., 35.], [200., -150., 100.])
        @test isapprox(u_singular, u, atol=1e-8)
    end
end


text:
function _interaction_type(N::T) where {T <: AbstractEcologicalNetwork}
    return eltype(N.edges)
end

function _species_type(N::T) where {T <: AbstractBipartiteNetwork}
    return eltype(N.B)
end

function _species_type(N::T) where {T <: AbstractUnipartiteNetwork}
    return eltype(N.S)
end

function _check_species_validity(::Type{T}) where {T <: Any}
  throw(ArgumentError("The type $(T) is not an allowed species type"))
end

_check_species_validity(::Type{T}) where {T <: Symbol} = nothing
_check_species_validity(::Type{T}) where {T <: AbstractString} = nothing
_check_species_validity(::Type{T}) where {T <: String} = nothing

"""
    richness(N::AbstractEcologicalNetwork, i::Int64)

Returns the number of species on either side of the network. The value of `i`
can be `1` (top-level species) or `2` (bottom-level species), as in the
`species` function.
"""
function richness(N::AbstractEcologicalNetwork; dims::Union{Nothing,Integer}=nothing)
  return length(species(N; dims=dims))
end

"""
    species(N::AbstractBipartiteNetwork)

Returns an array of all species in a bipartite network. The order of the species
corresponds to the order of rows (top level) and columns (bottom level) of the
adjacency matrix, in this order.
"""
function species(N::AbstractBipartiteNetwork; dims::Union{Nothing,Integer}=nothing)
  isnothing(dims) && return vcat(N.T, N.B)
  dims == 1 && return N.T
  dims == 2 && return N.B
  throw(ArgumentError("dims can only be 1 (top species) or 2 (bottom species), or `nothing` (all species), you used $(dims)"))
end

"""
    species(N::AbstractUnipartiteNetwork; dims::Int64=1)

Returns an array of species on either side of a unipartite network. In a
unipartite network, species are the same on either size. This function is
nevertheless useful when you want to write code that takes either side of the
network in a general way.
"""
function species(N::AbstractUnipartiteNetwork; dims::Union{Nothing,Integer}=nothing)
  return N.S
end

function _species_objects(N::AbstractUnipartiteNetwork)
  return (N.S,)
end

function _species_objects(N::AbstractBipartiteNetwork)
  return (N.T, N.B)
end

"""
    has_interaction{(N::AbstractEcologicalNetwork, i::NT, j::NT)

This function returns `true` if the interaction between `i` and `j` is not 0.
This refers to species by their names/values, and is the recommended way to test
for the presence of an interaction.

Use `N[i,j]` if you need to get the value of the interaction.
"""
function has_interaction(N::AbstractEcologicalNetwork, i::NT, j::NT) where {NT}
  _check_species_validity(NT)
  @assert i ∈ species(N; dims=1)
  @assert j ∈ species(N; dims=2)
  i_pos = something(findfirst(isequal(i), species(N; dims=1)),0)
  j_pos = something(findfirst(isequal(j), species(N; dims=2)),0)
  return has_interaction(N, i_pos, j_pos)
end

"""
    has_interaction(N::AbstractEcologicalNetwork, i::Int64, j::Int64)

This function returns `true` if the interaction between `i` and `j` is not 0.
This refers to species by their position instead of their name, and is not
recommended as the main solution. This is used internally by a few functions,
but is exported because it may be of general use.
"""
function has_interaction(N::AbstractEcologicalNetwork, i::Int64, j::Int64)
  # We need to make sure that the interaction is somewhere within the network
  @assert i <= size(N.edges, 1)
  @assert j <= size(N.edges, 2)
  # This should be reasonably general...
  return !iszero(N[i,j])
end

"""
    nodiagonal!(N::AbstractUnipartiteNetwork)

Modifies the network so that its diagonal is set to the appropriate zero.
"""
function nodiagonal!(N::AbstractUnipartiteNetwork)
  N.edges[diagind(N.edges)] .= zero(eltype(N.edges))
end

"""
    nodiagonal(N::AbstractUnipartiteNetwork)

Returns a copy of the network with its diagonal set to zero.
"""
function nodiagonal(N::AbstractUnipartiteNetwork)
  Y = copy(N)
  nodiagonal!(Y)
  return Y
end

"""
    nodiagonal(N::AbstractBipartiteNetwork)

Returns a *copy* of the network (because the diagonal of a bipartite network is
never a meaningful notion). This function is clearly useless, but allows to
write general code for all networks types when a step requires removing the
diagonal.
"""
function nodiagonal(N::AbstractBipartiteNetwork)
  return copy(N)
end

"""
    nodiagonal!(N::AbstractBipartiteNetwork)

Does nothing.
"""
function nodiagonal!(N::AbstractBipartiteNetwork)
end

"""
    interactions(N::AbstractEcologicalNetwork)

Returns the interactions in the ecological network. Interactions are returned as
an array of named tuples. *A minima*, these have fields `from` and `to`. For
networks that are probabilistic, there is a `probability` field. For networks
that are quantitative, there is a `strength` field. This functions allows to
iterate over interactions in a network in a convenient way.
"""
function interactions(N::BinaryNetwork)
    non_zero = count(!iszero, N.edges) # Number of non-zero entries in the matrix
    t_int, t_spe = _interaction_type(N), _species_type(N)
    edges_accumulator = Vector{NamedTuple{tuple(:from, :to),Tuple{t_spe, t_spe}}}(undef, non_zero)
    sp1 = species(N; dims=1)
    sp2 = species(N; dims=2)
    I, J, V = findnz(N.edges)
    for i in eachindex(I)
        edges_accumulator[i] = (from = sp1[I[i]], to = sp2[J[i]])
    end
    return unique(edges_accumulator)
end

function interactions(N::QuantitativeNetwork)
    non_zero = count(!iszero, N.edges) # Number of non-zero entries in the matrix
    t_int, t_spe = _interaction_type(N), _species_type(N)
    edges_accumulator = Vector{NamedTuple{tuple(:from, :to, :strength),Tuple{t_spe, t_spe, t_int}}}(undef, non_zero)
    sp1 = species(N; dims=1)
    sp2 = species(N; dims=2)
    I, J, V = findnz(N.edges)
    for i in eachindex(I)
        edges_accumulator[i] = (from = sp1[I[i]], to = sp2[J[i]], strength = V[i])
    end
    return unique(edges_accumulator)
end

function interactions(N::ProbabilisticNetwork)
    non_zero = count(!iszero, N.edges) # Number of non-zero entries in the matrix
    t_int, t_spe = _interaction_type(N), _species_type(N)
    edges_accumulator = Vector{NamedTuple{tuple(:from, :to, :probability),Tuple{t_spe, t_spe, t_int}}}(undef, non_zero)
    sp1 = species(N; dims=1)
    sp2 = species(N; dims=2)
    I, J, V = findnz(N.edges)
    for i in eachindex(I)
        edges_accumulator[i] = (from = sp1[I[i]], to = sp2[J[i]], probability = V[i])
    end
    return unique(edges_accumulator)
end

function adjacency(N::T) where {T <: AbstractEcologicalNetwork}
  return Array(N.edges)
end


text:
using OceanTurb

using OceanTurb: @add_standard_model_fields, AbstractParameters, AbstractModel

import OceanTurb.KPP: ∂B∂z

const nsol = 4
@solution U V T S

Base.@kwdef struct ModifiedPacanowskiPhilanderParameters{T} <: AbstractParameters
    Cν₀ :: T = 1e-4
    # Cν₋ :: T = 1e-1
    Cν₋ :: T = 0.1
    # Cν₋ :: T = 1.
    Pr  :: T = 1.0
    Riᶜ :: T = 0.25
    ΔRi :: T = 1e-1
end

struct ModifiedPacanowskiPhilanderModel{TS, G, T} <: AbstractModel{TS, G, T}
    @add_standard_model_fields
    parameters :: ModifiedPacanowskiPhilanderParameters{T}
    constants  :: Constants{T}
end

function ModifiedPacanowskiPhilanderModel(; N, L,
          grid = UniformGrid(N, L),
     constants = Constants(),
    parameters = Parameters(),
       stepper = :ForwardEuler,
           bcs = BoundaryConditions((ZeroFluxBoundaryConditions() for i=1:nsol)...)
    )

    solution = Solution((CellField(grid) for i=1:nsol)...)
    K = (U=KU, V=KV, T=KT, S=KS)
    R = (U=RU, V=RV, T=RT, S=RS)
    eqn = Equation(K=K, R=R)
    lhs = OceanTurb.build_lhs(solution)

    timestepper = Timestepper(stepper, eqn, solution, lhs)

    return ModifiedPacanowskiPhilanderModel(Clock(), grid, timestepper, solution, bcs, parameters, constants)
end

function local_richardson(U, V, T, S, g, α, β, i)
    Bz = ∂B∂z(T, S, g, α, β, i)
    S² = ∂z(U, i)^2 + ∂z(V, i)^2

    if S² == 0 && Bz == 0 # Alistair Adcroft's theorem
        return 0
    else
        return Bz / S²
    end
end

local_richardson(m, i) = local_richardson(m.solution.U, m.solution.V, m.solution.T,
                                          m.solution.S, m.constants.g, m.constants.α,
                                          m.constants.β, i)

tanh_step(x) = (1 - tanh(x)) / 2

# KU(Ri, ν₀, ν₋, Riᶜ, ΔRi) = ν₀ + ν₋ * tanh_step((Ri - Riᶜ) / ΔRi)
# KT() = 0
# KU(m, i) = KU(local_richardson(m, i), m.parameters.Cν₀, m.parameters.Cν₋,
#               m.parameters.Riᶜ, m.parameters.ΔRi)
# KT(m, i) = KT()

# KU() = 0
# KT(Ri, ν₀, ν₋, Riᶜ, ΔRi, Pr) = (ν₀ + ν₋ * tanh_step((Ri - Riᶜ) / ΔRi)) / Pr
# KU(m, i) = KU()
# KT(m, i) = KT(local_richardson(m, i), m.parameters.Cν₀, m.parameters.Cν₋,
#               m.parameters.Riᶜ, m.parameters.ΔRi, m.parameters.Pr)

KU(Ri, ν₀, ν₋, Riᶜ, ΔRi) = ν₀ + ν₋ * tanh_step((Ri - Riᶜ) / ΔRi)
KT(Ri, ν₀, ν₋, Riᶜ, ΔRi, Pr) = KU(Ri, ν₀, ν₋, Riᶜ, ΔRi) / Pr

KU(m, i) = KU(local_richardson(m, i), m.parameters.Cν₀, m.parameters.Cν₋, m.parameters.Riᶜ, m.parameters.ΔRi)
KT(m, i) = KT(local_richardson(m, i), m.parameters.Cν₀, m.parameters.Cν₋, m.parameters.Riᶜ, m.parameters.ΔRi, m.parameters.Pr)

const KV = KU
const KS = KT

RU(m, i) =   m.constants.f * m.solution.V[i]
RV(m, i) = - m.constants.f * m.solution.U[i]
RT(m, i) = 0
RS(m, i) = 0


text:
"""
  τ2λ(τ::Vector{Float64}, λ::Vector{Float64})

Suplement the vector of λ patrameters of Marshall-Olkin copula, given some of those
parameters and a vector fo Kendall's τ correlations. Wroks fo 2,3 variate MO copula
"""
function τ2λ(τ::Vector{Float64}, λ::Vector{Float64})
  if length(τ) == 1
    return [λ[1],λ[2],(λ[1]+λ[2])*τ[1]/(1-τ[1])]
  else
    t = τ./(1 .-τ)
    M = Matrix(1.0I, 3, 3) .-(ones(3,3) .-Matrix(1.0I, 3, 3)) .*t
    fm = hcat([1. 1. 0.; 0. 1. 1.; 1. 0. 1.].*t, .-[1.; 1.; 1.])
    ret = map(x -> (x > 0) ? (x) : (0.000000001), inv(M)*fm*λ)
    [λ[1:3]..., ret..., λ[end]]
  end
end

function moρ2τ(ρ::Float64)
  f(τ) = 1/2 .*sin.(τ[1]*pi/2)+τ[1]/2 - ρ[1]
  fzero(f, -0.999, 0.999)
end

#=
# some Marshall-Olkin copulas helpers

function getmoλ(λ::Vector{Float64}, ind::Vector{Int})
  n = floor(Int, log(2, length(λ)+1))
  s = collect(combinations(1:n))
  λ[[ind == a for a in s]]
end

function setmoλ!(λ::Vector{Float64}, ind::Vector{Int}, a::Float64)
  n = floor(Int, log(2, length(λ)+1))
  s = collect(combinations(1:n))
  λ[[ind == a for a in s]] = a
end
=#