R/rstream_mrg32k3a.R

In rstream: Streams of Random Numbers

#############################################################################
##                                                                         ##
##   Class: rstream.mrg32k3a                                               ##
##                                                                         ##
#############################################################################
##                                                                         ##
##   Interface to Pierre L'Ecuyers RngStreams Library                      ##
##                                                                         ##
#############################################################################


## Initialize global variables ----------------------------------------------

.rstream.mrg32k3a.init <- function () {

	## Indicate whether a package seed is already set
	if (!exists(".rstream.mrg32k3a.HasSeed", envir=.rstream.envir))
		assign(".rstream.mrg32k3a.HasSeed", FALSE, envir=.rstream.envir)

	## Default Package seed
	if (!exists(".rstream.mrg32k3a.DefaultSeed", envir=.rstream.envir))
		assign(".rstream.mrg32k3a.DefaultSeed", rep(as.double(Sys.time()),6), envir=.rstream.envir)
}


## Class --------------------------------------------------------------------

setClass( "rstream.mrg32k3a", representation(), contains = "rstream" )


## Initialize ---------------------------------------------------------------

setMethod( "initialize", "rstream.mrg32k3a",  
          function(.Object, name=NULL, seed=NULL, force.seed=FALSE, 
                   antithetic=FALSE, incprecision=FALSE) {

                  ## first increment counter for Rstreams
                  count <- get(".rstream.Count", envir=.rstream.envir) + 1
                  assign(".rstream.Count", count, envir=.rstream.envir)

                  ## name of the Rstream object.
                  ## by default we use type + number
                  if (is.null(name)) name <- paste("mrg32k3a", count, sep="")

                  ## type of Rstream object
                  .Object@type <- "mrg32k3a"

                  ## add info about Rstream type
                  .Object@info <- "RngStreams - library for multiple independent streams of Random Numbers"

                  ## at creation a Rstream is never packed
                  .Object@is.packed <- FALSE 

                  ## set (package) seed
                  hasSeed <- get(".rstream.mrg32k3a.HasSeed", envir=.rstream.envir)
                  defaultSeed <- get(".rstream.mrg32k3a.DefaultSeed", envir=.rstream.envir)
                  if (!is.null(seed) && !force.seed && hasSeed)
                          stop("rstream.mrg32k3a: Package already seeded! Seed ignored!\n",
                               "rstream.mrg32k3a: Set force.seed = TRUE if you really want to reseed.")
                  if (is.null(seed)) seed <- defaultSeed
                  seed <- .rstream.mrg32k3a.CheckSeed(seed)
                  .Call(C_R_RngStreams_SetPackageSeed, as.double(seed))

                  ## Create Rstream object
                  .Object@stream <- .Call(C_R_RngStreams_Init, .Object, as.character(name))
                  .Call(C_R_RngStreams_SetAntithetic, .Object@stream, as.integer(antithetic))
                  .Call(C_R_RngStreams_SetIncreasedPrecis, .Object@stream, as.integer(incprecision))

                  ## save seed as R variavble
                  assign(".rstream.mrg32k3a.DefaultSeed",
                         as.numeric(.Call(C_R_RngStreams_GetPackageSeed)), envir=.rstream.envir)
                  assign(".rstream.mrg32k3a.HasSeed", TRUE, envir=.rstream.envir)

                  ## return new Rstream object
                  .Object
          } )


## Validity -----------------------------------------------------------------

## .rstream.mrg32k3a.CheckSeed
##    make simple check on given seed
.rstream.mrg32k3a.CheckSeed <- function(seed) {
	ll <- length(seed)
	if (ll < 6)
		stop("rstream.mrg32k3a: seed too short; 6 numbers required")
	if (ll > 6) {
		warning("rstream.mrg32k3a: seed too long; truncated") 
		seed<-seed[1:6] }
###	for (i in 1:6)
###		if (seed[i] < 0) 
###			stop("invalid seed")
	seed
}


## Methods ------------------------------------------------------------------

## Access and Replacement methods ...........................................

## rstream.name
##    get and set name of Rstream object
setMethod("rstream.name", "rstream.mrg32k3a", 
          function(stream) { 
                  if (stream@is.packed) 
                          return (stream@pack$name)
                  else
                          return (.Call(C_R_RngStreams_GetName, stream@stream))
          } )

setReplaceMethod("rstream.name", "rstream.mrg32k3a",
                 function(stream, value) {
                         if (stream@is.packed) stop("Cannot change name for PACKED Rstream") 
                         .Call(C_R_RngStreams_SetName, stream@stream, as.character(value))
                         stream
                 } )


## rstream.antithetic
##   get and set flag for antithetic random numbers:  
setMethod("rstream.antithetic", "rstream.mrg32k3a", 
          function(stream) {
                  if (stream@is.packed) 
                          return (stream@pack$anti)
                  else 
                          return (as.logical(.Call(C_R_RngStreams_GetAntithetic, stream@stream)))
          } )

setReplaceMethod("rstream.antithetic", "rstream.mrg32k3a",
                 function(stream, value) { 
                         if (stream@is.packed) stop("Cannot change antithetic flag for PACKED Rstream") 
                         .Call(C_R_RngStreams_SetAntithetic, stream@stream, as.integer(value))
                         stream
                 } )


## rstream.incprecision
##    get and set flag for increased precision of random numbers
setMethod("rstream.incprecision", "rstream.mrg32k3a", 
          function(stream) {
                  if (stream@is.packed) 
                          return (stream@pack$incp)
                  else 
                          return (as.logical(.Call(C_R_RngStreams_GetIncreasedPrecis, stream@stream)))
          } )

