Source code for epyt_flow.data.benchmarks.leakdb

  1"""
  2LeakDB (Leakage Diagnosis Benchmark) by Vrachimis, S. G., Kyriakou, M. S., Eliades, D. G.,
  3and Polycarpou, M. M. (2018), is a realistic leakage dataset for water distribution networks.
  4The dataset is comprised of 1000 artificially created but realistic leakage
  5scenarios, on different water distribution networks, under varying conditions.
  6
  7See https://github.com/KIOS-Research/LeakDB/ for details.
  8
  9This module provides functions for loading the original LeakDB data set
 10:func:`~epyt_flow.data.benchmarks.leakdb.load_data`, as well as methods for loading the scenarios
 11:func:`~epyt_flow.data.benchmarks.leakdb.load_scenarios` and pre-generated SCADA data
 12:func:`~epyt_flow.data.benchmarks.leakdb.load_scada_data`.
 13The official scoring/evaluation is implemented in
 14:func:`~epyt_flow.data.benchmarks.leakdb.compute_evaluation_score` -- i.e. those results can be
 15directly compared to the official paper.
 16"""
 17import os
 18from typing import Union
 19import math
 20import json
 21import scipy
 22import numpy as np
 23import pandas as pd
 24from scipy.sparse import bsr_array
 25from sklearn.metrics import f1_score, recall_score as true_positive_rate
 26
 27from ..networks import load_net1, load_hanoi
 28from .leakdb_data import NET1_LEAKAGES, HANOI_LEAKAGES
 29from ...utils import get_temp_folder, to_seconds, unpack_zip_archive, create_path_if_not_exist, \
 30    download_if_necessary
 31from ...simulation import ScenarioSimulator, EpanetConstants
 32from ...simulation.events import AbruptLeakage, IncipientLeakage
 33from ...simulation import ScenarioConfig
 34from ...simulation.scada import ScadaData
 35from ...uncertainty import ModelUncertainty, UniformUncertainty
 36
 37
