Mercurial > repos > greg > insect_phenology_model
changeset 0:244c373f2a34 draft
Uploaded
| author | greg |
|---|---|
| date | Tue, 08 Aug 2017 13:14:39 -0400 |
| parents | |
| children | 3debfb664792 |
| files | .shed.yml insect_phenology_model.R insect_phenology_model.xml test-data/asheville2014.csv test-data/output.pdf tool-data/locations.txt.sample tool-data/years.txt.sample |
| diffstat | 7 files changed, 1231 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/.shed.yml Tue Aug 08 13:14:39 2017 -0400 @@ -0,0 +1,13 @@ +name: insect_phenology_model +owner: greg +description: | + Contains a tool that provides an agent-based stochastic model expressing stage-specific phenology and population dynamics + for an insect species across geographic regions. +homepage_url: https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/entomology/insect_phenology_model +long_description: | + Contains a tool that provides an agent-based stochastic model expressing stage-specific phenology and population dynamics + for an insect species across geographic regions. +remote_repository_url: https://github.com/gregvonkuster/galaxy_tools/tree/master/tools/entomology/insect_phenology_model +type: unrestricted +categories: +- Entomology
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/insect_phenology_model.R Tue Aug 08 13:14:39 2017 -0400 @@ -0,0 +1,661 @@ +#!/usr/bin/env Rscript + +suppressPackageStartupMessages(library("optparse")) + +option_list <- list( + make_option(c("-a", "--adult_mort"), action="store", dest="adult_mort", type="integer", help="Adjustment rate for adult mortality"), + make_option(c("-b", "--adult_accum"), action="store", dest="adult_accum", type="integer", help="Adjustment of DD accumulation (old nymph->adult)"), + make_option(c("-c", "--egg_mort"), action="store", dest="egg_mort", type="integer", help="Adjustment rate for egg mortality"), + make_option(c("-d", "--latitude"), action="store", dest="latitude", type="double", help="Latitude of selected location"), + 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("-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"), + make_option(c("-u", "--year"), action="store", dest="year", type="integer", help="Starting year"), + make_option(c("-v", "--temperature_dataset"), action="store", dest="temperature_dataset", help="Temperature data for selected location"), + make_option(c("-y", "--young_nymph_accum"), action="store", dest="young_nymph_accum", type="integer", help="Adjustment of DD accumulation (egg->young nymph)") +) + +parser <- OptionParser(usage="%prog [options] file", option_list=option_list) +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) + namedat <- paste(loc, year, ".Rdat", sep="") + temp.data <- read.csv(file=temperature.dataset, header=T) + + expdata[,1] <- c(1:365) + # Minimum + expdata[,2] <- temp.data[c(1:365), 3] + # Maximum + expdata[,3] <- temp.data[c(1:365), 2] + save(expdata, file=namedat) + namedat +} + +daylength=function(latitude) +{ + # from Forsythe 1995 + p=0.8333 + dl <- NULL + for (i in 1:365) { + 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 +} + +hourtemp=function(latitude, date, temperature_file_path) +{ + load(temperature_file_path) + threshold <- 14.17 # base development threshold for BMSB + 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 + T <- NULL # initialize hourly temperature + dh <- NULL #initialize degree hour vector + # date <- 200 + y <- dlprofile[date] # calculate daylength in given date + z <- 24 - y # night length + a <- 1.86 # lag coefficient + b <- 2.20 # night coefficient + #tempdata <- read.csv("tempdata.csv") #import raw data set + # Should be outside function otherwise its redundant + risetime <- 12 - y / 2 # sunrise time + settime <- 12 + y / 2 # sunset time + 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 + T[i]=(dxp - dnp) * sin(pi * m / (y + 2 * a)) + dnp + if (T[i]<8.4) { + dh[i] <- 0 + } + else { + dh[i] <- T[i] - 8.4 + } + } + else if (i > settime) { + n <- i - settime + T[i]=dnp + (ts - dnp) * exp( - b * n / z) + if (T[i]<8.4) { + dh[i] <- 0 + } + else { + dh[i] <- T[i] - 8.4 + } + } + else { + n <- i + 24 - settime + T[i]=dnp + (ts - dnp) * exp( - b * n / z) + if (T[i]<8.4) { + dh[i] <- 0 + } + else { + dh[i] <- T[i] - 8.4 + } + } + } + dd <- sum(dh) / 24 + } + return=c(dmean, dd) + return +} + +dev.egg = function(temperature) +{ + dev.rate= -0.9843 * temperature + 33.438 + return = dev.rate + return +} + +dev.young = function(temperature) +{ + n12 <- -0.3728 * temperature + 14.68 + n23 <- -0.6119 * temperature + 25.249 + dev.rate = mean(n12 + n23) + return = dev.rate + return +} + +dev.old = function(temperature) +{ + n34 <- -0.6119 * temperature + 17.602 + n45 <- -0.4408 * temperature + 19.036 + dev.rate = mean(n34 + n45) + return = dev.rate + return +} + +dev.emerg = function(temperature) +{ + emerg.rate <- -0.5332 * temperature + 24.147 + return = emerg.rate + return +} + +mortality.egg = function(temperature) +{ + if (temperature < 12.7) { + mort.prob = 0.8 + } + else { + mort.prob = 0.8 - temperature / 40.0 + if (mort.prob < 0) { + mort.prob = 0.01 + } + } + return = mort.prob + return +} + +mortality.nymph = function(temperature) +{ + if (temperature < 12.7) { + mort.prob = 0.03 + } + else { + mort.prob = temperature * 0.0008 + 0.03 + } + return = mort.prob + return +} + +mortality.adult = function(temperature) +{ + if (temperature < 12.7) { + mort.prob = 0.002 + } + else { + mort.prob = temperature * 0.0005 + 0.02 + } + return = mort.prob + return +} + +cat("Replications: ", opt$replications, "\n") +cat("Photoperiod: ", opt$photoperiod, "\n") +cat("Oviposition rate: ", opt$oviposition, "\n") +cat("Egg mortality rate: ", opt$egg_mort, "\n") +cat("Nymph mortality rate: ", opt$nymph_mort, "\n") +cat("Adult mortality rate: ", opt$adult_mort, "\n") +cat("Min clutch size: ", opt$min_clutch_size, "\n") +cat("Max clutch size: ", opt$max_clutch_size, "\n") +cat("(egg->young nymph): ", opt$young_nymph_accum, "\n") +cat("(young nymph->old nymph): ", opt$old_nymph_accum, "\n") +cat("(old nymph->adult): ", opt$adult_accum, "\n") + +# 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) +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) + +# loop through replications +for (N.rep in 1:opt$replications) { + # during each replication + # start with 1000 individuals -- user definable as well? + n <- 1000 + # Generation, Stage, DD, T, Diapause + vec.ini <- c(0, 3, 0, 0, 0) + # 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) + + # time series of population size + tot.pop <- NULL + # gen.0 pop size + gen0.pop <- rep(0, 365) + gen1.pop <- rep(0, 365) + gen2.pop <- rep(0, 365) + S0 <- S1 <- S2 <- S3 <- S4 <- S5 <- rep(0, 365) + g0.adult <- g1.adult <- g2.adult <- rep(0, 365) + N.newborn <- N.death <- N.adult <- rep(0, 365) + dd.day <- rep(0, 365) + + # start tick + ptm <- proc.time() + + # all the days + for (day in 1:365) { + # photoperiod in the day + photoperiod <- ph.p[day] + temp.profile <- hourtemp(opt$latitude, day, temperature_file_path) + mean.temp <- temp.profile[1] + dd.temp <- temp.profile[2] + dd.day[day] <- dd.temp + # trash bin for death + death.vec <- NULL + # new born + birth.vec <- NULL + + # all individuals + for (i in 1:n) { + # find individual record + vec.ind <- vec.mat[i,] + # first of all, still alive? + # adjustment for late season mortality rate + if (opt$latitude < 40.0) { + post.mort <- 1 + day.kill <- 300 + } + else { + post.mort <- 2 + day.kill <- 250 + } + if (vec.ind[2] == 0) { + # egg + death.prob = opt$egg_mort * mortality.egg(mean.temp) + } + else if (vec.ind[2] == 1 | vec.ind[2] == 2) { + death.prob = opt$nymph_mort * mortality.nymph(mean.temp) + } + else if (vec.ind[2] == 3 | vec.ind[2] == 4 | vec.ind[2] == 5) { + # for adult + if (day < day.kill) { + death.prob = opt$adult_mort * mortality.adult(mean.temp) + } + else { + # increase adult mortality after fall equinox + death.prob = opt$adult_mort * post.mort * mortality.adult(mean.temp) + } + } + # (or dependent on temperature and life stage?) + u.d <- runif(1) + if (u.d < death.prob) { + death.vec <- c(death.vec, i) + } + else { + # aggregrate index of dead bug + # event 1 end of diapause + if (vec.ind[1] == 0 && vec.ind[2] == 3) { + # overwintering adult (previttelogenic) + if (photoperiod > opt$photoperiod && vec.ind[3] > 68 && day < 180) { + # add 68C to become fully reproductively matured + # transfer to vittelogenic + vec.ind <- c(0, 4, 0, 0, 0) + vec.mat[i,] <- vec.ind + } + else { + # add to DD + vec.ind[3] <- vec.ind[3] + dd.temp + # add 1 day in current stage + vec.ind[4] <- vec.ind[4] + 1 + vec.mat[i,] <- vec.ind + } + } + if (vec.ind[1] != 0 && vec.ind[2] == 3) { + # NOT overwintering adult (previttelogenic) + current.gen <- vec.ind[1] + if (vec.ind[3] > 68) { + # add 68C to become fully reproductively matured + # transfer to vittelogenic + vec.ind <- c(current.gen, 4, 0, 0, 0) + vec.mat[i,] <- vec.ind + } + else { + # add to DD + vec.ind[3] <- vec.ind[3] + dd.temp + # add 1 day in current stage + vec.ind[4] <- vec.ind[4] + 1 + vec.mat[i,] <- vec.ind + } + } + + # event 2 oviposition -- where population dynamics comes from + if (vec.ind[2] == 4 && vec.ind[1] == 0 && mean.temp > 10) { + # vittelogenic stage, overwintering generation + if (vec.ind[4] == 0) { + # just turned in vittelogenic stage + n.birth=round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size)) + } + else { + # daily probability of birth + p.birth = opt$oviposition * 0.01 + u1 <- runif(1) + if (u1 < p.birth) { + n.birth=round(runif(1, 2, 8)) + } + } + # add to DD + vec.ind[3] <- vec.ind[3] + dd.temp + # add 1 day in current stage + vec.ind[4] <- vec.ind[4] + 1 + vec.mat[i,] <- vec.ind + if (n.birth > 0) { + # add new birth -- might be in different generations + # generation + 1 + new.gen <- vec.ind[1] + 1 + # egg profile + new.ind <- c(new.gen, 0, 0, 0, 0) + new.vec <- rep(new.ind, n.birth) + # update batch of egg profile + new.vec <- t(matrix(new.vec, nrow=5)) + # group with total eggs laid in that day + birth.vec <- rbind(birth.vec, new.vec) + } + } + + # event 2 oviposition -- for gen 1. + if (vec.ind[2] == 4 && vec.ind[1] == 1 && mean.temp > 12.5 && day < 222) { + # vittelogenic stage, 1st generation + if (vec.ind[4] == 0) { + # just turned in vittelogenic stage + n.birth=round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size)) + } + else { + # daily probability of birth + p.birth = opt$oviposition * 0.01 + u1 <- runif(1) + if (u1 < p.birth) { + n.birth = round(runif(1, 2, 8)) + } + } + # add to DD + vec.ind[3] <- vec.ind[3] + dd.temp + # add 1 day in current stage + vec.ind[4] <- vec.ind[4] + 1 + vec.mat[i,] <- vec.ind + if (n.birth > 0) { + # add new birth -- might be in different generations + # generation + 1 + new.gen <- vec.ind[1] + 1 + # egg profile + new.ind <- c(new.gen, 0, 0, 0, 0) + new.vec <- rep(new.ind, n.birth) + # update batch of egg profile + new.vec <- t(matrix(new.vec, nrow=5)) + # group with total eggs laid in that day + birth.vec <- rbind(birth.vec, new.vec) + } + } + + # event 3 development (with diapause determination) + # event 3.1 egg development to young nymph (vec.ind[2]=0 -> egg) + if (vec.ind[2] == 0) { + # egg stage + # add to DD + vec.ind[3] <- vec.ind[3] + dd.temp + if (vec.ind[3] >= (68 + opt$young_nymph_accum)) { + # from egg to young nymph, DD requirement met + current.gen <- vec.ind[1] + # transfer to young nym stage + vec.ind <- c(current.gen, 1, 0, 0, 0) + } + else { + # add 1 day in current stage + vec.ind[4] <- vec.ind[4] + 1 + } + vec.mat[i,] <- vec.ind + } + + # event 3.2 young nymph to old nymph (vec.ind[2]=1 -> young nymph: determines diapause) + if (vec.ind[2] == 1) { + # young nymph stage + # add to DD + vec.ind[3] <- vec.ind[3] + dd.temp + if (vec.ind[3] >= (250 + opt$old_nymph_accum)) { + # from young to old nymph, DD requirement met + current.gen <- vec.ind[1] + # transfer to old nym stage + vec.ind <- c(current.gen, 2, 0, 0, 0) + if (photoperiod < opt$photoperiod && day > 180) { + vec.ind[5] <- 1 + } # prepare for diapausing + } + else { + # add 1 day in current stage + vec.ind[4] <- vec.ind[4] + 1 + } + vec.mat[i,] <- vec.ind + } + + # event 3.3 old nymph to adult: previttelogenic or diapausing? + if (vec.ind[2] == 2) { + # old nymph stage + # add to DD + vec.ind[3] <- vec.ind[3] + dd.temp + if (vec.ind[3] >= (200 + opt$adult_accum)) { + # from old to adult, DD requirement met + current.gen <- vec.ind[1] + if (vec.ind[5] == 0) { + # non-diapausing adult -- previttelogenic + vec.ind <- c(current.gen, 3, 0, 0, 0) + } + else { + # diapausing + vec.ind <- c(current.gen, 5, 0, 0, 1) + } + } + else { + # add 1 day in current stage + vec.ind[4] <- vec.ind[4] + 1 + } + vec.mat[i,] <- vec.ind + } + + # event 4 growing of diapausing adult (unimportant, but still necessary)## + if (vec.ind[2] == 5) { + vec.ind[3] <- vec.ind[3] + dd.temp + vec.ind[4] <- vec.ind[4] + 1 + vec.mat[i,] <- vec.ind + } + } # else if it is still alive + } # end of the individual bug loop + + # find how many died + n.death <- length(death.vec) + if (n.death > 0) { + vec.mat <- vec.mat[-death.vec, ] + } + # remove record of dead + # find how many new born + n.newborn <- length(birth.vec[,1]) + vec.mat <- rbind(vec.mat, birth.vec) + # update population size for the next day + n <- n - n.death + n.newborn + + # aggregate results by day + tot.pop <- c(tot.pop, n) + # egg + s0 <- sum(vec.mat[,2] == 0) + # young nymph + s1 <- sum(vec.mat[,2] == 1) + # old nymph + s2 <- sum(vec.mat[,2] == 2) + # previtellogenic + s3 <- sum(vec.mat[,2] == 3) + # vitellogenic + s4 <- sum(vec.mat[,2] == 4) + # diapausing + s5 <- sum(vec.mat[,2] == 5) + # overwintering adult + gen0 <- sum(vec.mat[,1] == 0) + # first generation + gen1 <- sum(vec.mat[,1] == 1) + # second generation + gen2 <- sum(vec.mat[,1] == 2) + # sum of all adults + n.adult <- sum(vec.mat[,2] == 3) + sum(vec.mat[,2] == 4) + sum(vec.mat[,2] == 5) + # gen.0 pop size + gen0.pop[day] <- gen0 + gen1.pop[day] <- gen1 + gen2.pop[day] <- gen2 + S0[day] <- s0 + S1[day] <- s1 + S2[day] <- s2 + S3[day] <- s3 + S4[day] <- s4 + S5[day] <- s5 + g0.adult[day] <- sum(vec.mat[,1] == 0) + g1.adult[day] <- sum((vec.mat[,1] == 1 & vec.mat[,2] == 3) | (vec.mat[,1] == 1 & vec.mat[,2] == 4) | (vec.mat[,1] == 1 & vec.mat[,2] == 5)) + 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 365 days + + dd.cum <- cumsum(dd.day) + # collect all the outputs + S0.rep[,N.rep] <- S0 + S1.rep[,N.rep] <- S1 + S2.rep[,N.rep] <- S2 + S3.rep[,N.rep] <- S3 + S4.rep[,N.rep] <- S4 + S5.rep[,N.rep] <- S5 + newborn.rep[,N.rep] <- N.newborn + death.rep[,N.rep] <- N.death + adult.rep[,N.rep] <- N.adult + pop.rep[,N.rep] <- tot.pop + g0.rep[,N.rep] <- gen0.pop + g1.rep[,N.rep] <- gen1.pop + g2.rep[,N.rep] <- gen2.pop + 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) + +# 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) +# mean value for eggs +se <- apply(S0.rep, 1, mean) +# mean value for P +g0 <- apply(g0.rep, 1, mean) +# mean value for F1 +g1 <- apply(g1.rep, 1, mean) +# mean value for F2 +g2 <- apply(g2.rep, 1, mean) +# mean value for P adult +g0a <- apply(g0a.rep, 1, mean) +# mean value for F1 adult +g1a <- apply(g1a.rep, 1, mean) +# mean value for F2 adult +g2a <- apply(g2a.rep, 1, mean) + +# SE for adults +sa.se <- apply((S3.rep + S4.rep + S5.rep), 1, sd) / sqrt(opt$replications) +# SE for nymphs +sn.se <- apply((S1.rep + S2.rep) / sqrt(opt$replications), 1, sd) +# SE for eggs +se.se <- apply(S0.rep, 1, sd) / sqrt(opt$replications) +# SE value for P +g0.se <- apply(g0.rep, 1, sd) / sqrt(opt$replications) +# SE for F1 +g1.se <- apply(g1.rep, 1, sd) / sqrt(opt$replications) +# SE for F2 +g2.se <- apply(g2.rep, 1, sd) / sqrt(opt$replications) +# SE for P adult +g0a.se <- apply(g0a.rep, 1, sd) / sqrt(opt$replications) +# SE for F1 adult +g1a.se <- apply(g1a.rep, 1, sd) / sqrt(opt$replications) +# SE for F2 adult +g2a.se <- apply(g2a.rep, 1, sd) / sqrt(opt$replications) + +dev.new(width=20, height=20) + +# Start PDF device driver to save charts to output. +pdf(file=opt$output, height=20, width=20, bg="white") + +par(mar = c(5, 6, 4, 4), mfrow=c(3, 1)) + +# Subfigure 2: population size by life stage +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) +# Young and old nymphs +lines(day.all, sn, lwd = 2, lty = 1, col = 2) +# Eggs +lines(day.all, se, lwd = 2, lty = 1, col = 4) +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")) +axis(2, cex.axis = 2) +leg.text <- c("Egg", "Nymph", "Adult") +legend("topleft", leg.text, lty = c(1, 1, 1), col = c(4, 2, 1), cex = 2) +if (opt$se_plot == 1) { + # add SE lines to plot + # SE for adults + lines (day.all, sa + sa.se, lty = 2) + lines (day.all, sa - sa.se, lty = 2) + # SE for nymphs + lines (day.all, sn + sn.se, col = 2, lty = 2) + lines (day.all, sn - sn.se, col = 2, lty = 2) + # SE for eggs + lines (day.all, se + se.se, col = 4, lty = 2) + lines (day.all, se - se.se, col = 4, lty = 2) +} + +# Subfigure 3: population size by generation +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) +lines(day.all, g1, lwd = 2, lty = 1, col = 2) +lines(day.all, g2, lwd = 2, lty = 1, col = 4) +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")) +axis(2, cex.axis = 2) +leg.text <- c("P", "F1", "F2") +legend("topleft", leg.text, lty = c(1, 1, 1), col =c(1, 2, 4), cex = 2) +if (opt$se_plot == 1) { + # add SE lines to plot + # SE for adults + lines (day.all, g0 + g0.se, lty = 2) + lines (day.all, g0 - g0.se, lty = 2) + # SE for nymphs + lines (day.all, g1 + g1.se, col = 2, lty = 2) + lines (day.all, g1 - g1.se, col = 2, lty = 2) + # SE for eggs + lines (day.all, g2 + g2.se, col = 4, lty = 2) + lines (day.all, g2 - g2.se, col = 4, lty = 2) +} + +# Subfigure 4: adult population size by generation +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) +lines(day.all, g1a, lwd = 2, lty = 1, col = 2) +lines(day.all, g2a, lwd = 2, lty = 1, col = 4) +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")) +axis(2, cex.axis = 2) +leg.text <- c("P", "F1", "F2") +legend("topleft", leg.text, lty = c(1, 1, 1), col = c(1, 2, 4), cex = 2) +if (opt$se_plot == 1) { + # add SE lines to plot + # SE for adults + lines (day.all, g0a + g0a.se, lty = 2) + lines (day.all, g0a - g0a.se, lty = 2) + # SE for nymphs + lines (day.all, g1a + g1a.se, col = 2, lty = 2) + lines (day.all, g1a - g1a.se, col = 2, lty = 2) + # SE for eggs + lines (day.all, g2a + g2a.se, col = 4, lty = 2) + lines (day.all, g2a - g2a.se, col = 4, lty = 2) +} + +# Turn off device driver to flush output. +dev.off()
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/insect_phenology_model.xml Tue Aug 08 13:14:39 2017 -0400 @@ -0,0 +1,101 @@ +<?xml version="."?> +<tool id="insect_phenology_model" name="Insect phenology model" version="1.0.0"> + <description>expressing stage-specific phenology and population dynamics</description> + <requirements> + <requirement type="package" version="1.3.0">r-optparse</requirement> + </requirements> + <command><![CDATA[ +#set loc_lat_items = str($location).split(':') +#set loc = $loc_lat_items[0] +#set lat = $loc_lat_items[1] +Rscript '$__tool_directory__/insect_phenology_model.R' +-a $adult_mort +-b $adult_accum +-c $egg_mort +-d $lat +-e $loc +-f $max_clutch_size +-i $min_clutch_size +-j $nymph_mort +-k $old_nymph_accum +-o $output +-p $oviposition +-q $photoperiod +-s $replications +-t $se_plot +-u $year +-v '$temperature_data' +-y $young_nymph_accum + ]]></command> + <inputs> + <param name="location" type="select" label="Location"> + <options from_file="locations.txt"> + <column name="name" index="0"/> + <column name="value" index="1"/> + </options> + </param> + <param name="temperature_data" type="data" format="csv" label="Temperature data" /> + <param name="year" type="select" label="Temperature data year"> + <options from_file="years.txt"> + <column name="name" index="0"/> + <column name="value" index="1"/> + </options> + </param> + <param name="replications" type="integer" value="10" min="1" label="Number of replications" /> + <param name="photoperiod" type="float" value="13.5" min="0" label="Critical photoperiod for diapause induction/termination" /> + <param name="egg_mort" type="integer" value="1" min="0" label="Adjustment rate for egg mortality" /> + <param name="nymph_mort" type="integer" value="1" min="0" label="Adjustment rate for nymph mortality" /> + <param name="adult_mort" type="integer" value="1" min="0" label="Adjustment rate for adult mortality" /> + <param name="oviposition" type="integer" value="1" min="0" label="Adjustment oviposition rate" /> + <param name="min_clutch_size" type="integer" value="0" min="0" label="Adjustment of minimum clutch size" /> + <param name="max_clutch_size" type="integer" value="0" min="0" label="Adjustment of maximum clutch size" /> + <param name="young_nymph_accum" type="integer" value="0" min="0" label="Adjustment of DD accumulation (egg->young nymph)" /> + <param name="old_nymph_accum" type="integer" value="0" min="0" label="Adjustment of DD accumulation (young nymph->old nymph)" /> + <param name="adult_accum" type="integer" value="0" min="0" label="Adjustment of DD accumulation (old nymph->adult)" /> + <param name="se_plot" type="select" label="Plot SE?"> + <option value="1" selected="True">Yes</option> + <option value="0">No</option> + </param> + </inputs> + <outputs> + <data name="output" format="pdf" label="${tool.name} ${str($location).split(':')[0]}, ${year} on ${on_string}" /> + </outputs> + <tests> + <test> + <param name="temperature_data" value="asheville2014.csv" ftype="csv" /> + <param name="location" value="asheville:35.58" /> + <param name="year" value="2014" /> + <param name="replications" value="3" /> + <output name="output" file="output.pdf" ftype="pdf" compare="contains" /> + </test> + </tests> + <help> +**What it does** + +Provides an agent-based stochastic model expressing stage-specific phenology and population dynamics for an insect species across geographic regions. + +----- + +**Required options** + + * **Location** - the location associated with the selected temperature data. + * **Temperature data** - select the dataset from your history containing the temperature data. + * **Temperature data year** - the year during which the temperature data was recorded. + * **Number of replications** - number of replications. + * **Critical photoperiod for diapause induction/termination** - critical photoperiod for diapause induction/termination. + * **Adjustment rate for egg mortality** - adjustment rate for egg mortality. + * **Adjustment rate for nymph mortality** - adjustment rate for nymph mortality. + * **Adjustment rate for adult mortality** - adjustment rate for adult mortality. + * **Adjustment oviposition rate** - adjustment oviposition rate. + * **Adjustment of minimum clutch size** - adjustment of minimum clutch size. + * **Adjustment of maximum clutch size** - adjustment of maximum clutch size + * **Adjustment of DD accumulation (egg->young nymph)** - adjustment of DD accumulation (egg->young nymph). + * **Adjustment of DD accumulation (young nymph->old nymph)** - adjustment of DD accumulation (young nymph->old nymph). + * **Adjustment of DD accumulation (old nymph->adult)** - adjustment of DD accumulation (old nymph->adult). + * **Plot SE** - add SE lines to plot for eggs, nymphs and adults. + + </help> + <citations> + <citation type="doi">10.3389/fphys.2016.00165</citation> + </citations> +</tool>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/asheville2014.csv Tue Aug 08 13:14:39 2017 -0400 @@ -0,0 +1,366 @@ +"","TMAX","TMIN" +"100746",3.9,-1.7 +"100747",7.8,-3.3 +"100748",7.8,-7.8 +"100749",-3.9,-11.1 +"100750",6.7,-6.7 +"100751",11.1,5 +"100752",8.3,-17.2 +"100753",-6.7,-18.3 +"100754",8.3,-12.8 +"100755",10.6,-5 +"100756",12.2,1.7 +"100757",15,6.1 +"100758",12.8,3.3 +"100759",12.8,-0.6 +"100760",10.6,0.6 +"100761",3.3,-3.3 +"100762",6.1,-6.7 +"100763",9.4,-3.9 +"100764",2.8,-8.3 +"100765",7.2,-2.8 +"100766",16.7,-3.3 +"100767",6.1,-8.9 +"100768",-2.8,-12.8 +"100769",-2.8,-12.8 +"100770",-2.2,-16.7 +"100771",3.3,-10 +"100772",12.2,-3.9 +"100773",9.4,-6.1 +"100774",-5.6,-9.4 +"100775",-3.3,-13.3 +"100776",3.9,-14.4 +"100777",10.6,-7.2 +"100778",16.1,-3.9 +"100779",17.8,-1.1 +"100780",15,2.8 +"100781",7.2,2.8 +"100782",12.8,-2.8 +"100783",0.6,-4.4 +"100784",5.6,-3.9 +"100785",10,0 +"100786",11.7,-1.1 +"100787",3.3,-1.1 +"100788",1.1,-1.1 +"100789",-1.1,-5 +"100790",0.6,-3.3 +"100791",12.8,-2.2 +"100792",5,-3.3 +"100793",11.1,-3.9 +"100794",8.9,-1.1 +"100795",17.2,2.2 +"100796",20,7.2 +"100797",20.6,7.8 +"100798",17.2,3.3 +"100799",17.8,3.3 +"100800",20.6,2.2 +"100801",10,2.2 +"100802",10.6,0.6 +"100803",3.3,-2.8 +"100804",5.6,-6.1 +"100805",8.3,-5 +"100806",14.4,1.1 +"100807",19.4,1.1 +"100808",14.4,-2.2 +"100809",7.2,-2.8 +"100810",15.6,-0.6 +"100811",6.1,0.6 +"100812",12.2,1.1 +"100813",18.3,1.1 +"100814",16.1,7.2 +"100815",20,1.7 +"100816",24.4,3.9 +"100817",20.6,-1.7 +"100818",3.9,-3.3 +"100819",15,-3.9 +"100820",17.8,7.8 +"100821",10.6,5 +"100822",5,0 +"100823",10.6,0 +"100824",13.3,4.4 +"100825",13.9,6.7 +"100826",17.8,0 +"100827",20.6,7.2 +"100828",12.2,3.3 +"100829",11.7,-1.7 +"100830",5.6,-2.8 +"100831",6.1,-5 +"100832",14.4,-2.8 +"100833",14.4,9.4 +"100834",17.2,2.8 +"100835",10.6,0 +"100836",23.9,0.6 +"100837",26.1,7.2 +"100838",27.8,7.2 +"100839",27.2,12.2 +"100840",21.7,13.9 +"100841",16.7,5.6 +"100842",15.6,7.2 +"100843",13.3,6.1 +"100844",13.3,7.2 +"100845",15,5.6 +"100846",21.7,1.7 +"100847",23.9,7.2 +"100848",25,10 +"100849",25,13.3 +"100850",18.9,15.6 +"100851",16.1,0.6 +"100852",14.4,-2.2 +"100853",16.7,2.2 +"100854",11.7,5.6 +"100855",12.8,7.2 +"100856",22.2,5.6 +"100857",24.4,5.6 +"100858",18.9,11.7 +"100859",19.4,9.4 +"100860",22.8,5 +"100861",23.9,13.3 +"100862",26.1,8.3 +"100863",27.8,9.4 +"100864",25.6,16.1 +"100865",22.8,16.1 +"100866",21.7,15 +"100867",18.3,11.7 +"100868",20.6,8.3 +"100869",21.1,7.8 +"100870",26.7,11.1 +"100871",31.7,12.8 +"100872",30.6,13.3 +"100873",29.4,11.7 +"100874",30,12.2 +"100875",24.4,15 +"100876",21.1,14.4 +"100877",29.4,12.8 +"100878",28.9,15 +"100879",30,16.7 +"100880",26.7,15.6 +"100881",20,8.9 +"100882",16.7,6.1 +"100883",16.7,3.9 +"100884",16.1,7.8 +"100885",21.1,6.7 +"100886",25.6,10 +"100887",28.3,12.2 +"100888",27.8,15.6 +"100889",26.7,15 +"100890",26.1,12.2 +"100891",25.6,10.6 +"100892",28.3,14.4 +"100893",27.8,16.1 +"100894",28.9,15 +"100895",28.3,15.6 +"100896",27.8,17.8 +"100897",25.6,17.8 +"100898",22.8,16.7 +"100899",25.6,13.3 +"100900",26.1,13.9 +"100901",30,14.4 +"100902",27.2,17.8 +"100903",27.8,17.8 +"100904",28.9,17.2 +"100905",28.3,20 +"100906",28.3,15.6 +"100907",31.1,16.7 +"100908",27.2,17.8 +"100909",25,17.2 +"100910",27.2,16.7 +"100911",28.3,15.6 +"100912",28.9,16.1 +"100913",28.9,20 +"100914",31.7,17.2 +"100915",32.2,18.3 +"100916",31.7,17.8 +"100917",31.1,16.7 +"100918",29.4,18.3 +"100919",30,17.2 +"100920",28.3,20 +"100921",26.7,18.9 +"100922",27.2,18.9 +"100923",30,18.3 +"100924",28.3,20 +"100925",26.7,20 +"100926",29.4,19.4 +"100927",28.3,18.3 +"100928",32.2,19.4 +"100929",32.2,20.6 +"100930",30,18.3 +"100931",25.6,15 +"100932",27.2,11.7 +"100933",28.3,13.3 +"100934",30,15.6 +"100935",30.6,17.2 +"100936",26.7,18.3 +"100937",27.8,17.8 +"100938",29.4,18.3 +"100939",30,18.3 +"100940",30.6,19.4 +"100941",31.7,19.4 +"100942",26.1,18.3 +"100943",25,13.9 +"100944",26.1,13.3 +"100945",20.6,17.2 +"100946",18.9,16.7 +"100947",26.7,16.7 +"100948",25.6,18.9 +"100949",27.8,19.4 +"100950",28.3,18.9 +"100951",28.3,18.9 +"100952",28.9,17.8 +"100953",30.6,18.3 +"100954",31.1,19.4 +"100955",25,17.2 +"100956",25,15 +"100957",26.1,11.1 +"100958",20.6,13.9 +"100959",23.9,16.1 +"100960",26.7,17.2 +"100961",27.8,18.3 +"100962",26.1,16.1 +"100963",28.9,17.2 +"100964",29.4,15.6 +"100965",27.8,17.2 +"100966",23.9,18.9 +"100967",28.3,19.4 +"100968",23.9,17.8 +"100969",28.3,18.9 +"100970",27.8,19.4 +"100971",22.8,16.1 +"100972",27.2,12.8 +"100973",27.8,14.4 +"100974",27.2,15 +"100975",29.4,17.8 +"100976",28.9,20 +"100977",28.9,19.4 +"100978",32.2,17.8 +"100979",31.1,17.2 +"100980",31.7,19.4 +"100981",31.1,20 +"100982",26.1,20 +"100983",27.8,18.3 +"100984",28.9,14.4 +"100985",31.1,14.4 +"100986",32.8,16.1 +"100987",31.1,18.3 +"100988",30,19.4 +"100989",31.7,20 +"100990",32.2,20 +"100991",31.7,18.9 +"100992",31.7,17.8 +"100993",29.4,19.4 +"100994",31.1,18.9 +"100995",30.6,19.4 +"100996",30,18.3 +"100997",24.4,18.3 +"100998",29.4,17.8 +"100999",26.1,18.9 +"101000",29.4,17.8 +"101001",31.1,18.9 +"101002",25,17.2 +"101003",21.1,16.1 +"101004",27.2,16.1 +"101005",27.2,17.8 +"101006",25,16.7 +"101007",25.6,14.4 +"101008",24.4,15.6 +"101009",26.1,16.1 +"101010",28.3,12.8 +"101011",21.1,12.2 +"101012",21.1,8.9 +"101013",23.3,10 +"101014",23.9,15.6 +"101015",23.9,17.2 +"101016",24.4,15.6 +"101017",24.4,13.3 +"101018",19.4,16.1 +"101019",26.1,14.4 +"101020",27.8,11.1 +"101021",27.8,12.2 +"101022",23.3,16.1 +"101023",16.1,6.7 +"101024",18.3,1.1 +"101025",20,5 +"101026",20,12.2 +"101027",25,13.9 +"101028",27.2,10.6 +"101029",28.3,13.3 +"101030",25,16.7 +"101031",17.8,13.3 +"101032",21.1,12.8 +"101033",20,15 +"101034",19.4,10 +"101035",15.6,10 +"101036",24.4,7.2 +"101037",17.2,9.4 +"101038",18.9,7.8 +"101039",22.2,5 +"101040",17.8,7.2 +"101041",13.3,5.6 +"101042",16.7,5 +"101043",19.4,3.9 +"101044",21.1,3.3 +"101045",22.8,8.9 +"101046",26.7,7.2 +"101047",26.7,8.9 +"101048",16.7,7.2 +"101049",13.3,2.2 +"101050",13.9,2.2 +"101051",2.2,-0.6 +"101052",6.1,-1.1 +"101053",19.4,-2.8 +"101054",22.2,2.2 +"101055",17.8,6.7 +"101056",20.6,6.7 +"101057",8.3,1.7 +"101058",15,-1.1 +"101059",17.2,1.7 +"101060",18.3,1.1 +"101061",22.2,2.2 +"101062",17.8,6.1 +"101063",6.1,-1.7 +"101064",-1.1,-6.1 +"101065",6.7,-7.8 +"101066",12.8,1.7 +"101067",12.2,-5 +"101068",-2.8,-8.3 +"101069",6.1,-10 +"101070",11.1,-1.7 +"101071",10.6,0 +"101072",13.9,-1.7 +"101073",12.8,5 +"101074",22.2,11.7 +"101075",13.3,6.1 +"101076",10,3.3 +"101077",3.9,-2.2 +"101078",5,-3.9 +"101079",15.6,-2.8 +"101080",20,1.1 +"101081",21.1,4.4 +"101082",13.9,7.8 +"101083",16.7,10.6 +"101084",14.4,6.7 +"101085",11.1,7.2 +"101086",14.4,5.6 +"101087",10.6,2.2 +"101088",5,0.6 +"101089",4.4,0 +"101090",2.2,-0.6 +"101091",3.9,-2.2 +"101092",10,-4.4 +"101093",12.2,-1.7 +"101094",9.4,-1.7 +"101095",15.6,-3.9 +"101096",16.1,6.1 +"101097",8.3,2.8 +"101098",11.7,0.6 +"101099",10.6,2.8 +"101100",7.8,3.9 +"101101",9.4,3.9 +"101102",8.9,3.9 +"101103",11.1,3.9 +"101104",12.8,6.1 +"101105",8.3,1.7 +"101106",15.6,-2.2 +"101107",16.7,2.8 +"101108",14.4,10 +"101109",13.9,6.7 +"101110",6.7,0.6
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/output.pdf Tue Aug 08 13:14:39 2017 -0400 @@ -0,0 +1,51 @@ +%PDF-1.4 +%ρ\r +1 0 obj +<< +/Title (R Graphics Output) +/Creator (R) +>> +endobj +2 0 obj +<< /Type /Catalog /Pages 3 0 R >> +endobj +7 0 obj +<< /Type /Page /Parent 3 0 R /Contents 8 0 R /Resources 4 0 R >> +endobj +8 0 obj +endobj +3 0 obj +<< /Type /Pages /Kids [ 7 0 R ] /Count 1 /MediaBox [0 0 1440 1440] >> +endobj +4 0 obj +<< +/ProcSet [/PDF /Text] +/Font <</F2 10 0 R /F3 11 0 R >> +/ExtGState << >> +/ColorSpace << /sRGB 5 0 R >> +>> +endobj +5 0 obj +[/ICCBased 6 0 R] +endobj +6 0 obj +/Type /Encoding /BaseEncoding /WinAnsiEncoding +/Differences [ 45/minus 96/quoteleft +144/dotlessi /grave /acute /circumflex /tilde /macron /breve /dotaccent +/dieresis /.notdef /ring /cedilla /.notdef /hungarumlaut /ogonek /caron /space] +>> +endobj +10 0 obj +<< /Type /Font /Subtype /Type1 /Name /F2 /BaseFont /Helvetica +/Encoding 9 0 R >> +endobj +11 0 obj +<< /Type /Font /Subtype /Type1 /Name /F3 /BaseFont /Helvetica-Bold +/Encoding 9 0 R >> +endobj +xref +0 12 +trailer +<< /Size 12 /Info 1 0 R /Root 2 0 R >> +startxref +%%EOF
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tool-data/locations.txt.sample Tue Aug 08 13:14:39 2017 -0400 @@ -0,0 +1,8 @@ +Asheville NC asheville:35.58 +Bridgeton NJ bridgeton:39.43 +Davis CA davis:38.55 +Geneva NY geneva:42.88 +Homestead FL homestead:25.47 +Riverside CA riverside:33.95 +Salem OR salem:44.93 +Wenatchee WA wneatchee:47.42
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tool-data/years.txt.sample Tue Aug 08 13:14:39 2017 -0400 @@ -0,0 +1,31 @@ +1995 1995 +1996 1996 +1997 1997 +1998 1998 +1999 1999 +2000 2000 +2001 2001 +2002 2002 +2003 2003 +2004 2004 +2005 2005 +2006 2006 +2007 2007 +2008 2008 +2009 2009 +2010 2010 +2011 2011 +2012 2012 +2013 2013 +2014 2014 +2015 2015 +2016 2016 +2017 2017 +2018 2018 +2019 2019 +2020 2020 +2021 2021 +2022 2022 +2023 2023 +2024 2024 +2025 2025