Mercurial > repos > greg > insect_phenology_model
diff insect_phenology_model.R @ 112:bcb12b7e8563 draft
Uploaded
| author | greg |
|---|---|
| date | Tue, 29 May 2018 09:00:25 -0400 |
| parents | 37ac68b6ff10 |
| children | 9c998fd06628 |
line wrap: on
line diff
--- a/insect_phenology_model.R Tue Feb 13 13:47:32 2018 -0500 +++ b/insect_phenology_model.R Tue May 29 09:00:25 2018 -0400 @@ -6,20 +6,24 @@ 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)") ) @@ -27,9 +31,9 @@ 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]; @@ -45,52 +49,98 @@ 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) { - dev.rate = -0.9843 * temperature + 33.438; - return(dev.rate); +append_vector = function(data_frame, vec, new_column_name) { + num_columns = dim(data_frame)[2]; + current_column_names = colnames(data_frame); + # Append vector vec as a new column to data_frame. + data_frame[,num_columns+1] = vec; + # Reset the column names with the additional column for later access. + colnames(data_frame) = append(current_column_names, new_column_name); + return(data_frame); } -dev.emerg = function(temperature) { - emerg.rate = -0.5332 * temperature + 24.147; - return(emerg.rate); -} - -dev.old = function(temperature) { - n34 = -0.6119 * temperature + 17.602; - n45 = -0.4408 * temperature + 19.036; - dev.rate = mean(n34 + n45); - return(dev.rate); +get_file_path = function(life_stage, base_name, life_stage_nymph=NULL, life_stage_adult=NULL) { + if (!is.null(life_stage_nymph)) { + lsi = get_life_stage_index(life_stage, life_stage_nymph=life_stage_nymph); + file_name = paste(lsi, tolower(life_stage_nymph), base_name, sep="_"); + } 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="_"); + } else { + lsi = get_life_stage_index(life_stage); + file_name = paste(lsi, base_name, sep="_"); + } + file_path = paste("output_plots_dir", file_name, sep="/"); + return(file_path); } -dev.young = function(temperature) { - n12 = -0.3728 * temperature + 14.68; - n23 = -0.6119 * temperature + 25.249; - dev.rate = mean(n12 + n23); - return(dev.rate); +get_life_stage_index = function(life_stage, life_stage_nymph=NULL, life_stage_adult=NULL) { + # Name collection elements so that they + # are displayed in logical order. + if (life_stage=="Egg") { + lsi = "01"; + } else if (life_stage=="Nymph") { + if (life_stage_nymph=="Young") { + lsi = "02"; + } else if (life_stage_nymph=="Old") { + lsi = "03"; + } else if (life_stage_nymph=="Total") { + lsi="04"; + } + } else if (life_stage=="Adult") { + if (life_stage_adult=="Pre-vittelogenic") { + lsi = "05"; + } 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_date_labels = function(temperature_data_frame, num_rows) { - # Keep track of the years to see if spanning years. - month_labels = list(); - current_month_label = NULL; - for (i in 1:num_rows) { - # Get the year and month from the date which - # has the format YYYY-MM-DD. - date = format(temperature_data_frame$DATE[i]); - items = strsplit(date, "-")[[1]]; - month = items[2]; - month_label = month.abb[as.integer(month)]; - if (!identical(current_month_label, month_label)) { - month_labels[length(month_labels)+1] = month_label; - current_month_label = month_label; - } +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; } - return(c(unlist(month_labels))); + 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) { @@ -164,6 +214,112 @@ 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; @@ -197,125 +353,432 @@ return(mortality.probability); } -parse_input_data = function(input_file, num_rows) { - # Read in the input temperature datafile into a data frame. - temperature_data_frame = read.csv(file=input_file, header=T, strip.white=TRUE, sep=","); - num_columns = dim(temperature_data_frame)[2]; - if (num_columns == 6) { - # The input data has the following 6 columns: +parse_input_data = function(input_ytd, input_norm, num_days_ytd, location) { + if (is.null(input_ytd)) { + # We're analysing only the 30 year normals data, so create an empty + # data frame for containing temperature data after it is converted + # from the 30 year normals format to the year-to-date format. + 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.. - colnames(temperature_data_frame) = c("LATITUDE","LONGITUDE", "DATE", "DOY", "TMIN", "TMAX"); - # Add a column containing the daylight length for each day. - temperature_data_frame = add_daylight_length(temperature_data_frame, num_columns, num_rows); - # Reset the column names with the additional column for later access. - colnames(temperature_data_frame) = c("LATITUDE","LONGITUDE", "DATE", "DOY", "TMIN", "TMAX", "DAYLEN"); + temperature_data_frame = read.csv(file=input_ytd, header=T, strip.white=TRUE, stringsAsFactors=FALSE, sep=","); + # Set the temperature_data_frame column names for access. + colnames(temperature_data_frame) = c("LATITUDE", "LONGITUDE", "DATE", "DOY", "TMIN", "TMAX"); + # Get the start date. + start_date = temperature_data_frame$DATE[1]; + end_date = temperature_data_frame$DATE[num_days_ytd]; + # Extract the year from the start date. + 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. + start_doy_ytd = as.integer(temperature_data_frame$DOY[1]); + end_doy_ytd = as.integer(temperature_data_frame$DOY[num_days_ytd]); + } + # See if we're in a leap year. + is_leap_year = is_leap_year(start_date); + # Get the number of days in the year. + total_days = get_total_days(is_leap_year); + # 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),]; } - return(temperature_data_frame); + # 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(chart_type, insect, location, latitude, start_date, end_date, days, maxval, plot_std_error, - group1, group2, group3, group1_std_error, group2_std_error, group3_std_error, date_labels) { - if (chart_type == "pop_size_by_life_stage") { - title = paste(insect, ": Total pop. by life stage :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" "); - legend_text = c("Egg", "Nymph", "Adult"); - columns = c(4, 2, 1); - } else if (chart_type == "pop_size_by_generation") { - title = paste(insect, ": Total pop. by generation :", 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") { - title = paste(insect, ": Adult pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" "); +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); - } - 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); - 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(1, at=c(1:length(date_labels)) * 30 - 15, cex.axis=3, labels=date_labels); - axis(2, cex.axis=3); - if (plot_std_error==1) { - # Standard error for group1. - lines(days, group1+group1_std_error, lty=2); - lines(days, group1-group1_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); + 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); + } } } -temperature_data_frame = parse_input_data(opt$input, opt$num_days); -# All latitude values are the same, so get the value from the first row. +# Determine if we're plotting generations separately. +if (opt$plot_generations_separately=="yes") { + plot_generations_separately = TRUE; +} else { + plot_generations_separately = FALSE; +} +# 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); -YoungNymphs.replications = matrix(rep(0, opt$num_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); -Vitellogenic.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); -Diapausing.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); - -newborn.replications = matrix(rep(0, opt$num_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); - -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); -F1.replications = matrix(rep(0, opt$num_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); -F2_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); - -population.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); +if (process_eggs) { + Eggs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications); +} +if (process_young_nymphs | process_total_nymphs) { + YoungNymphs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications); +} +if (process_old_nymphs | process_total_nymphs) { + OldNymphs.replications = matrix(rep(0, total_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); +} +if (process_vittelogenic_adults | process_total_adults) { + Vittelogenic.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications); +} +if (process_diapausing_adults | process_total_adults) { + Diapausing.replications = matrix(rep(0, total_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); - YoungNymphs = rep(0, opt$num_days); - OldNymphs = rep(0, opt$num_days); - Previtellogenic = rep(0, opt$num_days); - Vitellogenic = rep(0, opt$num_days); - Diapausing = rep(0, opt$num_days); - - N.newborn = rep(0, opt$num_days); - N.adult = rep(0, opt$num_days); - N.death = rep(0, opt$num_days); - - overwintering_adult.population = rep(0, opt$num_days); - first_generation.population = rep(0, opt$num_days); - second_generation.population = rep(0, opt$num_days); - - P.adult = rep(0, opt$num_days); - F1.adult = rep(0, opt$num_days); - F2.adult = rep(0, opt$num_days); - + if (process_eggs) { + Eggs = rep(0, total_days); + } + if (process_young_nymphs | process_total_nymphs) { + YoungNymphs = rep(0, total_days); + } + if (process_old_nymphs | process_total_nymphs) { + OldNymphs = rep(0, total_days); + } + if (process_previttelogenic_adults | process_total_adults) { + Previttelogenic = rep(0, total_days); + } + if (process_vittelogenic_adults | process_total_adults) { + Vittelogenic = rep(0, total_days); + } + if (process_diapausing_adults | process_total_adults) { + Diapausing = rep(0, total_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, opt$num_days); - # All the days included in the input temperature dataset. - for (row in 1:opt$num_days) { + averages.day = rep(0, total_days); + # All the days included in the input_ytd temperature dataset. + for (row in 1:total_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; @@ -341,6 +804,7 @@ 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) { @@ -361,7 +825,7 @@ 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. @@ -369,7 +833,7 @@ 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; @@ -377,7 +841,7 @@ } } 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. @@ -492,7 +956,7 @@ } 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; @@ -535,19 +999,33 @@ # Update population size for the next day. num_insects = num_insects - num_insects.death + num_insects.newborn; - # Aggregate results by day. - # Egg population size. - Eggs[row] = sum(vector.matrix[,2]==0); - # Young nymph population size. - YoungNymphs[row] = sum(vector.matrix[,2]==1); - # Old nymph population size. - OldNymphs[row] = sum(vector.matrix[,2]==2); - # Previtellogenic population size. - Previtellogenic[row] = sum(vector.matrix[,2]==3); - # Vitellogenic population size. - Vitellogenic[row] = sum(vector.matrix[,2]==4); - # Diapausing population size. - Diapausing[row] = sum(vector.matrix[,2]==5); + # Aggregate results by day. Due to multiple transpose calls + # on vector.matrix above, the columns of vector.matrix + # are now Generation, Stage, degree-days, T, Diapause, + if (process_eggs) { + # For egg population size, column 2 (Stage), must be 0. + Eggs[row] = sum(vector.matrix[,2]==0); + } + if (process_young_nymphs | process_total_nymphs) { + # For young nymph population size, column 2 (Stage) must be 1. + YoungNymphs[row] = sum(vector.matrix[,2]==1); + } + if (process_old_nymphs | process_total_nymphs) { + # 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; @@ -558,117 +1036,637 @@ 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. - P.adult[row] = sum(vector.matrix[,1]==0); - # F1 adult 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)); - # F2 adult population size - 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)); - } # End of days specified in the input temperature data. + if (plot_generations_separately) { + if (process_eggs) { + # For egg life stage of generation P population size, + # column 1 (generation) is 0 and column 2 (Stage) is 0. + P.egg[row] = sum(vector.matrix[,1]==0 & vector.matrix[,2]==0); + # For egg life stage of generation F1 population size, + # 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; - YoungNymphs.replications[,N.replications] = YoungNymphs; - OldNymphs.replications[,N.replications] = OldNymphs; - Previtellogenic.replications[,N.replications] = Previtellogenic; - Vitellogenic.replications[,N.replications] = Vitellogenic; - Diapausing.replications[,N.replications] = Diapausing; + if (process_eggs) { + Eggs.replications[,current_replication] = Eggs; + } + if (process_young_nymphs | process_total_nymphs) { + YoungNymphs.replications[,current_replication] = YoungNymphs; + } + if (process_old_nymphs | process_total_nymphs) { + OldNymphs.replications[,current_replication] = OldNymphs; + } + if (process_previttelogenic_adults | process_total_adults) { + Previttelogenic.replications[,current_replication] = Previttelogenic; + } + if (process_vittelogenic_adults | process_total_adults) { + Vittelogenic.replications[,current_replication] = Vittelogenic; + } + if (process_diapausing_adults | process_total_adults) { + Diapausing.replications[,current_replication] = Diapausing; + } + newborn.replications[,current_replication] = N.newborn; + adult.replications[,current_replication] = N.adult; + death.replications[,current_replication] = N.death; + if (plot_generations_separately) { + # P is Parental, or overwintered adults. + P.replications[,current_replication] = overwintering_adult.population; + # F1 is the first field-produced generation. + F1.replications[,current_replication] = first_generation.population; + # F2 is the second field-produced generation. + F2.replications[,current_replication] = second_generation.population; + if (process_eggs) { + P_eggs.replications[,current_replication] = P.egg; + F1_eggs.replications[,current_replication] = F1.egg; + F2_eggs.replications[,current_replication] = F2.egg; + } + if (process_young_nymphs) { + P_young_nymphs.replications[,current_replication] = P.young_nymph; + F1_young_nymphs.replications[,current_replication] = F1.young_nymph; + F2_young_nymphs.replications[,current_replication] = F2.young_nymph; + } + if (process_old_nymphs) { + P_old_nymphs.replications[,current_replication] = P.old_nymph; + F1_old_nymphs.replications[,current_replication] = F1.old_nymph; + F2_old_nymphs.replications[,current_replication] = F2.old_nymph; + } + if (process_total_nymphs) { + P_total_nymphs.replications[,current_replication] = P.total_nymph; + F1_total_nymphs.replications[,current_replication] = F1.total_nymph; + F2_total_nymphs.replications[,current_replication] = F2.total_nymph; + } + if (process_previttelogenic_adults) { + P_previttelogenic_adults.replications[,current_replication] = P.previttelogenic_adult; + F1_previttelogenic_adults.replications[,current_replication] = F1.previttelogenic_adult; + F2_previttelogenic_adults.replications[,current_replication] = F2.previttelogenic_adult; + } + if (process_vittelogenic_adults) { + P_vittelogenic_adults.replications[,current_replication] = P.vittelogenic_adult; + F1_vittelogenic_adults.replications[,current_replication] = F1.vittelogenic_adult; + F2_vittelogenic_adults.replications[,current_replication] = F2.vittelogenic_adult; + } + if (process_diapausing_adults) { + P_diapausing_adults.replications[,current_replication] = P.diapausing_adult; + F1_diapausing_adults.replications[,current_replication] = F1.diapausing_adult; + F2_diapausing_adults.replications[,current_replication] = F2.diapausing_adult; + } + if (process_total_adults) { + P_total_adults.replications[,current_replication] = P.total_adult; + F1_total_adults.replications[,current_replication] = F1.total_adult; + F2_total_adults.replications[,current_replication] = F2.total_adult; + } + } + population.replications[,current_replication] = total.population; + # End processing replications. +} - newborn.replications[,N.replications] = N.newborn; - adult.replications[,N.replications] = N.adult; - death.replications[,N.replications] = N.death; - - P.replications[,N.replications] = overwintering_adult.population; - P_adults.replications[,N.replications] = P.adult; - F1.replications[,N.replications] = first_generation.population; - F1_adults.replications[,N.replications] = F1.adult; - F2.replications[,N.replications] = second_generation.population; - F2_adults.replications[,N.replications] = F2.adult; - - population.replications[,N.replications] = total.population; +if (process_eggs) { + # Mean value for eggs. + eggs = apply(Eggs.replications, 1, mean); + temperature_data_frame = append_vector(temperature_data_frame, eggs, "EGG"); + # Standard error for eggs. + eggs.std_error = apply(Eggs.replications, 1, sd) / sqrt(opt$replications); + temperature_data_frame = append_vector(temperature_data_frame, eggs.std_error, "EGGSE"); +} +if (process_nymphs) { + # Calculate nymph populations for selected life stage. + for (life_stage_nymph in life_stages_nymph) { + if (life_stage_nymph=="Young") { + # Mean value for young nymphs. + young_nymphs = apply(YoungNymphs.replications, 1, mean); + temperature_data_frame = append_vector(temperature_data_frame, young_nymphs, "YOUNGNYMPH"); + # Standard error for young nymphs. + young_nymphs.std_error = apply(YoungNymphs.replications / sqrt(opt$replications), 1, sd); + temperature_data_frame = append_vector(temperature_data_frame, young_nymphs.std_error, "YOUNGNYMPHSE"); + } else if (life_stage_nymph=="Old") { + # Mean value for old nymphs. + old_nymphs = apply(OldNymphs.replications, 1, mean); + temperature_data_frame = append_vector(temperature_data_frame, old_nymphs, "OLDNYMPH"); + # 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"); + } + } } -# Mean value for eggs. -eggs = apply(Eggs.replications, 1, mean); -# Standard error for eggs. -eggs.std_error = apply(Eggs.replications, 1, sd) / sqrt(opt$replications); - -# Mean value for nymphs. -nymphs = apply((YoungNymphs.replications+OldNymphs.replications), 1, mean); -# Standard error for nymphs. -nymphs.std_error = apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd); - -# Mean value for adults. -adults = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, mean); -# Standard error for adults. -adults.std_error = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications); +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"); + } +} -# Mean value for P. -P = apply(P.replications, 1, mean); -# Standard error for P. -P.std_error = apply(P.replications, 1, sd) / sqrt(opt$replications); - -# Mean value for P adults. -P_adults = apply(P_adults.replications, 1, mean); -# Standard error for P_adult. -P_adults.std_error = apply(P_adults.replications, 1, sd) / sqrt(opt$replications); - -# Mean value for F1. -F1 = apply(F1.replications, 1, mean); -# Standard error for F1. -F1.std_error = apply(F1.replications, 1, sd) / sqrt(opt$replications); - -# Mean value for F1 adults. -F1_adults = apply(F1_adults.replications, 1, mean); -# Standard error for F1 adult. -F1_adults.std_error = apply(F1_adults.replications, 1, sd) / sqrt(opt$replications); +# 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); +} -# Mean value for F2. -F2 = apply(F2.replications, 1, mean); -# Standard error for F2. -F2.std_error = apply(F2.replications, 1, sd) / sqrt(opt$replications); - -# Mean value for F2 adults. -F2_adults = apply(F2_adults.replications, 1, mean); -# Standard error for F2 adult. -F2_adults.std_error = apply(F2_adults.replications, 1, sd) / sqrt(opt$replications); - -# Display the total number of days in the Galaxy history item blurb. -cat("Number of days: ", opt$num_days, "\n"); - -dev.new(width=20, height=30); - -# 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)); - -# Data analysis and visualization plots only within a single calendar year. -days = c(1:opt$num_days); -start_date = temperature_data_frame$DATE[1]; -end_date = temperature_data_frame$DATE[opt$num_days]; - -# Subfigure 1: population size by life stage. -maxval = max(eggs+eggs.std_error, nymphs+nymphs.std_error, adults+adults.std_error); -render_chart("pop_size_by_life_stage", opt$insect, opt$location, latitude, start_date, end_date, days, maxval, - opt$std_error_plot, adults, nymphs, eggs, adults.std_error, nymphs.std_error, eggs.std_error, date_labels); -# Subfigure 2: population size by generation. -maxval = max(F2); -render_chart("pop_size_by_generation", opt$insect, opt$location, latitude, start_date, end_date, days, maxval, - opt$std_error_plot, P, F1, F2, P.std_error, F1.std_error, F2.std_error, date_labels); -# Subfigure 3: adult population size by generation. -maxval = max(F2_adults) + 100; -render_chart("adult_pop_size_by_generation", opt$insect, opt$location, latitude, start_date, end_date, days, maxval, - opt$std_error_plot, P_adults, F1_adults, F2_adults, P_adults.std_error, F1_adults.std_error, F2_adults.std_error, - date_labels); - -# Turn off device driver to flush output. -dev.off(); +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(); + } + } +}