#!/usr/bin/env python3
"""Chapitre 14 — Nettoyage, géocodage et qualité des données (cas 1, Gatineau, DONNÉES SYNTHÉTIQUES).

1. `generer_brut()` injecte des anomalies déclarées dans le jeu propre `cas1_gatineau_ventes.csv` (graine 20261031) et écrit
   `data/synthetic/cas1_gatineau_ventes_brut.csv` (618 inscriptions), `cas1_gatineau_anomalies_injectees.csv` (table de vérité)
   et `cas1_gatineau_role_extrait.csv` (extrait de rôle fictif pour l'appariement).
2. `pipeline()` applique les règles déterministes (R1-R7), le géocodage de repli (G1), l'imputation (I1), les scores robustes (S1),
   les diagnostics de régression (levier, résidus studentisés, distance de Cook ; exclusion E1) et écrit le journal d'exclusions
   `data/synthetic/cas1_gatineau_journal_exclusions.csv`.
3. `comparer_modeles()` estime `traite.cas1.hedonique` sur le jeu brut (coercé), après chaque étape et sur le jeu propre.
4. `simulation_geocodage()`, `appariement_role()`, `comparer_imputations()`, `sensibilite_seuils()` : sections 14.3 à 14.6.

Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch14_nettoyage.py
Toutes les lois d'injection sont documentées dans data/synthetic/README-brut.md."""
from __future__ import annotations
import sys, json, pathlib, datetime
import numpy as np, pandas as pd

sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
from traite import cas1, hedonique as h

ROOT = pathlib.Path(__file__).resolve().parents[3]
D = ROOT / "data/synthetic"
GRAINE = 20261031
PI2_PAR_M2 = 10.7639
BBOX = {"lat": (45.20, 45.80), "lon": (-76.10, -75.20)}          # emprise plausible du jeu (Gatineau et périphérie)
CENTRES = {"Hull": (45.428, -75.714), "Gatineau": (45.477, -75.650), "Aylmer": (45.395, -75.843),
           "Buckingham": (45.588, -75.418), "Masson-Angers": (45.545, -75.420)}  # centres de secteur (synth.py)
SEUIL_PRIX_MIN = 5000.0        # ventes < 5 000 $ écartées (convention empruntée au Règlement sur la proportion médiane)
SEUIL_S_MAX_M2 = 450.0         # superficie habitable implausible pour une unifamiliale du jeu → test d'unité
SEUIL_MAD = 3.5                # |z robuste| > 3,5
SEUIL_T = 3.0                  # |résidu studentisé externe| > 3
C_MAD = 1.4826                 # 1/Φ^{-1}(0,75)

# ----------------------------------------------------------------------------- 1. jeu brut

def _rues(rng, n):
    noms = ["rue des Érables", "rue Laurier", "boulevard Saint-Joseph", "rue Principale", "chemin de la Montagne", "rue Notre-Dame",
            "rue Jacques-Cartier", "avenue Gatineau", "rue de Lorraine", "chemin d'Aylmer", "rue Georges", "rue Maclaren"]
    return [f"{int(rng.integers(10, 1999))}, {noms[int(rng.integers(0, len(noms)))]} (fictive)" for _ in range(n)]

