
# ============================================================
# DESCRIPTION DU FICHIER
# Ce fichier permet de simuler les effets de la réforme "baisse uniforme du
# montant d'allocation" (baisse de 1% du taux de remplacement pour l'ensemble
# des demandeurs d'emploi) à partir des paramètres calibrés sur le régime de
# référence (cf. première partie du notebook simulations_reformes et calibration_modele_simule.py).

# Trois composantes peuvent être simulés via la fonction f_niveau_uniforme :
#   - contrefactuel (pas de réforme) : b_pre = b_post, élasticités à 0
#   - mécanique : b_post < b_pre, élasticités à 0 (comportements inchangés)
#   - total : b_post < b_pre, élasticités calibrées (réponses comportementales, il est possible de modifier les élasticités)
# L'effet comportemental s'obtient par différence : total - mécanique.
# La réforme s'applique au flux à partir du mois 240 et ne modifie pas les
# droits déjà ouverts (cf. Focus, section "Simulation des réformes").

# Utilisé par : simulations_reformes.ipynb (partie Réformes du montant)
# Fichiers produits : deux fichiers csv optionnels retraçant respectivement
#   - file1 : agrégats sur les 12 derniers mois et différences avec l'état
#             stationnaire pré-réforme (court terme et long terme)
#   - file2 : trajectoire mensuelle complète post-réforme (mois 240 à 480)
#             utilisée pour visualiser la montée en charge dans simulations_reforems.ipynb
# ============================================================

# packages
import numpy as np
from scipy.optimize import minimize
import pandas as pd
from concurrent.futures import ProcessPoolExecutor, as_completed




# -------------------------------------------------------------------
# CLASSE INDIVIDU AVEC PLUSIEURS FLUX DE NOMBRES ALÉATOIRES
# -------------------------------------------------------------------

from numpy.random import SeedSequence, default_rng

# ============================================================
# PARAMÈTRES FIXES ET STATISTIQUES CIBLES
# Ces paramètres sont définis directement ici car certains sont spécifiques à
# cette réforme (en particulier le salaire moyen et donc le niveau de
# cotisation, qui dépendent de la population affectée par la réforme).
# ============================================================

# Paramètres financiers (calibration MiDAS)
salaire_pre = 1968 # salaire mensuel moyen avant entrée au chômage pour la population affectée par la réforme baisse uniforme (calibration MiDAS)
tax = salaire_pre*(4+1.47*0.9825 + 0.89)/100 # cotisations mensuelles moyennes : 4% (cotisation employeur AC) + 1,47% × 0,9825 (CSG assise sur 98,25% du brut) + 0,89% (taux fictif imputant le financement des minima sociaux), soit 6,33% du salaire brut
b_rsa = 165.86 # niveau moyen des minima sociaux (RSA et ASS) en prenant en compte le taux de recours
rsa_eligibility_rate = 1 # fixé à 1 car le taux de recours est déjà pris en compte dans b_rsa

# Paramètres de la simulation
months_sim = 480 # durée de la simulation en mois (étendue par rapport à la calibration : 240 mois pour atteindre l'état stationnaire pré-réforme + 240 mois post-réforme pour observer la montée en charge complète)
economy_scale = 10 # 1 individu dans la simulation représente economy_scale individus dans l'économie réelle
nbr_threads = 10 # nombre de coeurs utilisés dans la parallélisation

# Paramètres de l'assurance chômage
# NB : ces paramètres sont structurés en listes [valeur_pré_réforme, valeur_post_réforme]
# pour garder une structure commune avec les codes des autres réformes (éligibilité,
# durée des droits). Dans le cas de la baisse uniforme du montant, ils sont
# identiques pré/post car la réforme ne porte que sur le niveau d'allocation b.

threshold_UI = [6,6] # seuil d'éligibilité (durée d'affiliation minimale) pré/post-réforme
pra_max = [24, 24] # période de référence affiliation pré/post-réforme
pbd_max = [18, 18] # durée potentielle d'indemnisation maximale pré/post-réforme
coef_conversion_UI = [0.75,0.75] # coefficient de conversion utilisé dans le calcul de la durée des droits pré/post-réforme

# Statistiques cibles (voir parametres_regime_reference.py pour le détail)
nbr_individuals_emp = 19341366/economy_scale # nombre d'individus en emploi
target_UI_stock = (2428884)/economy_scale # nombre d'individus au chômage indemnisé
tot_pop_guess = nbr_individuals_emp  + target_UI_stock+ (3279392)/economy_scale # Population totale : nombre d'individus en emploi + nombre d'individus au chômage indemnisé + nombre d'individus au chômage non indemnisé
target_reprise_inflow = (1345603/12)/economy_scale # Nombre de reprises de droits
target_recharg_inflow = (336752/12)/economy_scale # nombre de rechargements
target_UI_inflow = ((1895926)/12)/economy_scale # ouvertures de droits initiales et autres (Unedic, 2024) - ODD intermittents et 62+
target_UI_outflow = (2279654/12)/economy_scale # sorties courtes + longues
target_share68_inflow = 0.028192241	+ 0.026441572 # part des ouvertures de droits avec 6 ou 7 mois de durée d'affiliation
target_share68_inflow_odd = 73570/1302477 # part des ouvertures de droits avec 6 à 8 mois de durée d'affiliation (calibration MiDAS)
flow_system_tot = 1/12* 572727/economy_scale # nombre d'individus sortant/entrant sur le marché du travail (calibré pour un renouvellement complet en 44 ans, 62-18)
flow_system =  flow_system_tot/tot_pop_guess # taux de renouvellement de la population active

