# Performance problems with 1D cellular automata experiment in Python with Numpy

I'm an experienced programmer but relatively new to Python (a month or so) so it's quite likely that I've made some gauche errors. I very much appreciate your time spent taking a look at my code.

My intention is to do some experimentation with 1D cellular automata. For now, just totalistic ones (rules only based upon sum of neighbor cells of given symmetric distance). I started off with code from the book Think Complexity which is available online and have adapted it slightly. I find the performance to be quite dismal. I am using Numpy arrays, and I think I may be able to improve performance markedly if I can figure out how to use stride_tricks but thus far I've been unsuccessful at that. For reasons I do not understand, it takes a couple seconds to run even 100 generations on a field 200 cells wide on an i7 processor with 16GB of RAM.

Even more distressing is that after running for awhile this code eats astonishing amounts of RAM! I am using matplotlib to generate images but am not terribly concerned about its speed, though it is slow. My experiments will be numeric when I can run them and few images will need to be generated. The time of a few seconds I mentioned earlier is solely the duration of running the run() method and does not include any of the rendering. I will be sharing the simple loop used for rendering simple images of each ruleset though because I imagine I am doing something in there that prevents the Python garbage collector from working properly, leading to the 10GB+ memory usage after running for awhile.

The cellular automata class:

class CellularAutomata1DTotalistic(object):

def __init__(self, ruleset, width, radius=1, isRing=False):
self.cells = zeros((1, width + 1), dtype=int8)
self.width = width
self.rule = ruleset
self.isRing = isRing
self.table = self.buildTable(self.rule)
self.next = 0

def run(self, steps):
for i in xrange(steps):
self.step()

def step(self):
i = self.next
self.next += 1

if i == len(self.cells):
self.doubleCells()

self.cells[i, j] = self.table[self.cells[i - 1, j - self.radius:j + self.radius + 1].sum()]

def doubleCells(self):
newcells = zeros(self.cells.shape, dtype=int8)
self.cells = vstack((self.cells, newcells))

def reset(self):
self.cells = zeros((1, self.width + 1), dtype=int8)
self.next = 0

def startSingle(self):
self.reset()
self.cells[0, (self.width - 1) / 2] = 1
self.next = 1

def startWith(self, val):
self.reset()
# TODO: Center the provided pattern
self.cells[0] = val[:]
self.next = 1

def randomize(self, p):
for i, x in enumerate(random.random(width - 2)):
self.cells[0, i + 1] = int(x < p)
self.next = 1

def buildTable(self, rule):
table = {}
#bound = 2 ** (self.radius + 1)
bound = (self.radius * 2) + 2
for i, bit in enumerate(binary(rule, bound)):
table[bound - 1 - i] = bit
return table

def getLatest(self):
return self.cells[self.next - 1, :]

def getEntropy(self):
entropy = []
for i in xrange(1, self.next - 1):
p = float(self.cells[i].sum()) / float(self.width)
entropy.append(-(p * log(p, 2)))
return entropy

def getLiveCounts(self):
liveCounts = []
for i in xrange(0, self.next - 1):
liveCounts.append(self.cells[i].sum())
return liveCounts

def get_array(self, start=0, end=None):
"""Gets a slice of columns from the CA, with slice indices
(start, end).  Avoid copying if possible.
"""
if start == 0 and end == None:
return self.cells[0:self.next, :]
else:
return self.cells[0:self.next, start:end]

def binary(n, digits):
"""Returns a tuple of (digits) integers representing the
integer (n) in binary.  For example, binary(3,3) returns (0, 1, 1)"""
t = []
for i in range(digits):
n, r = divmod(n, 2)
t.append(r)

return tuple(reversed(t))


The simple program I am using to generate images of all rules run for 100 steps on a single black cell start state which has strictly increasing RAM usage:

import CellularAutomata1DTotalistic as CA1DT
import CellularAutomata1DRenderer as CA1DR

r = CA1DR.PyplotRenderer()

ca.startSingle()
ca.run(100)
r.draw(ca)
r.save('t[' + str(rule) + '][' + str(radius) + '].png')

for rule in xrange(2 ** ((radius * 2) + 2)):


Thanks very much for any tips! I should mention that I was worried the resizing of the numpy array by doubling its size when it needed more space might have been what was slowing it down so I created an alternative version which pre-created the cells array of the necessary size in the run function and I saw no difference whatsoever in the performance (using IPythons timeit functionality on runs of 1000 steps on fields 2000 cells wide).

-

This is the slow bit of your program:

for j in xrange(self.radius, self.width - self.radius):
self.cells[i, j] = self.table[self.cells[i - 1, j - self.radius:j + self.radius + 1].sum()]


The key to making faster numpy code is to vectorize. You want to get rid of explicit loops and instead use operations that contain implicit loops:

First thing, we can break this loop into two loops, one will do the summing and the other will do the table lookup:

    sums = zeros(self.width + 1, uint8)
sums[j] = self.cells[i - 1, j - self.radius:j + self.radius + 1].sum()

self.cells[i, j] = self.table[sums[j]]


The trick to vectorizing the first for loop is to realize that you can produce it by taking the cumulative sum and then subtracting. The sum of elements 10 - 16 is the sum of the elements up to 16 minus the sum of elements up to 10. We can ask numpy to calculate the sums once and then subtract the elements 2*self.radius elements apart like this:

    sums = zeros(self.width + 1, uint8)
cumsum = self.cells[i - 1].cumsum()


However, we aren't using the beginning or the end of this sums array. So let shift the indexes and get rid of that.

    cumsum = self.cells[i - 1].cumsum()

