Enumerate the isotropic neighbourhoods for 3D CA with Moore neighbourhood and 2D range-2 Moore neighbourhood

3DCA_iso_neighbourhoods.py

# Enumerate the isotropic neighbourhoods for 3D CA with Moore neighbourhood
# Determine equivalent neighbourhood configurations using permutations instead
# of matrix multiplication

from itertools import combinations, product
import numpy as np
from timeit import default_timer as timer

# Matrix representation of octahedral group (from Wikiversity)
# https://en.wikiversity.org/wiki/Full_octahedral_group#Code
octahedral = { (0, 0): [[ 1, 0, 0], [0, 1, 0], [0, 0, 1]], (0, 1): [[0, 1, 0], [ 1, 0, 0], [0, 0, 1]], 
               (0, 2): [[ 1, 0, 0], [0, 0, 1], [0, 1, 0]], (0, 3): [[0, 1, 0], [0, 0, 1], [ 1, 0, 0]],
               (0, 4): [[0, 0, 1], [ 1, 0, 0], [0, 1, 0]], (0, 5): [[0, 0, 1], [0, 1, 0], [ 1, 0, 0]],
               (1, 0): [[-1, 0, 0], [0, 1, 0], [0, 0, 1]], (1, 1): [[0,-1, 0], [ 1, 0, 0], [0, 0, 1]],
               (1, 2): [[-1, 0, 0], [0, 0, 1], [0, 1, 0]], (1, 3): [[0,-1, 0], [0, 0, 1], [ 1, 0, 0]],
               (1, 4): [[0, 0,-1], [ 1, 0, 0], [0, 1, 0]], (1, 5): [[0, 0,-1], [0, 1, 0], [ 1, 0, 0]],
               (2, 0): [[ 1, 0, 0], [0,-1, 0], [0, 0, 1]], (2, 1): [[0, 1, 0], [-1, 0, 0], [0, 0, 1]],
               (2, 2): [[ 1, 0, 0], [0, 0,-1], [0, 1, 0]], (2, 3): [[0, 1, 0], [0, 0,-1], [ 1, 0, 0]],
               (2, 4): [[0, 0, 1], [-1, 0, 0], [0, 1, 0]], (2, 5): [[0, 0, 1], [0,-1, 0], [ 1, 0, 0]],
               (3, 0): [[-1, 0, 0], [0,-1, 0], [0, 0, 1]], (3, 1): [[0,-1, 0], [-1, 0, 0], [0, 0, 1]],
               (3, 2): [[-1, 0, 0], [0, 0,-1], [0, 1, 0]], (3, 3): [[0,-1, 0], [0, 0,-1], [ 1, 0, 0]],
               (3, 4): [[0, 0,-1], [-1, 0, 0], [0, 1, 0]], (3, 5): [[0, 0,-1], [0,-1, 0], [ 1, 0, 0]],
               (4, 0): [[ 1, 0, 0], [0, 1, 0], [0, 0,-1]], (4, 1): [[0, 1, 0], [ 1, 0, 0], [0, 0,-1]],
               (4, 2): [[ 1, 0, 0], [0, 0, 1], [0,-1, 0]], (4, 3): [[0, 1, 0], [0, 0, 1], [-1, 0, 0]],
               (4, 4): [[0, 0, 1], [ 1, 0, 0], [0,-1, 0]], (4, 5): [[0, 0, 1], [0, 1, 0], [-1, 0, 0]],
               (5, 0): [[-1, 0, 0], [0, 1, 0], [0, 0,-1]], (5, 1): [[0,-1, 0], [ 1, 0, 0], [0, 0,-1]],
               (5, 2): [[-1, 0, 0], [0, 0, 1], [0,-1, 0]], (5, 3): [[0,-1, 0], [0, 0, 1], [-1, 0, 0]],
               (5, 4): [[0, 0,-1], [ 1, 0, 0], [0,-1, 0]], (5, 5): [[0, 0,-1], [0, 1, 0], [-1, 0, 0]],
               (6, 0): [[ 1, 0, 0], [0,-1, 0], [0, 0,-1]], (6, 1): [[0, 1, 0], [-1, 0, 0], [0, 0,-1]],
               (6, 2): [[ 1, 0, 0], [0, 0,-1], [0,-1, 0]], (6, 3): [[0, 1, 0], [0, 0,-1], [-1, 0, 0]],
               (6, 4): [[0, 0, 1], [-1, 0, 0], [0,-1, 0]], (6, 5): [[0, 0, 1], [0,-1, 0], [-1, 0, 0]],
               (7, 0): [[-1, 0, 0], [0,-1, 0], [0, 0,-1]], (7, 1): [[0,-1, 0], [-1, 0, 0], [0, 0,-1]],
               (7, 2): [[-1, 0, 0], [0, 0,-1], [0,-1, 0]], (7, 3): [[0,-1, 0], [0, 0,-1], [-1, 0, 0]],
               (7, 4): [[0, 0,-1], [-1, 0, 0], [0,-1, 0]], (7, 5): [[0, 0,-1], [0,-1, 0], [-1, 0, 0]] }
