Wolfram-style cellular automata (Python recipe) by Rick Muller

ActiveState Code (http://code.activestate.com/recipes/343386/)

Bitwise Versions. First, thanks for the recipe, Rick. While I have not read A New Kind of Science, I had the chance to see Wolfram lecture in support of the book and found myself amused and intrigued. I am eager to sit down and see what I can do with your code. That being said, here are bitwise arithmetic versions of the ca_data and pil_render functions.

<pre width="80"> def ca_data_bwa(height,width,dorandom,rulenum):

# Ensure that we are working with long integers, otherwise bitwise
# operations will not work properly for widths greater than 31 on most
# python implementations.
height = long(height)
width = long(width)
rulenum = long(rulenum)

# Create the first row, either randomly, or with a single 1
if dorandom:
    first_row = randint(0, (1 << width) - 1)
else:
    first_row = 1L << width/2L
results = [first_row]

rule = rulenum

for i in range(height-1L):
    data = results[-1]

    # Determine the new row based on the old data. We create appropriate
    # bitmasks to allow us to read/write to a particular bit and then use
    # bitwise operations to do the actual reading/writing.
    new = 0L
    for j in range(width):

        # Create a bitmask and then access the 'upper left' bit
        ul_bitmask = 1L << ((j - 1L) % width)
        ul_bit = (data & ul_bitmask) / ul_bitmask

        # Create a bitmask and then access the 'above' bit
        above_bitmask = 1L << j
        above_bit = (data & above_bitmask) / above_bitmask

        # Create a bitmask and then access the 'upper right' bit
        ur_bitmask = 1L << ((j + 1L) % width)
        ur_bit = (data & ur_bitmask) / ur_bitmask

        # Create a bitmask and then access the rule
        rule_bitmask = 1L << ((ul_bit << 2L) + (above_bit << 1L) + ur_bit)
        new_bit = (rule & rule_bitmask) / rule_bitmask

        # Create a bitmask and then adjust the value of new at j
        new_bitmask = new_bit << j
        new = new | new_bitmask

    results.append(new)
return results

def pil_render_bwa(data,height,width,fname="bs_bwa.png"): import Image, ImageDraw img = Image.new("RGB",(width,height),(255,255,255)) draw = ImageDraw.Draw(img)

for y in range(height):
    for x in range(width):
        if (data[y] & (1L << x)) / (1L << x): draw.point((x,y),(0,0,0))
img.save(fname,"PNG")
return

</pre>

(comment continued...)

(...continued from previous comment)

Unfortunately, bitwise ops are rather slow under Python, so the bitwise versions take roughly 6 times as long to run on my machine. So, no speed optimization for me, though your mileage may vary. I tried to optimize with various different approaches to bit shifting, but found that the << (left shift) operator was almost universally faster than pow and a number of other approaches.

On another dimension, these routines may provide an optimization. While I have not yet investigated the matter thoroughly, I suspect that the bitwise versions probably perform better than the originals with regard to memory consumption. While I am ignorant of the actual memory representation of Python long integers, it seems that, unless they are grossly inefficient, a list of arbitrarily long Python long integers must take up less space than an equivalent representation using a list of lists of 32-bit integers. Of course, the phrase "it seems reasonable that" often has very little to do with actual performance, so it would be nice if someone would check this out.

Regardless of their efficiency, these alternate versions should serve as basic examples of bit arithmetic in Python. With the exception of the randomization function (where using long integers opened up obvious optimizations too good to pass up) I have attempted to keep the flow of the functions as similar as possible. The comments are a bit sparse, so, should anyone be really puzzled, I can clarify further.

Enjoy,

LJ

profiled and speedy. I got fascinated by Wolfram's 1d automata experimentation as well. I went back to the math he presented in one of his papers, and came up with this implementation, which I profiled and tweaked until I got it as fast as I could make it without dropping into C:

import Image
import ImageOps
import random

def step(a, rule, k=2, r=1):
    nbrs = [a[c:] + a[:c] for c in range(-r, r+1, 1)]
    l = []
    for t in apply(zip, nbrs):
        result = 0
        for i in t:
            result = (result * k) + i
        l.append(result)
    return [((rule / (k ** v)) % k) for v in l]

def basicRun(rule, steps, stepper, seed=[1], k=2, r=1):
    seed = ([0] * steps) + seed + ([0] * steps)
    result = seed[:]
    for i in range(steps):
        seed = stepper(seed, rule, k=k, r=r)
        result += seed[:]
    return result, (len(seed), steps + 1)

def showResult(result, dims, k=2):
    i = Image.new("L", dims)
    i.putdata(result, (255 / (k - 1)))
    i = i.crop(i.getbbox())
    i = ImageOps.invert(i)
    i.load()
    i.show()

def runTest():
    lines = 400
    result, dims = basicRun(110, lines, step)
    showResult(result, dims)

if __name__ == "__main__":
    runTest()

step() is where all the real action occurs. That's where most of the tweaking happened. And while bitwise manipulation is slow, basic math seems to be fast.

Nice. Thanks, Dan. Go list comprehensions.

Er, Dave, sorry.

I came up with this shorter and faster implementation.

Time to compute Rule 110, 2500x2500:

1.9 s vs 8 s for David's solution (including image write)

0.9 s vs 6.8 s for David solution (without image write)

import png

def next_row(row, rule_bits):
    r = row[-1:] + row + row[:1]
    return [rule_bits[r[i]*4 + b*2 + r[i+2]] for (i, b) in enumerate(row)]

def save_png_data(rows, filename):
    ll = [[255 - b * 255 for b in row] for row in rows]
    png.from_array(ll, 'L').save(filename)

def run():

    rule = 110 # or 0b01101110

    (WIDTH, HEIGHT) = (2500, 2500)

    row = [0] * (WIDTH-1) + [1]

    s = bin(rule)[2:].zfill(8)
    rule_bits = [int(x) for x in s[::-1]]

    png_data = [row]

    for i in range(HEIGHT):
        row = next_row(row, rule_bits)
        png_data.append(row)

    save_png_data(png_data, "rule_%d_%d_%s.png" % (rule, WIDTH, HEIGHT))

if __name__ == "__main__":
    run()