source: misc/tools/event_selection/verify.py @ 6982

#!/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())