1#!/usr/bin/env @PYTHON_EXECUTABLE@
2"""
3Description: Numerically compare the contents of two Siconos mechanics-IO HDF5 simulation files.
4
5Sole output on standard output shall be the maximum difference between
6the columns selected for comparison.  Return code 0 will be returned
7if maximum difference is less than the threshold, otherwise 1 will be
8returned if greater, or 2 if the files could not be opened or do not
9contain the necessary tables.  Non-matching time indices shall be
10considered as failing the comparison.
11"""
12
13# Lighter imports before command line parsing
14from __future__ import print_function
15import os, sys, argparse, re
16
17parser = argparse.ArgumentParser(
18    description = __doc__)
19parser.add_argument('fns_in', metavar='<file>', type=str, nargs=2,
20                    help = 'input files (HDF5)')
21parser.add_argument('--columns','-c', metavar='<colspec>', type=str,
22                    help = """
23list of columns to include in the comparison (0-indexed,
24comma-separated, may be prefixed with table name, and multiple lists
25separated by semicolon, e.g., "dynamic:0,1,2;cf:0,1,2"; default
26table="dynamic").  See Siconos HDF5 file format docs for meaning of
27columns.""")
28parser.add_argument('--start', metavar='<time>', type=float,
29                    help = 'start time in seconds from beginning the comparison')
30parser.add_argument('--end', metavar='<time>', type=float,
31                    help = 'end time in seconds end the comparison')
32parser.add_argument('--interval', metavar='<time>', type=float,
33                    help = 'amount of time after start to compare')
34parser.add_argument('--threshold', metavar='<float>', type=float, default=1e-14,
35                    help = 'threshold for comparison')
36parser.add_argument('-V','--version', action='version',
37                    version='@SICONOS_VERSION@')
38
39if __name__ == '__main__':
40    args = parser.parse_args()
41    columns = args.columns
42    tablenames = set()
43    if args.columns is None:
44        columns = [('dynamic', None)]
45        tablenames.update(['dynamic'])
46    else:
47        columns = []
48        for spec in args.columns.split(';'):
49            if ':' in spec:
50                table, cols = spec.split(':')
51            else:
52                table = 'dynamic'
53                cols = spec
54            if ',' in cols:
55                columns += [(table, int(col)) for col in cols.split(',')]
56            else:
57                columns += [(table, None)]
58            tablenames.update([table])
59
60# Heavier imports after command line parsing
61import numpy as np
62import h5py
63
64def verify_tables(tablenames, columns, io1, io2):
65    """Verify files contain tables and tables have same columns."""
66    for table in tablenames:
67        if not table in io1['data']:
68            print('File "{}" does not have table "{}".'.format(
69                args.fns_in[0], table), file=sys.stderr)
70            sys.exit(2)
71        if not table in io2['data']:
72            print('File "{}" does not have table "{}".'.format(
73                args.fns_in[1], table), file=sys.stderr)
74            sys.exit(2)
75        for c in columns:
76            if c[1] is None:
77                continue
78            if c[0] == table and c[1] >= io1['data'][table].shape[1]:
79                print('Table "{}" in file "{}" does not have specified column {}.'
80                      .format(table, args.fns_in[0], c[1]), file=sys.stderr)
81                sys.exit(2)
82            if c[0] == table and c[1] >= io2['data'][table].shape[1]:
83                print('Table "{}" in file "{}" does not have specified column {}.'
84                      .format(table, args.fns_in[1], c[1]), file=sys.stderr)
85                sys.exit(2)
86        if io1['data'][table].shape[1] != io2['data'][table].shape[1]:
87            print('Tables "{}" do not have same number of columns in each file.'
88                  .format(table), file=sys.stderr)
89            sys.exit(2)
90
91def compare_tables(tablenames, columns, io1, io2):
92    """Compare tables and columns given."""
93    maxdiff = 0.0
94    for table in tablenames:
95        for c in columns:
96            if c[0]!=table:
97                continue
98            t1 = io1['data'][table]
99            t2 = io2['data'][table]
100
101            # start/end/interval times
102            S1, S2 = 0, 0
103            E1, E2 = t1.shape[0]-1, t2.shape[0]-1
104            if args.start is not None:
105                S1 = np.searchsorted(t1[:,0], args.start)
106                S2 = np.searchsorted(t2[:,0], args.start)
107            if args.end is not None:
108                E1 = np.searchsorted(t1[:,0], args.end)
109                E2 = np.searchsorted(t2[:,0], args.end)
110            elif args.interval is not None:
111                E1 = np.searchsorted(t1[:,0], args.start + args.interval)
112                E2 = np.searchsorted(t2[:,0], args.start + args.interval)
113            for t,n,s,T,fn in [(S1,args.start,'Start',t1,args.fns_in[0]),
114                               (S2,args.start,'Start',t2,args.fns_in[1]),
115                               (E1,args.end,'End',t1,args.fns_in[0]),
116                               (E2,args.end,'End',t2,args.fns_in[1])]:
117                if t >= T.shape[0]:
118                    print('{} time {} beyond the end of table "{}" for file "{}".'
119                          .format(s, n, table, fn), file=sys.stderr)
120                    sys.exit(2)
121
122            # TODO: we assume same sampling rate for now, later,
123            # support linear interpolation?
124            if (E1-S1) != (E2-S2):
125                print(('Tables "{}" do not have same number of rows in each file '
126                       +'for requested range.').format(table), file=sys.stderr)
127                sys.exit(2)
128
129            # Compare one column at a time
130            cols = [c[1]]
131            if c[1] is None:
132                cols = range(t1.shape[1])
133            for j in cols:
134                d = np.abs(t1[S1:E1,j] - t2[S2:E2,j]).max()
135                if d > maxdiff:
136                    maxdiff = d
137    print(maxdiff)
138    return int(maxdiff >= args.threshold)
139
140if __name__ == '__main__':
141    with h5py.File(args.fns_in[0], mode='r') as io1:
142        with h5py.File(args.fns_in[1], mode='r') as io2:
143            verify_tables(tablenames, columns, io1, io2)
144            sys.exit(compare_tables(tablenames, columns, io1, io2))
145