#!/usr/bin/env python3
# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca
"""LaTeX tables and number macros for WP23 from results/ (missing inputs skipped)."""
from __future__ import annotations

import json
import re
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))
import numpy as np
import polars as pl

from wp23 import config as C

T = C.ROOT / "paper" / "tables"; T.mkdir(parents=True, exist_ok=True); R = C.RESULTS
MACROS: dict[str, str] = {}
WORDS = {str(i): w for i, w in enumerate(["Zero", "One", "Two", "Three", "Four", "Five", "Six", "Seven", "Eight", "Nine"])}
LABELS = {"naive_mun_median": "Naive: municipal median", "hedonic_ols": "Hedonic OLS", "ridge": "Ridge", "elastic_net": "Elastic net", "random_forest": "Random forest", "xgboost": "XGBoost", "lightgbm": "LightGBM",
          "svr": "SVR (subsample)", "mlp": "Quantile MLP", "hshn_partial": "HSHN (partial pooling)", "hshn_full": "HSHN (full pooling)", "hshn_none": "HSHN (no pooling)",
          "no_change": "No change", "mean_growth": "Mean growth", "momentum": "Momentum", "rw_drift": "Random walk with drift", "seasonal_naive": "Seasonal naive", "ar2_cell": "AR(2), cell", "local_linear_trend": "Local linear trend",
          "ridge_panel": "Ridge (panel)", "lightgbm_panel": "LightGBM (panel)", "mlp_panel": "MLP (panel)", "lstm": "LSTM", "gru": "GRU", "tcn": "TCN", "tft_lite": "TFT-lite", "dsgfn": "DSGFN", "raref": "RAREF", "bhpf": "BHPF", "qanss": "QANSS", "qrefa": "QREFA"}
ORDER_A = ["naive_mun_median", "hedonic_ols", "ridge", "elastic_net", "svr", "random_forest", "xgboost", "lightgbm", "mlp", "hshn_none", "hshn_full", "hshn_partial"]
ORDER_B = ["no_change", "mean_growth", "momentum", "rw_drift", "seasonal_naive", "ar2_cell", "local_linear_trend", "ridge_panel", "lightgbm_panel", "mlp_panel", "lstm", "gru", "tcn", "tft_lite", "bhpf", "raref", "qanss", "dsgfn", "qrefa"]


def mname(n): return "".join(WORDS[c] if c.isdigit() else c for c in re.sub(r"[^A-Za-z0-9]", "", n))
def macro(n, v, fmt="{:,.1f}"):
    MACROS[mname(n)] = "n/a" if v is None or (isinstance(v, float) and np.isnan(v)) else fmt.format(v).replace(",", "{,}")
def esc(s): return str(s).replace("&", "\\&").replace("%", "\\%").replace("_", "\\_").replace("≤", "$\\le$").replace("≥", "$\\ge$").replace("×", "$\\times$")
def ok(*p):
    miss = [str(x) for x in p if not Path(x).exists()]
    if miss: print("skip:", miss)
    return not miss


def table(name, header, rows, caption, label, align=None, notes="", size="\\small", fit=None):
    align = align or "l" + "r" * (len(header) - 1); body = "\n".join(" & ".join(esc(c) for c in r) + " \\\\" for r in rows)
    fit = (len(header) > 7) if fit is None else fit
    tab = f"""\\begin{{tabular}}{{{align}}}
\\toprule
{' & '.join(esc(h) for h in header)} \\\\
\\midrule
{body}
\\bottomrule
\\end{{tabular}}"""
    if fit:  # wide table: scale the tabular to the text width; notes in a full-width parbox (threeparttable mis-measures scaled boxes)
        inner = "\\resizebox{\\textwidth}{!}{%\n" + tab + "}" + (("\n\\par\\vspace{3pt}\\parbox{\\textwidth}{\\footnotesize " + notes + "}") if notes else "")
    else:
        inner = "\\begin{threeparttable}\n" + tab + "\n" + (("\\begin{tablenotes}\\footnotesize\\item " + notes + "\\end{tablenotes}\n") if notes else "") + "\\end{threeparttable}"
    (T / f"{name}.tex").write_text(f"""\\begin{{table}}[htbp]
\\centering{size}
\\caption{{{caption}}}\\label{{{label}}}
{inner}
\\end{{table}}
"""); print("table", name)