target_UI_epuise = (745839/12)/economy_scale # Nombre d'épuisements de droits (non utilisé dans l'optimisation)


# ============================================================
# CLASSE INDIVIDU
# Représente un agent dans la simulation du marché du travail.
# Chaque individu possède un statut (employé, chômeur indemnisé,
# chômeur non indemnisé), une éligibilité à l'assurance chômage,
# des taux de transition emploi/chômage, et un historique de
# mois travaillés sur les 72 derniers mois.
# Différences avec la classe Individual du fichier de calibration :
#   - reçoit en argument le vecteur b = [b_pre, b_post] pour pouvoir
#     calculer la taille du choc relatif sur l'allocation
#   - possède un indicateur ref (0 = régime pré-réforme, 1 = régime
#     post-réforme) qui détermine si l'individu est affecté par la
#     réforme (ref passe à 1 quand l'individu n'a plus de droit ouvert/résiduel
#     dans l'ancien régime à un moment après le mois 240)
#   - possède un compteur month_simu pour repérer le mois courant
#     de la simulation et déclencher la réforme au mois 240
# ============================================================
class Individual:

    # ------------------------------------------------------------------
    # Fonction initialisant les caractéristiques de l'individu
    # ------------------------------------------------------------------
    def __init__(self, id, b, status='employed', eligible=False,
                 job_finding_rate=0.1, job_destruction_rate=0.01,
                 take_up_rate=0.5, months_worked=None, months_ui_left=0,
                 take_up=0, jdr_coef=1.0,
                 rngs=None):
        
        # Initialisation de l'historique de travail : 72 mois à 0 (= non travaillés) par défaut
        if months_worked is None:
            months_worked = [0] * 72 
        self.id = id
        self.status = status           # Statut courant : 'employed', 'ui_unemployed', 'noui_unemployed' (emploi, chômage indemnisé, chômage non indemnisé)
        self.eligible = eligible       # Éligibilité à l'assurance chômage
        self.months_worked = months_worked  # Historique binaire d'emploi des 72 derniers mois (1=travaillé, 0=non travaillé)
        self.tot_months_worked = sum(self.months_worked[0:pra_max[0]])  # Nombre total de mois travaillés sur la période de référence (régime pré-réforme à l'initialisation)
        self.months_ui_left = months_ui_left  # Mois de droits au chômage restants
        self.take_up = take_up         # Indicateur de recours effectif à l'assurance chômage (1=recours, 0=non-recours)
        self.take_up_rate = take_up_rate  # Taux de recours à l'assurance chômage (probabilité)
        self.jdr_coef = jdr_coef       # Coefficient de sensibilité du taux de destruction d'emploi à l'ancienneté

        # Marqueurs de flux mensuels pour calculer les stocks/flux agrégés (remis à 0 à chaque mise à jour)
        self.flow_into_ui = 0     # Indicateur d'une entrée en chômage indemnisé (ouverture de droits initiale)
        self.flow_out_ui = 0      # Indicateur d'une sortie du chômage indemnisé (retour à l'emploi)
        self.flow_recharg = 0     # Indicateur d'un rechargement des droits au chômage
        self.flow_reprise = 0     # Indicateur d'une reprise de droits préexistants après perte d'emploi
        self.odd = 0              # Indicateur d'une ouverture de droits (initiale ou rechargement)
        self.months_odd = 0       # Nombre de mois travaillés ayant donné lieu à l'ouverture de droits
        self.epuise_uni = 0       # Indicateur d'épuisement des droits
        self.delai_months_ui_left = 0  # Délai restant avant déchéance des droits
        self.month_simu = 0       # Compteur du mois courant de la simulation pour cet individu (utilisé pour déclencher la réforme au mois 240)
        self.ref = 0              # Indicateur du régime applicable à l'individu : 0 = régime pré-réforme (droits ouverts avant le mois 240), 1 = régime post-réforme (droits ouverts/rechargés à partir du mois 240) ; seuls les individus en ref=1 subissent le choc d'élasticité

        # Générateurs de nombres aléatoires indépendants par type d'événement
        # (permet la reproductibilité du code et l'indépendance des chocs)
        self.rngs = rngs or {
            'type': np.random.default_rng(),         # NB : ce stream 'type' est créé mais n'est pas utilisé dans cette réforme 
            'destruction': np.random.default_rng(),  # Destruction d'emploi
            'finding': np.random.default_rng(),       # Retour à l'emploi
            'system': np.random.default_rng(),        # Renouvellement de la population
            'take_up': np.random.default_rng(),       # Décision de recours à l'assurance chômage
        }

        
        self.job_finding_rate = job_finding_rate        # Taux mensuel de retour à l'emploi
        self.job_destruction_rate = job_destruction_rate  # Taux mensuel de destruction d'emploi

    # ------------------------------------------------------------------
    # Fonction mettant à jour l'état de l'individu et ses caractéristiques (historique d'emploi, droits restants) chaque mois       # en fonction des probabilités de transition, qui sont elles mêmes dépendantes de ses caractéristiques (état en début de       # mois, historique d'emploi, droits restants), et des règles de l'assurance chômage (éligibilité, durée des droits etc.).
    # Différences avec la calibration :
    #   - reçoit en argument b (vecteur [b_pre, b_post]) et les élasticités e_jdr, e_jfr
    #   - calcule la taille du choc relatif sur l'allocation et applique les
    #     élasticités aux taux de transition pour les individus soumis au régime
    #     post-réforme (ref=1)
    #   - passe l'individu en ref=1 à partir du mois 240 dès qu'il a consommé
    #     ses éventuels droits résiduels (réforme qui
    #     s'applique au flux et non au stock)
    # ------------------------------------------------------------------
    def update_status(self, b, e_jdr, e_jfr):
        # Calcul de la taille du choc d'allocation (variation relative entre b_pre et b_post)
        # choc[0] = 0 pour les individus en régime pré-réforme (non affectés)
        # choc[1] = (b_post - b_pre)/b_pre pour les individus en régime post-réforme (affectés)
        # Cette structure permet d'appliquer le choc uniquement aux affectés via l'indexation par self.ref
        choc_size = (b[1] - b[0])/b[0]
        choc = [0, choc_size]
        self.month_simu += 1
        # Passage en régime post-réforme : un individu bascule en ref=1 à partir
        # du mois 240 dès qu'il ne lui reste plus de droits ouverts pré-réforme
        # (la réforme s'applique au flux, pas au stock)
        if self.month_simu >= 240 and self.months_ui_left < 1:
            self.ref = 1

        # Remise à zéro des marqueurs de flux mensuels
        self.flow_into_ui = 0
        self.flow_recharg = 0
        self.flow_reprise = 0
        self.flow_out_ui = 0
        self.odd = 0
        self.months_odd = 0
        self.epuise_uni = 0

        # Calcul du nombre total de mois travaillés sur la période de référence applicable
        # (indexée par self.ref : ici identique en pré et post réforme du montant)
        self.tot_months_worked = sum(self.months_worked[0:pra_max[self.ref]])
        # Décompte du délai avant déchéance des droits (si non nul)
        if self.delai_months_ui_left >= 1:
            self.delai_months_ui_left -= 1

        # Renouvellement de la population : avec probabilité flow_system,l'individu est "réinitialisé" (simule les entrées/sorties du marché du travail, calibré de manière approximative avec un renouvellement en 44 ans (62-18)
        if self.rngs['system'].random() < flow_system:
            self.months_ui_left = 0
            self.months_worked = [0] * 72 
            self.take_up = 0
            # Un chômeur indemnisé affecté par la réinitialisation passe au chômage non indemnisé (car un nouvel entrant sur le marché du travail n'a pas de droits au chômage restants)
            if self.status == 'ui_unemployed':
                self.status = 'noui_unemployed'

        # ----------------------------------------------------------------
        # Cas 1 : Individu en emploi
        # ----------------------------------------------------------------
        if self.status == 'employed':
            self.take_up = 0
            # Si le délai de déchéance est écoulé, les droits restants sont annulés
            if self.delai_months_ui_left == 0:
                self.months_ui_left = 0
                # Si on est après la date d'application de la réforme cet individu devient soumis au régime réformé car il n'a plus de droits résiduels
                if self.month_simu >= 240:
                    self.ref = 1
            # Mise à jour de l'historique : marque le mois courant comme travaillé car l'individu est en emploi
            self.months_worked = [1] + self.months_worked[:-1]
            self.tot_months_worked = sum(self.months_worked[0:pra_max[self.ref]]) # Calcul du nombre total de mois travaillés sur la période de référence
            # Calcul de l'ancienneté continue (nombre de mois travaillés sans interruption)
            months_worked_continu = min(self.months_worked.index(0), 72) if 0 in self.months_worked else 72 # plafonne à la durée de la mémoire (= 5 ans) 

            # Vérification de l'éligibilité à une indemnisation de l'assurance chômage
            if self.months_ui_left > 0:
                # Si l'individu a encore des droits résiduels : éligible et son taux de séparation est affecté par la réforme (s'il est en régime post-réforme)
                self.eligible = True
                # Taux de destruction d'emploi : taux de base modulé par le choc et l'élasticité (1+choc*e_jdr) , taux décroissant avec l'ancienneté (jdr_coef < 0). e_jdr habituellement > 0 donc une baisse de l'allocation (choc<0) fait baisser le taux de destruction d'emploi à ancienneté donnée
                # choc[self.ref] vaut 0 pour les non-affectés (régime pré-réforme) donc le taux est inchangé pour eux
                jdrate = self.job_destruction_rate*(1+choc[self.ref]*e_jdr)*np.exp(self.jdr_coef*months_worked_continu)     
            else:
                if self.tot_months_worked >= threshold_UI[self.ref]:
                    # Si l'individu a travaillé suffisamment pour être éligible : son taux de séparation est affecté par la réforme (s'il est en régime post-réforme)
                    self.eligible = True
                    jdrate = self.job_destruction_rate*(1+choc[self.ref]*e_jdr)*np.exp(self.jdr_coef*months_worked_continu)
                else:
                    # Si l'individu n'a pas de droit restant et n'a pas assez travaillé pour être éligible : pas affecté par la réforme (ne peut pas ouvrir de droit donc le niveau d'allocation ne modifie pas ses incitations)
                    self.eligible = False
                    jdrate = self.job_destruction_rate*np.exp(self.jdr_coef*months_worked_continu) # pas affecté car ne peut pas ouvrir de droit

            # Tirage du choc de destruction d'emploi
            if self.rngs['destruction'].random() < jdrate: # L'individu perd son emploi avec une probabilité jdrate
                if self.eligible: # Si l'individu est éligible à l'assurance chômage
                    self.status = 'noui_unemployed'
                    if self.months_ui_left == 0 and self.rngs['take_up'].random() < self.take_up_rate:
                        # Ouverture de droits : l'individu n'avait pas de droits résiduels
                        # et décide de recourir à l'assurance chômage (probabilité self.take_up_rate)
                        self.flow_into_ui = 1
                        self.status = 'ui_unemployed'
                        self.odd = 1
                        self.take_up = 1
                        self.months_odd = self.tot_months_worked
                        months_worked_act = self.months_worked[0:pra_max[self.ref]] # vecteur des périodes travaillées sur la période de référence
                        # Longueur de la PRC (durée entre le premier et dernier mois travaillé sur la période de référence)
                        length_prc = (pra_max[self.ref] - months_worked_act[::-1].index(1)) if 1 in months_worked_act else 0
                        # Calcul des droits ouverts : entre 6 mois et pbd_max, selon la formule de conversion
                        rights = min(pbd_max[self.ref], max(round(min(self.tot_months_worked*1.7, length_prc)*coef_conversion_UI[self.ref]),6))
                        self.rights_opened = rights # NB : cet attribut est calculé mais n'est pas réutilisé dans la suite (héritage d'un code partagé avec d'autres réformes)
                        self.months_ui_left = rights
                        # Le délai de déchéance court sur 36 mois + durée des droits
                        self.delai_months_ui_left = 36 + rights
                        self.months_worked = [0]*72 
                    else:
                        # Reprise de droits résiduels si l'individu en avait encore
                        if self.months_ui_left > 0:
                            self.status = 'ui_unemployed'
                            self.flow_reprise = 1
                            self.take_up = 1
                else:
                    # Perte d'emploi sans éligibilité : chômage non indemnisé
                    self.status = 'noui_unemployed'
                    self.take_up = 0

        # ----------------------------------------------------------------
        # Cas 2 : Individu au chômage (indemnisé ou non)
        # ----------------------------------------------------------------
        else:
            self.months_worked = [0] + self.months_worked[:-1] # marque le mois courant comme non travaillé
            self.tot_months_worked = sum(self.months_worked[0:pra_max[self.ref]]) # calcule le nombre de mois travaillés sur la période de référence
            # Sous-cas 2a : Chômeur indemnisé (recours effectif à l'assurance chômage)
            if self.take_up > 0:
                self.months_ui_left = max(self.months_ui_left - 1, 0) # consommation d'un mois de droits
                # Taux de retour à l'emploi : taux de base modulé par l'ampleur du choc et l'élasticité (1+choc*e_jfr) pour les individus en régime post-réforme
                # Une baisse de l'allocation (choc < 0) couplée à e_jfr < 0 augmente le taux de retour à l'emploi (effet incitatif)
                jfr_ui = self.job_finding_rate*(1+choc[self.ref]*e_jfr) 
                if self.rngs['finding'].random() < jfr_ui: # si retrouve un emploi
                    self.status = 'employed'
                    self.take_up = 0
                    self.flow_out_ui = 1
                else: # si ne retrouve pas d'emploi
                    if self.months_ui_left == 0: # si est arrivé au bout de son droit
                        # Une fois les droits épuisés et le mois 240 atteint, bascule en régime post-réforme
                        if self.month_simu >= 240:
                            self.ref = 1
                        self.tot_months_worked = sum(self.months_worked[0:pra_max[self.ref]])
                        # Vérification de la durée retravaillée utilisable pour un rechargement :
                        if sum(self.months_worked) == 0 or self.months_worked.index(1) >= 48: # si pas retravaillé
                            self.tot_months_worked_recharg = 0
                        else:
                            # Calcul des mois travaillés sur la période de référence du rechargement
                            self.tot_months_worked_recharg = sum(self.months_worked[self.months_worked.index(1):self.months_worked.index(1)+pra_max[self.ref]]) # regarde le nombre de mois travaillés sur la PRA associée
                        if self.tot_months_worked_recharg >= threshold_UI[self.ref]: # si a retravaillé au moins 6 mois => ouvre un nouveau droit (on considère un taux de recours de 1 ici notamment grâce aux dispositifs automatiques de rappel)
                            # Rechargement des droits au chômage
                            self.flow_recharg = 1
                            self.odd = 1
                            self.months_odd = self.tot_months_worked
                            months_utiles = self.months_worked[self.months_worked.index(1):self.months_worked.index(1)+pra_max[self.ref]]
                            length_prc = (pra_max[self.ref] - months_utiles[::-1].index(1)) if 1 in self.months_worked else 0
                            # Calcul des nouveaux droits rechargés
                            rights = min(pbd_max[self.ref], max(round(min(self.tot_months_worked_recharg*1.7, length_prc)*coef_conversion_UI[self.ref]),6))
                            self.months_ui_left = rights
                            self.rights_opened = rights # NB : attribut non réutilisé dans la suite (cf. plus haut)
                            self.delai_months_ui_left = 36 + rights
                            self.months_worked = [0]*72 
                        else: # si n'a pas suffisamment retravaillé : sortie du chômage pour épuisement de droits, passe au chômage non indemnisé
                            self.epuise_uni = 1
                            self.status = 'noui_unemployed'
                            self.take_up = 0
                            self.eligible = False
            # Sous-cas 2b : Chômeur non indemnisé
            else:
                if self.rngs['finding'].random() < self.job_finding_rate: # Si il retrouve un emploi (probabilité supposée identique à celle d'un chômeur indemnisé dans ce modèle) - le taux n'est pas affecté par la réforme car cet individu n'est pas indemnisé
                    self.status = 'employed'
                    self.take_up = 0
                # sinon pas de changement d'état


