Archive for the ‘Python’ Category

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

Python arborealis

02/01/2013

I’ve been implementing a k-ary circularly linked tree (that is, a tree where each node could have a previous, a next, an up, or a down node) in Python.

The page Trees – How to Think Like a Computer Scientist is a nice introduction to some of the background concepts.

To __init__ each node of my TreeNode class you just need to provide the node’s contents and its links

class TreeNode:
    
    def __init__(self,
                 body=None,
                 Prev=None,
                 Next=None,
                 Up=None,
                 Down=None):
        self.body = body
        self.Prev = Prev
        self.Next = Next
        self.Up = Up
        self.Down = Down

I thought hard about getting __repr__ and __str__ right for this class (see Built-in Functions – repr, for example), but __repr__ is a bit tricky if you want to avoid too much recursion. I wrote a small getRepr function for doing this on any class. Members that are instances of the parent are recursed to a depth of one.

def getRepr(self,
            depth=0):
    keys = list(self.__dict__.keys())
    keys.sort()

    fields = []

    for key in keys:

        if (isinstance(self.__dict__[key], self.__class__)):

            if (depth < 1):
                value = getRepr(self.__dict__[key], depth=depth+1)
            else:
                value = ''.join(['<',
                                 self.__class__.__name__,
                                 '(...)>'])

        else:
            value = str(self.__dict__[key])

        fields.append(''.join([key,
                               '=',
                               value]))

    return ''.join([self.__class__.__name__,
                    '(',
                    ', '.join(fields),
                    ')'])

Hence my TreeNode.__repr__ is then really simple, and my class’s __str__ just returns the node’s body

    def __repr__(self):
        return getRepr(self)

    def __str__(self):
        return str(self.body)

Traversing a node goes all the way down to the tip of each branch and then backtracks to the point that a sibling exists. By default my traverse method processes each node in print mode.

    def isLeaf(self):
        return (self.Down is None)

    def traverse(self,
                 process_node=None):
        depth = 0
        node = self

        while (node is not None):

            if (process_node is None):
                import sys
                sys.stdout.write(' '*depth*2 + str(node) + '\n')
            else:
                process_node(node)

            if (node.isLeaf()):
            
                while (node is not None and
                       node != self and
                       node.Next is None):
                    node = node.Up
                    depth = depth - 1

                if (node is None or
                    node == self):
                    break

                node = node.Next
            else:
                node = node.Down
                depth = depth + 1

When it comes to actually populating a tree, my specific application is unusual in that I have a tree already, which is output by an external C program. So I currently have no functions for creating (or deleting) trees – I just read in my external tree from a file and set the links for each node accordingly. As an example (albeit somewhat clunky) of setting up a small tree directly, here’s the tree for the pseudocode 7 = 1 + 2 * 3; end (with the correct operator precedence!)

root = TreeNode()
asgn = TreeNode('=')
root.Down = asgn; asgn.Up = root
lhs = TreeNode(7)
asgn.Down = lhs; lhs.Up = asgn
plus = TreeNode('+')
lhs.Next = plus; plus.Prev = lhs; plus.Up = asgn
one = TreeNode(1)
plus.Down = one; one.Up = plus
times = TreeNode('*')
one.Next = times; times.Prev = one; times.Up = plus
two = TreeNode(2)
times.Down = two; two.Up = times
three = TreeNode(3)
two.Next = three; three.Prev = two; three.Up = times
end = TreeNode('end')
asgn.Next = end; end.Prev = asgn
print(repr(root))
root.traverse()

giving

TreeNode(Down=TreeNode(Down=<TreeNode(...)>, Next=<TreeNode(...)>, Prev=None, Up=<TreeNode(...)>, body==),
  Next=None, Prev=None, Up=None, body=None)
None
  =
    7
    +
      1
      *
        2
        3
  end

Source code: https://github.com/matcross/blog/tree/master/python-arborealis