--- a/insect_phenology_model.R	Thu Nov 09 12:02:37 2017 -0500
+++ b/insect_phenology_model.R	Thu Nov 09 13:25:06 2017 -0500
@@ -25,16 +25,19 @@
 args <- parse_args(parser, positional_arguments=TRUE)
 opt <- args$options
 
-convert_csv_to_rdata=function(loc, temperature_data)
+convert_csv_to_rdata=function(data_matrix)
 {
-    expdata <- matrix(rep(0, opt$num_days * 6), nrow=opt$num_days)
-    expdata[,1] <- c(1:opt$num_days)
+    # Also stores the current year.
+    date <- data_matrix[1, 3]
+    year <- substr(date, 1, 4)
+    # Integer day of the year.
+    data_matrix[,1] <- c(1:opt$num_days)
     # Minimum
-    expdata[,2] <- temperature_data[c(1:opt$num_days), 5]
+    data_matrix[,2] <- temperature_data[c(1:opt$num_days), 5]
     # Maximum
-    expdata[,3] <- temperature_data[c(1:opt$num_days), 6]
+    data_matrix[,3] <- temperature_data[c(1:opt$num_days), 6]
     namedat <- "tempdata.Rdat"
-    save(expdata, file=namedat)
+    save(data_matrix, file=namedat)
     namedat
 }
 
@@ -59,8 +62,8 @@
     # 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
+    dnp <- data_matrix[date, 2]  # daily minimum
+    dxp <- data_matrix[date, 3]  # daily maximum
     dmean <- 0.5 * (dnp + dxp)
     dd <- 0  # initialize degree day accumulation
 
@@ -200,7 +203,8 @@
 
 # Read in the input temperature datafile into a Data Frame object.
 temperature_data <- read.csv(file=opt$input, header=T, sep=",")
-temperature_file_path <- convert_csv_to_rdata(opt$location, temperature_data)
+raw_data_matrix <- matrix(rep(0, opt$num_days * 6), nrow=opt$num_days)
+temperature_file_path <- convert_csv_to_rdata(raw_data_matrix)
 latitude <- temperature_data[1, 1]
 
 cat("Number of days: ", opt$num_days, "\n")
@@ -578,7 +582,6 @@
 g2a.se <- apply(g2a.rep, 1, sd) / sqrt(opt$replications)
 
 dev.new(width=20, height=20)
-sink("/dev/null")
 
 # Start PDF device driver to save charts to output.
 pdf(file=opt$output, height=20, width=20, bg="white")
@@ -586,30 +589,32 @@
 par(mar = c(5, 6, 4, 4), mfrow=c(3, 1))
 
 # Subfigure 2: population size by life stage
-plot(day.all, sa, main = "BSMB Total Population Size by Life Stage", type = "l", ylim = c(0, max(se + se.se, sn + sn.se, sa + sa.se)), axes = F, lwd = 2, xlab = "", ylab = "Number", cex = 2, cex.lab = 2, cex.axis = 2, cex.main = 2)
-# Young and old nymphs
-lines(day.all, sn, lwd = 2, lty = 1, col = 2)
+title <- paste(opt$location, year, "BSMB Total Population Size by Life Stage", sep=" ")
+plot(day.all, sa, main=title, type="l", ylim=c(0, max(se + se.se, sn + sn.se, sa + sa.se)), axes=F, lwd=2, xlab="", ylab="Number", cex=2, cex.lab=2, cex.axis=2, cex.main=2)
+# Young and old nymphs.
+lines(day.all, sn, lwd=2, lty=1, col=2)
 # Eggs
-lines(day.all, se, lwd = 2, lty = 1, col = 4)
-axis(1, at = c(1:12) * 30 - 15, cex.axis = 2, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
+lines(day.all, se, lwd=2, lty=1, col=4)
+axis(1, at = c(1:12) * 30 - 15, cex.axis=2, labels=c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
 axis(2, cex.axis = 2)
 leg.text <- c("Egg", "Nymph", "Adult")
-legend("topleft", leg.text, lty = c(1, 1, 1), col = c(4, 2, 1), cex = 2)
+legend("topleft", leg.text, lty=c(1, 1, 1), col=c(4, 2, 1), cex=2)
 if (opt$se_plot == 1) {
     # add SE lines to plot
     # SE for adults
-    lines (day.all, sa + sa.se, lty = 2)
-    lines (day.all, sa - sa.se, lty = 2) 
+    lines (day.all, sa + sa.se, lty=2)
+    lines (day.all, sa - sa.se, lty=2) 
     # SE for nymphs
-    lines (day.all, sn + sn.se, col = 2, lty = 2)
-    lines (day.all, sn - sn.se, col = 2, lty = 2) 
+    lines (day.all, sn + sn.se, col=2, lty=2)
+    lines (day.all, sn - sn.se, col=2, lty=2) 
     # SE for eggs
-    lines (day.all, se + se.se, col = 4, lty = 2)
-    lines (day.all, se - se.se, col = 4, lty = 2) 
+    lines (day.all, se + se.se, col=4, lty=2)
+    lines (day.all, se - se.se, col=4, lty=2) 
 }
 
 # Subfigure 3: population size by generation
-plot(day.all, g0, main = "BSMB Total Population Size by Generation", type = "l", ylim = c(0, max(g2)), axes = F, lwd = 2, xlab = "", ylab = "Number", cex = 2, cex.lab = 2, cex.axis = 2, cex.main = 2)
+title <- paste(opt$location, year, "BSMB Total Population Size by Generation", sep=" ")
+plot(day.all, g0, main=title, type="l", ylim=c(0, max(g2)), axes=F, lwd=2, xlab="", ylab="Number", cex=2, cex.lab=2, cex.axis=2, cex.main=2)
 lines(day.all, g1, lwd = 2, lty = 1, col = 2)
 lines(day.all, g2, lwd = 2, lty = 1, col = 4)
 axis(1, at = c(1:12) * 30 - 15, cex.axis = 2, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
@@ -630,7 +635,8 @@
 }
 
 # Subfigure 4: adult population size by generation
-plot(day.all, g0a, ylim = c(0, max(g2a) + 100), main = "BSMB Adult Population Size by Generation", type = "l", axes = F, lwd = 2, xlab = "Year", ylab = "Number", cex = 2, cex.lab = 2, cex.axis = 2, cex.main = 2)
+title <- paste(opt$location, year, "BSMB Adult Population Size by Generation", sep=" ")
+plot(day.all, g0a, ylim=c(0, max(g2a) + 100), main=title, type="l", axes=F, lwd=2, xlab="Year", ylab="Number", cex=2, cex.lab=2, cex.axis=2, cex.main=2)
 lines(day.all, g1a, lwd = 2, lty = 1, col = 2)
 lines(day.all, g2a, lwd = 2, lty = 1, col = 4)
 axis(1, at = c(1:12) * 30 - 15, cex.axis = 2, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
@@ -650,6 +656,5 @@
     lines (day.all, g2a - g2a.se, col = 4, lty = 2) 
 }
 
-sink()
 # Turn off device driver to flush output.
 dev.off()