[Rcpp-devel] A question about Rcpp

[Romain Francois]

Everything Dirk said is fine. I just wanted to clarify why the initial 
example does not work.

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 first argument is of type "int", that is not good. every .Call 
function take between zero and 65 SEXP arguments, that's it.

What happens is that the SEXP (those are pointers) you pass in from R is 
converted to an int, as conversions from pointers to int is possible, 
and so you get unpredictable results.



In the example where you use :

NumericVector x (n,norm_rand);

You also need to deal with the RNG state. You can do this as Dirk did, 
with GetRNGstate(); and PutRNGstate(); or you can use RNGScope, 
something like this:

RNGScope scope ;
NumericVector x (n,norm_rand);


Romain

Le 07/11/10 15:52, Dirk Eddelbuettel a écrit :
>
> 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
> |
>


-- 
Romain Francois
Professional R Enthusiast
+33(0) 6 28 91 30 30
http://romainfrancois.blog.free.fr
|- http://bit.ly/czHPM7 : Rcpp Google tech talk on youtube
|- http://bit.ly/9P0eF9 : Google slides
`- http://bit.ly/cVPjPe : Chicago R Meetup slides