#!/usr/bin/env python3
"""Chapitre 20 — réconciliation et incertitude : combinaison des indications du cas 1 (Gatineau, données synthétiques).

Reprend sans les recalculer « à la main » : la sélection du chapitre 16 (traite.cas1.selection), la grille hédonique du
chapitre 17 (traite.cas1.ajustements_hedoniques), le modèle de référence et la prédiction du chapitre 18
(traite.cas1.hedonique / predire / intervalle_prediction) et le coût neuf du chapitre 21 (ch21_cout_neuf.cout_neuf_sujet).
Ajoute : bootstrap des comparables (candidats, ajustements, les deux), énumération exhaustive des sous-ensembles de cinq
comparables, indication par le coût sous hypothèses déclarées, pondération (égale, inverse de la variance, moindres carrés
généralisés avec corrélation, raisonnée), fourchette de la conclusion, et « vérité » de la loi de génération.
Imprime tous les chiffres cités dans le chapitre.
Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch20_reconciliation.py"""
from __future__ import annotations
import sys, pathlib, math, itertools, datetime
import numpy as np, pandas as pd
HERE = pathlib.Path(__file__).resolve()
sys.path.insert(0, str(HERE.parents[1]))
sys.path.insert(0, str(HERE.parent))
from traite import cas1, cout as ct, comparables as cp
import ch21_cout_neuf as c21

ROOT = HERE.parents[3]
GRAINE = 20261031
B = 2000                                   # réplications bootstrap
PLANCHER_BRUT = 0.01                       # plancher de l'inverse du brut (ch. 17, eq. 17.poids)
GROUPES = ["conditions_marche", "localisation", "superficie", "terrain", "age", "qualite", "garage", "sous_sol", "piscine", "type"]
COLS_GROUPES = {"conditions_marche": ["mois"], "localisation": ["dist_centre"] + [f"sect_{x}" for x in cas1.SECTEURS[1:]],
                "superficie": ["ln_S"], "terrain": ["ln_T"], "age": ["age", "age2_100"], "qualite": ["qualite"],
                "garage": ["garage"], "sous_sol": ["sous_sol"], "piscine": ["piscine"], "type": ["etages"]}
# Indication par le coût : reprise du chapitre 24 (ch24_depreciation.py : terrain extrait des ventes récentes de Hull, dépréciation
# extraite du marché à 28 ans + désuétude fonctionnelle). Les « hypothèses d'usage » ci-dessous ne servent que de contre-exemple.
ALPHA_TERRAIN_USAGE = 0.30  # allocation (part du terrain) — hypothèse d'usage, contre-exemple
VIE_ECONOMIQUE_USAGE = 60   # âge-vie — hypothèse d'usage, contre-exemple
AGE_EFFECTIF = 28
B_COUT = B                  # réplications du bootstrap de l'indication par le coût (mêmes tirages des 600 ventes que le modèle)
# Prédictions du sujet sous les quatre spécifications « en logarithme » du ch. 18 (tab. 18.variantes : semi-log, log-log,
# Box–Cox, indicatrices trimestrielles) : sert à l'incertitude de spécification du modèle hédonique.
PRED_SPEC_CH18 = [637953.0, 627260.0, 635171.0, 649980.0]
POIDS_RAISONNES = {"comparaison": 0.40, "hedonique": 0.45, "cout": 0.15}
ARRONDI_CONCLUSION = 5000


# ---------------------------------------------------------------- 1. indications héritées des chapitres 16 à 18
def base() -> dict:
    t = cas1.tout()
    m = t["modele"]; sel = t["selection"]; s = t["sujet"]; v = t["ventes"]
    cand = sel["candidats"].copy()                      # 16 candidats, triés par distance d'attributs croissante
    comps = sel["comparables"]
    adj16 = cas1.ajustements_hedoniques(m, cand, s)
    adj16["distance"] = cand["distance"].values
    res = pd.Series(m["residus"], index=v.index)
    adj16["residu_log"] = res.loc[cand.index].values
    adj16["retenu"] = np.arange(len(cand)) < len(comps)
    adj5 = adj16[adj16["retenu"]].copy()
    # indication par comparaison, exactement comme au chapitre 17 (bibliothèque traite.comparables)
    ca = [cp.ComparableAjuste(r["id_vente"], r["prix_vente"], {g: r[g] for g in GROUPES}) for _, r in adj5.iterrows()]
    w = cp.poids_inverse_brut(ca, PLANCHER_BRUT)
    rec = cp.reconciliation(ca, w)
    cols = list(m["X"].columns)
    Xs = cas1.matrice(cas1.sujet_df(s))[cols].values[0]
    Xc16 = cas1.matrice(cand)[cols].values
    G = np.zeros((len(cols), len(GROUPES)))
    for j, g in enumerate(GROUPES):
        for c in COLS_GROUPES[g]:
            G[cols.index(c), j] = 1.0
    return {**t, "cand": cand, "comps": comps, "adj16": adj16, "adj5": adj5, "poids_comp": np.array(w), "reconciliation": rec,
            "V_comp": rec["valeur_ponderee"], "sd_comp_pondere": rec["ecart_type_pondere"], "V_hed": t["prediction_hedonique"],
            "ip_hed": t["intervalle"], "cols": cols, "Xs": Xs, "Xc16": Xc16, "P16": cand["prix_vente"].values.astype(float), "G": G,
            "X": m["X"].values, "y": m["y"], "beta": m["beta"], "V_hc1": m["V_hc1"], "sigma": m["sigma"]}


