[Rcpp-devel] speed

[Douglas Bates]

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))
>> _______________________________________________
>> 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
>