function [y,counts] = blockedSummation(b,x,dim,omitnan,identity,vlen,vstart)
%MATLAB Code Generation Library Function

%   Performs blocked summation along dimension DIM of x.
%
%   B is the block size.
%
%   OMITNAN defaults to false;
%
%   DIM defaults to coder.internal.constNonSingletonDim(x).
%
%   The IDENTITY value determines both the output type and the accumulator
%   type. Note that IDENTITY must be complex if the output of F is complex.
%
%   VLEN is an optional input that defaults to size(X,DIM). VLEN must be
%   in the range 0 <= VLEN <= size(X,DIM). Only the first VLEN elements of
%   each vector are combined to form the output Y.
%
%   VSTART is an optional input that defaults to 1. Each vector processed
%   starts at index VSTART and ends at index VSTART + VLEN - 1.
%
%   The second output, COUNTS, is an array the same size as Y. COUNTS(k)
%   is the number of elements that contributed to the value of Y(k).

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

coder.internal.allowEnumInputs;
coder.internal.allowHalfInputs;
if nargin < 3
    dim = coder.internal.constNonSingletonDim(x);
end
if nargin < 4
    omitnan = false;
end
if nargin < 5
    % scalarEg(logical,logical) --> double
    % scalarEg(xchar,xchar) --> double
    identity = coder.internal.scalarEg(x,x);
end
if nargin < 6
    vlen = coder.internal.indexInt(size(x,dim));
end
if nargin < 7
    vstart = coder.internal.indexInt(1);
end
coder.internal.prefer_const(b,dim,omitnan,identity,vlen,vstart);
coder.internal.assert(coder.internal.isConst(dim), ...
    'Coder:toolbox:dimNotConst');
% Avoid counting during iterations unless counts are requested and there is
% a chance they will turn out to be anything but vlen.
DOCOUNTS = nargout == 2 && isfloat(x) && omitnan;
FLATTEN = ~coder.internal.avoidArrayFlattening;
IS2D = coder.const(coder.internal.ndims(x) == 2);
INCLUDENAN = coder.const(coder.internal.isConstTrue(~omitnan));
ND = coder.internal.indexInt(coder.internal.ndims(x));
SIM = coder.const(coder.target('Sfun') || coder.target('MEX'));
STRIDE1_ROWS = coder.const(dim == 2 && coder.isRowMajor);
STRIDE1_COLS = coder.const(dim == 1 && coder.isColumnMajor);
STRIDE1_SIMD_LOOPS = coder.const(SIM && INCLUDENAN && ( ...
    (STRIDE1_ROWS && IS2D) || ...
    (STRIDE1_COLS && (IS2D || FLATTEN))));
TRANSROW = coder.const(coder.isColumnMajor && SIM && INCLUDENAN && ...
    dim == 2 && coder.internal.isConstTrue(isrow(x)));
TRANSCOL = coder.const(coder.isRowMajor && SIM && INCLUDENAN && ...
    dim == 1 && coder.internal.isConstTrue(iscolumn(x)));
ONE = coder.internal.indexInt(1);
TWO = coder.internal.indexInt(2);
% Perform the operation.
if isempty(x) || vlen == 0
    % If size(x,dim) is zero, y can be non-empty when x is empty. It's
    % cleanest and safest to xdeal with this possibility as a special case.
    if dim >= 1 && dim <= ND
        sz = size(x);
        sz(dim) = 1;
        y = eml_expand(identity,sz);
    else
        y = eml_expand(identity,size(x));
    end
    if DOCOUNTS
        counts = zeros(size(y),'like',vlen);
    else
        counts = coder.nullcopy(zeros(size(y),'like',vlen));
    end
elseif dim > ND
    if DOCOUNTS
        [y,counts] = castLikeAndReplaceNaNs(x,identity);
    elseif omitnan
        y = castLikeAndReplaceNaNs(x,identity);
        counts = coder.nullcopy(zeros(size(y),'like',vlen));
    else
        y = cast(x,'like',identity);
        counts = coder.nullcopy(zeros(size(y),'like',vlen));
    end
