[Rcpp-devel] speed
Romain Francois
romain at r-enthusiasts.com
Thu Aug 5 11:42:17 CEST 2010
Le 04/08/10 21:56, Douglas Bates a écrit :
> 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]
Actually, what happens here is that the IntegerVector( int ) constructor
gets called which creates a vector of the requested size:
> fx <- cxxfunction( , 'IntegerVector x = 10 ; return x ;', plugin =
"Rcpp" )
> fx()
[1] 0 0 0 0 0 0 0 0 0 0
So Richard's code works because of two mistakes that sort of balance
each other.
Creating an IntegerVector of length one with an int can be done using
IntegerVector::create, as in:
> fx <- cxxfunction( , 'IntegerVector x = IntegerVector::create( 10 ) ;
return x ;', plugin = "Rcpp" )
> fx()
[1] 10
Romain
> 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 mailing list
> Rcpp-devel at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>
--
Romain Francois
Professional R Enthusiast
+33(0) 6 28 91 30 30
http://romainfrancois.blog.free.fr
|- http://bit.ly/aAyra4 : highlight 0.2-2
|- http://bit.ly/94EBKx : inline 0.3.6
`- http://bit.ly/aryfrk : useR! 2010
More information about the Rcpp-devel
mailing list