#!/usr/bin/env python3
# Auteur : Simon-Pierre Boucher — contact@spboucher.ai
"""ValoPlex — modèle hédonique avancé spécialisé plex (2+ logements, CUBF 1000).

Avancées vs modèle généraliste Vrai-Prix :
  · CIBLE EN RATIO : le modèle prédit log(prix / valeur au rôle), pas le prix.
    Invariant d'échelle, il évalue aussi justement le duplex de 400 k$ que la
    tour au rôle de 700 M$ (un modèle en prix brut compresse les gros
    immeubles vers la masse des duplex — bogue observé et corrigé) ;
  · domaine restreint aux plex, features « par porte » (aire/porte,
    valeur au rôle/porte, bâtiment $/m²) et signaux de mixité ;
  · fourchette P10-P90 par régression quantile PUIS calibration conforme
    (split conformal) : la couverture 80 % est garantie empiriquement sur un
    ensemble de calibration jamais vu ;
  · validation IAAO segmentée (par nombre de portes, par année, temporel 2026).

Usage : python3 hedonic_plex.py train | predict
"""
import json
import os
import sys

import joblib
import lightgbm as lgb
import numpy as np
import pandas as pd

QC = "/Users/simon-pierreboucher/Desktop/qc_house_eval"
BASE = "/Users/simon-pierreboucher/Desktop/valoplex"
MODELS = f"{BASE}/models"
ESTIM = f"{BASE}/data/estimations"
os.makedirs(MODELS, exist_ok=True)
os.makedirs(ESTIM, exist_ok=True)

T0 = pd.Timestamp("2021-01-01")

NUM = ["portes", "aire_etages_m2", "aire_par_porte", "superficie_terrain_m2",
       "front_terrain_m", "nb_etages", "nb_chambres_locatives",
       "nb_locaux_non_resid", "age", "lat", "lng", "n_adresses",
       "t_mois", "mois", "anciennete_role_mois",
       "log_val_terrain", "log_val_batiment", "log_val_immeuble",
       "val_par_porte", "bat_par_m2", "part_terrain", "croissance_role"]
CAT = ["lien_physique", "genre_construction", "code_mun_i", "region"]
FEATS = NUM + CAT
LGB = dict(objective="regression", metric="l1", num_leaves=192,
           learning_rate=0.045, n_estimators=4000, min_child_samples=30,
           colsample_bytree=0.85, subsample=0.85, subsample_freq=1,
           reg_lambda=1.5, n_jobs=16, verbose=-1)


def c_tx(name):  # colonnes côté transactions enrichies
    m = {"rl0308a": "role_aire_etages_m2", "rl0302a": "role_superficie_terrain_m2",
         "rl0301a": "role_front_terrain_m", "rl0306a": "role_nb_etages",
         "rl0311a": "role_nb_logements", "rl0312a": "role_nb_locaux_non_resid",
         "rl0313a": "role_nb_chambres_locatives", "lat": "role_lat",
         "lng": "role_lng", "n_adresses": "role_n_adresses",
         "rl0309a": "role_lien_physique", "rl0310a": "role_genre_construction",
         "code_mun": "role_code_mun", "rl0404a": "role_valeur_immeuble",
         "rl0402a": "role_valeur_terrain", "rl0403a": "role_valeur_batiment",
         "rl0405a": "role_valeur_role_anterieur", "rl0307a": "role_annee_construction",
         "dat_cond_mrche": "role_date_cond_marche"}
    return m[name]


