From edd at debian.org  Thu Nov 17 22:04:58 2011
From: edd at debian.org (Dirk Eddelbuettel)
Date: Thu, 17 Nov 2011 15:04:58 -0600
Subject: [Rd] .Call in R
In-Reply-To: <1321549770463-4080721.post@n4.nabble.com>
References: <1321549770463-4080721.post@n4.nabble.com>
Message-ID: <20165.30458.587277.817354@max.nulle.part>


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 ) )
[1] 5033
R> sum( fun( vecA ) )
[1] 5015
R> sum( fun( vecA ) )
[1] 66
R> sum( fun( vecA ) )
[1] 66
R> sum( fun( vecA ) )
[1] 5015
R> sum( fun( vecA ) )
[1] 5024
R> sum( fun( vecA ) )
[1] 5020
R>
R> sum(vecA)
[1] 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
compiles, links and loads your C or C++ snippets.

Dirk

-- 
"Outside of a dog, a book is a man's best friend. Inside of a dog, it is too
dark to read." -- Groucho Marx



[Rd] .Call in R
Dirk Eddelbuettel