#!/usr/bin/env python3 """EXP3 gộp theo LEAD WEEK: N mô hình AI cho từng (ô lưới, tuần), quyết định hai tầng bằng p = M/N. VÌ SAO PHẢI GỘP TUẦN. Bản exp3_multi.py trước đó đưa cả 46 lead day của một ô lưới vào CÙNG một mô hình và CÙNG một bảng contingency. Hệ quả: (1) "dung lượng mẫu" hiện ra là 19136 = 416 init x 46 lead day, tức mỗi init bị đếm 46 lần chứ không phải 19136 trường hợp độc lập; (2) mọi khác biệt giữa W1 và W6 bị trộn phẳng nên bảng đánh giá không còn ý nghĩa theo hạn dự báo. Ở đây một MẪU = một INITIALISATION, nên mỗi (ô lưới, tuần) có đúng 2080 mẫu cho 20 năm, tách thành 1664 mẫu huấn luyện (16 năm) và 416 mẫu kiểm định (4 năm). ĐỊNH NGHĨA yEHDE THEO TUẦN. Một tuần lead tại một ô lưới là yEHDE nếu có ÍT NHẤT MỘT ngày EHD trong tuần đó. Tương đương Sev_week = max(Sev_obs trong tuần) >= 0, nên target hồi quy vẫn là Sev_obs giữ liên tục và có dấu đúng như đặc tả EXP3. Base rate đo được: W1 0.430 (11 ngày), W2-W6 ~0.329 (7 ngày) -- so với 0.097 khi tính theo từng ngày. Tức tỷ lệ nEHDE là ~67%, không phải ~90%. MỖI (ĐIỂM, TUẦN) MỘT BỘ MÔ HÌNH RIÊNG -- 21 x 6 = 126 bộ, mỗi bộ N = 7 họ x 9 mức phân vị = 63 thành viên. Không gộp 6 tuần vào một mô hình chung. """ from __future__ import annotations import argparse import json import os import sys import time import numpy as np from joblib import Parallel, delayed sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import ehd_common as C import ehd_metrics as M import ai_data as D import build_ai_inputs as A import build_cube as B from exp3_multi import model_zoo, ZOO, NQ, ensemble_votes, oof_votes RES_DIR = os.path.join(C.OUT, "AI_Results") OUT_JSON = os.path.join(RES_DIR, "exp3_week.json") OUT_NPZ = os.path.join(RES_DIR, "exp3_week_p.npz") # Trường mà giá trị CỰC ĐẠI trong tuần mang thông tin khác hẳn giá trị trung bình: # một đợt nóng ba ngày kẹp giữa bốn ngày mát có trung bình tuần rất tầm thường. MAX_VARS = ["sv_mod", "ev_mod", "z_mx2t6", "z_t2m", "z_skt", "z_d2m", "z_thk_1000_500", "z_thk_925_700", "z_dt_850_500", "z_cape"] MIN_KEEP = 150 # tối thiểu số mẫu còn lại sau cổng tầng 1 MIN_KEEP_EVT = 20 # tối thiểu số mẫu dương trong tập đó # ------------------------------------------------------- gộp ngày -> tuần def weekly_frames(store, points, years, fset): """Ma trận đặc trưng và target ở mức (init, tuần, điểm). Trả về dict theo chỉ số tuần: (X, sev, evt, year, point, names). """ isel = np.nonzero(np.isin(store.init_year, years))[0] X, y, info = store.build_isel(points, isel, fset) ni, nk, npt = isel.size, C.NLT, points.size nf = X.shape[1] names = store.feature_names(fset) X = X.reshape(ni, nk, npt, nf) sev = y["sev"].reshape(ni, nk, npt) evt = y["evt"].reshape(ni, nk, npt) yr = store.init_year[isel] imax = [names.index(v) for v in MAX_VARS if v in names] isv = names.index("sv_mod") out = {} for w, (wname, lts) in enumerate(C.WEEKS): Xw = X[:, lts, :, :] # trung bình tuần của mọi đặc trưng, cực đại tuần của nhóm nhạy nhiệt, và # số ngày mô hình tự cho là nóng -- chính là "đếm EHD tuần" của EC cols = [np.nanmean(Xw, axis=1), np.nanmax(Xw[:, :, :, imax], axis=1), (Xw[:, :, :, isv] >= 0.0).sum(axis=1, dtype="f4")[:, :, None]] Xa = np.concatenate(cols, axis=2).reshape(ni * npt, -1) nm = (["m_" + s for s in names] + ["x_" + names[i] for i in imax] + ["n_hotday_mod"]) sw = np.nanmax(sev[:, lts, :], axis=1).reshape(ni * npt) ew = (np.nansum(evt[:, lts, :], axis=1) >= 1.0).reshape(ni * npt) yy = np.broadcast_to(yr[:, None], (ni, npt)).reshape(ni * npt) pp = np.broadcast_to(points[None, :], (ni, npt)).reshape(ni * npt) out[w] = (Xa, sw, ew.astype("f4"), yy, pp, nm, wname) return out # --------------------------------------------------------- một (điểm, tuần) def run_cell(ip, wname, Xtr, sev_tr, evt_tr, wtr, Xte, p1_grid): """Cả hai tầng EXP3 tại một ô lưới cho MỘT tuần lead.""" sd = Xtr.std(axis=0) col = sd > 1e-12 Xtr = np.ascontiguousarray(Xtr[:, col]); Xte = np.ascontiguousarray(Xte[:, col]) zoo = model_zoo() t0 = time.time() q_oof, _ = oof_votes(Xtr, sev_tr, wtr, zoo) q_val, _, mu_val = ensemble_votes(Xtr, sev_tr, Xte, zoo) out = {"stage1": (q_oof, q_val), "mu_frac_pos": float((mu_val >= 0).mean()), "nfeat": int(col.sum()), "stage2": {}} for p1 in p1_grid: keep = q_oof >= p1 if keep.sum() < MIN_KEEP or evt_tr[keep].sum() < MIN_KEEP_EVT: continue o2 = np.zeros(Xtr.shape[0]) a, _ = oof_votes(Xtr[keep], sev_tr[keep], wtr[keep], zoo) o2[keep] = a b, _, _ = ensemble_votes(Xtr[keep], sev_tr[keep], Xte, zoo) out["stage2"]["%.3f" % p1] = (o2, b, int(keep.sum())) print("[%s p%3d] train %d (%d dương), %d feature, %d ngưỡng p1, %.0fs" % (wname, ip, Xtr.shape[0], int(evt_tr.sum()), int(col.sum()), len(out["stage2"]), time.time() - t0), flush=True) return ip, wname, out def fbi1_threshold(score, obs): """Ngưỡng làm tần suất dự báo bằng đúng tần suất quan trắc, tức FBI = 1. Chốt ngưỡng bằng max-CSI luôn kéo về phía dự báo thừa khi mô hình còn yếu: CSI = a/(a+b+c) đánh đổi một false alarm lấy một hit với giá ngang nhau, nên ở W3-W6 nó chạy thẳng ra biên thấp nhất của lưới p2 và cho FBI ~ 2.5. Hiệu chỉnh theo tần suất khí hậu không có bậc tự do nào để chạy ra biên, và FBI ~ 1 vốn là một trong các chỉ tiêu phải đạt. """ r = float(obs.mean()) if r <= 0.0 or r >= 1.0: return np.inf return float(np.quantile(score, 1.0 - r)) def gated_score(q1, p1, s2): """Điểm số cuối: cổng tầng 1 đóng thì loại hẳn, mở thì lấy p của tầng 2.""" return np.where(q1 >= p1, s2, -1.0) def contingency(pred, obs): a, b, c, d = M.contingency(pred, obs) s = M.scores(a, b, c, d) s.update({"hit": a, "false_alarm": b, "miss": c, "correct_neg": d, "PC": (a + d) / max(a + b + c + d, 1), "n": int(obs.size), "base_rate": float(obs.mean())}) return s def main(): ap = argparse.ArgumentParser() ap.add_argument("--per-region", type=int, default=3) ap.add_argument("--seed", type=int, default=0) ap.add_argument("--fset", default="F3") ap.add_argument("--p1-grid", type=float, nargs="*", default=[0.05, 0.10, 0.20, 0.30, 0.40]) ap.add_argument("--p2-grid", type=float, nargs="*", default=[0.20, 0.30, 0.40, 0.50, 0.60, 0.70]) ap.add_argument("--jobs", type=int, default=20) ap.add_argument("--out", default=OUT_JSON) args = ap.parse_args() store = D.Store() pts = D.sample_points(store.grid, args.per_region, args.seed) print("[exp3w] %d điểm x %d tuần = %d bộ mô hình, N = %d họ x %d phân vị = %d " "thành viên" % (pts.size, len(C.WEEKS), pts.size * len(C.WEEKS), len(ZOO), NQ, len(ZOO) * NQ), flush=True) tr = weekly_frames(store, pts, C.TRAIN_YEARS, args.fset) te = weekly_frames(store, pts, C.VALID_YEARS, args.fset) for w, (_, _, ew, _, _, nm, wname) in tr.items(): print("[data] %s train %d mẫu/điểm, valid %d mẫu/điểm, %d feature tuần, " "base rate %.4f" % (wname, tr[w][0].shape[0] // pts.size, te[w][0].shape[0] // pts.size, len(nm), ew.mean()), flush=True) jobs = [] for w, (Xa, sw, ew, yy, pp, nm, wname) in tr.items(): Xv, sv_, ev_, yv, pv, _, _ = te[w] for p in pts: mt, mv = pp == p, pv == p ok_t = np.isfinite(Xa[mt]).all(axis=1) & np.isfinite(sw[mt]) ok_v = np.isfinite(Xv[mv]).all(axis=1) & np.isfinite(ev_[mv]) jobs.append((int(p), wname, Xa[mt][ok_t], sw[mt][ok_t], ew[mt][ok_t], yy[mt][ok_t], Xv[mv][ok_v], args.p1_grid)) got = Parallel(n_jobs=args.jobs, verbose=5)(delayed(run_cell)(*j) for j in jobs) # quan trắc của tập kiểm định, tra theo (điểm, tuần) obs = {} for w, (Xv, sv_, ev_, yv, pv, _, wname) in te.items(): for p in pts: mv = pv == p ok_v = np.isfinite(Xv[mv]).all(axis=1) & np.isfinite(ev_[mv]) obs[(int(p), wname)] = ev_[mv][ok_v] > 0.5 obs_tr = {} for w, (Xa, sw, ew, yy, pp, _, wname) in tr.items(): for p in pts: mt = pp == p ok_t = np.isfinite(Xa[mt]).all(axis=1) & np.isfinite(sw[mt]) obs_tr[(int(p), wname)] = ew[mt][ok_t] > 0.5 res = {"config": {"points": pts.tolist(), "weeks": C.WEEK_NAMES, "n_family": len(ZOO), "n_quantile": NQ, "n_members": len(ZOO) * NQ, "fset": args.fset, "target": "Sev_week = max(Sev_obs trong tuần); yEHDE nếu >= 0", "p1_grid": args.p1_grid, "p2_grid": args.p2_grid}, "models": [{"name": n, "desc": d} for n, d in ZOO], "by_cell": {}, "by_week": {}, "sweep": {}, "best": {}} # Ngưỡng (p1, p2) chốt RIÊNG TỪNG TUẦN, và chốt trên trường NGOÀI MẪU của các # năm huấn luyện -- dò ngưỡng theo CSI của tập kiểm định chính là dò tham số # trên tập test. byw = {} for ip, wname, r in got: byw.setdefault(wname, []).append((ip, r)) dump = {} for wname, cells in byw.items(): for ip, r in cells: dump["%d|%s|q_oof" % (ip, wname)] = r["stage1"][0] dump["%d|%s|q_val" % (ip, wname)] = r["stage1"][1] dump["%d|%s|obs_oof" % (ip, wname)] = obs_tr[(ip, wname)] dump["%d|%s|obs_val" % (ip, wname)] = obs[(ip, wname)] for k, v in r["stage2"].items(): dump["%d|%s|s2_oof|%s" % (ip, wname, k)] = v[0] dump["%d|%s|s2_val|%s" % (ip, wname, k)] = v[1] np.savez_compressed(OUT_NPZ, **dump) print("[dump] mảng p thô -> %s" % OUT_NPZ, flush=True) for wname, cells in byw.items(): # --- quy tac A: max-CSI tren truong ngoai mau (de doi chieu) best = None for p1 in args.p1_grid: k = "%.3f" % p1 for p2 in args.p2_grid: po, oo, pv, ov = [], [], [], [] for ip, r in cells: if k not in r["stage2"]: continue po.append((r["stage1"][0] >= p1) & (r["stage2"][k][0] > p2)) oo.append(obs_tr[(ip, wname)]) pv.append((r["stage1"][1] >= p1) & (r["stage2"][k][1] > p2)) ov.append(obs[(ip, wname)]) if not po: continue so = contingency(np.concatenate(po), np.concatenate(oo)) sv = contingency(np.concatenate(pv), np.concatenate(ov)) res["sweep"]["%s|p1=%.2f,p2=%.2f" % (wname, p1, p2)] = { "oof": so, "valid": sv} if best is None or so["CSI"] > best[0]: best = (so["CSI"], p1, p2) # --- quy tac B: p1 chon theo CSI ngoai mau, p2 hieu chinh RIENG TUNG DIEM # de FBI ~ 1. Moi o luoi co tan suat khi hau rieng nen mot p2 chung # cho ca 21 diem se thua o diem nay va thieu o diem kia. bestB = None for p1 in args.p1_grid: k = "%.3f" % p1 pv, ov = [], [] for ip, r in cells: if k not in r["stage2"]: continue so_ = gated_score(r["stage1"][0], p1, r["stage2"][k][0]) th = fbi1_threshold(so_, obs_tr[(ip, wname)]) sv_ = gated_score(r["stage1"][1], p1, r["stage2"][k][1]) pv.append(sv_ > th); ov.append(obs[(ip, wname)]) if not pv: continue oo, oob = [], [] for ip, r in cells: if k not in r["stage2"]: continue so_ = gated_score(r["stage1"][0], p1, r["stage2"][k][0]) th = fbi1_threshold(so_, obs_tr[(ip, wname)]) oo.append(so_ > th); oob.append(obs_tr[(ip, wname)]) sO = contingency(np.concatenate(oo), np.concatenate(oob)) sV = contingency(np.concatenate(pv), np.concatenate(ov)) res["sweep"]["%s|FBI1,p1=%.2f" % (wname, p1)] = {"oof": sO, "valid": sV} if bestB is None or sO["CSI"] > bestB[0]: bestB = (sO["CSI"], p1) if best is not None: _, p1b, p2b = best res["best"][wname] = {"rule": "maxCSI", "p1": p1b, "p2": p2b, "CSI_oof": best[0]} if bestB is None: continue _, p1c = bestB k = "%.3f" % p1c res.setdefault("best_fbi1", {})[wname] = {"rule": "FBI1", "p1": p1c, "CSI_oof": bestB[0]} pv, ov = [], [] for ip, r in cells: if k not in r["stage2"]: continue so_ = gated_score(r["stage1"][0], p1c, r["stage2"][k][0]) th = fbi1_threshold(so_, obs_tr[(ip, wname)]) sv_ = gated_score(r["stage1"][1], p1c, r["stage2"][k][1]) pred = sv_ > th o = obs[(ip, wname)] s = contingency(pred, o) s.update({"region": int(store.grid.reg[ip]), "lat": float(store.st["lat"][ip]), "lon": float(store.st["lon"][ip]), "gate_closed_frac": float((r["stage1"][1] < p1c).mean()), "n_stage2_train": r["stage2"][k][2], "nfeat": r["nfeat"], "p1": p1c, "p2_fbi1": th}) res["by_cell"]["%d|%s" % (ip, wname)] = s pv.append(pred); ov.append(o) if pv: res["by_week"][wname] = contingency(np.concatenate(pv), np.concatenate(ov)) b = res["by_week"][wname] print("[%s] FBI1 p1=%.2f -> POD %.3f POFD %.3f FAR %.3f FBI %.3f " "CSI %.3f (n=%d)" % (wname, p1c, b["POD"], b["POFD"], b["FAR"], b["FBI"], b["CSI"], b["n"]), flush=True) with open(args.out, "w") as fh: json.dump(res, fh, indent=1, default=float) print("\n[done] -> %s" % args.out, flush=True) if __name__ == "__main__": main()