"""Pressure sensor placement based on scalable sensitivity analysis.

The original implementation expanded a sparse water network into several dense
``node x node``, ``node x pipe``, and ``pipe x pipe`` matrices.  That made the
memory requirement quadratic and the explicit matrix inverse cubic in time.

This module keeps one algorithm for every network size:

* run EPANET once and reuse the first hydraulic state;
* keep incidence and hydraulic graphs sparse;
* estimate the row-wise L1 pressure sensitivity with deterministic Cauchy
  projections and one sparse factorization;
* estimate total directed hydraulic distance from a deterministic spatial
  coreset without materialising an all-pairs distance matrix;
* balance sensitivity score with geographic and pipe-network coverage without
  materialising candidate-to-candidate distances.

The random seed and sample counts are fixed, so the same model and request
produce the same placement on every run.
"""

from __future__ import annotations

import logging
from dataclasses import dataclass
from pathlib import Path
from tempfile import TemporaryDirectory
from time import perf_counter

import numpy as np
import wntr
from scipy.sparse import csr_matrix, eye
from scipy.sparse.csgraph import connected_components, dijkstra
from scipy.sparse.linalg import splu
from sklearn.cluster import MiniBatchKMeans


logger = logging.getLogger(__name__)

_RANDOM_SEED = 42
_SENSITIVITY_PROJECTIONS = 256
_HYDRAULIC_LANDMARKS = 256
_PROJECTION_BLOCK_SIZE = 16
_DIJKSTRA_BLOCK_SIZE = 16
_HEADLOSS_EPSILON = 1e-10
_DIAMETER_TOLERANCE_MM = 1e-9
_COVERAGE_ELIGIBILITY_RATIO = 0.70
_COVERAGE_EDGE_EPSILON = 1e-9


@dataclass(frozen=True)
class _PreparedNetwork:
    """Sparse data required by the placement pipeline."""

    node_names: tuple[str, ...]
    full_node_indices: np.ndarray
    candidate_indices: np.ndarray
    coordinates: np.ndarray
    incidence: csr_matrix
    conductance: np.ndarray
    roughness_response: np.ndarray
    distance_graph: csr_matrix
    coverage_graph: csr_matrix


@dataclass(frozen=True)
class _CandidatePool:
    """Aligned candidate arrays consumed by the placement stage."""

    full_indices: np.ndarray
    coordinates: np.ndarray
    names: np.ndarray
    scores: np.ndarray


def _run_hydraulic_simulation(
    wn: wntr.network.WaterNetworkModel,
):
    """Run only the initial EPANET state without shared ``temp.*`` files."""

    original_duration = wn.options.time.duration
    try:
        # Every downstream calculation reads ``iloc[0]``.  Running an extended
        # simulation only allocates unused time-series results, which is
        # especially expensive for daily models with tens of thousands of
        # nodes.  Restore the caller's model even when EPANET fails.
        wn.options.time.duration = 0
        with TemporaryDirectory(prefix="tjwater-sensitivity-") as temp_dir:
            file_prefix = str(Path(temp_dir) / "simulation")
            return wntr.sim.EpanetSimulator(wn).run_sim(file_prefix=file_prefix)
    finally:
        wn.options.time.duration = original_duration


def _excluded_elements(
    wn: wntr.network.WaterNetworkModel,
) -> tuple[set[str], set[str]]:
    """Return nodes that cannot host sensors and source-connected pipes.

    Reservoirs, tanks, pump/valve endpoints, and the junction immediately next
    to a reservoir or tank are treated as hydraulic boundary nodes.  Pipes
    connected directly to a source are removed from the perturbation set, as
    in the legacy algorithm.
    """

    source_nodes = set(wn.reservoir_name_list) | set(wn.tank_name_list)
    excluded_nodes = set(source_nodes)
    source_pipes: set[str] = set()

    for pipe_name, pipe in wn.pipes():
        endpoints = {pipe.start_node_name, pipe.end_node_name}
        if endpoints & source_nodes:
            source_pipes.add(pipe_name)
            excluded_nodes.update(endpoints)

    for _link_name, link in list(wn.pumps()) + list(wn.valves()):
        excluded_nodes.add(link.start_node_name)
        excluded_nodes.add(link.end_node_name)

    return excluded_nodes, source_pipes


