insect_phenology_model: insect_phenology

comparison insect_phenology_model.R @ 112:bcb12b7e8563 draft

Uploaded

author	greg
date	Tue, 29 May 2018 09:00:25 -0400
parents	37ac68b6ff10
children	9c998fd06628

comparison

equal deleted inserted replaced

-:37ac68b6ff10
+:bcb12b7e8563
 option_list <- list(
 make_option(c("--adult_mortality"), action="store", dest="adult_mortality", type="integer", help="Adjustment rate for adult mortality"),
 make_option(c("--adult_accumulation"), action="store", dest="adult_accumulation", type="integer", help="Adjustment of degree-days accumulation (old nymph->adult)"),
 make_option(c("--egg_mortality"), action="store", dest="egg_mortality", type="integer", help="Adjustment rate for egg mortality"),
-make_option(c("--input"), action="store", dest="input", help="Temperature data for selected location"),
+make_option(c("--input_norm"), action="store", dest="input_norm", help="30 year normals temperature data for selected station"),
+make_option(c("--input_ytd"), action="store", dest="input_ytd", default=NULL, help="Year-to-date temperature data for selected location"),
 make_option(c("--insect"), action="store", dest="insect", help="Insect name"),
 make_option(c("--insects_per_replication"), action="store", dest="insects_per_replication", type="integer", help="Number of insects with which to start each replication"),
-make_option(c("--location"), action="store", dest="location", help="Selected location"),
+make_option(c("--life_stages"), action="store", dest="life_stages", help="Selected life stages for plotting"),
+make_option(c("--life_stages_adult"), action="store", dest="life_stages_adult", default=NULL, help="Adult life stages for plotting"),
+make_option(c("--life_stages_nymph"), action="store", dest="life_stages_nymph", default=NULL, help="Nymph life stages for plotting"),
+make_option(c("--location"), action="store", dest="location", default=NULL, help="Selected location"),
 make_option(c("--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"),
 make_option(c("--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),
+make_option(c("--num_days_ytd"), action="store", dest="num_days_ytd", default=NULL, type="integer", help="Total number of days in the year-to-date temperature dataset"),
 make_option(c("--nymph_mortality"), action="store", dest="nymph_mortality", type="integer", help="Adjustment rate for nymph mortality"),
 make_option(c("--old_nymph_accumulation"), action="store", dest="old_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (young nymph->old nymph)"),
-make_option(c("--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"),
-make_option(c("--output"), action="store", dest="output", help="Output dataset"),
 make_option(c("--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"),
 make_option(c("--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"),
+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"),
+make_option(c("--plot_std_error"), action="store", dest="plot_std_error", help="Plot Standard error"),
 make_option(c("--replications"), action="store", dest="replications", type="integer", help="Number of replications"),
-make_option(c("--std_error_plot"), action="store", dest="std_error_plot", help="Plot Standard error"),
 make_option(c("--young_nymph_accumulation"), action="store", dest="young_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)")
 )
 parser <- OptionParser(usage="%prog [options] file", option_list=option_list);
 args <- parse_args(parser, positional_arguments=TRUE);
 opt <- args$options;
-add_daylight_length = function(temperature_data_frame, num_columns, num_rows) {
+add_daylight_length = function(temperature_data_frame, num_rows) {
 # Return a vector of daylight length (photoperido profile) for
-# the number of days specified in the input temperature data
+# the number of days specified in the input_ytd temperature data
 # (from Forsythe 1995).
 p = 0.8333;
 latitude = temperature_data_frame$LATITUDE[1];
 daylight_length_vector = NULL;
 for (i in 1:num_rows) {
 # Compute the length of daylight for the day of the year.
 darkness_length = 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi)));
 daylight_length_vector[i] = 24 - darkness_length;
 }
 # Append daylight_length_vector as a new column to temperature_data_frame.
-temperature_data_frame[, num_columns+1] = daylight_length_vector;
+temperature_data_frame = append_vector(temperature_data_frame, daylight_length_vector, "DAYLEN");
 return(temperature_data_frame);
 }
-dev.egg = function(temperature) {
+append_vector = function(data_frame, vec, new_column_name) {
-dev.rate = -0.9843 * temperature + 33.438;
+num_columns = dim(data_frame)[2];
-return(dev.rate);
+current_column_names = colnames(data_frame);
-}
+# Append vector vec as a new column to data_frame.
+data_frame[,num_columns+1] = vec;
-dev.emerg = function(temperature) {
+# Reset the column names with the additional column for later access.
-emerg.rate = -0.5332 * temperature + 24.147;
+colnames(data_frame) = append(current_column_names, new_column_name);
-return(emerg.rate);
+return(data_frame);
 }
-dev.old = function(temperature) {
+get_file_path = function(life_stage, base_name, life_stage_nymph=NULL, life_stage_adult=NULL) {
-n34 = -0.6119 * temperature + 17.602;
+if (!is.null(life_stage_nymph)) {
-n45 = -0.4408 * temperature + 19.036;
+lsi = get_life_stage_index(life_stage, life_stage_nymph=life_stage_nymph);
-dev.rate = mean(n34 + n45);
+file_name = paste(lsi, tolower(life_stage_nymph), base_name, sep="_");
-return(dev.rate);
+} else if (!is.null(life_stage_adult)) {
-}
+lsi = get_life_stage_index(life_stage, life_stage_adult=life_stage_adult);
+file_name = paste(lsi, tolower(life_stage_adult), base_name, sep="_");
-dev.young = function(temperature) {
+} else {
-n12 = -0.3728 * temperature + 14.68;
+lsi = get_life_stage_index(life_stage);
-n23 = -0.6119 * temperature + 25.249;
+file_name = paste(lsi, base_name, sep="_");
-dev.rate = mean(n12 + n23);
+}
-return(dev.rate);
+file_path = paste("output_plots_dir", file_name, sep="/");
-}
+return(file_path);
+}
-get_date_labels = function(temperature_data_frame, num_rows) {
+get_life_stage_index = function(life_stage, life_stage_nymph=NULL, life_stage_adult=NULL) {
-# Keep track of the years to see if spanning years.
+# Name collection elements so that they
-month_labels = list();
+# are displayed in logical order.
-current_month_label = NULL;
+if (life_stage=="Egg") {
-for (i in 1:num_rows) {
+lsi = "01";
-# Get the year and month from the date which
+} else if (life_stage=="Nymph") {
-# has the format YYYY-MM-DD.
+if (life_stage_nymph=="Young") {
-date = format(temperature_data_frame$DATE[i]);
+lsi = "02";
-items = strsplit(date, "-")[[1]];
+} else if (life_stage_nymph=="Old") {
-month = items[2];
+lsi = "03";
-month_label = month.abb[as.integer(month)];
+} else if (life_stage_nymph=="Total") {
-if (!identical(current_month_label, month_label)) {
+lsi="04";
-month_labels[length(month_labels)+1] = month_label;
+}
-current_month_label = month_label;
+} else if (life_stage=="Adult") {
-}
+if (life_stage_adult=="Pre-vittelogenic") {
-}
+lsi = "05";
-return(c(unlist(month_labels)));
+} else if (life_stage_adult=="Vittelogenic") {
+lsi = "06";
+} else if (life_stage_adult=="Diapausing") {
+lsi = "07";
+} else if (life_stage_adult=="Total") {
+lsi = "08";
+}
+} else if (life_stage=="Total") {
+lsi = "09";
+}
+return(lsi);
+}
+get_mean_and_std_error = function(p_replications, f1_replications, f2_replications) {
+# P mean.
+p_m = apply(p_replications, 1, mean);
+# P standard error.
+p_se = apply(p_replications, 1, sd) / sqrt(opt$replications);
+# F1 mean.
+f1_m = apply(f1_replications, 1, mean);
+# F1 standard error.
+f1_se = apply(f1_replications, 1, sd) / sqrt(opt$replications);
+# F2 mean.
+f2_m = apply(f2_replications, 1, mean);
+# F2 standard error.
+f2_se = apply(f2_replications, 1, sd) / sqrt(opt$replications);
+return(list(p_m, p_se, f1_m, f1_se, f2_m, f2_se))
+}
+get_next_normals_row = function(norm_data_frame, year, is_leap_year, index) {
+# Return the next 30 year normals row formatted
+# appropriately for the year-to-date data.
+latitude = norm_data_frame[index,"LATITUDE"][1];
+longitude = norm_data_frame[index,"LONGITUDE"][1];
+# Format the date.
+mmdd = norm_data_frame[index,"MMDD"][1];
+date_str = paste(year, mmdd, sep="-");
+doy = norm_data_frame[index,"DOY"][1];
+if (!is_leap_year) {
+# Since all normals data includes Feb 29, we have to
+# subtract 1 from DOY if we're not in a leap year since
+# we removed the Feb 29 row from the data frame above.
+doy = as.integer(doy) - 1;
+}
+tmin = norm_data_frame[index,"TMIN"][1];
+tmax = norm_data_frame[index,"TMAX"][1];
+return(list(latitude, longitude, date_str, doy, tmin, tmax));
 }
 get_temperature_at_hour = function(latitude, temperature_data_frame, row, num_days) {
 # Base development threshold for Brown Marmorated Stink Bug
 # insect phenology model.
 averages = sum(dh) / 24;
 }
 return(c(curr_mean_temp, averages))
 }
+get_tick_index = function(index, last_tick, ticks, month_labels) {
+# The R code tries hard not to draw overlapping tick labels, and so
+# will omit labels where they would abut or overlap previously drawn
+# labels. This can result in, for example, every other tick being
+# labelled.  We'll keep track of the last tick to make sure all of
+# the month labels are displayed, and missing ticks are restricted
+# to Sundays which have no labels anyway.
+if (last_tick==0) {
+return(length(ticks)+1);
+}
+last_saved_tick = ticks[[length(ticks)]];
+if (index-last_saved_tick<3) {
+last_saved_month = month_labels[[length(month_labels)]];
+if (last_saved_month=="") {
+# We're safe overwriting a tick
+# with no label (i.e., a Sunday tick).
+return(length(ticks));
+} else {
+# Don't eliminate a Month label.
+return(NULL);
+}
+}
+return(length(ticks)+1);
+}
+get_total_days = function(is_leap_year) {
+# Get the total number of days in the current year.
+if (is_leap_year) {
+return(366);
+} else {
+return(365);
+}
+}
+get_x_axis_ticks_and_labels = function(temperature_data_frame, num_rows, start_doy_ytd, end_doy_ytd) {
+# Keep track of the years to see if spanning years.
+month_labels = list();
+ticks = list();
+current_month_label = NULL;
+last_tick = 0;
+for (i in 1:num_rows) {
+if (start_doy_ytd > 1 & i==start_doy_ytd-1) {
+# Add a tick for the end of the 30 year normnals data
+# that was prepended to the year-to-date data.
+tick_index = get_tick_index(i, last_tick, ticks, month_labels)
+ticks[tick_index] = i;
+month_labels[tick_index] = "End prepended 30 year normals";
+last_tick = i;
+} else  if (end_doy_ytd > 0 & i==end_doy_ytd+1) {
+# Add a tick for the start of the 30 year normnals data
+# that was appended to the year-to-date data.
+tick_index = get_tick_index(i, last_tick, ticks, month_labels)
+ticks[tick_index] = i;
+month_labels[tick_index] = "Start appended 30 year normals";
+last_tick = i;
+} else if (i==num_rows) {
+# Add a tick for the last day of the year.
+tick_index = get_tick_index(i, last_tick, ticks, month_labels)
+ticks[tick_index] = i;
+month_labels[tick_index] = "";
+last_tick = i;
+} else {
+# Get the year and month from the date which
+# has the format YYYY-MM-DD.
+date = format(temperature_data_frame$DATE[i]);
+# Get the month label.
+items = strsplit(date, "-")[[1]];
+month = items[2];
+month_label = month.abb[as.integer(month)];
+if (!identical(current_month_label, month_label)) {
+# Add an x-axis tick for the month.
+tick_index = get_tick_index(i, last_tick, ticks, month_labels)
+ticks[tick_index] = i;
+month_labels[tick_index] = month_label;
+current_month_label = month_label;
+last_tick = i;
+}
+tick_index = get_tick_index(i, last_tick, ticks, month_labels)
+if (!is.null(tick_index)) {
+# Get the day.
+day = weekdays(as.Date(date));
+if (day=="Sunday") {
+# Add an x-axis tick if we're on a Sunday.
+ticks[tick_index] = i;
+# Add a blank month label so it is not displayed.
+month_labels[tick_index] = "";
+last_tick = i;
+}
+}
+}
+}
+return(list(ticks, month_labels));
+}
+is_leap_year = function(date_str) {
+# Extract the year from the date_str.
+date = format(date_str);
+items = strsplit(date, "-")[[1]];
+year = as.integer(items[1]);
+if (((year %% 4 == 0) & (year %% 100 != 0)) | (year %% 400 == 0)) {
+return(TRUE);
+} else {
+return(FALSE);
+}
+}
 mortality.adult = function(temperature) {
 if (temperature < 12.7) {
 mortality.probability = 0.002;
 }
 else {
 mortality.probability = temperature * 0.0008 + 0.03;
 }
 return(mortality.probability);
 }
-parse_input_data = function(input_file, num_rows) {
+parse_input_data = function(input_ytd, input_norm, num_days_ytd, location) {
-# Read in the input temperature datafile into a data frame.
+if (is.null(input_ytd)) {
-temperature_data_frame = read.csv(file=input_file, header=T, strip.white=TRUE, sep=",");
+# We're analysing only the 30 year normals data, so create an empty
-num_columns = dim(temperature_data_frame)[2];
+# data frame for containing temperature data after it is converted
-if (num_columns == 6) {
+# from the 30 year normals format to the year-to-date format.
-# The input data has the following 6 columns:
+temperature_data_frame = data.frame(matrix(ncol=6, nrow=0));
+colnames(temperature_data_frame) = c("LATITUDE", "LONGITUDE", "DATE", "DOY", "TMIN", "TMAX");
+# Base all dates on the current date since 30 year
+# normals data does not include any dates.
+year = format(Sys.Date(), "%Y");
+start_date = paste(year, "01", "01", sep="-");
+end_date = paste(year, "12", "31", sep="-");
+# Set invalid start and end DOY.
+start_doy_ytd = 0;
+end_doy_ytd = 0;
+} else {
+# Read the input_ytd temperature datafile into a data frame.
+# The input_ytd data has the following 6 columns:
 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
-# Set the column names for access when adding daylight length..
+temperature_data_frame = read.csv(file=input_ytd, header=T, strip.white=TRUE, stringsAsFactors=FALSE, sep=",");
-colnames(temperature_data_frame) = c("LATITUDE","LONGITUDE", "DATE", "DOY", "TMIN", "TMAX");
+# Set the temperature_data_frame column names for access.
-# Add a column containing the daylight length for each day.
+colnames(temperature_data_frame) = c("LATITUDE", "LONGITUDE", "DATE", "DOY", "TMIN", "TMAX");
-temperature_data_frame = add_daylight_length(temperature_data_frame, num_columns, num_rows);
+# Get the start date.
-# Reset the column names with the additional column for later access.
+start_date = temperature_data_frame$DATE[1];
-colnames(temperature_data_frame) = c("LATITUDE","LONGITUDE", "DATE", "DOY", "TMIN", "TMAX", "DAYLEN");
+end_date = temperature_data_frame$DATE[num_days_ytd];
-}
+# Extract the year from the start date.
-return(temperature_data_frame);
+date_str = format(start_date);
-}
+date_str_items = strsplit(date_str, "-")[[1]];
+year = date_str_items[1];
+# Save the first DOY to later check if start_date is Jan 1.
-render_chart = function(chart_type, insect, location, latitude, start_date, end_date, days, maxval, plot_std_error,
+start_doy_ytd = as.integer(temperature_data_frame$DOY[1]);
-group1, group2, group3, group1_std_error, group2_std_error, group3_std_error, date_labels) {
+end_doy_ytd = as.integer(temperature_data_frame$DOY[num_days_ytd]);
-if (chart_type == "pop_size_by_life_stage") {
+}
-title = paste(insect, ": Total pop. by life stage :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ");
+# See if we're in a leap year.
-legend_text = c("Egg", "Nymph", "Adult");
+is_leap_year = is_leap_year(start_date);
-columns = c(4, 2, 1);
+# Get the number of days in the year.
-} else if (chart_type == "pop_size_by_generation") {
+total_days = get_total_days(is_leap_year);
-title = paste(insect, ": Total pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ");
+# Read the input_norm temperature datafile into a data frame.
+# The input_norm data has the following 10 columns:
+# STATIONID, LATITUDE, LONGITUDE, ELEV_M, NAME, ST, MMDD, DOY, TMIN, TMAX
+norm_data_frame = read.csv(file=input_norm, header=T, strip.white=TRUE, stringsAsFactors=FALSE, sep=",");
+# Set the norm_data_frame column names for access.
+colnames(norm_data_frame) = c("STATIONID", "LATITUDE","LONGITUDE", "ELEV_M", "NAME", "ST", "MMDD", "DOY", "TMIN", "TMAX");
+# All normals data includes Feb 29 which is row 60 in
+# the data, so delete that row if we're not in a leap year.
+if (!is_leap_year) {
+norm_data_frame = norm_data_frame[-c(60),];
+}
+# Set the location to be the station name if the user elected no to enter it.
+if (is.null(location) | length(location)==0) {
+location = norm_data_frame$NAME[1];
+}
+if (is.null(input_ytd)) {
+# Convert the 30 year normals data to the year-to-date format.
+for (i in 1:total_days) {
+temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i);
+}
+} else {
+# Merge the year-to-date data with the 30 year normals data.
+if (start_doy_ytd > 1) {
+# The year-to-date data starts after Jan 1, so create a
+# temporary data frame to contain the 30 year normals data
+# from Jan 1 to the date immediately prior to start_date.
+tmp_data_frame = temperature_data_frame[FALSE,];
+for (i in 1:start_doy_ytd-1) {
+tmp_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i);
+}
+# Next merge the temporary data frame with the year-to-date data frame.
+temperature_data_frame = rbind(tmp_data_frame, temperature_data_frame);
+}
+# Define the next row for the year-to-date data from the 30 year normals data.
+first_normals_append_row = end_doy_ytd + 1;
+# Append the 30 year normals data to the year-to-date data.
+for (i in first_normals_append_row:total_days) {
+temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i);
+}
+}
+# Add a column containing the daylight length for each day.
+temperature_data_frame = add_daylight_length(temperature_data_frame, total_days);
+return(list(temperature_data_frame, start_date, end_date, start_doy_ytd, end_doy_ytd, is_leap_year, total_days, location));
+}
+render_chart = function(ticks, date_labels, chart_type, plot_std_error, insect, location, latitude, start_date, end_date, days, maxval,
+replications, life_stage, group, group_std_error, group2=NULL, group2_std_error=NULL, group3=NULL, group3_std_error=NULL,
+life_stages_adult=NULL, life_stages_nymph=NULL) {
+if (chart_type=="pop_size_by_life_stage") {
+if (life_stage=="Total") {
+title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
+legend_text = c("Egg", "Nymph", "Adult");
+columns = c(4, 2, 1);
+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);
+legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3);
+lines(days, group2, lwd=2, lty=1, col=2);
+lines(days, group3, lwd=2, lty=1, col=4);
+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);
+axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
+if (plot_std_error=="yes") {
+# Standard error for group.
+lines(days, group+group_std_error, lty=2);
+lines(days, group-group_std_error, lty=2);
+# Standard error for group2.
+lines(days, group2+group2_std_error, col=2, lty=2);
+lines(days, group2-group2_std_error, col=2, lty=2);
+# Standard error for group3.
+lines(days, group3+group3_std_error, col=4, lty=2);
+lines(days, group3-group3_std_error, col=4, lty=2);
+}
+} else {
+if (life_stage=="Egg") {
+title = paste(insect,  ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
+legend_text = c(life_stage);
+columns = c(4);
+} else if (life_stage=="Nymph") {
+stage = paste(life_stages_nymph, "Nymph Pop :", sep=" ");
+title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
+legend_text = c(paste(life_stages_nymph, life_stage, sep=" "));
+columns = c(2);
+} else if (life_stage=="Adult") {
+stage = paste(life_stages_adult, "Adult Pop", sep=" ");
+title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
+legend_text = c(paste(life_stages_adult, life_stage, sep=" "));
+columns = c(1);
+}
+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);
+legend("topleft", legend_text, lty=c(1), col="black", cex=3);
+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);
+axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
+if (plot_std_error=="yes") {
+# Standard error for group.
+lines(days, group+group_std_error, lty=2);
+lines(days, group-group_std_error, lty=2);
+}
+}
+} else if (chart_type=="pop_size_by_generation") {
+if (life_stage=="Total") {
+title_str = ": Total Pop by Gen :";
+} else if (life_stage=="Egg") {
+title_str = ": Egg Pop by Gen :";
+} else if (life_stage=="Nymph") {
+title_str = paste(":", life_stages_nymph, "Nymph Pop by Gen", ":", sep=" ");
+} else if (life_stage=="Adult") {
+title_str = paste(":", life_stages_adult, "Adult Pop by Gen", ":", sep=" ");
+}
+title = paste(insect, ": Reps", replications, title_str, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
 legend_text = c("P", "F1", "F2");
 columns = c(1, 2, 4);
-} else if (chart_type == "adult_pop_size_by_generation") {
+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);
-title = paste(insect, ": Adult pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ");
+legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3);
-legend_text = c("P", "F1", "F2");
+lines(days, group2, lwd=2, lty=1, col=2);
-columns = c(1, 2, 4);
+lines(days, group3, lwd=2, lty=1, col=4);
-}
+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);
-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);
+axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
-legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3);
+if (plot_std_error=="yes") {
-lines(days, group2, lwd=2, lty=1, col=2);
+# Standard error for group.
-lines(days, group3, lwd=2, lty=1, col=4);
+lines(days, group+group_std_error, lty=2);
-axis(1, at=c(1:length(date_labels)) * 30 - 15, cex.axis=3, labels=date_labels);
+lines(days, group-group_std_error, lty=2);
-axis(2, cex.axis=3);
+# Standard error for group2.
-if (plot_std_error==1) {
+lines(days, group2+group2_std_error, col=2, lty=2);
-# Standard error for group1.
+lines(days, group2-group2_std_error, col=2, lty=2);
-lines(days, group1+group1_std_error, lty=2);
+# Standard error for group3.
-lines(days, group1-group1_std_error, lty=2);
+lines(days, group3+group3_std_error, col=4, lty=2);
-# Standard error for group2.
+lines(days, group3-group3_std_error, col=4, lty=2);
-lines(days, group2+group2_std_error, col=2, lty=2);
+}
-lines(days, group2-group2_std_error, col=2, lty=2);
+}
-# Standard error for group3.
+}
-lines(days, group3+group3_std_error, col=4, lty=2);
-lines(days, group3-group3_std_error, col=4, lty=2);
+# Determine if we're plotting generations separately.
-}
+if (opt$plot_generations_separately=="yes") {
-}
+plot_generations_separately = TRUE;
+} else {
-temperature_data_frame = parse_input_data(opt$input, opt$num_days);
+plot_generations_separately = FALSE;
-# All latitude values are the same, so get the value from the first row.
+}
+# Display the total number of days in the Galaxy history item blurb.
+cat("Year-to-date number of days: ", opt$num_days_ytd, "\n");
+# Parse the inputs.
+data_list = parse_input_data(opt$input_ytd, opt$input_norm, opt$num_days_ytd, opt$location);
+temperature_data_frame = data_list[[1]];
+# Information needed for plots.
+start_date = data_list[[2]];
+end_date = data_list[[3]];
+start_doy_ytd = data_list[[4]];
+end_doy_ytd = data_list[[5]];
+is_leap_year = data_list[[6]];
+total_days = data_list[[7]];
+total_days_vector = c(1:total_days);
+location =  data_list[[8]];
+# Create copies of the temperature data for generations P, F1 and F2 if we're plotting generations separately.
+if (plot_generations_separately) {
+temperature_data_frame_P = data.frame(temperature_data_frame);
+temperature_data_frame_F1 = data.frame(temperature_data_frame);
+temperature_data_frame_F2 = data.frame(temperature_data_frame);
+}
+# Get the ticks date labels for plots.
+ticks_and_labels = get_x_axis_ticks_and_labels(temperature_data_frame, total_days, start_doy_ytd, end_doy_ytd);
+ticks = c(unlist(ticks_and_labels[1]));
+date_labels = c(unlist(ticks_and_labels[2]));
+# All latitude values are the same, so get the value for plots from the first row.
 latitude = temperature_data_frame$LATITUDE[1];
