Walking Randomly

Late last year I was asked to give a talk at my University to a group of academics studying financial mathematics (Black-Scholes, monte-carlo simulations, time series analysis – that sort of thing).  The talk was about the NAG Numerical Library, what it could do, what environments you can use it from, licensing and how they could get their hands on it via our site license.

I also spent a few minutes talking about Python including why I thought it was a good language for numerical computing and how to interface it with the NAG C library using my tutorials and PyNAG code.  I finished off with a simple graphical demonstration that minimised the Rosenbrock function using Python, Matplotlib and the NAG C-library running on Ubuntu Linux.

“So!”, I asked, “Any questions?”

The very first reply was “Does your Python-NAG interface work on Windows machines?” to which the answer at the time was “No!”  I took the opportunity to ask the audience how many of them did their numerical computing in Windows (Most of the room of around 50+ people), how many people did it using Mac OS X (A small group at the back), and how many people did it in Linux (about 3).

So, if I wanted the majority of that audience to use PyNAG then I had to get it working on Windows.  Fortunately, thanks to the portability of Python and the consistency of the NAG library across platforms, this proved to be rather easy to do and the result is now available for download on the main PyNAG webpage.

Let’s look at the main differences between the Linux and Windows versions

Loading the NAG Library

In the examples given in my Linux NAG-Python tutorials, the cllux08dgl NAG C-library was loaded using the line

libnag = cdll.LoadLibrary("/opt/NAG/cllux08dgl/lib/libnagc_nag.so.8")

In Windows the call is slightly different. Here it is for CLDLL084ZL (Mark 8 of the C library for Windows)

C_libnag = ctypes.windll.CLDLL084Z_mkl

If you are modifying any of the codes shown in my NAG-Python tutorials then you’ll need to change this line yourself. If you are using PyNAG, however, then it will already be done for you.

Callback functions

For my full introduction to NAG callback functions on Linux click here. The main difference between using callbacks on Windows compared to Linux is that on Linux we have

C_FUNC = CFUNCTYPE(c_double,c_double )

but on Windows we have

C_FUNC = WINFUNCTYPE(c_double,c_double )

There are several callback examples in the PyNAG distribution.

That’s pretty much it I think.  Working through all of the PyNAG examples and making sure that they ran on Windows uncovered one or two bugs in my codes that didn’t affect Linux for one reason or another so it was a useful exercise all in all.

So, now you head over to the main PyNAG page and download the Windows version of my Python/NAG interface which includes a set of example codes.  I also took the opportunity to throw in a couple of extra features and so upgraded PyNAG to version 0.16, check out the readme for more details.  Thanks to several employees at NAG for all of their help with this including Matt, John, David, Marcin and Sorin.

Over at Mathrecreation, Dan Mackinnon has been plotting two dimensional polygonal numbers on quadratic number spirals using a piece of software called Fathom. Today’s train-journey home project for me was to reproduce this work using Mathematica and Python.

I’ll leave the mathematics to Mathrecreation and Numberspiral.com and just get straight to the code. Here’s the Mathematica version:

• Download the Mathematica Notebook
• Download a fully interactive version that will run on the free Mathematica Player

p[k_, n_] := ((k - 2) n (n + 1))/2 - (k - 3) n;
Manipulate[
 nums = Table[p[k, n], {n, 1, numpoints}];
 If[showpoints,
  ListPolarPlot[Table[{2*Pi*Sqrt[n], Sqrt[n]}, {n, nums}]
   , Joined -> joinpoints, Axes -> False,
   PlotMarkers -> {Automatic, Small}, ImageSize -> {400, 400}],
  ListPolarPlot[Table[{2*Pi*Sqrt[n], Sqrt[n]}, {n, nums}]
   , Joined -> joinpoints, Axes -> False, ImageSize -> {400, 400}
   ]
  ]
 , {k, 2, 50, Appearance -> "Labeled"}
 , {{numpoints, 100, "Number of points"}, {100, 250, 500}}
 , {{showpoints, True, "Show points"}, {True, False}}
 , {{joinpoints, True, "Join points"}, {True, False}}
 , SaveDefinitions :> True
 ]

