From 7623f91709daca5373c8078892b3e4305e0bcfd5 Mon Sep 17 00:00:00 2001 From: grantbuster Date: Fri, 11 Oct 2024 16:41:23 -0600 Subject: [PATCH] fixing qdm again, this time with real test on real problematic data --- rex/bias_correction.py | 4 +-- rex/utilities/bc_utils.py | 24 ++++++++++++++- tests/test_bc.py | 63 ++++++++++++++++++++++++++++++++++++++- 3 files changed, 87 insertions(+), 4 deletions(-) diff --git a/rex/bias_correction.py b/rex/bias_correction.py index 93e555c7..d48f1eab 100644 --- a/rex/bias_correction.py +++ b/rex/bias_correction.py @@ -163,7 +163,7 @@ def qdm_irrad(ghi, dni, dhi, ghi_params_mf=None, dni_params_mf=None, dist='empirical', relative=True, sampling='linear', log_base=10, - delta_denom_min=None, delta_denom_zero=0.01): + delta_denom_min=0.01, delta_denom_zero=None): """Correct irradiance using the quantile delta mapping based on the method from Cannon et al., 2015 @@ -289,7 +289,7 @@ def qdm_irrad(ghi, dni, dhi, def qdm_ws(ws, params_oh, params_mh, params_mf=None, dist='empirical', relative=True, sampling='linear', log_base=10, - delta_denom_min=None, delta_denom_zero=0.01): + delta_denom_min=0.01, delta_denom_zero=None): """Correct windspeed using quantile delta mapping based on the method from Cannon et al., 2015 diff --git a/rex/utilities/bc_utils.py b/rex/utilities/bc_utils.py index f27c6010..daf32ca2 100644 --- a/rex/utilities/bc_utils.py +++ b/rex/utilities/bc_utils.py @@ -71,6 +71,28 @@ def sample_q_invlog(n_samples, log_base): return quantiles +def sample_cdf(quantiles, x_values, n_samples): + """Randomly draw a number of real values from a CDF. + + quantiles : np.ndarray + 1D array of quantile values from 0 to 1. Must be monotonic. + x_values : np.ndarray + Values on the x-axis of a CDF corresponding to quantiles. Must be + monotonic. + n_samples : int + Number of sample to draw + + Returns + ------- + samples : np.ndarray + 1D array of real values sampled from the CDF made up by quantiles and + x_values + """ + samples = np.random.uniform(0, 1, n_samples) + samples = np.interp(samples, quantiles, x_values) + return samples + + class QuantileDeltaMapping: """Class for quantile delta mapping based on the method from Cannon et al., 2015 @@ -376,7 +398,7 @@ def run_qdm(cls, arr, params_oh, params_mh, params_mf, if relative: if delta_denom_zero is not None: x_mh_mf[x_mh_mf == 0] = delta_denom_zero - elif delta_denom_min is not None: + if delta_denom_min is not None: x_mh_mf = np.maximum(x_mh_mf, delta_denom_min) delta = arr / x_mh_mf # Eq.4: x_m_p / F-1_m_h(Tau_m_p) arr_bc = x_oh * delta # Eq.6: x^_m_p = x^_o:m_h:p * delta diff --git a/tests/test_bc.py b/tests/test_bc.py index 4c2ac06a..2e210633 100644 --- a/tests/test_bc.py +++ b/tests/test_bc.py @@ -7,7 +7,8 @@ from flaky import flaky from rex.temporal_stats.temporal_stats import cdf -from rex.utilities.bc_utils import QuantileDeltaMapping +from rex.utilities.bc_utils import (QuantileDeltaMapping, sample_q_invlog, + sample_cdf) @flaky(max_runs=3, min_passes=1) @@ -98,3 +99,63 @@ def test_qdm_parallel(): arr_mf_bc_ser = qdm_rel_fut(arr_mf, max_workers=1) arr_mf_bc_par = qdm_rel_fut(arr_mf, max_workers=2) assert np.allclose(arr_mf_bc_ser, arr_mf_bc_par) + + +def test_difficult_qdm(): + """Relative QDM can blow up with extremely large deltas resulting from very + small values in the denominator. This test verifies the protection against + this using the ``delta_denom_min`` feature. The distributions are real + distributions from Sup3rCC EC-Earth3 SSP585 wind data""" + + params_mf = np.array([0.04, 1.51, 2.14, 2.64, 3.06, 3.44, 3.78, + 4.11, 4.43, 4.74, 5.03, 5.32, 5.6, 5.88, 6.15, + 6.42, 6.69, 6.95, 7.22, 7.48, 7.73, 7.99, 8.23, + 8.48, 8.72, 8.95, 9.2, 9.44, 9.67, 9.91, 10.15, + 10.39, 10.63, 10.88, 11.11, 11.35, 11.59, 11.84, + 12.1, 12.36, 12.63, 12.92, 13.23, 13.58, 13.96, + 14.42, 14.99, 15.75, 17.03, 25.48]) + + params_mh = np.array([0., 1.53, 2.19, 2.71, 3.15, 3.53, 3.89, 4.22, 4.53, + 4.83, 5.13, 5.42, 5.71, 5.98, 6.25, 6.53, 6.8, 7.06, + 7.3, 7.54, 7.79, 8.03, 8.27, 8.5, 8.74, 8.98, 9.21, + 9.46, 9.7, 9.94, 10.18, 10.41, 10.61, 10.87, 11.11, + 11.34, 11.56, 11.81, 12.06, 12.33, 12.61, 12.88, + 13.16, 13.49, 13.83, 14.23, 14.72, 15.41, 16.71, + 23.42]) + + params_oh = np.array([2.000e-02, 1.860e+00, 2.550e+00, 3.080e+00, + 3.530e+00, 3.950e+00, 4.330e+00, 4.690e+00, + 5.040e+00, 5.360e+00, 5.660e+00, 5.950e+00, + 6.220e+00, 6.490e+00, 6.740e+00, 6.990e+00, + 7.230e+00, 7.470e+00, 7.690e+00, 7.900e+00, + 8.120e+00, 8.340e+00, 8.550e+00, 8.770e+00, + 8.970e+00, 9.180e+00, 9.380e+00, 9.570e+00, + 9.770e+00, 9.960e+00, 1.016e+01, 1.035e+01, + 1.055e+01, 1.076e+01, 1.097e+01, 1.118e+01, + 1.140e+01, 1.163e+01, 1.188e+01, 1.211e+01, + 1.236e+01, 1.267e+01, 1.297e+01, 1.334e+01, + 1.377e+01, 1.427e+01, 1.479e+01, 1.547e+01, + 1.669e+01, 2.191e+01]) + + params_oh = params_oh[np.newaxis] # qdm expects (space, time) for params + params_mh = params_mh[np.newaxis] # qdm expects (space, time) for params + params_mf = params_mf[np.newaxis] # qdm expects (space, time) for params + + quantiles = sample_q_invlog(params_mf.shape[1], log_base=10) + arr = sample_cdf(quantiles, params_mf[0], int(1e4)) + arr[0] = 0.04 # trouble! + + # test bad result with large delta + qdm_rel_fut = QuantileDeltaMapping(params_oh, params_mh, params_mf, + dist='empirical', relative=True, + sampling='invlog', log_base=10) + arr_bc = qdm_rel_fut(arr[:, np.newaxis]) + assert arr_bc.max() > 1000 # bad result + + # test that delta_denom_min fixes this + qdm_rel_fut = QuantileDeltaMapping(params_oh, params_mh, params_mf, + dist='empirical', relative=True, + sampling='invlog', log_base=10, + delta_denom_min=0.01) + arr_bc = qdm_rel_fut(arr[:, np.newaxis]) + assert arr_bc.max() < 40 # fixed result