773 lines
30 KiB
Python
773 lines
30 KiB
Python
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<ts>.*?)\]
|
||
\s+UE\s+RNTI=(?P<rnti>\d+)
|
||
\s+DL_TBS=(?P<dl_tbs>\d+)
|
||
\s+UL_TBS=(?P<ul_tbs>\d+)
|
||
\s+DL_SDUs=(?P<dl_sdus>\d+)
|
||
\s+UL_SDUs=(?P<ul_sdus>\d+)
|
||
\s+DL_sched_RBs=(?P<dl_sched_rbs>\d+)
|
||
\s+UL_sched_RBs=(?P<ul_sched_rbs>\d+)
|
||
\s+DL_PRBs=(?P<dl_prbs>\d+)
|
||
\s+UL_PRBs=(?P<ul_prbs>\d+)
|
||
\s+DL_BLER=(?P<dl_bler>[\d.]+)
|
||
\s+UL_BLER=(?P<ul_bler>[\d.]+)
|
||
\s+PUSCH_SNR=(?P<pusch_snr>[\d.]+)
|
||
\s+PUCCH_SNR=(?P<pucch_snr>[\d.]+)
|
||
\s+WB_CQI=(?P<wb_cqi>\d+)
|
||
\s+DL_MCS1=(?P<dl_mcs1>\d+)
|
||
\s+UL_MCS1=(?P<ul_mcs1>\d+)
|
||
\s+DL_HARQ=(?P<dl_harq>\d+)
|
||
\s+UL_HARQ=(?P<ul_harq>\d+)
|
||
\s+BSR=(?P<bsr>\d+)
|
||
\s+PHR=(?P<phr>\d+)
|
||
""", re.VERBOSE,
|
||
)
|
||
RRC_PATTERNS: List[re.Pattern] = [
|
||
re.compile(r"\[(?P<ts>[^]]+)\].*RNTI=(?P<rnti>[0-9a-fA-F]+).*RRC[_ ]?MSG=(?P<msg>[A-Za-z_]+)", re.IGNORECASE),
|
||
re.compile(r"\[(?P<ts>[^]]+)\].*rrc[=:](?P<msg>[A-Za-z_]+).*rnti[=:](?P<rnti>[0-9a-fA-F]+)", re.IGNORECASE),
|
||
re.compile(r"\[(?P<ts>[^]]+)\].*(?P<msg>RRC[A-Za-z]+).*RNTI=(?P<rnti>[0-9a-fA-F]+)", re.IGNORECASE),
|
||
re.compile(r"\|\s*(?P<ts>\d{4}-\d{2}-\d{2}T\d{2}:\d{2}:\d{2}(?:\.\d+)?Z?).*?RNTI\s*(?P<rnti>[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()
|