#!/usr/bin/env python """Stage 17 — vote-fraction ensemble + regional event table (prompt section 6). Implements the four sub-requirements of prompt "6. Additional requirements": (1) Compare the AI members forecasting EHDE (priority: high POD, low POFD/FAR, esp. W3-W6) and assess each member's contribution to the ensemble (leave-one-member-out). (2) Vote-fraction probability: P(yEHDE) = m/n, n = number of AI members (here 23), m = number voting yEHDE. Each member votes with its OWN leakage-safe per-region x block decision threshold (stored in the per-fold npz), so the vote fraction is out-of-sample. (3) Deterministic forecast: EHDE = 1 if P(yEHDE) > P0. P0 is chosen to MAXIMISE PEIRCE (POD-POFD -> rewards high POD + low POFD), selected LEAVE-ONE-YEAR-OUT (train years only) so the headline scores stay leakage-safe; the SEDI-optimal point is reported too, and the full P0 sweep curve is written for the forecaster to pick an operating point. (4) Regional event table: for region R_j, week W_k, and each init, P = M_j/N_j (fraction of the region's cells forecast yEHDE). A region is alerted if P >= P_thres. The same M/N rule is applied to the OBSERVED field to define the regional truth; P_thres is swept and chosen per week (Peirce, leakage-safe per year), with a fixed 0.5 (majority-of-cells) reference. The 23 "members" are all stored, leakage-safe prediction columns: 7 sklearn families {logistic,hgb,rf,mlp,ridge_sev,hgb_sev,mlp_sev} x 3 cascade cutoffs {c1,c0.9,c0.99} + cnn_c1 + ann_torch_c1. Outputs (results/tables, results/figures, outputs): v6_member_scores_by_block.csv per-member binary scores per block (part 1) v6_member_contribution.csv leave-one-member-out dPeirce/dSEDI on the vote ensemble v6_vote_scores_by_block.csv vote ensemble at per-fold-optimal P0 (parts 2-3) v6_p0_sweep.csv block x P0 -> POD/POFD/FAR/CSI/Peirce/SEDI (diagnostic) v6_regional_event_table.csv region x week contingency + skill (part 4) v6_regional_by_week.csv pooled-over-regions skill per week + chosen P_thres outputs/v6_vote_maps.nc per-cell vote-ensemble skill + mean probability figures: v6_p0_sweep.png, v6_member_pod_pofd_far.png, v6_vote_skillmap.png, v6_regional_event_heatmap.png """ from __future__ import annotations import os, sys, glob import numpy as np import pandas as pd import xarray as xr 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 HERE = os.path.dirname(os.path.abspath(__file__)) ROOT = os.path.dirname(HERE) TABLES = os.path.join(ROOT, "results", "tables") FIGS = os.path.join(ROOT, "results", "figures") FAMILIES = ["logistic", "hgb", "rf", "mlp", "ridge_sev", "hgb_sev", "mlp_sev"] CUTOFFS = ["c1", "c0.9", "c0.99"] TORCH = ["cnn_c1", "ann_torch_c1"] # member tag list (order fixed): 21 sklearn config columns + 2 torch MEMBERS = [f"{fam}_{cut}" for fam in FAMILIES for cut in CUTOFFS] + TORCH NMEM = len(MEMBERS) PLOT_METRICS = ["POD", "POFD", "FAR", "CSI", "EDI", "SEDI"] # ------------------------------------------------------------------ helpers def _thr_by_region(decthr, tag, block, region_ids): """Map region id -> stored decision threshold for (tag, block); default 0.5.""" out = np.full(int(max(region_ids)) + 1, 0.5) sel = (decthr["tag"] == tag) & (decthr["block"] == block) for row in decthr[sel]: out[int(row["region"])] = float(row["thr"]) return out def _counts(y, flag): y = y.astype(int); flag = flag.astype(bool) a = int(np.sum(flag & (y == 1))); b = int(np.sum(flag & (y == 0))) c = int(np.sum(~flag & (y == 1))); d = int(np.sum(~flag & (y == 0))) return a, b, c, d def _scores(a, b, c, d): return vf.binary_scores_from_counts(a, b, c, d) def _best_frac_threshold(y, frac, metric="Peirce"): """Pick threshold t maximising `metric` for the rule frac > t (strict). Candidate cuts are midpoints between the distinct vote fractions, so every achievable split is tested. Returns (t, score).""" vals = np.unique(frac) cands = np.concatenate([[-1.0], (vals[:-1] + vals[1:]) / 2.0, [1.0]]) if vals.size > 1 \ else np.array([-1.0, 1.0]) best_t, best_s = -1.0, -np.inf yi = y.astype(int) for t in cands: a, b, c, d = _counts(yi, frac > t) s = _scores(a, b, c, d)[metric] if np.isfinite(s) and s > best_s: best_s, best_t = s, t return best_t, best_s # ------------------------------------------------------------------ load / vote def build_pooled(cfg, block_names, block_ids, region_ids, years): """Single pass over the per-year npz: build the 23-member flag matrix per block and pool. Returns per-block dict of pooled arrays (y, region, cell, m, flags, yr, init) plus the per-(member,block) contingency counts for the member comparison table.""" main = {int(os.path.basename(p).split("_")[-1][:-4]): p for p in glob.glob(os.path.join(ROOT, "outputs", "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, "outputs", "loyo_pred_v2_cnn", "pred_v2cnn_*.npz"))} years = [y for y in years if y in main and y in cnn] pool = {b: {k: [] for k in ("y", "region", "cell", "m", "flags", "yr", "init")} for b in block_names} mcounts = {(tag, b): np.zeros(4, dtype=np.int64) for tag in MEMBERS for b in block_names} for yi, yr in enumerate(years): d = np.load(main[yr], allow_pickle=True) de = np.load(cnn[yr], allow_pickle=True) dt_main, dt_cnn = d["decthr"], de["decthr"] bid = d["block_id"].astype(int); ebid = de["block_id"].astype(int) for b in block_names: bm = bid == block_ids[b]; ebm = ebid == block_ids[b] reg = d["region"][bm].astype(int); cell = d["cell"][bm].astype(int) y = d["y"][bm].astype(float) assert np.array_equal(cell, de["cell"][ebm].astype(int)), f"cell mismatch {yr} {b}" n = y.size init = np.arange(n) // np.unique(cell).size # init-major tiling (104 x 541) fin = np.isfinite(y) flags = np.zeros((n, NMEM), dtype=np.int8) for mi, tag in enumerate(MEMBERS): if tag in TORCH: p = de["p_" + tag][ebm].astype(float); thr = _thr_by_region(dt_cnn, tag, b, region_ids) else: p = d["p_" + tag][bm].astype(float); thr = _thr_by_region(dt_main, tag, b, region_ids) fl = (p >= thr[reg]).astype(np.int8) flags[:, mi] = fl a, bb, c, dd = _counts(y[fin], fl[fin].astype(bool)) mcounts[(tag, b)] += (a, bb, c, dd) m = flags.sum(axis=1).astype(np.int8) pool[b]["y"].append(y[fin]); pool[b]["region"].append(reg[fin]) pool[b]["cell"].append(cell[fin]); pool[b]["m"].append(m[fin]) pool[b]["flags"].append(flags[fin]); pool[b]["init"].append(init[fin].astype(np.int16)) pool[b]["yr"].append(np.full(int(fin.sum()), yi, dtype=np.int16)) for b in block_names: for k in pool[b]: pool[b][k] = np.concatenate(pool[b][k]) return pool, mcounts, years # ------------------------------------------------------------------ part 1 def member_table(mcounts, block_names): rows = [] for tag in MEMBERS: for b in block_names: a, bb, c, d = mcounts[(tag, b)] s = _scores(int(a), int(bb), int(c), int(d)) rows.append(dict(member=tag, block=b, n=int(a + bb + c + d), base_rate=round(s["base_rate"], 3), POD=round(s["POD"], 3), POFD=round(s["POFD"], 3), FAR=round(s["FAR"], 3), CSI=round(s["CSI"], 3), Peirce=round(s["Peirce"], 3), SEDI=round(s["SEDI"], 3), freq_bias=round(s["freq_bias"], 3))) df = pd.DataFrame(rows) df.to_csv(os.path.join(TABLES, "v6_member_scores_by_block.csv"), index=False) print(f"[m] wrote v6_member_scores_by_block.csv ({len(df)} rows)") return df def lomo_contribution(pool, block_names): """Leave-one-member-out: drop each member, re-optimise the vote P0 (pooled), and report the change in Peirce/SEDI. Positive dPeirce => the member HELPS the ensemble.""" rows = [] for b in block_names: y = pool[b]["y"]; m = pool[b]["m"].astype(np.int32); flags = pool[b]["flags"] _, full_pe = _best_frac_threshold(y, m / NMEM, "Peirce") _, full_se = _best_frac_threshold(y, m / NMEM, "SEDI") for mi, tag in enumerate(MEMBERS): mk = m - flags[:, mi] # drop member -> n-1 votes _, pe = _best_frac_threshold(y, mk / (NMEM - 1), "Peirce") _, se = _best_frac_threshold(y, mk / (NMEM - 1), "SEDI") rows.append(dict(member=tag, block=b, dPeirce=round(full_pe - pe, 4), dSEDI=round(full_se - se, 4))) df = pd.DataFrame(rows) df.to_csv(os.path.join(TABLES, "v6_member_contribution.csv"), index=False) print(f"[m] wrote v6_member_contribution.csv ({len(df)} rows)") return df # ------------------------------------------------------------------ parts 2-3 def p0_sweep_and_forecast(pool, block_names): """Diagnostic P0 sweep (pooled) + leakage-safe per-fold-optimal P0 headline forecast. Returns (sweep_df, vote_df, det_flag_by_block) where det_flag_by_block[b] is the per-sample deterministic EHDE from the per-fold P0 (used by the regional table).""" sweep_rows, vote_rows = [], [] det_flag = {} for b in block_names: y = pool[b]["y"]; m = pool[b]["m"].astype(np.int32); yr = pool[b]["yr"] frac = m / NMEM # --- diagnostic pooled sweep over the achievable P0 = k/NMEM --- for k in range(NMEM + 1): t = k / NMEM - 1e-9 # frac > t <=> m/N >= k/N a, bb, c, d = _counts(y, frac > t) s = _scores(a, bb, c, d) sweep_rows.append(dict(block=b, P0=round(k / NMEM, 4), m_min=k, POD=round(s["POD"], 3), POFD=round(s["POFD"], 3), FAR=round(s["FAR"], 3), CSI=round(s["CSI"], 3), Peirce=round(s["Peirce"], 3), SEDI=round(s["SEDI"], 3))) # --- pooled-optimal reference (Peirce & SEDI) --- t_pe, _ = _best_frac_threshold(y, frac, "Peirce") t_se, _ = _best_frac_threshold(y, frac, "SEDI") # --- leakage-safe per-fold P0 (Peirce on train years) -> honest deterministic forecast --- det = np.zeros(y.size, dtype=bool) p0_used = [] for yhold in np.unique(yr): tr = yr != yhold; te = yr == yhold t_fold, _ = _best_frac_threshold(y[tr], frac[tr], "Peirce") det[te] = frac[te] > t_fold p0_used.append((t_fold + 1e-9)) det_flag[b] = det a, bb, c, d = _counts(y, det) s = _scores(a, bb, c, d) vote_rows.append(dict( block=b, n=int(y.size), base_rate=round(float(np.mean(y)), 3), P0_peirce_pooled=round(max(t_pe, 0.0) + 1e-9, 3), P0_sedi_pooled=round(max(t_se, 0.0) + 1e-9, 3), P0_perfold_mean=round(float(np.mean(p0_used)), 3), POD=round(s["POD"], 3), POFD=round(s["POFD"], 3), FAR=round(s["FAR"], 3), CSI=round(s["CSI"], 3), Peirce=round(s["Peirce"], 3), SEDI=round(s["SEDI"], 3), HSS=round(s["HSS"], 3), freq_bias=round(s["freq_bias"], 3))) sweep_df = pd.DataFrame(sweep_rows); vote_df = pd.DataFrame(vote_rows) sweep_df.to_csv(os.path.join(TABLES, "v6_p0_sweep.csv"), index=False) vote_df.to_csv(os.path.join(TABLES, "v6_vote_scores_by_block.csv"), index=False) print(f"[p0] wrote v6_p0_sweep.csv ({len(sweep_df)}), v6_vote_scores_by_block.csv " f"({len(vote_df)})") return sweep_df, vote_df, det_flag # ------------------------------------------------------------------ part 4 def regional_table(pool, det_flag, block_names, region_ids): """Regional M_j/N_j event forecast, verified against the same rule on the observed field. Per (year, init) and region: P_fore = M_fore/N, P_obs = M_obs/N. A common P_thres is swept; the region is 'alerted' if P >= P_thres. P_thres is chosen per week (block) to maximise Peirce on train years (leakage-safe), and a fixed 0.5 reference is also reported.""" reg_rows, week_rows = [], [] for b in block_names: y = pool[b]["y"]; reg = pool[b]["region"]; yr = pool[b]["yr"]; init = pool[b]["init"] fore = det_flag[b].astype(float) g = yr.astype(np.int64) * 1000 + init # unique (year, init) group id gids, ginv = np.unique(g, return_inverse=True) gyr = np.zeros(gids.size, dtype=np.int64) gyr[ginv] = yr # year label per group # per region: coverage fractions per group Pf = {}; Po = {} for r in region_ids: sel = reg == r if sel.sum() == 0: continue N = np.bincount(ginv[sel], minlength=gids.size).astype(float) Mf = np.bincount(ginv[sel], weights=fore[sel], minlength=gids.size) Mo = np.bincount(ginv[sel], weights=y[sel], minlength=gids.size) ok = N > 0 Pf[r] = np.where(ok, Mf / np.maximum(N, 1), np.nan) Po[r] = np.where(ok, Mo / np.maximum(N, 1), np.nan) # --- choose P_thres per week (Peirce, leakage-safe) over pooled regions --- allf = np.concatenate([Pf[r] for r in Pf]); allo = np.concatenate([Po[r] for r in Po]) allyr = np.concatenate([gyr for _ in Pf]) fin = np.isfinite(allf) & np.isfinite(allo) allf, allo, allyr = allf[fin], allo[fin], allyr[fin] grid = np.round(np.arange(0.05, 1.0, 0.05), 2) def obs_event(po, t): # symmetric rule on observed field return po >= t # per-fold selection of a single P_thres det_alert = np.zeros(allf.size, dtype=bool) det_truth = np.zeros(allf.size, dtype=bool) pth_used = [] for yhold in np.unique(allyr): tr = allyr != yhold; te = allyr == yhold best_t, best_s = 0.5, -np.inf for t in grid: a, bb, c, d = _counts(obs_event(allo[tr], t), allf[tr] >= t) s = _scores(a, bb, c, d)["Peirce"] if np.isfinite(s) and s > best_s: best_s, best_t = s, t det_alert[te] = allf[te] >= best_t det_truth[te] = allo[te] >= best_t pth_used.append(best_t) a, bb, c, d = _counts(det_truth, det_alert) s = _scores(a, bb, c, d) # fixed 0.5 reference a5, b5, c5, d5 = _counts(allo >= 0.5, allf >= 0.5) s5 = _scores(a5, b5, c5, d5) week_rows.append(dict( block=b, n_region_inits=int(allf.size), Pthres_perfold_mean=round(float(np.mean(pth_used)), 3), POD=round(s["POD"], 3), POFD=round(s["POFD"], 3), FAR=round(s["FAR"], 3), CSI=round(s["CSI"], 3), Peirce=round(s["Peirce"], 3), SEDI=round(s["SEDI"], 3), POD_at0p5=round(s5["POD"], 3), FAR_at0p5=round(s5["FAR"], 3), Peirce_at0p5=round(s5["Peirce"], 3))) # per region x week at the mean chosen P_thres (descriptive) tstar = float(np.mean(pth_used)) for r in Pf: fmask = np.isfinite(Pf[r]) & np.isfinite(Po[r]) a, bb, c, d = _counts(Po[r][fmask] >= tstar, Pf[r][fmask] >= tstar) s = _scores(a, bb, c, d) reg_rows.append(dict(region=int(r), block=b, Pthres=round(tstar, 3), n_inits=int(fmask.sum()), hits=a, false_alarms=bb, misses=c, POD=round(s["POD"], 3), POFD=round(s["POFD"], 3), FAR=round(s["FAR"], 3), CSI=round(s["CSI"], 3), Peirce=round(s["Peirce"], 3), SEDI=round(s["SEDI"], 3))) reg_df = pd.DataFrame(reg_rows); week_df = pd.DataFrame(week_rows) reg_df.to_csv(os.path.join(TABLES, "v6_regional_event_table.csv"), index=False) week_df.to_csv(os.path.join(TABLES, "v6_regional_by_week.csv"), index=False) print(f"[reg] wrote v6_regional_event_table.csv ({len(reg_df)}), v6_regional_by_week.csv " f"({len(week_df)})") return reg_df, week_df # ------------------------------------------------------------------ maps (part 2) def _grid_geometry(cfg): ref = xr.open_dataarray(os.path.join(cfg["paths"]["outputs"], "forecast_daily_tx", f"fc_daily_tx_{cfg['forecast']['init_years'][0]}.nc")) reg = xr.open_dataset(cfg["paths"]["mask"])["reg"].reindex( lat=ref.lat, lon=ref.lon, method="nearest") iy, ix = np.where(np.isfinite(reg.values)) return ref.lat.values, ref.lon.values, iy, ix def vote_maps(cfg, pool, det_flag, block_names): lat, lon, iy, ix = _grid_geometry(cfg) ncell = iy.size metrics = ["POD", "POFD", "FAR", "CSI", "EDI", "SEDI"] arr = {mt: np.full((len(block_names), lat.size, lon.size), np.nan, np.float32) for mt in metrics} meanp = np.full((len(block_names), lat.size, lon.size), np.nan, np.float32) for bi, b in enumerate(block_names): y = pool[b]["y"]; cell = pool[b]["cell"]; det = det_flag[b].astype(bool) frac = pool[b]["m"].astype(float) / NMEM a = np.bincount(cell, weights=(det & (y == 1)), minlength=ncell)[:ncell] bb = np.bincount(cell, weights=(det & (y == 0)), minlength=ncell)[:ncell] c = np.bincount(cell, weights=(~det & (y == 1)), minlength=ncell)[:ncell] d = np.bincount(cell, weights=(~det & (y == 0)), minlength=ncell)[:ncell] mp = np.bincount(cell, weights=frac, minlength=ncell)[:ncell] cnt = np.bincount(cell, minlength=ncell)[:ncell] ev = a + c with np.errstate(divide="ignore", invalid="ignore"): POD = np.where(ev > 0, a / ev, np.nan) POFD = np.where((bb + d) > 0, bb / (bb + d), np.nan) FAR = np.where((a + bb) > 0, bb / (a + bb), np.nan) CSI = np.where((a + bb + c) > 0, a / (a + bb + c), np.nan) F = bb / (bb + d + 1e-12); H = a / (ev + 1e-12) e = 1e-12 edi = (np.log(F + e) - np.log(H + e)) / (np.log(F + e) + np.log(H + e) + e) sedi = ((np.log(F + e) - np.log(H + e) - np.log(1 - F + e) + np.log(1 - H + e)) / (np.log(F + e) + np.log(H + e) + np.log(1 - F + e) + np.log(1 - H + e) + e)) ok = (ev > 0) & ((bb + d) > 0) vals = dict(POD=POD, POFD=POFD, FAR=FAR, CSI=CSI, EDI=np.where(ok, edi, np.nan), SEDI=np.where(ok, sedi, np.nan)) for mt in metrics: g = np.full((lat.size, lon.size), np.nan); g[iy, ix] = vals[mt]; arr[mt][bi] = g gp = np.full((lat.size, lon.size), np.nan) gp[iy, ix] = np.where(cnt > 0, mp / np.maximum(cnt, 1), np.nan); meanp[bi] = gp ds = xr.Dataset({mt: (("block", "lat", "lon"), arr[mt]) for mt in metrics}, coords=dict(block=block_names, lat=lat, lon=lon)) ds["mean_prob"] = (("block", "lat", "lon"), meanp) ds.attrs["description"] = ("Vote-fraction ensemble (m/n, 23 members): per-cell binary skill " "of the deterministic forecast (per-fold P0) and mean P(yEHDE).") outnc = os.path.join(cfg["paths"]["outputs"], "v6_vote_maps.nc") ds.to_netcdf(outnc); print(f"[map] wrote {outnc}") _skillmap_fig(ds, block_names) # ------------------------------------------------------------------ figures def _skillmap_fig(ds, block_names): import matplotlib; matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.colors import BoundaryNorm, LinearSegmentedColormap, Normalize levels = [round(0.1 * i, 1) for i in range(11)] div_cmap = plt.get_cmap("RdYlGn_r"); div_norm = Normalize(0.0, 1.0) seq_cmap = LinearSegmentedColormap.from_list("yellow_red", ["#FFFF00", "#D00000"]) seq_cmap.set_under("#FFFF00"); seq_norm = BoundaryNorm(levels, ncolors=seq_cmap.N) styles = dict(POD=(div_cmap, div_norm), POFD=(div_cmap, div_norm), FAR=(div_cmap, div_norm), CSI=(seq_cmap, seq_norm), EDI=(seq_cmap, seq_norm), SEDI=(seq_cmap, seq_norm)) nr, nc = len(PLOT_METRICS), len(block_names) fig, axes = plt.subplots(nr, nc, figsize=(2.05 * nc, 1.95 * nr), squeeze=False) for ri, mt in enumerate(PLOT_METRICS): cmap, norm = styles[mt] for ci, b in enumerate(block_names): ax = axes[ri][ci] im = ax.pcolormesh(ds.lon, ds.lat, ds[mt].sel(block=b).values, cmap=cmap, norm=norm, shading="auto") ax.set_xticks([]); ax.set_yticks([]) if ri == 0: ax.set_title(b, fontsize=9) if ci == 0: ax.set_ylabel(mt, fontsize=9) cax = fig.add_axes([0.915, axes[ri][-1].get_position().y0, 0.012, axes[ri][-1].get_position().height]) fig.colorbar(im, cax=cax, ticks=levels) fig.suptitle("Vote-fraction ensemble (m/n, 23 members) — per-cell binary skill", fontsize=12) fig.subplots_adjust(left=0.05, right=0.9, top=0.93, bottom=0.03, wspace=0.05, hspace=0.12) out = os.path.join(FIGS, "v6_vote_skillmap.png"); fig.savefig(out, dpi=125); plt.close(fig) print(f"[map] wrote {out}") def _p0_sweep_fig(sweep_df, vote_df, block_names): import matplotlib; matplotlib.use("Agg") import matplotlib.pyplot as plt fig, axes = plt.subplots(2, 3, figsize=(13, 7), squeeze=False) for i, b in enumerate(block_names): ax = axes[i // 3][i % 3] s = sweep_df[sweep_df.block == b] for col, c in (("POD", "tab:green"), ("POFD", "tab:blue"), ("FAR", "tab:orange"), ("Peirce", "k"), ("SEDI", "tab:red")): ax.plot(s["P0"], s[col], marker=".", ms=4, lw=1.4, color=c, label=col) p0 = float(vote_df[vote_df.block == b]["P0_peirce_pooled"].iloc[0]) ax.axvline(p0, color="0.5", ls="--", lw=1) ax.set_title(f"{b} (P0*≈{p0:.2f})"); ax.set_xlabel("P0 = m/n threshold"); ax.grid(alpha=.3) ax.set_ylim(-0.1, 1.02) if i == 0: ax.legend(fontsize=8, ncol=2) fig.suptitle("Deterministic yEHDE vs vote threshold P0 (pooled LOYO)", fontsize=13) fig.tight_layout(rect=[0, 0, 1, 0.97]) out = os.path.join(FIGS, "v6_p0_sweep.png"); fig.savefig(out, dpi=130); plt.close(fig) print(f"[fig] wrote {out}") def _member_fig(mdf, block_names): import matplotlib; matplotlib.use("Agg") import matplotlib.pyplot as plt fig, axes = plt.subplots(1, 3, figsize=(15, 5), squeeze=False) for j, mt in enumerate(["POD", "POFD", "FAR"]): ax = axes[0][j] piv = mdf.pivot(index="member", columns="block", values=mt).reindex(MEMBERS)[block_names] im = ax.imshow(piv.values, aspect="auto", cmap="RdYlGn" if mt == "POD" else "RdYlGn_r", vmin=0, vmax=1) ax.set_xticks(range(len(block_names))); ax.set_xticklabels(block_names) ax.set_yticks(range(len(MEMBERS))); ax.set_yticklabels(MEMBERS, fontsize=7) ax.set_title(mt); fig.colorbar(im, ax=ax, fraction=0.046) fig.suptitle("Member skill by lead block (green=good): POD high-good, POFD/FAR low-good", fontsize=12) fig.tight_layout(rect=[0, 0, 1, 0.96]) out = os.path.join(FIGS, "v6_member_pod_pofd_far.png"); fig.savefig(out, dpi=125); plt.close(fig) print(f"[fig] wrote {out}") def _regional_fig(reg_df, block_names, region_ids): import matplotlib; matplotlib.use("Agg") import matplotlib.pyplot as plt piv = reg_df.pivot(index="region", columns="block", values="Peirce").reindex( sorted(region_ids))[block_names] fig, ax = plt.subplots(figsize=(7, 4.5)) im = ax.imshow(piv.values, aspect="auto", cmap="RdYlGn", vmin=-0.2, vmax=1.0) ax.set_xticks(range(len(block_names))); ax.set_xticklabels(block_names) ax.set_yticks(range(len(piv.index))); ax.set_yticklabels([f"R{r}" for r in piv.index]) for iy_ in range(piv.shape[0]): for ix_ in range(piv.shape[1]): v = piv.values[iy_, ix_] if np.isfinite(v): ax.text(ix_, iy_, f"{v:.2f}", ha="center", va="center", fontsize=8) ax.set_title("Regional event forecast — Peirce (region × week)") fig.colorbar(im, ax=ax, fraction=0.046, label="Peirce") fig.tight_layout() out = os.path.join(FIGS, "v6_regional_event_heatmap.png"); fig.savefig(out, dpi=130) plt.close(fig); print(f"[fig] wrote {out}") # ------------------------------------------------------------------ main def main(): cfg = load_config() block_names = list(lead_blocks(cfg).keys()) block_ids = {b: k + 1 for k, b in enumerate(block_names)} region_ids = list(cfg["regions"]["ids"]) years = list(cfg["cross_validation"]["years"]) os.makedirs(TABLES, exist_ok=True); os.makedirs(FIGS, exist_ok=True) print(f"[17] {NMEM} members: {MEMBERS}") pool, mcounts, used_years = build_pooled(cfg, block_names, block_ids, region_ids, years) print(f"[17] pooled {len(used_years)} folds; W1 samples={pool[block_names[0]]['y'].size}") mdf = member_table(mcounts, block_names) cdf = lomo_contribution(pool, block_names) sweep_df, vote_df, det_flag = p0_sweep_and_forecast(pool, block_names) reg_df, week_df = regional_table(pool, det_flag, block_names, region_ids) vote_maps(cfg, pool, det_flag, block_names) _p0_sweep_fig(sweep_df, vote_df, block_names) _member_fig(mdf, block_names) _regional_fig(reg_df, block_names, region_ids) # headline console summary print("\n[17] vote ensemble (per-fold P0) vs base rate, by block:") for _, r in vote_df.iterrows(): print(f" {r.block}: POD={r.POD:.3f} POFD={r.POFD:.3f} FAR={r.FAR:.3f} " f"Peirce={r.Peirce:.3f} SEDI={r.SEDI:.3f} (P0~{r.P0_perfold_mean:.2f})") print("\n[17] top LOMO contributors (mean dPeirce over W3-W6):") w = cdf[cdf.block.isin(["W3", "W4", "W5", "W6"])].groupby("member")["dPeirce"].mean() for tag, v in w.sort_values(ascending=False).head(6).items(): print(f" {tag:16s} dPeirce={v:+.4f}") print("[17] done.") if __name__ == "__main__": main()