function x = cumop(opstr,x,varargin)
%MATLAB Code Generation Private Function

%   Copyright 2002-2023 The MathWorks, Inc.
%#codegen

narginchk(2,5);
coder.inline('always');
coder.internal.prefer_const(opstr,varargin);
[dim,rev,omitnan] = processInputs(opstr,x,varargin{:});
if coder.const(coder.gpu.internal.isGpuCumSumSupported() && isreal(x) && strcmp(opstr,'cumsum'))...
        && ~coder.internal.isConstantFolding
    % We are calling the GPU Implementation of cumsum over here
    x = gpucoder.internal.cumsumImpl(x,dim,rev,omitnan);
else
    x = coder.internal.useConstantDim(3,coder.internal.ndims(x), ...
        @ndlooper,opstr,x,dim,rev,omitnan);
end

%--------------------------------------------------------------------------

function [dim,rev,omitnan] = processInputs(opstr,x,varargin)
coder.inline('always');
coder.internal.prefer_const(opstr,varargin);
eml_invariant(coder.internal.isBuiltInNumeric(x) || islogical(x), ...
    'Coder:toolbox:unsupportedClass',opstr,class(x));
dim = coder.internal.nonSingletonDim(x);
rev = false;
% Default value of omitnan is true for cummin and cummax, otherwise false.
omitnan = coder.const(strcmp(opstr,'cummin') || strcmp(opstr,'cummax'));
if nargin >= 3
    if coder.internal.isCharOrScalarString(varargin{1})
        firstopt = coder.internal.indexInt(1);
    else
        coder.internal.assertValidDim(varargin{1});
        dim = coder.internal.indexInt(varargin{1});
        firstopt = coder.internal.indexInt(2);
    end
    directionFlagSupplied = false;
    NaNFlagSupplied = false;
    matchParam = @coder.internal.matchStringParameter;
    for k = coder.unroll(firstopt:length(varargin))
        if matchParam(varargin{k},'reverse')
            coder.internal.assert(~directionFlagSupplied, ...
                'MATLAB:cumfun:duplicateDirection');
            directionFlagSupplied = true;
            rev = true;
        elseif matchParam(varargin{k},'forward')
            coder.internal.assert(~directionFlagSupplied, ...
                'MATLAB:cumfun:duplicateDirection');
            directionFlagSupplied = true;
            rev = false;
        elseif matchParam(varargin{k},'omitnan') || matchParam(varargin{k},'omitmissing')
            coder.internal.assert(~NaNFlagSupplied, ...
                'MATLAB:cumfun:duplicateNaNFlag');
            NaNFlagSupplied = true;
            omitnan = true;
        elseif matchParam(varargin{k},'includenan') || matchParam(varargin{k},'includemissing')
            coder.internal.assert(~NaNFlagSupplied, ...
                'MATLAB:cumfun:duplicateNaNFlag');
            NaNFlagSupplied = true;
            omitnan = false;
        else
            eml_invariant(false,'MATLAB:cumfun:wrongFlag');
        end
    end
end

%--------------------------------------------------------------------------

function p = minIncludeNaN(a,b)
p = min(a,b,'includenan');

function p = maxIncludeNaN(a,b)
p = max(a,b,'includenan');

%--------------------------------------------------------------------------

function x = ndlooper(opstr,x,dim,rev,omitnan)
coder.inline('always');
coder.internal.prefer_const(opstr,dim,rev,omitnan)
ND = coder.internal.indexInt(coder.internal.ndims(x));
ONE = coder.internal.indexInt(1);
if coder.isRowMajor
    thisDim = ONE;
else
    thisDim = ND;
end
% Specialize on options even if they aren't constants.
if rev
    if omitnan
        x = looper(opstr,x,dim,true,true,thisDim);
    else
        x = looper(opstr,x,dim,true,false,thisDim);
    end
else
    if omitnan
        x = looper(opstr,x,dim,false,true,thisDim);
    else
        x = looper(opstr,x,dim,false,false,thisDim);
    end
end

%--------------------------------------------------------------------------

function x = looper(opstr,x,dim,rev,omitnan,thisDim,varargin)
coder.inline('always');
coder.internal.prefer_const(opstr,dim,rev,omitnan,thisDim);
ND = coder.internal.indexInt(coder.internal.ndims(x));
ONE = coder.internal.indexInt(1);
XZERO = coder.internal.scalarEg(x);
[OP,NVAL,REPLACENANS] = coder.const(@constOptions,opstr,XZERO,omitnan);
FIRSTCALL = coder.const(nargin < 7);
UNITDIM = coder.internal.isConstTrue(size(x,dim) == 1);
if FIRSTCALL && isempty(x)
elseif coder.const(FIRSTCALL && UNITDIM)
    % With dim <= ND, this is an optimization.
    % With dim > ND it is required.
    if REPLACENANS
        x = replaceNaNs(x,NVAL,thisDim);
    end
elseif coder.const(FIRSTCALL && ...
        coder.internal.isConstTrue(isvector(x)))
    % Streamlined code treating vector inputs as 1D.
    x = vcumop(opstr,x,dim,rev,omitnan);
