Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 38:c0f76f4f84fc draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Tue, 10 Apr 2018 14:22:45 -0400
parents 5097cfeedc4f
children 169c8180205a
comparison
equal deleted inserted replaced
37:b7dcecf5476a 38:c0f76f4f84fc
4 4
5 option_list <- list( 5 option_list <- list(
6 make_option(c("--adult_mortality"), action="store", dest="adult_mortality", type="integer", help="Adjustment rate for adult mortality"), 6 make_option(c("--adult_mortality"), action="store", dest="adult_mortality", type="integer", help="Adjustment rate for adult mortality"),
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_norm"), action="store", dest="input_norm", help="30 year normals temperature data for selected station"),
10 make_option(c("--input_ytd"), action="store", dest="input_ytd", help="Year-to-date temperature data for selected location"),
10 make_option(c("--insect"), action="store", dest="insect", help="Insect name"), 11 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"), 12 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"), 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_stages_adult"), action="store", dest="life_stages_adult", default=NULL, help="Adult life stages for plotting"),
14 make_option(c("--life_stages_nymph"), action="store", dest="life_stages_nymph", default=NULL, help="Nymph life stages for plotting"), 15 make_option(c("--life_stages_nymph"), action="store", dest="life_stages_nymph", default=NULL, help="Nymph life stages for plotting"),
15 make_option(c("--location"), action="store", dest="location", help="Selected location"), 16 make_option(c("--location"), action="store", dest="location", help="Selected location"),
16 make_option(c("--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"), 17 make_option(c("--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"),
17 make_option(c("--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"), 18 make_option(c("--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),
18 make_option(c("--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"), 19 make_option(c("--num_days_ytd"), action="store", dest="num_days_ytd", type="integer", help="Total number of days in the temperature dataset"),
19 make_option(c("--nymph_mortality"), action="store", dest="nymph_mortality", type="integer", help="Adjustment rate for nymph mortality"), 20 make_option(c("--nymph_mortality"), action="store", dest="nymph_mortality", type="integer", help="Adjustment rate for nymph mortality"),
20 make_option(c("--old_nymph_accumulation"), action="store", dest="old_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (young nymph->old nymph)"), 21 make_option(c("--old_nymph_accumulation"), action="store", dest="old_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (young nymph->old nymph)"),
21 make_option(c("--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"), 22 make_option(c("--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"),
22 make_option(c("--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"), 23 make_option(c("--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"),
23 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"), 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"),
30 args <- parse_args(parser, positional_arguments=TRUE); 31 args <- parse_args(parser, positional_arguments=TRUE);
31 opt <- args$options; 32 opt <- args$options;
32 33
33 add_daylight_length = function(temperature_data_frame, num_rows) { 34 add_daylight_length = function(temperature_data_frame, num_rows) {
34 # Return a vector of daylight length (photoperido profile) for 35 # Return a vector of daylight length (photoperido profile) for
35 # the number of days specified in the input temperature data 36 # the number of days specified in the input_ytd temperature data
36 # (from Forsythe 1995). 37 # (from Forsythe 1995).
37 p = 0.8333; 38 p = 0.8333;
38 latitude = temperature_data_frame$LATITUDE[1]; 39 latitude = temperature_data_frame$LATITUDE[1];
39 daylight_length_vector = NULL; 40 daylight_length_vector = NULL;
40 for (i in 1:num_rows) { 41 for (i in 1:num_rows) {
216 } 217 }
217 } 218 }
218 return(length(ticks)+1); 219 return(length(ticks)+1);
219 } 220 }
220 221
221 get_x_axis_ticks_and_labels = function(temperature_data_frame, num_rows) { 222 get_total_days = function(is_leap_year) {
223 # Get the total number of days in the current year.
224 if (is_leap_year) {
225 return (366);
226 } else {
227 return (365);
228 }
229 }
230
231 get_x_axis_ticks_and_labels = function(temperature_data_frame, num_rows, num_days_ytd) {
222 # Keep track of the years to see if spanning years. 232 # Keep track of the years to see if spanning years.
223 month_labels = list(); 233 month_labels = list();
224 ticks = list(); 234 ticks = list();
225 current_month_label = NULL; 235 current_month_label = NULL;
226 last_tick = 0; 236 last_tick = 0;
227 for (i in 1:num_rows) { 237 for (i in 1:num_rows) {
238 if (i==num_days_ytd) {
239 # Add a tick for the start of the 30 year normnals data.
240 tick_index = get_tick_index(i, last_tick, ticks, month_labels)
241 ticks[tick_index] = i;
242 month_labels[tick_index] = "Start 30 year normals";
243 last_tick = i;
244 }
228 # Get the year and month from the date which 245 # Get the year and month from the date which
229 # has the format YYYY-MM-DD. 246 # has the format YYYY-MM-DD.
230 date = format(temperature_data_frame$DATE[i]); 247 date = format(temperature_data_frame$DATE[i]);
231 # Get the month label. 248 # Get the month label.
232 items = strsplit(date, "-")[[1]]; 249 items = strsplit(date, "-")[[1]];
233 month = items[2]; 250 month = items[2];
234 month_label = month.abb[as.integer(month)]; 251 month_label = month.abb[as.integer(month)];
235 tick_index = get_tick_index(i, last_tick, ticks, month_labels) 252 tick_index = get_tick_index(i, last_tick, ticks, month_labels)
253 if (!identical(current_month_label, month_label)) {
254 # Add an x-axis tick for the month.
255 ticks[tick_index] = i;
256 month_labels[tick_index] = month_label;
257 current_month_label = month_label;
258 last_tick = i;
259 }
260 tick_index = get_tick_index(i, last_tick, ticks, month_labels)
236 if (!is.null(tick_index)) { 261 if (!is.null(tick_index)) {
237 if (!identical(current_month_label, month_label)) {
238 # Add an x-axis tick for the month.
239 ticks[tick_index] = i;
240 month_labels[tick_index] = month_label;
241 current_month_label = month_label;
242 last_tick = i;
243 }
244 # Get the day. 262 # Get the day.
245 day = weekdays(as.Date(date)); 263 day = weekdays(as.Date(date));
246 if (day=="Sunday") { 264 if (day=="Sunday") {
247 # Add an x-axis tick if we're on a Sunday. 265 # Add an x-axis tick if we're on a Sunday.
248 ticks[tick_index] = i; 266 ticks[tick_index] = i;
249 # Add a blank month label so it is not displayed. 267 # Add a blank month label so it is not displayed.
250 month_labels[tick_index] = ""; 268 month_labels[tick_index] = "";
251 last_tick = i; 269 last_tick = i;
252 } 270 }
253 } 271 }
272 if (i==num_rows) {
273 # Add a tick for the last day of the year.
274 tick_index = get_tick_index(i, last_tick, ticks, month_labels)
275 ticks[tick_index] = i;
276 month_labels[tick_index] = "";
277 last_tick = i;
278 }
254 } 279 }
255 return(list(ticks, month_labels)); 280 return(list(ticks, month_labels));
281 }
282
283 is_leap_year = function(date_str) {
284 # Extract the year from the date_str.
285 date = format(date_str);
286 items = strsplit(date, "-")[[1]];
287 year = as.integer(items[1]);
288 if (((year %% 4 == 0) & (year %% 100 != 0)) | (year %% 400 == 0)) {
289 return (TRUE);
290 } else {
291 return (FALSE);
292 }
256 } 293 }
257 294
258 mortality.adult = function(temperature) { 295 mortality.adult = function(temperature) {
259 if (temperature < 12.7) { 296 if (temperature < 12.7) {
260 mortality.probability = 0.002; 297 mortality.probability = 0.002;
286 mortality.probability = temperature * 0.0008 + 0.03; 323 mortality.probability = temperature * 0.0008 + 0.03;
287 } 324 }
288 return(mortality.probability); 325 return(mortality.probability);
289 } 326 }
290 327
291 parse_input_data = function(input_file, num_rows) { 328 parse_input_data = function(input_ytd, input_norm, num_days_ytd) {
292 # Read in the input temperature datafile into a data frame. 329 # Read the input_ytd temperature datafile into a data frame.
293 temperature_data_frame = read.csv(file=input_file, header=T, strip.white=TRUE, sep=","); 330 # The input_ytd data has the following 6 columns:
294 num_columns = dim(temperature_data_frame)[2]; 331 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
295 if (num_columns == 6) { 332 temperature_data_frame = read.csv(file=input_ytd, header=T, strip.white=TRUE, stringsAsFactors=FALSE, sep=",");
296 # The input data has the following 6 columns: 333 # Set the temperature_data_frame column names for access.
297 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX 334 colnames(temperature_data_frame) = c("LATITUDE", "LONGITUDE", "DATE", "DOY", "TMIN", "TMAX");
298 # Set the column names for access when adding daylight length.. 335 # Get the start date.
299 colnames(temperature_data_frame) = c("LATITUDE","LONGITUDE", "DATE", "DOY", "TMIN", "TMAX"); 336 start_date = temperature_data_frame$DATE[1];
300 current_column_names = colnames(temperature_data_frame); 337 # See if we're in a leap year.
301 # Add a column containing the daylight length for each day. 338 is_leap_year = is_leap_year(start_date);
302 temperature_data_frame = add_daylight_length(temperature_data_frame, num_rows); 339 # get the number of days in the year.
303 } 340 total_days = get_total_days(is_leap_year);
341 # Extract the year from the start date.
342 date_str = format(start_date);
343 date_str_items = strsplit(date_str, "-")[[1]];
344 year = date_str_items[1];
345 # Read the input_norm temperature datafile into a data frame.
346 # The input_norm data has the following 10 columns:
347 # STATIONID, LATITUDE, LONGITUDE, ELEV_M, NAME, ST, MMDD, DOY, TMIN, TMAX
348 norm_data_frame = read.csv(file=input_norm, header=T, strip.white=TRUE, stringsAsFactors=FALSE, sep=",");
349 # Set the norm_data_frame column names for access.
350 colnames(norm_data_frame) = c("STATIONID", "LATITUDE","LONGITUDE", "ELEV_M", "NAME", "ST", "MMDD", "DOY", "TMIN", "TMAX");
351 # All normals data includes Feb 29 which is row 60 in
352 # the data, so delete that row if we're not in a leap year.
353 if (!is_leap_year) {
354 norm_data_frame = norm_data_frame[-c(60),];
355 }
356 # Define the next row for the temperature_data_frame from the 30 year normals data.
357 first_normals_row = num_days_ytd + 1;
358 for (i in first_normals_row:total_days) {
359 latitude = norm_data_frame[i,"LATITUDE"][1];
360 longitude = norm_data_frame[i,"LONGITUDE"][1];
361 # Format the date.
362 mmdd = norm_data_frame[i,"MMDD"][1];
363 date_str = paste(year, mmdd, sep="-");
364 doy = norm_data_frame[i,"DOY"][1];
365 if (!is_leap_year) {
366 # Since all normals data includes Feb 29, we have to
367 # subtract 1 from DOY if we're not in a leap year since
368 # we removed the Feb 29 row from the data frame above.
369 doy = as.integer(doy) - 1;
370 }
371 tmin = norm_data_frame[i,"TMIN"][1];
372 tmax = norm_data_frame[i,"TMAX"][1];
373 # Append the row to temperature_data_frame.
374 new_row = list(latitude, longitude, date_str, doy, tmin, tmax);
375 temperature_data_frame[i,] = new_row;
376 }
377 # Add a column containing the daylight length for each day.
378 temperature_data_frame = add_daylight_length(temperature_data_frame, total_days);
304 return(temperature_data_frame); 379 return(temperature_data_frame);
305 } 380 }
306 381
307 render_chart = function(ticks, date_labels, chart_type, plot_std_error, insect, location, latitude, start_date, end_date, days, maxval, 382 render_chart = function(ticks, date_labels, chart_type, plot_std_error, insect, location, latitude, start_date, end_date, days, maxval,
308 replications, life_stage, group, group_std_error, group2=NULL, group2_std_error=NULL, group3=NULL, group3_std_error=NULL, 383 replications, life_stage, group, group_std_error, group2=NULL, group2_std_error=NULL, group3=NULL, group3_std_error=NULL,
309 life_stages_adult=NULL, life_stages_nymph=NULL) { 384 life_stages_adult=NULL, life_stages_nymph=NULL) {
310 if (chart_type=="pop_size_by_life_stage") { 385 if (chart_type=="pop_size_by_life_stage") {
311 if (life_stage=="Total") { 386 if (life_stage=="Total") {
312 title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); 387 title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
313 legend_text = c("Egg", "Nymph", "Adult"); 388 legend_text = c("Egg", "Nymph", "Adult");
314 columns = c(4, 2, 1); 389 columns = c(4, 2, 1);
315 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=FALSE, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3); 390 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=FALSE, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3);
316 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3); 391 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3);
317 lines(days, group2, lwd=2, lty=1, col=2); 392 lines(days, group2, lwd=2, lty=1, col=2);
318 lines(days, group3, lwd=2, lty=1, col=4); 393 lines(days, group3, lwd=2, lty=1, col=4);
319 axis(side=1, at=ticks, labels=date_labels, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3); 394 axis(side=1, at=ticks, labels=date_labels, las=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
320 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3); 395 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
321 if (plot_std_error=="yes") { 396 if (plot_std_error=="yes") {
322 # Standard error for group. 397 # Standard error for group.
323 lines(days, group+group_std_error, lty=2); 398 lines(days, group+group_std_error, lty=2);
324 lines(days, group-group_std_error, lty=2); 399 lines(days, group-group_std_error, lty=2);
345 legend_text = c(paste(life_stages_adult, life_stage, sep=" ")); 420 legend_text = c(paste(life_stages_adult, life_stage, sep=" "));
346 columns = c(1); 421 columns = c(1);
347 } 422 }
348 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=FALSE, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3); 423 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=FALSE, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3);
349 legend("topleft", legend_text, lty=c(1), col="black", cex=3); 424 legend("topleft", legend_text, lty=c(1), col="black", cex=3);
350 axis(side=1, at=ticks, labels=date_labels, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3); 425 axis(side=1, at=ticks, labels=date_labels, las=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
351 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3); 426 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
352 if (plot_std_error=="yes") { 427 if (plot_std_error=="yes") {
353 # Standard error for group. 428 # Standard error for group.
354 lines(days, group+group_std_error, lty=2); 429 lines(days, group+group_std_error, lty=2);
355 lines(days, group-group_std_error, lty=2); 430 lines(days, group-group_std_error, lty=2);
370 columns = c(1, 2, 4); 445 columns = c(1, 2, 4);
371 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=FALSE, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3); 446 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=FALSE, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3);
372 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3); 447 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3);
373 lines(days, group2, lwd=2, lty=1, col=2); 448 lines(days, group2, lwd=2, lty=1, col=2);
374 lines(days, group3, lwd=2, lty=1, col=4); 449 lines(days, group3, lwd=2, lty=1, col=4);
375 axis(side=1, at=ticks, labels=date_labels, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3); 450 axis(side=1, at=ticks, labels=date_labels, las=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
376 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3); 451 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
377 if (plot_std_error=="yes") { 452 if (plot_std_error=="yes") {
378 # Standard error for group. 453 # Standard error for group.
379 lines(days, group+group_std_error, lty=2); 454 lines(days, group+group_std_error, lty=2);
380 lines(days, group-group_std_error, lty=2); 455 lines(days, group-group_std_error, lty=2);
392 if (opt$plot_generations_separately=="yes") { 467 if (opt$plot_generations_separately=="yes") {
393 plot_generations_separately = TRUE; 468 plot_generations_separately = TRUE;
394 } else { 469 } else {
395 plot_generations_separately = FALSE; 470 plot_generations_separately = FALSE;
396 } 471 }
472 # Display the total number of days in the Galaxy history item blurb.
473 cat("Year-to-date number of days: ", opt$num_days_ytd, "\n");
474
397 # Read the temperature data into a data frame. 475 # Read the temperature data into a data frame.
398 temperature_data_frame = parse_input_data(opt$input, opt$num_days); 476 temperature_data_frame = parse_input_data(opt$input_ytd, opt$input_norm, opt$num_days_ytd);
477
399 # Create copies of the temperature data for generations P, F1 and F2 if we're plotting generations separately. 478 # Create copies of the temperature data for generations P, F1 and F2 if we're plotting generations separately.
400 if (plot_generations_separately) { 479 if (plot_generations_separately) {
401 temperature_data_frame_P = data.frame(temperature_data_frame); 480 temperature_data_frame_P = data.frame(temperature_data_frame);
402 temperature_data_frame_F1 = data.frame(temperature_data_frame); 481 temperature_data_frame_F1 = data.frame(temperature_data_frame);
403 temperature_data_frame_F2 = data.frame(temperature_data_frame); 482 temperature_data_frame_F2 = data.frame(temperature_data_frame);
404 } 483 }
405 # Get the date labels for plots. 484
406 ticks_and_labels = get_x_axis_ticks_and_labels(temperature_data_frame, opt$num_days); 485 # Information needed for plots.
486 start_date = temperature_data_frame$DATE[1];
487 end_date = temperature_data_frame$DATE[opt$num_days_ytd];
488 # See if we're in a leap year.
489 is_leap_year = is_leap_year(start_date);
490 total_days = get_total_days(is_leap_year);
491 total_days_vector = c(1:total_days);
492
493 # Get the ticks date labels for plots.
494 ticks_and_labels = get_x_axis_ticks_and_labels(temperature_data_frame, total_days, opt$num_days_ytd);
407 ticks = c(unlist(ticks_and_labels[1])); 495 ticks = c(unlist(ticks_and_labels[1]));
408 date_labels = c(unlist(ticks_and_labels[2])); 496 date_labels = c(unlist(ticks_and_labels[2]));
409 # All latitude values are the same, so get the value for plots from the first row. 497 # All latitude values are the same, so get the value for plots from the first row.
410 latitude = temperature_data_frame$LATITUDE[1]; 498 latitude = temperature_data_frame$LATITUDE[1];
499
411 # Determine the specified life stages for processing. 500 # Determine the specified life stages for processing.
412 # Split life_stages into a list of strings for plots. 501 # Split life_stages into a list of strings for plots.
413 life_stages_str = as.character(opt$life_stages); 502 life_stages_str = as.character(opt$life_stages);
414 life_stages = strsplit(life_stages_str, ",")[[1]]; 503 life_stages = strsplit(life_stages_str, ",")[[1]];
504
415 # Determine the data we need to generate for plotting. 505 # Determine the data we need to generate for plotting.
416 process_eggs = FALSE; 506 process_eggs = FALSE;
417 process_nymphs = FALSE; 507 process_nymphs = FALSE;
418 process_young_nymphs = FALSE; 508 process_young_nymphs = FALSE;
419 process_old_nymphs = FALSE; 509 process_old_nymphs = FALSE;
466 } 556 }
467 } 557 }
468 } 558 }
469 # Initialize matrices. 559 # Initialize matrices.
470 if (process_eggs) { 560 if (process_eggs) {
471 Eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 561 Eggs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
472 } 562 }
473 if (process_young_nymphs | process_total_nymphs) { 563 if (process_young_nymphs | process_total_nymphs) {
474 YoungNymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 564 YoungNymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
475 } 565 }
476 if (process_old_nymphs | process_total_nymphs) { 566 if (process_old_nymphs | process_total_nymphs) {
477 OldNymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 567 OldNymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
478 } 568 }
479 if (process_previttelogenic_adults | process_total_adults) { 569 if (process_previttelogenic_adults | process_total_adults) {
480 Previttelogenic.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 570 Previttelogenic.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
481 } 571 }
482 if (process_vittelogenic_adults | process_total_adults) { 572 if (process_vittelogenic_adults | process_total_adults) {
483 Vittelogenic.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 573 Vittelogenic.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
484 } 574 }
485 if (process_diapausing_adults | process_total_adults) { 575 if (process_diapausing_adults | process_total_adults) {
486 Diapausing.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 576 Diapausing.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
487 } 577 }
488 newborn.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 578 newborn.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
489 adult.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 579 adult.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
490 death.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 580 death.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
491 if (plot_generations_separately) { 581 if (plot_generations_separately) {
492 # P is Parental, or overwintered adults. 582 # P is Parental, or overwintered adults.
493 P.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 583 P.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
494 # F1 is the first field-produced generation. 584 # F1 is the first field-produced generation.
495 F1.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 585 F1.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
496 # F2 is the second field-produced generation. 586 # F2 is the second field-produced generation.
497 F2.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 587 F2.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
498 if (process_eggs) { 588 if (process_eggs) {
499 P_eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 589 P_eggs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
500 F1_eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 590 F1_eggs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
501 F2_eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 591 F2_eggs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
502 } 592 }
503 if (process_young_nymphs) { 593 if (process_young_nymphs) {
504 P_young_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 594 P_young_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
505 F1_young_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 595 F1_young_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
506 F2_young_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 596 F2_young_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
507 } 597 }
508 if (process_old_nymphs) { 598 if (process_old_nymphs) {
509 P_old_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 599 P_old_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
510 F1_old_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 600 F1_old_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
511 F2_old_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 601 F2_old_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
512 } 602 }
513 if (process_total_nymphs) { 603 if (process_total_nymphs) {
514 P_total_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 604 P_total_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
515 F1_total_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 605 F1_total_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
516 F2_total_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 606 F2_total_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
517 } 607 }
518 if (process_previttelogenic_adults) { 608 if (process_previttelogenic_adults) {
519 P_previttelogenic_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 609 P_previttelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
520 F1_previttelogenic_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 610 F1_previttelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
521 F2_previttelogenic_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 611 F2_previttelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
522 } 612 }
523 if (process_vittelogenic_adults) { 613 if (process_vittelogenic_adults) {
524 P_vittelogenic_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 614 P_vittelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
525 F1_vittelogenic_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 615 F1_vittelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
526 F2_vittelogenic_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 616 F2_vittelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
527 } 617 }
528 if (process_diapausing_adults) { 618 if (process_diapausing_adults) {
529 P_diapausing_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 619 P_diapausing_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
530 F1_diapausing_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 620 F1_diapausing_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
531 F2_diapausing_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 621 F2_diapausing_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
532 } 622 }
533 if (process_total_adults) { 623 if (process_total_adults) {
534 P_total_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 624 P_total_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
535 F1_total_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 625 F1_total_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
536 F2_total_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 626 F2_total_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
537 } 627 }
538 } 628 }
539 # Total population. 629 # Total population.
540 population.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); 630 population.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
541 631
542 # Process replications. 632 # Process replications.
543 for (current_replication in 1:opt$replications) { 633 for (current_replication in 1:opt$replications) {
544 # Start with the user-defined number of insects per replication. 634 # Start with the user-defined number of insects per replication.
545 num_insects = opt$insects_per_replication; 635 num_insects = opt$insects_per_replication;
552 # Complete transposed matrix for the population, so now 642 # Complete transposed matrix for the population, so now
553 # the rows are Generation, Stage, degree-days, T, Diapause 643 # the rows are Generation, Stage, degree-days, T, Diapause
554 vector.matrix = base::t(matrix(vector.matrix, nrow=5)); 644 vector.matrix = base::t(matrix(vector.matrix, nrow=5));
555 # Time series of population size. 645 # Time series of population size.
556 if (process_eggs) { 646 if (process_eggs) {
557 Eggs = rep(0, opt$num_days); 647 Eggs = rep(0, total_days);
558 } 648 }
559 if (process_young_nymphs | process_total_nymphs) { 649 if (process_young_nymphs | process_total_nymphs) {
560 YoungNymphs = rep(0, opt$num_days); 650 YoungNymphs = rep(0, total_days);
561 } 651 }
562 if (process_old_nymphs | process_total_nymphs) { 652 if (process_old_nymphs | process_total_nymphs) {
563 OldNymphs = rep(0, opt$num_days); 653 OldNymphs = rep(0, total_days);
564 } 654 }
565 if (process_previttelogenic_adults | process_total_adults) { 655 if (process_previttelogenic_adults | process_total_adults) {
566 Previttelogenic = rep(0, opt$num_days); 656 Previttelogenic = rep(0, total_days);
567 } 657 }
568 if (process_vittelogenic_adults | process_total_adults) { 658 if (process_vittelogenic_adults | process_total_adults) {
569 Vittelogenic = rep(0, opt$num_days); 659 Vittelogenic = rep(0, total_days);
570 } 660 }
571 if (process_diapausing_adults | process_total_adults) { 661 if (process_diapausing_adults | process_total_adults) {
572 Diapausing = rep(0, opt$num_days); 662 Diapausing = rep(0, total_days);
573 } 663 }
574 N.newborn = rep(0, opt$num_days); 664 N.newborn = rep(0, total_days);
575 N.adult = rep(0, opt$num_days); 665 N.adult = rep(0, total_days);
576 N.death = rep(0, opt$num_days); 666 N.death = rep(0, total_days);
577 overwintering_adult.population = rep(0, opt$num_days); 667 overwintering_adult.population = rep(0, total_days);
578 first_generation.population = rep(0, opt$num_days); 668 first_generation.population = rep(0, total_days);
579 second_generation.population = rep(0, opt$num_days); 669 second_generation.population = rep(0, total_days);
580 if (plot_generations_separately) { 670 if (plot_generations_separately) {
581 # P is Parental, or overwintered adults. 671 # P is Parental, or overwintered adults.
582 # F1 is the first field-produced generation. 672 # F1 is the first field-produced generation.
583 # F2 is the second field-produced generation. 673 # F2 is the second field-produced generation.
584 if (process_eggs) { 674 if (process_eggs) {
585 P.egg = rep(0, opt$num_days); 675 P.egg = rep(0, total_days);
586 F1.egg = rep(0, opt$num_days); 676 F1.egg = rep(0, total_days);
587 F2.egg = rep(0, opt$num_days); 677 F2.egg = rep(0, total_days);
588 } 678 }
589 if (process_young_nymphs) { 679 if (process_young_nymphs) {
590 P.young_nymph = rep(0, opt$num_days); 680 P.young_nymph = rep(0, total_days);
591 F1.young_nymph = rep(0, opt$num_days); 681 F1.young_nymph = rep(0, total_days);
592 F2.young_nymph = rep(0, opt$num_days); 682 F2.young_nymph = rep(0, total_days);
593 } 683 }
594 if (process_old_nymphs) { 684 if (process_old_nymphs) {
595 P.old_nymph = rep(0, opt$num_days); 685 P.old_nymph = rep(0, total_days);
596 F1.old_nymph = rep(0, opt$num_days); 686 F1.old_nymph = rep(0, total_days);
597 F2.old_nymph = rep(0, opt$num_days); 687 F2.old_nymph = rep(0, total_days);
598 } 688 }
599 if (process_total_nymphs) { 689 if (process_total_nymphs) {
600 P.total_nymph = rep(0, opt$num_days); 690 P.total_nymph = rep(0, total_days);
601 F1.total_nymph = rep(0, opt$num_days); 691 F1.total_nymph = rep(0, total_days);
602 F2.total_nymph = rep(0, opt$num_days); 692 F2.total_nymph = rep(0, total_days);
603 } 693 }
604 if (process_previttelogenic_adults) { 694 if (process_previttelogenic_adults) {
605 P.previttelogenic_adult = rep(0, opt$num_days); 695 P.previttelogenic_adult = rep(0, total_days);
606 F1.previttelogenic_adult = rep(0, opt$num_days); 696 F1.previttelogenic_adult = rep(0, total_days);
607 F2.previttelogenic_adult = rep(0, opt$num_days); 697 F2.previttelogenic_adult = rep(0, total_days);
608 } 698 }
609 if (process_vittelogenic_adults) { 699 if (process_vittelogenic_adults) {
610 P.vittelogenic_adult = rep(0, opt$num_days); 700 P.vittelogenic_adult = rep(0, total_days);
611 F1.vittelogenic_adult = rep(0, opt$num_days); 701 F1.vittelogenic_adult = rep(0, total_days);
612 F2.vittelogenic_adult = rep(0, opt$num_days); 702 F2.vittelogenic_adult = rep(0, total_days);
613 } 703 }
614 if (process_diapausing_adults) { 704 if (process_diapausing_adults) {
615 P.diapausing_adult = rep(0, opt$num_days); 705 P.diapausing_adult = rep(0, total_days);
616 F1.diapausing_adult = rep(0, opt$num_days); 706 F1.diapausing_adult = rep(0, total_days);
617 F2.diapausing_adult = rep(0, opt$num_days); 707 F2.diapausing_adult = rep(0, total_days);
618 } 708 }
619 if (process_total_adults) { 709 if (process_total_adults) {
620 P.total_adult = rep(0, opt$num_days); 710 P.total_adult = rep(0, total_days);
621 F1.total_adult = rep(0, opt$num_days); 711 F1.total_adult = rep(0, total_days);
622 F2.total_adult = rep(0, opt$num_days); 712 F2.total_adult = rep(0, total_days);
623 } 713 }
624 } 714 }
625 total.population = NULL; 715 total.population = NULL;
626 averages.day = rep(0, opt$num_days); 716 averages.day = rep(0, total_days);
627 # All the days included in the input temperature dataset. 717 # All the days included in the input_ytd temperature dataset.
628 for (row in 1:opt$num_days) { 718 for (row in 1:total_days) {
629 # Get the integer day of the year for the current row. 719 # Get the integer day of the year for the current row.
630 doy = temperature_data_frame$DOY[row]; 720 doy = temperature_data_frame$DOY[row];
631 # Photoperiod in the day. 721 # Photoperiod in the day.
632 photoperiod = temperature_data_frame$DAYLEN[row]; 722 photoperiod = temperature_data_frame$DAYLEN[row];
633 temp.profile = get_temperature_at_hour(latitude, temperature_data_frame, row, opt$num_days); 723 temp.profile = get_temperature_at_hour(latitude, temperature_data_frame, row, total_days);
634 mean.temp = temp.profile[1]; 724 mean.temp = temp.profile[1];
635 averages.temp = temp.profile[2]; 725 averages.temp = temp.profile[2];
636 averages.day[row] = averages.temp; 726 averages.day[row] = averages.temp;
637 # Trash bin for death. 727 # Trash bin for death.
638 death.vector = NULL; 728 death.vector = NULL;
1016 # - column 1 (Generation) is 2 and column 2 (Stage) is 4 (Vittelogenic) 1106 # - column 1 (Generation) is 2 and column 2 (Stage) is 4 (Vittelogenic)
1017 # - column 1 (Generation) is 2 and column 2 (Stage) is 5 (Diapausing) 1107 # - column 1 (Generation) is 2 and column 2 (Stage) is 5 (Diapausing)
1018 F2.total_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)); 1108 F2.total_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));
1019 } 1109 }
1020 } 1110 }
1021 } # End of days specified in the input temperature data. 1111 } # End of days specified in the input_ytd temperature data.
1022 1112
1023 averages.cum = cumsum(averages.day); 1113 averages.cum = cumsum(averages.day);
1024 1114
1025 # Define the output values. 1115 # Define the output values.
1026 if (process_eggs) { 1116 if (process_eggs) {
1308 write.csv(temperature_data_frame_F1, file=file_path, row.names=F); 1398 write.csv(temperature_data_frame_F1, file=file_path, row.names=F);
1309 # Save the analyzed data for generation F2. 1399 # Save the analyzed data for generation F2.
1310 file_path = paste("output_data_dir", "03_generation_F2.csv", sep="/"); 1400 file_path = paste("output_data_dir", "03_generation_F2.csv", sep="/");
1311 write.csv(temperature_data_frame_F2, file=file_path, row.names=F); 1401 write.csv(temperature_data_frame_F2, file=file_path, row.names=F);
1312 } 1402 }
1313 # Display the total number of days in the Galaxy history item blurb.
1314 cat("Number of days: ", opt$num_days, "\n");
1315 # Information needed for plots plots.
1316 days = c(1:opt$num_days);
1317 start_date = temperature_data_frame$DATE[1];
1318 end_date = temperature_data_frame$DATE[opt$num_days];
1319 1403
1320 if (plot_generations_separately) { 1404 if (plot_generations_separately) {
1321 for (life_stage in life_stages) { 1405 for (life_stage in life_stages) {
1322 if (life_stage == "Egg") { 1406 if (life_stage == "Egg") {
1323 # Start PDF device driver. 1407 # Start PDF device driver.
1325 file_path = get_file_path(life_stage, "egg_pop_by_generation.pdf") 1409 file_path = get_file_path(life_stage, "egg_pop_by_generation.pdf")
1326 pdf(file=file_path, width=20, height=30, bg="white"); 1410 pdf(file=file_path, width=20, height=30, bg="white");
1327 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1411 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1328 # Egg population size by generation. 1412 # Egg population size by generation.
1329 maxval = max(P_eggs+F1_eggs+F2_eggs) + 100; 1413 maxval = max(P_eggs+F1_eggs+F2_eggs) + 100;
1330 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, 1414 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude,
1331 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, 1415 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=P_eggs, group_std_error=P_eggs.std_error,
1332 group3_std_error=F2_eggs.std_error); 1416 group2=F1_eggs, group2_std_error=F1_eggs.std_error, group3=F2_eggs, group3_std_error=F2_eggs.std_error);
1333 # Turn off device driver to flush output. 1417 # Turn off device driver to flush output.
1334 dev.off(); 1418 dev.off();
1335 } else if (life_stage == "Nymph") { 1419 } else if (life_stage == "Nymph") {
1336 for (life_stage_nymph in life_stages_nymph) { 1420 for (life_stage_nymph in life_stages_nymph) {
1337 # Start PDF device driver. 1421 # Start PDF device driver.
1365 group2 = F1_total_nymphs; 1449 group2 = F1_total_nymphs;
1366 group2_std_error = F1_total_nymphs.std_error; 1450 group2_std_error = F1_total_nymphs.std_error;
1367 group3 = F2_total_nymphs; 1451 group3 = F2_total_nymphs;
1368 group3_std_error = F2_total_nymphs.std_error; 1452 group3_std_error = F2_total_nymphs.std_error;
1369 } 1453 }
1370 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, 1454 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude,
1371 opt$replications, life_stage, group=group, group_std_error=group_std_error, group2=group2, group2_std_error=group2_std_error, 1455 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
1372 group3=group3, group3_std_error=group3_std_error, life_stages_nymph=life_stage_nymph); 1456 group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, life_stages_nymph=life_stage_nymph);
1373 # Turn off device driver to flush output. 1457 # Turn off device driver to flush output.
1374 dev.off(); 1458 dev.off();
1375 } 1459 }
1376 } else if (life_stage == "Adult") { 1460 } else if (life_stage == "Adult") {
1377 for (life_stage_adult in life_stages_adult) { 1461 for (life_stage_adult in life_stages_adult) {
1415 group2 = F1_total_adults; 1499 group2 = F1_total_adults;
1416 group2_std_error = F1_total_adults.std_error; 1500 group2_std_error = F1_total_adults.std_error;
1417 group3 = F2_total_adults; 1501 group3 = F2_total_adults;
1418 group3_std_error = F2_total_adults.std_error; 1502 group3_std_error = F2_total_adults.std_error;
1419 } 1503 }
1420 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, 1504 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude,
1421 opt$replications, life_stage, group=group, group_std_error=group_std_error, group2=group2, group2_std_error=group2_std_error, 1505 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
1422 group3=group3, group3_std_error=group3_std_error, life_stages_adult=life_stage_adult); 1506 group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, life_stages_adult=life_stage_adult);
1423 # Turn off device driver to flush output. 1507 # Turn off device driver to flush output.
1424 dev.off(); 1508 dev.off();
1425 } 1509 }
1426 } else if (life_stage == "Total") { 1510 } else if (life_stage == "Total") {
1427 # Start PDF device driver. 1511 # Start PDF device driver.
1431 file_path = get_file_path(life_stage, "total_pop_by_generation.pdf") 1515 file_path = get_file_path(life_stage, "total_pop_by_generation.pdf")
1432 pdf(file=file_path, width=20, height=30, bg="white"); 1516 pdf(file=file_path, width=20, height=30, bg="white");
1433 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1517 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1434 # Total population size by generation. 1518 # Total population size by generation.
1435 maxval = max(P+F1+F2) + 100; 1519 maxval = max(P+F1+F2) + 100;
1436 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, 1520 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude,
1437 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); 1521 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=P, group_std_error=P.std_error,
1522 group2=F1, group2_std_error=F1.std_error, group3=F2, group3_std_error=F2.std_error);
1438 # Turn off device driver to flush output. 1523 # Turn off device driver to flush output.
1439 dev.off(); 1524 dev.off();
1440 } 1525 }
1441 } 1526 }
1442 } else { 1527 } else {
1447 file_path = get_file_path(life_stage, "egg_pop.pdf") 1532 file_path = get_file_path(life_stage, "egg_pop.pdf")
1448 pdf(file=file_path, width=20, height=30, bg="white"); 1533 pdf(file=file_path, width=20, height=30, bg="white");
1449 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1534 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1450 # Egg population size. 1535 # Egg population size.
1451 maxval = max(eggs+eggs.std_error) + 100; 1536 maxval = max(eggs+eggs.std_error) + 100;
1452 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, 1537 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude,
1453 opt$replications, life_stage, group=eggs, group_std_error=eggs.std_error); 1538 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=eggs, group_std_error=eggs.std_error);
1454 # Turn off device driver to flush output. 1539 # Turn off device driver to flush output.
1455 dev.off(); 1540 dev.off();
1456 } else if (life_stage == "Nymph") { 1541 } else if (life_stage == "Nymph") {
1457 for (life_stage_nymph in life_stages_nymph) { 1542 for (life_stage_nymph in life_stages_nymph) {
1458 # Start PDF device driver. 1543 # Start PDF device driver.
1472 # Old nymph population size. 1557 # Old nymph population size.
1473 group = old_nymphs; 1558 group = old_nymphs;
1474 group_std_error = old_nymphs.std_error; 1559 group_std_error = old_nymphs.std_error;
1475 } 1560 }
1476 maxval = max(group+group_std_error) + 100; 1561 maxval = max(group+group_std_error) + 100;
1477 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, 1562 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude,
1478 opt$replications, life_stage, group=group, group_std_error=group_std_error, life_stages_nymph=life_stage_nymph); 1563 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
1564 life_stages_nymph=life_stage_nymph);
1479 # Turn off device driver to flush output. 1565 # Turn off device driver to flush output.
1480 dev.off(); 1566 dev.off();
1481 } 1567 }
1482 } else if (life_stage == "Adult") { 1568 } else if (life_stage == "Adult") {
1483 for (life_stage_adult in life_stages_adult) { 1569 for (life_stage_adult in life_stages_adult) {
1502 # Diapausing adult population size. 1588 # Diapausing adult population size.
1503 group = diapausing_adults; 1589 group = diapausing_adults;
1504 group_std_error = diapausing_adults.std_error 1590 group_std_error = diapausing_adults.std_error
1505 } 1591 }
1506 maxval = max(group+group_std_error) + 100; 1592 maxval = max(group+group_std_error) + 100;
1507 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, 1593 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude,
1508 opt$replications, life_stage, group=group, group_std_error=group_std_error, life_stages_adult=life_stage_adult); 1594 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
1595 life_stages_adult=life_stage_adult);
1509 # Turn off device driver to flush output. 1596 # Turn off device driver to flush output.
1510 dev.off(); 1597 dev.off();
1511 } 1598 }
1512 } else if (life_stage == "Total") { 1599 } else if (life_stage == "Total") {
1513 # Start PDF device driver. 1600 # Start PDF device driver.
1515 file_path = get_file_path(life_stage, "total_pop.pdf") 1602 file_path = get_file_path(life_stage, "total_pop.pdf")
1516 pdf(file=file_path, width=20, height=30, bg="white"); 1603 pdf(file=file_path, width=20, height=30, bg="white");
1517 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1604 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1518 # Total population size. 1605 # Total population size.
1519 maxval = max(eggs+eggs.std_error, total_nymphs+total_nymphs.std_error, total_adults+total_adults.std_error) + 100; 1606 maxval = max(eggs+eggs.std_error, total_nymphs+total_nymphs.std_error, total_adults+total_adults.std_error) + 100;
1520 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, 1607 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude,
1521 opt$replications, life_stage, group=total_adults, group_std_error=total_adults.std_error, group2=total_nymphs, group2_std_error=total_nymphs.std_error, group3=eggs, 1608 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=total_adults, group_std_error=total_adults.std_error,
1522 group3_std_error=eggs.std_error); 1609 group2=total_nymphs, group2_std_error=total_nymphs.std_error, group3=eggs, group3_std_error=eggs.std_error);
1523 # Turn off device driver to flush output. 1610 # Turn off device driver to flush output.
1524 dev.off(); 1611 dev.off();
1525 } 1612 }
1526 } 1613 }
1527 } 1614 }