# -*- coding: utf-8 -*-
"""05 — Analyse : métriques par configuration (validation 2024 → choix ; test 2025-26 → verdict),
surface Erreur(k, d, t), effets des facteurs (sélection, ajustement, pondération, réconciliation),
compétition de modèles M1-M13, hypothèses H1-H5, strates (région, densité, prix), intervalles,
règles conditionnelles. Sorties : tables/*.csv, results.json, fig/*.pdf"""
from scipy.optimize import nnls
import matplotlib; matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib import rcParams
from common import *

# ---------------------------------------------------------------- style (palette de référence, marques fines)
PAL = ["#2a78d6", "#eb6834", "#1baf7a", "#eda100", "#e87ba4", "#008300", "#6250d6", "#e34948"]
SEQ = ["#cde2fb", "#9ec5f4", "#6da7ec", "#3987e5", "#256abf", "#184f95", "#0d366b"]
rcParams.update({"font.family": "Helvetica Neue", "font.size": 9, "axes.edgecolor": "#c3c2b7", "axes.linewidth": 0.6,
                 "axes.grid": True, "grid.color": "#e8e7e3", "grid.linewidth": 0.5, "axes.spines.top": False,
                 "axes.spines.right": False, "legend.frameon": False, "figure.dpi": 150, "savefig.bbox": "tight",
                 "axes.titleweight": "bold", "axes.titlesize": 10, "axes.labelcolor": "#52514e", "xtick.color": "#52514e",
                 "ytick.color": "#52514e", "figure.facecolor": "#fcfcfb", "axes.facecolor": "#fcfcfb"})
from matplotlib.colors import LinearSegmentedColormap
CMAP = LinearSegmentedColormap.from_list("seqblue", SEQ)
TAB = OUT / "tables"; TAB.mkdir(exist_ok=True); FIG = OUT / "fig"
def savefig(fig, name): fig.savefig(FIG / f"{name}.pdf"); fig.savefig(FIG / f"{name}.png"); plt.close(fig)

# ---------------------------------------------------------------- chargement
T = load_parquet("targets.parquet"); S = load_parquet("sales.parquet").set_index("idx")
E1 = pd.read_csv(OUT / "configs_E1.csv"); E2 = pd.read_csv(OUT / "configs_E2.csv")
est1 = np.load(OUT / "est_E1.npy"); est2 = np.load(OUT / "est_E2.npy"); sd2 = np.load(OUT / "sd_E2.npy"); nav = np.load(OUT / "n_avail.npy")
meta = load_json("index_meta.json"); regions = meta["regions"]; M_EVAL0 = meta["M_EVAL0"]; IDX = np.load(OUT / "index_asof.npy")
rid = {r: i for i, r in enumerate(regions)}
lnP = T["ln_price"].to_numpy(np.float32)
VAL = (T["period"] == "val").to_numpy(); TEST = (T["period"] == "test").to_numpy()
log(f"{len(T):,} cibles (val {VAL.sum():,}, test {TEST.sum():,}) ; E1 {est1.shape[1]} ; E2 {est2.shape[1]}")

def err_of(est): return (np.exp(est - lnP[:, None]) - 1).astype(np.float32)
def metrics(err, sel):
    e = err[sel]; a = np.abs(e)
    ok = ~np.isnan(a); n = ok.sum(axis=0)
    with np.errstate(all="ignore"):
        return pd.DataFrame({"mdape": np.nanmedian(a, axis=0) * 100, "mape": np.nanmean(a, axis=0) * 100,
                             "biais": np.nanmedian(e, axis=0) * 100, "biais_moy": np.nanmean(e, axis=0) * 100,
                             "p5": np.nanmean(a <= 0.05, axis=0) * 100, "p10": np.nanmean(a <= 0.10, axis=0) * 100,
                             "p20": np.nanmean(a <= 0.20, axis=0) * 100, "couverture": n / sel.sum() * 100, "n": n})
err1 = err_of(est1); err2 = err_of(est2)
del est1
res = {}

# ---------------------------------------------------------------- densité & strates
T["dens"] = pd.cut(T["n_2km_12m"], [-1, 4, 19, 49, 199, 10 ** 9], labels=["< 5", "5-19", "20-49", "50-199", "≥ 200"])
T["prix_q"] = pd.qcut(T["amount"], 5, labels=["Q1 (bas)", "Q2", "Q3", "Q4", "Q5 (haut)"])
T["trimestre"] = pd.to_datetime(T["date"]).dt.to_period("Q").astype(str)
T["age_cl"] = pd.cut(T["age"], [-1, 3, 15, 40, 70, 1000], labels=["≤ 3 ans", "4-15", "16-40", "41-70", "> 70"])
T["cluster"] = pd.factorize(T["code_mun"].astype(str) + "|" + T["trimestre"])[0]

