8
\$\begingroup\$

I have written code to take a cellular automaton configuration and determine a state that could have existed one time-step prior, according to Game of Life rules.

The algorithm goes roughly as follows:

  1. Take a cell, and iterate through every possible 3x3 configuration of cells that could've evolved into this cell after one time step.
  2. If the configuration does not produce the correct state (i.e. live or dead), check the next.
  3. If the configuration does not overlap with the other proposed prior configurations in the area, check the next.
  4. If no configurations work, go back to the prior cell and recurse: continue iterating through possible cell configurations.
  5. If the configuration fits, put it on the grid. Go to the next cell.
from collections import deque
from itertools import product

class CellularAutomaton:
    configs = tuple(product((False, True), repeat=9))

    def __init__(self, rows, cols, fill=None):
        self.rows = rows
        self.cols = cols
        self.size = rows * cols

        """Cells can be False=dead, True=alive, None=undetermined. """
        self.g = [[fill] * rows for _ in range(cols)]

    def __getitem__(self, i):
        r, c = i
        return self.g[r % self.rows][c % self.cols]

    def __setitem__(self, i, v):
        r, c = i
        self.g[r % self.rows][c % self.cols] = v

    def __eq__(self, t):
        """
        Compare two automata.
        """
        return all(self[i] == t[i] for i in self)

    def __iter__(self):
        yield from product(range(self.rows), range(self.cols))

    def neighborhood(self, i):
        """
        Return the indices of all adjacent cells and cell itself.
        """
        r, c = i
        return (
            (r-1,c-1), (r-1,c), (r-1,c+1),
            (r  ,c-1), (r  ,c), (r  ,c+1),
            (r+1,c-1), (r+1,c), (r+1,c+1)
        )

    def index(self, i):
        """
        Converts 1-dimensional indices.
        """
        return (i // self.cols, i % self.cols)

    def reverse(self):
        """
        Return a new CellularAutomaton that evolves into this one after
        one evolution.
        """
        rows, cols, size = self.rows, self.cols, self.size
        ret = CellularAutomaton(rows, cols)

        """Hypothesis: a stack to keep track of which cell is using which
        configuration. """
        hypo = deque([iter(CellularAutomaton.configs)], maxlen=size+1)

        """A stack to keep track of which cells are changed with each
        configuration so if a configuration needs to be undone, we know which
        cells to revert. """
        undo = deque([]  , maxlen=size)

        while len(hypo)-1 < size:
            i = ret.index(len(hypo)-1)

            for cfg in hypo[-1]:
                nbhd = ret.neighborhood(i)

                """Does the configuration produce the right state? """
                if self[i] != ((cfg[4] and sum(cfg) in (3, 4)) or (not cfg[4] and sum(cfg) == 3)):
                    continue

                """Does the configuration fit on the automaton with previously-
                decided configurations? """
                if any(ret[n] != c and ret[n] is not None for n, c in zip(nbhd, cfg)):
                    continue

                """Only add undetermined cells to the undo stack because all
                other cells' states were determined by other nearby
                configurations. """
                undo.append(tuple(n for n in nbhd if ret[n] is None))

                """Update the return board with the configuration. """
                for n, c in zip(nbhd, cfg):
                    ret[n] = c

                hypo.append(iter(CellularAutomaton.configs))
                break
            else:
                """We iterated through every configuration and none of them
                worked. We need to go back to a previous cell. """

                hypo.pop()
                for c in undo.pop():
                    ret[c] = None

                continue

        return ret

    def forward(self):
        """
        Evolve one time-step forward.
        """
        rows, cols = self.rows, self.cols
        ret = CellularAutomaton(rows, cols)

        for i in self:
            nbhd = sum(self[j] for j in self.neighborhood(i)) - self[i]
            ret[i] = True if nbhd == 3 else (self[i] if nbhd == 2 else False)

        return ret

"""This is a smiley face. """
end = CellularAutomaton(10, 10, False)
end[2,2] = 1
end[2,3] = 1
end[7,2] = 1
end[7,3] = 1
end[1,5] = 1
end[2,6] = 1
end[3,7] = 1
end[4,7] = 1
end[5,7] = 1
end[6,7] = 1
end[7,6] = 1
end[8,5] = 1

assert end == end.reverse().forward(), "Not equal."

An important constraint of this problem is that I intend to use it on automata with dimensions potentially on the order of thousands of cells, so I don't believe I can use actual recursion here due to memory constraints. Instead, I am faking recursion with stacks.

I am looking for advice to make this code faster. It takes several minutes to reverse some configurations on my machine by just a single time-step. I am more concerned here about the result than the process, so I welcome algorithm changes, language/library recommendations, or even large-scale changes such as different automata systems, different boundary conditions, or other deviations from the project scope.

I am not very interested in style or organization improvements, but will welcome comments on these matters nonetheless.

\$\endgroup\$
1
  • 2
    \$\begingroup\$ Very interesting problem you have here. I instantly thought about combinatorics. See my comment on the accepted answer for how I would have approached the problem. I think you can mainly assume that dead areas were dead in the previous iteration, although it is possible with overpopulation, if you would want to go two steps back that would be weird :) \$\endgroup\$ Commented Apr 27, 2023 at 18:15

