function [coef, pval] = corr(x, varargin)
%CORR Linear or rank correlation for GPU data.
%   RHO = CORR(X)
%   RHO = CORR(X,Y,...)
%   [RHO,PVAL] = CORR(...)
%   [...] = CORR(...,'PARAM1',VAL1,'PARAM2',VAL2,...)
%
%   See also CORR, GPUARRAY.

%   See toolbox/stats/stats/corr for implementation details.

%   Copyright 2019-2023 The MathWorks, Inc.

if nargin > 1
    [varargin{:}] = convertStringsToChars(varargin{:});
end

if nargin < 1 || isempty(x)
    error(message('stats:corr:TooFewInputs'));
end
[n,p1] = size(x);

% Only x given, compute pairwise rank correlations
if (nargin < 2) || ischar(varargin{1})
    corrXX = true;
    y = x;
    p2 = p1;

% Both x and y given, compute the pairwise rank cross correlations
else
    y = varargin{1};
    varargin = varargin(2:end);
    if size(y,1) ~= n
        error(message('stats:corr:InputSizeMismatch'));
    end
    corrXX = false;
    p2 = size(y,2);
end

pnames = {'type'  'rows' 'tail' 'weights'};
dflts  = {'p'     'a'    'both' []};
[type,rows,tail,wts] = internal.stats.parseArgs(pnames,dflts,varargin{:});

% Validate the weights parameter
isweighted = ~isempty(wts);
if isweighted
    validateattributes(wts,{'double' 'single'}, {'size',[n 1],'nonnegative'},'corr','Weights');
end

% Validate the rows parameter.
rowsChoices = {'all' 'complete' 'pairwise'};
if ischar(rows)
    i = find(strncmpi(rows,rowsChoices,length(rows)));
    if isscalar(i)
        rows = rowsChoices{i}(1);
    end
end
switch rows
case 'a' % 'all'
    % Missing values are left in, so all cols have the same length.
case 'c' % 'complete'
    % Complete rows with missing values are removed, so all cols have the
    % same length.
    ok = ~any(isnan(x),2);
    if ~corrXX
        ok = ok & ~any(isnan(y),2);
    end
    if isweighted
        wts(wts==0) = nan("like", wts);
        ok = ok & ~isnan(wts);
    end
    n = sum(ok);
    if n < size(x,1)
        x = x(ok,:);
        if ~corrXX
            y = y(ok,:);
        end
        if isweighted
            wts = wts(ok);
        end
    end
case 'p' % 'pairwise'
    % Missing values are removed pairwise, so each column pair may have a
    % different length, we'll have to see. However, any NaN weights can be
    % removed up front.
    if isweighted
        ok = ~isnan(wts);
        n = sum(ok);
        if n<length(wts)
            wts = wts(ok);
            x = x(ok,:);
            if ~corrXX
                y = y(ok,:);
            end
        end
    end

otherwise
    error(message('stats:corr:UnknownRows'));
end

% Can't do anything without at least two observations.
if n < 2
    outType = internal.stats.dominantType(x,y);
    coef = NaN(p1,p2,"like",outType);
    pval = NaN(p1,p2,"like",outType);
    return
end

% Validate the type parameter, and set a handle to the correlation function.
typeChoices = {'pearson' 'kendall' 'spearman'};
if ischar(type)
    i = find(strncmpi(type,typeChoices,length(type)));
    if isscalar(i)
        type = typeChoices{i}(1);
    end
end
switch type
case 'p' % Pearson linear correlation
    corrFun = @internal.stats.corrPearson;
case 'k' % Kendall rank correlation (GPU-specific version)
    corrFun = @internal.stats.gpu.corrKendall;
case 's' % Spearman rank correlation
    corrFun = @internal.stats.corrSpearman;
otherwise
    error(message('stats:corr:UnknownType'));
end

% Validate the tail parameter.
tailChoices = {'left','both','right'};
if ischar(tail) && (size(tail,1)==1)
    i = find(strncmpi(tail,tailChoices,length(tail)));
    if isempty(i)
        i = find(strncmpi(tail,{'lt','ne','gt'},length(tail)));
    end
    if isscalar(i)
        tail = tailChoices{i}(1);
    elseif isempty(i)
        error(message('stats:corr:UnknownTail'));
    end
else
    error(message('stats:corr:UnknownTail'));
end

% Some things we just can't do for complex data.
complexdata = ~(isreal(x) || all(imag(x(:))==0));
if ~complexdata && ~corrXX
    complexdata = ~(isreal(y) || all(imag(y(:))==0));
end
if complexdata
    if nargout > 1
        error(message('stats:corr:PValComplexData'));
    elseif type ~= 'p'
        error(message('stats:corr:RankCorComplexData'));
    end
end

if corrXX
    if nargout<2
        coef = corrFun(rows,tail,wts,x);
    else
        [coef,pval] = corrFun(rows,tail,wts,x);
        pval = (pval + pval')/2;
    end

    coef = (coef + coef')/2;
else
    if nargout<2
        coef = corrFun(rows,tail,wts,x,y);
    else
        [coef,pval] = corrFun(rows,tail,wts,x,y);
    end
end
end