# 11.3: Simulating Cellular Automata

- Page ID
- 7831

\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)

\( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)

\( \newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\)

( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\)

\( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)

\( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\)

\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)

\( \newcommand{\Span}{\mathrm{span}}\)

\( \newcommand{\id}{\mathrm{id}}\)

\( \newcommand{\Span}{\mathrm{span}}\)

\( \newcommand{\kernel}{\mathrm{null}\,}\)

\( \newcommand{\range}{\mathrm{range}\,}\)

\( \newcommand{\RealPart}{\mathrm{Re}}\)

\( \newcommand{\ImaginaryPart}{\mathrm{Im}}\)

\( \newcommand{\Argument}{\mathrm{Arg}}\)

\( \newcommand{\norm}[1]{\| #1 \|}\)

\( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\)

\( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\AA}{\unicode[.8,0]{x212B}}\)

\( \newcommand{\vectorA}[1]{\vec{#1}} % arrow\)

\( \newcommand{\vectorAt}[1]{\vec{\text{#1}}} % arrow\)

\( \newcommand{\vectorB}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)

\( \newcommand{\vectorC}[1]{\textbf{#1}} \)

\( \newcommand{\vectorD}[1]{\overrightarrow{#1}} \)

\( \newcommand{\vectorDt}[1]{\overrightarrow{\text{#1}}} \)

\( \newcommand{\vectE}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash{\mathbf {#1}}}} \)

\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \)

\( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)

\(\newcommand{\avec}{\mathbf a}\) \(\newcommand{\bvec}{\mathbf b}\) \(\newcommand{\cvec}{\mathbf c}\) \(\newcommand{\dvec}{\mathbf d}\) \(\newcommand{\dtil}{\widetilde{\mathbf d}}\) \(\newcommand{\evec}{\mathbf e}\) \(\newcommand{\fvec}{\mathbf f}\) \(\newcommand{\nvec}{\mathbf n}\) \(\newcommand{\pvec}{\mathbf p}\) \(\newcommand{\qvec}{\mathbf q}\) \(\newcommand{\svec}{\mathbf s}\) \(\newcommand{\tvec}{\mathbf t}\) \(\newcommand{\uvec}{\mathbf u}\) \(\newcommand{\vvec}{\mathbf v}\) \(\newcommand{\wvec}{\mathbf w}\) \(\newcommand{\xvec}{\mathbf x}\) \(\newcommand{\yvec}{\mathbf y}\) \(\newcommand{\zvec}{\mathbf z}\) \(\newcommand{\rvec}{\mathbf r}\) \(\newcommand{\mvec}{\mathbf m}\) \(\newcommand{\zerovec}{\mathbf 0}\) \(\newcommand{\onevec}{\mathbf 1}\) \(\newcommand{\real}{\mathbb R}\) \(\newcommand{\twovec}[2]{\left[\begin{array}{r}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\ctwovec}[2]{\left[\begin{array}{c}#1 \\ #2 \end{array}\right]}\) \(\newcommand{\threevec}[3]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\cthreevec}[3]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \end{array}\right]}\) \(\newcommand{\fourvec}[4]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\cfourvec}[4]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \end{array}\right]}\) \(\newcommand{\fivevec}[5]{\left[\begin{array}{r}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\cfivevec}[5]{\left[\begin{array}{c}#1 \\ #2 \\ #3 \\ #4 \\ #5 \\ \end{array}\right]}\) \(\newcommand{\mattwo}[4]{\left[\begin{array}{rr}#1 \amp #2 \\ #3 \amp #4 \\ \end{array}\right]}\) \(\newcommand{\laspan}[1]{\text{Span}\{#1\}}\) \(\newcommand{\bcal}{\cal B}\) \(\newcommand{\ccal}{\cal C}\) \(\newcommand{\scal}{\cal S}\) \(\newcommand{\wcal}{\cal W}\) \(\newcommand{\ecal}{\cal E}\) \(\newcommand{\coords}[2]{\left\{#1\right\}_{#2}}\) \(\newcommand{\gray}[1]{\color{gray}{#1}}\) \(\newcommand{\lgray}[1]{\color{lightgray}{#1}}\) \(\newcommand{\rank}{\operatorname{rank}}\) \(\newcommand{\row}{\text{Row}}\) \(\newcommand{\col}{\text{Col}}\) \(\renewcommand{\row}{\text{Row}}\) \(\newcommand{\nul}{\text{Nul}}\) \(\newcommand{\var}{\text{Var}}\) \(\newcommand{\corr}{\text{corr}}\) \(\newcommand{\len}[1]{\left|#1\right|}\) \(\newcommand{\bbar}{\overline{\bvec}}\) \(\newcommand{\bhat}{\widehat{\bvec}}\) \(\newcommand{\bperp}{\bvec^\perp}\) \(\newcommand{\xhat}{\widehat{\xvec}}\) \(\newcommand{\vhat}{\widehat{\vvec}}\) \(\newcommand{\uhat}{\widehat{\uvec}}\) \(\newcommand{\what}{\widehat{\wvec}}\) \(\newcommand{\Sighat}{\widehat{\Sigma}}\) \(\newcommand{\lt}{<}\) \(\newcommand{\gt}{>}\) \(\newcommand{\amp}{&}\) \(\definecolor{fillinmathshade}{gray}{0.9}\)Despite their capability to represent various complex nonlinear phenomena, CA are relatively easy to implement and simulate because of their discreteness and homogeneity.

There are existing software tools^{2} and online interactive demonstrations^{3} already available for cellular automata simulation, but it is nonetheless helpful to learn how to develop a CA simulator by yourself. Let’s do so in Python, by working through the following example step by step. The CA model we plan to implement here is a binary CA model with the *droplet rule *[4]. Its state-transition function can be understood as a model of panic propagation among individuals sitting in a gym after a fire alarm goes off. Here is the rule (which uses the Moore neighborhoods):

- A normal individual will get panicky if he or she is surrounded by four or more panicky individuals.
- A panicky individual will remain panicky if he or she is surrounded by three or more panicky individuals. Otherwise he or she will revert back to normal.

Note that this rule can be further simplified to the following single rule:

- If there are four or more panicky individuals within the neighborhood, the central cell will become panicky; otherwise it will become normal.

Here are other model assumptions:

- Space: 2-D, \(n×n\) (\(n = 100\) for the time being)
- Boundary condition: periodic
- Initial condition: Random (panicky individuals with probability \(p\); normal ones with probability \(1−p\))

We will use pycxsimulator.py for dynamic CA simulations. Just as before, we need to design and implement the three essential components—initialization, observation, and updating.

To implement the initialization part, we have to decide how we are going to represent the states of the system. Since the configuration of this CA model is a regular two-dimensional grid, it is natural to use Python’s `array`

data structure. We can initialize it with randomly assigned states with probability \(p\). Moreover, we should actually prepare two such arrays, one for the current time step and the other for the next time step, in order to avoid any unwanted conﬂict during the state updating process. So here is an example of how to implement the initialization part:

Here,` zeros([a, b])`

is a function that generates an all-zero array with `a`

rows and `b`

columns. While `config `

is initialized with randomly assigned states (1 for panicky individuals, 0 for normal), `nextconfig `

is not, because it will be populated with the next states at the time of state updating.

The next thing is the observation. Fortunately, pylab has a built-in function` imshow`

that is perfect for our purpose to visualize the content of an array:

Note that the` cmap `

option in` imshow`

is to specify the color scheme used in the plot. Without it, pylab uses dark blue and red for binary values by default, which are rather hard to see.

The updating part requires some algorithmic thinking. The state-transition function needs to count the number of panicky individuals within a neighborhood. How can we do this? The positions of the neighbor cells within the Moore neighborhood around a position \((x,y)\) can be written as

\[{(x',y')|x-1 \leq x'\leq x+1, y-1 \leq y' \leq y+1}.\label{(11.3)} \]

This suggests that the neighbor cells can be swept through using nested `for`

loops for relative coordinate variables, say, \(dx\) and \(dy\), each ranging from `−1`

to `+1`

. Counting panicky individuals (1’s) using this idea can be implemented in Python as follows:

Here,` dx`

and` dy`

are relative coordinates around \((x,y)\), each varying from`−1`

to` +1`

. They are added to `x`

and` y`

, respectively, to look up the current state of the cell located at \((x + dx,y + dy)\) in `config`

. The expression (...) % n means that the value inside the parentheses is contained inside the `[0,n − 1]`

range by the mod operator (%). This is a useful coding technique to implement periodic boundary conditions in a very simple manner.

The counting code given above needs to be applied to all the cells in the space, so it should be included in another set of nested loops for `x `

and` y`

to sweep over the entire space. For each spatial location, the counting will be executed, and the next state of the cell will be determined based on the result. Here is the completed code for the updating:

Note the swapping of `config`

and `nextconfig`

at the end. This is precisely the same technique that we did for the simulation of multi-variable dynamical systems in the earlier chapters.

By putting all the above codes together, the completed simulator code in its entirety looks like this:

When you run this code, you see the results like those shown in Fig. 11.3.2. As you can see, with the initial configuration with \(p = 0.1\), most of the panicky states disappear quickly, leaving only a few self-sustaining clusters of panicky people. They may look like condensed water droplets, hence the name of the model (droplet rule).

Modify Code 11.5 to implement a simulator of the Game of Life CA. Simulate the dynamics from a random initial configuration. Measure the density of state 1’s in the configuration at each time step, and plot how the density changes over time. This can be done by creating an empty list in the` initialize `

function, and then making the measurement and appending the result to the list in the `observe`

function. The results stored in the list can be plotted manually after the simulation, or they could be plotted next to the visualization using pylab’s subplot function during the simulation.

Modify Code 11.5 to implement a simulator of the majority rule CA in a two-dimensional space. Then answer the following questions by conducting simulations:

• What happens if you change the ratio of binary states in the initial condition?

• What happens if you increase the radius of the neighborhoods?

• What happens if you increase the number of states?

The second model revision in the previous exercise (increasing the radius of neighborhoods) for the majority rule CA produces quite interesting spatial dynamics. Specifically, the boundaries between the two states tend to straighten out, and the characteristic scale of spatial features continuously becomes larger and larger over time (Fig. 11.3.3). This behavior is called *coarsening* (to be more specific, non-conserved coarsening). It can be seen in many real-world contexts, such as geographical distributions of distinct cultural/political/linguistic states in human populations, or incompatible genetic types in animal/plant populations.

What is most interesting about this coarsening behavior is that, once clear boundaries are established between domains of different states, the system’s macroscopic behavior can be described using emergent properties of those boundaries, such as their surface tensions and direction of movement[41,42]. The final fate of the system is determined not by the relative frequencies of the two states in the space, but by the topological features of the boundaries. For example, if a big “island” of white states is surrounded by a thin “ring” of black states, the latter minority will eventually dominate, no matter how small its initial proportion is (even though this is still driven by the majority rule!). This is an illustrative example that shows how important it is to consider emergent macroscopic properties of complex systems, and how counter-intuitive their behaviors can be sometimes.

Here is one more interesting fact about CA dynamics. The droplet/panic model discussed above has an interesting property: When you increase the initial density of panicky people (e.g., \(p = 0.3\)), the result changes dramatically. As seen in Fig. 11.3.4, the initially formed clusters tend to attach to each other, which makes their growth unstoppable. The whole space will eventually be filled up with all panicky people, which could be a disaster if this was a real situation.

You can explore the value of \(p\) to find out that the transition between these two distinct behaviors takes place at a rather narrow range of \(p\). This is an example of a* phase transition*, which is defined informally as follows:

A* phase transition* is a transition of macroscopic properties of a collective system that occurs when its environmental or internal conditions are varied

A familiar example of phase transitions is the transition between different phases of matter, i.e., solid, liquid, and gas, which occur when temperature and/or pressure are varied. Phase transitions can be characterized by measuring what physicists call order *parameters*^{4 }that represent how ordered a macroscopic state of the system is. A phase transition can be understood as a bifurcation observed in the macroscopic properties (i.e., order parameters) of a collective system.

Implement an interactive parameter setter for \(p\) in Code 11.5. Then conduct systematic simulations with varying \(p\), and identify its critical value below which isolated clusters are formed but above which the whole space is filled with panic.

__ __

^{2}Most notable is Golly (http://golly.sourceforge.net/).

^{3}For example, check out Wolfram Demonstrations Project(demonstrations.wolfram.com/)and Shodor.org’s interactive activities (http://www.shodor.org/interactivate/activities/).

^{4}Note that the word“parameter” in this context means an outcome of a measurement, and not a condition or input as in “model parameters.”