# SPDX-FileCopyrightText: 2026 Yating Zhang
# SPDX-License-Identifier: MIT

"""
numerical 6.1 compare four formulations 
name: CompareModel 
""" 

# system stay awake function 
function enable_stay_awake()
    if Sys.iswindows()
        # Windows: Call Kernel32 API
        # ES_CONTINUOUS (0x80000000) | ES_SYSTEM_REQUIRED (0x00000001)
        # Tells the system that this thread requires the system to be in the working state.
        try
            ccall((:SetThreadExecutionState, "kernel32"), UInt32, (UInt32,), 0x80000001)
            println(">> Windows: System sleep prevention enabled.")
        catch
            println(">> Windows: Failed to set sleep prevention via API.")
        end
    elseif Sys.isapple()
        # macOS: Use caffeinate
        # -i: prevent idle sleep
        # -w: wait for the process with the given PID to exit
        run(`caffeinate -i -w $(getpid())`, wait=false)
        println(">> macOS: System sleep prevention enabled (caffeinate).")
    elseif Sys.islinux()
        # Linux: Use systemd-inhibit
        # It creates a lock that prevents sleep/idle.
        try
            run(`systemd-inhibit --why="Numerical Experiments" --who="Julia" --what=idle:sleep sleep infinity`, wait=false)
            println(">> Linux: System sleep prevention enabled (systemd-inhibit).")
        catch
            println(">> Linux: systemd-inhibit not found. System may still sleep.")
        end
    end
end

enable_stay_awake()

# ================== set parameters ==================
mode_idx = findfirst(x -> x == "--mode", ARGS)
is_full = mode_idx !== nothing && mode_idx < length(ARGS) && ARGS[mode_idx+1] == "full"

if is_full
    @info "RUNNING IN FULL MODE: This will reproduce the full tables."
else
    @warn "RUNNING IN TEST MODE (Default): Only a subset of instances will be reproduced. To reproduce full results, add '--mode full' at the end of the command."
end

if is_full
    Mlist = [1,2] 
    Tlist = [5,10]
    stronglist = [0,1]  # weak and strong satiation 
    testn = 5
    Nlist = [30,50,70]
    time_bound = 60*60
else
    Mlist = [1] 
    Tlist = [5] 
    stronglist = [0,1]   
    testn = 5
    Nlist = [30]
    time_bound = 60*60 
end

# =======================================================

using JuMP,Pkg,CSV,DataFrames,Gurobi,BenchmarkTools,MosekTools, Dates, JSON , XLSX 
import LinearAlgebra,Random
Random.seed!(314)
import MathOptInterface as MOI  

# --- PATH SETUP ---
const SCRIPT_DIR = @__DIR__
const PROJECT_ROOT = dirname(SCRIPT_DIR) 
# ------------------

include(joinpath(PROJECT_ROOT, "src", "Gene_data.jl")) 
include(joinpath(PROJECT_ROOT, "src", "Multilinear.jl")) 
include(joinpath(PROJECT_ROOT, "src", "Conic_Env.jl"))
include(joinpath(PROJECT_ROOT, "src", "SCIP.jl"))


# -------------- Set parameters --------------
println("=" ^ 50)
println("Set optimization parameters ...") 

name = "CompareModel"
posi = 0 
mix = 0 
posi_ratio = 0   
# opt parameters 
v0 = 1  
stop_gap = 0.005
init = 0
init_X = Any[]
output = 1
IP = 1  # change IP to 0 to obtain continuous relaxation 

println("Numerical study: ", name)
println("Number of products N = ", Nlist)
println("Memory length M = ", Mlist)
println("Planning horizon T = ", Tlist ) 
println("Time limit = ", time_bound, " seconds ")
println("Stop optimality gap = ", stop_gap*100, "%")

println("Start solving HAP problems...")
println("=" ^ 50)
 
Result = [ ]
LP_result = [ ] 
results_folder = joinpath(PROJECT_ROOT, "results", "Table3_4_Sec6_1_CompareModel")
logs_folder = joinpath(results_folder, "logs")
mkpath(logs_folder)  
time_bound_opt = 2*60*60

