source: anuga_validation/automated_validation_tests/patong_beach_validation/compare_model_timeseries.py @ 7651

Last change on this file since 7651 was 7651, checked in by ole, 13 years ago

Got patong auto validation working and tested for TRIAL run

File size: 6.7 KB
Line 
1"""Compare selected timeseries from model outputs
2
3Usage:
4python compare_model_timeseries.py swwfile1 swwfile2
5
6Return 0 if timeseries are 'close enough', otherwise return non zero error code
7
8Typically swwfile1 would be model output from unit test and swwfile2 would
9be a reference model.
10"""
11
12from anuga.abstract_2d_finite_volumes.util import file_function
13from anuga.utilities.numerical_tools import cov, get_machine_precision
14
15
16import pylab
17import numpy as num
18import sys, os
19
20try:
21    import pylab
22    pylab.hold(False)  # Check if this command can be issued
23except:
24    print 'Could not import pylab'
25    plotting = False
26else:
27    # Create plots as png files
28    plotting = True
29
30
31
32# Model specific data
33gauges = {'g10a': [422233.4, 874380.2],
34          'g10b': [422107.9, 873556.8],
35          'g10c': [421966.8, 872890.3],
36          'g10d': [421606.1, 872106.2],
37          'g11a': [422628.2, 873933.2],
38          'g11b': [422716.2, 873420.6],
39          'g11c': [422689.1, 872859.8],
40          'g11d': [422408.7, 871940.3],
41          'north': [422572.4, 873992.6],
42          'south': [422004.4, 872014.4]}
43
44
45#---------------------------------------
46
47def usage():
48    print 'Usage:'
49    print 'python compare_model_timeseries.py swwfile1 swwfile2 epsilon'
50
51
52
53def get_timeseries(sww_filename, gauges):
54    """Generate time series for sww file based on gauges
55    """
56
57    gauge_locations = gauges.values()
58    gauge_names = gauges.keys()
59
60    tempfile = 'xyz1234tempfile.sww' # Has to end with sww
61
62    os.system('cp %s %s' % (sww_filename, tempfile))
63    f = file_function(tempfile, 
64                      quantities='stage',
65                      interpolation_points=gauge_locations,
66                      use_cache=True,
67                      verbose=True)
68
69    timevector = f.get_time()
70
71    timeseries = {}
72    for k, name in enumerate(gauge_names):
73        model = timeseries[name] = []
74       
75        for t in timevector:
76            model.append(f(t, point_id=k)[0])
77
78
79    return num.array(timevector), timeseries
80
81
82   
83def compare_timeseries(timevector, 
84                       timeseries1,
85                       timeseries2,
86                       name='',
87                       legend = ('Time series 1', 'Time series 2'),
88                       eps=1.0e-6):
89    """Compare and plot two timeseries
90    """
91
92   
93    timeseries1 = num.array(timeseries1)
94    timeseries2 = num.array(timeseries2) 
95    N = timevector.shape[0]
96    assert timeseries1.shape[0] == N   
97    assert timeseries2.shape[0] == N       
98
99    print 'Testing gauge "%s"' % name
100    print 'epsilon', eps
101    # Covariance measure   
102    res = cov(timeseries1-timeseries2)
103    print '2 norm diff', res       
104    msg = '2-norm of timeseries difference was too large: %e' % res
105    assert res < eps, msg
106     
107    # Maximum norm
108    nominator = max(abs(timeseries1-timeseries2))   
109    denominator = max(abs(timeseries1))
110    if denominator > 0:
111        # Relative measure
112        res = nominator/denominator
113    else:
114        # Absolute measure
115        res = nominator/N
116
117    print 'Max abs diff', res   
118    msg = '%s: Difference between timeseries was too large: %e' % (name, res)
119    assert res < eps, msg
120   
121   
122    nominator = sum(abs(timeseries1-timeseries2))   
123    denominator = sum(abs(timeseries1))
124    if denominator > 0:
125        # Relative measure
126        res = nominator/denominator
127    else:
128        # Absolute measure
129        res = nominator/N
130
131    print 'Sum abs diff', res   
132    msg = '%s: Difference between timeseries was too large: %e' % (name, res)
133    assert res < eps, msg
134
135    # Extrema
136    max1 = max(timeseries1)
137    max2 = max(timeseries2)
138    res = abs(max1-max2)
139    print 'Max diff', res
140    msg = '%s: Difference between maxima was too large: %e' % (name, res)
141    assert res < eps, msg   
142
143    min1 = min(timeseries1)
144    min2 = min(timeseries2)
145    res = abs(min1-min2)   
146    print 'Min diff', res
147    msg = '%s: Difference between minima was too large: %e' % (name, res)
148    assert res < eps, msg
149
150   
151    # Locations of extrema
152    i1 = num.argmax(timeseries1)
153    i2 = num.argmax(timeseries2)
154   
155    res = abs(timevector[i1]-timevector[i2])
156    print 'Max loc diff', res
157    msg = '%s: Difference between location of maxima was too large: %e' % (name, res)
158    assert res < eps, msg
159
160    i1 = num.argmin(timeseries1)
161    i2 = num.argmin(timeseries2)
162   
163    res = abs(timevector[i1]-timevector[i2])   
164    print 'Min loc diff', res   
165    msg = '%s: Difference between location of minima was too large: %e' % (name, res)
166    assert res < eps, msg
167   
168
169    if plotting:
170        # Generate plots
171        #pylab.ion() # No plotting on screen
172        pylab.hold(False)
173   
174        pylab.plot(timevector, timeseries1, 'r-',
175                   timevector, timeseries2, 'k-')
176
177        pylab.title('Gauge %s' % name)
178        pylab.xlabel('time(s)')
179        pylab.ylabel('stage (m)')   
180        pylab.legend(legend, shadow=True, loc='upper left')
181        pylab.savefig(name, dpi = 300)       
182
183        # Error vector
184        #pylab.ion() # No plotting on screen
185        pylab.hold(False)
186        pylab.plot(timevector, timeseries1-timeseries2, 'b-')
187        pylab.title('Gauge %s (difference)' % name)   
188        pylab.xlabel('time(s)')
189        pylab.ylabel('stage difference (m)') 
190        pylab.savefig(name + '_diff', dpi = 300)                 
191
192
193    print 'Gauge "%s" OK' % name
194    print
195   
196
197def compare_all_timeseries(swwfile1, swwfile2, eps=1.0e-6):
198   
199    timevector1, timeseries1 = get_timeseries(swwfile1, gauges)
200    timevector2, timeseries2 = get_timeseries(swwfile2, gauges)   
201
202    msg = 'Time vectors were different in models %s and %s' %\
203          (swwfile1, swwfile2)
204    assert num.allclose(timevector1, timevector2), msg
205    timevector = timevector1
206   
207    # Check that both timeseries exist for all gauges
208    for name in timeseries2:
209        assert name in timeseries1
210   
211    for name in timeseries1:
212        assert name in timeseries2
213           
214    # Compare all timeseries data individually
215    for name in timeseries1:
216        compare_timeseries(timevector, 
217                           timeseries1[name],
218                           timeseries2[name],
219                           name=name,
220                           eps=eps)
221       
222       
223       
224           
225   
226    print 'All gauges OK'
227   
228   
229
230if __name__ == '__main__':
231
232    if len(sys.argv) != 4:
233        usage()
234        sys.exit(-1) 
235
236    eps = float(sys.argv[3]) # Get tolerance to be used in comparisons
237    res = compare_all_timeseries(sys.argv[1], sys.argv[2], eps=eps)
238   
239    #try:
240    #    res = compare_all_timeseries(sys.argv[1], sys.argv[2], eps=eps)
241    #except:
242    #    print 'Failed'
243    #    sys.exit(-1)
244    #else:
245    #    print 'Good', res
246    #    sys.exit(0)
Note: See TracBrowser for help on using the repository browser.