| Previous changeset 7:ad26f07a7dd8 (2018-02-13) Next changeset 9:d9371485aaf5 (2018-02-27) |
|
Commit message:
Uploaded |
|
modified:
insect_phenology_model.R |
| b |
| diff -r ad26f07a7dd8 -r 37f1ad91a949 insect_phenology_model.R --- a/insect_phenology_model.R Tue Feb 13 13:53:29 2018 -0500 +++ b/insect_phenology_model.R Tue Feb 13 13:53:37 2018 -0500 |
| [ |
| b'@@ -23,145 +23,165 @@\n make_option(c("--young_nymph_accumulation"), action="store", dest="young_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)")\n )\n \n-parser <- OptionParser(usage="%prog [options] file", option_list=option_list)\n-args <- parse_args(parser, positional_arguments=TRUE)\n-opt <- args$options\n+parser <- OptionParser(usage="%prog [options] file", option_list=option_list);\n+args <- parse_args(parser, positional_arguments=TRUE);\n+opt <- args$options;\n \n add_daylight_length = function(temperature_data_frame, num_columns, num_rows) {\n # Return a vector of daylight length (photoperido profile) for\n # the number of days specified in the input temperature data\n # (from Forsythe 1995).\n- p = 0.8333\n- latitude <- temperature_data_frame$LATITUDE[1]\n- daylight_length_vector <- NULL\n+ p = 0.8333;\n+ latitude = temperature_data_frame$LATITUDE[1];\n+ daylight_length_vector = NULL;\n for (i in 1:num_rows) {\n # Get the day of the year from the current row\n # of the temperature data for computation.\n- doy <- temperature_data_frame$DOY[i]\n- theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (doy - 186)))\n- phi <- asin(0.39795 * cos(theta))\n+ doy = temperature_data_frame$DOY[i];\n+ theta = 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (doy - 186)));\n+ phi = asin(0.39795 * cos(theta));\n # Compute the length of daylight for the day of the year.\n- darkness_length <- 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi)))\n- daylight_length_vector[i] <- 24 - darkness_length\n+ darkness_length = 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi)));\n+ daylight_length_vector[i] = 24 - darkness_length;\n }\n # Append daylight_length_vector as a new column to temperature_data_frame.\n- temperature_data_frame[, num_columns+1] <- daylight_length_vector\n- return(temperature_data_frame)\n+ temperature_data_frame[, num_columns+1] = daylight_length_vector;\n+ return(temperature_data_frame);\n }\n \n dev.egg = function(temperature) {\n- dev.rate = -0.9843 * temperature + 33.438\n- return(dev.rate)\n+ dev.rate = -0.9843 * temperature + 33.438;\n+ return(dev.rate);\n }\n \n dev.emerg = function(temperature) {\n- emerg.rate <- -0.5332 * temperature + 24.147\n- return(emerg.rate)\n+ emerg.rate = -0.5332 * temperature + 24.147;\n+ return(emerg.rate);\n }\n \n dev.old = function(temperature) {\n- n34 <- -0.6119 * temperature + 17.602\n- n45 <- -0.4408 * temperature + 19.036\n- dev.rate = mean(n34 + n45)\n- return(dev.rate)\n+ n34 = -0.6119 * temperature + 17.602;\n+ n45 = -0.4408 * temperature + 19.036;\n+ dev.rate = mean(n34 + n45);\n+ return(dev.rate);\n }\n \n dev.young = function(temperature) {\n- n12 <- -0.3728 * temperature + 14.68\n- n23 <- -0.6119 * temperature + 25.249\n- dev.rate = mean(n12 + n23)\n- return(dev.rate)\n+ n12 = -0.3728 * temperature + 14.68;\n+ n23 = -0.6119 * temperature + 25.249;\n+ dev.rate = mean(n12 + n23);\n+ return(dev.rate);\n+}\n+\n+\n+get_date_labels = function(temperature_data_frame, num_rows) {\n+ # Keep track of the years to see if spanning years.\n+ month_labels = list();\n+ current_month_label = NULL;\n+ for (i in 1:num_rows) {\n+ # Get the year and month from the date which\n+ # has the format YYYY-MM-DD.\n+ date = format(temperature_data_frame$DATE[i]);\n+ items = strsplit(date, "-")[[1]];\n+ month = items[2];\n+ month_label = month.abb[as.integer(month)];\n+ if (!identical(current_month_label, month_label)) {\n+ month_labels[length(month_labels)+1] = month_label;\n+ current_month_label = month_label;\n+ }\n+ }\n+ return(c(unlist(month_labels)));\n }\n \n get_temperature_at_hour = function(lati'..b' 1, sd) / sqrt(opt$replications)\n+P.std_error = apply(P.replications, 1, sd) / sqrt(opt$replications);\n \n # Mean value for P adults.\n-P_adults <- apply(P_adults.replications, 1, mean)\n+P_adults = apply(P_adults.replications, 1, mean);\n # Standard error for P_adult.\n-P_adults.std_error <- apply(P_adults.replications, 1, sd) / sqrt(opt$replications)\n+P_adults.std_error = apply(P_adults.replications, 1, sd) / sqrt(opt$replications);\n \n # Mean value for F1.\n-F1 <- apply(F1.replications, 1, mean)\n+F1 = apply(F1.replications, 1, mean);\n # Standard error for F1.\n-F1.std_error <- apply(F1.replications, 1, sd) / sqrt(opt$replications)\n+F1.std_error = apply(F1.replications, 1, sd) / sqrt(opt$replications);\n \n # Mean value for F1 adults.\n-F1_adults <- apply(F1_adults.replications, 1, mean)\n+F1_adults = apply(F1_adults.replications, 1, mean);\n # Standard error for F1 adult.\n-F1_adults.std_error <- apply(F1_adults.replications, 1, sd) / sqrt(opt$replications)\n+F1_adults.std_error = apply(F1_adults.replications, 1, sd) / sqrt(opt$replications);\n \n # Mean value for F2.\n-F2 <- apply(F2.replications, 1, mean)\n+F2 = apply(F2.replications, 1, mean);\n # Standard error for F2.\n-F2.std_error <- apply(F2.replications, 1, sd) / sqrt(opt$replications)\n+F2.std_error = apply(F2.replications, 1, sd) / sqrt(opt$replications);\n \n # Mean value for F2 adults.\n-F2_adults <- apply(F2_adults.replications, 1, mean)\n+F2_adults = apply(F2_adults.replications, 1, mean);\n # Standard error for F2 adult.\n-F2_adults.std_error <- apply(F2_adults.replications, 1, sd) / sqrt(opt$replications)\n+F2_adults.std_error = apply(F2_adults.replications, 1, sd) / sqrt(opt$replications);\n \n # Display the total number of days in the Galaxy history item blurb.\n-cat("Number of days: ", opt$num_days, "\\n")\n+cat("Number of days: ", opt$num_days, "\\n");\n \n-dev.new(width=20, height=30)\n+dev.new(width=20, height=30);\n \n # Start PDF device driver to save charts to output.\n-pdf(file=opt$output, width=20, height=30, bg="white")\n-par(mar=c(5, 6, 4, 4), mfrow=c(3, 1))\n+pdf(file=opt$output, width=20, height=30, bg="white");\n+par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));\n \n # Data analysis and visualization plots only within a single calendar year.\n-days <- c(1:opt$num_days)\n-start_date <- temperature_data_frame$DATE[1]\n-end_date <- temperature_data_frame$DATE[opt$num_days]\n+days = c(1:opt$num_days);\n+start_date = temperature_data_frame$DATE[1];\n+end_date = temperature_data_frame$DATE[opt$num_days];\n \n # Subfigure 1: population size by life stage.\n-maxval <- max(eggs+eggs.std_error, nymphs+nymphs.std_error, adults+adults.std_error)\n+maxval = max(eggs+eggs.std_error, nymphs+nymphs.std_error, adults+adults.std_error);\n render_chart("pop_size_by_life_stage", opt$insect, opt$location, latitude, start_date, end_date, days, maxval,\n- opt$std_error_plot, adults, nymphs, eggs, adults.std_error, nymphs.std_error, eggs.std_error)\n+ opt$std_error_plot, adults, nymphs, eggs, adults.std_error, nymphs.std_error, eggs.std_error, date_labels);\n # Subfigure 2: population size by generation.\n-maxval <- max(F2)\n+maxval = max(F2);\n render_chart("pop_size_by_generation", opt$insect, opt$location, latitude, start_date, end_date, days, maxval,\n- opt$std_error_plot, P, F1, F2, P.std_error, F1.std_error, F2.std_error)\n+ opt$std_error_plot, P, F1, F2, P.std_error, F1.std_error, F2.std_error, date_labels);\n # Subfigure 3: adult population size by generation.\n-maxval <- max(F2_adults) + 100\n+maxval = max(F2_adults) + 100;\n render_chart("adult_pop_size_by_generation", opt$insect, opt$location, latitude, start_date, end_date, days, maxval,\n- opt$std_error_plot, P_adults, F1_adults, F2_adults, P_adults.std_error, F1_adults.std_error, F2_adults.std_error)\n+ opt$std_error_plot, P_adults, F1_adults, F2_adults, P_adults.std_error, F1_adults.std_error, F2_adults.std_error,\n+ date_labels);\n \n # Turn off device driver to flush output.\n-dev.off()\n+dev.off();\n' |