function y = prctile(varargin)
%
%   FOR INTERNAL USE ONLY -- This feature is intentionally undocumented.
%   Its behavior may change, or it may be removed in a future release.
%

%   Y = PRCTILE(X,P) returns percentiles of the values in X.  P is a scalar
%   or a vector of percent values.
%
%   Y = PRCTILE(X,P,'all') calculates percentiles of all the elements in X.
%   The smallest dimension index of Y has length LENGTH(P)
%
%   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 the vector VECDIM.
%
%   Y = PRCTILE(...,'Method','exact'/'approximate') does not support
%   'approximate'.
%
%   Examples:
%      y = matlab.internal.math.prctile(x,50); % the median of x
%      y = matlab.internal.math.prctile(x,[2.5 25 50 75 97.5]); % a useful summary of x

%   Copyright 2021 The MathWorks, Inc.

par = inputParser();
par.addRequired('x');
par.addRequired('p');
par.addOptional('dim',1,@(x) isnumeric(x) || validateDimAll(x));
par.addParameter('Method', 'exact');
par.addParameter('Delta',1e3);
par.addParameter('RandStream',[]);

par.parse(varargin{:});

x = par.Results.x;
p = par.Results.p;
dim = par.Results.dim;
method = par.Results.Method;

% Cannot have complex X
if ~isreal(x)
    error('MATLAB:prctile:NoComplexNamed','X');
end
if ~isvector(p) || numel(p) == 0 || any(p < 0 | p > 100) || ~isreal(p)
    error(message('MATLAB:prctile:BadPercents'));
end
validatestring(method,"exact");

% 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. Sometimes P will be data-like and X is table-like, and
% if someone decided to put P on the GPU we want to respect that.
yproto = dominantTypeHelper(x,p); % prototype output
calcproto = yproto;               % for internal calculations

% Make sure we are working in floating point to avoid rounding errors.
if isfloat(x)
    castOutput = true;
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;
    outType = underlyingType(x);
    if ismember(outType, ["int8" "uint8" "int16" "uint16"])
        % single precision is enough
        x = single(x);
    else
        % Needs double precision
        x = double(x);
    end
    calcproto = x([]); % prototype needed for calculation
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 like duration without double conversion, just proceed
    end
end

% Figure out which dimension prctile will work along.
sz = size(x);
if ismember('dim', par.UsingDefaults) 
    dim = find(sz ~= 1,1);
    if isempty(dim)
        dim = 1; 
    end
    dimArgGiven = false;
else
    if validateDimAll(dim)
        dim = 1:ndims(x);
    else
        % dim is numeric
        % Check the validity of dim
        if ~isreal(dim) || any(floor(dim) ~= ceil(dim)) || any(dim < 1) || any(~isfinite(dim))
            error(message('MATLAB:getdimarg:invalidDim'));
        end
        if ~isscalar(dim) && ~all(diff(sort(dim)))
            error(message('MATLAB:getdimarg:vecDimsMustBeUniquePositiveIntegers'));
        end
    end  
    % Permute the array so that the requested dimension is the first dim.
    if ~isequal(dim,1)
        nDimsX = ndims(x);
        dim = sort(dim);
        perm = [dim setdiff(1:max(nDimsX,max(dim)),dim)];
        x = permute(x, perm);
    end
    sz = size(x);
    dimArgGiven = true;
end

% If X is empty, return all NaNs.
if isempty(x)
    if isequal(x,[]) && ~dimArgGiven
        y = nan(size(p),'like',calcproto);
    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);
        y = nan(szout,'like',calcproto);
    end

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);
    
    % Algorithm for method = 'exact'
    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

    % 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 && isvector(x)
    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 bool = validateDimAll(dim)
bool = ((ischar(dim) && isrow(dim)) || ...
     (isstring(dim) && isscalar(dim) && (strlength(dim) > 0))) && ...
     strncmpi(dim,'all',max(strlength(dim), 1));

function proto = dominantTypeHelper(varargin)
%Copy of 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)

%   Copyright 2014-2019 The MathWorks, Inc.

% Add together empty arrays so that the types propagate but no actual
% calculation is required.
proto = varargin{1}([]);
for ii=2:nargin
    try
        proto = proto + varargin{ii}([]);
    catch
        % if types will not add, first type encountered wins    
    end
end
