`numeric`

`numarray`

David Mertz, Ph.D.

Idempotentate, Gnosis Software, Inc.

September, 2003

Numerical Python (often called NumPy) is a widely used extension library for fast operations on fixed-type arrays, of any dimensionality, in Python. Since the underlying code is well-optimized C, any speed limitations of Python's interpreter usually go away when major operations are performed in NumPy calls. As successful as NumPy has been, its developers have decided to supercede NumPy with a new module called Numarray that is mostly, but not quite entirely, compatible with NumPy. This installment looks both at the general features of NumPy, and at the specific improvements forthcoming with Numarray.

The first thing to know about the Numerical Python package is that it
doesn't let you do anything you cannot already do with standard
Python. What it does, is let you do many of the same things *a heck
of a lot faster*. Actually, that's not quite all; there are a number
of operations on arrays that are much more elegant to express in
* numeric* or

`numarray`

In essence, Numerical Python just gives you a new datatype, the array. Numerical arrays contain only elements of a homogeneous datatype, in contrast to to lists, tuples, and dictionaries which may contain heterogeneous elements. The other useful thing about Numerical arrays is that they may be multi-dimensional--but the dimensionality of arrays is a bit different that the simple nestability of lists. Numerical Python draws on the experience of programmers--particularly those from scientific backgrounds--who have abstracted the best features of arrays languages like APL, FORTRAN, MATLAB, and S, and created arrays whose shape and dimensionality is easily changeable. Back to this soon.

Operations on arrays in Numerical Python are performed *elementwise*.
Even though two dimensional arrays are similar to matrices from linear
algebra, operations (like multiply) have nothing to do with the
operations in linear algebra (like matrix multiplication).

Let us look at a quick example of all these points. In plain Python, you could create a "2-D list" with, e.g.:

>>> pyarr = [[1,2,3], ... [4,5,6], ... [7,8,9]] >>> print pyarr [[1, 2, 3], [4, 5, 6], [7, 8, 9]] >>> pyarr[1][1] = 0 >>> print pyarr [[1, 2, 3], [4, 0, 6], [7, 8, 9]]

Good enough, but about all you can do with this structure is set and retrieve elements by individual (multi-dimesional) index. In contrast, a Numerical array is a bit more flexible:

>>> from numarray import * >>> numarr = array(pyarr) >>> print numarr [[1 2 3] [4 0 6] [7 8 9]]

Not much has changed, but what about operations?

>>> numarr2 = numarr * 2 >>> print numarr2 [[ 2 4 6] [ 8 0 12] [14 16 18]] >>> print numarr2 + numarr [[ 3 6 9] [12 0 18] [21 24 27]]

And changing the shape:

>>> numarr2.shape = (9,) >>> print numarr2 [ 2 4 6 8 0 12 14 16 18]

For the most part, the newer package * numarray* is API compatible
with the earlier

`numeric`

`numeric`

`numeric`

`numarray`

`numarray`

`numeric`

Some of the improvements for * numarray* are:

