[Rcpp-devel] speed

[Romain Francois]

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