function v = ppValDiag(pp,xx)
%PPVALDIAG Evaluate a set of piecewise polynomials using each input only
%once.
%   V = PPVALDIAG(PP, XX) returns a value, for each entry in XX, of the
%   piecewise polynomial f contained in PP, as constructed by the spline
%   utility MKPP in INTERPSPLINE.
%
%   V is obtained by replacing each entry of XX by the value of f there
%   using a new piecewise polynomial for each evaluation.  It is required
%   that the dimension of the piecewise polynomial, pp.dim, be equal to the
%   number of elements in xx, which must be a vector.
%
%   Class support for input XX and the fields of PP:
%      float: double, single
%
%   See also PPVAL, CODISTRIBUTED/PRIVATE/INTERPSPLINE.

%   Copyright 2022 The Mathworks, Inc.

if ~isfloat(xx)
    error(message('MATLAB:chckxy:XQType'))
end
if ~isreal(xx)
    error(message('MATLAB:chckxy:XQComplex'))
end

if ~isrow(xx)
    assert(isvector(xx), 'The second input to PPVALDIAG must be a vector.')
    xx = xx';
end

lx = numel(xx);

% Take apart pp.
[b, c, l, k, dd] = unmkpp(pp);

assert(lx == dd, 'The length of the query input to PPVALDIAG must be equal to the dimension of the piecewise polynomial.')

% For each evaluation site, compute its breakpoint interval
% (mindful of the possibility that xx might be empty).
if lx
    [~, ~, index] = histcounts(xx, [-inf, b(2:l), inf]);
else
    index = ones(1, lx);
end

% Adjust for troubles, like evaluation sites that are NaN or +-inf.
infxx = (xx==inf);
if any(infxx)
    index(infxx) = l;
end
nogoodxx = (index==0);
if any(nogoodxx)
    xx(nogoodxx) = NaN;
    index(nogoodxx) = 1;
end

% Now go to local coordiates ...
xx = xx - b(index);
index = (1:dd) + dd*(index-1);

% ... and apply nested multiplication:
v = c(index, 1);
for ii = 2:k
    v = xx(:).*v + c(index, ii);
end
end