function [pks,locs,w,p] = findpeaks(y,varargin)
%FINDPEAKS Find local peaks in gpuArray data
%   PKS = FINDPEAKS(Y)
%   [PKS,LOCS] = FINDPEAKS(Y)
%   [PKS,LOCS] = FINDPEAKS(Y,X)
%   [PKS,LOCS] = FINDPEAKS(Y,Fs)
%   [PKS,LOCS,W] = FINDPEAKS(...)
%   [PKS,LOCS,W,P] = FINDPEAKS(...)
%   [...] = FINDPEAKS(...,'MinPeakHeight',MPH)
%   [...] = FINDPEAKS(...,'MinPeakProminence',MPP)
%   [...] = FINDPEAKS(...,'Threshold',TH)
%   [...] = FINDPEAKS(...,'MinPeakWidth',MINW)
%   [...] = FINDPEAKS(...,'MaxPeakWidth',MAXW)
%   [...] = FINDPEAKS(...,'MinPeakDistance',MPD)
%   [...] = FINDPEAKS(...,'SortStr',DIR)
%   [...] = FINDPEAKS(...,'NPeaks',NP)
%   FINDPEAKS(...,'Annotate',PLOTSTYLE)
%
%   See also FINDPEAKS, gpuArray.

% Copyright 2021-2023 The MathWorks, Inc.

% If Y is not on the GPU, run on CPU using the existing function.
if ~isgpuarray(y)
    [varargin{:}] = gather(varargin{:});
    if nargout > 2
        [pks,locs,w,p] = findpeaks(y,varargin{:});
    elseif nargout > 0
        [pks,locs] = findpeaks(y,varargin{:});
    else
        findpeaks(y,varargin{:});
    end
    return
end

% Extract the parameters from the input argument list.
[y,yIsRow,x,xIsRow,minH,minP,minW,maxW,minD,minT,maxN,sortDir,annotate,refW] ...
    = parse_inputs(y,varargin{:});

% Indicate if we need to compute the extent of a peak
needWidth = minW>0 || maxW<inf || minP>0 || nargout>2 || strcmp(annotate,'extents');

% Find indices of all finite and infinite peaks and the inflection points.
[iFinite,iInfinite,iInflect] = getAllPeaks(y);

% Keep only the indices of finite peaks that meet the required
% minimum height and threshold.
iPk = removePeaksBelowMinPeakHeight(y,iFinite,minH,refW);
iPk = removePeaksBelowThreshold(y,iPk,minT);

if needWidth
    % obtain the indices of each peak (iPk), the prominence base (bPk), and
    % the x- and y- coordinates of the peak base (bxPk, byPk) and the width
    % (wxPk)
    [iPk,bPk,bxPk,byPk,wxPk] = findExtents(y,x,iPk,iFinite,iInfinite,iInflect,minP,minW,maxW,refW);
else
    % combine finite and infinite peaks into one list
    [iPk,bPk,bxPk,byPk,wxPk] = combinePeaks(iPk,iInfinite);
end

% Find the indices of the largest peaks within the specified distance.
idx = findPeaksSeparatedByMoreThanMinPeakDistance(y,x,iPk,minD);

% Re-order and bound the number of peaks based upon the index vector.
idx = orderPeaks(y,iPk,idx,sortDir);
idx = keepAtMostNpPeaks(idx,maxN);

% Use the index vector to fetch the correct peaks.
iPk = iPk(idx);
if needWidth
    [bPk,bxPk,byPk,wxPk] = fetchPeakExtents(idx,bPk,bxPk,byPk,wxPk);
end

if nargout
    % Assign outputs.
    if needWidth
        [pks,locs,w,p] = assignFullOutputs(y,x,iPk,wxPk,bPk,yIsRow,xIsRow);
    else
        [pks,locs] = assignOutputs(y,x,iPk,yIsRow,xIsRow);
    end
else
    % No output arguments specified. Gather, plot and optionally annotate.
    [x,y,iPk,bPk,bxPk,byPk,wxPk,refW] = gather(x,y,iPk,bPk,bxPk,byPk,wxPk,refW);
    signal.internal.findpeaks.plot(x,y,iPk,bPk,bxPk,byPk,wxPk,refW,annotate)
end
end

%--------------------------------------------------------------------------
function [iPk,iInf,iInflect] = getAllPeaks(y)
% fetch indices all infinite peaks
iInf = find(isinf(y) & y>0);

% temporarily remove all +Inf values
yTemp = y;
yTemp(iInf) = NaN;

