Numpy tutorial

The Game of Life

The way of Python


In [1]:
Z = [[0,0,0,0,0,0],
     [0,0,0,1,0,0],
     [0,1,0,1,0,0],
     [0,0,1,1,0,0],
     [0,0,0,0,0,0],
     [0,0,0,0,0,0]]

First a neighbour counting function:


In [2]:
def compute_neighbours(Z):
    shape = len(Z), len(Z[0])
    N = [[0,]*(shape[0]) for i in range(shape[1])]
    for x in range(1,shape[0]-1):
        for y in range(1,shape[1]-1):
            N[x][y] = Z[x-1][y-1] + Z[x][y-1] + Z[x+1][y-1] \
                    + Z[x-1][y]               + Z[x+1][y]   \
                    + Z[x-1][y+1] + Z[x][y+1] + Z[x+1][y+1]
    return N

Then a function to perform this function on every cell:


In [3]:
def iterate(Z):
    N = compute_neighbours(Z)
    shape = len(Z), len(Z[0])
    for x in range(1,shape[0]-1):
        for y in range(1,shape[1]-1):
            if Z[x][y] == 1 and (N[x][y] < 2 or N[x][y] > 3):
                Z[x][y] = 0
            elif Z[x][y] == 0 and N[x][y] == 3:
                Z[x][y] = 1
    return Z

A function to show just the part of the board we are interested in:


In [4]:
def show(Z):
    for l in Z[1:-1]: print l[1:-1]
    print

Before any iterations:


In [5]:
show(Z)


[0, 0, 1, 0]
[1, 0, 1, 0]
[0, 1, 1, 0]
[0, 0, 0, 0]

After four iterations:


In [6]:
for i in range(4): iterate(Z)
show(Z)


[0, 0, 0, 0]
[0, 0, 0, 1]
[0, 1, 0, 1]
[0, 0, 1, 1]

The way of Numpy


In [7]:
import numpy as np
Z = np.array([[0,0,0,0,0,0],
              [0,0,0,1,0,0],
              [0,1,0,1,0,0],
              [0,0,1,1,0,0],
              [0,0,0,0,0,0],
              [0,0,0,0,0,0]])

Check the data type of the array:


In [8]:
print Z.dtype


int64

Check the shape of the array:


In [9]:
print Z.shape


(6, 6)

Accessing individual cells:

$\verb|Z[row_index, col_index]|$


In [10]:
print Z[0,5]


0

Slice notation:


In [11]:
print Z[1:5,1:5]


[[0 0 1 0]
 [1 0 1 0]
 [0 1 1 0]
 [0 0 0 0]]

Changing subparts impacts the whole array:


In [12]:
A = Z[1:5,1:5]
A[0,0] = 9
print A


[[9 0 1 0]
 [1 0 1 0]
 [0 1 1 0]
 [0 0 0 0]]

In [13]:
print Z


[[0 0 0 0 0 0]
 [0 9 0 1 0 0]
 [0 1 0 1 0 0]
 [0 0 1 1 0 0]
 [0 0 0 0 0 0]
 [0 0 0 0 0 0]]

We can check the influence of subparts by checking if the subpart is part of a larger array:


In [14]:
print Z.base is None


True

In [15]:
print A.base is None


False

In [16]:
print A.base is Z


True

Counting neighbours

Unlike in the Python case, we're going to act on the entire array (vectorisation).

We can manipulate Z as though it is a regular scalar:


In [17]:
print 1+(2*Z+3)


[[ 4  4  4  4  4  4]
 [ 4 22  4  6  4  4]
 [ 4  6  4  6  4  4]
 [ 4  4  6  6  4  4]
 [ 4  4  4  4  4  4]
 [ 4  4  4  4  4  4]]

Using vectorize computation:


In [18]:
Z = np.array([[0,0,0,0,0,0],
              [0,0,0,1,0,0],
              [0,1,0,1,0,0],
              [0,0,1,1,0,0],
              [0,0,0,0,0,0],
              [0,0,0,0,0,0]])

In [19]:
N = np.zeros(Z.shape, dtype=int)
N[1:-1,1:-1] += (Z[ :-2, :-2] + Z[ :-2,1:-1] + Z[ :-2,2:] + 
                 Z[1:-1, :-2]                + Z[1:-1,2:] + 
                 Z[2:  , :-2] + Z[2:  ,1:-1] + Z[2:  ,2:])

Iterate


