I recently found a solution to a problem that had been vexing me for about a year. In order to successfully import anything from scipy.integrate, I had edit the file scipy/integrate/__init__.py and comment out the line

from odepack import *

If not, I would get various import errors such as

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python2.5/site-packages/scipy/integrate/__init__.py", line 10,
in <module>
    from odepack import *
  File "/usr/lib/python2.5/site-packages/scipy/integrate/odepack.py", line 7,
in <module>
    import _odepack
ImportError: /usr/lib/python2.5/site-packages/scipy/integrate/_odepack.so:
undefined symbol: daxpy_

I found the solution in Gentoo bug 251165. The problem only occurs when Scipy is built with non-reference versions of BLAS and CBLAS. The solution is to install the reference implementation of BLAS and CBLAS, rebuild Scipy, and then use whichever implementation of BLAS or CBLAS you want. I don’t know why this happens, and I don’t know if it affects distributions other than Gentoo. By the way, Gentoo has a really useful system tool called eselect, which has various modules that are used to choose between different versions or implementations of tools on your system (BLAS, Java virtual machine, kernel sources, OpenGL, etc.)

I spent some time today profiling a Brownian dynamics simulation written in Python to see how I could make it faster before starting some long runs on a Linux cluster. In the sections below, I have attempted to quantify the speed improvements due to various changes. Keep in mind that the actual speed improvement in your code will vary, depending on where the actual bottlenecks occur. See my post about profiling Python code. Another caveat: I am running Python 2.4.4 because it’s installed on our cluster.

Numerical programming frequently requires the use of look-up tables. A look-up table is a collection of pre-computed values. When given an “x” value, the table returns a pre-computed “y” value. Look-up tables can be used to speed up numerical codes, when it is faster to look up a value in the table than it is to compute the value. They are also used when the data in the table cannot be computed–for example, experimental data or averaged results from an ensemble of Monte Carlo simulations. Another application is to compute a value when a function cannot be solved algebraically. Assume that you have a formula for a function q(h). You need the value of h for a given value of q, but the formula cannot be algebraically solved to get h(q). Instead, choose a range of h values, compute the function q(h), and store each value in a look-up table. Now you can get h(q) for any value stored in the table. The major limitation of a look-up table is that it cannot return valid results for any value of q which is outside the range of those stored in the table. Depending on its implementation, the table may be able to interpolate to return values between known points.

NOTE: these instructions are obsolete.  Please see Building NumPy on a 64-bit Red Hat Cluster with Intel MKL. In a previous post I described how I built Numpy with Intel compilers and the Math Kernel Library on a 64-bit cluster. Today I upgraded to Numpy-1.2.1 and I made a few improvements to my install process. Please read the previous post, since I will not duplicate some important information, and then read on. This time, I made use of a site.cfg file. Copy the file site.cfg.exampletosite.cfg and edit. At the end of the file, uncomment the [mkl] section and set the path to your library. Mine looks like:

Last night, in the shower, I realized that my collision detection routine could be even faster. Here is a representative snippet of code from my previous post:

    d2 = (x-self.x[0:i])*(x-self.x[0:i]) + (y-self.y[0:i])*(y-self.y[0:i]) + (z-self.z[0:i])*(z-self.z[0:i])

For some reason, I used the code (x-self.x)*(x-self.x) instead of (x-self.x)**2. Upon further reflection, I realized that (x-self.x)*(x-self.x) computes the difference between array elements twice, and then multiplies the results. Using a “power function” should enable the interpreter to compute the difference only once, and then multiply each element times itself. Here is the updated code, using Python’s power operator:

Today I installed numpy on a cluster. Normally, as a Gentoo admin, I just install things with emerge, and all the details are taken care of automagically. However, this cluster runs Red Hat Enterprise, and I don’t have admin privileges, so I had to install numpy in my home directory. I installed 1.0.4, to match the version used on another system. You may not need to do this for more recent versions of numpy, which may have an improved setup script. The overall process is:

My research group uses Labview 7.1 to write custom data acquisition (DAQ) software. I code everything else in Python, so I need to get data from Labview into Python for processing. Our DAQ program produces Labview binary files, so I had to find a way to read them with Python. Binary files are nice because they are a compact way to store numerical data as compared to ASCI or (heaven forbid) XML, but they are much harder to read. The binary format used by Labview is documented only indirectly, so I had to hack a little. The first thing to realize is that the Labview binary file is a direct dump of the data that was stored in RAM. How Labview stores data in memory is documented here. Indirectly, this documents how binary files are stored on disk. Our DAQ program writes a rather complex “cluster” (Labview’s version of a C structure) to disk. The elements of the cluster are stored contiguously as a sequence of bytes, and there’s no way to know which byte goes with which element, unless you know the size of each element and the order in which they are stored in the cluster. So, the first step is to document the cluster that’s being written to disk. You can use the context help in Labview to view the data type of the wire that leads to the VI that writes the file. With this in hand, you are ready to write Python code. First, make sure you open the file in binary mode: