#!/usr/bin/env python3
"""Chapitre 37 — Apprentissage automatique et explicabilité : exemple intégral sur les 600 ventes du cas 1
(Gatineau, DONNÉES SYNTHÉTIQUES, graine 20261031).

Le modèle hédonique de référence est `traite.cas1.hedonique` (semi-log, MCO, HC1) ; il n'est PAS redéfini ici. Ce script lui
oppose un arbre de régression, une forêt aléatoire (Breiman 2001), un gradient boosting (Friedman 2001 ; scikit-learn) et XGBoost
(Chen et Guestrin 2016), sous trois schémas de validation : plis aléatoires (les MÊMES dix plis que `ch18_hedonique.validation_croisee`,
graine 20261031), blocs spatiaux stratifiés par secteur (k-moyennes sur les coordonnées), et partition temporelle (ventes de 2026
prédites à partir de 2023-2025). Il ajoute la validation « un secteur exclu » (extrapolation hors support), un scénario où la loi de
génération est perturbée par une interaction et une rupture déclarées (le boosting y gagne), les valeurs SHAP (TreeSHAP) du boosting,
la dépendance partielle, l'importance par permutation, la prédiction du sujet par chaque modèle et les chiffres des exercices.
Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch37_ml.py"""
from __future__ import annotations
import sys, math, pathlib, warnings
import numpy as np, pandas as pd
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
from traite import cas1
from traite import hedonique as h

warnings.filterwarnings("ignore", category=FutureWarning); warnings.filterwarnings("ignore", category=UserWarning)
ROOT = pathlib.Path(__file__).resolve().parents[3]
GRAINE = 20261031
N_ARBRES_FORET = 500

# ----------------------------------------------------------------------------- paramètres de génération (README synthétique)
_g = (1 + 9.06 / 100) ** 0.5 - 1
VRAIS = {"ln_S": 0.55, "ln_T": 0.10, "age": -0.0055, "age2": 0.00003, "qualite": 0.06, "garage": 0.035, "sous_sol": 0.03, "piscine": 0.02,
         "etages": 0.02, "dist_centre": -0.006, "sigma": 0.085, "mois": math.log1p(_g) / 12,
         "sect": {"Hull": 0.0, "Gatineau": math.log(464_900 / 475_000), "Aylmer": math.log(530_000 / 475_000),
                  "Buckingham": math.log(395_000 / 475_000), "Masson-Angers": math.log(395_000 / 475_000) - 0.02}}
# scénario non linéaire déclaré (section 37.6) : prime « patrimoine rénové » et « mur de la mi-vie »
SCENARIO = {"prime_qualite_age": 0.12, "qualite_min": 4, "age_min_prime": 40, "mur": -0.10, "age_mur": (25, 45)}

# ----------------------------------------------------------------------------- régresseurs
VARIABLES_ML = ["superficie", "terrain", "age", "qualite", "garage", "sous_sol", "piscine", "etages", "dist_centre", "lat", "lon", "mois"] + \
               [f"sect_{s}" for s in cas1.SECTEURS]
ETIQUETTES = {"superficie": "Superficie habitable (m²)", "terrain": "Terrain (m²)", "age": "Âge (années)", "qualite": "Qualité (1-5)",
              "garage": "Garage (places)", "sous_sol": "Sous-sol fini", "piscine": "Piscine", "etages": "À étages", "dist_centre": "Distance au centre (km)",
              "lat": "Latitude", "lon": "Longitude", "mois": "Mois depuis 2024-07", "sect_Hull": "Secteur Hull", "sect_Gatineau": "Secteur Gatineau",
              "sect_Aylmer": "Secteur Aylmer", "sect_Buckingham": "Secteur Buckingham", "sect_Masson-Angers": "Secteur Masson-Angers"}

def matrice_ml(df: pd.DataFrame, s: dict | None = None) -> pd.DataFrame:
    """Régresseurs bruts pour les méthodes à base d'arbres : aucune transformation (les arbres sont invariants aux transformations
    monotones), secteur en indicatrices complètes, coordonnées brutes, temps en mois depuis la date de référence du jeu."""
    X = pd.DataFrame(index=df.index)
    X["superficie"] = df["superficie_habitable_m2"].astype(float); X["terrain"] = df["superficie_terrain_m2"].astype(float)
    X["age"] = df["age"].astype(float); X["qualite"] = df["qualite_1a5"].astype(float); X["garage"] = df["garage_places"].astype(float)
    X["sous_sol"] = df["sous_sol_fini"].astype(float); X["piscine"] = df["piscine"].astype(float)
    X["etages"] = (df["type"] == "à étages").astype(float); X["dist_centre"] = df["dist_centre_km"].astype(float)
    X["lat"] = df["lat"].astype(float); X["lon"] = df["lon"].astype(float)
    d = df["date_vente"] if "date_vente" in df else pd.Series([pd.Timestamp(s["date_evaluation"])] * len(df), index=df.index)
    X["mois"] = ((pd.to_datetime(d) - cas1.REF_DATE).dt.days / 30.4375).astype(float)
    for sct in cas1.SECTEURS:
        X[f"sect_{sct}"] = (df["secteur"] == sct).astype(float)
    return X[VARIABLES_ML]

