"""Fit causal replay thresholds in bounded numeric arrays, before publishing history.""" from datetime import datetime, timedelta import numpy as np import pandas as pd from app.domain.schemas.project_settings import ScadaMonitoringSettings from app.domain.schemas.scada_replay import ScadaReplayCreate from app.services.scada_normal_profile import CALIBRATION_VERSION, DIAGNOSTIC_POLICY def adaptive_score_threshold(scores: list[float]) -> dict: """Keep the cutoff below the normal lower tail, with a fixed noise margin.""" values = np.asarray(scores, dtype=float) if len(values) < 5 or not np.isfinite(values).all(): return {"status": "insufficient_data", "reason": "normal_score_history_incomplete"} center = float(np.median(values)) sigma = float(1.4826 * np.median(np.abs(values - center))) lower_tail = float(np.quantile(values, 0.05)) return { "status": "ready", "threshold": min(-0.04, center - 3 * sigma, lower_tail - 0.02), "normal_score_median": center, "normal_score_min": float(values.min()), "normal_score_max": float(values.max()), "normal_score_robust_sigma": sigma, "normal_score_lower_quantile": lower_tail, "normal_score_count": len(values), "method": "min(-0.04, median-3*robust_sigma, normal_q05-0.02)", } class MonitoringCalibrationBuilder: """Retain bounded numeric series, never full raw pages or device batches.""" def __init__(self, document: dict, config: ScadaReplayCreate): self.settings = ScadaMonitoringSettings.model_validate(document["monitoring_config"]) self.config = config self.cutoff = config.start_time self.start = self.cutoff - timedelta(days=15) self.interval = int(document["monitoring_interval_seconds"]) self.supported = self.interval % 60 == 0 and 86400 % self.interval == 0 self.count = 15 * (86400 // self.interval) if self.supported else 0 ids = dict.fromkeys(self.settings.pressure_device_ids + self.settings.outlet_flow_device_ids + [device for rule in config.rules if rule.kind == "relative_change" for device in rule.device_ids]) self.values = {device: np.full(self.count, np.nan) for device in ids} def add(self, rows: list[dict]) -> None: for row in rows: device = row["device_id"] if device not in self.values: continue at = datetime.fromisoformat(row["time"]) if at >= self.cutoff: raise ValueError("动态阈值拟合不能读取回放开始时刻及之后的观测") offset = (at - self.start).total_seconds() if (not self.supported or offset < 0 or offset % self.interval or row.get("quality") not in {"valid", "corrected"} or row.get("cleaned_value") is None): continue index = int(offset // self.interval) if np.isfinite(self.values[device][index]): raise ValueError("动态阈值拟合存在重复的设备时间点") self.values[device][index] = float(row["cleaned_value"]) def fit(self) -> dict: base = { "version": CALIBRATION_VERSION, "diagnostic_policy": DIAGNOSTIC_POLICY, "knowledge_cutoff": self.cutoff.isoformat(), "history_start": self.start.isoformat(), "required_samples_per_device": self.count, "pressure": {}, "outlet_flow": {}, "diagnostic": {}, } if not self.supported: return {**base, "status": "insufficient_data", "reason": "unsupported_sampling_interval"} available = {device: int(np.isfinite(values).sum()) for device, values in self.values.items()} # A regular diurnal rise/reversal can be a huge percentage of a small # preceding flow. Judge its absolute value against causal normal readings # at the same sampling phase before treating that percentage as abnormal. points_per_day = 86400 // self.interval for device, values in self.values.items(): if available[device] != self.count: continue days = values.reshape(15, points_per_day) center = np.median(days, axis=0) sigma = 1.4826 * np.median(np.abs(days - center), axis=0) margin = np.maximum(3 * sigma, np.maximum(np.abs(center) * 1e-6, 1e-9)) base["diagnostic"][device] = { "phase_origin": self.start.isoformat(), "interval_seconds": self.interval, "lower": (days.min(axis=0) - margin).tolist(), "upper": (days.max(axis=0) + margin).tolist(), "method": "normal_same_phase_range_plus_3_robust_sigma", "normal_day_count": 15, } for kind, ids, direction, floor in ( ("pressure", self.settings.pressure_device_ids, -1, 3.0), ("outlet_flow", self.settings.outlet_flow_device_ids, 1, 5.0), ): for device in ids: if available[device] != self.count: continue rule = next((rule for rule in self.config.rules if rule.kind == "relative_change" and device in rule.device_ids), None) if rule is None: continue series = pd.Series(self.values[device]) reference = series.rolling(rule.window_samples, min_periods=rule.minimum_samples).median().shift(1) changes = (direction * 100 * (series - reference) / reference.abs()).to_numpy() changes = changes[np.isfinite(changes)] if not len(changes): continue center = float(np.median(changes)) sigma = float(1.4826 * np.median(np.abs(changes - center))) base[kind][device] = { "threshold_percent": max(floor, center + 3 * sigma), "normal_change_median": center, "normal_change_robust_sigma": sigma, "normal_change_count": int(len(changes)), "method": "normal_recent_window_change_median_plus_3_robust_sigma", "normal_profile": base["diagnostic"].get(device), } pressure_ids = self.settings.pressure_device_ids incomplete = {device: available[device] for device in pressure_ids if available[device] != self.count} if len(pressure_ids) < 5 or incomplete: return {**base, "status": "insufficient_data", "reason": "pressure_history_incomplete", "available_samples": available} from app.services import scada_replay_monitoring as monitoring index = pd.date_range(start=self.start, periods=self.count, freq=pd.Timedelta(seconds=self.interval)) frame = pd.DataFrame({device: self.values[device] for device in pressure_ids}, index=index) result = monitoring.run_burst_detection( network="scada-replay-baseline", username="scada-replay", observed_pressure_data=frame, points_per_day=86400 // self.interval, mu=monitoring.TARGET_MU, iforest_params={"n_estimators": monitoring.TARGET_N_ESTIMATORS, "random_state": monitoring.TARGET_RANDOM_STATE, "contamination": "auto"}, sensor_nodes=pressure_ids, ) score = adaptive_score_threshold([float(row["Score"]) for row in result["rows"]]) return {**base, "status": score["status"], "burst_detection": score}