9

I have a numpy/pandas list of values:

a = np.random.randint(-100, 100, 10000)
b = a/100

I want to apply a custom cumsum function, but I haven't found a way to do it without loops. The custom function sets an upper limit of 1 and lower limit of -1 for the cumsum values, if the "add" to sum is beyond these limits the "add" becomes 0.

In the case that sum is between the limits of -1 and 1 but the "added" value would break beyond the limits, the "added" becomes the remainder to -1 or 1.

Here is the loop version:

def cumsum_with_limits(values):
    cumsum_values = []
    sum = 0
    for i in values:
        if sum+i <= 1 and sum+i >= -1: 
            sum += i
            cumsum_values.append(sum)
        elif sum+i >= 1:
            d = 1-sum # Remainder to 1
            sum += d
            cumsum_values.append(sum)
        elif sum+i <= -1:
            d = -1-sum # Remainder to -1
            sum += d
            cumsum_values.append(sum)

    return cumsum_values

Is there any way to vectorize this? I need to run this function on large datasets and performance is my current issue. Appreciate any help!


Update: Fixed the code a bit, and a little clarification for the outputs: Using np.random.seed(0), the first 6 values are:

b = [0.72, -0.53, 0.17, 0.92, -0.33, 0.95]

Expected output:

o = [0.72, 0.19, 0.36, 1, 0.67, 1]
9
  • If I understand correctly, cumsum_with_limits gives you the list of values such that their cumsum ever goes below -1 or above +1, right? So, what you want is that array of numbers, not the cumsum itself, is that correct?
    – javidcf
    Commented Nov 1, 2018 at 16:00
  • Yes correct, the output is the list of values not the cumsum itself e.g. [0, 0.3, 0.6, , 0.8, 1, 1, 1, 1, 1, 1, 0.6, 0.4, 0.1, -0.3, -0.6, -1, -1, -1, -1, -0.5], it can't go beyond 1 or -1 Commented Nov 1, 2018 at 16:03
  • Can you show a concrete example of an input and expected output that showcases the way you handle all the bounds? Commented Nov 1, 2018 at 16:11
  • I'm pretty sure there's some version of reduceat that would do this for you. Just have to figure out how to phrase it. Commented Nov 1, 2018 at 16:13
  • Yep, input [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2], output [0.2, 0.4, 0.6, 0.8, 1, 1, 1] Commented Nov 1, 2018 at 16:13

3 Answers 3

6

Loops aren't necessarily undesirable. If performance is an issue, consider numba. There's a ~330x improvement without materially changing your logic:

from numba import njit

np.random.seed(0)
a = np.random.randint(-100, 100, 10000)
b = a/100

@njit
def cumsum_with_limits_nb(values):
    n = len(values)
    res = np.empty(n)
    sum_val = 0
    for i in range(n):
        x = values[i]
        if (sum_val+x <= 1) and (sum_val+x >= -1):
            res[i] = x
            sum_val += x
        elif sum_val+x >= 1:
            d = 1-sum_val # Remainder to 1
            res[i] = d
            sum_val += d
        elif sum_val+x <= -1:
            d = -1-sum_val # Remainder to -1
            res[i] = d
            sum_val += d
    return res

assert np.isclose(cumsum_with_limits(b), cumsum_with_limits_nb(b)).all()

If you don't mind sacrificing some performance, you can rewrite this loop more succinctly:

@njit
def cumsum_with_limits_nb2(values):
    n = len(values)
    res = np.empty(n)
    sum_val = 0
    for i in range(n):
        x = values[i]
        next_sum = sum_val + x
        if np.abs(next_sum) >= 1:
            x = np.sign(next_sum) - sum_val
        res[i] = x
        sum_val += x
    return res

With similar performance to nb2, here's an alternative (thanks to @jdehesa):

@njit
def cumsum_with_limits_nb3(values):
    n = len(values)
    res = np.empty(n)
    sum_val = 0
    for i in range(n):
        x = min(max(sum_val + values[i], -1) , 1) - sum_val
        res[i] = x
        sum_val += x
    return res

Performance comparisons:

assert np.isclose(cumsum_with_limits(b), cumsum_with_limits_nb(b)).all()
assert np.isclose(cumsum_with_limits(b), cumsum_with_limits_nb2(b)).all()
assert np.isclose(cumsum_with_limits(b), cumsum_with_limits_nb3(b)).all()

