#!/usr/bin/env python3 """Screen alternative burst pipes near existing pressure meters using a frozen INP.""" from __future__ import annotations import argparse import asyncio import heapq import json import math import sys import os from collections import defaultdict, deque from pathlib import Path from uuid import uuid4 def nearest_pressure_distances( graph: dict[str, list[tuple[str, float]]], devices: list[dict], radius_m: float, ) -> dict[str, tuple[float, str]]: """Measure to the nearest pressure meter over the network, not just a selected meter.""" distances = {} pending = [] for device in devices: node, device_id = device["node_id"], device["device_id"] value = (0.0, device_id) if value < distances.get(node, (math.inf, "")): distances[node] = value heapq.heappush(pending, (0.0, device_id, node)) while pending: distance, device_id, node = heapq.heappop(pending) if distances[node] != (distance, device_id): continue for neighbor, length in graph.get(node, []): cost = distance + length value = (cost, device_id) if cost <= radius_m and value < distances.get(neighbor, (math.inf, "")): distances[neighbor] = value heapq.heappush(pending, (cost, device_id, neighbor)) return distances def main() -> int: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--calibration-dir", required=True, type=Path) parser.add_argument("--limit", type=int, default=24) parser.add_argument("--offset", type=int, default=0) parser.add_argument("--radius-m", type=float, default=400) parser.add_argument("--min-network-distance-m", type=float, default=0) parser.add_argument("--min-diameter", type=float, default=150) parser.add_argument("--require-isolation", action="store_true") parser.add_argument("--exclude-existing-km", type=float, default=0) parser.add_argument("--skip-screening", type=Path) args = parser.parse_args() os.umask(0o077) sys.path.insert(0, str(Path(__file__).resolve().parents[1])) from scripts.benchmark.calibration import configure_inp, digest, run_file, save_json from scripts.benchmark.simulation_design import BurstSimulationDesign from scripts.check_burst_simulation_design import readonly, resolve_project if (not 1 <= args.limit <= 100 or args.offset < 0 or not 100 <= args.radius_m <= 5000 or not 0 <= args.min_network_distance_m < args.radius_m or not 50 <= args.min_diameter <= 800 or not 0 <= args.exclude_existing_km <= 10): parser.error("invalid screen limit, offset, radius, diameter or exclusion distance") folder = args.calibration_dir.resolve() calibration = json.loads((folder / "calibration.json").read_text()) design = BurstSimulationDesign.model_validate(calibration["design"]) snapshot = folder / "source.inp" if digest(snapshot) != calibration["model_sha256"]: raise ValueError("model snapshot fingerprint changed") biz, _ = asyncio.run(resolve_project(design)) with readonly(biz) as connection: devices = connection.execute( "SELECT device_id,device_type,node_id,link_id FROM asset.scada_devices " "WHERE device_type IN ('pressure','pipe_flow') ORDER BY device_id" ).fetchall() pressure = [device for device in devices if device["device_type"] == "pressure"] flow = [device for device in devices if device["device_type"] == "pipe_flow" and device["device_id"] != "547063"] source = snapshot.read_text() graph = defaultdict(list) pipes = {} valve_links = {} source_nodes = set() closed_links = set() junctions = set() existing_emitters = set() section = "" for raw in source.splitlines(): line = raw.split(";", 1)[0].strip() if line.startswith("["): section = line.strip("[]").upper() continue tokens = line.split() if not tokens: continue if section == "JUNCTIONS": junctions.add(tokens[0]) elif section in {"RESERVOIRS", "TANKS"}: source_nodes.add(tokens[0]) elif section == "STATUS" and tokens[-1].upper() == "CLOSED": closed_links.add(tokens[0]) elif section == "EMITTERS": existing_emitters.add(tokens[0]) elif section in {"PIPES", "VALVES"}: link, left, right = tokens[:3] length = float(tokens[3]) if section == "PIPES" else 0.01 if section == "PIPES": if tokens[-1].upper() == "CLOSED": continue pipes[link] = {"pipe_id": link, "start_node_id": left, "emitter_node_id": right, "length_m": length, "diameter_mm": float(tokens[4])} else: valve_links[link] = (left, right) graph[left].append((right, max(length, 0.01))) graph[right].append((left, max(length, 0.01))) # Initial link overrides also apply to pipe and valve connectivity. graph.clear() pipe_graph = defaultdict(list) for link, pipe in pipes.items(): if link in closed_links: continue left, right = pipe["start_node_id"], pipe["emitter_node_id"] length = max(pipe["length_m"], 0.01) graph[left].append((right, length)) graph[right].append((left, length)) pipe_graph[left].append(right) pipe_graph[right].append(left) for link, (left, right) in valve_links.items(): if link not in closed_links: graph[left].append((right, 0.01)) graph[right].append((left, 0.01)) component_by_node = {} component_data = {} def isolation(node): if node not in component_by_node: seen = {node} pending = deque([node]) while pending: for neighbor in pipe_graph[pending.popleft()]: if neighbor not in seen: seen.add(neighbor) pending.append(neighbor) component_id = min(seen) for member in seen: component_by_node[member] = component_id boundaries = sorted(link for link, (left, right) in valve_links.items() if link not in closed_links and (left in seen) != (right in seen)) component_data[component_id] = {"region_id": component_id, "node_count": len(seen), "boundary_valve_ids": boundaries, "source_inside": bool(seen & source_nodes)} return component_data[component_by_node[node]] # Only real pressure-meter model nodes present in this frozen snapshot are used. pressure = [device for device in pressure if device["node_id"] in junctions] flow = [device for device in flow if device["link_id"] in pipes] distances = nearest_pressure_distances(graph, pressure, args.radius_m) candidates = [] for pipe in pipes.values(): node = pipe["emitter_node_id"] if (args.min_diameter <= pipe["diameter_mm"] <= 800 and pipe["length_m"] >= 50 and pipe["pipe_id"] not in closed_links and node in junctions and node not in existing_emitters and node in distances): distance, device = distances[node] if distance < args.min_network_distance_m: continue region = isolation(node) if args.require_isolation and (region["source_inside"] or not 1 <= len(region["boundary_valve_ids"]) <= 12): continue candidates.append({**pipe, "nearest_pressure_device": device, "network_distance_m": distance, "screen_priority": pipe["length_m"] / pipe["diameter_mm"] ** 2, "isolation": region}) skipped = set() if args.skip_screening: skipped = {entry["pipe_id"] for entry in json.loads(args.skip_screening.read_text())["results"]} with readonly(biz) as connection: coords = {row["node_id"]: row for row in connection.execute( "SELECT node_id,ST_X(ST_Transform(geom,4326)) longitude," "ST_Y(ST_Transform(geom,4326)) latitude FROM gis.node_geometries WHERE node_id=ANY(%s)", (list({pipe["emitter_node_id"] for pipe in candidates}),), )} existing_positions = connection.execute( "SELECT ST_X(ST_Transform(geom,4326)) longitude,ST_Y(ST_Transform(geom,4326)) latitude " "FROM gis.node_geometries WHERE node_id=ANY(%s)", ([pipe.emitter_node_id for pipe in design.pipes],), ).fetchall() filtered = [] for candidate in candidates: position = coords.get(candidate["emitter_node_id"]) if position is None or candidate["pipe_id"] in skipped: continue distances_km = [111.32 * math.hypot((position["longitude"] - old["longitude"]) * math.cos(math.radians(position["latitude"])), position["latitude"] - old["latitude"]) for old in existing_positions] if distances_km and min(distances_km) < args.exclude_existing_km: continue candidate.update(longitude=position["longitude"], latitude=position["latitude"], nearest_existing_km=min(distances_km) if distances_km else None) filtered.append(candidate) candidates = filtered # Cover different meter neighborhoods before taking additional pipes at the same meter. candidates.sort(key=lambda pipe: pipe["screen_priority"], reverse=True) used = set() unique, remaining = [], [] for candidate in candidates: device = candidate["nearest_pressure_device"] (remaining if device in used else unique).append(candidate) used.add(device) selected = (unique + remaining)[args.offset:args.offset + args.limit] node_ids = sorted({device["node_id"] for device in pressure} | {pipe["emitter_node_id"] for pipe in selected}) link_ids = sorted({device["link_id"] for device in flow}) at = max(design.event_times) output = folder / f"screening-{uuid4().hex[:8]}" output.mkdir(mode=0o700) report = {"version": "burst-pipe-screening@1", "model_sha256": digest(snapshot), "at": at.isoformat(), "candidate_count": len(candidates), "results": [], "screen_options": {key: str(value) if isinstance(value, Path) else value for key, value in vars(args).items()}, "device_bindings": devices, "scope": "daytime single-frame screening only; requires full paired calibration and isolation review"} def simulate(label, emitter=None): path = output / f"{label}.inp" path.write_text(configure_inp(source, at, 0, emitter=emitter, pda_options=calibration["pda_options"])) data = run_file(path, node_ids, link_ids, timeout=240) save_json(output / f"{label}.json", data) return data normal = simulate("normal") if normal["status"] != "completed" or normal["warnings"]: raise ValueError("normal screening baseline failed") for index, pipe in enumerate(selected): area = math.pi * (pipe["diameter_mm"] / 10) ** 2 / 4 / 8 burst = simulate(f"candidate-{index:02d}", (pipe["emitter_node_id"], area)) result = {**pipe, "orifice_area_cm2": area, "status": burst["status"], "output_file": f"candidate-{index:02d}.json"} if burst["status"] == "completed": drops = {device["device_id"]: 100 * (normal["nodes"][device["node_id"]][0]["pressure"] - burst["nodes"][device["node_id"]][0]["pressure"]) / max(abs(normal["nodes"][device["node_id"]][0]["pressure"]), 1e-6) for device in pressure} changes = {device["device_id"]: 100 * abs(normal["links"][device["link_id"]][0]["flow"] - burst["links"][device["link_id"]][0]["flow"]) / abs(normal["links"][device["link_id"]][0]["flow"]) for device in flow if abs(normal["links"][device["link_id"]][0]["flow"]) > 1e-6} ranked = sorted(drops.items(), key=lambda pair: pair[1], reverse=True) negative = any(row[0]["pressure"] < 0 for row in burst["nodes"].values()) result.update(pressure_drop_percent=dict(ranked), flow_change_percent=changes, possible_candidate=not burst["warnings"] and not negative and ranked[0][1] >= 30 and (ranked[1][1] >= 10 or max(changes.values(), default=0) >= 20)) print(f"{index+1}/{len(selected)} {pipe['pipe_id']}: max pressure drop={ranked[0][1]:.2f}%, " f"possible candidate={result['possible_candidate']}", flush=True) report["results"].append(result) save_json(output / "screening.json", report) print(f"screening complete; report: {output / 'screening.json'}", flush=True) return 0 if __name__ == "__main__": raise SystemExit(main())