#!/usr/bin/env python3 """Operational pack for one CFS init, in the SAME product set as the EC bulletin. The certified CFS chain (rt_forecast.py) issues one number per week per point: the downscaled weekly-max Tx and its margin over the local EHD threshold. The EC bulletin (../HW/ehd_forecast_ope1b.html) issues four things per point: accumulated severity, peak severity, a probability, and a binary field. This script produces the missing two WITHOUT re-fitting anything on the forecast itself: sev_acc sum over the days of the week of max(Sev, 0), Sev = Tx - threshold. Needs DAILY downscaled Tx, so the quantile map of the certified chain is refitted at daily resolution: same +-30 day window, same training years, same coarse cells, but the pairs are (daily model Tx at a lead inside the window, daily observed Tx of the valid day) instead of the weekly maxima. Nothing else changes. prob P(EHD in the week | margin), a logistic fit of the observed event on the margin of the certified chain, per week, TRAINING YEARS ONLY, pooled over points. This is where the CFS bulletin differs in kind from the EC one: EC has 72 AI members and reads p as their vote share M/N, CFS has one deterministic chain, so the probability is a calibration of its margin rather than an ensemble frequency. Both are probabilities of the same event and both are fitted out of sample; they are not the same estimator and the page says so. The weekly margin and the binary field are READ from the certified bulletin json, never recomputed, so the numbers on the page are the numbers that were certified. Also scored on the 726 validation inits, for the reference-skill and plausibility sections: coverage of the observed and forecast fields, POD/POFD/FAR/FBI/CSI/SEDI at the operational offset, and AUC of the margin. """ from __future__ import annotations import argparse import datetime as dt import json import os import sys import warnings import numpy as np warnings.filterwarnings("ignore") sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import cfs_common as C import build_obs as B import build_daily as D from bcsd import doy_window_mask, qmap_fit, qmap_apply from rt_forecast import read_rt, DELTA, REG_FRAC ITMAX = C.FLUX_NAMES.index("tmax") TXPTS = os.path.join(C.OUT, "Cube", "tx_daily_pts.f4") # (ninit, 46, 541) FIELDS = os.path.join(C.OUT, "AI_Results", "bcsd_fluxonly.json.w%d.fields.npz") # ------------------------------------------------------------------ daily cache def point_index(): z = np.load(B.TX_NPZ) g = C.flux_grid() olat, olon = z["lat"][z["jj"]], z["lon"][z["ii"]] pj = np.abs(olat[:, None] - g["wlat"][None, :]).argmin(1) pi = np.abs(olon[:, None] - g["wlon"][None, :]).argmin(1) _, cell = np.unique(pj * 1000 + pi, return_inverse=True) return pj, pi, cell def build_tx_daily_pts(pj, pi): """Daily model Tx at the 541 points, pulled once out of the 7.8 GB cube.""" if os.path.exists(TXPTS): return n, nlt, nv, ny, nx = D.shapes() src = np.memmap(D.DAILY, "float32", "r", shape=(n, nlt, nv, ny, nx)) dst = np.memmap(TXPTS, "float32", "w+", shape=(n, nlt, len(pj))) # Read each chunk CONTIGUOUSLY and slice in memory. Taking the strided view # src[:, :, ITMAX] straight off the memmap touches nearly every page of the # 7.8 GB cube anyway (the contiguous run per variable is 3456 B, smaller than # a page) but defeats Lustre readahead: measured >5 min and still unfinished, # against ~1 min for the sequential read below. for i0 in range(0, n, 64): i1 = min(i0 + 64, n) blk = np.asarray(src[i0:i1]) # sequential, ~163 MB dst[i0:i1] = blk[:, :, ITMAX, :, :][:, :, pj, pi] - 273.15 del blk dst.flush() print("[cache] %s" % TXPTS, flush=True) def obs_daily_at(inits, nlt): """Observed daily Tx on the valid day of every (init, lead).""" z = np.load(B.TX_NPZ) tx, oord = z["tx"], z["ordinal"] pos = {int(o): k for k, o in enumerate(oord)} out = np.full((len(inits), nlt, tx.shape[1]), np.nan, "f4") for i, s in enumerate(inits): s = str(s) o0 = dt.date(int(s[:4]), int(s[4:6]), int(s[6:8])).toordinal() for k in range(nlt): r = pos.get(o0 + k) if r is not None: out[i, k] = tx[r] return out # ------------------------------------------------------------------ scores def contingency(pred, obs): a = float(( pred & obs).sum()); b = float(( pred & ~obs).sum()) c = float((~pred & obs).sum()); d = float((~pred & ~obs).sum()) pod = a / max(a + c, 1); pofd = b / max(b + d, 1) far = b / max(a + b, 1); fbi = (a + b) / max(a + c, 1) csi = a / max(a + b + c, 1) # SEDI. The denominator is a SUM OF LOGS OF NUMBERS BELOW 1, i.e. negative; # clamping it away from zero with max(.., eps) turns it into +1e-12 and the # ratio explodes to ~1e12. Clamp the rates instead, never the denominator. h = min(max(pod, 1e-6), 1 - 1e-6) f = min(max(pofd, 1e-6), 1 - 1e-6) num = np.log(f) - np.log(h) - np.log(1 - f) + np.log(1 - h) den = np.log(f) + np.log(h) + np.log(1 - f) + np.log(1 - h) sedi = float(num / den) if den != 0 else float("nan") return dict(POD=pod, POFD=pofd, FAR=far, FBI=fbi, CSI=csi, SEDI=sedi) def region_events(field, reg, frac): """EXP1: a region has EHDE when more than `frac` of its points do.""" codes = np.unique(reg) return np.stack([field[:, reg == c].mean(1) > frac for c in codes], 1) def auc_score(score, obs): o = np.asarray(obs, bool); s = np.asarray(score, float) ok = np.isfinite(s) o, s = o[ok], s[ok] if o.all() or not o.any(): return float("nan") r = s.argsort().argsort().astype(float) + 1 n1 = float(o.sum()); n0 = float((~o).sum()) return float((r[o].sum() - n1 * (n1 + 1) / 2) / (n1 * n0)) def fit_prob(margin, event): """Logistic P(event | margin), fitted on the training rows handed in.""" from sklearn.linear_model import LogisticRegression ok = np.isfinite(margin) & np.isfinite(event.astype(float)) m = LogisticRegression(max_iter=1000) m.fit(margin[ok].reshape(-1, 1), event[ok].astype(int)) return m def main(): ap = argparse.ArgumentParser() ap.add_argument("--init", default="2026082900") ap.add_argument("--indir", default=None) ap.add_argument("--out", default=None) a = ap.parse_args() indir = a.indir or os.path.join(C.OUT, "RT", "cfs_%s" % a.init) out = a.out or os.path.join(C.OUT, "RT", "ope_%s.npz" % a.init) cert = json.load(open(os.path.join(C.OUT, "RT", "fcst_%s.json" % a.init))) pj, pi, cell = point_index() ncell = cell.max() + 1 npt = len(pj) build_tx_daily_pts(pj, pi) zt = np.load(B.TX_NPZ) reg = zt["reg"] thr_all = np.load(B.thr_path("train"))["thr"] # (365, npt) wk = np.load(os.path.join(C.OUT, "Cube", "weekly.npz")) inits, yrs, doy_w = wk["inits"], wk["yrs"], wk["doy"] train = np.isin(yrs, C.TRAIN_YEARS) valid = np.isin(yrs, C.VALID_YEARS) n = len(inits) Md = np.memmap(TXPTS, "float32", "r", shape=(n, C.NLT, npt)) Od = obs_daily_at(inits, C.NLT) # ------------------------------------------------- real-time daily model Tx step = read_rt(indir, a.init) # (184, ny, nx) K daily = np.stack([np.nanmax(step[4 * k:4 * k + 4], 0) for k in range(C.NLT)]) - 273.15 rt_daily = daily[:, pj, pi] # (46, npt) d0 = dt.date(int(a.init[:4]), int(a.init[4:6]), int(a.init[6:8])) sev_acc, sev_max, prob, pred, margin = {}, {}, {}, {}, {} vprob, vpred, vevent = {}, {}, {} skill, skill_reg, rel = {}, {}, {} for w, (wname, idx) in enumerate(C.WEEKS): cw = cert["weeks"][wname] margin[wname] = np.array(cw["margin"], "f4") pred[wname] = np.array(cw["event"], bool) sev_max[wname] = np.maximum(margin[wname], 0.0) q = cw["doy"] # ---- daily BCSD, same window and years as the certified weekly map ---- sel = np.nonzero(doy_window_mask(doy_w[:, w], q) & train)[0] md = np.asarray(Md[sel][:, idx, :]).reshape(-1, npt) # (nsel*nlead, npt) od = Od[sel][:, idx, :].reshape(-1, npt) ok = np.isfinite(md).all(1) & np.isfinite(od).all(1) md, od = md[ok], od[ok] mc = np.stack([md[:, cell == c].mean(1) for c in range(ncell)], 1) oc = np.stack([od[:, cell == c].mean(1) for c in range(ncell)], 1) rc = np.stack([rt_daily[:, cell == c][idx].mean(1) for c in range(ncell)], 1) fine_clim = np.nanmean(od, axis=0) # (npt) d1d = np.empty((len(idx), npt), "f4") for c in range(ncell): mq, oq = qmap_fit(mc[:, c], oc[:, c]) anom = qmap_apply(rc[:, c], mq, oq) - oc[:, c].mean() d1d[:, cell == c] = (fine_clim[cell == c][None, :] + anom[:, None]).astype("f4") days = [d0 + dt.timedelta(days=int(k)) for k in idx] thr_d = np.stack([thr_all[dd.timetuple().tm_yday - 1] for dd in days]) sev_acc[wname] = np.maximum(d1d - thr_d, 0.0).sum(0) # ---- probability, fitted on training years of the reforecast ---------- f = np.load(FIELDS % w) mg, ev, fy = f["d1"] - f["thr"], f["E"].astype(bool), f["yrs"] tr = np.isin(fy, C.TRAIN_YEARS); va = np.isin(fy, C.VALID_YEARS) mdl = fit_prob(mg[tr].ravel(), ev[tr].ravel()) prob[wname] = mdl.predict_proba(margin[wname].reshape(-1, 1))[:, 1] pv = mdl.predict_proba(mg[va].ravel().reshape(-1, 1))[:, 1] vprob[wname] = pv.reshape(mg[va].shape) vpred[wname] = (mg[va] >= DELTA) vevent[wname] = ev[va] s = contingency(vpred[wname].ravel(), vevent[wname].ravel()) s["AUC"] = auc_score(mg[va].ravel(), vevent[wname].ravel()) # EXP1 unit: the certified product is a REGIONAL yes/no, so score it there # too. The EC bulletin's table is per grid point, so the page carries both # and says which is which; the two units have different base rates (0.46 # against 0.65) and are not comparable to each other. rp = region_events(vpred[wname], reg, REG_FRAC) ro = region_events(vevent[wname], reg, REG_FRAC) rs = contingency(rp.ravel(), ro.ravel()) rs["AUC"] = auc_score( np.stack([mg[va][:, reg == c].mean(1) for c in np.unique(reg)], 1).ravel(), ro.ravel()) rs["base"] = float(ro.mean()) skill_reg[wname] = rs s["BRIER"] = float(np.mean((pv - vevent[wname].ravel().astype(float)) ** 2)) s["base"] = float(vevent[wname].mean()) skill[wname] = s # reliability, 10 bins, validation years bins = np.clip((pv * 10).astype(int), 0, 9) rel[wname] = [[int((bins == k).sum()), float(pv[bins == k].mean()) if (bins == k).any() else np.nan, float(vevent[wname].ravel()[bins == k].mean()) if (bins == k).any() else np.nan] for k in range(10)] print("%s sev_acc %.2f P %.3f | ô lưới POD %.3f FAR %.3f CSI %.3f AUC %.3f" " | vùng POD %.3f POFD %.3f FAR %.3f CSI %.3f" % (wname, sev_acc[wname].mean(), prob[wname].mean(), s["POD"], s["FAR"], s["CSI"], s["AUC"], rs["POD"], rs["POFD"], rs["FAR"], rs["CSI"]), flush=True) np.savez(out, reg=reg, delta=DELTA, reg_frac=REG_FRAC, init=a.init, skill=json.dumps(skill), skill_reg=json.dumps(skill_reg), rel=json.dumps(rel), **{"sev_acc_%s" % w: sev_acc[w] for w in C.WEEK_NAMES}, **{"sev_max_%s" % w: sev_max[w] for w in C.WEEK_NAMES}, **{"prob_%s" % w: prob[w] for w in C.WEEK_NAMES}, **{"pred_%s" % w: pred[w] for w in C.WEEK_NAMES}, **{"margin_%s" % w: margin[w] for w in C.WEEK_NAMES}, **{"cov_mod_%s" % w: vpred[w].mean(1) for w in C.WEEK_NAMES}, **{"cov_obs_%s" % w: vevent[w].mean(1) for w in C.WEEK_NAMES}, **{"cov_prob_%s" % w: (vprob[w] > 0.5).mean(1) for w in C.WEEK_NAMES}) print("\nsaved", out) if __name__ == "__main__": main()