New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DEoptim fails when fn calls another Rcpp-compiled function #14

Closed
kylebmetrum opened this Issue Jun 6, 2018 · 7 comments

Comments

Projects
None yet
3 participants
@kylebmetrum

kylebmetrum commented Jun 6, 2018

Using RcppDE::DEoptim to optimize a function that calls another Rcpp-compiled function. This seems to cause the optimization to fail (not an error, but it doesn't optimize properly).

A minimal example is shown below.

I saw this starting to happen when this was reported:
RcppCore/Rcpp#836

Up to that point, this sort of thing had been working for me. At the time I rolled back to 0.12.15 and that fixed it. Now running problems like this with 0.12.17 and still seeing the behavior.

Reporting here just because it's the only place where I can reproduce this.

library(Rcpp)
library(RcppDE)
## Now loading:
##   RcppDE: C++ Implementation of Differential Evolution Optimisation
##   Author: Dirk Eddelbuettel
## Based on:
##   DEoptim (version 2.0-7): Differential Evolution algorithm in R
##   Authors: David Ardia, Katharine Mullen, Brian Peterson and Joshua Ulrich
Rosenbrock <- function(x){
  x1 <- x[1]
  x2 <- x[2]
  100 * (x2 - x1 * x1)^2 + (1 - x1)^2
}

lower <- c(-10,-10)
upper <- -lower

set.seed(1234)
cont <- DEoptim.control(trace = 0)
fit0 <- DEoptim(Rosenbrock, lower, upper, control = cont)

fit0$optim
## $bestmem
## par1 par2 
##    1    1 
## 
## $bestval
## [1] 2.54782e-21
## 
## $nfeval
## [1] 10050
## 
## $iter
## [1] 200
sourceCpp(code='
          #include <Rcpp.h>
          
          // [[Rcpp::export]]
          int fibonacci(const int x) {
          if (x == 0) return(0);
          if (x == 1) return(1);
          return (fibonacci(x - 1)) + fibonacci(x - 2);
          }'
)

fibonacci(4)
## [1] 3
set.seed(1234)
Rosenbrock2 <- function(x){
  foo <- fibonacci(3)
  x1 <- x[1]
  x2 <- x[2]
  100 * (x2 - x1 * x1)^2 + (1 - x1)^2
}


set.seed(1234)
fit2 <- DEoptim(Rosenbrock2, lower, upper, control = cont)
fit2$optim
## $bestmem
##      par1      par2 
## -7.725932  2.445988 
## 
## $bestval
## [1] 327764.1
## 
## $nfeval
## [1] 10050
## 
## $iter
## [1] 200
library(DEoptim)
## Loading required package: parallel

## 
## DEoptim package
## Differential Evolution algorithm in R
## Authors: D. Ardia, K. Mullen, B. Peterson and J. Ulrich

## 
## Attaching package: 'DEoptim'

## The following objects are masked from 'package:RcppDE':
## 
##     DEoptim, DEoptim.control
set.seed(1234)
fit3 <- DEoptim(Rosenbrock2, lower, upper, control = DEoptim.control(trace = 0))
fit3$optim
## $bestmem
## par1 par2 
##    1    1 
## 
## $bestval
## [1] 3.010666e-16
## 
## $nfeval
## [1] 402
## 
## $iter
## [1] 200
sessionInfo()
## R version 3.5.0 (2018-04-23)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Sierra 10.12.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
## [1] DEoptim_2.2-4 RcppDE_0.1.5  Rcpp_0.12.17 
## 
## loaded via a namespace (and not attached):
##  [1] compiler_3.5.0  backports_1.1.2 magrittr_1.5    rprojroot_1.3-2
##  [5] tools_3.5.0     htmltools_0.3.6 yaml_2.1.19     stringi_1.2.2  
##  [9] rmarkdown_1.9   knitr_1.20      stringr_1.3.1   digest_0.6.15  
## [13] evaluate_0.10.1
@eddelbuettel

This comment has been minimized.

Owner

eddelbuettel commented Jun 6, 2018

That is confusing because the issue in #836 "only" had to do with whether / how the RNGstate was set and restored (which was something the Shiny team had noticed too).

I would not think it affects things here, but it might because each cppFunction() use sets / resets the RNG and what we did improces, I think, the general case. @kevinushey may have comments too.

For you, the need to add user-compiled functions may be an innocent bystand. I think this may just work if you turn off the RNG setting. I thought we had an option for this in sourceCpp() / cppFunction() but I do not see that now. @kevinushey didn't we add this?

In any event, it should be pretty easy for to cook up a cppFunction() clone that expands differently and does NOT pull these macros in. Then the RNG is left alone, and your (randomised, in the case of RcppDE) optimisation should proceed. Does that sounds reasonable?

@kevinushey

This comment has been minimized.

Contributor

kevinushey commented Jun 6, 2018

Just to be clear, the issue here is that this snippet of code:

library(Rcpp)
library(RcppDE)

sourceCpp(code='
          #include <Rcpp.h>
          
          // [[Rcpp::export]]
          int fibonacci(const int x) {
          if (x == 0) return(0);
          if (x == 1) return(1);
          return (fibonacci(x - 1)) + fibonacci(x - 2);
          }'
)


set.seed(1234)
Rosenbrock2 <- function(x){
    foo <- fibonacci(3)
    x1 <- x[1]
    x2 <- x[2]
    100 * (x2 - x1 * x1)^2 + (1 - x1)^2
}

lower <- c(-10,-10)
upper <- -lower
cont <- DEoptim.control(trace = 0)

set.seed(1234)
fit2 <- DEoptim(Rosenbrock2, lower, upper, control = cont)
fit2$optim

fails to optimize? And the other examples just illustrate situations wherein things still work fine; ie, the breakage is limited to these kinds of usages.

@eddelbuettel

This comment has been minimized.

Owner

eddelbuettel commented Jun 6, 2018

Yes, that's how I read it. Rosenbrock is a standard optimisation problem function, Kyle made it non-standard by having it call fibonacci. If everything was one "FibonacciRosenbrock" it would probably work -- as I understand it the nesting is now broken. And I think it is because we "cleaned up" the RNG states (driven by Joe/Winston).

@kevinushey

This comment has been minimized.

Contributor

kevinushey commented Jun 6, 2018

I think you are correct. With the older versions of Rcpp, GetRNGstate() and PutRNGstate() would only be called once each -- once when entering DEoptim() and once when exiting DEoptim(). Now, it gets called each time the fibonacci() function is invoked, and it looks like that ends up propagating a stale RNG state. I think this would affect all Rcpp functions exposed through attributes (not just sourceCpp()) so we should try to figure out a solution.

Some potential resolutions:

  1. This could potentially be accommodated for in the RcppDE package, by ensuring that the RNG state is synchronized before and after each call to the associated optimization function. This could be expensive depending on how many times the optimization function needs to be run to find the optimum, but at least we know it would be correct.

  2. We could provide some kind of API in Rcpp to 'suspend' updates to the RNG state, and RcppDE could use that to say "when running these functions, even if I'm calling an Rcpp function, please don't attempt to synchronize the RNG state". This could still potentially break callbacks to Rcpp functions that themselves mix usages of the RNG from the R and C sides, so this doesn't seem great, but it does seem unlikely that anyone would intentionally touch the RNG within the optimization function...

Other thoughts? If (1) feels okay, I can open a PR to RcppDE.

@kylebmetrum

This comment has been minimized.

kylebmetrum commented Jun 6, 2018

Apologies for the confusion on the report. Yeah, I was intending to show that the optimization works with no fibonacci and then doesn't work when you add it in; but works again with DEoptim::DEoptim.

I'm doing something similar in principle ... and just now was able to reproduce with the contrived example. In my application, the fibonacci function couldn't be reasonably wrapped into the same function.

Reporting now because this sort of thing used to work. If it isn't going to work going forward, that's cool. I thought it might be pointing to a deeper issue.

Here's what the iterations look like for Rosenbrock2 ... the estimates don't move anywhere from where they start.

Thanks,
Kyle

Iteration: 1 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 2 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 3 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 4 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 5 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 6 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 7 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 8 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 9 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 10 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 11 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 12 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 13 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 14 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 15 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 16 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 17 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 18 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 19 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
Iteration: 20 bestvalit: 327764.081065 bestmemit:   -7.725932    2.445988
@eddelbuettel

This comment has been minimized.

Owner

eddelbuettel commented Jun 6, 2018

@kylebmetrum Thanks for the follow-up
@kevinushey Good discussion, very helpful.

I think we should take a short break and think this through. I do think a case can be made for current behaviour being "as designed" so maybe what use to work for Kyle simply no longer does. Then again, our "technology" of compiling little snippets into R is attractive and has found its way into other places (stan, anyone?). Maybe we should support a general "support updates" feature as you had outlined. I am not sure yet how though.

I think if we make a change it should really be at the Rcpp level and not tucked away here in a somewhat obscure (yet useful ;-) package.

Thoughts?

@kevinushey

This comment has been minimized.

Contributor

kevinushey commented Jun 6, 2018

Here's an example that illustrates the issue can occur when calling any Rcpp Attributes function:

library(Rcpp)
library(RcppArmadillo)
library(RcppDE)

set.seed(1234)
Rosenbrock2 <- function(x) {
  RcppArmadillo:::armadillo_version(FALSE)  # this breaks the optimizer
  x1 <- x[1]
  x2 <- x[2]
  100 * (x2 - x1 * x1) ^ 2 + (1 - x1) ^ 2
}

lower <- c(-10, -10)
upper <- -lower
cont <- DEoptim.control(trace = 0)

set.seed(1234)
fit2 <- DEoptim(Rosenbrock2, lower, upper, control = cont)
fit2$optim

I chose armadillo_version(FALSE) here just because it's a super simple function exported with Rcpp Attributes. Try commenting out the call to RcppArmadillo:::armadillo_version(FALSE) -- you'll see that things work as expected again. Since the issue affects any calls to Rcpp functions (not just those compiled on the spot) I think it warrants getting a fix in.

eddelbuettel added a commit that referenced this issue Jun 9, 2018

Merge pull request #15 from kevinushey/bugfix/suspend-rng-sync
suspend RNG synchronization (fixes #14)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment