Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 10:61bc6bd8807d draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Tue, 27 Feb 2018 13:31:59 -0500
parents 37f1ad91a949
children d1f2c0634354
comparison
equal deleted inserted replaced
9:d9371485aaf5 10:61bc6bd8807d
7 make_option(c("--adult_accumulation"), action="store", dest="adult_accumulation", type="integer", help="Adjustment of degree-days accumulation (old nymph->adult)"), 7 make_option(c("--adult_accumulation"), action="store", dest="adult_accumulation", type="integer", help="Adjustment of degree-days accumulation (old nymph->adult)"),
8 make_option(c("--egg_mortality"), action="store", dest="egg_mortality", type="integer", help="Adjustment rate for egg mortality"), 8 make_option(c("--egg_mortality"), action="store", dest="egg_mortality", type="integer", help="Adjustment rate for egg mortality"),
9 make_option(c("--input"), action="store", dest="input", help="Temperature data for selected location"), 9 make_option(c("--input"), action="store", dest="input", help="Temperature data for selected location"),
10 make_option(c("--insect"), action="store", dest="insect", help="Insect name"), 10 make_option(c("--insect"), action="store", dest="insect", help="Insect name"),
11 make_option(c("--insects_per_replication"), action="store", dest="insects_per_replication", type="integer", help="Number of insects with which to start each replication"), 11 make_option(c("--insects_per_replication"), action="store", dest="insects_per_replication", type="integer", help="Number of insects with which to start each replication"),
12 make_option(c("--life_stages"), action="store", dest="life_stages", help="Selected life stages for plotting"),
13 make_option(c("--life_stages_adult"), action="store", dest="life_stages_adult", default=NULL, help="Adult life stages for plotting"),
14 make_option(c("--life_stage_nymph"), action="store", dest="life_stage_nymph", default=NULL, help="Nymph life stages for plotting"),
12 make_option(c("--location"), action="store", dest="location", help="Selected location"), 15 make_option(c("--location"), action="store", dest="location", help="Selected location"),
13 make_option(c("--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"), 16 make_option(c("--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"),
14 make_option(c("--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"), 17 make_option(c("--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),
15 make_option(c("--nymph_mortality"), action="store", dest="nymph_mortality", type="integer", help="Adjustment rate for nymph mortality"), 18 make_option(c("--nymph_mortality"), action="store", dest="nymph_mortality", type="integer", help="Adjustment rate for nymph mortality"),
16 make_option(c("--old_nymph_accumulation"), action="store", dest="old_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (young nymph->old nymph)"), 19 make_option(c("--old_nymph_accumulation"), action="store", dest="old_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (young nymph->old nymph)"),
17 make_option(c("--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"), 20 make_option(c("--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"),
18 make_option(c("--output"), action="store", dest="output", help="Output dataset"),
19 make_option(c("--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"), 21 make_option(c("--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"),
20 make_option(c("--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"), 22 make_option(c("--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"),
21 make_option(c("--replications"), action="store", dest="replications", type="integer", help="Number of replications"), 23 make_option(c("--replications"), action="store", dest="replications", type="integer", help="Number of replications"),
22 make_option(c("--std_error_plot"), action="store", dest="std_error_plot", help="Plot Standard error"), 24 make_option(c("--plot_generations_separately"), action="store", dest="plot_generations_separately", help="Plot Plot P, F1 and F2 as separate lines or pool across them"),
25 make_option(c("--plot_std_error"), action="store", dest="plot_std_error", help="Plot Standard error"),
23 make_option(c("--young_nymph_accumulation"), action="store", dest="young_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)") 26 make_option(c("--young_nymph_accumulation"), action="store", dest="young_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)")
24 ) 27 )
25 28
26 parser <- OptionParser(usage="%prog [options] file", option_list=option_list); 29 parser <- OptionParser(usage="%prog [options] file", option_list=option_list);
27 args <- parse_args(parser, positional_arguments=TRUE); 30 args <- parse_args(parser, positional_arguments=TRUE);
70 n12 = -0.3728 * temperature + 14.68; 73 n12 = -0.3728 * temperature + 14.68;
71 n23 = -0.6119 * temperature + 25.249; 74 n23 = -0.6119 * temperature + 25.249;
72 dev.rate = mean(n12 + n23); 75 dev.rate = mean(n12 + n23);
73 return(dev.rate); 76 return(dev.rate);
74 } 77 }
75
76 78
77 get_date_labels = function(temperature_data_frame, num_rows) { 79 get_date_labels = function(temperature_data_frame, num_rows) {
78 # Keep track of the years to see if spanning years. 80 # Keep track of the years to see if spanning years.
79 month_labels = list(); 81 month_labels = list();
80 current_month_label = NULL; 82 current_month_label = NULL;
213 } 215 }
214 return(temperature_data_frame); 216 return(temperature_data_frame);
215 } 217 }
216 218
217 219
218 render_chart = function(chart_type, insect, location, latitude, start_date, end_date, days, maxval, plot_std_error, 220 render_chart = function(date_labels, chart_type, plot_std_error, insect, location, latitude, start_date, end_date, days, maxval,
219 group1, group2, group3, group1_std_error, group2_std_error, group3_std_error, date_labels) { 221 replications, life_stage, group, group_std_error, group2=NULL, group2_std_error=NULL, group3=NULL, group3_std_error=NULL,
220 if (chart_type == "pop_size_by_life_stage") { 222 life_stages_adult=NULL, life_stage_nymph=NULL) {
221 title = paste(insect, ": Total pop. by life stage :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" "); 223 if (chart_type=="pop_size_by_life_stage") {
222 legend_text = c("Egg", "Nymph", "Adult"); 224 if (life_stage=="Total") {
223 columns = c(4, 2, 1); 225 title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
224 } else if (chart_type == "pop_size_by_generation") { 226 legend_text = c("Egg", "Nymph", "Adult");
225 title = paste(insect, ": Total pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" "); 227 columns = c(4, 2, 1);
228 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3);
229 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3);
230 lines(days, group2, lwd=2, lty=1, col=2);
231 lines(days, group3, lwd=2, lty=1, col=4);
232 axis(1, at=c(1:length(date_labels)) * 30 - 15, cex.axis=3, labels=date_labels);
233 axis(2, cex.axis=3);
234 if (plot_std_error=="yes") {
235 # Standard error for group.
236 lines(days, group+group_std_error, lty=2);
237 lines(days, group-group_std_error, lty=2);
238 # Standard error for group2.
239 lines(days, group2+group2_std_error, col=2, lty=2);
240 lines(days, group2-group2_std_error, col=2, lty=2);
241 # Standard error for group3.
242 lines(days, group3+group3_std_error, col=4, lty=2);
243 lines(days, group3-group3_std_error, col=4, lty=2);
244 }
245 } else {
246 if (life_stage=="Egg") {
247 title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
248 legend_text = c(life_stage);
249 columns = c(4);
250 } else if (life_stage=="Nymph") {
251 stage = paste(life_stage_nymph, "Nymph Pop :", sep=" ");
252 title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
253 legend_text = c(paste(life_stage_nymph, life_stage, sep=" "));
254 columns = c(2);
255 } else if (life_stage=="Adult") {
256 stage = paste(life_stages_adult, "Adult Pop", sep=" ");
257 title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
258 legend_text = c(paste(life_stages_adult, life_stage, sep=" "));
259 columns = c(1);
260 }
261 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3);
262 legend("topleft", legend_text, lty=c(1), col="black", cex=3);
263 axis(1, at=c(1:length(date_labels)) * 30 - 15, cex.axis=3, labels=date_labels);
264 axis(2, cex.axis=3);
265 if (plot_std_error=="yes") {
266 # Standard error for group.
267 lines(days, group+group_std_error, lty=2);
268 lines(days, group-group_std_error, lty=2);
269 }
270 }
271 } else if (chart_type=="pop_size_by_generation") {
272 if (life_stage=="Total") {
273 title_str = ": Total Pop by Gen :";
274 } else if (life_stage=="Egg") {
275 title_str = ": Egg Pop by Gen :";
276 } else if (life_stage=="Nymph") {
277 title_str = paste(":", life_stage_nymph, "Nymph Pop by Gen", ":", sep=" ");
278 } else if (life_stage=="Adult") {
279 title_str = paste(":", life_stages_adult, "Adult Pop by Gen", ":", sep=" ");
280 }
281 title = paste(insect, ": Reps", replications, title_str, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
226 legend_text = c("P", "F1", "F2"); 282 legend_text = c("P", "F1", "F2");
227 columns = c(1, 2, 4); 283 columns = c(1, 2, 4);
228 } else if (chart_type == "adult_pop_size_by_generation") { 284 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3);
229 title = paste(insect, ": Adult pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" "); 285 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3);
230 legend_text = c("P", "F1", "F2"); 286 lines(days, group2, lwd=2, lty=1, col=2);
231 columns = c(1, 2, 4); 287 lines(days, group3, lwd=2, lty=1, col=4);
232 } 288 axis(1, at=c(1:length(date_labels)) * 30 - 15, cex.axis=3, labels=date_labels);
233 plot(days, group1, main=title, type="l", ylim=c(0, maxval), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3); 289 axis(2, cex.axis=3);
234 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3); 290 if (plot_std_error=="yes") {
235 lines(days, group2, lwd=2, lty=1, col=2); 291 # Standard error for group.
236 lines(days, group3, lwd=2, lty=1, col=4); 292 lines(days, group+group_std_error, lty=2);
237 axis(1, at=c(1:length(date_labels)) * 30 - 15, cex.axis=3, labels=date_labels); 293 lines(days, group-group_std_error, lty=2);
238 axis(2, cex.axis=3); 294 # Standard error for group2.
239 if (plot_std_error==1) { 295 lines(days, group2+group2_std_error, col=2, lty=2);
240 # Standard error for group1. 296 lines(days, group2-group2_std_error, col=2, lty=2);
241 lines(days, group1+group1_std_error, lty=2); 297 # Standard error for group3.
242 lines(days, group1-group1_std_error, lty=2); 298 lines(days, group3+group3_std_error, col=4, lty=2);
243 # Standard error for group2. 299 lines(days, group3-group3_std_error, col=4, lty=2);
244 lines(days, group2+group2_std_error, col=2, lty=2); 300 }
245 lines(days, group2-group2_std_error, col=2, lty=2); 301 }
246 # Standard error for group3. 302 }
247 lines(days, group3+group3_std_error, col=4, lty=2); 303
248 lines(days, group3-group3_std_error, col=4, lty=2); 304 # Determine if we're plotting generations separately.
249 } 305 if (opt$plot_generations_separately=="yes") {
250 } 306 plot_generations_separately = TRUE;
251 307 } else {
308 plot_generations_separately = FALSE;
309 }
310 # Read the temperature data into a data frame.
252 temperature_data_frame = parse_input_data(opt$input, opt$num_days); 311 temperature_data_frame = parse_input_data(opt$input, opt$num_days);
253 # All latitude values are the same, so get the value from the first row. 312 output_dir = "output_dir";
313 # Get the date labels for plots.
314 date_labels = get_date_labels(temperature_data_frame, opt$num_days);
315 # All latitude values are the same, so get the value for plots from the first row.
254 latitude = temperature_data_frame$LATITUDE[1]; 316 latitude = temperature_data_frame$LATITUDE[1];
317 # Get the number of days for plots.
255 num_columns = dim(temperature_data_frame)[2]; 318 num_columns = dim(temperature_data_frame)[2];
256 date_labels = get_date_labels(temperature_data_frame, opt$num_days); 319 # Split life_stages into a list of strings for plots.
320 life_stages_str = as.character(opt$life_stages);
321 life_stages = strsplit(life_stages_str, ",")[[1]];
322 # Determine the data we need to generate for plotting.
323 process_eggs = FALSE;
324 process_nymphs = FALSE;
325 process_adults = FALSE;
326 for (life_stage in life_stages) {
327 if (life_stage=="Total") {
328 process_eggs = TRUE;
329 process_nymphs = TRUE;
330 life_stage_nymph = "Total";
331 process_adults = TRUE;
332 life_stages_adult = "Total";
333 } else if (life_stage=="Egg") {
334 process_eggs = TRUE;
335 } else if (life_stage=="Nymph") {
336 process_nymphs = TRUE;
337 life_stage_nymph = opt$life_stage_nymph;
338 } else if (life_stage=="Adult") {
339 process_adults = TRUE;
340 life_stages_adult = opt$life_stages_adult;
341 }
342 }
257 343
258 # Initialize matrices. 344 # Initialize matrices.
259 Eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 345 if (process_eggs) {
260 YoungNymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 346 Eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
261 OldNymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 347 }
262 Previtellogenic.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 348 if (process_nymphs) {
263 Vitellogenic.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 349 YoungNymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
264 Diapausing.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 350 OldNymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
265 351 }
352 if (process_adults) {
353 Previtellogenic.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
354 Vitellogenic.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
355 Diapausing.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
356 }
266 newborn.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 357 newborn.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
267 adult.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 358 adult.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
268 death.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 359 death.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
269 360 if (plot_generations_separately) {
270 P.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 361 # P is Parental, or overwintered adults.
271 P_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 362 P.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
272 F1.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 363 # F1 is the first field-produced generation.
273 F1_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 364 F1.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
274 F2.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 365 # F2 is the second field-produced generation.
275 F2_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 366 F2.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
276 367 if (process_eggs) {
368 P_eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
369 F1_eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
370 F2_eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
371 }
372 if (process_nymphs) {
373 P_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
374 F1_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
375 F2_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
376 }
377 if (process_adults) {
378 P_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
379 F1_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
380 F2_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
381 }
382 }
383 # Total population.
277 population.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 384 population.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
278 385
279 # Process replications. 386 # Process replications.
280 for (N.replications in 1:opt$replications) { 387 for (N.replications in 1:opt$replications) {
281 # Start with the user-defined number of insects per replication. 388 # Start with the user-defined number of insects per replication.
282 num_insects = opt$insects_per_replication; 389 num_insects = opt$insects_per_replication;
283 # Generation, Stage, degree-days, T, Diapause. 390 # Generation, Stage, degree-days, T, Diapause.
284 vector.ini = c(0, 3, 0, 0, 0); 391 vector.ini = c(0, 3, 0, 0, 0);
285 # Overwintering, previttelogenic, degree-days=0, T=0, no-diapause. 392 # Replicate to create a matrix where the columns are
393 # Generation, Stage, degree-days, T, Diapause and the
394 # rows are the initial number of insects per replication.
286 vector.matrix = rep(vector.ini, num_insects); 395 vector.matrix = rep(vector.ini, num_insects);
287 # Complete matrix for the population. 396 # Complete transposed matrix for the population, so now
397 # the rows are Generation, Stage, degree-days, T, Diapause
288 vector.matrix = base::t(matrix(vector.matrix, nrow=5)); 398 vector.matrix = base::t(matrix(vector.matrix, nrow=5));
289 # Time series of population size. 399 # Time series of population size.
290 Eggs = rep(0, opt$num_days); 400 if (process_eggs) {
291 YoungNymphs = rep(0, opt$num_days); 401 Eggs = rep(0, opt$num_days);
292 OldNymphs = rep(0, opt$num_days); 402 }
293 Previtellogenic = rep(0, opt$num_days); 403 if (process_nymphs) {
294 Vitellogenic = rep(0, opt$num_days); 404 YoungNymphs = rep(0, opt$num_days);
295 Diapausing = rep(0, opt$num_days); 405 OldNymphs = rep(0, opt$num_days);
296 406 }
407 if (process_adults) {
408 Previtellogenic = rep(0, opt$num_days);
409 Vitellogenic = rep(0, opt$num_days);
410 Diapausing = rep(0, opt$num_days);
411 }
297 N.newborn = rep(0, opt$num_days); 412 N.newborn = rep(0, opt$num_days);
298 N.adult = rep(0, opt$num_days); 413 N.adult = rep(0, opt$num_days);
299 N.death = rep(0, opt$num_days); 414 N.death = rep(0, opt$num_days);
300
301 overwintering_adult.population = rep(0, opt$num_days); 415 overwintering_adult.population = rep(0, opt$num_days);
302 first_generation.population = rep(0, opt$num_days); 416 first_generation.population = rep(0, opt$num_days);
303 second_generation.population = rep(0, opt$num_days); 417 second_generation.population = rep(0, opt$num_days);
304 418 if (plot_generations_separately) {
305 P.adult = rep(0, opt$num_days); 419 # P is Parental, or overwintered adults.
306 F1.adult = rep(0, opt$num_days); 420 # F1 is the first field-produced generation.
307 F2.adult = rep(0, opt$num_days); 421 # F2 is the second field-produced generation.
308 422 if (process_eggs) {
423 P.egg = rep(0, opt$num_days);
424 F1.egg = rep(0, opt$num_days);
425 F2.egg = rep(0, opt$num_days);
426 }
427 if (process_nymphs) {
428 P.nymph = rep(0, opt$num_days);
429 F1.nymph = rep(0, opt$num_days);
430 F2.nymph = rep(0, opt$num_days);
431 }
432 if (process_adults) {
433 P.adult = rep(0, opt$num_days);
434 F1.adult = rep(0, opt$num_days);
435 F2.adult = rep(0, opt$num_days);
436 }
437 }
309 total.population = NULL; 438 total.population = NULL;
310
311 averages.day = rep(0, opt$num_days); 439 averages.day = rep(0, opt$num_days);
312 # All the days included in the input temperature dataset. 440 # All the days included in the input temperature dataset.
313 for (row in 1:opt$num_days) { 441 for (row in 1:opt$num_days) {
314 # Get the integer day of the year for the current row. 442 # Get the integer day of the year for the current row.
315 doy = temperature_data_frame$DOY[row]; 443 doy = temperature_data_frame$DOY[row];
533 num_insects.newborn = length(birth.vector[,1]); 661 num_insects.newborn = length(birth.vector[,1]);
534 vector.matrix = rbind(vector.matrix, birth.vector); 662 vector.matrix = rbind(vector.matrix, birth.vector);
535 # Update population size for the next day. 663 # Update population size for the next day.
536 num_insects = num_insects - num_insects.death + num_insects.newborn; 664 num_insects = num_insects - num_insects.death + num_insects.newborn;
537 665
538 # Aggregate results by day. 666 # Aggregate results by day. Due to multiple transpose calls
539 # Egg population size. 667 # on vector.matrix above, the columns of vector.matrix
540 Eggs[row] = sum(vector.matrix[,2]==0); 668 # are now Generation, Stage, degree-days, T, Diapause,
541 # Young nymph population size. 669 if (process_eggs) {
542 YoungNymphs[row] = sum(vector.matrix[,2]==1); 670 # For egg population size, column 2 (Stage), must be 0.
543 # Old nymph population size. 671 Eggs[row] = sum(vector.matrix[,2]==0);
544 OldNymphs[row] = sum(vector.matrix[,2]==2); 672 }
545 # Previtellogenic population size. 673 if (process_nymphs) {
546 Previtellogenic[row] = sum(vector.matrix[,2]==3); 674 # For young nymph population size, column 2 (Stage) must be 1.
547 # Vitellogenic population size. 675 YoungNymphs[row] = sum(vector.matrix[,2]==1);
548 Vitellogenic[row] = sum(vector.matrix[,2]==4); 676 # For old nymph population size, column 2 (Stage) must be 2.
549 # Diapausing population size. 677 OldNymphs[row] = sum(vector.matrix[,2]==2);
550 Diapausing[row] = sum(vector.matrix[,2]==5); 678 }
679 if (process_adults) {
680 # For pre-vitellogenic population size, column 2 (Stage) must be 3.
681 Previtellogenic[row] = sum(vector.matrix[,2]==3);
682 # For vitellogenic population size, column 2 (Stage) must be 4.
683 Vitellogenic[row] = sum(vector.matrix[,2]==4);
684 # For diapausing population size, column 2 (Stage) must be 5.
685 Diapausing[row] = sum(vector.matrix[,2]==5);
686 }
551 687
552 # Newborn population size. 688 # Newborn population size.
553 N.newborn[row] = num_insects.newborn; 689 N.newborn[row] = num_insects.newborn;
554 # Adult population size. 690 # Adult population size.
555 N.adult[row] = sum(vector.matrix[,2]==3) + sum(vector.matrix[,2]==4) + sum(vector.matrix[,2]==5); 691 N.adult[row] = sum(vector.matrix[,2]==3) + sum(vector.matrix[,2]==4) + sum(vector.matrix[,2]==5);
556 # Dead population size. 692 # Dead population size.
557 N.death[row] = num_insects.death; 693 N.death[row] = num_insects.death;
558 694
559 total.population = c(total.population, num_insects); 695 total.population = c(total.population, num_insects);
560 696
561 # Overwintering adult population size. 697 # For overwintering adult (P) population
698 # size, column 1 (Generation) must be 0.
562 overwintering_adult.population[row] = sum(vector.matrix[,1]==0); 699 overwintering_adult.population[row] = sum(vector.matrix[,1]==0);
563 # First generation population size. 700 # For first field generation (F1) population
701 # size, column 1 (Generation) must be 1.
564 first_generation.population[row] = sum(vector.matrix[,1]==1); 702 first_generation.population[row] = sum(vector.matrix[,1]==1);
565 # Second generation population size. 703 # For second field generation (F2) population
704 # size, column 1 (Generation) must be 2.
566 second_generation.population[row] = sum(vector.matrix[,1]==2); 705 second_generation.population[row] = sum(vector.matrix[,1]==2);
567 706
568 # P adult population size. 707 if (plot_generations_separately) {
569 P.adult[row] = sum(vector.matrix[,1]==0); 708 if (process_eggs) {
570 # F1 adult population size. 709 # For egg life stage of generation F1 population size,
571 F1.adult[row] = sum((vector.matrix[,1]==1 & vector.matrix[,2]==3) | (vector.matrix[,1]==1 & vector.matrix[,2]==4) | (vector.matrix[,1]==1 & vector.matrix[,2]==5)); 710 # column 1 (generation) is 0 and column 2 (Stage) is 0.
572 # F2 adult population size 711 P.egg[row] = sum(vector.matrix[,1]==0 & vector.matrix[,2]==0);
573 F2.adult[row] = sum((vector.matrix[,1]==2 & vector.matrix[,2]==3) | (vector.matrix[,1]==2 & vector.matrix[,2]==4) | (vector.matrix[,1]==2 & vector.matrix[,2]==5)); 712 # For egg life stage of generation F1 population size,
713 # column 1 (generation) is 1 and column 2 (Stage) is 0.
714 F1.egg[row] = sum(vector.matrix[,1]==1 & vector.matrix[,2]==0);
715 # For egg life stage of generation F2 population size,
716 # column 1 (generation) is 2 and column 2 (Stage) is 0.
717 F2.egg[row] = sum(vector.matrix[,1]==2 & vector.matrix[,2]==0);
718 }
719 if (process_nymphs) {
720 # For nymph life stage of generation F1 population
721 # size, one of the following combinations is required:
722 # - column 1 (Generation) is 0 and column 2 (Stage) is 1 (Young nymph)
723 # - column 1 (Generation) is 0 and column 2 (Stage) is 2 (Old nymph)
724 P.nymph[row] = sum((vector.matrix[,1]==0 & vector.matrix[,2]==1) | (vector.matrix[,1]==0 & vector.matrix[,2]==2));
725 # For nymph life stage of generation F1 population
726 # size, one of the following combinations is required:
727 # - column 1 (Generation) is 1 and column 2 (Stage) is 1 (Young nymph)
728 # - column 1 (Generation) is 1 and column 2 (Stage) is 2 (Old nymph)
729 F1.nymph[row] = sum((vector.matrix[,1]==1 & vector.matrix[,2]==1) | (vector.matrix[,1]==1 & vector.matrix[,2]==2));
730 # For nymph life stage of generation F2 population
731 # size, one of the following combinations is required:
732 # - column 1 (Generation) is 2 and column 2 (Stage) is 1 (Young nymph)
733 # - column 1 (Generation) is 2 and column 2 (Stage) is 2 (Old nymph)
734 F2.nymph[row] = sum((vector.matrix[,1]==2 & vector.matrix[,2]==1) | (vector.matrix[,1]==2 & vector.matrix[,2]==2));
735 }
736 if (process_adults) {
737 # For adult life stage of generation P population
738 # size, one of the following combinations is required:
739 # - column 1 (Generation) is 0 and column 2 (Stage) is 3 (Pre-vitellogenic)
740 # - column 1 (Generation) is 0 and column 2 (Stage) is 4 (Vitellogenic)
741 # - column 1 (Generation) is 0 and column 2 (Stage) is 5 (Diapausing)
742 P.adult[row] = sum((vector.matrix[,1]==0 & vector.matrix[,2]==3) | (vector.matrix[,1]==0 & vector.matrix[,2]==4) | (vector.matrix[,1]==0 & vector.matrix[,2]==5));
743 # For adult life stage of generation F1 population
744 # size, one of the following combinations is required:
745 # - column 1 (Generation) is 1 and column 2 (Stage) is 3 (Pre-vitellogenic)
746 # - column 1 (Generation) is 1 and column 2 (Stage) is 4 (Vitellogenic)
747 # - column 1 (Generation) is 1 and column 2 (Stage) is 5 (Diapausing)
748 F1.adult[row] = sum((vector.matrix[,1]==1 & vector.matrix[,2]==3) | (vector.matrix[,1]==1 & vector.matrix[,2]==4) | (vector.matrix[,1]==1 & vector.matrix[,2]==5));
749 # For adult life stage of generation F2 population
750 # size, one of the following combinations is required:
751 # - column 1 (Generation) is 2 and column 2 (Stage) is 3 (Pre-vitellogenic)
752 # - column 1 (Generation) is 2 and column 2 (Stage) is 4 (Vitellogenic)
753 # - column 1 (Generation) is 2 and column 2 (Stage) is 5 (Diapausing)
754 F2.adult[row] = sum((vector.matrix[,1]==2 & vector.matrix[,2]==3) | (vector.matrix[,1]==2 & vector.matrix[,2]==4) | (vector.matrix[,1]==2 & vector.matrix[,2]==5));
755 }
756 }
574 } # End of days specified in the input temperature data. 757 } # End of days specified in the input temperature data.
575 758
576 averages.cum = cumsum(averages.day); 759 averages.cum = cumsum(averages.day);
577 760
578 # Define the output values. 761 # Define the output values.
579 Eggs.replications[,N.replications] = Eggs; 762 if (process_eggs) {
580 YoungNymphs.replications[,N.replications] = YoungNymphs; 763 Eggs.replications[,N.replications] = Eggs;
581 OldNymphs.replications[,N.replications] = OldNymphs; 764 }
582 Previtellogenic.replications[,N.replications] = Previtellogenic; 765 if (process_nymphs) {
583 Vitellogenic.replications[,N.replications] = Vitellogenic; 766 YoungNymphs.replications[,N.replications] = YoungNymphs;
584 Diapausing.replications[,N.replications] = Diapausing; 767 OldNymphs.replications[,N.replications] = OldNymphs;
585 768 }
769 if (process_adults) {
770 Previtellogenic.replications[,N.replications] = Previtellogenic;
771 Vitellogenic.replications[,N.replications] = Vitellogenic;
772 Diapausing.replications[,N.replications] = Diapausing;
773 }
586 newborn.replications[,N.replications] = N.newborn; 774 newborn.replications[,N.replications] = N.newborn;
587 adult.replications[,N.replications] = N.adult; 775 adult.replications[,N.replications] = N.adult;
588 death.replications[,N.replications] = N.death; 776 death.replications[,N.replications] = N.death;
589 777 if (plot_generations_separately) {
590 P.replications[,N.replications] = overwintering_adult.population; 778 # P is Parental, or overwintered adults.
591 P_adults.replications[,N.replications] = P.adult; 779 P.replications[,N.replications] = overwintering_adult.population;
592 F1.replications[,N.replications] = first_generation.population; 780 # F1 is the first field-produced generation.
593 F1_adults.replications[,N.replications] = F1.adult; 781 F1.replications[,N.replications] = first_generation.population;
594 F2.replications[,N.replications] = second_generation.population; 782 # F2 is the second field-produced generation.
595 F2_adults.replications[,N.replications] = F2.adult; 783 F2.replications[,N.replications] = second_generation.population;
596 784 if (process_eggs) {
785 P_eggs.replications[,N.replications] = P.egg;
786 F1_eggs.replications[,N.replications] = F1.egg;
787 F2_eggs.replications[,N.replications] = F2.egg;
788 }
789 if (process_nymphs) {
790 P_nymphs.replications[,N.replications] = P.nymph;
791 F1_nymphs.replications[,N.replications] = F1.nymph;
792 F2_nymphs.replications[,N.replications] = F2.nymph;
793 }
794 if (process_adults) {
795 P_adults.replications[,N.replications] = P.adult;
796 F1_adults.replications[,N.replications] = F1.adult;
797 F2_adults.replications[,N.replications] = F2.adult;
798 }
799 }
597 population.replications[,N.replications] = total.population; 800 population.replications[,N.replications] = total.population;
598 } 801 }
599 802
600 # Mean value for eggs. 803 if (process_eggs) {
601 eggs = apply(Eggs.replications, 1, mean); 804 # Mean value for eggs.
602 # Standard error for eggs. 805 eggs = apply(Eggs.replications, 1, mean);
603 eggs.std_error = apply(Eggs.replications, 1, sd) / sqrt(opt$replications); 806 # Standard error for eggs.
604 807 eggs.std_error = apply(Eggs.replications, 1, sd) / sqrt(opt$replications);
605 # Mean value for nymphs. 808 }
606 nymphs = apply((YoungNymphs.replications+OldNymphs.replications), 1, mean); 809 if (process_nymphs) {
607 # Standard error for nymphs. 810 # Calculate nymph populations for selected life stage.
608 nymphs.std_error = apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd); 811 if (life_stage_nymph=="Total") {
609 812 # Mean value for all nymphs.
610 # Mean value for adults. 813 nymphs = apply((YoungNymphs.replications+OldNymphs.replications), 1, mean);
611 adults = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, mean); 814 # Standard error for all nymphs.
612 # Standard error for adults. 815 nymphs.std_error = apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd);
613 adults.std_error = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications); 816 } else if (life_stage_nymph=="Young") {
614 817 # Mean value for young nymphs.
615 # Mean value for P. 818 nymphs = apply(YoungNymphs.replications, 1, mean);
616 P = apply(P.replications, 1, mean); 819 # Standard error for young nymphs.
617 # Standard error for P. 820 nymphs.std_error = apply(YoungNymphs.replications / sqrt(opt$replications), 1, sd);
618 P.std_error = apply(P.replications, 1, sd) / sqrt(opt$replications); 821 } else if (life_stage_nymph=="Old") {
619 822 # Mean value for old nymphs.
620 # Mean value for P adults. 823 nymphs = apply(OldNymphs.replications, 1, mean);
621 P_adults = apply(P_adults.replications, 1, mean); 824 # Standard error for old nymphs.
622 # Standard error for P_adult. 825 nymphs.std_error = apply(OldNymphs.replications / sqrt(opt$replications), 1, sd);
623 P_adults.std_error = apply(P_adults.replications, 1, sd) / sqrt(opt$replications); 826 }
624 827 }
625 # Mean value for F1. 828 if (process_adults) {
626 F1 = apply(F1.replications, 1, mean); 829 # Calculate adult populations for selected life stage.
627 # Standard error for F1. 830 if (life_stages_adult=="Total") {
628 F1.std_error = apply(F1.replications, 1, sd) / sqrt(opt$replications); 831 # Mean value for all adults.
629 832 adults = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, mean);
630 # Mean value for F1 adults. 833 # Standard error for all adults.
631 F1_adults = apply(F1_adults.replications, 1, mean); 834 adults.std_error = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications);
632 # Standard error for F1 adult. 835 } else if (life_stages_adult == "Pre-vittelogenic") {
633 F1_adults.std_error = apply(F1_adults.replications, 1, sd) / sqrt(opt$replications); 836 # Mean value for previtellogenic adults.
634 837 adults = apply(Previtellogenic.replications, 1, mean);
635 # Mean value for F2. 838 # Standard error for previtellogenic adults.
636 F2 = apply(F2.replications, 1, mean); 839 adults.std_error = apply(Previtellogenic.replications, 1, sd) / sqrt(opt$replications);
637 # Standard error for F2. 840 } else if (life_stages_adult == "Vittelogenic") {
638 F2.std_error = apply(F2.replications, 1, sd) / sqrt(opt$replications); 841 # Mean value for vitellogenic adults.
639 842 adults = apply(Vitellogenic.replications, 1, mean);
640 # Mean value for F2 adults. 843 # Standard error for vitellogenic adults.
641 F2_adults = apply(F2_adults.replications, 1, mean); 844 adults.std_error = apply(Vitellogenic.replications, 1, sd) / sqrt(opt$replications);
642 # Standard error for F2 adult. 845 } else if (life_stages_adult == "Diapausing") {
643 F2_adults.std_error = apply(F2_adults.replications, 1, sd) / sqrt(opt$replications); 846 # Mean value for vitellogenic adults.
847 adults = apply(Diapausing.replications, 1, mean);
848 # Standard error for vitellogenic adults.
849 adults.std_error = apply(Diapausing.replications, 1, sd) / sqrt(opt$replications);
850 }
851 }
852
853 if (plot_generations_separately) {
854 # Mean value for P which is Parental, or overwintered adults.
855 P = apply(P.replications, 1, mean);
856 # Standard error for P.
857 P.std_error = apply(P.replications, 1, sd) / sqrt(opt$replications);
858 # Mean value for F1, which is the first field-produced generation.
859 F1 = apply(F1.replications, 1, mean);
860 # Standard error for F1.
861 F1.std_error = apply(F1.replications, 1, sd) / sqrt(opt$replications);
862 # Mean value for F2, which is the second field-produced generation.
863 F2 = apply(F2.replications, 1, mean);
864 # Standard error for F2.
865 F2.std_error = apply(F2.replications, 1, sd) / sqrt(opt$replications);
866 if (process_eggs) {
867 # Mean value for P eggs.
868 P_eggs = apply(P_eggs.replications, 1, mean);
869 # Standard error for P_eggs.
870 P_eggs.std_error = apply(P_eggs.replications, 1, sd) / sqrt(opt$replications);
871 # Mean value for F1 eggs.
872 F1_eggs = apply(F1_eggs.replications, 1, mean);
873 # Standard error for F1 eggs.
874 F1_eggs.std_error = apply(F1_eggs.replications, 1, sd) / sqrt(opt$replications);
875 # Mean value for F2 eggs.
876 F2_eggs = apply(F2_eggs.replications, 1, mean);
877 # Standard error for F2 eggs.
878 F2_eggs.std_error = apply(F2_eggs.replications, 1, sd) / sqrt(opt$replications);
879 }
880 if (process_nymphs) {
881 # Mean value for P nymphs.
882 P_nymphs = apply(P_nymphs.replications, 1, mean);
883 # Standard error for P_nymphs.
884 P_nymphs.std_error = apply(P_nymphs.replications, 1, sd) / sqrt(opt$replications);
885 # Mean value for F1 nymphs.
886 F1_nymphs = apply(F1_nymphs.replications, 1, mean);
887 # Standard error for F1 nymphs.
888 F1_nymphs.std_error = apply(F1_nymphs.replications, 1, sd) / sqrt(opt$replications);
889 # Mean value for F2 nymphs.
890 F2_nymphs = apply(F2_nymphs.replications, 1, mean);
891 # Standard error for F2 eggs.
892 F2_nymphs.std_error = apply(F2_nymphs.replications, 1, sd) / sqrt(opt$replications);
893 }
894 if (process_adults) {
895 # Mean value for P adults.
896 P_adults = apply(P_adults.replications, 1, mean);
897 # Standard error for P_adults.
898 P_adults.std_error = apply(P_adults.replications, 1, sd) / sqrt(opt$replications);
899 # Mean value for F1 adults.
900 F1_adults = apply(F1_adults.replications, 1, mean);
901 # Standard error for F1 adults.
902 F1_adults.std_error = apply(F1_adults.replications, 1, sd) / sqrt(opt$replications);
903 # Mean value for F2 adults.
904 F2_adults = apply(F2_adults.replications, 1, mean);
905 # Standard error for F2 adults.
906 F2_adults.std_error = apply(F2_adults.replications, 1, sd) / sqrt(opt$replications);
907 }
908 }
644 909
645 # Display the total number of days in the Galaxy history item blurb. 910 # Display the total number of days in the Galaxy history item blurb.
646 cat("Number of days: ", opt$num_days, "\n"); 911 cat("Number of days: ", opt$num_days, "\n");
647 912
648 dev.new(width=20, height=30); 913 # Information needed for plots plots.
649
650 # Start PDF device driver to save charts to output.
651 pdf(file=opt$output, width=20, height=30, bg="white");
652 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
653
654 # Data analysis and visualization plots only within a single calendar year.
655 days = c(1:opt$num_days); 914 days = c(1:opt$num_days);
656 start_date = temperature_data_frame$DATE[1]; 915 start_date = temperature_data_frame$DATE[1];
657 end_date = temperature_data_frame$DATE[opt$num_days]; 916 end_date = temperature_data_frame$DATE[opt$num_days];
658 917
659 # Subfigure 1: population size by life stage. 918 if (plot_generations_separately) {
660 maxval = max(eggs+eggs.std_error, nymphs+nymphs.std_error, adults+adults.std_error); 919 for (life_stage in life_stages) {}
661 render_chart("pop_size_by_life_stage", opt$insect, opt$location, latitude, start_date, end_date, days, maxval, 920 if (life_stage == "Egg") {
662 opt$std_error_plot, adults, nymphs, eggs, adults.std_error, nymphs.std_error, eggs.std_error, date_labels); 921 # Start PDF device driver.
663 # Subfigure 2: population size by generation. 922 dev.new(width=20, height=30);
664 maxval = max(F2); 923 file_path = paste(output_dir, "egg_pop_by_generation.pdf", sep="/");
665 render_chart("pop_size_by_generation", opt$insect, opt$location, latitude, start_date, end_date, days, maxval, 924 pdf(file=file_path, width=20, height=30, bg="white");
666 opt$std_error_plot, P, F1, F2, P.std_error, F1.std_error, F2.std_error, date_labels); 925 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
667 # Subfigure 3: adult population size by generation. 926 # Egg population size by generation.
668 maxval = max(F2_adults) + 100; 927 maxval = max(F2_eggs) + 100;
669 render_chart("adult_pop_size_by_generation", opt$insect, opt$location, latitude, start_date, end_date, days, maxval, 928 render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
670 opt$std_error_plot, P_adults, F1_adults, F2_adults, P_adults.std_error, F1_adults.std_error, F2_adults.std_error, 929 opt$replications, life_stage, group=P_eggs, group_std_error=P_eggs.std_error, group2=F1_eggs, group2_std_error=F1_eggs.std_error, group3=F2_eggs,
671 date_labels); 930 group3_std_error=F2_eggs.std_error);
672 931 # Turn off device driver to flush output.
673 # Turn off device driver to flush output. 932 dev.off();
674 dev.off(); 933 } else if (life_stage == "Nymph") {
934 # Start PDF device driver.
935 dev.new(width=20, height=30);
936 file_name = paste(tolower(life_stage_nymph), "nymph_pop_by_generation.pdf", sep="_");
937 file_path = paste(output_dir, file_name, sep="/");
938 pdf(file=file_path, width=20, height=30, bg="white");
939 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
940 # Nymph population size by generation.
941 maxval = max(F2_nymphs) + 100;
942 render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
943 opt$replications, life_stage, group=P_nymphs, group_std_error=P_nymphs.std_error, group2=F1_nymphs, group2_std_error=F1_nymphs.std_error,
944 group3=F2_nymphs, group3_std_error=F2_nymphs.std_error, life_stage_nymph=life_stage_nymph);
945 # Turn off device driver to flush output.
946 dev.off();
947 } else if (life_stage == "Adult") {
948 # Start PDF device driver.
949 dev.new(width=20, height=30);
950 file_name = paste(tolower(life_stages_adult), "adult_pop_by_generation.pdf", sep="_");
951 file_path = paste(output_dir, file_name, sep="/");
952 pdf(file=file_path, width=20, height=30, bg="white");
953 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
954 # Adult population size by generation.
955 maxval = max(F2_adults) + 100;
956 render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
957 opt$replications, life_stage, group=P_adults, group_std_error=P_adults.std_error, group2=F1_adults, group2_std_error=F1_adults.std_error,
958 group3=F2_adults, group3_std_error=F2_adults.std_error, life_stages_adult=life_stages_adult);
959 # Turn off device driver to flush output.
960 dev.off();
961 } else if (life_stage == "Total") {
962 # Start PDF device driver.
963 dev.new(width=20, height=30);
964 file_path = paste(output_dir, "total_pop_by_generation.pdf", sep="/");
965 pdf(file=file_path, width=20, height=30, bg="white");
966 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
967 # Total population size by generation.
968 maxval = max(F2);
969 render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
970 opt$replications, life_stage, group=P, group_std_error=P.std_error, group2=F1, group2_std_error=F1.std_error, group3=F2, group3_std_error=F2.std_error);
971 # Turn off device driver to flush output.
972 dev.off();
973 }
974 } else {
975 for (life_stage in life_stages) {
976 if (life_stage == "Egg") {
977 # Start PDF device driver.
978 dev.new(width=20, height=30);
979 file_path = paste(output_dir, "egg_pop_size.pdf", sep="/");
980 pdf(file=file_path, width=20, height=30, bg="white");
981 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
982 # Egg population size.
983 maxval = max(eggs+eggs.std_error);
984 render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
985 opt$replications, life_stage, group=eggs, group_std_error=eggs.std_error);
986 # Turn off device driver to flush output.
987 dev.off();
988 } else if (life_stage == "Nymph") {
989 # Start PDF device driver.
990 dev.new(width=20, height=30);
991 file_name = paste(tolower(life_stage_nymph), "nymph_pop_size.pdf", sep="_");
992 file_path = paste(output_dir, file_name, sep="/");
993 pdf(file=file_path, width=20, height=30, bg="white");
994 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
995 # Nymph population size.
996 maxval = max(nymphs+nymphs.std_error);
997 render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
998 opt$replications, life_stage, group=nymphs, group_std_error=nymphs.std_error, life_stage_nymph=life_stage_nymph);
999 # Turn off device driver to flush output.
1000 dev.off();
1001 } else if (life_stage == "Adult") {
1002 # Start PDF device driver.
1003 dev.new(width=20, height=30);
1004 file_name = paste(tolower(life_stages_adult), "adult_pop_size.pdf", sep="_");
1005 file_path = paste(output_dir, file_name, sep="/");
1006 pdf(file=file_path, width=20, height=30, bg="white");
1007 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1008 # Adult population size.
1009 maxval = max(adults+adults.std_error);
1010 render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
1011 opt$replications, life_stage, group=adults, group_std_error=adults.std_error, life_stages_adult=life_stages_adult);
1012 # Turn off device driver to flush output.
1013 dev.off();
1014 } else if (life_stage == "Total") {
1015 # Start PDF device driver.
1016 dev.new(width=20, height=30);
1017 file_path = paste(output_dir, "total_pop_size.pdf", sep="/");
1018 pdf(file=file_path, width=20, height=30, bg="white");
1019 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1020 # Total population size.
1021 maxval = max(eggs+eggs.std_error, nymphs+nymphs.std_error, adults+adults.std_error);
1022 render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
1023 opt$replications, life_stage, group=adults, group_std_error=adults.std_error, group2=nymphs, group2_std_error=nymphs.std_error, group3=eggs,
1024 group3_std_error=eggs.std_error);
1025 # Turn off device driver to flush output.
1026 dev.off();
1027 }
1028 }
1029 }