spb/vrai-prix

1#!/usr/bin/env python32# Ré-estimation du modèle hédonique DIRECTEMENT SUR LE NŒUD.3#4# Le pipeline batch historique (laptop, hedonic.py/LightGBM) produit est_2021..5# est_2026 + p10/p90 par unité. Depuis que le connecteur ingest-jdm.mjs ajoute6# les ventes récentes en continu dans vraiprix.db, on peut ré-entraîner ici même7# (sklearn HistGradientBoostingRegressor ≈ LightGBM ; 28 cœurs / 96 Go).8#9# Deux modes :10#   --eval   Bench ÉQUITABLE : le modèle laptop a vu les ventes ≤ 2026-07-2711#            (build 2026-08-08). On entraîne le challenger sur ce même12#            périmètre, puis on compare les deux sur les ventes > 2026-07-2713#            (jamais vues par aucun des deux). Rien n'est écrit.14#   --apply  Ré-entraîne sur TOUTES les ventes, prédit les 3,7 M d'unités15#            (médiane + P10/P90 par modèles quantiles), sauvegarde les anciennes16#            colonnes dans `units_est_prev`, puis met à jour units.est_2026/p10/p90.17#18# Cible : log(prix). Enrichissements vs les features brutes du rôle :19#   - te_ratio : médiane locale (code_mun × unite_voisinage, replis municipalité20#     puis global) de log(prix/valeur_role) — « ratio d'évaluation » local ;21#   - ppm2_loc : médiane locale du prix/m² ;22#   - t : date de vente en mois depuis 2021-01 (tendance de marché) ;23#   - pondération de récence (demi-vie 36 mois) pour coller au marché actuel.24#25# Usage : python3 scripts/hedonic-retrain.py --eval26#         python3 scripts/hedonic-retrain.py --apply2728import argparse29import sqlite330import sys31import time3233import numpy as np34import pandas as pd35from sklearn.ensemble import HistGradientBoostingRegressor3637DB = "data/vraiprix.db"38BASELINE_CUTOFF = "2026-07-27"  # dernière vente vue par le build laptop 2026-08-0839T0 = pd.Timestamp("2021-01-01")40AMOUNT_MIN, AMOUNT_MAX = 50_000, 10_000_0004142NUM_FEATS = [43    "lat", "lng", "annee_construction", "aire_etages_m2", "superficie_terrain_m2",44    "front_terrain_m", "nb_etages", "nb_logements", "nb_locaux_non_resid",45    "nb_chambres_locatives", "n_adresses", "valeur_terrain", "valeur_batiment",46    "valeur_role", "valeur_anterieure",47]48CAT_FEATS = ["type_prop", "lien_physique", "genre_construction"]49# naive_role / naive_ppm2 : « estimés naïfs » locaux donnés explicitement au GBM50# (les arbres capturent mal les relations multiplicatives rôle × ratio local).51ENG_FEATS = ["t", "log_role", "te_ratio", "ppm2_loc",52             "naive_role", "naive_ppm2", "age", "log_aire", "log_terr"]535455def log(msg):56    print(f"[hedonic] {msg}", flush=True)575859def load_sales(con):60    q = f"""61    SELECT t.date, t.amount,62           u.id_provinc, u.code_mun, u.unite_voisinage, u.municipalite,63           u.lat, u.lng, u.annee_construction, u.aire_etages_m2,64           u.superficie_terrain_m2, u.front_terrain_m, u.nb_etages,65           u.nb_logements, u.nb_locaux_non_resid, u.nb_chambres_locatives,66           u.n_adresses, u.valeur_terrain, u.valeur_batiment, u.valeur_role,67           u.valeur_anterieure, u.type_prop, u.lien_physique, u.genre_construction68    FROM transactions t69    JOIN units u ON u.id_provinc = t.id_provinc70    WHERE t.amount BETWEEN {AMOUNT_MIN} AND {AMOUNT_MAX}71    """72    df = pd.read_sql_query(q, con)73    df["t"] = (pd.to_datetime(df["date"]) - T0).dt.days / 30.4474    return df757677# --- encodages locaux (calculés sur l'ENTRAÎNEMENT seulement, replis lissés) ---78class LocalEncoder:79    K = 8  # lissage : poids du repli8081    def fit(self, df):82        d = df[(df["valeur_role"] > 0)].copy()83        d["ratio"] = np.log(d["amount"] / d["valeur_role"])84        self.g_ratio = d["ratio"].median()85        m = d.groupby("municipalite")["ratio"].agg(["median", "size"])86        self.mun_ratio = ((m["median"] * m["size"] + self.g_ratio * self.K)87                          / (m["size"] + self.K)).to_dict()88        d["vk"] = d["code_mun"].astype(str) + "|" + d["unite_voisinage"].astype(str)89        v = d.groupby(["vk", "municipalite"])["ratio"].agg(["median", "size"]).reset_index()90        v["fb"] = v["municipalite"].map(self.mun_ratio).fillna(self.g_ratio)91        v["val"] = (v["median"] * v["size"] + v["fb"] * self.K) / (v["size"] + self.K)92        self.vois_ratio = dict(zip(v["vk"], v["val"]))9394        p = df[(df["aire_etages_m2"] > 20)].copy()95        p["ppm2"] = p["amount"] / p["aire_etages_m2"]96        self.g_ppm2 = p["ppm2"].median()97        m2 = p.groupby("municipalite")["ppm2"].agg(["median", "size"])98        self.mun_ppm2 = ((m2["median"] * m2["size"] + self.g_ppm2 * self.K)99                         / (m2["size"] + self.K)).to_dict()100        p["vk"] = p["code_mun"].astype(str) + "|" + p["unite_voisinage"].astype(str)101        v2 = p.groupby(["vk", "municipalite"])["ppm2"].agg(["median", "size"]).reset_index()102        v2["fb"] = v2["municipalite"].map(self.mun_ppm2).fillna(self.g_ppm2)103        v2["val"] = (v2["median"] * v2["size"] + v2["fb"] * self.K) / (v2["size"] + self.K)104        self.vois_ppm2 = dict(zip(v2["vk"], v2["val"]))105        return self106107    def transform(self, df):108        vk = df["code_mun"].astype(str) + "|" + df["unite_voisinage"].astype(str)109        mun_r = df["municipalite"].map(self.mun_ratio).fillna(self.g_ratio)110        df["te_ratio"] = vk.map(self.vois_ratio).fillna(mun_r)111        mun_p = df["municipalite"].map(self.mun_ppm2).fillna(self.g_ppm2)112        df["ppm2_loc"] = vk.map(self.vois_ppm2).fillna(mun_p)113        return df114115116def featurize(df, enc, cat_maps=None):117    df = enc.transform(df.copy())118    df["log_role"] = np.log1p(df["valeur_role"].clip(lower=0))119    df["naive_role"] = np.log1p((df["valeur_role"].clip(lower=0)120                                 * np.exp(df["te_ratio"])).clip(lower=0))121    df["naive_ppm2"] = np.log1p((df["aire_etages_m2"].clip(lower=0)122                                 * df["ppm2_loc"]).clip(lower=0))123    df["age"] = (2021 + df["t"] / 12.0) - df["annee_construction"]124    df["log_aire"] = np.log1p(df["aire_etages_m2"].clip(lower=0))125    df["log_terr"] = np.log1p(df["superficie_terrain_m2"].clip(lower=0))126    if cat_maps is None:127        cat_maps = {c: {v: i for i, v in enumerate(df[c].astype(str).fillna("NA").unique())}128                    for c in CAT_FEATS}129    for c in CAT_FEATS:130        df[c] = df[c].astype(str).fillna("NA").map(cat_maps[c]).fillna(-1).astype(int)131    cols = NUM_FEATS + ENG_FEATS + CAT_FEATS132    X = df[cols].astype(float).to_numpy()133    return X, cols, cat_maps134135136def make_model(loss="absolute_error", quantile=None, cols=None, light=False):137    # Config volontairement SIMPLE (≈ 8 min d'entraînement) : la variante lourde138    # (2500 it / 255 feuilles, ~40 min) ne gagnait que ~1-2 pts de MdAPE sur les139    # segments remplacés — pas rentable.140    cat_mask = [c in CAT_FEATS for c in cols]141    kw = dict(142        max_iter=600 if light else 900, learning_rate=0.06, max_leaf_nodes=127,143        min_samples_leaf=40, l2_regularization=0.1, max_bins=255,144        early_stopping=True, validation_fraction=0.05, n_iter_no_change=40,145        random_state=42, categorical_features=cat_mask,146    )147    if quantile is not None:148        return HistGradientBoostingRegressor(loss="quantile", quantile=quantile, **kw)149    return HistGradientBoostingRegressor(loss=loss, **kw)150151152def recency_weights(t, t_max, half_life=24.0):153    return np.power(0.5, (t_max - t) / half_life)154155156def metrics(y_true, y_pred, label):157    ape = np.abs(y_pred - y_true) / y_true158    return {159        "label": label, "n": len(y_true),160        "MdAPE": float(np.median(ape) * 100),161        "MAPE": float(np.mean(ape) * 100),162        "±10%": float(np.mean(ape <= 0.10) * 100),163        "±20%": float(np.mean(ape <= 0.20) * 100),164    }165166167def print_metrics(rows):168    hdr = f"{'modèle':<26}{'n':>7}{'MdAPE':>8}{'MAPE':>8}{'±10%':>7}{'±20%':>7}"169    print(hdr); print("-" * len(hdr))170    for r in rows:171        print(f"{r['label']:<26}{r['n']:>7}{r['MdAPE']:>7.2f}%{r['MAPE']:>7.2f}%"172              f"{r['±10%']:>6.1f}%{r['±20%']:>6.1f}%")173174175def run_eval(con):176    df = load_sales(con)177    log(f"ventes jointes exploitables : {len(df)}")178    train = df[df["date"] <= BASELINE_CUTOFF].copy()179    test = df[df["date"] > BASELINE_CUTOFF].copy()180    log(f"train (≤ {BASELINE_CUTOFF}) : {len(train)} | test (>) : {len(test)}")181182    enc = LocalEncoder().fit(train)183    Xtr, cols, cat_maps = featurize(train, enc)184    ytr = np.log(train["amount"].to_numpy())185    w = recency_weights(train["t"].to_numpy(), train["t"].max())186187    t0 = time.time()188    model = make_model(cols=cols)189    model.fit(Xtr, ytr, sample_weight=w)190    log(f"challenger entraîné ({model.n_iter_} itérations, {time.time()-t0:.0f}s)")191192    Xte, _, _ = featurize(test, enc, cat_maps)193    pred = np.exp(model.predict(Xte))194    y = test["amount"].to_numpy()195196    # baseline : est_2026 de l'unité jointe197    base = pd.read_sql_query(198        "SELECT id_provinc, est_2026, p10, p90 FROM units", con)199    test = test.merge(base, on="id_provinc", how="left")200    mask = test["est_2026"].notna().to_numpy()201202    # mélange géométrique laptop × challenger (moyenne sur l'échelle log)203    base_est = test["est_2026"].to_numpy()204    blend = np.where(mask, np.exp(0.5 * np.log(np.where(mask, base_est, 1))205                                  + 0.5 * np.log(pred)), pred)206207    rows = [208        metrics(y[mask], base_est[mask], "laptop est_2026"),209        metrics(y[mask], pred[mask], "challenger (nœud)"),210        metrics(y[mask], blend[mask], "blend 50/50 (géo)"),211    ]212    print(); print(f"=== Ventes jamais vues (> {BASELINE_CUTOFF}) — {mask.sum()} obs ===")213    print_metrics(rows)214215    # par type216    print("\n--- par type de propriété (MdAPE %) ---")217    tt = test[mask].copy(); tt["pred"] = pred[mask]; tt["y"] = y[mask]218    tt["blend"] = blend[mask]219    for tp, g in tt.groupby(test[mask]["type_prop"]):220        b = np.median(np.abs(g["est_2026"] - g["y"]) / g["y"]) * 100221        c = np.median(np.abs(g["pred"] - g["y"]) / g["y"]) * 100222        bl = np.median(np.abs(g["blend"] - g["y"]) / g["y"]) * 100223        print(f"  {tp:<16} n={len(g):>5}  laptop {b:6.2f}%  challenger {c:6.2f}%  blend {bl:6.2f}%")224225    # couverture des quantiles baseline226    cov10 = float(np.mean(y[mask] < test["p10"].to_numpy()[mask]) * 100)227    cov90 = float(np.mean(y[mask] > test["p90"].to_numpy()[mask]) * 100)228    print(f"\ncouverture P10/P90 laptop sur test : {cov10:.1f}% sous P10 (cible 10) | "229          f"{cov90:.1f}% au-dessus de P90 (cible 10)")230    return rows231232233# Types dont l'estimation laptop est remplacée par le challenger (gain validé234# sur ventes jamais vues : terrain MdAPE 62→36 %, autre 79→52 %).235REPLACE_TYPES = ("terrain", "autre")236237238def quantile_calibration(con, df):239    """Facteurs d'élargissement des P10/P90 laptop, calibrés sur les ventes240    jamais vues (> BASELINE_CUTOFF), hors types remplacés. Couverture observée241    18,5 %/16,8 % hors bornes → cible 10 %/10 % de chaque côté."""242    test = df[df["date"] > BASELINE_CUTOFF].copy()243    base = pd.read_sql_query(244        "SELECT id_provinc, est_2026, p10, p90, type_prop tp FROM units", con)245    test = test.merge(base, on="id_provinc", how="inner")246    test = test[test["est_2026"].notna() & (test["est_2026"] > 0)247                & (test["p10"] > 0) & (test["p90"] > 0)248                & ~test["tp"].isin(REPLACE_TYPES)]249    z = np.log(test["amount"] / test["est_2026"])250    r_lo = np.log(test["p10"] / test["est_2026"])   # < 0251    r_hi = np.log(test["p90"] / test["est_2026"])   # > 0252    ok = (r_lo < -1e-6) & (r_hi > 1e-6)253    s_lo = (z[ok] / r_lo[ok])   # >1 ⇒ sous P10254    s_hi = (z[ok] / r_hi[ok])   # >1 ⇒ au-dessus de P90255    a_lo = float(np.quantile(s_lo, 0.90))  # P(s_lo > a_lo) = 10 %256    a_hi = float(np.quantile(s_hi, 0.90))257    a_lo, a_hi = max(1.0, a_lo), max(1.0, a_hi)258    log(f"calibration quantiles laptop (n={ok.sum()}) : alpha_lo={a_lo:.3f}, "259        f"alpha_hi={a_hi:.3f}")260    return a_lo, a_hi261262263def run_apply(con):264    df = load_sales(con)265    a_lo, a_hi = quantile_calibration(con, df)266267    log(f"ré-entraînement final sur {len(df)} ventes (tout l'historique)")268    enc = LocalEncoder().fit(df)269    X, cols, cat_maps = featurize(df, enc)270    yl = np.log(df["amount"].to_numpy())271    w = recency_weights(df["t"].to_numpy(), df["t"].max())272273    t0 = time.time()274    med = make_model(cols=cols); med.fit(X, yl, sample_weight=w)275    log(f"modèle médian : {med.n_iter_} it, {time.time()-t0:.0f}s")276    t0 = time.time()277    q10 = make_model(cols=cols, quantile=0.10, light=True); q10.fit(X, yl, sample_weight=w)278    q90 = make_model(cols=cols, quantile=0.90, light=True); q90.fit(X, yl, sample_weight=w)279    log(f"modèles quantiles P10/P90 : {time.time()-t0:.0f}s")280281    # --- unités des types remplacés, prédites à la date « maintenant » ---282    t_now = float(df["t"].max())283    ph = ",".join("?" * len(REPLACE_TYPES))284    units = pd.read_sql_query(f"""285        SELECT id_provinc, code_mun, unite_voisinage, municipalite,286               {', '.join(NUM_FEATS)}, {', '.join(CAT_FEATS)}287        FROM units WHERE type_prop IN ({ph})""", con, params=REPLACE_TYPES)288    log(f"unités à ré-estimer ({' + '.join(REPLACE_TYPES)}) : {len(units)}")289    units["t"] = t_now290    Xu, _, _ = featurize(units, enc, cat_maps)291292    log("prédiction challenger (médiane + P10/P90)…")293    est = np.exp(med.predict(Xu))294    lo = np.exp(q10.predict(Xu))295    hi = np.exp(q90.predict(Xu))296    lo2 = np.minimum.reduce([lo, est, hi]); hi2 = np.maximum.reduce([lo, est, hi])297    out = pd.DataFrame({298        "id_provinc": units["id_provinc"],299        "est": np.clip(est, 1000, None).round(-2),300        "p10": np.clip(lo2, 1000, None).round(-2),301        "p90": np.clip(hi2, 1000, None).round(-2),302    })303304    log("écriture en base (sauvegarde units_est_prev puis UPDATE)…")305    cur = con.cursor()306    cur.execute("PRAGMA busy_timeout=30000")307    cur.execute("DROP TABLE IF EXISTS units_est_prev")308    cur.execute("""CREATE TABLE units_est_prev AS309                   SELECT id_provinc, est_2026, p10, p90 FROM units""")310    con.commit()311312    # 1) types remplacés → challenger313    cur.execute("DROP TABLE IF EXISTS units_est_new")314    cur.execute("""CREATE TABLE units_est_new315                   (id_provinc TEXT PRIMARY KEY, est REAL, p10 REAL, p90 REAL)""")316    cur.executemany("INSERT OR REPLACE INTO units_est_new VALUES (?,?,?,?)",317                    out.itertuples(index=False, name=None))318    con.commit()319    cur.execute("""UPDATE units SET320                     est_2026 = (SELECT est FROM units_est_new n WHERE n.id_provinc = units.id_provinc),321                     p10      = (SELECT p10 FROM units_est_new n WHERE n.id_provinc = units.id_provinc),322                     p90      = (SELECT p90 FROM units_est_new n WHERE n.id_provinc = units.id_provinc)323                   WHERE id_provinc IN (SELECT id_provinc FROM units_est_new)""")324    con.commit()325    cur.execute("DROP TABLE units_est_new")326    con.commit()327328    # 2) autres types → P10/P90 laptop élargis (p' = est × (p/est)^alpha)329    ph2 = ",".join("?" * len(REPLACE_TYPES))330    cur.execute(f"""UPDATE units SET331                      p10 = ROUND(est_2026 * POW(p10 / est_2026, ?), -2),332                      p90 = ROUND(est_2026 * POW(p90 / est_2026, ?), -2)333                    WHERE type_prop NOT IN ({ph2})334                      AND est_2026 > 0 AND p10 > 0 AND p90 > 0""",335                (a_lo, a_hi, *REPLACE_TYPES))336    con.commit()337    cur.execute("PRAGMA wal_checkpoint(TRUNCATE)")338    n = cur.execute(339        f"SELECT COUNT(*) FROM units WHERE type_prop IN ({ph2}) AND est_2026 IS NOT NULL",340        REPLACE_TYPES).fetchone()[0]341    log(f"terminé : {n} unités ré-estimées (challenger), P10/P90 recalibrés ailleurs "342        f"(alpha {a_lo:.3f}/{a_hi:.3f}). Ancienne version dans units_est_prev.")343344345if __name__ == "__main__":346    ap = argparse.ArgumentParser()347    ap.add_argument("--eval", action="store_true")348    ap.add_argument("--apply", action="store_true")349    ap.add_argument("--db", default=DB)350    a = ap.parse_args()351    con = sqlite3.connect(a.db)352    try:353        if a.eval:354            run_eval(con)355        elif a.apply:356            run_apply(con)357        else:358            ap.print_help(); sys.exit(1)359    finally:360        con.close()361