[Rd] .Call in R

Dirk Eddelbuettel edd at debian.org
Thu Nov 17 22:04:58 CET 2011

On 17 November 2011 at 09:09, Raymond wrote:
| Hi R developers,
|
|     I am new to this forum and hope someone can help me with .Call in R.
| Greatly appreciate any help!
|
|     Say, I have a vector called "vecA" of length 10000, I generate a vector
| called "vecR" with elements randomly generated from Uniform[0,1]. Both vecA
| and vecR are of double type. I want to replace elements vecA by elements in
| vecR only if sum of elements in vecR is greater than or equal to 5000.
| Otherwise, vecR remain unchanged. This is easy to do in R, which reads
|     vecA<-something;
|     vecR<-runif(10000);
|     if (sum(vecR)>=5000)){
|        vecA<-vecR;
|     }
|
|
|     Now my question is, if I am going to do the same thing in R using .Call.
| How can I achieve it in a more efficient way (i.e. less computation time
| compared with pure R code above.).  My c code (called "change_vecA.c") using
| .Call is like this:

Here is my take on it, using about the same number of commands in C++ thanks
to Rcpp and its vectorised sum() and runif() commands (which mimick the R
commands):

R> library(inline)
R> library(Rcpp)
R>
R> set.seed(42)         # fix RNG seed
R> vecA <- rt(10000, 6) # 'something' in vecA: t-dist with 6 df
R>
R> fun <- cxxfunction(signature(va="numeric"), # pass in a vector
+                    plugin="Rcpp",      # use Rcpp, and code below
+                    body='
+
+    Rcpp::NumericVector vA(va);
+    Rcpp::RNGScope tmp;                     // make sure RNG is set up
+    Rcpp::NumericVector vR = runif(10000);  // 10k of a U(0,1)
+
+    if (sum(vR) >= 5000) {                  // sum is an Rcpp sugar op.
+       vA = vR;                             // swap vR into vA
+    }
+    return(vA);                             // return vA
+ ')
R>
R> sum( fun( vecA ) )
 5033
R> sum( fun( vecA ) )
 5015
R> sum( fun( vecA ) )
 66
R> sum( fun( vecA ) )
 66
R> sum( fun( vecA ) )
 5015
R> sum( fun( vecA ) )
 5024
R> sum( fun( vecA ) )
 5020
R>
R> sum(vecA)
 66
R>

You can learn about Rcpp from the vignettes in the package, at my page at
http://dirk.eddelbuettel.com/code/rcpp.html as well different posts on my
blog, and of course the rcpp-devel mailing list.  The example above uses
cxxfunction() from the wonderful inline package you may find useful too as it