17 December 2006

From Bits to Cells: Simple Cellular Automata in Haskell, Part One

The little Haskell toy I wrote to simulate a dot-matrix printhead was useful in one sense: it helped me learn a little more Haskell! I couldn't help noticing that the results look like the transformations you get when applying simple cellular automata rules. So my next experiment is to implement those rules, and see what else I can learn from the results.

What follows is Literate Haskell, but there is no main program to run; it is for you to experiment with interactively, if you like.

Some Haskell code to implement simple cellular automata (CA) rules:
part one, round one.

The simplest type of one-dimensional CA behaves as follows:

- a cell may be on or off (one or zero).

- time ticks by in discrete intervals.

- at each tick, the state of each cell is updated using a function that
maps the current state of the cell to the left, the cell itself, and
the cell to the right to the new cell value.

Because this value-generating function depends on 3 bits of input, a
given CA rule actually consists of 2^3 = 8 possible mappings. The
mappings can be represented by the binary values 0b000 to 0b111
(decimal 0..7). By Wolfram's numbering convention for these simple
CA, we actually count down from 7 to 0. So, let's say we have a rule
that maps each possible set of left, center, and right values:

0b111, 0b110, 0b101, 0b100, 0b011, 0b010, 0b001, 0b000

to the new cell value:

0b0, 0b0, 0b0, 0b1, 0b1, 0b1, 0b1, 0b0

We can treat the resulting number, in this case 0b00011110, as the
rule number. This is rule 30 in Wolfram's schema; there are 2^8,
or 256, possible rules, numbered from 0 to 255.

First, let's make a function that, given a rule number and three
values (left, center, and right), returns the new cell state. We
turn the three values into an index to look up the appropriate bit
in the rule number.

>import Data.Bits

>genNextBit :: Int -> ( Bool, Bool, Bool ) -> Bool
>genNextBit rulenum ( left, center, right ) = rulenum `testBit` idx
>    where idx = ( if left then (4::Int) else (0::Int) ) .|.
>           ( if center then (2::Int) else (0::Int) ) .|.
>           ( if right then (1::Int) else (0::Int) )

Hmmm... it is lacking a certain elegance, and I don't quite like
the way the indentation rules work in this case, but let's test it
with a function that generates the 8 rule indices in the form of tuples:

>genRuleIndices :: [( Bool, Bool, Bool )]
>genRuleIndices = [ ( (val `testBit` 2), (val `testBit` 1),
>                   ( val `testBit` 0) ) | val <- [(7::Int),
>                   (6::Int)..(0::Int)] ]

Now if we write:


we get:


and if we write:

map (genNextBit 30) genRuleIndices

this expression curries us a function (mmm... curry!) which takes
a starting state and applies rule 30, then feeds it the possible
starting states. The result is:


That looks like the bit pattern for rule 30.

Just for fun, let's confirm by making a function that will
translate a list of output bit values back into a rule number.
The signature should look like this:

sumBitVals :: [Bool] -> Int

And we want the list

>test_bit_vals = [False,False,False,True,True,True,True,False]

to map back to 30. Take a shot at it yourself; you might find
that the result is instructional I'll wait.

(Musical interlude).

Did you try it? Let's look at some possible implementation strategies.
We could make a list of the powers of two and then do some list
manipulation to get the powers of two that correspond to our one-bits
summed up:

>powers_of_two_8_bits = reverse ( take 8 (iterate (2 *) 1) )


