rounding error fixes (#281)

Author: s-scherrer

CHANGELOG.rst

@@ -5,6 +5,10 @@ Changelog
 Unreleased
 ==========
 
+Version 0.14.3, 2022-12-26
+==========================
+- fixing rounding error issues in metric calculations (PR `#281 <https://github.com/TUW-GEO/pytesmo/pull/281>`_)
+
 Version 0.14.2, 2022-12-14
 ==========================
 - small bug fixes/doc updates (PRs `#273 <https://github.com/TUW-GEO/pytesmo/pull/273>`_, `#275 <https://github.com/TUW-GEO/pytesmo/pull/275>`_, `#276 <https://github.com/TUW-GEO/pytesmo/pull/276>`_,  `#278 <https://github.com/TUW-GEO/pytesmo/pull/278>`_)

src/pytesmo/metrics/_fast.c

@@ -132,7 +132,7 @@ cpdef mse_corr(floating [:] x, floating [:] y):
 cpdef _mse_corr_from_moments(
     floating mx, floating my, floating vx, floating vy, floating cov
 ):
-    return 2 * sqrt(vx) * sqrt(vy) - 2 * cov
+    return max(2 * sqrt(vx) * sqrt(vy) - 2 * cov, 0)
 
 
 
@@ -322,7 +322,12 @@ cpdef _ubrmsd(floating [:] x, floating [:] y):
 cpdef _pearsonr_from_moments(floating varx, floating vary, floating cov, int n):
     cdef floating R, p_R, t_squared, df, z
 
-    R = cov / sqrt(varx * vary)
+    # not defined in this case
+    if varx == 0 or vary == 0:
+        return np.nan, np.nan
+
+    # due to rounding errors, values slightly above 1 are possible
+    R = max(-1, min(cov / sqrt(varx * vary), 1))
 
     # p-value for R
     if fabs(R) == 1.0:

src/pytesmo/metrics/_fast_pairwise.c

@@ -20,7 +20,8 @@
 # Failure in temporal matching
 TEMPORAL_MATCHING_FAILED = 3
 # Temporal matching returned without error, but the data we need is not
-# available
+# available. This can happen if there is no temporal overlap of the required
+# datasets, or if there is no data at all for one of the required datasets.
 NO_TEMP_MATCHED_DATA = 4
 # the scaling failed
 SCALING_FAILED = 5

src/pytesmo/metrics/_fast_pairwise.pyx

@@ -1473,7 +1473,7 @@ def _calc_pairwise_metrics(
         mse_corr = _mse_corr_from_moments(mx, my, varx, vary, cov)
         mse_var = _mse_var_from_moments(mx, my, varx, vary, cov)
         mse_bias = _mse_bias_from_moments(mx, my, varx, vary, cov)
-        mse = mse_corr + mse_var + mse_bias
+        mse = max(mse_corr + mse_var + mse_bias, 0)
         result["mse_corr" + suffix][0] = mse_corr
         result["mse_var" + suffix][0] = mse_var
         result["mse_bias" + suffix][0] = mse_bias