Computing Aggregate Properties of Preimages for 2D Cellular Automata
AAggregate Properties of Preimages
Computing Aggregate Properties of Preimages for 2DCellular Automata
Randall D. Beer a) Cognitive Science Program, Indiana University, Bloomington,IN 47406 (Dated: 13 August 2018)
Computing properties of the set of precursors of a given configuration is a commonproblem underlying many important questions about cellular automata. Unfortu-nately, such computations quickly become intractable in dimension greater thanone. This paper presents an algorithm — incremental aggregation — that cancompute aggregate properties of the set of precursors exponentially faster thanna¨ıve approaches. The incremental aggregation algorithm is demonstrated on twoproblems from the two-dimensional binary Game of Life cellular automaton: pre-cursor count distributions and higher-order mean field theory coefficients. In bothcases, incremental aggregation allows us to obtain new results that were previouslybeyond reach.
As simple models of distributed dynamical systems, cellular automata (CAs)have found wide application across many areas of science. A CA consists ofa regular grid of cells, each of which can be in a finite number of states, thatevolves in discrete steps according to some rule that assigns to each cell anew state based on its previous state and that of its neighbors. Perhaps themost famous such model is Conway’s two-dimensional Game of Life cellularautomaton. Many interesting questions about CAs can be framed as questionsabout the set of precursor configurations that can evolve into a given targetconfiguration under the action of the update rule. Unfortunately, the set ofprecursors can be very expensive to compute. This paper presents a method forcomputing aggregate properties of the set of precursors of a configuration thatis exponentially faster than previous approaches. The method is demonstratedon two problems from the Game of Life.
I. INTRODUCTION
Cellular automata are paradigmatic examples of spatiotemporal complex systems. Theyhave a long history, dating from work by Ulam and von Neumann in the late 1940s and early1950s. Cellular automata have been used to model a wide range of systems, including gases,fluids, excitable media, morphogenesis, tumor growth, the spread of forest fires, traffic flow,urban sprawl, and parallel computation.
In addition, they have served as an importanttestbed for exploring theoretical concepts such as pattern formation, criticality, emergenceand computational universality.
Many interesting questions about a cellular automaton (CA) can be reduced to the com-putation of some property of the set of preimages of configurations of that automaton. TheGarden of Eden problem asks for the configurations of a CA which have no precursors, thatis, which can only occur as initial conditions. Finding the basin of an attractor involvestracing its precursors back to the Garden of Eden states that form its boundary. Statisticalmechanical calculations require partition functions, which involve counting the number ofprecursors with a given property. And so on. Although practical methods exist for com-puting such properties for one-dimensional CAs, there has been much less progress on a) Electronic mail: [email protected]. a r X i v : . [ n li n . C G ] N ov ggregate Properties of Preimages 2higher-dimensional CAs due to the O ( L n d ) scaling of na¨ıve approaches for an L -state size- n CA in d dimensions.This paper describes a new algorithm — incremental aggregation — that can achieve anexponential speedup on problems involving aggregate properties of the entire set of preim-ages of a configuration in two-dimensional cellular automata. After some mathematicalpreliminaries, we present the algorithm and describe its efficient implementation. The algo-rithm is then demonstrated on two problems in Conway’s Game of Life cellular automata. First, we use it to count the number of precursors of individual configurations and to com-pute the distribution of precursor counts over the complete set of configurations of a givengrid size. We show empirically that the algorithm achieves an O (2 n ) to O (2 n ) reduction inexecution time on precursor counting. Second, we compute the coefficients of the 3rd-ordermean field theory for the Game of Life. Along the way, various optimizations that canbe incorporated into the algorithm to further improve its performance are described. Thepaper concludes with a brief discussion of possible extensions to the algorithm and potentialfuture applications. II. PRELIMINARIES A cellular automaton is an autonomous dynamical system (cid:104) S L Σ , T, f t (cid:105) consisting of an L -dimensional state space S L Σ whose elements are drawn from an alphabet of symbols Σ,an ordered time set T (typically Z ∗ ) and an evolution operator f t (with t ∈ T ) definedas the t -times composition of a global transition function f : S L Σ → S L Σ . A key feature ofcellular automata is that f can be written element-wise in a uniform way, so that the globaldynamics unfolds according to the same local law at all points in the space. Note that f implicitly induces both a local and global topology on S L Σ by virtue of the way in which thenext state of one element of S L Σ depends on other elements. The three most common globaltopologies are finite/bounded (the space has an edge beyond which all states are assumedto have some fixed value), finite/unbounded (the space is periodic in all directions) andinfinite.In the Game of Life (GoL) cellular automaton, Σ = { , } and S L Σ has the local topologyof a rectangular grid of cells. The GoL update rule can be written element-wise as s t +1 x,y = δ M, + s tx,y δ M, , (1)where s tx,y is the state of the cell at location ( x, y ) at time t , M is the number of 1 cellsin that cell’s Moore neighborhood (the eight cells surrounding it) and δ i,j is the Kroneckerdelta function (which takes on the value 1 when i = j and the value 0 otherwise). In words,a cell will be ON in the next time step if either exactly three of its neighbors are currentlyON or it and exactly two of its neighbors are currently ON, otherwise it will be OFF. Fornotational simplicity only, we will assume in this paper that the rectangular grid is in factsquare, with L = n × n . We will also assume that the grid is sufficiently large that its globaltopology is irrelevant.Consider some n × n region R of a 2-dimensional grid at time t = 0. This regionsits at the intersection of a forward light cone of cells that it can influence, K + ( R ) = (cid:83) t K t ( R ), extending into the future and a backward light cone of cells that can influenceit, K − ( R ) = (cid:83) t K − t ( R ), extending into the past (Figure 1). A configuration of a region isan assignment of 0 or 1 to each of its constituent cells. If C R is a particular configurationof R , then its preimage f − ( C R ) is the set of all configurations C (cid:48)K − ( R ) of K − ( R ) suchthat f R ( C (cid:48)K − ( R ) ) = C R (where f R denotes f with its codomain restricted to R ). Thisdefinition generalizes straightforwardly to the t -preimage f − t ( C R ) = (cid:110) C (cid:48)K − t ( R ) (cid:12)(cid:12)(cid:12) f t R ( C (cid:48)K − t ( R ) ) = C R (cid:111) . (2)Here, f t denotes the t -times composition of f and f − t denotes the t -times inverse of f .Henceforth, we will drop the R subscripts from our notation, since the relevant regions areusually clear from context.ggregate Properties of Preimages 3Our task in this paper is to efficiently compute a given aggregate property of the set ofpreimages of a configuration. Because this is an exponential problem, the best that wecan hope to do is make the exponent as small as possible. Accordingly, it will be useful tohave some baseline preimage algorithms against which to judge our progress. The simplestsuch algorithm, which we will call the forward algorithm , is to iterate over all possibleconfigurations C (cid:48) of K − t ( R ), evolve each of them forward t times, and then select thosefor which f t ( C (cid:48) ) = C for further processing. Since, due to the Moore neighborhood inGoL, effects spread at a speed of one cell per time, step the number of cells in K − t ( R )is ( n + 2 t ) . Thus, the running time of this algorithm scales as O (2 ( n +2 t ) ). Despite thefact that forward iteration of GoL can be made extremely fast, the scaling behavior of thisforward algorithm severely limits its practical use.A second approach, which we will call the visitation algorithm , is to work backward froma given target pattern, doing a depth-first search through the possible inversions f − ( s x,y )of each cell state s x,y in turn and identifying the globally consistent ones. In the bestcase, the visitation algorithm should scale as O ( | f − t ( C ) | ), although in the worst case it candegenerate to the O (2 ( n +2 t ) ) scaling of the forward algorithm. The visitation algorithm isabout the best we can hope to do if we insist on “touching” each precursor individually. III. THE INCREMENTAL AGGREGATION ALGORITHM
In many applications, we are only interested in aggregate properties of the entire setof preimages of a given configuration. For example, the number of such precursors isan aggregate property. In that case, we can often compute the desired property withoutexplicitly visiting each precursor. The basic idea of the incremental aggregation algorithmis to work row-by-row, combining the desired aggregate property of each partial precursorwith that of each matching set of new row precursors. Henceforth, our notation will assumea radius-1 CA such as GoL for simplicity, although the algorithm can easily be applied tolarger radius CAs. Due to the level of abstraction at which the algorithm is formulated,while reading this section it may be useful to refer to the example of precursor countingdescribed in the next section for a concrete example. Graphical illustrations of the keyoperators and steps of the algorithm are provided in Figures 2 and 3, respectively.Let C be an n × n configuration and let C (cid:48) be one of its ( n + 2) × ( n + 2) precursors.Let C r ,r denote the partial configuration consisting of rows r through r inclusive of C ,with indices ranging from 1 to n . We will call C ,s the s -top of C and C n − s,n its s -bottom.As a special case, C r denotes the r th row of C . We will sometimes need to treat thisrestriction as an operation so that, for example, bottom ( C ,r ) extracts the last row of C ,r ,i.e. C r . A similar notation can be employed for C (cid:48) except that its row indices run from Aggregate Properties of Preimages 2ten element-wise in a uniform way, so that the globaldynamics unfolds according to the same local law at allpoints in the space. Note that f implicitly induces botha local and global topology on S L ⌃ by virtue of the wayin which the next state of one element of S L ⌃ depends onother elements. The three most common global topolo-gies are finite/bounded (the space has an edge beyondwhich all states are assumed to have some fixed value),finite/unbounded (the space is periodic in all directions)and infinite.In the Game of Life (GoL) cellular automaton, ⌃ = { , } and S L ⌃ has the local topology of a rectangular gridof cells. The GoL update rule can be written element-wise as s t +1 x,y = M, + s tx,y M, , (1)where s tx,y is the state of the cell at location ( x, y ) attime t , M is the number of 1 cells in that cell’s Mooreneighborhood (the eight cells surrounding it) and i,j isthe Kronecker delta function (which takes on the value 1when i = j and the value 0 otherwise). In words, a cellwill be ON in the next time step if either exactly threeof its neighbors are currently ON or it and exactly two ofits neighbors are currently ON, otherwise it will be OFF.For notational simplicity only, we will generally assumein this paper that the rectangular grid is in fact square,with L = n ⇥ n . We will also assume that the grid issu ciently large that its global topology is irrelevant.Consider some n ⇥ n region R of a 2-dimensionalgrid at time t = 0. This region sits at the intersec-tion of a forward light cone of cells that it can influ-ence, K + ( R ) = S t K t ( R ), extending into the futureand a backward light cone of cells that can influenceit, K ( R ) = S t K t ( R ), extending into the past (Fig-ure 1). A configuration of a region is an assignment of0 or 1 to each of its constituent cells. If C R is a par-ticular configuration of R , then its preimage f ( C R )is the set of all configurations C ( R ) of K ( R ) suchthat f R ( C ( R ) ) = C R (where f R denotes f with its R K ( R ) K + ( R ) K ( R ) K ( R ) K ( R ) K ( R ) K ( R ) K ( R ) FIG. 1. The forward and backward light cones of a region R (here a single cell) in the Game of Life. codomain restricted to R ). This definition generalizesstraightforwardly to the t -preimage f t ( C R ) = n C t ( R ) f t R ( C t ( R ) ) = C R o . (2)Here, f t denotes the t -times composition of f and f t denotes the t -times inverse of f . Henceforth, we will dropthe R subscripts from our notation, since the relevantregions are usually clear from context.Our task in this paper is to e ciently compute a givenaggregate property of the set of preimages of a configura-tion. Because this is an exponential problem, the bestthat we can hope to do is make the exponent as smallas possible. Accordingly, it will be useful to have somebaseline preimage algorithms against which to judge ourprogress. The simplest such algorithm, which we will callthe forward algorithm , is to iterate over all possible con-figurations C of K t ( R ), evolve each of them forward t times, and then select those for which f t ( C ) = C for fur-ther processing. Since, due to the Moore neighborhoodin GoL, e↵ects spread at a speed of one cell per timestep the number of cells in K t ( R ) is ( n + 2 t ) . Thus,the running time of this algorithm scales as O (2 ( n +2 t ) ).Despite the fact that forward iteration of GoL can bemade extremely fast, the scaling behavior of this forwardalgorithm severely limits its practical use.A second approach, which we will call the visitation al-gorithm , is to work backward from a given target pattern,doing a depth-first search through the possible inversions f ( s x,y ) of each cell state s x,y in turn and identifyingthe globally consistent ones. In the best case, the visita-tion algorithm should scale as O ( | f t ( C ) | ), although inthe worst case it can degenerate to the O (2 ( n +2 t ) ) scal-ing of the forward algorithm. The visitation algorithm isabout the best we can hope to do if we insist on “touch-ing” each precursor individually. III. THE INCREMENTAL AGGREGATION ALGORITHM
In many applications, we are only interested in aggre-gate properties of the entire set of preimages of a givenconfiguration. For example, the number of such precur-sors is an aggregate property. In that case, we can oftencompute the desired property without explicitly visitingeach precursor. The basic idea of the incremental aggre-gation algorithm is to work row-by-row, combining thedesired aggregate property of each partial precursor withthat of each matching set of new row precursors. Hence-forth, our notation will assume a radius-1 CA such asGoL for simplicity, although the algorithm can easily beapplied to larger radii CAs. Due to the level of abstrac-tion at which the algorithm is formulated, while readingthis section it may be useful to refer to the example ofprecursor counting described in the next section for aconcrete example.Let C be an n ⇥ n configuration and let C be oneof its ( n + 2) ⇥ ( n + 2) precursors. Let C r ,r denote FIG. 1. The forward and backward light cones of a region R (here a single cell) in the Game ofLife. ggregate Properties of Preimages 4 ) == ⨁ == bottom( = FIG. 2. A graphical illustration of the key operations used in the incremental aggregation algorithm.Shaded bars represent precursor rows (with identical shading denoting identical contents), androunded squares represent aggregate values (with brighter shades representing larger values). n + 1, so that a precursor of C would be denoted C (cid:48) , . Finally, we denote by tb theoverlapping vertical concatenation of a matching (i.e., bottom ( t ) = top ( b )) pair of a 2-top t and a 2-bottom b to form the corresponding 3-row configuration (Figure 2).Let (cid:104) , (cid:74) , ⊕(cid:105) be an aggregator for a particular domain of aggregation A , consisting ofan identity element , an extension operator (cid:74) which updates an ( r − ⊕ which combines twoaggregate values into one (Figure 2). For example, for precursor counting in Section IV,the identity element will be the integer 0, the extension operator will increment a partialcount each time a valid extension of a partial precursor is encountered, and the aggregationoperator will add two partial counts together.Let B → r A = { . . . , b i → a i , . . . } be a sequence of row-indexed maps from the 2-bottomsof the partial precursors of C ,r to A . These maps index the partial aggregate values wehave already computed up to row r by their 2-bottoms for easy extension. We use thenotation M [[ x ]] to refer to the value associated with the key x in the map M and thenotation M [[ x ]] ← y to denote an assignment to that value. Note that if b does not appearas a key in the bottom-to-aggregate map B → r A , the default value for M [[ x ]] is , i.e., theidentity element for the aggregator.We assume the existence of a set of C r -indexed maps T → C r Bs = { . . . , t j →{ . . . , b k , . . . } , . . . } from the 2-tops of all row precursors C (cid:48) r − ,r +1 of C r to the correspondingsets of matching 2-bottoms. These maps organize the precursors to row C r by their 2-tops,collecting all 2-bottoms that share the same 2-top. These top-to-bottom maps can beprecomputed for any given grid width by either the forward precursor algorithm or thevisitation algorithm. Note that a given 2-bottom b k can be associated with more than one2-top t j because both 2-tops and 2-bottoms individually underdetermine the C (cid:48) r − ,r +1 fromwhich they are derived.The central task of the incremental aggregation algorithm is to compute B → r A from B → r − A using T → C r Bs . This row extension step proceeds as follows. For each pairof entries ( b i → a i ) ∈ ( B → r − A ) and ( t j → { . . . , b (cid:48) k , . . . } ) ∈ ( T → C r Bs ) for which b i = t j , we extend the ( r − a i by computing a i (cid:74) bottom ( b (cid:48) k ).Then we aggregate this value with any previous value a k already associated with b (cid:48) k inthe new map B → r A we are constructing, so that the entry for b (cid:48) k in B → r A becomes a k ⊕ ( a i (cid:74) bottom ( b (cid:48) k )).In order to initialize this process, we need the bottom-to-aggregate map B → A for thefirst row. This initial map could be computed easily if the set of precursors to C wereexplicitly available, but it is not. However, we can reconstruct this set from the associated2-tops and 2-bottoms in T → C r Bs . Coupling this initialization step with a final aggre-gation over all partial results in B → n A produces the incremental aggregation algorithmggregate Properties of Preimages 5(Figure 3). AggregateOverPrecursors ( C, (cid:104) , (cid:74) , ⊕(cid:105) )—– Initialize first row —–( B → A ) ← ∅ for ( t → bs ) ∈ ( T → C Bs ) for b ∈ bs ( B → A )[[ b ]] ← ( B → A )[[ b ]] ⊕ ( (cid:74) tb )—– Process each subsequent row —– for r = 2 to n do ( B → r A ) ← ∅ for ( b → a ) ∈ ( B → r − A ) for b (cid:48) ∈ ( T → C r Bs (cid:48) )[[ b ]]( B → r A )[[ b (cid:48) ]] ← ( B → r A )[[ b (cid:48) ]] ⊕ ( a (cid:74) bottom ( b (cid:48) ))—– Return final aggregate —– result ← ( b → a ) ∈ ( B → n A ) result ← result ⊕ a return result In an efficient implementation of this algorithm, both the B → r A maps and the T → C r Bs maps can be represented as hash tables. Although these hash tables can become large,three optimizations are possible. First, for any given configuration C we need only constructthe T → C r Bs maps for the distinct C r entries that actually appear in C . Second, once B → r − A has been used to construct B → r A , the former map is no longer needed andits storage can be reclaimed. Third, we can avoid constructing the final B → n A map byaggregating directly into the result for the last row.We next examine two applications of the incremental aggregation algorithm. All timingspresented in this paper were performed using an optimized C++ implementation of theincremental aggregation algorithm running on a 3.0GHz 8-core/16-thread 2013 Apple MacPro (Xeon E5-1680 v2) with 64GB of RAM. IV. COUNTING PRECURSORS
One of the simplest aggregate properties that we might want to compute is the numberof preimages of a configuration, (cid:12)(cid:12) f − ( C ) (cid:12)(cid:12) . For precursor counting, the aggregation domain A is Z ∗ and we use the aggregator (cid:104) , , + (cid:105) , where the extension operator 1+ incrementsits aggregate value argument and ignores its row precursor argument. Here the propertywe are aggregating is simply the existence of a precursor and we just sum this property tocompute the total number of such precursors. If we have some number a of ( r − C (cid:48) ,r with a given 2-bottom and we have an allowable extension of that partialprecursor to an r -partial precursor C (cid:48) ,r +1 , then we increment a when we carry it over intothe new map. In addition, since such extensions can come from many different C (cid:48) ,r , wemust combine the new a with any a that already exists in the new map by summation. Ineffect, we are multiplying the number of each C (cid:48) ,r with a given 2-bottom by the number ofpossible extensions of that C (cid:48) ,r to C (cid:48) ,r +1 .This observation leads to an optimization that can be applied to processing the final rowof a target configuration. Not only can we avoid building B → n A , as mentioned above,but we can use multiplication to eliminate the innermost loop of processing for the lastrow. In particular, we terminate the row processing at row n − result ← for ( b → a ) ∈ ( B → n − A ) result ← result + a × length (( T → C n Bs (cid:48) )[[ b ]]) return result As an example of precursor counting, consider the family of n × n vacuum state (all0) configurations 0 n × n . Vacuum states are attractors of GoL and, since there is a strongtendency for a substantial fraction of a GoL grid to quickly decay to quiescence, the basinsof attraction of vacuum states must be quite large. We can use the incremental aggregationalgorithm to compute the 1-basin size of an n × n vacuum configuration by counting thenumber of ( n + 2) × ( n + 2) configurations that evolve to that configuration in one timestep. As can be seen in Figure 4, this number grows very quickly with grid size (blue curve),with (cid:12)(cid:12) f − (0 × ) (cid:12)(cid:12) = 3 965 375 048 845 134 539 385 175 457 630 019 267. In contrast, thenumber of preimages of the family of all 1 configurations 1 n × n (which are not attractorsbut rather 1-precursors of the corresponding vacuum state), grows much more slowly withgrid size (yellow curve), reaching only 1 829 325 441 at n = 10.These two families of configurations also provide an excellent vehicle for probing thescaling behavior of the incremental aggregation algorithm. As shown in Figure 5, thetime it takes this algorithm to compute such counts scales as O (2 cn ), with c ≈ .
79 for (cid:12)(cid:12) f − (0 n × n ) (cid:12)(cid:12) and c ≈ .
03 for (cid:12)(cid:12) f − (1 n × n ) (cid:12)(cid:12) . In both cases, the construction of the required T → C r Bs entries consumes a substantial fraction of this time. Fortunately, this cost can beamortized across multiple runs of the algorithm on configurations that share rows. Althoughstill exponential, this scaling behavior is exponentially better than the O (2 ( n +2) ) scaling ofthe forward algorithm or the best case O ( (cid:12)(cid:12) f − (0 n × n ) (cid:12)(cid:12) ) ≈ O (2 . n ) scaling of the visitationalgorithm. Indeed, given the size of (cid:12)(cid:12) f − (0 × ) (cid:12)(cid:12) , an execution time of about 36 minutesis remarkable; it amounts to an effective speed of almost 2 decillion (10 ) configurationsper second. Of course, the whole point of this algorithm is to avoid having to visit eachprecursor individually.Now we ask a more ambitious question: What form does the full distribution of precursorcounts take for a given grid size? Consider the 5 × × D symmetry of GoL to reduce the number ofconfigurations whose precursors must be counted by roughly a factor of 8. If configurationsare represented as n -bit integers, then these symmetry transformations can be efficientlyimplemented using bitwise techniques such as delta swaps. Second, we can interleave theprocess of scanning through configurations with row-by-row aggregation. That is, for eachpossible configuration of the first row, we construct the corresponding B → A map. Foreach of those, we then iterate over the possible configurations of the second row, constructingthe corresponding B → A maps, and so on until the last row. This ensures that eachbottom-to-aggregate map will be constructed only once and used to its maximum effectbefore being discarded. With these two optimizations, the execution time is reduced to justover 4 minutesA log-log plot of the 5 × × , thesmallest precursor count (40 528) does not correspond to 1 × . Note the small gaps betweenthe precursor counts for these two extreme configurations and the rest of the distribution,suggesting that they are somewhat exceptional. Smaller grids have a similar distribution,although it becomes considerably less smooth as grid size decreases. Monte Carlo samplesof larger grids also suggest that their precursor count distributions take a similar form.ggregate Properties of Preimages 7 V. MEAN FIELD THEORY
More refined aggregate properties of precursors can also be computed by the incrementalaggregation algorithm. Consider GoL mean field theory (MFT), which aims to predict thefuture grid density ρ t ( ρ ) from the initial density ρ . The first-order theory was derivedby Schulman and Seiden (blue curve in Figure 7). This theory predicts two stable fixedpoints at ρ ∞ = 0 and ρ ∞ ≈ . ρ ∞ ≈ . ρ ∞ ≈ . derivedan exact higher-order MFT that takes into account local correlations and they computedits second-order coefficients (yellow curve in Figure 7). This theory predicts only a singlestable fixed point at ρ ∞ = 0 and it is quite useful for explaining several other aspects of theshort-term evolution of GoL grids. Although even higher-order MFTs would provide a moreaccurate tool, Bagnoli et al. point out that computing their coefficients seems intractable.Next we show how the incremental aggregation algorithm can be used to efficiently computethe coefficients of the third-order theory MFT3.Higher-order MFTs can be derived as follows. Given the symmetries of GoL, the expecteddensity after t timesteps of a grid with initial density ρ is equivalent to the expected state ofan individual cell after t time steps averaged over the set of (2 t + 1) × (2 t + 1) configurationswith density ρ that determine that state. Thus, if we define (cid:12)(cid:12) f − (1) (cid:12)(cid:12) k to be the Hammingweight decomposition of (cid:12)(cid:12) f − (1) (cid:12)(cid:12) (that is, the number of distinct configurations in f − (1)that contain exactly k ρ t ( ρ ) = E (cid:34) (cid:12)(cid:12) f − (1) (cid:12)(cid:12) k (cid:0) (2 t +1) k (cid:1) (cid:35) , (3)where k ∼ B ((2 t + 1) , ρ ). Expanding this expectation using the PDF of the binomialdistribution and cancelling the normalization factor, we obtain ρ t ( ρ ) = (2 t +1) (cid:88) k =0 (cid:12)(cid:12) f − t (1) (cid:12)(cid:12) k ( ρ ) k (1 − ρ ) (2 t +1) − k . (4)The t h-order MFT is thus a polynomial in ρ defined by its coefficient vector c t withcomponents c tk = | f − t (1) | k . This vector can be computed in three steps. First, the visitationalgorithm is applied recursively to determine each symmetrically-distinct configuration C ∈ f − ( t − (1) and its multiplicity. Then we use the incremental aggregation algorithm tocompute the coefficient vector of each such C . Finally, we scale these vectors by theircorresponding multiplicities and sum them to obtain c t .For the computation of MFT coefficients, the aggregation domain is the space of coefficientvectors Z (2 t +1) +1 and the aggregator is (cid:104) , (cid:86) , + (cid:105) . Here, is the 0 vector of length (2 t +1) + 1 and + is vector addition. The extension operator (cid:86) updates the coefficient vectorfor the ( r − c i in c has a value v for an( r − w , then coefficient c i + w will have value v in the resulting r -partial precursor. Notethat, for this to work, (cid:86) x must be defined to set the component of corresponding tothe Hamming weight of x to 1.For MFT3, there are 1 065 921 symmetrically-distinct 5 × f − (1) and c can be calculated in just over two hours (Table I). As Bagnoli et al. argued, the firstthree and at least the last five coefficients are necessarily 0 due to the nature of the GoLupdate rule. The resulting ρ ( ρ ) curve is shown in green in Figure 7. It is similar to thesecond-order curve (yellow), but differs from it in two small but important ways. First,ggregate Properties of Preimages 8 k (cid:12)(cid:12) f − (1) (cid:12)(cid:12) k k (cid:12)(cid:12) f − (1) (cid:12)(cid:12) k (cid:12)(cid:12) f − (1) (cid:12)(cid:12) k of the third-order MFT ρ ( ρ ). its peak is lower. In fact, simulations suggest that, as the order increases, the maximumof ρ t ( ρ ) monotonically decreases toward its asymptotic value of ρ ∞ ≈ . Simulations suggest that the distance totangency steadily increases for higher orders, so MFT3 appears to be exceptional in thisregard.What are the prospects for computing even higher-order MFT coefficients? Even thenumber of symmetrically-distinct (2 t +1) × (2 t +1) configurations in f − t (1) grows extremelyquickly with t . Aside from this scaling, the most time-consuming portion of the computationis the dynamic memory allocation and deallocation of coefficient vectors in the B → r A hash tables. This could probably be reduced significantly by using a more specializedstatically-allocated data structure, by using a per-thread memory allocator, or by splittingthe computation across multiple processes rather than multiple threads in order to avoidmemory allocation contention. VI. CONCLUSIONS
This paper has presented a general method, the incremental aggregation algorithm, forcomputing aggregate properties of the set of preimages of a configuration in 2-dimensionalcellular automata. This algorithm is exponentially faster than standard approaches becauseit scales with the number of rows rather than grid area. The algorithm was demonstratedon two problems in the Game of Life: precursor counting and coefficient calculation forhigher-order mean field theories.The incremental aggregation algorithm can be extended in many ways. With the appro-ggregate Properties of Preimages 9priate aggregator it could be applied to other properties, such as finding Garden of Edenstates or computing the basins of attraction of a given configuration to a desired depth.In addition, since nothing in the algorithm is specific to the Game of Life, it could also beapplied easily to other cellular automata rules. For example, it would be interesting to seehow the mean field theory of other “Life-like” CAs compares to that for GoL. Finally, withsome tedious but straightforward work, the algorithm could be applied to non-square andeven non-rectangular grids and to CAs with more states, more dimensions, larger radii andother topologies. As always, the main limitation is the scaling of the problem relative tothe available computational resources. ACKNOWLEDGMENTS
This paper has benefited greatly from discussions with Alexander Gates. B. Chopard, “Cellular automata modeling of physical systems,” in
Encylopedia of Complexity and SystemsScience , edited by R. A. Meyers (Springer, 2009) pp. 865–892. L. B. Kier and P. G. Seybold, “Cellular automata modeling of complex biochemical systems,” in
Encylo-pedia of Complexity and Systems Science , edited by R. A. Meyers (Springer, 2009) pp. 848–865. T. Worsch, “Cellular automata as models of parallel computation,” in
Encylopedia of Complexity andSystems Science , edited by R. A. Meyers (Springer, 2009) pp. 741–755. J. E. Hanson, “Emergent phenomenon in cellular automata,” in
Encylopedia of Complexity and SystemsScience , edited by R. A. Meyers (Springer, 2009) pp. 768–778. J. Durand-Lose, “Universality of cellular automata,” in
Encylopedia of Complexity and Systems Science ,edited by R. A. Meyers (Springer, 2009) pp. 901–913. E. Jen, “Enumeration of preimages in cellular automata,” Complex Systems , 421–456 (1989). A. Wuensche and M. Lesser,
The Global Dynamics of Cellular Automata (Addison Wesley, 1994). E. R. Berlekamp, J. H. Conway, and R. K. Guy,
Winning Ways for Your Mathematical Plays, Volume2 (A. K. Peters, 1982). A. Adamatzky, ed.,
Game of Life Cellular Automata (Springer, 2010). K. Sutner, “On the computational complexity of finite cellular automata,” Journal of Computer andSystem Sciences , 87–97 (1995). D. E. Knuth,
The Art of Computer Programming, Volume 4A: Combinatorial Algorithms, Part 1 (Ad-dison Wesley, 2011). L. S. Schulman and P. Seiden, “Statistical mechanics of a dynamical system based on Conway’s Game ofLife,” Journal of Statistical Physics , 293–314 (1978). F. Bagnoli, R. Rechtman, and S. Ruffo, “Some facts of Life,” Physica A , 249–264 (1991). P. Bak, K. Chen, and M. Creutz, “Self-organized criticality in the Game of Life,” Nature , 780–782(1989). H. J. Blok and B. Bergersen, “Effect of boundary conditions on scaling in the ‘Game of Life’,” PhysicalReview E , 6249–6252 (1997). S. M. Reia and O. Kinouchi, “Conway’s Game of Life is a near-critical metastable state in the multiverseof cellular automata,” Physical Review E , 052123–052125 (2014). ggregate Properties of Preimages 10 { } { { } { } } { }( ) ⨁ ⨁ ( ) B ! r A : B ! r A : T ! C r Bs : { { } { } } { } ! ( ) ⨁ ⨁ ! ( ) T ! C Bs : B ! A : { } B ! n A : ⨁ ⨁ FIG. 3. A graphical illustration of of the operation of the incremental aggregation algorithm.Conventions are the same as in Figure 2. (Top) The initialization step. (Middle) The row extensionstep. (Bottom) The final aggregation step. ggregate Properties of Preimages 11 n N u m be r o f P r e c u r s o r s FIG. 4. A log plot of the scaling of the number of precursors with grid size for 0 n × n (blue)and 1 n × n (yellow) in the Game of Life. The fits are 2 ∧ (0 . n + 3 . n + 4 . ∧ (2 . n + 4 . - - - n t i m e ( s e cs ) FIG. 5. A log plot of the execution time scaling of precursor counting with grid size for 0 n × n (blue)and 1 n × n (yellow) in the Game of Life. The fits are 2 ∧ (2 . n − . ∧ (2 . n − . ggregate Properties of Preimages 12 - - - - - - - Number of Precursors P r obab ili t y FIG. 6. A log-log plot of the distribution of precursor counts for 5 × ρ ρ t FIG. 7. Exact first- (blue), second- (yellow) and third- (green) order mean field theories ρ t ( ρ0