SciPy India conference in Dec. 2009

I'm excited to announce that the first "Scientific Computing with Python" conference in India (http://scipy.in) will be held from December 12th to 17th, 2009 at the Technopark in Trivandrum, Kerala, India.

The theme of the conference will be "Scientific Python in Action" with respect to application and teaching. We are pleased to have Travis Oliphant, the creator and lead developer of numpy as the keynote speaker.

Here is a rough schedule of the conference:

Sat. Dec. 12 (conference)
Sun. Dec. 13 (conference)
Mon. Dec. 14 (tutorials)
Tues. Dec. 15 (tutorials)
Wed. Dec. 16 (sprint)
Thu. Dec. 17 (sprint)

The tutorial sessions will have two tracks, one specifically for teachers and one for the general public.

There are no registration fees.

Please register at:

http://scipy.in

The call for papers will be announced soon.

This conference is organized by the FOSSEE project funded by the Ministry of Human Resources and Development's National Mission on Education (NME) through Information and Communication Technology (ICT) jointly with SPACE-Kerala.

Sage notebook, Mayavi2 and 3D interaction on the web

Thanks to a cool open source x3d viewer along with its browser plugin -- freewrl, I can now use mayavi2 in a sage notebook and also get 3d interactivity! To see this in action see here. You should see something that looks like this:

The view on the bottom supports full 3d interactivity, which is really quite cool. This works seamlessly with the debian package on the freewrl site. I just installed the packge on my Ubuntu 8.04 machine and it works really nicely. The only trouble is that for more complex visualizations, the x3d file exported can be quite large but this can be partially alleviated by perhaps using mod_deflate and similar options on the server.

Now the x3d embedding above is a hack currently since I manually embedded the necessary HTML tag into the sage worksheet but as can be gleaned from this thread, this should be trivial to add to Sage proper. This opens up an exciting new possibility for mayavi2 on the web.

Mayavi2 in a Sage notebook on the web!

I've been pretty hopeless at blogging (or even responding to comments on my blog) but this was too good to not blog about. Ondrej Certik has setup a Sage notebook on the web with support for Mayavi2. This lets you generate pictures using mayavi2 offline, from anywhere on the web and see them as static images. No interactivity yet but that is technically possible and may happen sometime.

Ondrej blogged about this very briefly and mentioned his struggle on getting mayavi2 working in sage. Here is a link to his first attempt demonstrating a mayavi plot on the hpfem notebook:

http://nb.hpfem.org/home/pub/12/


I made a little worksheet with a few of the standard mayavi test cases here:

http://nb.hpfem.org/home/pub/16/


If your browser shows you a text rendering of the above, download the Force HTML bookmarklet from here and click on that, you should see the page nicely.

Quite a bit of work went into this before we got this far. I spent several hours figuring out how to get VTK built reliably so it can work for this. There was a minor bug in mlab that was fixed by Gael and me. Ondrej then figured out (after much pain) how to get mayavi2 and sage working nicely together. Jaap Spies' very cool mayavi2 packages for sage were also a big help I believe. Ondrej then ran into trouble with memory issues with the virtual machine he hosts this on. Once that was resolved this became possible.

Needless to say, this opens up some very interesting possibilities and is very cool.

Python list (in Cython) vs. NumPy

Taking my previous benchmark a little further I decided to see how well iterating over a Python list of doubles compares with using NumPy arrays. Here is an extremely simple example that implements the sum function in Cython and compares the result with NumPy's sum method.

Here is the Cython code:

# --- csum.pyx ---
def csum(list array):
  cdef int i, N=len(array)
  cdef double x, s=0.0
  for i in range(N):
      x = array[i]
      s += x
  return s

Here is a setup.py to build it:

# --- setup.py ---
from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
setup(cmdclass={'build_ext': build_ext},
    ext_modules = [Extension("csum", ["csum.pyx"])])

To time it in IPython I created a simple file called test.ipy like so:

# --- test.ipy ---
import csum
import numpy
for i in [10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 250000000]:
  print '-'*80
  print 'N =', i
  a = numpy.linspace(0, 1, i)
  b = a.tolist()
  print "Cython:",
  %timeit csum.csum(b)
  print "NumPy:",
  %timeit a.sum()

I run it using IPython and here are the results (formatted a little):

N         Cython       NumPy
10          534 ns       10.1 micros
100         1.76 micros  10.8 micros
1000        15.3 micros  19.3
10000       150 micros   101 micros
100000      1.75 ms      933 micros
1000000     19.7 ms       9.24 ms
10000000    198 ms        92.1 ms
25000000    499 ms       231 ms

This was done on a P4 3 Ghz machine and clearly lists do quite well. At small sizes they outperform NumPy and at really large sizes they are about twice as slow. This is very good considering how general lists are.

