diff --git a/doc/changes/dev/14199.apichange.rst b/doc/changes/dev/14199.apichange.rst new file mode 100644 index 00000000000..adc33299e3e --- /dev/null +++ b/doc/changes/dev/14199.apichange.rst @@ -0,0 +1 @@ +Add ``rng`` parameters to APIs that use randomness, while retaining keyword-only ``seed`` and ``random_state`` for compatibility (:gh:`9233` by `Bruno Aristimunha`_). diff --git a/examples/datasets/spm_faces_dataset.py b/examples/datasets/spm_faces_dataset.py index 32df7d1a9ed..7ad8f130e6b 100644 --- a/examples/datasets/spm_faces_dataset.py +++ b/examples/datasets/spm_faces_dataset.py @@ -36,7 +36,7 @@ raw.resample(100) raw.filter(1.0, None) # high-pass reject = dict(mag=5e-12) -ica = ICA(n_components=0.95, max_iter="auto", random_state=0) +ica = ICA(n_components=0.95, max_iter="auto", rng=97) ica.fit(raw, reject=reject) # compute correlation scores, get bad indices sorted by score eog_epochs = create_eog_epochs(raw, ch_name="MRT31-2908", reject=reject) diff --git a/examples/decoding/decoding_csp_eeg.py b/examples/decoding/decoding_csp_eeg.py index 48c0e7aaf19..6ff746c98c0 100644 --- a/examples/decoding/decoding_csp_eeg.py +++ b/examples/decoding/decoding_csp_eeg.py @@ -80,7 +80,7 @@ scores = [] epochs_data = epochs.get_data(copy=False) epochs_data_train = epochs_train.get_data(copy=False) -cv = ShuffleSplit(10, test_size=0.2, random_state=42) +cv = ShuffleSplit(10, test_size=0.2, random_state=103) cv_split = cv.split(epochs_data_train) # Assemble a classifier diff --git a/examples/decoding/decoding_csp_timefreq.py b/examples/decoding/decoding_csp_timefreq.py index 9c26bf05444..1601bc5e803 100644 --- a/examples/decoding/decoding_csp_timefreq.py +++ b/examples/decoding/decoding_csp_timefreq.py @@ -53,7 +53,7 @@ LinearDiscriminantAnalysis(), ) n_splits = 3 # for cross-validation, 5 is better, here we use 3 for speed -cv = StratifiedKFold(n_splits=n_splits, shuffle=True, random_state=42) +cv = StratifiedKFold(n_splits=n_splits, shuffle=True, random_state=191) # Classification & time-frequency parameters tmin, tmax = -0.200, 2.000 diff --git a/examples/decoding/decoding_rsa.py b/examples/decoding/decoding_rsa.py index 3412b87cd74..73c0124178e 100644 --- a/examples/decoding/decoding_rsa.py +++ b/examples/decoding/decoding_rsa.py @@ -125,13 +125,13 @@ # to focus the classifier on the time interval with best SNR. clf = make_pipeline( StandardScaler(), - OneVsRestClassifier(LogisticRegression(C=1)), + OneVsRestClassifier(LogisticRegression(C=1, random_state=79)), ) X = epochs.get_data(tmin=0.05, tmax=0.3).mean(axis=2) y = epochs.events[:, 2] classes = set(y) -cv = StratifiedKFold(n_splits=5, random_state=0, shuffle=True) +cv = StratifiedKFold(n_splits=5, random_state=83, shuffle=True) # Compute confusion matrix for each cross-validation fold y_pred = np.zeros((len(y), len(classes))) @@ -173,7 +173,7 @@ chance = 0.5 # TODO VERSION: this is MDS(2, n_init=4, init='random', metric='precomputed'), # but that spelling requires scikit-learn >= 1.8 -summary, _ = smacof(chance - confusion, n_components=2, n_init=4, random_state=0) +summary, _ = smacof(chance - confusion, n_components=2, n_init=4, random_state=89) cmap = plt.colormaps["rainbow"] colors = ["r", "b"] names = list(conds["condition"].values) diff --git a/examples/decoding/decoding_spatio_temporal_source.py b/examples/decoding/decoding_spatio_temporal_source.py index f724ea97b3b..e7fcf6f9ffd 100644 --- a/examples/decoding/decoding_spatio_temporal_source.py +++ b/examples/decoding/decoding_spatio_temporal_source.py @@ -103,7 +103,7 @@ clf = make_pipeline( StandardScaler(), # z-score normalization SelectKBest(f_classif, k=500), # select features for speed - LinearModel(LogisticRegression(C=1, solver="liblinear")), + LinearModel(LogisticRegression(C=1, solver="liblinear", random_state=107)), ) time_decod = SlidingEstimator(clf, scoring="roc_auc") diff --git a/examples/decoding/decoding_spoc_CMC.py b/examples/decoding/decoding_spoc_CMC.py index 3accd5b2cd6..f24519b76fa 100644 --- a/examples/decoding/decoding_spoc_CMC.py +++ b/examples/decoding/decoding_spoc_CMC.py @@ -60,7 +60,7 @@ # Classification pipeline with SPoC spatial filtering and Ridge Regression spoc = SPoC(n_components=2, log=True, reg="oas", rank="full") -clf = make_pipeline(spoc, Ridge()) +clf = make_pipeline(spoc, Ridge(random_state=127)) # Define a two fold cross-validation cv = KFold(n_splits=2, shuffle=False) diff --git a/examples/decoding/decoding_time_generalization_conditions.py b/examples/decoding/decoding_time_generalization_conditions.py index 5dcc9d6fcea..21c6e8e50eb 100644 --- a/examples/decoding/decoding_time_generalization_conditions.py +++ b/examples/decoding/decoding_time_generalization_conditions.py @@ -69,7 +69,8 @@ # and test on all right visual vs auditory trials. clf = make_pipeline( StandardScaler(), - LogisticRegression(solver="liblinear"), # liblinear is faster than lbfgs + # liblinear is faster than lbfgs + LogisticRegression(solver="liblinear", random_state=131), ) time_gen = GeneralizingEstimator(clf, scoring="roc_auc", n_jobs=None, verbose=True) diff --git a/examples/decoding/decoding_unsupervised_spatial_filter.py b/examples/decoding/decoding_unsupervised_spatial_filter.py index 2fb1a8fec46..764214adce3 100644 --- a/examples/decoding/decoding_unsupervised_spatial_filter.py +++ b/examples/decoding/decoding_unsupervised_spatial_filter.py @@ -63,7 +63,7 @@ ############################################################################## # Transform data with PCA computed on the average ie evoked response -pca = UnsupervisedSpatialFilter(PCA(30), average=False) +pca = UnsupervisedSpatialFilter(PCA(30, random_state=193), average=False) pca_data = pca.fit_transform(X) ev = mne.EvokedArray( np.mean(pca_data, axis=0), @@ -74,7 +74,9 @@ ############################################################################## # Transform data with ICA computed on the raw epochs (no averaging) -ica = UnsupervisedSpatialFilter(FastICA(30, whiten="unit-variance"), average=False) +ica = UnsupervisedSpatialFilter( + FastICA(30, whiten="unit-variance", random_state=197), average=False +) ica_data = ica.fit_transform(X) ev1 = mne.EvokedArray( np.mean(ica_data, axis=0), diff --git a/examples/decoding/decoding_xdawn_eeg.py b/examples/decoding/decoding_xdawn_eeg.py index 3749071baaa..cf08446ffc0 100644 --- a/examples/decoding/decoding_xdawn_eeg.py +++ b/examples/decoding/decoding_xdawn_eeg.py @@ -80,7 +80,9 @@ XdawnTransformer(n_components=n_filter), Vectorizer(), MinMaxScaler(), - OneVsRestClassifier(LogisticRegression(solver="liblinear", **kwargs)), + OneVsRestClassifier( + LogisticRegression(solver="liblinear", random_state=157, **kwargs) + ), ) # Get the data and labels @@ -89,7 +91,7 @@ y = epochs.events[:, -1] # Cross validator -cv = StratifiedKFold(n_splits=10, shuffle=True, random_state=42) +cv = StratifiedKFold(n_splits=10, shuffle=True, random_state=163) # Do cross-validation preds = np.empty(len(y)) diff --git a/examples/decoding/linear_model_patterns.py b/examples/decoding/linear_model_patterns.py index 48d679ed1fd..b36f769b3ba 100644 --- a/examples/decoding/linear_model_patterns.py +++ b/examples/decoding/linear_model_patterns.py @@ -73,7 +73,8 @@ # Decoding in sensor space using a LogisticRegression classifier # -------------------------------------------------------------- -clf = LogisticRegression(solver="liblinear") # liblinear is faster than lbfgs +# liblinear is faster than lbfgs +clf = LogisticRegression(solver="liblinear", random_state=137) scaler = StandardScaler() # create a linear model with LogisticRegression @@ -127,7 +128,7 @@ Vectorizer(), # 1) vectorize across time and channels StandardScaler(), # 2) normalize features across trials LinearModel( # 3) fits a logistic regression - LogisticRegression(solver="liblinear") + LogisticRegression(solver="liblinear", random_state=139) ), ) clf.fit(X, y) diff --git a/examples/inverse/mixed_norm_inverse.py b/examples/inverse/mixed_norm_inverse.py index 70764a53973..47f76f7f63d 100644 --- a/examples/inverse/mixed_norm_inverse.py +++ b/examples/inverse/mixed_norm_inverse.py @@ -81,10 +81,10 @@ return_residual=True, return_as_dipoles=True, verbose=True, - random_state=0, + rng=251, # for this dataset we know we should use a high alpha, so avoid some # of the slower (lower) alpha values - sure_alpha_grid=np.linspace(100, 40, 10), + sure_alpha_grid=np.linspace(90, 30, 10), ) t = 0.083 diff --git a/examples/preprocessing/find_ref_artifacts.py b/examples/preprocessing/find_ref_artifacts.py index d03701e1e5a..dc67fee8275 100644 --- a/examples/preprocessing/find_ref_artifacts.py +++ b/examples/preprocessing/find_ref_artifacts.py @@ -77,7 +77,7 @@ ica_kwargs = dict( method="picard", fit_params=dict(tol=1e-4), # use a high tol here for speed - random_state=99, + rng=99, ) all_picks = mne.pick_types(raw_tog.info, meg=True, ref_meg=True) ica_tog = ICA(n_components=60, max_iter="auto", allow_ref_meg=True, **ica_kwargs) diff --git a/examples/preprocessing/ica_comparison.py b/examples/preprocessing/ica_comparison.py index d4246b80362..b3e9bfb9468 100644 --- a/examples/preprocessing/ica_comparison.py +++ b/examples/preprocessing/ica_comparison.py @@ -50,7 +50,7 @@ def run_ica(method, fit_params=None): method=method, fit_params=fit_params, max_iter="auto", - random_state=0, + rng=29, ) t0 = time() ica.fit(raw, reject=reject) diff --git a/examples/preprocessing/muscle_ica.py b/examples/preprocessing/muscle_ica.py index 8ef1e451985..79382c23417 100644 --- a/examples/preprocessing/muscle_ica.py +++ b/examples/preprocessing/muscle_ica.py @@ -37,9 +37,7 @@ # %% # Run ICA -ica = mne.preprocessing.ICA( - n_components=15, method="picard", max_iter="auto", random_state=97 -) +ica = mne.preprocessing.ICA(n_components=15, method="picard", max_iter="auto", rng=97) ica.fit(raw) # %% @@ -104,7 +102,7 @@ # Run ICA ica = mne.preprocessing.ICA( - n_components=15, method="picard", max_iter="auto", random_state=97 + n_components=15, method="picard", max_iter="auto", rng=97 ) ica.fit(raw) ica.plot_sources(raw) diff --git a/examples/simulation/plot_stc_metrics.py b/examples/simulation/plot_stc_metrics.py index 8b481aed9e6..a86694bf090 100644 --- a/examples/simulation/plot_stc_metrics.py +++ b/examples/simulation/plot_stc_metrics.py @@ -32,7 +32,7 @@ spatial_deviation_error, ) -random_state = 42 # set random state to make this example deterministic +random_state = 229 # set random state to make this example deterministic # Import sample data data_path = sample.data_path() @@ -76,7 +76,7 @@ location=location, extent=extent, subjects_dir=subjects_dir, - random_state=random_state, + rng=random_state, ) # Dipole @@ -88,7 +88,7 @@ location=location, extent=extent, subjects_dir=subjects_dir, - random_state=random_state, + rng=random_state, ) # WHAT? @@ -128,7 +128,7 @@ raw_region = raw_region.pick(picks=["eeg", "stim"], exclude="bads") cov = mne.make_ad_hoc_cov(raw_region.info) mne.simulation.add_noise( - raw_region, cov, iir_filter=[0.2, -0.2, 0.04], random_state=random_state + raw_region, cov, iir_filter=[0.2, -0.2, 0.04], rng=random_state ) # Dipole @@ -136,7 +136,7 @@ raw_dipole = raw_dipole.pick(picks=["eeg", "stim"], exclude="bads") cov = mne.make_ad_hoc_cov(raw_dipole.info) mne.simulation.add_noise( - raw_dipole, cov, iir_filter=[0.2, -0.2, 0.04], random_state=random_state + raw_dipole, cov, iir_filter=[0.2, -0.2, 0.04], rng=random_state ) ############################################################################### diff --git a/examples/simulation/simulate_evoked_data.py b/examples/simulation/simulate_evoked_data.py index 330b7814481..fe7150d9d12 100644 --- a/examples/simulation/simulate_evoked_data.py +++ b/examples/simulation/simulate_evoked_data.py @@ -67,7 +67,7 @@ def data_fun(times): fwd["src"], n_dipoles=2, times=times, - random_state=42, + rng=239, labels=labels, data_fun=data_fun, ) diff --git a/examples/simulation/simulate_raw_data.py b/examples/simulation/simulate_raw_data.py index c5248c146dc..d2d641dd04e 100644 --- a/examples/simulation/simulate_raw_data.py +++ b/examples/simulation/simulate_raw_data.py @@ -67,7 +67,7 @@ def data_fun(times): fwd = mne.read_forward_solution(fwd_fname) src = fwd["src"] stc = simulate_sparse_stc( - src, n_dipoles=n_dipoles, times=times, data_fun=data_fun, random_state=rng + src, n_dipoles=n_dipoles, times=times, data_fun=data_fun, rng=rng ) # look at our source data fig, ax = plt.subplots(1) @@ -79,9 +79,9 @@ def data_fun(times): # Simulate raw data raw_sim = simulate_raw(raw.info, [stc] * 10, forward=fwd, verbose=True) cov = make_ad_hoc_cov(raw_sim.info) -add_noise(raw_sim, cov, iir_filter=[0.2, -0.2, 0.04], random_state=rng) -add_ecg(raw_sim, random_state=rng) -add_eog(raw_sim, random_state=rng) +add_noise(raw_sim, cov, iir_filter=[0.2, -0.2, 0.04], rng=rng) +add_ecg(raw_sim, rng=rng) +add_eog(raw_sim, rng=rng) raw_sim.plot() ############################################################################## diff --git a/examples/simulation/simulated_raw_data_using_subject_anatomy.py b/examples/simulation/simulated_raw_data_using_subject_anatomy.py index fa98684ff11..b9ae6e92c38 100644 --- a/examples/simulation/simulated_raw_data_using_subject_anatomy.py +++ b/examples/simulation/simulated_raw_data_using_subject_anatomy.py @@ -204,9 +204,9 @@ def data_fun(times, latency, duration): raw_sim = mne.simulation.simulate_raw(info, source_simulator, forward=fwd) raw_sim.set_eeg_reference(projection=True) -mne.simulation.add_noise(raw_sim, cov=noise_cov, random_state=0) -mne.simulation.add_eog(raw_sim, random_state=0) -mne.simulation.add_ecg(raw_sim, random_state=0) +mne.simulation.add_noise(raw_sim, cov=noise_cov, rng=211) +mne.simulation.add_eog(raw_sim, rng=223) +mne.simulation.add_ecg(raw_sim, rng=227) # Plot original and simulated raw data. raw_sim.plot(title="Simulated raw data") diff --git a/examples/simulation/source_simulator.py b/examples/simulation/source_simulator.py index 557e32afe40..e09834d72d8 100644 --- a/examples/simulation/source_simulator.py +++ b/examples/simulation/source_simulator.py @@ -85,7 +85,7 @@ class to generate source estimates and raw data. It is meant to be a brief # simulator can be given directly to the simulate_raw function. raw = mne.simulation.simulate_raw(info, source_simulator, forward=fwd) cov = mne.make_ad_hoc_cov(raw.info) -mne.simulation.add_noise(raw, cov, iir_filter=[0.2, -0.2, 0.04], random_state=97) +mne.simulation.add_noise(raw, cov, iir_filter=[0.2, -0.2, 0.04], rng=97) raw.plot() # %% diff --git a/examples/stats/cluster_stats_evoked.py b/examples/stats/cluster_stats_evoked.py index c0a2630cac1..597e8d6286a 100644 --- a/examples/stats/cluster_stats_evoked.py +++ b/examples/stats/cluster_stats_evoked.py @@ -70,7 +70,7 @@ threshold=threshold, tail=1, n_jobs=None, - seed=0, + rng=23, out_type="mask", ) diff --git a/examples/stats/sensor_permutation_test.py b/examples/stats/sensor_permutation_test.py index 9583d262166..8507b5fdc85 100644 --- a/examples/stats/sensor_permutation_test.py +++ b/examples/stats/sensor_permutation_test.py @@ -61,7 +61,7 @@ data = np.mean(data[:, :, temporal_mask], axis=2) n_permutations = 50000 -T0, p_values, H0 = permutation_t_test(data, n_permutations, n_jobs=None, seed=0) +T0, p_values, H0 = permutation_t_test(data, n_permutations, n_jobs=None, rng=17) significant_sensors = picks[p_values <= 0.05] significant_sensors_names = [raw.ch_names[k] for k in significant_sensors] diff --git a/examples/time_frequency/time_frequency_erds.py b/examples/time_frequency/time_frequency_erds.py index 93272eb7aa3..a9e6aaead27 100644 --- a/examples/time_frequency/time_frequency_erds.py +++ b/examples/time_frequency/time_frequency_erds.py @@ -90,7 +90,7 @@ cnorm = TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax) # min, center & max ERDS kwargs = dict( - n_permutations=100, step_down_p=0.05, seed=1, buffer_size=None, out_type="mask" + n_permutations=100, step_down_p=0.05, rng=1, buffer_size=None, out_type="mask" ) # for cluster test # %% diff --git a/examples/time_frequency/time_frequency_global_field_power.py b/examples/time_frequency/time_frequency_global_field_power.py index cc4ff14ce2a..efadad4ba03 100644 --- a/examples/time_frequency/time_frequency_global_field_power.py +++ b/examples/time_frequency/time_frequency_global_field_power.py @@ -136,7 +136,7 @@ def stat_fun(x): ax.plot(times, gfp, label=freq_name, color=color, linewidth=2.5) ax.axhline(0, linestyle="--", color="grey", linewidth=2) ci_low, ci_up = bootstrap_confidence_interval( - average.data, random_state=0, stat_fun=stat_fun + average.data, rng=109, stat_fun=stat_fun ) ci_low = rescale(ci_low, average.times, baseline=(None, 0)) ci_up = rescale(ci_up, average.times, baseline=(None, 0)) diff --git a/examples/time_frequency/time_frequency_simulated.py b/examples/time_frequency/time_frequency_simulated.py index e04dda3c7d7..f9eaaeb4f50 100644 --- a/examples/time_frequency/time_frequency_simulated.py +++ b/examples/time_frequency/time_frequency_simulated.py @@ -43,7 +43,7 @@ n_times = 1024 # Just over 1 second epochs n_epochs = 40 -seed = 42 +seed = 181 rng = np.random.default_rng(seed) data = rng.standard_normal((len(ch_names), n_times * n_epochs + 200)) # buffer diff --git a/examples/visualization/channel_epochs_image.py b/examples/visualization/channel_epochs_image.py index 9281270a8c1..cef054b84da 100644 --- a/examples/visualization/channel_epochs_image.py +++ b/examples/visualization/channel_epochs_image.py @@ -74,7 +74,7 @@ def order_func(times, data): this_data /= np.sqrt(np.sum(this_data**2, axis=1))[:, np.newaxis] return np.argsort( spectral_embedding( - rbf_kernel(this_data, gamma=1.0), n_components=1, random_state=0 + rbf_kernel(this_data, gamma=1.0), n_components=1, random_state=233 ).ravel() ) diff --git a/mne/_fiff/tests/test_what.py b/mne/_fiff/tests/test_what.py index 6f407b93eab..6430940fc93 100644 --- a/mne/_fiff/tests/test_what.py +++ b/mne/_fiff/tests/test_what.py @@ -23,7 +23,7 @@ def test_what(tmp_path, verbose_debug): """Test mne.what.""" pytest.importorskip("sklearn") # ICA - ica = ICA(max_iter=1, random_state=0) + ica = ICA(max_iter=1, rng=0) raw = RawArray( np.random.default_rng(0).standard_normal((3, 10)), create_info(3, 1000.0, "eeg") ) diff --git a/mne/beamformer/tests/test_rap_music.py b/mne/beamformer/tests/test_rap_music.py index 84d5fe220c2..20cbb027643 100644 --- a/mne/beamformer/tests/test_rap_music.py +++ b/mne/beamformer/tests/test_rap_music.py @@ -71,8 +71,7 @@ def simu_data(evoked, forward, noise_cov, n_dipoles, times, nave=1): tmin, tstep = times.min(), 1 / evoked.info["sfreq"] stc = mne.SourceEstimate(data, vertices=vertices, tmin=tmin, tstep=tstep) - # noise seed chosen to keep the explained-variance and gof values well - # inside the bounds asserted in the tests + # The bounds below were calibrated against this legacy noise stream. sim_evoked = mne.simulation.simulate_evoked( forward, stc, evoked.info, noise_cov, nave=nave, random_state=106 ) diff --git a/mne/cov.py b/mne/cov.py index 200b2cbf432..80f59811b86 100644 --- a/mne/cov.py +++ b/mne/cov.py @@ -1524,12 +1524,12 @@ def _auto_low_rank_model( iter_n_components = np.arange(5, data.shape[1], 5) from sklearn.decomposition import PCA, FactorAnalysis + random_state = method_params.pop("random_state", 0) if mode == "factor_analysis": - est = FactorAnalysis + est = FactorAnalysis(random_state=random_state, **method_params) else: assert mode == "pca" - est = PCA - est = est(**method_params) + est = PCA(random_state=random_state, **method_params) est.n_components = 1 scores = np.empty_like(iter_n_components, dtype=np.float64) scores.fill(np.nan) diff --git a/mne/decoding/base.py b/mne/decoding/base.py index 371888fedf6..5c4cc9a3b6f 100644 --- a/mne/decoding/base.py +++ b/mne/decoding/base.py @@ -469,7 +469,9 @@ def __init__(self, model=None): def __sklearn_tags__(self): """Get sklearn tags.""" tags = super().__sklearn_tags__() - model = self.model if self.model is not None else LogisticRegression() + model = ( + self.model if self.model is not None else LogisticRegression(random_state=0) + ) model_tags = model.__sklearn_tags__() tags.estimator_type = model_tags.estimator_type if tags.estimator_type is not None: @@ -538,7 +540,7 @@ def fit(self, X, y, **fit_params): self.model_ = ( clone(self.model) if self.model is not None - else LogisticRegression(solver="liblinear") + else LogisticRegression(solver="liblinear", random_state=0) ) self.model_.fit(X, y, **fit_params) diff --git a/mne/decoding/tests/test_base.py b/mne/decoding/tests/test_base.py index ca367a73ff0..adc439bd026 100644 --- a/mne/decoding/tests/test_base.py +++ b/mne/decoding/tests/test_base.py @@ -98,7 +98,9 @@ def _make_data(n_samples=1000, n_features=5, n_targets=3): def test_get_coef(): """Test getting linear coefficients (filters/patterns) from estimators.""" - lm_classification = LinearModel(LogisticRegression(solver="liblinear")) + lm_classification = LinearModel( + LogisticRegression(solver="liblinear", random_state=0) + ) assert hasattr(lm_classification, "__sklearn_tags__") if check_version("sklearn", "1.6"): print(lm_classification.__sklearn_tags__()) @@ -107,7 +109,7 @@ def test_get_coef(): assert not is_regressor(lm_classification.model) assert not is_regressor(lm_classification) - lm_regression = LinearModel(Ridge()) + lm_regression = LinearModel(Ridge(random_state=0)) assert is_regressor(lm_regression.model) assert is_regressor(lm_regression) assert not is_classifier(lm_regression.model) @@ -115,7 +117,7 @@ def test_get_coef(): parameters = {"kernel": ["linear"], "C": [1, 10]} lm_gs_classification = LinearModel( - GridSearchCV(svm.SVC(), parameters, cv=2, refit=True, n_jobs=None) + GridSearchCV(svm.SVC(random_state=0), parameters, cv=2, refit=True, n_jobs=None) ) assert is_classifier(lm_gs_classification) @@ -247,19 +249,21 @@ def transform(self, X): pytest.param( make_pipeline( Scaler(info=None, scalings="mean"), - SlidingEstimator(make_pipeline(LinearModel(Ridge()))), + SlidingEstimator(make_pipeline(LinearModel(Ridge(random_state=0)))), ), id="Scaler+SlidingEstimator", ), pytest.param( make_pipeline( _Noop(), - SlidingEstimator(make_pipeline(LinearModel(Ridge()))), + SlidingEstimator(make_pipeline(LinearModel(Ridge(random_state=0)))), ), id="Noop+SlidingEstimator", ), pytest.param( - SlidingEstimator(make_pipeline(StandardScaler(), LinearModel(Ridge()))), + SlidingEstimator( + make_pipeline(StandardScaler(), LinearModel(Ridge(random_state=0))) + ), id="SlidingEstimator+nested StandardScaler", ), ], @@ -293,7 +297,11 @@ def test_get_coef_inverse_step_name(): X, y, _ = _make_data(n_samples=100, n_features=5, n_targets=1) # Test with a simple pipeline - pipe = make_pipeline(StandardScaler(), PCA(n_components=3), LinearModel(Ridge())) + pipe = make_pipeline( + StandardScaler(), + PCA(n_components=3, random_state=0), + LinearModel(Ridge(random_state=0)), + ) pipe.fit(X, y) coef_inv_actual = get_coef( @@ -325,7 +333,9 @@ def test_get_coef_inverse_step_name(): ) # Test with a nested pipeline to check __ parsing - inner_pipe = make_pipeline(PCA(n_components=3), LinearModel(Ridge())) + inner_pipe = make_pipeline( + PCA(n_components=3, random_state=0), LinearModel(Ridge(random_state=0)) + ) nested_pipe = make_pipeline(StandardScaler(), inner_pipe) nested_pipe.fit(X, y) coef_nested_inv_actual = get_coef( @@ -360,7 +370,7 @@ def transform(self, X): # In a real scenario, this would modify X return X - pipe = make_pipeline(NonInvertibleTransformer(), LinearModel(Ridge())) + pipe = make_pipeline(NonInvertibleTransformer(), LinearModel(Ridge(random_state=0))) pipe.fit(X, y) with pytest.warns(RuntimeWarning, match="not invertible"): _ = get_coef( @@ -389,7 +399,7 @@ def test_get_coef_multiclass(n_features, n_targets): assert_array_equal(lm.filters_.shape, want_shape) if n_features > 1 and n_targets > 1: assert_array_almost_equal(A, lm.patterns_.T, decimal=2) - lm = LinearModel(Ridge(alpha=0)) + lm = LinearModel(Ridge(alpha=0, random_state=0)) clf = make_pipeline(lm) clf.fit(X, Y) if n_features > 1 and n_targets > 1: @@ -400,7 +410,7 @@ def test_get_coef_multiclass(n_features, n_targets): # With epochs, scaler, and vectorizer (typical use case) X_epo = X.reshape(X.shape + (1,)) info = create_info(n_features, 1000.0, "eeg") - lm = LinearModel(Ridge(alpha=1)) + lm = LinearModel(Ridge(alpha=1, random_state=0)) clf = make_pipeline( Scaler(info, scalings=dict(eeg=1.0)), # XXX adding this step breaks Vectorizer(), @@ -490,7 +500,7 @@ def test_linearmodel(): _ = clf.fit_transform(X, y) # check that model has to have coef_, RBF-SVM doesn't - clf = LinearModel(svm.SVC(kernel="rbf")) + clf = LinearModel(svm.SVC(kernel="rbf", random_state=0)) with pytest.raises(ValueError, match="does not have a `coef_`"): clf.fit(X, y) @@ -502,7 +512,7 @@ def test_linearmodel(): # check categorical target fit in standard linear model with GridSearchCV parameters = {"kernel": ["linear"], "C": [1, 10]} clf = LinearModel( - GridSearchCV(svm.SVC(), parameters, cv=2, refit=True, n_jobs=None) + GridSearchCV(svm.SVC(random_state=0), parameters, cv=2, refit=True, n_jobs=None) ) clf.fit(X, y) assert_equal(clf.filters_.shape, (n_features,)) @@ -596,7 +606,7 @@ def test_cross_val_multiscore(): assert_array_equal(manual, auto) -@parametrize_with_checks([LinearModel(LogisticRegression())]) +@parametrize_with_checks([LinearModel(LogisticRegression(random_state=0))]) def test_sklearn_compliance(estimator, check): """Test LinearModel compliance with sklearn.""" check(estimator) diff --git a/mne/decoding/tests/test_csp.py b/mne/decoding/tests/test_csp.py index 1110f0ddabb..04fdd6956fc 100644 --- a/mne/decoding/tests/test_csp.py +++ b/mne/decoding/tests/test_csp.py @@ -365,7 +365,7 @@ def test_regularized_csp(ch_type, rank, reg): clf = make_pipeline( sc, csp, - LinearModel(LogisticRegression(solver="liblinear")), + LinearModel(LogisticRegression(solver="liblinear", random_state=0)), ) score = cross_val_score(clf, epochs_data_orig, y, cv=cv, scoring="roc_auc").mean() assert 0.75 <= score <= 1.0 @@ -387,7 +387,7 @@ def test_regularized_csp(ch_type, rank, reg): def test_csp_pipeline(): """Test if CSP works in a pipeline.""" csp = CSP(reg=1, norm_trace=False) - svc = SVC() + svc = SVC(random_state=0) pipe = Pipeline([("CSP", csp), ("SVC", svc)]) pipe.set_params(CSP__reg=0.2) assert pipe.get_params()["CSP__reg"] == 0.2 diff --git a/mne/decoding/tests/test_receptive_field.py b/mne/decoding/tests/test_receptive_field.py index 01df9409e7e..bef99f12887 100644 --- a/mne/decoding/tests/test_receptive_field.py +++ b/mne/decoding/tests/test_receptive_field.py @@ -99,7 +99,7 @@ def test_rank_deficiency(): y = np.apply_along_axis(np.convolve, 0, eeg, win, mode="same") y += rng.normal(scale=100, size=y.shape) - for est in (Ridge(reg), reg): + for est in (Ridge(reg, random_state=0), reg): rf = ReceptiveField(tmin, tmax, fs, estimator=est, patterns=True) rf.fit(eeg, y) pred = rf.predict(eeg) @@ -179,7 +179,7 @@ def test_time_delay(): def test_receptive_field_basic(n_jobs): """Test model prep and fitting.""" # Make sure estimator pulling works - mod = Ridge() + mod = Ridge(random_state=0) rng = np.random.default_rng(1337) # Test the receptive field model @@ -401,7 +401,13 @@ def test_receptive_field_1d(n_jobs): fit_intercept=False, n_jobs=n_jobs, ) - for estimator in (Ridge(alpha=0.0), Ridge(alpha=0.1), 0.0, 0.1, lap): + for estimator in ( + Ridge(alpha=0.0, random_state=0), + Ridge(alpha=0.1, random_state=0), + 0.0, + 0.1, + lap, + ): for offset in (-100, 0, 100): model = ReceptiveField( smin, smax, 1.0, estimator=estimator, n_jobs=n_jobs @@ -456,7 +462,8 @@ def test_receptive_field_nd(n_jobs): smin, smax, 1.0, 0.1, n_jobs=n_jobs, edge_correction=False ) for estimator, atol in zip( - (Ridge(alpha=0.0), 0.0, 0.01, tdr_l, tdr_nc), (1e-3, 1e-3, 1e-3, 5e-3, 5e-2) + (Ridge(alpha=0.0, random_state=0), 0.0, 0.01, tdr_l, tdr_nc), + (1e-3, 1e-3, 1e-3, 5e-3, 5e-2), ): model = ReceptiveField(smin, smax, 1.0, estimator=estimator) model.fit(x, y) @@ -480,9 +487,9 @@ def test_receptive_field_nd(n_jobs): tdr = TimeDelayingRidge(smin, smax, 1.0, 0.0, n_jobs=n_jobs) tdr_no = TimeDelayingRidge(smin, smax, 1.0, 0.0, fit_intercept=False, n_jobs=n_jobs) for estimator in ( - Ridge(alpha=0.0), + Ridge(alpha=0.0, random_state=0), tdr, - Ridge(alpha=0.0, fit_intercept=False), + Ridge(alpha=0.0, fit_intercept=False, random_state=0), tdr_no, ): # first with no intercept in the data @@ -558,7 +565,7 @@ def test_inverse_coef(): # Check coefficient dims, for all estimator types X, y = _make_data(n_feats, n_targets, n_samples, tmin, tmax) tdr = TimeDelayingRidge(tmin, tmax, 1.0, 0.1, "laplacian") - for estimator in (0.0, 0.01, Ridge(alpha=0.0), tdr): + for estimator in (0.0, 0.01, Ridge(alpha=0.0, random_state=0), tdr): rf = ReceptiveField(tmin, tmax, 1.0, estimator=estimator, patterns=True) rf.fit(X, y) inv_rf = ReceptiveField(tmin, tmax, 1.0, estimator=estimator, patterns=True) @@ -581,7 +588,7 @@ def test_linalg_warning(): """Test that warnings are issued when no regularization is applied.""" n_feats, n_targets, n_samples = 5, 60, 50 X, y = _make_data(n_feats, n_targets, n_samples, tmin, tmax) - for estimator in (0.0, Ridge(alpha=0.0)): + for estimator in (0.0, Ridge(alpha=0.0, random_state=0)): rf = ReceptiveField(tmin, tmax, 1.0, estimator=estimator) with pytest.warns( (RuntimeWarning, UserWarning), match="[Singular|scipy.linalg.solve]" @@ -605,7 +612,9 @@ def test_tdr_sklearn_compliance(estimator, check): @pytest.mark.filterwarnings("ignore:.*invalid value encountered in subtract.*:") -@parametrize_with_checks([ReceptiveField(-1, 2, 1.0, estimator=Ridge(), patterns=True)]) +@parametrize_with_checks( + [ReceptiveField(-1, 2, 1.0, estimator=Ridge(random_state=0), patterns=True)] +) def test_rf_sklearn_compliance(estimator, check): """Test sklearn RF compliance.""" pytest.importorskip("sklearn", minversion="1.6") # TODO VERSION remove on 1.6+ diff --git a/mne/decoding/tests/test_search_light.py b/mne/decoding/tests/test_search_light.py index 162983a71da..934386b07a1 100644 --- a/mne/decoding/tests/test_search_light.py +++ b/mne/decoding/tests/test_search_light.py @@ -52,9 +52,9 @@ def test_search_light_basic(): sl = SlidingEstimator("foo") with pytest.raises(ValueError, match="must be"): sl.fit(X, y) - sl = SlidingEstimator(Ridge()) + sl = SlidingEstimator(Ridge(random_state=0)) assert not is_classifier(sl) - sl = SlidingEstimator(LogisticRegression(solver="liblinear")) + sl = SlidingEstimator(LogisticRegression(solver="liblinear", random_state=0)) assert is_classifier(sl.base_estimator) assert is_classifier(sl) # fit @@ -183,7 +183,9 @@ def transform(self, X): # Bagging classifiers X = rng.random((10, 3, 4)) for n_jobs in (1, 2): - pipe = SlidingEstimator(BaggingClassifier(None, 2), n_jobs=n_jobs) + pipe = SlidingEstimator( + BaggingClassifier(None, 2, random_state=0), n_jobs=n_jobs + ) pipe.fit(X, y) pipe.score(X, y) assert isinstance(pipe.estimators_[0], BaggingClassifier) @@ -324,9 +326,15 @@ def test_gl_score_branches(scoring, est_name, method): ) elif scoring == "accuracy_kwargs": # start from the default scorer but add a kwarg to prevent batching - acc_func = check_scoring(LogisticRegression(), "accuracy")._score_func + acc_func = check_scoring( + LogisticRegression(random_state=0), "accuracy" + )._score_func scoring = make_scorer(acc_func, normalize=False) - est = Ridge() if est_name == "ridge" else LogisticRegression(solver=solver) + est = ( + Ridge(random_state=0) + if est_name == "ridge" + else LogisticRegression(solver=solver, random_state=0) + ) gl = GeneralizingEstimator(est, scoring=scoring).fit(X, y) # Measure batching: count pred/call scores. Wraps `fn` calls so they append @@ -374,7 +382,7 @@ def force_fallback(e, X, y): def test_verbose_arg(capsys, n_jobs, verbose): """Test controlling output with the ``verbose`` argument.""" X, y = make_data() - clf = SVC() + clf = SVC(random_state=0) # shows progress bar and prints other messages to the console with use_log_level(True): @@ -424,8 +432,8 @@ def predict_proba(self, X): @pytest.mark.slowtest @parametrize_with_checks( [ - SlidingEstimator(LogisticRegression(), allow_2d=True), - GeneralizingEstimator(LogisticRegression(), allow_2d=True), + SlidingEstimator(LogisticRegression(random_state=0), allow_2d=True), + GeneralizingEstimator(LogisticRegression(random_state=0), allow_2d=True), ] ) def test_sklearn_compliance(estimator, check): diff --git a/mne/decoding/tests/test_transformer.py b/mne/decoding/tests/test_transformer.py index a91a4316521..f6d4bcc6031 100644 --- a/mne/decoding/tests/test_transformer.py +++ b/mne/decoding/tests/test_transformer.py @@ -262,9 +262,9 @@ def test_unsupervised_spatial_filter(): # Test fit n_components = 4 - usf = UnsupervisedSpatialFilter(PCA(n_components)) + usf = UnsupervisedSpatialFilter(PCA(n_components, random_state=0)) usf.fit(X) - usf1 = UnsupervisedSpatialFilter(PCA(n_components)) + usf1 = UnsupervisedSpatialFilter(PCA(n_components, random_state=0)) # test transform assert_equal(usf.transform(X).ndim, 3) @@ -274,9 +274,9 @@ def test_unsupervised_spatial_filter(): assert_array_almost_equal(usf.inverse_transform(usf.transform(X)), X) # Test with average param - usf = UnsupervisedSpatialFilter(PCA(4), average=True) + usf = UnsupervisedSpatialFilter(PCA(4, random_state=0), average=True) usf.fit_transform(X) - usf = UnsupervisedSpatialFilter(PCA(4), 2) + usf = UnsupervisedSpatialFilter(PCA(4, random_state=0), 2) with pytest.raises(TypeError, match="average must be"): usf.fit(X) @@ -335,7 +335,7 @@ def test_bad_triage(): Scaler(scalings="mean"), # Not easy to test Scaler(info) b/c number of channels must match TemporalFilter(), - UnsupervisedSpatialFilter(PCA()), + UnsupervisedSpatialFilter(PCA(random_state=0)), Vectorizer(), ] ) diff --git a/mne/decoding/transformer.py b/mne/decoding/transformer.py index 6fdb2f272e3..79656d40083 100644 --- a/mne/decoding/transformer.py +++ b/mne/decoding/transformer.py @@ -307,7 +307,9 @@ class Vectorizer(MNETransformerMixin, BaseEstimator): >>> from sklearn.linear_model import LogisticRegression >>> from sklearn.pipeline import make_pipeline >>> from sklearn.preprocessing import StandardScaler - >>> clf = make_pipeline(Vectorizer(), StandardScaler(), LogisticRegression()) + >>> clf = make_pipeline( + ... Vectorizer(), StandardScaler(), LogisticRegression(random_state=0) + ... ) """ def fit(self, X, y=None): diff --git a/mne/epochs.py b/mne/epochs.py index fa0292041b1..f313a94e2ac 100644 --- a/mne/epochs.py +++ b/mne/epochs.py @@ -90,6 +90,7 @@ _convert_times, _ensure_events, _gen_events, + _legacy_rng, _on_missing, _path_like, _pl, @@ -2488,12 +2489,14 @@ def export( export_epochs(fname, self, fmt, overwrite=overwrite, verbose=verbose) + @_legacy_rng("random_state") @fill_doc def equalize_event_counts( self, event_ids: list | dict | None = None, method: Literal["truncate", "mintime", "random"] = "mintime", *, + rng=None, random_state: int | RandomState | None = None, ) -> tuple: """Equalize the number of trials in each condition. @@ -2536,7 +2539,8 @@ def equalize_event_counts( The ``event_ids`` must identify non-overlapping subsets of the epochs. %(equalize_events_method)s - %(random_state)s Used only if ``method='random'``. + %(rng_method_random)s + %(random_state_rng_method_random)s Returns ------- @@ -2640,7 +2644,10 @@ def equalize_event_counts( eq_inds.append(self._keys_to_idx(eq)) sample_nums = [self.events[e, 0] for e in eq_inds] - indices = _get_drop_indices(sample_nums, method, random_state) + legacy_seed = ( + random_state if isinstance(random_state, int | np.integer) else None + ) + indices = _get_drop_indices(sample_nums, method, rng, legacy_seed=legacy_seed) # need to re-index indices indices = np.concatenate([e[idx] for e, idx in zip(eq_inds, indices)]) self.drop(indices, reason="EQUALIZED_COUNT") @@ -4011,11 +4018,13 @@ def combine_event_ids( return epochs +@_legacy_rng("random_state") @fill_doc def equalize_epoch_counts( epochs_list: list, method: Literal["truncate", "mintime", "random"] = "mintime", *, + rng=None, random_state: int | RandomState | None = None, ) -> None: """Equalize the number of trials in multiple Epochs or EpochsTFR instances. @@ -4025,7 +4034,8 @@ def equalize_epoch_counts( epochs_list : list of Epochs The Epochs instances to equalize trial counts for. %(equalize_events_method)s - %(random_state)s Used only if ``method='random'``. + %(rng_method_random)s + %(random_state_rng_method_random)s Notes ----- @@ -4052,12 +4062,13 @@ def equalize_epoch_counts( if not epoch._bad_dropped: epoch.drop_bad() sample_nums = [epoch.events[:, 0] for epoch in epochs_list] - indices = _get_drop_indices(sample_nums, method, random_state) + legacy_seed = random_state if isinstance(random_state, int | np.integer) else None + indices = _get_drop_indices(sample_nums, method, rng, legacy_seed=legacy_seed) for epoch, inds in zip(epochs_list, indices): epoch.drop(inds, reason="EQUALIZED_COUNT") -def _get_drop_indices(sample_nums, method, random_state): +def _get_drop_indices(sample_nums, method, rng, *, legacy_seed=None): """Get indices to drop from multiple event timing lists.""" small_idx = np.argmin([e.size for e in sample_nums]) small_epoch_indices = sample_nums[small_idx] @@ -4070,9 +4081,14 @@ def _get_drop_indices(sample_nums, method, random_state): mask = np.ones(event.size, dtype=bool) mask[small_epoch_indices.size :] = False elif method == "random": - rng = check_random_state(random_state) mask = np.zeros(event.size, dtype=bool) - idx = rng.choice( + # Historically an integer seed was normalized inside this loop, + # restarting the same stream for every event list. Preserve that + # behavior only for the deprecated parameter; ``rng`` advances. + this_rng = ( + check_random_state(legacy_seed) if legacy_seed is not None else rng + ) + idx = this_rng.choice( np.arange(event.size), size=small_epoch_indices.size, replace=False ) mask[idx] = True @@ -4653,15 +4669,17 @@ def _get_epoch_from_raw(self, idx, verbose=None): return data +@_legacy_rng("random_state") @fill_doc -def bootstrap(epochs, random_state=None): +def bootstrap(epochs, *, rng=None, random_state=None): """Compute epochs selected by bootstrapping. Parameters ---------- epochs : Epochs instance epochs data to be bootstrapped - %(random_state)s + %(rng)s + %(random_state_rng)s Returns ------- @@ -4675,7 +4693,6 @@ def bootstrap(epochs, random_state=None): "in the constructor." ) - rng = check_random_state(random_state) epochs_bootstrap = epochs.copy() n_events = len(epochs_bootstrap.events) idx = rng_uniform(rng)(0, n_events, n_events) diff --git a/mne/filter.py b/mne/filter.py index 6bc00186001..03f614169ce 100644 --- a/mne/filter.py +++ b/mne/filter.py @@ -2092,9 +2092,9 @@ def detrend(x, order=1, axis=-1): -------- As in :func:`scipy.signal.detrend`:: - >>> randgen = np.random.RandomState(9) + >>> rng = np.random.default_rng(11) >>> npoints = int(1e3) - >>> noise = randgen.randn(npoints) + >>> noise = rng.standard_normal(npoints) >>> x = 3 + 2*np.linspace(0, 1, npoints) + noise >>> bool((detrend(x) - noise).max() < 0.01) True diff --git a/mne/forward/tests/test_make_forward.py b/mne/forward/tests/test_make_forward.py index 42cd29ca425..7244aab406b 100644 --- a/mne/forward/tests/test_make_forward.py +++ b/mne/forward/tests/test_make_forward.py @@ -767,7 +767,7 @@ def test_make_forward_dipole(tmp_path): times, pos, amplitude, ori, gof = [], [], [], [], [] nave = 400 # add a tiny amount of noise to the simulated evokeds for s in stc: - evo_test = simulate_evoked(fwd, s, info, cov, nave=nave, random_state=rng) + evo_test = simulate_evoked(fwd, s, info, cov, nave=nave, rng=rng) # evo_test.add_proj(make_eeg_average_ref_proj(evo_test.info)) dfit, resid = fit_dipole(evo_test, cov, sphere, None) times += dfit.times.tolist() diff --git a/mne/inverse_sparse/mxne_debiasing.py b/mne/inverse_sparse/mxne_debiasing.py index 8dfc6054c37..c64dd43ce76 100644 --- a/mne/inverse_sparse/mxne_debiasing.py +++ b/mne/inverse_sparse/mxne_debiasing.py @@ -6,11 +6,11 @@ import numpy as np -from ..utils import check_random_state, fill_doc, logger, verbose +from ..utils import _check_rng, fill_doc, logger, verbose @fill_doc -def power_iteration_kron(A, C, max_iter=1000, tol=1e-3, random_state=0): +def power_iteration_kron(A, C, max_iter=1000, tol=1e-3, rng=0): """Find the largest singular value for the matrix kron(C.T, A). It uses power iterations. @@ -23,7 +23,7 @@ def power_iteration_kron(A, C, max_iter=1000, tol=1e-3, random_state=0): An array max_iter : int Maximum number of iterations - %(random_state)s + %(rng)s Returns ------- @@ -35,7 +35,7 @@ def power_iteration_kron(A, C, max_iter=1000, tol=1e-3, random_state=0): http://en.wikipedia.org/wiki/Power_iteration """ AS_size = C.shape[0] - rng = check_random_state(random_state) + rng = _check_rng(rng) B = rng.standard_normal((AS_size, AS_size)) B /= np.linalg.norm(B, "fro") ATA = np.dot(A.T, A) diff --git a/mne/inverse_sparse/mxne_inverse.py b/mne/inverse_sparse/mxne_inverse.py index 5c183ad6265..b6207444394 100644 --- a/mne/inverse_sparse/mxne_inverse.py +++ b/mne/inverse_sparse/mxne_inverse.py @@ -18,8 +18,8 @@ from ..utils import ( _check_depth, _check_option, + _legacy_rng, _validate_type, - check_random_state, logger, sum_squared, verbose, @@ -341,6 +341,7 @@ def make_stc_from_dipoles(dipoles, src, verbose=None): return stc +@_legacy_rng("random_state") @verbose def mixed_norm( evoked, @@ -364,8 +365,10 @@ def mixed_norm( rank=None, pick_ori=None, sure_alpha_grid="auto", - random_state=None, verbose=None, + *, + rng=None, + random_state=None, ): """Mixed-norm estimate (MxNE) and iterative reweighted MxNE (irMxNE). @@ -432,12 +435,12 @@ def mixed_norm( grid is directly specified. Ignored if alpha is not "sure". .. versionadded:: 0.24 - random_state : int | None - The random state used in a random number generator for delta and - epsilon used for the SURE computation. Defaults to None. + %(verbose)s + %(rng)s + %(random_state_rng)s + Used for the random delta and epsilon in the SURE computation. .. versionadded:: 0.24 - %(verbose)s Returns ------- @@ -541,7 +544,7 @@ def mixed_norm( gain, alpha_grid, sigma=1, - random_state=random_state, + rng=rng, n_mxne_iter=n_mxne_iter, maxit=maxit, tol=tol, @@ -926,7 +929,7 @@ def _compute_mxne_sure( debias, solver, dgap_freq, - random_state, + rng, verbose, ): """Stein Unbiased Risk Estimator (SURE). @@ -964,9 +967,9 @@ def _compute_mxne_sure( The algorithm to use for the optimization. dgap_freq : int or np.inf The duality gap is evaluated every dgap_freq iterations. - random_state : int | None - The random state used in a random number generator for delta and - epsilon used for the SURE computation. + rng : instance of numpy.random.Generator + The random number generator used for delta and epsilon in the SURE + computation. Returns ------- @@ -1072,7 +1075,6 @@ def _compute_sure_val(coef1, coef2, gain, M, sigma, delta, eps): sure_path = np.empty(len(alpha_grid)) - rng = check_random_state(random_state) # See Deledalle et al. 20214 Sec. 5.1 eps = 2 * sigma / (M.shape[0] ** 0.3) delta = rng.standard_normal(M.shape) diff --git a/mne/inverse_sparse/mxne_optim.py b/mne/inverse_sparse/mxne_optim.py index 5a4f1ccefae..b197be348ef 100644 --- a/mne/inverse_sparse/mxne_optim.py +++ b/mne/inverse_sparse/mxne_optim.py @@ -154,6 +154,7 @@ def _mixed_norm_solver_cd( tol=tol / sum_squared(M), fit_intercept=False, max_iter=maxit, + random_state=0, warm_start=True, ) if init is not None: diff --git a/mne/inverse_sparse/tests/test_mxne_inverse.py b/mne/inverse_sparse/tests/test_mxne_inverse.py index d83a5cbf8d8..7cce3daef02 100644 --- a/mne/inverse_sparse/tests/test_mxne_inverse.py +++ b/mne/inverse_sparse/tests/test_mxne_inverse.py @@ -370,7 +370,7 @@ def test_mxne_vol_sphere(): cov, nave=1e9, use_cps=True, - random_state=0, + rng=np.random.default_rng(0), ) dip_mxne = mixed_norm( @@ -547,7 +547,7 @@ def test_mxne_inverse_sure_synthetic( debias=True, solver="auto", dgap_freq=10, - random_state=0, + rng=np.random.default_rng(0), verbose=False, ) assert np.count_nonzero(active_set, axis=-1) == n_orient * nnz @@ -604,7 +604,13 @@ def data_fun(times): ) evoked = evoked.crop(tmin=0, tmax=10e-3) stc_ = mixed_norm( - evoked, forward, noise_cov, loose=0.9, n_mxne_iter=5, depth=0.9, random_state=1 + evoked, + forward, + noise_cov, + loose=0.9, + n_mxne_iter=5, + depth=0.9, + random_state=1, ) assert len(stc_.vertices) == len(stc.vertices) == 2 for si in range(len(stc_.vertices)): @@ -633,7 +639,7 @@ def test_mxne_inverse_empty(): n_mxne_iter=3, alpha=99, return_residual=True, - random_state=0, + rng=0, ) assert stc.data.size == 0 assert stc.vertices[0].size == 0 diff --git a/mne/label.py b/mne/label.py index df07faf1956..6907fac05c2 100644 --- a/mne/label.py +++ b/mne/label.py @@ -41,8 +41,8 @@ _check_option, _check_subject, _import_nibabel, + _legacy_rng, _validate_type, - check_random_state, fill_doc, get_subjects_dir, logger, @@ -1995,9 +1995,17 @@ def _grow_nonoverlapping_labels( return labels +@_legacy_rng("random_state") @fill_doc def random_parcellation( - subject, n_parcel, hemi, subjects_dir=None, surface="white", random_state=None + subject, + n_parcel, + hemi, + subjects_dir=None, + surface="white", + *, + rng=None, + random_state=None, ): """Generate random cortex parcellation by growing labels. @@ -2016,7 +2024,8 @@ def random_parcellation( parcels per hemisphere. %(subjects_dir)s %(surface)s - %(random_state)s + %(rng)s + %(random_state_rng)s Returns ------- @@ -2036,7 +2045,7 @@ def random_parcellation( dist[hemi] = mesh_dist(tris[hemi], vert[hemi]) # create the patches - labels = _cortex_parcellation(subject, n_parcel, hemis, vert, dist, random_state) + labels = _cortex_parcellation(subject, n_parcel, hemis, vert, dist, rng) # add a unique color to each label colors = _n_colors(len(labels)) @@ -2046,12 +2055,9 @@ def random_parcellation( return labels -def _cortex_parcellation( - subject, n_parcel, hemis, vertices_, graphs, random_state=None -): +def _cortex_parcellation(subject, n_parcel, hemis, vertices_, graphs, rng): """Random cortex parcellation.""" labels = [] - rng = check_random_state(random_state) for hemi in set(hemis): parcel_size = len(hemis) * len(vertices_[hemi]) // n_parcel graph = graphs[hemi] # distance graph @@ -2936,6 +2942,7 @@ def write_labels_to_annot( _write_annot(fname, annot, ctab, hemi_names, table_name) +@_legacy_rng("random_state") @fill_doc def select_sources( subject, @@ -2945,8 +2952,10 @@ def select_sources( grow_outside=True, subjects_dir=None, name=None, - random_state=None, surf="white", + *, + rng=None, + random_state=None, ): """Select sources from a label. @@ -2970,9 +2979,10 @@ def select_sources( %(subjects_dir)s name : None | str Assign name to the new label. - %(random_state)s surf : str The surface used to simulated the label, defaults to the white surface. + %(rng)s + %(random_state_rng)s Returns ------- @@ -3010,7 +3020,6 @@ def select_sources( subject, restrict_vertices=True, subjects_dir=subjects_dir, surf=surf ) else: - rng = check_random_state(random_state) seed = rng.choice(label.vertices) else: seed = label.vertices[location] diff --git a/mne/preprocessing/ica.py b/mne/preprocessing/ica.py index 435d88daf7e..1c706979460 100644 --- a/mne/preprocessing/ica.py +++ b/mne/preprocessing/ica.py @@ -64,6 +64,7 @@ _check_on_missing, _check_option, _check_preload, + _check_rng, _ensure_int, _get_inst_data, _limit_blas_threads, @@ -190,6 +191,13 @@ def _check_for_unsupported_ica_channels(picks, info, allow_ref_meg=False): _KNOWN_ICA_METHODS = ("fastica", "infomax", "picard") +def _rng_to_seed(rng): + """Adapt a Generator for third-party random_state parameters.""" + if isinstance(rng, np.random.Generator): + return int(rng.integers(np.iinfo(np.int32).max)) + return rng + + @fill_doc class ICA(ContainsMixin): """Data decomposition using Independent Component Analysis (ICA). @@ -239,7 +247,8 @@ class ICA(ContainsMixin): Noise covariance used for pre-whitening. If None (default), channels are scaled to unit variance ("z-standardized") as a group by channel type prior to the whitening by PCA. - %(random_state)s + %(rng)s + %(random_state_rng)s method : 'fastica' | 'infomax' | 'picard' The ICA method to use in the fit method. Use the ``fit_params`` argument to set additional parameters. Specifically, if you want Extended @@ -432,6 +441,7 @@ def __init__( n_components=None, *, noise_cov=None, + rng=None, random_state=None, method="fastica", fit_params=None, @@ -466,7 +476,14 @@ def __init__( self._max_pca_components = None self.n_pca_components = None self.ch_names = None + if rng is not None and random_state is not None: + raise TypeError("Specify only one of rng or random_state") + if random_state is not None: + logger.info("Use rng= instead of random_state= in new code") self.random_state = random_state + # stored un-normalized so that integer seeds stay intact for the + # third-party ``random_state`` parameters used during fitting + self.rng = rng if fit_params is None: fit_params = {} @@ -887,7 +904,14 @@ def _fit(self, data, fit_type): if not np.isfinite(data).all(): raise ValueError("Input data contains non-finite values (NaN/Inf). ") - random_state = check_random_state(self.random_state) + rng = getattr(self, "rng", None) + if self.random_state is None: + # Keep integer seeds intact for third-party ``random_state`` + # parameters, but use an independent Generator for the default. + if rng is None: + rng = _check_rng(None) + else: + rng = check_random_state(self.random_state) n_channels, n_samples = data.shape self._compute_pre_whitener(data) data = self._pre_whiten(data) @@ -955,14 +979,16 @@ def _fit(self, data, fit_type): if self.method == "fastica": from sklearn.decomposition import FastICA - ica = FastICA(whiten=False, random_state=random_state, **self.fit_params) + ica = FastICA( + whiten=False, random_state=_rng_to_seed(rng), **self.fit_params + ) ica.fit(data[:, sel]) self.unmixing_matrix_ = ica.components_ self.n_iter_ = ica.n_iter_ elif self.method in ("infomax", "extended-infomax"): unmixing_matrix, n_iter = infomax( data[:, sel], - random_state=random_state, + rng=_check_rng(rng), return_n_iter=True, **self.fit_params, ) @@ -976,7 +1002,7 @@ def _fit(self, data, fit_type): data[:, sel].T, whiten=False, return_n_iter=True, - random_state=random_state, + random_state=_rng_to_seed(rng), **self.fit_params, ) self.unmixing_matrix_ = W diff --git a/mne/preprocessing/infomax_.py b/mne/preprocessing/infomax_.py index 0e4f2a12ff3..64e4047eddb 100644 --- a/mne/preprocessing/infomax_.py +++ b/mne/preprocessing/infomax_.py @@ -7,9 +7,11 @@ import numpy as np from scipy.special import expit -from ..utils import check_random_state, logger, random_permutation, verbose +from ..utils import _legacy_rng, logger, verbose +from ..utils.numerics import _random_permutation +@_legacy_rng("random_state") @verbose def infomax( data, @@ -24,13 +26,15 @@ def infomax( kurt_size=6000, ext_blocks=1, max_iter=200, - random_state=None, blowup=1e4, blowup_fac=0.5, n_small_angle=20, use_bias=True, verbose=None, return_n_iter=False, + *, + rng=None, + random_state=None, ): """Run (extended) Infomax ICA decomposition on raw data. @@ -78,7 +82,6 @@ def infomax( Defaults to 1. max_iter : int The maximum number of iterations. Defaults to 200. - %(random_state)s blowup : float The maximum difference allowed between two successive estimations of the unmixing matrix. Defaults to 10000. @@ -98,6 +101,8 @@ def infomax( return_n_iter : bool Whether to return the number of iterations performed. Defaults to False. + %(rng)s + %(random_state_rng)s Returns ------- @@ -117,8 +122,6 @@ def infomax( """ from scipy.stats import kurtosis - rng = check_random_state(random_state) - # define some default parameters max_weight = 1e8 restart_fac = 0.9 @@ -182,7 +185,7 @@ def infomax( olddelta, oldchange = 1.0, 0.0 while step < max_iter: # shuffle data at each step - permute = random_permutation(n_samples, rng) + permute = _random_permutation(n_samples, rng) # ICA training block # loop across block samples diff --git a/mne/preprocessing/tests/test_eeglab_infomax.py b/mne/preprocessing/tests/test_eeglab_infomax.py index dfa8c9a748b..3e4d2537230 100644 --- a/mne/preprocessing/tests/test_eeglab_infomax.py +++ b/mne/preprocessing/tests/test_eeglab_infomax.py @@ -33,7 +33,7 @@ def generate_data_for_comparing_against_eeglab_infomax(ch_type, random_state): # select a small number of channels for the test number_of_channels_to_use = 5 - idx_perm = random_permutation(picks.shape[0], random_state) + idx_perm = random_permutation(picks.shape[0], random_state=random_state) picks = picks[idx_perm[:number_of_channels_to_use]] raw.filter( diff --git a/mne/preprocessing/tests/test_ica.py b/mne/preprocessing/tests/test_ica.py index b3487fc49ed..dacca6dccac 100644 --- a/mne/preprocessing/tests/test_ica.py +++ b/mne/preprocessing/tests/test_ica.py @@ -54,7 +54,7 @@ read_ica_eeglab, ) from mne.rank import _compute_rank_int -from mne.utils import _record_warnings, catch_logging, check_version +from mne.utils import _record_warnings, catch_logging, check_random_state, check_version data_dir = Path(__file__).parents[2] / "io" / "tests" / "data" raw_fname = data_dir / "test_raw.fif" @@ -82,8 +82,8 @@ def ICA(*args, **kwargs): """Fix the random state in tests.""" - if "random_state" not in kwargs: - kwargs["random_state"] = 0 + if "random_state" not in kwargs and "rng" not in kwargs: + kwargs["rng"] = 0 return _ICA(*args, **kwargs) @@ -117,9 +117,7 @@ def test_ica_full_data_recovery(method): for method in methods: stuff = [(2, n_channels, True), (2, n_channels // 2, False)] for n_components, n_pca_components, ok in stuff: - ica = ICA( - n_components=n_components, random_state=0, method=method, max_iter=1 - ) + ica = ICA(n_components=n_components, rng=0, method=method, max_iter=1) kwargs = dict(exclude=[], n_pca_components=n_pca_components) picks = list(range(n_channels)) with pytest.warns(UserWarning, match=None): # sometimes warns @@ -134,7 +132,7 @@ def test_ica_full_data_recovery(method): diff = np.abs(data[:n_channels] - raw2._data[:n_channels]) assert np.max(diff) > 1e-14 - ica = ICA(n_components=n_components, method=method, random_state=0) + ica = ICA(n_components=n_components, method=method, rng=0) with _record_warnings(): # sometimes warns ica.fit(epochs, picks=picks) _assert_ica_attributes(ica, epochs.get_data(picks)) @@ -171,7 +169,7 @@ def test_ica_simple(method): data = np.dot(A, S) info = create_info(data.shape[-2], 1000.0, "eeg") cov = make_ad_hoc_cov(info) - ica = ICA(n_components=n_components, method=method, random_state=0, noise_cov=cov) + ica = ICA(n_components=n_components, method=method, rng=0, noise_cov=cov) with ( pytest.warns(RuntimeWarning, match="high-pass filtered"), pytest.warns(RuntimeWarning, match="No average EEG.*"), @@ -191,7 +189,7 @@ def test_warnings(): epochs = Epochs( raw, events=events, baseline=None, preload=True, on_outside="ignore" ) - ica = ICA(n_components=2, max_iter=1, method="infomax", random_state=0) + ica = ICA(n_components=2, max_iter=1, method="infomax", rng=0) # not high-passed with epochs.info._unlock(): @@ -292,6 +290,36 @@ def test_ica_max_iter_(method, max_iter_default): ICA(max_iter=1.0) +def test_ica_rng_transition(): + """Test the transition from random_state to rng.""" + with pytest.raises(TypeError, match="only one"): + _ICA(random_state=0, rng=0) + + info = create_info(["Fz", "Cz", "Pz"], 100.0, "eeg") + with info._unlock(): + info["highpass"] = 1.0 + raw = RawArray(np.random.default_rng(0).standard_normal((3, 200)), info) + unmixings = [] + for kwargs in ( + dict(random_state=0), + dict(random_state=check_random_state(0)), + dict(rng=0), + ): + ica = _ICA( + n_components=2, + method="fastica", + max_iter=1000, + **kwargs, + ) + with _record_warnings(): # ICA does not necessarily converge + ica.fit(raw) + unmixings.append(ica.unmixing_matrix_) + assert_array_equal(unmixings[0], unmixings[1]) + # at the ICA/sklearn boundary an integer ``rng`` seed is forwarded verbatim, + # so it matches the same integer passed to the legacy parameter + assert_array_equal(unmixings[0], unmixings[2]) + + @pytest.mark.parametrize("method", ["infomax", "fastica", "picard"]) def test_ica_n_iter_(method, tmp_path): """Test that ICA.n_iter_ is set after fitting.""" @@ -300,9 +328,7 @@ def test_ica_n_iter_(method, tmp_path): raw = read_raw_fif(raw_fname).crop(0.5, stop).load_data() n_components = 3 max_iter = 1 - ica = ICA( - n_components=n_components, max_iter=max_iter, method=method, random_state=0 - ) + ica = ICA(n_components=n_components, max_iter=max_iter, method=method, rng=0) if method == "infomax": ica.fit(raw) @@ -1399,7 +1425,7 @@ def test_n_components_none(method, tmp_path): random_state = 12345 output_fname = tmp_path / "test_ica-ica.fif" - ica = ICA(method=method, n_components=n_components, random_state=random_state) + ica = ICA(method=method, n_components=n_components, rng=random_state) with _record_warnings(): ica.fit(epochs) _assert_ica_attributes(ica) @@ -1774,13 +1800,13 @@ def test_ica_rejects_nonfinite(): # Case 1: NaN raw = RawArray(data.copy(), info) raw._data[0, 25] = np.nan - ica = ICA(n_components=2, random_state=0, method="fastica", max_iter="auto") + ica = ICA(n_components=2, rng=0, method="fastica", max_iter="auto") with pytest.raises(ValueError, match=r"Input data contains non-finite values"): ica.fit(raw) # Case 2: Inf raw = RawArray(data.copy(), info) raw._data[1, 50] = np.inf - ica = ICA(n_components=2, random_state=0, method="fastica", max_iter="auto") + ica = ICA(n_components=2, rng=0, method="fastica", max_iter="auto") with pytest.raises(ValueError, match=r"Input data contains non-finite values"): ica.fit(raw) diff --git a/mne/preprocessing/tests/test_infomax.py b/mne/preprocessing/tests/test_infomax.py index 6ac98b85136..075a99f2d47 100644 --- a/mne/preprocessing/tests/test_infomax.py +++ b/mne/preprocessing/tests/test_infomax.py @@ -6,11 +6,11 @@ import numpy as np import pytest -from numpy.testing import assert_almost_equal +from numpy.testing import assert_almost_equal, assert_array_equal from scipy import stats from mne.preprocessing.infomax_ import infomax -from mne.utils import pinv +from mne.utils import check_random_state, pinv pytest.importorskip("sklearn") @@ -51,7 +51,7 @@ def test_infomax_blowup(): center_and_norm(m) X = _get_pca(0).fit_transform(m.T) - k_ = infomax(X, extended=True, l_rate=0.1, random_state=0) + k_ = infomax(X, extended=True, l_rate=0.1, rng=0) s_ = np.dot(k_, X.T) center_and_norm(s_) @@ -93,7 +93,7 @@ def test_infomax_simple(): algos = [True, False] for algo in algos: X = _get_pca(0).fit_transform(m.T) - k_ = infomax(X, extended=algo, random_state=0) + k_ = infomax(X, extended=algo, rng=0) s_ = np.dot(k_, X.T) center_and_norm(s_) @@ -120,8 +120,8 @@ def test_infomax_weights_ini(): X = rng.random((3, 100)) weights = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtype=np.float64) - w1 = infomax(X, max_iter=0, weights=weights, extended=True, random_state=0) - w2 = infomax(X, max_iter=0, weights=weights, extended=False, random_state=0) + w1 = infomax(X, max_iter=0, weights=weights, extended=True, rng=0) + w2 = infomax(X, max_iter=0, weights=weights, extended=False, rng=0) assert_almost_equal(w1, weights) assert_almost_equal(w2, weights) @@ -153,7 +153,7 @@ def test_non_square_infomax(): m = m.T m = _get_pca(0).fit_transform(m) # we need extended since input signals are sub-gaussian - unmixing_ = infomax(m, random_state=0, extended=True) + unmixing_ = infomax(m, rng=0, extended=True) s_ = np.dot(unmixing_, m.T) # Check that the mixing model described in the docstring holds: mixing_ = pinv(unmixing_.T) @@ -181,9 +181,7 @@ def test_infomax_n_iter(return_n_iter): rng = np.random.default_rng(0) X = rng.random((3, 100)) max_iter = 1 - r = infomax( - X, max_iter=max_iter, extended=True, return_n_iter=return_n_iter, random_state=0 - ) + r = infomax(X, max_iter=max_iter, extended=True, return_n_iter=return_n_iter, rng=0) if return_n_iter: assert isinstance(r, tuple) @@ -192,6 +190,23 @@ def test_infomax_n_iter(return_n_iter): assert isinstance(r, np.ndarray) +def test_infomax_legacy_rng_nested(): + """Test legacy RNGs survive Infomax's nested permutation path.""" + X = np.random.default_rng(0).standard_normal((20, 2)) + results = [] + for random_state in (0, check_random_state(0)): + results.append( + infomax( + X, + block=5, + extended=False, + max_iter=1, + random_state=random_state, + ) + ) + assert_array_equal(results[0], results[1]) + + def _get_pca(rng=None): from sklearn.decomposition import PCA diff --git a/mne/preprocessing/tests/test_xdawn.py b/mne/preprocessing/tests/test_xdawn.py index 86bae4f9997..0fd83247706 100644 --- a/mne/preprocessing/tests/test_xdawn.py +++ b/mne/preprocessing/tests/test_xdawn.py @@ -372,13 +372,13 @@ def test_xdawn_decoding_performance(): Xdawn(n_components=n_xdawn_comps), Vectorizer(), MinMaxScaler(), - LogisticRegression(solver="liblinear"), + LogisticRegression(solver="liblinear", random_state=0), ) xdawn_trans_pipe = make_pipeline( XdawnTransformer(n_components=n_xdawn_comps), Vectorizer(), MinMaxScaler(), - LogisticRegression(solver="liblinear"), + LogisticRegression(solver="liblinear", random_state=0), ) cv = KFold(n_splits=3, shuffle=False) diff --git a/mne/report/tests/test_report.py b/mne/report/tests/test_report.py index 4220c8aeea3..96af2bc37e5 100644 --- a/mne/report/tests/test_report.py +++ b/mne/report/tests/test_report.py @@ -957,9 +957,7 @@ def test_manual_report_2d(tmp_path, invisible_fig): evoked = evokeds[0].pick("eeg").decimate(10, verbose="error") with pytest.warns(ConvergenceWarning, match="did not converge"): - ica = ICA(n_components=3, max_iter=1, random_state=42).fit( - inst=raw.copy().crop(tmax=1) - ) + ica = ICA(n_components=3, max_iter=1, rng=42).fit(inst=raw.copy().crop(tmax=1)) ica_ecg_scores = ica_eog_scores = np.array([3, 0, 0]) ica_ecg_evoked = ica_eog_evoked = epochs_without_metadata.average() diff --git a/mne/simulation/evoked.py b/mne/simulation/evoked.py index 9805e520496..d1229c23a7b 100644 --- a/mne/simulation/evoked.py +++ b/mne/simulation/evoked.py @@ -13,9 +13,17 @@ from ..evoked import Evoked from ..forward import apply_forward from ..io import BaseRaw -from ..utils import _check_preload, _validate_type, check_random_state, logger, verbose +from ..utils import ( + _check_preload, + _legacy_rng, + _validate_type, + check_random_state, + logger, + verbose, +) +@_legacy_rng("random_state") @verbose def simulate_evoked( fwd, @@ -24,9 +32,11 @@ def simulate_evoked( cov=None, nave=30, iir_filter=None, - random_state=None, use_cps=True, verbose=None, + *, + rng=None, + random_state=None, ): """Generate noisy evoked data. @@ -51,11 +61,12 @@ def simulate_evoked( .. versionadded:: 0.15.0 iir_filter : None | array IIR filter coefficients (denominator) e.g. [1, -1, 0.2]. - %(random_state)s %(use_cps)s .. versionadded:: 0.15 %(verbose)s + %(rng)s + %(random_state_rng)s Returns ------- @@ -84,7 +95,7 @@ def simulate_evoked( return evoked if nave < np.inf: - noise = _simulate_noise_evoked(evoked, cov, iir_filter, random_state) + noise = _simulate_noise_evoked(evoked, cov, iir_filter, rng) evoked.data += noise.data / math.sqrt(nave) evoked.nave = np.int64(nave) if cov.get("projs", None): @@ -92,14 +103,15 @@ def simulate_evoked( return evoked -def _simulate_noise_evoked(evoked, cov, iir_filter, random_state): +def _simulate_noise_evoked(evoked, cov, iir_filter, rng): noise = evoked.copy() noise.data[:] = 0 - return _add_noise(noise, cov, iir_filter, random_state, allow_subselection=False) + return _add_noise(noise, cov, iir_filter, rng, allow_subselection=False) +@_legacy_rng("random_state") @verbose -def add_noise(inst, cov, iir_filter=None, random_state=None, verbose=None): +def add_noise(inst, cov, iir_filter=None, verbose=None, *, rng=None, random_state=None): """Create noise as a multivariate Gaussian. The spatial covariance of the noise is given from the cov matrix. @@ -112,8 +124,9 @@ def add_noise(inst, cov, iir_filter=None, random_state=None, verbose=None): The noise covariance. iir_filter : None | array-like IIR filter coefficients (denominator). - %(random_state)s %(verbose)s + %(rng)s + %(random_state_rng)s Returns ------- @@ -130,10 +143,13 @@ def add_noise(inst, cov, iir_filter=None, random_state=None, verbose=None): .. versionadded:: 0.18.0 """ # We always allow subselection here - return _add_noise(inst, cov, iir_filter, random_state) + legacy_seed = random_state if isinstance(random_state, int | np.integer) else None + return _add_noise(inst, cov, iir_filter, rng, legacy_seed=legacy_seed) -def _add_noise(inst, cov, iir_filter, random_state, allow_subselection=True): +def _add_noise( + inst, cov, iir_filter, rng, allow_subselection=True, *, legacy_seed=None +): """Add noise, possibly with channel subselection.""" _validate_type(cov, Covariance, "cov") _validate_type( @@ -162,17 +178,17 @@ def _add_noise(inst, cov, iir_filter, random_state, allow_subselection=True): gen_picks = np.arange(info["nchan"]) for epoch in data: + # An integer passed to the deprecated parameter historically restarted + # the same stream for each epoch. ``rng`` intentionally advances. + this_rng = check_random_state(legacy_seed) if legacy_seed is not None else rng epoch[picks] += _generate_noise( - info, cov, iir_filter, random_state, epoch.shape[1], picks=gen_picks + info, cov, iir_filter, this_rng, epoch.shape[1], picks=gen_picks )[0] return inst -def _generate_noise( - info, cov, iir_filter, random_state, n_samples, zi=None, picks=None -): +def _generate_noise(info, cov, iir_filter, rng, n_samples, zi=None, picks=None): """Create spatially colored and temporally IIR-filtered noise.""" - rng = check_random_state(random_state) _, _, colorer = compute_whitener( cov, info, pca=True, return_colorer=True, picks=picks, verbose=False ) diff --git a/mne/simulation/raw.py b/mne/simulation/raw.py index ccb32f13380..5d5ddbac0f3 100644 --- a/mne/simulation/raw.py +++ b/mne/simulation/raw.py @@ -44,10 +44,10 @@ from ..transforms import Transform, _get_trans, transform_surface_to from ..utils import ( _check_preload, + _legacy_rng, _pl, _validate_type, _verbose_safe_false, - check_random_state, logger, verbose, ) @@ -387,9 +387,17 @@ def simulate_raw( return raw +@_legacy_rng("random_state") @verbose def add_eog( - raw, head_pos=None, interp="cos2", n_jobs=None, random_state=None, verbose=None + raw, + head_pos=None, + interp="cos2", + n_jobs=None, + verbose=None, + *, + rng=None, + random_state=None, ): """Add blink noise to raw data. @@ -400,10 +408,11 @@ def add_eog( %(head_pos)s %(interp)s %(n_jobs)s - %(random_state)s + %(verbose)s + %(rng)s + %(random_state_rng)s The random generator state used for blink, ECG, and sensor noise randomization. - %(verbose)s Returns ------- @@ -439,12 +448,20 @@ def add_eog( ---------- .. footbibliography:: """ - return _add_exg(raw, "blink", head_pos, interp, n_jobs, random_state) + return _add_exg(raw, "blink", head_pos, interp, n_jobs, rng) +@_legacy_rng("random_state") @verbose def add_ecg( - raw, head_pos=None, interp="cos2", n_jobs=None, random_state=None, verbose=None + raw, + head_pos=None, + interp="cos2", + n_jobs=None, + verbose=None, + *, + rng=None, + random_state=None, ): """Add ECG noise to raw data. @@ -455,10 +472,11 @@ def add_ecg( %(head_pos)s %(interp)s %(n_jobs)s - %(random_state)s + %(verbose)s + %(rng)s + %(random_state_rng)s The random generator state used for blink, ECG, and sensor noise randomization. - %(verbose)s Returns ------- @@ -492,14 +510,13 @@ def add_ecg( .. versionadded:: 0.18 """ - return _add_exg(raw, "ecg", head_pos, interp, n_jobs, random_state) + return _add_exg(raw, "ecg", head_pos, interp, n_jobs, rng) -def _add_exg(raw, kind, head_pos, interp, n_jobs, random_state): +def _add_exg(raw, kind, head_pos, interp, n_jobs, rng): assert isinstance(kind, str) and kind in ("ecg", "blink") _validate_type(raw, BaseRaw, "raw") _check_preload(raw, f"Adding {kind} noise ") - rng = check_random_state(random_state) info, times, first_samp = raw.info, raw.times, raw.first_samp data = raw._data meg_picks = pick_types(info, meg=True, eeg=False, exclude=()) diff --git a/mne/simulation/source.py b/mne/simulation/source.py index a1353d9612d..25d1e90f749 100644 --- a/mne/simulation/source.py +++ b/mne/simulation/source.py @@ -13,22 +13,25 @@ _check_option, _ensure_events, _ensure_int, + _legacy_rng, _validate_type, - check_random_state, fill_doc, warn, ) +@_legacy_rng("random_state") @fill_doc def select_source_in_label( src, label, - random_state=None, location="random", subject=None, subjects_dir=None, surf="sphere", + *, + rng=None, + random_state=None, ): """Select source positions using a label. @@ -38,7 +41,6 @@ def select_source_in_label( The source space. label : Label The label. - %(random_state)s location : str The label location to choose. Can be 'random' (default) or 'center' to use :func:`mne.Label.center_of_mass` (restricting to vertices @@ -61,6 +63,8 @@ def select_source_in_label( with cortical folding. .. versionadded:: 0.13 + %(rng)s + %(random_state_rng)s Returns ------- @@ -69,11 +73,22 @@ def select_source_in_label( rh_vertno : list Selected source coefficients on the right hemisphere. """ + return _select_source_in_label( + src, + label, + location=location, + subject=subject, + subjects_dir=subjects_dir, + surf=surf, + rng=rng, + ) + + +def _select_source_in_label(src, label, *, location, subject, subjects_dir, surf, rng): + """Select a source in a label using an already normalized RNG.""" lh_vertno = list() rh_vertno = list() _check_option("location", location, ["random", "center"]) - - rng = check_random_state(random_state) if label.hemi == "lh": vertno = lh_vertno hemi_idx = 0 @@ -91,6 +106,7 @@ def select_source_in_label( return lh_vertno, rh_vertno +@_legacy_rng("random_state") @fill_doc def simulate_sparse_stc( src, @@ -98,11 +114,13 @@ def simulate_sparse_stc( times, data_fun=lambda t: 1e-7 * np.sin(20 * np.pi * t), labels=None, - random_state=None, location="random", subject=None, subjects_dir=None, surf="sphere", + *, + rng=None, + random_state=None, ): """Generate sparse (n_dipoles) sources time courses from data_fun. @@ -126,7 +144,6 @@ def simulate_sparse_stc( the same length containing the time courses. labels : None | list of Label The labels. The default is None, otherwise its size must be n_dipoles. - %(random_state)s location : str The label location to choose. Can be ``'random'`` (default) or ``'center'`` to use :func:`mne.Label.center_of_mass`. Note that for @@ -148,6 +165,8 @@ def simulate_sparse_stc( with cortical folding. .. versionadded:: 0.13 + %(rng)s + %(random_state_rng)s Returns ------- @@ -164,7 +183,6 @@ def simulate_sparse_stc( ----- .. versionadded:: 0.10.0 """ - rng = check_random_state(random_state) src = _ensure_src(src, verbose=False) subject_src = src._subject if subject is None: @@ -205,8 +223,14 @@ def simulate_sparse_stc( lh_data = [np.empty((0, data.shape[1]))] rh_data = [np.empty((0, data.shape[1]))] for i, label in enumerate(labels): - lh_vertno, rh_vertno = select_source_in_label( - src, label, rng, location, subject, subjects_dir, surf + lh_vertno, rh_vertno = _select_source_in_label( + src, + label, + location=location, + subject=subject, + subjects_dir=subjects_dir, + surf=surf, + rng=rng, ) vertno[0] += lh_vertno vertno[1] += rh_vertno diff --git a/mne/simulation/tests/test_evoked.py b/mne/simulation/tests/test_evoked.py index f8c64ef23a2..f805658fea3 100644 --- a/mne/simulation/tests/test_evoked.py +++ b/mne/simulation/tests/test_evoked.py @@ -14,10 +14,12 @@ ) from mne import ( + Covariance, EpochsArray, compute_covariance, compute_raw_covariance, convert_forward_solution, + create_info, pick_channels_cov, pick_types, pick_types_forward, @@ -60,7 +62,7 @@ def test_simulate_evoked(): times = np.linspace(tmin, tmin + n_samples * tstep, n_samples) # Generate times series for 2 dipoles - stc = simulate_sparse_stc(fwd["src"], n_dipoles=2, times=times, random_state=42) + stc = simulate_sparse_stc(fwd["src"], n_dipoles=2, times=times, rng=42) # Generate noisy evoked data iir_filter = [1, -0.9] @@ -71,7 +73,7 @@ def test_simulate_evoked(): cov, iir_filter=iir_filter, nave=nave, - random_state=0, + rng=0, ) assert_array_almost_equal(evoked.times, stc.times) assert len(evoked.data) == len(fwd["sol"]["data"]) @@ -128,7 +130,7 @@ def test_add_noise(): evoked = epochs.average(picks=np.arange(len(raw.ch_names))) for inst in (raw, epochs, evoked): with catch_logging() as log: - add_noise(inst, cov, random_state=rng, verbose=True) + add_noise(inst, cov, rng=rng, verbose=True) log = log.getvalue() want = "to {0}/{1} channels ({0}".format(len(cov["names"]), len(raw.ch_names)) assert want in log @@ -143,6 +145,28 @@ def test_add_noise(): r = np.corrcoef(cov["data"].ravel(), cov_new["data"].ravel())[0, 1] assert r > 0.99 + info = create_info(["EEG 001"], 100.0, "eeg") + small_cov = Covariance(np.ones(1), info["ch_names"], [], [], 1) + legacy = EpochsArray(np.zeros((2, 1, 5)), info, verbose=False) + add_noise(legacy, small_cov, random_state=0) + want = np.array( + [ + 1.764052345967664, + 0.4001572083672233, + 0.9787379841057392, + 2.240893199201458, + 1.8675579901499675, + ] + ) + assert_array_equal(legacy.get_data(copy=False)[0, 0], want) + assert_array_equal(legacy.get_data(copy=False)[1, 0], want) + + modern = EpochsArray(np.zeros((2, 1, 5)), info, verbose=False) + add_noise(modern, small_cov, rng=0) + assert not np.array_equal( + modern.get_data(copy=False)[0], modern.get_data(copy=False)[1] + ) + def test_rank_deficiency(): """Test adding noise from M/EEG float32 (I/O) cov with projectors.""" @@ -160,7 +184,7 @@ def test_rank_deficiency(): cov = regularize(cov, evoked.info, rank=None) cov = pick_channels_cov(cov, evoked.ch_names) evoked.data[:] = 0 - add_noise(evoked, cov, random_state=0) + add_noise(evoked, cov, rng=0) cov_new = compute_covariance( EpochsArray(evoked.data[np.newaxis], evoked.info), verbose="error" ) @@ -182,7 +206,7 @@ def test_order(): # MEG then EEG assert (eeg_picks > meg_picks.max()).all() times = np.arange(10) / 1000.0 - stc = simulate_sparse_stc(fwd["src"], 1, times=times, random_state=0) + stc = simulate_sparse_stc(fwd["src"], 1, times=times, rng=0) evoked_sim = simulate_evoked(fwd, stc, evoked.info, nave=np.inf) reorder = np.concatenate([eeg_picks, meg_picks]) evoked.reorder_channels([evoked.ch_names[pick] for pick in reorder]) diff --git a/mne/simulation/tests/test_metrics.py b/mne/simulation/tests/test_metrics.py index c2deca7816d..3b76264cc6f 100644 --- a/mne/simulation/tests/test_metrics.py +++ b/mne/simulation/tests/test_metrics.py @@ -21,8 +21,8 @@ def test_metrics(): src = read_source_spaces(src_fname) times = np.arange(600) / 1000.0 rng = np.random.default_rng(42) - stc1 = simulate_sparse_stc(src, n_dipoles=2, times=times, random_state=rng) - stc2 = simulate_sparse_stc(src, n_dipoles=2, times=times, random_state=rng) + stc1 = simulate_sparse_stc(src, n_dipoles=2, times=times, rng=rng) + stc2 = simulate_sparse_stc(src, n_dipoles=2, times=times, rng=rng) E1_rms = source_estimate_quantification(stc1, stc1, metric="rms") E2_rms = source_estimate_quantification(stc2, stc2, metric="rms") E1_cos = source_estimate_quantification(stc1, stc1, metric="cosine") diff --git a/mne/simulation/tests/test_raw.py b/mne/simulation/tests/test_raw.py index 4b41632b558..221d7ec1452 100644 --- a/mne/simulation/tests/test_raw.py +++ b/mne/simulation/tests/test_raw.py @@ -196,7 +196,7 @@ def _make_stc(raw, src): tstep = 1.0 / sfreq n_samples = len(raw.times) // 10 times = np.arange(0, n_samples) * tstep - stc = simulate_sparse_stc(src, 10, times, random_state=seed) + stc = simulate_sparse_stc(src, 10, times, rng=seed) return stc @@ -280,9 +280,9 @@ def test_simulate_raw_sphere(raw_data, tmp_path): raw.copy().pick(["meg", "eeg"]).info, stc, trans, src, sphere ) for this_raw in (raw_sim_meg, raw_sim_eeg, raw_sim_meeg): - add_eog(this_raw, random_state=seed) + add_eog(this_raw, rng=seed) for this_raw in (raw_sim_meg, raw_sim_meeg): - add_ecg(this_raw, random_state=seed) + add_ecg(this_raw, rng=seed) with pytest.raises(RuntimeError, match="only add ECG artifacts if MEG"): add_ecg(raw_sim_eeg) assert_allclose( @@ -590,24 +590,24 @@ def test_simulation_cascade(): # Calculate independent signal additions raw_eog = raw_null.copy() - add_eog(raw_eog, random_state=0) + add_eog(raw_eog, rng=0) raw_ecg = raw_null.copy() - add_ecg(raw_ecg, random_state=0) + add_ecg(raw_ecg, rng=0) raw_noise = raw_null.copy() cov = make_ad_hoc_cov(raw_null.info) - add_noise(raw_noise, cov, random_state=0) + add_noise(raw_noise, cov, rng=0) raw_chpi = raw_null.copy() add_chpi(raw_chpi) # Calculate Cascading signal additions raw_cascade = raw_null.copy() - add_eog(raw_cascade, random_state=0) - add_ecg(raw_cascade, random_state=0) + add_eog(raw_cascade, rng=0) + add_ecg(raw_cascade, rng=0) add_chpi(raw_cascade) - add_noise(raw_cascade, cov, random_state=0) + add_noise(raw_cascade, cov, rng=0) cascade_data = raw_cascade.get_data() serial_data = 0.0 diff --git a/mne/simulation/tests/test_source.py b/mne/simulation/tests/test_source.py index eb5eb1bd983..a2f35b2221e 100644 --- a/mne/simulation/tests/test_source.py +++ b/mne/simulation/tests/test_source.py @@ -155,7 +155,7 @@ def test_simulate_sparse_stc(_get_fwd_labels): len(mylabels), times, labels=mylabels, - random_state=random_state, + rng=random_state, location=location, subjects_dir=subjects_dir, ) @@ -170,7 +170,7 @@ def test_simulate_sparse_stc(_get_fwd_labels): len(mylabels), times, labels=mylabels, - random_state=random_state, + rng=random_state, location=location, subjects_dir=subjects_dir, ) @@ -218,7 +218,7 @@ def test_simulate_sparse_stc(_get_fwd_labels): len(mylabels) + 1, times, labels=mylabels, - random_state=random_state, + rng=random_state, location=location, subjects_dir=subjects_dir, ) @@ -291,7 +291,7 @@ def test_simulate_sparse_stc_single_hemi(_get_fwd_labels): len(labels_single_hemi), times, labels=labels_single_hemi, - random_state=0, + rng=0, ) assert stc_1.data.shape[0] == len(labels_single_hemi) @@ -303,7 +303,7 @@ def test_simulate_sparse_stc_single_hemi(_get_fwd_labels): len(labels_single_hemi), times, labels=labels_single_hemi, - random_state=0, + rng=0, ) assert_array_equal(stc_1.lh_vertno, stc_2.lh_vertno) diff --git a/mne/stats/cluster_level.py b/mne/stats/cluster_level.py index 0de6e009504..c98989cbdb8 100644 --- a/mne/stats/cluster_level.py +++ b/mne/stats/cluster_level.py @@ -12,9 +12,9 @@ from ..utils import ( ProgressBar, _check_option, + _legacy_rng, _pl, _validate_type, - check_random_state, logger, split_list, verbose, @@ -821,7 +821,7 @@ def _permutation_cluster_test( stat_fun, adjacency, n_jobs, - seed, + rng, max_step, exclude, step_down_p, @@ -952,12 +952,10 @@ def _permutation_cluster_test( if out_type == "indices": clusters = _cluster_mask_to_indices(clusters, t_obs.shape) - # convert our seed to orders + # Convert the RNG state to permutation orders. # check to see if we can do an exact test # (for a two-tailed test, we can exploit symmetry to just do half) extra = "" - rng = check_random_state(seed) - del seed if len(X) == 1: # 1-sample test do_perm_func = _do_1samp_permutations X_full = X[0] @@ -1102,6 +1100,7 @@ def _check_fun(X, stat_fun, threshold, tail=0, kind="within"): return stat_fun, threshold +@_legacy_rng("seed") @verbose def permutation_cluster_test( X, @@ -1111,7 +1110,6 @@ def permutation_cluster_test( stat_fun=None, adjacency=None, n_jobs=None, - seed=None, max_step=1, exclude=None, step_down_p=0, @@ -1120,6 +1118,9 @@ def permutation_cluster_test( check_disjoint=False, buffer_size=1000, verbose=None, + *, + rng=None, + seed=None, ): """Cluster-level statistical permutation test. @@ -1151,7 +1152,6 @@ def permutation_cluster_test( %(stat_fun_clust_f)s %(adjacency_clust_n)s %(n_jobs)s - %(seed)s %(max_step_clust)s %(exclude_clust)s %(step_down_p_clust)s @@ -1160,6 +1160,8 @@ def permutation_cluster_test( %(check_disjoint_clust)s %(buffer_size_clust)s %(verbose)s + %(rng)s + %(seed_rng)s Returns ------- @@ -1189,7 +1191,7 @@ def permutation_cluster_test( stat_fun=stat_fun, adjacency=adjacency, n_jobs=n_jobs, - seed=seed, + rng=rng, max_step=max_step, exclude=exclude, step_down_p=step_down_p, @@ -1200,6 +1202,7 @@ def permutation_cluster_test( ) +@_legacy_rng("seed") @verbose def permutation_cluster_1samp_test( X, @@ -1209,7 +1212,6 @@ def permutation_cluster_1samp_test( stat_fun=None, adjacency=None, n_jobs=None, - seed=None, max_step=1, exclude=None, step_down_p=0, @@ -1218,6 +1220,9 @@ def permutation_cluster_1samp_test( check_disjoint=False, buffer_size=1000, verbose=None, + *, + rng=None, + seed=None, ): """Non-parametric cluster-level paired t-test. @@ -1238,7 +1243,6 @@ def permutation_cluster_1samp_test( %(stat_fun_clust_t)s %(adjacency_clust_1)s %(n_jobs)s - %(seed)s %(max_step_clust)s %(exclude_clust)s %(step_down_p_clust)s @@ -1247,6 +1251,8 @@ def permutation_cluster_1samp_test( %(check_disjoint_clust)s %(buffer_size_clust)s %(verbose)s + %(rng)s + %(seed_rng)s Returns ------- @@ -1277,9 +1283,9 @@ def permutation_cluster_1samp_test( %(threshold_clust_t_notes)s If ``n_permutations`` exceeds the maximum number of possible permutations - given the number of observations, then ``n_permutations`` and ``seed`` - will be ignored since an exact test (full permutation test) will be - performed (this is the case when + given the number of observations, then ``n_permutations``, ``seed``, and + ``rng`` will be ignored since an exact test (full permutation test) will + be performed (this is the case when ``n_permutations >= 2 ** (n_observations - (tail == 0))``). If no initial clusters are found because all points in the true @@ -1299,7 +1305,7 @@ def permutation_cluster_1samp_test( stat_fun=stat_fun, adjacency=adjacency, n_jobs=n_jobs, - seed=seed, + rng=rng, max_step=max_step, exclude=exclude, step_down_p=step_down_p, @@ -1310,6 +1316,7 @@ def permutation_cluster_1samp_test( ) +@_legacy_rng("seed") @verbose def spatio_temporal_cluster_1samp_test( X, @@ -1319,7 +1326,6 @@ def spatio_temporal_cluster_1samp_test( stat_fun=None, adjacency=None, n_jobs=None, - seed=None, max_step=1, spatial_exclude=None, step_down_p=0, @@ -1328,6 +1334,9 @@ def spatio_temporal_cluster_1samp_test( check_disjoint=False, buffer_size=1000, verbose=None, + *, + rng=None, + seed=None, ): """Non-parametric cluster-level paired t-test for spatio-temporal data. @@ -1351,7 +1360,6 @@ def spatio_temporal_cluster_1samp_test( %(stat_fun_clust_t)s %(adjacency_clust_st1)s %(n_jobs)s - %(seed)s %(max_step_clust)s spatial_exclude : list of int or None List of spatial indices to exclude from clustering. @@ -1361,6 +1369,8 @@ def spatio_temporal_cluster_1samp_test( %(check_disjoint_clust)s %(buffer_size_clust)s %(verbose)s + %(rng)s + %(seed_rng)s Returns ------- @@ -1396,7 +1406,7 @@ def spatio_temporal_cluster_1samp_test( n_permutations=n_permutations, adjacency=adjacency, n_jobs=n_jobs, - seed=seed, + rng=rng, max_step=max_step, exclude=exclude, step_down_p=step_down_p, @@ -1407,6 +1417,7 @@ def spatio_temporal_cluster_1samp_test( ) +@_legacy_rng("seed") @verbose def spatio_temporal_cluster_test( X, @@ -1416,7 +1427,6 @@ def spatio_temporal_cluster_test( stat_fun=None, adjacency=None, n_jobs=None, - seed=None, max_step=1, spatial_exclude=None, step_down_p=0, @@ -1425,6 +1435,9 @@ def spatio_temporal_cluster_test( check_disjoint=False, buffer_size=1000, verbose=None, + *, + rng=None, + seed=None, ): """Non-parametric cluster-level test for spatio-temporal data. @@ -1450,7 +1463,6 @@ def spatio_temporal_cluster_test( %(stat_fun_clust_f)s %(adjacency_clust_stn)s %(n_jobs)s - %(seed)s %(max_step_clust)s spatial_exclude : list of int or None List of spatial indices to exclude from clustering. @@ -1460,6 +1472,8 @@ def spatio_temporal_cluster_test( %(check_disjoint_clust)s %(buffer_size_clust)s %(verbose)s + %(rng)s + %(seed_rng)s Returns ------- @@ -1495,7 +1509,7 @@ def spatio_temporal_cluster_test( n_permutations=n_permutations, adjacency=adjacency, n_jobs=n_jobs, - seed=seed, + rng=rng, max_step=max_step, exclude=exclude, step_down_p=step_down_p, diff --git a/mne/stats/permutations.py b/mne/stats/permutations.py index 903e1b56881..33ba199cac7 100644 --- a/mne/stats/permutations.py +++ b/mne/stats/permutations.py @@ -9,7 +9,13 @@ import numpy as np from ..parallel import parallel_func -from ..utils import _check_if_nan, check_random_state, logger, verbose +from ..utils import ( + _check_if_nan, + _legacy_rng, + fill_doc, + logger, + verbose, +) def _max_stat(X, X2, perms, dof_scaling): @@ -21,9 +27,17 @@ def _max_stat(X, X2, perms, dof_scaling): return max_abs +@_legacy_rng("seed") @verbose def permutation_t_test( - X, n_permutations=10000, tail=0, n_jobs=None, seed=None, verbose=None + X, + n_permutations=10000, + tail=0, + n_jobs=None, + verbose=None, + *, + rng=None, + seed=None, ): """One sample/paired sample permutation test based on a t-statistic. @@ -52,8 +66,9 @@ def permutation_t_test( than 0 (two tailed test). If tail is -1, the alternative hypothesis is that the mean of the data is less than 0 (lower tailed test). %(n_jobs)s - %(seed)s %(verbose)s + %(rng)s + %(seed_rng)s Returns ------- @@ -68,8 +83,8 @@ def permutation_t_test( Notes ----- If ``n_permutations >= 2 ** (n_samples - (tail == 0))``, - ``n_permutations`` and ``seed`` will be ignored since an exact test - (full permutation test) will be performed. + ``n_permutations``, ``seed``, and ``rng`` will be ignored since an exact + test (full permutation test) will be performed. References ---------- @@ -84,7 +99,6 @@ def permutation_t_test( dof_scaling = sqrt(n_samples / (n_samples - 1.0)) std0 = np.sqrt(X2 - mu0**2) * dof_scaling # get std with var splitting T_obs = np.mean(X, axis=0) / (std0 / sqrt(n_samples)) - rng = check_random_state(seed) orders, _, extra = _get_1samp_orders(n_samples, n_permutations, tail, rng) perms = 2 * np.array(orders) - 1 # from 0, 1 -> 1, -1 logger.info(f"Permuting {len(orders)} times{extra}...") @@ -105,8 +119,16 @@ def permutation_t_test( return T_obs, p_values, H0 +@_legacy_rng("random_state") +@fill_doc def bootstrap_confidence_interval( - arr, ci=0.95, n_bootstraps=2000, stat_fun="mean", random_state=None + arr, + ci=0.95, + n_bootstraps=2000, + stat_fun="mean", + *, + rng=None, + random_state=None, ): """Get confidence intervals from non-parametric bootstrap. @@ -120,8 +142,8 @@ def bootstrap_confidence_interval( Number of bootstraps. stat_fun : str | callable Can be "mean", "median", or a callable operating along ``axis=0``. - random_state : int | float | array_like | None - The seed at which to initialize the bootstrap. + %(rng)s + %(random_state_rng)s Returns ------- @@ -143,7 +165,6 @@ def stat_fun(x): raise ValueError("stat_fun must be 'mean', 'median' or callable.") n_trials = arr.shape[0] indices = np.arange(n_trials, dtype=int) # BCA would be cool to have too - rng = check_random_state(random_state) boot_indices = rng.choice(indices, replace=True, size=(n_bootstraps, len(indices))) stat = np.array([stat_fun(arr[inds]) for inds in boot_indices]) ci = (((1 - ci) / 2) * 100, (1 - ((1 - ci) / 2)) * 100) @@ -151,11 +172,11 @@ def stat_fun(x): return np.array([ci_low, ci_up]) -def _ci(arr, ci=0.95, method="bootstrap", n_bootstraps=2000, random_state=None): +def _ci(arr, ci=0.95, method="bootstrap", n_bootstraps=2000, rng=None): """Calculate confidence interval. Aux function for plot_compare_evokeds.""" if method == "bootstrap": return bootstrap_confidence_interval( - arr, ci=ci, n_bootstraps=n_bootstraps, random_state=random_state + arr, ci=ci, n_bootstraps=n_bootstraps, rng=rng ) else: from .parametric import _parametric_ci diff --git a/mne/stats/tests/test_cluster_level.py b/mne/stats/tests/test_cluster_level.py index 071861d85ff..e306de32c81 100644 --- a/mne/stats/tests/test_cluster_level.py +++ b/mne/stats/tests/test_cluster_level.py @@ -69,7 +69,7 @@ def test_thresholds(numba_conditional): with catch_logging() as log: with pytest.warns(RuntimeWarning, match="threshold is only valid"): out = permutation_cluster_1samp_test( - X, stat_fun=my_fun, seed=0, verbose=True, out_type="mask" + X, stat_fun=my_fun, rng=0, verbose=True, out_type="mask" ) log = log.getvalue() assert str(want_thresh)[:6] in log @@ -86,12 +86,12 @@ def test_thresholds(numba_conditional): with catch_logging() as log: with pytest.warns(RuntimeWarning, match="threshold is only valid"): out = permutation_cluster_test( - X, tail=1, stat_fun=my_fun, seed=0, verbose=True, out_type="mask" + X, tail=1, stat_fun=my_fun, rng=0, verbose=True, out_type="mask" ) log = log.getvalue() assert str(want_thresh)[:6] in log assert len(out[1]) == 1 # 1 cluster - assert_allclose(out[2], 0.031250, atol=1e-6) + assert_allclose(out[2], 0.03515625, atol=1e-6) with pytest.warns(RuntimeWarning, match='Ignoring argument "tail"'): permutation_cluster_test(X, tail=0, out_type="mask") @@ -103,7 +103,7 @@ def test_thresholds(numba_conditional): pytest.warns(RuntimeWarning, match="invalid value"), ): # NumPy out = permutation_cluster_1samp_test( - X, seed=0, threshold=dict(start=0, step=0.1), out_type="mask" + X, rng=0, threshold=dict(start=0, step=0.1), out_type="mask" ) assert (out[2] < 0.05).any() assert not (out[2] < 0.05).all() @@ -112,7 +112,7 @@ def test_thresholds(numba_conditional): with np.errstate(invalid="ignore"): permutation_cluster_1samp_test( X, - seed=0, + rng=0, threshold=dict(start=0, step=0.1), buffer_size=None, out_type="mask", @@ -137,7 +137,7 @@ def test_cache_dir(tmp_path, numba_conditional): buffer_size=None, n_jobs=2, n_permutations=1, - seed=0, + rng=0, stat_fun=ttest_1samp_no_p, verbose=False, out_type="mask", @@ -152,7 +152,7 @@ def test_cache_dir(tmp_path, numba_conditional): buffer_size=10, n_jobs=2, n_permutations=1, - seed=random_state, + rng=random_state, stat_fun=stat_fun, verbose=False, out_type="mask", @@ -175,7 +175,7 @@ def test_permutation_large_n_samples(numba_conditional): tails = (0, 1) if n_samples <= 20 else (0,) for tail in tails: H0 = permutation_cluster_1samp_test( - X[:n_samples], threshold=1e-4, tail=tail, seed=0, out_type="mask" + X[:n_samples], threshold=1e-4, tail=tail, rng=0, out_type="mask" )[-1] assert H0.shape == (1024,) assert len(np.unique(H0)) >= 1024 - (H0 == 0).sum() @@ -221,7 +221,7 @@ def test_cluster_permutation_test(numba_conditional): [condition1, condition2], n_permutations=100, tail=1, - seed=1, + rng=1, buffer_size=None, out_type="mask", ) @@ -235,7 +235,7 @@ def test_cluster_permutation_test(numba_conditional): [condition1, condition2], n_permutations=100, tail=1, - seed=1, + rng=1, n_jobs=2, buffer_size=buffer_size, out_type="mask", @@ -268,7 +268,7 @@ def test_cluster_permutation_t_test(numba_conditional, stat_fun): condition1, n_permutations=100, tail=0, - seed=1, + rng=1, out_type="mask", buffer_size=None, ) @@ -281,7 +281,7 @@ def test_cluster_permutation_t_test(numba_conditional, stat_fun): n_permutations=100, tail=1, threshold=1.67, - seed=1, + rng=1, stat_fun=stat_fun, out_type="mask", buffer_size=None, @@ -292,7 +292,7 @@ def test_cluster_permutation_t_test(numba_conditional, stat_fun): n_permutations=100, tail=-1, threshold=-1.67, - seed=1, + rng=1, stat_fun=stat_fun, buffer_size=None, out_type="mask", @@ -314,7 +314,7 @@ def test_cluster_permutation_t_test(numba_conditional, stat_fun): tail=-1, out_type="mask", threshold=-1.67, - seed=1, + rng=1, n_jobs=2, stat_fun=stat_fun, buffer_size=buffer_size, @@ -347,7 +347,7 @@ def test_cluster_permutation_with_adjacency(numba_conditional, monkeypatch): n_pts = condition1_1d.shape[1] # we don't care about p-values in any of these, so do fewer permutations args = dict( - seed=None, + rng=None, max_step=1, exclude=None, out_type="mask", @@ -610,7 +610,7 @@ def test_permutation_adjacency_equiv(numba_conditional): n_jobs=2, max_step=max_step, stat_fun=stat_fun, - seed=0, + rng=0, out_type="mask", ) # make sure our output datatype is correct @@ -677,7 +677,7 @@ def test_spatio_temporal_cluster_chain_merge(): max_step=1, n_permutations=20, out_type="indices", - seed=0, + rng=0, verbose=False, ) assert len(clusters) == 1 @@ -760,7 +760,7 @@ def test_spatio_temporal_cluster_adjacency(numba_conditional): adjacency=adj, n_permutations=50, tail=1, - seed=1, + rng=1, threshold=threshold, buffer_size=None, ) @@ -770,7 +770,7 @@ def test_spatio_temporal_cluster_adjacency(numba_conditional): [data1_2d, data2_2d], n_permutations=50, tail=1, - seed=1, + rng=1, threshold=threshold, n_jobs=2, buffer_size=buffer_size, @@ -783,7 +783,7 @@ def test_spatio_temporal_cluster_adjacency(numba_conditional): [data1_2d, data2_2d], n_permutations=50, tail=1, - seed=1, + rng=1, threshold=threshold, n_jobs=2, buffer_size=None, @@ -866,19 +866,19 @@ def test_permutation_test_H0(numba_conditional): data = rng.random((7, 10, 1)) - 0.5 with pytest.warns(RuntimeWarning, match="No clusters found"): t, clust, p, h0 = spatio_temporal_cluster_1samp_test( - data, threshold=100, n_permutations=1024, seed=rng + data, threshold=100, n_permutations=1024, rng=rng ) assert_equal(len(h0), 0) for n_permutations in (1024, 65, 64, 63): t, clust, p, h0 = spatio_temporal_cluster_1samp_test( - data, threshold=0.1, n_permutations=n_permutations, seed=rng + data, threshold=0.1, n_permutations=n_permutations, rng=rng ) assert_equal(len(h0), min(n_permutations, 64)) assert isinstance(clust[0], tuple) # sets of indices for tail, thresh in zip((-1, 0, 1), (-0.1, 0.1, 0.1)): t, clust, p, h0 = spatio_temporal_cluster_1samp_test( - data, threshold=thresh, seed=rng, tail=tail, out_type="mask" + data, threshold=thresh, rng=rng, tail=tail, out_type="mask" ) assert isinstance(clust[0], np.ndarray) # bool mask # same as "128 if tail else 64" diff --git a/mne/stats/tests/test_permutations.py b/mne/stats/tests/test_permutations.py index 24d4db6d9d4..981b99cb169 100644 --- a/mne/stats/tests/test_permutations.py +++ b/mne/stats/tests/test_permutations.py @@ -24,24 +24,24 @@ def test_permutation_t_test(): X = rng.standard_normal((n_samples, n_tests)) X[:, :2] += 1 - t_obs, p_values, H0 = permutation_t_test(X, n_permutations=999, tail=0, seed=0) + t_obs, p_values, H0 = permutation_t_test(X, n_permutations=999, tail=0, rng=0) assert (p_values > 0).all() assert len(H0) == 999 is_significant = p_values < 0.05 assert_array_equal(is_significant, [True, True, False, False, False]) - t_obs, p_values, H0 = permutation_t_test(X, n_permutations=999, tail=1, seed=0) + t_obs, p_values, H0 = permutation_t_test(X, n_permutations=999, tail=1, rng=0) assert (p_values > 0).all() assert len(H0) == 999 is_significant = p_values < 0.05 assert_array_equal(is_significant, [True, True, False, False, False]) - t_obs, p_values, H0 = permutation_t_test(X, n_permutations=999, tail=-1, seed=0) + t_obs, p_values, H0 = permutation_t_test(X, n_permutations=999, tail=-1, rng=0) is_significant = p_values < 0.05 assert_array_equal(is_significant, [False, False, False, False, False]) X *= -1 - t_obs, p_values, H0 = permutation_t_test(X, n_permutations=999, tail=-1, seed=0) + t_obs, p_values, H0 = permutation_t_test(X, n_permutations=999, tail=-1, rng=0) assert (p_values > 0).all() assert len(H0) == 999 is_significant = p_values < 0.05 @@ -50,7 +50,7 @@ def test_permutation_t_test(): # check equivalence with spatio_temporal_cluster_test for adjacency in (sparse.eye_array(n_tests), False): t_obs_clust, _, p_values_clust, _ = permutation_cluster_1samp_test( - X, n_permutations=999, seed=0, adjacency=adjacency, out_type="mask" + X, n_permutations=999, rng=0, adjacency=adjacency, out_type="mask" ) # the cluster tests drop any clusters that don't get thresholded keep = p_values < 1 @@ -90,7 +90,7 @@ def test_ci(): _ci(arr, method="parametric"), _ci(arr, method="bootstrap"), rtol=0.005 ) assert_allclose( - bootstrap_confidence_interval(arr, stat_fun="median", random_state=0), - bootstrap_confidence_interval(arr, stat_fun="mean", random_state=0), + bootstrap_confidence_interval(arr, stat_fun="median", rng=0), + bootstrap_confidence_interval(arr, stat_fun="mean", rng=0), rtol=0.1, ) diff --git a/mne/stats/tests/test_regression.py b/mne/stats/tests/test_regression.py index 84c23a33487..dd4444fde1f 100644 --- a/mne/stats/tests/test_regression.py +++ b/mne/stats/tests/test_regression.py @@ -141,7 +141,9 @@ def test_continuous_regression_with_overlap(): def solver(X, y): # Newer scikit-learn returns 1D array for ridge_regression, so ensure # 2D output - return np.atleast_2d(ridge_regression(X, y, alpha=0.0, solver="cholesky")) + return np.atleast_2d( + ridge_regression(X, y, alpha=0.0, solver="cholesky", random_state=0) + ) assert_allclose( effect, @@ -150,7 +152,7 @@ def solver(X, y): # test bad solvers def solT(X, y): - return ridge_regression(X, y, alpha=0.0, solver="cholesky").T + return ridge_regression(X, y, alpha=0.0, solver="cholesky", random_state=0).T pytest.raises(ValueError, linear_regression_raw, raw, events, solver=solT) pytest.raises(ValueError, linear_regression_raw, raw, events, solver="err") diff --git a/mne/tests/test_docstring_parameters.py b/mne/tests/test_docstring_parameters.py index 9df4cf8cdaa..dfb0bb4af93 100644 --- a/mne/tests/test_docstring_parameters.py +++ b/mne/tests/test_docstring_parameters.py @@ -16,7 +16,7 @@ import pytest import mne -from mne.utils import _pl, _record_warnings +from mne.utils import _pl, _record_warnings, _soft_import from mne.utils._typing import Color, FileLike from mne.utils.docs import _doc_special_members @@ -328,6 +328,46 @@ def _is_np_random(node): ) +def _sklearn_callables(tree): + """Resolve sklearn imports for signature inspection.""" + callables = {} + for node in ast.walk(tree): + if ( + isinstance(node, ast.ImportFrom) + and node.module + and node.module.startswith("sklearn") + ): + module = _soft_import(node.module, "checking RNG parameters", strict=False) + if module: + for alias in node.names: + callable_ = getattr(module, alias.name, None) + if callable_ is not None: + callables[alias.asname or alias.name] = callable_ + elif isinstance(node, ast.Import): + for alias in node.names: + if alias.name.startswith("sklearn"): + module = _soft_import( + alias.name, "checking RNG parameters", strict=False + ) + if module: + callables[alias.asname or alias.name] = module + return callables + + +def _rng_parameters(callable_, node): + """Get RNG parameters from a callable, if inspectable.""" + try: + parameters = inspect.signature(callable_).parameters + except (TypeError, ValueError): + return set() + shuffle = next((kw.value for kw in node.keywords if kw.arg == "shuffle"), None) + if "shuffle" in parameters and ( + shuffle is None or isinstance(shuffle, ast.Constant) and not shuffle.value + ): + return set() + return {"random_state", "rng", "seed"} & parameters.keys() + + def test_no_global_rng(): """Test that we use local generators and the modern numpy RNG API.""" root = pyproject_path.parent # only available in a dev/editable checkout @@ -338,26 +378,35 @@ def test_no_global_rng(): continue for path in sorted(base.rglob("*.py")): rel = path.relative_to(root).as_posix() - for node in ast.walk(ast.parse(path.read_text("utf-8"))): - # 1. the global RNG: ``np.random.`` / ``numpy.random.`` + tree = ast.parse(path.read_text("utf-8")) + callables = _sklearn_callables(tree) + for node in ast.walk(tree): if ( isinstance(node, ast.Attribute) and node.attr not in global_rng_ok and _is_np_random(node.value) ): - bad.append( - f"{rel}:{node.lineno}: np.random.{node.attr} " - "(use a local np.random.default_rng)" - ) - # 2. legacy RandomState-only methods, e.g. ``rng.randn(...)`` + bad.append(f"{rel}:{node.lineno}: global np.random.{node.attr}") elif ( isinstance(node, ast.Call) and isinstance(node.func, ast.Attribute) and node.func.attr in legacy_rng_methods and not _is_np_random(node.func.value) ): - want = legacy_rng_methods[node.func.attr] - bad.append(f"{rel}:{node.lineno}: .{node.func.attr}() (use {want})") + bad.append(f"{rel}:{node.lineno}: legacy .{node.func.attr}()") + elif isinstance(node, ast.Call): + name = getattr(node.func, "id", None) or getattr( + node.func, "attr", None + ) + callable_ = callables.get(name) + if isinstance(node.func, ast.Attribute): + module = callables.get(getattr(node.func.value, "id", None)) + if module is not None: + callable_ = getattr(module, node.func.attr, None) + parameters = _rng_parameters(callable_, node) + supplied = parameters & {kw.arg for kw in node.keywords} + if parameters and len(supplied) != 1: + bad.append(f"{rel}:{node.lineno}: {name}() needs one RNG") if bad: raise AssertionError( f"{len(bad)} outdated numpy RNG use{_pl(bad)} found:\n" + "\n".join(bad) diff --git a/mne/tests/test_epochs.py b/mne/tests/test_epochs.py index 922ee4e4d8a..d56b863d611 100644 --- a/mne/tests/test_epochs.py +++ b/mne/tests/test_epochs.py @@ -2723,12 +2723,14 @@ def test_bootstrap(): reject=reject, flat=flat, ) - random_states = [0, np.random.default_rng(0)] - for random_state in random_states: - epochs2 = bootstrap(epochs, random_state=random_state) + rngs = [0, np.random.default_rng(0)] + for rng in rngs: + epochs2 = bootstrap(epochs, rng=rng) assert len(epochs2.events) == len(epochs.events) assert epochs._data.shape == epochs2._data.shape + bootstrap(epochs, random_state=0) + def test_epochs_copy(): """Test copy epochs.""" diff --git a/mne/tests/test_label.py b/mne/tests/test_label.py index 19290bfd42d..eca3cc5ee07 100644 --- a/mne/tests/test_label.py +++ b/mne/tests/test_label.py @@ -1029,7 +1029,7 @@ def test_random_parcellation(): # Parcellation labels = random_parcellation( - subject, n_parcel, hemi, subjects_dir, surface=surface, random_state=rng + subject, n_parcel, hemi, subjects_dir, surface=surface, rng=rng ) # test number of labels @@ -1213,6 +1213,14 @@ def test_select_sources(): assert label.hemi == "rh" +def test_select_sources_rng_conflict_at_center(): + """Test RNG spelling conflicts when random selection is inactive.""" + with pytest.raises(TypeError, match="only one"): + select_sources(None, None, location="center", random_state=0, rng=1) + with pytest.raises(TypeError, match="only one"): + select_sources(None, None, location="center", random_state=None, rng=None) + + @testing.requires_testing_data @pytest.mark.parametrize( "fname, area", diff --git a/mne/utils/__init__.pyi b/mne/utils/__init__.pyi index 9d70ef0ba57..d607d763a1c 100644 --- a/mne/utils/__init__.pyi +++ b/mne/utils/__init__.pyi @@ -52,6 +52,8 @@ __all__ = [ "_check_qt_version", "_check_range", "_check_rank", + "_check_rng", + "_legacy_rng", "_check_sphere", "_check_src_normal", "_check_stc_units", @@ -259,6 +261,7 @@ from .check import ( _check_qt_version, _check_range, _check_rank, + _check_rng, _check_sphere, _check_src_normal, _check_stc_units, @@ -271,6 +274,7 @@ from .check import ( _import_nibabel, _import_pymatreader_funcs, _is_numeric, + _legacy_rng, _on_missing, _path_like, _require_version, diff --git a/mne/utils/check.py b/mne/utils/check.py index 5c41c628025..1c95b198f02 100644 --- a/mne/utils/check.py +++ b/mne/utils/check.py @@ -10,6 +10,7 @@ import re from builtins import input # noqa: A004, UP029 from difflib import get_close_matches +from functools import wraps from importlib import import_module from inspect import signature from pathlib import Path @@ -230,6 +231,43 @@ def check_random_state(seed): ) +def _check_rng(rng): + """Return a NumPy Generator, or a legacy RandomState unchanged. + + Legacy RandomState instances are accepted for interoperability with + third-party code such as scikit-learn that does not accept Generator + instances. + """ + if isinstance(rng, np.random.mtrand.RandomState): + return rng + return np.random.default_rng(rng) + + +def _legacy_rng(legacy_name): + """Handle presence-sensitive legacy RNG parameters at the call boundary. + + The decorated function must accept a keyword-only ``rng`` parameter. When + it is called, ``kwargs["rng"]`` is replaced by the normalized value before + the function runs, so the body only ever sees an already-normalized RNG. + """ + + def decorator(function): + @wraps(function) + def _legacy_rng_wrapper(*args, **kwargs): + if legacy_name not in kwargs: + kwargs["rng"] = _check_rng(kwargs.get("rng")) + return function(*args, **kwargs) + if "rng" in kwargs: + raise TypeError(f"Specify only one of rng or {legacy_name}") + logger.info(f"Use rng= instead of {legacy_name}= in new code") + kwargs["rng"] = check_random_state(kwargs[legacy_name]) + return function(*args, **kwargs) + + return _legacy_rng_wrapper + + return decorator + + def _check_event_id(event_id, events): """Check event_id and convert to default format.""" # check out event_id dict diff --git a/mne/utils/docs.py b/mne/utils/docs.py index 868ee4fe0e4..10eb70e4046 100644 --- a/mne/utils/docs.py +++ b/mne/utils/docs.py @@ -3760,13 +3760,22 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75): docdict["random_state"] = """ random_state : None | int | instance of ~numpy.random.RandomState A seed for the NumPy random number generator (RNG). If ``None`` (default), - the seed will be obtained from the operating system - (see :class:`~numpy.random.RandomState` for details), meaning it will most - likely produce different output every time this function or method is run. - To achieve reproducible results, pass a value here to explicitly initialize - the RNG with a defined state. + NumPy's global :class:`~numpy.random.RandomState` singleton is used. + Pass an int to use a new ``RandomState`` seeded with that value, or a + ``RandomState`` to control the random-number stream. """ +docdict["random_state_rng"] = """ +random_state : None | int | instance of ~numpy.random.RandomState + Supported for compatibility. New code should use ``rng``. If ``None``, + NumPy's global :class:`~numpy.random.RandomState` is used. +""" + +docdict["random_state_rng_method_random"] = ( + docdict["random_state_rng"].rstrip("\n") + + "\n Used only if ``method='random'``.\n" +) + _rank_base = """ rank : None | 'info' | 'full' | dict This controls the rank computation that can be read from the @@ -4034,6 +4043,28 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75): Default to False. """ +docdict["rng"] = """ +rng : None | int | instance of ~numpy.random.Generator | ~numpy.random.RandomState + The random number generator (RNG). If ``None`` (default), a new + :class:`numpy.random.Generator` seeded from entropy is used. Pass an int or + a :class:`numpy.random.Generator` for reproducible results, or a legacy + :class:`~numpy.random.RandomState` to control the random-number stream or + for interoperability with third-party code such as scikit-learn that does + not accept generators. An integer seed uses + :func:`numpy.random.default_rng` and therefore produces a different stream + than the same integer passed to a legacy ``random_state`` or ``seed`` + parameter. + + .. versionadded:: 1.13 +""" + +# The ``rng`` entry ends with a directive, so anything appended at the call site +# would be swallowed by it; make the ``method='random'`` variant here instead. +docdict["rng_method_random"] = docdict["rng"].replace( + "\n\n .. versionadded", + "\n Used only if ``method='random'``.\n\n .. versionadded", +) + docdict["roll"] = """ roll : float | None The roll of the camera rendering the view in degrees. @@ -4127,13 +4158,13 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75): docdict["seed"] = """ seed : None | int | instance of ~numpy.random.RandomState A seed for the NumPy random number generator (RNG). If ``None`` (default), - the seed will be obtained from the operating system - (see :class:`~numpy.random.RandomState` for details), meaning it will most - likely produce different output every time this function or method is run. - To achieve reproducible results, pass a value here to explicitly initialize - the RNG with a defined state. + NumPy's global :class:`~numpy.random.RandomState` singleton is used. + Pass an int to use a new ``RandomState`` seeded with that value, or a + ``RandomState`` to control the random-number stream. """ +docdict["seed_rng"] = docdict["random_state_rng"].replace("random_state", "seed") + docdict["seeg"] = """ seeg : bool If True (default), show sEEG electrodes. diff --git a/mne/utils/numerics.py b/mne/utils/numerics.py index 41fb91e5cf5..3c25ae2d745 100644 --- a/mne/utils/numerics.py +++ b/mne/utils/numerics.py @@ -27,8 +27,8 @@ from ._logging import logger, verbose, warn from .check import ( _ensure_int, + _legacy_rng, _validate_type, - check_random_state, ) from .docs import fill_doc from .misc import _empty_hash, _pl @@ -265,8 +265,9 @@ def compute_corr(x, y): return (np.dot(X.T, Y) / float(len(X) - 1)) / (x_sd * y_sd) +@_legacy_rng("random_state") @fill_doc -def random_permutation(n_samples, random_state=None): +def random_permutation(n_samples, *, rng=None, random_state=None): """Emulate the randperm matlab function. It returns a vector containing a random permutation of the @@ -288,14 +289,19 @@ def random_permutation(n_samples, random_state=None): n_samples : int End point of the sequence to be permuted (excluded, i.e., the end point is equal to n_samples-1) - %(random_state)s + %(rng)s + %(random_state_rng)s Returns ------- randperm : ndarray, int Randomly permuted sequence between 0 and n-1. """ - rng = check_random_state(random_state) + return _random_permutation(n_samples, rng) + + +def _random_permutation(n_samples, rng): + """Generate a MATLAB-compatible permutation with a normalized RNG.""" # This can't just be rng.permutation(n_samples) because it's not identical # to what MATLAB produces idx = rng.uniform(size=n_samples) diff --git a/mne/utils/tests/test_check.py b/mne/utils/tests/test_check.py index fbeea9b6136..0c455d64d26 100644 --- a/mne/utils/tests/test_check.py +++ b/mne/utils/tests/test_check.py @@ -10,7 +10,7 @@ import numpy as np import pytest -from numpy.testing import assert_allclose, assert_equal +from numpy.testing import assert_allclose, assert_array_equal, assert_equal import mne from mne import create_info, pick_channels_cov, read_vectorview_selection @@ -38,6 +38,7 @@ check_random_state, check_version, ) +from mne.utils.check import _check_rng, _legacy_rng data_path = testing.data_path(download=False) base_dir = data_path / "MEG" / "sample" @@ -48,6 +49,46 @@ reject = dict(grad=4000e-13, mag=4e-12) +def test_check_rng(): + """Test conversion to NumPy's modern random number generator.""" + assert isinstance(_check_rng(None), np.random.Generator) + + rng = _check_rng(0) + assert isinstance(rng, np.random.Generator) + assert_array_equal(rng.integers(10, size=3), _check_rng(0).integers(10, size=3)) + + assert _check_rng(rng) is rng + # legacy RandomState instances are passed through for scikit-learn interop + random_state = np.random.RandomState(0) + assert _check_rng(random_state) is random_state + with pytest.raises(TypeError, match="SeedSequence"): + _check_rng("foo") + + +@pytest.mark.parametrize("legacy_name", ("random_state", "seed")) +def test_legacy_rng_decorator(legacy_name): + """Test that the transition decorator normalizes and logs.""" + + @_legacy_rng(legacy_name) + def _func(*, rng=None, random_state=None, seed=None): + return rng + + # no argument: a fresh generator is created + assert isinstance(_func(), np.random.Generator) + assert isinstance(_func(rng=0), np.random.Generator) + # legacy RandomState passthrough + random_state = np.random.RandomState(0) + assert _func(rng=random_state) is random_state + # legacy int/None keep their RandomState semantics and log migration guidance + with catch_logging(verbose="info") as log: + assert isinstance(_func(**{legacy_name: 0}), np.random.RandomState) + assert f"Use rng= instead of {legacy_name}=" in log.getvalue() + assert isinstance(_func(**{legacy_name: None}), np.random.mtrand.RandomState) + # supplying both is an error + with pytest.raises(TypeError, match="Specify only one"): + _func(**{legacy_name: 0}, rng=0) + + @testing.requires_testing_data def test_check(tmp_path): """Test checking functions.""" diff --git a/mne/utils/tests/test_numerics.py b/mne/utils/tests/test_numerics.py index ea820cd9c75..a679e3d5dfd 100644 --- a/mne/utils/tests/test_numerics.py +++ b/mne/utils/tests/test_numerics.py @@ -213,15 +213,21 @@ def test_freq_mask(): def test_random_permutation(): - """Test random permutation function.""" + """Test random permutation function and its RNG transition.""" n_samples = 10 - random_state = 42 - python_randperm = random_permutation(n_samples, random_state) - # matlab output when we execute rng(42), randperm(10) - matlab_randperm = np.array([7, 6, 5, 1, 4, 9, 10, 3, 8, 2]) - - assert_array_equal(python_randperm, matlab_randperm - 1) + assert_array_equal( + random_permutation(n_samples, random_state=42), + np.array([7, 6, 5, 1, 4, 9, 10, 3, 8, 2]) - 1, + ) + # an integer ``rng`` seed is reproducible while a Generator instance advances + assert_array_equal( + random_permutation(n_samples, rng=42), random_permutation(n_samples, rng=42) + ) + rng = np.random.default_rng(42) + assert not np.array_equal( + random_permutation(n_samples, rng=rng), random_permutation(n_samples, rng=rng) + ) def test_cov_scaling(): @@ -456,7 +462,7 @@ def test_pca(n_components, whiten): X = np.random.default_rng(0).standard_normal((n_samples, n_dim)) X[:, -1] = np.mean(X[:, :-1], axis=-1) # true X dim is ndim - 1 X_orig = X.copy() - pca_skl = PCA(n_components, whiten=whiten, svd_solver="full") + pca_skl = PCA(n_components, whiten=whiten, svd_solver="full", random_state=0) pca_mne = _PCA(n_components, whiten=whiten) X_skl = pca_skl.fit_transform(X) assert_array_equal(X, X_orig) diff --git a/mne/viz/circle.py b/mne/viz/circle.py index 2b2b39705f9..7fc2e8d850c 100644 --- a/mne/viz/circle.py +++ b/mne/viz/circle.py @@ -268,7 +268,7 @@ def _plot_connectivity_circle( nodes_n_con[j] += 1 # initialize random number generator so plot is reproducible - rng = np.random.mtrand.RandomState(0) + rng = np.random.default_rng(0) n_con = len(indices[0]) noise_max = 0.25 * node_width diff --git a/mne/viz/tests/test_ica.py b/mne/viz/tests/test_ica.py index 787f2175785..332b9d10805 100644 --- a/mne/viz/tests/test_ica.py +++ b/mne/viz/tests/test_ica.py @@ -89,7 +89,7 @@ def test_plot_ica_components(): res = 8 fast_test = {"res": res, "contours": 0, "sensors": False} raw = _get_raw() - ica = ICA(noise_cov=read_cov(cov_fname), n_components=8, random_state=0) + ica = ICA(noise_cov=read_cov(cov_fname), n_components=8, rng=0) ica_picks = _get_picks(raw) with pytest.warns(RuntimeWarning, match="(projection)|(unstable mixing matrix)"): ica.fit(raw, picks=ica_picks) @@ -189,7 +189,7 @@ def test_plot_ica_properties_basic(): raw, events[:3], event_id, tmin, tmax, baseline=(None, 0), preload=True ) - ica = ICA(noise_cov=read_cov(cov_fname), n_components=2, max_iter=1, random_state=0) + ica = ICA(noise_cov=read_cov(cov_fname), n_components=2, max_iter=1, rng=0) with _record_warnings(), pytest.warns(RuntimeWarning, match="projection"): ica.fit(raw) @@ -264,14 +264,14 @@ def test_plot_ica_properties_basic(): raw = _get_raw(preload=True).pick(pick_names) raw.crop(0, 5) raw.info.normalize_proj() - ica = ICA(random_state=0, max_iter=1) + ica = ICA(rng=0, max_iter=1) with pytest.warns(UserWarning, match="did not converge"): ica.fit(raw) ica.plot_properties(raw) plt.close("all") # Test handling of zeros - ica = ICA(random_state=0, max_iter=1) + ica = ICA(rng=0, max_iter=1) epochs.pick(pick_names) with _record_warnings(), pytest.warns(UserWarning, match="did not converge"): ica.fit(epochs) @@ -326,7 +326,7 @@ def test_plot_ica_properties_reject(kind): raw.set_montage("spherical_1005") ica = ICA( n_components=2, - random_state=0, + rng=0, max_iter=1, ) with ( @@ -363,7 +363,7 @@ def test_plot_ica_sources(raw_orig, browser_backend, monkeypatch): ica_picks = pick_types( raw.info, meg=True, eeg=False, stim=False, ecg=False, eog=False, exclude="bads" ) - ica = ICA(n_components=2, random_state=0) + ica = ICA(n_components=2, rng=0) ica.fit(raw, picks=ica_picks) ica.exclude = [1] if sys.platform == "darwin": # unknown transformation bug @@ -544,7 +544,7 @@ def test_plot_ica_overlay(): with raw.info._unlock(): raw.info["highpass"] = 1.0 # fake high-pass filtering picks = _get_picks(raw) - ica = ICA(noise_cov=read_cov(cov_fname), n_components=2, random_state=0) + ica = ICA(noise_cov=read_cov(cov_fname), n_components=2, rng=0) # overlay plotting requires a fitted ICA with pytest.raises(RuntimeError, match="need to fit"): ica.plot_overlay(inst=raw) @@ -571,7 +571,7 @@ def test_plot_ica_overlay(): picks = pick_types(raw.info, meg=True, ref_meg=False) ica = ICA( n_components=2, - random_state=0, + rng=0, ) ica.fit(raw, picks=picks) with pytest.warns(RuntimeWarning, match="longer than"): @@ -590,7 +590,7 @@ def test_plot_ica_scores(): """Test plotting of ICA scores.""" raw = _get_raw() picks = _get_picks(raw) - ica = ICA(noise_cov=read_cov(cov_fname), n_components=2, random_state=0) + ica = ICA(noise_cov=read_cov(cov_fname), n_components=2, rng=0) with pytest.warns(RuntimeWarning, match="projection"): ica.fit(raw, picks=picks) ica.plot_scores([0.3, 0.2], axhline=[0.1, -0.1], figsize=(6.4, 2.7)) @@ -628,7 +628,7 @@ def test_plot_instance_components(browser_backend): """Test plotting of components as instances of raw and epochs.""" raw = _get_raw() picks = _get_picks(raw) - ica = ICA(noise_cov=read_cov(cov_fname), n_components=2, random_state=0) + ica = ICA(noise_cov=read_cov(cov_fname), n_components=2, rng=0) with pytest.warns(RuntimeWarning, match="projection"): ica.fit(raw, picks=picks) ica.exclude = [0] @@ -681,7 +681,7 @@ def test_plot_instance_components(browser_backend): def test_plot_components_opm(): """Test for gh-12934.""" evoked = read_evokeds(opm_fname, kind="average")[0] - ica = ICA(max_iter=1, random_state=0, n_components=10) + ica = ICA(max_iter=1, rng=0, n_components=10) ica.fit(RawArray(evoked.data, evoked.info), picks="mag", verbose="error") fig = ica.plot_components() # Biaxial OPM overlaps render grouped radial+tangential maps. @@ -694,7 +694,7 @@ def test_plot_components_opm(): ) def test_plot_components_opm_triaxial(triaxial_raw): """Test OPM component topomaps with colocated triaxial channels.""" - ica = ICA(max_iter=1, random_state=0, n_components=3) + ica = ICA(max_iter=1, rng=0, n_components=3) ica.fit(triaxial_raw, picks="mag", verbose="error") fig = ica.plot_components() assert len(fig.axes) == 6 diff --git a/mne/viz/tests/test_topomap.py b/mne/viz/tests/test_topomap.py index 072c1668eb7..9bd88662b2a 100644 --- a/mne/viz/tests/test_topomap.py +++ b/mne/viz/tests/test_topomap.py @@ -1006,7 +1006,7 @@ def test_plot_topomap_nirs_ica(fnirs_epochs): fnirs_epochs.info["highpass"] = 1.0 fnirs_epochs.baseline = None - ica = ICA(random_state=0).fit(fnirs_epochs) + ica = ICA(rng=0).fit(fnirs_epochs) fig = ica.plot_components() assert len(fig[0].axes) == 20 diff --git a/mne/viz/utils.py b/mne/viz/utils.py index 7e0b48f7e7e..80a5b146ec8 100644 --- a/mne/viz/utils.py +++ b/mne/viz/utils.py @@ -1403,7 +1403,7 @@ def _compute_scalings(scalings, inst, remove_dc=False, duration=10): # Load a random subset of epochs up to 100mb in size n_epochs = 1e8 // (len(inst.ch_names) * len(inst.times) * 8) n_epochs = int(np.clip(n_epochs, 1, len(inst))) - ixs_epochs = np.random.default_rng().choice( + ixs_epochs = np.random.default_rng(0).choice( len(inst), n_epochs, replace=False ) inst = inst.copy()[ixs_epochs].load_data() diff --git a/tools/vulture_allowlist.py b/tools/vulture_allowlist.py index a42a635ef20..811116d9adb 100644 --- a/tools/vulture_allowlist.py +++ b/tools/vulture_allowlist.py @@ -65,6 +65,7 @@ # Backward compat or rarely used RawFIF +select_source_in_label estimate_head_mri_t plot_epochs_psd_topomap plot_epochs_psd diff --git a/tutorials/clinical/60_sleep.py b/tutorials/clinical/60_sleep.py index e50b0740a7c..50755b5d411 100644 --- a/tutorials/clinical/60_sleep.py +++ b/tutorials/clinical/60_sleep.py @@ -301,7 +301,7 @@ def eeg_power_band(epochs): pipe = make_pipeline( FunctionTransformer(eeg_power_band, validate=False), - RandomForestClassifier(n_estimators=100, random_state=42), + RandomForestClassifier(n_estimators=100, random_state=241), ) # Train diff --git a/tutorials/intro/10_overview.py b/tutorials/intro/10_overview.py index d9fd1cef8d2..8d5f2acbcc5 100644 --- a/tutorials/intro/10_overview.py +++ b/tutorials/intro/10_overview.py @@ -98,7 +98,7 @@ # :ref:`tut-artifact-ica` for a detailed walk-through of that process). # set up and fit the ICA -ica = mne.preprocessing.ICA(n_components=20, random_state=97, max_iter=800) +ica = mne.preprocessing.ICA(n_components=20, rng=97, max_iter=800) ica.fit(raw) ica.exclude = [1, 2] # details on how we picked these are omitted here ica.plot_properties(raw, picks=ica.exclude) diff --git a/tutorials/intro/70_report.py b/tutorials/intro/70_report.py index f76a6a3cfbc..eb2bbc33f0a 100644 --- a/tutorials/intro/70_report.py +++ b/tutorials/intro/70_report.py @@ -246,7 +246,7 @@ ica = mne.preprocessing.ICA( n_components=5, # fit 5 ICA components fit_params=dict(tol=0.01), # assume very early on that ICA has converged - random_state=97, + rng=97, ) ica.fit(inst=raw) diff --git a/tutorials/machine-learning/50_decoding.py b/tutorials/machine-learning/50_decoding.py index 388b0f5a813..87bfffa15b7 100644 --- a/tutorials/machine-learning/50_decoding.py +++ b/tutorials/machine-learning/50_decoding.py @@ -141,7 +141,8 @@ clf = make_pipeline( Scaler(epochs.info), Vectorizer(), - LogisticRegression(solver="liblinear"), # liblinear is faster than lbfgs + # liblinear is faster than lbfgs + LogisticRegression(solver="liblinear", random_state=31), ) scores = cross_val_multiscore(clf, X, y, cv=5, n_jobs=None) @@ -225,7 +226,9 @@ # We can use CSP with these data with: csp = CSP(n_components=3, norm_trace=False) -clf_csp = make_pipeline(csp, LinearModel(LogisticRegression(solver="liblinear"))) +clf_csp = make_pipeline( + csp, LinearModel(LogisticRegression(solver="liblinear", random_state=37)) +) scores = cross_val_multiscore(clf_csp, X, y, cv=5, n_jobs=None) print(f"CSP: {100 * scores.mean():0.1f}%") @@ -324,7 +327,9 @@ # We will train the classifier on all left visual vs auditory trials on MEG -clf = make_pipeline(StandardScaler(), LogisticRegression(solver="liblinear")) +clf = make_pipeline( + StandardScaler(), LogisticRegression(solver="liblinear", random_state=41) +) time_decod = SlidingEstimator(clf, n_jobs=None, scoring="roc_auc", verbose=True) # here we use cv=3 just for speed @@ -347,7 +352,8 @@ # You can retrieve the spatial filters and spatial patterns if you explicitly # use a LinearModel clf = make_pipeline( - StandardScaler(), LinearModel(LogisticRegression(solver="liblinear")) + StandardScaler(), + LinearModel(LogisticRegression(solver="liblinear", random_state=43)), ) time_decod = SlidingEstimator(clf, n_jobs=None, scoring="roc_auc", verbose=True) time_decod.fit(X, y) diff --git a/tutorials/preprocessing/14_quality_control_report.py b/tutorials/preprocessing/14_quality_control_report.py index 1880ae014b2..708ec487b50 100644 --- a/tutorials/preprocessing/14_quality_control_report.py +++ b/tutorials/preprocessing/14_quality_control_report.py @@ -204,7 +204,7 @@ ica = ICA( n_components=15, - random_state=97, + rng=97, max_iter=50, # just for speed! ) diff --git a/tutorials/preprocessing/40_artifact_correction_ica.py b/tutorials/preprocessing/40_artifact_correction_ica.py index 670ba577316..4afc3d415ec 100644 --- a/tutorials/preprocessing/40_artifact_correction_ica.py +++ b/tutorials/preprocessing/40_artifact_correction_ica.py @@ -255,7 +255,7 @@ # **after** cleaning (and not before), should you require # baseline correction. -ica = ICA(n_components=15, max_iter="auto", random_state=97) +ica = ICA(n_components=15, max_iter="auto", rng=97) ica.fit(filt_raw, reject=dict(eeg=200e-6)) # avoid a couple of big artifacts ica @@ -470,7 +470,7 @@ # resolves out a little better: # refit the ICA with 30 components this time -new_ica = ICA(n_components=30, max_iter="auto", random_state=97) +new_ica = ICA(n_components=30, max_iter="auto", rng=97) new_ica.fit(filt_raw) # find which ICs match the ECG pattern @@ -542,7 +542,7 @@ # high-pass filter raw_filt = raw.copy().load_data().filter(l_freq=1.0, h_freq=None) # fit ICA, using low max_iter for speed - ica = ICA(n_components=30, max_iter=100, random_state=97) + ica = ICA(n_components=30, max_iter=100, rng=97) ica.fit(raw_filt, verbose="error") raws.append(raw) icas.append(ica) @@ -673,7 +673,7 @@ # Fit ICA model using the FastICA algorithm, detect and plot components # explaining ECG artifacts. -ica = ICA(n_components=15, method="fastica", max_iter="auto", random_state=97) +ica = ICA(n_components=15, method="fastica", max_iter="auto", rng=97) ica.fit(epochs) ecg_epochs = create_ecg_epochs(filt_raw, tmin=-0.5, tmax=0.5) diff --git a/tutorials/simulation/70_point_spread.py b/tutorials/simulation/70_point_spread.py index 485714e2c17..f803a398da6 100644 --- a/tutorials/simulation/70_point_spread.py +++ b/tutorials/simulation/70_point_spread.py @@ -26,7 +26,7 @@ # %% # First, we set some parameters. -seed = 42 +seed = 199 # parameters for inverse method method = "sLORETA" @@ -161,7 +161,7 @@ # (evoked) data from the known source-space signals. The amount of noise is # controlled by ``nave`` (higher values imply less noise). # -evoked_gen = simulate_evoked(fwd, stc_gen, evoked.info, cov, nave, random_state=seed) +evoked_gen = simulate_evoked(fwd, stc_gen, evoked.info, cov, nave, rng=seed) # Map the simulated sensor-space data to source-space using the inverse # operator. diff --git a/tutorials/simulation/80_dics.py b/tutorials/simulation/80_dics.py index a487db59752..5b00ead4fee 100644 --- a/tutorials/simulation/80_dics.py +++ b/tutorials/simulation/80_dics.py @@ -189,7 +189,7 @@ def coh_signal_gen(): ] # stacked in time duration = (len(stc_signal.times) * 2) / sfreq raw = simulate_raw(info, stcs, forward=fwd) -add_noise(raw, cov, iir_filter=[4, -4, 0.8], random_state=rand) +add_noise(raw, cov, iir_filter=[4, -4, 0.8], rng=rand) # %% diff --git a/tutorials/stats-sensor-space/10_background_stats.py b/tutorials/stats-sensor-space/10_background_stats.py index 7688794a195..1a25e741b7c 100644 --- a/tutorials/stats-sensor-space/10_background_stats.py +++ b/tutorials/stats-sensor-space/10_background_stats.py @@ -251,7 +251,7 @@ def plot_t_p(t, p, title, mcc, axes=None): ps.append(np.zeros(width * width)) mccs.append(False) for ii in range(n_src): - t, p = permutation_t_test(X[:, [ii]], verbose=False, seed=0)[:2] + t, p = permutation_t_test(X[:, [ii]], verbose=False, rng=47)[:2] ts[-1][ii], ps[-1][ii] = t[0], p[0] plot_t_p(ts[-1], ps[-1], titles[-1], mccs[-1]) @@ -370,7 +370,7 @@ def plot_t_p(t, p, title, mcc, axes=None): # of processed neuroimaging data). titles.append(r"$\mathbf{Perm_{max}}$") -out = permutation_t_test(X, verbose=False, seed=0)[:2] +out = permutation_t_test(X, verbose=False, rng=53)[:2] ts.append(out[0]) ps.append(out[1]) mccs.append(True) @@ -509,7 +509,7 @@ def plot_t_p(t, p, title, mcc, axes=None): # run the cluster test t_clust, clusters, p_values, H0 = permutation_cluster_1samp_test( X, - seed=0, + rng=59, n_jobs=None, threshold=t_thresh, adjacency=None, @@ -535,7 +535,7 @@ def plot_t_p(t, p, title, mcc, axes=None): stat_fun_hat = partial(ttest_1samp_no_p, sigma=sigma) t_hat, clusters, p_values, H0 = permutation_cluster_1samp_test( X, - seed=0, + rng=61, n_jobs=None, threshold=t_thresh, adjacency=None, @@ -579,7 +579,7 @@ def plot_t_p(t, p, title, mcc, axes=None): threshold_tfce = dict(start=0, step=0.2) t_tfce, _, p_tfce, H0 = permutation_cluster_1samp_test( X, - seed=0, + rng=67, n_jobs=None, threshold=threshold_tfce, adjacency=None, @@ -596,7 +596,7 @@ def plot_t_p(t, p, title, mcc, axes=None): titles.append(r"$\mathbf{C_{hat,TFCE}}$") t_tfce_hat, _, p_tfce_hat, H0 = permutation_cluster_1samp_test( X, - seed=0, + rng=71, n_jobs=None, threshold=threshold_tfce, adjacency=None, diff --git a/tutorials/stats-sensor-space/20_erp_stats.py b/tutorials/stats-sensor-space/20_erp_stats.py index 3b0b8af742b..17c3655845d 100644 --- a/tutorials/stats-sensor-space/20_erp_stats.py +++ b/tutorials/stats-sensor-space/20_erp_stats.py @@ -98,7 +98,7 @@ # Calculate statistical thresholds t_obs, clusters, cluster_pv, h0 = spatio_temporal_cluster_test( - X, tfce, adjacency=adjacency, n_permutations=100, seed=0 + X, tfce, adjacency=adjacency, n_permutations=100, rng=113 ) # a more standard number would be 1000+ significant_points = cluster_pv.reshape(t_obs.shape).T < 0.05 print(str(significant_points.sum()) + " points selected by TFCE ...") diff --git a/tutorials/stats-sensor-space/40_cluster_1samp_time_freq.py b/tutorials/stats-sensor-space/40_cluster_1samp_time_freq.py index c763a9af44f..3ad9b11e2e0 100644 --- a/tutorials/stats-sensor-space/40_cluster_1samp_time_freq.py +++ b/tutorials/stats-sensor-space/40_cluster_1samp_time_freq.py @@ -208,7 +208,7 @@ tail=tail, adjacency=adjacency, out_type="mask", - seed=0, + rng=73, verbose=True, ) diff --git a/tutorials/stats-sensor-space/50_cluster_between_time_freq.py b/tutorials/stats-sensor-space/50_cluster_between_time_freq.py index 0b4078ec883..7fff7d7197c 100644 --- a/tutorials/stats-sensor-space/50_cluster_between_time_freq.py +++ b/tutorials/stats-sensor-space/50_cluster_between_time_freq.py @@ -131,7 +131,7 @@ n_permutations=100, threshold=threshold, tail=0, - seed=np.random.default_rng(seed=8675309), + rng=np.random.default_rng(seed=8675309), ) # %% diff --git a/tutorials/stats-sensor-space/70_cluster_rmANOVA_time_freq.py b/tutorials/stats-sensor-space/70_cluster_rmANOVA_time_freq.py index c8ec0c5f0d5..febfdcb1a5e 100644 --- a/tutorials/stats-sensor-space/70_cluster_rmANOVA_time_freq.py +++ b/tutorials/stats-sensor-space/70_cluster_rmANOVA_time_freq.py @@ -239,7 +239,7 @@ def stat_fun(*args): n_permutations=n_permutations, buffer_size=None, out_type="mask", - seed=0, + rng=101, ) # %% diff --git a/tutorials/stats-sensor-space/75_cluster_ftest_spatiotemporal.py b/tutorials/stats-sensor-space/75_cluster_ftest_spatiotemporal.py index 44d4748fd88..9df6f787873 100644 --- a/tutorials/stats-sensor-space/75_cluster_ftest_spatiotemporal.py +++ b/tutorials/stats-sensor-space/75_cluster_ftest_spatiotemporal.py @@ -144,7 +144,7 @@ n_jobs=None, buffer_size=None, adjacency=adjacency, - seed=0, + rng=149, ) F_obs, clusters, p_values, _ = cluster_stats @@ -318,7 +318,7 @@ n_jobs=None, buffer_size=None, adjacency=tfr_adjacency, - seed=0, + rng=151, ) # %% diff --git a/tutorials/stats-source-space/20_cluster_1samp_spatiotemporal.py b/tutorials/stats-source-space/20_cluster_1samp_spatiotemporal.py index fc7ca962d44..f767cdf42c6 100644 --- a/tutorials/stats-source-space/20_cluster_1samp_spatiotemporal.py +++ b/tutorials/stats-source-space/20_cluster_1samp_spatiotemporal.py @@ -209,7 +209,7 @@ n_jobs=None, threshold=t_threshold, buffer_size=None, - seed=0, + rng=179, verbose=True, ) diff --git a/tutorials/stats-source-space/30_cluster_ftest_spatiotemporal.py b/tutorials/stats-source-space/30_cluster_ftest_spatiotemporal.py index 14a4488f7e6..890e2cbddf9 100644 --- a/tutorials/stats-source-space/30_cluster_ftest_spatiotemporal.py +++ b/tutorials/stats-source-space/30_cluster_ftest_spatiotemporal.py @@ -101,7 +101,7 @@ n_permutations=n_permutations, threshold=f_threshold, buffer_size=None, - seed=0, + rng=167, ) # Now select the clusters that are sig. at p < 0.05 (note that this value # is multiple-comparisons corrected). diff --git a/tutorials/stats-source-space/60_cluster_rmANOVA_spatiotemporal.py b/tutorials/stats-source-space/60_cluster_rmANOVA_spatiotemporal.py index 1f51ab95b1e..61ddb3cfbf5 100644 --- a/tutorials/stats-source-space/60_cluster_rmANOVA_spatiotemporal.py +++ b/tutorials/stats-source-space/60_cluster_rmANOVA_spatiotemporal.py @@ -245,7 +245,7 @@ def stat_fun(*args): stat_fun=stat_fun, n_permutations=n_permutations, buffer_size=None, - seed=0, + rng=173, ) # Now select the clusters that are sig. at p < 0.05 (note that this value # is multiple-comparisons corrected).