[WIP] Ensure that random number generation in tests is seeded & stand…

…ardized (#2574)

* Standardize random number generation for tests across glm, image, plotting, regions & tests (#2534)

* Standardize random number generation for tests across input_data, decoding + flake8 (#2534)

* Instantiate RandomState outside loops & list-comps (#2574)

nilearn/decoding/tests/test_searchlight.py

@@ -80,8 +80,10 @@ def test_searchlight():
         # the groups variable should have no effect.
         gcv = cv
 
-    groups = np.random.permutation(np.arange(frames, dtype=int) >
-                                   (frames // 2))
+    groups = np.random.RandomState(42).permutation(
+        np.arange(frames, dtype=int) > (frames // 2)
+    )
+
     sl = searchlight.SearchLight(mask_img, process_mask_img=mask_img, radius=1,
                                  n_jobs=n_jobs, scoring='accuracy', cv=gcv)
     sl.fit(data_img, cond, groups)

nilearn/glm/tests/test_contrasts.py

@@ -27,8 +27,9 @@ def test_expression_to_contrast_vector():
 
 
 def test_Tcontrast():
+    rng = np.random.RandomState(42)
     n, p, q = 100, 80, 10
-    X, Y = np.random.randn(p, q), np.random.randn(p, n)
+    X, Y = rng.standard_normal(size=(p, q)), rng.standard_normal(size=(p, n))
     labels, results = run_glm(Y, X, 'ar1')
     con_val = np.eye(q)[0]
     z_vals = compute_contrast(labels, results, con_val).z_score()
@@ -37,8 +38,9 @@ def test_Tcontrast():
 
 
 def test_Fcontrast():
+    rng = np.random.RandomState(42)
     n, p, q = 100, 80, 10
-    X, Y = np.random.randn(p, q), np.random.randn(p, n)
+    X, Y = rng.standard_normal(size=(p, q)), rng.standard_normal(size=(p, n))
     for model in ['ols', 'ar1']:
         labels, results = run_glm(Y, X, model)
         for con_val in [np.eye(q)[0], np.eye(q)[:3]]:
@@ -49,8 +51,9 @@ def test_Fcontrast():
 
 
 def test_t_contrast_add():
+    rng = np.random.RandomState(42)
     n, p, q = 100, 80, 10
-    X, Y = np.random.randn(p, q), np.random.randn(p, n)
+    X, Y = rng.standard_normal(size=(p, q)), rng.standard_normal(size=(p, n))
     lab, res = run_glm(Y, X, 'ols')
     c1, c2 = np.eye(q)[0], np.eye(q)[1]
     con = compute_contrast(lab, res, c1) + compute_contrast(lab, res, c2)
@@ -60,8 +63,9 @@ def test_t_contrast_add():
 
 
 def test_fixed_effect_contrast():
+    rng = np.random.RandomState(42)
     n, p, q = 100, 80, 10
-    X, Y = np.random.randn(p, q), np.random.randn(p, n)
+    X, Y = rng.standard_normal(size=(p, q)), rng.standard_normal(size=(p, n))
     lab, res = run_glm(Y, X, 'ols')
     c1, c2 = np.eye(q)[0], np.eye(q)[1]
     con = _compute_fixed_effect_contrast([lab, lab], [res, res], [c1, c2])
@@ -89,8 +93,9 @@ def test_fixed_effect_contrast_nonzero_effect():
 
 
 def test_F_contrast_add():
+    rng = np.random.RandomState(42)
     n, p, q = 100, 80, 10
-    X, Y = np.random.randn(p, q), np.random.randn(p, n)
+    X, Y = rng.standard_normal(size=(p, q)), rng.standard_normal(size=(p, n))
     lab, res = run_glm(Y, X, 'ar1')
     c1, c2 = np.eye(q)[:2], np.eye(q)[2:4]
     con = compute_contrast(lab, res, c1) + compute_contrast(lab, res, c2)
@@ -107,8 +112,9 @@ def test_F_contrast_add():
 
 
 def test_contrast_mul():
+    rng = np.random.RandomState(42)
     n, p, q = 100, 80, 10
-    X, Y = np.random.randn(p, q), np.random.randn(p, n)
+    X, Y = rng.standard_normal(size=(p, q)), rng.standard_normal(size=(p, n))
     lab, res = run_glm(Y, X, 'ar1')
     for c1 in [np.eye(q)[0], np.eye(q)[:3]]:
         con1 = compute_contrast(lab, res, c1)
@@ -119,8 +125,9 @@ def test_contrast_mul():
 
 def test_contrast_values():
     # but this test is circular and should be removed
+    rng = np.random.RandomState(42)
     n, p, q = 100, 80, 10
-    X, Y = np.random.randn(p, q), np.random.randn(p, n)
+    X, Y = rng.standard_normal(size=(p, q)), rng.standard_normal(size=(p, n))
     lab, res = run_glm(Y, X, 'ar1', bins=1)
     # t test
     cval = np.eye(q)[0]
@@ -137,9 +144,10 @@ def test_contrast_values():
 
 
 def test_low_level_fixed_effects():
+    rng = np.random.RandomState(42)
     p = 100
     # X1 is some effects estimate, V1 their variance for "session 1"
-    X1, V1 = np.random.randn(p), np.ones(p)
+    X1, V1 = rng.standard_normal(size=p), np.ones(p)
     # same thing for a "session 2"
     X2, V2 = 2 * X1, 4 * V1
     # compute the fixed effects estimate, Xf, their variance Vf,

nilearn/glm/tests/test_dmtx.py

@@ -54,10 +54,11 @@ def basic_paradigm():
 
 
 def modulated_block_paradigm():
+    rng = np.random.RandomState(42)
     conditions = ['c0', 'c0', 'c0', 'c1', 'c1', 'c1', 'c2', 'c2', 'c2']
     onsets = [30, 70, 100, 10, 30, 90, 30, 40, 60]
-    durations = 5 + 5 * np.random.rand(len(onsets))
-    values = 1 + np.random.rand(len(onsets))
+    durations = 5 + 5 * rng.uniform(size=len(onsets))
+    values = 1 + rng.uniform(size=len(onsets))
     events = pd.DataFrame({'trial_type': conditions,
                            'onset': onsets,
                            'duration': durations,
@@ -66,10 +67,11 @@ def modulated_block_paradigm():
 
 
 def modulated_event_paradigm():
+    rng = np.random.RandomState(42)
     conditions = ['c0', 'c0', 'c0', 'c1', 'c1', 'c1', 'c2', 'c2', 'c2']
     onsets = [30, 70, 100, 10, 30, 90, 30, 40, 60]
     durations = 1 * np.ones(9)
-    values = 1 + np.random.rand(len(onsets))
+    values = 1 + rng.uniform(size=len(onsets))
     events = pd.DataFrame({'trial_type': conditions,
                            'onset': onsets,
                            'duration': durations,
@@ -111,20 +113,21 @@ def test_design_matrix0():
 
 def test_design_matrix0c():
     # test design matrix creation when regressors are provided manually
+    rng = np.random.RandomState(42)
     tr = 1.0
     frame_times = np.linspace(0, 127 * tr, 128)
-    ax = np.random.randn(128, 4)
+    ax = rng.standard_normal(size=(128, 4))
     _, X, names = check_design_matrix(make_first_level_design_matrix(
         frame_times, drift_model='polynomial',
         drift_order=3, add_regs=ax))
     assert_almost_equal(X[:, 0], ax[:, 0])
-    ax = np.random.randn(127, 4)
+    ax = rng.standard_normal(size=(127, 4))
     with pytest.raises(
         AssertionError,
         match="Incorrect specification of additional regressors:."
     ):
         make_first_level_design_matrix(frame_times, add_regs=ax)
-    ax = np.random.randn(128, 4)
+    ax = rng.standard_normal(size=(128, 4))
     with pytest.raises(
         ValueError,
         match="Incorrect number of additional regressor names."
@@ -136,9 +139,10 @@ def test_design_matrix0c():
 
 def test_design_matrix0d():
     # test design matrix creation when regressors are provided manually
+    rng = np.random.RandomState(42)
     tr = 1.0
     frame_times = np.linspace(0, 127 * tr, 128)
-    ax = np.random.randn(128, 4)
+    ax = rng.standard_normal(size=(128, 4))
     _, X, names = check_design_matrix(make_first_level_design_matrix(
         frame_times, drift_model='polynomial', drift_order=3, add_regs=ax))
     assert len(names) == 8
@@ -338,11 +342,12 @@ def test_design_matrix14():
 
 def test_design_matrix15():
     # basic test based on basic_paradigm, plus user supplied regressors
+    rng = np.random.RandomState(42)
     tr = 1.0
     frame_times = np.linspace(0, 127 * tr, 128)
     events = basic_paradigm()
     hrf_model = 'glover'
-    ax = np.random.randn(128, 4)
+    ax = rng.standard_normal(size=(128, 4))
     X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
                                    drift_model='polynomial', drift_order=3,
                                    add_regs=ax)
@@ -352,11 +357,12 @@ def test_design_matrix15():
 
 def test_design_matrix16():
     # Check that additional regressors are put at the right place
+    rng = np.random.RandomState(42)
     tr = 1.0
     frame_times = np.linspace(0, 127 * tr, 128)
     events = basic_paradigm()
     hrf_model = 'glover'
-    ax = np.random.randn(128, 4)
+    ax = rng.standard_normal(size=(128, 4))
     X, names = design_matrix_light(frame_times, events, hrf_model=hrf_model,
                                    drift_model='polynomial', drift_order=3,
                                    add_regs=ax)
@@ -413,11 +419,12 @@ def test_design_matrix20():
 
 def test_design_matrix21():
     # basic test on repeated names of user supplied regressors
+    rng = np.random.RandomState(42)
     tr = 1.0
     frame_times = np.linspace(0, 127 * tr, 128)
     events = basic_paradigm()
     hrf_model = 'glover'
-    ax = np.random.randn(128, 4)
+    ax = rng.standard_normal(size=(128, 4))
     with pytest.raises(ValueError):
         design_matrix_light(frame_times, events,
                             hrf_model=hrf_model, drift_model='polynomial',

nilearn/glm/tests/test_first_level.py

@@ -217,8 +217,9 @@ def test_high_level_glm_different_design_matrices():
 
 
 def test_run_glm():
+    rng = np.random.RandomState(42)
     n, p, q = 100, 80, 10
-    X, Y = np.random.randn(p, q), np.random.randn(p, n)
+    X, Y = rng.standard_normal(size=(p, q)), rng.standard_normal(size=(p, n))
 
     # Ordinary Least Squares case
     labels, results = run_glm(Y, X, 'ols')
@@ -244,9 +245,10 @@ def test_run_glm():
 
 def test_scaling():
     """Test the scaling function"""
+    rng = np.random.RandomState(42)
     shape = (400, 10)
-    u = np.random.randn(*shape)
-    mean = 100 * np.random.rand(shape[1]) + 1
+    u = rng.standard_normal(size=shape)
+    mean = 100 * rng.uniform(size=shape[1]) + 1
     Y = u + mean
     Y_, mean_ = mean_scaling(Y)
     assert_almost_equal(Y_.mean(0), 0, 5)

nilearn/glm/tests/test_hemodynamic_models.py

@@ -80,7 +80,8 @@ def test_resample_regressor_nl():
 
 def test_orthogonalize():
     """ test that the orthogonalization is OK """
-    X = np.random.randn(100, 5)
+    rng = np.random.RandomState(42)
+    X = rng.standard_normal(size=(100, 5))
     X = _orthogonalize(X)
     K = np.dot(X.T, X)
     K -= np.diag(np.diag(K))
@@ -89,7 +90,8 @@ def test_orthogonalize():
 
 def test_orthogonalize_trivial():
     """ test that the orthogonalization is OK """
-    X = np.random.randn(100)
+    rng = np.random.RandomState(42)
+    X = rng.standard_normal(size=100)
     Y = X.copy()
     X = _orthogonalize(X)
     assert_array_equal(Y, X)

nilearn/glm/tests/test_paradigm.py

@@ -22,10 +22,11 @@ def basic_paradigm():
 
 
 def modulated_block_paradigm():
+    rng = np.random.RandomState(42)
     conditions = ['c0', 'c0', 'c0', 'c1', 'c1', 'c1', 'c2', 'c2', 'c2']
     onsets = [30, 70, 100, 10, 30, 90, 30, 40, 60]
-    durations = 5 + 5 * np.random.rand(len(onsets))
-    values = np.random.rand(len(onsets))
+    durations = 5 + 5 * rng.uniform(size=len(onsets))
+    values = rng.uniform(size=len(onsets))
     events = pd.DataFrame({'name': conditions,
                            'onset': onsets,
                            'duration': durations,
@@ -34,10 +35,11 @@ def modulated_block_paradigm():
 
 
 def modulated_event_paradigm():
+    rng = np.random.RandomState(42)
     conditions = ['c0', 'c0', 'c0', 'c1', 'c1', 'c1', 'c2', 'c2', 'c2']
     onsets = [30, 70, 100, 10, 30, 90, 30, 40, 60]
     durations = 1 * np.ones(9)
-    values = np.random.rand(len(onsets))
+    values = rng.uniform(size=len(onsets))
     events = pd.DataFrame({'name': conditions,
                            'onset': onsets,
                            'durations': durations,

nilearn/glm/tests/test_regression.py

@@ -13,9 +13,9 @@
 from nilearn.glm import OLSModel, ARModel
 
 
-RNG = np.random.RandomState(20110902)
-X = RNG.standard_normal((40, 10))
-Y = RNG.standard_normal((40,))
+RNG = np.random.RandomState(42)
+X = RNG.standard_normal(size=(40, 10))
+Y = RNG.standard_normal(size=(40,))
 
 
 def test_OLS():

nilearn/glm/tests/test_second_level.py

@@ -96,8 +96,9 @@ def test_fmri_inputs():
     # Test processing of FMRI inputs
     with InTemporaryDirectory():
         # prepare fake data
+        rng = np.random.RandomState(42)
         p, q = 80, 10
-        X = np.random.randn(p, q)
+        X = rng.standard_normal(size=(p, q))
         shapes = ((7, 8, 9, 10),)
         mask, FUNCFILE, _ = write_fake_fmri_data_and_design(shapes)
         FUNCFILE = FUNCFILE[0]
@@ -173,8 +174,9 @@ def test_fmri_inputs_for_non_parametric_inference():
     # Test processing of FMRI inputs
     with InTemporaryDirectory():
         # prepare fake data
+        rng = np.random.RandomState(42)
         p, q = 80, 10
-        X = np.random.randn(p, q)
+        X = rng.standard_normal(size=(p, q))
         shapes = ((7, 8, 9, 10),)
         mask, FUNCFILE, _ = write_fake_fmri_data_and_design(shapes)
         FUNCFILE = FUNCFILE[0]
@@ -344,7 +346,8 @@ def test_second_level_contrast_computation():
             model.compute_contrast(c1, None, None, '')
         # check that passing no explicit contrast when the design
         # matrix has more than one columns raises an error
-        X = pd.DataFrame(np.random.rand(4, 2), columns=['r1', 'r2'])
+        rng = np.random.RandomState(42)
+        X = pd.DataFrame(rng.uniform(size=(4, 2)), columns=["r1", "r2"])
         model = model.fit(Y, design_matrix=X)
         with pytest.raises(ValueError):
             model.compute_contrast(None)
@@ -389,7 +392,8 @@ def test_non_parametric_inference_contrast_computation():
             non_parametric_inference(Y, X, [], 'intercept', mask)
         # check that passing no explicit contrast when the design
         # matrix has more than one columns raises an error
-        X = pd.DataFrame(np.random.rand(4, 2), columns=['r1', 'r2'])
+        rng = np.random.RandomState(42)
+        X = pd.DataFrame(rng.uniform(size=(4, 2)), columns=["r1", "r2"])
         with pytest.raises(ValueError):
             non_parametric_inference(Y, X, None)
         del func_img, FUNCFILE, neg_log_pvals_img, X, Y

nilearn/glm/tests/test_thresholding.py

@@ -12,11 +12,12 @@
 
 
 def test_fdr():
+    rng = np.random.RandomState(42)
     n = 100
     x = np.linspace(.5 / n, 1. - .5 / n, n)
     x[:10] = .0005
     x = norm.isf(x)
-    np.random.shuffle(x)
+    rng.shuffle(x)
     assert_almost_equal(fdr_threshold(x, .1), norm.isf(.0005))
     assert fdr_threshold(x, .001) == np.infty
     with pytest.raises(ValueError):