diff --git a/inst/tests/froll.Rraw b/inst/tests/froll.Rraw index def2b3e9f..cdd54cad8 100644 --- a/inst/tests/froll.Rraw +++ b/inst/tests/froll.Rraw @@ -1202,8 +1202,8 @@ test(6000.931, frollprod(1:3, 2), c(NA, 2, 6), output="frollprodFast: running fo test(6000.932, frollprod(1:3, 2, align="left"), c(2, 6, NA), output="frollfun: align") test(6000.933, frollprod(c(1,2,NA), 2), c(NA, 2, NA), output="non-finite values are present in input, re-running with extra care for NFs") test(6000.934, frollprod(c(NA,2,3), 2), c(NA, NA, 6), output="non-finite values are present in input, skip non-finite inaware attempt and run with extra care for NFs straighaway") -test(6000.935, frollprod(1:3, c(2,2,2), adaptive=TRUE), c(NA, 2, 6), output="frolladaptiveprodFast: running for input length") -test(6000.936, frollprod(c(NA,2,3), c(2,2,2), adaptive=TRUE), c(NA, NA, 6), output="non-finite values are present in input, re-running with extra care for NFs") +test(6000.935, frollprod(1:3, c(2,2,2), adaptive=TRUE), c(NA, 2, 6), output="algo 0 not implemented, fall back to 1") +test(6000.936, frollprod(c(NA,2,3), c(2,2,2), adaptive=TRUE), c(NA, NA, 6), output="non-finite values are present in input, na.rm=FALSE and algo='exact' propagates NFs properply, no need to re-run") options(datatable.verbose=FALSE) # floating point overflow test(6000.941, frollprod(c(1e100, 1e100, 1e100, 1e100, 1e100), 5), c(NA,NA,NA,NA,Inf)) @@ -1222,6 +1222,15 @@ test(6000.953, frollprod(c(1e100, 1e100, 1e100, 1e100, 1e100), rep(5, 5), algo=" test(6000.954, frollprod(c(1e100, 1e100, 1e100, 1e100, 1e100), rep(4, 5), algo="exact", adaptive=TRUE), c(NA,NA,NA,Inf,Inf)) test(6000.955, frollprod(c(1e100, 1e100, 1e100, 1e100, -1e100), rep(5, 5), algo="exact", adaptive=TRUE), c(NA,NA,NA,NA,-Inf)) test(6000.956, frollprod(c(1e100, 1e100, 1e100, 1e100, -1e100), rep(4, 5), algo="exact", adaptive=TRUE), c(NA,NA,NA,Inf,-Inf)) +# rolling product and numerical stability #7349 +test(6000.9601, frollprod(c(2,2,0,2,2), 2), c(NA,4,0,0,4)) +test(6000.9602, frollprod(c(2,2,0,-2,2), 2), c(NA,4,0,0,-4)) +test(6000.9603, frollprod(c(2,2,0,-2,-2), 2), c(NA,4,0,0,4)) +test(6000.9604, frollprod(c(2,2,0,Inf,2), 2), c(NA,4,0,NaN,Inf)) +test(6000.9605, frollprod(c(2,2,0,-Inf,2), 2), c(NA,4,0,NaN,-Inf)) +test(6000.9606, frollprod(c(2,2,0,-Inf,-2), 2), c(NA,4,0,NaN,Inf)) +test(6000.9607, frollprod(c(0,2,2,2,2), 2), c(NA,0,4,4,4)) +test(6000.9608, frollprod(c(2,0,2,2,2), 2), c(NA,0,0,4,4)) # n==0, k==0, k[i]==0 test(6001.111, frollmean(1:3, 0), c(NaN,NaN,NaN), options=c("datatable.verbose"=TRUE), output="window width of size 0") @@ -2242,7 +2251,7 @@ rollfun = function(x, n, FUN, fill=NA_real_, na.rm=FALSE, nf.rm=FALSE, partial=F ans } base_compare = function(x, n, funs=c("mean","sum","max","min","prod","median"), algos=c("fast","exact")) { - num.step = 0.001 + num.step = 0.0001 for (fun in funs) { for (na.rm in c(FALSE, TRUE)) { for (fill in c(NA_real_, 0)) { @@ -2276,6 +2285,39 @@ base_compare = function(x, n, funs=c("mean","sum","max","min","prod","median"), } } } +num = 7000.0 +x = rnorm(1e3); n = 50 +base_compare(x, n) +x = rnorm(1e3+1); n = 50 ## uneven len +base_compare(x, n) +x = rnorm(1e3); n = 51 ## uneven window +base_compare(x, n) +x = rnorm(1e3+1); n = 51 +base_compare(x, n) +x = sort(rnorm(1e3)); n = 50 ## inc +base_compare(x, n) +x = sort(rnorm(1e3+1)); n = 50 +base_compare(x, n) +x = sort(rnorm(1e3)); n = 51 +base_compare(x, n) +x = sort(rnorm(1e3+1)); n = 51 +base_compare(x, n) +x = rev(sort(rnorm(1e3))); n = 50 ## desc +base_compare(x, n) +x = rev(sort(rnorm(1e3+1))); n = 50 +base_compare(x, n) +x = rev(sort(rnorm(1e3))); n = 51 +base_compare(x, n) +x = rev(sort(rnorm(1e3+1))); n = 51 +base_compare(x, n) +x = rep(rnorm(1), 1e3); n = 50 ## const +base_compare(x, n) +x = rep(rnorm(1), 1e3+1); n = 50 +base_compare(x, n) +x = rep(rnorm(1), 1e3); n = 51 +base_compare(x, n) +x = rep(rnorm(1), 1e3+1); n = 51 +base_compare(x, n) num = 7100.0 ## random NA non-finite x = makeNA(rnorm(1e3), nf=TRUE); n = 50 @@ -2433,6 +2475,39 @@ afun_compare = function(x, n, funs=c("mean","sum","max","min","prod","median"), } } num = 7300.0 +x = rnorm(1e3); n = sample(50, length(x), TRUE) +afun_compare(x, n) +x = rnorm(1e3+1); n = sample(50, length(x), TRUE) ## uneven len +afun_compare(x, n) +x = rnorm(1e3); n = sample(51, length(x), TRUE) ## uneven window +afun_compare(x, n) +x = rnorm(1e3+1); n = sample(51, length(x), TRUE) +afun_compare(x, n) +x = sort(rnorm(1e3)); n = sample(50, length(x), TRUE) ## inc +afun_compare(x, n) +x = sort(rnorm(1e3+1)); n = sample(50, length(x), TRUE) +afun_compare(x, n) +x = sort(rnorm(1e3)); n = sample(51, length(x), TRUE) +afun_compare(x, n) +x = sort(rnorm(1e3+1)); n = sample(51, length(x), TRUE) +afun_compare(x, n) +x = rev(sort(rnorm(1e3))); n = sample(50, length(x), TRUE) ## desc +afun_compare(x, n) +x = rev(sort(rnorm(1e3+1))); n = sample(50, length(x), TRUE) +afun_compare(x, n) +x = rev(sort(rnorm(1e3))); n = sample(51, length(x), TRUE) +afun_compare(x, n) +x = rev(sort(rnorm(1e3+1))); n = sample(51, length(x), TRUE) +afun_compare(x, n) +x = rep(rnorm(1), 1e3); n = sample(50, length(x), TRUE) ## const +afun_compare(x, n) +x = rep(rnorm(1), 1e3+1); n = sample(50, length(x), TRUE) +afun_compare(x, n) +x = rep(rnorm(1), 1e3); n = sample(51, length(x), TRUE) +afun_compare(x, n) +x = rep(rnorm(1), 1e3+1); n = sample(51, length(x), TRUE) +afun_compare(x, n) +num = 7400.0 #### no NA x = rnorm(1e3); n = sample(50, length(x), TRUE) # x even, n even afun_compare(x, n) diff --git a/src/data.table.h b/src/data.table.h index 6fb311c72..775d492aa 100644 --- a/src/data.table.h +++ b/src/data.table.h @@ -258,7 +258,7 @@ void frolladaptivesumExact(const double *x, uint64_t nx, ans_t *ans, const int * void frolladaptivemaxExact(const double *x, uint64_t nx, ans_t *ans, const int *k, double fill, bool narm, int hasnf, bool verbose); //void frolladaptiveminFast(const double *x, uint64_t nx, ans_t *ans, const int *k, double fill, bool narm, int hasnf, bool verbose); // does not exists as of now void frolladaptiveminExact(const double *x, uint64_t nx, ans_t *ans, const int *k, double fill, bool narm, int hasnf, bool verbose); -void frolladaptiveprodFast(const double *x, uint64_t nx, ans_t *ans, const int *k, double fill, bool narm, int hasnf, bool verbose); +//void frolladaptiveprodFast(const double *x, uint64_t nx, ans_t *ans, const int *k, double fill, bool narm, int hasnf, bool verbose); // does not exists as of now void frolladaptiveprodExact(const double *x, uint64_t nx, ans_t *ans, const int *k, double fill, bool narm, int hasnf, bool verbose); //void frolladaptivemedianFast(const double *x, uint64_t nx, ans_t *ans, const int *k, double fill, bool narm, int hasnf, bool verbose); // does not exists as of now void frolladaptivemedianExact(const double *x, uint64_t nx, ans_t *ans, const int *k, double fill, bool narm, int hasnf, bool verbose); diff --git a/src/froll.c b/src/froll.c index 532a3db13..6c47be246 100644 --- a/src/froll.c +++ b/src/froll.c @@ -895,7 +895,11 @@ void frollminExact(const double *x, uint64_t nx, ans_t *ans, int k, double fill, #undef PROD_WINDOW_STEP_FRONT #define PROD_WINDOW_STEP_FRONT \ if (R_FINITE(x[i])) { \ - w *= x[i]; \ + if (x[i] == 0.0) { \ + zerc++; \ + } else { \ + w *= x[i]; \ + } \ } else if (ISNAN(x[i])) { \ nc++; \ } else if (x[i]==R_PosInf) { \ @@ -906,7 +910,11 @@ void frollminExact(const double *x, uint64_t nx, ans_t *ans, int k, double fill, #undef PROD_WINDOW_STEP_BACK #define PROD_WINDOW_STEP_BACK \ if (R_FINITE(x[i-k])) { \ - w /= x[i-k]; \ + if (x[i-k] == 0.0) { \ + zerc--; \ + } else { \ + w /= x[i-k]; \ + } \ } else if (ISNAN(x[i-k])) { \ nc--; \ } else if (x[i-k]==R_PosInf) { \ @@ -918,18 +926,34 @@ void frollminExact(const double *x, uint64_t nx, ans_t *ans, int k, double fill, #define PROD_WINDOW_STEP_VALUE \ if (nc == 0) { \ if (pinf == 0 && ninf == 0) { \ - ans->dbl_v[i] = (double) w; \ + if (zerc) { \ + ans->dbl_v[i] = 0.0; \ + } else { \ + ans->dbl_v[i] = (double) w; \ + } \ } else { \ - ans->dbl_v[i] = (ninf+(w<0))%2 ? R_NegInf : R_PosInf; \ + if (zerc) { \ + ans->dbl_v[i] = R_NaN; \ + } else { \ + ans->dbl_v[i] = (ninf+(w<0))%2 ? R_NegInf : R_PosInf; \ + } \ } \ } else if (nc == k) { \ ans->dbl_v[i] = narm ? 1.0 : NA_REAL; \ } else { \ if (narm) { \ if (pinf == 0 && ninf == 0) { \ - ans->dbl_v[i] = (double) w; \ + if (zerc) { \ + ans->dbl_v[i] = 0.0; \ + } else { \ + ans->dbl_v[i] = (double) w; \ + } \ } else { \ - ans->dbl_v[i] = (ninf+(w<0))%2 ? R_NegInf : R_PosInf; \ + if (zerc) { \ + ans->dbl_v[i] = R_NaN; \ + } else { \ + ans->dbl_v[i] = (ninf+(w<0))%2 ? R_NegInf : R_PosInf; \ + } \ } \ } else { \ ans->dbl_v[i] = NA_REAL; \ @@ -951,34 +975,59 @@ void frollprodFast(const double *x, uint64_t nx, ans_t *ans, int k, double fill, return; } long double w = 1.0; + int zerc = 0; bool truehasnf = hasnf>0; if (!truehasnf) { int i; for (i=0; idbl_v[i] = fill; } - w *= x[i]; - ans->dbl_v[i] = (double) w; + if (x[i] == 0.0) { + zerc++; + } else { + w *= x[i]; + } + if (zerc) { + ans->dbl_v[i] = 0.0; + } else { + ans->dbl_v[i] = (double) w; + } if (R_FINITE((double) w)) { for (uint64_t i=k; idbl_v[i] = (double) w; + if (x[i-k] == 0.0) { + zerc--; + } else { + w /= x[i-k]; + } + if (x[i] == 0.0) { + zerc++; + } else { + w *= x[i]; + } + if (zerc) { + ans->dbl_v[i] = 0.0; + } else { + ans->dbl_v[i] = (double) w; + } } if (!R_FINITE((double) w)) { if (hasnf==-1) ansSetMsg(ans, 2, "%s: has.nf=FALSE used but non-finite values are present in input, use default has.nf=NA to avoid this warning", __func__); if (verbose) ansSetMsg(ans, 0, "%s: non-finite values are present in input, re-running with extra care for NFs\n", __func__); - w = 1.0; truehasnf = true; + w = 1.0; zerc = 0; truehasnf = true; } } else { if (hasnf==-1) ansSetMsg(ans, 2, "%s: has.nf=FALSE used but non-finite values are present in input, use default has.nf=NA to avoid this warning", __func__); if (verbose) ansSetMsg(ans, 0, "%s: non-finite values are present in input, skip non-finite inaware attempt and run with extra care for NFs straighaway\n", __func__); - w = 1.0; truehasnf = true; + w = 1.0; zerc = 0; truehasnf = true; } } if (truehasnf) { diff --git a/src/frolladaptive.c b/src/frolladaptive.c index 7019835d4..6e5bf8edd 100644 --- a/src/frolladaptive.c +++ b/src/frolladaptive.c @@ -41,11 +41,11 @@ void frolladaptivefun(rollfun_t rfun, unsigned int algo, const double *x, uint64 frolladaptiveminExact(x, nx, ans, k, fill, narm, hasnf, verbose); break; case PROD : - if (algo==0) { - frolladaptiveprodFast(x, nx, ans, k, fill, narm, hasnf, verbose); - } else if (algo==1) { - frolladaptiveprodExact(x, nx, ans, k, fill, narm, hasnf, verbose); + if (algo==0 && verbose) { + //frolladaptiveprodFast(x, nx, ans, k, fill, narm, hasnf, verbose); // frolladaptiveprodFast does not exists as of now + snprintf(end(ans->message[0]), 500, _("%s: algo %u not implemented, fall back to %u\n"), __func__, algo, (unsigned int) 1); } + frolladaptiveprodExact(x, nx, ans, k, fill, narm, hasnf, verbose); break; case MEDIAN : if (algo==0 && verbose) { @@ -642,123 +642,13 @@ void frolladaptiveminExact(const double *x, uint64_t nx, ans_t *ans, const int * } } -#undef PROD_WINDOW_STEP_VALUE -#define PROD_WINDOW_STEP_VALUE \ - if (wn>0) { \ - if (narm) { \ - if (wpinf == 0 && wninf == 0) { \ - int thisk = k[i] - ((int) wn); \ - ans->dbl_v[i] = thisk==0 ? 1.0 : (double) ws; \ - } else { \ - ans->dbl_v[i] = (wninf+(ws<0))%2 ? R_NegInf : R_PosInf; \ - } \ - } else { \ - ans->dbl_v[i] = NA_REAL; \ - } \ - } else { \ - if (wpinf == 0 && wninf == 0) { \ - ans->dbl_v[i] = (double) ws; /* k[i]==0 produces 1.0 */ \ - } else { \ - ans->dbl_v[i] = (wninf+(ws<0))%2 ? R_NegInf : R_PosInf; \ - } \ - } - -/* fast rolling adaptive prod - fast - * same as adaptive mean fast +/* fast rolling adaptive prod - fast - now redirects to exact! + * it was same as adaptive mean fast in 5975ffdc9839d507f4cd9cbba960198e301cf604 + * it was fast but not fully accurate and considering non-adaptive partial=TRUE falls back to adaptive we would like it to be as accurate as non-adaptive + * this could be implemented by tracking logsum rather than prod, together with zero and negative, as as non-adaptive fast roll prod + * for the current moment it will redirect to algo=exact, follow git hash for previous implementation if needed */ -void frolladaptiveprodFast(const double *x, uint64_t nx, ans_t *ans, const int *k, double fill, bool narm, int hasnf, bool verbose) { - if (verbose) - snprintf(end(ans->message[0]), 500, _("%s: running for input length %"PRIu64", hasnf %d, narm %d\n"), "frolladaptiveprodFast", (uint64_t)nx, hasnf, (int) narm); - bool truehasnf = hasnf>0; - long double w = 1.0; - double *cs = malloc(sizeof(*cs) * nx); - if (!cs) { // # nocov start - ansSetMsg(ans, 3, "%s: Unable to allocate memory for cumprod", __func__); // raise error - return; - } // # nocov end - if (!truehasnf) { - for (uint64_t i=0; idbl_v[i] = cs[i]; - } else if (i+1 > k[i]) { - ans->dbl_v[i] = cs[i]/cs[i-k[i]]; // for k[i]==0 it turns into cs[i]/cs[i] = 1, as expected - } else { - ans->dbl_v[i] = fill; - } - } - } else { - if (hasnf==-1) - ansSetMsg(ans, 2, "%s: has.nf=FALSE used but non-finite values are present in input, use default has.nf=NA to avoid this warning", __func__); - if (verbose) - ansSetMsg(ans, 0, "%s: non-finite values are present in input, re-running with extra care for NFs\n", __func__); - w = 1.0; truehasnf = true; - } - } - if (truehasnf) { - uint64_t nc = 0, pinf = 0, ninf = 0; // running NA counter - uint64_t *cn = malloc(sizeof(*cn) * nx); // cumulative NA counter, used the same way as cumprod, same as uint64_t cn[nx] but no segfault - if (!cn) { // # nocov start - ansSetMsg(ans, 3, "%s: Unable to allocate memory for cum NA counter", __func__); // raise error - free(cs); - return; - } // # nocov end - uint64_t *cpinf = malloc(sizeof(*cpinf) * nx); - if (!cpinf) { // # nocov start - ansSetMsg(ans, 3, "%s: Unable to allocate memory for cum Inf counter", __func__); // raise error - free(cs); free(cn); - return; - } // # nocov end - uint64_t *cninf = malloc(sizeof(*cninf) * nx); - if (!cninf) { // # nocov start - ansSetMsg(ans, 3, "%s: Unable to allocate memory for cum -Inf counter", __func__); // raise error - free(cs); free(cn); free(cpinf); - return; - } // # nocov end - for (uint64_t i=0; idbl_v[i] = fill; - } else if (i+1 == k[i]) { // first full window - wn = cn[i]; - wpinf = cpinf[i]; - wninf = cninf[i]; - ws = cs[i]; - PROD_WINDOW_STEP_VALUE - } else { // all the remaining full windows - wn = cn[i] - cn[i-k[i]]; // NAs in current window - wpinf = cpinf[i] - cpinf[i-k[i]]; // Inf in current window - wninf = cninf[i] - cninf[i-k[i]]; // -Inf in current window - ws = cs[i] / cs[i-k[i]]; // cumprod in current window - PROD_WINDOW_STEP_VALUE - } - } - free(cninf); free(cpinf); free(cn); - } - free(cs); -} + /* fast rolling adaptive prod - exact * same as adaptive mean exact */