Python/NAG Part 4 – Structures

April 16th, 2009 | Categories: NAG Library, programming, python | Tags:

This is part 4 of a series of articles devoted to demonstrating how to call the Numerical Algorithms Group (NAG) C library from Python. Click here for the index to this series.

The NAG C library makes extensive use of C structures throughout its code base and so we will need to know how to represent them in Python.  For example, the NAG routine d01ajc is a general purpose one dimensional numerical integrator and one of it’s arguments is a structure called Nag_QuadProgress with the following prototype (found in the nag_types.h file that comes with the NAG C library)

  typedef struct
{
Integer num_subint;
Integer fun_count;
double *sub_int_beg_pts;
double *sub_int_end_pts;
double *sub_int_result;
double *sub_int_error;
} Nag_QuadProgress;

following the ctypes documentation on Structures and Unions, I came up with the following Python representation of this structure and it seems to work just fine.

class Nag_QuadProgress(Structure):
_fields_ = [("num_subint", c_int),
("funcount", c_int),
("sub_int_beg_pts",POINTER(c_double) ),
("sub_int_end_pts",POINTER(c_double)),
("sub_int_result",POINTER(c_double)  ),
("sub_int_error",POINTER(c_double)  )
]


This is the only extra information we need in order to write a Python program that calculates the numerical approximation to an integral using the NAG routine d01ajc. Let’s calculate the integral of the function 4/(1.0+x*x) over the interval [0,1]. The code is as follows:

#!/usr/bin/env python
#Example 3 using NAG C library and ctypes
#Mike Croucher - April 2008
#Concept introduced : Creating and using Python versions of structures such as Nag_QuadProgress
from ctypes import *

d01ajc = libnag.d01ajc
d01ajc.restype=None

_fields_ = [("num_subint", c_int),
("funcount", c_int),
("sub_int_beg_pts",POINTER(c_double) ),
("sub_int_end_pts",POINTER(c_double)),
("sub_int_result",POINTER(c_double)  ),
("sub_int_error",POINTER(c_double)  )
]

#define python function
def py_func(x):
return 4./(1.0+x*x)

#create the prototype for the c callback function
C_FUNC = CFUNCTYPE(c_double,c_double )

#now create the c callback function
c_func = C_FUNC(py_func)

#set up the problem
a = c_double(0.0)
b = c_double(1.0)
nlimit = c_int(100)
epsr = c_double(0.00001)
ifail = c_int(0)

max_num_subint=c_int(200)
nlimit = c_int(0)
result=c_double(0.0)
epsabs= c_double(0.0)
epsrel=c_double(0.0001);
abserr=c_double(0.0);

d01ajc(c_func,a,b,epsabs,epsrel,max_num_subint,byref(result) ,byref(abserr),byref(qp),ifail )
print "Approximation of integral= %s" %(result.value)
print "Number of function evaluations=%s" % qp.funcount
print "Number of subintervals used= %s" % qp.num_subint


Copy and paste (or download) this to a file called d01ajc.py, make it executable and run it by typing the following at a bash prompt.
 chmod +x ./d01ajc.py ./d01ajc.py 
If all has gone well then you should get following result

Approximation of integral= 3.14159265359
Number of function evaluations=21
Number of subintervals used= 1


Note that in order to create an instance of the Nag_QuadProgress structure I just do
 qp=Nag_QuadProgress() 
It is also useful to compare the C prototype of d01ajc
 void nag_1d_quad_gen (double (*f)(double x), double a, double b, double epsabs, double epsrel, Integer max_num_subint, double *result, double *abserr, Nag_QuadProgress *qp, NagError *fail) 
with the way I called it in Python
 d01ajc(c_func,a,b,epsabs,epsrel,max_num_subint,byref(result) ,byref(abserr),byref(qp),ifail ) 
Note how pointer arguments such as
 double *abserr 
are called as follows
 byref(abserr) 

Yet again I have cheated with the NagError parameter and just used c_int(0) instead of a proper NagError structure. NagError will be treated properly in a future article.