From 4776c72b82660a6c3922783e32ad9576b5e37dd2 Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan <583379+adrn@users.noreply.github.com> Date: Sat, 24 Oct 2026 14:29:28 -0400 Subject: [PATCH 1/4] Add rng argument to fast_lyapunov_max The Cython DOP853 Lyapunov routines drew offset vectors from numpy's global RNG, so callers (and tests) could not seed them. Co-Authored-By: Claude Opus 5.5 (1M context) --- CHANGES.rst | 3 +++ src/gala/dynamics/lyapunov/dop853_lyapunov.pyx | 14 ++++++++++---- src/gala/dynamics/nonlinear.py | 7 +++++++ 3 files changed, 20 insertions(+), 4 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 7a547963a..c1979f5ae 100644 --- a/CHANGES.rst +++ b/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 --------- diff --git a/src/gala/dynamics/lyapunov/dop853_lyapunov.pyx b/src/gala/dynamics/lyapunov/dop853_lyapunov.pyx index c629a4e3e..32544e6eb 100644 --- a/src/gala/dynamics/lyapunov/dop853_lyapunov.pyx +++ b/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,7 +44,7 @@ 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 = ((hamiltonian.potential.c_instance)).cpotential @@ -52,6 +52,9 @@ cpdef dop853_lyapunov_max(hamiltonian, double[::1] w0, 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,7 +126,7 @@ 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 = ((hamiltonian.potential.c_instance)).cpotential @@ -131,6 +134,9 @@ cpdef dop853_lyapunov_max_dont_save(hamiltonian, double[::1] w0, 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 diff --git a/src/gala/dynamics/nonlinear.py b/src/gala/dynamics/nonlinear.py index 721dff918..3fe7132cc 100644 --- a/src/gala/dynamics/nonlinear.py +++ b/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: From e8dcfb58cd2462eed66c55869a08c636ad378a2f Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan <583379+adrn@users.noreply.github.com> Date: Sat, 24 Oct 2026 14:40:13 -0400 Subject: [PATCH 2/4] Seed all random number generation in tests Unseeded draws made tests nondeterministic; e.g. the chaotic forced pendulum Lyapunov test occasionally fell below its threshold in CI. - Replace legacy global np.random.* calls with local np.random.default_rng(42) generators. - Pass seeded rng/random_state to lyapunov_max, fast_lyapunov_max, and mockstream DFs, which otherwise draw from fresh entropy. - isotropic_w0 now takes an explicit RNG; TestActions passes RandomState(42) to keep the sample the Sanders reference is validated on. Co-Authored-By: Claude Opus 5.5 (1M context) --- .../actionangle/actionangle_helpers.py | 14 +- .../actionangle/test_actionangle_o2gf.py | 8 +- tests/dynamics/mockstream/test_df.py | 10 +- tests/dynamics/mockstream/test_mockstream.py | 8 +- .../mockstream/test_mockstream_class.py | 12 +- tests/dynamics/test_dynamics_core.py | 133 ++++++++++-------- tests/dynamics/test_dynamics_util.py | 21 +-- tests/dynamics/test_nonlinear.py | 6 + tests/dynamics/test_orbit.py | 118 +++++++++------- tests/dynamics/test_plot.py | 31 ++-- tests/dynamics/test_representation_nd.py | 10 +- tests/integrate/test_pyintegrators.py | 3 +- tests/potential/frame/test_transformations.py | 41 +++--- .../potential/test_potential_util.py | 4 +- tests/potential/scf/test_bfe.py | 5 +- .../test_psp_orbit_velocity_units.py | 10 +- 16 files changed, 240 insertions(+), 194 deletions(-) diff --git a/tests/dynamics/actionangle/actionangle_helpers.py b/tests/dynamics/actionangle/actionangle_helpers.py index bada61d40..35c8f3aa8 100644 --- a/tests/dynamics/actionangle/actionangle_helpers.py +++ b/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 diff --git a/tests/dynamics/actionangle/test_actionangle_o2gf.py b/tests/dynamics/actionangle/test_actionangle_o2gf.py index 41b1eae0d..971229808 100644 --- a/tests/dynamics/actionangle/test_actionangle_o2gf.py +++ b/tests/dynamics/actionangle/test_actionangle_o2gf.py @@ -193,8 +193,8 @@ def _setup_method(self, tmpdir): 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) + # RandomState reproduces the sample the Sanders reference code is validated on + w0 = isotropic_w0(np.random.RandomState(42), N=self.N) n_steps = 20000 # integrate orbits @@ -213,8 +213,10 @@ def test_compare_action_prepare(): _angle_prepare, ) + rng = np.random.default_rng(42) + 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) diff --git a/tests/dynamics/mockstream/test_df.py b/tests/dynamics/mockstream/test_df.py index 033a15d5d..196c61930 100644 --- a/tests/dynamics/mockstream/test_df.py +++ b/tests/dynamics/mockstream/test_df.py @@ -30,15 +30,15 @@ def test_init_sample(DF, DF_kwargs, pot): n_times = len(orbit.t) # Different ways to initialize successfully: - df = DF(**DF_kwargs) + df = DF(random_state=np.random.default_rng(42), **DF_kwargs) o = df.sample(orbit, 1e4 * u.Msun) assert len(o.x) == 2 * n_times - df = DF(lead=False, **DF_kwargs) + df = DF(lead=False, random_state=np.random.default_rng(42), **DF_kwargs) o = df.sample(orbit, 1e4 * u.Msun) assert len(o.x) == n_times - df = DF(trail=False, **DF_kwargs) + df = DF(trail=False, random_state=np.random.default_rng(42), **DF_kwargs) o = df.sample(orbit, 1e4 * u.Msun) assert len(o.x) == n_times @@ -82,10 +82,10 @@ def test_rotating_frame(): assert u.allclose(o.xyz, orbit_static.xyz, atol=1e-13 * u.kpc) assert u.allclose(o.v_xyz, orbit_static.v_xyz, atol=1e-13 * u.km / u.s) - df_static = DF(trail=False) + df_static = DF(trail=False, random_state=np.random.default_rng(42)) xvt_static = df_static.sample(orbit_static, 1e6 * u.Msun) - df_rotating = DF(trail=False) + df_rotating = DF(trail=False, random_state=np.random.default_rng(42)) xvt_rotating = df_rotating.sample(orbit_rotating, 1e6 * u.Msun) xvt_rotating_static = xvt_rotating.to_frame( H_static.frame, t=xvt_rotating.release_time diff --git a/tests/dynamics/mockstream/test_mockstream.py b/tests/dynamics/mockstream/test_mockstream.py index f3b7f3a16..b0a4c50bf 100644 --- a/tests/dynamics/mockstream/test_mockstream.py +++ b/tests/dynamics/mockstream/test_mockstream.py @@ -142,7 +142,7 @@ def test_run(rng, basic_hamiltonian, progenitor_w0, progenitor_mass): # Test custom n_particles array gen = MockStreamGenerator(df=df, hamiltonian=basic_hamiltonian) - n_particles = np.random.randint(0, 4, size=101) + n_particles = rng.integers(0, 4, size=101) stream3, _ = gen.run( progenitor_w0, progenitor_mass, @@ -677,7 +677,7 @@ def test_leapfrog_progenitor_final_position_uniform_release( pos=[13.0, 0.0, 20.0] * u.kpc, vel=[0, 130.0, 50] * u.km / u.s ) - df = ChenStreamDF() + df = ChenStreamDF(random_state=np.random.default_rng(42)) mw = MilkyWayPotential(version="latest") gen = MockStreamGenerator(df, mw) @@ -723,7 +723,7 @@ def test_leapfrog_progenitor_final_position_sparse_release( pos=[13.0, 0.0, 20.0] * u.kpc, vel=[0, 130.0, 50] * u.km / u.s ) - df = ChenStreamDF() + df = ChenStreamDF(random_state=np.random.default_rng(42)) mw = MilkyWayPotential(version="latest") gen = MockStreamGenerator(df, mw) @@ -765,7 +765,7 @@ def test_leapfrog_vs_dop853_consistency(basic_hamiltonian, progenitor_mass): pos=[13.0, 0.0, 20.0] * u.kpc, vel=[0, 130.0, 50] * u.km / u.s ) - df = ChenStreamDF() + df = ChenStreamDF(random_state=np.random.default_rng(42)) mw = MilkyWayPotential(version="latest") gen = MockStreamGenerator(df, mw) diff --git a/tests/dynamics/mockstream/test_mockstream_class.py b/tests/dynamics/mockstream/test_mockstream_class.py index 1d6a36bef..62429a9da 100644 --- a/tests/dynamics/mockstream/test_mockstream_class.py +++ b/tests/dynamics/mockstream/test_mockstream_class.py @@ -7,9 +7,10 @@ def test_init(): - xyz = np.random.random(size=(3, 100)) * u.kpc - vxyz = np.random.random(size=(3, 100)) * u.km / u.s - t1 = np.random.random(size=100) * u.Myr + rng = np.random.default_rng(42) + xyz = rng.random(size=(3, 100)) * u.kpc + vxyz = rng.random(size=(3, 100)) * u.km / u.s + t1 = rng.random(size=100) * u.Myr lead_trail = np.empty(100, dtype="U1") lead_trail[::2] = "t" @@ -27,8 +28,9 @@ def test_init(): def test_no_copy(): - xyz = np.random.random(size=(3, 100)) * u.kpc - vxyz = np.random.random(size=(3, 100)) * u.km / u.s + rng = np.random.default_rng(42) + xyz = rng.random(size=(3, 100)) * u.kpc + vxyz = rng.random(size=(3, 100)) * u.km / u.s s1 = MockStream(xyz, vxyz, copy=True) s2 = MockStream(xyz, vxyz, copy=False) diff --git a/tests/dynamics/test_dynamics_core.py b/tests/dynamics/test_dynamics_core.py index 66844ed62..44a6705d0 100644 --- a/tests/dynamics/test_dynamics_core.py +++ b/tests/dynamics/test_dynamics_core.py @@ -22,18 +22,19 @@ def test_initialize(): + rng = np.random.default_rng(42) with pytest.raises(ValueError): - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 8)) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 8)) PhaseSpacePosition(pos=x, vel=v) - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) o = PhaseSpacePosition(pos=x, vel=v) assert o.shape == (10,) - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.random(size=(3, 10)) * u.km / u.s + x = rng.random(size=(3, 10)) * u.kpc + v = rng.random(size=(3, 10)) * u.km / u.s o = PhaseSpacePosition(pos=x, vel=v) assert o.xyz.unit == u.kpc assert o.v_x.unit == u.km / u.s @@ -41,8 +42,8 @@ def test_initialize(): assert "s" in o.data.differentials # Not 3D - x = np.random.random(size=(2, 10)) - v = np.random.random(size=(2, 10)) + x = rng.random(size=(2, 10)) + v = rng.random(size=(2, 10)) o = PhaseSpacePosition(pos=x, vel=v) assert o.ndim == 2 @@ -51,14 +52,14 @@ def test_initialize(): assert o.frame is not None assert isinstance(o.frame, StaticFrame) - x = np.random.random(size=(4, 10)) - v = np.random.random(size=(4, 10)) + x = rng.random(size=(4, 10)) + v = rng.random(size=(4, 10)) o = PhaseSpacePosition(pos=x, vel=v) assert o.ndim == 4 # back to 3D - pos = CartesianRepresentation(np.random.random(size=(3, 10)) * u.one) - vel = CartesianDifferential(np.random.random(size=(3, 10)) * u.one) + pos = CartesianRepresentation(rng.random(size=(3, 10)) * u.one) + vel = CartesianDifferential(rng.random(size=(3, 10)) * u.one) o = PhaseSpacePosition(pos=pos, vel=vel) assert hasattr(o, "x") assert hasattr(o, "y") @@ -68,8 +69,8 @@ def test_initialize(): assert hasattr(o, "v_z") # passing a representation with a differential attached - pos = CartesianRepresentation(np.random.random(size=(3, 10)) * u.kpc) - vel = CartesianDifferential(np.random.random(size=(3, 10)) * u.km / u.s) + pos = CartesianRepresentation(rng.random(size=(3, 10)) * u.kpc) + vel = CartesianDifferential(rng.random(size=(3, 10)) * u.km / u.s) o = PhaseSpacePosition(pos.with_differentials({"s": vel})) assert hasattr(o, "x") assert hasattr(o, "y") @@ -91,8 +92,9 @@ def test_initialize(): def test_no_copy(): - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + rng = np.random.default_rng(42) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) o1 = PhaseSpacePosition(pos=x, vel=v, copy=True) o2 = PhaseSpacePosition(pos=x, vel=v, copy=False) assert np.all(o1.w() == o2.w()) @@ -102,7 +104,8 @@ def test_no_copy(): def test_from_w(): - w = np.random.random(size=(6, 10)) + rng = np.random.default_rng(42) + w = rng.random(size=(6, 10)) o = PhaseSpacePosition.from_w(w, galactic) assert o.x.unit == u.kpc assert o.v_x.unit == u.kpc / u.Myr @@ -111,15 +114,16 @@ def test_from_w(): def test_slice(): # simple - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + rng = np.random.default_rng(42) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) o = PhaseSpacePosition(pos=x, vel=v) new_o = o[:5] assert new_o.shape == (5,) # 1d slice on 3d - x = np.random.random(size=(3, 10, 8)) - v = np.random.random(size=(3, 10, 8)) + x = rng.random(size=(3, 10, 8)) + v = rng.random(size=(3, 10, 8)) o = PhaseSpacePosition(pos=x, vel=v) new_o = o[:5] assert new_o.shape == (5, 8) @@ -130,16 +134,16 @@ def test_slice(): assert new_o.shape == (5, 4) # boolean array - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) o = PhaseSpacePosition(pos=x, vel=v) ix = np.array([0, 0, 0, 0, 0, 1, 1, 1, 1, 1]).astype(bool) new_o = o[ix] assert new_o.shape == (sum(ix),) # integer array - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) o = PhaseSpacePosition(pos=x, vel=v) ix = np.array([0, 3, 5]) new_o = o[ix] @@ -147,8 +151,9 @@ def test_slice(): def test_reshape(): - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + rng = np.random.default_rng(42) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) o = PhaseSpacePosition(pos=x, vel=v) new_o = o.reshape((10, 1)) assert new_o.shape == (10, 1) @@ -157,8 +162,9 @@ def test_reshape(): def test_represent_as(): # simple / unitless - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + rng = np.random.default_rng(42) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) o = PhaseSpacePosition(pos=x, vel=v) new_o = o.represent_as(SphericalRepresentation) o.spherical @@ -169,8 +175,8 @@ def test_represent_as(): assert new_o.vel.d_distance.unit == u.one # simple / with units - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s + x = rng.random(size=(3, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = PhaseSpacePosition(pos=x, vel=v) sph = o.represent_as(SphericalRepresentation) assert sph.pos.distance.unit == u.kpc @@ -180,16 +186,17 @@ def test_represent_as(): assert u.allclose(getattr(sph.pos, c), getattr(sph2.pos, c), rtol=1e-12) # doesn't work for 2D - x = np.random.random(size=(2, 10)) - v = np.random.random(size=(2, 10)) + x = rng.random(size=(2, 10)) + v = rng.random(size=(2, 10)) o = PhaseSpacePosition(pos=x, vel=v) with pytest.raises(ValueError): o.represent_as(SphericalRepresentation) def test_represent_as_expected_attributes(): - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + rng = np.random.default_rng(42) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) o = PhaseSpacePosition(pos=x, vel=v) new_o = o.spherical @@ -235,8 +242,8 @@ def test_represent_as_expected_attributes(): assert hasattr(new_o, "v_xyz") # Check that this works with the NDCartesian classes too - x = np.random.random(size=(2, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(2, 10)) * u.km / u.s + x = rng.random(size=(2, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(2, 10)) * u.km / u.s new_o = PhaseSpacePosition(pos=x, vel=v) assert hasattr(new_o, "x1") @@ -249,8 +256,9 @@ def test_represent_as_expected_attributes(): def test_to_coord_frame(): # simple / unitless - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + rng = np.random.default_rng(42) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) o = PhaseSpacePosition(pos=x, vel=v) with ( @@ -260,16 +268,16 @@ def test_to_coord_frame(): o.to_coord_frame(Galactic()) # simple / with units - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s + x = rng.random(size=(3, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = PhaseSpacePosition(pos=x, vel=v) with coord.galactocentric_frame_defaults.set("v4.0"): coo = o.to_coord_frame(Galactic()) assert coo.name == "galactic" # doesn't work for 2D - x = np.random.random(size=(2, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(2, 10)) * u.km / u.s + x = rng.random(size=(2, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(2, 10)) * u.km / u.s o = PhaseSpacePosition(pos=x, vel=v) with coord.galactocentric_frame_defaults.set("v4.0"), pytest.raises(ValueError): o.to_coord_frame(Galactic()) @@ -277,34 +285,35 @@ def test_to_coord_frame(): def test_w(): # simple / unitless - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + rng = np.random.default_rng(42) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) o = PhaseSpacePosition(pos=x, vel=v) w = o.w() assert w.shape == (6, 10) - x = np.random.random(size=3) - v = np.random.random(size=3) + x = rng.random(size=3) + v = rng.random(size=3) o = PhaseSpacePosition(pos=x, vel=v) w = o.w() assert w.shape == (6, 1) # simple / unitless, 2D - x = np.random.random(size=(2, 10)) - v = np.random.random(size=(2, 10)) + x = rng.random(size=(2, 10)) + v = rng.random(size=(2, 10)) o = PhaseSpacePosition(pos=x, vel=v) w = o.w() assert w.shape == (4, 10) - x = np.random.random(size=2) - v = np.random.random(size=2) + x = rng.random(size=2) + v = rng.random(size=2) o = PhaseSpacePosition(pos=x, vel=v) w = o.w() assert w.shape == (4, 1) # simple / with units - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s + x = rng.random(size=(3, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = PhaseSpacePosition(pos=x, vel=v) with pytest.raises(ValueError): o.w() @@ -314,8 +323,8 @@ def test_w(): # simple / with units and potential p = HernquistPotential(units=galactic, m=1e11, c=0.25) - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s + x = rng.random(size=(3, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = PhaseSpacePosition(pos=x, vel=v) w = o.w(p.units) assert np.allclose(x.value, w[:3]) @@ -331,8 +340,9 @@ def test_w(): # ------------------------------------------------------------------------ def test_energy(): # with units - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s + rng = np.random.default_rng(42) + x = rng.random(size=(3, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = PhaseSpacePosition(pos=x, vel=v) KE = o.kinetic_energy() assert KE.unit == (o.v_x.unit) ** 2 @@ -341,14 +351,15 @@ def test_energy(): # with units and potential p = HernquistPotential(units=galactic, m=1e11, c=0.25) H = Hamiltonian(p) - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s + x = rng.random(size=(3, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = PhaseSpacePosition(pos=x, vel=v) PE = o.potential_energy(p) E = o.energy(H) def test_angular_momentum(): + rng = np.random.default_rng(42) w = PhaseSpacePosition([1.0, 0.0, 0.0], [0.0, 0.0, 1.0]) assert u.allclose(np.squeeze(w.angular_momentum()), [0.0, -1, 0] * u.one) @@ -370,8 +381,8 @@ def test_angular_momentum(): assert u.allclose(L, true_L) # multiple - random - q = np.random.uniform(size=(3, 128)) - p = np.random.uniform(size=(3, 128)) + q = rng.uniform(size=(3, 128)) + p = rng.uniform(size=(3, 128)) L = PhaseSpacePosition(q, p).angular_momentum() assert L.shape == (3, 128) diff --git a/tests/dynamics/test_dynamics_util.py b/tests/dynamics/test_dynamics_util.py index b8387b1f0..54017a6e0 100644 --- a/tests/dynamics/test_dynamics_util.py +++ b/tests/dynamics/test_dynamics_util.py @@ -48,21 +48,22 @@ def test_estimate_dt_n_steps(): class TestCombine: def setup_method(self): - x = np.random.random(size=(3,)) - v = np.random.random(size=(3,)) + rng = np.random.default_rng(42) + x = rng.random(size=(3,)) + v = rng.random(size=(3,)) p1 = PhaseSpacePosition(pos=x, vel=v) p2 = PhaseSpacePosition(pos=x, vel=v, frame=StaticFrame(galactic)) - x = np.random.random(size=(3, 5)) - v = np.random.random(size=(3, 5)) + x = rng.random(size=(3, 5)) + v = rng.random(size=(3, 5)) p3 = PhaseSpacePosition(pos=x, vel=v) p4 = PhaseSpacePosition(pos=x * u.kpc, vel=v * u.km / u.s) - x = np.random.random(size=(2, 5)) - v = np.random.random(size=(2, 5)) + x = rng.random(size=(2, 5)) + v = rng.random(size=(2, 5)) p5 = PhaseSpacePosition(pos=x, vel=v) self.psps = [p1, p2, p3, p4, p5] - x = np.random.random(size=(3, 8)) - v = np.random.random(size=(3, 8)) + x = rng.random(size=(3, 8)) + v = rng.random(size=(3, 8)) o1 = Orbit(pos=x, vel=v) o2 = Orbit(pos=x, vel=v, t=np.arange(8)) @@ -75,8 +76,8 @@ def setup_method(self): frame=StaticFrame(galactic), ) - x = np.random.random(size=(2, 8)) - v = np.random.random(size=(2, 8)) + x = rng.random(size=(2, 8)) + v = rng.random(size=(2, 8)) o4 = Orbit(pos=x, vel=v, t=np.arange(8)) self.orbs = [o1, o2, o3, o4] diff --git a/tests/dynamics/test_nonlinear.py b/tests/dynamics/test_nonlinear.py index d7175feac..8f6f0cc2a 100644 --- a/tests/dynamics/test_nonlinear.py +++ b/tests/dynamics/test_nonlinear.py @@ -36,6 +36,7 @@ def test_lyapunov_max(self, tmpdir): d0=d0, n_steps_per_pullback=n_steps_per_pullback, noffset_orbits=noffset, + rng=np.random.default_rng(42), ) regular_LEs = np.mean(regular_LEs, axis=1) @@ -49,6 +50,7 @@ def test_lyapunov_max(self, tmpdir): d0=d0, n_steps_per_pullback=n_steps_per_pullback, noffset_orbits=noffset, + rng=np.random.default_rng(42), ) chaotic_LEs = np.mean(chaotic_LEs, axis=1) assert chaotic_LEs[-1] > 1e-2 @@ -126,6 +128,7 @@ def test_lyapunov_max(self, tmpdir): n_steps=self.n_steps, d0=d0, noffset_orbits=noffset, + rng=np.random.default_rng(42), n_steps_per_pullback=n_steps_per_pullback, ) lyap = np.mean(lyap, axis=1) @@ -218,6 +221,7 @@ def test_fast_lyapunov_max(self, tmpdir): n_steps=self.n_steps, d0=d0, noffset_orbits=noffset, + rng=np.random.default_rng(42), n_steps_per_pullback=n_steps_per_pullback, ) lyap = np.mean(lyap, axis=1) @@ -246,6 +250,7 @@ def F(t, w): n_steps=self.n_steps // 8, d0=d0, noffset_orbits=noffset, + rng=np.random.default_rng(42), n_steps_per_pullback=n_steps_per_pullback, ) lyap1 = np.mean(lyap1, axis=1) @@ -262,6 +267,7 @@ def F(t, w): n_steps=self.n_steps // 8, d0=d0, noffset_orbits=noffset, + rng=np.random.default_rng(42), n_steps_per_pullback=n_steps_per_pullback, units=self.hamiltonian.units, ) diff --git a/tests/dynamics/test_orbit.py b/tests/dynamics/test_orbit.py index 2d369b774..9c76b9ec7 100644 --- a/tests/dynamics/test_orbit.py +++ b/tests/dynamics/test_orbit.py @@ -25,14 +25,15 @@ # Tests below should be cleaned up a bit... def test_initialize(): + rng = np.random.default_rng(42) with pytest.raises(ValueError): - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 8)) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 8)) Orbit(pos=x, vel=v) with pytest.raises(ValueError): - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) t = np.arange(8) Orbit(pos=x, vel=v, t=t) @@ -42,8 +43,8 @@ def test_initialize(): # o = Orbit(pos=x, vel=v) # assert o.ndim == 3 - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.random(size=(3, 10)) * u.km / u.s + x = rng.random(size=(3, 10)) * u.kpc + v = rng.random(size=(3, 10)) * u.km / u.s o = Orbit(pos=x, vel=v) assert o.xyz.unit == u.kpc assert o.v_x.unit == u.km / u.s @@ -56,8 +57,8 @@ def test_initialize(): # assert o.hamiltonian is None # Check that passing in frame and potential or Hamiltonian works - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.random(size=(3, 10)) * u.km / u.s + x = rng.random(size=(3, 10)) * u.kpc + v = rng.random(size=(3, 10)) * u.km / u.s frame = StaticFrame(galactic) potential = LogarithmicPotential( v_c=1.0, r_h=0.14, q1=1.0, q2=0.9, q3=1.0, units=galactic @@ -79,8 +80,9 @@ def test_initialize(): def test_no_copy(): - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + rng = np.random.default_rng(42) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) o1 = Orbit(pos=x, vel=v, copy=True) o2 = Orbit(pos=x, vel=v, copy=False) x[0, 0] = 9999.0 @@ -89,7 +91,8 @@ def test_no_copy(): def test_from_w(): - w = np.random.random(size=(6, 10)) + rng = np.random.default_rng(42) + w = rng.random(size=(6, 10)) o = Orbit.from_w(w, galactic) assert o.xyz.unit == u.kpc assert o.v_x.unit == u.kpc / u.Myr @@ -97,22 +100,23 @@ def test_from_w(): def test_slice(): # simple - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + rng = np.random.default_rng(42) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) o = Orbit(pos=x, vel=v) new_o = o[:5] assert new_o.shape == (5,) - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) t = np.linspace(0, 10, 10) o = Orbit(pos=x, vel=v, t=t) new_o = o[:5] assert new_o.shape == (5,) # 1d slice on 3d - x = np.random.random(size=(3, 10, 8)) - v = np.random.random(size=(3, 10, 8)) + x = rng.random(size=(3, 10, 8)) + v = rng.random(size=(3, 10, 8)) t = np.arange(x.shape[1]) o = Orbit(pos=x, vel=v, t=t) new_o = o[:5] @@ -142,8 +146,8 @@ def test_slice(): assert new_o.t.shape == (5,) # boolean array - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) t = np.arange(x.shape[1]) o = Orbit(pos=x, vel=v, t=t) ix = np.array([0, 0, 0, 0, 0, 1, 1, 1, 1, 1]).astype(bool) @@ -152,8 +156,8 @@ def test_slice(): assert new_o.t.shape == (5,) # boolean array - 3D - x = np.random.random(size=(3, 10, 4)) - v = np.random.random(size=(3, 10, 4)) + x = rng.random(size=(3, 10, 4)) + v = rng.random(size=(3, 10, 4)) t = np.arange(x.shape[1]) o = Orbit(pos=x, vel=v, t=t) ix = np.array([0, 0, 0, 0, 0, 1, 1, 1, 1, 1]).astype(bool) @@ -162,8 +166,8 @@ def test_slice(): assert new_o.t.shape == (5,) # integer array - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) t = np.arange(x.shape[1]) o = Orbit(pos=x, vel=v, t=t) ix = np.array([0, 3, 5]) @@ -174,8 +178,9 @@ def test_slice(): def test_reshape(): # 1d slice on 3d - x = np.random.random(size=(3, 10, 8)) - v = np.random.random(size=(3, 10, 8)) + rng = np.random.default_rng(42) + x = rng.random(size=(3, 10, 8)) + v = rng.random(size=(3, 10, 8)) t = np.arange(x.shape[1]) o = Orbit(pos=x, vel=v, t=t) new_o = o.reshape((10, 4, 2)) @@ -185,8 +190,9 @@ def test_reshape(): def test_represent_as(): # simple / unitless - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + rng = np.random.default_rng(42) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) o = Orbit(pos=x, vel=v) sph = o.represent_as(SphericalRepresentation) @@ -194,8 +200,8 @@ def test_represent_as(): assert sph.vel.d_distance.unit == u.one # simple / with units - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s + x = rng.random(size=(3, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = Orbit(pos=x, vel=v) sph = o.represent_as(SphericalRepresentation) assert sph.pos.distance.unit == u.kpc @@ -203,8 +209,9 @@ def test_represent_as(): def test_represent_as_expected_attributes(): - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s + rng = np.random.default_rng(42) + x = rng.random(size=(3, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = Orbit(pos=x, vel=v) new_o = o.spherical @@ -253,8 +260,8 @@ def test_represent_as_expected_attributes(): assert hasattr(new_o, "v_xyz") # Check that this works with the NDCartesian classes too - x = np.random.random(size=(2, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(2, 10)) * u.km / u.s + x = rng.random(size=(2, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(2, 10)) * u.km / u.s new_o = Orbit(pos=x, vel=v) assert hasattr(new_o, "x1") @@ -267,23 +274,24 @@ def test_represent_as_expected_attributes(): def test_to_coord_frame(): # simple / unitless - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + rng = np.random.default_rng(42) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) o = Orbit(pos=x, vel=v) with pytest.raises(u.UnitConversionError): o.to_coord_frame(Galactic()) # simple / with units - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s + x = rng.random(size=(3, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = Orbit(pos=x, vel=v) coo = o.to_coord_frame(Galactic()) assert coo.name == "galactic" # simple / with units and time - x = np.random.random(size=(3, 128, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(3, 128, 10)) * u.km / u.s + x = rng.random(size=(3, 128, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(3, 128, 10)) * u.km / u.s o = Orbit(pos=x, vel=v) coo = o.to_coord_frame(Galactic()) assert coo.name == "galactic" @@ -291,15 +299,16 @@ def test_to_coord_frame(): def test_w(): # simple / unitless - x = np.random.random(size=(3, 10)) - v = np.random.random(size=(3, 10)) + rng = np.random.default_rng(42) + x = rng.random(size=(3, 10)) + v = rng.random(size=(3, 10)) o = Orbit(pos=x, vel=v) w = o.w() assert w.shape == (6, 10) # simple / with units - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s + x = rng.random(size=(3, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = Orbit(pos=x, vel=v) with pytest.raises(ValueError): o.w() @@ -309,8 +318,8 @@ def test_w(): # simple / with units and potential p = HernquistPotential(units=galactic, m=1e11, c=0.25) - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s + x = rng.random(size=(3, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = Orbit(pos=x, vel=v, potential=p, frame=StaticFrame(galactic)) w = o.w() assert np.allclose(x.value, w[:3, :]) @@ -323,8 +332,9 @@ def test_w(): def test_energy(): # with units - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s + rng = np.random.default_rng(42) + x = rng.random(size=(3, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = Orbit(pos=x, vel=v) KE = o.kinetic_energy() assert KE.unit == (o.v_x.unit) ** 2 @@ -332,8 +342,8 @@ def test_energy(): # with units and potential p = HernquistPotential(units=galactic, m=1e11, c=0.25) - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s + x = rng.random(size=(3, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = Orbit(pos=x, vel=v, potential=p, frame=StaticFrame(galactic)) o.potential_energy() o.energy() @@ -341,8 +351,9 @@ def test_energy(): def test_angular_momentum(): # with units - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s + rng = np.random.default_rng(42) + x = rng.random(size=(3, 10)) * u.kpc + v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = Orbit(pos=x, vel=v) L = o.angular_momentum() assert L.unit == (o.v_x.unit * o.x.unit) @@ -575,13 +586,14 @@ def test_align_circulation(): def test_frame_transform(): + rng = np.random.default_rng(42) static = StaticFrame(galactic) rotating = ConstantRotatingFrame( Omega=[0.53, 1.241, 0.9394] * u.rad / u.Myr, units=galactic ) - x = np.random.random(size=(3, 10)) * u.kpc - v = np.random.random(size=(3, 10)) * u.km / u.s + x = rng.random(size=(3, 10)) * u.kpc + v = rng.random(size=(3, 10)) * u.km / u.s t = np.linspace(0, 1, 10) * u.Myr # no frame specified at init diff --git a/tests/dynamics/test_plot.py b/tests/dynamics/test_plot.py index 7df177462..e377d1453 100644 --- a/tests/dynamics/test_plot.py +++ b/tests/dynamics/test_plot.py @@ -18,6 +18,7 @@ def pytest_generate_tests(metafunc): + rng = np.random.default_rng(42) if "obj" not in metafunc.fixturenames: return @@ -25,40 +26,40 @@ def pytest_generate_tests(metafunc): norbits = 16 object_list.append( - PhaseSpacePosition(pos=np.random.random(size=3), vel=np.random.random(size=3)) + PhaseSpacePosition(pos=rng.random(size=3), vel=rng.random(size=3)) ) object_list.append( PhaseSpacePosition( - pos=np.random.random(size=(3, norbits)), - vel=np.random.random(size=(3, norbits)), + pos=rng.random(size=(3, norbits)), + vel=rng.random(size=(3, norbits)), ) ) object_list.append( PhaseSpacePosition( - pos=np.random.random(size=(3, norbits)) * u.kpc, - vel=np.random.random(size=(3, norbits)) * u.km / u.s, + pos=rng.random(size=(3, norbits)) * u.kpc, + vel=rng.random(size=(3, norbits)) * u.km / u.s, ) ) nsteps = 16 object_list.append( Orbit( - pos=np.random.random(size=(3, nsteps)), - vel=np.random.random(size=(3, nsteps)), + pos=rng.random(size=(3, nsteps)), + vel=rng.random(size=(3, nsteps)), t=np.linspace(0, 1, nsteps), ) ) object_list.append( Orbit( - pos=np.random.random(size=(3, nsteps, 2)), - vel=np.random.random(size=(3, nsteps, 2)), + pos=rng.random(size=(3, nsteps, 2)), + vel=rng.random(size=(3, nsteps, 2)), t=np.linspace(0, 1, nsteps), ) ) object_list.append( Orbit( - pos=np.random.random(size=(3, nsteps)) * u.kpc, - vel=np.random.random(size=(3, nsteps)) * u.km / u.s, + pos=rng.random(size=(3, nsteps)) * u.kpc, + vel=rng.random(size=(3, nsteps)) * u.km / u.s, t=np.linspace(0, 1, nsteps) * u.Myr, ) ) @@ -66,14 +67,14 @@ def pytest_generate_tests(metafunc): # 2D object_list.append( PhaseSpacePosition( - pos=np.random.random(size=(2, norbits)), - vel=np.random.random(size=(2, norbits)), + pos=rng.random(size=(2, norbits)), + vel=rng.random(size=(2, norbits)), ) ) object_list.append( Orbit( - pos=np.random.random(size=(2, nsteps)), - vel=np.random.random(size=(2, nsteps)), + pos=rng.random(size=(2, nsteps)), + vel=rng.random(size=(2, nsteps)), t=np.linspace(0, 1, nsteps), ) ) diff --git a/tests/dynamics/test_representation_nd.py b/tests/dynamics/test_representation_nd.py index 344c83700..19b70b37f 100644 --- a/tests/dynamics/test_representation_nd.py +++ b/tests/dynamics/test_representation_nd.py @@ -10,18 +10,19 @@ def test_init_repr(): # Passing in x1, x2 + rng = np.random.default_rng(42) rep = NDCartesianRepresentation([1.0, 1.0]) assert rep.xyz.shape == (2,) # Passing in x1, x2 - rep = NDCartesianRepresentation(np.random.random(size=(2, 8))) + rep = NDCartesianRepresentation(rng.random(size=(2, 8))) assert rep.xyz.shape == (2, 8) rep[:1] for n in range(1, 6 + 1): print("N: " + str(n)) - xs = np.random.uniform(size=(n, 16)) * u.one + xs = rng.uniform(size=(n, 16)) * u.one rep = NDCartesianRepresentation(xs) for i in range(1, n + 1): assert hasattr(rep, "x" + str(i)) @@ -37,20 +38,21 @@ def test_init_repr(): def test_init_diff(): # Passing in x1, x2 + rng = np.random.default_rng(42) rep = NDCartesianDifferential([1.0, 1.0]) assert rep.d_xyz.shape == (2,) with pytest.raises(TypeError): rep[:1] # Passing in x1, x2 - rep = NDCartesianDifferential(np.random.random(size=(2, 8))) + rep = NDCartesianDifferential(rng.random(size=(2, 8))) assert rep.d_xyz.shape == (2, 8) rep[:1] for n in range(1, 6 + 1): print("N: " + str(n)) - xs = np.random.uniform(size=(n, 16)) * u.one + xs = rng.uniform(size=(n, 16)) * u.one rep = NDCartesianDifferential(xs) for i in range(1, n + 1): assert hasattr(rep, "d_x" + str(i)) diff --git a/tests/integrate/test_pyintegrators.py b/tests/integrate/test_pyintegrators.py index dd0eecb61..5baed192e 100644 --- a/tests/integrate/test_pyintegrators.py +++ b/tests/integrate/test_pyintegrators.py @@ -138,6 +138,7 @@ def test_lorenz(Integrator): @pytest.mark.parametrize("Integrator", integrator_list) def test_memmap(tmpdir, Integrator): + rng = np.random.default_rng(42) dt = 0.1 n_steps = 1000 nw0 = 10000 @@ -145,7 +146,7 @@ def test_memmap(tmpdir, Integrator): filename = os.path.join(str(tmpdir), "test_memmap.npy") mmap = np.memmap(filename, mode="w+", shape=(2, n_steps + 1, nw0)) - w0 = np.random.uniform(-1, 1, size=(2, nw0)) + w0 = rng.uniform(-1, 1, size=(2, nw0)) integrator = Integrator(sho_F, func_args=(1.0,)) diff --git a/tests/potential/frame/test_transformations.py b/tests/potential/frame/test_transformations.py index 875981ec0..da6dfb78e 100644 --- a/tests/potential/frame/test_transformations.py +++ b/tests/potential/frame/test_transformations.py @@ -13,19 +13,20 @@ def test_axis_angle_rotate(): - for x in [np.random.random(size=(3, 32)), np.random.random(size=(3, 32, 8))]: - vec = np.random.random(size=(3, 32)) - theta = np.random.random(size=(32,)) + rng = np.random.default_rng(42) + for x in [rng.random(size=(3, 32)), rng.random(size=(3, 32, 8))]: + vec = rng.random(size=(3, 32)) + theta = rng.random(size=(32,)) out = rodrigues_axis_angle_rotate(x, vec, theta) assert out.shape == x.shape - vec = np.random.random(size=(3,)) - theta = np.random.random(size=(32,)) + vec = rng.random(size=(3,)) + theta = rng.random(size=(32,)) out = rodrigues_axis_angle_rotate(x, vec, theta) assert out.shape == x.shape - vec = np.random.random(size=(3,)) - theta = np.random.random(size=(1,)) + vec = rng.random(size=(3,)) + theta = rng.random(size=(1,)) out = rodrigues_axis_angle_rotate(x, vec, theta) assert out.shape == x.shape @@ -53,28 +54,29 @@ def _helper(fi, fr, w, t=None): def test_frame_transforms_3d(): + rng = np.random.default_rng(42) frame_i = StaticFrame(units=galactic) frame_r = ConstantRotatingFrame( Omega=[0.112, 1.235, 0.8656] * u.rad / u.Myr, units=galactic ) w = Orbit( - pos=np.random.random(size=(3, 32)) * u.kpc, - vel=np.random.random(size=(3, 32)) * u.kpc / u.Myr, + pos=rng.random(size=(3, 32)) * u.kpc, + vel=rng.random(size=(3, 32)) * u.kpc / u.Myr, t=np.linspace(0, 1, 32) * u.Myr, ) _helper(frame_i, frame_r, w, t=w.t) w = Orbit( - pos=np.random.random(size=(3, 32, 8)) * u.kpc, - vel=np.random.random(size=(3, 32, 8)) * u.kpc / u.Myr, + pos=rng.random(size=(3, 32, 8)) * u.kpc, + vel=rng.random(size=(3, 32, 8)) * u.kpc / u.Myr, t=np.linspace(0, 1, 32) * u.Myr, ) _helper(frame_i, frame_r, w, t=w.t) w = PhaseSpacePosition( - pos=np.random.random(size=3) * u.kpc, - vel=np.random.random(size=3) * u.kpc / u.Myr, + pos=rng.random(size=3) * u.kpc, + vel=rng.random(size=3) * u.kpc / u.Myr, ) with pytest.raises(ValueError): _helper(frame_i, frame_r, w) @@ -83,26 +85,27 @@ def test_frame_transforms_3d(): def test_frame_transforms_2d(): + rng = np.random.default_rng(42) frame_i = StaticFrame(units=galactic) frame_r = ConstantRotatingFrame(Omega=0.529 * u.rad / u.Myr, units=galactic) w = Orbit( - pos=np.random.random(size=(2, 32)) * u.kpc, - vel=np.random.random(size=(2, 32)) * u.kpc / u.Myr, + pos=rng.random(size=(2, 32)) * u.kpc, + vel=rng.random(size=(2, 32)) * u.kpc / u.Myr, t=np.linspace(0, 1, 32) * u.Myr, ) _helper(frame_i, frame_r, w, t=w.t) w = Orbit( - pos=np.random.random(size=(2, 32, 8)) * u.kpc, - vel=np.random.random(size=(2, 32, 8)) * u.kpc / u.Myr, + pos=rng.random(size=(2, 32, 8)) * u.kpc, + vel=rng.random(size=(2, 32, 8)) * u.kpc / u.Myr, t=np.linspace(0, 1, 32) * u.Myr, ) _helper(frame_i, frame_r, w, t=w.t) w = PhaseSpacePosition( - pos=np.random.random(size=2) * u.kpc, - vel=np.random.random(size=2) * u.kpc / u.Myr, + pos=rng.random(size=2) * u.kpc, + vel=rng.random(size=2) * u.kpc / u.Myr, ) with pytest.raises(ValueError): _helper(frame_i, frame_r, w) diff --git a/tests/potential/potential/test_potential_util.py b/tests/potential/potential/test_potential_util.py index 6d28a66ef..a59f6586e 100644 --- a/tests/potential/potential/test_potential_util.py +++ b/tests/potential/potential/test_potential_util.py @@ -36,7 +36,9 @@ class TestHarmonicOscillatorFromEquation(EquationBase): def test_derp(self): import numpy as np - self.potential.gradient(np.random.random(size=(1, 13))) + rng = np.random.default_rng(42) + + self.potential.gradient(rng.random(size=(1, 13))) @pytest.mark.skip(reason="to_sympy() not implemented") def test_against_sympy(self): diff --git a/tests/potential/scf/test_bfe.py b/tests/potential/scf/test_bfe.py index 54d8ce34e..c47cdae2c 100644 --- a/tests/potential/scf/test_bfe.py +++ b/tests/potential/scf/test_bfe.py @@ -185,13 +185,14 @@ def Ylmth(l, m, costh): def test_pure_py(): + rng = np.random.default_rng(42) nmax = 6 lmax = 4 # xyz = np.array([[1., 0., 1.], # [1., 1., 0.], # [0., 1., 1.]]) - xyz = np.random.uniform(-2.0, 2.0, size=(128, 3)) + xyz = rng.uniform(-2.0, 2.0, size=(128, 3)) # first try spherical: Snlm = np.zeros((nmax + 1, lmax + 1, lmax + 1)) @@ -209,7 +210,7 @@ def test_pure_py(): assert np.allclose(py_grd, cy_grd) # non-spherical: - Snlm = np.random.uniform(-1, 1, size=(nmax + 1, lmax + 1, lmax + 1)) + Snlm = rng.uniform(-1, 1, size=(nmax + 1, lmax + 1, lmax + 1)) Tnlm = np.zeros_like(Snlm) py_den, py_pot, py_grd = pure_py(xyz.T, Snlm, Tnlm, nmax, lmax) diff --git a/tests/regression/test_psp_orbit_velocity_units.py b/tests/regression/test_psp_orbit_velocity_units.py index 6b0b8293c..5e467eb01 100644 --- a/tests/regression/test_psp_orbit_velocity_units.py +++ b/tests/regression/test_psp_orbit_velocity_units.py @@ -28,8 +28,9 @@ @pytest.mark.parametrize(("pos_unit", "vel_unit", "should_raise"), unit_cases) def test_psp_velocity_units(pos_unit, vel_unit, should_raise): - pos = np.random.random(size=(3, 10)) * pos_unit - vel = np.random.random(size=(3, 10)) * vel_unit + rng = np.random.default_rng(42) + pos = rng.random(size=(3, 10)) * pos_unit + vel = rng.random(size=(3, 10)) * vel_unit if should_raise: with pytest.raises(u.UnitTypeError): @@ -40,8 +41,9 @@ def test_psp_velocity_units(pos_unit, vel_unit, should_raise): @pytest.mark.parametrize(("pos_unit", "vel_unit", "should_raise"), unit_cases) def test_orbit_velocity_units(pos_unit, vel_unit, should_raise): - pos = np.random.random(size=(3, 10)) * pos_unit - vel = np.random.random(size=(3, 10)) * vel_unit + rng = np.random.default_rng(42) + pos = rng.random(size=(3, 10)) * pos_unit + vel = rng.random(size=(3, 10)) * vel_unit t = np.linspace(0, 1, 10) * u.Myr if should_raise: From 7e38a1621bde9676bccbb0e842e0a94fe1b3bd5c Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan <583379+adrn@users.noreply.github.com> Date: Sat, 24 Oct 2026 14:58:36 -0400 Subject: [PATCH 3/4] Use a shared rng fixture; drop np.random.seed and RandomState np.random.seed and RandomState are legacy APIs. All test randomness now comes from np.random.default_rng: - Add a seeded `rng` fixture in conftest.py and use it in test functions instead of constructing generators inline; remove identical local copies. - setup_method/setup_class, helpers, and the Cython _coord test helpers (which can't take fixtures) use a local default_rng(42). This also seeds _test_cyl_to_car_roundtrip, whose seed was commented out. - TestActions overrides rng with seed 5: the vendored Sanders reference code crashes (IndexError / leastsq size errors) on most other samples. Co-Authored-By: Claude Opus 5.5 (1M context) --- conftest.py | 8 ++++ src/gala/dynamics/mockstream/_coord.pyx | 18 ++++---- tests/benchmarks/test_mockstream_benchmark.py | 7 ---- tests/coordinates/test_pm_cov_transform.py | 2 +- .../actionangle/test_actionangle_o2gf.py | 15 ++++--- .../actionangle/test_actionangle_staeckel.py | 3 +- .../actionangle/test_analyticactionangle.py | 6 +-- tests/dynamics/mockstream/test_df.py | 4 +- .../mockstream/test_mockstream_class.py | 6 +-- tests/dynamics/test_dynamics_core.py | 36 ++++++---------- tests/dynamics/test_orbit.py | 36 ++++++---------- tests/dynamics/test_representation_nd.py | 7 +--- tests/integrate/test_pyintegrators.py | 3 +- tests/potential/frame/test_transformations.py | 9 ++-- .../hamiltonian/hamiltonian_helpers.py | 42 +++++++++---------- .../potential/potential/potential_helpers.py | 4 -- tests/potential/potential/test_all_builtin.py | 3 +- .../potential/test_potential_util.py | 6 +-- tests/potential/scf/test_bfe.py | 3 +- tests/potential/scf/test_bfe_interp.py | 3 +- tests/potential/scf/test_computecoeff.py | 9 ++-- .../test_psp_orbit_velocity_units.py | 6 +-- 22 files changed, 96 insertions(+), 140 deletions(-) diff --git a/conftest.py b/conftest.py index 7b4b0db70..f9b390480 100644 --- a/conftest.py +++ b/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) diff --git a/src/gala/dynamics/mockstream/_coord.pyx b/src/gala/dynamics/mockstream/_coord.pyx index d4518a16f..5b8d433d7 100644 --- a/src/gala/dynamics/mockstream/_coord.pyx +++ b/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) diff --git a/tests/benchmarks/test_mockstream_benchmark.py b/tests/benchmarks/test_mockstream_benchmark.py index 37e97aed0..143f92fcb 100644 --- a/tests/benchmarks/test_mockstream_benchmark.py +++ b/tests/benchmarks/test_mockstream_benchmark.py @@ -1,5 +1,4 @@ import astropy.units as u -import numpy as np import pytest import gala.dynamics as gd @@ -7,12 +6,6 @@ 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.""" diff --git a/tests/coordinates/test_pm_cov_transform.py b/tests/coordinates/test_pm_cov_transform.py index cfc307784..2ab4c8630 100644 --- a/tests/coordinates/test_pm_cov_transform.py +++ b/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) diff --git a/tests/dynamics/actionangle/test_actionangle_o2gf.py b/tests/dynamics/actionangle/test_actionangle_o2gf.py index 971229808..d3961ebd4 100644 --- a/tests/dynamics/actionangle/test_actionangle_o2gf.py +++ b/tests/dynamics/actionangle/test_actionangle_o2gf.py @@ -184,8 +184,14 @@ 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 @@ -193,8 +199,7 @@ def _setup_method(self, tmpdir): v_c=0.2, r_s=20.0, a=1.0, b=0.77, c=0.55, units=galactic ) self.N = 8 - # RandomState reproduces the sample the Sanders reference code is validated on - w0 = isotropic_w0(np.random.RandomState(42), N=self.N) + w0 = isotropic_w0(rng, N=self.N) n_steps = 20000 # integrate orbits @@ -207,14 +212,12 @@ 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, ) - rng = np.random.default_rng(42) - logger.setLevel(logging.ERROR) AA = rng.uniform(0.0, 100.0, size=(1000, 6)) t = np.linspace(0.0, 100.0, 1000) diff --git a/tests/dynamics/actionangle/test_actionangle_staeckel.py b/tests/dynamics/actionangle/test_actionangle_staeckel.py index 8353720c4..6b2641f6b 100644 --- a/tests/dynamics/actionangle/test_actionangle_staeckel.py +++ b/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, diff --git a/tests/dynamics/actionangle/test_analyticactionangle.py b/tests/dynamics/actionangle/test_analyticactionangle.py index 847a50999..a356c8441 100644 --- a/tests/dynamics/actionangle/test_analyticactionangle.py +++ b/tests/dynamics/actionangle/test_analyticactionangle.py @@ -122,9 +122,9 @@ class TestHarmonicOscillator: def setup_method(self): logger.info("======== Harmonic Oscillator ========") self.N = 100 - np.random.seed(42) - x = np.random.uniform(-10.0, 10.0, size=(3, self.N)) - v = np.random.uniform(-1.0, 1.0, size=(3, self.N)) / 33.0 + rng = np.random.default_rng(42) + x = rng.uniform(-10.0, 10.0, size=(3, self.N)) + v = rng.uniform(-1.0, 1.0, size=(3, self.N)) / 33.0 w0 = np.vstack((x, v)) self.potential = HarmonicOscillatorPotential( diff --git a/tests/dynamics/mockstream/test_df.py b/tests/dynamics/mockstream/test_df.py index 196c61930..48549254f 100644 --- a/tests/dynamics/mockstream/test_df.py +++ b/tests/dynamics/mockstream/test_df.py @@ -42,9 +42,9 @@ def test_init_sample(DF, DF_kwargs, pot): o = df.sample(orbit, 1e4 * u.Msun) assert len(o.x) == n_times - df1 = DF(random_state=np.random.RandomState(42), **DF_kwargs) + df1 = DF(random_state=np.random.default_rng(42), **DF_kwargs) o1 = df1.sample(orbit, 1e4 * u.Msun) - df2 = DF(random_state=np.random.RandomState(42), **DF_kwargs) + df2 = DF(random_state=np.random.default_rng(42), **DF_kwargs) o2 = df2.sample(orbit, 1e4 * u.Msun) assert u.allclose(o1.xyz, o2.xyz) assert u.allclose(o1.v_xyz, o2.v_xyz) diff --git a/tests/dynamics/mockstream/test_mockstream_class.py b/tests/dynamics/mockstream/test_mockstream_class.py index 62429a9da..d28e4862d 100644 --- a/tests/dynamics/mockstream/test_mockstream_class.py +++ b/tests/dynamics/mockstream/test_mockstream_class.py @@ -6,8 +6,7 @@ from gala.dynamics.mockstream import MockStream -def test_init(): - rng = np.random.default_rng(42) +def test_init(rng): xyz = rng.random(size=(3, 100)) * u.kpc vxyz = rng.random(size=(3, 100)) * u.km / u.s t1 = rng.random(size=100) * u.Myr @@ -27,8 +26,7 @@ def test_init(): MockStream(xyz, vxyz, lead_trail=lead_trail[:-1]) -def test_no_copy(): - rng = np.random.default_rng(42) +def test_no_copy(rng): xyz = rng.random(size=(3, 100)) * u.kpc vxyz = rng.random(size=(3, 100)) * u.km / u.s diff --git a/tests/dynamics/test_dynamics_core.py b/tests/dynamics/test_dynamics_core.py index 44a6705d0..d202c3baa 100644 --- a/tests/dynamics/test_dynamics_core.py +++ b/tests/dynamics/test_dynamics_core.py @@ -21,8 +21,7 @@ from gala.units import galactic, solarsystem -def test_initialize(): - rng = np.random.default_rng(42) +def test_initialize(rng): with pytest.raises(ValueError): x = rng.random(size=(3, 10)) v = rng.random(size=(3, 8)) @@ -91,8 +90,7 @@ def test_initialize(): o = PhaseSpacePosition(pos=x, vel=v, frame="blah blah blah") -def test_no_copy(): - rng = np.random.default_rng(42) +def test_no_copy(rng): x = rng.random(size=(3, 10)) v = rng.random(size=(3, 10)) o1 = PhaseSpacePosition(pos=x, vel=v, copy=True) @@ -103,8 +101,7 @@ def test_no_copy(): assert o2.x[0].value == 9999.0 -def test_from_w(): - rng = np.random.default_rng(42) +def test_from_w(rng): w = rng.random(size=(6, 10)) o = PhaseSpacePosition.from_w(w, galactic) assert o.x.unit == u.kpc @@ -112,9 +109,8 @@ def test_from_w(): assert o.shape == (10,) -def test_slice(): +def test_slice(rng): # simple - rng = np.random.default_rng(42) x = rng.random(size=(3, 10)) v = rng.random(size=(3, 10)) o = PhaseSpacePosition(pos=x, vel=v) @@ -150,8 +146,7 @@ def test_slice(): assert new_o.shape == (len(ix),) -def test_reshape(): - rng = np.random.default_rng(42) +def test_reshape(rng): x = rng.random(size=(3, 10)) v = rng.random(size=(3, 10)) o = PhaseSpacePosition(pos=x, vel=v) @@ -160,9 +155,8 @@ def test_reshape(): assert new_o.x.shape == (10, 1) -def test_represent_as(): +def test_represent_as(rng): # simple / unitless - rng = np.random.default_rng(42) x = rng.random(size=(3, 10)) v = rng.random(size=(3, 10)) o = PhaseSpacePosition(pos=x, vel=v) @@ -193,8 +187,7 @@ def test_represent_as(): o.represent_as(SphericalRepresentation) -def test_represent_as_expected_attributes(): - rng = np.random.default_rng(42) +def test_represent_as_expected_attributes(rng): x = rng.random(size=(3, 10)) v = rng.random(size=(3, 10)) o = PhaseSpacePosition(pos=x, vel=v) @@ -254,9 +247,8 @@ def test_represent_as_expected_attributes(): assert hasattr(new_o, "v_xyz") -def test_to_coord_frame(): +def test_to_coord_frame(rng): # simple / unitless - rng = np.random.default_rng(42) x = rng.random(size=(3, 10)) v = rng.random(size=(3, 10)) o = PhaseSpacePosition(pos=x, vel=v) @@ -283,9 +275,8 @@ def test_to_coord_frame(): o.to_coord_frame(Galactic()) -def test_w(): +def test_w(rng): # simple / unitless - rng = np.random.default_rng(42) x = rng.random(size=(3, 10)) v = rng.random(size=(3, 10)) o = PhaseSpacePosition(pos=x, vel=v) @@ -338,9 +329,8 @@ def test_w(): # ------------------------------------------------------------------------ # Computed dynamical quantities # ------------------------------------------------------------------------ -def test_energy(): +def test_energy(rng): # with units - rng = np.random.default_rng(42) x = rng.random(size=(3, 10)) * u.kpc v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = PhaseSpacePosition(pos=x, vel=v) @@ -358,8 +348,7 @@ def test_energy(): E = o.energy(H) -def test_angular_momentum(): - rng = np.random.default_rng(42) +def test_angular_momentum(rng): w = PhaseSpacePosition([1.0, 0.0, 0.0], [0.0, 0.0, 1.0]) assert u.allclose(np.squeeze(w.angular_momentum()), [0.0, -1, 0] * u.one) @@ -387,8 +376,7 @@ def test_angular_momentum(): assert L.shape == (3, 128) -def test_guiding_radius(): - rng = np.random.default_rng(42) +def test_guiding_radius(rng): p = HernquistPotential(units=galactic, m=1e11, c=10.0) diff --git a/tests/dynamics/test_orbit.py b/tests/dynamics/test_orbit.py index 9c76b9ec7..a628cee34 100644 --- a/tests/dynamics/test_orbit.py +++ b/tests/dynamics/test_orbit.py @@ -24,8 +24,7 @@ # Tests below should be cleaned up a bit... -def test_initialize(): - rng = np.random.default_rng(42) +def test_initialize(rng): with pytest.raises(ValueError): x = rng.random(size=(3, 10)) v = rng.random(size=(3, 8)) @@ -79,8 +78,7 @@ def test_initialize(): assert isinstance(o.frame, StaticFrame) -def test_no_copy(): - rng = np.random.default_rng(42) +def test_no_copy(rng): x = rng.random(size=(3, 10)) v = rng.random(size=(3, 10)) o1 = Orbit(pos=x, vel=v, copy=True) @@ -90,17 +88,15 @@ def test_no_copy(): assert o2.x[0].value == 9999.0 -def test_from_w(): - rng = np.random.default_rng(42) +def test_from_w(rng): w = rng.random(size=(6, 10)) o = Orbit.from_w(w, galactic) assert o.xyz.unit == u.kpc assert o.v_x.unit == u.kpc / u.Myr -def test_slice(): +def test_slice(rng): # simple - rng = np.random.default_rng(42) x = rng.random(size=(3, 10)) v = rng.random(size=(3, 10)) o = Orbit(pos=x, vel=v) @@ -176,9 +172,8 @@ def test_slice(): assert new_o.t.shape == (len(ix),) -def test_reshape(): +def test_reshape(rng): # 1d slice on 3d - rng = np.random.default_rng(42) x = rng.random(size=(3, 10, 8)) v = rng.random(size=(3, 10, 8)) t = np.arange(x.shape[1]) @@ -188,9 +183,8 @@ def test_reshape(): assert new_o.x.shape == (10, 4, 2) -def test_represent_as(): +def test_represent_as(rng): # simple / unitless - rng = np.random.default_rng(42) x = rng.random(size=(3, 10)) v = rng.random(size=(3, 10)) o = Orbit(pos=x, vel=v) @@ -208,8 +202,7 @@ def test_represent_as(): assert sph.vel.d_distance.unit == u.km / u.s -def test_represent_as_expected_attributes(): - rng = np.random.default_rng(42) +def test_represent_as_expected_attributes(rng): x = rng.random(size=(3, 10)) * u.kpc v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = Orbit(pos=x, vel=v) @@ -272,9 +265,8 @@ def test_represent_as_expected_attributes(): assert hasattr(new_o, "v_xyz") -def test_to_coord_frame(): +def test_to_coord_frame(rng): # simple / unitless - rng = np.random.default_rng(42) x = rng.random(size=(3, 10)) v = rng.random(size=(3, 10)) o = Orbit(pos=x, vel=v) @@ -297,9 +289,8 @@ def test_to_coord_frame(): assert coo.name == "galactic" -def test_w(): +def test_w(rng): # simple / unitless - rng = np.random.default_rng(42) x = rng.random(size=(3, 10)) v = rng.random(size=(3, 10)) o = Orbit(pos=x, vel=v) @@ -330,9 +321,8 @@ def test_w(): assert np.allclose(v.value, (w[3:, :] * u.au / u.yr).to(u.km / u.s).value) -def test_energy(): +def test_energy(rng): # with units - rng = np.random.default_rng(42) x = rng.random(size=(3, 10)) * u.kpc v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = Orbit(pos=x, vel=v) @@ -349,9 +339,8 @@ def test_energy(): o.energy() -def test_angular_momentum(): +def test_angular_momentum(rng): # with units - rng = np.random.default_rng(42) x = rng.random(size=(3, 10)) * u.kpc v = rng.normal(0.0, 100.0, size=(3, 10)) * u.km / u.s o = Orbit(pos=x, vel=v) @@ -585,8 +574,7 @@ def test_align_circulation(): assert np.all(new_circ[:, 3] == 0.0) -def test_frame_transform(): - rng = np.random.default_rng(42) +def test_frame_transform(rng): static = StaticFrame(galactic) rotating = ConstantRotatingFrame( Omega=[0.53, 1.241, 0.9394] * u.rad / u.Myr, units=galactic diff --git a/tests/dynamics/test_representation_nd.py b/tests/dynamics/test_representation_nd.py index 19b70b37f..d3ccfb82f 100644 --- a/tests/dynamics/test_representation_nd.py +++ b/tests/dynamics/test_representation_nd.py @@ -1,5 +1,4 @@ import astropy.units as u -import numpy as np import pytest from gala.dynamics.representation_nd import ( @@ -8,9 +7,8 @@ ) -def test_init_repr(): +def test_init_repr(rng): # Passing in x1, x2 - rng = np.random.default_rng(42) rep = NDCartesianRepresentation([1.0, 1.0]) assert rep.xyz.shape == (2,) @@ -36,9 +34,8 @@ def test_init_repr(): assert rep2.shape == (8,) -def test_init_diff(): +def test_init_diff(rng): # Passing in x1, x2 - rng = np.random.default_rng(42) rep = NDCartesianDifferential([1.0, 1.0]) assert rep.d_xyz.shape == (2,) with pytest.raises(TypeError): diff --git a/tests/integrate/test_pyintegrators.py b/tests/integrate/test_pyintegrators.py index 5baed192e..a1ebb5099 100644 --- a/tests/integrate/test_pyintegrators.py +++ b/tests/integrate/test_pyintegrators.py @@ -137,8 +137,7 @@ def test_lorenz(Integrator): @pytest.mark.parametrize("Integrator", integrator_list) -def test_memmap(tmpdir, Integrator): - rng = np.random.default_rng(42) +def test_memmap(tmpdir, Integrator, rng): dt = 0.1 n_steps = 1000 nw0 = 10000 diff --git a/tests/potential/frame/test_transformations.py b/tests/potential/frame/test_transformations.py index da6dfb78e..8bfd90cfa 100644 --- a/tests/potential/frame/test_transformations.py +++ b/tests/potential/frame/test_transformations.py @@ -12,8 +12,7 @@ from gala.units import galactic -def test_axis_angle_rotate(): - rng = np.random.default_rng(42) +def test_axis_angle_rotate(rng): for x in [rng.random(size=(3, 32)), rng.random(size=(3, 32, 8))]: vec = rng.random(size=(3, 32)) theta = rng.random(size=(32,)) @@ -53,8 +52,7 @@ def _helper(fi, fr, w, t=None): assert u.allclose(vel_r, w.v_xyz) -def test_frame_transforms_3d(): - rng = np.random.default_rng(42) +def test_frame_transforms_3d(rng): frame_i = StaticFrame(units=galactic) frame_r = ConstantRotatingFrame( Omega=[0.112, 1.235, 0.8656] * u.rad / u.Myr, units=galactic @@ -84,8 +82,7 @@ def test_frame_transforms_3d(): _helper(frame_i, frame_r, w, t=0.0) -def test_frame_transforms_2d(): - rng = np.random.default_rng(42) +def test_frame_transforms_2d(rng): frame_i = StaticFrame(units=galactic) frame_r = ConstantRotatingFrame(Omega=0.529 * u.rad / u.Myr, units=galactic) diff --git a/tests/potential/hamiltonian/hamiltonian_helpers.py b/tests/potential/hamiltonian/hamiltonian_helpers.py index e802d2c51..e6aa7bb9f 100644 --- a/tests/potential/hamiltonian/hamiltonian_helpers.py +++ b/tests/potential/hamiltonian/hamiltonian_helpers.py @@ -14,7 +14,7 @@ class _TestBase: @classmethod def setup_class(cls): - np.random.seed(42) + rng = np.random.default_rng(42) ndim = 6 r_ndim = ndim # return ndim @@ -32,14 +32,14 @@ def setup_class(cls): # 1D - phase-space position cls.w0s.append( PSP( - pos=np.random.random(size=ndim // 2), - vel=np.random.random(size=ndim // 2), + pos=rng.random(size=ndim // 2), + vel=rng.random(size=ndim // 2), ) ) cls.w0s.append( PSP( - pos=np.random.random(size=ndim // 2) * u.kpc, - vel=np.random.random(size=ndim // 2) * u.km / u.s, + pos=rng.random(size=ndim // 2) * u.kpc, + vel=rng.random(size=ndim // 2) * u.km / u.s, ) ) cls.energy_return_shapes += [(1,)] * 2 @@ -49,14 +49,14 @@ def setup_class(cls): # 2D - phase-space position cls.w0s.append( PSP( - pos=np.random.random(size=(ndim // 2, norbits)), - vel=np.random.random(size=(ndim // 2, norbits)), + pos=rng.random(size=(ndim // 2, norbits)), + vel=rng.random(size=(ndim // 2, norbits)), ) ) cls.w0s.append( PSP( - pos=np.random.random(size=(ndim // 2, norbits)) * u.kpc, - vel=np.random.random(size=(ndim // 2, norbits)) * u.km / u.s, + pos=rng.random(size=(ndim // 2, norbits)) * u.kpc, + vel=rng.random(size=(ndim // 2, norbits)) * u.km / u.s, ) ) cls.energy_return_shapes += [(norbits,)] * 2 @@ -66,14 +66,14 @@ def setup_class(cls): # 3D - phase-space position cls.w0s.append( PSP( - pos=np.random.random(size=(ndim // 2, norbits, ntimes)), - vel=np.random.random(size=(ndim // 2, norbits, ntimes)), + pos=rng.random(size=(ndim // 2, norbits, ntimes)), + vel=rng.random(size=(ndim // 2, norbits, ntimes)), ) ) cls.w0s.append( PSP( - pos=np.random.random(size=(ndim // 2, norbits, ntimes)) * u.kpc, - vel=np.random.random(size=(ndim // 2, norbits, ntimes)) * u.km / u.s, + pos=rng.random(size=(ndim // 2, norbits, ntimes)) * u.kpc, + vel=rng.random(size=(ndim // 2, norbits, ntimes)) * u.km / u.s, ) ) cls.energy_return_shapes += [(norbits, ntimes)] * 2 @@ -83,14 +83,14 @@ def setup_class(cls): # 2D - orbit cls.w0s.append( ORB( - pos=np.random.random(size=(ndim // 2, ntimes)), - vel=np.random.random(size=(ndim // 2, ntimes)), + pos=rng.random(size=(ndim // 2, ntimes)), + vel=rng.random(size=(ndim // 2, ntimes)), ) ) cls.w0s.append( ORB( - pos=np.random.random(size=(ndim // 2, ntimes)) * u.kpc, - vel=np.random.random(size=(ndim // 2, ntimes)) * u.km / u.s, + pos=rng.random(size=(ndim // 2, ntimes)) * u.kpc, + vel=rng.random(size=(ndim // 2, ntimes)) * u.km / u.s, ) ) cls.energy_return_shapes += [(ntimes,)] * 2 @@ -111,14 +111,14 @@ def setup_class(cls): # 3D - orbit cls.w0s.append( ORB( - pos=np.random.random(size=(ndim // 2, ntimes, norbits)), - vel=np.random.random(size=(ndim // 2, ntimes, norbits)), + pos=rng.random(size=(ndim // 2, ntimes, norbits)), + vel=rng.random(size=(ndim // 2, ntimes, norbits)), ) ) cls.w0s.append( ORB( - pos=np.random.random(size=(ndim // 2, ntimes, norbits)) * u.kpc, - vel=np.random.random(size=(ndim // 2, ntimes, norbits)) * u.km / u.s, + pos=rng.random(size=(ndim // 2, ntimes, norbits)) * u.kpc, + vel=rng.random(size=(ndim // 2, ntimes, norbits)) * u.km / u.s, ) ) cls.energy_return_shapes += [(ntimes, norbits)] * 2 diff --git a/tests/potential/potential/potential_helpers.py b/tests/potential/potential/potential_helpers.py index 48c4f098b..67c6cb539 100644 --- a/tests/potential/potential/potential_helpers.py +++ b/tests/potential/potential/potential_helpers.py @@ -36,10 +36,6 @@ class PotentialTestBase: num_dx = None num_max_x = None - @pytest.fixture - def rng(self): - return np.random.default_rng(42) - def setup_method(self): # set up hamiltonian if self.frame is None: diff --git a/tests/potential/potential/test_all_builtin.py b/tests/potential/potential/test_all_builtin.py index 62323e792..9806c8b7b 100644 --- a/tests/potential/potential/test_all_builtin.py +++ b/tests/potential/potential/test_all_builtin.py @@ -649,7 +649,7 @@ def test_cylspline_against_agama_export(name): @pytest.mark.skipif(not GSL_ENABLED, reason="requires GSL to run this test") @pytest.mark.parametrize("name", ["mn", "triaxial", "tilted"]) -def test_cylspline_hessian_vs_gradient(name): +def test_cylspline_hessian_vs_gradient(name, rng): """ The Hessian should match finite differences of the gradient, including near and on the z axis and outside of the grid. @@ -657,7 +657,6 @@ def test_cylspline_hessian_vs_gradient(name): pot = p.CylSplinePotential.from_file( this_path / f"agama_cylspline_{name}.pot", units=galactic ) - rng = np.random.default_rng(42) xyz = np.hstack( ( rng.uniform(-20, 20, size=(3, 16)) * [[1], [1], [0.3]], diff --git a/tests/potential/potential/test_potential_util.py b/tests/potential/potential/test_potential_util.py index a59f6586e..616cc2a57 100644 --- a/tests/potential/potential/test_potential_util.py +++ b/tests/potential/potential/test_potential_util.py @@ -33,11 +33,7 @@ class TestHarmonicOscillatorFromEquation(EquationBase): potential = Potential(k=1.0) w0 = [1.0, 0.0] - def test_derp(self): - import numpy as np - - rng = np.random.default_rng(42) - + def test_derp(self, rng): self.potential.gradient(rng.random(size=(1, 13))) @pytest.mark.skip(reason="to_sympy() not implemented") diff --git a/tests/potential/scf/test_bfe.py b/tests/potential/scf/test_bfe.py index c47cdae2c..e482b1f6d 100644 --- a/tests/potential/scf/test_bfe.py +++ b/tests/potential/scf/test_bfe.py @@ -184,8 +184,7 @@ def Ylmth(l, m, costh): return density, potenti, gradien -def test_pure_py(): - rng = np.random.default_rng(42) +def test_pure_py(rng): nmax = 6 lmax = 4 diff --git a/tests/potential/scf/test_bfe_interp.py b/tests/potential/scf/test_bfe_interp.py index 12444ea5b..ba21be303 100644 --- a/tests/potential/scf/test_bfe_interp.py +++ b/tests/potential/scf/test_bfe_interp.py @@ -11,11 +11,10 @@ @pytest.mark.parametrize("func_name", ["energy", "density", "gradient"]) -def test_simple_compare_noninterp(func_name): +def test_simple_compare_noninterp(func_name, rng): """ Compare the interpolated to time-invariant versions for a trivial case """ - rng = np.random.default_rng(42) nmax = 5 lmax = 3 diff --git a/tests/potential/scf/test_computecoeff.py b/tests/potential/scf/test_computecoeff.py index 35de558fe..d677d3391 100644 --- a/tests/potential/scf/test_computecoeff.py +++ b/tests/potential/scf/test_computecoeff.py @@ -154,7 +154,7 @@ def flattened_hernquist_gradient(x, y, z, G, M, a, q): return np.array([gx, gy, gz]) -def test_flattened_hernquist(): +def test_flattened_hernquist(rng): """ This test compares the coefficients against some computed in the mathematica notebook 'flattened-hernquist.nb'. nmax and lmax here must match nmax and lmax @@ -193,11 +193,10 @@ def test_flattened_hernquist(): assert np.allclose(Snlm[0, ::2, 0], m_Snl0[0]) # check that random points match in gradient and density - np.random.seed(42) n_test = 1024 - r = 10.0 * np.cbrt(np.random.uniform(0.1**3, 1, size=n_test)) # 1 to 10 - t = np.arccos(2 * np.random.uniform(size=n_test) - 1) - ph = np.random.uniform(0, 2 * np.pi, size=n_test) + r = 10.0 * np.cbrt(rng.uniform(0.1**3, 1, size=n_test)) # 1 to 10 + t = np.arccos(2 * rng.uniform(size=n_test) - 1) + ph = rng.uniform(0, 2 * np.pi, size=n_test) x = r * np.cos(ph) * np.sin(t) y = r * np.sin(ph) * np.sin(t) z = r * np.cos(t) diff --git a/tests/regression/test_psp_orbit_velocity_units.py b/tests/regression/test_psp_orbit_velocity_units.py index 5e467eb01..b9231888f 100644 --- a/tests/regression/test_psp_orbit_velocity_units.py +++ b/tests/regression/test_psp_orbit_velocity_units.py @@ -27,8 +27,7 @@ @pytest.mark.parametrize(("pos_unit", "vel_unit", "should_raise"), unit_cases) -def test_psp_velocity_units(pos_unit, vel_unit, should_raise): - rng = np.random.default_rng(42) +def test_psp_velocity_units(pos_unit, vel_unit, should_raise, rng): pos = rng.random(size=(3, 10)) * pos_unit vel = rng.random(size=(3, 10)) * vel_unit @@ -40,8 +39,7 @@ def test_psp_velocity_units(pos_unit, vel_unit, should_raise): @pytest.mark.parametrize(("pos_unit", "vel_unit", "should_raise"), unit_cases) -def test_orbit_velocity_units(pos_unit, vel_unit, should_raise): - rng = np.random.default_rng(42) +def test_orbit_velocity_units(pos_unit, vel_unit, should_raise, rng): pos = rng.random(size=(3, 10)) * pos_unit vel = rng.random(size=(3, 10)) * vel_unit t = np.linspace(0, 1, 10) * u.Myr From 7eedbf59df3ca86d0055dd5ff17a3965d4659b8b Mon Sep 17 00:00:00 2001 From: Adrian Price-Whelan <583379+adrn@users.noreply.github.com> Date: Sat, 24 Oct 2026 15:29:07 -0400 Subject: [PATCH 4/4] Use numpy Generators in mockstream DFs; deprecate RandomState random_state now goes through np.random.default_rng, so it defaults to a Generator and also accepts an int seed. A legacy RandomState is still used as-is but raises GalaDeprecationWarning. Co-Authored-By: Claude Opus 5.5 (1M context) --- CHANGES.rst | 7 +++++ src/gala/dynamics/mockstream/df.pyx | 38 +++++++++++++++++++--------- tests/dynamics/mockstream/test_df.py | 17 +++++++++++++ 3 files changed, 50 insertions(+), 12 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index c1979f5ae..603b44cf1 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -91,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 ----- diff --git a/src/gala/dynamics/mockstream/df.pyx b/src/gala/dynamics/mockstream/df.pyx index b34759fb5..e4dea087e 100644 --- a/src/gala/dynamics/mockstream/df.pyx +++ b/src/gala/dynamics/mockstream/df.pyx @@ -41,8 +41,9 @@ 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): @@ -50,8 +51,17 @@ cdef class BaseStreamDF: 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 diff --git a/tests/dynamics/mockstream/test_df.py b/tests/dynamics/mockstream/test_df.py index 48549254f..774460039 100644 --- a/tests/dynamics/mockstream/test_df.py +++ b/tests/dynamics/mockstream/test_df.py @@ -58,6 +58,23 @@ def test_expected_failure(DF, DF_kwargs): DF(lead=False, trail=False, **DF_kwargs) +def test_random_state(): + orbit = gp.Hamiltonian(_TEST_POTENTIALS[0]).integrate_orbit( + [10.0, 0, 0, 0, 0.2, 0], dt=1.0, n_steps=100 + ) + + # an integer seed is passed through np.random.default_rng + df1 = ms.FardalStreamDF(gala_modified=True, random_state=42) + df2 = ms.FardalStreamDF(gala_modified=True, random_state=42) + assert isinstance(df1.random_state, np.random.Generator) + o1 = df1.sample(orbit, 1e4 * u.Msun) + o2 = df2.sample(orbit, 1e4 * u.Msun) + assert u.allclose(o1.xyz, o2.xyz) + + with pytest.warns(DeprecationWarning, match="RandomState"): + ms.FardalStreamDF(gala_modified=True, random_state=np.random.RandomState(42)) + + def test_rotating_frame(): DF = _DF_CLASSES[0] H_static = gp.Hamiltonian(_TEST_POTENTIALS[0])