# ----------------------------------------------------------------------------- modèles
class Hedonique:
    """Enveloppe du modèle hédonique de référence (traite.hedonique.mco sur traite.cas1.matrice) avec l'interface fit/predict."""
    def __init__(self, colonnes_exclues=()): self.excl = list(colonnes_exclues)
    def fit(self, X: pd.DataFrame, y):
        self.cols = [c for c in X.columns if c not in self.excl]
        self.r = h.mco(X[self.cols].values, np.asarray(y, float)); self.beta = self.r["beta"]; return self
    def predict(self, X: pd.DataFrame): return X[self.cols].values @ self.beta

def _foret(**k):
    from sklearn.ensemble import RandomForestRegressor
    return RandomForestRegressor(**{**dict(n_estimators=N_ARBRES_FORET, min_samples_leaf=2, max_features=0.5, random_state=GRAINE, n_jobs=-1), **k})

def _boosting(**k):
    from sklearn.ensemble import GradientBoostingRegressor
    return GradientBoostingRegressor(**{**dict(n_estimators=1200, learning_rate=0.03, max_depth=2, subsample=0.8, min_samples_leaf=5, random_state=GRAINE), **k})

def _xgboost(**k):
    import xgboost as xgb
    return xgb.XGBRegressor(**{**dict(n_estimators=1200, learning_rate=0.03, max_depth=2, subsample=0.8, colsample_bytree=0.8, min_child_weight=5,
                                      reg_lambda=1.0, random_state=GRAINE, n_jobs=4, verbosity=0), **k})

def _arbre(**k):
    from sklearn.tree import DecisionTreeRegressor
    return DecisionTreeRegressor(**{**dict(max_depth=3, min_samples_leaf=5, random_state=GRAINE), **k})

HYPERPARAMETRES = {"Arbre (profondeur 3)": "profondeur 3, feuilles ≥ 5", "Forêt aléatoire": f"{N_ARBRES_FORET} arbres, feuilles ≥ 2, 50 % des variables par division, bootstrap",
                   "Gradient boosting": "1 200 arbres de profondeur 2, ν = 0,03, sous-échantillon 80 %, feuilles ≥ 5",
                   "XGBoost": "1 200 arbres de profondeur 2, ν = 0,03, sous-échantillon 80 %, 80 % des colonnes, poids min. 5, λ = 1"}
MODELES = {"Hédonique (référence)": ("hedo", Hedonique), "Arbre (profondeur 3)": ("ml", _arbre), "Forêt aléatoire": ("ml", _foret),
           "Gradient boosting": ("ml", _boosting), "XGBoost": ("ml", _xgboost)}

# ----------------------------------------------------------------------------- plis
def plis_aleatoires(n: int, nfolds: int = 10, graine: int = GRAINE) -> list[np.ndarray]:
    """Les mêmes dix plis que ch18_hedonique.validation_croisee (même graine, même permutation, même découpage)."""
    rng = np.random.default_rng(graine); idx = rng.permutation(n); return list(np.array_split(idx, nfolds))

def plis_spatiaux(v: pd.DataFrame, par_secteur: int = 2, graine: int = GRAINE) -> tuple[list[np.ndarray], np.ndarray]:
    """Blocs spatiaux stratifiés : dans chaque secteur, k-moyennes (k = 2) sur les coordonnées en km → 10 blocs contigus.
    Chaque secteur reste représenté à l'estimation (par l'autre moitié), mais le voisinage immédiat d'une vente testée ne l'est jamais."""
    from sklearn.cluster import KMeans
    km = h.coordonnees_km(v["lat"].to_numpy(), v["lon"].to_numpy()); bloc = np.zeros(len(v), int); b = 0
    for sct in cas1.SECTEURS:
        ii = np.where(v["secteur"].to_numpy() == sct)[0]
        lab = KMeans(par_secteur, n_init=10, random_state=graine).fit_predict(km[ii]); bloc[ii] = b + lab; b += par_secteur
    return [np.where(bloc == j)[0] for j in range(b)], bloc

def plis_secteurs(v: pd.DataFrame) -> list[np.ndarray]:
    """Un secteur exclu à la fois (extrapolation hors support spatial)."""
    return [np.where(v["secteur"].to_numpy() == sct)[0] for sct in cas1.SECTEURS]

def partition_temporelle(v: pd.DataFrame, date: str = "2026-01-01") -> tuple[np.ndarray, np.ndarray]:
    test = np.where(v["date_vente"] >= pd.Timestamp(date))[0]; return np.setdiff1d(np.arange(len(v)), test), test

# ----------------------------------------------------------------------------- métriques
def metriques(pred_log: np.ndarray, y: np.ndarray, facteur: float = 1.0) -> dict:
    """RMSE du logarithme, R² hors échantillon, biais moyen (log), RMSE en dollars, MdAPE, statistiques de ratio IAAO
    (traite.hedonique.statistiques_ratio). `facteur` : correction de retransformation (Duan) appliquée aux prédictions en niveau."""
    m = ~np.isnan(pred_log); e = pred_log[m] - y[m]; P = np.exp(y[m]); V = np.exp(pred_log[m]) * facteur
    r = h.statistiques_ratio(V, P)
    return {"n": int(m.sum()), "rmse_log": float(np.sqrt(np.mean(e ** 2))), "biais_log": float(e.mean()),
            "r2_hors": float(1 - np.sum(e ** 2) / np.sum((y[m] - y[m].mean()) ** 2)), "rmse_dollars": float(np.sqrt(np.mean((V - P) ** 2))),
            "mdape": float(np.median(np.abs(V / P - 1)) * 100), "COD": r["COD"], "PRD": r["PRD"], "PRB": r["PRB"], "PRB_se": r["PRB_se"],
            "ratio_median": r["ratio_median"]}

