Use the seed in all rng in blending code

pysteps/blending/steps.py

@@ -176,9 +176,10 @@ class StepsBlendingState:
     # Mapping from NWP members to ensemble members
     mapping_list_NWP_member_to_ensemble_member: np.ndarray | None = None
 
-    # Random states for precipitation and motion
+    # Random states for precipitation, motion and probmatching
     randgen_precip: list[np.random.RandomState] | None = None
     randgen_motion: list[np.random.RandomState] | None = None
+    randgen_probmatching: list[np.random.RandomState] | None = None
 
     # Variables for final forecast computation
     previous_displacement: list[Any] | None = None
@@ -1095,19 +1096,28 @@ def __multiply_precip_cascade_to_match_ensemble_members(self):
     def __initialize_random_generators(self):
         # 6. Initialize all the random generators and prepare for the forecast loop
         """Initialize all the random generators."""
+        seed = self.__config.seed
         if self.__config.noise_method is not None:
             self.__state.randgen_precip = []
-            self.__state.randgen_motion = []
-            seed = self.__config.seed
             for j in range(self.__config.n_ens_members):
                 rs = np.random.RandomState(seed)
                 self.__state.randgen_precip.append(rs)
                 seed = rs.randint(0, high=1e9)
+
+        if self.__config.probmatching_method is not None:
+            self.__state.randgen_probmatching = []
+            for j in range(self.__config.n_ens_members):
                 rs = np.random.RandomState(seed)
-                self.__state.randgen_motion.append(rs)
+                self.__state.randgen_probmatching.append(rs)
                 seed = rs.randint(0, high=1e9)
 
         if self.__config.velocity_perturbation_method is not None:
+            self.__state.randgen_motion = []
+            for j in range(self.__config.n_ens_members):
+                rs = np.random.RandomState(seed)
+                self.__state.randgen_motion.append(rs)
+                seed = rs.randint(0, high=1e9)
+
             (
                 init_velocity_noise,
                 self.__params.generate_velocity_noise,
@@ -2373,6 +2383,7 @@ def __post_process_output(
                     first_array=arr1,
                     second_array=arr2,
                     probability_first_array=weights_probability_matching_normalized[0],
+                    randgen=self.__state.randgen_probmatching[j],
                 )
             )
         else:

pysteps/noise/utils.py

