From 8a8dff4cab320a5fefa299b443bc1be2563015d6 Mon Sep 17 00:00:00 2001 From: Marco Pagani Date: Sat, 17 Oct 2026 22:16:55 +0200 Subject: [PATCH] fnm: evaluate rupture jumps on the actual connections Distances and angles of a multifault rupture were read between consecutive indices of its rupture list, which is not ordered along the connections: unconnected pairs were evaluated (as 0 km) and real jumps skipped. Take the connections from the adjacency structure (minimum spanning tree of the rupture subgraph, new get_multifault_connection_pairs) and use the same pairs for the distance and angle plausibilities. --- openquake/fnm/rupture_connections.py | 131 ++++++++++++++++-- openquake/fnm/rupture_filtering.py | 40 ++++-- .../fnm/tests/test_rupture_connections.py | 27 ++++ openquake/fnm/tests/test_rupture_filtering.py | 23 +++ 4 files changed, 200 insertions(+), 21 deletions(-) diff --git a/openquake/fnm/rupture_connections.py b/openquake/fnm/rupture_connections.py index 93c94a3db..05be221ec 100644 --- a/openquake/fnm/rupture_connections.py +++ b/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, ) diff --git a/openquake/fnm/rupture_filtering.py b/openquake/fnm/rupture_filtering.py index 5af3fe620..93f9af52a 100644 --- a/openquake/fnm/rupture_filtering.py +++ b/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( diff --git a/openquake/fnm/tests/test_rupture_connections.py b/openquake/fnm/tests/test_rupture_connections.py index 348e2706d..58abc6f6f 100644 --- a/openquake/fnm/tests/test_rupture_connections.py +++ b/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], + ) diff --git a/openquake/fnm/tests/test_rupture_filtering.py b/openquake/fnm/tests/test_rupture_filtering.py index a3cc45bf1..3ae26e368 100644 --- a/openquake/fnm/tests/test_rupture_filtering.py +++ b/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(