"""Phase F — figures generated only from pipeline outputs (data/metadata, data/final, results/econometrics).
Writes results/figures/*.{pdf,png} and copies PDFs to paper/figures/. Run: python scripts/visualization/make_figures.py"""
import json, shutil
from pathlib import Path
import matplotlib; matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import FancyBboxPatch
import numpy as np, pandas as pd, yaml

ROOT = Path(__file__).resolve().parents[2]; META = ROOT / "data" / "metadata"; FINAL = ROOT / "data" / "final"; ECO = ROOT / "results" / "econometrics"
FIG = ROOT / "results" / "figures"; PFIG = ROOT / "paper" / "figures"; FIG.mkdir(parents=True, exist_ok=True); PFIG.mkdir(parents=True, exist_ok=True)
PAL = ["#2a78d6", "#eb6834", "#1baf7a", "#eda100", "#e87ba4", "#008300", "#6250d6", "#e34948"]; INK = "#0b0b0b"; INK2 = "#52514e"; GRID = "#e6e5e1"
plt.rcParams.update({"font.family": "serif", "font.size": 9, "axes.titlesize": 9.5, "axes.labelsize": 9, "legend.fontsize": 8, "xtick.labelsize": 8, "ytick.labelsize": 8,
                     "axes.edgecolor": INK2, "axes.grid": True, "grid.color": GRID, "grid.linewidth": 0.6, "axes.spines.top": False, "axes.spines.right": False,
                     "lines.linewidth": 1.6, "legend.frameon": False, "figure.dpi": 150, "savefig.dpi": 300, "axes.prop_cycle": matplotlib.cycler(color=PAL)})
LABELS = yaml.safe_load((ROOT / "scripts" / "harmonization" / "taxonomy.yaml").read_text())
LAB = {leaf: v["label"] for grp in LABELS.values() for leaf, v in grp.items()}; GRP = {leaf: g for g, grp in LABELS.items() for leaf in grp}
FAMLAB = {"pct_premium": "% premium (binary attribute, log price)", "semi_elasticity": "% per unit (log price, level X)", "elasticity": "elasticity (log-log)"}


def save(fig, name):
    fig.tight_layout(); fig.savefig(FIG / f"{name}.pdf", bbox_inches="tight"); fig.savefig(FIG / f"{name}.png", bbox_inches="tight"); shutil.copy(FIG / f"{name}.pdf", PFIG / f"{name}.pdf"); plt.close(fig); print("figure", name)


