#!/usr/bin/env python3
"""Chapitre 19 — modèles spatiaux et temporels sur le cas 1 (Gatineau, DONNÉES SYNTHÉTIQUES, graine 20261031).

Reproduit tous les chiffres cités dans manuscript/chapters/ch19.tex :
  1. matrices de poids (k plus proches voisins, distance inverse, spatio-temporelle de Dubé–Legros) ;
  2. I de Moran des résidus du modèle hédonique de référence (traite.cas1.hedonique) avec et sans indicatrices de secteur ;
  3. corrélogramme (I de Moran par bande de distance) ;
  4. tests du multiplicateur de Lagrange (LM-erreur, LM-décalage et versions robustes) ;
  5. MCO / SAR / SEM / SDM par maximum de vraisemblance concentrée, effets directs et indirects ;
  6. régression géographiquement pondérée (noyau gaussien fixe, bande passante par validation croisée) ;
  7. effets fixes de quartier (secteurs, grilles) et de temps (tendance, trimestres) ;
  8. variogramme empirique des résidus, modèle exponentiel et krigeage ordinaire ;
  9. décalage spatio-temporel (ventes antérieures seulement) ;
 10. validation croisée et indications pour le sujet ; sensibilité au choix de la matrice de poids.

Contrôle croisé (non requis pour l'exécution) : spreg (ML_Lag, ML_Error, LMtests), esda (Moran) et mgwr (GWR) donnent les
mêmes estimations à 1e-6 près ; voir --controle.
Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch19_spatial.py [--controle]"""
from __future__ import annotations
import sys, pathlib, warnings
import numpy as np, pandas as pd
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
from traite import cas1
from traite import hedonique as h

ROOT = pathlib.Path(__file__).resolve().parents[3]
K_VOISINS = 8            # matrice de référence : 8 plus proches voisins, standardisée en ligne
PORTEE_KM = 2.0          # matrice de distance inverse : portée 2 km
FENETRE_MOIS = 12        # matrice spatio-temporelle : ventes antérieures de moins de 12 mois
BANDES_CORRELOGRAMME = [0, 0.5, 1, 1.5, 2, 3, 4, 6, 8, 10, 12, 15]
BORNES_VARIOGRAMME = np.arange(0, 16.1, 1.0)
VOISINS_KRIGEAGE = 40
# primes de secteur du générateur (data/synthetic/README.md : log du rapport des médianes APCIQ T3 2026 à Hull ; Masson-Angers − 0,02)
PRIMES_DGP = {"Hull": 0.0, "Gatineau": float(np.log(464900 / 475000)), "Aylmer": float(np.log(530000 / 475000)),
              "Buckingham": float(np.log(395000 / 475000)), "Masson-Angers": float(np.log(395000 / 475000)) - 0.02}

_CACHE: dict = {}

# ----------------------------------------------------------------------------- base

def charger_base() -> dict:
    """Modèle de référence (traite.cas1.hedonique), spécification « sans secteur », coordonnées en km et sujet."""
    d = cas1.charger(); v = d["ventes"]; s = d["sujet"]
    m = cas1.hedonique(v)
    Xc = m["X"]; y = m["y"]
    Xr = Xc.drop(columns=[c for c in Xc.columns if c.startswith("sect_")])           # sans indicatrices de secteur
    Xn = Xr.drop(columns=["dist_centre"])                                              # ni secteur ni gradient de distance
    lat0, lon0 = float(v["lat"].mean()), float(v["lon"].mean())
    xy = h.coordonnees_km(v["lat"].values, v["lon"].values, lat0, lon0)
    cs = h.coordonnees_km([s["lat"]], [s["lon"]], lat0, lon0)
    xs_r = cas1.matrice(cas1.sujet_df(s))[Xr.columns].values[0]
    mco_r = h.mco(Xr.values, y); mco_n = h.mco(Xn.values, y)
    return {"ventes": v, "sujet": s, "modele": m, "Xc": Xc, "Xr": Xr, "Xn": Xn, "y": y, "xy": xy, "cs": cs, "xs_r": xs_r,
            "mois": Xc["mois"].values, "mco_r": mco_r, "mco_n": mco_n, "lat0": lat0, "lon0": lon0}