A set of 4 screenshots is shown below (click the image for a higher resolution version). The formula for polygonal numbers, p(kn,n), is only really valid for integer k and so I made a mistake when I allowed k to be any real value in the code above. However, I liked some of the resulting spirals so much that I decided to leave the mistake in.

When Dan showed the images he produced from Fathom, he mentioned that he would one day like to implement these systems in a more open platform. Well, you can’t get much more open than Python so here’s a piece of code that makes use of Python and the matplotlib package (download the file – polygonal_spiral.py).

#!/usr/bin/env python
import matplotlib.pyplot as plt
from math import *
 
def p(k,n):
   return(((k-2)*n*(n+1))/2 -(k-3)*n)
 
k=13
 
polygonal_nums = [p(k,n) for n in range(100)]
theta = [2*pi*sqrt(n) for n in polygonal_nums]
r = [sqrt(n) for n in polygonal_nums]
 
myplot = plt.polar(theta,r,'.-')
plt.gca().axis('off')
plt.show()

Now, something that I haven’t figured out yet is why the above python code works just fine for 100 points but if you increase it to 101 then you just get a single point. Does anyone have any ideas?

Unlike my Mathematica version, the above Python program does not give you an interactive GUI but I’ll leave that as an exercise for the reader.

If you enjoyed this post then you might also like:

• Simulating Harmonographs
• Quadraflakes,Pentaflakes,Hexaflakes and more

One of my plans for next year is to start giving short talks and tutorials about scientific Python to various groups around the University of Manchester.  Ever since I attended EuroSciPy in Leipzig earlier this year I have been excited about the possibilities offered by Python for students, researchers and teachers.  I genuinely believe that Python and SAGE will do for MATLAB and Mathematica users what R has done for users of SPSS and Stata.  It’s the future…I’ve tasted it!  Here are a few links to others who agree with me.

• Scientific Python tools use rising in education – a blog post from Fernando Perez
• A hands-on 2-day workshop at UC Berkeley – includes videotaped lectures
• Python for High Performance Computing (HPC)
• Should I switch to Python – Are you a MATLAB user considering the switch?
• Python in the Scientific World – From Guido van Rossum – Pythons benevolent dictator

Earlier this week I mentioned my plans to a colleague of mine over coffee and he said “Cool, I didn’t realise that you were a Python expert!”

Of course I’m NOT an expert on Python but it turns out that you just don’t have to be an expert in Python in order to get useful stuff done and THAT is my point!

Earlier today I was chatting to a lecturer over coffee about various mathematical packages that he might use for an upcoming Masters course  (note – offer me food or drink and I’m happy to talk about pretty much anything). He was mainly interested in Mathematica and so we spent most of our time discussing that but it is part of my job to make sure that he considers all of the alternatives – both commercial and open source. The course he was planning on running (which I’ll keep to myself out of respect for his confidentiality) was definitely a good fit for Mathematica but I felt that SAGE might suite him nicely as well.

“Does it have nice, interactive functionality like Mathematica’s Manipulate function?” he asked

Oh yes! Here is a toy example that I coded up in about the same amount of time that it took to write the introductory paragraph above (but hopefully it has no mistakes). With just a bit of effort pretty much anyone can make fully interactive mathematical demonstrations using completely free software. For more examples of SAGE’s interactive functionality check out their wiki.

Here’s the code:

def ftermSquare(n):
 return(1/n*sin(n*x*pi/3))

def ftermSawtooth(n):
 return(1/n*sin(n*x*pi/3))

def ftermParabola(n):
 return((-1)^n/n^2 * cos(n*x))

def fseriesSquare(n):
 return(4/pi*sum(ftermSquare(i) for i in range (1,2*n,2)))

def fseriesSawtooth(n):
 return(1/2-1/pi*sum(ftermSawtooth(i) for i in range (1,n)))

