#!/usr/bin/env python3
"""Chapitre 15 — Études de marché et segmentation. Calculs de l'exemple intégral et des tableaux du chapitre.

Trois blocs, tous reproductibles (graine 20261031) :
  A. indicateurs de marché (mois d'inventaire, taux d'absorption, délais) à partir des chiffres APCIQ du T3 2026 cités
     avec attribution (data/synthetic/calibration.json, sources/notes-E.md) — chiffres réels, calculs de l'auteur ;
  B. segmentation des 600 ventes SYNTHÉTIQUES du cas 1 (Gatineau) : a priori par secteur, puis statistique (k-moyennes,
     classification hiérarchique de Ward, silhouette) ; test de Chow (classique et Wald robuste HC1) de stabilité des
     coefficients hédoniques entre segments, via traite.hedonique.mco ; validation croisée pooled vs segmenté ;
  C. étude de marché locatif pour le cas 3 (12 logements, Rosemont) : loyers, inoccupation, univers, rotation et offre
     neuve de la zone SCHL 9 (Enquête sur les logements locatifs, octobre 2025) confrontés à l'état locatif synthétique
     (traite.cas3.revenus) — on ne recalcule aucun RNE ici.
Usage (racine du dépôt) : PYTHONPATH=analysis/python python3 analysis/python/cases/ch15_segmentation.py"""
from __future__ import annotations
import sys, math, pathlib
import numpy as np, pandas as pd
sys.path.insert(0, str(pathlib.Path(__file__).resolve().parents[1]))
from traite import cas1, cas3
from traite import hedonique as h
from scipy import stats
from scipy.cluster.hierarchy import linkage, fcluster
from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score, silhouette_samples, adjusted_rand_score
from sklearn.model_selection import KFold

ROOT = pathlib.Path(__file__).resolve().parents[3]
GRAINE = 20261031

# =============================================================================== A. indicateurs de marché (APCIQ T3 2026)
# Source : APCIQ par le système Centris, communiqué du 14 octobre 2026 et baromètres FSMI T3 2026 (consultés le 2026-10-31).
# Ventes résidentielles du trimestre, inscriptions en vigueur ; délais de vente médians (jours) unifamiliale / copropriété / plex.
APCIQ_T3_2026 = {
    "Province de Québec":    {"ventes": 20935, "inscriptions": 42408, "delais": (45, 58, 48)},
    "RMR de Montréal":       {"ventes": 9265,  "inscriptions": 20350, "delais": (40, 61, 51)},
    "RMR de Québec":         {"ventes": 2161,  "inscriptions": 2204,  "delais": (29, 34, 30)},
    "RMR de Gatineau":       {"ventes": 1044,  "inscriptions": 2256,  "delais": (34, 36, 38)},
    "RMR de Sherbrooke":     {"ventes": 569,   "inscriptions": 825,   "delais": None},
    "RMR de Trois-Rivières": {"ventes": 372,   "inscriptions": 447,   "delais": None},
    "RMR de Saguenay":       {"ventes": 371,   "inscriptions": 498,   "delais": None},
}
# Baromètre Gatineau T3 2026 : ventes résidentielles par secteur Centris et médiane unifamiliale ($).
GATINEAU_SECTEURS_T3_2026 = {
    "Gatineau": {"ventes": 346, "mediane_unif": 464900}, "Hull": {"ventes": 174, "mediane_unif": 475000},
    "Périphérie": {"ventes": 215, "mediane_unif": 530000}, "Buckingham/Masson-Angers": {"ventes": 108, "mediane_unif": 395000},
    "Aylmer": {"ventes": 201, "mediane_unif": 530000},
}
JOURS_PAR_MOIS = 365.25 / 12

def mois_inventaire(inscriptions: float, ventes_trimestre: float) -> float:
    """MI = inscriptions en vigueur / (ventes du trimestre / 3)."""
    return inscriptions / (ventes_trimestre / 3.0)

def taux_absorption(inscriptions: float, ventes_trimestre: float) -> float:
    """Part du stock écoulée par mois = ventes mensuelles / inscriptions = 1 / MI."""
    return (ventes_trimestre / 3.0) / inscriptions

