"""
prepare_data.py — Preparazione dati e generazione scenari

Script che fa da collante tra i dati grezzi e i modelli di ottimizzazione/ML.
Si occupa di:
  1. Caricare e preprocessare il catalogo prodotti da Excel.
  2. Inizializzare la simulazione di mercato (MarketSimulation + LatentChoiceModel + ClientArrivals).
  3. Generare e ridurre gli scenari di domanda stocastica.
  4. Creare i dataset di training per i modelli ML surrogati.

Funzioni presenti:

  import_dataset()          : legge TabellaProdotti.xlsx e separa qualità bassa/alta.
  market_sim_init()         : inizializza tutti gli oggetti di simulazione del mercato.
  prepare_input_data()      : funzione che genera l'intero input per l'ottimizzazione
                              e mette in cache il risultato su file .pkl per evitare ricalcoli.
  prepare_test_input_data() : variante per il test OOS con assortimento già fissato.
  prepare_NN_training_data(): genera dataset (a, y), calcola profitto_ricorso per la surrogate NN,
                              usando una logica exploration-exploitation guidata da Gurobi.
  prepare_arrivals_det()    : genera dataset assortimento, calcola domanda_media per il modello ibrido
                              deterministico, campionando assortimenti casuali.
  prepare_NN_arrivals_stoc(): (espansione futura) versione stocastica multi-scenario del dataset ibrido.

Salvataggio dati:

  prepare_input_data() salva l'output in un file .pkl (default: 'dati_progetto_test.pkl').
  Se il file esiste già, viene caricato senza ricalcolare. Usare force_recalculate=True
  per forzare il ricalcolo (es. dopo aver cambiato i parametri del mercato o il seed).

"""

import pandas as pd
import pickle
import time
import os
from models import *
from optimization import *
import numpy as np
import gurobipy as gp
from gurobipy import *

def import_dataset():
    """funzione per importare il dataset e fare le prime modifiche di base"""

    prodotti_tab=pd.read_excel(r'D:\POLI_Gaia\Materiale_Tesi_Magistrale\Codice\TabellaProdotti.xlsx')
    prodotti_tab.insert(0,'Qualità_bassa',0.0)
    prodotti_tab.insert(1,'Qualità_alta',0.0)
    #Dividiamo i due livelli di qualità:
    #Se la qualità è < 6.8, copia il valore in 'Qualità_bassa'
    prodotti_tab.loc[prodotti_tab['Qualità (0-10)']<6.8,'Qualità_bassa']=prodotti_tab['Qualità (0-10)']
    #Se la qualità è >= 6.8, copia il valore in 'Qualità_alta'
    prodotti_tab.loc[prodotti_tab['Qualità (0-10)']>=6.8,'Qualità_alta']=prodotti_tab['Qualità (0-10)']
    prodotti_tab.drop('Qualità (0-10)',axis=1,inplace=True)
    return prodotti_tab

def market_sim_init(df_prod, mean_arrival=None, seg_prob=None, coeff_rate_arrival=None, use_seed=1234):
    """
    funzione per inizializzare il mercato con le caratteristiche del modello di scelta e gli scenari di domanda

    Input: 
        1. il dataset dei prodotti
        2. il numero di clienti medi attesi giornalmente
        3. vettore con le probabilità dei segmenti di mercato
        4. vettore con i rate di arrivo gionralieri
        5. seed per riproducibilità

    Putput:
        1. istanza della classe MarketSimulation
        2. il dataframe pronto all'uso, con colonne dummy e suddivisioni necessarie 
        3. il modello di scelta
        4. istanza della classe ClientArrivals 
        5. gli arrivi dei clienti giornalieri medi
    """

    if mean_arrival is None:
        mean_arrival=330
    if coeff_rate_arrival is None:
        coeff_rate_arrival=[0.8,0.7,1,0.4,0.8,1.5,0] #come cambia la media nei giorni della settimana
    rate_arrival=[mean_arrival*coeff_rate_arrival[i] for i in range(len(coeff_rate_arrival))] 

    sim=MarketSimulation(df_prod,seg_probs=seg_prob)
    r=sim.sim_prod_prepare()

    #Modello di scelta
    choice_model=sim.use_model(MNLChoiceModel,LatentChoiceModel)
    arrivals=ClientArrivals(rate_arrival,seed=use_seed)
    return sim,r,choice_model,arrivals,rate_arrival

