insect_phenology_model: insect_phenology

comparison insect_phenology_model.R @ 13:80d38a26b2e8 draft

Uploaded

author	greg
date	Thu, 28 Sep 2017 10:45:53 -0400
parents	4114785abbe7
children	cd4d17bf3b9e

comparison

equal deleted inserted replaced

-:e82ea5f1a179
+:80d38a26b2e8
 make_option(c("-e", "--location"), action="store", dest="location", help="Selected location"),
 make_option(c("-f", "--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"),
 make_option(c("-i", "--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),
 make_option(c("-j", "--nymph_mort"), action="store", dest="nymph_mort", type="integer", help="Adjustment rate for nymph mortality"),
 make_option(c("-k", "--old_nymph_accum"), action="store", dest="old_nymph_accum", type="integer", help="Adjustment of DD accumulation (young nymph->old nymph)"),
+make_option(c("-n", "--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"),
 make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"),
 make_option(c("-p", "--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"),
 make_option(c("-q", "--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"),
 make_option(c("-s", "--replications"), action="store", dest="replications", type="integer", help="Number of replications"),
 make_option(c("-t", "--se_plot"), action="store", dest="se_plot", help="Plot SE"),
 args <- parse_args(parser, positional_arguments=TRUE)
 opt <- args$options
 data.input=function(loc, year, temperature.dataset)
 {
-expdata <- matrix(rep(0, 365 * 3), nrow=365)
+expdata <- matrix(rep(0, opt$num_days * 3), nrow=opt$num_days)
 namedat <- paste(loc,  year, ".Rdat", sep="")
 temp.data <- read.csv(file=temperature.dataset, header=T)
-expdata[,1] <- c(1:365)
+expdata[,1] <- c(1:opt$num_days)
 # Minimum
-expdata[,2] <- temp.data[c(1:365), 3]
+expdata[,2] <- temp.data[c(1:opt$num_days), 3]
 # Maximum
-expdata[,3] <- temp.data[c(1:365), 2]
+expdata[,3] <- temp.data[c(1:opt$num_days), 2]
 save(expdata, file=namedat)
 namedat
 }
-daylength=function(latitude)
+daylength=function(latitude, num_days)
 {
 # from Forsythe 1995
 p=0.8333
 dl <- NULL
-for (i in 1:365) {
+for (i in 1:num_days) {
 theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (i - 186)))
 phi <- asin(0.39795 * cos(theta))
 dl[i] <- 24 - 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi)))
 }
-dl   # return a vector of daylength in 365 days
+# return a vector of daylength for the number of
-}
+# days specifie din the input temperature data
+dl
-hourtemp=function(latitude, date, temperature_file_path)
+}
+hourtemp=function(latitude, date, temperature_file_path, num_days)
 {
 load(temperature_file_path)
-threshold <- 14.17  # base development threshold for BMSB
+# base development threshold for BMSB
+threshold <- 14.17
 dnp <- expdata[date, 2]  # daily minimum
 dxp <- expdata[date, 3]  # daily maximum
 dmean <- 0.5 * (dnp + dxp)
 dd <- 0  # initialize degree day accumulation
 if (dxp<threshold) {
 dd <- 0
 }
 else {
-dlprofile <- daylength(latitude)  # extract daylength data for entire year
+# extract daylength data for the number of
+# days specified in the input temperature data
+dlprofile <- daylength(latitude, num_days)
 T <- NULL  # initialize hourly temperature
 dh <- NULL #initialize degree hour vector
-# date <- 200
+# calculate daylength in given date
-y <- dlprofile[date]  # calculate daylength in given date
+y <- dlprofile[date]
-z <- 24 - y     # night length
+# night length
-a <- 1.86     # lag coefficient
+z <- 24 - y
-b <- 2.20     # night coefficient
+# lag coefficient
-#tempdata <- read.csv("tempdata.csv") #import raw data set
+a <- 1.86
-# Should be outside function otherwise its redundant
+# night coefficient
-risetime <- 12 - y / 2      # sunrise time
+b <- 2.20
-settime <- 12 + y / 2       # sunset time
+# sunrise time
+risetime <- 12 - y / 2
+# sunset time
+settime <- 12 + y / 2
 ts <- (dxp - dnp) * sin(pi * (settime - 5) / (y + 2 * a)) + dnp
 for (i in 1:24) {
 if (i > risetime && i<settime) {
-m <- i - 5  # number of hours after Tmin until sunset
+# number of hours after Tmin until sunset
+m <- i - 5
 T[i]=(dxp - dnp) * sin(pi * m / (y + 2 * a)) + dnp
 if (T[i]<8.4) {
 dh[i] <- 0
 }
 else {
 # Read in the input temperature datafile
 temperature_file_path <- data.input(opt$location, opt$year, opt$temperature_dataset)
 # Initialize matrix for results from all replications
-S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, 365 * opt$replications), ncol = opt$replications)
+S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol = opt$replications)
-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)
+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)
 # loop through replications
 for (N.rep in 1:opt$replications) {
 # during each replication
 # start with 1000 individuals -- user definable as well?
 # overwintering, previttelogenic, DD=0, T=0, no-diapause
 vec.mat <- rep(vec.ini, n)
 # complete matrix for the population
 vec.mat <- t(matrix(vec.mat, nrow=5))
 # complete photoperiod profile in a year, requires daylength function
-ph.p <- daylength(opt$latitude)
+ph.p <- daylength(opt$latitude, opt$num_days)
 # time series of population size
 tot.pop <- NULL
 # gen.0 pop size
-gen0.pop <- rep(0, 365)
+gen0.pop <- rep(0, opt$num_days)
-gen1.pop <- rep(0, 365)
+gen1.pop <- rep(0, opt$num_days)
-gen2.pop <- rep(0, 365)
+gen2.pop <- rep(0, opt$num_days)
-S0 <- S1 <- S2 <- S3 <- S4 <- S5 <- rep(0, 365)
+S0 <- S1 <- S2 <- S3 <- S4 <- S5 <- rep(0, opt$num_days)
-g0.adult <- g1.adult <- g2.adult <- rep(0, 365)
+g0.adult <- g1.adult <- g2.adult <- rep(0, opt$num_days)
-N.newborn <- N.death <- N.adult <- rep(0, 365)
+N.newborn <- N.death <- N.adult <- rep(0, opt$num_days)
-dd.day <- rep(0, 365)
+dd.day <- rep(0, opt$num_days)
 # start tick
 ptm <- proc.time()
 # all the days
-for (day in 1:365) {
+for (day in 1:opt$num_days) {
 # photoperiod in the day
 photoperiod <- ph.p[day]
-temp.profile <- hourtemp(opt$latitude, day, temperature_file_path)
+temp.profile <- hourtemp(opt$latitude, day, temperature_file_path, opt$num_days)
 mean.temp <- temp.profile[1]
 dd.temp <- temp.profile[2]
 dd.day[day] <- dd.temp
 # trash bin for death
 death.vec <- NULL
 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))
 N.newborn[day] <- n.newborn
 N.death[day] <- n.death
 N.adult[day] <- n.adult
-#print(c(N.rep, day, n, n.adult))
+}   # end of days specified in the input temperature data
-}   # end of 365 days
 dd.cum <- cumsum(dd.day)
 # collect all the outputs
 S0.rep[,N.rep] <- S0
 S1.rep[,N.rep] <- S1
 g0a.rep[,N.rep] <- g0.adult
 g1a.rep[,N.rep] <- g1.adult
 g2a.rep[,N.rep] <- g2.adult
 }
-# 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)
-# maybe do not need to export this bit, but for now just leave it as-is
-# do we need to export this Rdat file?
 # Data analysis and visualization
 # default: plot 1 year of result
 # but can be expanded to accommodate multiple years
 n.yr <- 1
-day.all <- c(1:365 * n.yr)
+day.all <- c(1:opt$num_days * n.yr)
 # mean value for adults
 sa <- apply((S3.rep + S4.rep + S5.rep), 1, mean)
 # mean value for nymphs
 sn <- apply((S1.rep + S2.rep), 1,mean)

Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 13:80d38a26b2e8 draft