function varargout = gradient(x,varargin)
%MATLAB Code Generation Library Function

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

coder.internal.assert(nargin > 0,'MATLAB:minrhs');
coder.internal.assert(isa(x,'float'), ...
    'Coder:toolbox:gradient_2');
NOUT = eml_max(1,nargout);
ONE = coder.internal.indexInt(1);
TWO = coder.internal.indexInt(2);
NESTED = coder.isRowMajor || ~coder.internal.get_eml_option('EnableMatrixFlattening');
coder.internal.assert(nargin - 1 <= ndims(x), ...
    'MATLAB:gradient:InvalidInputs');
% Implicit spacing is only supported under certain circumstances.
coder.internal.assert(nargin == 1 || ...
    (nargin == 2 && isscalar(varargin{1})) || ...
    (nargin == 2 && isvector(x)) || ...
    nargin == ndims(x) + 1, ...
    'MATLAB:gradient:InvalidInputs');
coder.internal.assert(nargout <= ndims(x), ...
    'Coder:toolbox:gradient_tooManyOutputArguments');
% We don't support returning two outputs for vector inputs.
coder.internal.assert(nargout <= 1 || ~isvector(x), ...
    'MATLAB:maxlhs');
if nargout <= 1 && nargin <= 2 && iscolumn(x)
    % For column vector input, varargout{1} is computed by differencing
    % along the column. Note that nargout <= 1 is redundant to the previous
    % error check, but including it here may help eliminate this case with
    % variable-size x that isn't a variable-length vector.
    varargout{1} = gradFlat(x,ONE,varargin{:});
else
    for k = coder.unroll(1:NOUT)
        % Except for k == 1 and k == 2, varargout{k} is computed by
        % differencing along dimension k. Output varargout{1} is computed by
        % differencing along the rows, and varargout{2} is computed by
        % differencing down the columns. Column vectors are a special case
        % that should already have been taken care of.
        if k == 1
            dim = TWO;
        elseif k == 2
            dim = ONE;
        else
            dim = k;
        end
        n = coder.internal.indexInt(size(x,dim));
        if n < 2
            if nargin > k
                ytype = outputType(x,varargin{1});
            else
                ytype = outputType(x);
            end
            varargout{k} = zeros(size(x),'like',ytype);
        elseif nargin > 1
            if nargin > k
                h = varargin{k};
            else
                h = varargin{1};
            end
            if NESTED
                varargout{k} = gradLoop(x,dim,h);
            else
                varargout{k} = gradFlat(x,dim,h);
            end
        else
            if NESTED
                varargout{k} = gradLoop(x,dim);
            else
                varargout{k} = gradFlat(x,dim);
            end
        end
    end
end

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

function ytype = outputType(x,h)
% Returns the output type and the type to which denominators should be
% cast.
xtype = coder.internal.scalarEg(x,x);
if nargin < 2
    ytype = xtype;
else
    % Spacing input does not change the class of the output type, but it
    % does determine whether or not it is complex.
    if ~isreal(xtype) || isreal(h)
        ytype = xtype;
    else
        ytype = complex(xtype);
    end
end

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

function y = gradFlat(x,dim,h)
% Perform gradient computations using flattened indexing. The computation
% is blocked, but in some cases it could be less cache friendly than the
% ND nested loop approach, which performs two passes per block to maximize
% data locality. The same technique could easily be incorporated here, but
% the result was slightly slower when tested.
coder.internal.prefer_const(dim);
DEFAULTDEN = nargin < 3;
n = coder.internal.indexInt(size(x,dim));
blocksize = coder.internal.prodsize(x,'below',dim);
pagesize = blocksize*n;
npages = coder.internal.prodsize(x,'above',dim);
if DEFAULTDEN
    ytype = outputType(x);
else
    den = h2den(n,h);
    ytype = outputType(x,h);