# ---------------------------------------------------------------- E1 : surface Erreur(k, d, t)
mv = metrics(err1, VAL).add_suffix("_val"); mt = metrics(err1, TEST).add_suffix("_test")
surf = pd.concat([E1, mv, mt], axis=1)
surf.to_csv(TAB / "surface_E1.csv", index=False)
ok = surf["couverture_val"] >= 97
best_e1 = surf[ok].sort_values("mdape_val").iloc[0]
res["E1_best_val"] = best_e1.to_dict()
log(f"E1 meilleure (val, couverture ≥ 97 %) : {best_e1[['sel','d','T','k','adj']].to_dict()} → val {best_e1.mdape_val:.2f} %, test {best_e1.mdape_test:.2f} %")
for (s, a), g in surf[ok].groupby(["sel", "adj"]):
    b = g.sort_values("mdape_val").iloc[0]
    res.setdefault("E1_best_by_sel_adj", {})[f"{s}|{a}"] = dict(k=int(b.k), d=float(b.d), T=int(b["T"]), mdape_val=float(b.mdape_val), mdape_test=float(b.mdape_test), couverture_test=float(b.couverture_test))

# figure : MdAPE vs k selon d (T = 12, geo, hedo, test)
fig, axes = plt.subplots(1, 2, figsize=(9, 3.4), sharey=True)
for ax, a in zip(axes, ["none", "hedo"]):
    for j, d in enumerate([0.5, 1, 2, 5, 10, 25]):
        g = surf[(surf.sel == "geo") & (surf["T"] == 12) & (surf.adj == a) & (surf.d == d)].sort_values("k")
        ax.plot(g.k, g.mdape_test, color=PAL[j], lw=2, marker="o", ms=4, label=f"d ≤ {d} km")
    ax.set_xscale("log"); ax.set_xticks([1, 3, 5, 10, 20, 50]); ax.set_xticklabels([1, 3, 5, 10, 20, 50])
    ax.set_xlabel("nombre de comparables k"); ax.set_title({"none": "sans ajustement (temps seulement)", "hedo": "grille hédonique"}[a])
axes[0].set_ylabel("MdAPE, test 2025-26 (%)"); axes[1].legend(ncol=2, fontsize=8)
savefig(fig, "f1_mdape_vs_k")
# figure : carte d × T au meilleur k (geo, hedo)
fig, axes = plt.subplots(1, 2, figsize=(8, 3.2))
for ax, a in zip(axes, ["none", "hedo"]):
    g = surf[(surf.sel == "geo") & (surf.adj == a)].groupby(["d", "T"])["mdape_test"].min().unstack("T")
    im = ax.imshow(g.to_numpy(), cmap=CMAP, aspect="auto")
    ax.set_xticks(range(len(g.columns))); ax.set_xticklabels([f"{t} mois" for t in g.columns]); ax.set_yticks(range(len(g.index))); ax.set_yticklabels([f"{d} km" for d in g.index])
    for i in range(g.shape[0]):
        for j in range(g.shape[1]):
            v = g.iloc[i, j]; ax.text(j, i, f"{v:.1f}", ha="center", va="center", fontsize=7.5, color="white" if v < np.nanpercentile(g.to_numpy(), 40) else "#0b0b0b")
    ax.grid(False); ax.set_title({"none": "sans ajustement", "hedo": "grille hédonique"}[a] + " — MdAPE au meilleur k")
savefig(fig, "f2_surface_d_T")

# ---------------------------------------------------------------- E2 : facteurs
mv2 = metrics(err2, VAL).add_suffix("_val"); mt2 = metrics(err2, TEST).add_suffix("_test")
fac = pd.concat([E2, mv2, mt2], axis=1); fac.to_csv(TAB / "facteurs_E2.csv", index=False)
anchor = (fac.k == 5) & (fac.d == 5) & (fac["T"] == 12)
base = anchor & (fac.sel == "geo") & (fac.w == "equal") & (fac.recon == "median")
res["E2_adjust_anchor"] = fac[base].set_index("adj")[["mdape_val", "mdape_test", "mape_test", "biais_test", "p10_test", "couverture_test"]].to_dict("index")
res["E2_sel_anchor"] = fac[anchor & (fac.adj == "hedo") & (fac.w == "equal") & (fac.recon == "median")].set_index("sel")[["mdape_val", "mdape_test", "p10_test"]].to_dict("index")
res["E2_wr_anchor"] = {f"{r.w}|{r.recon}": dict(mdape_val=r.mdape_val, mdape_test=r.mdape_test, p10_test=r.p10_test)
                       for r in fac[anchor & (fac.adj == "hedo") & (fac.sel == "geo")].itertuples()}