def matrices(b: dict) -> dict:
    xy, mois = b["xy"], b["mois"]
    return {"knn": h.matrice_poids_knn(xy, K_VOISINS), "dist": h.matrice_poids_distance_inverse(xy, PORTEE_KM),
            "st": h.matrice_poids_spatio_temporelle(xy, mois, PORTEE_KM, FENETRE_MOIS)}

def _loglik_mco(res: np.ndarray) -> float:
    n = len(res); return float(-n / 2 * (np.log(2 * np.pi * res @ res / n) + 1))

def _aic_mco(res: np.ndarray, k: int) -> float:
    return 2 * (k + 1) - 2 * _loglik_mco(res)

def _aicc(res: np.ndarray, k_eff: float) -> float:
    n = len(res); rss = float(res @ res)
    return n * np.log(rss / n) + n * np.log(2 * np.pi) + n * (n + k_eff) / (n - 2 - k_eff)

# ----------------------------------------------------------------------------- 2-3. Moran et corrélogramme

def tableau_moran(b: dict, Ws: dict) -> pd.DataFrame:
    series = {"ln P (sans modèle)": b["y"], "Résidus, sans secteur ni distance": b["mco_n"]["residus"],
              "Résidus, sans indicatrices de secteur": b["mco_r"]["residus"], "Résidus, modèle de référence": b["modele"]["residus"]}
    rows = []
    for nom, x in series.items():
        row = {"serie": nom}
        for wn, W in Ws.items():
            r = h.moran_i(x, W); row[f"I_{wn}"] = r["I"]; row[f"z_{wn}"] = r["z"]; row[f"p_{wn}"] = r["p"]
        rows.append(row)
    return pd.DataFrame(rows)

def correlogramme(b: dict) -> dict:
    return {"sans_secteur": h.moran_par_distance(b["mco_r"]["residus"], b["xy"], BANDES_CORRELOGRAMME),
            "reference": h.moran_par_distance(b["modele"]["residus"], b["xy"], BANDES_CORRELOGRAMME)}

# ----------------------------------------------------------------------------- 4-5. LM, SAR/SEM/SDM, effets

def tests_lm(b: dict, W: np.ndarray) -> dict:
    return {"reference": h.tests_lm_spatial(b["y"], b["Xc"].values, W), "sans_secteur": h.tests_lm_spatial(b["y"], b["Xr"].values, W)}

def modeles_spatiaux(b: dict, W: np.ndarray) -> dict:
    y, Xr, Xc = b["y"], b["Xr"].values, b["Xc"].values
    out = {"mco_r": b["mco_r"], "sar_r": h.sar_mv(y, Xr, W), "sem_r": h.sem_mv(y, Xr, W), "sdm_r": h.sar_mv(y, Xr, W, durbin=True),
           "sar_c": h.sar_mv(y, Xc, W), "sem_c": h.sem_mv(y, Xc, W)}
    for cle in ["sar_r", "sem_r", "sdm_r", "sar_c", "sem_c"]:
        res = out[cle]["innovations"] if "innovations" in out[cle] else out[cle]["residus"]
        out[cle]["moran_residus"] = h.moran_i(res, W)["I"]
    out["mco_r"]["loglik"] = _loglik_mco(b["mco_r"]["residus"]); out["mco_r"]["aic"] = _aic_mco(b["mco_r"]["residus"], Xr.shape[1])
    out["mco_r"]["moran_residus"] = h.moran_i(b["mco_r"]["residus"], W)["I"]
    out["mco_c"] = {"beta": b["modele"]["beta"], "se": b["modele"]["coef"]["se_hc1"].values, "residus": b["modele"]["residus"],
                    "loglik": _loglik_mco(b["modele"]["residus"]), "aic": _aic_mco(b["modele"]["residus"], Xc.shape[1]),
                    "sigma2": float(b["modele"]["sigma"] ** 2), "moran_residus": h.moran_i(b["modele"]["residus"], W)["I"]}
    return out

