Question about scientific calculations in Python

[Paul Hughett]

Martin Kaufmann <martinkaufmann at yahoo.com> wrote:

: I'm trying to write a program to calculate diffraction patterns for
: electron diffraction on clusters (100-10'000 atoms). So far I used
: Python only for CGI programming and system tools (on a rather low
: level...). Therefore I don't know whether it's reasonable to program
: such a project in Python (I really like the language). Presently my
: man concern is speed as the main function of the program is a nested
: loop, the inner over several 1000 bins of a histogram and the outer
: over several 1000 values of scattering factors. In the end I should be
: able to simulate diffraction patterns and compare them to the
: experimental data (i.e. run this function several times).

You might want to take a look at bblimage, which is an extension to
Python written in C; it emphasizes high performance, but the set of
operations implemented so far is oriented toward medical image
processing.  It supports pretty much all the basic operations (+,-, ...
sin, cos, exp, etc) on n-dimensional images in the usual integral and
floating point types (but not complex, alas; I haven't starting playing
with FFTs yet).  The URL is

         http://www.uphs.upenn.edu/bbl/bblimage

: Now my questions: Would it be best to

: (a) write the whole program in C/C++ (I know only the basics of C but
: it would be a good "excuse" to learn it...)?
: (b) write the main program in Python but the heavy calculations in C
: (I played today with scipy.weave -- the speed is much better but I
: didn't really understand what I was doing...)?
: (c) program it in Python and don't care about speed (or buy a new
: workstation...)?

It really depends on how many diffraction patterns you need to do.  It
doesn't make sense to spent two months learning C to save even 99% of
a ten-hour computation that you're going to do only once.  (Unless
your real purpose is to learn C.)  On the other hand. if you have to
do a few hundred patterns, then it might be a good investment of time.

For what it's worth, my own experience is that (b) using C for the
heavy lifting and Python to glue the pieces together can be very
effective if you implement the right piece in the right language.  I'm
using bblimage for registration and segmentation of hundreds of volume
images at half an hour apiece.  The current code runs almost as fast
as plain C (about 4 times faster than a competing C++ implementation
not particularly optimized for speed) but the Python code that glues
it all together is easier to modify and experiment with than a C++
implementation would be.

One general comment: The performance killer in many applications like
this is doing the innermost loop in the interpreted language (Python,
Matlab, etc).  If you can recast your algorithm as operations on
arrays rather than individual numbers, and your language does the
implicit loops over the elements efficiently, you typically get 10x or
so improvements in speed.  


Paul Hughett