392 lines
12 KiB
Python
392 lines
12 KiB
Python
from pathlib import Path
|
|
|
|
import numpy as np
|
|
import pytest
|
|
import wntr
|
|
from scipy.sparse import csr_matrix, isspmatrix_csr
|
|
from scipy.sparse.csgraph import dijkstra
|
|
|
|
from app.algorithms.sensor import sensitivity
|
|
|
|
|
|
def _build_test_network() -> wntr.network.WaterNetworkModel:
|
|
wn = wntr.network.WaterNetworkModel()
|
|
wn.options.time.duration = 0
|
|
wn.add_reservoir("R1", base_head=100.0, coordinates=(-1.0, 0.0))
|
|
wn.add_junction("J0", elevation=5.0, coordinates=(0.0, 0.0))
|
|
|
|
for index in range(1, 13):
|
|
wn.add_junction(
|
|
f"J{index}",
|
|
base_demand=0.001 + index * 0.00001,
|
|
elevation=5.0 + index * 0.05,
|
|
coordinates=(float(index % 4), float(index // 4)),
|
|
)
|
|
|
|
wn.add_pipe("P0", "R1", "J0", length=100.0, diameter=0.4, roughness=110)
|
|
wn.add_pipe("P1", "J0", "J1", length=100.0, diameter=0.4, roughness=110)
|
|
for index in range(1, 11):
|
|
wn.add_pipe(
|
|
f"P{index + 1}",
|
|
f"J{index}",
|
|
f"J{index + 1}",
|
|
length=80.0 + index,
|
|
diameter=0.3,
|
|
roughness=105,
|
|
)
|
|
# J12 is connected only through a small pipe, while J11 also touches P11.
|
|
wn.add_pipe("P12", "J11", "J12", length=90.0, diameter=0.1, roughness=105)
|
|
wn.add_pipe("PX1", "J2", "J6", length=120.0, diameter=0.3, roughness=105)
|
|
wn.add_pipe("PX2", "J5", "J9", length=120.0, diameter=0.3, roughness=105)
|
|
return wn
|
|
|
|
|
|
def _prepared_selection_network(
|
|
coordinates: np.ndarray,
|
|
edges: list[tuple[int, int, float]],
|
|
) -> sensitivity._PreparedNetwork:
|
|
node_count = len(coordinates)
|
|
rows: list[int] = []
|
|
columns: list[int] = []
|
|
weights: list[float] = []
|
|
for start, end, weight in edges:
|
|
rows.extend((start, end))
|
|
columns.extend((end, start))
|
|
weights.extend((weight, weight))
|
|
coverage_graph = csr_matrix(
|
|
(weights, (rows, columns)),
|
|
shape=(node_count, node_count),
|
|
)
|
|
return sensitivity._PreparedNetwork(
|
|
node_names=tuple(f"N{index:04d}" for index in range(node_count)),
|
|
full_node_indices=np.arange(node_count, dtype=np.int64),
|
|
candidate_indices=np.arange(node_count, dtype=np.int64),
|
|
coordinates=np.asarray(coordinates, dtype=np.float64),
|
|
incidence=csr_matrix((node_count, 1), dtype=np.float64),
|
|
conductance=np.ones(1, dtype=np.float64),
|
|
roughness_response=np.ones(1, dtype=np.float64),
|
|
distance_graph=csr_matrix((node_count, node_count), dtype=np.float64),
|
|
coverage_graph=coverage_graph,
|
|
)
|
|
|
|
|
|
def test_algorithm_is_deterministic_and_runs_epanet_once(monkeypatch, tmp_path):
|
|
wn = _build_test_network()
|
|
original_run_sim = wntr.sim.EpanetSimulator.run_sim
|
|
prefixes: list[str] = []
|
|
|
|
def counted_run_sim(simulator, *args, **kwargs):
|
|
prefixes.append(str(kwargs["file_prefix"]))
|
|
return original_run_sim(simulator, *args, **kwargs)
|
|
|
|
monkeypatch.setattr(wntr.sim.EpanetSimulator, "run_sim", counted_run_sim)
|
|
monkeypatch.chdir(tmp_path)
|
|
|
|
first = sensitivity.optimize_sensor_placement(wn, sensor_num=4, min_diameter=0)
|
|
second = sensitivity.optimize_sensor_placement(wn, sensor_num=4, min_diameter=0)
|
|
|
|
assert first == second
|
|
assert len(first) == len(set(first)) == 4
|
|
assert len(prefixes) == 2
|
|
assert all(not Path(prefix).parent.exists() for prefix in prefixes)
|
|
assert not list(tmp_path.glob("temp.*"))
|
|
|
|
|
|
def test_hydraulic_simulation_keeps_only_initial_state_and_restores_duration():
|
|
wn = _build_test_network()
|
|
wn.options.time.duration = 24 * 60 * 60
|
|
|
|
results = sensitivity._run_hydraulic_simulation(wn)
|
|
|
|
assert len(results.node["head"].index) == 1
|
|
assert wn.options.time.duration == 24 * 60 * 60
|
|
|
|
|
|
def test_preparation_keeps_network_matrices_sparse():
|
|
wn = _build_test_network()
|
|
results = sensitivity._run_hydraulic_simulation(wn)
|
|
|
|
prepared = sensitivity._prepare_network(wn, results, min_diameter=0)
|
|
|
|
assert isspmatrix_csr(prepared.incidence)
|
|
assert isspmatrix_csr(prepared.distance_graph)
|
|
assert isspmatrix_csr(prepared.coverage_graph)
|
|
assert prepared.incidence.nnz <= 2 * prepared.incidence.shape[1]
|
|
assert prepared.distance_graph.nnz <= wn.num_pipes
|
|
assert prepared.coverage_graph.nnz <= 2 * wn.num_links
|
|
assert (prepared.coverage_graph != prepared.coverage_graph.T).nnz == 0
|
|
dense_incidence_bytes = int(np.prod(prepared.incidence.shape)) * 8
|
|
sparse_payload_bytes = (
|
|
prepared.incidence.data.nbytes
|
|
+ prepared.incidence.indices.nbytes
|
|
+ prepared.incidence.indptr.nbytes
|
|
)
|
|
assert sparse_payload_bytes < dense_incidence_bytes
|
|
|
|
|
|
def test_minimum_diameter_filters_installation_candidates_in_millimetres():
|
|
wn = _build_test_network()
|
|
results = sensitivity._run_hydraulic_simulation(wn)
|
|
prepared = sensitivity._prepare_network(wn, results, min_diameter=300)
|
|
candidate_names = {
|
|
prepared.node_names[index] for index in prepared.candidate_indices
|
|
}
|
|
|
|
assert "J12" not in candidate_names
|
|
assert "J11" in candidate_names
|
|
|
|
selected = sensitivity.optimize_sensor_placement(
|
|
wn,
|
|
sensor_num=4,
|
|
min_diameter=300,
|
|
)
|
|
assert set(selected) <= candidate_names
|
|
|
|
with pytest.raises(ValueError, match="候选节点少于"):
|
|
sensitivity.optimize_sensor_placement(
|
|
wn,
|
|
sensor_num=len(candidate_names) + 1,
|
|
min_diameter=300,
|
|
)
|
|
|
|
|
|
def test_sparse_estimate_preserves_dense_reference_placement_quality():
|
|
wn = _build_test_network()
|
|
results = sensitivity._run_hydraulic_simulation(wn)
|
|
prepared = sensitivity._prepare_network(wn, results, min_diameter=0)
|
|
|
|
approximate_log_sensitivity = sensitivity._estimate_log_pressure_sensitivity(
|
|
prepared
|
|
)
|
|
approximate_distance = sensitivity._estimate_hydraulic_distance_sums(prepared)
|
|
approximate_selected = sensitivity._select_sensor_nodes(
|
|
prepared,
|
|
approximate_log_sensitivity,
|
|
approximate_distance,
|
|
sensor_num=4,
|
|
)
|
|
|
|
incidence = prepared.incidence.toarray()
|
|
laplacian = (
|
|
prepared.incidence.multiply(prepared.conductance)
|
|
@ prepared.incidence.T
|
|
).toarray()
|
|
diagonal_scale = float(np.max(np.abs(np.diag(laplacian))))
|
|
laplacian += np.eye(laplacian.shape[0]) * (
|
|
diagonal_scale * np.sqrt(np.finfo(np.float64).eps)
|
|
)
|
|
response = np.linalg.solve(
|
|
laplacian,
|
|
incidence * prepared.roughness_response,
|
|
)
|
|
exact_sensitivity = np.abs(response).sum(axis=1)
|
|
|
|
exact_distances = dijkstra(
|
|
prepared.distance_graph.transpose().tocsr(),
|
|
directed=True,
|
|
indices=prepared.full_node_indices,
|
|
return_predecessors=False,
|
|
)[:, prepared.full_node_indices]
|
|
exact_distances[~np.isfinite(exact_distances)] = 0.0
|
|
exact_distance = exact_distances.sum(axis=0)
|
|
exact_selected = sensitivity._select_sensor_nodes(
|
|
prepared,
|
|
np.log(np.maximum(exact_sensitivity, np.finfo(np.float64).tiny)),
|
|
exact_distance,
|
|
sensor_num=4,
|
|
)
|
|
|
|
exact_score = exact_sensitivity * exact_distance
|
|
node_index = {
|
|
node_name: index for index, node_name in enumerate(prepared.node_names)
|
|
}
|
|
approximate_objective = sum(
|
|
exact_score[node_index[node_name]] for node_name in approximate_selected
|
|
)
|
|
exact_objective = sum(
|
|
exact_score[node_index[node_name]] for node_name in exact_selected
|
|
)
|
|
|
|
assert approximate_objective / exact_objective >= 0.95
|
|
|
|
|
|
def test_mixed_coverage_avoids_candidate_density_bias():
|
|
dense_west = np.linspace(0.0, 2.0, 200)
|
|
sparse_east = np.linspace(3.0, 10.0, 20)
|
|
x_coordinates = np.concatenate((dense_west, sparse_east))
|
|
coordinates = np.column_stack(
|
|
(x_coordinates, np.zeros(len(x_coordinates), dtype=np.float64))
|
|
)
|
|
ordered = np.argsort(x_coordinates)
|
|
edges = [
|
|
(
|
|
int(start),
|
|
int(end),
|
|
float(x_coordinates[end] - x_coordinates[start]),
|
|
)
|
|
for start, end in zip(ordered[:-1], ordered[1:])
|
|
]
|
|
prepared = _prepared_selection_network(coordinates, edges)
|
|
log_scores = np.linspace(4.0, 0.0, len(coordinates))
|
|
distance_sums = np.ones(len(coordinates), dtype=np.float64)
|
|
|
|
selected = sensitivity._select_sensor_nodes(
|
|
prepared,
|
|
log_scores,
|
|
distance_sums,
|
|
sensor_num=6,
|
|
)
|
|
name_to_position = {
|
|
name: position for position, name in enumerate(prepared.node_names)
|
|
}
|
|
selected_positions = [name_to_position[name] for name in selected]
|
|
new_metrics = sensitivity._geographic_coverage_metrics(
|
|
coordinates,
|
|
selected_positions,
|
|
)
|
|
|
|
legacy_labels, _centers = sensitivity._cluster_labels(
|
|
coordinates,
|
|
6,
|
|
random_seed=sensitivity._RANDOM_SEED + 2,
|
|
)
|
|
legacy_positions: list[int] = []
|
|
represented: set[int] = set()
|
|
for position in np.argsort(-log_scores):
|
|
label = int(legacy_labels[position])
|
|
if label in represented:
|
|
continue
|
|
represented.add(label)
|
|
legacy_positions.append(int(position))
|
|
legacy_metrics = sensitivity._geographic_coverage_metrics(
|
|
coordinates,
|
|
legacy_positions,
|
|
)
|
|
|
|
assert new_metrics[0] <= legacy_metrics[0] * 0.6
|
|
assert new_metrics[2] >= legacy_metrics[2] * 1.5
|
|
assert max(x_coordinates[selected_positions]) >= 9.0
|
|
|
|
|
|
def test_disconnected_components_each_receive_a_sensor_when_slots_allow():
|
|
coordinates = np.asarray(
|
|
[
|
|
(0.0, 0.0),
|
|
(1.0, 0.0),
|
|
(0.0, 0.01),
|
|
(1.0, 0.01),
|
|
]
|
|
)
|
|
prepared = _prepared_selection_network(
|
|
coordinates,
|
|
[(0, 1, 1.0), (2, 3, 1.0)],
|
|
)
|
|
|
|
selected = sensitivity._select_sensor_nodes(
|
|
prepared,
|
|
np.asarray([10.0, 9.0, 8.0, 7.0]),
|
|
np.ones(4),
|
|
sensor_num=2,
|
|
)
|
|
|
|
assert len(set(selected) & {"N0000", "N0001"}) == 1
|
|
assert len(set(selected) & {"N0002", "N0003"}) == 1
|
|
|
|
|
|
def test_overlapping_components_respect_global_geographic_spacing():
|
|
coordinates = np.asarray(
|
|
[
|
|
(0.0, 0.0),
|
|
(10.0, 0.0),
|
|
(0.1, 0.0),
|
|
(10.1, 0.0),
|
|
]
|
|
)
|
|
prepared = _prepared_selection_network(
|
|
coordinates,
|
|
[(0, 1, 10.0), (2, 3, 10.0)],
|
|
)
|
|
|
|
selected = sensitivity._select_sensor_nodes(
|
|
prepared,
|
|
np.asarray([10.0, 1.0, 9.0, 0.0]),
|
|
np.ones(4),
|
|
sensor_num=2,
|
|
)
|
|
selected_positions = [prepared.node_names.index(name) for name in selected]
|
|
minimum_gap = sensitivity._geographic_coverage_metrics(
|
|
coordinates,
|
|
selected_positions,
|
|
)[2]
|
|
|
|
assert minimum_gap >= 0.9
|
|
|
|
|
|
def test_component_quota_prefers_longer_networks_when_slots_are_limited():
|
|
coordinates = np.asarray(
|
|
[
|
|
(0.0, 0.0),
|
|
(10.0, 0.0),
|
|
(20.0, 0.0),
|
|
(25.0, 0.0),
|
|
(30.0, 0.0),
|
|
(31.0, 0.0),
|
|
]
|
|
)
|
|
prepared = _prepared_selection_network(
|
|
coordinates,
|
|
[(0, 1, 10.0), (2, 3, 5.0), (4, 5, 1.0)],
|
|
)
|
|
|
|
selected = sensitivity._select_sensor_nodes(
|
|
prepared,
|
|
np.asarray([1.0, 1.0, 2.0, 2.0, 100.0, 100.0]),
|
|
np.ones(6),
|
|
sensor_num=2,
|
|
)
|
|
|
|
assert set(selected) <= {"N0000", "N0001", "N0002", "N0003"}
|
|
assert len(set(selected) & {"N0000", "N0001"}) == 1
|
|
assert len(set(selected) & {"N0002", "N0003"}) == 1
|
|
|
|
|
|
def test_duplicate_coordinates_use_topology_and_return_exact_count():
|
|
coordinates = np.zeros((6, 2), dtype=np.float64)
|
|
prepared = _prepared_selection_network(
|
|
coordinates,
|
|
[(index, index + 1, 1.0) for index in range(5)],
|
|
)
|
|
log_scores = np.linspace(6.0, 1.0, 6)
|
|
|
|
first = sensitivity._select_sensor_nodes(
|
|
prepared,
|
|
log_scores,
|
|
np.ones(6),
|
|
sensor_num=4,
|
|
)
|
|
second = sensitivity._select_sensor_nodes(
|
|
prepared,
|
|
log_scores,
|
|
np.ones(6),
|
|
sensor_num=4,
|
|
)
|
|
|
|
assert first == second
|
|
assert len(first) == len(set(first)) == 4
|
|
|
|
|
|
@pytest.mark.parametrize(
|
|
("sensor_num", "min_diameter", "message"),
|
|
[
|
|
(0, 0, "监测点数量必须大于 0"),
|
|
(1, -1, "最小管径不能小于 0"),
|
|
],
|
|
)
|
|
def test_algorithm_rejects_invalid_parameters(sensor_num, min_diameter, message):
|
|
with pytest.raises(ValueError, match=message):
|
|
sensitivity.optimize_sensor_placement(
|
|
_build_test_network(),
|
|
sensor_num=sensor_num,
|
|
min_diameter=min_diameter,
|
|
)
|