source: anuga_work/development/Hinwood_2008/calc_norm.py @ 5697

Last change on this file since 5697 was 5697, checked in by duncan, 16 years ago

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