Use seeded np.random.default_rng() random number generator in tests

CHANGES.rst

@@ -20,6 +20,9 @@ New Features
   apo/pericenter (Sanders & Binney 2016). This is much faster than the default
   energy-variance fit (``method="variance"``).
 
+- ``fast_lyapunov_max()`` now accepts an ``rng`` argument (a ``numpy.random.Generator``)
+  for the initial offset vectors, matching ``lyapunov_max()``.
+
 
 Bug fixes
 ---------
@@ -88,6 +91,13 @@ API changes
 - The ``shared`` and ``extra`` optional dependencies are merged into ``extra``
   (``pip install gala[extra]``).
 
+- Mock stream distribution functions (``FardalStreamDF``, etc.) now use a
+  ``numpy.random.Generator`` for ``random_state``. It defaults to
+  ``np.random.default_rng()`` and also accepts an integer seed. Passing a legacy
+  ``numpy.random.RandomState`` still works but is deprecated. The default generator
+  draws different random numbers than ``RandomState``, so unseeded streams will not
+  match those made with earlier versions.
+
 Other
 -----
 

conftest.py

@@ -2,6 +2,8 @@
 import sys
 from pathlib import Path
 
+import numpy as np
+import pytest
 from pytest_astropy_header.display import (
     PYTEST_HEADER_MODULES,
     TESTED_VERSIONS,
@@ -12,6 +14,12 @@
 sys.path.insert(0, str(tests_dir))
 
 
+@pytest.fixture
+def rng():
+    """Seeded random number generator; tests should draw from this, not np.random."""
+    return np.random.default_rng(42)
+
+
 def pytest_configure(config):
     config.option.astropy_header = True
     PYTEST_HEADER_MODULES.pop("Pandas", None)

src/gala/dynamics/lyapunov/dop853_lyapunov.pyx

@@ -24,7 +24,7 @@ cpdef dop853_lyapunov_max(hamiltonian, double[::1] w0,
                           double dt, int n_steps, double t0,
                           double d0, int n_steps_per_pullback, int noffset_orbits,
                           double atol=1E-10, double rtol=1E-10, int nmax=0,
-                          unsigned log_output=0):
+                          unsigned log_output=0, rng=None):
     cdef:
         int i, j, k, jiter
         int res
@@ -44,14 +44,17 @@ cpdef dop853_lyapunov_max(hamiltonian, double[::1] w0,
         double[:, :, ::1] all_w = np.zeros((n_steps, norbits, ndim))
 
         # temp stuff
-        double[:, ::1] d0_vec = np.random.uniform(size=(noffset_orbits, ndim))
+        double[:, ::1] d0_vec
 
         # whoa, so many dots
         CPotential* cp = (<CPotentialWrapper>(hamiltonian.potential.c_instance)).cpotential
         CFrameType cf = (<CFrameWrapper>(hamiltonian.frame.c_instance)).cframe
 
         void *args
 
+    rng = rng or np.random.default_rng()
+    d0_vec = rng.uniform(size=(noffset_orbits, ndim))
+
     # store initial conditions
     for i in range(norbits):
         if i == 0:  # store initial conditions for parent orbit
@@ -104,7 +107,7 @@ cpdef dop853_lyapunov_max_dont_save(hamiltonian, double[::1] w0,
                                     double dt, int n_steps, double t0,
                                     double d0, int n_steps_per_pullback, int noffset_orbits,
                                     double atol=1E-10, double rtol=1E-10, int nmax=0,
-                                    unsigned log_output=0):
+                                    unsigned log_output=0, rng=None):
     cdef:
         int i, j, k, jiter
         int res