def grille_rapide(beta: np.ndarray, Xs: np.ndarray, Xc: np.ndarray, P: np.ndarray, G: np.ndarray, plancher: float = PLANCHER_BRUT) -> dict:
    """Cascade hédonique vectorisée (identique à traite.cas1.ajustements_hedoniques) : Δ_k par groupe, bases cumulées,
    ajustements A_k = B_{k-1}(e^{Δ_k} − 1), prix ajusté P e^{ΣΔ}, brut Σ|A_k|/P, poids inverse du brut, valeur pondérée."""
    D = ((Xs[None, :] - Xc) * beta[None, :]) @ G           # n_c × n_groupes
    cum = np.cumsum(D, axis=1)
    bases = P[:, None] * np.exp(cum)
    bases_prec = np.hstack([P[:, None], bases[:, :-1]])
    A = bases_prec * (np.exp(D) - 1.0)
    P_aj = bases[:, -1]
    brut = np.abs(A).sum(axis=1) / P
    w = 1.0 / np.maximum(brut, plancher); w = w / w.sum()
    V = float(w @ P_aj)
    return {"P_aj": P_aj, "brut": brut, "poids": w, "valeur": V, "sd_pondere": float(math.sqrt(w @ (P_aj - V) ** 2))}


def valeur_ponderee(P_aj: np.ndarray, brut: np.ndarray, plancher: float = PLANCHER_BRUT) -> float:
    w = 1.0 / np.maximum(brut, plancher)
    return float((w @ P_aj) / w.sum())


# ---------------------------------------------------------------- 2. bootstrap des comparables
def bootstrap(b: dict, B: int = B, graine: int = GRAINE, k: int | None = None) -> dict:
    """Trois rééchantillonnages : (A) les 16 candidats avec remise puis sélection des 5 plus proches (modèle fixe) ;
    (B) les 600 ventes avec remise, réestimation du modèle, grille des 5 comparables fixes et prédiction du sujet ;
    (C) les deux à la fois (même réplication b pour la corrélation) ; variante paramétrique β* ~ N(β̂, V_HC1)."""
    rng = np.random.default_rng(graine)
    n_c = len(b["P16"]); k = k or len(b["comps"])
    g0 = grille_rapide(b["beta"], b["Xs"], b["Xc16"], b["P16"], b["G"])
    P_aj16, brut16 = g0["P_aj"], g0["brut"]
    # (A) candidats
    idxA = np.sort(rng.integers(0, n_c, size=(B, n_c)), axis=1)[:, :k]     # candidats triés par distance : rang = indice
    VA = np.array([valeur_ponderee(P_aj16[i], brut16[i]) for i in idxA])
    # (B) modèle (bootstrap par paires)
    X, y = b["X"], b["y"]; n = len(y)
    idxB = rng.integers(0, n, size=(B, n))
    betas = np.empty((B, X.shape[1])); duan = np.empty(B)
    for r in range(B):
        Xb, yb = X[idxB[r]], y[idxB[r]]
        beta_b, *_ = np.linalg.lstsq(Xb, yb, rcond=None)
        betas[r] = beta_b; duan[r] = float(np.mean(np.exp(yb - Xb @ beta_b)))
    pred = np.exp(betas @ b["Xs"]) * duan
    retenus = np.arange(k)
    VB = np.empty(B); VC = np.empty(B)
    for r in range(B):
        g = grille_rapide(betas[r], b["Xs"], b["Xc16"], b["P16"], b["G"])
        VB[r] = valeur_ponderee(g["P_aj"][retenus], g["brut"][retenus])
        VC[r] = valeur_ponderee(g["P_aj"][idxA[r]], g["brut"][idxA[r]])
    # variante paramétrique (coefficients tirés de la normale asymptotique, comparables fixes)
    betasP = rng.multivariate_normal(b["beta"], b["V_hc1"], size=B)
    VP = np.array([valeur_ponderee(*(lambda g: (g["P_aj"][retenus], g["brut"][retenus]))(grille_rapide(bp, b["Xs"], b["Xc16"], b["P16"], b["G"]))) for bp in betasP])
    def resume(x):
        q = np.percentile(x, [2.5, 5, 25, 50, 75, 95, 97.5])
        return {"moyenne": float(x.mean()), "ecart_type": float(x.std(ddof=1)), "q025": q[0], "q05": q[1], "q25": q[2], "q50": q[3], "q75": q[4], "q95": q[5], "q975": q[6]}
    rho = float(np.corrcoef(VC, pred)[0, 1])
    rhoB = float(np.corrcoef(VB, pred)[0, 1])
    return {"P_aj16": P_aj16, "brut16": brut16, "idxA": idxA, "idxB": idxB, "VA": VA, "VB": VB, "VC": VC, "VP": VP, "pred": pred, "betas": betas,
            "resume": {"candidats": resume(VA), "ajustements": resume(VB), "conjoint": resume(VC), "parametrique": resume(VP), "prediction": resume(pred)},
            "rho_conjoint_pred": rho, "rho_ajustements_pred": rhoB,
            "part_A_sous_hed": float((VA < b["V_hed"]).mean()), "part_C_sous_hed": float((VC < b["V_hed"]).mean()),
            "part_C_au_dessus_comp": float((VC > b["V_comp"]).mean())}