end
DENTYPE = real(coder.internal.scalarEg(x,x));
y = coder.nullcopy(zeros(size(x),'like',ytype));
if n >= 2
    for p = 1:npages
        pageoffset = (p - 1)*pagesize;
        lo = pageoffset;
        hi = pageoffset + blocksize;
        for k = 1:blocksize
            if DEFAULTDEN
                y(lo + k) = (x(hi + k) - x(lo + k));
            else
                y(lo + k) = (x(hi + k) - x(lo + k))/cast(den(1),'like',DENTYPE);
            end
        end
        coder.gpu.kernel;
        for j = 2:n-1
            lo = pageoffset + (j - 2)*blocksize;
            mid = pageoffset + (j - 1)*blocksize;
            hi = pageoffset + j*blocksize;
            coder.gpu.kernel;
            for k = 1:blocksize
                if DEFAULTDEN
                    y(mid + k) = (x(hi + k) - x(lo + k))/2;
                else
                    y(mid + k) = (x(hi + k) - x(lo + k))/ ...
                        cast(den(j),'like',DENTYPE);
                end
            end
        end
        lo = pageoffset + (n - 2)*blocksize;
        hi = pageoffset + (n - 1)*blocksize;
        for k = 1:blocksize
            if DEFAULTDEN
                y(hi + k) = (x(hi + k) - x(lo + k));
            else
                y(hi + k) = (x(hi + k) - x(lo + k))/cast(den(n),'like',DENTYPE);
            end
        end
    end
else
    y(:) = 0;
end

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

function den = h2den(n,h)
% Convert a spacing input to a vector of denominators in the difference
% calculations.
coder.internal.prefer_const(n);
nh = coder.internal.indexInt(numel(h));
htype = coder.internal.scalarEg(h,h); % converts logical to double.
den = coder.nullcopy(zeros(n,1,'like',htype));
if nh == 1
    % Scalar input implies constant spacing.
    den(1) = 1;
    for k = 2:n-1
        den(k) = 2;
    end
    den(n) = 1;
    den = den*h(1);
else
    den(1) = h(2) - h(1);
    for k = 2:nh-1
        den(k) = h(k + 1) - h(k - 1);
    end
    den(n) = h(n) - h(n - 1);
end

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

function y = gradLoop(x,dim,h)
% Generate nested loop structure to perform gradient computations with
% reasonable data locality.
coder.internal.prefer_const(dim);
if nargin < 3
    h = ones('like',real(x));
end
ONE = coder.internal.indexInt(1);
ND = coder.internal.indexInt(coder.internal.ndims(x));
ytype = outputType(x,h);
n = coder.internal.indexInt(size(x,dim));
den = h2den(n,h);
y = coder.nullcopy(zeros(size(x),'like',ytype));
if coder.isRowMajor
    y = initLoop(x,dim,y,ONE);
    y = diffLoop(x,dim,y,den,ONE);
else
    y = initLoop(x,dim,y,ND);
    y = diffLoop(x,dim,y,den,ND);
end

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

function y = initLoop(x,dim,y,loopdim,varargin)
% Copy the first term in the difference computation to the output array.
% The effect is like shifting the vectors x to the left one unit along
% dimension dim. The first element of each vector of x is not needed here,
% but we must make special provision to copy the last element of each
% vector over twice, since it is used in both the penultimate difference
% (x(n) - x(n - 2))/2 and in the last difference (x(n) - x(n - 1)).
coder.inline('always');
coder.internal.prefer_const(dim,loopdim);
ONE = coder.internal.indexInt(1);
ND = coder.internal.indexInt(coder.internal.ndims(x));
ROWMAJOR = coder.isRowMajor;
if coder.const(loopdim == dim)
    hi = coder.internal.indexInt(size(x,dim));
    if ROWMAJOR
        for k = 1:hi-1
            y = initLoop(x,dim,y,loopdim + 1,varargin{:},k);
        end
        y = initLoopLast(x,dim,y,loopdim + 1,varargin{:},hi);
    else
        for k = 1:hi-1
            y = initLoop(x,dim,y,loopdim - 1,k,varargin{:});
        end
        y = initLoopLast(x,dim,y,loopdim - 1,hi,varargin{:});
    end
elseif coder.const(loopdim >= ONE && loopdim <= ND)
    hi = coder.internal.indexInt(size(x,loopdim));
    if ROWMAJOR
        for k = 1:hi
            y = initLoop(x,dim,y,loopdim + 1,varargin{:},k);
        end
    else
        for k = 1:hi
            y = initLoop(x,dim,y,loopdim - 1,k,varargin{:});
        end
    end
