1 | """ |
---|
2 | Functions used to calculate the root mean square deviation. |
---|
3 | |
---|
4 | Duncan Gray, GA - 2007 |
---|
5 | |
---|
6 | """ |
---|
7 | |
---|
8 | |
---|
9 | #---------------------------------------------------------------------------- |
---|
10 | # Import necessary modules |
---|
11 | #---------------------------------------------------------------------------- |
---|
12 | |
---|
13 | # Standard modules |
---|
14 | import os |
---|
15 | from os import sep, path |
---|
16 | from csv import writer |
---|
17 | from time import localtime, strftime |
---|
18 | |
---|
19 | # Related major packages |
---|
20 | from Numeric import arange, array, zeros, Float, where, greater, less, \ |
---|
21 | compress, argmin, choose, searchsorted, sqrt, sum |
---|
22 | |
---|
23 | import project |
---|
24 | from os import sep |
---|
25 | from anuga.shallow_water.data_manager import csv2dict |
---|
26 | from anuga.utilities.numerical_tools import ensure_numeric, err, norm |
---|
27 | |
---|
28 | from slope import load_sensors |
---|
29 | from anuga.utilities.interp import interp |
---|
30 | |
---|
31 | def get_max_min_condition_array(min, max, vector): |
---|
32 | """ |
---|
33 | Given a vector of values, and minimum and maximum values, return a |
---|
34 | vector of 0/1's that can be used to cut arrays so only the times |
---|
35 | in the min max range are used. |
---|
36 | |
---|
37 | precondition: The vector values are ascending. |
---|
38 | |
---|
39 | """ |
---|
40 | |
---|
41 | SMALL_MIN = -1e10 # Not that small, but small enough |
---|
42 | vector = ensure_numeric(vector) |
---|
43 | assert min > SMALL_MIN |
---|
44 | no_maxs = where(less(vector,max), vector, SMALL_MIN) |
---|
45 | #print "no_maxs", no_maxs |
---|
46 | band_condition = greater(no_maxs, min) |
---|
47 | return band_condition |
---|
48 | |
---|
49 | def auto_rrms(outputdir_tag, scenarios, quantity, y_location_tag=':0.0'): |
---|
50 | """ |
---|
51 | Given a list of scenarios that have CSV guage files, calc the |
---|
52 | err, Number_of_samples and rmsd for all gauges in each scenario. |
---|
53 | Write this info to a file for each scenario. |
---|
54 | """ |
---|
55 | for run_data in scenarios: |
---|
56 | # Don't need to do this. The simulation .csv files are already |
---|
57 | # anuga_times |
---|
58 | #anuga_start_stop_times = [] |
---|
59 | #for time in run_data['wave_start_stop_times']: |
---|
60 | # anuga_start_stop_times.append( \ |
---|
61 | # time - run_data['ANUGA_start_time']) |
---|
62 | |
---|
63 | location_sims = [] |
---|
64 | location_exps = [] |
---|
65 | for gauge_x in run_data['gauge_x']: |
---|
66 | gauge_x = str(gauge_x) |
---|
67 | location_sims.append(gauge_x + y_location_tag) |
---|
68 | location_exps.append(gauge_x) |
---|
69 | |
---|
70 | id = run_data['scenario_id'] |
---|
71 | outputdir_name = id + outputdir_tag |
---|
72 | pro_instance = project.Project(['data','flumes','Hinwood_2008'], |
---|
73 | outputdir_name=outputdir_name) |
---|
74 | end = id + ".csv" |
---|
75 | file_sim = pro_instance.outputdir + quantity + "_" + end |
---|
76 | file_exp = pro_instance.raw_data_dir + sep + id + 'pressfilt_exp_' \ |
---|
77 | + quantity + '.csv' |
---|
78 | file_err = pro_instance.rmsd_dir + sep + outputdir_name + "_" \ |
---|
79 | + quantity + "_err.csv" |
---|
80 | |
---|
81 | |
---|
82 | simulation, _ = csv2dict(file_sim) |
---|
83 | experiment, _ = csv2dict(file_exp) |
---|
84 | |
---|
85 | time_sim = [float(x) for x in simulation['time']] |
---|
86 | time_exp = [float(x) for x in experiment['Time']] |
---|
87 | time_sim = ensure_numeric(time_sim) |
---|
88 | time_exp = ensure_numeric(time_exp) |
---|
89 | #print "min(time_exp)", min(time_exp) |
---|
90 | #print "max(time_exp)", max(time_exp) |
---|
91 | |
---|
92 | # Trim the experimental data |
---|
93 | condition_exp = get_max_min_condition_array(run_data['wave_times'][0], |
---|
94 | run_data['wave_times'][1], |
---|
95 | time_exp) |
---|
96 | time_exp_cut = compress(condition_exp, time_exp) #, axis=axis) |
---|
97 | #print "min(time_exp_cut)", min(time_exp_cut) |
---|
98 | #print "max(time_exp_cut)", max(time_exp_cut) |
---|
99 | |
---|
100 | #assert min(time_sim) < min(time_exp) |
---|
101 | |
---|
102 | # Trim the simulation data |
---|
103 | |
---|
104 | |
---|
105 | |
---|
106 | print "Writing to ", file_err |
---|
107 | |
---|
108 | err_list = [] |
---|
109 | points = [] |
---|
110 | rmsd_list = [] |
---|
111 | for location_sim, location_exp in map(None, location_sims, |
---|
112 | location_exps): |
---|
113 | quantity_sim = [float(x) for x in simulation[location_sim]] |
---|
114 | quantity_exp = [float(x) for x in experiment[location_exp]] |
---|
115 | |
---|
116 | quantity_exp_cut = compress(condition_exp, quantity_exp) |
---|
117 | |
---|
118 | # Now let's do interpolation |
---|
119 | quantity_sim_interp = interp(quantity_sim, time_sim, time_exp_cut) |
---|
120 | |
---|
121 | assert len(quantity_sim_interp) == len(quantity_exp_cut) |
---|
122 | norm = err(quantity_sim_interp, |
---|
123 | quantity_exp_cut, |
---|
124 | 2, relative = False) # 2nd norm (rel. RMS) |
---|
125 | err_list.append(norm) |
---|
126 | points.append(len(quantity_sim_interp)) |
---|
127 | rmsd_list.append(norm/sqrt(len(quantity_sim_interp))) |
---|
128 | #print "norm", norm |
---|
129 | #for i in range(len(quantity_sim_interp)): |
---|
130 | |
---|
131 | #print "quantity_sim_interp", quantity_sim_interp[i] |
---|
132 | #print "quantity_exp_cut", quantity_exp_cut[i] |
---|
133 | assert len(location_exps) == len(err_list) |
---|
134 | |
---|
135 | # Writing the file out for one scenario |
---|
136 | a_writer = writer(file(file_err, "wb")) |
---|
137 | a_writer.writerow(["x location", "err", "Number_of_samples", "rmsd"]) |
---|
138 | a_writer.writerows(map(None, |
---|
139 | location_exps, |
---|
140 | err_list, |
---|
141 | points, |
---|
142 | rmsd_list)) |
---|
143 | |
---|
144 | |
---|
145 | |
---|
146 | def plot_rrms_sensor_settings(run_data, outputdir_tags, quantity, |
---|
147 | save_as=True, |
---|
148 | is_interactive=False, max_rmsd=None): |
---|
149 | """ |
---|
150 | For a scenario, do |
---|
151 | """ |
---|
152 | from pylab import ion, plot, xlabel, ylabel, close, legend, \ |
---|
153 | savefig, title, axis, setp, subplot, grid, axvspan, figlegend |
---|
154 | from anuga.shallow_water.data_manager import csv2dict |
---|
155 | |
---|
156 | |
---|
157 | plot_type = ".pdf" |
---|
158 | |
---|
159 | id = run_data['scenario_id'] |
---|
160 | |
---|
161 | if is_interactive: |
---|
162 | ion() |
---|
163 | time_date = strftime('plot date: %d/%m/%Y Time: %H:%M:%S', |
---|
164 | localtime()) |
---|
165 | subplot(212) |
---|
166 | plot_title = id + " Root Mean Square Deviation comparison" + '\n' \ |
---|
167 | + time_date |
---|
168 | |
---|
169 | title(plot_title) |
---|
170 | y_label = "RMSD" |
---|
171 | ylabel(y_label) |
---|
172 | xlabel("x location, m") |
---|
173 | grid(True) |
---|
174 | |
---|
175 | lines = [] |
---|
176 | for outputdir_tag in outputdir_tags: |
---|
177 | |
---|
178 | outputdir_name = id + outputdir_tag |
---|
179 | pro_instance = project.Project(['data','flumes','Hinwood_2008'], |
---|
180 | outputdir_name=outputdir_name) |
---|
181 | |
---|
182 | file_err = pro_instance.rmsd_dir + sep + outputdir_name + "_" \ |
---|
183 | + quantity + "_err.csv" |
---|
184 | |
---|
185 | simulation, _ = csv2dict(file_err) |
---|
186 | locations = [float(x) for x in simulation['x location']] |
---|
187 | rmsd_list = [float(x) for x in simulation['rmsd']] |
---|
188 | lines.append(plot(locations, rmsd_list)) |
---|
189 | |
---|
190 | if max_rmsd is not None: |
---|
191 | #print "setting axis" |
---|
192 | axis(ymin=0, ymax=max_rmsd) |
---|
193 | |
---|
194 | for break_x in run_data['break_xs']: |
---|
195 | axvspan(break_x-0.001,break_x+0.001, facecolor='g') |
---|
196 | |
---|
197 | legend(lines, outputdir_tags,'upper left') |
---|
198 | |
---|
199 | if is_interactive: |
---|
200 | # Wait for enter pressed |
---|
201 | raw_input() |
---|
202 | |
---|
203 | save_as = pro_instance.plots_dir + sep + \ |
---|
204 | id + "_rmsd" + plot_type |
---|
205 | if save_as is not None: |
---|
206 | savefig(save_as) |
---|
207 | |
---|
208 | #Need to close this plot |
---|
209 | close() |
---|
210 | |
---|
211 | def auto_plot_test_rmsd(scenarios, outputdir_tag, |
---|
212 | to_publish_indexes, quantity="stage", |
---|
213 | save_as=True,add_run_info=False, |
---|
214 | is_interactive=False, max_rmsd=0.012): |
---|
215 | """ |
---|
216 | Produce the RMSD graphs for the anuga article. |
---|
217 | |
---|
218 | """ |
---|
219 | |
---|
220 | tests = [] |
---|
221 | file_names = [] |
---|
222 | |
---|
223 | for i in range(0,len(scenarios),2): |
---|
224 | if to_publish_indexes.has_key(scenarios[i]['scenario_id']): |
---|
225 | tests.append([scenarios[i], scenarios[i+1]]) |
---|
226 | else: |
---|
227 | tests.append([scenarios[i+1], scenarios[i]]) |
---|
228 | |
---|
229 | for test in tests: |
---|
230 | file_name = plot_test_rmsd(test, outputdir_tag, quantity, |
---|
231 | to_publish_indexes, |
---|
232 | save_as=save_as, add_run_info=add_run_info, |
---|
233 | is_interactive=is_interactive, |
---|
234 | max_rmsd=max_rmsd) |
---|
235 | file_names.append(file_name) |
---|
236 | return file_names |
---|
237 | |
---|
238 | def plot_test_rmsd(test, outputdir_tag, quantity, |
---|
239 | to_publish_indexes, |
---|
240 | save_as=None, add_run_info=False, |
---|
241 | is_interactive=False, max_rmsd=None): |
---|
242 | """ |
---|
243 | For a test, plot rmsd vs x location for both runs. |
---|
244 | Also, add the break measurement. |
---|
245 | """ |
---|
246 | from pylab import ion, plot, xlabel, ylabel, close, legend, \ |
---|
247 | savefig, title, axis, setp, subplot, grid, axvspan, figlegend |
---|
248 | from anuga.shallow_water.data_manager import csv2dict |
---|
249 | |
---|
250 | plot_type = ".pdf" |
---|
251 | id = 'S' + test[0]['scenario_id'][1:2] |
---|
252 | |
---|
253 | if add_run_info is True: |
---|
254 | plot_title = "Title will not be in final draft \n" + \ |
---|
255 | id + outputdir_tag |
---|
256 | |
---|
257 | title(plot_title) |
---|
258 | else: |
---|
259 | plot_title = "" |
---|
260 | |
---|
261 | if is_interactive: |
---|
262 | ion() |
---|
263 | time_date = strftime('plot date: %d/%m/%Y Time: %H:%M:%S', |
---|
264 | localtime()) |
---|
265 | #subplot(212) |
---|
266 | |
---|
267 | |
---|
268 | #title(plot_title) |
---|
269 | y_label = "RMSD" |
---|
270 | ylabel(y_label) |
---|
271 | xlabel("x location, m") |
---|
272 | grid(True) |
---|
273 | |
---|
274 | lines = [] |
---|
275 | legend_name = [] |
---|
276 | i = 0 |
---|
277 | for run in test: |
---|
278 | i += 1 |
---|
279 | legend_name.append("S" + test[0]['scenario_id'][1:2]+"R"+str(i)) |
---|
280 | outputdir_name = run['scenario_id'] + outputdir_tag |
---|
281 | pro_instance = project.Project(['data','flumes','Hinwood_2008'], |
---|
282 | outputdir_name=outputdir_name) |
---|
283 | |
---|
284 | file_err = pro_instance.rmsd_dir + sep + outputdir_name + "_" \ |
---|
285 | + quantity + "_err.csv" |
---|
286 | |
---|
287 | simulation, _ = csv2dict(file_err) |
---|
288 | locations = [float(x) for x in simulation['x location']] |
---|
289 | rmsd_list = [float(x) for x in simulation['rmsd']] |
---|
290 | |
---|
291 | |
---|
292 | lines.append(plot(locations, rmsd_list, 'k-x')) # k is black |
---|
293 | |
---|
294 | for i_gauge, (location, rmsd) in enumerate(map(None, |
---|
295 | locations, |
---|
296 | rmsd_list)): |
---|
297 | if to_publish_indexes.has_key(run['scenario_id']) and \ |
---|
298 | to_publish_indexes[run['scenario_id']].count(i_gauge) == 1: |
---|
299 | plot([location], [rmsd], 'ko', markeredgewidth=2, |
---|
300 | markeredgecolor='k') |
---|
301 | |
---|
302 | |
---|
303 | |
---|
304 | for break_x in run['break_xs']: |
---|
305 | axvspan(break_x-0.001,break_x+0.001, facecolor='g') |
---|
306 | |
---|
307 | if max_rmsd is not None: |
---|
308 | #print "setting axis" |
---|
309 | axis(ymin=0, ymax=max_rmsd) |
---|
310 | |
---|
311 | if id == 'S2': |
---|
312 | # Hack to get the desired width |
---|
313 | axis(xmax=7.0) |
---|
314 | |
---|
315 | setp(lines[1], linestyle='--') |
---|
316 | legend(lines, legend_name,'upper left') |
---|
317 | |
---|
318 | if is_interactive: |
---|
319 | raw_input() |
---|
320 | |
---|
321 | save_as = pro_instance.plots_dir + sep + \ |
---|
322 | id + "-rmsd" + plot_type |
---|
323 | print "save_as", save_as |
---|
324 | if save_as is not None: |
---|
325 | savefig(save_as) |
---|
326 | |
---|
327 | close() |
---|
328 | |
---|
329 | return save_as |
---|
330 | |
---|
331 | # Return a bunch of lists |
---|
332 | # The err files, for all scenarios |
---|
333 | def err_files(scenarios, outputdir_tag, rmsd_dir, quantity): |
---|
334 | """ |
---|
335 | Create a list of err files, for a list of scenarios. |
---|
336 | """ |
---|
337 | file_errs = [] |
---|
338 | for scenario in scenarios: |
---|
339 | id = scenario['scenario_id'] |
---|
340 | outputdir_name = id + outputdir_tag |
---|
341 | file_err = rmsd_dir + sep + outputdir_name + "_" \ |
---|
342 | + quantity + "_err.csv" |
---|
343 | file_errs.append(file_err) |
---|
344 | return file_errs |
---|
345 | |
---|
346 | |
---|
347 | def compare_different_settings(outputdir_tags, scenarios, quantity, |
---|
348 | save_as=None, |
---|
349 | is_interactive=False): |
---|
350 | """ |
---|
351 | Calculate the RMSD for all the tests in a scenario |
---|
352 | """ |
---|
353 | |
---|
354 | # A bit hacky. Getting a pro_instance to get the rmsd_dir. |
---|
355 | outputdir_name = scenarios[0]['scenario_id'] + outputdir_tags[0] |
---|
356 | pro_instance = project.Project(['data','flumes','Hinwood_2008'], |
---|
357 | outputdir_name=outputdir_name) |
---|
358 | |
---|
359 | |
---|
360 | settings = {} # keys are different settings. |
---|
361 | # For each setting there will be err and amount |
---|
362 | |
---|
363 | for outputdir_tag in outputdir_tags: |
---|
364 | files = err_files(scenarios, outputdir_tag, pro_instance.rmsd_dir, |
---|
365 | quantity) |
---|
366 | sim = {} |
---|
367 | for run_data, file in map(None, scenarios, files): |
---|
368 | |
---|
369 | simulation, _ = csv2dict(file) |
---|
370 | #locations = [float(x) for x in simulation['x location']] |
---|
371 | err_list = [float(x) for x in simulation['err']] |
---|
372 | amount_list = [float(x) for x in simulation['Number_of_samples']] |
---|
373 | |
---|
374 | if sim.has_key('err'): |
---|
375 | err_list.append(sim['err']) |
---|
376 | amount_list.append(sim['amount']) |
---|
377 | sim['err'], sim['amount'] = err_addition(err_list, |
---|
378 | amount_list) |
---|
379 | sim['rmsd'] = sim['err']/sqrt(sim['amount']) |
---|
380 | settings[outputdir_tag] = sim |
---|
381 | #print "settings", settings |
---|
382 | |
---|
383 | aux = [(settings[k]['rmsd'], k) for k in settings.keys()] |
---|
384 | aux.sort() |
---|
385 | for val, key in aux: |
---|
386 | print key + " " + str(val) |
---|
387 | |
---|
388 | |
---|
389 | |
---|
390 | def err_addition(err_list, amount_list): |
---|
391 | """ |
---|
392 | err1 is the err value (sqrt(sum_over_x&y((xi - yi)^2))) for a set of values |
---|
393 | amount1 is the number of values associated with the err. |
---|
394 | |
---|
395 | If this function gets used alot, maybe pull this out and make it an object |
---|
396 | """ |
---|
397 | err = norm(ensure_numeric(err_list)) |
---|
398 | amount = sum(ensure_numeric(amount_list)) |
---|
399 | |
---|
400 | return err, amount |
---|
401 | |
---|
402 | |
---|
403 | def calc_max_rmsd(scenarios, outputdir_tags, quantity): |
---|
404 | |
---|
405 | max_rmsd = 0 |
---|
406 | |
---|
407 | for run_data in scenarios: |
---|
408 | id = run_data['scenario_id'] |
---|
409 | for outputdir_tag in outputdir_tags: |
---|
410 | |
---|
411 | outputdir_name = id + outputdir_tag |
---|
412 | pro_instance = project.Project(['data','flumes','Hinwood_2008'], |
---|
413 | outputdir_name=outputdir_name) |
---|
414 | |
---|
415 | file_err = pro_instance.rmsd_dir + sep + outputdir_name + "_" \ |
---|
416 | + quantity + "_err.csv" |
---|
417 | |
---|
418 | simulation, _ = csv2dict(file_err) |
---|
419 | rmsd_list = [float(x) for x in simulation['rmsd']] |
---|
420 | max_rmsd = max(max(rmsd_list), max_rmsd) |
---|
421 | return max_rmsd |
---|
422 | |
---|
423 | |
---|
424 | def auto_plot_rrms_sensor_settings(outputdir_tags, scenarios, quantity, |
---|
425 | save_as=None, |
---|
426 | is_interactive=False): |
---|
427 | """ |
---|
428 | The general auto, for comparing things. |
---|
429 | |
---|
430 | """ |
---|
431 | max_rmsd = calc_max_rmsd(scenarios, outputdir_tags, quantity) |
---|
432 | print "max_rmsd", max_rmsd |
---|
433 | for run_data in scenarios: |
---|
434 | plot_rrms_sensor_settings(run_data, outputdir_tags, quantity, |
---|
435 | max_rmsd=max_rmsd) |
---|
436 | |
---|
437 | #------------------------------------------------------------- |
---|
438 | if __name__ == "__main__": |
---|
439 | """ |
---|
440 | """ |
---|
441 | from scenarios import scenarios |
---|
442 | |
---|
443 | outputdir_tags = [] |
---|
444 | outputdir_tags.append("_lmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_I") |
---|
445 | outputdir_tags.append("_no_velocity_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_I") |
---|
446 | outputdir_tags.append("_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.0001_I") |
---|
447 | outputdir_tags.append("_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.001_I") |
---|
448 | outputdir_tags.append("_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_I") |
---|
449 | outputdir_tags.append("_nolmts_wdth_0.1_z_0.012_ys_0.01_mta_0.01_I") |
---|
450 | outputdir_tags.append("_nolmts_wdth_1.0_z_0.0_ys_0.01_mta_0.01_I") |
---|
451 | outputdir_tags.append("_nolmts_wdth_0.1_z_0.012_ys_0.01_mta_0.0001_I") |
---|
452 | outputdir_tags.append("_lmts_wdth_0.1_z_0.0_ys_0.01_mta_0.0001_I") |
---|
453 | outputdir_tags.append("_no_velocity_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.0001_I") |
---|
454 | |
---|
455 | outputdir_tags = [] |
---|
456 | outputdir_tags.append("_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_I") |
---|
457 | #outputdir_tag = "_test_limiterC" |
---|
458 | #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!! |
---|
459 | #scenarios = scenarios[4:] # !!!!!!!!!!!!!!!!!!!!!! |
---|
460 | calc_rmsds = True |
---|
461 | #calc_rmsds = False |
---|
462 | if calc_rmsds: |
---|
463 | for outputdir_tag in outputdir_tags: |
---|
464 | auto_rrms(outputdir_tag, scenarios, "stage", y_location_tag=':0.0') |
---|
465 | |
---|
466 | #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!! |
---|
467 | #auto_plot_rrms_sensor_settings(outputdir_tags, scenarios, "stage") |
---|
468 | compare_different_settings(outputdir_tags, scenarios, "stage") |
---|
469 | |
---|
470 | outputdir_tag = "_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_I" |
---|
471 | #auto_plot_test_rmsd(scenarios, outputdir_tag) |
---|