function [LPY, LPI, codistr] = hMinMaxImpl(codistr, fcnMinMax, LP, dim, wantI, linearI)
%hMinMaxImpl Implementation for codistributor2dbc.

%   Copyright 2009-2019 The MathWorks, Inc.

% Should not get here if dim>2 or all dimensions being reduced are trivial
assert( all(dim<=2) && all(codistr.Cached.GlobalSize(dim) > 1), ...
    "Non-trivial MIN/MAX called for trivial reduction." );
orgCodistr = codistr;

% Min and max on a non-empty array lead to size in dimension being 1.
gsize = codistr.Cached.GlobalSize;
gsize(dim) = 1;

% Construct local parts for the workers that now have empty local parts.
codistr = codistributor2dbc(codistr.WorkerGrid, ...
                            codistr.BlockSize, ...
                            codistr.Orientation, gsize);

if wantI
    [LPY, LPI] = fcnMinMax(LP, [], dim);
    % Convert local indices into global indices
    if linearI
        LPI = orgCodistr.hLocalToGlobalLinearIndices(LPI, spmdIndex);
    else
        % One dimension only, so convert using globalIndices
        LPI = orgCodistr.hLocalToGlobalIndices(LPI, dim, spmdIndex);
    end
    % Make sure LPI and LPY are the same shape
    LPI = reshape(LPI, size(LPY));
else
    LPY = fcnMinMax(LP, [], dim);
    LPI = [];
end

if all(codistr.WorkerGrid(dim) == 1)
    % Each row/column resides in its entirety on a single worker, and all
    % workers are on the edge of the worker grid.
    return;
    
elseif all(ismember([1 2], dim))
    % Reducing in both dimensions, so result is scalar
    [LPY, LPI] = iReduceToScalar(codistr, fcnMinMax, LPY, LPI, dim, wantI, linearI);

elseif orgCodistr.Cached.GlobalSize(dim) <= codistr.BlockSize
    % We do not need to communicate because all the data already resides on the edge
    % of the worker grid.
    [LPY, LPI] = iHandleInteriorWorkers(codistr, LPY, LPI, dim, wantI);
    
else

    [LPY, LPI] = iMinMaxAcrossWorkers(codistr, fcnMinMax, LPY, LPI, dim, wantI, linearI);
end
end % End of hMinMaxImpl.

function [LPY, LPI] = iReduceToScalar(codistr, fcnMinMax, LPY, LPI, dim, wantI, linearI)
% We are taking MIN/MAX over both dimensions so all input local parts will
% be scalar or empty and result is scalar on worker 1.

if wantI
    % Reduce both values and indices. We avoid duplicate communication by
    % combining into a struct
    data = struct("Y", LPY, "I", LPI);
    assert(linearI, "Only linear indices can be returned when reducing in multiple dimensions");
    % Concatenate results in both dimensions rather than gathering because
    % we do not know the intermediate data size (required by gather).
    data = codistr.pCatToEdgeOfWorkerGrid(data, [1 2]);
    if spmdIndex == 1
        % To ensure we resolve ties to the lowest global index, put the local
        % results into global index order. We can treat as a vector since the
        % result is scalar anyway.
        LPY = cat(2,data.Y);
        LPI = cat(2,data.I);
        [LPI,idx] = sort(LPI(:));
        LPY =  LPY(idx);
    
        % Take the min/max and then find the corresponding index entries
        [LPY, idx] = fcnMinMax(LPY, [], dim);
        LPI = LPI(idx);
    else
        LPY = distributedutil.Allocator.create(codistr.hLocalSize(), LPY);
        LPI = distributedutil.Allocator.create(codistr.hLocalSize(), LPI);
    end
else
    % Value only, so order does not matter
    % Concatenate results in both dimensions rather than gathering because
    % we do not know the intermediate data size (required by gather).
    data = codistr.pCatToEdgeOfWorkerGrid(LPY, [1 2]);
    if spmdIndex == 1
        LPY = fcnMinMax(data, [], dim); 
    else
        LPY = distributedutil.Allocator.create(codistr.hLocalSize(), LPY);
    end
    LPI = [];    
end

end

function [LPY, LPI] = iMinMaxAcrossWorkers(codistr, fcnMinMax, LPY, LPI, dim, wantI, linearI)
% The input LPY and LPI are the data and the indices after calculating min/max
% on the local parts.
% We perform the second stage of min/max where we operate across the workers.
assert(numel(dim)==1, "General MIN/MAX called with more than one dimension");
if wantI
    % Get LPY and LPI into a cell array so that we only communicate once.  Odd
    % entries store the Y's, even entries store the indices.    
    data = cat(dim, {LPY}, {LPI});
    data = codistr.pCatToEdgeOfWorkerGrid(data, dim);
    if codistr.pIsOnWorkerGridEdge(dim)
        firstData = cat(dim, data{1:2:end-1});
        firstInd = cat(dim, data{2:2:end});
        % We want the min/max operation to be stable.  I.e. if there are ties, we want
        % to pick the entry with the lowest global index in the dim-th
        % dimension.  We achieve this by re-arranging firstData and firstInd so
        % that the global indices are in ascending order in the dim-th
        % dimension.  
        [firstInd, rearrangement] = sort(firstInd, dim);
        linInd = distributedutil.IndexManip.sortIndToLinear(dim, rearrangement);
        firstData = firstData(linInd);
 
        % Perform the second stage min/max across all the local parts.
        [LPY, indToFirst] = fcnMinMax(firstData,[], dim); 
        % Ensure the resulting index is linear
        if ~linearI
            indToFirst = distributedutil.IndexManip.minmaxIndToLinear( ...
                size(firstData), dim, indToFirst);
        end
        % Now firstData corresponds to firstInd, and LPY equals
        % firstData(linIndIntoFirst).  Therefore, Y corresponds to
        % firstInd(linIndIntoFirst).
        LPI = firstInd(indToFirst); 
    else
        % These workers don't store any data, so just allocate empty matrices of the
        % correct size and type.
        LPY = distributedutil.Allocator.create(codistr.hLocalSize(), LPY);
        LPI = distributedutil.Allocator.create(codistr.hLocalSize(), LPI);
    end
else
    data = codistr.pCatToEdgeOfWorkerGrid(LPY, dim);
    if codistr.pIsOnWorkerGridEdge(dim)
        LPY = fcnMinMax(data, [], dim); 
    else
        LPY = distributedutil.Allocator.create(codistr.hLocalSize(), LPY);
    end
    LPI = [];
end
end % End of iMinMaxAcrossWorkers.

function [LPY, LPI] = iHandleInteriorWorkers(codistr, LPY, LPI, dim, wantI)
% The input LP and codistr are the local part and codistributor after min/max
% along dim.  Handle the special case where the local part was non-trivial only
% on the edge of the worker grid, so we don't need to communicate to complete the
% min/max calculations.
if ~any(codistr.pIsOnWorkerGridEdge(dim))
    LPY = distributedutil.Allocator.create(codistr.hLocalSize(), LPY);
    if wantI
        LPI = distributedutil.Allocator.create(codistr.hLocalSize(), LPI);
    end
end
end % End of iHandleInteriorWorkers.