def generer_brut(seed: int = GRAINE) -> dict:
    """Jeu brut = jeu propre + anomalies déclarées. Retourne brut, table de vérité, extrait de rôle (DataFrames)."""
    rng = np.random.default_rng(seed)
    v = cas1.charger()["ventes"].copy()
    n = len(v)
    # --- identifiants : matricule fictif canonique « SS-FF-NNNN-C » et adresse fictive (jamais réelles)
    code_sect = {"Hull": "05", "Gatineau": "07", "Aylmer": "03", "Buckingham": "11", "Masson-Angers": "12"}
    feuillet = rng.integers(10, 99, n); numero = rng.choice(np.arange(1000, 9999), n, replace=False); verif = rng.integers(0, 9, n)
    v["matricule"] = [f"{code_sect[s]}-{f:02d}-{m:04d}-{c}" for s, f, m, c in zip(v["secteur"], feuillet, numero, verif)]
    v["adresse_fictive"] = _rues(rng, n)
    # --- extrait de rôle fictif (une ligne par unité d'évaluation, format canonique)
    role = v[["matricule", "adresse_fictive", "secteur", "superficie_terrain_m2", "annee_construction", "valeur_role_2025"]].copy()
    role = role.rename(columns={"superficie_terrain_m2": "role_superficie_terrain_m2", "annee_construction": "role_annee_construction"})
    # --- transactions : dates promesse / acte
    delai = rng.integers(21, 91, n)
    b = pd.DataFrame({
        "id_inscription": [f"B{i+1:04d}" for i in range(n)], "id_vente_origine": v["id_vente"], "matricule": v["matricule"],
        "adresse_fictive": v["adresse_fictive"], "date_promesse": v["date_vente"], "date_acte": v["date_vente"] + pd.to_timedelta(delai, unit="D"),
        "secteur": v["secteur"], "type": v["type"], "superficie_habitable_m2": v["superficie_habitable_m2"].astype(float),
        "superficie_terrain_m2": v["superficie_terrain_m2"].astype(float), "annee_construction": v["annee_construction"].astype(float),
        "qualite_1a5": v["qualite_1a5"].astype(float), "garage_places": v["garage_places"], "sous_sol_fini": v["sous_sol_fini"], "piscine": v["piscine"],
        "chambres": v["chambres"], "salles_de_bain": v["salles_de_bain"], "lat": v["lat"].astype(float), "lon": v["lon"].astype(float),
        "valeur_role_2025": v["valeur_role_2025"], "montant_meubles": 0.0, "type_transaction": "ordinaire", "prix_vente": v["prix_vente"].astype(float)})
    # --- tirage des lignes touchées (disjointes) : 109 lignes + 18 doublons ajoutés
    plan = [("prix_symbolique", 5), ("meubles_inclus", 8), ("superficie_pi2", 12), ("dates_inversees", 10), ("coord_nulles", 9),
            ("coord_hors_zone", 6), ("vente_liee_declaree", 8), ("vente_liee_non_declaree", 6), ("terrain_manquant", 20),
            ("qualite_manquante", 15), ("annee_manquante", 10)]
    idx = rng.permutation(n); pos = 0; verite = []
    def prendre(k):
        nonlocal pos; out = idx[pos:pos + k]; pos += k; return out
    for code, k in plan:
        rows = prendre(k)
        for i in rows:
            detail = ""
            if code == "prix_symbolique":
                detail = f"prix {b.at[i,'prix_vente']:.0f} → 1"; b.at[i, "prix_vente"] = 1.0
            elif code == "meubles_inclus":
                m = float(np.round(rng.uniform(8000, 35000) / 500) * 500); b.at[i, "montant_meubles"] = m; b.at[i, "prix_vente"] += m
                detail = f"meubles {m:.0f} ajoutés au prix"
            elif code == "superficie_pi2":
                s = b.at[i, "superficie_habitable_m2"]; b.at[i, "superficie_habitable_m2"] = float(np.round(s * PI2_PAR_M2)); detail = f"{s:.0f} m² saisis en pi²"
            elif code == "dates_inversees":
                b.at[i, "date_promesse"], b.at[i, "date_acte"] = b.at[i, "date_acte"], b.at[i, "date_promesse"]; detail = "promesse ↔ acte"
            elif code == "coord_nulles":
                b.at[i, "lat"] = 0.0; b.at[i, "lon"] = 0.0; detail = "lat/lon = 0"
            elif code == "coord_hors_zone":
                if rng.random() < 0.5:
                    b.at[i, "lon"] = -b.at[i, "lon"]; detail = "signe de la longitude inversé"
                else:
                    b.at[i, "lat"], b.at[i, "lon"] = b.at[i, "lon"], b.at[i, "lat"]; detail = "lat ↔ lon"
            elif code in ("vente_liee_declaree", "vente_liee_non_declaree"):
                f = float(rng.uniform(0.40, 0.60)); p = b.at[i, "prix_vente"]; b.at[i, "prix_vente"] = float(np.round(p * f / 500) * 500)
                if code == "vente_liee_declaree":
                    b.at[i, "type_transaction"] = "entre personnes liées" if rng.random() < 0.65 else "succession"
                detail = f"prix × {f:.2f}" + ("" if code == "vente_liee_declaree" else " (non déclarée)")
            elif code == "terrain_manquant":
                b.at[i, "superficie_terrain_m2"] = np.nan; detail = "terrain NA"
            elif code == "qualite_manquante":
                b.at[i, "qualite_1a5"] = np.nan; detail = "qualité NA"
            elif code == "annee_manquante":
                b.at[i, "annee_construction"] = np.nan; detail = "année NA"
            verite.append({"id_inscription": b.at[i, "id_inscription"], "id_vente_origine": b.at[i, "id_vente_origine"], "anomalie": code, "detail": detail})
    # --- doublons d'inscription : 18 copies (12 identiques, 6 avec acte décalé de 1 à 3 jours), ajoutées en fin de fichier
    dup_src = rng.choice(np.setdiff1d(np.arange(n), idx[:pos]), 18, replace=False)
    dups = b.loc[dup_src].copy(); dups["id_inscription"] = [f"B{n + j + 1:04d}" for j in range(18)]
    for j, (k, i) in enumerate(zip(dups.index, dup_src)):
        if j >= 12:
            dups.at[k, "date_acte"] = dups.at[k, "date_acte"] + pd.Timedelta(days=int(rng.integers(1, 4)))
        verite.append({"id_inscription": dups.at[k, "id_inscription"], "id_vente_origine": b.at[i, "id_vente_origine"], "anomalie": "doublon",
                       "detail": "copie de " + b.at[i, "id_inscription"] + ("" if j < 12 else " (acte décalé)")})
    b = pd.concat([b, dups], ignore_index=True)
    # --- bruit de format sur la clé matricule des transactions (appariement, §14.4) : la vérité reste dans id_vente_origine
    u = rng.random(len(b)); mat = b["matricule"].tolist()
    for i, (x, m) in enumerate(zip(u, mat)):
        if x < 0.04: mat[i] = np.nan                                   # clé absente
        elif x < 0.09: d = list(m.replace("-", "")); p = int(rng.integers(0, len(d))); d[p] = str((int(d[p]) + 1) % 10); mat[i] = "-".join(["".join(d[0:2]), "".join(d[2:4]), "".join(d[4:8]), d[8]])  # faute de frappe
        elif x < 0.19: mat[i] = m.replace("-", "")                     # sans séparateurs
        elif x < 0.44: mat[i] = m.replace("-", " ")                    # espaces
    b["matricule"] = mat
    b = b.sort_values(["date_acte", "id_inscription"]).reset_index(drop=True)
    for c in ("date_promesse", "date_acte"):
        b[c] = pd.to_datetime(b[c]).dt.strftime("%Y-%m-%d")
    verite = pd.DataFrame(verite)
    return {"brut": b, "verite": verite, "role": role.reset_index(drop=True)}

def ecrire_brut(seed: int = GRAINE) -> dict:
    g = generer_brut(seed)
    g["brut"].to_csv(D / "cas1_gatineau_ventes_brut.csv", index=False)
    g["verite"].to_csv(D / "cas1_gatineau_anomalies_injectees.csv", index=False)
    g["role"].to_csv(D / "cas1_gatineau_role_extrait.csv", index=False)
    return g

def charger_brut() -> dict:
    b = pd.read_csv(D / "cas1_gatineau_ventes_brut.csv", parse_dates=["date_promesse", "date_acte"])
    v = pd.read_csv(D / "cas1_gatineau_anomalies_injectees.csv")
    r = pd.read_csv(D / "cas1_gatineau_role_extrait.csv")
    return {"brut": b, "verite": v, "role": r}

# ----------------------------------------------------------------------------- 2. pipeline

def _dist_centre(lat, lon):
    return np.sqrt(((lat - 45.428) * 111) ** 2 + ((lon + 75.714) * 78) ** 2)