def valider(constructeur, genre: str, Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray, plis: list[np.ndarray]) -> np.ndarray:
    """Prédictions hors échantillon (chaque pli prédit par un modèle estimé sans lui)."""
    X = Xh if genre == "hedo" else Xm; n = len(y); pred = np.full(n, np.nan); tous = np.arange(n)
    for f in plis:
        tr = np.setdiff1d(tous, f); m = constructeur(); m.fit(X.iloc[tr], y[tr]); pred[f] = np.asarray(m.predict(X.iloc[f]), float)
    return pred

def facteur_duan(pred_log: np.ndarray, y: np.ndarray) -> float:
    """Facteur de retransformation estimé sur les résidus hors échantillon (honnête pour les modèles qui surajustent en échantillon)."""
    m = ~np.isnan(pred_log); return float(np.mean(np.exp(y[m] - pred_log[m])))

def comparaison(v: pd.DataFrame, Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray) -> dict:
    """Trois schémas de validation pour les cinq modèles. Retourne les prédictions hors échantillon et les métriques."""
    n = len(y); pa = plis_aleatoires(n); ps, bloc = plis_spatiaux(v); tr_t, te_t = partition_temporelle(v)
    out = {"plis": {"aleatoire": pa, "spatial": ps, "bloc_spatial": bloc, "temporel": (tr_t, te_t)}, "pred": {}, "metriques": {}, "duan": {}}
    for nom, (genre, ctor) in MODELES.items():
        p_a = valider(ctor, genre, Xh, Xm, y, pa); F = facteur_duan(p_a, y)
        p_s = valider(ctor, genre, Xh, Xm, y, ps)
        p_t = valider(ctor, genre, Xh, Xm, y, [te_t])            # un seul pli de test : les ventes de 2026
        out["pred"][nom] = {"aleatoire": p_a, "spatial": p_s, "temporel": p_t}; out["duan"][nom] = F
        out["metriques"][nom] = {"aleatoire": metriques(p_a, y, F), "spatial": metriques(p_s, y, F), "temporel": metriques(p_t, y, F)}
    # biais temporel par mois de 2026 pour le boosting et l'hédonique
    mois = Xm["mois"].to_numpy()
    out["biais_mensuel_2026"] = {nom: pd.Series(out["pred"][nom]["temporel"][te_t] - y[te_t]).groupby(np.floor(mois[te_t])).mean().to_dict()
                                 for nom in ("Hédonique (référence)", "Gradient boosting", "Forêt aléatoire")}
    return out

def boosting_detendance(v: pd.DataFrame, Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray) -> dict:
    """Remède à l'extrapolation temporelle : boosting estimé sur ln P − β̂_mois·mois (tendance de l'hédonique estimée sur l'entraînement),
    sans la variable mois, puis tendance rajoutée. Variante « feuilles ≥ 1 » pour montrer l'effet de la feuille extrême."""
    tr, te = partition_temporelle(v); n = len(y); cols = list(Xh.columns)
    hd = Hedonique().fit(Xh.iloc[tr], y[tr]); bm = hd.beta[cols.index("mois")]; mois = Xm["mois"].to_numpy()
    Xs = Xm.drop(columns=["mois"])
    g = _boosting().fit(Xs.iloc[tr], y[tr] - bm * mois[tr]); p = np.full(n, np.nan); p[te] = g.predict(Xs.iloc[te]) + bm * mois[te]
    g1 = _boosting(min_samples_leaf=1).fit(Xm.iloc[tr], y[tr]); p1 = np.full(n, np.nan); p1[te] = g1.predict(Xm.iloc[te])
    return {"beta_mois_entrainement": float(bm), "pred_detendance": p, "metriques_detendance": metriques(p, y), "pred_feuille1": p1,
            "metriques_feuille1": metriques(p1, y), "n_entrainement_apres_mois_17": int((mois[tr] > 17).sum())}

def validation_secteurs(v: pd.DataFrame, Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray) -> dict:
    """Un secteur exclu à la fois. L'hédonique n'a aucune information sur un secteur absent : son indicatrice est retirée et le secteur
    est prédit au niveau du secteur de référence (Hull) ; si Hull est exclu, Gatineau devient la référence. Ce choix est un JUGEMENT
    de l'analyste, explicite ; les arbres font implicitement le leur à partir des coordonnées et de la distance au centre."""
    n = len(y); plis = plis_secteurs(v); out = {"pred": {}, "metriques": {}, "biais_secteur": {}}
    for nom, (genre, ctor) in MODELES.items():
        pred = np.full(n, np.nan)
        for sct, f in zip(cas1.SECTEURS, plis):
            tr = np.setdiff1d(np.arange(n), f)
            if genre == "hedo":
                drop = f"sect_{sct}" if sct != "Hull" else "sect_Gatineau"; m = Hedonique([drop]).fit(Xh.iloc[tr], y[tr]); pred[f] = m.predict(Xh.iloc[f])
            else:
                m = ctor(); m.fit(Xm.iloc[tr], y[tr]); pred[f] = m.predict(Xm.iloc[f])
        out["pred"][nom] = pred; out["metriques"][nom] = metriques(pred, y)
        out["biais_secteur"][nom] = {sct: float(np.mean(pred[f] - y[f])) for sct, f in zip(cas1.SECTEURS, plis)}
    return out

