"""PHÉP THỬ QUYẾT ĐỊNH: trạng thái CHẬM tại init có nâng phân biệt EHD tuần 3-6 vượt trên trường DỰ BÁO đồng-thời của chính tuần đó không? Giả thuyết (khác trục mọi thứ đã làm): khả dự báo nội mùa của nắng nóng nằm ở trạng thái chậm biết từ lúc khởi tạo (độ ẩm đất, hoàn lưu quy mô lớn), KHÔNG ở quỹ đạo Tx. Mọi mô hình tuần hiện tại chỉ dùng trường của ĐÚNG tuần đích (đã bị sai số hỗn loạn chi phối ở hạn dài). Ở đây dùng trường tuần-1 (proxy cho init, gần quan trắc, dị thường dai dẳng) làm predictor CHO tuần đích 2..6. So trên 21 điểm (gộp), fit trên TRAIN_YEARS, chấm VALID_YEARS, pooled: CONCUR -- chỉ trường đồng-thời tuần đích (= thông tin các mô hình đang dùng). INIT -- chỉ trạng thái chậm tuần-1 (soil/hoàn lưu/EOF; BỎ Tx/severity nhanh). INIT_ALL -- toàn bộ trung bình tuần-1 (kể cả nhiệt ban đầu). BOTH -- đồng-thời + trạng thái chậm tuần-1 (câu hỏi chính: BOTH > CONCUR?). Đo AUC + SEDI tại điểm vận hành FBI=1. LogReg (co ngót) + LGBM. CĂN CHỈNH HÀNG: build_window trả hàng theo (init, point) chỉ phụ thuộc (years, points), KHÔNG phụ thuộc lts -> hai lần gọi với tuần-1 và tuần-đích căn chỉnh theo hàng. """ import json import os import sys import numpy as np from sklearn.linear_model import LogisticRegression import lightgbm as lgb sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import ehd_common as C from ai_data import Store from phase1_arch import build_window from ehd_metrics import auc_score, scores as sc_scores, contingency OUT = os.path.join(C.OUT, "AI_Results", "initstate") # Trạng thái CHẬM (dị thường dai dẳng) — cột z_ của F3, ở dạng trung bình tuần "m_*". # BỎ các biến nhiệt/severity nhanh (sv_mod, ev_mod, z_mx2t6, z_mn2t6, z_t2m, z_skt) # để cô lập tín hiệu trạng thái chậm khỏi chính dự báo nhiệt (vốn đã dùng, kém kỹ năng). SLOW = (["m_z_sm20", "m_z_sm100", "m_z_st20", "m_z_st100", "m_z_bowen", "m_z_slhf", "m_z_sshf", "m_z_tcw", "m_z_tcc", "m_z_cape", "m_z_msl", "m_z_u10", "m_z_v10", "m_z_uv10", "m_z_thk_1000_500", "m_z_thk_1000_850", "m_z_thk_850_700", "m_z_thk_925_500", "m_z_thk_925_700", "m_z_dt_850_500", "m_z_dt_925_850", "m_z_dt_700_500", "m_z_ivt", "m_z_ivt_u", "m_z_ivt_v", "m_z_w"] + ["m_h500%02d" % i for i in range(1, 11)] + ["m_t850%02d" % i for i in range(1, 11)]) def pick_points(store, per_region=3): reg = np.asarray(store.st["region"]).astype(int) pts = [] for r in range(1, 8): idx = np.nonzero(reg == r)[0] if idx.size: pts.extend(idx[np.linspace(0, idx.size - 1, per_region) .round().astype(int)].tolist()) return np.array(sorted(set(pts)), dtype=int) def prep(Xtr, Xte): med = np.nanmedian(Xtr, axis=0) Xtr = np.where(np.isnan(Xtr), med, Xtr) Xte = np.where(np.isnan(Xte), med, Xte) mu, sd = Xtr.mean(0), Xtr.std(0) sd[sd == 0] = 1.0 return (Xtr - mu) / sd, (Xte - mu) / sd def op(score, y): auc = auc_score(score, y) r = float(y.mean()) t = np.quantile(score, 1.0 - r) a, b, c, d = contingency(score >= t, y) s = sc_scores(a, b, c, d) return dict(AUC=auc, POD=s["POD"], FAR=s["FAR"], FBI=s["FBI"], CSI=s["CSI"], SEDI=s["SEDI"]) def fit_both(Xtr, ytr, Xte, yte): """AUC/SEDI của LogReg và LGBM; trả về (logreg_op, lgbm_op).""" m1 = LogisticRegression(C=0.3, max_iter=2000).fit(Xtr, ytr) p1 = m1.predict_proba(Xte)[:, 1] m2 = lgb.LGBMClassifier(n_estimators=300, learning_rate=0.05, num_leaves=31, min_child_samples=200, colsample_bytree=0.8, reg_lambda=5.0, n_jobs=8, verbose=-1).fit(Xtr, ytr) p2 = m2.predict_proba(Xte)[:, 1] return op(p1, yte), op(p2, yte) def main(): store = Store() pts = pick_points(store) os.makedirs(OUT, exist_ok=True) lt_w1 = C.WEEKS[0][1] # lead days của tuần 1 print("21-point testbed; init-proxy = tuần-1; đích = tuần 2..6", flush=True) res = {} print("\n%-4s %-9s | LogReg AUC/SEDI | LGBM AUC/SEDI" % ("W", "featset"), flush=True) for wi in range(1, 6): # W2..W6 wname, lt_w = C.WEEKS[wi] # đồng-thời tuần đích (train/valid) Ac_tr, _, N_tr, _, _, names = build_window( store, pts, C.TRAIN_YEARS, "F3", lt_w) Ac_te, _, N_te, _, _, _ = build_window( store, pts, C.VALID_YEARS, "F3", lt_w) # trạng thái tuần-1 (cùng hàng: cùng years+points) Ai_tr, _, _, _, _, _ = build_window(store, pts, C.TRAIN_YEARS, "F3", lt_w1) Ai_te, _, _, _, _, _ = build_window(store, pts, C.VALID_YEARS, "F3", lt_w1) islow = [names.index(c) for c in SLOW if c in names] nmean = sum(1 for c in names if c.startswith("m_")) # số cột trung bình imean = [i for i, c in enumerate(names) if c.startswith("m_")] Xi_tr, Xi_te = Ai_tr[:, islow], Ai_te[:, islow] # init slow Xia_tr, Xia_te = Ai_tr[:, imean], Ai_te[:, imean] # init all-mean ytr = (N_tr >= 1).astype(int) yte = (N_te >= 1).astype(int) feats = { "CONCUR": (Ac_tr, Ac_te), "INIT": (Xi_tr, Xi_te), "INIT_ALL": (Xia_tr, Xia_te), "BOTH": (np.column_stack([Ac_tr, Xi_tr]), np.column_stack([Ac_te, Xi_te])), } res[wname] = {"base": float(yte.mean()), "n_slow": len(islow)} for tag, (Xtr0, Xte0) in feats.items(): Xtr, Xte = prep(Xtr0, Xte0) lr, gb = fit_both(Xtr, ytr, Xte, yte) res[wname][tag] = {"logreg": lr, "lgbm": gb} print("%-4s %-9s | %.3f / %.3f | %.3f / %.3f" % (wname, tag, lr["AUC"], lr["SEDI"], gb["AUC"], gb["SEDI"]), flush=True) # kết luận nhanh cho tuần này dA_lr = res[wname]["BOTH"]["logreg"]["AUC"] - res[wname]["CONCUR"]["logreg"]["AUC"] dA_gb = res[wname]["BOTH"]["lgbm"]["AUC"] - res[wname]["CONCUR"]["lgbm"]["AUC"] print(" -> ΔAUC(BOTH-CONCUR) LogReg %+.3f LGBM %+.3f | INIT-alone AUC " "LR %.3f GB %.3f\n" % (dA_lr, dA_gb, res[wname]["INIT"]["logreg"]["AUC"], res[wname]["INIT"]["lgbm"]["AUC"]), flush=True) with open(os.path.join(OUT, "initstate_probe.json"), "w") as f: json.dump(res, f, indent=1, default=float) print("[done] ->", os.path.join(OUT, "initstate_probe.json"), flush=True) if __name__ == "__main__": main()