%timeit cumsum_with_limits(b)      # 12.5 ms per loop
%timeit cumsum_with_limits_nb(b)   # 40.9 µs per loop
%timeit cumsum_with_limits_nb2(b)  # 54.7 µs per loop
%timeit cumsum_with_limits_nb3(b)  # 54 µs per loop
7
  • @FrancWeser, No problem. It's probably because you mention you don't want a loop. But sometimes it's not necessarily a bad idea. For what it's worth, a specialized NumPy method is probably preferable, so please don't mark this as accepted until others have tried to figure it out.
    – jpp
    Commented Nov 1, 2018 at 16:15
  • 1
    I know I did misunderstood the question , so that downvote for me is fine , but why this answer get downvote as well ?
    – BENY
    Commented Nov 1, 2018 at 16:16
  • @W-B not sure, but I think your answer is the correct one. The user posted an example input output above that actually matches what you were doing :P. Very confusing Commented Nov 1, 2018 at 16:17
  • 1
    @user3483203 ummm, that is why I delete it :-) I am stepping away from unclear question
    – BENY
    Commented Nov 1, 2018 at 16:18
  • 1
    That's an amazing performance improvement. Thanks a lot! Curious if there are still ways to solve it without looping Commented Nov 1, 2018 at 17:06
5

Start with a regular cumsum:

b = ...
s = np.cumsum(b)

Find the first clip point:

i = np.argmax((s[0:] > 1) | (s[0:] < -1))

Adjust everything that follows:

s[i:] += (np.sign(s[i]) - s[i])

Rinse and repeat. This still requires a loop, but only over the adjustment points, which is generally expected to be much smaller than the total number of array size.

b = ...
s = np.cumsum(b)
while True:
    i = np.argmax((s[0:] > 1) | (s[0:] < -1))
    if np.abs(s[i]) <= 1:
        break
    s[i:] += (np.sign(s[i]) - s[i])

I still haven't found a way to completely pre-compute the adjustment points up front, so I would have to guess that the numba solution will be faster than this, even if it you compiled this with numba.

Starting with np.seed(0), your original example has 3090 adjustment points, which is approximately 1/3. Unfortunately, with all the temp arrays and extra sums, that makes the algorithmic complexity of my solution tend to O(n2). This is completely unacceptable.

0

I thought I had already answered the generic question of "cumulative sum with bounds" in the past, but I can't find it.

This solution also uses numba and is a bit more general (custom bounds) and concise than the ones given by @jpp.

It operates on the OP's problem (10K values, bounds at -1, 1) in 40 µs.

import numpy as np
from numba import njit

@njit
def cumsum_clip(a, xmin=-np.inf, xmax=np.inf):
    res = np.empty_like(a)
    c = 0
    for i in range(len(a)):
        c = min(max(c + a[i], xmin), xmax)
        res[i] = c
    return res

Example

np.random.seed(0)
x = np.random.randint(-100, 100, 10_000) / 100

>>> x[:6]
array([ 0.72, -0.53,  0.17,  0.92, -0.33,  0.95])

>>> cumsum_clip(x, -1, 1)[:6]
array([0.72, 0.19, 0.36, 1.  , 0.67, 1.  ])

%timeit cumsum_clip(x, -1, 1)
39.3 µs ± 31 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

Note: you can specify other bounds, e.g.:

>>> cumsum_clip(x, 0, 1)[:10]
array([0.72, 0.19, 0.36, 1.  , 0.67, 1.  , 1.  , 0.09, 0.  , 0.  ])

Or omit one of the bounds (for example here specifying only an upper bound):

>>> cumsum_clip(x, xmax=1)[:10]
array([ 0.72,  0.19,  0.36,  1.  ,  0.67,  1.  ,  1.  ,  0.09, -0.7 , -1.34])

Of course, it preserves the original dtype:

np.random.seed(0)
x = np.random.randint(-10, 10, 10)
>>> cumsum_clip(x, 0, 10)
array([ 2,  7,  0,  0,  0,  0,  0,  9, 10,  4])

>>> cumsum_clip(x, 0, 10).dtype
dtype('int64')

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.