def neighbours_find(a_dimension, a_depth, *args, **kwargs):
if (a_depth < <NUM_LIT:1>):<EOL><INDENT>l_D = []<EOL>l_I = {}<EOL>return l_I, l_D<EOL><DEDENT>b_returnUnion = False <EOL>b_wrapGridEdges = False <EOL>b_gridSize = False <EOL>for key, value in kwargs.items():<EOL><INDENT>if key == '<STR_LIT>':        b_returnUnion = value<EOL>if key == '<STR_LIT>':      b_wrapGridEdges = value<EOL>if key == '<STR_LIT>':<EOL><INDENT>a_gridSize = value<EOL>b_gridSize = True<EOL>if type(a_gridSize).__name__ != "<STR_LIT>":<EOL><INDENT>error_exit("<STR_LIT>",<EOL>"<STR_LIT>",<EOL>"<STR_LIT>", <NUM_LIT:2>)<EOL><DEDENT><DEDENT><DEDENT>if b_wrapGridEdges and not b_gridSize:<EOL><INDENT>error_exit("<STR_LIT>",<EOL>"<STR_LIT>",<EOL>"<STR_LIT>",<EOL><NUM_LIT:1>)<EOL><DEDENT>v_origin = zeros((a_dimension));<EOL>b_origin = <NUM_LIT:0>;<EOL>if len(args):<EOL><INDENT>av_origin = args[<NUM_LIT:0>];<EOL>if av_origin.size == a_dimension:<EOL><INDENT>v_origin = av_origin;<EOL>b_origin = <NUM_LIT:1>;<EOL><DEDENT><DEDENT>l_D = [zeros(a_dimension)] * a_depth<EOL>l_I = [zeros(a_dimension)] * a_depth<EOL>d = <NUM_LIT:1>;<EOL>hypercube = pow((<NUM_LIT:2> * d + <NUM_LIT:1>), a_dimension);<EOL>hypercubeInner = <NUM_LIT:1>;<EOL>orthogonals = <NUM_LIT:2> * a_dimension;<EOL>l_D[<NUM_LIT:0>] = zeros((orthogonals, a_dimension))<EOL>l_I[<NUM_LIT:0>] = zeros((hypercube - orthogonals - <NUM_LIT:1>, a_dimension));<EOL>for d in arange(<NUM_LIT:2>, a_depth + <NUM_LIT:1>):<EOL><INDENT>hypercubeInner = hypercube;<EOL>hypercube = pow((<NUM_LIT:2> * d + <NUM_LIT:1>), a_dimension);<EOL>l_D[d - <NUM_LIT:1>] = zeros((orthogonals, a_dimension))<EOL>l_I[d - <NUM_LIT:1>] = zeros((hypercube - orthogonals - hypercubeInner,<EOL>a_dimension))<EOL><DEDENT>M_bDoffset = ones((a_dimension)) * -a_depth;<EOL>M_current = ones((a_dimension));<EOL>M_currentAbs = ones((a_dimension));<EOL>M_ii = [<NUM_LIT:0>] * a_depth<EOL>M_dd = [<NUM_LIT:0>] * a_depth<EOL>for i in arange(<NUM_LIT:0>, hypercube):<EOL><INDENT>str_progress = '<STR_LIT>' %(i, hypercube - <NUM_LIT:1>, i / (hypercube - <NUM_LIT:1>) * <NUM_LIT:100>);<EOL>M_current = arr_base10toN(i, (<NUM_LIT:2> * a_depth + <NUM_LIT:1>), a_dimension);<EOL>M_current = M_current + M_bDoffset;<EOL>M_currentAbs = abs(M_current);<EOL>neighbour = max(M_currentAbs);<EOL>if b_origin:<EOL><INDENT>M_current += v_origin<EOL><DEDENT>if b_gridSize:<EOL><INDENT>if min(M_current) < <NUM_LIT:0>:<EOL><INDENT>if not b_wrapGridEdges: continue<EOL>else:<EOL><INDENT>Wt = where(M_current < <NUM_LIT:0>)<EOL>W = Wt[<NUM_LIT:0>]  <EOL>for el in arange(<NUM_LIT:0>, len(W)):<EOL><INDENT>M_current[W[el]] += a_gridSize[W[el]]<EOL><DEDENT><DEDENT><DEDENT>M_inGrid = a_gridSize - M_current<EOL>if min(M_inGrid) <= <NUM_LIT:0>:<EOL><INDENT>if not b_wrapGridEdges: continue<EOL>else:<EOL><INDENT>Wt = where(M_inGrid <= <NUM_LIT:0>)<EOL>W = Wt[<NUM_LIT:0>]  <EOL>for el in arange(<NUM_LIT:0>, len(W)):<EOL><INDENT>M_current[W[el]] = -M_inGrid[W[el]]<EOL><DEDENT>M_current = abs(M_current)<EOL><DEDENT><DEDENT><DEDENT>idx = int(neighbour) - <NUM_LIT:1><EOL>if(sum(M_currentAbs) > neighbour):<EOL><INDENT>l_I[idx][M_ii[idx]] = M_current<EOL>M_ii[idx] += <NUM_LIT:1>;<EOL><DEDENT>else:<EOL><INDENT>if(sum(M_currentAbs) == neighbour and neighbour):<EOL><INDENT>l_D[idx][M_dd[idx]] = M_current<EOL>M_dd[idx] += <NUM_LIT:1>;<EOL><DEDENT><DEDENT><DEDENT>if b_gridSize:<EOL><INDENT>for layer in arange(<NUM_LIT:0>, a_depth, dtype=int):<EOL><INDENT>l_I[layer] = l_I[layer][<NUM_LIT:0>:M_ii[layer], :]<EOL>l_D[layer] = l_D[layer][<NUM_LIT:0>:M_dd[layer], :]<EOL><DEDENT><DEDENT>if b_returnUnion:<EOL><INDENT>l_A = []<EOL>l_Al = v_origin<EOL>for layer in arange(<NUM_LIT:0>, a_depth, dtype=int):<EOL><INDENT>l_A.append(vstack((l_I[layer], l_D[layer])))<EOL>l_Al = vstack((l_Al, l_A[layer]))<EOL><DEDENT>return l_Al<EOL><DEDENT>return l_I, l_D<EOL>
SYNOPSIS

    [I, D] = neighbours_find(a_dimension, a_depth, *args, **kwargs)

