# Category Archives: Code

## Optimisation problems: Choosing a method

In this morning’s Code Coffee meeting, we discussed optimisation methods (for doing things like least-squares minimisation, not for optimising code). This is an extremely common problem in astrophysics – you have some data, and you want to find the model that best fits it. In fact, this is a very common problem in general, and so there are lots of methods out there for handling optimisation. Of course, depending on the type of problem you have, some methods will be better than others. So how do you choose between them? And how do you implement the methods? I’ve summarised some of the discussion we had on these topics below. There’s also a great website with a detailed optimisation guide here.

N.B. I should note that most problems that astrophysicists are interested in involve the optimisation of real functions given data (so-called continuous optimisation). There are vast branches of computer science research concerned with combinatorial optimisation, which I’m not so interested in here.

### Choosing which method to use

As I’ve already mentioned, optimisation is a really common problem, so there are very many different optimisation algorithms, most of them with a number of different implementations. Over the course of the meeting, we came up with a basic step-by-step guide to help you choose which method to use:

1. Is the problem linear? Use a linear method, since these are guaranteed to converge to the right answer (within some tolerance). The Conjugate Gradient method is a good first port of call for general linear problems. If you have a particularly simple problem, you could even directly invert the matrix for your linear operator.
2. Can the problem be rewritten to be linear? If you can change variables or otherwise re-cast the problem in such a way that it becomes linear, do that and use a linear solver!
3. Is simple brute-force good enough? Some of the simpler (or at least better-known) optimisation methods are very robust, and some are guaranteed to converge, even for pretty general non-linear problems. If you have a small number of parameters, or not too many data points, it’s often quicker and easier to use one of these. Ryan Houghton recommended four, in order of increasing complexity: Amoeba; Levenberg-Marquardt; NNLS (non-negative least-squares); and BVLS (bound variable least-squares). The first two are discussed in Numerical Recipes, and the latter two have pretty widely-used implementations.
4. Can you analytically calculate gradients? If you can analytically/cheaply calculate the Jacobian and/or Hessian matrix of your problem (i.e. the derivatives of the function you’re minimising with respect to the model parameters), you can use a gradient-based method. Providing this extra information can result in substantial increases in the speed and convergence properties of the optimiser, so it’s worth doing, especially for larger problems.
5. Is your problem computationally intensive/complex? Some problems are just too large or too messy to use any of the methods advocated so far, so you’ll have to look at more sophisticated methods with more complex implementations. Fortunately, there are some nice, general libraries that implement lots of fancy algorithms (e.g. NLopt, or some of the codes here).
6. Do you have access to a supercomputer? If more sophisticated algorithms are still unsuitable, you might have to bite the bullet and just use a brute-force method and a lot of computer time. This can be nasty, especially if you don’t know how long it will take for the optimiser to converge. In this case, you might want to spend some more time thinking about clever ways to simplify your problem.

### Interesting implementations

Over the course of the discussion, quite a few nice-looking optimisation libraries were mentioned. I’ve listed a few of them below.

• NLopt has interfaces for lots of programming languages, and implements a wide range of optimisation algorithms. You can switch between algorithms by simply changing a single line of code.
• L-BFGS is a limited-memory code for doing constrained and unconstrained large-scale optimisation problems. Consider using this if your problem has lots of parameters, or you have quite a lot of data. There’s also KNITRO from the same author, which looks more advanced.
• levmar is a nice code for doing Levenberg-Marquardt optimisation, with support for a number of programming languages.
• NNLS is a non-negative least-squares implementation, and is available for a number of programming languages commonly used in astrophysics. It’s supposed to have pretty good convergence properties.
• More algorithms are listed here, and particular implementations of them are either listed here, or are only a Google search away.

### Properties of optimisation algorithms

