Merge pull request #187 from NREL/gb/fix_qdm

fixing qdm again, this time with real test on real problematic data

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
 

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

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