def indicateurs_apciq() -> pd.DataFrame:
    rows = []
    for terr, d in APCIQ_T3_2026.items():
        mi = mois_inventaire(d["inscriptions"], d["ventes"])
        r = {"territoire": terr, "ventes_trim": d["ventes"], "inscriptions": d["inscriptions"], "ventes_mois": d["ventes"] / 3.0,
             "mois_inventaire": mi, "taux_absorption_mensuel": taux_absorption(d["inscriptions"], d["ventes"])}
        if d["delais"]:
            r.update({"delai_unif_j": d["delais"][0], "delai_copro_j": d["delais"][1], "delai_plex_j": d["delais"][2],
                      "delai_unif_mois": d["delais"][0] / JOURS_PAR_MOIS, "ratio_stock_flux": mi / (d["delais"][0] / JOURS_PAR_MOIS)})
        rows.append(r)
    return pd.DataFrame(rows)

def parts_secteurs_gatineau(v: pd.DataFrame) -> pd.DataFrame:
    """Parts des ventes par secteur : marché réel (APCIQ T3 2026) vs jeu synthétique du cas 1 (Buckingham + Masson-Angers regroupés ;
    la périphérie n'existe pas dans le jeu synthétique)."""
    reel = pd.Series({k: d["ventes"] for k, d in GATINEAU_SECTEURS_T3_2026.items()})
    syn = v["secteur"].replace({"Buckingham": "Buckingham/Masson-Angers", "Masson-Angers": "Buckingham/Masson-Angers"}).value_counts()
    out = pd.DataFrame({"ventes_apciq_t3_2026": reel, "part_apciq": reel / reel.sum(), "ventes_synthetiques": syn.reindex(reel.index).fillna(0).astype(int)})
    out["part_synthetique"] = out["ventes_synthetiques"] / out["ventes_synthetiques"].sum()
    return out

# =============================================================================== B. segmentation (cas 1, synthétique)
def matrice_base(df: pd.DataFrame) -> pd.DataFrame:
    """Régresseurs du modèle de référence SANS les indicatrices de secteur (les segments les remplacent)."""
    X = cas1.matrice(df)
    return X[[c for c in X.columns if not c.startswith("sect_")]]

def _ssr(X: np.ndarray, y: np.ndarray) -> float:
    r = h.mco(X, y); return float(r["residus"] @ r["residus"])

def chow(v: pd.DataFrame, groupes: pd.Series, nom: str = "groupe") -> dict:
    """Test de Chow (1960) de l'égalité des coefficients entre G segments.
    Version « pentes » : le modèle contraint garde une constante propre à chaque segment (indicatrices) et des pentes communes ;
    le modèle libre estime une équation par segment. F = [(SSR_c − SSR_l)/q] / [SSR_l/(n − G k0)], q = (G − 1)(k0 − 1).
    Version « classique » : modèle contraint entièrement commun (constante incluse), q = (G − 1) k0.
    Version robuste : Wald HC1 sur les interactions (régresseurs × indicatrices) du modèle empilé (χ² à q degrés de liberté)."""
    Xb = matrice_base(v); y = np.log(v["prix_vente"].astype(float)).values
    g = pd.Series(groupes, index=v.index).astype(str); niveaux = sorted(g.unique()); G = len(niveaux); n, k0 = Xb.shape
    ref = niveaux[0]
    D = pd.get_dummies(g, prefix="g").astype(float); D = D[[f"g_{l}" for l in niveaux[1:]]]
    # contraint (constantes par segment) et classique (tout commun)
    ssr_c = _ssr(np.column_stack([Xb.values, D.values]), y)
    ssr_classique = _ssr(Xb.values, y)
    # libre : une régression par segment
    par_segment = {}; ssr_l = 0.0
    for l in niveaux:
        m = (g == l).values; r = h.mco(Xb.values[m], y[m]); ssr_l += float(r["residus"] @ r["residus"])
        par_segment[l] = {"n": int(m.sum()), "coef": pd.Series(r["beta"], index=Xb.columns), "se_hc1": pd.Series(r["se_hc1"], index=Xb.columns),
                          "sigma": float(np.sqrt(r["sigma2"])), "r2": float(r["r2"])}
    q = (G - 1) * (k0 - 1); df2 = n - G * k0
    F = ((ssr_c - ssr_l) / q) / (ssr_l / df2); p = float(stats.f.sf(F, q, df2))
    q0 = (G - 1) * k0; F0 = ((ssr_classique - ssr_l) / q0) / (ssr_l / df2); p0 = float(stats.f.sf(F0, q0, df2))
    # Wald robuste sur le modèle empilé : X_b, puis X_b(sans const) × D_g pour g ≠ ref, puis D_g (constantes)
    pentes = [c for c in Xb.columns if c != "const"]
    blocs = [Xb.values]; noms = list(Xb.columns)
    for l in niveaux[1:]:
        d = D[f"g_{l}"].values[:, None]
        blocs.append(Xb[pentes].values * d); noms += [f"{c}×{l}" for c in pentes]
        blocs.append(d); noms.append(f"const×{l}")
    Xs = np.column_stack(blocs); r = h.mco(Xs, y)
    idx = [i for i, nme in enumerate(noms) if "×" in nme and not nme.startswith("const×")]
    theta = r["beta"][idx]; V = r["V_hc1"][np.ix_(idx, idx)]
    W = float(theta @ np.linalg.solve(V, theta)); pW = float(stats.chi2.sf(W, len(idx)))
    return {"nom": nom, "niveaux": niveaux, "G": G, "n": n, "k0": k0, "ssr_contraint": ssr_c, "ssr_classique": ssr_classique, "ssr_libre": ssr_l,
            "q": q, "df2": df2, "F": float(F), "p": p, "F_classique": float(F0), "q_classique": q0, "p_classique": p0,
            "wald_hc1": W, "p_wald": pW, "q_wald": len(idx), "par_segment": par_segment,
            "sigma_contraint": math.sqrt(ssr_c / (n - k0 - G + 1)), "sigma_libre": math.sqrt(ssr_l / df2)}

