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)

但在其他方面与接受的答案相同。

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接受非平面索引，或者如果有一个 clip ，这一切都会更容易mode='ignore'。