function [cp,ncp,optcost] = ischangePelt(cp,y,n,method,penalty, ...
    x,minsep,dopruning)
%MATLAB Code Generation Private Function

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


% x can be constant [] when the method is Variance or Mean.
coder.internal.prefer_const(n,method,dopruning);
ONE = coder.internal.indexInt(1);
yZERO = coder.internal.scalarEg(y);
methodstr = coder.const(validatestring(method, ...
    {'Linear','Variance','Mean'},'method'));
stat = coder.const(methodstr(1)); % L(inear), V(ariance), or M(ean).
nmax = coder.internal.indexInt(length(cp) - 1);
work = initWork(y,nmax,stat);
[ibest,pbest] = work.plist.front;
work.costs(pbest.istart) = -penalty;
work.dists(pbest.istop) = 0;
ix = ONE;
if n > 0
    cp(1) = ONE;
end
% initialize changepoint with distance to first legal changepoint
while ix < minsep && ix <= n
    work = updateDistsWithNewEndpoint(work,stat,getx(x,ix),y(ix),ix);
    cp(ix + 1) = ONE;
    ix = ix + 1;
end
while ix <= n
    % update all active changepoints with distance to new endpoint
    work = updateDistsWithNewEndpoint(work,stat,getx(x,ix),y(ix),ix);
    % add new point
    if ix > minsep
        istop = ix  + 1 - minsep;
        work.plist = work.plist.pushBack(makePartition(istop,istop));
        [d,work] = getDistToEnd(work,stat,x,y,istop,ix);
        work.dists(istop) = d;
    end
    % update ibest, splitting partitions if necessary.
    [ibest,best,work] = getBestPartition(work,stat,x,y,ix,ibest);
    % store result
    pbest = work.plist.getValue(ibest);
    cp(ix + 1) = pbest.istart;
    cost = best + penalty;
    work.costs(ix + 1) = cost;
    if dopruning
        work = prunePartitions(work,stat,x,y,ix,cost);
    end
    if needsMerge(stat)
        work = mergePartitions(work,best,ibest);
    end
    ix = ix + 1;
end
optcost = work.costs(n + 1);
[cp,ncp] = shiftIndices(cp);
optcost = optcost - cast(ncp,'like',yZERO)*cast(penalty,'like',yZERO);

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

function p = needsMerge(stat)
coder.internal.prefer_const(stat);
p = coder.const(stat == 'L' || stat == 'M');

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

function [cpidx,nchanges] = shiftIndices(cpidx)
np1 = coder.internal.indexInt(numel(cpidx));
ileft = np1;
iright = np1 + 1;
while ileft > 1 && cpidx(ileft) > 1 && cpidx(ileft) < ileft
    iright = iright - 1;
    cpidx(iright) = cpidx(ileft);
    ileft = cpidx(ileft);
end
iright = iright - 1;
nchanges = np1 - iright;
for ileft = 1:nchanges
    cpidx(ileft) = cpidx(iright + ileft);
end
for k = nchanges+1:np1
    cpidx(k) = 0;
end

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

function [ibest,best,work] = getBestPartition(work,stat,x,y,ix,ibest)
% initialize best new cost with the (updated) previous best location
pbest = work.plist.getValue(ibest);
best = work.costs(pbest.istart) + work.dists(pbest.istop);
[ip,part] = work.plist.front;
while ip ~= 0
    while work.costs(part.istart) + work.dists(part.istop) < best
        istart = part.istart;
        istop = part.istop;
        startcost = work.costs(istart);
        stopdist = work.dists(istop);
        if istart == istop
            % no partition needed.  just take the new best value
            best = startcost + stopdist;
            ibest = ip;
        elseif work.costs(istop) + stopdist < best
            % last point is best so far.  partition it separately
            [work.plist,ibest] = work.plist.insertAfter(ip, ...
                makePartition(istop,istop));
            best = work.costs(istop) + stopdist;
            istop = istop - 1;
            part.istop = istop;
            work.plist = work.plist.setValue(ip,part);
            [d,work] = growDist(work,stat,getx(x,istop),y(istop),ix,istop);
            work.dists(istop) = d;
        else
            imiddle = part.istop;
            while true
                imiddle = imiddle - 1;
                [middist,work] = growDist(work,stat, ...
                    getx(x,imiddle),y(imiddle),ix,imiddle);
                if ~( ...
                        startcost + middist < best && ...
                        work.costs(imiddle) <= work.costs(imiddle + 1) && ...
                        work.costs(imiddle) + stopdist >= best && ...
                        imiddle ~= istart)
                    break
                end
            end
            [work.plist,~] = work.plist.insertAfter( ...
                ip,makePartition(imiddle+1,istop));
            part.istop = imiddle;
            work.plist = work.plist.setValue(ip,part);
            work.dists(imiddle) = middist;
        end
    end
    [ip,part] = work.plist.next(ip);
