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

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.

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.

I have finally gotten back to programming in the last couple of days.  Our project has finally started to generate a lot of data, so I’ve been refactoring and improving my code that reads data stored in LabView binaries.  Today I spent a lot of time creating a GUI for browsing data.  Arguably, this wasn’t the best use of my time, but I learned a lot about multi-threaded Python GUI programming with wxPython.  You can find a gold mine of information on the multi-threaded wx programming at the wxPython wiki.  Because the LabView binary data has to be read sequentially, and the files are rather large, it takes a long time to read in a file.  I spawn a thread to handle the file reading,  while allowing the GUI to remain responsive.  The thread posts messages to the GUI window, which are used to update the user on the status of the file reading operation.  When the file is read, a final message containing the data is posted to the window.  It’s really pretty slick now that I’ve figured out how to do it. I will soon post a clever scheme to capture text output from the file-reading function, and display it in the GUI, without making substantial changes to the file-reading function.

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:

The site is back online after a server move.  Actually I can’t blame the server move for the downtime, because I went out of town and then had a lot of catching up to do, and as a result I didn’t switch the domain to point to the new server. I haven’t been blogging much because I haven’t done much software development, Linux admin, or lighting design lately.  I have been busy in the lab at work and I’ve been supervising a trainee lighting operator instead of running my own shows.  We haven’t been moving forward with plans to purchase a new lighting control console, so no update on that, either.  Hopefully, I’ll have something interesting to post soon.