def effets(mod: dict, W: np.ndarray, colonnes) -> pd.DataFrame:
    k = len(colonnes)
    e_sar = h.effets_sar(mod["sar_r"]["rho"], mod["sar_r"]["beta"], W)
    theta = np.concatenate([[0.0], mod["sdm_r"]["beta"][k:]])
    e_sdm = h.effets_sar(mod["sdm_r"]["rho"], mod["sdm_r"]["beta"][:k], W, theta)
    df = pd.DataFrame({"variable": colonnes, "beta_mco": mod["mco_r"]["beta"], "beta_sar": mod["sar_r"]["beta"],
                       "sar_direct": e_sar["direct"], "sar_indirect": e_sar["indirect"], "sar_total": e_sar["total"],
                       "beta_sdm": mod["sdm_r"]["beta"][:k], "theta_sdm": theta,
                       "sdm_direct": e_sdm["direct"], "sdm_indirect": e_sdm["indirect"], "sdm_total": e_sdm["total"]})
    return df.set_index("variable")

# ----------------------------------------------------------------------------- 6. GWR

def gwr_resultats(b: dict) -> dict:
    y, X, xy = b["y"], b["Xr"].values, b["xy"]
    sel = h.gwr_bande_cv(y, X, xy)
    g = h.gwr(y, X, xy, sel["bande"])
    gs = h.gwr(y, X, xy, sel["bande"], cibles=b["cs"])
    cols = list(b["Xr"].columns)
    beta = pd.DataFrame(g["beta"], columns=cols)
    quart = beta.describe(percentiles=[0.25, 0.5, 0.75]).T[["min", "25%", "50%", "75%", "max"]]
    # effet de localisation local = constante locale + coefficient local du gradient × distance (au-delà de 3 km), par secteur,
    # comparé à l'effet de localisation du modèle de référence (indicatrice + gradient) et à celui du générateur (prime − 0,006·max(d − 3, 0))
    sect = b["ventes"]["secteur"].values; dist = b["Xr"]["dist_centre"].values; j = cols.index("dist_centre")
    loc_gwr = g["beta"][:, 0] + g["beta"][:, j] * dist
    coef = b["modele"]["coef"]["coef"]
    loc_ref = np.array([coef.get(f"sect_{sc}", 0.0) for sc in sect]) + coef["dist_centre"] * dist
    loc_dgp = np.array([PRIMES_DGP[sc] for sc in sect]) - 0.006 * dist
    comp = pd.DataFrame({"gwr": loc_gwr, "reference": loc_ref, "generateur": loc_dgp, "secteur": sect}).groupby("secteur").median()
    comp["gwr"] -= float(np.median(g["beta"][sect == "Hull", 0]))          # écart à la constante locale médiane de Hull
    comp = comp.rename(columns={"gwr": "gwr_localisation_rel_hull", "reference": "reference_localisation", "generateur": "generateur_localisation"})
    comp.loc["Hull", "gwr_localisation_rel_hull"] = 0.0
    aicc_mco_r = _aicc(b["mco_r"]["residus"], X.shape[1]); aicc_mco_c = _aicc(b["modele"]["residus"], b["Xc"].shape[1])
    return {"selection": sel, "gwr": g, "beta": beta, "quartiles": quart, "comparaison_secteurs": comp, "sujet_beta": gs["beta"][0],
            "aicc_mco_sans_secteur": aicc_mco_r, "aicc_mco_reference": aicc_mco_c}

# ----------------------------------------------------------------------------- 7. effets fixes

def effets_fixes(b: dict, W: np.ndarray) -> pd.DataFrame:
    y, xy, v = b["y"], b["xy"], b["ventes"]; n = len(y)
    def ligne(nom, Xm, k_spatial, k_temps):
        r = h.mco(Xm, y); k = Xm.shape[1]
        return {"specification": nom, "parametres": k, "effets_spatiaux": k_spatial, "effets_temps": k_temps, "r2": r["r2"],
                "r2_ajuste": r["r2_ajuste"], "sigma": float(np.sqrt(r["sigma2"])), "aic": _aic_mco(r["residus"], k), "moran": h.moran_i(r["residus"], W)["I"]}
    rows = [ligne("Gradient de distance seul (sans secteur) ; tendance linéaire", b["Xr"].values, 1, 1),
            ligne("Secteurs et gradient (référence) ; tendance linéaire", b["Xc"].values, 5, 1)]
    for cell in [5.0, 3.0, 2.0]:
        gx = np.floor(xy[:, 0] / cell); gy = np.floor(xy[:, 1] / cell); cid = pd.factorize(gx * 1000 + gy)[0]
        D = pd.get_dummies(cid, drop_first=True).astype(float).values
        rows.append(ligne(f"Grille de {cell:g} km ({cid.max() + 1} cellules) ; tendance linéaire", np.column_stack([b["Xr"].values, D]), D.shape[1], 1))
    q = pd.PeriodIndex(v["date_vente"], freq="Q").astype(str)
    Dq = pd.get_dummies(q, drop_first=True).astype(float).values
    rows.append(ligne(f"Secteurs ; indicatrices trimestrielles ({Dq.shape[1]})", np.column_stack([b["Xc"].drop(columns=["mois"]).values, Dq]), 4, Dq.shape[1]))
    rt = h.mco(np.column_stack([b["Xc"].drop(columns=["mois"]).values, Dq]), y)
    df = pd.DataFrame(rows)
    df.attrs["indice_trimestriel"] = pd.Series(np.concatenate([[0.0], rt["beta"][-Dq.shape[1]:]]), index=sorted(q.unique()))
    df.attrs["tendance_mensuelle"] = float(b["modele"]["coef"].loc["mois", "coef"])
    return df

