根据另一个数组中包含的索引，高效地加总4D数组

Question

我有一个4D数组，基本上是一系列立方体。除了我知道其位置的值的子立方体之外，这些立方体大多是用零填充的。我需要把所有这些立方相加成一个立方。我只需要np.总和沿轴=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

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)

下面是我要总结的指数：

^{pr2}$

它只在点上求和：p1(0,0,1,1)，p2(2,1,5,6)，p3(4,3,8,7)和{}。在