Nothing
#' microsimulation
#'
#' Discrete event simulations in both R and C++ with Tools for Cost-Effectiveness Analysis.
#'
#' @section Introduction:
#'
#' Discrete event simulations in both R and C++ with Tools for Cost-Effectiveness Analysis.
#'
#' @name microsimulation-package
#' @aliases microsimulation
#' @author Mark Clements \email{mark.clements@ki.se}
#' @references \url{https://github.com/mclements/microsimulation}
#' @seealso \code{\link[Rcpp]{sourceCpp}}
#' @useDynLib microsimulation, .registration=TRUE
#' @import Rcpp
#' @import methods
#' @importFrom graphics lines plot
#' @importFrom stats predict rnorm sd
#' @importFrom ascii ascii
"_PACKAGE"
#' Cat a string for the library archive for use in loading the package
#'
#' @rdname Utilities
#' @return No return value, called for side effects
#' @export
LdFlags <- function()
cat(system.file("lib/libmicrosimulation.a", package="microsimulation", mustWork=TRUE))
#' On entry to the package, initialise the current stream in C.
#'
#' Is this function needed? We could define the current stream in open code.
#'
#' @param PACKAGE package from which this is called.
#'
#' @rdname Utilities
#' @return No return value, called for side effects
#' @export
microsimulation.init <- function(PACKAGE="microsimulation") {
.C("r_create_current_stream",PACKAGE=PACKAGE)
return(1)
}
#' On exit from the package, remove the current stream.
#'
#' Again, is this needed?
#'
#' @param PACKAGE package from which this is called.
#'
#' @return No return value, called for side effects
#' @rdname Utilities
#' @export
microsimulation.exit <- function(PACKAGE="microsimulation") {
.C("r_remove_current_stream",PACKAGE=PACKAGE)
return(1)
}
## http://r.789695.n4.nabble.com/How-to-construct-a-valid-seed-for-l-Ecuyer-s-method-with-given-Random-seed-td4656340.html
#' Convert from signed to unsigned
#'
#' @param seed signed seed (possibly a vector)
#' @return unsigned seed
#' @rdname Utilities
#' @export
unsigned <- function(seed) ifelse(seed < 0, seed + 2^32, seed)
#' Convert from unsigned to signed
#'
#' @param seed unsigned seed (possibly a vector)
#' @return signed seed
#'
#' @rdname Utilities
#' @export
signed <- function(seed) ifelse(seed>2^31, seed-2^32, seed)
#' Random draw for a positive (or other lower bound) random normal distribution
#'
#' @param n integer for the number of draws
#' @param mean numeric for the mean of the (untruncated) normal distribution (default=0)
#' @param sd numeric for the sd of the (untruncated) normal distribution (default=1)
#' @param lbound numeric for the lower bound (default=0)
#' @return numeric vector
#'
#' @importFrom stats rnorm sd
#' @rdname Utilities
#' @export
rnormPos <- function(n,mean=0,sd=1,lbound=0) {
if (length(mean)<n) mean <- rep(mean,length=n)
if (length(sd)<n) sd <- rep(sd,length=n)
x <- rnorm(n,mean,sd)
while(any(i <- which(x<lbound)))
x[i] <- rnorm(length(i),mean[i],sd[i])
x
}
#' Set the RngStream random number seed
#'
#' @param seed random number seed
#' @param PACKAGE package for the seed
#' @return invisibly returns the new seed
#' @rdname Utilities
#' @export
set.user.Random.seed <- function (seed,PACKAGE="microsimulation") {
seed <- as.double(unsigned(seed))
if (length(seed) == 1) seed <- rep(seed,6)
if (length(seed) == 7) seed <- seed[-1]
.C("r_set_user_random_seed",seed = seed,PACKAGE=PACKAGE)
return(invisible(seed))
}
#' Advance the RngStream random number seed
#'
#' @param seed random number seed
#' @param n number of sub-streams to advance
#' @param PACKAGE package for the seed
#' @return the advanced seed
#'
#' @rdname Utilities
#' @export
advance.substream <- function (seed,n,PACKAGE="microsimulation") {
seed <- as.double(unsigned(seed))
if (length(seed) == 1) seed <- rep(seed,6)
if (length(seed) == 7) seed <- seed[-1]
.C("r_rng_advance_substream", seed = seed, n = as.integer(n), PACKAGE=PACKAGE)$seed
}
#' Advance the RngStream random number seed to the next sub-stream
#'
#' @param seed random number seed
#' @param PACKAGE package for the seed
#' @return invisibly returns TRUE -- called for side effect
#'
#' @rdname Utilities
#' @export
next.user.Random.substream <- function(PACKAGE="microsimulation") {
.C("r_next_rng_substream", PACKAGE=PACKAGE)
return(invisible(TRUE))
}
#' Get the current RngStream random seed
#'
#' @param PACKAGE package for the seed
#' @return random seed
#' @rdname Utilities
#' @export
user.Random.seed <- function(PACKAGE="microsimulation") {
c(407L,
as.integer(signed(.C("r_get_user_random_seed", seed=as.double(rep(1,6)),
PACKAGE=PACKAGE)$seed)))
}
#' Add labels to a (default) zero-based integer to give a factor
#'
#' @param obj integer or logical for factor levels
#' @param labels labels for the factor levels
#' @param start first value of the levels
#' @return the new factor
#'
#' @rdname Utilities
#' @export
enum <- function(obj, labels, start=0) {
## stopifnot(is.logical(obj) || is.numeric(obj) || is.factor(obj))
obj <- as.integer(obj)
structure(factor(obj, levels=start + (0:(length(labels)-1)), labels=labels),start=start)
}
#' Assign labels to a (default) zero-based integer to give a factor
#'
#' @param obj integer or logical for factor levels
#' @param value labels for the factor levels
#' @return update the factor
#' @rdname Utilities
#' @export
"enum<-" <- function(obj, value) {
enum(if(is.factor(obj)) unclass(obj)-1+attr(obj,"start") else obj, value)
}
#' S3 class for maintain state for .Random.seed
#'
#' @return a list with oldseed (the old value of .Random.seed), and reset(), which resets .Random.seed
#'
#' @rdname Utilities
#' @export
RNGstate <- function() {
## house-keeping for random streams (see parallel::clusterSetRNGStream)
oldseed <- if (exists(".Random.seed", envir = .GlobalEnv,
inherits = FALSE))
get(".Random.seed", envir = .GlobalEnv, inherits = FALSE)
else NULL
oldkind <- RNGkind()
reset <- function() {
if (!is.null(oldseed))
assign(".Random.seed", oldseed, envir = .GlobalEnv)
else {
## clean up if we created a .Random.seed
if (exists(".Random.seed" ,envir = .GlobalEnv, inherits = FALSE))
rm(.Random.seed, envir = .GlobalEnv, inherits = FALSE)
do.call(RNGkind, as.list(oldkind))
}
}
list(oldseed = oldseed, reset = reset)
}
#' Utility to calculate the cost-efficiency frontier
#'
#' @param x vector of x coordinates
#' @param y vector y coordinates
#' @param concave logical for whether to calculate a concave frontier (default=TRUE)
#' @param convex logical for whether to calculate a convex frontier (default=NULL)
#' @return a list with components x and y for the frontier
#' @rdname Utilities
#' @export
frontier <- function (x, y, concave=TRUE, convex=NULL)
{
## check arguments
stopifnot(is.logical(concave),
is.null(convex) || is.logical(convex),
is.numeric(x),
is.numeric(y),
length(x) == length(y),
!any(is.na(x)),
!any(is.na(y)))
## Change concave if convex is defined
if (!is.null(convex))
concave <- !convex
if (concave) {
ichull <- grDevices::chull(cbind(x, y))
## case: length == 0
if ((n <- length(ichull)) == 0)
return(NULL)
## case: length == 1
if (n == 1)
return(list(x=x[ichull], y=y[ichull]))
## case: length > 1
## if min(x) value is not first in ichull, then re-order
if ((iminx <- which.min(x[ichull])) > 1) {
ichull <- ichull[c(iminx:n, 1:(iminx-1))]
}
## find those that are increasing for x and y
include <- c(TRUE,diff(x[ichull])>0 & diff(y[ichull])>0)
list(x=x[ichull][include], y=y[ichull][include])
} else # convex case
with(Recall(y,x,concave=TRUE), list(x=y, y=x))
}
#' plot lines for a frontier
#'
#' @param x vector of x coordinates
#' @param y vector of y coordinates
#' @param pch type of pch for the plotted symbols (default=19)
#' @param type join type (default="b")
#' @param ... other arguments to lines
#'
#' @return No return value, called for side effects
#'
#' @rdname Utilities
lines_frontier <- function(x,y,pch=19,type="b",...) {
index <- frontier(x,y)
lines(x[index],y[index],pch=pch,type=type,...)
}
#' Integrate a discounted value
#'
#' @param y the undiscounted value
#' @param start the start time
#' @param finish the finish time
#' @param dr discount rate, expressed as a percentage
#' @return numeric discounted value
#' @export
discountedInterval <- function(y, start, finish, dr) {
duration = finish - start
ifelse(dr<=0 | duration==0, y*duration,
y / (1 + dr)^ start / log(1 + dr) *
(1 - 1 / (1 + dr)^ duration))
}
#' Discounted value
#'
#' @param y the undiscounted value
#' @param time the time of the event
#' @param dr discount rate, expressed as a percentage
#' @return numeric vector
#' @rdname Utilities
#' @export
discountedPoint <- function(y, time, dr)
ifelse(dr<=0 | time==0, y, y / (1 + dr) ^ time)
#' Call the Person simulation from C++
#'
#' Example that uses the RngStream random number generator
#'
#' @param n number of simulations (default=20)
#' @param seed random number seed
#' @return data-frame
#' @rdname Examples
#' @export
callPersonSimulation <- function(n=20,seed=rep(12345,6)) {
state <- RNGstate(); on.exit(state$reset())
RNGkind("user")
set.user.Random.seed(seed)
stateT =c("Healthy","Localised","DxLocalised","LocallyAdvanced",
"DxLocallyAdvanced","Metastatic","DxMetastatic","Death")
eventT =c("toDeath", "toPCDeath", "toLocalised", "toDxLocalised",
"toDxLocallyAdvanced",
"toLocallyAdvanced", "toMetastatic", "toDxMetastatic")
out <- .Call(.callPersonSimulation,
as.integer(rep(seed,length=6)), # magic
list(n=as.integer(n)))
out$state = enum(out$state,stateT)
out$event = enum(out$event,eventT)
as.data.frame(out)
}
#' Call the SimplePerson simulation from C++
#'
#' Example that uses the Mersenne-Twister random number generator
#'
#' @param n number of simulations (default=10)
#' @return data-frame
#' @rdname Examples
#' @export
callSimplePerson <- function(n=10) {
state <- RNGstate(); on.exit(state$reset())
RNGkind("Mersenne-Twister")
set.seed(12345)
stateT <- c("Healthy","Cancer","Death")
eventT <- c("toOtherDeath", "toCancer", "toCancerDeath")
out <- .Call(.callSimplePerson,
parms=list(n=as.integer(n)))
out$state = enum(out$state,stateT)
out$event = enum(out$event,eventT)
as.data.frame(out)
}
#' Call the SimplePerson2 simulation from C++
#'
#' Example that uses the Mersenne-Twister random number generator
#'
#' @param n number of simulations (default=10)
#' @return data-frame
#' @rdname Examples
#' @export
callSimplePerson2 <- function(n=10) {
state <- RNGstate(); on.exit(state$reset())
RNGkind("Mersenne-Twister")
set.seed(12345)
stateT <- c("Healthy","Cancer","Death")
eventT <- c("toOtherDeath", "toCancer", "toCancerDeath")
out <- .Call(.callSimplePerson2,
parms=list(n=as.integer(n)))
reader <- function(obj)
cbind(data.frame(state=enum(obj[[1]],stateT)),
data.frame(obj[-1]))
out <- lapply(out,reader)
out$events <- with(out$events, data.frame(state=state,event=enum(event,eventT),age=age,number=number))
out$pt <- with(out$pt, data.frame(state=state,age=age,pt=pt))
out$prev <- with(out$prev, data.frame(state=state,age=age,number=number))
out
}
#' Call the IllnessDeath simulation from C++
#'
#' Example that uses the Mersenne-Twister random number generator
#'
#' @param n number of simulations (default=10)
#' @param cure probability of cure
#' @param zsd frailty standard deviation
#' @return data-frame
#' @rdname Examples
#' @export
callIllnessDeath <- function(n=10L,cure=0.1,zsd=0) {
state <- RNGstate(); on.exit(state$reset())
RNGkind("Mersenne-Twister")
set.seed(12345)
stateT <- c("Healthy","Cancer")
eventT <- c("toOtherDeath", "toCancer", "toCancerDeath")
out <- .Call(.callIllnessDeath,
parms=list(n=as.integer(n),cure=as.double(cure),zsd=as.double(zsd)))
reader <- function(obj)
cbind(data.frame(state=enum(obj[[1]],stateT)),
data.frame(obj[-1]))
out <- lapply(out,reader)
out$events <- with(out$events, data.frame(state=state,event=enum(event,eventT),age=age,number=number))
out$pt <- with(out$pt, data.frame(state=state,age=age,pt=pt))
out$prev <- with(out$prev, data.frame(state=state,age=age,number=number))
out
}
#' Calculate incremental cost-effectiveness ratios from two objects
#'
#' @param object1 first object
#' @param object2 second object
#' @param ... other arguments
#' @rdname Utilities
#' @export
ICER <- function(object1, object2, ...)
UseMethod("ICER")
#' Reference class implementation of an event queue in R
#'
#' This event queue is simple and useful for pedagogic purposes.
#'
#' The algorithm for pushing values into the queue is computationally
#' very simple: simply rank the times using \code{order()} and re-order
#' times and events. This approach is probably of acceptable performance
#' for smaller queue. A more computationally efficient approach for
#' pushing into larger queues would be to use a binary search (e.g. using
#' \code{findInterval()}).
#'
#' For faster alternatives, see \code{pqueue} and \code{PQueueRef}.
#'
#' @examples
#' pq = new("EventQueue")
#' pq$push(3,"Clear drains")
#' pq$push(4, "Feed cat")
#' pq$push(5, "Make tea")
#' pq$push(1, "Solve RC tasks")
#' pq$push(2, "Tax return")
#' while(!pq$empty())
#' print(pq$pop())
#'
#' @import methods
#' @exportClass EventQueue
#' @field times vector of times
#' @field events list of events
#' @rdname Classes
EventQueue <-
setRefClass("EventQueue",
fields = list(times = "numeric", events = "list"),
methods = list(
help = function() {
'Reference class implementation of an event queue. Fields for the event times and the events list.'
},
push = function(time,event) {
'Method to insert the event at the given time'
insert.ord <- findInterval(time,times)
times <<- append(times,time,insert.ord)
events <<- append(events,list(event),insert.ord)
},
pop = function() {
'Method to remove the head of the event queue and return its value'
head <- list(event=events[[1]], time=times[1])
times <<- times[-1]
events <<- events[-1]
return(head)
},
empty = function() {
'Method to check whether there are no events in the queue'
length(times) == 0
},
clear = function() {
'Method to clear the event queue'
times <<- numeric()
events <<- list()
},
cancel = function(predicate, ...) {
'Method to remove events that satisfy some predicate'
if (!empty()) {
i <- sapply(events, predicate, ...)
stopifnot(is.logical(i))
i[is.na(i)] <- TRUE
times <<- times[!i]
events <<- events[!i]
}
}))
#' Reference class implementation of a discrete event simulation
#'
#' Inherit from this class to represent a discrete event simulation. The
#' API is similar to that for Omnet++, where an \code{init} method sets up
#' the initial events using the \code{scheduleAt(time,event)} method, the
#' messages are handled using the \code{handleMessage(event)} method, the
#' simulation is run using the \code{run} method, and the \code{final}
#' method is called at the end of the simulation.
#'
#' @examples
#' DES = setRefClass("DES",
#' contains = "BaseDiscreteEventSimulation",
#' methods=list(
#' init=function() {
#' scheduleAt(3,"Clear drains")
#' scheduleAt(4, "Feed cat")
#' scheduleAt(5, "Make tea")
#' scheduleAt(1, "Solve RC tasks")
#' scheduleAt(2, "Tax return")
#' },
#' handleMessage=function(event) print(event)))
#'
#' des = new("DES")
#' des$run()
#' @examples
#' \dontrun{
#' testRsimulation1 <- function() {
#' ## A simple example
#' Simulation <-
#' setRefClass("Simulation",
#' contains = "BaseDiscreteEventSimulation")
#' Simulation$methods(
#' init = function() {
#' scheduleAt(rweibull(1,8,85), "Death due to other causes")
#' scheduleAt(rweibull(1,3,90), "Cancer diagnosis")
#' },
#' handleMessage = function(event) {
#' if (event %in% c("Death due to other causes", "Cancer death")) {
#' clear()
#' print(event)
#' }
#' else if (event == "Cancer diagnosis") {
#' if (runif(1) < 0.5)
#' scheduleAt(now() + rweibull(1,2,10), "Cancer death")
#' print(event)
#' }
#' })
#' Simulation$new()$run()
#' }
#'
#' ## An extension with individual life histories
#' testRsimulation2 <- function(n=100) {
#' Simulation <-
#' setRefClass("Simulation",
#' contains = "BaseDiscreteEventSimulation",
#' fields = list(state = "character", report = "data.frame"))
#' Simulation$methods(
#' init = function() {
#' report <<- data.frame()
#' state <<- "Healthy"
#' scheduleAt(rweibull(1,8,85), "Death due to other causes")
#' scheduleAt(rweibull(1,3,90), "Cancer diagnosis")
#' },
#' handleMessage = function(event) {
#' report <<- rbind(report, data.frame(state = state,
#' begin = attr(event,"sendingTime"),
#' end = currentTime,
#' event = event,
#' stringsAsFactors = FALSE))
#' if (event %in% c("Death due to other causes", "Cancer death")) {
#' clear()
#' }
#' else if (event == "Cancer diagnosis") {
#' state <<- "Cancer"
#' if (runif(1) < 0.5)
#' scheduleAt(now() + rweibull(1,2,10), "Cancer death")
#' }
#' },
#' final = function() report)
#' sim <- Simulation$new()
#' do.call("rbind", lapply(1:n, function(id) data.frame(id=id,sim$run())))
#' }
#'
#' ## reversible illness-death model
#' testRsimulation3 <- function(n=100) {
#' Simulation <-
#' setRefClass("Simulation",
#' contains = "BaseDiscreteEventSimulation",
#' fields = list(state = "character", everCancer = "logical",
#' report = "data.frame"))
#' Simulation$methods(
#' init = function() {
#' report <<- data.frame()
#' state <<- "Healthy"
#' everCancer <<- FALSE
#' scheduleAt(rweibull(1,8,85), "Death due to other causes")
#' scheduleAt(rweibull(1,3,90), "Cancer diagnosis")
#' },
#' handleMessage = function(event) {
#' report <<- rbind(report, data.frame(state = state,
#' everCancer = everCancer,
#' begin = attr(event,"sendingTime"),
#' end = currentTime,
#' event = event,
#' stringsAsFactors = FALSE))
#' if (event %in% c("Death due to other causes", "Cancer death")) {
#' clear()
#' }
#' else if (event == "Cancer diagnosis") {
#' state <<- "Cancer"
#' everCancer <<- TRUE
#' if (runif(1) < 0.5)
#' scheduleAt(now() + rweibull(1,2,10), "Cancer death")
#' scheduleAt(now() + 10, "Recovery")
#' }
#' else if (event == "Recovery") {
#' state <<- "Healthy"
#' scheduleAt(now() + rexp(1,10), "Cancer diagnosis")
#' }
#' },
#' final = function() report)
#' sim <- Simulation$new()
#' do.call("rbind", lapply(1:n, function(id) data.frame(id=id,sim$run())))
#' }
#'
#' ## cancer screening
#' testRsimulation4 <- function(n=1) {
#' Simulation <-
#' setRefClass("Simulation",
#' contains = "BaseDiscreteEventSimulation",
#' fields = list(state = "character", report = "data.frame"))
#' Simulation$methods(
#' init = function() {
#' report <<- data.frame()
#' state <<- "Healthy"
#' scheduleAt(rweibull(1,8,85), "Death due to other causes")
#' scheduleAt(rweibull(1,3,90), "Cancer onset")
#' scheduleAt(50,"Screening")
#' },
#' handleMessage = function(event) {
#' report <<- rbind(report, data.frame(state = state,
#' begin = attr(event,"sendingTime"),
#' end = currentTime,
#' event = event,
#' stringsAsFactors = FALSE))
#' if (event %in% c("Death due to other causes", "Cancer death")) {
#' clear()
#' }
#' else if (event == "Cancer onset") {
#' state <<- event
#' dx <- now() + rweibull(1,2,10)
#' scheduleAt(dx, "Clinical cancer diagnosis")
#' scheduleAt(dx + rweibull(1,1,10), "Cancer death")
#' scheduleAt(now() + rweibull(1,1,10), "Metastatic cancer")
#' }
#' else if (event == "Metastatic cancer") {
#' state <<- event
#' cancel(function(event) event %in%
#' c("Clinical cancer diagnosis","Cancer death")) # competing events
#' scheduleAt(now() + rweibull(1,2,5), "Cancer death")
#' }
#' else if (event == "Clinical cancer diagnosis") {
#' state <<- event
#' cancel(function(event) event == "Metastatic cancer")
#' }
#' else if (event == "Screening") {
#' switch(state,
#' "Cancer onset" = {
#' state <<- "Screen-detected cancer diagnosis"
#' cancel(function(event) event %in%
#' c("Clinical cancer diagnosis","Metastatic cancer"))
#' },
#' "Metastatic cancer" = {}, # ignore
#' "Clincal cancer diagnosis" = {}, # ignore
#' "Healthy" = {
#' if (now()<=68) scheduleAt(now()+2, "Screening")
#' })
#' }
#' else stop(event)
#' },
#' final = function() report)
#' sim <- Simulation$new()
#' do.call("rbind", lapply(1:n, function(id) data.frame(id=id,sim$run())))
#' }
#'
#' ## ticking bomb - toy example
#' testRsimulation5 <- function(n=1) {
#' Simulation <-
#' setRefClass("Simulation",
#' contains = "BaseDiscreteEventSimulation",
#' fields = list(report = "data.frame"))
#' Simulation$methods(
#' init = function() {
#' report <<- data.frame()
#' scheduleAt(rexp(1,1), "tick")
#' if (runif(1)<0.1)
#' scheduleAt(rexp(1,1), "explosion")
#' },
#' handleMessage = function(event) {
#' report <<- rbind(report, data.frame(begin = attr(event,"sendingTime"),
#' end = currentTime,
#' event = event,
#' stringsAsFactors = FALSE))
#' if (event == "explosion")
#' clear()
#' else {
#' clear() # queue
#' if (event == "tick") scheduleAt(currentTime+rexp(1,1), "tock")
#' else scheduleAt(currentTime+rexp(1,1), "tick")
#' if (runif(1)<0.1)
#' scheduleAt(currentTime+rexp(1,1), "explosion")
#' }
#' },
#' final = function() report)
#' sim <- Simulation$new()
#' do.call("rbind", lapply(1:n, function(id) data.frame(id=id,sim$run())))
#' }
#' }
#'
#' @import methods
#' @exportClass EventQueue
#' @field times vector of times
#' @field events list of events
#' @rdname Classes
BaseDiscreteEventSimulation <-
setRefClass("BaseDiscreteEventSimulation",
contains = "EventQueue",
## contains = "PQueueRef",
fields = list(currentTime = "numeric",
previousEventTime = "numeric"),
methods = list(
help = function() {
'Reference class implementation of an event-oriented discrete event simulation. Fields for the event times and the events list.'
},
scheduleAt = function(time, event) {
'Method that adds attributes for the event time and the sendingTime, and then insert the event into the event queue'
attr(event,"time") <- time
attr(event,"sendingTime") <- currentTime
push(time, event)
},
init = function() {
'Virtual method to initialise the event queue and attributes'
stop("VIRTUAL!")
},
handleMessage = function(event) {
'Virtual method to handle the messages as they arrive'
stop("VIRTUAL!")
},
final = function() {
'Method for finalising the simulation'
NULL
},
now = function() currentTime,
previous = function() previousEventTime,
reset = function(startTime = 0.0) {
'Method to reset the event queue'
clear()
previousEventTime <<- currentTime <<- startTime
},
run = function(startTime = 0.0) {
'Method to run the simulation'
reset(startTime)
init()
while(!empty()) {
event <- pop()$event
currentTime <<- attr(event, "time")
handleMessage(event)
previousEventTime <<- currentTime
}
final()
}))
#' S3 class to work with RngStream objects
#'
#' @param nextStream whether to move to the next stream (default=TRUE)
#' @param iseed set seed after changing RNG (otherwise keep the current seed)
#' @return list of class \code{RNGStream} with components:
#' \describe{
#' \item{resetRNGkind}{function to reset to the previous RNG and seed}
#' \item{seed}{function to return the current seed}
#' \item{open}{function to use the current seed}
#' \item{close}{function to make the current seed equal to .Random.seed}
#' \item{resetStream}{function to move back to start of stream}
#' \item{resetSubStream}{function to move back to start of sub-stream}
#' \item{nextSubStream}{function to move to next sub-stream}
#' \item{nextStream}{function to move to next stream}
#' }
#' @examples
#' ## set up one stream
#' s1 <- RNGStream()
#' s1$open()
#' rnorm(1)
#' s1$nextSubStream()
#' rnorm(1)
#' ## reset the stream
#' s1$resetStream()
#' rnorm(2)
#' s1$nextSubStream()
#' rnorm(2)
#'
#' ## now do with two streams
#' s1$resetStream()
#' s2 <- RNGStream()
#' with(s1,rnorm(1))
#' with(s2,rnorm(1))
#' s1$nextSubStream()
#' with(s1,rnorm(1))
#' ## now reset the streams and take two samples each time
#' s1$resetStream()
#' s2$resetStream()
#' with(s1,rnorm(2))
#' with(s2,rnorm(2))
#' s1$nextSubStream()
#' with(s1,rnorm(2))
#' @rdname RNGStream
#' @export
RNGStream <- function(nextStream = TRUE, iseed = NULL) {
stopifnot(exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE) || !is.null(iseed))
oldseed <- if (exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE))
get(".Random.seed", envir = .GlobalEnv, inherits = FALSE)
else NULL
if (RNGkind()[1] != "L'Ecuyer-CMRG") RNGkind("L'Ecuyer-CMRG")
if (!is.null(iseed)) set.seed(iseed)
current <- if (nextStream) parallel::nextRNGStream(.Random.seed) else .Random.seed
.Random.seed <<- startOfStream <- startOfSubStream <- current
structure(list(resetRNGkind = function()
{
if (!is.null(oldseed))
assign(".Random.seed", oldseed, envir = .GlobalEnv)
else rm(.Random.seed, envir = .GlobalEnv)
},
seed = function() current,
open = function() .Random.seed <<- current,
close = function() current <<- .Random.seed,
resetStream = function() .Random.seed <<- current <<- startOfSubStream <<- startOfStream,
resetSubStream = function() .Random.seed <<- current <<- startOfSubStream,
nextSubStream = function() .Random.seed <<- current <<- startOfSubStream <<- parallel::nextRNGSubStream(startOfSubStream),
nextStream = function() .Random.seed <<- current <<- startOfSubStream <<- startOfStream <<- parallel::nextRNGStream(startOfStream)),
class="RNGStream")
}
#' Use RNGStream as an old class
#'
#' @name RNGStream-class
#' @aliases RNGStream
#' @rdname RNGStream
#' @exportClass RNGStream
setOldClass("RNGStream")
#' With method for RNGStream S3 class
#'
#' @param data object of type RNGStream
#' @param expr expression using the RNGStream
#' @param ... other arguments passed to eval()
#' @return the value from the expression
#' @rdname RNGStream
with.RNGStream <- function(data,expr,...) {
data$open()
out <- eval(substitute(expr), enclos = parent.frame(), ...)
data$close()
out
}
#' Old data used in the prostata model
#'@rdname Data
"fhcrcData"
#' C++ function
#' @rdname Internal
#' @name callCalibrationSimulation
#' @return data-frame
NULL
#' C++ function
#' @rdname Internal
#' @name r_create_current_stream
#' @return No return value, called for side effects
#' @export
NULL
#' C++ function
#' @rdname Internal
#' @name r_remove_current_stream
#' @return No return value, called for side effects
#' @export
NULL
#' C++ function
#' @rdname Internal
#' @name r_set_user_random_seed
#' @return No return value, called for side effects
#' @export
NULL
#' C++ function
#' @rdname Internal
#' @name r_rng_advance_substream
#' @return No return value, called for side effects
#' @export
NULL
#' C++ function
#' @rdname Internal
#' @name r_next_rng_substream
#' @return No return value, called for side effects
#' @export
NULL
#' C++ function
#' @rdname Internal
#' @name r_get_user_random_seed
#' @return No return value, called for side effects
#' @export
NULL
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.