end

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

function work = initWork(y,n,stat)
coder.internal.prefer_const(n,stat);
coder.inline('always');
if coder.internal.isConst(n)
    % Constant allocation.
    np1 = coder.internal.indexInt(n) + 1;
else
    % Potentially bounded allocation.
    np1 = coder.internal.indexInt(numel(y)) + 1;
end
assert(np1 <= numel(y) + 1); %<HINT>
aPart = makePartition(1,1);
yZERO = zeros('like',y);
partitionList = coder.internal.list(aPart, ...
    'InitialCapacity',np1, ...
    'FixedCapacity',true);
partitionList = partitionList.pushBack(aPart);
if coder.const(stat == 'L')
    work = struct( ...
        'plist',partitionList, ...
        'costs',coder.internal.inf(np1,1,'like',yZERO), ...
        'dists',zeros(np1,1,'like',yZERO), ...
        'xOrigin',zeros(np1,1,'like',yZERO), ...
        'xBar',zeros(np1,1,'like',yZERO), ...
        'yBar',zeros(np1,1,'like',yZERO), ...
        'Sxx',zeros(np1,1,'like',yZERO), ...
        'Sxy',zeros(np1,1,'like',yZERO), ...
        'Syy',zeros(np1,1,'like',yZERO), ...
        'SxxSSE',zeros(np1,1,'like',yZERO));
elseif coder.const(stat == 'V')
    work = struct( ...
        'plist',partitionList, ...
        'costs',coder.internal.inf(np1,1,'like',yZERO), ...
        'dists',zeros(np1,1,'like',yZERO), ...
        'Syy',zeros(np1,1,'like',yZERO), ...
        'yBar',zeros(np1,1,'like',yZERO));
else % if stat == 'M'
    work = struct( ...
        'plist',partitionList, ...
        'costs',coder.internal.inf(np1,1,'like',yZERO), ...
        'dists',zeros(np1,1,'like',yZERO), ...
        'yBar',zeros(np1,1,'like',yZERO));
end

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

function work = updateDistWithNewEndpoint(work,stat,x,y,i,iend)
npoints = 1 + iend - i;
ydelta = y - work.yBar(i);
if coder.const(stat == 'L')
    if isa(x,'double')
        np = double(npoints);
    else
        np = cast(npoints,'like',y);
    end
    x = x - work.xOrigin(i);
    xdelta = x - work.xBar(i);
    work.xBar(i) = work.xBar(i) + xdelta/np;
    work.yBar(i) = work.yBar(i) + ydelta/np;
    dSyy = ydelta*(y - work.yBar(i));
    dSxx = xdelta*(x - work.xBar(i));
    dSxy = xdelta*ydelta*(np - 1)/np;
    work.Syy(i) = work.Syy(i) + dSyy;
    dSxxSSE = dSxx*work.Syy(i) + dSyy*work.Sxx(i) - ...
        dSxy*(2*work.Sxy(i) + dSxy);
    work.Sxx(i) = work.Sxx(i) + dSxx;
    work.Sxy(i) = work.Sxy(i) + dSxy;
    work.SxxSSE(i) = work.SxxSSE(i) + dSxxSSE;
    work.dists(i) = work.SxxSSE(i)/work.Sxx(i);
elseif coder.const(stat == 'V')
    np = cast(npoints,'like',y);
    work.yBar(i) = work.yBar(i) + ydelta/np;
    work.Syy(i) = work.Syy(i) + ydelta*(y - work.yBar(i));
    work.dists(i) = getVarDist(work.Syy(i),np);
else % if stat == 'M'
    np = cast(npoints,'like',y);
    work.yBar(i) = work.yBar(i) + ydelta/np;
    work.dists(i) = work.dists(i) + ydelta*(y - work.yBar(i));
