#!/usr/bin/env python3 """Thành viên NEO vào dự báo của chính mô hình EC (phương án b). Ý tưởng lấy từ Van Straaten et al. (2023, AIES 2): mạng KHÔNG dự báo lại từ đầu mà chỉ học phần HIỆU CHỈNH đặt trên nền tiên nghiệm của mô hình động lực. Ở đó log(p_raw) được cộng thẳng vào logit ngay trước softmax, nên đầu ra của mạng đọc được như logarit của một HỆ SỐ NHÂN phụ thuộc thời tiết lên xác suất của ECMWF. Bảng C1 của họ cho thấy phải làm đúng như vậy: nhét p_raw vào làm một biến đầu vào bình thường chỉ đạt BS 0.190, trong khi cộng vào logit đạt 0.168. KHÁC BIỆT BẮT BUỘC SO VỚI BÀI BÁO. Họ có 11 thành phần tổ hợp nên p_raw đếm được trực tiếp. Dữ liệu ở đây chỉ có control forecast, không có chiều thành phần, nên p_raw phải suy ra bằng cách hiệu chỉnh xác suất từ dự báo tất định. Việc hiệu chỉnh chỉ dùng đầu ra của chính EC tại thời điểm hiệu lực (không dùng quan trắc, theo quy ước vận hành ở AI17) và được ước lượng NGOÀI MẪU trên các năm huấn luyện. """ from __future__ import annotations import numpy as np EPS = 1e-4 def logit(p): p = np.clip(np.asarray(p, dtype="f8"), EPS, 1.0 - EPS) return np.log(p / (1.0 - p)) def sigmoid(x): return 1.0 / (1.0 + np.exp(-np.clip(np.asarray(x, dtype="f8"), -30, 30))) # ------------------------------------------------------------------- p_raw class RawProb: """Xác suất tiên nghiệm p_raw của EC cho yEHDE tuần, hiệu chỉnh từ EC-only. Đầu vào chỉ gồm các thống kê tuần của chính dự báo EC: mức nghiêm trọng trung bình và cực đại trong tuần, số ngày mô hình tự cho là nóng, và tần suất sự kiện theo mô hình. Dung lượng mô hình cố tình để rất thấp (một hồi quy logistic bốn hệ số) vì đây phải là TIÊN NGHIỆM chứ không phải một mô hình dự báo thứ hai -- nếu p_raw đã học nhiều thì phần "hiệu chỉnh" mất hết ý nghĩa. """ COLS = ("m_sv_mod", "x_sv_mod", "n_hotday_mod", "m_ev_mod") def __init__(self, names): self.idx = [names.index(c) for c in self.COLS if c in names] def fit(self, X, y): from sklearn.linear_model import LogisticRegression from sklearn.pipeline import make_pipeline from sklearn.preprocessing import StandardScaler y = np.asarray(y).astype(int) # không đủ điều kiện khớp thì lùi về tần suất khí hậu -- vẫn là một tiên # nghiệm hợp lệ, chỉ là không có thông tin thời tiết trong đó self.base_ = float(np.clip(y.mean(), EPS, 1.0 - EPS)) self.m_ = None if len(self.idx) and np.unique(y).size == 2: self.m_ = make_pipeline(StandardScaler(), LogisticRegression(C=1.0, max_iter=500)) self.m_.fit(X[:, self.idx], y) return self def predict_proba(self, X): if self.m_ is None: return np.full(X.shape[0], self.base_) return self.m_.predict_proba(X[:, self.idx])[:, 1] # ------------------------------------------------- mạng có offset (bài báo) class OffsetNet: """Phân loại yEHDE với logit neo: p = sigmoid(offset + net(X)). hidden = () cho đúng dòng "0 hidden layers" của bảng C1, tức hồi quy logistic có offset; hidden = (4,) là chính kiến trúc thắng cuộc của bài báo. """ def __init__(self, hidden=(4,), act="elu", dropout=0.0, epochs=200, lr=1.4e-3, wd=1e-4, batch=32, seed=0): self.hidden, self.act, self.dropout = hidden, act, dropout self.epochs, self.lr, self.wd, self.batch, self.seed = \ epochs, lr, wd, batch, seed def _act(self): import torch.nn as nn return {"relu": nn.ReLU, "tanh": nn.Tanh, "gelu": nn.GELU, "elu": nn.ELU}[self.act]() def fit(self, X, y, offset): import torch import torch.nn as nn torch.set_num_threads(1) torch.manual_seed(self.seed) X = np.asarray(X, dtype="f4") self.mu_, self.sd_ = X.mean(0), X.std(0) self.sd_[self.sd_ < 1e-6] = 1.0 layers, prev = [], X.shape[1] for h in self.hidden: layers += [nn.Linear(prev, h), self._act()] if self.dropout > 0: layers.append(nn.Dropout(self.dropout)) prev = h last = nn.Linear(prev, 1) nn.init.zeros_(last.weight); nn.init.zeros_(last.bias) # khởi động ở "không sửa gì" layers.append(last) net = nn.Sequential(*layers) opt = torch.optim.Adam(net.parameters(), lr=self.lr, weight_decay=self.wd) lossf = nn.BCEWithLogitsLoss() Xt = torch.from_numpy((X - self.mu_) / self.sd_) ot = torch.from_numpy(np.asarray(offset, dtype="f4")).unsqueeze(1) yt = torch.from_numpy(np.asarray(y, dtype="f4").ravel()).unsqueeze(1) n = X.shape[0] rng = np.random.default_rng(self.seed) net.train() for _ in range(self.epochs): order = rng.permutation(n) for a in range(0, n, self.batch): i = torch.from_numpy(order[a:a + self.batch]) opt.zero_grad() lossf(net(Xt[i]) + ot[i], yt[i]).backward() opt.step() net.eval() self.net_ = net return self def predict_proba(self, X, offset): import torch X = np.asarray(X, dtype="f4") with torch.no_grad(): z = self.net_(torch.from_numpy((X - self.mu_) / self.sd_)).numpy().ravel() return sigmoid(z + np.asarray(offset, dtype="f8")) class OffsetLGBM: """Cây tăng cường với init_score = offset -- cơ chế neo tương đương OffsetNet.""" def __init__(self, **kw): self.kw = kw def fit(self, X, y, offset): import lightgbm as lgb self.m_ = lgb.LGBMClassifier(n_jobs=1, verbose=-1, **self.kw) self.m_.fit(X, np.asarray(y).astype(int), init_score=np.asarray(offset, dtype="f8")) return self def predict_proba(self, X, offset): raw = self.m_.predict(X, raw_score=True) return sigmoid(raw + np.asarray(offset, dtype="f8")) def offset_zoo(): """Tám họ có thể neo THẬT (offset), phủ ba lớp mô hình khác nhau.""" return [ ("Off_Logit", "hồi quy logistic có offset (bảng C1: '0 hidden layers')", lambda: OffsetNet(hidden=(), epochs=200, seed=0)), ("Off_ANN_4", "ANN 1 lớp ẩn 4 nút, elu — đúng kiến trúc thắng cuộc của bài báo", lambda: OffsetNet(hidden=(4,), act="elu", epochs=200, seed=1)), ("Off_ANN_8_tanh", "ANN 1 lớp ẩn 8 nút, tanh", lambda: OffsetNet(hidden=(8,), act="tanh", epochs=200, seed=2)), ("Off_ANN_16", "ANN 1 lớp ẩn 16 nút, relu", lambda: OffsetNet(hidden=(16,), act="relu", epochs=200, seed=3)), ("Off_ANN_32_16", "ANN 2 lớp 32-16, relu, dropout 0.1", lambda: OffsetNet(hidden=(32, 16), act="relu", dropout=0.1, epochs=200, seed=4)), ("Off_LGBM_gbdt", "cây tăng cường độ dốc, init_score = offset", lambda: OffsetLGBM(n_estimators=200, num_leaves=15, learning_rate=0.05, min_child_samples=60, colsample_bytree=0.8, reg_lambda=5.0)), # LightGBM CẤM init_score với boosting_type="rf" (rf.hpp dòng 49), nên # biến thể rừng ngẫu nhiên không neo được; thay bằng một cấu hình gbdt sâu # hơn để vẫn có hai mức độ phức tạp cây khác nhau trong nhóm neo. ("Off_LGBM_deep", "cây tăng cường lá sâu (63 lá), init_score = offset", lambda: OffsetLGBM(n_estimators=150, num_leaves=63, learning_rate=0.08, min_child_samples=30, colsample_bytree=0.8, subsample=0.8, subsample_freq=1, reg_lambda=2.0)), ("Off_LGBM_slow", "cây tăng cường học chậm, 400 cây, lá nông", lambda: OffsetLGBM(n_estimators=400, num_leaves=7, learning_rate=0.02, min_child_samples=80, colsample_bytree=0.7, reg_lambda=10.0)), ] # ------------------------------------------------------- neo bằng Sev của EC def resid_families(zoo, names): """Mười họ hồi quy tốt nhất của lần chạy trước, dùng cho biến thể neo-Sev.""" keep = ["BayesRidge", "ElasticNet", "Huber", "Ridge", "PLS8", "RandomForest", "ExtraTrees", "LGBM_rf", "LGBM_dart", "LGBM_gbdt"] by = {z[0]: z for z in zoo} return [by[k] for k in keep if k in by]