RiDMC/R/basin.R

   1 Basin <- function(model, parms, xlim, xres=100, ylim,
   2   yres=100, attractorLimit, attractorIterations,
   3   method=c("fast", "slow"), ntries, seed) {
   4   checkModelParVar(model, parms)
   5   checkPositiveScalar(xres)
   6   checkPositiveScalar(yres)
   7   checkPositiveScalar(attractorLimit)
   8   checkPositiveScalar(attractorIterations)
   9   checkModelDiscrete(model)
  10   if(length(getModelVarNames(model))!=2)
  11     stop("'model' should have exactly 2 variables")
  12   method <- match.arg(method)
  13   ans <- list()
  14   ans$model <- model
  15   ans$xlim <- xlim
  16   ans$xres <- xres
  17   ans$ylim <- ylim
  18   ans$yres <- yres
  19   ans$attractorLimit <- attractorLimit
  20   ans$attractorIterations <- attractorIterations
  21   ans$method <- method
  22   if(method=="fast") {
  23     basin <- .Call("ridmc_basin_alloc", model$model, as.double(parms),
  24       as.double(xlim[1]), as.double(xlim[2]), as.integer(xres),
  25       as.double(ylim[1]), as.double(ylim[2]), as.integer(yres),
  26       as.integer(attractorLimit), as.integer(attractorIterations) ,
  27       PACKAGE='RiDMC')
  28     while(.Call("ridmc_basin_finished", basin, PACKAGE='RiDMC')==0)
  29       .Call("ridmc_basin_step", basin, PACKAGE='RiDMC')
  30     ans$data <- .Call("ridmc_basin_getData", basin, PACKAGE='RiDMC')
  31   } else {
  32     if(missing(ntries))
  33       stop("'ntries' is mandatory with 'slow' method")
  34     checkPositiveScalar(ntries)
  35     ans$ntries <- as.integer(ntries)
  36     basin <- .Call("ridmc_basin_slow_alloc", model$model, as.double(parms),
  37       as.double(xlim[1]), as.double(xlim[2]), as.integer(xres),
  38       as.double(ylim[1]), as.double(ylim[2]), as.integer(yres),
  39       as.integer(attractorLimit), as.integer(attractorIterations),
  40       as.integer(ntries), PACKAGE='RiDMC')
  41     if(!missing(seed))
  42       .Call("ridmc_basin_slow_setGslRngSeed", basin,
  43         as.integer(seed), PACKAGE='RiDMC')
  44     while(.Call("ridmc_basin_slow_finished", basin, PACKAGE='RiDMC')==0)
  45       .Call("ridmc_basin_slow_step", basin, PACKAGE='RiDMC')
  46     ans$data <- .Call("ridmc_basin_slow_getData", basin, PACKAGE='RiDMC')
  47   }
  48   class(ans) <- "idmc_basin"
  49   return(ans)
  50 }
  51 getBasinModel <- function(obj, ...)
  52   obj$model
  53 getBasinData <- function(obj, ...)
  54   obj$data
  55
  56 print.idmc_basin <- function(x, ...){
  57   mdl <- getBasinModel(x)
  58   cat('= iDMC basins of attraction =\n')
  59   cat('Model: ', getModelName(mdl), '\n')
  60   cat('x-range: ', paste(x$xlim, collapse=', '), '\n')
  61   cat('y-range: ', paste(x$ylim, collapse=', '), '\n')
  62   cat('resolution: ', x$xres, 'by', x$yres, '\n')
  63   cat('transient: ', x$attractorLimit, '\n')
  64   cat('attractor iterations: ', x$attractorIterations, '\n')
  65 }
  66
  67 makeBasinsPalette <- function(values, color.attractors, color.basins, color.infinity, default.palette=rainbow) {
  68   values <- unique(values)
  69   attrCodes <- values[(values %% 2)==0] ##attractor codes
  70   nl <- length(values) ##number of levels
  71   na <- length(attrCodes) ##number of attractors
  72   default.palette <- default.palette(na*2)
  73   if(missing(color.attractors))
  74     color.attractors <- default.palette[seq(1, by=2, length=na)]
  75   if(missing(color.basins))
  76     color.basins <- default.palette[seq(2, by=2, length=na)]
  77   if(missing(color.infinity))
  78     color.infinity <- 'black'
  79   if(length(color.attractors)<na)
  80     color.attractors <- c(color.attractors, default.palette[seq(1, by=2, length=na)])
  81   if(length(color.basins)<na)
  82     color.attractors <- c(color.attractors, default.palette[seq(2, by=2, length=na)])
  83   col <- numeric(2*na+1)
  84   col[1] <- color.infinity
  85   col[1+seq(1, by=2, length=na)] <- color.attractors[seq_len(na)]
  86   col[1+seq(2, by=2, length=na)] <- color.basins[seq_len(na)]
  87   col
  88 }
  89
  90 as.grob.idmc_basin <- function(x, color.attractors, color.basins,
  91   color.infinity, ...) {
  92   mat <- getBasinData(x)
  93   vals <- unique(as.vector(mat))
  94   vals <- vals[vals>1]
  95   attrCodes <- vals[(vals %% 2)==0]
  96   na <- length(attrCodes)
  97   mat1 <- matrix(1, NROW(mat), NCOL(mat))
  98   for(i in seq_along(attrCodes)) {
  99     mat1[mat==attrCodes[i]] <- (i-1)*2 + 2
 100     mat1[mat==(attrCodes[i]+1)] <- (i-1)*2 + 3
 101   }
 102   nc <- NCOL(mat)
 103   mat1 <- t(mat1[,nc:1])
 104   col <- makeBasinsPalette(values=vals, color.attractors, color.basins, color.infinity)
 105   ans <- imageGrob(matrix(col[as.vector(mat1)], NROW(mat1), NCOL(mat1)),
 106     xlim=x$xlim, ylim=x$ylim, respect = FALSE, name='image')
 107 }
 108
 109 plot.idmc_basin <- function(x, y, color.attractors, color.basins,
 110   color.infinity, labels.attr, labels.bas, label.infty='infinity',
 111   main = getModelName(getBasinModel(x)),
 112   xlab = getModelVarNames(getBasinModel(x))[1],
 113   ylab = getModelVarNames(getBasinModel(x))[2],
 114   axes=TRUE, legend=FALSE, attractorPoints=FALSE,
 115   pch=16, cex=0.2, ...) {
 116   imG <- as.grob(x, color.attractors=color.attractors, color.basins=color.basins, color.infinity=color.infinity)
 117   data <- getBasinData(x)
 118   vals <- unique(as.vector(data))
 119   vals <- vals[vals>1]
 120   attrCodes <- vals[(vals %% 2) == 0]
 121   na <- length(attrCodes)
 122   col <- makeBasinsPalette(values=vals, color.attractors, color.basins, color.infinity)
 123   if(legend) {
 124     if(missing(labels.attr))
 125       labels.attr <- paste('attractor', seq_len(na))
 126     if(missing(labels.bas))
 127       labels.bas <- paste('basin', seq_len(na))
 128     if(length(labels.attr)<na || length(labels.bas) < na)
 129       stop('there are ', na, 'attractors/basins pairs to plot: not enough labels provided')
 130     labels <- character(na*2+1)
 131     labels[1] <- label.infty
 132     labels[1+seq(1, by=2, length=na)] <- labels.attr
 133     labels[1+seq(2, by=2, length=na)] <- labels.bas
 134     yl <- unit(0, 'npc')
 135     xl <- unit(0, 'npc')
 136     clg <- colorLegendGrob(col, labels, y=yl, x=xl, name='legend')
 137     rightMargin <- convertWidth(widthDetails(clg), 'lines')
 138     mar <- c(4,4,4,rightMargin)
 139     mar[4] <- mar[4]*1.04
 140   } else
 141     mar <- NULL
 142   pG <- plotGrob(imG, axes=axes, main=main, xlab=xlab, ylab=ylab, mar=mar)
 143   grid.draw(pG)
 144   if(legend) {
 145     downViewport('rightMarginArea')
 146     grid.draw(clg)
 147     upViewport(0)
 148   }
 149   if(attractorPoints) {
 150     downViewport('plotArea')
 151     for(i in seq_along(attrCodes)) { ##for each attractor
 152       acd <- attrCodes[i]
 153       ids <- which(data==acd, arr.ind=TRUE)
 154       grid.points(ids[,1]/NCOL(data), 1-ids[,2]/NROW(data),
 155         pch=pch, size=unit(cex, 'char'),
 156         default.unit='npc', gp=gpar(col=col[i*2]))
 157     }
 158     upViewport(0)
 159   }
 160   invisible(NULL)
 161 }