end

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

function d = getVarDist(syy,npoints)
% std::max(
%   npoints*(logrealmin - log(npoints)),
%   npoints*log(syy/npoints)
%);
np = cast(npoints,'like',syy);
LOGREALMIN = log(realmin(class(syy)));
x1 = np*(LOGREALMIN - log(np));
x2 = np*log(syy/np);
d = eml_max(x1,x2);

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

function work = updateDistsWithNewEndpoint(work,stat,x,y,iend)
[ip,part] = work.plist.front;
while ip ~= 0
    work = updateDistWithNewEndpoint(work,stat,x,y,part.istop,iend);
    [ip,part] = work.plist.next(ip);
end

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

function [d,work] = getDistToEnd(work,stat,x,y,istart,istop)
ZERO = coder.internal.indexInt(0);
xZERO = zeros('like',x);
yZERO = zeros('like',y);
d = yZERO;
if coder.const(stat == 'L')
    if isa(x,'double')
        npType = xZERO;
    else
        npType = yZERO;
    end
    npoints = ZERO;
    xmean = xZERO;
    ymean = yZERO;
    sxx = xZERO;
    sxy = yZERO;
    syy = yZERO;
    sxxsse = yZERO;
    xorigin = x(istop);
    for i = istart:istop
        xi = x(i) - xorigin;
        xdelta = xi - xmean;
        ydelta = y(i) - ymean;
        npoints = npoints + 1;
        np = cast(npoints,'like',npType);
        xmean(1) = xmean + xdelta/np;
        ymean(1) = ymean + ydelta/np;
        dSyy = ydelta*(y(i) - ymean);
        dSxx = xdelta*(xi - xmean);
        dSxy = xdelta*ydelta*(np - 1)/np;
        syy(1) = syy + dSyy;
        dSxxSSE = dSxx*syy + dSyy*sxx - dSxy*(2*sxy + dSxy);
        sxx(1) = sxx + dSxx;
        sxy(1) = sxy + dSxy;
        sxxsse(1) = sxxsse + dSxxSSE;
    end
    work.xOrigin(istart) = xorigin;
    work.xBar(istart) = xmean;
    work.yBar(istart) = ymean;
    work.Sxx(istart) = sxx;
    work.Sxy(istart) = sxy;
    work.Syy(istart) = syy;
    work.SxxSSE(istart) = sxxsse;
    d(1) = sxxsse/sxx;
    work.dists(istart) = d;
elseif coder.const(stat == 'V')
    npoints = ZERO;
    ymean = yZERO;
    syy = yZERO;
    for i = istart:istop
        ydelta = y(i) - ymean;
        npoints = npoints + 1;
        ymean(1) = ymean + ydelta/cast(npoints,'like',yZERO);
        syy(1) = syy + ydelta*(y(i) - ymean);
    end
    work.yBar(istart) = ymean;
    work.Syy(istart) = syy;
    d(1) = getVarDist(syy,npoints);
    work.dists(istart) = d;
else % if stat == 'M'
    npoints = ZERO;
    ymean = yZERO;
    for i = istart:istop
        ydelta = y(i) - ymean;
        npoints = npoints + 1;
        ymean(1) = ymean + ydelta/cast(npoints,'like',yZERO);
        d(1) = d + ydelta*(y(i) - ymean);
    end
    work.yBar(istart) = ymean;
    work.dists(istart) = d;
end

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

function [d,work] = growDist(work,stat,x,y,ix,istart)
coder.internal.prefer_const(stat);
inext = istart + 1;
ymean = work.yBar(inext);
ydelta = y - ymean;
npoints = ix - istart + 1;
if coder.const(stat == 'L')
    if isa(x,'double')
        np = double(npoints);
    else
        np = cast(npoints,'like',y);
    end
    ymean = ymean + ydelta/np;
    work.yBar(istart) = ymean;
    x = x - work.xOrigin(inext);
    xmean = work.xBar(inext);
    xdelta = x - xmean;
    xmean = xmean + xdelta/np;
    dSyy = ydelta*(y - ymean);
    dSxx = xdelta*(x - xmean);
    dSxy = xdelta*ydelta*(np - 1)/np;
    work.Syy(istart) = work.Syy(inext) + dSyy;
    dSxxSSE = dSxx*work.Syy(istart) + dSyy*work.Sxx(inext) - ...
        dSxy*(2*work.Sxy(inext) + dSxy);
    work.xOrigin(istart) = work.xOrigin(inext);
    work.xBar(istart) = xmean;
    work.Sxx(istart) = work.Sxx(inext) + dSxx;
    work.Sxy(istart) = work.Sxy(inext) + dSxy;
    work.SxxSSE(istart) = work.SxxSSE(inext) + dSxxSSE;
    d = work.SxxSSE(istart)/work.Sxx(istart);
