

In [1]:

    
%matplotlib inline



In [2]:

    
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import convolve2d, fftconvolve
from numpy.fft import rfft2, irfft2



In [3]:

    
from matplotlib import colors
import matplotlib.animation as animation
from IPython.display import HTML

Numbers, Arrays and Life

Cellular automata are discrete-space, discrete-time models of spatio-temporal processes. In these models, signal propagation typically involves updating the state of each cell in a grid based on some function of its local neighbors. We will use a toy example, but the techniques you will implement generalize to useful contexts, e.g. spatial statistics or image processing.

We will use John Conway's Game of Life to practice our manipulation of numerical arrays. In particular, we will see how to use the vectorization, views, convolution and Fourier transform capabilities of numpy to solve the same problem in different ways.

See Wikipedia for detailed rules. Here is a summary:

Any live cell with fewer than two live neighbors dies, as if caused by under-population.
Any live cell with two or three live neighbors lives on to the next generation.
Any live cell with more than three live neighbors dies, as if by over-population.
Any dead cell with exactly three live neighbors becomes a live cell, as if by reproduction.

There are two common ways to count neighbors for cells that are at the edge of the grid. With periodic boundary conditions, the grid wraps around (so it is a torus), that is, the coordinate -1 is mapped to the max_coord, and max_coord + 1 is mapped to 0. Alternatively, just consider any neighbor outside the grid to be equal to zero. For simplicity, we have chosen a life configuration where both strategies give the same results, so you can implement either strategy.

Initial grid configuration



In [4]:

    
gun = np.load('gosper_gun.npy')
init_grid = np.zeros((64, 64)).astype('int')
for y, x in gun:
    init_grid[8+x, 8+y] = 1

Utility function to show animation



In [5]:

    
def animate(i):
    """Function to display animation for iteration i."""
    global grid
    grid = step(grid)
    im.set_array(grid)
    return im,

1. Write a function to play one step in the Game of Life. The function must be named step and take a single argument grid which is a numpy array containing the initial configuration, and return a new numpy array containing the updated configuration. Use for loops to implement this first version of the step function.



In [6]:

    
# Your solution here



In [ ]:

    
fig = plt.figure(figsize=(5, 5))
grid = init_grid.copy()
im = plt.imshow(grid, animated=True, interpolation='nearest', cmap='gray')
plt.close(fig)

anim = animation.FuncAnimation(fig, animate, 
                               frames=60, interval=50, blit=True)
HTML(anim.to_html5_video())

2. Rewrite the step function using vectorization. That is, use eight different views of the grid to calculate the neighbor sum instead of double for loops.



In [9]:

    
# Your solution here



In [ ]:

    
fig = plt.figure(figsize=(5, 5))
grid = init_grid.copy()
im = plt.imshow(grid, animated=True, interpolation='nearest', cmap='gray')
plt.close(fig)

anim = animation.FuncAnimation(fig, animate, 
                               frames=60, interval=50, blit=True)
HTML(anim.to_html5_video())

3a. A discrete 2D convolution generates a weighted sum of a 2D grid, with the weights given by a 2D kernel. For example, the kernel

0 1 0
1 1 1
0 1 0

would result in summing the current location and the N, E, S, W neighbors.

Use the convolve2d function from scipy.signal (already imported for you in this notebook) to implement the 3rd version of the step function with a suitable kernel and with mode=same to preserve the grid size across iterations.



In [65]:

    
# Your solution here



In [ ]:

    
fig = plt.figure(figsize=(5, 5))
grid = init_grid.copy()
im = plt.imshow(grid, animated=True, interpolation='nearest', cmap='gray')
plt.close(fig)

anim = animation.FuncAnimation(fig, animate, 
                               frames=60, interval=50, blit=True)
HTML(anim.to_html5_video())

3b. One way to multiply two numbers that you are familiar with is to move to the log domain where addition is equivalent to multiplication. For example, instead of calculating $100 \times 1000$, we can calculate $\exp(\log(100) + \log(1000))$ to get the same result to roundoff error. In the same way, a convolution is equivalent to multiplication in the Fourier domain.

Implement the fourth version of the step function using the fftconvolve function for real Fast Fourier Transforms from the scipy.signal package. Because of rounding errors, you need to round the results returned by fftconvolve.



In [14]:

    
# Your solution here



In [ ]:

    
fig = plt.figure(figsize=(5, 5))
grid = init_grid.copy()
im = plt.imshow(grid, animated=True, interpolation='nearest', cmap='gray')
plt.close(fig)

anim = animation.FuncAnimation(fig, animate, 
                               frames=60, interval=50, blit=True)
HTML(anim.to_html5_video())

4. Modify 3a to model the Forest Fire model instead of the Game of Life, with f=0.001 and p=0.1.



In [18]:

    
# Your solution here



In [ ]:

    
cmap = colors.ListedColormap(['black', 'green', 'red'])

fig = plt.figure(figsize=(5, 5))
grid = init_grid.copy()
im = plt.imshow(grid, animated=True, interpolation='nearest', cmap=cmap, vmax=2)
plt.close(fig)

anim = animation.FuncAnimation(fig, animate, 
                               frames=100, interval=50, blit=True)
HTML(anim.to_html5_video())