% determine the peaks and inflection points of the signal
[iPk,iInflect] = findLocalMaxima(yTemp);
end

%--------------------------------------------------------------------------
function [iPk, iInflect] = findLocalMaxima(yTemp)
% bookend Y by NaN and make index vector
% bookend Y by NaN and make index vector
yLeft = [NaN; yTemp]; % N+1 x 1
yRight = [yTemp; NaN]; % N+1 x 1

% Keep only the first of any adjacent pairs of equal values (including NaN).
iNeq = [true; arrayfun(@(x,y) ~isequaln(x,y), yLeft, yRight)];

% For performance keep the indices as uint32. gpuArrays are currently
% limited to 2^32 elements but we might have to go back to double or uint64
% if that changes.
iNeq = uint32(find(iNeq));

% Form the vector of values with equal sections compressed. We bookend with
% NaN so that when we get the second order diff we don't lose any values.
yNeq = [NaN;yTemp;NaN];
yNeq = yNeq(iNeq);

% Take the sign of the first sample derivative. Do this in single to
% preserve nans but allow faster processing.
s = single(sign(diff(yNeq)));

% Get the second derivative (changes of sign of derivative).
ds = [0;diff(s)];

% Find local maxima (2nd derivative is negative).
iMax = (ds<0);

% Find all transitions from rising to falling or to NaN.
iAny = ~(ds==0);

% Return indices into the original value vector without the NaN bookend.
iNeq = iNeq-1;
iInflect = iNeq(iAny);
iPk = iNeq(iMax);
end

%--------------------------------------------------------------------------
function iPk = removePeaksBelowMinPeakHeight(Y,iPk,Ph,widthRef)
if Ph>-inf && ~isempty(iPk)
    iPk = iPk(Y(iPk) > Ph);
    if isempty(iPk) && ~strcmp(widthRef,'halfheight')
        warning(message('signal:findpeaks:largeMinPeakHeight', 'MinPeakHeight', 'MinPeakHeight'));
    end
end
end

%--------------------------------------------------------------------------
function iPk = removePeaksBelowThreshold(Y,iPk,Th)
% Measure the peak against its neighbouring elements and remove those below
% the threshold.
if Th<=0
    % By definition peaks cannot be below their neighbors!
    return;
end
base = max(Y(iPk-1),Y(iPk+1));
iPk = iPk(Y(iPk)-base >= Th);
end


%--------------------------------------------------------------------------
function idx = findPeaksSeparatedByMoreThanMinPeakDistance(y,x,iPk,Pd)
% Start with the larger peaks to make sure we don't accidentally keep a
% small peak and remove a large peak in its neighborhood.