def tableau_coefficients(ch: dict, variables=("ln_S", "ln_T", "age", "qualite", "garage", "mois")) -> pd.DataFrame:
    rows = []
    for l, d in ch["par_segment"].items():
        row = {"segment": l, "n": d["n"], "sigma": d["sigma"], "r2": d["r2"]}
        for vv in variables:
            row[vv] = d["coef"][vv]; row[f"se_{vv}"] = d["se_hc1"][vv]
        rows.append(row)
    return pd.DataFrame(rows)

def caracteristiques_segmentation(v: pd.DataFrame, jeu: str = "attributs+coordonnees") -> pd.DataFrame:
    """Variables de classification, centrées-réduites. jeu ∈ {attributs, coordonnees, attributs+coordonnees}."""
    xy = h.coordonnees_km(v["lat"].values, v["lon"].values)
    F = pd.DataFrame(index=v.index)
    if "attributs" in jeu:
        F["ln_S"] = np.log(v["superficie_habitable_m2"]); F["ln_T"] = np.log(v["superficie_terrain_m2"])
        F["age"] = v["age"].astype(float); F["qualite"] = v["qualite_1a5"].astype(float); F["etages"] = (v["type"] == "à étages").astype(float)
    if "coordonnees" in jeu:
        F["x_km"] = xy[:, 0]; F["y_km"] = xy[:, 1]; F["dist_centre"] = v["dist_centre_km"].astype(float)
    Z = (F - F.mean()) / F.std(ddof=0)
    return Z

def silhouette_manuelle(Z: np.ndarray, labels: np.ndarray) -> np.ndarray:
    """s(i) = (b_i − a_i)/max(a_i, b_i) — implémentation directe (vérification de sklearn dans les tests)."""
    d = np.sqrt(((Z[:, None, :] - Z[None, :, :]) ** 2).sum(-1)); s = np.zeros(len(Z))
    for i in range(len(Z)):
        own = labels == labels[i]; own[i] = False
        if own.sum() == 0: s[i] = 0.0; continue
        a = d[i, own].mean()
        b = min(d[i, labels == l].mean() for l in np.unique(labels) if l != labels[i])
        s[i] = (b - a) / max(a, b)
    return s