def _minimum_weight_csr(
    rows: list[int],
    columns: list[int],
    weights: list[float],
    *,
    shape: tuple[int, int],
) -> csr_matrix:
    """Build a CSR graph while retaining the lightest parallel edge."""

    if not rows:
        return csr_matrix(shape, dtype=np.float64)

    row_array = np.asarray(rows, dtype=np.int64)
    column_array = np.asarray(columns, dtype=np.int64)
    weight_array = np.asarray(weights, dtype=np.float64)
    order = np.lexsort((column_array, row_array))
    row_array = row_array[order]
    column_array = column_array[order]
    weight_array = weight_array[order]

    group_start = np.empty(len(row_array), dtype=bool)
    group_start[0] = True
    group_start[1:] = (row_array[1:] != row_array[:-1]) | (
        column_array[1:] != column_array[:-1]
    )
    starts = np.flatnonzero(group_start)
    minimum_weights = np.minimum.reduceat(weight_array, starts)
    return csr_matrix(
        (minimum_weights, (row_array[starts], column_array[starts])),
        shape=shape,
    )


def _node_coordinates(
    wn: wntr.network.WaterNetworkModel,
    node_names: tuple[str, ...],
) -> np.ndarray:
    coordinate_series = wn.query_node_attribute("coordinates")
    coordinates = np.asarray(
        [coordinate_series.loc[node_name] for node_name in node_names],
        dtype=np.float64,
    )
    if coordinates.ndim != 2 or coordinates.shape[1] < 2:
        raise ValueError("管网节点缺少二维坐标，无法进行监测点空间布置")
    coordinates = coordinates[:, :2]
    if not np.isfinite(coordinates).all():
        raise ValueError("管网节点坐标包含非有限值，无法进行监测点空间布置")
    return coordinates


def _build_coverage_graph(
    wn: wntr.network.WaterNetworkModel,
    results,
    full_node_index: dict[str, int],
) -> csr_matrix:
    """Build the active undirected physical graph used to spread sensors."""

    status_series = results.link["status"].iloc[0]
    rows: list[int] = []
    columns: list[int] = []
    weights: list[float] = []

    for link_name, link in wn.links():
        if float(status_series.loc[link_name]) <= 0:
            continue

        start = full_node_index[link.start_node_name]
        end = full_node_index[link.end_node_name]
        # Pipes carry their physical length.  Pumps and valves are point
        # devices, so a tiny positive length preserves connectivity without
        # dominating shortest-path distance.
        weight = max(
            float(getattr(link, "length", 0.0)),
            _COVERAGE_EDGE_EPSILON,
        )
        rows.extend((start, end))
        columns.extend((end, start))
        weights.extend((weight, weight))

    return _minimum_weight_csr(
        rows,
        columns,
        weights,
        shape=(len(full_node_index), len(full_node_index)),
    )