def prepare_input_data(force_recalculate=False, SAA_flag=True, n_scenari_domanda_SAA=100, n_warm_scenarios=0, used_seed=None, nome_file='dati_progetto_test.pkl'):
    """ 
    funzione che prepare le costanti utili per l'ottimizzazione e produce gli scenari per il metodo SAA

    Input: 
        1. flag se si  vuole calcolare l'output o usare quello già esistente
        2. flag se si vogliono ottenere le specifiche degli scenari
        3. numero degli scenari
        4. numero degli scenari di warm-up utili nella generazione del dataset del surrogato (input nella funzione preparare_NN_training_data, chiamata in tarin_NN.py)
        5. seed per riproducibilità
        6. nome del file dove salvare output

    Output: un dizionario con tutte le costanti per il modello sotto la chiave 'base'; tutti gli elementi degli scenari sotto laa chiave 'saa_input'

    Dove viene chiamata:
        1. nei main_*.py dei modelli di ottimizzazione
        2. nei create_*_dataset.py per generare i dataset per il training
        3. in train_NN.py e active_learning.py per il training di NN surrogato
        
    """
    output_esistente=None
    if os.path.exists(nome_file):
        with open(nome_file,'rb') as f:
            output_esistente=pickle.load(f)
   
    if not force_recalculate and output_esistente is not None: 
        print(f"File '{nome_file}' caricato.")
        return output_esistente

    #Se il file non esiste (o force_recalculate=True), eseguiamo tutto il codice
    print("Aggiornamento file...")
    prodotti_tab=import_dataset()

    ########DATI SCELTI################
    prezzi=prodotti_tab['Prezzo €/kg'].values #p_i
    cost=prodotti_tab['Costi €/kg'].values #c_i

    holding=prodotti_tab['Holding €/kg'].values#h_i
    margin=prezzi-cost
    penalty_lost_sales=1.00*margin #la penalità sulle perdite è impostata sul valore del margine di profitto
    penale_SL=2.00 #*np.ones(len(margin)) #theta penalità sul livello di servizio, andiamo a penalizzare con due volte il margine se siamo sotto al livello di servizio
    SL_target=0.95 #SL_{i} livello di servizio legato ai prodotti che vanno out-of-service

    capacita_scaffale_dict={'Alto':480,'Occhi':480,'Basso':360} #C_l SCAFFALI NO FRIGO
    capacita_frigo_dict={'Alto':90,'Occhi':90,'Basso':40} #C_l SCAFFALI FRIGO

    costo_fisso_scaffale={'Alto':6.0,'Occhi':8.0,'Basso':5.0} #f_{i,l} [costo del mettere l'oggetto in assortimento, anche se poi non lo prendo]
    
    prodotti_tab.drop('Costi €/kg',axis=1,inplace=True)
    prodotti_tab.drop('Holding €/kg',axis=1,inplace=True)
    
    sim,r,choice_model,arrivals,rate_arrival=market_sim_init(df_prod=prodotti_tab, use_seed=used_seed)
    
    Prob=choice_model.weighted_Probabilities(r[1])
    sub_matrix=sim.substitution_Matrix(r[1],choice_model)
    sim.substitution_Matrix_bound(sub_matrix)


    if SAA_flag or output_esistente is None:
        n_tot=n_scenari_domanda_SAA+n_warm_scenarios
        #Generatore della domanda fisica (clienti con preferenze) SAA
        demand_weekly_tot,scenario_prob_tot=arrivals.reduced_scenarios(Prob,n_start=200000,n_target=n_tot)
        demand_weekly=demand_weekly_tot[:n_scenari_domanda_SAA,:]
        scenario_prob=scenario_prob_tot[:n_scenari_domanda_SAA]
        scenario_prob=scenario_prob/scenario_prob.sum()
        warm_demand=demand_weekly_tot[n_scenari_domanda_SAA:,:]
        max_revenue = np.zeros(len(margin))
        for i in range(len(margin)):
            max_d_prod = max(demand_weekly[s][i] for s in range(n_scenari_domanda_SAA))
            max_revenue[i] = max_d_prod * prezzi[i]
    else: 
        demand_weekly=output_esistente["saa_input"]["demand_weekly"]
        scenario_prob=output_esistente["saa_input"]["scenario_prob"]
        warm_demand=output_esistente["saa_input"]["warm_demand"]
        max_revenue=output_esistente["saa_input"]["max_revenue"]

    output = {
        "base":{"prodotti_tab":prodotti_tab,
                "rate_arrival":rate_arrival,
                "sub_matrix":sub_matrix,
                "prezzi":prezzi,
                "cost":cost,
                "holding":holding,
                "capacita_scaffale_dict":capacita_scaffale_dict,
                "capacita_frigo_dict":capacita_frigo_dict,
                "costo_fisso_scaffale":costo_fisso_scaffale,
                "penalty_lost_sales":penalty_lost_sales,
                "penale_SL":penale_SL,
                "SL_target":SL_target},
        "saa_input":{
                "demand_weekly":demand_weekly,
                "scenario_prob":scenario_prob,
                "warm_demand": warm_demand,
                "max_revenue": max_revenue}
    }

    with open(nome_file,'wb') as f:
        pickle.dump(output,f)
    
    print(f"Elaborazione completata e salvata in '{nome_file}'.")
    return output

