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.

Uniform deviates on the surface of a sphere

Here is a simple example of how to generate points randomly on the surface of the sphere that are uniformly distributed. The reason I'm posting it here is that it shows how convenient mayavi's mlab is for this sort of thing.

Anyone who has studied pseudo random numbers should know that the following will not generate points uniformly on the sphere. Here is a demonstration of the fact:

$ ipython -wthread
In [1]: from numpy import *
In [2]: from enthought.mayavi import mlab
In [3]: p = random.random(10000)*2*pi
In [4]: t = random.random(10000)*pi
In [5]: x, y, z = sin(t)*cos(p), sin(t)*sin(p), cos(t)
In [6]: g = mlab.points3d(x, y, z, z,  mode='point')
In [7]: g.actor.property.point_size = 2

The point size is made larger so you can see the points more clearly. Here is the picture it produces:


Notice the clustering at the poles. OTOH, the following works nicely,

In [3]: p = random.random(10000)*2*pi
In [10]: mlab.clf()
In [12]: t = arccos(random.random(10000)*2.0 - 1.0)
In [13]: x, y, z = sin(t)*cos(p), sin(t)*sin(p), cos(t)
In [14]: g1 = mlab.points3d(x, y, z, z,  mode='point')
In [15]: g1.actor.property.point_size = 2

This is the picture this one produces:



As you can see the transformation has worked and generates what appears to be a uniform distribution on the surface of a sphere. Visual testing isn't good enough uniformness but that isn't the point here.