from __future__ import annotations from collections import Counter from datetime import datetime, timedelta from typing import Any import numpy as np import pandas as pd from app.algorithms.burst_detection.burst_detector import BurstDetector from app.infra.db.timescaledb.internal_queries import InternalQueries from app.services.scheme_management import ( query_burst_detection_scheme_detail, query_burst_detection_schemes, scheme_name_exists, store_scheme_info, ) from app.services.tjnetwork import get_all_scada_info from app.services.time_api import extract_date, parse_utc_time, utc_now TARGET_DAY_COUNT = 15 DEFAULT_SAMPLE_INTERVAL_MINUTES = 15 TARGET_MU = 1 TARGET_N_ESTIMATORS = 50 TARGET_RANDOM_STATE = 42 TARGET_SCORE_THRESHOLD = -0.04 MIN_COMPLETE_SENSORS = 5 def run_burst_detection( *, network: str, username: str, observed_pressure_data: ( pd.DataFrame | dict[str, list[Any]] | list[dict[str, Any]] | list[list[Any]] | None ) = None, points_per_day: int = 1440, mu: int = 100, iforest_params: dict[str, Any] | None = None, target_time: datetime | str | None = None, sampling_interval_minutes: int | None = None, scada_start: datetime | str | None = None, scada_end: datetime | str | None = None, sensor_nodes: list[str] | None = None, scheme_name: str | None = None, data_source: str = "monitoring", simulation_scheme_name: str | None = None, simulation_scheme_type: str | None = None, ) -> dict[str, Any]: """ 运行爆管侦测服务入口。 调用方式三选一： - 不传数据时间窗，自动侦测最近完整时刻；可用 `target_time` 回放历史时刻 - 直接传 `observed_pressure_data` - 或传 `scada_start/scada_end` 让后端自动查询 SCADA 压力数据 `observed_pressure_data` 支持格式： - `pd.DataFrame` 行表示时间点，列表示传感器；列名应为传感器/节点 ID。 - `dict[str, list[Any]]` 键为传感器/节点 ID，值为按时间顺序排列的压力序列。 例如：`{"J1": [101.2, 101.0], "J2": [99.8, 99.7]}`。 - `list[dict[str, Any]]` 每个元素代表一个时间点的多传感器观测。 例如：`[{"J1": 101.2, "J2": 99.8}, {"J1": 101.0, "J2": 99.7}]`。 - `list[list[Any]]` 二维数组式 JSON，格式为 `(时间点数, 传感器数)`。 这是最接近原始 `burst_detector` 示例代码的调用方式。 数据约束： - 统一要求“行=时间点，列=传感器”。 - 总样本点数必须能被 `points_per_day` 整除。 - 至少要有 2 天数据，即 `sample_count >= 2 * points_per_day`。 - 若传入 `sensor_nodes`，输入数据必须包含这些列；SCADA 模式下也会只按这些节点取数。 """ if not network: raise ValueError("network is required.") selected_sensor_nodes = ( list(dict.fromkeys([node for node in (sensor_nodes or []) if node])) if sensor_nodes else None ) use_scada_source = scada_start is not None or scada_end is not None use_target_mode = ( observed_pressure_data is None and not use_scada_source and data_source != "simulation" ) or target_time is not None resolved_target_time: datetime | None = None requested_target_time: datetime | None = None excluded_sensors: list[dict[str, str]] = [] daily_times: list[datetime] | None = None resolved_sampling_interval_minutes: int | None = None if use_target_mode: if observed_pressure_data is not None or use_scada_source: raise ValueError( "target_time 不能与 observed_pressure_data 或 scada_start/scada_end 同时使用。" ) scada_sensor_nodes = ( selected_sensor_nodes if selected_sensor_nodes is not None else _get_pressure_sensor_nodes(network) ) requested_target_time = ( _to_datetime(target_time) if target_time is not None else None ) resolved_sampling_interval_minutes = _resolve_sampling_interval_minutes( network=network, sensor_nodes=scada_sensor_nodes, requested_interval=sampling_interval_minutes, ) target_points_per_day = 1440 // resolved_sampling_interval_minutes ( observed_input, resolved_target_time, excluded_sensors, ) = _build_target_pressure_from_scada( network=network, sensor_nodes=scada_sensor_nodes, requested_target_time=requested_target_time, sampling_interval_minutes=resolved_sampling_interval_minutes, points_per_day=target_points_per_day, ) selected_sensor_nodes = list(observed_input.columns) observed_source = ( "latest_monitoring" if target_time is None else "historical_monitoring" ) points_per_day = target_points_per_day mu = TARGET_MU iforest_params = { "n_estimators": TARGET_N_ESTIMATORS, "random_state": TARGET_RANDOM_STATE, "contamination": "auto", } daily_times = [ resolved_target_time - timedelta(days=offset) for offset in range(TARGET_DAY_COUNT - 1, -1, -1) ] elif use_scada_source: scada_sensor_nodes = ( selected_sensor_nodes if selected_sensor_nodes is not None else _get_pressure_sensor_nodes(network) ) if data_source == "simulation": if not simulation_scheme_name: raise ValueError("模拟方案模式必须提供 simulation_scheme_name。") observed_df = _build_observed_pressure_from_simulation( network=network, sensor_nodes=scada_sensor_nodes, scada_start=scada_start, scada_end=scada_end, simulation_scheme_name=simulation_scheme_name, simulation_scheme_type=simulation_scheme_type, ) observed_input = observed_df observed_source = "simulation_scheme_timerange" else: observed_df = _build_observed_pressure_from_scada( network=network, sensor_nodes=scada_sensor_nodes, scada_start=scada_start, scada_end=scada_end, ) observed_input = observed_df observed_source = "backend_timerange" else: if observed_pressure_data is None: raise ValueError( "未提供 observed_pressure_data，且未提供 scada_start/scada_end。" ) observed_input = observed_pressure_data observed_source = "request_payload" detector = BurstDetector( mu=mu, points_per_day=points_per_day, iforest_params=iforest_params, ) result_df = detector.run_detection( observed_input, sensor_nodes=selected_sensor_nodes, ) resolved_sensor_nodes = list(result_df.attrs.get("sensor_nodes", [])) rows = _serialize_result_rows( result_df, daily_times=daily_times, target_only=use_target_mode, ) summary = _build_detection_summary( result_df, daily_times=daily_times, target_only=use_target_mode, ) payload: dict[str, Any] = { "network": network, "sensor_nodes": resolved_sensor_nodes, "observed_source": observed_source, "sample_count": int(result_df.attrs.get("sample_count", 0)), "points_per_day": int(result_df.attrs.get("points_per_day", points_per_day)), "day_count": int(result_df.attrs.get("day_count", len(result_df))), "rows": rows, "summary": summary, "algorithm_params": { "mu": mu, "points_per_day": points_per_day, "iforest_params": detector.iforest_params, **( {"score_threshold": TARGET_SCORE_THRESHOLD} if use_target_mode else {} ), }, } if data_source == "simulation": payload["data_source"] = "simulation" payload["simulation_scheme"] = { "name": simulation_scheme_name, "type": simulation_scheme_type, } else: payload["data_source"] = "monitoring" if ( use_target_mode and resolved_target_time is not None and resolved_sampling_interval_minutes is not None ): sample_start = resolved_target_time - timedelta( days=TARGET_DAY_COUNT, minutes=-resolved_sampling_interval_minutes, ) payload.update( { "requested_target_time": ( requested_target_time.isoformat() if requested_target_time is not None else None ), "target_time": resolved_target_time.isoformat(), "reference_window": { "start": (resolved_target_time - timedelta(days=14)).isoformat(), "end": (resolved_target_time - timedelta(days=1)).isoformat(), "day_count": 14, }, "sampling_interval_minutes": resolved_sampling_interval_minutes, "daily_scores": [ { "timestamp": row["Timestamp"], "role": row["Role"], "score": row["Score"], "raw_prediction": row["Prediction"], } for row in rows ], "data_quality": { "included_sensors": resolved_sensor_nodes, "excluded_sensors": excluded_sensors, "minimum_required_sensors": MIN_COMPLETE_SENSORS, }, "scada_window": { "start": sample_start.isoformat(), "end": resolved_target_time.isoformat(), }, } ) elif use_scada_source: payload["scada_window"] = { "start": _to_datetime(scada_start).isoformat(), "end": _to_datetime(scada_end).isoformat(), } if scheme_name: _store_burst_detection_scheme( network=network, scheme_name=scheme_name, username=username, payload=payload, mu=mu, points_per_day=points_per_day, iforest_params=detector.iforest_params, ) payload["scheme_name"] = scheme_name return payload def _build_observed_pressure_from_simulation( *, network: str, sensor_nodes: list[str], scada_start: datetime | str | None, scada_end: datetime | str | None, simulation_scheme_name: str | None, simulation_scheme_type: str | None = None, ) -> pd.DataFrame: if scada_start is None or scada_end is None: raise ValueError("使用模拟方案查询时必须同时提供 scada_start 与 scada_end。") start_dt = _to_datetime(scada_start) end_dt = _to_datetime(scada_end) if start_dt >= end_dt: raise ValueError("SCADA 时间窗非法：scada_start 必须早于 scada_end。") # Reuse burst_location logic partially here or call internal queries directly # For burst detection, we need time series for all sensor nodes. # Check for missing nodes in simulation result if needed, but InternalQueries handles some of it. # We assume sensor_nodes are valid pressure nodes. scheme_type = simulation_scheme_type or "burst_analysis" simulation_data = InternalQueries.query_scheme_simulation_by_ids_timerange( db_name=network, scheme_type=scheme_type, scheme_name=simulation_scheme_name, element_ids=sensor_nodes, start_time=start_dt.isoformat(), end_time=end_dt.isoformat(), element_type="node", field="pressure", ) # simulation_data is {sensor_id: [{time, value}, ...]} # Convert to DataFrame: index=time, columns=sensor_ids data_dict = {} timestamps = set() for sensor_id in sensor_nodes: records = simulation_data.get(sensor_id, []) if not records: continue # Convert records to Series with time index ts_values = [] ts_index = [] for r in records: if r.get("value") is not None: ts_values.append(float(r["value"])) ts_index.append(pd.to_datetime(r["time"])) if ts_values: s = pd.Series(ts_values, index=ts_index) data_dict[sensor_id] = s timestamps.update(ts_index) if not data_dict: raise ValueError("指定时间窗内未查询到模拟压力数据。") observation_df = pd.DataFrame(data_dict) # Handle missing timestamps if any (though simulation usually has uniform steps) # Forward fill or interpolate might be needed if steps differ, but typically they align. observation_df = observation_df.sort_index() # Fill NaN if any missing points for some sensors observation_df = observation_df.fillna(method="ffill").fillna(method="bfill") if observation_df.empty: raise ValueError("模拟压力数据无法构建观测矩阵。") return observation_df def list_burst_detection_schemes( network: str, query_date: datetime | str | None = None, ) -> list[dict[str, Any]]: parsed_date = extract_date(query_date, field_name="query_date") if query_date is not None else None return query_burst_detection_schemes( name=network, network=network, query_date=parsed_date, ) def get_burst_detection_scheme_detail(network: str, scheme_name: str) -> dict[str, Any]: result = query_burst_detection_scheme_detail(network, scheme_name) if not result: raise ValueError(f"未找到爆管侦测方案: {scheme_name}") return result def _store_burst_detection_scheme( *, network: str, scheme_name: str, username: str, payload: dict[str, Any], mu: int, points_per_day: int, iforest_params: dict[str, Any], ) -> None: if scheme_name_exists(network, scheme_name): raise ValueError(f"方案名称已存在: {scheme_name}") now_iso = utc_now().isoformat() scheme_detail = { "network": network, "sensor_nodes": payload.get("sensor_nodes", []), "observed_source": payload.get("observed_source"), "scada_window": payload.get("scada_window"), "algorithm_params": { "mu": mu, "points_per_day": points_per_day, "iforest_params": iforest_params, }, "result_summary": payload.get("summary", {}), "result_payload": payload, } store_scheme_info( name=network, scheme_name=scheme_name, scheme_type="burst_detection", username=username, scheme_start_time=now_iso, scheme_detail=scheme_detail, ) def _serialize_result_rows( result_df: pd.DataFrame, *, daily_times: list[datetime] | None = None, target_only: bool = False, ) -> list[dict[str, Any]]: rows: list[dict[str, Any]] = [] raw_rows = result_df.to_dict(orient="records") for index, row in enumerate(raw_rows): is_target = index == len(raw_rows) - 1 is_burst = bool(row["IsBurst"]) if target_only: is_burst = is_target and float(row["Score"]) <= TARGET_SCORE_THRESHOLD rows.append( { "Day": int(row["Day"]), "Score": float(row["Score"]), "Prediction": int(row["Prediction"]), "IsBurst": is_burst, **( { "Timestamp": daily_times[index].isoformat(), "Role": "target" if is_target else "reference", } if daily_times is not None else {} ), } ) return rows def _build_detection_summary( result_df: pd.DataFrame, *, daily_times: list[datetime] | None = None, target_only: bool = False, ) -> dict[str, Any]: rows = _serialize_result_rows( result_df, daily_times=daily_times, target_only=target_only, ) if not rows: raise ValueError("爆管侦测结果为空。") score_series = result_df["Score"] most_anomalous_index = int(score_series.idxmin()) latest_row = rows[-1] anomaly_days = [row["Day"] for row in rows if row["IsBurst"]] summary = { "burst_detected": bool(latest_row["IsBurst"]), "latest_day": latest_row, "most_anomalous_day": int(result_df.iloc[most_anomalous_index]["Day"]), "anomaly_days": anomaly_days, "anomaly_day_count": len(anomaly_days), "latest_sensor_rankings": _build_latest_sensor_rankings(result_df), } if target_only: target_score = float(latest_row["Score"]) summary.update( { "target_score": target_score, "score_threshold": TARGET_SCORE_THRESHOLD, "target_rank": int(result_df["Score"].rank(method="min").iloc[-1]), "target_time": latest_row.get("Timestamp"), "reference_day_count": TARGET_DAY_COUNT - 1, } ) return summary def _build_latest_sensor_rankings(result_df: pd.DataFrame) -> list[dict[str, Any]]: feature_matrix = result_df.attrs.get("high_freq_features") sensor_nodes = list(result_df.attrs.get("sensor_nodes", [])) if feature_matrix is None or len(sensor_nodes) == 0: return [] latest_values = np.asarray(feature_matrix[-1], dtype=float) history = np.asarray(feature_matrix[:-1], dtype=float) history_means = history.mean(axis=0) history_stds = history.std(axis=0) safe_stds = np.where(history_stds > 1e-9, history_stds, 1e-9) deviations = (latest_values - history_means) / safe_stds ranking = sorted( zip( sensor_nodes, latest_values, history_means, history_stds, deviations, strict=False, ), key=lambda item: item[4], ) return [ { "sensor_node": sensor_id, "latest_high_frequency_value": float(value), "historical_mean": float(history_mean), "historical_std": float(history_std), "standardized_deviation": float(deviation), } for sensor_id, value, history_mean, history_std, deviation in ranking[ : min(10, len(ranking)) ] ] def _build_target_pressure_from_scada( *, network: str, sensor_nodes: list[str], requested_target_time: datetime | None, sampling_interval_minutes: int, points_per_day: int, ) -> tuple[pd.DataFrame, datetime, list[dict[str, str]]]: node_query_id = _get_pressure_sensor_mapping(network) mapped_nodes = [node for node in sensor_nodes if node in node_query_id] excluded_without_mapping = [ {"sensor_node": node, "reason": "missing_api_query_id"} for node in sensor_nodes if node not in node_query_id ] if len(mapped_nodes) < MIN_COMPLETE_SENSORS: raise ValueError( f"可查询的压力测点少于 {MIN_COMPLETE_SENSORS} 个，无法执行爆管侦测。" ) query_ids = [node_query_id[node] for node in mapped_nodes] candidate_before = requested_target_time last_excluded: list[dict[str, str]] = excluded_without_mapping for _ in range(4): resolved_target = InternalQueries.query_latest_scada_time( db_name=network, device_ids=query_ids, before_time=candidate_before, ) if resolved_target is None: break sample_start = resolved_target - timedelta( days=TARGET_DAY_COUNT, minutes=-sampling_interval_minutes, ) expected_index = pd.date_range( start=sample_start, end=resolved_target, freq=f"{sampling_interval_minutes}min", ) scada_data = InternalQueries.query_scada_by_ids_timerange( db_name=network, device_ids=query_ids, start_time=sample_start, end_time=resolved_target, ) complete_columns: dict[str, pd.Series] = {} excluded = list(excluded_without_mapping) for node_id in mapped_nodes: query_id = node_query_id[node_id] records = scada_data.get(query_id, []) if not records: excluded.append({"sensor_node": node_id, "reason": "no_data"}) continue record_frame = pd.DataFrame.from_records(records) record_frame["time"] = pd.to_datetime(record_frame["time"], utc=True) record_frame["value"] = pd.to_numeric( record_frame["value"], errors="coerce" ) series = ( record_frame.drop_duplicates(subset="time", keep="last") .set_index("time")["value"] .reindex(expected_index) ) if len(series) != TARGET_DAY_COUNT * points_per_day: excluded.append( {"sensor_node": node_id, "reason": "unexpected_sample_count"} ) continue if series.isna().any(): excluded.append( {"sensor_node": node_id, "reason": "missing_or_invalid_samples"} ) continue complete_columns[node_id] = series if len(complete_columns) >= MIN_COMPLETE_SENSORS: observation_df = pd.DataFrame(complete_columns, index=expected_index) return observation_df, resolved_target, excluded last_excluded = excluded candidate_before = resolved_target - timedelta(microseconds=1) excluded_preview = ", ".join( item["sensor_node"] for item in last_excluded[:10] ) raise ValueError( f"最近数据中完整压力测点少于 {MIN_COMPLETE_SENSORS} 个；" f"请检查 15 天数据完整性。排除测点: {excluded_preview or '无'}" ) def _resolve_sampling_interval_minutes( *, network: str, sensor_nodes: list[str], requested_interval: int | None, ) -> int: if requested_interval is not None: interval = int(requested_interval) else: selected_nodes = set(sensor_nodes) inferred_intervals = [ parsed for item in get_all_scada_info(network) if str(item.get("type", "")).lower() == "pressure" and str(item.get("associated_element_id", "")) in selected_nodes and ( parsed := _parse_sampling_interval_minutes( item.get("transmission_frequency") ) ) is not None ] interval = ( Counter(inferred_intervals).most_common(1)[0][0] if inferred_intervals else DEFAULT_SAMPLE_INTERVAL_MINUTES ) if interval <= 0 or 1440 % interval != 0: raise ValueError("采样间隔必须是能整除 1440 分钟的正整数。") return interval def _parse_sampling_interval_minutes(value: Any) -> int | None: if value is None: return None if isinstance(value, (int, float)): minutes = float(value) else: try: minutes = pd.to_timedelta(str(value)).total_seconds() / 60 except (TypeError, ValueError): return None rounded = round(minutes) if minutes <= 0 or abs(minutes - rounded) > 1e-6: return None return int(rounded) def _get_pressure_sensor_mapping(network: str) -> dict[str, str]: node_query_id: dict[str, str] = {} for item in get_all_scada_info(network): if str(item.get("type", "")).lower() != "pressure": continue node_id = item.get("associated_element_id") query_id = item.get("api_query_id") if node_id and query_id is not None: node_query_id[str(node_id)] = str(query_id) return node_query_id def _get_pressure_sensor_nodes(network: str) -> list[str]: sensor_nodes: list[str] = [] for item in get_all_scada_info(network): if str(item.get("type", "")).lower() != "pressure": continue node_id = item.get("associated_element_id") if isinstance(node_id, str) and node_id: sensor_nodes.append(node_id) sensor_nodes = list(dict.fromkeys(sensor_nodes)) if not sensor_nodes: raise ValueError("未找到压力传感器对应节点（scada_info.type=pressure）。") return sensor_nodes def _build_observed_pressure_from_scada( *, network: str, sensor_nodes: list[str], scada_start: datetime | str | None, scada_end: datetime | str | None, ) -> pd.DataFrame: if scada_start is None or scada_end is None: raise ValueError("使用后端 SCADA 查询时必须同时提供 scada_start 与 scada_end。") start_dt = _to_datetime(scada_start) end_dt = _to_datetime(scada_end) if start_dt >= end_dt: raise ValueError("SCADA 时间窗非法：scada_start 必须早于 scada_end。") node_query_id = _get_pressure_sensor_mapping(network) missing_nodes = [node_id for node_id in sensor_nodes if node_id not in node_query_id] if missing_nodes: preview = ", ".join(missing_nodes[:10]) raise ValueError(f"未找到可用于压力观测的 SCADA api_query_id: {preview}") query_ids = [node_query_id[node_id] for node_id in sensor_nodes] scada_data = InternalQueries.query_scada_by_ids_timerange( db_name=network, device_ids=query_ids, start_time=start_dt.isoformat(), end_time=end_dt.isoformat(), ) available_lengths = [ len(scada_data.get(query_id, [])) for query_id in query_ids if len(scada_data.get(query_id, [])) > 0 ] if not available_lengths: raise ValueError("指定时间窗内未查询到压力 SCADA 数据。") min_len = min(available_lengths) observation_df = pd.DataFrame() for node_id in sensor_nodes: query_id = node_query_id[node_id] records = scada_data.get(query_id, [])[:min_len] if len(records) < min_len: continue observation_df[node_id] = [float(item["value"]) for item in records] if observation_df.empty: raise ValueError("SCADA 压力数据无法构建观测矩阵。") return observation_df def _to_datetime(value: datetime | str) -> datetime: return parse_utc_time(value)