Mercurial > repos > greg > insect_phenology_model

changeset 99:51045fde125e draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Mon, 04 Dec 2017 13:04:11 -0500
parents d33a401c0153
children 52114847afff
files insect_phenology_model.R
diffstat 1 files changed, 49 insertions(+), 58 deletions(-) [+]
line wrap: on
line diff
--- a/insect_phenology_model.R	Mon Dec 04 11:54:10 2017 -0500
+++ b/insect_phenology_model.R	Mon Dec 04 13:04:11 2017 -0500
@@ -217,7 +217,7 @@
 g2.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
 g0a.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
 g1a.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
-g2a.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)
 
 # Loop through replications.
 for (N.replications in 1:opt$replications) {
@@ -467,13 +467,13 @@
                 }
             } # Else if it is still alive.
         } # End of the individual bug loop.
-        # Find how many died.
+        # Find the number of deaths.
         num_insects.death <- length(death.vector)
         if (num_insects.death > 0) {
+      # Remove record of dead.
             vector.matrix <- vector.matrix[-death.vector, ]
         }
-        # Remove record of dead.
-        # Find how many new born.
+        # Find the number of births.
         num_insects.newborn <- length(birth.vector[,1])
         vector.matrix <- rbind(vector.matrix, birth.vector)
         # Update population size for the next day.
@@ -481,38 +481,29 @@
 
         # Aggregate results by day.
         tot.pop <- c(tot.pop, num_insects)
-        # Egg.
-        sum_eggs <- sum(vector.matrix[,2] == 0)
-        # Young nymph.
-        sum_young_nymphs <- sum(vector.matrix[,2] == 1)
-        # Old nymph.
-        sum_old_nymphs <- sum(vector.matrix[,2] == 2)
-        # Previtellogenic.
-        sum_previtellogenic <- sum(vector.matrix[,2] == 3)
-        # Vitellogenic.
-        sum_vitellogenic <- sum(vector.matrix[,2] == 4)
-        # Diapausing.
-        sum_diapausing <- sum(vector.matrix[,2] == 5)
-        # Overwintering adult.
-        sum_overwintering_adults <- sum(vector.matrix[,1] == 0)
-        # First generation.
-        sum_first_generation <- sum(vector.matrix[,1] == 1)
-        # Second generation.
-        sum_second_generation <- sum(vector.matrix[,1] == 2)
-        # Sum of all adults.
+
+        # All adults population size.
         num_insects.adult <- sum(vector.matrix[,2] == 3) + sum(vector.matrix[,2] == 4) + sum(vector.matrix[,2] == 5)
 
-        # Population sizes.
-        gen0.pop[row] <- sum_overwintering_adults
-        gen1.pop[row] <- sum_first_generation
-        gen2.pop[row] <- sum_second_generation
+        # Overwintering adult population size.
+        gen0.pop[row] <- sum(vector.matrix[,1] == 0)
+        # First generation population size.
+        gen1.pop[row] <- sum(vector.matrix[,1] == 1)
+        # Second generation population size.
+        gen2.pop[row] <- sum(vector.matrix[,1] == 2)
 
-        S0[row] <- sum_eggs
-        S1[row] <- sum_young_nymphs
-        S2[row] <- sum_old_nymphs
-        S3[row] <- sum_previtellogenic
-        S4[row] <- sum_vitellogenic
-        S5[row] <- sum_diapausing
+        # Egg population size.
+        S0[row] <- sum(vector.matrix[,2] == 0)
+        # Young nymph population size.
+        S1[row] <- sum(vector.matrix[,2] == 1)
+        # Old nymph population size.
+        S2[row] <- sum(vector.matrix[,2] == 2)
+        # Previtellogenic population size.
+        S3[row] <- sum(vector.matrix[,2] == 3)
+        # Vitellogenic population size.
+        S4[row] <- sum(vector.matrix[,2] == 4)
+        # Diapausing population size.
+        S5[row] <- sum(vector.matrix[,2] == 5)
 
         g0.adult[row] <- sum(vector.matrix[,1] == 0)
         g1.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))
@@ -541,34 +532,34 @@
     g2.replications[,N.replications] <- gen2.pop
     g0a.replications[,N.replications] <- g0.adult
     g1a.replications[,N.replications] <- g1.adult
-    g2a.replications[,N.replications] <- g2.adult
+    F2_adults.replications[,N.replications] <- g2.adult
 }
 
 # Mean value for adults.
-adult_mean <- apply((S3.replications + S4.replications + S5.replications), 1, mean)
+adults <- apply((S3.replications + S4.replications + S5.replications), 1, mean)
 # Mean value for nymphs.
-nymph_mean <- apply((S1.replications + S2.replications), 1, mean)
+nymphs <- apply((S1.replications + S2.replications), 1, mean)
 # Mean value for eggs.
-egg_mean <- apply(S0.replications, 1, mean)
+eggs <- apply(S0.replications, 1, mean)
 # Mean value for P.
 g0 <- apply(g0.replications, 1, mean)
 # Mean value for F1.
 g1 <- apply(g1.replications, 1, mean)
 # Mean value for F2.
 g2 <- apply(g2.replications, 1, mean)
-# Mean value for P adult.
+# Mean value for P adults.
 g0a <- apply(g0a.replications, 1, mean)
-# Mean value for F1 adult.
+# Mean value for F1 adults.
 g1a <- apply(g1a.replications, 1, mean)