-num_columns = dim(temperature_data_frame)[2];
-date_labels = get_date_labels(temperature_data_frame, opt$num_days);
+# Determine the specified life stages for processing.
+# Split life_stages into a list of strings for plots.
+life_stages_str = as.character(opt$life_stages);
+life_stages = strsplit(life_stages_str, ",")[[1]];
+# Determine the data we need to generate for plotting.
+process_eggs = FALSE;
+process_nymphs = FALSE;
+process_young_nymphs = FALSE;
+process_old_nymphs = FALSE;
+process_total_nymphs = FALSE;
+process_adults = FALSE;
+process_previttelogenic_adults = FALSE;
+process_vittelogenic_adults = FALSE;
+process_diapausing_adults = FALSE;
+process_total_adults = FALSE;
+for (life_stage in life_stages) {
+if (life_stage=="Total") {
+process_eggs = TRUE;
+process_nymphs = TRUE;
+process_adults = TRUE;
+} else if (life_stage=="Egg") {
+process_eggs = TRUE;
+} else if (life_stage=="Nymph") {
+process_nymphs = TRUE;
+} else if (life_stage=="Adult") {
+process_adults = TRUE;
+}
+}
+if (process_nymphs) {
+# Split life_stages_nymph into a list of strings for plots.
+life_stages_nymph_str = as.character(opt$life_stages_nymph);
+life_stages_nymph = strsplit(life_stages_nymph_str, ",")[[1]];
+for (life_stage_nymph in life_stages_nymph) {
+if (life_stage_nymph=="Young") {
+process_young_nymphs = TRUE;
+} else if (life_stage_nymph=="Old") {
+process_old_nymphs = TRUE;
+} else if (life_stage_nymph=="Total") {
+process_total_nymphs = TRUE;
+}
+}
+}
+if (process_adults) {
+# Split life_stages_adult into a list of strings for plots.
+life_stages_adult_str = as.character(opt$life_stages_adult);
+life_stages_adult = strsplit(life_stages_adult_str, ",")[[1]];
+for (life_stage_adult in life_stages_adult) {
+if (life_stage_adult=="Pre-vittelogenic") {
+process_previttelogenic_adults = TRUE;
+} else if (life_stage_adult=="Vittelogenic") {
+process_vittelogenic_adults = TRUE;
+} else if (life_stage_adult=="Diapausing") {
+process_diapausing_adults = TRUE;
+} else if (life_stage_adult=="Total") {
+process_total_adults = TRUE;
+}
+}
+}
 # Initialize matrices.
-Eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+if (process_eggs) {
-YoungNymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+Eggs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
-OldNymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+}
-Previtellogenic.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+if (process_young_nymphs | process_total_nymphs) {
-Vitellogenic.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+YoungNymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
-Diapausing.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+}
+if (process_old_nymphs | process_total_nymphs) {
-newborn.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+OldNymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
-adult.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+}
-death.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+if (process_previttelogenic_adults | process_total_adults) {
+Previttelogenic.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
-P.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+}
-P_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+if (process_vittelogenic_adults | process_total_adults) {
-F1.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+Vittelogenic.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
-F1_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+}
-F2.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+if (process_diapausing_adults | process_total_adults) {
-F2_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+Diapausing.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+}
-population.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+newborn.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+adult.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+death.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+if (plot_generations_separately) {
+# P is Parental, or overwintered adults.
+P.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+# F1 is the first field-produced generation.
+F1.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+# F2 is the second field-produced generation.
+F2.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+if (process_eggs) {
+P_eggs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+F1_eggs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+F2_eggs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+}
+if (process_young_nymphs) {
+P_young_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+F1_young_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+F2_young_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+}
+if (process_old_nymphs) {
+P_old_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+F1_old_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+F2_old_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+}
+if (process_total_nymphs) {
+P_total_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+F1_total_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+F2_total_nymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+}
+if (process_previttelogenic_adults) {
+P_previttelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+F1_previttelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+F2_previttelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+}
+if (process_vittelogenic_adults) {
+P_vittelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+F1_vittelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+F2_vittelogenic_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+}
+if (process_diapausing_adults) {
+P_diapausing_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+F1_diapausing_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+F2_diapausing_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+}
+if (process_total_adults) {
+P_total_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+F1_total_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+F2_total_adults.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
+}
+}
+# Total population.
+population.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);
 # Process replications.
