Merge pull request #100 from ImperialCollegeLondon/91-design-output-m…

…etrics-as-in-metrics-for-the-networks-design-parameters

Metrics for the designed network

swmmanywhere/metric_utilities.py

@@ -517,6 +517,121 @@ def outlet_nse_flooding(synthetic_G: nx.Graph,
                          list(sg_syn.nodes),
                          list(sg_real.nodes))
 
+@metrics.register
+def outlet_kstest_diameters(real_G: nx.Graph,
+                            synthetic_G: nx.Graph,
+                            real_results: pd.DataFrame,
+                            real_subs: gpd.GeoDataFrame,
+                            **kwargs) -> float:
+    """Outlet KStest diameters.
+
+    Calculate the Kolmogorov-Smirnov statistic of the diameters in the subgraph
+    that drains to the dominant outlet node. The dominant outlet node of the
+    'real' network is calculated by dominant_outlet, while the dominant outlet
+    node of the 'synthetic' network is calculated by best_outlet_match.
+    """
+    # Identify synthetic and real outlet arcs
+    sg_syn, _ = best_outlet_match(synthetic_G, real_subs)
+    sg_real, _ = dominant_outlet(real_G, real_results)
+
+    # Extract the diameters
+    syn_diameters = nx.get_edge_attributes(sg_syn, 'diameter')
+    real_diameters = nx.get_edge_attributes(sg_real, 'diameter')
+    return stats.ks_2samp(list(syn_diameters.values()),
+                         list(real_diameters.values())).statistic
+
+@metrics.register
+def outlet_pbias_length(real_G: nx.Graph,
+                        synthetic_G: nx.Graph,
+                        real_results: pd.DataFrame,
+                        real_subs: gpd.GeoDataFrame,
+                        **kwargs) -> float:
+    r"""Outlet PBIAS length.
+
+    Calculate the percent bias of the total edge length in the subgraph that
+    drains to the dominant outlet node. The dominant outlet node of the 'real'
+    network is calculated by dominant_outlet, while the dominant outlet node of
+    the 'synthetic' network is calculated by best_outlet_match.
+
+    The percentage bias is calculated as:
+
+    .. math::
+
+        pbias = \\frac{{syn\_length - real\_length}}{{real\_length}}
+
+    where:
+    - :math:`syn\_length` is the synthetic length,
+    - :math:`real\_length` is the real length.
+    """
+    # Identify synthetic and real outlet arcs
+    sg_syn, _ = best_outlet_match(synthetic_G, real_subs)
+    sg_real, _ = dominant_outlet(real_G, real_results)
+
+    # Calculate the percent bias
+    syn_length = sum([d['length'] for u,v,d in sg_syn.edges(data=True)])
+    real_length = sum([d['length'] for u,v,d in sg_real.edges(data=True)])
+    return (syn_length - real_length) / real_length
+
+@metrics.register
+def outlet_pbias_nmanholes(real_G: nx.Graph,
+                           synthetic_G: nx.Graph,
+                           real_results: pd.DataFrame,
+                           real_subs: gpd.GeoDataFrame,
+                           **kwargs) -> float:
+    r"""Outlet PBIAS number of manholes (nodes).
+
+    Calculate the percent bias of the total number of nodes in the subgraph
+    that drains to the dominant outlet node. The dominant outlet node of the
+    'real' network is calculated by dominant_outlet, while the dominant outlet
+    node of the 'synthetic' network is calculated by best_outlet_match.
+
+    The percentage bias is calculated as:
+
+    .. math::
+
+        pbias = \\frac{{syn\_nnodes - real\_nnodes}}{{real\_nnodes}}
+
+    where:
+    - :math:`syn\_nnodes` is the number of synthetic nodes,
+    - :math:`real\_nnodes` is the real number of nodes.
+    """
+    # Identify synthetic and real outlet arcs
+    sg_syn, _ = best_outlet_match(synthetic_G, real_subs)
+    sg_real, _ = dominant_outlet(real_G, real_results)
+
+    return (sg_syn.number_of_nodes() - sg_real.number_of_nodes()) \
+        / sg_real.number_of_nodes()
+
+@metrics.register
+def outlet_pbias_npipes(real_G: nx.Graph,
+                        synthetic_G: nx.Graph,
+                        real_results: pd.DataFrame,
+                        real_subs: gpd.GeoDataFrame,
+                        **kwargs) -> float:
+    r"""Outlet PBIAS number of pipes (edges).
+
+    Calculate the percent bias of the total number of edges in the subgraph
+    that drains to the dominant outlet node. The dominant outlet node of the
+    'real' network is calculated by dominant_outlet, while the dominant outlet
+    node of the 'synthetic' network is calculated by best_outlet_match.
+
+    
+    The percentage bias is calculated as:
+
+    .. math::
+
+        pbias = \\frac{{syn\_nedges - real\_nedges}}{{real\_nedges}}
+
+    where:
+    - :math:`syn\_nedges` is the number of synthetic edges,
+    - :math:`real\_nedges` is the real number of edges.
+    """
+    # Identify synthetic and real outlet arcs
+    sg_syn, _ = best_outlet_match(synthetic_G, real_subs)
+    sg_real, _ = dominant_outlet(real_G, real_results)
+
+    return (sg_syn.number_of_edges() - sg_real.number_of_edges()) \
+        / sg_real.number_of_edges()
 
 
 @metrics.register