(An aside: since this is generated by a function that takes no
parameters at runtime, it would seem like a sufficiently smart
compiler could generate this list and stick it in the resulting
program so that no run-time calculation is required to generate
it. Does GHC do this? I don't know.)

Anyway, our implementation. We can tuple up our powers of two with
our bit values:

>tups = zip powers_of_two_8_bits test_bit_vals


Then we can turn this back into a list of only the powers of two that are
"on," and sum the results:

sum (map (\ tup -> if snd tup then fst tup else 0) tups )


It seems like this should be simpler still. I looked in the standard
library for a function to filter one list using a list of boolean
values as a comb. Let's write one:

>combListWithList :: [a] -> [Bool] -> [a]
>combListWithList [] [] = []
>combListWithList ls comb =
>   if (head comb) then (head ls) : combListWithList (tail ls) (tail comb)
>       else combListWithList (tail ls) (tail comb)

combListWithList powers_of_two_8_bits test_bit_vals


That seems good, although it doesn't express the idea that the lists
have to be the same length. It still seems amazing to me that I can
reel off functions like that in Haskell and have them work right the
first time! Now we can produce our final function:

>sumBitVals :: [Bool] -> Int
>sumBitVals ls = sum $ combListWithList powers_of_two_8_bits ls

sumBitVals test_bit_vals


There is probably an elementary way to implement the above function
using zipWith instead of my combList, but I'll leave that to you;
leave a comment if you figure it out!

As another aside, here is another function I came up with to generate
this bit vals; it doesn't rely on a specific length.

>foldBitVals :: [Bool] -> Int
>foldBitVals ls = snd (foldr
>   (\ flag tup -> if flag then ((fst tup) * 2, (snd tup) + (fst tup))
>                           else ((fst tup) * 2, snd tup) )
>   (1, 0) ls )

foldBitVals test_bit_vals


Take a moment to understand the above function; it is a little strange.
I perform a fold using a tuple. The first element is the power of two,
so it keeps track of our place in the list, and the second is an
accumulator as we add up our values. This is a trick for stashing
a state that is more complex than a single value into a fold.

Anyway, that was a long diversion, so let's end for now. Next time we'll look at ways to represent our CA universe as it evolves by applying our rules. As always, your comments are welcome!


Neil Mitchell said...

> An aside: since this is generated by a function that takes no...
> Does GHC do this? I don't know.)

I strongly suspect it doesn't. If you did, then you might (at runtime) have to do an unbounded amount of computation. In general its not worth doing - either the computation is small (low runtime cost anyway), or its big (high compile time cost, possibly large amount of result)

Another point, use the power of the pattern match!

combListWithList (lh:lt) (ch:ct) =
if ch then lh : combListWithList lt ct else combListWithList lt ct

(its also faster)

>foldBitVals :: [Bool] -> Int
>foldBitVals ls = snd (foldr
> (\ flag (t1,t2) -> if flag then (t1 * 2, t2 + t1)
> else (t1 * 2, t2) )
> (1, 0) ls )

Which, of course, you can rewrite as something like:

sum [(a,b) <- zip ls bit_values | a]

(where bit_values = [1,2,4,...] - similar to as defined previously)

The Alternate Moebyus said...


I guess the neatest way to do the base-conversion is something on the lines of (coded without ghc, so it might have a bit of handwaving):

b2num :: [Bool] -> Int
b2num = foldl (\ x y -> 2 * x + fromEnum y) 0

I use something similar for the Project Euler problems (incidentally, you might want to check them out - solving those problems is a nice way to learn a new language)

Unknown said...

"It still seems amazing to me that I can
reel off functions like that in Haskell and have them work right the
first time!"

I know! *Squee*

"There is probably an elementary way to implement the above function using zipWith instead of my combList, but I'll leave that to you; leave a comment if you figure it out!"

Indeed, that's the approach I took. Since zipWith is just zip then map, you practically wrote it yourself! zip is even equivalent to zipWith (,) (where (,) is the pair constructor, read as 'pair'). If you're just going to map the pairs together, you needn't construct them in the first place. Here's what I came up with:

> sumBitVals :: [Bool] -> Int
> sumBitVals bits = sum (zipWith intify (reverse bits) powersOfTwo)
> where powersOfTwo = map (^2) [0..]
> intify bit value = if bit then value else 0

I like the way you used iterate (*2) to generate the powersOfTwo. I'd have never thought of that. Then again, I'm not used to such a bizarre concept as "iterating" in Haskell. :) Recursion For The Win!

Paul R. Potts said...

Hi Alternate Moebyus,

That b2num function is crazy... I have to step through the expansion steps to see how it works!

On [], foldl doesn't ever call our lambda, so we just get the initial value.

On [False], fromEnum yields 0, and our initial x is zero, so we get 2 * 0 + 0, or 0.

On [True] fromEnum yields 1, so we get 2 * 0 + 1, or 1.

On [False, True] we are folding left, so we get 0 from the first pass and then 2 * 0 + 1 is generated.

On [True, False] we are folding left so we get 1 for the first pass, then do 2 * 1 + 0, for 2.

On [True, True] we get 1 from the first pass, then 2 * 1 + 1, for 3.

At first glance this function looks to me like it should not work for all lists... but it does seem to work!

The Alternate Moebyus said...

While I'd rather not derive it again :) it's supposed to be the optimized version of evaluating a polynomial (a finite power series, I guess) at 2.

a0 + 2(a1 + 2(a2 + ..)