From c72fa72877a2b92fd6cbfb5ce9fdcce182a49bb3 Mon Sep 17 00:00:00 2001 From: karp2601 Date: Tue, 24 Oct 2023 11:21:24 -0400 Subject: [PATCH 01/10] Adding support for btens to multi_shell_fiber_response function --- dipy/reconst/mcsd.py | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/dipy/reconst/mcsd.py b/dipy/reconst/mcsd.py index 2677d1fad6..fba1870bd6 100644 --- a/dipy/reconst/mcsd.py +++ b/dipy/reconst/mcsd.py @@ -437,7 +437,7 @@ def __call__(self, signal): def multi_shell_fiber_response(sh_order, bvals, wm_rf, gm_rf, csf_rf, - sphere=None, tol=20): + sphere=None, tol=20, btens=None): """Fiber response function estimation for multi-shell data. Parameters @@ -457,6 +457,17 @@ def multi_shell_fiber_response(sh_order, bvals, wm_rf, gm_rf, csf_rf, Sphere where the signal will be evaluated. tol : int, optional Tolerance gap for b-values clustering. + btens : can be any of two options, optional + + 1. an array of strings of shape (N,) specifying + encoding tensor shape associated with all unique b-values + separately. N corresponds to the number of unique b-values, + including the b0. Options for elements in array: 'LTE', + 'PTE', 'STE', 'CTE' corresponding to linear, planar, spherical, and + "cigar-shaped" tensor encoding. + 2. an array of shape (N,3,3) specifying the b-tensor of each unique + b-values exactly. N corresponds to the number of unique b-values, + including the b0. Returns ------- @@ -467,6 +478,11 @@ def multi_shell_fiber_response(sh_order, bvals, wm_rf, gm_rf, csf_rf, rcond_value = None if NUMPY_1_14_PLUS else -1 bvals = np.array(bvals, copy=True) + if btens is None: + btens = np.repeat(["LTE"], len(bvals)) + elif len(btens) != len(bvals): + msg = """bvals and btens parameters must have the same dimension.""" + raise ValueError(msg) evecs = np.zeros((3, 3)) z = np.array([0, 0, 1.]) evecs[:, 0] = z @@ -487,7 +503,7 @@ def multi_shell_fiber_response(sh_order, bvals, wm_rf, gm_rf, csf_rf, response = np.empty([len(bvals), len(n) + 2]) if bvals[0] < tol: - gtab = GradientTable(big_sphere.vertices * 0) + gtab = GradientTable(big_sphere.vertices * 0, btens=btens[0]) wm_response = single_tensor(gtab, wm_rf[0, 3], wm_rf[0, :3], evecs, snr=None) response[0, 2:] = np.linalg.lstsq(B, wm_response, rcond=rcond_value)[0] @@ -496,7 +512,8 @@ def multi_shell_fiber_response(sh_order, bvals, wm_rf, gm_rf, csf_rf, response[0, 0] = csf_rf[0, 3] / A for i, bvalue in enumerate(bvals[1:]): - gtab = GradientTable(big_sphere.vertices * bvalue) + gtab = GradientTable(big_sphere.vertices * bvalue, + btens=btens[i + 1]) wm_response = single_tensor(gtab, wm_rf[i, 3], wm_rf[i, :3], evecs, snr=None) response[i+1, 2:] = np.linalg.lstsq(B, wm_response, @@ -510,7 +527,8 @@ def multi_shell_fiber_response(sh_order, bvals, wm_rf, gm_rf, csf_rf, else: warnings.warn("""No b0 given. Proceeding either way.""", UserWarning) for i, bvalue in enumerate(bvals): - gtab = GradientTable(big_sphere.vertices * bvalue) + gtab = GradientTable(big_sphere.vertices * bvalue, + btens=btens[i]) wm_response = single_tensor(gtab, wm_rf[i, 3], wm_rf[i, :3], evecs, snr=None) response[i, 2:] = np.linalg.lstsq(B, wm_response, From 277fb8309fc59a4c7d73420832eb8f74a73a0605 Mon Sep 17 00:00:00 2001 From: karp2601 Date: Tue, 24 Oct 2023 15:11:44 -0400 Subject: [PATCH 02/10] Adding some basic test --- dipy/reconst/tests/test_mcsd.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/dipy/reconst/tests/test_mcsd.py b/dipy/reconst/tests/test_mcsd.py index d383c14e0a..7b2735d6ff 100644 --- a/dipy/reconst/tests/test_mcsd.py +++ b/dipy/reconst/tests/test_mcsd.py @@ -222,6 +222,19 @@ def test_multi_shell_fiber_response(): npt.assert_equal(response.response.shape, (4, 7)) + btens = ["LTE", "PTE", "STE", "CTE"] + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", message=shm.descoteaux07_legacy_msg, + category=PendingDeprecationWarning) + response = multi_shell_fiber_response(sh_order, [0, 1000, 2000, 3500], + wm_response, + gm_response, + csf_response, + btens=btens) + + npt.assert_equal(response.response.shape, (4, 7)) + with warnings.catch_warnings(record=True) as w: warnings.simplefilter("always", category=PendingDeprecationWarning) response = multi_shell_fiber_response(sh_order, [1000, 2000, 3500], From 8204e2ffb507ece5d6b764615bb19c9a6ce5f15c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jon=20Haitz=20Legarreta=20Gorro=C3=B1o?= Date: Thu, 26 Oct 2023 18:21:21 -0400 Subject: [PATCH 03/10] TEST: Filter legacy SH bases warnings in bootstrap direction getter test Filter legacy SH bases warnings in bootstrap direction getter test. Warnings are already being tested, and there is no further need to ensure that these warnings are raised. Filtering them avoids having them printed in the standard output, and thus, unnecessary clutter. Fixes: ``` direction/tests/test_bootstrap_direction_getter.py::test_bdg_get_direction /home/runner/work/dipy/dipy/venv/lib/python3.11/site-packages/dipy/reconst/shm.py:351: PendingDeprecationWarning: The legacy descoteaux07 SH basis uses absolute values for negative harmonic degrees. It is outdated and will be deprecated in a future DIPY release. Consider using the new descoteaux07 basis by setting the `legacy` parameter to `False`. warn(descoteaux07_legacy_msg, category=PendingDeprecationWarning) ``` Raised for example at: https://github.com/dipy/dipy/actions/runs/6657261609/job/18091595226?pr=2941#step:9:4028 --- .../tests/test_bootstrap_direction_getter.py | 49 +++++++++++++------ 1 file changed, 35 insertions(+), 14 deletions(-) diff --git a/dipy/direction/tests/test_bootstrap_direction_getter.py b/dipy/direction/tests/test_bootstrap_direction_getter.py index ca67e80bd6..6b6e6af1a1 100644 --- a/dipy/direction/tests/test_bootstrap_direction_getter.py +++ b/dipy/direction/tests/test_bootstrap_direction_getter.py @@ -82,30 +82,51 @@ def test_bdg_get_direction(): sphere = get_sphere('symmetric362') response = (np.array([0.0015, 0.0003, 0.0003]), 1) - csd_model = ConstrainedSphericalDeconvModel(gtab, response, sh_order=6) + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", message=shm.descoteaux07_legacy_msg, + category=PendingDeprecationWarning) + csd_model = ConstrainedSphericalDeconvModel(gtab, response, sh_order=6) + point = np.array([0., 0., 0.]) prev_direction = sphere.vertices[5] # test case in which no valid direction is found with default max attempts - boot_dg = BootDirectionGetter(data, model=csd_model, max_angle=10., - sphere=sphere) - npt.assert_equal(boot_dg.get_direction(point, prev_direction), 1) + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", message=shm.descoteaux07_legacy_msg, + category=PendingDeprecationWarning) + boot_dg = BootDirectionGetter(data, model=csd_model, max_angle=10., + sphere=sphere) + npt.assert_equal(boot_dg.get_direction(point, prev_direction), 1) # test case in which no valid direction is found with new max attempts - boot_dg = BootDirectionGetter(data, model=csd_model, max_angle=10, - sphere=sphere, max_attempts=3) - npt.assert_equal(boot_dg.get_direction(point, prev_direction), 1) + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", message=shm.descoteaux07_legacy_msg, + category=PendingDeprecationWarning) + boot_dg = BootDirectionGetter(data, model=csd_model, max_angle=10, + sphere=sphere, max_attempts=3) + npt.assert_equal(boot_dg.get_direction(point, prev_direction), 1) # test case in which a valid direction is found - boot_dg = BootDirectionGetter(data, model=csd_model, max_angle=60., - sphere=sphere, max_attempts=5) - npt.assert_equal(boot_dg.get_direction(point, prev_direction), 0) + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", message=shm.descoteaux07_legacy_msg, + category=PendingDeprecationWarning) + boot_dg = BootDirectionGetter(data, model=csd_model, max_angle=60., + sphere=sphere, max_attempts=5) + npt.assert_equal(boot_dg.get_direction(point, prev_direction), 0) # test invalid max_attempts parameters - npt.assert_raises( - ValueError, - lambda: BootDirectionGetter(data, csd_model, 60, sphere=sphere, - max_attempts=0)) + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", message=shm.descoteaux07_legacy_msg, + category=PendingDeprecationWarning) + npt.assert_raises( + ValueError, + lambda: BootDirectionGetter(data, csd_model, 60, sphere=sphere, + max_attempts=0)) def test_bdg_residual(): From c8c8b18e634e27bb5b7490f47c25c4f832145144 Mon Sep 17 00:00:00 2001 From: Maharshi Gor Date: Tue, 7 Nov 2023 12:34:54 -0500 Subject: [PATCH 04/10] volume slider fix --- dipy/viz/horizon/tab/slice.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/dipy/viz/horizon/tab/slice.py b/dipy/viz/horizon/tab/slice.py index 850008bc0c..c1862dbb2f 100644 --- a/dipy/viz/horizon/tab/slice.py +++ b/dipy/viz/horizon/tab/slice.py @@ -246,7 +246,7 @@ def _change_volume(self, slider): value = int(np.rint(slider.value)) if value != self._volume.selected_value: visible_slices = ( - self._slice_x.obj.value, self._slice_y.selected_value, + self._slice_x.selected_value, self._slice_y.selected_value, self._slice_z.selected_value) valid_vol = self._visualizer.change_volume( self._volume.selected_value, value, @@ -280,10 +280,10 @@ def _change_volume(self, slider): # Updating visibilities slices = [self._slice_x, self._slice_y, self._slice_z] - for s, i in enumerate(slices): + for i, s in enumerate(slices): self._update_slice_visibility( None, - s.obj, + s, self._visualizer.slice_actors[i], s.visibility ) From f19e6db1c8fd3a441d8199aef00810515af16ff1 Mon Sep 17 00:00:00 2001 From: Maharshi Gor Date: Tue, 7 Nov 2023 16:22:53 -0500 Subject: [PATCH 05/10] opacity and slider toggle fixed --- dipy/viz/horizon/tab/base.py | 5 +++++ dipy/viz/horizon/tab/slice.py | 21 +++++++++++++++++++++ 2 files changed, 26 insertions(+) diff --git a/dipy/viz/horizon/tab/base.py b/dipy/viz/horizon/tab/base.py index 3b02df48a2..47d03dddf4 100644 --- a/dipy/viz/horizon/tab/base.py +++ b/dipy/viz/horizon/tab/base.py @@ -122,11 +122,16 @@ def _render_tab_elements(self, tab_id, elements): self._tab_ui.add_element(tab_id, element.obj, element.position) def _tab_selected(self, tab_ui): + if self._active_tab_id == tab_ui.active_tab_idx: + self._active_tab_id = -1 + return + self._active_tab_id = tab_ui.active_tab_idx current_tab = self._tabs[self._active_tab_id] if current_tab.__class__.__name__ == 'SlicesTab': self.tab_changed(current_tab.actors) + current_tab.on_tab_selected() def reposition(self, win_size): """ diff --git a/dipy/viz/horizon/tab/slice.py b/dipy/viz/horizon/tab/slice.py index c1862dbb2f..370d6391aa 100644 --- a/dipy/viz/horizon/tab/slice.py +++ b/dipy/viz/horizon/tab/slice.py @@ -36,6 +36,11 @@ def __init__(self, slices_visualizer, slice_id=0, self.on_slice_change = lambda _tab_id, _x, _y, _z: None + self._opacity_toggle = build_checkbox( + labels=[''], + checked_labels=[''], + on_change=self._toggle_opacity) + self._slice_opacity_label, self._slice_opacity = build_slider( initial_value=1., max_value=1., @@ -167,6 +172,7 @@ def __init__(self, slices_visualizer, slice_id=0, self._visualizer.register_picker_callback(self._change_picked_voxel) self.register_elements( + self._opacity_toggle, self._slice_opacity_label, self._slice_opacity, self._slice_x_toggle, @@ -209,6 +215,15 @@ def _change_opacity(self, slider): self._slice_opacity.selected_value = slider.value self._update_opacities() + def _toggle_opacity(self, checkbox): + if '' in checkbox.checked_labels: + self._slice_opacity.selected_value = 1 + self._slice_opacity.obj.value = 1 + else: + self._slice_opacity.selected_value = 0 + self._slice_opacity.obj.value = 0 + self._update_opacities() + def _change_picked_voxel(self, message): self._voxel_data.obj.message = message self._voxel_data.selected_value = message @@ -321,6 +336,11 @@ def _update_opacities(self): slice_actor.GetProperty().SetOpacity( self._slice_opacity.selected_value) + def on_tab_selected(self): + self._slice_x.obj.set_visibility(self._slice_x.visibility) + self._slice_y.obj.set_visibility(self._slice_y.visibility) + self._slice_z.obj.set_visibility(self._slice_z.visibility) + def update_slices(self, x_slice, y_slice, z_slice): """Updates slicer positions. @@ -359,6 +379,7 @@ def build(self, tab_id, _tab_ui): self._tab_id = tab_id x_pos = .02 + self._opacity_toggle.position = (x_pos, .85) self._slice_x_toggle.position = (x_pos, .62) self._slice_y_toggle.position = (x_pos, .38) self._slice_z_toggle.position = (x_pos, .15) From bae47db23a01278d5def3dbdee7e59c117e47061 Mon Sep 17 00:00:00 2001 From: Maharshi Gor Date: Wed, 8 Nov 2023 14:17:37 -0500 Subject: [PATCH 06/10] bugfix for --force issue --- dipy/workflows/viz.py | 1 + 1 file changed, 1 insertion(+) diff --git a/dipy/workflows/viz.py b/dipy/workflows/viz.py index 18d3d1fd02..3575d7c6c9 100644 --- a/dipy/workflows/viz.py +++ b/dipy/workflows/viz.py @@ -106,6 +106,7 @@ def run(self, input_files, cluster=False, cluster_thr=15., adaptive visualization, Proceedings of: International Society of Magnetic Resonance in Medicine (ISMRM), Montreal, Canada, 2019. """ + super(HorizonFlow, self).__init__(force=True) verbose = True tractograms = [] images = [] From c5a21626cf6e7f7a51882adcd5b970427639879a Mon Sep 17 00:00:00 2001 From: Serge Koudoro Date: Fri, 10 Nov 2023 15:22:41 -0500 Subject: [PATCH 07/10] disable annotations --- .codecov.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.codecov.yml b/.codecov.yml index 6d85a429da..0e7e0e60de 100644 --- a/.codecov.yml +++ b/.codecov.yml @@ -1,3 +1,6 @@ +github_checks: + annotations: false + comment: layout: "reach, diff, files" behavior: default From b32aa21e74ef1c7047cb1d804915a815c87eb5ac Mon Sep 17 00:00:00 2001 From: Serge Koudoro Date: Fri, 10 Nov 2023 15:23:23 -0500 Subject: [PATCH 08/10] [MTN] indentation --- .codecov.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.codecov.yml b/.codecov.yml index 0e7e0e60de..136c265d0a 100644 --- a/.codecov.yml +++ b/.codecov.yml @@ -1,5 +1,5 @@ github_checks: - annotations: false + annotations: false comment: layout: "reach, diff, files" From 529aa01b328d9925fd291538fb8b4818d9ede613 Mon Sep 17 00:00:00 2001 From: Jong Sung Park Date: Fri, 3 Nov 2023 11:31:05 -0400 Subject: [PATCH 09/10] NF: Updating EVAC+ model and adding util function This commit changes the EVAC+ model and its default weights for better brain extraction output. It also adds an option to only select the largest foreground and remove its holes, so that most prediction noises can be removed. Slice wise filling holes were added as well for more desired results. Utils functions have changed to use voxsize, and tests have changed accordingly as well. --- dipy/data/fetcher.py | 8 +-- dipy/nn/evac.py | 44 ++++++++++------ dipy/nn/tests/test_evac.py | 18 +++++-- dipy/nn/tests/test_utils.py | 20 +++++--- dipy/nn/utils.py | 100 +++++++++++++++++++++++++++++++----- 5 files changed, 147 insertions(+), 43 deletions(-) diff --git a/dipy/data/fetcher.py b/dipy/data/fetcher.py index 7ddc10e4f9..7d27d2e49a 100644 --- a/dipy/data/fetcher.py +++ b/dipy/data/fetcher.py @@ -323,18 +323,18 @@ def fetcher(): "fetch_evac_weights", pjoin(dipy_home, 'evac'), 'https://ndownloader.figshare.com/files/', - ['40150867'], + ['43037191'], ['evac_default_weights.h5'], - ['998d5122e0aef1ccf7c0d58e41e978af'], + ['491cfa4f9a2860fad6c19f2b71b918e1'], doc="Download EVAC+ model weights for Park et. al 2022") fetch_evac_test = _make_fetcher( "fetch_evac_test", pjoin(dipy_home, 'evac'), 'https://ndownloader.figshare.com/files/', - ['40150894'], + ['43040074'], ['evac_test_data.npz'], - ['a2173e8f800ab7ab3b159b86bf3a8536'], + ['a5ad7116c9914a53ba891f9e73f3a132'], doc="Download EVAC+ test data for Park et. al 2022") fetch_stanford_t1 = _make_fetcher( diff --git a/dipy/nn/evac.py b/dipy/nn/evac.py index 5ba5a37bc1..515dbddcbd 100644 --- a/dipy/nn/evac.py +++ b/dipy/nn/evac.py @@ -15,13 +15,13 @@ from dipy.io.image import load_nifti from dipy.align.reslice import reslice from dipy.nn.utils import normalize, set_logger_level, transform_img -from dipy.nn.utils import recover_img +from dipy.nn.utils import recover_img, correct_minor_errors tf, have_tf, _ = optional_package('tensorflow') if have_tf: from tensorflow.keras.models import Model from tensorflow.keras.layers import Softmax, Conv3DTranspose - from tensorflow.keras.layers import Conv3D, LayerNormalization, LeakyReLU + from tensorflow.keras.layers import Conv3D, LayerNormalization, ReLU from tensorflow.keras.layers import Concatenate, Layer, Dropout, Add if Version(tf.__version__) < Version('2.0.0'): raise ImportError('Please upgrade to TensorFlow 2+') @@ -99,13 +99,13 @@ def __init__(self, out_channels, kernel_size, strides, padding=padding)) self.layer_list.append(Dropout(drop_r)) self.layer_list.append(LayerNormalization()) - self.layer_list.append(LeakyReLU()) + self.layer_list.append(ReLU()) if layer_type=='down': self.layer_list2.append(Conv3D(1, 2, strides=2, padding='same')) - self.layer_list2.append(LeakyReLU()) + self.layer_list2.append(ReLU()) elif layer_type=='up': self.layer_list2.append(Conv3DTranspose(1, 2, strides=2, padding='same')) - self.layer_list2.append(LeakyReLU()) + self.layer_list2.append(ReLU()) self.channel_sum = ChannelSum() self.add = Add() @@ -303,7 +303,8 @@ def __predict(self, x_test): def predict(self, T1, affine, voxsize=(1, 1, 1), batch_size=None, - return_affine=False, return_prob=False): + return_affine=False, return_prob=False, + largest_area=True): r""" Wrapper function to facilitate prediction of larger dataset. @@ -311,7 +312,7 @@ def predict(self, T1, affine, ---------- T1 : np.ndarray or list of np.ndarrys For a single image, input should be a 3D array. - If multiple images, it should be a 4D array or a list. + If multiple images, it should be a a list or tuple. affine : np.ndarray (4, 4) or (batch, 4, 4) or list of np.ndarrays with len of batch @@ -342,6 +343,11 @@ def predict(self, T1, affine, binary mask. Useful for testing. Default is False + largest_area : bool, optional + Whether to exclude only the largest background/foreground. + Useful for solving minor errors. + Default is True + Returns ------- pred_output : np.ndarray (...) or (batch, ...) @@ -350,6 +356,7 @@ def predict(self, T1, affine, affine : np.ndarray (...) or (batch, ...) affine matrix of mask only if return_affine is True + """ voxsize = np.array(voxsize) @@ -358,7 +365,7 @@ def predict(self, T1, affine, if isinstance(T1, (list, tuple)): dim = 4 T1 = np.array(T1) - elif len(T1.shape)==3: + elif len(T1.shape) == 3: dim = 3 if batch_size is not None: logger.warning('Batch size specified, but not used', @@ -370,24 +377,27 @@ def predict(self, T1, affine, affine = np.expand_dims(affine, 0) voxsize = np.expand_dims(voxsize, 0) else: - raise ValueError("T1 data should be a np.ndarray of dimension 3 or" - "a list/tuple of it") + raise ValueError("T1 data should be a np.ndarray of dimension 3 " + "or a list/tuple of it") input_data = np.zeros((128, 128, 128, len(T1))) rev_affine = np.zeros((len(T1), 4, 4)) + ori_shapes = np.zeros((len(T1), 3)).astype(int) # Normalize the data. n_T1 = np.zeros(T1.shape) for i, T1_img in enumerate(T1): n_T1[i] = normalize(T1_img, new_min=0, new_max=1) - t_img, t_affine = transform_img(n_T1[i], - affine[i]) + t_img, t_affine, ori_shape = transform_img(n_T1[i], + affine[i], + voxsize[i]) input_data[..., i] = t_img rev_affine[i] = t_affine + ori_shapes[i] = ori_shape # Prediction stage prediction = np.zeros((len(T1), 128, 128, 128), - dtype=np.float32) + dtype=np.float32) for batch_idx in range(batch_size, len(T1)+1, batch_size): batch = input_data[..., batch_idx-batch_size:batch_idx] temp_input = prepare_img(batch) @@ -403,9 +413,13 @@ def predict(self, T1, affine, for i, T1_img in enumerate(T1): output = recover_img(prediction[i], rev_affine[i], - T1_img.shape) - if return_prob == False: + voxsize=voxsize[i], + ori_shape=ori_shapes[i], + image_shape=T1_img.shape) + if not return_prob: output = np.where(output >= 0.5, 1, 0) + if largest_area: + output = correct_minor_errors(output) output_mask.append(output) if dim == 3: diff --git a/dipy/nn/tests/test_evac.py b/dipy/nn/tests/test_evac.py index 28f2787cd3..0c791d5117 100644 --- a/dipy/nn/tests/test_evac.py +++ b/dipy/nn/tests/test_evac.py @@ -1,10 +1,11 @@ import pytest from packaging.version import Version +import numpy as np +import numpy.testing as npt + from dipy.data import get_fnames from dipy.utils.optpkg import optional_package -import numpy as np -from numpy.testing import assert_almost_equal tf, have_tf, _ = optional_package('tensorflow') @@ -13,6 +14,7 @@ if Version(tf.__version__) < Version('2.0.0'): raise ImportError('Please upgrade to TensorFlow 2+') + @pytest.mark.skipif(not have_tf, reason='Requires TensorFlow') def test_default_weights(): file_path = get_fnames('evac_test_data') @@ -21,7 +23,7 @@ def test_default_weights(): evac_model = EVACPlus() results_arr = evac_model.predict(input_arr, np.eye(4), return_prob=True) - assert_almost_equal(results_arr, output_arr, decimal=4) + npt.assert_almost_equal(results_arr, output_arr, decimal=4) @pytest.mark.skipif(not have_tf, reason='Requires TensorFlow') @@ -35,4 +37,12 @@ def test_default_weights_batch(): fake_affine = np.array([np.eye(4), np.eye(4)]) fake_voxsize = np.ones((2, 3)) results_arr = evac_model.predict(input_arr, fake_affine, fake_voxsize, batch_size=2, return_prob=True) - assert_almost_equal(results_arr, output_arr, decimal=4) + npt.assert_almost_equal(results_arr, output_arr, decimal=4) + + +@pytest.mark.skipif(not have_tf, reason='Requires TensorFlow') +def test_T1_error(): + T1 = np.ones((3, 32, 32, 32)) + evac_model = EVACPlus() + fake_affine = np.array([np.eye(4), np.eye(4), np.eye(4)]) + npt.assert_raises(ValueError, evac_model.predict, T1, fake_affine) \ No newline at end of file diff --git a/dipy/nn/tests/test_utils.py b/dipy/nn/tests/test_utils.py index fc9a8a34fd..0c78cf035c 100644 --- a/dipy/nn/tests/test_utils.py +++ b/dipy/nn/tests/test_utils.py @@ -1,14 +1,22 @@ import numpy as np from dipy.nn.utils import normalize, unnormalize, transform_img, recover_img +from dipy.testing.decorators import set_random_number_generator -def test_norm(): - temp = np.random.random((8, 8, 8)) * 10 + +@set_random_number_generator() +def test_norm(rng=None): + temp = rng.random((8, 8, 8)) * 10 temp2 = normalize(temp) temp2 = unnormalize(temp2, -1, 1, 0, 10) np.testing.assert_almost_equal(temp, temp2, 1) -def test_transform(): - temp = np.random.random((30, 31, 32)) - temp2, new_affine = transform_img(temp, np.eye(4), (32, 32, 32)) - temp2 = recover_img(temp2, new_affine, temp.shape) + +@set_random_number_generator() +def test_transform(rng=None): + temp = rng.random((30, 31, 32)) + temp2, new_affine, ori_shape = transform_img(temp, np.eye(4), + init_shape=(32, 32, 32), + voxsize=np.ones(3)*2) + temp2 = recover_img(temp2, new_affine, ori_shape, temp.shape, + init_shape=(32, 32, 32), voxsize=np.ones(3)*2) np.testing.assert_almost_equal(np.array(temp.shape), np.array(temp2.shape)) \ No newline at end of file diff --git a/dipy/nn/utils.py b/dipy/nn/utils.py index 306fa6c862..0823d01032 100644 --- a/dipy/nn/utils.py +++ b/dipy/nn/utils.py @@ -1,7 +1,8 @@ import numpy as np -from scipy.ndimage import affine_transform +from scipy.ndimage import affine_transform, label from dipy.align.reslice import reslice + def normalize(image, min_v=None, max_v=None, new_min=-1, new_max=1): r""" normalization function @@ -9,20 +10,20 @@ def normalize(image, min_v=None, max_v=None, new_min=-1, new_max=1): Parameters ---------- image : np.ndarray - min_v : int or float (optional) + min_v : int or float, optional minimum value range for normalization intensities below min_v will be clipped if None it is set to min value of image Default : None - max_v : int or float (optional) + max_v : int or float, optional maximum value range for normalization intensities above max_v will be clipped if None it is set to max value of image Default : None - new_min : int or float (optional) + new_min : int or float, optional new minimum value after normalization Default : 0 - new_max : int or float (optional) + new_max : int or float, optional new maximum value after normalization Default : 1 @@ -37,6 +38,7 @@ def normalize(image, min_v=None, max_v=None, new_min=-1, new_max=1): max_v = np.max(image) return np.interp(image, (min_v, max_v), (new_min, new_max)) + def unnormalize(image, norm_min, norm_max, min_v, max_v): r""" unnormalization function @@ -73,7 +75,9 @@ def set_logger_level(log_level, logger): """ logger.setLevel(level=log_level) -def transform_img(image, affine, init_shape=(256, 256, 256), scale=2): + +def transform_img(image, affine, voxsize=None, + init_shape=(256, 256, 256), scale=2): r""" Function to reshape image as an input to the model @@ -83,10 +87,12 @@ def transform_img(image, affine, init_shape=(256, 256, 256), scale=2): Image to transform to voxelspace affine : np.ndarray Affine matrix provided by the file - init_shape : tuple + voxsize : np.ndarray (3,), optional + Voxel size of the image + init_shape : tuple, optional Initial shape to transform the image to Default is (256, 256, 256) - scale : float + scale : float, optional How much we want to scale the image Default is 2 @@ -94,7 +100,11 @@ def transform_img(image, affine, init_shape=(256, 256, 256), scale=2): ------- transformed_img : np.ndarray """ - affine2 = affine.copy() + if voxsize is not None and np.any(voxsize != np.ones(3)): + image, affine2 = reslice(image, affine, voxsize, (1, 1, 1)) + else: + affine2 = affine.copy() + ori_shape = np.array(image.shape) affine2[:3, 3] += np.array([init_shape[0]//2, init_shape[1]//2, init_shape[2]//2]) @@ -102,9 +112,11 @@ def transform_img(image, affine, init_shape=(256, 256, 256), scale=2): transformed_img = affine_transform(image, inv_affine, output_shape=init_shape) transformed_img, _ = reslice(transformed_img, np.eye(4), (1, 1, 1), (scale, scale, scale)) - return transformed_img, affine2 + return transformed_img, affine2, ori_shape -def recover_img(image, affine, ori_shape, scale=2): + +def recover_img(image, affine, ori_shape, image_shape, + init_shape=(256, 256, 256), voxsize=None, scale=2): r""" Function to recover image back to its original shape @@ -114,9 +126,16 @@ def recover_img(image, affine, ori_shape, scale=2): Image to recover affine : np.ndarray Affine matrix provided from transform_img - ori_shape : tuple - Original shape of image - scale : float + ori_shape : np.ndarray (3,) + Original shape of isotropic image + image_shape : tuple (3,) + Original shape of actual image + init_shape : tuple (3,), optional + Initial shape to transform the image to + Default is (256, 256, 256) + voxsize : np.ndarray (3,), optional + Voxel size of the original image + scale : float, optional Scale that was used in transform_img Default is 2 @@ -126,4 +145,57 @@ def recover_img(image, affine, ori_shape, scale=2): """ new_image, _ = reslice(image, np.eye(4), (scale, scale, scale), (1, 1, 1)) recovered_img = affine_transform(new_image, affine, output_shape=ori_shape) + affine[:3, 3] += np.array([init_shape[0]//2, + init_shape[1]//2, + init_shape[2]//2]) + if voxsize is not None and np.any(voxsize != np.ones(3)): + kwargs = {'order': 1, + 'mode': 'constant', + 'matrix': voxsize, + 'output_shape': image_shape} + recovered_img = np.round(affine_transform(recovered_img, **kwargs)) return recovered_img + + +def correct_minor_errors(binary_img): + """ + Remove any small mask chunks or holes + that could be in the segmentation output. + + Parameters + ---------- + binary_img : np.ndarray + Binary image + + Returns + ------- + largest_img : np.ndarray + """ + largest_img = np.zeros_like(binary_img) + chunks, n_chunk = label(np.abs(1-binary_img)) + u, c = np.unique(chunks[chunks != 0], return_counts=True) + target = u[np.argmax(c)] + largest_img = np.where(chunks == target, 0, 1) + + chunks, n_chunk = label(largest_img) + u, c = np.unique(chunks[chunks != 0], return_counts=True) + target = u[np.argmax(c)] + largest_img = np.where(chunks == target, 1, 0) + + for x in range(largest_img.shape[0]): + chunks, n_chunk = label(np.abs(1-largest_img[x])) + u, c = np.unique(chunks[chunks != 0], return_counts=True) + target = u[np.argmax(c)] + largest_img[x] = np.where(chunks == target, 0, 1) + for y in range(largest_img.shape[1]): + chunks, n_chunk = label(np.abs(1-largest_img[:, y])) + u, c = np.unique(chunks[chunks != 0], return_counts=True) + target = u[np.argmax(c)] + largest_img[:, y] = np.where(chunks == target, 0, 1) + for z in range(largest_img.shape[2]): + chunks, n_chunk = label(np.abs(1-largest_img[..., z])) + u, c = np.unique(chunks[chunks != 0], return_counts=True) + target = u[np.argmax(c)] + largest_img[..., z] = np.where(chunks == target, 0, 1) + + return largest_img From abbfbb618a71b1eb596b934ef52df41b8585f8be Mon Sep 17 00:00:00 2001 From: pjsjongsung Date: Fri, 3 Nov 2023 20:39:21 -0400 Subject: [PATCH 10/10] RF: Moving to numpy.random.Generator from numpy.random This commit is to move from numpy.random to numpy.random.Generator as numpy suggests. All tests that require random numbers are now using set_random_number_generator decorator from dipy.testing. Some functions were changed to accept rng (or replaced seed with rng) as an argument so that this complies with the changes made to the test. Randomstates were changed to Generators as well. Some seeds were modified for the test to work. test_mrf.py's global parameters are now all wrapped in a function. --- dipy/align/streamlinear.py | 12 +-- dipy/align/streamwarp.py | 4 +- dipy/align/tests/test_api.py | 8 +- dipy/align/tests/test_expectmax.py | 4 +- dipy/align/tests/test_imaffine.py | 22 ++--- dipy/align/tests/test_parzenhist.py | 37 ++++++--- dipy/align/tests/test_streamlinear.py | 42 +++++----- dipy/align/tests/test_transforms.py | 16 ++-- dipy/boots/resampling.py | 10 ++- dipy/core/gradients.py | 11 ++- dipy/core/sphere_stats.py | 5 +- dipy/core/tests/test_geometry.py | 25 +++--- dipy/core/tests/test_gradients.py | 20 +++-- dipy/core/tests/test_interpolation.py | 5 +- dipy/core/tests/test_optimize.py | 25 +++--- dipy/denoise/enhancement_kernel.pyx | 6 +- dipy/denoise/tests/test_ascm.py | 16 ++-- dipy/denoise/tests/test_lpca.py | 75 ++++++++++------- dipy/denoise/tests/test_nlmeans.py | 20 +++-- dipy/denoise/tests/test_noise_estimate.py | 7 +- dipy/denoise/tests/test_non_local_means.py | 16 ++-- dipy/denoise/tests/test_patch2self.py | 10 +-- .../tests/test_bootstrap_direction_getter.py | 6 +- dipy/direction/tests/test_peaks.py | 17 ++-- dipy/direction/tests/test_pmf.py | 8 +- dipy/io/tests/test_io_peaks.py | 11 ++- dipy/io/tests/test_utils.py | 18 +++-- dipy/nn/tests/test_cnn_1denoiser.py | 29 ++++--- dipy/reconst/tests/test_cross_validation.py | 5 +- dipy/reconst/tests/test_csdeconv.py | 11 ++- dipy/reconst/tests/test_cti.py | 5 +- dipy/reconst/tests/test_dti.py | 12 ++- dipy/reconst/tests/test_eudx_dg.py | 11 +-- dipy/reconst/tests/test_mcsd.py | 23 +++--- dipy/reconst/tests/test_multi_voxel.py | 6 +- dipy/reconst/tests/test_qti.py | 25 +++--- dipy/reconst/tests/test_shore_odf.py | 7 +- dipy/segment/bundles.py | 11 +-- dipy/segment/clustering.py | 6 +- dipy/segment/tests/test_bundles.py | 12 +-- dipy/segment/tests/test_clustering.py | 80 +++++++++++-------- dipy/segment/tests/test_feature.py | 19 +++-- dipy/segment/tests/test_metric.py | 28 ++++--- dipy/segment/tests/test_mrf.py | 75 ++++++++++------- dipy/segment/tests/test_qbx.py | 16 ++-- dipy/segment/tests/test_quickbundles.py | 4 +- dipy/testing/decorators.py | 9 ++- dipy/tracking/mesh.py | 8 +- dipy/tracking/stopping_criterion.pyx | 5 +- dipy/tracking/tests/test_distances.py | 12 ++- dipy/tracking/tests/test_metrics.py | 10 ++- .../tracking/tests/test_stopping_criterion.py | 22 ++--- dipy/tracking/tests/test_streamline.py | 63 ++++++++------- dipy/tracking/tests/test_tracking.py | 9 ++- dipy/tracking/tests/test_utils.py | 23 +++--- dipy/tracking/utils.py | 32 ++++---- dipy/utils/tests/test_arrfuncs.py | 6 +- dipy/viz/tests/test_apps.py | 11 +-- dipy/viz/tests/test_fury.py | 16 ++-- dipy/viz/tests/test_regtools.py | 8 +- dipy/workflows/stats.py | 2 +- dipy/workflows/tests/test_align.py | 8 +- dipy/workflows/tests/test_denoise.py | 6 +- dipy/workflows/tests/test_stats.py | 18 +++-- dipy/workflows/tests/test_viz.py | 10 ++- doc/examples/bundle_assignment_maps.py | 4 +- doc/examples/bundle_shape_similarity.py | 2 +- doc/examples/contextual_enhancement.py | 5 +- doc/examples/fiber_to_bundle_coherence.py | 2 +- doc/examples/gradients_spheres.py | 5 +- doc/examples/kfold_xval.py | 10 ++- doc/examples/reconst_msdki.py | 6 +- doc/examples/reconst_sh.py | 5 +- doc/examples/simulate_multi_tensor.py | 6 +- doc/examples/streamline_length.py | 26 +++--- doc/examples/viz_bundles.py | 8 +- 76 files changed, 709 insertions(+), 489 deletions(-) diff --git a/dipy/align/streamlinear.py b/dipy/align/streamlinear.py index d30fbb8695..44dc4421bf 100644 --- a/dipy/align/streamlinear.py +++ b/dipy/align/streamlinear.py @@ -1042,8 +1042,8 @@ def slr_with_qbx(static, moving, progressive : boolean, optional (default True) - rng : RandomState - If None creates RandomState in function. + rng : np.random.Generator + If None creates random generator in function. num_threads : int, optional Number of threads to be used for OpenMP parallelization. If None @@ -1072,7 +1072,7 @@ def slr_with_qbx(static, moving, """ if rng is None: - rng = np.random.RandomState() + rng = np.random.default_rng() if verbose: logger.info('Static streamlines size {}'.format(len(static))) @@ -1207,8 +1207,8 @@ def groupwise_slr(bundles, x0='affine', tol=0, max_iter=20, qbx_thr=[4], verbose : bool, optional If True, logs information. Default: False. - rng : RandomState - If None, creates RandomState in function. Default: None. + rng : np.random.Generator + If None, creates random generator in function. Default: None. References ---------- @@ -1233,7 +1233,7 @@ def group_distance(bundles, n_bundle): return d if rng is None: - rng = np.random.RandomState() + rng = np.random.default_rng() metric = JointBundleMinDistanceMetric() diff --git a/dipy/align/streamwarp.py b/dipy/align/streamwarp.py index 2aba246bdb..d3731f1b1a 100755 --- a/dipy/align/streamwarp.py +++ b/dipy/align/streamwarp.py @@ -265,7 +265,9 @@ def bundlewarp_shape_analysis(moving_aligned, deformed_bundle, no_disks=10, indx = assignment_map(deformed_bundle, deformed_bundle, n) indx = np.array(indx) - colors = np.random.rand(n, 3) + rng = np.random.default_rng() + + colors = rng.random((n, 3)) disks_color = [] for _, ind in enumerate(indx): diff --git a/dipy/align/tests/test_api.py b/dipy/align/tests/test_api.py index e2e667922a..7f313b7feb 100644 --- a/dipy/align/tests/test_api.py +++ b/dipy/align/tests/test_api.py @@ -23,6 +23,7 @@ from dipy.io.stateful_tractogram import StatefulTractogram, Space from dipy.io.gradients import read_bvals_bvecs from dipy.io.image import load_nifti +from dipy.testing.decorators import set_random_number_generator def setup_module(): @@ -246,11 +247,12 @@ def test_register_dwi_series_and_motion_correction(): npt.assert_array_equal(reg_affines, reg_affines_2) -def test_streamline_registration(): +@set_random_number_generator() +def test_streamline_registration(rng): sl1 = [np.array([[0, 0, 0], [0, 0, 0.5], [0, 0, 1], [0, 0, 1.5]]), np.array([[0, 0, 0], [0, 0.5, 0.5], [0, 1, 1]])] affine_mat = np.eye(4) - affine_mat[:3, 3] = np.random.randn(3) + affine_mat[:3, 3] = rng.standard_normal(3) sl2 = list(transform_tracking_output(sl1, affine_mat)) aligned, matrix = streamline_registration(sl2, sl1) npt.assert_almost_equal(matrix, np.linalg.inv(affine_mat)) @@ -259,7 +261,7 @@ def test_streamline_registration(): # We assume the two tracks come from the same space, but it might have # some affine associated with it: - base_aff = np.eye(4) * np.random.rand() + base_aff = np.eye(4) * rng.random() base_aff[:3, 3] = np.array([1, 2, 3]) base_aff[3, 3] = 1 diff --git a/dipy/align/tests/test_expectmax.py b/dipy/align/tests/test_expectmax.py index 1e24f7254d..bb3705f937 100644 --- a/dipy/align/tests/test_expectmax.py +++ b/dipy/align/tests/test_expectmax.py @@ -309,7 +309,7 @@ def test_quantize_positive_2d(rng): # make sure additive noise doesn't change the quantization result noise_amplitude = np.min([delta / 4.0, min_positive / 4.0]) sz = np.size(true_quantization) - noise = np.random.ranf(sz).reshape(img_shape) * noise_amplitude + noise = rng.random(sz).reshape(img_shape) * noise_amplitude noise = noise.astype(floating) input_image = np.ndarray(img_shape, dtype=floating) # assign intensities plus noise @@ -372,7 +372,7 @@ def test_quantize_positive_3d(rng): # make sure additive noise doesn't change the quantization result noise_amplitude = np.min([delta / 4.0, min_positive / 4.0]) sz = np.size(true_quantization) - noise = np.random.ranf(sz).reshape(img_shape) * noise_amplitude + noise = rng.random(sz).reshape(img_shape) * noise_amplitude noise = noise.astype(floating) input_image = np.ndarray(img_shape, dtype=floating) # assign intensities plus noise diff --git a/dipy/align/tests/test_imaffine.py b/dipy/align/tests/test_imaffine.py index e530f5046e..aaf608be5d 100644 --- a/dipy/align/tests/test_imaffine.py +++ b/dipy/align/tests/test_imaffine.py @@ -175,7 +175,8 @@ def test_transform_origins_3d(): assert_array_almost_equal(actual.affine, expected) -def test_affreg_all_transforms(): +@set_random_number_generator(202311) +def test_affreg_all_transforms(rng): # Test affine registration using all transforms with typical settings # Make sure dictionary entries are processed in the same order regardless @@ -198,7 +199,8 @@ def test_affreg_all_transforms(): trans, factor, nslices, - 1.0) + 1.0, + rng=rng) # Sum of absolute differences start_sad = np.abs(static - moving).sum() metric = imaffine.MutualInformationMetric(32, sampling_pc) @@ -243,7 +245,8 @@ def test_affreg_all_transforms(): np.zeros_like(smask), np.zeros_like(mmask)) -def test_affreg_defaults(): +@set_random_number_generator(202311) +def test_affreg_defaults(rng): # Test all default arguments with an arbitrary transform # Select an arbitrary transform (all of them are already tested # in test_affreg_all_transforms) @@ -260,7 +263,7 @@ def test_affreg_defaults(): factor = factors[ttype][0] transform = regtransforms[ttype] static, moving, static_grid2world, moving_grid2world, smask, mmask, T = \ - setup_random_transform(transform, factor, nslices, 1.0) + setup_random_transform(transform, factor, nslices, 1.0, rng=rng) # Sum of absolute differences start_sad = np.abs(static - moving).sum() @@ -297,7 +300,7 @@ def test_affreg_defaults(): end_sad = np.abs(moving - transformed_inv).sum() reduction = 1 - end_sad / start_sad print("%s>>%f" % (ttype, reduction)) - assert(reduction > 0.9) + assert(reduction > 0.89) @set_random_number_generator(2022966) @@ -326,7 +329,7 @@ def test_mi_gradient(rng): start.param_to_matrix(0.25 * rng.standard_normal(nrot)) # Get data (pair of images related to each other by an known transform) static, moving, static_g2w, moving_g2w, smask, mmask, M = \ - setup_random_transform(transform, factor, nslices, 2.0) + setup_random_transform(transform, factor, nslices, 2.0, rng=rng) # Prepare a MutualInformationMetric instance mi_metric = imaffine.MutualInformationMetric(32, sampling_proportion) @@ -608,10 +611,11 @@ def test_affine_map(rng): assert_raises(AffineInversionError, AffineMap, mat_large_dim) -def test_MIMetric_invalid_params(): +@set_random_number_generator() +def test_MIMetric_invalid_params(rng): transform = regtransforms[('AFFINE', 3)] - static = np.random.rand(20, 20, 20) - moving = np.random.rand(20, 20, 20) + static = rng.random((20, 20, 20)) + moving = rng.random((20, 20, 20)) n = transform.get_number_of_parameters() sampling_proportion = 0.3 theta_sing = np.zeros(n) diff --git a/dipy/align/tests/test_parzenhist.py b/dipy/align/tests/test_parzenhist.py index d7b7da1782..ab86bd63c5 100644 --- a/dipy/align/tests/test_parzenhist.py +++ b/dipy/align/tests/test_parzenhist.py @@ -31,7 +31,7 @@ ('AFFINE', 3): 0.1} -def create_random_image_pair(sh, nvals, seed): +def create_random_image_pair(sh, nvals, rng=None): r""" Create a pair of images with an arbitrary, non-uniform joint PDF Parameters @@ -42,6 +42,9 @@ def create_random_image_pair(sh, nvals, seed): maximum number of different values in the generated 2D images. The voxel intensities of the returned images will be in {0, 1, ..., nvals-1} + rng : numpy.random.generator + numpy's random generator. If None, it is set with a random seed. + Default is None Returns ------- @@ -50,7 +53,8 @@ def create_random_image_pair(sh, nvals, seed): moving : array, shape=sh second image in the image pair """ - rng = np.random.default_rng(seed) + if rng is None: + rng = np.random.default_rng() sz = reduce(mul, sh, 1) sh = tuple(sh) static = rng.integers(0, nvals, sz).reshape(sh) @@ -173,10 +177,10 @@ def test_parzen_joint_histogram(): assert_equal(index, P.padding + i) -def test_parzen_densities(): +@set_random_number_generator(1246592) +def test_parzen_densities(rng): # Test the computation of the joint intensity distribution # using a dense and a sparse set of values - seed = 1246592 nbins = 32 nr = 30 nc = 35 @@ -186,10 +190,10 @@ def test_parzen_densities(): for dim in [2, 3]: if dim == 2: shape = (nr, nc) - static, moving = create_random_image_pair(shape, nvals, seed) + static, moving = create_random_image_pair(shape, nvals, rng) else: shape = (ns, nr, nc) - static, moving = create_random_image_pair(shape, nvals, seed) + static, moving = create_random_image_pair(shape, nvals, rng) # Initialize parzen_hist = ParzenJointHistogram(nbins) @@ -257,7 +261,6 @@ def test_parzen_densities(): expected_smarginal_sparse) -@set_random_number_generator(3147702) def setup_random_transform(transform, rfactor, nslices=45, sigma=1, rng=None): r""" Creates a pair of images related to each other by an affine transform @@ -276,7 +279,14 @@ def setup_random_transform(transform, rfactor, nslices=45, sigma=1, rng=None): added to the identity parameters to create the random transform nslices: int number of slices to be stacked to form the volumes + sigma: float + standard deviation of the gaussian filter + rng: np.random.Generator + numpy's random generator. If None, it is set with a random seed. + Default is None """ + if rng is None: + rng = np.random.default_rng() dim = 2 if nslices == 1 else 3 if transform.get_dim() != dim: raise ValueError("Transform and requested volume have different dims.") @@ -316,7 +326,8 @@ def setup_random_transform(transform, rfactor, nslices=45, sigma=1, rng=None): return static, moving, static_g2w, moving_g2w, smask, mmask, M -def test_joint_pdf_gradients_dense(): +@set_random_number_generator(3147702) +def test_joint_pdf_gradients_dense(rng): # Compare the analytical and numerical (finite differences) gradient of # the joint distribution (i.e. derivatives of each histogram cell) w.r.t. # the transform parameters. Since the histograms are discrete partitions @@ -351,7 +362,8 @@ def test_joint_pdf_gradients_dense(): theta = transform.get_identity_parameters() static, moving, static_g2w, moving_g2w, smask, mmask, M = \ - setup_random_transform(transform, factor, nslices, 5.0) + setup_random_transform(transform, factor, nslices, 5.0, + rng=rng) parzen_hist = ParzenJointHistogram(32) parzen_hist.setup(static, moving, smask, mmask) @@ -412,7 +424,8 @@ def test_joint_pdf_gradients_dense(): assert(std_cosine < 0.25) -def test_joint_pdf_gradients_sparse(): +@set_random_number_generator(3147702) +def test_joint_pdf_gradients_sparse(rng): h = 1e-4 # Make sure dictionary entries are processed in the same order regardless @@ -435,14 +448,14 @@ def test_joint_pdf_gradients_sparse(): theta = transform.get_identity_parameters() static, moving, static_g2w, moving_g2w, smask, mmask, M = \ - setup_random_transform(transform, factor, nslices, 5.0) + setup_random_transform(transform, factor, nslices, 5.0, + rng=rng) parzen_hist = ParzenJointHistogram(32) parzen_hist.setup(static, moving, smask, mmask) # Sample the fixed-image domain k = 3 sigma = 0.25 - rng = np.random.default_rng(1234) shape = np.array(static.shape, dtype=np.int32) samples = sample_domain_regular(k, shape, static_g2w, sigma, rng) samples = np.array(samples) diff --git a/dipy/align/tests/test_streamlinear.py b/dipy/align/tests/test_streamlinear.py index 62f035159b..369885939b 100644 --- a/dipy/align/tests/test_streamlinear.py +++ b/dipy/align/tests/test_streamlinear.py @@ -27,6 +27,7 @@ from dipy.align.bundlemin import (_bundle_minimum_distance_matrix, _bundle_minimum_distance, distance_matrix_mdf) +from dipy.testing.decorators import set_random_number_generator def simulated_bundle(no_streamlines=10, waves=False, no_pts=12): @@ -198,12 +199,12 @@ def test_min_vs_min_fast_precision(): assert_equal(bmd.distance(x_test), bmdf.distance(x_test)) -def test_same_number_of_points(): - - A = [np.random.rand(10, 3), np.random.rand(20, 3)] - B = [np.random.rand(21, 3), np.random.rand(30, 3)] - C = [np.random.rand(10, 3), np.random.rand(10, 3)] - D = [np.random.rand(20, 3), np.random.rand(20, 3)] +@set_random_number_generator() +def test_same_number_of_points(rng): + A = [rng.random((10, 3)), rng.random((20, 3))] + B = [rng.random((21, 3)), rng.random((30, 3))] + C = [rng.random((10, 3)), rng.random((10, 3))] + D = [rng.random((20, 3)), rng.random((20, 3))] slr = StreamlineLinearRegistration() assert_raises(ValueError, slr.optimize, A, B) @@ -254,14 +255,14 @@ def test_efficient_bmd(): assert_almost_equal(dist, dist2) -def test_openmp_locks(): - +@set_random_number_generator() +def test_openmp_locks(rng): static = [] moving = [] pts = 20 for i in range(1000): - s = np.random.rand(pts, 3) + s = rng.random((pts, 3)) static.append(s) moving.append(s + 2) @@ -413,18 +414,18 @@ def test_vectorize_streamlines(): assert_equal(np.all(cb_subj1_pts_no == 10), True) -def test_x0_input(): - +@set_random_number_generator() +def test_x0_input(rng): for x0 in [6, 7, 12, "Rigid", 'rigid', "similarity", "Affine"]: StreamlineLinearRegistration(x0=x0) - for x0 in [np.random.rand(6), np.random.rand(7), np.random.rand(12)]: + for x0 in [rng.random(6), rng.random(7), rng.random(12)]: StreamlineLinearRegistration(x0=x0) - for x0 in [8, 20, "Whatever", np.random.rand(20), np.random.rand(20, 3)]: + for x0 in [8, 20, "Whatever", rng.random(20), rng.random((20, 3))]: assert_raises(ValueError, StreamlineLinearRegistration, x0=x0) - x0 = np.random.rand(4, 3) + x0 = rng.random((4, 3)) assert_raises(ValueError, StreamlineLinearRegistration, x0=x0) x0_6 = np.zeros(6) @@ -438,10 +439,10 @@ def test_x0_input(): assert_equal(slr.x0, x0_s[i]) -def test_compose_decompose_matrix44(): - +@set_random_number_generator() +def test_compose_decompose_matrix44(rng): for i in range(20): - x0 = np.random.rand(12) + x0 = rng.random(12) mat = compose_matrix44(x0[:6]) assert_array_almost_equal(x0[:6], decompose_matrix44(mat, size=6)) mat = compose_matrix44(x0[:7]) @@ -482,9 +483,10 @@ def test_cascade_of_optimizations_and_threading(): assert_(slm3.fopt < slm2.fopt) -def test_wrong_num_threads(): - A = [np.random.rand(10, 3), np.random.rand(10, 3)] - B = [np.random.rand(10, 3), np.random.rand(10, 3)] +@set_random_number_generator() +def test_wrong_num_threads(rng): + A = [rng.random((10, 3)), rng.random((10, 3))] + B = [rng.random((10, 3)), rng.random((10, 3))] slr = StreamlineLinearRegistration(num_threads=0) assert_raises(ValueError, slr.optimize, A, B) diff --git a/dipy/align/tests/test_transforms.py b/dipy/align/tests/test_transforms.py index db021ef7d5..00f6634f97 100644 --- a/dipy/align/tests/test_transforms.py +++ b/dipy/align/tests/test_transforms.py @@ -1,10 +1,12 @@ -from dipy.align.transforms import regtransforms, Transform import numpy as np from numpy.testing import (assert_array_equal, assert_array_almost_equal, assert_equal, assert_raises) +from dipy.align.transforms import regtransforms, Transform +from dipy.testing.decorators import set_random_number_generator + def test_number_of_parameters(): expected_params = {('TRANSLATION', 2): 2, @@ -28,8 +30,8 @@ def test_number_of_parameters(): expected_params[ttype]) -def test_param_to_matrix_2d(): - rng = np.random.RandomState() +@set_random_number_generator() +def test_param_to_matrix_2d(rng): # Test translation matrix 2D transform = regtransforms[('TRANSLATION', 2)] dx, dy = rng.uniform(size=(2,)) @@ -107,8 +109,8 @@ def test_param_to_matrix_2d(): assert_raises(ValueError, transform.param_to_matrix, theta) -def test_param_to_matrix_3d(): - rng = np.random.RandomState() +@set_random_number_generator() +def test_param_to_matrix_3d(rng): # Test translation matrix 3D transform = regtransforms[('TRANSLATION', 3)] dx, dy, dz = rng.uniform(size=(3,)) @@ -222,8 +224,8 @@ def test_identity_parameters(): assert_array_almost_equal(actual, expected) -def test_jacobian_functions(): - rng = np.random.RandomState() +@set_random_number_generator() +def test_jacobian_functions(rng): # Compare the analytical Jacobians with their numerical approximations h = 1e-8 nsamples = 50 diff --git a/dipy/boots/resampling.py b/dipy/boots/resampling.py index a83b046c13..216f244965 100644 --- a/dipy/boots/resampling.py +++ b/dipy/boots/resampling.py @@ -72,9 +72,11 @@ def bootstrap(x, statistic=bs_se, B=1000, alpha=0.95): N = len(x) bs_pdf = np.empty((B,)) + rng = np.random.default_rng() + for ii in range(0, B): # resample with replacement - rand_index = np.int16(np.round(np.random.random(N) * (N - 1))) + rand_index = np.int16(np.round(rng.random(N) * (N - 1))) bs_pdf[ii] = statistic(x[rand_index]) return bs_pdf, bs_se(bs_pdf), abc(x, statistic, alpha=alpha) @@ -272,11 +274,13 @@ def jackknife(pdf, statistic=np.std, M=None): M = np.minimum(M, N - 1) jk_pdf = np.empty((M,)) + rng = np.random.default_rng() + for ii in range(0, M): - rand_index = np.round(np.random.random() * (N - 1)) + rand_index = np.round(rng.random() * (N - 1)) # choose a unique random sample to remove while pdf_mask[rand_index] == 0: - rand_index = np.round(np.random.random() * (N - 1)) + rand_index = np.round(rng.random() * (N - 1)) # set mask to zero for chosen random index so not to choose again pdf_mask[rand_index] = 0 mask_index[rand_index] = 0 diff --git a/dipy/core/gradients.py b/dipy/core/gradients.py index ae0972ecc4..b7d80171f6 100644 --- a/dipy/core/gradients.py +++ b/dipy/core/gradients.py @@ -704,7 +704,7 @@ def reorient_bvecs(gtab, affines, atol=1e-2): b0_threshold=gtab.b0_threshold, atol=atol) -def generate_bvecs(N, iters=5000): +def generate_bvecs(N, iters=5000, rng=None): """Generates N bvectors. Uses dipy.core.sphere.disperse_charges to model electrostatic repulsion on @@ -717,6 +717,9 @@ def generate_bvecs(N, iters=5000): of bvals used. iters : int Number of iterations to run. + rng : numpy.random.Generator, optional + Numpy's random number generator. If None, the generator is created. + Default is None. Returns ------- @@ -724,8 +727,10 @@ def generate_bvecs(N, iters=5000): The generated directions, represented as a unit vector, of each gradient. """ - theta = np.pi * np.random.rand(N) - phi = 2 * np.pi * np.random.rand(N) + if rng is None: + rng = np.random.default_rng() + theta = np.pi * rng.random(N) + phi = 2 * np.pi * rng.random(N) hsph_initial = HemiSphere(theta=theta, phi=phi) hsph_updated, potential = disperse_charges(hsph_initial, iters) bvecs = hsph_updated.vertices diff --git a/dipy/core/sphere_stats.py b/dipy/core/sphere_stats.py index 209a1cf5dd..adbff95342 100644 --- a/dipy/core/sphere_stats.py +++ b/dipy/core/sphere_stats.py @@ -43,9 +43,10 @@ def random_uniform_on_sphere(n=1, coords='xyz'): >>> X.shape == (4, 3) True """ - z = np.random.uniform(-1, 1, n) + rng = np.random.default_rng() + z = rng.uniform(-1, 1, n) theta = np.arccos(z) - phi = np.random.uniform(0, 2*np.pi, n) + phi = rng.uniform(0, 2*np.pi, n) if coords == 'xyz': r = np.ones(n) return np.vstack(geometry.sphere2cart(r, theta, phi)).T diff --git a/dipy/core/tests/test_geometry.py b/dipy/core/tests/test_geometry.py index 03174875a3..f913217a93 100644 --- a/dipy/core/tests/test_geometry.py +++ b/dipy/core/tests/test_geometry.py @@ -26,6 +26,7 @@ assert_equal, assert_raises, assert_almost_equal,) from dipy.testing.spherepoints import sphere_points +from dipy.testing.decorators import set_random_number_generator from dipy.core.sphere_stats import random_uniform_on_sphere from itertools import permutations @@ -149,7 +150,8 @@ def test_sphere_distance(): assert_raises(ValueError, sphere_distance, [1, 0], [0, 1], 2.0) -def test_vector_cosine(): +@set_random_number_generator() +def test_vector_cosine(rng): a = [0, 1] b = [1, 0] assert_array_almost_equal(vector_cosine(a, b), 0) @@ -164,8 +166,8 @@ def test_vector_cosine(): assert_array_almost_equal(vector_cosine(pts1, pts2), [-1, 1]) # test relationship with correlation # not the same if non-zero vector mean - a = np.random.uniform(size=(100,)) - b = np.random.uniform(size=(100,)) + a = rng.uniform(size=(100,)) + b = rng.uniform(size=(100,)) cc = np.corrcoef(a, b)[0, 1] vcos = vector_cosine(a, b) assert not np.allclose(cc, vcos) @@ -267,12 +269,12 @@ def test_compose_transformations(): assert_raises(ValueError, compose_transformations, A) -def test_compose_decompose_matrix(): - - for translate in permutations(40 * np.random.rand(3), 3): - for angles in permutations(np.deg2rad(90 * np.random.rand(3)), 3): - for shears in permutations(3 * np.random.rand(3), 3): - for scale in permutations(3 * np.random.rand(3), 3): +@set_random_number_generator() +def test_compose_decompose_matrix(rng): + for translate in permutations(40 * rng.random(3), 3): + for angles in permutations(np.deg2rad(90 * rng.random(3)), 3): + for shears in permutations(3 * rng.random(3), 3): + for scale in permutations(3 * rng.random(3), 3): mat = compose_matrix(translate=translate, angles=angles, shear=shears, scale=scale) @@ -333,7 +335,8 @@ def _rotation_from_angles(r): return R -def test_dist_to_corner(): +@set_random_number_generator() +def test_dist_to_corner(rng): affine = np.eye(4) # Calculate the distance with the pythagorean theorem: pythagoras = np.sqrt(np.sum((np.diag(affine)[:-1] / 2) ** 2)) @@ -341,7 +344,7 @@ def test_dist_to_corner(): assert_array_almost_equal(dist_to_corner(affine), pythagoras) # Apply a rotation to the matrix, just to demonstrate the calculation is # robust to that: - R = _rotation_from_angles(np.random.randn(3) * np.pi) + R = _rotation_from_angles(rng.random(3) * np.pi) new_aff = np.vstack([np.dot(R, affine[:3, :]), [0, 0, 0, 1]]) assert_array_almost_equal(dist_to_corner(new_aff), pythagoras) diff --git a/dipy/core/tests/test_gradients.py b/dipy/core/tests/test_gradients.py index ac3a346588..a92719e3f9 100644 --- a/dipy/core/tests/test_gradients.py +++ b/dipy/core/tests/test_gradients.py @@ -20,6 +20,7 @@ from dipy.io.gradients import read_bvals_bvecs from dipy.utils.deprecator import ExpiredDeprecationError from dipy.testing import clear_and_catch_warnings +from dipy.testing.decorators import set_random_number_generator def test_unique_bvals_deprecated(): @@ -354,7 +355,8 @@ def test_qvalues(): npt.assert_almost_equal(bt.qvals, qvals) -def test_reorient_bvecs(): +@set_random_number_generator() +def test_reorient_bvecs(rng): sq2 = np.sqrt(2) / 2 bvals = np.concatenate([[0], np.ones(6) * 1000]) bvecs = np.array([[0, 0, 0], @@ -379,7 +381,7 @@ def test_reorient_bvecs(): rotation_affines = [] rotated_bvecs = bvecs[:] for i in np.where(~gt.b0s_mask)[0]: - rot_ang = np.random.rand() + rot_ang = rng.random() cos_rot = np.cos(rot_ang) sin_rot = np.sin(rot_ang) rotation_affines.append(np.array([[1, 0, 0, 0], @@ -449,16 +451,17 @@ def test_nan_bvecs(): npt.assert_(len(selected_w) == 0) -def test_generate_bvecs(): +@set_random_number_generator() +def test_generate_bvecs(rng): """Tests whether we have properly generated bvecs. """ # Test if the generated b-vectors are unit vectors - bvecs = generate_bvecs(100) + bvecs = generate_bvecs(100, rng=rng) norm = [np.linalg.norm(v) for v in bvecs] npt.assert_almost_equal(norm, np.ones(100)) # Test if two generated vectors are almost orthogonal - bvecs_2 = generate_bvecs(2) + bvecs_2 = generate_bvecs(2, rng=rng) cos_theta = np.dot(bvecs_2[0], bvecs_2[1]) npt.assert_almost_equal(cos_theta, 0., decimal=6) @@ -625,7 +628,8 @@ def test_check_multi_b(): npt.assert_(check_multi_b(gtab, 2, non_zero=False)) -def test_btens_to_params(): +@set_random_number_generator() +def test_btens_to_params(rng): """ Checks if bvals, bdeltas and b_etas are as expected for 4 b-tensor shapes (LTE, PTE, STE, CTE) as well as scaled and rotated versions of them @@ -687,14 +691,14 @@ def test_btens_to_params(): # Test function after rotating+scaling baseline tensors # ----------------------------------------------------- - scales = np.concatenate((np.array([1]), np.random.random(n_scales))) + scales = np.concatenate((np.array([1]), rng.random(n_scales))) for scale in scales: ebs = expected_bvals*scale # Generate `n_rotations` random 3-element vectors of norm 1 - v = np.random.random((n_rotations, 3)) - 0.5 + v = rng.random((n_rotations, 3)) - 0.5 u = np.apply_along_axis(lambda w: w/np.linalg.norm(w), axis=1, arr=v) for rot_idx in range(n_rotations): diff --git a/dipy/core/tests/test_interpolation.py b/dipy/core/tests/test_interpolation.py index 2dbb3e24b3..935737d653 100644 --- a/dipy/core/tests/test_interpolation.py +++ b/dipy/core/tests/test_interpolation.py @@ -20,10 +20,11 @@ from dipy.testing.decorators import set_random_number_generator -def test_trilinear_interpolate(): +@set_random_number_generator() +def test_trilinear_interpolate(rng): """This tests that the trilinear interpolation returns the correct values. """ - a, b, c = np.random.random(3) + a, b, c = rng.random(3) def linear_function(x, y, z): return a * x + b * y + c * z diff --git a/dipy/core/tests/test_optimize.py b/dipy/core/tests/test_optimize.py index 1841a520fe..47f0c8b065 100644 --- a/dipy/core/tests/test_optimize.py +++ b/dipy/core/tests/test_optimize.py @@ -4,6 +4,7 @@ import numpy.testing as npt from dipy.core.optimize import Optimizer, sparse_nnls, spdot import dipy.core.optimize as opt +from dipy.testing.decorators import set_random_number_generator def func(x): @@ -56,7 +57,8 @@ def test_optimize_new_scipy(): npt.assert_array_almost_equal(opt.xopt, np.array([0, 0, 0, 0.])) -def test_sklearn_linear_solver(): +@set_random_number_generator() +def test_sklearn_linear_solver(rng): class SillySolver(opt.SKLearnLinearSolver): def fit(self, X, y): self.coef_ = np.ones(X.shape[-1]) @@ -64,17 +66,18 @@ def fit(self, X, y): MySillySolver = SillySolver() n_samples = 100 n_features = 20 - y = np.random.rand(n_samples) + y = rng.random(n_samples) X = np.ones((n_samples, n_features)) MySillySolver.fit(X, y) npt.assert_equal(MySillySolver.coef_, np.ones(n_features)) npt.assert_equal(MySillySolver.predict(X), np.ones(n_samples) * 20) -def test_nonnegativeleastsquares(): +@set_random_number_generator() +def test_nonnegativeleastsquares(rng): n = 100 X = np.eye(n) - beta = np.random.rand(n) + beta = rng.random(n) y = np.dot(X, beta) my_nnls = opt.NonNegativeLeastSquares() my_nnls.fit(X, y) @@ -82,12 +85,13 @@ def test_nonnegativeleastsquares(): npt.assert_equal(my_nnls.predict(X), y) -def test_spdot(): +@set_random_number_generator() +def test_spdot(rng): n = 100 m = 20 k = 10 - A = np.random.randn(n, m) - B = np.random.randn(m, k) + A = rng.standard_normal((n, m)) + B = rng.standard_normal((m, k)) A_sparse = sps.csr_matrix(A) B_sparse = sps.csr_matrix(B) dense_dot = np.dot(A, B) @@ -98,10 +102,11 @@ def test_spdot(): npt.assert_array_almost_equal(dense_dot, spdot(A_sparse, B)) -def test_sparse_nnls(): +@set_random_number_generator() +def test_sparse_nnls(rng): # Set up the regression: - beta = np.random.rand(10) - X = np.random.randn(1000, 10) + beta = rng.random(10) + X = rng.standard_normal((1000, 10)) y = np.dot(X, beta) beta_hat = sparse_nnls(y, X) beta_hat_sparse = sparse_nnls(y, sps.csr_matrix(X)) diff --git a/dipy/denoise/enhancement_kernel.pyx b/dipy/denoise/enhancement_kernel.pyx index e1f369c8e1..fb5751f118 100644 --- a/dipy/denoise/enhancement_kernel.pyx +++ b/dipy/denoise/enhancement_kernel.pyx @@ -73,6 +73,8 @@ cdef class EnhancementKernel: self.t = t # define a sphere + rng = np.random.default_rng() + if isinstance(orientations, Sphere): # use the sphere defined by the user sphere = orientations @@ -82,8 +84,8 @@ cdef class EnhancementKernel: if n_pts == 0: sphere = None else: - theta = np.pi * np.random.rand(n_pts) - phi = 2 * np.pi * np.random.rand(n_pts) + theta = np.pi * rng.random(n_pts) + phi = 2 * np.pi * rng.random(n_pts) hsph_initial = HemiSphere(theta=theta, phi=phi) sphere, potential = disperse_charges(hsph_initial, 5000) else: diff --git a/dipy/denoise/tests/test_ascm.py b/dipy/denoise/tests/test_ascm.py index 06170e35b6..5faaee7ac1 100644 --- a/dipy/denoise/tests/test_ascm.py +++ b/dipy/denoise/tests/test_ascm.py @@ -7,6 +7,7 @@ from dipy.denoise.non_local_means import non_local_means from dipy.denoise.noise_estimate import estimate_sigma from dipy.denoise.adaptive_soft_matching import adaptive_soft_matching +from dipy.testing.decorators import set_random_number_generator def test_ascm_static(): @@ -19,8 +20,9 @@ def test_ascm_static(): assert_array_almost_equal(S0, S0n) -def test_ascm_random_noise(): - S0 = 100 + 2 * np.random.standard_normal((22, 23, 30)) +@set_random_number_generator() +def test_ascm_random_noise(rng): + S0 = 100 + 2 * rng.standard_normal((22, 23, 30)) S0n1 = non_local_means(S0, sigma=1, rician=False, patch_radius=1, block_radius=1) S0n2 = non_local_means(S0, sigma=1, rician=False, @@ -35,12 +37,13 @@ def test_ascm_random_noise(): assert_equal(np.round(S0n.mean()), 100) -def test_ascm_rmse_with_nlmeans(): +@set_random_number_generator() +def test_ascm_rmse_with_nlmeans(rng): # checks the smoothness S0 = np.ones((30, 30, 30)) * 100 S0[10:20, 10:20, 10:20] = 50 S0[20:30, 20:30, 20:30] = 0 - S0_noise = S0 + 20 * np.random.standard_normal((30, 30, 30)) + S0_noise = S0 + 20 * rng.standard_normal((30, 30, 30)) print("Original RMSE", np.sum(np.abs(S0 - S0_noise)) / np.sum(S0)) S0n1 = non_local_means( @@ -67,12 +70,13 @@ def test_ascm_rmse_with_nlmeans(): assert_(90 < np.mean(S0n) < 110) -def test_sharpness(): +@set_random_number_generator() +def test_sharpness(rng): # check the edge-preserving nature S0 = np.ones((30, 30, 30)) * 100 S0[10:20, 10:20, 10:20] = 50 S0[20:30, 20:30, 20:30] = 0 - S0_noise = S0 + 20 * np.random.standard_normal((30, 30, 30)) + S0_noise = S0 + 20 * rng.standard_normal((30, 30, 30)) S0n1 = non_local_means( S0_noise, sigma=400, diff --git a/dipy/denoise/tests/test_lpca.py b/dipy/denoise/tests/test_lpca.py index da965781f1..6832240d17 100644 --- a/dipy/denoise/tests/test_lpca.py +++ b/dipy/denoise/tests/test_lpca.py @@ -8,6 +8,7 @@ from dipy.denoise.localpca import (localpca, mppca, genpca, _pca_classifier) from dipy.sims.voxel import multi_tensor from dipy.core.gradients import gradient_table, generate_bvecs +from dipy.testing.decorators import set_random_number_generator def setup_module(): @@ -29,7 +30,7 @@ def setup_module(): gtab = gradient_table(bvals, bvecs) -def rfiw_phantom(gtab, snr=None): +def rfiw_phantom(gtab, snr=None, rng=None): """rectangle fiber immersed in water""" # define voxel index slice_ind = np.zeros((10, 10, 8)) @@ -87,9 +88,11 @@ def rfiw_phantom(gtab, snr=None): if snr is None: return DWI else: + if rng is None: + rng = np.random.default_rng() sigma = S2 * 1.0 / snr - n1 = np.random.normal(0, sigma, size=DWI.shape) - n2 = np.random.normal(0, sigma, size=DWI.shape) + n1 = rng.normal(0, sigma, size=DWI.shape) + n2 = rng.normal(0, sigma, size=DWI.shape) return [np.sqrt((DWI / np.sqrt(2) + n1)**2 + (DWI / np.sqrt(2) + n2)**2), sigma] @@ -100,8 +103,9 @@ def test_lpca_static(): assert_array_almost_equal(S0, S0ns) -def test_lpca_random_noise(): - S0 = 100 + 2 * np.random.standard_normal((22, 23, 30, 20)) +@set_random_number_generator() +def test_lpca_random_noise(rng): + S0 = 100 + 2 * rng.standard_normal((22, 23, 30, 20)) S0ns = localpca(S0, sigma=np.std(S0)) assert_(S0ns.min() > S0.min()) @@ -109,11 +113,12 @@ def test_lpca_random_noise(): assert_equal(np.round(S0ns.mean()), 100) -def test_lpca_boundary_behaviour(): +@set_random_number_generator() +def test_lpca_boundary_behaviour(rng): # check is first slice is getting denoised or not ? S0 = 100 * np.ones((20, 20, 20, 20), dtype='f8') S0[:, :, 0, :] = S0[:, :, 0, :] + 2 * \ - np.random.standard_normal((20, 20, 20)) + rng.standard_normal((20, 20, 20)) S0_first = S0[:, :, 0, :] S0ns = localpca(S0, sigma=np.std(S0)) S0ns_first = S0ns[:, :, 0, :] @@ -132,8 +137,9 @@ def test_lpca_boundary_behaviour(): assert_equal(np.round(S0ns_first.mean()), 100) -def test_lpca_rmse(): - S0_w_noise = 100 + 2 * np.random.standard_normal((22, 23, 30, 20)) +@set_random_number_generator() +def test_lpca_rmse(rng): + S0_w_noise = 100 + 2 * rng.standard_normal((22, 23, 30, 20)) rmse_w_noise = np.sqrt(np.mean((S0_w_noise - 100) ** 2)) S0_denoised = localpca(S0_w_noise, sigma=np.std(S0_w_noise)) rmse_denoised = np.sqrt(np.mean((S0_denoised - 100) ** 2)) @@ -141,11 +147,12 @@ def test_lpca_rmse(): assert_(rmse_denoised < rmse_w_noise) -def test_lpca_sharpness(): +@set_random_number_generator() +def test_lpca_sharpness(rng): S0 = np.ones((30, 30, 30, 20), dtype=np.float64) * 100 S0[10:20, 10:20, 10:20, :] = 50 S0[20:30, 20:30, 20:30, :] = 0 - S0 = S0 + 20 * np.random.standard_normal((30, 30, 30, 20)) + S0 = S0 + 20 * rng.standard_normal((30, 30, 30, 20)) S0ns = localpca(S0, sigma=20.0) # check the edge gradient edgs = np.abs(np.mean(S0ns[8, 10:20, 10:20] - S0ns[12, 10:20, 10:20]) - 50) @@ -183,9 +190,10 @@ def test_lpca_wrong(): assert_raises(ValueError, localpca, S0, sigma=1) -def test_phantom(): - DWI_clean = rfiw_phantom(gtab, snr=None) - DWI, sigma = rfiw_phantom(gtab, snr=30) +@set_random_number_generator() +def test_phantom(rng): + DWI_clean = rfiw_phantom(gtab, snr=None, rng=rng) + DWI, sigma = rfiw_phantom(gtab, snr=30, rng=rng) # To test without Rician correction temp = (DWI_clean / sigma)**2 DWI_clean_wrc = (sigma * np.sqrt(np.pi / 2) * np.exp(-0.5 * temp) * @@ -239,33 +247,38 @@ def test_phantom(): assert_(rmse_den_wrc < rmse_noisy_wrc) -def test_lpca_ill_conditioned(): - DWI, sigma = rfiw_phantom(gtab, snr=30) +@set_random_number_generator() +def test_lpca_ill_conditioned(rng): + DWI, sigma = rfiw_phantom(gtab, snr=30, rng=rng) for patch_radius in [1, [1, 1, 1]]: assert_warns(UserWarning, localpca, DWI, sigma, patch_radius=patch_radius) -def test_lpca_radius_wrong_shape(): - DWI, sigma = rfiw_phantom(gtab, snr=30) +@set_random_number_generator() +def test_lpca_radius_wrong_shape(rng): + DWI, sigma = rfiw_phantom(gtab, snr=30, rng=rng) for patch_radius in [[2, 2], [2, 2, 2, 2]]: assert_raises(ValueError, localpca, DWI, sigma, patch_radius=patch_radius) -def test_lpca_sigma_wrong_shape(): - DWI, sigma = rfiw_phantom(gtab, snr=30) +@set_random_number_generator() +def test_lpca_sigma_wrong_shape(rng): + DWI, sigma = rfiw_phantom(gtab, snr=30, rng=rng) # If sigma is 3D but shape is not like DWI.shape[:-1], an error is raised: sigma = np.zeros((DWI.shape[0], DWI.shape[1] + 1, DWI.shape[2])) assert_raises(ValueError, localpca, DWI, sigma) -def test_lpca_no_gtab_no_sigma(): - DWI, sigma = rfiw_phantom(gtab, snr=30) +@set_random_number_generator() +def test_lpca_no_gtab_no_sigma(rng): + DWI, sigma = rfiw_phantom(gtab, snr=30, rng=rng) assert_raises(ValueError, localpca, DWI, None, None) -def test_pca_classifier(): +@set_random_number_generator() +def test_pca_classifier(rng): # Produce small phantom with well aligned single voxels and ground truth # snr = 50, i.e signal std = 0.02 (Gaussian noise) std_gt = 0.02 @@ -277,7 +290,7 @@ def test_pca_classifier(): angles=[(0, 0, 1), (0, 0, 1)], fractions=(50, 50), snr=None) signal_test[..., :] = sig - noise = std_gt*np.random.standard_normal((5, 5, 5, ndir)) + noise = std_gt*rng.standard_normal((5, 5, 5, ndir)) dwi_test = signal_test + noise # Compute eigenvalues @@ -302,10 +315,11 @@ def test_pca_classifier(): assert_(std_error < 5) -def test_mppca_in_phantom(): - DWIgt = rfiw_phantom(gtab, snr=None) +@set_random_number_generator() +def test_mppca_in_phantom(rng): + DWIgt = rfiw_phantom(gtab, snr=None, rng=rng) std_gt = 0.02 - noise = std_gt*np.random.standard_normal(DWIgt.shape) + noise = std_gt*rng.standard_normal(DWIgt.shape) DWInoise = DWIgt + noise # patch radius (2: #samples > #features, 1: #samples < #features) @@ -319,10 +333,11 @@ def test_mppca_in_phantom(): assert_(rmse_den < rmse_noisy) -def test_mppca_returned_sigma(): - DWIgt = rfiw_phantom(gtab, snr=None) +@set_random_number_generator() +def test_mppca_returned_sigma(rng): + DWIgt = rfiw_phantom(gtab, snr=None, rng=rng) std_gt = 0.02 - noise = std_gt*np.random.standard_normal(DWIgt.shape) + noise = std_gt*rng.standard_normal(DWIgt.shape) DWInoise = DWIgt + noise # patch radius (2: #samples > #features, 1: #samples < #features) diff --git a/dipy/denoise/tests/test_nlmeans.py b/dipy/denoise/tests/test_nlmeans.py index 90cdabed3c..5fa2845ff6 100644 --- a/dipy/denoise/tests/test_nlmeans.py +++ b/dipy/denoise/tests/test_nlmeans.py @@ -1,17 +1,21 @@ +from time import time + import numpy as np from numpy.testing import (assert_, assert_equal, assert_array_almost_equal, assert_raises) import pytest + from dipy.denoise.nlmeans import nlmeans from dipy.denoise.denspeed import (add_padding_reflection, remove_padding) from dipy.utils.omp import cpu_count, have_openmp -from time import time +from dipy.testing.decorators import set_random_number_generator -def test_nlmeans_padding(): - S0 = 100 + 2 * np.random.standard_normal((50, 50, 50)) +@set_random_number_generator() +def test_nlmeans_padding(rng): + S0 = 100 + 2 * rng.standard_normal((50, 50, 50)) S0 = S0.astype('f8') S0n = add_padding_reflection(S0, 5) S0n2 = remove_padding(S0n, 5) @@ -34,8 +38,9 @@ def test_nlmeans_wrong(): assert_raises(ValueError, nlmeans, data, sigma, num_threads=0) -def test_nlmeans_random_noise(): - S0 = 100 + 2 * np.random.standard_normal((22, 23, 30)) +@set_random_number_generator() +def test_nlmeans_random_noise(rng): + S0 = 100 + 2 * rng.standard_normal((22, 23, 30)) S0n = nlmeans(S0, sigma=np.ones((22, 23, 30)) * np.std(S0), rician=False) @@ -47,12 +52,13 @@ def test_nlmeans_random_noise(): assert_equal(np.round(S0n.mean()), 100) -def test_nlmeans_boundary(): +@set_random_number_generator() +def test_nlmeans_boundary(rng): # nlmeans preserves boundaries S0 = 100 + np.zeros((20, 20, 20)) - noise = 2 * np.random.standard_normal((20, 20, 20)) + noise = 2 * rng.standard_normal((20, 20, 20)) S0 += noise diff --git a/dipy/denoise/tests/test_noise_estimate.py b/dipy/denoise/tests/test_noise_estimate.py index 6481977b47..d428c6941a 100644 --- a/dipy/denoise/tests/test_noise_estimate.py +++ b/dipy/denoise/tests/test_noise_estimate.py @@ -26,7 +26,8 @@ def test_inv_nchi(): assert_almost_equal(lambdaPlus, 9.722849086419043) -def test_piesno(): +@set_random_number_generator() +def test_piesno(rng): # Values taken from hispeed.OptimalPIESNO with the test data # in the package computed in matlab test_piesno_data = load_nifti_data(dpd.get_fnames("test_piesno")) @@ -34,8 +35,8 @@ def test_piesno(): return_mask=False) assert_almost_equal(sigma, 0.010749458025559) - noise1 = (np.random.randn(100, 100, 100) * 50) + 10 - noise2 = (np.random.randn(100, 100, 100) * 50) + 10 + noise1 = (rng.standard_normal((100, 100, 100)) * 50) + 10 + noise2 = (rng.standard_normal((100, 100, 100)) * 50) + 10 rician_noise = np.sqrt(noise1**2 + noise2**2) sigma, mask = piesno(rician_noise, N=1, alpha=0.01, l=1, eps=1e-10, return_mask=True) diff --git a/dipy/denoise/tests/test_non_local_means.py b/dipy/denoise/tests/test_non_local_means.py index 115cc8135a..d001e22cd5 100644 --- a/dipy/denoise/tests/test_non_local_means.py +++ b/dipy/denoise/tests/test_non_local_means.py @@ -4,6 +4,7 @@ assert_array_almost_equal, assert_raises) from dipy.denoise.non_local_means import non_local_means +from dipy.testing.decorators import set_random_number_generator def test_nlmeans_static(): @@ -12,8 +13,9 @@ def test_nlmeans_static(): assert_array_almost_equal(S0, S0nb) -def test_nlmeans_random_noise(): - S0 = 100 + 2 * np.random.standard_normal((22, 23, 30)) +@set_random_number_generator() +def test_nlmeans_random_noise(rng): + S0 = 100 + 2 * rng.standard_normal((22, 23, 30)) masker = np.zeros(S0.shape[:3]).astype(bool) masker[8:15, 8:15, 8:15] = 1 @@ -31,9 +33,10 @@ def test_nlmeans_random_noise(): assert_equal(np.round(S0nb[mask].mean()), 100) -def test_scalar_sigma(): +@set_random_number_generator() +def test_scalar_sigma(rng): S0 = 100 + np.zeros((20, 20, 20)) - noise = 2 * np.random.standard_normal((20, 20, 20)) + noise = 2 * rng.standard_normal((20, 20, 20)) S0 += noise S0[:10, :10, :10] = 300 + noise[:10, :10, :10] @@ -41,10 +44,11 @@ def test_scalar_sigma(): ValueError, non_local_means, S0, sigma=noise, rician=False) -def test_nlmeans_boundary(): +@set_random_number_generator() +def test_nlmeans_boundary(rng): # nlmeans preserves boundaries S0 = 100 + np.zeros((20, 20, 20)) - noise = 2 * np.random.standard_normal((20, 20, 20)) + noise = 2 * rng.standard_normal((20, 20, 20)) S0 += noise S0[:10, :10, :10] = 300 + noise[:10, :10, :10] non_local_means(S0, sigma=np.std(noise), rician=False) diff --git a/dipy/denoise/tests/test_patch2self.py b/dipy/denoise/tests/test_patch2self.py index 3c0d8bc243..b9c6d4f27a 100644 --- a/dipy/denoise/tests/test_patch2self.py +++ b/dipy/denoise/tests/test_patch2self.py @@ -18,7 +18,7 @@ @needs_sklearn @set_random_number_generator(1234) -def test_patch2self_random_noise(rng=None): +def test_patch2self_random_noise(rng): S0 = 30 + 2 * rng.standard_normal((20, 20, 20, 50)) bvals = np.repeat(30, 50) @@ -60,7 +60,7 @@ def test_patch2self_random_noise(rng=None): @needs_sklearn @set_random_number_generator(1234) -def test_patch2self_boundary(rng=None): +def test_patch2self_boundary(rng): # patch2self preserves boundaries S0 = 100 + np.zeros((20, 20, 20, 20)) noise = 2 * rng.standard_normal((20, 20, 20, 20)) @@ -74,7 +74,6 @@ def test_patch2self_boundary(rng=None): assert_less(S0[10, 10, 10, 10], 110) -@set_random_number_generator(4321) def rfiw_phantom(gtab, snr=None, rng=None): """rectangle fiber immersed in water""" # define voxel index @@ -140,7 +139,8 @@ def rfiw_phantom(gtab, snr=None, rng=None): @needs_sklearn -def test_phantom(): +@set_random_number_generator(4321) +def test_phantom(rng): # generate a gradient table for phantom data directions8 = generate_bvecs(8) @@ -157,7 +157,7 @@ def test_phantom(): directions8, directions30, directions60)) gtab = gradient_table(bvals, bvecs) - dwi, sigma = rfiw_phantom(gtab, snr=10) + dwi, sigma = rfiw_phantom(gtab, snr=10, rng=rng) dwi_den1 = p2s.patch2self(dwi, model='ridge', bvals=bvals, alpha=1.0) diff --git a/dipy/direction/tests/test_bootstrap_direction_getter.py b/dipy/direction/tests/test_bootstrap_direction_getter.py index 6b6e6af1a1..f4bcd57242 100644 --- a/dipy/direction/tests/test_bootstrap_direction_getter.py +++ b/dipy/direction/tests/test_bootstrap_direction_getter.py @@ -14,6 +14,7 @@ from dipy.reconst.csdeconv import (ConstrainedSphericalDeconvModel, TensorModel) from dipy.sims.voxel import multi_tensor, single_tensor +from dipy.testing.decorators import set_random_number_generator def test_bdg_initial_direction(): @@ -129,7 +130,8 @@ def test_bdg_get_direction(): max_attempts=0)) -def test_bdg_residual(): +@set_random_number_generator() +def test_bdg_residual(rng): """This tests the bootstrapping residual. """ @@ -146,7 +148,7 @@ def test_bdg_residual(): "ignore", message=shm.descoteaux07_legacy_msg, category=PendingDeprecationWarning) B, m, n = shm.real_sh_descoteaux(6, theta, phi) - shm_coeff = np.random.random(B.shape[1]) + shm_coeff = rng.random(B.shape[1]) # sphere_func is sampled of the spherical function for each point of # the sphere diff --git a/dipy/direction/tests/test_peaks.py b/dipy/direction/tests/test_peaks.py index 403f00031c..f12b9b853c 100644 --- a/dipy/direction/tests/test_peaks.py +++ b/dipy/direction/tests/test_peaks.py @@ -22,6 +22,7 @@ from dipy.core.sphere import HemiSphere from dipy.io.gradients import read_bvals_bvecs from dipy.reconst.shm import descoteaux07_legacy_msg, tournier07_legacy_msg +from dipy.testing.decorators import set_random_number_generator def test_peak_directions_nl(): @@ -403,7 +404,8 @@ def test_difference_with_minmax(): assert_almost_equal(values_1, values_4) -def test_degenerate_cases(): +@set_random_number_generator() +def test_degenerate_cases(rng): sphere = default_sphere @@ -440,13 +442,13 @@ def test_degenerate_cases(): assert_equal(values[0], 0.02) odf = np.ones(sphere.vertices.shape[0]) - odf += 0.1 * np.random.rand(odf.shape[0]) + odf += 0.1 * rng.random(odf.shape[0]) directions, values, indices = peak_directions(odf, sphere, .5, 25) assert_(all(values > values[0] * .5)) assert_array_equal(values, odf[indices]) odf = np.ones(sphere.vertices.shape[0]) - odf[1:] = np.finfo(float).eps * np.random.rand(odf.shape[0] - 1) + odf[1:] = np.finfo(float).eps * rng.random(odf.shape[0] - 1) directions, values, indices = peak_directions(odf, sphere, .5, 25) assert_equal(values[0], 1) @@ -670,11 +672,12 @@ def test_peaks_shm_coeff(): assert_array_almost_equal(pam.odf, odf2) -def test_reshape_peaks_for_visualization(): +@set_random_number_generator() +def test_reshape_peaks_for_visualization(rng): - data1 = np.random.randn(10, 5, 3).astype('float32') - data2 = np.random.randn(10, 2, 5, 3).astype('float32') - data3 = np.random.randn(10, 2, 12, 5, 3).astype('float32') + data1 = rng.standard_normal((10, 5, 3)).astype('float32') + data2 = rng.standard_normal((10, 2, 5, 3)).astype('float32') + data3 = rng.standard_normal((10, 2, 12, 5, 3)).astype('float32') data1_reshape = reshape_peaks_for_visualization(data1) data2_reshape = reshape_peaks_for_visualization(data2) diff --git a/dipy/direction/tests/test_pmf.py b/dipy/direction/tests/test_pmf.py index a99de775d1..93bad90a19 100644 --- a/dipy/direction/tests/test_pmf.py +++ b/dipy/direction/tests/test_pmf.py @@ -9,23 +9,25 @@ from dipy.direction.pmf import SimplePmfGen, SHCoeffPmfGen from dipy.reconst.csdeconv import ConstrainedSphericalDeconvModel from dipy.reconst.dti import TensorModel +from dipy.testing.decorators import set_random_number_generator response = (np.array([1.5e3, 0.3e3, 0.3e3]), 1) -def test_pmf_val(): +@set_random_number_generator() +def test_pmf_val(rng): sphere = get_sphere('symmetric724') with warnings.catch_warnings(): warnings.filterwarnings( 'ignore', message=shm.descoteaux07_legacy_msg, category=PendingDeprecationWarning) - pmfgen = SHCoeffPmfGen(np.random.random([2, 2, 2, 28]), sphere, None) + pmfgen = SHCoeffPmfGen(rng.random([2, 2, 2, 28]), sphere, None) point = np.array([1, 1, 1], dtype='float') for idx in [0, 5, 15, -1]: pmf = pmfgen.get_pmf(point) # Create a direction vector close to the vertex idx - xyz = sphere.vertices[idx] + np.random.random([3]) / 100 + xyz = sphere.vertices[idx] + rng.random([3]) / 100 pmf_idx = pmfgen.get_pmf_value(point, xyz) # Test that the pmf sampled for the direction xyz is correct npt.assert_array_almost_equal(pmf[idx], pmf_idx) diff --git a/dipy/io/tests/test_io_peaks.py b/dipy/io/tests/test_io_peaks.py index 49be372446..26f6bcd6bc 100644 --- a/dipy/io/tests/test_io_peaks.py +++ b/dipy/io/tests/test_io_peaks.py @@ -8,15 +8,17 @@ from dipy.direction.peaks import PeaksAndMetrics from dipy.data import default_sphere from dipy.io.peaks import load_peaks, save_peaks, peaks_to_niftis +from dipy.testing.decorators import set_random_number_generator -def test_io_peaks(): +@set_random_number_generator() +def test_io_peaks(rng): with TemporaryDirectory() as tmpdir: fname = 'test.pam5' pam = PeaksAndMetrics() pam.affine = np.eye(4) - pam.peak_dirs = np.random.rand(10, 10, 10, 5, 3) + pam.peak_dirs = rng.random((10, 10, 10, 5, 3)) pam.peak_values = np.zeros((10, 10, 10, 5)) pam.peak_indices = np.zeros((10, 10, 10, 5)) pam.shm_coeff = np.zeros((10, 10, 10, 45)) @@ -99,7 +101,8 @@ def test_io_peaks(): os.path.isfile(fname_gfa) -def test_io_save_peaks_error(): +@set_random_number_generator() +def test_io_save_peaks_error(rng): with TemporaryDirectory() as tmpdir: fname = 'test.pam5' @@ -110,7 +113,7 @@ def test_io_save_peaks_error(): pjoin(tmpdir, fname), pam) pam.affine = np.eye(4) - pam.peak_dirs = np.random.rand(10, 10, 10, 5, 3) + pam.peak_dirs = rng.random((10, 10, 10, 5, 3)) pam.peak_values = np.zeros((10, 10, 10, 5)) pam.peak_indices = np.zeros((10, 10, 10, 5)) pam.shm_coeff = np.zeros((10, 10, 10, 45)) diff --git a/dipy/io/tests/test_utils.py b/dipy/io/tests/test_utils.py index a486aad35d..e4ea9891ba 100644 --- a/dipy/io/tests/test_utils.py +++ b/dipy/io/tests/test_utils.py @@ -2,6 +2,11 @@ import tempfile import pytest +import nibabel as nib +import numpy as np +from numpy.testing import assert_allclose, assert_array_equal, assert_ +import trx.trx_file_memmap as tmm + from dipy.data import fetch_gold_standard_io from dipy.io.streamline import load_tractogram from dipy.io.utils import (create_nifti_header, @@ -9,11 +14,7 @@ get_reference_info, is_reference_info_valid, read_img_arr_or_path) - -import nibabel as nib -import numpy as np -from numpy.testing import assert_allclose, assert_array_equal, assert_ -import trx.trx_file_memmap as tmm +from dipy.testing.decorators import set_random_number_generator filepath_dix = {} files, folder = fetch_gold_standard_io() @@ -200,10 +201,11 @@ def test_all_zeros_affine(): msg='An all zeros affine should not be valid') -def test_read_img_arr_or_path(): - data = np.random.rand(4, 4, 4, 3) +@set_random_number_generator() +def test_read_img_arr_or_path(rng): + data = rng.random((4, 4, 4, 3)) aff = np.eye(4) - aff[:3, :] = np.random.randn(3, 4) + aff[:3, :] = rng.standard_normal((3, 4)) img = nib.Nifti1Image(data, aff) path = tempfile.NamedTemporaryFile().name + '.nii.gz' nib.save(img, path) diff --git a/dipy/nn/tests/test_cnn_1denoiser.py b/dipy/nn/tests/test_cnn_1denoiser.py index 204fdd4187..8bb6c37f67 100644 --- a/dipy/nn/tests/test_cnn_1denoiser.py +++ b/dipy/nn/tests/test_cnn_1denoiser.py @@ -1,18 +1,21 @@ import pytest -import numpy as np from dipy.utils.optpkg import optional_package +from dipy.testing.decorators import set_random_number_generator tf, have_tf, _ = optional_package('tensorflow') -if have_tf: +sklearn, have_sklearn, _ = optional_package('sklearn.model_selection') +if have_tf and have_sklearn: from dipy.nn.cnn_1d_denoising import Cnn1DDenoiser -@pytest.mark.skipif(not have_tf, reason='Requires TensorFlow') -def test_default_Cnn1DDenoiser_sequential(): +@pytest.mark.skipif(not have_tf or not have_sklearn, + reason='Requires TensorFlow and scikit-learn') +@set_random_number_generator() +def test_default_Cnn1DDenoiser_sequential(rng=None): # Create dummy data - normal_img = np.random.rand(10, 10, 10, 30) - nos_img = normal_img + np.random.normal(loc=0.0, scale=0.1, + normal_img = rng.random((10, 10, 10, 30)) + nos_img = normal_img + rng.normal(loc=0.0, scale=0.1, size=normal_img.shape) - x = np.random.rand(10, 10, 10, 30) + x = rng.random((10, 10, 10, 30)) # Test 1D denoiser model = Cnn1DDenoiser(30) model.compile(optimizer='adam', loss='mean_squared_error', @@ -25,13 +28,15 @@ def test_default_Cnn1DDenoiser_sequential(): assert data.shape == x.shape -@pytest.mark.skipif(not have_tf, reason='Requires TensorFlow') -def test_default_Cnn1DDenoiser_flow(pytestconfig): +@pytest.mark.skipif(not have_tf or not have_sklearn, + reason='Requires TensorFlow and scikit-learn') +@set_random_number_generator() +def test_default_Cnn1DDenoiser_flow(pytestconfig, rng): # Create dummy data - normal_img = np.random.rand(10, 10, 10, 30) - nos_img = normal_img + np.random.normal(loc=0.0, scale=0.1, + normal_img = rng.random((10, 10, 10, 30)) + nos_img = normal_img + rng.normal(loc=0.0, scale=0.1, size=normal_img.shape) - x = np.random.rand(10, 10, 10, 30) + x = rng.random((10, 10, 10, 30)) # Test 1D denoiser with flow API model = Cnn1DDenoiser(30) if pytestconfig.getoption('verbose') > 0: diff --git a/dipy/reconst/tests/test_cross_validation.py b/dipy/reconst/tests/test_cross_validation.py index 4b3542b6cb..f946daf1e6 100644 --- a/dipy/reconst/tests/test_cross_validation.py +++ b/dipy/reconst/tests/test_cross_validation.py @@ -20,8 +20,9 @@ fdata, fbval, fbvec = dpd.get_fnames('small_64D') -def test_coeff_of_determination(): - model = np.random.randn(10, 10, 10, 150) +@set_random_number_generator() +def test_coeff_of_determination(rng): + model = rng.standard_normal((10, 10, 10, 150)) data = np.copy(model) # If the model predicts the data perfectly, the COD is all 100s: cod = xval.coeff_of_determination(data, model) diff --git a/dipy/reconst/tests/test_csdeconv.py b/dipy/reconst/tests/test_csdeconv.py index b28c8e83fb..27dc65e5c6 100644 --- a/dipy/reconst/tests/test_csdeconv.py +++ b/dipy/reconst/tests/test_csdeconv.py @@ -36,6 +36,7 @@ from dipy.reconst.shm import lazy_index from dipy.utils.deprecator import ExpiredDeprecationError from dipy.io.gradients import read_bvals_bvecs +from dipy.testing.decorators import set_random_number_generator def get_test_data(): @@ -524,7 +525,8 @@ def test_r2_term_odf_sharp(): assert_equal(directions.shape[0], 2) -def test_csd_predict(): +@set_random_number_generator() +def test_csd_predict(rng): """ Test prediction API """ @@ -572,7 +574,7 @@ def test_csd_predict(): # For "well behaved" coefficients, the model should be able to find the # coefficients from the predicted signal. - coeff = np.random.random(csd_fit.shm_coeff.shape) - .5 + coeff = rng.random(csd_fit.shm_coeff.shape) - .5 coeff[..., 0] = 10. with warnings.catch_warnings(): warnings.filterwarnings( @@ -595,7 +597,8 @@ def test_csd_predict(): npt.assert_array_almost_equal(predict1, predict2) -def test_csd_predict_multi(): +@set_random_number_generator() +def test_csd_predict_multi(rng): """ Check that we can predict reasonably from multi-voxel fits: @@ -610,7 +613,7 @@ def test_csd_predict_multi(): "ignore", message=descoteaux07_legacy_msg, category=PendingDeprecationWarning) csd = ConstrainedSphericalDeconvModel(gtab, response) - coeff = np.random.random(45) - .5 + coeff = rng.random(45) - .5 coeff[..., 0] = 10. with warnings.catch_warnings(): warnings.filterwarnings( diff --git a/dipy/reconst/tests/test_cti.py b/dipy/reconst/tests/test_cti.py index 0df7663e81..6e30500da1 100644 --- a/dipy/reconst/tests/test_cti.py +++ b/dipy/reconst/tests/test_cti.py @@ -69,9 +69,10 @@ def _perpendicular_directions_temp(v, num=20, half=False): # compared to QTE, while Kmicro would be zero. +rng = np.random.default_rng(1234) n_pts = 20 -theta = np.pi * np.random.rand(n_pts) -phi = 2 * np.pi * np.random.rand(n_pts) +theta = np.pi * rng.random(n_pts) +phi = 2 * np.pi * rng.random(n_pts) hsph_initial = HemiSphere(theta=theta, phi=phi) hsph_updated, potential = disperse_charges(hsph_initial, 5000) diff --git a/dipy/reconst/tests/test_dti.py b/dipy/reconst/tests/test_dti.py index c44d37596a..bf82aba9a0 100644 --- a/dipy/reconst/tests/test_dti.py +++ b/dipy/reconst/tests/test_dti.py @@ -24,6 +24,8 @@ from dipy.sims.voxel import single_tensor +from dipy.testing.decorators import set_random_number_generator + def test_roll_evals(): # Just making sure this never passes through @@ -31,9 +33,10 @@ def test_roll_evals(): npt.assert_raises(ValueError, dti._roll_evals, weird_evals) -def test_tensor_algebra(): +@set_random_number_generator() +def test_tensor_algebra(rng): # Test that the computation of tensor determinant and norm is correct - test_arr = np.random.rand(10, 3, 3) + test_arr = rng.random((10, 3, 3)) t_det = dti.determinant(test_arr) t_norm = dti.norm(test_arr) for i, x in enumerate(test_arr): @@ -511,7 +514,8 @@ def test_mask(): dtifit.S0_hat[0, 0, 0]) -def test_nnls_jacobian_func(): +@set_random_number_generator() +def test_nnls_jacobian_func(rng): b0 = 1000. bval, bvecs = read_bvals_bvecs(*get_fnames('55dir_grad')) gtab = grad.gradient_table(bval, bvecs) @@ -526,7 +530,7 @@ def test_nnls_jacobian_func(): # Signals Y = np.exp(np.dot(X, D_orig)) scale = 10 - error = np.random.normal(scale=scale, size=Y.shape) + error = rng.normal(scale=scale, size=Y.shape) Y = Y + error # although sigma and gmm gradients seem correct from inspection, diff --git a/dipy/reconst/tests/test_eudx_dg.py b/dipy/reconst/tests/test_eudx_dg.py index ab3b972df4..afe9e57f0f 100644 --- a/dipy/reconst/tests/test_eudx_dg.py +++ b/dipy/reconst/tests/test_eudx_dg.py @@ -5,10 +5,11 @@ from dipy.direction.peaks import default_sphere, peaks_from_model from dipy.reconst.shm import descoteaux07_legacy_msg +from dipy.testing.decorators import set_random_number_generator -def test_EuDXDirectionGetter(): - +@set_random_number_generator() +def test_EuDXDirectionGetter(rng): class SillyModel: def fit(self, data, mask=None): return SillyFit(self) @@ -20,7 +21,7 @@ def __init__(self, model): def odf(self, sphere): odf = np.zeros(sphere.theta.shape) - r = np.random.randint(0, len(odf)) + r = rng.integers(0, len(odf)) odf[r] = 1 return odf @@ -29,7 +30,7 @@ def get_direction(dg, point, direction): state = dg.get_direction(point, newdir) return state, np.array(newdir) - data = np.random.random((3, 4, 5, 2)) + data = rng.random((3, 4, 5, 2)) with warnings.catch_warnings(): warnings.filterwarnings( "ignore", message=descoteaux07_legacy_msg, @@ -64,7 +65,7 @@ def get_direction(dg, point, direction): npt.assert_array_almost_equal(nd, -expected_dir) # Check that we can get directions at non-integer points - point += np.random.random(3) + point += rng.random(3) state, nd = get_direction(peaks, point, up) npt.assert_equal(state, 0) diff --git a/dipy/reconst/tests/test_mcsd.py b/dipy/reconst/tests/test_mcsd.py index b7d7d340d7..972feaeb74 100644 --- a/dipy/reconst/tests/test_mcsd.py +++ b/dipy/reconst/tests/test_mcsd.py @@ -32,8 +32,7 @@ [4.0E-4, 4.0E-4, 4.0E-4, 40.]]) -@set_random_number_generator(1234) -def get_test_data(rng=None): +def get_test_data(rng): gtab = get_3shell_gtab() evals_list = [np.array([1.7E-3, 0.4E-3, 0.4E-3]), np.array([6.0E-4, 4.0E-4, 4.0E-4]), @@ -255,8 +254,9 @@ def test_multi_shell_fiber_response(): npt.assert_equal(response.response.shape, (3, 7)) -def test_mask_for_response_msmt(): - gtab, data, masks_gt, _ = get_test_data() +@set_random_number_generator() +def test_mask_for_response_msmt(rng): + gtab, data, masks_gt, _ = get_test_data(rng) with warnings.catch_warnings(record=True) as w: wm_mask, gm_mask, csf_mask = mask_for_response_msmt(gtab, data, @@ -282,8 +282,9 @@ def test_mask_for_response_msmt(): npt.assert_array_almost_equal(masks_gt[2], csf_mask) -def test_mask_for_response_msmt_nvoxels(): - gtab, data, _, _ = get_test_data() +@set_random_number_generator() +def test_mask_for_response_msmt_nvoxels(rng): + gtab, data, _, _ = get_test_data(rng) with warnings.catch_warnings(record=True) as w: wm_mask, gm_mask, csf_mask = mask_for_response_msmt(gtab, data, @@ -339,8 +340,9 @@ def test_mask_for_response_msmt_nvoxels(): npt.assert_equal(csf_nvoxels, 0) -def test_response_from_mask_msmt(): - gtab, data, masks_gt, responses_gt = get_test_data() +@set_random_number_generator() +def test_response_from_mask_msmt(rng): + gtab, data, masks_gt, responses_gt = get_test_data(rng) response_wm, response_gm, response_csf \ = response_from_mask_msmt(gtab, data, masks_gt[0], @@ -364,8 +366,9 @@ def test_response_from_mask_msmt(): npt.assert_array_almost_equal(response_csf[0], responses_gt[2], 1) -def test_auto_response_msmt(): - gtab, data, _, _ = get_test_data() +@set_random_number_generator() +def test_auto_response_msmt(rng): + gtab, data, _, _ = get_test_data(rng) with warnings.catch_warnings(record=True) as w: response_auto_wm, response_auto_gm, response_auto_csf = \ diff --git a/dipy/reconst/tests/test_multi_voxel.py b/dipy/reconst/tests/test_multi_voxel.py index 3c8a79f28d..3c81be6df2 100644 --- a/dipy/reconst/tests/test_multi_voxel.py +++ b/dipy/reconst/tests/test_multi_voxel.py @@ -5,6 +5,7 @@ from dipy.reconst.multi_voxel import _squash, multi_voxel_fit, CallableArray from dipy.core.sphere import unit_icosahedron +from dipy.testing.decorators import set_random_number_generator def test_squash(): @@ -99,7 +100,8 @@ def test_CallableArray(): npt.assert_array_equal(callarray(4), expected) -def test_multi_voxel_fit(): +@set_random_number_generator() +def test_multi_voxel_fit(rng): class SillyModel: @@ -123,7 +125,7 @@ def odf(self, sphere): @property def directions(self): - n = np.random.randint(0, 10) + n = rng.integers(0, 10) return np.zeros((n, 3)) def predict(self, S0): diff --git a/dipy/reconst/tests/test_qti.py b/dipy/reconst/tests/test_qti.py index 527f9fb87d..a9957ffb81 100644 --- a/dipy/reconst/tests/test_qti.py +++ b/dipy/reconst/tests/test_qti.py @@ -8,6 +8,7 @@ from dipy.sims.voxel import vec2vec_rotmat from dipy.core.sphere import disperse_charges, HemiSphere from dipy.reconst.dti import fractional_anisotropy +from dipy.testing.decorators import set_random_number_generator from dipy.utils.optpkg import optional_package cp, have_cvxpy, _ = optional_package("cvxpy") @@ -255,14 +256,13 @@ def test_design_matrix(): [0., 0., 0.]]).T) -def _qti_gtab(): +def _qti_gtab(rng): """Return a gradient table with b0, 2 shells, 30 directions, and linear and planar tensor encoding for fitting QTI.""" - np.random.seed(123) n_dir = 30 hsph_initial = HemiSphere( - theta=np.pi * np.random.rand(n_dir), - phi=2 * np.pi * np.random.rand(n_dir)) + theta=np.pi * rng.random(n_dir), + phi=2 * np.pi * rng.random(n_dir)) hsph_updated, _ = disperse_charges(hsph_initial, 100) directions = hsph_updated.vertices bvecs = np.vstack([np.zeros(3)] + [directions for _ in range(4)]) @@ -276,10 +276,11 @@ def _qti_gtab(): return gtab -def test_ls_sdp_fits(): +@set_random_number_generator(123) +def test_ls_sdp_fits(rng): """Test ordinary and weighted least squares and semidefinite programming QTI fits by comparing the estimated parameters to the ground-truth values.""" - gtab = _qti_gtab() + gtab = _qti_gtab(rng) X = qti.design_matrix(gtab.btens) DTDs = [ _anisotropic_DTD(), @@ -307,7 +308,8 @@ def test_ls_sdp_fits(): cvxpy_solver='SCS'), params, decimal=2) -def test_qti_model(): +@set_random_number_generator(123) +def test_qti_model(rng): """Test the QTI model class.""" # Input validation @@ -315,15 +317,16 @@ def test_qti_model(): npt.assert_raises(ValueError, qti.QtiModel, gtab) gtab = gradient_table(np.ones(1), np.array([[1, 0, 0]]), btens='LTE') npt.assert_warns(UserWarning, qti.QtiModel, gtab) - npt.assert_raises(ValueError, qti.QtiModel, _qti_gtab(), 'non-linear') + npt.assert_raises(ValueError, qti.QtiModel, _qti_gtab(rng), 'non-linear') # Design matrix calculation - gtab = _qti_gtab() + gtab = _qti_gtab(rng) qtimodel = qti.QtiModel(gtab) npt.assert_almost_equal(qtimodel.X, qti.design_matrix(gtab.btens)) -def test_qti_fit(): +@set_random_number_generator(4321) +def test_qti_fit(rng): """Test the QTI fit class.""" # Generate a diffusion tensor distribution @@ -374,7 +377,7 @@ def test_qti_fit(): qti.from_6x6_to_21x1(qti.E_bulk)))[0, 0] # Fit QTI - gtab = _qti_gtab() + gtab = _qti_gtab(rng) if have_cvxpy: for fit_method in ['OLS', 'WLS', 'SDPdc']: diff --git a/dipy/reconst/tests/test_shore_odf.py b/dipy/reconst/tests/test_shore_odf.py index 53156ba5f0..8f6c24b0ba 100644 --- a/dipy/reconst/tests/test_shore_odf.py +++ b/dipy/reconst/tests/test_shore_odf.py @@ -2,6 +2,7 @@ import numpy as np import numpy.testing as npt + from dipy.data import default_sphere, get_isbi2013_2shell_gtab, get_3shell_gtab from dipy.reconst.shore import ShoreModel, shore_matrix from dipy.reconst.shm import sh_to_sf, descoteaux07_legacy_msg @@ -11,6 +12,7 @@ from dipy.core.subdivide_octahedron import create_unit_sphere from dipy.core.sphere_stats import angular_similarity from dipy.reconst.tests.test_dsi import sticks_and_ball_dummies +from dipy.testing.decorators import set_random_number_generator def test_shore_odf(): @@ -78,10 +80,11 @@ def test_shore_odf(): npt.assert_equal(gfa(odf) < 0.1, True) -def test_multivox_shore(): +@set_random_number_generator() +def test_multivox_shore(rng): gtab = get_3shell_gtab() - data = np.random.random([20, 30, 1, gtab.gradients.shape[0]]) + data = rng.random([20, 30, 1, gtab.gradients.shape[0]]) radial_order = 4 zeta = 700 asm = ShoreModel(gtab, radial_order=radial_order, diff --git a/dipy/segment/bundles.py b/dipy/segment/bundles.py index 627b65c784..eb5f424f6a 100644 --- a/dipy/segment/bundles.py +++ b/dipy/segment/bundles.py @@ -127,7 +127,8 @@ def cluster_bundle(bundle, clust_thr, rng, nb_pts=20, select_randomly=500000): White matter tract clust_thr : float clustering threshold used in quickbundlesX - rng : RandomState + rng : np.random.Generator + numpy's random generator for generating random values. nb_pts: integer (default 20) Discretizing streamlines to have nb_points number of points select_randomly: integer (default 500000) @@ -165,7 +166,7 @@ def bundle_shape_similarity(bundle1, bundle2, rng, clust_thr=(5, 3, 1.5), White matter tract from one subject (eg: AF_L) bundle2 : Streamlines White matter tract from another subject (eg: AF_L) - rng : RandomState + rng : np.random.Generator clust_thr : array-like, optional list of clustering thresholds used in quickbundlesX threshold : float, optional @@ -236,8 +237,8 @@ def __init__(self, streamlines, greater_than=50, less_than=1000000, Default: 15. nb_pts : int, optional. Number of points per streamline (default 20) - rng : RandomState - If None define RandomState in initialization function. + rng : np.random.Generator + If None define generator in initialization function. Default: None verbose: bool, optional. If True, log information. @@ -272,7 +273,7 @@ def __init__(self, streamlines, greater_than=50, less_than=1000000, self.start_thr = [40, 25, 20] if rng is None: - self.rng = np.random.RandomState() + self.rng = np.random.default_rng() else: self.rng = rng diff --git a/dipy/segment/clustering.py b/dipy/segment/clustering.py index 95a42d7c42..56f0dcaad2 100644 --- a/dipy/segment/clustering.py +++ b/dipy/segment/clustering.py @@ -691,8 +691,8 @@ def qbx_and_merge(streamlines, thresholds, select_randomly : int Randomly select a specific number of streamlines. If None all the streamlines are used. - rng : RandomState - If None then RandomState is initialized internally. + rng : numpy.random.Generator + If None then generator is initialized internally. verbose : bool, optional. If True, log information. Default False. @@ -719,7 +719,7 @@ def qbx_and_merge(streamlines, thresholds, select_randomly = len_s if rng is None: - rng = np.random.RandomState() + rng = np.random.default_rng() indices = rng.choice(len_s, min(select_randomly, len_s), replace=False) sample_streamlines = set_number_of_points(streamlines, nb_pts) diff --git a/dipy/segment/tests/test_bundles.py b/dipy/segment/tests/test_bundles.py index a10501e3e0..aa2c233dc0 100644 --- a/dipy/segment/tests/test_bundles.py +++ b/dipy/segment/tests/test_bundles.py @@ -1,15 +1,17 @@ import sys -import numpy as np import pytest import warnings +import numpy as np from numpy.testing import assert_equal, assert_almost_equal + from dipy.data import get_fnames from dipy.io.streamline import load_tractogram from dipy.segment.bundles import RecoBundles from dipy.tracking.distances import bundles_distances_mam from dipy.tracking.streamline import Streamlines from dipy.segment.clustering import qbx_and_merge +from dipy.testing.decorators import set_random_number_generator is_big_endian = 'big' in sys.byteorder.lower() @@ -105,11 +107,11 @@ def test_rb_disable_slr(): @pytest.mark.skipif(is_big_endian, reason="Little Endian architecture required") -def test_rb_slr_threads(): +@set_random_number_generator(42) +def test_rb_slr_threads(rng): - rng_multi = np.random.RandomState(42) rb_multi = RecoBundles(f, greater_than=0, clust_thr=10, - rng=np.random.RandomState(42)) + rng=rng) rec_trans_multi_threads, _ = rb_multi.recognize(model_bundle=f2, model_clust_thr=5., reduction_thr=10, @@ -117,7 +119,7 @@ def test_rb_slr_threads(): num_threads=None) rb_single = RecoBundles(f, greater_than=0, clust_thr=10, - rng=np.random.RandomState(42)) + rng=np.random.default_rng(42)) rec_trans_single_thread, _ = rb_single.recognize(model_bundle=f2, model_clust_thr=5., reduction_thr=10, diff --git a/dipy/segment/tests/test_clustering.py b/dipy/segment/tests/test_clustering.py index f17aca6601..4ac593a52b 100644 --- a/dipy/segment/tests/test_clustering.py +++ b/dipy/segment/tests/test_clustering.py @@ -1,13 +1,15 @@ -import numpy as np -import itertools import copy +import itertools + +import numpy as np +from numpy.testing import assert_array_equal, assert_raises, assert_equal from dipy.segment.clustering import Cluster, ClusterCentroid from dipy.segment.clustering import ClusterMap, ClusterMapCentroid from dipy.segment.clustering import Clustering from dipy.testing import assert_true, assert_false, assert_arrays_equal -from numpy.testing import assert_array_equal, assert_raises, assert_equal +from dipy.testing.decorators import set_random_number_generator features_shape = (1, 10) @@ -60,9 +62,10 @@ def test_cluster_assign(): assert_array_equal(cluster.indices, indices) -def test_cluster_iter(): +@set_random_number_generator() +def test_cluster_iter(rng): indices = list(range(len(data))) - np.random.shuffle(indices) # None trivial ordering + rng.shuffle(indices) # None trivial ordering # Test without specifying refdata cluster = Cluster() @@ -75,9 +78,10 @@ def test_cluster_iter(): assert_arrays_equal(list(cluster), [data[i] for i in indices]) -def test_cluster_getitem(): +@set_random_number_generator() +def test_cluster_getitem(rng): indices = list(range(len(data))) - np.random.shuffle(indices) # None trivial ordering + rng.shuffle(indices) # None trivial ordering advanced_indices = indices + [0, 1, 2, -1, -2, -3] # Test without specifying refdata in ClusterMap @@ -132,9 +136,10 @@ def test_cluster_getitem(): assert_raises(TypeError, cluster.__getitem__, "wrong") -def test_cluster_str_and_repr(): +@set_random_number_generator() +def test_cluster_str_and_repr(rng): indices = list(range(len(data))) - np.random.shuffle(indices) # None trivial ordering + rng.shuffle(indices) # None trivial ordering # Test without specifying refdata in ClusterMap cluster = Cluster() @@ -187,9 +192,10 @@ def test_cluster_centroid_assign(): assert_array_equal(cluster.centroid, centroid) -def test_cluster_centroid_iter(): +@set_random_number_generator() +def test_cluster_centroid_iter(rng): indices = list(range(len(data))) - np.random.shuffle(indices) # None trivial ordering + rng.shuffle(indices) # None trivial ordering # Test without specifying refdata in ClusterCentroid centroid = np.zeros(features_shape) @@ -205,9 +211,10 @@ def test_cluster_centroid_iter(): assert_arrays_equal(list(cluster), [data[i] for i in indices]) -def test_cluster_centroid_getitem(): +@set_random_number_generator() +def test_cluster_centroid_getitem(rng): indices = list(range(len(data))) - np.random.shuffle(indices) # None trivial ordering + rng.shuffle(indices) # None trivial ordering advanced_indices = indices + [0, 1, 2, -1, -2, -3] # Test without specifying refdata in ClusterCentroid @@ -355,15 +362,15 @@ def test_cluster_map_clear(): assert_array_equal(list(itertools.chain(*clusters)), []) -def test_cluster_map_iter(): - rng = np.random.RandomState(42) +@set_random_number_generator(42) +def test_cluster_map_iter(rng): nb_clusters = 11 # Test without specifying refdata in ClusterMap cluster_map = ClusterMap() clusters = [] for i in range(nb_clusters): - new_cluster = Cluster(indices=rng.randint(0, len(data), size=10)) + new_cluster = Cluster(indices=rng.integers(0, len(data), size=10)) cluster_map.add_cluster(new_cluster) clusters.append(new_cluster) @@ -383,10 +390,11 @@ def test_cluster_map_iter(): assert_array_equal(cluster_map, [cluster.indices for cluster in clusters]) -def test_cluster_map_getitem(): +@set_random_number_generator() +def test_cluster_map_getitem(rng): nb_clusters = 11 indices = list(range(nb_clusters)) - np.random.shuffle(indices) # None trivial ordering + rng.shuffle(indices) # None trivial ordering advanced_indices = indices + [0, 1, 2, -1, -2, -3] cluster_map = ClusterMap() @@ -441,11 +449,11 @@ def test_cluster_map_size(): assert_equal(cluster_map.size(), nb_clusters) -def test_cluster_map_clusters_sizes(): - rng = np.random.RandomState(42) +@set_random_number_generator(42) +def test_cluster_map_clusters_sizes(rng): nb_clusters = 11 # Generate random indices - indices = [range(rng.randint(1, 10)) for _ in range(nb_clusters)] + indices = [range(rng.integers(1, 10)) for _ in range(nb_clusters)] cluster_map = ClusterMap() clusters = [Cluster(indices=indices[i]) for i in range(nb_clusters)] @@ -454,18 +462,18 @@ def test_cluster_map_clusters_sizes(): assert_equal(cluster_map.clusters_sizes(), list(map(len, indices))) -def test_cluster_map_get_small_and_large_clusters(): - rng = np.random.RandomState(42) +@set_random_number_generator(42) +def test_cluster_map_get_small_and_large_clusters(rng): nb_clusters = 11 cluster_map = ClusterMap() # Randomly generate small clusters - indices = [rng.randint(0, 10, size=i) for i in range(1, nb_clusters+1)] + indices = [rng.integers(0, 10, size=i) for i in range(1, nb_clusters+1)] small_clusters = [Cluster(indices=indices[i]) for i in range(nb_clusters)] cluster_map.add_cluster(*small_clusters) # Randomly generate small clusters - indices = [rng.randint(0, 10, size=i) + indices = [rng.integers(0, 10, size=i) for i in range(nb_clusters+1, 2*nb_clusters+1)] large_clusters = [Cluster(indices=indices[i]) for i in range(nb_clusters)] cluster_map.add_cluster(*large_clusters) @@ -591,18 +599,19 @@ def test_cluster_map_centroid_add_cluster(): assert_raises(ValueError, cluster.assign, 123, features_too_long) -def test_cluster_map_centroid_remove_cluster(): +@set_random_number_generator() +def test_cluster_map_centroid_remove_cluster(rng): clusters = ClusterMapCentroid() - centroid1 = np.random.rand(*features_shape).astype(dtype) + centroid1 = rng.random(features_shape).astype(dtype) cluster1 = ClusterCentroid(centroid1, indices=[1]) clusters.add_cluster(cluster1) - centroid2 = np.random.rand(*features_shape).astype(dtype) + centroid2 = rng.random(features_shape).astype(dtype) cluster2 = ClusterCentroid(centroid2, indices=[1, 2]) clusters.add_cluster(cluster2) - centroid3 = np.random.rand(*features_shape).astype(dtype) + centroid3 = rng.random(features_shape).astype(dtype) cluster3 = ClusterCentroid(centroid3, indices=[1, 2, 3]) clusters.add_cluster(cluster3) @@ -628,8 +637,8 @@ def test_cluster_map_centroid_remove_cluster(): assert_array_equal(clusters.centroids, []) -def test_cluster_map_centroid_iter(): - rng = np.random.RandomState(42) +@set_random_number_generator(42) +def test_cluster_map_centroid_iter(rng): nb_clusters = 11 cluster_map = ClusterMapCentroid() @@ -637,8 +646,8 @@ def test_cluster_map_centroid_iter(): for i in range(nb_clusters): new_centroid = np.zeros_like(features) new_cluster = ClusterCentroid(new_centroid, - indices=rng.randint(0, len(data), - size=10)) + indices=rng.integers(0, len(data), + size=10)) cluster_map.add_cluster(new_cluster) clusters.append(new_cluster) @@ -654,10 +663,11 @@ def test_cluster_map_centroid_iter(): assert_arrays_equal(c1, [data[i] for i in c2.indices]) -def test_cluster_map_centroid_getitem(): +@set_random_number_generator() +def test_cluster_map_centroid_getitem(rng): nb_clusters = 11 indices = list(range(len(data))) - np.random.shuffle(indices) # None trivial ordering + rng.shuffle(indices) # None trivial ordering advanced_indices = indices + [0, 1, 2, -1, -2, -3] cluster_map = ClusterMapCentroid() diff --git a/dipy/segment/tests/test_feature.py b/dipy/segment/tests/test_feature.py index 729b708aaf..931743501e 100644 --- a/dipy/segment/tests/test_feature.py +++ b/dipy/segment/tests/test_feature.py @@ -1,20 +1,23 @@ import sys + import numpy as np +from numpy.testing import (assert_array_equal, assert_array_almost_equal, + assert_raises, assert_equal,) + import dipy.segment.featurespeed as dipysfeature import dipy.segment.metric as dipymetric import dipy.segment.metricspeed as dipysmetric from dipy.segment.featurespeed import extract - from dipy.testing import assert_true, assert_false -from numpy.testing import (assert_array_equal, assert_array_almost_equal, - assert_raises, assert_equal,) +from dipy.testing.decorators import set_random_number_generator dtype = "float32" +rng = np.random.default_rng() s1 = np.array([np.arange(10, dtype=dtype)]*3).T # 10x3 s2 = np.arange(3*10, dtype=dtype).reshape((-1, 3))[::-1] # 10x3 -s3 = np.random.rand(5, 4).astype(dtype) # 5x4 -s4 = np.random.rand(5, 3).astype(dtype) # 5x3 +s3 = rng.random((5, 4), dtype=dtype) # 5x4 +s4 = rng.random((5, 3), dtype=dtype) # 5x3 def test_identity_feature(): @@ -222,7 +225,8 @@ def extract(self, streamline): assert_array_almost_equal(features, -features_flip) -def test_feature_extract(): +@set_random_number_generator(1234) +def test_feature_extract(rng): # Test that features are automatically cast into float32 when # coming from Python space class CenterOfMass64bit(dipysfeature.Feature): @@ -232,12 +236,11 @@ def infer_shape(self, streamline): def extract(self, streamline): return np.mean(streamline.astype(np.float64), axis=0) - rng = np.random.RandomState(1234) nb_streamlines = 100 feature_shape = (1, 3) # One N-dimensional point feature = CenterOfMass64bit() - nb_points = rng.randint(20, 30, size=(nb_streamlines,)) * 3 + nb_points = rng.integers(20, 30, size=(nb_streamlines,)) * 3 streamlines = [np.arange(nb).reshape((-1, 3)).astype(np.float32) for nb in nb_points] features = extract(feature, streamlines) diff --git a/dipy/segment/tests/test_metric.py b/dipy/segment/tests/test_metric.py index d9e1492e21..63c39ccfa2 100644 --- a/dipy/segment/tests/test_metric.py +++ b/dipy/segment/tests/test_metric.py @@ -1,13 +1,15 @@ +import itertools + import numpy as np +from numpy.testing import (assert_array_equal, assert_raises, + assert_almost_equal, assert_equal) + import dipy.segment.featurespeed as dipysfeature import dipy.segment.metric as dipymetric import dipy.segment.metricspeed as dipysmetric -import itertools - from dipy.testing import (assert_true, assert_false, assert_greater_equal, assert_less_equal) -from numpy.testing import (assert_array_equal, assert_raises, - assert_almost_equal, assert_equal) +from dipy.testing.decorators import set_random_number_generator def norm(x, order=None, axis=None): @@ -22,9 +24,9 @@ def norm(x, order=None, axis=None): # Create wiggling streamline nb_points = 18 -rng = np.random.RandomState(42) +rng = np.random.default_rng(42) x = np.linspace(0, 10, nb_points) -y = rng.rand(nb_points) +y = rng.random(nb_points) z = np.sin(np.linspace(0, np.pi, nb_points)) # Bending s = np.array([x, y, z], dtype=dtype).T @@ -225,13 +227,14 @@ class EmptyMetric(dipysmetric.Metric): assert_raises(NotImplementedError, metric.dist, None, None) -def test_distance_matrix(): +@set_random_number_generator() +def test_distance_matrix(rng): metric = dipysmetric.SumPointwiseEuclideanMetric() for dtype in [np.int32, np.int64, np.float32, np.float64]: # Compute distances of all tuples spawn by the Cartesian product # of `data` with itself. - data = (np.random.rand(4, 10, 3)*10).astype(dtype) + data = (rng.random((4, 10, 3))*10).astype(dtype) D = dipysmetric.distance_matrix(metric, data) assert_equal(D.shape, (len(data), len(data))) assert_array_equal(np.diag(D), np.zeros(len(data))) @@ -247,7 +250,7 @@ def test_distance_matrix(): # Compute distances of all tuples spawn by the Cartesian product # of `data` with `data2`. - data2 = (np.random.rand(3, 10, 3)*10).astype(dtype) + data2 = (rng.random((3, 10, 3))*10).astype(dtype) D = dipysmetric.distance_matrix(metric, data, data2) assert_equal(D.shape, (len(data), len(data2))) @@ -257,12 +260,13 @@ def test_distance_matrix(): data2[j])) -def test_mean_distances(): +@set_random_number_generator() +def test_mean_distances(rng): nb_slines = 10 nb_pts = 22 dim = 3 - a = np.random.rand(nb_slines, nb_pts, dim) - b = np.random.rand(nb_slines, nb_pts, dim) + a = rng.random((nb_slines, nb_pts, dim)) + b = rng.random((nb_slines, nb_pts, dim)) diff = a - b # Test Euclidean distance (L2) diff --git a/dipy/segment/tests/test_mrf.py b/dipy/segment/tests/test_mrf.py index 54b9c5daa2..e08a96a254 100644 --- a/dipy/segment/tests/test_mrf.py +++ b/dipy/segment/tests/test_mrf.py @@ -8,34 +8,32 @@ from dipy.testing.decorators import set_random_number_generator -# Load a coronal slice from a T1-weighted MRI -fname = get_fnames('t1_coronal_slice') -single_slice = np.load(fname) +def create_image(): + # Load a coronal slice from a T1-weighted MRI + fname = get_fnames('t1_coronal_slice') + single_slice = np.load(fname) -# Stack a few copies to form a 3D volume -nslices = 5 -image = np.zeros(shape=single_slice.shape + (nslices,)) -image[..., :nslices] = single_slice[..., None] + # Stack a few copies to form a 3D volume + nslices = 5 + image = np.zeros(shape=single_slice.shape + (nslices,)) + image[..., :nslices] = single_slice[..., None] + return image -# Set up parameters -nclasses = 4 -beta = np.float64(0.0) -max_iter = 10 -background_noise = True # Making squares -square = np.zeros((256, 256, 3), dtype=np.int16) -square[42:213, 42:213, :] = 1 -square[71:185, 71:185, :] = 2 -square[99:157, 99:157, :] = 3 +def create_square(): + square = np.zeros((256, 256, 3), dtype=np.int16) + square[42:213, 42:213, :] = 1 + square[71:185, 71:185, :] = 2 + square[99:157, 99:157, :] = 3 + return square -@set_random_number_generator(1234) -def create_square_uniform(rng=None): + +def create_square_uniform(rng): square_1 = np.zeros((256, 256, 3)) + 0.001 square_1 = add_noise(square_1, 10000, 1, noise_type='gaussian', rng=rng) - rng = np.random.default_rng(1234) temp_1 = rng.integers(1, 21, size=(171, 171, 3)) temp_1 = np.where(temp_1 < 20, 1, 3) square_1[42:213, 42:213, :] = temp_1 @@ -49,13 +47,11 @@ def create_square_uniform(rng=None): return square_1 -@set_random_number_generator(1234) -def create_square_gauss(rng=None): +def create_square_gauss(rng): square_gauss = np.zeros((256, 256, 3)) + 0.001 square_gauss = add_noise(square_gauss, 10000, 1, noise_type='gaussian', rng=rng) square_gauss[42:213, 42:213, :] = 1 - rng = np.random.default_rng(1234) noise_1 = rng.normal(1.001, 0.0001, size=square_gauss[42:213, 42:213, :].shape) square_gauss[42:213, 42:213, :] = square_gauss[42:213, 42:213, :] + noise_1 @@ -73,6 +69,10 @@ def create_square_gauss(rng=None): def test_grayscale_image(): + nclasses = 4 + beta = np.float64(0.0) + + image = create_image() com = ConstantObservationModel() icm = IteratedConditionalModes() @@ -139,9 +139,12 @@ def test_grayscale_image(): def test_grayscale_iter(): - max_iter = 15 + nclasses = 4 beta = np.float64(0.1) + max_iter = 15 + background_noise = True + image = create_image() com = ConstantObservationModel() icm = IteratedConditionalModes() @@ -240,13 +243,18 @@ def test_grayscale_iter(): npt.assert_(np.abs(np.sum(difference_map)) != 0) -def test_square_iter(): +@set_random_number_generator() +def test_square_iter(rng): + + nclasses = 4 + beta = np.float64(0.0) + max_iter = 10 com = ConstantObservationModel() icm = IteratedConditionalModes() - initial_segmentation = square - square_gauss = create_square_gauss() + initial_segmentation = create_square() + square_gauss = create_square_gauss(rng) mu, sigmasq = com.seg_stats(square_gauss, initial_segmentation, nclasses) @@ -333,13 +341,17 @@ def test_square_iter(): npt.assert_(np.abs(np.sum(difference_map)) == 0.0) -def test_icm_square(): +@set_random_number_generator() +def test_icm_square(rng): + + nclasses = 4 + max_iter = 10 com = ConstantObservationModel() icm = IteratedConditionalModes() - initial_segmentation = square - square_1 = create_square_uniform() + initial_segmentation = create_square() + square_1 = create_square_uniform(rng) mu, sigma = com.seg_stats(square_1, initial_segmentation, nclasses) @@ -371,7 +383,7 @@ def test_icm_square(): initial_segmentation = final_segmentation_1 beta = 2 - initial_segmentation = square + initial_segmentation = create_square() for j in range(max_iter): @@ -391,10 +403,13 @@ def test_classify(): imgseg = TissueClassifierHMRF() + nclasses = 4 beta = 0.1 tolerance = 0.0001 max_iter = 10 + image = create_image() + npt.assert_(image.max() == 1.0) npt.assert_(image.min() == 0.0) diff --git a/dipy/segment/tests/test_qbx.py b/dipy/segment/tests/test_qbx.py index bd0eb2ed7f..83e5196238 100644 --- a/dipy/segment/tests/test_qbx.py +++ b/dipy/segment/tests/test_qbx.py @@ -1,4 +1,5 @@ import itertools + import numpy as np from numpy.testing import assert_array_equal, assert_equal, assert_raises @@ -9,18 +10,21 @@ ) from dipy.tracking.streamline import set_number_of_points from dipy.tracking.streamline import Streamlines +from dipy.testing.decorators import set_random_number_generator def straight_bundle(nb_streamlines=1, nb_pts=30, step_size=1, - radius=1, rng=np.random.RandomState(42)): + radius=1, rng=None): + if rng is None: + rng = np.random.default_rng(42) bundle = [] bundle_length = step_size * nb_pts Z = -np.linspace(0, bundle_length, nb_pts) for k in range(nb_streamlines): - theta = rng.rand() * (2*np.pi) - r = radius * rng.rand() + theta = rng.random() * (2*np.pi) + r = radius * rng.random() Xk = np.ones(nb_pts) * (r * np.cos(theta)) Yk = np.ones(nb_pts) * (r * np.sin(theta)) @@ -214,18 +218,18 @@ def test_raise_mdf(): assert_raises(ValueError, QuickBundles, thresholds[1], metric=metric) -def test_qbx_and_merge(): +@set_random_number_generator(42) +def test_qbx_and_merge(rng): # Generate synthetic streamlines bundles = bearing_bundles(4, 2) - bundles.append(straight_bundle(1)) + bundles.append(straight_bundle(1, rng=rng)) streamlines = Streamlines(list(itertools.chain(*bundles))) thresholds = [10, 2, 1] - rng = np.random.RandomState(seed=42) qbxm = qbx_and_merge(streamlines, thresholds, rng=rng) qbxm_centroids = qbxm.centroids diff --git a/dipy/segment/tests/test_quickbundles.py b/dipy/segment/tests/test_quickbundles.py index 3451e25b8d..a07e17e484 100644 --- a/dipy/segment/tests/test_quickbundles.py +++ b/dipy/segment/tests/test_quickbundles.py @@ -64,9 +64,9 @@ def test_quickbundles_shape_incompatibility(): list(itertools.chain(*clusters2))) -def test_quickbundles_2D(): +@set_random_number_generator(7) +def test_quickbundles_2D(rng): # Test quickbundles clustering using 2D points and the Eulidean metric. - rng = np.random.default_rng(7) data = [] data += \ [rng.standard_normal((1, 2)) + np.array([0, 0]) for _ in range(1)] diff --git a/dipy/testing/decorators.py b/dipy/testing/decorators.py index 7c355e644a..6fb08b149d 100644 --- a/dipy/testing/decorators.py +++ b/dipy/testing/decorators.py @@ -7,6 +7,7 @@ import re import os import platform +import inspect import numpy as np SKIP_RE = re.compile("(\s*>>>.*?)(\s*)#\s*skip\s+if\s+(.*)$") @@ -82,10 +83,14 @@ def set_random_number_generator(seed_v=1234): """ def _set_random_number_generator(func): - def _set_random_number_generator_wrapper(*args, **kwargs): + def _set_random_number_generator_wrapper(pytestconfig, *args, **kwargs): rng = np.random.default_rng(seed_v) kwargs['rng'] = rng - output = func(*args, **kwargs) + signature = inspect.signature(func) + if pytestconfig and 'pytestconfig' in signature.parameters.keys(): + output = func(pytestconfig, *args, **kwargs) + else: + output = func(*args, **kwargs) return output return _set_random_number_generator_wrapper return _set_random_number_generator diff --git a/dipy/tracking/mesh.py b/dipy/tracking/mesh.py index c289417bd8..4b2ae53231 100644 --- a/dipy/tracking/mesh.py +++ b/dipy/tracking/mesh.py @@ -18,7 +18,7 @@ def random_coordinates_from_surface(nb_triangles, nb_seed, triangles_mask=None, Specifies which triangles should be chosen (or not) triangles_weight : [n] numpy array Specifies the weight/probability of choosing each triangle - random_seed : int + rand_gen : int The seed for the random seed generator (numpy.random.seed). Returns @@ -44,15 +44,15 @@ def random_coordinates_from_surface(nb_triangles, nb_seed, triangles_mask=None, triangles_weight /= triangles_weight.sum() # Set the random seed generator - np.random.seed(rand_gen) + rng = np.random.default_rng(rand_gen) # Choose random triangles triangles_idx = \ - np.random.choice(nb_triangles, size=nb_seed, p=triangles_weight) + rng.choice(nb_triangles, size=nb_seed, p=triangles_weight) # Choose random trilinear coordinates # https://mathworld.wolfram.com/TrianglePointPicking.html - trilin_coord = np.random.rand(nb_seed, 3) + trilin_coord = rng.random((nb_seed, 3)) is_upper = (trilin_coord[:, 1:].sum(axis=-1) > 1.0) trilin_coord[is_upper] = 1.0 - trilin_coord[is_upper] trilin_coord[:, 0] = 1.0 - (trilin_coord[:, 1] + trilin_coord[:, 2]) diff --git a/dipy/tracking/stopping_criterion.pyx b/dipy/tracking/stopping_criterion.pyx index 58f6af591c..f08114ba00 100644 --- a/dipy/tracking/stopping_criterion.pyx +++ b/dipy/tracking/stopping_criterion.pyx @@ -277,18 +277,19 @@ cdef class CmcStoppingCriterion(AnatomicalStoppingCriterion): raise RuntimeError("Unexpected interpolation error " + "(exclude_map - code:%i)" % exclude_err) + rng = np.random.default_rng() # test if the tracking continues if include_result + exclude_result <= 0: return TRACKPOINT num = max(0, (1 - include_result - exclude_result)) den = num + include_result + exclude_result p = (num / den) ** self.correction_factor - if np.random.random() < p: + if rng.random() < p: return TRACKPOINT # test if the tracking stopped in the include tissue map p = (include_result / (include_result + exclude_result)) - if np.random.random() < p: + if rng.random() < p: return ENDPOINT # the tracking stopped in the exclude tissue map diff --git a/dipy/tracking/tests/test_distances.py b/dipy/tracking/tests/test_distances.py index 780aacbce4..ff79bbab6b 100644 --- a/dipy/tracking/tests/test_distances.py +++ b/dipy/tracking/tests/test_distances.py @@ -1,16 +1,20 @@ import warnings + import numpy as np -from dipy.testing import assert_true from numpy.testing import (assert_array_almost_equal, assert_equal, assert_almost_equal, assert_array_equal) + +from dipy.testing import assert_true +from dipy.testing.decorators import set_random_number_generator from dipy.tracking import distances as pf from dipy.tracking.streamline import set_number_of_points from dipy.data import get_fnames from dipy.io.streamline import load_tractogram -def test_LSCv2(verbose=False): +@set_random_number_generator() +def test_LSCv2(verbose=False, rng=None): xyz1 = np.array([[1, 0, 0], [2, 0, 0], [3, 0, 0]], dtype='float32') xyz2 = np.array([[1, 0, 0], [1, 2, 0], [1, 3, 0]], dtype='float32') xyz3 = np.array([[1.1, 0, 0], [1, 2, 0], [1, 3, 0]], dtype='float32') @@ -25,7 +29,7 @@ def test_LSCv2(verbose=False): pf.local_skeleton_clustering_3pts(T, 0.2) for i in range(40): - xyz = np.random.rand(3, 3).astype('f4') + xyz = rng.random((3, 3), dtype='f4') T.append(xyz) from time import time @@ -48,7 +52,7 @@ def test_LSCv2(verbose=False): T2 = [] for i in range(10**4): - xyz = np.random.rand(10, 3).astype('f4') + xyz = rng.random((10, 3), dtype='f4') T2.append(xyz) t1 = time() C5 = pf.local_skeleton_clustering(T2, .5) diff --git a/dipy/tracking/tests/test_metrics.py b/dipy/tracking/tests/test_metrics.py index 957ddd4e40..f4ec2c5c66 100644 --- a/dipy/tracking/tests/test_metrics.py +++ b/dipy/tracking/tests/test_metrics.py @@ -6,6 +6,7 @@ from dipy.tracking import metrics as tm from dipy.tracking import distances as pf from dipy.utils.deprecator import ExpiredDeprecationError +from dipy.testing.decorators import set_random_number_generator def test_downsample_deprecated(): @@ -13,16 +14,17 @@ def test_downsample_deprecated(): npt.assert_raises(ExpiredDeprecationError, tm.downsample, streamline, 12) -def test_splines(): +@set_random_number_generator() +def test_splines(rng): # create a helix t = np.linspace(0, 1.75*2*np.pi, 100) x = np.sin(t) y = np.cos(t) z = t # add noise - x += np.random.normal(scale=0.1, size=x.shape) - y += np.random.normal(scale=0.1, size=y.shape) - z += np.random.normal(scale=0.1, size=z.shape) + x += rng.normal(scale=0.1, size=x.shape) + y += rng.normal(scale=0.1, size=y.shape) + z += rng.normal(scale=0.1, size=z.shape) xyz = np.vstack((x, y, z)).T # get the B-splines smoothed result tm.spline(xyz, 3, 2, -1) diff --git a/dipy/tracking/tests/test_stopping_criterion.py b/dipy/tracking/tests/test_stopping_criterion.py index 115d877799..fad9492ba9 100644 --- a/dipy/tracking/tests/test_stopping_criterion.py +++ b/dipy/tracking/tests/test_stopping_criterion.py @@ -8,14 +8,16 @@ CmcStoppingCriterion, ThresholdStoppingCriterion, StreamlineStatus) +from dipy.testing.decorators import set_random_number_generator -def test_binary_stopping_criterion(): +@set_random_number_generator() +def test_binary_stopping_criterion(rng): """This tests that the binary stopping criterion returns expected streamline statuses. """ - mask = np.random.random((4, 4, 4)) + mask = rng.random((4, 4, 4)) mask[mask < 0.4] = 0.0 btc_boolean = BinaryStoppingCriterion(mask > 0) @@ -36,7 +38,7 @@ def test_binary_stopping_criterion(): # Test random points in voxel for ind in ndindex(mask.shape): for _ in range(50): - pts = np.array(ind, dtype='float64') + np.random.random(3) - 0.5 + pts = np.array(ind, dtype='float64') + rng.random(3) - 0.5 state_boolean = btc_boolean.check_point(pts) state_float64 = btc_float64.check_point(pts) if mask[ind] > 0: @@ -59,12 +61,13 @@ def test_binary_stopping_criterion(): npt.assert_equal(state_float64, int(StreamlineStatus.OUTSIDEIMAGE)) -def test_threshold_stopping_criterion(): +@set_random_number_generator() +def test_threshold_stopping_criterion(rng): """This tests that the threshold stopping criterion returns expected streamline statuses. """ - tissue_map = np.random.random((4, 4, 4)) + tissue_map = rng.random((4, 4, 4)) ttc = ThresholdStoppingCriterion(tissue_map.astype('float32'), 0.5) @@ -100,14 +103,15 @@ def test_threshold_stopping_criterion(): npt.assert_equal(state, int(StreamlineStatus.OUTSIDEIMAGE)) -def test_act_stopping_criterion(): +@set_random_number_generator() +def test_act_stopping_criterion(rng): """This tests that the act stopping criterion returns expected streamline statuses. """ - gm = np.random.random((4, 4, 4)) - wm = np.random.random((4, 4, 4)) - csf = np.random.random((4, 4, 4)) + gm = rng.random((4, 4, 4)) + wm = rng.random((4, 4, 4)) + csf = rng.random((4, 4, 4)) tissue_sum = gm + wm + csf gm /= tissue_sum wm /= tissue_sum diff --git a/dipy/tracking/tests/test_streamline.py b/dipy/tracking/tests/test_streamline.py index 071132a565..2ec7fe9e36 100644 --- a/dipy/tracking/tests/test_streamline.py +++ b/dipy/tracking/tests/test_streamline.py @@ -6,6 +6,7 @@ import numpy.testing as npt from dipy.testing.memory import get_type_refcount from dipy.testing import assert_arrays_equal +from dipy.testing.decorators import set_random_number_generator from dipy.testing import assert_true from numpy.testing import (assert_array_equal, assert_array_almost_equal, @@ -322,14 +323,15 @@ def test_set_number_of_points(): np.ones((10, 3)), np.ones((10, 3))], nb_points=1) -def test_set_number_of_points_memory_leaks(): +@set_random_number_generator(1234) +def test_set_number_of_points_memory_leaks(rng): # Test some dtypes dtypes = [np.float32, np.float64, np.int32, np.int64] for dtype in dtypes: - rng = np.random.default_rng(1234) + s_rng = np.random.default_rng(1234) NB_STREAMLINES = 10000 streamlines = \ - [rng.standard_normal((rng.integers(10, 100), 3)).astype(dtype) + [s_rng.standard_normal((s_rng.integers(10, 100), 3)).astype(dtype) for _ in range(NB_STREAMLINES)] list_refcount_before = get_type_refcount()["list"] @@ -344,12 +346,12 @@ def test_set_number_of_points_memory_leaks(): assert_equal(list_refcount_after, list_refcount_before+1) # Test mixed dtypes - rng = np.random.RandomState(1234) NB_STREAMLINES = 10000 streamlines = [] for i in range(NB_STREAMLINES): dtype = dtypes[i % len(dtypes)] - streamlines.append(rng.randn(rng.randint(10, 100), 3).astype(dtype)) + streamlines.append( + rng.standard_normal((rng.integers(10, 100), 3)).astype(dtype)) list_refcount_before = get_type_refcount()["list"] rstreamlines = set_number_of_points(streamlines, nb_points=2) @@ -467,14 +469,16 @@ def test_length(): [length_python(s) for s in streamlines_readonly]) -def test_length_memory_leaks(): +@set_random_number_generator(1234) +def test_length_memory_leaks(rng): # Test some dtypes dtypes = [np.float32, np.float64, np.int32, np.int64] for dtype in dtypes: - rng = np.random.RandomState(1234) + s_rng = np.random.default_rng(1234) NB_STREAMLINES = 10000 - streamlines = [rng.randn(rng.randint(10, 100), 3).astype(dtype) - for _ in range(NB_STREAMLINES)] + streamlines = \ + [s_rng.standard_normal((s_rng.integers(10, 100), 3)).astype(dtype) + for _ in range(NB_STREAMLINES)] list_refcount_before = get_type_refcount()["list"] @@ -486,12 +490,12 @@ def test_length_memory_leaks(): assert_equal(list_refcount_after, list_refcount_before) # Test mixed dtypes - rng = np.random.RandomState(1234) NB_STREAMLINES = 10000 streamlines = [] for i in range(NB_STREAMLINES): dtype = dtypes[i % len(dtypes)] - streamlines.append(rng.randn(rng.randint(10, 100), 3).astype(dtype)) + streamlines.append( + rng.standard_normal((rng.integers(10, 100), 3)).astype(dtype)) list_refcount_before = get_type_refcount()["list"] @@ -503,10 +507,11 @@ def test_length_memory_leaks(): assert_equal(list_refcount_after, list_refcount_before) -def test_unlist_relist_streamlines(): - streamlines = [np.random.rand(10, 3), - np.random.rand(20, 3), - np.random.rand(5, 3)] +@set_random_number_generator() +def test_unlist_relist_streamlines(rng): + streamlines = [rng.random((10, 3)), + rng.random((20, 3)), + rng.random((5, 3))] points, offsets = unlist_streamlines(streamlines) assert_equal(offsets.dtype, np.dtype('i8')) assert_equal(points.shape, (35, 3)) @@ -550,9 +555,10 @@ def test_transform_empty_streamlines(): assert_equal(len(streamlines), 0) -def test_deform_streamlines(): +@set_random_number_generator() +def test_deform_streamlines(rng): # Create Random deformation field - deformation_field = np.random.randn(200, 200, 200, 3) + deformation_field = rng.standard_normal((200, 200, 200, 3)) stream2grid = np.array([ [-0.13152201, -0.52553149, -0.06759869, -0.80014208], [1.01579851, 0.19840874, 0.18875411, 0.81826065], @@ -610,10 +616,11 @@ def test_center_and_transform(): assert_array_equal(streamlines3[0], B) -def test_select_random_streamlines(): - streamlines = [np.random.rand(10, 3), - np.random.rand(20, 3), - np.random.rand(5, 3)] +@set_random_number_generator() +def test_select_random_streamlines(rng): + streamlines = [rng.random((10, 3)), + rng.random((20, 3)), + rng.random((5, 3))] new_streamlines = select_random_set_of_streamlines(streamlines, 2) assert_equal(len(new_streamlines), 2) @@ -783,14 +790,16 @@ def test_compress_streamlines_identical_points(): [1, 1, 1]])) -def test_compress_streamlines_memory_leaks(): +@set_random_number_generator(1234) +def test_compress_streamlines_memory_leaks(rng): # Test some dtypes dtypes = [np.float32, np.float64, np.int32, np.int64] for dtype in dtypes: - rng = np.random.RandomState(1234) + s_rng = np.random.default_rng(1234) NB_STREAMLINES = 10000 - streamlines = [rng.randn(rng.randint(10, 100), 3).astype(dtype) - for _ in range(NB_STREAMLINES)] + streamlines = \ + [s_rng.standard_normal((s_rng.integers(10, 100), 3)).astype(dtype) + for _ in range(NB_STREAMLINES)] list_refcount_before = get_type_refcount()["list"] @@ -804,12 +813,12 @@ def test_compress_streamlines_memory_leaks(): assert_equal(list_refcount_after, list_refcount_before+1) # Test mixed dtypes - rng = np.random.RandomState(1234) NB_STREAMLINES = 10000 streamlines = [] for i in range(NB_STREAMLINES): dtype = dtypes[i % len(dtypes)] - streamlines.append(rng.randn(rng.randint(10, 100), 3).astype(dtype)) + streamlines.append( + rng.standard_normal((rng.integers(10, 100), 3)).astype(dtype)) list_refcount_before = get_type_refcount()["list"] cstreamlines = compress_streamlines(streamlines) diff --git a/dipy/tracking/tests/test_tracking.py b/dipy/tracking/tests/test_tracking.py index 311bd44884..f5836c41ab 100644 --- a/dipy/tracking/tests/test_tracking.py +++ b/dipy/tracking/tests/test_tracking.py @@ -959,7 +959,8 @@ def test_affine_transformations(): npt.assert_(np.allclose(streamlines_inv[1], expected[1], atol=0.3)) -def test_random_seed_initialization(): +@set_random_number_generator() +def test_random_seed_initialization(rng): """Test that the random generator can be initialized correctly with the tracking seeds. """ @@ -970,13 +971,13 @@ def test_random_seed_initialization(): z = np.array([1., 1, 1, 51.67881720942744]) seeds = np.row_stack([np.column_stack([x, y, z]), - np.random.random((10, 3))]) + rng.random((10, 3))]) sc = BinaryStoppingCriterion(np.ones((4, 4, 4))) dg = ProbabilisticDirectionGetter.from_pmf(pmf, 60, sphere) randoms_seeds = [None, 0, 1, -1, np.iinfo(np.uint32).max + 1] \ - + list(np.random.random(10)) \ - + list(np.random.randint(0, np.iinfo(np.int32).max, 10)) + + list(rng.random(10)) \ + + list(rng.integers(0, np.iinfo(np.int32).max, 10)) for rdm_seed in randoms_seeds: _ = Streamlines(LocalTracking(direction_getter=dg, diff --git a/dipy/tracking/tests/test_utils.py b/dipy/tracking/tests/test_utils.py index 67dd0b73a8..8aae942e11 100644 --- a/dipy/tracking/tests/test_utils.py +++ b/dipy/tracking/tests/test_utils.py @@ -17,6 +17,7 @@ from dipy.tracking.vox2track import streamline_mapping import numpy.testing as npt from dipy.testing import assert_true +from dipy.testing.decorators import set_random_number_generator def make_streamlines(return_seeds=False): @@ -464,7 +465,8 @@ def test_streamline_mapping(): npt.assert_equal(mapping, expected) -def test_length(): +@set_random_number_generator() +def test_length(rng): # Generate a simulated bundle of fibers: n_streamlines = 50 n_pts = 100 @@ -475,8 +477,8 @@ def test_length(): pts = np.vstack((np.cos(2 * t / np.pi), np.zeros(t.shape) + i, t)).T bundle.append(pts) - start = np.random.randint(10, 30, n_streamlines) - end = np.random.randint(60, 100, n_streamlines) + start = rng.integers(10, 30, n_streamlines) + end = rng.integers(60, 100, n_streamlines) bundle = [10 * streamline[start[i]:end[i]] for (i, streamline) in enumerate(bundle)] @@ -486,8 +488,9 @@ def test_length(): npt.assert_equal(this_length, metrics.length(bundle[idx])) -def test_seeds_from_mask(): - mask = np.random.randint(0, 1, size=(10, 10, 10)) +@set_random_number_generator() +def test_seeds_from_mask(rng): + mask = rng.integers(0, 1, size=(10, 10, 10)) seeds = seeds_from_mask(mask, np.eye(4), density=1) npt.assert_equal(mask.sum(), len(seeds)) npt.assert_array_equal(np.argwhere(mask), seeds) @@ -508,8 +511,9 @@ def test_seeds_from_mask(): npt.assert_equal(in_444.sum(), 3 * 4 * 5) -def test_random_seeds_from_mask(): - mask = np.random.randint(0, 1, size=(4, 6, 3)) +@set_random_number_generator() +def test_random_seeds_from_mask(rng): + mask = rng.integers(0, 1, size=(4, 6, 3)) seeds = random_seeds_from_mask(mask, np.eye(4), seeds_count=24, seed_count_per_voxel=True) @@ -624,7 +628,8 @@ def test_reduce_rois(): npt.assert_warns(UserWarning, reduce_rois, [roi2], [True]) -def test_path_length(): +@set_random_number_generator() +def test_path_length(rng): aoi = np.zeros((20, 20, 20), dtype=bool) aoi[0, 0, 0] = 1 @@ -660,7 +665,7 @@ def test_path_length(): aoi[0, 0, 0] = 1 streamlines = [] for i in range(1000): - rando = np.random.random(size=(100, 3)) * 19 + .5 + rando = rng.random(size=(100, 3)) * 19 + .5 assert (rando > .5).all() assert (rando < 19.5).all() streamlines.append(rando) diff --git a/dipy/tracking/utils.py b/dipy/tracking/utils.py index 7c7a69b132..8cbc08aaee 100644 --- a/dipy/tracking/utils.py +++ b/dipy/tracking/utils.py @@ -469,7 +469,7 @@ def random_seeds_from_mask(mask, affine, seeds_count=1, If True, seeds_count is per voxel, else seeds_count is the total number of seeds. random_seed : int - The seed for the random seed generator (numpy.random.seed). + The seed for the random seed generator (numpy.random.Generator). See Also -------- @@ -486,22 +486,22 @@ def random_seeds_from_mask(mask, affine, seeds_count=1, >>> mask[0,0,0] = 1 >>> random_seeds_from_mask(mask, np.eye(4), seeds_count=1, ... seed_count_per_voxel=True, random_seed=1) - array([[-0.0640051 , -0.47407377, 0.04966248]]) + array([[-0.23838787, -0.20150886, 0.31422574]]) >>> random_seeds_from_mask(mask, np.eye(4), seeds_count=6, ... seed_count_per_voxel=True, random_seed=1) - array([[-0.0640051 , -0.47407377, 0.04966248], - [ 0.0507979 , 0.20814782, -0.20909526], - [ 0.46702984, 0.04723225, 0.47268436], - [-0.27800683, 0.37073231, -0.29328084], - [ 0.39286015, -0.16802019, 0.32122912], - [-0.42369171, 0.27991879, -0.06159077]]) + array([[-0.23838787, -0.20150886, 0.31422574], + [-0.41435083, -0.26318949, 0.30127447], + [ 0.44305611, 0.01132755, 0.47624371], + [ 0.30500292, 0.30794079, 0.01532556], + [ 0.03816435, -0.15672913, -0.13093276], + [ 0.12509547, 0.3972138 , 0.27568569]]) >>> mask[0,1,2] = 1 >>> random_seeds_from_mask(mask, np.eye(4), ... seeds_count=2, seed_count_per_voxel=True, random_seed=1) - array([[-0.0640051 , -0.47407377, 0.04966248], - [-0.27800683, 1.37073231, 1.70671916], - [ 0.0507979 , 0.20814782, -0.20909526], - [-0.48962585, 1.00187459, 1.99577329]]) + array([[ 0.30500292, 1.30794079, 2.01532556], + [-0.23838787, -0.20150886, 0.31422574], + [ 0.3702492 , 0.78681721, 2.10314815], + [-0.41435083, -0.26318949, 0.30127447]]) """ mask = np.array(mask, dtype=bool, copy=False, ndmin=3) @@ -509,11 +509,11 @@ def random_seeds_from_mask(mask, affine, seeds_count=1, raise ValueError('mask cannot be more than 3d') # Randomize the voxels - np.random.seed(random_seed) + rng = np.random.default_rng(random_seed) shape = mask.shape mask = mask.flatten() indices = np.arange(len(mask)) - np.random.shuffle(indices) + rng.shuffle(indices) where = [np.unravel_index(i, shape) for i in indices if mask[i] == 1] num_voxels = len(where) @@ -532,9 +532,9 @@ def random_seeds_from_mask(mask, affine, seeds_count=1, if random_seed is not None: s_random_seed = hash((np.sum(s) + 1) * i + random_seed) \ % (2**32 - 1) - np.random.seed(s_random_seed) + rng = np.random.default_rng(s_random_seed) # Generate random triplet - grid = np.random.random(3) + grid = rng.random(3) seed = s + grid - .5 seeds.append(seed) seeds = np.asarray(seeds) diff --git a/dipy/utils/tests/test_arrfuncs.py b/dipy/utils/tests/test_arrfuncs.py index 0029321fbc..25759cf02d 100644 --- a/dipy/utils/tests/test_arrfuncs.py +++ b/dipy/utils/tests/test_arrfuncs.py @@ -9,6 +9,7 @@ from numpy.testing import (assert_array_almost_equal, assert_equal, assert_array_equal) from dipy.testing import assert_true, assert_false +from dipy.testing.decorators import set_random_number_generator NATIVE_ORDER = '<' if sys.byteorder == 'little' else '>' SWAPPED_ORDER = '>' if sys.byteorder == 'little' else '<' @@ -28,8 +29,9 @@ def test_as_native(): assert_equal(narr.dtype.byteorder, NATIVE_ORDER) -def test_pinv(): - arr = np.random.randn(4, 4, 4, 3, 7) +@set_random_number_generator() +def test_pinv(rng): + arr = rng.standard_normal((4, 4, 4, 3, 7)) _pinv = pinv(arr) for i in range(4): for j in range(4): diff --git a/dipy/viz/tests/test_apps.py b/dipy/viz/tests/test_apps.py index 8902e5eaa3..d8173e45e6 100644 --- a/dipy/viz/tests/test_apps.py +++ b/dipy/viz/tests/test_apps.py @@ -23,14 +23,15 @@ @pytest.mark.skipif(skip_it or not has_fury, reason="Needs xvfb") -def test_horizon_events(): +@set_random_number_generator() +def test_horizon_events(rng): # using here MNI template affine 2009a affine = np.array([[1., 0., 0., -98.], [0., 1., 0., -134.], [0., 0., 1., -72.], [0., 0., 0., 1.]]) - data = 255 * np.random.rand(197, 233, 189) + data = 255 * rng.random((197, 233, 189)) vox_size = (1., 1., 1.) images = [(data, affine)] @@ -60,8 +61,8 @@ def test_horizon_events(): @pytest.mark.skipif(skip_it or not has_fury, reason="Needs xvfb") -def test_horizon(): - +@set_random_number_generator() +def test_horizon(rng): s1 = 10 * np.array([[0, 0, 0], [1, 0, 0], [2, 0, 0], @@ -90,7 +91,7 @@ def test_horizon(): [0., 0., 1., -72.], [0., 0., 0., 1.]]) - data = 255 * np.random.rand(197, 233, 189) + data = 255 * rng.random((197, 233, 189)) vox_size = (1., 1., 1.) streamlines._data += np.array([-98., -134., -72.]) diff --git a/dipy/viz/tests/test_fury.py b/dipy/viz/tests/test_fury.py index a5503f0148..d737c63d57 100644 --- a/dipy/viz/tests/test_fury.py +++ b/dipy/viz/tests/test_fury.py @@ -3,7 +3,7 @@ import pytest import numpy as np import numpy.testing as npt -from dipy.testing.decorators import use_xvfb +from dipy.testing.decorators import use_xvfb, set_random_number_generator from dipy.utils.optpkg import optional_package from dipy.data import get_sphere from dipy.align.reslice import reslice @@ -30,9 +30,10 @@ @pytest.mark.skipif(skip_it or not has_fury, reason="Needs xvfb") -def test_slicer(): +@set_random_number_generator() +def test_slicer(rng): scene = window.Scene() - data = (255 * np.random.rand(50, 50, 50)) + data = (255 * rng.random((50, 50, 50))) affine = np.diag([1, 3, 2, 1]) data2, affine2 = reslice(data, affine, zooms=(1, 3, 2), new_zooms=(1, 1, 1)) @@ -113,7 +114,8 @@ def test_contour_from_roi(): @pytest.mark.skipif(skip_it or not has_fury, reason="Needs xvfb") -def test_bundle_maps(): +@set_random_number_generator() +def test_bundle_maps(rng): scene = window.Scene() bundle = fornix_streamlines() bundle, _ = center_streamlines(bundle) @@ -149,7 +151,7 @@ def test_bundle_maps(): scene.clear() nb_points = np.sum([len(b) for b in bundle]) - values = 100 * np.random.rand(nb_points) + values = 100 * rng.random((nb_points)) # values[:nb_points/2] = 0 line = actor.streamtube(bundle, values, linewidth=0.1, lookup_colormap=lut) @@ -161,7 +163,7 @@ def test_bundle_maps(): scene.clear() - colors = np.random.rand(nb_points, 3) + colors = rng.random((nb_points, 3)) # values[:nb_points/2] = 0 line = actor.line(bundle, colors, linewidth=2) @@ -181,7 +183,7 @@ def test_bundle_maps(): actor.line(bundle, (1., 0.5, 0)) actor.line(bundle, np.arange(len(bundle))) actor.line(bundle) - colors = [np.random.rand(*b.shape) for b in bundle] + colors = [rng.random(b.shape) for b in bundle] actor.line(bundle, colors=colors) diff --git a/dipy/viz/tests/test_regtools.py b/dipy/viz/tests/test_regtools.py index 26c0f4d790..a4da740885 100644 --- a/dipy/viz/tests/test_regtools.py +++ b/dipy/viz/tests/test_regtools.py @@ -4,6 +4,7 @@ import pytest from dipy.align.metrics import SSDMetric from dipy.align.imwarp import SymmetricDiffeomorphicRegistration +from dipy.testing.decorators import set_random_number_generator # Conditional import machinery for matplotlib from dipy.utils.optpkg import optional_package @@ -12,12 +13,13 @@ @pytest.mark.skipif(not have_matplotlib, reason='Requires Matplotlib') -def test_plot_2d_diffeomorphic_map(): +@set_random_number_generator() +def test_plot_2d_diffeomorphic_map(rng): # Test the regtools plotting interface (lightly). mv_shape = (11, 12) - moving = np.random.rand(*mv_shape) + moving = rng.random(mv_shape) st_shape = (13, 14) - static = np.random.rand(*st_shape) + static = rng.random(st_shape) dim = static.ndim metric = SSDMetric(dim) level_iters = [200, 100, 50, 25] diff --git a/dipy/workflows/stats.py b/dipy/workflows/stats.py index e524bf9c2c..e4eb7e6114 100755 --- a/dipy/workflows/stats.py +++ b/dipy/workflows/stats.py @@ -553,7 +553,7 @@ def run(self, subject_folder, clust_thr=(5, 3, 1.5), threshold=6, populations. Sci Rep 10, 17149 (2020) """ - rng = np.random.RandomState() + rng = np.random.default_rng() all_subjects = [] if os.path.isdir(subject_folder): groups = os.listdir(subject_folder) diff --git a/dipy/workflows/tests/test_align.py b/dipy/workflows/tests/test_align.py index 05feb10b82..58a129f623 100644 --- a/dipy/workflows/tests/test_align.py +++ b/dipy/workflows/tests/test_align.py @@ -19,6 +19,7 @@ ApplyTransformFlow, ResliceFlow, SlrWithQbxFlow, MotionCorrectionFlow, BundleWarpFlow) +from dipy.testing.decorators import set_random_number_generator _, have_pd, _ = optional_package("pandas") @@ -70,12 +71,13 @@ def test_slr_flow(): npt.assert_equal(os.path.isfile(out_path), True) -def test_image_registration(): +@set_random_number_generator(1234) +def test_image_registration(rng): with TemporaryDirectory() as temp_out_dir: static, moving, static_g2w, moving_g2w, smask, mmask, M\ = setup_random_transform(transform=regtransforms[('AFFINE', 3)], - rfactor=0.1) + rfactor=0.1, rng=rng) save_nifti(pjoin(temp_out_dir, 'b0.nii.gz'), data=static, affine=static_g2w) @@ -122,7 +124,7 @@ def test_translation(): out_quality='trans_q.txt') dist = read_distance('trans_q.txt') - npt.assert_almost_equal(float(dist), -0.3953547764454917, 1) + npt.assert_almost_equal(float(dist), -0.42097809101318934, 1) check_existence(out_moved, out_affine) def test_rigid(): diff --git a/dipy/workflows/tests/test_denoise.py b/dipy/workflows/tests/test_denoise.py index 87bfb81ba7..3a48b4c6d7 100644 --- a/dipy/workflows/tests/test_denoise.py +++ b/dipy/workflows/tests/test_denoise.py @@ -11,6 +11,7 @@ from dipy.testing import (assert_true, assert_false, assert_greater, assert_less) +from dipy.testing.decorators import set_random_number_generator from dipy.workflows.denoise import (NLMeansFlow, LPCAFlow, MPPCAFlow, GibbsRingingFlow, Patch2SelfFlow) @@ -60,9 +61,10 @@ def test_lpca_flow(): lpca_flow.last_generated_outputs['out_denoised'])) -def test_mppca_flow(): +@set_random_number_generator() +def test_mppca_flow(rng): with TemporaryDirectory() as out_dir: - S0 = 100 + 2 * np.random.standard_normal((22, 23, 30, 20)) + S0 = 100 + 2 * rng.standard_normal((22, 23, 30, 20)) data_path = os.path.join(out_dir, "random_noise.nii.gz") save_nifti(data_path, S0, np.eye(4)) diff --git a/dipy/workflows/tests/test_stats.py b/dipy/workflows/tests/test_stats.py index 1528f0f73d..b1193d90ed 100755 --- a/dipy/workflows/tests/test_stats.py +++ b/dipy/workflows/tests/test_stats.py @@ -4,15 +4,16 @@ from os.path import join from tempfile import TemporaryDirectory +import numpy as np import numpy.testing as npt import pytest -import numpy as np from dipy.data import get_fnames from dipy.io.image import save_nifti, load_nifti from dipy.io.stateful_tractogram import Space, StatefulTractogram from dipy.io.streamline import load_tractogram, save_tractogram from dipy.testing import assert_true +from dipy.testing.decorators import set_random_number_generator from dipy.tracking.streamline import Streamlines from dipy.utils.optpkg import optional_package from dipy.workflows.stats import SNRinCCFlow @@ -76,7 +77,8 @@ def test_stats(): @pytest.mark.skipif(not have_pandas or not have_statsmodels or not have_tables or not have_matplotlib, reason='Requires Pandas, StatsModels and PyTables') -def test_buan_bundle_profiles(): +@set_random_number_generator() +def test_buan_bundle_profiles(rng): with TemporaryDirectory() as dirpath: data_path = get_fnames('fornix') @@ -110,7 +112,7 @@ def test_buan_bundle_profiles(): dt = os.path.join(dirpath, "anatomical_measures") os.mkdir(dt) - fa = np.random.rand(255, 255, 255) + fa = rng.random((255, 255, 255)) save_nifti(os.path.join(dt, "fa.nii.gz"), fa, affine=np.eye(4)) @@ -128,7 +130,8 @@ def test_buan_bundle_profiles(): or not have_matplotlib, reason='Requires Pandas, StatsModels, PyTables, and ' 'matplotlib') -def test_bundle_analysis_tractometry_flow(): +@set_random_number_generator() +def test_bundle_analysis_tractometry_flow(rng): with TemporaryDirectory() as dirpath: data_path = get_fnames('fornix') @@ -171,7 +174,7 @@ def test_bundle_analysis_tractometry_flow(): bbox_valid_check=False) os.mkdir(os.path.join(pre, "anatomical_measures")) - fa = np.random.rand(255, 255, 255) + fa = rng.random((255, 255, 255)) save_nifti(os.path.join(pre, "anatomical_measures", "fa.nii.gz"), fa, affine=np.eye(4)) @@ -264,7 +267,8 @@ def test_linear_mixed_models_flow(): or not have_matplotlib, reason='Requires Pandas, StatsModels, PyTables, and ' 'matplotlib') -def test_bundle_shape_analysis_flow(): +@set_random_number_generator() +def test_bundle_shape_analysis_flow(rng): with TemporaryDirectory() as dirpath: data_path = get_fnames('fornix') @@ -307,7 +311,7 @@ def test_bundle_shape_analysis_flow(): bbox_valid_check=False) os.mkdir(os.path.join(pre, "anatomical_measures")) - fa = np.random.rand(255, 255, 255) + fa = rng.random((255, 255, 255)) save_nifti(os.path.join(pre, "anatomical_measures", "fa.nii.gz"), fa, affine=np.eye(4)) diff --git a/dipy/workflows/tests/test_viz.py b/dipy/workflows/tests/test_viz.py index 6742958530..d291d57885 100644 --- a/dipy/workflows/tests/test_viz.py +++ b/dipy/workflows/tests/test_viz.py @@ -12,6 +12,7 @@ from dipy.tracking.streamline import Streamlines from dipy.testing.decorators import use_xvfb from dipy.utils.optpkg import optional_package +from dipy.testing.decorators import set_random_number_generator fury, has_fury, setup_module = optional_package('fury') @@ -25,7 +26,8 @@ @pytest.mark.skipif(skip_it or not has_fury, reason='Requires FURY') -def test_horizon_flow(): +@set_random_number_generator() +def test_horizon_flow(rng): s1 = 10 * np.array([[0, 0, 0], [1, 0, 0], @@ -50,7 +52,7 @@ def test_horizon_flow(): [0., 0., 1., -72.], [0., 0., 0., 1.]]) - data = 255 * np.random.rand(197, 233, 189) + data = 255 * rng.random((197, 233, 189)) vox_size = (1., 1., 1.) streamlines = Streamlines() @@ -75,7 +77,7 @@ def test_horizon_flow(): world_coords=True, interactive=False) - data = 255 * np.random.rand(197, 233, 189) + data = 255 * rng.random((197, 233, 189)) images = [(data, affine)] @@ -96,7 +98,7 @@ def test_horizon_flow(): sft = StatefulTractogram(streamlines, nii_header, space=Space.RASMM) save_tractogram(sft, ftrk, bbox_valid_check=False) - pvalues = np.random.uniform(low=0, high=1, size=(10,)) + pvalues = rng.uniform(low=0, high=1, size=(10,)) np.save(fnpy, pvalues) input_files = [ftrk, fimg] diff --git a/doc/examples/bundle_assignment_maps.py b/doc/examples/bundle_assignment_maps.py index 3043ee21b1..86d6bed53f 100755 --- a/doc/examples/bundle_assignment_maps.py +++ b/doc/examples/bundle_assignment_maps.py @@ -57,11 +57,13 @@ # # Creating 100 bundle assignment maps on AF_L using BUAN [Chandio2020]_ +rng = np.random.default_rng() + n = 100 indx = assignment_map(model_af_l, model_af_l, n) indx = np.array(indx) -colors = [np.random.rand(3) for si in range(n)] +colors = [rng.random(3) for si in range(n)] disks_color = [] for i in range(len(indx)): diff --git a/doc/examples/bundle_shape_similarity.py b/doc/examples/bundle_shape_similarity.py index ccb767e946..adcc198221 100755 --- a/doc/examples/bundle_shape_similarity.py +++ b/doc/examples/bundle_shape_similarity.py @@ -33,7 +33,7 @@ # Let's create two streamline sets (bundles) from same bundle cb_subj1 by # randomly selecting 60 streamlines two times. -rng = np.random.RandomState() +rng = np.random.default_rng() bundle1 = select_random_set_of_streamlines(cb_subj1, 60, rng=None) bundle2 = select_random_set_of_streamlines(cb_subj1, 60, rng=None) diff --git a/doc/examples/contextual_enhancement.py b/doc/examples/contextual_enhancement.py index 75b1845daf..e19d065d86 100644 --- a/doc/examples/contextual_enhancement.py +++ b/doc/examples/contextual_enhancement.py @@ -108,8 +108,9 @@ # Add Rician noise b0_slice = data[:, :, :, 1] b0_mask, mask = median_otsu(b0_slice) -np.random.seed(1) -data_noisy = add_noise(data, 10.0, np.mean(b0_slice[mask]), noise_type='rician') +rng = np.random.default_rng(1) +data_noisy = add_noise(data, 10.0, np.mean(b0_slice[mask]), + noise_type='rician', rng=rng) # Select a small part of it. padding = 3 # Include a larger region to avoid boundary effects diff --git a/doc/examples/fiber_to_bundle_coherence.py b/doc/examples/fiber_to_bundle_coherence.py index e8be88a9c2..764f357868 100644 --- a/doc/examples/fiber_to_bundle_coherence.py +++ b/doc/examples/fiber_to_bundle_coherence.py @@ -60,7 +60,7 @@ # Enables/disables interactive visualization interactive = False # Fix seed -np.random.seed(1) +rng = np.random.default_rng(1) # Read data hardi_fname, hardi_bval_fname, hardi_bvec_fname = get_fnames('stanford_hardi') diff --git a/doc/examples/gradients_spheres.py b/doc/examples/gradients_spheres.py index e81dfa30eb..3580d9b258 100644 --- a/doc/examples/gradients_spheres.py +++ b/doc/examples/gradients_spheres.py @@ -33,9 +33,10 @@ # We can first create some random points on a ``HemiSphere`` using spherical # polar coordinates. +rng = np.random.default_rng() n_pts = 64 -theta = np.pi * np.random.rand(n_pts) -phi = 2 * np.pi * np.random.rand(n_pts) +theta = np.pi * rng.random(n_pts) +phi = 2 * np.pi * rng.random(n_pts) hsph_initial = HemiSphere(theta=theta, phi=phi) ############################################################################### diff --git a/doc/examples/kfold_xval.py b/doc/examples/kfold_xval.py index e6adee88ef..5824ec657a 100644 --- a/doc/examples/kfold_xval.py +++ b/doc/examples/kfold_xval.py @@ -80,10 +80,12 @@ # the model will be fit to half of the data, and used to predict the other # half. -dti_cc = xval.kfold_xval(dti_model, cc_vox, 2) -csd_cc = xval.kfold_xval(csd_model, cc_vox, 2, response) -dti_cso = xval.kfold_xval(dti_model, cso_vox, 2) -csd_cso = xval.kfold_xval(csd_model, cso_vox, 2, response) +rng = np.random.default_rng(2014) + +dti_cc = xval.kfold_xval(dti_model, cc_vox, 2, rng=rng) +csd_cc = xval.kfold_xval(csd_model, cc_vox, 2, response, rng=rng) +dti_cso = xval.kfold_xval(dti_model, cso_vox, 2, rng=rng) +csd_cso = xval.kfold_xval(csd_model, cso_vox, 2, response, rng=rng) ############################################################################### # We plot a scatter plot of the data with the model predictions in each of diff --git a/doc/examples/reconst_msdki.py b/doc/examples/reconst_msdki.py index a89e3f2c98..aa38fcc148 100644 --- a/doc/examples/reconst_msdki.py +++ b/doc/examples/reconst_msdki.py @@ -78,9 +78,11 @@ # three different b-values). # Sample the spherical coordinates of 60 random diffusion-weighted directions. +rng = np.random.default_rng() + n_pts = 60 -theta = np.pi * np.random.rand(n_pts) -phi = 2 * np.pi * np.random.rand(n_pts) +theta = np.pi * rng.random(n_pts) +phi = 2 * np.pi * rng.random(n_pts) # Convert direction to cartesian coordinates. hsph_initial = HemiSphere(theta=theta, phi=phi) diff --git a/doc/examples/reconst_sh.py b/doc/examples/reconst_sh.py index dbd15dec76..548dda4cec 100644 --- a/doc/examples/reconst_sh.py +++ b/doc/examples/reconst_sh.py @@ -21,9 +21,10 @@ # We can first create some random points on a ``HemiSphere`` using spherical # polar coordinates. +rng = np.random.default_rng() n_pts = 64 -theta = np.pi * np.random.rand(n_pts) -phi = 2 * np.pi * np.random.rand(n_pts) +theta = np.pi * rng.random(n_pts) +phi = 2 * np.pi * rng.random(n_pts) hsph_initial = HemiSphere(theta=theta, phi=phi) ############################################################################### diff --git a/doc/examples/simulate_multi_tensor.py b/doc/examples/simulate_multi_tensor.py index 3a6c17457a..10f91aa7a0 100644 --- a/doc/examples/simulate_multi_tensor.py +++ b/doc/examples/simulate_multi_tensor.py @@ -21,9 +21,11 @@ # b-vectors. To create one, we can first create some random points on a # ``HemiSphere`` using spherical polar coordinates. +rng = np.random.default_rng() + n_pts = 64 -theta = np.pi * np.random.rand(n_pts) -phi = 2 * np.pi * np.random.rand(n_pts) +theta = np.pi * rng.random(size=n_pts) +phi = 2 * np.pi * rng.random(size=n_pts) hsph_initial = HemiSphere(theta=theta, phi=phi) ############################################################################### diff --git a/doc/examples/streamline_length.py b/doc/examples/streamline_length.py index 553678da94..3e39a790bb 100644 --- a/doc/examples/streamline_length.py +++ b/doc/examples/streamline_length.py @@ -27,21 +27,23 @@ def simulated_bundles(no_streamlines=50, n_pts=100): - t = np.linspace(-10, 10, n_pts) + rng = np.random.default_rng() - bundle = [] - for i in np.linspace(3, 5, no_streamlines): - pts = np.vstack((np.cos(2 * t/np.pi), np.zeros(t.shape) + i, t)).T - bundle.append(pts) + t = np.linspace(-10, 10, n_pts) - start = np.random.randint(10, 30, no_streamlines) - end = np.random.randint(60, 100, no_streamlines) + bundle = [] + for i in np.linspace(3, 5, no_streamlines): + pts = np.vstack((np.cos(2 * t/np.pi), np.zeros(t.shape) + i, t )).T + bundle.append(pts) - bundle = [10 * streamline[start[i]:end[i]] - for (i, streamline) in enumerate(bundle)] - bundle = [np.ascontiguousarray(streamline) for streamline in bundle] + start = rng.integers(10, 30, no_streamlines) + end = rng.integers(60, 100, no_streamlines) - return bundle + bundle = [10 * streamline[start[i]:end[i]] + for (i, streamline) in enumerate(bundle)] + bundle = [np.ascontiguousarray(streamline) for streamline in bundle] + + return bundle bundle = simulated_bundles() @@ -128,7 +130,7 @@ def simulated_bundles(no_streamlines=50, n_pts=100): scene.set_camera(position=(0, 0, 0), focal_point=(30, 0, 0)) window.record(scene, out_path='simulated_cosine_bundle.png', size=(900, 900)) if interactive: - window.show(scene) + window.show(scene) ############################################################################### # .. rst-class:: centered small fst-italic fw-semibold diff --git a/doc/examples/viz_bundles.py b/doc/examples/viz_bundles.py index bfc8faa8db..954da37b5a 100644 --- a/doc/examples/viz_bundles.py +++ b/doc/examples/viz_bundles.py @@ -190,9 +190,11 @@ scene.clear() -colors = [np.random.rand(*streamline.shape) for streamline in bundle_native] -stream_actor6 = actor.line(bundle_native, np.array(colors, dtype=object), - linewidth=0.2) +rng = np.random.default_rng() + +colors = [rng.random(streamline.shape) for streamline in bundle_native] + +stream_actor6 = actor.line(bundle_native, colors, linewidth=0.2) scene.add(stream_actor6)