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
| [email protected]
|
--
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
_______________________________________________
Rcpp-devel mailing list
[email protected]
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