elseif TRANSROW
    % Since the input is a row vector, we inexpensively transpose it and
    % then use the optimized code for summing down columns in column-major
    % mode.
    y = sumMatrixIncludeNaN(b,identity,x(:),ONE,vlen,vstart);
elseif TRANSCOL
    % Since the input is a column vector, we can inexpensively transpose it
    % and then use the optimized code for summing across rows in row-major
    % mode.
    y = sumMatrixIncludeNaN(b,identity,x(:).',TWO,vlen,vstart);
elseif STRIDE1_SIMD_LOOPS
    % Faster code for stride = 1 sums.
    y = sumMatrixIncludeNaN(b,identity,x,dim,vlen,vstart);
elseif FLATTEN
    if DOCOUNTS
        [y,counts] = colMajorFlatIter(b,identity,x,dim,omitnan,vlen,vstart);
    else
        y = colMajorFlatIter(b,identity,x,dim,omitnan,vlen,vstart);
        counts = coder.nullcopy(zeros(size(y),'like',vlen));
    end
else
    if DOCOUNTS
        [y,counts] = nestedIter(b,identity,x,dim,omitnan,vlen,vstart);
    else
        y = nestedIter(b,identity,x,dim,omitnan,vlen,vstart);
        counts = coder.nullcopy(zeros(size(y),'like',vlen));
    end
end
if nargout == 2 && ~DOCOUNTS
    % We did not do counts during iterations, but the output is requested.
    counts = eml_expand(vlen,size(y));
end

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

function [y,counts] = colMajorFlatIter(b,identity,x,dim,omitnan,vlen,vstart)
coder.internal.prefer_const(b,identity,dim,omitnan,vlen,vstart);
ZERO = coder.internal.indexInt(0);
ONE = coder.internal.indexInt(1);
TWO = coder.internal.indexInt(2);
OMITNANS = isfloat(x) && omitnan;
DOCOUNTS = nargout == 2;
vstride = coder.internal.prodsize(x,'below',dim);
bvstride = b*vstride;
pagesize = vstride*coder.internal.indexInt(size(x,dim));
npages = coder.internal.prodsize(x,'above',dim);
voffset = (vstart - 1)*vstride;
% Preallocate outputs.
if coder.const(dim >= 1 && dim <= coder.internal.ndims(x))
    sz = size(x);
    sz(dim) = 1;
    y = coder.nullcopy(eml_expand(identity,sz));
else
    y = coder.nullcopy(eml_expand(identity,size(x)));
end
if DOCOUNTS
    counts = coder.nullcopy(zeros(size(y),coder.internal.indexIntClass));
end
bsum = coder.nullcopy(eml_expand(identity,[vstride,1]));
if vlen <= b
    firstBlockLength = vlen;
    lastBlockLength = ZERO;
    nblocks = ONE;
else
    firstBlockLength = b;
    nblocks = coder.internal.indexDivide(vlen,b);
    lastBlockLength = vlen - nblocks*b;
    if lastBlockLength > ZERO
        nblocks = nblocks + 1;
    else
        lastBlockLength = b;
    end
end
for xi = ONE:npages
    xpageoffset = (xi - 1)*pagesize + voffset;
    ypageoffset = (xi - 1)*vstride;
    % Initialize the current page of y with the first elements of
    % each vector to be processed.
    for xj = ONE:vstride
        ix = xpageoffset + xj; % index of the first element of x vector
        iy = ypageoffset + xj; % index of the output value
        if OMITNANS && isnan(x(ix))
            y(iy) = identity;
            if DOCOUNTS
                counts(iy) = 0;
            end
        else
            y(iy) = cast(x(ix),'like',identity);
            if DOCOUNTS
                counts(iy) = 1;
            end
        end
        bsum(xj) = 0;
    end
    % Accumulate the first block directly into Y.
    for k = TWO:firstBlockLength
        xoffset = xpageoffset + (k - 1)*vstride;
        for xj = ONE:vstride
            ix = xoffset + xj;
            iy = ypageoffset + xj;
            if ~(OMITNANS && isnan(x(ix)))
                y(iy) = y(iy) + cast(x(ix),'like',identity);
                if DOCOUNTS
                    counts(iy) = counts(iy) + 1;
                end
            end
        end
    end
    % Accumulate subsequent blocks into BSUM and then add to Y.
    for ib = TWO:nblocks
        xblockoffset = xpageoffset + (ib - 1)*bvstride;
        % Initialize BSUM.
        for xj = ONE:vstride
            ix = xblockoffset + xj;
            if OMITNANS && isnan(x(ix))
                bsum(xj) = identity;
            else
                bsum(xj) = cast(x(ix),'like',identity);
                if DOCOUNTS
                    iy = ypageoffset + xj;
                    counts(iy) = counts(iy) + 1;
                end
            end
        end
        % Add the remaining elements of the block.
        if ib == nblocks
            hi = lastBlockLength;
        else
            hi = b;
        end
        for k = TWO:hi
            xoffset = xblockoffset + (k - 1)*vstride;
            for xj = ONE:vstride
                ix = xoffset + xj;
                if ~(OMITNANS && isnan(x(ix)))
                    bsum(xj) = bsum(xj) + cast(x(ix),'like',identity);
                    if DOCOUNTS
                        iy = ypageoffset + xj;
                        counts(iy) = counts(iy) + 1;
                    end
                end
            end
        end
        for xj = ONE:vstride
            iy = ypageoffset + xj;
            y(iy) = y(iy) + bsum(xj);
        end
    end
end

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

function [y,counts] = nestedIter(b,identity,x,dim,omitnan,vlen,vstart)
% Assumes 1 <= dim <= length(size(x)).
coder.internal.prefer_const(b,identity,dim,omitnan,vlen,vstart);
ONE = coder.internal.indexInt(1);
DOCOUNTS = nargout == 2;
ND = coder.internal.indexInt(coder.internal.ndims(x));
sz = size(x);
sz(dim) = 1; % Would error here if dim were out of range.
y = coder.nullcopy(eml_expand(identity,sz));
szbsum = size(x);
if coder.isRowMajor
    for k = coder.unroll(1:dim)
        szbsum(k) = 1;
    end
else
    for k = coder.unroll(dim:ND)
        szbsum(k) = 1;
    end
end
bsum = coder.nullcopy(eml_expand(identity,szbsum));
if DOCOUNTS
    counts = coder.nullcopy(zeros(size(y),coder.internal.indexIntClass));
else
    counts = coder.internal.indexInt([]);
end
if coder.isRowMajor
    OuterDim = ONE;
else
    OuterDim = ND;
end
[y,bsum,counts] = genloops(b,identity,x,dim,omitnan,vlen,vstart, ...
    y,bsum,counts,OuterDim); %#ok<ASGLU>

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

function [y,bsum,counts] = genloops(b,identity,x,dim,omitnan,vlen,vstart, ...
    y,bsum,counts,loopdim,varargin)
% Iterate using xfull ND-indexing.
coder.inline('always');
coder.internal.prefer_const(b,identity,dim,omitnan,vlen,vstart,loopdim);
ROWMAJOR = coder.isRowMajor;
OMITNANS = isfloat(x) && omitnan;
EMPTYCOUNTS = coder.internal.isConstTrue(isempty(counts));
DOCOUNTS = nargout == 3 && ~EMPTYCOUNTS;
ZERO = coder.internal.indexInt(0);
ONE = coder.internal.indexInt(1);
TWO = coder.internal.indexInt(2);
ND = coder.internal.indexInt(coder.internal.ndims(x));
voffset = vstart - 1;
if coder.const(loopdim == dim)
    if vlen <= b
        firstBlockLength = vlen;
        lastBlockLength = ZERO;
        nblocks = ONE;
    else
        firstBlockLength = b;
        nblocks = coder.internal.indexDivide(vlen,b);
        lastBlockLength = vlen - nblocks*b;
        if lastBlockLength > ZERO
            nblocks = nblocks + 1;
        else
            lastBlockLength = b;
        end
    end
    % Process the first xmin(vlen,b) elements.
    if ROWMAJOR
        nextdim = loopdim + 1;
        [y,counts] = initPage(false,identity,x,dim,omitnan, ...
            y,counts,nextdim,varargin{:},vstart);
    else
        nextdim = loopdim - 1;
        [y,counts] = initPage(false,identity,x,dim,omitnan, ...
            y,counts,nextdim,vstart,varargin{:});
    end
    for k = TWO:firstBlockLength
        if ROWMAJOR
            [y,bsum,counts] = genloops(b,identity,x,dim,omitnan, ...
                vlen,vstart,y,bsum,counts,nextdim,varargin{:},voffset + k);
        else
            [y,bsum,counts] = genloops(b,identity,x,dim,omitnan, ...
                vlen,vstart,y,bsum,counts,nextdim,voffset + k,varargin{:});
        end
    end
    % Add remaining blocks.
    for ib = TWO:nblocks
        offset = (ib - 1)*b + voffset;
        if ROWMAJOR
            [bsum,counts] = initPage(true,identity,x,dim,omitnan, ...
                bsum,counts,nextdim,varargin{:},offset + ONE);
        else
            [bsum,counts] = initPage(true,identity,x,dim,omitnan, ...
                bsum,counts,nextdim,offset + ONE,varargin{:});
        end
        if ib == nblocks
            hi = lastBlockLength;
        else
            hi = b;
        end
        for k = TWO:hi
            if ROWMAJOR
                [bsum,counts] = genBlockSumLoops(b,identity,x,dim,omitnan, ...
                    bsum,counts,vlen,nextdim,varargin{:},offset + k);
            else
                [bsum,counts] = genBlockSumLoops(b,identity,x,dim,omitnan, ...
                    bsum,counts,vlen,nextdim,offset + k,varargin{:});
            end
        end
        if ROWMAJOR
            y = addBlockSum(y,bsum,dim,nextdim,varargin{:},ONE);
        else
            y = addBlockSum(y,bsum,dim,nextdim,ONE,varargin{:});
        end
    end
elseif coder.const(loopdim >= ONE && loopdim <= ND)
    for k = ONE:size(x,loopdim)
        if ROWMAJOR
            [y,bsum,counts] = genloops(b,identity,x,dim,omitnan, ...
                vlen,vstart,y,bsum,counts,loopdim + 1,varargin{:},k);
        else
            [y,bsum,counts] = genloops(b,identity,x,dim,omitnan, ...
                vlen,vstart,y,bsum,counts,loopdim - 1,k,varargin{:});
        end
    end
elseif vlen >= TWO
    % At this point we are in the innermost loop and need to add the
    % currently-indexed element of x to its accumulator in y.
    xk = x(varargin{:});
    if ~OMITNANS || ~isnan(xk)
        ysubs = varargin;
        ysubs{dim} = ONE;
        y(ysubs{:}) = y(ysubs{:}) + cast(xk,'like',identity);
        if DOCOUNTS
            csubs = varargin;
            csubs{dim} = ONE;
            counts(csubs{:}) = counts(csubs{:}) + 1;
        end
    end
end

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

function [bsum,counts] = genBlockSumLoops(b,identity,x,dim,omitnan,bsum, ...
    counts,vlen,loopdim,varargin)
% Iterate using xfull ND-indexing.
coder.inline('always');
coder.internal.prefer_const(b,identity,dim,omitnan,vlen,loopdim);
ROWMAJOR = coder.isRowMajor;
OMITNANS = isfloat(x) && omitnan;
EMPTYCOUNTS = coder.internal.isConstTrue(isempty(counts));
DOCOUNTS = nargout == 2 && ~EMPTYCOUNTS;
ONE = coder.internal.indexInt(1);
TWO = coder.internal.indexInt(2);
ND = coder.internal.indexInt(coder.internal.ndims(x));
if coder.const(loopdim >= ONE && loopdim <= ND)
    for k = ONE:size(x,loopdim)
        if ROWMAJOR
            [bsum,counts] = genBlockSumLoops(b,identity,x,dim,omitnan, ...
                bsum,counts,vlen,loopdim + 1,varargin{:},k);
        else
            [bsum,counts] = genBlockSumLoops(b,identity,x,dim,omitnan, ...
                bsum,counts,vlen,loopdim - 1,k,varargin{:});
        end
    end
elseif vlen >= TWO
    % At this point we are in the innermost loop and need to add the
    % currently-indexed element of x to its accumulator in y.
    xk = x(varargin{:});
    if ~OMITNANS || ~isnan(xk)
        bsubs = xSubsToBlockSubs(dim,varargin);
        bsum(bsubs{:}) = bsum(bsubs{:}) + cast(xk,'like',identity);
        if DOCOUNTS
            csubs = varargin;
            csubs{dim} = ONE;
            counts(csubs{:}) = counts(csubs{:}) + 1;
        end
    end
end

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

function [y,counts] = initPage(isblock,identity,x,dim,omitnan, ...
    y,counts,loopdim,varargin)
coder.inline('always');
coder.internal.prefer_const(isblock,identity,dim,omitnan,loopdim);
ROWMAJOR = coder.isRowMajor;
OMITNANS = isfloat(x) && omitnan;
EMPTYCOUNTS = coder.internal.isConstTrue(isempty(counts));
DOCOUNTS = nargout == 2 && ~EMPTYCOUNTS;
ONE = coder.internal.indexInt(1);
ND = coder.internal.indexInt(coder.internal.ndims(x));
if coder.const(loopdim >= ONE && loopdim <= ND)
    for k = ONE:size(x,loopdim)
        if ROWMAJOR
            [y,counts] = initPage(isblock,identity,x,dim,omitnan, ...
                y,counts,loopdim + 1,varargin{:},k);
        else
            [y,counts] = initPage(isblock,identity,x,dim,omitnan, ...
                y,counts,loopdim - 1,k,varargin{:});
        end
    end
elseif isblock
    % Y is just a block sum page. It must be indexed differently, and
    % COUNTS should not be initialized, only incremented, if appropriate.
    xsubs = varargin;
    xk = x(xsubs{:});
    bsubs = xSubsToBlockSubs(dim,xsubs);
    if OMITNANS && isnan(xk)
        y(bsubs{:}) = identity;
    else
        y(bsubs{:}) = xk;
        if DOCOUNTS
            countsubs = varargin;
            countsubs{dim} = ONE;
            counts(countsubs{:}) = counts(countsubs{:}) + 1;
        end
    end
else
    xk = x(varargin{:});
    ysubs = varargin;
    ysubs{dim} = ONE;
    if OMITNANS && isnan(xk)
        y(ysubs{:}) = identity;
        if DOCOUNTS
            counts(ysubs{:}) = 0;
        end
    else
        y(ysubs{:}) = xk;
        if DOCOUNTS
            counts(ysubs{:}) = 1;
        end
    end
end

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

function y = addBlockSum(y,bsum,dim,loopdim,varargin)
coder.inline('always');
coder.internal.prefer_const(dim,loopdim);
ONE = coder.internal.indexInt(1);
ND = coder.internal.indexInt(coder.internal.ndims(y));
ROWMAJOR = coder.isRowMajor;
if coder.const(loopdim >= ONE && loopdim <= ND)
    for k = ONE:size(y,loopdim)
        if ROWMAJOR
            y = addBlockSum(y,bsum,dim,loopdim + 1,varargin{:},k);
        else
            y = addBlockSum(y,bsum,dim,loopdim - 1,k,varargin{:});
        end
    end
else
    bsubs = xSubsToBlockSubs(dim,varargin);
    ysubs = varargin;
    ysubs{dim} = ONE;
    y(ysubs{:}) = y(ysubs{:}) + bsum(bsubs{:});
end

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

function bsubs = xSubsToBlockSubs(dim,xsubs)
coder.inline('always');
ONE = coder.internal.indexInt(1);
bsubs = xsubs;
coder.internal.prefer_const(dim);
ND = coder.internal.indexInt(length(xsubs));
if coder.isRowMajor
    for k = coder.unroll(1:dim)
        bsubs{k} = ONE;
    end
else
    for k = coder.unroll(dim:ND)
        bsubs{k} = ONE;
    end
end

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