# ----------------------------------------------------------------------------- scénario non linéaire déclaré
def perturbation_scenario(v: pd.DataFrame) -> np.ndarray:
    """ln P† = ln P + 0,12·1[qualité ≥ 4 et âge ≥ 40] − 0,10·1[25 ≤ âge ≤ 45] (hypothèse pédagogique déclarée)."""
    q = v["qualite_1a5"].to_numpy(float); a = v["age"].to_numpy(float); S = SCENARIO
    return S["prime_qualite_age"] * ((q >= S["qualite_min"]) & (a >= S["age_min_prime"])) + S["mur"] * ((a >= S["age_mur"][0]) & (a <= S["age_mur"][1]))

def scenario_nonlineaire(v: pd.DataFrame, Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray) -> dict:
    """Mêmes plis aléatoires ; hédonique de référence (mal spécifié désormais), forêt, boosting, XGBoost, et hédonique respécifié
    avec les deux indicatrices (ce que l'analyste ferait après avoir lu la dépendance partielle de l'âge et les interactions SHAP)."""
    delta = perturbation_scenario(v); y2 = y + delta; pa = plis_aleatoires(len(y))
    q = v["qualite_1a5"].to_numpy(float); a = v["age"].to_numpy(float); S = SCENARIO
    Xh2 = Xh.copy(); Xh2["prime"] = ((q >= S["qualite_min"]) & (a >= S["age_min_prime"])).astype(float); Xh2["mur"] = ((a >= S["age_mur"][0]) & (a <= S["age_mur"][1])).astype(float)
    out = {"delta": delta, "y": y2, "n_prime": int(Xh2["prime"].sum()), "n_mur": int(Xh2["mur"].sum()), "pred": {}, "metriques": {}, "Xh2": Xh2}
    for nom, (genre, ctor) in MODELES.items():
        if nom == "Arbre (profondeur 3)": continue
        p = valider(ctor, genre, Xh, Xm, y2, pa); out["pred"][nom] = p; out["metriques"][nom] = metriques(p, y2, facteur_duan(p, y2))
    p = valider(Hedonique, "hedo", Xh2, Xm, y2, pa); out["pred"]["Hédonique respécifié"] = p; out["metriques"]["Hédonique respécifié"] = metriques(p, y2, facteur_duan(p, y2))
    r = Hedonique().fit(Xh2, y2); cols = list(Xh2.columns)
    out["coef_respecifie"] = {"prime": float(r.beta[cols.index("prime")]), "mur": float(r.beta[cols.index("mur")]), "se_prime": float(r.r["se_hc1"][cols.index("prime")]),
                              "se_mur": float(r.r["se_hc1"][cols.index("mur")]), "sigma": float(np.sqrt(r.r["sigma2"])), "r2": float(r.r["r2"])}
    out["gbm_scenario"] = _boosting().fit(Xm, y2)
    out["hedo_scenario"] = Hedonique().fit(Xh, y2)
    return out

# ----------------------------------------------------------------------------- explicabilité
def shap_boosting(gbm, Xm: pd.DataFrame, xs: pd.DataFrame) -> dict:
    """TreeSHAP (Lundberg et al. 2020) du gradient boosting : valeurs SHAP des 600 ventes, du sujet, valeur de base, importance
    moyenne |φ|, interactions SHAP moyennes. Vérifie l'axiome d'efficacité (somme des φ + base = prédiction)."""
    import shap
    ex = shap.TreeExplainer(gbm); sv = np.asarray(ex.shap_values(Xm)); base = float(np.ravel(ex.expected_value)[0])
    svs = np.asarray(ex.shap_values(xs))[0]
    si = np.asarray(ex.shap_interaction_values(Xm)); inter = np.abs(si).mean(0).copy(); np.fill_diagonal(inter, 0.0)
    inter_df = pd.DataFrame(inter, index=Xm.columns, columns=Xm.columns)
    paires = inter_df.where(np.triu(np.ones_like(inter, bool), 1)).stack().sort_values(ascending=False)
    return {"valeurs": sv, "base": base, "sujet": pd.Series(svs, index=Xm.columns), "importance": pd.Series(np.abs(sv).mean(0), index=Xm.columns).sort_values(ascending=False),
            "interactions": paires, "efficacite_max_err": float(np.abs(sv.sum(1) + base - gbm.predict(Xm)).max()),
            "efficacite_sujet_err": float(abs(svs.sum() + base - gbm.predict(xs)[0]))}