def prepare_test_input_data(assortment, seg_probs=None, used_seed=None):
    """
    funzione variante per il test OOS con assortimento già fissato
    
    Input: 
        1. vettore di assortimento fissato
        2. vettore con la probabilità dei segmenti di mercato
        3. seed per riproducibilità

    Output:
        1. istanza della classe MarketSimulation
        2. le probabilità di acquisto
        3. la matrice di sostituzione
        4. l'etichette dei prodotti unici

    Dove viene chiamata: nei main_test_*.py per i test out-of-sample e negli scenari di shock
    """
    prodotti_tab=import_dataset()

    prodotti_tab.drop('Costi €/kg',axis=1,inplace=True)
    prodotti_tab.drop('Holding €/kg',axis=1,inplace=True)

    name_unique_prod = []

    for i, row in prodotti_tab.iterrows():
        etichetta = f"[{str(i).zfill(3)}] {row['Marca']} | {row['Prodotto']} | {row['Prezzo €/kg']}€/kg | {row['Scaffale']}"
        name_unique_prod.append(etichetta)

    sim,r,choice_model,arrivals,_=market_sim_init(df_prod=prodotti_tab, seg_prob=seg_probs, use_seed=used_seed)

    Prob=choice_model.weighted_Probabilities(r[1],available_prod=assortment)
    
    sub_matrix=sim.substitution_Matrix(r[1],choice_model, available_prod=assortment)

    return arrivals,Prob,sub_matrix,name_unique_prod

