python-logo

Wrapping an IMSL C function with a Python class using PyNL

on Mar 22, 16 • by Jennifer Locke • with No Comments

How to use the IMSL Python Numerical Library to access IMSL C functions from Python...

Home » Analytics and Data Mining » Wrapping an IMSL C function with a Python class using PyNL

As an IMSL user, you now have a clear path to access your favorite IMSL C functions from Python with the release of the IMSL Python Numerical Library (PyNL) 1.0. PyNL contains a set of Python wrappers for a subset of IMSL C functions and you can create additional wrappers to access all the other functions as needed. In this post, we’ll discuss the key steps in wrapping an IMSL C function using PyNL.

PyNL’s underlying architecture and integration with the IMSL C library makes adding new wrappers straightforward. PyNL also integrates with NumPy and SciPy, and is available with Python 2.7.8 and Python 3.4.2. For dependencies and available platforms, see our website.

Leveraging an existing PyNL wrapper

Since many of IMSL C’s 500+ APIs have common design features, the best way to create a new wrapper for an IMSL C function is to leverage an existing wrapper for a similar function, if possible. The current set of wrappers include a representative set of algorithms that serve as reference implementations for most API types. For example, to create a new wrapper for the function lin_sol_posdef(), we can reuse much of the existing wrapper for a similar function, lin_sol_gen(). The existing wrapper is a Python class named LU, so let’s walk through the steps.

Developing the wrapper

The recent white paper, How to access IMSL C functions from Python, describes wrapping an IMSL C function with a Python function. While Python functions work in most cases, sometimes it’s better to use a Python class so you can maintain state in the Python object, which we’ll show here. The complete wrapper code is provided at the end of this blog for closer analysis.

After you install PyNL, all wrapper code is available as Python source within the “imsl” package in the site-packages area of your Python installation. This source code is available to allow additional wrappers to be added to PyNL.

Setup

We recommend creating a new package within the PyNL area in your Python site-packages area to allow multiple user-written wrappers to coexist with the existing PyNL wrappers. Use the following steps to create a new package, where PYNL_AREA refers to the site-packages/imsl folder in your Python installation.

1. Create a new folder named PYNL_AREA/user_lib.

2. Add a file named __init__.py to the PYNL_AREA/user_lib folder, with the following contents:

"""Initialize imsl.user_lib package.""" 
from imsl.user_lib.cholesky import Cholesky 
__all__ = ('Cholesky',) 

3. Add the package name ‘user_lib’ to the file PYNL_AREA/__init__.py.

4. Expose the name of the IMSL C function to be wrapped in PYNL_AREA/_imsllib.py. For example, to expose imsl_d_lin_sol_posdef(), add the following line at the end of the set of other functions that are exposed:

 self.expose("imsl_d_lin_sol_posdef ", None, _err_check, "math") 

This sets up the exception handling of PyNL to throw Python exceptions when errors are detected within IMSL C.

Note: If you are wrapping a C/Stat function, the final argument to self.expose should be “stat” rather than “math“.

The wrapper

Again, the IMSL C API for lin_sol_posdef() is similar to the interface for lin_sol_gen(), so we’ll leverage the design and source code of the PyNL LU class to create a PyNL class for lin_sol_posdef(). Note that the existing PyNL LU class supports both real and complex versions of lin_sol_gen() but for the purposes of this posting, we’ll wrap just the real-valued IMSL C function imsl_d_lin_sol_posdef().

The PyNL LU class has the following characteristics:

• class LU(a)

Parameters

a: An NxN array containing the input matrix

Methods

cond() Compute the L1 norm condition number of the matrix
factor() Compute the pivoted LU factorization of the matrix
factor_full() Compute the pivoted LU factorization of the matrix
inv() Compute the inverse of the matrix
solve(b[, transpose])            Solve a system of linear equations using right-hand side b

These characteristics are similar to those required by the IMSL C function imsl_d_lin_sol_posdef(). Using the PyNL LU class as a basis, we can create a Python class with the following characteristics for imsl_d_lin_sol_posdef():

• class Cholesky(a)

Parameters

a: An NxN array containing the input matrix

Methods

cond() Compute the L1 norm condition number of the matrix
factor() Compute the Cholesky factorization of the matrix
inv() Compute the inverse of the matrix
solve(b)                                    Solve a system of linear equations using right-hand side b

This new class shares many of the methods of the LU class as well as the shape of the input array, a. It was created using the LU class source code (located in the site-packages/imsl/linalg folder in your Python installation) as a starting point. Note that the method LU.factor_full() is unique to the LU algorithm so is not needed for the Cholesky class.

