next-tjwater-server/app/services/scada_monitoring_calibration.py

"""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}