In [20]:
def iterate(Z):
    # Iterate the game of life : naive version
    # Count neighbours
    N = np.zeros(Z.shape, dtype=int)
    N[1:-1,1:-1] += (Z[0:-2,0:-2] + Z[0:-2,1:-1] + Z[0:-2,2:] + 
                     Z[1:-1,0:-2]                + Z[1:-1,2:] + 
                     Z[2:  ,0:-2] + Z[2:  ,1:-1] + Z[2:  ,2:])
    N_ = N.ravel()
    Z_ = Z.ravel()
    
    # Apply rules
    R1 = np.argwhere( (Z_==1) & (N_ < 2) )
    R2 = np.argwhere( (Z_==1) & (N_ > 3) )
    R3 = np.argwhere( (Z_==1) & ((N_==2) | (N_==3)))
    R4 = np.argwhere( (Z_==0) & (N_==3) )
    
    # Set new values
    Z_[R1] = 0
    Z_[R2] = 0
    Z_[R3] = Z_[R3]
    Z_[R4] = 1
    
    # Make sure borders stay at null
    Z[0,:] = Z[-1,:] = Z[:,0] = Z[:,-1] = 0
    
    return Z

In [21]:
def iterate_2(Z):
    # Iterate the game of life : version with broadcasting
    # Count neighbours
    N = (Z[0:-2,0:-2] + Z[0:-2,1:-1] + Z[0:-2,2:] + 
         Z[1:-1,0:-2]                + Z[1:-1,2:] + 
         Z[2:  ,0:-2] + Z[2:  ,1:-1] + Z[2:  ,2:])
    
    # Apply rules
    birth = (N==3) & (Z[1:-1,1:-1]==0)
    survive = ((N==2) | (N==3)) & (Z[1:-1,1:-1]==1)
    Z[...] = 0
    Z[1:-1,1:-1][birth|survive] = 1
    return Z

Testing:


In [22]:
print Z


[[0 0 0 0 0 0]
 [0 0 0 1 0 0]
 [0 1 0 1 0 0]
 [0 0 1 1 0 0]
 [0 0 0 0 0 0]
 [0 0 0 0 0 0]]

In [23]:
for i in range(4): iterate_2(Z)
print Z


[[0 0 0 0 0 0]
 [0 0 0 0 0 0]
 [0 0 0 0 1 0]
 [0 0 1 0 1 0]
 [0 0 0 1 1 0]
 [0 0 0 0 0 0]]

Getting bigger


In [24]:
Z = np.random.randint(0,2,(256,512))

In [25]:
for i in range(100): iterate(Z)

In [26]:
import matplotlib.pyplot as plt
% matplotlib inline

size = np.array(Z.shape)
dpi = 72.0
figsize = size[1]/float(dpi),size[0]/float(dpi)
fig = plt.figure(figsize=figsize, dpi=dpi, facecolor="white")
fig.add_axes([0.0, 0.0, 1.0, 1.0], frameon=True)
plt.imshow(Z,interpolation='nearest', cmap=plt.cm.gray_r)
plt.xticks([]), plt.yticks([])
plt.show()



In [27]:
def game_of_life(Z,iterations,step_size):
    
    ### NEEDS OPTIMISING ###
    
    # Display before calculations
    size = np.array(Z.shape)
    dpi = 72.0
    figsize = size[1]/float(dpi),size[0]/float(dpi)
    fig = plt.figure(figsize=figsize, dpi=dpi, facecolor="white")
    fig.add_axes([0.0, 0.0, 1.0, 1.0], frameon=True)
    plt.imshow(Z,interpolation='nearest', cmap=plt.cm.gray_r)
    plt.xticks([]), plt.yticks([])
    plt.show()
    
    # Display for each step
    for step in range(iterations/step_size):
    
        # Iterate
        for i in range(step_size): iterate(Z)

        # Display after calculations
        size = np.array(Z.shape)
        dpi = 72.0
        figsize = size[1]/float(dpi),size[0]/float(dpi)
        fig = plt.figure(figsize=figsize, dpi=dpi, facecolor="white")
        fig.add_axes([0.0, 0.0, 1.0, 1.0], frameon=True)
        plt.imshow(Z,interpolation='nearest', cmap=plt.cm.gray_r)
        plt.xticks([]), plt.yticks([])
        plt.show()
    
    return None

In [28]:
Z = np.random.randint(0,2,(128,256))

In [29]:
game_of_life(Z,500,50)