Mercurial > repos > greg > r_hello_world

comparison bmsb.R @ 11:887b76c1b7d1 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Deleted selected files
author greg
date Tue, 16 Aug 2016 09:51:53 -0400
parents a00ffd13739c
children
comparison
equal deleted inserted replaced
10:a00ffd13739c 11:887b76c1b7d1
1 #!/usr/bin/env Rscript
2
3 #suppressPackageStartupMessages(library("optparse"))
4 #suppressPackageStartupMessages(library("rjson"))
5
6 option_list <- list(
7 make_option(c("-s", "--save_log"), action="store_true", default=FALSE, help="Save R logs"),
8 make_option(c("-m", "--output_r_log"), action="store", help="Output dataset for R logs"),
9 make_option(c("-o", "--output"), action="store", help="Output dataset")
10 )
11
12 daylength = function(L){
13 # from Forsythe 1995
14 p = 0.8333
15 dl <- NULL
16 for (i in 1:365) {
17 theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (i - 186)))
18 phi <- asin(0.39795 * cos(theta))
19 dl[i] <- 24 - 24/pi * acos((sin(p * pi/180) + sin(L * pi/180) * sin(phi))/(cos(L * pi/180) * cos(phi)))
20 }
21 # return a vector of daylength in 365 days
22 dl
23 }
24
25
26 # source("daylength.R")
27 hourtemp = function(L,date){
28 # L = 37.5 specify this in main program
29 # base development threshold for BMSB
30 threshold <- 12.7
31 # threshold2 <- threshold/24 degree hour accumulation
32 #expdata <- tempdata[1:365,11:13] # Use daily max, min, mean
33 # daily minimum
34 dnp <- expdata[date,2]
35 # daily maximum
36 dxp <- expdata[date,3]
37 dmean <- 0.5 * (dnp + dxp)
38 #if (dmean>0) {
39 #dnp <- dnp - k1 * dmean
40 #dxp <- dxp + k2 * dmean
41 #} else {
42 #dnp <- dnp + k1 * dmean
43 #dxp <- dxp - k2 * dmean
44 #}
45 dd <- 0 # initialize degree day accumulation
46
47 if (dxp<threshold) {
48 dd <- 0
49 } else {
50 # extract daylength data for entire year
51 dlprofile <- daylength(L)
52 # initialize hourly temperature
53 T <- NULL
54 #initialize degree hour vector
55 dh <- NULL
56 # calculate daylength in given date
57 # date <- 200
58 y <- dlprofile[date]
59 # night length
60 z <- 24 - y
61 # lag coefficient
62 a <- 1.86
63 # night coefficient
64 b <- 2.20
65 #import raw data set
66 #tempdata <- read.csv("tempdata.csv")
67 # Should be outside function otherwise its redundant
68 # sunrise time
69 risetime <- 12 - y/2
70 # sunset time
71 settime <- 12 + y/2
72 ts <- (dxp - dnp) * sin(pi * (settime-5)/(y + 2 * a)) + dnp
73 for (i in 1:24) {
74 if (i > risetime && i < settime) {
75 # number of hours after Tmin until sunset
76 m <- i - 5
77 T[i] = (dxp - dnp) * sin(pi * m/(y + 2 * a)) + dnp
78 if (T[i]<8.4) {
79 dh[i] <- 0
80 } else {
81 dh[i] <- T[i] - 8.4
82 }
83 } else if (i>settime) {
84 n <- i - settime
85 T[i] = dnp + (ts - dnp) * exp(-b * n/z)
86 if (T[i]<8.4) {
87 dh[i] <- 0
88 } else {
89 dh[i] <- T[i] - 8.4
90 }
91 } else {
92 n <- i + 24 - settime
93 T[i] = dnp + (ts - dnp) * exp(-b * n / z)
94 if (T[i]<8.4) {
95 dh[i] <- 0
96 } else {
97 dh[i] <- T[i] - 8.4
98 }
99 }
100 }
101 dd <- sum(dh) / 24
102 }
103 return = c(dmean, dd)
104 return
105 }
106
107
108 mortality.egg = function(temperature) {
109 if (temperature < 12.7) {
110 mort.prob = 1
111 } else {
112 # 100% mortality if <12.7
113 mort.prob = 0.8 - temperature / 40
114 if (mort.prob<0) {
115 mort.prob = 0.01
116 }
117 }
118 return = mort.prob
119 return
120 }
121
122
123 mortality.nymph = function(temperature) {
124 if (temperature<12.7) {
125 mort.prob = 0.03
126 } else {
127 # at low temperature
128 mort.prob = -temperature * 0.0008 + 0.03
129 }
130 return = mort.prob
131 return
132 }
133
134
135 mortality.adult = function(temperature) {
136 if (temperature < 12.7) {
137 mort.prob = 0.002
138 } else {
139 mort.prob = -temperature * 0.0005 + 0.02
140 }
141 return = mort.prob
142 return
143 }
144
145
146 # model initialization
147 # TODO: add tool params for the following options.
148 # start with 1000 individuals
149 n <- 1000
150 # Generation, Stage, DD, T, Diapause
151 vec.ini <- c(0,3,0,0,0)
152 # overwintering, previttelogenic, DD = 0, T = 0, no-diapause
153 vec.mat <- rep(vec.ini,n)
154 # complete matrix for the population
155 vec.mat <- t(matrix(vec.mat, nrow = 5))
156 # latitude for Asheville NC
157 L <- 35.58
158 # complete photoperiod profile in a year, requires daylength function
159 ph.p <- daylength(L)
160
161 # load temperature data@location/year
162 load("asheville2014.Rdat")
163
164 # time series of population size
165 tot.pop <- NULL
166
167 # gen.0 pop size
168 gen0.pop <- rep(0, 365)
169 gen1.pop <- rep(0, 365)
170 gen2.pop <- rep(0, 365)
171
172 # aggregate
173 S0 <- S1 <- S2 <- S3 <- S4 <- S5 <- rep(0, 365)
174 g0.adult <- g1.adult <- g2.adult <- rep(0, 365)
175
176 # birth death adults
177 N.newborn <- N.death <- N.adult <- rep(0, 365)
178
179 # degree-day accumulation
180 dd.day <- rep(0, 365)
181
182 # start tick
183 ptm <- proc.time()
184
185 for (n.sim in 1:1000) {
186 # loop through 1000 simulations
187 for (day in 1:365) {
188 # loop through 365 day/yr
189 photoperiod <- ph.p[day]
190 # photoperiod in the day
191 temp.profile <- hourtemp(L,day)
192 # temperature profile
193 mean.temp <- temp.profile[1]
194 # mean temp
195 dd.temp <- temp.profile[2]
196 # degree-day
197 dd.day[day] <- dd.temp
198 death.vec <- NULL
199 # trash bin for death
200 birth.vec <- NULL
201 # record new born
202 for (i in 1:n) {
203 # loop through all individual
204 vec.ind <- vec.mat[i,]
205 # find individual record
206 # first of all, still alive?
207 if (vec.ind[2] == 0) {
208 # egg
209 death.prob = mortality.egg(mean.temp)
210 } else if (vec.ind[2] == 1 | vec.ind[2] == 2) {
211 # nymph
212 death.prob = mortality.nymph(mean.temp)
213 } else if (vec.ind[2] == 3 | vec.ind[2] == 4 | vec.ind[2] == 5) {
214 # for adult
215 death.prob = mortality.adult(mean.temp)
216 }
217 u.d <- runif(1)
218 if (u.d<death.prob) {
219 death.vec <- c(death.vec,i)
220 } else {
221 # aggregrate index of dead bug
222 # event 1 end of diapause
223 if (vec.ind[1] == 0 && vec.ind[2] == 3) {
224 # overwintering adult (previttelogenic)
225 if (photoperiod>13.5 && vec.ind[3] > 68 && day < 180) {
226 # add 68C to become fully reproductively matured
227 # transfer to vittelogenic
228 vec.ind <- c(0,4,0,0,0)
229 vec.mat[i,] <- vec.ind
230 } else {
231 # add to DD
232 vec.ind[3] <- vec.ind[3] + dd.temp
233 vec.ind[4] <- vec.ind[4] + 1 # add 1 day in current stage
234 vec.mat[i,] <- vec.ind
235 }
236 }
237 if (vec.ind[1]!=0 && vec.ind[2] == 3) {
238 # NOT overwintering adult (previttelogenic)
239 current.gen <- vec.ind[1]
240 if (vec.ind[3]>68) {
241 # add 68C to become fully reproductively matured
242 # transfer to vittelogenic
243 vec.ind <- c(current.gen,4,0,0,0)
244 vec.mat[i,] <- vec.ind
245 } else {
246 # add to DD
247 vec.ind[3] <- vec.ind[3] + dd.temp
248 # add 1 day in current stage
249 vec.ind[4] <- vec.ind[4] + 1
250 vec.mat[i,] <- vec.ind
251 }
252 }
253 # event 2 oviposition -- where population dynamics comes from
254 # vittelogenic stage, overwintering generation
255 if (vec.ind[2] == 4 && vec.ind[1] == 0 && mean.temp>10) {
256 if (vec.ind[4] == 0) {
257 # just turned in vittelogenic stage
258 n.birth = round(runif(1,10,20))
259 } else {
260 p.birth = 1/4/75
261 # prob of birth
262 u1 <- runif(1)
263 if (u1<p.birth) {
264 n.birth = n.birth
265 }
266 }
267 # add to DD
268 vec.ind[3] <- vec.ind[3] + dd.temp
269 # add 1 day in current stage
270 vec.ind[4] <- vec.ind[4] + 1
271 vec.mat[i,] <- vec.ind
272 if (n.birth>0) {
273 # add new birth -- might be in different generations
274 # generation + 1
275 new.gen <- vec.ind[1] + 1
276 # egg profile
277 new.ind <- c(new.gen,0,0,0,0)
278 new.vec <- rep(new.ind,n.birth)
279 # update batch of egg profile
280 new.vec <- t(matrix(new.vec,nrow = 5))
281 # group with total eggs laid in that day
282 birth.vec <- rbind(birth.vec,new.vec)
283 }
284 }
285 # event 2 oviposition -- for gen 1.
286 if (vec.ind[2] == 4 && vec.ind[1] == 1 && mean.temp>12.5 && day<222) {
287 # vittelogenic stage, 1st generation
288 if (vec.ind[4] == 0) {
289 # just turned in vittelogenic stage
290 n.birth = round(runif(1,10,20))
291 } else {
292 p.birth = 1/4/75
293 u1 <- runif(1)
294 if (u1<p.birth) {
295 n.birth = n.birth
296 }
297 }
298 # add to DD
299 vec.ind[3] <- vec.ind[3] + dd.temp
300 # add 1 day in current stage
301 vec.ind[4] <- vec.ind[4] + 1
302 vec.mat[i,] <- vec.ind
303 if (n.birth>0) {
304 # add new birth -- might be in different generations
305 # generation + 1
306 new.gen <- vec.ind[1] + 1
307 # egg profile
308 new.ind <- c(new.gen,0,0,0,0)
309 new.vec <- rep(new.ind,n.birth)
310 # update batch of egg profile
311 new.vec <- t(matrix(new.vec,nrow = 5))
312 # group with total eggs laid in that day
313 birth.vec <- rbind(birth.vec,new.vec)
314 }
315 }
316 # event 3 development (with diapause determination)
317 # event 3.1 egg development to young nymph (vec.ind[2] = 0 -> egg)
318 # egg stage
319 if (vec.ind[2] == 0) {
320 # add to DD
321 vec.ind[3] <- vec.ind[3] + dd.temp
322 # from egg to young nymph
323 if (vec.ind[3] >= 53.30 && -0.9843 * dd.temp + 33.438>0) {
324 current.gen <- vec.ind[1]
325 # transfer to young nym stage
326 vec.ind <- c(current.gen,1,0,0,0)
327 } else {
328 # add 1 day in current stage
329 vec.ind[4] <- vec.ind[4] + 1
330 }
331 vec.mat[i,] <- vec.ind
332 }
333 # event 3.2 young nymph to old nymph (vec.ind[2] = 1 -> young nymph: determines diapause)
334 # young nymph stage
335 if (vec.ind[2] == 1) {
336 # add to DD
337 vec.ind[3] <- vec.ind[3] + dd.temp
338 # from young to old nymph
339 if (vec.ind[3] >= 537/2 && -0.45 * dd.temp + 18>0) {
340 current.gen <- vec.ind[1]
341 # transfer to old nym stage
342 vec.ind <- c(current.gen,2,0,0,0)
343 # prepare for diapausing
344 if (photoperiod<13.5 && day > 180) {
345 vec.ind[5] <- 1
346 }
347 } else {
348 # add 1 day in current stage
349 vec.ind[4] <- vec.ind[4] + 1
350 }
351 vec.mat[i,] <- vec.ind
352 }
353 # event 3.3 old nymph to adult: previttelogenic or diapausing?
354 # old nymph stage
355 if (vec.ind[2] == 2) {
356 # add to DD
357 vec.ind[3] <- vec.ind[3] + dd.temp
358 # from old to adult
359 if (vec.ind[3] >= 537/2 && -0.50 * dd.temp + 22>0) {
360 current.gen <- vec.ind[1]
361 # non-diapausing adult -- previttelogenic
362 if (vec.ind[5] == 0) {
363 vec.ind <- c(current.gen,3,0,0,0)
364 # diapausing
365 } else {
366 vec.ind <- c(current.gen,5,0,0,1)
367 }
368 } else {
369 # add 1 day in current stage
370 vec.ind[4] <- vec.ind[4] + 1
371 }
372 vec.mat[i,] <- vec.ind
373 }
374 # event 4 growing of diapausing adult (unimportant, but still necessary)##
375 if (vec.ind[2] == 5) {
376 vec.ind[3] <- vec.ind[3] + dd.temp
377 vec.ind[4] <- vec.ind[4] + 1
378 vec.mat[i,] <- vec.ind
379 }
380 } # else if it is still alive
381 } # end of the individual bug loop
382 # find how many died
383 n.death <- length(death.vec)
384 if (n.death>0) {
385 vec.mat <- vec.mat[-death.vec, ]
386 }
387 # remove record of dead
388 # find how many new born
389 n.newborn <- length(birth.vec[,1])
390 vec.mat <- rbind(vec.mat,birth.vec)
391 # update population size for the next day
392 n <- n-n.death + n.newborn
393
394 # aggregate results by day
395 tot.pop <- c(tot.pop,n)
396 # egg
397 s0 <- sum(vec.mat[,2] == 0)
398 # young nymph
399 s1 <- sum(vec.mat[,2] == 1)
400 # old nymph
401 s2 <- sum(vec.mat[,2] == 2)
402 # previtellogenic
403 s3 <- sum(vec.mat[,2] == 3)
404 # vitellogenic
405 s4 <- sum(vec.mat[,2] == 4)
406 # diapausing
407 s5 <- sum(vec.mat[,2] == 5)
408 # overwintering adult
409 gen0 <- sum(vec.mat[,1] == 0)
410 # first generation
411 gen1 <- sum(vec.mat[,1] == 1)
412 # second generation
413 gen2 <- sum(vec.mat[,1] == 2)
414 # sum of all adults
415 n.adult <- sum(vec.mat[,2] == 3) + sum(vec.mat[,2] == 4) + sum(vec.mat[,2] == 5)
416 # gen.0 pop size
417 gen0.pop[day] <- gen0
418 gen1.pop[day] <- gen1
419 gen2.pop[day] <- gen2
420 S0[day] <- s0
421 S1[day] <- s1
422 S2[day] <- s2
423 S3[day] <- s3
424 S4[day] <- s4
425 S5[day] <- s5
426 g0.adult[day] <- sum(vec.mat[,1] == 0)
427 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))
428 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))
429 N.newborn[day] <- n.newborn
430 N.death[day] <- n.death
431 N.adult[day] <- n.adult
432 }
433 #print(n.sim)
434 }
435
436 proc.time() - ptm
437 dd.cum <- cumsum(dd.day)
438 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)