Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 3:24fa0d35a8bf draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Thu, 09 Nov 2017 14:20:42 -0500
parents 244c373f2a34
children e7b1fc0133bb
comparison
equal deleted inserted replaced
2:07444af6824f 3:24fa0d35a8bf
4 4
5 option_list <- list( 5 option_list <- list(
6 make_option(c("-a", "--adult_mort"), action="store", dest="adult_mort", type="integer", help="Adjustment rate for adult mortality"), 6 make_option(c("-a", "--adult_mort"), action="store", dest="adult_mort", type="integer", help="Adjustment rate for adult mortality"),
7 make_option(c("-b", "--adult_accum"), action="store", dest="adult_accum", type="integer", help="Adjustment of DD accumulation (old nymph->adult)"), 7 make_option(c("-b", "--adult_accum"), action="store", dest="adult_accum", type="integer", help="Adjustment of DD accumulation (old nymph->adult)"),
8 make_option(c("-c", "--egg_mort"), action="store", dest="egg_mort", type="integer", help="Adjustment rate for egg mortality"), 8 make_option(c("-c", "--egg_mort"), action="store", dest="egg_mort", type="integer", help="Adjustment rate for egg mortality"),
9 make_option(c("-d", "--latitude"), action="store", dest="latitude", type="double", help="Latitude of selected location"),
10 make_option(c("-e", "--location"), action="store", dest="location", help="Selected location"), 9 make_option(c("-e", "--location"), action="store", dest="location", help="Selected location"),
11 make_option(c("-f", "--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"), 10 make_option(c("-f", "--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"),
12 make_option(c("-i", "--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"), 11 make_option(c("-i", "--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),
13 make_option(c("-j", "--nymph_mort"), action="store", dest="nymph_mort", type="integer", help="Adjustment rate for nymph mortality"), 12 make_option(c("-j", "--nymph_mort"), action="store", dest="nymph_mort", type="integer", help="Adjustment rate for nymph mortality"),
14 make_option(c("-k", "--old_nymph_accum"), action="store", dest="old_nymph_accum", type="integer", help="Adjustment of DD accumulation (young nymph->old nymph)"), 13 make_option(c("-k", "--old_nymph_accum"), action="store", dest="old_nymph_accum", type="integer", help="Adjustment of DD accumulation (young nymph->old nymph)"),
14 make_option(c("-n", "--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"),
15 make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"), 15 make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"),
16 make_option(c("-p", "--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"), 16 make_option(c("-p", "--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"),
17 make_option(c("-q", "--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"), 17 make_option(c("-q", "--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"),
18 make_option(c("-s", "--replications"), action="store", dest="replications", type="integer", help="Number of replications"), 18 make_option(c("-s", "--replications"), action="store", dest="replications", type="integer", help="Number of replications"),
19 make_option(c("-t", "--se_plot"), action="store", dest="se_plot", help="Plot SE"), 19 make_option(c("-t", "--se_plot"), action="store", dest="se_plot", help="Plot SE"),
20 make_option(c("-u", "--year"), action="store", dest="year", type="integer", help="Starting year"), 20 make_option(c("-v", "--input"), action="store", dest="input", help="Temperature data for selected location"),
21 make_option(c("-v", "--temperature_dataset"), action="store", dest="temperature_dataset", help="Temperature data for selected location"),
22 make_option(c("-y", "--young_nymph_accum"), action="store", dest="young_nymph_accum", type="integer", help="Adjustment of DD accumulation (egg->young nymph)") 21 make_option(c("-y", "--young_nymph_accum"), action="store", dest="young_nymph_accum", type="integer", help="Adjustment of DD accumulation (egg->young nymph)")
23 ) 22 )
24 23
25 parser <- OptionParser(usage="%prog [options] file", option_list=option_list) 24 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
26 args <- parse_args(parser, positional_arguments=TRUE) 25 args <- parse_args(parser, positional_arguments=TRUE)
27 opt <- args$options 26 opt <- args$options
28 27
29 data.input=function(loc, year, temperature.dataset) 28 convert_csv_to_rdata=function(temperature_data, data_matrix)
30 { 29 {
31 expdata <- matrix(rep(0, 365 * 3), nrow=365) 30 # Integer day of the year.
32 namedat <- paste(loc, year, ".Rdat", sep="") 31 data_matrix[,1] <- c(1:opt$num_days)
33 temp.data <- read.csv(file=temperature.dataset, header=T)
34
35 expdata[,1] <- c(1:365)
36 # Minimum 32 # Minimum
37 expdata[,2] <- temp.data[c(1:365), 3] 33 data_matrix[,2] <- temperature_data[c(1:opt$num_days), 5]
38 # Maximum 34 # Maximum
39 expdata[,3] <- temp.data[c(1:365), 2] 35 data_matrix[,3] <- temperature_data[c(1:opt$num_days), 6]
40 save(expdata, file=namedat) 36 namedat <- "tempdata.Rdat"
37 save(data_matrix, file=namedat)
41 namedat 38 namedat
42 } 39 }
43 40
44 daylength=function(latitude) 41 daylength=function(latitude, num_days)
45 { 42 {
46 # from Forsythe 1995 43 # From Forsythe 1995.
47 p=0.8333 44 p=0.8333
48 dl <- NULL 45 dl <- NULL
49 for (i in 1:365) { 46 for (i in 1:num_days) {
50 theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (i - 186))) 47 theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (i - 186)))
51 phi <- asin(0.39795 * cos(theta)) 48 phi <- asin(0.39795 * cos(theta))
52 dl[i] <- 24 - 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi))) 49 dl[i] <- 24 - 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi)))
53 } 50 }
54 dl # return a vector of daylength in 365 days 51 # Return a vector of daylength for the number of
55 } 52 # days specified in the input temperature data.
56 53 dl
57 hourtemp=function(latitude, date, temperature_file_path) 54 }
55
56 hourtemp=function(latitude, date, temperature_file_path, num_days)
58 { 57 {
59 load(temperature_file_path) 58 load(temperature_file_path)
60 threshold <- 14.17 # base development threshold for BMSB 59 # Base development threshold for Brown Marmolated Stink Bug
61 dnp <- expdata[date, 2] # daily minimum 60 # insect phenology model.
62 dxp <- expdata[date, 3] # daily maximum 61 threshold <- 14.17
62 dnp <- data_matrix[date, 2] # daily minimum
63 dxp <- data_matrix[date, 3] # daily maximum
63 dmean <- 0.5 * (dnp + dxp) 64 dmean <- 0.5 * (dnp + dxp)
64 dd <- 0 # initialize degree day accumulation 65 dd <- 0 # initialize degree day accumulation
65 66
66 if (dxp<threshold) { 67 if (dxp<threshold) {
67 dd <- 0 68 dd <- 0
68 } 69 }
69 else { 70 else {
70 dlprofile <- daylength(latitude) # extract daylength data for entire year 71 # Extract daylength data for the number of
71 T <- NULL # initialize hourly temperature 72 # days specified in the input temperature data.
72 dh <- NULL #initialize degree hour vector 73 dlprofile <- daylength(latitude, num_days)
73 # date <- 200 74 # Initialize hourly temperature.
74 y <- dlprofile[date] # calculate daylength in given date 75 T <- NULL
75 z <- 24 - y # night length 76 # Initialize degree hour vector.
76 a <- 1.86 # lag coefficient 77 dh <- NULL
77 b <- 2.20 # night coefficient 78 # Calculate daylength in given date.
78 #tempdata <- read.csv("tempdata.csv") #import raw data set 79 y <- dlprofile[date]
79 # Should be outside function otherwise its redundant 80 # Night length.
80 risetime <- 12 - y / 2 # sunrise time 81 z <- 24 - y
81 settime <- 12 + y / 2 # sunset time 82 # Lag coefficient.
83 a <- 1.86
84 # Night coefficient.
85 b <- 2.20
86 # Sunrise time.
87 risetime <- 12 - y / 2
88 # Sunset time.
89 settime <- 12 + y / 2
82 ts <- (dxp - dnp) * sin(pi * (settime - 5) / (y + 2 * a)) + dnp 90 ts <- (dxp - dnp) * sin(pi * (settime - 5) / (y + 2 * a)) + dnp
83 for (i in 1:24) { 91 for (i in 1:24) {
84 if (i > risetime && i<settime) { 92 if (i > risetime && i<settime) {
85 m <- i - 5 # number of hours after Tmin until sunset 93 # Number of hours after Tmin until sunset.
94 m <- i - 5
86 T[i]=(dxp - dnp) * sin(pi * m / (y + 2 * a)) + dnp 95 T[i]=(dxp - dnp) * sin(pi * m / (y + 2 * a)) + dnp
87 if (T[i]<8.4) { 96 if (T[i]<8.4) {
88 dh[i] <- 0 97 dh[i] <- 0
89 } 98 }
90 else { 99 else {
187 } 196 }
188 return = mort.prob 197 return = mort.prob
189 return 198 return
190 } 199 }
191 200
192 cat("Replications: ", opt$replications, "\n") 201 # Read in the input temperature datafile into a Data Frame object.
193 cat("Photoperiod: ", opt$photoperiod, "\n") 202 temperature_data <- read.csv(file=opt$input, header=T, sep=",")
194 cat("Oviposition rate: ", opt$oviposition, "\n") 203 start_date <- temperature_data[c(1:1), 3]
195 cat("Egg mortality rate: ", opt$egg_mort, "\n") 204 end_date <- temperature_data[c(opt$num_days:opt$num_days), 3]
196 cat("Nymph mortality rate: ", opt$nymph_mort, "\n") 205 raw_data_matrix <- matrix(rep(0, opt$num_days * 6), nrow=opt$num_days)
197 cat("Adult mortality rate: ", opt$adult_mort, "\n") 206 temperature_file_path <- convert_csv_to_rdata(temperature_data, raw_data_matrix)
198 cat("Min clutch size: ", opt$min_clutch_size, "\n") 207 latitude <- temperature_data[1, 1]
199 cat("Max clutch size: ", opt$max_clutch_size, "\n") 208
200 cat("(egg->young nymph): ", opt$young_nymph_accum, "\n") 209 cat("Number of days: ", opt$num_days, "\n")
201 cat("(young nymph->old nymph): ", opt$old_nymph_accum, "\n") 210
202 cat("(old nymph->adult): ", opt$adult_accum, "\n") 211 # Initialize matrix for results from all replications.
203 212 S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol = opt$replications)
204 # Read in the input temperature datafile 213 newborn.rep <- death.rep <- adult.rep <- pop.rep <- g0.rep <- g1.rep <- g2.rep <- g0a.rep <- g1a.rep <- g2a.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol=opt$replications)
205 temperature_file_path <- data.input(opt$location, opt$year, opt$temperature_dataset)
206
207 # Initialize matrix for results from all replications
208 S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, 365 * opt$replications), ncol = opt$replications)
209 newborn.rep <- death.rep <- adult.rep <- pop.rep <- g0.rep <- g1.rep <- g2.rep <- g0a.rep <- g1a.rep <- g2a.rep <- matrix(rep(0, 365 * opt$replications), ncol=opt$replications)
210 214
211 # loop through replications 215 # loop through replications
212 for (N.rep in 1:opt$replications) { 216 for (N.rep in 1:opt$replications) {
213 # during each replication 217 # During each replication start with 1000 individuals.
214 # start with 1000 individuals -- user definable as well? 218 # TODO: user definable as well?
215 n <- 1000 219 n <- 1000
216 # Generation, Stage, DD, T, Diapause 220 # Generation, Stage, DD, T, Diapause.
217 vec.ini <- c(0, 3, 0, 0, 0) 221 vec.ini <- c(0, 3, 0, 0, 0)
218 # overwintering, previttelogenic, DD=0, T=0, no-diapause 222 # Overwintering, previttelogenic, DD=0, T=0, no-diapause.
219 vec.mat <- rep(vec.ini, n) 223 vec.mat <- rep(vec.ini, n)
220 # complete matrix for the population 224 # Complete matrix for the population.
221 vec.mat <- t(matrix(vec.mat, nrow=5)) 225 vec.mat <- base::t(matrix(vec.mat, nrow=5))
222 # complete photoperiod profile in a year, requires daylength function 226 # Complete photoperiod profile in a year, requires daylength function.
223 ph.p <- daylength(opt$latitude) 227 ph.p <- daylength(latitude, opt$num_days)
224 228
225 # time series of population size 229 # Time series of population size.
226 tot.pop <- NULL 230 tot.pop <- NULL
227 # gen.0 pop size 231 gen0.pop <- rep(0, opt$num_days)
228 gen0.pop <- rep(0, 365) 232 gen1.pop <- rep(0, opt$num_days)
229 gen1.pop <- rep(0, 365) 233 gen2.pop <- rep(0, opt$num_days)
230 gen2.pop <- rep(0, 365) 234 S0 <- S1 <- S2 <- S3 <- S4 <- S5 <- rep(0, opt$num_days)
231 S0 <- S1 <- S2 <- S3 <- S4 <- S5 <- rep(0, 365) 235 g0.adult <- g1.adult <- g2.adult <- rep(0, opt$num_days)
232 g0.adult <- g1.adult <- g2.adult <- rep(0, 365) 236 N.newborn <- N.death <- N.adult <- rep(0, opt$num_days)
233 N.newborn <- N.death <- N.adult <- rep(0, 365) 237 dd.day <- rep(0, opt$num_days)
234 dd.day <- rep(0, 365) 238
235 239 # All the days included in the input temperature dataset.
236 # start tick 240 for (day in 1:opt$num_days) {
237 ptm <- proc.time() 241 # Photoperiod in the day.
238
239 # all the days
240 for (day in 1:365) {
241 # photoperiod in the day
242 photoperiod <- ph.p[day] 242 photoperiod <- ph.p[day]
243 temp.profile <- hourtemp(opt$latitude, day, temperature_file_path) 243 temp.profile <- hourtemp(latitude, day, temperature_file_path, opt$num_days)
244 mean.temp <- temp.profile[1] 244 mean.temp <- temp.profile[1]
245 dd.temp <- temp.profile[2] 245 dd.temp <- temp.profile[2]
246 dd.day[day] <- dd.temp 246 dd.day[day] <- dd.temp
247 # trash bin for death 247 # Trash bin for death.
248 death.vec <- NULL 248 death.vec <- NULL
249 # new born 249 # Newborn.
250 birth.vec <- NULL 250 birth.vec <- NULL
251 251
252 # all individuals 252 # All individuals.
253 for (i in 1:n) { 253 for (i in 1:n) {
254 # find individual record 254 # Find individual record.
255 vec.ind <- vec.mat[i,] 255 vec.ind <- vec.mat[i,]
256 # first of all, still alive? 256 # First of all, still alive?
257 # adjustment for late season mortality rate 257 # Adjustment for late season mortality rate.
258 if (opt$latitude < 40.0) { 258 if (latitude < 40.0) {
259 post.mort <- 1 259 post.mort <- 1
260 day.kill <- 300 260 day.kill <- 300
261 } 261 }
262 else { 262 else {
263 post.mort <- 2 263 post.mort <- 2
264 day.kill <- 250 264 day.kill <- 250
265 } 265 }
266 if (vec.ind[2] == 0) { 266 if (vec.ind[2] == 0) {
267 # egg 267 # Egg.
268 death.prob = opt$egg_mort * mortality.egg(mean.temp) 268 death.prob = opt$egg_mort * mortality.egg(mean.temp)
269 } 269 }
270 else if (vec.ind[2] == 1 | vec.ind[2] == 2) { 270 else if (vec.ind[2] == 1 | vec.ind[2] == 2) {
271 death.prob = opt$nymph_mort * mortality.nymph(mean.temp) 271 death.prob = opt$nymph_mort * mortality.nymph(mean.temp)
272 } 272 }
273 else if (vec.ind[2] == 3 | vec.ind[2] == 4 | vec.ind[2] == 5) { 273 else if (vec.ind[2] == 3 | vec.ind[2] == 4 | vec.ind[2] == 5) {
274 # for adult 274 # For adult.
275 if (day < day.kill) { 275 if (day < day.kill) {
276 death.prob = opt$adult_mort * mortality.adult(mean.temp) 276 death.prob = opt$adult_mort * mortality.adult(mean.temp)
277 } 277 }
278 else { 278 else {
279 # increase adult mortality after fall equinox 279 # Increase adult mortality after fall equinox.
280 death.prob = opt$adult_mort * post.mort * mortality.adult(mean.temp) 280 death.prob = opt$adult_mort * post.mort * mortality.adult(mean.temp)
281 } 281 }
282 } 282 }
283 # (or dependent on temperature and life stage?) 283 # (or dependent on temperature and life stage?)
284 u.d <- runif(1) 284 u.d <- runif(1)
285 if (u.d < death.prob) { 285 if (u.d < death.prob) {
286 death.vec <- c(death.vec, i) 286 death.vec <- c(death.vec, i)
287 } 287 }
288 else { 288 else {
289 # aggregrate index of dead bug 289 # Aggregrate index of dead bug.
290 # event 1 end of diapause 290 # Event 1 end of diapause.
291 if (vec.ind[1] == 0 && vec.ind[2] == 3) { 291 if (vec.ind[1] == 0 && vec.ind[2] == 3) {
292 # overwintering adult (previttelogenic) 292 # Overwintering adult (previttelogenic).
293 if (photoperiod > opt$photoperiod && vec.ind[3] > 68 && day < 180) { 293 if (photoperiod > opt$photoperiod && vec.ind[3] > 68 && day < 180) {
294 # add 68C to become fully reproductively matured 294 # Add 68C to become fully reproductively matured.
295 # transfer to vittelogenic 295 # Transfer to vittelogenic.
296 vec.ind <- c(0, 4, 0, 0, 0) 296 vec.ind <- c(0, 4, 0, 0, 0)
297 vec.mat[i,] <- vec.ind 297 vec.mat[i,] <- vec.ind
298 } 298 }
299 else { 299 else {
300 # add to DD 300 # Add to dd.
301 vec.ind[3] <- vec.ind[3] + dd.temp 301 vec.ind[3] <- vec.ind[3] + dd.temp
302 # add 1 day in current stage 302 # Add 1 day in current stage.
303 vec.ind[4] <- vec.ind[4] + 1 303 vec.ind[4] <- vec.ind[4] + 1
304 vec.mat[i,] <- vec.ind 304 vec.mat[i,] <- vec.ind
305 } 305 }
306 } 306 }
307 if (vec.ind[1] != 0 && vec.ind[2] == 3) { 307 if (vec.ind[1] != 0 && vec.ind[2] == 3) {
308 # NOT overwintering adult (previttelogenic) 308 # Not overwintering adult (previttelogenic).
309 current.gen <- vec.ind[1] 309 current.gen <- vec.ind[1]
310 if (vec.ind[3] > 68) { 310 if (vec.ind[3] > 68) {
311 # add 68C to become fully reproductively matured 311 # Add 68C to become fully reproductively matured.
312 # transfer to vittelogenic 312 # Transfer to vittelogenic.
313 vec.ind <- c(current.gen, 4, 0, 0, 0) 313 vec.ind <- c(current.gen, 4, 0, 0, 0)
314 vec.mat[i,] <- vec.ind 314 vec.mat[i,] <- vec.ind
315 } 315 }
316 else { 316 else {
317 # add to DD 317 # Add to dd.
318 vec.ind[3] <- vec.ind[3] + dd.temp 318 vec.ind[3] <- vec.ind[3] + dd.temp
319 # add 1 day in current stage 319 # Add 1 day in current stage.
320 vec.ind[4] <- vec.ind[4] + 1 320 vec.ind[4] <- vec.ind[4] + 1
321 vec.mat[i,] <- vec.ind 321 vec.mat[i,] <- vec.ind
322 } 322 }
323 } 323 }
324 324
325 # event 2 oviposition -- where population dynamics comes from 325 # Event 2 oviposition -- where population dynamics comes from.
326 if (vec.ind[2] == 4 && vec.ind[1] == 0 && mean.temp > 10) { 326 if (vec.ind[2] == 4 && vec.ind[1] == 0 && mean.temp > 10) {
327 # vittelogenic stage, overwintering generation 327 # Vittelogenic stage, overwintering generation.
328 if (vec.ind[4] == 0) { 328 if (vec.ind[4] == 0) {
329 # just turned in vittelogenic stage 329 # Just turned in vittelogenic stage.
330 n.birth=round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size)) 330 n.birth=round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size))
331 } 331 }
332 else { 332 else {
333 # daily probability of birth 333 # Daily probability of birth.
334 p.birth = opt$oviposition * 0.01 334 p.birth = opt$oviposition * 0.01
335 u1 <- runif(1) 335 u1 <- runif(1)
336 if (u1 < p.birth) { 336 if (u1 < p.birth) {
337 n.birth=round(runif(1, 2, 8)) 337 n.birth=round(runif(1, 2, 8))
338 } 338 }
339 } 339 }
340 # add to DD 340 # Add to dd.
341 vec.ind[3] <- vec.ind[3] + dd.temp 341 vec.ind[3] <- vec.ind[3] + dd.temp
342 # add 1 day in current stage 342 # Add 1 day in current stage.
343 vec.ind[4] <- vec.ind[4] + 1 343 vec.ind[4] <- vec.ind[4] + 1
344 vec.mat[i,] <- vec.ind 344 vec.mat[i,] <- vec.ind
345 if (n.birth > 0) { 345 if (n.birth > 0) {
346 # add new birth -- might be in different generations 346 # Add new birth -- might be in different generations.
347 # generation + 1
348 new.gen <- vec.ind[1] + 1 347 new.gen <- vec.ind[1] + 1
349 # egg profile 348 # Egg profile.
350 new.ind <- c(new.gen, 0, 0, 0, 0) 349 new.ind <- c(new.gen, 0, 0, 0, 0)
351 new.vec <- rep(new.ind, n.birth) 350 new.vec <- rep(new.ind, n.birth)
352 # update batch of egg profile 351 # Update batch of egg profile.
353 new.vec <- t(matrix(new.vec, nrow=5)) 352 new.vec <- t(matrix(new.vec, nrow=5))
354 # group with total eggs laid in that day 353 # Group with total eggs laid in that day.
355 birth.vec <- rbind(birth.vec, new.vec) 354 birth.vec <- rbind(birth.vec, new.vec)
356 } 355 }
357 } 356 }
358 357
359 # event 2 oviposition -- for gen 1. 358 # Event 2 oviposition -- for gen 1.
360 if (vec.ind[2] == 4 && vec.ind[1] == 1 && mean.temp > 12.5 && day < 222) { 359 if (vec.ind[2] == 4 && vec.ind[1] == 1 && mean.temp > 12.5 && day < 222) {
361 # vittelogenic stage, 1st generation 360 # Vittelogenic stage, 1st generation
362 if (vec.ind[4] == 0) { 361 if (vec.ind[4] == 0) {
363 # just turned in vittelogenic stage 362 # Just turned in vittelogenic stage.
364 n.birth=round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size)) 363 n.birth=round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size))
365 } 364 }
366 else { 365 else {
367 # daily probability of birth 366 # Daily probability of birth.
368 p.birth = opt$oviposition * 0.01 367 p.birth = opt$oviposition * 0.01
369 u1 <- runif(1) 368 u1 <- runif(1)
370 if (u1 < p.birth) { 369 if (u1 < p.birth) {
371 n.birth = round(runif(1, 2, 8)) 370 n.birth = round(runif(1, 2, 8))
372 } 371 }
373 } 372 }
374 # add to DD 373 # Add to dd.
375 vec.ind[3] <- vec.ind[3] + dd.temp 374 vec.ind[3] <- vec.ind[3] + dd.temp
376 # add 1 day in current stage 375 # Add 1 day in current stage.
377 vec.ind[4] <- vec.ind[4] + 1 376 vec.ind[4] <- vec.ind[4] + 1
378 vec.mat[i,] <- vec.ind 377 vec.mat[i,] <- vec.ind
379 if (n.birth > 0) { 378 if (n.birth > 0) {
380 # add new birth -- might be in different generations 379 # Add new birth -- might be in different generations.
381 # generation + 1
382 new.gen <- vec.ind[1] + 1 380 new.gen <- vec.ind[1] + 1
383 # egg profile 381 # Egg profile.
384 new.ind <- c(new.gen, 0, 0, 0, 0) 382 new.ind <- c(new.gen, 0, 0, 0, 0)
385 new.vec <- rep(new.ind, n.birth) 383 new.vec <- rep(new.ind, n.birth)
386 # update batch of egg profile 384 # Update batch of egg profile.
387 new.vec <- t(matrix(new.vec, nrow=5)) 385 new.vec <- t(matrix(new.vec, nrow=5))
388 # group with total eggs laid in that day 386 # Group with total eggs laid in that day.
389 birth.vec <- rbind(birth.vec, new.vec) 387 birth.vec <- rbind(birth.vec, new.vec)
390 } 388 }
391 } 389 }
392 390
393 # event 3 development (with diapause determination) 391 # Event 3 development (with diapause determination).
394 # event 3.1 egg development to young nymph (vec.ind[2]=0 -> egg) 392 # Event 3.1 egg development to young nymph (vec.ind[2]=0 -> egg).
395 if (vec.ind[2] == 0) { 393 if (vec.ind[2] == 0) {
396 # egg stage 394 # Egg stage.
397 # add to DD 395 # Add to dd.
398 vec.ind[3] <- vec.ind[3] + dd.temp 396 vec.ind[3] <- vec.ind[3] + dd.temp
399 if (vec.ind[3] >= (68 + opt$young_nymph_accum)) { 397 if (vec.ind[3] >= (68 + opt$young_nymph_accum)) {
400 # from egg to young nymph, DD requirement met 398 # From egg to young nymph, DD requirement met.
401 current.gen <- vec.ind[1] 399 current.gen <- vec.ind[1]
402 # transfer to young nym stage 400 # Transfer to young nymph stage.
403 vec.ind <- c(current.gen, 1, 0, 0, 0) 401 vec.ind <- c(current.gen, 1, 0, 0, 0)
404 } 402 }
405 else { 403 else {
406 # add 1 day in current stage 404 # Add 1 day in current stage.
407 vec.ind[4] <- vec.ind[4] + 1 405 vec.ind[4] <- vec.ind[4] + 1
408 } 406 }
409 vec.mat[i,] <- vec.ind 407 vec.mat[i,] <- vec.ind
410 } 408 }
411 409
412 # event 3.2 young nymph to old nymph (vec.ind[2]=1 -> young nymph: determines diapause) 410 # Event 3.2 young nymph to old nymph (vec.ind[2]=1 -> young nymph: determines diapause).
413 if (vec.ind[2] == 1) { 411 if (vec.ind[2] == 1) {
414 # young nymph stage 412 # young nymph stage.
415 # add to DD 413 # add to dd.
416 vec.ind[3] <- vec.ind[3] + dd.temp 414 vec.ind[3] <- vec.ind[3] + dd.temp
417 if (vec.ind[3] >= (250 + opt$old_nymph_accum)) { 415 if (vec.ind[3] >= (250 + opt$old_nymph_accum)) {
418 # from young to old nymph, DD requirement met 416 # From young to old nymph, dd requirement met.
419 current.gen <- vec.ind[1] 417 current.gen <- vec.ind[1]
420 # transfer to old nym stage 418 # Transfer to old nym stage.
421 vec.ind <- c(current.gen, 2, 0, 0, 0) 419 vec.ind <- c(current.gen, 2, 0, 0, 0)
422 if (photoperiod < opt$photoperiod && day > 180) { 420 if (photoperiod < opt$photoperiod && day > 180) {
423 vec.ind[5] <- 1 421 vec.ind[5] <- 1
424 } # prepare for diapausing 422 } # Prepare for diapausing.
425 } 423 }
426 else { 424 else {
427 # add 1 day in current stage 425 # Add 1 day in current stage.
428 vec.ind[4] <- vec.ind[4] + 1 426 vec.ind[4] <- vec.ind[4] + 1
429 } 427 }
430 vec.mat[i,] <- vec.ind 428 vec.mat[i,] <- vec.ind
431 } 429 }
432 430
433 # event 3.3 old nymph to adult: previttelogenic or diapausing? 431 # Event 3.3 old nymph to adult: previttelogenic or diapausing?
434 if (vec.ind[2] == 2) { 432 if (vec.ind[2] == 2) {
435 # old nymph stage 433 # Old nymph stage.
436 # add to DD 434 # add to dd.
437 vec.ind[3] <- vec.ind[3] + dd.temp 435 vec.ind[3] <- vec.ind[3] + dd.temp
438 if (vec.ind[3] >= (200 + opt$adult_accum)) { 436 if (vec.ind[3] >= (200 + opt$adult_accum)) {
439 # from old to adult, DD requirement met 437 # From old to adult, dd requirement met.
440 current.gen <- vec.ind[1] 438 current.gen <- vec.ind[1]
441 if (vec.ind[5] == 0) { 439 if (vec.ind[5] == 0) {
442 # non-diapausing adult -- previttelogenic 440 # Non-diapausing adult -- previttelogenic.
443 vec.ind <- c(current.gen, 3, 0, 0, 0) 441 vec.ind <- c(current.gen, 3, 0, 0, 0)
444 } 442 }
445 else { 443 else {
446 # diapausing 444 # Diapausing.
447 vec.ind <- c(current.gen, 5, 0, 0, 1) 445 vec.ind <- c(current.gen, 5, 0, 0, 1)
448 } 446 }
449 } 447 }
450 else { 448 else {
451 # add 1 day in current stage 449 # Add 1 day in current stage.
452 vec.ind[4] <- vec.ind[4] + 1 450 vec.ind[4] <- vec.ind[4] + 1
453 } 451 }
454 vec.mat[i,] <- vec.ind 452 vec.mat[i,] <- vec.ind
455 } 453 }
456 454
457 # event 4 growing of diapausing adult (unimportant, but still necessary)## 455 # Event 4 growing of diapausing adult (unimportant, but still necessary).
458 if (vec.ind[2] == 5) { 456 if (vec.ind[2] == 5) {
459 vec.ind[3] <- vec.ind[3] + dd.temp 457 vec.ind[3] <- vec.ind[3] + dd.temp
460 vec.ind[4] <- vec.ind[4] + 1 458 vec.ind[4] <- vec.ind[4] + 1
461 vec.mat[i,] <- vec.ind 459 vec.mat[i,] <- vec.ind
462 } 460 }
463 } # else if it is still alive 461 } # Else if it is still alive.
464 } # end of the individual bug loop 462 } # End of the individual bug loop.
465 463
466 # find how many died 464 # Find how many died.
467 n.death <- length(death.vec) 465 n.death <- length(death.vec)
468 if (n.death > 0) { 466 if (n.death > 0) {
469 vec.mat <- vec.mat[-death.vec, ] 467 vec.mat <- vec.mat[-death.vec, ]
470 } 468 }
471 # remove record of dead 469 # Remove record of dead.
472 # find how many new born 470 # Find how many new born.
473 n.newborn <- length(birth.vec[,1]) 471 n.newborn <- length(birth.vec[,1])
474 vec.mat <- rbind(vec.mat, birth.vec) 472 vec.mat <- rbind(vec.mat, birth.vec)
475 # update population size for the next day 473 # Update population size for the next day.
476 n <- n - n.death + n.newborn 474 n <- n - n.death + n.newborn
477 475
478 # aggregate results by day 476 # Aggregate results by day.
479 tot.pop <- c(tot.pop, n) 477 tot.pop <- c(tot.pop, n)
480 # egg 478 # Egg.
481 s0 <- sum(vec.mat[,2] == 0) 479 s0 <- sum(vec.mat[,2] == 0)
482 # young nymph 480 # Young nymph.
483 s1 <- sum(vec.mat[,2] == 1) 481 s1 <- sum(vec.mat[,2] == 1)
484 # old nymph 482 # Old nymph.
485 s2 <- sum(vec.mat[,2] == 2) 483 s2 <- sum(vec.mat[,2] == 2)
486 # previtellogenic 484 # Previtellogenic.
487 s3 <- sum(vec.mat[,2] == 3) 485 s3 <- sum(vec.mat[,2] == 3)
488 # vitellogenic 486 # Vitellogenic.
489 s4 <- sum(vec.mat[,2] == 4) 487 s4 <- sum(vec.mat[,2] == 4)
490 # diapausing 488 # Diapausing.
491 s5 <- sum(vec.mat[,2] == 5) 489 s5 <- sum(vec.mat[,2] == 5)
492 # overwintering adult 490 # Overwintering adult.
493 gen0 <- sum(vec.mat[,1] == 0) 491 gen0 <- sum(vec.mat[,1] == 0)
494 # first generation 492 # First generation.
495 gen1 <- sum(vec.mat[,1] == 1) 493 gen1 <- sum(vec.mat[,1] == 1)
496 # second generation 494 # Second generation.
497 gen2 <- sum(vec.mat[,1] == 2) 495 gen2 <- sum(vec.mat[,1] == 2)
498 # sum of all adults 496 # Sum of all adults.
499 n.adult <- sum(vec.mat[,2] == 3) + sum(vec.mat[,2] == 4) + sum(vec.mat[,2] == 5) 497 n.adult <- sum(vec.mat[,2] == 3) + sum(vec.mat[,2] == 4) + sum(vec.mat[,2] == 5)
500 # gen.0 pop size 498 # Gen eration 0 pop size.
501 gen0.pop[day] <- gen0 499 gen0.pop[day] <- gen0
502 gen1.pop[day] <- gen1 500 gen1.pop[day] <- gen1
503 gen2.pop[day] <- gen2 501 gen2.pop[day] <- gen2
504 S0[day] <- s0 502 S0[day] <- s0
505 S1[day] <- s1 503 S1[day] <- s1
512 g2.adult[day] <- sum((vec.mat[,1]== 2 & vec.mat[,2] == 3) | (vec.mat[,1] == 2 & vec.mat[,2] == 4) | (vec.mat[,1] == 2 & vec.mat[,2] == 5)) 510 g2.adult[day] <- sum((vec.mat[,1]== 2 & vec.mat[,2] == 3) | (vec.mat[,1] == 2 & vec.mat[,2] == 4) | (vec.mat[,1] == 2 & vec.mat[,2] == 5))
513 511
514 N.newborn[day] <- n.newborn 512 N.newborn[day] <- n.newborn
515 N.death[day] <- n.death 513 N.death[day] <- n.death
516 N.adult[day] <- n.adult 514 N.adult[day] <- n.adult
517 #print(c(N.rep, day, n, n.adult)) 515 } # end of days specified in the input temperature data
518 } # end of 365 days
519 516
520 dd.cum <- cumsum(dd.day) 517 dd.cum <- cumsum(dd.day)
521 # collect all the outputs 518 # Collect all the outputs.
522 S0.rep[,N.rep] <- S0 519 S0.rep[,N.rep] <- S0
523 S1.rep[,N.rep] <- S1 520 S1.rep[,N.rep] <- S1
524 S2.rep[,N.rep] <- S2 521 S2.rep[,N.rep] <- S2
525 S3.rep[,N.rep] <- S3 522 S3.rep[,N.rep] <- S3
526 S4.rep[,N.rep] <- S4 523 S4.rep[,N.rep] <- S4
535 g0a.rep[,N.rep] <- g0.adult 532 g0a.rep[,N.rep] <- g0.adult
536 g1a.rep[,N.rep] <- g1.adult 533 g1a.rep[,N.rep] <- g1.adult
537 g2a.rep[,N.rep] <- g2.adult 534 g2a.rep[,N.rep] <- g2.adult
538 } 535 }
539 536
540 # save(dd.day, dd.cum, S0.rep, S1.rep, S2.rep, S3.rep, S4.rep, S5.rep, newborn.rep, death.rep, adult.rep, pop.rep, g0.rep, g1.rep, g2.rep, g0a.rep, g1a.rep, g2a.rep, file=opt$output)
541 # maybe do not need to export this bit, but for now just leave it as-is
542 # do we need to export this Rdat file?
543
544 # Data analysis and visualization 537 # Data analysis and visualization
545 # default: plot 1 year of result 538 # default: plot 1 year of result
546 # but can be expanded to accommodate multiple years 539 # but can be expanded to accommodate multiple years
547 n.yr <- 1 540 n.yr <- 1
548 day.all <- c(1:365 * n.yr) 541 day.all <- c(1:opt$num_days * n.yr)
549 542
550 # mean value for adults 543 # mean value for adults
551 sa <- apply((S3.rep + S4.rep + S5.rep), 1, mean) 544 sa <- apply((S3.rep + S4.rep + S5.rep), 1, mean)
552 # mean value for nymphs 545 # mean value for nymphs
553 sn <- apply((S1.rep + S2.rep), 1,mean) 546 sn <- apply((S1.rep + S2.rep), 1,mean)
591 pdf(file=opt$output, height=20, width=20, bg="white") 584 pdf(file=opt$output, height=20, width=20, bg="white")
592 585
593 par(mar = c(5, 6, 4, 4), mfrow=c(3, 1)) 586 par(mar = c(5, 6, 4, 4), mfrow=c(3, 1))
594 587
595 # Subfigure 2: population size by life stage 588 # Subfigure 2: population size by life stage
596 plot(day.all, sa, main = "BSMB Total Population Size by Life Stage", type = "l", ylim = c(0, max(se + se.se, sn + sn.se, sa + sa.se)), axes = F, lwd = 2, xlab = "", ylab = "Number", cex = 2, cex.lab = 2, cex.axis = 2, cex.main = 2) 589 title <- paste("BSMB Total Population Size by Life Stage:", opt$location, ", Latitude:", latitude, ", Temperature Dates:", start_date, "to", end_date, sep=" ")
597 # Young and old nymphs 590 plot(day.all, sa, main=title, type="l", ylim=c(0, max(se + se.se, sn + sn.se, sa + sa.se)), axes=F, lwd=2, xlab="", ylab="Number", cex=2, cex.lab=2, cex.axis=2, cex.main=2)
598 lines(day.all, sn, lwd = 2, lty = 1, col = 2) 591 # Young and old nymphs.
592 lines(day.all, sn, lwd=2, lty=1, col=2)
599 # Eggs 593 # Eggs
600 lines(day.all, se, lwd = 2, lty = 1, col = 4) 594 lines(day.all, se, lwd=2, lty=1, col=4)
601 axis(1, at = c(1:12) * 30 - 15, cex.axis = 2, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")) 595 axis(1, at = c(1:12) * 30 - 15, cex.axis=2, labels=c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
602 axis(2, cex.axis = 2) 596 axis(2, cex.axis = 2)
603 leg.text <- c("Egg", "Nymph", "Adult") 597 leg.text <- c("Egg", "Nymph", "Adult")
604 legend("topleft", leg.text, lty = c(1, 1, 1), col = c(4, 2, 1), cex = 2) 598 legend("topleft", leg.text, lty=c(1, 1, 1), col=c(4, 2, 1), cex=2)
605 if (opt$se_plot == 1) { 599 if (opt$se_plot == 1) {
606 # add SE lines to plot 600 # add SE lines to plot
607 # SE for adults 601 # SE for adults
608 lines (day.all, sa + sa.se, lty = 2) 602 lines (day.all, sa + sa.se, lty=2)
609 lines (day.all, sa - sa.se, lty = 2) 603 lines (day.all, sa - sa.se, lty=2)
610 # SE for nymphs 604 # SE for nymphs
611 lines (day.all, sn + sn.se, col = 2, lty = 2) 605 lines (day.all, sn + sn.se, col=2, lty=2)
612 lines (day.all, sn - sn.se, col = 2, lty = 2) 606 lines (day.all, sn - sn.se, col=2, lty=2)
613 # SE for eggs 607 # SE for eggs
614 lines (day.all, se + se.se, col = 4, lty = 2) 608 lines (day.all, se + se.se, col=4, lty=2)
615 lines (day.all, se - se.se, col = 4, lty = 2) 609 lines (day.all, se - se.se, col=4, lty=2)
616 } 610 }
617 611
618 # Subfigure 3: population size by generation 612 # Subfigure 3: population size by generation
619 plot(day.all, g0, main = "BSMB Total Population Size by Generation", type = "l", ylim = c(0, max(g2)), axes = F, lwd = 2, xlab = "", ylab = "Number", cex = 2, cex.lab = 2, cex.axis = 2, cex.main = 2) 613 title <- paste("BSMB Total Population Size by Generation:", opt$location, ", Latitude:", latitude, ", Temperature Dates:", start_date, "to", end_date, sep=" ")
614 plot(day.all, g0, main=title, type="l", ylim=c(0, max(g2)), axes=F, lwd=2, xlab="", ylab="Number", cex=2, cex.lab=2, cex.axis=2, cex.main=2)
620 lines(day.all, g1, lwd = 2, lty = 1, col = 2) 615 lines(day.all, g1, lwd = 2, lty = 1, col = 2)
621 lines(day.all, g2, lwd = 2, lty = 1, col = 4) 616 lines(day.all, g2, lwd = 2, lty = 1, col = 4)
622 axis(1, at = c(1:12) * 30 - 15, cex.axis = 2, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")) 617 axis(1, at = c(1:12) * 30 - 15, cex.axis = 2, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
623 axis(2, cex.axis = 2) 618 axis(2, cex.axis = 2)
624 leg.text <- c("P", "F1", "F2") 619 leg.text <- c("P", "F1", "F2")
635 lines (day.all, g2 + g2.se, col = 4, lty = 2) 630 lines (day.all, g2 + g2.se, col = 4, lty = 2)
636 lines (day.all, g2 - g2.se, col = 4, lty = 2) 631 lines (day.all, g2 - g2.se, col = 4, lty = 2)
637 } 632 }
638 633
639 # Subfigure 4: adult population size by generation 634 # Subfigure 4: adult population size by generation
640 plot(day.all, g0a, ylim = c(0, max(g2a) + 100), main = "BSMB Adult Population Size by Generation", type = "l", axes = F, lwd = 2, xlab = "Year", ylab = "Number", cex = 2, cex.lab = 2, cex.axis = 2, cex.main = 2) 635 title <- paste("BSMB Adult Population Size by Generation:", opt$location, ", Latitude:", latitude, ", Temperature Dates:", start_date, "to", end_date, sep=" ")
636 plot(day.all, g0a, ylim=c(0, max(g2a) + 100), main=title, type="l", axes=F, lwd=2, xlab="Year", ylab="Number", cex=2, cex.lab=2, cex.axis=2, cex.main=2)
641 lines(day.all, g1a, lwd = 2, lty = 1, col = 2) 637 lines(day.all, g1a, lwd = 2, lty = 1, col = 2)
642 lines(day.all, g2a, lwd = 2, lty = 1, col = 4) 638 lines(day.all, g2a, lwd = 2, lty = 1, col = 4)
643 axis(1, at = c(1:12) * 30 - 15, cex.axis = 2, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")) 639 axis(1, at = c(1:12) * 30 - 15, cex.axis = 2, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
644 axis(2, cex.axis = 2) 640 axis(2, cex.axis = 2)
645 leg.text <- c("P", "F1", "F2") 641 leg.text <- c("P", "F1", "F2")