"""
optimization.py — Modelli di Ottimizzazione (Gurobi)

Questo script implementa tutti i modelli di ottimizzazione del progetto, raccolti nella
classe AssortmentOptimizer. Il problema da risolvere è un'ottimizzazione stocastica a
due stadi: decidere quali prodotti inserire in assortimento (a) e quante unità ordinare (y)
in modo da massimizzare il profitto atteso, tenendo conto dell'incertezza sulla domanda.

Struttura della classe AssortmentOptimizer:

  Inizializzazione e strutture interne:
    __init__()              : riceve tutti i parametri economici e fisici del problema.
    _prepare_mappings()     : inizializza strutture sparse utili nei modelli di ottimizzazione.
    _mean_parameters()      : precalcola coefficienti pesati per la funzione obiettivo stocastica.

  Costruzione modelli Gurobi:
    build_model_SAA()           : modello stocastico completo a due stadi.
    build_model_NN()            : surrogata NN: la rete neurale stima E[Q(a,y)].
    build_model_hybrid_det()    : ibrido deterministico: Ridge predice la domanda media.
    build_model_hybrid_NN_stoc(): ibrido stocastico: NN predice la domanda per ogni scenario.
    build_recurse_only()        : solo il modello di ricorso (secondo stadio), per valutare soluzioni.

  Soluzione e metodi di supporto:
    solve()                         : lancia Gurobi e stampa il risultato.
    solve_recurse()                 : fissa (a,y) e risolve il solo ricorso.
    get_single_scenario_solution()  : risolve SAA su un singolo scenario (usato per generare il dataset NN).
    debug_valore_assortimento()     : valuta il profitto totale di un assortimento dato.
    close_model()                   : libera la memoria del modello Gurobi.

  Export risultati:
    export_to_text_SAA()           : report per il modello SAA.
    export_to_text_NN()            : report per il modello NN surrogata.
    export_to_text_hybrid_det()    : report per il modello ibrido deterministico.
    export_to_text_hybrid_NN_stoc(): report per il modello ibrido stocastico.


Nota sui prodotti:

  I prodotti sono organizzati in triplette: ogni gruppo di 3 indici consecutivi
  (3*j, 3*j+1, 3*j+2) rappresenta le 3 opzioni di scaffale (Alto/Occhi/Basso) per
  lo stesso prodotto unico. Ci sono dei vincoli SOS1 nei modelli che impediscono 
  che vengano scelti Alto e Basso contemporaneamente (ma è ammesso Occhi + uno degli altri).

"""

import numpy as np
import pandas as pd
from scipy.sparse import coo_array
import gurobipy as gp
from gurobipy import GRB
from gurobi_ml import add_predictor_constr
import joblib