def t_data():
    macro("NBenchB", 14, "{:.0f}"); macro("NBenchA", 8, "{:.0f}")
    s = pl.read_parquet(C.WP22_ANALYTIC, columns=["tx_id", "property_class", "region_code", "mun_code", "month", "price", "floor_area_m2", "land_area_m2", "age", "storeys"])
    macro("NSales", s.height, "{:,.0f}"); macro("NMun", s["mun_code"].n_unique(), "{:,.0f}")
    rows = [["Floor area of storeys (m²)", "roll", "continuous; cubic in centred log", "all tasks"], ["Lot area (m²)", "roll", "quadratic in centred log", "all tasks"], ["Year built / age", "roll", "nine age bins; missingness indicator", "all tasks"],
            ["Storeys, building genre, physical link, dwellings (plex)", "roll", "categorical", "Task A"], ["Coordinates, municipality, MRC, region", "feed / SDA V2026-07", "location; nested metric cells (WP22)", "all tasks"],
            ["Assessment in force at sale", "feed", "not used as a predictor (post-dated for Task A test sales would be fine, but excluded for comparability)", "—"],
            ["Municipal / regional market features at origin", "computed from sales ≤ origin", "median log price and quality-stripped level (3/6/12/24 m), momentum, volumes, dispersion", "Task A"],
            ["10 nearest past sales (≤ 12 m)", "computed", "mean quality-stripped log price, mean distance, dispersion", "Task A"],
            ["Real-time SHARP-HPI vintage (13-month history)", "WP22 index_vintages (smoothed mode)", "1/3/6/12-m changes, volatilities, acceleration, counts, s.e.; parent (region, province) momentum", "Task B"],
            ["Bedrooms, bathrooms, parking, basement, condition, renovations, quality, listing, images", "—", "not available", "—"]]
    table("tab_variables", ["Variable", "Source", "Treatment", "Used in"], rows, "Variables available and their use", "tab:variables", align="p{4.2cm}p{3cm}p{5.2cm}p{1.6cm}", size="\\footnotesize")
    cov = s.group_by(["region_code", "property_class"]).agg(pl.len().alias("n")).pivot(on="property_class", index="region_code", values="n").fill_null(0).sort("region_code")
    rows = [[r["region_code"], f"{r.get('house',0):,}".replace(",", "{,}"), f"{r.get('condo',0):,}".replace(",", "{,}"), f"{r.get('plex',0):,}".replace(",", "{,}")] for r in cov.iter_rows(named=True)]
    if ok(R / "task_a" / "timing_trees.csv"):
        ta = pl.concat([pl.read_csv(f) for f in (R / "task_a").glob("timing_*.csv")]); macro("NOriginsA", ta["origin"].n_unique(), "{:.0f}"); macro("NTestA", ta.group_by("origin").agg(pl.col("n_test").first())["n_test"].sum(), "{:,.0f}")
    if ok(C.PROCESSED / "task_b_panel.parquet"):
        pb = pl.read_parquet(C.PROCESSED / "task_b_panel.parquet"); macro("NCellsB", pb.select(["class", "level", "geo_id"]).unique().height, "{:.0f}"); macro("NRowsB", pb.height, "{:,.0f}"); macro("NOriginsB", pb["origin"].n_unique(), "{:.0f}")
        for h in C.HORIZONS: macro(f"NOriginsBh{h}", pb.filter(pl.col("h") == h)["origin"].n_unique(), "{:.0f}")
    table("tab_coverage", ["Region", "Houses", "Condominiums", "Plexes"], rows, "Headline sales by administrative region and class (WP22 sample reused)", "tab:coverage", size="\\footnotesize")
    # selection: sold vs stock
    if ok(C.WP22_STOCK):
        st = pl.read_parquet(C.WP22_STOCK, columns=["cubf", "nb_logements", "lien_physique", "annee_construction", "aire_etages_m2", "superficie_terrain_m2"]).filter(pl.col("cubf") == 1000)
        st = st.with_columns(pl.when(pl.col("nb_logements") >= 6).then(None).when(pl.col("nb_logements").is_between(2, 5)).then(pl.lit("plex")).when((pl.col("nb_logements") <= 1) & (pl.col("lien_physique") == "Intégré")).then(pl.lit("condo")).when(pl.col("nb_logements") <= 1).then(pl.lit("house")).otherwise(None).alias("cls")).filter(pl.col("cls").is_not_null())
        rows = []
        for cls in ["house", "condo", "plex"]:
            a = s.filter(pl.col("property_class") == cls); b = st.filter(pl.col("cls") == cls)
            rows.append([cls, f"{a['floor_area_m2'].median():.0f}", f"{b['aire_etages_m2'].median():.0f}", f"{a['age'].median():.0f}", f"{2026 - b['annee_construction'].median():.0f}", f"{a['land_area_m2'].median():.0f}", f"{b['superficie_terrain_m2'].median():.0f}", f"{a.height:,}".replace(",", "{,}"), f"{b.height:,}".replace(",", "{,}")])
            macro(f"SelFa{cls}Sold", a["floor_area_m2"].median(), "{:.0f}"); macro(f"SelFa{cls}Stock", b["aire_etages_m2"].median(), "{:.0f}"); macro(f"SelAge{cls}Sold", a["age"].median(), "{:.0f}"); macro(f"SelAge{cls}Stock", 2026 - b["annee_construction"].median(), "{:.0f}")
        table("tab_selection", ["Class", "Floor area sold", "Floor area stock", "Age sold", "Age stock", "Lot sold", "Lot stock", "Sales", "Stock units"], rows, "Transaction selection: medians of sold dwellings (2021--2026) versus the assessment-roll stock", "tab:selection", size="\\footnotesize",
              notes="Stock = residential use code 1000 in the 2026 roll, classified with the WP22 rules; age of the stock as of 2026.")
    rows = [["HSHN", "embeddings 8 (geo) / 4 (h); encoder 256-128 GELU, dropout 0.1; AdamW, one-cycle lr $2\\times10^{-3}$, 8 epochs, batch 4,096; municipality weight decay $5\\times10^{-3}$"],
            ["MLP (A)", "256-128 GELU, dropout 0.1; pinball loss on 5 quantiles; same optimiser"], ["LightGBM (A)", "800 trees, lr 0.04, 63 leaves, min child 30, subsample 0.8; quantile models 400 trees"], ["XGBoost (A)", "600 trees, lr 0.05, depth 8, subsample 0.8"], ["Random forest (A)", "150 trees, depth 18, min leaf 5, max features 0.5"], ["SVR (A)", "RBF, C=3, $\\epsilon$=0.05, 15,000-sale subsample"],
            ["Panel MLP / LSTM / GRU / TCN / TFT-lite (B)", "hidden 64 (MLP) / 32 (RNN, attention) / 16 channels (TCN); 60 epochs, batch 512, AdamW lr $2\\times10^{-3}$, weight decay $10^{-4}$; pinball loss"],
            ["DSGFN (B)", "node encoder 32; one attention layer; prior adjacency: Gaussian distance kernel (60 km) $\\times$ co-movement; 60 epochs"], ["RAREF (B)", "3-state Gaussian HMM (60 EM iterations) on provincial real-time growth; 3 ridge experts ($\\alpha$=3)"],
            ["BHPF (B)", "conjugate Gaussian hierarchy; 15 variance-component iterations; priors $\\tau_b^2$=1 (slopes)"], ["QANSS (B)", "32 tanh units; 200 Adam steps, lr $5\\times10^{-3}$; Gaussian likelihood"], ["QREFA (B)", "NNLS stacking of QANSS, BHPF, DSGFN (40 epochs), RAREF on training origins"],
            ["LightGBM (B)", "300 trees, lr 0.03, 15 leaves, min child 40; quantile models"], ["Ridge (B)", "$\\alpha$=3 on standardised features"]]
    table("tab_hyper", ["Model", "Hyper-parameters (fixed in advance, never tuned on test origins)"], rows, "Hyper-parameters", "tab:hyper", align="p{3.4cm}p{11.5cm}", size="\\footnotesize")
    rows = [["Hedonic coefficients", "WP22 stage 1 re-estimated at every vintage; Task A features never use coefficients"], ["Geographic clustering / location cells", "WP22 nested grid rebuilt from data ≤ origin at every vintage; Task A uses raw coordinates"],
            ["Normalisation", "standardisation constants from training rows only; index levels enter as changes"], ["Feature selection", "fixed feature lists, no selection on test data"], ["Imputation", "WP22 conditional medians (no price information); indicators time-stable"],
            ["Target encoding", "none; one-hot with training categories"], ["Spatial graph", "DSGFN adjacency recomputed at each origin from centroids and the real-time history ≤ origin"], ["Hyper-parameter tuning", "fixed a priori (Table~\\ref{tab:hyper})"],
            ["Index revisions", "predictors = vintage published for the origin; target = final series (revision component reported)"], ["Partial last month", "October 2026 excluded; last target September 2026"]]
    table("tab_leakage", ["Item", "Treatment"], rows, "Leakage audit", "tab:leakage", align="p{4cm}p{11cm}", size="\\footnotesize")


