Simulating Experiences: Excursions in Programming

"Modelling Nature with Cellular Automata"
Part I.  One-Dimensional CA's

Mathematica in Education
Vol.2 No.3
(Spring 1993)

by
Richard J. Gaylord
Department of Materials Science
University of Illinois at Urbana-Champaign
gaylord@ux1.cso.uiuc.edu

Introduction

A cellular automaton (CA) is a discrete dynamic system of lattice sites
that evolve in discrete time steps as each site assumes a value from a
finite set of values according to local and uniform rules determined by the
values of a neighborhood of sites around it.

Stephen Wolfram [Wolfram, 1984] has expressed the view that "any physical
process can be described by an algorithm and therefore can be represented
by a computational process" and that "cellular automata can be regarded as
computers and as models of physical systems".

One cellular automaton which was discussed earlier [Gaylord, 1992], is the
sandpile model which has been used to describe self-organized criticality
(SOC) and catastrophes in complex systems. Now I want to take a  look at
the two most famous cellular automata. In Part I, I'll discuss
one-dimensional cellular automata In Part II, I'll turn to the
two-dimensional cellular automaton known as the "Game of Life".

One-Dimensional Cellular Automata

In landmark investigations [Wolfram, 1993] which introduced computation as
a primary tool or methodology for scientific investigation, equal in
importance to both theory and experiment, Wolfram studied the evolution of
one-dimensional cellular automata (1D CA).

A 1D CA starts with a row of lattice sites (or cells) each of which has
some value. The row either consists of a single seed site with a non-zero
value, located in the center of the row, or is a disordered state in which
seed sites are randomly distributed in the row.

The CA "evolves" by producing successive rows of sites, each new row
representing a time step. The value of each site in a new row is determined
by applying a specific rule to the site in the preceding row which lies
directly above the new site.

The rule is homogeneous, which means that the same rule is used in
determining the value of each new site, and the rule is local, which means
that it is based the values of a "neighborhood" of sites, consisting of the
site and/or its nearest-neighbor (nn) sites. The rule is applied to all of
the sites in a row simultaneously. For a finite CA, the nn sites of the
first and last sites in a row are taken from the other end of the row (eg.,
the value of the next nn site on the left-hand side of the first site is
given by the last site). The use of this "wrap-around" reflecting boundary
condition essentially turns the line into a circle.

As the updating of the current state of the CA proceeds and new rows
appear, a pattern emerges.

Wolfram recognized that 1D CA's evolve into patterns falling into one of
four general classes (analogies to physical and computing systems are noted
in parentheses):

#1 - homogeneous state [analogous to limit points] 
#2 - periodic structures [analogous to limit cycles]
#3 - chaotic pattern [analogous to strange attractors]
#4 - complex, localized structures [possibly analogous to universal computers]

The analogies of the behaviors of CA's in the first three classes to the
behaviors of physical systems is quite interesting. Even more intriguing is
the behavior of a CA in the fourth class where there is no simple algorithm
for predicting behavior and it is necessary to actually run the CA
evolution simulation to determine the outcome.

The fact that 1D CA's show such great "complexity" in their behaviors has
rather profound implications (eg., some people have postulated that the
universe itself, is a cellular automaton). Here, I will develop a program
for creating 1D CA's which can be used to join in the exploration of these
fascinating systems.

CA Rules

I can write down the rules that govern the evolution of lattice sites based
on the number of lattice sites that interact and the values that a lattice
site can have.

I'll first calculate the total number of rules for a given set of k and r
values. When 2r+1 sites interact (ie., a given site and its r nn sites on
each side) and each site can have k values (the integer values 0 through
k-1), the total number of possible rules is given by:

totalrules[k_, r_] := k^(k^(2r+1))

Applying this formula to some specific examples:

(1) When each site can have two values, 0 or 1, and sites evolve
independently of the state of their neighbors, there are four possible
rules,

totalrules[2, 0]

which I can write as