def exhaustif(b: dict, P_aj16: np.ndarray, brut16: np.ndarray, k: int = 5) -> dict:
    """Tous les sous-ensembles de k candidats parmi les 16 (C(16,5) = 4368), modèle fixe."""
    n_c = len(P_aj16)
    vals = []; n_retenus = []
    for sub in itertools.combinations(range(n_c), k):
        i = np.array(sub)
        vals.append(valeur_ponderee(P_aj16[i], brut16[i])); n_retenus.append(int((i < k).sum()))
    vals = np.array(vals); n_retenus = np.array(n_retenus)
    q = np.percentile(vals, [2.5, 5, 50, 95, 97.5])
    par_n = pd.DataFrame({"n_retenus": n_retenus, "valeur": vals}).groupby("n_retenus")["valeur"].agg(["count", "mean", "min", "max"])
    return {"valeurs": vals, "n": len(vals), "moyenne": float(vals.mean()), "ecart_type": float(vals.std(ddof=1)), "min": float(vals.min()),
            "max": float(vals.max()), "q025": q[0], "q05": q[1], "q50": q[2], "q95": q[3], "q975": q[4], "par_n_retenus": par_n,
            "rang_selection": int((vals < b["V_comp"]).sum()), "part_sous_hed": float((vals < b["V_hed"]).mean())}


def elargissement(b: dict, P_aj16: np.ndarray, brut16: np.ndarray) -> pd.DataFrame:
    """Indication par comparaison lorsque l'on retient les k premiers candidats (k = 5 … 16) : valeur pondérée (inverse du brut),
    moyenne simple, résidu moyen (exp(e) − 1) des comparables retenus. Le modèle est fixe (arbitrage biais–variance du ch. 16)."""
    res = b["adj16"]["residu_log"].values
    rows = []
    for k in range(5, len(P_aj16) + 1):
        i = np.arange(k)
        rows.append({"k": k, "valeur_ponderee": valeur_ponderee(P_aj16[i], brut16[i]), "moyenne": float(P_aj16[i].mean()),
                     "residu_moyen_pct": float(np.mean(np.exp(res[i]) - 1)), "brut_moyen": float(brut16[i].mean()),
                     "distance_max": float(b["adj16"]["distance"].values[k - 1])})
    return pd.DataFrame(rows)


# ---------------------------------------------------------------- 3. indication par le coût (reprise du chapitre 24)
def _ch24() -> dict:
    import io, contextlib
    import ch24_depreciation as c24
    with contextlib.redirect_stdout(io.StringIO()):
        return c24.main()


def cout(b: dict) -> dict:
    """Indication par le coût du chapitre 24 : V = V_L + C_N − D_extraite(28 ans) − D_fonctionnelle. Ajoute les diagnostics
    d'extraction aux deux autres indications, la lecture « hypothèses d'usage » (allocation 30 %, âge-vie 60 ans) comme
    contre-exemple, et l'incertitude de spécification (variantes d'extraction du ch. 24)."""
    import ch24_depreciation as c24
    r24 = _ch24()
    CN = r24["cn"]["total"]; vs = r24["valeur"]; V_L = r24["vl_retenue"]["retenue"]; fonc = r24["salle_eau"]["desuetude"]
    D_ext = vs["d_extraite"]; D_tot = vs["B_extraction"]["depreciation_totale"]
    # variantes d'extraction du ch. 24 (fenêtre des ventes récentes, élasticité, qualité) → valeur sous chaque variante
    sx = r24["sens_extraction"].copy()
    sx["valeur"] = [ct.valeur_par_cout(float(vl), CN, float(t) * CN, fonc)["valeur"] for vl, t in zip(sx["vl_hull"], sx["taux_28"])]
    # contre-exemple : hypothèses d'usage (allocation sur l'indication par comparaison, âge-vie)
    V_L_usage = ALPHA_TERRAIN_USAGE * b["V_comp"]; D_usage = ct.depreciation_age_vie(CN, AGE_EFFECTIF, VIE_ECONOMIQUE_USAGE)
    usage = ct.valeur_par_cout(V_L_usage, CN, D_usage)
    return {"r24": r24, "CN": CN, "V_L": V_L, "V_L_mediane": r24["vl_retenue"]["mediane"], "n_terrain": r24["vl_retenue"]["n"],
            "part_terrain": r24["vl_retenue"]["part"], "D_ext": D_ext, "taux_ext": vs["taux_extrait_28"], "D_fonc": fonc, "D_tot": D_tot,
            "taux_tot": vs["B_extraction"]["taux_depreciation"], "contribution_batiment": vs["B_extraction"]["contribution_batiment"],
            "valeur": vs["valeur_retenue"], "arrondie": vs["valeur_arrondie"], "ventilee": vs["A_ventilee"]["valeur"],
            "taux_ventilee": vs["A_ventilee"]["taux_depreciation"], "tornade": r24["tornade"], "sens_extraction": sx,
            "sd_spec": float(sx["valeur"].std(ddof=1)), "cout_neuf_plus_terrain": CN + V_L,
            "extraction_hed": ct.extraction_depreciation_marche(b["V_hed"], V_L, CN, AGE_EFFECTIF),
            "extraction_comp": ct.extraction_depreciation_marche(b["V_comp"], V_L, CN, AGE_EFFECTIF),
            "usage": {"alpha": ALPHA_TERRAIN_USAGE, "vie": VIE_ECONOMIQUE_USAGE, "V_L": V_L_usage, "D": D_usage, "valeur": usage["valeur"],
                      "extraction_hed": ct.extraction_depreciation_marche(b["V_hed"], V_L_usage, CN, AGE_EFFECTIF)},
            "fonc_fn": fonc}