class AssortmentOptimizer:
    def __init__(self,prodotti_df,sub_matrix,prices,costs,holding_costs,capacity,fridge_capacity,fixed_cost,penalty_lost_sales,penalty_service_level,service_level_target,flag_shelf=False):
        self.df=prodotti_df
        self.sub_matrix=sub_matrix
        self.prices=prices
        self.costs=costs
        self.h_costs=holding_costs
        self.cap_shelf=capacity
        self.cap_fridge=fridge_capacity
        self.f_cost=fixed_cost
        self.pen_lost=penalty_lost_sales
        self.pen_SL=penalty_service_level
        self.SL_target=service_level_target #percentuale di livello di servizio da mantenere
        self.shelf=flag_shelf #Se True, i vincoli di capacità vengono applicati a livello di tipologia di scaffale; 
                            #se False, vengono aggregati a livello globale.
        #dizionario per memorizzare le istanze dei diversi modelli di ottimizzazione
        self.models={"saa":None,
                     "recurse":None,
                     "hybrid_det":None,
                     "hybrid_nn_stoc":None,
                     "nn":None}
        self.results={}
    
    def _prepare_mappings(self):
        """metodo per inizializzare delle matrici e delle strutture utilizzate nei modelli di ottimizzazione"""
        #evita ricalcolo se già chiamato
        if hasattr(self,'_mappings_ready'):
            return
        self._mappings_ready=True
        self.num_opzioni=len(self.df)
        self.I=range(self.num_opzioni)
        self.num_unique_prod=self.num_opzioni//3 #struttura del dataset con prodotti unici su tre scaffali diversi

        #Logica scaffali e prodotti (Gnocchi)
        self.opzioni_per_scaffale={} 
        self.gnocchi_id=[]
        for idx, row in self.df.iterrows():
            s_key=row['Scaffale']
            if s_key not in self.opzioni_per_scaffale:
                self.opzioni_per_scaffale[s_key]=[]
            self.opzioni_per_scaffale[s_key].append(idx)
            if "Gnocchi" in str(row['Prodotto']):
                self.gnocchi_id.append(idx)
        
        self.capacity_shelf_tot=sum(self.cap_shelf.values())
        self.capacity_frdge_tot=sum(self.cap_fridge.values())

        #Matricidi aggregazione e sostituzione
        #identifica le coordinate (riga, colonna) dove esiste una probabilità di sostituzione > 0, escludendo l'ultima riga/colonna (no-buy)
        self.rows_sost,self.cols_sost=np.where(self.sub_matrix[:-1,:-1]>0)
        self.num_couples=len(self.rows_sost)
        #memorizzazione come coo_array per efficienza di memoria
        #matrix_in: mappa chi riceve la domanda sostituita, matrix_out: da chi esce la domamnda
        self.matrix_in=coo_array((np.ones(self.num_couples),(self.cols_sost,np.arange(self.num_couples))),shape=(self.num_opzioni,self.num_couples))
        self.matrix_out=coo_array((np.ones(self.num_couples),(self.rows_sost,np.arange(self.num_couples))),shape=(self.num_opzioni,self.num_couples))

        col_indices=np.arange(self.num_opzioni)
        row_indices=np.repeat(np.arange(self.num_unique_prod),3) #crea array [0,0,0, 1,1,1, 2,2,2...] per collegare 3 opzioni a 1 prodotto
        self.matrix_sum_prod=coo_array((np.ones(self.num_opzioni),(row_indices,col_indices)),shape=(self.num_unique_prod,self.num_opzioni))

    def _mean_parameters(self,demand,prob_scen):
        """
        calcola valori medi per i coefficienti nel problema stocastico SAA.

        Moltiplica i vettori dei parametri economici per le probabilità degli scenari, generando matrici pronte per essere 
        inserite direttamente nella funzione obiettivo del modello SAA.
        
        Input:
            1. matrice degli scenari di domanda di dimensione (n_scenari, n_opzioni + 1);
            2. vettore probabilità scenario.

        Dove viene usato: nel metodo build_model_SAA e debug_valore_assortimento; e negli script prepare_data.py e active_learning.py per il dataset del surrogato
        
        """
        self._prepare_mappings() 
        self.p=np.array(prob_scen)
        demand=np.array(demand)
        self.S=range(demand.shape[0])

        #Calcolo target ServiceLevel Fisso (basato sulla media degli scenari)
        self.d=demand[:,:-1].T
        domanda_tot_prod=self.d.sum(axis=0)
        self.domanda_tot_prod_SL=domanda_tot_prod*(1-self.SL_target)

        #Calcolo parametri mediati 
        self.price_medio=np.outer(self.prices,self.p)
        self.pen_lost_medio=np.outer(self.pen_lost,self.p)
        self.pen_inventory_medio=np.outer(self.h_costs,self.p)
        self.pen_sost_medio=np.outer((0.2*np.ones(self.pen_lost.size)),self.p)

        self.pen_SL_medio=self.pen_SL*self.p    
        

    def build_model_SAA(self,demand,probability,gap=False,write=False):
        """costruisce il modello di ottimizzazione SAA"""

        self.models["saa"]=gp.Model("SAA_Esteso")
        model=self.models["saa"]

        self._mean_parameters(demand,probability)
        demand_T=self.d
        num_p=self.num_opzioni
        num_s=demand_T.shape[1]
        num_unique_prod=self.num_unique_prod

        model.ModelSense=GRB.MAXIMIZE

        model.setParam('Method',4)
        model.setParam('Heuristics',0.3)
        model.setParam('MIPFocus',2)
        model.setParam('Presolve',2)
        model.setParam('Symmetry',2)
        if gap:
            model.setParam('MIPGap',0.1)
        else:
            model.setParam('MIPGap',0.05)
        model.setParam('Cuts',2)
        if len(probability)>10:
            model.setParam('NoRelHeurTime',300)

        #Varibili primo stadio
        a=model.addMVar((num_p,1),lb=0.0,vtype=GRB.BINARY,name="a")#Variabile binaria di appartenenza allo scaffale e all'altezza stabilita
        y=model.addMVar((num_p,1),lb=0.0,vtype=GRB.INTEGER,name="y")#quantità di prodotto (unità=1kg)
        a.BranchPriority=100
        y.BranchPriority=50


        #Variabili secondo stadio
        x=model.addMVar((num_p,num_s),lb=0.0,ub=demand_T,vtype=GRB.CONTINUOUS,name="x")#quantità primaria
        sost=model.addMVar((self.num_couples,num_s),lb=0.0,vtype=GRB.CONTINUOUS,name="sost")#sostituita
        r=model.addMVar((num_p,num_s),lb=0.0,ub=300.0,vtype=GRB.CONTINUOUS,name="r")#rimanenza
        lost=model.addMVar((num_p,num_s),lb=0.0,ub=demand_T,vtype=GRB.CONTINUOUS,name="lost_demand")#domanda persa
        service_level=model.addMVar(num_s,lb=0.0,ub=demand_T.sum(axis=0).reshape(-1),name="service_level")#domanda legata al livello diservizio sul prodotto totale

        #sostituzioni
        sost_out=self.matrix_out@sost
        sost_in=self.matrix_in@sost
        lost_tot=lost.sum(axis=0)

        ricavo_atteso=(self.price_medio*x).sum()+(self.price_medio*sost_in).sum()
        penali_attese=(self.pen_lost_medio*lost).sum()+(self.pen_sost_medio*sost_out).sum()+(self.pen_inventory_medio*r).sum()+(self.pen_SL_medio*service_level).sum()

        secondo_stadio_atteso=ricavo_atteso-penali_attese

        f_vector=np.array([self.f_cost[self.df.iloc[i]['Scaffale']]for i in self.I])
        costo_primo_stadio=(self.costs@y)+(f_vector@a)

        #f obj
        model.setObjective(secondo_stadio_atteso-costo_primo_stadio)
        
        #Mutua esclusione prodotti
        for i in range(num_unique_prod):
            model.addSOS(GRB.SOS_TYPE1,[a[3*i],a[3*i+1]])

        
        if self.shelf:
            for s_name,indici in self.opzioni_per_scaffale.items():

                indici_freschi=[i for i in indici if i in self.gnocchi_id]
                indici_secchi=[i for i in indici if i not in self.gnocchi_id]        
            
                #Capacità frigo
                if indici_freschi:
                    cap_f=self.cap_fridge[s_name]
                    model.addConstr(y[indici_freschi].sum()<=cap_f,name=f"Capacita_fresco_{s_name}")
                    
                #Capacità stagera
                if indici_secchi:
                    cap_s=self.cap_shelf[s_name]
                    model.addConstr(y[indici_secchi].sum()<=cap_s,name=f"Capacita_secco_{s_name}")
        else:
            tutti_freschi=self.gnocchi_id
            tutti_secchi= [i for i in range(self.num_opzioni) if i not in self.gnocchi_id]
            
            if tutti_freschi:
                model.addConstr(y[tutti_freschi].sum()<=self.capacity_frdge_tot,name="Capacita_fresco_globale")
            if tutti_secchi:
                model.addConstr(y[tutti_secchi].sum()<=self.capacity_shelf_tot,name="Capacita_secco_globale")


        model.addConstr(y<=200*a,name=f"Link_ay") ##200 è LA GRANDE M

        #Vincoli del ricorso
        
        model.addConstr(x<=demand_T*a,name="Primary_sell")

        model.addConstr(sost<=100*a[self.cols_sost,0].reshape(-1,1),name="Link_asost")

        M_c_saa=self.sub_matrix[self.rows_sost,self.cols_sost]
        model.addConstr(sost+(M_c_saa[:,np.newaxis]*x[self.rows_sost,:])<=(M_c_saa[:,np.newaxis]*demand_T[self.rows_sost,:]),name="Subst_constrs")

        #per ogni prodotto che riceve sostituzione, questo può effettivamente riceverla SE è in assortimento
        #in sost io ho [indice di coppia,scenario]. Io posso dire sost[k dove k sono tutti quegli indici in cui si arriva in j]<=big_M*a[j,0]

        model.addConstr(x+sost_in+r==y,name="Stock_balance")

        model.addConstr(x+lost+sost_out==demand_T,name="Demand_balance_all")

        model.addConstr(lost_tot-service_level<=self.domanda_tot_prod_SL,name="ServiceLevel_goal")

        model.update()

        if write:
            model.write("il_mio_modello.lp")

        model._a_saa=a
        model._y_saa=y
        model._x_saa=x
        model._sost_saa=sost
        model._r_saa=r
        model._lost_saa=lost
        model._service_level_saa=service_level

    def build_model_NN(self, mlp_path='surrogate_mlp.pkl', scaler_path_input='scaler_input.pkl',scaler_path_target='scaler_target.pkl', write=False):
        """costruisce il modello di ottimizzazione con ricorso appreso dalla rete neurale"""

        self.models["nn"]=gp.Model("NN_model")
        model=self.models["nn"]
        self._prepare_mappings()#Solo struttura
        num_p=self.num_opzioni

        model.ModelSense=GRB.MAXIMIZE

        #Variabili di Primo Stadio
        a=model.addMVar((num_p,1),vtype=GRB.BINARY,name="a")
        y=model.addMVar((num_p,1),vtype=GRB.INTEGER,name="y")
        
        scaler_X=joblib.load(scaler_path_input)
        X_min=scaler_X.data_min_
        X_range=scaler_X.data_range_

        nn_input_scaled=model.addMVar(num_p,lb=0.0,ub=1.0,name="nn_input_scaled")
        for i in range(num_p):
            if X_range[i] > 0: 
                model.addConstr(nn_input_scaled[i]==(y[i]-X_min[i])/X_range[i], name=f"scale_input_{i}")
            else:
                model.addConstr(nn_input_scaled[i]==0.0, name=f"scale_input_zero_{i}")

        
        scaler_y=joblib.load(scaler_path_target)
        mu=scaler_y.mean_[0]#La media calcolata dallo scaler
        sigma=scaler_y.scale_[0]#La deviazione standard calcolata dallo scaler

        #Rappresenta E[Q(a,y,xi)] predetto dalla rete
        expected_recourse_aus=model.addVar(lb=-GRB.INFINITY,name="recourse_value_standardized")
        mlp_model = joblib.load(mlp_path)

        #La NN predice il profitto del secondo stadio dati'a'e'y'
        add_predictor_constr(model,mlp_model,nn_input_scaled,expected_recourse_aus)

        expected_recourse=model.addVar(lb=-GRB.INFINITY,name="expected_recourse")


        model.addConstr(expected_recourse==expected_recourse_aus*sigma+mu,name="denormalize_output")

        f_vector=np.array([self.f_cost[self.df.iloc[i]['Scaffale']] for i in self.I])
        costo_primo_stadio=(self.costs@y)+(f_vector@a)

        #f obj
        model.setObjective(expected_recourse-costo_primo_stadio)

        #5.Vincoli(Solo quelli fisici di primo stadio)
        model.addConstr(y<=200*a)#Link stock-assortimento
        model.addConstr(y>=a, name="link_a_y_min")
       
        #Mutua esclusione prodotti
        for i in range(self.num_unique_prod):
            model.addSOS(GRB.SOS_TYPE1,[a[3*i], a[3*i+1]])

        
        if self.shelf:
            for s_name, indici in self.opzioni_per_scaffale.items():

                indici_freschi = [i for i in indici if i in self.gnocchi_id]
                indici_secchi  = [i for i in indici if i not in self.gnocchi_id]        
            
                #Capacità frigo
                if indici_freschi:
                    cap_f=self.cap_fridge[s_name]
                    model.addConstr(y[indici_freschi].sum()<=cap_f,name=f"Capacita_fresco_{s_name}")
                    
                #Capacità stagera
                if indici_secchi:
                    cap_s=self.cap_shelf[s_name]
                    model.addConstr(y[indici_secchi].sum()<=cap_s,name=f"Capacita_secco_{s_name}")
        else:
            tutti_freschi=self.gnocchi_id
            tutti_secchi=[i for i in range(self.num_opzioni) if i not in self.gnocchi_id]
            
            if tutti_freschi:
                model.addConstr(y[tutti_freschi].sum()<=self.capacity_frdge_tot,name="Capacita_fresco_globale")
            if tutti_secchi:
                model.addConstr(y[tutti_secchi].sum()<=self.capacity_shelf_tot,name="Capacita_secco_globale")

                model.update()

        if write:
            model.write("il_mio_modello.lp")

        model._a_nn = a
        model._y_nn = y

    #L'idea è generare solo domanda attesa dato l'assortimento, poi Gurobi ottimizza avendo una 'domanda deterministica'
    #NB il predittore impara una domanda più complessa: infatti capisce che domanda si ottiene variando l'assortimento iniziale, 
    #quindi riesce a valutare una logica di sostituzione più completa. 
    #Nella fase di ottimizzazione, vediamo invece come, data questa domanda più 'vera', si comporta la sostituzione in-store.
    #Svantaggio: passiamo un'informazione aggragata media, si perde l'info sulla variabilità, quindi l'assortimento è più fittato sulla domanda, poche sostituzioni
    def build_model_hybrid_det(self,demand_tot_mean, mlp_path='best_linear_models_list.pkl', scaler_path_target='scaler_y.pkl',write=False):
        """costruisce il modello di ottimizzazione ibrido deterministico"""

        self.models["hybrid_det"]=gp.Model("Hybrid_det")
        model=self.models["hybrid_det"]
        self._prepare_mappings()#struttura
        num_p=len(self.df)
        num_unique_prod=self.num_unique_prod

        model.ModelSense=GRB.MAXIMIZE

        model.setParam('MIPFocus',2)
        model.setParam("Presolve",2)
        model.setParam('NoRelHeurTime',60)
        model.setParam('Cuts',3)
        #Variabili
        a=model.addMVar(num_p,vtype=GRB.BINARY,name="a")
        y=model.addMVar(num_p,lb=0.0,vtype=GRB.INTEGER,name="y")
        a.BranchPriority=100
        y.BranchPriority=50

        scaler_y=joblib.load(scaler_path_target)
        y_min=scaler_y.data_min_
        y_range=scaler_y.data_range_

        #d_pred è l'output del predittore:rappresenta la DOMANDA TOTALE(considerando primaria+sost a catalogo)
        models_list = joblib.load(mlp_path)
        d_pred_raw_list = []
        
        for i in range(num_p):
            v = model.addVar(lb=0.0, vtype=GRB.CONTINUOUS, name=f"d_pred_ml_{i}")         
            add_predictor_constr(model, models_list[i], a, [v])    
            d_pred_raw_list.append(v)
        delta_no_buy_raw=model.addVar(lb=0.0, vtype=GRB.CONTINUOUS,name=f"delta_no_buy_raw")
        add_predictor_constr(model,models_list[num_p],a,[delta_no_buy_raw])

        #Riscaliamo l'output
        d_pred = model.addMVar(num_p, lb=0.0, ub=200 ,vtype=GRB.CONTINUOUS, name="d_pred")
        for i in range(num_p):
            model.addConstr(d_pred[i]==(d_pred_raw_list[i] * y_range[i]) + y_min[i], name=f"scale_output_{i}")
        delta_no_buy=model.addVar(lb=0.0, vtype=GRB.CONTINUOUS, name=f"delta_no_buy")
        model.addConstr(delta_no_buy==(delta_no_buy_raw*y_range[num_p])+y_min[num_p], name=f"scale_nobuy")

        d_pred_real=model.addMVar(num_p, lb=0.0, ub=200, name="real_demand")
        model.addConstr(d_pred_real<=d_pred)
        model.addConstr(d_pred_real<=200*a, name="Liank_a_to_d")
        model.addConstr(d_pred_real>=d_pred-200*(1-a), name="Link_a_to_d_min")

        model.addConstr(y>=a)

        #Variabili per coerenza con SAA
        x_prim=model.addMVar(num_p,lb=0.0,name="x_primaria")
        sost=model.addMVar(self.num_couples,lb=0.0,name="sost")
        r=model.addMVar(num_p,lb=0.0,name="r")
        lost=model.addMVar(num_p,lb=0.0,name="lost_demand")
        service_level=model.addVar(lb=0.0,name="service_level")

        #Precalcolati utili
        sost_in=self.matrix_in@sost
        sost_out=self.matrix_out@sost
        pen_lost_mean=self.pen_lost.mean()
        lost_tot=lost.sum()

        tolleranza_lost=demand_tot_mean*(1-self.SL_target)

        #f obj
        ricavo=(self.prices@x_prim)+(self.prices@sost_in)

        #Penali e Costi
        penali=(self.pen_lost@lost)+(self.h_costs@r)+(pen_lost_mean*delta_no_buy)+(self.pen_SL*service_level)+(0.2*sost_out.sum())
        f_vector=np.array([self.f_cost[self.df.iloc[i]['Scaffale']] for i in self.I])
        costi_fissi_e_acquisto=(self.costs@y)+(f_vector@a)

        model.setObjective(ricavo-penali-costi_fissi_e_acquisto)
 
        #Mutua esclusione
        for i in range(num_unique_prod):
            model.addSOS(GRB.SOS_TYPE1, [a[3*i], a[3*i+1]])

        if self.shelf:
            for s_name, indici in self.opzioni_per_scaffale.items():

                indici_freschi = [i for i in indici if i in self.gnocchi_id]
                indici_secchi  = [i for i in indici if i not in self.gnocchi_id]        
            
                #Capacità frigo
                if indici_freschi:
                    cap_f=self.cap_fridge[s_name]
                    model.addConstr(y[indici_freschi].sum()<=cap_f,name=f"Capacita_fresco_{s_name}")
                    
                #Capacità stagera
                if indici_secchi:
                    cap_s=self.cap_shelf[s_name]
                    model.addConstr(y[indici_secchi].sum()<=cap_s,name=f"Capacita_secco_{s_name}")
        else:
            tutti_freschi=self.gnocchi_id
            tutti_secchi=[i for i in range(self.num_opzioni) if i not in self.gnocchi_id]
            
            if tutti_freschi:
                model.addConstr(y[tutti_freschi].sum()<=self.capacity_frdge_tot,name="Capacita_fresco_globale")
            if tutti_secchi:
                model.addConstr(y[tutti_secchi].sum()<=self.capacity_shelf_tot,name="Capacita_secco_globale")


        model.addConstr(y<=200*a,name=f"Link_ay") ##200 è LA GRANDE M
        
        model.addConstr(x_prim<=200*a,name="Primary_sell")
        
        model.addConstr(sost<=15*a[self.cols_sost],name=f"Link_asost") ##100 è LA GRANDE M

        M_c_saa=self.sub_matrix[self.rows_sost,self.cols_sost]
        model.addConstr(sost+(M_c_saa*x_prim[self.rows_sost])-(M_c_saa*d_pred_real[self.rows_sost])<=0,name="Subst_constrs")

        #per ogni prodotto che riceve sostituzione, questo può effettivamente riceverla SE è in assortimento
        #in sost io ho [indice di coppia, scenario]. Io posso dire sost[k dove k sono tutti quegli indici in cui si arriva in j]<= big_M * a[j,0]  

        model.addConstr(x_prim+sost_in+r==y,name="Stock_balance")

        model.addConstr(x_prim+lost+sost_out-d_pred_real==0,name="Demand_balance_all")

        #Le vendite perse, meno lo slack della penale, devono essere sotto la tolleranza
        model.addConstr(lost_tot-service_level<=tolleranza_lost,name="ServiceLevel_goal")

        model.update()
        if write:
            model.write("modello_hybrid_confronto.lp")

        model._a_hnn=a
        model._y_hnn=y
        model._x_hnn=x_prim
        model._sost_hnn=sost
        model._r_hnn=r
        model._lost_hnn=lost
        model._service_level_hnn=service_level
        model._d_pred_hnn=d_pred_real
        model._demand_saa=demand_tot_mean

    def solve(self,model):
        """Lancia l'ottimizzazione e salva i risultati internamente"""
        print("Optimization starts:")
        model.optimize()
        if model.status == GRB.OPTIMAL:
            print(f"Optimization done. Expected profit: {model.ObjVal:.2f}")
            return True
        else:
            print(f"Optimization status: {model.status}")
            return True
        


    def get_single_scenario_solution(self, single_demand):
        """
        funzione che risolve SAA su un singolo scenario

        Usanto in prepare_data.py per stabilire gli assortimenti più scelti a scenario singolo per il dataset del surrogato
        
        """
        
        self.build_model_SAA(demand=single_demand,probability=[1.0],gap=True)
        mini_model=self.models["saa"]
        
        mini_model.setParam('OutputFlag',0)
        mini_model.optimize()
        
        if mini_model.Status in [GRB.OPTIMAL,GRB.TIME_LIMIT]:
            a_sol=(mini_model._a_saa.X > 0.5).astype(int).flatten() 
            y_sol=mini_model._y_saa.X.flatten()
            return a_sol,y_sol
        return None,None

    def build_recurse_only(self,n_scenarios=50,write=False):
        """
        Questo metodo crea il modello del ricorso per generare il dataset del modello surrogato e per verificare il valore del ricorso per alcuni controlli effettuati con altri metodi.
        Calcola il profitto atteso del secondo stadio, fissando le decisioni del primo (a, y).
        a_fixed: vettore binario (num_p,)
        y_fixed: vettore int (num_p,)

        Usato in: prepare_data.py e active_learning.py per il dataset del surrogato, per ottenere il valore target per la rete
        """
        
        num_p=self.num_opzioni
        num_s=n_scenarios

        self.models["recurse"]=gp.Model("Recourse_Evaluation")
        m=self.models["recurse"]
        m.setParam('OutputFlag',0)
        m.setParam('Method',1)

        x=m.addMVar((num_p,num_s),lb=0.0,name="x")
        sost=m.addMVar((self.num_couples,num_s),lb=0.0,name="sost")
        r=m.addMVar((num_p,num_s),lb=0.0,name="r")
        lost=m.addMVar((num_p,num_s),lb=0.0,name="lost")
        sl=m.addMVar(num_s,lb=0.0,name="sl")

        sost_in=self.matrix_in@sost
        sost_out=self.matrix_out@sost
        lost_tot=lost.sum(axis=0)

        a_param=m.addMVar(num_p, vtype=GRB.BINARY, name="a_param")
        y_param=m.addMVar(num_p, vtype=GRB.INTEGER, name="y_param")


        m.addConstr(x+sost_in+r==y_param.reshape(-1,1),name="Stock_Balance")

        m.addConstr(x+lost+sost_out==self.d,name="Demand_Balance")

        m.addConstr(x<=self.d*a_param.reshape(-1,1),name="Primary_sell")

        m.addConstr(sost<=100*a_param[self.cols_sost].reshape(-1,1),name="Sost_Limit")

        M_c_saa=self.sub_matrix[self.rows_sost,self.cols_sost]
        m.addConstr(sost+(M_c_saa[:,np.newaxis]*x[self.rows_sost,:])<=(M_c_saa[:,np.newaxis]*self.d[self.rows_sost,:]),name="Subst_constrs")

        m.addConstr(lost_tot-sl<=self.domanda_tot_prod_SL,name="SL_Goal")
        

        ricavo=(self.price_medio*x).sum()+(self.price_medio*sost_in).sum()
        penali=(self.pen_lost_medio*lost).sum()+(self.pen_sost_medio*sost_out).sum()+(self.pen_inventory_medio*r).sum()+(self.pen_SL_medio*sl).sum()

        m.setObjective(ricavo-penali,GRB.MAXIMIZE)
        m.update()

        if write:
            m.write("recourse_model.lp")

        m._a_rec=a_param
        m._y_rec=y_param
        m._x_rec=x
        m._sost_rec=sost
        m._r_rec=r
        m._lost_rec=lost
        m._service_level_rec=sl
    
    def solve_recurse(self, model,a_fixed,y_fixed):
            
        model._a_rec.LB=a_fixed
        model._y_rec.LB=y_fixed
        model._a_rec.UB=a_fixed
        model._y_rec.UB=y_fixed

        #model.update()
        model.reset()
        model.optimize()
        if model.Status==GRB.OPTIMAL:            
            return model.ObjVal
        else:
            return -999999 #To be sure not infeasibility

    def debug_valore_assortimento(self, vettore_totale, demand, probability, write_model=False):
        """metodo per valutare il profitto totale dato un assortimento fissato."""
        self._mean_parameters(demand,probability)
        n_scen=demand.shape[0]
        num_p=len(vettore_totale)//2
        a_test=vettore_totale[:num_p]
        y_test=vettore_totale[num_p:]
        
        #Calcolo recourse
        self.build_recurse_only(n_scenarios=n_scen,write=write_model)
        model_recurse=self.models["recurse"]
        profitto_secondo_stadio=self.solve_recurse(model=model_recurse,a_fixed=a_test,y_fixed=y_test)
        
   
        f_vector=np.array([self.f_cost[self.df.iloc[i]['Scaffale']] for i in self.I])
        costi_fissi=(f_vector@a_test)
        costi_primo_stadio=(self.costs@ y_test)+costi_fissi
        
        valore_totale=profitto_secondo_stadio-costi_primo_stadio
        return valore_totale,profitto_secondo_stadio,costi_primo_stadio, costi_fissi


    def export_to_text_SAA(self,demand_d,file_mean="analisi_risultati.txt"):
        """Metodo per report su ottimizzazione SAA"""
        model=self.models["saa"]
        a_values=(model._a_saa.X > 0.5).astype(int)
        y_values=model._y_saa.X
        r_values=model._r_saa.X
        lost_values=model._lost_saa.X
        service_level_values=model._service_level_saa.X
        x_values=model._x_saa.X
        sost_values=model._sost_saa.X

        tot_lost=lost_values.mean(axis=1)
        no_buy_medio=demand_d[:,self.num_opzioni].mean()
        domanda_potenziale=demand_d.sum()/demand_d.shape[0]
        tot_y=y_values.sum()
        media_prodotti=demand_d[:,:self.num_opzioni].mean(axis=0)
        r_medie=r_values.mean(axis=1)
        service_level_medio=service_level_values.mean()

        x_medie=x_values.mean(axis=1)
        y_vals=y_values.flatten()
        a_vals=a_values.flatten()

        media_sost=sost_values.mean(axis=1)
        tot_sost_per_prodotto=np.zeros(self.num_opzioni)
        for k, j in enumerate(self.cols_sost):
            tot_sost_per_prodotto[j]+=media_sost[k]

    
        df_report=self.df.copy()
        df_report['STOCK (y)']=y_vals
        df_report['REST (r)']=r_medie
        df_report['LOST (l)']=tot_lost
        df_report['PRIMARY DEMAND (x)']=x_medie
        df_report['SOST_RECIVED']=tot_sost_per_prodotto

        indici_assortimento=np.where(a_vals>0.5)[0].tolist()

        var_df=pd.DataFrame({'y_values':y_vals,'a_values':a_vals})
        var_df.to_csv(r'test_res\assortimento.csv',header=True,index=False)
        with open(file_mean, "w") as f_r:
                f_r.write("\n--- DETAILED ANALYSIS (mean values over the scenarios) ---\n")
                cols = ['Prodotto', 'Marca', 'Prezzo €/kg', 'Scaffale', 'STOCK (y)', 'REST (r)', 'LOST (l)', 'PRIMARY DEMAND (x)', 'SOST_RECIVED']
                f_r.write(df_report[cols].to_string(index=False, justify='left', float_format="%.2f"))

                f_r.write("\n\n" + "="*60 + "\n")
                f_r.write(f"Mean Value from the sistem \n")
                f_r.write("="*60 + "\n")
                f_r.write(f"{'Total Stock (tot_y):':<40} {tot_y:.2f} unità\n")
                f_r.write(f"{'Total Lost Sales:':<40} {tot_lost.sum():.2f} unità\n")
                f_r.write(f"{'Client No-Buy (a priori):':<40} {no_buy_medio:.2f}\n")
                f_r.write(f"{'Lost Market Share (a priori) (%):':<40} {(no_buy_medio/domanda_potenziale)*100:.2f}%\n")
                f_r.write(f"{'Total Lost Market Share (%):':<40} {((no_buy_medio+tot_lost.sum())/domanda_potenziale)*100:.2f}%\n")
                f_r.write(f"{'Optimal Profit (Objective Value):':<40} {model.ObjVal:.2f} €\n")
                
                f_r.write("-" * 60 + "\n")
                f_r.write(f"{'Total potential client:':<40} {domanda_potenziale:.2f}\n")
                f_r.write(f"Mean DEMAND for product:\n{media_prodotti}\n")
                f_r.write(f"{'Number of products in the assortment:':<40} {len(indici_assortimento)}\n")
                f_r.write(f"INDICES VECTOR IN THE ASSORTMENT: {indici_assortimento}\n")
                f_r.write(f"SERVICE LEVEL: {service_level_medio}\n")

                f_r.write("\nDETAILED SUBTITUTIONS (mean on the scenarios):\n")
                found = False
                for k in range(len(media_sost)):
                    if media_sost[k]>0.01:
                        found=True
                        i=self.rows_sost[k]
                        j=self.cols_sost[k]
                        p_da=self.df.iloc[i]['Prodotto']
                        p_a=self.df.iloc[j]['Prodotto']
                        f_r.write(f"- [{i:>3}] {p_da:<15} --> [{j:>3}] {p_a:<15}: {media_sost[k]:.2f} unità\n")
                
                if not found:
                    f_r.write("Nessuna sostituzione rilevata.\n")

    def export_to_text_hybrid_det(self,no_buy_a_priori, file_nn_det="analisi_risultati_hybrid_det.txt"):
        """Metodo per report su ottimizzazione ibrido deterministico"""
        model=self.models["hybrid_det"]
       
        y_vals=model._y_hnn.X.flatten()
        a_vals=(model._a_hnn.X > 0.5).astype(int).flatten()

        r_vals=model._r_hnn.X.flatten()
        lost_vals=model._lost_hnn.X.flatten()
        x_vals=model._x_hnn.X.flatten()
        sost_vals=model._sost_hnn.X.flatten()
        SL_vals=model._service_level_hnn.X
        d_pred_vals=model._d_pred_hnn.X.flatten()

        demand_arr=model._demand_saa

        tot_sost_per_prodotto=np.zeros(self.num_opzioni)
        for k, j in enumerate(self.cols_sost):
            tot_sost_per_prodotto[j]+=sost_vals[k]


        df_report=self.df.copy()
        df_report['STOCK (y)']=y_vals
        df_report['REST (r)']=r_vals
        df_report['LOST (l)']=lost_vals
        df_report['PRIMARY DEMAND (x)']=x_vals
        df_report['PREDICTED DEMAND (d)']=d_pred_vals
        df_report['SOST_RECEIVED']=tot_sost_per_prodotto

        indici_assortimento=np.where(a_vals > 0.5)[0].tolist()
        domanda_totale_predetta=d_pred_vals.sum()

        import os
        os.makedirs('test_res', exist_ok=True)
        var_df=pd.DataFrame({'y_values': y_vals, 'a_values': a_vals})
        var_df.to_csv(r'test_res\assortimento_ibrido_det.csv', header=True, index=False)

        with open(file_nn_det, "w") as f_r:
            f_r.write("\n--- HYBRID DETERMINISTIC MODEL ANALYSIS ---\n")
            cols = ['Prodotto', 'Marca', 'Prezzo €/kg', 'Scaffale', 'STOCK (y)', 'REST (r)', 'LOST (l)', 'PRIMARY DEMAND (x)', 'PREDICTED DEMAND (d)', 'SOST_RECEIVED']
            f_r.write(df_report[cols].to_string(index=False, justify='left', float_format="%.2f"))

            f_r.write("\n\n" + "="*60 + "\n")
            f_r.write(f"Hybrid Deterministic System Performance\n")
            f_r.write("="*60 + "\n")
            f_r.write(f"{'Total Stock Allocated (tot_y):':<40} {y_vals.sum():.2f} unità\n")
            f_r.write(f"{'Total Lost Sales:':<40} {lost_vals.sum():.2f} unità\n")
            f_r.write(f"{'Optimal Profit (Objective Value):':<40} {model.ObjVal:.2f} €\n")
            f_r.write(f"{'Total Predicted Market Demand:':<40} {domanda_totale_predetta:.2f} unità\n")
            f_r.write(f"{'Total actual mean Demand:':<40} {demand_arr:.2f} unità\n")
            f_r.write(f"{'Client No-Buy (a priori):':<40} {no_buy_a_priori:.2f}\n")
            f_r.write(f"{'Lost Market Share (a priori) (%):':<40} {(no_buy_a_priori/demand_arr)*100:.2f}%\n")
            f_r.write(f"{'Total Lost Market Share (%):':<40} {((no_buy_a_priori+lost_vals.sum())/demand_arr)*100:.2f}%\n")
            f_r.write(f"{'Optimal Profit (Objective Value):':<40} {model.ObjVal:.2f} €\n")
            f_r.write(f"{'Number of products in the assortment:':<40} {len(indici_assortimento)}\n")
            f_r.write(f"INDICES VECTOR IN THE ASSORTMENT: {indici_assortimento}\n")
            f_r.write(f"SERVICE LEVEL SLACK VARIABLE: {SL_vals:.2f}\n")
            f_r.write(f"{'Gurobi Execution Time (Runtime):':<40} {model.Runtime:.4f} secondi\n")
            
            f_r.write("\nDETAILED SUBSTITUTIONS:\n")
            found = False
            for k in range(len(sost_vals)):
                if sost_vals[k]>0.01:
                    found=True
                    i=self.rows_sost[k]
                    j=self.cols_sost[k]
                    p_da=self.df.iloc[i]['Prodotto']
                    p_a=self.df.iloc[j]['Prodotto']
                    f_r.write(f"- [{i:>3}] {p_da:<15} --> [{j:>3}] {p_a:<15}: {sost_vals[k]:.2f} unità\n")
            
            if not found:
                f_r.write("Nessuna sostituzione rilevata.\n")

    def export_to_text_NN(self,file_nn='analisi_risultati_NN.txt'):
        """Metodo per report su ottimizzazione Surrogato"""
        model=self.models["nn"]
        y_vals=model._y_nn.X.flatten()
        a_vals=(model._a_nn.X > 0.5).astype(int).flatten()
        indici_assortimento=np.where(a_vals>0.5)[0].tolist()
        y_vals_ass=[0]*len(indici_assortimento)
        j=0
        for i in range(len(y_vals)):
            if i in indici_assortimento:
                y_vals_ass[j]=int(y_vals[i])
                j+=1

        f_vector=np.array([self.f_cost[self.df.iloc[i]['Scaffale']] for i in self.I])
        costo_fisso_tot=(f_vector@a_vals)
        costo_acquisto_tot=(self.costs@y_vals)

        df_report=self.df.copy()
        df_report['STOCK (y)']=y_vals
        
        with open(file_nn, "w") as f:
            f.write("--- REPORT OPTIMIZATION SURROGATE MODEL ---\n")
            f.write(f"Data: {pd.Timestamp.now()}\n\n")
            
            f.write("DETAIL ASSORTMENT AND CHOSEN STOCK:\n")
            cols = ['Prodotto', 'Marca', 'Prezzo €/kg', 'Scaffale', 'STOCK (y)']
            f.write(df_report[cols].to_string(index=False, justify='left'))

            f.write("\n\n" + "="*60 + "\n")
            f.write("ECONOMIC BEHAVIOUR ESTIMATE BY THE MODEL\n")
            f.write("="*60 + "\n")

            f.write(f"{'TOTAL PROFIT (Objective Value):':<40} {model.ObjVal:.2f} €\n")
            f.write(f"{'Fixed Shelfs Cost:':<40} {costo_fisso_tot:.2f} €\n")
            f.write(f"{'Stock Purchase Cost:':<40} {costo_acquisto_tot:.2f} €\n")
            f.write(f"{'Total Ordered Stock (units):':<40} {y_vals.sum():.2f}\n")
            
            f.write("-" * 60 + "\n")
            f.write(f"NUMBER OF PRODUCTS IN THE ASSORTIMENT: {len(indici_assortimento)}\n")
            f.write(f"INDICES VECTOR: {indici_assortimento}\n")
            f.write(f"QUANTITIS VECTOR:{y_vals_ass}\n")
  

    def close_model(self,model):
        self.models[model].dispose()