def _journal(j, df, ids, etape, action, motif, champ="", avant=None, apres=None):
    for k, i in enumerate(ids):
        j.append({"etape": etape, "id_inscription": df.at[i, "id_inscription"], "action": action, "motif": motif, "champ": champ,
                  "valeur_avant": None if avant is None else avant[k], "valeur_apres": None if apres is None else apres[k]})

def regles(b: pd.DataFrame) -> dict:
    """R1 doublons, R2 dates inversées, R3 meubles, R4 prix symbolique, R5 unité de superficie, R6 coordonnées, R7 transaction non représentative."""
    df = b.copy(); j = []; etapes = [("0 lecture", len(df))]
    df["date_promesse"] = pd.to_datetime(df["date_promesse"]); df["date_acte"] = pd.to_datetime(df["date_acte"])
    # R2 — dates : une promesse postérieure à l'acte est impossible → permutation (correction)
    m = df["date_acte"] < df["date_promesse"]
    ids = df.index[m]
    _journal(j, df, ids, "R2", "correction", "date de promesse postérieure à l'acte : champs permutés", "date_promesse",
             df.loc[ids, "date_promesse"].dt.strftime("%Y-%m-%d").tolist(), df.loc[ids, "date_acte"].dt.strftime("%Y-%m-%d").tolist())
    df.loc[m, ["date_promesse", "date_acte"]] = df.loc[m, ["date_acte", "date_promesse"]].values
    df["delai_acte_jours"] = (df["date_acte"] - df["date_promesse"]).dt.days
    # R1 — doublons : même unité (clé normalisée ou, à défaut, attributs), même promesse, même prix → on garde l'inscription la plus récente
    df["cle_unite"] = df["matricule"].astype(str).str.replace(r"\D", "", regex=True)
    sans = df["cle_unite"].isin(["", "nan"])
    df.loc[sans, "cle_unite"] = ("ATTR" + df.loc[sans, "secteur"] + df.loc[sans, "superficie_habitable_m2"].astype(str) + df.loc[sans, "superficie_terrain_m2"].astype(str) + df.loc[sans, "annee_construction"].astype(str))
    df = df.sort_values(["date_acte", "id_inscription"])
    dup = df.duplicated(subset=["cle_unite", "date_promesse", "prix_vente"], keep="last")
    _journal(j, df, df.index[dup], "R1", "exclusion", "doublon d'inscription (même clé d'unité, même promesse, même prix) : inscription la plus récente conservée")
    df = df[~dup]
    # R1b — doublon par attributs (la clé peut être absente ou mal saisie sur l'une des deux inscriptions)
    attr = ["secteur", "type", "superficie_habitable_m2", "superficie_terrain_m2", "annee_construction", "date_promesse", "prix_vente"]
    dup2 = df.duplicated(subset=attr, keep="last")
    _journal(j, df, df.index[dup2], "R1", "exclusion", "doublon d'inscription repéré par les attributs (clé d'unité absente ou discordante) : inscription la plus récente conservée")
    df = df[~dup2]; etapes.append(("R1 doublons", len(df)))
    # R3 — meubles : le prix de l'immeuble est le prix déclaré moins la contrepartie des meubles (correction)
    m = df["montant_meubles"].fillna(0) > 0; ids = df.index[m]
    _journal(j, df, ids, "R3", "correction", "contrepartie des meubles retranchée du prix", "prix_vente", df.loc[ids, "prix_vente"].tolist(),
             (df.loc[ids, "prix_vente"] - df.loc[ids, "montant_meubles"]).tolist())
    df["prix_vente"] = df["prix_vente"] - df["montant_meubles"].fillna(0)
    # R4 — prix symbolique
    m = df["prix_vente"] < SEUIL_PRIX_MIN
    _journal(j, df, df.index[m], "R4", "exclusion", f"prix inférieur à {SEUIL_PRIX_MIN:.0f} $ (prix symbolique ou donation)", "prix_vente", df.loc[m, "prix_vente"].tolist())
    df = df[~m]; etapes.append(("R4 prix symbolique", len(df)))
    # R5 — unité de superficie : valeur implausible dont la conversion pi² → m² est plausible
    s = df["superficie_habitable_m2"]; conv = s / PI2_PAR_M2
    m = (s > SEUIL_S_MAX_M2) & conv.between(40, 400); ids = df.index[m]
    _journal(j, df, ids, "R5", "correction", "superficie implausible (> 450 m²) et plausible une fois divisée par 10,7639 : saisie en pi²", "superficie_habitable_m2",
             s[m].tolist(), conv[m].round(0).tolist())
    df.loc[m, "superficie_habitable_m2"] = conv[m].round(0)
    m2 = (df["superficie_habitable_m2"] > SEUIL_S_MAX_M2)   # implausible et non convertible → exclusion (aucun cas attendu)
    _journal(j, df, df.index[m2], "R5", "exclusion", "superficie implausible non convertible", "superficie_habitable_m2", df.loc[m2, "superficie_habitable_m2"].tolist())
    df = df[~m2]
    # R6 — coordonnées nulles ou hors emprise → invalidées (géocodage de repli à l'étape G1)
    m = (df["lat"].fillna(0) == 0) | (df["lon"].fillna(0) == 0) | ~df["lat"].between(*BBOX["lat"]) | ~df["lon"].between(*BBOX["lon"])
    ids = df.index[m]
    _journal(j, df, ids, "R6", "correction", "coordonnées nulles ou hors de l'emprise du marché : invalidées, à géocoder", "lat;lon",
             [f"{a:.5f};{o:.5f}" for a, o in zip(df.loc[ids, "lat"], df.loc[ids, "lon"])])
    df["coord_valide"] = ~m; df.loc[m, ["lat", "lon"]] = np.nan
    # R7 — conditions de vente non représentatives déclarées
    m = df["type_transaction"].fillna("ordinaire").str.lower() != "ordinaire"
    _journal(j, df, df.index[m], "R7", "exclusion", "transaction déclarée non représentative (personnes liées, succession, reprise)", "type_transaction", df.loc[m, "type_transaction"].tolist())
    df = df[~m]; etapes.append(("R7 non représentatives", len(df)))
    df["date_vente"] = df["date_promesse"]
    df["age"] = 2026 - df["annee_construction"]
    return {"df": df, "journal": j, "etapes": etapes}