# ----------------------------------------------------------------------------- 8. variogramme et krigeage

def variogramme_krigeage(b: dict) -> dict:
    e = b["mco_r"]["residus"]; xy = b["xy"]
    vg = h.variogramme_empirique(xy, e, BORNES_VARIOGRAMME)
    vg_ref = h.variogramme_empirique(xy, b["modele"]["residus"], BORNES_VARIOGRAMME)
    p = h.ajuster_variogramme(vg["h"], vg["gamma"], vg["n"], portee_max=5.0)
    loo = h.krigeage_validation_croisee(xy, e, p, VOISINS_KRIGEAGE)
    sujet = h.krigeage_ordinaire(xy, e, b["cs"], p, VOISINS_KRIGEAGE)
    # surface pour la figure : grille régulière (km)
    gx = np.linspace(xy[:, 0].min() - 1, xy[:, 0].max() + 1, 90); gy = np.linspace(xy[:, 1].min() - 1, xy[:, 1].max() + 1, 70)
    GX, GY = np.meshgrid(gx, gy); grille = np.column_stack([GX.ravel(), GY.ravel()])
    kg = h.krigeage_ordinaire(xy, e, grille, p, VOISINS_KRIGEAGE)
    return {"variogramme": vg, "variogramme_reference": vg_ref, "parametres": p, "loo": loo, "sujet": sujet,
            "grille": {"x": gx, "y": gy, "z": kg["prediction"].reshape(GX.shape), "var": kg["variance"].reshape(GX.shape)},
            "variance_residus": float(e.var()), "sigma2_reference": float(b["modele"]["sigma"] ** 2)}

# ----------------------------------------------------------------------------- 9. spatio-temporel

def spatio_temporel(b: dict, Ws: dict) -> dict:
    y, Xr, Wst = b["y"], b["Xr"].values, Ws["st"]
    e = b["mco_r"]["residus"]
    ok = Wst.sum(1) > 0                                   # ventes ayant au moins une vente antérieure voisine
    Wy = Wst @ y
    # décalage spatio-temporel : prédéterminé (matrice triangulaire après tri par date) → MV = MCO sur [W_st y, X]
    Z = np.column_stack([Wy[ok], Xr[ok]]); r = h.mco(Z, y[ok])
    sem_st = h.sem_mv(y, Xr, Wst)
    return {"n_retenues": int(ok.sum()), "n_sans_voisin": int((~ok).sum()), "voisins_moyen": float((Wst[ok] > 0).sum(1).mean()),
            "rho_st": float(r["beta"][0]), "se_rho_st": float(r["se_hc1"][0]), "beta_st": r["beta"][1:], "r2_st": r["r2"],
            "moran_e_st": h.moran_i(e, Wst), "moran_e_knn": h.moran_i(e, Ws["knn"]), "moran_e_dist": h.moran_i(e, Ws["dist"]),
            "sem_st": sem_st, "max_abs_valeur_propre": float(np.abs(np.linalg.eigvals(Wst)).max()),
            "moran_e_apres_lag": h.moran_i(r["residus"], Wst[np.ix_(ok, ok)])["I"]}

# ----------------------------------------------------------------------------- 10. validation croisée et sujet

