3

Lets say I use an 2D/3D numpy array to model a box consisting of cells. The cells are labeled by increasing numbers starting from 0. Until here, this could be done by

box = np.arange(np.prod(n_cells))
box = box.reshape(n_cells)

where n_cells is a np.array which stores the number of cells the box should have in each dimensions, basically the dimensions of the box counted in cells.

Now, the hard part for me is to find an intelligent way to find the neighbors for each cell. My goal is to set up an array or a list that contains the neighbors to each cell - ideally lets look at a small 2D example.

0   1   2   3
4   5   6   7
8   9   10 11
12  13  14 15

Here, without periodic boundary conditions the neighbors are

0: 1, 4, 5
1: 0, 2, 4, 5, 6
2: 1, 5, 6, 7, 3
...

But I would like to have it with periodic boundary conditions, such that

0: 1, 4, 5, 3, 7, 12, 13, 15
1: 0, 2, 4, 5, 6, 12, 13, 14

This way, every element has 8 neighbors in 2D. Ideally, I would like to be able do create such a list/array for any dimensions, but I am especially interested in 2D/3D if a generalized solution is not available.

Another issue, I will have to adress is that this way, I count all the pairs twice, i.e. 0 is a neighbor of 1 and 1 is a neighbor of 0. This is something, I need to get rid of as well but it is not the main problem.

Marsl
  • 285
  • 1
  • 10
  • Like this? https://stackoverflow.com/a/34908879/8881141 Or this? https://stackoverflow.com/a/2373689/8881141 – Mr. T Dec 30 '17 at 14:26
  • Possible duplicate of [How to Find the Neighbors of a Cell in an ndarray?](https://stackoverflow.com/questions/34905274/how-to-find-the-neighbors-of-a-cell-in-an-ndarray) – sam-pyt Dec 30 '17 at 14:32
  • Thx, this goes in the direction of want I need but it does not adress periodicity – Marsl Dec 30 '17 at 14:37
  • You can adapt these methods with `x %= len(dimension)` and `if x < 0: x += len(dimension)` – Mr. T Dec 30 '17 at 15:21

1 Answers1

5

This can be done with numpy.roll which "rolls" the array along given axes, with exactly the kind of wrap-around that you want. For example, rolling by (-1, -1) shifts everything to the left and up, so the array becomes

  [[ 5,  6,  7,  4],
   [ 9, 10, 11,  8],
   [13, 14, 15, 12],
   [ 1,  2,  3,  0]]

In this way, we just found the south-east neighbor for every point. It remains to flatten this list (ravel), repeat the process for each of 9 offsets (including (0, 0) which means the number itself), and stack the results. The solution works for arbitrary dimension of array b:

dim = len(b.shape)       # number of dimensions
offsets = [0, -1, 1]     # offsets, 0 first so the original entry is first 
columns = []
for shift in itertools.product(offsets, repeat=dim):   # equivalent to dim nested loops over offsets
    columns.append(np.roll(b, shift, np.arange(dim)).ravel())
neighbors = np.stack(columns, axis=-1)

Output (the value of neighbors):

  [[ 0,  1,  3,  4,  5,  7, 12, 13, 15],
   [ 1,  2,  0,  5,  6,  4, 13, 14, 12],
   [ 2,  3,  1,  6,  7,  5, 14, 15, 13],
   [ 3,  0,  2,  7,  4,  6, 15, 12, 14],
   [ 4,  5,  7,  8,  9, 11,  0,  1,  3],
   [ 5,  6,  4,  9, 10,  8,  1,  2,  0],
   [ 6,  7,  5, 10, 11,  9,  2,  3,  1],
   [ 7,  4,  6, 11,  8, 10,  3,  0,  2],
   [ 8,  9, 11, 12, 13, 15,  4,  5,  7],
   [ 9, 10,  8, 13, 14, 12,  5,  6,  4],
   [10, 11,  9, 14, 15, 13,  6,  7,  5],
   [11,  8, 10, 15, 12, 14,  7,  4,  6],
   [12, 13, 15,  0,  1,  3,  8,  9, 11],
   [13, 14, 12,  1,  2,  0,  9, 10,  8],
   [14, 15, 13,  2,  3,  1, 10, 11,  9],
   [15, 12, 14,  3,  0,  2, 11,  8, 10]]

In each row, the first entry is the original number, the others are its neighbors.

To have each entry-neighbor pair listed only once, you can mask the redundant entries, for example with NaN:

np.where(neighbors >= neighbors[:, [0]], neighbors, np.nan)

  [[  0.,   1.,   3.,   4.,   5.,   7.,  12.,  13.,  15.],
   [  1.,   2.,  nan,   5.,   6.,   4.,  13.,  14.,  12.],
   [  2.,   3.,  nan,   6.,   7.,   5.,  14.,  15.,  13.],
   [  3.,  nan,  nan,   7.,   4.,   6.,  15.,  12.,  14.],
   [  4.,   5.,   7.,   8.,   9.,  11.,  nan,  nan,  nan],
   [  5.,   6.,  nan,   9.,  10.,   8.,  nan,  nan,  nan],
   [  6.,   7.,  nan,  10.,  11.,   9.,  nan,  nan,  nan],
   [  7.,  nan,  nan,  11.,   8.,  10.,  nan,  nan,  nan],
   [  8.,   9.,  11.,  12.,  13.,  15.,  nan,  nan,  nan],
   [  9.,  10.,  nan,  13.,  14.,  12.,  nan,  nan,  nan],
   [ 10.,  11.,  nan,  14.,  15.,  13.,  nan,  nan,  nan],
   [ 11.,  nan,  nan,  15.,  12.,  14.,  nan,  nan,  nan],
   [ 12.,  13.,  15.,  nan,  nan,  nan,  nan,  nan,  nan],
   [ 13.,  14.,  nan,  nan,  nan,  nan,  nan,  nan,  nan],
   [ 14.,  15.,  nan,  nan,  nan,  nan,  nan,  nan,  nan],
   [ 15.,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan]])

The idea is that neighbors >= neighbors[:, [0]] lists only those who have bigger number than the cell itself.