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

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

Current Hinwood scenario

File size: 10.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       
130def plot_rrms_sensor_settings(run_data, outputdir_tags, quantity,
131              save_as=None,
132              is_interactive=False):
133    """
134    For a scenario, do
135    """
136    from pylab import ion, plot, xlabel, ylabel, close, legend, \
137         savefig, title, axis, setp, subplot, grid, axvspan, figlegend
138    from anuga.shallow_water.data_manager import csv2dict
139   
140    plot_type = ".pdf"
141   
142    id = run_data['scenario_id']
143   
144    if is_interactive:
145        ion()
146    time_date = strftime('plot date: %d/%m/%Y Time: %H:%M:%S',
147                                      localtime()) 
148    subplot(212)   
149    plot_title = id + " Root Mean Square Deviation comparison" + '\n' \
150                 + time_date
151                 
152    title(plot_title)
153    y_label = "RMSD"
154    ylabel(y_label)
155    xlabel("x location, m")
156    grid(True)
157
158    lines = []
159    for outputdir_tag in outputdir_tags:
160           
161        outputdir_name = id + outputdir_tag
162        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
163                                       outputdir_name=outputdir_name)
164       
165        file_err = pro_instance.rmsd_dir + sep + outputdir_name + "_" \
166                   + quantity + "_err.csv"
167           
168        simulation, _ = csv2dict(file_err)
169        locations = [float(x) for x in simulation['x location']]
170        rmsd_list = [float(x) for x in simulation['rmsd']]
171        lines.append(plot(locations, rmsd_list))
172       
173
174    for break_x in run_data['break_xs']:
175        axvspan(break_x-0.001,break_x+0.001, facecolor='g')
176       
177    figlegend(lines, outputdir_tags,'upper left')
178   
179    if is_interactive:
180        # Wait for enter pressed
181        raw_input()
182
183    save_as = pro_instance.plots_dir + sep + \
184              id + "_rmsd" + plot_type
185    if save_as is not None:
186        savefig(save_as)
187   
188    #Need to close this plot
189    close()
190
191# Return a bunch of lists
192# The err files, for all scenarios
193def err_files(scenarios, outputdir_tag, rmsd_dir, quantity):
194    """
195    The err files, for a list of scenarios
196    """
197    file_errs = []
198    for scenario in scenarios:
199        id = scenario['scenario_id']
200        outputdir_name = id + outputdir_tag
201        file_err = rmsd_dir + sep + outputdir_name + "_" \
202                   + quantity + "_err.csv"
203        file_errs.append(file_err)
204    return file_errs
205   
206
207def compare_different_settings(outputdir_tags, scenarios, quantity,
208                   save_as=None,
209                   is_interactive=False):
210
211    # A bit hacky.  Getting a pro_instance to get the rmsd_dir.
212    outputdir_name = scenarios[0]['scenario_id'] + outputdir_tags[0]
213    pro_instance = project.Project(['data','flumes','Hinwood_2008'],
214                                   outputdir_name=outputdir_name)
215   
216     
217    settings = {} # keys are different settings.
218    # For each setting there will be err and amount
219       
220    for outputdir_tag in outputdir_tags:
221        files = err_files(scenarios, outputdir_tag, pro_instance.rmsd_dir,
222                          quantity)
223        sim = {}
224        for run_data, file in map(None, scenarios, files):
225           
226            simulation, _ = csv2dict(file)
227            #locations = [float(x) for x in simulation['x location']]
228            err_list = [float(x) for x in simulation['err']]
229            amount_list = [float(x) for x in simulation['Number_of_samples']]
230           
231            if sim.has_key('err'):
232                err_list.append(sim['err'])
233                amount_list.append(sim['amount'])
234            sim['err'], sim['amount'] = err_addition(err_list,
235                                                     amount_list)
236        sim['rmsd'] = sim['err']/sqrt(sim['amount'])
237        settings[outputdir_tag] = sim
238    #print "settings", settings
239
240    aux = [(settings[k]['rmsd'], k) for k in settings.keys()]
241    aux.sort()
242    for val, key in aux:
243        print key + "   " + str(val)
244   
245   
246   
247def err_addition(err_list, amount_list):
248    """
249    err1 is the err value (sqrt(sum_over_x&y((xi - yi)^2))) for a set of values
250    amount1 is the number of values associated with the err.
251   
252    If this function gets used alot, maybe pull this out and make it an object
253    """
254    err = norm(ensure_numeric(err_list))
255    amount = sum(ensure_numeric(amount_list))
256
257    return err, amount
258
259   
260       
261def auto_plot_rrms_sensor_settings(outputdir_tags, scenarios, quantity,
262                   save_as=None,
263                   is_interactive=False):
264       
265    for run_data in scenarios:
266        plot_rrms_sensor_settings(run_data, outputdir_tags, quantity) 
267                 
268#-------------------------------------------------------------
269if __name__ == "__main__":
270    """
271    """
272    from scenarios import scenarios
273   
274    outputdir_tags = []
275    #outputdir_tags.append("_nolmts_wdth_0.1_z_0.012_ys_0.01_mta_0.01_H")
276    outputdir_tags.append("_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_G")
277    #outputdir_tag = "_good_tri_area_0.01_D"
278    outputdir_tags.append("_nolmts_wdth_0.1_z_0.012_ys_0.01_mta_0.01_G")
279    outputdir_tags.append("_lmts_wdth_0.1_z_0.012_ys_0.01_mta_0.01_G")
280    outputdir_tags.append("_nolmts_wdth_0.1_z_0.012_ys_0.05_mta_1e-05_G")
281    outputdir_tags.append("_nolmts_wdth_0.1_z_0.012_ys_0.01_mta_0.001_G")
282    outputdir_tags.append("_nolmts_wdth_0.01_z_0.012_ys_0.01_mta_1e-05_G")
283    outputdir_tags = [] # _no_velocity_nolmts_wdth_0.1_z_0.012_ys_0.01_mta_0.01_G
284    outputdir_tags.append("_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_G")
285   
286    #outputdir_tag = "_good_lmts_wdth_1.0_z_0.012_ys_0.01_mta_0.001_H"
287    #outputdir_tag = "_test_limiterC"
288    #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
289    #scenarios = scenarios[4:] # !!!!!!!!!!!!!!!!!!!!!!
290    #for outputdir_tag in outputdir_tags:
291     #   auto_rrms(outputdir_tag, scenarios, "stage", y_location_tag=':0.0')
292   
293    #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
294    auto_plot_rrms_sensor_settings(outputdir_tags, scenarios, "stage")
295    #compare_different_settings(outputdir_tags, scenarios, "stage")
Note: See TracBrowser for help on using the repository browser.