Probably to give an idea of origin of the problem, I should start with that this is a surface diffusion problem.
For example I have a matrix of values:
| 0 0 0 0 3 2 |
| 0 0 0 1 5 5 |
| 0 1 0 2 4 5 |
| 2 5 1 0 0 1 |
| 1 5 5 1 0 0 |
| 0 1 2 0 0 0 |
And for every cell I need to calculate a Laplace operator(simply sum of the neighbours minus cell value*4) for diffusion modelling. What I started with was Numpy's np.roll. It works very efficiently and it can be speeded up by shifting axes manually:
grid = np.array((6,6))
grid out *= -4
grid_out[1:, :] += grid[:-1, :]
grid_out[:-1, :] += grid[1:, :]
grid_out[:, 1:] += grid[:, :-1]
grid_out[:, :-1] += grid[:, 1:]
grid += grid_out
The problem here is that cells without values (<0) should 'mirror' the value of the neighbouring cell, so that they don't influence the outcoming value.
A seemingly possible solution could be pre-filling the matrix with mirrored values and running the code above. A bug arrises here, because, for example, for a cell grid[2,2] it is not clear which value to mirror 1 or 2. Thus 'mirroring' could be pre-determined only for one axis (which demands filling and clearing before and after calculation respectively) at a time and of course for every cell individually.
I would like to utilise Numpy's versatility and efficiency with arrays, but so far I found only how to do simple arithmetics on arrays. The closest function in Numpy I have found is np.where which is basically conditional cell-wise array processing and is seemingly exactly what I need, but I didn't figure out how to make it look into the neighbouring cells.
Performance here plays a big role, as this will be used for larger grids of 3d arrays, thus classical iterating is a not anticipated. The question is if this could be implemented with Numpy, because its faster than straight forward iterating with a bunch of if's even with Numba on, which I already tried.
