# ============================================================
# DESCRIPTION DU FICHIER
# Ce fichier permet de simuler les effets de la réforme "raccourcissement de
# la période de référence d'éligibilité pour les seniors" (par exemple
# alignement de la PRA d'éligibilité des 55-57 ans sur celle du régime
# général, soit un passage de 36 mois à 24 mois) à partir des paramètres
# calibrés sur le régime général.
# (pas de calibration spécifique ici pour les séniors mais cela serait souhaitable).
#
# Différence majeure avec les réformes du montant et de la durée :
# Pour les réformes d'éligibilité, le code ne simule QUE L'EFFET MÉCANIQUE.
# Aucune élasticité n'est mobilisée : les fonctions update_status,
# objective_function et f_PRA_elig_senior ne reçoivent pas d'arguments e_jfr
# ou e_jdr, et les taux de transition restent au niveau calibré sur le
# régime de référence senior quelle que soit la situation de l'individu.
#
# Spécificité-clé 1 : PRA_elig vs pra_max
# Le paramètre qui varie est PRA_elig, qui désigne la période de référence
# utilisée pour le TEST D'ÉLIGIBILITÉ. Il est conceptuellement distinct de
# pra_max, qui désigne la période de référence utilisée pour le CALCUL DE LA
# DURÉE DES DROITS OUVERTS (qui reste inchangée par cette réforme).
# - PRA_elig[self.ref] : utilisé en dynamique d'emploi pour calculer
#   tot_months_worked (puis comparé à DA_min pour l'éligibilité) et lors
#   du test de rechargement
# - pra_max[self.ref] : utilisé partout ailleurs (dynamique de chômage,
#   calcul des droits ouverts, etc.)
#
# Spécificités-clé 2 : population senior
#   - La simulation ne porte que sur la population senior 55-57 ans :
#     l'ensemble des cibles de calage et la population totale sont multipliés
#     par share_seniors = 2,65% (part des 55-57 ans dans la population des
#     demandeurs d'emploi indemnisables, calibration MiDAS)
#   - La période de référence pour le calcul des droits (pra_max) est de
#     36 mois (régime senior, vs 24 mois pour le régime général)
#   - Le plafond de durée est de 22 mois (régime senior)
#   - economy_scale est fixé à 1 (et non 10) en raison de la petite taille
#     de la population senior, pour conserver suffisamment d'individus dans
#     la simulation
#
# Deux composantes peuvent être simulées via la fonction f_PRA_elig_senior :
#   - contrefactuel (pas de réforme) : PRA_elig_pre = PRA_elig_post
#   - mécanique : PRA_elig_post < PRA_elig_pre (raccourcissement, comportements inchangés)
# L'effet total = effet mécanique pour ce code (pas de réponse comportementale).
# La réforme s'applique au flux à partir du mois 240 et ne modifie pas les
# droits déjà ouverts.
#
# Utilisé par : simulations_reformes.ipynb (partie Réformes seniors)
# 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_reformes.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 senior frontalière affectée).
# Les statistiques cibles sont également mises à l'échelle de la population
# senior 55-57 ans via le coefficient share_seniors.
# ============================================================

# Paramètres financiers (calibration MiDAS sur la population senior 55-57 ans affectée)
b = [935.7, 935.7]  # niveau d'allocation moyen mensuel des individus senior affectés par la réforme ; structure liste conservée pour homogénéité, le montant ne varie pas pour une réforme d'éligibilité. 
salaire_pre = 1389 # salaire mensuel moyen avant entrée au chômage pour la population senior affectée (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 = 122 # niveau moyen des minima sociaux pour la population frontalière (122€ vs 165,86€ utilisé dans les réformes de montant/durée) ; la population à la frontière d'éligibilité a une calibration différente car la population affectée n'est pas éligible à l'ASS.
rsa_eligibility_rate = 1 # fixé à 1 car le taux de recours est déjà pris en compte dans b_rsa

