September 4th, 2015
A common exercise in programming is making an instance of Conway's game of life. It's a great example of a simple set of rules that generate incredibly complex behavior. And if you're like me, you're sick of it. I mean, it's such an overused example and it got really old by now, which is a shame because the subject of cellular automata (of which the game of life is but a single example) is very fascinating. This inspired me to make a series of posts exploring the subject and it's applications. And I decided to do it in python, because python is so hot right now, you guys.
Simply put, a (1D) cellular automata is composed of an array of cells distributed in a circular pattern, where they can be found in one of two states: 0 or 1. In the beginning of each time step, every cell looks at its own state and the state of its nearest neighbors, and according to a deterministic rule decide whether to change its state or keep the current one.
Simple right? Taking by the name "cellular automata" the uninitiated reader might be mislead to think that the cells are the interesting part of this whole thing. Wrong! The rules that the cells use to update their state is where the money is at. Here is what a rule looks like
The input row represents all the possible neighborhoods of a cell. Depending on its own state and the states of its first neighbors to the left and right, the cell will turn to the state given in the output row. In cellular automata jargon, we say this is a rule with \(k=3\), because each cell look at its own state and the state of its direct neighbors to the left and to the right. If each cell looked instead to the two nearest neighbors in each direction this would be a \(k=5\) rule, rules with \(k=7\) look at its three nearest neighbors and so on. Each rule must take \(2^k\) inputs.
To deal with rules with more ease, we can treat each input as a 3 bit integer. Doing that and converting them to base 10, you should notice that they are ordered from 7 to 0. If, in this order, you take the all of the outputs as an 8 bit integer, the resulting number is the "name" of the rule. The rule above is therefore rule 126.
To get our cellular automata going, the first
thing we gotta do is to write a function that takes the rule's name and
generates a python list out of it. This should be quite simple to do with
some bit shifting and whatnot. I'm feeling lazy though, so I'll just use the
built-in bin function
print(bin(126))
#[stdout] 0b1111110
The problem is that it returns a string prefixed with "0b". No big deal, well strip the two first characters, transform each character into an integer, left-pad with zeros and reverse it (because the rules are organized in descending order, but lists are indexed in ascending order)
def make_rule(rule_name, k):
rule_name = int(rule_name)
rule_len = 2 ** k
rule = [*map(int, bin(rule_name)[2:])]
rule = [0] * (rule_len - len(rule)) + rule
rule = rule[::-1]
return rule
print(make_rule(126, 3))
#[stdout] [0, 1, 1, 1, 1, 1, 1, 0]
Next we need to make a function that take a list representing the state and
turn into an integer. It's just the reverse of what make_rule
does
def state_id(state):
# Thank god for the `base` argument in python's build-in `int` =)
return int(''.join(map(str, state)), base=2)
assert state_id([0,1,1]) == 3
assert state_id([1,0,1]) == 5
I know, this looks clunky, converting to and from strings just to convert an integer from one base to another, but it's a perfectly valid approach in my opinion, and not much slower than the other pure python alternatives (but more on that on a later post).
We'll also need a simple rolling window function (and make it a generator for extra python points)
from pprint import pprint as print
def rolling_window(arr, wsize):
arr = arr[-wsize//2 + 1:] + arr + arr[:wsize//2]
for i in range(len(arr) - wsize + 1):
yield arr[i:i + wsize]
print(list(rolling_window([1,2,3,4,5,6], 5)))
#[stdout] [[5, 6, 1, 2, 3], [6, 1, 2, 3, 4], [1, 2, 3, 4, 5], [2, 3, 4, 5, 6], [3, 4, 5, 6, 1], [4, 5, 6, 1, 2]]
Now applying a rule is just a matter of putting it all together
def apply_rule(cells, rule, k):
assert len(rule) == 2 ** k
return [rule[state_id(w)] for w in rolling_window(cells, k)]
assert apply_rule([0,0,1,1,0,0], make_rule(126, 3), 3) == [0, 1, 1, 1, 1, 0]
And it works! But not very interesting right now, right? Before trying to do something fancier, let's first put together all we did in a nice little class
import matplotlib.pyplot as plt
from random import randint
import numpy as np # for now we'll only use numpy to help with plotting
class CAutomaton:
def __init__(self, ncells, k, init='center'):
self.ncells = ncells
self.k = k
if init == 'center':
self.cells = [0] * self.ncells
self.cells[self.ncells // 2] = 1
elif init == 'random':
self.cells = [randint(0, 1) for _ in range(self.ncells)]
else:
self.cells = list(init)
assert len(self.cells) == self.ncells
@staticmethod
def make_rule(rule_id, k):
rule_len = 2 ** k
rule = [*map(int, bin(rule_id)[2:])]
rule = [0] * (rule_len - len(rule)) + rule
rule = np.array(rule)[::-1]
return rule
def state_id(self, state):
assert len(state) == self.k
return int(''.join(map(str, state)), base=2)
def apply_rule(self, rule):
assert len(rule) == 2 ** self.k
self.cells = [rule[self.state_id(w)]
for w in rolling_window(self.cells, self.k)]
return self
def iterate_and_plot(self, rule_id, niter=None, ax=None):
if niter is None:
niter = self.ncells // 2
rule = CAutomaton.make_rule(rule_id, self.k)
states = [self.cells] + [self.apply_rule(rule).cells for _ in range(niter)]
states = np.array(states)[::-1]
if ax is None:
fig, ax = plt.subplots(figsize=(10,10))
ax.pcolormesh(states, cmap="Greys")
ax.set_xlim(0, self.ncells)
ax.set_ylim(0, niter+1)
ax.set_title(r'$Rule\;%i$' % rule_id, fontsize=20)
ax.set_axis_off()
ax.set_aspect("equal")
return ax
I made a practical iterate_and_plot method to repeatedly apply
the rule to the cells and plot the history (time flowing from top to bottom).
I also added two possible initial conditions: 'center' (where
all cells start at state 0 except for the middle one) and
'random' (which should be pretty obvious).
Let's see what happens when we apply rule 126 a bunch of times with a simple initial condition
CAutomaton(ncells=520, k=3, init='center').iterate_and_plot(rule_id=126, niter=255) plt.show()
Whoa! Look at that nice Sierpinski triangle! If you're not sold on cellular automata by now, there's not much I can do for you, buddy. As you play with different rules you'll notice they show a variety of different behaviors
fig, ax = plt.subplots(ncols=2, nrows=2, figsize=(8,5)) CAutomaton(ncells=63, k=3).iterate_and_plot(rule_id=12, ax=ax[0,0]) CAutomaton(ncells=63, k=3).iterate_and_plot(rule_id=1, ax=ax[0,1]) CAutomaton(ncells=63, k=3).iterate_and_plot(rule_id=45, ax=ax[1,0]) CAutomaton(ncells=63, k=3).iterate_and_plot(rule_id=150, ax=ax[1,1]) fig.tight_layout() plt.show()
These patterns didn't passed unnoticed to Stephen Wolfram (of Mathematica© fame). He created an interesting taxonomy for cellular automata rules
fig, ax = plt.subplots(ncols=4, figsize=(10, 10)) CAutomaton(63, 3, 'random').iterate_and_plot(rule_id=250, niter=100, ax=ax[0]) CAutomaton(63, 3, 'random').iterate_and_plot(rule_id=1, niter=100, ax=ax[1]) CAutomaton(63, 3, 'random').iterate_and_plot(rule_id=30, niter=100, ax=ax[2]) CAutomaton(63, 3, 'random').iterate_and_plot(rule_id=126, niter=100, ax=ax[3]) fig.tight_layout() plt.show()
Class 1 rules normally converge to uniform states when subject to random initial conditions, like rule 250 above;
Class 2 rules never really converge, but get stuck in simple periodic patterns, like rule 1 above;
Class 3 rules are completely aperiodic and chaotic. Looking from a distance they usually look like white noise, like rule 30 above;
Class 4 rules are also aperiodic, but unlike class 3 rules they show some sort of organized structure even when looking at a distance, like rule 126. These are the most interesting rules, as they lie somewhere between order and chaos, which Wolfram referred as "interesting complexity".
You can take a look at all the 256 existing rules with \(k=3\) here, see if you can tell to which class each rule belongs.
I totally encourage you to play with this simple code. Try seeing what happens when you use \(k=5\) or \(k=7\), but keep in mind that the number of available rules grows with \(2^{2^k}\) (with \(k=9\) you have more possible rules than atoms in the universe :O), and most of them look like they belong to class 3. It's also very simple to generalize this code to cellular automata with more than two states, you just have to remember to codify the rules in base \(N_{states}\) instead of base 2. You might also wanna try extending it to more then one dimension. Conway's game of life is probably the most famous cellular automata, and it's a 2-dimensional one.
If you're interested to know more about cellular automata, Stanford Encyclopedia of Philosophy has an awesome article about it.
❦