def _prepare_network(
    wn: wntr.network.WaterNetworkModel,
    results,
    *,
    min_diameter: int,
) -> _PreparedNetwork:
    excluded_nodes, source_pipes = _excluded_elements(wn)
    full_node_names = tuple(wn.node_name_list)
    full_node_index = {
        node_name: index for index, node_name in enumerate(full_node_names)
    }
    node_names = tuple(
        node_name for node_name in full_node_names if node_name not in excluded_nodes
    )
    if not node_names:
        raise ValueError("管网中没有可参与灵敏度分析的节点")

    node_index = {node_name: index for index, node_name in enumerate(node_names)}
    full_node_indices = np.asarray(
        [full_node_index[node_name] for node_name in node_names],
        dtype=np.int64,
    )
    coordinates = _node_coordinates(wn, node_names)

    flow_series = results.link["flowrate"].iloc[0]
    headloss_series = results.link["headloss"].iloc[0]
    head_series = results.node["head"].iloc[0]

    candidate_nodes: set[str] = set()
    for _pipe_name, pipe in wn.pipes():
        diameter_mm = float(pipe.diameter) * 1000.0
        if diameter_mm + _DIAMETER_TOLERANCE_MM < min_diameter:
            continue
        if pipe.start_node_name in node_index:
            candidate_nodes.add(pipe.start_node_name)
        if pipe.end_node_name in node_index:
            candidate_nodes.add(pipe.end_node_name)

    incidence_rows: list[int] = []
    incidence_columns: list[int] = []
    incidence_values: list[float] = []
    conductance: list[float] = []
    roughness_response: list[float] = []
    distance_rows: list[int] = []
    distance_columns: list[int] = []
    distance_weights: list[float] = []

    kept_pipe_count = 0
    for pipe_name, pipe in wn.pipes():
        if pipe_name in source_pipes:
            continue

        start_name = pipe.start_node_name
        end_name = pipe.end_node_name
        if start_name not in node_index and end_name not in node_index:
            continue

        flow = float(flow_series.loc[pipe_name])
        absolute_flow = abs(flow)
        headloss = abs(float(headloss_series.loc[pipe_name]))
        roughness = float(pipe.roughness)
        if roughness <= 0:
            raise ValueError(f"管道 {pipe_name} 的粗糙度必须大于 0")

        orientation = -1.0 if flow < 0 else 1.0
        if start_name in node_index:
            incidence_rows.append(node_index[start_name])
            incidence_columns.append(kept_pipe_count)
            incidence_values.append(-orientation)
        if end_name in node_index:
            incidence_rows.append(node_index[end_name])
            incidence_columns.append(kept_pipe_count)
            incidence_values.append(orientation)

        conductance.append(
            absolute_flow / (1.852 * headloss + _HEADLOSS_EPSILON)
        )
        roughness_response.append(absolute_flow / roughness)

        if flow > 0:
            upstream_name, downstream_name = start_name, end_name
        else:
            upstream_name, downstream_name = end_name, start_name
        hydraulic_weight = (
            abs(float(head_series.loc[start_name]) - float(head_series.loc[end_name]))
            * float(pipe.length)
        )
        distance_rows.append(full_node_index[upstream_name])
        distance_columns.append(full_node_index[downstream_name])
        distance_weights.append(hydraulic_weight)
        kept_pipe_count += 1

    if kept_pipe_count == 0:
        raise ValueError("管网中没有可用于灵敏度分析的管道")

    incidence = csr_matrix(
        (
            np.asarray(incidence_values, dtype=np.float64),
            (
                np.asarray(incidence_rows, dtype=np.int64),
                np.asarray(incidence_columns, dtype=np.int64),
            ),
        ),
        shape=(len(node_names), kept_pipe_count),
    )
    conductance_array = np.asarray(conductance, dtype=np.float64)
    response_array = np.asarray(roughness_response, dtype=np.float64)
    if not np.isfinite(conductance_array).all() or not np.isfinite(
        response_array
    ).all():
        raise ValueError("水力结果产生了非有限灵敏度系数")

    distance_graph = _minimum_weight_csr(
        distance_rows,
        distance_columns,
        distance_weights,
        shape=(len(full_node_names), len(full_node_names)),
    )
    coverage_graph = _build_coverage_graph(wn, results, full_node_index)
    candidate_indices = np.asarray(
        [
            index
            for index, node_name in enumerate(node_names)
            if node_name in candidate_nodes
        ],
        dtype=np.int64,
    )

    return _PreparedNetwork(
        node_names=node_names,
        full_node_indices=full_node_indices,
        candidate_indices=candidate_indices,
        coordinates=coordinates,
        incidence=incidence,
        conductance=conductance_array,
        roughness_response=response_array,
        distance_graph=distance_graph,
        coverage_graph=coverage_graph,
    )


def _axis_normalized_coordinates(coordinates: np.ndarray) -> np.ndarray:
    """Scale each axis independently for MiniBatchKMeans."""

    minimum = coordinates.min(axis=0)
    span = np.ptp(coordinates, axis=0)
    span[span == 0] = 1.0
    return (coordinates - minimum) / span


def _isotropic_coordinates(coordinates: np.ndarray) -> np.ndarray:
    """Normalize coordinates without distorting the network aspect ratio."""

    minimum = coordinates.min(axis=0)
    scale = float(np.max(np.ptp(coordinates, axis=0), initial=0.0))
    if scale == 0:
        scale = 1.0
    return (coordinates - minimum) / scale


def _cluster_labels(
    coordinates: np.ndarray,
    cluster_count: int,
    *,
    random_seed: int,
) -> tuple[np.ndarray, np.ndarray]:
    """Cluster coordinates deterministically with one implementation at all sizes."""

    normalized = _axis_normalized_coordinates(coordinates)
    if cluster_count == 1:
        return np.zeros(len(coordinates), dtype=np.int64), normalized[[0]]
    if cluster_count >= len(coordinates):
        return np.arange(len(coordinates), dtype=np.int64), normalized.copy()

    model = MiniBatchKMeans(
        n_clusters=cluster_count,
        random_state=random_seed,
        n_init=3,
        batch_size=min(len(coordinates), max(1024, cluster_count * 3)),
        max_iter=100,
        max_no_improvement=20,
        reassignment_ratio=0.0,
    )
    labels = model.fit_predict(normalized).astype(np.int64, copy=False)
    return labels, np.asarray(model.cluster_centers_, dtype=np.float64)


