| blob | 4cd0022a8e0415102f021b06ac5e559b60d56f30 |
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
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 }
67 grid.draw.idmc_basin <- function(x, y, color.attractors, color.basins,
68 color.infinity, axes=TRUE, ...) {
69 mat <- getBasinData(x)
70 vals <- unique(as.vector(mat))
71 vals <- vals[vals>1]
72 attrCodes <- vals[(vals %% 2)==0]
73 mat1 <- matrix(1, NROW(mat), NCOL(mat))
74 for(i in seq_along(attrCodes)) {
75 mat1[mat==attrCodes[i]] <- (i-1)*2 + 2
76 mat1[mat==(attrCodes[i]+1)] <- (i-1)*2 + 3
77 }
78 nl <- length(vals)
79 na <- length(attrCodes) ##number of attractors
80 nc <- NCOL(mat)
81 mat1 <- t(mat1[,nc:1])
82 mdl <- getBasinModel(x)
83 default.palette <- rainbow(na*2)
84 if(missing(color.attractors))
85 color.attractors <- default.palette[seq(1, by=2, length=na)]
86 if(missing(color.basins))
87 color.basins <- default.palette[seq(2, by=2, length=na)]
88 if(missing(color.infinity))
89 color.infinity <- 'black'
90 if(length(color.attractors)<na)
91 color.attractors <- c(color.attractors, default.palette[seq(1, by=2, length=na)])
92 if(length(color.basins)<na)
93 color.attractors <- c(color.attractors, default.palette[seq(2, by=2, length=na)])
94 col <- numeric(2*na+1)
95 col[1] <- color.infinity
96 col[1+seq(1, by=2, length=na)] <- color.attractors[seq_len(na)]
97 col[1+seq(2, by=2, length=na)] <- color.basins[seq_len(na)]
98 imG <- imageDiscretePlotGrob(mat1, col, xlim=x$xlim, ylim=x$ylim, axes=axes,
99 name='image', gp=NULL, vp=NULL)
100 grid.draw(imG)
101 }
103 plot.idmc_basin <- function(x, y, color.attractors, color.basins,
104 color.infinity, main, xlab, ylab, ...) {
105 mat <- getBasinData(x)
106 vals <- unique(as.vector(mat))
107 vals <- vals[vals>1]
108 attrCodes <- vals[(vals %% 2)==0]
109 mat1 <- matrix(1, NROW(mat), NCOL(mat))
110 for(i in seq_along(attrCodes)) {
111 mat1[mat==attrCodes[i]] <- (i-1)*2 + 2
112 mat1[mat==(attrCodes[i]+1)] <- (i-1)*2 + 3
113 }
114 nl <- length(vals)
115 na <- length(attrCodes) ##number of attractors
116 nc <- NCOL(mat)
117 mdl <- getBasinModel(x)
118 default.palette <- rainbow(na*2)
119 if(missing(color.attractors))
120 color.attractors <- default.palette[seq(1, by=2, length=na)]
121 if(missing(color.basins))
122 color.basins <- default.palette[seq(2, by=2, length=na)]
123 if(missing(color.infinity))
124 color.infinity <- 'black'
125 if(length(color.attractors)<na)
126 color.attractors <- c(color.attractors, default.palette[seq(1, by=2, length=na)])
127 if(length(color.basins)<na)
128 color.attractors <- c(color.attractors, default.palette[seq(2, by=2, length=na)])
129 if(missing(main))
130 main <- getModelName(mdl)
131 if(missing(xlab))
132 xlab <- getModelVarNames(mdl)[1]
133 if(missing(ylab))
134 ylab <- getModelVarNames(mdl)[2]
135 col <- numeric(2*na+1)
136 col[1] <- color.infinity
137 col[1+seq(1, by=2, length=na)] <- color.attractors[seq_len(na)]
138 col[1+seq(2, by=2, length=na)] <- color.basins[seq_len(na)]
139 image(x=seq(x$xlim[1], x$xlim[2], length=x$xres),
140 y=seq(x$ylim[1], x$ylim[2], length=x$yres),
141 z=mat1[,nc:1],
142 breaks=c( 0:(na*2+1) + 0.5 ) , col=col,
143 main=main, xlab=xlab, ylab=ylab, ... )
144 }