def geocodage_repli(df: pd.DataFrame, j: list) -> pd.DataFrame:
    """G1 : sans adresse réelle à géocoder, les coordonnées invalides reçoivent le centre de leur secteur (précision « secteur »)."""
    df = df.copy(); m = ~df["coord_valide"]; ids = df.index[m]
    df.loc[m, "lat"] = [CENTRES[s][0] for s in df.loc[m, "secteur"]]; df.loc[m, "lon"] = [CENTRES[s][1] for s in df.loc[m, "secteur"]]
    df["precision_geocodage"] = np.where(m, "centre de secteur", "point d'origine")
    _journal(j, df, ids, "G1", "imputation", "géocodage de repli au centre du secteur (précision kilométrique) ; variables de distance à utiliser avec prudence", "lat;lon", None,
             [f"{a:.5f};{o:.5f}" for a, o in zip(df.loc[ids, "lat"], df.loc[ids, "lon"])])
    df["dist_centre_km"] = _dist_centre(df["lat"], df["lon"]).round(2)
    return df

def imputer(df: pd.DataFrame, j: list | None = None, methode: str = "mediane") -> pd.DataFrame:
    """I1 : valeurs manquantes. methode ∈ {mediane (conditionnelle secteur × type), regression, indicateur, suppression}."""
    df = df.copy()
    champs = {"superficie_terrain_m2": "terrain", "qualite_1a5": "qualite", "annee_construction": "annee"}
    if methode == "suppression":
        m = df[list(champs)].isna().any(axis=1)
        if j is not None: _journal(j, df, df.index[m], "I1", "exclusion", "valeur manquante (suppression par liste)")
        df = df[~m].copy(); df["age"] = 2026 - df["annee_construction"]; return df
    for c in champs:
        m = df[c].isna()
        if not m.any(): continue
        if methode in ("mediane", "indicateur"):
            med = df.groupby(["secteur", "type"])[c].transform("median")
            val = med[m].fillna(df[c].median())
        elif methode == "regression":
            comp = df[~m]
            def X_of(d):
                X = pd.DataFrame({"const": 1.0, "ln_S": np.log(d["superficie_habitable_m2"]), "etages": (d["type"] == "à étages").astype(float)}, index=d.index)
                for s in cas1.SECTEURS[1:]: X[f"s_{s}"] = (d["secteur"] == s).astype(float)
                if c != "annee_construction": X["annee"] = d["annee_construction"].fillna(d["annee_construction"].median())
                if c != "superficie_terrain_m2": X["ln_T"] = np.log(d["superficie_terrain_m2"].fillna(d["superficie_terrain_m2"].median()))
                return X
            y = np.log(comp[c]) if c == "superficie_terrain_m2" else comp[c]
            r = h.mco(X_of(comp).values, y.values); pred = X_of(df[m]).values @ r["beta"]
            val = pd.Series(np.exp(pred) if c == "superficie_terrain_m2" else pred, index=df.index[m])
            if c == "qualite_1a5": val = val.round().clip(1, 5)
            if c == "annee_construction": val = val.round()
        if j is not None:
            _journal(j, df, df.index[m], "I1", "imputation", f"{champs[c]} manquant : imputation ({methode})", c, None, val.round(1).tolist())
        df.loc[m, c] = val.values
        if methode == "indicateur": df[f"manq_{champs[c]}"] = m.astype(float)
    df["age"] = 2026 - df["annee_construction"]
    return df

def z_robuste(x: np.ndarray) -> np.ndarray:
    x = np.asarray(x, float); med = np.nanmedian(x); mad = np.nanmedian(np.abs(x - med))
    return (x - med) / (C_MAD * mad)

def scores_robustes(df: pd.DataFrame, j: list | None = None, seuil: float = SEUIL_MAD) -> dict:
    """S1 : scores robustes (médiane / MAD) du ln du prix au m² par secteur, du ln de la superficie habitable et du ln du terrain."""
    df = df.copy()
    df["ln_prix_m2"] = np.log(df["prix_vente"] / df["superficie_habitable_m2"])
    df["z_prix_m2"] = df.groupby("secteur")["ln_prix_m2"].transform(lambda s: z_robuste(s.values))
    df["z_superficie"] = z_robuste(np.log(df["superficie_habitable_m2"].values))
    df["z_terrain"] = df.groupby("secteur")["superficie_terrain_m2"].transform(lambda s: z_robuste(np.log(s.values)))
    df["z_classique_prix_m2"] = df.groupby("secteur")["ln_prix_m2"].transform(lambda s: (s - s.mean()) / s.std(ddof=1))
    m = (df["z_prix_m2"].abs() > seuil) | (df["z_superficie"].abs() > seuil) | (df["z_terrain"].abs() > seuil)
    ids = df.index[m]
    if j is not None:
        _journal(j, df, ids, "S1", "exclusion", f"score robuste |z| > {seuil} (prix au m² par secteur, superficie ou terrain)", "z_prix_m2;z_superficie;z_terrain", None,
                 [f"{a:.2f};{b_:.2f};{c:.2f}" for a, b_, c in zip(df.loc[ids, "z_prix_m2"], df.loc[ids, "z_superficie"], df.loc[ids, "z_terrain"])])
    return {"df": df, "exclus": df[m], "garde": df[~m]}

