"""Stage 19 — prune the weakest members with NESTED (leakage-safe) selection -> registry v2. Stage 2 : pruning is a fitted decision, so it must live INSIDE the outer LOYO. For each held-out year y the prune/keep decision is computed from the OTHER 19 years only (their stored out-of-sample votes), then applied to y. The member set may differ by lead week (a member weak at W2 may be kept at W5) — allowed and expected. Criterion (STAGE2_PLAN §0 default, from the user's stated priority = high POD, low POFD/FAR): Peirce skill score PSS = POD - POFD per lead week W2-W6, pooled over regions, subject to a FAR ceiling; drop the bottom tier ONLY IF the vote-ensemble PSS on the training years does not fall (ensemble-effect check — a weak-but-decorrelated member that helps the ensemble is kept). Consensus registry v2: a member is marked `pruned` if it is dropped in >= 50% of the outer folds (aggregated over W2-W6 fold-week decisions). We also report the OPTIMISM GAP: the score of a set chosen on all 20 years and scored on all 20 years, minus the nested honest score. Each member's vote for year t is reconstructed leakage-safely from the stored artefacts: vote = (stored P(yEHDE) for year t) >= (stored per-region x block threshold fit on years != t). Inputs : outputs/loyo_pred_v2/pred_v2_*.npz, outputs/loyo_pred_v2_cnn/pred_v2cnn_*.npz, outputs/member_registry/registry_v1.json Outputs: outputs/member_registry/registry_v2.json + .csv results/tables/v2_pruning_ranking.csv (per fold x week x member: train PSS/FAR, rank) results/tables/v2_pruning_consensus.csv (per week x member: times dropped / 20, status) results/tables/v2_pruning_ensemble_effect.csv (full vs pruned ensemble, nested + optimism gap) results/figures/v2_pruning_ranking.png """ from __future__ import annotations import os, sys, json, glob, logging 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 from ehd import verify as vf log = logging.getLogger("stage19") logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s") FAR_CEILING = 0.62 # members above this FAR at a week are prune candidates regardless of PSS BOTTOM_FRAC = 0.34 # bottom third by train PSS are prune candidates EPS = 1e-3 # ensemble PSS may not fall by more than this to accept a prune MIN_KEEP = 10 # never prune below this many members CONSENSUS = 0.50 # dropped in >= this fraction of folds -> pruned in v2 def _peirce_from_votefrac(y, frac): """Max Peirce (POD-POFD) over thresholds on a member-vote fraction, and the argmax threshold.""" grid = np.linspace(0.05, 0.95, 19) best, bt = -np.inf, 0.5 for t in grid: s = vf.binary_scores(y, frac, t)["Peirce"] if np.isfinite(s) and s > best: best, bt = s, t return best, bt def load_votes(cfg): """Return dict year -> dict(block -> dict(votes[n,M] uint8, y[n], region[n])) and member list.""" root = cfg["paths"]["outputs"] reg = json.load(open(os.path.join(root, "member_registry", "registry_v1.json"))) members = [m["member_id"] for m in reg["members"]] prob_col = {m["member_id"]: m["provenance"]["prob_column"] for m in reg["members"]} is_torch = {m["member_id"]: (m["family"] in ("cnn", "ann_torch")) for m in reg["members"]} blocks = list(cfg["lead_blocks"].keys()) bid_of = {b: i + 1 for i, b in enumerate(blocks)} # W1->1..W6->6 (matches stored block_id) main = {int(os.path.basename(p).split("_")[-1][:-4]): p for p in glob.glob(os.path.join(root, "loyo_pred_v2", "pred_v2_*.npz"))} cnn = {int(os.path.basename(p).split("_")[-1][:-4]): p for p in glob.glob(os.path.join(root, "loyo_pred_v2_cnn", "pred_v2cnn_*.npz"))} years = sorted(main) def thr_lookup(decthr): d = {} for row in decthr: d[(str(row["tag"]), str(row["block"]), int(row["region"]))] = float(row["thr"]) return d data = {} for y in years: dm = np.load(main[y], allow_pickle=True) dc = np.load(cnn[y], allow_pickle=True) thr = thr_lookup(dm["decthr"]); thr.update(thr_lookup(dc["decthr"])) bid = dm["block_id"].astype(int); reg_arr = dm["region"].astype(int) yv = dm["y"].astype(float) cbid = dc["block_id"].astype(int) assert np.array_equal(dm["cell"].astype(int), dc["cell"].astype(int)), f"cell mismatch {y}" data[y] = {} for b in blocks: bcode = bid_of[b] m = bid == bcode reg_b = reg_arr[m]; y_b = yv[m] V = np.zeros((m.sum(), len(members)), dtype=np.uint8) for j, mid in enumerate(members): src = dc if is_torch[mid] else dm p = src[prob_col[mid]][(cbid == bcode) if is_torch[mid] else m].astype(float) tvec = np.array([thr.get((mid, b, int(r)), 0.5) for r in reg_b]) V[:, j] = (p >= tvec).astype(np.uint8) data[y][b] = dict(votes=V, y=y_b, region=reg_b) return data, members, blocks def member_scores(votes_col, y): a = int(np.sum((votes_col == 1) & (y == 1))) b = int(np.sum((votes_col == 1) & (y == 0))) c = int(np.sum((votes_col == 0) & (y == 1))) d = int(np.sum((votes_col == 0) & (y == 0))) s = vf.binary_scores_from_counts(a, b, c, d) return s["Peirce"], s["FAR"] def prune_decision(Vtr, ytr, members): """Return (drop_set indices, full_pss, pruned_pss) on the TRAIN pool for one week.""" M = Vtr.shape[1] pss = np.array([member_scores(Vtr[:, j], ytr)[0] for j in range(M)]) far = np.array([member_scores(Vtr[:, j], ytr)[1] for j in range(M)]) order = np.argsort(pss) # weakest first n_bottom = max(1, int(np.floor(BOTTOM_FRAC * M))) cand = set(order[:n_bottom].tolist()) cand |= set(np.where(far > FAR_CEILING)[0].tolist()) # never prune below MIN_KEEP cand = set(sorted(cand, key=lambda j: pss[j])[:max(0, M - MIN_KEEP)]) keep_full = np.arange(M) keep_pruned = np.array([j for j in range(M) if j not in cand]) full_pss, _ = _peirce_from_votefrac(ytr, Vtr[:, keep_full].mean(1)) pruned_pss, _ = _peirce_from_votefrac(ytr, Vtr[:, keep_pruned].mean(1)) if len(keep_pruned) else (np.nan, .5) if not (pruned_pss >= full_pss - EPS): # ensemble-effect gate: reject the prune cand = set() return cand, pss, far, full_pss, pruned_pss def main(): cfg = load_config() outdir = os.path.join(cfg["paths"]["outputs"], "member_registry") tdir = os.path.join(cfg["paths"]["results"], "tables") fdir = os.path.join(cfg["paths"]["results"], "figures") os.makedirs(fdir, exist_ok=True) log.info("loading stored votes ...") data, members, blocks = load_votes(cfg) years = sorted(data) weeks = [b for b in blocks if b != "W1"] # prune on operational range W2-W6 M = len(members) rank_rows, eff_rows = [], [] drop_count = {(w, mid): 0 for w in weeks for mid in members} # nested consensus # nested honest: for each week, pooled held-out ensemble PSS (full and fold-specific pruned) honest = {w: {"full": [], "pruned": []} for w in weeks} for w in weeks: for y in years: tr = [t for t in years if t != y] Vtr = np.concatenate([data[t][w]["votes"] for t in tr]) ytr = np.concatenate([data[t][w]["y"] for t in tr]) cand, pss, far, full_pss_tr, pruned_pss_tr = prune_decision(Vtr, ytr, members) # record ranking (train) order = np.argsort(-pss) for rk, j in enumerate(order): rank_rows.append(dict(week=w, held_out=y, member_id=members[j], train_PSS=round(float(pss[j]), 4), train_FAR=round(float(far[j]), 4), rank=rk + 1, dropped=int(j in cand))) for j in cand: drop_count[(w, members[j])] += 1 # honest OUT-OF-SAMPLE eval on held-out y with the train-derived P0 Vte = data[y][w]["votes"]; yte = data[y][w]["y"] keep_pruned = np.array([j for j in range(M) if j not in cand]) _, p0_full = _peirce_from_votefrac(ytr, Vtr.mean(1)) _, p0_prun = _peirce_from_votefrac(ytr, Vtr[:, keep_pruned].mean(1)) honest[w]["full"].append( (yte, (Vte.mean(1) >= p0_full).astype(int))) honest[w]["pruned"].append( (yte, (Vte[:, keep_pruned].mean(1) >= p0_prun).astype(int))) # ---- consensus registry v2 (per week) + overall ---- cons_rows = [] overall_drop = {mid: 0 for mid in members} # counts over all week-fold decisions for w in weeks: for mid in members: times = drop_count[(w, mid)] cons_rows.append(dict(week=w, member_id=mid, times_dropped=times, n_folds=len(years), frac_dropped=round(times / len(years), 3), consensus_pruned=int(times / len(years) >= CONSENSUS))) overall_drop[mid] += times cons = pd.DataFrame(cons_rows) n_decisions = len(years) * len(weeks) # per-week consensus: member pruned AT week w if dropped in >= CONSENSUS of that week's folds pruned_at = {mid: sorted([w for w in weeks if drop_count[(w, mid)] / len(years) >= CONSENSUS]) for mid in members} # active-by-week mask (W1 never pruned; keep all members there) active_rows = [] for mid in members: for b in blocks: act = 0 if (b in pruned_at[mid]) else 1 active_rows.append(dict(member_id=mid, block=b, active=act)) pd.DataFrame(active_rows).to_csv(os.path.join(tdir, "v2_active_members_by_week.csv"), index=False) overall_status = {mid: ("pruned@" + ",".join(pruned_at[mid]) if pruned_at[mid] else "active") for mid in members} # ---- nested honest vs optimistic (optimism gap) ---- def pooled_peirce(pairs): a=b=c=d=0 for yv, fl in pairs: a += int(np.sum((fl == 1) & (yv == 1))); b += int(np.sum((fl == 1) & (yv == 0))) c += int(np.sum((fl == 0) & (yv == 1))); d += int(np.sum((fl == 0) & (yv == 0))) return vf.binary_scores_from_counts(a, b, c, d)["Peirce"] for w in weeks: # optimistic: one consensus set chosen on ALL years, scored on ALL years Vall = np.concatenate([data[t][w]["votes"] for t in years]) yall = np.concatenate([data[t][w]["y"] for t in years]) cand_all, _, _, _, _ = prune_decision(Vall, yall, members) keep_all = np.array([j for j in range(M) if j not in cand_all]) opt_full, p0 = _peirce_from_votefrac(yall, Vall.mean(1)) opt_prun, p0p = _peirce_from_votefrac(yall, Vall[:, keep_all].mean(1)) opt_prun_score = vf.binary_scores(yall, Vall[:, keep_all].mean(1), p0p)["Peirce"] h_full = pooled_peirce(honest[w]["full"]) h_prun = pooled_peirce(honest[w]["pruned"]) eff_rows.append(dict( week=w, n_members=M, full_pool_PSS_nested=round(h_full, 4), pruned_pool_PSS_nested=round(h_prun, 4), nested_delta=round(h_prun - h_full, 4), optimistic_pruned_PSS=round(opt_prun_score, 4), optimism_gap=round(opt_prun_score - h_prun, 4), n_dropped_allyear=int(len(cand_all)))) pd.DataFrame(rank_rows).to_csv(os.path.join(tdir, "v2_pruning_ranking.csv"), index=False) cons.to_csv(os.path.join(tdir, "v2_pruning_consensus.csv"), index=False) eff = pd.DataFrame(eff_rows) eff.to_csv(os.path.join(tdir, "v2_pruning_ensemble_effect.csv"), index=False) # ---- write registry v2 ---- reg1 = json.load(open(os.path.join(outdir, "registry_v1.json"))) for m in reg1["members"]: mid = m["member_id"] st = overall_status[mid] m["status"] = st wk = cons[cons.member_id == mid] m["status_reason"] = ( f"Stage-19 nested per-week prune: consensus-pruned at [{','.join(pruned_at[mid]) or 'none'}] " f"(>= {int(CONSENSUS*100)}% of folds); active at the other weeks.") m["registry_version"] = "v2" m["active_by_week"] = {b: (0 if b in pruned_at[mid] else 1) for b in blocks} m["pruned_at_weeks"] = pruned_at[mid] reg1["registry_version"] = "v2" reg1["description"] = ("Registry v2 = v1 after nested leakage-safe pruning (stage 19). Selection " "is inside the outer LOYO and PER LEAD WEEK; a member may be pruned at one " "week and kept at another. Downstream items use v2_active_members_by_week.csv.") reg1["n_member_week_prunes"] = int(sum(len(pruned_at[mid]) for mid in members)) reg1["n_members_pruned_any_week"] = sum(1 for mid in members if pruned_at[mid]) reg1["pruning_criterion"] = (f"PSS(POD-POFD) per week W2-W6 pooled over regions, FAR ceiling " f"{FAR_CEILING}, bottom {int(BOTTOM_FRAC*100)}% candidates, " f"ensemble-effect gate eps={EPS}, min_keep={MIN_KEEP}.") with open(os.path.join(outdir, "registry_v2.json"), "w") as f: json.dump(reg1, f, indent=2) pd.DataFrame([dict(member_id=m["member_id"], family=m["family"], cascade_cutoff=m["configuration"]["cascade_cutoff"], status=m["status"], status_reason=m["status_reason"]) for m in reg1["members"]]).to_csv( os.path.join(outdir, "registry_v2.csv"), index=False) # ---- figure: mean train-PSS ranking with drop consensus ---- import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt rr = pd.DataFrame(rank_rows) order_mid = (rr.groupby("member_id")["train_PSS"].mean().sort_values(ascending=False).index.tolist()) mean_pss = rr.groupby("member_id")["train_PSS"].mean().reindex(order_mid) frac_drop = (cons.groupby("member_id")["frac_dropped"].mean().reindex(order_mid)) fig, ax = plt.subplots(figsize=(10, 6)) colors = ["#b42318" if pruned_at[m] else "#1a7f37" for m in order_mid] ax.barh(range(len(order_mid)), mean_pss.values, color=colors) ax.set_yticks(range(len(order_mid))); ax.set_yticklabels(order_mid, fontsize=8) ax.invert_yaxis(); ax.set_xlabel("mean train PSS (W2-W6, nested)") ax.set_title("Stage-19 member ranking — red = pruned at >=1 lead week (>=50% fold consensus)") for i, m in enumerate(order_mid): ax.text(mean_pss[m] + 0.001, i, f"{frac_drop[m]*100:.0f}%", va="center", fontsize=7) fig.tight_layout(); fig.savefig(os.path.join(fdir, "v2_pruning_ranking.png"), dpi=130) log.info("registry v2: %d member-week prunes across %d members", reg1["n_member_week_prunes"], reg1["n_members_pruned_any_week"]) print("\n=== per-week consensus prunes (>=50% of folds) ===") for mid in members: if pruned_at[mid]: print(f" {mid}: pruned at {pruned_at[mid]}") print("\n=== ensemble effect (nested honest) + optimism gap ===") print(eff.to_string(index=False)) if __name__ == "__main__": main()