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

Last change on this file since 5664 was 5664, checked in by duncan, 15 years ago

Current Hinwood - calc norm set up for _G end tag comparison
File size: 14.3 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=True,
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
143    plot_type = ".pdf"
144   
145    id = run_data['scenario_id']
146   
147    if is_interactive:
148        ion()
149    time_date = strftime('plot date: %d/%m/%Y Time: %H:%M:%S',
150                                      localtime()) 
151    subplot(212)   
152    plot_title = id + " Root Mean Square Deviation comparison" + '\n' \
153                 + time_date
154                 
155    title(plot_title)
156    y_label = "RMSD"
157    ylabel(y_label)
158    xlabel("x location, m")
159    grid(True)
160
161    lines = []
162    for outputdir_tag in outputdir_tags:
163           
164        outputdir_name = id + outputdir_tag
165        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
166                                       outputdir_name=outputdir_name)
167       
168        file_err = pro_instance.rmsd_dir + sep + outputdir_name + "_" \
169                   + quantity + "_err.csv"
170           
171        simulation, _ = csv2dict(file_err)
172        locations = [float(x) for x in simulation['x location']]
173        rmsd_list = [float(x) for x in simulation['rmsd']]
174        lines.append(plot(locations, rmsd_list))
175       
176        if max_rmsd is not None:
177            #print "setting axis"
178            axis(ymin=0, ymax=max_rmsd)
179
180    for break_x in run_data['break_xs']:
181        axvspan(break_x-0.001,break_x+0.001, facecolor='g')
182       
183    figlegend(lines, outputdir_tags,'upper left')
184   
185    if is_interactive:
186        # Wait for enter pressed
187        raw_input()
188
189    save_as = pro_instance.plots_dir + sep + \
190              id + "_rmsd" + plot_type
191    if save_as is not None:
192        savefig(save_as)
193   
194    #Need to close this plot
195    close()
196
197def auto_plot_test_rmsd(scenarios, outputdir_tag, quantity="stage",
198              save_as=True,
199              is_interactive=False, max_rmsd=1.0):
200   
201    tests = []
202   
203    for i in range(0,len(scenarios),2):
204        tests.append([scenarios[i], scenarios[i+1]])
205
206    for test in tests:
207        plot_test_rmsd(test, outputdir_tag, quantity,
208              save_as=save_as,
209              is_interactive=is_interactive, max_rmsd=max_rmsd)
210       
211def plot_test_rmsd(test, outputdir_tag, quantity,
212              save_as=None,
213              is_interactive=False, max_rmsd=None):
214    """
215    For a test, plot rmsd vs x location for both runs.
216    Also, add the break measurement.
217    """
218    from pylab import ion, plot, xlabel, ylabel, close, legend, \
219         savefig, title, axis, setp, subplot, grid, axvspan, figlegend
220    from anuga.shallow_water.data_manager import csv2dict
221       
222    plot_type = ".pdf"   
223    id = test[0]['scenario_id'][:2]
224   
225    if is_interactive:
226        ion()
227    time_date = strftime('plot date: %d/%m/%Y Time: %H:%M:%S',
228                                      localtime()) 
229    #subplot(212)   
230    plot_title = id + " Root Mean Square Deviation comparison" + '\n' \
231                 + time_date
232                 
233    #title(plot_title)
234    y_label = "RMSD"
235    ylabel(y_label)
236    xlabel("x location, m")
237    grid(True)
238
239    lines = []
240    for run in test:
241           
242        outputdir_name = run['scenario_id'] + outputdir_tag
243        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
244                                       outputdir_name=outputdir_name)
245       
246        file_err = pro_instance.rmsd_dir + sep + outputdir_name + "_" \
247                   + quantity + "_err.csv"
248           
249        simulation, _ = csv2dict(file_err)
250        locations = [float(x) for x in simulation['x location']]
251        rmsd_list = [float(x) for x in simulation['rmsd']]
252        lines.append(plot(locations, rmsd_list))
253       
254        for break_x in run['break_xs']:
255            axvspan(break_x-0.001,break_x+0.001, facecolor='g')
256       
257    if max_rmsd is not None:
258        #print "setting axis"
259        axis(ymin=0, ymax=max_rmsd)
260       
261    figlegend(lines, outputdir_tags,'upper left')
262   
263    if is_interactive:
264        raw_input()
265
266    save_as = pro_instance.plots_dir + sep + \
267              id + "_rmsd" + plot_type
268    if save_as is not None:
269        savefig(save_as)
270   
271    close()
272
273# Return a bunch of lists
274# The err files, for all scenarios
275def err_files(scenarios, outputdir_tag, rmsd_dir, quantity):
276    """
277    The err files, for a list of scenarios
278    """
279    file_errs = []
280    for scenario in scenarios:
281        id = scenario['scenario_id']
282        outputdir_name = id + outputdir_tag
283        file_err = rmsd_dir + sep + outputdir_name + "_" \
284                   + quantity + "_err.csv"
285        file_errs.append(file_err)
286    return file_errs
287   
288
289def compare_different_settings(outputdir_tags, scenarios, quantity,
290                   save_as=None,
291                   is_interactive=False):
292
293    # A bit hacky.  Getting a pro_instance to get the rmsd_dir.
294    outputdir_name = scenarios[0]['scenario_id'] + outputdir_tags[0]
295    pro_instance = project.Project(['data','flumes','Hinwood_2008'],
296                                   outputdir_name=outputdir_name)
297   
298     
299    settings = {} # keys are different settings.
300    # For each setting there will be err and amount
301       
302    for outputdir_tag in outputdir_tags:
303        files = err_files(scenarios, outputdir_tag, pro_instance.rmsd_dir,
304                          quantity)
305        sim = {}
306        for run_data, file in map(None, scenarios, files):
307           
308            simulation, _ = csv2dict(file)
309            #locations = [float(x) for x in simulation['x location']]
310            err_list = [float(x) for x in simulation['err']]
311            amount_list = [float(x) for x in simulation['Number_of_samples']]
312           
313            if sim.has_key('err'):
314                err_list.append(sim['err'])
315                amount_list.append(sim['amount'])
316            sim['err'], sim['amount'] = err_addition(err_list,
317                                                     amount_list)
318        sim['rmsd'] = sim['err']/sqrt(sim['amount'])
319        settings[outputdir_tag] = sim
320    #print "settings", settings
321
322    aux = [(settings[k]['rmsd'], k) for k in settings.keys()]
323    aux.sort()
324    for val, key in aux:
325        print key + "   " + str(val)
326   
327   
328   
329def err_addition(err_list, amount_list):
330    """
331    err1 is the err value (sqrt(sum_over_x&y((xi - yi)^2))) for a set of values
332    amount1 is the number of values associated with the err.
333   
334    If this function gets used alot, maybe pull this out and make it an object
335    """
336    err = norm(ensure_numeric(err_list))
337    amount = sum(ensure_numeric(amount_list))
338
339    return err, amount
340
341
342def calc_max_rmsd(scenarios, outputdir_tags, quantity):
343   
344    max_rmsd = 0
345   
346    for run_data in scenarios:
347        id = run_data['scenario_id']
348        for outputdir_tag in outputdir_tags:
349           
350            outputdir_name = id + outputdir_tag
351            pro_instance = project.Project(['data','flumes','Hinwood_2008'],
352                                           outputdir_name=outputdir_name)
353           
354            file_err = pro_instance.rmsd_dir + sep + outputdir_name + "_" \
355                       + quantity + "_err.csv"
356           
357            simulation, _ = csv2dict(file_err)
358            rmsd_list = [float(x) for x in simulation['rmsd']]
359            max_rmsd = max(max(rmsd_list), max_rmsd)
360    return max_rmsd
361
362       
363def auto_plot_rrms_sensor_settings(outputdir_tags, scenarios, quantity,
364                   save_as=None,
365                   is_interactive=False):
366    max_rmsd = calc_max_rmsd(scenarios, outputdir_tags, quantity)
367    print "max_rmsd", max_rmsd
368    for run_data in scenarios:
369        plot_rrms_sensor_settings(run_data, outputdir_tags, quantity,
370                                  max_rmsd=max_rmsd) 
371                 
372#-------------------------------------------------------------
373if __name__ == "__main__":
374    """
375    """
376    from scenarios import scenarios
377   
378    outputdir_tags = []
379    #outputdir_tags.append("_nolmts_wdth_0.1_z_0.012_ys_0.01_mta_0.01_H")
380    outputdir_tags.append("_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_G")
381    #outputdir_tag = "_good_tri_area_0.01_D"
382    outputdir_tags.append("_nolmts_wdth_0.1_z_0.012_ys_0.01_mta_0.01_G")
383    outputdir_tags.append("_lmts_wdth_0.1_z_0.012_ys_0.01_mta_0.01_G")
384    outputdir_tags.append("_nolmts_wdth_0.1_z_0.012_ys_0.05_mta_1e-05_G")
385    outputdir_tags.append("_nolmts_wdth_0.1_z_0.012_ys_0.01_mta_0.001_G")
386    outputdir_tags.append("_nolmts_wdth_0.01_z_0.012_ys_0.01_mta_1e-05_G")
387    outputdir_tags = []
388    outputdir_tags.append("_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_G")
389   
390    #outputdir_tag = "_good_lmts_wdth_1.0_z_0.012_ys_0.01_mta_0.001_H"
391    #outputdir_tag = "_test_limiterC"
392    #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
393    #scenarios = scenarios[4:] # !!!!!!!!!!!!!!!!!!!!!!
394    if False:
395        for outputdir_tag in outputdir_tags:
396            auto_rrms(outputdir_tag, scenarios, "stage", y_location_tag=':0.0')
397   
398    #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
399    #auto_plot_rrms_sensor_settings(outputdir_tags, scenarios, "stage")
400    #compare_different_settings(outputdir_tags, scenarios, "stage")
401
402    outputdir_tag = "_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_G" 
403    auto_plot_test_rmsd(scenarios, outputdir_tag)
Note: See TracBrowser for help on using the repository browser.