def _estimate_log_pressure_sensitivity(prepared: _PreparedNetwork) -> np.ndarray:
    """Estimate each row's L1 sensitivity using streaming Cauchy projections."""

    weighted_incidence = prepared.incidence.multiply(prepared.conductance)
    laplacian = (weighted_incidence @ prepared.incidence.T).tocsc()
    diagonal = np.asarray(laplacian.diagonal(), dtype=np.float64)
    diagonal_scale = float(np.max(np.abs(diagonal), initial=0.0))
    if diagonal_scale == 0:
        raise ValueError("水力雅可比矩阵为空，无法计算压力灵敏度")

    regularization = diagonal_scale * np.sqrt(np.finfo(np.float64).eps)
    laplacian = laplacian + eye(
        laplacian.shape[0], format="csc", dtype=np.float64
    ) * regularization
    factor = splu(
        laplacian,
        permc_spec="MMD_AT_PLUS_A",
        diag_pivot_thresh=0.0,
        options={"SymmetricMode": True},
    )

    random = np.random.default_rng(_RANDOM_SEED)
    log_absolute_sum = np.zeros(len(prepared.node_names), dtype=np.float64)
    projection_count = 0
    float_epsilon = np.finfo(np.float64).eps
    float_tiny = np.finfo(np.float64).tiny

    while projection_count < _SENSITIVITY_PROJECTIONS:
        block_size = min(
            _PROJECTION_BLOCK_SIZE,
            _SENSITIVITY_PROJECTIONS - projection_count,
        )
        uniform = random.random((prepared.incidence.shape[1], block_size))
        np.clip(uniform, float_epsilon, 1.0 - float_epsilon, out=uniform)
        cauchy_projection = np.tan(np.pi * (uniform - 0.5))
        projected_response = prepared.incidence @ (
            prepared.roughness_response[:, None] * cauchy_projection
        )
        solution = factor.solve(np.asarray(projected_response, dtype=np.float64))
        log_absolute_sum += np.log(
            np.maximum(np.abs(solution), float_tiny)
        ).sum(axis=1)
        projection_count += block_size

    # For a standard Cauchy variable E[log(abs(X))] is zero.  Therefore this
    # streaming geometric mean estimates log(||row||_1) without retaining the
    # node-by-projection matrix.  A finite-sample bias is common to all rows and
    # does not affect ranking.
    return log_absolute_sum / _SENSITIVITY_PROJECTIONS


def _landmark_coreset(prepared: _PreparedNetwork) -> tuple[np.ndarray, np.ndarray]:
    landmark_count = min(_HYDRAULIC_LANDMARKS, len(prepared.node_names))
    labels, centers = _cluster_labels(
        prepared.coordinates,
        landmark_count,
        random_seed=_RANDOM_SEED + 1,
    )
    normalized = _axis_normalized_coordinates(prepared.coordinates)
    landmarks: list[int] = []
    weights: list[float] = []

    for label in np.unique(labels):
        members = np.flatnonzero(labels == label)
        center = centers[int(label)]
        squared_distance = np.square(normalized[members] - center).sum(axis=1)
        landmarks.append(int(members[int(np.argmin(squared_distance))]))
        weights.append(float(len(members)))

    return (
        np.asarray(landmarks, dtype=np.int64),
        np.asarray(weights, dtype=np.float64),
    )


