Besides its obvious scientific uses, NumPy can also be used as an efficient
multi-dimensional container of generic data. Arbitrary data-types can be
defined and this allows NumPy to seamlessly and speedily integrate with a wide
variety of projects. We are going to explore numpy through a simple example,
implementing the Game of Life.

Numpy is slanted toward scientific computing and we'll consider in this section
the game of life by
John Conway which is one of the earliest example of cellular automata (see
figure below). Those cellular automaton can be conveniently considered as array
of cells that are connected together through the notion of neighbours. We'll
show in the following sections implementation of this game using pure python
and numpy in order to illustrate main differences with python and numpy.

The Game of Life, also known simply as Life, is a cellular automaton devised
by the British mathematician John Horton Conway in 1970. It is the
best-known example of a cellular automaton. The "game" is actually a
zero-player game, meaning that its evolution is determined by its initial
state, needing no input from human players. One interacts with the Game of
Life by creating an initial configuration and observing how it evolves.

The universe of the Game of Life is an infinite two-dimensional orthogonal grid
of square cells, each of which is in one of two possible states, live or
dead. Every cell interacts with its eight neighbours, which are the cells that
are directly horizontally, vertically, or diagonally adjacent. At each step in
time, the following transitions occur:

Any live cell with fewer than two live neighbours dies, as if by needs caused
by underpopulation.

Any live cell with more than three live neighbours dies, as if by
overcrowding.

Any live cell with two or three live neighbours lives, unchanged, to the next
generation.

Any dead cell with exactly three live neighbours becomes a live cell.

The initial pattern constitutes the 'seed' of the system. The first generation
is created by applying the above rules simultaneously to every cell in the seed
– births and deaths happen simultaneously, and the discrete moment at which
this happens is sometimes called a tick. (In other words, each generation is a
pure function of the one before.) The rules continue to be applied repeatedly
to create further generations.

We'll first use a very simple setup and more precisely, we'll use the glider pattern that is known to
move one step diagonally in 4 iterations as illustrated below:

t = 0

t = 1

t = 2

t = 3

t = 4

This property will help us debug our scripts.

The way of python

Note

We could have used the more efficient python array interface but people may be more
familiar with the list object.

In pure python, we can code the Game of Life using a list of lists representing
the board where cells are supposed to evolve:

Note that we did not specify the data type of the array and thus, numpy will
choose one for us. Since all elements are integers, numpy will then choose an
integer data type. This can be easily checked using:

>>> print Z.dtype
int64

We can also check the shape of the array to make sure it is 6x6:

>>> print Z.shape
(6, 6)

Each element of Z can be accessed using a row and a column
index (in that order):

>>> print Z[0,5]
0

Note

This element access is actually called indexing and this
is very powerful tool for vectorized computation.

But even better, we can also access a subpart of the array using the slice
notation:

>>> print Z[1:5,1:5]
[[0 0 1 0]
[1 0 1 0]
[0 1 1 0]
[0 0 0 0]]

In the example above, we actually extract a subpart of Z ranging from rows 1 to
5 and columns 1 to 5. It is important to understand at this point that this is
really a subpart of Z in the sense that any change to this subpart will
have immediate impact on Z:

We set the value of A[0,0] to 9 and we see immediate change in Z[1,1]
because A[0,0] actually corresponds to Z[1,1]. This may seem trivial
with such simple arrays, but things can become much more complex (we'll see
that later). If in doubt, you can check easily if an array is part of another
one:

>>> print Z.base is None
True
>>> print A.base is Z
True

Counting neighbours

Note

It is not always possible to vectorize computations and it requires
generally some experience. You'll acquire this experience by using numpy (of
course) but also by asking questions on the mailing list

We now need a function to count the neighbours. We could do it the same way as
for the python version, but this would make things very slow because of the
nested loops. We would prefer to act on the whole array whenever possible, this
is called vectorization.

Ok, let's start then...

First, you need to know that you can manipulate Zas if (and only as
if) it was a regular scalar:

If you look carefully at the output, you may realize that the ouptut
corresponds to the formula above applied individually to each element. Said
differently, we have (1+(2*Z+3))[i,j] == (1+(2*Z[i,j]+3)) for any i,j.

Ok, so far, so good. Now what happens if we add Z with one of its subpart,
let's say Z[-1:1,-1:1] ?

This raises a Value Error, but more interestingly, numpy complains about
the impossibility of broadcasting the two arrays together. Broadcasting is
a very powerful feature of numpy and most of the time, it saves you a lot of
hassle. Let's consider for example the following code:

How can a matrix and a scalar be added together ? Well, they can't. But numpy
is smart enough to guess that you actually want to add 1 to each of the element
of Z. This concept of broadcasting is quite powerful and it will take you
some time before masterizing it fully (if even possible).

However, in the present case (counting neighbours if you remember), we won't
use broadcasting (uh ?). But we'll use vectorize computation using the
following code:

What we actually did with the above code is to add all the darker blue squares
together. Since they have been chosen carefully, the result will be exactly
what we expected. If you want to convince yourself, consider a cell in the
lighter blue area of the central sub-figure and check what will the result for
a given cell.

Iterate

Note

Note the use of the ravel function that flatten an array. This is necessary since the argwhere function returns flattened indices.

In a first approach, we can write the iterate function using the argwhere
method that will give us the indices where a given condition is True.

Even if this first version does not use nested loops, it is far from optimal
because of the use of the 4 argwhere calls that may be quite slow. We can
instead take advantages of numpy features the following way.

If you look at the birth and survive lines, you'll see that these two
variables are indeed arrays. The right-hand side of these two expressions are
in fact logical expressions that will result in boolean arrays (just print them
to check). We then set all Z values to 0 (all cells become dead) and we use
the birth and survive arrays to conditionally set Z values
to 1. And we're done ! Let's test this:

Getting bigger

While numpy works perfectly with very small arrays, you'll really benefit from
numpy power with big to very big arrays. So let us reconsider the game of life
with a bigger array. First, we won't initalize the array by hand but initalize
it randomly:

A step further

Reaction and diffusion of chemical species can produce a variety of patterns,
reminiscent of those often seen in nature. The Gray Scott equations model such
a reaction. For more information on this chemical system see the article
Complex Patterns in a Simple System, John E. Pearson, Science, Volume 261,
9 July 1993.

Let's consider two chemical species U and V with respective
concentrations u and v and diffusion rates ru and
rv. V is converted into P with a rate of conversion
k. f represents the rate of the process that feeds U
and drains U, V and P. We can now write:

Chemical reaction

Equations

U + 2V⟶3V

V⟶P

(∂u)/(∂t) = ru∇2u − uv2 + f(1 − u)

(∂v)/(∂t) = rv∇2v + uv2 − (f + k)v

Examples

(click figure to see movie)

Obviously, you may think we need two arrays, one for U and for V. But
since U and V are tighly linked, it may be indeed better to use a
single array. Numpy allows to do that with the notion of structured array:

The size of the array is (n+2,n+2) since we need the borders when computing the
neighbours. However, we'll compute differential equation only in the center
part, so we can already creates some useful views of this array:

>>> U,V = Z['U'], Z['V']
>>> u,v = U[1:-1,1:-1], V[1:-1,1:-1]

Next, we need to compute the Laplacian and we'll use a discrete approximation
obtained via the finite difference method using the same vectorization as for the Game of Life:

Other Tutorials

The SciPy Lecture notes offers a teaching material on the scientific Python
ecosystem as well as quick introduction to central tools and techniques. The
different chapters each correspond to a 1 to 2 hours course with increasing
level of expertise, from beginner to expert.