1########################################################################
2##
3## Copyright (C) 2001-2021 The Octave Project Developers
4##
5## See the file COPYRIGHT.md in the top-level directory of this
6## distribution or <https://octave.org/copyright/>.
7##
8## This file is part of Octave.
9##
10## Octave is free software: you can redistribute it and/or modify it
11## under the terms of the GNU General Public License as published by
12## the Free Software Foundation, either version 3 of the License, or
13## (at your option) any later version.
14##
15## Octave is distributed in the hope that it will be useful, but
16## WITHOUT ANY WARRANTY; without even the implied warranty of
17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18## GNU General Public License for more details.
19##
20## You should have received a copy of the GNU General Public License
21## along with Octave; see the file COPYING.  If not, see
22## <https://www.gnu.org/licenses/>.
23##
24########################################################################
25
26## -*- texinfo -*-
27## @deftypefn  {} {@var{pp} =} pchip (@var{x}, @var{y})
28## @deftypefnx {} {@var{yi} =} pchip (@var{x}, @var{y}, @var{xi})
29## Return the Piecewise Cubic Hermite Interpolating Polynomial (pchip) of
30## points @var{x} and @var{y}.
31##
32## If called with two arguments, return the piecewise polynomial @var{pp}
33## that may be used with @code{ppval} to evaluate the polynomial at specific
34## points.
35##
36## When called with a third input argument, @code{pchip} evaluates the pchip
37## polynomial at the points @var{xi}.  The third calling form is equivalent to
38## @code{ppval (pchip (@var{x}, @var{y}), @var{xi})}.
39##
40## The variable @var{x} must be a strictly monotonic vector (either increasing
41## or decreasing) of length @var{n}.
42##
43## @var{y} can be either a vector or array.  If @var{y} is a vector then it
44## must be the same length @var{n} as @var{x}.  If @var{y} is an array then
45## the size of @var{y} must have the form
46## @tex
47## $$[s_1, s_2, \cdots, s_k, n]$$
48## @end tex
49## @ifnottex
50## @code{[@var{s1}, @var{s2}, @dots{}, @var{sk}, @var{n}]}
51## @end ifnottex
52## The array is reshaped internally to a matrix where the leading dimension is
53## given by
54## @tex
55## $$s_1 s_2 \cdots s_k$$
56## @end tex
57## @ifnottex
58## @code{@var{s1} * @var{s2} * @dots{} * @var{sk}}
59## @end ifnottex
60## and each row of this matrix is then treated separately.  Note that this is
61## exactly opposite to @code{interp1} but is done for @sc{matlab}
62## compatibility.
63##
64## @seealso{spline, ppval, mkpp, unmkpp}
65## @end deftypefn
66
67## Algorithm:
68## S_k = a_k + b_k*x + c_k*x^2 + d_k*x^3; (spline polynomial)
69##
70## 4 conditions:
71## S_k(x_k) = y_k;
72## S_k(x_k+1) = y_k+1;
73## S_k'(x_k) = y_k';
74## S_k'(x_k+1) = y_k+1';
75
76function ret = pchip (x, y, xi)
77
78  if (nargin < 2 || nargin > 3)
79    print_usage ();
80  endif
81
82  ## make row vector
83  x = x(:).';
84  n = length (x);
85
86  ## Check the size and shape of y
87  if (isvector (y))
88    y = y(:).';  # force row vector
89    szy = size (y);
90    if (! size_equal (x, y))
91      error ("pchip: length of X and Y must match");
92    endif
93  else
94    szy = size (y);
95    if (n != szy(end))
96      error ("pchip: length of X and last dimension of Y must match");
97    endif
98    y = reshape (y, [prod(szy(1:end-1)), szy(end)]);
99  endif
100
101  h = diff (x);
102  if (all (h < 0))
103    x = fliplr (x);
104    h = diff (x);
105    y = fliplr (y);
106  elseif (any (h <= 0))
107    error ("pchip: X must be strictly monotonic");
108  endif
109
110  f1 = y(:, 1:n-1);
111
112  ## Compute derivatives.
113  d = __pchip_deriv__ (x, y, 2);
114  d1 = d(:, 1:n-1);
115  d2 = d(:, 2:n);
116
117  ## This is taken from SLATEC.
118  h = diag (h);
119
120  delta = diff (y, 1, 2) / h;
121  del1 = (d1 - delta) / h;
122  del2 = (d2 - delta) / h;
123  c3 = del1 + del2;
124  c2 = -c3 - del1;
125  c3 /= h;
126  coeffs = cat (3, c3, c2, d1, f1);
127
128  ret = mkpp (x, coeffs, szy(1:end-1));
129
130  if (nargin == 3)
131    ret = ppval (ret, xi);
132  endif
133
134endfunction
135
136
137%!demo
138%! x = 0:8;
139%! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0];
140%! xi = 0:0.01:8;
141%! yspline = spline (x,y,xi);
142%! ypchip = pchip (x,y,xi);
143%! title ("pchip and spline fit to discontinuous function");
144%! plot (xi,yspline, xi,ypchip,"-", x,y,"+");
145%! legend ("spline", "pchip", "data");
146%! %-------------------------------------------------------------------
147%! % confirm that pchip agreed better to discontinuous data than spline
148
149%!shared x, y, y2, pp, yi1, yi2, yi3
150%! x = 0:8;
151%! y = [1, 1, 1, 1, 0.5, 0, 0, 0, 0];
152%!assert (pchip (x,y,x), y)
153%!assert (pchip (x,y,x'), y')
154%!assert (pchip (x',y',x'), y')
155%!assert (pchip (x',y',x), y)
156%!assert (isempty (pchip (x',y',[])))
157%!assert (isempty (pchip (x,y,[])))
158%!assert (pchip (x,[y;y],x), [pchip(x,y,x);pchip(x,y,x)])
159%!assert (pchip (x,[y;y],x'), [pchip(x,y,x);pchip(x,y,x)])
160%!assert (pchip (x',[y;y],x), [pchip(x,y,x);pchip(x,y,x)])
161%!assert (pchip (x',[y;y],x'), [pchip(x,y,x);pchip(x,y,x)])
162%!test
163%! x = (0:8)*pi/4; y = [sin(x); cos(x)];
164%! y2(:,:,1) = y; y2(:,:,2) = y+1; y2(:,:,3) = y-1;
165%! pp = pchip (x, shiftdim (y2,2));
166%! yi1 = ppval (pp, (1:4)*pi/4);
167%! yi2 = ppval (pp, repmat ((1:4)*pi/4, [5,1]));
168%! yi3 = ppval (pp, [pi/2,pi]);
169%!assert (size (pp.coefs), [48,4])
170%!assert (pp.pieces, 8)
171%!assert (pp.order, 4)
172%!assert (pp.dim, [3,2])
173%!assert (ppval (pp,pi), [0,-1;1,0;-1,-2], 1e-14)
174%!assert (yi3(:,:,2), ppval (pp,pi), 1e-14)
175%!assert (yi3(:,:,1), [1,0;2,1;0,-1], 1e-14)
176%!assert (squeeze (yi1(1,2,:)), [1/sqrt(2); 0; -1/sqrt(2);-1], 1e-14)
177%!assert (size (yi2), [3,2,5,4])
178%!assert (squeeze (yi2(1,2,3,:)), [1/sqrt(2); 0; -1/sqrt(2);-1], 1e-14)
179
180%!error (pchip (1,2))
181%!error (pchip (1,2,3))
182