Tuesday, December 19, 2006

Evaluating cellular automata is comonadic

Paul Potts's post inspired me to say something about cellular automata too.

So here's the deal: whenever you see large datastructures pieced together from lots of small but similar computations there's a good chance that we're dealing with a comonad. In cellular automata we compute the value of each cell in the next generation by performing a local computation based on the neighbourhood of that cell. So cellular automata look like they might form a good candidate for comonadic evaluation.

I want to work on 'universes' that extend to infinity in both directions. And I want this universe to be constructed lazily on demand. One way of doing that is to represent a 'universe' as a centre point, a list of all elements to the left of that centre point and a list of all elements to the right. Here's a suitable
type:


> data U x = U [x] x [x]


For example U [-1,-2..] 0 [1,2..] can be thought of as representing all of the integers in sequence.

But this actually contains slightly more information than a list that extends to infinity both ways. The centre point forms a kind of focus of attention. We could shift that focus of attention left or right. For example consider

U [-2,-3..] (-1) [0,1..]

This represents the same sequence of integers but the focus has been shifted left. So think of the type U x as being a doubly infinite sequence with a cursor. (In fact, this makes it a kind of zipper.)

We can formalise the notion of shifting left and right as follows:


> right (U a b (c:cs)) = U (b:a) c cs
> left (U (a:as) b c) = U as a (b:c)


An object of type U is semantically like a C pointer into a const block of memory. You can increment it, decrement it and dereference it using the function I'll call coreturn below.

As U is a kind of list structure, it needs a map. In fact, we can define fmap for it:


> instance Functor U where
> fmap f (U a b c) = U (map f a) (f b) (map f c)


Now the fun starts. First I'll bemoan the fact that Comonads aren't in the standard Haskell libraries (at least I don't think they are). So I have to define them myself:


> class Functor w => Comonad w where
> (=>>) :: w a -> (w a -> b) -> w b
> coreturn :: w a -> a
> cojoin :: w a -> w (w a)
> x =>> f = fmap f (cojoin x)


cojoin is the dual to the usual join function. I've chosen to do things the category theoretical way and define =>> in terms of cojoin.

And here's why U forms a Comonad:


> instance Comonad U where
> cojoin a = U (tail $ iterate left a) a (tail $ iterate right a)
> coreturn (U _ b _) = b


Look closely at cojoin. It turns a into a 'universe' of 'universes' where each element is a copy of a shifted left or right a number of times. This is where all the work is happening. The reason we want to do this is as follows: we want to write rules that work on the local neighbourhoods of our universe. We can think of a universe with the cursor pointing at a particular element as being an element with a neighbourhood on each side. For example, we can write a cellular automaton rule like this:


> rule (U (a:_) b (c:_)) = not (a && b && not c || (a==b))


In order to apply this everywhere in the universe we need to apply the rule to each possible shift of the universe. And that's what cojoin does, it constructs a universe of all possible shifts of a. Compare with what I said here. So believe it or not, we've already written the code to evaluate cellular automata. u =>> rule applies the rule to u. The rest is just boring IO:


> shift i u = (iterate (if i<0 then left else right) u) !! abs i
>
> toList i j u = take (j-i) $ half $ shift i u where
> half (U _ b c) = [b] ++ c
>
> test = let u = U (repeat False) True (repeat False)
> in putStr $
> unlines $
> take 20 $
> map (map (\x -> if x then '#' else ' ') . toList (-20) 20) $
> iterate (=>> rule) u


Lazy infinite structures, comonads, zippers. I think I'm just beginning to get the hang of this functional programming lark! Over Xmas I might feel ready to try writing a piece of code longer than a dozen or so lines.

Anyway, I must end with a credit. I probably wouldn't have come up with this if I hadn't read this paper by Uustalu and Vene.

16 comments:

Saizan said...

I like very much this lazy solution, but how would you generalize it to spaces with more dimensions?
(Like in Life where you have 8 other cells in your neighborhood.)
I can imagine increasing the number of the lists, but not very well how to travel this new universe, especially how to reach the cells in the "oblique" directions

sigfpe said...

I think there may be an easy but inefficient solution to the problem of working in 2D using the type U (U Bool). I'd imagine 'up' and 'down' to be defined like 'fmap left' and 'fmap right' or something like that. I think you can then construct a new 'cojoin' with a little bit of work. But this might be a bit inefficient. But maybe not - I'd have to think about it harder but all those 'fmap left's might still only do O(1) work per cell in total. When I next have a moment...

sigfpe said...

Here are some more details:

> data U2 x = U2 (U (U x)) deriving Show

> instance Functor U2 where

> fmap f (U2 u) = U2 $ fmap (fmap f) u

> instance Comonad U2 where

> coreturn (U2 u) = coreturn (coreturn u)

> cojoin (U2 u) = fmap U2 $ U2 $ roll $ roll u where

> iterate1 f = tail . iterate f

> roll a = U (iterate1 (fmap left) a) a (iterate1 (fmap right) a)

(Sorry about lack of indentation.) I have the game of life running fine.

The strategy for U2 is much the same as that for U. cojoin makes a 2D grid of 2D grids where each inner grid is shifted an amount corresponding to its location within the outer grid. A similar strategy can be used to implement many other comonads. In fact, that was the ulterior motive in writing this CA code: getting intuition about what cojoin typically looks like.

Saizan said...

well I'm quite lost on the intuition of your roll, and why you use fmap U2 in cojoin.
However, a rule would be written like this?
rule (U2 (U ((U (a:_) b (c:_)):_) (U (d:_) x (f:_)) ((U (g:_) h (i:_)):_))) = ...
where the letters have this alignment?
a b c
d x f
g h i

sigfpe said...

Here is the exact rule I used in my code:

rule' (U2 (U
(U (u0:_) u1 (u2:_):_)
(U (u3:_) u4 (u5:_))
(U (u6:_) u7 (u8:_):_))) =
let n = length $ filter id [u0,u1,u2,u3,u5,u6,u7,u8] in
u4 && (n==2 || n==3) || (not u4) && n==3

