1# Copyright (C) 2009-2012 Renaud Gaujoux 2# 3# This file is part of the rngtools package for R. 4# This program is free software; you can redistribute it and/or 5# modify it under the terms of the GNU General Public License 6# as published by the Free Software Foundation; either version 2 7# of the License, or (at your option) any later version. 8# 9# This program is distributed in the hope that it will be useful, 10# but WITHOUT ANY WARRANTY; without even the implied warranty of 11# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12# GNU General Public License for more details. 13# 14# You should have received a copy of the GNU General Public License 15# along with this program; if not, write to the Free Software 16# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 17# 18# Creation: 08 Nov 2011 19############################################################################### 20 21###% Returns all the libraries that provides a user-supplied RNG 22###% 23###% The library that provides the wrapper hooks for the management multiple 24###% user-supplied RNG is removed from the output list. 25###% 26#' @importFrom utils tail 27RNGlibs <- function(n=0, full=FALSE, hook="user_unif_rand", unlist=TRUE){ 28 dlls <- getLoadedDLLs() 29 res <- lapply(dlls, function(d){ 30 dname <- d[['name']] 31 if( dname=='' ) 32 return(NA) 33 34 symb.unif_rand <- RNGlib(PACKAGE=dname, hook=hook) 35 if( is.null(symb.unif_rand) ) 36 NA 37 else 38 symb.unif_rand 39 }) 40 41 res <- res[!is.na(res)] 42 if( !full ) 43 res <- names(res) 44 45 # limit the results if requested 46 if( n>0 ) 47 res <- tail(res, n) 48 49 # return result 50 if( unlist && length(res) == 1 ) 51 res[[1]] 52 else 53 res 54} 55 56###% Returns the library that provides the current user-supplied RNG hooks. 57###% 58###% This is the library that is first called by runif when using setting RNG 59###% kind to "user-supplied". 60###% In general this will be rstream, except if a package providing the RNG hook 61###% 'user_unif_rand' is loaded after rstream, and no call to RNGkind or getRNG 62###% were done thereafter. 63###% 64###% @return an object of class NativeSymbolInfo or NULL if no hook were found 65###% 66RNGlib <- function(PACKAGE='', full=FALSE, hook="user_unif_rand", ...){ 67 68 if( !missing(PACKAGE) ) 69 full = TRUE 70 if( !missing(hook) ) 71 hook <- match.arg(hook, c('user_unif_rand', 'user_unif_init', 'user_unif_nseed', 'user_unif_seedloc')) 72 73 # lookup for the hook "user_unif_rand" in all the loaded libraries 74 symb.unif_rand <- try( getNativeSymbolInfo(hook, PACKAGE=PACKAGE, ...), silent=TRUE) 75 if( is(symb.unif_rand, 'try-error') ){ 76 77 if( !full ) '' else NULL 78 79 }else if( PACKAGE=='' && is.null(symb.unif_rand$package) ){ 80 #special case for MS Windows when PACKAGE is not specified: if two 81 # RNGlibs are loaded, the first one is seen, not the last one as on Unix 82 libs <- RNGlibs(full=TRUE, unlist=FALSE, hook=hook) 83 w <- which(sapply(libs, function(l) identical(l$address, symb.unif_rand$address))) 84 85 # returns full info or just the name 86 if( full ) libs[[w]] 87 else names(libs)[w] 88 89 }else if( full ) symb.unif_rand 90 else symb.unif_rand$package[['name']] 91} 92 93###% Returns the package that provides the current RNG managed by rstream 94###% 95###% It returns the name of the package to which are currently passed the RNG 96###% calls (runif, set.seed). 97###% This is either 'base' if core RNG is in use (e.g. Mersenne-Twister, Marsaglia-Multicarry, etc...) 98###% or the package that provides the actual RNG hooks called by the rstream 99###% wrapper hooks. This one was set either explicitly via RNGkind or implicitly 100###% when rstream was first loaded. In this latter case, the provider was identified 101###% at loading time as 'base' if core RNGs were in use or as the package that was 102###% providing the RNG hook 'user_unif_rand' if the RNG in used was "user-supplied". 103###% 104RNGprovider <- function(kind=RNGkind(), user.supplied=FALSE){ 105 106 if( kind[1L] == 'user-supplied' || user.supplied ) RNGlib() 107 else 'base' 108} 109 110#' Directly Getting or Setting the RNG Seed 111#' 112#' These functions provide a direct access to the RNG seed object `.Random.seed`. 113#' 114#' @name RNGseed 115NULL 116 117#' @describeIn RNGseed directly gets/sets the current RNG seed \code{.Random.seed}. 118#' It can typically be used to backup and restore the RNG state on exit of 119#' functions, enabling local RNG changes. 120#' 121#' @param seed an RNG seed, i.e. an integer vector. 122#' No validity check is performed, so it \strong{must} be a 123#' valid seed. 124#' 125#' @return invisibly the current RNG seed when called with no arguments, 126#' or the -- old -- value of the seed before changing it to 127#' \code{seed}. 128#' 129#' @export 130#' @examples 131#' 132#' # get current seed 133#' RNGseed() 134#' # directly set seed 135#' old <- RNGseed(c(401L, 1L, 1L)) 136#' # show old/new seed description 137#' showRNG(old) 138#' showRNG() 139#' 140#' # set bad seed 141#' RNGseed(2:3) 142#' try( showRNG() ) 143#' # recover from bad state 144#' RNGrecovery() 145#' showRNG() 146#' 147#' # example of backup/restore of RNG in functions 148#' f <- function(){ 149#' orng <- RNGseed() 150#' on.exit(RNGseed(orng)) 151#' RNGkind('Marsaglia') 152#' runif(10) 153#' } 154#' 155#' sample(NA) 156#' s <- .Random.seed 157#' f() 158#' identical(s, .Random.seed) 159#' \dontshow{ stopifnot(identical(s, .Random.seed)) } 160#' 161RNGseed <- function(seed){ 162 163 res <- if( missing(seed) ){ 164 if( exists('.Random.seed', where = .GlobalEnv) ) 165 get('.Random.seed', envir = .GlobalEnv) 166 }else if( is.null(seed) ){ 167 if( exists('.Random.seed', where = .GlobalEnv) ) 168 rm('.Random.seed', envir = .GlobalEnv) 169 }else{ 170 old <- RNGseed() 171 assign('.Random.seed', seed, envir = .GlobalEnv) 172 old 173 } 174 invisible(res) 175} 176 177#' @describeIn RNGseed recovers from a broken state of \code{.Random.seed}, 178#' and reset the RNG settings to defaults. 179#' 180#' @export 181RNGrecovery <- function(){ 182 s <- as.integer(c(401,0,0)) 183 assign(".Random.seed", s, envir=.GlobalEnv) 184 do.call(RNGkind, as.list(rep("default", .RNGkind_length()))) 185 186} 187 188.getRNGattribute <- function(object){ 189 if( .hasSlot(object, 'rng') ) slot(object, 'rng') 190 else if( .hasSlot(object, 'rng.seed') ) slot(object, 'rng.seed') # for back compatibility 191 else if( !is.null(attr(object, 'rng')) ) attr(object, 'rng') 192 else if( is.list(object) ){ # compatibility with package setRNG 193 if( !is.null(object[['rng']]) ) object[['rng']] 194 else if( is.list(object[['noise']]) && !is.null(object[['noise']][['rng']]) ) 195 object[['noise']][['rng']] 196 }else NULL 197} 198 199#' Getting/Setting RNGs 200#' 201#' \code{getRNG} returns the Random Number Generator (RNG) settings used for 202#' computing an object, using a suitable \code{.getRNG} S4 method to extract 203#' these settings. 204#' For example, in the case of objects that result from multiple model fits, 205#' it would return the RNG settings used to compute the best fit. 206#' 207#' This function handles single number RNG specifications in the following way: 208#' \describe{ 209#' \item{integers}{Return them unchanged, considering them as encoded RNG kind 210#' specification (see \code{\link{RNG}}). No validity check is performed.} 211#' \item{real numbers}{If \code{num.ok=TRUE} return them unchanged. 212#' Otherwise, consider them as (pre-)seeds and pass them to \code{\link{set.seed}} 213#' to get a proper RNG seed. 214#' Hence calling \code{getRNG(1234)} is equivalent to \code{set.seed(1234); getRNG()} 215#' (See examples). 216#' } 217#' } 218#' 219#' @param object an R object from which RNG settings can be extracted, e.g. an 220#' integer vector containing a suitable value for \code{.Random.seed} or embedded 221#' RNG data, e.g., in S3/S4 slot \code{rng} or \code{rng$noise}. 222#' @param ... extra arguments to allow extension and passed to a suitable S4 method 223#' \code{.getRNG} or \code{.setRNG}. 224#' @param num.ok logical that indicates if single numeric (not integer) RNG data should be 225#' considered as a valid RNG seed (\code{TRUE}) or passed to \code{\link{set.seed}} 226#' into a proper RNG seed (\code{FALSE}) (See details and examples). 227#' @param extract logical that indicates if embedded RNG data should be looked for and 228#' extracted (\code{TRUE}) or if the object itself should be considered as an 229#' RNG specification. 230#' @param recursive logical that indicates if embedded RNG data should be extracted 231#' recursively (\code{TRUE}) or only once (\code{FASE}). 232#' 233#' @return \code{getRNG}, \code{getRNG1}, \code{nextRNG} and \code{setRNG} 234#' usually return an integer vector of length > 2L, like \code{\link{.Random.seed}}. 235#' 236#' \code{getRNG} and \code{getRNG1} return \code{NULL} if no RNG data was found. 237#' 238#' @rdname rng 239#' @seealso \code{\link{.Random.seed}}, \code{\link{showRNG}} 240#' @export 241#' 242#' @examples 243#' # get current RNG settings 244#' s <- getRNG() 245#' head(s) 246#' showRNG(s) 247#' 248#' # get RNG from a given single numeric seed 249#' s1234 <- getRNG(1234) 250#' head(s1234) 251#' showRNG(s1234) 252#' # this is identical to the RNG seed as after set.seed() 253#' set.seed(1234) 254#' identical(s1234, .Random.seed) 255#' # but if num.ok=TRUE the object is returned unchanged 256#' getRNG(1234, num.ok=TRUE) 257#' 258#' # single integer RNG data = encoded kind 259#' head(getRNG(1L)) 260#' 261#' # embedded RNG data 262#' s <- getRNG(list(1L, rng=1234)) 263#' identical(s, s1234) 264#' 265getRNG <- function(object, ..., num.ok=FALSE, extract=TRUE, recursive=TRUE){ 266 267 if( missing(object) || is.null(object) ) return( .getRNG() ) 268 269 # use RNG data from object if available 270 if( extract && !is.null(rng <- .getRNGattribute(object)) ){ 271 if( recursive && hasRNG(rng) ) getRNG(rng, ..., num.ok=num.ok) 272 else rng 273 }else if( isNumber(object) ){ 274 if( num.ok && isReal(object) ) object 275 else if( isInteger(object) ) object 276 else nextRNG(object, ...) # return RNG as if after setting seed 277 }else .getRNG(object, ...) # call S4 method on object 278 279} 280 281#' \code{hasRNG} tells if an object has embedded RNG data. 282#' @export 283#' @rdname rng 284#' 285#' @examples 286#' # test for embedded RNG data 287#' hasRNG(1) 288#' hasRNG( structure(1, rng=1:3) ) 289#' hasRNG( list(1, 2, 3) ) 290#' hasRNG( list(1, 2, 3, rng=1:3) ) 291#' hasRNG( list(1, 2, 3, noise=list(1:3, rng=1)) ) 292#' 293hasRNG <- function(object){ 294 !is.null(.getRNGattribute(object)) 295} 296 297#' Getting RNG Seeds 298#' 299#' \code{.getRNG} is an S4 generic that extract RNG settings from a variety of 300#' object types. 301#' Its methods define the workhorse functions that are called by \code{getRNG}. 302#' 303#' @inheritParams getRNG 304#' @export 305setGeneric('.getRNG', function(object, ...) standardGeneric('.getRNG') ) 306 307#' @describeIn .getRNG Default method that tries to extract RNG information from \code{object}, by 308#' looking sequentially to a slot named \code{'rng'}, a slot named \code{'rng.seed'} 309#' or an attribute names \code{'rng'}. 310#' 311#' It returns \code{NULL} if no RNG data was found. 312setMethod('.getRNG', 'ANY', 313 function(object, ...){ 314 if( missing(object) ) .get_Random.seed() # safe-guard 315 else .getRNGattribute(object) 316 } 317) 318 319.get_Random.seed <- function(){ 320 # return current value of .Random.seed 321 # ensuring it exists first 322 if( !exists('.Random.seed', envir = .GlobalEnv) ) 323 sample(NA) 324 325 return( get('.Random.seed', envir = .GlobalEnv) ) 326 327} 328#' @describeIn .getRNG Returns the current RNG settings. 329setMethod('.getRNG', 'missing', 330 function(object){ 331 .get_Random.seed() 332 333 } 334) 335 336#' @describeIn .getRNG Method for S3 objects, that aims at reproducing the behaviour of the function 337#' \code{getRNG} of the package \code{getRNG}. 338#' 339#' It sequentially looks for RNG data in elements \code{'rng'}, \code{noise$rng} 340#' if element \code{'noise'} exists and is a \code{list}, or in attribute \code{'rng'}. 341#' 342setMethod('.getRNG', 'list', 343 function(object){ 344 # lookup for some specific elements 345 if( !is.null(object$rng) ) object$rng 346 else if( is.list(object$noise) ) object$noise$rng 347 else attr(object, 'rng') 348 } 349) 350#setMethod('.getRNG', 'rstream', 351# function(object){ 352# object 353# } 354#) 355#' @describeIn .getRNG Method for numeric vectors, which returns the object itself, coerced into an integer 356#' vector if necessary, as it is assumed to already represent a value for 357#' \code{\link{.Random.seed}}. 358#' 359setMethod('.getRNG', 'numeric', 360 function(object, ...){ 361 as.integer(object) 362 } 363) 364 365#' Extracting RNG Settings from Computation Result Objects 366#' 367#' \code{getRNG1} is an S4 generic that returns the \strong{initial} RNG settings 368#' used for computing an object. 369#' For example, in the case of results from multiple model fits, it would 370#' return the RNG settings used to compute the \emph{first} fit. 371#' 372#' \code{getRNG1} is defined to provide separate access to the RNG settings as 373#' they were at the very beginning of a whole computation, which might differ 374#' from the RNG settings returned by \code{getRNG}, that allows to reproduce the 375#' result only. 376#' 377#' Think of a sequence of separate computations, from which only one result is 378#' used for the result (e.g. the one that maximizes a likelihood): 379#' \code{getRNG1} would return the RNG settings to reproduce the complete sequence 380#' of computations, while \code{getRNG} would return the RNG settings necessary to 381#' reproduce only the computation whose result has maximum likelihood. 382#' 383#' @param object an R object. 384#' @param ... extra arguments to allow extension. 385#' 386#' @export 387setGeneric('getRNG1', function(object, ...) standardGeneric('getRNG1') ) 388#' @describeIn getRNG1 Default method that is identical to \code{getRNG(object, ...)}. 389setMethod('getRNG1', 'ANY', 390 function(object, ...){ 391 getRNG(object, ...) 392 } 393) 394 395#' \code{nextRNG} returns the RNG settings as they would be after seeding with 396#' \code{object}. 397#' 398#' @param ndraw number of draws to perform before returning the RNG seed. 399#' 400#' @rdname rng 401#' @export 402#' @examples 403#' head(nextRNG()) 404#' head(nextRNG(1234)) 405#' head(nextRNG(1234, ndraw=10)) 406#' 407nextRNG <- function(object, ..., ndraw=0L){ 408 409 # get/restore .Random.seed on.exit 410 orseed <- RNGseed() 411 on.exit(RNGseed(orseed)) 412 413 # return next state of current RNG if object is missing 414 if( missing(object) ){ 415 runif(1) 416 return( getRNG() ) 417 } 418 419 # extract RNG from object 420 rng <- .getRNGattribute(object) 421 if( !is.null(rng) ){ 422 on.exit() 423 return( nextRNG(rng, ...) ) 424 } 425 426 # only work for numeric seeds 427 if( !is.numeric(object) ) 428 stop("Invalid seed: expecting a numeric seed.") 429 430 # set RNG 431 .setRNG(object, ...) 432 433 # perform some draws 434 if( ndraw > 0 ){ 435 if( !isNumber(ndraw) ) 436 stop("Invalid value in argument `ndraw`: single numeric value expected.") 437 runif(ndraw) 438 } 439 # return new RNG settings 440 res <- RNGseed() 441 res 442} 443 444#' @importFrom utils head 445.collapse <- function(x, sep=', ', n=7L){ 446 447 res <- paste(head(x, n), collapse=', ') 448 if( length(x) > n ) 449 res <- paste(res, '...', sep=', ') 450 res 451} 452 453#' \code{setRNG} set the current RNG with a seed, 454#' using a suitable \code{.setRNG} method to set these settings. 455#' 456#' @param verbose a logical that indicates if the new RNG settings should 457#' be displayed. 458#' 459#' @param check logical that indicates if only valid RNG kinds should be 460#' accepted, or if invalid values should just throw a warning. 461#' Note that this argument is used only on R >= 3.0.2. 462#' 463#' @return \code{setRNG} invisibly returns the old RNG settings as 464#' they were before changing them. 465#' 466#' @export 467#' @rdname rng 468#' @examples 469#' 470#' obj <- list(x=1000, rng=123) 471#' setRNG(obj) 472#' rng <- getRNG() 473#' runif(10) 474#' set.seed(123) 475#' rng.equal(rng) 476#' 477setRNG <- function(object, ..., verbose=FALSE, check = TRUE){ 478 479 # do nothing if null 480 if( is.null(object) ) return() 481 482 # use RNG data from object if available 483 rng <- getRNG(object, ...) 484 if( !is.null(rng) && !identical(rng, object) ) return( setRNG(rng, ...) ) 485 486 # get/restore .Random.seed on.exit in case of errors 487 orseed <- getRNG() 488 on.exit({ 489 message("Restoring RNG settings probably due to an error in setRNG") 490 RNGseed(orseed) 491 }) 492 493 # call S4 method on object 494 # check validity of the seed 495 invalid_seed <- NULL 496 withCallingHandlers(.setRNG(object, ...) 497 , warning = function(err){ 498 invalid_seed <<- grepl("\\.Random\\.seed.* is not a valid", err$message) 499 if( check && testRversion('> 3.0.1') && invalid_seed ){ 500 stop("setRNG - Invalid RNG kind [", str_out(object), "]: " 501 , err$message, '.' 502 , call.=FALSE) 503 }else if( !invalid_seed ){ # if the seed is not invalid then we show the warning and continue 504 warning(err) 505 invokeRestart("muffleWarning") 506 507 } 508 } 509 ) 510 if( isTRUE(invalid_seed) ) return() 511 # cancel RNG restoration 512 on.exit() 513 if( verbose ) showRNG() 514 515 invisible(orseed) 516} 517 518#' Setting RNG Seeds 519#' 520#' \code{.setRNG} is an S4 generic that sets the current RNG settings, from a 521#' variety of specifications. 522#' Its methods define the workhorse functions that are called by \code{setRNG}. 523#' 524#' @inheritParams setRNG 525#' @export 526setGeneric('.setRNG', function(object, ...) standardGeneric('.setRNG') ) 527#' @describeIn .setRNG Sets the RNG to the kind(s) specified in \code{object}. 528#' If object is a single string that is a valid RNG kind, then this method is equivalent to \code{RNGkind(object, ...}. 529#' Otherwise, each element is assumed to be a valid argument for [RNGkind]. 530#' Note that this latter case the names of `object`, if any, are used as argument names in the call to [RNGkind]. 531#' 532#' @examples 533#' # set RNG kind 534#' old <- setRNG('Marsaglia') 535#' # restore 536#' setRNG(old) 537setMethod('.setRNG', 'character', 538 function(object, ...){ 539 if( length(object) == 1L ) 540 RNGkind(kind=object, ...) 541 else{ 542 n0 <- .RNGkind_length() 543 if( length(object) > n0 ){ 544 warning("RNG character specification is too long: discarding elements ", paste0(tail(object, -n0), collapse = ", ")) 545 546 } 547 do.call(RNGkind, as.list(head(object, n0))) 548 549 } 550 } 551) 552 553#' @describeIn .setRNG Sets the RNG settings using \code{object} directly the new value for 554#' \code{.Random.seed} or to initialise it with \code{\link{set.seed}}. 555#' 556#' @examples 557#' 558#' # directly set .Random.seed 559#' rng <- getRNG() 560#' r <- runif(10) 561#' setRNG(rng) 562#' rng.equal(rng) 563#' 564#' # initialise from a single number (<=> set.seed) 565#' setRNG(123) 566#' rng <- getRNG() 567#' runif(10) 568#' set.seed(123) 569#' rng.equal(rng) 570#' 571setMethod('.setRNG', 'numeric', 572 function(object, ...){ 573 574 if( length(object) == 1L ){ 575 if( is.integer(object) ){ # set kind and draw once to fix seed 576 RNGseed(object) 577 # check validity of the seed 578 tryCatch(runif(1) 579 , error = function(err){ 580 stop("setRNG - Invalid RNG kind [", object, "]: " 581 , err$message, '.' 582 , call.=FALSE) 583 } 584 ) 585 RNGseed() 586 }else{ # pass to set.seed 587 set.seed(object, ...) 588 } 589 }else{ 590 seed <- as.integer(object) 591 RNGseed(seed) 592 # check validity of the seed 593 tryCatch(runif(1) 594 , error=function(err){ 595 stop("setRNG - Invalid numeric seed [" 596 , .collapse(seed, n=5), "]: ", err$message, '.' 597 , call.=FALSE) 598 } 599 , warning = function(w){ 600 stop("setRNG - Invalid numeric seed [" 601 , .collapse(seed, n=5), "]: ", w$message, '.' 602 , call.=FALSE) 603 } 604 , finally = { 605 if( !identical(seed[1L], RNGseed()[1L]) ){ 606 msg <- "detected that the RNG kind would change after frist draw." 607 stop("setRNG - Invalid numeric seed [" 608 , .collapse(seed, n=5), "]: ", msg, '.' 609 , call.=FALSE) 610 } 611 }) 612 # re-force setting the seed if no error happened 613 RNGseed(seed) 614 } 615 } 616) 617 618#' \code{RNGdigest} computes a hash from the RNG settings associated with an 619#' object. 620#' 621#' @import digest 622#' @rdname RNGstr 623#' @export 624#' 625#' @examples 626#' # compute digest hash from RNG settings 627#' RNGdigest() 628#' RNGdigest(1234) 629#' # no validity check is performed 630#' RNGdigest(2:3) 631#' 632RNGdigest <- function(object=getRNG()){ 633 634 x <- object 635 object <- getRNG(x) 636 637 # exit if no RNG was extracted 638 if( is.null(object) ){ 639 warning("Found no embedded RNG data in object [", class(x),"]: returned NULL digest [", digest(NULL), '].') 640 return(digest(NULL)) # TODO: return NULL 641 } 642 643 digest(object) 644 645} 646 647#' Comparing RNG Settings 648#' 649#' \code{rng.equal} compares the RNG settings associated with two objects. 650#' 651#' These functions return \code{TRUE} if the RNG settings are identical, 652#' and \code{FALSE} otherwise. 653#' The comparison is made between the hashes returned by \code{RNGdigest}. 654#' 655#' @param x objects from which RNG settings are extracted 656#' @param y object from which RNG settings are extracted 657#' 658#' @return \code{rng.equal} and \code{rng.equal1} return a \code{TRUE} or 659#' \code{FALSE}. 660#' 661#' @rdname rngcmp 662#' @export 663rng.equal <- function(x, y){ 664 if( missing(y) ) 665 y <- getRNG() 666 identical(RNGdigest(x), RNGdigest(y)) 667} 668 669#' \code{rng1.equal} tests whether two objects have identical 670#' \strong{initial} RNG settings. 671#' 672#' @rdname rngcmp 673#' @export 674rng1.equal <- function(x, y){ 675 if( missing(y) ) 676 y <- getRNG() 677 rng.equal(getRNG1(x), getRNG1(y)) 678} 679