function y = dispatchSP(op,x,k,dim,nanflag,nrm,spIn)
%MATLAB Code Generation Library Function

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

coder.internal.allowEnumInputs;
coder.internal.prefer_const(op,k,dim,nanflag,spIn,nrm,spIn);

ONE = coder.internal.indexInt(1);
omitnan = nanflag(1) == 'o';

[kleft,kright] = toKleftKright(k,spIn);
inclusive = (numel(k) == 2);

ytype = getYType(x,op);
y = coder.nullcopy(zeros(size(x),'like',ytype));

sizeX = coder.internal.indexInt(size(x));

[stride,numelPerStride] = getSizeInfo(sizeX,dim);

xInfo = struct('size',sizeX,'stride',stride, ...
    'numelPerStride',numelPerStride,'numel',coder.internal.indexInt(numel(x)));

windowStart = ONE;
windowEnd = ONE;

if isa(spIn,'int64')
    % Convert signed integer sample points to unsigned without altering the
    % difference between values. This prevents overflow when subtracting
    % sample points values
    sp = typecast(bitxor(spIn,intmin('int64')),'uint64');
elseif isnumeric(spIn) && ~isa(spIn,'uint64')
    sp = double(spIn);
else
    sp = spIn;
end

for ii = ONE:numel(sp)
    [windowStart,windowEnd] = ...
        getWindowBounds(windowStart,windowEnd,sp,ii,kleft,kright,inclusive);  
    inds = getAllInds(ii,xInfo,dim);
    y(inds) = dispatchToFun(op,x,omitnan,ytype,windowStart,windowEnd,xInfo,inds,ii,nrm,dim);
end
end

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

function ytype = getYType(x,op)
coder.internal.prefer_const(op);
coder.inline('always');
switch coder.const(op)
    case {'var','std'}
        ytype = coder.internal.scalarEg(x,x);
    case {'sum','prod'}
        ytype = coder.internal.scalarEg(sum(x));
    case 'mean'
        ytype = coder.internal.scalarEg(mean(x));
    case {'median','min','max','mad'}
        ytype = coder.internal.scalarEg(x);
end
end

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

function [windowStart,windowEnd] = ...
    getWindowBounds(windowStart,windowEnd,sp,k,kleft,kright,inclusive)
coder.internal.prefer_const(sp);

while(windowStart < k && sp(k) - sp(windowStart) > kleft)
    windowStart = windowStart+1;
end

windowEnd = max(windowEnd,k);
if inclusive
    while(windowEnd < numel(sp) && sp(windowEnd+1) - sp(k) <= kright)
        windowEnd = windowEnd+1;
    end
else
    while(windowEnd < numel(sp) && sp(windowEnd+1) - sp(k) < kright)
        windowEnd = windowEnd+1;
    end
end
end

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

function out = getAllInds(indexInDim,xInfo,dim)
coder.internal.prefer_const(xInfo,dim);

ONE = coder.internal.indexInt(1);
currentDimValue = ONE;
indexOut = ONE;

out = coder.nullcopy(zeros(1,xInfo.numelPerStride,coder.internal.indexIntClass));

indexInDimCount = ONE;
for ii = ONE:xInfo.numel
    if indexInDimCount > xInfo.stride
        if currentDimValue == coder.internal.indexInt(xInfo.size(dim))
            currentDimValue = ONE;
        else
            currentDimValue = currentDimValue + ONE;
        end
        indexInDimCount = ONE;
    end
    if currentDimValue == indexInDim
        out(indexOut) = ii;
        indexOut = indexOut + ONE;
        if indexOut > xInfo.numelPerStride
            break
        end
    end
    indexInDimCount = indexInDimCount + ONE;
end
end

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

function y = dispatchToFun(op,x,omitnan,ytype,windowStart,windowEnd,xInfo,inds,indInDim,nrm,dim)
coder.internal.prefer_const(op,omitnan,ytype,xInfo,dim);
ONE = coder.internal.indexInt(1);
nx = windowEnd - windowStart + ONE;

switch op
    case 'sum'
        identity = zeros('like',ytype);
        y = combineVectorElements(@plus,x,identity,dim,omitnan,nx,windowStart);
    case 'mean'
        y = meanInPlace(x,ytype,xInfo,omitnan,nx,windowStart,ONE,dim);
    case 'prod'
        identity = ones('like',ytype);
        y = combineVectorElements(@times,x,identity,dim,omitnan,nx,windowStart);
    case 'max'
        fun = @(windowStartInd,windowEndInd) ...
            coder.internal.vectorMinOrMaxInPlace('max',x,omitnan, ...
            windowStartInd,windowEndInd,xInfo.stride);
        y = expandVectorFun(fun,ytype,windowStart,windowEnd,inds,xInfo,indInDim);
    case 'min'
        fun = @(windowStartInd,windowEndInd) ...
            coder.internal.vectorMinOrMaxInPlace('min',x,omitnan, ...
            windowStartInd,windowEndInd,xInfo.stride);
        y = expandVectorFun(fun,ytype,windowStart,windowEnd,inds,xInfo,indInDim);
    case 'median'
        y = medianInPlace(x,ytype,xInfo,omitnan,nx,windowStart,inds,indInDim);
    case {'var','std','mad'}
        y = varStdMad(x,op,ytype,xInfo,omitnan,nx,windowStart,indInDim,nrm,inds,dim);
end
end

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

function [kleft,kright] = toKleftKright(k,sp)
% Split K into kleft and kright. The SP input is only needed to determine
% the output type. If either K or SP is double, the results are double.
coder.internal.allowEnumInputs;
coder.internal.prefer_const(k,sp);
SCALARK = coder.internal.isConstTrue(isscalar(k));
if isnumeric(sp) && isinteger(k) && SCALARK
    kright = k/2;
    kleft = k/2 - mod(k,cast(2,'like',k)); % Subtract 1 when k is odd