def validation_croisee(b: dict, mod: dict, g: dict, vk: dict, W: np.ndarray) -> pd.DataFrame:
    y = b["y"]; n = len(y)
    def loo_mco(Xm):
        r = h.mco(Xm, y); H = Xm @ np.linalg.pinv(Xm.T @ Xm) @ Xm.T; return r["residus"] / (1 - np.diag(H))
    sem = mod["sem_r"]; u = sem["residus"]; B = np.eye(n) - sem["lambda"] * W; Q = B.T @ B
    Eu = np.array([-(Q[i] @ u - Q[i, i] * u[i]) / Q[i, i] for i in range(n)])
    series = {"MCO, modèle de référence (secteurs)": loo_mco(b["Xc"].values),
              "MCO, sans indicatrices de secteur": loo_mco(b["Xr"].values),
              "MCO sans secteur + krigeage des résidus": b["mco_r"]["residus"] - vk["loo"],
              "SEM sans secteur (espérance conditionnelle)": u - Eu,
              "SAR sans secteur (prédiction conditionnelle)": y - (mod["sar_r"]["rho"] * (W @ y) + b["Xr"].values @ mod["sar_r"]["beta"]),
              "GWR sans secteur (validation croisée)": None}
    rows = []
    for nom, r in series.items():
        if r is None:
            rmse = float(np.sqrt(g["gwr"]["cv"] / n)); rows.append({"modele": nom, "rmse_log": rmse, "erreur_abs_moy_pct": np.nan, "moran_erreurs": h.moran_i(g["gwr"]["residus"], W)["I"]})
        else:
            rows.append({"modele": nom, "rmse_log": float(np.sqrt(np.mean(r ** 2))), "erreur_abs_moy_pct": float(np.mean(np.abs(np.exp(r) - 1))), "moran_erreurs": h.moran_i(r, W)["I"]})
    return pd.DataFrame(rows)

def predictions_sujet(b: dict, mod: dict, g: dict, vk: dict, W: np.ndarray) -> pd.DataFrame:
    s, xs = b["sujet"], b["xs_r"]; y = b["y"]; n = len(y)
    duan = lambda res: float(np.mean(np.exp(res)))
    ref = float(cas1.predire(b["modele"], cas1.sujet_df(s))[0])
    mco_r = float(np.exp(xs @ b["mco_r"]["beta"]) * duan(b["mco_r"]["residus"]))
    krig = float(np.exp(xs @ b["mco_r"]["beta"] + vk["sujet"]["prediction"][0]) * duan(b["mco_r"]["residus"] - vk["loo"]))
    # SEM : espérance conditionnelle de u au sujet (matrice augmentée : le sujet n'est voisin d'aucune vente)
    xy_aug = np.vstack([b["xy"], b["cs"]]); Wa = np.zeros((n + 1, n + 1)); Wa[:n, :n] = W
    Wa[n, :n] = h.matrice_poids_knn(xy_aug, K_VOISINS)[n, :n]; Wa[n] /= Wa[n].sum()
    sem = mod["sem_r"]; u_s = float(h.sem_esperance_conditionnelle(sem["lambda"], Wa, sem["residus"], [n])[0])
    sem_pred = float(np.exp(xs @ sem["beta"] + u_s) * duan(sem["innovations"]))
    # SAR : prédiction conditionnelle avec le décalage observé (moyenne des ln P des k ventes voisines)
    idx = np.nonzero(Wa[n, :n])[0]; Wy_s = float(Wa[n, :n] @ y)
    sar = mod["sar_r"]; sar_pred = float(np.exp(sar["rho"] * Wy_s + xs @ sar["beta"]) * duan(sar["residus"]))
    gwr_pred = float(np.exp(xs @ g["sujet_beta"]) * duan(g["gwr"]["residus"]))
    df = pd.DataFrame([
        {"modele": "MCO, modèle de référence (indicatrices de secteur)", "ln_correction": 0.0, "indication": ref},
        {"modele": "MCO sans indicatrices de secteur", "ln_correction": 0.0, "indication": mco_r},
        {"modele": "MCO sans secteur + krigeage du résidu au sujet", "ln_correction": float(vk["sujet"]["prediction"][0]), "indication": krig},
        {"modele": "SEM sans secteur (espérance conditionnelle du résidu)", "ln_correction": u_s, "indication": sem_pred},
        {"modele": "SAR sans secteur (décalage des 8 ventes voisines)", "ln_correction": float(sar["rho"] * Wy_s), "indication": sar_pred},
        {"modele": "GWR sans secteur (coefficients locaux au sujet)", "ln_correction": float(g["sujet_beta"][0] - b["mco_r"]["beta"][0]), "indication": gwr_pred}])
    df["ecart_log_vs_mco_sans_secteur"] = np.log(df["indication"] / mco_r)
    df.attrs["voisins_sujet"] = list(b["ventes"]["id_vente"].values[idx]); df.attrs["secteurs_voisins"] = list(b["ventes"]["secteur"].values[idx])
    df.attrs["variance_krigeage_sujet"] = float(vk["sujet"]["variance"][0])
    return df

