Numpy将较小的矩阵与较大的矩阵相加

Question

我有一个很大的3D矩阵来表示代理在3D空间中的位置。如果矩阵上没有代理，则矩阵的值为0；如果矩阵上有代理，则矩阵的值为1。然后，我的问题是，我希望代理在某种意义上“增长”，即我希望它们由比方(3x3x3)的1确定。如果已经有办法做到这一点，但是当智能体靠近边界时，我遇到了麻烦。例如，我有一个位置矩阵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))

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 2

一种稍微更程序化的解决问题的方法：

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 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'，一切都会更容易。