fnm: evaluate rupture jumps on the actual connections

openquake/fnm/rupture_connections.py

@@ -1653,9 +1653,103 @@ def rdist_to_dist_matrix(
     return dist_matrix
 
 
+def adj_matrix_to_nbr_dict(adj_matrix) -> Dict[int, Dict[int, float]]:
+    """
+    Converts a (dense or sparse) adjacency matrix into a nested dictionary
+    `{i: {j: value}}` holding only the non-zero entries.
+    """
+    nbr_dict = {}
+    if issparse(adj_matrix):
+        mat = adj_matrix.tocsr()
+        indptr, indices, data = mat.indptr, mat.indices, mat.data
+        for i in range(mat.shape[0]):
+            row = {
+                int(j): float(v)
+                for j, v in zip(
+                    indices[indptr[i] : indptr[i + 1]],
+                    data[indptr[i] : indptr[i + 1]],
+                )
+                if v != 0.0
+            }
+            if row:
+                nbr_dict[i] = row
+    else:
+        arr = np.asarray(adj_matrix)
+        for i, j in zip(*np.nonzero(arr)):
+            nbr_dict.setdefault(int(i), {})[int(j)] = float(arr[i, j])
+    return nbr_dict
+
+
+def get_mf_connection_pairs(mf, nbr_dict) -> list[tuple[int, int]]:
+    """
+    Get the connections (pairs of single-fault ruptures) that join the
+    ruptures of a multifault rupture. The rupture indices in `mf` are not
+    in any particular order, so the connections are taken from the
+    adjacency structure: they are the edges of the minimum spanning tree
+    (by adjacency value, i.e. distance) of the subgraph induced by `mf`.
+
+    Parameters
+    ----------
+    mf : list of int
+        List of rupture indices.
+    nbr_dict : dict
+        Nested dictionary `{i: {j: value}}`, see `adj_matrix_to_nbr_dict`.
+
+    Returns
+    -------
+    list of tuples
+        The `len(mf) - 1` connections for a connected rupture.
+    """
+    if len(mf) < 2:
+        return []
+    members = set(int(r) for r in mf)
+    in_tree = [int(mf[0])]
+    in_tree_set = set(in_tree)
+    pairs = []
+    while len(in_tree) < len(members):
+        best = None
+        for u in in_tree:
+            for w, val in nbr_dict.get(u, {}).items():
+                if w in members and w not in in_tree_set:
+                    if best is None or val < best[2]:
+                        best = (u, w, val)
+        if best is None:
+            logging.warning("Multifault rupture %s is not connected", mf)
+            break
+        pairs.append((best[0], best[1]))
+        in_tree.append(best[1])
+        in_tree_set.add(best[1])
+    return pairs
+
+
+def get_multifault_connection_pairs(rupture_df, adj_matrix) -> pd.Series:
+    """
+    Get the connections for every rupture in `rupture_df` (see
+    `get_mf_connection_pairs`) as a Pandas Series indexed like `rupture_df`.
+    """
+    nbr_dict = adj_matrix_to_nbr_dict(adj_matrix)
+    return pd.Series(
+        data=[
+            get_mf_connection_pairs(row.ruptures, nbr_dict)
+            for row in rupture_df.itertuples()
+        ],
+        index=rupture_df.index,
+        dtype=object,
+    )
+
+
+def _distances_from_pairs(pairs, nbr_dict) -> np.ndarray:
+    if len(pairs) == 0:
+        return np.zeros(1)
+    return np.array(
+        [nbr_dict.get(i, {}).get(j, 0.0) for i, j in pairs], dtype=float
+    )
+
+
 def get_mf_distances_from_adj_matrix(mf, dist_adj_matrix) -> np.ndarray:
     """
-    Get the distances between the ruptures in a multifault rupture.
+    Get the distances between the ruptures in a multifault rupture, i.e.
+    the distances of the connections returned by `get_mf_connection_pairs`.
 
     Parameters
     ----------
@@ -1671,13 +1765,19 @@ def get_mf_distances_from_adj_matrix(mf, dist_adj_matrix) -> np.ndarray:
     """
     if len(mf) == 1:
         return np.zeros(1)
-    distances = np.empty(len(mf) - 1)
-    for i in range(len(distances)):
-        distances[i] = dist_adj_matrix[mf[i], mf[i + 1]]
-    return distances
+    idx = np.asarray(mf, dtype=int)
+    if issparse(dist_adj_matrix):
+        sub = dist_adj_matrix.tocsr()[idx][:, idx].toarray()
+    else:
+        sub = np.asarray(dist_adj_matrix)[np.ix_(idx, idx)]
+    local = adj_matrix_to_nbr_dict(sub)
+    local_pairs = get_mf_connection_pairs(list(range(len(idx))), local)
+    return _distances_from_pairs(local_pairs, local)
 
 
-def get_multifault_rupture_distances(rupture_df, distance_matrix) -> pd.Series:
+def get_multifault_rupture_distances(
+    rupture_df, distance_matrix, pairs: Optional[pd.Series] = None
+) -> pd.Series:
     """
     Get the distances between the ruptures in a multifault rupture.
     The distances are given as a Pandas Series, where the index is the rupture
@@ -1693,17 +1793,28 @@ def get_multifault_rupture_distances(rupture_df, distance_matrix) -> pd.Series:
             fault: fault identifier
     distance_matrix : np.ndarray
         Rupture adjacency matrix.
+    pairs : pd.Series, optional
+        Connections for each rupture, as returned by
+        `get_multifault_connection_pairs`. Computed from `distance_matrix`
+        if not given.
 
     Returns
     -------
     pd.Series
         Series of distances between the ruptures in the multifault rupture.
     """
+    nbr_dict = adj_matrix_to_nbr_dict(distance_matrix)
+    if pairs is None:
+        pairs = pd.Series(
+            data=[
+                get_mf_connection_pairs(row.ruptures, nbr_dict)
+                for row in rupture_df.itertuples()
+            ],
+            index=rupture_df.index,
+            dtype=object,
+        )
     distances = pd.Series(
-        data=[
-            get_mf_distances_from_adj_matrix(row.ruptures, distance_matrix)
-            for row in rupture_df.itertuples()
-        ],
+        data=[_distances_from_pairs(p, nbr_dict) for p in pairs.values],
         index=rupture_df.index,
     )
 

openquake/fnm/rupture_filtering.py

@@ -49,6 +49,7 @@
 
 from openquake.fnm.rupture_connections import (
     get_multifault_rupture_distances,
+    get_multifault_connection_pairs,
     get_proximal_rup_angles,
 )
 
@@ -206,23 +207,29 @@ def _padded_array_from_sequences(
     return arr
 
 
-def _matrix_to_rupture_series(rupture_df, matrix):
+def _matrix_to_rupture_series(rupture_df, matrix, pairs=None):
+    """
+    Reads the values of `matrix` at the connections of each rupture.
+    `pairs` holds the connections of each rupture (see
+    `get_multifault_connection_pairs`); if not given they are derived from
+    the non-zero structure of `matrix` itself.
+    """
     if matrix is None:
         return None
 
+    if pairs is None:
+        pairs = get_multifault_connection_pairs(rupture_df, matrix)
+
     mat = matrix.todok() if issparse(matrix) else np.asarray(matrix)
     is_sparse = issparse(mat)
     seqs = []
-    for row in rupture_df.itertuples():
-        rup_indices = getattr(row, "ruptures")
-        if len(rup_indices) == 1:
+    for rup_pairs in pairs.values:
+        if len(rup_pairs) == 0:
             seqs.append(np.zeros(1, dtype=np.float64))
             continue
-
-        vals = []
-        for idx in range(len(rup_indices) - 1):
-            i, j = rup_indices[idx], rup_indices[idx + 1]
-            vals.append(float(mat[i, j] if is_sparse else mat[i][j]))
+        vals = [
+            float(mat[i, j] if is_sparse else mat[i][j]) for i, j in rup_pairs
+        ]
         seqs.append(np.asarray(vals, dtype=np.float64))
 
     return pd.Series(seqs, index=rupture_df.index)
@@ -278,6 +285,7 @@ def _compute_angle_plausibility(
     connection_angle_function,
     connection_angle_midpoint,
     no_connection_val=-1.0,
+    pairs=None,
 ):
     n_rups = rupture_df.shape[0]
     if connection_angle_midpoint is None:
@@ -306,7 +314,7 @@ def _compute_angle_plausibility(
                 np.asarray(angle_matrix_input), dtype=np.float64
             )
         angle_series = _matrix_to_rupture_series(
-            rupture_df, angle_matrix_input
+            rupture_df, angle_matrix_input, pairs=pairs
         )
 
     if angle_series is None:
@@ -440,13 +448,22 @@ def get_rupture_plausibilities(
     if rupture_df.shape[0] == 0:
         return pd.DataFrame(index=rupture_df.index, columns=columns)
 
+    # The connections of each multifault rupture, shared by the distance
+    # and angle plausibilities so both are evaluated on the same jumps.
+    pair_matrix = distance_matrix if distance_matrix is not None else bin_adj_mat
+    pairs = (
+        get_multifault_connection_pairs(rupture_df, pair_matrix)
+        if pair_matrix is not None
+        else None
+    )
+
     if distances is None:
         if distance_matrix is None:
             raise ValueError(
                 "Either distances or distance_matrix must be provided."
             )
         distances = get_multifault_rupture_distances(
-            rupture_df, distance_matrix
+            rupture_df, distance_matrix, pairs=pairs
         )
 
     conn_total = _compute_distance_plausibility(
@@ -466,6 +483,7 @@ def get_rupture_plausibilities(
         subfaults,
         connection_angle_function,
         connection_angle_midpoint,
+        pairs=pairs,
     )
 
     slip_plaus = _compute_slip_az_plausibility(

openquake/fnm/tests/test_rupture_connections.py

@@ -93,6 +93,8 @@
     get_multifault_rupture_distances,
     subgraphs_from_connected_components,
     get_rupture_grouping,
+    get_mf_connection_pairs,
+    adj_matrix_to_nbr_dict,
 )
 
 HERE = pathlib.Path(__file__).parent.absolute()
@@ -1959,3 +1961,28 @@ def test__lonlat_to_xyz_roundtrip():
 
 def test__geog_vec_to_xyz():
     pass
+
+
+class TestConnectionPairs(unittest.TestCase):
+    def test_pairs_follow_connections_not_list_order(self):
+        # path 2 - 9 - 5; the rupture list order puts 2 and 5 side by side
+        dist = dok_array((10, 10), dtype=np.float32)
+        for a, b in [(2, 9), (9, 5)]:
+            dist[a, b] = dist[b, a] = 8.0
+        nbr = adj_matrix_to_nbr_dict(dist)
+        pairs = get_mf_connection_pairs([2, 5, 9], nbr)
+        self.assertEqual(
+            {frozenset(p) for p in pairs}, {frozenset({2, 9}), frozenset({9, 5})}
+        )
+        np.testing.assert_array_almost_equal(
+            get_mf_distances_from_adj_matrix([2, 5, 9], dist), [8.0, 8.0]
+        )
+
+    def test_pairs_use_shortest_jumps_in_cycle(self):
+        dist = np.array(
+            [[0.0, 1.0, 9.0], [1.0, 0.0, 2.0], [9.0, 2.0, 0.0]]
+        )
+        np.testing.assert_array_almost_equal(
+            sorted(get_mf_distances_from_adj_matrix([0, 1, 2], dist)),
+            [1.0, 2.0],
+        )

openquake/fnm/tests/test_rupture_filtering.py

@@ -104,6 +104,29 @@ def test_connection_angle_plausibility_from_sparse_angles(self):
         )
 
 
+class TestPlausibilityRegressions(unittest.TestCase):
+    def test_distance_plausibility_counts_every_jump(self):
+        # path 2 - 9 - 5 with two 8 km jumps; the rupture list order puts
+        # the unconnected ruptures 2 and 5 next to each other
+        dist = dok_array((10, 10), dtype=np.float32)
+        for a, b in [(2, 9), (9, 5)]:
+            dist[a, b] = dist[b, a] = 8.0
+        rup_df = pd.DataFrame(
+            {"ruptures": [[2, 5, 9], [2, 9]], "slip_azimuth": [[0.0] * 3, [0.0] * 2]}
+        )
+        plaus = get_rupture_plausibilities(
+            rup_df,
+            distance_matrix=dist,
+            connection_distance_midpoint=5.0,
+            connection_angle_midpoint=None,
+            slip_azimuth_midpoint=None,
+        )
+        one_jump = np.exp(-np.log(2.0) / 5.0 * 8.0)
+        np.testing.assert_allclose(
+            plaus["connection_distance"].values, [one_jump**2, one_jump]
+        )
+
+
 class TestHelperComputations(unittest.TestCase):
     def setUp(self):
         self.rupture_df = pd.DataFrame(