## Copyright (C) 2001-2012 Kai Habel
##
## 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{pp} =} pchip (@var{x}, @var{y})
## @deftypefnx {Function File} {@var{yi} =} pchip (@var{x}, @var{y}, @var{xi})
## Return the Piecewise Cubic Hermite Interpolating Polynomial (pchip) of
## points @var{x} and @var{y}.
##
## If called with two arguments, return the piecewise polynomial @var{pp}
## that may be used with @code{ppval} to evaluate the polynomial at specific
## points.  When called with a third input argument, @code{pchip} evaluates
## the pchip polynomial at the points @var{xi}.  The third calling form is
## equivalent to @code{ppval (pchip (@var{x}, @var{y}), @var{xi})}.
##
## The variable @var{x} must be a strictly monotonic vector (either
## increasing or decreasing) of length @var{n}.  @var{y} can be either a
## vector or array.  If @var{y} is a vector then it must be the same length
## @var{n} as @var{x}.  If @var{y} is an array then the size of @var{y} must
## have the form
## @tex
## $$[s_1, s_2, \cdots, s_k, n]$$
## @end tex
## @ifnottex
## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n}]}
## @end ifnottex
## The array is reshaped internally to a matrix where the leading
## dimension is given by
## @tex
## $$s_1 s_2 \cdots s_k$$
## @end tex
## @ifnottex
## @code{@var{s1} * @var{s2} * @dots{} * @var{sk}}
## @end ifnottex
## and each row of this matrix is then treated separately.  Note that this
## is exactly opposite to @code{interp1} but is done for @sc{matlab}
## compatibility.
##
## @seealso{spline, ppval, mkpp, unmkpp}
## @end deftypefn

## Author:  Kai Habel <kai.habel@gmx.de>
## Date: 9. mar 2001
##
## S_k = a_k + b_k*x + c_k*x^2 + d_k*x^3; (spline polynom)
##
## 4 conditions:
## S_k(x_k) = y_k;
## S_k(x_k+1) = y_k+1;
## S_k'(x_k) = y_k';
## S_k'(x_k+1) = y_k+1';

function ret = pchip (x, y, xi)

  if (nargin < 2 || nargin > 3)
    print_usage ();
  endif

  ## make row vector
  x = x(:).';
  n = length (x);

  ## Check the size and shape of y
  if (isvector (y))
    y = y(:).'; ##row vector
    szy = size (y);
    if !(size_equal (x, y))
      error ("pchip: length of X and Y must match")
    endif
  else
    szy = size (y);
    if (n != szy(end))
      error ("pchip: length of X and last dimension of Y must match")
    endif
    y = reshape (y, [prod(szy(1:end-1)), szy(end)]);
  endif

  h = diff (x);
  if (all (h < 0))
    x = fliplr (x);
    h = diff (x);
    y = fliplr (y);
  elseif (any (h <= 0))
    error("pchip: X must be strictly monotonic");
  endif

  f1 = y(:, 1:n-1);

  ## Compute derivatives.
  d = __pchip_deriv__ (x, y, 2);
  d1 = d(:, 1:n-1);
  d2 = d(:, 2:n);

  ## This is taken from SLATEC.
  h = diag (h);

  delta = diff (y, 1, 2) / h;
  del1 = (d1 - delta) / h;
  del2 = (d2 - delta) / h;
  c3 = del1 + del2;
  c2 = -c3 - del1;
  c3 = c3 / h;
  coeffs = cat (3, c3, c2, d1, f1);

  ret = mkpp (x, coeffs, szy(1:end-1));

  if (nargin == 3)
    ret = ppval (ret, xi);
  endif

endfunction

%!demo
%! x = 0:8;
%! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0];
%! xi = 0:0.01:8;
%! yspline = spline(x,y,xi);
%! ypchip = pchip(x,y,xi);
%! title("pchip and spline fit to discontinuous function");
%! plot(xi,yspline,xi,ypchip,"-",x,y,"+");
%! legend ("spline","pchip","data");
%! %-------------------------------------------------------------------
%! % confirm that pchip agreed better to discontinuous data than spline

%!shared x,y,y2,pp,yi1,yi2,yi3
%! x = 0:8;
%! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0];
%!assert (pchip(x,y,x), y);
%!assert (pchip(x,y,x'), y');
%!assert (pchip(x',y',x'), y');
%!assert (pchip(x',y',x), y);
%!assert (isempty(pchip(x',y',[])));
%!assert (isempty(pchip(x,y,[])));
%!assert (pchip(x,[y;y],x), [pchip(x,y,x);pchip(x,y,x)])
%!assert (pchip(x,[y;y],x'), [pchip(x,y,x);pchip(x,y,x)])
%!assert (pchip(x',[y;y],x), [pchip(x,y,x);pchip(x,y,x)])
%!assert (pchip(x',[y;y],x'), [pchip(x,y,x);pchip(x,y,x)])
%!test
%! x=(0:8)*pi/4;y=[sin(x);cos(x)];
%! y2(:,:,1)=y;y2(:,:,2)=y+1;y2(:,:,3)=y-1;
%! pp=pchip(x,shiftdim(y2,2));
%! yi1=ppval(pp,(1:4)*pi/4);
%! yi2=ppval(pp,repmat((1:4)*pi/4,[5,1]));
%! yi3=ppval(pp,[pi/2,pi]);
%!assert(size(pp.coefs),[48,4]);
%!assert(pp.pieces,8);
%!assert(pp.order,4);
%!assert(pp.dim,[3,2]);
%!assert(ppval(pp,pi),[0,-1;1,0;-1,-2],1e-14);
%!assert(yi3(:,:,2),ppval(pp,pi),1e-14);
%!assert(yi3(:,:,1),[1,0;2,1;0,-1],1e-14);
%!assert(squeeze(yi1(1,2,:)),[1/sqrt(2); 0; -1/sqrt(2);-1],1e-14);
%!assert(size(yi2),[3,2,5,4]);
%!assert(squeeze(yi2(1,2,3,:)),[1/sqrt(2); 0; -1/sqrt(2);-1],1e-14);