source: trunk/misc/tools/event_selection/verify.py @ 8523

Last change on this file since 8523 was 6251, checked in by rwilson, 16 years ago

Added Event Selection tool.
File size: 7.0 KB
Line 
1#!/bin/env python
2
3"""
4Program to verify that the output of Dave Burbidge's FORTRAN programs
5(list_quakes.f90 and get_multimux.f90) are the same as from the python
6replacement programs.
7
8This program will *not* run under Windows!
9"""
10
11import sys
12if sys.platform == 'win32':
13    print 'Sorry, this program cannot run under Windows'
14    sys.exit(10)
15import os
16import os.path
17import glob
18import re
19import EventSelection as es
20
21# Regions we are going to check in
22RegionNames = ['Australia', 'Indian', 'SW_Pacific']
23
24# where the data files live (Unix/Linux)
25BaseDirectory = '/nas/gemd/georisk_models/tsunami/models/PTHA'
26
27# where the FORTRAN programs live and full program pathnames
28FTNDir = '/nas/gemd/georisk_models/tsunami/models/template/event_selection'
29FTN_list_quakes = os.path.join(FTNDir, 'list_quakes_solaris')
30FTN_get_multimux = os.path.join(FTNDir, 'get_multimux_solaris')
31
32# 'T' files directory
33TFilesDir = 'Tfiles'
34
35# 'multimux' directory
36MultimuxDir = 'multimux'
37
38# invall filename
39InvallFilename = 'i_invall'
40
41# file containing FORTRAN list_quakes input data
42FTN_list_quakes_inputfile = 'lq_input'
43
44# min and max wave heights used, along with step value
45MinWaveHeight = 0.3
46MaxWaveHeight = 1.6
47StepWaveHeight = 0.2
48
49# accuracy to decide of two numbers are equal
50Epsilon = 0.0001
51
52# regular expression pattern to split a line on multiple spaces
53SpacesPattern = re.compile(' +')
54
55
56##
57# @brief Compare the list_quake output files, FORTRAN and python.
58# @param gn The guage number.
59# @param min_wh The minimum wave height.
60# @param max_wh The maximum wave height.
61# @note Writes to stdout if there is a difference.
62def list_quake_compare(gn, min_wh, max_wh):
63    # first, check for outut files identical
64    faultxy_OK = (os.system('diff fault.xy.py fault.xy >/dev/null 2>&1') == 0)
65    quake_prob_OK = (os.system('diff quake_prob.txt.py quake_prob.txt >/dev/null 2>&1') == 0)
66    if faultxy_OK and quake_prob_OK:
67        return
68
69    # set error message
70    error = None
71
72    # now check for logical equality
73    if not faultxy_OK:
74        # if no output from FORTRAN list_quakes, python must be 1 line exactly (header)
75        if not os.path.exists('fault.xy'):
76            fd = open('fault.xy.py', 'r')
77            lines = fd.readlines()
78            fd.close()
79            if len(lines) != 1:
80                error = error + ', ' + 'fault.xy'
81        else:
82            # same number of lines?
83            fd = open('fault.xy.py', 'r')
84            py_lines = fd.readlines()
85            fd.close()
86            fd = open('fault.xy', 'r')
87            ftn_lines = fd.readlines()
88            fd.close()
89            if len(py_lines) != len(ftn_lines):
90                error = 'fault.xy'
91            else:
92                # check logical difference - should never be different
93                raise RuntimeError, "Can't get here!?"
94
95    if not quake_prob_OK:
96        # if no output from FORTRAN list_quakes, python must be 1 line exactly (header)
97        if not os.path.exists('quake_prob.txt'):
98            fd = open('quake_prob.txt.py', 'r')
99            lines = fd.readlines()
100            fd.close()
101            if len(lines) != 1:
102                error = error + ', ' + 'quake_prob.txt (1)'
103        else:
104            # same number of lines?
105            fd = open('quake_prob.txt.py', 'r')
106            py_lines = fd.readlines()
107            fd.close()
108            fd = open('quake_prob.txt', 'r')
109            ftn_lines = fd.readlines()
110            fd.close()
111            if len(py_lines) != len(ftn_lines):
112                if error:
113                    error += ', quake_prob.txt (2)'
114                else:
115                    error = 'quake_prob.txt (2)'
116            else:
117                # files are textually different, check for logical difference
118                py_data = read_qp_data(py_lines)
119                ftn_data = read_qp_data(ftn_lines)
120                if not cmp_qp_data(py_data, ftn_data, Epsilon):
121                    if error:
122                        error += ', quake_prob.txt (3)'
123                    else:
124                        error = 'quake_prob.txt (3)'
125
126    if error:
127        print 'Gauge=%d, min=%.1f, max=%.1f, file differ: %s' % (gn, min_wh, max_wh, error)
128
129def read_qp_data(lines):
130    data = [lines[0]]
131    for l in lines[1:]:
132        (_, v1, v2, v3, v4) = re.split(SpacesPattern, l)
133        data.append((int(v1), float(v2), float(v3), float(v4)))
134    return data
135
136
137def cmp_qp_data(d1, d2, eps):
138    # now check each numeric tuple - (%d, %f, %f, %f)
139    # headers should be identical
140    z = zip(d1[1:], d2[1:])
141    for (t1, t2) in z:
142        (d1, f11, f12, f13) = t1
143        (d2, f21, f22, f23) = t2
144        if (d1 != d2) or (f11 != f21) or (f12 != f22) or (f13 != f23):
145            print 'd1=%d, f11=%f, F12=%f, f13=%f' % (d1, f11, f12, f13)
146            print 'd2=%d, f21=%f, F22=%f, f23=%f' % (d2, f21, f22, f23)
147            sys.exit(0)
148    return True
149       
150
151def main():
152    for r in RegionNames:
153        r_dir = os.path.join(BaseDirectory, r)
154        if not os.path.exists(r_dir):
155            print "Sorry, region '%s' doesn't exist" % r
156            continue
157        # link to i_invall file
158        invall_file = os.path.join(r_dir, MultimuxDir, InvallFilename)
159        os.system('rm -f %s' % InvallFilename)
160        os.system('ln -s %s %s' % (invall_file, InvallFilename))
161
162        # get tide gauge numbers, if none, goto next region
163        tf_dir = os.path.join(r_dir, TFilesDir)
164        tf_dir_mask = os.path.join(tf_dir, 'T-[0-9]*')
165        tgn_files = [os.path.basename(x) for x in glob.glob(tf_dir_mask)]
166        if not tgn_files:
167            print "No T files for region '%s'" % r_dir
168            continue
169
170        # now cycle through all T files
171        for t in tgn_files:
172            t_file = os.path.join(tf_dir, t)
173            # link to T file
174            os.system('rm -f %s' % t)
175            os.system('ln -s %s %s' % (t_file, t))
176
177            gn = int(t[2:])
178            min_wh = MinWaveHeight
179            while min_wh <= MaxWaveHeight:
180                max_wh = min_wh
181                while max_wh <= MaxWaveHeight:
182                    # run the python list_quakes
183                    es.list_quakes(gn, min_wh, max_wh, invall_file, t_file, 'fault.xy.py', 'quake_prob.txt.py')
184                    # run the FORTRAN list_quakes
185                    fin = open(FTN_list_quakes_inputfile, 'w')
186                    fin.write('%d\n' % gn)
187                    fin.write('%.1f\n' % min_wh)
188                    fin.write('%.1f\n' % max_wh)
189                    fin.close()
190#                    print '%s < %s >/dev/null' % (FTN_list_quakes, FTN_list_quakes_inputfile)
191                    os.system('rm -f fault.xy quake_prob.txt')
192                    os.system('%s < %s >/dev/null' % (FTN_list_quakes, FTN_list_quakes_inputfile))
193                    # compare list_quakes output
194                    print '%s %d, %.1f, %.1f' % (r, gn, min_wh, max_wh)
195                    list_quake_compare(gn, min_wh, max_wh)
196                    max_wh += StepWaveHeight
197                min_wh += StepWaveHeight
198            os.system('rm -f %s' % t)
199
200if __name__ == '__main__':
201    sys.exit(main())
Note: See TracBrowser for help on using the repository browser.