text:
#Abstract runtime policy
abstract AbstractPolicy

#Factory function, currently supported type is :mh
policy(s::Symbol,args...;keyargs...) = _policy(Val{s},args...;keyargs...)


text:
struct Batch{S, A<:AbstractVector, M} <: TokenStages
    x::A
    meta::M
end

Batch{S}(x, meta = nothing) where S = Batch{S, typeof(x), typeof(meta)}(x, meta)

setmeta(x::Batch{S}, meta) where S = Batch{S}(x.x, meta)
setvalue(x::Batch{S}, y) where S = Batch{S}(y, x.meta)

function Base.show(io::IO, x::Batch{S}) where S
    print(io, "Batch{", S, "}(", x.x)
    isnothing(x.meta) || print(io, ", ", x.meta)
    print(io, ')')
end

splittability(::BaseTokenization, ::Batch) = Splittable()

splitting(::BaseTokenization, s::Batch) = s.x

wrap(::BaseTokenization, b::Batch{S},  x) where S = S(x, getmeta(b))

# nested

tokenize(tkr::NestedTokenizer, p::ParentStages, t::AbstractTokenization, x::Batch{Document}) = tokenize_procedure!(push!, Vector{Vector{TokenStage}}[], tkr, p, t, x)

tokenize(tkr::NestedTokenizer, p::ParentStages, t::AbstractTokenization, x::Batch{Sentence}) = tokenize_procedure!(push!, Vector{TokenStage}[], tkr, p, t, x)

# indexed

splitting(p::ParentStages, t::IndexedTokenization, b::Batch{Sentence}, x) = enumerate(splitting(p, t.base, b, x))

wrap(p::ParentStages, t::IndexedTokenization, b::Batch{Sentence}, (i, x)) = updatemeta(wrap(p, t.base, b, x), (sentence_id = i,))

splitting(p::ParentStages, t::IndexedTokenization, b::Batch{Document}, x) = enumerate(splitting(p, t.base, b, x))
wrap(p::ParentStages, t::IndexedTokenization, b::Batch{Document}, (i, x)) = updatemeta(wrap(p, t.base, b, x), (document_id = i,))

ext: jl


text:
## -------  CES in A and [h, l, a]  OR  A and a.
# Elasticity = 1 / (1 - cesCoeff)
# Complementarity requires cesCoeff < 0

# How does the CES aggregator work?
abstract type AbstractHCesAggr end
# Aggregate A_q and ability
struct hCesAggrA <: AbstractHCesAggr end
# Aggregate A_q and (a, h, l)
struct hCesAggrAhl <: AbstractHCesAggr end

description(::Type{hCesAggrA}) = "CES aggregator of TFP and ability";
description(::Type{hCesAggrAhl}) =
    "CES aggregator of TFP and (ability, h, study time)";
description(h :: T) where T <: AbstractHCesAggr = description(T);

"""
    $(SIGNATURES)

For a single college (not a ModelObject)
"""
mutable struct HcProdCES{T <: AbstractHCesAggr} <: AbstractHcProdFct
    # College specific tfp
    tfp  ::  Double
    # Neutral tfp (premultiplies everything)
    tfpNeutral :: Double
    timeExp  ::  Double
    hExp  ::  Double
    # Depreciation rate
    deltaH :: Double
    # Ability scale
    aScale  ::  Double
    # Elasticity of substitution
    substElast :: Double
    # Fixed time cost per course
    timePerCourse :: Double
    # Study time per course minimum (this is assigned when study time very low)
    minTimePerCourse :: Double
end

# -----  All colleges

Base.@kwdef mutable struct HcProdCesSwitches{T <: AbstractHCesAggr} <: AbstractHcProdSwitches

    # Same exponents on time and h. If yes, ignore `hExp`.
    sameExponents :: Bool = true
    timeExp :: Double = 0.6
    calTimeExp :: Bool = true
    hExp :: Double = 0.6
    hExpLb :: Double = -1.0
    calHExp :: Bool = false
    deltaH :: Double = 0.0
    calDeltaH :: Bool = false
    aScale :: Double = 0.2
    calAScale :: Bool = true
    substElast :: Double = 0.8
    calSubstElast :: Bool = true
    # TFP for first college
    tfp0 :: Double = 1.0
    tfp0Lb :: Double = 0.02
    tfp0Ub :: Double = 5.0
end

"""
    $(SIGNATURES)

ModelObject that specifies human capital production functions for all colleges. Ben-Porath functional form.

Since all colleges share some parameters, we need a model object that keeps
track of parameters that are common or differ by college.
"""
mutable struct HcProdCesSet{T <: AbstractHCesAggr} <: AbstractHcProdSet
    objId :: ObjectId
    switches :: HcProdCesSwitches{T}
    nc :: CollInt

    # Calibrated parameters
    tfpV  ::  IncreasingVector
    tfpNeutral :: Double
    timeExp  ::  Double
    hExp  ::  Double
    deltaH :: Double
    # Ability scale
    aScale  ::  Double
    substElast :: Double

    # Fixed time cost per course
    timePerCourse :: Double
    # Study time per course minimum (this is assigned when study time very low)
    minTimePerCourse :: Double

    pvec :: ParamVector
end

## -------------- CES: For one college
# For numerical reasons, the substitution elasticity is bounded away from 1.
# If it is set close to 1, the value is simply changed away from 1 when computing dh.

function make_test_hprod_ces(T)
    deltaH = 0.05;
    tfpNeutral = 0.2;
    h = HcProdCES{T}(1.2, tfpNeutral, 0.3, 0.45,
        deltaH, 1.2, 0.9, 0.02, 0.005);
    @assert validate_hprod(h)
    return h
end

function validate_hprod(hS :: HcProdCES{T}) where T
    isValid = (tfp(hS) > 0.01)  &&  (0.0 < time_exp(hS) ≤ 1.0);
    if (subst_elast(hS) < 0.1)  ||  (subst_elast(hS) > 10.0)
        @warn "Substitution elasticity out of bounds: $(subst_elast(hS))";
        isValid = false;
    end
    return isValid
end

tfp_neutral(h :: HcProdCES{T}) where T = h.tfpNeutral;
tfp_neutral(h :: HcProdCesSet{T}) where T = h.tfpNeutral;
subst_elast(h :: HcProdCES{T}) where T = h.substElast;
subst_elast(h :: HcProdCesSet{T}) where T = h.substElast;
ces_coeff(h :: HcProdCES{T}) where T = 1.0 - 1.0 / h.substElast;

# Ces coefficient, bounded away from 0 for computational reasons.
function ces_coeff_bounded(h :: HcProdCES{T}) where T
    cc = ces_coeff(h);
    if (cc >= 0.0)  &&  (cc < 0.05)
        cc = 0.05;
    elseif (cc <= 0.0)  &&  (cc > -0.05)
        cc = -0.05;
    end
    return cc
end

"""
    $(SIGNATURES)

H produced (before shock is realized). Nonnegative.

# Arguments
- nTriedV
    number of courses attempted this period.
- h0V
    h endowments, so that `hV - h0V` is learning.
"""
function dh(hS :: HcProdCES{hCesAggrAhl}, abilV,  hV, h0V, timeV, nTriedV)
    sTimeV = study_time_per_course(hS, timeV, nTriedV);
    xV = (hV .^ h_exp(hS)) .* (sTimeV .^ hS.timeExp) .* exp.(hS.aScale .* abilV);
    pCES = ces_coeff_bounded(hS);
    # The 0.5 ensures that this becomes Cobb Douglas as the elasticity -> 1
    # Numerically it is probably best to pull out tfp
    dhV = nTriedV .* tfp_neutral(hS) .* tfp(hS) .*
        (0.5 .* ((xV ./ tfp(hS)) .^ pCES) .+ 0.5) .^ (1.0 / pCES);
    return dhV
end

# Same with different aggregator
function dh(hS :: HcProdCES{hCesAggrA}, abilV,  hV, h0V, timeV, nTriedV)
    sTimeV = study_time_per_course(hS, timeV, nTriedV);
    xV = (hV .^ h_exp(hS)) .* (sTimeV .^ hS.timeExp);
    pCES = ces_coeff_bounded(hS);
    dhV = nTriedV .* tfp_neutral(hS) .* xV .* tfp(hS) .*
        (0.5 .* ((exp.(hS.aScale .* abilV) ./ tfp(hS)) .^ pCES) .+
        0.5) .^ (1.0 / pCES);
    return dhV
end

show_string(hS :: HcProdCES) =
    "CES with elasticity $(subst_elast(hS))";

Base.show(io :: IO, hS :: HcProdCES) =
    print(io,  "H prod fct:  $(show_string(hS))");

## ---------------------  For all colleges

# Also used for `describe`.
function settings_table(h :: HcProdCesSwitches{T}) where T
    same_exponents(h)  ?  expStr = "the same"  :  expStr = "different";
    ddh = delta_h(h);
    cal_delta_h(h)  ?  deprecStr = "fixed at $ddh"  :  deprecStr = "calibrated";
    return [
        "H production function" "CES";
        "Productivity"  description(T);
        "Exponents on time and h"  expStr;
        "Depreciation"  deprecStr
    ]
end

function StructLH.describe(h :: HcProdCesSet{T}) where T
    tfpString = format_vector("",  [tfp(h, ic)  for  ic = 1 : n_colleges(h)], 1);
    return [
        "H production function"  "CES";
        "Productivity"  description(T);
        "TFP by college"  tfpString
    ]
end

function settings_list(h :: HcProdCesSwitches{T}, st) where T
    eqnHChange = ["H production function", "eqnHChange",  eqn_hchange(h)];
    return [eqnHChange]
end

function validate_hprod(s :: HcProdCesSwitches)
    isValid = true;
    if s.sameExponents && s.calHExp
        @warn "Cannot calibrate hExp if it has the same value as timeExp"
        isValid = false;
    end
    return isValid
end

function validate_hprod_set(h :: HcProdCesSet)
    isValid = (h.nc > 1)  &&  (h.timeExp > 0.0)  &&
        (h.aScale > 0.0)  &&  (h.timePerCourse > 0.0);
    isValid = isValid  &&  (1.0 > delta_h(h) >= 0.0);
    return isValid
end

function eqn_hchange(h :: HcProdCesSwitches{hCesAggrAhl})
    "\\hTfpNeutral (\\hTfp ^ \\hSeCoeff + [\\sTimePerCourse^{\\hTimeExp}  \\hc^{\\hHExp} e^{\\hAScale \\abil}] ^ \\hSeCoeff) ^ (1/\\hSeCoeff)"
end

function eqn_hchange(h :: HcProdCesSwitches{hCesAggrA})
    "\\hTfpNeutral \\sTimePerCourse^{\\hTimeExp}  \\hc^{\\hHExp} (\\hTfp ^ \\hSeCoeff + [e^{\\hAScale \\abil}] ^ \\hSeCoeff) ^ (1/\\hSeCoeff)"
end

# Initialize with defaults
function make_hc_prod_set(objId :: ObjectId, nc :: Integer,
    switches :: HcProdCesSwitches{T}) where T

    @assert validate_hprod(switches);

    pTimePerCourse = init_time_per_course();
    pTfp = init_tfp(objId, nc;
        tfp0 = switches.tfp0, tfp0Lb = switches.tfp0Lb, tfp0Ub = switches.tfp0Ub);

    timeExp = switches.timeExp;
    pTimeExp = Param(:timeExp,  ldescription(:hTimeExp), lsymbol(:hTimeExp),
        timeExp, timeExp, 0.2, 0.9, switches.calTimeExp);

    tfpNeutral = 0.5;
    pTfpNeutral = Param(:tfpNeutral,
         ldescription(:hTfpNeutral), lsymbol(:hTfpNeutral),
        tfpNeutral, tfpNeutral, 0.02, 2.0, true);

    pDeltaH = init_h_depreciation(switches);

    if switches.sameExponents
        hExp = switches.timeExp;
        calHExp = false;
    else
        hExp = switches.hExp;
        calHExp = switches.calHExp;
    end
    # Allowing negative hExp for concavity.
    pHExp = Param(:hExp,  ldescription(:hHExp), lsymbol(:hHExp),
        hExp, hExp, switches.hExpLb, 0.9, calHExp);

    aScale = switches.aScale;
    pAScale = Param(:aScale,  ldescription(:hAScale), lsymbol(:hAScale),
        aScale, aScale, 0.02, 2.0, switches.calAScale);

    # Elasticity is 1/(1 - seCoeff).
    # Elasticity of 0.5 is pretty close to Leontief. Avoid.
    substElast = switches.substElast;
    pSe = Param(:substElast,  ldescription(:hSubstElast), lsymbol(:hSubstElast),
        substElast, substElast, 0.5, 4.0, switches.calSubstElast);

    pvec = ParamVector(objId,
        [pTfpNeutral, pTimeExp, pHExp, pDeltaH, pAScale, pSe, pTimePerCourse]);
    # Min study time required per course. Should never bind.
    minTimePerCourse =
        hours_per_week_to_mtu(0.1 / data_to_model_courses(1));

    return HcProdCesSet(objId, switches, CollInt(nc), pTfp, tfpNeutral,
        timeExp, hExp, ModelParams.value(pDeltaH), aScale,
        substElast, pTimePerCourse.value, minTimePerCourse, pvec)
end

make_test_hc_ces_set(T) =
    make_hc_prod_set(ObjectId(:HProd), 4, HcProdCesSwitches{T}(deltaH = 0.05));

# Make h production function for one college
function make_h_prod(hs :: HcProdCesSet{T}, iCollege :: Integer) where T
    return HcProdCES{T}(tfp(hs, iCollege), tfp_neutral(hs),
        time_exp(hs), h_exp(hs),
        delta_h(hs), hs.aScale, subst_elast(hs),
        hs.timePerCourse,  hs.minTimePerCourse);
end

# ----------------


text:
using BSplineKit
using Test

# TODO
# - test recombined basis (test function should satisfy BCs...)

# We make sure that the function is specialised on `f`, to fully avoid
# allocations.
# See https://docs.julialang.org/en/v1/manual/performance-tips/#Be-aware-of-when-Julia-avoids-specializing
function check_approximate_in_place!(f::F, fapprox) where {F}
    S = copy(spline(fapprox))  # assumed to already approximate `f`
    @test 0 == @allocated approximate!(f, fapprox)  # no allocations!
    @test spline(fapprox) == S  # check that nothing changed
    nothing
end

