Logical indexing with Numpy

I just discovered one nuance of Python's Numpy package which, given my Matlab background, is a bit unintuitive.

Suppose I have an array of numeric data and I would like to filter out only elements whose values lie between, say, 10 and 60, inclusive. In Matlab, I would do this:

filterData = data(data >= 10 & data <=60)

However, logical operations with Numpy sometimes require special functions. The equivalent expression in Python is

filterData = data[np.logical_and(data >= 10, data <= 60)].

np.logical_and() computes the element-wise truth values of (data >= 10) AND(data <= 60), allowing me to logical index (or fancy index, as pythonistas say) into data.

Numpy's memory maps for increased file read/write efficiency

I mainly use Matlab and Python when I need to analyze and visualize data and run simulations. As I mentioned in my previous post, I like to document when I learn something new that makes my programs written with these tools cleaner or more efficient.

I've just discovered Numpy's memory map object for efficiently reading and writing to very large arrays in files. A memory map is used to access only part of an array that is saved in a binary file on a hard disk. The memory map may be indexed like any other array that exists in memory, but only the part of the array that you are working on will be loaded into memory. This is similar to the idea of generators that I wrote about in my previous post.

As an example, I create a memory map object called inFile which interacts with a file whose location is stored in the variable fName:

inFile = numpy.memmap(filename, dtype = 'float32', mode = 'w+', shape = (3,4))

The mode argument 'w+' tells the program that I want to create the file if it doesn't exist or overwrite it in case that it does. Printing the contents of inFile at this point will result in a 3x4 array of zeros being displayed on the screen. When I want to write to this memory map, I can write to a portion of it like:

inFile[0,3] = 3.2

which writes 3.2 to element in the first row, fourth column of the array. When I delete inFile, the contents will be flushed to memory:

del inFile

Reading parts of this array from fName is done by expressions like the following:

inFile = numpy.memmap(filename, dtype = 'float32', mode = 'r', shape = (3,4))

 Now, to move some of the file's contents to an array in memory, I only need to type

myNewVariable = inFile[:,2]

which, for example, moves the third column of the file to the array myNewVariable.

I've noticed an increase in speed of my programs by using memory maps, though I haven't tested exactly how much time is saved.  

 

Increase a Python program's memory efficiency with generators

I've recently been working with very large data sets (more than a million data points) and have encountered a serious reduction in the efficiency of my Python programs' computations. One reason for this is that I have been reading a large data file into memory all at once, performing my computation on this data file, and then moving onto the next file, all the while keeping my calculation results in memory so that I can plot them at the end.

One tool I've discovered for increasing the efficiency of these types of operations is the idea of a generator. Now, I realize that these are well-known in Python circles. However, I am not primarily a programmer, so I am not fully aware of the tools available to me when I write programs. Hence, I sometimes use this blog as a notebook so I can easily find these tools again in the future.

Simply put, a generator is an object that executes commands and yet may be iterated over like the elements in the list. One example, shown here, is the readline() method. The code looks like this:

inFile = open('sampledata.txt')
w = inFile.readline()

w is a generator function that reads only one line in the file into memory. Each time readline() is called, w is assigned the string representing the next line in the file sampledata.txt. According to the above link, an object that offers the abilities to load just parts of an object into memory and yet still iterate over all the parts is a generator.

Another nice description of the use of generators is given here. And, because I am having some trouble plotting all this data, a discussion of possible solutions to plotting large data sets is presented here at Stack Overflow.

The importance of knowing your random number generator

I recently gained access to the Stokes computing cluster at UCF and ran some long-time simulations that I had written in Matlab. These were not in parallel, but rather utilized each node to run a separate simulation. Each simulation took a few days to run.

I was a bit irritated to find the output from each simulation was exactly the same, despite having run on separate nodes. The reason, I discovered, is that Matlab will generate the same string of random numbers every time you start a new instance.

I tried this experiment with Matlab R2011a. After opening a fresh instance of Matlab, I entered and received the following in the Command Window:

>> randn(1,5)

ans =

    0.5377    1.8339   -2.2588    0.8622    0.3188

I then closed Matlab and reopened it. Entering the same command gives the same output as before:

>> randn(1,5)

ans =

    0.5377    1.8339   -2.2588    0.8622    0.3188

Oh, boy. I performed the same procedure on a different model of computer in our lab and received the same output. I also tried this with Python 2.7 and Numpy's random.rand() method. The output is different each time I start an iPython session, so the seed must be automatically shuffled at the start.

Lesson learned. If you're working with random numbers in Matlab (at least with R2011a), be sure to randomize the seed with the command rng shuffle. Now to go redo a few week's worth of work...

A resource for volume visualization in Matlab

As many engineers may already know, visualizing scalar or vector volume information is difficult. Examples include the heat distribution in a room (scalar) or the electromagnetic field around a complex-shaped scatterer (vector). There is a good set of resources on how to do so in Matlab here: http://blogs.mathworks.com/videos/2009/10/30/basics-volume-visualization-29-examples-of-scalar-and-vector-fields/. I haven't tried volume visualization with any Python packages (see Python x,y for a collection of useful scientific packages), but if I find the time to play with it, I'll compare and contrast the two.