(1) Element types are arranged in a hierarchical class structure to
support isinstance() checks. * numeric* used only character type codes
when specifying data types (but initializers in

`numarray`

(2) The type coercion rules have changed to preserve types in arrays (more often), rather than converting to the type of a Python scalar.

(3) Added array attributes (not just getters and setters).

(4) More flexible exception handling.

New users probably do not to worry about any of these changes; it is
best to start out with * numarray* at this point, in most cases.

Let us try to get a sense of just how much faster operations in Numerical Python can be than standard Python. As a toy task, we will create a sequence of numbers, then double each one. First, a few variants on standard Python approaches:

def timer(fun, n, comment=""): from time import clock start = clock() print comment, len(fun(n)), "elements", print "in %.2f seconds" % (clock()-start) def double1(n): return map(lambda n: 2*n, xrange(n)) timer(double1, 5000000, "Running map() on xrange iterator:") def double2(n): return [2*n for n in xrange(n)] timer(double2, 5000000, "Running listcomp on xrange iter: ") def double3(n): double = [] for n in xrange(n): double.append(2*n) return double timer(double3, 5000000, "Building new list from iterator: ")

We can check whether any speed differences exist between `map()`

, list
comprehensions, and traditional loops. What about standard module
* array*, which also requires homogeneous element types, perhaps it
will run faster:

import array def double4(n): return [2*n for n in array.array('i',range(n))] timer(double4, 5000000, "Running listcomp on array.array: ")

Finally, let us look at how * numarray* does. As an extra check, we
should see whether we still benefit if it is necessary to convert an
array back into a standard list:

from numarray import * def double5(n): return 2*arange(n) timer(double5, 5000000, "Numarray scalar multiplication: ") def double6(n): return (2*arange(n)).tolist() timer(double6, 5000000, "Numarray mult, returning list: ")

Now running it:

$ python2.3 timing.py Running map() on xrange iterator: 5000000 elements in 13.61 seconds Running listcomp on xrange iter: 5000000 elements in 16.46 seconds Building new list from iterator: 5000000 elements in 20.13 seconds Running listcomp on array.array: 5000000 elements in 25.58 seconds Numarray scalar multiplication: 5000000 elements in 0.61 seconds Numarray mult, returning list: 5000000 elements in 3.70 seconds

The minor differences between techniques for working with lists is
perhaps worth noticing, as is the misstep in trying standard module
* array*. But the general moral is that

`numarray`

Do not put too much weight in such a simple comparison, but the speedup is probably about typical. For big scientific calculations, moving from months to days--or days to hours--become quite worthwhile.

The typical use case for Numerical Python is in scientific modeling,
or perhaps in related areas like graphic manipulations and
convolutions or signal processing. To illustrate a number of features
of * numarray*, I will use a moderately realistic problem. Suppose you
have a 3D physical space over which a quantity varies. Abstractly, any
parametric space, with any number of dimensions, is equally applicable
to

`numarray`

While this point is perhaps obvious, it is worth mentioning
explicitly. My below examples work with smallish arrays because it is
easier to illustrate with them. But * numarray* remains fast for
arrays with millions of elements, rather than merely dozens; the
former are probably more common for real scientific models.

As a start, let us create the "room". There are a number of ways to
do this--the most general is with the callable `array()`

. Using that,
we can generate Numerical arrays with a variety of initialization
parameters, including initial data from any Python sequences. But for
our room, the convenience function `zeros()`

produces a nice--evenly
chilly--room:

>>> from numarray import * >>> room = zeros((4,3,5),Float) >>> print room [[[ 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0.]] [[ 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0.]] [[ 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0.]] [[ 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0.]]]

Each 2-D "matrix" printed represents on horizontal slice of the 3-D room--we will call it top-to-bottom. First thing, let us raise the whole room temperature to a comfortable 70 degree Farenheit:

>>> room += 70 >>> print room [[[ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.]] [[ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.]] [[ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.]] [[ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.]]]

There is a very important difference between Numerical arrays and Python lists to keep in mind for the rest of our operations. When you take slices of arrays--and we will see that slicing is quite flexible and powerful for multiple dimensions--what you get is not a copy, but a "view." There are multiple ways to point to the same data. Let us take a look. Suppose our room has a draft that comes in and cools the floor by four degrees:

>>> floor = room[3] >>> floor -= 4 >>> print room [[[ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.]] [[ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.]] [[ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.]] [[ 66. 66. 66. 66. 66.] [ 66. 66. 66. 66. 66.] [ 66. 66. 66. 66. 66.]]]

In contrast, the fireplace on the north wall heats each position next to it by 8 degrees, while itself being an ambient 90 degrees:

>>> north = room[:,0] >>> near_fireplace = north[2:4,2:5] >>> near_fireplace += 8 >>> north[3,2] = 90 # the fireplace cell itself >>> print room [[[ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.]] [[ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.]] [[ 70. 78. 78. 78. 70.] [ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.]] [[ 66. 74. 90. 74. 66.] [ 66. 66. 66. 66. 66.] [ 66. 66. 66. 66. 66.]]]

Here we get into some tricker indexing, where slices can be indicated along multiple dimensions. These views are nice to keep around. For example, you might wonder what the current temperatures are over the north wall:

>>> print north [[ 70. 70. 70. 70. 70.] [ 70. 70. 70. 70. 70.] [ 70. 78. 78. 78. 70.] [ 66. 74. 90. 74. 66.]]

This introductory look will only a small part of all the convenience
functions and array methods/attributes that are available in
* numarray*. My hope is only to give reader a general feel; the

`numarray`

Now that our room has been subject to various temperature fluctuation, we might want to judge its overall state. For example, what is the current average temperature in the room:

>>> add.reduce(room.flat)/len(room.flat) 70.066666666666663

This merits a bit of explanation. All the operations you can perform
on arrays have corresponding *Universal Functions* (ufuncs). So where
we had done something like `floor-= 4`

in the above code, we could
have equally well used `subtract(floor,4,floor)`

. With three
arguments, to `substract()`

, the operation is performed in-place. You
could also create a copy of `floor`

using `floor=subtract(floor,4)`

;
but this is probably not what you want, since the change will be in a
new array, not a subsection of `room`

.

However, ufuncs are not merely functions. They are callable objects
that also have some methods of their own: `.reduce()`

is probably the
most useful of these. It works just like the Python built-in function
`reduce()`

, with each operation being the underlying ufunc (the
methods are much faster though when applied to Numerical arrays). In
other words `add.reduce()`

means `sum()`

and `multiply.reduce()`

means
`product()`

.(those shortcut names are defined too).

Before you can sum the room cell temperatures, you need to take a flat, 1-D, view of the data. Otherwise, you add along the primary dimension, producing an new array of reduced dimension, e.g.:

>>> add.reduce(room) array([[ 276., 292., 308., 292., 276.], [ 276., 276., 276., 276., 276.], [ 276., 276., 276., 276., 276.]])

Such dimensional summary can be useful, but it is not what we want in this case.

Since we are modeling a physical system, let us add a bit more realism. Rooms have air micro-currents in them that mix up temperature. We might model that by supposing that every once in a while, a cell adjusts to the temperatures around it:

>>> def equalize(room): ... z,y,x = map(randint, (1,1,1), room.shape) ... zmin,ymin,xmin = maximum([z-2,y-2,x-2],[0,0,0]).tolist() ... zmax,ymax,xmax = [z+1,y+1,x+1] ... region = room[zmin:zmax,ymin:ymax,xmin:xmax].copy() ... room[z-1,y-1,x-1] = sum(region.flat)/len(region.flat) ... return room

This model is a bit unrealistic, of course: cells pick up on the
temperature around them, but do not impart a differential back to
their neighbors. Still, let us look at what it does. First we pick a
random cell--or actually we pick one more than the cell itself in each
dimension, since `.shape`

gives us the length rather than the maximum
index. The value `zmin`

and friends make sure we only index down to
zero, not in negative ranges; `zmax`

and friends are not really
needed, since indexing past the size of an array dimension is just
treated as the maximum value (much as with Python lists).

Next thing, we need to define the region of adjacent cells. Given our
small room dimensions, picking a cell at a room surface, edge, or
corner is rather common--the `region`

of a cell is likely to be
smaller than the maximum 27-elements of a 3x3x3 subsection. That is
not a problem, we just need to be sure to use the right denominator
for our average. This new average temperature is assigned to the
randomly chosen cell.

You can perform as many equalization events as you wish to in your
model. Each call just adjusts one cell. Asymptotically, a lot of
calls will lead to equalization of the temperature in parts of the
room. Even though arrays are mutable, the function `equalize()`

returns its array. This is useful if you only want to equalize a
*copy* of the model:

>>> print equalize(room.copy()) [[[ 70. 70. 70. 70. 70. ] [ 70. 70. 70. 70. 70. ] [ 70. 70. 70. 70. 70. ]] [[ 70. 70. 71.333333 70. 70. ] [ 70. 70. 70. 70. 70. ] [ 70. 70. 70. 70. 70. ]] [[ 70. 78. 78. 78. 70. ] [ 70. 70. 70. 70. 70. ] [ 70. 70. 70. 70. 70. ]] [[ 66. 74. 90. 74. 66. ] [ 66. 66. 66. 66. 66. ] [ 66. 66. 66. 68. 66. ]]]

This introduction touched on only a bit of what * numarry* can do. The
full capabilities are quite fancy. For example, you can populate
arrays based on populating functions, which is likely to be useful for
physical models. You can specify subsections of arrays not only with
slices, but also with index arrays--this not only lets you operate on
non-contiguous segments of an array, it also lets you--via the

`take()`

function, redimension and reform arrays in interesting ways.
Moreover, for the most part, I have presented only operations between
scalars and arrays; you can also perform operations between arrays,
including ones of different dimensionality. There is a lot to it, but
the API manages to remain fairly intuitive throughout.
I encourage readers to install * numarray* and/or

`numeric`

The Numerical Python homepage is at the below URL. This has good
links to documentation on both * numarray* and

`numeric`

http://www.pfdubois.com/numpy/

The project itself is hosted by Sourceforge, and you can find the Numerical Python project page at:

http://sourceforge.net/projects/numpy

While David Mertz also likes laziness and impatience, this installment is about hubris. David may be reached at mertz@gnosis.cx; his life pored over athttp://gnosis.cx/publish/. Suggestions and recommendations on this, past, or future, columns are welcomed. Check out his book, Text Processing in Python at <http://gnosis.cx/TPiP>.