def scores_table(tag, order, caption, label, hs, ycol_desc):
    f = R / f"eval_{tag}" / "scores.csv"
    if not ok(f): return
    d = pl.read_csv(f, infer_schema_length=None).filter(pl.col("group") == "all"); mcs = pl.read_csv(R / f"eval_{tag}" / "mcs.csv", infer_schema_length=None) if (R / f"eval_{tag}" / "mcs.csv").exists() else None
    tests = pl.read_csv(R / f"eval_{tag}" / "tests_vs_reference.csv", infer_schema_length=None) if (R / f"eval_{tag}" / "tests_vs_reference.csv").exists() else None
    rows = []
    for m in order:
        if m not in d["model"].to_list(): continue
        r = [LABELS.get(m, m)]
        for h in hs:
            x = d.filter((pl.col("model") == m) & (pl.col("h").cast(pl.Utf8) == str(h)))
            if x.height == 0: r += ["—", "—"]; continue
            mae = 100 * x["mae"][0]; star = ""
            if tests is not None:
                tt = tests.filter((pl.col("model") == m) & (pl.col("h").cast(pl.Utf8) == str(h)))
                if tt.height and tt["dm_p"][0] == tt["dm_p"][0]:
                    star = "***" if tt["dm_p"][0] < 0.01 else "**" if tt["dm_p"][0] < 0.05 else "*" if tt["dm_p"][0] < 0.10 else ""
            inm = ""
            if mcs is not None:
                mm = mcs.filter((pl.col("model") == m) & (pl.col("h").cast(pl.Utf8) == str(h)))
                if mm.height and mm["in_mcs_90"][0]: inm = "$^{\\dagger}$"
            r += [f"{mae:.2f}{star}{inm}", f"{x['mase'][0]:.2f}" if x["mase"][0] == x["mase"][0] else "—"]
            macro(f"{tag}MAE{m}h{h}", mae, "{:.2f}"); macro(f"{tag}MASE{m}h{h}", x["mase"][0], "{:.2f}"); macro(f"{tag}CRPS{m}h{h}", 100 * x["crps"][0] if x["crps"][0] == x["crps"][0] else None, "{:.2f}"); macro(f"{tag}Cov{m}h{h}", 100 * x["cov90"][0] if x["cov90"][0] == x["cov90"][0] else None, "{:.0f}")
        rows.append(r)
    header = ["Model"] + [f"{h}: MAE" if i % 2 == 0 else f"{h}: MASE" for h in hs for i in range(2)]
    table(f"tab_scores_{tag}", header, rows, caption, label, size="\\scriptsize",
          notes=f"{ycol_desc} MAE in log points ×100; MASE relative to the naive reference. Stars: Diebold--Mariano test of equal squared-error loss against the reference on per-origin mean losses (* 10\\%, ** 5\\%, *** 1\\%); $\\dagger$ = in the 90\\% model confidence set of the horizon.")
    rows = []
    for m in order:
        if m not in d["model"].to_list(): continue
        r = [LABELS.get(m, m)]
        for h in hs:
            x = d.filter((pl.col("model") == m) & (pl.col("h").cast(pl.Utf8) == str(h)))
            r += [f"{100*x['crps'][0]:.2f}" if x.height and x["crps"][0] == x["crps"][0] else "—", f"{100*x['cov90'][0]:.0f}" if x.height and x["cov90"][0] == x["cov90"][0] else "—", f"{100*x['width90'][0]:.1f}" if x.height and x["width90"][0] == x["width90"][0] else "—"]
        rows.append(r)
    header = ["Model"] + [f"{h}: {k}" for h in hs for k in ("CRPS", "cov90", "width90")]
    table(f"tab_prob_{tag}", header, rows, caption.replace("Point", "Probabilistic").replace("accuracy", "accuracy (CRPS ×100, 90\\% interval coverage \\% and width ×100)"), label.replace("scores", "prob"), size="\\scriptsize")


