高效地按照另一个数组中包含的索引对4D数组进行求和

Question

我有一个4D数组，本质上是一系列立方体。这些立方体大多数都填充了零值，除了我知道位置的值子立方体。我需要将所有这些立方体相加成一个立方体。我可以简单地使用np.sum沿axis=3完成此操作，但这是蒙特卡罗过程的一部分，并且会重复进行很多次。我想知道，由于我知道立方体中的子立方体在哪里，我能否更有效地将它们相加，因为大多数求和操作将是添加零-这些立方体的大小很大（>100 ^ 3），因此如果我可以节省时间就最好了。我对Python / Numpy还比较新，发现很难摆脱循环思维方式！总之，我正在寻找一种方法，可以仅在某些部分上操作大型n维数组。我意识到遮罩数组在这里很合适，但我已经尝试过了，并不认为它在这种情况下有任何加速作用；除非我完全错了！ 编辑：这里有三个糟糕的版本，其中包含我正在尝试做的事情-可能在上下文之外不太清楚-基本上涉及计算彼此之间距离的多个电荷的影响，但每个电荷不会影响自己。其中一个在运行时确定距离，另外两个使用预先计算的数组，其中包含信息，但必须“对齐”并相加-同样可能在上下文之外没有意义，但这是我到目前为止所做的

def coloumbicForces3(carrierArray, cubeEnergeticDisorderArray, array):
cubeLen=100

offsetArray=np.array([[0,0,0],[0,0,+1],[0,0,-1],[0,+1,0],[0,-1,0],[+1,0,0],[-1,0,0]])
indices=np.zeros((7,3,len(carrierArray)), dtype=np.int32)
indices1=a=indices.reshape((7*len(carrierArray),3))
tIndices=cubeLen-indices
superimposedArray=np.zeros((cubeLen,cubeLen,cubeLen,2+2), dtype=myFloat)
sumArray=np.zeros((cubeLen,cubeLen,cubeLen))
for i in range(len(carrierArray)):
    indices[:,:,i]=offsetArray+carrierArray[i,1:4]

for c, carrierC in enumerate(carrierArray[:,0]):
    for k, carrierK in enumerate(carrierArray[:,0]):
        if c==k:
            continue
        for (x,y,z) in indices[:,:,k]:
        #print c, indices[:,:,c]
            if(carrierC==1):
                superimposedArray[x,y,z,c]=cubeEnergeticDisorderArray[x,y,z,c]=-1*array[cubeLen-x,cubeLen-y, cubeLen-z]
            else:
                superimposedArray[x,y,z,c]=cubeEnergeticDisorderArray[x,y,z,c]=array[cubeLen-x,cubeLen-y, cubeLen-z]



b = np.ascontiguousarray(a).view(np.dtype((np.void, a.dtype.itemsize * a.shape[1])))
_,idx = np.unique(b, return_index=True)
aUnique=a[idx]

for (i,j,k) in aUnique:
    sumArray[i,j,k]=np.sum(superimposedArray[i,j,k])

for c, carrierC in enumerate(carrierArray[:,0]):
    for (i,j,k) in aUnique:
        cubeEnergeticDisorderArray[i,j,k,c]=cubeEnergeticDisorderArray[i,j,k,-2]+sumArray[i,j,k]

return cubeEnergeticDisorderArray


def coloumbicForces(carrierArray, cubeEnergeticDisorderArray, array):
cubeLen= len(cubeEnergeticDisorderArray[0,:,:,0])
superimposedArray=np.zeros((cubeLen,cubeLen,cubeLen,2+2), dtype=myFloat)

for k, carrier in enumerate(carrierArray[:,0]):
    superimposedArray[:,:,:,k]=cubeEnergeticDisorderArray[:,:,:,k]=array[cubeLen-carrierArray[k,1]:2*cubeLen-carrierArray[k,1],cubeLen-carrierArray[k,2]:2*cubeLen-carrierArray[k,2],cubeLen-carrierArray[k,3]:2*cubeLen-carrierArray[k,3]]

    if (carrier==1):
        a=superimposedArray[:,:,:,k]
        b=cubeEnergeticDisorderArray[:,:,:,k]
        superimposedArray[:,:,:,k]=ne.evaluate("a*-1")
        cubeEnergeticDisorderArray[:,:,:,k]=ne.evaluate("b*-1")

sumArray=ne.evaluate("sum(superimposedArray, axis=3)")

for k, carrier in enumerate(carrierArray[:,0]):
    a=cubeEnergeticDisorderArray[:,:,:,k]
    b=cubeEnergeticDisorderArray[:,:,:,-2]
    cubeEnergeticDisorderArray[:,:,:,k]=ne.evaluate("sumArray-a+b")

return cubeEnergeticDisorderArray

def coloumbicForces2(carrierArray, cubeEnergeticDisorderArray, array):
x0=carrierArray[0,1]
y0=carrierArray[0,2]
z0=carrierArray[0,3]
x1=carrierArray[1,1]
y1=carrierArray[1,2]
z1=carrierArray[1,3]

