"""Multi-model ensemble / stacking of the v2 zoo -> EHD_ens forecast. The reforecast is control-only (one physical member), so the "members" combined here are the AI model families of the v2 zoo, each already an out-of-sample, leakage-safe, isotonic- calibrated P(yEHDE). Two combinations are produced: ens_mean unweighted mean of the member probabilities (no fitting; robust baseline). ens_stack a logistic meta-learner on the member probabilities, trained LEAVE-ONE-YEAR-OUT at the meta level: for held-out year y the stacker is fit on the OTHER 19 years' (member_probs -> y) and applied to y. Because the member probabilities for y were themselves produced by base models that never saw y (stage 12), and the meta model for y is fit without y, the stacked prediction for y stays leakage-safe. Members are taken at cascade cutoff 1.0 (single-level, no screening) because the cascade screen was shown ~neutral (MEMORY L-12); this keeps the ensemble members clean and comparable. The per-region decision threshold for an ensemble is fit, per held-out year, on the ensemble probabilities of the fit years (the 19 non-held-out years) — the analogue of the zoo's train-only threshold — then applied to the held-out year (leakage-safe). """ from __future__ import annotations import numpy as np from sklearn.linear_model import LogisticRegression from sklearn.preprocessing import StandardScaler from sklearn.pipeline import make_pipeline from . import verify as vf MEMBERS = ["logistic", "hgb", "rf", "mlp", "ridge_sev", "hgb_sev", "mlp_sev"] EXTRA_MEMBERS = ["cnn", "ann_torch"] # torch families (Phase B), merged if their npz exist MEMBER_TAG = "c1" # cascade cutoff 1.0 columns (no screening) ENS_TAGS = ["ensmean_c1", "ensstack_c1"] def member_cols(members=MEMBERS): return [f"p_{m}_{MEMBER_TAG}" for m in members] def load_by_year(files, block_ids, members=MEMBERS, extra_files=None, extra_members=None): """Read per-year npz into a nested dict: data[year][block] = dict(M, y, cell, region). M is [n_samples, n_members] of member probabilities; block is the block *name*. `files` maps year -> path; `block_ids` maps block name -> integer id. If `extra_files` (year -> path) is given, the `extra_members` columns from those npz are appended to M after asserting the (block, cell) ordering matches the main file (stages 12/12b assemble identically). Returns (data, members_used). """ cols = member_cols(members) extra_members = extra_members or [] ecols = member_cols(extra_members) out = {} used = list(members) have_extra = bool(extra_files) and extra_members for yr, path in files.items(): d = np.load(path, allow_pickle=True) bid = d["block_id"].astype(int) de = None if have_extra and yr in extra_files: de = np.load(extra_files[yr], allow_pickle=True) ebid = de["block_id"].astype(int) out[yr] = {} for bname, bcode in block_ids.items(): bm = bid == bcode M = np.column_stack([d[c][bm].astype(float) for c in cols]) cell = d["cell"][bm].astype(int) if de is not None: ebm = ebid == bcode assert np.array_equal(cell, de["cell"][ebm].astype(int)), \ f"cnn/main cell ordering mismatch {yr} {bname}" Mx = np.column_stack([de[c][ebm].astype(float) for c in ecols]) M = np.column_stack([M, Mx]) out[yr][bname] = dict(M=M, y=d["y"][bm].astype(float), cell=cell, region=d["region"][bm].astype(int)) if have_extra: used = list(members) + list(extra_members) return out, used def _fit_region_thresholds(y, p, region, region_ids, metric="Peirce"): g_thr, _ = vf.best_threshold(y, p, metric) out = {int(r): g_thr for r in region_ids} for r in region_ids: m = region == r if m.sum() > 200 and np.nansum(y[m]) > 10: out[int(r)], _ = vf.best_threshold(y[m], p[m], metric) return out def build_ensembles(data, block_names, region_ids, years): """LOYO mean + stacked ensemble per block. Returns per (block) a dict with pooled arrays across all years (in `years` order): y, cell, region, p_mean, p_stack, flag_mean, flag_stack plus decthr_rows: list of (tag, block, region, thr) fit per held-out year (averaged? no — each year uses its own fold threshold; we store the per-fold thresholds only for provenance). Flags are already assigned per sample using that sample's held-out-year threshold. """ result = {} decthr_rows = [] for b in block_names: ys, cells, regs = [], [], [] pmean, pstack, fmean, fstack = [], [], [], [] for yr in years: te = data[yr][b] tr_years = [t for t in years if t != yr] # --- mean ensemble --- pm_te = np.nanmean(te["M"], axis=1) # --- stacked ensemble: fit meta on other years --- Mtr = np.concatenate([data[t][b]["M"] for t in tr_years]) ytr = np.concatenate([data[t][b]["y"] for t in tr_years]) mtr = np.isfinite(ytr) & np.all(np.isfinite(Mtr), axis=1) meta = make_pipeline(StandardScaler(), LogisticRegression(C=1.0, max_iter=1000)) meta.fit(Mtr[mtr], ytr[mtr].astype(int)) ps_te = meta.predict_proba(np.nan_to_num(te["M"], nan=0.0))[:, 1] # --- decision thresholds: fit on the fit-years' ensemble probs (per region) --- pm_tr = np.nanmean(Mtr, axis=1) ps_tr = meta.predict_proba(np.nan_to_num(Mtr, nan=0.0))[:, 1] rtr = np.concatenate([data[t][b]["region"] for t in tr_years]) thr_mean = _fit_region_thresholds(ytr[mtr], pm_tr[mtr], rtr[mtr], region_ids) thr_stack = _fit_region_thresholds(ytr[mtr], ps_tr[mtr], rtr[mtr], region_ids) for r in region_ids: decthr_rows.append(("ensmean_c1", b, int(r), thr_mean[int(r)])) decthr_rows.append(("ensstack_c1", b, int(r), thr_stack[int(r)])) # apply per-sample threshold by region tvec_m = np.array([thr_mean[int(r)] for r in te["region"]]) tvec_s = np.array([thr_stack[int(r)] for r in te["region"]]) ys.append(te["y"]); cells.append(te["cell"]); regs.append(te["region"]) pmean.append(pm_te); pstack.append(ps_te) fmean.append(pm_te >= tvec_m); fstack.append(ps_te >= tvec_s) result[b] = dict( y=np.concatenate(ys), cell=np.concatenate(cells), region=np.concatenate(regs), p_mean=np.concatenate(pmean), p_stack=np.concatenate(pstack), flag_mean=np.concatenate(fmean), flag_stack=np.concatenate(fstack)) return result, decthr_rows