RiDMC/R/lexp.R

   1 lexp <- function(idmc_model, par, var, time, eps) {
   2   m <- idmc_model
   3   checkModelParVar(m, par, var)
   4   checkPositiveScalar(time)
   5   if(getModelType(m)=='C') {
   6     if(missing(eps)) {
   7       eps <- getOption('ts.eps')
   8       message('using eps = ', eps)
   9     }
  10     checkPositiveScalar(eps)
  11     ans <- .Call("ridmc_lexp_ode", m$model, as.double(par), as.double(var),
  12       as.double(time), as.double(eps), PACKAGE='RiDMC')
  13   } else {
  14     ans <- .Call("ridmc_lexp", m$model,
  15       as.double(par), as.double(var), as.integer(time), PACKAGE='RiDMC')
  16   }
  17   return(ans)
  18 }
  19
  20 LyapunovExponents <- function(idmc_model, par, var, time, eps,
  21   which.par.vary=1, par.min, par.max, par.howMany=100) {
  22   checkModelParVar(idmc_model, par, var)
  23   checkPositiveScalar(time)
  24   modelType <- getModelType(idmc_model)
  25   if(modelType=='C') {
  26     if(missing(eps)) {
  27       eps <- getOption('ts.eps')
  28       message('using eps = ', eps)
  29     }
  30     checkPositiveScalar(eps)
  31   }
  32   m <- idmc_model
  33   np <- par.howMany
  34   nv <- length(var)
  35   parValues <- seq(par.min, par.max, length=np)
  36   val <- matrix(,np,nv)
  37   if(modelType=='D') {
  38     for(i in seq_along(parValues)) {
  39       par[which.par.vary] <- parValues[i]
  40       val[i,] <- .Call("ridmc_lexp", m$model,
  41         as.double(par), as.double(var),
  42         as.integer(time), PACKAGE='RiDMC')
  43     }
  44   } else if(modelType=='C'){
  45     for(i in seq_along(parValues)) {
  46       par[which.par.vary] <- parValues[i]
  47       val[i,] <- .Call("ridmc_lexp_ode",
  48         m$model, as.double(par), as.double(var),
  49         as.double(time), as.double(eps), PACKAGE='RiDMC')
  50     }
  51   } else
  52     stop('invalid model type')
  53   ans <- list()
  54   ans$model <- idmc_model
  55   ans$var <- var
  56   ans$par <- par
  57   ans$values <- val
  58   ans$par.values <- parValues
  59   ans$which.par <- getModelParNames(m)[which.par.vary]
  60   class(ans) <- 'idmc_lexp_diagram'
  61   return(ans)
  62 }
  63 as.matrix.idmc_lexp_diagram <- function(x, ...)
  64   x$values
  65
  66 print.idmc_lexp_diagram <- function(x, ...) {
  67   m <- x$model
  68   cat('=iDMC Lyapunov exponents diagram=\n')
  69   cat('Model: ', getModelName(m), '\n')
  70   cat('Starting point: ')
  71     tmp <- getModelVarNames(m)
  72     cat(paste(tmp, x$var, sep=' = ', collapse=', '), '\n')
  73   cat('Parameter values: ')
  74     tmp <- getModelParNames(m)
  75     cat(paste(tmp, x$par, sep=' = ', collapse=', '), '\n')
  76   cat('Varying par.: ', x$which.par, '\n')
  77   cat('Varying par. range: [', paste(range(x$par.values), collapse=', '), ']\n')
  78   mle <- range(apply(x$values, 1, max, na.rm=TRUE), na.rm=TRUE)
  79   cat('MLE range: [', paste(format(mle, digits=4), collapse=', '), ']\n')
  80 }
  81
  82 as.grob.idmc_lexp_diagram <- function(x, col, lty, ylim = NULL, ...) {
  83   y <- x$values
  84   x1 <- x$par.values
  85   nc <- NCOL(y)
  86   if(missing(col))
  87     col <- seq_len(nc)
  88   if(missing(lty))
  89     lty <- rep(1,nc)
  90   xlim <- range(x1)
  91   if(is.null(ylim))
  92     ylim <- range(as.vector(y))
  93   lines <- list()
  94   for(j in seq_len(nc))
  95     lines[[j]] <- linesGrob(x1, y[,j], gp=gpar(col=col[j], lty=lty[j], ...), default.units = "native")
  96   lines <- do.call(gList, lines)
  97   cG <- gTree(children=lines, name='lines')
  98   contentsGrob(cG, xlim=xlim, ylim=ylim)
  99 }
 100
 101 plot.idmc_lexp_diagram <- function(x, y, col, lty,
 102   main = getModelName(x$model),
 103   xlab = x$which.par,
 104   ylab = 'Lyapunov exponent',
 105   ylim = NULL,
 106   mar = NULL,
 107   axes=TRUE,
 108   bty=TRUE,
 109   ...) {
 110   cG <- as.grob(x, col, lty, ylim = NULL, ...)
 111   pG <- plotGrob(cG, axes=axes, main=main, xlab=xlab, ylab=ylab, mar=mar)
 112   grid.draw(pG)
 113 }