#!/usr/bin/env python3 """EXP3 tầng 1, phương án (b): thành viên NEO vào dự báo của chính EC. Ba biến thể được chấm trên cùng một tập kiểm định, cộng một mốc tham chiếu: p_raw -- tiên nghiệm của chính EC, hiệu chỉnh xác suất từ bốn thống kê tuần của dự báo tất định. Đây là MỐC THAM CHIẾU mà mọi mô hình AI phải vượt; các bảng trước đây của dự án chưa hề có mốc này. A -- KHÔNG neo: thành viên hồi quy thẳng Sev_obs (bản đang dùng ở mục 3f). B1 -- NEO BẰNG Sev: thành viên hồi quy PHẦN DƯ Sev_obs - Sev_mod, dự báo cuối là Sev_mod + hiệu chỉnh. Mô hình chỉ phải học sai số của EC. B2 -- NEO BẰNG XÁC SUẤT, offset thật: p = sigmoid(logit(p_raw) + net(X)), đúng cơ chế của Van Straaten et al. (2023). Công thức p = M/N giữ nguyên ở cả ba biến thể. Với thành viên hồi quy, thành viên (họ f, phân vị q) bỏ phiếu yEHDE khi Sev dự báo >= 0. Với thành viên phân loại có offset, thành viên (f, q) bỏ phiếu khi p_f >= tau_q, nên M/N vẫn là một con số trong [0, 1] và xấp xỉ trung bình xác suất của các họ. """ from __future__ import annotations import argparse import json import os import sys import time import numpy as np from joblib import Parallel, delayed from scipy.stats import norm 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 exp3_zoo as Z import exp3_offset as OFF from exp3_multi import conditional_spread, apply_spread from exp3_week import MAX_VARS from exp3_week_s1 import gate_scores, pick_p1, P1_GRID RES_DIR = os.path.join(C.OUT, "AI_Results") OUT_JSON = os.path.join(RES_DIR, "exp3_week_s1b.json") OUT_NPZ = os.path.join(RES_DIR, "exp3_week_s1b_votes.npz") NFOLD = 4 NQ = 9 QLEV = (np.arange(NQ) + 0.5) / NQ ZQ = norm.ppf(QLEV) def agg_names(store, fset): names = store.feature_names(fset) imax = [names.index(v) for v in MAX_VARS if v in names] return (["m_" + s for s in names] + ["x_" + names[i] for i in imax] + ["n_hotday_mod"]), imax def weekly_data(store, points, years, fset): 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, :, :] agg = np.concatenate( [np.nanmean(Xw, axis=1), np.nanmax(Xw[:, :, :, imax], axis=1), (Xw[:, :, :, isv] >= 0.0).sum(axis=1, dtype="f4")[:, :, None]], axis=2).reshape(ni * npt, -1) out[wname] = dict( agg=agg, sev=np.nanmax(sev[:, lts, :], axis=1).reshape(ni * npt), # mức nghiêm trọng tuần THEO MÔ HÌNH -- tiên nghiệm để neo biến thể B1 sev_mod=np.nanmax(Xw[:, :, :, isv], axis=1).reshape(ni * npt), evt=(np.nansum(evt[:, lts, :], axis=1) >= 1.0).reshape(ni * npt), year=np.broadcast_to(yr[:, None], (ni, npt)).reshape(ni * npt), point=np.broadcast_to(points[None, :], (ni, npt)).reshape(ni * npt)) return out def folds_of(year): ys = np.unique(year) return np.array_split(np.random.default_rng(0).permutation(ys), NFOLD) # ------------------------------------------------------------- các biến thể def votes_resid(zoo, Xtr, res_tr, mod_tr, Xte, mod_te): """Thành viên hồi quy PHẦN DƯ; dự báo cuối = Sev_mod + hiệu chỉnh.""" out = np.zeros((len(zoo), Xte.shape[0]), dtype="i1") for i, fam in enumerate(zoo): make = fam[-1] try: m = make() m.fit(Xtr, res_tr) mu_tr = np.asarray(m.predict(Xtr)).ravel() edges, sd = conditional_spread(mu_tr, res_tr - mu_tr) mu_te = np.asarray(m.predict(Xte)).ravel() sg_te = apply_spread(mu_te, edges, sd) pred = mod_te[:, None] + mu_te[:, None] + sg_te[:, None] * ZQ[None, :] out[i] = (pred >= 0.0).sum(axis=1) except Exception as e: print(" [warn] resid %s: %s" % (fam[0], e), flush=True) out[i] = -1 return out def votes_offset(zoo, Xtr, ytr, off_tr, Xte, off_te): """Thành viên phân loại có offset; (f, q) bỏ phiếu khi p_f >= tau_q.""" out = np.zeros((len(zoo), Xte.shape[0]), dtype="i1") for i, (nm, _, make) in enumerate(zoo): try: m = make() m.fit(Xtr, ytr, off_tr) p = m.predict_proba(Xte, off_te) out[i] = (p[:, None] >= QLEV[None, :]).sum(axis=1) except Exception as e: print(" [warn] offset %s: %s" % (nm, e), flush=True) out[i] = -1 return out def run_cell(ip, wname, Xtr, sev_tr, mod_tr, evt_tr, yr_tr, Xte, mod_te, names): t0 = time.time() sd = Xtr.std(axis=0) col = sd > 1e-12 nm = [n for n, k in zip(names, col) if k] Xtr = np.ascontiguousarray(Xtr[:, col]); Xte = np.ascontiguousarray(Xte[:, col]) zr = OFF.resid_families(Z.zoo_full(), nm) zo = OFF.offset_zoo() res_tr = sev_tr - mod_tr fl = folds_of(yr_tr) # --- p_raw: ngoài mẫu trên tập huấn luyện, khớp đầy đủ cho tập kiểm định praw_oof = np.zeros(Xtr.shape[0]) for f in fl: te = np.isin(yr_tr, f) if te.all() or not te.any(): continue praw_oof[te] = OFF.RawProb(nm).fit(Xtr[~te], evt_tr[~te] ).predict_proba(Xtr[te]) praw_val = OFF.RawProb(nm).fit(Xtr, evt_tr).predict_proba(Xte) # --- B1 neo bằng Sev, và B2 neo bằng xác suất b1_oof = np.zeros((len(zr), Xtr.shape[0]), dtype="i1") b2_oof = np.zeros((len(zo), Xtr.shape[0]), dtype="i1") for f in fl: te = np.isin(yr_tr, f) if te.all() or not te.any(): continue b1_oof[:, te] = votes_resid(zr, Xtr[~te], res_tr[~te], mod_tr[~te], Xtr[te], mod_tr[te]) b2_oof[:, te] = votes_offset(zo, Xtr[~te], evt_tr[~te], OFF.logit(praw_oof[~te]), Xtr[te], OFF.logit(praw_oof[te])) b1_val = votes_resid(zr, Xtr, res_tr, mod_tr, Xte, mod_te) b2_val = votes_offset(zo, Xtr, evt_tr, OFF.logit(praw_oof), Xte, OFF.logit(praw_val)) print("[%s p%3d] %d họ neo-Sev + %d họ offset, %.0fs" % (wname, ip, len(zr), len(zo), time.time() - t0), flush=True) return ip, wname, praw_oof, praw_val, b1_oof, b1_val, b2_oof, b2_val 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("--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) names, _ = agg_names(store, args.fset) zr = OFF.resid_families(Z.zoo_full(), names) zo = OFF.offset_zoo() print("[s1b] %d ô; B1 neo-Sev %d họ = %d thành viên; B2 offset %d họ = %d " "thành viên" % (pts.size * len(C.WEEKS), len(zr), len(zr) * NQ, len(zo), len(zo) * NQ), flush=True) tr = weekly_data(store, pts, C.TRAIN_YEARS, args.fset) te = weekly_data(store, pts, C.VALID_YEARS, args.fset) jobs, meta = [], [] for wname in C.WEEK_NAMES: a, b = tr[wname], te[wname] for p in pts: mt, mv = a["point"] == p, b["point"] == p ok_t = (np.isfinite(a["agg"][mt]).all(axis=1) & np.isfinite(a["sev"][mt]) & np.isfinite(a["sev_mod"][mt])) ok_v = (np.isfinite(b["agg"][mv]).all(axis=1) & np.isfinite(b["sev_mod"][mv]) & np.isfinite(b["evt"][mv])) jobs.append((int(p), wname, a["agg"][mt][ok_t], a["sev"][mt][ok_t], a["sev_mod"][mt][ok_t], a["evt"][mt][ok_t].astype("f4"), a["year"][mt][ok_t], b["agg"][mv][ok_v], b["sev_mod"][mv][ok_v], names)) meta.append((int(p), wname, a["evt"][mt][ok_t] > 0.5, b["evt"][mv][ok_v] > 0.5)) got = Parallel(n_jobs=args.jobs, verbose=5)(delayed(run_cell)(*j) for j in jobs) obs_tr = {("%d|%s" % (p, w)): o for p, w, o, _ in meta} obs_va = {("%d|%s" % (p, w)): o for p, w, _, o in meta} store_ = {} dump = {} for ip, w, po, pv, b1o, b1v, b2o, b2v in got: k = "%d|%s" % (ip, w) store_[k] = dict(praw_oof=po, praw_val=pv, b1_oof=b1o, b1_val=b1v, b2_oof=b2o, b2_val=b2v) for nmk, arr in store_[k].items(): dump["%s|%s" % (k, nmk)] = arr dump[k + "|obs_oof"] = obs_tr[k]; dump[k + "|obs_val"] = obs_va[k] dump["fam_b1"] = np.array([z[0] for z in zr]) dump["fam_b2"] = np.array([z[0] for z in zo]) np.savez_compressed(OUT_NPZ, **dump) print("[dump] -> %s" % OUT_NPZ, flush=True) def p_of(v): good = v >= 0 return np.where(good.any(axis=0), np.where(good, v, 0).sum(axis=0) / np.maximum(good.sum(axis=0) * NQ, 1), np.nan) # biến thể A (không neo) lấy lại từ lần chạy trước, cùng 10 họ prev = os.path.join(RES_DIR, "exp3_week_s1_votes.npz") A = None if os.path.exists(prev): za = np.load(prev, allow_pickle=True) fa = [str(x) for x in za["families"]] sel = np.array([fa.index(n) for n in [z[0] for z in zr] if n in fa]) A = (za, sel) res = {"config": {"points": pts.tolist(), "weeks": C.WEEK_NAMES, "fset": args.fset, "n_quantile": NQ, "n_family_b1": len(zr), "n_family_b2": len(zo)}, "models_b1": [{"name": z[0], "desc": z[2]} for z in zr], "models_b2": [{"name": n, "desc": d} for n, d, _ in zo], "by_week": {}, "family_b2": {}} for w in C.WEEK_NAMES: ks = [k for k in store_ if k.split("|")[1] == w] if not ks: continue oo = np.concatenate([obs_tr[k] for k in ks]) ov = np.concatenate([obs_va[k] for k in ks]) blk = {} cands = { "p_raw (EC gốc)": (np.concatenate([store_[k]["praw_oof"] for k in ks]), np.concatenate([store_[k]["praw_val"] for k in ks])), "B1 neo bằng Sev": (np.concatenate([p_of(store_[k]["b1_oof"]) for k in ks]), np.concatenate([p_of(store_[k]["b1_val"]) for k in ks])), "B2 neo bằng xác suất": ( np.concatenate([p_of(store_[k]["b2_oof"]) for k in ks]), np.concatenate([p_of(store_[k]["b2_val"]) for k in ks])), } if A is not None: za, sel = A def pa(kind): out = [] for k in ks: v = za["%s|%s" % (k, kind)][sel] good = v >= 0 out.append(np.where(good.any(axis=0), np.where(good, v, 0).sum(axis=0) / np.maximum(good.sum(axis=0) * NQ, 1), np.nan)) return np.concatenate(out) # lần chạy trước lọc NaN kèm cả chuỗi lead day nên số hàng lệch vài # đơn vị; phải lấy luôn quan trắc CỦA CHÍNH lần chạy đó cands["A không neo"] = ( pa("oof"), pa("val"), np.concatenate([za[k + "|obs_oof"] for k in ks]), np.concatenate([za[k + "|obs_val"] for k in ks])) for nm, val in cands.items(): po, pv, oo_, ov_ = (val if len(val) == 4 else (val[0], val[1], oo, ov)) p1 = pick_p1(po, oo_) g = gate_scores(pv, ov_, p1) g["AUC"] = float(M.auc_score(pv, ov_)) g["AUC_oof"] = float(M.auc_score(po, oo_)) blk[nm] = g res["by_week"][w] = blk print("[%s] " % w + " | ".join( "%s AUC %.3f đóng %.0f%% NPV %.3f rò %.0f%%" % (nm, g["AUC"], 100 * g["closure"], g["NPV"], 100 * g["leakage"]) for nm, g in blk.items()), flush=True) for i, (nm, desc, _) in enumerate(zo): pv = np.concatenate([store_[k]["b2_val"][i].astype("f8") for k in store_]) ov = np.concatenate([obs_va[k] for k in store_]) good = pv >= 0 res["family_b2"][nm] = { "desc": desc, "AUC_val": float(M.auc_score(pv[good], ov[good])) if good.any() else np.nan, "fail_frac": float(1.0 - good.mean())} 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()