def fig01_prisma():
    c = json.load(open(META / "prisma_counts.json"))
    fig, ax = plt.subplots(figsize=(7.4, 7.2)); ax.axis("off"); ax.set_xlim(0, 10); ax.set_ylim(0, 12)
    def box(x, y, w, h, text, col="#eef3fb"):
        ax.add_patch(FancyBboxPatch((x, y), w, h, boxstyle="round,pad=0.02,rounding_size=0.08", fc=col, ec=INK2, lw=0.8)); ax.text(x + w / 2, y + h / 2, text, ha="center", va="center", fontsize=7.6, color=INK)
    def arrow(x1, y1, x2, y2): ax.annotate("", (x2, y2), (x1, y1), arrowprops=dict(arrowstyle="->", color=INK2, lw=0.8))
    src = c["by_source"]
    box(0.3, 10.6, 5.4, 1.1, f"Records identified from databases (n = {c['records_identified_raw']:,})\nSemantic Scholar {src.get('semanticscholar_bulk',0):,} · citations {src.get('semanticscholar_citations',0):,} · Crossref {src.get('crossref',0):,}\narXiv {src.get('arxiv',0):,} · OpenAlex {src.get('openalex',0):,}")
    box(6.0, 10.6, 3.7, 1.1, f"Duplicates removed\n(n = {c['records_identified_raw'] - c['records_after_dedup']:,})", "#f6f6f4")
    arrow(3.0, 10.6, 3.0, 9.6)
    box(0.3, 8.5, 5.4, 1.1, f"Unique records screened by rules\n(title/abstract keyword scoring; n = {c['records_after_dedup']:,})")
    box(6.0, 8.5, 3.7, 1.1, f"Excluded by rules\n(no residential price outcome or no\nhedonic/regression signal; n = {c['excluded_by_rules']:,})", "#f6f6f4"); arrow(5.7, 9.05, 6.0, 9.05)
    arrow(3.0, 8.5, 3.0, 7.5)
    box(0.3, 6.4, 5.4, 1.1, f"Stage-1 candidates screened by independent agents\n(title/abstract; n = {c['screened_by_llm']:,})")
    ex = c["llm_exclusion_reasons"]; top = "; ".join(f"{k.split('_')[0]} {v}" for k, v in sorted(ex.items(), key=lambda kv: -kv[1])[:5])
    box(6.0, 6.0, 3.7, 1.5, f"Excluded with reason (n = {c['llm_excluded']:,})\n{top}\nunclear, flagged for review: {c['llm_unclear']:,}", "#f6f6f4"); arrow(5.7, 6.95, 6.0, 6.95)
    arrow(3.0, 6.4, 3.0, 5.4)
    box(0.3, 4.3, 5.4, 1.1, f"Included at stage 1: hedonic coefficient studies n = {c['included_hedonic_stage1']:,}\nML/AVM predictive studies n = {c['included_ml_predictive_stage1']:,} (AI synthesis)")
    arrow(3.0, 4.3, 3.0, 3.3)
    box(0.3, 2.2, 5.4, 1.1, f"Open-access full texts sought (hedonic + unclear): n = {c['fulltext_sought_hedonic_or_unclear']:,}\nretrieved n = {c['fulltext_retrieved']:,}")
    box(6.0, 2.2, 3.7, 1.1, f"Not retrievable legally\n(paywalled / no OA copy / download failed;\nn = {c['fulltext_not_retrieved']:,})", "#f6f6f4"); arrow(5.7, 2.75, 6.0, 2.75)
    arrow(3.0, 2.2, 3.0, 1.2)
    box(0.3, 0.1, 5.4, 1.1, f"Full texts assessed and extracted: n = {c['fulltext_assessed_extracted']:,}\nstudies eligible for coefficient synthesis: n = {c['fulltext_eligible']:,}", "#e3f1ea")
    box(6.0, 0.1, 3.7, 1.1, f"Excluded at full text\n(n = {c['fulltext_ineligible']:,})", "#f6f6f4"); arrow(5.7, 0.65, 6.0, 0.65)
    save(fig, "fig01_prisma_flow")


def fig02_trends():
    rec = pd.read_csv(META / "records_dedup.csv", low_memory=False, usecols=["record_id", "year", "stage1_candidate"]); llm = pd.read_csv(META / "screening_stage1_llm.csv")
    m = rec.merge(llm[["record_id", "reason_code", "ai_ml_flag"]], on="record_id", how="left")
    yrs = np.arange(1970, 2027)
    inc = m[m.reason_code == "ELIGIBLE_HEDONIC"].groupby("year").size().reindex(yrs, fill_value=0); ml = m[m.reason_code == "ELIGIBLE_ML_PREDICTIVE_ONLY"].groupby("year").size().reindex(yrs, fill_value=0)
    fig, ax = plt.subplots(1, 2, figsize=(8.2, 3.0))
    ax[0].bar(yrs, inc.values, color=PAL[0], width=0.8, label="hedonic coefficient studies (stage-1 include)"); ax[0].bar(yrs, ml.values, bottom=inc.values, color=PAL[1], width=0.8, label="ML / AVM predictive studies")
    ax[0].set_xlabel("publication year"); ax[0].set_ylabel("records included at stage 1"); ax[0].set_title("(a) Annual publication trend of included records"); ax[0].legend(fontsize=7); ax[0].set_xlim(1969, 2027)
    share = (ml / (inc + ml)).rolling(3, center=True, min_periods=1).mean()
    ax[1].plot(yrs, 100 * share.values, color=PAL[1]); ax[1].set_xlabel("publication year"); ax[1].set_ylabel("share of ML/AVM studies (%)"); ax[1].set_title("(b) Rise of machine-learning valuation studies"); ax[1].set_xlim(1995, 2027)
    save(fig, "fig02_publication_trends")


