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

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

Work on Patong validation regression test

File size: 5.3 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
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:'
38    print 'python compare_model_timeseries.py swwfile1 swwfile2'
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()
48
49    tempfile = 'xyz1234tempfile.sww' # Has to end with sww
50
51    os.system('cp %s %s' % (sww_filename, tempfile))
52    f = file_function(tempfile, 
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
68    return num.array(timevector), timeseries
69
70
71   
72def compare_timeseries(timevector, 
73                       timeseries1,
74                       timeseries2,
75                       name='',
76                       legend = ('Time series 1', 'Time series 2')):
77    """Compare and plot two timeseries
78    """
79
80   
81    timeseries1 = num.array(timeseries1)
82    timeseries2 = num.array(timeseries2) 
83    N = timevector.shape[0]
84    assert timeseries1.shape[0] == N   
85    assert timeseries2.shape[0] == N       
86   
87    eps = 1.0e-3
88   
89    # Covariance measure   
90    #res = cov(timeseries1, timeseries2)
91    #msg = 'Covariance between timeseries was too large: %e' % res
92    #assert res < 0.01, msg
93     
94    # Difference measures
95    print timeseries1
96    print timeseries1-timeseries2
97   
98   
99    # Maximum norm
100    nominator = max(abs(timeseries1-timeseries2))   
101    denominator = max(abs(timeseries1))
102    print nominator, denominator, N
103    if denominator > 0:
104        # Relative measure
105        res = nominator/denominator
106    else:
107        # Absolute measure
108        res = nominator/N
109               
110    msg = '%s: Difference between timeseries was too large: %e' % (name, res)
111    #assert res < eps, msg
112   
113   
114    nominator = sum(abs(timeseries1-timeseries2))   
115    denominator = sum(abs(timeseries1))
116    print nominator, denominator, N
117    if denominator > 0:
118        # Relative measure
119        res = nominator/denominator
120    else:
121        # Absolute measure
122        res = nominator/N
123               
124    msg = '%s: Difference between timeseries was too large: %e' % (name, res)
125    #assert res < eps, msg
126
127    # Extrema
128    #res = max(model)
129    #report_difference('Maximum', res, expected_maximum[name], rtol, atol)
130   
131   
132
133    # Generate plots
134    pylab.ion() # No plotting on screen
135    pylab.hold(False)
136   
137    pylab.plot(timevector, timeseries1, 'r-',
138               timevector, timeseries2, 'k-')
139
140    pylab.title('Gauge %s' % name)
141    pylab.xlabel('time(s)')
142    pylab.ylabel('stage (m)')   
143    pylab.legend(legend, shadow=True, loc='upper left')
144    pylab.savefig(name, dpi = 300)       
145
146    # Error vector
147    pylab.ion() # No plotting on screen
148    pylab.hold(False)
149    pylab.plot(timevector, timeseries1-timeseries2, 'b-')
150    pylab.title('Gauge %s (difference)' % name)   
151    pylab.xlabel('time(s)')
152    pylab.ylabel('stage difference (m)') 
153    pylab.savefig(name + '_diff', dpi = 300)                 
154   
155   
156
157def compare_all_timeseries(swwfile1, swwfile2):
158   
159    timevector1, timeseries1 = get_timeseries(swwfile1, gauges)
160    timevector2, timeseries2 = get_timeseries(swwfile2, gauges)   
161
162    msg = 'Time vectors were different in models %s and %s' %\
163          (swwfile1, swwfile2)
164    assert num.allclose(timevector1, timevector2), msg
165    timevector = timevector1
166   
167    # Check that both timeseries exist for all gauges
168    for name in timeseries2:
169        assert name in timeseries1
170   
171    for name in timeseries1:
172        assert name in timeseries2
173           
174    # Compare all timeseries data individually
175    for name in timeseries1:
176        compare_timeseries(timevector, 
177                           timeseries1[name],
178                           timeseries2[name],
179                           name=name)
180       
181       
182       
183           
184   
185    print 'OK'
186   
187   
188
189if __name__ == '__main__':
190
191    if len(sys.argv) != 3:
192        usage()
193        import sys; sys.exit(-1) 
194
195    #res = compare_all_timeseries(sys.argv[1], sys.argv[2])
196    try:
197        res = compare_all_timeseries(sys.argv[1], sys.argv[2])
198    except:
199        sys.exit(-1)
200    else:
201        sys.exit(0) 
Note: See TracBrowser for help on using the repository browser.