# effets marginaux moyens (moyenne de MdAPE test sur toutes les autres modalités)
marg = {}
for f in ["sel", "adj", "w", "recon"]:
    marg[f] = fac.groupby(f)["mdape_test"].mean().round(2).to_dict()
res["E2_marginal"] = marg
anc_tab = fac.groupby(["k", "d", "T"])["mdape_test"].min().round(2).to_dict()
res["E2_anchor_best"] = {f"{k}|{d}|{t}": v for (k, d, t), v in anc_tab.items()}
# figure : ajustements à l'ancre (val et test)
fig, ax = plt.subplots(figsize=(7, 3.2))
g = fac[base].set_index("adj").loc[["raw", "none", "conv", "hedo", "hedo_reg", "local", "gbm", "role"]]
lab = {"raw": "aucun", "none": "temps seul", "conv": "grille conv.", "hedo": "hédonique", "hedo_reg": "hédo. régionale", "local": "locale", "gbm": "GBM", "role": "ratio au rôle"}
x = np.arange(len(g)); ax.bar(x - 0.2, g.mdape_val, 0.38, color=PAL[0], label="validation 2024"); ax.bar(x + 0.2, g.mdape_test, 0.38, color=PAL[1], label="test 2025-26")
for i, (v1, v2) in enumerate(zip(g.mdape_val, g.mdape_test)):
    ax.text(i - 0.2, v1 + 0.2, f"{v1:.1f}", ha="center", fontsize=7); ax.text(i + 0.2, v2 + 0.2, f"{v2:.1f}", ha="center", fontsize=7)
ax.set_xticks(x); ax.set_xticklabels([lab[a] for a in g.index]); ax.set_ylabel("MdAPE (%)"); ax.legend(); ax.set_title("Type d'ajustement — 5 comparables les plus proches, ≤ 5 km, ≤ 12 mois, médiane")
savefig(fig, "f3_ajustements")

# ---------------------------------------------------------------- modèles de référence (M10, M11) et compétition
def level_shift(m_asof):
    mt_ = T["month"].to_numpy(); r = T["region"].map(rid).to_numpy()
    return IDX[mt_ - M_EVAL0, r, mt_ - 1] - IDX[mt_ - M_EVAL0, r, m_asof - 1]
Sx = S.loc[T["idx"]]
XF = np.column_stack([Sx.ln_floor, Sx.ln_land, Sx.age, Sx.age ** 2 / 100, Sx.lien == 1, Sx.lien >= 2, Sx.genre == 1, Sx.genre == 2, Sx.genre == 3, Sx.etages, Sx.neuf]).astype(float)
M10 = np.full(len(T), np.nan); M11 = np.full(len(T), np.nan); M11b = np.full(len(T), np.nan)
for a in load_json("asof.json")["asof"]:
    mask = (T["asof"] == a).to_numpy(); period = str(a)
    if not mask.any(): continue
    mod = load_json(f"models_{period}.json"); beta = np.array([mod["beta_global"][f] for f in mod["feats"]])
    shift = level_shift(mod["m_asof"])
    fe = (T["region"].map(mod["region_fe"]).fillna(0) + T["mrc"].map(mod["mrc_eff"]).fillna(0) + T["code_mun"].map(mod["mun_eff"]).fillna(0)
          + (T["code_mun"].astype(str) + "|" + T["unite_voisinage"].astype(str)).map(mod["uv_eff"]).fillna(0)).to_numpy()
    M10[mask] = (XF @ beta + fe + shift)[mask]
    M11[mask] = (np.load(OUT / f"gbm_role_pred_{period}.npy")[T["idx"].to_numpy()] + shift)[mask]
    M11b[mask] = (np.load(OUT / f"gbm_pred_{period}.npy")[T["idx"].to_numpy()] + shift)[mask]

# M10 : recentrage médian calibré sur la validation 2024 (la MCO vise la moyenne du log, distribution asymétrique)
M10 = M10 + np.nanmedian((lnP - M10)[VAL])
def e2col(sel, k, d, Tm, adj, w, recon):
    m = (E2.sel == sel) & (E2.k == k) & (E2.d == d) & (E2["T"] == Tm) & (E2.adj == adj) & (E2.w == w) & (E2.recon == recon)
    return int(E2.index[m][0])
