# HG changeset patch
# User greg
# Date 1519756319 18000
# Node ID 61bc6bd8807d91b24680cbefd48f88156ec6af7d
# Parent  d9371485aaf519244c5ced8b00d69c05fd4c2889
Uploaded

diff -r d9371485aaf5 -r 61bc6bd8807d insect_phenology_model.R
--- a/insect_phenology_model.R	Tue Feb 27 13:31:49 2018 -0500
+++ b/insect_phenology_model.R	Tue Feb 27 13:31:59 2018 -0500
@@ -9,17 +9,20 @@
     make_option(c("--input"), action="store", dest="input", help="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("--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_stage_nymph"), action="store", dest="life_stage_nymph", default=NULL, help="Nymph life stages for plotting"),
     make_option(c("--location"), action="store", dest="location", 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("--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("--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("--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("--young_nymph_accumulation"), action="store", dest="young_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)")
 )
 
@@ -73,7 +76,6 @@
     return(dev.rate);
 }
 
-
 get_date_labels = function(temperature_data_frame, num_rows) {
     # Keep track of the years to see if spanning years.
     month_labels = list();
@@ -215,65 +217,170 @@
 }
 
 
-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(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_stage_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=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=="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_stage_nymph, "Nymph Pop :", sep=" ");
+                title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
+                legend_text = c(paste(life_stage_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=F, 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(1, at=c(1:length(date_labels)) * 30 - 15, cex.axis=3, labels=date_labels);
+            axis(2, cex.axis=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_stage_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=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=="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);
+        }
     }
 }
 
+# Determine if we're plotting generations separately.
+if (opt$plot_generations_separately=="yes") {
+    plot_generations_separately = TRUE;
+} else {
+    plot_generations_separately = FALSE;
+}
+# Read the temperature data into a data frame.
 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.
+output_dir = "output_dir";
+# Get the date labels for plots.
+date_labels = get_date_labels(temperature_data_frame, opt$num_days);
+# All latitude values are the same, so get the value for plots from the first row.
 latitude = temperature_data_frame$LATITUDE[1];
+# Get the number of days for plots.
 num_columns = dim(temperature_data_frame)[2];
-date_labels = get_date_labels(temperature_data_frame, opt$num_days);
+# 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_adults = FALSE;
+for (life_stage in life_stages) {
+    if (life_stage=="Total") {
+        process_eggs = TRUE;
+        process_nymphs = TRUE;
+        life_stage_nymph = "Total";
+        process_adults = TRUE;
+        life_stages_adult = "Total";
+    } else if (life_stage=="Egg") {
+        process_eggs = TRUE;
+    } else if (life_stage=="Nymph") {
+        process_nymphs = TRUE;
+        life_stage_nymph = opt$life_stage_nymph;
+    } else if (life_stage=="Adult") {
+        process_adults = TRUE;
+        life_stages_adult = opt$life_stages_adult;
+    }
+}
 
 # 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);
-
+if (process_eggs) {
+    Eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+}
+if (process_nymphs) {
+    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);
+}
+if (process_adults) {
+    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);
-
+if (plot_generations_separately) {
+    # P is Parental, or overwintered adults.
+    P.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+    # F1 is the first field-produced generation.
+    F1.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+    # F2 is the second field-produced generation.
+    F2.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+    if (process_eggs) {
+        P_eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+        F1_eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+        F2_eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+    }
+    if (process_nymphs) {
+        P_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+        F1_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+        F2_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+    }
+    if (process_adults) {
+        P_adults.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_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
+    }
+}
+# Total population.
 population.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications);
 
 # Process replications.
@@ -282,32 +389,53 @@
     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);
