comparison utils.R @ 0:1868b7913590 draft

Uploaded
author greg
date Tue, 07 Aug 2018 13:05:10 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:1868b7913590
1 #!/usr/bin/env Rscript
2
3 get_file_path = function(life_stage, base_name, sub_life_stage=NULL) {
4 if (is.null(sub_life_stage)) {
5 lsi = get_life_stage_index(life_stage);
6 file_name = paste(lsi, base_name, sep="_");
7 } else {
8 lsi = get_life_stage_index(life_stage, sub_life_stage=sub_life_stage);
9 file_name = paste(lsi, tolower(sub_life_stage), base_name, sep="_");
10 }
11 file_path = paste("output_plots_dir", file_name, sep="/");
12 return(file_path);
13 }
14
15 get_year_from_date = function(date_str) {
16 date_str_items = strsplit(date_str, "-")[[1]];
17 return (date_str_items[1]);
18 }
19
20 get_life_stage_index = function(life_stage, sub_life_stage=NULL) {
21 # Name collection elements so that they
22 # are displayed in logical order.
23 if (life_stage=="Egg") {
24 lsi = "01";
25 } else if (life_stage=="Nymph") {
26 if (sub_life_stage=="Young") {
27 lsi = "02";
28 } else if (sub_life_stage=="Old") {
29 lsi = "03";
30 } else if (sub_life_stage=="Total") {
31 lsi="04";
32 }
33 } else if (life_stage=="Adult") {
34 if (sub_life_stage=="Pre-vittelogenic") {
35 lsi = "05";
36 } else if (sub_life_stage=="Vittelogenic") {
37 lsi = "06";
38 } else if (sub_life_stage=="Diapausing") {
39 lsi = "07";
40 } else if (sub_life_stage=="Total") {
41 lsi = "08";
42 }
43 } else if (life_stage=="Total") {
44 lsi = "09";
45 }
46 return(lsi);
47 }
48
49 get_mean_and_std_error = function(p_replications, f1_replications, f2_replications) {
50 # P mean.
51 p_m = apply(p_replications, 1, mean);
52 # P standard error.
53 p_se = apply(p_replications, 1, sd) / sqrt(opt$replications);
54 # F1 mean.
55 f1_m = apply(f1_replications, 1, mean);
56 # F1 standard error.
57 f1_se = apply(f1_replications, 1, sd) / sqrt(opt$replications);
58 # F2 mean.
59 f2_m = apply(f2_replications, 1, mean);
60 # F2 standard error.
61 f2_se = apply(f2_replications, 1, sd) / sqrt(opt$replications);
62 return(list(p_m, p_se, f1_m, f1_se, f2_m, f2_se))
63 }
64
65 get_tick_index = function(index, last_tick, ticks, tick_labels, tick_sep) {
66 # The R code tries hard not to draw overlapping tick labels, and so
67 # will omit labels where they would abut or overlap previously drawn
68 # labels. This can result in, for example, every other tick being
69 # labelled. We'll keep track of the last tick to make sure all of
70 # the month labels are displayed, and missing ticks are restricted
71 # to Sundays which have no labels anyway.
72 if (last_tick==0) {
73 return(length(ticks)+1);
74 }
75 last_saved_tick = ticks[[length(ticks)]];
76 if (index-last_saved_tick<tick_sep) {
77 last_saved_month = tick_labels[[length(tick_labels)]];
78 if (last_saved_month=="") {
79 # We're safe overwriting a tick
80 # with no label (i.e., a Sunday tick).
81 return(length(ticks));
82 } else {
83 # Don't eliminate a Month label.
84 return(NULL);
85 }
86 }
87 return(length(ticks)+1);
88 }
89
90 get_total_days = function(is_leap_year) {
91 # Get the total number of days in the current year.
92 if (is_leap_year) {
93 return(366);
94 } else {
95 return(365);
96 }
97 }
98
99 get_x_axis_ticks_and_labels = function(temperature_data_frame, prepend_end_doy_norm=0, append_start_doy_norm=0, date_interval=FALSE) {
100 # Generate a list of ticks and labels for plotting the x axis.
101 if (prepend_end_doy_norm > 0) {
102 prepend_end_norm_row = which(temperature_data_frame$DOY==prepend_end_doy_norm);
103 } else {
104 prepend_end_norm_row = 0;
105 }
106 if (append_start_doy_norm > 0) {
107 append_start_norm_row = which(temperature_data_frame$DOY==append_start_doy_norm);
108 } else {
109 append_start_norm_row = 0;
110 }
111 num_rows = dim(temperature_data_frame)[1];
112 tick_labels = list();
113 ticks = list();
114 current_month_label = NULL;
115 last_tick = 0;
116 if (date_interval) {
117 tick_sep = 0;
118 } else {
119 tick_sep = 3;
120 }
121 for (i in 1:num_rows) {
122 # Get the year and month from the date which
123 # has the format YYYY-MM-DD.
124 date = format(temperature_data_frame$DATE[i]);
125 # Get the month label.
126 items = strsplit(date, "-")[[1]];
127 month = items[2];
128 month_label = month.abb[as.integer(month)];
129 day = as.integer(items[3]);
130 doy = as.integer(temperature_data_frame$DOY[i]);
131 # We're plotting the entire year, so ticks will
132 # occur on Sundays and the first of each month.
133 if (i == prepend_end_norm_row) {
134 # Add a tick for the end of the 30 year normnals data
135 # that was prepended to the year-to-date data.
136 label_str = "End prepended 30 year normals";
137 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep)
138 ticks[tick_index] = i;
139 if (date_interval) {
140 # Append the day to label_str
141 tick_labels[tick_index] = paste(label_str, day, sep=" ");
142 } else {
143 tick_labels[tick_index] = label_str;
144 }
145 last_tick = i;
146 } else if (doy == append_start_doy_norm) {
147 # Add a tick for the start of the 30 year normnals data
148 # that was appended to the year-to-date data.
149 label_str = "Start appended 30 year normals";
150 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep)
151 ticks[tick_index] = i;
152 if (!identical(current_month_label, month_label)) {
153 # Append the month to label_str.
154 label_str = paste(label_str, month_label, spe=" ");
155 current_month_label = month_label;
156 }
157 if (date_interval) {
158 # Append the day to label_str
159 label_str = paste(label_str, day, sep=" ");
160 }
161 tick_labels[tick_index] = label_str;
162 last_tick = i;
163 } else if (i==num_rows) {
164 # Add a tick for the last day of the year.
165 label_str = "";
166 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep)
167 ticks[tick_index] = i;
168 if (!identical(current_month_label, month_label)) {
169 # Append the month to label_str.
170 label_str = month_label;
171 current_month_label = month_label;
172 }
173 if (date_interval) {
174 # Append the day to label_str
175 label_str = paste(label_str, day, sep=" ");
176 }
177 tick_labels[tick_index] = label_str;
178 } else {
179 if (!identical(current_month_label, month_label)) {
180 # Add a tick for the month.
181 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep)
182 ticks[tick_index] = i;
183 if (date_interval) {
184 # Append the day to the month.
185 tick_labels[tick_index] = paste(month_label, day, sep=" ");
186 } else {
187 tick_labels[tick_index] = month_label;
188 }
189 current_month_label = month_label;
190 last_tick = i;
191 }
192 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep)
193 if (!is.null(tick_index)) {
194 if (date_interval) {
195 # Add a tick for every day. The first tick is the
196 # month label, so add a tick only if i is not 1
197 if (i>1 & day>1) {
198 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep)
199 ticks[tick_index] = i;
200 # Add the day as the label.
201 tick_labels[tick_index] = day;
202 last_tick = i;
203 }
204 } else {
205 # Get the day.
206 day = weekdays(as.Date(date));
207 if (day=="Sunday") {
208 # Add a tick if we're on a Sunday.
209 ticks[tick_index] = i;
210 # Add a blank month label so it is not displayed.
211 tick_labels[tick_index] = "";
212 last_tick = i;
213 }
214 }
215 }
216 }
217 }
218 return(list(ticks, tick_labels));
219 }
220
221 render_chart = function(ticks, date_labels, chart_type, plot_std_error, insect, location, latitude, start_date, end_date, days, maxval,
222 replications, life_stage, group, group_std_error, group2=NULL, group2_std_error=NULL, group3=NULL, group3_std_error=NULL,
223 sub_life_stage=NULL) {
224 if (chart_type=="pop_size_by_life_stage") {
225 if (life_stage=="Total") {
226 title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
227 legend_text = c("Egg", "Nymph", "Adult");
228 columns = c(4, 2, 1);
229 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=FALSE, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3);
230 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3);
231 lines(days, group2, lwd=2, lty=1, col=2);
232 lines(days, group3, lwd=2, lty=1, col=4);
233 axis(side=1, at=ticks, labels=date_labels, las=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
234 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
235 if (plot_std_error=="yes") {
236 # Standard error for group.
237 lines(days, group+group_std_error, lty=2);
238 lines(days, group-group_std_error, lty=2);
239 # Standard error for group2.
240 lines(days, group2+group2_std_error, col=2, lty=2);
241 lines(days, group2-group2_std_error, col=2, lty=2);
242 # Standard error for group3.
243 lines(days, group3+group3_std_error, col=4, lty=2);
244 lines(days, group3-group3_std_error, col=4, lty=2);
245 }
246 } else {
247 if (life_stage=="Egg") {
248 title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
249 legend_text = c(life_stage);
250 columns = c(4);
251 } else if (life_stage=="Nymph") {
252 stage = paste(sub_life_stage, "Nymph Pop :", sep=" ");
253 title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
254 legend_text = c(paste(sub_life_stage, life_stage, sep=" "));
255 columns = c(2);
256 } else if (life_stage=="Adult") {
257 stage = paste(sub_life_stage, "Adult Pop", sep=" ");
258 title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
259 legend_text = c(paste(sub_life_stage, life_stage, sep=" "));
260 columns = c(1);
261 }
262 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=FALSE, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3);
263 legend("topleft", legend_text, lty=c(1), col="black", cex=3);
264 axis(side=1, at=ticks, labels=date_labels, las=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
265 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
266 if (plot_std_error=="yes") {
267 # Standard error for group.
268 lines(days, group+group_std_error, lty=2);
269 lines(days, group-group_std_error, lty=2);
270 }
271 }
272 } else if (chart_type=="pop_size_by_generation") {
273 if (life_stage=="Total") {
274 title_str = ": Total Pop by Gen :";
275 } else if (life_stage=="Egg") {
276 title_str = ": Egg Pop by Gen :";
277 } else if (life_stage=="Nymph") {
278 title_str = paste(":", sub_life_stage, "Nymph Pop by Gen", ":", sep=" ");
279 } else if (life_stage=="Adult") {
280 title_str = paste(":", sub_life_stage, "Adult Pop by Gen", ":", sep=" ");
281 }
282 title = paste(insect, ": Reps", replications, title_str, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
283 legend_text = c("P", "F1", "F2");
284 columns = c(1, 2, 4);
285 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=FALSE, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3);
286 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3);
287 lines(days, group2, lwd=2, lty=1, col=2);
288 lines(days, group3, lwd=2, lty=1, col=4);
289 axis(side=1, at=ticks, labels=date_labels, las=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
290 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
291 if (plot_std_error=="yes") {
292 # Standard error for group.
293 lines(days, group+group_std_error, lty=2);
294 lines(days, group-group_std_error, lty=2);
295 # Standard error for group2.
296 lines(days, group2+group2_std_error, col=2, lty=2);
297 lines(days, group2-group2_std_error, col=2, lty=2);
298 # Standard error for group3.
299 lines(days, group3+group3_std_error, col=4, lty=2);
300 lines(days, group3-group3_std_error, col=4, lty=2);
301 }
302 }
303 }
304
305 stop_err = function(msg) {
306 cat(msg, file=stderr());
307 quit(save="no", status=1);
308 }
309
310 validate_date = function(date_str) {
311 valid_date = as.Date(date_str, format="%Y-%m-%d");
312 if( class(valid_date)=="try-error" || is.na(valid_date)) {
313 msg = paste("Invalid date: ", date_str, ", valid date format is yyyy-mm-dd.", sep="");
314 stop_err(msg);
315 }
316 return(valid_date);
317 }