Mainly Tech projects on Python and Electronic Design Automation.

Sunday, January 31, 2021

Counting types of neighbour - in Python.

(I tried using the Spyder IDE to write a more literal Python program. Output samples needed to be pasted in to the source, but I didn't particularly want many of Jupyters strengths, so tried Spyder together with http://hilite.me/)

Best viewed on larger than 7" screens.


# -*- coding: utf-8 -*-
"""
Created on Sat Jan 30 10:24:54 2021

@author: Paddy3118
"""

#from conway_cubes_2_plus_dims import *
#from conway_cubes_2_plus_dims import neighbours
from itertools import product


# %% #Full neighbourhood
"""
In my [previous post][1] I showed how the number of neighbours to a point in
N-Dimensional space increases exponentially:

[1]: https://paddy3118.blogspot.com/2021/01/conways-game-of-life-in-7-dimensions.html

"""

# %% ##Function

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

for dim in range(1,10):
    n = neighbours(tuple([0] * dim))
    print(f"A {dim}-D grid point has {len(n):_} neighbours.")
    assert len(n) == 3**dim - 1

# %% ##Output
"""
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.
"""

# %% ##Visuals
"""
Here we define neighbours as all cells, n whose indices differ by at most 1,
i.e. by -1, 0, or +1 from the point X itself; apart from the point X.

In 1-D:

    .nXn.

In 2-D:

    .....
    .nnn.
    .nXn.
    .nnn.
    .....

# 'Axial' neighbours

 There is another definition of neighbour that "cuts out the diagonals" in 2-D
 to form:

     .....
     ..n..
     .nXn.
     ..n..
     .....

In 3-D this would add two extra cells, one directly above and one below X.

In 1-D this would be the two cells one either side of X.

I call this axial because if X is translated to the origin then axial
neighbours are the 2N cells in which only one coordinate can differ by either
+/- 1.

# Let's calclate them.
"""

# %% #Axial neighbourhood
from typing import Tuple, Set, Dict

Point = Tuple[int, int, int]
PointSet = Set[Point]


def axial_neighbours(point: Point) -> PointSet:
    dim, pt = len(point), list(point)
    return {tuple(pt[:d] + [pt[d] + delta] + pt[d+1:])
            for d in range(dim) for delta in (-1, +1)}

print("\n" + '='*60 + '\n')
for dim in range(1,10):
    n = axial_neighbours(tuple([0] * dim))
    text = ':\n  ' + str(sorted(n)) if dim <= 2 else ''
    print(f"A {dim}-D grid point has {len(n):_} axial neighbours{text}")
    assert len(n) == 2*dim

# %% ##Output
"""
============================================================

A 1-D grid point has 2 axial neighbours:
  [(-1,), (1,)]
A 2-D grid point has 4 axial neighbours:
  [(-1, 0), (0, -1), (0, 1), (1, 0)]
A 3-D grid point has 6 axial neighbours
A 4-D grid point has 8 axial neighbours
A 5-D grid point has 10 axial neighbours
A 6-D grid point has 12 axial neighbours
A 7-D grid point has 14 axial neighbours
A 8-D grid point has 16 axial neighbours
A 9-D grid point has 18 axial neighbours
"""

# %% #Another way of classifying...
"""
The 3-D axial neighbourhood around X can be described as:
    Think of X as a cube with six sides oriented with the axes of the 3
    dimensions. The axial neighbourhood is six similar cubes with one face
    aligned with each of X's faces. A kind of 3-D cross of cubes.

In 3-D, the "full" neighbourhood around point X describes a 3x3x3 units cube
centered on X.

In 3-D: Thinking of that 3x3x3 cube around X:

    What if we excluded the six corners of that cube?
    * Those excluded corner points have all their three coordinates different
      from those of X,
      i.e. if excluded e = Point(e_x, e_y, e_z) and X = Point(x_x, x_y, x_z):
          e_x != x_x and e_y != x_y and e_z != x_z
      Also:
          |e_x - x_x| == 1 and |e_y - x_y| == 1 and |e_z - x_z| == 1

    * We _Keep_ all points in the cube that have 1 and 2 different coords to
      those of X

    Another definition of the _axial_ neighbourhood case is
    * We _Keep_ all points in the cube that have only 1 different coord to
      those of X

This can be generalised to thinking in N dimensions of neighbourhood types
keeping from 1 to _up to_ N differences in coordinates (DinC).
"""
def axial_neighbours(point: Point) -> PointSet:
    dim, pt = len(point), list(point)
    return {tuple(pt[:d] + [pt[d] + delta] + pt[d+1:])
            for d in range(dim) for delta in (-1, +1)}

#%% #Neighbourhood by Differences in Coordinates

from collections import defaultdict
from pprint import pformat
from math import comb, factorial as fact


def d_in_c_neighbourhood(dim: int) -> Dict[int, PointSet]:
    """
    Split neighbourhood cube around origin point in dim-dimensions mapping
    count of coords not-equal to origin -> those neighbours
    """
    origin = tuple([0] * dim)
    cube = {n
            for n in product(*((n-1, n, n+1)
                               for n in origin))
            if n != origin}
    d_in_c = defaultdict(set)
    for n in cube:
        d_in_c[dim - n.count(0)].add(n)
    return dict(d_in_c)

def _check_counts(d: int, c: int, n: int) -> None:
    """
    Some checks on counts

    Parameters
    ----------
    d : int
        Dimensionality.
    c : int
        Count of neighbour coords UNequal to origin.
    n : int
        Number of neighbour points with exactly c off-origin coords.

    Returns
    -------
    None.

    """
    # # From OEIS
    # if   c == 1: assert n == d* 2
    # elif c == 2: assert n == d*(d-1)* 2
    # elif c == 3: assert n == d*(d-1)*(d-2)* 4 / 3
    # elif c == 4: assert n == d*(d-1)*(d-2)*(d-3)* 2 / 3
    # elif c == 5: assert n == d*(d-1)*(d-2)*(d-3)*(d-4)* 4 / 15
    # # Getting the hang of it
    # elif c == 6: assert n == d*(d-1)*(d-2)*(d-3)*(d-4)*(d-5)* 4 / 45

    # Noticed some powers of 2 leading to
    # if c == 6: assert n == d*(d-1)*(d-2)*(d-3)*(d-4)*(d-5)* 2**6 / fact(6)
    # if c == 6: assert n == fact(d) / fact(c) / fact(d - c) * 2**c

    # Finally:
    assert n == fact(d) / fact(c) / fact(d - c) * 2**c


print("\n" + '='*60 + '\n')
for dim in range(1,10):
    d = d_in_c_neighbourhood(dim)
    tot = sum(len(n_set) for n_set in d.values())
    print(f"A {dim}-D point has a total of {tot:_} full neighbours of which:")
    for diff_count in sorted(d):
        n_count = len(d[diff_count])
        print(f"  {n_count:4_} have exactly {diff_count:2} different coords.")
        _check_counts(dim, diff_count, n_count)

# %% ##Output
"""
============================================================

A 1-D point has a total of 2 full neighbours of which:
     2 have exactly  1 different coords.
A 2-D point has a total of 8 full neighbours of which:
     4 have exactly  1 different coords.
     4 have exactly  2 different coords.
A 3-D point has a total of 26 full neighbours of which:
     6 have exactly  1 different coords.
    12 have exactly  2 different coords.
     8 have exactly  3 different coords.
A 4-D point has a total of 80 full neighbours of which:
     8 have exactly  1 different coords.
    24 have exactly  2 different coords.
    32 have exactly  3 different coords.
    16 have exactly  4 different coords.
A 5-D point has a total of 242 full neighbours of which:
    10 have exactly  1 different coords.
    40 have exactly  2 different coords.
    80 have exactly  3 different coords.
    80 have exactly  4 different coords.
    32 have exactly  5 different coords.
A 6-D point has a total of 728 full neighbours of which:
    12 have exactly  1 different coords.
    60 have exactly  2 different coords.
   160 have exactly  3 different coords.
   240 have exactly  4 different coords.
   192 have exactly  5 different coords.
    64 have exactly  6 different coords.
A 7-D point has a total of 2_186 full neighbours of which:
    14 have exactly  1 different coords.
    84 have exactly  2 different coords.
   280 have exactly  3 different coords.
   560 have exactly  4 different coords.
   672 have exactly  5 different coords.
   448 have exactly  6 different coords.
   128 have exactly  7 different coords.
A 8-D point has a total of 6_560 full neighbours of which:
    16 have exactly  1 different coords.
   112 have exactly  2 different coords.
   448 have exactly  3 different coords.
  1_120 have exactly  4 different coords.
  1_792 have exactly  5 different coords.
  1_792 have exactly  6 different coords.
  1_024 have exactly  7 different coords.
   256 have exactly  8 different coords.
A 9-D point has a total of 19_682 full neighbours of which:
    18 have exactly  1 different coords.
   144 have exactly  2 different coords.
   672 have exactly  3 different coords.
  2_016 have exactly  4 different coords.
  4_032 have exactly  5 different coords.
  5_376 have exactly  6 different coords.
  4_608 have exactly  7 different coords.
  2_304 have exactly  8 different coords.
   512 have exactly  9 different coords.
"""

# %% #Manhattan
"""
What, up until now, I have called the Differences in Coordinates is better
known as the Manhattan distance, or [Taxicab geometry][2].


[2]: https://en.wikipedia.org/wiki/Taxicab_geometry
"""

# %% #Formula for Taxicab counts of cell neighbours
"""
Function _check_counts was originally set up as a crude check of taxicab
distance of 1 as the formula was obvious.

The formulea for taxcab distances 2, 3, and 4 were got from a search on [OEIS][3]

The need for factorials is obvious. The count for a taxicab distance of N in
N-D space is 2**N which allowed me to work out a final formula:

If:
    d is the number of spatial dimiensions.
    t is the taxicab distance of a neighbouring point to the origin.
Then:
    n, the count of neighbours with exactly this taxicab distance is

        n = f(d, t)
          = d! / t! / (d-t)! * 2**t


We can use the assertion from the exercising of function neighbours() at the
beginning to state that:

        sum(f(d, t)) for t from 0 .. d
          = g(d)
          = 3**d - 1


[3]: https://oeis.org/
"""

# %% ##Taxicab visualisation
"""
Lets create a visualisation of the difference in coords for neighbours in <= 4-D.
(The function is general, but I'm lost after 3-D)!

The origin will show as zero, 0; and neighbours surround it as digits which are
the taxicab distance from 0.
"""

def to_str(taxi: Dict[int, PointSet], indent: int=4) -> str:
    """


    Parameters
    ----------
    taxi : Dict[int, PointSet]
        Map taxicab distance from origin
        -> set off neighbours with that difference.
    indent : int
        indent output with spaces

    Returns
    -------
    str
        Visusalisation of region.

    """

    if not taxi:
        return ''

    ap = set()              # all points
    neighbour2taxi = {}
    for taxi_distance, nbrs in taxi.items():
        ap |= nbrs
        for n in nbrs:
            neighbour2taxi[n] = taxi_distance

    # Dimensionality
    dims = len(n)

    # Add in origin showing as 0
    origin = tuple([0] * dims)
    ap.add(origin)
    neighbour2taxi[origin] = 0

    # Plots neighbourhood of origin (plus extra space)
    space = 1
    minmax = [range(-(1 + space), (1 + space) + 1)
              for dim in range(dims)]

    txt = []
    indent_txt = ' ' * indent
    for plane_coords in product(*minmax[2:]):
        ptxt = ['\n' + indent_txt
                + ', '.join(f"dim{dim}={val}"
                            for dim, val in enumerate(plane_coords, 2))]

        ptxt += [''.join((str(neighbour2taxi[tuple([x, y] + list(plane_coords))])
                          if tuple([x, y] + list(plane_coords)) in ap
                          else '.')
                         for x in minmax[0])
                 for y in minmax[1]]

        # Don't plot planes with no neighbours (due to extra space).
        if ''.join(ptxt).count('.') < (3 + space*2)**2:
            txt += ptxt

    return '\n'.join(indent_txt + t for t in txt)


print("\n" + '='*60 + '\n')
for dim in range(2,5):
    d = d_in_c_neighbourhood(dim)
    tot = sum(len(n_set) for n_set in d.values())
    print(f"\nA {dim}-D point has a total of {tot:_} full neighbours of which:")
    for diff_count in sorted(d):
        n_count = len(d[diff_count])
        print(f"  {n_count:4_} have taxicab distance {diff_count:2} from the origin.")
        _check_counts(dim, diff_count, n_count)
    print(to_str(d))

# %% ##Output
"""
============================================================


A 2-D point has a total of 8 full neighbours of which:
     4 have taxicab distance  1 from the origin.
     4 have taxicab distance  2 from the origin.


    .....
    .212.
    .101.
    .212.
    .....

A 3-D point has a total of 26 full neighbours of which:
     6 have taxicab distance  1 from the origin.
    12 have taxicab distance  2 from the origin.
     8 have taxicab distance  3 from the origin.

    dim2=-1
    .....
    .323.
    .212.
    .323.
    .....

    dim2=0
    .....
    .212.
    .101.
    .212.
    .....

    dim2=1
    .....
    .323.
    .212.
    .323.
    .....

A 4-D point has a total of 80 full neighbours of which:
     8 have taxicab distance  1 from the origin.
    24 have taxicab distance  2 from the origin.
    32 have taxicab distance  3 from the origin.
    16 have taxicab distance  4 from the origin.

    dim2=-1, dim3=-1
    .....
    .434.
    .323.
    .434.
    .....

    dim2=-1, dim3=0
    .....
    .323.
    .212.
    .323.
    .....

    dim2=-1, dim3=1
    .....
    .434.
    .323.
    .434.
    .....

    dim2=0, dim3=-1
    .....
    .323.
    .212.
    .323.
    .....

    dim2=0, dim3=0
    .....
    .212.
    .101.
    .212.
    .....

    dim2=0, dim3=1
    .....
    .323.
    .212.
    .323.
    .....

    dim2=1, dim3=-1
    .....
    .434.
    .323.
    .434.
    .....

    dim2=1, dim3=0
    .....
    .323.
    .212.
    .323.
    .....

    dim2=1, dim3=1
    .....
    .434.
    .323.
    .434.
    .....

"""

No comments:

Post a Comment

Followers

Subscribe Now: google

Add to Google Reader or Homepage

Go deh too!

whos.amung.us

Blog Archive