self.cells[i, j] = self.table[sums[j - self.radius]]


Furthermore, that second loop can be rewritten to start from zero by shifting the indexes

    for j in xrange(self.width - 2*self.radius):
self.cells[i, j + self.radius] = self.table[sums[j]]


Next, self.table is a dictionary. But it'll be easier to vectorize if we make it a numpy array.

    bound = (self.radius * 2) + 2
table = zeros(bound, uint8)
for i, bit in enumerate(binary(rule, bound)):
table[bound - 1 - i] = bit
return table


Now the loop looking up the table can be rewritten as:

    self.cells[i, self.radius:self.width - self.radius] = self.table[sums]


This gives me:

    cumsum = self.cells[i - 1].cumsum()


General code review:

class CellularAutomata1DTotalistic(object):

def __init__(self, ruleset, width, radius=1, isRing=False):
self.cells = zeros((1, width + 1), dtype=int8)


Why are you doing width + 1? That makes the code more complicated because you've got this extra cell.

        self.width = width
self.rule = ruleset
self.isRing = isRing


Python standard is to use lowercase_with_underscores for local variables, parameter, and attributes

        self.table = self.buildTable(self.rule)
self.next = 0

def run(self, steps):
for i in xrange(steps):
self.step()

def step(self):
i = self.next


I recommend against using i here as its not immeadiately obvious what you are doing with it.

        self.next += 1


I recommend having self.counter = itertools.counter() in __init__. Then self.counter.next() will return the next number.

        if i == len(self.cells):
self.doubleCells()

self.cells[i, j] = self.table[self.cells[i - 1, j - self.radius:j + self.radius + 1].sum()]

def doubleCells(self):


Python standard convention is lowercase_with_underscores for method names

        newcells = zeros(self.cells.shape, dtype=int8)
self.cells = vstack((self.cells, newcells))


It might actually make more sense for self.cells to be a python list. Then you just append onto it. You may find it more useful to have it an array, but its something to consider

    def reset(self):
self.cells = zeros((1, self.width + 1), dtype=int8)
self.next = 0


Consider whether you really want to have reset instead of just making a new object

    def startSingle(self):
self.reset()
self.cells[0, (self.width - 1) / 2] = 1
self.next = 1

def startWith(self, val):
self.reset()
# TODO: Center the provided pattern
self.cells[0] = val[:]
self.next = 1

def randomize(self, p):
for i, x in enumerate(random.random(width - 2)):
self.cells[0, i + 1] = int(x < p)
self.next = 1


Use numpy's random functions. self.cells[0, 1:-1] = numpy.random.random(width - 2) < p. Actually self is undefined here suggesting you haven't done any testing on this function.

    def buildTable(self, rule):
table = {}
#bound = 2 ** (self.radius + 1)


Don't keep dead code, just delete it

        bound = (self.radius * 2) + 2
for i, bit in enumerate(binary(rule, bound)):
table[bound - 1 - i] = bit
return table


Lookup numpy.unpackbits. It'll extract the bits from the number for you.

    def getLatest(self):
return self.cells[self.next - 1, :]

def getEntropy(self):
entropy = []
for i in xrange(1, self.next - 1):
p = float(self.cells[i].sum()) / float(self.width)
entropy.append(-(p * log(p, 2)))
return entropy


This could be readily vectorized, but its probably not performance critical.

    def getLiveCounts(self):
liveCounts = []
for i in xrange(0, self.next - 1):
liveCounts.append(self.cells[i].sum())
return liveCounts

def get_array(self, start=0, end=None):
"""Gets a slice of columns from the CA, with slice indices
(start, end).  Avoid copying if possible.
"""
if start == 0 and end == None:
return self.cells[0:self.next, :]
else:
return self.cells[0:self.next, start:end]


What if start is supplied but end is not?

def binary(n, digits):
"""Returns a tuple of (digits) integers representing the
integer (n) in binary.  For example, binary(3,3) returns (0, 1, 1)"""
t = []
for i in range(digits):
n, r = divmod(n, 2)
t.append(r)

return tuple(reversed(t))


This isn't a good place for a tuple, as it is a list of digits, not a collection of hetrogenous items.

-
First off, thanks very much, your feedback is golden! Some of the things you questioned, such as various bits dealing with the unused cells on each end of the array were there to make it easier to adapt this code to allow the cells array to wraparound, that comes once I've got this down. You were right about randomize, I have never used it, should just remove it. You mentioned using a Python list rather than numpy array. I was using a numpy array because I thought it was the way to go for performance with this kind of thing, am I incorrect in that assumption? Again, thanks very much! –  otakucode Feb 28 '13 at 0:36
@otakucode, numpy arrays are slower than python lists if used the same way. numpy arrays are faster only if you can use vector operations. If you are explicitly looping over the array you aren't gaining any performance. (Memory consumption will be down, but speed will not improve) –  Winston Ewert Feb 28 '13 at 0:53
@otakucode, regarding the wraparound, I suspect that unused cells aren't a good way of doing that. But I'll leave that up to you. –  Winston Ewert Feb 28 '13 at 0:58
I implemented the changes you suggested and it significantly improved the speed of my code. However, it still very rapidly consumes memory when run in the loop I posted. I'm very puzzled by this as I don't see how I could be leaving dangling references. Do you have any tips about that? –  otakucode Feb 28 '13 at 2:15
@otakucode, when I run your code (with my modifications) I'm not seeing the memory consumption issue. Perhaps you can supply a runnable example of your code in its current state? –  Winston Ewert Feb 28 '13 at 2:20