function test_approximation(B)
    k = order(B)

    ftest = sin

    f_fast = @inferred approximate(ftest, B, VariationDiminishing())
    @test startswith(repr(f_fast), "SplineApproximation containing")
    check_approximate_in_place!(ftest, f_fast)

    # With these two methods, the approximation should exactly match a
    # polynomial of degree k - 1.
    f_interp = @inferred approximate(ftest, B, ApproxByInterpolation(B))
    @test startswith(repr(f_interp), "SplineApproximation containing")
    check_approximate_in_place!(ftest, f_interp)

    f_opt = @inferred approximate(ftest, B, MinimiseL2Error())
    @test startswith(repr(f_opt), "SplineApproximation containing")
    check_approximate_in_place!(ftest, f_opt)

    # Polynomial of degree k - 1.
    # The interpolation and minimisation methods should exactly approximate such
    # a polynomial.
    @testset "Polynomial" begin
        poly_coefs = ntuple(d -> (-d)^d, Val(k))
        fpoly(x) = evalpoly(x, poly_coefs)

        xfine = range(boundaries(B)...; length = 3 * length(B))

        approximate!(fpoly, f_interp)
        @test fpoly.(xfine) ≈ f_interp.(xfine)

        approximate!(fpoly, f_opt)
        @test fpoly.(xfine) ≈ f_opt.(xfine)
    end

    nothing
end

@testset "Approximation" begin
    N = 16
    breaks = [-cos(π * n / N) for n = 0:N]
    @testset "Order $k" for k = 4:5
        B = BSplineBasis(BSplineOrder(k), copy(breaks))
        # R = RecombinedBSplineBasis(B, Derivative(1))  # Neumann BCs
        test_approximation(B)
    end
end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
using Cairo, Gtk, Colors

struct AtariCairoCtx <: AbstractCtx end
struct AtariRGBCtx <: AbstractCtx end
struct NoCtx <: AbstractCtx end
struct AtariGtkCtx <: AbstractCtx
    canvas::GtkCanvas
    win::GtkWindowLeaf
end

function Ctx(env::AtariEnv, mode::Symbol=:no_render)
    if mode == :human_pane
        AtariCairoCtx()

    elseif mode == :human_window
        canvas = @GtkCanvas()
        win = GtkWindow(canvas, "AtariEnv")
        show(canvas)
        visible(win, false)
        signal_connect(win, "delete-event") do widget, event
            ccall((:gtk_widget_hide_on_delete, Gtk.libgtk), Bool, (Ptr{GObject},), win)
        end

        AtariGtkCtx(canvas, win)
    elseif mode == :rgb

        AtariRGBCtx()
    elseif mode == :no_render
        NoCtx()
    else
        error("Unrecognized mode in Ctx(): $(mode)")
    end
end

function drawcanvas(env::AtariEnv)
    screen_grab = getScreenRGB(env.ale)
    w, h = getScreenWidth(env.ale), getScreenHeight(env.ale)

    # Converting RGB values to 32 bit color...
    r_screen_grab = screen_grab[1:3:end] * 0x00010000
    g_screen_grab = screen_grab[2:3:end] * 0x00000100
    b_screen_grab = screen_grab[3:3:end]

    rgb_array = reshape(r_screen_grab .+ b_screen_grab .+ g_screen_grab, w, h) |> transpose |> Array
end

render!(env::AbstractEnv, ctx::NoCtx) = nothing
render!(env::AtariEnv, ctx::AtariCairoCtx) = env |> drawcanvas |> CairoRGBSurface |> display

function render!(env::AtariEnv, ctx::AtariGtkCtx)
    !visible(ctx.win) && visible(ctx.win, true)
    @guarded draw(ctx.canvas) do widget
        context = getgc(ctx.canvas)
        image(context, CairoRGBSurface(drawcanvas(env)),
            0, 0, getScreenWidth(env.ale), getScreenHeight(env.ale))
    end
end

render!(env::AtariEnv, ctx::AtariRGBCtx) = get_preprocessed_RGB(env)


text:
#  Copyright 2021, Author: Ali Daraeepour (a.daraeepour@duke.edu)
#                  Contributors: Mauricio Hernandez (mmh54@duke.edu)
#  This Source Code Form is subject to the terms of the MIT  License. If a copy of the MIT was
# not distributed with this file, You can obtain one at https://mit-license.org/.
#############################################################################
# GRACE BAU UC
# This program solves the Business as Usual Unit Commitment problem of
# Duke Energy Power Sytem
# See https://github.com/power-systems-op/grace-arpae/
#############################################################################

"""
File: BAU_OPM.jl
Version: 7.7
...
# Arguments: None
# Outputs: dataframe
# Examples: N/A
"""

#NOTES: This file was originally labeled as BAU_OPM_V5_Nuc_Cogen_Imp_Exp_MustRun_DCC
# Changes in this version: Nuclear data was included
#= Changes
1. Specify data type of the vectors that store the data read from CSV input files.
2. Replace several global variables by local variables.
3. Replace several auxiliary variables by data structures.
4. Replace the variables that store the results of the first UC by data structures.
5. Divide the main script in different Julia scripts:
 5.1 constants.jl: Store the constants used across Julia files
 5.2 data_structure.jl: Define the data structures used across Julia files
 5.3 import_data.jl: Read input files and store data in dataframes and vectors
 5.4 fucr_model.jl: Solve first unit commitment model
 5.5 BAU_OPM.jl: Main program
 5.6 bucr_model.jl: Solve Balancing UC problem
 5.6 Embed main program in a function
=#

#=
---------------------------
Suggestions to test the code:
- Don't force the nature of the variables.
- Compare results, objective functions. For 2 days.
- If values are not the same.

- Objective function is similar or smaller
- Cost, if they are higher.
- Optimality gap could.

- Run the code v7 in different computers.
=#

#using Queryverse
using CPLEX
using CpuId
using CSV
using DataFrames
using Dates
using DelimitedFiles
using JuMP
using Logging
using Parameters

t1 = time_ns()
include("constants.jl")
include("data_structure.jl")
t1_read_data = time_ns()
include("import_data.jl")
t2_read_data = time_ns()
time_read_data = (t2_read_data -t1_read_data)/1.0e9;

include("fucr_model.jl")
include("bucr1_model.jl")
include("sucr_model.jl")
include("bucr2_model.jl")
## Enabling debugging code, use ENV["JULIA_DEBUG"] = "" to desable ging code
ENV["JULIA_DEBUG"] = "all"

# Log file
io_log = open( string(".//outputs//logs//UC_BAU_",
        Dates.format(now(), "yyyy-mm-dd_HH-MM-SS"), ".txt"), FILE_ACCESS_APPEND)

logger = SimpleLogger(io_log)
flush(io_log)
global_logger(logger)

@info "Hardware Features: " cpuinfo()

write(io_log, "Running model from day $INITIAL_DAY to day $FINAL_DAY with the following parameters:\n")
write(io_log, "Load-Shedding Penalty: $LOAD_SHED_PENALTY, Over-generation Penalty: $OVERGEN_PENALTY\n")
write(io_log, "Max Load-Shedding Penalty $LOAD_SHED_MAX, Max Over-generation Penalty: $OVERGEN_MAX\n")
write(io_log, "MaxGenLimit Viol Penalty: $VIOLATION_PENALTY, OptimalityGap: $OVERGEN_PENALTY\n")