if isempty(iPk) || Pd==0
    idx = cast((1:length(iPk)).', "like", iPk);
    return
end

% copy peak values and locations to a temporary place
pks = y(iPk);
locs = x(iPk);

% Order peaks from large to small. Note that we gather the location indices
% since the search operation is faster on CPU.
[~, sortIdx] = sort(pks,'descend');
locs_temp = gather(locs(sortIdx));

idelete = true(size(locs_temp))<0;
for ii = 1:length(locs_temp)
    if ~idelete(ii)
        % If the peak is not in the neighborhood of a larger peak, find
        % secondary peaks to eliminate.
        idelete = idelete | (locs_temp>=locs_temp(ii)-Pd) & (locs_temp<=locs_temp(ii)+Pd);
        idelete(ii) = false; % Keep current peak
    end
end

% report back indices in consecutive order
idx = sort(sortIdx(~idelete));
end

%--------------------------------------------------------------------------
function [iPk,bPk,bxPk,byPk,wxPk] = combinePeaks(iPk,iInf)
% Combine normal peaks and infinite peaks into a single sorted list.

% Avoid calling UNION if iInf is empty since it will re-sort iPk!
if ~isempty(iInf)
    iPk = union(iPk,iInf);
end
bPk = zeros(0,1,"like",iPk);
bxPk = zeros(0,2,"like",iPk);
byPk = zeros(0,2,"like",iPk);
wxPk = zeros(0,2,"like",iPk);
end

%--------------------------------------------------------------------------
function idx = orderPeaks(Y,iPk,idx,Str)
% Order the peaks in height (Y) order
if isempty(idx) || strcmp(Str,'none')
    return
end

[~,s]  = sort(Y(iPk(idx)),Str);
idx = idx(s);
end

%--------------------------------------------------------------------------
function idx = keepAtMostNpPeaks(idx,Np)

if length(idx)>Np
    idx = idx(matlab.internal.ColonDescriptor(1,1,Np));
end
end

%--------------------------------------------------------------------------
function [bPk,bxPk,byPk,wxPk] = fetchPeakExtents(idx,bPk,bxPk,byPk,wxPk)
bPk = bPk(idx);
bxPk = bxPk(idx,:);
byPk = byPk(idx,:);
wxPk = wxPk(idx,:);
end


%--------------------------------------------------------------------------
function [iPk,bPk,bxPk,byPk,wxPk] = findExtents(y,x,iPk,iFin,iInf,iInflect,minP,minW,maxW,refW)
% findExtents is a highly iterative algorithm that loops over elements near
% each peak. It is therefore not very suitable for running on the GPU and we
% gather and call the CPU version to avoid dramatic slow-downs.

% Only the data args can be on the GPU so gather them.
[y,x,iPk,iFin,iInf,iInflect] = gather(y,x,iPk,iFin,iInf,iInflect);

% Call the CPU version.
[iPk,bPk,bxPk,byPk,wxPk] = signal.internal.findpeaks.findExtents( ...
    y,x,iPk,iFin,iInf,iInflect,minP,minW,maxW,refW);

% Put the data args back on GPU. values that depend on X only go on the GPU
% if X is numeric (i.e. not a datetime, duration, or similar).
iPk = gpuArray(iPk);
bPk = gpuArray(bPk);
byPk = gpuArray(byPk);
if isnumeric(x)
    bxPk = gpuArray(bxPk);
    wxPk = gpuArray(wxPk);
end

end

%--------------------------------------------------------------------------
function [Ypk,Xpk] = assignOutputs(y,x,iPk,yIsRow,xIsRow)

% fetch the coordinates of the peak
Ypk = y(iPk);
Xpk = x(iPk);

% preserve orientation of Y
if yIsRow
    Ypk = reshape(Ypk,1,[]);
end

% preserve orientation of X
if xIsRow
    Xpk = reshape(Xpk,1,[]);
end

end

%--------------------------------------------------------------------------
function [Ypk,Xpk,Wpk,Ppk] = assignFullOutputs(y,x,iPk,wxPk,bPk,yIsRow,xIsRow)

% fetch the coordinates of the peak
Ypk = y(iPk);
Xpk = x(iPk);

% compute the width and prominence
Wpk = diff(wxPk,1,2);
Ppk = Ypk-bPk;

% preserve orientation of Y (and P)
if yIsRow
    Ypk = reshape(Ypk,1,[]);
    Ppk = reshape(Ppk,1,[]);
end

% preserve orientation of X (and W)
if xIsRow
    Xpk = reshape(Xpk,1,[]);
    Wpk = reshape(Wpk,1,[]);
end
end
%--------------------------------------------------------------------------
function [y,yIsRow,x,xIsRow,Ph,Pp,Wmin,Wmax,Pd,Th,Np,Str,Ann,Ref] = parse_inputs(Yin,varargin)

% Validate input signal
validateattributes(Yin,{'double','single'},{'nonempty','real','vector'},...
    'findpeaks','Y');
y = reshape(Yin,[],1);
M = length(y);

if M < 3
    error(message('signal:findpeaks:emptyDataSet'));
end
yIsRow = isrow(Yin);

% indicate if the user specified an Fs or X
hasX = ~isempty(varargin) && (isnumeric(varargin{1}) || ...
    isdatetime(varargin{1}) && length(varargin{1})>1);

if hasX
    startArg = 2;
    FsSupplied = isscalar(varargin{1});
    if FsSupplied
        % Fs
        Fs = varargin{1};
        validateattributes(Fs,{'numeric'},{'real','finite','positive'},'findpeaks','Fs');
        Fs = double(Fs(1));
        x = gpuArray.colon(0,M-1).'/Fs;
        xRange = M-1./Fs;
        xIsRow = yIsRow;
    else
        % X
        Xin = varargin{1};

        if isnumeric(Xin)
            validateattributes(Xin,{'numeric'},{'real','finite','vector','increasing'},'findpeaks','X');
            xRange = gather(double(Xin(M) - Xin(1)));
            % Convert to gpuArray of double if not already.
            Xin = gpuArray(double(Xin));
        else % isdatetime(Xin)
            validateattributes(seconds(Xin-Xin(1)),{'double'},{'real','finite','vector','increasing'},'findpeaks','X');
        end

        if length(Xin) ~= M
            throwAsCaller(MException(message('signal:findpeaks:mismatchYX')));
        end
        xIsRow = isrow(Xin);
        x = reshape(Xin,[],1);
    end
else
    startArg = 1;
    % unspecified, use index vector
    x = gpuArray.colon(1,M).';
    xRange = M-1;
    xIsRow = yIsRow;
end

% All params after Y, X must be in-memory
if numel(varargin) >= startArg
    [varargin{startArg:end}] = gather(varargin{startArg:end});
end

% Parse the inputs
parser = signal.internal.findpeaks.getParser();
parse(parser,varargin{startArg:end});
Ph = parser.Results.MinPeakHeight;
Pp = parser.Results.MinPeakProminence;
Wmin = parser.Results.MinPeakWidth;
Wmax = parser.Results.MaxPeakWidth;
Pd = parser.Results.MinPeakDistance;
Th = parser.Results.Threshold;
Np = parser.Results.NPeaks;
Str = parser.Results.SortStr;
Ann = parser.Results.Annotate;
Ref = parser.Results.WidthReference;

% Make sure we stop checking as soon as we can
numToCheck = numel(fieldnames(parser.Results)) - numel(parser.UsingDefaults);
userSpecified = iIsUserDefined(parser);

% Limit the number of peaks to the number of input samples if not supplied.
if isempty(Np)
    Np = M;
else
    validateattributes(Np,{'numeric'},{'real','scalar','nonempty','integer','positive'},'findpeaks','NPeaks');
    numToCheck = numToCheck - 1;
end

% If no other options were supplied we can stop now and use the defaults.
if numToCheck == 0
    return
end

% Only validate inputs passed explicitly by the user.
if userSpecified.MinPeakHeight
    validateattributes(Ph,{'numeric'},{'real','scalar','nonempty'},'findpeaks','MinPeakHeight');
end
if userSpecified.MinPeakDistance
    if isnumeric(x)
        validateattributes(Pd,{'numeric'},{'real','scalar','nonempty','nonnegative','<',xRange},'findpeaks','MinPeakDistance');
    else
        if isduration(Pd)
            validateattributes(seconds(Pd),{'numeric'},{'real','scalar','nonempty','nonnegative'},'findpeaks','MinPeakDistance');
        else
            validateattributes(Pd,{'numeric'},{'real','scalar','nonempty','nonnegative'},'findpeaks','MinPeakDistance');
        end
    end
end
if userSpecified.MinPeakProminence
    validateattributes(Pp,{'numeric'},{'real','scalar','nonempty','nonnegative'},'findpeaks','MinPeakProminence');
end
if userSpecified.MinPeakWidth
    if isduration(Wmin)
        validateattributes(seconds(Wmin),{'numeric'},{'real','scalar','finite','nonempty','nonnegative'},'findpeaks','MinPeakWidth');
    else
        validateattributes(Wmin,{'numeric'},{'real','scalar','finite','nonempty','nonnegative'},'findpeaks','MinPeakWidth');
    end
end
if userSpecified.MaxPeakWidth
    if isduration(Wmax)
        validateattributes(seconds(Wmax),{'numeric'},{'real','scalar','nonnan','nonempty','nonnegative'},'findpeaks','MaxPeakWidth');
    else
        validateattributes(Wmax,{'numeric'},{'real','scalar','nonnan','nonempty','nonnegative'},'findpeaks','MaxPeakWidth');
    end
end
if userSpecified.Threshold
    validateattributes(Th,{'numeric'},{'real','scalar','nonempty','nonnegative'},'findpeaks','Threshold');
end
if userSpecified.SortStr
    Str = validatestring(Str,{'ascend','none','descend'},'findpeaks','SortStr');
end
if userSpecified.Annotate
    Ann = validatestring(Ann,{'peaks','extents'},'findpeaks','SortStr');
end
if userSpecified.WidthReference
    Ref = validatestring(Ref,{'halfprom','halfheight'},'findpeaks','WidthReference');
    % ignore peaks below zero when using halfheight width reference
    if strcmp(Ref,'halfheight')
        Ph = max(Ph,0);
    end
end
end

%--------------------------------------------------------------------------
% GPU-specific helpers
%--------------------------------------------------------------------------
function s = iIsUserDefined(parser)
% Return a struct indicating which parameters were specified by the user.
% The struct has one field per parameter with true indicating
% user-specified and false indicating default value.
allParams = fieldnames(parser.Results);
userSpecified = ~ismember(allParams, parser.UsingDefaults);
s = cell2struct(num2cell(userSpecified), allParams, 1);
end