def fig03_map():
    S = pd.read_csv(FINAL / "studies.csv"); S = S[S.eligible == True]
    llm = pd.read_csv(META / "screening_stage1_llm.csv"); cc = yaml.safe_load((ROOT / "scripts" / "harmonization" / "countries.yaml").read_text())
    g1 = llm[llm.reason_code == "ELIGIBLE_HEDONIC"].country_guess.fillna("unknown").map(lambda c: str(c).split(",")[0].strip()).value_counts()
    fig, ax = plt.subplots(figsize=(8.4, 4.4)); ax.set_xlim(-170, 180); ax.set_ylim(-58, 80); ax.set_xlabel("longitude"); ax.set_ylabel("latitude"); ax.set_aspect("equal")
    for name, (cont, inc, iso, lat, lon) in cc.items():
        n1 = int(g1.get(name, 0)); n2 = int((S.country == name).sum())
        if n1 > 0: ax.scatter(lon, lat, s=6 + 2.0 * n1, color=PAL[0], alpha=0.35, lw=0)
        if n2 > 0: ax.scatter(lon, lat, s=6 + 14 * n2, color=PAL[1], alpha=0.85, lw=0.4, edgecolor=INK)
        if n1 >= 20 or n2 >= 3: ax.text(lon + 3, lat + 2, f"{iso} {n1}/{n2}", fontsize=6.5, color=INK2)
    ax.scatter([], [], color=PAL[0], alpha=0.35, s=60, label="stage-1 included hedonic records (abstract-based country)"); ax.scatter([], [], color=PAL[1], alpha=0.85, s=60, label="studies extracted for coefficient synthesis")
    ax.legend(loc="lower left", fontsize=7); ax.set_title("Geographic coverage (bubble area proportional to counts; label = ISO3 stage-1/extracted)")
    save(fig, "fig03_world_coverage")


def fig04_heatmap():
    C = pd.read_csv(FINAL / "coefficients_raw.csv"); S = pd.read_csv(FINAL / "studies.csv")
    C = C[C.taxonomy_id != "other"].merge(S[["record_id", "continent"]], on="record_id", how="left"); C["continent"] = C.continent.fillna("unknown")
    st = C.groupby(["taxonomy_id", "continent"]).record_id.nunique().unstack(fill_value=0)
    st = st.loc[st.sum(1).sort_values(ascending=False).index[:30]]
    fig, ax = plt.subplots(figsize=(7.2, 8.0)); im = ax.imshow(st.values, cmap="Blues", aspect="auto")
    ax.set_xticks(range(len(st.columns))); ax.set_xticklabels(st.columns, rotation=30, ha="right"); ax.set_yticks(range(len(st))); ax.set_yticklabels([LAB.get(i, i) for i in st.index]); ax.grid(False)
    for i in range(st.shape[0]):
        for j in range(st.shape[1]):
            v = st.values[i, j]
            if v: ax.text(j, i, str(v), ha="center", va="center", fontsize=6.5, color="white" if v > st.values.max() * 0.6 else INK)
    fig.colorbar(im, ax=ax, label="number of extracted studies reporting the characteristic"); ax.set_title("Variable frequency by continent (extracted studies)")
    save(fig, "fig04_variable_frequency_heatmap")