def prepare_NN_training_data(optimizer, demand_scenario_structure, demand_scenarios_eval, prob_scen, n_samples=10000, n_extreme=30,seed=1234):
    """
    funzione che crea il dataset di training per il modello surrogato, campionando gli elementi con metodo di exploration-exploitation framework e calcola l'output target usando la funzione di ricorso esatta
    
    Input:
        1. istanza della classe AssortmentOptimizer 
        2. gli scenari di domanda di warm up, ottenuti da prepare_input_data, utilizzati per calcolare le frequenze di selezione dei preodotti più comuni 
        3. gli scenari di domanda per la valutazione esatta della funzione di ricorso, per il target della rete
        4. vettore delle probabilità associate agli scenari di valutazione
        5. numero di campioni standard da generare
        6. numero di campioni estremi da generarre 
        7. seed per la riproducibilità

    Output: dataset finale con dimensione (n_samples + n_extreme, 2*num_p + 1): con 'a' assortimento + 'y' quantità + valore ricorso

    Dove viene usato: in create_NN_dataset.py

    """
    np.random.seed(seed)

    #parametri utili per costruire il modello     
    num_p=demand_scenario_structure.shape[1]-1
    num_unique=num_p//3
    freq_a=np.zeros(num_p)
    n_struct_scenario=demand_scenario_structure.shape[0]

    for i in range(n_struct_scenario): # Risolve i singoli scenari di warm-up per capire quali prodotti sono scelti più spesso, considerando un modello semplice a singolo scenraio
        d_s=demand_scenario_structure[i:i+1,:]
        sol_a,_=optimizer.get_single_scenario_solution(d_s)
        freq_a+=sol_a
    freq_a=freq_a/n_struct_scenario

    optimizer._mean_parameters(demand_scenarios_eval,prob_scen) 

    n_eval_scenarios=demand_scenarios_eval.shape[0]

    row_size=2*num_p+1
    dataset=np.zeros((n_samples+n_extreme,row_size))

    #Costruzione del modello 
    gen_model=gp.Model("Generator")
    gen_model.setParam('OutputFlag',0)
    gen_model.setParam('MIPGap',0.05)
    
    #variabili
    a=gen_model.addMVar(num_p,vtype=GRB.BINARY,name="a")
    y=gen_model.addMVar(num_p,lb=0.0,vtype=GRB.INTEGER,name="y")

    #variabili ausiliarie per linearizzare il valore assoluto: la funzione obiettivo vuole minimizzare la distanza tra a/y con il target_A/y ottenuto dal campionamento, in modo tle da ottenere soluzioni feasibile partendo da valori random
    z_y=gen_model.addMVar(num_p, vtype=GRB.CONTINUOUS, name="z_y") 
    z_a=gen_model.addMVar(num_p, vtype=GRB.CONTINUOUS, name="z_a")
    #variabili ausiliarie per evitare di ricostruire il modello da zero ad ogni nuovo campione: si setta UB e LB uguali al nuovo emlemento campionato
    target_y_placeholder=gen_model.addMVar(num_p, lb=0, name="target_y_placeholder")
    target_a_placeholder=gen_model.addMVar(num_p, lb=0, vtype=GRB.BINARY, name="target_a_placeholder")


    #vincoli
    for i in range(optimizer.num_unique_prod):
        gen_model.addSOS(GRB.SOS_TYPE1, [a[3*i], a[3*i+1]])
    if optimizer.shelf:
        for s_name, indici in optimizer.opzioni_per_scaffale.items():

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


    gen_model.addConstr(y<=200*a) ##200 è LA GRANDE M
  
    #linearizzazione del valore assoluto
    gen_model.addConstr(z_y>=y-target_y_placeholder)
    gen_model.addConstr(z_y>=target_y_placeholder-y)
    gen_model.addConstr(z_a>=a-target_a_placeholder)
    gen_model.addConstr(z_a>=target_a_placeholder-a)

    #f obj
    gen_model.setObjective(z_y.sum()+100*z_a.sum(), GRB.MINIMIZE)

    #ottimizzazione del ricorso
    start_time=time.time()
    optimizer.build_recurse_only(n_scenarios=n_eval_scenarios, write=False)
    model_recurse=optimizer.models["recurse"]

    
    print(f"Inizio generazione di {n_samples+n_extreme} campioni...")

    idx=0
    #casi estremi
    if n_extreme>0:
        #primo campione estremo
        a_null=np.zeros(num_p)
        y_null=np.zeros(num_p)
        val_null=optimizer.solve_recurse(model=model_recurse,a_fixed=a_null,y_fixed=y_null) #valore del target con niente in assortimento
        sample=np.concatenate([a_null,y_null,[val_null]])
        dataset[idx]=sample  

        idx+=1
        #altri campioni estremi (densi)
        for _ in range(n_extreme-1):
            a_full=np.zeros(num_p)

            for i in range(num_unique): #si cerca di generare qualcosa di vicino all'ammissibilità anmche campionando casualmente tenendo conto dei vincoli di mutua esclusione
                a_full[3*i+2]=1

                a_next=np.random.choice([0,1])
                a_full[3*i+a_next]=1
            
            intensity=np.random.uniform(0.3,1)
            max_y_full_campione = int(200*intensity)+1 #si fissa la quantità su livelli alti
            y_full=np.random.randint(0, max_y_full_campione, size=num_p)
            #dati i valori campionati, si usano per fissare i valori delle variabili ausiliarie, così da usare il modello esistente
            target_y_placeholder.setAttr(GRB.Attr.LB, y_full) 
            target_y_placeholder.setAttr(GRB.Attr.UB, y_full)
            target_a_placeholder.setAttr(GRB.Attr.LB, a_full)
            target_a_placeholder.setAttr(GRB.Attr.UB, a_full)

            gen_model.optimize()
            if gen_model.SolCount>0: #trovati dei valori di a e y feasible
                current_a=(a.X>0.5).astype(int)
                current_y=y.X
                
                profitto_ricorso=optimizer.solve_recurse(model=model_recurse, a_fixed=current_a, y_fixed=current_y) #si trovare il valore del ricorso con le variabili trovare
                
                sample=np.concatenate([current_a,current_y,[profitto_ricorso]])
                dataset[idx]=sample
                idx+=1
        print(f"Campioni estremi fatti")   

    #campioni standard: per inserire i valori delle varibili random si fa nello stesso modo, cambia il campionamento
    for i in range(n_samples):
        alpha=np.random.uniform(0,1) #exploration/exploitation
        prob_explore=alpha*freq_a+(1-alpha)*0.5 #con prob aplpha exploit la sol SAA,  altrimenti generazione casuale
        target_a=(np.random.rand(num_p)<prob_explore).astype(int)

        intensita_target=np.random.uniform(0.1, 1.0)
        max_y_campione=int(200*intensita_target)+1 
        target_y=np.random.randint(0, max_y_campione, size=num_p)

        target_y=target_y*target_a #per mettere y=0 quanto a=0

        if np.random.rand()<0.2:
            j=np.random.randint(0,num_unique) 
            
            idx_stress=3*j+np.random.randint(0,3)
            
            target_y[idx_stress]=np.random.randint(190,200)

            target_a[idx_stress]=1
            
            if idx_stress%3==0: 
                target_a[idx_stress+1]=0
            elif idx_stress%3==1: 
                target_a[idx_stress-1]=0
            
            target_y[idx_stress]=target_y[idx_stress]
 
        target_y_placeholder.setAttr(GRB.Attr.LB, target_y)
        target_y_placeholder.setAttr(GRB.Attr.UB, target_y)        
        target_a_placeholder.setAttr(GRB.Attr.LB, target_a)
        target_a_placeholder.setAttr(GRB.Attr.UB, target_a)

        gen_model.optimize()
        
        if gen_model.SolCount>0: #trovati dei valori di a e y feasible
            current_a=(a.X>0.5).astype(int)
            current_y=y.X
            
            profitto_ricorso=optimizer.solve_recurse(model=model_recurse, a_fixed=current_a, y_fixed=current_y)

            sample=np.concatenate([current_a,current_y,[profitto_ricorso]])
            dataset[idx]=sample
            idx+=1
            
        if i%1000==0:
            print(f"Campione {i+1} fatto")

    end_time=time.time()
    elapsed_time=end_time-start_time

    hours,rem=divmod(elapsed_time,3600)
    minutes,seconds=divmod(rem,60)

    print(f"\nGenerazione completata in: {int(hours)}h {int(minutes)}m {seconds:.2f}s")

    return pd.DataFrame(dataset)