@@ -539,4 +654,4 @@ def subcatchment_nse_flooding(synthetic_G: nx.Graph,
                                     synthetic_G = synthetic_G,
                                     real_G = real_G)
 
-    return median_nse_by_group(results, 'sub_id')
+    return median_nse_by_group(results, 'sub_id')

tests/test_graph_utilities.py

@@ -4,6 +4,7 @@
 @author: Barney
 """
 import math
+import os
 import tempfile
 from pathlib import Path
 
@@ -15,6 +16,7 @@
 from swmmanywhere.graph_utilities import graphfcns as gu
 from swmmanywhere.graph_utilities import iterate_graphfcns, load_graph, save_graph
 
+os.environ['SWMMANYWHERE_VERBOSE'] = "false"
 
 def load_street_network():
     """Load a street network."""

tests/test_metric_utilities.py

@@ -279,6 +279,58 @@ def test_outlet_nse_flooding():
                                     real_subs = subs)
     assert val == 0.0
 
+def test_design_params():
+    """Test the design param related metrics."""
+    G = load_graph(Path(__file__).parent / 'test_data' / 'graph_topo_derived.json')
+    nx.set_edge_attributes(G, 0.15, 'diameter')
+    subs = get_subs()
+
+    # Mock results (only needed for dominant outlet)
+    results = pd.DataFrame([{'id' : 4253560,
+                             'variable' : 'flow',
+                             'value' : 10,
+                             'date' : pd.to_datetime('2021-01-01 00:00:00')},
+                             {'id' : 4253560,
+                             'variable' : 'flow',
+                             'value' : 5,
+                             'date' : pd.to_datetime('2021-01-01 00:00:05')},
+                             ])
+
+    # Target results
+    design_results = {'outlet_kstest_diameters' : 0.0625,
+               'outlet_pbias_length' : -0.15088965,
+               'outlet_pbias_nmanholes' : -0.05,
+               'outlet_pbias_npipes' : -0.15789473}
+
+    # Iterate for G = G, i.e., perfect results
+    metrics = mu.iterate_metrics(synthetic_G = G,
+                                 synthetic_subs = None,
+                                 synthetic_results = None,
+                                 real_G = G,
+                                 real_subs = subs,
+                                 real_results = results,
+                                 metric_list = design_results.keys())
+    for metric, val in metrics.items():
+        assert metric in design_results
+        assert np.isclose(val, 0)
+
+    # edit the graph for target results
+    G_ = G.copy()
+    G_.remove_node(list(G.nodes)[0])
+    G_.edges[list(G_.edges)[0]]['diameter'] = 0.3
+
+    metrics = mu.iterate_metrics(synthetic_G = G_,
+                                 synthetic_subs = None,
+                                 synthetic_results = None,
+                                 real_G = G,
+                                 real_subs = subs,
+                                 real_results = results,
+                                 metric_list = design_results.keys())
+
+    for metric, val in metrics.items():
+        assert metric in design_results
+        assert np.isclose(val, design_results[metric]), metric
+
 def test_netcomp_iterate():
     """Test the netcomp metrics and iterate_metrics."""
     netcomp_results = {'nc_deltacon0' : 0.00129408,