--- a/insect_phenology_model.R	Thu Nov 09 11:30:20 2017 -0500
+++ b/insect_phenology_model.R	Thu Nov 09 11:30:27 2017 -0500
@@ -25,24 +25,22 @@
 args <- parse_args(parser, positional_arguments=TRUE)
 opt <- args$options
 
-get_temperature_file_path=function(loc, temperature_data)
+convert_csv_to_rdata=function(loc, temperature_data)
 {
-    expdata <- matrix(rep(0, opt$num_days * 3), nrow=opt$num_days)
+    expdata <- matrix(rep(0, opt$num_days * 6), nrow=opt$num_days)
     expdata[,1] <- c(1:opt$num_days)
     # Minimum
-    expdata[,2] <- temperature_data[c(1:opt$num_days), 3]
+    expdata[,2] <- temperature_data[c(1:opt$num_days), 5]
     # Maximum
-    expdata[,3] <- temperature_data[c(1:opt$num_days), 2]
-    date <- temperature_data[1, 3]
-    year <- substr(date, 1, 4)
-    namedat <- paste(loc,  year, ".Rdat", sep="")
+    expdata[,3] <- temperature_data[c(1:opt$num_days), 6]
+    namedat <- "tempdata.Rdat"
     save(expdata, file=namedat)
     namedat
 }
 
 daylength=function(latitude, num_days)
 {
-    # from Forsythe 1995
+    # From Forsythe 1995.
     p=0.8333
     dl <- NULL
     for (i in 1:num_days) {
@@ -50,15 +48,16 @@
         phi <- asin(0.39795 * cos(theta))
         dl[i] <- 24 - 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi)))
     }
