function y = detrend(x,varargin)
%MATLAB Code Generation Library Function

%   Limitations:
%       * If BP supplied and is not logical, it must be sorted in strictly
%       ascending order.
%       * Uses an inexpensive method of detecting non-unique or badly
%       conditioned problems for the purpose of issuing a warning. May not
%       match MATLAB in this respect.
%       * Variable-sizing is generally required if SamplePoints are
%       supplied.
%       * Variable-sizing is generally required if N is not a constant, unless
%       N is logical. There must be an upper bound on the value of N to avoid
%       dynamic memory allocation.
%       * If variable-sizing is disabled and BP is supplied, it must
%       contain only integer values in the range 1:n-2, where n = length(x)
%       if x is a vector or n = size(x,1), otherwise.

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

coder.internal.prefer_const(varargin);
ISROW = coder.internal.isConstTrue(isrow(x));
ISSCALAR = coder.internal.isConstTrue(isscalar(x));
ISMATRIX = coder.internal.isConstTrue(ismatrix(x));
ZERO = coder.internal.indexInt(0);
ONE = coder.internal.indexInt(1);
TWO = coder.internal.indexInt(2);
if ISROW && ~ISSCALAR
    y1 = detrend(x(1,:).',varargin{:});
    y = y1.';
    return
end
coder.internal.assert(isscalar(x) || ~isrow(x), ...
    'Coder:toolbox:detrend_rowVecSpecialCase');
if ~ISMATRIX
    y1 = detrend(x(:,:),varargin{:});
    y = reshape(y1,size(x));
    return
end
N = coder.internal.indexInt(size(x,1));
BPPossible = nargin >= 3 && (isnumeric(varargin{2}) || ...
    islogical(varargin{2}) || isduration(varargin{2}) || ...
    isdatetime(varargin{2}));

matchParam = @coder.internal.matchStringParameter;
if nargin == 1
    polyDeg = ONE;
    nextArg = ONE;
elseif coder.internal.isTextRow(varargin{1})
    coder.internal.assert(coder.internal.isConst(varargin{1}), ...
        'Coder:toolbox:OptionStringsMustBeConstant');
    if matchParam(varargin{1},'linear')
        polyDeg = ONE;
        nextArg = TWO;
    elseif matchParam(varargin{1},'constant')
        polyDeg = ZERO;
        nextArg = TWO;
    else
        polyDeg = ONE;
        nextArg = ONE;
        BPPossible = false;
    end
else
    pdt = varargin{1};
    coder.internal.assert(coder.internal.isConstTrue(isscalar(pdt)), ...
        'MATLAB:detrend:InvalidTrendInputType');
    if ~islogical(pdt)
        coder.internal.assert(isnumeric(pdt) && pdt >= 0 && ...
            floor(pdt) == pdt, ...
            'MATLAB:detrend:InvalidTrendInputType');
    end
    polyDeg = coder.internal.indexInt(pdt);
    nextArg = TWO;
end
if BPPossible && ~coder.internal.isTextRow(varargin{2})
    bpInputOriginalType = varargin{2};
    nextArg = nextArg + 1;
    BPSUPPLIED = true;
else
    bpInputOriginalType = zeros(0,'like',real(x));
    BPSUPPLIED = false;
end
OMITNAN = false;
if nextArg<=(nargin-1) && coder.internal.isTextRow(varargin{nextArg})
    if matchParam(varargin{nextArg},'omitnan') || matchParam(varargin{nextArg},'omitmissing')
        OMITNAN = coder.internal.get_eml_option('NonFinitesSupport'); %Bypass NaN handling if NaNs are not supported
        nextArg = nextArg + 1;
    elseif matchParam(varargin{nextArg},'includenan') || matchParam(varargin{nextArg},'includemissing')
        nextArg = nextArg + 1;
    end
end
parms = {'Continuous','SamplePoints'};
pop = struct( ...
    'CaseSensitivity',false, ...
    'PartialMatching','unique', ...
    'StructExpand',false, ...
    'IgnoreNulls',true);
pstruct = coder.internal.parseParameterInputs(parms,pop,varargin{nextArg:end});
CONTINUOUS = coder.internal.getParameterValue(pstruct.Continuous,true, ...
    varargin{nextArg:end});
sOriginalType = coder.internal.getParameterValue(pstruct.SamplePoints,double(0:N-1)', ...
    varargin{nextArg:end});
SPSUPPLIED = pstruct.SamplePoints > 0;
coder.internal.assert( ...
    isdatetime(bpInputOriginalType)  || isduration(bpInputOriginalType) || ...  breakpoints are durations or datetimes OR
    (coder.internal.isBuiltInNumeric(bpInputOriginalType) || islogical(bpInputOriginalType)) && ... breakpoints are numeric or logical AND are real and not sparse
    isreal(bpInputOriginalType) && ~issparse(bpInputOriginalType), ...
    'MATLAB:detrend:BreakpointsInvalid');
if SPSUPPLIED
    if isdatetime(sOriginalType)
        if isempty(sOriginalType)
            s = milliseconds([]);
            minS = s;
        else
            s = milliseconds(sOriginalType-sOriginalType(1));
            minS = s(1);
        end
    else
        if isduration(sOriginalType)
            s = milliseconds(sOriginalType);
        else
            s = sOriginalType;
        end
        minS = s(1);
        s = s - minS;
    end
     coder.internal.assert(coder.internal.isBuiltInNumeric(s) && ...
         (islogical(bpInputOriginalType) || isa(bpInputOriginalType,class(sOriginalType)) || isempty(bpInputOriginalType)), ...
         'MATLAB:detrend:BreakpointsInvalid');
    checkSamplePoints(s,x);
else
    minS = 1;
    s = sOriginalType;
end
if N == 0
    y = x;
    return
end
if N == 1
    y = x*0;
    return
end
VSIZE = coder.const(~coder.internal.eml_option_eq('VariableSizing','DisableInInference'));
if isduration(bpInputOriginalType)
    bpInput = milliseconds(bpInputOriginalType);
elseif isdatetime(bpInputOriginalType)
    bpInput = milliseconds(bpInputOriginalType - sOriginalType(1));
else
    bpInput = bpInputOriginalType;
end
if BPSUPPLIED && ~SPSUPPLIED && CONTINUOUS && ~VSIZE
    % Enforce legacy limitations to make fixed-size compilation possible.
    rtAssertValidLegacyBPInput(bpInput,N-2);
    bp = [0;bpInput(:)-1];
    lbp = coder.internal.indexInt(size(bp,1));
elseif SPSUPPLIED || BPSUPPLIED
    [bp,lbp] = processBP(bpInput(:),s,minS,N);
else
    bp = s(1);
    lbp = ONE;
end
if OMITNAN
    y = coder.nullcopy(x);
    if isvector(x)
        numDataPoints = coder.internal.indexInt(length(s));
        y = NanColumnDetrend(y,x,ONE,s,numDataPoints,polyDeg,bp,CONTINUOUS,VSIZE,lbp);
    else
        numColumns = coder.internal.indexInt(size(x,2));
        numDataPoints = coder.internal.indexInt(size(x,1));
        columnsWithoutNansInd = coder.nullcopy((ONE:numColumns)');
        numColumnsWithoutNans = ZERO;
        % Check columns of x for NaNs
        for ii = ONE:numColumns
            hasNan = false;
            for jj = ONE:numDataPoints
                if isnan(x(jj,ii))
                    hasNan = true;
                    break
                end
            end
            % Process columns with and without NaNs separately
            if hasNan
                y = NanColumnDetrend(y,x,ii,s,numDataPoints,polyDeg,bp,CONTINUOUS,VSIZE,lbp);
            else
                numColumnsWithoutNans = numColumnsWithoutNans + 1;
                columnsWithoutNansInd(numColumnsWithoutNans) = ii;
            end
        end
        if numColumnsWithoutNans > 0
            xNoNans = x(:,columnsWithoutNansInd(1:numColumnsWithoutNans));
            yTmp = chooseDetrendMethod(CONTINUOUS,polyDeg,VSIZE,xNoNans,bp,lbp,s);
            y(:,columnsWithoutNansInd(1:numColumnsWithoutNans)) = yTmp;
        end
    end
    return;
end
y = chooseDetrendMethod(CONTINUOUS,polyDeg,VSIZE,x,bp,lbp,s);

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

function y = chooseDetrendMethod(CONTINUOUS,polyDeg,VSIZE,x,bp,lbp,s)
coder.inline('always');
coder.internal.prefer_const(CONTINUOUS,polyDeg,VSIZE,lbp);
% Apply method
if CONTINUOUS
    if polyDeg == 0
        y = CCDetrend(x,bp,lbp,s);
    elseif VSIZE || coder.internal.isConst(polyDeg)
        y = CPPDetrend(x,bp,lbp,polyDeg,s(:));
    else
        % Enable fixed-size code generation to work as it did before when
        % polyDeg could only have been 0 or 1.
        coder.internal.assert(polyDeg == 1, ...
            'MATLAB:maxlhs');
        y = CPPDetrend(x,bp,lbp,ONE,s(:));
    end
else
    y = PPDetrend(x,bp,lbp,polyDeg,s(:));
end

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

function checkSamplePoints(s,x)
% Validate SamplePoints value
coder.internal.assert((isfloat(s) || isduration(s) || isdatetime(s)) && ...
    (isvector(s) || isempty(s)), ...
    'MATLAB:detrend:SamplePointsInvalidDatatype');
ns = coder.internal.indexInt(numel(s));
coder.internal.assert(ns == size(x,1)*~isempty(x), ...
    'MATLAB:detrend:SamplePointsLength');
if isfloat(s)
    coder.internal.assert(isreal(s),'MATLAB:detrend:SamplePointsComplex');
    coder.internal.assert(~issparse(s),'MATLAB:detrend:SamplePointsSparse');
end
coder.internal.assert(~coder.internal.anyNonFinite(s), ...
    'MATLAB:detrend:SamplePointsFinite');
DUPLICATE = int8(1);
NOTASCENDING = int8(2);
ec = int8(0);
for k = 1:ns-1
    if s(k) < s(k + 1)
    elseif ~(s(k) > s(k + 1))
        ec = DUPLICATE;
        break
    else
        ec = NOTASCENDING;
        break
    end
end
coder.internal.assert(ec ~= DUPLICATE,'MATLAB:detrend:SamplePointsDuplicate');
coder.internal.assert(ec ~= NOTASCENDING,'MATLAB:detrend:SamplePointsSorted');

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

function xs = subsum(x,begSeg,endSeg)
% xs = sum(x(begSeg:endSeg),1) with blocking for single precision.
nxsub = endSeg - begSeg + 1;
xsZERO = coder.internal.scalarEg(x,x);
xs = xsZERO;
if isa(x,'double')
    xs = sumab(x,begSeg,endSeg);
else
    % Blocked summation.
    BLOCKLENGTH = coder.internal.indexInt(1024);
    if nxsub > BLOCKLENGTH
        nblocks = coder.internal.indexDivide(nxsub,BLOCKLENGTH);
        for b = 1:nblocks
            begBlock = begSeg + (b - 1)*BLOCKLENGTH;
            endBlock = begBlock + BLOCKLENGTH - 1;
            xs = xs + sumab(x,begBlock,endBlock);
        end
        begBlock = begSeg + nblocks*BLOCKLENGTH;
    else
        begBlock = begSeg;
    end
    if endSeg >= begBlock
        xs = xs + sumab(x,begBlock,endSeg);
    end
end

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

function xs = sumab(x,begSeg,endSeg)
% xs = sum(x(begSeg:endSeg),1) without blocking.
xsZERO = coder.internal.scalarEg(x,x);
xs = xsZERO;
for k = begSeg:endSeg
    xs = xs + x(k);
end

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

function p = RealStrictAscending(bp1,lbp)
p = true;
for k = 2:lbp
    p = p && bp1(k) > bp1(k - 1);
end

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

function [bp,lbp] = processBP(bpInput,s,minS,N)
ONE = coder.internal.indexInt(1);
nbp = coder.internal.indexInt(length(bpInput));
if islogical(bpInput)
    bp = coder.nullcopy(zeros(nbp + 1,1,'like',s));
    % bp(1) = true;
    % bp1 = s(bp);
    lbp = ONE;
    bp(1) = 0;
    lastBp = min(nbp,N-1); %Prevents the last point being a breakpoint, trims excess breakpoints
    for k = 2:lastBp
        if bpInput(k)
            lbp = lbp + 1;
            bp(lbp) = s(k);
        end
    end
else
    bp = zeros(nbp + 1,1,'like',s);
    lbp = ONE;
    smax = s(end);
    for k = 1:nbp
        bpk = bpInput(k) - minS;
        if bpk > 0 && bpk < smax
            lbp = lbp + 1;
            bp(lbp) = bpk;
        end
    end
    coder.internal.assert(RealStrictAscending(bp,lbp), ...
        'Coder:MATLAB:detrend_BPnotsorted');
end
assert(lbp <= nbp + 1); %#<HINT>

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

function rtAssertValidLegacyBPInput(bp,maxbp)
for k = 1:numel(bp)-1
    bpk = bp(k);
    if floor(bpk) ~= bpk
        coder.internal.error('Coder:MATLAB:detrend_BPnonfinite');
    elseif bpk == bp(k+1)
        coder.internal.error('Coder:MATLAB:detrend_BPnonunique');
    elseif ~(bpk < bp(k+1))
        coder.internal.error('Coder:MATLAB:detrend_BPnotsorted');
    elseif bpk < 1 || bpk > maxbp
        coder.internal.error('Coder:MATLAB:detrend_BPoutofrange');
    end
end

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

function y = CCDetrend(x,bp,lbp,s)
% Continuous constant subtracts the mean of the first segment
ONE = coder.internal.indexInt(1);
XTYPE = coder.internal.scalarEg(x,x);
if lbp == 1
    endSeg = coder.internal.indexInt(numel(s));
else
    % [~,endSeg] = min(abs(s - bp1(2)));
    endSeg = coder.internal.bsearch(s,bp(2));
    d1 = bp(2) - s(endSeg);
    d2 = s(endSeg + 1) - bp(2);
    if d2 < d1
        endSeg = endSeg + 1;
    end
end
% segMean = mean(x(begSeg:endSeg,:),1);
partsum = @(x)subsum(x,ONE,endSeg);
segsum = coder.internal.applyVectorFunction(partsum,x,ONE,ONE,XTYPE);
segsum = segsum/cast(endSeg - ONE + 1,'like',segsum);
y = bsxfun(@minus,x,segsum);

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

function y = CPPDetrend(x,bp,lbp,polyDeg,s)
% Continuous, piecewise polynomial trend
% Assumes bp and s are column-vectors.
coder.internal.prefer_const(lbp,polyDeg);
if isa(x,'double') || ~isfloat(s)
    % Do everything in double precision.
    atype = 0;
else
    % Use the supplied types of s and bp, which might be double
    % while x is single.
    atype = zeros('like',s);
end
ns = coder.internal.indexInt(size(s,1));
scaleS = cast(s(end),'like',atype);
a = coder.nullcopy(zeros(ns,lbp,'like',atype));
for j = 1:lbp
    for i = 1:ns
        a(i,j) = cast(s(i),'like',atype) - cast(bp(j),'like',atype);
    end
end
a = a/scaleS;
a = max(a,0);
W = coder.nullcopy(zeros(ns,lbp*polyDeg + 1,'like',real(x)));
for panel = 0:polyDeg - 1
    b = cast(polyDeg - panel,'like',a);
    for j = 1:lbp
        for i = 1:ns
            W(i,j + panel*lbp) = a(i,j)^b;
        end
    end
end
for i = 1:ns
    W(i,end) = 1;
end

% Remove best fit
[p,rr] = coder.internal.qrsolve(W,x,'InferRank','NoWarn');

if size(W,1) < size(W,2) || rr < size(W,2)
    coder.internal.warning('MATLAB:detrend:PolyNotUnique');
end
y = x - W*p;

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

function y = PPDetrend(x,bp,lbp,polyDeg,s)
% Assumes s is column-vector.
coder.internal.prefer_const(polyDeg,lbp);
ZERO = coder.internal.indexInt(0);
y = zeros(size(x),'like',x);
ns = coder.internal.indexInt(size(s,1));
segments = coder.nullcopy(zeros(ns,1,coder.internal.indexIntClass));
% segments = sum(s(:) >= bp1',2);
for k = 1:ns
    segk = ZERO;
    sk = s(k);
    for j = 1:lbp
        if sk >= bp(j)
            segk = segk + 1;
        end
    end
    segments(k) = segk;
end
notUniquePoly = false;
for k = 1:lbp
    segidx = segments == k;
    nseg = nnz(segidx);
    if polyDeg == 0 || nseg == 1
        % Remove mean from each segment
        y(segidx,:) = bsxfun(@minus,x(segidx,:),mean(x(segidx,:),1));
    elseif nseg >= 2
        % Normalize before fitting polynomial
        if isa(x,'double') || ~isfloat(s)
            % Do everything in double precision.
            atype = 0;
        else
            % Use the supplied types of s and bp, which might be double
            % while x is single.
            atype = zeros('like',s);
        end
        a = cast(s(segidx),'like',atype);
        % Normalize a = (a - mean(a))/std(a);
        a = a - mean(a(:));
        % Since a is already centered, std(a) = norm(a)/sqrt(length(a)-1);
        stda = norm(a)/sqrt(nseg - 1);
        a = a/stda;
        % Construct the Vandermonde matrix
        V = coder.internal.vander(a,polyDeg);
        % Solve least squares problem to get polynomial coefficients
        [p,rr] = coder.internal.qrsolve(V,x(segidx,:),'InferRank','NoWarn');

        if size(V,1) < size(V,2) || rr < size(V,2)
            notUniquePoly = true;
        end
        % Remove best fit
        y(segidx,:) = bsxfun(@minus,x(segidx,:),V*p);
    end
end
if notUniquePoly
    coder.internal.warning('MATLAB:detrend:PolyNotUnique');
end

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

function y = NanColumnDetrend(y,x,column,s,numDataPoints,polyDeg,bp,CONTINUOUS,VSIZE,lbp)
%Remove NaNs, detrend, then replace NaNs
coder.internal.prefer_const(s,numDataPoints,polyDeg,bp,CONTINUOUS,VSIZE,lbp);
[xNoNans, sNoNans] = RemoveNans(x,s,column);
if(~isempty(xNoNans))
    yNoNans = chooseDetrendMethod(CONTINUOUS,polyDeg,VSIZE,xNoNans,bp,lbp,sNoNans);
else
    yNoNans = cast([],class(x));
end
y = ReplaceNans(y,x,column,yNoNans,numDataPoints);

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

function [xNoNans, sNoNans] = RemoveNans(x,s,column)
ZERO = coder.internal.indexInt(0);
ONE = coder.internal.indexInt(1);
numDataPoints = coder.internal.indexInt(length(s));
indNotNan = coder.nullcopy((ONE:numDataPoints)');
numNotNan = ZERO;
currentInd = ONE;
%Simultaneously locate and count NaNs in data
for ii = 1:numDataPoints
    if ~isnan(x(ii,column))
        indNotNan(currentInd) = ii;
        currentInd = currentInd+1;
        numNotNan = numNotNan+ONE;
    end
end

sNoNans = coder.nullcopy(s(1:numNotNan));
xNoNans = coder.nullcopy(x(1:numNotNan,column));

for ii = ONE:numNotNan
    sNoNans(ii) = s(indNotNan(ii));
    xNoNans(ii) = x(indNotNan(ii),column);
end

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

function y = ReplaceNans(y,x,column,yNoNans,numDataPoints)
%Repopulate with yNoNans, filling in with NaNs anywhere that x is also a
%NaN
coder.internal.prefer_const(numDataPoints);
NAN = nan(1,'like',x);
currentInd = coder.internal.indexInt(1);
for ii = coder.internal.indexInt(1):numDataPoints
    if ~isnan(x(ii,column))
        y(ii,column) = yNoNans(currentInd);
        currentInd = currentInd+1;
    else
        y(ii,column) = NAN;
    end
end

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