# ===================================================================
# CLASSE SIMULATION
# Orchestre la simulation sur un ensemble d'individus et un nombre
# de mois donné. Agrège les statistiques mensuelles dans un DataFrame.
# Différences avec la calibration :
#   - reçoit en argument b (vecteur [b_pre, b_post]) et les élasticités
#     e_jdr, e_jfr, qui sont propagés aux mises à jour individuelles
#   - agrège les statistiques supplémentaires liées à la réforme :
#     ui_unemployed_regime_new (chômeurs indemnisés sous régime post-réforme),
#     ind_regime_new (tous individus en régime post-réforme)
#   - calcule le coût de l'assurance chômage en distinguant les individus
#     indemnisés à l'ancien montant b_pre (régime pré-réforme) et au nouveau
#     montant b_post (régime post-réforme)
# ===================================================================
  

class Simulation:
    # ------------------------------------------------------------------
    # Fonction initialisant les individus
    # ------------------------------------------------------------------
    def __init__(self, b, e_jdr, e_jfr, num_individuals=100, months=months_sim,
                 job_finding_rate=0.1, take_up_rate=0.3, job_destruction_rate= 0.01,
                 jdr_coef=1.0, base_seed=None):

        # Utilisation d'une SeedSequence pour garantir la reproductibilité
        # et l'indépendance des flux aléatoires entre individus
        ss_master = SeedSequence(base_seed)
        individual_seeds = ss_master.spawn(num_individuals)

        self.num_individuals = num_individuals
        self.months = months
        self.individuals = []

        for i in range(num_individuals):
            # Création de 5 flux RNG indépendants par individu
            # NB : le flux 'type' est créé mais non utilisé dans cette réforme (cf. classe Individual)
            ss_ind = individual_seeds[i].spawn(5)
            rngs = {
                'type': default_rng(ss_ind[0]),
                'destruction': default_rng(ss_ind[1]),
                'finding': default_rng(ss_ind[2]),
                'system': default_rng(ss_ind[3]),
                'take_up': default_rng(ss_ind[4]),
            }

            ind = Individual(
                id=i,
                b = b,
                job_finding_rate=job_finding_rate,
                take_up_rate=take_up_rate,
                job_destruction_rate=job_destruction_rate,
                jdr_coef=jdr_coef,
                rngs=rngs
            )
            self.individuals.append(ind)

        # DataFrame de stockage de l'historique de la simulation
        self.history = pd.DataFrame(columns=['month', 'employed', 'ui_unemployed', 'noui_unemployed'])
    # ------------------------------------------------------------------
    # Fonction exécutant la simulation chaque mois et enregistrant les statistiques agrégées
    # ------------------------------------------------------------------
    def run(self, b, e_jdr, e_jfr):
        records = []
        # Exécute la simulation sur tous les mois
        for month in range(self.months):
            # ----------------------------------------------------------
            # Calcul des stocks et flux agrégés en début de mois
            # ----------------------------------------------------------
            employed = sum(1 for ind in self.individuals if ind.status == 'employed')
            employed_eligible = sum(1 for ind in self.individuals if ind.status == 'employed' and ind.eligible)
            ui_unemployed = sum(1 for ind in self.individuals if ind.status == 'ui_unemployed')
            # Chômeurs indemnisés sous régime post-réforme (statistique supplémentaire spécifique aux simulations de réforme : permet de calculer la part de bénéficiaires affectés par la réforme et de suivre la montée en charge)
            ui_unemployed_regime_new = sum(1 for ind in self.individuals if ind.status == 'ui_unemployed' and ind.ref == 1)
            # Ensemble des individus passés en régime post-réforme (tous statuts confondus)
            ind_regime_new = sum(1 for ind in self.individuals if ind.ref == 1)
            noui_unemployed = sum(1 for ind in self.individuals if ind.status == 'noui_unemployed')
            noui_eligible = sum(1 for ind in self.individuals if ind.status == 'noui_unemployed' and ind.eligible)
            # Employés ayant entre 6 et 8 mois d'ancienneté (segment d'intérêt pour le calage)
            employed_6_8_history = sum(1 for ind in self.individuals if ind.status == 'employed' and 6 <= ind.tot_months_worked < 8)
            # Agrégation des flux mensuels
            flow_out_ui = sum(ind.flow_out_ui for ind in self.individuals)
            flow_into_ui = sum(ind.flow_into_ui for ind in self.individuals)
            flow_recharg = sum(ind.flow_recharg for ind in self.individuals)
            flow_reprise = sum(ind.flow_reprise for ind in self.individuals)
            epuise_uni = sum(ind.epuise_uni for ind in self.individuals)
            odd = sum(ind.odd for ind in self.individuals) # Total des ouvertures/rechargements de droits
            # Part des ouvertures de droits pour les individus ayant 6-8 mois travaillés
            odd_6_8 = sum(1 for ind in self.individuals if ind.odd == 1 and 6 <= ind.months_odd < 8)
            tax_rev = tax * employed # cotisations
            # Coût d'assurance chômage : les indemnisés du régime post-réforme perçoivent b_post=b[1], les autres b_pre=b[0]
            ui_cost = ui_unemployed_regime_new*b[1]+ (ui_unemployed-ui_unemployed_regime_new)*b[0]
            rsa_cost = b_rsa * noui_unemployed * rsa_eligibility_rate # minima sociaux
            net_budget = tax_rev - ui_cost - rsa_cost # dépense mensuelle nette du système

            # Enregistrement des statistiques du mois
            records.append({
                'month': month,
                'employed': employed,
                'employed_eligible': employed_eligible,
                'ui_unemployed': ui_unemployed,
                'ui_unemployed_regime_new': ui_unemployed_regime_new,
                'ind_regime_new': ind_regime_new,
                'noui_unemployed': noui_unemployed,
                'noui_eligible': noui_eligible,
                'employed_6_8_history': employed_6_8_history,
                'flow_out_ui': flow_out_ui,
                'flow_into_ui': flow_into_ui,
                'flow_recharg': flow_recharg,
                'flow_reprise': flow_reprise,
                'odd': odd,
                'odd_6_8': odd_6_8,
                'epuise_uni': epuise_uni,
                'tax_rev': tax_rev,
                'ui_cost': ui_cost,
                'rsa_cost': rsa_cost,
                'net_budget': net_budget
            })
            # Mise à jour des individus (propagation de b, e_jdr, e_jfr à chaque update individuel)
            for ind in self.individuals:
                ind.update_status(b, e_jdr, e_jfr)
        self.history = pd.concat([self.history, pd.DataFrame(records)], ignore_index=True)

    def get_history(self):
        return self.history