# Paramètres de l'assurance chômage - régime senior 55-57 ans
# 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.
# Le paramètre PRA_elig (distinct de pra_max) n'apparaît pas ici car c'est lui
# qui varie pré/post pour cette réforme (passé en argument à f_PRA_elig_senior).
pra_max = [36, 36] # période de référence d'affiliation senior = 36 mois (vs 24 mois pour le régime général) ; utilisée pour le CALCUL des droits ouverts, inchangée par cette réforme
DA_min = [6, 6] # seuil de durée minimale d'affiliation (inchangé par cette 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 (inchangé)
pbd_max = [22, 22] # durée potentielle d'indemnisation maximale senior = 22 mois (vs 18 mois pour le régime général), inchangée par cette réforme

# Paramètres de la simulation
nbr_threads = 10 # nombre de coeurs utilisés dans la parallélisation
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 = 1 # 1 individu dans la simulation représente economy_scale individus dans l'économie réelle, mis à 1 ici (vs 10 dans les autres codes) en raison de la petite population senior, pour garder suffisamment d'individus dans la simulation 

# Statistiques cibles, mises à l'échelle de la population senior par share_seniors
# (voir parametres_regime_reference.py pour le détail des valeurs absolues utilisées)
share_seniors = 0.0265 # part des 55-57 ans dans la population des demandeurs d'emploi indemnisables (calibration MiDAS) - multiplie l'ensemble des cibles ci-dessous pour redimensionner les statistiques à la population senior
nbr_individuals_emp = 19341366/economy_scale # nombre d'individus en emploi (population totale, le scaling senior est appliqué dans tot_pop_guess)
target_UI_stock = (2428884)/economy_scale # nombre d'individus au chômage indemnisé
tot_pop_guess = share_seniors*(nbr_individuals_emp  + target_UI_stock+ (3279392)/economy_scale) # Population totale senior : (emploi + UI + non-UI) × share_seniors
target_reprise_inflow = share_seniors*(1345603/12)/economy_scale # Nombre de reprises de droits chez les seniors
target_recharg_inflow = share_seniors*(336752/12)/economy_scale # nombre de rechargements chez les seniors
target_UI_inflow = (share_seniors*(1895926)/12)/economy_scale # ouvertures de droits initiales et autres chez les seniors
target_UI_outflow = share_seniors*(2279654/12)/economy_scale # sorties courtes + longues chez les seniors
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 = share_seniors*1/12* 572727/economy_scale # nombre d'individus sortant/entrant sur le marché du travail chez les seniors (calibré pour un renouvellement complet en 44 ans, 62-18, à l'échelle senior)
flow_system =  flow_system_tot/tot_pop_guess # taux de renouvellement de la population senior

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