def importance_permutation(Xh: pd.DataFrame, Xm: pd.DataFrame, y: np.ndarray, plis: list[np.ndarray], repetitions: int = 10, graine: int = GRAINE) -> pd.DataFrame:
    """Importance par permutation (Breiman 2001) mesurée HORS échantillon : pour chaque pli, la variable j est permutée dans les
    ventes testées et l'accroissement de l'erreur quadratique (log) est moyenné. Boosting et forêt."""
    rng = np.random.default_rng(graine); n = len(y); res = {}
    for nom, ctor in (("Gradient boosting", _boosting), ("Forêt aléatoire", _foret)):
        gain = np.zeros(Xm.shape[1]); base_err = 0.0
        for f in plis:
            tr = np.setdiff1d(np.arange(n), f); m = ctor(); m.fit(Xm.iloc[tr].to_numpy(), y[tr]); Xt = Xm.iloc[f].to_numpy(); e0 = np.mean((m.predict(Xt) - y[f]) ** 2); base_err += e0 * len(f) / n
            for j in range(Xm.shape[1]):
                acc = 0.0
                for _ in range(repetitions):
                    Xp = Xt.copy(); Xp[:, j] = rng.permutation(Xp[:, j]); acc += np.mean((m.predict(Xp) - y[f]) ** 2)
                gain[j] += (acc / repetitions - e0) * len(f) / n
        res[nom] = pd.Series(gain, index=Xm.columns); res[nom + " (EQM de base)"] = base_err
    return res

def dependance_partielle(modele, X: pd.DataFrame, variable: str, grille: np.ndarray) -> np.ndarray:
    """PDP (Friedman 2001) : moyenne des prédictions lorsque la variable est fixée à chaque valeur de la grille, les autres variables
    gardant leurs valeurs observées. Retourne la courbe en logarithme du prix."""
    out = np.empty(len(grille))
    for i, g in enumerate(grille):
        Xg = X.copy(); Xg[variable] = g; out[i] = float(np.mean(modele.predict(Xg)))
    return out

def pdp_vraie(v: pd.DataFrame, variable: str, grille: np.ndarray, scenario: bool = False) -> np.ndarray:
    """Dépendance partielle de la LOI DE GÉNÉRATION (connue) : utile pour juger les modèles sur données synthétiques."""
    a0 = v["age"].to_numpy(float); S0 = v["superficie_habitable_m2"].to_numpy(float); q = v["qualite_1a5"].to_numpy(float); out = np.empty(len(grille))
    for i, g in enumerate(grille):
        a = np.full_like(a0, g) if variable == "age" else a0; S = np.full_like(S0, g) if variable == "superficie" else S0
        f = VRAIS["ln_S"] * np.log(S / 140.0) + VRAIS["age"] * a + VRAIS["age2"] * a ** 2
        if scenario:
            f = f + SCENARIO["prime_qualite_age"] * ((q >= SCENARIO["qualite_min"]) & (a >= SCENARIO["age_min_prime"])) + SCENARIO["mur"] * ((a >= SCENARIO["age_mur"][0]) & (a <= SCENARIO["age_mur"][1]))
        out[i] = float(np.mean(f))
    return out

def courbes_apprentissage(Xm: pd.DataFrame, y: np.ndarray, plis: list[np.ndarray], profondeurs=(1, 2, 4)) -> dict:
    """Erreur quadratique (log) d'entraînement et de validation croisée selon le nombre d'arbres du boosting (staged_predict), pour
    plusieurs profondeurs : surapprentissage et arrêt précoce."""
    n = len(y); out = {}
    for d in profondeurs:
        tr_err = np.zeros(1200); te_err = np.zeros(1200)
        for f in plis:
            tr = np.setdiff1d(np.arange(n), f); m = _boosting(max_depth=d).fit(Xm.iloc[tr], y[tr])
            for k, (pt, pv) in enumerate(zip(m.staged_predict(Xm.iloc[tr]), m.staged_predict(Xm.iloc[f]))):
                tr_err[k] += np.mean((pt - y[tr]) ** 2) / len(plis); te_err[k] += np.mean((pv - y[f]) ** 2) * len(f) / n
        out[d] = {"entrainement": np.sqrt(tr_err), "validation": np.sqrt(te_err), "optimum": int(np.argmin(te_err)) + 1, "rmse_min": float(np.sqrt(te_err.min()))}
    return out

def sensibilite_hyperparametres(Xm: pd.DataFrame, y: np.ndarray, plis: list[np.ndarray]) -> pd.DataFrame:
    """Grille d'hyperparamètres du boosting (profondeur, taux d'apprentissage, nombre d'arbres, taille minimale des feuilles) et de la
    forêt (feuilles, variables par division) : RMSE log et COD en validation croisée aléatoire."""
    lignes = []
    grille_b = [dict(max_depth=1, n_estimators=1200), dict(max_depth=2, n_estimators=300), dict(max_depth=2, n_estimators=1200), dict(max_depth=2, n_estimators=3000),
                dict(max_depth=3, n_estimators=1200), dict(max_depth=4, n_estimators=1200), dict(max_depth=2, n_estimators=1200, learning_rate=0.1),
                dict(max_depth=2, n_estimators=1200, min_samples_leaf=1), dict(max_depth=2, n_estimators=1200, subsample=1.0)]
    for k in grille_b:
        p = valider(lambda: _boosting(**k), "ml", None, Xm, y, plis); m = metriques(p, y, facteur_duan(p, y))
        lignes.append({"modele": "Gradient boosting", "profondeur": k.get("max_depth", 2), "arbres": k.get("n_estimators", 1200), "nu": k.get("learning_rate", 0.03),
                       "feuille_min": k.get("min_samples_leaf", 5), "sous_ech": k.get("subsample", 0.8), "variables": "toutes", **m})
    for k in [dict(min_samples_leaf=1, max_features=1.0), dict(min_samples_leaf=2, max_features=0.5), dict(min_samples_leaf=5, max_features=0.33), dict(min_samples_leaf=2, max_features=0.5, n_estimators=100)]:
        p = valider(lambda: _foret(**k), "ml", None, Xm, y, plis); m = metriques(p, y, facteur_duan(p, y))
        lignes.append({"modele": "Forêt aléatoire", "profondeur": "max", "arbres": k.get("n_estimators", N_ARBRES_FORET), "nu": np.nan, "feuille_min": k["min_samples_leaf"],
                       "sous_ech": "bootstrap", "variables": f"{int(round(k['max_features'] * 100))} %", **m})
    return pd.DataFrame(lignes)