def _estimate_hydraulic_distance_sums(prepared: _PreparedNetwork) -> np.ndarray:
    """Estimate outbound distance sums without an all-pairs distance matrix."""

    landmark_indices, landmark_weights = _landmark_coreset(prepared)
    full_landmark_indices = prepared.full_node_indices[landmark_indices]
    reversed_graph = prepared.distance_graph.transpose().tocsr()
    distance_sums = np.zeros(len(prepared.node_names), dtype=np.float64)

    for start in range(0, len(landmark_indices), _DIJKSTRA_BLOCK_SIZE):
        stop = min(start + _DIJKSTRA_BLOCK_SIZE, len(landmark_indices))
        distances = dijkstra(
            reversed_graph,
            directed=True,
            indices=full_landmark_indices[start:stop],
            return_predecessors=False,
        )
        distances = np.atleast_2d(distances)[:, prepared.full_node_indices]
        # The legacy matrix represented unreachable pairs as zero.  Retaining
        # that convention prevents disconnected branches from receiving an
        # artificial infinite score.
        distances[~np.isfinite(distances)] = 0.0
        distance_sums += landmark_weights[start:stop] @ distances

    return distance_sums


def _build_candidate_pool(
    prepared: _PreparedNetwork,
    log_sensitivity: np.ndarray,
    hydraulic_distance_sums: np.ndarray,
) -> _CandidatePool:
    candidate_indices = prepared.candidate_indices
    candidate_distance = hydraulic_distance_sums[candidate_indices]
    with np.errstate(divide="ignore", invalid="ignore"):
        scores = log_sensitivity[candidate_indices] + np.log(candidate_distance)
    scores = np.nan_to_num(
        scores,
        nan=-np.inf,
        neginf=-np.inf,
        posinf=np.finfo(np.float64).max,
    )
    return _CandidatePool(
        full_indices=prepared.full_node_indices[candidate_indices],
        coordinates=_isotropic_coordinates(prepared.coordinates[candidate_indices]),
        names=np.asarray(
            [prepared.node_names[index] for index in candidate_indices],
            dtype=str,
        ),
        scores=scores,
    )


def _highest_scoring_position(
    names: np.ndarray,
    scores: np.ndarray,
    positions: np.ndarray,
) -> int:
    """Return the best position, breaking score ties by node name."""

    order = np.lexsort((names[positions], -scores[positions]))
    return int(positions[order[0]])


def _relative_gap(
    distances: np.ndarray,
    available: np.ndarray,
) -> np.ndarray:
    """Normalize available distances to their current finite maximum."""

    maximum = float(np.max(distances[available], initial=0.0))
    if not np.isfinite(maximum) or maximum <= np.finfo(np.float64).eps:
        return np.zeros(len(distances), dtype=np.float64)
    return distances / maximum


def _eligible_gap_positions(
    relative_gap: np.ndarray,
    available: np.ndarray,
) -> np.ndarray:
    """Return positions within the configured fraction of the largest gap."""

    maximum = float(np.max(relative_gap[available], initial=0.0))
    if maximum <= np.finfo(np.float64).eps:
        return np.flatnonzero(available)
    threshold = _COVERAGE_ELIGIBILITY_RATIO * maximum
    return np.flatnonzero(
        available & (relative_gap >= threshold - np.finfo(np.float64).eps)
    )


def _allocate_component_quotas(
    coverage_graph: csr_matrix,
    candidate_full_indices: np.ndarray,
    candidate_scores: np.ndarray,
    *,
    sensor_num: int,
) -> tuple[np.ndarray, np.ndarray]:
    """Allocate sensor counts by active pipe length with candidate caps."""

    component_count, node_components = connected_components(
        coverage_graph,
        directed=False,
        return_labels=True,
    )
    candidate_components = node_components[candidate_full_indices]
    capacities = np.bincount(
        candidate_components,
        minlength=component_count,
    ).astype(np.int64, copy=False)

    # The graph is symmetric.  Summed row weights count every physical edge
    # twice, hence the division by two after aggregation by component.
    node_lengths = np.asarray(coverage_graph.sum(axis=1)).ravel()
    component_lengths = np.bincount(
        node_components,
        weights=node_lengths,
        minlength=component_count,
    ) / 2.0
    component_best_scores = np.full(component_count, -np.inf, dtype=np.float64)
    np.maximum.at(
        component_best_scores,
        candidate_components,
        candidate_scores,
    )

    active_components = np.flatnonzero(capacities)
    quotas = np.zeros(component_count, dtype=np.int64)
    if len(active_components) > sensor_num:
        order = np.lexsort(
            (
                active_components,
                -component_best_scores[active_components],
                -component_lengths[active_components],
            )
        )
        quotas[active_components[order[:sensor_num]]] = 1
        return candidate_components, quotas

    quotas[active_components] = 1
    remaining = sensor_num - len(active_components)
    while remaining > 0:
        available = active_components[
            quotas[active_components] < capacities[active_components]
        ]
        if len(available) == 0:
            raise ValueError("连通区域中的候选节点不足，无法分配监测点名额")

        weights = component_lengths[available]
        if float(weights.sum()) <= 0:
            weights = (capacities[available] - quotas[available]).astype(
                np.float64,
                copy=False,
            )
        ideal = remaining * weights / float(weights.sum())
        whole = np.minimum(
            np.floor(ideal).astype(np.int64),
            capacities[available] - quotas[available],
        )
        whole_count = int(whole.sum())
        if whole_count:
            quotas[available] += whole
            remaining -= whole_count
            continue

        fractional = ideal - np.floor(ideal)
        order = np.lexsort(
            (
                available,
                -component_best_scores[available],
                -weights,
                -fractional,
            )
        )
        for component in available[order]:
            quotas[component] += 1
            remaining -= 1
            if remaining == 0:
                break

    return candidate_components, quotas