# -------------------------------------------------------------------
# Fonction exécutant la simulation pour une population d'individus donnée
# Différence avec la calibration : propage en plus les paramètres de réforme
# (b, e_jdr, e_jfr) à la classe Simulation
# -------------------------------------------------------------------
def run_one(sim_id, b, e_jdr, e_jfr, num_individuals=10000, months=months_sim, seed=None,
            job_finding_rate=0.1, take_up_rate=0.9, job_destruction_rate=0.01,
            jdr_coef=1.0):
    # Lance une simulation unique identifiée par sim_id et retourne son historique.
    sim = Simulation(b = b, e_jdr = e_jdr, e_jfr=e_jfr, num_individuals=num_individuals, months=months,
                     job_finding_rate=job_finding_rate, take_up_rate=take_up_rate,
                     job_destruction_rate=job_destruction_rate,
                     jdr_coef=jdr_coef, base_seed=seed)
    sim.run(b, e_jdr, e_jfr)
    hist = sim.get_history().copy()
    hist["run_id"] = sim_id # Ajout de l'identifiant du thread pour le suivi
    return hist


# -------------------------------------------------------------------
# Fonction de simulation parallélisée
# Permet de paralléliser la simulation sur plusieurs coeurs
# pour réduire le temps de calcul sur de grandes populations.
# Différence avec la calibration : propage en plus les paramètres de réforme
# (b, e_jdr, e_jfr) à chaque thread
# -------------------------------------------------------------------
def run_all(b, e_jdr, e_jfr, job_finding_rate= 0.1, job_destruction_rate= 0.01, take_up_rate=0.3, months=months_sim,
            jdr_coef=1.0):
    # Lance nbr_threads simulations en parallèle et agrège les résultats par mois.

    # Répartition de la population totale entre les threads
    tot_pop = tot_pop_guess
    X = nbr_threads
    num_individuals_thread = int(tot_pop / X)
    # Graines différentes par thread pour garantir l'indépendance des simulations
    base_seed = 12345
    seeds = [base_seed + i * 1000 for i in range(X)]

    futures = []
    results = []
    # Soumission des tâches au pool de processus
    with ProcessPoolExecutor(max_workers=X) as ex:
        for i in range(X):
            futures.append(ex.submit(run_one, i, b, e_jdr, e_jfr, num_individuals_thread, months, seeds[i],
                                     job_finding_rate, take_up_rate, job_destruction_rate,
                                     jdr_coef))
        # Collecte des résultats au fur et à mesure de leur complétion
        for f in as_completed(futures):
            results.append(f.result())

    # Agrégation par mois : somme des stocks/flux sur tous les threads
    avg_history = pd.concat(results).groupby('month').sum().reset_index()
    return avg_history