time_performance_header    = ["Section", "Time", "Note1", "Note2", "Note3", "Note4"]
open(".//outputs//time_performance.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(time_performance_header), ',')
end; # closes file

@info "Time to read input data (s): $time_read_data";
open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
        writedlm(io, hcat("Read Input Data", time_read_data, "",
                "", "", "Read CSV files"), ',')
end;

objective_values_header    = ["Section", "Time", "Time2", "Note1", "Value"]
open(".//outputs//objective_values_v76.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(objective_values_header), ',')
end; # closes file

## Headers of output files
@time begin
FUCR_GenOutputs_header      = ["Day", "Hour", "GeneratorID", "UNIT_NAME", "MinPowerOut", "MaxPowerOut", "Output", "On/off", "ShutDown", "Startup", "UpSpinRes", "Non_SpinRes", "DownSpinRes", "TotalGenPosV", "TotalGenNegV", "MaxGenVio", "MinGenVio"]
FUCR_PeakerOutputs_header   = ["Day", "Hour", "PeakerID", "UNIT_NAME", "MinPowerOut", "MaxPowerOut", "Output", "On/off", "ShutDown", "Startup", "UpSpinRes", "Non_SpinRes", "DownSpinRes", "TotalGenPosV", "TotalGenNegV", "MaxGenVio", "MinGenVio"]
FUCR_StorageOutputs_header  = ["Day", "Hour", "StorageUniID", "Power", "EnergyLimit", "Charge_St", "Discharge_St", "Idle_St", "storgChrgPwr", "storgDiscPwr", "storgSOC", "storgResUp", "storgResDn"]
FUCR_TranFlowOutputs_header = ["Day", "Time period", "Source", "Sink", "Flow", "TransCap"]
FUCR_Curtail_header         = ["Day", "Hour", "Zone", "OverGeneration", "DemandCurtailment", "WindCrtailment", "SolarCurtailment", "HydroSpillage"]
SUCR_GenOutputs_header      = ["Day", "Hour", "GeneratorID", "UNIT_NAME", "MinPowerOut", "MaxPowerOut", "Output", "On/off", "ShutDown", "Startup", "UpSpinRes", "Non_SpinRes", "DownSpinRes", "TotalGenPosV", "TotalGenNegV", "MaxGenVio", "MinGenVio"]
SUCR_PeakerOutputs_header   = ["Day", "Hour", "GeneratorID", "UNIT_NAME", "MinPowerOut", "MaxPowerOut", "Output", "On/off", "ShutDown", "Startup", "UpSpinRes", "Non_SpinRes", "DownSpinRes", "TotalGenPosV", "TotalGenNegV", "MaxGenVio", "MinGenVio"]
SUCR_StorageOutputs_header  = ["Day", "Hour", "StorageUniID", "Power", "EnergyLimit", "Charge_St", "Discharge_St", "Idle_St", "storgChrgPwr", "storgDiscPwr", "storgSOC", "storgResUp", "storgResDn"]
SUCR_TranFlowOutputs_header = ["Day", "Time period", "Source", "Sink", "Flow", "TransCap"]
SUCR_Curtail_header         = ["Day", "Hour", "Zone", "OverGeneration", "DemandCurtailment", "WindCrtailment", "SolarCurtailment", "HydroSpillage"]
BUCR_GenOutputs_header      = ["Day", "Hour", "GeneratorID", "UNIT_NAME", "MinPowerOut", "MaxPowerOut", "Output", "On/off", "ShutDown", "Startup", "TotalGenPosV", "TotalGenNegV", "MaxGenVio", "MinGenVio"]
BUCR_PeakerOutputs_header   = ["Day", "Hour", "GeneratorID", "UNIT_NAME", "MinPowerOut", "MaxPowerOut", "Output", "On/off", "ShutDown", "Startup", "TotalGenPosV", "TotalGenNegV", "MaxGenVio", "MinGenVio"]
BUCR_StorageOutputs_header  = ["Day", "Hour", "StorageUniID", "Power", "EnergyLimit", "Charge_St", "Discharge_St", "Idle_St", "storgChrgPwr", "storgDiscPwr", "storgSOC", "storgResUp", "storgResDn"]
BUCR_TranFlowOutputs_header = ["Day", "Time period", "Source", "Sink", "Flow", "TransCap"]
BUCR_Curtail_header         = ["Day", "Hour", "Zone", "OverGeneration", "DemandCurtailment", "WindCrtailment", "SolarCurtailment", "HydroSpillage"]

end #timer

## Spreadsheets for the first unit commitment run
# Creating conventional generating units' schedules in the 1st unit commitment run
open(".//outputs//FUCR_GenOutputs.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(FUCR_GenOutputs_header), ',')
end; # closes file

open(".//outputs//FUCR_PeakerOutputs.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(FUCR_PeakerOutputs_header), ',')
end; # closes file

open(".//outputs//FUCR_StorageOutputs.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(FUCR_StorageOutputs_header), ',')
end;

open(".//outputs//FUCR_TranFlowOutputs.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(FUCR_TranFlowOutputs_header), ',')
end; # closes file

open(".//outputs//FUCR_Curtail.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(FUCR_Curtail_header), ',')
end; # closes file

# Spreadsheets for the second unit commitment run
# Creating conventional generating units' schedules in the 2nd unit commitment run
open(".//outputs//SUCR_GenOutputs.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(SUCR_GenOutputs_header), ',')
end; # closes file

open(".//outputs//SUCR_PeakerOutputs.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(SUCR_PeakerOutputs_header), ',')
end; # closes file

open(".//outputs//SUCR_StorageOutputs.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(SUCR_StorageOutputs_header), ',')
end; # closes file

open(".//outputs//SUCR_TranFlowOutputs.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(SUCR_TranFlowOutputs_header), ',')
end; # closes file

open(".//outputs//SUCR_Curtail.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(SUCR_Curtail_header), ',')
end; # closes file

# Spreadsheets for the balancing unit commitment run
# Write the conventional generators' schedules
open(".//outputs//BUCR_GenOutputs.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(BUCR_GenOutputs_header), ',')
end; # closes file

open(".//outputs//BUCR_PeakerOutputs.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(BUCR_PeakerOutputs_header), ',')
end; # closes file

# Writing storage units' optimal schedules into CSV file
open(".//outputs//BUCR_StorageOutputs.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(BUCR_StorageOutputs_header), ',')
end; # closes file

# Writing the transmission flow schedules in CSV file
open(".//outputs//BUCR_TranFlowOutputs.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(BUCR_TranFlowOutputs_header), ',')
end; # closes file

open(".//outputs//BUCR_Curtail.csv", FILE_ACCESS_OVER) do io
    writedlm(io, permutedims(BUCR_Curtail_header), ',')
end; # closes file

function run_main_code()
	# Setting initial values
	fucr_gens = GensResults(
		onoff_init = copy(DF_Generators.StatusInit),
		power_out_init = copy(DF_Generators.PowerInit),
		uptime_init = copy(DF_Generators.UpTimeInit),
		dntime_init = copy(DF_Generators.DownTimeInit))

	fucr_peakers = PeakersResults(
		onoff_init = copy(DF_Peakers.StatusInit),
		power_out_init = copy(DF_Peakers.PowerInit),
		uptime_init = copy(DF_Peakers.UpTimeInit),
		dntime_init = copy(DF_Peakers.DownTimeInit))

	fucr_storg = StorageResults(soc_init =copy(DF_Storage.SOCInit))

	FUCR = UC_Results(fucr_gens, fucr_peakers, fucr_storg)

	sucr_gens = GensResults(
		onoff_init = copy(DF_Generators.StatusInit),
		power_out_init = copy(DF_Generators.PowerInit),
		uptime_init = copy(DF_Generators.UpTimeInit),
		dntime_init = copy(DF_Generators.DownTimeInit))

	sucr_peakers = PeakersResults(
		onoff_init = copy(DF_Peakers.StatusInit),
		power_out_init = copy(DF_Peakers.PowerInit),
		uptime_init = copy(DF_Peakers.UpTimeInit),
		dntime_init = copy(DF_Peakers.DownTimeInit),
		genout_block = zeros(Float64, PEAKERS,BLOCKS,24-INIT_HR_SUCR+INIT_HR_FUCR))

	sucr_storg = StorageResults(soc_init =copy(DF_Storage.SOCInit))

	SUCR = UC_Results(sucr_gens, sucr_peakers, sucr_storg)

	bucr1_gens = GensResults(
		onoff_init = copy(DF_Generators.StatusInit),
		power_out_init = copy(DF_Generators.PowerInit),
		uptime_init = copy(DF_Generators.UpTimeInit),
		dntime_init = copy(DF_Generators.DownTimeInit))

	bucr1_peakers = PeakersResults(
		onoff_init = copy(DF_Peakers.StatusInit),
		power_out_init = copy(DF_Peakers.PowerInit),
		uptime_init = copy(DF_Peakers.UpTimeInit),
		dntime_init = copy(DF_Peakers.DownTimeInit))

	bucr1_storg = StorageResults(soc_init =copy(DF_Storage.SOCInit))

	BUCR1 = UC_Results(bucr1_gens, bucr1_peakers, bucr1_storg)

	bucr2_gens = GensResults()
	bucr2_peakers = PeakersResults()
	bucr2_storg = StorageResults(soc_init = zeros(STORG_UNITS))

	BUCR2 = UC_Results(bucr2_gens, bucr2_peakers, bucr2_storg)

	fucr_to_bucr1_gens = GensResults()
	fucr_to_bucr1_peakers = PeakersResults()
	fucr_to_bucr1_storg = StorageResults()

	FUCRtoBUCR1 = UC_Results(fucr_to_bucr1_gens, fucr_to_bucr1_peakers,
						fucr_to_bucr1_storg)

	# Creating variables that transfer optimal schedules between the Models
	#TODO: Some of the below variables may be unneccsary and can be deleted. Check at the end

	sucr_to_bucr2_gens = GensResults()
	sucr_to_bucr2_peakers = PeakersResults()
	sucr_to_bucr2_storg = StorageResults()

	SUCRtoBUCR2 = UC_Results(sucr_to_bucr2_gens, sucr_to_bucr2_peakers,
						sucr_to_bucr2_storg)

	# Auxiliary variables for enforcing commitment of slow and fast-start units
	# in BUCRs
	#TODO: Include these variables in UC_Results structures
	BUCR1_Commit_LB = zeros(Int64, GENS)
	BUCR1_Commit_UB = zeros(Int64, GENS)
	BUCR2_Commit_LB = zeros(Int64, GENS)
	BUCR2_Commit_UB = zeros(Int64, GENS)

	BUCR1_Commit_Peaker_LB = zeros(Int64, PEAKERS)
	BUCR1_Commit_Peaker_UB = zeros(Int64, PEAKERS)
	BUCR2_Commit_Peaker_LB = zeros(Int64, PEAKERS)
	BUCR2_Commit_Peaker_UB = zeros(Int64, PEAKERS)

	## Pre-processing the data to calculate the model inputs:
	# The time range lower-bound for min up constraint using the alternative approach

	# LB for slow-start conventional generators
	lb_MUT = zeros(Int64, GENS, HRS_FUCR);
	lb_MDT = zeros(Int64, GENS, HRS_FUCR);
	for g in 1:GENS, t in 1:HRS_FUCR
	    lb_MUT[g,t] =  replace_lower_limit(t-DF_Generators.MinUpTime[g]+1);
	    lb_MDT[g,t] =  replace_lower_limit(t-DF_Generators.MinDownTime[g]+1);
	end;

	# LB for peakers
	lb_MUT_Peaker=zeros(Int64, PEAKERS, HRS_FUCR)
	lb_MDT_Peaker=zeros(Int64, PEAKERS, HRS_FUCR)
	for k in 1:PEAKERS , t in 1:HRS_FUCR
	    lb_MUT_Peaker[k,t] =  replace_lower_limit(t-DF_Peakers.MinUpTime[k]+1);
	    lb_MDT_Peaker[k,t] =  replace_lower_limit(t-DF_Peakers.MinDownTime[k]+1);
	end;
	## TODO: We need two more loops for calculating the lb_MUT and lb_MDT for the second UC runs
	#************************************************************************************
	## The foor loop runs two WAUC models and the RTUC models every day
	for day = INITIAL_DAY:FINAL_DAY
	   t1_day_execution = time_ns()

	   # Bottom and upper cells of the demand data needed for running the first
	   # WA-UC run at 6 am with 7-day look-ahead horizon
	   D_Rng_Dn_FUCR = ((day-1)*(INIT_HR_FUCR+HRS_FUCR)) + INIT_HR_FUCR + 1;
	   D_Rng_Up_FUCR = day*(INIT_HR_FUCR+HRS_FUCR);
	   #rngdn_fucr
	   R_Rng_Dn_FUCR = ((day-1)*(24))+INIT_HR_FUCR+1
	   R_Rng_Up_FUCR = ((day+6)*(24))

	   # Demand Data Pre-Processing for FUCR and SUCR
	   fucr_prepoc = DemandWAPreprocGens(
	   		copy(FUCR_Demands[D_Rng_Dn_FUCR:D_Rng_Up_FUCR, :]),
	    	copy(FUCR_SolarGs[R_Rng_Dn_FUCR:R_Rng_Up_FUCR, :]),
	        copy(FUCR_WindGs[R_Rng_Dn_FUCR:R_Rng_Up_FUCR, :]),
	        copy(FUCR_HydroGs[R_Rng_Dn_FUCR:R_Rng_Up_FUCR, :]),
	        copy(FUCR_NuclearGs[R_Rng_Dn_FUCR:R_Rng_Up_FUCR, :]),
	        copy(FUCR_CogenGs[R_Rng_Dn_FUCR:R_Rng_Up_FUCR, :]) )

	    # Bottom and upper cells of the demand data needed for running the second
	    # WA-UC run at 5 pm with 7-day look-ahead horizon
	    D_Rng_Dn_SUCR = ((day-1)*(INIT_HR_FUCR+HRS_FUCR))+INIT_HR_SUCR+1;
	    D_Rng_Up_SUCR = day*(INIT_HR_FUCR+HRS_FUCR);

	    R_Rng_Dn_SUCR = ((day-1)*(24))+INIT_HR_SUCR+1;
	    R_Rng_Up_SUCR = ((day+6)*(24));

	    sucr_prepoc = DemandWAPreprocGens(
	    	copy(SUCR_Demands[D_Rng_Dn_SUCR:D_Rng_Up_SUCR, :]),
	    	copy(SUCR_SolarGs[R_Rng_Dn_SUCR:R_Rng_Up_SUCR, :]),
	    	copy(SUCR_WindGs[R_Rng_Dn_SUCR:R_Rng_Up_SUCR, :]),
	    	copy(SUCR_HydroGs[R_Rng_Dn_SUCR:R_Rng_Up_SUCR, :]),
	    	copy(SUCR_NuclearGs[R_Rng_Dn_SUCR:R_Rng_Up_SUCR, :]),
	    	copy(SUCR_CogenGs[R_Rng_Dn_SUCR:R_Rng_Up_SUCR, :]) )

		# This block models the first UC optimization that is run in the morning
	    t1_FUCRmodel = time_ns()

		FUCRtoBUCR1 = fucr_model(day, DF_Generators, DF_Peakers, DF_Storage,
			FuelPrice, FuelPricePeakers, Map_Gens, Map_Peakers, Map_Storage,
            FUCR, fucr_prepoc, lb_MDT, lb_MUT, lb_MDT_Peaker, lb_MUT_Peaker)

		# This block models the Balancing Unit Commitment Runs between the
		# morning and evening UC Runs
	    for h=1:INIT_HR_SUCR-INIT_HR_FUCR # number of BUCR periods between FUCR and SUCR
	        # Pre-processing demand variables
	        t1_BUCR_SUCR_data_hand = time_ns()
			# Bottom cell of the demand data needed for running the first WA-UC
			# run at 6 am with 7-day look-ahead horizon
	        D_Rng_BUCR1 = ((day-1)*24)+INIT_HR_FUCR+h

			bucr1_prepoc = DemandHrPreprocGens(
				copy(BUCR_Demands[D_Rng_BUCR1, :]),
				copy(BUCR_SolarGs[D_Rng_BUCR1, :]),
				copy(BUCR_WindGs[D_Rng_BUCR1, :]),
				copy(BUCR_HydroGs[D_Rng_BUCR1, :]),
				copy(BUCR_NuclearGs[D_Rng_BUCR1, :]),
				copy(BUCR_CogenGs[D_Rng_BUCR1, :]) )

	        # Preprocessing module that fixes the commitment of slow-start units
			# to their FUCR's outcome and determines the binary commitment bounds
			# for fast-start units dependent to their initial up/down time and
			# min up/down time limits if the units are slow their BAUC's
			# commitment is fixed to their FUCR's schedule

	        for g=1:GENS
				BUCR1_Commit_LB[g] = fix_uc_slow_start(FUCRtoBUCR1.gens.onoff[g,h])
				BUCR1_Commit_UB[g] = fix_uc_slow_start(FUCRtoBUCR1.gens.onoff[g,h])
	        end;

	        # If the units are fast, their BAUC's commitment could be fixed
	        # to 0 or 1 or vary between 0 or 1 dependent to their initial
			# up/down time and min up/down time
	        for k=1:PEAKERS
	            if BUCR1.peakers.dntime_init[k] == 0
	                if BUCR1.peakers.uptime_init[k] < DF_Peakers.MinUpTime[k]
	                    BUCR1_Commit_Peaker_LB[k] = 1;
	                    BUCR1_Commit_Peaker_UB[k] = 1;
	                else
	                    BUCR1_Commit_Peaker_LB[k] = 0;
	                    BUCR1_Commit_Peaker_UB[k] = 1;
	                end
	            elseif BUCR1.peakers.dntime_init[k] < DF_Peakers.MinDownTime[k]
	                BUCR1_Commit_Peaker_LB[k] = 0;
	                BUCR1_Commit_Peaker_UB[k] = 0;
	            else
	                BUCR1_Commit_Peaker_LB[k] = 0;
	                BUCR1_Commit_Peaker_UB[k] = 1;
	            end
	        end

	        t2_BUCR_SUCR_data_hand = time_ns();

	        time_BUCR_SUCR_data_hand = (t2_BUCR_SUCR_data_hand -t1_BUCR_SUCR_data_hand)/1.0e9;
	        @info "BUCR_SUCR data handling for day $day executed in (s): $time_BUCR_SUCR_data_hand";

	        open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
	                writedlm(io, hcat("BUCR1model", time_BUCR_SUCR_data_hand,
						"day: $day", "hour $(h+INIT_HR_FUCR)",
						"Pre-processing variables", "Data Manipulation"), ',')
	        end; #closes file

			# This function models the 2nd UC optimization that is run in the evening
			BUCR1, SUCR, BUCR2 = bucr1_model(day, h, DF_Generators, DF_Peakers,
				DF_Storage, FuelPrice, FuelPricePeakers, Map_Gens, Map_Peakers,
				Map_Storage, BUCR1, SUCR, BUCR2, FUCRtoBUCR1,
				bucr1_prepoc, BUCR1_Commit_LB, BUCR1_Commit_UB,
				BUCR1_Commit_Peaker_LB, BUCR1_Commit_Peaker_UB)

	    end # ends loop that runs hourly BUCR1 between INIT_HR_FUCR & INIT_HR_SUCR

		# Run second unit commitment model
		SUCRtoBUCR2 = sucr_model(day, DF_Generators, DF_Peakers, DF_Storage,
					FuelPrice, FuelPricePeakers, Map_Gens,
					Map_Storage, SUCR, sucr_prepoc, lb_MDT, lb_MUT,
					lb_MDT_Peaker, lb_MUT_Peaker);

		## Initilization of the next UC Run
		#=
		    # This must be updated later when we run two WAUCs every day and then RTUCs
		        for g=1:GENS
		            DF_Generators.StatusInit[g]=JuMP.value.(SUCR_genOnOff[g,24-INIT_HR_SUCR+INIT_HR_FUCR]);
		            DF_Generators.PowerInit[g]=JuMP.value.(SUCR_genOut[g,24-INIT_HR_SUCR+INIT_HR_FUCR]);
		        end
		        for p=1:STORG_UNITS
		            DF_Storage.SOCInit[p]=JuMP.value.(SUCR_storgSOC[p,24-INIT_HR_SUCR+INIT_HR_FUCR]);
		        end
		=#

		# This block models the second Balancing Unit Commitment Run, which is
		# for the time range between the morning and evening UC Runs
	    for h= 1:(24-INIT_HR_SUCR+INIT_HR_FUCR) # number of BUCR periods between FUCR and SUCR
	        if h==24-INIT_HR_SUCR+INIT_HR_FUCR
	            println("this period is hour: ", h)
	        end

	        # Preprocessing the demand data
			# Cell of the demand data needed for running the first WAUC run
			# at 6 am with 7-day look-ahead horizon
	        D_Rng_BUCR2 = ((day-1)*24)+INIT_HR_SUCR+h

			bucr2_prepoc = DemandHrPreprocGens(
				copy(BUCR_Demands[D_Rng_BUCR2, :]),
				copy(BUCR_SolarGs[D_Rng_BUCR2, :]),
				copy(BUCR_WindGs[D_Rng_BUCR2, :]),
				copy(BUCR_HydroGs[D_Rng_BUCR2, :]),
				copy(BUCR_NuclearGs[D_Rng_BUCR2, :]),
				copy(BUCR_CogenGs[D_Rng_BUCR2, :]) )

			#TODO: Delete this
			#=BUCR2_Hr_Demand = BUCR_Demands[D_Rng_BUCR2, :]
	        BUCR2_Hr_SolarG = BUCR_SolarGs[D_Rng_BUCR2, :]
	        BUCR2_Hr_WindG = BUCR_WindGs[D_Rng_BUCR2, :]
	        BUCR2_Hr_HydroG = BUCR_HydroGs[D_Rng_BUCR2, :]
	        BUCR2_Hr_NuclearG = BUCR_NuclearGs[D_Rng_BUCR2, :]
	        BUCR2_Hr_CogenG = BUCR_CogenGs[D_Rng_BUCR2, :]
			=#

	        # Preprocessing module that fixes the commitment of slow-start units
			# to their FUCR's outcome and determines the binary commitment bounds
			# for fast-start units dependent to their initial up/down time and
			# minimum up/down time limits. If the units are slow their BAUC's
			# commitment is fixed to their FUCR's schedule
	        for g=1:GENS
	            # If DF_Generators.FastStart[g]==0 # if the units are slow their
				# BAUC's commitment is fixed to their FUCR's schedule
				BUCR2_Commit_LB[g] =  fix_uc_slow_start(SUCRtoBUCR2.gens.onoff[g,h]);
				BUCR2_Commit_UB[g] =  fix_uc_slow_start(SUCRtoBUCR2.gens.onoff[g,h]);
	            #= TODO: Delete this
				if SUCRtoBUCR2.gens.onoff[g,h]==0
	                BUCR2_Commit_LB[g] = 0;
	                BUCR2_Commit_UB[g] = 0;
	            else
	                BUCR2_Commit_LB[g] = 1;
	                BUCR2_Commit_UB[g] = 1;
	            end=#
	        end

	        # if the units are fast their BAUC's commitment could be fixed to 0
			# or 1 or vary between 0 or 1 dependent to their initial up/down
			# time and minimum up/down time
	        for k=1:PEAKERS
	            if BUCR2.peakers.dntime_init[k]==0
	                if BUCR2.peakers.uptime_init[k] < DF_Peakers.MinUpTime[k]
	                    BUCR2_Commit_Peaker_LB[k] = 1;
	                    BUCR2_Commit_Peaker_UB[k] = 1;
	                else
	                    BUCR2_Commit_Peaker_LB[k] = 0;
	                    BUCR2_Commit_Peaker_UB[k] = 1;
	                end
	            elseif BUCR2.peakers.dntime_init[k] < DF_Peakers.MinDownTime[k]
	                    BUCR2_Commit_Peaker_LB[k] = 0;
	                    BUCR2_Commit_Peaker_UB[k] = 0;
	            else
	                    BUCR2_Commit_Peaker_LB[k] = 0;
	                    BUCR2_Commit_Peaker_UB[k] = 1;
	            end
	        end

			# Run Second Balancing Unit Commitment Model
			FUCR, BUCR1, BUCR2 = bucr2_model(day, h, DF_Generators, DF_Peakers,
				DF_Storage, FuelPrice, FuelPricePeakers, Map_Gens, Map_Peakers,
				Map_Storage, FUCR, BUCR1, BUCR2, SUCRtoBUCR2, bucr2_prepoc,
				BUCR2_Commit_LB, BUCR2_Commit_UB, BUCR2_Commit_Peaker_LB,
				BUCR2_Commit_Peaker_UB)

		end # loop that runs hourly BUCR between 1 to (24 - INIT_HR_SUCR + INIT_HR_FUCR)

		t2_day_execution = time_ns()
		time_day_execution = (t2_day_execution - t1_day_execution)/1.0e9;

		open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
				writedlm(io, hcat("Whole Day", time_day_execution,
						"day: $day", " ", "Whole Day Execution"), ',')
		end; #closes file

	end # foor loop that runs the UC model on  a daily basis

end; # end of run_main_code

##

# Function to "replace" values equal to zero by 0, 1 otherwise
fix_uc_slow_start(x::Int64) = x == 0 ? 0 : 1

# Function to replace values below 1 by 1
replace_lower_limit(x::Int64) = x < 1 ? 1 : x

run_main_code()

t2 = time_ns()
elapsedTime = (t2 -t1)/1.0e9;

write(io_log, "Whole program time execution (s):\t $elapsedTime\n")
@info "Whole Program setup executed in (s):" elapsedTime;

open(".//outputs//time_performance.csv", FILE_ACCESS_APPEND) do io
        writedlm(io, hcat("Whole Program", elapsedTime,
                "", "", "Whole Execution"), ',')
end; #closes file

close(io_log);

max_stars_repo_name: power-systems-op/grace-arpae

max_issues_repo_name: power-systems-op/grace-arpae

max_forks_repo_name: power-systems-op/grace-arpae


text:
using Plots, Dates

a=[
"21.693180629997933" "
15:09:05.389";
"14.849758469419463" "15:09:05.58";
"18.871073102595346" "15:09:05.581";
"17.705181430195292" "15:09:05.583";
"15.210479632070335" "15:09:05.584";
"15.271612129811938" "15:09:05.586";
"15.239694384551788" "15:09:05.588";
"14.305067502483697" "15:09:05.589";
"13.482562423502614" "15:09:05.591";
"12.876001722397385" "15:09:05.593";
"12.805594833628685" "15:09:05.594";
"12.457738102830431" "15:09:05.596";
"11.605904207061823" "15:09:05.598";
"11.194548534361212" "15:09:05.599";
"11.240457357737089" "15:09:05.601";
"11.190018463507247" "15:09:05.603";
"10.54234617836919" "15:09:05.604";
"9.771247999681453" "15:09:05.606";
"9.971748447964599" "15:09:05.608";
"9.961127631863747" "15:09:05.609";
"9.585688886522933" "15:09:05.611";
"8.831398877508802" "15:09:05.612";
"8.772680290843674" "15:09:05.614";
"8.892681182576261" "15:09:05.615";
"8.47483656117191" "15:09:05.617";
"8.051257038356551" "15:09:05.619";
"7.836416224615945" "15:09:05.621";
"7.6973318939204525" "15:09:05.623";
"7.56861396088807" "15:09:05.638";
"7.21708085332945" "15:09:05.639";
"7.041396300087717" "15:09:05.641";
"6.976721507394556" "15:09:05.643";
"6.774189584976573" "15:09:05.645";
"6.606433801017014" "15:09:05.646";
"6.491238479787147" "15:09:05.648";
"6.465407956324111" "15:09:05.65";
"6.33759611310141" "15:09:05.652";
"6.2663196306267235" "15:09:05.654";
"6.243419793892277" "15:09:05.657";
"6.32328832395769" "15:09:05.658";
"6.619533222380913" "15:09:05.66";
"6.683464169573939" "15:09:05.662";
"6.226088942942654" "15:09:05.663";
"5.837544256661132" "15:09:05.665";
"6.001552472892467" "15:09:05.667";
"6.142284234786329" "15:09:05.669";
"5.713154173254792" "15:09:05.67";
"5.627224369412082" "15:09:05.672";
"5.798986357032883" "15:09:05.674";
"5.5638203596149545" "15:09:05.675";
"77.91969" "15:11:53.607";
"88.45319" "15:11:54.244";
"94.22918" "15:11:54.709";
"96.63985" "15:11:55.078";
"97.462326" "15:11:55.439";
"97.221504" "15:11:55.793";
"96.14967" "15:11:56.172";
"94.2463" "15:11:56.536";
"91.18372" "15:11:56.845";
"85.45838" "15:11:57.249";
"84.35536" "15:11:57.636";
"86.03393" "15:11:58.014";
"89.09288" "15:11:58.443";
"90.427055" "15:11:58.786";
"90.72829" "15:11:59.117";
"90.3256" "15:11:59.496";
"89.34307" "15:11:59.824";
"87.82792" "15:12:00.2";
"85.00899" "15:12:00.591";
"80.95103" "15:12:01.045";
"78.994484" "15:12:01.431";
"83.956985" "15:12:01.817";
"81.82451" "15:12:02.181";
"78.20917" "15:12:02.474";
"79.43976" "15:12:02.809";
"80.80074" "15:12:03.23";
"80.69117" "15:12:03.623";
"78.97905" "15:12:04.01";
"82.39736" "15:12:04.375";
"77.42406" "15:12:04.68";
"77.60327" "15:12:05.089";
"77.83742" "15:12:05.447";
"78.21078" "15:12:05.801";
"78.11332" "15:12:06.13";
"76.9447" "15:12:06.588";
"76.84605" "15:12:06.885";
"75.77428" "15:12:07.186";
"76.31214" "15:12:07.566";
"76.18572" "15:12:07.958";
"76.10903" "15:12:08.369";
"75.55867" "15:12:08.703";
"75.36168" "15:12:09.108";
"74.78529" "15:12:09.51";
"74.687485" "15:12:09.894";
"74.37139" "15:12:10.237";
"74.36885" "15:12:10.681";
"74.062386" "15:12:11.04";
"74.10397" "15:12:11.398";
"73.71966" "15:12:11.787";
"73.50538" "15:12:12.149";
"103.78020466853134" "15:12:12.151";
"71.39917545838638" "15:12:12.152";
"43.835507904825384" "15:12:12.154";
"26.037818480506466" "15:12:12.155";
"22.466109727815276" "15:12:12.157";
"25.573757293590923" "15:12:12.158";
"29.028966350800705" "15:12:12.16";
"30.835888518567454" "15:12:12.162";
"30.104994695305347" "15:12:12.164";
"27.536860817441028" "15:12:12.165";
"24.383022026604706" "15:12:12.167";
"21.29960591113438" "15:12:12.169";
"17.747236100309646" "15:12:12.17";
"14.206853881666227" "15:12:12.172";
"11.992404667204356" "15:12:12.173";
"11.588198002152236" "15:12:12.175";
"11.7388534104683" "15:12:12.176";
"11.980444916622721" "15:12:12.178";
"12.587588881517497" "15:12:12.18";
"12.769639450593779" "15:12:12.181";
"11.583519151470709" "15:12:12.183";
"9.501164335292207" "15:12:12.185";
"8.654295580582641" "15:12:12.187";
"9.01548765256084" "15:12:12.189";
"8.594953734742683" "15:12:12.191";
"8.147538233592984" "15:12:12.192";
"8.471366362882927" "15:12:12.194";
"8.24035846217022" "15:12:12.196";
"7.523414724306299" "15:12:12.198";
"7.189982395776787" "15:12:12.2";
"6.508849669389" "15:12:12.202";
"5.851533928189284" "15:12:12.203";
"6.4037789888217125" "15:12:12.205";
"6.783906378051663" "15:12:12.207";
"6.618028982643132" "15:12:12.208";
"6.232803966905912" "15:12:12.21";
"5.8461338103875065" "15:12:12.212";
"5.4375080999554735" "15:12:12.214";
"5.300198820924965" "15:12:12.216";
"5.587699213538676" "15:12:12.217";
"5.40344208807284" "15:12:12.219";
"5.291385393476899" "15:12:12.227";
"4.984261057651273" "15:12:12.229";
"4.881643671665893" "15:12:12.231";
"4.93980303827467" "15:12:12.232";
"4.802080297751322" "15:12:12.234";
"4.744227552063704" "15:12:12.236";
"4.524137289118702" "15:12:12.238";
"4.389026696913442" "15:12:12.239";
"4.450599357746215" "15:12:12.241";
"64.41709" "15:12:13.458";
"63.632217" "15:12:14.301";
"46.92427" "15:12:15.097";
"32.96313" "15:12:15.942";
"29.291273" "15:12:16.558";
"25.644917" "15:12:17.104";
"23.561989" "15:12:17.76";
"52.707478" "15:12:18.683";
"22.207323" "15:12:19.651";
"23.706646" "15:12:20.61";
"14.20956" "15:12:21.684";
"25.883623" "15:12:22.636";
"18.306864" "15:12:23.996";
"6.875953" "15:12:25.484";
"11.9977" "15:12:26.321";
"7.160552" "15:12:27.511";
"14.344087" "15:12:28.546";
"8.279809" "15:12:29.848";
"7.8758035" "15:12:30.775";
"6.038576" "15:12:31.813";
"2.2966294" "15:12:32.899";
"4.332708" "15:12:33.531";
"6.4547825" "15:12:34.351";
"5.7979975" "15:12:35.247";
"4.0804954" "15:12:36.187";
"2.7464151" "15:12:36.936";
"2.3627076" "15:12:37.706";
"3.2782743" "15:12:38.453";
"2.6265779" "15:12:39.358";
"3.2288523" "15:12:40.17";
"2.7752945" "15:12:40.948";
"1.186291" "15:12:41.789";
"1.6154346" "15:12:42.472";
"2.494385" "15:12:43.135";
"1.9619559" "15:12:43.949";
"1.1906984" "15:12:44.74";
"1.4041933" "15:12:45.398";
"1.347243" "15:12:46.142";
"0.81995565" "15:12:46.901";
"1.5716455" "15:12:47.564";
"0.86126065" "15:12:48.435";
"0.6891742" "15:12:49.155";
"1.0051694" "15:12:49.865";
"0.7354157" "15:12:50.641";
"0.67939365" "15:12:51.28";
"0.8814287" "15:12:51.875";
"0.53287417" "15:12:52.639";
"0.4200542" "15:12:53.325";
"0.66252524" "15:12:53.93";
"0.49660772" "15:12:54.602";
"17.675070570669615" "15:12:54.604";
"14.570600872069344" "15:12:54.606";
"12.050048514403779" "15:12:54.607";
"9.678179637714395" "15:12:54.608";
"10.152522356746942" "15:12:54.61";
"10.22044639649735" "15:12:54.612";
"8.04914110373006" "15:12:54.613";
"7.032138587621534" "15:12:54.615";
"8.230037639094844" "15:12:54.616";
"8.112301095722152" "15:12:54.618";
"6.82777864639771" "15:12:54.62";
"6.175164073087483" "15:12:54.621";
"5.709938230935337" "15:12:54.623";
"5.843172529033399" "15:12:54.624";
"5.980040523007635" "15:12:54.626";
"5.1580553635125685" "15:12:54.627";
"5.092255952432414" "15:12:54.629";
"5.459202269494638" "15:12:54.631";
"5.049653865780248" "15:12:54.632";
"4.467758722562376" "15:12:54.634";
"4.324186417796585" "15:12:54.636";
"4.752285597590552" "15:12:54.637";
"4.531041848894491" "15:12:54.639";
"4.09737535292966" "15:12:54.641";
"4.204933133294741" "15:12:54.647";
"3.9613821548000945" "15:12:54.649";
"3.9099247196684535" "15:12:54.651";
"3.6989380596096066" "15:12:54.652";
"3.6227483875970634" "15:12:54.654";
"3.5715589688151126" "15:12:54.656";
"3.532335251351003" "15:12:54.658";
"3.447700079042611" "15:12:54.659";
"3.296995317481246" "15:12:54.661";
"3.2295508862437052" "15:12:54.663";
"3.259050918846149" "15:12:54.665";
"3.006368540790016" "15:12:54.666";
"3.1330418058276033" "15:12:54.668";
"2.97521086419771" "15:12:54.671";
"2.8645573897595584" "15:12:54.673";
"2.8758895799615374" "15:12:54.675";
"2.827000684053438" "15:12:54.677";
"2.894350970873364" "15:12:54.678";
"3.0104642745408046" "15:12:54.68";
"2.9511782418400787" "15:12:54.682";
"2.6087699579605776" "15:12:54.684";
"2.4392316464460313" "15:12:54.685";
"2.4049713120039122" "15:12:54.687";
"2.503992372666395" "15:12:54.689";
"2.632122117518553" "15:12:54.691";
"2.6686821823890137" "15:12:54.692";
"124.01136" "15:12:55.768";
"120.530655" "15:12:56.904";
"130.76521" "15:12:57.876";
"128.29715" "15:12:59.114";
"132.95764" "15:13:02.378";
"133.61787" "15:13:02.963";
"131.75017" "15:13:03.371";
"128.18472" "15:13:03.796";
"123.39702" "15:13:04.245";
"118.72645" "15:13:04.737";
"115.100555" "15:13:05.239";
"111.79445" "15:13:05.738";
"109.15439" "15:13:06.218";
"112.42572" "15:13:07.253";
"111.7254" "15:13:08.565";
"107.00366" "15:13:09.416";
"108.18689" "15:13:09.945";
"108.58419" "15:13:10.436";
"108.5795" "15:13:11.015";
"108.48891" "15:13:11.459";
"108.35627" "15:13:11.848";
"108.194405" "15:13:12.256";
"108.10015" "15:13:12.729";
"107.76345" "15:13:13.326";
"107.23125" "15:13:13.933";
"106.53842" "15:13:14.611";
"105.88434" "15:13:15.007";
"105.334015" "15:13:15.438";
"104.64102" "15:13:16.08";
"103.829285" "15:13:16.664";
"102.915405" "15:13:17.288";
"101.906334" "15:13:17.909";
"100.79671" "15:13:18.449";
"99.87119" "15:13:18.911";
"98.782524" "15:13:19.475";
"97.74554" "15:13:19.908";
"96.41364" "15:13:20.404";
"94.55979" "15:13:20.859";
"91.341" "15:13:21.272";
"95.19458" "15:13:21.788";
"96.5668" "15:13:22.447";
"95.09491" "15:13:22.994";
"94.69362" "15:13:23.822";
"93.513855" "15:13:24.176";
"87.31163" "15:13:24.688";
"91.13651" "15:13:25.186";
"92.91186" "15:13:25.539";
"93.84257" "15:13:25.932";
"94.31615" "15:13:26.388";
"94.47222" "15:13:26.814";
"158.3530547348604" "15:13:26.816";
"133.21145896256962" "15:13:26.817";
"109.95818779174513" "15:13:26.818";
"89.37608376298566" "15:13:26.82";
"71.49050970650565" "15:13:26.821";
"57.846327114885334" "15:13:26.822";
"47.12060542162359" "15:13:26.823";
"37.89888007616456" "15:13:26.825";
"31.234619085237597" "15:13:26.826";
"27.780180162198516" "15:13:26.827";
"26.97730221566408" "15:13:26.828";
"28.624742707782985" "15:13:26.833";
"31.09468259053384" "15:13:26.834";
"32.250395519855594" "15:13:26.836";
"30.987321283873293" "15:13:26.837";
"27.367452840014725" "15:13:26.839";
"22.172456814405077" "15:13:26.84";
"16.92889659112651" "15:13:26.842";
"14.416349034228634" "15:13:26.843";
"14.881528572265081" "15:13:26.845";
"15.034067016651338" "15:13:26.847";
"14.532645649427373" "15:13:26.848";
"14.12980102068133" "15:13:26.85";
"13.331876850835851" "15:13:26.852";
"11.734434891357704" "15:13:26.853";
"9.79347304696053" "15:13:26.855";
"8.923428367225746" "15:13:26.857";
"9.500625495678467" "15:13:26.858";
"9.889106080223089" "15:13:26.86";
"9.576257519493089" "15:13:26.862";
"9.022203282268617" "15:13:26.864";
"8.230390385103814" "15:13:26.865";
"7.010404205800276" "15:13:26.867";
"6.027325078165178" "15:13:26.868";
"6.367077396523484" "15:13:26.87";
"6.500906299478687" "15:13:26.872";
"6.185294860217844" "15:13:26.873";
"5.905033135434491" "15:13:26.875";
"5.147103971289184" "15:13:26.876";
"4.755633322769082" "15:13:26.878";
"5.22930110660714" "15:13:26.879";
"5.097657354332236" "15:13:26.881";
"4.534309628634492" "15:13:26.883";
"3.993775575867627" "15:13:26.884";
"4.193965508165071" "15:13:26.886";
"4.521357031684894" "15:13:26.888";
"4.183303506044945" "15:13:26.89";
"3.7814488560591353" "15:13:26.891";
"3.3334368323626533" "15:13:26.893";
"3.623635121442286" "15:13:26.895";
"110.57353" "15:13:51.7";
"86.8244" "15:13:53.062";
"75.86707" "15:13:54.028";
"80.52017" "15:13:55.181";
"85.348114" "15:13:55.98";
"89.67025" "15:13:56.805";
"87.89375" "15:13:57.68";
"82.94136" "15:13:58.481";
"80.68053" "15:13:59.305";
"79.259125" "15:14:00.198";
"77.48437" "15:14:01.046";
"75.9717" "15:14:09.98";
"75.012375" "15:14:16.305";
"74.05318" "15:14:20.039";
"73.15254" "15:14:21.313";
"73.09031" "15:14:22.013";
"73.19003" "15:14:22.772";
"70.64611" "15:14:23.539";
"72.2913" "15:14:24.406";
"69.248245" "15:14:25.049";
"70.624664" "15:14:25.979";
"74.0138" "15:14:27.701";
"75.61718" "15:14:28.513";
"76.62279" "15:14:29.01";
"77.35976" "15:14:29.795";
"77.80881" "15:14:30.175";
"78.013214" "15:14:30.579";
"77.79789" "15:14:31.101";
"77.46091" "15:14:32.717";
"76.46782" "15:14:33.256";
"74.392006" "15:14:34.05";
"71.76419" "15:14:34.821";
"70.16536" "15:14:35.617";
"74.769485" "15:14:36.476";
"69.258385" "15:14:37.203";
"68.651634" "15:14:38.079";
"69.92104" "15:14:38.893";
"69.79346" "15:14:39.805";
"70.064224" "15:14:40.852";
"71.08431" "15:14:41.71";
"71.802895" "15:14:42.492";
"73.09411" "15:14:52.021";
"73.53909" "15:14:53.08";
"73.58596" "15:14:54.119";
"73.637695" "15:14:55.077";
"73.821686" "15:14:55.634";
"73.91312" "15:14:56.178";
"73.79058" "15:15:00.104";
"73.53084" "15:15:01.445";
"73.30578" "15:15:01.951";
"96.27300961370284" "15:15:01.953";
"78.43964693255776" "15:15:01.955";
"62.89331744201945" "15:15:01.957";
"47.19305295496835" "15:15:01.958";
"33.2131786101106" "15:15:01.96";
"23.211070527480693" "15:15:01.961";
"17.833322001818093" "15:15:01.963";
"18.78514599416313" "15:15:01.965";
"22.42366990923824" "15:15:01.967";
"24.148420827834173" "15:15:01.968";
"23.93792049010538" "15:15:01.97";
"22.849654311851463" "15:15:01.972";
"21.26536706481854" "15:15:01.974";
"18.80688075798271" "15:15:01.975";
"15.685917335595885" "15:15:01.977";
"13.196725287320543" "15:15:01.978";
"11.256649177207132" "15:15:01.98";
"9.661792943746024" "15:15:01.982";
"10.26180257958145" "15:15:01.983";
"10.852587373452849" "15:15:01.985";
"11.396074140007112" "15:15:01.987";
"11.894937219323069" "15:15:01.988";
"11.35520112096304" "15:15:01.99";
"10.427744469785017" "15:15:01.992";
"9.526417877752873" "15:15:01.994";
"8.28665452707891" "15:15:01.996";
"7.849091548214682" "15:15:01.998";
"7.669730159660486" "15:15:02";
"7.861044694877725" "15:15:02.001";
"8.04182667626384" "15:15:02.003";
"7.871354809346516" "15:15:02.005";
"7.627199538480839" "15:15:02.011";
"6.819161819138998" "15:15:02.012";
"6.3060252420365" "15:15:02.014";
"5.771065768577108" "15:15:02.016";
"5.784233780445008" "15:15:02.017";
"5.845757087879443" "15:15:02.019";
"5.927953817799761" "15:15:02.021";
"5.692619285524822" "15:15:02.023";
"5.263192174206345" "15:15:02.025";
"4.914356844013468" "15:15:02.027";
"4.6623280580615" "15:15:02.029";
"4.7387990784480545" "15:15:02.031";
"4.780357490230598" "15:15:02.033";
"4.669946580162062" "15:15:02.035";
"4.569636303179074" "15:15:02.037";
"4.456541512195911" "15:15:02.039";
"4.302276110299467" "15:15:02.04";
"4.165433772364541" "15:15:02.042";
"4.036534754423933" "15:15:02.044";
"66.81997" "15:15:14.322";
"59.855667" "15:15:15.771";
"37.67311" "15:15:17.079";
"59.809708" "15:15:18.092";
"83.349525" "15:15:21.05";
"859.32135" "15:15:24.092";
"1505.0326" "15:15:27.17";
"1782.3123" "15:15:49.881";
"1940.5017" "15:15:50.962";
"2057.049" "15:15:51.849";
"2008.4354" "15:15:54.444";
"1695.337" "15:16:12.839";
"1547.6293" "15:16:14.358";
"1142.6104" "15:16:16.373";
"741.1111" "15:16:19.662";
"821.3444" "15:16:33.699";
"187.79512" "15:16:51.611";
"172.14047" "15:16:54.041";
"164.72462" "15:16:55.532";
"159.62753" "15:16:57.177";
"154.88857" "15:16:59.282";
"150.65578" "15:17:05.735";
"146.64275" "15:17:12.948";
"143.15862" "15:17:13.432";
"140.07317" "15:17:14.029";
"137.26106" "15:17:17.153";
"134.77591" "15:17:17.761";
"132.59354" "15:17:21.48";
"130.64209" "15:17:22.095";
"128.69514" "15:17:33.931";
"126.85361" "15:17:34.528";
"125.0247" "15:17:35.128";
"123.163956" "15:17:35.754";
"121.26326" "15:17:36.403";
"119.48151" "15:17:40.147";
"117.63254" "15:17:40.809";
"115.721245" "15:17:41.692";
"113.977936" "15:17:42.349";
"112.719696" "15:17:48.318";
"113.039085" "15:17:49.479";
"114.776535" "15:17:52.55";
"117.75608" "15:17:54.216";
"120.59915" "15:17:55.607";
"119.04658" "15:17:58.143";
"118.475395" "15:18:00.489";
"115.34969" "15:18:03.021";
"97.04737" "15:18:05.972";
"94.84225" "15:18:06.501";
"97.21866" "15:18:06.945";
"99.53592" "15:18:07.431";
]

function parse_dates(x)
    t=split(x,(':','.'))
    w=map(x->parse(Float64,x),t)
    Dates.Time(w...)
end
function norm_time(x)
    abs_times = x
    start_time = abs_times[1]
    norm_times = Vector(undef,length(abs_times))
    for i in 1:length(abs_times)
        norm_times[i] = abs_times[i]-start_time
    end
    norm_times
end

loss_values = map(x->parse(Float64,x),a[:,1])
time_values = map(x->parse_dates(x),a[:,2])
norm_time_vals = norm_time(time_values)

plot(time_values[1:50],loss_values[1:50],linecolor="blue")
plot!(time_values[51:100],loss_values[51:100],linecolor="red")
plot!(time_values[101:150],loss_values[101:150],linecolor="blue")
plot!(time_values[151:150],loss_values[151:150],linecolor="red")
plot!(time_values[201:250],loss_values[201:250],linecolor="blue")
plot!(time_values[251:300],loss_values[251:300],linecolor="red")
plot!(time_values[301:350],loss_values[301:350],linecolor="blue")
plot!(time_values[351:400],loss_values[351:400],linecolor="red")
plot!(time_values[401:450],loss_values[401:450],linecolor="blue")
plot!(time_values[451:500],loss_values[451:500],linecolor="red")
Int64(norm_time_vals[1])

print(loss_values[500])


text:
function welcome_message()

    println("\n \n")
    println("___________________________________________________________")
    println(" _  __ _____    _____   ____   _  __      __ ______  _____ ")
    println("| |/ // ____|  / ____| / __ \\ | | \\ \\    / /|  ____||  __ \\ ")
    println("| ' /| |  __  | (___  | |  | || |  \\ \\  / / | |__   | |__) |")
    println("|  < | | |_ |  \\___ \\ | |  | || |   \\ \\/ /  |  __|  |  _  / ")
    println("| . \\| |__| |  ____) || |__| || |____\\  /   | |____ | | \\ \\ ")
    println("|_|\\_\\\\_____| |_____/  \\____/ |______|\\/    |______||_|  \\_\\")
    println("___________________________________________________________")
    println(" ")
    println(" Welcome to the Klein gordon field Solver")
    println("___________________________________________________________")

end


text:
# def normal_boundary_weights(nobjs, divisions_outer, divisions_inner=0):
#     def generate_recursive(weights, weight, left, total, index):
#         if index == nobjs - 1:
#             weight[index] = float(left) / float(total)
#             weights.append(copy.copy(weight))
#         else:
#             for i in range(left+1):
#                 weight[index] = float(i) / float(total)
#                 generate_recursive(weights, weight, left-i, total, index+1)
#
#     def generate_weights(divisions):
#         weights = []
#         generate_recursive(weights, [0.0]*nobjs, divisions, divisions, 0)
#         return weights
#
#     weights = generate_weights(divisions_outer)
#
#     if divisions_inner > 0:
#         inner_weights = generate_weights(divisions_inner)
#
#         for i in range(len(inner_weights)):
#             weight = inner_weights[i]
#
#             for j in range(len(weight)):
#                 weight[j] = (1.0 / nobjs + weight[j]) / 2.0
#
#             weights.append(weight)
#
#     return weights

"""
    generate_regular_reference_points(M, p)

Generates reference points following the Das and Dennis systematic approach.
Notation follows Deb and Jain (NSGA-III paper)
M = number of objectives
p = number of divisions
"""
function generate_regular_reference_points(M, p)
    weights = zeros(0,M)
    weights = generate_recursive(weights, zeros(1, M), p, p, M, 1)
    return weights
end

function generate_recursive(weights, weight, left, total, M, index)
    if index == M
        # termination case
        weight[index] = left / total
        weights = [weights; weight]
    else
        for i in 0:left
            weight[index] = i / total
            weights = generate_recursive(weights, weight, left-i, total, M, index+1)
        end
    end
    return weights
end

max_stars_repo_name: phawthorne/echidna

max_issues_repo_name: phawthorne/echidna

max_forks_repo_name: phawthorne/echidna


text:
#Global Setting
use_gpu = false
use_darrays = false #Only relevant if use_gpu = false

random_seed = 1234567
random_seed = nothing

#Data Loading specifics
data_path = "data/multi-d-gaussians//"
data_prefix = "g"
groups_count = 9
global_preprocessing = nothing
local_preprocessing = nothing

#Model Parameters
iterations = 100
hard_clustering = false

total_dim = 3
local_dim = 3

α = 500.0
γ = 10000.0
global_weight = 1.0
local_weight= 1.0

initial_global_clusters = 1
initial_local_clusters = 1 #this is per group

use_dict_for_global = false

split_delays = false

ignore_local = false

split_stop = 10
argmax_sample_stop = 10

use_mean_for_merge = false

global_hyper_params = niw_hyperparams(1.0,
    zeros(2),
    20.0,
    Matrix{Float64}(I, 2, 2)*1)

local_hyper_params = niw_hyperparams(1.0,
    zeros(1),
    20.0,
    Matrix{Float64}(I, 1, 1)*1)

# #Model Parameters
# iterations = 10
#
# total_dim = 3
# local_dim = 3
#
# α = 4.0
# γ = 2.0
# global_weight = 1.0
# local_weight= 1.0
#
#
# global_hyper_params = niw_hyperparams(5.0,
#     ones(3)*2.5,
#     10.0,
#     [[0.8662817 0.78323282 0.41225376];[0.78323282 0.74170384 0.50340258];[0.41225376 0.50340258 0.79185577]])
#
# local_hyper_params = niw_hyperparams(1.0,
#     [217.0,510.0],
#     10.0,
#     Matrix{Float64}(I, 2, 2)*0.5)


text:
export viz_policy_drl
function viz_policy_drl(neuralNetworkPath::AbstractString, alpha::Matrix{Float64}=[1.0 1.0]; batch_size::Int64=500,XandY::Bool = false)

    nnet = NNet(neuralNetworkPath);
    if length(alpha)!=2
        if !XandY
            grid = RectangleGrid(Ranges, Thetas, Bearings, Speeds, Speeds)
        else
            grid = RectangleGrid(Xs, Ys, Bearings, Speeds, Speeds)
        end
        policy = read_policy(Actions', alpha)
    end

    c = RGB{U8}(1.,1.,1.) # white
    e = RGB{U8}(.0,.0,.5) # pink
    a = RGB{U8}(.0,.600,.0) # green
    d = RGB{U8}(.5,.5,.5) # grey
    b = RGB{U8}(.7,.9,.0) # neon green
    f = RGB{U8}(0.94,1.0,.7) # pale yellow
    colors =[a; b; f; c; d; e]

    @manipulate for psi_int  = round(rad2deg(psis_drl)),
        v_own = sos_drl,
        v_int = sis_drl,
        zoom = [4, 3, 2, 1.5,1],
        nbin = [100,150,200,250],
        interp = [true, false]

        #mat =  ccall((:load_network,LIB_BLAS),Ptr{Void},(Ptr{UInt8},),neuralNetworkPath)
        inputsNet= zeros(nbin*nbin,STATE_DIM)
        ind = 1
        for i=linspace(round(Int,-1*RangeMax/zoom),round(Int,RangeMax/zoom),nbin)
            for j=linspace(round(Int,-1*RangeMax/zoom),round(Int,RangeMax/zoom),nbin)
                r = sqrt(i^2+j^2)
                th = atan2(j,i)
                inputsNet[ind,:] = [r,th,deg2rad(psi_int),v_own,v_int];
                ind = ind+1
            end
        end

        q_nnet = zeros(nbin*nbin,ACTION_DIM);
        ind = 1
        while ind+batch_size<nbin*nbin
            input = inputsNet[ind:(ind+batch_size-1),:]'
            output = evaluate_network_multiple(nnet,input)
            q_nnet = [q_nnet[1:(ind-1),:];output';q_nnet[ind+batch_size:end,:]]
            ind=ind+batch_size
        end
        input = inputsNet[ind:end,:]'
        output = evaluate_network_multiple(nnet,input)
        q_nnet = [q_nnet[1:(ind-1),:];output']

        ind = 1
        function get_heat1(x::Float64, y::Float64)
            if !XandY
                r = norm([x,y])
                th = atan2(y,x)
                qvals = evaluate(policy, get_belief(
                [r, th, deg2rad(psi_int), v_own, v_int], grid,interp,true))
            else
                qvals = evaluate(policy, get_belief(
                [x, y, deg2rad(psi_int), v_own, v_int], grid,interp,true,true))
            end
            return Actions[indmax(qvals)]
        end # function get_heat1

        function get_heat2(x::Float64, y::Float64)
            r = sqrt(x^2+y^2)
            th = atan2(y,x)
            action  = Actions[indmax(q_nnet[ind,:])]
            ind = ind+1
            return action
        end # function get_heat2

        g = GroupPlot(2, 1, groupStyle = "horizontal sep=3cm")
        if length(alpha)!=2
            push!(g, Axis([
                Plots.Image(get_heat1, (round(Int,-1*RangeMax/zoom), round(Int,RangeMax/zoom)),
                (round(Int,-1*RangeMax/zoom), round(Int,RangeMax/zoom)),
                           zmin = -20, zmax = 20,
                           xbins = nbin, ybins = nbin,
                           colormap = ColorMaps.RGBArray(colors), colorbar = false),
               Plots.Node(L">", 0, 0, style=string("font=\\Huge")),
               Plots.Node(L">", 2500/zoom, 2500/zoom, style=string("rotate=", psi_int, ",font=\\Huge"))
                ], width="10cm", height="10cm", xlabel="x (ft)", ylabel="y (ft)", title="Discrete Value Iteration Policy"))
        end
        ind = 1
        push!(g, Axis([
            Plots.Image(get_heat2, (round(Int,-1*RangeMax/zoom), round(Int,RangeMax/zoom)),
            (round(Int,-1*RangeMax/zoom), round(Int,RangeMax/zoom)),
                       zmin = -20, zmax = 20,
                       xbins = nbin, ybins = nbin,
                       colormap = ColorMaps.RGBArray(colors)),
           Plots.Node(L">", 0, 0, style=string("font=\\Huge")),
           Plots.Node(L">", 2500/zoom, 2500/zoom, style=string("rotate=", psi_int, ",font=\\Huge"))
            ], width="10cm", height="10cm", xlabel="x (ft)", ylabel="y (ft)", title="Neural Net Policy"))
        g
    end # for p_int, v0, v1, pa, ta
end # function viz_policy_drl


text:
#=

# 9: Kepler vortex (explicit, symplectic, and reversible)

```@raw html
	<img src='../assets/Kepler.png' alt='missing' width="50%" height="50%" /><br>
```
Kepler vortex, where a fluid forms a Gaussian ring rotating around the center of gravity. The ring is centered around ``r0=10`` with standard deviation ``2.5`` (approximately nine thousand particles). Strength of the gravity is ``GM=1000``, and speed of sound is chosen as ``c=0.001 v_\varphi``, where ``v_\varphi`` is the velocity at ``r_0``. The final time is approximately 63, which corresponds to ten revolutions of the ring.
=#

module Kepler_vortex

using Printf
using SmoothedParticles
using Parameters
using Plots
using DataFrames # to store the csv file
using CSV # to store the csv file
using QuadGK #To calculate the integral over Σ
using Interpolations
using Roots
include("utils/FixPA.jl")
using .FixPA

#=
Declare constant parameters
=#
##graviational
const r0 = 10.0 #central ring radius
const GM = 1000. #graviational foce

#=
Calculate the Gaussian distribution of particles in the vortex, according to Cartwright et al.
=#
function vphi_r(r::Float64)::Float64
	return sqrt(GM)/sqrt(r)
end
const vphi0 = vphi_r(r0)
@show vphi0
const omega0= vphi0/r0
@show omega0

function Σ(r::Float64)::Float64
	return 2*pi*r*exp(-30*(1-r/r0)^2)
end
const denominator = quadgk(Σ, 0, 40, rtol=1.0e-06)[1]
@show denominator
function f(r::Float64)::Float64
	return quadgk(Σ, 0, r, rtol=1.0e-03)[1]/denominator
end

const rs = 0.:0.5:25.0 # collection of radii
const f_table = [f(r) for r in rs]
const f_interpolation = interpolate(f_table, BSpline(Cubic(Line(OnGrid()))))
const f_interpolation_scaled = scale(f_interpolation, rs)
function r_f(F::Float64)::Float64 #radius as a function of the uniformly distributed random variable
	return find_zero(x -> (f_interpolation_scaled(x)-F), r0)
end
@show r_f(0.5)

## computational
#const rs_in_vortex = 7.:dr:14. # equidistant rings in the vortex. Boundaries chosen as the edges of the histogram for 10^4 points over the interval [0,1]
const N_rings = 25
@show N_rings
const rs_in_vortex = [r_f(u) for u in 0.01:(0.99-0.01)/N_rings:0.99] # Gaussian rings
dr = r_f(0.25+1/N_rings)-r_f(0.25)
@show dr

const h = 3.0*dr           #size of kernel support
const rho0 = 1.   	   #fluid density
const m = rho0*dr^2        #particle mass
#const mu = 0.0#8.4e-4          #dynamic viscosity of water

##geometrical
const box_width = 4*r0
const wall_width = 2.5*dr

##artificial
const c = 0.01             #numerical speed of sound
const dr_wall = 0.95*dr
const E_wall = GM/r0
const eps = 1e-16

##temporal
const dt = 0.0001*h/c
@show dt
const t_end = 10 * 2 * pi / omega0 #ten revolutions
@show t_end
const dt_frame = t_end/200

##particle types
const FLUID = 0.
const WALL = 1.

@with_kw mutable struct Particle <: AbstractParticle
	x::RealVector #position
	v::RealVector = VEC0 #velocity
	a::RealVector = VEC0 #acceleration
	P::Float64 = 0. #pressure
	rho::Float64 = 0. #density
    rho0::Float64 = 0.
	type::Float64 #particle_type
end

#=
Define geometry and make particles
=#

function generate_circle!(sys::ParticleSystem, r::Float64, dφ::Float64; φ_start=0.0, x_center=0.0, y_center=0.0, vφ= 0.0)
	φ = φ_start
	while φ < 2*pi+φ_start
		x = cos(φ) #x position on unit circle
		y = sin(φ) #y position on unit circle
		particle_position = RealVector(x_center+r*x, y_center+r*y, 0.)
		particle_velocity = RealVector(-vφ*y, vφ*x, 0.)
		push!(sys.particles, Particle(x=particle_position, v=particle_velocity, type = FLUID))
		φ += dφ
	end
end

function make_system()
    domain = Rectangle(-box_width, -box_width, box_width, box_width)
    sys = ParticleSystem(Particle, domain, h)

	dφ = rs_in_vortex[2]/rs_in_vortex[1]-1.0 # increment in the angle
	for i in 1:length(rs_in_vortex)-1 # not the last one because dphi would be unknown
		r = rs_in_vortex[i]
		generate_circle!(sys, r, dφ; vφ = vphi_r(r))
		dφ = (rs_in_vortex[i+1]-r)/r
	end

	return sys
end

#=
Define particle interactions
=#

@inbounds function find_rho!(p::Particle, q::Particle, r::Float64)
    if p.type == FLUID && q.type == FLUID
		p.rho += m*wendland2(h,r)
	end
end

@inbounds function find_rho0!(p::Particle, q::Particle, r::Float64)
    if p.type == FLUID && q.type == FLUID
		p.rho0 += m*wendland2(h,r)
	end
end

function find_pressure!(p::Particle)
	p.P = c^2*(p.rho - p.rho0)
end

@inbounds function internal_force!(p::Particle, q::Particle, r::Float64)
	if p.type == FLUID && q.type == FLUID
		ker = m*rDwendland2(h,r)
		p.a += -ker*(p.P/rho0^2 + q.P/rho0^2)*(p.x - q.x)
		#p.a += +2*ker*mu/rho0^2*(p.v - q.v)
	elseif p.type == FLUID && q.type == WALL && r < dr_wall
		s = dr_wall/(r + eps)
		p.a += -E_wall/(r + eps)^2*(s^2 - s^4)*(p.x - q.x)
	end
end

function reset_a!(p::Particle)
    p.a = zero(RealVector)
end

function reset_rho!(p::Particle)
    p.rho = 0.0
end

function move!(p::Particle)
	if p.type == FLUID
		p.x = rev_add(p.x, dt*p.v)
	end
end

function accelerate!(p::Particle)
	if p.type == FLUID
		p.v = rev_add(p.v, 0.5*dt*(rev_add(p.a, -GM/(norm(p.x)^3)*p.x)))
	end
end

function LJ_potential(p::Particle, q::Particle, r::Float64)::Float64
	if q.type == WALL && p.type == FLUID && r < dr_wall
		s = dr_wall/(r + eps)
		return m*E_wall*(0.5s^2 - 0.25s^4 -0.25)
	else
		return 0.0
	end
end

function energy_kinetic(sys::ParticleSystem)::Float64
	return sum(p -> 0.5*m*dot(p.v, p.v), sys.particles)
end

function energy_potential(sys::ParticleSystem)::Float64
	return sum(p -> -GM*m/(norm(p.x)), sys.particles)
end

function energy_internal(sys::ParticleSystem)::Float64
	return sum(p -> 0.5*m*c^2*(p.rho - p.rho0)^2/rho0^2, sys.particles)
end

function energy(sys::ParticleSystem, p::Particle)::Float64
	kinetic = 0.5*m*dot(p.v, p.v)
	internal =  0.5*m*c^2*(p.rho - p.rho0)^2/rho0^2
	gravity_potential = - GM*m/(norm(p.x))
	wall_potential = SmoothedParticles.sum(sys, LJ_potential, p)
	return kinetic + internal + gravity_potential + wall_potential
end

#=
Put everything into a time loop
=#

function verlet_step!(sys::ParticleSystem)
    apply!(sys, accelerate!)
    apply!(sys, move!)
    create_cell_list!(sys)
    apply!(sys, reset_rho!)
    apply!(sys, find_rho!, self = true)
    apply!(sys, find_pressure!)
    apply!(sys, reset_a!)
    apply!(sys, internal_force!)
    apply!(sys, accelerate!)
end

function save_results!(out::SmoothedParticles.DataStorage, sys::ParticleSystem, k::Int64)
    if (k %  Int64(round(dt_frame/dt)) == 0)
        @printf("t = %.6e\n", k*dt)
        #energy
        Etot = sum(p -> energy(sys,p), sys.particles)
        @show Etot
        Ekin = energy_kinetic(sys)
        @show Ekin
        Epot = energy_potential(sys)
        @show Epot
        Eint = energy_internal(sys)
        @show Eint
        println("# of part. = ", length(sys.particles))
        println()
        save_frame!(out, sys, :v, :a, :P, :rho, :rho0)
    end
end

function main()
	sys = make_system()
	out = new_pvd_file("results/Kepler_vortex")

    #initialization
    create_cell_list!(sys)
    apply!(sys, find_rho0!, self = true)
    apply!(sys, find_rho!, self = true)
    apply!(sys, find_pressure!)
    apply!(sys, internal_force!)

	N_of_particles = length(sys.particles)
	@show(N_of_particles)
	@show(m)

	step_final = Int64(round(t_end/dt))
	for k = 0 : step_final
        verlet_step!(sys)
        save_results!(out, sys, k)
	end

	save_pvd_file(out)

end ## function main

end ## module


text:
include("Potts.jl")
include("/home/abhijith/Dropbox/Code/GRISE_code/Potts_learning.jl")
using IterTools:product

p =15
q = 2
L = 6

Terms = Dict{}()
#The interaction strengths used in the paper are given below.
#J = [0.2598910820966015  -0.653715969944616  0.07787582451255748 -0.3704026001241134 0.0397234544123 0.23441234]
#For now we will generate them at random
J =  0.15*randn(L)
#The larger the size of these numbers , more smaples will be required to leran the model. This can be seen from the IT bound
println(J)

for r in 1:L
    M = pth_order_tensor(r)
    for i in 1:p-L+1
        index =  Tuple(i:i+r-1)
        Terms[index] =  J[r]*M
    end
end

F = FactorGraph([1:L...],p,q,Terms);

#Sample from model and save at path
n_samples = 10^6

s = raw_sampler_potts(F, n_samples, false)
samples=[2*a .- 3  for a in s]
sample_path ="samples.csv"
writedlm( sample_path,hcat(samples...))
#println("exact", countmap(samples))

#Generating another set of samples for testing the trained model.

s = raw_sampler_potts(F, 10^6, false)
samples2=[2*a .- 3  for a in s]
sample_path ="samples_test.csv"
writedlm( sample_path,hcat(samples2...))


text:
"""
    diff_optimizer(optimizer_constructor)::Optimizer

Creates a `DiffOpt.Optimizer`, which is an MOI layer with an internal optimizer
and other utility methods. Results (primal, dual and slack values) are obtained
by querying the internal optimizer instantiated using the
`optimizer_constructor`. These values are required for find jacobians with respect to problem data.

One define a differentiable model by using any solver of choice. Example:

```julia
julia> import DiffOpt, HiGHS

julia> model = DiffOpt.diff_optimizer(HiGHS.Optimizer)
julia> x = model.add_variable(model)
julia> model.add_constraint(model, ...)
```
"""
function diff_optimizer(optimizer_constructor)::Optimizer
    optimizer = MOI.instantiate(optimizer_constructor, with_bridge_type=Float64)
    # When we do `MOI.copy_to(diff, optimizer)` we need to efficiently `MOI.get`
    # the model information from `optimizer`. However, 1) `optimizer` may not
    # implement some getters or it may be inefficient and 2) the getters may be
    # unimplemented or inefficient through some bridges.
    # For this reason we add a cache layer, the same cache JuMP adds.
    caching_opt = MOIU.CachingOptimizer(
        MOIU.UniversalFallback(MOIU.Model{Float64}()),
        optimizer,
    )
    return Optimizer(caching_opt)
end

mutable struct Optimizer{OT <: MOI.ModelLike} <: MOI.AbstractOptimizer
    optimizer::OT

    model_constructors::Vector{Any}
    model_constructor::Any

    diff::Any

    index_map::Union{Nothing,MOI.Utilities.IndexMap}

    # sensitivity input cache using MOI like sparse format
    input_cache::InputCache

    function Optimizer(optimizer::OT) where {OT <: MOI.ModelLike}
        output = new{OT}(
            optimizer,
            Any[],
            nothing,
            nothing,
            nothing,
            InputCache(),
        )
        add_all_model_constructors(output)
        return output
    end
end

"""
   add_model_constructor(optimizer::Optimizer, model_constructor)

Add the constructor of [`AbstractModel`](@ref) for `optimizer` to choose
from when trying to differentiate.
"""
function add_model_constructor(optimizer::Optimizer, model_constructor)
    push!(optimizer.model_constructors, model_constructor)
    return
end

function MOI.add_variable(model::Optimizer)
    model.diff = nothing
    vi = MOI.add_variable(model.optimizer)
    return vi
end

function MOI.add_variables(model::Optimizer, N::Int)
    model.diff = nothing
    return VI[MOI.add_variable(model) for i in 1:N]
end

function MOI.add_constraint(model::Optimizer, f::MOI.AbstractFunction, s::MOI.AbstractSet)
    model.diff = nothing
    return MOI.add_constraint(model.optimizer, f, s)
end

function MOI.add_constraints(model::Optimizer, f::AbstractVector{F}, s::AbstractVector{S}) where {F<:MOI.AbstractFunction, S<: MOI.AbstractSet}
    model.diff = nothing
    return CI{F, S}[MOI.add_constraint(model, f[i], s[i]) for i in eachindex(f)]
end

function MOI.set(model::Optimizer, attr::MOI.ObjectiveFunction{F}, f::F) where {F<:MOI.AbstractFunction}
    model.diff = nothing
    MOI.set(model.optimizer, attr, f)
end

MOI.supports(model::Optimizer, attr::MOI.ObjectiveSense) = MOI.supports(model.optimizer, attr)

function MOI.set(model::Optimizer, attr::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    model.diff = nothing
    return MOI.set(model.optimizer, attr, sense)
end

function MOI.get(model::Optimizer, attr::MOI.AbstractModelAttribute)
    return MOI.get(model.optimizer, attr)
end

function MOI.get(model::Optimizer, attr::MOI.ListOfConstraintIndices)
    return MOI.get(model.optimizer, attr)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintSet, ci::CI)
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(
    model::Optimizer, attr::MOI.ConstraintSet, ci::CI{F, S}, s::S
) where {F,S}
    model.diff = nothing
    return MOI.set(model.optimizer, attr, ci, s)
end

function MOI.get(
    model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}
) where {F,S}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.set(
    model::Optimizer, attr::MOI.ConstraintFunction, ci::CI{F, S}, f::F
) where {F,S}
    model.diff = nothing
    return MOI.set(model.optimizer, attr, ci, f)
end

# `MOI.supports` methods

function MOI.supports(model::Optimizer, attr::MOI.AbstractModelAttribute)
    return MOI.supports(model.optimizer, attr)
end

function MOI.supports(model::Optimizer, attr::MOI.ObjectiveFunction)
    return MOI.supports(model.optimizer, attr)
end

function MOI.supports_constraint(
    model::Optimizer, ::Type{F}, ::Type{S},
) where {F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
    MOI.supports_constraint(model.optimizer, F, S)
end
function MOI.supports(
    model::Optimizer, attr::MOI.ConstraintName, ::Type{CI{F, S}}
) where {F,S}
    return MOI.supports(model.optimizer, attr, CI{F,S})
end

MOI.get(model::Optimizer, attr::MOI.SolveTimeSec) = MOI.get(model.optimizer, attr)

function MOI.empty!(model::Optimizer)
    MOI.empty!(model.optimizer)
    model.diff = nothing
    model.index_map = nothing
    empty!(model.input_cache)
    return
end

function MOI.is_empty(model::Optimizer)
    return MOI.is_empty(model.optimizer) &&
           model.diff === nothing
end

function MOI.supports_incremental_interface(model::Optimizer)
    if !MOI.supports_incremental_interface(model.optimizer)
        error(
            "DiffOpt requires a solver that " *
            "`MOI.supports_incremental_interface`, which is not the case for " *
            "$(model.optimizer)")
    end
    return true
end

function MOI.copy_to(model::Optimizer, src::MOI.ModelLike)
    model.diff = nothing
    return MOIU.default_copy_to(model.optimizer, src)
end

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    return MOI.get(model.optimizer, MOI.TerminationStatus())
end

function MOI.set(model::Optimizer, ::MOI.VariablePrimalStart,
                 vi::VI, value::Float64)
    model.diff = nothing
    MOI.set(model.optimizer, MOI.VariablePrimalStart(), vi, value)
end

function MOI.supports(model::Optimizer, attr::MOI.AbstractVariableAttribute,
                      ::Type{MOI.VariableIndex})
    return MOI.supports(model.optimizer, attr, MOI.VariableIndex)
end

function MOI.set(model::Optimizer, attr::MOI.AbstractVariableAttribute, v::MOI.VariableIndex, value)
    MOI.set(model.optimizer, attr, v, value)
end

function MOI.get(model::Optimizer, attr::MOI.AbstractVariableAttribute, vi::MOI.VariableIndex)
    return MOI.get(model.optimizer, attr, vi)
end

function MOI.delete(model::Optimizer, ci::CI{F,S}) where {F,S}
    model.diff = nothing
    MOI.delete(model.optimizer, ci)
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal, ci::CI{F,S}) where {F,S}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.is_valid(model::Optimizer, v::VI)
    return MOI.is_valid(model.optimizer, v::VI)
end

MOI.is_valid(model::Optimizer, con::CI) = MOI.is_valid(model.optimizer, con)

function MOI.get(model::Optimizer, attr::MOI.ConstraintDual, ci::CI{F,S}) where {F,S}
    return MOI.get(model.optimizer, attr, ci)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintBasisStatus, ci::CI{F,S}) where {F,S}
    return MOI.get(model.optimizer, MOI.ConstraintBasisStatus(), ci)
end

# helper methods to check if a constraint contains a Variable
function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.VariableIndex})
    return v == MOI.get(model, MOI.ConstraintFunction(), ci)
end

function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.ScalarAffineFunction{Float64}})
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return any(term -> v == term.variable, func.terms)
end

function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.VectorOfVariables})
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return v in func.variables
end

function _constraint_contains(model::Optimizer, v::VI, ci::CI{MOI.VectorAffineFunction{Float64}})
    func = MOI.get(model, MOI.ConstraintFunction(), ci)
    return any(term -> v == term.scalar_term.variable, func.terms)
end

function MOI.delete(model::Optimizer, v::VI)
    model.diff = nothing
    MOI.delete(model.optimizer, v)
end

# for array deletion
function MOI.delete(model::Optimizer, indices::Vector{VI})
    model.diff = nothing
    for i in indices
        MOI.delete(model, i)
    end
end

function MOI.modify(
    model::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
    chg::MOI.AbstractFunctionModification
)
    model.diff = nothing
    MOI.modify(
        model.optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        chg
    )
end

function MOI.modify(model::Optimizer, ci::CI, chg::MOI.AbstractFunctionModification)
    model.diff = nothing
    MOI.modify(model.optimizer, ci, chg)
end

function MOI.get(model::Optimizer, ::Type{MOI.VariableIndex}, name::String)
    return MOI.get(model.optimizer, MOI.VariableIndex, name)
end

function MOI.set(model::Optimizer, ::MOI.ConstraintName, con::CI, name::String)
    MOI.set(model.optimizer, MOI.ConstraintName(), con, name)
end

function MOI.get(model::Optimizer, ::Type{MOI.ConstraintIndex}, name::String)
    return MOI.get(model.optimizer, MOI.ConstraintIndex, name)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintName, con::CI)
    return MOI.get(model.optimizer, MOI.ConstraintName(), con)
end

function MOI.get(model::Optimizer, ::MOI.ConstraintName, ::Type{CI{F, S}}) where {F,S}
    return MOI.get(model.optimizer, MOI.ConstraintName(), CI{F,S})
end

function MOI.set(model::Optimizer, ::MOI.Name, name::String)
    MOI.set(model.optimizer, MOI.Name(), name)
end

function MOI.get(model::Optimizer, ::Type{CI{F, S}}, name::String) where {F,S}
    return MOI.get(model.optimizer, CI{F,S}, name)
end

function MOI.supports(model::Optimizer, attr::MOI.TimeLimitSec)
    MOI.supports(model.optimizer, attr)
end
function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value::Union{Real, Nothing})
    MOI.set(model.optimizer, MOI.TimeLimitSec(), value)
end
function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    return MOI.get(model.optimizer, MOI.TimeLimitSec())
end

MOI.supports(model::Optimizer, ::MOI.Silent) = MOI.supports(model.optimizer, MOI.Silent())

function MOI.set(model::Optimizer, ::MOI.Silent, value)
    MOI.set(model.optimizer, MOI.Silent(), value)
end

function MOI.get(model::Optimizer, ::MOI.Silent)
    return MOI.get(model.optimizer, MOI.Silent())
end

function MOI.get(model::Optimizer, attr::A) where {A <: Union{MOI.SolverName, MOI.SolverVersion}}
    return MOI.get(model.optimizer, attr)
end

function MOI.optimize!(model::Optimizer)
    model.diff = nothing
    MOI.optimize!(model.optimizer)
    return
end

"""
    ModelConstructor <: MOI.AbstractOptimizerAttribute

Determines which subtype of [`DiffOpt.AbstractModel`](@ref) to use for
differentiation. When set to `nothing`, the first one out of
`model.model_constructors` that support the problem is used.
"""
struct ModelConstructor <: MOI.AbstractOptimizerAttribute end

MOI.supports(::Optimizer, ::ModelConstructor) = true
MOI.get(model::Optimizer, ::ModelConstructor) = model.model_constructor
function MOI.set(model::Optimizer, ::ModelConstructor, model_constructor)
    model.model_constructor = model_constructor
end

function reverse_differentiate!(model::Optimizer)
    st = MOI.get(model.optimizer, MOI.TerminationStatus())
    if !in(st,  (MOI.LOCALLY_SOLVED, MOI.OPTIMAL))
        error("Trying to compute the reverse differentiation on a model with termination status $(st)")
    end
    diff = _diff(model)
    for (vi, value) in model.input_cache.dx
        MOI.set(diff, ReverseVariablePrimal(), model.index_map[vi], value)
    end
    return reverse_differentiate!(diff)
end

function _copy_forward_in_constraint(diff, index_map, con_map, constraints)
    for (index, value) in constraints
        MOI.set(diff, ForwardConstraintFunction(), con_map[index], MOI.Utilities.map_indices(index_map, value))
    end
end

function forward_differentiate!(model::Optimizer)
    st = MOI.get(model.optimizer, MOI.TerminationStatus())
    if !in(st,  (MOI.LOCALLY_SOLVED, MOI.OPTIMAL))
        error("Trying to compute the forward differentiation on a model with termination status $(st)")
    end
    diff = _diff(model)
    if model.input_cache.objective !== nothing
        MOI.set(diff, ForwardObjectiveFunction(), MOI.Utilities.map_indices(model.index_map, model.input_cache.objective))
    end
    for (F, S) in keys(model.input_cache.scalar_constraints.dict)
        _copy_forward_in_constraint(diff, model.index_map, model.index_map.con_map[F, S], model.input_cache.scalar_constraints[F, S])
    end
    for (F, S) in keys(model.input_cache.vector_constraints.dict)
        _copy_forward_in_constraint(diff, model.index_map, model.index_map.con_map[F, S], model.input_cache.vector_constraints[F, S])
    end
    forward_differentiate!(diff)
end

function _copy_constraint_start(dest, src, index_map::MOIU.DoubleDicts.IndexDoubleDictInner{F, S}, dest_attr, src_attr) where {F, S}
    for ci in MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
        value = MOI.get(src, src_attr, ci)
        MOI.set(dest, dest_attr, index_map[ci], value)
    end
end

function _instantiate_with_bridges(model_constructor)
    model = MOI.Bridges.LazyBridgeOptimizer(
        MOI.instantiate(model_constructor)
    )
    # We don't add any variable bridge here because:
    # 1) If `ZerosBridge` is used, `MOI.Bridges.unbridged_function` does not work.
    #    This is in fact expected: since `ZerosBridge` drops the variable, we dont
    #    compute the derivative of the value of this variable as a function of its fixed value.
    #    This could be easily determined as the same as the derivative of the value but
    #    since the variable was also dropped from other constraints, we would ignore its impact on the other constraints.
    # 2) For affine variable bridges, `bridged_function` and `unbridged_function` don't treat the function as a derivative hence they will add constants
    MOI.Bridges.Constraint.add_all_bridges(model, Float64)
    MOI.Bridges.Objective.add_all_bridges(model, Float64)
    return model
end

function _diff(model::Optimizer)
    if model.diff === nothing
        _check_termination_status(model)
        model_constructor = MOI.get(model, ModelConstructor())
        if isnothing(model_constructor)
            model.diff = nothing
            for constructor in model.model_constructors
                model.diff = _instantiate_with_bridges(constructor)
                try
                    model.index_map = MOI.copy_to(model.diff, model.optimizer)
                catch err
                    if err isa MOI.UnsupportedConstraint || err isa MOI.UnsupportedAttribute
                        model.diff = nothing
                    else
                        rethrow(err)
                    end
                end
                if !isnothing(model.diff)
                    break
                end
            end
            if isnothing(model.diff)
                error("No differentiation model supports the problem. If you believe it should be supported, say by `DiffOpt.QuadraticProgram.Model`, use `MOI.set(model, DiffOpt.ModelConstructor, DiffOpt.QuadraticProgram.Model)` and try again to see an error indicating why it is not supported.")
            end
        else
            model.diff = _instantiate_with_bridges(model_constructor)
            model.index_map = MOI.copy_to(model.diff, model.optimizer)
        end
        vis_src = MOI.get(model.optimizer, MOI.ListOfVariableIndices())
        MOI.set(model.diff, MOI.VariablePrimalStart(), getindex.(Ref(model.index_map), vis_src), MOI.get(model.optimizer, MOI.VariablePrimal(), vis_src))
        for (F, S) in MOI.get(model.optimizer, MOI.ListOfConstraintTypesPresent())
            _copy_constraint_start(model.diff, model.optimizer, model.index_map.con_map[F, S], MOI.ConstraintPrimalStart(), MOI.ConstraintPrimal())
            _copy_constraint_start(model.diff, model.optimizer, model.index_map.con_map[F, S], MOI.ConstraintDualStart(), MOI.ConstraintDual())
        end
    end
    return model.diff
end

function _check_termination_status(model::Optimizer)
    if !in(
        MOI.get(model, MOI.TerminationStatus()), (MOI.LOCALLY_SOLVED, MOI.OPTIMAL)
        )
        error("problem status: ", MOI.get(model.optimizer, MOI.TerminationStatus()))
    end
end

# DiffOpt attributes redirected to `diff`

function _checked_diff(model::Optimizer, ::MOI.AnyAttribute, call)
    if model.diff === nothing
        error("Cannot get attribute `attr`. First call `DiffOpt.$call`.")
    end
    return model.diff
end

function MOI.get(model::Optimizer, attr::ReverseObjectiveFunction)
    return IndexMappedFunction(
        MOI.get(_checked_diff(model, attr, :reverse), attr),
        model.index_map,
    )
end
MOI.supports(::Optimizer, ::ForwardObjectiveFunction) = true
function MOI.get(model::Optimizer, ::ForwardObjectiveFunction)
    return model.input_cache.objective
end
function MOI.set(model::Optimizer, ::ForwardObjectiveFunction, objective)
    model.input_cache.objective = objective
    return
end

function MOI.get(model::Optimizer, attr::ForwardVariablePrimal, vi::MOI.VariableIndex)
    return MOI.get(_checked_diff(model, attr, :forward), attr, model.index_map[vi])
end
MOI.supports(::Optimizer, ::ReverseVariablePrimal, ::Type{VI}) = true
function MOI.get(model::Optimizer, ::ReverseVariablePrimal, vi::VI)
    return get(model.input_cache.dx, vi, 0.0)
end
function MOI.set(model::Optimizer, ::ReverseVariablePrimal, vi::VI, val)
    model.input_cache.dx[vi] = val
    return
end

function MOI.get(model::Optimizer, attr::ReverseConstraintFunction, ci::MOI.ConstraintIndex)
    return IndexMappedFunction(
        MOI.get(_checked_diff(model, attr, :reverse), attr, model.index_map[ci]),
        model.index_map,
    )
end
function MOI.get(model::Optimizer,
    ::ForwardConstraintFunction, ci::CI{MOI.ScalarAffineFunction{T},S}
) where {T,S}
    return get(model.input_cache.scalar_constraints, ci, zero(MOI.ScalarAffineFunction{T}))
end
function MOI.supports(
    ::Optimizer,
    ::ForwardConstraintFunction,
    ::Type{CI{F,S}}
) where {F<:Union{MOI.ScalarAffineFunction,MOI.VectorAffineFunction},S}
    return true
end
function MOI.get(model::Optimizer,
    ::ForwardConstraintFunction, ci::CI{MOI.VectorAffineFunction{T},S}
) where {T,S}
    func = get(model.input_cache.vector_constraints, ci, nothing)
    if func === nothing
        set = MOI.get(model, MOI.ConstraintSet(), ci)
        dim = MOI.dimension(set)
        return MOI.Utilities.zero_with_output_dimension(MOI.VectorAffineFunction{T}, dim)
    else
        return func
    end
end
function MOI.set(model::Optimizer,
    ::ForwardConstraintFunction,
    ci::CI{MOI.ScalarAffineFunction{T},S},
    func::MOI.ScalarAffineFunction{T},
) where {T,S}
    model.input_cache.scalar_constraints[ci] = func
    return
end
function MOI.set(model::Optimizer,
    ::ForwardConstraintFunction,
    ci::CI{MOI.VectorAffineFunction{T},S},
    func::MOI.VectorAffineFunction{T},
) where {T,S}
    model.input_cache.vector_constraints[ci] = func
    return
end


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# type to contain a probability distribution sampler and corresponding random number generator
mutable struct DistrRng{T<:Sampleable}
	d::T
	rng::MersenneTwister

	function DistrRng(d::T, rng::MersenneTwister) where T <: Sampleable
		return new{T}(d, deepcopy(rng))
	end
	function DistrRng(d::T; seed::Int = nullIndex) where T <: Sampleable
		rng = (seed >= 0 ? MersenneTwister(seed) : MersenneTwister(rand(UInt32)))
		return new{T}(d, rng)
	end
end

global GlobalRngBackup = MersenneTwister(0); # global variable to sometimes store GLOBAL_RNG state

function copyRng!(dest::MersenneTwister, src::MersenneTwister)
	# faster than using the copy! function
	dest.seed = src.seed
	dest.state = src.state
	dest.vals = src.vals
	dest.ints = src.ints
	dest.idxF = src.idxF
	dest.idxI = src.idxI
	# # slower, by several times:
	# for fname in [:seed, :state, :vals, :ints, :idxF, :idxI] # = fieldnames(MersenneTwister)
		# setfield!(dest, fname, getfield(src, fname))
	# end
end

# Generate next random value from distribution in distrRng, with RNG in distrRng.
# Base.rand(::AbstractRNG, ::Distribution) does not exist for all distributions, so have to use
# Base.rand(::Distribution) and after setting GLOBAL_RNG.
function Base.rand(distrRng::DistrRng, n::Int)
	if isa(GLOBAL_RNG, MersenneTwister) # julia <= v1.2
		rng = distrRng.rng # shorthand
		# store GLOBAL_RNG state in backup, set GLOBAL_RNG to rng
		copyRng!(GlobalRngBackup, GLOBAL_RNG)
		copyRng!(GLOBAL_RNG, rng)
		value = rand(distrRng.d, n)
		# set rng to state of GLOBAL_RNG, restore GLOBAL_RNG from backup
		copyRng!(rng, GLOBAL_RNG)
		copyRng!(GLOBAL_RNG, GlobalRngBackup) # for safety, so rng is no longer tied to GLOBAL_RNG
		return value
	else # julia v1.3
		return rand(distrRng.rng, distrRng.d, n)
	end
end
function Base.rand(distrRng::DistrRng)
	return rand(distrRng, 1)[1]
end


text:
using Base: @deprecate_binding

# Types ...f0 renamed to ...f
@deprecate_binding Vecf0 Vecf
@deprecate_binding Pointf0 Pointf
for i in 1:4
    for T in [:Point, :Vec]
        oldname = Symbol("$T$(i)f0")
        newname = Symbol("$T$(i)f")
        @eval begin
            @deprecate_binding $oldname $newname
        end
    end
    oldname = Symbol("Mat$(i)f0")
    newname = Symbol("Mat$(i)f")
    @eval begin
        @deprecate_binding $oldname $newname
    end
end

# Rect types
@deprecate_binding Rect2D  Rect2
@deprecate_binding Rect3D  Rect3
@deprecate_binding FRect   Rectf
@deprecate_binding FRect2D Rect2f
@deprecate_binding FRect3D Rect3f
@deprecate_binding IRect   Recti
@deprecate_binding IRect2D Rect2i
@deprecate_binding IRect3D Rect3i
@deprecate_binding TRect   RectT


text:
# ....................
# @license: %MIT License%:~ http://www.opensource.org/licenses/MIT
# @project: julia
# @file: /try.jl
# @created: Thursday, 26th March 2020
# @author: brentian (chuwzhang@126.com)
# @modified: brentian (chuwzhang@126.com>)
#    Thursday, 26th March 2020 3:31:03 pm

# @description:
__precompile__()
module lr_with_outlier

using JuMP
using Gurobi
using Plots
using LinearAlgebra
using Printf

# export module functions
export f

function f(x, y0, n, m, N, x_min, x_max, y_min, y_max)
    model = Model()
    # Declaring variables
    @variable(model, b[1:m] >= 0)
    # @variable(model, 0 <= q[1:n] <= 1)
    @variable(model, q[1:n], Bin)

    @variable(model, res[1:n])
    # == absolute value ==
    @constraint(model, RES[i = 1:n], res[i] + dot(b, x[i, :]) >= y0[i])
    @constraint(model, RES_NEG[i = 1:n], res[i] - dot(b, x[i, :]) >= -y0[i])
    # == quadratic ==
    # @constraint(model,
    #   RES_QUAD[i = 1:n],
    #   [res[i], dot(b, x[i,:]) - y0[i]] in SecondOrderCone())
    block_m = [
        1 * diagm(ones(n)) 0.5 * diagm(ones(n))
        0.5 * diagm(ones(n)) 1 * diagm(ones(n))
    ]
    @constraint(model, sum(q) >= N)
    @objective(model, Min, dot([res; q], block_m, [res; q]))
    set_optimizer(model, Gurobi.Optimizer)

    # Printing the prepared optimization model
    print(model)

    # Solving the optimization problem
    optimize!(model)

    print(value.(q), '\n')
    print(value.(b), '\n')

    # === summarizing results ===
    q1 = value.(q)

    # metrics
    y_bar = sum(y0) / length(y0)
    s_tot = sum((y0 .- y_bar).^2 .* q1)
    r = value.(res)
    s_res = sum(r.^2 .* q1)

    # make inferences
    y_pred = [dot(value.(b), x[i, :]) * (y_max - y_min) + y_min for i in 1:n]

    selected_index = [i for i = 1:n if q1[i] == 1]
    dropped_index = [i for i = 1:n if q1[i] != 1]
    s_x = x[selected_index] .* (x_max - x_min) .+ x_min
    d_x = x[dropped_index] .* (x_max - x_min) .+ x_min
    s_y0 = y0[selected_index] .* (y_max - y_min) .+ y_min
    d_y0 = y0[dropped_index] .* (y_max - y_min) .+ y_min
    s_y_pred = y_pred[selected_index]
    d_y_pred = y_pred[dropped_index]
    r_sqr = (1 - s_res / s_tot)
    p = plot(
        s_x,
        s_y_pred,
        lab = ["selected_regression"],
        title = "drop $(n - N) =>\nR^2: $r_sqr",
        legend = false,
    )
    scatter!(s_x, s_y0, lab = ["selected_samples"])
    scatter!(d_x, d_y0, lab = ["outliers"])
    return [p, value.(b), N, q1, r_sqr, s_x, s_y_pred]
end

end

max_stars_repo_name: brentian/artoopt

max_issues_repo_name: brentian/artoopt

max_forks_repo_name: brentian/artoopt