def t_results():
    f = R / "eval_a" / "scores.csv"
    if f.exists():
        d = pl.read_csv(f, infer_schema_length=None).filter(pl.col("group") == "all")
        for m in ["lightgbm", "hshn_partial", "hedonic_ols"]:
            a = d.filter((pl.col("model") == m) & (pl.col("h") == "1-3")); b = d.filter((pl.col("model") == m) & (pl.col("h") == "7-12"))
            if a.height and b.height: macro(f"aMAEDelta{m}", 100 * (b["mae"][0] - a["mae"][0]), "{:.1f}")
    scores_table("a", ORDER_A, "Task A — Point accuracy of future sale-price forecasts by horizon band (months after the origin)", "tab:scores_a", ["1-3", "4-6", "7-12"], "Target: log sale price.")
    scores_table("b", ORDER_B, "Task B — Point accuracy of quality-adjusted index-change forecasts by horizon (months)", "tab:scores_b", [1, 3, 6, 12], "Target: final SHARP-HPI log change.")
    for tag, grp, lab in [("b", "level", "level"), ("b", "class", "class"), ("b", "size_level", "cell size"), ("a", "mun_size", "municipality size"), ("a", "property_class", "class"), ("a", "segment", "price segment")]:
        f = R / f"eval_{tag}" / "scores.csv"
        if not ok(f): continue
        d = pl.read_csv(f, infer_schema_length=None).filter(pl.col("group") == grp)
        if d.height == 0: continue
        groups = sorted(d[grp].unique().to_list()); order = ORDER_B if tag == "b" else ORDER_A; hs = [3] if tag == "b" else ["4-6"]
        rows = []
        for m in order:
            if m not in d["model"].to_list(): continue
            r = [LABELS.get(m, m)]
            for g in groups:
                x = d.filter((pl.col("model") == m) & (pl.col(grp) == g) & (pl.col("h").cast(pl.Utf8) == str(hs[0])))
                r.append(f"{100*x['mae'][0]:.2f}" if x.height else "—")
            rows.append(r)
        table(f"tab_het_{tag}_{grp}", ["Model"] + [str(g) for g in groups], rows, f"Task {tag.upper()} — MAE (×100) by {lab}, horizon {hs[0]}", f"tab:het_{tag}_{grp}", size="\\scriptsize")


def t_ablation():
    f = R / "ablation_b" / "ablation_summary.csv"
    if not ok(f): return
    d = pl.read_csv(f, infer_schema_length=None).filter(pl.col("level") == "municipality")
    lab = {"dsgfn": "DSGFN (full)", "dsgfn_no_distance": "DSGFN without distance kernel", "dsgfn_no_comovement": "DSGFN without co-movement", "dsgfn_no_graph": "DSGFN with uninformative prior graph", "qanss": "QANSS (full)", "qanss_no_parent": "QANSS without parent drifts",
           "qrefa": "QREFA (NNLS stack)", "equal_weight_combo": "Equal-weight combination of the four", "qrefa_without_qanss": "QREFA without QANSS", "qrefa_without_bhpf": "QREFA without BHPF", "qrefa_without_dsgfn": "QREFA without DSGFN", "qrefa_without_raref": "QREFA without RAREF",
           "ridge_panel": "Ridge panel (full features)", "ridge_no_momentum": "Ridge without own momentum", "ridge_no_parent": "Ridge without parent (region/province) features", "ridge_no_volume": "Ridge without volume features"}
    rows = []
    for m in lab:
        r = [lab[m]]
        for h in (1, 3, 6):
            x = d.filter((pl.col("model") == m) & (pl.col("h") == h))
            if x.height == 0: r += ["—", "—"]; continue
            diff = x["mae_diff_vs_ref"][0]; lo, hi = x["diff_lo"][0], x["diff_hi"][0]
            sig = "" if diff is None or diff != diff else ("*" if (lo > 0 or hi < 0) else "")
            r += [f"{100*x['mae'][0]:.2f}", (f"{100*diff:+.2f}{sig}" if diff == diff and diff is not None else "ref")]
            macro(f"Abl{m}h{h}", 100 * x["mae"][0], "{:.2f}")
            if diff == diff and diff is not None: macro(f"AblDiff{m}h{h}", 100 * diff, "{:+.2f}")
        rows.append(r)
    table("tab_ablation_b", ["Variant", "h=1: MAE", "vs ref", "h=3: MAE", "vs ref", "h=6: MAE", "vs ref"], rows, "Task B ablations at the municipality level: MAE (×100) and difference to the reference variant (per-origin mean, * = 95\\% stationary-bootstrap interval excludes zero)", "tab:ablation_b", size="\\scriptsize",
          notes="Reference variants: DSGFN, QANSS, QREFA and the full ridge panel. Positive differences mean the ablated variant is worse.", fit=True)
    a = R / "ablation_b" / "dsgfn_attention_summary.csv"
    if a.exists():
        x = pl.read_csv(a, infer_schema_length=None)
        for r in x.iter_rows(named=True):
            macro(f"AttTop{r['h']}", r["top_weight"], "{:.2f}"); macro(f"AttDist{r['h']}", r["median_dist_top_km"], "{:.0f}"); macro(f"AttSameReg{r['h']}", 100 * r["same_region"], "{:.0f}"); macro(f"AttSameMrc{r['h']}", 100 * r["same_mrc"], "{:.0f}")