def diagnostics(X: np.ndarray, y: np.ndarray) -> dict:
    """Levier h_ii, résidus studentisés internes r_i et externes t_i, distance de Cook D_i, DFFITS ; vérifie Σh = k."""
    X = np.asarray(X, float); y = np.asarray(y, float); n, k = X.shape
    r = h.mco(X, y); e = r["residus"]; s2 = r["sigma2"]
    XtX_inv = np.linalg.pinv(X.T @ X)
    hii = np.einsum("ij,jk,ik->i", X, XtX_inv, X)
    r_int = e / np.sqrt(s2 * (1 - hii))
    s2_i = ((n - k) * s2 - e ** 2 / (1 - hii)) / (n - k - 1)
    t_ext = e / np.sqrt(s2_i * (1 - hii))
    cook = r_int ** 2 / k * hii / (1 - hii)
    dffits = t_ext * np.sqrt(hii / (1 - hii))
    return {"h": hii, "r_int": r_int, "t_ext": t_ext, "cook": cook, "dffits": dffits, "n": n, "k": k, "somme_h": float(hii.sum()), "mco": r}

def cook_direct(X: np.ndarray, y: np.ndarray, i: int) -> float:
    """Distance de Cook par réestimation sans l'observation i (contrôle de l'identité algébrique)."""
    X = np.asarray(X, float); y = np.asarray(y, float); n, k = X.shape
    r = h.mco(X, y); m = np.ones(n, bool); m[i] = False
    b_i = h.mco(X[m], y[m])["beta"]; d = r["beta"] - b_i
    return float(d @ (X.T @ X) @ d / (k * r["sigma2"]))

def exclusions_regression(df: pd.DataFrame, j: list | None = None, seuil_t: float = SEUIL_T, max_iter: int = 2) -> dict:
    """E1 : ajuste le modèle de référence, exclut |t_i| > seuil, réitère (au plus max_iter fois)."""
    df = df.copy(); exclus = []; iterations = []
    for it in range(1, max_iter + 1):
        X = cas1.matrice(df); y = np.log(df["prix_vente"].astype(float))
        dg = diagnostics(X.values, y.values)
        df["levier"] = dg["h"]; df["t_ext"] = dg["t_ext"]; df["cook"] = dg["cook"]
        m = np.abs(dg["t_ext"]) > seuil_t
        iterations.append({"iteration": it, "n": len(df), "exclus": int(m.sum()), "sigma": float(np.sqrt(dg["mco"]["sigma2"])), "cook_max": float(dg["cook"].max()),
                           "levier_max": float(dg["h"].max()), "n_cook_4n": int((dg["cook"] > 4 / len(df)).sum()), "n_levier_2kn": int((dg["h"] > 2 * dg["k"] / len(df)).sum())})
        if not m.any(): break
        ids = df.index[m]
        if j is not None:
            _journal(j, df, ids, "E1", "exclusion", f"résidu studentisé externe |t| > {seuil_t} (itération {it})", "t_ext;cook;levier", None,
                     [f"{a:.2f};{c:.3f};{l:.3f}" for a, c, l in zip(df.loc[ids, "t_ext"], df.loc[ids, "cook"], df.loc[ids, "levier"])])
        exclus.append(df[m]); df = df[~m]
    return {"df": df, "exclus": pd.concat(exclus) if exclus else df.iloc[0:0], "iterations": iterations, "diag": dg}

def pipeline(b: pd.DataFrame | None = None, seuil_mad: float = SEUIL_MAD, seuil_t: float = SEUIL_T, imputation: str = "mediane", ecrire: bool = False) -> dict:
    """Chaîne complète : règles → géocodage de repli → imputation → scores robustes → régression → jeu nettoyé + journal."""
    if b is None: b = charger_brut()["brut"]
    R = regles(b); j = R["journal"]; etapes = R["etapes"]
    df = geocodage_repli(R["df"], j)
    df = imputer(df, j, imputation); etapes.append(("I1 imputation", len(df)))
    S = scores_robustes(df, j, seuil_mad); df = S["garde"]; etapes.append(("S1 scores robustes", len(df)))
    E = exclusions_regression(df, j, seuil_t); df = E["df"]; etapes.append(("E1 régression", len(df)))
    jour = pd.DataFrame(j, columns=["etape", "id_inscription", "action", "motif", "champ", "valeur_avant", "valeur_apres"])
    jour.insert(0, "version_pipeline", "ch14_nettoyage.py 1.0 (2026-10-31)")
    if ecrire: jour.to_csv(D / "cas1_gatineau_journal_exclusions.csv", index=False)
    return {"propre": df, "journal": jour, "etapes": etapes, "regles": R, "scores": S, "regression": E, "apres_regles": R["df"]}

# ----------------------------------------------------------------------------- 3. effet sur le modèle

def _modele(df: pd.DataFrame) -> dict:
    d = df.copy()
    if "date_vente" not in d: d["date_vente"] = pd.to_datetime(d["date_promesse"])
    if "age" not in d or d["age"].isna().any(): d["age"] = 2026 - d["annee_construction"]
    if "dist_centre_km" not in d: d["dist_centre_km"] = _dist_centre(d["lat"].fillna(45.428), d["lon"].fillna(-75.714))
    m = cas1.hedonique(d)
    return m

def brut_coerce(b: pd.DataFrame) -> pd.DataFrame:
    """Le jeu brut tel qu'un analyste pressé le passerait au modèle : lignes complètes seulement, aucune autre correction."""
    d = b.copy(); d["date_vente"] = pd.to_datetime(d["date_promesse"]); d["age"] = 2026 - d["annee_construction"]
    d["dist_centre_km"] = _dist_centre(d["lat"], d["lon"])
    d = d.dropna(subset=["superficie_terrain_m2", "qualite_1a5", "annee_construction"])
    return d