Your code is amazingly close!

Maybe I'll tidy up my code and fill it with useful comments so I can post the entire thing here.

U2 converts a 1D universe of 1D universes into a 2D universe. 'fmap U2' converts...take a deep breath...a 2D universe of 1D universes of 1D universes into a 2D universe of 2D universes. 'fmap U2 . fmap' converts a 1D universe of 1D universes of 1D universes of 1D universes into a 2D universe of 2D universes. Conceptually, a 2D universe is the same thing as a 1D universe of 1D universes. The guts of the code is written in terms of 1D universes, and 'fmap U2 . U2' is just a bit of fluff at the end to convert everything back to 2D universes so cojoin has the right type.

The 'roll' is just taking advantage of the thing I said in an earlier comment. If you think of the inner U's as columns and the outer one as rows, then fmap left and right shift up and down. The slightly weird thing about 'roll' is that exactly the same piece of code does rows and columns whereas you might expect one function with fmaps for columns and one function without for rows. You can think of fmap as 'ducking down' a layer of U-flavoured onion skin. The two rolls do different things because each time it is used there are different layers of onion skin in place. The best way to make sense of it is to do what I did - make U and U2 instances of Show that can print stuff out with a pretty layout and then go through cojoin one step at a time.

kowey said...

When I first saw this, I was frightened off by the seeming abundance of scary new concepts: cellular automata, comonads, functors (for no good reason I still get nervous when I see that word), universes. Then I read it a second time, read it instead of skimming it, and realised that oh... there really is nothing to it.

The comonad stuff will take a while to sink in, but reading your post, I get the impression that they are not fundamentally deeper than monads. No cobrain is required to understand them, after all. Note: fellow weaklings might do well to have Cale's Monads as containers open in a separate window or browser tab.

Thanks for this post. For starters, the literate Haskell helps, pasting and running is useful. More importantly, the post had this wonderful effect of mutual reinforcement, in which several things one does not know well get tied together, and as a result all become clearer. Not only do I now have a slight taste of comonad, I also have now more than heard about cellular automata, and my tenuous grasp of zippers is now strenghtened. Makes me wish zippers were first presented in list terms from the very beginning.

What might have been helpful is some means for the easily frightened to realise that the post contains easy stuff, but short of inventing a Haskell coloured belts system, [this post is rated yellow belt; if you understand monads in terms of >>= and return, etc], I don't know what such a mechanism would consist of. Another potentially useful aid might be to put the example automaton rule into words, or even pictures. Incidentally, this particular example rule does not seem to make use of c. Is that right?

sigfpe said...

I think you're right, my rule doesn't seem to use c. I just made up random rules until I liked the result.

there really is nothing to it
Yes! There are some incredibly difficult looking papers out there but when you actually get down to implementing examples they can suddenly seem close to trivial. Check out the comonad I define here for an example that's reminiscent of the kinds of adding-up-lists-of-prices examples given in books like "Learn C in 21 days". And yet the paper I drew on is scary as hell!

Anonymous said...

You should try implementing the "Cellular Potts" model developed by Glazer & Granier. It's more physically meaningful than CA.

kathy green said...

Whoa, and I thought this site was about the iPhone)!

geophf said...

Are all Zippers that are instances of Functor Comonads? Conversely, and simply, are all Comonads Zippers that are Functors? I've done a bit of reading on (co)monads, and the above is the current straw I'm grasping at for understanding comonads better. Am I missing the boat entirely?

As for the cellular automata implementation, Mathworld, etc, has a list of standardized rules. Rule 30 is demonstrated graphically as the banner of my company's webpage: http://www.logicaltypes.com/.

Here is a way to pattern-match to build these standardized rules:

