[Rcpp-devel] speed

[Dirk Eddelbuettel]

Richard,

Great question. 

One obvious consideration for improvement is to not call dbinom and dpois
from R, but rather using the C entry points provided by Rmath.h. With that, I
get these performances:

> require(rbenchmark)
Loading required package: rbenchmark
> benchmark(fx(0:100),
+           fx2(0:100),
+           fxR(0:100))
        test replications elapsed relative user.self 
1  fx(0:100)          100  89.426 3.811850     89.42 
2 fx2(0:100)          100  23.460 1.000000     23.46 
3 fxR(0:100)          100  41.699 1.777451     41.69 
> # deleted empty sys.self, user,child, sys.child columns here

The bad news is that I don't get the same results:

> all.equal(fx(0:10), fxR(0:10))
[1] TRUE
> all.equal(fx2(0:10), fxR(0:10))
[1] "Mean relative difference: 0.2535309"
> 

But as Doug suggested, there may be a need for reworking things and checks
anyway. Or maybe I just introduced a bug -- dunno. I merely meant to help on
pointing out that Rmath.h is there too, and I did this as a quick and dirty
wrap around the 'atomistic' Rmath functions. Maybe someone would want to
contribute vectorised versions of these ?  Patches welcome, as they say...

Dirk


inc2 <- '
#include <Rmath.h>

NumericVector wrap_dbinom(IntegerVector c, int k, double f) {
  NumericVector x(c.size());
  for (int i=0; i<c.size(); i++) x[i] = dbinom(c[i], k, f, false);
  return x;
}

NumericVector wrap_dpois(IntegerVector c, int l) {
  NumericVector x(c.size());
  for (int i=0; i<c.size(); i++) x[i] = dpois(c[i], l, false);
  return x;
}
'

fx2 <- cxxfunction(signature(Kr="integer"), '
    IntegerVector K(Kr);
    NumericMatrix bpsum(K.size(), K.size());
    for(int i=0; i<K.size(); i++) {
        for(int j=0; j<K.size(); j++) {
            IntegerVector Ki = K[i];
            IntegerVector cmin = seq_len(Ki.size()+1);
            IntegerVector cmin0 = cmin-1;
            NumericVector bin = wrap_dbinom(cmin0, K[j], 0.5);
            NumericVector pois = wrap_dpois(K[j]-cmin0, 1.5);
            NumericVector bp = bin * pois;
            bpsum(i, j) = std::accumulate(bp.begin(), bp.end(), 0.0);
            }
        }
    return bpsum;
    ', plugin="Rcpp", include=inc2)




-- 
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com