Fix overflow in point-set distance norms

src/pyrecest/utils/metrics.py

@@ -661,7 +661,23 @@ def _pairwise_distances(
 ) -> np.ndarray:
     if distance_fn is None:
         differences = estimated[:, None, :] - reference[None, :, :]
-        return np.linalg.norm(differences, axis=-1)
+        scales = np.max(np.abs(differences), axis=-1)
+        distances = scales.copy()
+        finite_nonzero = np.isfinite(scales) & (scales > 0.0)
+        if np.any(finite_nonzero):
+            normalized = np.zeros_like(differences)
+            np.divide(
+                differences,
+                scales[..., None],
+                out=normalized,
+                where=finite_nonzero[..., None],
+            )
+            normalized_norms = np.linalg.norm(normalized, axis=-1)
+            distances[finite_nonzero] = (
+                scales[finite_nonzero] * normalized_norms[finite_nonzero]
+            )
+        distances[scales == 0.0] = 0.0
+        return distances
     distances = np.empty((estimated.shape[0], reference.shape[0]), dtype=float)
     for row, estimate in enumerate(estimated):
         for col, truth in enumerate(reference):

tests/test_metrics_extreme_coordinates.py

@@ -0,0 +1,52 @@
+import numpy as np
+
+from pyrecest.utils.metrics import gospa_distance, ospa_distance
+
+
+def _large_coordinate_case():
+    component = 1.0e308
+    estimated = np.array([[component, component]])
+    reference = np.zeros((1, 2))
+    expected = np.hypot(component, component)
+    cutoff = 1.5e308
+    return estimated, reference, expected, cutoff
+
+
+def test_ospa_preserves_representable_large_coordinate_distance():
+    estimated, reference, expected, cutoff = _large_coordinate_case()
+
+    with np.errstate(over="raise", invalid="raise"):
+        distance = ospa_distance(
+            estimated,
+            reference,
+            cutoff=cutoff,
+            order=1.0,
+        )
+
+    assert np.isfinite(distance)
+    np.testing.assert_allclose(
+        distance,
+        expected,
+        rtol=4.0 * np.finfo(float).eps,
+        atol=0.0,
+    )
+
+
+def test_gospa_preserves_representable_large_coordinate_distance():
+    estimated, reference, expected, cutoff = _large_coordinate_case()
+
+    with np.errstate(over="raise", invalid="raise"):
+        distance = gospa_distance(
+            estimated,
+            reference,
+            cutoff=cutoff,
+            order=1.0,
+        )
+
+    assert np.isfinite(distance)
+    np.testing.assert_allclose(
+        distance,
+        expected,
+        rtol=4.0 * np.finfo(float).eps,
+        atol=0.0,
+    )