-for (N.replications in 1:opt$replications) {
+for (current_replication in 1:opt$replications) {
 # Start with the user-defined number of insects per replication.
 num_insects = opt$insects_per_replication;
 # Generation, Stage, degree-days, T, Diapause.
 vector.ini = c(0, 3, 0, 0, 0);
-# Overwintering, previttelogenic, degree-days=0, T=0, no-diapause.
+# Replicate to create a matrix where the columns are
+# Generation, Stage, degree-days, T, Diapause and the
+# rows are the initial number of insects per replication.
 vector.matrix = rep(vector.ini, num_insects);
-# Complete matrix for the population.
+# Complete transposed matrix for the population, so now
+# the rows are Generation, Stage, degree-days, T, Diapause
 vector.matrix = base::t(matrix(vector.matrix, nrow=5));
 # Time series of population size.
-Eggs = rep(0, opt$num_days);
+if (process_eggs) {
-YoungNymphs = rep(0, opt$num_days);
+Eggs = rep(0, total_days);
-OldNymphs = rep(0, opt$num_days);
+}
-Previtellogenic = rep(0, opt$num_days);
+if (process_young_nymphs | process_total_nymphs) {
-Vitellogenic = rep(0, opt$num_days);
+YoungNymphs = rep(0, total_days);
-Diapausing = rep(0, opt$num_days);
+}
+if (process_old_nymphs | process_total_nymphs) {
-N.newborn = rep(0, opt$num_days);
+OldNymphs = rep(0, total_days);
-N.adult = rep(0, opt$num_days);
+}
-N.death = rep(0, opt$num_days);
+if (process_previttelogenic_adults | process_total_adults) {
+Previttelogenic = rep(0, total_days);
-overwintering_adult.population = rep(0, opt$num_days);
+}
-first_generation.population = rep(0, opt$num_days);
+if (process_vittelogenic_adults | process_total_adults) {
-second_generation.population = rep(0, opt$num_days);
+Vittelogenic = rep(0, total_days);
+}
-P.adult = rep(0, opt$num_days);
+if (process_diapausing_adults | process_total_adults) {
-F1.adult = rep(0, opt$num_days);
+Diapausing = rep(0, total_days);
-F2.adult = rep(0, opt$num_days);
+}
+N.newborn = rep(0, total_days);
+N.adult = rep(0, total_days);
+N.death = rep(0, total_days);
+overwintering_adult.population = rep(0, total_days);
+first_generation.population = rep(0, total_days);
+second_generation.population = rep(0, total_days);
+if (plot_generations_separately) {
+# P is Parental, or overwintered adults.
+# F1 is the first field-produced generation.
+# F2 is the second field-produced generation.
+if (process_eggs) {
+P.egg = rep(0, total_days);
+F1.egg = rep(0, total_days);
+F2.egg = rep(0, total_days);
+}
+if (process_young_nymphs) {
+P.young_nymph = rep(0, total_days);
+F1.young_nymph = rep(0, total_days);
+F2.young_nymph = rep(0, total_days);
+}
+if (process_old_nymphs) {
+P.old_nymph = rep(0, total_days);
+F1.old_nymph = rep(0, total_days);
+F2.old_nymph = rep(0, total_days);
+}
+if (process_total_nymphs) {
+P.total_nymph = rep(0, total_days);
+F1.total_nymph = rep(0, total_days);
+F2.total_nymph = rep(0, total_days);
+}
+if (process_previttelogenic_adults) {
+P.previttelogenic_adult = rep(0, total_days);
+F1.previttelogenic_adult = rep(0, total_days);
+F2.previttelogenic_adult = rep(0, total_days);
+}
+if (process_vittelogenic_adults) {
+P.vittelogenic_adult = rep(0, total_days);
+F1.vittelogenic_adult = rep(0, total_days);
+F2.vittelogenic_adult = rep(0, total_days);
+}
+if (process_diapausing_adults) {
+P.diapausing_adult = rep(0, total_days);
+F1.diapausing_adult = rep(0, total_days);
+F2.diapausing_adult = rep(0, total_days);
+}
+if (process_total_adults) {
+P.total_adult = rep(0, total_days);
+F1.total_adult = rep(0, total_days);
+F2.total_adult = rep(0, total_days);
+}
+}
 total.population = NULL;
+averages.day = rep(0, total_days);
-averages.day = rep(0, opt$num_days);
+# All the days included in the input_ytd temperature dataset.
-# All the days included in the input temperature dataset.
+for (row in 1:total_days) {
-for (row in 1:opt$num_days) {
 # Get the integer day of the year for the current row.
 doy = temperature_data_frame$DOY[row];
 # Photoperiod in the day.
 photoperiod = temperature_data_frame$DAYLEN[row];
-temp.profile = get_temperature_at_hour(latitude, temperature_data_frame, row, opt$num_days);
+temp.profile = get_temperature_at_hour(latitude, temperature_data_frame, row, total_days);
 mean.temp = temp.profile[1];
 averages.temp = temp.profile[2];
 averages.day[row] = averages.temp;
 # Trash bin for death.
 death.vector = NULL;
 if (vector.individual[2] == 0) {
 # Egg.
 death.probability = opt$egg_mortality * mortality.egg(mean.temp);
 }
 else if (vector.individual[2] == 1 | vector.individual[2] == 2) {
+# Nymph.
 death.probability = opt$nymph_mortality * mortality.nymph(mean.temp);
 }
 else if (vector.individual[2] == 3 | vector.individual[2] == 4 | vector.individual[2] == 5) {
 # Adult.
 if (doy < day.kill) {
 death.vector = c(death.vector, i);
 }
 else {
 # End of diapause.
 if (vector.individual[1] == 0 && vector.individual[2] == 3) {
-# Overwintering adult (previttelogenic).
+# Overwintering adult (pre-vittelogenic).
 if (photoperiod > opt$photoperiod && vector.individual[3] > 68 && doy < 180) {
 # Add 68C to become fully reproductively matured.
 # Transfer to vittelogenic.
 vector.individual = c(0, 4, 0, 0, 0);
 vector.matrix[i,] = vector.individual;
 }
 else {
-# Add to # Add average temperature for current day.
+# Add average temperature for current day.
 vector.individual[3] = vector.individual[3] + averages.temp;
 # Add 1 day in current stage.
 vector.individual[4] = vector.individual[4] + 1;
 vector.matrix[i,] = vector.individual;
 }
 }
 if (vector.individual[1] != 0 && vector.individual[2] == 3) {
-# Not overwintering adult (previttelogenic).
+# Not overwintering adult (pre-vittelogenic).
 current.gen = vector.individual[1];
 if (vector.individual[3] > 68) {
 # Add 68C to become fully reproductively matured.
 # Transfer to vittelogenic.
 vector.individual = c(current.gen, 4, 0, 0, 0);
 # Add 1 day in current stage.
 vector.individual[4] = vector.individual[4] + 1;
 }
 vector.matrix[i,] = vector.individual;
 }
-# Old nymph to adult: previttelogenic or diapausing?
+# Old nymph to adult: pre-vittelogenic or diapausing?
 if (vector.individual[2] == 2) {
 # Add average temperature for current day.
 vector.individual[3] = vector.individual[3] + averages.temp;
 if (vector.individual[3] >= (200+opt$adult_accumulation)) {
 # From old to adult, degree_days requirement met.
 num_insects.newborn = length(birth.vector[,1]);
 vector.matrix = rbind(vector.matrix, birth.vector);
 # Update population size for the next day.
 num_insects = num_insects - num_insects.death + num_insects.newborn;
-# Aggregate results by day.
+# Aggregate results by day.  Due to multiple transpose calls
-# Egg population size.
+# on vector.matrix above, the columns of vector.matrix
-Eggs[row] = sum(vector.matrix[,2]==0);
+# are now Generation, Stage, degree-days, T, Diapause,
-# Young nymph population size.
+if (process_eggs) {
-YoungNymphs[row] = sum(vector.matrix[,2]==1);
+# For egg population size, column 2 (Stage), must be 0.
-# Old nymph population size.
+Eggs[row] = sum(vector.matrix[,2]==0);
-OldNymphs[row] = sum(vector.matrix[,2]==2);
+}
-# Previtellogenic population size.
+if (process_young_nymphs | process_total_nymphs) {
-Previtellogenic[row] = sum(vector.matrix[,2]==3);
+# For young nymph population size, column 2 (Stage) must be 1.
-# Vitellogenic population size.
+YoungNymphs[row] = sum(vector.matrix[,2]==1);
-Vitellogenic[row] = sum(vector.matrix[,2]==4);
+}
-# Diapausing population size.
+if (process_old_nymphs | process_total_nymphs) {
-Diapausing[row] = sum(vector.matrix[,2]==5);
+# For old nymph population size, column 2 (Stage) must be 2.
+OldNymphs[row] = sum(vector.matrix[,2]==2);
+}
+if (process_previttelogenic_adults | process_total_adults) {
+# For pre-vittelogenic population size, column 2 (Stage) must be 3.
+Previttelogenic[row] = sum(vector.matrix[,2]==3);
+}
+if (process_vittelogenic_adults | process_total_adults) {
+# For vittelogenic population size, column 2 (Stage) must be 4.
+Vittelogenic[row] = sum(vector.matrix[,2]==4);
+}
+if (process_diapausing_adults | process_total_adults) {
+# For diapausing population size, column 2 (Stage) must be 5.
+Diapausing[row] = sum(vector.matrix[,2]==5);
+}
 # Newborn population size.
 N.newborn[row] = num_insects.newborn;
 # Adult population size.
 N.adult[row] = sum(vector.matrix[,2]==3) + sum(vector.matrix[,2]==4) + sum(vector.matrix[,2]==5);
 # Dead population size.
 N.death[row] = num_insects.death;
 total.population = c(total.population, num_insects);
-# Overwintering adult population size.
+# For overwintering adult (P) population
+# size, column 1 (Generation) must be 0.
 overwintering_adult.population[row] = sum(vector.matrix[,1]==0);
-# First generation population size.
+# For first field generation (F1) population
+# size, column 1 (Generation) must be 1.
 first_generation.population[row] = sum(vector.matrix[,1]==1);
-# Second generation population size.
+# For second field generation (F2) population
+# size, column 1 (Generation) must be 2.
 second_generation.population[row] = sum(vector.matrix[,1]==2);
-# P adult population size.
+if (plot_generations_separately) {
-P.adult[row] = sum(vector.matrix[,1]==0);
+if (process_eggs) {
-# F1 adult population size.
+# For egg life stage of generation P population size,
-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));
+# column 1 (generation) is 0 and column 2 (Stage) is 0.
-# F2 adult population size
+P.egg[row] = sum(vector.matrix[,1]==0 & vector.matrix[,2]==0);
-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));
+# For egg life stage of generation F1 population size,
-}   # End of days specified in the input temperature data.
+# column 1 (generation) is 1 and column 2 (Stage) is 0.
+F1.egg[row] = sum(vector.matrix[,1]==1 & vector.matrix[,2]==0);
+# For egg life stage of generation F2 population size,
+# column 1 (generation) is 2 and column 2 (Stage) is 0.
+F2.egg[row] = sum(vector.matrix[,1]==2 & vector.matrix[,2]==0);
+}
+if (process_young_nymphs) {
+# For young nymph life stage of generation P population
+# size, the following combination is required:
+# - column 1 (Generation) is 0 and column 2 (Stage) is 1 (Young nymph)
+P.young_nymph[row] = sum(vector.matrix[,1]==0 & vector.matrix[,2]==1);
+# For young nymph life stage of generation F1 population
+# size, the following combination is required:
+# - column 1 (Generation) is 1 and column 2 (Stage) is 1 (Young nymph)
+F1.young_nymph[row] = sum(vector.matrix[,1]==1 & vector.matrix[,2]==1);
+# For young nymph life stage of generation F2 population
+# size, the following combination is required:
+# - column 1 (Generation) is 2 and column 2 (Stage) is 1 (Young nymph)
+F2.young_nymph[row] = sum(vector.matrix[,1]==2 & vector.matrix[,2]==1);
+}
+if (process_old_nymphs) {
+# For old nymph life stage of generation P population
+# size, the following combination is required:
+# - column 1 (Generation) is 0 and column 2 (Stage) is 2 (Old nymph)
+P.old_nymph[row] = sum(vector.matrix[,1]==0 & vector.matrix[,2]==2);
+# For old nymph life stage of generation F1 population
+# size, the following combination is required:
+# - column 1 (Generation) is 1 and column 2 (Stage) is 2 (Old nymph)
+F1.old_nymph[row] = sum(vector.matrix[,1]==1 & vector.matrix[,2]==2);
+# For old nymph life stage of generation F2 population
+# size, the following combination is required:
+# - column 1 (Generation) is 2 and column 2 (Stage) is 2 (Old nymph)
+F2.old_nymph[row] = sum(vector.matrix[,1]==2 & vector.matrix[,2]==2);
+}
+if (process_total_nymphs) {
+# For total nymph life stage of generation P population
+# size, one of the following combinations is required:
+# - column 1 (Generation) is 0 and column 2 (Stage) is 1 (Young nymph)
+# - column 1 (Generation) is 0 and column 2 (Stage) is 2 (Old nymph)
+P.total_nymph[row] = sum((vector.matrix[,1]==0 & vector.matrix[,2]==1) | (vector.matrix[,1]==0 & vector.matrix[,2]==2));
+# For total nymph life stage of generation F1 population
+# size, one of the following combinations is required:
+# - column 1 (Generation) is 1 and column 2 (Stage) is 1 (Young nymph)
+# - column 1 (Generation) is 1 and column 2 (Stage) is 2 (Old nymph)
+F1.total_nymph[row] = sum((vector.matrix[,1]==1 & vector.matrix[,2]==1) | (vector.matrix[,1]==1 & vector.matrix[,2]==2));
+# For total nymph life stage of generation F2 population
+# size, one of the following combinations is required:
+# - column 1 (Generation) is 2 and column 2 (Stage) is 1 (Young nymph)
+# - column 1 (Generation) is 2 and column 2 (Stage) is 2 (Old nymph)
+F2.total_nymph[row] = sum((vector.matrix[,1]==2 & vector.matrix[,2]==1) | (vector.matrix[,1]==2 & vector.matrix[,2]==2));
+}
+if (process_previttelogenic_adults) {
+# For previttelogenic adult life stage of generation P population
+# size, the following combination is required:
+# - column 1 (Generation) is 0 and column 2 (Stage) is 3 (Pre-vittelogenic)
+P.previttelogenic_adult[row] = sum(vector.matrix[,1]==0 & vector.matrix[,2]==3);
+# For previttelogenic adult life stage of generation F1 population
+# size, the following combination is required:
+# - column 1 (Generation) is 1 and column 2 (Stage) is 3 (Pre-vittelogenic)
+F1.previttelogenic_adult[row] = sum(vector.matrix[,1]==1 & vector.matrix[,2]==3);
+# For previttelogenic adult life stage of generation F2 population
+# size, the following combination is required:
+# - column 1 (Generation) is 2 and column 2 (Stage) is 3 (Pre-vittelogenic)
+F2.previttelogenic_adult[row] = sum(vector.matrix[,1]==2 & vector.matrix[,2]==3);
+}
+if (process_vittelogenic_adults) {
+# For vittelogenic adult life stage of generation P population
+# size, the following combination is required:
+# - column 1 (Generation) is 0 and column 2 (Stage) is 4 (Vittelogenic)
+P.vittelogenic_adult[row] = sum(vector.matrix[,1]==0 & vector.matrix[,2]==4);
+# For vittelogenic adult life stage of generation F1 population
+# size, the following combination is required:
+# - column 1 (Generation) is 1 and column 2 (Stage) is 4 (Vittelogenic)
+F1.vittelogenic_adult[row] = sum(vector.matrix[,1]==1 & vector.matrix[,2]==4);
+# For vittelogenic adult life stage of generation F2 population
+# size, the following combination is required:
+# - column 1 (Generation) is 2 and column 2 (Stage) is 4 (Vittelogenic)
+F2.vittelogenic_adult[row] = sum(vector.matrix[,1]==2 & vector.matrix[,2]==4);
+}
+if (process_diapausing_adults) {
+# For diapausing adult life stage of generation P population
+# size, the following combination is required:
+# - column 1 (Generation) is 0 and column 2 (Stage) is 5 (Diapausing)
+P.diapausing_adult[row] = sum(vector.matrix[,1]==0 & vector.matrix[,2]==5);
+# For diapausing adult life stage of generation F1 population
+# size, the following combination is required:
+# - column 1 (Generation) is 1 and column 2 (Stage) is 5 (Diapausing)
+F1.diapausing_adult[row] = sum(vector.matrix[,1]==1 & vector.matrix[,2]==5);
+# For diapausing adult life stage of generation F2 population
+# size, the following combination is required:
+# - column 1 (Generation) is 2 and column 2 (Stage) is 5 (Diapausing)
+F2.diapausing_adult[row] = sum(vector.matrix[,1]==2 & vector.matrix[,2]==5);
+}
+if (process_total_adults) {
+# For total adult life stage of generation P population
+# size, one of the following combinations is required:
+# - column 1 (Generation) is 0 and column 2 (Stage) is 3 (Pre-vittelogenic)
+# - column 1 (Generation) is 0 and column 2 (Stage) is 4 (Vittelogenic)
+# - column 1 (Generation) is 0 and column 2 (Stage) is 5 (Diapausing)
+P.total_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));
+# For total adult life stage of generation F1 population
+# size, one of the following combinations is required:
+# - column 1 (Generation) is 1 and column 2 (Stage) is 3 (Pre-vittelogenic)
+# - column 1 (Generation) is 1 and column 2 (Stage) is 4 (Vittelogenic)
+# - column 1 (Generation) is 1 and column 2 (Stage) is 5 (Diapausing)
+F1.total_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));
+# For total adult life stage of generation F2 population
+# size, one of the following combinations is required:
+# - column 1 (Generation) is 2 and column 2 (Stage) is 3 (Pre-vittelogenic)
+# - column 1 (Generation) is 2 and column 2 (Stage) is 4 (Vittelogenic)
+# - column 1 (Generation) is 2 and column 2 (Stage) is 5 (Diapausing)
+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));
+}
+}
+}   # End of days specified in the input_ytd temperature data.
 averages.cum = cumsum(averages.day);
 # Define the output values.
