from __future__ import annotations from dataclasses import dataclass from typing import Any import numpy as np CHANNELS = ( "demand_norm", "stack_load_norm", "feed_flow_norm", "coolant_flow_norm", "stack_temperature_norm", "discharge_pressure_norm", "auxiliary_current_norm", "valve_position_norm", ) @dataclass(frozen=True) class TelemetryRun: run_id: str values: np.ndarray labels: np.ndarray anomaly_type: str event_start: int | None cadence_seconds: float = 1.0 def generate_telemetry_run( seed: int, anomaly_type: str = "none", length: int = 360, ) -> TelemetryRun: """Create a normalized synthetic startup/run/shutdown telemetry episode.""" supported = { "none", "sensor_drift", "valve_stiction", "cooling_degradation", "auxiliary_trip", } if anomaly_type not in supported: raise ValueError(f"unsupported anomaly type: {anomaly_type}") if length < 240: raise ValueError("length must be at least 240 samples") rng = np.random.default_rng(seed) demand = np.zeros(length, dtype=np.float64) ramp_up_end = min(105, length // 3) run_end = max(ramp_up_end + 40, length - 80) ramp_down_end = min(length, run_end + 60) demand[30:ramp_up_end] = np.linspace(0.0, 0.72, ramp_up_end - 30, endpoint=False) demand[ramp_up_end:run_end] = 0.72 + 0.035 * np.sin( np.linspace(0.0, 8.0 * np.pi, run_end - ramp_up_end) ) demand[run_end:ramp_down_end] = np.linspace( demand[run_end - 1], 0.0, ramp_down_end - run_end, endpoint=False ) values = np.zeros((length, len(CHANNELS)), dtype=np.float64) load = flow = coolant = temp = pressure = aux = valve = 0.0 for tick in range(length): load += 0.18 * (demand[tick] - load) + rng.normal(0.0, 0.003) valve += 0.25 * ((0.04 + 0.88 * demand[tick]) - valve) + rng.normal(0.0, 0.004) flow += 0.20 * ((0.03 + 0.86 * load) - flow) + rng.normal(0.0, 0.005) coolant += 0.16 * ((0.08 + 0.76 * load) - coolant) + rng.normal(0.0, 0.004) temp += 0.055 * ((0.22 + 0.52 * load) - temp) + rng.normal(0.0, 0.002) pressure += 0.12 * ((0.12 + 0.58 * flow) - pressure) + rng.normal(0.0, 0.003) aux += 0.20 * ((0.10 + 0.56 * load) - aux) + rng.normal(0.0, 0.004) values[tick] = demand[tick], load, flow, coolant, temp, pressure, aux, valve labels = np.zeros(length, dtype=np.int8) event_start: int | None = None if anomaly_type != "none": event_start = min(max(150, length // 2), length - 70) labels[event_start:] = 1 horizon = length - event_start # Reach the full injected deviation within 90 seconds and then hold. # This keeps event-delay metrics interpretable at the 1 Hz cadence. severity = np.minimum(1.0, np.arange(horizon, dtype=np.float64) / 90.0) if anomaly_type == "sensor_drift": values[event_start:, 5] += 0.34 * severity elif anomaly_type == "valve_stiction": stuck = values[event_start, 7] values[event_start:, 7] = stuck values[event_start:, 2] -= 0.30 * severity values[event_start:, 5] -= 0.12 * severity elif anomaly_type == "cooling_degradation": values[event_start:, 3] -= 0.42 * severity values[event_start:, 4] += 0.38 * severity elif anomaly_type == "auxiliary_trip": values[event_start:, 6] *= np.exp(-np.linspace(0.0, 8.0, horizon)) values[event_start:, 3] -= 0.24 * severity values[event_start:, 4] += 0.20 * severity return TelemetryRun( run_id=f"synthetic-{seed:04d}", values=np.clip(values, 0.0, 1.2), labels=labels, anomaly_type=anomaly_type, event_start=event_start, ) def causal_window_features( values: np.ndarray, window_size: int = 24 ) -> tuple[np.ndarray, np.ndarray]: if values.ndim != 2 or values.shape[1] != len(CHANNELS): raise ValueError(f"values must have shape [time, {len(CHANNELS)}]") if len(values) < window_size: raise ValueError("telemetry run is shorter than the causal window") x_axis = np.linspace(-1.0, 1.0, window_size) x_energy = float(np.sum(x_axis**2)) rows = [] ticks = [] for end in range(window_size - 1, len(values)): window = values[end - window_size + 1 : end + 1] centered = window - np.mean(window, axis=0) slope = (x_axis[:, None] * centered).sum(axis=0) / x_energy rows.append( np.concatenate( [ window[-1], np.mean(window, axis=0), np.std(window, axis=0), slope, window[-1] - window[-2], ] ) ) ticks.append(end) return np.asarray(rows, dtype=np.float64), np.asarray(ticks, dtype=np.int64) class RobustPCADetector: """Compact unsupervised baseline with robust scaling and PCA residuals.""" def __init__( self, center: np.ndarray, scale: np.ndarray, components: np.ndarray, score_center: float, score_scale: float, threshold: float, window_size: int, ) -> None: self.center = np.asarray(center, dtype=np.float64) self.scale = np.asarray(scale, dtype=np.float64) self.components = np.asarray(components, dtype=np.float64) self.score_center = float(score_center) self.score_scale = float(score_scale) self.threshold = float(threshold) self.window_size = int(window_size) @staticmethod def _raw_score(z: np.ndarray, components: np.ndarray) -> np.ndarray: projection = z @ components.T reconstruction = projection @ components return np.mean((z - reconstruction) ** 2, axis=1) @classmethod def fit( cls, train_features: np.ndarray, calibration_features: np.ndarray, window_size: int = 24, explained_variance: float = 0.92, calibration_quantile: float = 0.995, ) -> RobustPCADetector: center = np.median(train_features, axis=0) mad = np.median(np.abs(train_features - center), axis=0) scale = np.maximum(1.4826 * mad, 1e-4) z_train = np.clip((train_features - center) / scale, -12.0, 12.0) _, singular, vh = np.linalg.svd(z_train, full_matrices=False) variance = singular**2 cumulative = np.cumsum(variance) / np.sum(variance) component_count = int(np.searchsorted(cumulative, explained_variance) + 1) component_count = min(max(component_count, 2), train_features.shape[1] - 1) components = vh[:component_count] raw_train = cls._raw_score(z_train, components) score_center = float(np.median(raw_train)) score_mad = float(np.median(np.abs(raw_train - score_center))) score_scale = max(1.4826 * score_mad, 1e-6) z_cal = np.clip((calibration_features - center) / scale, -12.0, 12.0) calibrated = (cls._raw_score(z_cal, components) - score_center) / score_scale threshold = float(np.quantile(calibrated, calibration_quantile)) return cls( center, scale, components, score_center, score_scale, threshold, window_size, ) def score_features(self, features: np.ndarray) -> np.ndarray: z = np.clip((features - self.center) / self.scale, -12.0, 12.0) raw = self._raw_score(z, self.components) return (raw - self.score_center) / self.score_scale def score_run(self, run: TelemetryRun) -> tuple[np.ndarray, np.ndarray]: features, ticks = causal_window_features(run.values, self.window_size) return ticks, self.score_features(features) @staticmethod def persistent_alerts(scores: np.ndarray, threshold: float, samples: int = 3) -> np.ndarray: """Require consecutive exceedances to suppress isolated nuisance alerts.""" if samples < 1: raise ValueError("samples must be positive") raw = np.asarray(scores) > threshold if samples == 1: return raw counts = np.convolve(raw.astype(np.int8), np.ones(samples, dtype=np.int8), mode="full") return counts[: len(raw)] >= samples def to_dict(self) -> dict[str, Any]: return { "schema_version": "forge.anomaly.robust-pca.v1", "channels": list(CHANNELS), "window_size": self.window_size, "center": self.center.tolist(), "scale": self.scale.tolist(), "components": self.components.tolist(), "score_center": self.score_center, "score_scale": self.score_scale, "threshold": self.threshold, "authority": "L0 evidence only; never an interlock or control input", } @classmethod def from_dict(cls, payload: dict[str, Any]) -> RobustPCADetector: if payload["schema_version"] != "forge.anomaly.robust-pca.v1": raise ValueError("unsupported anomaly detector artifact") if tuple(payload["channels"]) != CHANNELS: raise ValueError("channel contract mismatch") return cls( center=np.asarray(payload["center"]), scale=np.asarray(payload["scale"]), components=np.asarray(payload["components"]), score_center=payload["score_center"], score_scale=payload["score_scale"], threshold=payload["threshold"], window_size=payload["window_size"], )