Skip to content

[Rcpp-devel] RNGScope, importing function, and if-else statement

10 messages · Marie Auger-Methe, Davor Cubranic, Dirk Eddelbuettel

#
Hi,

I want to use a random generator that is not available in the Rcpp: rvm 
from R's CircStats package and I want the number generated to be the 
same as the equivalent function in R. Importing the function within a 
Rcpp inline function works fine except if you have an if-else statement. 
The strange thing is that it works as long at the if-statement is TRUE, 
and then it doens't work anymore and some numbers are repeated but not 
in the same order.

Here is an example:

suppressMessages(library(CircStats))
suppressMessages(library(inline))


ft <- '
   Rcpp::NumericVector xc(clone(x));
   Rcpp::Function vmc(vm);
   int n = xc.size();
   NumericVector y(n);

   RNGScope Scope;

   for(int j=0; j<n; j++){
     if(xc[j] == 0){
       y[j] = as<double>(runif(1,0,2*PI));
     }else{
       y[j] = as<double>(vmc(1,0,1));
     }
   }

   return y;'

ftt <- cxxfunction(signature( x ="numeric", vm = "function"),
             body = ft, plugin = "Rcpp")

ftr <- function(x){
     n <- length(x)
     y <- numeric(n)
     for(i in 1:n){
           if(x[i] == 0){
             y[i] = runif(1,0,2*pi)
           }else{
             y[i] = rvm(1, 0, 1)
           }
     }
     return(y)
}

x <- c(0,0,0,1,1,1,0,0,0)

set.seed(9)
r1 <- ftt(x,rvm)
set.seed(9)
r2 <- ftr(x)
cbind(x,r1,r2)

Any idea what my be causing the difference between R and Rcpp? Many 
thanks in advance!

Marie
#
Hi Marie,
On 15 April 2012 at 19:01, Marie Auger-Methe wrote:
| Hi,
| 
| I want to use a random generator that is not available in the Rcpp: rvm 
| from R's CircStats package and I want the number generated to be the 
| same as the equivalent function in R. Importing the function within a 
| Rcpp inline function works fine except if you have an if-else statement. 
| The strange thing is that it works as long at the if-statement is TRUE, 
| and then it doens't work anymore and some numbers are repeated but not 
| in the same order.
| 
| Here is an example:
| 
| suppressMessages(library(CircStats))
| suppressMessages(library(inline))
| 
| 
| ft <- '
|    Rcpp::NumericVector xc(clone(x));
|    Rcpp::Function vmc(vm);
|    int n = xc.size();
|    NumericVector y(n);
| 
|    RNGScope Scope;
| 
|    for(int j=0; j<n; j++){
|      if(xc[j] == 0){
|        y[j] = as<double>(runif(1,0,2*PI));
|      }else{
|        y[j] = as<double>(vmc(1,0,1));
|      }
|    }
| 
|    return y;'
| 
| ftt <- cxxfunction(signature( x ="numeric", vm = "function"),
|              body = ft, plugin = "Rcpp")
| 
| ftr <- function(x){
|      n <- length(x)
|      y <- numeric(n)
|      for(i in 1:n){
|            if(x[i] == 0){
|              y[i] = runif(1,0,2*pi)
|            }else{
|              y[i] = rvm(1, 0, 1)
|            }
|      }
|      return(y)
| }
| 
| x <- c(0,0,0,1,1,1,0,0,0)
| 
| set.seed(9)
| r1 <- ftt(x,rvm)
| set.seed(9)
| r2 <- ftr(x)
| cbind(x,r1,r2)
| 
| Any idea what my be causing the difference between R and Rcpp? Many 
| thanks in advance!

I won't have a chance to disentangle this in detail, but could you possibly
just "reimplement" what rvm does using the Rcpp sugar functions?  

That way you would be guaranteed the "nice" behaviour of properly replicating
RNG draws from R in C++.

Dirk

