Numpy 将较小的矩阵添加到较大的矩阵

Question

我有很大的 3D 矩阵来指示代理在 3D 空间中的位置。如果上面没有代理，则矩阵的值为 0，如果上面有代理，则为 1。 然后，我的问题是我希望代理“成长”，因为我希望它们由一个立方体（3x3x3）决定。如果已经有办法做到这一点，但是当代理靠近边界时我遇到了麻烦。 例如，我有一个位置矩阵 100x100x100，如果我知道我的代理在位置 (x, y, z) 我会这样做：

positions_matrix = numpy.zeros((100, 100, 100))
positions_matrix[x - 1: x + 2, y - 1: y + 2, z - 1: z + 2] +=  numpy.ones((3, 3, 3))

当然，在我的真实代码中，我会循环更多位置，但基本上就是这样。这是可行的，但是当代理靠近无法求和的边界时，问题就来了，因为切片得到的矩阵会小于那些矩阵。

知道如何解决它，或者 numpy 或任何其他包是否对此有实现？虽然我很确定我不是第一个面对这个问题的人，但我找不到它。

Answer 1

解决问题的更程序化的方式：

import numpy as np


m = np.zeros((100, 100, 100))

slicing = tuple(
    slice(max(0, x_i - 1), min(x_i + 2, d - 1))
    for x_i, d in zip((x, y, z), m.shape))
ones_shape = tuple(s.stop - s.start for s in slicing)

m[slicing] += np.ones(ones_shape)

但在其他方面与accepted answer 相同。

Answer 2

你应该在上下界进行切割，使用类似的东西：

import numpy as np
m = np.zeros((100, 100, 100))

x_min, x_max = np.max([0, x-1]), np.min([x+2, m.shape[0]-1])
y_min, y_max = np.max([0, y-1]), np.min([y+2, m.shape[1]-1])
z_min, z_max = np.max([0, z-1]), np.min([z+2, m.shape[2]-1])

m[x_min:x_max, y_min:y_max, z_min:z_max] += np.ones((x_max-x_min, y_max-y_min, z_max-z_min))

Answer 3

有一个使用np.put 及其'clip' 选项的解决方案。 它只需要一点技巧，因为该函数需要扁平矩阵中的索引；幸运的是，np.ravel_multi_index 函数完成了这项工作：

import itertools
import numpy as np
x, y, z = 2, 0, 4
positions_matrix = np.zeros((100,100,100))

indices = np.array( list( itertools.product( (x-1, x, x+1), (y-1, y, y+1), (z-1, z, z+1)) ))
flat_indices = np.ravel_multi_index(indices.T, positions_matrix.shape, mode='clip')
positions_matrix.put(flat_indices, 1+positions_matrix.take(flat_indices))
# positions_matrix[2,1,4] is now 1.0

这个解决方案的好处是您可以使用其他模式，例如'wrap'（如果您的代理住在甜甜圈上;-) 或周期性空间）。

我将解释它如何在较小的 2D 矩阵上工作：

import itertools
import numpy as np

positions_matrix = np.zeros((8,8))
ones = np.ones((3,3))
x, y = 0, 4

indices = np.array( list( itertools.product( (x-1, x, x+1), (y-1, y, y+1) )))
# array([[-1,  3],
#        [-1,  4],
#        [-1,  5],
#        [ 0,  3],
#        [ 0,  4],
#        [ 0,  5],
#        [ 1,  3],
#        [ 1,  4],
#        [ 1,  5]])

flat_indices = np.ravel_multi_index(indices.T, positions_matrix.shape, mode='clip')
# array([ 3,  4,  5,  3,  4,  5, 11, 12, 13])

positions_matrix.put(flat_indices, ones, mode='clip')
# positions_matrix is now:
# array([[0., 0., 0., 1., 1., 1., 0., 0.],
#        [0., 0., 0., 1., 1., 1., 0., 0.],
#        [0., 0., 0., 0., 0., 0., 0., 0.],
#        [ ...

顺便说一句，在这种情况下，mode='clip' 对于put 来说是多余的。

好吧，我刚刚作弊put 完成了一项任务。 +=1 需要 take 和 put：

positions_matrix.put(flat_indices, ones.flat + positions_matrix.take(flat_indices))
# notice that ones has to be flattened, or alternatively the result of take could be reshaped (3,3)
# positions_matrix is now: 
# array([[0., 0., 0., 2., 2., 2., 0., 0.],
#        [0., 0., 0., 2., 2., 2., 0., 0.],
#        [0., 0., 0., 0., 0., 0., 0., 0.],
#        [ ...

与其他解决方案相比，此解决方案有一个重要区别：ones 矩阵始终为 (3,3)， 这可能是也可能不是优势。 诀窍就在这个 flat_indices 列表中，它有重复的条目（剪辑的结果）。

因此，如果您在最大索引处添加非常数子矩阵，则可能需要采取一些预防措施：

x, y = 1, 7
values = 1 + np.arange(9)
indices = np.array( list( itertools.product( (x-1, x, x+1), (y-1, y, y+1) )))
flat_indices = np.ravel_multi_index(indices.T, positions_matrix.shape, mode='clip')
positions_matrix.put(flat_indices, values, mode='clip')
# positions_matrix is now:
# array([[0., 0., 0., 2., 2., 2., 1., 3.],
#        [0., 0., 0., 2., 2., 2., 4., 6.],
#        [0., 0., 0., 0., 0., 0., 7., 9.],

...您可能期望最后一列是 2 5 8。 目前，您可以处理flat_indices，例如将-1 放在越界位置。 但是如果np.put 接受非平面索引，或者如果有一个剪辑mode='ignore'，这一切都会更容易。