#!/bin/env python """ Program to verify that the output of Dave Burbidge's FORTRAN programs (list_quakes.f90 and get_multimux.f90) are the same as from the python replacement programs. This program will *not* run under Windows! """ import sys if sys.platform == 'win32': print 'Sorry, this program cannot run under Windows' sys.exit(10) import os import os.path import glob import re import EventSelection as es # Regions we are going to check in RegionNames = ['Australia', 'Indian', 'SW_Pacific'] # where the data files live (Unix/Linux) BaseDirectory = '/nas/gemd/georisk_models/tsunami/models/PTHA' # where the FORTRAN programs live and full program pathnames FTNDir = '/nas/gemd/georisk_models/tsunami/models/template/event_selection' FTN_list_quakes = os.path.join(FTNDir, 'list_quakes_solaris') FTN_get_multimux = os.path.join(FTNDir, 'get_multimux_solaris') # 'T' files directory TFilesDir = 'Tfiles' # 'multimux' directory MultimuxDir = 'multimux' # invall filename InvallFilename = 'i_invall' # file containing FORTRAN list_quakes input data FTN_list_quakes_inputfile = 'lq_input' # min and max wave heights used, along with step value MinWaveHeight = 0.3 MaxWaveHeight = 1.6 StepWaveHeight = 0.2 # accuracy to decide of two numbers are equal Epsilon = 0.0001 # regular expression pattern to split a line on multiple spaces SpacesPattern = re.compile(' +') ## # @brief Compare the list_quake output files, FORTRAN and python. # @param gn The guage number. # @param min_wh The minimum wave height. # @param max_wh The maximum wave height. # @note Writes to stdout if there is a difference. def list_quake_compare(gn, min_wh, max_wh): # first, check for outut files identical faultxy_OK = (os.system('diff fault.xy.py fault.xy >/dev/null 2>&1') == 0) quake_prob_OK = (os.system('diff quake_prob.txt.py quake_prob.txt >/dev/null 2>&1') == 0) if faultxy_OK and quake_prob_OK: return # set error message error = None # now check for logical equality if not faultxy_OK: # if no output from FORTRAN list_quakes, python must be 1 line exactly (header) if not os.path.exists('fault.xy'): fd = open('fault.xy.py', 'r') lines = fd.readlines() fd.close() if len(lines) != 1: error = error + ', ' + 'fault.xy' else: # same number of lines? fd = open('fault.xy.py', 'r') py_lines = fd.readlines() fd.close() fd = open('fault.xy', 'r') ftn_lines = fd.readlines() fd.close() if len(py_lines) != len(ftn_lines): error = 'fault.xy' else: # check logical difference - should never be different raise RuntimeError, "Can't get here!?" if not quake_prob_OK: # if no output from FORTRAN list_quakes, python must be 1 line exactly (header) if not os.path.exists('quake_prob.txt'): fd = open('quake_prob.txt.py', 'r') lines = fd.readlines() fd.close() if len(lines) != 1: error = error + ', ' + 'quake_prob.txt (1)' else: # same number of lines? fd = open('quake_prob.txt.py', 'r') py_lines = fd.readlines() fd.close() fd = open('quake_prob.txt', 'r') ftn_lines = fd.readlines() fd.close() if len(py_lines) != len(ftn_lines): if error: error += ', quake_prob.txt (2)' else: error = 'quake_prob.txt (2)' else: # files are textually different, check for logical difference py_data = read_qp_data(py_lines) ftn_data = read_qp_data(ftn_lines) if not cmp_qp_data(py_data, ftn_data, Epsilon): if error: error += ', quake_prob.txt (3)' else: error = 'quake_prob.txt (3)' if error: print 'Gauge=%d, min=%.1f, max=%.1f, file differ: %s' % (gn, min_wh, max_wh, error) def read_qp_data(lines): data = [lines[0]] for l in lines[1:]: (_, v1, v2, v3, v4) = re.split(SpacesPattern, l) data.append((int(v1), float(v2), float(v3), float(v4))) return data def cmp_qp_data(d1, d2, eps): # now check each numeric tuple - (%d, %f, %f, %f) # headers should be identical z = zip(d1[1:], d2[1:]) for (t1, t2) in z: (d1, f11, f12, f13) = t1 (d2, f21, f22, f23) = t2 if (d1 != d2) or (f11 != f21) or (f12 != f22) or (f13 != f23): print 'd1=%d, f11=%f, F12=%f, f13=%f' % (d1, f11, f12, f13) print 'd2=%d, f21=%f, F22=%f, f23=%f' % (d2, f21, f22, f23) sys.exit(0) return True def main(): for r in RegionNames: r_dir = os.path.join(BaseDirectory, r) if not os.path.exists(r_dir): print "Sorry, region '%s' doesn't exist" % r continue # link to i_invall file invall_file = os.path.join(r_dir, MultimuxDir, InvallFilename) os.system('rm -f %s' % InvallFilename) os.system('ln -s %s %s' % (invall_file, InvallFilename)) # get tide gauge numbers, if none, goto next region tf_dir = os.path.join(r_dir, TFilesDir) tf_dir_mask = os.path.join(tf_dir, 'T-[0-9]*') tgn_files = [os.path.basename(x) for x in glob.glob(tf_dir_mask)] if not tgn_files: print "No T files for region '%s'" % r_dir continue # now cycle through all T files for t in tgn_files: t_file = os.path.join(tf_dir, t) # link to T file os.system('rm -f %s' % t) os.system('ln -s %s %s' % (t_file, t)) gn = int(t[2:]) min_wh = MinWaveHeight while min_wh <= MaxWaveHeight: max_wh = min_wh while max_wh <= MaxWaveHeight: # run the python list_quakes es.list_quakes(gn, min_wh, max_wh, invall_file, t_file, 'fault.xy.py', 'quake_prob.txt.py') # run the FORTRAN list_quakes fin = open(FTN_list_quakes_inputfile, 'w') fin.write('%d\n' % gn) fin.write('%.1f\n' % min_wh) fin.write('%.1f\n' % max_wh) fin.close() # print '%s < %s >/dev/null' % (FTN_list_quakes, FTN_list_quakes_inputfile) os.system('rm -f fault.xy quake_prob.txt') os.system('%s < %s >/dev/null' % (FTN_list_quakes, FTN_list_quakes_inputfile)) # compare list_quakes output print '%s %d, %.1f, %.1f' % (r, gn, min_wh, max_wh) list_quake_compare(gn, min_wh, max_wh) max_wh += StepWaveHeight min_wh += StepWaveHeight os.system('rm -f %s' % t) if __name__ == '__main__': sys.exit(main())