Introduce resampler that checks if resampling required

stonesoup/resampler/particle.py

@@ -2,6 +2,7 @@
 import numpy as np
 
 from .base import Resampler
+from ..base import Property
 from ..types.numeric import Probability
 from ..types.particle import Particles
 
@@ -25,7 +26,7 @@ def resample(self, particles):
         if not isinstance(particles, Particles):
             particles = Particles(particle_list=particles)
         n_particles = len(particles)
-        weight = Probability(1/n_particles)
+        weight = Probability(1 / n_particles)
 
         log_weights = np.array([weight.log_value for weight in particles.weight])
         weight_order = np.argsort(log_weights, kind='stable')
@@ -42,7 +43,51 @@ def resample(self, particles):
         u_j = u_i + (1 / n_particles) * np.arange(n_particles)
         index = weight_order[np.searchsorted(cdf, np.log(u_j))]
         new_particles = Particles(state_vector=particles.state_vector[:, index],
-                                  weight=[weight]*n_particles,
+                                  weight=[weight] * n_particles,
                                   parent=Particles(state_vector=particles.state_vector[:, index],
                                                    weight=particles.weight[index]))
         return new_particles
+
+
+class ESSResampler(Resampler):
+    """ESSResampler class
+        This wrapper uses a :class:`~.Resampler` to resample the particles inside
+        an instant of :class:`~.Particles`, but only after checking if this is necessary
+        by comparing Effective Sample Size (ESS) with a supplied threshold (numeric).
+        Kish's ESS is used, as in Section 3.5 of [1]_.   Doucet-Johansen, 2008.
+
+        References
+        ----------
+        .. [1] Doucet A., Johansen A.M., 2009, Tutorial on Particle Filtering \
+        and Smoothing: Fifteen years later, Handbook of Nonlinear Filtering, Vol. 12.
+
+        """
+
+    threshold: float = Property(default=None,
+                                doc='Threshold compared with ESS to decide whether to resample. \
+                                    Default is number of particles divided by 2, \
+                                        set in resample method')
+    resampler: Resampler = Property(default=SystematicResampler,
+                                    doc='Resampler to wrap, which is called \
+                                        when ESS below threshold')
+
+    def resample(self, particles):
+        """
+        Parameters
+        ----------
+        particles : list of :class:`~.Particle`
+            The particles to be resampled according to their weight
+
+        Returns
+        -------
+        particles : list of :class:`~.Particle`
+            The particles, either unchanged or resampled, depending on weight degeneracy
+        """
+        if not isinstance(particles, Particles):
+            particles = Particles(particle_list=particles)
+        if self.threshold is None:
+            self.threshold = len(particles) / 2
+        if 1 / np.sum(np.square(particles.weight)) < self.threshold:  # If ESS too small, resample
+            return self.resampler.resample(self.resampler, particles)
+        else:
+            return particles

stonesoup/resampler/tests/test_particle.py

@@ -1,16 +1,18 @@
 # -*- coding: utf-8 -*-
 import numpy as np
 
-from ...types.particle import Particle
+from ...types.particle import Particle, \
+    Particles
 from ..particle import SystematicResampler
+from ..particle import ESSResampler
 
 
 def test_systematic_equal():
     particles = [Particle(np.array([[i]]), weight=1/20) for i in range(20)]
 
-    resamlper = SystematicResampler()
+    resampler = SystematicResampler()
 
-    new_particles = resamlper.resample(particles)
+    new_particles = resampler.resample(particles)
 
     # Particles equal weight, so should end up with similar particles.
     assert all(np.array_equal(np.array([[i]]), new_particle.state_vector)
@@ -21,9 +23,9 @@ def test_systematic_single():
     particles = [Particle(np.array([[i]]), weight=1 if i == 10 else 0)
                  for i in range(20)]
 
-    resamlper = SystematicResampler()
+    resampler = SystematicResampler()
 
-    new_particles = resamlper.resample(particles)
+    new_particles = resampler.resample(particles)
 
     # Weight all at single particle at index/vector 10, so all should be
     # at same point
@@ -35,10 +37,55 @@ def test_systematic_even():
     particles = [Particle(np.array([[i]]), weight=1/10 if i % 2 == 0 else 0)
                  for i in range(20)]
 
-    resamlper = SystematicResampler()
+    resampler = SystematicResampler()
 
-    new_particles = resamlper.resample(particles)
+    new_particles = resampler.resample(particles)
 
     # Weight all at even particles, so new should all be at even vector
     assert all(np.array_equal(np.array([[i//2*2]]), new_particle.state_vector)
                for i, new_particle in enumerate(new_particles))
+
+
+def test_ess_zero():
+    particles = [Particle(np.array([[i]]), weight=(i+1) / 55)
+                 for i in range(10)]
+    particles = Particles(particle_list=particles)  # Needed for assert line
+
+    resampler = ESSResampler(0)  # Should never resample
+
+    new_particles = resampler.resample(particles)  # Should not have changed
+
+    assert new_particles == particles
+
+
+def test_ess_n():
+    particles = [Particle(np.array([[i]]), weight=(i+1) / 55)
+                 for i in range(10)]
+
+    resampler = ESSResampler(10)  # This resampler should always resample
+
+    new_particles = resampler.resample(particles)  # All weights should be equal due to resampling
+
+    assert all([w == 1/10 for w in new_particles.weight])
+
+
+def test_ess_inequality():
+    particles = [Particle(np.array([[i]]), weight=(i + 1) / 55)
+                 for i in range(10)]
+    particles = Particles(particle_list=particles)
+
+    resampler1 = ESSResampler(3025/385)  # This resampler should not resample
+    resampler2 = ESSResampler(3025/385 + 0.01)  # This resampler should resample
+
+    new_particles1 = resampler1.resample(particles)
+    new_particles2 = resampler2.resample(particles)
+
+    assert new_particles1 == particles and all([w == 1 / 10 for w in new_particles2.weight])
+
+
+def test_ess_empty_threshold():
+    particles = [Particle(np.array([[i]]), weight=(i + 1) / 55)
+                 for i in range(10)]
+    resampler = ESSResampler()
+    resampler.resample(particles)
+    assert resampler.threshold == 5