R/RNGseq.R

In rngtools: Utility functions for working with Random Number Generators

# Generate a sequence of RNGs suitable for parallel computation 
# using L'Ecuyer's RNG
# 
# Author: Renaud Gaujoux
###############################################################################

# or-NULL operator (borrowed from Hadley Wickham)
'%||%' <- function(x, y) if( !is.null(x) ) x else y

#' Generate Sequence of Random Streams
#' 
#' Create a given number of seeds for L'Ecuyer's RNG, that can be used to seed 
#' parallel computation, making them fully reproducible.
#' 
#' This ensures complete reproducibility of the set of run. 
#' The streams are created using L'Ecuyer's RNG, implemented in R core since
#' version 2.14.0 under the name \code{"L'Ecuyer-CMRG"} (see \code{\link{RNG}}).
#' 
#' Generating a sequence without specifying a seed uses a single draw of the 
#' current RNG. The generation of a sequence using seed (a single or 6-length 
#' numeric) a should not affect the current RNG state. 
#' 
#' @param n Number of streams to be created
#' @param seed seed specification used to initialise the set of streams 
#' using \code{\link{RNGseq_seed}}.
#' @param simplify a logical that specifies if sequences of length 1 should be
#' unlisted and returned as a single vector.
#' @param ... extra arguments passed to \code{\link{RNGseq_seed}}.  
#' 
#' @return a list of integer vectors (or a single integer  vector if 
#' \code{n=1} and \code{unlist=TRUE}).
#' 
#' @export 
#' @examples
#' 
#' RNGseq(3)
#' RNGseq(3)
#' RNGseq(3, seed=123)
#' # or identically
#' set.seed(123)
#' identical(RNGseq(3), RNGseq(3, seed=123))
#' \dontshow{
#' set.seed(123)
#' stopifnot( identical(RNGseq(3), RNGseq(3, seed=123)) ) 
#' }
#' 
#' RNGseq(3, seed=1:6, verbose=TRUE)
#' # select Normal kind
#' RNGseq(3, seed=123, normal.kind="Ahrens")
#' 
RNGseq <- function(n, seed=NULL, ..., simplify=TRUE, version=2){
	
	library(parallel)
	# check parameters
	if( n <= 0 )
		stop("NMF::createStream - invalid value for 'n' [positive value expected]")
	
	# extract RNG setting from object if possible
	if( !is.null(seed) ) seed <- getRNG(seed, num.ok=TRUE) %||% seed
	
	# convert matrix into a list of seed
	if( is.matrix(seed) )
		seed <- lapply(seq(ncol(seed)), function(i) seed[,i])
	
	# if already a sequence of seeds: use directly
	#print(seed)
	if( is.list(seed) ){
		# check length
		if( length(seed) > n ){
			warning("Reference seed sequence is longer than the required number of seed: only using the ", n, " first seeds.")
			seed <- seed[1:n]
		}else if( length(seed) < n )
			stop("Reference seed sequence is shorter [",length(seed),"] than the required number of seed [", n, "].")
		
		res <- lapply(seed, as.integer)
	}else{ # otherwise: get initial seed for the CMRG stream sequence
		
		orng <- RNGseed()
		.s <- RNGseq_seed(seed, ..., version=version)
		
		res <- lapply(1:n, function(i){
					if( i == 1 ) .s
					else .s <<- nextRNGStream(.s)				
				})
		
		# if not seeded and current RNG is L'Ecuyer-CMRG => move to stream after last stream
		if( is.null(seed) && RNGkind()[1L] == "L'Ecuyer-CMRG" && version>=2 ){
			# ensure old normal kind is used 
			RNGseed(c(orng[1L], nextRNGStream(.s)[2:7]))
		}
	}
	
	# return list or single RNG
	if( n==1 && simplify )
		res[[1]]
	else
		res
	
}

