function varargout = ischange(A,varargin)
%MATLAB Code Generation Library Function

%   Limitations:
%   * Results with single-precision inputs may not match MATLAB. Reason:
%     with the 'mean' and 'variance' methods, the generated code does not
%     use double precision for intermediate computations. With the 'linear'
%     method, double precision is used only for SamplePoints, but more of
%     the intermediate computations and temporary arrays are still
%     performed in single precision.

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

narginchk(1,inf);
nargoutchk(0,3);
coder.internal.prefer_const(varargin);
ONE = coder.internal.indexInt(1);
nextArg = ONE;
nVarArgs = coder.internal.indexInt(nargin - 1);
% Parse second input, if supplied. It could be the method string, a dim
% argument, or the first parameter name.
if nVarArgs == 1 && coder.internal.isTextRow(varargin{1})
    % Supplying this code just makes for a slightly better error message
    % when there are only two inputs, since the second input must be a
    % method input, not a parameter name.
    method = parseMethod(varargin{1});
    coder.internal.assert(method ~= UNRECOGNIZED, ...
        'MATLAB:ischange:MethodInvalid')
    LOOKFORDIM = false;
    nextArg = nextArg + 1;
elseif nVarArgs >= 1 && coder.internal.isTextRow(varargin{1})
    method1 = parseMethod(varargin{1});
    if method1 ~= UNRECOGNIZED
        method = method1;
        LOOKFORDIM = true;
        nextArg = nextArg + 1;
    else
        method = MEAN;
        LOOKFORDIM = false;
    end
else
    LOOKFORDIM = true;
    method = MEAN;
end
% Default minimum separation is 1 for 'mean' and 2 for the others
if method == MEAN
    separation = coder.internal.indexInt(1);
else
    separation = coder.internal.indexInt(2);
end
ISTABLE = isa(A,'table');
ISTIMETABLE = isa(A,'timetable');
ISTABULAR = ISTABLE || ISTIMETABLE;
if LOOKFORDIM && nextArg <= nVarArgs && isnumeric(varargin{nextArg})
    coder.internal.assert(~ISTABULAR,'MATLAB:ischange:DimensionTable');
    coder.internal.assertValidDim(varargin{nextArg});
    dim = coder.internal.indexInt(varargin{nextArg});
    nextArg = nextArg + 1;
elseif ISTABULAR
    dim = ONE;
else
    dim = coder.internal.preferConstNonSingletonDim(A);
end
parms = {'SamplePoints','DataVariables','MaxNumChanges','Threshold','OutputFormat'};
poptions = struct( ...
    'CaseSensitivity',false, ...
    'PartialMatching','unique', ...
    'StructExpand',false, ...
    'IgnoreNulls',false);
pstruct = coder.internal.parseParameterInputs( ...
    parms,poptions,varargin{nextArg:end});
coder.internal.assert(ISTABULAR || pstruct.DataVariables == 0, ...
    'MATLAB:ischange:DataVariablesArray');
% Process Threshold input, if present.
coder.internal.assert(pstruct.Threshold == 0 || ...
    pstruct.MaxNumChanges == 0, ...
    'MATLAB:ischange:MaxNumChangesThreshold');
threshold = coder.internal.getParameterValue(pstruct.Threshold, ...
    1,varargin{nextArg:end});
coder.internal.assert((isnumeric(threshold) || islogical(threshold)) && ...
    isscalar(threshold) && ...
    isreal(threshold) && ...
    threshold >= 0 && ...
    isfinite(threshold), ...
    'MATLAB:ischange:Threshold');
% Process MaxNumChanges input, if present.
if pstruct.MaxNumChanges == 0
    maxNumChanges = [];
else
    nchg = coder.internal.getParameterValue(pstruct.MaxNumChanges, ...
        [],varargin{nextArg:end});
    coder.internal.assert((isnumeric(nchg) || islogical(nchg)) && ...
        isscalar(nchg) && isreal(nchg) && ...
        nchg == floor(nchg) && ...
        nchg > 0, ...
        'MATLAB:ischange:MaxNumChanges');
    maxNumChanges = coder.internal.indexInt(nchg(1));
end
% Process SamplePoints input, if present.
if pstruct.SamplePoints == 0
    if ISTIMETABLE
        spts = A.Properties.RowTimes;
    else
        spts = [];
    end
else
    spts = coder.internal.getParameterValue(pstruct.SamplePoints, ...
        [],varargin{nextArg:end});
    coder.internal.assert(~ISTIMETABLE, ...
        'MATLAB:ischange:SamplePointsTimeTable');