def t_robust():
    f = R / "robustness_b" / "summary.csv"
    if not ok(f): return
    d = pl.read_csv(f, infer_schema_length=None).filter(pl.col("level") == "municipality")
    variants = ["expanding", "rolling_24", "min_n_10", "winsorised", "seed_1", "seed_2", "bandwidth"]; labv = {"expanding": "Expanding (baseline)", "rolling_24": "Rolling 24 origins", "min_n_10": "Cells with $n \\ge 10$", "winsorised": "Winsorised target", "seed_1": "Seed 1", "seed_2": "Seed 2", "bandwidth": "Bandwidth 30 / 120 km"}
    models = ["no_change", "momentum", "ar2_cell", "ridge_panel", "qanss", "tft_lite", "dsgfn", "dsgfn_bw30", "dsgfn_bw120"]
    rows = []
    for m in models:
        r = [LABELS.get(m, m.replace("dsgfn_bw", "DSGFN, bandwidth ").replace("30", "30 km").replace("120", "120 km") if m.startswith("dsgfn_bw") else LABELS.get(m, m))]
        for v in variants:
            x = d.filter((pl.col("model") == m) & (pl.col("variant") == v)); r.append(f"{100*x['mae'][0]:.2f}" if x.height else "—")
            if x.height: macro(f"Rob{m}{v}", 100 * x["mae"][0], "{:.2f}"); macro(f"RobCov{m}{v}", 100 * x["cov90"][0], "{:.0f}")
        rows.append(r)
    table("tab_robust_b", ["Model"] + [labv[v] for v in variants], rows, "Task B robustness at the municipality level, horizon 3: MAE (×100) under alternative designs", "tab:robust_b", size="\\scriptsize",
          notes="Rolling: models estimated on the 24 most recent origins only. Seeds: neural models re-estimated with different random seeds (the baseline is seed 0). Winsorised: targets clipped at the 1st and 99th percentiles in training and test.", fit=True)


def t_interpret():
    f = R / "interpret" / "shap_importance_wide.csv"
    if not ok(f): return
    d = pl.read_csv(f, infer_schema_length=None).head(15); cols = [c for c in d.columns if c not in ("feature", "mean")]
    lab = {"knn_q": "10 nearest past sales: quality-stripped log price", "lfa": "log floor area (centred)", "age_f": "age", "origin_t": "origin month (time index)", "lla": "log lot area (centred)", "mun_q_12": "municipality: quality-stripped level, 12 m",
           "y_m": "northing", "x_m": "easting", "mun_medlp_3": "municipality: median log price, 3 m", "knn_dist_km": "distance to the 10 nearest sales", "knn_q_sd": "dispersion of the 10 nearest sales", "mun_n_3": "municipality: sales, 3 m", "mun_n_12": "municipality: sales, 12 m",
           "reg_q_12": "region: quality-stripped level, 12 m", "reg_q_3": "region: level, 3 m", "lfa2": "log floor area squared", "mun_mom_3_12": "municipality: momentum 3 vs 12 m", "reg_mom_3_12": "region: momentum", "mun_vol_ratio": "municipality: volume ratio"}
    rows = [[lab.get(r["feature"], r["feature"])] + [f"{100*r[c]:.1f}" for c in cols] + [f"{100*r['mean']:.1f}"] for r in d.iter_rows(named=True)]
    table("tab_shap", ["Feature"] + cols + ["Mean"], rows, "Task A — mean absolute SHAP value (×100, log-price units) of the LightGBM model at three origins, 4,000 test sales each", "tab:shap", size="\\footnotesize",
          notes="Importance is predictive, not causal. Origin-month importance declines as the training history lengthens.", fit=True)
    for r in d.iter_rows(named=True):
        macro(f"Shap{r['feature']}", 100 * r["mean"], "{:.1f}")
        for c in cols: macro(f"Shap{r['feature']}{c}", 100 * r[c], "{:.1f}")
    a = pl.read_csv(R / "interpret" / "ale.csv", infer_schema_length=None).filter(pl.col("origin") == "2024-03")
    for feat in ["mun_mom_3_12", "knn_q", "mun_vol_ratio"]:
        x = a.filter(pl.col("feature") == feat).sort("bin"); macro(f"ALE{feat}Range", 100 * (x["ale"].max() - x["ale"].min()), "{:.1f}"); macro(f"ALE{feat}Top", 100 * x["ale"][-1], "{:+.1f}"); macro(f"ALE{feat}Bottom", 100 * x["ale"][0], "{:+.1f}")


