Mercurial > repos > greg > bmsb
changeset 18:c6668285a216 draft
Uploaded
| author | greg |
|---|---|
| date | Tue, 16 Aug 2016 11:38:43 -0400 |
| parents | f6977dcfbec3 |
| children | d965e188feab |
| files | bmsb.R |
| diffstat | 1 files changed, 4 insertions(+), 430 deletions(-) [+] |
line wrap: on
line diff
--- a/bmsb.R Tue Aug 16 11:36:22 2016 -0400 +++ b/bmsb.R Tue Aug 16 11:38:43 2016 -0400 @@ -1,9 +1,7 @@ #!/usr/bin/env Rscript -#suppressPackageStartupMessages(library("optparse")) - options_list <- list( - make_option(c("-o", "--output"), help="Output dataset") + make_option(c("-o", "--output"), action="store", help="Output dataset") ) parser <- OptionParser(usage="%prog [options] file", options_list) @@ -11,430 +9,6 @@ opt <- args$options -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) +fileConn<-file(opt$output) +writeLines(c("Hello World!"), fileConn) +close(fileConn)