from __future__ import annotations """ Memory-safe, drop-in refined version of prediction.py - Streams very large log files (no giant raw DataFrames) - Computes observed flags without large Python sets - Zero-copy feature audit (mask features on the fly) You can run this directly. Key differences from the original are marked with "STREAMING" and "ZERO-COPY" comments. """ import glob, os, re from collections import defaultdict from dataclasses import dataclass from datetime import datetime from pathlib import Path from typing import List, Tuple, Dict, Optional import numpy as np import pandas as pd import matplotlib.pyplot as plt import torch import torch.nn as nn from torch.utils.data import Dataset, DataLoader # ---------------------------- # Utils # ---------------------------- def set_seed(seed=42): import random random.seed(seed) np.random.seed(seed) torch.manual_seed(seed) torch.cuda.manual_seed_all(seed) set_seed(42) def ensure_dir(p: Path): p.mkdir(parents=True, exist_ok=True) return p def _run_lengths(mask: np.ndarray): """Return number of runs and their lengths (consecutive True segments).""" if mask.sum() == 0: return 0, np.array([], dtype=int) m = mask.astype(bool) dm = np.diff(np.concatenate(([False], m, [False])).astype(np.int8)) starts = np.where(dm == 1)[0] ends = np.where(dm == -1)[0] lens = ends - starts return len(lens), lens def report_underpredictions( timestamps: pd.Series, actual: np.ndarray, pred: np.ndarray, out_csv_path ): """ Save a per-minute CSV of under-predictions (pred < actual) and return a dict of summary metrics. """ actual = actual.astype("float64") pred = pred.astype("float64") finite = np.isfinite(actual) & np.isfinite(pred) if finite.sum() == 0: return { "n_pred_points": 0, "n_under": 0, "under_rate": np.nan, "mae_under": np.nan, "sum_under": np.nan, "max_under": np.nan, "p95_under": np.nan, "n_runs": 0, "avg_run": np.nan, "max_run": 0 } under = finite & (pred < actual) gap = (actual - pred) # positive when under rel_gap = np.divide(gap, np.maximum(actual, 1e-9)) # summary metrics n_pred = int(finite.sum()) n_under = int(under.sum()) under_rate = n_under / n_pred if n_under > 0: mae_under = float(gap[under].mean()) sum_under = float(gap[under].sum()) max_under = float(gap[under].max()) p95_under = float(np.percentile(gap[under], 95)) runs, lens = _run_lengths(under) avg_run = float(lens.mean()) if runs > 0 else np.nan max_run = int(lens.max()) if runs > 0 else 0 else: mae_under = sum_under = max_under = p95_under = np.nan runs = 0; avg_run = np.nan; max_run = 0 # write detailed CSV (only the under-predicted rows) df = pd.DataFrame({ "timestamp": np.asarray(timestamps)[finite], "actual": actual[finite], "pred": pred[finite], "under": (pred[finite] < actual[finite]), "gap": gap[finite], "rel_gap": rel_gap[finite], }) df_under = df[df["under"]].drop(columns=["under"]) df_under.to_csv(out_csv_path, index=False) return { "n_pred_points": n_pred, "n_under": n_under, "under_rate": under_rate, "mae_under": mae_under, "sum_under": sum_under, "max_under": max_under, "p95_under": p95_under, "n_runs": runs, "avg_run": avg_run, "max_run": max_run, } # ---------------------------- # Config # ---------------------------- # NOTE: fixed missing comma in DATE_STRS DATE_STRS = ["20250828", "20250829"] # train/eval days handled independently LOG_DIR = "artifacts" RESULT_DIR = Path("results") SEQ_LEN = 45 # history length (minutes) PRED_LEN = 10 # forecast horizon (minutes) BATCH = 256 EPOCHS_FIRST = 4 # per-target epochs on the first day (you can raise) EPOCHS_NEXT = 2 # per-target epochs on subsequent days LR = 1e-3 ADD_TIME_OF_DAY = True # add sin/cos time-of-day features # === Feature selection === ENABLE_FEATURE_SELECTION = True AUDIT_DAY = DATE_STRS[0] FS_EPOCHS = max(2, min(3, EPOCHS_FIRST)) # quick audit training FS_IMPROVE_THRESH = 0.01 # drop if masking improves val loss by >=1% FS_NOISE_BAND = 0.002 # ignore micro-changes (<0.2%) DEVICE = torch.device("cuda" if torch.cuda.is_available() else "cpu") # Model targets TARGETS = ["num_ues", "traffic_load"] # Columns we aim to have per minute (some may be zero if not present) FOCUS_COLS = [ "num_ues","traffic_load","dl_prbs","ul_prbs","dl_sched_rbs","ul_sched_rbs", "wb_cqi","dl_mcs1","ul_mcs1","pusch_snr","pucch_snr", "dl_bler","ul_bler","bsr","phr", "rrc_total","rrc_unique_rnti","rrc_setup_cpl","rrc_release", "rrc_recfg","rrc_recfg_cpl","rrc_resume","rrc_paging" ] # ---------------------------- # Regex / Parsing # ---------------------------- LOG_PATTERN_FULL = re.compile( r"""\[ (?P.*?)\] \s+UE\s+RNTI=(?P\d+) \s+DL_TBS=(?P\d+) \s+UL_TBS=(?P\d+) \s+DL_SDUs=(?P\d+) \s+UL_SDUs=(?P\d+) \s+DL_sched_RBs=(?P\d+) \s+UL_sched_RBs=(?P\d+) \s+DL_PRBs=(?P\d+) \s+UL_PRBs=(?P\d+) \s+DL_BLER=(?P[\d.]+) \s+UL_BLER=(?P[\d.]+) \s+PUSCH_SNR=(?P[\d.]+) \s+PUCCH_SNR=(?P[\d.]+) \s+WB_CQI=(?P\d+) \s+DL_MCS1=(?P\d+) \s+UL_MCS1=(?P\d+) \s+DL_HARQ=(?P\d+) \s+UL_HARQ=(?P\d+) \s+BSR=(?P\d+) \s+PHR=(?P\d+) """, re.VERBOSE, ) RRC_PATTERNS: List[re.Pattern] = [ re.compile(r"\[(?P[^]]+)\].*RNTI=(?P[0-9a-fA-F]+).*RRC[_ ]?MSG=(?P[A-Za-z_]+)", re.IGNORECASE), re.compile(r"\[(?P[^]]+)\].*rrc[=:](?P[A-Za-z_]+).*rnti[=:](?P[0-9a-fA-F]+)", re.IGNORECASE), re.compile(r"\[(?P[^]]+)\].*(?PRRC[A-Za-z]+).*RNTI=(?P[0-9a-fA-F]+)", re.IGNORECASE), re.compile(r"\|\s*(?P\d{4}-\d{2}-\d{2}T\d{2}:\d{2}:\d{2}(?:\.\d+)?Z?).*?RNTI\s*(?P[0-9a-fA-F]+)", re.IGNORECASE), ] # ---------------------------- # STREAMING: Per-minute aggregators without raw-frame concatenation # ---------------------------- C_RNTI_MIN, C_RNTI_MAX = 0x003D, 0xFFFC @dataclass class _MACAgg: dl_sdus_sum: float = 0.0; ul_sdus_sum: float = 0.0 dl_sched_rbs_sum: float = 0.0; ul_sched_rbs_sum: float = 0.0 dl_prbs_sum: float = 0.0; ul_prbs_sum: float = 0.0 dl_bler_sum: float = 0.0; ul_bler_sum: float = 0.0 pusch_snr_sum: float = 0.0; pucch_snr_sum: float = 0.0 wb_cqi_sum: float = 0.0; dl_mcs1_sum: float = 0.0; ul_mcs1_sum: float = 0.0 dl_harq_sum: float = 0.0; ul_harq_sum: float = 0.0 bsr_sum: float = 0.0; phr_sum: float = 0.0 count: int = 0 dl_bytes: int = 0; ul_bytes: int = 0 ues: Optional[set] = None def _floor_min(ts: datetime) -> pd.Timestamp: return pd.Timestamp(ts).floor("min") def stream_mac_minutes(files: List[str]) -> pd.DataFrame: if not files: return pd.DataFrame() acc: Dict[pd.Timestamp, _MACAgg] = {} for path in files: with open(path, "r", errors="ignore") as f: for line in f: m = LOG_PATTERN_FULL.search(line) if not m: continue d = m.groupdict() try: ts = datetime.fromisoformat(d["ts"]) # may raise except Exception: continue t = _floor_min(ts) A = acc.get(t) if A is None: A = _MACAgg(ues=set()) acc[t] = A # unique UE when real RNTI with any traffic try: rnti_val = int(d["rnti"]) # MAC logs use decimal here except Exception: rnti_val = -1 dl_tbs = int(d["dl_tbs"]); ul_tbs = int(d["ul_tbs"]) if (C_RNTI_MIN <= rnti_val <= C_RNTI_MAX) and (dl_tbs > 0 or ul_tbs > 0): A.ues.add(rnti_val) # sums & bytes A.dl_bytes += dl_tbs; A.ul_bytes += ul_tbs A.dl_sdus_sum += int(d["dl_sdus"]); A.ul_sdus_sum += int(d["ul_sdus"]) A.dl_sched_rbs_sum += int(d["dl_sched_rbs"]); A.ul_sched_rbs_sum += int(d["ul_sched_rbs"]) A.dl_prbs_sum += int(d["dl_prbs"]); A.ul_prbs_sum += int(d["ul_prbs"]) A.dl_bler_sum += float(d["dl_bler"]); A.ul_bler_sum += float(d["ul_bler"]) A.pusch_snr_sum += float(d["pusch_snr"]); A.pucch_snr_sum += float(d["pucch_snr"]) A.wb_cqi_sum += int(d["wb_cqi"]); A.dl_mcs1_sum += int(d["dl_mcs1"]); A.ul_mcs1_sum += int(d["ul_mcs1"]) A.dl_harq_sum += int(d["dl_harq"]); A.ul_harq_sum += int(d["ul_harq"]) A.bsr_sum += int(d["bsr"]); A.phr_sum += int(d["phr"]) A.count += 1 if not acc: return pd.DataFrame() rows = [] for t in sorted(acc.keys()): A = acc[t] c = max(1, A.count) row = { "timestamp": t, "num_ues": len(A.ues) if A.ues is not None else 0, "dl_bytes": A.dl_bytes, "ul_bytes": A.ul_bytes, "dl_sdus": A.dl_sdus_sum / c, "ul_sdus": A.ul_sdus_sum / c, "dl_sched_rbs": A.dl_sched_rbs_sum / c, "ul_sched_rbs": A.ul_sched_rbs_sum / c, "dl_prbs": A.dl_prbs_sum / c, "ul_prbs": A.ul_prbs_sum / c, "dl_bler": A.dl_bler_sum / c, "ul_bler": A.ul_bler_sum / c, "pusch_snr": A.pusch_snr_sum / c, "pucch_snr": A.pucch_snr_sum / c, "wb_cqi": A.wb_cqi_sum / c, "dl_mcs1": A.dl_mcs1_sum / c, "ul_mcs1": A.ul_mcs1_sum / c, "dl_harq": A.dl_harq_sum / c, "ul_harq": A.ul_harq_sum / c, "bsr": A.bsr_sum / c, "phr": A.phr_sum / c, "traffic_load": A.dl_bytes + A.ul_bytes, "mac_present": True, } rows.append(row) out = pd.DataFrame(rows) out.sort_values("timestamp", inplace=True) return out RRC_FAM_COLS = ["rrc_conn_req","rrc_setup_cpl","rrc_recfg","rrc_recfg_cpl","rrc_reestab", "rrc_release","rrc_resume","rrc_suspend","rrc_paging","rrc_other"] def canonical_rrc_family(msg: str | None) -> str: if not msg: return "rrc_other" m = msg.lower() if "connectionrequest" in m or "connrequest" in m or "setupreq" in m: return "rrc_conn_req" if "connectionsetupcomplete" in m or "setupcomplete" in m: return "rrc_setup_cpl" if "reconfigurationcomplete" in m or "reconfigcomplete" in m: return "rrc_recfg_cpl" if "reconfiguration" in m or "reconfig" in m: return "rrc_recfg" if "reestablish" in m: return "rrc_reestab" if "release" in m: return "rrc_release" if "pag" in m: return "rrc_paging" if "resume" in m: return "rrc_resume" if "suspend" in m: return "rrc_suspend" return "rrc_other" def stream_rrc_minutes(files: List[str]) -> pd.DataFrame: if not files: return pd.DataFrame() total = defaultdict(int) unique = defaultdict(set) fam = {c: defaultdict(int) for c in RRC_FAM_COLS} for path in files: with open(path, "r", errors="ignore") as f: for line in f: ts = rnti = msg = None for pat in RRC_PATTERNS: m = pat.search(line) if m: gd = m.groupdict(); ts = gd.get("ts"); rnti = gd.get("rnti"); msg = gd.get("msg"); break if not ts: continue ts = ts.rstrip("Z") try: t = _floor_min(datetime.fromisoformat(ts)) except ValueError: continue total[t] += 1 fam_key = canonical_rrc_family(msg) fam[fam_key][t] += 1 # unique per minute if we can parse rnti rid = None if rnti is not None: try: rid = int(rnti, 16) if any(c.isalpha() for c in rnti) else int(rnti, 10) except ValueError: pass if rid is not None: unique[t].add(rid) if not total: return pd.DataFrame() minutes = sorted(total.keys()) data = { "timestamp": minutes, "rrc_total": [total[t] for t in minutes], "rrc_unique_rnti": [len(unique[t]) if t in unique else 0 for t in minutes], "rrc_present": [True] * len(minutes), } for c in RRC_FAM_COLS: data[c] = [fam[c][t] if t in fam[c] else 0 for t in minutes] out = pd.DataFrame(data) return out # ---------------------------- # Per-minute features (compat helpers kept for reference) # ---------------------------- def full_minute_index(start, end): return pd.date_range(start.floor("min"), end.ceil("min"), freq="1min") def add_time_of_day_feats(df: pd.DataFrame) -> pd.DataFrame: if not ADD_TIME_OF_DAY: return df out = df.copy() m = out["timestamp"].dt.hour * 60 + out["timestamp"].dt.minute out["tod_sin"] = np.sin(2*np.pi*m/1440.0).astype("float32") out["tod_cos"] = np.cos(2*np.pi*m/1440.0).astype("float32") return out # ---------------------------- # STREAMING: Memory-safe build_day_features # ---------------------------- def build_day_features(day: str) -> pd.DataFrame: mac_files = glob.glob(os.path.join(LOG_DIR, f"mac_stats_{day}.log")) \ + glob.glob(os.path.join(LOG_DIR, f"mac_stats_{day}_*.log")) rrc_files = glob.glob(os.path.join(LOG_DIR, f"rrc_messages_{day}.log")) \ + glob.glob(os.path.join(LOG_DIR, f"rrc_messages_{day}_*.log")) mac_feat = stream_mac_minutes(mac_files) rrc_feat = stream_rrc_minutes(rrc_files) if mac_feat.empty and rrc_feat.empty: raise SystemExit(f"[{day}] no parsable lines found") # build the full per-minute index from whichever side exists frames = [df for df in (mac_feat, rrc_feat) if not df.empty] t_start = min(df["timestamp"].min() for df in frames) t_end = max(df["timestamp"].max() for df in frames) idx = pd.date_range(t_start.floor("min"), t_end.ceil("min"), freq="1min") # --- MAC side: reindex or fabricate skeleton if empty if not mac_feat.empty: mac_feat = (mac_feat.set_index("timestamp") .reindex(idx).fillna(0.0) .rename_axis("timestamp").reset_index()) mac_feat["mac_present"] = True else: mac_feat = pd.DataFrame({"timestamp": idx}) mac_feat["mac_present"] = False # no MAC rows this minute # --- RRC side: reindex or fabricate skeleton if empty if not rrc_feat.empty: rrc_feat = (rrc_feat.set_index("timestamp") .reindex(idx).fillna(0.0) .rename_axis("timestamp").reset_index()) rrc_feat["rrc_present"] = True else: rrc_feat = pd.DataFrame({"timestamp": idx}) # minimal set the rest of the code expects: rrc_feat["rrc_total"] = 0.0 rrc_feat["rrc_unique_rnti"] = 0.0 for c in RRC_FAM_COLS: rrc_feat[c] = 0.0 rrc_feat["rrc_present"] = False # merge & flags df = (mac_feat.merge(rrc_feat, on="timestamp", how="left") .fillna(0.0).reset_index(drop=True)) present_cols = [] if "mac_present" in df: present_cols.append(df["mac_present"].astype(bool)) if "rrc_present" in df: present_cols.append(df["rrc_present"].astype(bool)) df["observed"] = False if not present_cols else np.logical_or.reduce(present_cols) # time-of-day & compact dtypes if ADD_TIME_OF_DAY: m = df["timestamp"].dt.hour * 60 + df["timestamp"].dt.minute df["tod_sin"] = np.sin(2*np.pi*m/1440.0).astype("float32") df["tod_cos"] = np.cos(2*np.pi*m/1440.0).astype("float32") float64_cols = df.select_dtypes(include=["float64"]).columns if len(float64_cols): df[float64_cols] = df[float64_cols].astype("float32") return df # ---------------------------- # Sequence dataset (streaming) # ---------------------------- def valid_t_indices(observed: np.ndarray, seq_len: int, pred_len: int) -> np.ndarray: """Return t indices where [t-seq_len, t+pred_len-1] are all observed=True.""" obs = observed.astype(np.uint8) w = seq_len + pred_len s = pd.Series(obs).rolling(w, min_periods=w).sum().to_numpy() idx = [] for i in range(w - 1, len(obs)): if s[i] == w: t = i - pred_len + 1 if t >= seq_len: idx.append(t) return np.array(idx, dtype=np.int64) class TSWindowDataset(Dataset): def __init__(self, values_norm: np.ndarray, t_indices: np.ndarray, seq_len: int, pred_len: int, tgt_idx: int): self.V = values_norm self.t_idx = t_indices self.seq = seq_len self.h = pred_len self.tgt = tgt_idx def __len__(self): return len(self.t_idx) def __getitem__(self, i): t = self.t_idx[i] x = self.V[t - self.seq : t, :] y = self.V[t : t + self.h, self.tgt] return torch.from_numpy(x), torch.from_numpy(y) # ---------------------------- # LSTM model # ---------------------------- class Seq2PointLSTM(nn.Module): def __init__(self, input_size, hidden=128, layers=1, pred_len=PRED_LEN, dropout=0.0): super().__init__() self.lstm = nn.LSTM(input_size, hidden, layers, batch_first=True, dropout=dropout if layers>1 else 0.0) self.fc = nn.Linear(hidden, pred_len) def forward(self, x): enc, _ = self.lstm(x) return self.fc(enc[:, -1, :]) # ---------------------------- # Train / Forecast per day # ---------------------------- def prepare_day_matrix(day_df: pd.DataFrame, feat_cols_override: List[str] | None = None) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, list]: """Return values_norm, mu, sigma, observed, feature_names in a fixed order.""" base_cols = FOCUS_COLS + (["tod_sin","tod_cos"] if ADD_TIME_OF_DAY else []) feat_cols = feat_cols_override if feat_cols_override is not None else base_cols # ensure presence; if missing, create zeros df = day_df.copy() for c in feat_cols: if c not in df.columns: df[c] = 0.0 Xdf = df[feat_cols].astype("float32").copy() observed = df["observed"].to_numpy() obs_mask = observed.astype(bool) mu = Xdf[obs_mask].mean(0).to_numpy(dtype=np.float32) sigma = Xdf[obs_mask].std(0).to_numpy(dtype=np.float32) + 1e-8 values_norm = ((Xdf.to_numpy(dtype=np.float32) - mu) / sigma).astype(np.float32) return values_norm, mu, sigma, observed.astype(np.bool_), feat_cols def build_dataloaders(values_norm: np.ndarray, observed: np.ndarray, feat_cols: list, target: str): tgt_idx = feat_cols.index(target) t_idx_all = valid_t_indices(observed, SEQ_LEN, PRED_LEN) if len(t_idx_all) < 100: raise RuntimeError(f"Not enough valid windows for {target} on this day.") split = int(0.9 * len(t_idx_all)) t_tr, t_va = t_idx_all[:split], t_idx_all[split:] ds_tr = TSWindowDataset(values_norm, t_tr, SEQ_LEN, PRED_LEN, tgt_idx) ds_va = TSWindowDataset(values_norm, t_va, SEQ_LEN, PRED_LEN, tgt_idx) dl_tr = DataLoader(ds_tr, batch_size=BATCH, shuffle=True, num_workers=0, pin_memory=True) dl_va = DataLoader(ds_va, batch_size=BATCH, shuffle=False, num_workers=0, pin_memory=True) return dl_tr, dl_va, tgt_idx def eval_val_loss(model: nn.Module, dl_va: DataLoader) -> float: crit = nn.MSELoss() was_training = model.training # remember mode model.eval() with torch.no_grad(): tot = 0.0; n = 0 for xb, yb in dl_va: xb = xb.to(DEVICE); yb = yb.to(DEVICE) pred = model(xb) tot += crit(pred, yb).item() * xb.size(0); n += xb.size(0) if was_training: model.train(True) # restore mode return tot / max(1, n) def train_one_day(model: nn.Module, dl_tr: DataLoader, dl_va: DataLoader, epochs: int) -> float: opt = torch.optim.Adam(model.parameters(), lr=LR) crit = nn.MSELoss() best_va = float("inf") for _ in range(epochs): model.train(True) for xb, yb in dl_tr: xb = xb.to(DEVICE); yb = yb.to(DEVICE) opt.zero_grad() pred = model(xb) loss = crit(pred, yb) loss.backward() opt.step() va = eval_val_loss(model, dl_va) best_va = min(best_va, va) return best_va def rolling_forecast_day(model: nn.Module, values_norm: np.ndarray, observed: np.ndarray, feat_cols: list, target: str, mu: np.ndarray, sigma: np.ndarray) -> np.ndarray: tgt_idx = feat_cols.index(target) N = values_norm.shape[0] merged = np.full(N, np.nan, dtype=np.float32) counts = np.zeros(N, dtype=np.float32) t_idx_all = valid_t_indices(observed, SEQ_LEN, PRED_LEN) model.eval() with torch.no_grad(): for t in t_idx_all: x = torch.from_numpy(values_norm[t-SEQ_LEN:t, :]).unsqueeze(0).to(DEVICE) yhat = model(x).squeeze(0).cpu().numpy() # normalized yhat_dn = yhat * sigma[tgt_idx] + mu[tgt_idx] for j, v in enumerate(yhat_dn): i = t + j if i < N: if np.isnan(merged[i]): merged[i] = 0.0 merged[i] += v; counts[i] += 1.0 with np.errstate(invalid="ignore", divide="ignore"): merged = np.divide(merged, counts, out=merged, where=counts > 0) return merged # ---------------------------- # ZERO-COPY: Feature selection (mask-on-the-fly) # ---------------------------- def eval_val_loss_masked(model: nn.Module, dl_va: DataLoader, col_to_zero: Optional[int]) -> float: crit = nn.MSELoss() was_training = model.training model.eval() tot = 0.0; n = 0 with torch.no_grad(): for xb, yb in dl_va: if col_to_zero is not None: xb = xb.clone() xb[:, :, col_to_zero] = 0.0 xb = xb.to(DEVICE); yb = yb.to(DEVICE) pred = model(xb) tot += crit(pred, yb).item() * xb.size(0); n += xb.size(0) if was_training: model.train(True) return tot / max(1, n) def feature_audit_and_select(day_df: pd.DataFrame, target: str, base_cols: List[str]) -> Tuple[List[str], Dict[str, float]]: """ Memory-safe audit: - Train once on AUDIT_DAY (fast) - For each feature, evaluate validation while zeroing that column on-the-fly (no full matrix copies) Returns (kept_features, deltas) where deltas[c] = (masked - base)/base. """ print(f"[{AUDIT_DAY}] auditing features for target={target}...") values_norm, mu, sigma, observed, feat_cols = prepare_day_matrix(day_df, feat_cols_override=base_cols) dl_tr, dl_va, _ = build_dataloaders(values_norm, observed, feat_cols, target) # Baseline model model = Seq2PointLSTM(input_size=len(feat_cols), hidden=128, layers=1, pred_len=PRED_LEN).to(DEVICE) base_va = train_one_day(model, dl_tr, dl_va, epochs=FS_EPOCHS) print(f"[{AUDIT_DAY}] baseline val MSE (audit) ~ {base_va:.6f}") deltas: Dict[str, float] = {} for i, col in enumerate(feat_cols): if col == target: deltas[col] = 0.0 continue va_mask = eval_val_loss_masked(model, dl_va, i) delta = (va_mask - base_va) / max(1e-12, base_va) deltas[col] = float(delta) harmful = [c for c, d in deltas.items() if d <= -FS_IMPROVE_THRESH and c != target] ignored = [c for c, d in deltas.items() if abs(d) < FS_NOISE_BAND and c != target] if harmful: print(f"[{AUDIT_DAY}] dropping harmful features (improve ≥{FS_IMPROVE_THRESH*100:.1f}%): {harmful}") if ignored: print(f"[{AUDIT_DAY}] (near-zero impact features): {ignored}") kept = [c for c in feat_cols if (c == target or c not in harmful)] print(f"[{AUDIT_DAY}] kept {len(kept)}/{len(feat_cols)} features") top_help = sorted([(c, d) for c, d in deltas.items() if d < 0], key=lambda x: x[1])[:8] top_hurt = sorted([(c, d) for c, d in deltas.items() if d > 0], key=lambda x: x[1], reverse=True)[:8] def fmt(entries): return ", ".join([f"{c}:{d*100:+.1f}%" for c,d in entries]) if entries else "-" print(f"[{AUDIT_DAY}] best removals (↓val): {fmt(top_help)}") print(f"[{AUDIT_DAY}] worst removals (↑val): {fmt(top_hurt)}") return kept, deltas # ---------------------------- # Orchestration # ---------------------------- def run_pipeline_for_target(target: str, feat_cols: List[str], run_tag: str, day_tables: Dict[str, pd.DataFrame]) -> Dict[str, float]: """ Train across days (inheriting weights). Save plots & underpred CSVs. Return dict of per-day best validation losses. """ print(f"\n=== [{run_tag}] Target: {target} with {len(feat_cols)} features ===") run_dir = ensure_dir(RESULT_DIR / run_tag) val_losses: Dict[str, float] = {} model = None for k, day in enumerate(DATE_STRS): print(f"[{day}] preparing data...") df_day = day_tables[day] values_norm, mu, sigma, observed, used_cols = prepare_day_matrix(df_day, feat_cols_override=feat_cols) input_size = values_norm.shape[1] if model is None: model = Seq2PointLSTM(input_size=input_size, hidden=128, layers=1, pred_len=PRED_LEN).to(DEVICE) dl_tr, dl_va, _ = build_dataloaders(values_norm, observed, used_cols, target) epochs = EPOCHS_FIRST if k == 0 else EPOCHS_NEXT print(f"[{day}] training ({epochs} epochs)...") va = train_one_day(model, dl_tr, dl_va, epochs) val_losses[day] = va print(f"[{day}] best val MSE ~ {va:.6f}") print(f"[{day}] forecasting...") pred = rolling_forecast_day(model, values_norm, observed, used_cols, target, mu, sigma) actual = df_day[target].astype("float32").to_numpy() # Under-prediction report csv_path = run_dir / f"{target}_underpred_{day}.csv" metrics = report_underpredictions( timestamps=df_day["timestamp"], actual=actual, pred=pred, out_csv_path=csv_path ) print( f"[{day}] {target} under-predictions: " f"{metrics['n_under']}/{metrics['n_pred_points']} " f"({metrics['under_rate']*100:.1f}%) | " f"MAE_under={metrics['mae_under']:.3f}, " f"sum_under={metrics['sum_under']:.3f}, " f"max_under={metrics['max_under']:.3f}, " f"p95_under={metrics['p95_under']:.3f}, " f"runs={metrics['n_runs']}, avg_run={metrics['avg_run']:.2f} min, " f"max_run={metrics['max_run']} min" ) print("saved under-prediction CSV →", csv_path) # Plot with under-prediction highlighting plt.figure(figsize=(12,4)) plt.plot(actual, label=f"actual {target}", lw=1) if np.isfinite(pred).any(): plt.plot(pred, label=f"forecast {target}", lw=2) finite = np.isfinite(actual) & np.isfinite(pred) under = finite & (pred < actual) if under.any(): x = np.arange(len(actual)) pred_f = np.where(np.isfinite(pred), pred, actual) plt.fill_between(x, pred_f, actual, where=under, alpha=0.15, color="red", label="under-predicted") plt.title(f"{target} – {PRED_LEN}-min horizon (LSTM, day {day}) [{run_tag}]") plt.xlabel("minutes (index within day)") plt.ylabel(target); plt.grid(True); plt.legend() outp = run_dir / f"{target}_forecast_{day}_lstm.png" plt.tight_layout(); plt.savefig(outp); plt.close() print("saved", outp) return val_losses def summarize_comparison(target: str, losses_base: Dict[str, float], losses_pruned: Dict[str, float]): days = sorted(losses_base.keys()) print(f"\n=== Comparison (target={target}) ===") for d in days: b, p = losses_base[d], losses_pruned[d] delta = (p - b) / max(1e-12, b) print(f"{d}: base={b:.6f} | pruned={p:.6f} | Δ={delta*100:+.2f}%") b_mean = np.mean([losses_base[d] for d in days]) p_mean = np.mean([losses_pruned[d] for d in days]) d_mean = (p_mean - b_mean) / max(1e-12, b_mean) print(f"Avg: base={b_mean:.6f} | pruned={p_mean:.6f} | Δ={d_mean*100:+.2f}%") def main(): ensure_dir(RESULT_DIR) # Build per-day feature tables once (independently per day) day_tables: Dict[str, pd.DataFrame] = {d: build_day_features(d) for d in DATE_STRS} base_cols = FOCUS_COLS + (["tod_sin","tod_cos"] if ADD_TIME_OF_DAY else []) for target in TARGETS: # 1) Baseline full-feature run losses_base = run_pipeline_for_target(target, feat_cols=base_cols, run_tag="baseline", day_tables=day_tables) # 2) Audit and select (on AUDIT_DAY) then re-run with pruned set if ENABLE_FEATURE_SELECTION: kept_cols = base_cols kept_map = {} try: kept_cols, deltas = feature_audit_and_select(day_tables[AUDIT_DAY], target, base_cols) kept_map = deltas except Exception as e: print(f"[WARN] feature audit failed ({e}); proceeding without pruning.") losses_pruned = run_pipeline_for_target(target, feat_cols=kept_cols, run_tag="pruned", day_tables=day_tables) # 3) Summarize before/after summarize_comparison(target, losses_base, losses_pruned) # Save audit table audit_path = ensure_dir(RESULT_DIR / "pruned") / f"{target}_feature_audit_{AUDIT_DAY}.csv" if kept_map: pd.DataFrame( [{"feature": c, "delta_val_pct": 100.0 * kept_map[c], "dropped": (c not in kept_cols)} for c in kept_map] ).sort_values("delta_val_pct").to_csv(audit_path, index=False) print("saved feature audit →", audit_path) else: print("[INFO] Feature selection disabled; only baseline run executed.") if __name__ == "__main__": main()