• Convergence: Some algorithms are guaranteed to converge, but that doesn’t mean that they will converge to the correct answer!
• Constrained optimisation: Sometimes it’s necessary to put constraints into the optimisation problem (e.g. you don’t want the mass of your black hole to go negative). You might be able to avoid doing this if you have a good guess for the starting point of the optimisation, so it’s less likely to wander off into a disallowed region of parameter space. If not, there are a number of methods for performing constrained optimisations, but they tend to be reliable only for linear constraints. Word on the street is that codes capable of handling more complicated non-linear constraints are jealously-guarded Wall Street secrets…

Thanks go to Mark Walker, Joe Zuntz, Min-Su Shin, Sigurd Naess and Ryan Houghton for contributing greatly to this discussion.

## Code Coffee: One-line guides to getting started in astro computing

I was away on my travels last week, so I missed the first Code Coffee of the new term here in Oxford. Joe Zuntz was present to officiate, though, and led a discussion (allegedly) for the benefit of our new doctoral students, who’d just arrived at that point. The result was a collection of “one-line guides” on getting started with astrophysics-style computing and programming, which I’ve copied below. (Be warned: It’s a little Oxford-specific in parts.)

## Don’t code maths!

### Files

• Don’t invent your own file formats.
• FITS tables are not a great format. Consider using HDF instead.

### Compilers

• GNU: gcc, g++, gfortran (freely-available, cross-platform)
• Intel: icc, ifort (tend to be faster for some applications)
• Avoid mixing Intel/GNU compilers, and different versions of the compilers, when using Fortran. The compiled Fortran modules that are produced may not be compatible between different versions.

## Handling Fortran binary files with Python

Today’s top tip, straight from the front lines of, err, cosmology, concerns the handling of Fortran binary files with Python. If you have to deal with legacy Fortran code, you might find that it has a penchant for outputting files in an “unformatted” binary format. Exactly how the data is stored in this format can apparently depend on a lot of factors, such as the endian-ness of your machine. Hurrah! The files are non-portable. The necessity to figure what exactly is being dumped into the file can also cause headaches if you’re trying to read Fortran output into a program written in another language (for example, C/C++).

And it gets worse – what if you’re trying to write a file to be read in by the Fortran code, using a different language? Well, this particular problem cropped up for me, and of course I wanted to use Python to solve it.

For Python, fortunately, Neil Martinsen-Burrell has written FortranFile, which is subclassed off the standard Python file object, and which provides some convenient functions for reading and writing Fortran-formatted files. (There’s also an entry in the SciPy Cookbook that is relevant for Fortran I/O.) The pertinent function for me was `writeReals()`, which writes a NumPy array to a file in the appropriate format. The only subtlety I came across while briefly fiddling with this was that, when constructing a new `FortranFile` object for writing, you have to pass the `mode` keyword (e.g. `mode="w"`), which is inherited from file. (There’s a thread with a couple of tips on using FortranFile here.)

There are probably better ways of doing this than downloading someone else’s code off the internet, but like I said, I didn’t want to spend all afternoon on it.

## Drawing random numbers from a multivariate Gaussian distribution

I like software libraries. A software library means that someone has already done some of the heavy lifting for me. A software library means that I won’t have to reinvent the wheel. It means that I might not have to spend my afternoon fixing bugs and stringently validating my own code*.

It also means that I don’t necessarily have a good idea of what the computer is doing “behind the scenes”; that I’m potentially unaware of the clever (or not) algorithm that is being implemented, and what it’s limitations (or bugs) are.

For a long time now, I’ve been blithely accepting that numpy.random.multivariate_normal applies some sensible, but mysterious, procedure in order to generate random variables from a multivariate Gaussian (normal) distribution. And it does – I just never tried to unravel said mystery. So, prompted by Dag Sverre, and courtesy of Wikipedia,  I present a Python snippet that draws multivariate Gaussians given (a) some way of generating a set of standard univariate Gaussian random numbers, and (b) a Cholesky decomposition routine. In all, it’s about 10% faster than the NumPy routine, which applies a bunch of validity tests that I don’t.

