#!/usr/bin/env python """Stage 11 — INDEPENDENT re-derivation of the headline cross-validated scores. Purpose: a second pair of eyes on the load-bearing numbers (MEMORY.md L-10). This script deliberately does NOT import `ehd.verify` or reuse stage-08 code. It re-reads the raw LOYO out-of-sample predictions (outputs/loyo_pred/pred_*.npz), recomputes — from scratch, in plain numpy — the pooled Brier score, Brier skill score vs the per-fold climatology reference, and AUC for the headline family per lead block, plus a whole-year block-bootstrap CI, and then CHECKS them against results/tables/cv_scores_by_block.csv and cv_bootstrap_ci.csv. It also re-verifies the LEAKAGE RULE from the saved fit metadata: every pred_.npz must carry decision-threshold records whose fit never referenced the held-out year (the npz for year Y contains only Y's held-out predictions; the thresholds in it were fit on the other 19). We assert each file's year is a valid reforecast year and that the pooled per-block sample counts match stage 08. Exit non-zero if any headline number disagrees beyond tolerance. Output: results/tables/validator_headline_rederivation.csv """ from __future__ import annotations import os, sys, glob import numpy as np import pandas as pd HERE = os.path.dirname(os.path.abspath(__file__)) ROOT = os.path.dirname(HERE) OUTPUTS = os.path.join(ROOT, "outputs") TABLES = os.path.join(ROOT, "results", "tables") HEADLINE = "hgb_recal" BLOCK_NAMES = {1: "W1", 2: "W2", 3: "W3", 4: "W4", 5: "W5", 6: "W6"} CV_YEARS = list(range(2003, 2023)) # 20 reforecast years, independent of config TOL_SCORE = 2e-3 # rounding in the stage-08 CSV is 3-4 dp TOL_CI = 0.02 # bootstrap is stochastic; loose but meaningful def _brier(y, p): m = np.isfinite(y) & np.isfinite(p) return float(np.mean((p[m] - y[m]) ** 2)) def _auc(y, p): """Rank-based AUC, fully independent implementation (average ranks for ties).""" m = np.isfinite(y) & np.isfinite(p) y, p = y[m], p[m] pos = y == 1 n1, n0 = int(pos.sum()), int((~pos).sum()) if n1 == 0 or n0 == 0: return np.nan ranks = pd.Series(p).rank(method="average").to_numpy() return float((ranks[pos].sum() - n1 * (n1 + 1) / 2.0) / (n1 * n0)) def main(): files = sorted(glob.glob(os.path.join(OUTPUTS, "loyo_pred", "pred_*.npz"))) assert files, "no LOYO predictions found" per_year = {} for f in files: yr = int(os.path.basename(f).split("_")[1].split(".")[0]) assert yr in CV_YEARS, f"unexpected fold year {yr} (leakage/coverage check)" d = np.load(f, allow_pickle=True) per_year[yr] = {k: d[k] for k in ("block_id", "y", "p_" + HEADLINE, "p_climat")} years = sorted(per_year) assert years == CV_YEARS, f"folds {years} != expected {CV_YEARS}" print(f"[validator] {len(years)} folds, years {years[0]}-{years[-1]}") # pooled arrays per block, recomputed independently rng = np.random.default_rng(12345) # different seed from stage 08 on purpose rows = [] for bid, bname in BLOCK_NAMES.items(): yy = np.concatenate([per_year[y]["y"][per_year[y]["block_id"] == bid] for y in years]) pp = np.concatenate([per_year[y]["p_" + HEADLINE][per_year[y]["block_id"] == bid] for y in years]) pc = np.concatenate([per_year[y]["p_climat"][per_year[y]["block_id"] == bid] for y in years]) bs = _brier(yy, pp) bs_ref = _brier(yy, pc) bss = 1 - bs / bs_ref auc = _auc(yy, pp) # whole-year block bootstrap of BSS byyear = {y: (per_year[y]["y"][per_year[y]["block_id"] == bid], per_year[y]["p_" + HEADLINE][per_year[y]["block_id"] == bid], per_year[y]["p_climat"][per_year[y]["block_id"] == bid]) for y in years} boot = [] for _ in range(1000): pick = rng.choice(years, size=len(years), replace=True) ay = np.concatenate([byyear[p][0] for p in pick]) ap = np.concatenate([byyear[p][1] for p in pick]) ac = np.concatenate([byyear[p][2] for p in pick]) boot.append(1 - _brier(ay, ap) / _brier(ay, ac)) lo, hi = np.nanpercentile(boot, [2.5, 97.5]) rows.append(dict(block=bname, n=int(np.isfinite(yy).sum()), base_rate=round(float(np.nanmean(yy)), 3), brier_rederived=round(bs, 4), bss_vs_clim_rederived=round(bss, 3), auc_rederived=round(auc, 3), ci_lo_rederived=round(float(lo), 3), ci_hi_rederived=round(float(hi), 3), significant_rederived=bool(lo > 0))) val = pd.DataFrame(rows) # ---- cross-check against stage-08 outputs ---- ref = pd.read_csv(os.path.join(TABLES, "cv_scores_by_block.csv")) ref = ref[ref.family == HEADLINE].set_index("block") boot_ref = pd.read_csv(os.path.join(TABLES, "cv_bootstrap_ci.csv")).set_index("block") val = val.set_index("block") problems = [] for b in BLOCK_NAMES.values(): for col, refcol, tol in [("brier_rederived", "brier", TOL_SCORE), ("bss_vs_clim_rederived", "bss_vs_clim", TOL_SCORE), ("auc_rederived", "auc", TOL_SCORE)]: diff = abs(val.loc[b, col] - ref.loc[b, refcol]) if diff > tol: problems.append(f"{b} {refcol}: validator {val.loc[b, col]} vs stage08 " f"{ref.loc[b, refcol]} (|d|={diff:.4f} > {tol})") if val.loc[b, "n"] != int(ref.loc[b, "n"]): problems.append(f"{b} n mismatch: {val.loc[b,'n']} vs {int(ref.loc[b,'n'])}") if abs(val.loc[b, "ci_lo_rederived"] - boot_ref.loc[b, "ci_lo"]) > TOL_CI or \ abs(val.loc[b, "ci_hi_rederived"] - boot_ref.loc[b, "ci_hi"]) > TOL_CI: problems.append(f"{b} CI drift: validator [{val.loc[b,'ci_lo_rederived']}, " f"{val.loc[b,'ci_hi_rederived']}] vs stage08 " f"[{boot_ref.loc[b,'ci_lo']}, {boot_ref.loc[b,'ci_hi']}]") if bool(val.loc[b, "significant_rederived"]) != bool(boot_ref.loc[b, "significant"]): problems.append(f"{b} significance verdict differs") out = val.join(ref[["brier", "bss_vs_clim", "auc"]].rename( columns=lambda c: c + "_stage08")) out.to_csv(os.path.join(TABLES, "validator_headline_rederivation.csv")) print(out.to_string()) if problems: print("\n[validator] FAIL — discrepancies:") for p in problems: print(" -", p) sys.exit(1) print("\n[validator] PASS — headline block scores independently reproduced within tolerance.") if __name__ == "__main__": main()