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

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

moving Honwood files around

File size: 15.9 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 anuga.utilities.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    legend(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,
198                        to_publish_indexes, quantity="stage",
199              save_as=True,add_run_info=False,
200              is_interactive=False, max_rmsd=0.012):
201    """
202    Produce the RMSD graphs for the anuga article.
203
204    """
205   
206    tests = []
207    file_names = []
208   
209    for i in range(0,len(scenarios),2):
210        if to_publish_indexes.has_key(scenarios[i]['scenario_id']):           
211            tests.append([scenarios[i], scenarios[i+1]])
212        else:           
213            tests.append([scenarios[i+1], scenarios[i]])
214
215    for test in tests:
216        file_name = plot_test_rmsd(test, outputdir_tag, quantity,
217                                   to_publish_indexes,
218                                   save_as=save_as, add_run_info=add_run_info,
219                                   is_interactive=is_interactive,
220                                   max_rmsd=max_rmsd)
221        file_names.append(file_name)
222    return file_names
223       
224def plot_test_rmsd(test, outputdir_tag, quantity,
225                   to_publish_indexes,
226                   save_as=None, add_run_info=False,
227                   is_interactive=False, max_rmsd=None):
228    """
229    For a test, plot rmsd vs x location for both runs.
230    Also, add the break measurement.
231    """
232    from pylab import ion, plot, xlabel, ylabel, close, legend, \
233         savefig, title, axis, setp, subplot, grid, axvspan, figlegend
234    from anuga.shallow_water.data_manager import csv2dict
235       
236    plot_type = ".pdf"   
237    id = 'S' + test[0]['scenario_id'][1:2]
238   
239    if add_run_info is True:
240        plot_title = "Title will not be in final draft \n" + \
241                     id + outputdir_tag
242       
243        title(plot_title)
244    else:
245        plot_title = ""
246           
247    if is_interactive:
248        ion()
249    time_date = strftime('plot date: %d/%m/%Y Time: %H:%M:%S',
250                                      localtime()) 
251    #subplot(212)   
252   
253                 
254    #title(plot_title)
255    y_label = "RMSD"
256    ylabel(y_label)
257    xlabel("x location, m")
258    grid(True)
259
260    lines = []
261    legend_name = []
262    i = 0
263    for run in test:
264        i += 1
265        legend_name.append("S" + test[0]['scenario_id'][1:2]+"R"+str(i))   
266        outputdir_name = run['scenario_id'] + outputdir_tag
267        pro_instance = project.Project(['data','flumes','Hinwood_2008'],
268                                       outputdir_name=outputdir_name)
269       
270        file_err = pro_instance.rmsd_dir + sep + outputdir_name + "_" \
271                   + quantity + "_err.csv"
272           
273        simulation, _ = csv2dict(file_err)
274        locations = [float(x) for x in simulation['x location']]
275        rmsd_list = [float(x) for x in simulation['rmsd']]
276
277       
278        lines.append(plot(locations, rmsd_list, 'k-x')) # k is black
279
280        for i_gauge, (location, rmsd) in enumerate(map(None,
281                                                     locations,
282                                                     rmsd_list)):
283            if to_publish_indexes.has_key(run['scenario_id']) and \
284                   to_publish_indexes[run['scenario_id']].count(i_gauge) == 1:
285                plot([location], [rmsd], 'ko', markeredgewidth=2,
286                     markeredgecolor='k')
287               
288
289       
290        for break_x in run['break_xs']:
291            axvspan(break_x-0.001,break_x+0.001, facecolor='g')
292       
293    if max_rmsd is not None:
294        #print "setting axis"
295        axis(ymin=0, ymax=max_rmsd)
296       
297    if id == 'S2':
298        # Hack to get the desired width
299        axis(xmax=7.0)
300       
301    setp(lines[1], linestyle='--')   
302    legend(lines, legend_name,'upper left')
303   
304    if is_interactive:
305        raw_input()
306
307    save_as = pro_instance.plots_dir + sep + \
308              id + "-rmsd" + plot_type
309    print "save_as", save_as
310    if save_as is not None:
311        savefig(save_as)
312   
313    close()
314
315    return save_as
316
317# Return a bunch of lists
318# The err files, for all scenarios
319def err_files(scenarios, outputdir_tag, rmsd_dir, quantity):
320    """
321    The err files, for a list of scenarios
322    """
323    file_errs = []
324    for scenario in scenarios:
325        id = scenario['scenario_id']
326        outputdir_name = id + outputdir_tag
327        file_err = rmsd_dir + sep + outputdir_name + "_" \
328                   + quantity + "_err.csv"
329        file_errs.append(file_err)
330    return file_errs
331   
332
333def compare_different_settings(outputdir_tags, scenarios, quantity,
334                   save_as=None,
335                   is_interactive=False):
336
337    # A bit hacky.  Getting a pro_instance to get the rmsd_dir.
338    outputdir_name = scenarios[0]['scenario_id'] + outputdir_tags[0]
339    pro_instance = project.Project(['data','flumes','Hinwood_2008'],
340                                   outputdir_name=outputdir_name)
341   
342     
343    settings = {} # keys are different settings.
344    # For each setting there will be err and amount
345       
346    for outputdir_tag in outputdir_tags:
347        files = err_files(scenarios, outputdir_tag, pro_instance.rmsd_dir,
348                          quantity)
349        sim = {}
350        for run_data, file in map(None, scenarios, files):
351           
352            simulation, _ = csv2dict(file)
353            #locations = [float(x) for x in simulation['x location']]
354            err_list = [float(x) for x in simulation['err']]
355            amount_list = [float(x) for x in simulation['Number_of_samples']]
356           
357            if sim.has_key('err'):
358                err_list.append(sim['err'])
359                amount_list.append(sim['amount'])
360            sim['err'], sim['amount'] = err_addition(err_list,
361                                                     amount_list)
362        sim['rmsd'] = sim['err']/sqrt(sim['amount'])
363        settings[outputdir_tag] = sim
364    #print "settings", settings
365
366    aux = [(settings[k]['rmsd'], k) for k in settings.keys()]
367    aux.sort()
368    for val, key in aux:
369        print key + "   " + str(val)
370   
371   
372   
373def err_addition(err_list, amount_list):
374    """
375    err1 is the err value (sqrt(sum_over_x&y((xi - yi)^2))) for a set of values
376    amount1 is the number of values associated with the err.
377   
378    If this function gets used alot, maybe pull this out and make it an object
379    """
380    err = norm(ensure_numeric(err_list))
381    amount = sum(ensure_numeric(amount_list))
382
383    return err, amount
384
385
386def calc_max_rmsd(scenarios, outputdir_tags, quantity):
387   
388    max_rmsd = 0
389   
390    for run_data in scenarios:
391        id = run_data['scenario_id']
392        for outputdir_tag in outputdir_tags:
393           
394            outputdir_name = id + outputdir_tag
395            pro_instance = project.Project(['data','flumes','Hinwood_2008'],
396                                           outputdir_name=outputdir_name)
397           
398            file_err = pro_instance.rmsd_dir + sep + outputdir_name + "_" \
399                       + quantity + "_err.csv"
400           
401            simulation, _ = csv2dict(file_err)
402            rmsd_list = [float(x) for x in simulation['rmsd']]
403            max_rmsd = max(max(rmsd_list), max_rmsd)
404    return max_rmsd
405
406       
407def auto_plot_rrms_sensor_settings(outputdir_tags, scenarios, quantity,
408                   save_as=None,
409                   is_interactive=False):
410    """
411    The general auto, for comparing things.
412
413    """
414    max_rmsd = calc_max_rmsd(scenarios, outputdir_tags, quantity)
415    print "max_rmsd", max_rmsd
416    for run_data in scenarios:
417        plot_rrms_sensor_settings(run_data, outputdir_tags, quantity,
418                                  max_rmsd=max_rmsd) 
419                 
420#-------------------------------------------------------------
421if __name__ == "__main__":
422    """
423    """
424    from scenarios import scenarios
425   
426    outputdir_tags = []
427    outputdir_tags.append("_lmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_I")
428    outputdir_tags.append("_no_velocity_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_I")
429    outputdir_tags.append("_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.0001_I")
430    outputdir_tags.append("_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.001_I")
431    outputdir_tags.append("_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_I")
432    outputdir_tags.append("_nolmts_wdth_0.1_z_0.012_ys_0.01_mta_0.01_I")
433    outputdir_tags.append("_nolmts_wdth_1.0_z_0.0_ys_0.01_mta_0.01_I")
434    outputdir_tags.append("_nolmts_wdth_0.1_z_0.012_ys_0.01_mta_0.0001_I")
435    outputdir_tags.append("_lmts_wdth_0.1_z_0.0_ys_0.01_mta_0.0001_I")
436    outputdir_tags.append("_no_velocity_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.0001_I")
437   
438    #outputdir_tag = "_test_limiterC"
439    #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
440    #scenarios = scenarios[4:] # !!!!!!!!!!!!!!!!!!!!!!
441    calc_norms = True
442    #calc_norms = False
443    if calc_norms:
444        for outputdir_tag in outputdir_tags:
445            auto_rrms(outputdir_tag, scenarios, "stage", y_location_tag=':0.0')
446   
447    #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
448    #auto_plot_rrms_sensor_settings(outputdir_tags, scenarios, "stage")
449    compare_different_settings(outputdir_tags, scenarios, "stage")
450
451    outputdir_tag = "_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01_I" 
452    #auto_plot_test_rmsd(scenarios, outputdir_tag)
Note: See TracBrowser for help on using the repository browser.