def bootstrap_cout(b: dict, bs: dict, co: dict, B: int = B_COUT) -> dict:
    """Bootstrap de l'indication par le coût du ch. 24 : mêmes rééchantillons des 600 ventes que le modèle (idxB) ; pour chaque
    réplication, prix ramenés à la date par la tendance du modèle rééstimé, terrain extrait des ventes récentes de Hull,
    dépréciation extraite à 28 ans (ajustement quadratique), désuétude fonctionnelle fixe. Les ventes d'un secteur sans vente
    récente dans la réplication (terrain non extractible) sont exclues de l'ajustement du profil."""
    import ch24_depreciation as c24
    v = b["ventes"]; CN = co["CN"]; fonc = co["D_fonc"]; i_mois = b["cols"].index("mois")
    VK = np.empty(B); VL = np.empty(B); T28 = np.empty(B)
    for r in range(B):
        vb = v.iloc[bs["idxB"][r]].reset_index(drop=True)
        mb = {"croissance_mensuelle": float(np.exp(bs["betas"][r][i_mois]) - 1.0)}
        w, vl = c24.extraction(vb, mb, CN)
        w = w.dropna(subset=["taux"])          # secteurs sans vente récente dans la réplication : pas de terrain, exclus du profil
        q = c24.ajuster_profil(w)["quadratique"]; t28 = float(c24.taux_quadratique(q, AGE_EFFECTIF))
        VL[r] = float(vl.loc["Hull", "vl_mediane"]); T28[r] = t28
        VK[r] = ct.valeur_par_cout(VL[r], CN, t28 * CN, fonc)["valeur"]
    qv = np.percentile(VK, [2.5, 5, 50, 95, 97.5])
    return {"VK": VK, "VL": VL, "T28": T28,
            "resume": {"moyenne": float(VK.mean()), "ecart_type": float(VK.std(ddof=1)), "q025": qv[0], "q05": qv[1], "q50": qv[2], "q95": qv[3], "q975": qv[4]},
            "sd_VL": float(VL.std(ddof=1)), "sd_T28": float(T28.std(ddof=1)),
            "rho_cout_pred": float(np.corrcoef(VK, bs["pred"][:B])[0, 1]), "rho_cout_comp": float(np.corrcoef(VK, bs["VC"][:B])[0, 1])}


# ---------------------------------------------------------------- 4. combinaison
def poids_optimaux_deux(s1: float, s2: float, rho: float) -> float:
    """Poids optimal de l'indication 1 parmi deux indications corrélées (Bates–Granger) : w1 = (σ2² − ρσ1σ2)/(σ1² + σ2² − 2ρσ1σ2)."""
    return (s2 ** 2 - rho * s1 * s2) / (s1 ** 2 + s2 ** 2 - 2 * rho * s1 * s2)


def poids_gls(Sigma: np.ndarray) -> np.ndarray:
    """w = Σ⁻¹1 / 1'Σ⁻¹1 (variance minimale sous Σw = 1, poids non contraints)."""
    inv = np.linalg.inv(Sigma); un = np.ones(Sigma.shape[0])
    w = inv @ un
    return w / w.sum()