def fseriesParabola(n):
 return(pi^2/3 + 4*sum(ftermParabola(i) for i in range(1,n)))

@interact
def plotFourier(n=slider(1, 30,1,10,'Number of terms')
,plotpoints=('Value of plot_points',[100,500,1000]),Function=['Saw Tooth','Square Wave','Periodic Parabola']):
    if Function=='Saw Tooth':
     show(plot(fseriesSawtooth(n),x,-6,6,plot_points=plotpoints))
    if Function=='Square Wave':
     show(plot(fseriesSquare(n),x,-6,6,plot_points=plotpoints))
    if Function=='Periodic Parabola':
     show(plot(fseriesParabola(n),x,-6,6,plot_points=plotpoints))

Just a quick note to say that I’ve released version 0.15 of PyNAG – my Python interface to the NAG C library and given it its own dedicated page on Walking Randomly.  New features include

• Basic support for 165 functions from Mark 8 of the NAG C library.
• 29 simple but fully functional examples.
• Added initial support for NAG’s Complex structure.

If you have learned how to program then the odds are that you have written a program that displays a list of prime numbers.  A prime number generator is a programming rite of passage that is almost up there with the ubiquitous ‘Hello World‘ application (for programmers who are going on to do numerical work at least).

There are various ways you might choose to do it but I bet that your first (or even your second or third) thought wouldn’t be ‘I know – I’ll use a regular expression‘. Someone known only as ‘Abigail’ obviously thought that the world desperately needed such a regex since he / she came up with the following incantation

^1?$|^(11+?)\1+$

Which is a regular expression that I assure you can be used to generate a list of all the prime numbers.  If you are in the mood for a puzzle then sit back, grab a cup of coffee and try to work out HOW it can generate the primes.

If you give up then head over to Python Programming for a full explanation along with Python code that demonstrates that it really does work.

Can anyone come up with regular expressions for other interesting series of numbers?

Update (27th August 2009): Someone pointed me to an alternative explanation here.

I’ll be attending the 2009 EuroSciPy scientific Python conference in Leipzig next month – July 25th and 26th to be precise – and would love to hear from anyone else who’ll be attending.

Back in 1997 I was a 2nd year undergraduate of Physics and I was taught how to program in Fortran, a language that has survived over 40 years due to several facts including

• It is very good at what it does.  Well written Fortran code, pushed through the right compiler is screamingly fast.
• There are millions of lines of legacy code still being used in the wild.  If you end up doing research in subjects such as Chemistry, Physics or Engineering then you will almost certainly bump into Fortran code (I did!).
• A beginner’s course in Fortran has been part of the staple diet in degrees in Physics, Chemistry and various engineering disciplines (among others) for decades.
• It constantly re-invents itself to include new features.  I was taught Fortran 77 (despite it being 1997) but you can now also have your pick of Fortran 90, 95, 2003, and soon 2008.

Almost everyone I knew hated that 1997 Fortran course and the reasons for the hatred essentially boiled down to one of two points depending on your past experience.

• Fortran was too hard!  So much work for such small gains (First time – programmers)
• The course was far too easy.  It was just a matter of learning Fortran syntax and blitzing through the exercises. (People with prior experience)

The course was followed by a numerical methods course which culminated in a set of projects that had to be solved in Fortran.  People hated the follow on course for one of two reasons

• They didn’t have a clue what was going on in the first course and now they were completely lost.
• The problems given were very dull and could be solved too easily.  In Excel!  Fortran was then used to pass the course.

Do you see a pattern here?

Fast forward to 2009 and I see that Fortran is still being taught to many undergraduates all over the world as their first ever introduction to programming.  Bear in mind that these students are used to being able to get interactive 3D plots from the likes of Mathematica, Maple or MATLAB and can solve complex differential equations simply by typing them into Wolfram Alpha on the web.  They can solve problems infinitely more complicated then the ones I was faced with in even my most advanced Fortran courses with just a couple of lines of code.