# ----------------------------------------------------------------------------- 11. sensibilité

def sensibilite(b: dict) -> pd.DataFrame:
    y, Xr, xy, e = b["y"], b["Xr"].values, b["xy"], b["mco_r"]["residus"]
    specs = [("k plus proches voisins", k, h.matrice_poids_knn(xy, k)) for k in [4, 6, 8, 12, 16]]
    specs += [("distance inverse", d, h.matrice_poids_distance_inverse(xy, d)) for d in [2.0, 3.0, 4.0]]
    rows = []
    for fam, par, W in specs:
        mi = h.moran_i(e, W); sar = h.sar_mv(y, Xr, W); sem = h.sem_mv(y, Xr, W)
        rows.append({"famille": fam, "parametre": par, "voisins_moyen": float((W > 0).sum(1).mean()), "lignes_vides": int((W.sum(1) == 0).sum()), "I": mi["I"], "z": mi["z"],
                     "rho": sar["rho"], "se_rho": sar["se_rho"], "aic_sar": sar["aic"], "lambda": sem["lambda"], "se_lambda": sem["se_lambda"], "aic_sem": sem["aic"]})
    return pd.DataFrame(rows)

# ----------------------------------------------------------------------------- contrôle croisé facultatif

def controle_croise(b: dict, W: np.ndarray, g: dict) -> dict:
    """Compare les estimations maison à spreg/esda/mgwr (installés sur la machine de rédaction)."""
    import spreg
    from libpysal.weights import W as Wcls
    from esda.moran import Moran
    from mgwr.sel_bw import Sel_BW
    from mgwr.gwr import GWR
    y, X = b["y"], b["Xr"]
    neigh = {i: list(np.nonzero(W[i])[0]) for i in range(len(W))}; wts = {i: [W[i, j] for j in neigh[i]] for i in neigh}
    w = Wcls(neigh, wts); w.transform = "r"; Xs = X.drop(columns=["const"]).values; yy = y.reshape(-1, 1)
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        ml = spreg.ML_Lag(yy, Xs, w, method="full"); me = spreg.ML_Error(yy, Xs, w, method="full")
        ols = spreg.OLS(yy, Xs, w=w, spat_diag=True, moran=True)
        mo = Moran(b["mco_r"]["residus"], w, permutations=0)
        bw = Sel_BW(b["xy"], yy, Xs, fixed=True, kernel="gaussian").search(criterion="CV")
        gm = GWR(b["xy"], yy, Xs, bw, fixed=True, kernel="gaussian").fit()
    mine_sar = h.sar_mv(y, X.values, W); mine_sem = h.sem_mv(y, X.values, W); lm = h.tests_lm_spatial(y, X.values, W)
    return {"rho_spreg": float(ml.rho), "rho_maison": mine_sar["rho"], "lambda_spreg": float(me.lam), "lambda_maison": mine_sem["lambda"],
            "lm_err_spreg": float(ols.lm_error[0]), "lm_err_maison": lm["LM_erreur"][0], "rlm_lag_spreg": float(ols.rlm_lag[0]), "rlm_lag_maison": lm["RLM_decalage"][0],
            "moran_esda_z": float(mo.z_norm), "moran_maison_z": h.moran_i(b["mco_r"]["residus"], W)["z"],
            "bande_mgwr": float(bw), "bande_maison": g["selection"]["bande"], "aicc_mgwr": float(gm.aicc), "aicc_maison_meme_bande": h.gwr(y, X.values, b["xy"], float(bw))["aicc"]}

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

