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

Question

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 = []
    #add the origin cell
    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)

        #discard cells with edge coordinates that are already in the shell
        smaller_cube_cells = filtered_cube_cells

    #add result to cache
    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]]

Answer 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.