source: anuga_validation/Hinwood_2008/calc_norm.py @ 5691

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

Adding a header to the boundary files

File size: 6.9 KB
Line 
1"""
2
3All err results are going into the same dir, and it can't really be changed.
4"""
5
6
7#----------------------------------------------------------------------------
8# Import necessary modules
9#----------------------------------------------------------------------------
10
11# Standard modules
12import os
13from csv import writer
14from time import localtime, strftime
15
16# Related major packages
17from Numeric import zeros, Float, where, greater, less, compress, sqrt, sum
18
19from anuga.shallow_water.data_manager import csv2dict
20from anuga.utilities.numerical_tools import ensure_numeric, err, norm
21from anuga.utilities.interp import interp
22
23def get_max_min_condition_array(min, max, vector):
24   
25    SMALL_MIN = -1e10  # Not that small, but small enough
26    vector = ensure_numeric(vector)
27    assert min > SMALL_MIN
28    no_maxs = where(less(vector,max), vector, SMALL_MIN)
29    band_condition = greater(no_maxs, min)
30    return band_condition
31
32   
33def auto_rrms(outputdir_tag, scenarios, quantity='stage',
34              y_location_tag=':0.0'):
35    """
36    Given a list of scenarios that have CSV guage files, calc the
37    err, Number_of_samples and rmsd for all gauges in each scenario.
38    Write this info to a file for each scenario.
39    """
40    for run_data in scenarios:           
41        location_sims = []
42        location_exps = []
43        for gauge_x in run_data['gauge_x']:
44            gauge_x = str(gauge_x)
45            location_sims.append(gauge_x + y_location_tag)
46            location_exps.append(gauge_x)
47       
48        id = run_data['scenario_id']
49        outputdir_name = id + outputdir_tag
50        file_sim = outputdir_name + '_' +  quantity + ".csv"
51        file_exp = id + '_exp_' + quantity + '.csv'
52        file_err = outputdir_name + "_" + quantity + "_err.csv"
53       
54
55        simulation, _ = csv2dict(file_sim)
56        experiment, _ = csv2dict(file_exp)
57       
58        time_sim = [float(x) for x in simulation['time']]
59        time_exp = [float(x) for x in experiment['Time']]
60        time_sim = ensure_numeric(time_sim)
61        time_exp = ensure_numeric(time_exp)
62       
63        condition = get_max_min_condition_array(run_data['wave_times'][0],
64                                                run_data['wave_times'][1],
65                                                time_exp)
66        time_exp_cut = compress(condition, time_exp)
67       
68        print "Writing to ", file_err
69       
70        err_list = []
71        points = []
72        rmsd_list = []
73        for location_sim, location_exp in map(None, location_sims,
74                                              location_exps):
75            quantity_sim = [float(x) for x in simulation[location_sim]]
76            quantity_exp = [float(x) for x in experiment[location_exp]]
77
78            quantity_exp_cut = compress(condition, quantity_exp)
79
80            # Now let's do interpolation
81            quantity_sim_interp = interp(quantity_sim, time_sim, time_exp_cut)
82
83            assert len(quantity_sim_interp) == len(quantity_exp_cut)
84            norm = err(quantity_sim_interp,
85                       quantity_exp_cut,
86                       2, relative = False)  # 2nd norm (rel. RMS)
87            err_list.append(norm)
88            points.append(len(quantity_sim_interp))
89            rmsd_list.append(norm/sqrt(len(quantity_sim_interp))) 
90        assert len(location_exps) == len(err_list)
91
92        # Writing the file out for one scenario
93        a_writer = writer(file(file_err, "wb"))
94        a_writer.writerow(["x location", "err", "Number_of_samples", "rmsd"])
95        a_writer.writerows(map(None,
96                               location_exps,
97                               err_list,
98                               points,
99                               rmsd_list))
100
101
102
103def load_sensors(quantity_file):
104    """
105    Load a csv file, where the first row is the column header and
106    the first colum explains the rows.
107    """
108   
109    # Read the depth file
110    dfid = open(quantity_file)
111    lines = dfid.readlines()
112    dfid.close()
113
114    title = lines.pop(0)
115    n_time = len(lines)
116    n_sensors = len(lines[0].split(','))-1  # -1 to remove time
117    times = zeros(n_time, Float)  #Time
118    depths = zeros(n_time, Float)  #
119    sensors = zeros((n_time,n_sensors), Float)
120    quantity_locations = title.split(',')
121    quantity_locations.pop(0) # remove 'time'
122
123    # Doing j.split(':')[0] drops the y location
124    locations = [float(j.split(':')[0]) for j in quantity_locations]
125   
126    for i, line in enumerate(lines):
127        fields = line.split(',')
128        fields = [float(j) for j in fields]
129        times[i] = fields[0]
130        sensors[i] = fields[1:] # 1: to remove time
131
132    return times, locations, sensors                 
133   
134       
135       
136
137# Return a bunch of lists
138# The err files, for all scenarios
139def err_files(scenarios, outputdir_tag, quantity='stage'):
140    """
141    Create a list of err files, for a list of scenarios.
142    """
143    file_errs = []
144    for scenario in scenarios:
145        id = scenario['scenario_id']
146        outputdir_name = id + outputdir_tag
147        file_err =  outputdir_name + "_" + quantity + "_err.csv"
148        file_errs.append(file_err)
149    return file_errs
150   
151
152def compare_different_settings(outputdir_tag, scenarios, quantity='stage'):
153    """
154    Calculate the RMSD for all the tests in a scenario
155    """
156    files = err_files(scenarios, outputdir_tag, quantity=quantity)
157    err = 0.0
158    number_of_samples = 0
159    for run_data, file in map(None, scenarios, files):
160       
161        simulation, _ = csv2dict(file)
162        err_list = [float(x) for x in simulation['err']]
163        number_of_samples_list = [float(x) for x in \
164                                  simulation['Number_of_samples']]
165       
166        if number_of_samples is not 0:
167            err_list.append(err)
168            number_of_samples_list.append(number_of_samples)
169        err, number_of_samples = err_addition(err_list, number_of_samples_list)
170    rmsd = err/sqrt(number_of_samples)
171    print outputdir_tag + "   " + str(rmsd)
172   
173   
174   
175def err_addition(err_list, number_of_samples_list):
176    """
177    This function 'sums' a list of errs and sums a list of samples
178   
179    err is the err value (sqrt(sum_over_x&y((xi - yi)^2))) for a set of values.
180    number_of_samples is the number of values associated with the err.
181   
182    If this function gets used alot, maybe pull this out and make it an object
183    """
184    err = norm(ensure_numeric(err_list))
185    number_of_samples = sum(ensure_numeric(number_of_samples_list))
186
187    return err, number_of_samples
188
189                 
190#-------------------------------------------------------------
191if __name__ == "__main__":
192    """
193    """
194    from scenarios import scenarios
195   
196    #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
197
198    outputdir_tag = "_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01"
199    calc_norms = True
200    #calc_norms = False
201    if calc_norms:
202        auto_rrms(outputdir_tag, scenarios, "stage", y_location_tag=':0.0')
203    compare_different_settings(outputdir_tag, scenarios, "stage")
204   
Note: See TracBrowser for help on using the repository browser.