| 
| Marie
| 
| _______________________________________________
| 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
#
At Sun, 15 Apr 2012 19:01:30 +0100, Marie Auger-Methe wrote:
In this example, it looks to me that the key is that as soon as you
start using 'vmc', the numbers get out of order. If that's the case,
then perhaps 'vmc' and 'runif' somehow both use the random seed, but
RNGScope keeps track only of 'runif's?

Davor
#
Hi Dirk,

Thank you for the prompt reply!
By reimplement rvm using Rcpp sugar functions, do you mean re-write the 
rvm function using sugar functions?

Marie
On 15/04/2012 9:39 PM, Dirk Eddelbuettel wrote:
#
Yes, there is something weird about the RNG. If you replace vmc by a 
function known to Rcpp such as rnorm this problem doesn't occur. If you 
create a function that uses only the imported vmc the problem doesn't 
occur either.  The problem resides in calling two functions, of which 
one is imported, since the same problem occur without the if-else statement:

suppressMessages(library(CircStats))
suppressMessages(library(inline))

ft <- '
   Rcpp::Function vmc(vm);
   Rcpp::NumericVector y(2);

   RNGScope Scope;

   y[0] = as<double>(runif(1,0,2*PI));
   y[1] = as<double>(vmc(1,0,1));

return y;'

ftt <- cxxfunction(signature( vm = "function"),
                    body = ft, plugin = "Rcpp")

ftr <- function(){
   y <- numeric(2)
     y[1] = runif(1,0,2*pi)
     y[2] = rvm(1, 0, 1)
   return(y)
}


set.seed(11)
r1 <- ftt(rvm)
set.seed(11)
r2 <- ftr()
cbind(r1,r2)

Marie
On 16/04/2012 6:11 AM, Davor Cubranic wrote:
#
Marie,
On 16 April 2012 at 08:56, Marie Auger-Methe wrote:
| Hi Dirk,
| 
| Thank you for the prompt reply!
| By reimplement rvm using Rcpp sugar functions, do you mean re-write the 
| rvm function using sugar functions?

Yes. 

(But as I said I have not yet looked at the code for rvm at all...)

Dirk
#
Hi Dirk and list,

As suggested by Dirk I have re-written the function rvm in Rcpp (see 
below) and it now gives the same results as the rvm function in R. So my 
problem is fixed (although it would nice if you could just import an R 
RNG function in Rcpp... maybe a future development ;) ).
I'm still not sure what was the problem with my original method, my only 
clue is that the function rvm itself randomly generates multiple numbers 
using runif?

Thanks for all the help,

Marie

suppressMessages(library(CircStats))
suppressMessages(library(inline))
suppressMessages(library(sugaR))

rvmscr <- '
   // Inputs
   int nc = as<int>(n);
   double meanc = as<double>(mean);
   double kc = as<double>(k);

   // Assign memory
   Rcpp::NumericVector vm(nc);
   double U1;
   double U2;
   double U3;
   double z;
   double f;
   double c;

   double a = 1 + pow((1 + 4 * pow(kc,2)),0.5);
   double b = (a - pow(2*a,0.5))/(2*kc);
   double r = (1 + pow(b,2))/(2 * b);

   int obs = 0;

   RNGScope scope;

   while(obs < nc){
     U1 = as<double>(runif(1, 0, 1));
     z = cos(PI * U1);
     f = (1 + r * z)/(r + z);
     c = kc * (r - f);
     U2 = as<double>(runif(1, 0, 1));
     if(c * (2 - c) - U2 > 0){
       U3 = as<double>(runif(1, 0, 1));
         if(U3 < 0.5){
           vm[obs] = 2 * PI - acos(f) + meanc;
         }else{
           vm[obs] = acos(f) + meanc;
         }
       vm[obs] = fmod(vm[obs], 2*PI);
       obs += 1;
     }else{
       if(log(c/U2) + 1 - c >= 0){
         U3 = as<double>(runif(1, 0, 1));
         if(U3 < 0.5){
           vm[obs] = 2 * PI - acos(f) + meanc;
         }else{
           vm[obs] = acos(f) + meanc;
         }
         vm[obs] = fmod(vm[obs], 2*PI);
         obs += 1;
       }
     }
   }

   return vm;
