Exploring Computational Irreducibility 
and the Predictability of Complex Systems Using Mathematica

Seth J. Chandler
University of Houston Law Center

SChandler@uh.edu

Introduction

The principle of computational irreducibility set forth by Stephen Wolfram in A New Kind of Science  [Wolfram 2003]  suggests that the consequences of very simple programs are often impossible to foretell without actually running the program. Thus, to the extent, as Wolfram and others claim, that physical, biological or social phenomena may be emulated by simple programs, the principle of computational irreducibility asserts that their long run behavior is simply not susceptible to intelligent prediction. A butterfly flapping its wings in Brazil may result deterministically and yet unpredictably in a Texas tornado two weeks hence. A recent article by Navot Israeli and Nigel Goldenfeld, Computational Irreducibility and the Predictability of Complex Physical Systems, publised in the February 20, 2004 volume of Physical Review Letters [Israeli & Goldenfeld 2004], scrutinizes this proposition. It shows that, if certain "microscopic" details are disregarded through a process known as "coarse graining", the long run behavior of many elementary cellular automata can be emulated by other cellular automaton that are often less complex than that of the original. And this is true even with respect to some (but apparently not all) "Class 3" cellular automata that exhibit complex behavior. Their paper is thus important in several regards. First, from a theoretical perspective, the paper suggests that, although some automata may be computationally irreducible, others hitherto placed in that category may have their behavior approximated (and thus "reduced" in some sense) using other cellular automata. Accordingly, the paper speculates as to whether Wolfram's classifications of automata and the understanding of computational irreducibility he presented in A New Kind of Science may need some refinement. Second, from a practical perspective, to the extent that cellular automata can be used to model physical, biological or social phenomena, systems hitherto thought to be completely unpredictable in the long run (such as the weather), might now be subject to some macroscopic predictability. While we might still not have an ability to forecast the weather two weeks hence in a particular location, we might be able to say that at least it will not be raining in both that location and a neighboring town.

This article exposits in some detail the methodology developed by Israeli and Goldenfeld and shows how Mathematica may be used to explore and extend their work. The main purpose of this effort is to make their work more accessible and to permit a broader audience to conduct at least many of the relevant experiments. This "translational" scholarship is important because cellular automata can help explicate the behavior not only of physical systems, but biological and social systems as well. And perhaps because Mathematica is such a natural language for the study of cellular automata, the code for accomplishing the translation proves relatively simple. I proceed by explaining how Mathematica can be used to "coarse grain" a system. I then discuss an algorithm for determining the identity of any cellular automaton that emulates the behavior of such a coarse grained system. I then use this algorithm to replicate all but the most complex experiments conducted by Israeli and Goldenfeld. These experiments show the reduction in complexity that can be obtained if microscopic details of the automaton's behavior are disregarded. They also show, however, that the reduction in complexity may not be terribly large. Further analysis of the experiments using Mathematica exposes what appear to be previously unreported symmetries between characteristics of the elementary cellular automaton under study, the method of coarse graining employed, and the emulating elementary cellular automaton. These symmetries suggest, though they do not prove, that there may be somewhat more accelerated methods for finding emulating cellular automata. I conclude the paper by extending Israeli and Goldenfeld's methodology to the domain of 3-color totalistic automata. I develop sampling methods that accelerate discovery of successful emulations. I find that complex 3-color totalistic automata can, like their 2-color elementary cousins, frequently be emulated by simpler systems, although, again, the enhancement in the ability to predict the future evolution of the system may be small.

Speed issues associated with Mathematica make it difficult as a practical matter to replicate the most complex experiments conducted by Israeli and Goldenfeld. To some extent, these limitations may be fought back through sampling methods developed in this article. They may also be softened to some extent by faster hardware, grid computing, patience, or perhaps algorithms that exploit Mathematica more effectively than has been possible in this initial foray. The analysis here suggests, however, that, absent the development of truly clever algorithms, there may be more fundamental problems of computational explosion and irreducibility with trying to find emulating cellular automata that rely on larger block size. Moreover, the reduction in complexity obtained by the sort of emulation described here is often not terribly substantial. Many cellular automata may retain their secrets.

Finding Emulating Cellular Automata

Accelerating The Computation of a Cellular Automaton