Once we have defined the design of the new class, the next steps include:

• Validating the input data types
• Performing data conversion in preparation for calling the IMSL C function
• Setting up arguments and calling the IMSL C function

Let’s look at some code snippets to see how to accomplish these steps.

Validating input data types

Since IMSL C is a C library that expects data to be of a specific type, the Python wrapper determines whether the input data can be converted to the target precision. This excerpt is from the wrapper method Cholesky.__init__:

 
        # attempt to promote matrix a to a compatible type.
        common_type = _numpy.promote_types(_numpy.float64, self._a.dtype)
        self._a = _numpy.asarray(self._a, dtype=common_type)

        if (not _numpy.issubdtype(self._a.dtype, _numpy.float64)):
            raise ValueError("array type {} not supported".format(
                             self._a.dtype.name))

Performing data conversion in preparation for calling the IMSL C function

Again, since IMSL C expects specific datatypes, all data passed to the IMSL C function is converted to the desired type using the NumPy method asarray. For example, in the wrapper method Cholesky.solve:

      b = _numpy.asarray(b, dtype=self._a.dtype, order='C')

Setting up arguments to pass to the IMSL C function

The Python wrapper creates an argument list and passes it to the IMSL C function. This argument list corresponds to the set of required and optional arguments necessary to perform the desired action. Note that IMSL C requires a terminating zero at the end of the optional argument list. Consider the method Cholesky.inv:

            args = []

            row_dim = self._a.shape[0]

            args.append(row_dim)
            args.append(self._a.ctypes.data_as(_ctypes.c_void_p))
            args.append(0)  # b is ignored

            self._inverse = _numpy.empty([row_dim, row_dim],
                                         dtype=self._a.dtype)

            args.append(_constants.IMSL_INVERSE_ONLY)
            args.append(_constants.IMSL_INVERSE_USER)
            args.append(self._inverse.ctypes.data_as(_ctypes.c_void_p))

            args.append(0)

            func = _imsllib.imsllib.imsl_d_lin_sol_posdef
            func(*args)

Sample usage

Let’s see how you can use the Python interface to the IMSL C function imsl_d_lin_sol_posdef() to solve a system of three linear equations with a symmetric positive definite coefficient matrix. We’ll also demonstrate how to compute the inverse of the coefficient matrix, the L1 norm condition number, and Cholesky factors. The equations are listed below:

x1 − 3x2 + 2x3 = 27

−3x1 + 10x2 − 5x3 = −78

2x1 − 5x2 + 6x3 = 64

 """ Cholesky class example. """
import imsl.user_lib.cholesky as chol
import numpy as np

a = [[1.0, -3.0, 2.0], [-3.0, 10.0, -5.0], [2.0, -5.0, 6.0]]
b = [27.0, -78.0, 64.0]
chol = chol.Cholesky(a)
x = chol.solve(b)
print("Solution:")
print(x)
print("")

inv = chol.inv()
print("Inverse:")
print(inv)
print("")

# Compute a*inv(a) (=I)
prod = np.dot(a, inv)
print("a*inv(a):")
print(prod)
print("")

print("L1 condition number:")
print(chol.cond())
print("")

L = chol.factor()
print("Cholesky factors:")
print(L)
print("")

# Compute a-L*trans(L) (=0)
for i in range(x.size):
    L[i, i+1:] = 0

print("a - L*trans(L):")
print(a - np.dot(L, np.transpose(L)))

Output

Solution:
[ 1. -4.  7.]

Inverse:
[[ 35.   8.  -5.]
 [  8.   2.  -1.]
 [ -5.  -1.   1.]]

a*inv(a):
[[ 1.  0.  0.]
 [ 0.  1.  0.]
 [ 0.  0.  1.]]

L1 condition number:
673.838709677

Cholesky factors:
[[ 1. -3.  2.]
 [-3.  1.  1.]
 [ 2.  1.  1.]]

a - L*trans(L):
[[ 0.  0.  0.]
 [ 0.  0.  0.]
 [ 0.  0.  0.]]

Complete wrapper code for imsl_d_lin_sol_posdef()

The complete wrapper for the IMSL C function imsl_d_lin_sol_posdef(), using a Python class named Cholesky is provided below.

"""Cholesky factorization related class, methods."""
import ctypes as _ctypes
import numpy as _numpy

import imsl._constants as _constants
import imsl._imsllib as _imsllib