# ------------------------------------------------------------------
# Fonction de calcul des statistiques agrégées issues de la simulation.
# Différence majeure avec la calibration : calcule plusieurs fenêtres
# temporelles pour pouvoir analyser à la fois l'état stationnaire pré-réforme
# (mois 228-240), le court terme post-réforme (mois 240-252) et la
# trajectoire mensuelle complète post-réforme (mois 240-480, pour visualiser
# la montée en charge).
# Exportation optionnelle vers deux fichiers CSV :
#   - file1 : agrégats moyens et différences pré/post (état stationnaire et court terme)
#   - file2 : trajectoire mensuelle complète post-réforme (montée en charge)
# ------------------------------------------------------------------
def moments_computation(avg_history, file1=None, file2=None):
    # Fenêtres temporelles utilisées pour l'analyse :
    avg_history_last_12 = avg_history[avg_history['month'] >= (months_sim - 12)].copy() # 12 derniers mois (nouvel état stationnaire après montée en charge complète)
    avg_history_steady_pre = avg_history[avg_history['month'].isin(range(240-12, 240))].copy() # 12 mois précédant la réforme (état stationnaire pré-réforme, sert de référence pour les différences)
    avg_history_12_post = avg_history[avg_history['month'].isin(range(240, 240+12))].copy() # 12 premiers mois post-réforme (effet de court terme, première année de montée en charge)
    months_following = avg_history[avg_history['month'].isin(range(240, 480))].copy() # ensemble de la période post-réforme (trajectoire complète pour visualiser la montée en charge)

    # stocks/flux (moyennes sur les 12 derniers mois = état stationnaire post-réforme)
    employed_avg = avg_history_last_12['employed'].mean()
    employed_eligible_avg = avg_history_last_12['employed_eligible'].mean()
    ui_unemployed_avg = avg_history_last_12['ui_unemployed'].mean()
    noui_unemployed_avg = avg_history_last_12['noui_unemployed'].mean()
    noui_eligible_avg = avg_history_last_12['noui_eligible'].mean()
    employed_6_8_history = avg_history_last_12['employed_6_8_history'].mean()
    flow_out_ui_avg = avg_history_last_12['flow_out_ui'].mean()
    flow_into_ui_avg = avg_history_last_12['flow_into_ui'].mean()
    flow_recharg_avg = avg_history_last_12['flow_recharg'].mean()
    flow_reprise_avg = avg_history_last_12['flow_reprise'].mean()
    odd_6_8_hist = avg_history_last_12['odd_6_8'].mean()
    odd_hist = avg_history_last_12['odd'].mean()
    epuise_uni_avg = avg_history_last_12['epuise_uni'].mean()
    # Part des chômeurs indemnisés sous régime post-réforme parmi l'ensemble des chômeurs indemnisés (mesure de l'avancement de la montée en charge ; tend vers 1 à long terme)
    share_regime_new = avg_history_last_12['ui_unemployed_regime_new'].mean()/max(avg_history_last_12['ui_unemployed'].mean(), 1e-12)

    # Liste ordonnée des moments retournés
    stocks = [
        employed_avg, ui_unemployed_avg, noui_unemployed_avg, noui_eligible_avg, employed_eligible_avg,
        employed_6_8_history, flow_out_ui_avg, flow_into_ui_avg, flow_recharg_avg, flow_reprise_avg,
        odd_6_8_hist, odd_hist, epuise_uni_avg, share_regime_new
    ]

    # Exports optionnels vers des fichiers csv externes
    if file1 is not None:
        # Trois lignes exportées (à l'échelle de l'économie réelle) :
        # - avg_last_12 : moyennes sur les 12 derniers mois (= post-réforme stabilisé)
        # - diff_pre_post_stable : différence entre l'état stationnaire post-réforme et l'état stationnaire pré-réforme (effet à long terme)
        # - diff_pre_post_12 : différence entre les 12 premiers mois post-réforme et l'état stationnaire pré-réforme (effet à court terme)
        s_last = (avg_history_last_12.mean(numeric_only=True) * economy_scale).rename('avg_last_12')
        s_diff_stable = ((avg_history_last_12.mean(numeric_only=True) - avg_history_steady_pre.mean(numeric_only=True)) * economy_scale).rename('diff_pre_post_stable')
        s_diff_12 = ((avg_history_12_post.mean(numeric_only=True) - avg_history_steady_pre.mean(numeric_only=True)) * economy_scale).rename('diff_pre_post_12')
        df_concat = pd.concat([s_last, s_diff_stable, s_diff_12], axis=1).T
        df_concat.to_csv(file1, index=False)

    if file2 is not None:
        # Exportation de la trajectoire mensuelle complète post-réforme (mois 240 à 480) à l'échelle de l'économie réelle
        # Utilisée par le notebook pour tracer les courbes de montée en charge et calculer les économies annuelles à différents horizons (1 an, 5 ans, 10 ans)
        scaled = months_following * economy_scale
        pre = avg_history_steady_pre.mean(numeric_only=True)*economy_scale
        diff = scaled.sub(pre, axis='columns')  # disponible si l'on souhaite exporter directement les différences pré/post (non utilisé ici, la différence étant calculée a posteriori dans le notebook)
        scaled.to_csv(file2, index=False)

    return stocks