'

rvmc <- cxxfunction(signature(n ="integer", mean = "numeric", k = 
"numeric"),
                     body = rvmscr, plugin = "Rcpp")

n <- 1000
mu <- 0
k <- 3
set.seed(89)
r1 <- rvmc(n,mu,k)
set.seed(89)
r2 <- rvm(n,mu,k)

max(r1 - r2)
cbind(r1[1:20],r2[1:20])
On 16/04/2012 10:29 AM, Dirk Eddelbuettel wrote:
#
On 16 April 2012 at 12:43, Marie Auger-Methe wrote:
| Hi Dirk and list,
| 
| As suggested by Dirk I have re-written the function rvm in Rcpp (see 
| below) and it now gives the same results as the rvm function in R. So my 

Super!

| problem is fixed (although it would nice if you could just import an R 
| RNG function in Rcpp... maybe a future development ;) ).
| I'm still not sure what was the problem with my original method, my only 
| clue is that the function rvm itself randomly generates multiple numbers 
| using runif?

I am still travelling and have not had a chance to look at this, but I
suspect that maybe something odd happen with PutRNGState / GetRNGState on the
side of rvm (and RNGScope deals with this for us).  There is no reason why it
shouldn't work, maybe we can work out a fix to be passed to CircStats. Else
we can always provide a suitable rvm in Rcpp too I suppose...

Dirk

 
| Thanks for all the help,
| 
| Marie
| 
| suppressMessages(library(CircStats))
| suppressMessages(library(inline))
| suppressMessages(library(sugaR))
| 
| rvmscr <- '
|    // Inputs
|    int nc = as<int>(n);
|    double meanc = as<double>(mean);
|    double kc = as<double>(k);
| 
|    // Assign memory
|    Rcpp::NumericVector vm(nc);
|    double U1;
|    double U2;
|    double U3;
|    double z;
|    double f;
|    double c;
| 
|    double a = 1 + pow((1 + 4 * pow(kc,2)),0.5);
|    double b = (a - pow(2*a,0.5))/(2*kc);
|    double r = (1 + pow(b,2))/(2 * b);
| 
|    int obs = 0;
| 
|    RNGScope scope;
| 
|    while(obs < nc){
|      U1 = as<double>(runif(1, 0, 1));
|      z = cos(PI * U1);
|      f = (1 + r * z)/(r + z);
|      c = kc * (r - f);
|      U2 = as<double>(runif(1, 0, 1));
|      if(c * (2 - c) - U2 > 0){
|        U3 = as<double>(runif(1, 0, 1));
|          if(U3 < 0.5){
|            vm[obs] = 2 * PI - acos(f) + meanc;
|          }else{
|            vm[obs] = acos(f) + meanc;
|          }
|        vm[obs] = fmod(vm[obs], 2*PI);
|        obs += 1;
|      }else{
|        if(log(c/U2) + 1 - c >= 0){
|          U3 = as<double>(runif(1, 0, 1));
|          if(U3 < 0.5){
|            vm[obs] = 2 * PI - acos(f) + meanc;
|          }else{
|            vm[obs] = acos(f) + meanc;
|          }
|          vm[obs] = fmod(vm[obs], 2*PI);
|          obs += 1;
|        }
|      }
|    }
| 
|    return vm;
| '
| 
| rvmc <- cxxfunction(signature(n ="integer", mean = "numeric", k = 
| "numeric"),
|                      body = rvmscr, plugin = "Rcpp")
| 
| n <- 1000
| mu <- 0
| k <- 3
| set.seed(89)
| r1 <- rvmc(n,mu,k)
| set.seed(89)
| r2 <- rvm(n,mu,k)
| 
| max(r1 - r2)
| cbind(r1[1:20],r2[1:20])
|
| On 16/04/2012 10:29 AM, Dirk Eddelbuettel wrote:
| > Marie,
| >
| > On 16 April 2012 at 08:56, Marie Auger-Methe wrote:
| > | Hi Dirk,
| > |
| > | Thank you for the prompt reply!
| > | By reimplement rvm using Rcpp sugar functions, do you mean re-write the
| > | rvm function using sugar functions?
| >
| > Yes.
| >
| > (But as I said I have not yet looked at the code for rvm at all...)
| >
| > Dirk
| >
|
2 days later
#
Hi list,

