insect_phenology_model: insect_phenology

comparison insect_phenology_model.R @ 109:8871797b1be4 draft

Uploaded

author	greg
date	Tue, 05 Dec 2017 15:16:14 -0500
parents	2c2903e7ce68
children	37ac68b6ff10

comparison

equal deleted inserted replaced

-:9d5824e4b5d0
+:8871797b1be4
 #!/usr/bin/env Rscript
 suppressPackageStartupMessages(library("optparse"))
 option_list <- list(
-make_option(c("-a", "--adult_mortality"), action="store", dest="adult_mortality", type="integer", help="Adjustment rate for adult mortality"),
+make_option(c("--adult_mortality"), action="store", dest="adult_mortality", type="integer", help="Adjustment rate for adult mortality"),
-make_option(c("-b", "--adult_accumulation"), action="store", dest="adult_accumulation", type="integer", help="Adjustment of degree-days accumulation (old nymph->adult)"),
+make_option(c("--adult_accumulation"), action="store", dest="adult_accumulation", type="integer", help="Adjustment of degree-days accumulation (old nymph->adult)"),
-make_option(c("-c", "--egg_mortality"), action="store", dest="egg_mortality", type="integer", help="Adjustment rate for egg mortality"),
+make_option(c("--egg_mortality"), action="store", dest="egg_mortality", type="integer", help="Adjustment rate for egg mortality"),
-make_option(c("-e", "--location"), action="store", dest="location", help="Selected location"),
+make_option(c("--input"), action="store", dest="input", help="Temperature data for selected location"),
-make_option(c("-f", "--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"),
+make_option(c("--insect"), action="store", dest="insect", help="Insect name"),
-make_option(c("-i", "--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),
+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("-j", "--nymph_mortality"), action="store", dest="nymph_mortality", type="integer", help="Adjustment rate for nymph mortality"),
+make_option(c("--location"), action="store", dest="location", help="Selected location"),
-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)"),
+make_option(c("--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"),
-make_option(c("-n", "--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"),
+make_option(c("--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),
-make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"),
+make_option(c("--nymph_mortality"), action="store", dest="nymph_mortality", type="integer", help="Adjustment rate for nymph mortality"),
-make_option(c("-p", "--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"),
+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("-q", "--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"),
+make_option(c("--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"),
-make_option(c("-s", "--replications"), action="store", dest="replications", type="integer", help="Number of replications"),
+make_option(c("--output"), action="store", dest="output", help="Output dataset"),
-make_option(c("-t", "--std_error_plot"), action="store", dest="std_error_plot", help="Plot Standard error"),
+make_option(c("--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"),
-make_option(c("-v", "--input"), action="store", dest="input", help="Temperature data for selected location"),
+make_option(c("--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"),
-make_option(c("-y", "--young_nymph_accumulation"), action="store", dest="young_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)"),
+make_option(c("--replications"), action="store", dest="replications", type="integer", help="Number of replications"),
-make_option(c("-x", "--insect"), action="store", dest="insect", help="Insect name")
+make_option(c("--std_error_plot"), action="store", dest="std_error_plot", 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)")
 )
 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
 args <- parse_args(parser, positional_arguments=TRUE)
 opt <- args$options
 lines(days, group3-group3_std_error, col=4, lty=2)
 }
 }
 temperature_data_frame <- parse_input_data(opt$input, opt$num_days)
-# All latitude values are the same,
+# All latitude values are the same, so get the value from the first row.
-# so get the value from the first row.
 latitude <- temperature_data_frame$LATITUDE[1]
 # 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)
 population.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
 # Process replications.
 for (N.replications in 1:opt$replications) {
-# During each replication start with 1000 individuals.
+# Start with the user-defined number of insects per replication.
-# TODO: user definable as well?
+num_insects <- opt$insects_per_replication
-num_insects <- 1000
 # Generation, Stage, degree-days, T, Diapause.
 vector.ini <- c(0, 3, 0, 0, 0)
 # Overwintering, previttelogenic, degree-days=0, T=0, no-diapause.
 vector.matrix <- rep(vector.ini, num_insects)
 # Complete matrix for the population.
 u.d <- runif(1)
 if (u.d < death.probability) {
 death.vector <- c(death.vector, i)
 }
 else {
-# Aggregrate index of dead bug.
 # End of diapause.
 if (vector.individual[1] == 0 && vector.individual[2] == 3) {
 # Overwintering adult (previttelogenic).
 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
 }
 else {
-# Add to # Add average temperature for current day..
+# Add to # Add average temperature for current day.
 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
 }
 # Add 1 day in current stage.
 vector.individual[4] <- vector.individual[4] + 1
 vector.matrix[i,] <- vector.individual
 }
 }
-# Event 2 oviposition -- where population dynamics comes from.
+# Oviposition -- where population dynamics comes from.
 if (vector.individual[2] == 4 && vector.individual[1] == 0 && mean.temp > 10) {
 # 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))
 new.vector <- t(matrix(new.vector, nrow=5))
 # Group with total eggs laid in that day.
 birth.vector <- rbind(birth.vector, new.vector)
 }
 }
-# Event 2 oviposition -- for generation 1.
+# Oviposition -- for generation 1.
 if (vector.individual[2] == 4 && vector.individual[1] == 1 && mean.temp > 12.5 && doy < 222) {
 # 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))
 new.vector <- t(matrix(new.vector, nrow=5))
 # Group with total eggs laid in that day.
 birth.vector <- rbind(birth.vector, new.vector)
 }
 }
-# Event 3 development (with diapause determination).
+# Egg to young nymph.
-# Event 3.1 egg development to young nymph (vector.individual[2]=0 -> egg).
 if (vector.individual[2] == 0) {
-# Egg stage.
 # Add average temperature for current day.
 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]
 # Add 1 day in current stage.
 vector.individual[4] <- vector.individual[4] + 1
 }
 vector.matrix[i,] <- vector.individual
 }
-# Event 3.2 young nymph to old nymph (vector.individual[2]=1 -> young nymph: determines diapause).
+# Young nymph to old nymph.
 if (vector.individual[2] == 1) {
-# Young nymph stage.
 # Add average temperature for current day.
 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]
 # Add 1 day in current stage.
 vector.individual[4] <- vector.individual[4] + 1
 }
 vector.matrix[i,] <- vector.individual
 }
-# Event 3.3 old nymph to adult: previttelogenic or diapausing?
+# Old nymph to adult: previttelogenic or diapausing?
 if (vector.individual[2] == 2) {
-# Old nymph stage.
 # Add average temperature for current day.
 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]
 if (vector.individual[5] == 0) {
-# Non-diapausing adult -- previttelogenic.
+# Previttelogenic.
 vector.individual <- c(current.gen, 3, 0, 0, 0)
 }
 else {
 # Diapausing.
 vector.individual <- c(current.gen, 5, 0, 0, 1)
 # Add 1 day in current stage.
 vector.individual[4] <- vector.individual[4] + 1
 }
 vector.matrix[i,] <- vector.individual
 }
-# Event 4 growing of diapausing adult (unimportant, but still necessary).
+# 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
 }
 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
-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
 }
 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=40)
+dev.new(width=20, height=30)
 # Start PDF device driver to save charts to output.
-pdf(file=opt$output, width=20, height=40, bg="white")
+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
+# Data analysis and visualization plots only within a single calendar year.
-# 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]
-# Subfigure 1: population size by life stage
+# 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)
-# Subfigure 2: population size by generation
+# 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)
-# Subfigure 3: adult population size by generation
+# 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)
 # Turn off device driver to flush output.

Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 109:8871797b1be4 draft