-
Notifications
You must be signed in to change notification settings - Fork 14
/
example_POM.R
167 lines (145 loc) · 6.84 KB
/
example_POM.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
if (interactive()) {
set.seed(89462)
library(parallel)
library(raster)
mySim <- simInit(
times = list(start = 0.0, end = 2.0, timeunit = "year"),
params = list(
.globals = list(stackName = "landscape", burnStats = "nPixelsBurned"),
fireSpread = list(nfires = 5),
randomLandscapes = list(nx = 300, ny = 300)
),
modules = list("randomLandscapes", "fireSpread", "caribouMovement"),
paths = list(modulePath = system.file("sampleModules", package = "SpaDES.core"))
)
# Since this is a made up example, we don't have real data
# to run POM against. Instead, we will run the model once,
# take the values at the end of the simulation as if they
# are real data, then rerun the POM function next,
# comparing these "data" with the simulated values
# using Mean Absolute Deviation
outData <- spades(reproducible::Copy(mySim), .plotInitialTime = NA)
# Extract the "true" data, in this case, the "proportion of cells burned"
# Function defined that will use landscape$Fires map from simList,
# i.e., sim$landscape$Fires
# the return value being compared via MAD with propCellBurnedData
propCellBurnedFn <- function(landscape) {
sum(getValues(landscape$Fires) > 0) / ncell(landscape$Fires)
}
# visualize the burned maps of true "data"
propCellBurnedData <- propCellBurnedFn(outData$landscape)
clearPlot()
if (interactive()) {
fires <- outData$landscape$Fires # Plot doesn't do well with many nested layers
Plot(fires)
}
# Example 1 - 1 parameter
# In words, this says, "find the best value of spreadprob such that
# the proportion of the area burned in the simulation
# is as close as possible to the proportion area burned in
# the "data", using \code{DEoptim()}.
# Can use cluster if computer is multi-threaded.
# This example can use parallelType = 1 in DEoptim. For this, you must manually
# pass all packages and variables as character strings.
# cl <- makeCluster(detectCores() - 1) # not implemented yet in DEoptim
out1 <- POM(mySim, "spreadprob",
list(propCellBurnedData = propCellBurnedFn), # data = pattern pair
#optimControl = list(parallelType = 1),
logObjFnVals = TRUE)
## Once cl arg is available from DEoptim, this will work:
# out1 <- POM(mySim, "spreadprob", cl = cl,
# list(propCellBurnedData = propCellBurnedFn)) # data = pattern pair
# Example 2 - 2 parameters
# Function defined that will use caribou from sim$caribou, with
# the return value being compared via MAD with nPattern
# module, parameter N, is from 10 to 1000)
caribouFn <- function(caribou) length(caribou)
# Extract "data" from simList object (normally, this would be actual data)
nPattern <- caribouFn(outData$caribou)
aTime <- Sys.time()
parsToVary <- c("spreadprob", "N")
out2 <- POM(mySim, parsToVary,
list(propCellBurnedData = propCellBurnedFn,
nPattern = caribouFn), logObjFnVals = TRUE)
#optimControl = list(parallelType = 1))
#cl = cl) # not yet implemented, waiting for DEoptim
bTime <- Sys.time()
# check that population overlaps known values (0.225 and 100)
apply(out2$member$pop, 2, quantile, c(0.025, 0.975))
hists <- apply(out2$member$pop, 2, hist, plot = FALSE)
clearPlot()
for (i in seq_along(hists)) Plot(hists[[i]], addTo = parsToVary[i],
title = parsToVary[i], axes = TRUE)
print(paste("DEoptim", format(bTime - aTime)))
#stopCluster(cl) # not yet implemented, waiting for DEoptim
# Example 3 - using objFn instead of objects
# list all the parameters in the simList, from these, we select to vary
params(mySim)
# Objective Function Example:
# objective function must have several elements
# - first argument must be parameter vector, passed to and used by DEoptim
# - likely needs to take sim object, likely needs a copy
# because of pass-by-reference semantics of sim objects
# - pass data that will be used internally for objective function
objFnEx <- function(pars, # param values
sim, # simList object
nPattern, propCellBurnedData, caribouFn, propCellBurnedFn) {
### data
# make a copy of simList because it will possibly be altered by spades call
sim1 <- reproducible::Copy(sim)
# take the parameters and assign them to simList
params(sim1)$fireSpread$spreadprob <- pars[1]
params(sim1)$caribouMovement$N <- pars[2]
# run spades, without plotting
out <- spades(sim1, .plotInitialTime = NA)
# calculate outputs
propCellBurnedOut <- propCellBurnedFn(out$landscape)
nPattern_Out <- caribouFn(out$caribou)
minimizeFn <- abs(nPattern_Out - nPattern) +
abs(propCellBurnedOut - propCellBurnedData)
# have more info reported to console, if desired
# cat(minimizeFn)
# cat(" ")
# cat(pars)
# cat("\n")
return(minimizeFn)
}
# Run DEoptim with custom objFn, identifying 2 parameters to allow
# to vary, and pass all necessary objects required for the
# objFn
# choose 2 of them to vary. Need to identify them in params & inside objFn
# Change optimization parameters to alter how convergence is achieved
out5 <- POM(mySim, params = c("spreadprob", "N"),
objFn = objFnEx,
nPattern = nPattern,
propCellBurnedData = propCellBurnedData,
caribouFn = caribouFn,
propCellBurnedFn = propCellBurnedFn,
#cl = cl, # uncomment for cluster # not yet implemented, waiting for DEoptim
# see ?DEoptim.control for explanation of these options
optimControl = list(
NP = 100, # run 100 populations, allowing quantiles to be calculated
initialpop = matrix(c(runif(100, 0.2, 0.24), runif(100, 80, 120)), ncol = 2),
parallelType = 1
)
)
# Can also use an optimizer directly -- miss automatic parameter bounds,
# and automatic objective function using option 2
library(DEoptim)
out7 <- DEoptim(fn = objFnEx,
sim = mySim,
nPattern = nPattern,
propCellBurnedData = propCellBurnedData,
caribouFn = caribouFn,
propCellBurnedFn = propCellBurnedFn,
# cl = cl, # uncomment for cluster
# see ?DEoptim.control for explanation of these options
control = DEoptim.control(
steptol = 3,
parallelType = 1, # parallelType = 3,
packages = list("raster", "SpaDES.core", "RColorBrewer"),
parVar = list("objFnEx"),
initialpop = matrix(c(runif(40, 0.2, 0.24),
runif(40, 80, 120)), ncol = 2)),
lower = c(0.2, 80), upper = c(0.24, 120))
}