Solve closest points once per fault pair in rupture propagation

source_modelling/rupture_propagation.py

@@ -12,6 +12,7 @@
 [0]: https://github.com/ucgmsim/source_modelling/wiki/Rupture-Propagation
 """
 
+import itertools
 import random
 import warnings
 from collections import defaultdict, namedtuple
@@ -459,18 +460,18 @@ def sample_rupture_propagation(
         raise ValueError(
             "Cannot specify an initial source and initial source distribution at the same time."
         )
-    distance_graph = {
-        source_a_name: {
-            source_b_name: distance_between(
-                source_a,
-                source_b,
-                *sources.closest_point_between_sources(source_a, source_b),
-            )
-            for source_b_name, source_b in sources_map.items()
-            if source_a_name != source_b_name
-        }
-        for source_a_name, source_a in sources_map.items()
-    }
+    # Distance is symmetric, so solve for the closest points once per pair.
+    distance_graph: DistanceGraph = {source_name: {} for source_name in sources_map}
+    for (source_a_name, source_a), (source_b_name, source_b) in itertools.combinations(
+        sources_map.items(), 2
+    ):
+        distance = distance_between(
+            source_a,
+            source_b,
+            *sources.closest_point_between_sources(source_a, source_b),
+        )
+        distance_graph[source_a_name][source_b_name] = distance
+        distance_graph[source_b_name][source_a_name] = distance
     # Prune the distance graph to remove physically impossible jumps.
     distance_graph = prune_distance_graph(
         distance_graph, jump_impossibility_limit_distance

tests/test_rupture_propagation.py

@@ -501,6 +501,37 @@ def test_sample_rupture_propagation(
                 assert distance <= jump_impossibility_limit_distance
 
 
+def test_sample_rupture_propagation_solves_each_pair_once(
+    monkeypatch: pytest.MonkeyPatch,
+):
+    """Closest points should be solved once per unordered pair of sources."""
+    sources_map = {
+        name: sources.Point(np.array([-41.2865, lon, 0]), 1000, 1000, 0, 0, 0)
+        for name, lon in zip("ABCD", [174.7762, 174.7862, 174.7962, 174.8062])
+    }
+    closest_point_between_sources = sources.closest_point_between_sources
+    calls = []
+
+    def counting_closest_point_between_sources(
+        source_a: sources.IsSource, source_b: sources.IsSource
+    ) -> tuple[np.ndarray, np.ndarray]:
+        calls.append((source_a, source_b))
+        return closest_point_between_sources(source_a, source_b)
+
+    monkeypatch.setattr(
+        sources,
+        "closest_point_between_sources",
+        counting_closest_point_between_sources,
+    )
+    tree = rupture_propagation.sample_rupture_propagation(
+        sources_map,  # ty: ignore[invalid-argument-type]
+        initial_source="A",
+    )
+    n = len(sources_map)
+    assert len(calls) == n * (n - 1) // 2
+    assert set(tree) == set(sources_map)
+
+
 @pytest.mark.parametrize(
     "source_map, rupture_causality_tree, min_depth, expected_jump_points",
     [