Mercurial > repos > greg > bmsb
changeset 34:7c43bc0ce69c draft
Deleted selected files
| author | greg |
|---|---|
| date | Fri, 16 Dec 2016 08:51:06 -0500 |
| parents | 390ed5192839 |
| children | 7c40c2b303f1 |
| files | BMSB.R |
| diffstat | 1 files changed, 0 insertions(+), 578 deletions(-) [+] |
line wrap: on
line diff
--- a/BMSB.R Fri Dec 16 08:50:54 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,578 +0,0 @@ -#!/usr/bin/env Rscript - -suppressPackageStartupMessages(library("optparse")) - -options_list <- list( - make_option(c("-adult_mort", "--adult_mort"), action="store", help="Adjustment rate for adult mortality"), - make_option(c("-adult_nymph_accum", "--adult_nymph_accum"), action="store", help="Adjustment of DD accumulation (old nymph->adult)"), - make_option(c("-egg_mort", "--egg_mort"), action="store", help="Adjustment rate for egg mortality"), - make_option(c("-latitude", "--latitude"), action="store", help="Latitude of selected location"), - make_option(c("-location", "--location"), action="store", help="Selected location"), - make_option(c("-min_clutch_size", "--min_clutch_size"), action="store", help="Adjustment of minimum clutch size"), - make_option(c("-max_clutch_size", "--max_clutch_size"), action="store", help="Adjustment of maximum clutch size"), - make_option(c("-nymph_mort", "--nymph_mort"), action="store", help="Adjustment rate for nymph mortality"), - make_option(c("-old_nymph_accum", "--old_nymph_accum"), action="store", help="Adjustment of DD accumulation (young nymph->old nymph)"), - make_option(c("-output", "--output"), action="store", help="Output dataset"), - make_option(c("-oviposition", "--oviposition"), action="store", help="Adjustment for oviposition rate"), - make_option(c("-photoperiod", "--photoperiod"), action="store", help="Critical photoperiod for diapause induction/termination"), - make_option(c("-replications", "--replications"), action="store", help="Number of replications"), - make_option(c("-se_plot", "--se_plot"), action="store", help="Plot SE"), - make_option(c("-start_year", "--start_year"), action="store", help="Starting year"), - make_option(c("-sim_year", "--sim_year"), action="store", help="Simulation year"), - make_option(c("-temperature_datasets", "--temperature_datasets"), action="store", help="Temperature data for selected location"), - make_option(c("-young_nymph_accum", "--young_nymph_accum"), action="store", help="Adjustment of DD accumulation (egg->young nymph)"), -) - -parser <- OptionParser(usage="%prog [options] file", options_list) -arguements <- parse_args(parser, positional_arguments=TRUE) -opt <- args$options -args <- arguments$args - -temperature_datasets <- strsplit(opt$temperature_datasets, ",") - -# read in the input temperature datafile -data.input(opt$location, opt$start_year, temperature_datasets) - -input.name<-paste(opt$location, opt$sim_year, ".Rdat" ,sep="") -output.name<-paste(opt$location, opt$sim_year, "sim.Rdat", sep="") -load(input.name) - -# initialize matrix for results from all replications -S0.rep<-S1.rep<-S2.rep<-S3.rep<-S4.rep<-S5.rep<-matrix(rep(0,365*n.rep),ncol=n.rep) -newborn.rep<-death.rep<-adult.rep<-pop.rep<-g0.rep<-g1.rep<-g2.rep<-g0a.rep<-g1a.rep<-g2a.rep<-matrix(rep(0,365*n.rep),ncol=n.rep) - -# loop through replications -for (N.rep in 1:n.rep) -{ - # during each replication - n<-1000 # start with 1000 individuals -- user definable as well? - # 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) - vec.mat<-t(matrix(vec.mat,nrow=5)) # complete matrix for the population - ph.p<-daylength(L) # complete photoperiod profile in a year, requires daylength function - - tot.pop<-NULL # time series of population size - gen0.pop<-rep(0,365) # gen.0 pop size - 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) - - ptm <- proc.time() # start tick - - # all the days - for (day in 1:365) - { - photoperiod<-ph.p[day] # photoperiod in the day - temp.profile<-hourtemp(L,day) - mean.temp<-temp.profile[1] - dd.temp<-temp.profile[2] - dd.day[day]<-dd.temp - death.vec<-NULL # trash bin for death - birth.vec<-NULL # new born - #n<-length(vec.mat[,1]) # population size at previous day - - # all individuals - for (i in 1:n) - { - vec.ind<-vec.mat[i,] # find individual record - # first of all, still alive? - # adjustment for late season mortality rate - if (L<40) - { - post.mort<-1 - day.kill<-300 - } - else - { - post.mort<-2 - day.kill<-250 - } - # egg - if(vec.ind[2]==0) - { - death.prob=ar.em*mortality.egg(mean.temp) - } - else if (vec.ind[2]==1 | vec.ind[2]==2) - { - death.prob=ar.nm*mortality.nymph(mean.temp) - } - # for adult - else if (vec.ind[2]==3 | vec.ind[2]==4 | vec.ind[2]==5) - { - if (day<day.kill) - { - death.prob=ar.am*mortality.adult(mean.temp) - } - else - { - death.prob=ar.am*post.mort*mortality.adult(mean.temp)} # increase adult mortality after fall equinox - } - #(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>ph.cr && vec.ind[3]>68 && day<180) - { - # add 68C to become fully reproductively matured - vec.ind<-c(0,4,0,0,0) # transfer to vittelogenic - vec.mat[i,]<-vec.ind - } - else - { - vec.ind[3]<-vec.ind[3]+dd.temp # add to DD - 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 - vec.ind<-c(current.gen,4,0,0,0) # transfer to vittelogenic - vec.mat[i,]<-vec.ind - } - else - { - vec.ind[3]<-vec.ind[3]+dd.temp # add to DD - vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage - 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+min.ovi.adj,8+max.ovi.adj)) - } - else - { - p.birth=ar.ovi*0.01 # daily probability of birth - u1<-runif(1) - if (u1<p.birth) - { - n.birth=round(runif(1,2,8)) - } - } - vec.ind[3]<-vec.ind[3]+dd.temp # add to DD - vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage - vec.mat[i,]<-vec.ind - if (n.birth>0) - { - # add new birth -- might be in different generations - new.gen<-vec.ind[1]+1 # generation +1 - new.ind<-c(new.gen,0,0,0,0) # egg profile - new.vec<-rep(new.ind,n.birth) - new.vec<-t(matrix(new.vec,nrow=5)) # update batch of egg profile - birth.vec<-rbind(birth.vec,new.vec) # group with total eggs laid in that day - } - } - - # 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+min.ovi.adj,8+max.ovi.adj)) - } - else - { - p.birth=ar.ovi*0.01 # daily probability of birth - u1<-runif(1) - if (u1<p.birth) - { - n.birth=round(runif(1,2,8)) - } - } - vec.ind[3]<-vec.ind[3]+dd.temp # add to DD - vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage - vec.mat[i,]<-vec.ind - if (n.birth>0) - { - # add new birth -- might be in different generations - new.gen<-vec.ind[1]+1 # generation +1 - new.ind<-c(new.gen,0,0,0,0) # egg profile - new.vec<-rep(new.ind,n.birth) - new.vec<-t(matrix(new.vec,nrow=5)) # update batch of egg profile - birth.vec<-rbind(birth.vec,new.vec) # group with total eggs laid in that day - } - } - - # 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 - vec.ind[3]<-vec.ind[3]+dd.temp # add to DD - if (vec.ind[3]>=(68+dd.adj1)) - { - # from egg to young nymph, DD requirement met - current.gen<-vec.ind[1] - vec.ind<-c(current.gen,1,0,0,0) # transfer to young nym stage - } - else - { - vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage - } - 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 - vec.ind[3]<-vec.ind[3]+dd.temp # add to DD - if (vec.ind[3]>=(250+dd.adj2)) - { - # from young to old nymph, DD requirement met - current.gen<-vec.ind[1] - vec.ind<-c(current.gen,2,0,0,0) # transfer to old nym stage - if (photoperiod<ph.cr && day > 180) - { - vec.ind[5]<-1 - } # prepare for diapausing - } - else - { - vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage - } - vec.mat[i,]<-vec.ind - } - - # event 3.3 old nymph to adult: previttelogenic or diapausing? - if (vec.ind[2]==2) - { - # old nymph stage - vec.ind[3]<-vec.ind[3]+dd.temp # add to DD - if (vec.ind[3]>=(200+dd.adj3)) - { - # 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 - { - vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage - } - 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) - s0<-sum(vec.mat[,2]==0) #egg - s1<-sum(vec.mat[,2]==1) # young nymph - s2<-sum(vec.mat[,2]==2) # old nymph - s3<-sum(vec.mat[,2]==3) # previtellogenic - s4<-sum(vec.mat[,2]==4) # vitellogenic - s5<-sum(vec.mat[,2]==5) # diapausing - gen0<-sum(vec.mat[,1]==0) # overwintering adult - gen1<-sum(vec.mat[,1]==1) # first generation - gen2<-sum(vec.mat[,1]==2) # second generation - n.adult<-sum(vec.mat[,2]==3)+sum(vec.mat[,2]==4)+sum(vec.mat[,2]==5) # sum of all adults - gen0.pop[day]<-gen0 # gen.0 pop size - 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 - - #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="wenatchee2013sim.Rdat") - #newborn.rep<-death.rep<-adult.rep<-pop.rep<-g0.rep<-g1.rep<-g2.rep<-g0a.rep<-g1a.rep<-g2a.rep<-matrix(rep(0,365*n.rep),ncol=n.rep) - # 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? - - -######################################### -# input starting year and how many years -# n.yr and start.yr needs to be integer -# loc.name needs to be CHARACTER and matches exactly the name in the csv file!!! -data.input=function(loc, start.yr, temperature.datasets) -{ - n.yr <- length(temperature_datasets) - for (i in 1:n.yr) - { - expdata<-matrix(rep(0,365*3),nrow=365) - yr<-start.yr+i # replace 2004 with start. yr - name.input<-paste(temperature.datasets[i], ".csv", sep="") - namedat<-paste(loc, yr,".Rdat",sep="") - temp.data<-read.csv(file=name.input, header=T) - - expdata[,1]<-c(1:365) - expdata[,2]<-temp.data[c(1:365),3] #minimum - expdata[,3]<-temp.data[c(1:365),2] #maximum - save(expdata,file=namedat) - } -} -######################################### - -######################################### -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))) - } - dl # return a vector of daylength in 365 days -} -######################################### - -######################################### -hourtemp=function(L,date) -{ - 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(L) # 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 - 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 -} -######################################### \ No newline at end of file