2 Answers 2

7
\$\begingroup\$

Reduce the search space

For each cell, it looks like the code always considers all 512 possible 9x9 configurations for the previous generation. However, if the cell is currently alive, you know that the previous configuration must have had 3 or 4 live cells. There are only 126 of those configurations, so only considering those would reduce the search space by about 75%. Similarly, if the cell is dead, there are only 386 possible configurations for the previous generation, a 25% reduction.

Here are the changes. Organize configs so that the configurations that would result in a live cell are in configs[1] and those that would result in a dead cell are in configs[0]:

class CellularAutomaton:
    configs = ([],[])
    
    for t in product((False, True), repeat=9):
        s = sum(t)
        if s == 3 or s == 4 and t[4]:
            configs[1].append(t)
        else:
            configs[0].append(t)

Then, in reverse(), the first thing to push on the hypo stack, are the configs that could result in cell[0,0]:

hypo = deque([iter(CellularAutomaton_A.configs[self[0, 0]])], maxlen=size+1)

The test to make sure the chosen config produces the right state can be eliminated, because it's precomputed in the way configs was built.

Lastly, after updating the return board with the configuration, we push on the stack the possible configurations for the next cell based on it's current state:

nr, nc = ret.index(len(hypo))
hypo.append(iter(CellularAutomaton_A.configs[self[nr, nc]]))

On my machine, the original code takes about 3.6 seconds to run the example provided. With the simple changes above, the code takes about 2.3 seconds...about 33% faster!

Reduce it some more.

The neighborhoods of adjacent cells overlap. This can be used to further prune the search space. Consider cells 6 and 7 below:

1  2  3  4
5  6  7  8
9 10 11 12

When cell 6 has a possible configuration, it has proposed values for cells 2, 3, 6, 7, 10 and 11 which are in the neighborhood of cell 7. When pushing possible configurations for cell 7, we know the current state of cell 7 and the proposed configuration for cell 6. Given what is known, there are at most 8 possibilities for cells 4, 8, and 12. That's 8 possible configurations instead of 512, a 98% reduction.

Let's call the tuple for cells 2, 3, 6, 7, 10, 11 the "right_key" of cell 6. It is also the "left_key" of cell 7. The only valid configs for cell 7 are configs that have a left_key that matches the right_key of cell 6.

As before, configs is first indexed by the current state of a cell, but each part is now a dict. The keys to the dict are the left_key of of a configuration and the values are list of configurations having the corresponding left_key.

One tricky part is that cells in column 0, don't have a cell on the left to get a key from, so in that case we push an iterable over all possibilities that can result in the current value of that cell. This is done using itertools.chain.from_iterable on the values() of one of the dicts in configs.

Here is code based on this idea:

from collections import defaultdict, deque
from itertools import chain, product

def right_key(t):
    return t[1], t[2], t[4], t[5], t[7], t[8]

def left_key(t):
    return t[0], t[1], t[3], t[4], t[6], t[7]