def comparer_modeles(P: dict, b: pd.DataFrame) -> dict:
    propre_ref = cas1.charger()["ventes"]
    jeux = {"brut (lignes complètes)": brut_coerce(b), "après règles R1–R7 (+G1, I1)": imputer(geocodage_repli(P["regles"]["df"], []), None, "mediane"),
            "après scores robustes S1": P["scores"]["garde"], "après régression E1 (nettoyé)": P["propre"], "jeu propre de référence": propre_ref}
    rows = []; modeles = {}
    for nom, d in jeux.items():
        m = _modele(d); modeles[nom] = m; c = m["coef"]
        rows.append({"jeu": nom, "n": m["n"], "R2": m["r2"], "sigma": m["sigma"], "ln_S": c.loc["ln_S", "coef"], "se_ln_S": c.loc["ln_S", "se_hc1"], "ln_T": c.loc["ln_T", "coef"],
                     "age": c.loc["age", "coef"], "qualite": c.loc["qualite", "coef"], "garage": c.loc["garage", "coef"], "mois_pct": 100 * m["croissance_mensuelle"],
                     "sect_Aylmer": c.loc["sect_Aylmer", "coef"], "const": c.loc["const", "coef"]})
    tab = pd.DataFrame(rows).set_index("jeu")
    # prédiction pour le sujet
    s = cas1.charger()["sujet"]; sd = cas1.sujet_df(s)
    tab["prediction_sujet"] = [float(cas1.predire(modeles[n], sd)[0]) for n in tab.index]
    # écart nettoyé vs référence en erreurs-types
    ref = modeles["jeu propre de référence"]["coef"]; net = modeles["après régression E1 (nettoyé)"]["coef"]
    ecarts = pd.DataFrame({"ref": ref["coef"], "nettoye": net["coef"], "se_ref": ref["se_hc1"]}); ecarts["ecart_en_se"] = (ecarts["nettoye"] - ecarts["ref"]) / ecarts["se_ref"]
    return {"tableau": tab, "modeles": modeles, "ecarts": ecarts}

# ----------------------------------------------------------------------------- 4. géocodage (simulation)

TECHNIQUES = {"point d'adresse": 0.005, "centroïde de parcelle": 0.015, "interpolation le long de la rue": 0.040, "centroïde de code postal (urbain)": 0.150, "centroïde de code postal (périurbain)": 1.000}  # écart-type radial, km — HYPOTHÈSES pédagogiques

def simulation_geocodage(seed: int = GRAINE, k_voisins: int = 8) -> dict:
    """Erreur de position gaussienne isotrope d'écart-type σ (par technique) appliquée aux coordonnées propres (projetées en km) :
    erreur médiane, 95e centile, part des ventes dont l'appartenance au rayon de 3 km du sujet change, part des k voisins conservés,
    variation du coefficient de distance au centre du modèle hédonique."""
    rng = np.random.default_rng(seed)
    d = cas1.charger(); v = d["ventes"].copy(); s = d["sujet"]
    xy = h.coordonnees_km(v["lat"], v["lon"], 45.428, -75.714); xs = h.coordonnees_km([s["lat"]], [s["lon"]], 45.428, -75.714)[0]
    dist_s = np.sqrt(((xy - xs) ** 2).sum(1)); dans_rayon = dist_s <= 3.0
    W0 = h.matrice_poids_knn(xy, k_voisins, standardiser=False)
    m_ref = cas1.hedonique(v); beta_ref = m_ref["coef"].loc["dist_centre", "coef"]
    rows = []
    for nom, sig in TECHNIQUES.items():
        err = rng.normal(0, sig, size=xy.shape); xy2 = xy + err
        e = np.sqrt((err ** 2).sum(1))
        d2 = np.sqrt(((xy2 - xs) ** 2).sum(1)); change = (d2 <= 3.0) != dans_rayon
        W1 = h.matrice_poids_knn(xy2, k_voisins, standardiser=False)
        conserves = (W0 * W1).sum(1) / k_voisins
        v2 = v.copy(); v2["dist_centre_km"] = np.sqrt((xy2 ** 2).sum(1))
        m2 = cas1.hedonique(v2)
        rows.append({"technique": nom, "sigma_km": sig, "erreur_mediane_m": 1000 * np.median(e), "erreur_p95_m": 1000 * np.quantile(e, 0.95),
                     "pct_rayon_change": 100 * change.mean(), "pct_voisins_conserves": 100 * conserves.mean(), "beta_dist": m2["coef"].loc["dist_centre", "coef"],
                     "sigma_modele": m2["sigma"], "erreurs_m": 1000 * e})
    tab = pd.DataFrame(rows).set_index("technique")
    return {"tableau": tab.drop(columns="erreurs_m"), "erreurs": {r["technique"]: r["erreurs_m"] for r in rows}, "beta_dist_ref": beta_ref, "n_dans_rayon": int(dans_rayon.sum())}

# ----------------------------------------------------------------------------- 5. appariement rôle – transactions

def normaliser_matricule(x) -> str:
    if pd.isna(x): return ""
    return "".join(ch for ch in str(x) if ch.isdigit())