def classifier(Z: pd.DataFrame, k_max: int = 8) -> dict:
    """k-moyennes (20 initialisations, graine fixée) et Ward pour k = 2..k_max ; silhouette moyenne pour chaque k."""
    Zv = Z.values; L = linkage(Zv, method="ward")
    rows = []; labels_km = {}; labels_w = {}; inertie = {}
    for k in range(2, k_max + 1):
        km = KMeans(n_clusters=k, n_init=20, random_state=GRAINE).fit(Zv)
        lw = fcluster(L, t=k, criterion="maxclust")
        labels_km[k] = km.labels_; labels_w[k] = lw; inertie[k] = float(km.inertia_)
        rows.append({"k": k, "silhouette_kmeans": float(silhouette_score(Zv, km.labels_)), "silhouette_ward": float(silhouette_score(Zv, lw)),
                     "inertie_kmeans": float(km.inertia_), "ari_kmeans_ward": float(adjusted_rand_score(km.labels_, lw))})
    tab = pd.DataFrame(rows)
    k_opt = int(tab.loc[tab["silhouette_kmeans"].idxmax(), "k"])
    return {"linkage": L, "tableau": tab, "labels_kmeans": labels_km, "labels_ward": labels_w, "k_opt": k_opt, "inertie": inertie}

def profil_segments(v: pd.DataFrame, labels: np.ndarray) -> pd.DataFrame:
    d = v.assign(segment=labels)
    g = d.groupby("segment")
    out = pd.DataFrame({"n": g.size(), "superficie_med": g["superficie_habitable_m2"].median(), "age_med": g["age"].median(),
                        "qualite_moy": g["qualite_1a5"].mean(), "dist_centre_med": g["dist_centre_km"].median(),
                        "prix_med": g["prix_vente"].median(), "prix_m2_med": (d["prix_vente"] / d["superficie_habitable_m2"]).groupby(d["segment"]).median(),
                        "secteur_dominant": g["secteur"].agg(lambda s: s.value_counts().index[0]),
                        "part_secteur_dominant": g["secteur"].agg(lambda s: s.value_counts().iloc[0] / len(s))})
    return out

def validation_croisee(v: pd.DataFrame, segments: dict, n_plis: int = 10) -> pd.DataFrame:
    """RMSE (log) et erreur absolue médiane (%) hors échantillon, 10 plis : modèle de référence (pentes communes + indicatrices de
    secteur), pentes communes + indicatrices de segment, et une équation par segment (variance accrue sur les petits segments)."""
    y = np.log(v["prix_vente"].astype(float)).values; Xb = matrice_base(v)
    kf = KFold(n_splits=n_plis, shuffle=True, random_state=GRAINE)
    def ajuster_predire(Xtr, ytr, Xte):
        b = h.mco(Xtr, ytr)["beta"]; return Xte @ b
    resultats = []
    specs = {"Référence : pentes communes + indicatrices de secteur": ("dummies", v["secteur"].astype(str))}
    for nom, lab in segments.items():
        specs[f"Pentes communes + indicatrices ({nom})"] = ("dummies", pd.Series(lab, index=v.index).astype(str))
        specs[f"Une équation par segment ({nom})"] = ("separe", pd.Series(lab, index=v.index).astype(str))
    for nom, (mode, lab) in specs.items():
        err = np.zeros(len(y))
        D = pd.get_dummies(lab).astype(float).values
        for tr, te in kf.split(Xb):
            if mode == "dummies":
                Xd = np.column_stack([Xb.values, D[:, 1:]])
                err[te] = y[te] - ajuster_predire(Xd[tr], y[tr], Xd[te])
            else:
                for l in np.unique(lab):
                    mtr = tr[(lab.values[tr] == l)]; mte = te[(lab.values[te] == l)]
                    if len(mte) == 0: continue
                    if len(mtr) <= Xb.shape[1] + 2:   # segment trop petit : repli sur le modèle commun
                        err[mte] = y[mte] - ajuster_predire(Xb.values[tr], y[tr], Xb.values[mte])
                    else:
                        err[mte] = y[mte] - ajuster_predire(Xb.values[mtr], y[mtr], Xb.values[mte])
        resultats.append({"specification": nom, "rmse_log": float(np.sqrt(np.mean(err ** 2))), "mape_mediane_pct": float(100 * np.median(np.abs(np.exp(err) - 1))),
                          "erreur_abs_moy_pct": float(100 * np.mean(np.abs(np.exp(err) - 1)))})
    return pd.DataFrame(resultats)

