#!/usr/bin/env python3 """GĐ1: chốt kiến trúc triển khai toàn miền -- per-point, gộp vùng, hay gộp toàn quốc. Ba biến thể được chấm trên CÙNG một tập điểm đánh giá và CÙNG một danh mục họ mô hình, nên chênh lệch đọc được là chênh lệch của kiến trúc chứ không phải của mô hình. point -- mỗi ô lưới một bộ mô hình riêng, khớp trên 1664 mẫu của chính nó. Đây là kiến trúc đang dùng ở mục 3h. region -- mỗi vùng khí hậu một bộ, khớp trên toàn bộ ô lưới của vùng đó (85 000 - 175 000 mẫu). national -- một bộ cho cả nước, khớp trên 541 ô (900 000 mẫu). Ở CẢ BA biến thể, NGƯỠNG vẫn chốt RIÊNG TỪNG Ô LƯỚI trên trường ngoài mẫu. Đó là điểm mấu chốt: mọi kết quả từ AI11 tới AI25 cho thấy thông tin cục bộ nằm ở ngưỡng chứ không nằm ở hệ số mô hình, nên gộp mô hình mà giữ ngưỡng riêng là phép thử đúng. TÍNH KHÁI QUÁT KHÔNG GIAN. Mọi số liệu trước đây đến từ 21 điểm lấy với seed 0. Tập đánh giá ở đây gồm 21 điểm seed 0 CỘNG 21 điểm seed 1, và kết quả được tách riêng cho hai nhóm, để biết điểm vận hành có chuyển sang tập điểm khác hay không. BỘ NHỚ. Ma trận toàn quốc TRƯỚC khi gộp tuần là 1664 x 46 x 541 x 70 x 4 B = 11,6 GB, không giữ nổi. Vì thế initialisation được duyệt theo chunk và gộp tuần ngay trong từng chunk; kết quả sau gộp chỉ còn 292 MB. """ 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 ai_data as D import exp3_offset as OFF from exp3_multi import conditional_spread, apply_spread from exp3_week import MAX_VARS from exp3_week_s1 import pick_p1 RES_DIR = os.path.join(C.OUT, "AI_Results", "phase1") NFOLD = 4 NQ = 9 QLEV = (np.arange(NQ) + 0.5) / NQ ZQ = norm.ppf(QLEV) CHUNK = 64 # Từng tuần riêng là sản phẩm nghiệp vụ người dùng cần (dự báo "tuần thứ 4", không # phải "ba tuần tới gộp lại"); các cửa sổ cộng dồn giữ lại để đối chiếu. WINDOWS = {"W1": [0], "W2": [1], "W3": [2], "W4": [3], "W5": [4], "W6": [5], "W1-2": [0, 1], "W1-3": [0, 1, 2], "W1-6": [0, 1, 2, 3, 4, 5]} # ------------------------------------------------- danh mục dùng chung 3 biến thể class PLS1D: """PLS trả về véc-tơ một chiều, cho khớp giao diện của các họ còn lại. Định nghĩa ở CẤP MODULE chứ không lồng trong hàm: bộ tham số vận hành phải pickle được, mà pickle không lưu nổi một lớp khai báo bên trong hàm. """ def __init__(self, n_components=8): from sklearn.cross_decomposition import PLSRegression from sklearn.pipeline import make_pipeline from sklearn.preprocessing import StandardScaler self.m_ = make_pipeline(StandardScaler(), PLSRegression(n_components=n_components)) def fit(self, X, y): self.m_.fit(X, y) return self def predict(self, X): return np.asarray(self.m_.predict(X)).ravel() def zoo_s1(n): """Tám họ hồi quy Sev. Siêu tham số của cây co giãn theo cỡ mẫu -- giữ nguyên min_samples_leaf khi n tăng 500 lần là để mô hình học thuộc lòng.""" from sklearn.linear_model import Ridge, ElasticNet, BayesianRidge, HuberRegressor from sklearn.ensemble import HistGradientBoostingRegressor from sklearn.pipeline import make_pipeline from sklearn.preprocessing import StandardScaler import lightgbm as lgb leaf = max(20, int(n / 80)) return [ ("Ridge", lambda: make_pipeline(StandardScaler(), Ridge(alpha=10.0))), ("ElasticNet", lambda: make_pipeline( StandardScaler(), ElasticNet(alpha=0.01, l1_ratio=0.5, max_iter=2000))), ("BayesRidge", lambda: make_pipeline(StandardScaler(), BayesianRidge())), ("Huber", lambda: make_pipeline( StandardScaler(), HuberRegressor(alpha=1e-3, max_iter=300))), ("PLS8", lambda: PLS1D(n_components=8)), ("LGBM_gbdt", lambda: lgb.LGBMRegressor( n_estimators=300, num_leaves=31, learning_rate=0.05, min_child_samples=leaf, subsample=0.8, subsample_freq=1, colsample_bytree=0.8, reg_lambda=5.0, n_jobs=1, verbose=-1)), ("LGBM_dart", lambda: lgb.LGBMRegressor( boosting_type="dart", n_estimators=200, num_leaves=31, learning_rate=0.05, min_child_samples=leaf, colsample_bytree=0.8, n_jobs=1, verbose=-1)), ("HistGB", lambda: HistGradientBoostingRegressor( max_iter=250, max_leaf_nodes=31, learning_rate=0.06, min_samples_leaf=leaf, l2_regularization=1.0, random_state=0)), ] def zoo_s2(n): """Sáu họ phân loại trực tiếp yEHDE cho tầng 2.""" from sklearn.linear_model import LogisticRegression from sklearn.ensemble import HistGradientBoostingClassifier from sklearn.pipeline import make_pipeline from sklearn.preprocessing import StandardScaler import lightgbm as lgb leaf = max(20, int(n / 80)) return [ ("LogReg_L2", lambda: make_pipeline( StandardScaler(), LogisticRegression(C=0.3, max_iter=800))), ("LogReg_EN", lambda: make_pipeline( StandardScaler(), LogisticRegression( penalty="elasticnet", C=0.5, l1_ratio=0.5, solver="saga", max_iter=400))), ("LGBM_clf", lambda: lgb.LGBMClassifier( n_estimators=250, num_leaves=31, learning_rate=0.05, min_child_samples=leaf, colsample_bytree=0.8, reg_lambda=5.0, n_jobs=1, verbose=-1)), ("LGBM_clf_slow", lambda: lgb.LGBMClassifier( n_estimators=400, num_leaves=15, learning_rate=0.02, min_child_samples=leaf * 2, colsample_bytree=0.7, reg_lambda=10.0, n_jobs=1, verbose=-1)), ("HistGB_clf", lambda: HistGradientBoostingClassifier( max_iter=250, max_leaf_nodes=31, learning_rate=0.06, min_samples_leaf=leaf, l2_regularization=1.0, random_state=0)), ("ANN_16", "ann"), ] # --------------------------------------------------------------- dữ liệu def build_window(store, points, years, fset, lts, chunk=CHUNK): """Gộp tuần ngay trong từng chunk initialisation để không giữ ma trận 11,6 GB.""" names = store.feature_names(fset) imax = [names.index(v) for v in MAX_VARS if v in names] isv = names.index("sv_mod") isel_all = np.nonzero(np.isin(store.init_year, years))[0] A, S, N, Y, P = [], [], [], [], [] for a in range(0, isel_all.size, chunk): isel = isel_all[a:a + chunk] X, y, _ = store.build_isel(points, isel, fset) ni, nk, npt, nf = isel.size, C.NLT, points.size, X.shape[1] X = X.reshape(ni, nk, npt, nf)[:, lts, :, :] sev = y["sev"].reshape(ni, nk, npt)[:, lts, :] evt = y["evt"].reshape(ni, nk, npt)[:, lts, :] A.append(np.concatenate( [np.nanmean(X, axis=1), np.nanmax(X[:, :, :, imax], axis=1), (X[:, :, :, isv] >= 0.0).sum(axis=1, dtype="f4")[:, :, None]], axis=2).reshape(ni * npt, -1).astype("f4")) S.append(np.nanmax(sev, axis=1).reshape(ni * npt).astype("f4")) N.append(np.nansum(evt, axis=1).reshape(ni * npt).astype("f4")) Y.append(np.broadcast_to(store.init_year[isel][:, None], (ni, npt)).reshape(ni * npt)) P.append(np.broadcast_to(points[None, :], (ni, npt)).reshape(ni * npt)) del X, sev, evt agg_names = (["m_" + s for s in names] + ["x_" + names[i] for i in imax] + ["n_hotday_mod"]) return (np.concatenate(A), np.concatenate(S), np.concatenate(N), np.concatenate(Y), np.concatenate(P), agg_names) # ------------------------------------------------------------------- khớp mô hình def _fit_s1(mk, nm, Xtr, ztr, Xte, thr_te): """Một họ, một lần khớp -> số phiếu của 9 thành viên phân vị. `thr_te` là ngưỡng bỏ phiếu THEO TỪNG HÀNG: với target ">= 1 ngày" nó là 0 (dự báo Sev cực đại vượt ngưỡng EHD), còn với target ">= k ngày" nó là k của chính ô lưới đó, vì k được chọn riêng cho mỗi ô theo khí hậu của ô. """ try: m = mk(); m.fit(Xtr, ztr) mu = np.asarray(m.predict(Xtr)).ravel() edges, sd = conditional_spread(mu, ztr - mu) mt = np.asarray(m.predict(Xte)).ravel() sg = apply_spread(mt, edges, sd) return ((mt[:, None] + sg[:, None] * ZQ[None, :] ) >= np.asarray(thr_te, "f8")[:, None]).sum(axis=1).astype("i1") except Exception as e: print(" [warn] s1 %s: %s" % (nm, e), flush=True) return np.full(Xte.shape[0], -1, dtype="i1") def _fit_s2(mk, nm, Xtr, ytr, Xte): try: if mk == "ann": big = Xtr.shape[0] > 20000 m = OFF.OffsetNet(hidden=(16,), act="relu", epochs=30 if big else 150, batch=2048 if big else 32, seed=0) m.fit(Xtr, ytr, np.zeros(Xtr.shape[0])) return np.clip(m.predict_proba(Xte, np.zeros(Xte.shape[0])), 0, 1) m = mk(); m.fit(Xtr, np.asarray(ytr).astype(int)) return np.clip(np.asarray(m.predict_proba(Xte)[:, 1], "f8"), 0, 1) except Exception as e: print(" [warn] s2 %s: %s" % (nm, e), flush=True) return None 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) def run_group(gname, Atr, ztr, ytr, yrtr, ptr, thr_tr, Ate, yte, pte, thr_te, eval_pts, jobs): """Cả hai tầng cho một nhóm huấn luyện. Song song hoá ở mức (fold x họ mô hình) chứ không ở mức nhóm: biến thể toàn quốc chỉ có MỘT nhóm, nếu song song theo nhóm thì 19 trong 20 lõi ngồi không. """ t0 = time.time() n = Atr.shape[0] z1, z2 = zoo_s1(n), zoo_s2(n) ys = np.unique(yrtr) folds = [f for f in np.array_split( np.random.default_rng(0).permutation(ys), NFOLD) if f.size] par = Parallel(n_jobs=jobs) # ---- tầng 1, ngoài mẫu: 4 fold x 8 họ, cộng 8 lần khớp đầy đủ masks = [np.isin(yrtr, f) for f in folds] masks = [m for m in masks if m.any() and not m.all()] tasks = [(j, i) for j in range(len(masks)) for i in range(len(z1))] outs = par(delayed(_fit_s1)(z1[i][1], z1[i][0], Atr[~masks[j]], ztr[~masks[j]], Atr[masks[j]], thr_tr[masks[j]]) for j, i in tasks) v_oof = np.zeros((len(z1), n), dtype="i1") for (j, i), o in zip(tasks, outs): v_oof[i, masks[j]] = o v_val = np.vstack(par(delayed(_fit_s1)(mk, nm, Atr, ztr, Ate, thr_te) for nm, mk in z1)) q_oof, q_val = p_of(v_oof), p_of(v_val) print(" [%s] tầng 1 xong: %d mẫu, %.0fs" % (gname, n, time.time() - t0), flush=True) # ---- cổng: ngưỡng p1 chốt RIÊNG TỪNG Ô trên trường ngoài mẫu của chính ô đó info = {} keep = np.zeros(n, bool) for p in eval_pts: mo, mv = ptr == p, pte == p if not mo.any() or not mv.any(): continue p1 = pick_p1(q_oof[mo], ytr[mo]) info[int(p)] = (mo, mv, p1) keep[mo] = q_oof[mo] >= p1 # ô lưới KHÔNG nằm trong tập đánh giá vẫn phải qua cổng để tầng 2 có đủ mẫu; # dùng ngưỡng trung vị của các ô đánh giá cùng nhóm if info: p1m = float(np.median([v[2] for v in info.values()])) rest = ~np.isin(ptr, list(info)) keep[rest] = q_oof[rest] >= p1m # ---- tầng 2 trên tập cổng giữ lại s2o = np.zeros(n); s2v = np.zeros(Ate.shape[0]) if keep.sum() >= 100 and 0 < ytr[keep].sum() < keep.sum(): t2 = [(j, i) for j in range(len(masks)) for i in range(len(z2))] o2 = par(delayed(_fit_s2)( z2[i][1], z2[i][0], Atr[(~masks[j]) & keep], ytr[(~masks[j]) & keep], Atr[masks[j] & keep]) for j, i in t2) acc = {} for (j, i), o in zip(t2, o2): if o is not None: acc.setdefault(j, []).append(o) for j, lst in acc.items(): s2o[masks[j] & keep] = np.mean(lst, axis=0) ov = par(delayed(_fit_s2)(mk, nm, Atr[keep], ytr[keep], Ate) for nm, mk in z2) ov = [o for o in ov if o is not None] if ov: s2v = np.mean(ov, axis=0) print(" [%s] tầng 2 xong: giữ %d/%d mẫu, tổng %.0fs" % (gname, int(keep.sum()), n, time.time() - t0), flush=True) res = {} for p, (mo, mv, p1) in info.items(): res[p] = dict(s_oof=np.where(q_oof[mo] >= p1, s2o[mo], -1.0), s_val=np.where(q_val[mv] >= p1, s2v[mv], -1.0), obs_oof=ytr[mo] > 0.5, obs_val=yte[mv] > 0.5, p1=float(p1)) return res def main(): ap = argparse.ArgumentParser() ap.add_argument("--variant", required=True, choices=["point", "region", "national"]) ap.add_argument("--window", required=True, choices=list(WINDOWS)) ap.add_argument("--fset", default="F3") ap.add_argument("--jobs", type=int, default=20) ap.add_argument("--target", default="so1", choices=["so1", "geq_k"], help="'so1' = có >= 1 ngày EHD trong cửa sổ (tiêu chí đang " "dùng); 'geq_k' = có >= k ngày, k chọn riêng cho mỗi ô " "lưới trên các năm huấn luyện để tần suất gần 1/3 nhất") ap.add_argument("--eval", default="sample", choices=["sample", "all"], help="'sample' = 42 điểm kiểm thử của GĐ1; " "'all' = toàn bộ 541 ô lưới đất, dùng cho GĐ2") args = ap.parse_args() os.makedirs(RES_DIR, exist_ok=True) store = D.Store() e0 = D.sample_points(store.grid, 3, 0) e1 = D.sample_points(store.grid, 3, 1) allp = np.arange(store.grid.reg.size) eval_pts = (allp if args.eval == "all" else np.unique(np.concatenate([e0, e1]))) train_pts = allp if args.variant != "point" else eval_pts lts = np.concatenate([C.WEEKS[i][1] for i in WINDOWS[args.window]]) print("[p1] %s / %s: %d điểm huấn luyện, %d điểm đánh giá" % (args.variant, args.window, train_pts.size, eval_pts.size), flush=True) t0 = time.time() Atr, Str, Ctr, Ytr, Ptr, _ = build_window(store, train_pts, C.TRAIN_YEARS, args.fset, lts) Ate, Ste, Cte, Yte, Pte, _ = build_window(store, eval_pts, C.VALID_YEARS, args.fset, lts) ok_t = np.isfinite(Atr).all(axis=1) & np.isfinite(Str) & np.isfinite(Ctr) ok_v = np.isfinite(Ate).all(axis=1) & np.isfinite(Ste) & np.isfinite(Cte) Atr, Str, Ctr, Ytr, Ptr = (Atr[ok_t], Str[ok_t], Ctr[ok_t], Ytr[ok_t], Ptr[ok_t]) Ate, Cte, Pte = Ate[ok_v], Cte[ok_v], Pte[ok_v] if args.target == "so1": ztr = Str.astype("f8") # hồi quy Sev cực đại tuần thr_tr = np.zeros(Atr.shape[0]); thr_te = np.zeros(Ate.shape[0]) ytr = (Ctr >= 1).astype("f8"); yte = (Cte >= 1).astype("f8") else: # k riêng cho từng ô lưới, chốt TRÊN CÁC NĂM HUẤN LUYỆN cho tần suất ~ 1/3 kmap = {} for p in np.unique(Ptr): c = Ctr[Ptr == p] ks = np.arange(1, int(max(c.max(), 1)) + 1) rates = np.array([(c >= kk).mean() for kk in ks]) kmap[int(p)] = float(ks[np.argmin(np.abs(rates - 1.0 / 3.0))]) thr_tr = np.array([kmap[int(p)] for p in Ptr]) thr_te = np.array([kmap.get(int(p), 1.0) for p in Pte]) ztr = Ctr.astype("f8") # hồi quy SỐ NGÀY EHD ytr = (Ctr >= thr_tr).astype("f8"); yte = (Cte >= thr_te).astype("f8") kk = np.array(list(kmap.values())) print("[p1] target >= k ngày: k trung vị %.0f, khoảng %d-%d" % (np.median(kk), kk.min(), kk.max()), flush=True) print("[p1] dữ liệu: train %d hàng, valid %d hàng, %.0fs" % (Atr.shape[0], Ate.shape[0], time.time() - t0), flush=True) reg = store.grid.reg if args.variant == "national": groups = [("national", np.ones(Atr.shape[0], bool), np.ones(Ate.shape[0], bool), eval_pts)] elif args.variant == "region": groups = [] for r in range(1, 8): gp = eval_pts[reg[eval_pts] == r] if gp.size: groups.append(("S%d" % r, reg[Ptr] == r, reg[Pte] == r, gp)) else: groups = [("p%d" % p, Ptr == p, Pte == p, np.array([p])) for p in eval_pts] # nhóm ít thì dồn lõi vào bên trong mỗi nhóm; nhóm nhiều thì song song theo nhóm if len(groups) >= args.jobs: got = Parallel(n_jobs=args.jobs, verbose=5)( delayed(run_group)(g, Atr[mt], ztr[mt], ytr[mt], Ytr[mt], Ptr[mt], thr_tr[mt], Ate[mv], yte[mv], Pte[mv], thr_te[mv], gp, 1) for g, mt, mv, gp in groups) else: got = [run_group(g, Atr[mt], ztr[mt], ytr[mt], Ytr[mt], Ptr[mt], thr_tr[mt], Ate[mv], yte[mv], Pte[mv], thr_te[mv], gp, args.jobs) for g, mt, mv, gp in groups] dump = {} for r in got: for p, d in r.items(): for k in ("s_oof", "s_val", "obs_oof", "obs_val"): dump["%d|%s" % (p, k)] = d[k] dump["%d|p1" % p] = np.array([d["p1"]]) dump["eval_seed0"] = e0 dump["eval_seed1"] = e1 tag = "%s_%s" % (args.variant, args.window) if args.target != "so1": tag = tag + "_" + args.target if args.eval == "all": tag = "full_" + tag out = os.path.join(RES_DIR, "%s.npz" % tag) np.savez_compressed(out, **dump) print("[done] %.0f phút -> %s" % ((time.time() - t0) / 60, out), flush=True) if __name__ == "__main__": main()