function Y = pdist(X,dist,varargin)
%PDIST Pairwise distance between observations.
%   D = PDIST(X)
%   D = PDIST(X, DISTANCE)
%
%   The DISTANCE argument must be a string.
%
%   Example:
%      % Compute the ordinary Euclidean distance
%      X = gpuArray.randn(100, 5);        % some random points
%      D = pdist(X, 'euclidean');         % euclidean distance
%
%
%   See also pdist, gpuArray.

%   Copyright 2015-2023 The MathWorks, Inc.

if nargin < 2
    dist = "euc";
else
    if matlab.internal.datatypes.isScalarText(dist, false)
        allMethods = ["euclidean"; "seuclidean"; "cityblock"; "chebychev"; ...
            "fasteuclidean"; "fastseuclidean"; "fastsquaredeuclidean"; ...
            "mahalanobis"; "minkowski"; "cosine"; "correlation"; ...
            "spearman"; "hamming"; "jaccard"; "squaredeuclidean"];
        noGpuSupportMethods = ["fasteuclidean", "fastseuclidean", "fastsquaredeuclidean"];

        i = find(strncmpi(dist,allMethods,strlength(dist)));
        if isempty(i)
            error(message("stats:gpu:PdistDistanceNotSupported", dist));
        elseif length(i) > 1
            error(message("stats:pdist:AmbiguousDistance", dist));
        end

        dist = allMethods(i);
        
        if any(dist == noGpuSupportMethods)
            error(message("stats:gpu:PdistDistanceNotSupported", dist));
        end

        if dist == "squaredeuclidean"
            dist = "sqe"; % represent squared Euclidean
        else
            dist = extractBetween(dist,1,3);
        end

    else
        error(message("stats:gpu:PdistBadDistance"));
    end
end

% Integer/logical/char/anything data may be handled by a caller-defined
% distance function, otherwise it is converted to double.  Complex floating
% point data must also be handled by a caller-defined distance function.
if ~isfloat(X)
    warning(message("stats:pdist:DataConversion", underlyingType(X)));
    X = double(X);
elseif ~isreal(X)
    error(message("stats:pdist:InvalidData"));
end

if ~ismatrix(X)
    error(message("stats:pdist:UnsupportedND"));
end

[n,p] = size(X);

% Degenerate case, just return an empty of the proper size.
if n < 2
    Y = zeros(1,0,"like",X); % X was single/double, or cast to double
    return;
end

additionalArg = [];
if dist == "seu" % Standardized Euclidean weights by coordinate variance
    if nargin < 3
        additionalArg = var(X, "omitnan");
        if any(additionalArg == 0)
            warning(message("stats:pdist:ConstantColumns"));
        end
        additionalArg = 1./ additionalArg;
    else
        additionalArg = varargin{1};
        if ~(isvector(additionalArg) && length(additionalArg) == p...
                && all(additionalArg >= 0))
            error(message("stats:pdist:InvalidWeights"));
        end
        if any(additionalArg == 0)
            warning(message("stats:pdist:ZeroInverseWeights"));
        end
        %We will apply the inverse weight on each coordinate in the sum
        %of squares.
        additionalArg = 1./ (additionalArg.^2);
    end
elseif dist == "mah" % Mahalanobis
    if nargin < 3
        additionalArg = cov(X, "omitrows");
        [T,flag] = chol(additionalArg);
    else
        additionalArg = varargin{1};
        if ~isequal(size(additionalArg),[p,p])
            error(message("stats:pdist:InvalidCov"));
        end
        %use cholcov because we also need to check whether the matrix is symmetric
        [T,flag] = cholcov(additionalArg,0);
    end
    if flag ~= 0
        error(message("stats:pdist:SingularCov"));
    end
    additionalArg = inv(T);
elseif dist == "min" % Minkowski distance needs a third argument
    if nargin < 3  % use default value for exponent
        additionalArg = 2;
    elseif ( isscalar(varargin{1}) && varargin{1} > 0)
        additionalArg = gather(varargin{1}); % get exponent from input args (it must be a scalar)
    else
        error(message("stats:pdist:InvalidExponent"));
    end
elseif dist == "cos" % Cosine
    [X,flag] = normalizeX(X);
    if flag
        warning(message("stats:pdist:ZeroPoints"));
    end
elseif dist == "cor" % Correlation
    X = X - mean(X,2);
    [X, flag] = normalizeX(X);
    if flag
        warning(message("stats:pdist:ConstantPoints"));
    end
elseif dist == "spe" %Spearman
    X = tiedrank(X.').'; % treat rows as a series
    X = X - (p+1)/2; % subtract off the (constant) mean
    [X,flag] = normalizeX(X);
    if flag
        warning(message("stats:pdist:TiedPoints"));
    end
end

% Note that if there is any code for case 'che','euc' or 'cit' in the
% above switch dist block, the some code need to be repeated in the
% corresponding block below.
if strcmp(dist,"min") % Minkowski distance
    if isinf(additionalArg) %the exponent is inf
        dist = "che";
        additionalArg = [];
    elseif additionalArg == 2 %the exponent is 2
        dist = "euc";
        additionalArg = [];
    elseif additionalArg == 1 %the exponent is 1
        dist = "cit";
        additionalArg = [];
    end
end

% Call a mex file to compute distances for the standard distance measures
% and full real double or single data.
additionalArg = cast(additionalArg,underlyingType(X));
Y = internal.stats.gpu.pdist(X.',dist,additionalArg); % X is a real matrix

%---------------------------------------------
% Normalize the rows in X to have unit norm.
function [X,flag] = normalizeX(X)
% Rescale each row by largest element to prevent over/underflow, and
% compute the 2-norm.
Xmax = max(abs(X),[],2);
X2 = X ./ Xmax;
Xnorm = sqrt(sum(X2.^2, 2));

% The norm will be NaN for rows that are all zeros, fix that for the test
% below.
Xnorm(Xmax==0) = 0;

% The norm will be NaN for rows of X that have any +/-Inf. Those should be
% Inf, but leave them as is so those rows will not affect the test below.
% The points can't be normalized, so any distances from them will be NaN
% anyway.

% Find points that are effectively zero relative to the point with largest norm.
flag = any(Xnorm <= eps(max(Xnorm)));

% Back out the rescaling, and normalize rows of X to have unit 2-norm.
% Rows can't be normalized become all NaN.
Xnorm = Xnorm .* Xmax;
X = X ./ Xnorm;

