Skip to content

[Rcpp-devel] Rcpp bug with rgamma

6 messages · Colin Rundel, Dirk Eddelbuettel, Romain Francois

Original
Colin Rundel
# 
I am working on a project using RcppArmadillo and I've run into an issue with the rgamma function in Rcpp. When calling rgamma the function pegs R's cpu utilization and the process continues to churn forcing me to kill it. I've let things run for around 5 mins with no end in sight.

I can replicate the problem with the following code:

library(Rcpp)
library(inline)

foo = "return(rgamma(1,1,1));"
cxxfunction(signature(), foo, plugin = "Rcpp" )()


Other r* functions seem to work with out issue (ie rnorm, rbeta, etc).

I am using of Rcpp 0.86 on OSX with R 2.11.1. I have not taken the time to look at the Rcpp code to see what the issue might be, but hopefully this is something that is an easy fix.

-Colin
Dirk Eddelbuettel
#
On 2 October 2010 at 22:17, Colin Rundel wrote:
| I am working on a project using RcppArmadillo and I've run into an issue with the rgamma function in Rcpp. When calling rgamma the function pegs R's cpu utilization and the process continues to churn forcing me to kill it. I've let things run for around 5 mins with no end in sight.
| 
| I can replicate the problem with the following code:
| 
| library(Rcpp)
| library(inline)
| 
| foo = "return(rgamma(1,1,1));"
| cxxfunction(signature(), foo, plugin = "Rcpp" )()

I can confirm this on Ubuntu 10.04. 

The code for rgamma is pretty simple -- if all args are well we call R's
rgamma n time with argument a and scale:

	return NumericVector( n, ::Rf_rgamma, a, scale ) ;

Changing your test program to 

        foo <- "return(Rcpp::wrap(Rf_rgamma(1,1)));"

so that it calls the function directly also goes astray.  There must now be
something we define or set that sends the R engine over the edge here.

I'll try to take another peak.  Thanks for the bug report! 

Dirk
 

| Other r* functions seem to work with out issue (ie rnorm, rbeta, etc).
| 
| I am using of Rcpp 0.86 on OSX with R 2.11.1. I have not taken the time to look at the Rcpp code to see what the issue might be, but hopefully this is something that is an easy fix.
| 
| -Colin
| 
| 
| _______________________________________________
| Rcpp-devel mailing list
| Rcpp-devel at lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
Dirk Eddelbuettel
# 
A standalone program linked against libRmath (the R stats functions compiled
into their own library) works:

#define MATHLIB_STANDALONE 1
#include <Rmath.h>

int
main(int argc, char** argv)
{
    int i;
    set_seed(123, 456);
    for (i=0; i<10; i++) {
	printf(" %d rgamma(1,1) %f\n", i, rgamma(1,1));
    }
    return 0;
}

edd at max:~$ gcc -o /tmp/colin /tmp/colin.c -lRmath -lm
edd at max:~$ /tmp/colin 
 0 rgamma(1,1) 0.314001
 1 rgamma(1,1) 1.498724
 2 rgamma(1,1) 0.154940
 3 rgamma(1,1) 0.053388
 4 rgamma(1,1) 0.404130
 5 rgamma(1,1) 0.167460
 6 rgamma(1,1) 0.260743
 7 rgamma(1,1) 3.204146
 8 rgamma(1,1) 0.271189
 9 rgamma(1,1) 1.129695
edd at max:~$ 

So the ball in our court.  Rcpp must be defining something that upsets the
inner loops in rgamma.

Dirk
Dirk Eddelbuettel
# 
Colin,

The problem goes away as soon as you set the seeds to the RNGs:

   edd at max:/tmp$ cat colin.r
   library(Rcpp)
   library(inline)
   
   set.seed(42)  # this makes all the difference
   foo <- "return(Rcpp::wrap(rgamma(3,1,1)));"
   fun <- cxxfunction(signature(), foo, plugin = "Rcpp")
   
   print(fun())
   cat("Bye\n");
   edd at max:/tmp$ r colin.r
   Loading required package: methods
   [1] 1.93930 0.18042 0.53443
   Bye
   edd at max:/tmp$ 

I am not quite sure why that is needed for rgamma() when it isn't for all the
other distributions. We'll take a look -- but your code should now be fine.

Cheers, Dirk
Dirk Eddelbuettel
# 
Colin,

Just as one last comment, the issue at hand is described in Section 6.3 of
Writing R Extensions:

   6.3 Random number generation
   ============================
   
   The interface to R's internal random number generation routines is
   
        double unif_rand();
        double norm_rand();
        double exp_rand();
   
   giving one uniform, normal or exponential pseudo-random variate.
   However, before these are used, the user must call
   
        GetRNGstate();
   
   and after all the required variates have been generated, call
   
        PutRNGstate();
   
   These essentially read in (or create) `.Random.seed' and write it out
   after use.

   [...]

We had not made that sufficiently clear. What you encountered with rgamma was
an internal algorithm using a rejection sampleing scheme calling unif_rand()
and exp_rand() ... but creating an infinite loop at it would always end up
with the same draws.

I amended the 'Rcpp-sugar' vignette as well as the header files to point
users to this Section 6.3 to make sure GetRNGstate() and PutRNGstate() are
used.

Cheers,  Dirk
Romain Francois
# 
Le 03/10/10 23:29, Dirk Eddelbuettel a ?crit :
In addition, we have the RNGScope class, whose constructor calls 
GetRNGstate and destruvctor calls PutRNGstate, so that you can do :

fx <- cxxfunction( , '
	RNGScope scope ;
	NumericVector x = rgamma( 10, 1, 1 ) ;
	return x ;
', plugin = "Rcpp" )

fx()

Romain