Mercurial > repos > greg > insect_phenology_model
changeset 8:37f1ad91a949 draft
Uploaded
| author | greg |
|---|---|
| date | Tue, 13 Feb 2018 13:53:37 -0500 |
| parents | ad26f07a7dd8 |
| children | d9371485aaf5 |
| files | insect_phenology_model.R |
| diffstat | 1 files changed, 297 insertions(+), 273 deletions(-) [+] |
line wrap: on
line diff
--- a/insect_phenology_model.R Tue Feb 13 13:53:29 2018 -0500 +++ b/insect_phenology_model.R Tue Feb 13 13:53:37 2018 -0500 @@ -23,145 +23,165 @@ 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) -args <- parse_args(parser, positional_arguments=TRUE) -opt <- args$options +parser <- OptionParser(usage="%prog [options] file", option_list=option_list); +args <- parse_args(parser, positional_arguments=TRUE); +opt <- args$options; add_daylight_length = function(temperature_data_frame, num_columns, num_rows) { # Return a vector of daylight length (photoperido profile) for # the number of days specified in the input temperature data # (from Forsythe 1995). - p = 0.8333 - latitude <- temperature_data_frame$LATITUDE[1] - daylight_length_vector <- NULL + p = 0.8333; + latitude = temperature_data_frame$LATITUDE[1]; + daylight_length_vector = NULL; for (i in 1:num_rows) { # Get the day of the year from the current row # of the temperature data for computation. - doy <- temperature_data_frame$DOY[i] - theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (doy - 186))) - phi <- asin(0.39795 * cos(theta)) + doy = temperature_data_frame$DOY[i]; + theta = 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (doy - 186))); + phi = asin(0.39795 * cos(theta)); # Compute the length of daylight for the day of the year. - darkness_length <- 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi))) - daylight_length_vector[i] <- 24 - darkness_length + darkness_length = 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi))); + daylight_length_vector[i] = 24 - darkness_length; } # Append daylight_length_vector as a new column to temperature_data_frame. - temperature_data_frame[, num_columns+1] <- daylight_length_vector - return(temperature_data_frame) + temperature_data_frame[, num_columns+1] = daylight_length_vector; + return(temperature_data_frame); } dev.egg = function(temperature) { - dev.rate = -0.9843 * temperature + 33.438 - return(dev.rate) + dev.rate = -0.9843 * temperature + 33.438; + return(dev.rate); } dev.emerg = function(temperature) { - emerg.rate <- -0.5332 * temperature + 24.147 - return(emerg.rate) + emerg.rate = -0.5332 * temperature + 24.147; + return(emerg.rate); } dev.old = function(temperature) { - n34 <- -0.6119 * temperature + 17.602 - n45 <- -0.4408 * temperature + 19.036 - dev.rate = mean(n34 + n45) - return(dev.rate) + n34 = -0.6119 * temperature + 17.602; + n45 = -0.4408 * temperature + 19.036; + dev.rate = mean(n34 + n45); + return(dev.rate); } dev.young = function(temperature) { - n12 <- -0.3728 * temperature + 14.68 - n23 <- -0.6119 * temperature + 25.249 - dev.rate = mean(n12 + n23) - return(dev.rate) + n12 = -0.3728 * temperature + 14.68; + n23 = -0.6119 * temperature + 25.249; + dev.rate = mean(n12 + n23); + return(dev.rate); +} + + +get_date_labels = function(temperature_data_frame, num_rows) { + # Keep track of the years to see if spanning years. + month_labels = list(); + current_month_label = NULL; + for (i in 1:num_rows) { + # Get the year and month from the date which + # has the format YYYY-MM-DD. + date = format(temperature_data_frame$DATE[i]); + items = strsplit(date, "-")[[1]]; + month = items[2]; + month_label = month.abb[as.integer(month)]; + if (!identical(current_month_label, month_label)) { + month_labels[length(month_labels)+1] = month_label; + current_month_label = month_label; + } + } + return(c(unlist(month_labels))); } get_temperature_at_hour = function(latitude, temperature_data_frame, row, num_days) { - # Base development threshold for Brown Marmolated Stink Bug + # Base development threshold for Brown Marmorated Stink Bug # insect phenology model. - threshold <- 14.17 + threshold = 14.17; # Minimum temperature for current row. - curr_min_temp <- temperature_data_frame$TMIN[row] + curr_min_temp = temperature_data_frame$TMIN[row]; # Maximum temperature for current row. - curr_max_temp <- temperature_data_frame$TMAX[row] + curr_max_temp = temperature_data_frame$TMAX[row]; # Mean temperature for current row. - curr_mean_temp <- 0.5 * (curr_min_temp + curr_max_temp) + curr_mean_temp = 0.5 * (curr_min_temp + curr_max_temp); # Initialize degree day accumulation - averages <- 0 + averages = 0; if (curr_max_temp < threshold) { - averages <- 0 + averages = 0; } else { # Initialize hourly temperature. - T <- NULL + T = NULL; # Initialize degree hour vector. - dh <- NULL + dh = NULL; # Daylight length for current row. - y <- temperature_data_frame$DAYLEN[row] + y = temperature_data_frame$DAYLEN[row]; # Darkness length. - z <- 24 - y + z = 24 - y; # Lag coefficient. - a <- 1.86 + a = 1.86; # Darkness coefficient. - b <- 2.20 + b = 2.20; # Sunrise time. - risetime <- 12 - y / 2 + risetime = 12 - y / 2; # Sunset time. - settime <- 12 + y / 2 - ts <- (curr_max_temp - curr_min_temp) * sin(pi * (settime - 5) / (y + 2 * a)) + curr_min_temp + settime = 12 + y / 2; + ts = (curr_max_temp - curr_min_temp) * sin(pi * (settime - 5) / (y + 2 * a)) + curr_min_temp; for (i in 1:24) { if (i > risetime && i < settime) { # Number of hours after Tmin until sunset. - m <- i - 5 - T[i] = (curr_max_temp - curr_min_temp) * sin(pi * m / (y + 2 * a)) + curr_min_temp + m = i - 5; + T[i] = (curr_max_temp - curr_min_temp) * sin(pi * m / (y + 2 * a)) + curr_min_temp; if (T[i] < 8.4) { - dh[i] <- 0 + dh[i] = 0; } else { - dh[i] <- T[i] - 8.4 + dh[i] = T[i] - 8.4; } } else if (i > settime) { - n <- i - settime - T[i] = curr_min_temp + (ts - curr_min_temp) * exp( - b * n / z) + n = i - settime; + T[i] = curr_min_temp + (ts - curr_min_temp) * exp( - b * n / z); if (T[i] < 8.4) { - dh[i] <- 0 + dh[i] = 0; } else { - dh[i] <- T[i] - 8.4 + dh[i] = T[i] - 8.4; } } else { - n <- i + 24 - settime - T[i] = curr_min_temp + (ts - curr_min_temp) * exp( - b * n / z) + n = i + 24 - settime; + T[i] = curr_min_temp + (ts - curr_min_temp) * exp( - b * n / z); if (T[i] < 8.4) { - dh[i] <- 0 + dh[i] = 0; } else { - dh[i] <- T[i] - 8.4 + dh[i] = T[i] - 8.4; } } } - averages <- sum(dh) / 24 + averages = sum(dh) / 24; } return(c(curr_mean_temp, averages)) } mortality.adult = function(temperature) { if (temperature < 12.7) { - mortality.probability = 0.002 + mortality.probability = 0.002; } else { - mortality.probability = temperature * 0.0005 + 0.02 + mortality.probability = temperature * 0.0005 + 0.02; } return(mortality.probability) } mortality.egg = function(temperature) { if (temperature < 12.7) { - mortality.probability = 0.8 + mortality.probability = 0.8; } else { - mortality.probability = 0.8 - temperature / 40.0 + mortality.probability = 0.8 - temperature / 40.0; if (mortality.probability < 0) { - mortality.probability = 0.01 + mortality.probability = 0.01; } } return(mortality.probability) @@ -169,171 +189,174 @@ mortality.nymph = function(temperature) { if (temperature < 12.7) { - mortality.probability = 0.03 + mortality.probability = 0.03; } else { - mortality.probability = temperature * 0.0008 + 0.03 + mortality.probability = temperature * 0.0008 + 0.03; } - return(mortality.probability) + return(mortality.probability); } parse_input_data = function(input_file, num_rows) { # Read in the input temperature datafile into a data frame. - temperature_data_frame <- read.csv(file=input_file, header=T, strip.white=TRUE, sep=",") - num_columns <- dim(temperature_data_frame)[2] + temperature_data_frame = read.csv(file=input_file, header=T, strip.white=TRUE, sep=","); + num_columns = dim(temperature_data_frame)[2]; if (num_columns == 6) { # The input data has the following 6 columns: # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX # Set the column names for access when adding daylight length.. - colnames(temperature_data_frame) <- c("LATITUDE","LONGITUDE", "DATE", "DOY", "TMIN", "TMAX") + colnames(temperature_data_frame) = c("LATITUDE","LONGITUDE", "DATE", "DOY", "TMIN", "TMAX"); # Add a column containing the daylight length for each day. - temperature_data_frame <- add_daylight_length(temperature_data_frame, num_columns, num_rows) + temperature_data_frame = add_daylight_length(temperature_data_frame, num_columns, num_rows); # Reset the column names with the additional column for later access. - colnames(temperature_data_frame) <- c("LATITUDE","LONGITUDE", "DATE", "DOY", "TMIN", "TMAX", "DAYLEN") + colnames(temperature_data_frame) = c("LATITUDE","LONGITUDE", "DATE", "DOY", "TMIN", "TMAX", "DAYLEN"); } - return(temperature_data_frame) + return(temperature_data_frame); } + render_chart = function(chart_type, insect, location, latitude, start_date, end_date, days, maxval, plot_std_error, - group1, group2, group3, group1_std_error, group2_std_error, group3_std_error) { + group1, group2, group3, group1_std_error, group2_std_error, group3_std_error, date_labels) { if (chart_type == "pop_size_by_life_stage") { - title <- paste(insect, ": Total pop. by life stage :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ") - legend_text <- c("Egg", "Nymph", "Adult") - columns <- c(4, 2, 1) + title = paste(insect, ": Total pop. by life stage :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" "); + legend_text = c("Egg", "Nymph", "Adult"); + columns = c(4, 2, 1); } else if (chart_type == "pop_size_by_generation") { - title <- paste(insect, ": Total pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ") - legend_text <- c("P", "F1", "F2") - columns <- c(1, 2, 4) + title = paste(insect, ": Total pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" "); + legend_text = c("P", "F1", "F2"); + columns = c(1, 2, 4); } else if (chart_type == "adult_pop_size_by_generation") { - title <- paste(insect, ": Adult pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ") - legend_text <- c("P", "F1", "F2") - columns <- c(1, 2, 4) + title = paste(insect, ": Adult pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" "); + legend_text = c("P", "F1", "F2"); + columns = c(1, 2, 4); } - 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) - legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3) - lines(days, group2, lwd=2, lty=1, col=2) - lines(days, group3, 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) + 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); + legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3); + lines(days, group2, lwd=2, lty=1, col=2); + lines(days, group3, lwd=2, lty=1, col=4); + axis(1, at=c(1:length(date_labels)) * 30 - 15, cex.axis=3, labels=date_labels); + axis(2, cex.axis=3); if (plot_std_error==1) { # Standard error for group1. - lines(days, group1+group1_std_error, lty=2) - lines (days, group1-group1_std_error, lty=2) + lines(days, group1+group1_std_error, lty=2); + lines(days, group1-group1_std_error, lty=2); # Standard error for group2. - lines(days, group2+group2_std_error, col=2, lty=2) - lines(days, group2-group2_std_error, col=2, lty=2) + lines(days, group2+group2_std_error, col=2, lty=2); + lines(days, group2-group2_std_error, col=2, lty=2); # Standard error for group3. - lines(days, group3+group3_std_error, col=4, lty=2) - lines(days, group3-group3_std_error, col=4, lty=2) + lines(days, group3+group3_std_error, col=4, lty=2); + lines(days, group3-group3_std_error, col=4, lty=2); } } -temperature_data_frame <- parse_input_data(opt$input, opt$num_days) +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. -latitude <- temperature_data_frame$LATITUDE[1] +latitude = temperature_data_frame$LATITUDE[1]; +num_columns = dim(temperature_data_frame)[2]; +date_labels = get_date_labels(temperature_data_frame, opt$num_days); # Initialize matrices. -Eggs.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) -YoungNymphs.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) -OldNymphs.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) -Previtellogenic.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) -Vitellogenic.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) -Diapausing.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +Eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); +YoungNymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); +OldNymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); +Previtellogenic.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); +Vitellogenic.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); +Diapausing.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); -newborn.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) -adult.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) -death.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +newborn.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); +adult.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); +death.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); -P.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) -P_adults.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) -F1.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) -F1_adults.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) -F2.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) +P.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); +P_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); +F1.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); +F1_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); +F2.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); -population.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) +population.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); # Process replications. for (N.replications in 1:opt$replications) { # Start with the user-defined number of insects per replication. - num_insects <- opt$insects_per_replication + num_insects = opt$insects_per_replication; # Generation, Stage, degree-days, T, Diapause. - vector.ini <- c(0, 3, 0, 0, 0) + vector.ini = c(0, 3, 0, 0, 0); # Overwintering, previttelogenic, degree-days=0, T=0, no-diapause. - vector.matrix <- rep(vector.ini, num_insects) + vector.matrix = rep(vector.ini, num_insects); # Complete matrix for the population. - vector.matrix <- base::t(matrix(vector.matrix, nrow=5)) + vector.matrix = base::t(matrix(vector.matrix, nrow=5)); # Time series of population size. - Eggs <- rep(0, opt$num_days) - YoungNymphs <- rep(0, opt$num_days) - OldNymphs <- rep(0, opt$num_days) - Previtellogenic <- rep(0, opt$num_days) - Vitellogenic <- rep(0, opt$num_days) - Diapausing <- rep(0, opt$num_days) + Eggs = rep(0, opt$num_days); + YoungNymphs = rep(0, opt$num_days); + OldNymphs = rep(0, opt$num_days); + Previtellogenic = rep(0, opt$num_days); + Vitellogenic = rep(0, opt$num_days); + Diapausing = rep(0, opt$num_days); - N.newborn <- rep(0, opt$num_days) - N.adult <- rep(0, opt$num_days) - N.death <- rep(0, opt$num_days) + N.newborn = rep(0, opt$num_days); + N.adult = rep(0, opt$num_days); + N.death = rep(0, opt$num_days); - overwintering_adult.population <- rep(0, opt$num_days) - first_generation.population <- rep(0, opt$num_days) - second_generation.population <- rep(0, opt$num_days) + overwintering_adult.population = rep(0, opt$num_days); + first_generation.population = rep(0, opt$num_days); + second_generation.population = rep(0, opt$num_days); - P.adult <- rep(0, opt$num_days) - F1.adult <- rep(0, opt$num_days) - F2.adult <- rep(0, opt$num_days) + P.adult = rep(0, opt$num_days); + F1.adult = rep(0, opt$num_days); + F2.adult = rep(0, opt$num_days); - total.population <- NULL + total.population = NULL; - averages.day <- rep(0, opt$num_days) + averages.day = rep(0, opt$num_days); # All the days included in the input temperature dataset. for (row in 1:opt$num_days) { # Get the integer day of the year for the current row. - doy <- temperature_data_frame$DOY[row] + doy = temperature_data_frame$DOY[row]; # Photoperiod in the day. - photoperiod <- temperature_data_frame$DAYLEN[row] - temp.profile <- get_temperature_at_hour(latitude, temperature_data_frame, row, opt$num_days) - mean.temp <- temp.profile[1] - averages.temp <- temp.profile[2] - averages.day[row] <- averages.temp + photoperiod = temperature_data_frame$DAYLEN[row]; + temp.profile = get_temperature_at_hour(latitude, temperature_data_frame, row, opt$num_days); + mean.temp = temp.profile[1]; + averages.temp = temp.profile[2]; + averages.day[row] = averages.temp; # Trash bin for death. - death.vector <- NULL + death.vector = NULL; # Newborn. - birth.vector <- NULL + birth.vector = NULL; # All individuals. for (i in 1:num_insects) { # Individual record. - vector.individual <- vector.matrix[i,] + vector.individual = vector.matrix[i,]; # Adjustment for late season mortality rate (still alive?). if (latitude < 40.0) { - post.mortality <- 1 - day.kill <- 300 + post.mortality = 1; + day.kill = 300; } else { - post.mortality <- 2 - day.kill <- 250 + post.mortality = 2; + day.kill = 250; } if (vector.individual[2] == 0) { # Egg. - death.probability = opt$egg_mortality * mortality.egg(mean.temp) + death.probability = opt$egg_mortality * mortality.egg(mean.temp); } else if (vector.individual[2] == 1 | vector.individual[2] == 2) { - death.probability = opt$nymph_mortality * mortality.nymph(mean.temp) + death.probability = opt$nymph_mortality * mortality.nymph(mean.temp); } else if (vector.individual[2] == 3 | vector.individual[2] == 4 | vector.individual[2] == 5) { # Adult. if (doy < day.kill) { - death.probability = opt$adult_mortality * mortality.adult(mean.temp) + death.probability = opt$adult_mortality * mortality.adult(mean.temp); } else { # Increase adult mortality after fall equinox. - death.probability = opt$adult_mortality * post.mortality * mortality.adult(mean.temp) + death.probability = opt$adult_mortality * post.mortality * mortality.adult(mean.temp); } } # Dependent on temperature and life stage? - u.d <- runif(1) + u.d = runif(1); if (u.d < death.probability) { - death.vector <- c(death.vector, i) + death.vector = c(death.vector, i); } else { # End of diapause. @@ -342,32 +365,32 @@ if (photoperiod > opt$photoperiod && vector.individual[3] > 68 && doy < 180) { # Add 68C to become fully reproductively matured. # Transfer to vittelogenic. - vector.individual <- c(0, 4, 0, 0, 0) - vector.matrix[i,] <- vector.individual + vector.individual = c(0, 4, 0, 0, 0); + vector.matrix[i,] = vector.individual; } else { # Add to # Add average temperature for current day. - vector.individual[3] <- vector.individual[3] + averages.temp + vector.individual[3] = vector.individual[3] + averages.temp; # Add 1 day in current stage. - vector.individual[4] <- vector.individual[4] + 1 - vector.matrix[i,] <- vector.individual + vector.individual[4] = vector.individual[4] + 1; + vector.matrix[i,] = vector.individual; } } if (vector.individual[1] != 0 && vector.individual[2] == 3) { # Not overwintering adult (previttelogenic). - current.gen <- vector.individual[1] + current.gen = vector.individual[1]; if (vector.individual[3] > 68) { # Add 68C to become fully reproductively matured. # Transfer to vittelogenic. - vector.individual <- c(current.gen, 4, 0, 0, 0) - vector.matrix[i,] <- vector.individual + vector.individual = c(current.gen, 4, 0, 0, 0); + vector.matrix[i,] = vector.individual; } else { # Add average temperature for current day. - vector.individual[3] <- vector.individual[3] + averages.temp + vector.individual[3] = vector.individual[3] + averages.temp; # Add 1 day in current stage. - vector.individual[4] <- vector.individual[4] + 1 - vector.matrix[i,] <- vector.individual + vector.individual[4] = vector.individual[4] + 1; + vector.matrix[i,] = vector.individual; } } # Oviposition -- where population dynamics comes from. @@ -375,31 +398,31 @@ # Vittelogenic stage, overwintering generation. if (vector.individual[4] == 0) { # Just turned in vittelogenic stage. - num_insects.birth = round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size)) + num_insects.birth = round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size)); } else { # Daily probability of birth. - p.birth = opt$oviposition * 0.01 - u1 <- runif(1) + p.birth = opt$oviposition * 0.01; + u1 = runif(1); if (u1 < p.birth) { - num_insects.birth = round(runif(1, 2, 8)) + num_insects.birth = round(runif(1, 2, 8)); } } # Add average temperature for current day. - vector.individual[3] <- vector.individual[3] + averages.temp + vector.individual[3] = vector.individual[3] + averages.temp; # Add 1 day in current stage. - vector.individual[4] <- vector.individual[4] + 1 - vector.matrix[i,] <- vector.individual + vector.individual[4] = vector.individual[4] + 1; + vector.matrix[i,] = vector.individual; if (num_insects.birth > 0) { # Add new birth -- might be in different generations. - new.gen <- vector.individual[1] + 1 + new.gen = vector.individual[1] + 1; # Egg profile. - new.individual <- c(new.gen, 0, 0, 0, 0) - new.vector <- rep(new.individual, num_insects.birth) + new.individual = c(new.gen, 0, 0, 0, 0); + new.vector = rep(new.individual, num_insects.birth); # Update batch of egg profile. - new.vector <- t(matrix(new.vector, nrow=5)) + new.vector = t(matrix(new.vector, nrow=5)); # Group with total eggs laid in that day. - birth.vector <- rbind(birth.vector, new.vector) + birth.vector = rbind(birth.vector, new.vector); } } # Oviposition -- for generation 1. @@ -407,244 +430,245 @@ # Vittelogenic stage, 1st generation if (vector.individual[4] == 0) { # Just turned in vittelogenic stage. - num_insects.birth = round(runif(1, 2+opt$min_clutch_size, 8+opt$max_clutch_size)) + num_insects.birth = round(runif(1, 2+opt$min_clutch_size, 8+opt$max_clutch_size)); } else { # Daily probability of birth. - p.birth = opt$oviposition * 0.01 - u1 <- runif(1) + p.birth = opt$oviposition * 0.01; + u1 = runif(1); if (u1 < p.birth) { - num_insects.birth = round(runif(1, 2, 8)) + num_insects.birth = round(runif(1, 2, 8)); } } # Add average temperature for current day. - vector.individual[3] <- vector.individual[3] + averages.temp + vector.individual[3] = vector.individual[3] + averages.temp; # Add 1 day in current stage. - vector.individual[4] <- vector.individual[4] + 1 - vector.matrix[i,] <- vector.individual + vector.individual[4] = vector.individual[4] + 1; + vector.matrix[i,] = vector.individual; if (num_insects.birth > 0) { # Add new birth -- might be in different generations. - new.gen <- vector.individual[1] + 1 + new.gen = vector.individual[1] + 1; # Egg profile. - new.individual <- c(new.gen, 0, 0, 0, 0) - new.vector <- rep(new.individual, num_insects.birth) + new.individual = c(new.gen, 0, 0, 0, 0); + new.vector = rep(new.individual, num_insects.birth); # Update batch of egg profile. - new.vector <- t(matrix(new.vector, nrow=5)) + new.vector = t(matrix(new.vector, nrow=5)); # Group with total eggs laid in that day. - birth.vector <- rbind(birth.vector, new.vector) + birth.vector = rbind(birth.vector, new.vector); } } # Egg to young nymph. if (vector.individual[2] == 0) { # Add average temperature for current day. - vector.individual[3] <- vector.individual[3] + averages.temp + vector.individual[3] = vector.individual[3] + averages.temp; if (vector.individual[3] >= (68+opt$young_nymph_accumulation)) { # From egg to young nymph, degree-days requirement met. - current.gen <- vector.individual[1] + current.gen = vector.individual[1]; # Transfer to young nymph stage. - vector.individual <- c(current.gen, 1, 0, 0, 0) + vector.individual = c(current.gen, 1, 0, 0, 0); } else { # Add 1 day in current stage. - vector.individual[4] <- vector.individual[4] + 1 + vector.individual[4] = vector.individual[4] + 1; } - vector.matrix[i,] <- vector.individual + vector.matrix[i,] = vector.individual; } # Young nymph to old nymph. if (vector.individual[2] == 1) { # Add average temperature for current day. - vector.individual[3] <- vector.individual[3] + averages.temp + vector.individual[3] = vector.individual[3] + averages.temp; if (vector.individual[3] >= (250+opt$old_nymph_accumulation)) { # From young to old nymph, degree_days requirement met. - current.gen <- vector.individual[1] + current.gen = vector.individual[1]; # Transfer to old nym stage. - vector.individual <- c(current.gen, 2, 0, 0, 0) + vector.individual = c(current.gen, 2, 0, 0, 0); if (photoperiod < opt$photoperiod && doy > 180) { - vector.individual[5] <- 1 + vector.individual[5] = 1; } # Prepare for diapausing. } else { # Add 1 day in current stage. - vector.individual[4] <- vector.individual[4] + 1 + vector.individual[4] = vector.individual[4] + 1; } - vector.matrix[i,] <- vector.individual + vector.matrix[i,] = vector.individual; } # Old nymph to adult: previttelogenic or diapausing? if (vector.individual[2] == 2) { # Add average temperature for current day. - vector.individual[3] <- vector.individual[3] + averages.temp + 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] + current.gen = vector.individual[1]; if (vector.individual[5] == 0) { # Previttelogenic. - vector.individual <- c(current.gen, 3, 0, 0, 0) + vector.individual = c(current.gen, 3, 0, 0, 0); } else { # Diapausing. - vector.individual <- c(current.gen, 5, 0, 0, 1) + vector.individual = c(current.gen, 5, 0, 0, 1); } } else { # Add 1 day in current stage. - vector.individual[4] <- vector.individual[4] + 1 + vector.individual[4] = vector.individual[4] + 1; } - vector.matrix[i,] <- vector.individual + vector.matrix[i,] = vector.individual; } # 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 - vector.matrix[i,] <- vector.individual + vector.individual[3] = vector.individual[3] + averages.temp; + vector.individual[4] = vector.individual[4] + 1; + vector.matrix[i,] = vector.individual; } } # Else if it is still alive. } # End of the individual bug loop. # Number of deaths. - num_insects.death <- length(death.vector) + num_insects.death = length(death.vector); if (num_insects.death > 0) { # Remove record of dead. - vector.matrix <- vector.matrix[-death.vector, ] + vector.matrix = vector.matrix[-death.vector,]; } # Number of births. - num_insects.newborn <- length(birth.vector[,1]) - vector.matrix <- rbind(vector.matrix, birth.vector) + num_insects.newborn = length(birth.vector[,1]); + vector.matrix = rbind(vector.matrix, birth.vector); # Update population size for the next day. - num_insects <- num_insects - num_insects.death + num_insects.newborn + num_insects = num_insects - num_insects.death + num_insects.newborn; # Aggregate results by day. # Egg population size. - Eggs[row] <- sum(vector.matrix[,2]==0) + Eggs[row] = sum(vector.matrix[,2]==0); # Young nymph population size. - YoungNymphs[row] <- sum(vector.matrix[,2]==1) + YoungNymphs[row] = sum(vector.matrix[,2]==1); # Old nymph population size. - OldNymphs[row] <- sum(vector.matrix[,2]==2) + OldNymphs[row] = sum(vector.matrix[,2]==2); # Previtellogenic population size. - Previtellogenic[row] <- sum(vector.matrix[,2]==3) + Previtellogenic[row] = sum(vector.matrix[,2]==3); # Vitellogenic population size. - Vitellogenic[row] <- sum(vector.matrix[,2]==4) + Vitellogenic[row] = sum(vector.matrix[,2]==4); # Diapausing population size. - Diapausing[row] <- sum(vector.matrix[,2]==5) + Diapausing[row] = sum(vector.matrix[,2]==5); # Newborn population size. - N.newborn[row] <- num_insects.newborn + N.newborn[row] = num_insects.newborn; # Adult population size. - N.adult[row] <- sum(vector.matrix[,2]==3) + sum(vector.matrix[,2]==4) + sum(vector.matrix[,2]==5) + N.adult[row] = sum(vector.matrix[,2]==3) + sum(vector.matrix[,2]==4) + sum(vector.matrix[,2]==5); # Dead population size. - N.death[row] <- num_insects.death + N.death[row] = num_insects.death; - total.population <- c(total.population, num_insects) + total.population = c(total.population, num_insects); # Overwintering adult population size. - overwintering_adult.population[row] <- sum(vector.matrix[,1]==0) + overwintering_adult.population[row] = sum(vector.matrix[,1]==0); # First generation population size. - first_generation.population[row] <- sum(vector.matrix[,1]==1) + first_generation.population[row] = sum(vector.matrix[,1]==1); # Second generation population size. - second_generation.population[row] <- sum(vector.matrix[,1]==2) + second_generation.population[row] = sum(vector.matrix[,1]==2); # P adult population size. - P.adult[row] <- sum(vector.matrix[,1]==0) + P.adult[row] = sum(vector.matrix[,1]==0); # F1 adult population size. - 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)) + 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)); # F2 adult population size - 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)) + 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)); } # End of days specified in the input temperature data. - averages.cum <- cumsum(averages.day) + averages.cum = cumsum(averages.day); # Define the output values. - Eggs.replications[,N.replications] <- Eggs - YoungNymphs.replications[,N.replications] <- YoungNymphs - OldNymphs.replications[,N.replications] <- OldNymphs - Previtellogenic.replications[,N.replications] <- Previtellogenic - Vitellogenic.replications[,N.replications] <- Vitellogenic - Diapausing.replications[,N.replications] <- Diapausing + Eggs.replications[,N.replications] = Eggs; + YoungNymphs.replications[,N.replications] = YoungNymphs; + OldNymphs.replications[,N.replications] = OldNymphs; + Previtellogenic.replications[,N.replications] = Previtellogenic; + Vitellogenic.replications[,N.replications] = Vitellogenic; + Diapausing.replications[,N.replications] = Diapausing; - newborn.replications[,N.replications] <- N.newborn - adult.replications[,N.replications] <- N.adult - death.replications[,N.replications] <- N.death + newborn.replications[,N.replications] = N.newborn; + adult.replications[,N.replications] = N.adult; + 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 - F2_adults.replications[,N.replications] <- F2.adult + 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; + F2_adults.replications[,N.replications] = F2.adult; - population.replications[,N.replications] <- total.population + population.replications[,N.replications] = total.population; } # Mean value for eggs. -eggs <- apply(Eggs.replications, 1, mean) +eggs = apply(Eggs.replications, 1, mean); # Standard error for eggs. -eggs.std_error <- apply(Eggs.replications, 1, sd) / sqrt(opt$replications) +eggs.std_error = apply(Eggs.replications, 1, sd) / sqrt(opt$replications); # Mean value for nymphs. -nymphs <- apply((YoungNymphs.replications+OldNymphs.replications), 1, mean) +nymphs = apply((YoungNymphs.replications+OldNymphs.replications), 1, mean); # Standard error for nymphs. -nymphs.std_error <- apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd) +nymphs.std_error = apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd); # Mean value for adults. -adults <- apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, mean) +adults = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, mean); # Standard error for adults. -adults.std_error <- apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications) +adults.std_error = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications); # Mean value for P. -P <- apply(P.replications, 1, mean) +P = apply(P.replications, 1, mean); # Standard error for P. -P.std_error <- apply(P.replications, 1, sd) / sqrt(opt$replications) +P.std_error = apply(P.replications, 1, sd) / sqrt(opt$replications); # Mean value for P adults. -P_adults <- apply(P_adults.replications, 1, mean) +P_adults = apply(P_adults.replications, 1, mean); # Standard error for P_adult. -P_adults.std_error <- apply(P_adults.replications, 1, sd) / sqrt(opt$replications) +P_adults.std_error = apply(P_adults.replications, 1, sd) / sqrt(opt$replications); # Mean value for F1. -F1 <- apply(F1.replications, 1, mean) +F1 = apply(F1.replications, 1, mean); # Standard error for F1. -F1.std_error <- apply(F1.replications, 1, sd) / sqrt(opt$replications) +F1.std_error = apply(F1.replications, 1, sd) / sqrt(opt$replications); # Mean value for F1 adults. -F1_adults <- apply(F1_adults.replications, 1, mean) +F1_adults = apply(F1_adults.replications, 1, mean); # Standard error for F1 adult. -F1_adults.std_error <- apply(F1_adults.replications, 1, sd) / sqrt(opt$replications) +F1_adults.std_error = apply(F1_adults.replications, 1, sd) / sqrt(opt$replications); # Mean value for F2. -F2 <- apply(F2.replications, 1, mean) +F2 = apply(F2.replications, 1, mean); # Standard error for F2. -F2.std_error <- apply(F2.replications, 1, sd) / sqrt(opt$replications) +F2.std_error = apply(F2.replications, 1, sd) / sqrt(opt$replications); # Mean value for F2 adults. -F2_adults <- apply(F2_adults.replications, 1, mean) +F2_adults = apply(F2_adults.replications, 1, mean); # Standard error for F2 adult. -F2_adults.std_error <- apply(F2_adults.replications, 1, sd) / sqrt(opt$replications) +F2_adults.std_error = apply(F2_adults.replications, 1, sd) / sqrt(opt$replications); # Display the total number of days in the Galaxy history item blurb. -cat("Number of days: ", opt$num_days, "\n") +cat("Number of days: ", opt$num_days, "\n"); -dev.new(width=20, height=30) +dev.new(width=20, height=30); # Start PDF device driver to save charts to output. -pdf(file=opt$output, width=20, height=30, bg="white") -par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)) +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. -days <- c(1:opt$num_days) -start_date <- temperature_data_frame$DATE[1] -end_date <- temperature_data_frame$DATE[opt$num_days] +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. -maxval <- max(eggs+eggs.std_error, nymphs+nymphs.std_error, adults+adults.std_error) +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) + opt$std_error_plot, adults, nymphs, eggs, adults.std_error, nymphs.std_error, eggs.std_error, date_labels); # Subfigure 2: population size by generation. -maxval <- max(F2) +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) + opt$std_error_plot, P, F1, F2, P.std_error, F1.std_error, F2.std_error, date_labels); # Subfigure 3: adult population size by generation. -maxval <- max(F2_adults) + 100 +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) + opt$std_error_plot, P_adults, F1_adults, F2_adults, P_adults.std_error, F1_adults.std_error, F2_adults.std_error, + date_labels); # Turn off device driver to flush output. -dev.off() +dev.off();