Cellular automatons are great at illustrating the power of functional programming. They also make for very pretty Clojure programs.
This article talks about Rule 30, what's so interesting about it, and how to implement it in Clojure.
There will be a second part to this article as well, covering ways in which we can make this automaton go faster using Clojure's concurrency and parallel programming tools.
But first, a pretty picture:
Oooooh!.. Now that we have that out of the way, on with the article: A Whole Lotta Automata.
Rule 30
(run-steps 10 stepper-3 rule30 state)
([0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 1 1 0 0 1 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 1 1 0 1 1 1 1 0 0 0 0 0 0]
[0 0 0 0 0 0 1 1 0 0 1 0 0 0 1 0 0 0 0 0]
[0 0 0 0 0 1 1 0 1 1 1 1 0 1 1 1 0 0 0 0]
[0 0 0 0 1 1 0 0 1 0 0 0 0 1 0 0 1 0 0 0]
[0 0 0 1 1 0 1 1 1 1 0 0 1 1 1 1 1 1 0 0]
[0 0 1 1 0 0 1 0 0 0 1 1 1 0 0 0 0 0 1 0]
[0 1 1 0 1 1 1 1 0 1 1 0 0 1 0 0 0 1 1 1])
If all goes well, the lines you see above should look something like the picture you see below:
Okay, so the ones and zeros don't convey quite the contrast that this picture does... You just have to look at them from far away.
What's going on here? It's a two one dimensional cellular automaton, where the state exists as a horizontal line and time exists on the vertical Y axis. The function scans across the horizontal line, creating a new line beneath it. Whether a block is black or white (on or off, 1 or 0, whatever you like) depends on the blocks above it -- specifically the block to it's Northeast, it's North and it's Northwest.
Observe, at the top of the picture above, the mappings that specify Rule 30. Like the one on the far left, for instance: If the function comes across three black boxes, it places a white box below them. It does this until it reaches the end of the line and then starts over again, using the new line as the input to the function that creates another subsequent line, and so forth.
Now look at another picture of Rule 30 from really far away:
Unbounded Complexity *eagle screams*
What's interesting is how such complexity can grow from such a simple set of rules. For a three bit sequence there are eight possible states:
(defn mkrule [v1 v2 v3 v4 v5 v6 v7 v8]
{ [1 1 1] v1
[1 1 0] v2
[1 0 1] v3
[1 0 0] v4
[0 1 1] v5
[0 1 0] v6
[0 0 1] v7
[0 0 0] v8 })
If we map each possible state to an output of either 1 or 0, we have 256 possible mappings. Well, it turns out that half of those mappings are simply mirror images of the other half of those mappings, resulting in only 128 possible unique mappings.
There are a few of those rules, however, that exhibit what some call "unbounded complexity," such as Rule 30 above. The complexity of our output is bounded only by the width of our automaton and the length of time we have to compute our new lines. If our automaton was only 10 cells wide (like in my first example above), then it would not take long for the output to begin repeating itself. Perhaps we might see a 100 or 1000 generation cycle that just loops infinitely.
Given unlimited time and space, though, it is believed that Rule 30 grows in complexity forever -- meaning it can constantly produce new patterns that are unique to the whole output. That such a simple function -- such a minimal amount of initial information -- can produce infinite amounts of information is sure to stoke the philosophical embers in anyone's camp fire.
The definition for Rule 30 in our program is simply:
(def rule30 (mkrule 0 0 0 1 1 1 1 0))
The initial conditions are:
(def state (assoc (vec (take 20 (iterate identity 0))) 10 1))
This produces:
[0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0]
... A vector of zeros with a one in the center.
This isn't a whole lot of data to start off with, but we end up with a very rich output. The magic is in the logic of the functions.
How it works
So, we need to provide our function run-steps with the number of time steps we want to run, the-steps; the function we use to transform the current state, the-stepper; the rule that we will be supplying to your transformer function, the-rule; and the actual state we will use for our initial conditions, the-state.
(defn run-steps [the-steps the-stepper the-rule the-state]
(take the-steps
(iterate #(the-stepper the-rule %) the-state)))
iterate simply feeds the output of the-stepper back into its input, dropping a copy off with us each time, over and over until we stop asking.
the-stepper gets a little bit interesting because we have to "wrap-around" the state. On the very first cell, the cell to our Northwest is actually the one at the other end of the line. And when we get to the end of the line the cell to our Northeast is the one at the beginning of the line. So, we want a stepper function that not only returns a consecutive sequence of vectors that each have three elements, incrementing by an index of one, but also that has the behavior of wrap-around when we are at either end of the line.
In other words, for the sequence [0 1 2 3 4 5 6 7 8 9], we want the following output:
[9 0 1] [0 1 2] [1 2 3] ... [7 8 9] [8 9 0]
Of course, our sequences will have strictly ones and zeros. Notice that, although this looks like a partitioning operation, we are only indexing by one.
So, on to our stepper function:
(defn stepper-3 [the-rule the-state]
(->> (range (count the-state))
(map #(the-rule (wrap-sub the-state % 3)))
(vec)))
Again, the tricky wrapping stuff is offloaded to another helper function, wrap-sub. stepper-3 just creates a range from zero to however long the sequence is: (range (count [1 0 1])) => (0 1 2) ... and maps the range values into the wrap-sub function (along with the state), handing the resulting three-element vector to the-rule, and then finally returns the whole result as a vector.
The wrap-sub function requires a collection, the-state; a location so it knows where to grab from in the collection, the-index; and the number of elements you would like it to grab, the-width:
(defn wrap-sub [the-state the-index the-width]
(let [left-index (- (int (/ the-width 2)))
right-index (+ the-width left-index)]
(->> (range left-index right-index)
(map #(wrap-nth the-state the-index %) )
(vec))))
the-width, or the number of elements we want to grab are always three (for the purposes of this kind of automaton). In other automatons or for other reasons (as you will see below), we may want to grab five or 15 or 1000. Either way, wrap-sub assumes that you want your index to be at the center of the group of elements it grabs. So, using again our sequence from above, suppose we want the sixth element of [0 1 2 3 4 5 6 7 8 9] with a grouping width of five. The group returned will be [4 5 6 7 8], with six in the center.
The first let statement halves the-width and makes it a negative int, giving us our left-index. The second let statement returns the right-index for the positive half of our grouping. Once again we create a range, this time starting with a negative number, and map that range into another helper function, returning a vector.
(defn wrap-nth [the-state absolute-index relative-index]
(nth the-state
(rem (+ absolute-index size relative-index)
(count the-state))))
Using our example sequence from above, with an index of six and a desired group size of five:
(wrap-nth [0 1 2 3 4 5 6 7 8 9] 6 -2)
This adds up the size (10), the absolute-index (6) and the relative-index (-2) to get 14.
Then it returns the remainder of 14 divided by the size of the collection (10), resulting in four. It then gets and returns the 4th element in the collection.
The next relative-index that gets mapped in from wrap-sub:
(wrap-nth [0 1 2 3 4 5 6 7 8 9] 6 -1)
... results in the 5th element. And so forth.
Is this the most efficient way to get the result we want? Probably not. There is probably some incantation of partition-by or some other new function that makes this much easier. Hopefully someone will add some wisdom in the comments. For the time being though, the method used above works.
Automata in Practice
There are a lot of reasons why you would want to learn to program in this way. This program is purely functional. We could have just as easily created a stateful grid of values and 'banged on them in place.' Instead, this program pipelines the state through a number of functions and produces a result. An automaton epitomizes the functional model because the same input should _always_ produce the same output. When thinking about a program, a useful exercise is to see if you can think of it in terms of a cellular automaton.
Stephen Wolfram, the guy who invented this two one mensional two-state automaton, believes that the whole universe is a cellular automaton. The application of small programs like this are, indeed, seen throughout nature. Observe this sea-shell and notice the patterns it creates:
Ooooh!.. Pretty shell.
Some people believe that the most complex biological systems we see in nature are simply many layers of extremely simple programs, like Rule 30.
Cellular automatons have also been used in cryptography and Rule 30 has been used as a pseudo-random number generator. Cellular automatons have been used to simulate weather patterns, traffic patterns, social patterns like ant colonies, etc.
Rich Hickey's ant colony simulation, designed to teach concurrency modelling in Clojure, is also a cellular automaton.
The Program
Here is our Rule 30 program in it's entirety:
(ns Automata.core)
(def state (assoc (vec (take 20 (iterate identity 0))) 10 1))
(defn mkrule [v1 v2 v3 v4 v5 v6 v7 v8]
{ [1 1 1] v1
[1 1 0] v2
[1 0 1] v3
[1 0 0] v4
[0 1 1] v5
[0 1 0] v6
[0 0 1] v7
[0 0 0] v8 })
(def rule30 (mkrule 0 0 0 1 1 1 1 0))
(defn wrap-nth [the-state absolute-index relative-index]
(let [size (count the-state)]
(nth the-state
(rem (+ absolute-index size relative-index) size))))
(defn wrap-sub [the-state the-index the-width]
(let [left-index (- (int (/ the-width 2)))
right-index (+ the-width left-index)]
(->> (range left-index right-index)
(map #(wrap-nth the-state the-index %) )
(vec))))
(defn stepper-3 [the-rule the-state]
(->> (range (count the-state))
(map #(the-rule (wrap-sub the-state % 3)))
(vec)))
(defn run-steps [the-steps the-stepper the-rule the-state]
(take the-steps
(iterate #(the-stepper the-rule %) the-state)))
(run-steps 10 stepper-3 rule45 state)
And, if you're feeling daring, run these as well and let me know what your benchmarks are in the comments:
(def state (assoc (vec (take 1000000 (iterate identity 0))) 50000 1))
(time (take 20 (drop 49990 (last (run-steps 10 stepper-3 rule30 state)))))
On my quad-core i7, clojure 1.2, I get about 22 seconds:
"Elapsed time: 22114.679595 msecs"
(0 1 1 0 1 1 1 1 0 1 1 0 0 1 0 0 0 1 1 1)
Also, feel free to optimize my code. I'd love to see this go faster!
Concurrency and Parallelism
In my next blog post, I'm going to explore some of Clojure's concurrency and parallel programming tools to see if I can get this automaton to go any faster.
Here is some code I'm going to be working on if you want to get a head start:
;; Expects pre-wrapped/ear-muffed sequences
(defn chunk-stepper-3 [the-rule the-state]
(->> (drop 1 (drop-last (range (count the-state))))
(map #(the-rule (wrap-sub the-state % 3)))
(vec)))
;; Don't look at it unless you want to make it prettier ;)
;; Makes chunks of size given and conjs smaller chunk if
;; there are any remaining elements
(defn chunkify [the-state chunk-size]
(let [window (+ 2 chunk-size)
first-index (- (int (/ window 2)) 1)
the-range (range first-index
(* chunk-size
(int (/ (count the-state) chunk-size)))
chunk-size)
size (count the-state)
rem-size (rem size chunk-size)
tail-sub (if (not= 0 rem-size)
(let [range-tail
(vec (range (- size rem-size) size))
tail-size (count range-tail)
tail-index
(nth range-tail
(int (/ tail-size 2)))]
(wrap-sub the-state tail-index
(+ tail-size 2))))]
(if tail-sub
(conj (vec (pmap #(wrap-sub the-state % window)
the-range))
tail-sub)
(vec (pmap #(wrap-sub the-state % window)
the-range)))))
;; Chunks and un-chunks using the stepper
(defn chunkerator [the-stepper the-rule the-state chunk-size]
(vec (flatten (pmap #(the-stepper the-rule %) (chunkify the-state chunk-size)))))
;; iterates the chunkerator
(defn run-chunks [the-steps the-stepper the-rule the-state chunk-size]
(take the-steps (iterate #(chunkerator the-stepper the-rule % chunk-size) the-state)))
(time (take 20 (drop 49990 (last (run-chunks 10 chunk-stepper-3 rule30 state 1000)))))
"Elapsed time: 19441.61972 msecs"
(0 1 1 0 1 1 1 1 0 1 1 0 0 1 0 0 0 1 1 1)
As you can see pmap only wins me about 2 seconds (10%) over map. That's with a chunk size of 1000. In my next post I am going to explore the possibility of delegating out chunks to agents. Does anyone have any other ideas?
Conclusion
So that's it. Clojure obviously makes for a very expressive language when working with cellular automatons. Stay tuned for part two on concurrency and parallelism on automatons.
Pretty expression:
((((((((((((defn λ☯☮[]λ☯☮))))))))))))