-    # return a vector of daylength for the number of
-    # days specifie din the input temperature data
+    # Return a vector of daylength for the number of
+    # days specified in the input temperature data.
     dl
 }
 
 hourtemp=function(latitude, date, temperature_file_path, num_days)
 {
     load(temperature_file_path)
-    # base development threshold for BMSB
+    # Base development threshold for Brown Marmolated Stink Bug
+    # insect phenology model.
     threshold <- 14.17
     dnp <- expdata[date, 2]  # daily minimum
     dxp <- expdata[date, 3]  # daily maximum
@@ -69,27 +68,29 @@
         dd <- 0
     }
     else {
-        # extract daylength data for the number of
-        # days specified in the input temperature data
+        # Extract daylength data for the number of
+        # days specified in the input temperature data.
         dlprofile <- daylength(latitude, num_days)
-        T <- NULL  # initialize hourly temperature
-        dh <- NULL #initialize degree hour vector
-        # calculate daylength in given date
+        # Initialize hourly temperature.
+        T <- NULL
+        # Initialize degree hour vector.
+        dh <- NULL
+        # Calculate daylength in given date.
         y <- dlprofile[date]
-        # night length
+        # Night length.
         z <- 24 - y
-        # lag coefficient
+        # Lag coefficient.
         a <- 1.86
-        # night coefficient
+        # Night coefficient.
         b <- 2.20
-        # sunrise time
+        # Sunrise time.
         risetime <- 12 - y / 2
-        # sunset time
+        # Sunset time.
         settime <- 12 + y / 2
         ts <- (dxp - dnp) * sin(pi * (settime - 5) / (y + 2 * a)) + dnp
         for (i in 1:24) {
             if (i > risetime && i<settime) {
-                # number of hours after Tmin until sunset
+                # Number of hours after Tmin until sunset.
                 m <- i - 5
                 T[i]=(dxp - dnp) * sin(pi * m / (y + 2 * a)) + dnp
                 if (T[i]<8.4) {
@@ -197,44 +198,45 @@
     return
 }
 
-cat("Replications: ", opt$replications, "\n")
-cat("Photoperiod: ", opt$photoperiod, "\n")
-cat("Oviposition rate: ", opt$oviposition, "\n")
-cat("Egg mortality rate: ", opt$egg_mort, "\n")
-cat("Nymph mortality rate: ", opt$nymph_mort, "\n")
-cat("Adult mortality rate: ", opt$adult_mort, "\n")
-cat("Min clutch size: ", opt$min_clutch_size, "\n")
-cat("Max clutch size: ", opt$max_clutch_size, "\n")
-cat("(egg->young nymph): ", opt$young_nymph_accum, "\n")
-cat("(young nymph->old nymph): ", opt$old_nymph_accum, "\n")
-cat("(old nymph->adult): ", opt$adult_accum, "\n")
-
 # Read in the input temperature datafile into a Data Frame object.
 temperature_data <- read.csv(file=opt$input, header=T, sep=",")
-temperature_file_path <- get_temperature_file_path(opt$location, temperature_data)
+temperature_file_path <- convert_csv_to_rdata(opt$location, temperature_data)
 latitude <- temperature_data[1, 1]
 
+cat("Number of days: ", opt$num_days, ", ")
+cat("Latitude: ", latitude, "\n")
+cat("Replications: ", opt$replications, ", ")
+cat("Photoperiod: ", opt$photoperiod, "\n")
+cat("Oviposition rate: ", opt$oviposition, ", ")
+cat("Egg mortality rate: ", opt$egg_mort, "\n")
+cat("Nymph mortality rate: ", opt$nymph_mort, ", ")
+cat("Adult mortality rate: ", opt$adult_mort, "\n")
+cat("Min clutch size: ", opt$min_clutch_size, ", ")
+cat("Max clutch size: ", opt$max_clutch_size, "\n")
+cat("(egg->young nymph): ", opt$young_nymph_accum, ", ")
+cat("(young nymph->old nymph): ", opt$old_nymph_accum, "\n")
+cat("(old nymph->adult): ", opt$adult_accum
+
 # Initialize matrix for results from all replications
 S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol = opt$replications)
 newborn.rep <- death.rep <- adult.rep <- pop.rep <- g0.rep <- g1.rep <- g2.rep <- g0a.rep <- g1a.rep <- g2a.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol=opt$replications)
 
 # loop through replications
 for (N.rep in 1:opt$replications) {
-    # during each replication
-    # start with 1000 individuals -- user definable as well?
+    # During each replication start with 1000 individuals.
+    # TODO: user definable as well?
     n <- 1000
-    # Generation, Stage, DD, T, Diapause
+    # Generation, Stage, DD, T, Diapause.
     vec.ini <- c(0, 3, 0, 0, 0)
-    # overwintering, previttelogenic, DD=0, T=0, no-diapause
+    # Overwintering, previttelogenic, DD=0, T=0, no-diapause.
     vec.mat <- rep(vec.ini, n)
-    # complete matrix for the population
-    vec.mat <- t(matrix(vec.mat, nrow=5))
-    # complete photoperiod profile in a year, requires daylength function
+    # Complete matrix for the population.
+    vec.mat <- base::t(matrix(vec.mat, nrow=5))
+    # Complete photoperiod profile in a year, requires daylength function.
     ph.p <- daylength(latitude, opt$num_days)
 
-    # time series of population size
+    # Time series of population size.
     tot.pop <- NULL
-    # gen.0 pop size
     gen0.pop <- rep(0, opt$num_days)
     gen1.pop <- rep(0, opt$num_days)
     gen2.pop <- rep(0, opt$num_days)
@@ -243,28 +245,25 @@
     N.newborn <- N.death <- N.adult <- rep(0, opt$num_days)
     dd.day <- rep(0, opt$num_days)
 
-    # start tick
-    ptm <- proc.time()
-
     # All the days included in the input temperature dataset.
     for (day in 1:opt$num_days) {
-        # photoperiod in the day
+        # Photoperiod in the day.
         photoperiod <- ph.p[day]
         temp.profile <- hourtemp(latitude, day, temperature_file_path, opt$num_days)
         mean.temp <- temp.profile[1]
         dd.temp <- temp.profile[2]
         dd.day[day] <- dd.temp
-        # trash bin for death
+        # Trash bin for death.
         death.vec <- NULL
-        # new born
+        # Newborn.
         birth.vec <- NULL
 
-        # all individuals
+        # All individuals.
         for (i in 1:n) {
-            # find individual record
+            # Find individual record.
             vec.ind <- vec.mat[i,]
-            # first of all, still alive?  
-            # adjustment for late season mortality rate
+            # First of all, still alive?
+            # Adjustment for late season mortality rate.
             if (latitude < 40.0) {
                 post.mort <- 1
                 day.kill <- 300
@@ -274,19 +273,19 @@
                 day.kill <- 250
             }
             if (vec.ind[2] == 0) {
-                # egg
+                # Egg.
                 death.prob = opt$egg_mort * mortality.egg(mean.temp)
             }
             else if (vec.ind[2] == 1 | vec.ind[2] == 2) {
                 death.prob = opt$nymph_mort * mortality.nymph(mean.temp)
             }
             else if (vec.ind[2] == 3 | vec.ind[2] == 4 | vec.ind[2] == 5) {
-                # for adult
+                # For adult.
                 if (day < day.kill) {
                     death.prob = opt$adult_mort * mortality.adult(mean.temp)
                 }
                 else {
-                    # increase adult mortality after fall equinox
+                    # Increase adult mortality after fall equinox.
                     death.prob = opt$adult_mort * post.mort * mortality.adult(mean.temp)
                 }
             }
@@ -296,218 +295,216 @@
                 death.vec <- c(death.vec, i)
             } 
             else {
-                # aggregrate index of dead bug
-                # event 1 end of diapause
+                # Aggregrate index of dead bug.
+                # Event 1 end of diapause.
                 if (vec.ind[1] == 0 && vec.ind[2] == 3) {
-                    # overwintering adult (previttelogenic)
+                    # Overwintering adult (previttelogenic).
                     if (photoperiod > opt$photoperiod && vec.ind[3] > 68 && day < 180) {
-                        # add 68C to become fully reproductively matured
-                        # transfer to vittelogenic
+                        # Add 68C to become fully reproductively matured.
+                        # Transfer to vittelogenic.
                         vec.ind <- c(0, 4, 0, 0, 0)
                         vec.mat[i,] <- vec.ind
                     }
                     else {
-                        # add to DD
+                        # Add to dd.
                         vec.ind[3] <- vec.ind[3] + dd.temp
-                        # add 1 day in current stage
+                        # Add 1 day in current stage.
                         vec.ind[4] <- vec.ind[4] + 1
                         vec.mat[i,] <- vec.ind
                     }
                 }
                 if (vec.ind[1] != 0 && vec.ind[2] == 3) {
-                    # NOT overwintering adult (previttelogenic)
+                    # Not overwintering adult (previttelogenic).
                     current.gen <- vec.ind[1]
                     if (vec.ind[3] > 68) {
-                        # add 68C to become fully reproductively matured
-                        # transfer to vittelogenic
+                        # Add 68C to become fully reproductively matured.
+                        # Transfer to vittelogenic.
                         vec.ind <- c(current.gen, 4, 0, 0, 0)
                         vec.mat[i,] <- vec.ind
                     }
                     else {
-                        # add to DD
+                        # Add to dd.
                         vec.ind[3] <- vec.ind[3] + dd.temp
-                        # add 1 day in current stage
+                        # Add 1 day in current stage.
                         vec.ind[4] <- vec.ind[4] + 1
                         vec.mat[i,] <- vec.ind
                     }
                 }
 
-                # event 2 oviposition -- where population dynamics comes from
+                # Event 2 oviposition -- where population dynamics comes from.
                 if (vec.ind[2] == 4 && vec.ind[1] == 0 && mean.temp > 10) {
-                    # vittelogenic stage, overwintering generation
+                    # Vittelogenic stage, overwintering generation.
                     if (vec.ind[4] == 0) {
-                        # just turned in vittelogenic stage
+                        # Just turned in vittelogenic stage.
                         n.birth=round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size))
                     }
                     else {
-                        # daily probability of birth
+                        # Daily probability of birth.
                         p.birth = opt$oviposition * 0.01
                         u1 <- runif(1)
                         if (u1 < p.birth) {
                             n.birth=round(runif(1, 2, 8))
                         }
                     }
-                    # add to DD
+                    # Add to dd.
                     vec.ind[3] <- vec.ind[3] + dd.temp
-                    # add 1 day in current stage
+                    # Add 1 day in current stage.
                     vec.ind[4] <- vec.ind[4] + 1
                     vec.mat[i,] <- vec.ind
                     if (n.birth > 0) {
-                        # add new birth -- might be in different generations
-                        # generation + 1
+                        # Add new birth -- might be in different generations.
                         new.gen <- vec.ind[1] + 1
-                        # egg profile
+                        # Egg profile.
                         new.ind <- c(new.gen, 0, 0, 0, 0)
                         new.vec <- rep(new.ind, n.birth)
-                        # update batch of egg profile
+                        # Update batch of egg profile.
                         new.vec <- t(matrix(new.vec, nrow=5))
-                        # group with total eggs laid in that day
+                        # Group with total eggs laid in that day.
                         birth.vec <- rbind(birth.vec, new.vec)
                     }
                 }
 
-                # event 2 oviposition -- for gen 1.
+                # Event 2 oviposition -- for gen 1.
                 if (vec.ind[2] == 4 && vec.ind[1] == 1 && mean.temp > 12.5 && day < 222) {
-                    # vittelogenic stage, 1st generation
+                    # Vittelogenic stage, 1st generation
                     if (vec.ind[4] == 0) {
-                        # just turned in vittelogenic stage
+                        # Just turned in vittelogenic stage.
                         n.birth=round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size))
                     }
                     else {
-                        # daily probability of birth
+                        # Daily probability of birth.
                         p.birth = opt$oviposition * 0.01
                         u1 <- runif(1)
                         if (u1 < p.birth) {
                             n.birth = round(runif(1, 2, 8))
                         }
                     }
-                    # add to DD
+                    # Add to dd.
                     vec.ind[3] <- vec.ind[3] + dd.temp
-                    # add 1 day in current stage
+                    # Add 1 day in current stage.
                     vec.ind[4] <- vec.ind[4] + 1
                     vec.mat[i,] <- vec.ind
                     if (n.birth > 0) {
-                        # add new birth -- might be in different generations
-                        # generation + 1
+                        # Add new birth -- might be in different generations.
                         new.gen <- vec.ind[1] + 1
-                        # egg profile
+                        # Egg profile.
                         new.ind <- c(new.gen, 0, 0, 0, 0)
                         new.vec <- rep(new.ind, n.birth)
-                        # update batch of egg profile
+                        # Update batch of egg profile.
                         new.vec <- t(matrix(new.vec, nrow=5))
-                        # group with total eggs laid in that day
+                        # Group with total eggs laid in that day.
                         birth.vec <- rbind(birth.vec, new.vec)
                     }
                 }
 
-                # event 3 development (with diapause determination)
-                # event 3.1 egg development to young nymph (vec.ind[2]=0 -> egg)
+                # Event 3 development (with diapause determination).
+                # Event 3.1 egg development to young nymph (vec.ind[2]=0 -> egg).
                 if (vec.ind[2] == 0) {
-                    # egg stage
-                    # add to DD
+                    # Egg stage.
+                    # Add to dd.
                     vec.ind[3] <- vec.ind[3] + dd.temp
                     if (vec.ind[3] >= (68 + opt$young_nymph_accum)) {
-                        # from egg to young nymph, DD requirement met
+                        # From egg to young nymph, DD requirement met.
                         current.gen <- vec.ind[1]
-                        # transfer to young nym stage
+                        # Transfer to young nymph stage.
                         vec.ind <- c(current.gen, 1, 0, 0, 0)
                     }
                     else {
-                        # add 1 day in current stage
+                        # Add 1 day in current stage.
                         vec.ind[4] <- vec.ind[4] + 1
                     }
                     vec.mat[i,] <- vec.ind
                 }
 
-                # event 3.2 young nymph to old nymph (vec.ind[2]=1 -> young nymph: determines diapause)
+                # Event 3.2 young nymph to old nymph (vec.ind[2]=1 -> young nymph: determines diapause).
                 if (vec.ind[2] == 1) {
-                    # young nymph stage
-                    # add to DD
+                    # young nymph stage.
+                    # add to dd.
                     vec.ind[3] <- vec.ind[3] + dd.temp
                     if (vec.ind[3] >= (250 + opt$old_nymph_accum)) {
-                        # from young to old nymph, DD requirement met
+                        # From young to old nymph, dd requirement met.
                         current.gen <- vec.ind[1]
-                        # transfer to old nym stage
+                        # Transfer to old nym stage.
                         vec.ind <- c(current.gen, 2, 0, 0, 0)
                         if (photoperiod < opt$photoperiod && day > 180) {
                             vec.ind[5] <- 1
-                        } # prepare for diapausing
+                        } # Prepare for diapausing.
                     }
                     else {
-                        # add 1 day in current stage
+                        # Add 1 day in current stage.
                         vec.ind[4] <- vec.ind[4] + 1
                     }
                     vec.mat[i,] <- vec.ind
                 }  
 
