1#! /usr/bin/env python
2
3"""
4Checks computed forces from an apbs run
5"""
6
7import sys, re
8from apbs_logger import Logger
9
10error_tolerance = 1e-6
11
12class PolarForce:
13    """
14    Exctracts and compares computations of polar forces
15    """
16
17    # A crazy regex pattern used to match label/value sets
18    pattern=r'\s+(?P<label>[a-zA-Z]+)\s+(?P<x>[+-]?\d\.\d+E[+-]\d+)\s+(?P<y>[+-]?\d\.\d+E[+-]\d+)\s+(?P<z>[+-]?\d\.\d+E[+-]\d+)'
19
20    def __init__( self, label, x, y, z ):
21        """
22        Constructs a polar force result from supplied values
23        """
24        self.label = label
25        self.x = x
26        self.y = y
27        self.z = z
28
29    def __init__( self, line ):
30        """
31        Extracts ploar force results from a file at a given line
32        """
33        m = re.search( self.pattern, line )
34        self.label = m.group( 'label' )
35        self.x = float( m.group( 'x' ) )
36        self.y = float( m.group( 'y' ) )
37        self.z = float( m.group( 'z' ) )
38
39    def diff( self, other ):
40        """
41        Compares the value of two polar force field results
42        """
43
44        diff_dict = {}
45        for d in ( 'x', 'y', 'z' ):
46            diff_dict[ d ] = abs( getattr( self, d ) - getattr( other, d ) )
47        return diff_dict
48
49    def __repr__( self ):
50        return "PolarForce{ label:%s, x:%g, y:%g, z:%g}\n" % ( self.label, self.x, self.y, self.z )
51
52    def short( self ):
53        return self.label
54
55
56
57class ApolarForce( PolarForce ):
58    """
59    Exctracts and compares computations of apolar forces
60    """
61
62    # A crazy regex pattern used to match label/value sets
63    pattern=r'\s+(?P<label>[a-zA-Z]+)\s+(?P<atom>\w+)\s+(?P<x>[+-]?\d\.\d+E[+-]\d+)\s+(?P<y>[+-]?\d\.\d+E[+-]\d+)\s+(?P<z>[+-]?\d\.\d+E[+-]\d+)'
64
65    def __init__( self, label, atom, x, y, z ):
66        """
67        Constructs an apolar force result from supplied values
68        """
69        super( ApolarForce, self ).__init__( self, x, y, z )
70        self.label = label
71
72    def __init__( self, line ):
73        """
74        Extracts aploar force results from a file at a given line
75        """
76        m = re.search( self.pattern, line )
77        self.label = m.group( 'label' )
78        self.atom = m.group( 'atom' )
79        self.x = float( m.group( 'x' ) )
80        self.y = float( m.group( 'y' ) )
81        self.z = float( m.group( 'z' ) )
82
83    def __rpr__( self ):
84        return "ApolarForce{ label:%s, atom:%s, x:%g, y:%g, z:%g}\n" % ( self.label, self.atom, self.x, self.y, self.z )
85
86    def short( self ):
87        return '%s for %s' % ( self.label, self.atom )
88
89
90
91
92def extract_forces( force_class, lines, start_pattern, ):
93    """
94    Extracts force results
95    """
96    force_dict = {}
97    in_section = False
98    in_forces = False
99    start_line = -1
100    end_line = -1
101    for line_number, line_text in enumerate(lines):
102        if not in_section:
103            if line_text.startswith( start_pattern ):
104                in_section = True
105        if in_section and not in_forces:
106            if re.search( force_class.pattern, line_text ):
107                in_forces = True
108                start_line = line_number
109        if in_section and in_forces:
110            if not re.search( force_class.pattern, line_text ):
111                end_line = line_number
112                break
113    return parse_forces( force_class, lines[ start_line : end_line ] )
114
115
116
117def parse_forces( force_class, lines ):
118    force_dict = {}
119    for line in lines:
120        force_item = force_class( line )
121        force_dict[ force_item.label ] = force_item
122    return force_dict
123
124
125
126def compare_force_dicts( test_force_dict, true_force_dict ):
127    """
128    Compares force dictionaries
129    """
130
131    for force_key in test_force_dict.keys():
132        test_force = test_force_dict[ force_key ]
133        true_force = true_force_dict[ force_key ]
134        diff_dict = test_force.diff( true_force )
135
136        for ( diff_key, diff_value ) in diff_dict.items():
137            test_value = getattr( test_force, diff_key )
138            true_value = getattr( true_force, diff_key )
139
140            if diff_value == 0.0:
141                logger.message( '*** Comparison %s in %s PASSED ***' % ( test_force.short(), diff_key ) )
142                logger.log( 'Comparison %s in %s PASSED (%g)' % ( test_force.short(), diff_key, test_value ) )
143
144            elif diff_value < error_tolerance:
145                logger.message( '*** Comparison %s in %s PASSED (with rounding error - see log)***' % ( test_force.short(), diff_key ) )
146                logger.log( 'Comparison %s in %s PASSED within error (%g; expected %g)' %( test_force.short(), diff_key, test_value, true_value ) )
147
148            else:
149                logger.message( '*** Comparison %s in %s FAILED ***' % ( test_force.short(), diff_key ) )
150                logger.message( '   APBS returned %g' % test_value )
151                logger.message( '   Expected result is %g (difference of: %g)' % ( true_value, diff_value ) )
152                logger.log( 'Comparison %s in %s FAILED (%g; expected %g)' %( test_force.short(), diff_key, test_value, true_value ) )
153
154
155
156def check_forces( input_file, polar_file, apolar_file, logger ):
157
158    logger.both( "Checking forces for input file %s" % input_file )
159
160    f = None
161    try:
162        f = open( input_file, 'r' )
163    except IOError:
164        print >> sys.stderr, "Couldn't read from forces file %s" % input_file
165    input_lines = f.readlines()
166
167    test_polar_force_dict = extract_forces( PolarForce, input_lines, 'print force' )
168    test_apolar_force_dict = extract_forces( ApolarForce, input_lines, 'print APOL force' )
169
170    try:
171        f = open( polar_file, 'r' )
172    except IOError:
173        print >> sys.stderr, "Couldn't read from forces file %s" % input_file
174    input_lines = f.readlines()
175    true_polar_force_dict = parse_forces( PolarForce, input_lines )
176
177    try:
178        f = open( apolar_file, 'r' )
179    except IOError:
180        print >> sys.stderr, "Couldn't read from forces file %s" % input_file
181    input_lines = f.readlines()
182    true_apolar_force_dict = parse_forces( ApolarForce, input_lines )
183
184    logger.both( "Checking Polar Forces" )
185    compare_force_dicts( test_polar_force_dict, true_polar_force_dict )
186
187    logger.both( "Checking Apolar Forces" )
188    compare_force_dicts( test_apolar_force_dict, true_apolar_force_dict )
189
190
191
192def test():
193    l = open( 'forces.log', 'w' )
194    logger = Logger( sys.stderr, l )
195    check_forces( 'apbs-forces.out', 'polarforces', 'apolarforces', logger )
196
197
198
199if __name__ == '__main__':
200    print >> sys.stderr, "The python source file %s is a module and not runnable" % sys.argv[ 0 ]
201    sys.exit( 1 )
202