[MAINT] Use Numpy Random Generator to replace legacy RandomState (#4084)

* Update random generator

* Fix failures

* Revert some changes for sklearn compat

* Update rng in modules where possible

* Fix some tests possibly related to test pollution

* Revert changes to _permuted_ols_on_chunk

* Fix tests

* Remove unused fixture

nilearn/_utils/data_gen.py

@@ -12,7 +12,6 @@
 import scipy.signal
 from nibabel import Nifti1Image
 from scipy.ndimage import binary_dilation
-from sklearn.utils import check_random_state
 
 from nilearn import datasets, image, maskers, masking
 from nilearn._utils import as_ndarray, logger
@@ -47,11 +46,11 @@ def generate_mni_space_img(n_scans=1, res=30, random_state=0, mask_dilation=2):
         Generated mask in MNI space.
 
     """
-    rand_gen = check_random_state(random_state)
+    rand_gen = np.random.default_rng(random_state)
     mask_img = datasets.load_mni152_brain_mask(resolution=res)
     masker = maskers.NiftiMasker(mask_img).fit()
     n_voxels = image.get_data(mask_img).sum()
-    data = rand_gen.randn(n_scans, n_voxels)
+    data = rand_gen.standard_normal((n_scans, n_voxels))
     if mask_dilation is not None and mask_dilation > 0:
         mask_img = image.new_img_like(
             mask_img,
@@ -82,8 +81,8 @@ def generate_timeseries(n_timepoints, n_features, random_state=0):
         Generated time series.
 
     """
-    rand_gen = check_random_state(random_state)
-    return rand_gen.randn(n_timepoints, n_features)
+    rand_gen = np.random.default_rng(random_state)
+    return rand_gen.standard_normal((n_timepoints, n_features))
 
 
 def generate_regions_ts(n_features,
@@ -118,7 +117,7 @@ def generate_regions_ts(n_features,
         shape (n_features, n_regions)
 
     """
-    rand_gen = check_random_state(random_state)
+    rand_gen = np.random.default_rng(random_state)
     if window is None:
         window = "boxcar"
 
@@ -318,9 +317,9 @@ def generate_fake_fmri(shape=(10, 11, 12),
     width = [s // 2 for s in shape]
     shift = [s // 4 for s in shape]
 
-    rand_gen = check_random_state(random_state)
+    rand_gen = np.random.default_rng(random_state)
     if kind == "noise":
-        signals = rand_gen.randint(256, size=(width + [length]))
+        signals = rand_gen.integers(256, size=(width + [length]))
     elif kind == "step":
         signals = np.ones(width + [length])
         signals[..., :length // 2] = 0.5
@@ -348,20 +347,20 @@ def generate_fake_fmri(shape=(10, 11, 12),
                          f'to put {n_blocks} blocks of size {block_size}')
     t_start = 0
     if rest_max_size > 0:
-        t_start = rand_gen.randint(0, rest_max_size, 1)[0]
+        t_start = rand_gen.integers(0, rest_max_size, 1)[0]
     for block in range(n_blocks):
         if block_type == 'classification':
             # Select a random voxel and add some signal to the background
-            voxel_idx = rand_gen.randint(0, flat_fmri.shape[0], 1)[0]
-            trials_effect = (rand_gen.random_sample(block_size) + 1) * 3.
+            voxel_idx = rand_gen.integers(0, flat_fmri.shape[0], 1)[0]
+            trials_effect = (rand_gen.random(block_size) + 1) * 3.
         else:
             # Select the voxel in the image center and add some signal
             # that increases with each block
             voxel_idx = flat_fmri.shape[0] // 2
-            trials_effect = (rand_gen.random_sample(block_size) + 1) * block
+            trials_effect = (rand_gen.random(block_size) + 1) * block
         t_rest = 0
         if rest_max_size > 0:
-            t_rest = rand_gen.randint(0, rest_max_size, 1)[0]
+            t_rest = rand_gen.integers(0, rest_max_size, 1)[0]
         flat_fmri[voxel_idx, t_start:t_start + block_size] += trials_effect
         target[t_start:t_start + block_size] = block + 1
         t_start += t_rest + block_size
@@ -409,17 +408,20 @@ def generate_fake_fmri_data_and_design(shapes,
     """
     fmri_data = []
     design_matrices = []
-    rand_gen = check_random_state(random_state)
+    rand_gen = np.random.default_rng(random_state)
     for i, shape in enumerate(shapes):
-        data = rand_gen.randn(*shape)
+        data = rand_gen.standard_normal(shape)
         data[1:-1, 1:-1, 1:-1] += 100
         fmri_data.append(Nifti1Image(data, affine))
         columns = rand_gen.choice(list(string.ascii_lowercase),
                                   size=rk,
                                   replace=False)
         design_matrices.append(
-            pd.DataFrame(rand_gen.randn(shape[3], rk), columns=columns))
-    mask = Nifti1Image((rand_gen.rand(*shape[:3]) > .5).astype(np.int8),
+            pd.DataFrame(
+                rand_gen.standard_normal((shape[3], rk)), columns=columns
+            )
+        )
+    mask = Nifti1Image((rand_gen.random(shape[:3]) > .5).astype(np.int8),
                        affine)
     return mask, fmri_data, design_matrices
 
@@ -471,19 +473,19 @@ def write_fake_fmri_data_and_design(shapes,
 
     mask_file, fmri_files, design_files = file_path / 'mask.nii', [], []
 
-    rand_gen = check_random_state(random_state)
+    rand_gen = np.random.default_rng(random_state)
     for i, shape in enumerate(shapes):
 
-        data = rand_gen.randn(*shape)
+        data = rand_gen.standard_normal(shape)
         data[1:-1, 1:-1, 1:-1] += 100
         fmri_files.append(str(file_path / f'fmri_run{i:d}.nii'))
         Nifti1Image(data, affine).to_filename(fmri_files[-1])
 
         design_files.append(str(file_path / f'dmtx_{i:d}.csv'))
-        pd.DataFrame(rand_gen.randn(shape[3], rk),
+        pd.DataFrame(rand_gen.standard_normal((shape[3], rk)),
                      columns=['', '', '']).to_csv(design_files[-1])
 
-    Nifti1Image((rand_gen.rand(*shape[:3]) > .5).astype(np.int8),
+    Nifti1Image((rand_gen.random(shape[:3]) > .5).astype(np.int8),
                 affine).to_filename(mask_file)
 
     return mask_file, fmri_files, design_files
@@ -537,8 +539,8 @@ def write_fake_bold_img(file_path,
         Output file path.
 
     """
-    rand_gen = check_random_state(random_state)
-    data = rand_gen.randn(*shape)
+    rand_gen = np.random.default_rng(random_state)
+    data = rand_gen.standard_normal(shape)
     data[1:-1, 1:-1, 1:-1] += 100
     Nifti1Image(data, affine).to_filename(file_path)
     return file_path
@@ -574,12 +576,12 @@ def _generate_signals_from_precisions(precisions,
         (sample number, precisions[n].shape[0]).
 
     """
-    rand_gen = check_random_state(random_state)
+    rand_gen = np.random.default_rng(random_state)
 
     signals = []
-    n_samples = rand_gen.randint(min_n_samples,
-                                 high=max_n_samples,
-                                 size=len(precisions))
+    n_samples = rand_gen.integers(min_n_samples,
+                                  high=max_n_samples,
+                                  size=len(precisions))
 
     mean = np.zeros(precisions[0].shape[0])
     for n, prec in zip(n_samples, precisions):
@@ -635,16 +637,15 @@ def generate_group_sparse_gaussian_graphs(n_subjects=5,
         and signals.
 
     """
-    rand_gen = check_random_state(random_state)
+    rand_gen = np.random.default_rng(random_state)
     # Generate topology (upper triangular binary matrix, with zeros on the
     # diagonal)
     topology = np.empty((n_features, n_features))
     topology[:, :] = np.triu(
-        (rand_gen.randint(0,
-                          high=int(1. / density),
-                          size=n_features * n_features)).reshape(
-                              n_features, n_features) == 0,
-        k=1)
+        (rand_gen.integers(0,
+                           high=int(1. / density),
+                           size=n_features * n_features)).reshape(
+                               n_features, n_features) == 0, k=1)
 
     # Generate edges weights on topology
     precisions = []
@@ -734,11 +735,11 @@ def _basic_confounds(length, random_state=0):
         'trans_x', 'trans_y', 'trans_z'.
 
     """
-    rand_gen = check_random_state(random_state)
+    rand_gen = np.random.default_rng(random_state)
     columns = ['csf', 'white_matter', 'global_signal',
                'rot_x', 'rot_y', 'rot_z',
                'trans_x', 'trans_y', 'trans_z']
-    data = rand_gen.rand(length, len(columns))
+    data = rand_gen.random((length, len(columns)))
     confounds = pd.DataFrame(data, columns=columns)
     return confounds
 
@@ -792,7 +793,7 @@ def add_metadata_to_bids_dataset(bids_path,
 def generate_random_img(
     shape,
     affine=np.eye(4),
-    random_state=np.random.RandomState(0),
+    random_state=0,
 ):
     """Create a random 3D or 4D image with a given shape and affine.
 
@@ -805,8 +806,8 @@ def generate_random_img(
     affine : 4x4 numpy.ndarray
         The affine of the image
 
-    random_state : numpy.random.RandomState instance, optional
-        random number generator.
+    random_state : int, optional
+        Seed for random number generator.
 
     Returns
     -------
@@ -816,7 +817,8 @@ def generate_random_img(
     mask_img : 3D niimg
         The mask image.
     """
-    data = random_state.standard_normal(size=shape)
+    rng = np.random.default_rng(random_state)
+    data = rng.standard_normal(size=shape)
     data_img = Nifti1Image(data, affine)
     if len(shape) == 4:
         mask_data = as_ndarray(data[..., 0] > 0.2, dtype=np.int8)
@@ -910,7 +912,7 @@ def create_fake_bids_dataset(
     """
     n_voxels = 4
 
-    rand_gen = check_random_state(random_state)
+    rand_gen = np.random.default_rng(random_state)
 
     bids_dataset_dir = "bids_dataset"
     bids_path = Path(base_dir) / bids_dataset_dir

nilearn/conftest.py

@@ -132,8 +132,8 @@ def suppress_specific_warning():
 # ------------------------   RNG   ------------------------#
 
 
-def _rng():
-    return np.random.RandomState(42)
+def _rng(seed=42):
+    return np.random.default_rng(seed)
 
 
 @pytest.fixture()
@@ -238,7 +238,7 @@ def _img_3d_rand(affine=_affine_eye()):
 
     Mostly used for set up in other fixtures in other testing modules.
     """
-    data = _rng().rand(*_shape_3d_default())
+    data = _rng().random(_shape_3d_default())
     return Nifti1Image(data, affine)
 
 
@@ -251,7 +251,7 @@ def img_3d_rand_eye():
 def _img_3d_mni(affine=_affine_mni()):
     data_positive = np.zeros((7, 7, 3))
     rng = _rng()
-    data_rng = rng.rand(7, 7, 3)
+    data_rng = rng.random((7, 7, 3))
     data_positive[1:-1, 2:-1, 1:] = data_rng[1:-1, 2:-1, 1:]
     return Nifti1Image(data_positive, affine)
 
@@ -326,7 +326,7 @@ def img_4d_ones_eye():
 @pytest.fixture
 def img_4d_rand_eye():
     """Return a default random filled 4D Nifti1Image (identity affine)."""
-    data = _rng().rand(*_shape_4d_default())
+    data = _rng().random(_shape_4d_default())
     return Nifti1Image(data, _affine_eye())
 
 

nilearn/connectome/tests/test_connectivity_matrices.py

@@ -8,7 +8,6 @@
 from pandas import DataFrame
 from scipy import linalg
 from sklearn.covariance import EmpiricalCovariance, LedoitWolf
-from sklearn.utils import check_random_state
 
 from nilearn._utils.extmath import is_spd
 from nilearn.connectome.connectivity_matrices import (
@@ -60,8 +59,8 @@ def random_diagonal(p, v_min=1.0, v_max=2.0, random_state=0):
         A diagonal matrix with the given minimal and maximal elements.
 
     """
-    random_state = check_random_state(random_state)
-    diag = random_state.rand(p) * (v_max - v_min) + v_min
+    random_state = np.random.default_rng(random_state)
+    diag = random_state.random(p) * (v_max - v_min) + v_min
     diag[diag == np.amax(diag)] = v_max
     diag[diag == np.amin(diag)] = v_min
     return np.diag(diag)
@@ -91,8 +90,8 @@ def random_spd(p, eig_min, cond, random_state=0):
         A symmetric positive definite matrix with the given minimal eigenvalue
         and condition number.
     """
-    random_state = check_random_state(random_state)
-    mat = random_state.randn(p, p)
+    rand_gen = np.random.default_rng(random_state)
+    mat = rand_gen.standard_normal((p, p))
     unitary, _ = linalg.qr(mat)
     diag = random_diagonal(
         p, v_min=eig_min, v_max=cond * eig_min, random_state=random_state
@@ -273,12 +272,12 @@ def random_non_singular(p, sing_min=1.0, sing_max=2.0, random_state=0):
         A nonsingular matrix with the given minimal and maximal singular
         values.
     """
-    random_state = check_random_state(random_state)
+    rand_gen = np.random.default_rng(random_state)
     diag = random_diagonal(
         p, v_min=sing_min, v_max=sing_max, random_state=random_state
     )
-    mat1 = random_state.randn(p, p)
-    mat2 = random_state.randn(p, p)
+    mat1 = rand_gen.standard_normal((p, p))
+    mat2 = rand_gen.standard_normal((p, p))
     unitary1, _ = linalg.qr(mat1)
     unitary2, _ = linalg.qr(mat2)
     return unitary1.dot(diag).dot(unitary2.T)
@@ -468,13 +467,13 @@ def test_sym_matrix_to_vec_is_the_inverse_of_vec_to_sym_matrix(rng):
     p = n * (n + 1) // 2
 
     # when diagonal is included
-    vec = rng.rand(p)
+    vec = rng.random(p)
     sym = vec_to_sym_matrix(vec)
 
     assert_array_almost_equal(sym_matrix_to_vec(sym), vec)
 
     # when diagonal given separately
-    diagonal = rng.rand(n + 1)
+    diagonal = rng.random(n + 1)
     sym = vec_to_sym_matrix(vec, diagonal=diagonal)
 
     assert_array_almost_equal(

nilearn/datasets/tests/test_atlas.py

@@ -334,7 +334,7 @@ def test_fetch_coords_seitzman_2018():
 def _destrieux_data():
     """Function mocking the download of the destrieux atlas."""
     data = {"destrieux2009.rst": "readme"}
-    atlas = _rng().randint(0, 10, (10, 10, 10), dtype="int32")
+    atlas = _rng().integers(0, 10, (10, 10, 10), dtype="int32")
     atlas_img = nibabel.Nifti1Image(atlas, np.eye(4))
     labels = "\n".join([f"{idx},label {idx}" for idx in range(10)])
     labels = f"index,name\n{labels}"

nilearn/datasets/tests/test_func.py

@@ -390,7 +390,9 @@ def test__load_mixed_gambles(rng, affine_eye):
     n_trials = 48
     for n_subjects in [1, 5, 16]:
         zmaps = [
-            nibabel.Nifti1Image(rng.randn(3, 4, 5, n_trials), affine_eye)
+            nibabel.Nifti1Image(
+                rng.standard_normal((3, 4, 5, n_trials)), affine_eye
+            )
             for _ in range(n_subjects)
         ]
         zmaps, gain, _ = func._load_mixed_gambles(zmaps)

nilearn/datasets/tests/test_neurovault.py

@@ -53,7 +53,7 @@ def _get_neurovault_data():
     ).values
     collections["number_of_images"] = collections[
         "true_number_of_images"
-    ] + rng.binomial(1, 0.1, n_collections) * rng.randint(
+    ] + rng.binomial(1, 0.1, n_collections) * rng.integers(
         0, 100, n_collections
     )
     images["not_mni"] = rng.binomial(1, 0.1, size=n_images).astype(bool)
@@ -77,7 +77,7 @@ def _get_neurovault_data():
         p=[0.4, 0.4, 0.2],
     )
     images["some_key"] = "some_value"
-    images[13] = rng.randn(n_images)
+    images[13] = rng.standard_normal(n_images)
     url = "https://neurovault.org/media/images/{}/{}.nii.gz"
     image_names = [
         hashlib.sha1(bytes(img_id)).hexdigest()[:4] for img_id in image_ids

nilearn/decoding/decoder.py

@@ -40,7 +40,6 @@
 )
 from sklearn.preprocessing import LabelBinarizer
 from sklearn.svm import SVR, LinearSVC, l1_min_c
-from sklearn.utils import check_random_state
 from sklearn.utils.extmath import safe_sparse_dot
 from sklearn.utils.validation import check_is_fitted, check_X_y
 
@@ -1051,7 +1050,7 @@ def _predict_dummy(self, n_samples):
             if strategy in ["most_frequent", "prior"]:
                 scores = np.tile(dummy_output, reps=(n_samples, 1))
             elif strategy == "stratified":
-                rs = check_random_state(0)
+                rs = np.random.default_rng(0)
                 scores = rs.multinomial(1, dummy_output, size=n_samples)
 
         elif isinstance(self.estimator, DummyRegressor):