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

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

Current Hinwood scenario calculating RMSD's

File size: 10.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 Numeric import arange, array, zeros, Float, where, greater, less, \
17     compress, argmin, choose, searchsorted, sqrt, sum
18
19# Related major packages
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
26from interp import interp
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   
38def auto_rrms(outputdir_tag, scenarios, quantity):
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)
56            location_sims.append(gauge_x + ':0.0')
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)
80       
81        condition = get_max_min_condition_array(run_data['wave_times'][0],
82                                                run_data['wave_times'][1],
83                                                time_exp)
84        time_exp_cut = compress(condition, time_exp) #, axis=axis)
85        #print "min(time_exp_cut)", min(time_exp_cut)
86        #print "max(time_exp_cut)", max(time_exp_cut)
87       
88        #assert min(time_sim) < min(time_exp)
89       
90        print "Writing to ", file_err
91       
92        err_list = []
93        points = []
94        rmsd_list = []
95        for location_sim, location_exp in map(None, location_sims,
96                                              location_exps):
97            quantity_sim = [float(x) for x in simulation[location_sim]]
98            quantity_exp = [float(x) for x in experiment[location_exp]]
99
100            quantity_exp_cut = compress(condition, quantity_exp)
101
102            # Now let's do interpolation
103            quantity_sim_interp = interp(quantity_sim, time_sim, time_exp_cut)
104
105            assert len(quantity_sim_interp) == len(quantity_exp_cut)
106            norm = err(quantity_sim_interp,
107                       quantity_exp_cut,
108                       2, relative = False)  # 2nd norm (rel. RMS)
109            err_list.append(norm)
110            points.append(len(quantity_sim_interp))
111            rmsd_list.append(norm/sqrt(len(quantity_sim_interp)))             
112            #print "norm", norm
113            #for i in range(len(quantity_sim_interp)):
114               
115                #print "quantity_sim_interp", quantity_sim_interp[i]
116                #print "quantity_exp_cut", quantity_exp_cut[i]
117        assert len(location_exps) == len(err_list)
118
119        # Writing the file out for one scenario
120        a_writer = writer(file(file_err, "wb"))
121        a_writer.writerow(["x location", "err", "Number_of_samples", "rmsd"])
122        a_writer.writerows(map(None,
123                               location_exps,
124                               err_list,
125                               points,
126                               rmsd_list))
127       
128def plot_rrms_sensor_settings(run_data, outputdir_tags, quantity,
129              save_as=None,
130              is_interactive=False):
131    """
132    For a scenario, do
133    """
134    from pylab import ion, plot, xlabel, ylabel, close, legend, \
135         savefig, title, axis, setp, subplot, grid, axvspan, figlegend
136    from anuga.shallow_water.data_manager import csv2dict
137   
138    plot_type = ".pdf"
139   
140    id = run_data['scenario_id']
141   
142    if is_interactive:
143        ion()
144       
145    subplot(212)   
146    plot_title = id + " Root Mean Square Deviation comparison"
147    title(plot_title)
148    y_label = "RMSD"
149    ylabel(y_label)
150    grid(True)
151
152    lines = []
153    for outputdir_tag in outputdir_tags:
154           
155        outputdir_name = id + outputdir_tag
156        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
157                                       outputdir_name=outputdir_name)
158       
159        file_err = pro_instance.rmsd_dir + sep + outputdir_name + "_" \
160                   + quantity + "_err.csv"
161           
162        simulation, _ = csv2dict(file_err)
163        locations = [float(x) for x in simulation['x location']]
164        rmsd_list = [float(x) for x in simulation['rmsd']]
165        lines.append(plot(locations, rmsd_list))
166       
167    figlegend(lines, outputdir_tags,'upper left')
168   
169    if is_interactive:
170        # Wait for enter pressed
171        raw_input()
172
173    save_as = pro_instance.plots_dir + sep + \
174              id + "_rmsd" + plot_type
175    if save_as is not None:
176        savefig(save_as)
177   
178    #Need to close this plot
179    close()
180
181# Return a bunch of lists
182# The err files, for all scenarios
183def err_files(scenarios, outputdir_tag, rmsd_dir, quantity):
184    """
185    The err files, for a list of scenarios
186    """
187    file_errs = []
188    for scenario in scenarios:
189        id = scenario['scenario_id']
190        outputdir_name = id + outputdir_tag
191        file_err = rmsd_dir + sep + outputdir_name + "_" \
192                   + quantity + "_err.csv"
193        file_errs.append(file_err)
194    return file_errs
195   
196
197def plot_settings(outputdir_tags, scenarios, quantity,
198                   save_as=None,
199                   is_interactive=False):
200
201    # A bit hacky.  Getting a pro_instance to get the rmsd_dir.
202    outputdir_name = scenarios[0]['scenario_id'] + outputdir_tags[0]
203    pro_instance = project.Project(['data','flumes','Hinwood_2008'],
204                                   outputdir_name=outputdir_name)
205   
206     
207    settings = {} # keys are different settings.
208    # For each setting there will be err and amount
209       
210    for outputdir_tag in outputdir_tags:
211        files = err_files(scenarios, outputdir_tag, pro_instance.rmsd_dir,
212                          quantity)
213        sim = {}
214        for run_data, file in map(None, scenarios, files):
215           
216            simulation, _ = csv2dict(file)
217            #locations = [float(x) for x in simulation['x location']]
218            err_list = [float(x) for x in simulation['err']]
219            amount_list = [float(x) for x in simulation['Number_of_samples']]
220           
221            if sim.has_key('err'):
222                err_list.append(sim['err'])
223                amount_list.append(sim['amount'])
224            sim['err'], sim['amount'] = err_addition(err_list,
225                                                     amount_list)
226        sim['rmsd'] = sim['err']/sqrt(sim['amount'])
227        settings[outputdir_tag] = sim
228    print "settings", settings
229
230    aux = [(settings[k]['rmsd'], k) for k in settings.keys()]
231    aux.sort()
232    for val, key in aux:
233        print key + "   " + str(val)
234   
235   
236   
237def err_addition(err_list, amount_list):
238    """
239    err1 is the err value (sqrt(sum_over_x&y((xi - yi)^2))) for a set of values
240    amount1 is the number of values associated with the err.
241   
242    If this function gets used alot, maybe pull this out and make it an object
243    """
244    err = norm(ensure_numeric(err_list))
245    amount = sum(ensure_numeric(amount_list))
246
247    return err, amount
248
249   
250       
251def auto_plot_rrms_sensor_settings(outputdir_tags, scenarios, quantity,
252                   save_as=None,
253                   is_interactive=False):
254       
255    for run_data in scenarios:
256        plot_rrms_sensor_settings(run_data, outputdir_tags, quantity) 
257                 
258#-------------------------------------------------------------
259if __name__ == "__main__":
260    """
261    """
262    from scenarios import scenarios
263   
264    outputdir_tags = []
265    outputdir_tags.append("_good_lmts_wdth_0.1_z_0.012_ys_0.01_mta_0.01_F")
266    outputdir_tags.append("_good_lmts_wdth_0.1_z_0.012_ys_0.01_mta_0.001_F")
267    outputdir_tags.append("_good_lmts_wdth_0.1_z_0.012_ys_0.01_mta_0.0001_F")
268    outputdir_tags.append("_good_lmts_wdth_1.0_z_0.012_ys_0.01_mta_0.001_F")
269    outputdir_tags.append("_good_nolmts_wdth_0.1_z_0.012_ys_0.01_mta_0.01_F")
270    outputdir_tags.append("_good_lmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_F")
271    outputdir_tag = "_good_lmts_wdth_1.0_z_0.012_ys_0.01_mta_0.001_F"
272    #outputdir_tag = "_test_limiterC"
273    #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
274    #scenarios = scenarios[4:] # !!!!!!!!!!!!!!!!!!!!!!
275    for outputdir_tag in outputdir_tags:
276        auto_rrms(outputdir_tag, scenarios, "stage")
277   
278    #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
279    #auto_plot_rrms(outputdir_tags, scenarios, "stage")
280    plot_settings(outputdir_tags, scenarios, "stage")
Note: See TracBrowser for help on using the repository browser.