{"id":1058,"date":"2009-04-24T16:22:38","date_gmt":"2009-04-24T15:22:38","guid":{"rendered":"http:\/\/www.walkingrandomly.com\/?p=1058"},"modified":"2009-05-13T19:32:03","modified_gmt":"2009-05-13T18:32:03","slug":"pythonnag-part-5-callback-functions-with-communication","status":"publish","type":"post","link":"https:\/\/walkingrandomly.com\/?p=1058","title":{"rendered":"Python\/NAG Part 5 – Callback functions with communication."},"content":{"rendered":"

This is part 5 of a series of articles devoted to demonstrating how to call the Numerical Algorithms Group (NAG) C library<\/a> from Python.\u00a0 Click here for the index<\/a> to this series.<\/p>\n

The Optimisation chapter of the NAG C library<\/a> is one of the most popular chapters in my experience.\u00a0 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.\u00a0 From this point on it’s just a matter of increased complexity.<\/p>\n

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<\/a>.\u00a0 This routine minimizes a function of one variable, using function values only.\u00a0 In other words – you don’t need to provide it with derivative information which simplifies the programming at the cost of performance.<\/p>\n

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.<\/p>\n

\"Minimum<\/p>\n

The NAG\/Python code to do this is given below (or you can download it here<\/a>).<\/p>\n

#!\/usr\/bin\/env python\r\n#Example 4 using NAG C library and ctypes - e04abc\r\n#Mike Croucher - April 2008\r\n#Concepts : Communication callback functions that use the C prototype void funct(double, double *,Nag_comm *)<\/strong><\/pre>\n
from ctypes import *\r\nfrom math import *\r\nlibnag = cdll.LoadLibrary(\"\/opt\/NAG\/cllux08dgl\/lib\/libnagc_nag.so.8\")\r\ne04abc = libnag.e04abc\r\ne04abc.restype=None<\/strong><\/pre>\n
class Nag_Comm(Structure):\r\n\t_fields_ = [(\"flag\", c_int),\r\n                (\"first\", c_int),\r\n \t\t(\"last\", c_int),\r\n\t\t(\"nf\",c_int),\r\n\t\t(\"it_prt\",c_int),\r\n\t\t(\"it_maj_prt\",c_int),\r\n\t\t(\"sol_sqp_prt\",c_int),\r\n\t\t(\"sol_prt\",c_int),\r\n\t\t(\"rootnode_sol_prt\",c_int),\r\n\t\t(\"node_prt\",c_int),\r\n\t\t(\"rootnode_prt\",c_int),\r\n\t\t(\"g_prt\",c_int),\r\n\t\t(\"new_lm\",c_int),\r\n\t\t(\"needf\",c_int),\r\n\t\t(\"p\",c_void_p),\r\n\t\t(\"iuser\",POINTER(c_int) ),\r\n\t\t(\"user\",POINTER(c_double) ),\r\n\t\t(\"nag_w\",POINTER(c_double) ),\r\n\t\t(\"nag_iw\",POINTER(c_int) ),\r\n\t\t(\"nag_p\",c_void_p),\r\n\t\t(\"nag_print_w\",POINTER(c_double) ),\r\n\t\t(\"nag_print_iw\",POINTER(c_int) ),\r\n\t\t(\"nrealloc\",c_int)\r\n\t\t]\r\n<\/strong><\/pre>\n
#define python objective function\r\ndef py_obj_func(x,ans,comm):\r\n  ans[0] = sin(x) \/ x\r\n  return\r\n<\/strong><\/pre>\n
#create the type for the c callback function\r\nC_FUNC = CFUNCTYPE(None,c_double,POINTER(c_double),POINTER(Nag_Comm))\r\n<\/strong><\/pre>\n
#now create the c callback function\r\nc_func = C_FUNC(py_obj_func)\r\n<\/strong><\/pre>\n
#Set up the problem parameters\r\ncomm = Nag_Comm()\r\n#e1 and e2 control the solution accuracy.  Setting them to 0 chooses the default values.\r\n#See the documentation for e04abc for more details\r\ne1 = c_double(0.0)\r\ne2 = c_double(0.0)\r\n#a and b define our starting interval\r\na = c_double(3.5)\r\nb = c_double(5.0)\r\n#xt will contain the x value of the minimum\r\nxt=c_double(0.0)\r\n#f will contain the function value at the minimum\r\nf=c_double(0.0)\r\n#max_fun is the maximum number of function evaluations we will allow\r\nmax_fun=c_int(100)\r\n#Yet again I cheat for fail parameter.\r\nfail=c_long(0)\r\n<\/strong><\/pre>\n
e04abc(c_func,e1,e2,byref(a),byref(b),max_fun,byref(xt),byref(f),byref(comm),fail)\r\n<\/strong><\/pre>\n
print \"The minimum is in the interval  [%s,%s]\" %(a.value,b.value)\r\nprint \"The estimated result = %s\" % xt.value\r\nprint \"The function value at the minimum is %s \" % f.value\r\nprint \"%s evaluations were required\" % comm.nf\r\n<\/strong><\/pre>\n
If you run this code then you should get the following output<\/p>\n
The minimum is in the interval  [4.4934093959,4.49340951173]\r\nThe estimated result = 4.49340945381\r\nThe function value at the minimum is -0.217233628211\r\n10 evaluations were required<\/pre>\n
IFAIL oddness<\/strong><\/p>\n
My original version of this code used the trick of setting the NAG IFAIL parameter to c_int(0)<\/strong> to save me from worrying about the NagError structure (this is first explained in part 1<\/a>) but when this example was tested it didn’t work on 64bit machines although it worked just fine on 32bit machines.\u00a0 Changing it to c_long(0)<\/strong> (as I have done in this example)\u00a0 made it work on both 32 and 64bit machines.<\/p>\n
What I find odd is that in all my previous examples, c_int(0) worked just fine on both architectures.\u00a0 Eventually I will cover how to use the IFAIL parameter ‘properly’ but for now it is sufficient to be aware of this issue.<\/p>\n
Scary Structures<\/strong><\/p>\n
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<\/strong> structure correct by carefully working though the C specification of it in the nag_types.h<\/strong> file (included with the NAG C library installation) and transposing it to Python.<\/p>\n
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<\/strong>.  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.<\/p>\n
Callback functions with Communication<\/strong><\/p>\n
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<\/a>.\u00a0 This NAG routine expects the C prototype of your objective function to look like this:<\/p>\n
void funct(double xc, double *fc, Nag_Comm *comm)<\/strong><\/p>\n