Generators are Python functions that have a yield statement instead of a return statement. Or rather, calling a function that has a yield statement returns a generator, which is an iterable object that uses the underlying function to produce values. Each time you call next() on the generator, it yields the next value in the sequence, or throws StopIteration, which indicates that it's out of values. Python's famous list comprehensions are actually generators. For example, you can create a list of square numbers with a command like [x**2 for x in range(5)]. If you leave out the brackets you get a generator:



Python 2.5 added some more sophistication to generators that allow them to be used as coroutines (PEP 342). For instance, generators now have a close() method that causes the generating function to throw a GeneratorExit, which can act like a signal to tell the generator to clean itself up. It's also possible to pass an arbitrary exception to the generator, or to pass values back to the generator.

The advantage of using coroutines is that it allows you to modularize your code without having to break out a lot of heavy object-oriented features. For example, I often have to write scripts that plot a bunch of little figures. I often want to break the figure across multiple pages, so I generate axes until the figure is full, then write the figure to disk and start a new one. There's some bookkeeping involved in this; I have to create the initial figure, then check after every plot if the figure is full and take the appropriate actions if it is. All of that is extraneous to what the script is actually doing, so we want to factor that out. There are many ways of solving this, but a generator is remarkably simple:



Using this generator is fairly simple. The following code plots 4 figures with 5 panels each, with only a single line of extra code to initialize the coroutine.



Several things to note about this function:

• The grid function is itself a generator that spits out Axes objects. If the user just specifies the number of columns and rows, I create a generator that uses matplotlib's add_subplot() method. But I have other generators I can pass in that create irregularly spaced grids, so this is highly modular.
  
• Remember that because of Python's closure rules, the grid() generator uses the values of nx and ny defined in the main function. Thus, even though the signature for grid is grid(fig), it can still access these values. Note that if you define the gridding function in a different scope, it binds to the variables in that context. In other words, you could modify the grid dynamically by changing those variables in the calling scope, although this would probably be an undesirable side effect since the generator doesn't tell you when it's started on a new figure.
  
• The grid generator stops yielding axes when the figure is full, at which point the post-processing function is called and a new figure is generated.
  
• The try/except block catches calls to the close() or throw() methods of the generator and makes sure figfun() gets run on the final figure.
  

The motivation for using scons to automate analysis workflows is to formally specify dependencies in how the analysis gets done. This is the idea behind "reproducible research" (a term clearly coined by non-biologists; I prefer "reproducible analysis"). You run one command and it turns your raw data into figures that are nearly ready for publication, removing a lot of potential sources of error. I've used Sweave for a couple of projects. It's a system that allows you to embed R code chunks in a LaTeX document, and the calculations and figures generated from the embedded code get placed directly into the resulting PDF file. It's a powerful idea, but it doesn't lend itself well to the development of multi-stage analysis workflows where the output of one step in data analysis flows into the next step. Change one thing and you have to rerun the entire script. It makes a lot more sense to design each step as a modular component and formally specify the dependencies between the steps.

In the previous post, I described how to use a custom builder to add python scripts to a scons workflow. You can do the same thing for R scripts and Sweave documents. With a little regular expression kung-fu you can even get scons to recognize which R scripts and data sources are being imported.



You still have to manually specify what each script outputs (except for the Sweave builder, which knows the output will be a .tex file), for instance:

I've started using scons to manage my data analysis workflows. A lot of the work is done by python scripts that read in a bunch of data from one more sources, crunch it, and spit out a new table. So there is a dependency not only on the input data but on the script and any modules it imports. You can get scons to run a script using a simple Command builder. For instance, you have some "script.py" that expects a command line like "script.py input1 input2 output".



This works fairly well, but if script.py imports another module (e.g. with functions common to a bunch of scripts) you have to manually specify that dependency. But with a little extra code you can get scons to automatically scan python scripts for import statements and include any imports from the local directory as dependencies. I also like to use an Emitter that will move the script to the front of the list of dependencies, so I don't have to worry about what order I specify them in.



Now you can use the custom builder as follows, and scons will recognize any modules script.py depends on.

R doesn't have great string or set manipulation functions, but you can accomplish a lot by using factors and apply(). For instance, I had a column in a data frame that consisted of 3 two-character tokens concatenated together (e.g. 'abdqbk','cbdabb') that represented sets of 3 sounds played in order. I needed to recode this as a set of two factors, one indicating which set of symbols had been used in a particular trial, and another indicating the order they had been played. Solving the problem turned out to be an interesting illustration of how your tools constrain your thinking, and in this case the result wound up being fairly elegant (IMHO). Whereas if I had been working in Python I probably would have written a couple of loops to run through the strings, one to find the unique sets of symbols, and the other to assign unique values to each permutation.



strmsplit() is another little bit of apply() magic that will split a string at a set of fixed cut points. The code is from one R tip a day.