GC.gc() 
for M in Mlist 
    if M ==1 
        rank_set = [1]
    else 
        rank_set = [ [1,2], [2,1]]
    end 
    for strong in stronglist   
        json_path = joinpath(PROJECT_ROOT, "data", "Satiation_M"*string(M)*"_strong"*string(strong)*".json")
        test_data = open(json_path, "r") do io
        JSON.parse(read(io, String))
        end 
        for tn in 1:length(Nlist)
            N = Nlist[tn] 
            s = zeros(N,M) # only consider satiation cases 

            for k in 1:testn 
                print(M, ",", N , ",", strong,   "\n")

                r = test_data[(tn-1)*testn + k ][1]
                u = test_data[(tn-1)*testn + k ][2]
                beta0 = test_data[(tn-1)*testn + k ][3] 
                print(length(r), ",", length(u) , ",", length(beta0[1]), "\n")

                beta = zeros(N,M)
                # set up parameters 
                for i in 1:N 
                    for m in 1:M 
                        beta[i,m] = beta0[m][i]
                    end 
                end  
            
                betanew = zeros(N,M,M+1)
                for t in 2:M+1
                    for i in 1:N  
                        if (t<=M) 
                            for m in 1:t-1 
                                betanew[i,m,t] = beta[i,m]
                            end 
                        else 
                            for m in 1:M 
                                betanew[i,m,t] =  beta[i,m] 
                            end 
                        end 
                    end 
                end  
                
                # optimization for different periods 
                for T in Tlist 
                    # unconstrained 
                    cardinal = N 
                    showtime = T 

                    # four models 
                    print("----------solve conic model \n")  
                    logname1 = joinpath(logs_folder,"conic_M"*string(M)*"_N"*string(N)*"_T"*string(T)*"_inst"*string(k)*"_strong"*string(strong)*".txt" )
                    conic = MIP(T,N,M,r,beta,u,v0,IP, time_bound,stop_gap,init,init_X, output, logname1, s, betanew,rank_set, cardinal, showtime , 1 )  
                    
                    print("-----------solve envelope tight \n")
                    logname2 = joinpath(logs_folder,"env_M"*string(M)*"_N"*string(N)*"_T"*string(T)*"_inst"*string(k)*"_strong"*string(strong)*".txt")
                    if M == 2 
                        env = MIP_M2_Satiation_Strong(T, N, M, r, beta, u, v0, time_bound, IP, cardinal, showtime, stop_gap,
                    init, init_X, logname2, s, output)  
                    else
                        env = MIP_M1_Strong( T,N,M,r,beta,u,v0,IP, time_bound,stop_gap,init,init_X, output, logname2, s,  cardinal, showtime ) 
                    end  
                    
                    print("------------solve multilinear \n")
                    logname3 = joinpath(logs_folder,"ml_M"*string(M)*"_N"*string(N)*"_T"*string(T)*"_inst"*string(k)*"_strong"*string(strong)*".txt")
                    ml = Multilinear(T,N,M, r,u,beta,v0,IP,stop_gap,time_bound, output,logname3, s, cardinal,showtime )
                    
                    print("------------solve scip \n")
                    scipres = FracBinary(T,N,M,r,u,beta,v0,IP,stop_gap,time_bound, cardinal,showtime)

                    push!(Result,[M, strong, N, k , T, cardinal, showtime, conic[1], conic[2], conic[3],  
                    env[1], env[3], env[4],
                    ml[1], ml[2], ml[3], 
                    scipres[1], scipres[2], scipres[3]    ])

                    # lp relaxation
                    coniclp = MIP(T,N,M,r,beta,u,v0,0, time_bound,stop_gap,init,init_X, output, "", s, betanew,rank_set, cardinal, showtime , 1 )
                    if M == 2 
                        envlp = MIP_M2_Satiation_Strong(T, N, M, r, beta, u, v0, time_bound,0, cardinal, showtime, stop_gap,init, init_X,logname2, s, output)  
                    else
                        envlp = MIP_M1_Strong( T,N,M,r,beta,u,v0,0, time_bound,stop_gap,init,init_X, output, logname2, s,  cardinal, showtime   ) 
                    end 
                    ml_lp = Multilinear(T,N,M, r,u,beta,v0,0,stop_gap,time_bound, output,"", s, cardinal,showtime )

                    # find the optimal IP solution 
                    bmk = 0 
                    if env[4] <= stop_gap+0.0001  # already optimal 
                        bmk = env[1]
                    else 
                        if M == 2 
                            envopt = MIP_M2_Satiation_Strong(T, N, M, r, beta, u, v0, time_bound_opt, 1, cardinal, showtime, stop_gap, 1, env[2], "", s, output)  
                        else
                            envopt = MIP_M1_Strong( T,N,M,r,beta,u,v0,IP, time_bound_opt,stop_gap,1,env[2], output, "", s,  cardinal, showtime ) 
                        end   
                        bmk = envopt[1]
                    end 
                    push!(LP_result,[M, strong, N, k , T, cardinal, showtime, (coniclp[1] - bmk)/bmk, 
                    (envlp[1] - bmk)/bmk,  (ml_lp[1] - bmk)/bmk  ])
                    GC.gc()
                end 
            end 
        end
    end 