class CellularAutomaton_B:
    configs = (defaultdict(list), defaultdict(list))
    
    for t in product((False, True), repeat=9):
        s = sum(t)
        k = left_key(t)
        if s == 3 or s == 4 and t[4]:
            configs[1][k].append(t)
        else:
            configs[0][k].append(t)
            
    def __init__(self, rows, cols, fill=None):
        self.rows = rows
        self.cols = cols
        self.size = rows * cols

        """Cells can be False=dead, True=alive, None=undetermined. """
        self.g = [[fill] * rows for _ in range(cols)]

    def __getitem__(self, i):
        r, c = i
        return self.g[r % self.rows][c % self.cols]

    def __setitem__(self, i, v):
        r, c = i
        self.g[r % self.rows][c % self.cols] = v

    def __eq__(self, other):
        """
        Compare two automata.
        """
        return self.g == other.g

    def __iter__(self):
        yield from product(range(self.rows), range(self.cols))

    def neighborhood(self, i):
        """
        Return the indices of all adjacent cells and cell itself.
        """
        r, c = i
        return (
            (r-1,c-1), (r-1,c), (r-1,c+1),
            (r  ,c-1), (r  ,c), (r  ,c+1),
            (r+1,c-1), (r+1,c), (r+1,c+1)
        )

    def index(self, i):
        """
        Converts 1-dimensional indices.
        """
        return (i // self.cols, i % self.cols)

    def reverse(self):
        """
        Return a new CellularAutomaton that evolves into this one after
        one evolution.
        """
        rows, cols, size = self.rows, self.cols, self.size
        ret = CellularAutomaton_B(rows, cols)

        """Hypothesis: a stack to keep track of which cell is using which
        configuration. """
        initial_hypos = chain.from_iterable(CellularAutomaton_B.configs[self[0, 0]].values())
        hypo = deque([initial_hypos], maxlen=size+1)
        
        """A stack to keep track of which cells are changed with each
        configuration so if a configuration needs to be undone, we know which
        cells to revert. """
        undo = deque([], maxlen=size)

        while len(hypo) <= size:
            i = ret.index(len(hypo)-1)

            for cfg in hypo[-1]:
                nbhd = ret.neighborhood(i)

                """Does the configuration fit on the automaton with previously-
                decided configurations? """
                if any(ret[n] != c and ret[n] is not None for n, c in zip(nbhd, cfg)):
                    continue

                """Only add undetermined cells to the undo stack because all
                other cells' states were determined by other nearby
                configurations. """
                undo.append(tuple(n for n in nbhd if ret[n] is None))

                """Update the return board with the configuration. """
                for n, c in zip(nbhd, cfg):
                    ret[n] = c

                nr, nc = ret.index(len(hypo))
                if nc:
                    key = right_key(cfg)
                    hypo.append(iter(CellularAutomaton_B.configs[self[nr, nc]][key]))
                    
                else:
                    first_col_hypo = chain.from_iterable(CellularAutomaton_B.configs[self[nr, nc]].values())
                    hypo.append(iter(first_col_hypo))
                break
            else:
                """We iterated through every configuration and none of them
                worked. We need to go back to a previous cell. """

                hypo.pop()
                for c in undo.pop():
                    ret[c] = None

                continue

        return ret

    def forward(self):
        """
        Evolve one time-step forward.
        """
        rows, cols = self.rows, self.cols
        ret = CellularAutomaton_B(rows, cols)

        for i in self:
            nbhd = sum(self[j] for j in self.neighborhood(i)) - self[i]
            ret[i] = True if nbhd == 3 else (self[i] if nbhd == 2 else False)

        return ret

This code runs the example in about 0.6 seconds, or 1/10 the time of the original code.

More?

That first cell needs to consider all possible configurations. The other cells in the first row have information about the configuration of the cell to their right, which information is used in the program above. The cells in the first column have information about the cell above it, which could be used in a similar manner. The cells in the rest of the grid have information about the cell above it and the cell to it's left. Using both leaves only 2 configurations to consider (the to possible values of the bottom, right cell). Extending the code is left as an exercise for the reader.

\$\endgroup\$
3
  • \$\begingroup\$ This is great stuff. \$\endgroup\$
    – SirPython
    Commented Apr 24, 2023 at 20:06
  • 1
    \$\begingroup\$ When I first read this question, I thought about solving it by playing Minesweeper. Although you haven't formulized it as such here, it feels like the same concepts are being applied. This problem is just a more generalized version where you know that e.g. 2 <= a + b + c + d + e <= 3 \$\endgroup\$ Commented Apr 27, 2023 at 18:06
  • \$\begingroup\$ @SimonForsberg the original code is basically a depth first search of potential solutions. However, it didn't use all the information at hand to reduce the search space. I showed that the current value of a cell and a proposed configuration of nearby cells limit the possible choices at each point in the search. Data can be pre-calculated, so that it isn't necessary to test if configurations of adjacent cells are compatible because the pre-calculated data structure only returns compatible configurations. But in the end it's still a DFS. \$\endgroup\$
    – RootTwo
    Commented Apr 28, 2023 at 1:32
3
\$\begingroup\$

representation

Target CA's are square grids with sides of e.g. 10², 10³, 10⁴ cells.

        self.g = [[fill] * rows for _ in range(cols)]

This works. But it uses lots of memory, and all those object pointers induce random reads.

Prefer to use a python array, or perhaps a numpy ndarray. Then we can store a million uint8's in approximately a megabyte, with the overhead of just a single object rather than a million of them. And we can do cache-friendly sequential scans of DRAM instead of random reads that are hard for Intel to predict.


algorithm

This is essentially about representation. Let's assume you adopted the ndarray approach and are able to easily switch between dense and sparse datastructures.

Consider those three grid sizes, the largest having a hundred million cells. We will stick to simulating in just the .forward() direction. Put a glider gun on each. Should we suffer N^2 quadratic performance? No!

Position a Gosper glider gun on each of those three grids at, IDK, (80, 80), so that gliders are launched toward the origin where they disappear. Is the N=1000 CA really a harder problem than for N=100 ? No. Most cells are dead and will stay that way. We need to avoid visiting them. Similarly for N=10,000.

You need a way of identifying active / inactive bounding boxes, changing the bbox sizes, and creating / merging them. There is a robust literature on this topic, some of it motivated by the desire to run simulations forward at high speed.

Here's another motivating example. You want to be able to efficiently handle a pair of guns positioned near the midpoint of each side, launching streams of gliders which almost touch before disappearing off the edge of the world. Initially you will have a pair of active zones plus the majority of the grid being dead zones. At some point the initial gliders will nearly meet and you'll want the ability to merge two active zones into one.


early pruning

Your .reverse() exploration produces many dead ends. Some of that is isolated "sparks" which expire in the subsequent generation. But some of that is blocks, beehives, and other still lifes which usually won't get you closer to your goal.

Recognizing such dead ends sooner than later will reduce your β branching factor, to enable faster solutions.


solve subproblems

Here is one more motivating example. Suppose you're searching for the 8 gliders that combine to form a glider gun. You know the answer already so you can place them and simulate forward.

Now place them and erase one of the gliders. Ask your code to search for the solution. Place the gliders further away and solve that related problem. Verify the running time scales sensibly -- this will require keeping track of solved / unsolved regions.


division

The .index() method performs more / division and % modulo operations than I would expect to see in code that's trying to be fast. With the python bytecode interpreter, maybe it's lost in the noise, and will only become apparent when you translate to Rust. Nonetheless, consider keeping (x, y) separate so you don't need to recover them with mod. (In general the row length won't be a nice power of two.)

