Sunday, September 21, 2008

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.

Saturday, September 20, 2008

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:


Friday, September 19, 2008

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.

Thursday, September 18, 2008

Announcing the Mayavi2-3.x series

This is a long overdue announcement. ETS-3.0.0 was released just before the SciPy conference in August 2008. Mayavi-3.0.0 was released as part of this. There are a huge number of significant changes to Mayavi as compared to 2.x. Note that we are still calling this Mayavi2-3.x.y since the Mayavi2 represents a departure from the older Mayavi-1.x series and Mayavi2-3.x is simply the next version of Mayavi2. The full details of the changes from the 2.x series are documented in the CHANGES.txt file in the mayavi documentation directory. Here is a summary of the major changes to this series thus far. The current release is 3.0.3.

Core Mayavi:
  1. I've added all the modules and filters that were available in Mayavi1 into Mayavi2. The only module I didn't port is the Locator module which didn't seem very useful. Mayavi2 now has more modules and filters than Mayavi1 had. Now there isn't an excuse to continue using Mayavi1.5.
  2. Users can now right-click on the nodes on the tree view to create new sources, filters and modules.
  3. The menu entries for the modules and filters (on the app and on right-click) are all context sensitive. So if your data doesn't support a particular module you shouldn't be able to add it from the UI.
  4. The file->open menu is far cleaner and exposes just one "Open" item that automatically lets you open any supported data.
  5. Added a toolbar to the engine view that offer icons to make it easy to add new sources/filters and modules. Special "Adder nodes" are added to the tree view when a scene/source is empty that makes it easy for new users to use the mayavi pipeline.
  6. Added a -o/--offscreen option to the mayavi2 application so you can run mayavi offscreen if your VTK version supports it. This in combination with the -x command line option makes for a powerful combination.
  7. New and much easier extension mechanism for the mayavi library and app via a user_mayavi.py and site_mayavi.py.
  8. Added a tvtk_doc module/script that lets you search through the TVTK classes (with and/or keyword support), this is similar to Mayavi1's vtk_doc script.
  9. Added a GenericModule that makes it very easy to put together a bunch of components/filters to create a new module.
  10. Added Optional, Collection filters that let you easily build filters out of combinations of existing components or filters.
  11. Added a new SetActiveAttribute filter that lets you choose the active scalar/vector/tensor attribute, this lets you do neat things like plot iso-contours of one scalar on top of the iso-contour of another, see examples/mayavi/contour_contour.py for an example.
  12. Gaël sphinxified the documentation to make it look much nicer and fully searchable.
  13. Better and more complete testing, these are unfortunately integration tests currently and will slowly be made into proper unit tests.
  14. The mayavi2 application and plugins are now ported to use Envisage3 which is much cleaner and nicer to work with than Envisage2.
  15. There is now a full-fledged preferences framework for Mayavi (to access the preferences use, from enthought.mayavi.preferences.api import preference_manager).
  16. Some parts of the API and file organization has been cleaned up. This is mostly related to the location of some modules, the core scripting API hasn't really changed.
  17. The project is now called Mayavi and not MayaVi as before. This avoids unnecessary confusion on how to pronounce the name and avoids any comparison with either Maya or Vi.
  18. ETS itself is reorganized into a much smaller set of packages unlike the 40 odd packages in the ETS-2.x series. This makes dependency handling, packaging and installing much easier.

Mlab:

  1. The enthought.mayavi.mlab.pipeline is complete and can be used to fully script mayavi.
  2. The mlab API has changed to be more consistent with the naming style used in ETS, for example isosurface has become iso_surface, extractedges becomes extract_edges etc.
  3. Added a show() function and decorator to allow users to easily create standalone scripts.
  4. mlab.pipeline.open lets you open any supported data.
  5. The mlab API can now take either engine or figure keyword arguments. This allows to avoid the use of the global sate set in the mlab engine. Mlab also now exposes a set_engine function.
  6. It is easy to change visualized data using the .mlab_source attribute on objects created from mlab. This makes it very efficient and easy to create animations from mlab. See here for more details.
  7. Mlab by default uses a MayaviScene that features a convenient Mayavi icon which brings up the engine view using which you can edit the pipeline from the UI (using the toolbar or right clicks). This gives mlab the full power of mayavi.
Apart from these significant feature additions there have been the usual round of bug fixes and new bugs introduced.

As you can see this is a very significant release that marks a very important phase for mayavi2. All the additions made at the sprint went into this release.

Currently it is probably easiest to install mayavi via either enthought's EPD or Python(x,y). Gaël has made available Ubuntu packages that are available at a link he mentions here. Dave Peterson has been making all of ETS available from PyPI, however you'll have to get all the dependencies installed (numpy, VTK, wxPython or Qt4). The best place to look for installation instructions is here, https://svn.enthought.com/enthought/wiki/Install

Enjoy.