`import numpy as np`

``` # Define mean, covariance matrix, and no. of samples required mean = ... # some array of length N cov = ... # some positive-definite NxN matrix Ndraws = 1000 # Do factorisation (can store this and use it again later) L = np.linalg.cholesky(cov) # Get 3*Ndraws Gaussian random variables (mean=0, variance=1) norm = np.random.normal(size=Ndraws*3).reshape(3, Ndraws) ```

```# Construct final set of random numbers (with correct mean) rand = mean + np.dot(L, norm) ```

Of course, this now introduces two new mysteries: How do you efficiently implement Cholesky decomposition, and how do you generate univariate Gaussian random numbers? Those are exercises for another day, I think.

* It depends on how sloppy I’m feeling; of course, you should validate everything that matters, regardless of who wrote it.

## Healpix map2gif options

I’m doing a lot of work with CMB maps at the moment, most of which are pixelised using Healpix. It’s good to be able to visualise the maps, and the Healpix distribution comes with a utility to do just that: `map2gif`. Unfortunately, map2gif has somewhat sparse documentation – there are a few options that aren’t documented and that sort of thing. I ended up poking through the source code to find some options, so I figured I should share them here to make other peoples’ lives a little easier.

Rough full-sky map of the Sunyaev-Zel’dovich effect, extracted from WMAP data.

(You can access much of this information by typing `map2gif -hlp` to print usage instructions, but they’re quite basic.)

#### Basic usage

At its simplest, you need to supply filenames for an input FITS file (in the Healpix format) and an output GIF:

`\$ map2gif -inp input.fits -out map.gif`

If you have a CMB map, you’re probably most interested in seeing the CMB anisotropies. Unfortunately, residual foregrounds can mess up the scaling of the map, so you have to add a scale yourself. For a map in units of micro-Kelvin, the following scale should produce good results for the anisotropies:

`\$ map2gif -inp input.fits -out map.gif -min -300 -max 300`

For most uses, you probably want a colour scale to be displayed on the output image too. Use the `-bar` option for this (it’s a boolean value, so specify True):

`\$ map2gif -inp input.fits -out map.gif -bar True`

A rough (Gnomonic-projected) map of the Sunyaev-Zel’dovich effect centred around the Virgo cluster, extracted from WMAP data.

#### Projection and centering

For a lot of applications, you’ll be wanting to zoom in on some region of the map rather than seeing the full sky. Healpix supports the Mollweide (default) and Gnomonic projections, which you can switch between using the `-proj` switch. You only need to type the first three letters of each projection. If you want to zoom in on some specific area, use Gnomonic and specify a different map origin using the (FITS file coordinate system-dependent) `-lat` and `-lon` switches:

`\$ map2gif -inp input.fits -out map.gif -proj gno -lat 45 -lon 87`

By default this will look quite zoomed-in, so you might want to use the resolution (`-res`) switch to zoom out a little (higher numbers zoom out more):

`\$ map2gif -inp inp.fits -out map.gif -proj gno -lat 45 -lon 87 -res 4`

#### Presentation (colours, labels, image size)

There are a number of switches for controlling the presentation of the output image, most of which are explained sufficiently well in the Healpix documentation. For example, you can add a title to the plot (`-ttl`), increase its width in pixels (`-xsz`; the aspect ratio is kept the same) or use a different colormap (`-col`; try different numbers to get different colour maps).

#### Other options

There are a few other options that are a bit more specialised; for example, you can scale the signal by arcsinh (`-ash`) for some reason. A particularly interesting one for CMB applications is the map selection option (`-sig`); it’s typical to include maps for several Stokes parameters (I, Q, and U) inside a single FITS file for experiments that do polarisation mapping, so you can use this switch to select between them. For example, to pick the Q map out of a WMAP IQU map, use:

`\$ map2gif -inp input.fits -out map.gif -sig 2`

