handle var < 0 spotted by Ben

Author: jangorecki

inst/tests/froll.Rraw

@@ -1439,6 +1439,11 @@ test(6001.727, frollvar(adaptive=TRUE, c(1:2,NA), c(2,0,2), algo="exact"), c(NA_
 test(6001.728, frollvar(adaptive=TRUE, c(1:2,NA), c(2,0,2), algo="exact", na.rm=TRUE), c(NA_real_,NA_real_,NA_real_))
 test(6001.729, frollvar(adaptive=TRUE, c(1:2,NA), c(2,0,2), algo="exact", na.rm=TRUE, partial=TRUE), c(NA_real_,NA_real_,NA_real_))
 test(6001.730, frollvar(adaptive=TRUE, c(1:2,NA), c(2,0,2), fill=99, algo="exact", na.rm=TRUE), c(99,NA,NA))
+y = c(1e8+2.980232e-8, 1e8, 1e8, 1e8) # CLAMP0 test
+test(6001.731, frollvar(y, 3)[4L], 0)
+test(6001.732, frollsd(y, 3)[4L], 0)
+test(6001.733, frollvar(y, c(3,3,3,3), adaptive=TRUE)[4L], 0)
+test(6001.734, frollsd(y, c(3,3,3,3), adaptive=TRUE)[4L], 0)
 test(6001.781, frollapply(FUN=var, 1:3, 0), c(NA_real_,NA_real_,NA_real_))
 test(6001.782, frollapply(FUN=var, 1:3, 0, fill=99), c(NA_real_,NA_real_,NA_real_))
 test(6001.783, frollapply(FUN=var, c(1:2,NA), 0), c(NA_real_,NA_real_,NA_real_))

src/froll.c

@@ -1137,6 +1137,8 @@ void frollprodExact(const double *x, uint64_t nx, ans_t *ans, int k, double fill
   }
 }
 
+#define CLAMP0(x) (((x) < 0) ? 0 : (x))
+
 /* fast rolling var - fast
   Welford's online algorithm
   numerically stable
@@ -1169,7 +1171,8 @@ void frollvarFast(const double *x, uint64_t nx, ans_t *ans, int k, double fill,
     long double delta = x[i] - wmean; // i==k-1
     wmean += delta / (i + 1);
     m2 += delta * (x[i] - wmean);
-    ans->dbl_v[i] = (double)(m2 / k0);
+    double ans_i = m2 / k0;
+    ans->dbl_v[i] = CLAMP0(ans_i);
     if (R_FINITE((double) m2)) {
       for (uint64_t i=k; i<nx; i++) {
           double x_out = x[i-k];
@@ -1180,7 +1183,8 @@ void frollvarFast(const double *x, uint64_t nx, ans_t *ans, int k, double fill,
           long double delta_in = x_in - wmean;
           wmean += delta_in / k;
           m2 += delta_in * (x_in - wmean);
-          ans->dbl_v[i] = (double)(m2 / k0);
+          double ans_i = m2 / k0;
+          ans->dbl_v[i] = CLAMP0(ans_i);
       }
       if (!R_FINITE((double) m2)) {
         if (hasnf==-1)
@@ -1235,7 +1239,7 @@ void frollvarExact(const double *x, uint64_t nx, ans_t *ans, int k, double fill,
       if (!R_FINITE((double) wsum)) {
         if (ISNAN((double) wsum)) {
           if (!narm) {
-            ans->dbl_v[i] = (double) wsum;
+            ans->dbl_v[i] = (double) wsum; // propagate NAs
           }
           truehasnf = true;
         } else {
@@ -1249,7 +1253,8 @@ void frollvarExact(const double *x, uint64_t nx, ans_t *ans, int k, double fill,
           xi = x[i+j] - wmean;
           wsumxi += (xi * xi);
         }
-        ans->dbl_v[i] = (double) (wsumxi / (k - 1));
+        double ans_i = wsumxi / (k - 1);
+        ans->dbl_v[i] = CLAMP0(ans_i);
       }
     }
     if (truehasnf) {
@@ -1284,7 +1289,8 @@ void frollvarExact(const double *x, uint64_t nx, ans_t *ans, int k, double fill,
             xi = x[i+j] - wmean;
             wsumxi += (xi * xi);
           }
-          ans->dbl_v[i] = (double) (wsumxi / (k - 1));
+          double ans_i = wsumxi / (k - 1);
+          ans->dbl_v[i] = CLAMP0(ans_i);
         } else if (nc < (k - 1)) { // var(scalar) is also NA thats why k-1 so at least 2 numbers must be there
           long double wmean = wsum / (k - nc);
           long double xi = 0.0;
@@ -1295,7 +1301,8 @@ void frollvarExact(const double *x, uint64_t nx, ans_t *ans, int k, double fill,
               wsumxi += (xi * xi);
             }
           }
-          ans->dbl_v[i] = (double) (wsumxi / (k - nc - 1));
+          double ans_i = wsumxi / (k - nc - 1);
+          ans->dbl_v[i] = CLAMP0(ans_i);
         } else {
           ans->dbl_v[i] = NA_REAL;
         }

src/frolladaptive.c

@@ -737,6 +737,8 @@ void frolladaptiveprodExact(const double *x, uint64_t nx, ans_t *ans, const int
   }
 }
 
+#define CLAMP0(x) (((x) < 0) ? 0 : (x))
+
 /* fast rolling adaptive var - exact
  */
 void frolladaptivevarExact(const double *x, uint64_t nx, ans_t *ans, const int *k, double fill, bool narm, int hasnf, bool verbose) {
@@ -761,7 +763,7 @@ void frolladaptivevarExact(const double *x, uint64_t nx, ans_t *ans, const int *
         if (!R_FINITE((double) wsum)) {
           if (ISNAN((double) wsum)) {
             if (!narm) {
-              ans->dbl_v[i] = (double) wsum;
+              ans->dbl_v[i] = (double) wsum; // propagate NAs
             }
             truehasnf = true;
           } else {
@@ -775,7 +777,8 @@ void frolladaptivevarExact(const double *x, uint64_t nx, ans_t *ans, const int *
             xi = x[i+j] - wmean;
             wsumxi += (xi * xi);
           }
-          ans->dbl_v[i] = (double) (wsumxi / (k[i] - 1));
+          double ans_i = (wsumxi / (k[i] - 1));
+          ans->dbl_v[i] = CLAMP0(ans_i);
         }
       }
     }