-# Mean value for F2 adult.
-F2_adult_mean <- apply(g2a.replications, 1, mean)
+# Mean value for F2 adults.
+F2_adults <- apply(F2_adults.replications, 1, mean)
 
 # Standard error for adults.
-adult_mean.std_error <- apply((S3.replications + S4.replications + S5.replications), 1, sd) / sqrt(opt$replications)
+adults.std_error <- apply((S3.replications + S4.replications + S5.replications), 1, sd) / sqrt(opt$replications)
 # Standard error for nymphs.
-nymph_mean.std_error <- apply((S1.replications + S2.replications) / sqrt(opt$replications), 1, sd)
+nymphs.std_error <- apply((S1.replications + S2.replications) / sqrt(opt$replications), 1, sd)
 # Standard error for eggs.
-egg_mean.std_error <- apply(S0.replications, 1, sd) / sqrt(opt$replications)
+eggs.std_error <- apply(S0.replications, 1, sd) / sqrt(opt$replications)
 # Standard error value for P.
 g0.std_error <- apply(g0.replications, 1, sd) / sqrt(opt$replications)
 # Standard error for F1.
@@ -580,7 +571,7 @@
 # Standard error for F1 adult.
 g1a.std_error <- apply(g1a.replications, 1, sd) / sqrt(opt$replications)
 # Standard error for F2 adult.
-g2a.std_error <- apply(g2a.replications, 1, sd) / sqrt(opt$replications)
+F2_adults.std_error <- apply(F2_adults.replications, 1, sd) / sqrt(opt$replications)
 
 dev.new(width=20, height=30)
 
@@ -597,25 +588,25 @@
 
 # Subfigure 1: population size by life stage
 title <- paste(opt$insect, ": Total pop. by life stage :", opt$location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ")
-plot(day.all, adult_mean, main=title, type="l", ylim=c(0, max(egg_mean + egg_mean.std_error, nymph_mean + nymph_mean.std_error, adult_mean + adult_mean.std_error)), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)
+plot(day.all, adults, main=title, type="l", ylim=c(0, max(eggs + eggs.std_error, nymphs + nymphs.std_error, adults + adults.std_error)), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)
 # Young and old nymphs.
-lines(day.all, nymph_mean, lwd=2, lty=1, col=2)
+lines(day.all, nymphs, lwd=2, lty=1, col=2)
 # Eggs
-lines(day.all, egg_mean, lwd=2, lty=1, col=4)
+lines(day.all, eggs, lwd=2, lty=1, col=4)
 axis(1, at=c(1:12) * 30 - 15, cex.axis=3, labels=c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
 axis(2, cex.axis=3)
 legend("topleft", c("Egg", "Nymph", "Adult"), lty=c(1, 1, 1), col=c(4, 2, 1), cex=3)
 if (opt$std_error_plot == 1) {
     # Add Standard error lines to plot
     # Standard error for adults
-    lines (day.all, adult_mean+adult_mean.std_error, lty=2)
-    lines (day.all, adult_mean-adult_mean.std_error, lty=2)
+    lines (day.all, adults+adults.std_error, lty=2)
+    lines (day.all, adults-adults.std_error, lty=2)
     # Standard error for nymphs
-    lines (day.all, nymph_mean+nymph_mean.std_error, col=2, lty=2)
-    lines (day.all, nymph_mean-nymph_mean.std_error, col=2, lty=2)
+    lines (day.all, nymphs+nymphs.std_error, col=2, lty=2)
+    lines (day.all, nymphs-nymphs.std_error, col=2, lty=2)
     # Standard error for eggs
-    lines (day.all, egg_mean+egg_mean.std_error, col=4, lty=2)
-    lines (day.all, egg_mean-egg_mean.std_error, col=4, lty=2)
+    lines (day.all, eggs+eggs.std_error, col=4, lty=2)
+    lines (day.all, eggs-eggs.std_error, col=4, lty=2)
 }
 
 # Subfigure 2: population size by generation
@@ -641,9 +632,9 @@
 
 # Subfigure 3: adult population size by generation
 title <- paste(opt$insect, ": Adult pop. by generation :", opt$location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ")
-plot(day.all, g0a, ylim=c(0, max(F2_adult_mean) + 100), main=title, type="l", axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)
+plot(day.all, g0a, ylim=c(0, max(F2_adults) + 100), main=title, type="l", axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)
 lines(day.all, g1a, lwd = 2, lty = 1, col=2)
-lines(day.all, F2_adult_mean, lwd = 2, lty = 1, col=4)
+lines(day.all, F2_adults, lwd = 2, lty = 1, col=4)
 axis(1, at=c(1:12) * 30 - 15, cex.axis=3, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
 axis(2, cex.axis=3)
 legend("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
@@ -656,8 +647,8 @@
     lines (day.all, g1a+g1a.std_error, col=2, lty=2)
     lines (day.all, g1a-g1a.std_error, col=2, lty=2)
     # Standard error for eggs
-    lines (day.all, F2_adult_mean+g2a.std_error, col=4, lty=2)
-    lines (day.all, F2_adult_mean-g2a.std_error, col=4, lty=2)
+    lines (day.all, F2_adults+F2_adults.std_error, col=4, lty=2)
+    lines (day.all, F2_adults-F2_adults.std_error, col=4, lty=2)
 }
 
 # Turn off device driver to flush output.