def e1col(sel, k, d, Tm, adj):
    m = (E1.sel == sel) & (E1.k == k) & (E1.d == d) & (E1["T"] == Tm) & (E1.adj == adj)
    return int(E1.index[m][0])
MODELS = {
    "M1 médiane des 3 plus proches (temps seul)": ("E2", e2col("geo", 3, 2, 6, "none", "equal", "median")),
    "M2 moyenne des 5 plus proches (temps seul)": ("E2", e2col("geo", 5, 5, 12, "none", "equal", "mean")),
    "M3 sélection par similarité physique": ("E2", e2col("sim", 5, 5, 12, "none", "equal", "median")),
    "M4 pondération par distance": ("E2", e2col("geo", 5, 5, 12, "none", "dist", "mean")),
    "M5 pondération par récence": ("E2", e2col("geo", 5, 5, 12, "none", "rec", "mean")),
    "M6 grille d'ajustement conventionnelle": ("E2", e2col("geo", 5, 5, 12, "conv", "equal", "median")),
    "M7 grille d'ajustement hédonique": ("E2", e2col("geo", 5, 5, 12, "hedo", "equal", "median")),
    "M8a grille locale (régression sur le voisinage)": ("E2", e2col("geo", 5, 5, 12, "local", "equal", "median")),
    "M8b ajustements non linéaires (GBM)": ("E2", e2col("geo", 5, 5, 12, "gbm", "equal", "median")),
    "M13 ratio au rôle (5 plus proches, même rôle)": ("E2", e2col("geo", 5, 5, 12, "role", "equal", "median")),
    "M10 hédonique globale directe": ("X", M10),
    "M11a gradient boosting (sans rôle)": ("X", M11b),
    "M11b gradient boosting (avec rôle)": ("X", M11),
}
def est_of(spec):
    kind, c = spec
    return est2[:, c].astype(float) if kind == "E2" else c
# M0 : meilleure configuration empirique choisie sur la validation (E1 ∪ E2, couverture ≥ 90 %)
cand_all = pd.concat([surf.assign(src="E1"), fac.assign(src="E2")], ignore_index=True)
cand_ok = cand_all[cand_all.couverture_val >= 97].sort_values("mdape_val")
b0 = cand_ok.iloc[0]
M0_est = (est2[:, int(b0.cfg)] if b0.src == "E2" else np.load(OUT / "est_E1.npy", mmap_mode="r")[:, int(b0.cfg)]).astype(float)
res["M0"] = b0[["src", "sel", "d", "T", "k", "adj", "w", "recon", "mdape_val", "mdape_test", "couverture_test"]].to_dict()
res["M0_sans_contrainte"] = cand_all[cand_all.couverture_val >= 50].sort_values("mdape_val").iloc[0][["src", "sel", "d", "T", "k", "adj", "w", "recon", "mdape_val", "mdape_test", "couverture_val", "couverture_test"]].to_dict()
MODELS["M0 meilleure configuration empirique (choisie en 2024, couverture ≥ 97 %)"] = ("X", M0_est)
# M9 : réconciliation optimisée = combinaison convexe (NNLS sur la validation) d'estimations de comparables
stack_cols = {a: e2col("geo", 5, 5, 12, a, "equal", "median") for a in ["none", "conv", "hedo", "hedo_reg", "local", "gbm", "role"]}
def stack(names, mats, fit_mask):
    """Combinaison convexe apprise par NNLS sur la validation ; si une composante manque pour une cible,
    les poids des composantes disponibles sont renormalisés."""
    A = np.column_stack(mats); okr = ~np.isnan(A).any(axis=1)
    w, _ = nnls(A[okr & fit_mask], lnP[okr & fit_mask].astype(float)); w = w / w.sum()
    W = np.where(np.isnan(A), 0.0, w[None, :]); den = W.sum(axis=1)
    out = np.nansum(A * W, axis=1) / np.where(den > 0, den, np.nan)
    return out, dict(zip(names, w.round(3)))
M9_est, w9 = stack(list(stack_cols), [est2[:, c].astype(float) for c in stack_cols.values()], VAL)
res["M9_weights"] = w9; MODELS["M9 réconciliation optimisée (poids appris en 2024)"] = ("X", M9_est)
M12_est, w12 = stack(["gbm_role", "comps_gbm", "comps_hedo"], [M11, est2[:, stack_cols["gbm"]].astype(float), est2[:, stack_cols["hedo"]].astype(float)], VAL)
res["M12_weights"] = w12; MODELS["M12 hybride comparables + GBM (poids appris en 2024)"] = ("X", M12_est)