# ----------------------------------------------------------------------------- sujet
def prediction_sujet(modeles_ajustes: dict, duan: dict, Xh_s: pd.DataFrame, Xm_s: pd.DataFrame, ref: dict, s: dict) -> pd.DataFrame:
    """Prédiction du sujet par chaque modèle estimé sur les 600 ventes : logarithme, niveau avec correction de retransformation
    (facteur hors échantillon ; pour l'hédonique, traite.cas1.predire donne la valeur de référence 637 953 $)."""
    rows = []
    for nom, m in modeles_ajustes.items():
        lp = float(m.predict(Xh_s if nom.startswith("Hédonique") else Xm_s)[0]); F = duan[nom]
        rows.append({"modele": nom, "ln_pred": lp, "naif": math.exp(lp), "facteur": F, "prediction": math.exp(lp) * F})
    df = pd.DataFrame(rows).set_index("modele")
    df.loc["Hédonique (référence)", "prediction_cas1"] = float(cas1.predire(ref, cas1.sujet_df(s))[0])
    return df

# ----------------------------------------------------------------------------- exercices
def exercices(v: pd.DataFrame, comp: dict, secteurs: dict, scen: dict, sh: dict) -> dict:
    e = {}
    # 37.1 : division d'un nœud sur huit ventes (ln P, superficie) — réduction de la somme des carrés
    sub = v.sort_values("id_vente").iloc[:8]; lnP = np.log(sub["prix_vente"].to_numpy(float)); S = sub["superficie_habitable_m2"].to_numpy(float)
    ordre = np.argsort(S); S, lnP = S[ordre], lnP[ordre]; sse0 = float(np.sum((lnP - lnP.mean()) ** 2)); best = None; essais = []
    for k in range(1, 8):
        g, d = lnP[:k], lnP[k:]; sse = float(np.sum((g - g.mean()) ** 2) + np.sum((d - d.mean()) ** 2)); seuil = float((S[k - 1] + S[k]) / 2)
        essais.append({"seuil": seuil, "sse": sse, "gain": sse0 - sse})
        if best is None or sse < best["sse"]: best = {"seuil": seuil, "sse": sse, "gain": sse0 - sse, "k": k, "c_gauche": float(g.mean()), "c_droite": float(d.mean())}
    e["1"] = {"ids": list(sub.sort_values("superficie_habitable_m2")["id_vente"]), "S": S.tolist(), "lnP": lnP.tolist(), "prix": sub.sort_values("superficie_habitable_m2")["prix_vente"].tolist(),
              "sse0": sse0, "essais": essais, "best": best}
    # 37.2 : variance d'une moyenne d'arbres corrélés
    s2, rho = 0.12 ** 2, 0.30; e["2"] = {"sigma": 0.12, "rho": rho, **{f"B{B}": rho * s2 + (1 - rho) * s2 / B for B in (1, 10, 100, 500)}, "limite": rho * s2, "ecart_type_limite": math.sqrt(rho * s2),
                                        "B_pour_1pct": (1 - rho) / 0.01 / rho if rho else np.inf}
    # 37.3 : boosting à la main
    nu = 0.1; F = 12.80; yv = 13.00; traj = [F]
    for _ in range(3): F = F + nu * (yv - F); traj.append(F)
    e["3"] = {"nu": nu, "F": traj, "residus": [yv - t for t in traj], "m_1pct": math.log(0.01) / math.log(1 - nu), "F_nu1": yv}
    # 37.4 : valeurs de Shapley d'un jeu à trois joueurs (S, Q, L) ; v(S) en milliers de dollars au-dessus de la base
    vS = {frozenset(): 0.0, frozenset("S"): 60.0, frozenset("Q"): 30.0, frozenset("L"): 20.0, frozenset("SQ"): 100.0, frozenset("SL"): 85.0, frozenset("QL"): 45.0, frozenset("SQL"): 130.0}
    import itertools
    joueurs = ["S", "Q", "L"]; phi = {}
    for j in joueurs:
        tot = 0.0
        for perm in itertools.permutations(joueurs):
            avant = frozenset(perm[:perm.index(j)]); tot += vS[avant | {j}] - vS[avant]
        phi[j] = tot / 6
    e["4"] = {"v": {"+".join(sorted(k)) if k else "∅": val for k, val in vS.items()}, "phi": phi, "somme": sum(phi.values()), "total": vS[frozenset("SQL")]}
    # 37.5 : dix prédictions d'un AVM fictif
    prix = np.array([312000, 358500, 401000, 445000, 489500, 532000, 598000, 655000, 742000, 895000], float)
    pred = np.array([335000, 349000, 388000, 468000, 471000, 561000, 575000, 702000, 718000, 839000], float)
    r = h.statistiques_ratio(pred, prix); el = np.log(pred / prix)
    e["5"] = {"prix": prix.tolist(), "pred": pred.tolist(), "ratios": (pred / prix).tolist(), "rmse_log": float(np.sqrt(np.mean(el ** 2))), "mdape": float(np.median(np.abs(pred / prix - 1)) * 100),
              "mape": float(np.mean(np.abs(pred / prix - 1)) * 100), **{k: float(val) for k, val in r.items()}}
    # 37.6 : validation par secteur vs aléatoire (chiffres du script)
    e["6"] = {nom: {"aleatoire": comp["metriques"][nom]["aleatoire"]["COD"], "spatial": comp["metriques"][nom]["spatial"]["COD"], "secteur": secteurs["metriques"][nom]["COD"]}
              for nom in ("Hédonique (référence)", "Forêt aléatoire", "Gradient boosting")}
    e["6"]["biais_secteur"] = secteurs["biais_secteur"]
    # 37.7 : scénario — respécification
    e["7"] = {nom: scen["metriques"][nom]["COD"] for nom in scen["metriques"]}; e["7"]["coef"] = scen["coef_respecifie"]; e["7"]["interactions_top"] = sh["scenario_interactions"].head(3).to_dict()
    return e

