--- a/insect_phenology_model.R	Tue Dec 05 15:15:49 2017 -0500
+++ b/insect_phenology_model.R	Tue Dec 05 15:16:14 2017 -0500
@@ -3,23 +3,24 @@
 suppressPackageStartupMessages(library("optparse"))
 
 option_list <- list(
-    make_option(c("-a", "--adult_mortality"), action="store", dest="adult_mortality", type="integer", help="Adjustment rate for adult mortality"),
-    make_option(c("-b", "--adult_accumulation"), action="store", dest="adult_accumulation", type="integer", help="Adjustment of degree-days accumulation (old nymph->adult)"),
-    make_option(c("-c", "--egg_mortality"), action="store", dest="egg_mortality", type="integer", help="Adjustment rate for egg mortality"),
-    make_option(c("-e", "--location"), action="store", dest="location", help="Selected location"),
-    make_option(c("-f", "--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"),
-    make_option(c("-i", "--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),
-    make_option(c("-j", "--nymph_mortality"), action="store", dest="nymph_mortality", type="integer", help="Adjustment rate for nymph mortality"),
-    make_option(c("-k", "--old_nymph_accumulation"), action="store", dest="old_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (young nymph->old nymph)"),
-    make_option(c("-n", "--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"),
-    make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"),
-    make_option(c("-p", "--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"),
-    make_option(c("-q", "--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"),
-    make_option(c("-s", "--replications"), action="store", dest="replications", type="integer", help="Number of replications"),
-    make_option(c("-t", "--std_error_plot"), action="store", dest="std_error_plot", help="Plot Standard error"),
-    make_option(c("-v", "--input"), action="store", dest="input", help="Temperature data for selected location"),
-    make_option(c("-y", "--young_nymph_accumulation"), action="store", dest="young_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)"),
-    make_option(c("-x", "--insect"), action="store", dest="insect", help="Insect name")
+    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("--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("--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("--young_nymph_accumulation"), action="store", dest="young_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)")
 )
 
 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
@@ -228,8 +229,7 @@
 }
 
 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.
+# All latitude values are the same, so get the value from the first row.
 latitude <- temperature_data_frame$LATITUDE[1]
 
 # Initialize matrices.
@@ -255,9 +255,8 @@
 
 # Process replications.
 for (N.replications in 1:opt$replications) {
-    # During each replication start with 1000 individuals.
-    # TODO: user definable as well?
-    num_insects <- 1000
+    # 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.
@@ -337,7 +336,6 @@
                 death.vector <- c(death.vector, i)
             }
             else {
-                # Aggregrate index of dead bug.
                 # End of diapause.
                 if (vector.individual[1] == 0 && vector.individual[2] == 3) {
                     # Overwintering adult (previttelogenic).
@@ -348,7 +346,7 @@
                         vector.matrix[i,] <- vector.individual
                     }
                     else {
-                        # Add to # Add average temperature for current day..
+                        # Add to # 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
@@ -372,7 +370,7 @@
                         vector.matrix[i,] <- vector.individual
                     }
                 }
-                # Event 2 oviposition -- where population dynamics comes from.
+                # Oviposition -- where population dynamics comes from.
                 if (vector.individual[2] == 4 && vector.individual[1] == 0 && mean.temp > 10) {
                     # Vittelogenic stage, overwintering generation.
                     if (vector.individual[4] == 0) {
@@ -404,7 +402,7 @@
                         birth.vector <- rbind(birth.vector, new.vector)
                     }
                 }
-                # Event 2 oviposition -- for generation 1.
+                # Oviposition -- for generation 1.
                 if (vector.individual[2] == 4 && vector.individual[1] == 1 && mean.temp > 12.5 && doy < 222) {
                     # Vittelogenic stage, 1st generation
                     if (vector.individual[4] == 0) {
@@ -436,10 +434,8 @@
                         birth.vector <- rbind(birth.vector, new.vector)
                     }
                 }
-                # Event 3 development (with diapause determination).
-                # Event 3.1 egg development to young nymph (vector.individual[2]=0 -> egg).
+                # Egg to young nymph.
                 if (vector.individual[2] == 0) {
-                    # Egg stage.
                     # Add average temperature for current day.
                     vector.individual[3] <- vector.individual[3] + averages.temp
                     if (vector.individual[3] >= (68+opt$young_nymph_accumulation)) {
@@ -454,9 +450,8 @@
                     }
                     vector.matrix[i,] <- vector.individual
                 }
-                # Event 3.2 young nymph to old nymph (vector.individual[2]=1 -> young nymph: determines diapause).
+                # Young nymph to old nymph.
                 if (vector.individual[2] == 1) {
-                    # Young nymph stage.
                     # Add average temperature for current day.
                     vector.individual[3] <- vector.individual[3] + averages.temp
                     if (vector.individual[3] >= (250+opt$old_nymph_accumulation)) {
@@ -474,16 +469,15 @@
                     }
                     vector.matrix[i,] <- vector.individual
                 }
-                # Event 3.3 old nymph to adult: previttelogenic or diapausing?
+                # Old nymph to adult: previttelogenic or diapausing?
                 if (vector.individual[2] == 2) {
-                    # Old nymph stage.
                     # Add average temperature for current day.
                     vector.individual[3] <- vector.individual[3] + averages.temp
                     if (vector.individual[3] >= (200+opt$adult_accumulation)) {
                         # From old to adult, degree_days requirement met.
                         current.gen <- vector.individual[1]
                         if (vector.individual[5] == 0) {
-                            # Non-diapausing adult -- previttelogenic.
+                            # Previttelogenic.
                             vector.individual <- c(current.gen, 3, 0, 0, 0)
                         }
                         else {
@@ -497,7 +491,7 @@
                     }
                     vector.matrix[i,] <- vector.individual
                 }
-                # Event 4 growing of diapausing adult (unimportant, but still necessary).
+                # Growing of diapausing adult (unimportant, but still necessary).
                 if (vector.individual[2] == 5) {
                     vector.individual[3] <- vector.individual[3] + averages.temp
                     vector.individual[4] <- vector.individual[4] + 1
@@ -571,10 +565,10 @@
     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
-    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
@@ -628,27 +622,26 @@
 # 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=40)
+dev.new(width=20, height=30)
 
 # Start PDF device driver to save charts to output.
-pdf(file=opt$output, width=20, height=40, bg="white")
+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.
+# 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
+# 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)
-# Subfigure 2: population size by generation
+# 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)
-# Subfigure 3: adult population size by generation
+# 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)