def t_extra():
    X = R / "extra"
    if not ok(X / "tests_holm_b.csv"): return
    t = pl.read_csv(X / "tests_holm_b.csv", infer_schema_length=None)
    rows, rows2 = [], []
    for m in ORDER_B:
        if m == "no_change": continue
        r, r2 = [LABELS.get(m, m)], [LABELS.get(m, m)]
        for h in (1, 3, 6, 12):
            x = t.filter((pl.col("model") == m) & (pl.col("h") == h))
            has = x.height and x["dm_stat"][0] is not None
            r += [f"{x['dm_stat'][0]:+.2f}" if has else "—", f"{x['dm_p_holm'][0]:.2f}" if has else "—"]
            r2 += [f"{x['dm_p'][0]:.3f}" if has else "—", f"{x['gw_p'][0]:.3f}" if has and x["gw_p"][0] is not None else "—"]
            if has: macro(f"HolmP{m}h{h}", x["dm_p_holm"][0], "{:.2f}"); macro(f"DMP{m}h{h}", x["dm_p"][0], "{:.3f}"); macro(f"GWP{m}h{h}", x["gw_p"][0], "{:.3f}")
        rows.append(r); rows2.append(r2)
    table("tab_tests_b", ["Model"] + [f"h={h}: {k}" for h in (1, 3, 6, 12) for k in ("DM", "Holm $p$")], rows, "Task B — Diebold--Mariano statistics against the no-change forecast (negative = better than no change) and Holm-adjusted $p$-values across the 18 models of each horizon", "tab:tests_b", size="\\scriptsize", fit=True,
          notes="HAC variance with $h-1$ Bartlett lags and the Harvey--Leybourne--Newbold correction; $t_{n-1}$ reference distribution; Holm step-down over the 18 two-sided $p$-values of each horizon.")
    table("tab_gw_b", ["Model"] + [f"h={h}: {k}" for h in (1, 3, 6, 12) for k in ("DM $p$", "GW $p$")], rows2, "Task B — unadjusted Diebold--Mariano $p$-values and Giacomini--White conditional-predictive-ability $p$-values against the no-change forecast", "tab:gw_b", size="\\scriptsize", fit=True,
          notes="The Giacomini--White test regresses the loss differential on a constant and its own lag; a rejection means the relative accuracy of the two forecasts is predictable from the past.")
    c = pl.read_csv(X / "computational_cost.csv", infer_schema_length=None)
    rows = [[r["task"], LABELS.get(r["model"], r["model"]), f"{r['mean_seconds_per_fit']:.1f}", f"{r['total_seconds']:.0f}", str(r["fits"]), f"{r['mean_train_rows']:,.0f}".replace(",", "{,}")] for r in c.iter_rows(named=True)]
    table("tab_cost", ["Task", "Model", "Seconds per fit", "Total seconds", "Fits", "Mean training rows"], rows, "Computational cost (wall-clock on a 32-core Mac Studio M3 Ultra; tree models 10 threads, neural models 8 threads)", "tab:cost", size="\\scriptsize")
    for r in c.iter_rows(named=True): macro(f"Cost{r['task']}{r['model']}", r["mean_seconds_per_fit"], "{:.1f}")
    y = pl.read_csv(X / "error_by_year_b.csv", infer_schema_length=None).filter(pl.col("h") == 3)
    yrs = sorted(y["year"].unique().to_list()); rows = []
    for m in KEYS_B:
        x = y.filter(pl.col("model") == m)
        if x.height == 0: continue
        rows.append([LABELS.get(m, m)] + [f"{100*x.filter(pl.col('year') == yr)['mae'][0]:.2f} ({x.filter(pl.col('year') == yr)['mase'][0]:.2f})" if x.filter(pl.col("year") == yr).height else "—" for yr in yrs])
        for yr in yrs:
            z = x.filter(pl.col("year") == yr)
            if z.height: macro(f"Year{m}{yr}", z["mase"][0], "{:.2f}")
    table("tab_year_b", ["Model"] + [f"Origins in {yr}" for yr in yrs], rows, "Task B, horizon 3 — MAE (×100) and MASE (in parentheses) by year of the forecast origin", "tab:year_b", size="\\scriptsize", notes="2023 covers origins from October to December only; 2026 covers January to June.")
    g = pl.read_csv(X / "error_by_regime_b.csv", infer_schema_length=None).filter(pl.col("h") == 3); regs = ["rising", "flat", "falling"]; rows = []
    for m in KEYS_B:
        x = g.filter(pl.col("model") == m)
        if x.height == 0: continue
        rows.append([LABELS.get(m, m)] + [f"{100*x.filter(pl.col('regime') == rg)['mae'][0]:.2f} ({x.filter(pl.col('regime') == rg)['mase'][0]:.2f})" if x.filter(pl.col("regime") == rg).height else "—" for rg in regs])
        for rg in regs:
            z = x.filter(pl.col("regime") == rg)
            if z.height: macro(f"Regime{m}{rg}", z["mase"][0], "{:.2f}")
    ns = {rg: (int(g.filter((pl.col("regime") == rg) & (pl.col("model") == "no_change"))["origins"][0]) if g.filter((pl.col("regime") == rg) & (pl.col("model") == "no_change")).height else 0) for rg in regs}
    table("tab_regime_b", ["Model"] + [f"{rg} ({ns[rg]} origins)" for rg in regs], rows, "Task B, horizon 3 — MAE (×100) and MASE (in parentheses) by ex-post provincial regime of the target window (rising: provincial house index change above +2 log points over the horizon; falling: below $-2$; flat otherwise)", "tab:regime_b", size="\\scriptsize")
    for rg in regs: macro(f"RegimeOrigins{rg}", ns[rg], "{:.0f}")
    d = pl.read_csv(X / "directional_b.csv", infer_schema_length=None); rows = []
    for m in KEYS_B:
        if m == "no_change": continue
        r = [LABELS.get(m, m)]
        for h in (1, 3, 6, 12):
            x = d.filter((pl.col("model") == m) & (pl.col("h") == h))
            r += [f"{100*x['directional_accuracy'][0]:.0f}" if x.height else "—", f"{100*x['directional_accuracy_at_turning_points'][0]:.0f}" if x.height and x["directional_accuracy_at_turning_points"][0] is not None else "—"]
            if x.height: macro(f"Dir{m}h{h}", 100 * x["directional_accuracy"][0], "{:.0f}"); macro(f"DirTP{m}h{h}", 100 * x["directional_accuracy_at_turning_points"][0] if x["directional_accuracy_at_turning_points"][0] is not None else None, "{:.0f}")
        rows.append(r)
    table("tab_direction_b", ["Model"] + [f"h={h}: {k}" for h in (1, 3, 6, 12) for k in ("all", "turning pts")], rows, "Task B — directional accuracy (\\% of forecasts with the sign of the realised change, changes above 0.5 log points) overall and at turning points (origins where the realised change reverses sign)", "tab:direction_b", size="\\scriptsize",
          notes="The no-change forecast has no sign; momentum has the sign of the last three-month change.", fit=True)
    dc = pl.read_csv(X / "decline_b.csv", infer_schema_length=None); rows = []
    for m in KEYS_B:
        r = [LABELS.get(m, m)]
        for h in (1, 3, 6, 12):
            x = dc.filter((pl.col("model") == m) & (pl.col("h") == h))
            r += [f"{x['brier_decline'][0]:.3f}" if x.height else "—", f"{x['brier_skill_decline'][0]:+.2f}" if x.height else "—"]
            if x.height: macro(f"BrierSkill{m}h{h}", x["brier_skill_decline"][0], "{:+.2f}"); macro(f"ShareDeclinesh{h}", 100 * x["share_declines"][0], "{:.0f}")
        rows.append(r)
    table("tab_decline_b", ["Model"] + [f"h={h}: {k}" for h in (1, 3, 6, 12) for k in ("Brier", "skill")], rows, "Task C — probability of an index decline over the horizon: Brier score and Brier skill score relative to the climatological frequency of declines", "tab:decline_b", size="\\scriptsize",
          notes="Probabilities are read off the five-quantile forecasts by linear interpolation of the quantile function. Skill = $1 - $ Brier / Brier of the constant climatological probability.")
    cb = pl.read_csv(X / "combinations_b.csv", infer_schema_length=None); lab = {"eq_lgbm_qanss": "Equal weights: LightGBM + QANSS", "eq_lgbm_qanss_ridge_ar2": "Equal weights: LightGBM + QANSS + ridge + AR(2)", "eq_lgbm_ar2": "Equal weights: LightGBM + AR(2)", "eq_all_with_nochange": "Equal weights: four models + no change"}
    rows = []
    for name in ["lightgbm_panel", "qanss", "ridge_panel", "ar2_cell", "qrefa", "eq_lgbm_qanss", "eq_lgbm_ar2", "eq_lgbm_qanss_ridge_ar2", "eq_all_with_nochange"]:
        x = cb.filter(pl.col("combination") == name); r = [lab.get(name, LABELS.get(name, name))]
        for h in (1, 3, 6, 12):
            z = x.filter(pl.col("h") == h); r.append(f"{z['mase'][0]:.3f}" if z.height else "—")
            if z.height: macro(f"Combo{name}h{h}", z["mase"][0], "{:.3f}")
        rows.append(r)
    table("tab_combo_b", ["Forecast", "h=1", "h=3", "h=6", "h=12"], rows, "Task B — MASE of equal-weight combinations versus their components (same cells and origins)", "tab:combo_b", size="\\footnotesize")
    md = pl.read_csv(X / "municipal_mase_summary_b.csv", infer_schema_length=None)
    rows = [[LABELS.get(r["model"], r["model"]), f"{r['median_mase']:.2f}", f"{r['p10']:.2f}", f"{r['p90']:.2f}", f"{100*r['share_beating_no_change']:.0f}", str(r["cells"])] for r in md.iter_rows(named=True)]
    table("tab_mun_dist_b", ["Model", "Median MASE", "P10", "P90", "% of cells beating no change", "Cells"], rows, "Task B, horizon 3 — distribution of the MASE across the 192 municipal cells (64 municipalities × 3 classes)", "tab:mun_dist_b", size="\\footnotesize")
    for r in md.iter_rows(named=True): macro(f"MunShare{r['model']}", 100 * r["share_beating_no_change"], "{:.0f}"); macro(f"MunMedMASE{r['model']}", r["median_mase"], "{:.2f}"); macro(f"MunPninety{r['model']}", r["p90"], "{:.2f}")
    de = pl.read_csv(X / "dollar_errors_a.csv", infer_schema_length=None)
    rows = [[r["property_class"], LABELS.get(r["model"], r["model"]), f"{100*r['mae_log']:.1f}", f"{r['median_price']:,.0f}".replace(",", "{,}"), f"{r['median_abs_dollar_error']:,.0f}".replace(",", "{,}")] for r in de.iter_rows(named=True)]
    table("tab_dollar_a", ["Class", "Model", "MAE (log ×100)", "Median price (\\$)", "Median absolute error (\\$)"], rows, "Task A, horizons 1--3 months — economic magnitude of forecast errors", "tab:dollar_a", size="\\footnotesize")
    for r in de.iter_rows(named=True): macro(f"Dollar{r['model']}{r['property_class']}", r["median_abs_dollar_error"], "{:,.0f}")
    sr = pl.read_csv(X / "selection_rates.csv", infer_schema_length=None).filter(pl.col("quintile").is_not_null()); rows = []
    for cls in ["house", "condo", "plex"]:
        for dim, lab2 in [("fa_q", "floor area"), ("age_q", "age")]:
            x = sr.filter((pl.col("class") == cls) & (pl.col("dimension") == dim)).sort("quintile"); rows.append([f"{cls}, by {lab2} quintile (1 = smallest / newest)"] + [f"{v:.1f}" for v in x["annual_sale_rate_pct"].to_list()])
            for q, v in zip(x["quintile"].to_list(), x["annual_sale_rate_pct"].to_list()): macro(f"SaleRate{cls}{q}", v, "{:.1f}")
    table("tab_selection_rates", ["Class and dimension", "Q1", "Q2", "Q3", "Q4", "Q5"], rows, "Transaction selection — annual sale rate of the roll stock (\\% of units sold per year, 2021--2026) by quintile of floor area and of age", "tab:selection_rates", size="\\footnotesize")


