Mercurial > repos > greg > insect_phenology_model
changeset 10:61bc6bd8807d draft
Uploaded
| author | greg |
|---|---|
| date | Tue, 27 Feb 2018 13:31:59 -0500 |
| parents | d9371485aaf5 |
| children | d1f2c0634354 |
| files | insect_phenology_model.R |
| diffstat | 1 files changed, 523 insertions(+), 168 deletions(-) [+] |
line wrap: on
line diff
--- a/insect_phenology_model.R Tue Feb 27 13:31:49 2018 -0500 +++ b/insect_phenology_model.R Tue Feb 27 13:31:59 2018 -0500 @@ -9,17 +9,20 @@ make_option(c("--input"), action="store", dest="input", help="Temperature data for selected location"), make_option(c("--insect"), action="store", dest="insect", help="Insect name"), make_option(c("--insects_per_replication"), action="store", dest="insects_per_replication", type="integer", help="Number of insects with which to start each replication"), + make_option(c("--life_stages"), action="store", dest="life_stages", help="Selected life stages for plotting"), + make_option(c("--life_stages_adult"), action="store", dest="life_stages_adult", default=NULL, help="Adult life stages for plotting"), + make_option(c("--life_stage_nymph"), action="store", dest="life_stage_nymph", default=NULL, help="Nymph life stages for plotting"), make_option(c("--location"), action="store", dest="location", help="Selected location"), make_option(c("--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"), make_option(c("--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"), make_option(c("--nymph_mortality"), action="store", dest="nymph_mortality", type="integer", help="Adjustment rate for nymph mortality"), make_option(c("--old_nymph_accumulation"), action="store", dest="old_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (young nymph->old nymph)"), make_option(c("--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"), - make_option(c("--output"), action="store", dest="output", help="Output dataset"), make_option(c("--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"), make_option(c("--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"), make_option(c("--replications"), action="store", dest="replications", type="integer", help="Number of replications"), - make_option(c("--std_error_plot"), action="store", dest="std_error_plot", help="Plot Standard error"), + make_option(c("--plot_generations_separately"), action="store", dest="plot_generations_separately", help="Plot Plot P, F1 and F2 as separate lines or pool across them"), + make_option(c("--plot_std_error"), action="store", dest="plot_std_error", help="Plot Standard error"), make_option(c("--young_nymph_accumulation"), action="store", dest="young_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)") ) @@ -73,7 +76,6 @@ 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(); @@ -215,65 +217,170 @@ } -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, 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); - } 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); - } else if (chart_type == "adult_pop_size_by_generation") { - title = paste(insect, ": Adult pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" "); +render_chart = function(date_labels, chart_type, plot_std_error, insect, location, latitude, start_date, end_date, days, maxval, + replications, life_stage, group, group_std_error, group2=NULL, group2_std_error=NULL, group3=NULL, group3_std_error=NULL, + life_stages_adult=NULL, life_stage_nymph=NULL) { + if (chart_type=="pop_size_by_life_stage") { + if (life_stage=="Total") { + title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); + legend_text = c("Egg", "Nymph", "Adult"); + columns = c(4, 2, 1); + plot(days, group, 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=="yes") { + # Standard error for group. + lines(days, group+group_std_error, lty=2); + lines(days, group-group_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); + # Standard error for group3. + lines(days, group3+group3_std_error, col=4, lty=2); + lines(days, group3-group3_std_error, col=4, lty=2); + } + } else { + if (life_stage=="Egg") { + title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); + legend_text = c(life_stage); + columns = c(4); + } else if (life_stage=="Nymph") { + stage = paste(life_stage_nymph, "Nymph Pop :", sep=" "); + title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); + legend_text = c(paste(life_stage_nymph, life_stage, sep=" ")); + columns = c(2); + } else if (life_stage=="Adult") { + stage = paste(life_stages_adult, "Adult Pop", sep=" "); + title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); + legend_text = c(paste(life_stages_adult, life_stage, sep=" ")); + columns = c(1); + } + plot(days, group, 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), col="black", cex=3); + 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=="yes") { + # Standard error for group. + lines(days, group+group_std_error, lty=2); + lines(days, group-group_std_error, lty=2); + } + } + } else if (chart_type=="pop_size_by_generation") { + if (life_stage=="Total") { + title_str = ": Total Pop by Gen :"; + } else if (life_stage=="Egg") { + title_str = ": Egg Pop by Gen :"; + } else if (life_stage=="Nymph") { + title_str = paste(":", life_stage_nymph, "Nymph Pop by Gen", ":", sep=" "); + } else if (life_stage=="Adult") { + title_str = paste(":", life_stages_adult, "Adult Pop by Gen", ":", sep=" "); + } + title = paste(insect, ": Reps", replications, title_str, 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: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); - # Standard error for group2. - 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); + plot(days, group, 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=="yes") { + # Standard error for group. + lines(days, group+group_std_error, lty=2); + lines(days, group-group_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); + # Standard error for group3. + lines(days, group3+group3_std_error, col=4, lty=2); + lines(days, group3-group3_std_error, col=4, lty=2); + } } } +# Determine if we're plotting generations separately. +if (opt$plot_generations_separately=="yes") { + plot_generations_separately = TRUE; +} else { + plot_generations_separately = FALSE; +} +# Read the temperature data into a data frame. 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. +output_dir = "output_dir"; +# Get the date labels for plots. +date_labels = get_date_labels(temperature_data_frame, opt$num_days); +# All latitude values are the same, so get the value for plots from the first row. latitude = temperature_data_frame$LATITUDE[1]; +# Get the number of days for plots. num_columns = dim(temperature_data_frame)[2]; -date_labels = get_date_labels(temperature_data_frame, opt$num_days); +# Split life_stages into a list of strings for plots. +life_stages_str = as.character(opt$life_stages); +life_stages = strsplit(life_stages_str, ",")[[1]]; +# Determine the data we need to generate for plotting. +process_eggs = FALSE; +process_nymphs = FALSE; +process_adults = FALSE; +for (life_stage in life_stages) { + if (life_stage=="Total") { + process_eggs = TRUE; + process_nymphs = TRUE; + life_stage_nymph = "Total"; + process_adults = TRUE; + life_stages_adult = "Total"; + } else if (life_stage=="Egg") { + process_eggs = TRUE; + } else if (life_stage=="Nymph") { + process_nymphs = TRUE; + life_stage_nymph = opt$life_stage_nymph; + } else if (life_stage=="Adult") { + process_adults = TRUE; + life_stages_adult = opt$life_stages_adult; + } +} # 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); - +if (process_eggs) { + Eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); +} +if (process_nymphs) { + 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); +} +if (process_adults) { + 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); - -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); - +if (plot_generations_separately) { + # P is Parental, or overwintered adults. + P.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); + # F1 is the first field-produced generation. + F1.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); + # F2 is the second field-produced generation. + F2.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); + if (process_eggs) { + P_eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); + F1_eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); + F2_eggs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); + } + if (process_nymphs) { + P_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); + F1_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); + F2_nymphs.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); + } + if (process_adults) { + P_adults.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_adults.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); + } +} +# Total population. population.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); # Process replications. @@ -282,32 +389,53 @@ num_insects = opt$insects_per_replication; # Generation, Stage, degree-days, T, Diapause. vector.ini = c(0, 3, 0, 0, 0); - # Overwintering, previttelogenic, degree-days=0, T=0, no-diapause. + # Replicate to create a matrix where the columns are + # Generation, Stage, degree-days, T, Diapause and the + # rows are the initial number of insects per replication. vector.matrix = rep(vector.ini, num_insects); - # Complete matrix for the population. + # Complete transposed matrix for the population, so now + # the rows are Generation, Stage, degree-days, T, Diapause 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); - + if (process_eggs) { + Eggs = rep(0, opt$num_days); + } + if (process_nymphs) { + YoungNymphs = rep(0, opt$num_days); + OldNymphs = rep(0, opt$num_days); + } + if (process_adults) { + 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); - 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); - + if (plot_generations_separately) { + # P is Parental, or overwintered adults. + # F1 is the first field-produced generation. + # F2 is the second field-produced generation. + if (process_eggs) { + P.egg = rep(0, opt$num_days); + F1.egg = rep(0, opt$num_days); + F2.egg = rep(0, opt$num_days); + } + if (process_nymphs) { + P.nymph = rep(0, opt$num_days); + F1.nymph = rep(0, opt$num_days); + F2.nymph = rep(0, opt$num_days); + } + if (process_adults) { + P.adult = rep(0, opt$num_days); + F1.adult = rep(0, opt$num_days); + F2.adult = rep(0, opt$num_days); + } + } total.population = NULL; - averages.day = rep(0, opt$num_days); # All the days included in the input temperature dataset. for (row in 1:opt$num_days) { @@ -535,19 +663,27 @@ # Update population size for the next day. num_insects = num_insects - num_insects.death + num_insects.newborn; - # Aggregate results by day. - # Egg population size. - Eggs[row] = sum(vector.matrix[,2]==0); - # Young nymph population size. - YoungNymphs[row] = sum(vector.matrix[,2]==1); - # Old nymph population size. - OldNymphs[row] = sum(vector.matrix[,2]==2); - # Previtellogenic population size. - Previtellogenic[row] = sum(vector.matrix[,2]==3); - # Vitellogenic population size. - Vitellogenic[row] = sum(vector.matrix[,2]==4); - # Diapausing population size. - Diapausing[row] = sum(vector.matrix[,2]==5); + # Aggregate results by day. Due to multiple transpose calls + # on vector.matrix above, the columns of vector.matrix + # are now Generation, Stage, degree-days, T, Diapause, + if (process_eggs) { + # For egg population size, column 2 (Stage), must be 0. + Eggs[row] = sum(vector.matrix[,2]==0); + } + if (process_nymphs) { + # For young nymph population size, column 2 (Stage) must be 1. + YoungNymphs[row] = sum(vector.matrix[,2]==1); + # For old nymph population size, column 2 (Stage) must be 2. + OldNymphs[row] = sum(vector.matrix[,2]==2); + } + if (process_adults) { + # For pre-vitellogenic population size, column 2 (Stage) must be 3. + Previtellogenic[row] = sum(vector.matrix[,2]==3); + # For vitellogenic population size, column 2 (Stage) must be 4. + Vitellogenic[row] = sum(vector.matrix[,2]==4); + # For diapausing population size, column 2 (Stage) must be 5. + Diapausing[row] = sum(vector.matrix[,2]==5); + } # Newborn population size. N.newborn[row] = num_insects.newborn; @@ -558,117 +694,336 @@ total.population = c(total.population, num_insects); - # Overwintering adult population size. + # For overwintering adult (P) population + # size, column 1 (Generation) must be 0. overwintering_adult.population[row] = sum(vector.matrix[,1]==0); - # First generation population size. + # For first field generation (F1) population + # size, column 1 (Generation) must be 1. first_generation.population[row] = sum(vector.matrix[,1]==1); - # Second generation population size. + # For second field generation (F2) population + # size, column 1 (Generation) must be 2. second_generation.population[row] = sum(vector.matrix[,1]==2); - # P adult population size. - 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)); - # 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)); + if (plot_generations_separately) { + if (process_eggs) { + # For egg life stage of generation F1 population size, + # column 1 (generation) is 0 and column 2 (Stage) is 0. + P.egg[row] = sum(vector.matrix[,1]==0 & vector.matrix[,2]==0); + # For egg life stage of generation F1 population size, + # column 1 (generation) is 1 and column 2 (Stage) is 0. + F1.egg[row] = sum(vector.matrix[,1]==1 & vector.matrix[,2]==0); + # For egg life stage of generation F2 population size, + # column 1 (generation) is 2 and column 2 (Stage) is 0. + F2.egg[row] = sum(vector.matrix[,1]==2 & vector.matrix[,2]==0); + } + if (process_nymphs) { + # For nymph life stage of generation F1 population + # size, one of the following combinations is required: + # - column 1 (Generation) is 0 and column 2 (Stage) is 1 (Young nymph) + # - column 1 (Generation) is 0 and column 2 (Stage) is 2 (Old nymph) + P.nymph[row] = sum((vector.matrix[,1]==0 & vector.matrix[,2]==1) | (vector.matrix[,1]==0 & vector.matrix[,2]==2)); + # For nymph life stage of generation F1 population + # size, one of the following combinations is required: + # - column 1 (Generation) is 1 and column 2 (Stage) is 1 (Young nymph) + # - column 1 (Generation) is 1 and column 2 (Stage) is 2 (Old nymph) + F1.nymph[row] = sum((vector.matrix[,1]==1 & vector.matrix[,2]==1) | (vector.matrix[,1]==1 & vector.matrix[,2]==2)); + # For nymph life stage of generation F2 population + # size, one of the following combinations is required: + # - column 1 (Generation) is 2 and column 2 (Stage) is 1 (Young nymph) + # - column 1 (Generation) is 2 and column 2 (Stage) is 2 (Old nymph) + F2.nymph[row] = sum((vector.matrix[,1]==2 & vector.matrix[,2]==1) | (vector.matrix[,1]==2 & vector.matrix[,2]==2)); + } + if (process_adults) { + # For adult life stage of generation P population + # size, one of the following combinations is required: + # - column 1 (Generation) is 0 and column 2 (Stage) is 3 (Pre-vitellogenic) + # - column 1 (Generation) is 0 and column 2 (Stage) is 4 (Vitellogenic) + # - column 1 (Generation) is 0 and column 2 (Stage) is 5 (Diapausing) + P.adult[row] = sum((vector.matrix[,1]==0 & vector.matrix[,2]==3) | (vector.matrix[,1]==0 & vector.matrix[,2]==4) | (vector.matrix[,1]==0 & vector.matrix[,2]==5)); + # For adult life stage of generation F1 population + # size, one of the following combinations is required: + # - column 1 (Generation) is 1 and column 2 (Stage) is 3 (Pre-vitellogenic) + # - column 1 (Generation) is 1 and column 2 (Stage) is 4 (Vitellogenic) + # - column 1 (Generation) is 1 and column 2 (Stage) is 5 (Diapausing) + 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)); + # For adult life stage of generation F2 population + # size, one of the following combinations is required: + # - column 1 (Generation) is 2 and column 2 (Stage) is 3 (Pre-vitellogenic) + # - column 1 (Generation) is 2 and column 2 (Stage) is 4 (Vitellogenic) + # - column 1 (Generation) is 2 and column 2 (Stage) is 5 (Diapausing) + 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); # 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; - + if (process_eggs) { + Eggs.replications[,N.replications] = Eggs; + } + if (process_nymphs) { + YoungNymphs.replications[,N.replications] = YoungNymphs; + OldNymphs.replications[,N.replications] = OldNymphs; + } + if (process_adults) { + 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; - - 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; - + if (plot_generations_separately) { + # P is Parental, or overwintered adults. + P.replications[,N.replications] = overwintering_adult.population; + # F1 is the first field-produced generation. + F1.replications[,N.replications] = first_generation.population; + # F2 is the second field-produced generation. + F2.replications[,N.replications] = second_generation.population; + if (process_eggs) { + P_eggs.replications[,N.replications] = P.egg; + F1_eggs.replications[,N.replications] = F1.egg; + F2_eggs.replications[,N.replications] = F2.egg; + } + if (process_nymphs) { + P_nymphs.replications[,N.replications] = P.nymph; + F1_nymphs.replications[,N.replications] = F1.nymph; + F2_nymphs.replications[,N.replications] = F2.nymph; + } + if (process_adults) { + P_adults.replications[,N.replications] = P.adult; + F1_adults.replications[,N.replications] = F1.adult; + F2_adults.replications[,N.replications] = F2.adult; + } + } population.replications[,N.replications] = total.population; } -# Mean value for eggs. -eggs = apply(Eggs.replications, 1, mean); -# Standard error for eggs. -eggs.std_error = apply(Eggs.replications, 1, sd) / sqrt(opt$replications); - -# Mean value for nymphs. -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); - -# Mean value for adults. -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); - -# Mean value for P. -P = apply(P.replications, 1, mean); -# Standard error for P. -P.std_error = apply(P.replications, 1, sd) / sqrt(opt$replications); +if (process_eggs) { + # Mean value for eggs. + eggs = apply(Eggs.replications, 1, mean); + # Standard error for eggs. + eggs.std_error = apply(Eggs.replications, 1, sd) / sqrt(opt$replications); +} +if (process_nymphs) { + # Calculate nymph populations for selected life stage. + if (life_stage_nymph=="Total") { + # Mean value for all nymphs. + nymphs = apply((YoungNymphs.replications+OldNymphs.replications), 1, mean); + # Standard error for all nymphs. + nymphs.std_error = apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd); + } else if (life_stage_nymph=="Young") { + # Mean value for young nymphs. + nymphs = apply(YoungNymphs.replications, 1, mean); + # Standard error for young nymphs. + nymphs.std_error = apply(YoungNymphs.replications / sqrt(opt$replications), 1, sd); + } else if (life_stage_nymph=="Old") { + # Mean value for old nymphs. + nymphs = apply(OldNymphs.replications, 1, mean); + # Standard error for old nymphs. + nymphs.std_error = apply(OldNymphs.replications / sqrt(opt$replications), 1, sd); + } +} +if (process_adults) { + # Calculate adult populations for selected life stage. + if (life_stages_adult=="Total") { + # Mean value for all adults. + adults = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, mean); + # Standard error for all adults. + adults.std_error = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications); + } else if (life_stages_adult == "Pre-vittelogenic") { + # Mean value for previtellogenic adults. + adults = apply(Previtellogenic.replications, 1, mean); + # Standard error for previtellogenic adults. + adults.std_error = apply(Previtellogenic.replications, 1, sd) / sqrt(opt$replications); + } else if (life_stages_adult == "Vittelogenic") { + # Mean value for vitellogenic adults. + adults = apply(Vitellogenic.replications, 1, mean); + # Standard error for vitellogenic adults. + adults.std_error = apply(Vitellogenic.replications, 1, sd) / sqrt(opt$replications); + } else if (life_stages_adult == "Diapausing") { + # Mean value for vitellogenic adults. + adults = apply(Diapausing.replications, 1, mean); + # Standard error for vitellogenic adults. + adults.std_error = apply(Diapausing.replications, 1, sd) / sqrt(opt$replications); + } +} -# Mean value for P adults. -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); - -# Mean value for F1. -F1 = apply(F1.replications, 1, mean); -# Standard error for F1. -F1.std_error = apply(F1.replications, 1, sd) / sqrt(opt$replications); - -# Mean value for F1 adults. -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); - -# Mean value for F2. -F2 = apply(F2.replications, 1, mean); -# Standard error for F2. -F2.std_error = apply(F2.replications, 1, sd) / sqrt(opt$replications); - -# Mean value for F2 adults. -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); +if (plot_generations_separately) { + # Mean value for P which is Parental, or overwintered adults. + P = apply(P.replications, 1, mean); + # Standard error for P. + P.std_error = apply(P.replications, 1, sd) / sqrt(opt$replications); + # Mean value for F1, which is the first field-produced generation. + F1 = apply(F1.replications, 1, mean); + # Standard error for F1. + F1.std_error = apply(F1.replications, 1, sd) / sqrt(opt$replications); + # Mean value for F2, which is the second field-produced generation. + F2 = apply(F2.replications, 1, mean); + # Standard error for F2. + F2.std_error = apply(F2.replications, 1, sd) / sqrt(opt$replications); + if (process_eggs) { + # Mean value for P eggs. + P_eggs = apply(P_eggs.replications, 1, mean); + # Standard error for P_eggs. + P_eggs.std_error = apply(P_eggs.replications, 1, sd) / sqrt(opt$replications); + # Mean value for F1 eggs. + F1_eggs = apply(F1_eggs.replications, 1, mean); + # Standard error for F1 eggs. + F1_eggs.std_error = apply(F1_eggs.replications, 1, sd) / sqrt(opt$replications); + # Mean value for F2 eggs. + F2_eggs = apply(F2_eggs.replications, 1, mean); + # Standard error for F2 eggs. + F2_eggs.std_error = apply(F2_eggs.replications, 1, sd) / sqrt(opt$replications); + } + if (process_nymphs) { + # Mean value for P nymphs. + P_nymphs = apply(P_nymphs.replications, 1, mean); + # Standard error for P_nymphs. + P_nymphs.std_error = apply(P_nymphs.replications, 1, sd) / sqrt(opt$replications); + # Mean value for F1 nymphs. + F1_nymphs = apply(F1_nymphs.replications, 1, mean); + # Standard error for F1 nymphs. + F1_nymphs.std_error = apply(F1_nymphs.replications, 1, sd) / sqrt(opt$replications); + # Mean value for F2 nymphs. + F2_nymphs = apply(F2_nymphs.replications, 1, mean); + # Standard error for F2 eggs. + F2_nymphs.std_error = apply(F2_nymphs.replications, 1, sd) / sqrt(opt$replications); + } + if (process_adults) { + # Mean value for P adults. + P_adults = apply(P_adults.replications, 1, mean); + # Standard error for P_adults. + P_adults.std_error = apply(P_adults.replications, 1, sd) / sqrt(opt$replications); + # Mean value for F1 adults. + F1_adults = apply(F1_adults.replications, 1, mean); + # Standard error for F1 adults. + F1_adults.std_error = apply(F1_adults.replications, 1, sd) / sqrt(opt$replications); + # Mean value for F2 adults. + F2_adults = apply(F2_adults.replications, 1, mean); + # Standard error for F2 adults. + 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"); -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)); - -# Data analysis and visualization plots only within a single calendar year. +# Information needed for plots plots. 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); -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, date_labels); -# Subfigure 2: population size by generation. -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, date_labels); -# Subfigure 3: adult population size by generation. -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, - date_labels); - -# Turn off device driver to flush output. -dev.off(); +if (plot_generations_separately) { + for (life_stage in life_stages) {} + if (life_stage == "Egg") { + # Start PDF device driver. + dev.new(width=20, height=30); + file_path = paste(output_dir, "egg_pop_by_generation.pdf", sep="/"); + pdf(file=file_path, width=20, height=30, bg="white"); + par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); + # Egg population size by generation. + maxval = max(F2_eggs) + 100; + render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, + opt$replications, life_stage, group=P_eggs, group_std_error=P_eggs.std_error, group2=F1_eggs, group2_std_error=F1_eggs.std_error, group3=F2_eggs, + group3_std_error=F2_eggs.std_error); + # Turn off device driver to flush output. + dev.off(); + } else if (life_stage == "Nymph") { + # Start PDF device driver. + dev.new(width=20, height=30); + file_name = paste(tolower(life_stage_nymph), "nymph_pop_by_generation.pdf", sep="_"); + file_path = paste(output_dir, file_name, sep="/"); + pdf(file=file_path, width=20, height=30, bg="white"); + par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); + # Nymph population size by generation. + maxval = max(F2_nymphs) + 100; + render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, + opt$replications, life_stage, group=P_nymphs, group_std_error=P_nymphs.std_error, group2=F1_nymphs, group2_std_error=F1_nymphs.std_error, + group3=F2_nymphs, group3_std_error=F2_nymphs.std_error, life_stage_nymph=life_stage_nymph); + # Turn off device driver to flush output. + dev.off(); + } else if (life_stage == "Adult") { + # Start PDF device driver. + dev.new(width=20, height=30); + file_name = paste(tolower(life_stages_adult), "adult_pop_by_generation.pdf", sep="_"); + file_path = paste(output_dir, file_name, sep="/"); + pdf(file=file_path, width=20, height=30, bg="white"); + par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); + # Adult population size by generation. + maxval = max(F2_adults) + 100; + render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, + opt$replications, life_stage, group=P_adults, group_std_error=P_adults.std_error, group2=F1_adults, group2_std_error=F1_adults.std_error, + group3=F2_adults, group3_std_error=F2_adults.std_error, life_stages_adult=life_stages_adult); + # Turn off device driver to flush output. + dev.off(); + } else if (life_stage == "Total") { + # Start PDF device driver. + dev.new(width=20, height=30); + file_path = paste(output_dir, "total_pop_by_generation.pdf", sep="/"); + pdf(file=file_path, width=20, height=30, bg="white"); + par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); + # Total population size by generation. + maxval = max(F2); + render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, + opt$replications, life_stage, group=P, group_std_error=P.std_error, group2=F1, group2_std_error=F1.std_error, group3=F2, group3_std_error=F2.std_error); + # Turn off device driver to flush output. + dev.off(); + } +} else { + for (life_stage in life_stages) { + if (life_stage == "Egg") { + # Start PDF device driver. + dev.new(width=20, height=30); + file_path = paste(output_dir, "egg_pop_size.pdf", sep="/"); + pdf(file=file_path, width=20, height=30, bg="white"); + par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); + # Egg population size. + maxval = max(eggs+eggs.std_error); + render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, + opt$replications, life_stage, group=eggs, group_std_error=eggs.std_error); + # Turn off device driver to flush output. + dev.off(); + } else if (life_stage == "Nymph") { + # Start PDF device driver. + dev.new(width=20, height=30); + file_name = paste(tolower(life_stage_nymph), "nymph_pop_size.pdf", sep="_"); + file_path = paste(output_dir, file_name, sep="/"); + pdf(file=file_path, width=20, height=30, bg="white"); + par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); + # Nymph population size. + maxval = max(nymphs+nymphs.std_error); + render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, + opt$replications, life_stage, group=nymphs, group_std_error=nymphs.std_error, life_stage_nymph=life_stage_nymph); + # Turn off device driver to flush output. + dev.off(); + } else if (life_stage == "Adult") { + # Start PDF device driver. + dev.new(width=20, height=30); + file_name = paste(tolower(life_stages_adult), "adult_pop_size.pdf", sep="_"); + file_path = paste(output_dir, file_name, sep="/"); + pdf(file=file_path, width=20, height=30, bg="white"); + par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); + # Adult population size. + maxval = max(adults+adults.std_error); + render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, + opt$replications, life_stage, group=adults, group_std_error=adults.std_error, life_stages_adult=life_stages_adult); + # Turn off device driver to flush output. + dev.off(); + } else if (life_stage == "Total") { + # Start PDF device driver. + dev.new(width=20, height=30); + file_path = paste(output_dir, "total_pop_size.pdf", sep="/"); + pdf(file=file_path, width=20, height=30, bg="white"); + par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); + # Total population size. + maxval = max(eggs+eggs.std_error, nymphs+nymphs.std_error, adults+adults.std_error); + render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, + opt$replications, life_stage, group=adults, group_std_error=adults.std_error, group2=nymphs, group2_std_error=nymphs.std_error, group3=eggs, + group3_std_error=eggs.std_error); + # Turn off device driver to flush output. + dev.off(); + } + } +}