-                # event 3.3 old nymph to adult: previttelogenic or diapausing?
+                # Event 3.3 old nymph to adult: previttelogenic or diapausing?
                 if (vec.ind[2] == 2) {
-                    # old nymph stage
-                    # add to DD
+                    # Old nymph stage.
+                    # add to dd.
                     vec.ind[3] <- vec.ind[3] + dd.temp
                     if (vec.ind[3] >= (200 + opt$adult_accum)) {
-                        # from old to adult, DD requirement met
+                        # From old to adult, dd requirement met.
                         current.gen <- vec.ind[1]
                         if (vec.ind[5] == 0) {
-                            # non-diapausing adult -- previttelogenic
+                            # Non-diapausing adult -- previttelogenic.
                             vec.ind <- c(current.gen, 3, 0, 0, 0)
                         }
                         else {
-                            # diapausing 
+                            # Diapausing.
                             vec.ind <- c(current.gen, 5, 0, 0, 1)
                         }
                     }
                     else {
-                        # add 1 day in current stage
+                        # Add 1 day in current stage.
                         vec.ind[4] <- vec.ind[4] + 1
                     }
                     vec.mat[i,] <- vec.ind
                 }
 
-                # event 4 growing of diapausing adult (unimportant, but still necessary)## 
+                # Event 4 growing of diapausing adult (unimportant, but still necessary).
                 if (vec.ind[2] == 5) {
                     vec.ind[3] <- vec.ind[3] + dd.temp
                     vec.ind[4] <- vec.ind[4] + 1
                     vec.mat[i,] <- vec.ind
                 }