@@ -101,8 +101,9 @@ def compute_noise_stddev_adjs(
     randstates = []
 
     for k in range(num_iter):
-        randstates.append(np.random.RandomState(seed=seed))
-        seed = np.random.randint(0, high=1e9)
+        rs = np.random.RandomState(seed=seed)
+        randstates.append(rs)
+        seed = rs.randint(0, high=1e9)
 
     def worker(k):
         # generate Gaussian white noise field, filter it using the chosen

pysteps/postprocessing/probmatching.py

@@ -274,7 +274,7 @@ def _get_error(scale):
     return shift, scale, R.reshape(shape)
 
 
-def resample_distributions(first_array, second_array, probability_first_array):
+def resample_distributions(first_array, second_array, probability_first_array, randgen):
     """
     Merges two distributions (e.g., from the extrapolation nowcast and NWP in the blending module)
     to effectively combine two distributions for probability matching without losing extremes.
@@ -324,7 +324,7 @@ def resample_distributions(first_array, second_array, probability_first_array):
     n = asort.shape[0]
 
     # Resample the distributions
-    idxsamples = np.random.binomial(1, probability_first_array, n).astype(bool)
+    idxsamples = randgen.binomial(1, probability_first_array, n).astype(bool)
     csort = np.where(idxsamples, asort, bsort)
     csort = np.sort(csort)[::-1]
 

pysteps/tests/test_blending_steps.py

@@ -8,48 +8,48 @@
 import pysteps
 from pysteps import blending, cascade
 
+# fmt:off
 steps_arg_values = [
-    # Test the case where both the radar image and the NWP fields contain no rain.
-    (1, 3, 4, 8, None, None, False, "spn", True, 4, False, False, 0, False),
-    (1, 3, 4, 8, "obs", None, False, "spn", True, 4, False, False, 0, False),
-    (1, 3, 4, 8, "incremental", None, False, "spn", True, 4, False, False, 0, False),
-    (1, 3, 4, 8, None, "mean", False, "spn", True, 4, False, False, 0, False),
-    (1, 3, 4, 8, None, "mean", False, "spn", True, 4, False, False, 0, True),
-    (1, 3, 4, 8, None, "cdf", False, "spn", True, 4, False, False, 0, False),
-    (1, [1, 2, 3], 4, 8, None, "cdf", False, "spn", True, 4, False, False, 0, False),
-    (1, 3, 4, 8, "incremental", "cdf", False, "spn", True, 4, False, False, 0, False),
-    (1, 3, 4, 6, "incremental", "cdf", False, "bps", True, 4, False, False, 0, False),
-    (1, 3, 4, 6, "incremental", "cdf", False, "bps", False, 4, False, False, 0, False),
-    (1, 3, 4, 6, "incremental", "cdf", False, "bps", False, 4, False, False, 0, True),
-    (1, 3, 4, 9, "incremental", "cdf", False, "spn", True, 4, False, False, 0, False),
-    (2, 3, 10, 8, "incremental", "cdf", False, "spn", True, 10, False, False, 0, False),
-    (5, 3, 5, 8, "incremental", "cdf", False, "spn", True, 5, False, False, 0, False),
-    (1, 10, 1, 8, "incremental", "cdf", False, "spn", True, 1, False, False, 0, False),
-    (2, 3, 2, 8, "incremental", "cdf", True, "spn", True, 2, False, False, 0, False),
-    # TODO: make next test work! This is currently not working on the main branch
-    # (2, 3, 4, 8, "incremental", "cdf", True, "spn", True, 2, False, False, 0, False),
-    # (2, 3, 4, 8, "incremental", "cdf", False, "spn", True, 2, False, False, 0, False),
-    (1, 3, 6, 8, None, None, False, "spn", True, 6, False, False, 0, False),
+    (1, 3, 4, 8, None, None, False, "spn", True, 4, False, False, 0, False, None),
+    (1, 3, 4, 8, "obs", None, False, "spn", True, 4, False, False, 0, False, None),
+    (1, 3, 4, 8, "incremental", None, False, "spn", True, 4, False, False, 0, False, None),
+    (1, 3, 4, 8, None, "mean", False, "spn", True, 4, False, False, 0, False, None),
+    (1, 3, 4, 8, None, "mean", False, "spn", True, 4, False, False, 0, True, None),
+    (1, 3, 4, 8, None, "cdf", False, "spn", True, 4, False, False, 0, False, None),
+    (1, [1, 2, 3], 4, 8, None, "cdf", False, "spn", True, 4, False, False, 0, False, None),
+    (1, 3, 4, 8, "incremental", "cdf", False, "spn", True, 4, False, False, 0, False, None),
+    (1, 3, 4, 6, "incremental", "cdf", False, "bps", True, 4, False, False, 0, False, None),
+    (1, 3, 4, 6, "incremental", "cdf", False, "bps", False, 4, False, False, 0, False, None),
+    (1, 3, 4, 6, "incremental", "cdf", False, "bps", False, 4, False, False, 0, True, None),
+    (1, 3, 4, 9, "incremental", "cdf", False, "spn", True, 4, False, False, 0, False, None),
+    (2, 3, 10, 8, "incremental", "cdf", False, "spn", True, 10, False, False, 0, False, None),
+    (5, 3, 5, 8, "incremental", "cdf", False, "spn", True, 5, False, False, 0, False, None),
+    (1, 10, 1, 8, "incremental", "cdf", False, "spn", True, 1, False, False, 0, False, None),
+    (2, 3, 2, 8, "incremental", "cdf", True, "spn", True, 2, False, False, 0, False, None),
+    (1, 3, 6, 8, None, None, False, "spn", True, 6, False, False, 0, False, None),
+    (1, 3, 6, 8, None, None, False, "spn", True, 6, False, False, 0, False, "bps"),
     #    Test the case where the radar image contains no rain.
-    (1, 3, 6, 8, None, None, False, "spn", True, 6, True, False, 0, False),
-    (5, 3, 5, 6, "incremental", "cdf", False, "spn", False, 5, True, False, 0, False),
-    (5, 3, 5, 6, "incremental", "cdf", False, "spn", False, 5, True, False, 0, True),
+    (1, 3, 6, 8, None, None, False, "spn", True, 6, True, False, 0, False, None),
+    (5, 3, 5, 6, "incremental", "cdf", False, "spn", False, 5, True, False, 0, False, None),
+    (5, 3, 5, 6, "incremental", "cdf", False, "spn", False, 5, True, False, 0, True, None),
     #   Test the case where the NWP fields contain no rain.
-    (1, 3, 6, 8, None, None, False, "spn", True, 6, False, True, 0, False),
-    (5, 3, 5, 6, "incremental", "cdf", False, "spn", False, 5, False, True, 0, True),
-    (1, 3, 6, 8, None, None, False, "spn", True, 6, True, True, 0, False),
-    (5, 3, 5, 6, "incremental", "cdf", False, "spn", False, 5, True, True, 0, False),
-    (5, 3, 5, 6, "obs", "mean", True, "spn", True, 5, True, True, 0, False),
+    (1, 3, 6, 8, None, None, False, "spn", True, 6, False, True, 0, False, None),
+    (5, 3, 5, 6, "incremental", "cdf", False, "spn", False, 5, False, True, 0, True, None),
+    # Test the case where both the radar image and the NWP fields contain no rain.
+    (1, 3, 6, 8, None, None, False, "spn", True, 6, True, True, 0, False, None),
+    (5, 3, 5, 6, "incremental", "cdf", False, "spn", False, 5, True, True, 0, False, None),
+    (5, 3, 5, 6, "obs", "mean", True, "spn", True, 5, True, True, 0, False, None),
     # Test for smooth radar mask
-    (1, 3, 6, 8, None, None, False, "spn", True, 6, False, False, 80, False),
-    (5, 3, 5, 6, "incremental", "cdf", False, "spn", False, 5, False, False, 80, False),
-    (5, 3, 5, 6, "obs", "mean", False, "spn", False, 5, False, False, 80, False),
-    (1, 3, 6, 8, None, None, False, "spn", True, 6, False, True, 80, False),
-    (5, 3, 5, 6, "incremental", "cdf", False, "spn", False, 5, True, False, 80, True),
-    (5, 3, 5, 6, "obs", "mean", False, "spn", False, 5, True, True, 80, False),
-    (5, [1, 2, 3], 5, 6, "obs", "mean", False, "spn", False, 5, True, True, 80, False),
-    (5, [1, 3], 5, 6, "obs", "mean", False, "spn", False, 5, True, True, 80, False),
+    (1, 3, 6, 8, None, None, False, "spn", True, 6, False, False, 80, False, None),
+    (5, 3, 5, 6, "incremental", "cdf", False, "spn", False, 5, False, False, 80, False, None),
+    (5, 3, 5, 6, "obs", "mean", False, "spn", False, 5, False, False, 80, False, None),
+    (1, 3, 6, 8, None, None, False, "spn", True, 6, False, True, 80, False, None),
+    (5, 3, 5, 6, "incremental", "cdf", False, "spn", False, 5, True, False, 80, True, None),
+    (5, 3, 5, 6, "obs", "mean", False, "spn", False, 5, True, True, 80, False, None),
+    (5, [1, 2, 3], 5, 6, "obs", "mean", False, "spn", False, 5, True, True, 80, False, None),
+    (5, [1, 3], 5, 6, "obs", "mean", False, "spn", False, 5, True, True, 80, False, None),
 ]
+# fmt:on
 
 steps_arg_names = (
     "n_models",
@@ -66,6 +66,7 @@
     "zero_nwp",
     "smooth_radar_mask_range",
     "resample_distribution",
+    "vel_pert_method",
 )
 
 
@@ -85,6 +86,7 @@ def test_steps_blending(
     zero_nwp,
     smooth_radar_mask_range,
     resample_distribution,
+    vel_pert_method,
 ):
     pytest.importorskip("cv2")
 
@@ -278,7 +280,7 @@ def test_steps_blending(
         noise_method="nonparametric",
         noise_stddev_adj="auto",
         ar_order=2,
-        vel_pert_method=None,
+        vel_pert_method=vel_pert_method,
         weights_method=weights_method,
         conditional=False,
         probmatching_method=probmatching_method,

pysteps/tests/test_postprocessing_probmatching.py

@@ -1,7 +1,10 @@
 import numpy as np
 import pytest
-from pysteps.postprocessing.probmatching import resample_distributions
-from pysteps.postprocessing.probmatching import nonparam_match_empirical_cdf
+
+from pysteps.postprocessing.probmatching import (
+    nonparam_match_empirical_cdf,
+    resample_distributions,
+)
 
 
 class TestResampleDistributions:
@@ -16,7 +19,7 @@ def test_valid_inputs(self):
         second_array = np.array([2, 4, 6, 8, 10])
         probability_first_array = 0.6
         result = resample_distributions(
-            first_array, second_array, probability_first_array
+            first_array, second_array, probability_first_array, np.random
         )
         expected_result = np.array([9, 8, 6, 3, 1])  # Expected result based on the seed
         assert result.shape == first_array.shape
@@ -27,7 +30,7 @@ def test_probability_zero(self):
         second_array = np.array([2, 4, 6, 8, 10])
         probability_first_array = 0.0
         result = resample_distributions(
-            first_array, second_array, probability_first_array
+            first_array, second_array, probability_first_array, np.random
         )
         assert np.array_equal(result, np.sort(second_array)[::-1])
 
@@ -36,7 +39,7 @@ def test_probability_one(self):
         second_array = np.array([2, 4, 6, 8, 10])
         probability_first_array = 1.0
         result = resample_distributions(
-            first_array, second_array, probability_first_array
+            first_array, second_array, probability_first_array, np.random
         )
         assert np.array_equal(result, np.sort(first_array)[::-1])
 
@@ -45,7 +48,7 @@ def test_nan_in_arr1_prob_1(self):
         array_without_nan = np.array([2.0, 4, 6, 8, 10])
         probability_first_array = 1.0
         result = resample_distributions(
-            array_with_nan, array_without_nan, probability_first_array
+            array_with_nan, array_without_nan, probability_first_array, np.random
         )
         expected_result = np.array([np.nan, 9, 7, 3, 1], dtype=float)
         assert np.allclose(result, expected_result, equal_nan=True)
@@ -55,7 +58,7 @@ def test_nan_in_arr1_prob_0(self):
         array_without_nan = np.array([2, 4, 6, 8, 10])
         probability_first_array = 0.0
         result = resample_distributions(
-            array_with_nan, array_without_nan, probability_first_array
+            array_with_nan, array_without_nan, probability_first_array, np.random
         )
         expected_result = np.array([np.nan, 10, 8, 4, 2], dtype=float)
         assert np.allclose(result, expected_result, equal_nan=True)
@@ -65,7 +68,7 @@ def test_nan_in_arr2_prob_1(self):
         array_with_nan = np.array([2.0, 4, 6, np.nan, 10])
         probability_first_array = 1.0
         result = resample_distributions(
-            array_without_nan, array_with_nan, probability_first_array
+            array_without_nan, array_with_nan, probability_first_array, np.random
         )
         expected_result = np.array([np.nan, 9, 5, 3, 1], dtype=float)
         assert np.allclose(result, expected_result, equal_nan=True)
@@ -75,7 +78,7 @@ def test_nan_in_arr2_prob_0(self):
         array_with_nan = np.array([2, 4, 6, np.nan, 10])
         probability_first_array = 0.0
         result = resample_distributions(
-            array_without_nan, array_with_nan, probability_first_array
+            array_without_nan, array_with_nan, probability_first_array, np.random
         )
         expected_result = np.array([np.nan, 10, 6, 4, 2], dtype=float)
         assert np.allclose(result, expected_result, equal_nan=True)
@@ -85,7 +88,7 @@ def test_nan_in_both_prob_1(self):
         array2_with_nan = np.array([2.0, 4, np.nan, np.nan, 10])
         probability_first_array = 1.0
         result = resample_distributions(
-            array1_with_nan, array2_with_nan, probability_first_array
+            array1_with_nan, array2_with_nan, probability_first_array, np.random
         )
         expected_result = np.array([np.nan, np.nan, np.nan, 9, 1], dtype=float)
         assert np.allclose(result, expected_result, equal_nan=True)
@@ -95,7 +98,7 @@ def test_nan_in_both_prob_0(self):
         array2_with_nan = np.array([2.0, 4, np.nan, np.nan, 10])
         probability_first_array = 0.0
         result = resample_distributions(
-            array1_with_nan, array2_with_nan, probability_first_array
+            array1_with_nan, array2_with_nan, probability_first_array, np.random
         )
         expected_result = np.array([np.nan, np.nan, np.nan, 10, 2], dtype=float)
         assert np.allclose(result, expected_result, equal_nan=True)