#' \code{RNGseq_seed} generates the -- next -- random seed used as the first seed in 
#' the sequence generated by \code{\link{RNGseq}}. 
#' 
#' @param normal.kind Type of Normal random generator. See \code{\link{RNG}}.
#' @param verbose logical to toggle verbose messages
#' @param version version of the function to use, to reproduce old behaviours.
#' Version 1 had a bug which made the generated stream sequences share most of their
#' seeds (!), as well as being not equivalent to calling \code{set.seed(seed); RNGseq_seed(NULL)}. 
#' Version 2 fixes this bug.
#' 
#' @return a 7-length numeric vector.
#' @seealso \code{\link{RNGseq}}
#' 
#' @rdname RNGseq
#' @export
#' @examples 
#' 
#' ## generate a seed for RNGseq
#' # random  
#' RNGseq_seed() 
#' RNGseq_seed()
#' RNGseq_seed(NULL)
#' # fixed
#' RNGseq_seed(1)
#' RNGseq_seed(1:6)
#' 
#' # `RNGseq_seed(1)` is identical to 
#' set.seed(1)
#' s <- RNGseq_seed()
#' identical(s, RNGseq_seed(1))
#' \dontshow{ stopifnot(identical(s, RNGseq_seed(1))) }
#'  
RNGseq_seed <- function(seed=NULL, normal.kind=NULL, verbose=FALSE, version=2){
	
	# retrieve current seed
	orng <- RNGseed()
	# setup RNG restoration in case of an error
	on.exit({
		RNGseed(orng)
		if( verbose ) message("# Restoring RNG: ", paste(RNGkind(), collapse=' - '), ' [', .collapse(orng), ']')
	})
	
	rkind_not_CMRG <- RNGkind()[1L] != "L'Ecuyer-CMRG"
	
	if( verbose ) message("# Original RNG: ", paste(RNGkind(), collapse=' - '), ' [', .collapse(orng), ']')
	# seed with numeric seed
	if( is.numeric(seed) ){
		if( length(seed) == 1L ){
			
			if( verbose ) message("# Generate RNGstream random seed from ", seed, " ... ", appendLF=FALSE)
			if( version<2 || rkind_not_CMRG ){ # behaviour prior 1.4
				set.seed(seed)
				RNGkind(kind="L'Ecuyer-CMRG", normal.kind=normal.kind)
			}else{ # fix seed after switching to CMRG: ensure result independence from the current RNG
				set.seed(seed, kind="L'Ecuyer-CMRG", normal.kind=normal.kind)
			}
			if( verbose ) message("OK")
		}
		else if( length(seed) == 6L ){
			if( verbose ) message("# Directly use 6-long seed: ",  paste(seed, collapse=', '), " ... ", appendLF=FALSE)
			RNGkind("L'Ecuyer-CMRG", normal.kind=normal.kind)
			s <- RNGseed()
			s[2:7] <- as.integer(seed)
			RNGseed(s)
			if( verbose ) message("OK")
		}else if ( length(seed) == 7L ){
			if( seed[1] %% 100 != 7L )
				stop("RNGseq_seed - Invalid 7-long numeric seed: RNG code should be '7', i.e. of type \"L'Ecuyer-CMRG\"")
			if( verbose ) message("# Directly use CMRG seed: ",  paste(seed, collapse=', '), " ... ", appendLF=FALSE)
			RNGseed(seed)
			if( verbose ) message("OK")
		}else
			stop("RNGseq_seed - Invalid numeric seed: should be a numeric of length 1, 6 or 7")		
	}else if( is.null(seed) ){
		if( rkind_not_CMRG ){ # seed with random seed
			
			# draw once from the current calling RNG to ensure different seeds
			# for separate loops, but to ensure identical results as with set.seed
			# one must still use the current RNG before changing RNG kind
			runif(1)
			orng1 <- RNGseed()
			RNGseed(orng)
			orng <- orng1
			
			if( verbose ) message("# Generate random RNGstream seed: ", appendLF=FALSE)
			RNGkind(kind="L'Ecuyer", normal.kind=normal.kind)
			if( verbose ) message("OK")
		}else{ # seed with next RNG stream
			if( version < 2 ){
				on.exit() # cancel RNG restoration
				s <- nextRNGStream(orng)
				if( verbose ) message("# Use next active RNG stream: ", .collapse(s[2:7]))
				RNGseed(s)	
			}else{
				# only change normal kind if necessary and use current stream state
				if( !is.null(normal.kind) ) RNGkind(normal.kind=normal.kind)
				if( verbose ) message("# Use current active RNG stream: ", .collapse(RNGseed()[2:7]))
			}
		}
	}else
		stop("RNGseq_seed - Invalid seed value: should be a numeric or NULL")
	
	s <- RNGseed()
	if( verbose ) message("# Seed RNGkind is: ", paste(RNGkind(), collapse=' - '), ' [', .collapse(s), ']')	
	s	
}