# ------------------------------------------------------------------
# Fonction lançant la simulation pour un jeu de paramètres donné et exportant
# les statistiques agrégées.
# Différence majeure avec la calibration : cette fonction ne calcule pas
# d'écart à des cibles (pas de minimisation : les paramètres structurels sont
# déjà calibrés) ; elle se contente de lancer la simulation et d'écrire les
# fichiers de sortie. Le nom "objective_function" est conservé par cohérence
# avec le code de calibration.
# ------------------------------------------------------------------
def objective_function(params, b, e_jdr, e_jfr, file1=None, file2=None):
    # params: [jfr, jdr, jdr_coef, take_up_rate] - paramètres structurels calibrés sur le régime de référence
    jfr, jdr, jdr_coef, take_up_r = params

    # Conversion en types numériques et application des transformations
    job_finding_rate = float(np.exp(jfr)) # taux de retour à l'emploi
    job_destruction_rate = float(np.exp(jdr)) # delta 0 (partie constante du taux de destruction d'emploi)
    jdr_coef = float(jdr_coef) # delta 1 (sensibilité du taux de destruction d'emploi à l'ancienneté)
    take_up_rate = float(take_up_r) # taux de recours
    

    # Lancement de la simulation avec les paramètres structurels calibrés et les paramètres de réforme (b, e_jdr, e_jfr)
    avg_history = run_all(b= b, e_jdr=e_jdr, e_jfr =e_jfr, job_finding_rate=job_finding_rate, job_destruction_rate=job_destruction_rate, jdr_coef=jdr_coef, take_up_rate=take_up_rate)

    # Calcul et exportation des statistiques agrégées (aucune valeur retournée car pas de minimisation)
    moments_computation(avg_history, file1, file2)
    return 

