
####################################################################################################################
#### Load packages ####
using Pkg
Pkg.activate(".")
Pkg.add("JuMP")
Pkg.add("Combinatorics")
Pkg.add("DataFrames")
Pkg.add("CSV")
Pkg.add("ExcelFiles")
Pkg.add("CSVFiles")
Pkg.add("DelimitedFiles")
Pkg.add("XLSX")
Pkg.add("BenchmarkTools")
Pkg.add("GLPK")
Pkg.add("Plots")
Pkg.add("StatsPlots")
Pkg.add(Pkg.PackageSpec(name = "Gurobi", version = v"0.9.4"))

## Precompile project
Pkg.precompile()

using JuMP
using Combinatorics
using DataFrames
using CSV
using ExcelFiles
using DelimitedFiles
using XLSX
using BenchmarkTools
using GLPK
using Gurobi




####################################################################################################################
# #### Data Prep ####
model_first = 0;
model_last = 4;

for i in (model_first:model_last)

    println("Scenario $i initiated")


    ##### Sets #####
    # read sets
    sets_df = load("Scenario$i.xlsx", "sets") |> DataFrame;
    sets_df;
    set_K = sets_df[!,:K]; 
    set_S = sets_df[!,:S];
    set_L = sets_df[!,:L];
    set_N = sets_df[!,:N];

    set_K = findall(!ismissing, set_K);
    set_S = findall(!ismissing, set_S);
    set_L = findall(!ismissing, set_L);
    set_N = findall(!ismissing, set_N);

    set_A = [(i,j) for (i,j) in permutations(set_N,2) if i != j];
    set_IL = [(i,l) for i in set_N for l in set_L];
    set_KN = [(k,i) for k in set_K for i in set_N];
    set_KSA = [(k,s,i,j) for k in set_K for s in set_S for (i,j) in set_A];
    set_KSN = [(k,s,i) for k in set_K for s in set_S for i in set_N];

    ##### Parameters with multiple indices #####
    # Read from input data file
    param_c_df = load("Scenario$i.xlsx", "param_cc") |> DataFrame;
    # Convert and create parameters
    param_c_df[!,[:k, :s, :i, :j]] = convert.(Int,param_c_df[!,[:k, :s, :i, :j]]) #convert to int
    param_c_df[!,:index] = zip(param_c_df[!,:k], param_c_df[!,:s],param_c_df[!,:i],param_c_df[!,:j]) |> collect #create index col
    param_c_df = DataFrames.select!(param_c_df,[:index,:value]);
    param_c = eachrow(param_c_df) |> Dict # convert to dict type

    param_U_df = load("Scenario$i.xlsx", "param_U") |> DataFrame;
    param_U_df[!,[:s, :i, :j]] = convert.(Int,param_U_df[!,[:s, :i, :j]]);
    param_U_df[!,:index] = zip(param_U_df[!,:s],param_U_df[!,:i],param_U_df[!,:j]) |> collect;
    param_U_df = DataFrames.select!(param_U_df,[:index,:value]);
    param_U = eachrow(param_U_df) |> Dict;

    param_v_df = load("Scenario$i.xlsx", "param_vv") |> DataFrame;
    param_v_df[!,[:k,:s,:i]] = convert.(Int,param_v_df[!,[:k,:s,:i]]);
    param_v_df[!,:index] = zip(param_v_df[!,:k], param_v_df[!,:s],param_v_df[!,:i]) |> collect;
    param_v_df = DataFrames.select!(param_v_df,[:index,:value]);
    param_v = eachrow(param_v_df) |> Dict;

    param_prop_df = load("Scenario$i.xlsx", "param_prop") |> DataFrame;
    param_prop_df[!,[:k,:s,:i]] = convert.(Int,param_prop_df[!,[:k,:s,:i]]);
    param_prop_df[!,:index] = zip(param_prop_df[!,:k], param_prop_df[!,:s],param_prop_df[!,:i]) |> collect; 
    param_prop_df = DataFrames.select!(param_prop_df,[:index,:value]);
    param_prop = eachrow(param_prop_df) |> Dict;

    param_F_df = load("Scenario$i.xlsx", "param_F") |> DataFrame;
    param_F_df[!,[:i,:l]] = convert.(Int,param_F_df[!,[:i,:l]]);
    param_F_df[!,:index] = zip(param_F_df[!,:i], param_F_df[!,:l]) |> collect;
    param_F_df = DataFrames.select!(param_F_df,[:index,:value]); 
    param_F = eachrow(param_F_df) |> Dict;

    param_d_df = load("Scenario$i.xlsx", "param_dd") |> DataFrame;
    param_d_df[!,[:i,:j]] = convert.(Int,param_d_df[!,[:i,:j]]);
    param_d_df[!,:index] = zip(param_d_df[!,:i], param_d_df[!,:j]) |> collect;
    param_d_df = DataFrames.select!(param_d_df,[:index,:value]);
    param_d = eachrow(param_d_df) |> Dict;

    ##### Parameters with single index #####
    param_M_df = load("Scenario$i.xlsx", "param_M") |> DataFrame;
    param_M_df[!,[:index]] = convert.(Int,param_M_df[!,[:index]]);
    param_M = eachrow(param_M_df) |> Dict;

    param_b_df = load("Scenario$i.xlsx", "param_bb") |> DataFrame;
    param_b_df[!,[:index]] = convert.(Int,param_b_df[!,[:index]]);
    param_b = eachrow(param_b_df) |> Dict;

    param_q_df = load("Scenario$i.xlsx", "param_qq") |> DataFrame;
    param_q_df[!,[:index]] = convert.(Int,param_q_df[!,[:index]]);
    param_q = eachrow(param_q_df) |> Dict;

    param_h_df = load("Scenario$i.xlsx", "param_hh") |> DataFrame;
    param_h_df[!,[:index]] = convert.(Int,param_h_df[!,[:index]]);
    param_h = eachrow(param_h_df) |> Dict;

    param_p_df = load("Scenario$i.xlsx", "param_pp") |> DataFrame;
    param_p_df[!,[:index]] = convert.(Int,param_p_df[!,[:index]]);
    param_p = eachrow(param_p_df) |> Dict;

    param_P_df = load("Scenario$i.xlsx", "param_P") |> DataFrame;
    param_P_df[!,[:index]] = convert.(Int,param_P_df[!,[:index]]);
    param_P = eachrow(param_P_df) |> Dict;

    param_u_df = load("Scenario$i.xlsx", "param_uu") |> DataFrame;
    param_u_df[!,[:index]] = convert.(Int,param_u_df[!,[:index]]);
    param_u = eachrow(param_u_df) |> Dict;

    param_bigM_df = load("Scenario$i.xlsx", "param_bigM") |> DataFrame;
    param_bigM_df[!,[:index]] = convert.(Int,param_bigM_df[!,[:index]]);
    param_bigM = eachrow(param_bigM_df) |> Dict;

    # Adjusting parameters
    param_D_df = load("Scenario$i.xlsx", "param_D") |> DataFrame;
    param_D_df[!,[:index]] = convert.(Int,param_D_df[!,[:index]]);
    param_D = eachrow(param_D_df) |> Dict;

    param_alpha_df = load("Scenario$i.xlsx", "param_alpha") |> DataFrame;
    param_alpha_df[!,[:index]] = convert.(Int,param_alpha_df[!,[:index]]);
    param_alpha = eachrow(param_alpha_df) |> Dict;

    println("Parameters set up successfully for scenario $i")
    ####################################################################################################################
    #### Model Prep ####
    # Time the model
    TimeTrack_start = time();

    # Preparing an optimization model
    m = Model(Gurobi.Optimizer);
    # m = Model(GLPK.Optimizer);

    ##Declaring variables

    # variables
    @variable(m, y[index in set_IL], Bin);
    @variable(m, r[index in set_KN] >= 0,Int);
    @variable(m, x[index in set_KSA] >= 0,Int);
    @variable(m, z[index in set_KSN] >= 0,Int);
    @variable(m, w[index in set_KSN] >= 0,Int);
    @variable(m, gamma[s in set_S], Bin);


    ##Objective function
    @objective(m, Min, sum(param_F[(i,l)] * y[(i,l)] for (i,l) in set_IL) # Fix Cost
                        + sum(param_q[k] * r[(k,n)] for (k,n) in set_KN)  # Preposition Cost
                        + sum(param_P[s] * param_c[(k,s,i,j)] * x[(k,s,i,j)] for (k,s,i,j) in set_KSA) # Scenario Specific Cost 1
                        + sum(param_P[s] * param_h[k] * z[(k,s,i)] + param_P[s] * param_p[k] * w[(k,s,i)] for (k,s,i) in set_KSN) # Scenario Specific Cost 2
                ); 

    ##Constraints 
    # Flow Conservation
    fc = @constraint(m, fc[k in set_K, s in set_S, i in set_N],
                            sum(x[(k,s,j,i)] for j in set_N if j != i) + param_prop[(k,s,i)]* r[(k,i)] - z[(k,s,i)]
                            ==
                            sum(x[(k,s,i,j)] for j in set_N if j != i) + param_v[(k,s,i)] - w[(k,s,i)]
                );

    # Open Facility Capacity
    ofc = @constraint(m, ofc[i in set_N],
                        sum(param_b[k] * r[(k,i)] for k in set_K)
                        <= 
                        sum(param_M[l] * y[(i,l)] for l in set_L)
                );
    # Facility Number per Node
    fpn = @constraint(m, fpn[i in set_N],
                    sum(y[(i,l)] for l in set_L)
                    <=
                    1
                );
    # Arc Capacity
    ac = @constraint(m, ac[s in set_S, (i,j) in set_A],
                    sum(param_u[k] * x[(k,s,i,j)] for k in set_K)
                    <=
                    param_U[(s,i,j)]
                );


    # Average Distance Limit
    adl = @constraint(m, adl[k in set_K, s in set_S],
                        sum(param_d[(i,j)] * x[(k,s,i,j)] for (i,j) in set_A)
                        <=
                        param_D[k] * sum(param_v[(k,s,i)] for i in set_N) + param_bigM[1] * (1 - gamma[s])            
                );

    # Reliability Set Definition
    rsd = @constraint(m, rsd,
                    sum(param_P[s] * gamma[s] for s in set_S)
                    >=
                    param_alpha[1]
                );

    # Demand requirements for scenarios in reliable set
    dmrs = @constraint(m, dmrs[k in set_K, s in set_S, i in set_N],
                        w[(k,s,i)]
                        <=
                        param_v[(k,s,i)] * (1 - gamma[s])               
                );

    println("Model successfully set up for scenario $i")
    ####################################################################################################################
    #### Model solving ####

    # Printing the prepared optimization model
    # print(m);

    # Solving the optimization problem
    JuMP.optimize!(m);
    TimeTrack_elapsed = time() - TimeTrack_start;

    ####################################################################################################################
    #### Post-processing ####

    # Print result
    println("Model running time is: ", round(TimeTrack_elapsed,digits = 5),"s")
    println("Model solving time:",round(JuMP.solve_time(m),digits = 5), "s")

    # save the optimized variables to *_star
    y_star = JuMP.value.(y);
    r_star = JuMP.value.(r);
    x_star = JuMP.value.(x);
    z_star = JuMP.value.(z);
    w_star = JuMP.value.(w);
    gamma_star = JuMP.value.(gamma);

    # Print objective value
    println("Objective value: ", JuMP.objective_value(m))

    # Print the non-zero variables
    println("If a variable is not printed, the value is 0!")
    for (i,l) in set_IL
        if y_star[(i,l)] == 1
            println("facility $l located in node $i")
        else end
    end;

    for (k,i) in set_KN
        if r_star[(k,i)] > 0
            println("Amount of commodity $k pre-position at node $i: ", r_star[(k,i)])
        else end
    end;

    for (k,s,i,j) in set_KSA
        if x_star[(k,s,i,j)] > 0 
        println("Amount of commodity $k shipped across link ($i,$j) in scenario $s: ", x_star[(k,s,i,j)])
        else end
    end;

    for (k,s,i) in set_KSN
        if z_star[(k,s,i)] > 0 
            println("Commodity $k not used in scenario $s at node $i: ", z_star[(k,s,i)] )
        else end
    end;

    for (k,s,i) in set_KSN
        if w_star[(k,s,i)] > 0
            println("Shortage of commodity $k at node $i in scenario $s: ", w_star[(k,s,i)])
        else end
    end;

    for s in set_S
        if gamma_star[s] > 0
            println("Scenario $s belong to the reliable set")
        else end
    end;


    ####################################################################################################################
    #### Export Results ####

    # Write result to data frames
    y_df = DataFrame(variable_type =[], index=[],value=[]);
    for (i,l) in set_IL
        val = y_star[(i,l)]
        push!(y_df, vcat("y",(i,l),val))
    end;


    r_df = DataFrame(variable_type =[],index=[],value=[]);
    for (k,i) in set_KN
        val = r_star[(k,i)]
        push!(r_df, vcat("r",(k,i),val))
    end;

    x_df = DataFrame(variable_type =[],index=[],value=[]);
    for (k,s,i,j) in set_KSA
        val = x_star[(k,s,i,j)]
        push!(x_df, vcat("x",(k,s,i,j),val))
    end;

    z_df = DataFrame(variable_type =[],index=[],value=[]);
    for (k,s,i) in set_KSN
        val = z_star[(k,s,i)]
        push!(z_df, vcat("z",(k,s,i),val))
    end;

    w_df = DataFrame(variable_type =[],index=[],value=[]);
    for (k,s,i) in set_KSN
        val = w_star[(k,s,i)]
        push!(w_df, vcat("w",(k,s,i),val))
    end;

    gamma_df = DataFrame(variable_type =[],index=[],value=[]);
    for s in set_S
        val = gamma_star[s]
        push!(gamma_df, vcat("gamma",s,val))
    end;

    # Merge all the data frames together
    result_df = DataFrame(variable_type=[], index=[], value=[]);
    result_df = DataFrames.append!(result_df,y_df);
    result_df = DataFrames.append!(result_df,r_df);
    result_df = DataFrames.append!(result_df,x_df);
    result_df = DataFrames.append!(result_df,z_df);
    result_df = DataFrames.append!(result_df,w_df);
    result_df = DataFrames.append!(result_df,gamma_df);

    # Export the merged data frame
    CSV.write("output_var$i.csv",result_df);

    # Export Objective value and Running Time (Add another column for model no)
    vars_no = length(y_star) + length(r_star) + length(x_star) + length(z_star) + length(w_star)
    constraint_no = length(fc) + length(ofc) + length(fpn) + length(ac) + length(adl) + 1 + length(dmrs) 
    stat_df = DataFrame(scenario=[], objective_value = [], modeling_time = [], solve_time = [], 
                        vars_no = [], y_no = [], r_no = [] , x_no = [], z_no =[] , w_no = [],
                        constraint_no = [], fc_no = [], ofc_no = [], fpn_no = [],ac_no = [], adl_no = [], rsd_no = [], dmrs_no = []
                        );
    DataFrames.push!(stat_df, vcat(i, JuMP.objective_value(m), round(TimeTrack_elapsed,digits = 5), round(JuMP.solve_time(m),digits = 5),
    vars_no,length(y_star), length(r_star), length(x_star), length(z_star), length(w_star),
    constraint_no,length(fc),length(ofc),length(fpn),length(ac),length(adl),1,length(dmrs)
    ));

    CSV.write("output_stat$i.csv",stat_df);

    println("Results successfully exported for scenario $i");



end