Learn Fortran – spend a semester achieving not very much

These students study subjects such as physics and chemistry because they are interested in the subject matter and computers are just a way of crunching through numbers (and the algebra for that matter) as far as they are concerned.  Despite having access to untold amounts of computational and visualisation power coupled with easy programming languages thanks to languages such as MATLAB and Python, these enthusiastic, young potential programmers get forced to bend their minds around the foibles of Fortran.

For many it’s their first ever introduction to programming and they get forced to work with one of the most painful programming languages in existence (in my opinion at least).   Looking at a typical one-semester course it seems that by the time they have finished they will be able to produce command line only programs that do things like multiply matrices together (slowly), solve the quadratic equation and find the mean and standard deviation of a list of numbers.

Not particularity impressive for an introduction to the power of computation in their subject is it?

Fortran syntax is rather unforgiving compared to something like Python and it takes many lines of code to achieve even relatively simple results.  Don’t believe me?  OK, write a program in pure Fortran that gives a plot of a Sin(n*x) for integer n and x ranging from -2*pi to 2*pi.  Now connect that plot up to a slider control which will control the value of n.  Done?  Ok – now get it working on another operating system (eg if you originally used Windows, get it to work on Mac OS X).

Now try the same exercise in Python or Mathematica.

This is still a very basic program but it would give a much greater sense of achievement compared to finding the mean of list of numbers and could easily be extended for more able students (Fourier Series perhaps).  I believe that many introductory Fortran programming courses end up teaching students that programming means ‘Calculating things the hard way’ when they should be leaving an introductory course with the opposite impression.

Learn Fortran – and never use it again.

Did you learn Fortran at University?  Are you still using it?  If you answer yes to both of those questions then there is a high probability that you are still involved in research or that advanced numerical analysis is the mainstay of your job.  I know a lot of people in the (non-academic) IT industry – many of them were undergraduates in Physics or Chemistry and so they learned Fortran.  They don’t use it anymore.  In fact they never used it since completing their 1st semester, 2nd year exam and that includes the computational projects they chose to do as part of their degrees.

The fact of the matter is that most undergraduates in subjects such as Physics end up in careers that have nothing to do with their degree subject and so most of them will never use Fortran ever again.  The ‘programming concepts’ they learned might be useful if they end up learning Java, Python or something along those lines but that’s about it.  Would it not be more sensible to teach a language that can support the computational concepts required in the underlying subject that also has an outside chance of being used outside of academia?

Teach Fortran – and spend a fortune on compilers

There are free Fortran compilers available but in my experience these are not the ones that get used the most for teaching or research and this is because the commercial Fortran compilers tend to be (or at least,they are perceived to be) much better.  The problem is that when you are involved with looking after the software portfolio of a large University (and I am) then no one will agree on what ‘the best’ compiler is.  (Very) roughly paraphrased, here are some comments I have received from Fortran programmers and teachers over the last four years or so.

• When you teach Fortran, you MUST use the NAG Fortran compiler since it is the most standards compliant.  The site license allows students to have it on their own machines which is useful.
• When you teach Fortran, you MUST use the Silverfrost compiler because it supports Windows GUI programming via Clearwin.
• We MUST avoid the Silverfrost compiler for teaching because it is not available on platforms such as Mac and Linux.
• We MUST avoid the NAG compiler for teaching because although it has a great user interface for Windows, it is command line only for Linux and Mac.  This confuses students.
• If we are going to teach Fortran then we simply MUST use the Intel Fortran Compiler.  It’s the best and the fastest.
• We need the Intel Fortran Compiler for research because it is the only one compatible with Abaqus.
• A site license for Absoft Fortran is essential.  It’s the only one that will compile <insert application here> which is essential for my group’s work.
• We should only ever use free Fortran compilers – relying on commercial offerings is wrong.

I’m only getting started!  Choose a compiler, any compiler – I’ll get back to you within the week with an advocate who thinks we should get a site license for it and another who hates it with a passion.  Let’s say you had a blank cheque book and you gave everyone exactly what they wanted – your institution would be spending tens of thousands of pounds on Fortran compilers and you’d still be annoying the free-software advocate.

