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

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

bug fixes

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