function window = chooseWindowSize(A, tau, sp)
%chooseWindowSizeUniform Determine a heuristic tuning size based on a tuning factor.
%
%   FOR INTERNAL USE ONLY -- This feature is intentionally undocumented.
%   Its behavior may change, or it may be removed in a future release.
%

% Copyright 2019-2020 The Mathworks, Inc.

%#codegen

coder.internal.prefer_const(tau);

coder.internal.assert(~isa(sp,'datetime'),'MATLAB:smoothdata:noWindowDatetimeSP');

[uniformSP,minDiff] = isuniform(sp);
if uniformSP
    % Tune on uniform nodes, scale the window size
    window = chooseWindowSizeNoSP(A,tau);
    window = window * (sp(2) - sp(1));
    window = clipWindowValue(window, class(sp));
    return;
end

if size(A,1) > 1
    if ~isempty(sp)
        % Limiting cases.
        if tau == 0
            window = (sp(end) - sp(1));
            window = clipWindowValue(window, class(sp));
            return;
        elseif tau == 1
            window = minDiff/2;
            window = clipWindowValue(window, class(sp));
            return;
        end
        
        % Interpolate onto a uniform grid
        Ai = A;
        if isinteger(sp)
            scaleFactor = getIntegerScaleFactor(sp);
            if (sp(end) > flintmax) || (sp(1) < -flintmax)
                spScale = rescaleIntegerData(sp, true);
            else
                spScale = double(sp);
            end
        else
            scaleFactor = (sp(end) - sp(1))/(length(sp)-1);
            spScale = sp;
        end
        if isa(spScale,'duration')
            tgrid = linspace(milliseconds(spScale(1)), ...
                milliseconds(spScale(end)), length(spScale))';
        else
            tgrid = linspace(spScale(1), spScale(end), length(spScale))';
        end
        
        % Have to convert integer data to double
        if isinteger(Ai) || islogical(Ai)
            if (max(Ai(:)) > flintmax)
                Ad = rescaleIntegerData(Ai, false);
            else
                Ad = double(Ai);
            end
        else
            Ad = A;
        end
        Ad = fillmissing(Ad, 'pchip', 'SamplePoints', spScale, ...
            'EndValues', 'extrap');
        if isa(spScale,'duration')
            Ad = interp1(milliseconds(spScale), Ad, tgrid, 'pchip');
        else
            Ad = interp1(spScale, Ad, tgrid, 'pchip');
        end
        
        % Use uniform tuning and scale
        window = chooseWindowSizeNoSP(Ad,tau);
        window = window * scaleFactor;
        window = clipWindowValue(window, class(spScale));
    else
        window = chooseWindowSizeNoSP(A,tau);
    end
else
    % Trivial case.
    if isa(sp, 'datetime') || isa(sp, 'duration')
        window = milliseconds(1);
    else
        window = ones('like',coder.internal.scalarEg(sp));
    end
end

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

function window = chooseWindowSizeNoSP(A,tau)
coder.noImplicitExpansionInFunction; % use of - when re-centering with same size operands: Ab - (same_size_as_Ab)
ONE = coder.internal.indexInt(1);

% Non-trivial case
if tau == 0
    window = size(A, 1);
    return
elseif tau == 1
    window = 1;
    return
end

% Get an estimate of the average cutoff frequency below which
% most of the input's energy is contained
if ~isfloat(A)
    Ab = double(A);
else
    Ab = fillmissing(A, 'pchip');
end

% Re-center the values
Ab = Ab - repmat(mean(Ab,1),size(Ab,1),1);
N = size(Ab, 1);
Ac = abs(fft(Ab, 2*N)).^2 / (2*N);
numRowsAc = coder.internal.indexInt(size(Ac,1));

if isreal(A)
    Ab(:,:) = Ac(1:N,:);
else
    for jj = ONE:size(Ab,2)
        for ii = ONE:size(Ab,1)
            if ii > 1
                Ab(ii,jj) = 0.5*Ac(ii,jj)+Ac(numRowsAc-ii+ONE,:);
            else
                Ab(ii,jj) = Ac(ii,jj);
            end
        end
    end