def forest(fam, name, title):
    p = pd.read_csv(ECO / "pooled_effects.csv"); p = p[p.family == fam].sort_values("est_re")
    if len(p) == 0: return
    fig, ax = plt.subplots(figsize=(7.0, 0.32 * len(p) + 1.6)); y = np.arange(len(p))
    ax.errorbar(p.est_re, y, xerr=[p.est_re - p.ci_lo, p.ci_hi - p.est_re], fmt="o", color=PAL[0], ms=4, capsize=2, lw=1, label="REML random effects, 95% CI")
    ax.hlines(y, p.pi_lo, p.pi_hi, color=PAL[0], alpha=0.25, lw=5, label="95% prediction interval")
    if p.est_ml.notna().any(): ax.scatter(p.est_ml, y + 0.18, marker="s", s=14, color=PAL[1], label="three-level (study/model) estimate")
    ax.axvline(0, color=INK2, lw=0.8); ax.set_yticks(y); ax.set_yticklabels([f"{LAB.get(t, t)}  (k={k}, n={n})" for t, k, n in zip(p.taxonomy_id, p.k_studies, p.n_estimates)])
    ax.set_xlabel(FAMLAB.get(fam, fam)); ax.set_title(title); ax.legend(fontsize=7, loc="lower right")
    lim = np.nanpercentile(np.r_[p.pi_lo, p.pi_hi], [2, 98]); ax.set_xlim(lim[0] - 0.1 * abs(lim[0]) - 1, lim[1] + 0.1 * abs(lim[1]) + 1)
    save(fig, name)


def fig05_forests():
    forest("pct_premium", "fig05_forest_pct_premium", "Pooled percentage premiums of binary attributes"); forest("semi_elasticity", "fig06_forest_semi_elasticity", "Pooled semi-elasticities of continuous attributes"); forest("elasticity", "fig07_forest_elasticity", "Pooled elasticities (log-log models)")


def fig08_geo_heterogeneity():
    E = pd.read_csv(FINAL / "effects.csv"); top = E.groupby(["family", "taxonomy_id"]).study_id.nunique().sort_values(ascending=False)
    top = [k for k in top.index if top[k] >= 8][:6]
    if not top: return
    fig, axes = plt.subplots(2, 3, figsize=(9.6, 5.6)); axes = axes.ravel()
    for a_, (fam, tx) in zip(axes, top):
        s = E[(E.family == fam) & (E.taxonomy_id == tx)]; s = s[s.continent.notna()]
        order = s.groupby("continent").study_id.nunique().sort_values(ascending=False).index
        data = [s[s.continent == c].es.values for c in order]
        a_.boxplot(data, vert=False, labels=[f"{c} (k={s[s.continent == c].study_id.nunique()})" for c in order], widths=0.6, showfliers=False, medianprops=dict(color=PAL[1]))
        a_.axvline(0, color=INK2, lw=0.8); a_.set_title(f"{LAB.get(tx, tx)} — {FAMLAB.get(fam, fam).split(' (')[0]}", fontsize=8); a_.set_xlabel("harmonised effect size")
    for a_ in axes[len(top):]: a_.axis("off")
    save(fig, "fig08_geographic_heterogeneity")


def fig09_temporal():
    E = pd.read_csv(FINAL / "effects.csv"); E = E[E.data_mid_year.notna()]
    top = E.groupby(["family", "taxonomy_id"]).study_id.nunique().sort_values(ascending=False); top = [k for k in top.index if top[k] >= 8][:6]
    if not top: return
    fig, axes = plt.subplots(2, 3, figsize=(9.6, 5.6)); axes = axes.ravel()
    for a_, (fam, tx) in zip(axes, top):
        s = E[(E.family == fam) & (E.taxonomy_id == tx)]
        a_.scatter(s.data_mid_year, s.es, s=np.clip(8 / np.maximum(s.es_se, 1e-6), 4, 60), color=PAL[0], alpha=0.5, lw=0)
        if len(s) > 5:
            w = (1 / s.es_se ** 2).values; X = np.column_stack([np.ones(len(s)), s.data_mid_year.values - 2000]); b = np.linalg.lstsq(X * np.sqrt(w)[:, None], s.es.values * np.sqrt(w), rcond=None)[0]
            xs = np.linspace(s.data_mid_year.min(), s.data_mid_year.max(), 50); a_.plot(xs, b[0] + b[1] * (xs - 2000), color=PAL[1], label=f"WLS trend {b[1]:+.3f}/yr")
            a_.legend(fontsize=7)
        a_.axhline(0, color=INK2, lw=0.8); a_.set_title(f"{LAB.get(tx, tx)}", fontsize=8); a_.set_xlabel("midpoint of the transaction sample period"); a_.set_ylabel("harmonised effect")
    for a_ in axes[len(top):]: a_.axis("off")
    save(fig, "fig09_temporal_evolution")


