RiDMC/R/lexpMap.R

   1 LyapunovExponentsMap <- function(idmc_model, par, var, time, eps,
   2   par.x = 1, par.x.range, par.x.howMany=100, par.y = 2, par.y.range, par.y.howMany=100,
   3   eps.zero=sqrt(.Machine$double.eps)) {
   4   checkModelParVar(idmc_model, par, var)
   5   checkPositiveScalar(time)
   6   modelType <- getModelType(idmc_model)
   7   var <- as.double(var)
   8   nvar <- length(var)
   9   ##Lookup matrix:
  10   lm <- enumerateExponents(nvar)
  11   m <- idmc_model$model
  12   recodeExp <- function(x) {
  13     codes <- integer(nvar)
  14     codes[is.finite(x) & x>0] <- 3
  15     codes[is.finite(x) & x<0] <- 1
  16     codes[is.finite(x) & abs(x)<eps.zero] <- 2
  17     codes[!is.finite(x)] <- 4
  18     cd <- c(sum(codes==1), sum(codes==2), sum(codes==3), sum(codes==4))
  19     which(apply(lm, 1, function(x) all(x == cd)))
  20   }
  21   if(modelType=='C') {
  22     if(missing(eps)) {
  23       eps <- getOption('ts.eps')
  24       message('using eps = ', eps)
  25     }
  26     checkPositiveScalar(eps)
  27     eps <- as.double(eps)
  28     time <- as.double(time)
  29     lexpLocal <- function(par)
  30       recodeExp( .Call("ridmc_lexp_ode", m, as.double(par), var, time, eps, PACKAGE='RiDMC') )
  31   } else {
  32     time <- as.integer(time)
  33     lexpLocal <- function(par)
  34       recodeExp( .Call("ridmc_lexp", m, as.double(par), var, time, PACKAGE='RiDMC') )
  35   }
  36
  37   f <- function(a, b) {
  38     par[par.x] <- a
  39     par[par.y] <- b
  40     lexpLocal(par)
  41   }
  42
  43   p1seq <- seq(par.x.range[1], par.x.range[2], len=par.x.howMany)
  44   p2seq <- seq(par.y.range[1], par.y.range[2], len=par.y.howMany)
  45   mat <- outer(p1seq, p2seq, Vectorize(f))
  46   ans <- list()
  47   ans$model <- idmc_model
  48   ans$var <- var
  49   ans$par <- par
  50   ans$par.x <- par.x
  51   ans$par.x.range <- par.x.range
  52   ans$par.x.howMany <- par.x.howMany
  53   ans$par.y <- par.y
  54   ans$par.y.range <- par.y.range
  55   ans$par.y.howMany <- par.y.howMany
  56   ans$mat <- mat
  57   ans$labels <- rownames(lm)
  58   class(ans) <- 'idmc_lexp_map'
  59   ans
  60 }
  61
  62 enumerateExponents <- function(nvar) {
  63   size <- 4^nvar ##each exponent can be <0, =0, >0, non-finite
  64   permutations <- function(a) {
  65     a <- as.matrix(a)
  66     id <- rownames(a) <- seq_len(NROW(a))
  67     tmp <- expand.grid(seq_along(id), 1:4)
  68     cbind(a[tmp[,1],], tmp[,2])
  69   }
  70
  71   combine <- function(nvar, current=1:4) {
  72     if(nvar==1)
  73       current
  74     else
  75       combine(nvar-1, permutations(current))
  76   }
  77
  78   ##compute all permutations:
  79   ans2 <- combine(nvar)
  80   ans2 <- as.matrix(ans2)
  81   ##set row labels:
  82   howManyPos <- function(x) sum(x==3)
  83   howManyNeg <- function(x) sum(x==1)
  84   howManyNull <- function(x) sum(x==2)
  85   howManyDiverging <- function(x) sum(x==4)
  86   ans <- matrix(,NROW(ans2), 4)
  87   ans[,1] <- apply(ans2, 1, howManyNeg)
  88   ans[,2] <- apply(ans2, 1, howManyNull)
  89   ans[,3] <- apply(ans2, 1, howManyPos)
  90   ans[,4] <- apply(ans2, 1, howManyDiverging)
  91   ans <- unique(ans)
  92   nms <- apply(ans, 1, function(x) {
  93       paste(if(x[3]) paste(x[3], "positive, ") else "",
  94         if(x[1]) paste(x[1], 'negative, ') else "",
  95         if(x[2]) paste(x[2], 'zero') else "",
  96         if(x[4]) paste(x[4], 'diverging') else "", sep="")
  97     })
  98   rownames(ans) <- nms
  99   return(ans)
 100 }
 101
 102 print.idmc_lexp_map <- function(x, ...) {
 103   m <- x$model
 104   cat('=iDMC Lyapunov exponents map=\n')
 105   cat('Model: ', getModelName(m), '\n')
 106   cat('Starting point: ')
 107     tmp <- getModelVarNames(m)
 108     cat(paste(tmp, signif(x$var), sep=' = ', collapse=', '), '\n')
 109   cat('Parameter values: ')
 110     tmp <- getModelParNames(m)
 111     cat(paste(tmp, x$par, sep=' = ', collapse=', '), '\n')
 112   pn <- getModelParNames(m)
 113   cat('Varying x-axis par.: ', pn[x$par.x], '\n')
 114   cat('Varying x-axis par. range: [', paste(x$par.x.range, collapse=', '), ']\n')
 115   cat('Varying y-axis par.: ', pn[x$par.y], '\n')
 116   cat('Varying y-axis par. range: [', paste(x$par.y.range, collapse=', '), ']\n')
 117 }
 118
 119 as.grob.idmc_lexp_map <- function(x, colors, ...) {
 120   mat <- x$mat
 121   if(missing(colors))
 122     colors <- seq_len(max(as.vector(mat)))
 123   colors <- colors[as.vector(mat)]
 124   colors <- matrix(colors, NROW(mat))
 125   imageGrob(t(colors), xlim=x$par.x.range, ylim=x$par.y.range, respect = FALSE, name='image')
 126 }
 127
 128 plot.idmc_lexp_map <- function(x, y, colors, labels,
 129   main = getModelName(x$model),
 130   xlab, ylab, axes=TRUE, mar=NULL, legend=TRUE, ...) {
 131   m <- x$model
 132   pn <- getModelParNames(m)
 133   if(missing(xlab))
 134     xlab <- pn[x$par.x]
 135   if(missing(ylab))
 136     ylab <- pn[x$par.y]
 137   mat <- x$mat
 138   levels <- seq_along(x$labels)
 139   colors.all <- seq_along(levels)
 140   ids <- unique(as.vector(mat))
 141   if(missing(colors))
 142     colors <- colors.all[ids]
 143   else
 144     if(length(colors) < length(ids))
 145       stop(length(ids), 'colors must be specified')
 146   colors.all[ids] <- colors
 147   cG <- as.grob(x, colors=colors.all)
 148   if(legend) {
 149     labels.all <- x$labels
 150     if(missing(labels))
 151       labels <- labels.all[ids]
 152     if(length(labels) < length(ids))
 153       stop('wrong number of labels. There are ', length(ids), 'levels to be labelled')
 154     cl <- colors.all[ids]
 155     lb <- labels
 156     clg <- colorLegendGrob(cl, lb, y=unit(0, 'npc'), x=unit(0, 'npc'), name='legend')
 157     rightMargin <- convertWidth(widthDetails(clg), 'lines')
 158     mar <- c(4,4,4,rightMargin)
 159     mar[4] <- mar[4]*1.04
 160   } else {
 161     mar <- NULL
 162   }
 163   pG <- plotGrob(cG, main=main, xlab=xlab, ylab=ylab, axes=axes, mar=mar, ...)
 164   grid.draw(pG)
 165   if(legend) {
 166     downViewport('rightMarginArea')
 167     grid.draw(clg)
 168     upViewport(0)
 169   }
 170 }