# ----------------------------------------------------------------------------- orchestration
def main(verbose: bool = True, complet: bool = True) -> dict:
    d = cas1.charger(); v = d["ventes"].reset_index(drop=True); s = d["sujet"]
    ref = cas1.hedonique(v); Xh = cas1.matrice(v); Xm = matrice_ml(v); y = np.log(v["prix_vente"].astype(float).to_numpy())
    sd = cas1.sujet_df(s); Xh_s = cas1.matrice(sd)[Xh.columns]; Xm_s = matrice_ml(sd, s)
    comp = comparaison(v, Xh, Xm, y)
    det = boosting_detendance(v, Xh, Xm, y)
    sect = validation_secteurs(v, Xh, Xm, y)
    scen = scenario_nonlineaire(v, Xh, Xm, y)
    ajustes = {nom: (ctor().fit(Xh if genre == "hedo" else Xm, y)) for nom, (genre, ctor) in MODELES.items()}
    sujet = prediction_sujet(ajustes, comp["duan"], Xh_s, Xm_s, ref, s)
    sh = shap_boosting(ajustes["Gradient boosting"], Xm, Xm_s)
    sh_sc = shap_boosting(scen["gbm_scenario"], Xm, Xm_s); sh["scenario_interactions"] = sh_sc["interactions"]; sh["scenario_importance"] = sh_sc["importance"]
    pa = comp["plis"]["aleatoire"]
    perm = importance_permutation(Xh, Xm, y, pa) if complet else None
    grille_S = np.linspace(80, 260, 37); grille_A = np.linspace(1, 95, 48)
    pdp = {"superficie": {"grille": grille_S, "Gradient boosting": dependance_partielle(ajustes["Gradient boosting"], Xm, "superficie", grille_S),
                          "Forêt aléatoire": dependance_partielle(ajustes["Forêt aléatoire"], Xm, "superficie", grille_S),
                          "Hédonique (référence)": dependance_partielle(ajustes["Hédonique (référence)"], Xh, "ln_S", np.log(grille_S / 140.0)), "vraie": pdp_vraie(v, "superficie", grille_S)},
           "age": {"grille": grille_A, "Gradient boosting": dependance_partielle(ajustes["Gradient boosting"], Xm, "age", grille_A),
                   "Forêt aléatoire": dependance_partielle(ajustes["Forêt aléatoire"], Xm, "age", grille_A), "vraie": pdp_vraie(v, "age", grille_A)}}
    # hédonique : l'âge entre par deux colonnes (age, age2_100) → PDP par substitution conjointe
    def pdp_hedo_age(modele, Xh_, grille):
        out = np.empty(len(grille))
        for i, g in enumerate(grille):
            Xg = Xh_.copy(); Xg["age"] = g; Xg["age2_100"] = g ** 2 / 100.0; out[i] = float(np.mean(modele.predict(Xg)))
        return out
    pdp["age"]["Hédonique (référence)"] = pdp_hedo_age(ajustes["Hédonique (référence)"], Xh, grille_A)
    pdp["age_scenario"] = {"grille": grille_A, "Gradient boosting": dependance_partielle(scen["gbm_scenario"], Xm, "age", grille_A),
                           "Hédonique (référence)": pdp_hedo_age(scen["hedo_scenario"], Xh, grille_A), "vraie": pdp_vraie(v, "age", grille_A, scenario=True)}
    courbes = courbes_apprentissage(Xm, y, pa) if complet else None
    sens = sensibilite_hyperparametres(Xm, y, pa) if complet else None
    ex = exercices(v, comp, sect, scen, sh)
    out = {"ventes": v, "sujet": s, "modele_ref": ref, "Xh": Xh, "Xm": Xm, "y": y, "Xh_s": Xh_s, "Xm_s": Xm_s, "comparaison": comp, "detendance": det, "secteurs": sect,
           "scenario": scen, "ajustes": ajustes, "sujet_pred": sujet, "shap": sh, "permutation": perm, "pdp": pdp, "courbes": courbes, "sensibilite": sens, "exercices": ex}
    if verbose:
        pd.set_option("display.width", 220); pd.set_option("display.max_columns", 30)
        print("== Modèle hédonique de référence : R² %.3f, σ %.4f, n %d, k %d" % (ref["r2"], ref["sigma"], ref["n"], ref["k"]))
        for schema in ("aleatoire", "spatial", "temporel"):
            print(f"== Validation {schema} :")
            for nom, m in comp["metriques"].items():
                mm = m[schema]; print(f"   {nom:24s} n={mm['n']:3d} RMSElog={mm['rmse_log']:.4f} R²hors={mm['r2_hors']:.3f} biais={mm['biais_log']:+.4f} RMSE$={mm['rmse_dollars']:,.0f} MdAPE={mm['mdape']:.2f} COD={mm['COD']:.2f} PRD={mm['PRD']:.3f} PRB={mm['PRB']:.3f} méd={mm['ratio_median']:.3f} Duan={comp['duan'][nom]:.4f}")
        print("   biais mensuel 2026 :", {k: {int(kk): round(vv, 3) for kk, vv in vals.items()} for k, vals in comp["biais_mensuel_2026"].items()})
        md = det["metriques_detendance"]; m1 = det["metriques_feuille1"]
        print(f"== Boosting détendancé (test 2026) : biais={md['biais_log']:+.4f} RMSElog={md['rmse_log']:.4f} COD={md['COD']:.2f} MdAPE={md['mdape']:.2f} ; β_mois entraînement {det['beta_mois_entrainement']:.5f} ; n(mois>17)={det['n_entrainement_apres_mois_17']}")
        print(f"   Boosting feuilles ≥ 1 (test 2026) : biais={m1['biais_log']:+.4f} RMSElog={m1['rmse_log']:.4f} COD={m1['COD']:.2f}")
        print("== Un secteur exclu :")
        for nom, m in sect["metriques"].items():
            print(f"   {nom:24s} RMSElog={m['rmse_log']:.4f} MdAPE={m['mdape']:.2f} COD={m['COD']:.2f} biais/secteur=", {k: round(vv, 3) for k, vv in sect["biais_secteur"][nom].items()})
        print(f"== Scénario non linéaire (n prime={scen['n_prime']}, n mur={scen['n_mur']}) :")
        for nom, m in scen["metriques"].items():
            print(f"   {nom:24s} RMSElog={m['rmse_log']:.4f} MdAPE={m['mdape']:.2f} COD={m['COD']:.2f} PRB={m['PRB']:.3f}")
        print("   coefficients respécifiés :", {k: round(vv, 4) for k, vv in scen["coef_respecifie"].items()})
        print("== Sujet :"); print(sujet.round(4).to_string())
        print(f"== SHAP (boosting) : base {sh['base']:.4f} (exp {math.exp(sh['base']):,.0f} $) ; efficacité max err {sh['efficacite_max_err']:.2e}")
        print("   importance moyenne |φ| :", sh["importance"].round(4).to_dict())
        print("   sujet φ :", sh["sujet"].sort_values().round(4).to_dict())
        print("   interactions (base) :", sh["interactions"].head(4).round(4).to_dict()); print("   interactions (scénario) :", sh["scenario_interactions"].head(4).round(4).to_dict())
        if perm is not None:
            print("== Importance par permutation (ΔEQM log, hors échantillon) :")
            for nom in ("Gradient boosting", "Forêt aléatoire"): print(f"   {nom} (EQM base {perm[nom + ' (EQM de base)']:.5f}) :", perm[nom].sort_values(ascending=False).round(5).head(8).to_dict())
        if courbes is not None:
            print("== Courbes d'apprentissage :", {d_: (c["optimum"], round(c["rmse_min"], 4), round(float(c["entrainement"][-1]), 4)) for d_, c in courbes.items()})
        if sens is not None:
            print("== Sensibilité aux hyperparamètres :"); print(sens[["modele", "profondeur", "arbres", "nu", "feuille_min", "sous_ech", "variables", "rmse_log", "mdape", "COD", "PRB"]].round(4).to_string(index=False))
        print("== Exercices :")
        print("   37.1", {k: (round(vv, 4) if isinstance(vv, float) else vv) for k, vv in ex["1"]["best"].items()}, "SSE0", round(ex["1"]["sse0"], 4)); print("        ventes", ex["1"]["ids"], ex["1"]["S"], [round(x, 4) for x in ex["1"]["lnP"]])
        print("   37.2", {k: (round(vv, 6) if isinstance(vv, float) else vv) for k, vv in ex["2"].items()})
        print("   37.3", {k: (np.round(vv, 4).tolist() if isinstance(vv, list) else round(vv, 4)) for k, vv in ex["3"].items()})
        print("   37.4", {k: (round(vv, 3) if isinstance(vv, float) else vv) for k, vv in ex["4"].items()})
        print("   37.5", {k: (round(vv, 4) if isinstance(vv, float) else vv) for k, vv in ex["5"].items() if k not in ("prix", "pred", "ratios")}, [round(r_, 3) for r_ in ex["5"]["ratios"]])
        print("   37.6", ex["6"]); print("   37.7", ex["7"])
    return out

if __name__ == "__main__":
    main()