> rule30 (U (a:_) b (c:_)) = bitfrom (rule30' (bitof a b c))
> where rule30' 111 = 0
> rule30' 110 = 0
> rule30' 101 = 0
> rule30' 100 = 1
> rule30' 011 = 1
> rule30' 010 = 1
> rule30' 001 = 1
> rule30' 000 = 0
> bitof x y z = 100 * t x + 10 * t y + t z
> t = fromEnum
> bitfrom = toEnum

(Ah, well, the indentation is not preserved, but you get the idea.)

Since all the standardized rules follow the above pattern (where the rule number is the bit-pattern result from the three-bit-pattern input), it is trivial in Haskell to write a generalized

rule :: Enum a => Int -> U a -> a

for all the rules that follow that pattern.

scc said...

Pondering rewriting this but representing the 'universe' as a function from an integer leads one to realise that we get a comonad from functions from any monoid as follows:

coreturn a = a mempty
cojoin a = curry (a . uncurry mappend)
fmap f = (f.)

Which allows us to handle automata on 1d, square grids, hex grids, and so on. Albeit perhaps inefficiently.

I haven't tried this so it's probably embarrassingly wrong.

Alex S. said...

So I've gone forward and implemented the game of life, without any fancy printing, but still, it seems to work:

http://hpaste.org/fastcgi/hpaste.fcgi/view?id=11012#a11019

Comments and suggestions are very welcome!

Alex.

Alex S. said...

So I went ahead and wrote the code. It's on hpaste under http://hpaste.org/fastcgi/hpaste.fcgi/view?id=11012#a11019

Try to run it, seems to work for me!

Alex.

gereeter said...

scc,

That's an interesting little duality. The product type can be a comonad with no restrictions, but only a monad if its "extra" type is a monoid. Similarly, the reader type (or function type, whichever you prefer) can be a monad with no restrictions, but can only be a comonad if the "extra" type is a monoid.

Also, you were expressing a lack of confidence in your definitions, so let's verify the comonad laws:

1.

coreturn . cojoin
= \x -> coreturn (cojoin x)
= \x -> cojoin x mempty
= \x -> (curry (x . uncurry mappend)) mempty
= \x y -> (x . uncurry mappend) (mempty, y)
= \x y -> x (uncurry mappend (mempty, y))
= \x y -> x (mappend mempty y)
= \x y -> x y
= \x -> x
= id

2.

fmap coreturn . cojoin
= \x -> fmap coreturn (cojoin x)
= \x -> coreturn . cojoin x
= \x y -> coreturn (cojoin x y)
= \x y -> cojoin x y mempty
= \x y -> curry (x . uncurry mappend) y mempty
= \x y -> (x . uncurry mappend) (y, mempty)
= \x y -> x (uncurry mappend (y, mempty))
= \x y -> x (mappend y mempty)
= \x y -> x y
= \x -> x
= id

3.

cojoin . cojoin
= \x -> cojoin (cojoin x)
= \x -> cojoin (curry (x . uncurry mappend))
= \x -> curry (curry (x . uncurry mappend) . uncurry mappend)
= \x y z -> curry (curry (x . uncurry mappend) . uncurry mappend) y z
= \x y z -> (curry (x . uncurry mappend) . uncurry mappend) (y, z)
= \x y z -> curry (x . uncurry mappend) (uncurry mappend (y, z))
= \x y z -> curry (x . uncurry mappend) (mappend y z)
= \x y z w -> curry (x . uncurry mappend) (mappend y z) w
= \x y z w -> (x . uncurry mappend) (mappend y z, w)
= \x y z w -> x (uncurry mappend (mappend y z, w))
= \x y z w -> x (mappend (mappend y z) w)
= \x y z w -> x (mappend y (mappend z w))
= \x y z w -> x (uncurry mappend (y, mappend z w))
= \x y z w -> (x . uncurry mappend) (y, mappend z w)
= \x y z w -> curry (x . uncurry mappend) y (mappend z w)
= \x y z w -> cojoin x y (mappend z w)
= \x y z w -> cojoin x y (uncurry mappend (z, w))
= \x y z w -> (cojoin x y . uncurry mappend) (z, w)
= \x y z w -> curry (cojoin x y . uncurry mappend) z w
= \x y z w -> cojoin (cojoin x y) z w
= \x y -> cojoin (cojoin x y)
= \x y -> (cojoin . cojoin x) y
= \x -> cojoin . cojoin x
= \x -> fmap cojoin (cojoin x)
= \x -> (fmap cojoin . cojoin) x
= fmap cojoin . cojoin


Phew. That was long. But yes, your definitions are correct.

Jeremy List said...

Toroidal universes are a pretty simple extension (you could just use circular lists but that would lead to recalculating a lot of values)

> right (U a b []) = right (U [] b (reverse a))
> right (U a b (c:cs)) = U (b:a) c cs
> left (U [] b c) = left (U (reverse c) b [])
> left (U (a:as) b c) = U as a (b:c)

Jeremy List said...

Toroidal universes are a pretty simple extension (you could just use circular lists but that would lead to recalculating a lot of values)

> right (U a b []) = right (U [] b (reverse a))
> right (U a b (c:cs)) = U (b:a) c cs
> left (U [] b c) = left (U (reverse c) b [])
> left (U (a:as) b c) = U as a (b:c)

Blog Archive