A Cellular Automata
Programming Toolkit

Rule-based programming with Mathematica provides efficient templates for creating cellular automata, including the famous Game of Life.

by Richard J. Gaylord and Kazume Nishidate

Cellular automata can be used to model a broad range of natural phenomena [Gaylord and Nishidate 1996]. Here we exploit Mathematica's rule-based programming paradigm and its APL-like ability to manipulate all of the elements of a matrix simultaneously. The result is a simple approach to developing CA programs that minimizes the programming needed for CA creation.

A scientific definition of a cellular automaton

Scientists generally define a cellular automaton (CA) to be a discrete dynamical system, where space, time, and the states of the system are all discrete and have the following properties:

- Space is represented by a regular lattice in one, two, or three dimensions.

- Each site, or cell, in (or on) the CA lattice can be in one of a finite number of states. The states are represented by integer number values.

Note: A variant of the CA, known as a coupled map-lattice (CML) or cell-dynamic scheme (CDS), uses continuous lattice site values. In general, lattice site values can be integers, reals, or symbols, as well as tuplets (lists) with elements of any type, where the choice of value is determined by the kinds of system components that are being represented (eg., when objects in the system have more than one changeable attribute or property, the values of these quantities are grouped in a n-tuple).

- The CA system evolves over a succession of time steps. The values of all of the sites in the lattice are updated synchronously in each time step.

Note: While the values of all the sites are updated in each time step, the value of a site need not change in a given time step.

- Site values are updated using a set of rules known as a lookup table  which take the values of the site and its neighboring sites into account.

A computational definition of a cellular automaton

Since a CA model is always expressed as an algorithm and is invariably implemented as a computer program, it's useful to give a computational definition of a CA.

A CA is a computer program in which the following computations are performed, in order:

- A matrix is created with specific element values (integer, real, symbol, or list).

- A function, or a set of functions, that can be used to change the value of a matrix element, based on the values of the element and elements in its vicinity, is defined.

- The function is repeatedly applied to the matrix, each time  simultaneously changing the values of all of the matrix elements.

The CA lattice

We usually work with two-dimensional cellular automata, employing a rectangular lattice. A simple square lattice is shown below.

Table["+", {5}, {5}] // MatrixForm

Note: The use of a rectangular lattice does not limit us to CA's based on this type of lattice. CA's based on other two-dimensional lattice systems (eg., trigonal and hexagonal) can also be handled with this lattice [Gaylord and Nishidate 1996].

Each lattice site has a value. Here is a random Boolean lattice where sites have values of 0 or 1.

(randomBoolean = Table[Random[Integer], {5}, {5}])//MatrixForm

A graphical representation of a CA is especially useful for observing patterns of structure and behavior that occur during the evolution of the CA. One way to visualize a CA lattice configuration is by using a RasterArray graphics where each lattice site value is transformed into an RGBColor.

The RasterArray representation of the above randomBoolean CA lattice is shown here:

Show[ Graphics[ RasterArray[Reverse[randomBoolean] /. {1 -> RGBColor[1, 0, 0], 
0 -> RGBColor[0, 1, 0]}], AspectRatio -> Automatic, GridLines -> {Range[0, 5], Range[0, 5]}] ]

(*Joe, can this be rendered as black and white*)

Note: The Reverse and GridLines commands were used to create the above figure in order to make it easier to see the correspondence between the sites in the randomBoolean matrix and the cells in the RasterArray graphic.

The neighborhood

In updating the values of a site in a CA lattice, it is necessary to consider the site's value and the values of the sites in its vicinity. The sites considered in a particular CA will depend on what interactions are being modeled in the system. 

The two most commonly used CA neighborhoods are shown below:

The von Neumann neighborhood

N  
          |  
      W - X - E 
          |  
          S

The set of five sites consisting of the X site and its four nearest neighbor sites lying north (above), east (right), south (below), and west (left) of X are known as the von Neumann neighborhood of X.

The Moore neighborhood

NW   N   NE
        \ | /
     W  - X -  E 
        / | \
     SW   S   SE

The set of nine sites consisting of the X site, the four nearest neighbor sites in its von Neumann neighborhood, and the four nearest neighbor sites lying northeast, southeast,  southwest, and northwest of X are known as the Moore neighborhood of X.

