# Recursive algorithm to obtain grid points inside a n-cube volume or surface

I wrote two functions to return point coordinates spanning the n-dimensional hypercubic lattice of half size K (meaning that the coordinates along each axis go from -K,.. 0.. K). The n-dimensional cube volume should have (2*K+1)**N grid points.

The functionality is split across two functions: one for returning the grid points inside the n-dimensional cube volume get_k_full_cube_of_dim(K, N) and the other is for returning the grid points in the surface get_k_cube_shell_of_dim(K, N). The two functions are mutually recursive, that is, each one calls each other.

The volume function is implemented by taking the (disjoint) union of all the surfaces of inner cubes of same dimension N but smaller K, including the origin grid point. The surface function is implemented by taking the volumes of cubes of one lower dimensionality (for example, for the surface of a 3-cube, I would get the volume of squares, which are 2-cubes) and use them to build the faces along each axis, filtering cells that are already included in previous iterations.

Right now the code seems to be working correctly, but I wanted to elicit some comments on readability and general implementation.

Here is the listing:

import math
import copy as cp

enable_caching = True
#cache of results
full_K_cube_of_dim_N = {}

def get_k_full_cube_of_dim(K, N):
if enable_caching:
if N in full_K_cube_of_dim_N:
K_cubes = full_K_cube_of_dim_N[N]
if K in K_cubes:
return cp.deepcopy(K_cubes[K])
if N == 1:
prev_K_cubes = full_K_cube_of_dim_N.get(N, {})
v = [[i-K] for i in range(2*K+1)]
prev_K_cubes[K] = v
full_K_cube_of_dim_N[N] = prev_K_cubes #not sure if needed, just in case
return list(v)
cube = []
cube.append([0 for i in range(N)])
for k in range(1, K+1):
subshell = get_k_cube_shell_of_dim(k,N)
print "get_k_full_cube_of_dim: retrieve shell of (k={0}, N={1}) = {2}".format(k, N, str(subshell))
cube.extend(subshell)
prev_K_cubes = full_K_cube_of_dim_N.get(N, {})
prev_K_cubes[K] = cube
full_K_cube_of_dim_N[N] = prev_K_cubes
return cp.deepcopy(cube)

#cache of results
K_cube_shell_of_dim_N = {}

def get_k_cube_shell_of_dim(K, N):
if enable_caching:
if N in K_cube_shell_of_dim_N:
K_cshell = K_cube_shell_of_dim_N[N]
if K in K_cshell:
return cp.deepcopy(K_cshell[K])
if N == 0:
return []
if N == 1:
prev_K_cshell = K_cube_shell_of_dim_N.get(N, {})
v = [[-K], [K]]
prev_K_cshell[K] = v
K_cube_shell_of_dim_N[N] = prev_K_cshell #not sure if needed, just in case
return list(v)
values = [-K,K]
shell = []
smaller_cube_cells = get_k_full_cube_of_dim(K,N-1)
if len(smaller_cube_cells) != (2*K+1)**(N-1):
ValueError( "smaller cube cells has {0} items, expected {1}".format(len(smaller_cube_cells), (2*K+1)**(N-1)) )
print "smaller cubes for building shell: " + str(smaller_cube_cells)
for axis in range(N):
empty_entry = [0 for i in range(N)]
filtered_cube_cells = []
for cube_cell in smaller_cube_cells:
entry = list(empty_entry)
iter_cube = 0
skip_cell = False
for iter in range(N):
#the axis coordinate will be set at the end
if iter == axis:
continue
#we filter the edge cells overlapping with the planes of the shell added on previous axes
print "axis: {0} iter: {1} iter_cube: {2} cube_cell: {3}".format(axis,iter,iter_cube,str(cube_cell))
cube_cell_coordinate = cube_cell[ iter_cube ]
if iter < axis and math.fabs(cube_cell_coordinate) == K:
skip_cell = True
print "skipped cell: " + str(cube_cell)
break
entry[ iter ] = cube_cell_coordinate
iter_cube += 1

if skip_cell:
continue
filtered_cube_cells.append(cube_cell)

print "cell not skipped: "+ str(cube_cell) + " update cube cells in subshell: " + str(filtered_cube_cells)
for signs in values:
finished_entry = list(entry)
finished_entry[axis] = signs
shell.append(finished_entry)
print "top shell updated: "+ str(shell)

smaller_cube_cells = filtered_cube_cells

prev_K_cshell = K_cube_shell_of_dim_N.get(N, {})
prev_K_cshell[K] = shell
K_cube_shell_of_dim_N[N] = prev_K_cshell
return cp.deepcopy(shell)


Some test results:

>>> get_k_cube_shell_of_dim(1,1)
[[-1], [1]]
>>> get_k_cube_shell_of_dim(1,2)
[[-1, -1], [1, -1], [-1, 0], [1, 0], [-1, 1], [1, 1], [0, -1], [0, 1]]
>>> get_k_full_cube_of_dim(1,3)
[[0, 0, 0], [-1, 0, 0], [1, 0, 0], [-1, -1, -1], [1, -1, -1], [-1, 1, -1], [1, 1, -1], [-1, -1, 0], [1, -1, 0], [-1, 1, 0], [1, 1, 0], [-1, -1, 1], [1, -1, 1], [-1, 1, 1], [1, 1, 1], [-1, 0, -1], [1, 0, -1], [-1, 0, 1], [1, 0, 1], [0, -1, 0], [0, 1, 0], [0, -1, -1], [0, 1, -1], [0, -1, 1], [0, 1, 1], [0, 0, -1], [0, 0, 1]]


• I find your code hard to understand. You should write docstrings and more comments.
• The function get_k_cube_shell_of_dim is long and complicated. We can't even view it all at once on this page. I did not go through all of it. Break it into smaller functions.
• The caching you are doing is commonly called memoization, and can be conveniently implemented via a decorator. That would simplify the function itself. See for example http://wiki.python.org/moin/PythonDecoratorLibrary#Memoize

• Disregarding the previous point, full_K_cube_of_dim_N would be better as an instance of collections.defaultdict(dict).

That would allow you to replace

prev_K_cubes = full_K_cube_of_dim_N.get(N, {})
prev_K_cubes[K] = cube
full_K_cube_of_dim_N[N] = prev_K_cubes


with

full_K_cube_of_dim_N[N][K] = cube

• [0 for i in range(N)] is equivalent to N*[0]

• get_k_cube_shell_of_dim does not trust get_k_full_cube_of_dim to return a valid result (there's a len check to raise a ValueError). An assertion before the return statement in get_k_cube_shell_of_dim would be more appropriate.