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

[Ron Wehrens]

We are in the process of updating the kohonen package and are 
considering to move to Rcpp.

As part of this process we are running timing comparisons between the 
three ways to call C/C++ code from R. 1) using the .C interface, 2) 
using the .Call interface, and 3) using Rcpp. Here we notice something 
strange that we can capture in the following minimalistic example:

Consider the task of computing the minimal distance (sum of squares) 
between the columns of two data matrices. For this, we have an 
implementation for .C, .Call, and Rcpp (see below for the code of .Call 
and Rcpp). If we compare the timings for these three, we notice that 
Rcpp is much slower than the .C and .Call version (see the graphics at 
https://github.com/jwkruisselbrink/rcpp-timings). When checking this 
further, we notice that if the ISNAN checks are removed from all three 
versions, Rcpp catches up with the .Call version.



The questions are therefor: How can this be? What is going on with ISNAN 
in Rcpp? What can we do about it? Or are we doing something very wrong here?



Best regards,



Johannes Kruisselbrink and Ron Wehrens



Source is available here: https://github.com/jwkruisselbrink/rcpp-timings



/****************** Start Rcpp ***********************/ #include 
<Rcpp.h> using namespace Rcpp;
// [[Rcpp::export]]

double RcppDist(NumericMatrix sData1,
   NumericMatrix sData2,
   int numObjects,
   int numVars,
   int numCodes) {
   int i, j, k;
   double dist, tmp, mindist,
     *data = REAL(sData1),
     *codes = REAL(sData2);
   mindist = 10000.0;

   for (i = 0; i < numObjects; i++) {
     for (j = 0; j < numCodes; j++) {
       dist = 0;
       for (k = 0; k < numVars; k++) {
         if (!ISNAN(data[i * numVars + k])) {
           tmp = data[i * numVars + k] - codes[j * numVars + k];
           dist += tmp * tmp;
         }
       }
       if (dist < mindist) {
         mindist = dist;
       }
     }
   }

   return mindist;
}

/****************** End Rcpp ***********************/



/****************** Start Call ***********************/


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

SEXP CallDist(
   SEXP sData,
   SEXP sCodes,
   SEXP sNumObjects,
   SEXP sNumVars,
   SEXP sNumCodes
) {
   int
     i, j, k,
     numObjects = *INTEGER(sNumObjects),
     numCodes = *INTEGER(sNumCodes),
     numVars = *INTEGER(sNumVars);
   double
     dist, tmp, *mindist,
     *data = REAL(sData),
     *codes = REAL(sCodes);
   SEXP output = PROTECT(allocVector(REALSXP, 1));
   mindist = REAL(output);
   *mindist = 10000.0;

   for (i = 0; i < numObjects; i++) {
     for (j = 0; j < numCodes; j++) {
       dist = 0;
       for (k = 0; k < numVars; k++) {
         if (!ISNAN(data[i * numVars + k])) {
           tmp = data[i * numVars + k] - codes[j * numVars + k];
           dist += tmp * tmp;
         }
       }
       if (dist < *mindist) {
         *mindist = dist;
       }
     }
   }

   UNPROTECT(1);
   return output;
}

/****************** End Call ***********************/