def combinaison(b: dict, bs: dict, co: dict, bk: dict, graine: int = GRAINE) -> dict:
    V = np.array([b["V_comp"], b["V_hed"], co["valeur"]])
    sd_spec = float(np.std(PRED_SPEC_CH18, ddof=1))
    sd_pred_boot = bs["resume"]["prediction"]["ecart_type"]
    s = np.array([bs["resume"]["conjoint"]["ecart_type"], math.sqrt(sd_pred_boot ** 2 + sd_spec ** 2), math.sqrt(bk["resume"]["ecart_type"] ** 2 + co["sd_spec"] ** 2)])
    rho = bs["rho_conjoint_pred"]; rho_ck = bk["rho_cout_comp"]; rho_hk = bk["rho_cout_pred"]
    R = np.array([[1.0, rho, rho_ck], [rho, 1.0, rho_hk], [rho_ck, rho_hk, 1.0]])
    Sigma = R * np.outer(s, s)
    schemas = {
        "comparaison seule": np.array([1.0, 0.0, 0.0]),
        "hédonique seule": np.array([0.0, 1.0, 0.0]),
        "égale (un tiers chacune)": np.ones(3) / 3,
        "inverse de la variance (indépendance)": (1 / s ** 2) / (1 / s ** 2).sum(),
        "variance minimale avec corrélation (MCG)": poids_gls(Sigma),
        "raisonnée (retenue)": np.array([POIDS_RAISONNES["comparaison"], POIDS_RAISONNES["hedonique"], POIDS_RAISONNES["cout"]]),
    }
    tab = pd.DataFrame([{"schema": k, "w_comp": w[0], "w_hed": w[1], "w_cout": w[2], "valeur": float(w @ V), "ecart_type": float(math.sqrt(w @ Sigma @ w))}
                        for k, w in schemas.items()])
    # distribution bootstrap de la conclusion retenue (réplication conjointe : comparaison et modèle issus du même tirage ; coût tiré de la grille)
    n_b = len(bk["VK"]); cout_b = bk["VK"]
    w = schemas["raisonnée (retenue)"]
    Vb = w[0] * bs["VC"][:n_b] + w[1] * bs["pred"][:n_b] + w[2] * cout_b
    q = np.percentile(Vb, [2.5, 5, 50, 95, 97.5])
    V_raisonnee = float(w @ V)
    arrondi = ARRONDI_CONCLUSION
    conclusion = round(V_raisonnee / arrondi) * arrondi
    fourchette90 = (math.floor(q[1] / arrondi) * arrondi, math.ceil(q[3] / arrondi) * arrondi)
    # fourchette de valeur conclue : encadrée par les deux indications crédibles (comparaison, hédonique), arrondies vers l'extérieur
    fourchette_conclue = (math.floor(min(V[0], V[1]) / arrondi) * arrondi, math.ceil(max(V[0], V[1]) / arrondi) * arrondi)
    # courbes de sensibilité aux poids : w_hed de 0 à 1 − w_cout, pour trois poids du coût
    wh = np.linspace(0, 1, 101)
    courbes = {wc: np.array([ (1 - h - wc) * V[0] + h * V[1] + wc * V[2] if h <= 1 - wc + 1e-12 else np.nan for h in wh]) for wc in (0.0, 0.15, 0.30)}
    # moyenne simple sous les hypothèses d'usage du coût (contre-exemple)
    V_usage = np.array([b["V_comp"], b["V_hed"], co["usage"]["valeur"]])
    # variance de la combinaison de deux indications (comparaison, hédonique) selon ρ
    var_deux = {r: np.array([((1 - h) ** 2) * s[0] ** 2 + h ** 2 * s[1] ** 2 + 2 * h * (1 - h) * r * s[0] * s[1] for h in wh]) for r in (0.0, 0.3, 0.6, 0.9)}
    w_opt_deux = {r: poids_optimaux_deux(s[0], s[1], r) for r in (0.0, rho, 0.3, 0.6, 0.9)}
    return {"V": V, "s": s, "sd_spec_hed": sd_spec, "sd_pred_boot": sd_pred_boot, "rho": rho, "R": R, "Sigma": Sigma, "schemas": schemas, "tableau": tab,
            "moyenne_simple_usage": float(V_usage.mean()), "inverse_variance_usage": None,
            "Vb": Vb, "cout_b": cout_b, "q": {"q025": q[0], "q05": q[1], "q50": q[2], "q95": q[3], "q975": q[4]}, "V_raisonnee": V_raisonnee,
            "conclusion": conclusion, "fourchette90": fourchette90, "fourchette_conclue": fourchette_conclue,
            "largeur_relative_conclue": (fourchette_conclue[1] - fourchette_conclue[0]) / conclusion, "wh": wh, "courbes": courbes, "var_deux": var_deux, "w_opt_deux": w_opt_deux,
            "ecart_conclusion_comp": conclusion / b["V_comp"] - 1, "ecart_conclusion_hed": conclusion / b["V_hed"] - 1}


# ---------------------------------------------------------------- 5. ce que révèle la loi de génération
def verite(b: dict) -> dict:
    """Espérance conditionnelle « vraie » du prix du sujet selon la loi déclarée dans data/synthetic/README.md ;
    la constante (normalisation de la médiane) est estimée sur les 600 ventes (erreur-type σ/√n ≈ 0,35 %)."""
    from traite.synth import _mois_depuis
    v = b["ventes"]; s = b["sujet"]
    med = 475000.0; med_sect = {"Gatineau": 464900.0, "Aylmer": 530000.0, "Buckingham": 395000.0, "Masson-Angers": 395000.0}
    prime = {"Hull": 0.0, **{k: math.log(m / med) for k, m in med_sect.items()}}; prime["Masson-Angers"] -= 0.02
    g = 1.0906 ** 0.5 - 1
    ref = _mois_depuis([datetime.date(2024, 7, 1)])[0]
    def f(df, dates):
        mois = _mois_depuis([d.date() for d in pd.to_datetime(dates)])
        return (0.55 * np.log(df["superficie_habitable_m2"].values / 140) + 0.10 * np.log(df["superficie_terrain_m2"].values / 560)
                - 0.0055 * df["age"].values + 0.00003 * df["age"].values ** 2 + 0.06 * (df["qualite_1a5"].values - 3)
                + 0.035 * df["garage_places"].values + 0.03 * df["sous_sol_fini"].values + 0.02 * df["piscine"].values
                + np.array([prime[x] for x in df["secteur"].values]) - 0.006 * np.clip(df["dist_centre_km"].values - 3, 0, None)
                + 0.02 * (df["type"].values == "à étages") + math.log1p(g) * (mois - ref) / 12)
    fv = f(v, v["date_vente"]); e = np.log(v["prix_vente"].values) - fv
    c_hat = float(e.mean()); se_c = float(e.std(ddof=1) / math.sqrt(len(e)))
    sd = cas1.sujet_df(s); fs = float(f(sd, sd["date_vente"])[0])
    V_true = math.exp(fs + c_hat) * math.exp(0.085 ** 2 / 2)
    return {"c_hat": c_hat, "se_c": se_c, "sigma_e": float(e.std(ddof=1)), "V_true": V_true,
            "V_true_bas": V_true * math.exp(-1.96 * se_c), "V_true_haut": V_true * math.exp(1.96 * se_c),
            "ecart_comp": b["V_comp"] / V_true - 1, "ecart_hed": b["V_hed"] / V_true - 1}


