R - vectorised conditional replace

Question

Hi I'm trying manipulate a list of numbers and I would like to do so without a for loop, using fast native operation in R. The pseudocode for the manipulation is :

By default the starting total is 100 (for every block within zeros)

From the first zero to next zero, the moment the cumulative total falls by more than 2% replace all subsequent numbers with zero.

Do this far all blocks of numbers within zeros

The cumulative sums resets to 100 every time

For example if following were my data :

d <- c(0,0,0,1,3,4,5,-1,2,3,-5,8,0,0,-2,-3,3,5,0,0,0,-1,-1,-1,-1);

Results would be :

0 0 0 1 3 4 5 -1 2 3 -5 0 0 0 -2 -3 0 0 0 0 0 -1 -1 -1 0

Currently I have an implementation with a for loop, but since my vector is really long, the performance is terrible.

Thanks in advance.

Here is a running sample code :

d <- c(0,0,0,1,3,4,5,-1,2,3,-5,8,0,0,-2,-3,3,5,0,0,0,-1,-1,-1,-1);
ans <- d;
running_total <- 100;
count <- 1;
max <- 100;
toggle <- FALSE;
processing <- FALSE;

for(i in d){
  if( i != 0 ){  
       processing <- TRUE; 
       if(toggle == TRUE){
          ans[count] = 0;  
       }
       else{
         running_total = running_total + i;

          if( running_total > max ){ max = running_total;}
          else if ( 0.98*max > running_total){
              toggle <- TRUE;  
          }
      }
   }

   if( i == 0 && processing == TRUE )
   { 
       running_total = 100; 
       max = 100;
       toggle <- FALSE;
   }
   count <- count + 1;
}
cat(ans)

Answer 1

I am not sure how to translate your loop into vectorized operations. However, there are two fairly easy options for large performance improvements. The first is to simply put your loop into an R function, and use the compiler package to precompile it. The second slightly more complicated option is to translate your R loop into a c++ loop and use the Rcpp package to link it to an R function. Then you call an R function that passes it to c++ code which is fast. I show both these options and timings. I do want to gratefully acknowledge the help of Alexandre Bujard from the Rcpp listserv, who helped me with a pointer issue I did not understand.

First, here is your R loop as a function, foo.r.

## Your R loop as a function
foo.r <- function(d) {
  ans <- d
  running_total <- 100
  count <- 1
  max <- 100
  toggle <- FALSE
  processing <- FALSE

  for(i in d){
    if(i != 0 ){
      processing <- TRUE
      if(toggle == TRUE){
        ans[count] <- 0
      } else {
        running_total = running_total + i;
        if (running_total > max) {
          max <- running_total
        } else if (0.98*max > running_total) {
          toggle <- TRUE
        }
      }
    }
    if(i == 0 && processing == TRUE) {
      running_total <- 100
      max <- 100
      toggle <- FALSE
    }
    count <- count + 1
  }
  return(ans)
}

Now we can load the compiler package and compile the function and call it foo.rcomp.

## load compiler package and compile your R loop
require(compiler)
foo.rcomp <- cmpfun(foo.r)

That is all it takes for the compilation route. It is all R and obviously very easy. Now for the c++ approach, we use the Rcpp package as well as the inline package which allows us to "inline" the c++ code. That is, we do not have to make a source file and compile it, we just include it in the R code and the compilation is handled for us.

## load Rcpp package and inline for ease of linking
require(Rcpp)
require(inline)

## Rcpp version
src <- '
  const NumericVector xx(x);
  int n = xx.size();
  NumericVector res = clone(xx);
  int toggle = 0;
  int processing = 0;
  int tot = 100;
  int max = 100;

  typedef NumericVector::iterator vec_iterator;
  vec_iterator ixx = xx.begin();
  vec_iterator ires = res.begin();
  for (int i = 0; i < n; i++) {
    if (ixx[i] != 0) {
      processing = 1;
      if (toggle == 1) {
        ires[i] = 0;
      } else {
        tot += ixx[i];
        if (tot > max) {
          max = tot;
        } else if (.98 * max > tot) {
            toggle = 1;
          }
      }
    }

   if (ixx[i] == 0 && processing == 1) {
     tot = 100;
     max = 100;
     toggle = 0;
   }
  }
  return res;
'

foo.rcpp <- cxxfunction(signature(x = "numeric"), src, plugin = "Rcpp")

Now we can test that we get the expected results:

## demonstrate equivalence
d <- c(0,0,0,1,3,4,5,-1,2,3,-5,8,0,0,-2,-3,3,5,0,0,0,-1,-1,-1,-1)
all.equal(foo.r(d), foo.rcpp(d))

Finally, create a much larger version of d by repeating it 10e4 times. Then we can run the three different functions, pure R code, compiled R code, and R function linked to c++ code.

## make larger vector to test performance
dbig <- rep(d, 10^5)

system.time(res.r <- foo.r(dbig))
system.time(res.rcomp <- foo.rcomp(dbig))
system.time(res.rcpp <- foo.rcpp(dbig))

Which on my system, gives:

> system.time(res.r <- foo.r(dbig))
   user  system elapsed 
  12.55    0.02   12.61 
> system.time(res.rcomp <- foo.rcomp(dbig))
   user  system elapsed 
   2.17    0.01    2.19 
> system.time(res.rcpp <- foo.rcpp(dbig))
   user  system elapsed 
   0.01    0.00    0.02 

The compiled R code takes about 1/6 the time the uncompiled R code taking only 2 seconds to operate on the vector of 2.5 million. The c++ code is orders of magnitude faster even then the compiled R code requiring just .02 seconds to complete. Aside from the initial setup, the syntax for the basic loop is nearly identical in R and c++ so you do not even lose clarity. I suspect that even if parts or all of your loop could be vectorized in R, you would be sore pressed to beat the performance of the R function linked to c++. Lastly, just for proof:

> all.equal(res.r, res.rcomp)
[1] TRUE
> all.equal(res.r, res.rcpp)
[1] TRUE

The different functions return the same results.