#!/usr/bin/env python3 """EXP2: a two-step forecast -- first ask whether an EHDE is possible at all. The user's proposal. Split the sample into a group where the region has an EHDE (yEHD) and a group where it does not (nEHD), using the EXP1 rule that a region with P_i = m_i/n_i below 5 % counts as nEHD for that week. Operationally the forecast then runs in two steps: 1. a gate answers "can this region have an EHDE in this week?" If the probability is near zero, no grid-point forecast is issued at all; 2. only where the gate opens does the grid-point EHD forecast run. Two things have to be said plainly about this design. First, step 1 already exists. "Region with P_i < 5 %" is exactly the regional EHDE label at coverage 5 % that section AA models, so the gate is the stage-2 logistic regression, not a new object. What is new is using it to *suppress* the point forecast rather than to score the region. Second, the evaluation must not be conditioned on the gate. Scoring the point forecast only where the gate opened is the mistake section S records for the `daily_capable` frame: it conditions verification on the forecast itself and flatters the result. Every table here is computed over ALL validation rows; a gate that closes on a day that did have events simply turns those into misses. That makes the tables comparable with every other point-level table in the project, and it makes the cost of gating visible instead of hidden. What the cascade can and cannot do: P(EHD at point) = P(EHDE in region) x P(EHD at point | EHDE in region) is a factorisation, not extra information, so with perfectly calibrated stages the product is what a single correct model would already give. The gain, if there is one, comes from the second factor being an easier problem to fit -- a homogeneous subpopulation with a much higher base rate -- which is the same effect that made a linear model beat the trees in section AA. """ from __future__ import annotations import argparse import json import os import sys import numpy as np sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import ehd_common as C import ehd_metrics as M from train_national import RES_DIR import regional_model as RM RM.KEYS = RM.KEYS + ("point",) # the cascade joins rows on the grid point too load = RM.load def gate_key(init, lt, region): return (np.asarray(init, "i8") * C.NLT + np.asarray(lt, "i8")) * 8 + region def join_rows(src, dst): """Score of `src` for each row of `dst`, joined on (init, lead, region, point).""" ks = gate_key(src["init"], src["lt"], src["region"]) * 1000 + src["point"] kd = gate_key(dst["init"], dst["lt"], dst["region"]) * 1000 + dst["point"] o = np.argsort(ks, kind="stable") pos = np.clip(np.searchsorted(ks[o], kd), 0, ks.size - 1) return np.where(ks[o][pos] == kd, src["score"][o][pos], np.nan) def join_gate(d, keys, vals): """Gate probability of each point row, looked up by its (init, lead, region).""" o = np.argsort(keys, kind="stable") ks, vs = keys[o], vals[o] pos = np.searchsorted(ks, gate_key(d["init"], d["lt"], d["region"])) pos = np.clip(pos, 0, ks.size - 1) hit = ks[pos] == gate_key(d["init"], d["lt"], d["region"]) out = np.full(d["score"].size, np.nan) out[hit] = vs[pos[hit]] return out, hit def weekly_fbi1(score, obs, week): """Threshold per week at the climatological event frequency of that week.""" return {w: (float(np.quantile(score[week == w], 1.0 - obs[week == w].mean())) if (week == w).sum() and obs[week == w].sum() else np.inf) for w in range(len(C.WEEKS))} def apply_weekly(score, week, th): p = np.zeros(score.size, bool) for w, t in th.items(): m = week == w p[m] = score[m] >= t return p def table(pred, obs, week, score, title, weeks=("W4", "W5", "W6", "ALL")): rows = {} print("\n %s" % title) print(" %9s %8s %7s %7s %7s %7s %7s %7s" % ("base_rate", "n_evt", "POD", "POFD", "FAR", "FBI", "CSI", "AUC")) for w in list(range(len(C.WEEKS))) + ["ALL"]: name = C.WEEK_NAMES[w] if w != "ALL" else "ALL" m = (week == w) if w != "ALL" else np.ones(obs.size, bool) if m.sum() == 0 or obs[m].sum() == 0: continue s = M.score_binary(pred[m], obs[m]) s["AUC"] = M.auc_score(score[m], obs[m]) s["n_events"] = int(obs[m].sum()) rows[name] = s if name in weeks: print(" %-4s %9.4f %8d %7.3f %7.3f %7.3f %7.3f %7.3f %7.3f" % (name, s["base_rate"], s["n_events"], s["POD"], s["POFD"], s["FAR"], s["FBI"], s["CSI"], s["AUC"])) return rows def combine(point, gate, mode, g): """Turn the point probability and the gate probability into one score.""" if mode == "point": return point if mode == "prod": return point * gate if mode == "hard": # step 1 of the proposal, literally return np.where(gate >= g, point, 0.0) raise ValueError(mode) def main(): ap = argparse.ArgumentParser() ap.add_argument("--tag", default="_ens") ap.add_argument("--gate", default=os.path.join(RES_DIR, "gate_ens.npz")) ap.add_argument("--gates", type=float, nargs="+", default=[0.05, 0.10, 0.20, 0.30, 0.50]) ap.add_argument("--spec-n-tag", default="", help="tag of the nEHD expert; with it the cascade becomes the " "complete mixture P(g)*P(.|yEHD) + (1-P(g))*P(.|nEHD)") ap.add_argument("--spec-tag", default="", help="tag of the step-2 specialist, fitted on the yEHD region-days " "only; when given, the proper factorisation " "P(gate) x P(point | yEHD) is scored as well") ap.add_argument("--out", default=os.path.join(RES_DIR, "cascade_ens.json")) args = ap.parse_args() G = np.load(args.gate) cover = float(G["cover"]) tr = load(args.tag, "oofpred_*.npz") va = load(args.tag, "pred_*.npz") print("[cascade] %d hàng điểm huấn luyện, %d kiểm định; cổng dùng ngưỡng phủ %.0f%%" % (tr["score"].size, va["score"].size, 100 * cover), flush=True) gt, ok_t = join_gate(tr, gate_key(G["tr_init"], G["tr_lead"], G["tr_region"]), G["tr_gate"]) gv, ok_v = join_gate(va, gate_key(G["va_init"], G["va_lead"], G["va_region"]), G["va_gate"]) ok_t &= np.isfinite(gt) ok_v &= np.isfinite(gv) tr = {k: v[ok_t] for k, v in tr.items()}; gt = gt[ok_t] va = {k: v[ok_v] for k, v in va.items()}; gv = gv[ok_v] print("[cascade] ghép được cổng cho %.1f%% hàng huấn luyện, %.1f%% hàng kiểm định" % (100 * ok_t.mean(), 100 * ok_v.mean()), flush=True) et, ev = tr["evt"] > 0.5, va["evt"] > 0.5 out = {"cover": cover, "results": {}, "gate_cost": {}} # What the gate throws away, before any model is applied to what is left. print("\n[chi phí của cổng] trên các năm kiểm định, toàn bộ hàng:") print(" %-8s %10s %10s %10s" % ("ngưỡng", "% hàng bỏ", "% sự kiện mất", "base rate còn lại")) for g in [0.0] + list(args.gates): closed = gv < g cost = {"rows_dropped": float(closed.mean()), "events_lost": float(ev[closed].sum()) / max(float(ev.sum()), 1.0), "base_rate_kept": float(ev[~closed].mean()) if (~closed).sum() else 0.0} out["gate_cost"]["%.2f" % g] = cost print(" %-8.2f %9.1f%% %9.1f%% %13.4f" % (g, 100 * cost["rows_dropped"], 100 * cost["events_lost"], cost["base_rate_kept"])) sp_t = sp_v = None if args.spec_tag: st_ = load(args.spec_tag, "oofpred_*.npz") sv_ = load(args.spec_tag, "pred_*.npz") # the specialist should write the same rows as the base model, but join on # (init, lead, region, point) rather than assume the order matches kt = gate_key(st_["init"], st_["lt"], st_["region"]) * 1000 + st_["point"] kb = gate_key(tr["init"], tr["lt"], tr["region"]) * 1000 + tr["point"] o = np.argsort(kt, kind="stable") pos = np.clip(np.searchsorted(kt[o], kb), 0, kt.size - 1) sp_t = np.where(kt[o][pos] == kb, st_["score"][o][pos], np.nan) kt = gate_key(sv_["init"], sv_["lt"], sv_["region"]) * 1000 + sv_["point"] kb = gate_key(va["init"], va["lt"], va["region"]) * 1000 + va["point"] o = np.argsort(kt, kind="stable") pos = np.clip(np.searchsorted(kt[o], kb), 0, kt.size - 1) sp_v = np.where(kt[o][pos] == kb, sv_["score"][o][pos], np.nan) print("[cascade] ghép được mô hình chuyên biệt cho %.1f%% hàng huấn luyện, " "%.1f%% hàng kiểm định" % (100 * np.isfinite(sp_t).mean(), 100 * np.isfinite(sp_v).mean()), flush=True) bad = ~np.isfinite(sp_t) sp_t[bad] = tr["score"][bad] bad = ~np.isfinite(sp_v) sp_v[bad] = va["score"][bad] sn_t = sn_v = None if args.spec_n_tag: sn_t = join_rows(load(args.spec_n_tag, "oofpred_*.npz"), tr) sn_v = join_rows(load(args.spec_n_tag, "pred_*.npz"), va) print("[cascade] ghép được chuyên gia nEHD cho %.1f%% / %.1f%% hàng" % (100 * np.isfinite(sn_t).mean(), 100 * np.isfinite(sn_v).mean()), flush=True) sn_t[~np.isfinite(sn_t)] = 0.0 sn_v[~np.isfinite(sn_v)] = 0.0 for frame, mt, mv in (("daily_all", np.ones(et.size, bool), np.ones(ev.size, bool)), ("daily_capable", tr["capable"] > 0.5, va["capable"] > 0.5)): print("\n" + "=" * 78 + "\n KHUNG %s (%d hàng kiểm định, base rate %.4f)" % (frame, int(mv.sum()), float(ev[mv].mean()))) wt, wv = tr["week"][mt], va["week"][mv] yt, yv = et[mt], ev[mv] specs = [("point", 0.0, "chỉ mô hình điểm (mốc so sánh hiện tại)"), ("prod", 0.0, "nhân xác suất cổng với xác suất điểm")] specs += [("hard", g, "cổng cứng: bỏ dự báo khi P(EHDE vùng) < %.2f" % g) for g in args.gates] for mode, g, lbl in specs: st = combine(tr["score"][mt], gt[mt], mode, g) sv = combine(va["score"][mv], gv[mv], mode, g) th = weekly_fbi1(st, yt, wt) key = "%s_%s%s" % (frame, mode, ("_%02d" % int(100 * g)) if mode == "hard" else "") out["results"][key] = table(apply_weekly(sv, wv, th), yv, wv, sv, lbl) if sp_t is None: continue pairs = [ ("spec_alone", sp_t[mt], sp_v[mv], "CHỈ mô hình chuyên biệt (huấn luyện trên yEHD), áp cho mọi hàng"), ("gate_x_spec", gt[mt] * sp_t[mt], gv[mv] * sp_v[mv], "EXP2 hai bước: P(EHDE vùng) x P(EHD ô | yEHD)")] if sn_t is not None: pairs.append(( "mixture", gt[mt] * sp_t[mt] + (1 - gt[mt]) * sn_t[mt], gv[mv] * sp_v[mv] + (1 - gv[mv]) * sn_v[mv], "EXP2 ĐẦY ĐỦ (hỗn hợp): P(g)xP(ô|yEHD) + (1-P(g))xP(ô|nEHD)")) for key, st, sv, lbl in (pairs): th = weekly_fbi1(st, yt, wt) out["results"]["%s_%s" % (frame, key)] = table( apply_weekly(sv, wv, th), yv, wv, sv, lbl) if sp_t is not None: # Diagnostic, NOT an operational score: restrict to the region-days that # truly were yEHD. Conditioning on truth is legitimate for asking whether # the specialist is better at its own job; conditioning on the gate would # not be. yes_t = np.isin(gate_key(tr["init"], tr["lt"], tr["region"]), gate_key(G["tr_init"], G["tr_lead"], G["tr_region"])[G["tr_Pobs"] > cover]) yes_v = np.isin(gate_key(va["init"], va["lt"], va["region"]), gate_key(G["va_init"], G["va_lead"], G["va_region"])[G["va_Pobs"] > cover]) print("\n" + "=" * 78 + "\n CHẨN ĐOÁN — chỉ các region-day THỰC SỰ là yEHD (%d hàng kiểm định, " "base rate %.4f). Không phải điểm vận hành." % (int(yes_v.sum()), float(ev[yes_v].mean()))) for key, st, sv, lbl in ( ("cond_base", tr["score"][yes_t], va["score"][yes_v], "mô hình điểm gốc (không chuyên biệt)"), ("cond_spec", sp_t[yes_t], sp_v[yes_v], "mô hình chuyên biệt huấn luyện trên đúng nhóm này")): th = weekly_fbi1(st, et[yes_t], tr["week"][yes_t]) out["results"][key] = table( apply_weekly(sv, va["week"][yes_v], th), ev[yes_v], va["week"][yes_v], sv, lbl) with open(args.out, "w") as fh: json.dump(out, fh, indent=1, default=float) print("\n[done] -> %s" % args.out, flush=True) if __name__ == "__main__": main()