NumPy

The Numerical Python Extensions, also called NumPy or Numeric, turn Python into an ideal tool for experimental numerical and scientific computing (better than specialized programs like MatLab, Octave, RLab or SciLab). NumPy is useful for everybody who analyzes data with the help of a spreadsheet program like Microsoft Excel—it is not just for mathematicians and scientists who crunch lots of data.

NumPy defines a new data type, NumPy array, and a very complete set of operators and functions to manipulate NumPy arrays. All the functionality of NumPy can be obtained in pure Python, but NumPy gives you speed and elegance.

In the following, I assume that you have installed NumPy on your system. Doing so is not really difficult. There are binary packages for Windows, or source packages for all platforms. A source package is installed using distutils (see Chapter 26), by typing

root@calcifer:/home/boud# python setup_all.py install
    

Once numpy is installed, you can start Python (or open the Interpreter window in BlackAdder) and import the NumPy extension:

[packer@slow packer]$ python
Python 2.1.1 (#1, Aug 20 2001, 08:17:33) 
[GCC 2.95.3 19991030 (prerelease)] on linux2
Type "copyright", "credits" or "license" for more information.
>>> from Numeric import *
>>>
    

A NumPy array looks like a list and can be created from a list (in fact, from any sequency type: list, tuple or string).

Let's create and print a 1-dimensional NumPy array of Python floats:

>>> a = array([1.0, 4.0, 9.0, 16.0, 25.0, 36.0, 49.0, 64.0, 81.0, 100.0])
>>> print a
[   1.    4.    9.   16.   25.   36.   49.   64.   81.  100.]
>>>
    

This creates a 1-dimensional NumPy array. All elements in the list should have the same data type.

A 2-dimensional NumPy array is created from a list of sub-lists:

>>> b = array([[0.0, 1.0], [2.0, 3.0]])
>>> print b
[[ 0.  1.]
 [ 2.  3.]]
>>>
    

The sub-lists should have the same length, and all the elements in all the sub-lists should have the same data type.

You can show off with NumPy arrays of even higher dimensions (up to 40, by default). For example, a 3-dimensional NumPy array is created from a list of sub-lists of sub-sub-lists:

>>> c = array([[[0.0, 1.0], [2.0, 3.0]], [[4.0, 5.0], [6.0, 7.0]]])
>>> print c
[[[ 0.  1.]
  [ 2.  3.]]
 [[ 4.  5.]
  [ 6.  7.]]]
>>>
    

The sub-lists should have the same length, the sub-sub-lists should have the same length, and all elements of all sub-sub-lists should have the same data type.

In the following, I am going to compare the functionality of NumPy arrays and lists. Here is an easier method to create a NumPy array:

>>> ax = arange(0.0, 5.0, 0.5)
>>> print ax
[ 0.   0.5  1.   1.5  2.   2.5  3.   3.5  4.   4.5]
>>>
    

The function call arange(0.0, 5.0, 0.5) returns an array with elements ranging from 0.0 to 5.0 (non-inclusive) in steps of 0.5. Here is a similiar function to return a list with the same properties:

def lrange(start, stop, step):
    start, stop, step = float(start), float(stop), float(step)
    size = int(round((stop-start)/step))
    result = [start] * size
    for i in xrange(size):
        result[i] += i * step
    return result
    

After copying and pasting the function definition in your Python interpreter, do:

>>> lx = lrange(0.0, 5.0, 0.5)
>>> print lx
[0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5]
    

Why are NumPy arrays better than lists? The full answer is speed and elegance. To compare lists and NumPy arrays with respect to elegance, lets use a simple function:

def lorentzian(x):
    return 1.0/(1.0+(x-2.5)**2)
    

To calculate a list, ly, containing the function values for each element of ly, we can do:

>>> ly = [0.0]*len(lx)
>>> for i in range(len(lx)): ly[i] = lorentzian(lx[i])
...
    

Do you know that you can get rid of the loop? The following is more elegant and slightly faster:

>>> ly = map(lorentzian, lx)
    

NumPy arrays are even more elegant, and they allow:

>>> ay = lorentzian(ax)
    

Almost magic, isn't it? I wrote the function lorentzian(x) assuming that x is Python float. If you call lorentzian with a NumPy array as argument, it returns a NumPy array. This does not work with lists:

>>> ly = lorentzian(lx)
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
  File "<stdin>", line 2, in lorentzian
TypeError: unsupported operand type(s) for -
>>>
    

To compare speed, we create a list, xl, and a NumPy array, xa, with 100000 elements and use the profile module to time the statements yl = map(lorentzian, xl) and ya = lorentzian(xa):

>>> import profile
>>> xl = lrange(0, 10, 0.0001)
>>> xa = arange(0, 10, 0.0001)
>>> profile.run('yl = map(lorentzian, xl)')
         100002 function calls in 2.200 CPU seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
   100000    1.000    0.000    1.000    0.000 <stdin>:1(lorentzian)
        1    1.200    1.200    2.200    2.200 <string>:1(?)
        0    0.000             0.000          profile:0(profiler)
        1    0.000    0.000    2.200    2.200 profile:0(yl = map(lorentzian, xl))


>>> profile.run('ya = lorentzian(xa)')
         3 function calls in 0.090 CPU seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.090    0.090    0.090    0.090 <stdin>:1(lorentzian)
        1    0.000    0.000    0.090    0.090 <string>:1(?)
        0    0.000             0.000          profile:0(profiler)
        1    0.000    0.000    0.090    0.090 profile:0(ya = lorentzian(xa))


>>>
    

On my computer, the Numerical Python extensions are almost 25 times faster than pure Python!

There exists a scientific plotting program, SciGraphica (http://scigraphica.sourceforge.net), which allows you to manipulate your data in a spreadsheet. The underlying engine is a python interpreter with the NumPy. Each column in the spreadsheet is in reality a NumPy array. This clearly demostrates the power of this extension. If you want to know more about NumPy, you can consult the excellent documentation at http://www.pfdubois.com/numpy, the homepage of NumPy.