-
+    if (process_eggs) {
+        Eggs = rep(0, opt$num_days);
+    }
+    if (process_nymphs) {
+        YoungNymphs = rep(0, opt$num_days);
+        OldNymphs = rep(0, opt$num_days);
+    }
+    if (process_adults) {
+        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 (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, opt$num_days);
+            F1.egg = rep(0, opt$num_days);
+            F2.egg = rep(0, opt$num_days);
+        }
+        if (process_nymphs) {
+            P.nymph = rep(0, opt$num_days);
+            F1.nymph = rep(0, opt$num_days);
+            F2.nymph = rep(0, opt$num_days);
+        }
+        if (process_adults) {
+            P.adult = rep(0, opt$num_days);
+            F1.adult = rep(0, opt$num_days);
+            F2.adult = rep(0, opt$num_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) {
@@ -535,19 +663,27 @@
         # 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_nymphs) {
+            # For young nymph population size, column 2 (Stage) must be 1.
+            YoungNymphs[row] = sum(vector.matrix[,2]==1);
+            # For old nymph population size, column 2 (Stage) must be 2.
+            OldNymphs[row] = sum(vector.matrix[,2]==2);
+        }
+        if (process_adults) {
+            # For pre-vitellogenic population size, column 2 (Stage) must be 3.
+            Previtellogenic[row] = sum(vector.matrix[,2]==3);
+            # For vitellogenic population size, column 2 (Stage) must be 4.
+            Vitellogenic[row] = sum(vector.matrix[,2]==4);
+            # 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 +694,336 @@
 
         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));
+        if (plot_generations_separately) {
+            if (process_eggs) {
+                # For egg life stage of generation F1 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_nymphs) {
+                # For nymph life stage of generation F1 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.nymph[row] = sum((vector.matrix[,1]==0 & vector.matrix[,2]==1) | (vector.matrix[,1]==0 & vector.matrix[,2]==2));
+                # For 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.nymph[row] = sum((vector.matrix[,1]==1 & vector.matrix[,2]==1) | (vector.matrix[,1]==1 & vector.matrix[,2]==2));
+                # For 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.nymph[row] = sum((vector.matrix[,1]==2 & vector.matrix[,2]==1) | (vector.matrix[,1]==2 & vector.matrix[,2]==2));
+            }
+            if (process_adults) {
+                # For 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-vitellogenic)
+                # - column 1 (Generation) is 0 and column 2 (Stage) is 4 (Vitellogenic)
+                # - column 1 (Generation) is 0 and column 2 (Stage) is 5 (Diapausing)
+                P.adult[row] = sum((vector.matrix[,1]==0 & vector.matrix[,2]==3) | (vector.matrix[,1]==0 & vector.matrix[,2]==4) | (vector.matrix[,1]==0 & vector.matrix[,2]==5));
+                # For 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-vitellogenic)
+                # - column 1 (Generation) is 1 and column 2 (Stage) is 4 (Vitellogenic)
+                # - column 1 (Generation) is 1 and column 2 (Stage) is 5 (Diapausing)
+                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));
+                # For 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-vitellogenic)
+                # - column 1 (Generation) is 2 and column 2 (Stage) is 4 (Vitellogenic)
+                # - column 1 (Generation) is 2 and column 2 (Stage) is 5 (Diapausing)
+                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.
 
     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[,N.replications] = Eggs;
+    }
+    if (process_nymphs) {
+        YoungNymphs.replications[,N.replications] = YoungNymphs;
+        OldNymphs.replications[,N.replications] = OldNymphs;
+    }
+    if (process_adults) {
+        Previtellogenic.replications[,N.replications] = Previtellogenic;
+        Vitellogenic.replications[,N.replications] = Vitellogenic;
+        Diapausing.replications[,N.replications] = Diapausing;
+    }
     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;
-
+    if (plot_generations_separately) {
+        # P is Parental, or overwintered adults.
+        P.replications[,N.replications] = overwintering_adult.population;
+        # F1 is the first field-produced generation.
+        F1.replications[,N.replications] = first_generation.population;
+        # F2 is the second field-produced generation.
+        F2.replications[,N.replications] = second_generation.population;
+        if (process_eggs) {
+            P_eggs.replications[,N.replications] = P.egg;
+            F1_eggs.replications[,N.replications] = F1.egg;
+            F2_eggs.replications[,N.replications] = F2.egg;
+        }
+        if (process_nymphs) {
+            P_nymphs.replications[,N.replications] = P.nymph;
+            F1_nymphs.replications[,N.replications] = F1.nymph;
+            F2_nymphs.replications[,N.replications] = F2.nymph;
+        }
+        if (process_adults) {
+            P_adults.replications[,N.replications] = P.adult;
+            F1_adults.replications[,N.replications] = F1.adult;
+            F2_adults.replications[,N.replications] = F2.adult;
+        }
+    }
     population.replications[,N.replications] = total.population;
 }
 
