function y = prctile(x,p,varargin)
%PRCTILE Percentiles of a sample.
%   Y = PRCTILE(X,P) returns percentiles of the values in X. P is a scalar
%   or a vector of percentage values for which to compute percentiles. When
%   X is a vector, Y is the same size as P, and Y(i) contains the P(i)-th
%   percentile. When X is a matrix, the i-th row of Y contains the P(i)-th
%   percentiles of each column of X. For N-D arrays, PRCTILE(X,P) operates
%   along the first non-singleton dimension.
%
%   Y = PRCTILE(X,P,"all") calculates percentiles of all the elements in X.
%
%   Y = PRCTILE(X,P,DIM) calculates percentiles along dimension DIM.
%
%   Y = PRCTILE(X,P,VECDIM) calculates percentiles of elements of X based
%   on the dimensions specified in vector VECDIM. For example,
%   PRCTILE(X,P,[1 2]) operates on the elements contained in the first and
%   second dimensions of X.
%
%   Y = PRCTILE(...,"Method",METHOD) specifies the method for calculating
%   percentiles. The value of METHOD must be:
%     "exact" - (default) Computes using sorting as explained below.
%     "approximate" - Computes using an approximation algorithm based on
%     t-digests.
%
%   Percentiles are specified using percentages, from 0 to 100.  For an N
%   element vector X, PRCTILE computes percentiles as follows:
%      1) The sorted values in X are taken as the 100*(0.5/N), 100*(1.5/N),
%         ..., 100*((N-0.5)/N) percentiles.
%      2) Linear interpolation is used to compute percentiles for percent
%         values between 100*(0.5/N) and 100*((N-0.5)/N).
%      3) The minimum or maximum values in X are assigned to percentiles
%         for percent values outside that range.
%
%   PRCTILE treats NaN values as missing values, and removes them.
%
%   Examples:
%      y = prctile(x,50); % the median of x
%      y = prctile(x,[2.5 25 50 75 97.5]); % a useful summary of x
%
%   See also IQR, MEDIAN, QUANTILE.

%   Copyright 1993-2024 The MathWorks, Inc.

if ~isreal(x)
    error(message('MATLAB:prctile:InvalidData'));
end
if ~(isreal(p) && (isnumeric(p) || islogical(p)) && isvector(p) && ~isempty(p)) || any(p < 0 | p > 100)
    error(message('MATLAB:prctile:BadPercents'));
end

% Make sure we are working in floating point to avoid rounding errors.
if isfloat(x)
    castOutput = true;
    yproto = dominantTypeHelper(x,p); % prototype output
    calcproto = yproto;               % for internal calculations
elseif isinteger(x)
    % integer types are up-cast to either double or single and the result
    % is down-cast back to the input type
    castOutput = true;
    yproto = dominantTypeHelper(x,p); % prototype output
    if contains(underlyingType(x), ["8" "16"])
        % single precision is enough
        x = single(x);
    else
        % Needs double precision
        x = double(x);
    end
    calcproto = x([]); % prototype needed for calculation
elseif isduration(x)
    castOutput = false;
    calcproto = x([]);
else
    % All other types (e.g. char, logical) are cast to double and the result is
    % double, ignoring the type of p
    castOutput = false;
    try
        x = double(x);
    catch
        % For types without double conversion, just proceed
    end
    calcproto = x([]);
end

% Parse dim
dimArgGiven = false;
if nargin > 2
    dim = varargin{1};
    if matlab.internal.math.checkInputName(dim,'all')
        dim = 1:ndims(x);
        dimArgGiven = true;
    elseif ~ischar(dim) && ~(isstring(dim) && isscalar(dim))
        if ~isvector(dim) || ~isreal(dim) || any(floor(dim) ~= dim) || any(dim < 1) || ~allfinite(dim) || isempty(dim) || ~isnumeric(dim)
            error(message('MATLAB:getdimarg:invalidDim'));
        end
        if ~isscalar(dim) && ~all(diff(sort(dim)))
            error(message('MATLAB:getdimarg:vecDimsMustBeUniquePositiveIntegers'));
        end
        dimArgGiven = true;
    end
end
sz = size(x);
if ~dimArgGiven
    dim = matlab.internal.math.firstNonSingletonDim(x);
end

% Parse name-value arguments
exactMethod = true;
delta = 1e3;
rs = [];
deltaRSGiven = false;
indStart = dimArgGiven + 1;
num = numel(varargin);
if rem(num - indStart,2) == 0
    error(message('MATLAB:prctile:KeyWithoutValue'));
end
for j = indStart:2:num
    % Allow undocumented parameters 'Delta' and 'RandStream' for backward compatibility.
    indName = matlab.internal.math.checkInputName(varargin{j},{'Method','Delta','RandStream'});
    if ~any(indName)
        error(message('MATLAB:prctile:ParseNames'));
    end
    if indName(1)
        indMethod = matlab.internal.math.checkInputName(varargin{j+1},{'exact', 'approximate'});
        if ~any(indMethod)
            error(message('MATLAB:prctile:InvalidMethod'))
        end
        exactMethod = indMethod(1);
    elseif indName(2)
        % Delta is undocumented, and TDigest does error checking.
        delta = varargin{j+1};
        deltaRSGiven = true;
    else
        % RandStream is undocumented, and TDigest does error checking.
        rs = varargin{j+1};
        deltaRSGiven = true;
    end
end
if deltaRSGiven && exactMethod
    % Delta and RandStream only apply to the approximate method.
    error(message('MATLAB:prctile:DeltaRSApproxOnly'));
end
if ~exactMethod
    checkSupportedForApproximateMethod(x);
end

% Permute the array so that the requested dimension is the first dim.
if dimArgGiven
    if all(dim > numel(sz))
        % Avoid unecessary calculations and return quickly.
        np = numel(p);
        r = ones(1,min(dim));
        r(end) = np;
        if castOutput
            y = cast(x, 'like', yproto);
        else
            y = x;
        end
        y = repmat(y,r);
        return;
    end
    if ~isequal(dim,1)
        dim = sort(dim);
        if issparse(x) && any(dim > 2)
            % permuting beyond second dimension not supported for sparse
            dim(dim > 2) = [];
        end
        perm = [dim setdiff(1:max(ndims(x),max(dim)),dim)];
        x = permute(x, perm);
        sz = size(x);
    end
end

% Special behavior for vectors needed after percentile computation. Store
% whether x is a vector now before x is reshaped during the computation.
XisVector = isvector(x);

% If X is empty, return all NaNs.
if isempty(x)
    if isequal(size(x),[0 0]) && ~dimArgGiven
        szout = size(p);
    else
        if dimArgGiven
            work_dim = 1:numel(dim);
        else
            work_dim = dim;
        end
        szout = sz; 
        szout(work_dim) = 1;
        szout(work_dim(1)) = numel(p);
    end
    y = createArray(szout,Like=calcproto,FillValue=missing);
else
    % Drop X's leading singleton dims, and combine its trailing dims.  This
    % leaves a matrix, and we can work along columns.
    if dimArgGiven
        work_dim = 1:numel(dim);
    else
        work_dim = dim;
    end
    work_dim = work_dim(work_dim <= numel(sz));
    nrows = prod(sz(work_dim));
    ncols = numel(x) ./ nrows;
    x = reshape(x, nrows, ncols);

    if ~isfloat(p)
        p = double(p);
    end
    
    if exactMethod
        if isequal(p,50) && ~issparse(x) && isnumeric(x) && ~isobject(x)
            % Utilize internal fast median
            y = matlab.internal.math.columnmedian(x,true);
        else
            x = sort(x,1);
            n = sum(~isnan(x), 1); % Number of non-NaN values in each column
            
            % For columns with no valid data, set n=1 to get nan in the result
            n(n==0) = 1;
            
            % If the number of non-nans in each column is the same, do all cols at once.
            if all(n == n(1))
                n = n(1);
                if isequal(p,50) % make the median fast
                    if rem(n,2) % n is odd
                        y = x((n+1)/2,:);
                    else        % n is even
                        y = (x(n/2,:) + x(n/2+1,:))/2;
                    end
                else
                    y = interpColsSame(x,p,n);
                end
                
            else
                % Get percentiles of the non-NaN values in each column.
                y = interpColsDiffer(x,p,n);
            end
        end
    else
        if issparse(x)
            numCols = size(x,2);
            y = zeros(numel(p),numCols,'like',calcproto);
            for k = 1:numCols
                v = x(:,k);
                td = matlab.internal.math.TDigestArray(full(v), delta, rs);
                y(:,k) = matlab.internal.math.tdigestICDF(td,p,min(v,[],1),max(v,[],1));
            end
        else
            td = matlab.internal.math.TDigestArray(x, delta, rs);
            y = matlab.internal.math.tdigestICDF(td,p,min(x,[],1),max(x,[],1));
        end
    end

    % Reshape Y to conform to X's original shape and size.
    szout = sz; 
    szout(work_dim) = 1;
    szout(work_dim(1)) = numel(p);
    y = reshape(y,szout);
end
% undo the DIM permutation
if dimArgGiven && ~isequal(dim,1)
     y = ipermute(y,perm);  
end

% If X is a vector, the shape of Y should follow that of P, unless an
% explicit DIM arg was given.
if ~dimArgGiven && XisVector
    y = reshape(y,size(p)); 
end

if castOutput
    y = cast(y, 'like', yproto);
end

function y = interpColsSame(x, p, n)
%INTERPCOLSSAME An aternative approach of 1-D linear interpolation which is
%   faster than using INTERP1Q and can deal with invalid data so long as
%   all columns have the same number of valid entries (scalar n).

% Make p a column vector. Note that n is assumed to be scalar.
if isrow(p)
    p = p';
end

% Form the vector of index values (numel(p) x 1)
r = (p/100)*n;
k = floor(r+0.5); % K gives the index for the row just before r
kp1 = k + 1;      % K+1 gives the index for the row just after r
r = r - k;        % R is the ratio between the K and K+1 rows

% Find indices that are out of the range 1 to n and cap them
k(k<1 | isnan(k)) = 1;
kp1 = min( kp1, n );

% Use simple linear interpolation for the valid percentages
y = (0.5+r).*x(kp1,:)+(0.5-r).*x(k,:);

% Make sure that values we hit exactly are copied rather than interpolated
exact = (r==-0.5);
if any(exact)
    y(exact,:) = x(k(exact),:);
end

% Make sure that identical values are copied rather than interpolated
same = (x(k,:)==x(kp1,:));
if any(same(:))
    x = x(k,:); % expand x
    y(same) = x(same);
end

function y = interpColsDiffer(x, p, n)
%INTERPCOLSDIFFER A simple 1-D linear interpolation of columns that can
%deal with columns with differing numbers of valid entries (vector n).

[nrows, ncols] = size(x);

% Make p a column vector. n is already a row vector with ncols columns.
if isrow(p)
    p = p';
end

% Form the grid of index values (numel(p) x numel(n))
r = (p/100)*n;
k = floor(r+0.5); % K gives the index for the row just before r
kp1 = k + 1;      % K+1 gives the index for the row just after r
r = r - k;        % R is the ratio between the K and K+1 rows

% Find indices that are out of the range 1 to n and cap them
k(k<1 | isnan(k)) = 1;
kp1 = min( kp1, n );

% Convert K and Kp1 into linear indices
offset = nrows*(0:ncols-1);
k =  k + offset;
kp1 = kp1 + offset;

% Use simple linear interpolation for the valid percentages.
% Note that NaNs in r produce NaN rows.
y = (0.5-r).*x(k) + (0.5+r).*x(kp1);

% Make sure that values we hit exactly are copied rather than interpolated
exact = (r==-0.5);
if any(exact(:))
    y(exact) = x(k(exact));
end

% Make sure that identical values are copied rather than interpolated
same = (x(k)==x(kp1));
if any(same(:))
    x = x(k); % expand x
    y(same) = x(same);
end

function proto = dominantTypeHelper(x,p)
% Adapted from internal.stats.dominantType
%DOMINANTTYPE return a prototype of the dominant input type
%   DOMINANTTYPE uses standard arithmetic rules to determine a prototype of
%   the dominant type from the inputs.
%
%   Example:
%   dominantType( single(1), double(2) ) => single
%   dominantType( single(1), gpuArray(2) ) => gpuArray(single)
%
% Get a prototype output to determine the type. Mostly X is in the nature
% of data here, and P is something that determine how we index into X, but
% we're going to treat them both as data for the purposes of determining
% the output type. If someone decided to put P on the GPU we want to
% respect that.

% Add together empty arrays so that the types propagate but no actual
% calculation is required.
x = x([]);
if isa(p,'double')
    proto = x;
else
    try
        proto = x + p([]);
        if issparse(x)
            proto = sparse(proto);
        end
    catch
        % if types will not add, first type encountered wins
        proto = x;
    end
end

function checkSupportedForApproximateMethod(x)
m = [];

if ~isnumeric(x)
    m = message('MATLAB:prctile:NumericRequiredForApprox');
elseif isobject(x)
    m = message('MATLAB:prctile:NoObjectsForApprox');
end
if ~isempty(m)
    throwAsCaller(MException(m.Identifier, '%s', getString(m)));
end
