Mainly Tech projects on Python and Electronic Design Automation.

Friday, January 29, 2021

Conways Game of Life in 7 Dimensions!

 

I stumbled across the site "Advent of Code" and its puzzle Conways cube by accident. I really liked the idea of the problem, even though it took me some time to understand the given examples, but then I really got into it.

I wrote a 3D answer; updated it to a 4D answer, then set about changing the Python source to handle any dimension >= 2.

The problem 

This is the gist of what is requested

  • Conways Game of life, (GOL), is usually played on a 2-D rectangular grid of cells that are in one of two states: inactive or active at the end of any cycle
  • We map a coordinate system on the grid where grid intersections are atall whole  integer coordinates. (Negative too).
  • The neighbours of a cell on the grid are defined as any cell with all coordinates that differs from the target cell by either -1, 0, or +1 independently - apart from the target cell in the center of course.
    For the original 2-D GOL, each cell has 8 neighbours.
  • The grid, (called a  Pocket Domain in the puzzle), is first seeded/initialised with an arrangement of activated  cells in the 2-D plane coresponding to cycle(0).
  • The configuration of active cells in the next cycle cycle(n+1) is calculated from the configuration in the current cycle cycle(n) by applying the following two simple rules to all cells:
    1. If a cell is active in cycle(n) it will stay active if it has 2 or 3 active neighbours; otherwise it becomes inactive in cycle(n+1).
    2. If a cell is inactive in cycle(n) it will become active if it has exactly 3 active neighbours; otherwise it becomes inactive in cycle(n+1).
  •  Note that although the pocket domain/grid is infinite, the rules need only be applied in the neighbourhood of any active cedlls as all cells outside that are inactive and can only remain inactive in the next cycle.
  • The puzzle defines a starting configuration/format and asks how many cells are active in some future cycle.

The textual format

For 2-D it is a rectangular plane of full stops '.' and hashes '#' where hashes represent active cells on the plane, and full-stops inactive cells. Lets call this the x-y coordinate plane

  • The rectangle shown is only large enough to show all active cells.
 
.#.
..#
###

A 3-D pocket domain is visualised as one or more x-y slices through the cubic x,y,z extent of the domain,The z coordinate for each slice in the cycle is different and stated above each slice like so:

z=-1
#..
..#
.#.

z=0
#.#
.##
.#.

z=1
#..
..#
.#.

The range of x and y coords shown, of all x-y slices for a cycle stay the same; but may be different between cycles as the extent of all active cells changes between cycles.

Their 4-D textual format introduces w as the forth dimension and now x-y planes through the encompassing hyper-rectangle of all active cells need both the z and w coords specified for each plane.

z=-1, w=-1
#..
..#
.#.

z=0, w=-1
#..
..#
.#.

z=1, w=-1
#..
..#
.#.

z=-1, w=0
#..
..#
.#.

N-Dimensional Solution

My N-Dimensional extension to the Text Format

My format is a slight change in that it will read the above, but names the dimensions greater than two, (assuming the x-y plane  as dim0 and dim1), dim2, dim3, ... for example:
dim2=0, dim3=0, dim4=0
.#.
..#
###

The internal grid data-structure

I started with no idea if there was to be a limit to the extent of the 3-D pocket dimension, and new that 2-D GOL looked pretty sparse w.r.t. active cells.
 
So I opted for a sparse grid.

A bit of thought and I confirmed I only needed to keep active cells, from which I could calculate all neighbours that could be affected then calculate what the active cells of the next cycle become.

If I keep the coordinate of cells as a tuple then I will be doing a lot of lookup on whether a coord is active or not. I made the pocket domain/grid/pd a set of active points where points were tuples. Checking if a neighbour is active amounts to checking if it is in the set.

Reading the textual format

My initial 3-D reader only reads the 3-d format. It uses a regexp to extract the Z and XY sections and turns them into a set of the active points defined.

from itertools import product
import re

#%% to/From textual format

# Help: https://regex101.com/codegen?language=python
pocket_plane_re = r"""
	# NOTE: add an extra \n\n to end of string being matched

	^z=(?P<Z>[-+0-9]+?)    # Single z-plane coord
	\s*?\n
	(?P<XY>.*?\n)\n        # Following xy-plane description
	"""