else
    % y(...,k,...) = x(...,k + 1,...), k < size(x,dim)
    subsr = varargin;
    subsr{dim} = varargin{dim} + 1;
    y(varargin{:}) = x(subsr{:});
end

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

function y = initLoopLast(x,dim,y,loopdim,varargin)
% Copy the last element along dimension dim to y.
coder.inline('always');
coder.internal.prefer_const(dim,loopdim);
ONE = coder.internal.indexInt(1);
ND = coder.internal.indexInt(coder.internal.ndims(x));
ROWMAJOR = coder.isRowMajor;
if coder.const(loopdim >= ONE && loopdim <= ND)
    hi = coder.internal.indexInt(size(x,loopdim));
    if ROWMAJOR
        for k = 1:hi
            y = initLoopLast(x,dim,y,loopdim + 1,varargin{:},k);
        end
    else
        for k = 1:hi
            y = initLoopLast(x,dim,y,loopdim - 1,k,varargin{:});
        end
    end
else
    % y(...,n,...) = x(...,n,...), n = size(x,dim)
    y(varargin{:}) = x(varargin{:});
end

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

function y = diffLoop(x,dim,y,den,loopdim,varargin)
% Peform the difference calculations. What we subtract from each vector of
% y along dimension dim is the corresponding vector of x along dimension
% dim shifted to the right one unit. Obviously special provision must be
% made for the first element of each vector of x, since it is used to
% compute both first and second elements of the vectors of y: (x(2) - x(1))
% and (x(3) - x(1))/2, respectively.
coder.inline('always');
coder.internal.prefer_const(dim,loopdim);
ONE = coder.internal.indexInt(1);
TWO = coder.internal.indexInt(2);
ND = coder.internal.indexInt(coder.internal.ndims(x));
ROWMAJOR = coder.isRowMajor;
DENTYPE = real(coder.internal.scalarEg(x,x));
if coder.const(loopdim == dim)
    hi = coder.internal.indexInt(size(x,dim));
    if ROWMAJOR
        y = diffLoopFirst(x,dim,y,den,loopdim + 1,varargin{:},ONE);
        for k = TWO:hi
            y = diffLoop(x,dim,y,den,loopdim + 1,varargin{:},k);
        end
    else
        y = diffLoopFirst(x,dim,y,den,loopdim - 1,ONE,varargin{:});
        for k = TWO:hi
            y = diffLoop(x,dim,y,den,loopdim - 1,k,varargin{:});
        end
    end
elseif coder.const(loopdim >= ONE && loopdim <= ND)
    hi = coder.internal.indexInt(size(x,loopdim));
    if ROWMAJOR
        for k = 1:hi
            y = diffLoop(x,dim,y,den,loopdim + 1,varargin{:},k);
        end
    else
        for k = 1:hi
            y = diffLoop(x,dim,y,den,loopdim - 1,k,varargin{:});
        end
    end
else
    % y(...,k,...) = (y(...,k,...) - x(...,k - 1,...))/den(k), 2<=k<=size(x,dim)
    subsr = varargin;
    subsr{dim} = varargin{dim} - 1;
    y(varargin{:}) = (y(varargin{:}) - x(subsr{:}))/ ...
        cast(den(varargin{dim}),'like',DENTYPE);
end

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

function y = diffLoopFirst(x,dim,y,den,loopdim,varargin)
% Compute just the first element of y along dimension dim.
coder.inline('always');
coder.internal.prefer_const(dim,loopdim);
ONE = coder.internal.indexInt(1);
ND = coder.internal.indexInt(coder.internal.ndims(x));
ROWMAJOR = coder.isRowMajor;
DENTYPE = real(coder.internal.scalarEg(x,x));
if coder.const(loopdim >= ONE && loopdim <= ND)
    hi = coder.internal.indexInt(size(x,loopdim));
    if ROWMAJOR
        for k = 1:hi
            y = diffLoopFirst(x,dim,y,den,loopdim + 1,varargin{:},k);
        end
    else
        for k = 1:hi
            y = diffLoopFirst(x,dim,y,den,loopdim - 1,k,varargin{:});
        end
    end
else
    % y(...,1,...) = (y(...,1,...) - x(...,1,...))/den(1)
    y(varargin{:}) = (y(varargin{:}) - x(varargin{:}))/ ...
        cast(den(1),'like',DENTYPE);
end

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