# ============================================================
# CLASSE INDIVIDU
# Représente un agent senior (55-57 ans) 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 PRA_elig = [PRA_elig_pre, PRA_elig_post]
#     qui désigne la période de référence utilisée UNIQUEMENT pour le test
#     d'éligibilité (et le test de rechargement). Le calcul de la durée
#     des droits ouverts reste fondé sur pra_max global (36 mois senior),
#     qui ne change pas
#   - possède un indicateur ref (0 = régime pré-réforme, 1 = régime
#     post-réforme) qui détermine quelle PRA_elig appliquer
#   - possède un compteur month_simu pour repérer le mois courant
#     de la simulation et déclencher la réforme au mois 240
# NB : à la différence des réformes du montant ou de la durée, AUCUN effet comportemental
# n'est calculé ici (pas de paramètres e_jdr/e_jfr) ; les taux
# de transition restent à leur valeur calibrée sur le régime de référence senior
# ============================================================
class Individual:
    def __init__(self, id, PRA_elig, 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 PRA senior (36 mois)
        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, 1 = régime post-réforme

        # 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 (population déjà restreinte aux seniors)
            '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 senior (plus faible que pour la population générale, ajusté dans le notebook)
        self.job_destruction_rate = job_destruction_rate  # Taux mensuel de destruction d'emploi (non modulé par la réforme dans ce code)

    # ------------------------------------------------------------------
    # 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 PRA_elig (vecteur [PRA_elig_pre, PRA_elig_post])
    #   - utilise PRA_elig[self.ref] pour calculer tot_months_worked dans la
    #     dynamique d'emploi (test d'éligibilité) et dans le test de rechargement
    #   - utilise pra_max[self.ref] (global = 36, senior) pour le calcul de la
    #     durée des droits ouverts (inchangée par la réforme)
    #   - passe l'individu en ref=1 à partir du mois 240 dès qu'il a consommé
    #     ses éventuels droits résiduels (réforme appliquée au flux et non au stock)
    # ------------------------------------------------------------------
    def update_status(self, PRA_elig):
        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
        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 PRA senior "calcul des droits" (pra_max = 36)
        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é"
        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é
            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
                if self.month_simu >= 240:
                    self.ref = 1
            # Mise à jour de l'historique : marque le mois courant comme travaillé
            self.months_worked = [1] + self.months_worked[:-1]
            # CHANGEMENT-CLÉ : tot_months_worked est ici recalculé sur PRA_elig (test d'éligibilité), et non sur pra_max
            self.tot_months_worked = sum(self.months_worked[0:PRA_elig[self.ref]])
            # 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:
                # L'individu a encore des droits résiduels : éligible (la PRA_elig ne s'applique pas aux droits déjà ouverts)
                self.eligible = True
                # Taux de destruction d'emploi non modulé par la réforme (pas d'élasticité dans ce code)
                jdrate = self.job_destruction_rate*np.exp(self.jdr_coef*months_worked_continu)     
            else:
                # Test d'éligibilité : nombre de mois travaillés sur PRA_elig comparé à DA_min (global, inchangé pour cette réforme)
                if self.tot_months_worked >= DA_min[self.ref]:
                    self.eligible = True
                    jdrate = self.job_destruction_rate*np.exp(self.jdr_coef*months_worked_continu) 
                else:
                    self.eligible = False
                    jdrate = self.job_destruction_rate*np.exp(self.jdr_coef*months_worked_continu)

            # 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
                        # IMPORTANT : on RECALCULE tot_months_worked sur pra_max (36, senior) pour le calcul des droits ouverts ; la durée des droits est donc fondée sur la PRA "calcul des droits" senior, inchangée par la réforme
                        self.tot_months_worked = sum(self.months_worked[0:pra_max[self.ref]])
                        self.months_odd = self.tot_months_worked
                        months_worked_act = self.months_worked[0:pra_max[self.ref]]
                        # Longueur de la PRC (durée entre le premier et dernier mois travaillé sur la période de référence pra_max senior)
                        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 senior (22 mois)
                        rights = min(pbd_max[self.ref], max(round(min(self.tot_months_worked*1.7, length_prc)*coef_conversion_UI[self.ref]),6))
                        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 pra_max senior
            # 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 non modulé par la réforme (pas d'élasticité dans ce code)
                jfr_ui = self.job_finding_rate
                    
                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 :
                        # NB : pour les seniors, le seuil sur l'index du dernier emploi est 36 mois (vs 48 dans le régime général), pour des contraintes de durée de mémoire
                        if sum(self.months_worked) == 0 or self.months_worked.index(1) >= 36: # si pas retravaillé 
                            self.tot_months_worked_recharg = 0
                        else:
                            # CHANGEMENT-CLÉ : le calcul des mois utilisables pour le rechargement utilise PRA_elig (et non pra_max) car le rechargement est aussi un test d'éligibilité
                            self.tot_months_worked_recharg = sum(self.months_worked[self.months_worked.index(1):self.months_worked.index(1)+PRA_elig[self.ref]])
                        if self.tot_months_worked_recharg >= DA_min[self.ref]: # si a retravaillé au moins DA_min mois sur PRA_elig => ouvre un nouveau droit
                            # Rechargement des droits au chômage
                            self.flow_recharg = 1
                            self.odd = 1
                            self.months_odd = self.tot_months_worked
                            # Pour le CALCUL de la durée des droits rechargés, on revient à pra_max senior
                            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
                            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.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é)
                    self.status = 'employed'
                    self.take_up = 0
                # sinon pas de changement d'état


