diff --git a/pymskt/__init__.py b/pymskt/__init__.py index c392fc4..84f74d6 100644 --- a/pymskt/__init__.py +++ b/pymskt/__init__.py @@ -12,4 +12,4 @@ RTOL = 1e-4 ATOL = 1e-5 -__version__ = "0.1.20" +__version__ = "0.1.21" diff --git a/pymskt/mesh/meshTools.py b/pymskt/mesh/meshTools.py index be93774..4dc6877 100644 --- a/pymskt/mesh/meshTools.py +++ b/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) diff --git a/pymskt/mesh/meshes.py b/pymskt/mesh/meshes.py index 6faf901..8da294b 100644 --- a/pymskt/mesh/meshes.py +++ b/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): diff --git a/testing/mesh/meshTools/rand_sample_pts_mesh_seed_test.py b/testing/mesh/meshTools/rand_sample_pts_mesh_seed_test.py new file mode 100644 index 0000000..248d876 --- /dev/null +++ b/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