CODEWITHSUNDEEP
pythonnumpymatrixcumsum
fact
fact

Reputation: 557

Inverse cumsum for numpy

A is a ((d,e)) numpy array. I compute a ((d,e)) numpy array B where I compute the entry B[i,j] as follows

b=0
for k in range(i+1,d):
    for l in range(j+1,e):
        b=b+A[k,l]
B[i,j]=b

In other words, B[i,j] is the sum of A[k,l] taken over all indices k>i, l>j; this is sort of the opposite of the usual cumsum applied to both axis. I am wondering if there is a more elegant and faster way to do this (e.g. using np.cumsum)?

Upvotes: 1

Views: 920

Answers (1)

farenorth
farenorth

Reputation: 10781

Assuming you're trying to do this:

A = np.arange(15).reshape((5, -1))

def cumsum2_reverse(arr):
    out = np.empty_like(arr)
    d, e = arr.shape
    for i in xrange(d):
        for j in xrange(e):
            b = 0
            for k in xrange(i + 1, d):
                for l in xrange(j + 1, e):
                    b += arr[k, l]
            out[i, j] = b
    return out

Then if you do,

In [1]: A_revsum = cumsum2_reverse(A)

In [2]: A_revsum
Out[2]: 
array([[72, 38,  0],
      [63, 33,  0],
      [48, 25,  0],
      [27, 14,  0],
      [ 0,  0,  0]])

You could use np.cumsum on the reverse-ordered arrays to compute the sum. For example, at first you might try something similar to what @Jaime suggested:

In [3]: np.cumsum(np.cumsum(A[::-1, ::-1], 0), 1)[::-1, ::-1]
Out[3]:
array([[105,  75,  40],
       [102,  72,  38],
       [ 90,  63,  33],
       [ 69,  48,  25],
       [ 39,  27,  14]])

Here we remember that np.cumsum starts with the value in the first column (in this case last column), so to ensure zeros there, you could shift the output of this operation. This might look like:

def cumsum2_reverse_alt(arr):
    out = np.zeros_like(arr)
    out[:-1, :-1] = np.cumsum(np.cumsum(arr[:0:-1, :0:-1], 0), 1)[::-1, ::-1]
    return out

This gives the same values as above.

In [4]: (cumsum2_reverse(A) == cumsum2_reverse_alt(A)).all()
Out[4]: True

Note, that the one that utilizes np.cumsum is much faster for large arrays. For example:

In [5]: A=np.arange(3000).reshape((50, -1))

In [6]: %timeit cumsum2_reverse(A)
1 loops, best of 3: 453 ms per loop

In [7]: %timeit cumsum2_reverse_alt(A)
10000 loops, best of 3: 24.7 us per loop

Upvotes: 2

Related Questions