octahedral = [np.array(octahedral[k]).transpose() for k in sorted(octahedral)]

# Neighbours in 3D range-1 neighbourhood
neighbours = list(product([-1, 0, 1], repeat=3))
neighbours.remove((0, 0, 0))
nbhrIdx = list(range(0, len(neighbours)))

# Permutation representation of octahedral group
octahedral_perms = []
for trans in octahedral:
    perm = [neighbours.index(tuple(np.array(nbhr) @ trans)) for nbhr in neighbours]
    octahedral_perms.append(perm)

# Generate all neighbourhood configurations (up to population count > Num_neighbours/2)
# and count all unique configurations under action of the octehdral group
isoCounts = [0]*14
isoCounts[0] = 1
for N in range(1, 14):
    time0 = timer()
    isoNbhds = set()
    for nbhd in combinations(nbhrIdx, N):
        inSet = False
        for perm in octahedral_perms[1:]:
            tr_nbhd = frozenset([perm[nb] for nb in nbhd])
            if tr_nbhd in isoNbhds:
                inSet = True
                break
        if not inSet:
            isoNbhds.add(frozenset(nbhd))
    isoCounts[N] = len(isoNbhds)
    time1 = timer()
    print(int(time1-time0), len(isoNbhds))
    time0 = time1

print(isoCounts)
print(sum(isoCounts + isoCounts[:-1]))

R2_iso_neighbourhoods.py

from itertools import combinations, product
import numpy as np
from timeit import default_timer as timer

# Matrix representation of dihedral_4 group (cubic symmetries)
dihedral4 = { 'r0': [[1, 0], [0, 1]], 'r1':[[0, -1], [1, 0]],
              'r2': [[-1, 0], [0, -1]], 'r3': [[0, 1], [-1, 0]],
              's0': [[1, 0], [0, -1]], 's1':[[0, 1], [1, 0]],
              's2': [[-1, 0], [0, 1]], 's3': [[0, -1], [-1, 0]] }

dihedral4 = [np.array(dihedral4[k]).transpose() for k in  sorted(dihedral4)]

# Neighbours in 2D range-2 neighbourhood
neighbours = list(product([-2, -1, 0, 1, 2], repeat=2))
neighbours.remove((0, 0))

# Generate and count all the combinations of N neighbours under action of the dihedral_4 group:
isoCounts = [0]*13
isoCounts[0] = 1
for N in range(1, 13):
    time0 = timer()
    isoNbhds = set()
    for neighbourhood in combinations(neighbours, N):
        points = np.array(neighbourhood)
        inSet = False
        for trans in dihedral4[1:]:
            trPoints = frozenset((tuple(row) for row in points @ trans))
            if trPoints in isoNbhds:
                inSet = True
                break
        if not inSet:
            fsPoints = frozenset((tuple(row) for row in points))
            isoNbhds.add(fsPoints)
    isoCounts[N] = len(isoNbhds)
    time1 = timer()
    print(int(time1-time0), len(isoNbhds))
    time0 = time1

print(isoCounts)
print(sum(isoCounts + isoCounts[:-1]))