end
% Process OutputFormat input, if present.
if pstruct.OutputFormat == 0
    tabularOutput = false;
else
    coder.internal.assert(ISTABULAR,'MATLAB:ischange:OutputFormatArray');
    outputFormat = coder.internal.getParameterValue( ...
        pstruct.OutputFormat,'logical',varargin{nextArg:end});
    fmt = validatestring(outputFormat,{'logical','tabular'},'isLocalExtrema','OutputFormat');
    coder.internal.assert(coder.internal.isConst(fmt), ...
        'Coder:toolbox:OptionStringsMustBeConstant');
    tabularOutput = strcmp(fmt,'tabular');
end
if method == LINEAR
    coder.internal.errorIf(isdatetime(spts), ...
        'MATLAB:ischange:SamplePointsDatetimeLinear');
    % Always use a double abscissa x.
    if isempty(spts)
        nv = coder.internal.indexInt(size(A,dim));
        x = zeros(nv,1);
        for k = 1:nv
            x(k) = k;
        end
    elseif isduration(spts)
        x = milliseconds(spts);
        checkSamplePoints(x,A,dim);
    else
        x = double(spts);
        checkSamplePoints(x,A,dim);
    end
else
    if ~isempty(spts)
        % We ignore the sample points, but we check them anyway.
        checkSamplePoints(spts,A,dim);
    end
    x = []; % Don't create any x
end
if ISTABULAR
    dataVars = coder.internal.getParameterValue( ...
        pstruct.DataVariables,1:width(A),varargin{nextArg:end});
    [varargout{1:nargout}] = ischangeTable( ...
        A,method,x,maxNumChanges,separation,threshold,dataVars,tabularOutput);
else
    coder.internal.assert(isfloat(A),'MATLAB:ischange:FirstInputInvalid');
    coder.internal.assert(isreal(A),'MATLAB:ischange:ComplexInput');
    penalty = cast(threshold,'like',A);
    [varargout{1:nargout}] = ischangeArray( ...
        A,method,x,maxNumChanges,separation,penalty,dim);
end

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

function [TF,S1,S2] = ischangeTable( ...
    A,method,x,numchanges,separation,thresh,dataVars,tabularOutput)
coder.internal.prefer_const(method,numchanges,separation,thresh, ...
    dataVars,tabularOutput);
ONE = coder.internal.indexInt(1);
TFTmp = false(size(A));
if nargout > 1
    S1 = A(:,dataVars);
    if nargout > 2
        S2 = S1;
    end
end
coder.internal.assert(coder.internal.isConst(dataVars), ...
    'Coder:toolbox:DataVarsMustBeConstant');
if coder.internal.isTextRow(dataVars)
    nv = ONE;
else
    nv = coder.const(coder.internal.indexInt(length(dataVars)));
end
varIdx = coder.nullcopy(coder.internal.indexInt(zeros(1,nv)));
coder.unroll;
for j = 1:nv
    if coder.internal.isTextRow(dataVars)
        vj = dataVars;
        Avj = A.(vj);
        tfcolidx = varnameToColumnIndex(vj,A.Properties.VariableNames);
    elseif iscell(dataVars)
        vj = dataVars{j};
        Avj = A.(vj);
        tfcolidx = varnameToColumnIndex(vj,A.Properties.VariableNames);
    else
        vj = dataVars(j);
        Avj = A.(vj);
        tfcolidx = vj;
    end
    varIdx(j) = tfcolidx;
    penalty = cast(thresh,'like',Avj);
    coder.internal.assert(isfloat(Avj), ...
        'MATLAB:ischange:UnsupportedTableVariable');
    coder.internal.assert(isreal(Avj), ...
        'MATLAB:ischange:ComplexInput');
    coder.internal.assert(iscolumn(Avj) && ~isempty(Avj), ...
        'MATLAB:ischange:NonVectorTableVariable');
    if nargout < 2
        TFTmp(:,tfcolidx) = ischangeArray( ...
            Avj,method,x,numchanges,separation,penalty,ONE);
    elseif nargout < 3
        [TFTmp(:,tfcolidx),S1.(j)] = ischangeArray( ...
            Avj,method,x,numchanges,separation,penalty,ONE);
    else
        [TFTmp(:,tfcolidx),S1.(j),S2.(j)] = ischangeArray( ...
            Avj,method,x,numchanges,separation,penalty,ONE);
    end
