#!/usr/bin/env python """Stage 20g — verify the calibration-recovery models (stage 20f): base / sel / stack, out-of-sample. Pools LOYO predictions from outputs/loyo_pred_v3cal/ and scores each model per block (BSS vs per-fold climatology, AUC, full binary inventory), with a whole-year block bootstrap on the BSS gain of stack and sel over base — the test of whether proper OOF calibration + stacking converts the discrimination gain (stage 20d/20e) into calibrated probability skill. Outputs: results/tables/v3cal_scores_by_block.csv base/sel/stack: base_rate, Brier, BSS_clim, AUC, SEDI results/tables/v3cal_bss_gain_bootstrap.csv per block: BSS gain vs base for sel and stack + 95% CI results/tables/v3cal_binary_by_block.csv full 10-metric inventory per model results/figures/v3cal_skill_by_lead.png BSS & AUC by lead: base/sel/stack """ from __future__ import annotations import os, sys, glob, argparse import numpy as np import pandas as pd sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..", "src")) from ehd.config import load_config, lead_blocks from ehd import verify as vf ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) TABLES = os.path.join(ROOT, "results", "tables") FIGS = os.path.join(ROOT, "results", "figures") MODELS = ["base", "sel", "stack"] def load(cfg, pat): paths = sorted(glob.glob(os.path.join(cfg["paths"]["outputs"], "loyo_pred_v3cal", pat))) if not paths: raise SystemExit(f"no outputs/loyo_pred_v3cal/{pat} — run stage 20f first") frames = [] for p in paths: yr = int(os.path.basename(p).split("_")[-1][:-4]) d = np.load(p, allow_pickle=True) df = pd.DataFrame({k: d[k] for k in d.files}); df["year"] = yr frames.append(df) return pd.concat(frames, ignore_index=True) def main(): ap = argparse.ArgumentParser() ap.add_argument("--pat", default="pred_v3cal_*.npz", help="glob (use pred_v3cal_nocap_*.npz for no-cap)") ap.add_argument("--suffix", default="", help="output filename suffix, e.g. _nocap") args = ap.parse_args() cfg = load_config(); blocks = lead_blocks(cfg) bid = {k + 1: b for k, b in enumerate(blocks)} years = list(cfg["cross_validation"]["years"]) df = load(cfg, args.pat) rows, irows, brows = [], [], [] rng = np.random.default_rng(cfg["seed"]); B = 1000 for b_id, bname in bid.items(): sub = df[df.block_id == b_id] y = sub.y.values; pc = np.full(len(sub), float(np.nanmean(y))) for tag in MODELS: p = sub[f"p_{tag}"].values; fl = sub[f"flag_{tag}"].values.astype(bool) m = np.isfinite(y) & np.isfinite(p); yv, pv, fv = y[m], p[m], fl[m] a = int(np.sum(fv & (yv == 1))); bb = int(np.sum(fv & (yv == 0))) c = int(np.sum(~fv & (yv == 1))); d = int(np.sum(~fv & (yv == 0))) bs = vf.binary_scores_from_counts(a, bb, c, d) rows.append(dict(block=bname, model=tag, base_rate=round(float(np.nanmean(yv)), 3), brier=round(float(np.nanmean((pv - yv) ** 2)), 4), bss_vs_clim=round(vf.brier_skill_score(y, p, pc), 3), auc=round(vf.auc(yv, pv), 3), SEDI=round(bs["SEDI"], 3))) irows.append(dict(block=bname, model=tag, POD=round(bs["POD"], 3), POFD=round(bs["POFD"], 3), FAR=round(bs["FAR"], 3), CSI=round(bs["CSI"], 3), ETS=round(bs["ETS"], 3), HSS=round(bs["HSS"], 3), Peirce=round(bs["Peirce"], 3), EDI=round(bs["EDI"], 3), SEDI=round(bs["SEDI"], 3), freq_bias=round(bs["freq_bias"], 2))) per_year = {yv: sub[sub.year == yv] for yv in years} def bss_of(col, pick): yy = np.concatenate([per_year[q].y.values for q in pick]) pp = np.concatenate([per_year[q][col].values for q in pick]) return vf.brier_skill_score(yy, pp, np.full(len(yy), float(np.nanmean(yy)))) gains = {"sel": [], "stack": []} for _ in range(B): pick = rng.choice(years, size=len(years), replace=True) gb = bss_of("p_base", pick) gains["sel"].append(bss_of("p_sel", pick) - gb) gains["stack"].append(bss_of("p_stack", pick) - gb) base_bss = {r["model"]: r["bss_vs_clim"] for r in rows if r["block"] == bname} for tag in ("sel", "stack"): lo, hi = np.nanpercentile(gains[tag], [2.5, 97.5]) brows.append(dict(block=bname, model=tag, bss_base=base_bss["base"], bss_model=base_bss[tag], bss_gain=round(base_bss[tag] - base_bss["base"], 3), gain_ci_lo=round(float(lo), 3), gain_ci_hi=round(float(hi), 3), significant=bool(lo > 0))) pd.DataFrame(rows).to_csv(os.path.join(TABLES, f"v3cal_scores_by_block{args.suffix}.csv"), index=False) pd.DataFrame(brows).to_csv(os.path.join(TABLES, f"v3cal_bss_gain_bootstrap{args.suffix}.csv"), index=False) pd.DataFrame(irows).to_csv(os.path.join(TABLES, f"v3cal_binary_by_block{args.suffix}.csv"), index=False) print(f"[20g{args.suffix}] scores:\n", pd.DataFrame(rows).to_string(index=False)) print(f"[20g{args.suffix}] BSS gain vs base:\n", pd.DataFrame(brows).to_string(index=False)) _fig(pd.DataFrame(rows), list(blocks), args.suffix) def _fig(df, block_names, suffix): import matplotlib; matplotlib.use("Agg"); import matplotlib.pyplot as plt x = np.arange(len(block_names)) fig, ax = plt.subplots(1, 2, figsize=(11, 4.3)) for tag, col in (("base", "0.5"), ("sel", "tab:orange"), ("stack", "tab:red")): s = df[df.model == tag].set_index("block").reindex(block_names) ax[0].plot(x, s.bss_vs_clim, marker="o", color=col, label=tag) ax[1].plot(x, s.auc, marker="o", color=col, label=tag) for a, t in zip(ax, ["BSS vs climatology", "AUC"]): a.set_xticks(x); a.set_xticklabels(block_names) a.axhline(0 if t.startswith("BSS") else 0.5, color="k", lw=.6) a.set_title(f"Stage-2 calibration recovery — {t}"); a.grid(alpha=.3); a.legend() os.makedirs(FIGS, exist_ok=True); fig.tight_layout() fig.savefig(os.path.join(FIGS, f"v3cal_skill_by_lead{suffix}.png"), dpi=130) print(f"[20g] wrote results/figures/v3cal_skill_by_lead{suffix}.png") if __name__ == "__main__": main()