Having recently upgraded emacs to the latest version (I was running the last release of CarbonEmacs), I took the opportunity to clean up my .emacs file and add/enable a few new packages:

• smex: provides ido-like completion of M-x commands
  
• ropemacs: this looks really useful for python programming. Docstring lookup, name completion, refactoring.
  
• darcsum: interactive darcs recording
  

Packages I've been using for some time:

• ESS: Run R in an emacs buffer. Easily run bits of code in the running R process.
  
• AucTeX: edit and compile LaTeX documents
  
• doxymacs: package to easily insert doxygen comments
  

Deprecated packages:

swbuff: cycles between buffers with a single keystroke, similar to how you can move between tabs in a browser with C-pagedown and C-pageup. I now use ido for this, using the following lines in my .emacs file:

Why are (neuro)scientists such terrible programmers? I just spent most of the morning trying to add some simple functionality to some code from a fairly prominent computational neuroscience lab, and the code was so bad, I began to entertain serious doubts about the method and a whole string of papers that have been published using the method. I briefly debated rewriting the whole thing from scratch, but I have about twelve other projects on my plate now, and the benefits of challenging a fairly accepted method are pretty low, unless you can also upset a few cherished hypothesis at the same time. So I crammed my patches into the spaghetti, prayed it didn't break anything, and started the analysis, which only makes a rather minor contribution to my results anyway. No doubt better programmers than me have felt and done similar things, and in all likelihood the problem is much larger: science is so competitive now, with everyone looking for big, new, flashy results, that verification has become a really low priority. There's not a lot I can do about that, but at least in the realm of computational techniques it may be worthwhile to outline the problem.

What makes code good, and why is it important? Good code is first and foremost, easy to interpret by someone who didn't write it. Anyone with a basic understanding of programming, and the scientific question being addressed, should be able to tell what's going on in every line. This goes far beyond commenting; in fact, really good code doesn't need to have a lot of comments, because the expressions are largely self-explanatory. But regardless of how it's achieved, when code is readable the methodology is clear. It should be obvious that this is important to good science, because it allows your colleagues to evaluate what you did and extend your work. But let's be honest: not everyone shares those goals. A cynical explanation for why there's so much bad scientific code is that the people releasing it want to maintain control over the method, for commercial or competitive reasons. I don't think this motive is all that common in neuroscience, so I prefer to think that it stems from a lack of education, which I am willing to do my part to try to dispel. Toward that end, I'll be posting articles here detailing some common mistakes and mistaken attitudes.

In the long run, as computational techniques become more complex and more common, I would propose that code needs to be treated as a part of a paper's Methods section. If a program hasn't been used previously, it should be subject to review, and even if the method isn't new, authors should be required to deposit the code they used in a given paper somewhere online. Many journals make similar requirements for gene sequences. I think this would do a lot to improve the quality of programming, and it would keep researchers from purposefully obfuscating code to keep one foot in the door for commercial exploitation. But of course there is a much larger discussion about how well scientific and pecuniary motives can coexist in research.

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.