| Previous changeset 89:1615e60cf61c (2017-12-01) Next changeset 91:b6b42e12e173 (2017-12-01) |
|
Commit message:
Uploaded |
|
modified:
insect_phenology_model.R |
| b |
| diff -r 1615e60cf61c -r 7433448e313d insect_phenology_model.R --- a/insect_phenology_model.R Fri Dec 01 09:20:49 2017 -0500 +++ b/insect_phenology_model.R Fri Dec 01 09:20:56 2017 -0500 |
| b |
| b'@@ -4,21 +4,21 @@\n \n option_list <- list(\n make_option(c("-a", "--adult_mort"), action="store", dest="adult_mort", type="integer", help="Adjustment rate for adult mortality"),\n- make_option(c("-b", "--adult_accum"), action="store", dest="adult_accum", type="integer", help="Adjustment of DD accumulation (old nymph->adult)"),\n+ make_option(c("-b", "--adult_accum"), action="store", dest="adult_accum", type="integer", help="Adjustment of degree-days accumulation (old nymph->adult)"),\n make_option(c("-c", "--egg_mort"), action="store", dest="egg_mort", type="integer", help="Adjustment rate for egg mortality"),\n make_option(c("-e", "--location"), action="store", dest="location", help="Selected location"),\n make_option(c("-f", "--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"),\n make_option(c("-i", "--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),\n make_option(c("-j", "--nymph_mort"), action="store", dest="nymph_mort", type="integer", help="Adjustment rate for nymph mortality"),\n- make_option(c("-k", "--old_nymph_accum"), action="store", dest="old_nymph_accum", type="integer", help="Adjustment of DD accumulation (young nymph->old nymph)"),\n+ 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)"),\n make_option(c("-n", "--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"),\n make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"),\n make_option(c("-p", "--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"),\n make_option(c("-q", "--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"),\n make_option(c("-s", "--replications"), action="store", dest="replications", type="integer", help="Number of replications"),\n- make_option(c("-t", "--se_plot"), action="store", dest="se_plot", help="Plot SE"),\n+ make_option(c("-t", "--std_error_plot"), action="store", dest="std_error_plot", help="Plot Standard error"),\n make_option(c("-v", "--input"), action="store", dest="input", help="Temperature data for selected location"),\n- make_option(c("-y", "--young_nymph_accum"), action="store", dest="young_nymph_accum", type="integer", help="Adjustment of DD accumulation (egg->young nymph)"),\n+ make_option(c("-y", "--young_nymph_accum"), action="store", dest="young_nymph_accum", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)"),\n make_option(c("-x", "--insect"), action="store", dest="insect", help="Insect name")\n )\n \n@@ -77,9 +77,9 @@\n # Mean temperature for current row.\n dmean <- 0.5 * (dnp + dxp)\n # Initialize degree day accumulation\n- dd <- 0\n+ degree_days <- 0\n if (dxp < threshold) {\n- dd <- 0\n+ degree_days <- 0\n }\n else {\n # Initialize hourly temperature.\n@@ -132,9 +132,9 @@\n }\n }\n }\n- dd <- sum(dh) / 24\n+ degree_days <- sum(dh) / 24\n }\n- return(c(dmean, dd))\n+ return(c(dmean, degree_days))\n }\n \n dev.egg = function(temperature) {\n@@ -209,10 +209,10 @@\n for (N.rep in 1:opt$replications) {\n # During each replication start with 1000 individuals.\n # TODO: user definable as well?\n- n <- 1000\n- # Generation, Stage, DD, T, Diapause.\n+ num_insects <- 1000\n+ # Generation, Stage, degree-days, T, Diapause.\n vec.ini <- c(0, 3, 0, 0, 0)\n- # Overwintering, previttelogenic, DD=0, T=0, no-diapause.\n+ # Overwintering, previttelogenic, degree-days=0, T=0, no-diapause.\n vec.mat <- rep(vec.ini, n)\n # Complete matrix for the population.\n vec.mat <'..b' latitude, ":", start_date, "-", end_date, sep=" ")\n-plot(day.all, sa, main=title, type="l", ylim=c(0, max(se + se.se, sn + sn.se, sa + sa.se)), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)\n+plot(day.all, mean_value_adult, main=title, type="l", ylim=c(0, max(mean_value_egg + mean_value_egg.se, mean_value_nymph + mean_value_nymph.se, mean_value_adult + mean_value_adult.se)), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)\n # Young and old nymphs.\n-lines(day.all, sn, lwd=2, lty=1, col=2)\n+lines(day.all, mean_value_nymph, lwd=2, lty=1, col=2)\n # Eggs\n-lines(day.all, se, lwd=2, lty=1, col=4)\n+lines(day.all, mean_value_egg, lwd=2, lty=1, col=4)\n 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"))\n axis(2, cex.axis=3)\n leg.text <- c("Egg", "Nymph", "Adult")\n legend("topleft", leg.text, lty=c(1, 1, 1), col=c(4, 2, 1), cex=3)\n if (opt$se_plot == 1) {\n- # Add SE lines to plot\n- # SE for adults\n- lines (day.all, sa + sa.se, lty=2)\n- lines (day.all, sa - sa.se, lty=2) \n- # SE for nymphs\n- lines (day.all, sn + sn.se, col=2, lty=2)\n- lines (day.all, sn - sn.se, col=2, lty=2)\n- # SE for eggs\n- lines (day.all, se + se.se, col=4, lty=2)\n- lines (day.all, se - se.se, col=4, lty=2)\n+ # Add Standard error lines to plot\n+ # Standard error for adults\n+ lines (day.all, mean_value_adult+mean_value_adult.se, lty=2)\n+ lines (day.all, mean_value_adult-mean_value_adult.se, lty=2) \n+ # Standard error for nymphs\n+ lines (day.all, mean_value_nymph+mean_value_nymph.se, col=2, lty=2)\n+ lines (day.all, mean_value_nymph-mean_value_nymph.se, col=2, lty=2)\n+ # Standard error for eggs\n+ lines (day.all, mean_value_egg+mean_value_egg.se, col=4, lty=2)\n+ lines (day.all, mean_value_egg-mean_value_egg.se, col=4, lty=2)\n }\n \n # Subfigure 2: population size by generation\n@@ -609,16 +606,16 @@\n leg.text <- c("P", "F1", "F2")\n legend("topleft", leg.text, lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)\n if (opt$se_plot == 1) {\n- # Add SE lines to plot\n- # SE for adults\n- lines (day.all, g0+g0.se, lty=2)\n- lines (day.all, g0-g0.se, lty=2) \n- # SE for nymphs\n- lines (day.all, g1+g1.se, col=2, lty=2)\n- lines (day.all, g1-g1.se, col=2, lty=2)\n- # SE for eggs\n- lines (day.all, g2+g2.se, col=4, lty=2)\n- lines (day.all, g2-g2.se, col=4, lty=2)\n+ # Add Standard error lines to plot\n+ # Standard error for adults\n+ lines (day.all, g0+g0.std_error, lty=2)\n+ lines (day.all, g0-g0.std_error, lty=2) \n+ # Standard error for nymphs\n+ lines (day.all, g1+g1.std_error, col=2, lty=2)\n+ lines (day.all, g1-g1.std_error, col=2, lty=2)\n+ # Standard error for eggs\n+ lines (day.all, g2+g2.std_error, col=4, lty=2)\n+ lines (day.all, g2-g2.std_error, col=4, lty=2)\n }\n \n # Subfigure 3: adult population size by generation\n@@ -630,17 +627,17 @@\n axis(2, cex.axis=3)\n leg.text <- c("P", "F1", "F2")\n legend("topleft", leg.text, lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)\n-if (opt$se_plot == 1) {\n- # Add SE lines to plot\n- # SE for adults\n- lines (day.all, g0a+g0a.se, lty=2)\n- lines (day.all, g0a-g0a.se, lty=2) \n- # SE for nymphs\n- lines (day.all, g1a+g1a.se, col=2, lty=2)\n- lines (day.all, g1a-g1a.se, col=2, lty=2)\n- # SE for eggs\n- lines (day.all, g2a+g2a.se, col=4, lty=2)\n- lines (day.all, g2a-g2a.se, col=4, lty=2)\n+if (opt$std_error_plot == 1) {\n+ # Add Standard error lines to plot\n+ # Standard error for adults\n+ lines (day.all, g0a+g0a.std_error, lty=2)\n+ lines (day.all, g0a-g0a.std_error, lty=2) \n+ # Standard error for nymphs\n+ lines (day.all, g1a+g1a.std_error, col=2, lty=2)\n+ lines (day.all, g1a-g1a.std_error, col=2, lty=2)\n+ # Standard error for eggs\n+ lines (day.all, g2a+g2a.std_error, col=4, lty=2)\n+ lines (day.all, g2a-g2a.std_error, col=4, lty=2)\n }\n \n # Turn off device driver to flush output.\n' |