def _plane_from_str(text, z):
    """single pocket dim z-plane to set of active cell coords"""
    pd = {(x, y, z)         # Pocket dimension
          for y, row in enumerate(text.strip().split('\n'))
          for x, state in enumerate(row)
          if state == '#'}
    return pd

def from_str(text):
    """Pocket dim text to set of active cell coords, the internal format"""
    pd = set()              # pocket dim active cells accumulator
    matches = re.finditer(pocket_plane_re, text + '\n\n',
                          re.MULTILINE | re.DOTALL | re.VERBOSE)
    for match in matches:
        z = int(match.groupdict()['Z'])
        xy =    match.groupdict()['XY']
        pd |= _plane_from_str(xy, z)
    return pd

function from_str separates and extracts the z coord and the xy plane from the text then passes them to function _plane_from_str to return a set of active cell coordinates to add to the pocket domain.

The N-Dimensional reader 


from itertools import product
import re

#%% to/From textual format

# Help: https://regex101.com/codegen?language=python

# Parses textaul version of state.
pocket_plane_re = r"""
	# Note: add an extra \n\n to end of string being matched,
    #       and a \n to the beginning.

	# Coordinates of plane, (optional for 2D).
	^(?P<PLANE>(?: [a-z][a-z0-9]*=(?P<NUM1>-?\d+))
	    (?:,\s+ (?: [a-z][a-z0-9]*=(?P<NUM2>-?\d+)))*)?

	\s*?\n
	(?P<XY>[.#\n]+?\n)\n        # xy-plane description
	"""

def _plane_from_str(text, plane_coords):
    """single pocket dim z-plane to set of active cell coords"""
    pcoords = list(plane_coords)
    pd = {tuple([x, y] + pcoords)         # Pocket dimension
          for y, row in enumerate(text.strip().split('\n'))
          for x, state in enumerate(row)
          if state == '#'}
    return pd

def from_str(text: str) -> set:
    """Pocket dim text to set of active cell coords, the internal format"""
    pd = set()              # pocket dim active cells accumulator
    matches = re.finditer(pocket_plane_re, '\n' + text + '\n\n',
                          re.MULTILINE | re.DOTALL | re.VERBOSE)
    for match in matches:
        # Extra coords for the xy plane
        plane_coords_txt = match.groupdict()['PLANE'] or ''
        plane_coords = [int(matchp.group()[1:])
                        for matchp in
                        re.finditer(r"=-?\d+", plane_coords_txt)]
        # xy plane
        xy = match.groupdict()['XY']
        # Accumulate active cells in plane
        pd |= _plane_from_str(xy, plane_coords)
    return pd
 
An extended regexp handles their being no extra coord needed for the x-y plane, when working in 2-D, as well as the unbounded list of extra dimension coords in more than 2-D.
 

Writing the textual format.

Creates a textual representation of the internal format of a cycle.
 
In the 3-D case I first finds the minimum/maximum extents of activity on all three axis in minmax. Stepping through the range of the z axis it then finds and prints all the activity in the x-y plane for that value of z

I had had difficulty in understanding that the extend of x and y changes between cycles so argument print_extent was used in debugging.

def to_str(pd, print_extent=False):
    "From set of active cell coords to output textual format"
    # Extent of pocket dimenson
    minmax = [range(min(pd, key=lambda pt:pt[dim])[dim],
                    max(pd, key=lambda pt:pt[dim])[dim] + 1)
              for dim in range(3)]

    txt = [f"\n// FROM: {tuple(r.start for r in minmax)}"
           f" TO: {tuple(r.stop -1 for r in minmax)}"] if print_extent else []

    append = txt.append
    for z in minmax[2]:
        append(f"\nz={z}")
        for y in minmax[1]:
            append(''.join(('#' if (x, y, z) in pd else '.')
                           for x in minmax[0]))

    return '\n'.join(txt)

The N-Dimensional writer

