Make the surface samplers seedable (closes #54)

pymskt/__init__.py

@@ -12,4 +12,4 @@
 
 RTOL = 1e-4
 ATOL = 1e-5
-__version__ = "0.1.20"
+__version__ = "0.1.21"

pymskt/mesh/meshTools.py

@@ -1530,19 +1530,52 @@ def consistent_normals(mesh, points_dtype=np.float64, faces_dtype=np.int32):
     return new_mesh
 
 
+def pcu_random_seed(seed):
+    """
+    Derive a seed `point_cloud_utils` will honour, from one a user would write.
+
+    pcu's `random_seed` is not an ordinary seed argument, so a user seed cannot be passed
+    through. Three distinct problems, which is why this is a derivation and not a special
+    case:
+
+    - `random_seed=0` means "seed from the current time", not "seed 0". Passing a user's
+      0 through would silently mean "unseeded" -- and 0 is the first seed most people try.
+    - Remapping only 0 to some constant would make that constant collide with a real seed:
+      `if seed == 0: seed = 1` gives seeds 0 and 1 the same run.
+    - The binding is a 32-bit C++ signature and raises TypeError outside that range, so a
+      timestamp or a hash used as a seed fails.
+
+    Deriving a fresh non-zero 32-bit value from `seed` handles all three under one rule.
+    Negative seeds raise, inheriting numpy's contract rather than inventing another.
+
+    `None` maps to pcu's own 0 sentinel, leaving the unseeded default untouched.
+    """
+    if seed is None:
+        return 0
+    return int(np.random.default_rng(seed).integers(1, 2**31 - 1))
+
+
 def rand_sample_pts_mesh(
-    mesh, n_pts, method="bluenoise", points_dtype=np.float64, faces_dtype=np.int32
+    mesh, n_pts, method="bluenoise", points_dtype=np.float64, faces_dtype=np.int32, seed=None
 ):
     """
     Randomly sample points from a mesh
+
+    Args:
+        seed (int, optional): Makes the sampling reproducible. Defaults to None, which
+            samples differently on every call -- the historical behaviour.
     """
     # get faces and points
     faces, points = get_faces_vertices(mesh, points_dtype=points_dtype, faces_dtype=faces_dtype)
 
+    random_seed = pcu_random_seed(seed)
+
     if method == "random":
-        fid, bc = pcu.sample_mesh_random(points, faces, n_pts)
+        fid, bc = pcu.sample_mesh_random(points, faces, n_pts, random_seed=random_seed)
     elif method == "bluenoise":
-        fid, bc = pcu.sample_mesh_poisson_disk(points, faces, num_samples=n_pts)
+        fid, bc = pcu.sample_mesh_poisson_disk(
+            points, faces, num_samples=n_pts, random_seed=random_seed
+        )
 
     rand_pts = pcu.interpolate_barycentric_coords(faces, fid, bc, points)
 

pymskt/mesh/meshes.py

@@ -346,24 +346,43 @@ def get_largest(self):
         print(mesh_.n_points)
         self.deep_copy(mesh_)
 
-    def rand_surface_pts(self, n_pts=100_000, method="bluenoise"):
+    def rand_surface_pts(self, n_pts=100_000, method="bluenoise", seed=None):
         """
         Sample points from the surface of the mesh.
+
+        Args:
+            seed (int, optional): Makes the sampling reproducible. Defaults to None, which
+                samples differently on every call -- the historical behaviour.
         """
-        return rand_sample_pts_mesh(self, n_pts=n_pts, method=method)
+        return rand_sample_pts_mesh(self, n_pts=n_pts, method=method, seed=seed)
 
     def rand_pts_around_surface(
-        self, n_pts=100_000, surface_method="bluenoise", distribution="normal", sigma=1.0
+        self,
+        n_pts=100_000,
+        surface_method="bluenoise",
+        distribution="normal",
+        sigma=1.0,
+        seed=None,
     ):
         """
         Sample points around the surface of the mesh. For SDF sampling & neural implicit representation models.
+
+        Args:
+            seed (int, optional): Makes the sampling reproducible. Defaults to None, which
+                samples differently on every call -- the historical behaviour.
+
+        Notes:
+            Both draws here are seeded: the base surface points (via `rand_surface_pts`)
+            and the offsets applied to them. Seeding only one leaves the result random.
         """
+        rng = np.random.default_rng(seed)
+
         if distribution == "normal":
-            rand_gen = np.random.default_rng().multivariate_normal
+            rand_gen = rng.multivariate_normal
         elif distribution == "laplace":
-            rand_gen = np.random.default_rng().laplace
+            rand_gen = rng.laplace
 
-        base_pts = self.rand_surface_pts(n_pts=n_pts, method=surface_method)
+        base_pts = self.rand_surface_pts(n_pts=n_pts, method=surface_method, seed=seed)
         mean = [0, 0, 0]
 
         if (distribution == "normal") and (sigma is not None):

testing/mesh/meshTools/rand_sample_pts_mesh_seed_test.py

@@ -0,0 +1,114 @@
+"""
+Reproducibility of the random surface samplers.
+
+Regression tests for gattia/pymskt#54: ``Mesh.rand_pts_around_surface`` had two
+independent random draws that a caller could not reach, so identical inputs produced
+different point clouds on every call.
+
+* the base surface points went to ``pcu.sample_mesh_random`` / ``sample_mesh_poisson_disk``
+  with no ``random_seed``, and pcu's ``random_seed=0`` default means "seed from the current
+  time", not "seed 0";
+* the offsets used ``np.random.default_rng()`` with no argument, which seeds itself from OS
+  entropy and ignores ``np.random.seed()``.
+
+Seeding only one of the two leaves the result random, so every test here checks the whole
+call rather than a single draw.
+"""
+
+import numpy as np
+import pytest
+import pyvista as pv
+
+from pymskt.mesh import Mesh
+from pymskt.mesh.meshTools import pcu_random_seed, rand_sample_pts_mesh
+
+
+@pytest.fixture
+def sphere(tmp_path):
+    path = str(tmp_path / "sphere.vtk")
+    pv.Sphere(radius=1.0, theta_resolution=24, phi_resolution=24).triangulate().save(path)
+    return path
+
+
+def around(path, seed, method="random", distribution="normal", n_pts=500):
+    return Mesh(path).rand_pts_around_surface(
+        n_pts=n_pts, surface_method=method, distribution=distribution, sigma=0.01, seed=seed
+    )
+
+
+class TestRandPtsAroundSurface:
+    def test_the_same_seed_gives_the_same_points(self, sphere):
+        assert np.array_equal(around(sphere, 42), around(sphere, 42))
+
+    def test_different_seeds_give_different_points(self, sphere):
+        assert not np.array_equal(around(sphere, 42), around(sphere, 7))
+
+    def test_seed_zero_is_reproducible(self, sphere):
+        """
+        The one that would slip through. pcu reads ``random_seed=0`` as "use the clock", so
+        passing a user's seed straight through would make ``seed=0`` silently mean
+        "unseeded" -- and 0 is the first seed most people try.
+        """
+        assert np.array_equal(around(sphere, 0), around(sphere, 0))
+
+    def test_the_default_is_still_unseeded(self, sphere):
+        """
+        Backwards compatibility: ``seed=None`` must behave exactly as before, including
+        being immune to ``np.random.seed()``, which never reached either draw.
+        """
+        np.random.seed(0)
+        first = around(sphere, None)
+        np.random.seed(0)
+        second = around(sphere, None)
+        assert not np.array_equal(first, second)
+
+    def test_the_laplace_distribution_is_seeded_too(self, sphere):
+        assert np.array_equal(
+            around(sphere, 5, distribution="laplace"),
+            around(sphere, 5, distribution="laplace"),
+        )
+
+
+class TestRandSamplePtsMesh:
+    """
+    The surface-point draw on its own. Both methods are covered here; ``bluenoise`` is not
+    covered end-to-end through ``rand_pts_around_surface`` because that path is broken for
+    an unrelated reason -- ``sample_mesh_poisson_disk`` returns approximately, not exactly,
+    ``num_samples`` points, so adding the offsets raises a broadcast error. That is a
+    separate defect and is not addressed here.
+    """
+
+    @pytest.mark.parametrize("method", ["random", "bluenoise"])
+    def test_the_same_seed_gives_the_same_points(self, sphere, method):
+        mesh = Mesh(sphere)
+        first = rand_sample_pts_mesh(mesh, n_pts=400, method=method, seed=11)
+        second = rand_sample_pts_mesh(mesh, n_pts=400, method=method, seed=11)
+        assert np.array_equal(first, second)
+
+    @pytest.mark.parametrize("method", ["random", "bluenoise"])
+    def test_different_seeds_give_different_points(self, sphere, method):
+        mesh = Mesh(sphere)
+        first = rand_sample_pts_mesh(mesh, n_pts=400, method=method, seed=11)
+        second = rand_sample_pts_mesh(mesh, n_pts=400, method=method, seed=12)
+        assert first.shape != second.shape or not np.array_equal(first, second)
+
+
+class TestPcuRandomSeed:
+    def test_none_maps_to_pcus_unseeded_sentinel(self):
+        assert pcu_random_seed(None) == 0
+
+    def test_zero_does_not_map_to_zero(self):
+        """Otherwise ``seed=0`` would mean "seed from the clock" inside pcu."""
+        assert pcu_random_seed(0) != 0
+
+    def test_it_is_deterministic(self):
+        assert pcu_random_seed(3) == pcu_random_seed(3)
+
+    def test_distinct_seeds_map_to_distinct_values(self):
+        assert len({pcu_random_seed(i) for i in range(25)}) == 25
+
+    def test_the_result_is_a_positive_32_bit_int(self):
+        for seed in (0, 1, 12345, 2**31):
+            value = pcu_random_seed(seed)
+            assert isinstance(value, int)
+            assert 0 < value < 2**31