rows = []
ERR = {}
for name, spec in MODELS.items():
    e = (np.exp(est_of(spec) - lnP) - 1).astype(np.float32); ERR[name] = e
    for per, mask in [("val", VAL), ("test", TEST)]:
        m = metrics(e[:, None], mask).iloc[0].to_dict(); m.update(modele=name, periode=per); rows.append(m)
comp = pd.DataFrame(rows); comp.to_csv(TAB / "competition_modeles.csv", index=False)
pd.DataFrame({n.split(" ")[0]: est_of(s_) for n, s_ in MODELS.items()}).assign(idx=T["idx"].to_numpy()).to_parquet(OUT / "models_est.parquet", index=False)
res["competition"] = comp[comp.periode == "test"].set_index("modele")[["mdape", "mape", "biais", "p5", "p10", "p20", "couverture"]].round(2).to_dict("index")
res["competition_val"] = comp[comp.periode == "val"].set_index("modele")[["mdape", "p10", "couverture"]].round(2).to_dict("index")
log("compétition (test) :\n" + comp[comp.periode == "test"][["modele", "mdape", "mape", "biais", "p10", "couverture"]].round(2).to_string(index=False))
fig, ax = plt.subplots(figsize=(7.5, 5))
g = comp[comp.periode == "test"].sort_values("mdape", ascending=False)
cols = [PAL[0] if ("M1" in n and "M1 " not in n) or n.startswith(("M10", "M11", "M12")) else PAL[2] if n.startswith(("M0", "M9")) else PAL[1] for n in g.modele]
cols = [PAL[0] if n.startswith(("M10", "M11", "M12")) else (PAL[2] if n.startswith(("M0", "M9")) else PAL[1]) for n in g.modele]
ax.barh(g.modele, g.mdape, color=cols, height=0.65)
for i, v in enumerate(g.mdape): ax.text(v + 0.15, i, f"{v:.1f} %", va="center", fontsize=8)
ax.set_xlabel("MdAPE, test 2025-26 (%)"); ax.set_title("Compétition de modèles (orange : comparables ; vert : optimisés ; bleu : références externes)")
savefig(fig, "f4_competition")

# ---------------------------------------------------------------- bootstrap par grappes (municipalité × trimestre) des différences appariées
rng = np.random.default_rng(7)
def paired(eA, eB, mask, B=400):
    a = np.abs(eA[mask]); b = np.abs(eB[mask]); cl = T["cluster"].to_numpy()[mask]
    okp = ~np.isnan(a) & ~np.isnan(b); a, b, cl = a[okp], b[okp], cl[okp]
    d = a - b
    cl_ids, inv = np.unique(cl, return_inverse=True); nc = len(cl_ids)
    sums = np.bincount(inv, d, nc); cnts = np.bincount(inv, minlength=nc)
    boots_mean = []; boots_med = []
    order = np.argsort(inv); d_sorted = d[order]; starts = np.searchsorted(inv[order], np.arange(nc)); ends = np.append(starts[1:], len(d))
    for _ in range(B):
        pick = rng.integers(0, nc, nc)
        boots_mean.append(sums[pick].sum() / cnts[pick].sum())
        if _ < 150:
            idxs = np.concatenate([np.arange(starts[p], ends[p]) for p in pick])
            boots_med.append(np.median(a[order][idxs]) - np.median(b[order][idxs]))
    return dict(n=int(len(d)), diff_mape=float(d.mean() * 100), ci_mape=[float(np.percentile(boots_mean, 2.5) * 100), float(np.percentile(boots_mean, 97.5) * 100)],
                diff_mdape=float((np.median(a) - np.median(b)) * 100), ci_mdape=[float(np.percentile(boots_med, 2.5) * 100), float(np.percentile(boots_med, 97.5) * 100)],
                part_A_meilleur=float((a < b).mean() * 100))
