## Copyright (C) 2000-2013 Paul Kienzle
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} {@var{yi} =} ppval (@var{pp}, @var{xi})
## Evaluate the piecewise polynomial structure @var{pp} at the points @var{xi}.
## If @var{pp} describes a scalar polynomial function, the result is an
## array of the same shape as @var{xi}.
## Otherwise, the size of the result is @code{[pp.dim, length(@var{xi})]} if
## @var{xi} is a vector, or @code{[pp.dim, size(@var{xi})]} if it is a
## multi-dimensional array.
## @seealso{mkpp, unmkpp, spline, pchip}
## @end deftypefn

function yi = ppval (pp, xi)

  if (nargin != 2)
    print_usage ();
  endif
  if (! (isstruct (pp) && strcmp (pp.form, "pp")))
    error ("ppval: first argument must be a pp-form structure");
  endif

  ## Extract info.
  [x, P, n, k, d] = unmkpp (pp);

  ## dimension checks
  sxi = size (xi);
  if (isvector (xi))
    xi = xi(:).';
  endif

  nd = length (d);

  ## Determine intervals.
  xn = numel (xi);
  idx = lookup (x, xi, "lr");

  P = reshape (P, [d, n * k]);
  P = shiftdim (P, nd);
  P = reshape (P, [n, k, d]);
  Pidx = P(idx(:), :);#2d matrix size x: coefs*prod(d) y: prod(sxi)

  if (isvector (xi))
    Pidx = reshape (Pidx, [xn, k, d]);
    Pidx = shiftdim (Pidx, 1);
    dimvec = [d, xn];
  else
    Pidx = reshape (Pidx, [sxi, k, d]);
    Pidx = shiftdim (Pidx, length (sxi));
    dimvec = [d, sxi];
  endif
  ndv = length (dimvec);

  ## Offsets.
  dx = (xi - x(idx))(:)';
  dx = repmat (dx, [prod(d), 1]);
  dx = reshape (dx, dimvec);
  dx = shiftdim (dx, ndv - 1);

  ## Use Horner scheme.
  if (k > 1)
    yi = shiftdim (reshape (Pidx(1,:), dimvec), ndv - 1);
  else
    yi = shiftdim (reshape (Pidx, dimvec), ndv - 1);
  endif

  for i = 2 : k;
    yi .*= dx;
    yi += shiftdim (reshape (Pidx(i,:), dimvec), ndv - 1);
  endfor

  ## Adjust shape.
  if ((numel (xi) > 1) || (length (d) == 1))
    yi = reshape (shiftdim (yi, 1), dimvec);
  endif

  if (isvector (xi) && (d == 1))
    yi = reshape (yi, sxi);
  elseif (isfield (pp, "orient") && strcmp (pp.orient, "first"))
    yi = shiftdim (yi, nd);
  endif

  ## FIXME: Why is this commented out, rather than just removed?
  #if (d == 1)
  #  yi = reshape (yi, sxi);
  #endif

endfunction


%!shared b, c, pp, pp2, xi, abserr
%! b = 1:3;
%! c = ones (2);
%! pp = mkpp (b, c);
%! abserr = 1e-14;
%! pp2 = mkpp (b, [c;c], 2);
%! xi = [1.1 1.3 1.9 2.1];
%!assert (ppval (pp, 1.1), 1.1, abserr)
%!assert (ppval (pp, 2.1), 1.1, abserr)
%!assert (ppval (pp, xi), [1.1 1.3 1.9 1.1], abserr)
%!assert (ppval (pp, xi.'), [1.1 1.3 1.9 1.1].', abserr)
%!assert (ppval (pp2, 1.1), [1.1;1.1], abserr)
%!assert (ppval (pp2, 2.1), [1.1;1.1], abserr)
%!assert (ppval (pp2, xi), [1.1 1.3 1.9 1.1;1.1 1.3 1.9 1.1], abserr)
%!assert (ppval (pp2, xi'), [1.1 1.3 1.9 1.1;1.1 1.3 1.9 1.1], abserr)
%!assert (size (ppval (pp2, [xi;xi])), [2 2 4])
%!test
%! breaks = [0, 1, 2, 3];
%! coefs = rand (6, 4);
%! pp = mkpp (breaks, coefs, 2);
%! ret = zeros (2, 4, 2);
%! ret(:,:,1) = ppval (pp, breaks');
%! ret(:,:,2) = ppval (pp, breaks');
%! assert (ppval (pp, [breaks',breaks']), ret)