def _select_component_positions(
    coverage_graph: csr_matrix,
    candidates: _CandidatePool,
    component_positions: np.ndarray,
    existing_positions: list[int],
    *,
    quota: int,
) -> list[int]:
    """Select one component's sensors with score-aware farthest-first search."""

    local_coordinates = candidates.coordinates[component_positions]
    local_names = candidates.names[component_positions]
    local_scores = candidates.scores[component_positions]
    nearest_geographic = np.full(len(component_positions), np.inf)
    nearest_topological = np.full(len(component_positions), np.inf)

    for position in existing_positions:
        nearest_geographic = np.minimum(
            nearest_geographic,
            np.linalg.norm(
                local_coordinates - candidates.coordinates[position],
                axis=1,
            ),
        )

    if existing_positions:
        all_local = np.ones(len(component_positions), dtype=bool)
        seed_eligible = _eligible_gap_positions(
            _relative_gap(nearest_geographic, all_local),
            all_local,
        )
    else:
        seed_eligible = np.arange(len(component_positions), dtype=np.int64)

    seed = _highest_scoring_position(
        local_names,
        local_scores,
        seed_eligible,
    )
    selected_local = [seed]
    remaining = np.ones(len(component_positions), dtype=bool)
    remaining[seed] = False

    while len(selected_local) < quota:
        newest = selected_local[-1]
        geographic_distance = np.linalg.norm(
            local_coordinates - local_coordinates[newest],
            axis=1,
        )
        nearest_geographic = np.minimum(
            nearest_geographic,
            geographic_distance,
        )

        source = int(candidates.full_indices[component_positions[newest]])
        topological_distance = dijkstra(
            coverage_graph,
            directed=False,
            indices=source,
            return_predecessors=False,
        )[candidates.full_indices[component_positions]]
        nearest_topological = np.minimum(
            nearest_topological,
            topological_distance,
        )

        coverage_gap = np.maximum(
            _relative_gap(nearest_geographic, remaining),
            _relative_gap(nearest_topological, remaining),
        )
        eligible_local = _eligible_gap_positions(coverage_gap, remaining)
        next_local = _highest_scoring_position(
            local_names,
            local_scores,
            eligible_local,
        )
        selected_local.append(next_local)
        remaining[next_local] = False

    return [int(component_positions[position]) for position in selected_local]


def _geographic_coverage_metrics(
    candidate_coordinates: np.ndarray,
    selected_positions: list[int],
) -> tuple[float, float, float]:
    normalized = _isotropic_coordinates(candidate_coordinates)
    nearest = np.full(len(normalized), np.inf)
    for position in selected_positions:
        nearest = np.minimum(
            nearest,
            np.linalg.norm(normalized - normalized[position], axis=1),
        )

    selected_coordinates = normalized[selected_positions]
    if len(selected_positions) < 2:
        minimum_gap = 0.0
    else:
        pairwise = np.linalg.norm(
            selected_coordinates[:, None, :] - selected_coordinates[None, :, :],
            axis=2,
        )
        np.fill_diagonal(pairwise, np.inf)
        minimum_gap = float(pairwise.min())
    return (
        float(nearest.max()),
        float(np.quantile(nearest, 0.95)),
        minimum_gap,
    )