@@ -123,14 +126,17 @@ cpdef dop853_lyapunov_max_dont_save(hamiltonian, double[::1] w0,
         double[:, ::1] LEs = np.zeros((niter, noffset_orbits))
 
         # temp stuff
-        double[:, ::1] d0_vec = np.random.uniform(size=(noffset_orbits, ndim))
+        double[:, ::1] d0_vec
 
         # whoa, so many dots
         CPotential* cp = (<CPotentialWrapper>(hamiltonian.potential.c_instance)).cpotential
         CFrameType cf = (<CFrameWrapper>(hamiltonian.frame.c_instance)).cframe
 
         void *args
 
+    rng = rng or np.random.default_rng()
+    d0_vec = rng.uniform(size=(noffset_orbits, ndim))
+
     # store initial conditions
     for i in range(norbits):
         if i == 0:  # store initial conditions for parent orbit

src/gala/dynamics/mockstream/_coord.pyx

@@ -151,7 +151,7 @@ cdef void cyl_to_car(double *cyl, # in
 
 cpdef _test_sat_rotation_matrix():
     import numpy as np
-    np.random.seed(42)
+    rng = np.random.default_rng(42)
     n = 1024
 
     cdef:
@@ -162,7 +162,7 @@ cpdef _test_sat_rotation_matrix():
         unsigned int i, j
 
     for i in range(n):
-        w = np.random.uniform(size=6)
+        w = rng.uniform(size=6)
         sat_rotation_matrix(&w[0], &R[0, 0])
 
         x = np.array(R).dot(np.array(w)[:3])
@@ -187,12 +187,12 @@ cpdef _test_sat_rotation_matrix():
 
 cpdef _test_to_sat_coords_roundtrip():
     import numpy as np
-    np.random.seed(42)
+    rng = np.random.default_rng(42)
     n = 1024
 
     cdef:
-        double[:, ::1] w = np.random.uniform(size=(n, 6))
-        double[:, ::1] w_sat = np.random.uniform(size=(n, 6))
+        double[:, ::1] w = rng.uniform(size=(n, 6))
+        double[:, ::1] w_sat = rng.uniform(size=(n, 6))
         double[:, ::1] R = np.zeros((3, 3))
 
         double[::1] w_prime = np.zeros(6)
@@ -210,11 +210,11 @@ cpdef _test_to_sat_coords_roundtrip():
 
 cpdef _test_car_to_cyl_roundtrip():
     import numpy as np
-    np.random.seed(42)
+    rng = np.random.default_rng(42)
     n = 1024
 
     cdef:
-        double[:, ::1] w = np.random.uniform(-10, 10, size=(n, 6))
+        double[:, ::1] w = rng.uniform(-10, 10, size=(n, 6))
         double[::1] cyl = np.zeros(6)
         double[::1] w2 = np.zeros(6)
 
@@ -228,11 +228,11 @@ cpdef _test_car_to_cyl_roundtrip():
 
 cpdef _test_cyl_to_car_roundtrip():
     import numpy as np
-    # np.random.seed(42)
+    rng = np.random.default_rng(42)
     n = 1024
 
     cdef:
-        double[:, ::1] cyl = np.random.uniform(0, 2*np.pi, size=(n, 6))
+        double[:, ::1] cyl = rng.uniform(0, 2*np.pi, size=(n, 6))
         double[::1] w = np.zeros(6)
         double[::1] cyl2 = np.zeros(6)
 

src/gala/dynamics/mockstream/df.pyx

@@ -41,17 +41,27 @@ cdef class BaseStreamDF:
         Generate a leading tail. Default: True.
     trail : bool (optional)
         Generate a trailing tail. Default: True.
-    random_state : `~numpy.random.RandomState` (optional)
-        To control random number generation.
+    random_state : `~numpy.random.Generator`, int, or None (optional)
+        To control random number generation. Passed to
+        :func:`numpy.random.default_rng`, so an integer seed also works.
 
     """
     def __init__(self, lead=True, trail=True, random_state=None):
 
         self._lead = int(lead)
         self._trail = int(trail)
 
-        if random_state is None:
-            random_state = np.random.RandomState()
+        if isinstance(random_state, np.random.RandomState):
+            import warnings
+            from gala.util import GalaDeprecationWarning
+            warnings.warn(
+                "Passing a numpy.random.RandomState as random_state is deprecated; "
+                "pass a numpy.random.Generator instead, e.g. "
+                "np.random.default_rng(seed).",
+                GalaDeprecationWarning,
+            )
+        else:
+            random_state = np.random.default_rng(random_state)
         self.random_state = random_state
 
         if not self.lead and not self.trail:
@@ -250,8 +260,9 @@ cdef class StreaklineStreamDF(BaseStreamDF):
         Generate a leading tail. Default: True.
     trail : bool (optional)
         Generate a trailing tail. Default: True.
-    random_state : `~numpy.random.RandomState` (optional)
-        To control random number generation.
+    random_state : `~numpy.random.Generator`, int, or None (optional)
+        To control random number generation. Passed to
+        :func:`numpy.random.default_rng`, so an integer seed also works.
     """
 
     cpdef _sample(self, potential,
@@ -332,8 +343,9 @@ cdef class FardalStreamDF(BaseStreamDF):
         Generate a leading tail. Default: True.
     trail : bool (optional)
         Generate a trailing tail. Default: True.
-    random_state : `~numpy.random.RandomState` (optional)
-        To control random number generation.
+    random_state : `~numpy.random.Generator`, int, or None (optional)
+        To control random number generation. Passed to
+        :func:`numpy.random.default_rng`, so an integer seed also works.
     """
     def __init__(
         self, gala_modified=None, lead=True, trail=True, random_state=None
@@ -472,8 +484,9 @@ cdef class LagrangeCloudStreamDF(BaseStreamDF):
         Generate a leading tail. Default: True.
     trail : bool (optional)
         Generate a trailing tail. Default: True.
-    random_state : `~numpy.random.RandomState` (optional)
-        To control random number generation.
+    random_state : `~numpy.random.Generator`, int, or None (optional)
+        To control random number generation. Passed to
+        :func:`numpy.random.default_rng`, so an integer seed also works.
     """
 
     cdef public object v_disp
@@ -565,8 +578,9 @@ cdef class ChenStreamDF(BaseStreamDF):
         Generate a leading tail. Default: True.
     trail : bool (optional)
         Generate a trailing tail. Default: True.
-    random_state : `~numpy.random.RandomState` (optional)
-        To control random number generation.
+    random_state : `~numpy.random.Generator`, int, or None (optional)
+        To control random number generation. Passed to
+        :func:`numpy.random.default_rng`, so an integer seed also works.
     """
     def __init__(
         self, lead=True, trail=True, random_state=None

src/gala/dynamics/nonlinear.py

@@ -21,6 +21,7 @@ def fast_lyapunov_max(
     rtol=1e-10,
     nmax=0,
     return_orbit=True,
+    rng=None,
 ):
     """
     Compute the maximum Lyapunov exponent using a fast C-implemented method.
@@ -59,6 +60,10 @@ def fast_lyapunov_max(
     return_orbit : bool, optional
         Whether to return the integrated orbit along with the Lyapunov
         exponent. Default is True.
+    rng : `numpy.random.Generator` (optional)
+        If provided, will use this random number generator to generate the
+        initial offset vectors. If not provided, will use
+        ``numpy.random.default_rng()``.
 
     Returns
     -------
@@ -121,6 +126,7 @@ def fast_lyapunov_max(
             atol,
             rtol,
             nmax,
+            rng=rng,
         )
         w = np.rollaxis(w, -1)
 
@@ -150,6 +156,7 @@ def fast_lyapunov_max(
         atol,
         rtol,
         nmax,
+        rng=rng,
     )
 
     try:

tests/benchmarks/test_mockstream_benchmark.py

@@ -1,18 +1,11 @@
 import astropy.units as u
-import numpy as np
 import pytest
 
 import gala.dynamics as gd
 import gala.potential as gp
 from gala.dynamics import mockstream as ms
 
 
-@pytest.fixture
-def rng():
-    """Random number generator with fixed seed for reproducibility."""
-    return np.random.default_rng(42)
-
-
 @pytest.fixture
 def hamiltonian():
     """Standard NFW potential for benchmarks."""

tests/coordinates/test_pm_cov_transform.py

@@ -77,7 +77,7 @@ def test_transform_correctness(to_frame):
 
     # generate proper motion samples and transform the samples:
     pm = np.vstack((c.pm_ra_cosdec.value, c.pm_dec.value)).T
-    rnd = np.random.RandomState(42)
+    rnd = np.random.default_rng(42)
 
     for i in range(len(c)):
         pm_samples = rnd.multivariate_normal(pm[i], cov[i], size=2**16)

tests/dynamics/actionangle/actionangle_helpers.py

@@ -116,18 +116,18 @@ def plot_angles(t, angles, freqs, subsample_factor=1000):
     # axes[1].scatter(theta[0, ix], theta[2], alpha=0.5, marker='o', c=t)
 
 
-def isotropic_w0(N=100):
+def isotropic_w0(rng, N=100):
     # positions
-    d = np.random.lognormal(mean=np.log(25), sigma=0.5, size=N)
-    phi = np.random.uniform(0, 2 * np.pi, size=N)
-    theta = np.arccos(np.random.uniform(size=N) - 0.5)
+    d = rng.lognormal(mean=np.log(25), sigma=0.5, size=N)
+    phi = rng.uniform(0, 2 * np.pi, size=N)
+    theta = np.arccos(rng.uniform(size=N) - 0.5)
 
-    vr = np.random.normal(150.0, 40.0, size=N) * u.km / u.s
-    vt = np.random.normal(100.0, 40.0, size=N)
+    vr = rng.normal(150.0, 40.0, size=N) * u.km / u.s
+    vt = rng.normal(100.0, 40.0, size=N)
     vt = np.vstack((vt, np.zeros_like(vt))).T
 
     # rotate to be random position angle
-    pa = np.random.uniform(0, 2 * np.pi, size=N)
+    pa = rng.uniform(0, 2 * np.pi, size=N)
     M = np.array([[np.cos(pa), -np.sin(pa)], [np.sin(pa), np.cos(pa)]]).T
     vt = np.array([vv.dot(MM) for (vv, MM) in zip(vt, M)]) * u.km / u.s
     vphi, vtheta = vt.T

tests/dynamics/actionangle/test_actionangle_o2gf.py

@@ -184,17 +184,22 @@ def test_actions(self):
 
 
 class TestActions(ActionsBase):
+    @pytest.fixture
+    def rng(self):
+        # The vendored Sanders reference code (_genfunc) crashes on many orbit
+        # samples (most seeds); seed 5 gives 8 orbits it can handle.
+        return np.random.default_rng(5)
+
     @pytest.fixture(autouse=True)
-    def _setup_method(self, tmpdir):
+    def _setup_method(self, tmpdir, rng):
         self.plot_path = tmpdir.mkdir("normal")
 
         self.units = galactic
         self.potential = LeeSutoTriaxialNFWPotential(
             v_c=0.2, r_s=20.0, a=1.0, b=0.77, c=0.55, units=galactic
         )
         self.N = 8
-        np.random.seed(42)
-        w0 = isotropic_w0(N=self.N)
+        w0 = isotropic_w0(rng, N=self.N)
         n_steps = 20000
 
         # integrate orbits
@@ -207,14 +212,14 @@ def _setup_method(self, tmpdir):
         self.w = orbit.w()
 
 
-def test_compare_action_prepare():
+def test_compare_action_prepare(rng):
     from gala.dynamics.actionangle.actionangle_o2gf import (
         _action_prepare,
         _angle_prepare,
     )
 
     logger.setLevel(logging.ERROR)
-    AA = np.random.uniform(0.0, 100.0, size=(1000, 6))
+    AA = rng.uniform(0.0, 100.0, size=(1000, 6))
     t = np.linspace(0.0, 100.0, 1000)
 
     act_san, _n_vectors = solver.solver(AA, N_max=6, symNx=2)

tests/dynamics/actionangle/test_actionangle_staeckel.py

@@ -152,10 +152,9 @@ def test_compare_staeckel_o2gf(find_actions):
 
 
 @requires_gsl
-def test_staeckel_isochrone():
+def test_staeckel_isochrone(rng):
     """In a spherical potential with Delta -> 0, the fudge is exact up to quadrature."""
     pot = gp.IsochronePotential(m=1e11, b=3.0, units=galactic)
-    rng = np.random.default_rng(42)
     N = 256
     w = PhaseSpacePosition(
         pos=rng.uniform(-8, 8, size=(3, N)) * u.kpc,