def eng(df, col, when):
    """Features plex. `col` : fonction nom générique -> nom de colonne."""
    o = pd.DataFrame(index=df.index)
    num = lambda n: pd.to_numeric(df[col(n)], errors="coerce")
    portes = num("rl0311a")
    aire = num("rl0308a")
    val = num("rl0404a")
    terr = num("rl0402a")
    bat = num("rl0403a")
    prev = num("rl0405a")
    yb = num("rl0307a")

    o["portes"] = portes
    o["aire_etages_m2"] = aire
    o["aire_par_porte"] = (aire / portes.replace(0, np.nan)).clip(10, 500)
    o["superficie_terrain_m2"] = num("rl0302a")
    o["front_terrain_m"] = num("rl0301a")
    o["nb_etages"] = num("rl0306a")
    o["nb_chambres_locatives"] = num("rl0313a")
    o["nb_locaux_non_resid"] = num("rl0312a")
    o["age"] = (when.dt.year - yb).where((yb > 1600) & (yb <= when.dt.year + 2))
    o["lat"] = num("lat")
    o["lng"] = num("lng")
    o["n_adresses"] = num("n_adresses")
    o["t_mois"] = (when - T0).dt.days / 30.44
    o["mois"] = when.dt.month.astype(float)
    # ancienneté du rôle : mois entre la date des conditions du marché du rôle
    # (dat_cond_mrche) et la date d'évaluation — LA variable d'un modèle en ratio :
    # plus le rôle est vieux, plus le marché s'en est éloigné
    dcm = pd.to_datetime(df[col("dat_cond_mrche")].astype(str).str[:10], errors="coerce")
    o["anciennete_role_mois"] = ((when - dcm).dt.days / 30.44).clip(-6, 72)
    o["log_val_terrain"] = np.log1p(terr.clip(lower=0))
    o["log_val_batiment"] = np.log1p(bat.clip(lower=0))
    o["log_val_immeuble"] = np.log1p(val.clip(lower=0))
    o["val_par_porte"] = (val / portes.replace(0, np.nan)).clip(5_000, 2_000_000)
    o["bat_par_m2"] = (bat / aire.replace(0, np.nan)).clip(50, 20_000)
    o["part_terrain"] = (terr / val.replace(0, np.nan)).clip(0, 1)
    o["croissance_role"] = (val / prev.replace(0, np.nan)).clip(0.2, 5)
    o["lien_physique"] = df[col("rl0309a")].astype(str).astype("category")
    o["genre_construction"] = df[col("rl0310a")].astype(str).astype("category")
    cm = pd.to_numeric(df[col("code_mun")].astype(str).str.strip(), errors="coerce")
    o["code_mun_i"] = cm.fillna(-1).astype(int).astype("category")
    o["region"] = (cm // 1000).fillna(-1).astype(int).astype("category")
    return o


def align(X, levels=None):
    if levels is None:
        return {k: X[k].cat.categories for k in CAT}
    for k in CAT:
        X[k] = pd.Categorical(X[k], categories=levels[k])
    return levels


def build_training():
    df = pd.read_parquet(f"{QC}/province_transactions_enrichi.parquet")
    lg_ = pd.to_numeric(df["role_nb_logements"], errors="coerce")
    cu = pd.to_numeric(df["role_cubf"], errors="coerce")
    val = pd.to_numeric(df["role_valeur_immeuble"], errors="coerce")
    ratio = df["amount"] / val.replace(0, np.nan)
    keep = (
        (cu == 1000) & (lg_ >= 2) & df["role_id_provinc"].notna() &
        (df["match_valeur_exacte"] | (df["match_dist_m"] <= 50)) &
        df["amount"].between(60_000, 30_000_000) & ratio.between(0.25, 4.0)
    )
    df = df[keep].copy()
    df["when"] = pd.to_datetime(df["date"])
    print(f"transactions plex retenues : {len(df)}", flush=True)
    val_ok = pd.to_numeric(df["role_valeur_immeuble"], errors="coerce")
    df = df[val_ok > 10_000].copy()  # cible en ratio : rôle > 0 requis
    X = eng(df, c_tx, df["when"])
    base = pd.to_numeric(df["role_valeur_immeuble"], errors="coerce").values
    y = np.log(df["amount"].values / base)  # log(prix / valeur au rôle)
    meta = df[["id", "date", "amount", "tx_year"]].reset_index(drop=True)
    meta["portes"] = pd.to_numeric(df["role_nb_logements"], errors="coerce").values
    meta["base"] = base
    meta["valeur_role"] = base
    return X.reset_index(drop=True), y, meta


def iaao(pred, price):
    r = pred / price
    med = np.median(r)
    cod = 100 * np.mean(np.abs(r - med)) / med
    prd = np.mean(r) / (np.sum(pred) / np.sum(price))
    proxy = 0.5 * (pred / med + price)
    z = np.log(proxy) / np.log(2)
    beta = np.linalg.lstsq(np.c_[np.ones(len(z)), z], r / med - 1, rcond=None)[0]
    return dict(ratio_median=round(float(med), 4), COD=round(float(cod), 2),
                PRD=round(float(prd), 4), PRB=round(float(beta[1]), 4))


def metrics(pred, price):
    if len(price) == 0:
        return {}
    ape = np.abs(pred - price) / price
    return dict(n=int(len(price)), MdAPE=round(float(np.median(ape)) * 100, 2),
                MAPE=round(float(np.mean(ape)) * 100, 2),
                dans_10pct=round(float((ape <= .10).mean()) * 100, 1),
                dans_20pct=round(float((ape <= .20).mean()) * 100, 1),
                R2_log=round(float(1 - np.var(np.log(pred) - np.log(price)) /
                                   np.var(np.log(price))), 4),
                **iaao(pred, price))


BANDS = [(2, 2), (3, 3), (4, 5), (6, 12), (13, 10**9)]


def band_of(doors):
    import numpy as _np
    b = _np.full(len(doors), len(BANDS) - 1)
    for i, (lo, hi) in enumerate(BANDS):
        b[(doors >= lo) & (doors <= hi)] = i
    return b


def shrink_weight(doors):
    """Poids du modèle selon le support des ventes : 1 dans le domaine,
    décroissant hors domaine (13+ portes = extrapolation)."""
    import numpy as _np
    w = _np.ones(len(doors))
    w[(doors >= 13) & (doors <= 24)] = 0.55
    w[(doors >= 25) & (doors <= 49)] = 0.35
    w[doors >= 50] = 0.20
    return w


def train():
    X, y, meta = build_training()
    rng = np.random.default_rng(20260809)
    n = len(X)
    u = rng.random(n)
    test = u < 0.15                    # évaluation finale
    calib = (u >= 0.15) & (u < 0.30)   # calibration conforme
    fit = u >= 0.30                    # entraînement
    report = {"n_total": int(n), "n_fit": int(fit.sum()),
              "n_calib": int(calib.sum()), "n_test": int(test.sum())}

    print("=== modèle central (monotone) ===", flush=True)
    m = lgb.LGBMRegressor(**LGB)
    m.fit(X[FEATS][fit], y[fit], eval_set=[(X[FEATS][test], y[test])],
          callbacks=[lgb.early_stopping(200, verbose=False)])
    best = int(m.best_iteration_ or LGB["n_estimators"])
    report["best_iter"] = best
    sm = float(np.mean(np.exp(y[fit] - m.predict(X[FEATS][fit]))))
    report["smearing"] = round(sm, 4)
    base_test = meta.loc[test, "base"].values
    pred_test = base_test * np.exp(m.predict(X[FEATS][test])) * sm
    price_test = base_test * np.exp(y[test])
    report["holdout"] = metrics(pred_test, price_test)

    # segmentation par portes
    pt = meta.loc[test, "portes"].values
    vr = meta.loc[test, "valeur_role"].values
    for lab, mask in [("2 portes", pt == 2), ("3 portes", pt == 3),
                      ("4-5 portes", (pt >= 4) & (pt <= 5)),
                      ("6-12 portes", (pt >= 6) & (pt <= 12)),
                      ("13+ portes", pt >= 13),
                      ("rôle > 2M$", vr > 2_000_000),
                      ("rôle > 5M$", vr > 5_000_000)]:
        report[f"seg_{lab}"] = metrics(pred_test[mask], price_test[mask])
    # temporel
    t26 = (meta["tx_year"] == 2026).values
    m26 = lgb.LGBMRegressor(**{**LGB, "n_estimators": best})
    m26.fit(X[FEATS][~t26], y[~t26])
    sm26 = float(np.mean(np.exp(y[~t26] - m26.predict(X[FEATS][~t26]))))
    b26 = meta.loc[t26, "base"].values
    report["test_2026"] = metrics(b26 * np.exp(m26.predict(X[FEATS][t26])) * sm26,
                                  b26 * np.exp(y[t26]))

    print("=== quantiles + calibration conforme ===", flush=True)
    qmods = {}
    for a, nm in [(0.1, "p10"), (0.9, "p90")]:
        mq = lgb.LGBMRegressor(**{**LGB, "objective": "quantile", "alpha": a,
                                  "metric": "quantile", "n_estimators": best})
        mq.fit(X[FEATS][fit], y[fit])
        qmods[nm] = mq
    lo_c = qmods["p10"].predict(X[FEATS][calib])
    hi_c = qmods["p90"].predict(X[FEATS][calib])
    # score de non-conformité (intervalle symétrisé) et quantile 80 %
    s = np.maximum(lo_c - y[calib], y[calib] - hi_c)
    k = int(np.ceil(0.80 * (calib.sum() + 1))) - 1
    qhat = float(np.sort(s)[min(k, len(s) - 1)])
    report["conformal_qhat_log"] = round(qhat, 4)
    # couverture avant/après sur le test
    lo_t = qmods["p10"].predict(X[FEATS][test])
    hi_t = qmods["p90"].predict(X[FEATS][test])
    report["couverture_brute_pct"] = round(float(((y[test] >= lo_t) & (y[test] <= hi_t)).mean()) * 100, 1)
    report["couverture_conforme_pct"] = round(float(((y[test] >= lo_t - qhat) & (y[test] <= hi_t + qhat)).mean()) * 100, 1)

    # --- calibration Mondrian (par bande de portes) + ancre de rétrécissement ---
    doors_cal = meta.loc[calib, "portes"].values
    bc = band_of(doors_cal)
    band_stats = []
    for i, (blo, bhi) in enumerate(BANDS):
        m_b = bc == i
        yb = y[calib][m_b]
        if m_b.sum() >= 30:
            sb = np.maximum(lo_c[m_b] - yb, yb - hi_c[m_b])
            kb = int(np.ceil(0.80 * (m_b.sum() + 1))) - 1
            qhat_b = float(np.sort(sb)[min(kb, len(sb) - 1)])
        else:
            qhat_b = qhat  # repli global si bande trop mince
        band_stats.append({
            "lo": blo, "hi": bhi, "n_calib": int(m_b.sum()),
            "qhat": qhat_b,
            "ratio_median": float(np.exp(np.median(yb))) if m_b.sum() >= 10 else 1.0,
            "ratio_q10": float(np.exp(np.quantile(yb, 0.10))) if m_b.sum() >= 10 else None,
            "ratio_q90": float(np.exp(np.quantile(yb, 0.90))) if m_b.sum() >= 10 else None,
        })
    report["band_stats"] = band_stats
    # couverture Mondrian sur le test
    bt = band_of(meta.loc[test, "portes"].values)
    q_t = np.array([band_stats[i]["qhat"] for i in bt])
    report["couverture_mondrian_pct"] = round(float(((y[test] >= lo_t - q_t) & (y[test] <= hi_t + q_t)).mean()) * 100, 1)

    # modèle final : toutes les données
    print("=== réentraînement final ===", flush=True)
    mf = lgb.LGBMRegressor(**{**LGB, "n_estimators": int(best * 1.15)})
    mf.fit(X[FEATS], y)
    sm_f = float(np.mean(np.exp(y - mf.predict(X[FEATS]))))
    qf = {}
    for a, nm in [(0.1, "p10"), (0.9, "p90")]:
        mq = lgb.LGBMRegressor(**{**LGB, "objective": "quantile", "alpha": a,
                                  "metric": "quantile", "n_estimators": best})
        mq.fit(X[FEATS], y)
        qf[nm] = mq

    imp = pd.Series(mf.feature_importances_, index=FEATS).sort_values(ascending=False)
    report["importance"] = {k: int(v) for k, v in imp.head(12).items()}
    joblib.dump({"model": mf, "smearing": sm_f, "p10": qf["p10"], "p90": qf["p90"],
                 "qhat": qhat, "band_stats": band_stats, "levels": align(X),
                 "feats": FEATS},
                f"{MODELS}/valoplex_lgb.joblib")
    with open(f"{MODELS}/rapport_validation.json", "w") as f:
        json.dump(report, f, indent=2, ensure_ascii=False, default=str)
    print(json.dumps({k: v for k, v in report.items() if k != "importance"},
                     indent=2, ensure_ascii=False))
    print("TRAIN_DONE")


def predict():
    b = joblib.load(f"{MODELS}/valoplex_lgb.joblib")
    mf, sm, qhat = b["model"], b["smearing"], b["qhat"]
    idx = pd.read_csv(f"{QC}/data/indexRole2026.csv")
    muns = dict(zip(idx.iloc[:, 0].astype(str).str.lstrip("0"), idx.iloc[:, 1]))
    ident = lambda n: n
    for Y in range(2021, 2027):
        df = pd.read_parquet(f"{QC}/data/parquet/role_{Y}.parquet")
        cu = pd.to_numeric(df["rl0105a"], errors="coerce")
        lg_ = pd.to_numeric(df["rl0311a"], errors="coerce")
        df = df[(cu == 1000) & (lg_ >= 2)].reset_index(drop=True)
        when = pd.Series(pd.Timestamp(f"{Y}-06-15"), index=df.index)
        X = eng(df, ident, when)
        align(X, b["levels"])
        base = pd.to_numeric(df["rl0404a"], errors="coerce").values
        ok = base > 10_000
        doors_arr = pd.to_numeric(df["rl0311a"], errors="coerce").fillna(2).values
        bands_arr = band_of(doors_arr)
        bs = b["band_stats"]
        pl_raw = mf.predict(X[b["feats"]])
        # rétrécissement vers l'ancre empirique de la bande hors domaine
        w = shrink_weight(doors_arr)
        anchor = np.array([np.log(max(bs[i]["ratio_median"], 0.05)) for i in bands_arr])
        pl = w * pl_raw + (1 - w) * anchor
        # intervalles : quantiles + qhat de bande dans le domaine ;
        # dispersion empirique de bande recentrée hors domaine (13+ portes)
        q_b = np.array([bs[i]["qhat"] for i in bands_arr])
        lo = b["p10"].predict(X[b["feats"]]) - q_b
        hi = b["p90"].predict(X[b["feats"]]) + q_b
        big = doors_arr >= 13
        if big.any():
            med_b = np.array([np.log(max(bs[i]["ratio_median"], 0.05)) for i in bands_arr])
            q10_b = np.array([np.log(bs[i]["ratio_q10"]) if bs[i]["ratio_q10"] else med_b[k] - 0.4
                              for k, i in enumerate(bands_arr)])
            q90_b = np.array([np.log(bs[i]["ratio_q90"]) if bs[i]["ratio_q90"] else med_b[k] + 0.4
                              for k, i in enumerate(bands_arr)])
            lo[big] = pl[big] + (q10_b[big] - med_b[big])
            hi[big] = pl[big] + (q90_b[big] - med_b[big])
        out = pd.DataFrame({
            "id_provinc": df["id_provinc"], "annee": Y,
            "code_mun": df["code_mun"], "arrond": df["arrond"],
            "adresse": df["adresse"], "lat": df["lat"], "lng": df["lng"],
            "portes": lg_[ (cu == 1000) & (lg_ >= 2) ].values,
            "valeur_role": base,
            "valeur_estimee": np.where(ok, np.round(base * np.exp(pl) * sm, -2), np.nan),
            "valeur_estimee_p10": np.where(ok, np.round(base * np.exp(lo) * sm, -2), np.nan),
            "valeur_estimee_p90": np.where(ok, np.round(base * np.exp(hi) * sm, -2), np.nan),
        })
        out["municipalite"] = out["code_mun"].astype(str).str.strip() \
            .str.lstrip("0").map(muns)
        out.to_parquet(f"{ESTIM}/plex_estimes_{Y}.parquet", index=False)
        print(f"[{Y}] {len(out)} plex | médiane {out['valeur_estimee'].median():,.0f} $",
              flush=True)
    print("PREDICT_DONE")


if __name__ == "__main__":
    {"train": train, "predict": predict}[sys.argv[1]]()