From d633f103d9c94ccda8e2cb7dd0887ff9a1eafb53 Mon Sep 17 00:00:00 2001
From: Daniel Fremont <dfremont@ucsc.edu>
Date: Sat, 11 Jan 2025 15:26:28 -0800
Subject: [PATCH] fix reproducibility issue

---
 src/scenic/core/utils.py           |  7 +++-
 tests/syntax/test_distributions.py | 54 +++++++++++++++++++++++++-----
 2 files changed, 51 insertions(+), 10 deletions(-)

diff --git a/src/scenic/core/utils.py b/src/scenic/core/utils.py
index 68ab3c374..9817843f5 100644
--- a/src/scenic/core/utils.py
+++ b/src/scenic/core/utils.py
@@ -320,7 +320,12 @@ def findMeshInteriorPoint(mesh, num_samples=None):
             num_samples = 1
         else:
             num_samples = math.ceil(min(1e6, max(1, math.log(0.01, 1 - p_volume))))
-    samples = trimesh.sample.volume_mesh(mesh, num_samples)
+
+    # Do the "random" number generation ourselves so that it's deterministic
+    # (this helps debugging and reproducibility)
+    rng = numpy.random.default_rng(49493130352093220597973654454967996892)
+    pts = (rng.random((num_samples, 3)) * mesh.extents) + mesh.bounds[0]
+    samples = pts[mesh.contains(pts)]
     if samples.size > 0:
         return samples[0]
 
diff --git a/tests/syntax/test_distributions.py b/tests/syntax/test_distributions.py
index c45572988..6b983c5d2 100644
--- a/tests/syntax/test_distributions.py
+++ b/tests/syntax/test_distributions.py
@@ -668,18 +668,54 @@ def test_reproducibility():
             assert iterations == baseIterations
 
 
+def test_reproducibility_lazy_interior():
+    """Regression test for a reproducibility issue involving lazy region computations.
+
+    In this test, an interior point of the objects' shape is computed on demand
+    during the first sample, then cached. The code for doing so used NumPy's PRNG,
+    meaning that a second sample with the same random seed could differ.
+    """
+    scenario = compileScenic(
+        """
+        import numpy
+        from scenic.core.distributions import distributionFunction
+        @distributionFunction
+        def gen(arg):
+            return numpy.random.random()
+
+        region = BoxRegion().difference(BoxRegion(dimensions=(0.1, 0.1, 0.1)))
+        shape = MeshShape(region.mesh)  # Shape which does not contain its center
+        other = new Object with shape shape
+        ego = new Object at (Range(0.9, 1.1), 0), with shape shape
+        param foo = ego intersects other  # trigger computation of interior point
+        param bar = gen(ego)  # generate number using NumPy's PRNG
+        """
+    )
+    seed = random.randint(0, 1000000000)
+    random.seed(seed)
+    numpy.random.seed(seed)
+    s1 = sampleScene(scenario, maxIterations=60)
+    random.seed(seed)
+    numpy.random.seed(seed)
+    s2 = sampleScene(scenario, maxIterations=60)
+    assert s1.params["bar"] == s2.params["bar"]
+    assert s1.egoObject.x == s2.egoObject.x
+
+
 @pytest.mark.slow
 def test_reproducibility_3d():
     scenario = compileScenic(
-        "ego = new Object\n"
-        "workspace = Workspace(SpheroidRegion(dimensions=(25,15,10)))\n"
-        "region = BoxRegion(dimensions=(25,15,0.1))\n"
-        "obj_1 = new Object in workspace, facing Range(0, 360) deg, with width Range(0.5, 1), with length Range(0.5,1)\n"
-        "obj_2 = new Object in workspace, facing (Range(0, 360) deg, Range(0, 360) deg, Range(0, 360) deg)\n"
-        "obj_3 = new Object in workspace, on region\n"
-        "param foo = Uniform(1, 4, 9, 16, 25, 36)\n"
-        "x = Range(0, 1)\n"
-        "require x > 0.8"
+        """
+        ego = new Object
+        workspace = Workspace(SpheroidRegion(dimensions=(5,5,5)))
+        region = BoxRegion(dimensions=(25,15,0.1))
+        #obj_1 = new Object in workspace, facing Range(0, 360) deg, with width Range(0.5, 1), with length Range(0.5,1)
+        obj_2 = new Object in workspace, facing (Range(0, 360) deg, Range(0, 360) deg, Range(0, 360) deg)
+        #obj_3 = new Object in workspace, on region
+        param foo = ego intersects obj_2
+        x = Range(0, 1)
+        require x > 0.8
+        """
     )
     seeds = [random.randint(0, 100) for i in range(10)]
     for seed in seeds: