Code covered by the BSD License

# RCUMSUMC

by

### Matt Fig (view profile)

RCUMSUM cumulative sum of elements, restarted after every zero.

rcumprod(A,D)
```function [A] = rcumprod(A,D)
%RCUMPROD cumulative sum of elements, restarted after every zero.
% For vectors, RCUMPROD(X) is a vector containing the cumulative sum of
% the elements of X, with the summation restarting after a zero is
% encountered. For matrices, RCUMPROD(X) is a matrix the same size
% as X containing the cumulative restarted sums over each column.  For N-D
% arrays, RCUMPROD(X) operates along the first non-singleton dimension.
%
% RCUMPROD(X,DIM) works along the dimension DIM.
%
%     Example: If X = [0 1 2
%                      3 0 4
%                      1 1 3]
%
%     then rcumprod(X,1) is [0 1 2  and rcumprod(X,2) is [0 1 2
%                            3 0 8                        3 0 4
%                            3 1 24]                      1 1 3]
%
% This function has a MEX file helper function which will run only if
% compiled.  If the MEX function is not compiled, then a MATLAB
% implementation will be used.  The MEX function is RCUMPRODC.CPP
%
%    Class support: double, single.
%
%
%
% Author:  Matt Fig
% Date:  9/2/2010

S = size(A);

if nargin<2
D = find(S>1,1,'first'); % Default is the same as CUMPROD.
end

if D>ndims(A) || ~isnumeric(D) || ~isscalar(D) || isempty(A) || D<1
return  % Just return A as it is passed.
end

CLA = strcmp(class(A),{'single';'double';'logical'}); % Get input class.

if ~any(CLA)
error('Input A must be a single, double or logical array.')
end

try
if CLA(3)
A = double(A);  % Turn a logical and return.
return
end

A = rcumprodc(A,double(D));  % Call helper function.

catch
% Note this is slow!  Of course there is an analogous vectorized
% code to the one found in RCUMSUM, but it has overflow problems even
% with midsize arrays.  It also has numerical problems.  Thus this code
% is used, although one really should just compile the MEX function
% in both cases!
L = 1:length(S);
L(D) = [];
A = permute(A,[D L]);  % Work down rows.
H = 2:numel(A);
H = H(logical(mod(H-1,S(D))));  % These are the changed indexes.

for ii = H  % O.K., this will be slow.
if A(ii-1)
A(ii) = A(ii-1) * A(ii);
end
end

A = ipermute(A,[D L]);  % Original dimensions.
end```