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

Last change on this file since 7649 was 7649, checked in by ole, 14 years ago

Added variable tolerance and more measures to Patong auto validation

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