end

% Compute normalized cumulative sum, average it over columns
Ab = Ab ./ repmat(sum(Ab),size(Ab,1),1);
Ab = cumsum(Ab, 'omitnan');
Ab = mean(Ab(:, :), 2);

% Determine cutoff bandwidth
% rho = find(Ab > tau, 1, 'first');
rho = [];
for ii = ONE:numel(Ab)
    if Ab(ii)>tau
        rho = cast(ii,class(0));
        break
    end
end

if isempty(rho)
    % Columns are constant -- nothing to do
    winsz = 1;
else
    % Convert to moving average filter width
    winsz = ceil(sqrt((0.44294*2*N/max(rho(1)-1,1))^2+1));
    % rho can only be a scalar, but specifying the first element
    % permits better constant folding.
end
window = winsz;

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

function [tf,minDiff] = isuniform(sp)
coder.internal.prefer_const(sp);
coder.inline('always');
if numel(sp) <= 1
    tf = false;
    minDiff = [];
else
    if coder.internal.isConst(sp)
        SPACING = sp(2)-sp(1);
        DIFFSP = coder.const(feval('diff',sp));
        minDiff = coder.const(feval('min',DIFFSP));
        DIFFEQSPACING = coder.const(feval('eq',DIFFSP,SPACING));
        tf = coder.const(feval('all',DIFFEQSPACING));
    else
        spacing = sp(2)-sp(1);
        minDiff = spacing;
        tf = true;
        for ii = 3:numel(sp)
            currentDiff = sp(ii)-sp(ii-1);
            if currentDiff ~= spacing
                tf = false;
                if currentDiff < minDiff
                    minDiff = currentDiff;
                end
            end
        end
    end
end

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

function window = clipWindowValue(window, cls)
coder.inline('always')
coder.internal.prefer_const(window, cls)
if ~isinteger(window)
    if window == 0
        if strcmp(cls, 'duration') || strcmp(cls, 'datetime')
            for ii = 1:numel(window)
                window(ii) = milliseconds(realmin);
            end
        else
            window = realmin(cls);
        end
        return;
    end
    if isfinite(window)
        return;
    end
    if strcmp(cls, 'duration') || strcmp(cls, 'datetime')
        for ii = 1:numel(window)
            window(ii) = milliseconds(realmax);
        end
    else
        window = realmax(cls);
    end
end

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

function sf = getIntegerScaleFactor(sp)
% Return the scaling factor for the window size used for auto-tuned selection
N = uint64(numel(sp) - 1);
if (sp(1) < 0) && (sp(end) < 0)
    padBit = uint64(sp(1) == intmin('int64'));
    sf = cast((padBit + uint64(-sp(1)) - uint64(-sp(end))) / N, 'like', sp);
elseif (sp(1) > 0) && (sp(end) > 0)
    sf = cast((uint64(sp(end)) - uint64(sp(1))) / N, 'like', sp);
else
    padBit = uint64(sp(1) == intmin('int64'));
    sf = cast((uint64(sp(end)) + padBit + uint64(-sp(1))) / N, 'like', sp);
end

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

function td = rescaleIntegerData(t, mustBeSorted)
% Rescale int64/uint64 sample points

if isa(t, 'int64')
    % Convert the int64 to be positive and cast to uint64
    minVal = uint64(0);
    if min(t) < 0
        minVal = cast(-min(t), 'uint64');
    end
    idx = t > 0;
    tU64 = zeros(size(t), 'uint64');
    tU64(idx) = cast(t(idx), 'uint64') + minVal;
    tU64(~idx) = minVal - cast(-t(~idx), 'uint64');
    tAdjusted = tU64;
else
    tAdjusted = t;
end


tShifted = tAdjusted - min(tAdjusted);

if max(tShifted) > flintmax
    % Approximate sample points if we can't get things aligned
    % exactly
    td = (double(tShifted) / double(max(tShifted)));
    if mustBeSorted
        for i = 2:length(td)
            if (td(i) <= td(i-1))
                td(i) = td(i-1) + eps(td(i-1));
            end
        end
    end
else
    td = double(tShifted);
end