elseif coder.const(stat == 'V')
    np = cast(ix - istart + 1,'like',y);
    ymean = ymean + ydelta/np;
    work.yBar(istart) = ymean;
    work.Syy(istart) = work.Syy(inext) + ydelta*(y - ymean);
    d = getVarDist(work.Syy(istart),np);
else % if stat == 'M'
    np = cast(ix - istart + 1,'like',y);
    ymean = ymean + ydelta/np;
    work.yBar(istart) = ymean;
    d = work.dists(inext) + ydelta*(y - ymean);
end
work.dists(istart) = d;

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

function p = makePartition(istart,istop)
p = struct( ...
    'istart',coder.internal.indexInt(istart), ...
    'istop',coder.internal.indexInt(istop));

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

function [pruned,work,ip,part] = prunePartition(work,ip,part,thresh)
pruned = work.costs(part.istart) + work.dists(part.istop) > thresh;
if pruned
    ip0 = ip;
    [ip,part] = work.plist.next(ip);
    work.plist = work.plist.erase(ip0);
end

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

function work = prunePartitions(work,stat,x,y,ix,thresh)
coder.internal.prefer_const(stat);
[ip,part] = work.plist.front;
while ip ~= 0
    [pruned,work,ip,part] = prunePartition(work,ip,part,thresh);
    % If pruned == true, ip has been incremented.
    if ~pruned
        istart = part.istart;
        istop = part.istop;
        startcost = work.costs(istart);
        stopdist = work.dists(istop);
        while istart < istop && ...
                startcost + stopdist < thresh && ...
                work.costs(istop) + stopdist > thresh
            istop = istop - 1;
            [stopdist,work] = growDist(work,stat, ...
                getx(x,istop),y(istop),ix,istop);
        end
        part.istop = istop;
        work.plist = work.plist.setValue(ip,part);
        [pruned,work,ip,part] = prunePartition(work,ip,part,thresh);
        if ~pruned
            [ip,part] = work.plist.next(ip);
        end
    end
end

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

function work = mergePartitions(work,best,ibest)
[ileft,pleft] = work.plist.front;
if ileft ~= 0
    [imid,pmid] = work.plist.next(ileft);
    if imid ~= 0
        [iright,pright] = work.plist.next(imid);
        while iright ~= 0
            if  imid ~= ibest && ileft ~= ibest && iright ~= ibest && ...
                    work.costs(pleft.istart) + work.dists(pright.istop) >= best && ...
                    work.costs(pleft.istart) <= work.costs(pmid.istart) && ...
                    work.costs(pmid.istart) <= work.costs(pright.istart)
                if work.costs(pleft.istart) + work.dists(pleft.istop) > ...
                        work.costs(pright.istart) + work.dists(pright.istop)
                    pleft.istop = pmid.istop;
                    work.plist = work.plist.setValue(ileft,pleft);
                else
                    pright.istart = pmid.istart;
                    work.plist = work.plist.setValue(iright,pright);
                end
                work.plist = work.plist.erase(imid);
                imid = iright;
                pmid = pright;
                [iright,pright] = work.plist.next(iright);
            else
                [ileft,pleft] = work.plist.next(ileft);
                [imid,pmid] = work.plist.next(imid);
                [iright,pright] = work.plist.next(iright);
            end
        end
    end
end

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

function xk = getx(x,k)
coder.inline('always');
if coder.internal.isConstTrue(isempty(x))
    xk = cast(k,'like',x);
else
    xk = x(k);
end

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

% LocalWords:  inear ariance ean changepoint changepoints ibest plist dists Sxx Sxy Syy npoints
% LocalWords:  logrealmin syy istart istop ip
