Mercurial > repos > greg > bmsb
changeset 19:d965e188feab draft
Uploaded
| author | greg |
|---|---|
| date | Tue, 16 Aug 2016 13:05:09 -0400 |
| parents | c6668285a216 |
| children | b477344d20f1 |
| files | bmsb.R |
| diffstat | 1 files changed, 432 insertions(+), 4 deletions(-) [+] |
line wrap: on
line diff
--- a/bmsb.R Tue Aug 16 11:38:43 2016 -0400 +++ b/bmsb.R Tue Aug 16 13:05:09 2016 -0400 @@ -1,7 +1,11 @@ #!/usr/bin/env Rscript +#suppressPackageStartupMessages(library("optparse")) + options_list <- list( - make_option(c("-o", "--output"), action="store", help="Output dataset") + make_option(c("-s", "--save_log"), action="store_true", default=FALSE, help="Save R logs"), + make_option(c("-m", "--output_r_log"), help="Output dataset for R logs"), + make_option(c("-o", "--output"), help="Output dataset") ) parser <- OptionParser(usage="%prog [options] file", options_list) @@ -9,6 +13,430 @@ opt <- args$options -fileConn<-file(opt$output) -writeLines(c("Hello World!"), fileConn) -close(fileConn) +daylength = function(L){ + # 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(L * pi/180) * sin(phi))/(cos(L * pi/180) * cos(phi))) + } + # return a vector of daylength in 365 days + dl +} + + +# source("daylength.R") +hourtemp = function(L,date){ + # L = 37.5 specify this in main program + # base development threshold for BMSB + threshold <- 12.7 + # threshold2 <- threshold/24 degree hour accumulation + #expdata <- tempdata[1:365,11:13] # Use daily max, min, mean + # daily minimum + dnp <- expdata[date,2] + # daily maximum + dxp <- expdata[date,3] + dmean <- 0.5 * (dnp + dxp) + #if (dmean>0) { + #dnp <- dnp - k1 * dmean + #dxp <- dxp + k2 * dmean + #} else { + #dnp <- dnp + k1 * dmean + #dxp <- dxp - k2 * dmean + #} + dd <- 0 # initialize degree day accumulation + + if (dxp<threshold) { + dd <- 0 + } else { + # extract daylength data for entire year + dlprofile <- daylength(L) + # initialize hourly temperature + T <- NULL + #initialize degree hour vector + dh <- NULL + # calculate daylength in given date + # date <- 200 + y <- dlprofile[date] + # night length + z <- 24 - y + # lag coefficient + a <- 1.86 + # night coefficient + b <- 2.20 + #import raw data set + #tempdata <- read.csv("tempdata.csv") + # Should be outside function otherwise its redundant + # 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) { + # 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 { + 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 +} + + +mortality.egg = function(temperature) { + if (temperature < 12.7) { + mort.prob = 1 + } else { + # 100% mortality if <12.7 + mort.prob = 0.8 - temperature / 40 + 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 { + # at low temperature + 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 +} + + +# model initialization +# TODO: add tool params for the following options. +# start with 1000 individuals +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)) +# latitude for Asheville NC +L <- 35.58 +# complete photoperiod profile in a year, requires daylength function +ph.p <- daylength(L) + +# load temperature data@location/year +load("asheville2014.Rdat") + +# 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) + +# aggregate +S0 <- S1 <- S2 <- S3 <- S4 <- S5 <- rep(0, 365) +g0.adult <- g1.adult <- g2.adult <- rep(0, 365) + +# birth death adults +N.newborn <- N.death <- N.adult <- rep(0, 365) + +# degree-day accumulation +dd.day <- rep(0, 365) + +# start tick +ptm <- proc.time() + +for (n.sim in 1:1000) { + # loop through 1000 simulations + for (day in 1:365) { + # loop through 365 day/yr + photoperiod <- ph.p[day] + # photoperiod in the day + temp.profile <- hourtemp(L,day) + # temperature profile + mean.temp <- temp.profile[1] + # mean temp + dd.temp <- temp.profile[2] + # degree-day + dd.day[day] <- dd.temp + death.vec <- NULL + # trash bin for death + birth.vec <- NULL + # record new born + for (i in 1:n) { + # loop through all individual + vec.ind <- vec.mat[i,] + # find individual record + # first of all, still alive? + if (vec.ind[2] == 0) { + # egg + death.prob = mortality.egg(mean.temp) + } else if (vec.ind[2] == 1 | vec.ind[2] == 2) { + # nymph + death.prob = mortality.nymph(mean.temp) + } else if (vec.ind[2] == 3 | vec.ind[2] == 4 | vec.ind[2] == 5) { + # for adult + death.prob = mortality.adult(mean.temp) + } + 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>13.5 && 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 + vec.ind[4] <- vec.ind[4] + 1 # add 1 day in current stage + 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 + # vittelogenic stage, overwintering generation + if (vec.ind[2] == 4 && vec.ind[1] == 0 && mean.temp>10) { + if (vec.ind[4] == 0) { + # just turned in vittelogenic stage + n.birth = round(runif(1,10,20)) + } else { + p.birth = 1/4/75 + # prob of birth + u1 <- runif(1) + if (u1<p.birth) { + n.birth = n.birth + } + } + # 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,10,20)) + } else { + p.birth = 1/4/75 + u1 <- runif(1) + if (u1<p.birth) { + n.birth = n.birth + } + } + # 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) + # egg stage + if (vec.ind[2] == 0) { + # add to DD + vec.ind[3] <- vec.ind[3] + dd.temp + # from egg to young nymph + if (vec.ind[3] >= 53.30 && -0.9843 * dd.temp + 33.438>0) { + 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) + # young nymph stage + if (vec.ind[2] == 1) { + # add to DD + vec.ind[3] <- vec.ind[3] + dd.temp + # from young to old nymph + if (vec.ind[3] >= 537/2 && -0.45 * dd.temp + 18>0) { + current.gen <- vec.ind[1] + # transfer to old nym stage + vec.ind <- c(current.gen,2,0,0,0) + # prepare for diapausing + if (photoperiod<13.5 && day > 180) { + vec.ind[5] <- 1 + } + } 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? + # old nymph stage + if (vec.ind[2] == 2) { + # add to DD + vec.ind[3] <- vec.ind[3] + dd.temp + # from old to adult + if (vec.ind[3] >= 537/2 && -0.50 * dd.temp + 22>0) { + current.gen <- vec.ind[1] + # non-diapausing adult -- previttelogenic + if (vec.ind[5] == 0) { + vec.ind <- c(current.gen,3,0,0,0) + # diapausing + } else { + 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(n.sim) +} + +proc.time() - ptm +dd.cum <- cumsum(dd.day) +save(dd.day, dd.cum, S0, S1, S2, S3, S4, S5, N.newborn, N.death, N.adult, tot.pop, gen0.pop, gen1.pop, gen2.pop, g0.adult, g1.adult, g2.adult, file=opt$output)