Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 102:ce996e90f5a7 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Tue, 05 Dec 2017 10:32:04 -0500
parents 691b677d8f83
children 4d9e40d8af0f
comparison
equal deleted inserted replaced
101:691b677d8f83 102:ce996e90f5a7
1 #!/usr/bin/env Rscript 1 #!/usr/bin/env Rscript
2 2
3 suppressPackageStartupMessages(library("optparse")) 3 suppressPackageStartupMessages(library("optparse"))
4 4
5 option_list <- list( 5 option_list <- list(
6 make_option(c("-a", "--adult_mort"), action="store", dest="adult_mort", type="integer", help="Adjustment rate for adult mortality"), 6 make_option(c("-a", "--adult_mortality"), action="store", dest="adult_mortality", type="integer", help="Adjustment rate for adult mortality"),
7 make_option(c("-b", "--adult_accum"), action="store", dest="adult_accum", type="integer", help="Adjustment of degree-days accumulation (old nymph->adult)"), 7 make_option(c("-b", "--adult_accumulation"), action="store", dest="adult_accumulation", type="integer", help="Adjustment of degree-days accumulation (old nymph->adult)"),
8 make_option(c("-c", "--egg_mort"), action="store", dest="egg_mort", type="integer", help="Adjustment rate for egg mortality"), 8 make_option(c("-c", "--egg_mortality"), action="store", dest="egg_mortality", type="integer", help="Adjustment rate for egg mortality"),
9 make_option(c("-e", "--location"), action="store", dest="location", help="Selected location"), 9 make_option(c("-e", "--location"), action="store", dest="location", help="Selected location"),
10 make_option(c("-f", "--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"), 10 make_option(c("-f", "--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"),
11 make_option(c("-i", "--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"), 11 make_option(c("-i", "--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),
12 make_option(c("-j", "--nymph_mort"), action="store", dest="nymph_mort", type="integer", help="Adjustment rate for nymph mortality"), 12 make_option(c("-j", "--nymph_mortality"), action="store", dest="nymph_mortality", type="integer", help="Adjustment rate for nymph mortality"),
13 make_option(c("-k", "--old_nymph_accum"), action="store", dest="old_nymph_accum", type="integer", help="Adjustment of degree-days accumulation (young nymph->old nymph)"), 13 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)"),
14 make_option(c("-n", "--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"), 14 make_option(c("-n", "--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"),
15 make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"), 15 make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"),
16 make_option(c("-p", "--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"), 16 make_option(c("-p", "--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"),
17 make_option(c("-q", "--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"), 17 make_option(c("-q", "--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"),
18 make_option(c("-s", "--replications"), action="store", dest="replications", type="integer", help="Number of replications"), 18 make_option(c("-s", "--replications"), action="store", dest="replications", type="integer", help="Number of replications"),
19 make_option(c("-t", "--std_error_plot"), action="store", dest="std_error_plot", help="Plot Standard error"), 19 make_option(c("-t", "--std_error_plot"), action="store", dest="std_error_plot", help="Plot Standard error"),
20 make_option(c("-v", "--input"), action="store", dest="input", help="Temperature data for selected location"), 20 make_option(c("-v", "--input"), action="store", dest="input", help="Temperature data for selected location"),
21 make_option(c("-y", "--young_nymph_accum"), action="store", dest="young_nymph_accum", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)"), 21 make_option(c("-y", "--young_nymph_accumulation"), action="store", dest="young_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)"),
22 make_option(c("-x", "--insect"), action="store", dest="insect", help="Insect name") 22 make_option(c("-x", "--insect"), action="store", dest="insect", help="Insect name")
23 ) 23 )
24 24
25 parser <- OptionParser(usage="%prog [options] file", option_list=option_list) 25 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
26 args <- parse_args(parser, positional_arguments=TRUE) 26 args <- parse_args(parser, positional_arguments=TRUE)
27 opt <- args$options 27 opt <- args$options
28
29 parse_input_data = function(input_file, num_rows) {
30 # Read in the input temperature datafile into a data frame.
31 temperature_data_frame <- read.csv(file=input_file, header=T, strip.white=TRUE, sep=",")
32 num_columns <- dim(temperature_data_frame)[2]
33 if (num_columns == 6) {
34 # The input data has the following 6 columns:
35 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
36 # Set the column names for access when adding daylight length..
37 colnames(temperature_data_frame) <- c("LATITUDE","LONGITUDE", "DATE", "DOY", "TMIN", "TMAX")
38 # Add a column containing the daylight length for each day.
39 temperature_data_frame <- add_daylight_length(temperature_data_frame, num_columns, num_rows)
40 # Reset the column names with the additional column for later access.
41 colnames(temperature_data_frame) <- c("LATITUDE","LONGITUDE", "DATE", "DOY", "TMIN", "TMAX", "DAYLEN")
42 }
43 return(temperature_data_frame)
44 }
45 28
46 add_daylight_length = function(temperature_data_frame, num_columns, num_rows) { 29 add_daylight_length = function(temperature_data_frame, num_columns, num_rows) {
47 # Return a vector of daylight length (photoperido profile) for 30 # Return a vector of daylight length (photoperido profile) for
48 # the number of days specified in the input temperature data 31 # the number of days specified in the input temperature data
49 # (from Forsythe 1995). 32 # (from Forsythe 1995).
55 # of the temperature data for computation. 38 # of the temperature data for computation.
56 doy <- temperature_data_frame$DOY[i] 39 doy <- temperature_data_frame$DOY[i]
57 theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (doy - 186))) 40 theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (doy - 186)))
58 phi <- asin(0.39795 * cos(theta)) 41 phi <- asin(0.39795 * cos(theta))
59 # Compute the length of daylight for the day of the year. 42 # Compute the length of daylight for the day of the year.
60 daylight_length_vector[i] <- 24 - (24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi)))) 43 darkness_length <- 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi)))
44 daylight_length_vector[i] <- 24 - darkness_length
61 } 45 }
62 # Append daylight_length_vector as a new column to temperature_data_frame. 46 # Append daylight_length_vector as a new column to temperature_data_frame.
63 temperature_data_frame[, num_columns+1] <- daylight_length_vector 47 temperature_data_frame[, num_columns+1] <- daylight_length_vector
64 return(temperature_data_frame) 48 return(temperature_data_frame)
65 } 49 }
66 50
51 dev.egg = function(temperature) {
52 dev.rate = -0.9843 * temperature + 33.438
53 return(dev.rate)
54 }
55
56 dev.emerg = function(temperature) {
57 emerg.rate <- -0.5332 * temperature + 24.147
58 return(emerg.rate)
59 }
60
61 dev.old = function(temperature) {
62 n34 <- -0.6119 * temperature + 17.602
63 n45 <- -0.4408 * temperature + 19.036
64 dev.rate = mean(n34 + n45)
65 return(dev.rate)
66 }
67
68 dev.young = function(temperature) {
69 n12 <- -0.3728 * temperature + 14.68
70 n23 <- -0.6119 * temperature + 25.249
71 dev.rate = mean(n12 + n23)
72 return(dev.rate)
73 }
74
67 get_temperature_at_hour = function(latitude, temperature_data_frame, row, num_days) { 75 get_temperature_at_hour = function(latitude, temperature_data_frame, row, num_days) {
68 # Base development threshold for Brown Marmolated Stink Bug 76 # Base development threshold for Brown Marmolated Stink Bug
69 # insect phenology model. 77 # insect phenology model.
70 # TODO: Pass insect on the command line to accomodate more
71 # the just the Brown Marmolated Stink Bub.
72 threshold <- 14.17 78 threshold <- 14.17
73 # Minimum temperature for current row. 79 # Minimum temperature for current row.
74 dnp <- temperature_data_frame$TMIN[row] 80 curr_min_temp <- temperature_data_frame$TMIN[row]
75 # Maximum temperature for current row. 81 # Maximum temperature for current row.
76 dxp <- temperature_data_frame$TMAX[row] 82 curr_max_temp <- temperature_data_frame$TMAX[row]
77 # Mean temperature for current row. 83 # Mean temperature for current row.
78 dmean <- 0.5 * (dnp + dxp) 84 curr_mean_temp <- 0.5 * (curr_min_temp + curr_max_temp)
79 # Initialize degree day accumulation 85 # Initialize degree day accumulation
80 degree_days <- 0 86 degree_days <- 0
81 if (dxp < threshold) { 87 if (curr_max_temp < threshold) {
82 degree_days <- 0 88 degree_days <- 0
83 } 89 }
84 else { 90 else {
85 # Initialize hourly temperature. 91 # Initialize hourly temperature.
86 T <- NULL 92 T <- NULL
96 b <- 2.20 102 b <- 2.20
97 # Sunrise time. 103 # Sunrise time.
98 risetime <- 12 - y / 2 104 risetime <- 12 - y / 2
99 # Sunset time. 105 # Sunset time.
100 settime <- 12 + y / 2 106 settime <- 12 + y / 2
101 ts <- (dxp - dnp) * sin(pi * (settime - 5) / (y + 2 * a)) + dnp 107 ts <- (curr_max_temp - curr_min_temp) * sin(pi * (settime - 5) / (y + 2 * a)) + curr_min_temp
102 for (i in 1:24) { 108 for (i in 1:24) {
103 if (i > risetime && i < settime) { 109 if (i > risetime && i < settime) {
104 # Number of hours after Tmin until sunset. 110 # Number of hours after Tmin until sunset.
105 m <- i - 5 111 m <- i - 5
106 T[i] = (dxp - dnp) * sin(pi * m / (y + 2 * a)) + dnp 112 T[i] = (curr_max_temp - curr_min_temp) * sin(pi * m / (y + 2 * a)) + curr_min_temp
107 if (T[i] < 8.4) { 113 if (T[i] < 8.4) {
108 dh[i] <- 0 114 dh[i] <- 0
109 } 115 }
110 else { 116 else {
111 dh[i] <- T[i] - 8.4 117 dh[i] <- T[i] - 8.4
112 } 118 }
113 } 119 }
114 else if (i > settime) { 120 else if (i > settime) {
115 n <- i - settime 121 n <- i - settime
116 T[i] = dnp + (ts - dnp) * exp( - b * n / z) 122 T[i] = curr_min_temp + (ts - curr_min_temp) * exp( - b * n / z)
117 if (T[i] < 8.4) { 123 if (T[i] < 8.4) {
118 dh[i] <- 0 124 dh[i] <- 0
119 } 125 }
120 else { 126 else {
121 dh[i] <- T[i] - 8.4 127 dh[i] <- T[i] - 8.4
122 } 128 }
123 } 129 }
124 else { 130 else {
125 n <- i + 24 - settime 131 n <- i + 24 - settime
126 T[i]=dnp + (ts - dnp) * exp( - b * n / z) 132 T[i] = curr_min_temp + (ts - curr_min_temp) * exp( - b * n / z)
127 if (T[i] < 8.4) { 133 if (T[i] < 8.4) {
128 dh[i] <- 0 134 dh[i] <- 0
129 } 135 }
130 else { 136 else {
131 dh[i] <- T[i] - 8.4 137 dh[i] <- T[i] - 8.4
132 } 138 }
133 } 139 }
134 } 140 }
135 degree_days <- sum(dh) / 24 141 degree_days <- sum(dh) / 24
136 } 142 }
137 return(c(dmean, degree_days)) 143 return(c(curr_mean_temp, degree_days))
138 } 144 }
139 145
140 dev.egg = function(temperature) { 146 mortality.adult = function(temperature) {
141 dev.rate= -0.9843 * temperature + 33.438 147 if (temperature < 12.7) {
142 return(dev.rate) 148 mortality.probability = 0.002
143 } 149 }
144 150 else {
145 dev.young = function(temperature) { 151 mortality.probability = temperature * 0.0005 + 0.02
146 n12 <- -0.3728 * temperature + 14.68 152 }
147 n23 <- -0.6119 * temperature + 25.249 153 return(mortality.probability)
148 dev.rate = mean(n12 + n23)
149 return(dev.rate)
150 }
151
152 dev.old = function(temperature) {
153 n34 <- -0.6119 * temperature + 17.602
154 n45 <- -0.4408 * temperature + 19.036
155 dev.rate = mean(n34 + n45)
156 return(dev.rate)
157 }
158
159 dev.emerg = function(temperature) {
160 emerg.rate <- -0.5332 * temperature + 24.147
161 return(emerg.rate)
162 } 154 }
163 155
164 mortality.egg = function(temperature) { 156 mortality.egg = function(temperature) {
165 if (temperature < 12.7) { 157 if (temperature < 12.7) {
166 mort.prob = 0.8 158 mortality.probability = 0.8
167 } 159 }
168 else { 160 else {
169 mort.prob = 0.8 - temperature / 40.0 161 mortality.probability = 0.8 - temperature / 40.0
170 if (mort.prob < 0) { 162 if (mortality.probability < 0) {
171 mort.prob = 0.01 163 mortality.probability = 0.01
172 } 164 }
173 } 165 }
174 return(mort.prob) 166 return(mortality.probability)
175 } 167 }
176 168
177 mortality.nymph = function(temperature) { 169 mortality.nymph = function(temperature) {
178 if (temperature < 12.7) { 170 if (temperature < 12.7) {
179 mort.prob = 0.03 171 mortality.probability = 0.03
180 } 172 }
181 else { 173 else {
182 mort.prob = temperature * 0.0008 + 0.03 174 mortality.probability = temperature * 0.0008 + 0.03
183 } 175 }
184 return(mort.prob) 176 return(mortality.probability)
185 } 177 }
186 178
187 mortality.adult = function(temperature) { 179 parse_input_data = function(input_file, num_rows) {
188 if (temperature < 12.7) { 180 # Read in the input temperature datafile into a data frame.
189 mort.prob = 0.002 181 temperature_data_frame <- read.csv(file=input_file, header=T, strip.white=TRUE, sep=",")
190 } 182 num_columns <- dim(temperature_data_frame)[2]
191 else { 183 if (num_columns == 6) {
192 mort.prob = temperature * 0.0005 + 0.02 184 # The input data has the following 6 columns:
193 } 185 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
194 return(mort.prob) 186 # Set the column names for access when adding daylight length..
187 colnames(temperature_data_frame) <- c("LATITUDE","LONGITUDE", "DATE", "DOY", "TMIN", "TMAX")
188 # Add a column containing the daylight length for each day.
189 temperature_data_frame <- add_daylight_length(temperature_data_frame, num_columns, num_rows)
190 # Reset the column names with the additional column for later access.
191 colnames(temperature_data_frame) <- c("LATITUDE","LONGITUDE", "DATE", "DOY", "TMIN", "TMAX", "DAYLEN")
192 }
193 return(temperature_data_frame)
194 }
195
196 render_chart = function(chart_type, insect, location, latitude, start_date, end_date, days, maxval, plot_std_error,
197 group1, group2, group3, group1_std_error, group2_std_error, group3_std_error) {
198 if (chart_type == "pop_size_by_life_stage") {
199 title <- paste(insect, ": Total pop. by life stage :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ")
200 legend("topleft", c("Egg", "Nymph", "Adult"), lty=c(1, 1, 1), col=c(4, 2, 1), cex=3)
201 } else if (chart_type == "pop_size_by_generation") {
202 title <- paste(insect, ": Total pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ")
203 legend("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
204 } else if (chart_type == "adult_pop_size_by_generation") {
205 title <- paste(insect, ": Adult pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ")
206 legend("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
207 }
208 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)
209 lines(days, group2, lwd=2, lty=1, col=2)
210 lines(days, group3, lwd=2, lty=1, col=4)
211 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"))
212 axis(2, cex.axis=3)
213 if (plot_std_error==1) {
214 # Standard error for group1.
215 lines(days, group1+group1_std_error, lty=2)
216 lines (days, group1-group1_std_error, lty=2)
217 # Standard error for group2.
218 lines(days, group2+group2_std_error, col=2, lty=2)
219 lines(days, group2-group2_std_error, col=2, lty=2)
220 # Standard error for group3.
221 lines(days, group3+group3_std_error, col=4, lty=2)
222 lines(days, group3-group3_std_error, col=4, lty=2)
223 }
195 } 224 }
196 225
197 temperature_data_frame <- parse_input_data(opt$input, opt$num_days) 226 temperature_data_frame <- parse_input_data(opt$input, opt$num_days)
198 # All latitude values are the same, 227 # All latitude values are the same,
199 # so get the value from the first row. 228 # so get the value from the first row.
200 latitude <- temperature_data_frame$LATITUDE[1] 229 latitude <- temperature_data_frame$LATITUDE[1]
201 230
202 cat("Number of days: ", opt$num_days, "\n") 231 cat("Number of days: ", opt$num_days, "\n")
203 232
204 # Initialize matrices. 233 # Initialize matrices.
205 S0.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 234 Eggs.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
206 S1.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 235 YoungNymphs.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
207 S2.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 236 OldNymphs.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
208 S3.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 237 Previtellogenic.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
209 S4.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 238 Vitellogenic.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
210 S5.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 239 Diapausing.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
211 newborn.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 240 newborn.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
212 death.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 241 death.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
213 adult.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 242 adult.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
214 pop.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 243 population.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
215 P.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 244 P.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
216 F1.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 245 F1.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
217 F2.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 246 F2.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
218 P_adults.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 247 P_adults.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
219 F1_adults.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 248 F1_adults.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
220 F2_adults.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 249 F2_adults.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
221 250
222 # Loop through replications. 251 # Process replications.
223 for (N.replications in 1:opt$replications) { 252 for (N.replications in 1:opt$replications) {
224 # During each replication start with 1000 individuals. 253 # During each replication start with 1000 individuals.
225 # TODO: user definable as well? 254 # TODO: user definable as well?
226 num_insects <- 1000 255 num_insects <- 1000
227 # Generation, Stage, degree-days, T, Diapause. 256 # Generation, Stage, degree-days, T, Diapause.
233 # Time series of population size. 262 # Time series of population size.
234 total.population <- NULL 263 total.population <- NULL
235 overwintering_adult.population <- rep(0, opt$num_days) 264 overwintering_adult.population <- rep(0, opt$num_days)
236 first_generation.population <- rep(0, opt$num_days) 265 first_generation.population <- rep(0, opt$num_days)
237 second_generation.population <- rep(0, opt$num_days) 266 second_generation.population <- rep(0, opt$num_days)
238 S0 <- rep(0, opt$num_days) 267 Eggs <- rep(0, opt$num_days)
239 S1 <- rep(0, opt$num_days) 268 YoungNymphs <- rep(0, opt$num_days)
240 S2 <- rep(0, opt$num_days) 269 OldNymphs <- rep(0, opt$num_days)
241 S3 <- rep(0, opt$num_days) 270 Previtellogenic <- rep(0, opt$num_days)
242 S4 <- rep(0, opt$num_days) 271 Vitellogenic <- rep(0, opt$num_days)
243 S5 <- rep(0, opt$num_days) 272 Diapausing <- rep(0, opt$num_days)
244 P.adult <- rep(0, opt$num_days) 273 P.adult <- rep(0, opt$num_days)
245 F1.adult <- rep(0, opt$num_days) 274 F1.adult <- rep(0, opt$num_days)
246 F2.adult <- rep(0, opt$num_days) 275 F2.adult <- rep(0, opt$num_days)
247 N.newborn <- rep(0, opt$num_days) 276 N.newborn <- rep(0, opt$num_days)
248 N.death <- rep(0, opt$num_days) 277 N.death <- rep(0, opt$num_days)
263 # Newborn. 292 # Newborn.
264 birth.vector <- NULL 293 birth.vector <- NULL
265 # All individuals. 294 # All individuals.
266 for (i in 1:num_insects) { 295 for (i in 1:num_insects) {
267 # Find individual record. 296 # Find individual record.
268 vector.ind <- vector.matrix[i,] 297 vector.individual <- vector.matrix[i,]
269 # First of all, still alive? 298 # First of all, still alive?
270 # Adjustment for late season mortality rate. 299 # Adjustment for late season mortality rate.
271 if (latitude < 40.0) { 300 if (latitude < 40.0) {
272 post.mort <- 1 301 post.mortality <- 1
273 day.kill <- 300 302 day.kill <- 300
274 } 303 }
275 else { 304 else {
276 post.mort <- 2 305 post.mortality <- 2
277 day.kill <- 250 306 day.kill <- 250
278 } 307 }
279 if (vector.ind[2] == 0) { 308 if (vector.individual[2] == 0) {
280 # Egg. 309 # Egg.
281 death.prob = opt$egg_mort * mortality.egg(mean.temp) 310 death.probability = opt$egg_mortality * mortality.egg(mean.temp)
282 } 311 }
283 else if (vector.ind[2] == 1 | vector.ind[2] == 2) { 312 else if (vector.individual[2] == 1 | vector.individual[2] == 2) {
284 death.prob = opt$nymph_mort * mortality.nymph(mean.temp) 313 death.probability = opt$nymph_mortality * mortality.nymph(mean.temp)
285 } 314 }
286 else if (vector.ind[2] == 3 | vector.ind[2] == 4 | vector.ind[2] == 5) { 315 else if (vector.individual[2] == 3 | vector.individual[2] == 4 | vector.individual[2] == 5) {
287 # For adult. 316 # For adult.
288 if (doy < day.kill) { 317 if (doy < day.kill) {
289 death.prob = opt$adult_mort * mortality.adult(mean.temp) 318 death.probability = opt$adult_mortality * mortality.adult(mean.temp)
290 } 319 }
291 else { 320 else {
292 # Increase adult mortality after fall equinox. 321 # Increase adult mortality after fall equinox.
293 death.prob = opt$adult_mort * post.mort * mortality.adult(mean.temp) 322 death.probability = opt$adult_mortality * post.mortality * mortality.adult(mean.temp)
294 } 323 }
295 } 324 }
296 # (or dependent on temperature and life stage?) 325 # (or dependent on temperature and life stage?)
297 u.d <- runif(1) 326 u.d <- runif(1)
298 if (u.d < death.prob) { 327 if (u.d < death.probability) {
299 death.vector <- c(death.vector, i) 328 death.vector <- c(death.vector, i)
300 } 329 }
301 else { 330 else {
302 # Aggregrate index of dead bug. 331 # Aggregrate index of dead bug.
303 # Event 1 end of diapause. 332 # Event 1 end of diapause.
304 if (vector.ind[1] == 0 && vector.ind[2] == 3) { 333 if (vector.individual[1] == 0 && vector.individual[2] == 3) {
305 # Overwintering adult (previttelogenic). 334 # Overwintering adult (previttelogenic).
306 if (photoperiod > opt$photoperiod && vector.ind[3] > 68 && doy < 180) { 335 if (photoperiod > opt$photoperiod && vector.individual[3] > 68 && doy < 180) {
307 # Add 68C to become fully reproductively matured. 336 # Add 68C to become fully reproductively matured.
308 # Transfer to vittelogenic. 337 # Transfer to vittelogenic.
309 vector.ind <- c(0, 4, 0, 0, 0) 338 vector.individual <- c(0, 4, 0, 0, 0)
310 vector.matrix[i,] <- vector.ind 339 vector.matrix[i,] <- vector.individual
311 } 340 }
312 else { 341 else {
313 # Add to degree_days. 342 # Add to degree_days.
314 vector.ind[3] <- vector.ind[3] + degree_days.temp 343 vector.individual[3] <- vector.individual[3] + degree_days.temp
315 # Add 1 day in current stage. 344 # Add 1 day in current stage.
316 vector.ind[4] <- vector.ind[4] + 1 345 vector.individual[4] <- vector.individual[4] + 1
317 vector.matrix[i,] <- vector.ind 346 vector.matrix[i,] <- vector.individual
318 } 347 }
319 } 348 }
320 if (vector.ind[1] != 0 && vector.ind[2] == 3) { 349 if (vector.individual[1] != 0 && vector.individual[2] == 3) {
321 # Not overwintering adult (previttelogenic). 350 # Not overwintering adult (previttelogenic).
322 current.gen <- vector.ind[1] 351 current.gen <- vector.individual[1]
323 if (vector.ind[3] > 68) { 352 if (vector.individual[3] > 68) {
324 # Add 68C to become fully reproductively matured. 353 # Add 68C to become fully reproductively matured.
325 # Transfer to vittelogenic. 354 # Transfer to vittelogenic.
326 vector.ind <- c(current.gen, 4, 0, 0, 0) 355 vector.individual <- c(current.gen, 4, 0, 0, 0)
327 vector.matrix[i,] <- vector.ind 356 vector.matrix[i,] <- vector.individual
328 } 357 }
329 else { 358 else {
330 # Add to degree_days. 359 # Add to degree_days.
331 vector.ind[3] <- vector.ind[3] + degree_days.temp 360 vector.individual[3] <- vector.individual[3] + degree_days.temp
332 # Add 1 day in current stage. 361 # Add 1 day in current stage.
333 vector.ind[4] <- vector.ind[4] + 1 362 vector.individual[4] <- vector.individual[4] + 1
334 vector.matrix[i,] <- vector.ind 363 vector.matrix[i,] <- vector.individual
335 } 364 }
336 } 365 }
337 # Event 2 oviposition -- where population dynamics comes from. 366 # Event 2 oviposition -- where population dynamics comes from.
338 if (vector.ind[2] == 4 && vector.ind[1] == 0 && mean.temp > 10) { 367 if (vector.individual[2] == 4 && vector.individual[1] == 0 && mean.temp > 10) {
339 # Vittelogenic stage, overwintering generation. 368 # Vittelogenic stage, overwintering generation.
340 if (vector.ind[4] == 0) { 369 if (vector.individual[4] == 0) {
341 # Just turned in vittelogenic stage. 370 # Just turned in vittelogenic stage.
342 num_insects.birth = round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size)) 371 num_insects.birth = round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size))
343 } 372 }
344 else { 373 else {
345 # Daily probability of birth. 374 # Daily probability of birth.
348 if (u1 < p.birth) { 377 if (u1 < p.birth) {
349 num_insects.birth = round(runif(1, 2, 8)) 378 num_insects.birth = round(runif(1, 2, 8))
350 } 379 }
351 } 380 }
352 # Add to degree_days. 381 # Add to degree_days.
353 vector.ind[3] <- vector.ind[3] + degree_days.temp 382 vector.individual[3] <- vector.individual[3] + degree_days.temp
354 # Add 1 day in current stage. 383 # Add 1 day in current stage.
355 vector.ind[4] <- vector.ind[4] + 1 384 vector.individual[4] <- vector.individual[4] + 1
356 vector.matrix[i,] <- vector.ind 385 vector.matrix[i,] <- vector.individual
357 if (num_insects.birth > 0) { 386 if (num_insects.birth > 0) {
358 # Add new birth -- might be in different generations. 387 # Add new birth -- might be in different generations.
359 new.gen <- vector.ind[1] + 1 388 new.gen <- vector.individual[1] + 1
360 # Egg profile. 389 # Egg profile.
361 new.ind <- c(new.gen, 0, 0, 0, 0) 390 new.individual <- c(new.gen, 0, 0, 0, 0)
362 new.vector <- rep(new.ind, num_insects.birth) 391 new.vector <- rep(new.individual, num_insects.birth)
363 # Update batch of egg profile. 392 # Update batch of egg profile.
364 new.vector <- t(matrix(new.vector, nrow=5)) 393 new.vector <- t(matrix(new.vector, nrow=5))
365 # Group with total eggs laid in that day. 394 # Group with total eggs laid in that day.
366 birth.vector <- rbind(birth.vector, new.vector) 395 birth.vector <- rbind(birth.vector, new.vector)
367 } 396 }
368 } 397 }
369 # Event 2 oviposition -- for generation 1. 398 # Event 2 oviposition -- for generation 1.
370 if (vector.ind[2] == 4 && vector.ind[1] == 1 && mean.temp > 12.5 && doy < 222) { 399 if (vector.individual[2] == 4 && vector.individual[1] == 1 && mean.temp > 12.5 && doy < 222) {
371 # Vittelogenic stage, 1st generation 400 # Vittelogenic stage, 1st generation
372 if (vector.ind[4] == 0) { 401 if (vector.individual[4] == 0) {
373 # Just turned in vittelogenic stage. 402 # Just turned in vittelogenic stage.
374 num_insects.birth=round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size)) 403 num_insects.birth = round(runif(1, 2+opt$min_clutch_size, 8+opt$max_clutch_size))
375 } 404 }
376 else { 405 else {
377 # Daily probability of birth. 406 # Daily probability of birth.
378 p.birth = opt$oviposition * 0.01 407 p.birth = opt$oviposition * 0.01
379 u1 <- runif(1) 408 u1 <- runif(1)
380 if (u1 < p.birth) { 409 if (u1 < p.birth) {
381 num_insects.birth = round(runif(1, 2, 8)) 410 num_insects.birth = round(runif(1, 2, 8))
382 } 411 }
383 } 412 }
384 # Add to degree_days. 413 # Add to degree_days.
385 vector.ind[3] <- vector.ind[3] + degree_days.temp 414 vector.individual[3] <- vector.individual[3] + degree_days.temp
386 # Add 1 day in current stage. 415 # Add 1 day in current stage.
387 vector.ind[4] <- vector.ind[4] + 1 416 vector.individual[4] <- vector.individual[4] + 1
388 vector.matrix[i,] <- vector.ind 417 vector.matrix[i,] <- vector.individual
389 if (num_insects.birth > 0) { 418 if (num_insects.birth > 0) {
390 # Add new birth -- might be in different generations. 419 # Add new birth -- might be in different generations.
391 new.gen <- vector.ind[1] + 1 420 new.gen <- vector.individual[1] + 1
392 # Egg profile. 421 # Egg profile.
393 new.ind <- c(new.gen, 0, 0, 0, 0) 422 new.individual <- c(new.gen, 0, 0, 0, 0)
394 new.vector <- rep(new.ind, num_insects.birth) 423 new.vector <- rep(new.individual, num_insects.birth)
395 # Update batch of egg profile. 424 # Update batch of egg profile.
396 new.vector <- t(matrix(new.vector, nrow=5)) 425 new.vector <- t(matrix(new.vector, nrow=5))
397 # Group with total eggs laid in that day. 426 # Group with total eggs laid in that day.
398 birth.vector <- rbind(birth.vector, new.vector) 427 birth.vector <- rbind(birth.vector, new.vector)
399 } 428 }
400 } 429 }
401 # Event 3 development (with diapause determination). 430 # Event 3 development (with diapause determination).
402 # Event 3.1 egg development to young nymph (vector.ind[2]=0 -> egg). 431 # Event 3.1 egg development to young nymph (vector.individual[2]=0 -> egg).
403 if (vector.ind[2] == 0) { 432 if (vector.individual[2] == 0) {
404 # Egg stage. 433 # Egg stage.
405 # Add to degree_days. 434 # Add to degree_days.
406 vector.ind[3] <- vector.ind[3] + degree_days.temp 435 vector.individual[3] <- vector.individual[3] + degree_days.temp
407 if (vector.ind[3] >= (68 + opt$young_nymph_accum)) { 436 if (vector.individual[3] >= (68+opt$young_nymph_accumulation)) {
408 # From egg to young nymph, degree-days requirement met. 437 # From egg to young nymph, degree-days requirement met.
409 current.gen <- vector.ind[1] 438 current.gen <- vector.individual[1]
410 # Transfer to young nymph stage. 439 # Transfer to young nymph stage.
411 vector.ind <- c(current.gen, 1, 0, 0, 0) 440 vector.individual <- c(current.gen, 1, 0, 0, 0)
412 } 441 }
413 else { 442 else {
414 # Add 1 day in current stage. 443 # Add 1 day in current stage.
415 vector.ind[4] <- vector.ind[4] + 1 444 vector.individual[4] <- vector.individual[4] + 1
416 } 445 }
417 vector.matrix[i,] <- vector.ind 446 vector.matrix[i,] <- vector.individual
418 } 447 }
419 # Event 3.2 young nymph to old nymph (vector.ind[2]=1 -> young nymph: determines diapause). 448 # Event 3.2 young nymph to old nymph (vector.individual[2]=1 -> young nymph: determines diapause).
420 if (vector.ind[2] == 1) { 449 if (vector.individual[2] == 1) {
421 # Young nymph stage. 450 # Young nymph stage.
422 # Add to degree_days. 451 # Add to degree_days.
423 vector.ind[3] <- vector.ind[3] + degree_days.temp 452 vector.individual[3] <- vector.individual[3] + degree_days.temp
424 if (vector.ind[3] >= (250 + opt$old_nymph_accum)) { 453 if (vector.individual[3] >= (250+opt$old_nymph_accumulation)) {
425 # From young to old nymph, degree_days requirement met. 454 # From young to old nymph, degree_days requirement met.
426 current.gen <- vector.ind[1] 455 current.gen <- vector.individual[1]
427 # Transfer to old nym stage. 456 # Transfer to old nym stage.
428 vector.ind <- c(current.gen, 2, 0, 0, 0) 457 vector.individual <- c(current.gen, 2, 0, 0, 0)
429 if (photoperiod < opt$photoperiod && doy > 180) { 458 if (photoperiod < opt$photoperiod && doy > 180) {
430 vector.ind[5] <- 1 459 vector.individual[5] <- 1
431 } # Prepare for diapausing. 460 } # Prepare for diapausing.
432 } 461 }
433 else { 462 else {
434 # Add 1 day in current stage. 463 # Add 1 day in current stage.
435 vector.ind[4] <- vector.ind[4] + 1 464 vector.individual[4] <- vector.individual[4] + 1
436 } 465 }
437 vector.matrix[i,] <- vector.ind 466 vector.matrix[i,] <- vector.individual
438 } 467 }
439 # Event 3.3 old nymph to adult: previttelogenic or diapausing? 468 # Event 3.3 old nymph to adult: previttelogenic or diapausing?
440 if (vector.ind[2] == 2) { 469 if (vector.individual[2] == 2) {
441 # Old nymph stage. 470 # Old nymph stage.
442 # Add to degree_days. 471 # Add to degree_days.
443 vector.ind[3] <- vector.ind[3] + degree_days.temp 472 vector.individual[3] <- vector.individual[3] + degree_days.temp
444 if (vector.ind[3] >= (200 + opt$adult_accum)) { 473 if (vector.individual[3] >= (200+opt$adult_accumulation)) {
445 # From old to adult, degree_days requirement met. 474 # From old to adult, degree_days requirement met.
446 current.gen <- vector.ind[1] 475 current.gen <- vector.individual[1]
447 if (vector.ind[5] == 0) { 476 if (vector.individual[5] == 0) {
448 # Non-diapausing adult -- previttelogenic. 477 # Non-diapausing adult -- previttelogenic.
449 vector.ind <- c(current.gen, 3, 0, 0, 0) 478 vector.individual <- c(current.gen, 3, 0, 0, 0)
450 } 479 }
451 else { 480 else {
452 # Diapausing. 481 # Diapausing.
453 vector.ind <- c(current.gen, 5, 0, 0, 1) 482 vector.individual <- c(current.gen, 5, 0, 0, 1)
454 } 483 }
455 } 484 }
456 else { 485 else {
457 # Add 1 day in current stage. 486 # Add 1 day in current stage.
458 vector.ind[4] <- vector.ind[4] + 1 487 vector.individual[4] <- vector.individual[4] + 1
459 } 488 }
460 vector.matrix[i,] <- vector.ind 489 vector.matrix[i,] <- vector.individual
461 } 490 }
462 # Event 4 growing of diapausing adult (unimportant, but still necessary). 491 # Event 4 growing of diapausing adult (unimportant, but still necessary).
463 if (vector.ind[2] == 5) { 492 if (vector.individual[2] == 5) {
464 vector.ind[3] <- vector.ind[3] + degree_days.temp 493 vector.individual[3] <- vector.individual[3] + degree_days.temp
465 vector.ind[4] <- vector.ind[4] + 1 494 vector.individual[4] <- vector.individual[4] + 1
466 vector.matrix[i,] <- vector.ind 495 vector.matrix[i,] <- vector.individual
467 } 496 }
468 } # Else if it is still alive. 497 } # Else if it is still alive.
469 } # End of the individual bug loop. 498 } # End of the individual bug loop.
470 # Find the number of deaths. 499 # Find the number of deaths.
471 num_insects.death <- length(death.vector) 500 num_insects.death <- length(death.vector)
472 if (num_insects.death > 0) { 501 if (num_insects.death > 0) {
473 # Remove record of dead. 502 # Remove record of dead.
474 vector.matrix <- vector.matrix[-death.vector, ] 503 vector.matrix <- vector.matrix[-death.vector, ]
475 } 504 }
476 # Find the number of births. 505 # Find the number of births.
477 num_insects.newborn <- length(birth.vector[,1]) 506 num_insects.newborn <- length(birth.vector[,1])
478 vector.matrix <- rbind(vector.matrix, birth.vector) 507 vector.matrix <- rbind(vector.matrix, birth.vector)
479 # Update population size for the next day. 508 # Update population size for the next day.
480 num_insects <- num_insects - num_insects.death + num_insects.newborn 509 num_insects <- num_insects - num_insects.death + num_insects.newborn
481 510
482 # Aggregate results by day. 511 # Aggregate results by day.
483 total.population <- c(total.population, num_insects) 512 total.population <- c(total.population, num_insects)
484
485 # All adults population size. 513 # All adults population size.
486 num_insects.adult <- sum(vector.matrix[,2] == 3) + sum(vector.matrix[,2] == 4) + sum(vector.matrix[,2] == 5) 514 num_insects.adult <- sum(vector.matrix[,2]==3)+sum(vector.matrix[,2]==4)+sum(vector.matrix[,2]==5)
487
488 # Overwintering adult population size. 515 # Overwintering adult population size.
489 overwintering_adult.population[row] <- sum(vector.matrix[,1] == 0) 516 overwintering_adult.population[row] <- sum(vector.matrix[,1]==0)
490 # First generation population size. 517 # First generation population size.
491 first_generation.population[row] <- sum(vector.matrix[,1] == 1) 518 first_generation.population[row] <- sum(vector.matrix[,1]==1)
492 # Second generation population size. 519 # Second generation population size.
493 second_generation.population[row] <- sum(vector.matrix[,1] == 2) 520 second_generation.population[row] <- sum(vector.matrix[,1]==2)
494
495 # Egg population size. 521 # Egg population size.
496 S0[row] <- sum(vector.matrix[,2]==0) 522 Eggs[row] <- sum(vector.matrix[,2]==0)
497 # Young nymph population size. 523 # Young nymph population size.
498 S1[row] <- sum(vector.matrix[,2]==1) 524 YoungNymphs[row] <- sum(vector.matrix[,2]==1)
499 # Old nymph population size. 525 # Old nymph population size.
500 S2[row] <- sum(vector.matrix[,2]==2) 526 OldNymphs[row] <- sum(vector.matrix[,2]==2)
501 # Previtellogenic population size. 527 # Previtellogenic population size.
502 S3[row] <- sum(vector.matrix[,2]==3) 528 Previtellogenic[row] <- sum(vector.matrix[,2]==3)
503 # Vitellogenic population size. 529 # Vitellogenic population size.
504 S4[row] <- sum(vector.matrix[,2]==4) 530 Vitellogenic[row] <- sum(vector.matrix[,2]==4)
505 # Diapausing population size. 531 # Diapausing population size.
506 S5[row] <- sum(vector.matrix[,2]==5) 532 Diapausing[row] <- sum(vector.matrix[,2]==5)
507 533 # P adult population size.
508 P.adult[row] <- sum(vector.matrix[,1]==0) 534 P.adult[row] <- sum(vector.matrix[,1]==0)
535 # F1 adult population size.
509 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)) 536 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))
537 # F2 adult population size
510 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)) 538 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))
511 539 # Newborn population size.
512 N.newborn[row] <- num_insects.newborn 540 N.newborn[row] <- num_insects.newborn
541 # Dead population size.
513 N.death[row] <- num_insects.death 542 N.death[row] <- num_insects.death
543 # Adult population size.
514 N.adult[row] <- num_insects.adult 544 N.adult[row] <- num_insects.adult
515 } # end of days specified in the input temperature data 545 } # End of days specified in the input temperature data.
516 546
517 degree_days.cum <- cumsum(degree_days.day) 547 degree_days.cum <- cumsum(degree_days.day)
518 548
519 # Collect all the outputs. 549 # Define the output values.
520 S0.replications[,N.replications] <- S0 550 Eggs.replications[,N.replications] <- Eggs
521 S1.replications[,N.replications] <- S1 551 YoungNymphs.replications[,N.replications] <- YoungNymphs
522 S2.replications[,N.replications] <- S2 552 OldNymphs.replications[,N.replications] <- OldNymphs
523 S3.replications[,N.replications] <- S3 553 Previtellogenic.replications[,N.replications] <- Previtellogenic
524 S4.replications[,N.replications] <- S4 554 Vitellogenic.replications[,N.replications] <- Vitellogenic
525 S5.replications[,N.replications] <- S5 555 Diapausing.replications[,N.replications] <- Diapausing
526 newborn.replications[,N.replications] <- N.newborn 556 newborn.replications[,N.replications] <- N.newborn
527 death.replications[,N.replications] <- N.death 557 death.replications[,N.replications] <- N.death
528 adult.replications[,N.replications] <- N.adult 558 adult.replications[,N.replications] <- N.adult
529 pop.replications[,N.replications] <- total.population 559 population.replications[,N.replications] <- total.population
530 P.replications[,N.replications] <- overwintering_adult.population 560 P.replications[,N.replications] <- overwintering_adult.population
531 F1.replications[,N.replications] <- first_generation.population 561 F1.replications[,N.replications] <- first_generation.population
532 F2.replications[,N.replications] <- second_generation.population 562 F2.replications[,N.replications] <- second_generation.population
533 P_adults.replications[,N.replications] <- P.adult 563 P_adults.replications[,N.replications] <- P.adult
534 F1_adults.replications[,N.replications] <- F1.adult 564 F1_adults.replications[,N.replications] <- F1.adult
535 F2_adults.replications[,N.replications] <- F2.adult 565 F2_adults.replications[,N.replications] <- F2.adult
536 } 566 }
537 567
538 # Mean value for adults. 568 # Mean value for adults.
539 adults <- apply((S3.replications + S4.replications + S5.replications), 1, mean) 569 adults <- apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, mean)
540 # Mean value for nymphs. 570 # Mean value for nymphs.
541 nymphs <- apply((S1.replications + S2.replications), 1, mean) 571 nymphs <- apply((YoungNymphs.replications+OldNymphs.replications), 1, mean)
542 # Mean value for eggs. 572 # Mean value for eggs.
543 eggs <- apply(S0.replications, 1, mean) 573 eggs <- apply(Eggs.replications, 1, mean)
544 # Mean value for P. 574 # Mean value for P.
545 P <- apply(P.replications, 1, mean) 575 P <- apply(P.replications, 1, mean)
546 # Mean value for F1. 576 # Mean value for F1.
547 F1 <- apply(F1.replications, 1, mean) 577 F1 <- apply(F1.replications, 1, mean)
548 # Mean value for F2. 578 # Mean value for F2.
551 P_adults <- apply(P_adults.replications, 1, mean) 581 P_adults <- apply(P_adults.replications, 1, mean)
552 # Mean value for F1 adults. 582 # Mean value for F1 adults.
553 F1_adults <- apply(F1_adults.replications, 1, mean) 583 F1_adults <- apply(F1_adults.replications, 1, mean)
554 # Mean value for F2 adults. 584 # Mean value for F2 adults.
555 F2_adults <- apply(F2_adults.replications, 1, mean) 585 F2_adults <- apply(F2_adults.replications, 1, mean)
556
557 # Standard error for adults. 586 # Standard error for adults.
558 adults.std_error <- apply((S3.replications + S4.replications + S5.replications), 1, sd) / sqrt(opt$replications) 587 adults.std_error <- apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications)
559 # Standard error for nymphs. 588 # Standard error for nymphs.
560 nymphs.std_error <- apply((S1.replications + S2.replications) / sqrt(opt$replications), 1, sd) 589 nymphs.std_error <- apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd)
561 # Standard error for eggs. 590 # Standard error for eggs.
562 eggs.std_error <- apply(S0.replications, 1, sd) / sqrt(opt$replications) 591 eggs.std_error <- apply(Eggs.replications, 1, sd) / sqrt(opt$replications)
563 # Standard error value for P. 592 # Standard error value for P.
564 P.std_error <- apply(P.replications, 1, sd) / sqrt(opt$replications) 593 P.std_error <- apply(P.replications, 1, sd) / sqrt(opt$replications)
565 # Standard error for F1. 594 # Standard error for F1.
566 F1.std_error <- apply(F1.replications, 1, sd) / sqrt(opt$replications) 595 F1.std_error <- apply(F1.replications, 1, sd) / sqrt(opt$replications)
567 # Standard error for F2. 596 # Standard error for F2.
571 # Standard error for F1 adult. 600 # Standard error for F1 adult.
572 F1_adults.std_error <- apply(F1_adults.replications, 1, sd) / sqrt(opt$replications) 601 F1_adults.std_error <- apply(F1_adults.replications, 1, sd) / sqrt(opt$replications)
573 # Standard error for F2 adult. 602 # Standard error for F2 adult.
574 F2_adults.std_error <- apply(F2_adults.replications, 1, sd) / sqrt(opt$replications) 603 F2_adults.std_error <- apply(F2_adults.replications, 1, sd) / sqrt(opt$replications)
575 604
576 dev.new(width=20, height=30) 605 dev.new(width=20, height=40)
577 606
578 # Start PDF device driver to save charts to output. 607 # Start PDF device driver to save charts to output.
579 pdf(file=opt$output, width=20, height=30, bg="white") 608 pdf(file=opt$output, width=20, height=40, bg="white")
580
581 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)) 609 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1))
582 610
583 # Data analysis and visualization plots 611 # Data analysis and visualization plots
584 # only within a single calendar year. 612 # only within a single calendar year.
585 day.all <- c(1:opt$num_days) 613 day.all <- c(1:opt$num_days)
586 start_date <- temperature_data_frame$DATE[1] 614 start_date <- temperature_data_frame$DATE[1]
587 end_date <- temperature_data_frame$DATE[opt$num_days] 615 end_date <- temperature_data_frame$DATE[opt$num_days]
588 616
589 # Subfigure 1: population size by life stage 617 # Subfigure 1: population size by life stage
590 title <- paste(opt$insect, ": Total pop. by life stage :", opt$location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ") 618 maxval <- max(eggs+eggs.std_error, nymphs+nymphs.std_error, adults+adults.std_error)
591 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) 619 render_chart("pop_size_by_life_stage", opt$insect, opt$location, latitude, start_date, end_date, days.all, maxval,
592 # Young and old nymphs. 620 opt$std_error_plot, adults, nymphs, eggs, adults.std_error, nymphs.std_error, eggs.std_error)
593 lines(day.all, nymphs, lwd=2, lty=1, col=2)
594 # Eggs
595 lines(day.all, eggs, lwd=2, lty=1, col=4)
596 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"))
597 axis(2, cex.axis=3)
598 legend("topleft", c("Egg", "Nymph", "Adult"), lty=c(1, 1, 1), col=c(4, 2, 1), cex=3)
599 if (opt$std_error_plot == 1) {
600 # Add Standard error lines to plot
601 # Standard error for adults
602 lines (day.all, adults+adults.std_error, lty=2)
603 lines (day.all, adults-adults.std_error, lty=2)
604 # Standard error for nymphs
605 lines (day.all, nymphs+nymphs.std_error, col=2, lty=2)
606 lines (day.all, nymphs-nymphs.std_error, col=2, lty=2)
607 # Standard error for eggs
608 lines (day.all, eggs+eggs.std_error, col=4, lty=2)
609 lines (day.all, eggs-eggs.std_error, col=4, lty=2)
610 }
611 621
612 # Subfigure 2: population size by generation 622 # Subfigure 2: population size by generation
613 title <- paste(opt$insect, ": Total pop. by generation :", opt$location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ") 623 maxval <- max(F2)
614 plot(day.all, P, main=title, type="l", ylim=c(0, max(F2)), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3) 624 render_chart("pop_size_by_generation", opt$insect, opt$location, latitude, start_date, end_date, days.all, maxval,
615 lines(day.all, F1, lwd = 2, lty = 1, col=2) 625 opt$std_error_plot, P, F1, F2, P.std_error, F1.std_error, F2.std_error)
616 lines(day.all, F2, lwd = 2, lty = 1, col=4) 626
617 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"))
618 axis(2, cex.axis=3)
619 legend("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
620 if (opt$std_error_plot == 1) {
621 # Add Standard error lines to plot
622 # Standard error for adults
623 lines (day.all, P+P.std_error, lty=2)
624 lines (day.all, P-P.std_error, lty=2)
625 # Standard error for nymphs
626 lines (day.all, F1+F1.std_error, col=2, lty=2)
627 lines (day.all, F1-F1.std_error, col=2, lty=2)
628 # Standard error for eggs
629 lines (day.all, F2+F2.std_error, col=4, lty=2)
630 lines (day.all, F2-F2.std_error, col=4, lty=2)
631 }
632 627
633 # Subfigure 3: adult population size by generation 628 # Subfigure 3: adult population size by generation
634 title <- paste(opt$insect, ": Adult pop. by generation :", opt$location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ") 629 maxval <- max(F2_adults) + 100
635 plot(day.all, P_adults, 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) 630 render_chart("adult_pop_size_by_generation", opt$insect, opt$location, latitude, start_date, end_date, days.all, maxval,
636 lines(day.all, F1_adults, lwd = 2, lty = 1, col=2) 631 opt$std_error_plot, P_adults, F1_adults, F2_adults, P_adults.std_error, F1_adults.std_error, F2_adults2.std_error)
637 lines(day.all, F2_adults, lwd = 2, lty = 1, col=4) 632
638 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"))
639 axis(2, cex.axis=3)
640 legend("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
641 if (opt$std_error_plot == 1) {
642 # Add Standard error lines to plot
643 # Standard error for adults
644 lines (day.all, P_adults+P_adults.std_error, lty=2)
645 lines (day.all, P_adults-P_adults.std_error, lty=2)
646 # Standard error for nymphs
647 lines (day.all, F1_adults+F1_adults.std_error, col=2, lty=2)
648 lines (day.all, F1_adults-F1_adults.std_error, col=2, lty=2)
649 # Standard error for eggs
650 lines (day.all, F2_adults+F2_adults.std_error, col=4, lty=2)
651 lines (day.all, F2_adults-F2_adults.std_error, col=4, lty=2)
652 }
653 633
654 # Turn off device driver to flush output. 634 # Turn off device driver to flush output.
655 dev.off() 635 dev.off()