Needless to say, most people don’t have a blank cheque book.  At Manchester University (my workplace) we support 2 site-wide Fortran compilers

• The NAG Compiler – recommended for teaching.
• The Intel Fortran Compiler – recommended for research

We have a truly unlimited site license for NAG (every student can have a copy on their own machine if they wish) which makes it perfect from a licensing point of view and many people like to use it to teach.  The Widows version has a nice GUI and help system for example – perfect for beginners.  There are Mac and Linux versions too and although these are command-line only, they are better than nothing.

The Intel Fortran Compiler licenses we have are in the form of network licenses and we have a very limited number of these – enough to support the research of every (Windows and Linux) Fortran programmer on campus but nowhere near enough for teaching.  To get enough for teaching would cost a lot.

There are also pockets of usage of various other compilers but nothing on a site-wide basis.

Of course, not everyone is happy with this set-up, and I was only recently on the receiving end of a very nasty email from someone because I had to inform him that we didn’t have the money to buy his compiler of choice.  It’s not my fault that he couldn’t have it (I don’t have a budget and never have had) but he felt that I was due an ear lashing I guess and blamed all of his entire department’s woes on me personally.  Anway, I digress….

The upshot is that Fortran compilers are expensive and no one can agree on which one(s) you should get.  If you get them all to please (almost) everybody then you will be very cash poor and it won’t be long before the C++ guys start knocking on your door to discuss the issue of a commercial C++ compiler or three….

Learn Fortran – because it is still very useful

Fortran has been around in one form or another for significantly longer than I have been alive for a very good reason. It’s good at what it does. I know a few High Performance Computing specialists and have been reliably informed that if every second counts in your code execution – if it absolutely, positively, definitely has to run as fast as humanly possible then it is hard to beat Fortran. I take their word for it because they know what they are talking about.

Organisations such as the Numerical Algorithm’s Group (NAG) seem to agree with this stance since their core product is a Fortran Library and NAG have a reputation that is hard to ignore.  When it comes to numerics – they know their stuff and they work in Fortran so it MUST be a good choice for certain applications.

If you find yourself using research-level applications such as Abaqus or Gaussian then you’ll probably end up needing to use Fortran, just as you will if you end up having to modify one of the thousands of legacy applications out there.  Fortran is a fact of life for many graduate students.  It was a fact of life for me too once and I hated it.

I know that it is used a lot at Manchester for research because ,as mentioned earlier, we have network licenses for Intel Fortran and when I am particularly bored I grep the usage logs to see how many unique user names have used it to compile something.  There are many.

I am not for a second suggesting that Fortran is irrelevant because it so obviously isn’t.  It is heavily used in certain, specialist areas.  What I am suggesting is that it is not the ideal language to use as a first exposure to programming.  All of the good reasons for using Fortran seem to come up at a time in your career when you are a reasonably advanced programmer not when you first come across the concepts loops and arrays.

Put bluntly I feel that the correct place for learning Fortran is in grad school and only if it is needed.

If not Introductory Fortran then Introductory what?

If you have got this far then you can probably guess what I am going to suggest – Python.  Python is infinitely more suitable for beginner programmers in my opinion and with even just a smattering of the language it is possible to achieve a great deal.  Standard python modules such as matplotlib and numpy help take care of plotting and heavy-duty numerics with ease and there are no licensing problems to speak of .  It’s free!

You can use it interactively which helps with learning the basics and even a beginner can produce impressive looking results with relatively little effort.  Not only is it more fun than Fortran but it is probably a lot more useful for an undergraduate because they might actually use it.  Its use in the computational sciences is growing very quickly and it also has a lot of applications in areas such as the web, games and general task automation.

If and when a student is ready to move onto graduate level problems then he/she may well find that Python still fulfills their needs but if they really do need raw speed then they can either hook into existing Fortran libraries using Python modules such as ctypes and f2py or they can roll up their sleeves and learn Fortran syntax.