sorry for bombarding the list with questions.
I'm trying to use Rcpp via the inline package and to use parallel 
computing at the same time. But I get this error:

first error: NULL value passed as symbol address
I've read a few posts (not related to parallel computing) discussing this error but I couldn't link it back to my specific problem.

Here is an example:

library(inline)
library(parallel)

# A silly Rcpp function
sillyscr<- '
   int x = as<int>(i);
   double y = 5.6;
   NumericVector j(1);
   j[0] =  x + y;

   return j;
'
silly<- cxxfunction(signature(i = "int"), body = sillyscr, plugin = "Rcpp")
silly(1) # Works!

# Equivalent function in R to make sure the problem is not the parallel computing
silly_notcpp<- function(i){
   j<- i +5.6
   return(j)
}
silly_notcpp(1) # Works !


M<- detectCores()
cl<- makeCluster(M)
clusterExport(cl, varlist=list("silly","silly_notcpp"))
res1<- parSapply(cl, 1:10,silly_notcpp) # Works!
res1
res2<- parSapply(cl, 1:10,silly) # Does not work !
res2
stopCluster(cl)

Many thanks!

Marie

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20120419/4e07e85a/attachment.html>
#
(Resending, forgot to CC list on first reply, Sorry --Dirk)

Marie,
On 19 April 2012 at 16:11, Marie Auger-Methe wrote:
| Hi list,
| 
| sorry for bombarding the list with questions.
| I'm trying to use Rcpp via the inline package and to use parallel computing at
| the same time. But I get this error:
| 
| first error: NULL value passed as symbol address
| I've read a few posts (not related to parallel computing) discussing this error but I couldn't link it back to my specific problem.
| 
| Here is an example:
| 
| library(inline)
| library(parallel)
| 
| # A silly Rcpp function
| sillyscr <- '
|   int x = as<int>(i);
|   double y = 5.6;
|   NumericVector j(1);
|   j[0] =  x + y;
| 
|   return j;
| '
| silly <- cxxfunction(signature(i = "int"), body = sillyscr, plugin = "Rcpp")
| silly(1) # Works!
| 
| # Equivalent function in R to make sure the problem is not the parallel computing
| silly_notcpp <- function(i){
|   j <- i +5.6
|   return(j)
| }
| silly_notcpp(1) # Works !
| 
| 
| M <- detectCores()
| cl <- makeCluster(M)
| clusterExport(cl, varlist=list("silly","silly_notcpp"))
| res1 <- parSapply(cl, 1:10,silly_notcpp) # Works!
| res1
| res2 <- parSapply(cl, 1:10,silly) # Does not work !
| res2
| stopCluster(cl)
| 
| Many thanks!

inline and cxxfunction create what looks to you like a normal function in the
current session, here this function is called 'silly'.  It really does a
.Call() of some object code compiled in the session temp directory.

While you send 'silly' to the nodes, I suspect that what gets called does
not. Recall that R really is single threaded, and that parallel (just like
snow) creates new sessions for the nodes, and new session temp directories.
Those have 'silly', but not the code it calls (eg the object code etc).

For 'real applications', we generally suggest a package.
Rcpp.package.skeleton() can take silly as an argument to aid you in building
package.  

Else, you may have to have the cxxfunction() call, and hence the 'sillyscr'
text variable, on each node.

The r-sig-hpc lists discusses some of these issues (about how these simple
parSapply calls can fail for lack of local resources) every now and then.

Dirk

-- 
R/Finance 2012 Conference on May 11 and 12, 2012 at UIC in Chicago, IL
See agenda, registration details and more at http://www.RinFinance.com