Avoid diagonal-only hypertoroidal particle initialization

src/pyrecest/filters/hypertoroidal_particle_filter.py

@@ -2,17 +2,16 @@
 from typing import Union
 
 import numpy as np
+from scipy.stats import qmc
 
 # pylint: disable=redefined-builtin,no-name-in-module,no-member
 # pylint: disable=no-name-in-module,no-member
 from pyrecest.backend import (
-    arange,
     int32,
     int64,
     linspace,
     mod,
     pi,
-    tile,
 )
 from pyrecest.distributions import (
     AbstractHypertoroidalDistribution,
@@ -57,6 +56,22 @@ def _validate_positive_integer(value, name: str) -> int:
     return integer
 
 
+def _initial_particle_locations(n_particles: int, dim: int):
+    if dim == 1:
+        return linspace(0.0, 2.0 * pi, num=n_particles, endpoint=False)
+
+    # Tiling one circular grid across dimensions confines every particle to the
+    # diagonal x_1 = ... = x_dim. A deterministic Latin hypercube preserves one
+    # evenly spaced angular sample per marginal bin without imposing that
+    # artificial perfect dependence between coordinates.
+    unit_hypercube_points = qmc.LatinHypercube(
+        d=dim,
+        scramble=False,
+        seed=0,
+    ).random(n=n_particles)
+    return 2.0 * np.pi * unit_hypercube_points
+
+
 class HypertoroidalParticleFilter(AbstractParticleFilter, HypertoroidalFilterMixin):
     def __init__(
         self,
@@ -65,12 +80,7 @@ def __init__(
     ):
         n_particles = _validate_positive_integer(n_particles, "n_particles")
         dim = _validate_positive_integer(dim, "dim")
-        if dim == 1:
-            points = linspace(0.0, 2.0 * pi, num=n_particles, endpoint=False)
-        else:
-            points = tile(
-                arange(0.0, 2.0 * pi, 2.0 * pi / n_particles), (dim, 1)
-            ).T.squeeze()
+        points = _initial_particle_locations(n_particles, dim)
         filter_state = HypertoroidalDiracDistribution(points, dim=dim)
         HypertoroidalFilterMixin.__init__(self)
         AbstractParticleFilter.__init__(self, filter_state)

tests/filters/test_hypertoroidal_particle_filter.py

@@ -1,5 +1,6 @@
 import unittest
 
+import numpy as np
 import numpy.testing as npt
 import pyrecest.backend
 
@@ -41,6 +42,36 @@ def test_constructor_rejects_invalid_dimension(self):
                 with self.assertRaisesRegex(ValueError, "dim"):
                     HypertoroidalParticleFilter(5, dim)
 
+    def test_multidimensional_constructor_does_not_collapse_to_diagonal(self):
+        n_particles = 32
+        particles = np.asarray(
+            HypertoroidalParticleFilter(n_particles, 3).filter_state.d
+        )
+
+        self.assertEqual(particles.shape, (n_particles, 3))
+        self.assertFalse(np.allclose(particles[:, 0], particles[:, 1]))
+
+        expected_spacing = 2.0 * np.pi / n_particles
+        for dimension in range(particles.shape[1]):
+            sorted_angles = np.sort(particles[:, dimension])
+            wrapped_angles = np.concatenate(
+                (sorted_angles, [sorted_angles[0] + 2.0 * np.pi])
+            )
+            npt.assert_allclose(np.diff(wrapped_angles), expected_spacing)
+
+    def test_constructor_preserves_requested_particle_count(self):
+        hpf = HypertoroidalParticleFilter(61, 2)
+
+        self.assertEqual(hpf.filter_state.d.shape, (61, 2))
+        self.assertEqual(hpf.filter_state.w.shape, (61,))
+
+    def test_constructor_preserves_singleton_particle_axis(self):
+        hpf = HypertoroidalParticleFilter(1, 3)
+
+        self.assertEqual(hpf.filter_state.d.shape, (1, 3))
+        self.assertEqual(hpf.filter_state.w.shape, (1,))
+        self.assertEqual(hpf.get_point_estimate().shape, (3,))
+
     @unittest.skipIf(
         pyrecest.backend.__backend_name__ == "jax", reason="Backend not supported'"
     )