Game of Life in NumPy (Preliminaries)

You might be familiar with Conway’s Game of Life. It is a simple game based on the idea of modeling a dynamic system progressing over time. We progress from discrete time points \(t_i \rightarrow t_i+1\) and apply very simple rules at each time point. The game progresses on a square grid of discrete points. Each cell on the grid has eight neighbors (the four cardinal directions, plus diagonals). The progression to the next time interval is given by the following:

  • If I’m alive and I have fewer than two alive neighors, I die of loneliness.
  • If I’m alive and I have two or three alive neighbors, I live.
  • If I’m alive and I have more than three alive neighbors, I die of starvation.
  • If I’m dead and I have three (exactly) live neighbors, I become alive by spontaneous combustion.

Typically, these rules are applied by loops and conditionals. I want to show-off a lesser known corner of NumPy by looking at as_strided from NumPy’s dark corner of stride_tricks. We’ll get to know as_strided today and, in a future post, we’ll play the game of life using stride tricks – woot!

Introducing Array Strides

Before we get to anything fancy, we should review the role of strides in NumPy. Memory for all the elements in a NumPy array are layed out as one long contiguous block of (linear) memory. We use clever tricks to give this linear block of memory the structure we need to deal with rows and columns (or more dimensions). Consider the following:

In [1]:
import numpy as np
In [2]:
arr_1d = np.arange(12)
print arr_1d

[ 0  1  2  3  4  5  6  7  8  9 10 11]
In [3]:
arr_2d = np.arange(12).reshape(3,4)
print arr_2d

[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]

Both of these arrays have exactly the same elements. Perhaps surprisingly, the elements can be – and are! – stored in the exactly same flat sequence of values. The following gives us a look at that.

In [4]:
print arr_1d.flatten()
print arr_2d.flatten()

[ 0  1  2  3  4  5  6  7  8  9 10 11]
[ 0  1  2  3  4  5  6  7  8  9 10 11]

It turns out that the real difference between the two arrays is not how they are stored in memory, but what we do with that flat block of data. And that is determined by the shape of the array.

In [5]:
print arr_1d.shape
print arr_2d.shape

(12,)
(3, 4)

Let’s look at the size of the 2-D array. The value of 4 tells us that the inner most dimension (it’s the last thing in the shape) has four elements in it. This has an interesting side-effect: it determines the step sized need for the outer dimension. To get from outer dimension 0 (the first visual row) to outer dimension 1 (the second visual row), we have to walk across the 4 elements of the inner dimension.

We can apply this to give a formula to go from an multi-dimensional index arr_2d[1,2] to the specific element we need in the flat data as it really exists in memory \([1,2]{\rightarrow}1*4+2\) More generally, that is \([row,col]{\rightarrow}(row*rowSize)+col\) or perhaps \([1,2]{\rightarrow}(row*numCols)+col\). The same sort of reasoning applies as we add more dimensions. Here, \(rowSize\) (or, equivalently, \(numCols\)) determines the "number of steps" or the stride size (in the flat array) for the outer dimension. We can look at these explicitly:

In [6]:
print arr_1d.strides
print arr_2d.strides

(8,)
(32, 8)

Whoa! Those aren’t at all what I led you to expect.

One other factoid: we’re looking at the lower levels (inner workings and implementation details) of the NumPy library. This means we have to start worrying about the machine representations of our data. Our arrays have dtype=np.int_ which is 64 bits per int on modern 64-bit CPUs (aka, np.int64). So, moving forward one data entry means moving forward 64 bits which is equivalent to 8 bytes. Dividing those values given by arr.strides by 8 gives us "element strides" (as opposed to the "memory strides" we’ll be working with) of (1,) for arr_1d and of (4,1) for arr_2d. We can rewrite \([1,2]{\rightarrow}1*4+2\) as \([1,2]{\rightarrow}1*\color{red}{4}+2*\color{red}{1}\) and the "coefficient" on the indices are now \(\color{red}{(4,1)}\) which is the same as our "element strides".

In the remainder of this post, we will be working with the proper "memory strides", not these slightly more convenient "element strides".

Introducing as_strided

In [7]:
from numpy.lib.stride_tricks import as_strided

We can use as_strided to customize the way we walk across an array. That is, we can specify strides other than those implied by the shape of the array. Because of the flexibility in specifying the strides, we also have to tell NumPy what the resulting output shape will be. Let’s revisit the 1-D array:

In [8]:
print arr_1d
shapeShifted = as_strided(arr_1d,            # array to monkey with
                          strides=(4*8,1*8), # custom step sizes
                          shape=(3,4))       # output shape

