CODEWITHSUNDEEP
pythonnumpycumsum
Touki
Touki

Reputation: 843

"np.cumsum() like" - iterative, on actual values

I'm searching a way to implement with numpy this piece of python code:

  N = np.arange(5)

  for i in range(5):
     for k in range(i):
        N[i] += N[k]

Assuming that I work in fact on big 2-D arrays (1300*1300).

np.cumsum() provide a good way, on one axis N[0][i] or N[i][0], except that it only sums the values of the original array, not of the evolving array.

I can't figure a way to do that. Any Idea ?

@Edit :

To make things clear :

With 1-D array the loop give

Out[89]: array([ 0,  1,  3,  7, 15])

With cumsum :

array([ 0,  1,  3,  6, 10])

With a 2-D, it would give something like :

N = np.arange(25).reshape(5,5)
for i in range(len(N)):
    N = np.cumsum(N, axis=i)

Upvotes: 1

Views: 395

Answers (2)

Jaime
Jaime

Reputation: 67457

If you work out the result of your loop, starting with a sequence a[n], and what comes out of the two for loops is a sequence b[n], then:

b[n] = a[n] + a[n-1] + 2*a[n-2] + 4*a[n-3] + ... + 2**(n-2)*a[0] =
       a[n] + c[n-1]

where I am defining:

c[n-1] = a[n-1] + 2*a[n-2] + 4*a[n-3] + ... + 2**(n-2)*a[0]

With this last expression in mind, there are ways to vectorize your double loop. But note the very large factors(2**(n-2)) by which you have to multiply the items in your sequence. If your sequence has positive and negative terms, these may cancel out and return reasonable numbers. But if you have an array of over 1000 positive elements, then you are going to overflow any numpy dtype.

So for short sequences of less than 30 items, perhaps 60 if you force using int64's, the following will work faster than your for loops:

def evolving_cumsum(arr):
    arr = np.array(arr) # makes a copy of the data
    pows = 2**np.arange(len(arr))[::-1]
    c = np.cumsum(arr*pows)
    c /= pows
    arr[1:] += c[:-1]
    return arr

>>> a = np.arange(10)
>>> evolving_cumsum(a)
array([  0,   1,   3,   7,  15,  31,  63, 127, 255, 511])
>>> for i in range(len(a)):
...     for k in range(i):
...         a[i] += a[k]
... 
>>> a
array([  0,   1,   3,   7,  15,  31,  63, 127, 255, 511])

But in general I am afraid you are going to have to keep your loops.

Upvotes: 2

Shankar
Shankar

Reputation: 3624

import numpy as np
N = np.arange(5)
print np.cumsum(N) # outputs [ 0  1  3  6 10]

for i in range(5):
    for k in range(i):
        N[i] += N[k]

print np.cumsum(N) # outputs [ 0  1  4 11 26]

What do you mean by evolving array?

Upvotes: 0

Related Questions