KEYS_B = ["no_change", "mean_growth", "momentum", "ar2_cell", "ridge_panel", "lightgbm_panel", "tft_lite", "bhpf", "raref", "qanss", "dsgfn", "qrefa"]


def t_simulation():
    f = R / "simulation" / "summary.csv"
    if not ok(f): return
    d = pl.read_csv(f, infer_schema_length=None); scen = ["rw", "momentum", "regimes", "spillover", "sparse"]; scen = [s for s in scen if s in d["scenario"].to_list()]
    for lvl in ["municipality", "region", "province"]:
        rows = []
        for m in ORDER_B:
            x = d.filter((pl.col("model") == m) & (pl.col("level") == lvl) & (pl.col("h") == 3))
            if x.height == 0: continue
            r = [LABELS.get(m, m)]
            for s in scen:
                y = x.filter(pl.col("scenario") == s); r.append(f"{y['rel_mae'][0]:.2f} ({y['rel_mae_mcse'][0]:.2f})" if y.height else "—")
                if y.height: macro(f"Sim{lvl}{m}{s}", 100 * y["mae"][0], "{:.2f}"); macro(f"SimRel{lvl}{m}{s}", y["rel_mae"][0], "{:.2f}"); macro(f"SimCov{lvl}{m}{s}", 100 * y["cov90"][0], "{:.0f}")
            rows.append(r)
        table(f"tab_sim_{lvl}", ["Model"] + scen, rows, f"Simulation — MAE of 3-month index-change forecasts relative to the no-change forecast at the {lvl} level, by scenario (Monte Carlo s.e. over replications)", f"tab:sim_{lvl}", size="\\scriptsize",
              notes="rw: random walk with a constant drift; momentum: AR(1) drifts at every level; regimes: Markov-switching provincial drift; spillover: municipal drift follows the lagged regional drift; sparse: momentum with 35\\% of the sales. Values below 1 beat the no-change forecast. Classic and tree models: 10 replications; neural models: 6.")