print shapeShifted

[ 0  1  2  3  4  5  6  7  8  9 10 11]
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]

To recap, we told NumPy to walk across arr_1d using steps of 1 element on the inner-most axis (walking across the visual rows) and steps of 4 elements to kick off a new row.

We can pull some other neat tricks. For example, look at what happens if we never "start a new row" (i.e., our outer stride is 0).

In [9]:
print arr_1d.itemsize # number of bytes per element (aka, 8)
print arr_1d.shape[0]
duplicated = as_strided(arr_1d, 
                        strides=(0, 1*arr_1d.itemsize), 
                        shape = (2,arr_1d.shape[0]))
print duplicated

8
12
[[ 0  1  2  3  4  5  6  7  8  9 10 11]
 [ 0  1  2  3  4  5  6  7  8  9 10 11]]

Really? This is exactly the sort of replication that happens in broadcasting. How did that output array get filled in? Consider how we filled in duplicated[0,5]. To fill in item [0,5] we have to take stride length for the outer dimension times the index of the outer dimension (\(0\cdot0=0\)). We then have to take the stride length for the inner dimension times the index of the inner dimension (\(8\cdot5=40\)). Adding the two gives us a flat offset of 40 (from \(0\cdot0+8\cdot5\). When we try to fill in the second row of duplicated, we do something similar. For duplicated[1,5] the values are: \(0\cdot1 + 8\cdot5=40\) (note, the 0 is the same in both computations because are stride for the outer dimension is 0).

We can do something analogous to the other axis and broadcast column-wise, instead of row-wise. We use both styles (outer/column and inner/row) if we want to compute an outer product.

In [10]:
arr1 = np.arange(1,5)
arr2 = np.arange(5,8)

outerProdShape = arr1.shape + arr2.shape

x = as_strided(arr1, 
               strides=(arr1.itemsize, 0), 
               shape=outerProdShape) # copy oPS[0] times (arr2 shape)
y = as_strided(arr2, 
               strides=(0, arr2.itemsize), 
               shape=outerProdShape) # copy oPS[1] times (arr1 shape)

print "x:"
print x

print "y:"
print y

print "x*y:"
print x * y

x:
[[1 1 1]
 [2 2 2]
 [3 3 3]
 [4 4 4]]
y:
[[5 6 7]
 [5 6 7]
 [5 6 7]
 [5 6 7]]
x*y:
[[ 5  6  7]
 [10 12 14]
 [15 18 21]
 [20 24 28]]

This is the same answer that we get using NumPy’s built-in outer function:

In [11]:
print np.outer(arr1, arr2)

[[ 5  6  7]
 [10 12 14]
 [15 18 21]
 [20 24 28]]

You might know that we can do some outer product like tasks with broadcasting by judicious use of a new axis (np.newaxis). Here’s an example:

In [12]:
print arr1[:,np.newaxis] * arr2[np.newaxis,:]

[[ 5  6  7]
 [10 12 14]
 [15 18 21]
 [20 24 28]]

Indeed, it gives the same answer!

Now, compare arr1[:, np.newaxis] with the custom strides we applied to arr1: strides took the value (arr1.itemsize, 0). On the inner dimension, the np.newaxis performs the same task as the 0 stride size. No kidding!

Neighborhoods

Let’s see how we can apply these tricks to pick out sub-arrays.

In [13]:
board = np.arange(16).reshape((4,4))
print board
print board.shape

[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]
 [12 13 14 15]]
(4, 4)

We want to pick out the four "complete" neighborhoods centered around 5, 6, 7, and 8. Let’s look at the neighborhood for 5. What is the shape of the result? 3×3. What are the strides? Well, to walk across a row is still just walking one element at a time and to get to the next row is still 4 elements at a time. These are the same as the strides in the original. The difference is we don’t take "everything", we just take a selection. Let’s see if that actually works:

In [14]:
neighbors = as_strided(board, shape=(3,3), strides=board.strides)
print neighbors

[[ 0  1  2]
 [ 4  5  6]
 [ 8  9 10]]

Ok, nice. Now, if we want all four neighborhoods, what is the output shape? We have several 3×3 results. How many? In this case, we have 2×2 of them (for each of the "center" cells). This gives a shape of (2,2,3,3) – the neighborhoods are the inner dimensions and the organization of the neighborhoods is the outer dimensions.

We know how to get one neighborhood out. That moves us across the inner two dimensions. Now we need to figure out how to get from one neighborhood to the next. Once we have the complete neighborhood for 5, we want to start the neighborhood for 6. It’s top-left neighbor is the number 1 – but this is just one element beyond where we started for the previous neighborhood (which means our stride is one element). And we can keep doing that the whole way across the board (we’ll stop based on the output shape).

