[Rcpp-devel] speed

Sanjog Misra misra at simon.rochester.edu
Wed Aug 4 21:27:46 CEST 2010


Richard,

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

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



More information about the Rcpp-devel mailing list