end 



temp_matrix  = tupe2mat(Result) 
f1 = joinpath(results_folder,"result_detail.csv")
temp_data = DataFrame(temp_matrix, :auto) 
rename!(temp_data, ["M", "beta", "N", "k","T","cardsize", "showtime","conic_obj","conic_time","conic_gap", "env_obj", "env_time", "env_gap", "ml_obj", "ml_time","ml_gap", "scip_obj","scip_time","scip_gap"]) 
CSV.write(f1,  temp_data, writeheader=true)

temp_matrix2  = tupe2mat(LP_result) 
f2 = joinpath(results_folder,"relaxation_detail.csv")
relax_data = DataFrame(temp_matrix2, :auto) 
rename!(relax_data, ["M", "beta", "N", "k","T","cardsize", "showtime","conic_gap","env_gap","ml_gap"]) 
CSV.write(f2,  relax_data, writeheader=true)

# ========================== compute table 3 ========================= 
println("Compare models ...")
println("=" ^ 50)

# Get current time
t = Dates.now() 
timestr = Dates.format(t, "HHMM")

function compute_res(testn, time_bound, soludata)
    solved = 0 
    total_time = 0 
    total_endgap = 0 
    rename!(soludata, ["obj","runtime","gap"])
    for (index, row) in enumerate(eachrow(soludata))
        if row.gap <= stop_gap 
            solved +=1 
            total_time += row.runtime 
        else 
            total_endgap += row.gap
        end
    end 
    if solved ==0 
        return [solved, time_bound, round(100*total_endgap/testn, digits=4)] 
    elseif solved == testn 
        return [solved, round(total_time/testn, digits=4), round(100*total_endgap/testn ) ] 
    else
        return [solved, round(total_time/solved, digits=4), round(100*total_endgap/(testn-solved), digits=4)] 
    end
end 

summary = [ ] 
for T in Tlist
    for M in Mlist   
        for N in Nlist 
            for stv in stronglist
                filterdata = filter(row -> row.M == M && row.N == N && row.T ==T && row.beta==stv, temp_data)
                res = zeros(16)
                res[1:4] = [T,M,N,stv] 
                # compute result for each model. 
                res1 = compute_res(testn, time_bound, select(filterdata, :conic_obj, :conic_time, :conic_gap))
                res[5:7] = res1 
                res2 = compute_res(testn, time_bound, select(filterdata, :env_obj, :env_time, :env_gap)) 
                res[8:10] = res2 
                res3 = compute_res(testn, time_bound, select(filterdata, :ml_obj, :ml_time, :ml_gap)) 
                res[11:13] = res3 
                res4 = compute_res(testn, time_bound, select(filterdata, :scip_obj, :scip_time, :scip_gap)) 
                res[14:16] = res4 
                push!(summary, res)
            end 
        end 
    end 
end 

summary = tupe2mat(summary) 
summary = DataFrame(summary, :auto)   
rename!(summary, ["T", "M", "N", "beta", "conic_solved", "conic_time", "conic_gap", "env_solved", "env_time", "env_gap", "ml_solved", "ml_time", "ml_gap", "scip_solved", "scip_time", "scip_gap"  ])
summary.beta = replace(summary.beta, 0 => "W", 1 => "S")
table3 = joinpath(results_folder,"Table3_IPsolu_"*timestr*".xlsx")
XLSX.writetable(table3, summary)



# ================== compute table 4 ================= 

sumdata = [ ] 
for T in Tlist
    for M in Mlist   
        for N in Nlist 
            for stv in stronglist
                filterdata = filter(row -> row.M == M && row.N == N && row.T ==T && row.beta==stv, relax_data)
                push!(sumdata, [T,M,N,stv, round(100*mean(filterdata.conic_gap),digits=4), round(100*mean(filterdata.env_gap),digits=4), round(100*mean(filterdata.ml_gap),digits=4) ])
            end 
        end 
    end 
end 

sumdata = tupe2mat(sumdata) 
sumdata = DataFrame(sumdata, :auto)   
rename!(sumdata, ["T", "M", "N", "beta", "conic_root_gap","env_root_gap", "ml_root_gap" ])
sumdata.beta = replace(sumdata.beta, 0 => "W", 1 => "S")


table4 =joinpath(results_folder, "Table4_LPsolu_"*timestr*".xlsx"   )
XLSX.writetable(table4, sumdata)

println("Finish numerical study ", name, "!")
 