rcannood / GillespieSSA2

Gillespie’s Stochastic Simulation Algorithm for Impatient People ⌛
https://rcannood.github.io/GillespieSSA2/
7 stars 0 forks source link

issues with seeds/repetition #8

Closed bbolker closed 3 years ago

bbolker commented 3 years ago

I get identical dynamics when I run a GillespieSSA2 model twice in a row. I can't find any discussion of random seeds (setting/resetting/etc.) in the package/repository, despite the fact that set.seed() is used in the examples. The example below is adapted from the Lotka-Volterra vignette, but I think this is general ... (is some environment being stored/reset in a way that accidentally restores the previous state of the random-number seed?)

library(GillespieSSA2)
sim_name <- "Lotka Predator-Prey model"
params <- c(c1 = 10, c2 = .01, c3 = 10)
final_time <- 2
initial_state <- c(Y1 = 1000, Y2 = 1000)
## Define reactions
reactions <- list(
  reaction("c1 * Y1", c(Y1 = +1)),
  reaction("c2 * Y1 * Y2", c(Y1 = -1, Y2 = +1)),
  reaction("c3 * Y2", c(Y2 = -1))
)
## Run simulations with the Exact method
set.seed(1)
out1 <- ssa(
  initial_state = initial_state,
  reactions = reactions,
  params = params,
  final_time = final_time,
  method = ssa_exact(),
  sim_name = sim_name
) 
out2 <- ssa(
  initial_state = initial_state,
  reactions = reactions,
  params = params,
  final_time = final_time,
  method = ssa_exact(),
  sim_name = sim_name
) 
identical(out1$state,out2$state)  # TRUE

Using set.seed(2) before the second run does make the results different, but unless I'm very confused (always a possibility) this feels like something that needs to be documented ??

rcannood commented 3 years ago

Thanks for bringing this to my attention!

It seemed I was calling R's RNG functions without calling all of the required functions as documented in section 6.3 of Writing R Extensions. I used Rcpp's RNGScope (as discussed by Dirk Eddelbuettel here) to solve the problem.

Installing the latest release of GillespieSSA2 should resolve this issue.