fix: NPY002

random number generator stuff.

This adds possibility to use rng properly in ModOpt.

src/modopt/math/matrix.py

@@ -274,6 +274,8 @@ class PowerMethod:
         initialisation (default is ``True``)
     verbose : bool, optional
         Optional verbosity (default is ``False``)
+    rng: int, xp.random.Generator or None (default is ``None``)
+        Random number generator or seed.
 
     Examples
     --------
@@ -300,6 +302,7 @@ def __init__(
         auto_run=True,
         compute_backend="numpy",
         verbose=False,
+        rng=None,
     ):
 
         self._operator = operator
@@ -308,6 +311,7 @@ def __init__(
         self._verbose = verbose
         xp, compute_backend = get_backend(compute_backend)
         self.xp = xp
+        self.rng = None
         self.compute_backend = compute_backend
         if auto_run:
             self.get_spec_rad()
@@ -324,7 +328,8 @@ def _set_initial_x(self):
             Random values of the same shape as the input data
 
         """
-        return self.xp.random.random(self._data_shape).astype(self._data_type)
+        rng = self.xp.random.default_rng(self.rng)
+        return rng.random(self._data_shape).astype(self._data_type)
 
     def get_spec_rad(self, tolerance=1e-6, max_iter=20, extra_factor=1.0):
         """Get spectral radius.

src/modopt/signal/noise.py

@@ -11,7 +11,7 @@
 from modopt.base.backend import get_array_module
 
 
-def add_noise(input_data, sigma=1.0, noise_type="gauss"):
+def add_noise(input_data, sigma=1.0, noise_type="gauss", rng=None):
     """Add noise to data.
 
     This method adds Gaussian or Poisson noise to the input data.
@@ -25,6 +25,9 @@ def add_noise(input_data, sigma=1.0, noise_type="gauss"):
         default is ``1.0``)
     noise_type : {'gauss', 'poisson'}
         Type of noise to be added (default is ``'gauss'``)
+    rng: np.random.Generator or int
+        A Random number generator or a seed to initialize one.
+
 
     Returns
     -------
@@ -64,6 +67,9 @@ def add_noise(input_data, sigma=1.0, noise_type="gauss"):
     array([ 3.24869073, -1.22351283, -1.0563435 , -2.14593724,  1.73081526])
 
     """
+    if not isinstance(rng, np.random.Generator):
+        rng = np.random.default_rng(rng)
+
     input_data = np.array(input_data)
 
     if noise_type not in {"gauss", "poisson"}:
@@ -78,10 +84,10 @@ def add_noise(input_data, sigma=1.0, noise_type="gauss"):
             )
 
     if noise_type == "gauss":
-        random = np.random.randn(*input_data.shape)
+        random = rng.standard_normal(input_data.shape)
 
     elif noise_type == "poisson":
-        random = np.random.poisson(np.abs(input_data))
+        random = rng.poisson(np.abs(input_data))
 
     if isinstance(sigma, (int, float)):
         return input_data + sigma * random

src/modopt/signal/validation.py

@@ -16,6 +16,7 @@ def transpose_test(
     x_args=None,
     y_shape=None,
     y_args=None,
+    rng=None,
 ):
     """Transpose test.
 
@@ -36,6 +37,8 @@ def transpose_test(
         Shape of transpose operator input data (default is ``None``)
     y_args : tuple, optional
         Arguments to be passed to transpose operator (default is ``None``)
+    rng: np.random.Generator or int or None (default is ``None``)
+        Initialized random number generator or seed.
 
     Raises
     ------
@@ -62,9 +65,11 @@ def transpose_test(
     if isinstance(y_args, type(None)):
         y_args = x_args
 
+    if not isinstance(rng, np.random.Generator):
+        rng = np.random.default_rng(rng)
     # Generate random arrays.
-    x_val = np.random.ranf(x_shape)
-    y_val = np.random.ranf(y_shape)
+    x_val = rng.random(x_shape)
+    y_val = rng.random(y_shape)
 
     # Calculate <MX, Y>
     mx_y = np.sum(np.multiply(operator(x_val, x_args), y_val))

tests/test_algorithms.py

@@ -24,6 +24,9 @@
     SKLEARN_AVAILABLE = False
 
 
+rng = np.random.default_rng()
+
+
 @fixture
 def idty():
     """Identity function."""
@@ -84,7 +87,7 @@ class AlgoCases:
     """
 
     data1 = np.arange(9).reshape(3, 3).astype(float)
-    data2 = data1 + np.random.randn(*data1.shape) * 1e-6
+    data2 = data1 + rng.standard_normal(data1.shape) * 1e-6
     max_iter = 20
 
     @parametrize(

tests/test_math.py

@@ -29,6 +29,8 @@
 else:
     SKIMAGE_AVAILABLE = True
 
+rng = np.random.default_rng(1)
+
 
 class TestConvolve:
     """Test convolve functions."""
@@ -136,18 +138,15 @@ class TestMatrix:
         ),
     )
 
-    @pytest.fixture
+    @pytest.fixture(scope="module")
     def pm_instance(self, request):
         """Power Method instance."""
-        np.random.seed(1)
         pm = matrix.PowerMethod(
             lambda x_val: x_val.dot(x_val.T),
             self.array33.shape,
-            auto_run=request.param,
             verbose=True,
+            rng=np.random.default_rng(0),
         )
-        if not request.param:
-            pm.get_spec_rad(max_iter=1)
         return pm
 
     @pytest.mark.parametrize(
@@ -195,12 +194,7 @@ def test_rotate(self):
 
         npt.assert_raises(ValueError, matrix.rotate, self.array23, np.pi / 2)
 
-    @pytest.mark.parametrize(
-        ("pm_instance", "value"),
-        [(True, 1.0), (False, 0.8675467477372257)],
-        indirect=["pm_instance"],
-    )
-    def test_power_method(self, pm_instance, value):
+    def test_power_method(self, pm_instance, value=1):
         """Test power method."""
         npt.assert_almost_equal(pm_instance.spec_rad, value)
         npt.assert_almost_equal(pm_instance.inv_spec_rad, 1 / value)

tests/test_signal.py

@@ -45,13 +45,13 @@ class TestNoise:
 
     data1 = np.arange(9).reshape(3, 3).astype(float)
     data2 = np.array(
-        [[0, 2.0, 2.0], [4.0, 5.0, 10], [11.0, 15.0, 18.0]],
+        [[0, 3.0, 4.0], [6.0, 9.0, 8.0], [14.0, 14.0, 17.0]],
     )
     data3 = np.array(
         [
-            [1.62434536, 0.38824359, 1.47182825],
-            [1.92703138, 4.86540763, 2.6984613],
-            [7.74481176, 6.2387931, 8.3190391],
+            [0.3455842, 1.8216181, 2.3304371],
+            [1.6968428, 4.9053559, 5.4463746],
+            [5.4630468, 7.5811181, 8.3645724],
         ]
     )
     data4 = np.array([[0, 0, 0], [0, 0, 5.0], [6.0, 7.0, 8.0]])
@@ -70,9 +70,10 @@ class TestNoise:
     )
     def test_add_noise(self, data, noise_type, sigma, data_noise):
         """Test add_noise."""
-        np.random.seed(1)
+        rng = np.random.default_rng(1)
         npt.assert_almost_equal(
-            noise.add_noise(data, sigma=sigma, noise_type=noise_type), data_noise
+            noise.add_noise(data, sigma=sigma, noise_type=noise_type, rng=rng),
+            data_noise,
         )
 
     @pytest.mark.parametrize(
@@ -241,13 +242,13 @@ class TestValidation:
 
     def test_transpose_test(self):
         """Test transpose_test."""
-        np.random.seed(2)
         npt.assert_equal(
             validation.transpose_test(
                 lambda x_val, y_val: x_val.dot(y_val),
                 lambda x_val, y_val: x_val.dot(y_val.T),
                 self.array33.shape,
                 x_args=self.array33,
+                rng=2,
             ),
             None,
         )