function xy = cov(x,varargin)
%MATLAB Code Generation Library Function

%   Notes:
%   1. Output on the diagonal is always real (as it should be).
%      Example:  cov([1,2i;nan,4;7,2])
%   2. If the input is variably sized and is [] at run-time, returns []
%      rather than NaN.

%   Copyright 1984-2023 The MathWorks, Inc.
%#codegen

narginchk(1,4);
eml_invariant(isfloat(x) || islogical(x) || ischar(x), ...
    'Coder:toolbox:unsupportedClass','cov',class(x));
eml_invariant(ismatrix(x),'MATLAB:cov:InputDim');
coder.internal.prefer_const(varargin);
if islogical(x) || ischar(x) % Cast logical and char inputs to double.
    xy = cov(double(x),varargin{:});
    return
elseif nargin == 1
    domle = false;
    opt = INCLUDENAN;
elseif nargin == 2 && coder.internal.isCharOrScalarString(varargin{1})
    domle = false;
    opt = parseTrailingStringOption(varargin{1});
    coder.internal.assert(opt ~= INVALID,'MATLAB:cov:unknownFlag');
elseif nargin == 2 && coder.internal.isConstTrue(isscalar(varargin{1}))
    domle = isequal(varargin{1},1);
    coder.internal.assert(domle || isequal(varargin{1},0), ...
        'Coder:toolbox:invalidCovNormOption');
    opt = INCLUDENAN;
elseif nargin == 3 && coder.internal.isCharOrScalarString(varargin{2}) && ...
        coder.internal.isConstTrue(isscalar(varargin{1}))
    domle = isequal(varargin{1},1);
    coder.internal.assert(domle || isequal(varargin{1},0), ...
        'Coder:toolbox:invalidCovNormOption');
    opt = parseTrailingStringOption(varargin{2});
    coder.internal.assert(opt ~= INVALID,'MATLAB:cov:unknownFlag');
else
    y = varargin{1};
    coder.internal.assert(numel(x) == numel(y), ...
        'MATLAB:cov:XYlengthMismatch');
    coder.internal.assert(isfloat(y) || islogical(y) || ischar(y), ...
        'Coder:toolbox:unsupportedClass','cov',class(y));
    xy0 = coder.internal.scalarEg(x,y);
    xy = zeros(2,'like',xy0);
    % We intercept the scalar case because in MATLAB cov(1,3)-->zeros(2),
    % but cov([1,3])-->2. Apart from that edge case, we can calculate
    % cov(x,y) as cov([x(:),y(:)]).
    if isscalar(x) && isscalar(y)
        NAN = coder.internal.nan('like',xy);
        if ~isfinite(x(1))
            xy(1,1) = NAN;
            xy(2,1) = NAN;
            xy(1,2) = NAN;
        end
        if ~isfinite(y(1))
            xy(2,1) = NAN;
            xy(1,2) = NAN;
            xy(2,2) = NAN;
        end
    else
        xy(:,:) = cov([x(:),cast(y(:),'like',xy0)],varargin{2:end});
    end
    return
end
if coder.internal.isConstTrue(isvector(x))
    xy = local_cov(x(:),domle,opt);