PAIRS = [
    ("M7 grille d'ajustement hédonique", "M2 moyenne des 5 plus proches (temps seul)"),
    ("M7 grille d'ajustement hédonique", "M6 grille d'ajustement conventionnelle"),
    ("M8a grille locale (régression sur le voisinage)", "M7 grille d'ajustement hédonique"),
    ("M8b ajustements non linéaires (GBM)", "M7 grille d'ajustement hédonique"),
    ("M11b gradient boosting (avec rôle)", "M0 meilleure configuration empirique (choisie en 2024, couverture ≥ 97 %)"),
    ("M11a gradient boosting (sans rôle)", "M0 meilleure configuration empirique (choisie en 2024, couverture ≥ 97 %)"),
    ("M12 hybride comparables + GBM (poids appris en 2024)", "M11b gradient boosting (avec rôle)"),
    ("M13 ratio au rôle (5 plus proches, même rôle)", "M7 grille d'ajustement hédonique"),
    ("M3 sélection par similarité physique", "M2 moyenne des 5 plus proches (temps seul)"),
    ("M9 réconciliation optimisée (poids appris en 2024)", "M7 grille d'ajustement hédonique"),
    ("M1 médiane des 3 plus proches (temps seul)", "M2 moyenne des 5 plus proches (temps seul)"),
]
res["paires"] = {f"{a} vs {b}": paired(ERR[a], ERR[b], TEST) for a, b in PAIRS}
pd.DataFrame([dict(A=a, B=b, **res["paires"][f"{a} vs {b}"]) for a, b in PAIRS]).to_csv(TAB / "paires_bootstrap.csv", index=False)

# ---------------------------------------------------------------- strates
def strata(err, by):
    out = []
    for lev, g in T[TEST].groupby(by, observed=True):
        m = np.zeros(len(T), bool); m[g.index] = True
        r = metrics(err[:, None], m).iloc[0].to_dict(); r[by] = str(lev); out.append(r)
    return pd.DataFrame(out)
KEY = ["M2 moyenne des 5 plus proches (temps seul)", "M7 grille d'ajustement hédonique", "M8a grille locale (régression sur le voisinage)",
       "M13 ratio au rôle (5 plus proches, même rôle)", "M0 meilleure configuration empirique (choisie en 2024, couverture ≥ 97 %)", "M10 hédonique globale directe", "M11b gradient boosting (avec rôle)", "M12 hybride comparables + GBM (poids appris en 2024)"]
for by in ["region_nom", "dens", "prix_q", "trimestre", "age_cl", "lien"]:
    tabs = []
    for name in KEY:
        s = strata(ERR[name], by); s["modele"] = name.split(" ")[0]; tabs.append(s)
    st = pd.concat(tabs); st.to_csv(TAB / f"strates_{by}.csv", index=False)
    res[f"strates_{by}"] = st.pivot(index=by, columns="modele", values="mdape").round(2).to_dict("index")
    res[f"strates_{by}_n"] = st[st.modele == "M7"].set_index(by)["n"].to_dict()
# figure : densité
fig, ax = plt.subplots(figsize=(6.5, 3.4))
st = pd.read_csv(TAB / "strates_dens.csv")
order = ["< 5", "5-19", "20-49", "50-199", "≥ 200"]
for j, (name, lab_) in enumerate([("M2", "M2 5 plus proches"), ("M7", "M7 grille hédonique"), ("M0", "M0 meilleure config."), ("M11b", "M11b GBM"), ("M12", "M12 hybride")]):
    g = st[st.modele == name].set_index("dens").loc[order]
    ax.plot(order, g.mdape, color=PAL[j], lw=2, marker="o", ms=4, label=lab_)
ax.set_xlabel("ventes disponibles dans un rayon de 2 km au cours des 12 mois précédents"); ax.set_ylabel("MdAPE, test (%)"); ax.legend(fontsize=8); ax.set_title("Précision selon la liquidité du sous-marché")
savefig(fig, "f5_densite")

# ---------------------------------------------------------------- hypothèses
H = {}
# H1 : proche mais dissemblable vs loin mais semblable (k = 5, T = 12, sans ajustement)
H["H1"] = {"geo_d1": surf.loc[e1col("geo", 5, 1, 12, "none"), ["mdape_test", "couverture_test"]].to_dict(),
           "sim_d10": surf.loc[e1col("sim", 5, 10, 12, "none"), ["mdape_test", "couverture_test"]].to_dict(),
           "geo_d10": surf.loc[e1col("geo", 5, 10, 12, "none"), ["mdape_test", "couverture_test"]].to_dict(),
           "mix_d10": surf.loc[e1col("mix", 5, 10, 12, "none"), ["mdape_test", "couverture_test"]].to_dict(),
           "sim_d10_hedo": surf.loc[e1col("sim", 5, 10, 12, "hedo"), ["mdape_test", "couverture_test"]].to_dict(),
           "geo_d10_hedo": surf.loc[e1col("geo", 5, 10, 12, "hedo"), ["mdape_test", "couverture_test"]].to_dict(),
           "mix_d10_hedo": surf.loc[e1col("mix", 5, 10, 12, "hedo"), ["mdape_test", "couverture_test"]].to_dict()}