One of the features I appreciate most about numpy is array broadcasting, probably because I hated having to use repmat() in my MATLAB days. With MATLAB, if you want to add the values in a vector to each column of a matrix, you have to replicate the vector to the same size as the matrix first, which is a colossal waste of memory (MATLAB does support something called lazy evaluation, but I'm not sure it applies to these sorts of expressions). In numpy you can simply add the vector to the matrix; however, the rules are a little tricky. There's an excellent FAQ on broadcasting here - how to specify which dimension should be matched between the two arrays, and when broadcasting is actually more inefficient.

Numpy supports structured arrays, which are the nearest thing to R's data.frame class. Data are organized into fields and records. Each field (column) has a name and data type, and each record (row) has a value for all the fields. Columns are indexed by name, and rows are indexed by integers. Recarray objects can be generated from nested Python iterable objects using numpy.rec.fromrecords:



Note that the 'price' field has a float data type because one of the records has a float value, and the field is promoted to the most general data type. For more precise control over field data types, fromrecords() takes a format argument, which is a comma-delimited list of format strings. For instance, to force 'price' to be an integer, call D = numpy.rec.fromrecords(D, names='quality,price,size', formats='S4,i4,i4')

For reading and writing recarrays, use matplotlib.mlab.rec2csv() and matplotlib.mlab.csv2rec(). The format of each field can be specified using a dictionary. There are a number of arguments to both functions that can be used to control how the data is read in (e.g. delimiter, is the first row a list of field names, etc), most of which are documented. The rec2csv() function always outputs field names as headers. To avoid this behavior, or to avoid having a dependency on matplotlib, use numpy.savetxt()

A useful and underappreciated MATLAB function is accumarray(), which aggregates a vector of data elements into a multi-dimensional array based on a separate set of index vectors. It's essentially the same as the R function xtabs(), except that the indices in accumarray() have to be numeric. I found myself needing something like this function in python when I was writing an algorithm for calculating reassigned spectrograms (more on this later; I'm working on a little monograph on time-frequency representations of birdsong). I couldn't find anything, and wound up writing the algorithm in weave/blitz -- which I probably would have done anyway to tweak the performance.

Later I discovered that there is an equivalent to accumarray() in scipy: it's the scipy.sparse.coo_matrix() constructor. Here's the example from the docstring:



The main downside to this method is that it only works with 2D arrays. I also don't know anything about the performance. Probably if your matrix is fairly sparse (i.e. lots of zeros) this is the way to go. If you've got more or fewer dimensions, or the output is dense, you're probably better off just writing the loop in C++.

Contour plots are useful for visualizing 2D arrays of numbers where the values vary relatively smoothly over space. A less common use is to annotate regions in a plot, in which case the identity of each region is typically coded by an integer. If you have multiple regions, the values aren't going to vary smoothly, and the built-in matplot contour command won't look quite right: higher-valued regions will have multiple borders around them. The solution is to make a separate plot call for each unique value of the label. An optional gaussian filter can be used to smooth off region edges.





The following code, or something equivalent, is necessary to ensure that the contour colors cycle through different values for each level. I find these particular colors pleasing and fairly easy to distinguish, but your mileage may vary, especially if you need to annotate more than 15 different regions.

Let's say you're making a nice class and you want to expose some constant property of the class so that code that uses your class will know, for example, what kind of data type the class deals in. In Java or C++ you create a constant static member and you're set. In python you do the same thing, except these variables are called class variables. But how do you do this when you're defining a Python class in C or C++?

The answer is probably totally obvious to anyone who uses python a lot, but it took me more than 15 minutes to find the answer, so it goes here: you put it in the class dictionary. With normal python objects both the instance and the class have a dictionary that contains the attributes of the instance. The instance dictionary takes precedence, so if you define a class variable and then subsequently assign some value to that variable for an instance, it gets stored in the instance dictionary, and that's the value you get if you access that attribute in the future. Types defined in C extension modules generally don't have instance dictionaries, so anything you put in the type dictionary winds up being read-only.

Here's the line I used (in the initialization function) to create the class variable _dtype and assign it the numpy dtype associated with the data the class processes.

I've been using functions in the scipy.io module to read and write binary and ASCII data, specifically the fread(), fwrite() and write_array() functions. These have been deprecated due to some new functionality in numpy. I'm not sure how bleeding-edge this is, since I have been pretty erratic about installing numpy/scipy from svn and from official releases, but since there's no easily locatable documentation on the web at present, here are my notes. In each example, fp represents a file object (opened in the correct mode, of course), and the first line is the old scipy call that should be replaced with the new numpy call in the second line.

Writing binary data:



Writing ASCII (i.e. tabular) data. I do this all the time when exporting large data sets to R. Write a header
line to the file first, and adjust the format string for your data.



Reading binary data (PCM sound data in this case):



An alternative is to use numpy's builtin memory-mapped file access, described here. This method supports complex and structured data types. For instance, I could replace my function for reading PCM data with the following line:



Python should handle the garbage collection on the object (data) when its reference count drops to zero, but I haven't tested this extensively. And of course you have to be careful about writing to the array, because it will affect what's on the disk. Setting the object to read-only mode (as above) will prevent this from happening, but a RuntimeError will be thrown if you try to modify the array.