A big "Game of Life"

Question

Our quest: Create a big simulation for Conway's Game of Life, and record the entire simulation history.

Current Approach: Cython is used for an iterate method. The history of life is placed in an HDF store, with help from pandas, and matplotlib is used to visually check results.

What I could use help with: While I'd love advice on anything, I'm especially interested in improving the pandas implementation in this code. Also, how could I make this file easier for others to read?

---

Explanation of definitions

• Z (ndarray, int32): 2D array of life (1 or 0).
• N (ndarray, int32): Temporary array used to tally the count of neighbors at each point on Z. N has the shape of Z minus one on each axis.
• Z_chunk (DataFrame, int32): History of the past few iterations of Z

iterate.pyx

#cython: wraparound=False, boundscheck=False, cdivision=True
#cython: profile=False, nonecheck=False, overflowcheck=False
#cython: cdivision_warnings=False, unraisable_tracebacks=False
import numpy as np
cimport numpy as np

cpdef iterate(Z, c):
    '''Element by elemenent iteration with optimized Cython.

    Args:
        Z (ndarray - int32) - Represents 2D space
        c (namedtuple) - Container for constants
    Returns:
        Z (ndarray - int32)
    '''

    N = np.zeros((c.rows-1, c.cols-1), dtype=np.int32)

    cdef int rows = c.rows
    cdef int cols = c.cols
    cdef int [:, :] N_ = N
    cdef int [:, :] Z_ = Z
    cdef int x, y

    with nogil:
        # Count neighbors
        for x in range(1, rows-1):
            for y in range(1, cols-1):
                N_[x, y] = (Z_[x-1, y-1] + Z_[x-1, y] + Z_[x-1, y+1] +
                            Z_[x,   y-1]              + Z_[x,   y+1] +
                            Z_[x+1, y-1] + Z_[x+1, y] + Z_[x+1, y+1])
        # Apply rules
        for x in range(1, rows-1):
            for y in range(1, cols-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 np.array(Z_)

Main file

import time
import os, sys
from matplotlib import pylab as plt
import numpy as np
import pandas as pd
from pandas import Series, DataFrame
from collections import namedtuple
import pyximport; pyximport.install()
from iterate import iterate

# Configure constants
#####################

# WARNING! Do not make these values large.
# The pandas part of this code is poorly written and does not scale well.
n_chunks = 3  
chunk_size = 10
n_iterations = n_chunks * chunk_size
rows = 128
cols = 128

lbls_row = ['row%05i'%(i) for i in range(rows)]
lbls_col = ['col%05i'%(i) for i in range(cols)]
lbls_iter_all = ['iter%05i'%(i) for i in range(n_iterations)]

# `c` is passed to all functions
Const = namedtuple('c', ['rows', 'cols', 'n_iterations', 'chunk_size', 'n_chunks',
                         'store_path', 'save_path', 'lbls_row', 'lbls_col',
                         'lbls_iter_all'])
c = Const(rows=rows, cols=cols, n_iterations=n_iterations, chunk_size=chunk_size,
          n_chunks=n_chunks, store_path='life_store.h5', save_path='life_results',
         lbls_row=lbls_row, lbls_col=lbls_col, lbls_iter_all=lbls_iter_all)

# Use `m` only for matplotlib parameters
dpi = 72.0
figsize = c.cols / float(dpi), c.rows / float(dpi)
Config_mpl = namedtuple('m', ['figsize', 'dpi'])
m = Config_mpl(figsize=figsize, dpi=dpi)


# Run the game
##############

start = time.time()

# Clear previous results if they exist
if os.path.exists(c.store_path):
   os.system('rm ' + c.store_path)
# Create new store
store = pd.HDFStore(c.store_path)

# Randomly place `1`s throughout the map
Z = np.random.randint(0, 2, (c.rows, c.cols)).astype(np.int32)
Z[0, :] = 0  # Clear boarders-- will act as a boundary
Z[-1, :] = 0
Z[:, 0] = 0
Z[:, -1] = 0

for i in range(c.n_chunks):
    print 'Chunk', i

    # Initialize Z_chunk
    lbls_iter = ['iter%05i'%(k) for k in range(c.chunk_size * i, 
                                               c.chunk_size * i + c.chunk_size)]
    columns = pd.MultiIndex.from_product([lbls_iter, c.lbls_col], names=['iter', 'col'])
    Z_chunk = DataFrame(index=c.lbls_row, columns=columns, dtype=np.int32)

    for j in range(c.chunk_size):
        Z_chunk.loc[:, lbls_iter[j]] = iterate(Z, c)

    store['chunk%03i'%i] = Z_chunk

print 'Life lasted: ', time.time() - start


# View results
##############

# Create a folder in the current directory
if not os.path.exists(c.save_path):
    os.makedirs(c.save_path)
else:  # Clear previous images if the directory exists
    cmd = 'rm ./' + c.save_path + '/*.png'
    os.system(cmd)

# Create figure
fig = plt.figure(figsize=m.figsize, dpi=m.dpi, facecolor="white")

ax = fig.add_axes([0.0, 0.0, 1.0, 1.0], frameon=False)
ax.get_xaxis().set_visible(False)
ax.get_yaxis().set_visible(False)

# Because of my sloppy HDF writing, we will need nested for loops here
for i in range(c.n_chunks):
    Z_chunk = store['chunk%03i'%i]
    lbls_iter = ['iter%05i'%(k) for k in range(c.chunk_size * i, 
                                               c.chunk_size * i + c.chunk_size)]
    for j, lbl in enumerate(lbls_iter):
        # Extract
        Z = Z_chunk.loc[:, lbl]
        # Show
        ax.imshow(Z, interpolation='nearest', cmap=plt.cm.gray_r)
        # Save
        fig.savefig('./{}/img_{}.png'.format(c.save_path, (1+i)*(j+1)))
        # Clear
        ax.cla()

# Close
plt.close(fig)
store.close()

References: Python Tutorial

---

Update!

As suggested by Curt, h5py is used instead of pandas to remove index support, and NumPy is used instead of Cython for clarity.

With the following values

n_chunks = 10
chunk_size = 100
n_iterations = n_chunks * chunk_size
rows = 1024
cols = 1024

Comparison:

• Cython - 31.8s
• scipy.signal.convolve - 229.5s
• NumPy vectorization - 49.1s

(Note: The NumPy vectorization method is from the NumPy Tutorial linked above in the reference. I don't understand the mathematics behind convolution well enough to explain why it's slower.)

The following code seems to work very well for big values. Even though the Cython version is slightly faster, both methods are plenty fast enough when compared to disk writing and image processing.

import time
import os, sys
from matplotlib import pylab as plt
import numpy as np
import h5py
from scipy.signal import convolve
from collections import namedtuple


def iterate(Z):
    # 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


# Configure constants
#####################

n_chunks = 10
chunk_size = 100
n_iterations = n_chunks * chunk_size
rows = 1024
cols = 1024

# `c` is passed to all functions
Const = namedtuple('c', ['rows', 'cols', 'n_iterations', 'chunk_size', 'n_chunks',
                         'h5_path', 'save_path'])
c = Const(rows=rows, cols=cols, n_iterations=n_iterations, chunk_size=chunk_size,
          n_chunks=n_chunks, h5_path='life.hdf5', save_path='life_results')

# Use `m` only for matplotlib parameters
dpi = 72.0
figsize = c.cols / float(dpi), c.rows / float(dpi)
Config_mpl = namedtuple('m', ['figsize', 'dpi'])
m = Config_mpl(figsize=figsize, dpi=dpi)


# Run the game
##############

start = time.time()

# Clear previous results if they exist
if os.path.exists(c.h5_path):
   os.system('rm ' + c.h5_path)
# Create new store
f = h5py.File(c.h5_path, 'w')

dset = f.create_dataset("Results", (c.rows, c.cols, c.n_iterations), dtype=np.int32,
                        chunks=(c.rows, c.cols, c.chunk_size))

# Randomly place `1`s throughout the map
Z = np.random.randint(0, 2, (c.rows, c.cols)).astype(np.int32)
Z[0, :] = 0  # Clear boarders-- will act as a boundary
Z[-1, :] = 0
Z[:, 0] = 0
Z[:, -1] = 0

for i in range(c.n_chunks):
    print 'Chunk', i

    # Initialize Z_chunk
    Z_chunk = np.zeros((c.rows, c.cols, c.chunk_size), dtype=np.int32)

    for j in range(c.chunk_size):
        Z_chunk[:, :, j] = iterate(Z)

    dset[:, :, i*c.chunk_size : i*c.chunk_size+chunk_size] = Z_chunk

print 'Time elapsed: ', time.time() - start


# View the results
##################

# Create a folder in the current directory
if not os.path.exists(c.h5_path):
    os.makedirs(c.h5_path)
else:  # Clear previous images if the directory exists
    cmd = 'rm ./' + c.h5_path + '/*.png'
    os.system(cmd)

# Create figure
fig = plt.figure(figsize=m.figsize, dpi=m.dpi, facecolor="white")

ax = fig.add_axes([0.0, 0.0, 1.0, 1.0], frameon=False)
ax.get_xaxis().set_visible(False)
ax.get_yaxis().set_visible(False)

# Iterate through `dset` to create images
for i in range(c.n_iterations):

    # Extract
    Z = dset[:, :, i]
    # Show
    ax.imshow(Z, cmap=plt.cm.gray_r)
    # Save
    fig.savefig('./{}/img_{}.png'.format(c.save_path, i))
    # Clear
    ax.cla()

# Close
plt.close(fig)

Example frame

Answer 1

DISCLAIMER: I don't know cython and have never used it, so if any of my advice doesn't apply because of cython limitations, feel free to disregard it.

1. Counting neighbors in a game board is very easy to do via convolution with the proper kernel. Below, I used SciPy's convolve function, not the much faster fftconvolve, because the latter only works on floats and so rounding back to integers / bools would be required. Nonetheless, I suspect for large game boards, the FFT method and rounding may still be faster. I doubt that Cython would speed up NumPy very much, but I would love to be proven wrong! (A Cython guide for NumPy users says Typical Python numerical programs would tend to gain very little as most time is spent in lower-level C that is used in a high-level fashion.) Anyway, to my eye at least this function is a bit easier to read and interpret.

---

import numpy as np
from scipy.signal import convolve

def iterate_game_of_life(board, neighbors_kernel=None):
    """Performs one iteration of Conway's game of life on a 2d numpy array of bools"""
    if neighbors_kernel is None:
        neighbors_kernel = np.ones(shape=(3, 3), dtype='int')
        neighbors_kernel[1, 1] = 0

    neighbors_count = convolve(board, neighbors_kernel, mode='same')
    has_three, has_two = neighbors_count == 3, neighbors_count == 2
    return np.logical_or(has_three, np.logical_and(board, has_two))

2. Why do you need pandas at all? Are you using it just to store HDF5 files? If so, I'd recommend looking at h5py instead. You may have to learn a bit more about how HDF5 file formats work, but I think it would be more efficient than Pandas. In particular, in an application like this where the game board is naturally a matrix, Pandas' DataFrame model that rows and columns need to have labels seems like annoying overhead instead of a feature. You could even store the iterations of the game as a single datacube instead of as separate HDF5 datasets.

3. Probably defining your own chunking will improve HDF5 performance at the largest sizes, but it depends on how you want to view the data. Do you want to view the entire board for a single iteration easily? Or do you want to view a small region of the board across multiple iterations easily? Unless you have particular views in mind, a good starting place would be to just use h5py's autochunking, i.e., don't define your own.

4. This is minor, but why are you using any interpolation at all in ax.imshow(Z, interpolation='nearest', cmap=plt.cm.gray_r)? I'd recommend ax.imshow(Z, interpolation='none', cmap=plt.cm.gray_r).