setReplaceMethod("rstream.incprecision", "rstream.mrg32k3a",
                 function(stream, value) { 
                         if (stream@is.packed) stop("Cannot change increased precision flag for PACKED Rstream") 
                         .Call(C_R_RngStreams_SetIncreasedPrecis, stream@stream, as.integer(value))
                         stream
                 } )


## Sampling methods .........................................................

## rstream.sample
##    make a random sample
setMethod("rstream.sample", "rstream.mrg32k3a",
          function(stream,n=1) { 
                  if (stream@is.packed) stop("Cannot sample from PACKED Rstream") 
                  .Call(C_R_RngStreams_Sample, stream@stream, as.integer(n)) } )

setMethod("r", "rstream.mrg32k3a",
          function(stream,n=1) { 
                  if (stream@is.packed) stop("Cannot sample from PACKED Rstream") 
                  .Call(C_R_RngStreams_Sample, stream@stream, as.integer(n)) } )


## Jump methods .............................................................

## rstream.resetsubstream
##   reset current substream of Rstream object
if(!isGeneric("rstream.resetsubstream"))
        setGeneric("rstream.resetsubstream", function(stream) standardGeneric("rstream.resetsubstream"))
setMethod("rstream.resetsubstream", "rstream.mrg32k3a", 
          function(stream) { 
                  if (stream@is.packed) stop("Cannot reset PACKED Rstream") 
                  dummy <- .Call(C_R_RngStreams_ResetStartSubstream, stream@stream)
          } )


## rstream.nextsubstream
##   skip to next substream of Rstream object
if(!isGeneric("rstream.nextsubstream"))
        setGeneric("rstream.nextsubstream", function(stream) standardGeneric("rstream.nextsubstream"))
setMethod("rstream.nextsubstream", "rstream.mrg32k3a", 
          function(stream) { 
                  if (stream@is.packed) stop("Cannot skip substream of PACKED Rstream") 
                  dummy <- .Call(C_R_RngStreams_ResetNextSubstream, stream@stream)
          } )


## Reset and copy methods ...................................................

## rstream.reset
##   reset Rstream object
setMethod("rstream.reset", "rstream.mrg32k3a", 
          function(stream) { 
                  if (stream@is.packed) stop("Cannot reset PACKED Rstream") 
                  dummy <- .Call(C_R_RngStreams_ResetStartStream, stream@stream)
          } )


## rstream.clone
##    clone (copy) Rstream object
setMethod("rstream.clone", "rstream.mrg32k3a", 
          function(stream) { 
                  if (stream@is.packed) stop("Cannot clone PACKED Rstream") 
                  clone <- stream
                  name <- paste(rstream.name(stream),".",sep="")
                  clone@stream <- .Call(C_R_RngStreams_Clone, clone, stream@stream, name)
                  clone
          } )


## rstream.pack, rstream.unpack
##    pack and unpack Rstream object such that all data are contained in R object
##    (and can be easily copied within R)
setReplaceMethod("rstream.packed", "rstream.mrg32k3a", 
                 function(stream, value) {
                         value <- as.logical(value)
                         if (value && stream@is.packed)   return(stream)
                         if (!value && !stream@is.packed) return(stream)
                         if (value) {	# pack
                                 name <- rstream.name(stream)
                                 anti <- rstream.antithetic(stream)
                                 incp <- rstream.incprecision(stream)
                                 stream@is.packed <- TRUE
                                 stream@pack <- list()
                                 stream@pack$state <- as.double(.Call(C_R_RngStreams_GetData, stream@stream))
                                 stream@pack$name <- name 
                                 stream@pack$anti <- anti
                                 stream@pack$incp <- incp
                                 .Call(C_R_RngStreams_Free, stream@stream)
                         }
                         else {		# unpack
                                 stream@is.packed <- FALSE
                                 .Call(C_R_RngStreams_SetData, stream,
                                       stream@stream, stream@pack$state, stream@pack$name)
                         }
                         stream
                 } )


## Printing et al. ..........................................................

## print:
##    print all data of a Rstream object
setMethod( "print", "rstream.mrg32k3a",
          function(x, ...) { 
                  .rstream.PrintData(x)
                  if (!x@is.packed) {
                          state <- .Call(C_R_RngStreams_GetData, x@stream)
                          cat("\n\tInternal state:\n",
                              "\t  initial state:\n",
                              "\t\t", state[13:18], "\n",
                              "\t  starting point of current substream:\n",
                              "\t\t", state[7:12], "\n",
                              "\t  current state:\n",
                              "\t\t", state[1:6], "\n\n" ) }
                  else
                          cat("\n")
          } )


## Rstream objects <-> R generators -----------------------------------------

## .rstream.getRNG
##    get Rstream object for current R generator
##    (internal method; not exported)
setMethod(".rstream.getRNG", "rstream.mrg32k3a", 
          function(stream) { rstream.clone(stream) } )


## .rstream.setRNG
##    set R generator to given Rstream object
##    (internal method; not exported)
setMethod(".rstream.setRNG", "rstream.mrg32k3a", 
          function(stream) {
            if (isTRUE(get(".rstream.version", envir=.rstream.envir) < 1003)) {
              ## pre 1.3:
              .Call(C_R_RngStreams_setRNG, stream@stream)
              RNGkind(kind="user-supplied")
            } else {
              ## new version:
              RNGkind(kind="user-supplied")
              .Call(C_R_RngStreams_setRNG, stream@stream)
            }
            stream
          } )

## End ----------------------------------------------------------------------