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

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

Current Hinwood - Removed sensor 8 from last 4 runs

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