source: anuga_work/development/Hinwood_2008/calc_rmsd.py @ 5710

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

comments

File size: 16.4 KB
RevLine 
[5577]1"""
[5710]2Functions used to calculate the root mean square deviation.
[5577]3
[5710]4Duncan Gray, GA - 2007
5
[5577]6"""
7
8
9#----------------------------------------------------------------------------
10# Import necessary modules
11#----------------------------------------------------------------------------
12
13# Standard modules
14import os
15from os import sep, path
16from csv import writer
[5590]17from time import localtime, strftime
18
19# Related major packages
[5577]20from Numeric import arange, array, zeros, Float, where, greater, less, \
21     compress, argmin, choose, searchsorted, sqrt, sum
22
23import project         
24from os import sep
25from anuga.shallow_water.data_manager import csv2dict
26from anuga.utilities.numerical_tools import ensure_numeric, err, norm
27
28from slope import load_sensors
[5681]29from anuga.utilities.interp import interp
[5577]30
31def get_max_min_condition_array(min, max, vector):
[5710]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.
[5577]36   
[5710]37    precondition: The vector values are ascending.
38   
39    """
40   
[5577]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   
[5616]49def auto_rrms(outputdir_tag, scenarios, quantity, y_location_tag=':0.0'):
[5577]50    """
[5710]51    Given a list of scenarios that have CSV guage files, calc the
[5577]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)
[5616]67            location_sims.append(gauge_x + y_location_tag)
[5577]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)
[5697]91
92        # Trim the experimental data
93        condition_exp = get_max_min_condition_array(run_data['wave_times'][0],
[5577]94                                                run_data['wave_times'][1],
95                                                time_exp)
[5697]96        time_exp_cut = compress(condition_exp, time_exp) #, axis=axis)
[5577]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)
[5697]101
102        # Trim the simulation data
103
104
[5577]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
[5697]116            quantity_exp_cut = compress(condition_exp, quantity_exp)
[5577]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))
[5659]143
144
145                 
[5577]146def plot_rrms_sensor_settings(run_data, outputdir_tags, quantity,
[5664]147              save_as=True,
[5659]148              is_interactive=False, max_rmsd=None):
[5577]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
[5659]155
[5664]156
[5577]157    plot_type = ".pdf"
158   
159    id = run_data['scenario_id']
160   
161    if is_interactive:
162        ion()
[5590]163    time_date = strftime('plot date: %d/%m/%Y Time: %H:%M:%S',
164                                      localtime()) 
[5577]165    subplot(212)   
[5590]166    plot_title = id + " Root Mean Square Deviation comparison" + '\n' \
167                 + time_date
168                 
[5577]169    title(plot_title)
170    y_label = "RMSD"
171    ylabel(y_label)
[5590]172    xlabel("x location, m")
[5577]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       
[5659]190        if max_rmsd is not None:
191            #print "setting axis"
192            axis(ymin=0, ymax=max_rmsd)
[5590]193
194    for break_x in run_data['break_xs']:
195        axvspan(break_x-0.001,break_x+0.001, facecolor='g')
196       
[5670]197    legend(lines, outputdir_tags,'upper left')
[5577]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
[5681]211def 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):
[5670]215    """
216    Produce the RMSD graphs for the anuga article.
217
218    """
[5664]219   
220    tests = []
[5670]221    file_names = []
[5664]222   
223    for i in range(0,len(scenarios),2):
[5681]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]])
[5664]228
229    for test in tests:
[5670]230        file_name = plot_test_rmsd(test, outputdir_tag, quantity,
[5681]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)
[5670]235        file_names.append(file_name)
236    return file_names
[5664]237       
238def plot_test_rmsd(test, outputdir_tag, quantity,
[5681]239                   to_publish_indexes,
240                   save_as=None, add_run_info=False,
241                   is_interactive=False, max_rmsd=None):
[5664]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"   
[5681]251    id = 'S' + test[0]['scenario_id'][1:2]
[5664]252   
[5681]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           
[5664]261    if is_interactive:
262        ion()
263    time_date = strftime('plot date: %d/%m/%Y Time: %H:%M:%S',
264                                      localtime()) 
265    #subplot(212)   
[5681]266   
[5664]267                 
268    #title(plot_title)
269    y_label = "RMSD"
270    ylabel(y_label)
271    xlabel("x location, m")
272    grid(True)
273
274    lines = []
[5670]275    legend_name = []
276    i = 0
[5664]277    for run in test:
[5670]278        i += 1
[5681]279        legend_name.append("S" + test[0]['scenario_id'][1:2]+"R"+str(i))   
[5664]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']]
[5681]290
[5664]291       
[5681]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       
[5664]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)
[5670]310       
[5681]311    if id == 'S2':
312        # Hack to get the desired width
313        axis(xmax=7.0)
314       
315    setp(lines[1], linestyle='--')   
[5670]316    legend(lines, legend_name,'upper left')
[5664]317   
318    if is_interactive:
319        raw_input()
320
321    save_as = pro_instance.plots_dir + sep + \
[5681]322              id + "-rmsd" + plot_type
[5670]323    print "save_as", save_as
[5664]324    if save_as is not None:
325        savefig(save_as)
326   
327    close()
328
[5670]329    return save_as
330
[5577]331# Return a bunch of lists
332# The err files, for all scenarios
333def err_files(scenarios, outputdir_tag, rmsd_dir, quantity):
334    """
[5710]335    Create a list of err files, for a list of scenarios.
[5577]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
[5594]347def compare_different_settings(outputdir_tags, scenarios, quantity,
[5577]348                   save_as=None,
349                   is_interactive=False):
[5710]350    """
351    Calculate the RMSD for all the tests in a scenario
352    """
[5577]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
[5616]381    #print "settings", settings
[5577]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   
390def 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
[5659]402
403def calc_max_rmsd(scenarios, outputdir_tags, quantity):
[5577]404   
[5659]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
[5577]423       
424def auto_plot_rrms_sensor_settings(outputdir_tags, scenarios, quantity,
425                   save_as=None,
426                   is_interactive=False):
[5681]427    """
428    The general auto, for comparing things.
429
430    """
[5659]431    max_rmsd = calc_max_rmsd(scenarios, outputdir_tags, quantity)
432    print "max_rmsd", max_rmsd
[5577]433    for run_data in scenarios:
[5659]434        plot_rrms_sensor_settings(run_data, outputdir_tags, quantity,
435                                  max_rmsd=max_rmsd) 
[5577]436                 
437#-------------------------------------------------------------
438if __name__ == "__main__":
439    """
440    """
441    from scenarios import scenarios
442   
443    outputdir_tags = []
[5670]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")
[5681]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")
[5616]454   
[5697]455    outputdir_tags = []
[5698]456    outputdir_tags.append("_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_I")
[5577]457    #outputdir_tag = "_test_limiterC"
458    #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
459    #scenarios = scenarios[4:] # !!!!!!!!!!!!!!!!!!!!!!
[5709]460    calc_rmsds = True
461    #calc_rmsds = False
462    if calc_rmsds:
[5664]463        for outputdir_tag in outputdir_tags:
464            auto_rrms(outputdir_tag, scenarios, "stage", y_location_tag=':0.0')
[5577]465   
466    #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
[5664]467    #auto_plot_rrms_sensor_settings(outputdir_tags, scenarios, "stage")
[5681]468    compare_different_settings(outputdir_tags, scenarios, "stage")
[5664]469
[5670]470    outputdir_tag = "_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_I" 
[5681]471    #auto_plot_test_rmsd(scenarios, outputdir_tag)
Note: See TracBrowser for help on using the repository browser.