def appariement_role(b: pd.DataFrame, role: pd.DataFrame, verite: pd.DataFrame | None = None) -> dict:
    """Trois paliers : (1) clé exacte, (2) clé normalisée (chiffres seulement), (3) blocage par secteur + année + terrain ± 2 % avec score
    de concordance (un seul candidat exigé). Taux par palier et, si la vérité est connue (id_vente_origine ↔ matricule canonique), faux appariements."""
    t = b.copy(); r = role.copy()
    r["cle"] = r["matricule"].map(normaliser_matricule)
    t["cle"] = t["matricule"].map(normaliser_matricule)
    t["palier"] = ""; t["matricule_apparie"] = pd.Series([None] * len(t), index=t.index, dtype=object)
    # palier 1 : exact
    exact = t["matricule"].isin(r["matricule"]); t.loc[exact, "palier"] = "1 exact"; t.loc[exact, "matricule_apparie"] = t.loc[exact, "matricule"]
    # palier 2 : normalisé
    m2 = (~exact) & t["cle"].isin(r["cle"]) & (t["cle"] != "")
    cle_to_mat = dict(zip(r["cle"], r["matricule"]))
    t.loc[m2, "palier"] = "2 normalisé"; t.loc[m2, "matricule_apparie"] = t.loc[m2, "cle"].map(cle_to_mat)
    # palier 3 : blocage + concordance (secteur, année, terrain ±2 %) — un seul candidat
    reste = t.index[t["palier"] == ""]
    n3 = 0; ambigus = 0
    for i in reste:
        row = t.loc[i]
        cand = r[(r["secteur"] == row["secteur"]) & (r["role_annee_construction"] == row["annee_construction"])]
        if pd.notna(row["superficie_terrain_m2"]):
            cand = cand[(cand["role_superficie_terrain_m2"] - row["superficie_terrain_m2"]).abs() <= 0.02 * row["superficie_terrain_m2"]]
        else:
            cand = cand.iloc[0:0]
        if len(cand) == 1:
            t.at[i, "palier"] = "3 attributs"; t.at[i, "matricule_apparie"] = cand["matricule"].iloc[0]; n3 += 1
        elif len(cand) > 1:
            ambigus += 1
    t.loc[t["palier"] == "", "palier"] = "non apparié"
    taux = t["palier"].value_counts().reindex(["1 exact", "2 normalisé", "3 attributs", "non apparié"]).fillna(0).astype(int)
    out = {"table": t, "taux": taux, "pct": (100 * taux / len(t)).round(1), "ambigus": ambigus}
    if verite is not None:
        # vérité : matricule canonique de l'unité d'origine
        vrai = dict(zip(role["matricule"], role["matricule"]))
        orig = cas1.charger()["ventes"][["id_vente"]].copy()
        orig["matricule_vrai"] = role["matricule"].values     # même ordre que le jeu propre (generer_brut)
        t = t.merge(orig.rename(columns={"id_vente": "id_vente_origine"}), on="id_vente_origine", how="left")
        app = t[t["palier"] != "non apparié"]
        out["faux_appariements"] = int((app["matricule_apparie"] != app["matricule_vrai"]).sum())
        out["faux_par_palier"] = app.assign(faux=app["matricule_apparie"] != app["matricule_vrai"]).groupby("palier")["faux"].sum().to_dict()
        out["table"] = t
    return out

# ----------------------------------------------------------------------------- 5b. effet de masquage (illustration)

def illustration_masquage(b: pd.DataFrame, seuil: float = SEUIL_MAD) -> dict:
    """Sur le jeu après R1-R3 seulement (prix symboliques encore présents) : combien d'observations dépassent le seuil selon
    le score classique (moyenne, écart-type) et selon le score robuste (médiane, MAD) du ln du prix au m² par secteur ?"""
    df = b.copy(); df["date_promesse"] = pd.to_datetime(df["date_promesse"]); df["date_acte"] = pd.to_datetime(df["date_acte"])
    df["prix_net"] = df["prix_vente"] - df["montant_meubles"].fillna(0)
    df["lpm"] = np.log(df["prix_net"].clip(lower=1) / df["superficie_habitable_m2"])
    out = []
    for sct, g in df.groupby("secteur"):
        x = g["lpm"].values; zc = (x - x.mean()) / x.std(ddof=1); zr = z_robuste(x)
        out.append({"secteur": sct, "n": len(x), "moyenne": x.mean(), "ecart_type": x.std(ddof=1), "mediane": np.median(x), "mad_x1.4826": C_MAD * np.median(np.abs(x - np.median(x))),
                    "n_classique": int((np.abs(zc) > seuil).sum()), "n_robuste": int((np.abs(zr) > seuil).sum()), "n_prix_1": int((g["prix_net"] < SEUIL_PRIX_MIN).sum()),
                    "z_classique_min": zc.min(), "z_robuste_min": zr.min()})
    return {"tableau": pd.DataFrame(out).set_index("secteur"), "df": df}

# ----------------------------------------------------------------------------- 6. imputation comparée et sensibilité

def _modele_indicateurs(df: pd.DataFrame) -> dict:
    """Modèle de référence augmenté des indicatrices de valeur manquante (méthode de l'indicateur) ; prédiction du sujet avec indicatrices nulles."""
    X = cas1.matrice(df); extra = [c for c in df.columns if c.startswith("manq_")]
    for c in extra: X[c] = df[c].values
    y = np.log(df["prix_vente"].astype(float)); r = h.mco(X.values, y.values)
    coef = pd.DataFrame({"coef": r["beta"], "se_hc1": r["se_hc1"]}, index=X.columns)
    s = cas1.charger()["sujet"]; Xs = cas1.matrice(cas1.sujet_df(s)); 
    for c in extra: Xs[c] = 0.0
    pred = float(np.exp(Xs[X.columns].values @ r["beta"])[0] * np.mean(np.exp(r["residus"])))
    return {"coef": coef, "sigma": float(np.sqrt(r["sigma2"])), "n": r["n"], "r2": r["r2"], "prediction": pred}

def comparer_imputations(P: dict) -> pd.DataFrame:
    """Quatre traitements des 45 valeurs manquantes, à pipeline égal par ailleurs (S1 puis E1) : suppression par liste, médiane conditionnelle
    (secteur × type), régression auxiliaire, médiane + indicatrices de manquant dans le modèle."""
    base = geocodage_repli(P["regles"]["df"], [])
    rows = []
    for meth in ["suppression", "mediane", "regression", "indicateur"]:
        d = imputer(base, None, meth)
        S = scores_robustes(d, None); E = exclusions_regression(S["garde"], None)
        if meth == "indicateur":
            m = _modele_indicateurs(E["df"]); c = m["coef"]; pred = m["prediction"]
            ind = {k: c.loc[k, "coef"] for k in c.index if k.startswith("manq_")}
        else:
            m = _modele(E["df"]); c = m["coef"]; pred = float(cas1.predire(m, cas1.sujet_df(cas1.charger()["sujet"]))[0]); ind = {}
        rows.append({"methode": meth, "n": m["n"], "sigma": m["sigma"], "ln_T": c.loc["ln_T", "coef"], "se_ln_T": c.loc["ln_T", "se_hc1"], "qualite": c.loc["qualite", "coef"],
                     "age": c.loc["age", "coef"], "prediction_sujet": pred, **{f"coef_{k}": v for k, v in ind.items()}})
    return pd.DataFrame(rows).set_index("methode")