def _select_sensor_nodes(
    prepared: _PreparedNetwork,
    log_sensitivity: np.ndarray,
    hydraulic_distance_sums: np.ndarray,
    *,
    sensor_num: int,
) -> list[str]:
    candidate_indices = prepared.candidate_indices
    if len(candidate_indices) < sensor_num:
        raise ValueError(
            "满足最小管径要求的候选节点少于请求的监测点数量："
            f"候选 {len(candidate_indices)} 个，请求 {sensor_num} 个"
        )

    candidates = _build_candidate_pool(
        prepared,
        log_sensitivity,
        hydraulic_distance_sums,
    )
    candidate_components, component_quotas = _allocate_component_quotas(
        prepared.coverage_graph,
        candidates.full_indices,
        candidates.scores,
        sensor_num=sensor_num,
    )
    selected_positions: list[int] = []
    quota_components = np.flatnonzero(component_quotas)
    component_order = np.lexsort(
        (quota_components, -component_quotas[quota_components])
    )
    for component in quota_components[component_order]:
        component_positions = np.flatnonzero(candidate_components == component)
        selected_positions.extend(
            _select_component_positions(
                prepared.coverage_graph,
                candidates,
                component_positions,
                selected_positions,
                quota=int(component_quotas[component]),
            )
        )

    selected_array = np.asarray(selected_positions, dtype=np.int64)
    selected_order = np.lexsort(
        (
            candidates.names[selected_array],
            -candidates.scores[selected_array],
        )
    )
    selected_positions = selected_array[selected_order].tolist()
    maximum_radius, p95_radius, minimum_gap = _geographic_coverage_metrics(
        prepared.coordinates[candidate_indices],
        selected_positions,
    )
    logger.info(
        "Sensitivity placement coverage: components=%d max_radius=%.6f "
        "p95_radius=%.6f min_sensor_gap=%.6f",
        int(np.count_nonzero(component_quotas)),
        maximum_radius,
        p95_radius,
        minimum_gap,
    )
    return [str(candidates.names[position]) for position in selected_positions]


def optimize_sensor_placement(
    wn: wntr.network.WaterNetworkModel,
    sensor_num: int,
    min_diameter: int,
) -> list[str]:
    """Return deterministic pressure monitoring nodes for a loaded network.

    ``min_diameter`` is expressed in millimetres, matching the HTTP contract.
    A node is a valid installation candidate when at least one incident pipe
    meets the threshold.  All valid hydraulic nodes still participate in the
    sensitivity calculation so small pipes continue to influence the result.
    """

    if sensor_num <= 0:
        raise ValueError("监测点数量必须大于 0")
    if min_diameter < 0:
        raise ValueError("最小管径不能小于 0")

    total_started = perf_counter()
    simulation_started = total_started
    results = _run_hydraulic_simulation(wn)
    simulation_seconds = perf_counter() - simulation_started

    preparation_started = perf_counter()
    prepared = _prepare_network(wn, results, min_diameter=min_diameter)
    preparation_seconds = perf_counter() - preparation_started

    sensitivity_started = perf_counter()
    log_sensitivity = _estimate_log_pressure_sensitivity(prepared)
    sensitivity_seconds = perf_counter() - sensitivity_started

    distance_started = perf_counter()
    hydraulic_distance_sums = _estimate_hydraulic_distance_sums(prepared)
    distance_seconds = perf_counter() - distance_started

    selection_started = perf_counter()
    selected = _select_sensor_nodes(
        prepared,
        log_sensitivity,
        hydraulic_distance_sums,
        sensor_num=sensor_num,
    )
    selection_seconds = perf_counter() - selection_started

    logger.info(
        "Sensitivity placement completed: nodes=%d pipes=%d candidates=%d "
        "sensors=%d seconds=%.3f "
        "(simulation=%.3f preparation=%.3f sensitivity=%.3f "
        "distance=%.3f selection=%.3f)",
        len(prepared.node_names),
        prepared.incidence.shape[1],
        len(prepared.candidate_indices),
        len(selected),
        perf_counter() - total_started,
        simulation_seconds,
        preparation_seconds,
        sensitivity_seconds,
        distance_seconds,
        selection_seconds,
    )
    return selected


def optimize_sensor_placement_from_inp(
    inp_path: str | Path,
    sensor_num: int,
    min_diameter: int,
) -> list[str]:
    """Load an EPANET INP model and run the unified placement algorithm."""

    wn = wntr.network.WaterNetworkModel(str(inp_path))
    return optimize_sensor_placement(
        wn,
        sensor_num=sensor_num,
        min_diameter=min_diameter,
    )