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 ***********************/
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
