SciPy - the embarrassing way to code

After spending most of the last night trying to model and fit a bunch of data with scipy, my main bottleneck was my understanding; working out what I needed to do. Once that was done, everything could be implemented in a single call. I don't think that's a particularly bad bottleneck at all.

There was some good discussion about this last time this was posted, and also current discussion over on reddit.

http://news.ycombinator.com/item?id=182676 http://www.reddit.com/r/Python/comments/a2xzj/the_embarrassi...

I had the exact same experience recently while writing a dynamic programming simulation in Python. The naive (and elegant) solution to the problem consists of a 20-line recursive function which took me about an hour to write. Unfortunately function calls are very expensive in Python, and I my sim required a couple billion of them, so I started searching for ways to eliminate the recursion. 3 weeks and a deep foray into Numpy later and I had it.

The second implementation is about the same length and replaces all of the simplicity of the recursive solution with a bunch of five-dimensional arrays. Never before did I have to think as deeply about vectorized functions, array broadcasting, dimension matching, etc. The non-recursive version is about an order of magnitude faster, and what's better, using the numpy operators means nearly all of the processing is offloaded into low-level C routines, so the Python overhead is effectively nil. On the downside the code is extremely difficult to follow now... trying to visualize what's going on in 5 dimensions is a real pain.

I did some SciPy work in computer vision applications, and it was both a blessing and a curse. About 80% of the code could be reduced into matrix operations, which provided the 'embarrassing' moments that the article mentions. However, I'd hit something like Hysteresis Thresholding in the Canny Edge Detector, which traces through the image looking for faint edges, and spend a day figuring out how to write a Python extension against the Numpy API. This was necessary, because any code that iterated over 1000 x 1000 matrices in Python that I wrote took well over 10 seconds, instead of the < 100 milliseconds it should take.

As a nublet, I probably didn't end up saving much time over C++, but I could see that in the right hands, this tool would be a lifesaver

scipy still doesn't have a good submodule for doing manipulations with multivariate gaussians, which is a HUGE part of machine learning applications.

I'm referring to the analog of R's mvtnorm library.

Yes, you can use rpy2. Yes, you can roll your own dmvnorm (but be careful about the degenerate situations with zero eigenvalues that always arise in high dimensional problems with real data). It is significantly more of a pain to roll your own pmvnorm, because you are now talking about efficient computation of high dimensional integrals, which always seems to involve a trip to Numerical Recipes or the like (see the Genz and Miwa references when you do "?pmvnorm" in R).

Basically this is a huge thing to lack.

PS: yes, if you grep the codebase, you can find this: http://www.scipy.org/doc/api_docs/SciPy.stats.kde.gaussian_k...

But you need to write wrappers because it doesn't provide a basic interface to the multivariate gaussian density.

I could go on about other parts of scipy -- it has some advanced stuff but lacks a lot of basics...

"It's just a matrix" and all the comments sound like SciPy is some kind of APL for the rest of us, is that impression correct? There you too can achieve incredible conciseness once you manage to formulate you problem using some matrices and vectors with funny dimension. (And you get the benefit that the resulting code doesn't look too simple, so you get something to show off witth in return for your effort.)

This is precisely why people are afraid of math - the more work you do, the smaller the end result. E = mc^2 is actually an approximation, but even the full set isn't very complicated, and imagine how difficult it is to get there. The other truth is that less code is generally faster, but not as small or as fast as a better algorithm, which is in turn not as small or as fast as a mathematical equation. The poster talks about moving across this space, which is actually a great accomplishment... but will probably get the best reception among great programmers or mathematicians.

I once spent around ten weeks writing a C program to numerically integrate the Schrödinger equation for an electron in a particular magnetic field. When I started looking at results, they looked eerily like trig functions. A little digging revealed that the problem was stated and solved (exactly) in one of the exercises to Quantum Mechanics by Landau and Lifshitz.

I sheepishly told my advisor at the time what had happened, worried he would think me an idiot, but he was thrilled: "You solved the problem! That's what's important."

I think the author is taking the wrong thing away from this, all those little realisations are things that if he comes across in the future he will have a better understanding of what to do.

In other words rather than become embarrassed after realising a months work boiled down to just a few lines think about how much you learned to get to those few lines.

Just curious, can SciPy also automatically take advantage of cell processors or graphics processors? If it is all just matrix operations, it seems doable?