elseif thisDim == dim
    n = coder.internal.indexInt(size(x,dim));
    if coder.const(rev)
        if coder.isRowMajor
            if REPLACENANS
                x = replaceNaNs(x,NVAL,thisDim + 1,varargin{:},n);
            end
            for k = n-1:-1:ONE
                x = looper(opstr,x,dim,rev,omitnan,thisDim+1,varargin{:},k);
            end
        else
            if REPLACENANS
                x = replaceNaNs(x,NVAL,thisDim - 1,n,varargin{:});
            end
            for k = n-1:-1:ONE
                x = looper(opstr,x,dim,rev,omitnan,thisDim-1,k,varargin{:});
            end
        end
    else
        if coder.isRowMajor
            if REPLACENANS
                x = replaceNaNs(x,NVAL,thisDim + 1,varargin{:},ONE);
            end
            for k = ONE:n-1
                x = looper(opstr,x,dim,rev,omitnan,thisDim+1,varargin{:},k);
            end
        else
            if REPLACENANS
                x = replaceNaNs(x,NVAL,thisDim - 1,ONE,varargin{:});
            end
            for k = ONE:n-1
                x = looper(opstr,x,dim,rev,omitnan,thisDim-1,k,varargin{:});
            end
        end
    end
elseif ONE <= thisDim && thisDim <= ND
    if coder.isRowMajor
        for k = ONE:size(x,thisDim)
            x = looper(opstr,x,dim,rev,omitnan,thisDim+1,varargin{:},k);
        end
    else
        for k = ONE:size(x,thisDim)
            x = looper(opstr,x,dim,rev,omitnan,thisDim-1,k,varargin{:});
        end
    end
else
    subsb = varargin;
    subsb{dim} = subsb{dim} + 1;
    if coder.const(rev)
        if REPLACENANS && isnan(x(varargin{:}))
            x(varargin{:}) = OP(NVAL,x(subsb{:}));
        else
            x(varargin{:}) = OP(x(varargin{:}),x(subsb{:}));            
        end
    else
        if REPLACENANS && isnan(x(subsb{:}))
            x(subsb{:}) = OP(x(varargin{:}),NVAL);
        else
            x(subsb{:}) = OP(x(varargin{:}),x(subsb{:}));            
        end
    end
end

%--------------------------------------------------------------------------

function [OP,NVAL,REPLACENANS] = constOptions(opstr,xzero,omitnan)
coder.internal.prefer_const(opstr,xzero,omitnan);
if coder.const(strcmp(opstr,'cummin'))
    if coder.const(isfloat(xzero) && ~omitnan)
        NVAL = xzero;
        OP = @minIncludeNaN;
        REPLACENANS = false;
    else
        NVAL = xzero;
        OP = @min;
        REPLACENANS = false;
    end
elseif coder.const(strcmp(opstr,'cummax'))
    if coder.const(isfloat(xzero) && ~omitnan)
        NVAL = xzero;
        OP = @maxIncludeNaN;
        REPLACENANS = false;
    else
        NVAL = xzero;
        OP = @max;
        REPLACENANS = false;
    end
elseif coder.const(strcmp(opstr,'cumprod'))
    if islogical(xzero)
        OP = @and;
    else
        OP = @times;
    end
    NVAL = cast(1,'like',xzero);
    REPLACENANS = coder.const(isfloat(xzero) && omitnan);
else % cumsum
    if islogical(xzero)
        OP = @or;
    else
        OP = @plus;
    end
    NVAL = xzero;
    REPLACENANS = coder.const(isfloat(xzero) && omitnan);
end

%--------------------------------------------------------------------------

function x = replaceNaNs(x,nval,dim,varargin)
% Replaces NaNs in the first or last element in the dimension to operate
% along with nval.
coder.inline('always');
coder.internal.prefer_const(nval,dim)
ND = coder.internal.indexInt(coder.internal.ndims(x));
ONE = coder.internal.indexInt(1);
if ONE <= dim && dim <= ND
    if coder.isRowMajor
        for k = ONE:size(x,dim)
            x = replaceNaNs(x,nval,dim + 1,varargin{:},k);
        end
    else
        for k = ONE:size(x,dim)
            x = replaceNaNs(x,nval,dim - 1,k,varargin{:});
        end
    end
else
    if isnan(x(varargin{:}))
        x(varargin{:}) = nval;
    end
end

%--------------------------------------------------------------------------

function x = vcumop(opstr,x,dim,rev,omitnan)
% Vector cumop treating x as 1D. Handles the case where dim doesn't
% correspond to the orientation of the x vector.
coder.inline('always');
coder.internal.prefer_const(opstr,dim,rev,omitnan);
nx = coder.internal.indexInt(length(x));
ONE = coder.internal.indexInt(1);
XZERO = coder.internal.scalarEg(x);
[OP,NVAL,REPLACENANS] = coder.const(@constOptions,opstr,XZERO,omitnan);
if isempty(x) || size(x,dim) == 1
    if REPLACENANS
        for k = 1:nx
            if isnan(x(k))
                x(k) = NVAL;
            end
        end
    end
elseif coder.const(rev)
    if REPLACENANS && isnan(x(nx))
        x(nx) = NVAL;
    end
    for k = nx - 1:-1:ONE
        if REPLACENANS && isnan(x(k))
            x(k) = OP(NVAL,x(k + 1));
        else
            x(k) = OP(x(k),x(k + 1));
        end
    end
else
    if REPLACENANS && isnan(x(1))
        x(1) = NVAL;
    end
    for k = ONE:nx - 1
        if REPLACENANS && isnan(x(k + 1))
            x(k + 1) = OP(x(k),NVAL);
        else
            x(k + 1) = OP(x(k),x(k + 1));
        end
    end
end

%--------------------------------------------------------------------------
