227 lines
11 KiB
Python
227 lines
11 KiB
Python
"""File-isolated hydraulic calibration; database connections are read-only."""
|
|
from __future__ import annotations
|
|
|
|
import ctypes
|
|
import hashlib
|
|
import json
|
|
import math
|
|
import os
|
|
import re
|
|
import subprocess
|
|
from datetime import datetime, timedelta
|
|
from pathlib import Path
|
|
from uuid import NAMESPACE_URL, uuid5
|
|
from zoneinfo import ZoneInfo
|
|
|
|
from scripts.benchmark.simulation_design import PipeCalibration, SelectionRule
|
|
|
|
|
|
def digest(path: Path) -> str:
|
|
return hashlib.sha256(path.read_bytes()).hexdigest()
|
|
|
|
|
|
def save_json(path: Path, value: dict) -> None:
|
|
temporary = path.with_suffix(".pending")
|
|
with temporary.open("x", encoding="utf-8") as output:
|
|
os.chmod(temporary, 0o600)
|
|
json.dump(value, output, ensure_ascii=False, indent=2, allow_nan=False, default=str)
|
|
output.write("\n")
|
|
temporary.replace(path)
|
|
|
|
|
|
def select_calibrated_area(
|
|
pipe: PipeCalibration, passed_levels: list[str], case_results: list[dict], retain_weak: bool,
|
|
) -> dict:
|
|
selected = passed_levels[-1] if passed_levels else None
|
|
soft_blockers = {"insufficient_consecutive_pressure_and_corroboration", "event_rule_did_not_trigger"}
|
|
physically_valid = bool(case_results) and all(
|
|
not set(case["blockers"]) - soft_blockers for case in case_results
|
|
)
|
|
status = "calibrated_candidate" if selected else (
|
|
"weak_signal_candidate" if retain_weak and physically_valid else "unsuitable"
|
|
)
|
|
return {"status": status, "severity": selected,
|
|
"area_cm2": pipe.candidate_orifice_areas_cm2[selected] if selected else None}
|
|
|
|
|
|
def configure_inp(
|
|
source: str,
|
|
start: datetime,
|
|
duration: int,
|
|
*,
|
|
emitter: tuple[str, float] | None = None,
|
|
pda_options: dict[str, str] | None = None,
|
|
) -> str:
|
|
"""Apply the burst service's options to a file, including the zero-leak control."""
|
|
if start.tzinfo is None:
|
|
raise ValueError("simulation start must include a timezone")
|
|
start = start.astimezone(ZoneInfo("Asia/Shanghai"))
|
|
updates = {
|
|
"TIMES": {"DURATION": f"{duration // 3600}:{duration % 3600 // 60:02d}:00",
|
|
"PATTERN START": start.strftime("%H:%M:%S"),
|
|
"REPORT START": "0:00:00", "STATISTIC": "NONE"},
|
|
"OPTIONS": {**(pda_options or {}), "DEMAND MODEL": "PDA", "REQUIRED PRESSURE": "10.0000"},
|
|
"REPORT": {"NODES": "NONE", "LINKS": "NONE"},
|
|
}
|
|
section = None
|
|
output = []
|
|
pending = {}
|
|
for line in source.splitlines():
|
|
match = re.fullmatch(r"\s*\[([^]]+)\]\s*", line)
|
|
if match:
|
|
output.extend(f"{key} {value}" for key, value in pending.items())
|
|
section = match[1].upper()
|
|
pending = dict(updates.get(section, {}))
|
|
output.append(line)
|
|
if section == "EMITTERS" and emitter:
|
|
node, area = emitter
|
|
coefficient = 0.65 * area / 10000 * math.sqrt(19.6) * 1000
|
|
output.append(f"{node} {coefficient:.12g}")
|
|
continue
|
|
tokens = line.split(";", 1)[0].upper().split()
|
|
if any(tokens[:len(key.split())] == key.split() for key in updates.get(section, {})):
|
|
continue
|
|
output.append(line)
|
|
output.extend(f"{key} {value}" for key, value in pending.items())
|
|
if "[EMITTERS]" not in source.upper() and emitter:
|
|
raise ValueError("exported model has no emitter section")
|
|
return "\n".join(output) + "\n"
|
|
|
|
|
|
def run_file(inp: Path, node_ids: list[str], link_ids: list[str], *, timeout: int) -> dict:
|
|
from app.infra.epanet.epanet import Output
|
|
|
|
executable = Path(__file__).resolve().parents[2] / "app/infra/epanet/linux/runepanet"
|
|
env = os.environ.copy()
|
|
env["LD_LIBRARY_PATH"] = f"{executable.parent}:{env.get('LD_LIBRARY_PATH', '')}"
|
|
report_path, binary = inp.with_suffix(".rpt"), inp.with_suffix(".bin")
|
|
result = subprocess.run(
|
|
[str(executable), str(inp), str(report_path), str(binary)],
|
|
env=env, capture_output=True, text=True, timeout=timeout, check=False,
|
|
)
|
|
report = report_path.read_text(encoding="utf-8", errors="replace") if report_path.exists() else ""
|
|
warnings = [line.strip() for line in report.splitlines()
|
|
if re.search(r"^\s*(?:warning\b|error\s+\d+\b)|system\s+(?:is\s+)?unbalanced|cannot solve hydraulic",
|
|
line, re.IGNORECASE)]
|
|
summary = {"exit_code": result.returncode, "warnings": warnings,
|
|
"inp_sha256": digest(inp), "report_sha256": digest(report_path) if report_path.exists() else None}
|
|
if result.returncode or not binary.exists():
|
|
return {**summary, "status": "failed"}
|
|
reader = Output(str(binary))
|
|
names = reader.element_name()
|
|
summary.update(status="completed", units=reader.units(), times=reader.times(), nodes={}, links={})
|
|
for kind, selected, attributes, function in (
|
|
("nodes", node_ids, ["demand", "head", "pressure", "quality"], reader._lib.ENR_getNodeResult),
|
|
("links", link_ids, ["flow", "velocity", "headloss", "quality", "status", "setting", "reaction", "friction"],
|
|
reader._lib.ENR_getLinkResult),
|
|
):
|
|
indices = {name: index + 1 for index, name in enumerate(names[kind])}
|
|
for name in selected:
|
|
index = indices.get(name)
|
|
if index is None:
|
|
raise ValueError("requested hydraulic output is missing")
|
|
series = []
|
|
for period in range(summary["times"]["num_periods"]):
|
|
values = ctypes.POINTER(ctypes.c_float)()
|
|
length = ctypes.c_int()
|
|
reader._check(function(reader._handle, period, index, ctypes.byref(values), ctypes.byref(length)))
|
|
try:
|
|
if length.value != len(attributes):
|
|
raise ValueError("unexpected hydraulic result attributes")
|
|
row = {key: float(values[i]) for i, key in enumerate(attributes)}
|
|
if not all(math.isfinite(value) for value in row.values()):
|
|
raise ValueError("non-finite hydraulic result")
|
|
series.append(row)
|
|
finally:
|
|
reader._lib.ENR_free(ctypes.byref(values))
|
|
summary[kind][name] = series
|
|
return summary
|
|
|
|
|
|
def pressure_preview(
|
|
pipe: PipeCalibration, start: datetime, duration: int, step: int,
|
|
history: dict, current: dict, rule: SelectionRule,
|
|
) -> dict:
|
|
"""Exercise the actual replay cleaner/detector without database mutations."""
|
|
from app.domain.schemas.scada_replay import ScadaReplayCreate
|
|
from app.services.scada_replay_processing import detect_frame, prepare_cleaning
|
|
|
|
ids = [binding.device_id for binding in pipe.pressure_bindings]
|
|
config = ScadaReplayCreate(
|
|
request_id=uuid5(NAMESPACE_URL, f"calibration:{pipe.pipe_id}:{start.isoformat()}"),
|
|
start_time=start, end_time=start + timedelta(seconds=duration + step), device_ids=ids,
|
|
rules=[{"rule_id": "calibration-pressure", "device_ids": ids,
|
|
"kind": "relative_change", "direction": "decrease", "unit": "m",
|
|
"threshold": rule.primary_pressure_drop_percent,
|
|
"window_samples": 4, "minimum_samples": 4,
|
|
"consecutive_samples": rule.consecutive_samples, "minimum_devices": 1}],
|
|
)
|
|
document = {"replay_id": str(config.request_id), "devices": [
|
|
{"device_id": key, "transmission_frequency": "0:15:00"} for key in ids],
|
|
"histories": {binding.device_id: [row["pressure"] for row in history[binding.node_id]]
|
|
for binding in pipe.pressure_bindings}}
|
|
raised = []
|
|
for index in range(duration // step + 1):
|
|
at = start + timedelta(seconds=index * step)
|
|
rows = [{"device_id": binding.device_id, "time": at,
|
|
"monitored_value": current[binding.node_id][index]["pressure"],
|
|
"source_hash": hashlib.sha256(f"{binding.device_id}:{at}".encode()).hexdigest()}
|
|
for binding in pipe.pressure_bindings]
|
|
document, changes = detect_frame(document, prepare_cleaning(rows, config), config, at)
|
|
raised.extend(change for change in changes
|
|
if change["change"] == "raised" and change.get("category") != "data_quality")
|
|
return {"raised": bool(raised), "first_signal_at": raised[0]["occurred_at"] if raised else None,
|
|
"scope": "pure replay rule preview; not monitoring worker or Agent acceptance"}
|
|
|
|
|
|
def assess(
|
|
pipe: PipeCalibration, normal: dict, burst: dict, rule: SelectionRule,
|
|
*, duration: int, step: int,
|
|
) -> dict:
|
|
expected_times = {"report_start": 0, "report_step": step,
|
|
"sim_duration": duration, "num_periods": duration // step + 1}
|
|
blockers = []
|
|
for label, result in (("normal", normal), ("burst", burst)):
|
|
if result.get("status") != "completed" or result.get("times") != expected_times:
|
|
blockers.append(f"{label}_incomplete_or_failed")
|
|
if result.get("warnings"):
|
|
blockers.append(f"{label}_solver_warning")
|
|
if result.get("units") != {"flow": "LPS", "pressure": "MTR", "quality": "NONE"}:
|
|
blockers.append(f"{label}_unexpected_units")
|
|
if blockers:
|
|
return {"passed": False, "blockers": blockers}
|
|
pressures, flows = {}, {}
|
|
for binding in pipe.pressure_bindings:
|
|
before = [row["pressure"] for row in normal["nodes"][binding.node_id]]
|
|
after = [row["pressure"] for row in burst["nodes"][binding.node_id]]
|
|
if min(before) <= rule.relative_change_epsilon:
|
|
blockers.append(f"nonpositive_normal_pressure:{binding.device_id}")
|
|
if min(after) < 0:
|
|
blockers.append(f"negative_monitor_pressure:{binding.device_id}")
|
|
pressures[binding.device_id] = [100 * (a - b) / max(abs(a), rule.relative_change_epsilon)
|
|
for a, b in zip(before, after)]
|
|
for binding in pipe.flow_bindings:
|
|
flows[binding.device_id] = [100 * abs(a["flow"] - b["flow"]) /
|
|
max(abs(a["flow"]), rule.relative_change_epsilon)
|
|
for a, b in zip(normal["links"][binding.link_id], burst["links"][binding.link_id])]
|
|
witnesses = []
|
|
count = duration // step + 1
|
|
for primary, drops in pressures.items():
|
|
for start in range(count - rule.consecutive_samples + 1):
|
|
window = range(start, start + rule.consecutive_samples)
|
|
if not all(drops[i] >= rule.primary_pressure_drop_percent for i in window):
|
|
continue
|
|
second = [key for key, values in pressures.items() if key != primary and
|
|
all(values[i] >= rule.corroborating_pressure_drop_percent for i in window)]
|
|
flow = [key for key, values in flows.items() if
|
|
all(values[i] >= rule.flow_change_percent for i in window) and
|
|
all(abs(normal["links"][next(b.link_id for b in pipe.flow_bindings if b.device_id == key)][i]["flow"])
|
|
> rule.relative_change_epsilon for i in window)]
|
|
if second or flow:
|
|
witnesses.append({"primary_device": primary, "second_pressure_devices": second,
|
|
"flow_devices": flow, "start_offset_seconds": start * step})
|
|
if not witnesses:
|
|
blockers.append("insufficient_consecutive_pressure_and_corroboration")
|
|
return {"passed": not blockers, "blockers": blockers, "witnesses": witnesses,
|
|
"pressure_drop_percent": pressures, "flow_change_percent": flows}
|