if __name__ == "__main__":
    for fn in [t_data, t_results, t_ablation, t_robust, t_interpret, t_extra, t_simulation]:
        try: fn()
        except Exception as ex:
            import traceback; traceback.print_exc(); print("FAILED", fn.__name__)
    (T / "numbers.tex").write_text("% generated by scripts/50_tables.py\n" + "\n".join(f"\\newcommand{{\\{k}}}{{{v}}}" for k, v in sorted(MACROS.items())) + "\n")
    used = set()
    for f in list((C.ROOT / "paper" / "sections").glob("*.tex")) + list((C.ROOT / "paper" / "appendix").glob("*.tex")):
        used |= set(re.findall(r"\\([A-Z][A-Za-z]+)(?=[^A-Za-z]|$)", f.read_text()))
    skip = {"E", "Var", "DeclareMathOperator", "LARGE", "Large", "Delta", "WPabstract", "N"}
    missing = sorted(u for u in used - set(MACROS) - skip if not u.startswith("WP"))
    (T / "numbers_missing.tex").write_text("% must be empty in the final version\n" + "\n".join(f"\\providecommand{{\\{k}}}{{\\textbf{{??{k}??}}}}" for k in missing) + "\n")
    print(len(MACROS), "macros;", len(missing), "missing:", missing)