Using Fortran to teach raw beginners is a bit like teaching complex numbers to kindergarten children before they can count to ten.  Sure, some of them will need complex numbers one day but if it was complex numbers or nothing then a lot of people would really struggle to count.

So, I am (finally) done.  Do I have it all wrong?  Is Fortran so essential that Universities would be short changing students of numerate degrees if they didn’t teach it or are people just teaching what they were taught themselves because that’s the way it has always been done?  I am particularly interested in hearing from either students or teachers of Fortran but, as always, comments from everyone are welcome – even if you disagree with me.

If you are happy to talk in public then please use the comments section but if you would prefer a private discussion then feel free to email me.

(I am skating closer to my job than I usually do on this blog so here is the hopefully unnecessary disclaimer.  These are my opinions alone and do not necessarily reflect the policy of my employer, The Univeristy of Manchester.  If you find yourself coming to me for Fortran compiler support there then I will do the very best I can for you – no matter what my personal programming prejudices may be.  However, if you are up for a friendly chat over coffee concerning this then feel free to drop me a line. )

This is part 6 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.

Over the last few articles I have demonstrated that it is relatively straightforward to interface Python to the NAG C library.  The resulting code, however, looks a bit ugly and overly verbose since you have to include a lot of boiler plate code such as telling Python where the C library is, defining pythonic versions of C structures, defining various constants and so on.  I found myself doing a lot of copying and pasting in the early stages of working all this stuff out.

Just recently I made a small step forward by moving all of this boilerplate into a Python package which I have called PyNAG.   It’s a relatively simple package since all it does is define Pythonic versions of various C-structures, link to the NAG library, perform restypes on the function names, define a few often used NAG variables such as NAGERR_DEFAULT and so on.

The practical upshot is that you need to write a lot less code when interfacing with the C library.  For example here is the Python code to minimise a 2D function via the simplex method using e04ccc of the NAG library.

#!/usr/bin/env python
#Mike Croucher - April 2008
#A minimisation problem solved using e04ccc
from PyNAG import e04ccc,NAGERR_DEFAULT,E04_DEFAULT,Nag_Comm
from ctypes import *
from math import *

#Create the communication structure
comm = Nag_Comm() 

#set up the objective function
def py_obj_func(n,xc,objf,comm):
	objf[0] = exp(xc[0]) * (xc[0] * 4.0 * ( xc[0]+xc[1] ) + xc[1] * 2.0 * (xc[1] +1.0) +1.0 )
C_FUNC = CFUNCTYPE(None,c_int,POINTER(c_double),POINTER(c_double),POINTER(Nag_Comm) )
c_obj_func = C_FUNC(py_obj_func)

#Set up the starting point
n=2
x=(c_double*2)()
x[0]=c_double(0.4)
x[1]=c_double(-0.8)

#Set up other variables
objf = c_double()

e04ccc(n, c_obj_func, byref(x), byref(objf),E04_DEFAULT,comm,NAGERR_DEFAULT)

It’s still not as pretty as it would it if it were written in pure python but it’s a (small) step in the right direction I think.

The good

• PyNAG is free and always will be.  You’ll obviously need to buy a copy of the NAG C library though.
• It includes a definition of one of the largest C-Structures I have ever seen – Nag_E04_Opt – which has over 130 members.  This took me quite a while to get working and now you don’t need to go through the same pain I did.  All the C structures mentioned in these tutorials have been included to save your typing fingers.
• PyNAG has been tested on 32bit and 64bit flavours of Linux and seems to work fine.
• It is relatively easy to add definitions for NAG functions that are not already currently included.  In the best case scenario all you need to do is add two lines of code per NAG function. (They’re not all this easy though!)
• It comes with a set of fully functional examples (15 so far with more on the way) along with expected output.  This output has been checked against the original C versions of the code to ensure that the results are correct

The bad