def prepare_arrivals_det(class_arr, sim_model, choice_model, data_set, num_samples):

    """
    simula l'arrivo di num_samples e restituisce la domanda prevista che include la domanda primaria
    e la domanda di sostituzione per tutti i prodotti in num_samples in diversi casi di assortimento. Aggiungiamo anche il delta_no_buy
    che è la differenza tra il valore effettivo di no-buy e il valore di riferimento di no-buy per l'intero assortimento.

    Input:
        1. instanza della classe ClientArrivals 
        2. istanza della classe MarketSimulation
        3. instanza della classe LatentChoiceModel
        4. il dataset dei prodotti ripulito, dummizzato pronto all'uso
        5. il numero di campioni da generare

    Output: una tupla con i valori di input per il dataset del training e quelli di target 

    Dove si usa: in create_hybrid_dataset.py
    """
    num_products=data_set.shape[0]-1  #non contanto la no_buy option
    num_unique_prod=num_products//3
    tot_arrival=sum(class_arr.rate_arrival) 

    #calcoliamo il valore della no-buy a priori consideraando il negozio completo
    full_assortment=np.ones(num_products, dtype=int).tolist()
    prob_full=choice_model.weighted_Probabilities(data_set, full_assortment)
    baseline_no_buy=tot_arrival*prob_full[-1]

    #inizializzazione strutture 
    X_det=np.zeros((num_samples,num_products))
    y_det=np.zeros((num_samples,num_products))
    delta_no_buy_det=np.zeros((num_samples,1))


    for s in range(num_samples):
        p=class_arr.rng.uniform(0.05,0.95) #probabilità casuale per variare la densità dell'assortimento ad ogni ciclo
        
        #struttura del vettore che rispetti il vincolo di mutua esclusione (per ogni prodotto unico)
        a=class_arr.rng.binomial(1,p,size=num_unique_prod) #campioniamo la presenza del primo elemento (Alto)
        b=np.zeros(num_unique_prod, dtype=int)
        mask_unique_prod=(a==0) #vediamo se il primo è zero o no
        if np.sum(mask_unique_prod)>0:
            b[mask_unique_prod]=class_arr.rng.binomial(1,p,size=np.sum(mask_unique_prod)) #se il primo è zero campioniamo la presenza del secondo (Basso)
        c=class_arr.rng.binomial(1,p,size=num_unique_prod) #campioniamo la presenza del terzo (Occhi)
        #costruiamo l'assortimento
        actual_assortment=np.column_stack((a,b,c)).flatten()
        assortment_list=actual_assortment.tolist()
        probabilities=choice_model.weighted_Probabilities(data_set, assortment_list)
        subs_matrix=sim_model.substitution_Matrix(data_set, choice_model, assortment_list)
        sim_model.substitution_Matrix_bound(subs_matrix)
        prob_choice=class_arr._total_demand_flow(probabilities, subs_matrix,assortment_list) #valutiamo probabilità complessiva

        exp_demand=tot_arrival*prob_choice[:-1] #domanda attesa per ogni prodotto
        X_det[s,:]=actual_assortment #assortimento valutato
        y_det[s,:]=exp_demand
        no_buy_actual=tot_arrival*prob_choice[-1]
        delta_no_buy_det[s,0]=no_buy_actual-baseline_no_buy #no_buy aggiunta rispetto a quella a catalogo pieno

    return (X_det,y_det,delta_no_buy_det) 