{{{1} -> 0, {0} -> 0}, 
 {{1} -> 0, {0} -> 1}, 
 {{1} -> 1, {0} -> 0}, 
 {{1} -> 1, {0} -> 1}}

Let me explain the notation being used here.  The outermost list is a
nested list consisting of lists of transformation rules for determining the
value of a new site. Each inner list (eg., {{1} -> 1, {0} -> 1}) is called
a rule set.  The list on the left-hand side of a transformation rule in a
rule set is a tuplet whose elements are the values of the neighborhood of a
site (the middle element of the tuplet is the value of the site). The
integer on the right-hand side of the transformation rule is the value of
the new site created directly beneath the site during a time step of CA
evolution (eg., the transformation rule {0,1,1} -> 0 says that a new site
created beneath a site having a value 1 and a left next nn with value 0 and
a right next nn with value 1, will have a value 0).

(2) As the number of nearest-neighbor sites involved in the interaction
increases, the number of rules rapidly increases. 
When three (2r+1 = 3) sites (ie., a given site and its next nn sites on
each side (vwx)) interact and each site can have k = 2 values, there are
the following possible neighborhood values:

{{1, 1, 1}, {1, 1, 0}, {1, 0, 1},  {1, 0, 0}, {0, 1, 1}, 
 {0, 1, 0}, {0, 0, 1}, {0, 0, 0}}

Thus, there are eight or (2 2 2) = 2^3 = 8 possible triplets that can be
formed.  Each of these triplets can result a site value of 1 or 0. For
example, one rule set is:

{{1,1,1} -> 1, {1,1,0} -> 1, {1,0,1} -> 1,  {1,0,0} -> 0, 
  {0,1,1} -> 1, {0,1,0} -> 0, {0,0,1} -> 0, {0,0,0} -> 1}

Since each triplet can result in either a 1 or 0, there are a total of 
(2 2 2 2 2 2 2 2) = 2^8 = 256  rule sets.

totalrules[2, 1]

(3) When each site can have two values and interacts with its two nn sites
on each side (uvwxy) there are 2^5 = 32 possible quintuplets and 2^32 =
4,294,967,296 rules sets!

totalrules[2, 2]

(4) Increasing the number of values that a site can have beyond two 
increases the number of rule sets even more dramatically. When each site
can have three values and sites interact with their next nn site on each
side, there are 3^(3^3) =  7,625,597,484,987 rules sets.

totalrules[3, 1]

I would like to have a shorthand notation for expressing a rule set. I can
illustrate how this can be done using example 2.

{{1,1,1} -> 1, {1,1,0} -> 1, {1,0,1} -> 1,  
  {1,0,0} -> 0, {0,1,1} -> 1, {0,1,0} -> 0, 
  {0,0,1} -> 0, {0,0,0} -> 1}

I take the right-hand sides of the transformation rules in this rule set
and form the binary number 11101001. The base ten equivalent of this binary
number is 233.

2^^11101001

So, I can say that the rule number 233 corresponds to the rule set

{{1,1,1} -> 1, {1,1,0} -> 1, {1,0,1} -> 1,  
  {1,0,0} -> 0, {0,1,1} -> 1, {0,1,0} -> 0, 
  {0,0,1} -> 0, {0,0,0} -> 1}

Each of the k^(k^(2r+1)) rules in a 1D CA can be given its own rule number
in this way, using a base k number to represent the right-hand sides of the
transformation rules in the rule set and taking its base 10 equivalent.

Playing By The Rules

I now turn to creating a rule set for any of the k^(k^(2r+1)) rules in a 1D CA.

I'll start by looking at the rule number 233 example used above. I can
easily create the required rule set by hand:

MapThread[Rule, 
      {{{1,1,1}, {1,1,0}, {1,0,1}, {1,0,0}, 
       {0,1,1}, {0,1,0}, {0,0,1},{0,0,0}},
       {1, 1, 1, 0, 1, 0, 0, 1}}]

Now I'll write a general program that generates the rule set for a given
value of k, r and the rule number w.

