[Rcpp-devel] A question about Rcpp

[Dirk Eddelbuettel]

Hi Chunhua and Siddharta,

I have taken the liberty of redirecting this to the rcpp-devel which is the
best forum for questions about Rcpp as it exposes them to a number of clever
readers.

On 7 November 2010 at 06:08, WU, Chunhua wrote:
| Hi, Dirk,
| 
| We are quite interested in your Rcpp and RcppArmadillo packages. We are using these packages to create our own package. We have a question regarding how to call the internal random number generators in R. After searching for a while we tried the following codes:
| 
| extern "C" SEXP test_normal(int n, SEXP mu, SEXP sd){
|   Function rnorm("rnorm");
|   NumericVector x = rnorm(n, ["mean"]=mu,_["sd"]=sd);
|   return x ;
| }
| 
| The above codes compiled successfully under both Linux (Ubuntu 10.10 64bit) and also Windows XP 32bit. It worked under Linux but did not work on Windows. On windows, it outputted a full screen of numbers and stating that still 43309713 entries omitted.
| 
| We also tried the code :
| 
| NumericVector x (n,norm_rand);
| 
| It did not compile under Linux, while compiled on Windows with same problem of output too many numbers.
| Could you help what is wrong with that?

It uses the "pre-sugar" form of calling the random number generator. You can
do a little better now.

One main problem (besides a few small syntactic ones) is that you did NOT set
and restore the random number generator state.  So a quick test I just wrote
would be

-----------------------------------------------------------------------------
library(inline)

src <- '
    int n = Rcpp::as<int>(ns);
    double m = Rcpp::as<double>(ms);
    double sd = Rcpp::as<double>(sds);
    GetRNGstate();
    Rcpp::NumericVector x = rnorm(n, m, sd);
    PutRNGstate();
    return x;'

fun <- cxxfunction(signature(ns="integer", ms="numeric", sds="numeric"),
                   body=src, plugin="Rcpp")
-----------------------------------------------------------------------------


which works just fine:


> src <- '
+     int n = Rcpp::as<int>(ns);
+     double m = Rcpp::as<double>(ms);
+     double sd = Rcpp::as<double>(sds);
+     GetRNGstate();
+     Rcpp::NumericVector x = rnorm(n, m, sd);
+     PutRNGstate();
+     return x;'
> 
> fun <- cxxfunction(signature(ns="integer", ms="numeric", sds="numeric"),
+                    body=src, plugin="Rcpp")
> set.seed(42); fun(3,1,1);
[1] 2.3710 0.4353 1.3631
> set.seed(42); rnorm(3,1,1);
[1] 2.3710 0.4353 1.3631
> 

and we get the same numbers rnorm gives us in R which is rather nice for
verification.

There are a few different 'patterns' for how you can call rnorm() et al from
C++ -- have a look at the existing example, the sugar vignette and the unit
tests.

 
| Also we read somewhere that the first method (define the function) is not
|efficient, is that true? 

Yes, there is a cost to the lookup from R.

Sugar as used here should be faster.

| Our ultimate goal is to use RcppArmadillo and call random number generators
| extensively. What would be an efficient way you would suggest? 

That's a good route. Armadillo and therefore RcppArmadillo have their own
RNGs.  Which ones are 'more efficient' I have not tested.  Going with the
ones from R had advantages for comparison and debugging.

Hope this helps, keep bringing questions to the list.

Cheers, Dirk
 
| Thanks very much!
| 
| Best,
| Chunhua
| 
| 
| __________________________
| Chunhua Wu
| Ph.D. Candidate in Marketing
| Olin Business School,
| Washington University in St. Louis
| chunhuawu at wustl.edu
| 

-- 
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com