[docs] 38def true_negative_rate(y_pred: np.ndarray, y: np.ndarray) -> float: 39 """ 40 Computes the true negative rate (also called specificity). 41 42 Parameters 43 ---------- 44 y_pred : `numpy.ndarray <https://numpy.org/doc/stable/reference/generated/numpy.ndarray.html>`_ 45 Predicted labels. 46 y : `numpy.ndarray <https://numpy.org/doc/stable/reference/generated/numpy.ndarray.html>`_ 47 Ground truth labels. 48 49 Returns 50 ------- 51 `float` 52 True negative rate. 53 """ 54 if not isinstance(y_pred, np.ndarray): 55 raise TypeError("'y_pred' must be an instance of 'numpy.ndarray' " + 56 f"but not of '{type(y_pred)}'") 57 if not isinstance(y, np.ndarray): 58 raise TypeError("'y' must be an instance of 'numpy.ndarray' " + 59 f"but not of '{type(y)}'") 60 if y_pred.shape != y.shape: 61 raise ValueError(f"Shape mismatch: {y_pred.shape} vs. {y.shape}") 62 if len(y_pred.shape) > 1: 63 raise ValueError("'y_pred' must be a 1d array") 64 if len(y.shape) > 1: 65 raise ValueError("'y' must be a 1d array") 66 if set(np.unique(y_pred)) != set([0, 1]): 67 raise ValueError("Labels must be either '0' or '1'") 68 69 tn = np.sum((y == 0) & (y_pred == 0)) 70 fp = np.sum((y == 0) & (y_pred == 1)) 71 72 return tn / (tn + fp)
73 74 75def __leak_time_to_idx(t: int, round_up: bool = False, hydraulic_time_step: int = 1800): 76 if round_up is False: 77 return math.floor(t / hydraulic_time_step) 78 else: 79 return math.ceil(t / hydraulic_time_step) 80 81 82def __get_leak_time_windows(s_id: int, leaks_info: dict, 83 hydraulic_time_step: int = 1800) -> list[tuple[int, int]]: 84 time_windows = [] 85 if str(s_id) in leaks_info: 86 for leak in leaks_info[str(s_id)]: 87 t_idx_start = __leak_time_to_idx(leak["leak_start_time"] * hydraulic_time_step) 88 t_idx_end = __leak_time_to_idx(leak["leak_end_time"] * hydraulic_time_step, 89 round_up=True) 90 91 time_windows.append((t_idx_start, t_idx_end)) 92 93 return time_windows 94 95 96def __create_labels(s_id: int, n_time_steps: int, nodes: list[str], 97 leaks_info: dict, 98 hydraulic_time_step: int = 1800) -> tuple[np.ndarray, scipy.sparse.bsr_array]: 99 y = np.zeros(n_time_steps) 100 101 leak_locations_row = [] 102 leak_locations_col = [] 103 if str(s_id) in leaks_info: 104 for leak in leaks_info[str(s_id)]: 105 t_idx_start = __leak_time_to_idx(leak["leak_start_time"] * hydraulic_time_step) 106 t_idx_end = __leak_time_to_idx(leak["leak_end_time"] * hydraulic_time_step, 107 round_up=True) 108 109 leak_node_idx = nodes.index(leak["node_id"]) 110 111 for t in range(t_idx_end - t_idx_start): 112 leak_locations_row.append(t_idx_start + t) 113 leak_locations_col.append(leak_node_idx) 114 115 y[t_idx_start:t_idx_end] = 1 116 117 y_leak_locations = bsr_array( 118 (np.ones(len(leak_locations_row)), (leak_locations_row, leak_locations_col)), 119 shape=(n_time_steps, len(nodes))) 120 121 return y, y_leak_locations 122 123
[docs] 124def compute_evaluation_score(scenarios_id: list[int], use_net1: bool, 125 y_pred_labels_per_scenario: list[np.ndarray]) -> dict: 126 """ 127 Evaluates the predictions (leakage detection) for a list of given scenarios. 128 129 Parameters 130 ---------- 131 scenarios_id : `list[int]` 132 List of scenarios ID that are to be evaluated -- there is a total number of 1000 scenarios. 133 use_net1 : `bool` 134 If True, Net1 LeakDB will be used for evaluation, otherwise the Hanoi LeakDB will be used. 135 y_pred_labels_per_scenario : `list[numpy.ndarray] <https://numpy.org/doc/stable/reference/generated/numpy.ndarray.html>`_ 136 Predicted binary labels (over time) for each scenario in `scenarios_id`. 137 138 Returns 139 ------- 140 `dict` 141 Dictionary containing the f1-score, true positive rate, true negative rate, 142 and early detection score. 143 """ 144 # Original MATLAB implementation: https://github.com/KIOS-Research/LeakDB/blob/master/CCWI-WDSA2018/Scoring%20Function/scoring_algorithm.m 145 if len(scenarios_id) != len(y_pred_labels_per_scenario): 146 raise ValueError("Number of scenarios does not match number of predictions -- " + 147 f"expected {len(scenarios_id)} but got {len(y_pred_labels_per_scenario)}") 148 149 # Load ground truth 150 if use_net1 is True: 151 leaks_info = json.loads(NET1_LEAKAGES) 152 else: 153 leaks_info = json.loads(HANOI_LEAKAGES) 154 155 network_config = load_net1() if use_net1 is True \ 156 else load_hanoi() 157 nodes = network_config.sensor_config.nodes 158 159 y_true = [] 160 for i, s_id in enumerate(scenarios_id): 161 y, _ = __create_labels(s_id, len(y_pred_labels_per_scenario[i]), nodes, leaks_info) 162 if len(y) != len(y_pred_labels_per_scenario[i]): 163 raise ValueError("A prediction must be provided for each time step -- " + 164 f"mismatch for scenario {i}, expected {len(y)} but got " + 165 f"{y_pred_labels_per_scenario[i]}") 166 y_true.append(y) 167 168 y_true = np.stack(y_true, axis=0) 169 y_pred = np.stack(y_pred_labels_per_scenario, axis=0) 170 171 # Evaluate predictions 172 f1 = f1_score(y_true, y_pred) 173 tpr = true_positive_rate(y_true, y_pred) 174 tnr = true_negative_rate(y_pred, y_true) 175 176 early_detection_score = 0 177 normalizing = [] 178 n_time_steps_tolerance = 10 179 detection_threshold = .75 180 for i, s_id in enumerate(scenarios_id): 181 y_pred_i = y_pred_labels_per_scenario[i] 182 leaks_time_window = __get_leak_time_windows(s_id, leaks_info) 183 184 scores = [] 185 for t0, _ in leaks_time_window: 186 normalizing.append(1.) 187 188 y_pred_window = y_pred_i[t0:t0+n_time_steps_tolerance] 189 if 1 in y_pred_window and \ 190 np.sum(y_pred_window) / len(y_pred_window) > detection_threshold: 191 t_idx = np.argwhere(y_pred_window)[0] + 1 192 scores.append(2. / (1 + np.exp((5. / n_time_steps_tolerance) * t_idx))) 193 else: 194 scores.append(0.) 195 196 early_detection_score += np.sum(scores) 197 198 early_detection_score = early_detection_score / np.sum(normalizing) 199 200 return {"f1_score": f1, "true_positive_rate": tpr, 201 "true_negative_rate": tnr, "early_detection_score": early_detection_score}
202 203
[docs] 204def load_data(scenarios_id: list[int], use_net1: bool, download_dir: str = None, 205 return_X_y: bool = False, return_features_desc: bool = False, 206 return_leak_locations: bool = False, verbose: bool = True) -> dict: 207 """ 208 Loads the original LeakDB benchmark data set. 209 210 .. warning:: 211 212 All scenarios together are a huge data set -- approx. 8GB for Net1 and 25GB for Hanoi. 213 Downloading and loading might take some time! Also, a sufficient amount of hard disk 214 memory is required. 215 216 Parameters 217 ---------- 218 scenarios_id : `list[int]` 219 List of scenarios ID that are to be loaded -- there are a total number of 1000 scenarios. 220 use_net1 : `bool` 221 If True, Net1 LeakDB will be loaded, otherwise the Hanoi LeakDB will be loaded. 222 download_dir : `str`, optional 223 Path to the data files -- if None, the temp folder will be used. 224 If the path does not exist, the data files will be downloaded to the given path. 225 226 The default is None. 227 return_X_y : `bool`, optional 228 If True, the data is returned together with the labels (presence of a leakage) as 229 two Numpy arrays, otherwise, the data is returned as Pandas data frames. 230 231 The default is False. 232 return_features_desc : `bool`, optional 233 If True and if `return_X_y` is True, the returned dictionary contains the 234 features' descriptions (i.e. names) under the key "features_desc". 235 236 The default is False. 237 return_leak_locations : `bool` 238 If True and if `return_X_y` is True, the leak locations are returned as well -- 239 as an instance of `scipy.sparse.bsr_array <https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.bsr_array.html>`_. 240 241 The default is False. 242 verbose : `bool`, optional 243 If True, a progress bar is shown while downloading files. 244 245 The default is True. 246 247 Returns 248 ------- 249 `dict` 250 Dictionary containing the scenario data sets. Data of each requested scenario 251 can be accessed by using the scenario ID as a key. 252 """ 253 url_data = "https://filedn.com/lumBFq2P9S74PNoLPWtzxG4/EPyT-Flow/LeakDB-Original/" +\ 254 f"{'Net1_CMH/' if use_net1 is True else 'Hanoi_CMH/'}" 255 256 if use_net1 is True: 257 network_desc = "Net1" 258 leaks_info = json.loads(NET1_LEAKAGES) 259 else: 260 network_desc = "Hanoi" 261 leaks_info = json.loads(HANOI_LEAKAGES) 262 263 download_dir = download_dir if download_dir is not None else get_temp_folder() 264 download_dir = os.path.join(download_dir, network_desc) 265 create_path_if_not_exist(download_dir) 266 267 results = {} 268 for s_id in scenarios_id: 269 scenario_data = f"Scenario-{s_id}.zip" 270 scenario_data_url = url_data + scenario_data 271 scenario_data_file_in = os.path.join(download_dir, scenario_data) 272 scenario_data_folder_in = os.path.join(download_dir, f"Scenario-{s_id}") 273 274 download_if_necessary(scenario_data_file_in, scenario_data_url, verbose) 275 create_path_if_not_exist(scenario_data_folder_in) 276 unpack_zip_archive(scenario_data_file_in, scenario_data_folder_in) 277 278 scenario_data_folder_in = os.path.join(scenario_data_folder_in, f"Scenario-{s_id}") 279 280 # Load and parse data 281 pressure_files = list(filter(lambda d: d.endswith(".csv"), 282 os.listdir(os.path.join(scenario_data_folder_in, 283 "Pressures")))) 284 pressure_readings = {} 285 all_nodes = [] 286 for f_in in pressure_files: 287 df = pd.read_csv(os.path.join(scenario_data_folder_in, "Pressures", f_in)) 288 node_id = f_in.replace(".csv", "") 289 all_nodes.append(node_id) 290 pressure_readings[f"Pressure-{node_id}"] = df["Value"] 291 292 flow_files = list(filter(lambda d: d.endswith(".csv"), 293 os.listdir(os.path.join(scenario_data_folder_in, "Flows")))) 294 flow_readings = {} 295 for f_in in flow_files: 296 df = pd.read_csv(os.path.join(scenario_data_folder_in, "Flows", f_in)) 297 flow_readings[f"Flow-{f_in.replace('.csv', '')}"] = df["Value"] 298 299 df_labels = pd.read_csv(os.path.join(scenario_data_folder_in, "Labels.csv")) 300 labels = df_labels["Label"] 301 302 df_timestamps = pd.read_csv(os.path.join(scenario_data_folder_in, "Timestamps.csv")) 303 sensor_reading_times = df_timestamps["Timestamp"] 304 305 df_final = pd.DataFrame(pressure_readings | flow_readings | 306 {"labels": labels, "timestamps": sensor_reading_times}) 307 308 # Prepare final data 309 if return_X_y is True: 310 X = df_final[list(pressure_readings.keys()) + list(flow_readings.keys())].to_numpy() 311 y = labels.to_numpy() 312 313 network_config = load_net1(download_dir) if use_net1 is True \ 314 else load_hanoi(download_dir) 315 nodes = network_config.sensor_config.nodes 316 _, y_leak_locations = __create_labels(s_id, X.shape[0], nodes, leaks_info) 317 318 if return_features_desc is True and "features_desc" not in results: 319 results["features_desc"] = list(pressure_readings.keys()) + \ 320 list(flow_readings.keys()) 321 322 if return_leak_locations is True: 323 results[s_id] = (X, y, y_leak_locations) 324 else: 325 results[s_id] = (X, y) 326 else: 327 results[s_id] = df_final 328 329 return results
330 331
[docs] 332def load_scada_data(scenarios_id: list[int], use_net1: bool = True, download_dir: str = None, 333 return_X_y: bool = False, return_leak_locations: bool = False, 334 verbose: bool = True 335 ) -> Union[list[ScadaData], list[tuple[np.ndarray, np.ndarray]]]: 336 """ 337 Loads the SCADA data of the simulated LeakDB benchmark scenarios -- see 338 :func:`~epyt_flow.data.benchmarks.leakdb.load_scenarios`. 339 340 .. note:: 341 Note that due to the randomness in the demand creation as well as in the model 342 uncertainties, the SCADA data differs from the original data set 343 which can be loaded by calling :func:`~epyt_flow.data.benchmarks.leakdb.load_data`. 344 However, the leakages (i.e. location and profile) are consistent with the original data set. 345 346 Parameters 347 ---------- 348 scenarios_id : `list[int]` 349 List of scenarios ID that are to be loaded -- there are a total number of 1000 scenarios. 350 use_net1 : `bool`, optional 351 If True, Net1 LeakDB will be loaded, otherwise the Hanoi LeakDB will be loaded. 352 353 The default is True. 354 download_dir : `str`, optional 355 Path to the data files -- if None, the temp folder will be used. 356 If the path does not exist, the data files will be downloaded to the given path. 357 358 The default is None. 359 return_X_y : `bool`, optional 360 If True, the data is returned together with the labels (presence of a leakage) as 361 two Numpy arrays, otherwise, the data is returned as 362 :class:`~epyt_flow.simulation.scada.scada_data.ScadaData` instances. 363 364 The default is False. 365 return_leak_locations : `bool` 366 If True, the leak locations are returned as well -- 367 as an instance of `scipy.sparse.bsr_array <https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.bsr_array.html>`_. 368 369 The default is False. 370 verbose : `bool`, optional 371 If True, a progress bar is shown while downloading files. 372 373 The default is True. 374 375 Returns 376 ------- 377 list[`:class:`~epyt_flow.simulation.scada.scada_data.ScadaData`] or `list[tuple[numpy.ndarray, numpy.ndarray]]` 378 The simulated benchmark scenarios as either a list of 379 :class:`~epyt_flow.simulation.scada.scada_data.ScadaData` instances or as a list of 380 (X, y) Numpy arrays. If 'return_leak_locations' is True, the leak locations are included 381 as an instance of `scipy.sparse.bsr_array <https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.bsr_array.html>`_ as well. 382 """ 383 download_dir = download_dir if download_dir is not None else get_temp_folder() 384 385 url_data = "https://filedn.com/lumBFq2P9S74PNoLPWtzxG4/EPyT-Flow/LeakDB/" +\ 386 f"{'Net1/' if use_net1 is True else 'Hanoi/'}" 387 388 if use_net1 is True: 389 leaks_info = json.loads(NET1_LEAKAGES) 390 else: 391 leaks_info = json.loads(HANOI_LEAKAGES) 392 393 r = [] 394 395 for s_id in scenarios_id: 396 f_in = f"{'Net1_ID' if use_net1 is True else 'Hanoi_ID'}={s_id}.epytflow_scada_data" 397 download_if_necessary(os.path.join(download_dir, f_in), url_data + f_in, verbose) 398 399 data = ScadaData.load_from_file(os.path.join(download_dir, f_in)) 400 401 X = data.get_data() 402 y, y_leak_locations = __create_labels(s_id, X.shape[0], data.sensor_config.nodes, 403 leaks_info) 404 405 if return_X_y is True: 406 if return_leak_locations is True: 407 r.append((X, y, y_leak_locations)) 408 else: 409 r.append((X, y)) 410 else: 411 if return_leak_locations is True: 412 r.append((data, y_leak_locations)) 413 else: 414 r.append(data) 415 416 return r
417 418
[docs] 419def load_scenarios(scenarios_id: list[int], use_net1: bool = True, 420 download_dir: str = None, verbose: bool = True) -> list[ScenarioConfig]: 421 """ 422 Creates and returns the LeakDB scenarios -- they can be either modified or 423 passed directly to the simulator 424 :class:`~epyt_flow.simulation.scenario_simulator.ScenarioSimulator`. 425 426 .. note:: 427 Note that due to the randomness in the demand creation as well as in the model 428 uncertainties, the simulation results will differ between different runs, and 429 will also differ from the original data set 430 (see :func:`~epyt_flow.data.benchmarks.leakdb.load_data`). 431 However, the leakages (i.e. location and profile) will be always the same and be 432 consistent with the original data set. 433 434 Parameters 435 ---------- 436 scenarios_id : `list[int]` 437 List of scenarios ID that are to be loaded -- there is a total number of 1000 scenarios. 438 use_net1 : `bool`, optional 439 If True, Net1 network will be used, otherwise the Hanoi network will be used. 440 441 The default is True. 442 download_dir : `str`, optional 443 Path to the Net1.inp or Hanoi.inp file -- if None, the temp folder will be used. 444 If the path does not exist, the .inp will be downloaded to the give path. 445 446 The default is None. 447 verbose : `bool`, optional 448 If True, a progress bar is shown while downloading files. 449 450 The default is True. 451 452 Returns 453 ------- 454 list[:class:`~epyt_flow.simulation.scenario_config.ScenarioConfig`] 455 LeakDB scenarios. 456 """ 457 scenarios_inp = [] 458 459 # Load the network 460 load_network = load_net1 if use_net1 is True else load_hanoi 461 download_dir = download_dir if download_dir is not None else get_temp_folder() 462 network_config = load_network(download_dir) 463 464 # Set simulation duration and other general parameters such as the demand model and flow units 465 hydraulic_time_step = to_seconds(minutes=30) # 30min time steps 466 general_params = {"simulation_duration": to_seconds(days=365), # One year 467 "hydraulic_time_step": hydraulic_time_step, 468 "reporting_time_step": hydraulic_time_step, 469 "flow_units_id": EpanetConstants.EN_CMH, 470 "demand_model": {"type": EpanetConstants.EN_PDA, "pressure_min": 0, 471 "pressure_required": 0.1, 472 "pressure_exponent": 0.5} 473 } | network_config.general_params 474 475 # Add demand patterns 476 def gen_dem(download_dir): 477 # Taken from https://github.com/KIOS-Research/LeakDB/blob/master/CCWI-WDSA2018/Dataset_Generator_Py3/demandGenerator.py 478 week_pat = scipy.io.loadmat(os.path.join(download_dir, "weekPat_30min.mat")) 479 a_w = week_pat['Aw'] 480 nw = week_pat['nw'] 481 year_offset = scipy.io.loadmat(os.path.join(download_dir, "yearOffset_30min.mat")) 482 a_y = year_offset['Ay'] 483 ny = year_offset['ny'] 484 485 # Create yearly component 486 days = 365 487 488 t = (288/6)*days # one year period in five minute intervals 489 w = 2*np.pi/t 490 k = np.arange(1, days*288/6+1, 1) # number of time steps in time series 491 n = ny[0][0] # number of fourier coefficients 492 h_y = [1]*len(k) 493 494 for i in range(1, n+1): 495 h_y = np.column_stack((h_y, np.sin(i*w*k), np.cos(i*w*k))) 496 497 unc_y = 0.1 498 a_y_r = a_y*(1-unc_y + 2*unc_y*np.random.rand(int(a_y.shape[0]), int(a_y.shape[1]))) 499 year_offset = np.dot(h_y, a_y_r) 500 501 # Create weekly component 502 t = (288/6)*7 # one week period in five minute intervals 503 w = 2*np.pi/t 504 k = np.arange(1, days*288/6+1, 1) # number of time steps in time series 505 n = nw[0][0] # number of fourier coefficients 506 h_w = [1]*len(k) 507 for i in range(1, n+1): 508 h_w = np.column_stack((h_w, np.sin(i*w*k), np.cos(i*w*k))) 509 510 unc_w = 0.1 511 a_w_r = a_w*(1-unc_w + 2*unc_w*np.random.rand(int(a_w.shape[0]), int(a_w.shape[1]))) 512 week_year_pat = np.dot(h_w, a_w_r) 513 514 # Create random component 515 unc_r = 0.05 516 random = np.random.normal(0, (-unc_r+2*unc_r), 517 (int(week_year_pat.shape[0]), int(week_year_pat.shape[1]))) 518 519 # Create demand 520 base = 1 521 variation = 0.75 + np.random.normal(0, 0.07) # from 0 to 1 522 dem = base * (year_offset+1) * (week_year_pat*variation+1) * (random+1) 523 524 dem = dem.tolist() 525 dem_final = [] 526 for d in dem: 527 dem_final.append(d[0]) 528 529 return dem_final 530 531 week_pattern_url = "https://filedn.com/lumBFq2P9S74PNoLPWtzxG4/EPyT-Flow/Networks/CCWI-WDSA2018/" +\ 532 "Dataset_Generator_Py3/weekPat_30min.mat" 533 year_offset_url = "https://filedn.com/lumBFq2P9S74PNoLPWtzxG4/EPyT-Flow/Networks/CCWI-WDSA2018/" +\ 534 "Dataset_Generator_Py3/yearOffset_30min.mat" 535 536 download_if_necessary(os.path.join(download_dir, "weekPat_30min.mat"), 537 week_pattern_url, verbose) 538 download_if_necessary(os.path.join(download_dir, "yearOffset_30min.mat"), 539 year_offset_url, verbose) 540 541 for s_id in scenarios_id: # Create new .inp files with demands if necessary 542 f_inp_in = os.path.join(download_dir, 543 f"{'Net1' if use_net1 is True else 'Hanoi'}_LeakDB_ID={s_id}.inp") 544 scenarios_inp.append(f_inp_in) 545 546 if not os.path.exists(f_inp_in): 547 with ScenarioSimulator(f_inp_in=network_config.f_inp_in) as wdn: 548 wdn.set_general_parameters(**general_params) 549 wdn.epanet_api.set_hydraulic_time_step(hydraulic_time_step) 550 wdn.epanet_api.settimeparam(EpanetConstants.EN_PATTERNSTEP, hydraulic_time_step) 551 552 for idx in range(1, wdn.epanet_api.getcount(EpanetConstants.EN_PATCOUNT) + 1): 553 wdn.epanet_api.deletepattern(idx) 554 555 reservoir_nodes_id = wdn.epanet_api.get_all_reservoirs_id() 556 for node_id in network_config.sensor_config.nodes: 557 if node_id in network_config.sensor_config.tanks or\ 558 node_id in reservoir_nodes_id: 559 continue 560 561 node_idx = wdn.epanet_api.get_node_idx(node_id) 562 base_demand = wdn.epanet_api.get_node_base_demand(node_idx) 563 564 my_demand_pattern = np.array(gen_dem(download_dir)) 565 566 wdn.set_node_demand_pattern(node_id=node_id, base_demand=base_demand, 567 demand_pattern_id=f"demand_{node_id}", 568 demand_pattern=my_demand_pattern) 569 570 wdn.epanet_api.saveinpfile(f_inp_in) 571 572 # Create uncertainties 573 class MyUniformUncertainty(UniformUncertainty): 574 """ 575 Custom uniform uncertainty for LeakDB scenarios. 576 """ 577 def __init__(self, **kwds): 578 super().__init__(**kwds) 579 580 def apply(self, data: float) -> float: 581 z = data * np.random.uniform(low=self.low, high=self.high) 582 lower = data - z 583 upper = data + z 584 return lower + np.random.uniform() * (upper - lower) 585 586 my_uncertainties = {"global_pipe_length_uncertainty": MyUniformUncertainty(low=0, high=0.25), 587 "global_pipe_roughness_uncertainty": MyUniformUncertainty(low=0, high=0.25), 588 "global_base_demand_uncertainty": MyUniformUncertainty(low=0, high=0.25)} 589 model_uncertainty = ModelUncertainty(**my_uncertainties) 590 591 # Create sensor config (place pressure and flow sensors everywhere) 592 sensor_config = network_config.sensor_config 593 sensor_config.pressure_sensors = sensor_config.nodes 594 sensor_config.flow_sensors = sensor_config.links 595 596 # Add leakages 597 leaks_all = [] 598 599 if use_net1 is True: 600 leaks_info = json.loads(NET1_LEAKAGES) 601 else: 602 leaks_info = json.loads(HANOI_LEAKAGES) 603 604 for s_id in scenarios_id: 605 leaks_data = [] 606 607 if str(s_id) in leaks_info: 608 for leak in leaks_info[str(s_id)]: 609 if leak["leak_type"] == "incipient": 610 leaks_data.append( 611 IncipientLeakage(node_id=leak["node_id"], link_id=None, 612 diameter=leak["leak_diameter"], 613 start_time=leak["leak_start_time"] * hydraulic_time_step, 614 end_time=leak["leak_end_time"] * hydraulic_time_step, 615 peak_time=leak["leak_peak_time"] * hydraulic_time_step)) 616 else: 617 leaks_data.append( 618 AbruptLeakage(node_id=leak["node_id"], link_id=None, 619 diameter=leak["leak_diameter"], 620 start_time=leak["leak_start_time"] * hydraulic_time_step, 621 end_time=leak["leak_end_time"] * hydraulic_time_step)) 622 623 leaks_all.append(leaks_data) 624 625 # Build final scenarios 626 return [ScenarioConfig(f_inp_in=f_inp_in, general_params=general_params, 627 sensor_config=sensor_config, model_uncertainty=model_uncertainty, 628 system_events=leaks) 629 for f_inp_in, leaks in zip(scenarios_inp, leaks_all)]