# H2 : k optimal selon la densité (geo, d = 10, T = 12, hedo et none)
h2 = {}
for a in ["none", "hedo"]:
    for lev in order:
        m = TEST & (T["dens"] == lev).to_numpy()
        cols_ = [e1col("geo", k, 10, 12, a) for k in [1, 3, 5, 7, 10, 15, 20, 30, 50]]
        md = np.nanmedian(np.abs(err1[m][:, cols_]), axis=0) * 100
        h2[f"{a}|{lev}"] = dict(zip([1, 3, 5, 7, 10, 15, 20, 30, 50], md.round(2)))
H["H2"] = h2
fig, axes = plt.subplots(1, 2, figsize=(9, 3.3), sharey=True)
for ax, a in zip(axes, ["none", "hedo"]):
    for j, lev in enumerate(order):
        v = h2[f"{a}|{lev}"]; ax.plot(list(v), list(v.values()), color=PAL[j], lw=2, marker="o", ms=3.5, label=f"densité {lev}")
    ax.set_xscale("log"); ax.set_xticks([1, 3, 5, 10, 20, 50]); ax.set_xticklabels([1, 3, 5, 10, 20, 50]); ax.set_xlabel("k"); ax.set_title({"none": "temps seul", "hedo": "grille hédonique"}[a] + " — d ≤ 10 km, 12 mois")
axes[0].set_ylabel("MdAPE, test (%)"); axes[1].legend(fontsize=8)
savefig(fig, "f6_k_par_densite")
# H3 : local vs hédonique fixe selon la densité, aux trois ancres
h3 = {}
for (k, d, Tm) in [(3, 2, 6), (5, 5, 12), (10, 10, 24)]:
    for lev in order:
        m = TEST & (T["dens"] == lev).to_numpy()
        h3[f"{k}|{d}|{Tm}|{lev}"] = {a: float(np.nanmedian(np.abs(err2[m][:, e2col("geo", k, d, Tm, a, "equal", "median")])) * 100) for a in ["conv", "hedo", "hedo_reg", "local", "gbm"]}
H["H3"] = h3
# H4 : pondération selon la densité (ancre 5/5/12, hedo)
h4 = {}
for lev in order:
    m = TEST & (T["dens"] == lev).to_numpy()
    h4[lev] = {f"{w}|{r}": float(np.nanmedian(np.abs(err2[m][:, e2col("geo", 5, 5, 12, "hedo", w, r)])) * 100) for (w, r) in [("equal", "median"), ("equal", "mean"), ("dist", "mean"), ("rec", "mean"), ("sim", "mean"), ("all", "mean"), ("all", "median")]}
H["H4"] = h4
# H5 : méthodes simples vs ML selon l'homogénéité (dispersion des prix ajustés des 10 voisins à 10 km / 24 mois)
het = sd2[:, e2col("geo", 10, 10, 24, "hedo", "equal", "median")].astype(float)
T["heterog"] = pd.qcut(pd.Series(np.where(np.isnan(het), np.nanmax(het), het)), 3, labels=["homogène", "moyen", "hétérogène"])
h5 = {}
for lev in ["homogène", "moyen", "hétérogène"]:
    m = TEST & (T["heterog"] == lev).to_numpy()
    h5[lev] = {n.split(" ")[0]: float(np.nanmedian(np.abs(ERR[n][m])) * 100) for n in KEY}
    h5[lev]["sd_median"] = float(np.nanmedian(het[m]))
H["H5"] = h5
res["hypotheses"] = H

# ---------------------------------------------------------------- intervalles d'incertitude (couverture)
intervals = {}
for a in ["none", "hedo", "gbm"]:
    for (k, d, Tm) in [(5, 5, 12), (10, 10, 24)]:
        c = e2col("geo", k, d, Tm, a, "equal", "median"); e = est2[:, c].astype(float); s = sd2[:, c].astype(float)
        okv = VAL & ~np.isnan(e) & ~np.isnan(s) & (s > 0)
        z_needed = np.abs(lnP[okv] - e[okv]) / s[okv]; z80 = float(np.percentile(z_needed, 80))
        okt = TEST & ~np.isnan(e) & ~np.isnan(s) & (s > 0)
        cov_nominal = float(np.mean(np.abs(lnP[okt] - e[okt]) <= 1.2816 * s[okt]) * 100)
        cov_cal = float(np.mean(np.abs(lnP[okt] - e[okt]) <= z80 * s[okt]) * 100)
        width = float(np.median(2 * z80 * s[okt]) * 100)
        intervals[f"{a}|{k}|{d}|{Tm}"] = dict(z80_val=z80, couverture_nominale_80=cov_nominal, couverture_calibree=cov_cal, largeur_mediane_pct=width, n_test=int(okt.sum()))