Another neighborhood which we often use (eg., to model diffusive processes) is shown below

The MvonN neighborhood

Nn
          |          
     NW   N   NE
        \ | /
Ww - W  - X -  E - Ee 
        / | \
     SW   S   SE
          |
          Ss

This neighborhood consists of the nine sites in the Moore neighborhood of X, and the four next nearest neighbors lying in the north, east, south, and west direction. For lack of a better name, we refer to this as the MvonN neighborhood of X.

Note: While the neighborhoods shown above are regularly-shaped, CA's can use irregularly-shaped neighborhoods [Gaylord and Nishidate 1996].

To illustrate which sites lie in the Moore neighborhood of a lattice site, the Moore neighborhoods of the sites on a small rectangular  lattice are illustrated below.

(Lat = Partition[Range[25], 5]) // MatrixForm

Wrap = Join[{Last[#]}, #, {First[#]}]&;
(LatNeighborhoods = Partition[Wrap[Map[Wrap, Lat]], {3, 3}, {1, 1}]) // MatrixForm

Looking the outermost neighborhoods in LatNeighborhoods, the nearest neighbors of sites lying at various locations along the borders of Lat can be identified as follows:

- the nearest neighbor west of a site on the left border is the site in the same row on the right border.

- the nearest neighbor northwest of a site on the left border is the site in the previous row on the right border.

- the nearest neighbor southwest of a site on the left border is the site in the following row on the right border.

- the nearest neighbor east of a site on the right border is the site in the same row on the left border. 

- the nearest neighbor northeast of a site on the right border is the site in the previous row on the left border.

- the nearest neighbor southeast of a site on the right border is the site in the following row on the left border.

- the nearest neighbor north of a site on the top border is the site in the same column on the bottom border.

- the nearest neighbor northwest of a site on the top border is the site in the previous column on the bottom border.

- the nearest neighbor northeast of a site on the top border is the site in the following column on the bottom border.

- the nearest neighbor below a site on the bottom border is the site in the same column on the top border. 

- the nearest neighbor southwest of a site on the bottom border is the site in the previous column on the top border.

- the nearest neighbor southeast of a site on the bottom border is the site in the following column on the top border.

- the nearest neighbor northwest of the site in the top left corner is the site in the bottom right corner.

- the nearest neighbor northeast of the site in the top right corner is the site in the bottom left corner.

- the nearest neighbor southwest of the site in the bottom left corner is the site in the top right corner.

- the nearest neighbor southeast of the site in the bottom right corner is the site in the top left corner.

Dealing with boundaries:

The boundary conditions that are used in a CA program depend on the specific physics of the system being modeled. Traditionally, the values of the lattice sites in the neighborhood of a border site are different for different boundary conditions. For example, under periodic boundary conditions, border sites have sites on the opposing border as neighbors (a site on the left border has the site in the same row on the right border as its left neighbor). Border sites under absorbing boundary conditions have no off-lattice neighbors (a site on the left border has no left neighbor). 

In our computations, the values of the lattice sites in the neighborhood of a border site are computed as described in the previous section under all boundary conditions. Different boundary conditions are dealt with by using appropriate values for sites along the lattice border.

Periodic boundaries

This boundary, commonly used in CA models of infinite systems, is the default boundary condition.

Note: In a CA model of mobile objects moving in space, the wrap-around boundary condition results in  objects that leave the system at one border location and re-enter at the opposing border location.

Absorbing (open) and reflecting (closed) boundaries and sinks and sources

These boundary conditions are implemented by using a unique value for border sites, different from the values used for interior sites, and writing rules that specify what happens to the value of a site in the neighborhood of a border site and that leave the value of a border site unchanged.

Note: By using unique values for border sites, rules can be written for a CA model of mobile objects moving in space, which take an object occupying a site adjacent to a border site and freezes the object in place, turns the object around, or removes the object  from the system or, if the site is empty, puts an object on the site.

Moving boundaries

One way to avoid having to deal with specific boundary conditions is to adjust the size of the CA lattice on each time step so that nothing interesting happens on the border sites.

Update rules for the CA lattice sites

Rules that update the value of a lattice site, based on the values of the sites in the neighborhood of the site, have one of the following forms, depending on the type of neighborhood.

Rules that are applied to a site and the neighbor sites in its von Neumann neighborhood have the form:

update[site, N, E, S, W]

where the five arguments represent the value of the site and the values of the four nearest neighbors in the N, E, S, W directions.

Rules that are applied to a site and the neighbor sites in its Moore neighborhood have the form:

update[site, N, E, S, W, NE, SE, SW, NW]

where the nine arguments represent the value of the site, the values of the four nearest neighbors in the N, E, S, W directions, and the values of the four nearest neighbors in the NE, SE, SW, NW directions.

Rules that are applied to a site and the neighbor sites in its Moore neighborhood and the neighbor sites in its von Neumann neighborhood have the form:

update[site, N, E, S, W, NE, SE, SW, NW, Nn, Ee, Ss, Ww]

where the thirteen arguments represent the value of the site, the values of the four nearest neighbors in the N, E, S, W directions, the values of the four nearest neighbors in the NE, SE, SW, NW directions, and the values of its four next nearest neighbors in the N, E, S, W directions.

There is no formula for creating CA update rules.The update rules in some CA's may depend on the values of the neighborhood sites while the update rules in other CA's may depend on some quantity that is computed from the values of the neighborhood sites. Additionally, some CA's may use more than two or more sets of rules (which need not depend on the same set of neighborhood sites) that are applied consecutively. Many examples of update rule sets are given in [Gaylord and Nishidate, 1996]

Applying update rules to the CA lattice

The following functions can be used to apply update rules to all of the sites of a CA lattice.

For a rule that depends on the values of sites in the von Neumann neighborhood of a site, we define

VonNeumann[func__, lat_] :=   
        MapThread[func, Map[RotateRight[lat, #]&, {{0, 0}, {1, 0}, {0, -1}, {-1, 0}, {0, 1}}], 2]

(*Joe, Richard's formatting looked better than this in Mma3.0. So maybe you can align these a little better. Also in code below. David *)

and apply this function to the CA lattice using the anonymous function

VonNeumann[update, #]&

For a rule that depends on the values of sites in the Moore neighborhood of a site, we define

Moore[func__, lat_] :=   
  MapThread[func, 
                     Map[RotateRight[lat, #]&, 
                             {{0, 0}, {1, 0}, {0, -1}, {-1, 0}, {0, 1}, {1,  -1}, {-1, -1}, {-1, 1}, {1, 1}}], 2]

and apply this function to the CA lattice using the anonymous function

Moore[update, #]&

For a rule that depends on the values of sites in the MvonN neighborhood of a site, we define

MvonN[func__, lat_] :=   
     MapThread[func, Map[RotateRight[lat, #]&, {{0, 0}, {1, 0}, {0, -1}, {-1, 0}, {0, 1}, {1,  -1}, {-1, -1}, 
                                                                               {-1, 1}, {1, 1},  {2, 0}, {0, -2}, {-2, 0}, {0, 2}}], 2]

and apply this function to the CA lattice using the anonymous function

MvonN[update, #]&

The effect of applying these functions to a CA lattice can be seen by using a simple 3 by 3 matrix.

(mat3 = Partition[Range[9], 3]) // MatrixForm

Applying an undefined function, g, to mat3 using the VonNeumann function gives

VonNeumann[g, #]&[mat3]

The result is a 3-by-3 matrix in which each element is a function call  to five arguments. To identify the arguments in these function calls, we can focus on the center site in the following 5-by-5 matrix:

(mat5 = Partition[Range[25], 5] /.{13 -> site, 8 ->N, 14 -> E, 18 -> S, 12 -> W, 9 -> NE,  19 -> SE, 17 -> SW, 7 -> NW, 3 -> Nn, 15 -> Ee, 23 ->Ss, 11 -> Ww}) // MatrixForm

Applying each of the three anonymous functions defined above and the undefined function, g, to mat5 gives the following values for the center site in the resulting matrix:

VonNeumann[g, #]&[mat5][[3, 3]]

Moore[g, #]& [mat5][[3, 3]]

MvonN[g, #]&[mat5][[3, 3]]

CA Evolution

The CA can be made to evolve over a specified number of time steps by repeatedly applying one of the anonymous functions given above to the CA lattice. This is done using the general form

NestList[VonNeumann[update, #]&, initConfig, t],

where t represents the specified number of time steps. If  we expect that the CA may, at some point, stop evolving permanently (as opposed to a temporary pause), we can use the FixedPoint function with three argument values to halt the evolution of the CA when the system becomes invariant. This is done using the general form

FixedPointList[Moore[update, #]&, initConfig, t]

If we wish the CA program to stop running when either a certain condition is met, the system becomes static, or a specified number of time steps have taken place, we can use the FixedPoint operation with four arguments, the last of which is the SameTest option. This is done using the general form

FixedPointList[MvonN[update, #]&,  initconfig, t, SameTest -> ...]

Note: If we are only interested in the final configuration of the CA lattice, and not the intermediate configurations, we can use the Nest and FixedPoint functions.

Watching the CA evolve

We can create a flip-card animation of the evolution of a CA using the following commands:

ShowCA[list_, opts___]:=
 Map[Show[Graphics[RasterArray[list[[#]] /. 
                         {val1 -> RGBColor[x1, y1, z1], 
                          val2 -> RGBColor[x2, y2, z2], ...}] ], 
                   opts]&,
         Range[Length[list]] ]

Example: The Game of Life CA

The Game of Life CA is the most well-known of all cellular automata [Gardner 1983]. It is played on an n-by-n square lattice with periodic boundaries. Each lattice site has a value of either 0 or 1 and a lattice site whose value is 1 is said to be alive and a lattice site whose value is 0 is said to be dead. The initial lattice configuration is usually a random distribution of dead and living sites (although the lattice may be seeded with specific life-forms [see Gaylord and Nishidate, 1996]). The values of the sites of the lattice are updated for a fixed number of time steps or until the lattice configuration stops changing, using update rules that are based on values of a site and the nearest neighbor sites in its Moore neighborhood.

A general template for the Game of Life CA program, lacking only the specification of the update rules, can be written as:

OblaDeOblaDa[n_, t_]:=  
  Module[{initConfig, Moore, update},
	
	         initConfig = Table[Random[Integer],{n},{n}] ;
	
          update...
          update...
          .
          .
          .
         update...
               
         Moore[func__, lat_] :=   
                 MapThread[func, Map[RotateRight[lat, #]&, 
                                    {{0, 0}, {1, 0}, {0, -1}, {-1, 0}, {0, 1}, {1,  -1}, {-1, -1}, {-1, 1}, {1, 1}}],  2];
          
          FixedPointList[Moore[update, #]&, initConfig, t] ]

The Game of Life update rules

The life and death update rules for updating a lattice site in the Game of Life CA are as follows:

- any site (dead or alive) with three living nearest neighbor sites stays alive or is born.
- a living site with two living nearest neighbor sites remains alive.  
- all other sites either remain dead or die.

Traditionally, the implementation of these update rules involves determining the number of occurrences of 1's in the Moore neighborhood of a site [Gaylord and Nishidate 1996].  However, it is much more efficient, in terms of the runtime speed of the  program, to write life and death update rules specific to each of the 29 possible Moore neighborhood configurations (eg., update[1, 0, 1, 0, 1, 1, 0, 1, 1] = 0). Since typing in 512 update rules requires too much work, we will have Mathematica generate a look-up table of these rules for us.

Note: The code used in the next section to generate the 512 rewrite rules of the Game of Life, is a bit advanced (and demonstrates some of Mathematica's power). In most CA's [Gaylord and Nishidate 1996], the rule sets are considerably smaller and can be entered manually by the programmer rather than generated by the program itself.

Creating the Game of Life look-up rule table

The set of all possible Moore neighborhood configurations is computed using

LiveConfigs = Join[Map[Join[{0}, #]&, Permutations[{1, 1, 1, 0, 0, 0, 0, 0}]], 
                                Map[Join[{1}, #]&, Permutations[{1, 1, 1, 0, 0, 0, 0, 0}]], 
                               Map[Join[{1}, #]&, Permutations[{1, 1, 0, 0, 0, 0, 0, 0}]] ];

DieConfigs = Complement[Flatten[Map[Permutations, 
             Map[Join[Table[1, {#}], Table[0, {(9 - #)}] ]&, Range[0, 9]]], 1], LiveConfigs];

LiveConfigs and DieConfigs are nested lists. LiveConfigs contains lists of those Moore neighborhood configurations whose center site value will be updated to a value of 1 and DieConfigs contains lists of those Moore neighborhood configurations whose center site value will be updated to a value of 0 (in all of these lists, the center site value is given by the first element in the list). We can calculate the number of life-giving configurations and death-taking configurations, respectively.

Map[Length, {LiveConfigs, DieConfigs}]

We now create the update rules for each of the LiveConfigs and DieConfigs configurations as follows:

Apply[(update[##] = 1)&, LiveConfigs, 1];

Apply[(update[##] = 0)&, DieConfigs, 1];

We show how this works below for a set of nine DieConfigs configuration.

Apply[(update[##] = 0)&, Permutations[{1, 1, 1, 1, 1, 1, 1, 1, 0}], 1];

?update

It would take too much space to print  (and be too boring to read) all 512 update rules, but we can check that we have, in fact, generated all of the update rules.

Length[DownValues[update]]

Overall then, the following set of commands creates the entire set of 512 update rules for the Game of Life CA:

(*Make WRICourier, or what we are using. I'll write you about this. David *)

LiveConfigs = Join[Map[Join[{0}, #]&, Permutations[{1, 1, 1, 0, 0, 0, 0, 0}]], 
                                Map[Join[{1}, #]&, Permutations[{1, 1, 1, 0, 0, 0, 0, 0}]], 
                                Map[Join[{1}, #]&, Permutations[{1, 1, 0, 0, 0, 0, 0, 0}]] ];
DieConfigs = Complement[Flatten[Map[Permutations, 
                                     Map[Join[Table[1, {#}], Table[0, {(9 - #)}] ]&, Range[0, 9]]], 1], LiveConfigs]; 
Apply[(update[##] = 1)&, LiveConfigs, 1];  
Apply[(update[##] = 0)&, DieConfigs, 1];

To create the Game of Life CA, all we have to do is to insert these commands into the program template.

The Program

OblaDeOblaDa[n_, t_]:=  
  Module[{initConfig, Moore, update, LiveConfigs, DieConfigs},
	
   	initConfig = Table[Random[Integer],{n},{n}] ;
	
    LiveConfigs = Join[Map[Join[{0}, #]&, Permutations[{1, 1, 1, 0, 0, 0, 0, 0}]], 
                                    Map[Join[{1}, #]&, Permutations[{1, 1, 1, 0, 0, 0, 0, 0}]], 
                                    Map[Join[{1}, #]&, Permutations[{1, 1, 0, 0, 0, 0, 0, 0}]] ];
    DieConfigs = Complement[Flatten[Map[Permutations, 
                 Map[Join[Table[1, {#}], Table[0, {(9 - #)}] ]&, Range[0, 9]]], 1], LiveConfigs]; 
    Apply[(update[##] = 1)&, LiveConfigs, 1];  
    Apply[(update[##] = 0)&, DieConfigs, 1]; 
               
    Moore[func__, lat_] :=   
        MapThread[func, Map[RotateRight[lat, #]&, 
                                   {{0, 0}, {1, 0}, {0, -1}, {-1, 0}, {0, 1}, {1,  -1}, {-1, -1}, {-1, 1}, {1, 1}}], 2];

    FixedPointList[Moore[update, #]&, initConfig, t] ]

The following graphics show the state of a Life game at times t = 1, t = 64, and t = 128, played on a 32-by-32 grid.

Note: This material is adopted from our new book "Modeling Nature: Cellular Automata Simulations using Mathematica" [Gaylord and Nishidate, 1996].

References

Gardner, M. "The fantastic combinations of John Conway's new solitary game 'Life.'" Scientific American, 223, 120-123.

Gaylord, Richard J and Nishidate, Kazume. "Modeling Nature: Cellular Automata Simulations with Mathematica". TELOS/Springer-Verlag (1996).

ABOUT THE EDITORS. (*Standard for this column.*)

ELECTRONIC SUBSCRIPTIONS

Included in the distribution for each electronic subscription is the file ToolKit.nb containing Mathematica code for the material described in this article.