#!/usr/bin/env python3 """Test one thing: does training a probabilistic classifier on the EVENT beat thresholding a downscaled deterministic Tx, given exactly the same information? The deterministic route (downscale Tx, then flag Tx >= Tx_thres) reaches only CSI ~0.14 at grid points and ~0.20 for regional EHDE. This script keeps the same inputs -- the BCSD-downscaled Tx and the fields derived from it -- but * predicts P(EHD) with a binary classifier instead of thresholding, and * tunes the decision probability per (region, forecast week). Both the classifier and the decision threshold are fitted on the training years of each fold only; the threshold search uses an inner split of those training years, so no test year influences either the model or its operating point. Usage: ehd_classifier.py [--pred outputs/downscale/pred/pred_bcsd.npy] [--folds 5] """ import argparse import json import os import sys import numpy as np import xarray as xr import lightgbm as lgb ROOT = "/home/tanpv/Progs/claude/HW_EC_10" OBS = f"{ROOT}/outputs/OBS_4EC15/tx_obs_aligned.nc" STATIC = f"{ROOT}/outputs/Cli/static_fields_0.25.nc" CLI = f"{ROOT}/outputs/Cli" WEEKS = [("W1", 0, 11), ("W2", 11, 18), ("W3", 18, 25), ("W4", 25, 32), ("W5", 32, 39), ("W6", 39, 46)] FOLD_START = [2003, 2007, 2011, 2015, 2019] NDOY = 365 TARGET_POD = 0.60 # user-set operating point, see CLAUDE.md TARGET_FAR = 0.20 sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) from evaluate_ehd import contingency # noqa: E402 def main(): ap = argparse.ArgumentParser() ap.add_argument("--pred", default=f"{ROOT}/outputs/downscale/pred/pred_bcsd.npy") ap.add_argument("--folds", type=int, default=5) ap.add_argument("--subsample", type=int, default=6_000_000) ap.add_argument("--cover", type=float, default=0.20) ap.add_argument("--out", default=f"{ROOT}/outputs/downscale/ehd_classifier.json") ap.add_argument("--seed", type=int, default=0) args = ap.parse_args() rng = np.random.default_rng(args.seed) o = xr.open_dataset(OBS) st = xr.open_dataset(STATIC) land = (np.isfinite(st.vn.values) & (st.vn.values > 0)).ravel() land_idx = np.where(land)[0] nland = len(land_idx) reg = st.reg.values.ravel()[land_idx] ny, ni, nl = o.sizes["year"], o.sizes["init"], o.sizes["lead"] obs = o.tx_obs.values.reshape(ny, ni, nl, -1)[..., land_idx] doy = o.doy.values txp = np.load(args.pred) print(f"pred {txp.shape}, {nland} land points", flush=True) thr = np.empty_like(obs) for k, y0 in enumerate(FOLD_START): f = np.arange(k * 4, k * 4 + 4) t = xr.open_dataset(f"{CLI}/Qp_tx_obs_fold_{y0}.nc" ).tx_thres_obs.values.reshape(NDOY, -1)[:, land_idx] thr[f] = t[doy[f] - 1] ehd_obs = obs >= thr sev_pred = txp - thr # the deterministic decision variable wk_of_lead = np.zeros(nl, dtype=int) for w, (_, l0, l1) in enumerate(WEEKS): wk_of_lead[l0:l1] = w stat = np.stack([st.dz.values.ravel(), st.orog.values.ravel(), st.orog_std.values.ravel(), st.lat2d.values.ravel(), st.lon2d.values.ravel()], axis=1)[land_idx].astype("float32") FEAT = ["sev_pred", "tx_pred", "thr", "dz", "orog", "orog_std", "lat", "lon", "lead", "doy_sin", "doy_cos", "reg"] def build(years, rows): sv = sev_pred[years].reshape(-1, nland)[rows] tp = txp[years].reshape(-1, nland)[rows] th = thr[years].reshape(-1, nland)[rows] dd = doy[years].reshape(-1)[rows] ll = np.tile(np.arange(nl), len(years) * ni)[rows] n = len(rows) * nland X = np.empty((n, len(FEAT)), dtype="float32") X[:, 0] = sv.ravel(); X[:, 1] = tp.ravel(); X[:, 2] = th.ravel() for k in range(3): X[:, 3 + k] = np.tile(stat[:, k], len(rows)) X[:, 6] = np.tile(stat[:, 3], len(rows)) X[:, 7] = np.tile(stat[:, 4], len(rows)) X[:, 8] = np.repeat(ll / 45.0, nland) X[:, 9] = np.repeat(np.sin(2 * np.pi * dd / NDOY), nland) X[:, 10] = np.repeat(np.cos(2 * np.pi * dd / NDOY), nland) X[:, 11] = np.tile(reg, len(rows)) y = ehd_obs[years].reshape(-1, nland)[rows].ravel() return X, y cuts, cuts_maxcsi = {}, {} per = ny // args.folds # Two separate arrays on purpose. A year that is a TEST year of one fold can be an # inner-validation year of another; sharing one array would let a later fold # overwrite an earlier fold's held-out prediction with a different model's output. prob = np.full((ny, ni, nl, nland), np.nan, dtype="float32") # test years only prob_val = np.full((ny, ni, nl, nland), np.nan, dtype="float32") # tuning only for fi in range(args.folds): te = np.arange(fi * per, (fi + 1) * per) tr_all = np.array([i for i in range(ny) if i not in te]) # inner split of the training years: 3/4 fit, 1/4 to pick the thresholds inner_val = tr_all[rng.permutation(len(tr_all))[:max(1, len(tr_all) // 4)]] inner_fit = np.array([i for i in tr_all if i not in inner_val]) ntime = len(inner_fit) * ni * nl want = max(1, args.subsample // nland) rows = np.sort(rng.choice(ntime, min(want, ntime), replace=False)) Xtr, ytr = build(inner_fit, rows) ok = np.isfinite(Xtr).all(axis=1) & np.isfinite(ytr.astype("float32")) print(f"fold {fi+1}: fit on {ok.sum():,} rows, base rate {ytr[ok].mean():.3f}", flush=True) clf = lgb.LGBMClassifier(n_estimators=500, learning_rate=0.05, num_leaves=127, min_child_samples=200, subsample=0.8, subsample_freq=1, colsample_bytree=0.8, reg_lambda=1.0, n_jobs=16, verbose=-1) clf.fit(Xtr[ok], ytr[ok], feature_name=FEAT) del Xtr, ytr for years, dest in ((inner_val, prob_val), (te, prob)): nt = len(years) * ni * nl out = np.empty((nt, nland), dtype="float32") for a in range(0, nt, 4000): r = np.arange(a, min(a + 4000, nt)) X, _ = build(years, r) m = np.isfinite(X).all(axis=1) p = np.full(len(X), np.nan, dtype="float32") if m.any(): p[m] = clf.predict_proba(X[m])[:, 1] out[a:a + len(r)] = p.reshape(len(r), nland) dest[years] = out.reshape(len(years), ni, nl, nland) print(f"fold {fi+1}/{args.folds} scored", flush=True) # decision thresholds per (region, week) from the inner validation years only for w, (wn, l0, l1) in enumerate(WEEKS): for r in range(1, 8): m = reg == r if m.sum() == 0: continue pv = prob_val[inner_val][:, :, l0:l1][..., m].ravel() ov = ehd_obs[inner_val][:, :, l0:l1][..., m].ravel() good = np.isfinite(pv) pv, ov = pv[good], ov[good] # Operating point per user decision (2026-09-03): satisfy # POD >= 0.6 and FAR <= 0.2 first, maximising CSI among the feasible # cuts; if no cut is feasible, fall back to the best CSI overall. best_c, best_csi = 0.5, -1.0 feas_c, feas_csi = None, -1.0 for cc in np.arange(0.02, 0.99, 0.01): hit = ((pv >= cc) & ov).sum() miss = ((pv < cc) & ov).sum() fa = ((pv >= cc) & ~ov).sum() csi = hit / max(hit + miss + fa, 1) pod = hit / max(hit + miss, 1) far = fa / max(hit + fa, 1) if csi > best_csi: best_csi, best_c = csi, cc if pod >= TARGET_POD and far <= TARGET_FAR and csi > feas_csi: feas_csi, feas_c = csi, cc cuts[(fi, wn, r)] = feas_c if feas_c is not None else best_c cuts_maxcsi[(fi, wn, r)] = best_c # apply the tuned cut to the held-out years ehd_pred = np.zeros_like(ehd_obs) for fi in range(args.folds): te = np.arange(fi * per, (fi + 1) * per) for wn, l0, l1 in WEEKS: for r in range(1, 8): m = reg == r if m.sum() == 0: continue c = cuts[(fi, wn, r)] blk = prob[te][:, :, l0:l1][..., m] >= c idx = np.where(m)[0] ehd_pred[np.ix_(te, np.arange(ni), np.arange(l0, l1), idx)] = blk valid = np.isfinite(prob) & np.isfinite(obs) res = {"gridpoint": {}, "regional_EHDE": {}} for wn, l0, l1 in WEEKS + [("ALL", 0, nl)]: v = valid[:, :, l0:l1] res["gridpoint"][wn] = contingency(ehd_pred[:, :, l0:l1][v], ehd_obs[:, :, l0:l1][v]) po, pp = [], [] for r in range(1, 8): m = reg == r if m.sum() == 0: continue po.append(ehd_obs[:, :, l0:l1][..., m].mean(axis=-1).ravel()) pp.append(ehd_pred[:, :, l0:l1][..., m].mean(axis=-1).ravel()) res["regional_EHDE"][wn] = contingency(np.concatenate(pp) > args.cover, np.concatenate(po) > args.cover) for lvl in ("gridpoint", "regional_EHDE"): v = res[lvl]["ALL"] print(f"{lvl:15s} base={v['base']:.3f} POD={v['POD']:.3f} POFD={v['POFD']:.3f} " f"FAR={v['FAR']:.3f} FBI={v['FBI']:.3f} CSI={v['CSI']:.3f}", flush=True) os.makedirs(os.path.dirname(args.out), exist_ok=True) with open(args.out, "w") as fh: json.dump({"classifier": res, "cuts": {f"{k[0]}|{k[1]}|S{k[2]}": float(v) for k, v in cuts.items()}}, fh, indent=2) np.save(f"{ROOT}/outputs/downscale/pred/prob_ehd_bcsd.npy", prob) np.save(f"{ROOT}/outputs/downscale/pred/prob_ehd_bcsd_innerval.npy", prob_val) np.save(f"{ROOT}/outputs/downscale/pred/ehd_obs.npy", ehd_obs) print(f"wrote {args.out}") return 0 if __name__ == "__main__": sys.exit(main())