Use seed for all rng in blending to make a test run completely determ…

…inistic

pysteps/blending/steps.py

@@ -734,16 +734,19 @@ def forecast(
         )
 
         # 6. Initialize all the random generators and prepare for the forecast loop
-        randgen_prec, vps, generate_vel_noise = _init_random_generators(
-            velocity,
-            noise_method,
-            vel_pert_method,
-            vp_par,
-            vp_perp,
-            seed,
-            n_ens_members,
-            kmperpixel,
-            timestep,
+        randgen_prec, vps, generate_vel_noise, randgen_probmatching = (
+            _init_random_generators(
+                velocity,
+                noise_method,
+                probmatching_method,
+                vel_pert_method,
+                vp_par,
+                vp_perp,
+                seed,
+                n_ens_members,
+                kmperpixel,
+                timestep,
+            )
         )
         D, D_Yn, D_pb, R_f, R_m, mask_rim, struct, fft_objs = _prepare_forecast_loop(
             precip_cascade,
@@ -1621,6 +1624,7 @@ def worker(j):
                                 first_array=arr1,
                                 second_array=arr2,
                                 probability_first_array=weights_pm_normalized[0],
+                                randgen=randgen_probmatching[j],
                             )
                         else:
                             R_pm_resampled = R_pm_blended.copy()
@@ -2290,6 +2294,7 @@ def _find_nwp_combination(
 def _init_random_generators(
     velocity,
     noise_method,
+    probmatching_method,
     vel_pert_method,
     vp_par,
     vp_perp,
@@ -2301,16 +2306,28 @@ def _init_random_generators(
     """Initialize all the random generators."""
     if noise_method is not None:
         randgen_prec = []
-        randgen_motion = []
         for j in range(n_ens_members):
             rs = np.random.RandomState(seed)
             randgen_prec.append(rs)
             seed = rs.randint(0, high=1e9)
+    else:
+        randgen_prec = None
+
+    if probmatching_method is not None:
+        randgen_probmatching = []
+        for j in range(n_ens_members):
             rs = np.random.RandomState(seed)
-            randgen_motion.append(rs)
+            randgen_probmatching.append(rs)
             seed = rs.randint(0, high=1e9)
+    else:
+        randgen_probmatching = None
 
     if vel_pert_method is not None:
+        randgen_motion = []
+        for j in range(n_ens_members):
+            rs = np.random.RandomState(seed)
+            randgen_motion.append(rs)
+            seed = rs.randint(0, high=1e9)
         init_vel_noise, generate_vel_noise = noise.get_method(vel_pert_method)
 
         # initialize the perturbation generators for the motion field
@@ -2326,7 +2343,7 @@ def _init_random_generators(
     else:
         vps, generate_vel_noise = None, None
 
-    return randgen_prec, vps, generate_vel_noise
+    return randgen_prec, vps, generate_vel_noise, randgen_probmatching
 
 
 def _prepare_forecast_loop(

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_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)