def to_str(pd, print_extent=False) -> str:
    "From set of active cell coords to output textual format"
    if not pd:
        return ''
    # Dimensionality of pocket dimension
    point = pd.pop()
    dims = len(point)
    pd |= {point}           # Put that point back

    # Extent of pocket dimenson
    minmax = [range(min(pd, key=lambda pt:pt[dim])[dim],
                    max(pd, key=lambda pt:pt[dim])[dim] + 1)
              for dim in range(dims)]

    txt = [f"\n// FROM: {tuple(r.start for r in minmax)}"
           f" TO: {tuple(r.stop -1 for r in minmax)}"] if print_extent else []

    for plane_coords in product(*minmax[2:]):

        ptxt = ['\n' + ', '.join(f"dim{dim}={val}"
                                 for dim, val in enumerate(plane_coords, 2))]
        ptxt += [''.join(('#'
                          if tuple([x, y] + list(plane_coords)) in pd
                          else '.')
                         for x in minmax[0])
                 for y in minmax[1]]
        if '#' in ''.join(ptxt):
            txt += ptxt

    return '\n'.join(txt)

To find out how many dimensions there are, a point is poped from, (and later added back to), the set of active points, pd, and stored in dims. minmax and assignments are similar,

The creation of the text for each x-y plane is separated out. many of the x-y planes can be of all inactive cells - this modification will not display x-y planes if they have no set cells.

Generating all neighbours to a cell

An important part of the program that is repeated many times, as well as many times for the same cell coordinates.

3-D neighbour function exploration

Knowing that it is calculated many times for the same cell I decided to first have a cached function. I skipped using lru_cache for a simple, fast homebrew of using a dict _neighbour_cache.

_neighbour_cache = {}
def neighbours(point):
    "return the points 26 neighbour coords"
    if point not in _neighbour_cache:
        x, y, z = point
        _neighbour_cache[point] = {xyz
                                   for xyz in product((x-1, x, x+1),
                                                      (y-1, y, y+1),
                                                      (z-1, z, z+1))
                                   if xyz != point}
    return _neighbour_cache[point]

That cache was large and in higher dimensions would get huge or inefficient, so I decided to time some alternate functions for the 3-D case.

Alternate implementations and timings

def neighboursX(point):
    x, y, z = point
    return {
         (x-1, y-1, z-1), (x  , y-1, z-1), (x+1, y-1, z-1),
         (x-1, y  , z-1), (x  , y  , z-1), (x+1, y  , z-1),
         (x-1, y+1, z-1), (x  , y+1, z-1), (x+1, y+1, z-1),

         (x-1, y-1, z  ), (x  , y-1, z  ), (x+1, y-1, z  ),
         (x-1, y  , z  ),                  (x+1, y  , z  ),
         (x-1, y+1, z  ), (x  , y+1, z  ), (x+1, y+1, z  ),

         (x-1, y-1, z+1), (x  , y-1, z+1), (x+1, y-1, z+1),
         (x-1, y  , z+1), (x  , y  , z+1), (x+1, y  , z+1),
         (x-1, y+1, z+1), (x  , y+1, z+1), (x+1, y+1, z+1),
        }

def neighboursY(point):
    x, y, z = point
    return {xyz
            for xyz in product((x-1, x, x+1),
                               (y-1, y, y+1),
                               (z-1, z, z+1))
            if xyz != point}

"""

In [80]: p = (1, 1, 1)

In [81]: neighbours(p) == neighboursX(p) == neighboursY(p)
Out[81]: True

In [82]: %timeit neighbours(p)
719 ns ± 20.2 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

In [83]: %timeit neighboursX(p)
9.75 µs ± 82 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

In [84]: %timeit neighboursY(p)
11.2 µs ± 265 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

In [85]:

"""

Caching using a simple global dict in neighbours is an order of magnitude
faster than calculations, but calculations take no extra memory

neighboursY seems to hit the sweet spot of maintainability speed and memory use and was adopted for use in the 4-D version.

4-D neighbour function

def neighbours(point):
    x, y, z, w = point
    return {xyzw
            for xyzw in product((x-1, x, x+1),
                                (y-1, y, y+1),
                                (z-1, z, z+1),
                                (w-1, w, w+1))
            if xyzw != point}

It is just a 4-D extension of the 3-D neighboursY function.

The N-D neighbour function

def neighbours(point):
    return {n
            for n in product(*((n-1, n, n+1)
                               for n in point))
            if n != point}
 
Working that out was sweet :-)
 

