diff --git a/pysteps/blending/steps.py b/pysteps/blending/steps.py index a6ecc09f0..d15208064 100644 --- a/pysteps/blending/steps.py +++ b/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( diff --git a/pysteps/noise/utils.py b/pysteps/noise/utils.py index 58495b8ad..aaae82f9b 100644 --- a/pysteps/noise/utils.py +++ b/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 diff --git a/pysteps/postprocessing/probmatching.py b/pysteps/postprocessing/probmatching.py index afc951542..efda2b103 100644 --- a/pysteps/postprocessing/probmatching.py +++ b/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] diff --git a/pysteps/tests/test_postprocessing_probmatching.py b/pysteps/tests/test_postprocessing_probmatching.py index 8c7c12f4f..b42a95e7e 100644 --- a/pysteps/tests/test_postprocessing_probmatching.py +++ b/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)