bmsb: bmsb.R comparison

comparison bmsb.R @ 0:ff341ead2c11 draft

Uploaded

author	greg
date	Mon, 01 Aug 2016 11:16:51 -0400
parents
children	99a386ac1f5b

comparison

equal deleted inserted replaced

--1:000000000000
+:ff341ead2c11
+#!/usr/bin/env Rscript
+suppressPackageStartupMessages(require("hash"))
+suppressPackageStartupMessages(require("KernSmooth"))
+suppressPackageStartupMessages(require("methods"))
+suppressPackageStartupMessages(library("optparse"))
+suppressPackageStartupMessages(library("rjson"))
+suppressPackageStartupMessages(require("utils"))
+option_list <- list(
+make_option(c("-s", "--save_log"), action="store_true", default=FALSE, help="Save R logs"),
+make_option(c("-m", "--output_r_log"), action="store", help="Output dataset for R logs"),
+make_option(c("-o", "--output"), action="store", help="Output dataset")
+)
+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: setwd(“/home/lunarmouse/Nielsen's project/”)
+# 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="asheville2014sim.Rdat")

Mercurial > repos > greg > bmsb

comparison bmsb.R @ 0:ff341ead2c11 draft