Sum 4D array efficiently according to indices contained in another array

Question

I have a 4D array, a series of cubes essentially. These cubes are mostly filled with zeroes apart from sub-cubes of values of which I know the locations. I need to sum all these cubes together into one cube. I can do this simply with np.sum along axis=3 but this is part of a monte carlo process and is done a huge number of times. I was wondering since I know where the sub-cubes are in the cubes can I sum them more efficiently since the majority of the summing operation would be adding zeros- these cube are of significant size (>100^3) so if I can it will be a huge saving. I'm newish to Python/Numpy and am finding it hard to get out of the loop mindset! In general am looking for a way to manipulate large n-dimensional arrays together but only in certain parts. I realise masked arrays come to mind here but I have tried and don't think it offers any speed up in this case; unless I'm way off there!

Edit: Here are three terrible versions of what I am trying to do - probably won't make much sense out of context - basically involves calculating the effects of multiple charges at a distance from each other but each charge does not affect itself. One of these determines distance on the fly, two use precalculated arrays with the information in but it must be 'aligned' and summed - again probably won't make sense out of context but this is what I have so far

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

Did you try to use sparse matrix, since your cubes contain a lot of zeros? — JPG
– JPG, Commented Jan 28, 2014 at 18:00

Saullo G. P. Castro · Accepted Answer · 2014-01-28 19:42:56Z

1

If you know were the sub-cubes are you can use fancy indexing to sum only where you need, for example:

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)

And here are the indices where I want to sum:

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]

Which will sum only at points: p1(0,0,1,1), p2(2,1,5,6), p3(4,3,8,7) and p4(6,7,9,8).

answered Jan 28, 2014 at 19:42

Saullo G. P. Castro

59.4k28 gold badges191 silver badges244 bronze badges

Sign up to request clarification or add additional context in comments.

Collectives™ on Stack Overflow

Sum 4D array efficiently according to indices contained in another array

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related