I first create the list of the number of tuplets that are involved when a
1D CA has 2r+1 nearest neightbor interactions and k site values. There are
(as usual) several ways to create the tuplets. The following function
appears to be the fastest, by far:

tuplets[k_, r_] := 
      Reverse[Flatten[Apply[Outer,
          Prepend[Table[Range[0,k-1],{(2r+1)}],List]],2r]];

Trying the definition out on the k = 2, r = 1 CA,

tuplets[2, 1]

I next create the list of the right hand sides of the transformation rule
in the rule set.  This list is called the update value list.

I'll need to convert the rule number w into a list whose elements represent
the number in base k. There is a built-in function, IntegerDigits, which
does this. For example,

IntegerDigits[233, 2]

I now have everything I need: a list of the tuplets, a list of update
values, and the MapThread function. However, there's a slight complication.
Consider rule 24 for a k = 2, r = 1 CA.

IntegerDigits[24, 2]

Length[tuplets[2, 1]] == Length[IntegerDigits[24, 2]]

I can't thread together lists having different lengths (and taking
2^^10000000 indicates that half of the rule numbers, those below 128, have
less than the prerequisite eight elements for threading with the list of
eight triplets for the k = 2, r = 1 CA). What is needed is to add 0's to
the front of IntegerDigits[w, k] to bring the list length up to k^(2r+1). I
can do this using

