Walking Randomly

This is part 5 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 Optimisation chapter of the NAG C library is one of the most popular chapters in my experience.  It is also one of the most complicated from a programming point of view but if you have been working through this series from the beginning then you already know all of the Python techniques you need.  From this point on it’s just a matter of increased complexity.

As always, however, I like to build the complexity up slowly so let’s take a look at possibly the simplest NAG optimisation routine – e04abc.  This routine minimizes a function of one variable, using function values only.  In other words – you don’t need to provide it with derivative information which simplifies the programming at the cost of performance.

The example Python code I am going to give will solve a very simple local optimisation problem – The minimisation of the function sin(x)/x in the range [3.5,5.0] or, put another way, what are the x and y values of the centre of the red dot in the picture below.

The NAG/Python code to do this is given below (or you can download it here).

#!/usr/bin/env python
#Example 4 using NAG C library and ctypes - e04abc
#Mike Croucher - April 2008
#Concepts : Communication callback functions that use the C prototype void funct(double, double *,Nag_comm *)

from ctypes import *
from math import *
libnag = cdll.LoadLibrary("/opt/NAG/cllux08dgl/lib/libnagc_nag.so.8")
e04abc = libnag.e04abc
e04abc.restype=None

class Nag_Comm(Structure):
	_fields_ = [("flag", c_int),
                ("first", c_int),
 		("last", c_int),
		("nf",c_int),
		("it_prt",c_int),
		("it_maj_prt",c_int),
		("sol_sqp_prt",c_int),
		("sol_prt",c_int),
		("rootnode_sol_prt",c_int),
		("node_prt",c_int),
		("rootnode_prt",c_int),
		("g_prt",c_int),
		("new_lm",c_int),
		("needf",c_int),
		("p",c_void_p),
		("iuser",POINTER(c_int) ),
		("user",POINTER(c_double) ),
		("nag_w",POINTER(c_double) ),
		("nag_iw",POINTER(c_int) ),
		("nag_p",c_void_p),
		("nag_print_w",POINTER(c_double) ),
		("nag_print_iw",POINTER(c_int) ),
		("nrealloc",c_int)
		]

#define python objective function
def py_obj_func(x,ans,comm):
  ans[0] = sin(x) / x
  return

#create the type for the c callback function
C_FUNC = CFUNCTYPE(None,c_double,POINTER(c_double),POINTER(Nag_Comm))

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

#Set up the problem parameters
comm = Nag_Comm()
#e1 and e2 control the solution accuracy.  Setting them to 0 chooses the default values.
#See the documentation for e04abc for more details
e1 = c_double(0.0)
e2 = c_double(0.0)
#a and b define our starting interval
a = c_double(3.5)
b = c_double(5.0)
#xt will contain the x value of the minimum
xt=c_double(0.0)
#f will contain the function value at the minimum
f=c_double(0.0)
#max_fun is the maximum number of function evaluations we will allow
max_fun=c_int(100)
#Yet again I cheat for fail parameter.
fail=c_long(0)

e04abc(c_func,e1,e2,byref(a),byref(b),max_fun,byref(xt),byref(f),byref(comm),fail)

print "The minimum is in the interval  [%s,%s]" %(a.value,b.value)
print "The estimated result = %s" % xt.value
print "The function value at the minimum is %s " % f.value
print "%s evaluations were required" % comm.nf

If you run this code then you should get the following output

The minimum is in the interval  [4.4934093959,4.49340951173]
The estimated result = 4.49340945381
The function value at the minimum is -0.217233628211
10 evaluations were required

IFAIL oddness

My original version of this code used the trick of setting the NAG IFAIL parameter to c_int(0) to save me from worrying about the NagError structure (this is first explained in part 1) but when this example was tested it didn’t work on 64bit machines although it worked just fine on 32bit machines.  Changing it to c_long(0) (as I have done in this example)  made it work on both 32 and 64bit machines.

What I find odd is that in all my previous examples, c_int(0) worked just fine on both architectures.  Eventually I will cover how to use the IFAIL parameter ‘properly’ but for now it is sufficient to be aware of this issue.

Scary Structures

Although this is the longest piece of code we have seen so far in this series, it contains no new Python/ctypes techniques at all. The only major difficulty I faced when writing it for the first time was making sure that I got the specification of the Nag_Comm structure correct by carefully working though the C specification of it in the nag_types.h file (included with the NAG C library installation) and transposing it to Python.

One of the biggest hurdles you will face when working with the NAG C library and Python is getting these structure definitions correct and some of them are truly HUGE. If you can get away with it then you are advised to copy and paste someone else’s definitions rather than create your own from scratch. In a later part of this series I will be releasing a Python module called PyNAG which will predefine some of the more common structures for you.

Callback functions with Communication

The objective function (that is, the function we want to minimise) is implemented as a callback function but it is a bit more complicated than the one we saw in part 3.  This NAG routine expects the C prototype of your objective function to look like this:

void funct(double xc, double *fc, Nag_Comm *comm)

• xc is an input argument and is the point at which the value of f(x) is required.
• fc is an output argument and is the value of the function f at the current point x.  Note that it is a pointer to double.
• Is a pointer to a structure of type Nag_Comm.  Your function may or may not use this directly but the NAG routine almost certainly will.  It is through this structure that the NAG routine can communicate various things to and from your function as and when is necessary.

so our C-types function prototype looks like this

C_FUNC = CFUNCTYPE(None,c_double,POINTER(c_double),POINTER(Nag_Comm))

whereas the python function itself is

#define python objective function
def py_obj_func(x,ans,comm):
  ans[0] = sin(x) / x
  return

Note that we need to do

ans[0] = sin(x) / x

rather than simply

ans = sin(x) / x

since ans is a pointer to double rather than just a double.

In this particular example I have only used the simplest feature provided by the Nag_Comm structure and that is to ask it how many times the objective function was called.

print "%s evaluations were required" % comm.nf

That’s it for this part of the series. Next time I’ll be looking at a more advanced optimisation routine. Feel free to post any comments, questions or suggestions.