def main(verbose: bool = True, controle: bool = False) -> dict:
    if "res" in _CACHE and not controle:
        return _CACHE["res"]
    b = charger_base(); Ws = matrices(b); W = Ws["knn"]
    res = {"base": b, "W": Ws, "moran": tableau_moran(b, Ws), "correlogramme": correlogramme(b), "lm": tests_lm(b, W)}
    res["modeles"] = modeles_spatiaux(b, W)
    res["effets"] = effets(res["modeles"], W, list(b["Xr"].columns))
    res["gwr"] = gwr_resultats(b)
    res["effets_fixes"] = effets_fixes(b, W)
    res["variogramme"] = variogramme_krigeage(b)
    res["spatio_temporel"] = spatio_temporel(b, Ws)
    res["validation"] = validation_croisee(b, res["modeles"], res["gwr"], res["variogramme"], W)
    res["sujet"] = predictions_sujet(b, res["modeles"], res["gwr"], res["variogramme"], W)
    res["sensibilite"] = sensibilite(b)
    if controle:
        res["controle"] = controle_croise(b, W, res["gwr"])
    _CACHE["res"] = res
    if verbose:
        imprimer(res)
    return res

def imprimer(res: dict):
    pd.set_option("display.width", 220); pd.set_option("display.max_columns", 30)
    b = res["base"]; m = b["modele"]
    print("== Modèle de référence : R² %.3f, σ %.4f, n %d, k %d" % (m["r2"], m["sigma"], m["n"], m["k"]))
    print("== Sans secteur : R² %.3f, σ %.4f ; sans secteur ni distance : R² %.3f" % (b["mco_r"]["r2"], np.sqrt(b["mco_r"]["sigma2"]), b["mco_n"]["r2"]))
    print("== Étendue des coordonnées (km) :", b["xy"].min(0).round(1), b["xy"].max(0).round(1), "; sujet :", b["cs"].round(2))
    for wn, W in res["W"].items():
        print(f"   matrice {wn} : voisins moyens {(W > 0).sum(1).mean():.1f}, lignes vides {(W.sum(1) == 0).sum()}")
    print("== I de Moran :\n", res["moran"].round(4).to_string(index=False))
    print("== Corrélogramme (sans secteur) :", [(c["d_min"], c["d_max"], round(c["I"], 3), round(c["z"], 1), c["paires"]) for c in res["correlogramme"]["sans_secteur"]])
    print("== Corrélogramme (référence) :", [(c["d_min"], c["d_max"], round(c["I"], 3), round(c["z"], 1)) for c in res["correlogramme"]["reference"]])
    for k, v in res["lm"].items():
        print(f"== Tests LM ({k}) :", {a: (round(s, 2), round(p, 4)) for a, (s, p) in v.items()})
    mod = res["modeles"]
    print("== SAR sans secteur : ρ %.4f (%.4f), logL %.2f, AIC %.1f, I résidus %.3f" % (mod["sar_r"]["rho"], mod["sar_r"]["se_rho"], mod["sar_r"]["loglik"], mod["sar_r"]["aic"], mod["sar_r"]["moran_residus"]))
    print("== SEM sans secteur : λ %.4f (%.4f), logL %.2f, AIC %.1f, I innovations %.3f" % (mod["sem_r"]["lambda"], mod["sem_r"]["se_lambda"], mod["sem_r"]["loglik"], mod["sem_r"]["aic"], mod["sem_r"]["moran_residus"]))
    print("== SDM sans secteur : ρ %.4f (%.4f), logL %.2f, AIC %.1f, I résidus %.3f" % (mod["sdm_r"]["rho"], mod["sdm_r"]["se_rho"], mod["sdm_r"]["loglik"], mod["sdm_r"]["aic"], mod["sdm_r"]["moran_residus"]))
    print("== MCO sans secteur : logL %.2f, AIC %.1f ; MCO référence : logL %.2f, AIC %.1f" % (mod["mco_r"]["loglik"], mod["mco_r"]["aic"], mod["mco_c"]["loglik"], mod["mco_c"]["aic"]))
    print("== SAR référence : ρ %.4f (%.4f) ; SEM référence : λ %.4f (%.4f)" % (mod["sar_c"]["rho"], mod["sar_c"]["se_rho"], mod["sem_c"]["lambda"], mod["sem_c"]["se_lambda"]))
    cols = list(b["Xr"].columns)
    tab = pd.DataFrame({"MCO réf.": pd.Series(m["beta"], index=b["Xc"].columns).reindex(cols), "MCO": mod["mco_r"]["beta"], "SAR": mod["sar_r"]["beta"], "SEM": mod["sem_r"]["beta"], "SDM β": mod["sdm_r"]["beta"][:len(cols)],
                        "SDM θ": np.concatenate([[np.nan], mod["sdm_r"]["beta"][len(cols):]]),
                        "se MCO": mod["mco_r"]["se_hc1"], "se SAR": mod["sar_r"]["se"], "se SEM": mod["sem_r"]["se"], "se SDM": mod["sdm_r"]["se"][:len(cols)]}, index=cols)
    print("== Coefficients :\n", tab.round(4).to_string())
    print("== Effets directs/indirects/totaux :\n", res["effets"].round(4).to_string())
    g = res["gwr"]
    print("== GWR : bande %.3f km (CV %.4f), tr(S) %.1f, AICc %.1f (MCO sans secteur %.1f ; référence %.1f), R² %.3f, σ %.4f" % (g["selection"]["bande"], g["selection"]["cv"], g["gwr"]["trace_S"], g["gwr"]["aicc"], g["aicc_mco_sans_secteur"], g["aicc_mco_reference"], g["gwr"]["r2"], g["gwr"]["sigma"]))
    print("   quartiles des coefficients locaux :\n", g["quartiles"].round(4).to_string())
    print("   constante locale médiane par secteur (écart à Hull) vs indicatrice vs générateur :\n", g["comparaison_secteurs"].round(4).to_string())
    print("   coefficients locaux au sujet :", dict(zip(cols, g["sujet_beta"].round(4))))
    print("== Effets fixes :\n", res["effets_fixes"].round(4).to_string(index=False))
    print("   indice trimestriel (log) :", res["effets_fixes"].attrs["indice_trimestriel"].round(3).to_dict())
    vk = res["variogramme"]
    print("== Variogramme (sans secteur) : h", vk["variogramme"]["h"].round(2), "\n   γ", vk["variogramme"]["gamma"].round(5), "\n   N", vk["variogramme"]["n"])
    print("   référence γ", vk["variogramme_reference"]["gamma"].round(5))
    print("   paramètres :", {k: round(v, 5) for k, v in vk["parametres"].items()}, "; var(e) %.5f ; σ² réf %.5f" % (vk["variance_residus"], vk["sigma2_reference"]))
    print("   krigeage au sujet : %.4f (variance %.5f)" % (vk["sujet"]["prediction"][0], vk["sujet"]["variance"][0]))
    st = res["spatio_temporel"]
    print("== Spatio-temporel : %d ventes retenues (%d sans voisin antérieur), %.1f voisins en moyenne, max |valeur propre| %.1e" % (st["n_retenues"], st["n_sans_voisin"], st["voisins_moyen"], st["max_abs_valeur_propre"]))
    print("   I(e) : W_st %.3f (z %.1f) ; knn %.3f ; distance %.3f" % (st["moran_e_st"]["I"], st["moran_e_st"]["z"], st["moran_e_knn"]["I"], st["moran_e_dist"]["I"]))
    print("   ρ_st (MCO sur ventes antérieures) %.4f (%.4f), R² %.3f, I résidus après décalage %.3f ; SEM-st λ %.4f (%.4f), AIC %.1f" % (st["rho_st"], st["se_rho_st"], st["r2_st"], st["moran_e_apres_lag"], st["sem_st"]["lambda"], st["sem_st"]["se_lambda"], st["sem_st"]["aic"]))
    print("== Validation croisée :\n", res["validation"].round(4).to_string(index=False))
    print("== Sujet :\n", res["sujet"].round(4).to_string(index=False)); print("   voisins du sujet :", res["sujet"].attrs["voisins_sujet"], res["sujet"].attrs["secteurs_voisins"])
    print("== Sensibilité :\n", res["sensibilite"].round(4).to_string(index=False))
    if "controle" in res:
        print("== Contrôle croisé spreg/esda/mgwr :", {k: round(v, 5) for k, v in res["controle"].items()})

if __name__ == "__main__":
    main(controle="--controle" in sys.argv)