fillcode[k_, r_, w_] :=
  Join[Table[0,{k^(2r+1)-Length[#]}],#]&[IntegerDigits[w, k]]

which works for all rule numbers

fillcode[2, 1, 24]

fillcode[2, 1, 233]

(note: I am taking advantage here of the fact that Table returns an empty
list when the iterator is zero or negative).

Overall then, I can create a rule set with the following formula:

ruleset1[k_, r_, w_] :=
      Module[{fillcode1, tuplets1},
           fillcode1[x_, y_, z_] :=
                Join[Table[0,{x^(2y+1)-Length[#]}],#]&[IntegerDigits[z, x]];
           tuplets1[p_, s_] :=
                 Reverse[Flatten[Apply[Outer,
                         Prepend[Table[Range[0,p-1],{(2s+1)}],List]],2s]];
            MapThread[Rule, {tuplets1[k, r], fillcode1[k, r, w]}]
                 ]

or in the form of a totally inscrutable one-liner,

ruleset2[k_, r_, w_] :=
    MapThread[Rule, 
       {Reverse[Flatten[Apply[Outer,
           Prepend[Table[Range[0,k-1],{(2r+1)}],List]],2r]],
               Join[Table[0,{k^(2r+1)-Length[#]}],#]&[IntegerDigits[w, k]]}]

Just to check that this definition works, I'll generate the list of rule
sets for the simplest case of non-interacting sites having values of 1 and
0 which I discussed at the beginning of the previous section.

{ruleset2[2, 0, 0], ruleset2[2, 0, 1], ruleset2[2, 0, 2], ruleset2[2, 0, 3]}

(note that if I try to calculate ruleset2[2, 0, 4], I get an error message.)

ruleset2[2,0,4]

Now, I'm ready to develop the program for evolving a one-dimensional
cellular automaton.

The 1D CA Program

I'll begin by defining the symbols that appear in the program:

m = number of sites per row
n = maximum number of rows (generations)
w = rule number (this number must be less than k^(k^(2r+1)))
k = number of values that a site can assume 
2r+1 = number of interacting sites

I'll first explain the individual parts that go into the program and then
put them together.

The CA starts from row, which is either a disordered state in which there
are randomly distributed seed sites

row = Table[Random[Integer, {0, k-1}],{m}]

or a single seed site

row = ReplacePart[Table[0, {m}],
                  Random[Integer, {0, k-1}],Ceiling[m/2]]

(note: the seed sites here have randomly determined values. It is also
possible to specify the seed site values).

I'll use CArows to evolve the CA.

CArows = FixedPointList[(Partition[
      Join[Take[#,-r],#,Take[#,r]],(2r+1),1]/.ruleset)&, 
       row, n]

FixedPointList applies an anonymous function to the row of seed site(s) to
produce a new row of sites, and then applies the anonymous function to the
new row of sites to produce another row of sites and so on until there is
no change in the values in successive rows or to a maximum of n rows,
whichever occurs first. The anonymous function

(Partition[Join[Take[#,-r],#,Take[#,r]],(2r+1),1]/.ruleset)&

performs the following sequence of operations:

(1) a reflecting boundary condition is applied to a row of m sites (denoted
here by {...}) to create a list of m+2r values.

Join[Take[#,-r],#,Take[#,r]]& [{...}]

(2) the above list is partitioned to form a list of m tuplets of size 2r+1
(the elements of successive tuplets will overlap when r > 0) containing the
values of the neighborhoods of each of the m sites in the row.

Partition[%,(2r+1),1]

(3) a new row is created using the rule set to determine the value of the
new sites beneath the old sites.

% /. ruleset

Putting all of the code fragments together and adding some code to: (1)
print out the number of generations produced by the CA, and (2) create a
two-dimensional RasterArray graphics with randomly chosen RGBColors
representing different site values, I have

CAdisordered[w_, k_, r_, m_, n_] :=
 Module[{ruleset, row, CArows},
 ruleset =
    MapThread[Rule, 
       {Reverse[Flatten[Apply[Outer,
            Prepend[Table[Range[0,k-1],{(2r+1)}],List]],2r]],
               Join[Table[0,{k^(2r+1)-Length[#]}],#]&[IntegerDigits[w, k]]}];
    row = Table[Random[Integer,{0, k-1}],{m}];
    CArows = FixedPointList[(Partition[
        Join[Take[#,-r],#,Take[#,r]],(2r+1),1]/.ruleset)&, row, n];
    Print[Length[CArows]];
    Show[Graphics[RasterArray[
      Reverse[CArows]/.Thread[Range[0, k-1]->
            (Map[RGBColor, Table[Random[],{k},{3}]] /.
                 RGBColor[{x__}]->RGBColor[x]) ]]], AspectRatio->Automatic]
             ]

Running the Program

We show a few illustrative k = 2, r = 1 1D CA's below.

GaylordFig.1.eps

SeedRandom[6] 
CAdisordered[78, 2, 1, 50, 15];

GaylordFig.2.eps

SeedRandom[9]
CAdisordered[558430, 4, 1, 50, 15];

CA Animation

While a static figure can be used to detect CA patterns, it would be nice
to watch the evolutionary process unfolding. To do this, an animation can
be created.

In creating an animation using RasterArray graphics, the main thing I have
to take into account is that the rectangles, representing site values, that
re-appear in successive figures in the series that constitute the flip
cards of the animation must be the same size and color. This requires two
things:

(1) each figure must have the same number of rows.  This is accomplished by
adding enough rows of 0's to the CArows used in each figure to bring the
total number of lines in the figure up to the Length of the final CA (so I
have to complete the entire CA computation before I create the graphics).

(2) the transformation rule for converting site values to RGBColors must be
applied after all of the site values in all of the CA figures in the
animation are calculated and before the graphics command is applied.

The CA animation program is:

CAdisorderedAnimation[w_, k_, r_, m_, n_] :=
 Module[{ruleset, row},
    ruleset = 
       MapThread[Rule, 
           {Reverse[Flatten[Apply[Outer,
             Prepend[Table[Range[0,k-1],{(2r+1)}],List]],2r]],
               Join[Table[0,{k^(2r+1)-Length[#]}],#]&[IntegerDigits[w, k]]}];
    row = Table[Random[Integer,{0, k-1}],{m}];
    CArows = FixedPointList[(Partition[
         Join[Take[#,-r],#,Take[#,r]],(2r+1),1]/.ruleset)&, row, n];
    Print[Length[CArows]];
    Map[Show[Graphics[RasterArray[#]],AspectRatio->Automatic]&,
     Map[Join[Table[Table[0,{m}],{Length[CArows]-#}], 
              Reverse[Take[CArows,#]]]&,
                             Range[Length[CArows]] ]/.Thread[Range[0, k-1]->
             (Map[RGBColor, Table[Random[],{k},{3}]]/.
                         RGBColor[{x__}]->RGBColor[x])] ];
           ]

SeedRandom[9] 
CAdisorderedAnimation[1111111113,4,1,100,20]

(note: after creating the animation in a Mathematica notebook, double-click
on the bracket enclosing all of the graphics cells (this bracket is created
when Automatic Grouping has been selected in the Cell menu). This will
close up the graphics cells, leaving only the first cell exposed.
Double-clicking on the exposed cell starts the animation. The speed of the
animation can be adjusted using the number keys (1 is the slowest speed and
9 is the fastest speed)).

Conclusion

Happy pattern hunting. We will next look (finally) at various programs for
the Game of Life [unlike politicians, I always fufill my promises
eventually].

References

Wolfram, Stephen. 1984. "Computer Software in Science and Mathematics",
Scientific American, 251(Sept.), pp. 188-203

Wolfram, Stephen. 1993, "Wolfram on Cellular Automata and Complexity: A
Reprint Volume", Addison-Wesley, Reading, Mass. (the original research
articles in this collection, which is due to be published in September
1993, are largely accessible to the non-specialist).

Gaylord, Richard J. 1992. "Catastrophes in Complex Systems", Mathematica in
Education, 2(1), pp. 19-23.

Algorithmic Housekeeping

Below are an improved version of one program (Spreading) and a correction
to another program (Percolation Clustering) that I  presented previously.

Correction to the Percolation Clustering Algorithm:

The clusters in ufinal in Out[1] in the percolation clustering article
[MiE, 2(2), 21-24 (1993)] are supposed to be consecutively numbered but
aren't. This is because the code for relabelling in the clusterLabel
program is wrong.

relabelrules2 = Thread[new->correctlabels];
ufinal = ul /. relabelrules2;

should be replaced by

relabelrules2 = 
   Thread[Union[correctlabels]->
      Range[Length[Union[correctlabels]]]];
ufinal = ul //. relabelrules1 /. relabelrules2;

Sorry about that (I might note that the incorrect code was run a number
times without producing an incorrect output; this demonstrates the need to
perform a great many runs to ensure, as much as possible, that a program
works in all cases).

Improved Spreading Algorithm

By changing the way in which the perimeter site with the lowest random
number is selected for the invasion percolation case in the spreading
program [MiE, 1(4), 17-19 (1992)], the program becomes more efficient (and
elegant). This requires a number of changes in the original program and so
I simply give the complete improved program below:

spreading[n_Integer,p_Real:-1.0]:= 
  Module[
   {cluster = {{0, 0}}, 
    choices = {{1, 0},{0, 1},{-1, 0},{0, -1}},  
    per = Transpose[{Table[Random[],{4}], {{1,0},{0,1},{-1,0},{0,-1}}}], 
    perTocluster, nn, new, newper, testsite},  
     perTocluster := 
        (Module[{nn, new, newper},
            AppendTo[cluster, Flatten[#, 1][[2]] ];
        per = Complement[per, #];
        nn = Map[Function[x, x + Last[cluster]], choices]; 
        new = Complement[nn, cluster, Transpose[per][[2]] ];
        newper = Transpose[{Table[Random[],{Length[new]}],new}];
        per = Join[per,newper]])&;
     While[Length[cluster] < n + 1 && Length[per] >0,
               testsite = Random[Integer,{1, Length[per]}];
           Which[p < 0, perTocluster[First[Sort[per]]],
                  per[[testsite, 2]] <= p, perTocluster[{per[[testsite]]}], 
                  True, per = Complement[per, {per[[testsite]]}] ] 
                 ];
     Show[Graphics[{RGBColor[0.3,0.3,0.6],
                    Map[(Rectangle[#-{0.5,0.5},#+{0.5,0.5}])&,cluster]}],
                AspectRatio->1, 
                PlotRange->Map[({Min[#], Max[#]})&,Transpose[cluster]]] 
     ]