[Rcpp-devel] speed

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

_______________________________________________
Rcpp-devel mailing list
Rcpp-devel at lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel

[Rcpp-devel] speed

Thread (14 messages)