elseif isrow(x) % Variable-size matrix becomes a row at run-time.
    xy = local_cov(x(1,:).',domle,opt);
else
    xy = local_cov(x,domle,opt);
end

%--------------------------------------------------------------------------

function opt = parseTrailingStringOption(s)
coder.inline('always');
coder.internal.prefer_const(s);
matchParam = @coder.internal.matchStringParameter;
if matchParam(s,'partialrows')
    opt = PAIRWISE;
elseif matchParam(s,'omitrows')
    opt = COMPLETE;
elseif matchParam(s,'includenan') || matchParam(s,'includemissing')
    opt = INCLUDENAN;
else
    opt = INVALID;
end

%--------------------------------------------------------------------------

function c = local_cov(x,domle,missingOption)
coder.internal.prefer_const(domle,missingOption);
coder.internal.assert(ismatrix(x),'Coder:toolbox:InputsMustBe2D');
m = coder.internal.indexInt(size(x,1));
n = coder.internal.indexInt(size(x,2));
if eml_is_const(size(x)) && isequal(size(x),[0,0])
    c = coder.internal.nan(class(x));
    return
end
c = zeros(n,'like',x);
if isempty(x)
    c(:) = coder.internal.nan(class(c));
    return
end
XZERO = zeros('like',x);
IZERO = coder.internal.indexInt(0);
IONE = coder.internal.indexInt(1);
if missingOption == INCLUDENAN
    xnan = false(size(x)); % unreferenced
    anynans = false;
else
    xnan = isnan(x);
    anynans = vany(xnan);
end
if missingOption ~= PAIRWISE || ~anynans
    if anynans
        % Remove rows with nans.
        nrows = IZERO;
        for k = 1:m
            keeprow = true;
            for j = 1:n
                if xnan(k,j)
                    keeprow = false;
                    break
                end
            end
            if keeprow
                nrows = nrows + 1;
                if nrows < k
                    for j = 1:n
                        x(nrows,j) = x(k,j);
                    end
                end
            end
        end
    else
        nrows = m;
    end
    XONE = XZERO + 1;
    LDA = coder.internal.indexInt(size(x,1));
    % if nrows == 1, c = zeros(n) is already the correct answer.
    if nrows == 0
        c(:) = coder.internal.nan(class(x));
    elseif nrows >= 2
        nx = cast(nrows,'like',real(x));
        % Center the data.
        for j = 1:n
            % Calculate the mean.
            muj = XZERO;
            for i = 1:nrows
                muj = muj + x(i,j);
            end
            muj = muj/nx;
            % Center the column.
            for i = 1:nrows
                x(i,j) = x(i,j) - muj;
            end
        end
        if domle
            alpha = XONE/nx;
        else
            alpha = XONE/(nx - 1);
        end
        % c = xc' * xc / denom;
        % c = alpha*x(1:nrows,:)'*x(1:nrows,:);
        if coder.isRowMajor
            for i = nrows+1:m
                for j = 1:n
                    x(i,j) = 0;
                end
            end
            x = x*sqrt(alpha);
            c = x'*x;
        else
            c(:) = 0;
            c = coder.internal.blas.xgemm('C','N',n,n,nrows, ...
                alpha,x,IONE,LDA,x,IONE,LDA,XZERO,c,IONE,n);
        end
    end
else % pairwise with some NaNs
    for j = 1:n
        for i = 1:j
            mui = XZERO;
            muj = XZERO;
            nrows = IZERO;
            for k = 1:m
                if ~xnan(k,i) && ~xnan(k,j)
                    nrows = nrows + 1;
                    mui = mui + x(k,i);
                    muj = muj + x(k,j);
                end
            end
            if nrows == 0
                c(i,j) = coder.internal.nan(class(x));
            elseif nrows >= 2
                nx = cast(nrows,'like',real(x));
                mui = mui/nx;
                muj = muj/nx;
                if domle
                    denom = nx;
                else
                    denom = nx - 1;
                end
                cij = XZERO;
                for k = 1:m
                    if ~xnan(k,i) && ~xnan(k,j)
                        cij = cij + (x(k,i) - mui)*(x(k,j) - muj);
                    end
                end
                cij = cij/denom;
                c(i,j) = cij;
            end
            if i < j
                c(j,i) = c(i,j);
            end
        end
    end
end

%--------------------------------------------------------------------------

function p = vany(x)
% Efficiently evaluate p = any(x(:)).
for k = 1:numel(x)
    if x(k)
        p = true;
        return
    end
end
p = false;

%--------------------------------------------------------------------------

function n = INCLUDENAN
coder.inline('always');
n = int8(0);

function n = COMPLETE
coder.inline('always');
n = int8(1);

function n = PAIRWISE
coder.inline('always');
n = int8(2);

function n = INVALID
coder.inline('always');
n = int8(-1);

%--------------------------------------------------------------------------
