#!/usr/bin/env python """Stage 25 — whole-dataset (2003-2022) out-of-sample verification of the v3 zoo (stage 21). "Whole dataset" = the 20 pooled LEAVE-ONE-YEAR-OUT held-out years, so the leakage rule holds (every prediction was made by a model that never saw its year). For each family x lead block we report the full binary contingency inventory the hybrid-design note asks for (doc §4/§8): probabilistic (threshold-free): Brier, BSS vs climatology, ROC-AUC, PR-AUC (rare-event aware). binary (at the leakage-safe per-region x block decision threshold stored per fold): POD, POFD, FAR, CSI, ETS, HSS, Peirce, EDI, SEDI, frequency bias. Binary counts are accumulated across folds using EACH fold's own decision threshold (fit on the other 19 years), then scored once on the pooled contingency table — no threshold ever sees its own test year. Climatology reference = per-block base rate of the pooled held-out labels. Outputs: results/tables/v3zoo_scores_by_model_block.csv family x block: all scores results/tables/v3zoo_best_model_by_block.csv best family per block (by SEDI and by BSS) results/figures/v3zoo_skill_by_lead.png BSS + SEDI vs lead for the leading families Consumes outputs/loyo_pred_v3zoo/pred_v3zoo_*.npz (+ optional stack from stage 23 if present). """ 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 from sklearn.metrics import average_precision_score 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") def load_folds(indir, pattern): files = sorted(glob.glob(os.path.join(indir, pattern))) if not files: raise SystemExit(f"no {pattern} in {indir} — run stage 21 first") return files def family_cols(d): return sorted(c for c in d.files if c.startswith("p_")) def main(): ap = argparse.ArgumentParser() ap.add_argument("--indir", default=None) ap.add_argument("--pattern", default="pred_v3zoo_*.npz") ap.add_argument("--out-prefix", default="v3zoo", help="output-table stem, e.g. v4ff -> v4ff_scores_by_model_block.csv") args = ap.parse_args() pref = args.out_prefix cfg = load_config(); blocks = lead_blocks(cfg); block_names = list(blocks) block_ids = {b: k + 1 for k, b in enumerate(block_names)} indir = args.indir or os.path.join(cfg["paths"]["outputs"], "loyo_pred_v3zoo") files = load_folds(indir, args.pattern) # pool per (family, block): y, p, and a leakage-safe per-fold per-region flag d0 = np.load(files[0], allow_pickle=True) fams = [c[2:] for c in family_cols(d0)] print(f"[25] {len(files)} folds | {len(fams)} families: {fams}") pool = {(f, b): {"y": [], "p": [], "a": 0, "b": 0, "c": 0, "d": 0} for f in fams for b in block_names} for path in files: d = np.load(path, allow_pickle=True) bid = d["block_id"].astype(int); reg = d["region"].astype(int); y = d["y"].astype(float) dec = d["decthr"] if "decthr" in d.files else None for b in block_names: mb = bid == block_ids[b] for f in fams: p = d["p_" + f][mb].astype(float); yb = y[mb]; rb = reg[mb] ok = np.isfinite(yb) & np.isfinite(p) pool[(f, b)]["y"].append(yb[ok]); pool[(f, b)]["p"].append(p[ok]) if dec is not None: sel = (dec["tag"] == f) & (dec["block"] == b) thr_by_reg = {int(r): float(t) for r, t in zip(dec["region"][sel], dec["thr"][sel])} if thr_by_reg: gthr = float(np.median(list(thr_by_reg.values()))) tvec = np.array([thr_by_reg.get(int(r), gthr) for r in rb])[ok] a, bb, c, dd = vf.contingency(yb[ok], p[ok], tvec) pool[(f, b)]["a"] += a; pool[(f, b)]["b"] += bb pool[(f, b)]["c"] += c; pool[(f, b)]["d"] += dd rows = [] for f in fams: for b in block_names: P = pool[(f, b)] y = np.concatenate(P["y"]) if P["y"] else np.array([]) p = np.concatenate(P["p"]) if P["p"] else np.array([]) if y.size == 0: continue base = float(np.mean(y)) rec = dict(family=f, block=b, base_rate=round(base, 4), n=int(y.size), brier=round(vf.brier(y, p), 5), bss=round(vf.brier_skill_score(y, p, np.full(y.size, base)), 4), auc=round(vf.auc(y, p), 4), pr_auc=round(float(average_precision_score(y, p)) if y.sum() > 0 else np.nan, 4)) if (P["a"] + P["b"] + P["c"] + P["d"]) > 0: bs = vf.binary_scores_from_counts(P["a"], P["b"], P["c"], P["d"]) for k in ("POD", "POFD", "FAR", "CSI", "ETS", "HSS", "Peirce", "EDI", "SEDI", "freq_bias"): rec[k] = round(float(bs[k]), 4) rows.append(rec) df = pd.DataFrame(rows) os.makedirs(TABLES, exist_ok=True) df.to_csv(os.path.join(TABLES, f"{pref}_scores_by_model_block.csv"), index=False) print(f"[25] wrote {pref}_scores_by_model_block.csv") best = [] for b in block_names: sub = df[df.block == b] if sub.empty: continue best.append(dict(block=b, best_by_SEDI=sub.loc[sub.SEDI.idxmax(), "family"] if "SEDI" in sub else "", SEDI=float(sub.SEDI.max()) if "SEDI" in sub else np.nan, best_by_BSS=sub.loc[sub.bss.idxmax(), "family"], BSS=float(sub.bss.max()))) pd.DataFrame(best).to_csv(os.path.join(TABLES, f"{pref}_best_model_by_block.csv"), index=False) print(f"[25] wrote {pref}_best_model_by_block.csv") for r in best: print(f" {r['block']}: best BSS {r['best_by_BSS']}={r['BSS']:.3f} | " f"best SEDI {r['best_by_SEDI']}={r['SEDI']:.3f}") _plot(df, block_names, pref) def _plot(df, block_names, pref="v3zoo"): 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)) # show the leading families by mean BSS lead = (df.groupby("family").bss.mean().sort_values(ascending=False).head(5).index.tolist()) for f in lead: s = df[df.family == f].set_index("block").reindex(block_names) ax[0].plot(x, s.bss.values, marker="o", label=f) if "SEDI" in s: ax[1].plot(x, s.SEDI.values, marker="o", label=f) for a, t in zip(ax, ("BSS vs climatology", "SEDI")): a.set_xticks(x); a.set_xticklabels(block_names); a.set_title(t) a.axhline(0, color="k", lw=0.6); a.grid(alpha=0.3); a.legend(fontsize=8) fig.suptitle(f"{pref} — out-of-sample skill by lead block (2003-2022 pooled)") os.makedirs(FIGS, exist_ok=True) fig.tight_layout(); fig.savefig(os.path.join(FIGS, f"{pref}_skill_by_lead.png"), dpi=130) print(f"[25] wrote results/figures/{pref}_skill_by_lead.png") if __name__ == "__main__": main()