GCANs: Global Cellular Automaton Networks

© Seth J. Chandler 2004
SChandler@uh.edu

Introduction

Many social, economic, biological and legal interactions can be thought of as a situation in which actors have private information, public signals that they give off, and various algorithms for figuring out which signals to emit. Often those algorithms involve introspection – a look at one's own condition and values – and an examination of the behavior of "neighboring" members of society. This article conceptualizes something known as a Global Cellular Automaton Network (GCAN) as a general tool for exploring systems that feature distributed information (with varying degrees of privacy), potentially heterogeneous values, and complex patterns of connection between agents. It does so using the Mathematica programming language and drawing on some of the insights of Stephen Wolfram's book, A New Kind of Science. The GCAN derives from the somewhat simpler construct of a global cellular automaton ("GCA"). A GCA is a single cellular automaton in which the evolution of each site is a function not only of the values of the site's neighbors – which is the distinguishing feature of conventional cellular automata – but also of certain global features of the automaton. A GCAN preserves the notion of evolution that depends in some way on global features of the system but further extends from both a cellular automaton and a GCA in at least two ways: (1) there are multiple cellular automata that are networked with each other; (2) each site within each cellular automaton within the GCAN evolves based not only the values of its neighbors and on its own global characteristics, but also on certain characteristics of the cellular automata to which it is connected. This extension creates a tool capable of great flexibility in the systematic study of complex systems. Moreover, GCANs can be used to explore "Random Cellular Networks," in which the update rule to be applied by each site on each iteration is determined by application of a node-specific cellular automaton to the signals being emitted by its neighboring node. They provide a vehicle for studying composition of cellular automata. This notebook begins the study of GCANs and shows how they can be implemented in Mathematica.Third,  The code for this versatile mechanism appears to be relatively speedy and almost as general as Mathematica's underlying CellularAutomaton function.

Implementing Global Cellular Automata

Introductory Code

I adopt a convention in this article of "graying out" code that is inessential, such as cosmetic directives to graphics routines.

Needs["DiscreteMath`Combinatorica`"]; 
    Needs["Statistics`DescriptiveStatistics`"]; 
    Needs["Graphics`Graphics`"]; Off[General::spell1];

RasterGraphics[state_?MatrixQ, colors_Integer : 2, opts___] := 
      Graphics[Raster[Reverse[1 - state/(colors - 1)]], 
        AspectRatio -> (AspectRatio /. {opts} /. 
              AspectRatio -> Automatic), opts]

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

SeedRandom[1];

Discussion

I first create a function called GCA that takes a list of potential cellular automaton rules that might be applied, a system against which a particular cellular automaton rule is to be applied, and a switching function for choosing, based on the state, which of the potential rules to apply.  I then use the built-in Mathematica CellularAutomaton command to evolve the system to its next iteration. This implements a globally controlled cellular automaton: the rule that is to be applied is a function of the state of the entire system. (This code derives from an October 14, 2003  post by Stephen Wolfram to the "Pure NKS" forum at www.wolframscience.com .)

\(GCA[rlist_List, \ init_List, \ 
          f_]\  := \(CellularAutomaton[rlist[\([f[init]]\)], \ init, \ 

            1, \ {\(-1\)}]\)[\([1]\)]\)

Here are some examples of how it operates. The expression below takes a 100 site random array and applies elementary rule 110 or elementary rule 250 to all the sites. The system determines which rule to apply by examining the length of runs of 0s or 1s. If there are runs of 0s or 1s at least one fifth as long as the number of sites, then rule 110 is applied. Otherwise rule 250 is applied.  The size of the runs is determined by the auxiliary longestcircularrun function.

longestcircularrun[x_] := 
      Max[Length /@ 
          With[{s = Split[x]}, 
            If[s[[1, 1]] === 
                s[[-1, 1]], 
              s /. {α_, β___, ω_} :> {Join[α, ω], β}, s]]]

