source: anuga_validation/Hinwood_2008/calc_rmsd.py @ 5695

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

Get the output directory going

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