There are also some basic rescaling operations that you can apply using the addition (`-add`) and multiplication (`-mul`) switches. Read more about those in the Healpix docs.

#### Listing of command line options

Here’s a full listing of command-line options to map2gif, taken directly from the output of the `-hlp` switch:

```usage: MAP2GIF [-inp input_file] [-out output_file] [-max usermax] [-min usermin] [-add offset] [-mul factor] [-pro projection] [-res gridresn] [-lon lon0] [-lat lat0] [-col color_table] [-sig signal] [-xsz xsize] [-log logflg] [-ash ashflg] [-bar add color bar] [-ttl title] [-hlp]```

## Weird gfortran compiler error of the week

In response to the following compiler error from gfortran, the correct response is apparently to run `make clean` and then `make` again:

`Error: More actual than formal arguments in procedure call at (1)`

God knows what was actually going wrong.

## Getting codes in different languages to interact

In a Code Coffee meeting last month (9th July), we discussed ways of getting codes in different languages to talk to one another. This is quite a common task; for example, if you need your modern data analysis pipeline to interact with some legacy code (normally Fortran!), or if you have a mostly Python code that has a computationally-intensive bit that would be faster written in C. Our focus was on Python, C, and Fortran, since they seem to be the most commonly used languages in the group.

People brought a number of different solutions with them. Some were very easy to implement, just as a simple few-line block in your Python code that could be used to speed up a loop. Others were quite a bit more complicated, requiring extra work on building modules and getting your hands dirty with memory management, but having the virtue of being much more flexible. The methods we discussed are summarised below, in increasing order of complexity. If you skip to the bottom, there’s also a link to a tarball, with some examples.

(Edit: Be sure to check out some of the interesting comments, below, too. They mostly concern Cython.)

## Weave

Weave is a way of writing C code “inline” – directly into a Python script.  Python variables are automatically made available to the C code (you sometimes need to do a bit of extra work to cast them to a specific type, but not much), which makes it nice and easy to use, with little boilerplate required in the C code. Essentially, you can just stick a fast, few-line C loop in the middle of a Python function, then call Weave and it’ll handle compilation and shuffling data between Python and C for you automagically, without requiring any complicated wrapper code or makefile magic. There’s also a “blitz” mode, which takes a Python/NumPy expression and converts it into C++ in the background, which can also help speed things up (although it’s generally slower than the inline mode).

Unfortunately, the documentation isn’t great, but there are plenty of examples out there (SageMath have a couple of very clear ones). It should definitely be your first port of call if all you want to do is speed up part of your Python code. It’s not so good if you want to interface with legacy code, or have something more complicated that you need to do in C (e.g. something that’s split into different functions). See also Cython, which seems to be a bit more flexible, but not too much more difficult to use. There’s a nice speed comparison with NumPy, Weave blitz and inline, and MATLAB here, and a Cython/Weave example (with benchmarks) here. There’s also an interesting blog post on getting GSL to work with Weave/Cython.

## f2py

f2py is a way of automatically generating Python modules from Fortran code. It seems pretty easy to use – all you need to do is run it over your Fortran source code, and it will produce a compiled module that you can simply import into your Python code and use as you would any other. The f2py program scans the Fortran code and figures out how to deal with function arguments and convert between types itself, which is very handy. It even generates function documentation (docstrings) for the Python module, telling you what format the function arguments need to be in. One problem with f2py, pointed out by Ryan, is that it tends to fail silently, making it difficult to debug! It’s apparently quite robust though, so this probably won’t be necessary that often. The documentation looks pretty good, and you can find other resources here. (N.B. They ask you to cite their paper if you use f2py in research code.)

## Python ctypes