Also, consider using if or max() for bounds checks, rather than supporting a toroidal wrap-around world. One trick is to create a dead region by unconditionally writing zeros around the edges on each tick.


JIT

You may find that a numba @jit decorator can speed things up, at least in a simple helper method.

Start small. Use type hints. Benchmark as you go.


style

        rows, cols, size = self.rows, self.cols, self.size
        ret = CellularAutomaton(rows, cols)

        """Hypothesis: a stack to keep track of which cell is using which
        configuration. """

This is a peculiar way to express your ideas. Embrace self already -- this is python! Rather than pushing a str constant onto the bytecode stack and then discarding it, or generating NOPs, prefer to write a comment:

        ret = CellularAutomaton(self.rows, self.cols)

        # Hypothesis: a stack to keep track of which cell is using which
        # configuration.

This code achieves its design goals, at least for small problems.

With some simple style changes I would be willing to delegate or accept maintenance tasks on this codebase.

\$\endgroup\$
1
  • \$\begingroup\$ I appreciate this feedback. If I had more time for this project and could spend more time researching the ins-and-outs of Game of Life, these optimizations would be very helpful. But unfortunately I'm on a tighter schedule and can't pursue these, so I've accepted the other answer. Thank you for your time nonetheless. \$\endgroup\$
    – SirPython
    Commented Apr 24, 2023 at 20:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.