#!/usr/bin/env python3 """EXP3 tầng 1 mở rộng: cổng nEHDE với N lớn, gộp theo lead week. CHỈ LÀM TẦNG 1. Tầng 2 để sau vì nó nhiều khả năng dùng bộ mô hình khác hẳn. Câu hỏi duy nhất của tầng 1: "mẫu này gần như chắc chắn nEHDE?" Trả lời bằng p = M/N. Bản trước dùng N = 7 họ x 9 phân vị = 63 thành viên; bản này dùng N = 28 họ x 9 phân vị = 252 thành viên, thêm bảy cấu hình ANN (số lớp, số nút, hàm kích hoạt, dropout, phạt trọng số khác nhau), ba CNN 1-D chạy dọc chuỗi lead day trong tuần, cùng các họ tuyến tính bền vững, cây khác cơ chế và SVR. Đầu ra được thiết kế để trả lời được câu "cổng loại bỏ nEHDE chính xác đến đâu": * ĐỘ ĐÓNG (closure) -- bao nhiêu phần trăm mẫu bị kết luận nEHDE; * ĐỘ TINH KHIẾT (NPV) -- trong số bị đóng, bao nhiêu phần trăm thật sự nEHDE, luôn đặt cạnh mốc tham chiếu 1 - base_rate của việc đóng bừa; * RÒ RỈ (leakage) -- bao nhiêu phần trăm SỰ KIỆN thật bị cổng ném đi, con số quyết định cổng có dùng được hay không; * AUC / EDI / SEDI của chính điểm số p, không phụ thuộc ngưỡng và base rate. Ma trận số phiếu theo TỪNG HỌ được lưu ra .npz, nên mọi tổ hợp con của 28 họ đều đánh giá lại được mà không phải khớp lại mô hình. """ 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 from exp3_multi import conditional_spread, apply_spread from exp3_week import MAX_VARS RES_DIR = os.path.join(C.OUT, "AI_Results") OUT_JSON = os.path.join(RES_DIR, "exp3_week_s1.json") OUT_NPZ = os.path.join(RES_DIR, "exp3_week_s1_votes.npz") NFOLD = 4 NQ = 9 QLEVELS = (np.arange(NQ) + 0.5) / NQ ZQ = norm.ppf(QLEVELS) P1_GRID = np.round(np.arange(0.02, 0.62, 0.02), 3) MAX_LEAK = 0.05 # cổng không được ném đi quá 5 % số sự kiện thật # ------------------------------------------------------ dữ liệu gộp theo tuần def weekly_data(store, points, years, fset): """Trả về, theo từng tuần: đặc trưng đã gộp, CHUỖI LEAD DAY, và target.""" 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, :, :] # (ni, T, npt, nf) 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) # chuỗi lead day cho CNN: (mẫu, ngày, biến) seq = Xw.transpose(0, 2, 1, 3).reshape(ni * npt, len(lts), nf) 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[wname] = dict(agg=agg, seq=seq, sev=sw, evt=ew.astype("f4"), year=yy, point=pp) return out # ------------------------------------------------------------ phiếu theo họ def family_votes(zoo, Atr, Str, ytr, Ate, Ste): """Số thành viên (trong 9 mức phân vị) cho yEHDE, TÁCH RIÊNG theo từng họ.""" out = np.zeros((len(zoo), Ate.shape[0]), dtype="i1") for i, (_, kind, _, make) in enumerate(zoo): Xtr, Xte = (Str, Ste) if kind == "seq" else (Atr, Ate) try: m = make() m.fit(Xtr, ytr) mu_tr = np.asarray(m.predict(Xtr)).ravel() edges, sd = conditional_spread(mu_tr, ytr - mu_tr) mu_te = np.asarray(m.predict(Xte)).ravel() sg_te = apply_spread(mu_te, edges, sd) out[i] = ((mu_te[:, None] + sg_te[:, None] * ZQ[None, :]) >= 0.0 ).sum(axis=1) except Exception as e: # một họ hỏng không được giết cả ô print(" [warn] %s: %s" % (zoo[i][0], e), flush=True) out[i] = -1 return out def oof_family_votes(zoo, A, S, y, year, nfold=NFOLD): years = np.unique(year) folds = np.array_split(np.random.default_rng(0).permutation(years), nfold) out = np.zeros((len(zoo), A.shape[0]), dtype="i1") for f in folds: te = np.isin(year, f) if te.all() or not te.any(): continue out[:, te] = family_votes(zoo, A[~te], S[~te], y[~te], A[te], S[te]) return out def run_cell(ip, wname, Atr, Str, sev_tr, evt_tr, yr_tr, Ate, Ste): zoo = Z.zoo_full() t0 = time.time() sd = Atr.std(axis=0) col = sd > 1e-12 Atr = np.ascontiguousarray(Atr[:, col]); Ate = np.ascontiguousarray(Ate[:, col]) v_oof = oof_family_votes(zoo, Atr, Str, sev_tr, yr_tr) v_val = family_votes(zoo, Atr, Str, sev_tr, Ate, Ste) print("[%s p%3d] %d họ, train %d (%d dương), %.0fs" % (wname, ip, len(zoo), Atr.shape[0], int(evt_tr.sum()), time.time() - t0), flush=True) return ip, wname, v_oof, v_val # ------------------------------------------------------------- chẩn đoán cổng def gate_scores(p, obs, p1): """Chất lượng của cổng 'p < p1 => nEHDE' tại một ngưỡng.""" closed = p < p1 n = obs.size n_closed = int(closed.sum()) lost = int(obs[closed].sum()) # sự kiện thật bị ném đi npv = 1.0 - lost / n_closed if n_closed else np.nan a, b, c, d = M.contingency(~closed, obs) # 'cổng mở' như một dự báo yEHDE s = M.scores(a, b, c, d) return {"p1": float(p1), "closure": n_closed / n, "n_closed": n_closed, "NPV": npv, "NPV_ref": 1.0 - float(obs.mean()), "leakage": lost / max(int(obs.sum()), 1), "n_lost": lost, "POD": s["POD"], "POFD": s["POFD"], "FAR": s["FAR"], "CSI": s["CSI"], "EDI": s["EDI"], "SEDI": s["SEDI"], "n": n, "base_rate": float(obs.mean())} def pick_p1(p, obs, max_leak=MAX_LEAK): """Ngưỡng đóng được NHIỀU mẫu nhất mà vẫn không ném đi quá max_leak sự kiện.""" best = 0.0 for p1 in P1_GRID: if p[obs.astype(bool)].size == 0: break leak = float((p[obs.astype(bool)] < p1).mean()) if leak <= max_leak: best = float(p1) return best 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() zoo = Z.zoo_full() nbase = len(Z.zoo_base()) store = D.Store() pts = D.sample_points(store.grid, args.per_region, args.seed) print("[s1] %d điểm x %d tuần = %d ô, N = %d họ x %d phân vị = %d thành viên " "(bản gốc: %d họ = %d thành viên)" % (pts.size, len(C.WEEKS), pts.size * len(C.WEEKS), len(zoo), NQ, len(zoo) * NQ, nbase, nbase * 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["seq"][mt]).all(axis=(1, 2)) & np.isfinite(a["sev"][mt])) ok_v = (np.isfinite(b["agg"][mv]).all(axis=1) & np.isfinite(b["seq"][mv]).all(axis=(1, 2)) & np.isfinite(b["evt"][mv])) jobs.append((int(p), wname, a["agg"][mt][ok_t], a["seq"][mt][ok_t], a["sev"][mt][ok_t], a["evt"][mt][ok_t], a["year"][mt][ok_t], b["agg"][mv][ok_v], b["seq"][mv][ok_v])) meta.append((int(p), wname, a["evt"][mt][ok_t] > 0.5, b["evt"][mv][ok_v] > 0.5)) print("[data] %d ô, mỗi ô %d mẫu huấn luyện / %d mẫu kiểm định" % (len(jobs), jobs[0][2].shape[0], jobs[0][7].shape[0]), flush=True) got = Parallel(n_jobs=args.jobs, verbose=5)(delayed(run_cell)(*j) for j in jobs) votes = {("%d|%s" % (ip, w)): (vo, vv) for ip, w, vo, vv in got} obs_tr = {("%d|%s" % (ip, w)): o for ip, w, o, _ in meta} obs_va = {("%d|%s" % (ip, w)): o for ip, w, _, o in meta} dump = {} for k, (vo, vv) in votes.items(): dump[k + "|oof"] = vo; dump[k + "|val"] = vv dump[k + "|obs_oof"] = obs_tr[k]; dump[k + "|obs_val"] = obs_va[k] dump["families"] = np.array([z[0] for z in zoo]) np.savez_compressed(OUT_NPZ, **dump) print("[dump] phiếu theo từng họ -> %s" % OUT_NPZ, flush=True) def p_of(v, sel): """p = M/N cho một tổ hợp họ; họ hỏng (mã -1) bị loại khỏi mẫu số.""" sub = v[sel] good = sub >= 0 return np.where(good.any(axis=0), np.where(good, sub, 0).sum(axis=0) / np.maximum(good.sum(axis=0) * NQ, 1), np.nan) sel_base = np.arange(nbase) sel_full = np.arange(len(zoo)) res = {"config": {"points": pts.tolist(), "weeks": C.WEEK_NAMES, "fset": args.fset, "n_quantile": NQ, "n_family_base": nbase, "n_family_full": len(zoo), "n_member_base": nbase * NQ, "n_member_full": len(zoo) * NQ, "max_leak": MAX_LEAK}, "models": [{"name": n, "kind": k, "desc": d, "new": i >= nbase} for i, (n, k, d, _) in enumerate(zoo)], "by_week": {}, "curve": {}, "by_cell": {}, "family": {}} for wname in C.WEEK_NAMES: ks = [k for k in votes if k.split("|")[1] == wname] 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 = {} for tag, sel in (("base", sel_base), ("full", sel_full)): po = np.concatenate([p_of(votes[k][0], sel) for k in ks]) pv = np.concatenate([p_of(votes[k][1], sel) for k in ks]) 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)) g["p1_from_oof"] = p1 g["closure_oof"] = float((po < p1).mean()) g["leakage_oof"] = float((po[oo] < p1).mean()) if oo.any() else np.nan blk[tag] = g res["curve"]["%s|%s" % (wname, tag)] = [ gate_scores(pv, ov, q) for q in P1_GRID] res["by_week"][wname] = blk b, f = blk["base"], blk["full"] print("[%s] base: AUC %.3f p1 %.2f đóng %.1f%% NPV %.3f rò %.1f%% | " "full: AUC %.3f p1 %.2f đóng %.1f%% NPV %.3f rò %.1f%% (mốc %.3f)" % (wname, b["AUC"], b["p1"], 100 * b["closure"], b["NPV"], 100 * b["leakage"], f["AUC"], f["p1"], 100 * f["closure"], f["NPV"], 100 * f["leakage"], f["NPV_ref"]), flush=True) # từng ô lưới tại ngưỡng của tổ hợp đầy đủ p1f = blk["full"]["p1_from_oof"] for k in ks: ip = int(k.split("|")[0]) pv = p_of(votes[k][1], sel_full) g = gate_scores(pv, obs_va[k], p1f) g.update({"region": int(store.grid.reg[ip]), "lat": float(store.st["lat"][ip]), "lon": float(store.st["lon"][ip]), "AUC": float(M.auc_score(pv, obs_va[k]))}) res["by_cell"][k] = g # đóng góp riêng của từng họ: AUC của phiếu chính họ đó, gộp toàn bộ 126 ô for i, (nm, kind, desc, _) in enumerate(zoo): po = np.concatenate([votes[k][0][i].astype("f8") for k in votes]) pv = np.concatenate([votes[k][1][i].astype("f8") for k in votes]) oo = np.concatenate([obs_tr[k] for k in votes]) ov = np.concatenate([obs_va[k] for k in votes]) good_o, good_v = po >= 0, pv >= 0 res["family"][nm] = { "kind": kind, "desc": desc, "new": i >= nbase, "AUC_oof": float(M.auc_score(po[good_o], oo[good_o])) if good_o.any() else np.nan, "AUC_val": float(M.auc_score(pv[good_v], ov[good_v])) if good_v.any() else np.nan, "vote_mean": float(np.mean(pv[good_v]) / NQ) if good_v.any() else np.nan, "fail_frac": float(1.0 - good_v.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()