Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 22:cd4d17bf3b9e 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 08:09:12 -0500
parents 80d38a26b2e8
children 5354a82ed295
comparison
equal deleted inserted replaced
21:a61f74f6d165 22:cd4d17bf3b9e
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)"),
16 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"),
17 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"),
18 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"),
19 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"),
20 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"),
21 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"),
22 make_option(c("-v", "--temperature_dataset"), action="store", dest="temperature_dataset", help="Temperature data for selected location"),
23 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)")
24 ) 22 )
25 23
26 parser <- OptionParser(usage="%prog [options] file", option_list=option_list) 24 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
27 args <- parse_args(parser, positional_arguments=TRUE) 25 args <- parse_args(parser, positional_arguments=TRUE)
28 opt <- args$options 26 opt <- args$options
29 27
30 data.input=function(loc, year, temperature.dataset) 28 get_temperature_file_path=function(loc, temperature_data)
31 { 29 {
32 expdata <- matrix(rep(0, opt$num_days * 3), nrow=opt$num_days) 30 expdata <- matrix(rep(0, opt$num_days * 3), nrow=opt$num_days)
33 namedat <- paste(loc, year, ".Rdat", sep="")
34 temp.data <- read.csv(file=temperature.dataset, header=T)
35
36 expdata[,1] <- c(1:opt$num_days) 31 expdata[,1] <- c(1:opt$num_days)
37 # Minimum 32 # Minimum
38 expdata[,2] <- temp.data[c(1:opt$num_days), 3] 33 expdata[,2] <- temperature_data[c(1:opt$num_days), 3]
39 # Maximum 34 # Maximum
40 expdata[,3] <- temp.data[c(1:opt$num_days), 2] 35 expdata[,3] <- temperature_data[c(1:opt$num_days), 2]
36 date <- temperature_data[1, 3]
37 year <- substr(date, 1, 4)
38 namedat <- paste(loc, year, ".Rdat", sep="")
41 save(expdata, file=namedat) 39 save(expdata, file=namedat)
42 namedat 40 namedat
43 } 41 }
44 42
45 daylength=function(latitude, num_days) 43 daylength=function(latitude, num_days)
209 cat("Max clutch size: ", opt$max_clutch_size, "\n") 207 cat("Max clutch size: ", opt$max_clutch_size, "\n")
210 cat("(egg->young nymph): ", opt$young_nymph_accum, "\n") 208 cat("(egg->young nymph): ", opt$young_nymph_accum, "\n")
211 cat("(young nymph->old nymph): ", opt$old_nymph_accum, "\n") 209 cat("(young nymph->old nymph): ", opt$old_nymph_accum, "\n")
212 cat("(old nymph->adult): ", opt$adult_accum, "\n") 210 cat("(old nymph->adult): ", opt$adult_accum, "\n")
213 211
214 # Read in the input temperature datafile 212 # Read in the input temperature datafile into a Data Frame object.
215 temperature_file_path <- data.input(opt$location, opt$year, opt$temperature_dataset) 213 temperature_data <= temp.data <- read.csv(file=opt$input, header=T, sep=",")
214 temperature_file_path <- get_temperature_file_path(opt$location, temperature_data)
215 latitude <- temperature_data[1, 1]
216 216
217 # Initialize matrix for results from all replications 217 # Initialize matrix for results from all replications
218 S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol = opt$replications) 218 S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol = opt$replications)
219 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) 219 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)
220 220
228 # overwintering, previttelogenic, DD=0, T=0, no-diapause 228 # overwintering, previttelogenic, DD=0, T=0, no-diapause
229 vec.mat <- rep(vec.ini, n) 229 vec.mat <- rep(vec.ini, n)
230 # complete matrix for the population 230 # complete matrix for the population
231 vec.mat <- t(matrix(vec.mat, nrow=5)) 231 vec.mat <- t(matrix(vec.mat, nrow=5))
232 # complete photoperiod profile in a year, requires daylength function 232 # complete photoperiod profile in a year, requires daylength function
233 ph.p <- daylength(opt$latitude, opt$num_days) 233 ph.p <- daylength(latitude, opt$num_days)
234 234
235 # time series of population size 235 # time series of population size
236 tot.pop <- NULL 236 tot.pop <- NULL
237 # gen.0 pop size 237 # gen.0 pop size
238 gen0.pop <- rep(0, opt$num_days) 238 gen0.pop <- rep(0, opt$num_days)
244 dd.day <- rep(0, opt$num_days) 244 dd.day <- rep(0, opt$num_days)
245 245
246 # start tick 246 # start tick
247 ptm <- proc.time() 247 ptm <- proc.time()
248 248
249 # all the days 249 # All the days included in the input temperature dataset.
250 for (day in 1:opt$num_days) { 250 for (day in 1:opt$num_days) {
251 # photoperiod in the day 251 # photoperiod in the day
252 photoperiod <- ph.p[day] 252 photoperiod <- ph.p[day]
253 temp.profile <- hourtemp(opt$latitude, day, temperature_file_path, opt$num_days) 253 temp.profile <- hourtemp(latitude, day, temperature_file_path, opt$num_days)
254 mean.temp <- temp.profile[1] 254 mean.temp <- temp.profile[1]
255 dd.temp <- temp.profile[2] 255 dd.temp <- temp.profile[2]
256 dd.day[day] <- dd.temp 256 dd.day[day] <- dd.temp
257 # trash bin for death 257 # trash bin for death
258 death.vec <- NULL 258 death.vec <- NULL
263 for (i in 1:n) { 263 for (i in 1:n) {
264 # find individual record 264 # find individual record
265 vec.ind <- vec.mat[i,] 265 vec.ind <- vec.mat[i,]
266 # first of all, still alive? 266 # first of all, still alive?
267 # adjustment for late season mortality rate 267 # adjustment for late season mortality rate
268 if (opt$latitude < 40.0) { 268 if (latitude < 40.0) {
269 post.mort <- 1 269 post.mort <- 1
270 day.kill <- 300 270 day.kill <- 300
271 } 271 }
272 else { 272 else {
273 post.mort <- 2 273 post.mort <- 2