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

using JuMP,Gurobi,BenchmarkTools, LinearAlgebra  
import MathOptInterface as MOI   
using SCIP 
 
# ==================== Optimization Model ====================================
# This module solves HAP problem via the globle solver SCIP. 
# ================================================================================

 
function FracBinary(T,N,M,r,u,beta,v0,IP,stop_gap,time_bound, cardinal,showtime)
    # Create model with SCIP
    model = Model(SCIP.Optimizer)

    # Variables  
    @variable(model, 0<= x[1:N,1:T]<=1)  
    @variable(model, w[1:T]) # revenue in each period 
    @variable(model, a[1:N, 1:T]) # attraction value 
    # objective
    @objective(model, Max, sum(w[t] for t in  1:T)  )

    # Constraints
    for t in 1:T
        if t == 1
            @NLconstraint(model, [i in 1:N], a[i,t] == exp(u[i])  )  
        elseif t==2
            @NLconstraint(model, [i in 1:N], a[i,t] == exp(u[i]+beta[i,1]*x[i,t-1])  )
        elseif (M==1)
            @NLconstraint(model, [i in 1:N], a[i,t] == exp(u[i]+beta[i,1]*x[i,t-1])  )
        else
            @NLconstraint(model, [i in 1:N], a[i,t] == exp(u[i]+beta[i,1]*x[i,t-1] + beta[i,2]*x[i,t-2])  )
        end 
    end 

    @NLconstraint(model, [t in 1:T], w[t] <= sum(r[i]*x[i,t]*a[i,t] for i in 1:N)/(v0+ sum(x[i,t]*a[i,t] for i in 1:N) ))

    @constraint(model, [i in 1:N], sum(x[i,t] for t in 1:T)<= showtime ) 
    @constraint(model, [t in 1:T], sum(x[i,t] for i in 1:N)<= cardinal ) 

    
    if IP == 1 
        set_binary.(x)
        set_optimizer_attribute(model, "limits/time", time_bound) 
        set_optimizer_attribute(model, "limits/gap", stop_gap  )   
    end 
    # Solve
    optimize!(model)

    # Results
    return objective_value(model), solve_time(model), MOI.get(model, MOI.RelativeGap())
end 
 