After we’ve walked across the board (neighborhood by neighborhood), how do we get to the next row? There are two ways to look at it. The first is by comparing the neighborhood for 5 to the neighborhood for 9. They differ by one row. So, our stride is four elements (the number of elements per row). The other way to think about this is where we are at after we walked the neighborhoods across the previous row. We will have gone through all the neighborhoods (inner most dimensions) and the first of the outer two dimensions. We now need to move forward one row – which again, is four elements. I prefer the first explanation!

So, our strides (in terms of elements) end up being (4,0) within one neighborhood and (4,0) for progressing neighborhood to neighborhood. The total stride (element wise) is: (4,0,4,0). But, the component strides (our outer two dimensions) are the same as the strides of the board. This means that our neighborhood strides are board.strides + board.strides.

In [15]:
print board.strides + board.strides

neighborhoods = as_strided(board, 
                           shape=(2,2,3,3), 
                           strides=board.strides+board.strides)

print neighborhoods[0,0]
print neighborhoods[-1, -1]

(32, 8, 32, 8)
[[ 0  1  2]
 [ 4  5  6]
 [ 8  9 10]]
[[ 5  6  7]
 [ 9 10 11]
 [13 14 15]]

Generalizing

We can generalize this a bit to work for more dimensions. For example, we could compute the neighbors in 3-D space or in 4-D space.

The key points are that: (1) taking the neighborhood has the same number of dimensions as the input, (2) the number of neighborhoods is the length of the axis minus 2 (we don’t use the first or last elements on an axis as a neighborhood center), and (3) the new strides are just the original strides concatenated to the original strides.

In [ ]:
def grid_nD(arr):
    assert all(_len>2 for _len in arr.shape)
    
    nDims = len(arr.shape)
    newShape = [_len-2 for _len in arr.shape]
    newShape.extend([3] * nDims)
    print "shape:", newShape
    
    newStrides = arr.strides + arr.strides
    print "stride:", newStrides
    return as_strided(arr, shape=newShape, strides=newStrides)

Let’s test that on a 2D grid, since we’re pretty familiar with that case already.

In [16]:
x = np.arange(25).reshape((5,5))
print x
neighborhoods = grid_nD(x)
print neighborhoods[0,0]
print neighborhoods[-1,-1]

[[ 0  1  2  3  4]
 [ 5  6  7  8  9]
 [10 11 12 13 14]
 [15 16 17 18 19]
 [20 21 22 23 24]]
shape: [3, 3, 3, 3]
stride: (40, 8, 40, 8)
[[ 0  1  2]
 [ 5  6  7]
 [10 11 12]]
[[12 13 14]
 [17 18 19]
 [22 23 24]]

And now let’s see how it works on a 3D grid. Here, we’ll use 64 total grid points, but write it as \(4*4*4\) to emphasize the similarity between the size of the dimensions and this cube-shaped grid.

In [17]:
x = np.arange(4*4*4).reshape((4,4,4))
neighborhoods = grid_nD(x)
print neighborhoods[0,0,0]
print neighborhoods[-1,-1,-1]

shape: [2, 2, 2, 3, 3, 3]
stride: (128, 32, 8, 128, 32, 8)
[[[ 0  1  2]
  [ 4  5  6]
  [ 8  9 10]]

 [[16 17 18]
  [20 21 22]
  [24 25 26]]

 [[32 33 34]
  [36 37 38]
  [40 41 42]]]
[[[21 22 23]
  [25 26 27]
  [29 30 31]]

 [[37 38 39]
  [41 42 43]
  [45 46 47]]

 [[53 54 55]
  [57 58 59]
  [61 62 63]]]

References for this Material

I made varying use of the following material to put this together:

  • http://scipy-lectures.github.io/advanced/advanced_numpy/#indexing-scheme-strides
  • http://chintaksheth.wordpress.com/2013/07/31/numpy-the-tricks-of-the-trade-part-ii/
  • http://wiki.scipy.org/Cookbook/GameOfLifeStrides
  • http://www.rigtorp.se/2011/01/01/rolling-statistics-numpy.html

Additional Resources

You can grab a copy of this notebook.

Even better, you can view it using nbviewer.

License

Unless otherwise noted, the contents of this notebook are under the following license. The code in the notebook should be considered part of the text (i.e., licensed and treated as as follows).

Creative Commons License
DrsFenner.org Blog And Notebooks by Mark and Barbara Fenner is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.
Permissions beyond the scope of this license may be available at drsfenner.org/blog/about-and-contacts.