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
RevLine 
[5687]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
[5691]17from Numeric import zeros, Float, where, greater, less, compress, sqrt, sum
[5687]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    """
[5691]36    Given a list of scenarios that have CSV guage files, calc the
[5687]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
[5689]50        file_sim = outputdir_name + '_' +  quantity + ".csv"
51        file_exp = id + '_exp_' + quantity + '.csv'
52        file_err = outputdir_name + "_" + quantity + "_err.csv"
[5687]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)
[5691]66        time_exp_cut = compress(condition, time_exp)
[5687]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))
[5691]89            rmsd_list.append(norm/sqrt(len(quantity_sim_interp))) 
[5687]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
[5689]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    """
[5687]108   
[5689]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)
[5691]120    quantity_locations = title.split(',')
[5689]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):
[5691]127        fields = line.split(',')
[5689]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   
[5687]134       
135       
136
137# Return a bunch of lists
138# The err files, for all scenarios
[5689]139def err_files(scenarios, outputdir_tag, quantity='stage'):
[5687]140    """
[5691]141    Create a list of err files, for a list of scenarios.
[5687]142    """
143    file_errs = []
144    for scenario in scenarios:
145        id = scenario['scenario_id']
146        outputdir_name = id + outputdir_tag
[5689]147        file_err =  outputdir_name + "_" + quantity + "_err.csv"
[5687]148        file_errs.append(file_err)
149    return file_errs
150   
151
[5689]152def compare_different_settings(outputdir_tag, scenarios, quantity='stage'):
[5691]153    """
154    Calculate the RMSD for all the tests in a scenario
155    """
[5689]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):
[5687]160       
[5689]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)
[5687]172   
173   
174   
[5689]175def err_addition(err_list, number_of_samples_list):
[5687]176    """
[5691]177    This function 'sums' a list of errs and sums a list of samples
[5687]178   
[5691]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   
[5687]182    If this function gets used alot, maybe pull this out and make it an object
183    """
184    err = norm(ensure_numeric(err_list))
[5689]185    number_of_samples = sum(ensure_numeric(number_of_samples_list))
[5687]186
[5689]187    return err, number_of_samples
[5687]188
189                 
190#-------------------------------------------------------------
191if __name__ == "__main__":
192    """
193    """
194    from scenarios import scenarios
195   
[5691]196    #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!!
[5689]197
[5691]198    outputdir_tag = "_nolmts_wdth_0.1_z_0.0_ys_0.01_mta_0.01"
[5687]199    calc_norms = True
200    #calc_norms = False
201    if calc_norms:
[5689]202        auto_rrms(outputdir_tag, scenarios, "stage", y_location_tag=':0.0')
203    compare_different_settings(outputdir_tag, scenarios, "stage")
[5687]204   
Note: See TracBrowser for help on using the repository browser.