cumsum with upper & lower limits?

Question

I'd like to find a vectorized way to calculate the cumulative sums of a vector, but with upper and lower limits.

In my case, the input only contains 1's and -1's. You can use this assumption in your answer. Of course, a more general solution is also welcome.

For example:

x     = [1 1 1 1 -1 -1 -1 -1 -1 -1];
upper = 3;
lower = 0;
s     = cumsum(x)                    %// Ordinary cumsum.
s =
     1     2     3     4     3     2     1     0    -1    -2

y     = cumsumlim(x, upper, lower)   %// Cumsum with limits.
y =
     1     2     3     3     2     1     0     0     0     0
                 ^                       ^
                 |                       |
            upper limit             lower limit

When the cumulative sum reaches the upper limit (at the 3rd element), it won't increase anymore. Likewise, when the cumulative sum reaches the lower limit (at the 7th element), it won't decrease anymore. A for-loop version would be like this:

function y = cumsumlim(x, upper, lower)

y = zeros(size(x));
y(1) = x(1);

for i = 2 : numel(x)
    y(i) = y(i-1) + x(i);
    y(i) = min(y(i), upper);
    y(i) = max(y(i), lower);
end

end

Do you have any ideas?

Trimeter

I don't quite understand the example output you have showed. Can you be more verbose and explain how you arrived at your desired output? How exactly do the upper and lower limits come into play with your "cumsum" function? – Trafficator

Does x contain 1's and -1's only? – Sarcocarp

@BenW either specify limits on what x can contain or choose a more representative x. Firstly have it go more than 1 above and/or below the limits, have it breach the same limit more than once and most importantly as Divakar mentions, if x can contain other numbers please include some – Kurt

@Trafficator , I've added more explanation. Hope that will be helpful. – Trimeter

@Divakar, in my case, the input only contains 1's and -1's. And also, I've added a for-loop version code of what I wanted to do. Hope that will help to explain. – Trimeter

Should this cumsum be able reach max more than once? – Barbitone

@HamtaroWarrior, yes, of course. – Trimeter

Just asking because it'd be quite quick otherwise ;) – Barbitone

@Ben.W, I suggest you accept Luis' answer (it's a nice gesture when someone answers one of your questions). It's the tick mark on the left side of the answer. (You could also accept mine if you for some reason didn't like his solution). – Hendrik

@StewieGriffin, thanks for the reminder. I found both answers useful, so I up-voted both. However, I haven't found a solution to my problem yet, that's why I didn't accept any answer. Anyway, it seems like no one else is going to answer, so I'm happy to accept Luis' answer. – Trimeter

Answer 1

This is a somewhat hackish solution, but perhaps worth mentioning.

You can do the sum using a signed integer data type, and exploit the inherent limits of that data type. For this to work, the input needs to be converted to that integer type and multiplied by the appropiate factor, and an initial offset needs to be applied. The factor and offset are chosen as a function of lower and upper. After cumsum, the multiplication and offset are undone to obtain the desired result.

In your example, data type int8 suffices; and the required factor and offset are 85 and -128 respectively:

x = [1 1 1 1 -1 -1 -1 -1 -1 -1];
result = cumsum([-128 int8(x)*85]); %// integer sum, with factor and initial offset
result = (double(result(2:end))+128)/85; %// undo factor and offset

which gives

result =
     1     2     3     3     2     1     0     0     0     0

Enrol

Answer 2

I won't provide you with a magic vectorized way to do this, but I'll provide you with some data that probably will help you get on with your work.

Your cumsumlim function is very fast!

tic
for ii = 1:100
    y = cumsumlim(x,3,0);
end
t = toc;
disp(['Length of vector: ' num2str(numel(x))])
disp(['Total time for one execution: ' num2str(t*10), ' ms.'])
Length of vector: 65000
Total time for one execution: 1.7965 ms.

I really doubt this is your bottleneck. Have you tried profiling the code?

Hendrik