• It only works on Linux at the moment.  Depending on interest levels (mine and yours) I may get it working on Mac OS X and Windows.
• It only includes all the boilerplate needed to interface with 22 NAG functions at the time of writing.  This is a very small percentage of the library as a whole.
• It doesn’t include all of the structures you might need to use throughout the entire library at the moment.
• It doesn’t include any of the enumerations that the NAG C library uses.   I plan to include some of the most commonly used ones but will simply have to define them as variables since ctypes doesn’t support enumerations as far as I know.
• This is my first ever Python package so things may be a bit rough around the edges. I welcome advice from more experienced Pythonistas concerning things like package deployment.
• The Python version of the NagError structure has a member called pprint rather than the original print. This is because print is a reserved keyword in Python so I couldn’t use it.
• It has only been tested on Python versions 2.5.x and 2.6.x.  I have no idea if it will work on other versions at the moment.
• Documentation is rather sparse at the moment.  I haven’t had time to do a better job.

TheUgly

• All PyNAG really does is save you from having to type in a load of boilerplate.  It doesn’t properly ‘pythonise’ the NAG C library and so it is painfully obvious that your code is talking to a C library and not a Python library.  Ideally I would like to change this and make the user interface a lot more user-friendly.  Any input on how best to go about this would be gratefully received.

The Installation

The first public release of PyNAG is version 0.12 and is currently available here (Linux only).  Installation instructions are

• Download the tar gzip file above and run the following commands
• tar -xvzf ./PyNAG-0.12.tar.gz
  cd PyNAG-0.12/src
• edit the following line in the file __init__.py so that it points to your installation of the library
  C_libnag= ctypes.cdll.LoadLibrary(“/opt/NAG/cllux08dgl/lib/libnagc_nag.so.8″)
• Move back into the root PyNAG-0.12 directory:
  cd..
• run the following command with administrator priviledges
  python setup.py install

The installer will scatter the PyNAG files where ever your machine is set up to put them.  On my machine it was put in the following location but your mileage may vary

/usr/local/lib/python2.6/dist-packages/PyNAG

The Examples

I have written a set of example scripts that show PyNAG in use and these can be found in the src/examples directory of the PyNAG distribution.   Mathematically speaking they are all rather dull but they do show how you can call the routines I have included so far.  More interesting examples will be appearing soon.

The Future

Getting the NAG C library to work with Python has always been a side project for me and I never expected to take it even this far.  Most of this work has been done on the train while commuting to and from work and, although it is far from perfect, I am happy with it so far.  There is much to do though and the rate it gets done depends on user-interest.  Things that spring to mind (in no particular order) are

• Get PyNAG on a proper CVS system such as github so that other people can easily contribute.
• Add support for Numpy matrices.
• Add support for more of the most commonly used NAG functions.
• Get a windows version up and running.
• Write some more interesting demos.
• Write some sort of automated test script that runs all the examples and compares the results with a set of expected results.
• See if I get any users and take requests

Final points.

I do not work for NAG – I just like their products. My professional involvement with them is via the University of Manchester where I am the site-license administrator for NAG products (among many other things).  I have a vested interest in getting the best out of their products as well as getting the best out of the likes of Wolfram, Mathworks and a shed-load of other software vendors.

One of my hobbies is pointing out to NAG where their library could do with extra functionality.  I am aided in this venture by a load of academics who use the library in it’s various guises and who contact me to say ‘I like NAG but it would be much better if it had the ability to do XYZ’.   I am very grateful to NAG for accepting my steady stream of emails with good grace and also for implementing some of our suggestions in the latest version of their library.  Feel free to point out any of your favourite algorithms that are missing from the library and why it should have them.

Finally, I’d like to thank the various people at NAG who helped me get my head around their library and for helping me test my initial tentative steps into the NAG/Python world.  I owe many of you several beverages of your choice.

As always, comments are welcomed.

Update (5th November 2009): PyNAG has been updated several times since this post was written. The latest version will always be available from http://www.walkingrandomly.com/?page_id=1662

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.