def variante_mnar_terrain(seed: int = GRAINE) -> pd.DataFrame:
    """Exercice 14.7 : les 20 terrains manquants ne sont plus tirés au hasard (MCAR) mais sont les 20 plus petits terrains du jeu
    parmi les lignes non altérées autrement (MNAR). Retourne le tableau de `comparer_imputations` pour cette variante (+ e.-t. des indicatrices)."""
    g = generer_brut(seed); b = g["brut"].copy(); ver = g["verite"]
    v = cas1.charger()["ventes"].set_index("id_vente")
    ids = ver.loc[ver["anomalie"] == "terrain_manquant", "id_vente_origine"]
    m = b["id_vente_origine"].isin(ids); b.loc[m, "superficie_terrain_m2"] = v.loc[b.loc[m, "id_vente_origine"], "superficie_terrain_m2"].values
    cand = b[~b["id_inscription"].isin(set(ver["id_inscription"]))].nsmallest(20, "superficie_terrain_m2")
    b.loc[cand.index, "superficie_terrain_m2"] = np.nan
    P = pipeline(b); tab = comparer_imputations(P)
    # erreur-type de l'indicatrice de terrain manquant
    base = geocodage_repli(P["regles"]["df"], []); d = imputer(base, None, "indicateur")
    S = scores_robustes(d, None); E = exclusions_regression(S["garde"], None); mi = _modele_indicateurs(E["df"])
    tab.loc["indicateur", "se_manq_terrain"] = mi["coef"].loc["manq_terrain", "se_hc1"]
    return tab

def sensibilite_seuils(b: pd.DataFrame, verite: pd.DataFrame) -> pd.DataFrame:
    """Grille seuil MAD × seuil t : exclusions, σ, β_ln S, et détection des anomalies « à exclure » (ventes liées non déclarées)."""
    a_detecter = set(verite.loc[verite["anomalie"] == "vente_liee_non_declaree", "id_inscription"])
    propres = set(b["id_inscription"]) - set(verite["id_inscription"])
    rows = []
    for zm in [2.5, 3.0, 3.5, 4.0]:
        for zt in [2.5, 3.0, 3.5]:
            P = pipeline(b, zm, zt)
            excl = set(P["journal"].loc[P["journal"]["etape"].isin(["S1", "E1"]), "id_inscription"])
            m = _modele(P["propre"])
            rows.append({"seuil_mad": zm, "seuil_t": zt, "n_final": m["n"], "exclus_S1": int((P["journal"]["etape"] == "S1").sum()), "exclus_E1": int((P["journal"]["etape"] == "E1").sum()),
                         "liees_detectees": len(excl & a_detecter), "faux_positifs": len(excl & propres), "sigma": m["sigma"], "ln_S": m["coef"].loc["ln_S", "coef"],
                         "prediction_sujet": float(cas1.predire(m, cas1.sujet_df(cas1.charger()["sujet"]))[0])})
    return pd.DataFrame(rows)

# ----------------------------------------------------------------------------- main

def main(ecrire: bool = True, verbose: bool = True) -> dict:
    pd.set_option("display.width", 220); pd.set_option("display.max_columns", 30)
    P0 = lambda *a: print(*a) if verbose else None
    g = ecrire_brut() if ecrire else generer_brut()
    b, verite, role = g["brut"], g["verite"], g["role"]
    P0("== Jeu brut :", b.shape, "| anomalies injectées :", len(verite)); P0(verite["anomalie"].value_counts().to_string())
    P = pipeline(b, ecrire=ecrire)
    P0("== Entonnoir :", P["etapes"])
    P0("== Journal par étape/action :\n", P["journal"].groupby(["etape", "action"]).size().to_string())
    P0("== Itérations E1 :", P["regression"]["iterations"])
    # détection : quelles anomalies le pipeline a-t-il traitées ?
    jv = P["journal"].merge(verite, on="id_inscription", how="left")
    P0("== Étape × anomalie :\n", pd.crosstab(jv["etape"], jv["anomalie"].fillna("(ligne propre)")).to_string())
    C = comparer_modeles(P, b)
    P0("== Modèles :\n", C["tableau"].round(4).to_string())
    P0("== Écarts nettoyé − référence (en erreurs-types de la référence) :\n", C["ecarts"].round(3).to_string())
    S = P["scores"]["df"]
    P0("== Scores robustes : max |z_prix_m2| =", round(S["z_prix_m2"].abs().max(), 2), "; exclus S1 :", len(P["scores"]["exclus"]))
    P0(P["scores"]["exclus"][["id_inscription", "secteur", "superficie_habitable_m2", "prix_vente", "z_prix_m2", "z_superficie", "z_terrain"]].round(2).to_string(index=False))
    E = P["regression"]["exclus"]
    P0("== Exclus E1 :\n", E[["id_inscription", "secteur", "superficie_habitable_m2", "prix_vente", "t_ext", "cook", "levier"]].round(3).to_string(index=False))
    M = illustration_masquage(b); P0("== Masquage (après R1-R3, prix symboliques présents) :\n", M["tableau"].round(3).to_string())
    G = simulation_geocodage(); P0("== Géocodage (simulation) :\n", G["tableau"].round(3).to_string(), "\nβ_dist référence :", round(G["beta_dist_ref"], 4), "; ventes dans le rayon de 3 km :", G["n_dans_rayon"])
    A = appariement_role(b, role, verite); P0("== Appariement :\n", A["taux"].to_string(), "\n", A["pct"].to_string(), "\nambigus :", A["ambigus"], "; faux appariements :", A.get("faux_appariements"), A.get("faux_par_palier"))
    I = comparer_imputations(P); P0("== Imputations :\n", I.round(4).to_string())
    T = sensibilite_seuils(b, verite); P0("== Sensibilité aux seuils :\n", T.round(4).to_string(index=False))
    V = variante_mnar_terrain(); P0("== Variante MNAR (ex. 14.7) :\n", V.round(4).to_string())
    return {"brut": b, "verite": verite, "role": role, "pipeline": P, "modeles": C, "geocodage": G, "appariement": A, "imputations": I, "sensibilite": T, "masquage": M, "mnar": V}

if __name__ == "__main__":
    main(ecrire="--sans-ecriture" not in sys.argv)