-            } # else if it is still alive
-        } # end of the individual bug loop
+            } # Else if it is still alive.
+        } # End of the individual bug loop.
 
-        # find how many died
+        # Find how many died.
         n.death <- length(death.vec)
         if (n.death > 0) {
             vec.mat <- vec.mat[-death.vec, ]
         }
-        # remove record of dead
-        # find how many new born  
+        # Remove record of dead.
+        # Find how many new born.
         n.newborn <- length(birth.vec[,1])
         vec.mat <- rbind(vec.mat, birth.vec)
-        # update population size for the next day
+        # Update population size for the next day.
         n <- n - n.death + n.newborn 
 
-        # aggregate results by day
+        # Aggregate results by day.
         tot.pop <- c(tot.pop, n) 
-        # egg
+        # Egg.
         s0 <- sum(vec.mat[,2] == 0)
-        # young nymph
+        # Young nymph.
         s1 <- sum(vec.mat[,2] == 1)
-        # old nymph
+        # Old nymph.
         s2 <- sum(vec.mat[,2] == 2)
-        # previtellogenic
+        # Previtellogenic.
         s3 <- sum(vec.mat[,2] == 3)
-        # vitellogenic
+        # Vitellogenic.
         s4 <- sum(vec.mat[,2] == 4)
-        # diapausing
+        # Diapausing.
         s5 <- sum(vec.mat[,2] == 5)
-        # overwintering adult
+        # Overwintering adult.
         gen0 <- sum(vec.mat[,1] == 0)
-        # first generation
+        # First generation.
         gen1 <- sum(vec.mat[,1] == 1)
-        # second generation
+        # Second generation.
         gen2 <- sum(vec.mat[,1] == 2)
-        # sum of all adults
+        # Sum of all adults.
         n.adult <- sum(vec.mat[,2] == 3) + sum(vec.mat[,2] == 4) + sum(vec.mat[,2] == 5)
-        # gen.0 pop size
+        # Gen eration 0 pop size.
         gen0.pop[day] <- gen0
         gen1.pop[day] <- gen1
         gen2.pop[day] <- gen2
@@ -527,7 +524,7 @@
     }   # end of days specified in the input temperature data
 
     dd.cum <- cumsum(dd.day)
-    # collect all the outputs
+    # Collect all the outputs.
     S0.rep[,N.rep] <- S0
     S1.rep[,N.rep] <- S1
     S2.rep[,N.rep] <- S2