res["intervalles"] = intervals

# ---------------------------------------------------------------- règles conditionnelles (système de recommandation)
cond = {}; est_cond = np.full(len(T), np.nan)
for lev in order:
    mv_ = VAL & (T["dens"] == lev).to_numpy(); mt_ = TEST & (T["dens"] == lev).to_numpy()
    m1 = metrics(err1, mv_); m2 = metrics(err2, mv_)
    allm = pd.concat([pd.concat([E1, m1], axis=1).assign(src="E1"), pd.concat([E2, m2], axis=1).assign(src="E2")], ignore_index=True)
    b = allm[allm.couverture >= 95].sort_values("mdape").iloc[0]
    e = (err1 if b.src == "E1" else err2)[:, int(b.cfg)]
    est_cond[mt_] = e[mt_]
    cond[lev] = dict(b[["src", "sel", "d", "T", "k", "adj", "w", "recon", "mdape", "couverture"]].to_dict(), mdape_test=float(np.nanmedian(np.abs(e[mt_])) * 100), n_test=int(mt_.sum()))
res["regles_conditionnelles"] = cond
res["regles_conditionnelles_global_test"] = dict(mdape=float(np.nanmedian(np.abs(est_cond[TEST])) * 100), p10=float(np.nanmean(np.abs(est_cond[TEST]) <= 0.1) * 100), couverture=float(np.mean(~np.isnan(est_cond[TEST])) * 100))
pd.DataFrame(cond).T.to_csv(TAB / "regles_conditionnelles.csv")

# ---------------------------------------------------------------- stabilité temporelle, indice, divers
tq = []
for q, g in T[TEST].groupby("trimestre"):
    m = np.zeros(len(T), bool); m[g.index] = True
    tq.append(dict(trimestre=q, n=int(m.sum()), **{n.split(" ")[0]: float(np.nanmedian(np.abs(ERR[n][m])) * 100) for n in KEY},
                   **{"biais_" + n.split(" ")[0]: float(np.nanmedian(ERR[n][m]) * 100) for n in KEY}))
pd.DataFrame(tq).to_csv(TAB / "stabilite_trimestres.csv", index=False); res["stabilite"] = tq
fig, ax = plt.subplots(figsize=(7, 3.2))
idxf = pd.read_csv(OUT / "index_full.csv")
for j, (r, lab_) in enumerate([("QC", "Québec"), ("06", "Montréal"), ("03", "Capitale-Nationale"), ("07", "Outaouais"), ("05", "Estrie")]):
    g = idxf[idxf.region.astype(str).str.zfill(2) == r] if r != "QC" else idxf[idxf.region == "QC"]
    if len(g): ax.plot(pd.to_datetime(g.ym), np.exp(g.log_idx) * 100, color=PAL[j], lw=2 if r == "QC" else 1.4, label=lab_)
ax.set_ylabel("indice hédonique (2021-01 = 100)"); ax.legend(fontsize=8); ax.set_title("Indice de prix mensuel (plein échantillon, descriptif)")
savefig(fig, "f7_indice")
res["dataset"] = load_json("dataset_stats.json")
desc = {}
for per in ["train", "val", "test"]:
    d = S[S["period"] == per]
    desc[per] = dict(n=int(len(d)), prix_med=float(d.amount.median()), prix_moy=float(d.amount.mean()), aire_med=float(d.floor_area.median()),
                     terrain_med=float(d.land_area.median()), age_med=float(d.age.median()), detache=float((d.lien == 0).mean() * 100),
                     etages_moy=float(d.etages.mean()), ratio_role_med=float(d.ratio_role.median()))
for per, m in [("val", VAL), ("test", TEST)]:
    desc[per]["dens_med"] = float(np.median(T[m].n_2km_12m)); desc[per]["n_cand_med"] = float(np.median(T[m].n_cand))
res["descriptives"] = desc
res["n_targets"] = dict(val=int(VAL.sum()), test=int(TEST.sum()))
res["densite_distribution"] = T[TEST]["dens"].value_counts().sort_index().to_dict()
res["n_cand_median"] = dict(val=float(np.median(T[VAL].n_cand)), test=float(np.median(T[TEST].n_cand)))
save_json(res, "results.json")
log("analyse terminée")