# ---------------------------------------------------------------- 6. exercices (chiffres des solutions)
def exercices(b: dict, bs: dict, co: dict, cb: dict) -> dict:
    out = {}
    # ex. 1 : deux indications indépendantes
    V1, s1, V2, s2 = 582000.0, 30000.0, 636000.0, 12000.0
    w1 = (1 / s1 ** 2) / (1 / s1 ** 2 + 1 / s2 ** 2)
    out["ex1"] = {"w1": w1, "w2": 1 - w1, "V": w1 * V1 + (1 - w1) * V2, "sd": math.sqrt(1 / (1 / s1 ** 2 + 1 / s2 ** 2)),
                  "V_egal": 0.5 * (V1 + V2), "sd_egal": math.sqrt(0.25 * s1 ** 2 + 0.25 * s2 ** 2)}
    # ex. 2 : intervalle de confiance de la valeur (E[P|x]) contre intervalle de prédiction du prix (ch. 18)
    lnp, var_est, sig2 = 13.3627, 0.00010, 0.0830 ** 2
    duan = 1.00335
    out["ex2"] = {"ic_bas": math.exp(lnp - 1.96 * math.sqrt(var_est)) * duan, "ic_haut": math.exp(lnp + 1.96 * math.sqrt(var_est)) * duan,
                  "ip_bas": math.exp(lnp - 1.96 * math.sqrt(var_est + sig2)), "ip_haut": math.exp(lnp + 1.96 * math.sqrt(var_est + sig2)),
                  "demi_largeur_ic_pct": 1.96 * math.sqrt(var_est), "demi_largeur_ip_pct": 1.96 * math.sqrt(var_est + sig2),
                  "rapport_largeurs": math.sqrt(var_est + sig2) / math.sqrt(var_est)}
    # ex. 3 : poids optimaux avec corrélation
    out["ex3"] = {r: {"w1": poids_optimaux_deux(s1, s2, r),
                      "sd": math.sqrt(max(0.0, (lambda w: w ** 2 * s1 ** 2 + (1 - w) ** 2 * s2 ** 2 + 2 * w * (1 - w) * r * s1 * s2)(poids_optimaux_deux(s1, s2, r))))}
                  for r in (0.0, 0.3, 0.5, 0.6)}
    # ex. 4 : bootstrap à la main (percentiles de 20 réplications données)
    rep = np.array([561, 568, 571, 574, 576, 577, 579, 580, 582, 583, 585, 586, 588, 590, 592, 595, 598, 603, 611, 624]) * 1000.0
    out["ex4"] = {"moyenne": float(rep.mean()), "ecart_type": float(rep.std(ddof=1)), "q05": float(np.percentile(rep, 5)), "q95": float(np.percentile(rep, 95)),
                  "q10": float(np.percentile(rep, 10)), "q90": float(np.percentile(rep, 90))}
    # ex. 5 : coût sous hypothèses d'usage, puis extraction avec le terrain du ch. 24
    CN = co["CN"]
    out["ex5"] = {"a": ct.valeur_par_cout(0.35 * b["V_comp"], CN, ct.depreciation_age_vie(CN, 22, 60))["valeur"],
                  "b": ct.valeur_par_cout(0.35 * b["V_comp"], CN, ct.depreciation_age_vie(CN, 28, 75))["valeur"],
                  "extraction_hed": ct.extraction_depreciation_marche(b["V_hed"], co["V_L"], CN, 28),
                  "extraction_comp": ct.extraction_depreciation_marche(b["V_comp"], co["V_L"], CN, 28),
                  "V_L_pour_hed": b["V_hed"] - (CN - co["D_tot"])}
    # ex. 6 : autres pondérations et comparaison élargie (k = 16)
    V16 = float(r_el["valeur_ponderee"].iloc[-1]) if (r_el := _CACHE.get("el")) is not None else float(elargissement(b, bs["P_aj16"], bs["brut16"])["valeur_ponderee"].iloc[-1])
    w2 = np.array([0.30, 0.55, 0.15]); w1 = np.array([POIDS_RAISONNES["comparaison"], POIDS_RAISONNES["hedonique"], POIDS_RAISONNES["cout"]])
    V = np.array([b["V_comp"], b["V_hed"], co["valeur"]]); V16v = np.array([V16, b["V_hed"], co["valeur"]])
    a = ARRONDI_CONCLUSION
    bs8 = bootstrap(b, B=B, k=8)
    out["ex6"] = {"sd_conjoint_k8": bs8["resume"]["conjoint"]["ecart_type"], "V_k8": float(valeur_ponderee(bs["P_aj16"][:8], bs["brut16"][:8])),
                  "V_30_55_15": float(w2 @ V), "arrondi_30_55_15": round(float(w2 @ V) / a) * a, "V16": V16,
                  "V_k16_raisonnee": float(w1 @ V16v), "arrondi_k16": round(float(w1 @ V16v) / a) * a,
                  "V_55_30_15": float(np.array([0.55, 0.30, 0.15]) @ V), "arrondi_55_30_15": round(float(np.array([0.55, 0.30, 0.15]) @ V) / a) * a}
    # ex. 7 : incidence d'une hypothèse extraordinaire (eq. 10.incidence / 10.esperance)
    out["ex7"] = {"V_H": cb["conclusion"], "V_nonH": cb["conclusion"] - 14000.0, "p": 0.7}
    out["ex7"]["delta"] = out["ex7"]["V_H"] - out["ex7"]["V_nonH"]
    out["ex7"]["esperance"] = out["ex7"]["V_nonH"] + out["ex7"]["p"] * out["ex7"]["delta"]
    return out