elseif SCALARK
    kleft = k/2;
    kright = kleft;
else
    kleft = k(1);
    kright = k(2);
end
end

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

function [stride,numelPerStride] = getSizeInfo(sizeX,dim)
ONE = coder.internal.indexInt(1);
if coder.internal.isConst(sizeX)
    % Evaluate the product of dimensions extrinsically to ensure const-ness
    strideTmp = coder.const(feval('prod',sizeX((1:numel(sizeX)<dim))));
    stride = coder.internal.indexInt(strideTmp);
    numelOutTmp = coder.const(feval('prod',sizeX((1:numel(sizeX)~=dim))));
    numelPerStride = coder.internal.indexInt(numelOutTmp);
else
    stride = ONE;
    numelPerStride = ONE;
    for ii = 1:numel(sizeX)
        if ii < dim
            stride = stride*sizeX(ii);
        end
        if ii ~= dim
            numelPerStride = numelPerStride*sizeX(ii);
        end
    end
end
end

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

function [y] = meanInPlace(x,ytype,xInfo,omitnan,nx,windowStart,nrm,dim)
coder.internal.prefer_const(ytype,xInfo,omitnan,nrm,dim);
ONE = coder.internal.indexInt(1);

identity = zeros('like',ytype);
[y,count] = combineVectorElements(@plus,x,identity,dim,omitnan,nx,windowStart);
for ii = ONE:numel(y)
    if count(ii) == 0
        y(ii) = coder.internal.nan();
    else
        if nrm == 0
            current_count = eml_max(ones('like',count),count(ii)-1);
        else
            current_count = count(ii);
        end
        y(ii) = y(ii)./cast(current_count,'like',y);
    end
end
end

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

function y = medianInPlace(x,ytype,xInfo,omitnan,nx,windowStart,inds,indInDim)
coder.internal.prefer_const(ytype,xInfo,omitnan)
ONE = coder.internal.indexInt(1);
coder.inline('always');
if omitnan
    nanflag = 'o';
else
    nanflag = 'i';
end
windowEnd = windowStart+nx-1;
fun = @(windowStartInd,windowEndInd) ...
    vmedian(x,nx,nanflag,windowStartInd-ONE,xInfo.stride,true);
y = expandVectorFun(fun,ytype,windowStart,windowEnd,inds,xInfo,indInDim);
end

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

function y = expandVectorFun(fun,ytype,windowStart,windowEnd,inds,xInfo,indInDim)
% Take a function that collapsese a vector and apply it to a dimension of
% a matrix (i.e. use it to collapse a matrix from n dims to n-1 dims)
%
% inds is a vector containing the linear index of the centers of each
% window. These all must have the same index in the dimension that is being
% collapsed - e.g. if we are collapsing the 2nd dimension of a 2D matrix,
% inds might equal sub2ind(size(x),1:size(x,1),repmat(3,[1 size(x,1)]))
%
% windowStart and windowEnd are the boundries of each window, defined by
% a number of strided before and after the center point
%
% indInDim is the current index in the dimension that is being collapsed
coder.internal.prefer_const(ytype,xInfo);
y = coder.nullcopy(zeros(1,xInfo.numelPerStride,'like',ytype));
ONE = coder.internal.indexInt(1);
startOffset = (windowStart-indInDim)*xInfo.stride;
endOffset = (windowEnd-indInDim)*xInfo.stride;
count = ONE;
for ii = 1:numel(inds)
    current_ind = inds(ii);
    y(count) = fun(startOffset+current_ind,endOffset+current_ind);
    count = count + ONE;
end
end

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

function y = varStdMad(x,op,ytype,xInfo,omitnan,nx,windowStart,indInDim,nrm,inds,dim)
coder.internal.prefer_const(op,ytype,xInfo,omitnan,nrm,dim);
ONE = coder.internal.indexInt(1);

if strcmp(op,'var') || strcmp(op,'std')
    NOTMAD = true;
    fun = @(x,nrm)meanInPlace(x,ytype,xInfo,omitnan,nx,windowStart,nrm,dim);
else
    NOTMAD = false;
    fun = @(x,~)medianInPlace(x,ytype,xInfo,omitnan,nx,windowStart,inds,indInDim);
end
meanOrMedian = fun(x,1);
x1 = cast(x,'like',ytype);
% Explicitly handing the dim stuff rather than using combineVectorElements
% because this allows for better control of numerical accuracy for the std
% and var cases -- this lets us do sum(abs(x-mean)^2) rather than 
% sum(abs(x)^2) - n*abs(mean)^2
for ii = ONE:numel(inds)
    offset = inds(ii);
    for jj = windowStart:windowStart+nx-1
        currentInd = (jj-indInDim)*xInfo.stride+offset;
        % Avoiding doing abs(x-y) in case x is unsigned
        if isfloat(x1)
            x1(currentInd) = abs(x1(currentInd)-meanOrMedian(ii));
        else
            % Special handing for ints in case they are unsigned
            % This should onlyg get hit by movmad because movvar and movstd
            % cannot have integer ytypes
            if x1(currentInd) >= meanOrMedian(ii)
                x1(currentInd) = x1(currentInd)-meanOrMedian(ii);
            else
                x1(currentInd) = meanOrMedian(ii) - x1(currentInd);
            end
        end
        if NOTMAD
            x1(currentInd) = x1(currentInd)^2;
        end
    end
end
y = fun(x1,nrm);
if strcmp(op,'std')
    y = sqrt(y);
end
end
