#!/usr/bin/env python3 """EXP3 (bản cập nhật): N mô hình AI tại mỗi ô lưới, quyết định hai tầng bằng p = M/N. Đặc tả của người dùng, thực hiện đúng cấu trúc: * 21 điểm (3 điểm x 7 vùng), mỗi điểm có bộ mô hình riêng; * target là Sev_obs, giữ liên tục và có dấu (cả âm lẫn không âm); * M trong N thành viên cho yEHDE -> p = M/N, q = 1 - p; * tầng 1: p < p1 => nEHDE, bỏ mẫu. p >= p1 => giữ lại; * tầng 2: khớp lại trên tập giữ lại, p2 > nguong => yEHDE; * đánh giá trên TOÀN BỘ tập mẫu kiểm định, gộp kết quả của cả hai tầng -- mẫu bị tầng 1 kết luận nEHDE thì chỉ cần một lần dự báo; * xuất bảng contingency cho TỪNG điểm trong 21 điểm, kèm danh mục mô hình. VÌ SAO THÀNH VIÊN LÀ CẶP (HỌ MÔ HÌNH, MỨC PHÂN VỊ). Đo trực tiếp trên một ô lưới: Sev_obs quan trắc >= 0 ở 12.65% số hàng, nhưng dự báo điểm của TÁM họ mô hình khác nhau đều chỉ vượt 0 ở dưới 0.9% (Ridge 0.27, ElasticNet 0.06, LGBM 0.84, LGBM-RF 0.00, ExtraTrees 0.24, RandomForest 0.45, kNN 0.19, MLP 0.86). Nguyên nhân không phải các mô hình giống nhau mà là mọi mô hình đều ước lượng KỲ VỌNG CÓ ĐIỀU KIỆN của Sev, vốn gần như luôn âm. Đếm phiếu trên các kỳ vọng vì thế luôn cho M/N ~ 0 và không ngưỡng p nào cứu được. Mỗi họ ở đây do đó sinh ra một chùm thành viên theo mức phân vị: thành viên (f, q) dự báo mu_f + sigma_f * z_q, với sigma ước lượng theo thập phân vị của mu nên bắt được tính không đồng nhất phương sai. Công thức p = M/N giữ nguyên; chỉ có bản chất một "thành viên" là khác. Biến thể `mean` (đếm phiếu trên đúng N kỳ vọng của N họ) vẫn được chạy song song để mức suy biến hiện ra bằng số chứ không phải bằng lập luận. """ 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 train_ai as T import pointwise as PW RES_DIR = os.path.join(C.OUT, "AI_Results") OUT_JSON = os.path.join(RES_DIR, "exp3_multi.json") NFOLD = 4 NQ = 9 # số mức phân vị cho mỗi họ mô hình SPREAD_BINS = 10 # số thập phân vị của mu để ước lượng sigma # ------------------------------------------------------------ danh mục mô hình def model_zoo(): """Bảy họ mô hình AI, mỗi họ một cơ chế học khác hẳn nhau.""" from sklearn.linear_model import Ridge, ElasticNet from sklearn.ensemble import ExtraTreesRegressor from sklearn.neighbors import KNeighborsRegressor from sklearn.neural_network import MLPRegressor from sklearn.pipeline import make_pipeline from sklearn.preprocessing import StandardScaler import lightgbm as lgb return [ ("Ridge", "hồi quy tuyến tính, phạt L2", lambda: make_pipeline(StandardScaler(), Ridge(alpha=10.0))), ("ElasticNet", "hồi quy tuyến tính thưa, phạt L1+L2", lambda: make_pipeline(StandardScaler(), ElasticNet(alpha=0.01, l1_ratio=0.5, max_iter=2000))), ("LGBM_gbdt", "cây tăng cường độ dốc (gradient boosting)", lambda: lgb.LGBMRegressor(n_estimators=300, num_leaves=15, learning_rate=0.05, min_child_samples=100, subsample=0.8, subsample_freq=1, colsample_bytree=0.8, reg_lambda=5.0, n_jobs=1, verbose=-1)), ("LGBM_rf", "rừng ngẫu nhiên dựng bằng LightGBM", lambda: lgb.LGBMRegressor(boosting_type="rf", n_estimators=200, num_leaves=63, subsample=0.7, subsample_freq=1, colsample_bytree=0.7, min_child_samples=50, n_jobs=1, verbose=-1)), ("ExtraTrees", "cây cực kỳ ngẫu nhiên (extremely randomised trees)", lambda: ExtraTreesRegressor(n_estimators=100, min_samples_leaf=20, max_features=0.5, n_jobs=1, random_state=0)), ("kNN", "láng giềng gần nhất, k = 25", lambda: make_pipeline(StandardScaler(), KNeighborsRegressor(n_neighbors=25, n_jobs=1))), ("MLP", "mạng nơ-ron truyền thẳng, 64-32", lambda: make_pipeline(StandardScaler(), MLPRegressor(hidden_layer_sizes=(64, 32), max_iter=60, early_stopping=True, random_state=0))), ] ZOO = [(n, d) for n, d, _ in model_zoo()] QLEVELS = (np.arange(NQ) + 0.5) / NQ # ------------------------------------------------------------------ ước lượng def conditional_spread(mu, resid): """Độ lệch chuẩn của phần dư theo thập phân vị của mu (bắt heteroscedasticity).""" edges = np.quantile(mu, np.linspace(0, 1, SPREAD_BINS + 1)) edges[0] -= 1e-6; edges[-1] += 1e-6 b = np.clip(np.digitize(mu, edges) - 1, 0, SPREAD_BINS - 1) sd = np.full(SPREAD_BINS, np.nanstd(resid) if resid.size else 1.0) for k in range(SPREAD_BINS): m = b == k if m.sum() >= 50: sd[k] = float(np.std(resid[m])) return edges, np.maximum(sd, 1e-3) def apply_spread(mu, edges, sd): b = np.clip(np.digitize(mu, edges) - 1, 0, SPREAD_BINS - 1) return sd[b] def fit_family(make, Xtr, ytr, Xte): """Kỳ vọng dự báo của một họ, cùng độ tán phần dư ước lượng trong tập huấn luyện.""" m = make() m.fit(Xtr, ytr) mu_tr = m.predict(Xtr) edges, sd = conditional_spread(mu_tr, ytr - mu_tr) mu_te = m.predict(Xte) return mu_te, apply_spread(mu_te, edges, sd) def ensemble_votes(Xtr, ytr, Xte, zoo): """p = M/N cho hai cách đọc 'thành viên', từ cùng một lần khớp mô hình. Trả về (p_quantile, p_mean, mu_all): thành viên theo cặp (họ, mức phân vị), thành viên chỉ là kỳ vọng của mỗi họ, và ma trận kỳ vọng để chẩn đoán. """ nf = len(zoo) mu = np.empty((nf, Xte.shape[0]), dtype="f4") sg = np.empty((nf, Xte.shape[0]), dtype="f4") for i, (_, _, make) in enumerate(zoo): mu[i], sg[i] = fit_family(make, Xtr, ytr, Xte) z = norm.ppf(QLEVELS).astype("f4") # (NQ,) # thành viên (f, q) dự báo mu_f + sigma_f * z_q -> đếm số thành viên >= 0 hit = (mu[:, :, None] + sg[:, :, None] * z[None, None, :]) >= 0.0 return hit.mean(axis=(0, 2)), (mu >= 0.0).mean(axis=0), mu def oof_votes(X, y, year, zoo, nfold=NFOLD): """p ngoài mẫu trên toàn bộ tập huấn luyện, chia theo năm.""" years = np.unique(year) folds = np.array_split(np.random.default_rng(0).permutation(years), nfold) pq = np.zeros(X.shape[0]); pm = np.zeros(X.shape[0]) for f in folds: te = np.isin(year, f) if te.all() or not te.any(): continue a, b, _ = ensemble_votes(X[~te], y[~te], X[te], zoo) pq[te], pm[te] = a, b return pq, pm # ------------------------------------------------------------------ một điểm def run_point(ip, Xtr, sev_tr, evt_tr, wtr, Xte, p1_grid, seed0): """Cả hai tầng EXP3 tại một ô lưới, cho mọi ngưỡng p1 trong lưới quét.""" 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, m_oof = oof_votes(Xtr, sev_tr, wtr, zoo) q_val, m_val, mu_val = ensemble_votes(Xtr, sev_tr, Xte, zoo) out = {"stage1": (q_oof, q_val, m_oof, m_val), "mu_frac_pos": float((mu_val >= 0).mean()), "stage2": {}} for p1 in p1_grid: keep = q_oof >= p1 if keep.sum() < 500 or evt_tr[keep].sum() < 30: continue o2 = np.zeros(Xtr.shape[0]); v2 = np.zeros(Xte.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) v2 = b out["stage2"]["%.3f" % p1] = (o2, v2) print("[point %3d] %d hàng, %d ngưỡng p1, %.0fs" % (ip, Xtr.shape[0], len(out["stage2"]), time.time() - t0), flush=True) return ip, out # ---------------------------------------------------------------- báo cáo 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)}) 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]) 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("[exp3] %d điểm, N = %d họ x %d mức phân vị = %d thành viên" % (pts.size, len(ZOO), NQ, len(ZOO) * NQ), flush=True) for n, d in ZOO: print(" %-12s %s" % (n, d), flush=True) Xtr, ytr_, itr = store.build(pts, C.TRAIN_YEARS, args.fset) Xte, yte_, ite = store.build(pts, C.VALID_YEARS, args.fset) ok_tr = (np.isfinite(Xtr).all(axis=1) & np.isfinite(ytr_["sev"]) & np.isfinite(ytr_["evt"])) ok_te = np.isfinite(Xte).all(axis=1) & np.isfinite(yte_["evt"]) print("[data] train %d hàng, valid %d hàng, base rate %.4f" % (ok_tr.sum(), ok_te.sum(), ytr_["evt"][ok_tr].mean()), flush=True) got = Parallel(n_jobs=args.jobs, verbose=5)( delayed(run_point)( int(p), Xtr[ok_tr & (itr["point"] == p)], ytr_["sev"][ok_tr & (itr["point"] == p)], ytr_["evt"][ok_tr & (itr["point"] == p)], itr["year"][ok_tr & (itr["point"] == p)], Xte[ok_te & (ite["point"] == p)], args.p1_grid, args.seed + int(p)) for p in pts) res = {"config": {"points": pts.tolist(), "n_family": len(ZOO), "n_quantile": NQ, "n_members": len(ZOO) * NQ, "fset": args.fset, "p1_grid": args.p1_grid, "p2_grid": args.p2_grid}, "models": [{"name": n, "desc": d} for n, d in ZOO], "by_point": {}, "sweep": {}, "degenerate": {}} idx_te = {int(p): (ok_te & (ite["point"] == p)) for p in pts} for ip, r in got: m = idx_te[ip] obs = yte_["evt"][m] > 0.5 res["degenerate"][str(ip)] = { "mean_vote_frac_pos": r["mu_frac_pos"], "obs_frac_pos": float(obs.mean()), "quantile_vote_mean_p": float(np.mean(r["stage1"][1])), "mean_vote_mean_p": float(np.mean(r["stage1"][3]))} # Quét (p1, p2). Cấu hình được CHỐT TRÊN TRƯỜNG NGOÀI MẪU CỦA CÁC NĂM HUẤN # LUYỆN, không bao giờ trên tập kiểm định: chọn ngưỡng theo CSI của tập kiểm # định chính là dò tham số trên tập test và sẽ cho con số lạc quan giả. # Bảng trên tập kiểm định vẫn được xuất, nhưng chỉ để đối chiếu, và nó là # CẬN TRÊN chứ không phải kết quả có thể đạt được trong vận hành. idx_tr = {int(p): (ok_tr & (itr["point"] == p)) for p in pts} best = None for p1 in args.p1_grid: for p2 in args.p2_grid: k = "%.3f" % p1 pr_o, ob_o, pr_v, ob_v = [], [], [], [] for ip, r in got: if k not in r["stage2"]: continue pr_o.append((r["stage1"][0] >= p1) & (r["stage2"][k][0] > p2)) ob_o.append(ytr_["evt"][idx_tr[ip]] > 0.5) pr_v.append((r["stage1"][1] >= p1) & (r["stage2"][k][1] > p2)) ob_v.append(yte_["evt"][idx_te[ip]] > 0.5) if not pr_o: continue so = contingency(np.concatenate(pr_o), np.concatenate(ob_o)) sv = contingency(np.concatenate(pr_v), np.concatenate(ob_v)) tag = "p1=%.2f,p2=%.2f" % (p1, p2) res["sweep"][tag] = sv res.setdefault("sweep_oof", {})[tag] = so if best is None or so["CSI"] > best[0]: best = (so["CSI"], p1, p2) if best: _, p1b, p2b = best res["best"] = {"p1": p1b, "p2": p2b} print("\n[chốt] cấu hình chọn trên trường NGOÀI MẪU của năm huấn luyện: " "p1 = %.2f, p2 = %.2f (CSI ngoài mẫu %.3f)" % (p1b, p2b, best[0]), flush=True) sv = res["sweep"]["p1=%.2f,p2=%.2f" % (p1b, p2b)] print("[kiểm định] POD %.3f POFD %.3f FAR %.3f FBI %.3f CSI %.3f PC %.3f" % (sv["POD"], sv["POFD"], sv["FAR"], sv["FBI"], sv["CSI"], sv["PC"]), flush=True) for ip, r in got: k = "%.3f" % p1b if k not in r["stage2"]: continue m = idx_te[ip] q_val = r["stage1"][1] pred = (q_val >= p1b) & (r["stage2"][k][1] > p2b) obs = yte_["evt"][m] > 0.5 s = contingency(pred, obs) s["region"] = int(store.grid.reg[ip]) s["lat"] = float(store.st["lat"][ip]); s["lon"] = float(store.st["lon"][ip]) s["gate_closed_frac"] = float((q_val < p1b).mean()) res["by_point"][str(ip)] = s 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()