ARGS

   INPUT
   a_dimension         int32           number of dimensions
   a_depth             int32           depth of neighbours to find

   OPTIONAL
   av_origin           array          row vector that defines the
                                      + origin in the <a_dimension>
                                      + space. Neighbours' locations
                                      + are returned relative to this
                                      + origin. Default is the zero
                                      + origin.

                                      If specified, this *must* be
                                      a 1xa_dimension nparray, i.e.
                                      np.array( [(x,y)] ) in the
                                      2D case.

   OUTPUT
   I                   list           a python list (per depth) of
                                           "indirect" neighbour
                                           information. Each index
                                           of the list is a row - order
                                           matrix of "index" distant
                                           indirect neighbours
   D                   list           a python list (per depth) of
                                           "direct" neighbour
                                           information. Each index
                                           of the list is a row - order
                                           matrix of "index" distant
                                           direct neighbours
                - OR -
   A                   list           a single array that is the union
                                           of I and D across all layers

Named keyword args

  "returnUnion"      If True, will return A = I union D, flattened
                     including the origin point. In effect, this is
                     all combinations of the coordinates along the
                     axes.
  "wrapGridEdges"    If True, wrap around the edges of grid.
                     If False, do not wrap around grid edges.
  "gridSize"         Specifies the gridSize for 'wrapEdges'.

  If "gridSize" is set, and "wrapEdges" is not, then the neighbors will
  not include coordinates "outside" of the grid.

DESC
       This method determines the neighbours of a point in an
       n - dimensional discrete space. The "depth" (or ply) to
       calculate is `a_depth'.

       Indirect neighbours are non-orthogonal.
       Direct neighbours are orthogonal.

       Note: this operation is identical to finding all the
       combinations of a given set of elements.

       The neighbours_findFast() method is several orders
       of magnitude faster, and should probably always be
       used over this method. This will remain mostly for
       historical/testing purposes. Also, this method does
       return neighbours as "indirect" and "direct" and
       per "layer" depth, which the ..._findFast does not.


 PRECONDITIONS
       o The underlying problem is discrete.

 POSTCONDITIONS
       o I cell array contains the Indirect neighbours
       o D cell array contains the Direct neighbours
       o size{I} = size{D}
       o size{I} = a_depth

 HISTORY
 20 December 2011
 o Translated from MatLAB script of same name.