# ---------------------------------------------------------------- tout
_CACHE: dict = {}

def tout(verbose: bool = False) -> dict:
    if "r" in _CACHE:
        return _CACHE["r"]
    b = base()
    bs = bootstrap(b)
    ex = exhaustif(b, bs["P_aj16"], bs["brut16"])
    el = elargissement(b, bs["P_aj16"], bs["brut16"]); _CACHE["el"] = el
    co = cout(b)
    bk = bootstrap_cout(b, bs, co)
    cb = combinaison(b, bs, co, bk)
    vt = verite(b)
    exo = exercices(b, bs, co, cb)
    r = {"base": b, "bootstrap": bs, "exhaustif": ex, "elargissement": el, "cout": co, "bootstrap_cout": bk, "combinaison": cb, "verite": vt, "exercices": exo}
    _CACHE["r"] = r
    if verbose:
        imprimer(r)
    return r


def imprimer(r: dict) -> None:
    b, bs, ex, co, cb, vt = r["base"], r["bootstrap"], r["exhaustif"], r["cout"], r["combinaison"], r["verite"]
    f = lambda x: f"{x:,.0f}".replace(",", " ")
    print("=== Indications héritées ===")
    print(f"Comparaison (ch. 17) : pondérée {f(b['V_comp'])} $, écart-type pondéré {f(b['sd_comp_pondere'])} $, poids {np.round(b['poids_comp'], 3)}")
    ip = b["ip_hed"]
    print(f"Hédonique (ch. 18) : {f(b['V_hed'])} $ (Duan) ; IP 95 % [{f(ip['bas'])} ; {f(ip['haut'])}]")
    print("\n=== 16 candidats : prix ajusté, brut, résidu (%) ===")
    a = b["adj16"]
    print(a[["id_vente", "distance", "prix_vente", "prix_ajuste", "ajust_brut_pct", "residu_log", "retenu"]].assign(residu_pct=lambda d: np.exp(d["residu_log"]) - 1).round(4).to_string(index=False))
    print(f"résidu moyen (exp−1) des 5 retenus : {np.mean(np.exp(a.loc[a['retenu'], 'residu_log']) - 1):.4f} ; des 11 autres : {np.mean(np.exp(a.loc[~a['retenu'], 'residu_log']) - 1):.4f} ; des 16 : {np.mean(np.exp(a['residu_log']) - 1):.4f}")
    print("\n=== Bootstrap (B = %d) ===" % B)
    for k, v in bs["resume"].items():
        print(f"{k:14s} moy {f(v['moyenne'])} é.-t. {f(v['ecart_type'])} q2.5 {f(v['q025'])} q5 {f(v['q05'])} q50 {f(v['q50'])} q95 {f(v['q95'])} q97.5 {f(v['q975'])}")
    print(f"corrélation (conjoint, prédiction) = {bs['rho_conjoint_pred']:.3f} ; (ajustements, prédiction) = {bs['rho_ajustements_pred']:.3f}")
    print(f"part des rééchantillons de candidats sous la prédiction hédonique : {bs['part_A_sous_hed']:.3f} ; conjoint : {bs['part_C_sous_hed']:.3f}")
    print("\n=== Énumération exhaustive C(16,5) ===")
    print(f"n = {ex['n']} ; moyenne {f(ex['moyenne'])} ; é.-t. {f(ex['ecart_type'])} ; min {f(ex['min'])} ; max {f(ex['max'])} ; q5 {f(ex['q05'])} ; q50 {f(ex['q50'])} ; q95 {f(ex['q95'])}")
    print(f"rang de la sélection du ch. 16 : {ex['rang_selection']} sous-ensembles sur {ex['n']} donnent moins ; part sous la prédiction hédonique {ex['part_sous_hed']:.3f}")
    print(ex["par_n_retenus"].round(0).to_string())
    print("\n=== Élargissement k = 5 … 16 (modèle fixe) ===")
    print(r["elargissement"].round(4).to_string(index=False))
    bk = r["bootstrap_cout"]
    print("\n=== Coût (chapitre 24) ===")
    print(f"C_N {f(co['CN'])} ; V_L {f(co['V_L'])} (médiane {f(co['V_L_mediane'])} de {co['n_terrain']} ventes récentes de Hull, part {co['part_terrain']:.3f}) ; D extraite {f(co['D_ext'])} ({co['taux_ext']:.4f}) ; D fonctionnelle {f(co['D_fonc'])} ; D totale {f(co['D_tot'])} ({co['taux_tot']:.4f}) ; contribution bâtiment {f(co['contribution_batiment'])} ; valeur {f(co['valeur'])} → {f(co['arrondie'])} ; lecture ventilée {f(co['ventilee'])} ({co['taux_ventilee']:.4f})")
    print(f"C_N + V_L = {f(co['cout_neuf_plus_terrain'])} ; extraction à la prédiction hédonique : D = {f(co['extraction_hed']['depreciation'])} (taux {co['extraction_hed']['taux']:.4f}) ; à la comparaison : D = {f(co['extraction_comp']['depreciation'])} (taux {co['extraction_comp']['taux']:.4f})")
    u = co["usage"]
    print(f"hypothèses d'usage (allocation {u['alpha']:.2f} × comparaison = {f(u['V_L'])}, âge-vie {u['vie']} ans : D = {f(u['D'])}) : valeur {f(u['valeur'])} ; extraction au modèle avec ce terrain : D = {f(u['extraction_hed']['depreciation'])} ({u['extraction_hed']['taux']:.3f})")
    print("variantes d'extraction (ch. 24) :"); print(co["sens_extraction"][["variante", "vl_hull", "taux_28", "valeur"]].round(4).to_string(index=False)); print(f"écart-type de spécification {f(co['sd_spec'])}")
    print("tornade ch. 24 :"); [print(f"  {t['param']:50s} {f(t['bas'])} – {f(t['haut'])}") for t in co["tornade"]]
    rk = bk["resume"]
    print(f"bootstrap du coût (B = {len(bk['VK'])}) : moy {f(rk['moyenne'])} é.-t. {f(rk['ecart_type'])} q5 {f(rk['q05'])} q50 {f(rk['q50'])} q95 {f(rk['q95'])} ; é.-t. V_L {f(bk['sd_VL'])} ; é.-t. taux 28 ans {bk['sd_T28']:.4f} ; ρ(coût, prédiction) {bk['rho_cout_pred']:.3f} ; ρ(coût, comparaison) {bk['rho_cout_comp']:.3f}")
    print("\n=== Combinaison ===")
    print(f"indications : {[f(x) for x in cb['V']]} ; écarts-types : {[f(x) for x in cb['s']]} (hédonique : boot {f(cb['sd_pred_boot'])} + spécification {f(cb['sd_spec_hed'])} ; coût : boot {f(bk['resume']['ecart_type'])} + spécification {f(co['sd_spec'])})")
    print("matrice de corrélation :"); print(np.round(cb["R"], 3))
    print(f"moyenne simple avec le coût sous hypothèses d'usage : {f(cb['moyenne_simple_usage'])}")
    print(cb["tableau"].round(3).to_string(index=False))
    print(f"poids optimal de la comparaison à deux indications selon ρ : { {round(k,3): round(v,3) for k, v in cb['w_opt_deux'].items()} }")
    q = cb["q"]
    print(f"conclusion raisonnée {f(cb['V_raisonnee'])} → arrondie {f(cb['conclusion'])} ; bootstrap de la combinaison : q2.5 {f(q['q025'])} q5 {f(q['q05'])} q50 {f(q['q50'])} q95 {f(q['q95'])} q97.5 {f(q['q975'])} ; fourchette 90 % arrondie {f(cb['fourchette90'][0])}–{f(cb['fourchette90'][1])} ; fourchette conclue {f(cb['fourchette_conclue'][0])}–{f(cb['fourchette_conclue'][1])} ({cb['largeur_relative_conclue']:.3f})")
    print(f"écart conclusion / comparaison {cb['ecart_conclusion_comp']:.4f} ; / hédonique {cb['ecart_conclusion_hed']:.4f}")
    print("\n=== Loi de génération ===")
    print(f"c_hat {vt['c_hat']:.4f} (e.-t. {vt['se_c']:.4f}, σ_e {vt['sigma_e']:.4f}) ; V_vraie {f(vt['V_true'])} [{f(vt['V_true_bas'])} ; {f(vt['V_true_haut'])}] ; écart comparaison {vt['ecart_comp']:.4f} ; hédonique {vt['ecart_hed']:.4f} ; conclusion {cb['conclusion']/vt['V_true']-1:.4f}")
    print("\n=== Exercices ===")
    for k, v in r["exercices"].items():
        print(k, {kk: (round(vv, 4) if isinstance(vv, float) else vv) for kk, vv in v.items()} if isinstance(v, dict) else v)


if __name__ == "__main__":
    tout(verbose=True)