class Cholesky():

    """Solve a real symmetric positive definite system of linear equations."""

    def __init__(self, a):
        """Instantiate Cholesky class."""
        self._a = a
        self._factor_factor = None
        self._cond = None
        self._inverse = None

        if self._a is None:
            raise TypeError("None not supported")

        self._a = _numpy.array(a, order='C')

        # attempt to promote matrix a to a compatible type.
        common_type = _numpy.promote_types(_numpy.float64, self._a.dtype)
        self._a = _numpy.asarray(self._a, dtype=common_type)

        if (not _numpy.issubdtype(self._a.dtype, _numpy.float64)):
            raise ValueError("array type {} not supported".format(
                             self._a.dtype.name))

        if self._a.ndim != 2:
            raise ValueError("array of dimension {} not"
                             " supported".format(self._a.ndim))
        if self._a.size == 0:
            raise ValueError("empty array not supported")
        if (self._a.shape[0] != self._a.shape[1]):
            raise ValueError("input matrix must be square")

    def solve(self, b):
        r"""Solve a real symmetric positive definite system of linear equations 
            using right-hand side `b`."""
        if b is None:
            raise TypeError("None not supported")

        b = _numpy.asarray(b, dtype=self._a.dtype, order='C')

        if b.ndim != 1:
            raise ValueError("array of dimension {} not"
                             " supported".format(b.ndim))

        if b.size != self._a.shape[0]:
            raise ValueError("dependent variable values length ({}) does not"
                             " match coefficient matrix row count"
                             " ({})".format(b.size, self._a.shape[0]))

        args = []

        row_dim = self._a.shape[0]

        args.append(row_dim)
        args.append(self._a.ctypes.data_as(_ctypes.c_void_p))
        args.append(b.ctypes.data_as(_ctypes.c_void_p))

        solution = _numpy.empty_like(b)
        args.append(_constants.IMSL_RETURN_USER)
        args.append(solution.ctypes.data_as(_ctypes.c_void_p))

        if self._factor_factor is None:
            self._factor_factor = _numpy.empty([row_dim, row_dim],
                                               dtype=self._a.dtype)
        else:
            args.append(_constants.IMSL_SOLVE_ONLY)

        args.append(_constants.IMSL_FACTOR_USER)
        args.append(self._factor_factor.ctypes.data_as(_ctypes.c_void_p))

        if self._cond is None:
            args.append(_constants.IMSL_CONDITION)
            cond = _ctypes.c_double()
            args.append(_ctypes.byref(cond))

        args.append(0)

        func = _imsllib.imsllib.imsl_d_lin_sol_posdef
        func(*args)

        if self._cond is None:
            self._cond = cond.value

        return solution

    def factor(self):
        """Compute the Cholesky factorization of the matrix."""
        if (self._factor_factor is None or self._cond is None):
            args = []

            row_dim = self._a.shape[0]

            args.append(row_dim)
            args.append(self._a.ctypes.data_as(_ctypes.c_void_p))
            args.append(0)  # b is ignored

            self._factor_factor = _numpy.empty([row_dim, row_dim],
                                               dtype=self._a.dtype)

            args.append(_constants.IMSL_FACTOR_ONLY)
            args.append(_constants.IMSL_FACTOR_USER)
            args.append(self._factor_factor.ctypes.data_as(_ctypes.c_void_p))

            args.append(_constants.IMSL_CONDITION)
            cond = _ctypes.c_double()
            args.append(_ctypes.byref(cond))

            args.append(0)

            func = _imsllib.imsllib.imsl_d_lin_sol_posdef
            func(*args)

            self._cond = cond.value

        return (_numpy.copy(self._factor_factor))


    def inv(self):
        """Compute the inverse of the matrix."""
        if self._inverse is None:
            args = []

            row_dim = self._a.shape[0]

            args.append(row_dim)
            args.append(self._a.ctypes.data_as(_ctypes.c_void_p))
            args.append(0)  # b is ignored

            self._inverse = _numpy.empty([row_dim, row_dim],
                                         dtype=self._a.dtype)

            args.append(_constants.IMSL_INVERSE_ONLY)
            args.append(_constants.IMSL_INVERSE_USER)
            args.append(self._inverse.ctypes.data_as(_ctypes.c_void_p))

            args.append(0)

            func = _imsllib.imsllib.imsl_d_lin_sol_posdef
            func(*args)

        return _numpy.copy(self._inverse)

    def cond(self):
        """Compute the L_1 norm condition number of the matrix."""
        if self._cond is None:
            self.factor()

        return self._cond

Related Posts

Leave a Reply

Your email address will not be published. Required fields are marked *

Scroll to top