[Rcpp-devel] First foray into Rcpp and comparison of speed with R

[Dirk Eddelbuettel]

Hi Gregor,

And welcome!

On 1 September 2011 at 10:06, Gregor GORJANC wrote:
| I wrote a simple application in R. Application is not very complex, but since I
| have to deal with large data I have devoted a lot of attention to vectorize R
| computations as much as possible. This lead me to a very decent performance,
| but my computation involves a large outer loop, due to the recursive nature of
| a problem. I looked into Rcpp and I was amazed that I was able to come up with
| its use in my application even though I do not had *any* practical experience
| with C++ before!!! I managed to half the execution speed of my application,
| which is great. I kind of hoped for more, but given that I vectorized much of R
| code and that computations are simple additions and multiplications I wonder if
| I could get much more gains with Rcpp foray.

Hard to say in general terms.  Rcpp sugar can help, Romain has hidden some
loop unrolling inside it and that helps in the benchmark example in the
Rcpp-introduction vignette (and JSS article) but it doesn't imply all sugar
code will be faster.

| Does anyone have time/will to quickly look at my code bellow to see if I am
| missing something very obvious. I am basically doing one walk from top to
| bottom of a rather large matrix (most of it is dense). In R I was able to
| vectorize all within line computations, but in C++ I have to use explicit
| looping, right?

Look at the Rcpp-sugar vignette. Some things can be more compact, but the
devil is as always in the details.  Starting with explicit loops is a good
idea. 
 
| Thx!
| 
| Gregor Gorjanc
| 
| SEXP fx(SEXP c1_, SEXP c2_, SEXP nI_, SEXP nP_, SEXP nT_, SEXP y_, SEXP P_,
| SEXP Px_) {
|   using namespace Rcpp ;
|   int i, j, k, t, p;
|   double c1 = Rcpp::as<double>(c1_);  
|   double c2 = Rcpp::as<double>(c2_); 
|   int nI = Rcpp::as<int>(nI_); 
|   int nP = Rcpp::as<int>(nP_);
|   int nT = Rcpp::as<int>(nT_);
|   Rcpp::NumericMatrix ped(y_);
|   Rcpp::IntegerVector P(P_);  
|   Rcpp::IntegerVector Px(Px_);
|   Rcpp::NumericMatrix pa(nI+1, nT);
|   Rcpp::NumericMatrix  w(nI+1, nT);
|   Rcpp::NumericMatrix xa(nI+1, nP*nT);
|   
|   // --- Compute ---
|       
|   for(i = 1; i < nI+1; i++) {

Careful about i=1:  C/C++ indices run from 0 to n-1.

|     for(t = 0; t < nT; t++) {
|       pa(i, t) = c1 * ped(ped(i, 1), t+3) +
|                  c2 * ped(ped(i, 2), t+3);
|     
|       w(i, t) = ped(i, t+3) - pa(i, t);

The recursive nature may make it tricky to vectorise this.

|     
|       for(p = 0; p < nP; p++) {
|         j = Px[t] + p;
|         xa(i, j) = c1 * xa(ped(i, 1), j) +
|                    c2 * xa(ped(i, 2), j);   
|       }
|       k = Px[t] + P[i];
|       xa(i, k) += w(i, t);
|     }
|   }
|   
|   // --- Return ---
|       
|   Rcpp::List ret;
|   ret["pa"] = pa;
|   ret["w"]  = w;
|   ret["xa"] = xa;
|   return(ret);

We usually do 

  Rcpp::List ret = Rcpp::List::create(Rcpp::Named("pa", pa),
                                      Rcpp::Named("w", w);
                                      Rcpp::Named("xa", xa));
  return ret;

to avoid re-allocation.  But that is minuscule, and your code is arguably
as easy or easier to read :)

Hope this helps,  Dirk

-- 
Two new Rcpp master classes for R and C++ integration scheduled for 
New York (Sep 24) and San Francisco (Oct 8), more details are at
http://dirk.eddelbuettel.com/blog/2011/08/04#rcpp_classes_2011-09_and_2011-10
http://www.revolutionanalytics.com/products/training/public/rcpp-master-class.php