Archive for October, 2013

Object-oriented Python Interfaces to NAG Routines

30/10/2013

There have been a number of tutorials about how to call a NAG Library from Python: for example, Mike Croucher’s pages over at Walking Randomly or some of the Python posts on the NAG blog. For this post here I’ve dabbled with a complementary example which I’ve developed mainly as a proof-of-concept for how object-oriented wrappers to a NAG Library might look. Potentially this is the kind of interface design that you might want to put on top of Mike Croucher’s ctypes set up. You could then customize it yourself to provide the kind of functionality you need in your own application.

Here’s the code:

#!/usr/bin/env python
"""
An example Python interface to a NAG (Fortran) Library, using ctypes.
"""
class NagFunctionError(Exception):
    "A base NAG exception class."
    def __init__(self, function_name, code, message=None):
        "Here, message - if supplied - should be list."
        self.code = code
        self.function_name = function_name
        self.message = message

    def __str__(self):
        if (self.message is None):
            error_message = []
        else:
            error_message = self.message
        error_message.append('ABNORMAL EXIT, code ' + str(self.code) +
                             ', from NAG function ' + self.function_name + '.')
        return '\n** ' + '\n** '.join(error_message)

class NagLicenceError(NagFunctionError):
    "A NAG exception class for a missing licence."
    def __init__(self, function_name, code):
        super(self.__class__, self).__init__(function_name, code,
                                             ['Function call not licensed.'])

class NagTypeError(NagFunctionError):
    "A NAG exception class for an invalid argument instance."
    def __init__(self, function_name, msg):
        super(self.__class__, self).__init__(function_name, -199, msg)

class NagFunction(object):
    "A root NAG function class."
    from ctypes import c_int, Structure
    # The default error-handling mode for each (Fortran) Library call:
    ifail = c_int(1)

    class Complex(Structure):
        "A class to emulate a complex value using ctypes."
        from ctypes import c_double
        _fields_ = [("re", c_double), ("im", c_double)]

    def check_fortran_ifail(self):
        """
        Takes action based on the Fortran ifail value returned by a NAG Fortran
        Library call.
        """
        if (self.ifail.value == 0):
            return
        elif (self.ifail.value == -399):
            raise NagLicenceError(self.__class__.__name__, self.ifail.value)
        else:
            raise NagFunctionError(self.__class__.__name__, self.ifail.value)

    def check_type(self,
                   argument,
                   argument_name,
                   expected_type):
        "Verifies that argument is an instance of expected_type."
        if (isinstance(argument, expected_type)):
            return
        raise NagTypeError(self.__class__.__name__,
                           ['Invalid type ' + argument.__class__.__name__ +
                            ' for ' + argument_name + '.',
                            'Expected type ' + expected_type.__name__ + '.'])

class NagRootsLambertWComplex(NagFunction):
    "A NAG class for the complex Lambert W function, c05bb."
    def __init__(self, naglib_h):
        super(self.__class__, self).__init__()
        self.fcn_h = naglib_h.c05bbf_
        self.fcn_h.restype = None

    def evaluate(self, branch, offset, z):
        "The wrapper to the NAG Library call."
        from ctypes import byref, c_double, c_int
        self.check_type(branch, 'branch', int)
        self.check_type(offset, 'offset', bool)
        self.check_type(z, 'z', complex)
        branch_f = c_int(branch)
        if (offset):
            offset_f = c_int(1)
        else:
            offset_f = c_int(0)
        z_f = self.Complex(z.real, z.imag)
        w_f = self.Complex()
        resid_f = c_double()
        self.fcn_h(byref(branch_f),
                   byref(offset_f),
                   byref(z_f),
                   byref(w_f),
                   byref(resid_f),
                   byref(self.ifail))
        self.check_fortran_ifail()
        return complex(w_f.re, w_f.im), resid_f.value