# =============================================================================== C. marché locatif, cas 3 (SCHL zone 9 Rosemont)
# Source : SCHL, Enquête sur les logements locatifs, octobre 2025 (Rapport sur le marché locatif publié le 11 décembre 2025),
# tableaux de données RMR de Montréal : 1.1.1 (inoccupation), 1.1.2 (loyer moyen), 1.1.3 (univers), 1.1.5 (variation même échantillon),
# 1.1.6 (rotation), 1.1.9 (loyer des logements inoccupés vs occupés), 1.3.1-1.3.4 (taille d'immeuble), 3.1.7 (parc neuf).
# Lettres : indicateur de qualité SCHL (a excellent, b très bon, c bon, d passable) ; None = donnée supprimée (**).
SCHL_ZONE9 = {
    "zone": "Zone 9 — Rosemont/La Petite-Patrie",
    "loyer_2025": {"studio": 1010, "1 ch.": 1055, "2 ch.": 1348, "3 ch. +": 1750, "tous": 1237},
    "qualite_loyer_2025": {"studio": "c", "1 ch.": "d", "2 ch.": "c", "3 ch. +": "c", "tous": "c"},
    "loyer_2024": {"studio": 819, "1 ch.": 1093, "2 ch.": None, "3 ch. +": 1782, "tous": 1237},
    "univers_2025": {"studio": 5553, "1 ch.": 12405, "2 ch.": 22520, "3 ch. +": 4076, "tous": 44554},
    "inocc_2025": {"studio": None, "1 ch.": 0.5, "2 ch.": None, "3 ch. +": None, "tous": None},
    "inocc_2024": {"studio": 2.0, "1 ch.": 0.4, "2 ch.": None, "3 ch. +": 0.0, "tous": 0.7},
    "variation_2025": {"studio": None, "1 ch.": 8.8, "2 ch.": 6.3, "3 ch. +": 4.2, "tous": 7.5},
    "inocc_taille_2025": {"3 à 5": None, "6 à 19": None, "20 à 49": 4.1, "50 à 99": 0.0, "100 +": 2.8},
    "loyer_inoccupe_vs_occupe_2025": {"studio": (905, 1012), "1 ch.": (1295, 1054), "2 ch.": (None, 1285), "3 ch. +": (None, 1750), "tous": (None, 1213)},
    "neuf_2022_2025": {"univers": 627, "loyer_tous": 1920, "loyer_2ch": 2226, "univers_2ch": 306},
}
SCHL_ILE_6_19 = {"inocc_2025": 2.2, "inocc_2024": 1.7, "loyer_2ch_2025": 1177, "loyer_tous_2025": 1162, "rotation_2025": 8.5, "rotation_2024": 9.4}
SCHL_RMR = {"inocc_2025": 2.9, "inocc_2024": 2.1, "loyer_2ch_2025": 1346, "loyer_tous_2025": 1290, "rotation_2025": 9.4, "rotation_2024": 11.1,
            "loyer_inoccupe_tous": 1628, "loyer_occupe_tous": 1281, "loyer_inoccupe_2ch": 1817, "loyer_occupe_2ch": 1333, "univers_2025": 667266}

def etude_locative_cas3() -> dict:
    c = cas3.charger(); etat = c["etat"]; hyp = c["hyp"]; rev = cas3.revenus(etat, hyp)
    corr = {"3½": "1 ch.", "4½": "2 ch.", "5½": "3 ch. +"}
    g = etat.groupby("pieces").agg(n=("logement", "size"), contrat=("loyer_mensuel_contrat", "mean"), marche=("loyer_mensuel_marche", "mean"),
                                   superficie=("superficie_m2", "mean"))
    g["type_schl"] = [corr[p] for p in g.index]
    g["loyer_schl_zone9_2025"] = [SCHL_ZONE9["loyer_2025"][t] for t in g["type_schl"]]
    g["ratio_marche_schl"] = g["marche"] / g["loyer_schl_zone9_2025"]
    g["ecart_contrat_marche_pct"] = 100 * (g["marche"] / g["contrat"] - 1)
    u = SCHL_ZONE9["univers_2025"]
    parts = {k: u[k] / u["tous"] for k in ["studio", "1 ch.", "2 ch.", "3 ch. +"]}
    ecart_1ch = SCHL_ZONE9["loyer_inoccupe_vs_occupe_2025"]["1 ch."]
    ecart_rmr = (SCHL_RMR["loyer_inoccupe_tous"] / SCHL_RMR["loyer_occupe_tous"] - 1)
    rot = {"hypothese_cas3": hyp["rotation_annuelle"], "schl_rmr_2025": SCHL_RMR["rotation_2025"] / 100, "schl_ile_6_19_2025": SCHL_ILE_6_19["rotation_2025"] / 100}
    rattrapage = {k: r * rev["ecart_marche"] for k, r in rot.items()}          # part de l'écart récupérable la première année
    horizon = {k: 1.0 / r for k, r in rot.items()}                              # années pour renouveler tous les baux au rythme r
    return {"par_type": g, "revenus": rev, "ecart_global_pct": 100 * (rev["rbp_marche"] / rev["rbp_contrat"] - 1), "parts_univers": parts,
            "ecart_inoccupe_occupe_1ch_zone9_pct": 100 * (ecart_1ch[0] / ecart_1ch[1] - 1), "ecart_inoccupe_occupe_rmr_pct": 100 * ecart_rmr,
            "rotation": rot, "rattrapage_an1": rattrapage, "horizon_renouvellement": horizon,
            "neuf_prime_pct": 100 * (SCHL_ZONE9["neuf_2022_2025"]["loyer_tous"] / SCHL_ZONE9["loyer_2025"]["tous"] - 1),
            "neuf_part_univers_pct": 100 * SCHL_ZONE9["neuf_2022_2025"]["univers"] / u["tous"],
            "inocc_hyp_cas3_pct": 100 * hyp["taux_inoccupation_pertes"]}