Python vs. Cython vs. D (PyD) vs. C++ (SWIG)

In April 2008 there was a thread on the scipy-dev list regarding the inclusion of Cython code in SciPy. In that thread, I mentioned a particular use case of interest to me -- creating and manipulating an array of objects (rather than arrays of elementary data types) and being able to do that with Cython easily and efficiently.

The problem I was considering is a simple one. I create a list of "Vortex" objects and compute (naively) the velocity of a collection of these particles on one another. This is an O(N^2) computation since every particle influences every other. The idea was to create simple OO code to be able to perform these computations. Here is an outline of the Python code for doing this:

class Vortex(object):
 def __init__(self, pos=0.0, strength=1.0):
    # ...
 def eval_velocity(self, pos):
    return -1j*self.strength/(2*pi*(pos - self.position))

class VortexManager(object):
 def __init__(self, vortices=None):
    # vortices is a list of vortex objects.
    self.vortices = vortices

 def set(self, pos, str):
    # ...

 def velocity(self):
    for va in self.vortices:
        vel = complex(0, 0)
        for vb in self.vortices:
            if va is vb:
                continue
            else:
                vel += vb.eval_velocity(va.position)
        va.velocity = vel

Very straightforward code. Now, back in April I implemented this in pure Python, C++ and wrapped the C++ code to Python using SWIG. I also implemented it in D and wrapped that using PyD. I found that D was about 1.7 times slower than C++. C++ was about 300-400 times faster than the pure Python version.

I attended Robert Bradshaw's Cython tutorial at SciPy08 and really liked it. About 10 days ago I finally found the time to create a Cython version and the winner is ...

Cython!

I've put up all of the code here. To use the code, untar the tarball and do the following:

$ cd cpython_d_cpp
$ python setup.py build_ext --inplace

This requires SWIG, numpy and Cython to build. If you have PyD installed it will build the PyD extension also. To test the performance do this:

$ python test_perf.py

This produces the following output for me on a P4 (3 Ghz):

dee(4000): 1.87730193138
(1411.53285812+1411.53285812j) (945.091286479+945.091286479j)
swig(4000): 1.10782289505
(1411.53285812+1411.53285812j) (945.091286479+945.091286479j)
cython(4000): 1.15034103394
(1411.53285812+1411.53285812j) (945.091286479+945.091286479j)
Pure Python(200): 1.14771318436
# N SWIG Cython Ratio
1000 0.071 0.069 0.967
2000 0.283 0.274 0.968
3000 0.638 0.619 0.970
4000 1.135 1.100 0.970
5000 1.767 1.720 0.973
6000 2.517 2.473 0.983
7000 3.474 3.370 0.970
8000 4.541 4.403 0.970
9000 5.698 5.575 0.978
10000 7.000 6.879 0.983

The first few numbers just test one single case of 4000 particles. D is slower than both C++ and Cython. Python is dog slow (or donkey slow as I like to say it)! For some reason I was getting segfaults when I tried to test the D wrapper for more than 3000 particles. On my Macbook the Cython version is actually 30% faster than the C++ version and on a Sempron 2800 Cython is about 25% slower. So different machines produce different numbers. However, C++ and Cython are both in the same ballpark.

What I loved about the Cython code is that I use a Python list to manage the Vortex objects. This shows that we can use the normal Python containers to manage objects. This is extremely convenient. This isn't very surprising either since Python containers are also heavily optimized. "cython -a" was a huge help when getting things done. For more details on the Cython code look at the tarball.

Clearly, if you are building code from scratch and need speed, Cython is an excellent option. For this I really must congratulate the Cython and Pyrex developers.

Mayavi Screencast now on blogger

Finally, I've managed to upload the screencast I posted last about on blogger.
If you are on a Mac or on Windows and downloaded the file I put up earlier, you might need to install the theora component from here.

I have no idea how it will eventually show up with the processing that goes into it. Depending on how this goes I'll put up new screencasts either here or perhaps at Showmedo. Anyways, here it is:

Mayavi2-3.x screencast

Here is a screencast (16 Mb Ogg Theora) of some of the new Mayavi2 features in the 3.x series.

http://www.aero.iitb.ac.in/~prabhu/tmp/videos/mayavi2-3-screencast.avi


Its an ogg theora video. This is my first screencast so its going to be a little rough. There is too much noise in the sound recording and I'll try and fix that next time. BTW, I used the excellent istanbul for the video recording. I recorded the sound track separately and mixed the two with a gst script I got from here.

I was unable to upload this video on this blog because the IITB network at this time is unbelievably slow. If you know of another place I should consider uploading the video or would like to host the video elsewhere please let me know. Feel free to download it and host it elsewhere.