insect_phenology_model: insect_phenology

comparison insect_phenology_model.R @ 103:4d9e40d8af0f draft

Uploaded

author	greg
date	Tue, 05 Dec 2017 11:01:23 -0500
parents	ce996e90f5a7
children	d6e37828b094

comparison

equal deleted inserted replaced

-:ce996e90f5a7
+:4d9e40d8af0f
 # Maximum temperature for current 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)
 # Initialize degree day accumulation
-degree_days <- 0
+averages <- 0
 if (curr_max_temp < threshold) {
-degree_days <- 0
+averages <- 0
 }
 else {
 # Initialize hourly temperature.
 T <- NULL
 # Initialize degree hour vector.
 else {
 dh[i] <- T[i] - 8.4
 }
 }
 }
-degree_days <- sum(dh) / 24
+averages <- sum(dh) / 24
 }
-return(c(curr_mean_temp, degree_days))
+return(c(curr_mean_temp, averages))
 }
 mortality.adult = function(temperature) {
 if (temperature < 12.7) {
 mortality.probability = 0.002
 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) {
 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("topleft", c("Egg", "Nymph", "Adult"), lty=c(1, 1, 1), col=c(4, 2, 1), cex=3)
+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("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
+legend_text <- c("P", "F1", "F2")
+columns <- col=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("topleft", c("P", "F1", "F2"), lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
+legend_text <- c("P", "F1", "F2")
+columns <- col=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), 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)
 if (plot_std_error==1) {
 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)
-adult.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)
-P.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)
-F2.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_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)
 # Process replications.
 for (N.replications in 1:opt$replications) {
 # During each replication start with 1000 individuals.
 # TODO: user definable as well?
 # Overwintering, previttelogenic, degree-days=0, T=0, no-diapause.
 vector.matrix <- rep(vector.ini, num_insects)
 # Complete matrix for the population.
 vector.matrix <- base::t(matrix(vector.matrix, nrow=5))
 # Time series of population size.
-total.population <- NULL
-overwintering_adult.population <- rep(0, opt$num_days)
-first_generation.population <- rep(0, opt$num_days)
-second_generation.population <- 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)
+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)
-N.newborn <- rep(0, opt$num_days)
-N.death <- rep(0, opt$num_days)
+total.population <- NULL
-N.adult <- rep(0, opt$num_days)
-degree_days.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]
 # 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]
-degree_days.temp <- temp.profile[2]
+averages.temp <- temp.profile[2]
-degree_days.day[row] <- degree_days.temp
+averages.day[row] <- averages.temp
 # Trash bin for death.
 death.vector <- NULL
 # Newborn.
 birth.vector <- NULL
 # All individuals.
 for (i in 1:num_insects) {
-# Find individual record.
+# Individual record.
 vector.individual <- vector.matrix[i,]
-# First of all, still alive?
+# Adjustment for late season mortality rate (still alive?).
-# Adjustment for late season mortality rate.
 if (latitude < 40.0) {
 post.mortality <- 1
 day.kill <- 300
 }
 else {
 }
 else if (vector.individual[2] == 1 | vector.individual[2] == 2) {
 death.probability = opt$nymph_mortality * mortality.nymph(mean.temp)
 }
 else if (vector.individual[2] == 3 | vector.individual[2] == 4 | vector.individual[2] == 5) {
-# For adult.
+# Adult.
 if (doy < day.kill) {
 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)
 }
 }
-# (or dependent on temperature and life stage?)
+# Dependent on temperature and life stage?
 u.d <- runif(1)
 if (u.d < death.probability) {
 death.vector <- c(death.vector, i)
 }
 else {
 # Aggregrate index of dead bug.
-# Event 1 end of diapause.
+# 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 degree_days.
+# Add to # Add average temperature for current day..
-vector.individual[3] <- vector.individual[3] + degree_days.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
 }
 }
 # Transfer to vittelogenic.
 vector.individual <- c(current.gen, 4, 0, 0, 0)
 vector.matrix[i,] <- vector.individual
 }
 else {
-# Add to degree_days.
+# Add average temperature for current day.
-vector.individual[3] <- vector.individual[3] + degree_days.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
 }
 }
 u1 <- runif(1)
 if (u1 < p.birth) {
 num_insects.birth = round(runif(1, 2, 8))
 }
 }
-# Add to degree_days.
+# Add average temperature for current day.
-vector.individual[3] <- vector.individual[3] + degree_days.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
 if (num_insects.birth > 0) {
 # Add new birth -- might be in different generations.
 u1 <- runif(1)
 if (u1 < p.birth) {
 num_insects.birth = round(runif(1, 2, 8))
 }
 }
-# Add to degree_days.
+# Add average temperature for current day.
-vector.individual[3] <- vector.individual[3] + degree_days.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
 if (num_insects.birth > 0) {
 # Add new birth -- might be in different generations.
 }
 # Event 3 development (with diapause determination).
 # Event 3.1 egg development to young nymph (vector.individual[2]=0 -> egg).
 if (vector.individual[2] == 0) {
 # Egg stage.
-# Add to degree_days.
+# Add average temperature for current day.
-vector.individual[3] <- vector.individual[3] + degree_days.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]
 # Transfer to young nymph stage.
 vector.individual <- c(current.gen, 1, 0, 0, 0)
 vector.matrix[i,] <- vector.individual
 }
 # Event 3.2 young nymph to old nymph (vector.individual[2]=1 -> young nymph: determines diapause).
 if (vector.individual[2] == 1) {
 # Young nymph stage.
-# Add to degree_days.
+# Add average temperature for current day.
-vector.individual[3] <- vector.individual[3] + degree_days.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]
 # Transfer to old nym stage.
 vector.individual <- c(current.gen, 2, 0, 0, 0)
 vector.matrix[i,] <- vector.individual
 }
 # Event 3.3 old nymph to adult: previttelogenic or diapausing?
 if (vector.individual[2] == 2) {
 # Old nymph stage.
-# Add to degree_days.
+# Add average temperature for current day.
-vector.individual[3] <- vector.individual[3] + degree_days.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]
 if (vector.individual[5] == 0) {
 # Non-diapausing adult -- previttelogenic.
 }
 vector.matrix[i,] <- vector.individual
 }
 # Event 4 growing of diapausing adult (unimportant, but still necessary).
 if (vector.individual[2] == 5) {
-vector.individual[3] <- vector.individual[3] + degree_days.temp
+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.
 N.death[row] <- num_insects.death
 # Adult population size.
 N.adult[row] <- num_insects.adult
 }   # End of days specified in the input temperature data.
-degree_days.cum <- cumsum(degree_days.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
 P_adults.replications[,N.replications] <- P.adult
 F1_adults.replications[,N.replications] <- F1.adult
 F2_adults.replications[,N.replications] <- F2.adult
 }
+# 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)
-# Mean value for nymphs.
+# Standard error for adults.
-nymphs <- apply((YoungNymphs.replications+OldNymphs.replications), 1, mean)
+adults.std_error <- apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications)
-# Mean value for eggs.
-eggs <- apply(Eggs.replications, 1, mean)
 # 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)
+# 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)
-# Mean value for P adults.
+# Standard error for F2.
-P_adults <- apply(P_adults.replications, 1, mean)
+F2.std_error <- apply(F2.replications, 1, sd) / sqrt(opt$replications)
-# Mean value for F1 adults.
-F1_adults <- apply(F1_adults.replications, 1, mean)
 # Mean value for F2 adults.
 F2_adults <- apply(F2_adults.replications, 1, mean)
-# Standard error for adults.
-adults.std_error <- apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications)
-# Standard error for nymphs.
-nymphs.std_error <- apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd)
-# Standard error for eggs.
-eggs.std_error <- apply(Eggs.replications, 1, sd) / sqrt(opt$replications)
-# Standard error value for P.
-P.std_error <- apply(P.replications, 1, sd) / sqrt(opt$replications)
-# Standard error for F1.
-F1.std_error <- apply(F1.replications, 1, sd) / sqrt(opt$replications)
-# Standard error for F2.
-F2.std_error <- apply(F2.replications, 1, sd) / sqrt(opt$replications)
-# Standard error for P adult.
-P_adults.std_error <- apply(P_adults.replications, 1, sd) / sqrt(opt$replications)
-# Standard error for F1 adult.
-F1_adults.std_error <- apply(F1_adults.replications, 1, sd) / sqrt(opt$replications)
 # Standard error for F2 adult.
 F2_adults.std_error <- apply(F2_adults.replications, 1, sd) / sqrt(opt$replications)
 dev.new(width=20, height=40)

Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 103:4d9e40d8af0f draft