As luck would have it, my day job will require me to teach a mini-class on writing computer simulations for scientific purposes next week, and one of the examples I thought I’d include turned out to be just the thing I needed to provide an illustration for the review of John Allen Paulos’ Irreligion (2007) I’ve got in the works. Since the details are a little beyond the scope of that review, I might as well make another post out of them.

First, the punchline:

This is a fifty-frame animation of Conway’s Game of Life which begins with a 32×32 grid of random values and successively advances by applying the transition rules of that famous automaton. The purpose of the exercise is twofold: first, to demonstrate some basic programming techniques and show how much you can do in Python with half an hour of free time; and second, to see how persistent and cyclic features arise from random configurations.

I chose to establish periodic boundary conditions, so that patterns leaving one side wrap around to the other (readers of my generation will recall the game Asteroids). This means we’re playing Life on a torus — is there anything donuts can’t do? As before, we’ll be using the Python language with a few extra modules for stuff like visual output.

Now, the code!

#!/usr/bin/python
# conway.py
# Blake Stacey
# (bstacey at alum.removethis.mit.andthis.edu)
import math, scipy, random, pylab

def iTest(x, p):
	if(x <= p):
		return 0
	else:
		return 1

fP = 0.5            # probability of dead square
iGridEdge = 32      # size of grid
iTrials = 50        # number of times to loop

random.seed()
	
imThis = scipy.array(scipy.rand(iGridEdge, iGridEdge))
imThat = scipy.array(scipy.zeros((iGridEdge, iGridEdge)))

for i in range(iGridEdge):
	for j in range(iGridEdge):
		imThis[i][j] = iTest(imThis[i][j], fP)
		
pylab.ion()
pylab.hold(False)

for t in range(iTrials):
	for i in range(iGridEdge):
		for j in range(iGridEdge):
			# bounds checking
			iPlus = (i + 1) % iGridEdge
			iMinus = (i - 1) % iGridEdge
			jPlus = (j + 1) % iGridEdge
			jMinus = (j - 1) % iGridEdge
			
			# count number of neighbors
			iNext = (imThis[i][jPlus]
					 + imThis[iPlus][j]
					 + imThis[iMinus][j]
					 + imThis[i][jMinus]
					 + imThis[iPlus][jPlus]
					 + imThis[iMinus][jMinus]
					 + imThis[iPlus][jMinus]
					 + imThis[iMinus][jPlus])
			iCurrent = imThis[i][j]
			
			# update rule
			if(iCurrent == 0):
				if(iNext == 3):
					imThat[i][j] = 1
				else:
					imThat[i][j] = 0
			else:
				if(iNext <= 1):
					imThat[i][j] = 0
				elif(iNext > 3):
					imThat[i][j] = 0
				elif((iNext == 2) or (iNext == 3)):
					imThat[i][j] = 1					
	# end loop over elements

	# swap buffers
	imOther = imThis
	imThis = imThat
	imThat = imOther
	
	# print summary statistic and display image
	print sum(sum(imThis)) / float((iGridEdge * iGridEdge))
	pylab.imshow(imThis, cmap = pylab.cm.prism)
	pylab.draw()

By sticking a pylab.savefig command after the pylab.imshow call, the images can be written to files:

	pylab.imshow(imThis, cmap = pylab.cm.prism)
	pylab.savefig("frame_" + str(t) + ".png")
	pylab.draw()

There’s an awful lot of whitespace in the resulting images, and the quickest hack for, well, hacking it off is to use ImageMagick:

 convert -crop 682x682+272+83 frame_0.png trimmed_0.png

I actually did all fifty frames with a bash loop:

for i in `seq 0 49`; do convert -crop 682x682+272+83
  frame_$i.png $i.png; done

Noting that the images had more resolution than was really necessary, I sized them down:

for i in `seq 0 49`; do convert -resize 100x100 $i.png
  resize_$i.png; done

Then I made the frames into an animated GIF, realized that they were out of order because the single-digit names were not prefixed with 0, fixed that problem, and made a new animation:

convert -delay 20 -loop 0 resize*.png animation.gif

The result was the image seen above. Here, I used a random initial configuration, because I had allotted myself half an hour to do the whole thing and I didn’t feel like specifying anything more detailed. However, you can put in a feature like a glider with code like this:

imThis = scipy.array(scipy.zeros((iGridEdge, iGridEdge)))
imThat = scipy.array(scipy.zeros((iGridEdge, iGridEdge)))
imThis[5][5] = 1
imThis[5][6] = 1
imThis[5][7] = 1
imThis[6][5] = 1
imThis[7][6] = 1

Lately, I’ve found myself growing tired of the usual creationist-bashing and snarky atheology which is the easiest type of science blogging to produce. Consequently, I plan to transition Sunclipse over to more math and physics posts, more “breakfast experiments” with programming and data mining, and more happy ramblings in the “bibliophilia” category.