Mercurial > repos > greg > bmsb

comparison bmsb.R @ 29:be7c61620bb1 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Sun, 21 Aug 2016 10:18:33 -0400
parents 79cef4a790cc
children
comparison
equal deleted inserted replaced
28:22043dadfd2c 29:be7c61620bb1
1 #!/usr/bin/env Rscript 1 #!/usr/bin/env Rscript
2 2
3 suppressPackageStartupMessages(library("optparse")) 3 suppressPackageStartupMessages(library("optparse"))
4 4
5 options_list <- list( 5 options_list <- list(
6 make_option(c("-n", "--input"), action="store", help="Input dataset") 6 make_option(c("-i", "--input"), action="store", help="Input dataset"),
7 make_option(c("-o", "--output"), action="store", help="Output dataset") 7 make_option(c("-o", "--output"), action="store", help="Output dataset")
8 ) 8 )
9 9
10 parser <- OptionParser(usage="%prog [options] file", options_list) 10 parser <- OptionParser(usage="%prog [options] file", options_list)
11 args <- parse_args(parser, positional_arguments=TRUE) 11 args <- parse_args(parser, positional_arguments=TRUE)
179 { 179 {
180 # event 1 end of diapause 180 # event 1 end of diapause
181 if (vec.ind[1]==0 && vec.ind[2]==3) { # overwintering adult (previttelogenic) 181 if (vec.ind[1]==0 && vec.ind[2]==3) { # overwintering adult (previttelogenic)
182 if (photoperiod>13.5 && vec.ind[3]>77 && day<180) { # add 77C to become fully reproductively matured 182 if (photoperiod>13.5 && vec.ind[3]>77 && day<180) { # add 77C to become fully reproductively matured
183 vec.ind<-c(0,4,0,0,0) # transfer to vittelogenic 183 vec.ind<-c(0,4,0,0,0) # transfer to vittelogenic
184 vec.mat[i,]<-vec.ind 184 vec.mat[i,]<-vec.ind
185 185
186 } else { 186 } else {
187 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD 187 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
188 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage 188 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
189 vec.mat[i,]<-vec.ind 189 vec.mat[i,]<-vec.ind
190 } 190 }
191 } 191 }
192 192
193 if (vec.ind[1]!=0 && vec.ind[2]==3) { # NOT overwintering adult (previttelogenic) 193 if (vec.ind[1]!=0 && vec.ind[2]==3) { # NOT overwintering adult (previttelogenic)
194 current.gen<-vec.ind[1] 194 current.gen<-vec.ind[1]
195 if (vec.ind[3]>77) { # add 77C to become fully reproductively matured 195 if (vec.ind[3]>77) { # add 77C to become fully reproductively matured
196 vec.ind<-c(current.gen,4,0,0,0) # transfer to vittelogenic 196 vec.ind<-c(current.gen,4,0,0,0) # transfer to vittelogenic
197 vec.mat[i,]<-vec.ind 197 vec.mat[i,]<-vec.ind
198 } else { 198 } else {
199 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD 199 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
200 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage 200 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
201 vec.mat[i,]<-vec.ind 201 vec.mat[i,]<-vec.ind
202 } 202 }
203 } 203 }
204 204
205 205
206 206
212 u1<-runif(1) 212 u1<-runif(1)
213 if (u1<p.birth) {n.birth=round(runif(1,2,8))} 213 if (u1<p.birth) {n.birth=round(runif(1,2,8))}
214 } 214 }
215 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD 215 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
216 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage 216 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
217 vec.mat[i,]<-vec.ind 217 vec.mat[i,]<-vec.ind
218 if (n.birth>0) { # add new birth -- might be in different generations 218 if (n.birth>0) { # add new birth -- might be in different generations
219 new.gen<-vec.ind[1]+1 # generation +1 219 new.gen<-vec.ind[1]+1 # generation +1
220 new.ind<-c(new.gen,0,0,0,0) # egg profile 220 new.ind<-c(new.gen,0,0,0,0) # egg profile
221 new.vec<-rep(new.ind,n.birth) 221 new.vec<-rep(new.ind,n.birth)
222 new.vec<-t(matrix(new.vec,nrow=5)) # update batch of egg profile 222 new.vec<-t(matrix(new.vec,nrow=5)) # update batch of egg profile
223 birth.vec<-rbind(birth.vec,new.vec) # group with total eggs laid in that day 223 birth.vec<-rbind(birth.vec,new.vec) # group with total eggs laid in that day
224 } 224 }
225 } 225 }
226 226
227 # event 2 oviposition -- for gen 1. 227 # event 2 oviposition -- for gen 1.
228 if (vec.ind[2]==4 && vec.ind[1]==1 && mean.temp>12.5 && day<222) { # vittelogenic stage, 1st generation 228 if (vec.ind[2]==4 && vec.ind[1]==1 && mean.temp>12.5 && day<222) { # vittelogenic stage, 1st generation
229 if (vec.ind[4]==0) { # just turned in vittelogenic stage 229 if (vec.ind[4]==0) { # just turned in vittelogenic stage
232 u1<-runif(1) 232 u1<-runif(1)
233 if (u1<p.birth) {n.birth=round(runif(1,2,8))} 233 if (u1<p.birth) {n.birth=round(runif(1,2,8))}
234 } 234 }
235 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD 235 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
236 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage 236 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
237 vec.mat[i,]<-vec.ind 237 vec.mat[i,]<-vec.ind
238 if (n.birth>0) { # add new birth -- might be in different generations 238 if (n.birth>0) { # add new birth -- might be in different generations
239 new.gen<-vec.ind[1]+1 # generation +1 239 new.gen<-vec.ind[1]+1 # generation +1
240 new.ind<-c(new.gen,0,0,0,0) # egg profile 240 new.ind<-c(new.gen,0,0,0,0) # egg profile
241 new.vec<-rep(new.ind,n.birth) 241 new.vec<-rep(new.ind,n.birth)
242 new.vec<-t(matrix(new.vec,nrow=5)) # update batch of egg profile 242 new.vec<-t(matrix(new.vec,nrow=5)) # update batch of egg profile
243 birth.vec<-rbind(birth.vec,new.vec) # group with total eggs laid in that day 243 birth.vec<-rbind(birth.vec,new.vec) # group with total eggs laid in that day
244 } 244 }
245 } 245 }
246 246
247 247
248 248
249 # event 3 development (with diapause determination) 249 # event 3 development (with diapause determination)
250 # event 3.1 egg development to young nymph (vec.ind[2]=0 -> egg) 250 # event 3.1 egg development to young nymph (vec.ind[2]=0 -> egg)
251 if (vec.ind[2]==0) { # egg stage 251 if (vec.ind[2]==0) { # egg stage
252 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD 252 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
253 if (vec.ind[3]>=68) { # from egg to young nymph, DD requirement met 253 if (vec.ind[3]>=68) { # from egg to young nymph, DD requirement met
254 current.gen<-vec.ind[1] 254 current.gen<-vec.ind[1]
255 vec.ind<-c(current.gen,1,0,0,0) # transfer to young nym stage 255 vec.ind<-c(current.gen,1,0,0,0) # transfer to young nym stage
256 } else { 256 } else {
257 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage 257 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
258 } 258 }
259 vec.mat[i,]<-vec.ind 259 vec.mat[i,]<-vec.ind
260 } 260 }
261 261
262 # event 3.2 young nymph to old nymph (vec.ind[2]=1 -> young nymph: determines diapause) 262 # event 3.2 young nymph to old nymph (vec.ind[2]=1 -> young nymph: determines diapause)
263 if (vec.ind[2]==1) { # young nymph stage 263 if (vec.ind[2]==1) { # young nymph stage
264 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD 264 vec.ind[3]<-vec.ind[3]+dd.temp # add to DD
265 if (vec.ind[3]>=250) { # from young to old nymph, DD requirement met 265 if (vec.ind[3]>=250) { # from young to old nymph, DD requirement met
266 current.gen<-vec.ind[1] 266 current.gen<-vec.ind[1]
267 vec.ind<-c(current.gen,2,0,0,0) # transfer to old nym stage 267 vec.ind<-c(current.gen,2,0,0,0) # transfer to old nym stage
268 if (photoperiod<13.5 && day > 180) {vec.ind[5]<-1} # prepare for diapausing 268 if (photoperiod<13.5 && day > 180) {vec.ind[5]<-1} # prepare for diapausing
269 } else { 269 } else {
270 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage 270 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
271 } 271 }
272 vec.mat[i,]<-vec.ind 272 vec.mat[i,]<-vec.ind
273 } 273 }
274 274
275 275
276 # event 3.3 old nymph to adult: previttelogenic or diapausing? 276 # event 3.3 old nymph to adult: previttelogenic or diapausing?
280 current.gen<-vec.ind[1] 280 current.gen<-vec.ind[1]
281 if (vec.ind[5]==0) { # non-diapausing adult -- previttelogenic 281 if (vec.ind[5]==0) { # non-diapausing adult -- previttelogenic
282 vec.ind<-c(current.gen,3,0,0,0) 282 vec.ind<-c(current.gen,3,0,0,0)
283 } else { # diapausing 283 } else { # diapausing
284 vec.ind<-c(current.gen,5,0,0,1) 284 vec.ind<-c(current.gen,5,0,0,1)
285 } 285 }
286 } else { 286 } else {
287 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage 287 vec.ind[4]<-vec.ind[4]+1 # add 1 day in current stage
288 } 288 }
289 vec.mat[i,]<-vec.ind 289 vec.mat[i,]<-vec.ind
290 } 290 }
339 339
340 340
341 N.newborn[day]<-n.newborn 341 N.newborn[day]<-n.newborn
342 N.death[day]<-n.death 342 N.death[day]<-n.death
343 N.adult[day]<-n.adult 343 N.adult[day]<-n.adult
344 print(c(day,n,n.adult)) 344 #print(c(day,n,n.adult))
345 } 345 }
346 346
347 proc.time() - ptm 347 proc.time() - ptm
348 dd.cum<-cumsum(dd.day) 348 dd.cum<-cumsum(dd.day)
349 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") 349 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)
350
351