Ordinarily one computes the evolution of a cellular automaton one step at a time. One applies a set of rules to an input, produces an output, and iteratively applies the same rules over and over again. Thus, to determine the outcome of a cellular automaton after t iterations often requires t invocations of the cellular automaton rule transformation set. One can, however, accelerate the computation of this evolutionary path by putting the computation into a kind of "warp drive." The idea is simple. To compute the value of an 𝒩-length sequence of an infinite length cellular automaton whose rule operates over radius r, one must know, at t steps in the past, the value of the 𝒩 + 2 r t  sites that are centered "above" the original sequence.  This is so because, with each step backwards, one must know (1) the values of the preceding states of the r sites to the left of the left boundary of the 𝒩-length sequence, (2) the values of the preceding states of the 𝒩-length sequence, and (3) the values of the preceding state for the r sites to the right of the right boundary of the 𝒩-length sequence. Inversely, therefore, if σ is the length of some (non-periodic) sequence of an infinite length cellular automaton whose rule operates over radius r, and one knows the values of σ, then, after t iterations, one can also compute the values of the σ - 2 r t length sequence centered below the original sequence. Thus, by some simple algebra, if one wants to know the values of an 𝒩-length sequence derived from an q length sequence, one must run the automaton \(TraditionalForm\`Floor[\(σ - 𝒩\)\/\(2\ r\)]\) times.

This computation can be done in Mathematica using DeBruijn Sequences. A Debruijn Sequence is the shortest circular sequence of length (σ^a) such that every string of length n on the alphabet a of size σ occurs as a contiguous subrange of the sequence described by a. [Pemmaraju 2003] For example, the DeBruijn Sequence of length 3 on the alphabet {x,y,z} is {z, x, z, y, x, y, x, x, z, x, x, x, y, z, x, y, y, 
        y, z, z, z, y, z, y, y, x, z}. (Mathematica's Combinatorica package contains a routine for swiftly finding DeBruijn Sequences: DeBruijnSequence[
        a, (2  r + 1) 𝒩]. Thus, if one runs a cellular automaton α just once (over 𝒩 iterations) on initial conditions that are the DeBruijn Sequence of length σ on the alphabet {0,…,k-1} one has all the information needed to compute the accelerated cellular automaton. To save time, I load in some pre-computed DeBruijn Sequences: binary alphabets of lengths 6, 9 and 12 and a ternary alphabet on length 6. The DeBruijn package also loads in some graphics and housekeeping code.

Needs["DiscreteMath`DeBruijn`"]

DeBruijnSequence[{x, y, z}, 3]

The βrules function implementing this "warp drive" computation takes as input (1) virtually any specification for a one-dimensional cellular automaton ("the source cellular automaton"), (2) the DeBruijn Sequence (a) over an alphabet a with the same number of  symbols as there are colors in the source cellular automaton and (b) of a length σ equal to the product of (i)  the block size 𝒩 and (ii)  the diameter 2r+1 of the source cellular automaton , (3) σ itself, and (4) the block size 𝒩. It produces a list each element of which is a list containing the σ-length input and the 𝒩-length output.

βrules[α : {_, _, r_}, 
        debruijnsequence_List, σ_Integer, 𝒩_Integer] := 
      Module[{σ2, 𝒩2, t, cadb}, 
        t = Floor[(σ - 𝒩)/2  
              r]; 
σ2 = 
          Ceiling[σ/2]; 
𝒩2 = 
          Ceiling[𝒩/2]; 
cadb = 
          CellularAutomaton[α, debruijnsequence, 
            t, {{0, t, t}}]; 
MapThread[Join, 
          List[Partition[cadb[[1]], σ, 1, {σ2, σ2}], 
            Partition[cadb[[-1]], 𝒩, 
              1, {𝒩2, 𝒩2}]]]]

To illustrate, here is a short version of the transformation rules for elementary rule 105 with a block size of 2.  The function output2rules presents the output in a familiar rule-based format.

output2rules[x_List, 𝒩_] := 
      Map[Rule[Drop[#, -𝒩], 
            Take[#, -𝒩]] &, x]

Short[
      output2rules[β105 = βrules[{105, 2, 1}, db26, 6, 2], 2], 
      2]

Here is another example showing the transformation rules for elementary rule 146 with a block size of 3.

Short[{output2rules[β146 = βrules[{146, 2, 1}, db29, 9, 3], 
          3]}, 2]

The code works for automata with more than 2 colors and for arbitrarily weighted automata.

Short[{output2rules[β4625597484987 = βrules[{4625597484987, 
                3, 1}, db36, 6, 2], 2], 
        output2rules[βweighted = βrules[{101, {2, {4, 2, 3}}, 1}, 
              db26, 6, 2], 2]}, 2]

Coarse Graining is Supercelling followed by Projection

Basic Code

I can now take our list of rules and see if they describe a cellular automaton when the input and output are "coarse grained." By coarse grained, Israeli and Goldenfeld refer to a two step process of "supercelling" and then "projecting."  "Supercelling" means to take the base k digit sequence of the input ("the raw input") and output ("the raw output") for each rule and convert each of them into their equivalent base (k^𝒩) digit sequences. Thus, for a block size of 2 (𝒩=2), the rule {0, 1, 0, 1, 1, 1} -> {1, 1} would be supercelled into the base 4 digit rule {1, 1, 3} -> {3}. So, if the original input has σ digits, the supercelled version on the input has (σ/𝒩)digits and the supercelled version of the output has 1 digit. Here is an implementation of this concept in Mathematica. The supercell function takes as its first input argument the output from the βrules function. The remaining arguments simply tell supercell the size of the block size 𝒩 and k, the base in which the digits are to be understood.

frombase2base[x_List, from_Integer, to_Integer, len_] := 
      IntegerDigits[FromDigits[x, from], to, len]

supercell[x : List[__Integer], 𝒩_Integer, k_Integer] := 
        With[{kn = (k^𝒩), σ = 
              Length[x] - 𝒩}, 
          Join[frombase2base[Drop[x, -𝒩], k, 
              kn, σ/𝒩], 
            frombase2base[Take[x, -𝒩], k, kn, 
              1]]];

I then have the option of projecting each digit in the supercell onto a target alphabet of size (Subscript[k, β]). (Generally, though not necessarily, I want (Subscript[k, β]) to be the same as k). Notice that, because there can be (k^𝒩) different digits in the supercells, the projection operation needs to know how to handle (k^𝒩) different inputs. I implement projection using Mathematica's ReplaceAll command.

projection[proj : List[__Rule] | proj : List[__List] | proj_Rule, 
        x_List] := ReplaceAll[x, proj]

I can combine the supercelling and projecting operations into a single coarsegrain function.

coarsegrain[proj : List[__Rule] | proj : List[__List] | proj_Rule, 
        x : List[___Integer], 𝒩_Integer, k_Integer] := 
      projection[proj, supercell[x, 𝒩, k]]

Representation of Projections

Because a projection maps (k^𝒩) possible inputs onto (Subscript[k, β]) possible outputs, there are  ((Subscript[k, β])^((k^𝒩))) possible projections; I can accordingly represent the projection operation by a number between 0 and  ((Subscript[k, β])^((k^𝒩))) - 1. I thus number projections in a fashion akin to the way that Wolfram numbers cellular automata.  The proj2rule function translates between projection number and rule for a given k, (Subscript[k, β]) and 𝒩. It is also capable of producing the universe of possible projections for a given k, (Subscript[k, β]) and 𝒩.

SetAttributes[proj2rule, Listable];
proj2rule[proj_Integer, 𝒩_Integer, k_Integer, 
          kβ_Integer] := 
        Thread[Rule[Range[(k^𝒩) - 1, 0, -1], 
            IntegerDigits[proj, kβ, 
              (k^𝒩)]]];
proj2rule[𝒩_Integer, k_Integer, kβ_Integer] := 
        proj2rule[
          Range[0, 
            k(β^((k^𝒩))) - 1], 𝒩, k, 
          kβ];

Thus, projection 73 for k = 2, 
      (Subscript[k, β]) = 2 and 𝒩 = 3, projection 54321 for k = 2, kβ = 2, 𝒩 = 4, and projection 12345678 for  k = 2, kβ = 3, 𝒩 = 4 are respectively computed as follows:

{proj2rule[73, 2, 2, 3], proj2rule[54321, 2, 2, 4], 
      proj2rule[123456789, 2, 3, 4]}

I can also invert the process.

rule2proj[x : List[__Rule], k_Integer, kβ_Integer] := 
      With[{𝒩 = Log[k, Length[x]]}, {FromDigits[Last /@ x, 
            kβ], k, kβ, 𝒩}]

{rule2proj[proj2rule[73, 3, 2, 2], 2, 2], 
      rule2proj[proj2rule[54321, 4, 2, 2], 2, 2], 
      rule2proj[proj2rule[12345678, 4, 2, 3], 2, 3]}

A Timed Example Using Elementary Rule 146

I can now show the process step by step using elementary rule 146 with a block size of 3. I use a Map operation to show how supercelling and projecting works on all possible inputs.

Timing[With[{p73 = proj2rule[73, 3, 2, 2]}, 
          Map[coarsegrain[p73, #, 3, 2] &, βrules[{146, 2, 1}, db29, 9, 
              3]]];]

I can also do this without loss of much efficiency by using a FoldList operation in which the first argument to FoldList proves irrelevant. Recall that FoldList[f, 0, {a, b, c}] is {0, f[0, a], f[f[0, a], b], 
        f[f[f[0, a], b], c]}. (Mathematica's Map command may generally be thought of as a "subset" of FoldList).

Timing[
      With[{p73 = proj2rule[73, 3, 2, 2]}, 
          Rest[FoldList[
              coarsegrain[p73, #2, 3, 
                  2] &, "irrelevant", βrules[{146, 2, 1}, db29, 9, 
                3]]]] === With[{p73 = proj2rule[73, 3, 2, 2]}, 
          Map[coarsegrain[p73, #, 3, 2] &, βrules[{146, 2, 1}, db29, 9, 
              3]]]]

Determining If a Cellular Automaton Can Emulate A Coarse Grained System

With this background, I can now determine if the coarse grained system amounts to a cellular automaton. This is important because, if the emulating system is an automaton that is "simpler" than the automaton that generated the coarse grained system, I can say that, if one is willing to surrender information about certain details, the system behaves in a simpler fashion than one might have originally supposed. Macro level predictions about computationally irreducible systems may thus prove possible.

The algorithm for determining whether one of these systems is a cellular automaton is relatively straightforward. I create a throwvalue, something that both represents a "clearly wrong guess" as to the Wolfram number of any cellular automaton that our system might be emulating (the "hypothesized emulating cellular automaton" or "HECA") and that serves as an output if the system proves not to be a cellular automaton. Specifically, I set the throwvalue to a list digits that are each equal to -1.  The length of the list is the number of possible inputs that the HECA might encounter.  For each raw input that I encounter, I coarse grain it and see whether I already know how the HECA is supposed to handle it. If I don't know, I say that the output of the cellular automaton is equal to the output of the coarse graining. If I do know, and that already-known output of the cellular automaton is not equal to the output from the coarse graining, I abort the process. I abort because the HECA would then be multivalued. Otherwise, I basically do nothing. I implement this with compiled function modifyωc to enhance speed. The arguments are an array ω containing the best current guess as the digits of the HECA, cgc, a list of integers that joins the list of input integers from a warp drive computation and the list of output integers, z, the number of digits in the HECA, u, the possible values the HECA must know how to handle, kβ, the number of colors in the HECA, wts, the weights accorded to the input digits that the HECA confronts, and throwvalue, the previously discussed clearly wrong guess as to the HECA. The function outputs a new guess as to the digits of the HECA.

modifyωc = 
        Compile[{{ω, _Integer, 1}, {cgc, _Integer, 
              1}, {z, _Integer}, {u, _Integer, 
              1}, {kβ, _Integer}, {wts, _Integer, 
              1}, {throwvalue, _Integer, 
              1}}, 
Module[{index, 
              value, ωindex}, 
index = 
              z + 1 - Position[u, Most[cgc] . wts][[1, 
                    1]]; 
value = 
              Last[cgc]; 
ωindex = ω[[index]]; 
If[ωindex == 
                throwvalue[[index]], ReplacePart[ω, value, index], 
              If[ωindex ≠ value, 
                throwvalue, ω]]], {{index, _Integer}, {value, _Integer}, {ωindex, _Integer}, {_FromDigits, _Integer}}];

I can then put the process together as follows in the function findemulatingca. Notice that as in the code above, the Fold command proves the key. Basically I take our starting guess for the Wolfram number of the automaton and refine it as the Fold folds in additional pieces of information about the coarse graining results. If one of those pieces of information shows the system cannot be a cellular automaton, I abort. If at the end of the process, I don't know how the system would handle certain types of input, I call findemulatingca recursively, but this time assuming that the HECA has one fewer colors.

Options[
          findemulatingca] = {EmulatingWeights -> Elementary};

\(findemulatingca[
        i_Integer, β : List[__List], 𝒩_Integer, 
        k_Integer, kβ_Integer?\((# >= 2 &)\), 
        opts___] := 
Catch[
        With[{w = 
              Module[{emulwts, u, z, tv, 
                  p}, 
emulwts = \(EmulatingWeights /. {opts}\) /. Options[findemulatingca]; 
If[
                  emulwts === Elementary, 
                  emulwts = kβ\^{2, 1, 0}]; 
u = 
                  Union[Strings[Range[0, kβ - 1], Length[emulwts]] . 
                      emulwts]; 
z = 
                  Length[u]; 
tv = 
                  Table[\(-1\), {z}]; 
p = 
                  proj2rule[i, 𝒩, k, 
                    kβ]; 
Fold[
                  With[{m = 
                          modifyωc[#1, 
                            coarsegrain[p, #2, 𝒩, k], z, u, 
                            kβ, emulwts, tv]}, 
                      If[m === tv, Throw["multivalued"], m]] &, 
                  tv, β]]}, 
If[
            FreeQ[w, \(-1\)], {FromDigits[w, kβ], kβ, 1}, 
            findemulatingca[i, β, 𝒩, k, kβ - 1, 
              opts]]]]\)

findemulatingca[
        i_Integer, β : List[__List], 𝒩_Integer, 
        k_Integer, kβ_?(# < 2 &), 
        opts___] := "Not a Cellular Automaton"

This implementation allows us to find which, if any projections, create an emulating cellular automaton. The whole thing works quite swiftly. Learning about a particular projection takes very little time. Learning about all possible projections takes on the order of a second, depending on the speed of one's machine. I do this through the function f below.

cadomain[k_, r_] := 
      Which[Head[k] === Integer, 
Range[0, 
          (k^(2  r + 1)) - 1], 
Head[k] === 
          List, 
With[{kv = k[[1]], wts = k[[2]]}, 
          Union[Strings[Range[0, kv - 1], Length[wts]] . wts]]]

Options[f] = {EmulatingWeights -> Elementary, 
          Projections -> All};

f[α : {n_List, k_, r_}, db_List, 𝒩_, kβ_, 
        opts___] := 
Module[{kv, 
          proj}, 
Which[Head[k] === Integer, kv = k, 
          Head[k] === List, 
          kv = k[[1]]]; 
proj = Projections /. {opts} /. Options[f]; 
Which[proj === All, 
          proj = Range[0, 
              k(β^((kv^𝒩))) - 
                1]]; 
Map[
          Function[Α, {Α, 
              Cases[With[{β = βrules[{Α, k, r}, 
                        db, (2  r + 
                              1) 𝒩, 𝒩]}, 
                  Map[{#, 
                        findemulatingca[#, β, 𝒩, kv, 
                          kβ, opts]} &, 
                    proj]], {_Integer, {__Integer}}]}], n]]

f[α : {All, k_, r_}, db_List, 𝒩_, kβ_, 
        opts___] := 
Module[{kv}, 
Which[
          Head[k] === Integer, kv = k, Head[k] === List, kv = k[[1]]]; 
        With[{n = Range[0, kv^Length[cadomain[k, r]] - 1]}, 
          f[{n, k, r}, db, 𝒩, kβ, opts]]]

The following code shows how the f function operates, as well as providing a sense of computation time.

{Timing[
        Short[f[{{146, 147}, 2, 1}, db29, 3, 2, 
            Projections -> Range[0, 127]], 5]], 
      Timing[Short[f[{{105}, 2, 1}, db26, 2, 2]]], 
      Timing[Short[
          f[{{146, 147}, 2, 1}, db29, 3, 2, Projections -> Range[120, 130]], 
          5]]}

The results above show, for example, that the Class 3 elementary cellular automaton 146 can be emulated under projection 128 (proj2rule[128, 3, 2, 
          2] ={7 -> 1, 6 -> 0, 5 -> 0, 
          4 -> 0, 3 -> 0, 2 -> 0, 1 -> 0, 
          0 -> 0} ) by the Class 1 elementary cellular automaton 128, which takes all inputs and makes them zero. Thus, one knows that, in the long run, Rule 146 creates subsequences of 3 whose "supercell" values are 0 through 6, but not 7. We have thus gained some information about the future states of systems produced by Rule 146.

SetOptions[Graphics, ImageSize -> 12*72]; 
    SetOptions[GraphicsArray, ImageSize -> 18*72];

Show[
        GraphicsArray[{RasterGraphics[
              CellularAutomaton[{146, 2, 1}, db29, 300], 
              PlotLabel -> "Rule 146"], 
            RasterGraphics[CellularAutomaton[{128, 2, 1}, db29, 300], 
              PlotLabel -> "Rule 128"]}]];

Similarly, the results above show that the Class 3 elementary cellular automaton 105 can be emulated under projection 3 (and others) ({3 -> 0, 2 -> 0, 1 -> 1, 
        0 -> 1} by the Class 3 elementary cellular automaton 150. Thus one chaotic pattern can be emulated by another chaotic pattern in which supercells of 0 and 1 are replaced by 1s and supercells of 2 and 3 are replaced by zeros. We have not gained much ability to predict the long run behavior of Rule 105 systems (because Rule 150 is about as complex), but we have perhaps gained insight into the relationship between two previously unconnected rules.

Show[
        GraphicsArray[{RasterGraphics[
              CellularAutomaton[{105, 2, 1}, db29, 300], 
              PlotLabel -> "Rule 105"], 
            RasterGraphics[CellularAutomaton[{150, 2, 1}, db29, 300], 
              PlotLabel -> "Rule 150"]}]];

Replicating and Extending Israeli and Goldenfeld's Experiments

𝒩 (block size) of 2

I can now use this code to determine all elementary rules that are emulated by other rules under the universe of projections for 𝒩 (block size) of 2.

Short[
      t2 = Timing[
          f[{All, 2, 1}, 2, 2] = 
            DeleteCases[f[{All, 2, 1}, db26, 2, 2], {_Integer, {}}]], 
      3]

It would be exciting to find a relationship among elementary cellular automaton, projection, and emulating cellular automaton, because knowledge of such a relationship might accelerate the process of finding emulating cellular automata. The experiments shown below contain some attempts to find a relationship. I find no obvious relationship. I do find a possible hint, however, that there may in fact sometimes be some relationship.

triplets𝒩2 = 
        Flatten[Thread /@ 
              f[{All, 2, 1}, 2, 
                2] /. {e_Integer, {proj_, {n_, k_, r_}}} :> {e, proj, n}, 
          1];

{ca2eca𝒩2OP, ca2proj𝒩2OP, 
          proj2eca𝒩2OP} = {Union[{#[[1]], #[[3]]} & /@ triplets𝒩2], 
          Union[{#[[1]], #[[2]]} & /@ 
              triplets𝒩2], 
          Union[{#[[2]], #[[3]]} & /@ 
              triplets𝒩2]};

As shown in the left panel of the graphic below, there appear to be some intriguing symmetries (vertical around the number 4 and horizontal around the number 4) between the number of 1s in the base 2 representation of the elementary cellular automaton rule and the number of 1s in the base 2 representation of the corresponding emulating cellular automaton rules. As shown in the central panel of the graphic below, there are identical symmetries between the number of 1s in the base 2 representation of the elementary cellular automaton rule and the number of 1s in the base 2 representation of the projections that successfully generate emulating cellular automaton rules. There is only a diagonal symmetry, however, between the number of 1s in the projection and the number of 1 in the emulating cellular automaton rule.

SetOptions[Graphics, Prolog -> PointSize[0.04], AspectRatio -> 1, 
        PlotLabel -> "Block size of 2", Frame -> True];

Show[
        GraphicsArray[{Graphics[
              Point /@ Union[
                  Map[DigitCount[#, 2, 1] &, 
                    ca2eca𝒩2OP, {2}]], 
              FrameLabel -> {"# of 1s in Elem. CA", "# of 1s in Emul. \
CA"}], Graphics[
              Point /@ Union[
                  Map[DigitCount[#, 2, 1] &, 
                    ca2proj𝒩2OP, {2}]], 
              FrameLabel -> {"# of 1s in Elem. CA", "# of 1s in \
Proj."}], 
            Graphics[
              Point /@ Union[
                  Map[DigitCount[#, 2, 1] &, 
                    proj2eca𝒩2OP, {2}]], 
              FrameLabel -> {"# of 1s in Proj.", "# of 1s in Emul. \
CA"}]}]];

I can also find the distribution of the number of 1's in the projection rules that work. The distribution, as might be expected, is symmetric. Also, it appears that projections in which half the inputs are projected onto 1 and half are projected onto 0 have a greater rate of success. The projections are also unbiased in that transforming a 3 to a 1 succeeds as frequently as transforming a 2 to a 1, a 1 to a 1 or a zero to a 1.

Module[{m, dc, tid}, 
m = 
        Map[#[[2]] &, 
          triplets𝒩2]; 
dc = DigitCount[#, 2, 1] & /@ m; 
tid = 
        Transpose[IntegerDigits[#, 2, 4] & /@ 
            m]; 
{Table[Count[dc, i], {i, 0, 4}], 
        Total /@ tid}]

Only 13 of the 10480142147 possible subsets of the 16 possible transformations ever appear to work for 𝒩=2.

Frequencies[Map[First, Last /@ f[{All, 2, 1}, 2, 2], {2}]]

This fact suggests that one might be able to find many (though certainly not all) of the emulations by sampling the successful projections for a subset of source cellular automata and then restricting the search for emulation to a subset of projections that appears particularly promising. Here I draw on a small sample.

{{t2[[1]], Length[t2[[2]]]}, 
      Timing[With[{p = 
              With[{promising = 
                    Map[First, 
                      Last /@ DeleteCases[
                          f[{RandomKSubset[256, 20] - 1, 2, 1}, db26, 2, 
                            2], {_Integer, {}}], {2}]}, Sort[promising, 
                    Length[#1] < Length[#2] &][[1]]]}, 
          Length[DeleteCases[
              f[{All, 2, 1}, db26, 2, 2, 
                Projections -> p], {_Integer, {}}]]]]}

Here I draw on a bigger sample, take a little more time, and find more emulations.

{{t2[[1]], Length[t2[[2]]]}, 
      Timing[With[{p = 
              With[{promising = 
                    Map[First, 
                      Last /@ DeleteCases[
                          f[{RandomKSubset[256, 20] - 1, 2, 1}, db26, 2, 
                            2], {_Integer, {}}], {2}]}, 
                Union[Flatten[
                    Join[Sort[promising, 
                          Length[#1] < Length[#2] &][[{1, 
                            2}]]]]]]}, 
Length[
            DeleteCases[
              f[{All, 2, 1}, db26, 2, 2, 
                Projections -> p], {_Integer, {}}]]]]}

Using Coarse Graining to Emulate Two Color Rules with Three Color Rules

I can see that there are also k=3 rules that emulate the behavior coarse grained version of a k=2 rule. But the set of rules that can be successfully emulated for a block size of two using three color rules is identical to the set of rules that can be successfully emulated for a block size of two using two color rules. So increasing the number of colors does not seem to advance knowledge about the possibilities of emulation.

Timing[
      Short[f[{All, 2, 1}, 2, 3] = 
          DeleteCases[f[{All, 2, 1}, db26, 2, 3], {_Integer, {}}], 
        3]]

First /@ f[{All, 2, 1}, 2, 3] === 
      First /@ f[{All, 2, 1}, 2, 3]

Again, only 25 of the  2.4761288718465865`*^88possible subsets of the 81 possible projections ever appear to work in this scenario.

Frequencies[Map[First, Last /@ f[{All, 2, 1}, 2, 3], {2}]]

𝒩 (block size) of 3

I can use the same code to determine all elementary rules that are emulated by other rules under the universe of projections for 𝒩 (block size) of 3.

Short[
      t3 = Timing[
          f[{All, 2, 1}, 3, 2] = 
            DeleteCases[f[{All, 2, 1}, db29, 3, 2], {_Integer, {}}]], 
      2]

triplets𝒩3 = 
        Flatten[Thread /@ 
              f[{All, 2, 1}, 3, 
                2] /. {e_Integer, {proj_, {n_, k_, r_}}} :> {e, proj, n}, 
          1];

{ca2eca𝒩3OP, ca2proj𝒩3OP, 
          proj2eca𝒩3OP} = {Union[{#[[1]], #[[3]]} & /@ triplets𝒩3], 
          Union[{#[[1]], #[[2]]} & /@ 
              triplets𝒩3], 
          Union[{#[[2]], #[[3]]} & /@ 
              triplets𝒩3]};

As shown in the graphic below, the symmetry patterns observed with block sizes of 2 recur with block sizes of 3.

SetOptions[Graphics, 
        PlotLabel -> "Block size of 3"];

Show[
        GraphicsArray[{Graphics[
              Point /@ Union[
                  Map[DigitCount[#, 2, 1] &, 
                    ca2eca𝒩3OP, {2}]], 
              FrameLabel -> {"# of 1s in Elem. CA", "# of 1s in Emul. \
CA"}], Graphics[
              Point /@ Union[
                  Map[DigitCount[#, 2, 1] &, 
                    ca2proj𝒩3OP, {2}]], 
              FrameLabel -> {"# of 1s in Elem. CA", "# of 1s in \
Proj."}], 
            Graphics[
              Point /@ Union[
                  Map[DigitCount[#, 2, 1] &, 
                    proj2eca𝒩3OP, {2}]], 
              FrameLabel -> {"# of 1s in Proj.", "# of 1s in Emul. \
CA"}]}]];

I can again find the distribution of the number of 1's in the projection rules that work. The distribution, as might be expected, is symmetric. Also, it appears that projections in which half the inputs are projected onto 1 and half are projected onto 0 have a greater rate of success. The projections are also unbiased in that transforming a 3 to a 1 succeeds as frequently as transforming a 2 to a 1, a 1 to a 1 or a zero to a 1.

Module[{m, dc, tid}, 
m = 
        Map[#[[2]] &, 
          triplets𝒩3]; 
dc = DigitCount[#, 2, 1] & /@ m; 
tid = 
        Transpose[IntegerDigits[#, 2, 8] & /@ 
            m]; 
{Table[Count[dc, i], {i, 0, 8}], 
        Total /@ tid}]

Only 48 of the approximately \(TraditionalForm\`7.25027830432137117274515927`15.954589770191005*^373\)possible sets of 256 possible projections ever appear to work for 𝒩=3.

Short[Frequencies[Map[First, Last /@ f[{All, 2, 1}, 3, 2], {2}]], 
      2]

I can again exploit this fact to conduct an accelerated search for at least some emulations.

{{t3[[1]], Length[t3[[2]]]}, 
      Timing[With[{p = 
              With[{promising = 
                    Map[First, 
                      Last /@ DeleteCases[
                          f[{RandomKSubset[256, 20] - 1, 2, 1}, db29, 3, 
                            2], {_Integer, {}}], {2}]}, Sort[promising, 
                    Length[#1] < Length[#2] &][[1]]]}, 
          Length[DeleteCases[
              f[{All, 2, 1}, db29, 3, 2, 
                Projections -> p], {_Integer, {}}]]]]}

𝒩 (block size) of 4

Finally, I can explore what is involved in having Mathematica search for blocks of size 4. The main problem is that with blocks of size 4, there are 4096 possible inputs that need to be explored. And each of the results from these 4096 inputs can be projected in one of (2^((2^4))) = 65,536  ways. And there are (2^8) = 256 possible source rules. So, even if each projection for each rule could be done in an average of 1/100 of a second, it still would take two days to uncover the entire emulation system While I might cut this time perhaps by a factor of 10 with a faster personal computer or a few tricks such as (a) deletion of "isomorphic" projections, (b) exclusion of cellular automata for which an emulating cellular automaton has already been found on smaller block sizes, examination of this large block size does not appear to be well mated to the convenience and elegance of the Mathematica programming environment running on even a fast single machine. Unless a speedier algorithm can be found, this task is better suited either for grid-based computing or for linking to programming languages that gain speed through tighter coupling to the underlying hardware.

In the interim, however, I can use some intelligent sampling to get a feel for the problem. Unfortunately, the sampling method frequently fails to find much of interest.

Timing[
      With[{unfound = 
            Complement[Range[0, 255], 
              Union[Join[First /@ f[{All, 2, 1}, 2, 2], 
                  First /@ f[{All, 2, 1}, 3, 
                      2]]]]}, 
f[{RandomKSubset[unfound, 
              2], 2, 1}, db212, 4, 2, 
          Projections -> Range[0, 32767]]]]

We can use the projections that are found to be successful, however, to find a few other elementary rules that work with block size 4 emulations.

Timing[
      With[{unfound = 
            Complement[Range[0, 255], 
              Union[Join[First /@ f[{All, 2, 1}, 2, 2], 
                  First /@ f[{All, 2, 1}, 3, 
                      2]]]]}, 
DeleteCases[
          f[{unfound, 2, 1}, db212, 4, 2, 
            Projections -> {8, 32767}], {_Integer, {}}]]]

Extending the Code to Totalistic Rules

Two color totalistic cellular automata

I can now see if totalistic k=2 cellular automata can be emulated by other totalistic k=2 cellular automata by projection on block sizes of 2.

Timing[
      Short[f[{All, {2, {1, 1, 1}}, 1}, 2, 2, 
            EmulatingWeights -> {1, 1, 1}] = 
          DeleteCases[
            f[{All, {2, {1, 1, 1}}, 1}, db26, 2, 2, 
              EmulatingWeights -> {1, 1, 1}], {_Integer, {}}], 2]]

totalistictriplets𝒩2 = 
        Flatten[Thread /@ 
              f[{All, {2, {1, 1, 1}}, 1}, 2, 2, 
                EmulatingWeights -> {1, 1, 1}] /. {e_Integer, {proj_, {n_, 
                    k_, r_}}} :> {e, proj, n}, 1];

Union[
      Map[DigitCount[#, 2, 1] &, 
        totalistictriplets𝒩2, {2}]]

{totalisticca2eca𝒩2OP, 
          totalisticca2proj𝒩2OP, 
          totalisticproj2eca𝒩2OP} = {Union[{#[[1]], #[[3]]} & /@ totalistictriplets𝒩2], 
          Union[{#[[1]], #[[2]]} & /@ 
              triplets𝒩2], 
          Union[{#[[2]], #[[3]]} & /@ 
              triplets𝒩2]};

The same symmetries that existed when elementary automata were emulated by elementary automata exist when totalistic automata are emulated by totalistic automata.

SetOptions[Graphics, Prolog -> PointSize[0.08], ImageSize -> 4*72, 
        AspectRatio -> 1, 
        PlotLabel -> "Block size of 2\nTotalistic Automata", 
        Frame -> True];

Show[
        GraphicsArray[{Graphics[
              Point /@ Union[
                  Map[DigitCount[#, 2, 1] &, 
                    totalisticca2eca𝒩2OP, {2}]], 
              FrameLabel -> {"# of 1s in Elem. CA", "# of 1s in Emul. \
CA"}], Graphics[
              Point /@ Union[
                  Map[DigitCount[#, 2, 1] &, 
                    totalisticca2proj𝒩2OP, {2}]], 
              FrameLabel -> {"# of 1s in Elem. CA", "# of 1s in \
Proj."}], 
            Graphics[
              Point /@ Union[
                  Map[DigitCount[#, 2, 1] &, 
                    totalisticproj2eca𝒩2OP, {2}]], 
              FrameLabel -> {"# of 1s in Proj.", "# of 1s in Emul. \
CA"}]}]];

I can again find the distribution of the number of 1's in the projection rules that work. The distribution, as before, is symmetric. Also, it appears that projections in which half the inputs are projected onto 1 and half are projected onto 0 have a greater rate of success. The projections are also unbiased in that transforming a 3 to a 1 succeeds as frequently as transforming a 2 to a 1, a 1 to a 1 or a zero to a 1.

Module[{m, dc, tid}, 
m = 
        Map[#[[2]] &, 
          totalistictriplets𝒩2]; 
dc = DigitCount[#, 2, 1] & /@ m; 
tid = 
        Transpose[IntegerDigits[#, 2, 4] & /@ 
            m]; 
{Table[Count[dc, i], {i, 0, 4}], 
        Total /@ tid}]

Only some sets of transformations ever appear to work for 𝒩=2.

Union[
      Map[First, 
        Last /@ f[{All, {2, {1, 1, 1}}, 1}, 2, 2, 
            EmulatingWeights -> {1, 1, 1}], {2}]]

Length[
      Union[Map[First, 
          Last /@ f[{All, {2, {1, 1, 1}}, 1}, 2, 2, 
              EmulatingWeights -> {1, 1, 1}], {2}]]]

Three color totalistic cellular automata

The Mathematica code developed here also permits exploration of whether 3-color totalistic automata, some of which produce beautiful and complex patterns, can be emulated by other 2-color or 3-color totalistic automata, including perhaps those that produce simpler patterns. I check below, for example to see if totalistic rule 2043 can be emulated with block sizes of 2 by 2-color or 3-color totalistic automata or by 2 or 3 color elementary automata. I find that no emulation exists.

Timing[
      Short[f[{{2043}, {3, {1, 1, 1}}, 1}, 2, 3, 
            EmulatingWeights -> {1, 1, 1}] = 
          DeleteCases[
            f[{{2043}, {3, {1, 1, 1}}, 1}, db36, 2, 3, 
              EmulatingWeights -> {1, 1, 1}], {_Integer, {}}], 3]]

Timing[
      Short[f[{{2043}, {3, {1, 1, 1}}, 1}, 2, 3] = 
          DeleteCases[
            f[{{2043}, {3, {1, 1, 1}}, 1}, db36, 2, 3], {_Integer, {}}], 
        3]]

Other 3-color totalistic automata can be emulated by 2-color or 3-color automata. Rule 689 provides an example. The graphic below shows that the cellular automaton appears to be Class 3.

Show[
        GraphicsArray[{RasterGraphics[
              CellularAutomaton[{689, {3, {1, 1, 1}}, 1}, db36, 500], 3, 
              PlotLabel -> "Rule 689"], 
            RasterGraphics[
              CellularAutomaton[{1093, {3, {1, 1, 1}}, 1}, db36, 500], 3, 
              PlotLabel -> "Rule 1093"]}]];

I find, however, that 3-color totalistic rule 689 can be emulated by several Class 1 cellular automata. For example, take a projection (9817) that maps everything to 1 except for 3, which it maps to zero: proj2rule[9865, 2, 3, 3] = {8 -> 1, 
          7 -> 1, 6 -> 1, 5 -> 1, 4 -> 1, 3 -> 2, 
          2 -> 1, 1 -> 0, 0 -> 1}. I find that Class 1 3-color rule 3812798742493 or Class 1 totalistic rule 1093, which maps everything to 1, emulate rule 689 under this and other projections. From this we can determine that no supercell of 1 will exist in the long run behavior of Rule 689.

Timing[
      Short[f[{{689}, {3, {1, 1, 1}}, 1}, 2, 3] = 
          DeleteCases[
            f[{{689}, {3, {1, 1, 1}}, 1}, db36, 2, 3], {_Integer, {}}], 
        3]]

Timing[
      Short[f[{{689}, {3, {1, 1, 1}}, 1}, 2, 3] = 
          DeleteCases[
            f[{{689}, {3, {1, 1, 1}}, 1}, db36, 2, 3, 
              EmulatingWeights -> {1, 1, 1}], {_Integer, {}}], 3]]

Indeed, I find that the projections that work with Rule 689 succeed in finding 376  3-color totalistic automata that can be emulated by other 3-color totalistic automata. The automata susceptible to successful emulation appear to all be Class 1, 2 or 3.

With[{promising3colortotalisticprojections = 
          First /@ f[{{689}, {3, {1, 1, 1}}, 1}, 2, 3][[1, 2]]}, 
      Length[Map[First, 
          DeleteCases[
            f[{All, {3, {1, 1, 1}}, 1}, 2, 3, EmulatingWeights -> {1, 1, 1}, 
                Projections -> promising3colortotalisticprojections] = 
              f[{All, {3, {1, 1, 1}}, 1}, db36, 2, 3, 
                EmulatingWeights -> {1, 1, 1}, 
                Projections -> 
                  promising3colortotalisticprojections], {_Integer, {}}]]]]

Brute force checking for emulation by other automata under projection with block sizes of 3 proves impracticable. I would have to check (3^(((3 )^3))) = 7625587484987 projections. Even the sampling process described above and even grid computing would have difficulty making a dent in a problem of that magnitude.

Conclusion

The ideas of Israeli and Goldenfeld are potentially important for both physical and other complex systems. If it were possible, for example, easily to find emulating cellular automata that substantially reduced the complexity of source cellular automata, the ability to make useful predictions about at least some complex or chaotic phenomena shown to be "irreducible" by Wolfram might be restored. And this analysis shows Mathematica to prove quite adept as an instrument for presenting and exploring these ideas. Its ability to compute cellular automata rapidly, its ability to create functions of cellular automata, its ability to transform values using projection rules, its support for functional programming constructs such as Fold and FoldList and its unification of ways in which graphical and numeric objects can be programmed, make Mathematica a magnificent programming environment in which to probe these cutting edge ideas. Indeed, the explorations done here suggest (though they certainly do not prove) that there may be certain useful symmetries between characteristics of source rules subject to emulation, projection rules and emulating rules.

Mathematica is not without its limitations, however. As some of the timings above suggest, it takes a long time to examine systems with color counts greater than 2 and block sizes greater than 3. Some disadvantages of programming languages that are abstracted far away from the underlying computer hardware are made apparent by this exercise. Perhaps some of this can be cured with speedier computers, grid computing, more clever Mathematica programming, or development of better algorithms for suggesting the most promising projections for examination, but methods that retain use of Mathematica while accelerating computation by the several orders of magnitude that would be desirable are not readily apparent.

The difficulties in this enterprise may be more than a software issue, however. The reduction in complexity actually achieved by this process of emulation has not yet been shown to be substantial in most cases. Moreover, the process of finding emulating cellular automata may itself be computationally irreducible. On its face, it is certainly a problem subject to computational explosion. The number of projections that need to be studied for any rule grows as ((Subscript[k, β])^((k^𝒩))). The size of the DeBruijn Sequence on which the cellular automata must first be run grows as (2^(3  𝒩)). Thus, even for binary cellular automata emulations of other binary cellular automata, one is starting to look at computations that appear to grow as \(TraditionalForm\`2\^\(2\^𝒩 + 3  𝒩\)\). So, research in this field may ultimately be stymied by one of the issues discussed by Stephen Wolfram in a New Kind of Science: there is indeed a solution – an emulating cellular automaton – it is just essentially impossible to find it. Although the symmetries shown in the relationships among the source rule, the projection and the emulating rule described above may provide tantalizing leads into an algorithm for accelerated calculation of emulating rules, at the moment, the method for compacting or reducing many cellular automata remains hidden, as does any order in the physical, economic, legal, or other system they describe.

The author thanks the Professors Israeli and Goldenfeld for their support and correction in this project. The conclusions and programs are those of the author alone.

References

[Israeli & Goldenfeld 2004] N. Israeli & N. Goldenfeld, Phys. Rev. Lett. 92, 7

[Pemmaraju 2003] S. Pemmaraju, Computational Discrete Mathematics: Combinatorics and Graph Theory with Mathematica, (Cambridge University Press, Cambridge, UK, 2003)

[Wolfram 2003] S. Wolfram, A New Kind of Science (Wolfram Media, Champaign, IL, 2003)