# =============================================================================== tout
def tout() -> dict:
    d = cas1.charger(); v = d["ventes"]
    ind = indicateurs_apciq(); parts = parts_secteurs_gatineau(v)
    ref = cas1.hedonique(v)
    # segmentation a priori : secteurs ; autres partitions a priori (type, période)
    ch_secteur = chow(v, v["secteur"], "secteur")
    est = v["secteur"].replace({"Buckingham": "Est", "Masson-Angers": "Est"})
    ch_secteur4 = chow(v, est, "secteur (4 groupes, Est regroupé)")
    ch_type = chow(v, v["type"], "type de bâtiment")
    periode = np.where(v["date_vente"] >= pd.Timestamp("2025-01-01"), "2025-2026", "2023-2024")
    ch_periode = chow(v, pd.Series(periode, index=v.index), "période")
    # segmentation statistique
    Z = caracteristiques_segmentation(v, "attributs+coordonnees"); cl = classifier(Z)
    k = cl["k_opt"]; lab_km = cl["labels_kmeans"][k]; lab_w = cl["labels_ward"][k]
    sil_i = silhouette_samples(Z.values, lab_km)
    prof = profil_segments(v, lab_km + 1)
    croise = pd.crosstab(pd.Series(lab_km + 1, index=v.index, name="segment"), v["secteur"])
    ari_secteur = float(adjusted_rand_score(v["secteur"], lab_km)); ari_ward = float(adjusted_rand_score(lab_km, lab_w))
    jeux = {}
    for jeu in ["attributs", "coordonnees", "attributs+coordonnees"]:
        Zj = caracteristiques_segmentation(v, jeu); cj = classifier(Zj); kj = cj["k_opt"]
        jeux[jeu] = {"k_opt": kj, "silhouette": float(cj["tableau"].loc[cj["tableau"]["k"] == kj, "silhouette_kmeans"].iloc[0]),
                     "ari_secteur": float(adjusted_rand_score(v["secteur"], cj["labels_kmeans"][kj])),
                     "silhouette_k5": float(cj["tableau"].loc[cj["tableau"]["k"] == 5, "silhouette_kmeans"].iloc[0]),
                     "ari_secteur_k5": float(adjusted_rand_score(v["secteur"], cj["labels_kmeans"][5]))}
    ch_kmeans = chow(v, pd.Series(lab_km + 1, index=v.index), f"k-moyennes (k = {k})")
    ch_ward = chow(v, pd.Series(lab_w, index=v.index), f"Ward (k = {k})")
    cv = validation_croisee(v, {"k-moyennes": lab_km + 1, "Ward": lab_w})
    loc = etude_locative_cas3()
    return {"ventes": v, "sujet": d["sujet"], "indicateurs": ind, "parts_gatineau": parts, "reference": ref,
            "chow": {"secteur": ch_secteur, "secteur4": ch_secteur4, "type": ch_type, "periode": ch_periode, "kmeans": ch_kmeans, "ward": ch_ward},
            "Z": Z, "classification": cl, "k_opt": k, "labels_kmeans": lab_km, "labels_ward": lab_w, "silhouette_i": sil_i, "profil": prof,
            "croise": croise, "ari_secteur": ari_secteur, "ari_kmeans_ward": ari_ward, "jeux": jeux, "cv": cv, "locatif": loc}