Checking if a cell becomes active in the next cycle

 Take a point and the set of all active points in the current cycle then checks the points neighbours to see if this point will be active or not in the next cycle

def is_now_active(point, pd):
    "Is point active w.r.t. pd"
    ncount = 0              # active neighbours count
    for npoint in neighbours(point):
        ncount += npoint in pd
        if ncount >3:
            break
    if point in pd:         # Currently active
        return 2 <= ncount <= 3
    else:                   # Currently inactive
        return ncount == 3
 
The above didn't need changing from being written for the 3-D case and being reused in the 4-D and N-D versions.
 
in higher dimensions points can have many neighbours and many possible active neighbours but we only need a count of up to 4 for the following checks to work.
 

 Generating the next cycle from the current.

This again did not need changing from the 3-D case

def next_cycle(pd: set) -> set:
    "generate next cycle from current Pocket Dimension pd"
    possibles = set()           # cells that could become active
    for active in pd:
        possibles |= neighbours(active)
    possibles |= pd

    return {point for point in possibles if is_now_active(point, pd)}

it checks all current points and their neighbours to see if they will be active in the next cycle.


Examples of N-Dimensional results

if __name__ == '__main__':
    cycle0 = """
.#.
..#
###
"""

    pd2 =  from_str(cycle0)
    start_pd_for_dim = {d: {tuple(list(pt) + [0]*(d-2)) for pt in pd2}
                        for d in range(2, 8)}
    for dim in range(2, 8):
        print('='*32)
        n, pd = 0, start_pd_for_dim[dim]
        lines, active = 0, 1
        while True:
            text = to_str(pd)
            active = text.count('#')
            active_planes = ('\n' + text).count('\n\n')
            lines = text.count('\n')
            print(f"\n{dim}-D: After {n} cycles, {active:_} cells on "
                  f"{active_planes:_} active plane(s):\n"
                  f"{text if lines < 60 else ''}")
            if not (lines <2**12 and 0 < active < 2_000 and n < 10):
                break
            n, pd = n + 1, next_cycle(pd)

Just using the one 2-D initial cycle to create the same initial x-y plane in 2 to seven pocket domains in dict start_pd_for_dim.

Then for increasing pocket domain dimensions run the Conway cycles. cycle printing is stopped for large outputs. and whole dimensions are truncated when they get "large"

Output

A sample from the end of the output for the __main__ above

================================

6-D: After 0 cycles, 5 cells on 1 active plane(s):

dim2=0, dim3=0, dim4=0, dim5=0
.#.
..#
###

6-D: After 1 cycles, 245 cells on 81 active plane(s):


6-D: After 2 cycles, 464 cells on 48 active plane(s):


6-D: After 3 cycles, 15_744 cells on 1_504 active plane(s):

================================

7-D: After 0 cycles, 5 cells on 1 active plane(s):

dim2=0, dim3=0, dim4=0, dim5=0, dim6=0
.#.
..#
###

7-D: After 1 cycles, 731 cells on 243 active plane(s):


7-D: After 2 cycles, 1_152 cells on 112 active plane(s):


7-D: After 3 cycles, 106_400 cells on 10_320 active plane(s):


Extra: Just how many neighbours does a point have?

>>>  for d in range(1,15):
         n = len(neighbours(tuple([0] * d)))
         print(f"A {d}-D grid point has {n:_} neighbours.")
         assert n == 3**d - 1
A 1-D grid point has 2 neighbours.
A 2-D grid point has 8 neighbours.
A 3-D grid point has 26 neighbours.
A 4-D grid point has 80 neighbours.
A 5-D grid point has 242 neighbours.
A 6-D grid point has 728 neighbours.
A 7-D grid point has 2_186 neighbours.
A 8-D grid point has 6_560 neighbours.
A 9-D grid point has 19_682 neighbours.
A 10-D grid point has 59_048 neighbours.
A 11-D grid point has 177_146 neighbours.
A 12-D grid point has 531_440 neighbours.
A 13-D grid point has 1_594_322 neighbours.
A 14-D grid point has 4_782_968 neighbours.

In [3]: 
 
There's a lot of neighbours in higher dimensions!
 

 END.

No comments:

Post a Comment

Followers

Subscribe Now: google

Add to Google Reader or Homepage

Go deh too!

whos.amung.us

Blog Archive