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.
%
%
% See also cumsum, cumprod, sum, prod.
%
% 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