cubeEnergeticDisorderArray[x0,y0,z0,0]=cubeEnergeticDisorderArray[x0,y0,z0,-2]-(1.60217657e-19)*2995850595.79/(distance([x0,y0,z0], carrierArray[1,1:4])*1e-9)
cubeEnergeticDisorderArray[x0-1,y0,z0,0]=cubeEnergeticDisorderArray[x0-1,y0,z0,-2]-(1.60217657e-19)*2995850595.79/(distance([x0-1,y0,z0], carrierArray[1,1:4])*1e-9)
cubeEnergeticDisorderArray[x0+1,y0,z0,0]=cubeEnergeticDisorderArray[x0+1,y0,z0,-2]-(1.60217657e-19)*2995850595.79/(distance([x0+1,y0,z0], carrierArray[1,1:4])*1e-9)
cubeEnergeticDisorderArray[x0,y0-1,z0,0]=cubeEnergeticDisorderArray[x0,y0-1,z0,-2]-(1.60217657e-19)*2995850595.79/(distance([x0,y0-1,z0], carrierArray[1,1:4])*1e-9)
cubeEnergeticDisorderArray[x0,y0+1,z0,0]=cubeEnergeticDisorderArray[x0,y0+1,z0,-2]-(1.60217657e-19)*2995850595.79/(distance([x0,y0+1,z0], carrierArray[1,1:4])*1e-9)
cubeEnergeticDisorderArray[x0,y0,z0-1,0]=cubeEnergeticDisorderArray[x0,y0,z0-1,-2]-(1.60217657e-19)*2995850595.79/(distance([x0,y0,z0-1], carrierArray[1,1:4])*1e-9)
cubeEnergeticDisorderArray[x0,y0,z0+1,0]=cubeEnergeticDisorderArray[x0,y0,z0+1,-2]-(1.60217657e-19)*2995850595.79/(distance([x0,y0,z0+1], carrierArray[1,1:4])*1e-9)

cubeEnergeticDisorderArray[x1,y1,z1,1]=cubeEnergeticDisorderArray[x1,y1,z1,-2]+(1.60217657e-19)*2995850595.79/(distance([x1,y1,z1], carrierArray[0,1:4])*1e-9)
cubeEnergeticDisorderArray[x1-1,y1,z1,1]=cubeEnergeticDisorderArray[x1-1,y1,z1,-2]+(1.60217657e-19)*2995850595.79/(distance([x1-1,y1,z1], carrierArray[0,1:4])*1e-9)
cubeEnergeticDisorderArray[x1+1,y1,z1,1]=cubeEnergeticDisorderArray[x1+1,y1,z1,-2]+(1.60217657e-19)*2995850595.79/(distance([x1+1,y1,z1], carrierArray[0,1:4])*1e-9)
cubeEnergeticDisorderArray[x1,y1-1,z1,1]=cubeEnergeticDisorderArray[x1,y1-1,z1,-2]+(1.60217657e-19)*2995850595.79/(distance([x1,y1-1,z1], carrierArray[0,1:4])*1e-9)
cubeEnergeticDisorderArray[x1,y1+1,z1,1]=cubeEnergeticDisorderArray[x1,y1+1,z1,-2]+(1.60217657e-19)*2995850595.79/(distance([x1,y1+1,z1], carrierArray[0,1:4])*1e-9)
cubeEnergeticDisorderArray[x1,y1,z1-1,1]=cubeEnergeticDisorderArray[x1,y1,z1-1,-2]+(1.60217657e-19)*2995850595.79/(distance([x1,y1,z1-1], carrierArray[0,1:4])*1e-9)
cubeEnergeticDisorderArray[x1,y1,z1+1,1]=cubeEnergeticDisorderArray[x1,y1,z1+1,-2]+(1.60217657e-19)*2995850595.79/(distance([x1,y1,z1+1], carrierArray[0,1:4])*1e-9)
return cubeEnergeticDisorderArray

- user1857560

---

你尝试使用稀疏矩阵了吗？因为你的立方体包含很多零。 - JPG

稀疏矩阵仅限于2D。 - user1857560

Answer 1

如果您知道子立方体的位置，您可以使用高级索引仅在需要的位置进行求和，例如：

import numpy as np
from numpy.random import random

c1 = random((10, 10, 10, 10))
c2 = random((10, 10, 10, 10))
c3 = np.zeros_like(c2)

这里是我想要求和的索引：

i1 = [0, 2, 4, 6]
i2 = [0, 1, 3, 7]
i3 = [1, 5, 8, 9]
i4 = [1, 6, 7, 8]

c3[i1,i2,i3,i4] = c1[i1,i2,i3,i4] + c2[i1,i2,i3,i4]

这将仅在以下点求和：p1(0,0,1,1)、p2(2,1,5,6)、p3(4,3,8,7) 和 p4(6,7,9,8)。

- Saullo G. P. Castro