Code covered by the BSD License  

Highlights from
RCUMSUMC

image thumbnail

RCUMSUMC

by

 

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.
%
%
%    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

Contact us