# ------------------------------------------------------------------
# Fonction principale appelée par le notebook pour simuler la réforme
# "baisse uniforme du montant d'allocation".
# C'est cette fonction qui orchestre les 3 scénarios à simuler successivement :
#   1) contrefactuel : b_pre = b_post (pas de réforme), élasticités à 0
#      → trajectoire de référence sans réforme
#   2) effet mécanique : b_post < b_pre, élasticités à 0
#      → effet à comportements inchangés
#   3) effet total : b_post < b_pre, élasticités calibrées (e_jfr=-0,3, e_jdr=0,1)
#      → effet total intégrant les réponses comportementales
# L'effet comportemental s'obtient ensuite par différence : (3) - (2).
#
# Paramètres :
#   - params : paramètres structurels calibrés sur le régime de référence
#   - e_jfr  : élasticité du taux de retour à l'emploi au montant d'allocation
#              (valeur calibrée -0,3 pour les réformes du montant ; mettre à 0
#              pour simuler l'effet mécanique seul)
#   - e_jdr  : élasticité du taux de séparation au montant d'allocation
#              (valeur calibrée 0,1 pour les réformes du montant ; mettre à 0
#              pour simuler l'effet mécanique seul)
#   - b_pre  : montant d'allocation pré-réforme (par défaut : 1212.9€,
#              montant moyen d'allocation des individus affectés par la
#              réforme baisse uniforme - calibration MiDAS)
#   - b_post : montant d'allocation post-réforme (par défaut : 1200.984€,
#              correspondant à une baisse de ~1% appliquée à b_pre, ajustée
#              pour tenir compte de la dégressivité - calibration MiDAS)
#   - file1, file2 : chemins des fichiers csv de sortie (cf. moments_computation)
# ------------------------------------------------------------------
def f_niveau_uniforme(params, e_jfr, e_jdr, b_post=1200.984, b_pre=1212.9,  file1 = None, file2=None): # attention à la calibration de b_pre et b_post : la dégressivité est intégrée dans les valeurs (calibration MiDAS)
    
    # Affichage du type de simulation lancée (diagnostic)
    if b_pre == b_post:
        print("simulation contrefactuelle")
    else:
        if e_jfr== 0:
            if e_jdr == 0:
                print("simulation après application de l'effet mecanique")
            else:
                print("simulation après application de l'effet séparation (+ mécanique)")
        else:
            if e_jdr == 0:
                print("simulation après application de l'effet sortie (+ mécanique) sans application de l'effet séparation")
            else:
                print("simulation après application de l'effet total (mécanique +séparation+sortie)")
    e_jfr = e_jfr # -0.3 # passer à 0 pour avoir l'effet mécanique seul
    e_jdr = e_jdr # 0.1 # passer à 0 pour avoir l'effet mécanique seul

    # Construction du vecteur b = [pré-réforme, post-réforme] propagé à toute la simulation
    b=[b_pre, b_post]
   
    return objective_function(params = params, b = b, e_jdr = e_jdr, e_jfr = e_jfr, file1 = file1, file2=file2)

