RiDMC/R/trajectory.R

   1 Trajectory <- function(idmc_model, par, var, time=1, transient=0,
   2   seed, eps=getOption("ts.eps"), integrator=2) {
   3   m <- idmc_model
   4   checkModelParVar(m, par, var)
   5   ans <- list()
   6   ans$transient <- transient
   7   ans$time <- time
   8   ans$par <- par
   9   ans$var <- var
  10   ans$model <- m
  11   if(getModelType(m)=='C') { ##Continuous model
  12     integrator <- as.integer(integrator)
  13     if((integrator<0)||(integrator>8))
  14       stop('\'integrator\' should be an integer code between 0 and 8')
  15     checkPositiveScalar(eps)
  16     ans$eps <- eps
  17     ans$integrator <- integrator
  18     trajectory <- .Call("ridmc_ctrajectory_alloc",
  19       m$model, as.double(par), as.double(var),
  20       as.double(eps), as.integer(integrator), PACKAGE='RiDMC')
  21     ans$trajectory <- trajectory
  22     ans$step <- function()
  23       .Call("ridmc_ctrajectory_step", trajectory, PACKAGE='RiDMC')
  24     ans$getValue <- function()
  25       .Call("ridmc_ctrajectory_getValue", trajectory, PACKAGE='RiDMC')
  26     ans$getModel <- function() {
  27       pp <- .Call("ridmc_ctrajectory_getModel", trajectory, PACKAGE='RiDMC')
  28       buildModel(pp, getModelText(m))
  29     }
  30     class(ans) <- c("idmc_ctrajectory","idmc_trajectory")
  31   } else { ##Discrete model
  32     eps <- 1
  33     trajectory <- .Call("ridmc_dtrajectory_alloc",
  34       m$model, as.double(par), as.double(var), PACKAGE='RiDMC')
  35     ans$eps <- eps
  36     ans$trajectory <- trajectory
  37     ans$step <- function()
  38       .Call("ridmc_dtrajectory_step", trajectory, PACKAGE='RiDMC')
  39     ans$getValue <- function()
  40       .Call("ridmc_dtrajectory_getValue", trajectory, PACKAGE='RiDMC')
  41     ans$getModel <- function() {
  42       pp <- .Call("ridmc_dtrajectory_getModel", trajectory, PACKAGE='RiDMC')
  43       buildModel(pp, m$text)
  44     }
  45     class(ans) <- c("idmc_dtrajectory","idmc_trajectory")
  46   }
  47   vnames <- getModelVarNames(m)
  48   values <- matrix(var, 1, length(vnames))
  49   colnames(values) <- vnames
  50   ans$values <- values
  51   if(!missing(seed))
  52     setTrajectorySeed(ans, seed)
  53   if(transient>0) {
  54     ans <- stepTrajectory(ans, transient)
  55     ans$values <- ans$values[NROW(ans$values),,drop=FALSE]
  56   }
  57   stepTrajectory(ans, time)
  58 }
  59 getTrajectoryModel <- function(idmc_trajectory)
  60   idmc_trajectory$getModel()
  61 getTrajectoryValues <- function(idmc_trajectory)
  62   idmc_trajectory$values
  63 setTrajectorySeed <- function(idmc_trajectory, seed)
  64   getTrajectoryModel(idmc_trajectory)$set.seed(seed)
  65 as.matrix.idmc_trajectory <- function(x, ...)
  66   getTrajectoryValues(x)
  67 as.ts.idmc_trajectory <- function(x, ...)
  68   ts(as.matrix(x), frequency = 1/x$eps,
  69     start=x$transient, ...)
  70
  71 stepTrajectory <- function(idmc_trajectory, time=1) {
  72   tr <- idmc_trajectory
  73   eps <- tr$eps
  74   nsteps <- floor(time/eps)
  75   vars <- getModelVarNames(getTrajectoryModel(tr))
  76   oldTr <- getTrajectoryValues(idmc_trajectory)
  77   newTr <- matrix(,nsteps, length(vars))
  78   for(i in seq_len(nsteps)) {
  79     tr$step()
  80     newTr[i,] <- tr$getValue()
  81   }
  82   tr$values <- rbind(oldTr, newTr)
  83   tr
  84 }
  85
  86 print.idmc_ctrajectory <- function(x, ...) {
  87   modelInfo <- x$model$infos
  88   cat('= iDMC model continuous trajectory =\n')
  89   cat('model: ', getModelName(getTrajectoryModel(x)), '\n')
  90   cat('parameter values: ', paste(x$par, sep=','), '\n')
  91   cat('starting point: ', paste(x$var, sep=','),'\n')
  92   cat('transient length: ', x$transient, '\n')
  93   cat('time span: ', x$time, '\n')
  94   cat('step size: ', x$eps, '\n')
  95   cat('step function: ', x$integrator, '\n')
  96 }
  97 print.idmc_dtrajectory <- function(x, ...) {
  98   cat('= iDMC model discrete trajectory =\n')
  99   cat('model: ', getModelName(getTrajectoryModel(x)), '\n')
 100   cat('parameter values: ', paste(x$par, sep=','), '\n')
 101   cat('starting point: ', paste(x$var, sep=','),'\n')
 102   cat('transient length: ', x$transient, '\n')
 103   cat('time span: ', x$time, '\n')
 104 }
 105
 106 as.grob.idmc_trajectory <- function(x, vars=1:2, type='l', ...) {
 107   mdl <- getTrajectoryModel(x)
 108   varNames <- getModelVarNames(mdl)
 109   names(varNames) <- varNames
 110   vars <- vars[1:2]
 111   vars <- varNames[vars]
 112   if(length(varNames)<2) {
 113     y <- as.ts(x)
 114     xyGrob(time(y), y, type=type, name='xy', ...)
 115   } else {
 116     xx <- as.matrix(x)[,vars]
 117     xyGrob(xx[,1], xx[,2], type=type, name='xy', ...)
 118   }
 119 }
 120
 121 plot.idmc_trajectory <- function(x, y, vars=1:2, type='l',
 122   main = getModelName(getTrajectoryModel(x)), xlab, ylab,
 123   mar = NULL, axes=TRUE, bty=TRUE, ...) {
 124   mdl <- getTrajectoryModel(x)
 125   varNames <- getModelVarNames(mdl)
 126   names(varNames) <- varNames
 127   vars <- vars[1:2]
 128   vars <- varNames[vars]
 129   if(missing(xlab))
 130     xlab <- vars[1]
 131   if(missing(ylab))
 132     ylab <- vars[2]
 133   if(length(varNames)<2) {
 134     ylab <- varNames
 135     xlab <- 'time'
 136   }
 137   cG <- as.grob(x, vars=vars, type=type, ...)
 138   pG <- plotGrob(cG, axes=axes, main=main, xlab=xlab, ylab=ylab, mar=mar, bty=bty)
 139   grid.draw(pG)
 140 }
 141
 142 trajectoryList <- function(idmc_model, n=2, par, var, time=1, transient=0,
 143   seed, eps=getOption("ts.eps"), integrator=2) {
 144   argList <- expandArgList(n=n, par=par, var=var, time=time, transient=transient, eps=eps, integrator=integrator)
 145   argList <- lapply(argList, function(x) cons(idmc_model=idmc_model, tail=x))
 146   if(!missing(seed))
 147     argList <- lapply(argList, function(x) cons(seed=seed, tail=x))
 148   trList <- lapply(argList, do.call, what=Trajectory)
 149   class(trList) <- 'idmc_trajectoryList'
 150   trList
 151 }
 152
 153 as.grob.idmc_trajectoryList <- function(x, colors, ...) {
 154   as.grob2 <- function(obj, color, ...)
 155     as.grob(obj, gp=gpar(col=color), ...)
 156   childs <- mapply(as.grob2, x, colors, ..., SIMPLIFY=FALSE)
 157   xmin <- min(sapply(childs, function(x) min(Xlim(x))))
 158   xmax <- min(sapply(childs, function(x) max(Xlim(x))))
 159   ymin <- min(sapply(childs, function(x) min(Ylim(x))))
 160   ymax <- min(sapply(childs, function(x) max(Ylim(x))))
 161   childs <- do.call(gList, childs)
 162   contentsGrob(gTree(children=childs), xlim=c(xmin, xmax), ylim=c(ymin, ymax), respect=FALSE)
 163 }
 164
 165 plot.idmc_trajectoryList <- function(x, y, vars=1:2, type='l', colors,
 166   main = getModelName(getTrajectoryModel(x[[1]])), xlab, ylab,
 167   mar = NULL, axes=TRUE, bty=TRUE, ...) {
 168   if(missing(colors))
 169     colors <- seq_along(x)
 170   mdl <- getTrajectoryModel(x[[1]])
 171   varNames <- getModelVarNames(mdl)
 172   names(varNames) <- varNames
 173   vars <- vars[1:2]
 174   vars <- varNames[vars]
 175   if(missing(xlab))
 176     xlab <- vars[1]
 177   if(missing(ylab))
 178     ylab <- vars[2]
 179   if(length(varNames)<2) {
 180     ylab <- varNames
 181     xlab <- 'time'
 182   }
 183   cG <- as.grob(x, colors=as.list(colors), vars=vars, type=type, ...)
 184   grid.draw(plotGrob(cG, main=main, xlab=xlab, ylab=ylab, axes=axes, mar=mar, bty=bty))
 185 }