VTK is an open-source software system for “3D computer graphics, image processing, and visualization” developed by by Kitware. VTK is the foundation of Paraview, an industrial-strength CFD visualization tool that I have found to be very useful. I generate “second generation” XML-based files from my Python code and import them into Paraview for visualization. I am in the process of creating some Python classes to do, and I hope to publish them soon. Until then, I want to share some useful resources. The VTK file formats are specified in this document. It’s a pretty good specification, but it lacks some examples. Soon I will post an example of a valid unstructured, serial .vtu file. Each VTK file includes data from only one time step, so you have to keep track of time yourself (the filename is an easy solution). Paraview can read in data from multiple time steps, but you have to specify them in a .pvd file. This is also an XML file, with the following format: (reference)

I wasted some time today tracking down a bug in one of my programs.  It turned out to be “unexpected behavior” rather than a bug.  I was aware of this aspect of the language, but I made an assumption and got bit.  Read on for a valuable lesson. Python handles integer math differently than floating point math.  If you type a number without a decimal point, Python treats it as an integer.  All math performed only with integers results in integers.For example, 1/2 evaluates to 0 while 1./2. evaluates to 0.5.  If you mix integers and floats, Python will  produce a floating point result (1/2.=0.5), but you must be very careful.  For example, you might expect the expression 4/3*3.14159 to yield a floating point result.  It does yield a floating point number, butnotthe one you were expecting!  4/3*3.14159 yields 3.14159.  What happened?  Python works from left to right.  4/3 evaluates to the integer “1”.  1*3.14159 evaluates to  3.14159.  For comparison, 4./3.*3.14159 evaluates to 4.1887866.  Here’s the problem with this particular aspect of Python: according to the rules of math, 4/3*3.14159 is exactly the same expression as 4*3.14159/3, but in Python they yield different results if you forget the decimal points!  4*3.14159 evaluates to a floating point, so (4*3.14159)/3 yields the “correct” floating point value. Lesson Learned: be explicit about specifyingall floats if you are doing floating-point math!  Sometimes I get lazy and leave a trailing decimal point off of a number when doing a floating point calculation, knowing that the results are “upcast” into floats. Not any more! Note: this unexpected behavior goes away in Python 3.0

“Speed” is a complicated term when used in the context of software.  Does it mean raw speed of execution, or reducing the amount of time until a correct result is obtained?  Python is not the first language that comes to mind when people think of “fast software.”  It is true that pure Python will usually not execute as quickly as the same algorithm directly coded in C or Fortran.  However, when you define speed as “least amount of time until you get the right answer,” then Python is pretty fast.  It is so easy to develop correct code in Python, when compared to low-level compiled languages, that Python is often the fastet route to a correct answer, even if the execution time is longer.  Having said that, there aretimes when code has to execute quickly, and that’s why I will introduce you to profiling Python code. Like all things Python, profiling is easier than you think.  Python 2.4 has the profile module, and Python 2.5 has both profile and cProfile.  cProfile is written in C for lower overhead, and it’s the recommended version.  I am stuck with profile, because the cluster that I am working with still uses Python 2.4.  You can read the docs for more details, but I will quickly outline what I find to be the most helpful usage.  For example, say I want to profile the filerun_sim.py.  I use the following command line:

I had to re-build Numpy because our cluster was upgraded and the Intel compilers and libraries were moved to a different directory.  This turned out to be a half-day affair of trial-and-error.  I learned a few important things, which I will try to list here:

*Delete the numpy-1.0.4/build directory after every build attempt.  Doing “python setup.py clean” isnoteffective.  I kept getting errors about undefined symbols when I tried to “import numpy” on the Python command line.  It was looking for symbols in the old locations, even though I had just rebuilt the code using the new library locations.  It turned out that I needed to delete the build directory in order to force a complete bottom-up rebuild.*The use of “setup.py” from distutils is not well documented online.  The best thing to do is run “python setup.py –help-commands” to get a list of available commands.  Then run “python setup.py <cmd> –help” to get help for that specific command.  You can string commands together on the command line, as I will show in the example below.When you test the new numpy, make sure you arenotin the numpy-1.0.4 directory!  If you are in the numpy source directory, when you import numpy, you will get the message “Running from numpy source directory.” and you will not be able to load any symbols from numpy. On 64-bit architectures, you need to compile position-independent library code.  For some reason, distutils does not do this automatically, and the compilation will fail with an error similar to the following:

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:

Python is a very-high-level language.  That makes it easy to write code quickly, but the program may not be as fast as a program compiled from a lower-level language.  For this reason, many scientific programs are written in Fortran or C++.  However, it has always been my experience that the majority of time on a project is spent in writing, modifiying, and debugging code, rather than executing.  Fortunately, if written correctly, the time-critical parts of Python code can execute almost as fast as compiled software.  Here is an example of a collision-detection algorithm which achieved almost a ten-fold increase in speed when written to use Numpy.

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:

The pickle module provided with Python is so useful that I’m surprised I haven’t used it before. Pickle allows you to save an entire data structure (such as an object) to disk as a binary file in a effortless (and fairly efficient) manner. For example, in my latest project I have created a Monte Carlo simulation that can take quite a bit of time to run. I also need to make multiple runs to get statistics on the results. At the end of each run, I need to dump the resulting data to disk so that it can be read in later by an analysis program. If I had to write data in a format that could be interchanged with other scientific software, I’d use the hdf5 format with the pytables package. However, right now I just need to get something working, and the pickle module is perfect. Here is how I save an object called box:

It’s remarkably easy to spawn a Python thread.  However, before doing so, I caution you that a Python thread is not the same thing as an OS thread.  Python threads run within the Python interpreter, but the Python interpreter always executes in a single process.  The reasons why have already been explained elsewhere, so I refer you to the thread module documentation to learn about the Global Interpreter Lock.  You probably have objections to this state of affairs, and I assure you they have already been voiced by Juergen Brendel and responded to by Guido van Rossum (creator of Python).  Anyway, the upshot is that Python can only utilize one core of a multi-core CPU.  This isn’t such a big deal for me because I’m a scientific programmer, and if I really need to write parallel code it’s going to have to run on a cluster or a grid.  Threads don’t help with that. Having said all that, threads in Python are still useful.  I will detail one example in which I spawn a thread to load a large binary file.  While this doesn’t spread the work across multiple CPU cores, it does enable the GUI to remain interactive while the file loads. All you have to do to create a Python thread is create a class that is derived from Thread. In the example below, I derived a class called Loader, which “wraps” the function that actually reads the binary files. The __init__ method accepts the filename and other options as arguments. The run() method is required. Don’t call run() directly–instead, call the start() method (inherited from the base class) to start the thread.