| Previous changeset 25:3852133fb058 (2017-11-09) Next changeset 27:24949e72f7ec (2017-11-09) |
|
Commit message:
Uploaded |
|
modified:
insect_phenology_model.R |
| b |
| diff -r 3852133fb058 -r 3c06cab3db2c insect_phenology_model.R --- a/insect_phenology_model.R Thu Nov 09 11:30:20 2017 -0500 +++ b/insect_phenology_model.R Thu Nov 09 11:30:27 2017 -0500 |
| [ |
| b'@@ -25,24 +25,22 @@\n args <- parse_args(parser, positional_arguments=TRUE)\r\n opt <- args$options\r\n \r\n-get_temperature_file_path=function(loc, temperature_data)\r\n+convert_csv_to_rdata=function(loc, temperature_data)\r\n {\r\n- expdata <- matrix(rep(0, opt$num_days * 3), nrow=opt$num_days)\r\n+ expdata <- matrix(rep(0, opt$num_days * 6), nrow=opt$num_days)\r\n expdata[,1] <- c(1:opt$num_days)\r\n # Minimum\r\n- expdata[,2] <- temperature_data[c(1:opt$num_days), 3]\r\n+ expdata[,2] <- temperature_data[c(1:opt$num_days), 5]\r\n # Maximum\r\n- expdata[,3] <- temperature_data[c(1:opt$num_days), 2]\r\n- date <- temperature_data[1, 3]\r\n- year <- substr(date, 1, 4)\r\n- namedat <- paste(loc, year, ".Rdat", sep="")\r\n+ expdata[,3] <- temperature_data[c(1:opt$num_days), 6]\r\n+ namedat <- "tempdata.Rdat"\r\n save(expdata, file=namedat)\r\n namedat\r\n }\r\n \r\n daylength=function(latitude, num_days)\r\n {\r\n- # from Forsythe 1995\r\n+ # From Forsythe 1995.\r\n p=0.8333\r\n dl <- NULL\r\n for (i in 1:num_days) {\r\n@@ -50,15 +48,16 @@\n phi <- asin(0.39795 * cos(theta))\r\n dl[i] <- 24 - 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi)))\r\n }\r\n- # return a vector of daylength for the number of\r\n- # days specifie din the input temperature data\r\n+ # Return a vector of daylength for the number of\r\n+ # days specified in the input temperature data.\r\n dl\r\n }\r\n \r\n hourtemp=function(latitude, date, temperature_file_path, num_days)\r\n {\r\n load(temperature_file_path)\r\n- # base development threshold for BMSB\r\n+ # Base development threshold for Brown Marmolated Stink Bug\r\n+ # insect phenology model.\r\n threshold <- 14.17\r\n dnp <- expdata[date, 2] # daily minimum\r\n dxp <- expdata[date, 3] # daily maximum\r\n@@ -69,27 +68,29 @@\n dd <- 0\r\n }\r\n else {\r\n- # extract daylength data for the number of\r\n- # days specified in the input temperature data\r\n+ # Extract daylength data for the number of\r\n+ # days specified in the input temperature data.\r\n dlprofile <- daylength(latitude, num_days)\r\n- T <- NULL # initialize hourly temperature\r\n- dh <- NULL #initialize degree hour vector\r\n- # calculate daylength in given date\r\n+ # Initialize hourly temperature.\r\n+ T <- NULL\r\n+ # Initialize degree hour vector.\r\n+ dh <- NULL\r\n+ # Calculate daylength in given date.\r\n y <- dlprofile[date]\r\n- # night length\r\n+ # Night length.\r\n z <- 24 - y\r\n- # lag coefficient\r\n+ # Lag coefficient.\r\n a <- 1.86\r\n- # night coefficient\r\n+ # Night coefficient.\r\n b <- 2.20\r\n- # sunrise time\r\n+ # Sunrise time.\r\n risetime <- 12 - y / 2\r\n- # sunset time\r\n+ # Sunset time.\r\n settime <- 12 + y / 2\r\n ts <- (dxp - dnp) * sin(pi * (settime - 5) / (y + 2 * a)) + dnp\r\n for (i in 1:24) {\r\n if (i > risetime && i<settime) {\r\n- # number of hours after Tmin until sunset\r\n+ # Number of hours after Tmin until sunset.\r\n m <- i - 5\r\n T[i]=(dxp - dnp) * sin(pi * m / (y + 2 * a)) + dnp\r\n if (T[i]<8.4) {\r\n@@ -197,44 +198,45 @@\n return\r\n }\r\n \r\n-cat("Replications: ", opt$replications, "\\n")\r\n-cat("Photoperiod: ", opt$photoperiod, "\\n")\r\n-cat("Oviposition rate: ", opt$oviposition, "\\n")\r\n-cat("Egg mortality rate: ", opt$egg_mort, "\\n")\r\n-cat("Nymph mortality rate: ", opt$nymph_mort, "\\n")\r\n-cat("Adult mortality rate: ", opt$adult_mort, "\\n")\r\n-cat("Min clutch size: ", opt$min_clutch_size, "\\n")\r\n-cat("Max clutch size: ", opt$max_clutch_size, "\\n")\r\n-cat("(egg->young nymph): ", opt$young_nymph_accum, "\\n")\r\n-cat("(young nymph->old nymph): ", opt$old_nymph_accum, "\\n")\r\n-cat("(old nymph->adult): ", opt$adult_accum, "\\n")\r\n-\r\n # Read in the input temperat'..b'elogenic or diapausing?\r\n if (vec.ind[2] == 2) {\r\n- # old nymph stage\r\n- # add to DD\r\n+ # Old nymph stage.\r\n+ # add to dd.\r\n vec.ind[3] <- vec.ind[3] + dd.temp\r\n if (vec.ind[3] >= (200 + opt$adult_accum)) {\r\n- # from old to adult, DD requirement met\r\n+ # From old to adult, dd requirement met.\r\n current.gen <- vec.ind[1]\r\n if (vec.ind[5] == 0) {\r\n- # non-diapausing adult -- previttelogenic\r\n+ # Non-diapausing adult -- previttelogenic.\r\n vec.ind <- c(current.gen, 3, 0, 0, 0)\r\n }\r\n else {\r\n- # diapausing \r\n+ # Diapausing.\r\n vec.ind <- c(current.gen, 5, 0, 0, 1)\r\n }\r\n }\r\n else {\r\n- # add 1 day in current stage\r\n+ # Add 1 day in current stage.\r\n vec.ind[4] <- vec.ind[4] + 1\r\n }\r\n vec.mat[i,] <- vec.ind\r\n }\r\n \r\n- # event 4 growing of diapausing adult (unimportant, but still necessary)## \r\n+ # Event 4 growing of diapausing adult (unimportant, but still necessary).\r\n if (vec.ind[2] == 5) {\r\n vec.ind[3] <- vec.ind[3] + dd.temp\r\n vec.ind[4] <- vec.ind[4] + 1\r\n vec.mat[i,] <- vec.ind\r\n }\r\n- } # else if it is still alive\r\n- } # end of the individual bug loop\r\n+ } # Else if it is still alive.\r\n+ } # End of the individual bug loop.\r\n \r\n- # find how many died\r\n+ # Find how many died.\r\n n.death <- length(death.vec)\r\n if (n.death > 0) {\r\n vec.mat <- vec.mat[-death.vec, ]\r\n }\r\n- # remove record of dead\r\n- # find how many new born \r\n+ # Remove record of dead.\r\n+ # Find how many new born.\r\n n.newborn <- length(birth.vec[,1])\r\n vec.mat <- rbind(vec.mat, birth.vec)\r\n- # update population size for the next day\r\n+ # Update population size for the next day.\r\n n <- n - n.death + n.newborn \r\n \r\n- # aggregate results by day\r\n+ # Aggregate results by day.\r\n tot.pop <- c(tot.pop, n) \r\n- # egg\r\n+ # Egg.\r\n s0 <- sum(vec.mat[,2] == 0)\r\n- # young nymph\r\n+ # Young nymph.\r\n s1 <- sum(vec.mat[,2] == 1)\r\n- # old nymph\r\n+ # Old nymph.\r\n s2 <- sum(vec.mat[,2] == 2)\r\n- # previtellogenic\r\n+ # Previtellogenic.\r\n s3 <- sum(vec.mat[,2] == 3)\r\n- # vitellogenic\r\n+ # Vitellogenic.\r\n s4 <- sum(vec.mat[,2] == 4)\r\n- # diapausing\r\n+ # Diapausing.\r\n s5 <- sum(vec.mat[,2] == 5)\r\n- # overwintering adult\r\n+ # Overwintering adult.\r\n gen0 <- sum(vec.mat[,1] == 0)\r\n- # first generation\r\n+ # First generation.\r\n gen1 <- sum(vec.mat[,1] == 1)\r\n- # second generation\r\n+ # Second generation.\r\n gen2 <- sum(vec.mat[,1] == 2)\r\n- # sum of all adults\r\n+ # Sum of all adults.\r\n n.adult <- sum(vec.mat[,2] == 3) + sum(vec.mat[,2] == 4) + sum(vec.mat[,2] == 5)\r\n- # gen.0 pop size\r\n+ # Gen eration 0 pop size.\r\n gen0.pop[day] <- gen0\r\n gen1.pop[day] <- gen1\r\n gen2.pop[day] <- gen2\r\n@@ -527,7 +524,7 @@\n } # end of days specified in the input temperature data\r\n \r\n dd.cum <- cumsum(dd.day)\r\n- # collect all the outputs\r\n+ # Collect all the outputs.\r\n S0.rep[,N.rep] <- S0\r\n S1.rep[,N.rep] <- S1\r\n S2.rep[,N.rep] <- S2\r\n' |