ctypes is a Python module for handling calls to shared libraries (whether they are written in C/C++ or Fortran). The syntax for loading libraries and calling functions is exceptionally simple – little different from importing a native Python module and calling functions in that. However, there’s some added complexity that comes as a result of having to convert between C and Python data types in function arguments. For the most part, all you need to do is call a simple ctypes helper function to make the conversion; e.g. `c_double(42.3)` will give you something that can be passed to a function that expects a C double as an argument. In more recent versions of NumPy, ndarrays also have ctypes conversion methods (see the `ndarray.ctypes.data_as()` method; you probably want to get a pointer to the array, for which you should use `ctypes.POINTER`), although there are also some subtleties to be aware of in terms of how Python/NumPy and C/C++ or Fortran store multidimensional arrays in memory (see the order property of NumPy arrays, for example). The Sage website has a good tutorial on working with ndarrays in ctypes. A slightly nicer way of handling types conversion is to write a small Python wrapper for the C function (or whatever), and use the argtypes attribute to specify what arguments the function accepts. A comprehensive tutorial can be found here, which explains most of the workings of ctypes, and which has good examples.

The tutorial on the Python website mostly concerns calling existing system libraries, but of course you can write your own libraries too. For C/C++ code compiled with GCC, this is pretty simple – all you need to do is add a couple of compiler flags (see Section 3.4 of the TDLP Library Program tutorial). Of course, your C code should be set out like a library first – just compiling something that’s meant to be a standalone program as a library won’t work too well, and you’ll need to write header files etc. That’s easy to do, and existing “non-library” code can often be refactored into a working library with very little effort. There’s a very nice, simple tutorial on writing and compiling libraries here. Note that the libraries don’t need any “special treatment” to interface with Python using ctypes – they’re just standard C/Fortran libraries, and can happily be called by other C/Fortran programs too.

## Python C API

Big chunks of Python are written in C, so it should come as no surprise that there is a C library which provides access to Python data types and functions. Access it by including Python.h in your C code. I’ll keep discussion of this one to a deliberately high-level, because it’s significantly more complicated than the others, and doesn’t seem to be fantastically well documented (although see the official Python docs for it here and here). This method allows you to write new Python modules directly in C; the end result will be something that looks and behaves exactly like a standard Python module (of the sort that you may have written before, in Python), which requires no fiddling about with data types or what have you in the Python code that uses it. That’s because all of the fiddling around is done in the C code! I find this quite a bit trickier to write – your C files need some boilerplate to get themselves (and their constituent functions) recognised by the Python interpreter, and there’s a bit of makefile magic required to actually build and install the module too. Plus, you have to figure out how to handle type conversions between native C types and Python types (mostly PyObjects). It’s not super-difficult, but it is fiddly in places. Take a look at this nice tutorial by Dan Foreman-Mackey (and another one from JPL which looks more specifically at handling NumPy arrays).

Some of the more confusing issues that you’re likely to run into are to do with reference counting. This is a reasonably technical computer science concept, which has to do with how Python stores and labels variables. (Recall that, in Python, the same data can be referred to by multiple variable names.) It’s important to make sure your C code properly tracks references to Python objects as it manipulates them; otherwise, you’re likely to run into a host of problems, from memory leaks to weird crashes. If you’re going to use the Python C API, I strongly recommend that you invest a good amount of time in understanding how this works.

All in all, I think the C API is only the way forward if you’re specifically setting out to write a new Python module for use elsewhere; if all you have is some existing C code that you want to quickly plug in to some Python code, it’s going to be a lot of hassle. Still, once you’ve set up a project once, and used the C API a little bit, it’s a lot quicker to get up and running with the next project.

## Summary

All in all, ctypes is probably the best route to go for most people, unless you have something particularly trivial you want to speed-up. It’s fast (native C function calls), easy enough to set up, and conceptually simple to understand. It’s not quite as “automatic” as Weave/Cython or f2py, but you’ll probably end up with cleaner, more flexible, and more robust code for your troubles. And it’s much easier to work with than the Python C API.

The Performance Python page has a nice comparison of some of the methods mentioned above. A tarball containing example code for most of the above methods can be downloaded from here; thanks to Ryan Houghton, Neale Gibson and Joe Zuntz for providing these.