Use numpy's random Generator interface.

* sambal/sambal.py: Import numpy.random.default_rng. Do not import
numpy.random.randn, numpy.random.random.
(random_on_sphere, random_on_disk, random_on_cap): Accept optional rng
argument. Use rng to generate random numbers.

1 files changed, 9 insertions, 9 deletions

diff --git a/sambal/sambal.py b/sambal/sambal.py
index f4d64f4..836d0b6 100644
--- a/sambal/sambal.py
+++ b/sambal/sambal.py
@@ -17,12 +17,12 @@
 # <https://www.gnu.org/licenses/>.
 
 from numpy import cos, dot, empty, log, sin, sqrt, pi
-from numpy.random import randn, random
+from numpy.random import default_rng
 from numpy.linalg import norm
 
-def random_on_sphere(dim):
+def random_on_sphere(dim, rng=default_rng()):
     """Return a random vector uniformly distributed on the unit sphere."""
-    x = randn(dim)
+    x = rng.standard_normal(dim)
     return x / norm(x)
 
 def rotate_from_nth_canonical(x, axis):
@@ -40,18 +40,18 @@ def rotate_from_nth_canonical(x, axis):
             - axisn*(xn*s + a*(axisn - 1))/b
     return x
 
-def random_on_disk(axis, planar_angle):
+def random_on_disk(axis, planar_angle, rng=default_rng()):
     """Return a random vector uniformly distributed on the periphery of a
 disk.
 
     """
     dim = axis.size
     x = empty(dim)
-    x[:-1] = sin(planar_angle) * random_on_sphere(dim - 1)
+    x[:-1] = sin(planar_angle) * random_on_sphere(dim - 1, rng)
     x[-1] = cos(planar_angle)
     return rotate_from_nth_canonical(x, axis)
 
-def random_on_cap(axis, maximum_planar_angle):
+def random_on_cap(axis, maximum_planar_angle, rng=default_rng()):
     """Return a random vector uniformly distributed on a spherical
 cap. The random planar angle is generated using rejection sampling.
 
@@ -61,8 +61,8 @@ cap. The random planar angle is generated using rejection sampling.
     dim = axis.size
     box_height = (dim-2)*log(sin(min(maximum_planar_angle, pi/2)))
     while True:
-        theta = maximum_planar_angle*random()
-        f = box_height + log(random())
+        theta = maximum_planar_angle*rng.random()
+        f = box_height + log(rng.random())
         if f < (dim-2)*log(sin(theta)):
             break
-    return random_on_disk(axis, theta)
+    return random_on_disk(axis, theta, rng)