#PER ESTENSIONE FUTURA
def prepare_NN_arrivals_stoc(class_arr, sim_model, choice_model, data_set, num_samples,num_scenarios=None):

    """
    Espansione futura: Modello ibrido stocastico
    Questa funzione simula num_samples x num_s arrivi e restituisce le probabilità di domanda, per tutti gli scenari, che contengono la domanda primaria
    e di sostituzione per tutti i prodotti in num_samples in diversi casi di assortimento.

    Utilizza la distribuzione multinomiale per ottenere il numero effettivo di prodotti necessari in ogni scenario. Uno scenario è descritto dal
    numero di clienti che entrano nel sistema.

    È utile creare il dataset per addestrare un predittore per tutti gli scenari al fine di
    calcolare il modello stocastico in Gurobi. In questo caso, l'obiettivo è prevedere la domanda dello scenario
    e utilizzarla nel problema di ottimizzazione come domanda effettiva per quello scenario.
    """
    num_products=data_set.shape[0]-1  #non contanto la no_buy option
    num_unique_prod=num_products//3
    num_days=len(class_arr.rate_arrival)
    n_arri=class_arr._only_arrivals(num_scenarios=num_scenarios)

    n_arri_week=n_arri.sum(axis=1)
    full_assortment=np.ones(num_products, dtype=int).tolist()
    prob_full=choice_model.weighted_Probabilities(data_set, full_assortment)  

    X_stoch=[]
    y_stoch=[]
    delta_no_buy_stoch=[]

    for _ in range(num_samples):
        p=class_arr.rng.uniform(0.05,0.95)
        a=class_arr.rng.binomial(1,p,size=num_unique_prod)
        b=np.zeros(num_unique_prod, dtype=int)
        mask_unique_prod=(a==0)
        if np.sum(mask_unique_prod)>0:
            b[mask_unique_prod]=class_arr.rng.binomial(1,p,size=np.sum(mask_unique_prod))
        c=class_arr.rng.binomial(1,p,size=num_unique_prod)

        actual_assortment=np.column_stack((a,b,c)).flatten()
        assortment_list=actual_assortment.tolist()
        probabilities=choice_model.weighted_Probabilities(data_set, assortment_list)
        subs_matrix=sim_model.substitution_Matrix(data_set, choice_model, assortment_list)
        sim_model.substitution_Matrix_bound(subs_matrix)
        prob_choice=class_arr._total_demand_flow(probabilities, subs_matrix,assortment_list)

        if num_scenarios is not None:
            demand_scenarios=np.zeros((num_scenarios, num_days, num_products+1), dtype=int)
            for sc in range(num_scenarios):
                baseline_no_buy=n_arri_week[sc]*prob_full[-1]
                for d in range(num_days):
                    demand_scenarios[sc, d, :]=class_arr.rng.multinomial(n_arri[sc,d],prob_choice)
                weekly_demand_scenario=demand_scenarios[sc,:,:-1].sum(axis=0)
                weekly_no_buy_scenario=demand_scenarios[sc,:,-1].sum()
                arrivi_totali = int(n_arri[sc,:].sum())
                X_comb=actual_assortment.tolist()+[arrivi_totali]

                X_stoch.append(X_comb)
                y_stoch.append(weekly_demand_scenario)
                delta_no_buy_stoch.append([weekly_no_buy_scenario-baseline_no_buy])
            
    X_stoch_fin=np.array(X_stoch)
    y_stoch_fin=np.array(y_stoch)
    delta_no_buy_stoch_fin=np.array(delta_no_buy_stoch)


    return (X_stoch_fin, y_stoch_fin, delta_no_buy_stoch_fin)   