def fig10_funnel():
    E = pd.read_csv(FINAL / "effects.csv"); top = E.groupby(["family", "taxonomy_id"]).size().sort_values(ascending=False); top = [k for k in top.index if top[k] >= 15][:6]
    if not top: return
    fig, axes = plt.subplots(2, 3, figsize=(9.6, 5.6)); axes = axes.ravel()
    for a_, (fam, tx) in zip(axes, top):
        s = E[(E.family == fam) & (E.taxonomy_id == tx)]; w = 1 / s.es_se ** 2; mu = (w * s.es).sum() / w.sum()
        a_.scatter(s.es, s.es_se, s=10, color=PAL[0], alpha=0.6, lw=0); se = np.linspace(0, s.es_se.max(), 50)
        a_.plot(mu - 1.96 * se, se, color=INK2, lw=0.7, ls="--"); a_.plot(mu + 1.96 * se, se, color=INK2, lw=0.7, ls="--"); a_.axvline(mu, color=PAL[1], lw=0.9)
        a_.invert_yaxis(); a_.set_title(f"{LAB.get(tx, tx)} (n={len(s)})", fontsize=8); a_.set_xlabel("harmonised effect"); a_.set_ylabel("standard error")
    for a_ in axes[len(top):]: a_.axis("off")
    save(fig, "fig10_funnel_plots")


def fig11_moderators():
    f = ECO / "meta_regression.csv"
    if not f.exists(): return
    m = pd.read_csv(f); m = m[m.term != "intrcpt"]
    if len(m) == 0: return
    keys = m.groupby(["family", "taxonomy_id"]).size().index[:6]
    fig, axes = plt.subplots(1, len(keys), figsize=(2.6 * len(keys) + 1, 3.6), sharey=True); axes = np.atleast_1d(axes)
    for a_, (fam, tx) in zip(axes, keys):
        s = m[(m.family == fam) & (m.taxonomy_id == tx)]; y = np.arange(len(s))
        a_.errorbar(s.est, y, xerr=1.96 * s.se_rve.fillna(s.se), fmt="o", color=PAL[0], ms=4, capsize=2, lw=1); a_.axvline(0, color=INK2, lw=0.8)
        a_.set_yticks(y); a_.set_yticklabels(s.term); a_.set_title(f"{LAB.get(tx, tx)}", fontsize=8); a_.set_xlabel("moderator effect (95% RVE CI)")
    save(fig, "fig11_moderator_effects")