def c_example():
    """
    Calls the Lambert W wrapper for a range of input values and prints the
    results.
    """
    def format_complex(z):
        "Formats a complex z for output as an (a, b) pair."
        return ''.join(['(',
                        format(z.real, '15.8e'),
                        ', ',
                        format(z.imag, '15.8e'),
                        ')'])

    from ctypes import cdll
    import sys
    naglib = '/path/to/libnag_nag.so'
    naglib_h = cdll.LoadLibrary(naglib)
    branch = 0
    offset = False
    sys.stdout.write('Lambert W function values for branch = ' + str(branch) +
                     ', offset = ' + str(offset) + ':\n')
    for z in [
            complex(0.5, -1.0),
            complex(1.0, 2.3),
            complex(4.5, -0.1),
            complex(6.0, 6.0)
    ]:
        w, resid = NagRootsLambertWComplex(naglib_h).evaluate(branch, offset, z)
        sys.stdout.write('z = ' + format_complex(z) +
                         ', W(z) = ' + format_complex(w) +
                         '. residual = ' + format(resid, '15.8e') + '\n')

if (__name__=='__main__'):
    c_example()

The Nag*Error classes, as you can see, are independent of the platform and NAG Library language being used. They serve as suggestions for how you could map the mechanism of the Fortran ifail integer or C NagError structure to something Pythonic.

The core NagFunction class would presumably contain all the utility methods used when wrapping the NAG calls. Since the NAG routines are statically typed you may want to provide a unified NAG look-and-feel by checking for correctly-typed input with the check_type method. If not you could just rely on Python’s own checking in the ctypes constructors.

At this point you could perhaps break the provided Library functionality into its conceptual Chapters by using a class hierarchy to reflect the grouping. I’ve chosen to subclass from NagFunction without any further indirection.

The work for wrapping a Library call is done in an evaluate method in a class for the Library routine. For the Lambert W function in this example the wrapping amounts simply to massaging the method’s input to be more C like, calling the Library routine, checking its exit status and then returning the relevant output data as Python objects.

Source code: https://github.com/matcross/blog/tree/master/object-oriented-python-interfaces-to-nag-routines

Advertisements

Customizing svn diff

26/10/2013

A few years ago I wrote about how my discovery of colordiff had improved my svn diff experience. Since then I’ve made a number of other tweaks to better customize my diffing:

#!/bin/sh -u
exclude_file=/path/to/exclude/file
svn diff --no-diff-deleted --diff-cmd diff -x "--ignore-all-space --text --unified=0" $* | \
  filterdiff --verbose --exclude-from-file ${exclude_file} ${tmpfile} | \
  grep -v --file=${exclude_file} | colordiff | less -RS

I never want to see full diffs for deleted files, hence I give the --no-diff-deleted option to svn diff. (I would also use --no-diff-added if it existed!)

I also regularly diff artifacts archived from automatic jobs which build and test NAG Fortran and C Libraries. There are some differences there that, although they need to be archived for possible future reference, I never need to see on a day-to-day basis — diffs for non-repeatable RNGs, for example. I use filterdiff (from patchutils) to remove these from the svn diff output. I have the patterns to exclude listed in a separate file (pointed to by the exclude_file in the script). These patterns are shell wildcards, so look something like *g05kg*e.x for ignoring the results differences from the Fortran (respectively C) example program for the NAG non-repeatable RNG initializer g05kgf (respectively g05kgc).

A block of differences excluded by filterdiff leaves behind the separator Index: filename from svn diff, so in the script above a call to grep filters these out too. Unfortunately that means each exclusion needs to appear twice in the exclude_file: once as a shell wildcard for filterdiff and then once as a basic regular expression for grep: so as both *g05kg*e.x and g05kg.e\.x say. I haven’t yet worked out how to unify this. You also end up with orphaned ==== separator lines as well, but I don’t feel too inconvenienced by this: I just jump between the resulting diffs according to the presence of the ^Index separator for the diffs that haven’t been excluded. These diffs will come out colorized via colordiff and less -R as discussed in the older post.

Note also that Index lines for deleted files, filtered by svn diff --no-diff-deleted, still appear in the final output. I don’t exclude these at the grep step because I still like to see that a file has been deleted, without needing to see what has been deleted. And although the facility exists in the [helpers] section of ~/.subversion/config to set the script above to override Subversion’s diff implementation, I like to run my diffs by invoking this script explicitly and then I can still use plain svn diff as a sanity check or fallback, especially if filtering is not required or whitespace in the diff is significant.

Source code: https://github.com/matcross/blog/tree/master/customizing-svn-diff