def main():
    pd.set_option("display.width", 220); pd.set_option("display.max_columns", 30)
    r = tout()
    print("=== A. Indicateurs APCIQ T3 2026 (chiffres cités ; calculs de l'auteur)")
    print(r["indicateurs"].round(3).to_string(index=False))
    print("\nParts des ventes par secteur, Gatineau (APCIQ vs synthétique) :\n", r["parts_gatineau"].round(3).to_string())
    print("\n=== B. Modèle de référence (traite.cas1.hedonique) : R² = %.3f, σ = %.4f" % (r["reference"]["r2"], r["reference"]["sigma"]))
    for key, ch in r["chow"].items():
        print(f"\n--- Chow [{ch['nom']}] : G = {ch['G']}, k0 = {ch['k0']}, n = {ch['n']}")
        print(f"    SSR contraint (constantes par segment) = {ch['ssr_contraint']:.4f} ; SSR classique = {ch['ssr_classique']:.4f} ; SSR libre = {ch['ssr_libre']:.4f}")
        print(f"    F pentes = {ch['F']:.3f} (q = {ch['q']}, df2 = {ch['df2']}), p = {ch['p']:.4f} ; F classique = {ch['F_classique']:.3f} (q = {ch['q_classique']}), p = {ch['p_classique']:.4g}")
        print(f"    Wald HC1 = {ch['wald_hc1']:.2f} (q = {ch['q_wald']}), p = {ch['p_wald']:.4f} ; σ contraint {ch['sigma_contraint']:.4f} vs libre {ch['sigma_libre']:.4f}")
        if key in ("secteur", "kmeans"):
            print(tableau_coefficients(ch).round(4).to_string(index=False))
    print("\n--- Classification (attributs + coordonnées) :\n", r["classification"]["tableau"].round(3).to_string(index=False))
    print(f"    k optimal (silhouette k-moyennes) = {r['k_opt']} ; ARI(k-moyennes, secteurs) = {r['ari_secteur']:.3f} ; ARI(k-moyennes, Ward) = {r['ari_kmeans_ward']:.3f}")
    print("    silhouette moyenne par segment :", pd.Series(r["silhouette_i"]).groupby(r["labels_kmeans"] + 1).mean().round(3).to_dict())
    print("\n--- Profil des segments k-moyennes :\n", r["profil"].round(2).to_string())
    print("\n--- Segments × secteurs :\n", r["croise"].to_string())
    print("\n--- Jeux de variables :\n", pd.DataFrame(r["jeux"]).T.round(3).to_string())
    print("\n--- Validation croisée (10 plis) :\n", r["cv"].round(4).to_string(index=False))
    L = r["locatif"]
    print("\n=== C. Étude de marché locatif, cas 3 (SCHL zone 9, oct. 2025 ; état locatif synthétique)")
    print(L["par_type"].round(3).to_string())
    print("    RBP contrat %.0f $, RBP marché %.0f $, écart %.0f $ (%.1f %%)" % (L["revenus"]["rbp_contrat"], L["revenus"]["rbp_marche"], L["revenus"]["ecart_marche"], L["ecart_global_pct"]))
    print("    parts de l'univers zone 9 :", {k: round(x, 3) for k, x in L["parts_univers"].items()})
    print("    écart loyer inoccupé/occupé : 1 ch. zone 9 %.1f %% ; tous RMR %.1f %%" % (L["ecart_inoccupe_occupe_1ch_zone9_pct"], L["ecart_inoccupe_occupe_rmr_pct"]))
    print("    rotation :", {k: round(x, 3) for k, x in L["rotation"].items()}, "; rattrapage an 1 :", {k: round(x) for k, x in L["rattrapage_an1"].items()},
          "; horizon (ans) :", {k: round(x, 1) for k, x in L["horizon_renouvellement"].items()})
    print("    parc neuf 2022-2025 zone 9 : prime de loyer %.1f %%, part de l'univers %.2f %%" % (L["neuf_prime_pct"], L["neuf_part_univers_pct"]))
    return r

if __name__ == "__main__":
    main()