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
    


Answer: 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)`.