# ===================================================================
# CLASSE SIMULATION
# Orchestre la simulation sur un ensemble d'individus seniors et un nombre
# de mois donné. Agrège les statistiques mensuelles dans un DataFrame.
# Différences avec la calibration :
#   - reçoit en argument PRA_elig (vecteur [PRA_elig_pre, PRA_elig_post]),
#     propagé 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 deux régimes
#     (ici b[0]=b[1] donc le coût ne dépend pas de la répartition ; la
#     structure est conservée pour homogénéité avec les réformes du montant)
# ===================================================================
  

class Simulation:
    # ------------------------------------------------------------------
    # Fonction initialisant les individus
    # ------------------------------------------------------------------
    def __init__(self, PRA_elig, 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
            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,
                PRA_elig = PRA_elig,
                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, PRA_elig):
        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 seniors sous régime post-réforme (statistique spécifique aux simulations de réforme : permet 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 seniors 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é sur leur PRA_elig courante (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 : b[0] = b[1] pour cette réforme (le montant ne change pas), donc le coût total ne dépend pas de la répartition régime pré/post ; la formule est conservée pour homogénéité avec les codes des réformes du montant
            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 PRA_elig à chaque update individuel)
            for ind in self.individuals:
                ind.update_status(PRA_elig)
        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 de seniors donnée
# Différence avec la calibration : propage en plus le paramètre de réforme
# PRA_elig à la classe Simulation
# -------------------------------------------------------------------
def run_one(sim_id, PRA_elig, 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(PRA_elig = PRA_elig, 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(PRA_elig)
    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 la population senior.
# Différence avec la calibration : propage en plus le paramètre de réforme
# PRA_elig à chaque thread
# -------------------------------------------------------------------
def run_all(PRA_elig, 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 senior 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, PRA_elig, 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 ; ici economy_scale=1) :
        # - 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.
# Différence avec les réformes du montant/durée : ne reçoit PAS d'élasticités
# e_jdr/e_jfr (effet mécanique uniquement)
# ------------------------------------------------------------------
def objective_function(params, PRA_elig, file1=None, file2=None):
    # params: [jfr, jdr, jdr_coef, take_up_rate] - paramètres structurels calibrés sur la population senior (notamment job_finding_rate plus faible que pour la population générale, ajusté par le notebook)
    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 senior
    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 le paramètre de réforme PRA_elig
    avg_history = run_all(PRA_elig = PRA_elig, 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
# "raccourcissement de la période de référence d'éligibilité pour les
# seniors" (par exemple alignement de la PRA d'éligibilité des 55-57 ans
# sur celle du régime général, soit un passage typique de 36 à 24 mois).
# C'est cette fonction qui orchestre les 2 scénarios à simuler successivement :
#   1) contrefactuel : PRA_elig_pre = PRA_elig_post (pas de réforme)
#      → trajectoire de référence sans réforme
#   2) effet mécanique : PRA_elig_post < PRA_elig_pre (raccourcissement)
#      → effet à comportements inchangés (seul effet disponible, pas d'élasticité)
#
# NB : contrairement aux réformes du montant et de la durée des droits,
# il n'y a pas de scénario "effet total" intégrant des réponses
# comportementales.
#
# Paramètres :
#   - params : paramètres structurels calibrés sur la population senior
#   - PRA_elig_pre  : période de référence d'éligibilité senior pré-réforme
#                     (par défaut : 36 mois, régime senior actuel)
#   - PRA_elig_post : période de référence d'éligibilité senior post-réforme
#                     (typiquement 24 mois, alignement sur le régime général)
#   - file1, file2 : chemins des fichiers csv de sortie (cf. moments_computation)
# ------------------------------------------------------------------
def f_PRA_elig_senior(params, PRA_elig_post, PRA_elig_pre=36,  file1 = None, file2=None):
    # Affichage du type de simulation lancée (diagnostic)
    if PRA_elig_pre == PRA_elig_post:
        print("simulation contrefactuelle")
    else:
        print("simulation après application de l'effet mecanique")

    # Construction du vecteur PRA_elig = [pré-réforme, post-réforme] propagé à toute la simulation
    PRA_elig = [PRA_elig_pre, PRA_elig_post] 
   
    return objective_function(params = params, PRA_elig = PRA_elig, file1 = file1, file2=file2)