def fig12_evidence_density():
    C = pd.read_csv(FINAL / "coefficients_raw.csv"); S = pd.read_csv(FINAL / "studies.csv")
    C = C[C.taxonomy_id != "other"].merge(S[["record_id", "sample_start", "sample_end", "year"]], on="record_id", how="left")
    for c in ("sample_start", "sample_end", "year"):
        C[c] = pd.to_numeric(C[c], errors="coerce"); C.loc[(C[c] < 1900) | (C[c] > 2027), c] = np.nan
    C["mid"] = C[["sample_start", "sample_end"]].mean(axis=1).fillna(C.year - 3); C = C[C.mid.notna()]; C["decade"] = (C.mid // 10 * 10).astype("Int64")
    st = C.groupby(["taxonomy_id", "decade"]).record_id.nunique().unstack(fill_value=0); st = st.loc[st.sum(1).sort_values(ascending=False).index[:30]]
    fig, ax = plt.subplots(figsize=(6.4, 8.0)); im = ax.imshow(st.values, cmap="Oranges", aspect="auto"); ax.grid(False)
    ax.set_xticks(range(len(st.columns))); ax.set_xticklabels([f"{int(d)}s" for d in st.columns]); ax.set_yticks(range(len(st))); ax.set_yticklabels([LAB.get(i, i) for i in st.index])
    for i in range(st.shape[0]):
        for j in range(st.shape[1]):
            if st.values[i, j]: ax.text(j, i, str(st.values[i, j]), ha="center", va="center", fontsize=6.5, color=INK)
    fig.colorbar(im, ax=ax, label="studies"); ax.set_title("Evidence density: characteristic × decade of the transaction data")
    save(fig, "fig12_evidence_density")


def fig13_ai_evolution():
    llm = pd.read_csv(META / "screening_stage1_llm.csv"); rec = pd.read_csv(META / "records_dedup.csv", low_memory=False, usecols=["record_id", "year", "title", "abstract"])
    m = llm.merge(rec, on="record_id"); m = m[m.decision == "include"]; t = (m.title.fillna("") + " " + m.abstract.fillna("")).str.lower()
    fams = {"random forest / boosting": r"random forest|gradient boost|xgboost|lightgbm|catboost|boosting", "neural networks / deep learning": r"neural network|deep learning|cnn|lstm|transformer", "computer vision / images": r"street view|satellite|image|computer vision|photo", "text / NLP": r"text|nlp|natural language|sentiment|listing description|word", "spatial econometrics": r"spatial lag|spatial error|geographically weighted|gwr|spatial autoregress|spatial durbin|spatial econometric"}
    yrs = np.arange(2000, 2027); fig, ax = plt.subplots(figsize=(7.2, 3.2))
    for i, (lab, pat) in enumerate(fams.items()):
        cnt = m[t.str.contains(pat)].groupby("year").size().reindex(yrs, fill_value=0).rolling(3, center=True, min_periods=1).mean(); ax.plot(yrs, cnt.values, color=PAL[i], label=lab)
    ax.set_xlabel("publication year"); ax.set_ylabel("included records mentioning the method (3-yr moving average)"); ax.set_title("Methodological evolution in the included literature (title/abstract mentions)"); ax.legend(fontsize=7)
    save(fig, "fig13_ai_methodology_evolution")


def fig14_summary():
    p = pd.read_csv(ECO / "pooled_effects.csv"); p = p[p.family.isin(["pct_premium", "semi_elasticity"])].copy(); p["label"] = p.taxonomy_id.map(lambda t: LAB.get(t, t)) + p.family.map({"pct_premium": " (premium %)", "semi_elasticity": " (% per unit)"})
    p = p.sort_values("est_re"); fig, ax = plt.subplots(figsize=(7.4, 0.26 * len(p) + 1.5)); y = np.arange(len(p))
    cols = [PAL[GRP.get(t, "A")[0] in "A" and 0 or (1 if GRP.get(t, "")[0] in "BC" else 2 if GRP.get(t, "")[0] == "D" else 3 if GRP.get(t, "")[0] == "E" else 4)] for t in p.taxonomy_id]
    ax.barh(y, p.est_re, color=cols, height=0.6); ax.errorbar(p.est_re, y, xerr=[p.est_re - p.ci_lo, p.ci_hi - p.est_re], fmt="none", ecolor=INK2, capsize=2, lw=0.8)
    ax.set_yticks(y); ax.set_yticklabels(p.label); ax.axvline(0, color=INK2, lw=0.8); ax.set_xlabel("pooled estimate (REML) with 95% CI"); ax.set_title("Summary of the pooled evidence (colour = taxonomy group)")
    save(fig, "fig14_summary_evidence")


ALL = [fig01_prisma, fig02_trends, fig03_map, fig04_heatmap, fig05_forests, fig08_geo_heterogeneity, fig09_temporal, fig10_funnel, fig11_moderators, fig12_evidence_density, fig13_ai_evolution, fig14_summary]
if __name__ == "__main__":
    import sys, traceback
    bad = []
    for f in ALL:
        try: f()
        except Exception as e: bad.append(f.__name__); print("FAILED", f.__name__, repr(e)); traceback.print_exc()
    sys.exit(1 if bad else 0)