在 2d 中显示 4d 矩阵的更多“pythonic”方式

Question

我想将 4d 矩阵绘制为带有索引的 2d 矩阵：

[i][j][k][l] --> [i * nj + j][ k * nl + l]

我有一个工作版本here。

这是我想要的工作，但它不是很优雅。我研究了“重塑”，但这并不是我想要的，或者我使用不正确。

给定一个形状为 (100000,4) 的 4d 数组“r”，我要替换的相关 sn-p 是：

def transform(i,j,k,l, s1, s2):
    return [i * s1 + j, k * s2 + l]

nx = 5
ny = 11
iedges=np.linspace(0,100, nx)
jedges=np.linspace(0, 20, ny)
bins = ( iedges,jedges,iedges,jedges )

H, edges = np.histogramdd(r, bins=bins )

H2 = np.zeros(( (nx-1)*(ny-1),(nx-1)*(ny-1)))
for i in range(nx-1):
    for j in range(ny-1):
        for k in range(nx-1):
            for l in range(ny-1):
                x,y = transform(i,j,k,l,ny-1,ny-1)                
                H2[x][y] = H[i][j][k][l]

在这种情况下，H2 的值将对应于 H 的值，但是条目

i,j,k,l

将显示为

i*ny + j, k * ny + l

。

示例图：

Answer 1

看起来你可以从x,y 计算出i,j,k,l？这应该是这样的：

from functools import partial

def get_ijkl(x, y, s1, s2):
    # "Reverse" of `transform`
    i, j = divmod(x, s1)
    k, l = divmod(y, s2)
    return (i, j, k, l)


def get_2d_val(x, y, s1, s2, four_dim_array):
    return four_dim_array[get_ijkl(x, y, s1, s2)]


smaller_shape = ((nx-1)*(ny-1), (nx-1)*(ny-1))

知道这一点有几种可能的方法：

• numpy.fromfunction：

H3 = np.fromfunction(
    partial(get_2d_val, s1=ny-1, s2=ny-1, four_dim_array=H),
    shape=smaller_shape,
    dtype=int,
)
assert np.all(H2 == H3)

• indexing:

indices_to_take = np.array([
    [list(get_ijkl(x, y, ny-1, ny-1)) for x in range(smaller_shape[0])] for y in range(smaller_shape[1])
]).transpose()
H4 = H[tuple(indices_to_take)]
assert np.all(H2 == H4)

• 正如@hpaulj 回答的那样，您可以简单地重塑数组，它会更快。但是如果你有一些不同的transform 并且可以计算适当的“反向”函数，那么使用fromfunction 或自定义索引将变得有用

Answer 2

您确定reshape 不起作用吗？

我在一个小的随机 r 上运行了您的代码。 H 的非零项是：

In [13]: np.argwhere(H)                                                                                  
Out[13]: 
array([[0, 9, 3, 1],
       [1, 1, 1, 2],
       [1, 2, 1, 3],
       [2, 2, 2, 3],
       [3, 1, 1, 8]])

对于转换后的H2：

In [14]: np.argwhere(H2)                                                                                 
Out[14]: 
array([[ 9, 31],
       [11, 12],
       [12, 13],
       [22, 23],
       [31, 18]])

其中一个H 索引转换为H2 索引：

In [16]: transform(0,9,3,1,4,10)                                                                         
Out[16]: [9, 31]

如果我只是重塑H，我会得到与H2 相同的数组：

In [17]: H3=H.reshape(40,40)                                                                             
In [18]: np.argwhere(H3)                                                                                 
Out[18]: 
array([[ 9, 31],
       [11, 12],
       [12, 13],
       [22, 23],
       [31, 18]])
In [19]: np.allclose(H2,H3)                                                                              
Out[19]: True

因此，无需深入研究您的代码细节，它在我看来就像一个简单的重塑。