Skip to content

[Rcpp-devel] Rcpp ISNAN slower than C ISNAN?

2 messages · Johannes Kruisselbrink, Romain Francois

Johannes Kruisselbrink
# 
Here is a further reduced example (see below). Now it is a function to 
count NaNs in a vector, and it shows the same behaviour. Code is also 
available at 
https://github.com/jwkruisselbrink/rcpp-timings/tree/master/minimal.

Regarding your question:
| Why not drop data and codes and use  sData1(i,k) - sData2(j,k)

Well, I wanted to rule out that Rcpp sugar was causing the slowdown. If 
it isn't we'll put it back in after having a fix for this issue.

/=============== call.c ===============/

#include <R.h>
#include <Rinternals.h>

SEXP CountNans(SEXP sX, SEXP sLength) {
   int i, n = *INTEGER(sLength);
   int *nanCount;
   double *x = REAL(sX);

   SEXP output = PROTECT(allocVector(INTSXP, 1));
   nanCount = INTEGER(output);
   *nanCount = 0;
   for (i = 0; i < n; i++) {
     if (ISNAN(x[i])) {
       *nanCount += 1;
     }
   }
   UNPROTECT(1);
   return output;
}

/=============== rcpp.cpp ===============/

#include <R.h>
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
int CountNans(NumericVector x) {
   int n = x.length();
   int nanCount = 0;
   for (int i = 0; i < n; i++) {
     if (ISNAN(x[i])) {
       nanCount++;
     }
   }
   return nanCount;
}

/=============== R code ===============/

library(Rcpp)
library(microbenchmark)
library(ggplot2)

sourceCpp('rcpp.cpp')

if (is.loaded(paste("call", .Platform$dynlib.ext, sep=""))) {
   dyn.unload(paste("call", .Platform$dynlib.ext, sep=""))
}
system(paste("R CMD SHLIB call.c", sep=""))
dyn.load(paste("call", .Platform$dynlib.ext, sep=""))

n <- as.integer(100000)
x <- rnorm(n)
mb <- microbenchmark(
   rcpp <- CountNans(x),
   call <- .Call("CountNans", x, n),
   times = 10000
)

autoplot(mb)
Romain Francois
# 
I?d be interesting to see what this more C++idomatic version would perform 

nanCount = std::count_if( x.begin(), x.end(), ISNAN ) ;

Because despite all the inlining efforts that has been put in the implementation of NumericVector, its operator[] might not perform as well as using [] on a double*. 

Romain