[JOSS review] Seeding

[mllg]

This is a part of the JOSS review outlined in openjournals/joss-reviews#1493.

The seeds currently default to 1:njobs. I've followed a similar approach in batchtools where I determine a start seed randomly, and then increment the seed for each job. Even with a random initialization, I'm not sure if this approach is feasible (see https://github.com/mllg/batchtools/issues/81). Defaulting to 1:njobs seems to be even worse, there is a reason RNGs are not initialized with a constant value.

[gvegayon]

I really liked this post https://rpubs.com/Jouni_Helske/225931 (mentioned in the issue you pointed out). Following one of the (rather simple) tests, which as you suggested the simple 1:njobs does not pass, I wrote this simple (sort of arbitrary) function:

# A naive seeding function
new_seeds <- function(nseeds) {

  sum_of_seed <- sum(.Random.seed)
  sum_of_seed <- sum_of_seed - (sum_of_seed %/% 100L)*100L
  seeds <- (1L:nseeds) * 10L
  (seeds + sum_of_seed ^ (1L:nseeds %% 2L))[order(stats::runif(nseeds))]

}

# "Test" from https://rpubs.com/Jouni_Helske/225931
last_rand_seq <- list()
random <- function(kind) {
  set.seed(123)
  RNGkind(kind = kind)
  x <- numeric(1e5)
  seeds <- new_seeds(1e5)
  for(i in seq_along(x)){
    set.seed(seeds[i])
    x[i] <- runif(1)
  }

  last_rand_seq <<- c(last_rand_seq, list(x))

  acf(x, main = kind)
}
kinds <- c("Mersenne-Twister", "Marsaglia-Multicarry", "L'Ecuyer-CMRG")

invisible(sapply(kinds, random)) # Looks nice, no autocorrelation

This function is reproducible, i.e., seeds can be reproduced. Will including this function to generate the random seeds instead of using 1:njobs be enough for JOSS? From what I read in the discussion this is still a problem without a clear solution.

[mllg]

I'm not sure I understand how the function new_seeds() helps to solve the problem

Your seeds are determined by the current RNG state on which you apply some transformations. Wouldn't you get the same quality of randomness with as.integer(runif(nseeds, 1, .Machine$integer.max))?

Either way, I think that the proposed solution (or just calling runif()) would be sufficient for me -- the arguably better solution is the one proposed by @HenrikBengtsson in the linked issue.

Note that there are two special cases which you might also need to address to get the seeding right:

1) .Random.seed is NULL until the RNG is initiated (usually by drawing the first random number). 2) sum(.Random.seed) might overflow.