-Eggs.replications[,N.replications] = Eggs;
+if (process_eggs) {
-YoungNymphs.replications[,N.replications] = YoungNymphs;
+Eggs.replications[,current_replication] = Eggs;
-OldNymphs.replications[,N.replications] = OldNymphs;
+}
-Previtellogenic.replications[,N.replications] = Previtellogenic;
+if (process_young_nymphs | process_total_nymphs) {
-Vitellogenic.replications[,N.replications] = Vitellogenic;
+YoungNymphs.replications[,current_replication] = YoungNymphs;
-Diapausing.replications[,N.replications] = Diapausing;
+}
+if (process_old_nymphs | process_total_nymphs) {
-newborn.replications[,N.replications] = N.newborn;
+OldNymphs.replications[,current_replication] = OldNymphs;
-adult.replications[,N.replications] = N.adult;
+}
-death.replications[,N.replications] = N.death;
+if (process_previttelogenic_adults | process_total_adults) {
+Previttelogenic.replications[,current_replication] = Previttelogenic;
-P.replications[,N.replications] = overwintering_adult.population;
+}
-P_adults.replications[,N.replications] = P.adult;
+if (process_vittelogenic_adults | process_total_adults) {
-F1.replications[,N.replications] = first_generation.population;
+Vittelogenic.replications[,current_replication] = Vittelogenic;
-F1_adults.replications[,N.replications] = F1.adult;
+}
-F2.replications[,N.replications] = second_generation.population;
+if (process_diapausing_adults | process_total_adults) {
-F2_adults.replications[,N.replications] = F2.adult;
+Diapausing.replications[,current_replication] = Diapausing;
+}
-population.replications[,N.replications] = total.population;
+newborn.replications[,current_replication] = N.newborn;
-}
+adult.replications[,current_replication] = N.adult;
+death.replications[,current_replication] = N.death;
-# Mean value for eggs.
+if (plot_generations_separately) {
-eggs = apply(Eggs.replications, 1, mean);
+# P is Parental, or overwintered adults.
-# Standard error for eggs.
+P.replications[,current_replication] = overwintering_adult.population;
-eggs.std_error = apply(Eggs.replications, 1, sd) / sqrt(opt$replications);
+# F1 is the first field-produced generation.
+F1.replications[,current_replication] = first_generation.population;
-# Mean value for nymphs.
+# F2 is the second field-produced generation.
-nymphs = apply((YoungNymphs.replications+OldNymphs.replications), 1, mean);
+F2.replications[,current_replication] = second_generation.population;
-# Standard error for nymphs.
+if (process_eggs) {
-nymphs.std_error = apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd);
+P_eggs.replications[,current_replication] = P.egg;
+F1_eggs.replications[,current_replication] = F1.egg;
-# Mean value for adults.
+F2_eggs.replications[,current_replication] = F2.egg;
-adults = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, mean);
+}
-# Standard error for adults.
+if (process_young_nymphs) {
-adults.std_error = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications);
+P_young_nymphs.replications[,current_replication] = P.young_nymph;
+F1_young_nymphs.replications[,current_replication] = F1.young_nymph;
-# Mean value for P.
+F2_young_nymphs.replications[,current_replication] = F2.young_nymph;
-P = apply(P.replications, 1, mean);
+}
-# Standard error for P.
+if (process_old_nymphs) {
-P.std_error = apply(P.replications, 1, sd) / sqrt(opt$replications);
+P_old_nymphs.replications[,current_replication] = P.old_nymph;
+F1_old_nymphs.replications[,current_replication] = F1.old_nymph;
-# Mean value for P adults.
+F2_old_nymphs.replications[,current_replication] = F2.old_nymph;
-P_adults = apply(P_adults.replications, 1, mean);
+}
-# Standard error for P_adult.
+if (process_total_nymphs) {
-P_adults.std_error = apply(P_adults.replications, 1, sd) / sqrt(opt$replications);
+P_total_nymphs.replications[,current_replication] = P.total_nymph;
+F1_total_nymphs.replications[,current_replication] = F1.total_nymph;
-# Mean value for F1.
+F2_total_nymphs.replications[,current_replication] = F2.total_nymph;
-F1 = apply(F1.replications, 1, mean);
+}
-# Standard error for F1.
+if (process_previttelogenic_adults) {
-F1.std_error = apply(F1.replications, 1, sd) / sqrt(opt$replications);
+P_previttelogenic_adults.replications[,current_replication] = P.previttelogenic_adult;
+F1_previttelogenic_adults.replications[,current_replication] = F1.previttelogenic_adult;
-# Mean value for F1 adults.
+F2_previttelogenic_adults.replications[,current_replication] = F2.previttelogenic_adult;
-F1_adults = apply(F1_adults.replications, 1, mean);
+}
-# Standard error for F1 adult.
+if (process_vittelogenic_adults) {
-F1_adults.std_error = apply(F1_adults.replications, 1, sd) / sqrt(opt$replications);
+P_vittelogenic_adults.replications[,current_replication] = P.vittelogenic_adult;
+F1_vittelogenic_adults.replications[,current_replication] = F1.vittelogenic_adult;
-# Mean value for F2.
+F2_vittelogenic_adults.replications[,current_replication] = F2.vittelogenic_adult;
-F2 = apply(F2.replications, 1, mean);
+}
-# Standard error for F2.
+if (process_diapausing_adults) {
-F2.std_error = apply(F2.replications, 1, sd) / sqrt(opt$replications);
+P_diapausing_adults.replications[,current_replication] = P.diapausing_adult;
+F1_diapausing_adults.replications[,current_replication] = F1.diapausing_adult;
-# Mean value for F2 adults.
+F2_diapausing_adults.replications[,current_replication] = F2.diapausing_adult;
-F2_adults = apply(F2_adults.replications, 1, mean);
+}
-# Standard error for F2 adult.
+if (process_total_adults) {
-F2_adults.std_error = apply(F2_adults.replications, 1, sd) / sqrt(opt$replications);
+P_total_adults.replications[,current_replication] = P.total_adult;
+F1_total_adults.replications[,current_replication] = F1.total_adult;
-# Display the total number of days in the Galaxy history item blurb.
+F2_total_adults.replications[,current_replication] = F2.total_adult;
-cat("Number of days: ", opt$num_days, "\n");
+}
+}
-dev.new(width=20, height=30);
+population.replications[,current_replication] = total.population;
+# End processing replications.
-# Start PDF device driver to save charts to output.
+}
-pdf(file=opt$output, width=20, height=30, bg="white");
-par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+if (process_eggs) {
+# Mean value for eggs.
-# Data analysis and visualization plots only within a single calendar year.
+eggs = apply(Eggs.replications, 1, mean);
-days = c(1:opt$num_days);
+temperature_data_frame = append_vector(temperature_data_frame, eggs, "EGG");
-start_date = temperature_data_frame$DATE[1];
+# Standard error for eggs.
-end_date = temperature_data_frame$DATE[opt$num_days];
+eggs.std_error = apply(Eggs.replications, 1, sd) / sqrt(opt$replications);
+temperature_data_frame = append_vector(temperature_data_frame, eggs.std_error, "EGGSE");
-# Subfigure 1: population size by life stage.
+}
-maxval = max(eggs+eggs.std_error, nymphs+nymphs.std_error, adults+adults.std_error);
+if (process_nymphs) {
-render_chart("pop_size_by_life_stage", opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
+# Calculate nymph populations for selected life stage.
-opt$std_error_plot, adults, nymphs, eggs, adults.std_error, nymphs.std_error, eggs.std_error, date_labels);
+for (life_stage_nymph in life_stages_nymph) {
-# Subfigure 2: population size by generation.
+if (life_stage_nymph=="Young") {
-maxval = max(F2);
+# Mean value for young nymphs.
-render_chart("pop_size_by_generation", opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
+young_nymphs = apply(YoungNymphs.replications, 1, mean);
-opt$std_error_plot, P, F1, F2, P.std_error, F1.std_error, F2.std_error, date_labels);
+temperature_data_frame = append_vector(temperature_data_frame, young_nymphs, "YOUNGNYMPH");
-# Subfigure 3: adult population size by generation.
+# Standard error for young nymphs.
-maxval = max(F2_adults) + 100;
+young_nymphs.std_error = apply(YoungNymphs.replications / sqrt(opt$replications), 1, sd);
-render_chart("adult_pop_size_by_generation", opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
+temperature_data_frame = append_vector(temperature_data_frame, young_nymphs.std_error, "YOUNGNYMPHSE");
-opt$std_error_plot, P_adults, F1_adults, F2_adults, P_adults.std_error, F1_adults.std_error, F2_adults.std_error,
+} else if (life_stage_nymph=="Old") {
-date_labels);
+# Mean value for old nymphs.
+old_nymphs = apply(OldNymphs.replications, 1, mean);
-# Turn off device driver to flush output.
+temperature_data_frame = append_vector(temperature_data_frame, old_nymphs, "OLDNYMPH");
-dev.off();
+# Standard error for old nymphs.
+old_nymphs.std_error = apply(OldNymphs.replications / sqrt(opt$replications), 1, sd);
+temperature_data_frame = append_vector(temperature_data_frame, old_nymphs.std_error, "OLDNYMPHSE");
+} else if (life_stage_nymph=="Total") {
+# Mean value for all nymphs.
+total_nymphs = apply((YoungNymphs.replications+OldNymphs.replications), 1, mean);
+temperature_data_frame = append_vector(temperature_data_frame, total_nymphs, "TOTALNYMPH");
+# Standard error for all nymphs.
+total_nymphs.std_error = apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd);
+temperature_data_frame = append_vector(temperature_data_frame, total_nymphs.std_error, "TOTALNYMPHSE");
+}
+}
+}
+if (process_adults) {
+# Calculate adult populations for selected life stage.
+for (life_stage_adult in life_stages_adult) {
+if (life_stage_adult == "Pre-vittelogenic") {
+# Mean value for previttelogenic adults.
+previttelogenic_adults = apply(Previttelogenic.replications, 1, mean);
+temperature_data_frame = append_vector(temperature_data_frame, previttelogenic_adults, "PRE-VITADULT");
+# Standard error for previttelogenic adults.
+previttelogenic_adults.std_error = apply(Previttelogenic.replications, 1, sd) / sqrt(opt$replications);
+temperature_data_frame = append_vector(temperature_data_frame, previttelogenic_adults.std_error, "PRE-VITADULTSE");
+} else if (life_stage_adult == "Vittelogenic") {
+# Mean value for vittelogenic adults.
+vittelogenic_adults = apply(Vittelogenic.replications, 1, mean);
+temperature_data_frame = append_vector(temperature_data_frame, vittelogenic_adults, "VITADULT");
+# Standard error for vittelogenic adults.
+vittelogenic_adults.std_error = apply(Vittelogenic.replications, 1, sd) / sqrt(opt$replications);
+temperature_data_frame = append_vector(temperature_data_frame, vittelogenic_adults.std_error, "VITADULTSE");
+} else if (life_stage_adult == "Diapausing") {
+# Mean value for vittelogenic adults.
+diapausing_adults = apply(Diapausing.replications, 1, mean);
+temperature_data_frame = append_vector(temperature_data_frame, diapausing_adults, "DIAPAUSINGADULT");
+# Standard error for vittelogenic adults.
+diapausing_adults.std_error = apply(Diapausing.replications, 1, sd) / sqrt(opt$replications);
+temperature_data_frame = append_vector(temperature_data_frame, diapausing_adults.std_error, "DIAPAUSINGADULTSE");
+} else if (life_stage_adult=="Total") {
+# Mean value for all adults.
+total_adults = apply((Previttelogenic.replications+Vittelogenic.replications+Diapausing.replications), 1, mean);
+temperature_data_frame = append_vector(temperature_data_frame, total_adults, "TOTALADULT");
+# Standard error for all adults.
+total_adults.std_error = apply((Previttelogenic.replications+Vittelogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications);
+temperature_data_frame = append_vector(temperature_data_frame, total_adults.std_error, "TOTALADULTSE");
+}
+}
+}
+if (plot_generations_separately) {
+m_se = get_mean_and_std_error(P.replications, F1.replications, F2.replications);
+P = m_se[[1]];
+P.std_error = m_se[[2]];
+F1 = m_se[[3]];
+F1.std_error = m_se[[4]];
+F2 = m_se[[5]];
+F2.std_error = m_se[[6]];
+if (process_eggs) {
+m_se = get_mean_and_std_error(P_eggs.replications, F1_eggs.replications, F2_eggs.replications);
+P_eggs = m_se[[1]];
+P_eggs.std_error = m_se[[2]];
+temperature_data_frame_P = append_vector(temperature_data_frame_P, P_eggs, "EGG-P");
+temperature_data_frame_P = append_vector(temperature_data_frame_P, P_eggs.std_error, "EGG-P-SE");
+F1_eggs = m_se[[3]];
+F1_eggs.std_error = m_se[[4]];
+temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_eggs, "EGG-F1");
+temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_eggs.std_error, "EGG-F1-SE");
+F2_eggs = m_se[[5]];
+F2_eggs.std_error = m_se[[6]];
+temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_eggs, "EGG-F2");
+temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_eggs.std_error, "EGG-F2-SE");
+}
+if (process_young_nymphs) {
+m_se = get_mean_and_std_error(P_young_nymphs.replications, F1_young_nymphs.replications, F2_young_nymphs.replications);
+P_young_nymphs = m_se[[1]];
+P_young_nymphs.std_error = m_se[[2]];
+temperature_data_frame_P = append_vector(temperature_data_frame_P, P_young_nymphs, "YOUNGNYMPH-P");
+temperature_data_frame_P = append_vector(temperature_data_frame_P, P_young_nymphs.std_error, "YOUNGNYMPH-P-SE");
+F1_young_nymphs = m_se[[3]];
+F1_young_nymphs.std_error = m_se[[4]];
+temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_young_nymphs, "YOUNGNYMPH-F1");
+temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_young_nymphs.std_error, "YOUNGNYMPH-F1-SE");
+F2_young_nymphs = m_se[[5]];
+F2_young_nymphs.std_error = m_se[[6]];
+temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_young_nymphs, "YOUNGNYMPH-F2");
+temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_young_nymphs.std_error, "YOUNGNYMPH-F2-SE");
+}
+if (process_old_nymphs) {
+m_se = get_mean_and_std_error(P_old_nymphs.replications, F1_old_nymphs.replications, F2_old_nymphs.replications);
+P_old_nymphs = m_se[[1]];
+P_old_nymphs.std_error = m_se[[2]];
+temperature_data_frame_P = append_vector(temperature_data_frame_P, P_old_nymphs, "OLDNYMPH-P");
+temperature_data_frame_P = append_vector(temperature_data_frame_P, P_old_nymphs.std_error, "OLDNYMPH-P-SE");
+F1_old_nymphs = m_se[[3]];
+F1_old_nymphs.std_error = m_se[[4]];
+temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_old_nymphs, "OLDNYMPH-F1");
+temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_old_nymphs.std_error, "OLDNYMPH-F1-SE");
+F2_old_nymphs = m_se[[5]];
+F2_old_nymphs.std_error = m_se[[6]];
+temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_old_nymphs, "OLDNYMPH-F2");
+temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_old_nymphs.std_error, "OLDNYMPH-F2-SE");
+}
+if (process_total_nymphs) {
+m_se = get_mean_and_std_error(P_total_nymphs.replications, F1_total_nymphs.replications, F2_total_nymphs.replications);
+P_total_nymphs = m_se[[1]];
+P_total_nymphs.std_error = m_se[[2]];
+temperature_data_frame_P = append_vector(temperature_data_frame_P, P_total_nymphs, "TOTALNYMPH-P");
+temperature_data_frame_P = append_vector(temperature_data_frame_P, P_total_nymphs.std_error, "TOTALNYMPH-P-SE");
+F1_total_nymphs = m_se[[3]];
+F1_total_nymphs.std_error = m_se[[4]];
+temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_total_nymphs, "TOTALNYMPH-F1");
+temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_total_nymphs.std_error, "TOTALNYMPH-F1-SE");
+F2_total_nymphs = m_se[[5]];
+F2_total_nymphs.std_error = m_se[[6]];
+temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_total_nymphs, "TOTALNYMPH-F2");
+temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_total_nymphs.std_error, "TOTALNYMPH-F2-SE");
+}
+if (process_previttelogenic_adults) {
+m_se = get_mean_and_std_error(P_previttelogenic_adults.replications, F1_previttelogenic_adults.replications, F2_previttelogenic_adults.replications);
+P_previttelogenic_adults = m_se[[1]];
+P_previttelogenic_adults.std_error = m_se[[2]];
+temperature_data_frame_P = append_vector(temperature_data_frame_P, P_previttelogenic_adults, "PRE-VITADULT-P");
+temperature_data_frame_P = append_vector(temperature_data_frame_P, P_previttelogenic_adults.std_error, "PRE-VITADULT-P-SE");
+F1_previttelogenic_adults = m_se[[3]];
+F1_previttelogenic_adults.std_error = m_se[[4]];
+temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_previttelogenic_adults, "PRE-VITADULT-F1");
+temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_previttelogenic_adults.std_error, "PRE-VITADULT-F1-SE");
+F2_previttelogenic_adults = m_se[[5]];
+F2_previttelogenic_adults.std_error = m_se[[6]];
+temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_previttelogenic_adults, "PRE-VITADULT-F2");
+temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_previttelogenic_adults.std_error, "PRE-VITADULT-F2-SE");
+}
+if (process_vittelogenic_adults) {
+m_se = get_mean_and_std_error(P_vittelogenic_adults.replications, F1_vittelogenic_adults.replications, F2_vittelogenic_adults.replications);
+P_vittelogenic_adults = m_se[[1]];
+P_vittelogenic_adults.std_error = m_se[[2]];
+temperature_data_frame_P = append_vector(temperature_data_frame_P, P_vittelogenic_adults, "VITADULT-P");
+temperature_data_frame_P = append_vector(temperature_data_frame_P, P_vittelogenic_adults.std_error, "VITADULT-P-SE");
+F1_vittelogenic_adults = m_se[[3]];
+F1_vittelogenic_adults.std_error = m_se[[4]];
+temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_vittelogenic_adults, "VITADULT-F1");
+temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_vittelogenic_adults.std_error, "VITADULT-F1-SE");
+F2_vittelogenic_adults = m_se[[5]];
+F2_vittelogenic_adults.std_error = m_se[[6]];
+temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_vittelogenic_adults, "VITADULT-F2");
+temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_vittelogenic_adults.std_error, "VITADULT-F2-SE");
+}
+if (process_diapausing_adults) {
+m_se = get_mean_and_std_error(P_diapausing_adults.replications, F1_diapausing_adults.replications, F2_diapausing_adults.replications);
+P_diapausing_adults = m_se[[1]];
+P_diapausing_adults.std_error = m_se[[2]];
+temperature_data_frame_P = append_vector(temperature_data_frame_P, P_diapausing_adults, "DIAPAUSINGADULT-P");
+temperature_data_frame_P = append_vector(temperature_data_frame_P, P_diapausing_adults.std_error, "DIAPAUSINGADULT-P-SE");
+F1_diapausing_adults = m_se[[3]];
+F1_diapausing_adults.std_error = m_se[[4]];
+temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_diapausing_adults, "DIAPAUSINGADULT-F1");
+temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_diapausing_adults.std_error, "DIAPAUSINGADULT-F1-SE");
+F2_diapausing_adults = m_se[[5]];
+F2_diapausing_adults.std_error = m_se[[6]];
+temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_diapausing_adults, "DIAPAUSINGADULT-F2");
+temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_diapausing_adults.std_error, "DIAPAUSINGADULT-F2-SE");
+}
+if (process_total_adults) {
+m_se = get_mean_and_std_error(P_total_adults.replications, F1_total_adults.replications, F2_total_adults.replications);
+P_total_adults = m_se[[1]];
+P_total_adults.std_error = m_se[[2]];
+temperature_data_frame_P = append_vector(temperature_data_frame_P, P_total_adults, "TOTALADULT-P");
+temperature_data_frame_P = append_vector(temperature_data_frame_P, P_total_adults.std_error, "TOTALADULT-P-SE");
+F1_total_adults = m_se[[3]];
+F1_total_adults.std_error = m_se[[4]];
+temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_total_adults, "TOTALADULT-F1");
+temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_total_adults.std_error, "TOTALADULT-F1-SE");
+F2_total_adults = m_se[[5]];
+F2_total_adults.std_error = m_se[[6]];
+temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_total_adults, "TOTALADULT-F2");
+temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_total_adults.std_error, "TOTALADULT-F2-SE");
+}
+}
+# Save the analyzed data for combined generations.
+file_path = paste("output_data_dir", "04_combined_generations.csv", sep="/");
+write.csv(temperature_data_frame, file=file_path, row.names=F);
+if (plot_generations_separately) {
+# Save the analyzed data for generation P.
+file_path = paste("output_data_dir", "01_generation_P.csv", sep="/");
+write.csv(temperature_data_frame_P, file=file_path, row.names=F);
+# Save the analyzed data for generation F1.
+file_path = paste("output_data_dir", "02_generation_F1.csv", sep="/");
+write.csv(temperature_data_frame_F1, file=file_path, row.names=F);
+# Save the analyzed data for generation F2.
+file_path = paste("output_data_dir", "03_generation_F2.csv", sep="/");
+write.csv(temperature_data_frame_F2, file=file_path, row.names=F);
+}
+if (plot_generations_separately) {
+for (life_stage in life_stages) {
+if (life_stage == "Egg") {
+# Start PDF device driver.
+dev.new(width=20, height=30);
+file_path = get_file_path(life_stage, "egg_pop_by_generation.pdf")
+pdf(file=file_path, width=20, height=30, bg="white");
+par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+# Egg population size by generation.
+maxval = max(P_eggs+F1_eggs+F2_eggs) + 100;
+render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, location, latitude,
+start_date, end_date, total_days_vector, maxval, 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, group3_std_error=F2_eggs.std_error);
+# Turn off device driver to flush output.
+dev.off();
+} else if (life_stage == "Nymph") {
+for (life_stage_nymph in life_stages_nymph) {
+# Start PDF device driver.
+dev.new(width=20, height=30);
+file_path = get_file_path(life_stage, "nymph_pop_by_generation.pdf", life_stage_nymph=life_stage_nymph)
+pdf(file=file_path, width=20, height=30, bg="white");
+par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+if (life_stage_nymph=="Young") {
+# Young nymph population size by generation.
+maxval = max(P_young_nymphs+F1_young_nymphs+F2_young_nymphs) + 100;
+group = P_young_nymphs;
+group_std_error = P_young_nymphs.std_error;
+group2 = F1_young_nymphs;
+group2_std_error = F1_young_nymphs.std_error;
+group3 = F2_young_nymphs;
+group3_std_error = F2_young_nymphs.std_error;
+} else if (life_stage_nymph=="Old") {
+# Total nymph population size by generation.
+maxval = max(P_old_nymphs+F1_old_nymphs+F2_old_nymphs) + 100;
+group = P_old_nymphs;
+group_std_error = P_old_nymphs.std_error;
+group2 = F1_old_nymphs;
+group2_std_error = F1_old_nymphs.std_error;
+group3 = F2_old_nymphs;
+group3_std_error = F2_old_nymphs.std_error;
+} else if (life_stage_nymph=="Total") {
+# Total nymph population size by generation.
+maxval = max(P_total_nymphs+F1_total_nymphs+F2_total_nymphs) + 100;
+group = P_total_nymphs;
+group_std_error = P_total_nymphs.std_error;
+group2 = F1_total_nymphs;
+group2_std_error = F1_total_nymphs.std_error;
+group3 = F2_total_nymphs;
+group3_std_error = F2_total_nymphs.std_error;
+}
+render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, location, latitude,
+start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
+group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, life_stages_nymph=life_stage_nymph);
+# Turn off device driver to flush output.
+dev.off();
+}
+} else if (life_stage == "Adult") {
+for (life_stage_adult in life_stages_adult) {
+# Start PDF device driver.
+dev.new(width=20, height=30);
+file_path = get_file_path(life_stage, "adult_pop_by_generation.pdf", life_stage_adult=life_stage_adult)
+pdf(file=file_path, width=20, height=30, bg="white");
+par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+if (life_stage_adult=="Pre-vittelogenic") {
+# Pre-vittelogenic adult population size by generation.
+maxval = max(P_previttelogenic_adults+F1_previttelogenic_adults+F2_previttelogenic_adults) + 100;
+group = P_previttelogenic_adults;
+group_std_error = P_previttelogenic_adults.std_error;
+group2 = F1_previttelogenic_adults;
+group2_std_error = F1_previttelogenic_adults.std_error;
+group3 = F2_previttelogenic_adults;
+group3_std_error = F2_previttelogenic_adults.std_error;
+} else if (life_stage_adult=="Vittelogenic") {
+# Vittelogenic adult population size by generation.
+maxval = max(P_vittelogenic_adults+F1_vittelogenic_adults+F2_vittelogenic_adults) + 100;
+group = P_vittelogenic_adults;
+group_std_error = P_vittelogenic_adults.std_error;
+group2 = F1_vittelogenic_adults;
+group2_std_error = F1_vittelogenic_adults.std_error;
+group3 = F2_vittelogenic_adults;
+group3_std_error = F2_vittelogenic_adults.std_error;
+} else if (life_stage_adult=="Diapausing") {
+# Diapausing adult population size by generation.
+maxval = max(P_diapausing_adults+F1_diapausing_adults+F2_diapausing_adults) + 100;
+group = P_diapausing_adults;
+group_std_error = P_diapausing_adults.std_error;
+group2 = F1_diapausing_adults;
+group2_std_error = F1_diapausing_adults.std_error;
+group3 = F2_diapausing_adults;
+group3_std_error = F2_diapausing_adults.std_error;
+} else if (life_stage_adult=="Total") {
+# Total adult population size by generation.
+maxval = max(P_total_adults+F1_total_adults+F2_total_adults) + 100;
+group = P_total_adults;
+group_std_error = P_total_adults.std_error;
+group2 = F1_total_adults;
+group2_std_error = F1_total_adults.std_error;
+group3 = F2_total_adults;
+group3_std_error = F2_total_adults.std_error;
+}
+render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, location, latitude,
+start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
+group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, life_stages_adult=life_stage_adult);
+# Turn off device driver to flush output.
+dev.off();
+}
+} else if (life_stage == "Total") {
+# Start PDF device driver.
+# Name collection elements so that they
+# are displayed in logical order.
+dev.new(width=20, height=30);
+file_path = get_file_path(life_stage, "total_pop_by_generation.pdf")
+pdf(file=file_path, width=20, height=30, bg="white");
+par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+# Total population size by generation.
+maxval = max(P+F1+F2) + 100;
+render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, location, latitude,
+start_date, end_date, total_days_vector, maxval, 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);
+# Turn off device driver to flush output.
+dev.off();
+}
+}
+} else {
+for (life_stage in life_stages) {
+if (life_stage == "Egg") {
+# Start PDF device driver.
+dev.new(width=20, height=30);
+file_path = get_file_path(life_stage, "egg_pop.pdf")
+pdf(file=file_path, width=20, height=30, bg="white");
+par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+# Egg population size.
+maxval = max(eggs+eggs.std_error) + 100;
+render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, location, latitude,
+start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=eggs, group_std_error=eggs.std_error);
+# Turn off device driver to flush output.
+dev.off();
+} else if (life_stage == "Nymph") {
+for (life_stage_nymph in life_stages_nymph) {
+# Start PDF device driver.
+dev.new(width=20, height=30);
+file_path = get_file_path(life_stage, "nymph_pop.pdf", life_stage_nymph=life_stage_nymph)
+pdf(file=file_path, width=20, height=30, bg="white");
+par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+if (life_stage_nymph=="Total") {
+# Total nymph population size.
+group = total_nymphs;
+group_std_error = total_nymphs.std_error;
+} else if (life_stage_nymph=="Young") {
+# Young nymph population size.
+group = young_nymphs;
+group_std_error = young_nymphs.std_error;
+} else if (life_stage_nymph=="Old") {
+# Old nymph population size.
+group = old_nymphs;
+group_std_error = old_nymphs.std_error;
+}
+maxval = max(group+group_std_error) + 100;
+render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, location, latitude,
+start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
+life_stages_nymph=life_stage_nymph);
+# Turn off device driver to flush output.
+dev.off();
+}
+} else if (life_stage == "Adult") {
+for (life_stage_adult in life_stages_adult) {
+# Start PDF device driver.
+dev.new(width=20, height=30);
+file_path = get_file_path(life_stage, "adult_pop.pdf", life_stage_adult=life_stage_adult)
+pdf(file=file_path, width=20, height=30, bg="white");
+par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+if (life_stage_adult=="Total") {
+# Total adult population size.
+group = total_adults;
+group_std_error = total_adults.std_error
+} else if (life_stage_adult=="Pre-vittelogenic") {
+# Pre-vittelogenic adult population size.
+group = previttelogenic_adults;
+group_std_error = previttelogenic_adults.std_error
+} else if (life_stage_adult=="Vittelogenic") {
+# Vittelogenic adult population size.
+group = vittelogenic_adults;
+group_std_error = vittelogenic_adults.std_error
+} else if (life_stage_adult=="Diapausing") {
+# Diapausing adult population size.
+group = diapausing_adults;
+group_std_error = diapausing_adults.std_error
+}
+maxval = max(group+group_std_error) + 100;
+render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, location, latitude,
+start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
+life_stages_adult=life_stage_adult);
+# Turn off device driver to flush output.
+dev.off();
+}
+} else if (life_stage == "Total") {
+# Start PDF device driver.
+dev.new(width=20, height=30);
+file_path = get_file_path(life_stage, "total_pop.pdf")
+pdf(file=file_path, width=20, height=30, bg="white");
+par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+# Total population size.
+maxval = max(eggs+eggs.std_error, total_nymphs+total_nymphs.std_error, total_adults+total_adults.std_error) + 100;
+render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, location, latitude,
+start_date, end_date, total_days_vector, maxval, 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, group3_std_error=eggs.std_error);
+# Turn off device driver to flush output.
+dev.off();
+}
+}
+}

Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 112:bcb12b7e8563 draft