end
if tabularOutput
    if istimetable(A)
        TF = array2timetable(TFTmp(:,varIdx), ...
            'VariableNames', {A.Properties.VariableNames{varIdx}}, ...
            'RowTimes', A.Properties.RowTimes, ...
            'DimensionNames', A.Properties.DimensionNames);
    else
        TF = array2table(TFTmp(:,varIdx), ...
            'VariableNames',{A.Properties.VariableNames{varIdx}}, ...
            'DimensionNames',A.Properties.DimensionNames, ...
            'RowNames',A.Properties.RowNames);
    end
else
    TF = TFTmp;
end

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

function j = varnameToColumnIndex(name,varnames)
% Convert a table variable name to its corresponding column index.
coder.internal.prefer_const(name,varnames);
j = coder.internal.indexInt(0);
coder.unroll;
for k = 1:coder.internal.indexInt(length(varnames))
    if strcmp(name,varnames{k})
        j = k;
        return
    end
end

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

function [TF,S1,S2] = ischangeArray( ...
    A,method,x,numchanges,separation,penalty,dim)
% Changepoint detection for arrays
coder.internal.prefer_const(method,numchanges,separation,penalty,dim)
vlen = coder.internal.indexInt(size(A,dim));
FITSEGMENTS = nargout > 1;
% Process a couple of special cases, sparse and constant dim with only
% one output.
if issparse(A)
    if FITSEGMENTS
        [TF,S1,S2] = isChangeArraySparse(A, ...
            method,x,numchanges,separation,penalty,dim);
    else
        TF = isChangeArraySparse(A, ...
            method,x,numchanges,separation,penalty,dim);
    end
    return
elseif ~FITSEGMENTS && coder.internal.isConst(dim)
    % coder.internal.applyVectorFunction is row-major friendly, so use it
    % if we only need one output and dim is constant.
    isChangeFun = @(v)ischangeArrayColumn(v, ...
        method,x,numchanges,separation,penalty);
    TF = coder.internal.applyVectorFunction(isChangeFun,A,dim,vlen,false);
    return
end
% Process full inputs when dim is not constant or at least the second
% output is required.
TF = coder.nullcopy(false(size(A)));
if FITSEGMENTS
    S1 = coder.nullcopy(A);
    S2 = coder.nullcopy(A);
end
if coder.internal.isConstTrue(isrow(A)) && dim == 2
    if FITSEGMENTS
        [TF(1,:),S1(1,:),S2(1,:)] = ischangeArrayColumn(A(1,:).', ...
            method,x,numchanges,separation,penalty);
    else
        TF(1,:) = ischangeArrayColumn(A(1,:).', ...
            method,x,numchanges,separation,penalty);
    end
elseif dim == 1
    ncols = coder.internal.prodsize(A,'above',1);
    for k = 1:ncols
        if FITSEGMENTS
            [TF(:,k),S1(:,k),S2(:,k)] = ischangeArrayColumn(A(:,k), ...
                method,x,numchanges,separation,penalty);
        else
            TF(:,k) = ischangeArrayColumn(A(:,k), ...
                method,x,numchanges,separation,penalty);
        end
    end
else
    npages = coder.internal.prodsize(A,'above',dim);
    vlen = coder.internal.indexInt(size(A,dim));
    stride = coder.internal.prodsize(A,'below',dim);
    pagesize = stride*npages;
    for page = 1:npages
        pageOffset = (page - 1)*pagesize;
        for j = 1:stride
            vOffset = pageOffset + (j - 1);
            v = coder.nullcopy(zeros(vlen,1,'like',A));
            for k = 1:vlen
                v(k) = A(vOffset + (k - 1)*stride + 1);
            end
            if FITSEGMENTS
                [tfv,s1v,s2v] = ischangeArrayColumn(v,method,x, ...
                    numchanges,separation,penalty);
            else
                tfv = ischangeArrayColumn(v,method,x,numchanges, ...
                    separation,penalty);
            end
            for k = 1:vlen
                TF(vOffset + (k - 1)*stride + 1) = tfv(k);
            end
            if FITSEGMENTS
                for k = 1:vlen
                    S1(vOffset + (k - 1)*stride + 1) = s1v(k);
                end
                for k = 1:vlen
                    S2(vOffset + (k - 1)*stride + 1) = s2v(k);
                end
            end
        end
    end
end

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

function [TF,S1,S2] = isChangeArraySparse( ...
    A,method,x,numchanges,separation,penalty,dim)
