Stabilize nonadditive particle-noise weight normalization

src/pyrecest/filters/abstract_particle_filter.py

@@ -8,6 +8,7 @@
     array,
     hstack,
     isfinite,
+    max as backend_max,
     ndim,
     ones_like,
     random,
@@ -230,10 +231,18 @@ def predict_nonlinear_nonadditive(self, f, samples, weights):
             raise ValueError("Noise weights must be finite.")
         if not bool(all(weights >= 0.0)):
             raise ValueError("Noise weights must be nonnegative.")
-        weight_sum = sum(weights)
-        if not bool(isfinite(weight_sum)) or not bool(weight_sum > 0.0):
+        if weights.shape[0] == 0:
             raise ValueError("Noise weights must have positive finite total mass.")
-        weights = weights / weight_sum
+        weight_scale = backend_max(weights)
+        if not bool(isfinite(weight_scale)) or not bool(weight_scale > 0.0):
+            raise ValueError("Noise weights must have positive finite total mass.")
+        scaled_weights = weights / weight_scale
+        scaled_weight_sum = sum(scaled_weights)
+        if not bool(isfinite(scaled_weight_sum)) or not bool(
+            scaled_weight_sum > 0.0
+        ):
+            raise ValueError("Noise weights must have positive finite total mass.")
+        weights = scaled_weights / scaled_weight_sum
         n_particles = self.filter_state.w.shape[0]
         noise_samples = random.choice(samples, n_particles, p=weights)
 
@@ -275,7 +284,7 @@ def update_model(
         """Update using a reusable particle measurement model."""
         if not hasattr(measurement_model, "likelihood"):
             raise TypeError(
-                "Particle-filter measurement models must expose a likelihood callable."
+                "Particle measurement models must expose a likelihood callable."
             )
 
         return self.update_nonlinear_using_likelihood(

tests/filters/test_particle_filter_extreme_nonadditive_weights.py

@@ -0,0 +1,32 @@
+from unittest import mock
+
+import numpy as np
+import numpy.testing as npt
+
+from pyrecest.backend import array, to_numpy
+from pyrecest.distributions import LinearDiracDistribution
+from pyrecest.filters.euclidean_particle_filter import EuclideanParticleFilter
+
+
+def test_nonadditive_prediction_normalizes_extreme_finite_noise_weights():
+    particles = array([[10.0], [20.0]])
+    samples = array([[1.0], [2.0]])
+    backend_dtype = to_numpy(array([1.0])).dtype
+    max_weight = np.finfo(backend_dtype).max
+    weights = array([max_weight, max_weight / 2.0])
+    particle_filter = EuclideanParticleFilter(n_particles=2, dim=1)
+    particle_filter.filter_state = LinearDiracDistribution(particles)
+
+    with mock.patch(
+        "pyrecest.filters.abstract_particle_filter.random.choice",
+        return_value=samples,
+    ) as choice_mock:
+        particle_filter.predict_nonlinear_nonadditive(
+            lambda particle, noise: particle + noise,
+            samples,
+            weights,
+        )
+
+    normalized_weights = choice_mock.call_args.kwargs["p"]
+    npt.assert_allclose(to_numpy(normalized_weights), [2.0 / 3.0, 1.0 / 3.0])
+    npt.assert_allclose(to_numpy(particle_filter.filter_state.d), [[11.0], [22.0]])