Short[
      GCA[{110, 250}, Array[Random[Integer] &, {400}], 
        If[longestcircularrun[#] ≥  Length[#]/5, 1, 2] &], 
      2]

Notice that one can then NestList the GCA command to evolve the system from some starting configuation. The output can be visualized in the customary way.

gcaexample1 = 
        NestList[
          Function[system, 
            GCA[{110, 250}, system, 
              If[longestcircularrun[#] ≥  Length[#]/5, 1, 
                  2] &]], Array[Random[Integer] &, {400}], 400];

Show[RasterGraphics[gcaexample1], 
        AspectRatio -> 1];

The basic GCA function is quite flexible. There can be any number of rules, the rules can be of different types: they can take different number of colors and different neighborhoods; they can even be functions. The switching function that is the last argument to GCA can be any function that produces an integer output bound between 1 and the length of list of possible rules.

Implementing Global Cellular Automata Networks

I next create an extension of GCA that undertakes parallel processing of networked globally controlled cellular automata (GCAN). The idea is to have a network of cellular automata each of which can apply varying rules depending on the environment. The choice of each cellular automaton as to which rule to apply depends not only on the global state of its own system, but also on the state of neighboring systems. These neighboring systems themselves co-evolve based on their own state and that of their neighbors. The Mathematica code for implementing this proves relatively simple. The first argument is a list containing lists of rules that each global cellular automaton can potentially applies. Thus if the first GCA can apply rules 23, 45 or 52 depending on the output of its switching function and the second GCA can apply rules 110, 150, and 223 depending on the output of its switching function, the first argument would read {{23,45,52},{110,150,223}}. The second argument is a list containing the initial state of each system, what I will call the "meta-state" of the system.  So, if the first GCA started in configuaration {0,0,1,0,0,0,1,0,0} and the second system started in configuration {0,1,0}, the second argument (the meta-state) would read {{0,0,1,0,0,0,1,0,0},{0,1,0}}. The third argument is a list of dependencies. So, in a system with five parallel cellular automata in which the first system depended on the third and fourth systems, the second system depended only on the second system, the third system depended on the first and fifth systems, the fourth system depended on the fifth system and the fifth system depended on the second system, the "adjacency list" would read {{3,4},{2},{1,5},{5},{2}}. (This representation of the dependencies facilitates interfacing with Mathematica's powerful Combinatorica package, which can likewise represent graphs in this fashion.) The fourth argument is a list of "meta-functions" each of which produces a switching function. The meta-functions do so by applying themselves to the metastate of the system (i.e. the state of all the cellular automata in the system) and the list of all neighbors of the node associated with the meta-function. That process produces a list of functions (one for each node) which can then be applied to the state of each of the corresponding nodes. By way of example, consider the meta-function metaf set forth below, which will employ the auxiliary restrictknowledge and middlecell functions set forth below. For a particular node, the metaf function looks at the sum of the middle cells of its neighbors. If that sum is even, the value of the switch is one plus the value of that node's first cell. If that sum is odd, the value of the switch is one minus the value of that node's first cell.

restrictknowledge[metastate_, alsublist_, restrictionfunction_] := 
      restrictionfunction /@ Part[metastate, alsublist]

middlecell[state_] := 
      state[[Quotient[Length[state], 2, -2]]]

Those concerned about the detailed operations of metafunctions are urged to look at the example below. When we  execute the command metaf below, we obtain 1 + s1[1]. We do so because our "restrict knowledge" of neighboring nodes 2 and 3 is {1,1}. We will call that resulting list s. We then use the ca23 auxilliary function to run CellularAutomaton with elementary rule 23 on s over Max[1, 
        Floor[((Length[s] - 𝒩)/(2 r))]] iterations, where r is the radius of the cellular automaton rule being applied and where 𝒩 is the number of digits we want output. This formula for the number of iterations insures that all sites of s have an opportunity to influence the the value of the 𝒩 digits extracted. When we do this in our example, we find that the middle cell of the last iteration of the cellular automaton is 0. We thus perform the computation 1 + First[self]. First[self] is s1[1], so the result is 1+s1[1]. This method lets the GCAN look at properties of neighboring nodes (including itself) and then permits that data to effect the evaluation of its own global state. And all of this is done to determine which update rule that node is to apply on its succeeding turn.

ca23[s_, 𝒩_, r_] := 
      CellularAutomaton[23, s, 
        Max[1, Floor[(Length[s] - 𝒩)/(2  
                  r)]], {-1, 0}]

metaf = 
        Function[{init, alsublist}, 
          Function[self, 
            With[{s = 
                  restrictknowledge[init, alsublist, 
                    middlecell], 𝒩 = 1, 
                r = 1}, 
If[
                ca23[s, 𝒩, r] === {{0}}, 1 + First[self], 
                1 + (1 - First[self])]]]];

system = {Array[s1, {7}], {1, 0, 1, 0, 0}, {1, 1, 1}, {0, 0, 0, 0, 0, 
            0, 1}};

metaf[system, {2, 3}][system[[1]]]

The piece of code below generalizes the GCAN procedure. GCAN and GCA thus become the two pieces of code crucial to implement the concept of a Global Cellular Automaton Network.

GCAN[rlist : List[__List], metastate : List[__List], 
        adjacencylist_List, f_List] := 
      MapThread[
        GCA[#1, #2, #4[metastate, #3]] &, {rlist, metastate, adjacencylist, 
          f}]

Using GCANs To Model Real World Complex Systems

The Concept

The most obvious use of GCANs is to study models of phenomena such as political or judicial systems in which individuals have their own private beliefs regarding candidates or doctrine but make tentative voting decisions based on the signals emitted by those they regard as "neighbors." GCANs can also be used to study an investment market composed of individuals who have private information but who give off public signals (such as purchasing decisions). The GCAN unifies these areas by permitting exploration of the co-evolution of systems in which the rule each cellular automaton (agent) uses to update itself depends in part on the signals emitted by its neighboring nodes. The idea is that each node looks at its own data and the signals being emitted by its neighbor and then uses those signals to determine which update rule to apply to that data. The application of the update rule determines the new state of each node and the signal it is emitting.

In the example below, the signal given off by each node (voter, market participant, computer) is arbitrarily determined to be the value of its middle cell. The middle cell is simply a window into the otherwise private system of the node. Each node then determines the modal value of the middle cells of the nodes to which it is connected. (Ties are arbitrarily resolved in favor of the highest valued state).  Each node then applies the rule in its rule set corresponding to the modal value of the signal being emitted by its neighbors. Thus, if the middle cells of the systems in the neighborhood of a particular system are {0,1,1,2,2,1,1,0}, then the system will apply its second rule to its data on the succeeding iteration. Why? It does so because 1 is the modal value of the middle cells. (For those confused by the reference to the second rule, remember that, although the values of the system start with 0, the rules are indexed starting from 1. ) A system that has no neighbors applies its first rule to its data. The GCAN mechanism we have developed makes this system fairly simple to implement.

mostvotes[
          votes : List[__]] := Sort[
            Frequencies[
              votes]][[-1, -1]];
mostvotes[{}] := 0

The randomnetwork defined as illustrated below shows a possible set of connections between a small set of voters or judges or other decision makers.

ShowLabeledGraph[FromAdjacencyLists[
            randomnetwork = ToAdjacencyLists[RandomGraph[10, 0.4]]],AspectRatio -> 1/2];

Each node can hold a different amount of data. Here the nodes each hold between 5 and 15 bits or information.

randominitialstates = 
        Table[With[{q = Random[Integer, {5, 15}]}, 
            Table[Random[Integer, {0, 1}], {q}]], {Length[
              randomnetwork]}];

Homogeneous Rule Sets

A Sample Homogeneous Rule Set

I can start with a random GCAN in which all sites (systems) choose from the same set of rules.

homogeneousrulesets = 
        Table[{30, 45}, {Length[randomnetwork]}];

I now evolve the systems over 100 iterations and hold the results in the variable networksystem.

networksystem = 
        NestList[
          Function[metastate, 
            GCAN[homogeneousrulesets, metastate, randomnetwork, 
              Table[Function[{init, al}, 
                  Function[self, 
                    With[{votes = Map[middlecell, init[[al]]]}, 
                      1 + mostvotes[votes]]]], {Length[randomnetwork]}]]], 
          randominitialstates, 400];

I can examine the evolution of the private data of all of the systems.

Show[GraphicsArray[MapIndexed[RasterGraphics[#1,PlotLabel -> "System\n" <> 
                          ToString[#2[[1]]],AspectRatio -> 8]&,Transpose[networksystem, {2, 1}]]]];

There is an issue about how to determine the periodicity of a GCAN. One way would be to determine the the periodicity of the set of signals given off by the agents in the GCAN. This method can sometimes be somewhat misleading, however, because the return of the set of signals to a given earlier condition does not absolutely guarantee that the future path of the signals will replicate those of the path: the "private" condition of each of the constituent automata may influence future signals too. Nonetheless, for simplicity, I examine "signal periodicity" here.

period[z_] := 
      Replace[z, {{a___, b_, c___, b_, e___} :> 
            Length[{b, c}], {a__, b_} :> 
            Interval[{1 + Length[{a, b}], ∞}]}]

period[networksystem]

signalperiodicity[gcan_] := 
      period[Transpose[
          Map[middlecell, Transpose[gcan, {2, 1}], {2}]]]

signalperiodicity[networksystem]

Or I can examine aggregations concerning the signals the system is emitting. (One might think of this as how the market or electorate is "voting" over time.) Here, I look at how many agents are voting for "1." Notice how this relatively simple system can create deterministic economic or political markets with highly complex and essentially unpredicatable dynamics.

Show[ListPlot[Total /@ 
            Transpose[
              Map[middlecell, Transpose[networksystem, {2, 1}], {2}]],PlotJoined -> True,AspectRatio -> 1/3,AxesLabel -> {"Time", "Votes for 1"}]]

The Set of All Homogeneous Rule Sets

For a given random network, one can look at issues of periodicity over all possible homogeneous rule sets. Here, for example is the distribution of the periods of the signals emitted by the agents for all 256 homogeneous elementary totalistic rule sets. That computation takes on the order of a minute. The graphic below shows the distribution of periodicity over all homogeneous elementary totalistic rule sets for our random network but over four different random initial conditions. I develop some auxilliary functions to aid in the graphic presentation.

ns[network_, initialstates_, iters_, ruleset_] := 
        NestList[
          Function[metastate, 
            GCAN[ruleset, metastate, network, 
              Table[Function[{init, al}, 
                  Function[self, 
                    With[{votes = Map[middlecell, init[[al]]]}, 
                      1 + mostvotes[votes]]]], {Length[network]}]]], 
          initialstates, iters];

freq[x_, set_] := 
      Module[{fx, c}, 
        fx = Frequencies[x]; 
c = 
          Partition[Complement[set, Last /@ fx], 1, 
            1]; 
Replace[
          Sort[Join[c, fx], 
            Or[And[Not[NumberQ[#2[[-1]]]]], 
                Last[#1] < Last[#2]] &], {z_} -> {0, 
              z}, {1}]]

periodicityfrequency[pt_, iters_] := 
      freq[Flatten[Map[Last, pt, {2}]], 
        Append[Range[0, iters + 1], 
          Interval[{iters + 2, ∞}]]]

barlabels[x_List, blankout_List] := 
      ReplacePart[Last /@ x, "", Partition[blankout, 1]]

Show[GraphicsArray[Table[With[{pf = 
                        periodicityfrequency[
                          Table[{i, j, 
                              signalperiodicity[
                                ns[randomnetwork, 
                                  Table[
                                    With[{q = Random[Integer, {5, 15}]}, 
                                      Table[
                                        Random[
                                        Integer, {0, 1}], {q}]], {Length[
                                        randomnetwork]}], 19, 
                                  Table[{{i, {2, {1, 1, 1}}, 
                                        1}, {j, {2, {1, 1, 1}}, 1}}, {Length[
                                        randomnetwork]}]]]}, {i, 0, 15}, {j, 
                              0, 15}], 19]},BarChart[pf,DisplayFunction -> Identity,PlotRange -> {Automatic, {0, 150}},\(AxesLabel -> {"Period", "Frequency"}\),TextStyle -> {FontSize -> 16, 
                              FontFamily -> "Times"},BarLabels -> 
                            barlabels[pf, 
                              Drop[Complement[Range[2, Length[pf]], 
                                  Range[2, Length[pf], 10]], -1]]]],{2},{2}],PlotLabel -> 
                StyleForm["Distribution of Periods\nFor Different Random \
Initial States", FontSize -> 18, FontFamily -> "Times"],GraphicsSpacing -> 0]];

One can also see whether the pattern of connections in the random network significantly affects the distribution of periods. It turns out also generally to have little effect. Most random networks generate of this size end up creating signals with short periods, though a few create signals with long periods.

Show[GraphicsArray[Table[With[{pf = 
                        periodicityfrequency[
                          Table[{i, j, 
                              signalperiodicity[
                                ns[ToAdjacencyLists[RandomGraph[10, 0.4]], 
                                  randominitialstates, 19, 
                                  Table[{{i, {2, {1, 1, 1}}, 
                                        1}, {j, {2, {1, 1, 1}}, 1}}, {Length[
                                        randomnetwork]}]]]}, {i, 0, 15}, {j, 
                              0, 15}], 19]},BarChart[pf,DisplayFunction -> Identity,PlotRange -> {Automatic, {0, 150}},\(AxesLabel -> {"Period", "Frequency"}\),TextStyle -> {FontSize -> 16, 
                              FontFamily -> "Times"},BarLabels -> 
                            barlabels[pf, 
                              Drop[Complement[Range[2, Length[pf]], 
                                  Range[2, Length[pf], 10]], -1]]]],{2},{2}],PlotLabel -> 
                StyleForm["Distribution of Periods\nFor Different Random \
Networks", FontSize -> 18, FontFamily -> "Times"],GraphicsSpacing -> 0]];

Finally, one can vary the probability of connectedness among the agents and determine what effect it has on the distribution of periods of signals emitted by the agents. Again, the distribution of periodicity remains quite similar over the specified range.

Show[GraphicsArray[
Table[
With[{pf = 
                        periodicityfrequency[
                          Table[
{i, j, 
                              With[
{rn = 
                                    ToAdjacencyLists[
                                      RandomGraph[10, 
                                        FromDigits[{x, y}, 3]/10. ]]}, 
                                signalperiodicity[
                                  ns[rn, randominitialstates, 19, 
                                    Table[{{i, {2, {1, 1, 1}}, 
                                        1}, {j, {2, {1, 1, 1}}, 1}}, {Length[
                                        rn]}]]]]}, 
{i, 0, 
                              15}, {j, 0, 15}], 19]},
BarChart[pf,DisplayFunction -> Identity,PlotRange -> {Automatic, {0, 150}},\(AxesLabel -> {"Period", "Frequency"}\),PlotLabel -> 
                            ToString[FromDigits[{x, y}, 3]/10. ],TextStyle -> {FontSize -> 16, 
                              FontFamily -> "Times"},BarLabels -> 
                            barlabels[pf, 
                              Drop[Complement[Range[2, Length[pf]], 
                                  Range[2, Length[pf], 10]], -1]]]],{x, 2},{y, 2}],
PlotLabel -> 
                StyleForm["Distribution of Periods\nFor Different Prob. of \
Connection", FontSize -> 18, FontFamily -> "Times"],GraphicsSpacing -> 0]];

Heterogeneous Rule Sets

Obviously, one can have more complex systems in which each node not only has its own private data, but also has its own set of update rules. Consistent with Wolfram's notion of a principle of computational equivalence, this increase in the complexity of the model often results in no increase in the complexity of the system's outputs. Indeed, sometimes the output can be less complex. Here the variable ruleheterogeneousnetworksystem holds the results of evolving this GCAN.

ruleheterogeneousnetworksystem = 
        With[{heterogeneousrulesets = 
              Array[Random[Integer, {0, 255}] &, {Length[randomnetwork], 
                  2}]}, NestList[
            Function[metastate, 
              GCAN[heterogeneousrulesets, metastate, randomnetwork, 
                Table[Function[{init, al}, 
                    Function[self, 
                      With[{votes = Map[middlecell, init[[al]]]}, 
                        1 + mostvotes[votes]]]], {Length[randomnetwork]}]]], 
            randominitialstates, 150]];

The graphic below shows the evolution of each of the systems.

Show[
        GraphicsArray[
          MapIndexed[
            RasterGraphics[#1, 
                PlotLabel -> "System\n" <> ToString[#2[[1]]], 
                AspectRatio -> 8] &, 
            Transpose[ruleheterogeneousnetworksystem, {2, 1}]]], 
        ImageSize -> 432];

The graphic below shows the evolution of the number of nodes "voting" for one over time.

ListPlot[
        Total /@ Transpose[
            Map[middlecell, 
              Transpose[ruleheterogeneousnetworksystem, {2, 1}], {2}]], 
        PlotJoined -> True, AspectRatio -> 1/3, 
        AxesLabel -> {"Time", "Votes for 1"}];

One can see whether adding heterogeneity to the model affects the distribution of periods. If heterogeneity in the rule sets adds to "complexity," I would expect to see that fact reflected in a significant, rightwards shift in the distribution of periods. And, indeed, that is precisely what we see. The graphic below examines 256 randomly selected rule sets. A substantial number develop long periods.

SetOptions[BarChart,DisplayFunction -> Identity,PlotRange -> {Automatic, {0, 150}},AxesLabel -> {"Period", "Frequency"},TextStyle -> {FontSize -> 16, 
                FontFamily -> "Times"}];

Show[GraphicsArray[Table[With[{pf = 
                        periodicityfrequency[
                          Table[{i, j, 
                              signalperiodicity[
                                ns[randomnetwork, randominitialstates, 19, 
                                  Array[
                                    Random[Integer, {0, 15}] &, {Length[
                                        randomnetwork], 2}]]]}, {i, 0, 
                              15}, {j, 0, 15}], 19]},BarChart[pf,BarLabels -> 
                            barlabels[pf, 
                              Drop[Complement[Range[2, Length[pf]], 
                                  Range[2, Length[pf], 10]], -1]]]],{1},{1}],PlotLabel -> 
                StyleForm["Distribution of Periods\nFor Heterogeneous Rule \
Sets", FontSize -> 10, FontFamily -> "Times"],GraphicsSpacing -> 0]];

Past the Complexity Threshhold

One can further "complexify" the model by permitting the nodes to give off two-bit signals instead of single bit signals and then giving each node four rules to choose from instead of two. Again, however, this increase in realism, ends up resulting in little additional complexity in the overall system. Here, I show the evolution of a system in which each node has seven neighbors but different amounts of private data.

bigregularnetwork = 
        ToAdjacencyLists[RegularGraph[7, 50]];

Each node has its own set of update rules to use on its private data.

verycomplexrulesets = 
        Array[Random[Integer, {0, 255}] &, {Length[bigregularnetwork], 
            4}];

Each site has different amounts of private data.

complexinitialstates = 
        Table[With[{q = Random[Integer, {5, 15}]}, 
            Table[Random[Integer, {0, 1}], {q}]], {Length[
              bigregularnetwork]}];

middletwocells[state_] := 
        With[{q = Quotient[Length[state], 2, -2]}, 
          Take[state, {q - 1, q}]];

We now evolve this more complex system.

verycomplexnetworksystem = 
        NestList[
          Function[metastate, 
            GCAN[verycomplexrulesets, metastate, bigregularnetwork, 
              Table[Function[{init, al}, 
                  Function[self, 
                    With[{votes = 
                          Map[FromDigits[middletwocells[#], 2] &, 
                            init[[al]]]}, 
                      1 + mostvotes[votes]]]], {Length[
                    bigregularnetwork]}]]], complexinitialstates, 
          150];

The graphic below shows the evolving percentage of the nodes issuing each of the four signals.

PercentileBarChart[
        Sequence @@ (Transpose[
              Map[First, 
                Map[freq[#, {0, 1, 2, 3}] &, 
                  Map[FromDigits[middletwocells[#], 2] &, 
                    Take[verycomplexnetworksystem, 40], {2}]], {2}]]), 
        Axes -> None, Frame -> True, 
        FrameLabel -> {"Time", "Voting"}, 
        FrameTicks -> {None, Automatic}, 
        BarEdgeStyle -> Thickness[0.00001], BarSpacing -> 0.4, 
        ImageSize -> 432];

Although the system is deterministic, there is no particular pattern that can be observed. The graphic below confirms this by showing the number of nodes issuing "signal 3" for all 150 iterations of the system.

ListPlot[
        Last[(Transpose[
              Map[First, 
                Map[freq[#, {0, 1, 2, 3}] &, 
                  Map[FromDigits[middletwocells[#], 2] &, 
                    verycomplexnetworksystem, {2}]], {2}]])], 
        PlotJoined -> True, AspectRatio -> 1/3, 
        AxesLabel -> {"Time", "Votes for 1"}];

Random Cellular Networks

We can also use the GCAN methodology to evolve systems of "Random Cellular Networks" in which the switching function is simply the output of a random cellular automaton on the "signals" emitted by the neighboring cellular automata. The cellular automaton iterates the number of times necessary to ensure that all cells have the opportunity to influence the middle cell. The name and methodology is similar to Stuart Kauffman's "Random Boolean Network" except that here a cellular automaton, rather than an arbitrary boolean function, updates the system.

cametaf[ca : {n_, k_, r_}, 𝒩_:  1] := 
        Function[{init, al}, 
          Function[self, 
            With[{s = restrictknowledge[init, al, middlecell]}, 
              If[CellularAutomaton[ca, s, 
                    Max[1, Floor[(Length[s] - 𝒩)/(2  
                              r)]], {-1, 0}] === {{0}}, 1, 
                2]]]];

RCNnetwork = ToAdjacencyLists[RandomGraph[100, 0.1]];

RCNrulesets = 
        Array[Random[Integer, {0, 255}] &, {Length[RCNnetwork], 
            2}];

RCNfunctions = 
        Table[cametaf[{Random[Integer, {0, 255}], 2, 1}], {Length[
              RCNnetwork]}];

RCNinitialstates = 
        Table[With[{q = Random[Integer, {10, 100}]}, 
            Table[Random[Integer, {0, 1}], {q}]], {Length[
              RCNnetwork]}];

Timing[RCNnetworksystem = 
          NestList[
            Function[metastate, 
              GCAN[RCNrulesets, metastate, RCNnetwork, RCNfunctions]], 
            RCNinitialstates, 150];]

The following graphic shows the internal state of a few of the nodes over time.

Show[
        GraphicsArray[
          MapIndexed[
            RasterGraphics[#1, 
                PlotLabel -> "System\n" <> ToString[#2[[1]]], 
                AspectRatio -> 8] &, 
            Take[Transpose[RCNnetworksystem, {2, 1}], 10]]]];

And the graphic below shows the number of "one" signals being emitted by the 100 nodes in the system.

ListPlot[
        Total /@ Transpose[
            Map[middlecell, Transpose[RCNnetworksystem, {2, 1}], {2}]], 
        PlotJoined -> True, AspectRatio -> 1/3, 
        AxesLabel -> {"Time", "Votes for 1"}];

We have great flexibility in the sort of cellular automata we use to determine the value of the "switch." Here is an example in which the switches are k=2, r=2 automata.

RCNfunctionsr2 = 
        Table[cametaf[{Random[Integer, {0, (2^((2^5))) - 1}], 2, 2}], {Length[
              RCNnetwork]}];

Timing[RCNnetworksystemr2 = 
          NestList[
            Function[metastate, 
              GCAN[RCNrulesets, metastate, RCNnetwork, RCNfunctionsr2]], 
            RCNinitialstates, 150];]

Show[
        GraphicsArray[
          MapIndexed[
            RasterGraphics[#1, 
                PlotLabel -> "System\n" <> ToString[#2[[1]]], 
                AspectRatio -> 8] &, 
            Take[Transpose[RCNnetworksystemr2, {2, 1}], 10]]], 
        ImageSize -> 432];

Moreover, instead of assigning the cellular automaton randomly at the outset, one could make the update rule endogenous. In the example below, the cellular automaton function used to "reduce" the set of signals from the neighboring nodes to a single signal is determined by the base 255 sum of each site's own values. Again, notwithstanding the complexity of the overall process, the GCAN code works quite swiftly.

endogenousfunction = 
        Function[{init, al}, 
          Function[self, 
            With[{s = restrictknowledge[init, al, middlecell]}, 
              If[CellularAutomaton[Mod[Total[self], 255], s, 
                    Max[1, Floor[(Length[s] - 1)/(2 1)]], {-1, 
                      0}] === {{0}}, 1, 2]]]];

endogenousfunctiontable = 
        Table[endogenousfunction, {Length[RCNnetwork]}];

Timing[endogenousRCNnetworksystem = 
          NestList[
            Function[metastate, 
              GCAN[RCNrulesets, metastate, RCNnetwork, 
                endogenousfunctiontable]], RCNinitialstates, 
            150];]

Here is a sample of the evolution of the nodes of the resulting system.

Show[
        GraphicsArray[
          MapIndexed[
            RasterGraphics[#1, 
                PlotLabel -> "System\n" <> ToString[#2[[1]]], 
                AspectRatio -> 8] &, 
            Take[Transpose[endogenousRCNnetworksystem, {2, 1}], 
              10]]]];

Emulating Composed Systems with GCANs

GCANs have many uses. A GCAN can emulate a composed system in which the update rules proceed in some fixed periodic sequence. The system may be said to be temporally inhomogeneous and periodic. The emulation occurs simply by having one of the systems implement a counter. Here is one in which the first system applies Rule 110 or Rule 250 depending on the state of the middle cell of the second system (the onlyneighbor function is used to implement this), and the second system applies Rule {1,2,0} (a two-color automaton of radius 0) regardless of the state of the first system (the onenomatterwhat function is used to implement this). Since the middle column of a system is initiated as {1} and updated using Rule {1,2,0}, it value simply flips back and forth between 0 and 1. The first system thus alternates application of Rule 110 and Rule 250.

onlyneighbor = 
        Function[{init, al}, 
          Function[self, 
            1 + First[
                restrictknowledge[init, al, middlecell]]]];

onenomatterwhat = 
        Function[{init, al}, Function[self, 1]];

emulatedcomposedsystem = 
        With[{r100 = Table[Random[Integer, {0, 1}], {100}]}, 
          NestList[
            Function[metastate, 
              GCAN[{{110, 250}, {{1, 2, 0}}}, 
                metastate, {{2}, {1, 1}}, {onlyneighbor, 
                  onenomatterwhat}]], {r100, {1}}, 200]];

Show[
        GraphicsArray[
          MapIndexed[
            RasterGraphics[#1, 
                PlotLabel -> "System " <> ToString[#2[[1]]], 
                AspectRatio -> 1] &, 
            Transpose[emulatedcomposedsystem, {2, 1}]]], 
        AspectRatio -> 1/2, ImageSize -> 432];

"Composing" and "Decomposing" Cellular Automata Using GCANs

The ability of GCANs to readily compose cellular automata can be used to explore (a) whether composition can create greater complexity than that existing in the constituent rules,and (b) the relationship between the complexity generated by a cellular automaton and the ability of that cellular automatona to be decomposed into constituent simpler rules.

The function compositionca, which works somewhat analogously to Mathematica's built-in Composition function, takes a sequence of cellular automaton specifications and the initial condition of the system and produces the terminal condition of the system.

compositionca[ca__][ic_List] := 
      With[{z = Length[{ca}]}, 
        First[Nest[
            Function[metastate, 
              GCAN[{{ca}, {{FromDigits[RotateLeft[Range[z - 1, 0, -1]], 
                        z], z, 0}}}, 
                metastate, {{2}, {1}}, {onlyneighbor, 
                  onenomatterwhat}]], {ic, {z - 1}}, z]]]

Thus, if we (1) construct the DeBruijn Sequence (see http://mathworld.wolfram.com/deBruijnSequence.html) (a) over the alphabet employed by the cellular automata being composed and (b) of a length equal to the diameter of a cellular automaton that might emulate the behavior of the composed sequence, and then (2) run the composed cellular automata on that sequence, we can determine the relationship between all possible beginning states of the hypothesized emulating cellular automaton and the resulting end states.  I can then write some Mathematica code that takes these starting and ending values and determines the identity of any cellular automaton whose behavior that transformation might emulate. The code handles elementary and weighted cellular automata, such as totalistic automata. The Mathematica compiler is used to increase speed.

startendform2βform[{start_List, end_List}, r_:  1] := 
      MapThread[Join, 
        List[Partition[start, 2  r + 1, 1, {r + 1, r + 1}], 
          Partition[end, 1, 1, {1, 1}]]]

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}}];

Options[βform2CAform] = {EmulatingWeights -> Elementary, 
          EmulatingRadius -> 1};

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

Composition of k=2, r= 1/2 Cellular Automata Using GCANs

A little experimentation with Mathematica shows that, of the 16 k=2, r=1/2 cellular automata, few are capable of generating more than the simplest behaviors. When started from seed on an infinite lattice, none generate anything more complex than nesting behavior. When started from random conditions, only two, Rule 6 and its mirror Rule 9, are generally capable of preserving randomness. The rest move swiftly to homogeneous states or systems with very short periods. And when started from seed on some finite lattices, only the same two, Rules 6 and 9, are capable of generating complexity. So, we wish to see whether the composition of rules excluding Rule 6 and Rule 9 is capable of generating anything interesting. We also wish to see whether any of the elementary  (k=2, r=1) cellular automata known to produce particularly complex behavior can be represented as the composition of two simpler k=2, r=1/2 cellular automata. GCANs are a simple way of exploring these issues.

The table below shows the way in which two k=2, r=1/2 automata can be composed to produce a k=2, r=1 automaton. The composition of 0 and 1, for example, yields elementary rule 255.

Short[
      Timing[k2r12table = 
          With[{db23 = DeBruijnSequence[{0, 1}, 3]}, 
            Table[{i, 
                j, βform2CAform[
                    startendform2βform[{db23, compositionca[{j, 
                              2, {{0}, {-1}}}, {i, 2, {{0}, {1}}}][
                          db23]}, 1], 2, 2][[1]]}, {i, 0, 15}, {j, 0, 
                15}]]], 3]

DisplayForm[
      GridBox[Join[{Flatten[{" ", ToString /@ Range[0, 15]}]}, 
          MapIndexed[Prepend[#1, #2[[1]] - 1] &, 
            Map[StyleForm[ToString[Last[#]], FontSize -> 16] &, 
              k2r12table, {2}]]], ColumnWidths -> 3, 
        ColumnAlignments -> Right]]

Notice that not all elementary cellular automata can be "factored" into these r=1/2 automata. Here is a list of the "unfactorable" elementary automata.

Complement[Range[0, 255], 
      Union[Flatten[Map[Last, k2r12table, {2}]]]]

If one looks at the composed systems that do not have Rule 6 or 9 as one of their constituents, one finds that none are capable of generating interesting behavior from seed on an infinite or finite lattice; none are capable of preserving any randomness in the initial conditions; and none are even capable of generating nesting. Moreover, even among the composed systems that do have Rule 6 or 9 as one of their constituent parts, none are capable of producing behavior more complex than nesting when started from seed on an infinite lattice. Thus, the elementary cellular automata with truly complex behavior (such as Rule 30 or Rule 110) appear not to be subject to decomposition into smaller k=2, r=1/2 cellular automata. At least for k=2, r=1/2 cellular automata, composition does not appear to bring Wolfram's boundary of complex behavior any closer.

k2r12excluding6or9 = 
      Union[Last /@ 
          DeleteCases[
            Flatten[k2r12table, 1], {___, 6  |  9, __}]]

Composition of k=2, r=1 Totalistic Cellular Automata using GCANs

I can do the same sort of experiment on the totalistic cellular automata. None of the 16 k=2, r=1 totalistic cellular automata are capable of generating anything more complex than nesting when started from seed on an infinite lattice. Five rules (2, 5, 6, 9, 10, and 11), however, are capable of (a) generating complex patterns when started from seed on a finite lattice and (b) capable of preserving randomness in the initial conditions. GCANs can help study the k=2, r=2 cellular automata that result from the composition of these simpler totalistic CAs. Generally, these composed systems will not themselves be totalistic.

Short[
      Timing[k2r1ttable = 
          With[{db25 = DeBruijnSequence[{0, 1}, 5]}, 
            Table[{i, 
                j, βform2CAform[
                    startendform2βform[{db25, compositionca[{j, {2, {1, 1, 1}}, 1}, {i, {2, {1, 1, 1}}, 1}][db25]}, 2], 2, 2, 
                    EmulatingRadius -> 2][[1]]}, {i, 0, 15}, {j, 0, 
                15}]]], 3]

Although the entire table is too large to fit conveniently on a conventional sheet of paper, we can provide a flavor of the results by displaying part of the resulting composition table.

DisplayForm[
      GridBox[Join[{Flatten[{" ", ToString /@ Range[0, 8]}]}, 
          MapIndexed[Prepend[#1, #2[[1]] - 1] &, 
            Map[StyleForm[ToString[Last[#]], 
                  FontSize -> 16] &, Take[#, 9] & /@ 
                k2r1ttable, {2}]]], ColumnWidths -> 9, 
        ColumnAlignments -> Right]]

Wecan also determine the identity of any totalistic cellular automata that are the composition of two k=2, r=1 totalistic cellular automata. In the table below "x" means that the system is not a totalistic cellular automaton.

\(Short[
      Timing[k2r1tttable = 
          With[{db25 = DeBruijnSequence[{0, 1}, 5]}, 
            Table[{i, j, 
                With[{β = βform2CAform[
                        startendform2βform[{db25, \(compositionca[{j, {2, {1, 1, 1}}, 1}, {i, {2, {1, 1, 1}}, 1}]\)[db25]}, 2], 2, 2, 
                        EmulatingRadius -> 2, 
                        EmulatingWeights -> {1, 1, 1, 1, 1}]}, 
                  If[Head[β] === 
                      String, "x", β[\([1]\)]]]}, {i, 0, 15}, {j, 
                0, 15}]]], 3]\)

DisplayForm[
      GridBox[Join[{Flatten[{" ", ToString /@ Range[0, 15]}]}, 
          MapIndexed[Prepend[#1, #2[[1]] - 1] &, 
            Map[StyleForm[ToString[Last[#]], FontSize -> 16] &, 
              k2r1tttable, {2}]]], ColumnWidths -> 5, 
        ColumnAlignments -> Right]]

To see if complexity can emerge from the composition of two automata not themselves capable of much complexity, I temporarily restrict attention to those 48 composed systems no part of which contains a k=2, r=1 totalistic system that is capable of producing more than Wolfram Class 1 or Wolfram Class 2 behavior, i.e. I delete systems in which one of the constituent rules is Rule 2, 5, 6, 9, 10 or 11. The excluded rules have the ability to preserve randomness.  I find that, when started from seed, none of the composed systems produce anything more complex than nesting. When started from a random configuration, however, many of the systems are Class 3; some are possibly Class 4.  An example (the composition of totalistic 14 with totalistic 4) is shown below. Thus, even though neither of the constituent cellular automata are generally capable of preserving randomness, the composed system is capable of doing so.

simplecomposedtotalistic = 
        With[{a = Alternatives[2, 5, 6, 9, 10, 11]}, 
          Flatten[DeleteCases[k2r1ttable, {a, _, _}  |  {_, a, _}, {2}], 
            1]];

Show[
        RasterGraphics[
          CellularAutomaton[{Cases[
                  simplecomposedtotalistic, {14, 4, x_} :> x][[1]], 2, 
              2}, Table[Random[Integer, {0, 1}], {400}], 400]]];

complexcomposedtotalistic = 
        With[{a = Alternatives[2, 5, 6, 9, 10, 11]}, 
          Cases[k2r1ttable, {a, _, _}  |  {_, a, _}, {2}]];

I can also see if the composition of k=2, r=1 Class 3 or Class 4 totalistic automata ever lead to behavior of a different class. There are many instances in which the composition of these more complex automata lead to simple behavior. The composition of Rule 5 with Rule 11 leads, for example, to Class 1 behavior. And, as shown by the perturbation results displayed below, the composition of Rule 10 with Rule 9 results in Rule 2286166545, which produces behavior that may will be Class 4.

Show[
        RasterGraphics[
          With[{ra = Array[Random[Integer, {0, 1}] &, {301}]}, 
            Abs[CellularAutomaton[{Cases[
                        k2r1ttable, {10, 9, _}, {2}][[1, -1]], 2, 
                    2}, ra, 300, {All, All}] - 
                CellularAutomaton[{Cases[
                        k2r1ttable, {10, 9, _}, {2}][[1, -1]], 2, 
                    2}, MapAt[1 - # &, ra, {{151}}], 
                  300, {All, All}]]]]];

Conclusion

This paper implements the concept of a Global Cellular Automaton Network in Mathematica. The mechanism developed is capable of studying both composition of cellular automata in which some of the nodes essentialy serve only as "counters," and more complex systems in which the nodes interact with each other. The latter system provides a parsimonious yet sophisticated account of how complex human and computational systems may interact with each other. The versatility of the GCAN mechanism also permits it to be used in the study of "Random Cellular Networks," in which cellular automata themselves determine the update rule to be applied at each site. Further research is needed to elucidate issues such as (a) the relationship between the connection patterns of the nodes and the complexity of its evolution; (b) the actual GCANs that might emulate a particular system of human interaction; and (c) whether GCANs might be useful in conducting practical computations.