[Rcpp-devel] speed

On Wed, Aug 4, 2010 at 2:27 PM, Sanjog Misra<misra at simon.rochester.edu>  wrote:

Richard,

I did a quick check and the Rcpp/C++ version is at least twice as fast as
the pure R code.

I get a similar result.

But a more important question is whether your C++ code is producing
the same values are the R function.  In R your sequence is 0:K[i] but
in C++ your code is

             IntegerVector Ki = K[i];
             IntegerVector cmin = seq_len(Ki.size()+1);
             IntegerVector cmin0 = cmin-1;

Because K is an Integervector, K[i] will be an int and Ki will be an
IntegerVector of length 1.  That is Ki.size() will always be 1 and
cmin0 will be equivalent to the R expression 0:1, not 0:K[i]

There are several other places that the code could be tightened.  For
example, you can avoid the call to the R functions because Rf_dbinom
and Rf_dpois are part of the R API.  Furthermore, the inner loop can
work with scalars instead of vectors if you do the accumulation in a
C++ loop.  Dirk and Romain have done a magnificent job making
different types of R vectors available within C++ but, even so,
operations on C++ native types will always be faster.

But before pursuing that it would be best to straighten out what cmin
and cmin0 should be.  I think they are incorrectly defined.  (Also,
they are loop invariants and their definition, however it is
reconciled, can be moved outside the loop on j.  A smart compiler may
do this but you shouldn't count on that.)

Best
Sanjog


On 8/4/10 3:21 PM, "Richard Chandler"<richard.chandlers at gmail.com>  wrote:

Hi,

I wrote my first Rcpp/C++ program in hopes of speeding up an R
function, but alas it is slower. I think the C++ program is slower
because I have made heavy use of dbinom and dpois from R. Is there a
way to do this without calling back to R? Are there any other obvious
ways to speed up the C++ program? I realize that I can vectorize the R
function and avoid zero probabilities, but I have showed it this way
for simplicity.

Thanks,
Richard



fx<- cxxfunction(signature(Kr="integer"), '
     Environment stats("package:stats");
     Function dbinom = stats["dbinom"];
     Function dpois = stats["dpois"];

     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 = dbinom(cmin0, K[j], 0.5);
             NumericVector pois = 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")



fxR<- function(K) {
     lk<- length(K)
     bpsum<- matrix(NA, lk, lk)
     for(i in 1:lk)
         for(j in 1:lk)
             bpsum[i, j]<- sum(dbinom(0:K[i], K[j], 0.5) *
dpois(K[j]-(0:K[i]), 1.5))
     return(bpsum)
     }


all.equal(fx(0:10), fxR(0:10))
system.time(fx(0:100))
system.time(fxR(0:100))