-# 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);
-
-# 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);
+if (process_eggs) {
+    # 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);
+}
+if (process_nymphs) {
+    # Calculate nymph populations for selected life stage.
+    if (life_stage_nymph=="Total") {
+        # Mean value for all nymphs.
+        nymphs = apply((YoungNymphs.replications+OldNymphs.replications), 1, mean);
+        # Standard error for all nymphs.
+        nymphs.std_error = apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd);
+    } else if (life_stage_nymph=="Young") {
+        # Mean value for young nymphs.
+        nymphs = apply(YoungNymphs.replications, 1, mean);
+        # Standard error for young nymphs.
+        nymphs.std_error = apply(YoungNymphs.replications / sqrt(opt$replications), 1, sd);
+    } else if (life_stage_nymph=="Old") {
+        # Mean value for old nymphs.
+        nymphs = apply(OldNymphs.replications, 1, mean);
+        # Standard error for old nymphs.
+        nymphs.std_error = apply(OldNymphs.replications / sqrt(opt$replications), 1, sd);
+    }
+}
+if (process_adults) {
+    # Calculate adult populations for selected life stage.
+    if (life_stages_adult=="Total") {
+        # Mean value for all adults.
+        adults = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, mean);
+        # Standard error for all adults.
+        adults.std_error = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications);
+    } else if (life_stages_adult == "Pre-vittelogenic") {
+        # Mean value for previtellogenic adults.
+        adults = apply(Previtellogenic.replications, 1, mean);
+        # Standard error for previtellogenic adults.
+        adults.std_error = apply(Previtellogenic.replications, 1, sd) / sqrt(opt$replications);
+    } else if (life_stages_adult == "Vittelogenic") {
+        # Mean value for vitellogenic adults.
+        adults = apply(Vitellogenic.replications, 1, mean);
+        # Standard error for vitellogenic adults.
+        adults.std_error = apply(Vitellogenic.replications, 1, sd) / sqrt(opt$replications);
+    } else if (life_stages_adult == "Diapausing") {
+        # Mean value for vitellogenic adults.
+        adults = apply(Diapausing.replications, 1, mean);
+        # Standard error for vitellogenic adults.
+        adults.std_error = apply(Diapausing.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);
-
-# 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);
+if (plot_generations_separately) {
+    # Mean value for P which is Parental, or overwintered adults.
+    P = apply(P.replications, 1, mean);
+    # Standard error for P.
+    P.std_error = apply(P.replications, 1, sd) / sqrt(opt$replications);
+    # Mean value for F1, which is the first field-produced generation.
+    F1 = apply(F1.replications, 1, mean);
+    # Standard error for F1.
+    F1.std_error = apply(F1.replications, 1, sd) / sqrt(opt$replications);
+    # Mean value for F2, which is the second field-produced generation.
+    F2 = apply(F2.replications, 1, mean);
+    # Standard error for F2.
+    F2.std_error = apply(F2.replications, 1, sd) / sqrt(opt$replications);
+    if (process_eggs) {
+        # Mean value for P eggs.
+        P_eggs = apply(P_eggs.replications, 1, mean);
+        # Standard error for P_eggs.
+        P_eggs.std_error = apply(P_eggs.replications, 1, sd) / sqrt(opt$replications);
+        # Mean value for F1 eggs.
+        F1_eggs = apply(F1_eggs.replications, 1, mean);
+        # Standard error for F1 eggs.
+        F1_eggs.std_error = apply(F1_eggs.replications, 1, sd) / sqrt(opt$replications);
+        # Mean value for F2 eggs.
+        F2_eggs = apply(F2_eggs.replications, 1, mean);
+        # Standard error for F2 eggs.
+        F2_eggs.std_error = apply(F2_eggs.replications, 1, sd) / sqrt(opt$replications);
+    }
+    if (process_nymphs) {
+        # Mean value for P nymphs.
+        P_nymphs = apply(P_nymphs.replications, 1, mean);
+        # Standard error for P_nymphs.
+        P_nymphs.std_error = apply(P_nymphs.replications, 1, sd) / sqrt(opt$replications);
+        # Mean value for F1 nymphs.
+        F1_nymphs = apply(F1_nymphs.replications, 1, mean);
+        # Standard error for F1 nymphs.
+        F1_nymphs.std_error = apply(F1_nymphs.replications, 1, sd) / sqrt(opt$replications);
+        # Mean value for F2 nymphs.
+        F2_nymphs = apply(F2_nymphs.replications, 1, mean);
+        # Standard error for F2 eggs.
+        F2_nymphs.std_error = apply(F2_nymphs.replications, 1, sd) / sqrt(opt$replications);
+    }
+    if (process_adults) {
+        # Mean value for P adults.
+        P_adults = apply(P_adults.replications, 1, mean);
+        # Standard error for P_adults.
+        P_adults.std_error = apply(P_adults.replications, 1, sd) / sqrt(opt$replications);
+        # Mean value for F1 adults.
+        F1_adults = apply(F1_adults.replications, 1, mean);
+        # Standard error for F1 adults.
+        F1_adults.std_error = apply(F1_adults.replications, 1, sd) / sqrt(opt$replications);
+        # Mean value for F2 adults.
+        F2_adults = apply(F2_adults.replications, 1, mean);
+        # Standard error for F2 adults.
+        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.
+# Information needed for plots plots.
 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 = paste(output_dir, "egg_pop_by_generation.pdf", sep="/");
+            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(F2_eggs) + 100;
+            render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, 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") {
+            # Start PDF device driver.
+            dev.new(width=20, height=30);
+            file_name = paste(tolower(life_stage_nymph), "nymph_pop_by_generation.pdf", sep="_");
+            file_path = paste(output_dir, file_name, sep="/");
+            pdf(file=file_path, width=20, height=30, bg="white");
+            par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+            # Nymph population size by generation.
+            maxval = max(F2_nymphs) + 100;
+            render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
+                opt$replications, life_stage, group=P_nymphs, group_std_error=P_nymphs.std_error, group2=F1_nymphs, group2_std_error=F1_nymphs.std_error,
+                group3=F2_nymphs, group3_std_error=F2_nymphs.std_error, life_stage_nymph=life_stage_nymph);
+            # Turn off device driver to flush output.
+            dev.off();
+        } else if (life_stage == "Adult") {
+            # Start PDF device driver.
+            dev.new(width=20, height=30);
+            file_name = paste(tolower(life_stages_adult), "adult_pop_by_generation.pdf", sep="_");
+            file_path = paste(output_dir, file_name, sep="/");
+            pdf(file=file_path, width=20, height=30, bg="white");
+            par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+            # Adult population size by generation.
+            maxval = max(F2_adults) + 100;
+            render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
+                opt$replications, life_stage, group=P_adults, group_std_error=P_adults.std_error, group2=F1_adults, group2_std_error=F1_adults.std_error,
+                group3=F2_adults, group3_std_error=F2_adults.std_error, life_stages_adult=life_stages_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 = paste(output_dir, "total_pop_by_generation.pdf", sep="/");
+            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(F2);
+            render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, 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 = paste(output_dir, "egg_pop_size.pdf", sep="/");
+            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);
+            render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, 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") {
+            # Start PDF device driver.
+            dev.new(width=20, height=30);
+            file_name = paste(tolower(life_stage_nymph), "nymph_pop_size.pdf", sep="_");
+            file_path = paste(output_dir, file_name, sep="/");
+            pdf(file=file_path, width=20, height=30, bg="white");
+            par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+            # Nymph population size.
+            maxval = max(nymphs+nymphs.std_error);
+            render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
+                opt$replications, life_stage, group=nymphs, group_std_error=nymphs.std_error, life_stage_nymph=life_stage_nymph);
+            # Turn off device driver to flush output.
+            dev.off();
+        } else if (life_stage == "Adult") {
+            # Start PDF device driver.
+            dev.new(width=20, height=30);
+            file_name = paste(tolower(life_stages_adult), "adult_pop_size.pdf", sep="_");
+            file_path = paste(output_dir, file_name, sep="/");
+            pdf(file=file_path, width=20, height=30, bg="white");
+            par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
+            # Adult population size.
+            maxval = max(adults+adults.std_error);
+            render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
+                opt$replications, life_stage, group=adults, group_std_error=adults.std_error, life_stages_adult=life_stages_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 = paste(output_dir, "total_pop_size.pdf", sep="/");
+            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, nymphs+nymphs.std_error, adults+adults.std_error);
+            render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
+                opt$replications, life_stage, group=adults, group_std_error=adults.std_error, group2=nymphs, group2_std_error=nymphs.std_error, group3=eggs,
+                group3_std_error=eggs.std_error);
+            # Turn off device driver to flush output.
+            dev.off();
+        }
+    }
+}