function sp = spapi(knots,x,y)
%SPAPI Spline interpolation.
%
%   SPAPI(KNOTS,X,Y)      (with both KNOTS and X vectors)
%   returns the spline  f  (if any) of order
%         k := length(KNOTS) - length(X)
%   with knot sequence KNOTS for which  Y(:,j) = f(X(j)) , all j.
%   This is taken in the osculatory sense in case some X are repeated,
%   i.e., in the sense that  D^m(i) f(X(i)) = Y(:,i)  in case the X are
%   in nondecreasing order, with  m(i) := #{ j<i : X(j) = X(i) }.
%   For this, it is important to realize the following: if the input X is
%   not nondecreasing, it will be reordered to make it so, and the input Y 
%   will be reordered in the same way. See the cubic Hermite interpolation 
%   example below.
%
%   SPAPI(K,X,Y)          (with K a positive integer)
%   uses APTKNT to obtain from K and X an apt knot sequence. Precisely,
%   it has the same effect as  SPAPI(APTKNT(sort(X),K),X,Y).
%
%   SPAPI({KNOTS1,...,KNOTSm},{X1,...,Xm},Y)
%   returns in SP the m-variate tensor-product spline of coordinate order
%      ki := length(KNOTSi) - length(Xi)
%   with knot sequence KNOTSi in the i-th variable, i=1,...,m, for which
%   Y(:,j1,...,jm) = f(X1(j1),...,Xm(jm)),  all j := (j1,...,jm) 
%   As in the univariate case, KNOTSi may also be a positive integer, in which
%   case the i-th knot sequence is obtained from Xi via APTKNT.
%   Note the possibility of interpolating to vector-valued data. However,
%   in contrast to the univariate case, if the data to be interpolated are
%   scalar-valued, then the input array Y is permitted to be m-dimensional,
%   in which case
%   Y(j1,...,jm) = f(X1(j1),...,Xm(jm)),  all j := (j1,...,jm) 
%   Multiplicities in the sequences Xi lead to osculatory interpolation just
%   as in the univariate case.
%
%   For example, if the points in the vector  t  are all distinct, then
%
%      sp = spapi(augknt(t,4,2),[t t],[cos(t) -sin(t)]);
%
%   provides the C^1 piecewise cubic Hermite interpolant to the function 
%   f(x) = cos(x)  at the points in  t . If matching of slopes is only
%   required at some subsequence  s  of  t  but that includes the leftmost
%   and the rightmost point in  t , one would use instead
%
%      sp = spapi( augknt([t s],4), [t s], [cos(t) -sin(s)] );
%
%   or even just
%
%      sp = spapi( 4, [t s], [cos(t) -sin(s)] );
%
%   and the last works even if s fails to include the extreme points of t,
%   and produces even a C^2 piecewise cubic interpolant to these Hermite data.
%
%   As another example,
%
%      sp = spapi( {[0 0 1 1],[0 0 1 1]}, {[0 1],[0 1]}, [0 0;0 1] );
%
%   constructs the bilinear interpolant to values at the corners of a square,
%   as would the statement
%
%      sp = spapi({2,2},{[0 1],[0 1]},[0 0;0 1]);
%
%   As a final example, the statements
%
%      x = -2:.5:2; y=-1:.25:1; [xx, yy] = ndgrid(x,y); 
%      z = exp(-(xx.^2+yy.^2)); 
%      sp = spapi({3,4},{x,y},z);
%      fnplt(sp)
%
%   produce the picture of an interpolant (piecewise quadratic in x, 
%   piecewise cubic in y) to a bivariate function.
%   Use of MESHGRID instead of NDGRID here would produce an error.
%
%   See also SPAPS, SPAP2.

%   Carl de Boor 24 dec 89
%   cb :  22 apr 95
%   cb :  9 may '95 (use .' instead of ')
%   cb :  9 mar '96 (flip  X  in case it has more than one row)
%   mathworks: 19mar96 (clean up error(...))
%   cb : 12sep96 (string argument for SPCOL)
%   cb : 18oct97 (include tensor product)
%   cb : 28nov98 (handle sparse Y, add bivariate example, mention NDGRID)
%   cb : 14feb99 (with KNOTS a positive integer, take it as the order and use
%                 APTKNT to construct a knot sequence apt for the given X)
%   cb : 21nov99, 20feb00 (improve the help, permit k=1)
%   cb : 10apr00 (introduce optional third argument for multidimensional case;
%                 also handle just one data point in case k=1;
%                 also use the actual order used in APTKNT in case the user
%                 specified the order)
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   $Revision: 1.15 $

if iscell(knots) % gridded data are to be interpolated by tensor product splines
   if ~iscell(x)
      error('If KNOTS is a cell-array, then also X must be one.'), end
   m = length(knots);
   if m~=length(x)
      error('If KNOTS is a cell-array, then X must be one of the same length.')
   end
   sizey = size(y);
   switch length(sizey)
     case m  % grid values of a scalar-valued function
        if issparse(y), y = full(y); end 
        sizey = [1 sizey]; y = reshape(y, sizey); 
     case m+1
     otherwise
        error(['If KNOTS is a cell-array of length m, then Y must be an ', ...
               'm- or (m+1)-dimensional array.'])
   end
   
   v = y; sizev = sizey;
   for i=m:-1:1   % carry out coordinatewise interpolation
      temp = spapi1(knots{i},x{i}, reshape(v,prod(sizev(1:m)),sizev(m+1)));
      v = reshape(spbrk(temp,'c'), sizev);
      if length(knots{i})==1, knots{i} = spbrk(temp,'knots'); end
      if m>1
         v = permute(v,[1,m+1,2:m]); sizev(2:m+1) = sizev([m+1,2:m]);
      end
   end
   % At this point, V contains the tensor-product B-spline coefficients;
   % also, the various knot sequences will have been checked. 
   % It remains to return information:
   sp = spmak(knots, v, sizev);

else             % univariate spline interpolation
   sp = spapi1(knots,x,y);
end

function sp = spapi1(knots,x,y)
%SPAPI1 univariate spline interpolation

if length(knots)==1 % the order is being specified; get the knots via APTKNT
   k = knots;
   if k~=fix(k), error('The order must be a natural number.'), end
   if k<1, error('The order must be at least 1.'), end
   if k==1&length(x)==1
      knots = [x, x+1];
   else
      knots = x; if any(diff(knots)<0), knots = sort(knots); end
      [knots,k] = aptknt(knots,k);
   end
else

   if ~isempty(find(diff(knots)<0))
      error('The knot sequence should be nondecreasing.')
   end
   
   n = length(x); npk = length(knots); k = npk-n;
   if k<1
     error(sprintf(['The number of data points, %.0f, should be less than\n',...
             'the number, %.0f, of knots.'],n, npk));
   end
end

% It is assumed that  y  has the same sense as  x  which is important in case
% the entries of  y  are vectors rather than scalars
[rx,cx] = size(x); [ry,cy] = size(y);
if rx>1
   if ry~=rx
     error(sprintf(['Since  x  has more than one row,\n',...
            '   y  should have the same number of rows as  x .']));
   else, x = x.'; y = y.'; end
else,
   if cy~=cx error('Number of data sites and data values should match.')
   end
end

%  sort the given abscissae.
[x,index] = sort(x); y = y(:,index);

%  Generate the collocation matrix and divide it into the possibly reordered
%  sequence of given ordinates to generate the B-spline coefficients of the
%  interpolant, then put it all together into SP.

sp = spmak(knots,slvblk(spcol(knots,k,x,'slvblk'),y.').');