coder.internal.prefer_const(method,numchanges,separation,penalty,dim);
ONE = coder.internal.indexInt(1);
if dim == 1
    ncols = coder.internal.indexInt(size(A,2));
    TF = logical(spalloc(size(A,1),size(A,2),nnz(A)));
    if nargout > 1
        S1 = coder.nullcopy(A);
        S2 = coder.nullcopy(A);
        for k = 1:ncols
            [TF(:,k),S1(:,k),S2(:,k)] = ischangeArrayColumn( ...
                A(:,k),method,x,numchanges,separation,penalty);
        end
    else
        for k = 1:ncols
            TF(:,k) = ischangeArrayColumn( ...
                A(:,k),method,x,numchanges,separation,penalty);
        end
    end
elseif dim == 2
    if nargout > 1
        [TFt,S1t,S2t] = isChangeArraySparse( ...
            A.',method,x,numchanges,separation,penalty,ONE);
        TF = TFt.';
        S1 = S1t.';
        S2 = S2t.';
    else
        TFt = isChangeArraySparse( ...
            A.',method,x,numchanges,separation,penalty,ONE);
        TF = TFt.';
    end
else
    TF = logical(spalloc(size(A,1),size(A,2),0));
    if nargout > 1
        S1 = A;
        S2 = spalloc(size(A,1),size(A,2),0);
    end
end

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

function [tf,s1,s2] = ischangeArrayColumn(y,method,x,Kmax,Lmin,penalty)
% Changepoint detection for one column
if issparse(y)
    y = full(y);
end
ZERO = coder.internal.indexInt(0);
nmax = preferConstN(x,y);
n = ZERO;
hasInfs = false;
nNaNs = ZERO;
locNaN = zeros(nmax,1,coder.internal.indexIntClass);
XEMPTY = coder.internal.isConstTrue(isempty(x));
% Omit NaNs by taking them out, computing with no NaNs, and then accounting
% for the NaN locations in the original data.
for k = 1:nmax
    yk = y(k);
    if isinf(yk)
        hasInfs = true;
    elseif isnan(yk)
        nNaNs = nNaNs + 1;
        locNaN(nNaNs) = k;
    else
        n = n + 1;
        y(n) = y(k);
        if ~XEMPTY
            x(n) = x(k);
        end
    end
end
coder.internal.assert(~hasInfs,'MATLAB:ischange:InfInput');
% Allocate space for cp vector.
cp = ones(nmax + 1,1,coder.internal.indexIntClass);
if nmax < 2 % in general: numel(a) < Lmin
    % No change points if data is empty or scalar
    ncp = ZERO;
elseif isempty(Kmax)
    % ischange(a) and ischange(a,...,'Threshold',penalty)
    [cp,ncp] = cpmanual(cp,y,n,method,Lmin,penalty,x);
else
    % ischange(a,...,'MaxNumChanges',Kmax)
    [cp,ncp] = cpmulti(cp,y,n,method,Lmin,Kmax,x);
end
tf = false(size(y));
for k = 1:ncp
    tf(cp(k)) = true;
end
if nargout > 1
    [s1,s2] = fitsegments(y,n,cp,ncp,method,x);
    if nNaNs > 0
        % Account for the omitted NaNs.
        iNaN = locNaN(nNaNs);
        for k = nmax:-1:iNaN+1
            s1(k) = s1(k - nNaNs);
            s2(k) = s2(k - nNaNs);
        end
        s1(iNaN) = nan;
        s2(iNaN) = nan;
        for j = nNaNs-1:-1:1
            khi = locNaN(j + 1) - 1;
            iNaN = locNaN(j);
            for k = khi:-1:iNaN + 1
                s1(k) = s1(k - j);
                s2(k) = s2(k - j);
            end
            s1(iNaN) = nan;
            s2(iNaN) = nan;
        end
    end
end

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

function [cp,ncp] = cpmanual(cp,a,n,method,Lmin,penalty,x)
% Multiple changepoints using 'Threshold' = penalty
if Lmin < 1 || Lmin > n || n <= 1
    ncp = coder.internal.indexInt(0);
else
    cost0 = cpnochange(a,n,method,x);
    [cost1,emptyFlag] = cpsingle(a,n,method,Lmin,x);
    if emptyFlag || cost1 + penalty > cost0
        ncp = coder.internal.indexInt(0);
    else
        [cp,ncp] = ischangePelt(cp,a,n,method,penalty,x,Lmin,true);
    end
end

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

function cost = cpnochange(a,n,statistic,x)
% Evaluate cost function of segment.
ONE = coder.internal.indexInt(1);
[abar,avar] = segmentMeanAndVariance(a,ONE,n,true);
na = cast(n,'like',a);
if statistic == MEAN
    % cost = Syy
    cost = na*avar;
elseif statistic == VARIANCE
    cost = na*log(avar);
else % statistic = 'linear'
    % cost = Syy - Sxy^2 / Sxx
    [xbar,xvar] = segmentMeanAndVariance(x,ONE,n,true);
    cost = coder.internal.scalarEg(a,x);
    for k = 1:n
        cost = cost + (a(k) - abar)*(x(k) - xbar);
    end
    cost = na*avar - (cost/na)*(cost/xvar);
end

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

function [cost,emptyFlag] = cpsingle(a,n,statistic,Lmin,x)
% Cost function if data had only one changepoint
% Evaluate costs for all possible combinations of two segments
coder.internal.prefer_const(statistic);
if n < 2*Lmin
    emptyFlag = true;
    cost = zeros('like',a);
    return
end
emptyFlag = false;
if statistic == MEAN
    [costSegment1,costSegment2] = costTwoSegmentsMean(a,n);
elseif statistic == VARIANCE
    [costSegment1,costSegment2] = costTwoSegmentsVariance(a,n);
else
    [costSegment1,costSegment2] = costTwoSegmentsLinear(a,x,n);
end
% Find the optimal cost of partitioning the data into two segments by
% finding the minimum of cost(1:i-1) + cost(i:n) for i = 2:n
cost = costSegment1(Lmin) + costSegment2(Lmin + 1);
for k = Lmin+1:n-Lmin
    tk = costSegment1(k) + costSegment2(k + 1);
    if tk < cost
        cost = tk;
    end
end

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

function [costSegment1,costSegment2] = costTwoSegmentsMean(y,n)
% Evaluate cost function for all possible combinations of two segments
ny = coder.internal.indexInt(numel(y));
if coder.internal.isConst(ny)
    costSegment1 = coder.nullcopy(zeros(ny,1,'like',y));
    costSegment2 = coder.nullcopy(zeros(ny,1,'like',y));
    for k = n+1:ny
        costSegment1(k) = 0;
        costSegment2(k) = 0;
    end
else
    assert(n <= ny); %#HINT
    costSegment1 = coder.nullcopy(zeros(n,1,'like',y));
    costSegment2 = coder.nullcopy(zeros(n,1,'like',y));
end
% Cost of all possible first segments:
% same iterative formula as the builtin instead of cpnochange on a(1:ix)
ymean = zeros('like',y);
Syy = ymean;
for k = 1:n
    ydelta = y(k) - ymean;
    ymean = ymean + ydelta/cast(k,'like',y);
    Syy = Syy + ydelta*(y(k) - ymean);
    costSegment1(k) = Syy;
end
% Cost of all possible second segments:
% same iterative formula as the builtin instead of cpnochange on a(ix:n)
ymean = zeros('like',y);
Syy = ymean;
for k = n:-1:1
    ydelta = y(k) - ymean;
    ymean = ymean + ydelta/cast(n - k + 1,'like',y);
    Syy = Syy + ydelta*(y(k) - ymean);
    costSegment2(k) = Syy;
end

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

function [costSegment1,costSegment2] = costTwoSegmentsVariance(y,n)
% Evaluate cost function for all possible combinations of two segments
ny = coder.internal.indexInt(numel(y));
if coder.internal.isConst(ny)
    costSegment1 = coder.nullcopy(zeros(ny,1,'like',y));
    costSegment2 = coder.nullcopy(zeros(ny,1,'like',y));
    for k = n+1:ny
        costSegment1(k) = 0;
        costSegment2(k) = 0;
    end
else
    assert(n <= ny); %#HINT
    costSegment1 = coder.nullcopy(zeros(n,1,'like',y));
    costSegment2 = coder.nullcopy(zeros(n,1,'like',y));
end
LOGREALMIN = log(realmin(class(y))); % MATLAB has log(realmin)
% Cost of all possible first segments:
% same iterative formula as the builtin instead of cpnochange on a(1:ix)
ymean = zeros('like',y);
Syy = ymean;
for k = 1:n
    npoints = cast(k,'like',y);
    ydelta = y(k) - ymean;
    ymean = ymean + ydelta/npoints;
    Syy = Syy + ydelta*(y(k) - ymean);
    costSegment1(k) = npoints* ...
        max(LOGREALMIN - log(npoints),log(Syy/npoints));
end
% Cost of all possible second segments:
% same iterative formula as the builtin instead of cpnochange on a(ix:n)
ymean = zeros('like',y);
Syy = ymean;
for k = n:-1:1
    ydelta = y(k) - ymean;
    npoints = cast(n + 1 - k,'like',y);
    ymean = ymean + ydelta/npoints;
    Syy = Syy + ydelta*(y(k) - ymean);
    costSegment2(k) = npoints* ...
        max(LOGREALMIN - log(npoints),log(Syy/npoints));
end

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

function [costSegment1,costSegment2] = costTwoSegmentsLinear(y,x,n)
% Evaluate cost function for all possible combinations of two segments
yZERO = zeros('like',y);
ny = coder.internal.indexInt(numel(y));
if coder.internal.isConst(ny)
    costSegment1 = coder.nullcopy(zeros(ny,1,'like',y));
    costSegment2 = coder.nullcopy(zeros(ny,1,'like',y));
    for k = n+1:ny
        costSegment1(k) = 0;
        costSegment2(k) = 0;
    end
else
    assert(n <= ny); %#HINT
    costSegment1 = coder.nullcopy(zeros(n,1,'like',y));
    costSegment2 = coder.nullcopy(zeros(n,1,'like',y));
end
if n == 0
    return
end
% Cost of all possible first segments:
% same iterative formula as the builtin instead of cpnochange on a(1:ix)
xZERO = zeros('like',x);
Sxx = xZERO;
Syy = yZERO;
Sxy = yZERO;
SxxSSE = yZERO;
if isa(x,'double')
    npType = 0;
else
    npType = yZERO;
end
ymean = y(1);
xmean = x(1);
costSegment1(1) = coder.internal.nan('like',yZERO);
for k = 2:n
    npoints = cast(k,'like',npType);
    ydelta = y(k) - ymean;
    xdelta = x(k) - xmean;
    ymean(1) = ymean + ydelta/npoints;
    xmean(1) = xmean + xdelta/npoints;
    dSyy = ydelta*(y(k) - ymean);
    dSxx = xdelta*(x(k) - xmean);
    c = (npoints - 1)/npoints;
    dSxy = xdelta*ydelta*c;
    Syy = Syy + dSyy;
    dSxxSSE = dSxx*Syy + dSyy*Sxx - dSxy*(2*Sxy + dSxy);
    Sxx = Sxx + dSxx;
    Sxy = Sxy + dSxy;
    SxxSSE = SxxSSE + dSxxSSE;
    costSegment1(k) = SxxSSE/Sxx;
end
% Cost of all possible second segments:
% same iterative formula as the builtin instead of cpnochange on a(ix:n)
Sxx = xZERO;
Syy = yZERO;
Sxy = yZERO;
ymean = y(n);
xmean = x(n);
costSegment2(n) = coder.internal.nan('like',yZERO);
for k = n-1:-1:1
    npoints = cast(n + 1 - k,'like',npType);
    ydelta = y(k) - ymean;
    xdelta = x(k) - xmean;
    ymean = ymean + ydelta/npoints;
    xmean = xmean + xdelta/npoints;
    Syy = Syy + ydelta*(y(k) - ymean);
    Sxx = Sxx + xdelta*(x(k) - xmean);
    c = (npoints - 1)/npoints;
    Sxy = Sxy + xdelta*ydelta*c;
    costSegment2(k) = Syy - Sxy^2/Sxx;
end

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

function [cp,ncp] = cpmulti(cp,a,n,method,Lmin,Kmax,x)
% Multiple changepoints using 'MaxNumChanges' = Kmax
% Try one change
AZERO = zeros('like',a);
AEPS = eps(class(a));
% Try one change
[cost1,emptyFlag] = cpsingle(a,n,method,Lmin,x);
if emptyFlag
    ncp = coder.internal.indexInt(0);
    return
end
% Cost lower bound: overfit with penalty set to zero
[cp,ncp,mincost] = ischangePelt(cp,a,n,method,0,x,Lmin,true);
if Kmax >= ncp
    return
end
% Cost upper bound: cost function if data had no changepoints
maxcost = cpnochange(a,n,method,x);
maxpenalty = max(maxcost - mincost,AEPS);
% Early return if the cost of having no changepoints is the same as the
% cost of having the largest possible number of changepoints (for example,
% data on an almost straight line for 'linear'). 
tol = 4*max(AEPS,eps(maxpenalty));
if ~(maxpenalty > tol)
    ncp = coder.internal.indexInt(0);
    return
end
% Initial penalty: cost1 + penalty = maxcost
penalty = max(maxcost - cost1,AEPS);
[cp,ncp] = ischangePelt(cp,a,n,method,penalty,x,Lmin,true);
% Seek lower bound:
% decrease the penalty until it produces more changepoints than Kmax
minpenalty = AZERO;
cpmax = cp;
ncpmax = ncp;
seekUpper = true;
while ncp < Kmax && abs(penalty - minpenalty) >= tol
    seekUpper = false;
    maxpenalty = penalty;
    ncpmax = ncp;
    for k = 1:ncpmax
        cpmax(k) = cp(k);
    end
    penalty = penalty/2;
    [cp,ncp] = ischangePelt(cp,a,n,method,penalty,x,Lmin,true);
    if ncp > Kmax
        minpenalty = penalty;
    end
end
% Seek upper bound:
% increase the penalty until it produces less changepoints than Kmax
if seekUpper
    while ncp > Kmax && abs(penalty - maxpenalty) >= tol
        minpenalty = penalty;
        penalty = 2*penalty;
        [cp,ncp] = ischangePelt(cp,a,n,method,penalty,x,Lmin,true);
        if ncp < Kmax
            maxpenalty = penalty;
            ncpmax = ncp;
            for k = 1:ncpmax
                cpmax(k) = cp(k);
            end
        end
    end
end
if ncp ~= Kmax
    % Search for specified number of changes
    penalty = (maxpenalty + minpenalty)/2;
    while ncp ~= Kmax && abs(minpenalty - penalty) > tol && abs(penalty - maxpenalty) > tol
        [cp,ncp] = ischangePelt(cp,a,n,method,penalty,x,Lmin,true);
        if ncp < Kmax
            maxpenalty = penalty;
            ncpmax = ncp;
            for k = 1:ncpmax
                cpmax(k) = cp(k);
            end
        else
            minpenalty = penalty;
        end
        penalty = (maxpenalty + minpenalty)/2;
    end
    % Don't exceed Kmax when reporting penalties
    if ncp ~= Kmax
        ncp = ncpmax;
        for k = 1:ncp
            cp(k) = cpmax(k);
        end
    end
end

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

function [s1,s2] = fitsegments(a,n,cp,ncp,method,x)
% Useful information about each segment delimited by the changepoints
coder.internal.prefer_const(method);
s1 = coder.nullcopy(a);
s2 = coder.nullcopy(a);
ONE = coder.internal.indexInt(1);
if isempty(a)
    return
end
if method == LINEAR
    fitfun = @(istart,istop)linfit(x,a,istart,istop);
else
    fitfun = @(istart,istop)segmentMeanAndVariance(a,istart,istop,false);
end
if ncp == 0
    [m,c] = fitfun(ONE,n);
    s1(:) = m;
    s2(:) = c;
    return
end
% The calculations over each segment are independent and could be
% performed in parallel.
% First segment
istop = cp(1) - 1;
[m,c] = fitfun(ONE,istop);
for j = ONE:istop
    s1(j) = m;
    s2(j) = c;
end
% Middle segments
for k = 2:ncp
    istart = cp(k - 1);
    istop = cp(k) - 1;
    [m,c] = fitfun(istart,istop);
    for j = istart:istop
        s1(j) = m;
        s2(j) = c;
    end
end
% Last segment
istart = cp(ncp);
[m,c] = fitfun(istart,n);
for j = istart:n
    s1(j) = m;
    s2(j) = c;
end

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

function [slope,yintercept] = linfit(x,y,istart,istop)
coder.internal.assert(coder.const(~coder.internal.eml_option_eq('VariableSizing','DisableInInference')), ...
    'Coder:toolbox:IschangeVarSizeNeeded');
n = istop - istart + 1;
assert(n <= numel(y)); %<HINT>
V = ones(n,2,'like',y);
b = zeros(n,1,'like',y);
for j = istart:istop
    V(j - istart + 1,1) = x(j);
    b(j - istart + 1) = y(j);
end
p = coder.internal.qrsolve(V,b,'FullRank','NoWarn');
slope = p(1);
yintercept = p(2);

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

function s = vsubsum(x,ilo,ihi)
% xs = sum(x(begSeg:endSeg)) with blocking for single precision.
s = zeros(class(x));
if coder.internal.isConstTrue(isempty(x))
    % s = cast(sum(ilo:ihi),'like',x)
    % Assume ihi >= ilo - 1.
    s = cast((ihi - ilo + 1)*(ilo + ihi)/2,'like',x);
elseif isempty(x)
elseif issparse(x)
    s = full(sum(x(ilo:ihi)));
elseif isa(x,'double')
    s = sumab(x,ilo,ihi);
else
    % Blocked summation.
    BLOCKLENGTH = coder.internal.indexInt(1024);
    nxsub = ihi - ilo + 1;
    if nxsub > BLOCKLENGTH
        nblocks = coder.internal.indexDivide(nxsub,BLOCKLENGTH);
        for b = 1:nblocks
            begBlock = ilo + (b - 1)*BLOCKLENGTH;
            endBlock = begBlock + BLOCKLENGTH - 1;
            s = s + sumab(x,begBlock,endBlock);
        end
        begBlock = ilo + nblocks*BLOCKLENGTH;
    else
        begBlock = ilo;
    end
    if ihi >= begBlock
        s = s + sumab(x,begBlock,ihi);
    end
end

function s = sumab(x,a,b)
coder.inline('always');
s = zeros(class(x));
for k = a:b
    s = s + x(k);
end

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

function [m,v] = segmentMeanAndVariance(x,begSeg,endSeg,bias)
if issparse(x)
    seg = x(begSeg:endSeg);
    m = mean(seg);
    v = var(seg);
else
    ONE = coder.internal.indexInt(1);
    nxsub = endSeg - begSeg + 1;
    dnxsub = cast(nxsub,'like',real(x));
    m = vsubsum(x,begSeg,endSeg);
    m = m/dnxsub;
    work = coder.nullcopy(real(x));
    for k = begSeg:endSeg
        work(k) = abs(x(k) - m);
    end
    v = coder.internal.blas.xnrm2(nxsub,work,begSeg,ONE).^2;
    if nxsub >= 2
        if bias
            v = v/dnxsub;
        else
            v = v/(dnxsub - 1);
        end
    end
end

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

function y = LINEAR
coder.inline('always');
y = 'L';

function y = MEAN
coder.inline('always');
y = 'M';

function y = VARIANCE
coder.inline('always');
y = 'V';

function y = UNRECOGNIZED
coder.inline('always');
y = 'U';

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

function checkSamplePoints(x,A,dim)
% Validate SamplePoints value.
coder.internal.prefer_const(dim);
coder.internal.assert(isfloat(x) || isduration(x) || isdatetime(x), ...
    'MATLAB:ischange:SamplePointsInvalidDatatype');
EMPTYX = isempty(x);
EMPTYA = isempty(A);
n = coder.internal.indexInt(numel(x));
if isa(A,'timetable')
    MUSTBEFINITEMSG = 'MATLAB:ischange:RowTimesFinite';
    MUSTBESORTEDMSG = 'MATLAB:ischange:RowTimesSorted';
    HASDUPLICATEMSG = 'MATLAB:ischange:RowTimesDuplicate';
else
    coder.internal.assert(isvector(x) || EMPTYX, ...
        'MATLAB:ischange:SamplePointsInvalidDatatype');
    coder.internal.assert(n == size(A,dim) || (EMPTYX && EMPTYA), ...
        'MATLAB:ischange:SamplePointsLength');
    if isfloat(x)
        coder.internal.assert(isreal(x), ...
            'MATLAB:ischange:SamplePointsComplex');
        coder.internal.assert(~issparse(x), ...
            'MATLAB:ischange:SamplePointsSparse');
    end
    MUSTBEFINITEMSG = 'MATLAB:ischange:SamplePointsFinite';
    MUSTBESORTEDMSG = 'MATLAB:ischange:SamplePointsSorted';
    HASDUPLICATEMSG = 'MATLAB:ischange:SamplePointsDuplicate';
end
coder.internal.errorIf(coder.internal.anyNonFinite(x),MUSTBEFINITEMSG);
for k = 2:n
    if x(k) <= x(k - 1)
        coder.internal.assert(x(k) == x(k - 1),MUSTBESORTEDMSG);
        coder.internal.assert(x(k) ~= x(k - 1),HASDUPLICATEMSG);
    end
end

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

function n = preferConstN(x,y)
XEMPTY = coder.internal.isConstTrue(isempty(x));
nx = coder.internal.indexInt(numel(x));
ny = coder.internal.indexInt(numel(y));
CONSTNX = coder.internal.isConst(nx);
CONSTNY = coder.internal.isConst(ny);
if XEMPTY
    n = ny;
elseif CONSTNX && CONSTNY
    n = coder.const(eml_min(nx,ny));
elseif CONSTNX
    n = coder.const(nx);
elseif CONSTNY
    n = coder.const(ny);
else
    n = eml_min(nx,ny);
end

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

function method = parseMethod(s)
coder.internal.prefer_const(s);
coder.internal.assert(coder.internal.isConst(s), ...
    'Coder:toolbox:OptionStringsMustBeConstant');
len = max(1,coder.internal.indexInt(strlength(s)));
if coder.const(strncmpi(s,'mean',len))
    method = MEAN;
elseif strncmpi(s,'variance',len)
    method = VARIANCE;
elseif strncmpi(s,'linear',len)
    method = LINEAR;
else
    method = UNRECOGNIZED;
end

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