function sp = spap2(knots,k,x,y,w)
%SPAP2 Least squares spline approximation.
%
%   SPAP2(KNOTS,K,X,Y) returns the spline  f  of order K with knot sequence 
%   KNOTS for which
%
%   (*)    Y = f(X)
%
%   in the mean-square sense, i.e.,  the one that minimizes
%
%   (**)   sum_j W(j) ( Y(j) - f(X(j)) )^2 .
%
%   with  W = ones(size(X)). Other weights can be specified by an optional
%   additional argument, i.e., by using
%
%   SPAP2(KNOTS,K,X,Y,W) which returns the spline  f  of order K with knot
%   sequence KNOTS that minimizes (**). A better choice than the default
%   W = ones(size(X))  would be the composite trapezoid rule weights 
%   W = ([dx;0]+[0;dx]).'/2 , with dx = diff(X(:)) .
%
%   If the data sites satisfy the Schoenberg-Whitney conditions
%
%   (***)   KNOTS(j) < X(j) < KNOTS(j+K) ,
%                               j=1:length(X)=length(KNOTS)-K ,
%
%   (with equality permitted at knots of multiplicity K), then  f  is
%   the unique spline of that order satisfying  (***)  exactly.  No
%   spline is returned unless (***) is satisfied for some subsequence
%   of X.
%   
%   Since it might be difficult to supply such a knot sequence for the given
%   data sites, it is also possible to specify KNOTS as a positive integer,
%   in which case, if possible, a knot sequence will be supplied that satisfies
%   (***) for some subsequence of X and results in a spline consisting of KNOTS 
%   polynomial pieces.
%
%   It is also possible to fit to gridded data:
%
%   SPAP2( {KNOTS1,...,KNOTSm}, [K1,...,Km], {X1,...,Xm}, Y ) returns
%   the m-variate tensor-product spline of coordinate order Ki and with knot 
%   sequence KNOTSi in the i-th variable, i=1,...,m, for which
%
%   Y(:,i1,...,im) = f(X1(i1),...,Xm(im)),  all i := (i1,...,im) 
%
%   in the (possibly weighted) mean-square sense.
%   Note the possibility of fitting to vector-valued data. However,
%   in contrast to the univariate case, if the data to be fitted are
%   scalar-valued, then the input array Y is permitted to be m-dimensional,
%   in which case
%   Y(i1,...,im) = f(X1(i1),...,Xm(im)),  all i := (i1,...,im) 
%   in the (possibly weighted) mean-square sense.
%
%   For example, 
%
%      spap2(augknt(x([1 end]),2),2,x,y);
%
%   provides the least-squares straight-line fit to data x,y, assuming that 
%   x(1) <= x(i) <= x(end), while
%
%      spap2(1,2,x,y);
%
%   will accomplish this without that assumption, and, with that assumption,
%
%      w = ones(size(x)); w([1 end]) = 100;
%      spap2(1,2,x,y,w);
%
%   forces that fit to come very close to the first and last data point.
%
%   For example, the statements
%
%      x = -2:.2:2; y=-1:.25:1; [xx, yy] = ndgrid(x,y); 
%      z = exp(-(xx.^2+yy.^2)); 
%      sp = spap2({augknt([-2:2],3),2},[3 4],{x,y},z);
%      fnplt(sp)
%
%   produce the picture of an approximant to a bivariate function. 
%   Use of MESHGRID instead of NDGRID here would produce an error.
%
%   See also SPAPI, SPAPS.

%   Carl de Boor 24 dec 89
%   cb : 7 jun 94; make use of SPCOL facility to suppress
%                  interpretation of multiplicities as derivative matching
%   cb : 15 sep 94; permit the optional use of weights.
%   cb :  9 may '95 (use .' instead of ')
%   cb : 25nov95 (update help), 12sep96 (update help; string args for SPCOL)
%        26jan97 (missing paren in help)
%   cb : 26oct97 (also handle gridded data)
%   cb : 28nov98 (handle sparse Y)
%   cb : 10may98 (standardize the help)
%   cb : 21jan99 (add bivariate example and mention NDGRID in the help)
%   cb : 06feb00 (trap errors from SLVBLK; allow KNOTS to be just a number) 
%   cb : 20feb00 (improve the help; make sure latest change ok for multivar.) 
%   cb : 10apr00 (introduce optional third argument for multidimensional case;
%                 handle 1 piece of order 1 correctly)
%   cb : 14may00 (reduce number of error returns, handle suppression of
%                                    trailing singleton dimensions correctly)
%   cb : 13jul00 (switch from fprintf(...),error(' '), to error(sprintf(...)) )
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   $Revision: 1.14 $

if nargin<5, w = []; end

if iscell(knots) % gridded data are to be fitted 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
       if length(sizey)<m % interpret missing dimensions as singletons
          sizey = [1 sizey ones(1,m-length(sizey))]; y = reshape(y,sizey);
       else
          error(['If KNOTS is a cell-array of length m, then Y must be an ', ...
               'm- or (m+1)-dimensional array.'])
       end
   end

   if iscell(k), k = cat(2,k{:}); end
   if length(k)==1, k = repmat(k,1,m); end
   if isempty(w), w = cell(1,m); end
   
   v = y; sizev = sizey;
   for i=m:-1:1   % carry out coordinatewise least-squares fitting
      [knots{i},v,sizev(m+1),k(i)] = spbrk(spap21(knots{i}, k(i), x{i}, ...
                      reshape(v,prod(sizev(1:m)),sizev(m+1)),w{i}));
      v = reshape(v,sizev);
      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.
   % It remains to put together the spline:
   sp = spmak(knots, v, sizev);

else             % univariate spline interpolation
   sp = spap21(knots,k,x,y,w);
end

function sp = spap21(knots,k,x,y,w)
%SPAP21 Univariate least squares spline approximation.

% 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(' y  should contain exactly as many entries as  x .')
   else, y = y.'; end
else
   if (cy~=cx), error(' y  should contain exactly as many entries as  x .')
   end
end

%  sort the given data sites.
if any(diff(x)<0), [x,index] = sort(x); y = y(:,index); end, nx = length(x);

if length(knots)==1 % we are to use a spline with  KNOTS  pieces
   k = min(k,nx); maxpieces = nx-k+1;
   if knots<1|knots>maxpieces
      warning(...
      sprintf(['The number of polynomial pieces to be used must be', ...
            ' positive\n',...
            'but not larger than length(x)-k+1 = ',num2str(maxpieces),'.\n']))
      knots = max(1,min(maxpieces,knots));
   end
   if knots==1&k==1
      if nx==1, knots = [x(1) x(1)+1];
      else,     knots = x([1 end]);
      end
   else
      knots = aptknt(x(round(linspace(1,nx,knots-1+k))),k);
   end
%else  % adjust the input to make sure that  0 < length(knots)-k <= length(x)
%   k = max(1, min(length(knots)-1,k));
%   nxpk = nx+k;
%   if length(knots)>nxpk % delete some knots
%     knots = knots(round(linspace(1,end,nxpk)));
%   end
end

%  Generate the collocation matrix and divide it into the possibly reordered
%  sequence of given values to generate the B-spline coefficients of the
%  interpolant, then put it all together into SP. But trap any error from
%  SLVBLK, in order to provide a more helpful error message.

lasterr('')
try
   if ~isempty(w)
      sp = spmak(knots,slvblk(spcol(knots,k,x,'slvblk','noderiv'),y.',w).');
   else
      sp = spmak(knots,slvblk(spcol(knots,k,x,'slvblk','noderiv'),y.').');
   end
catch
  if ~isempty(findstr('slvblk',lasterr))
     error(...
     sprintf(['\nThe given knots and order are incompatible with the given', ...
             ' data sites.\n', ...
       'No subsequence of the data sites satisfies the',...
       ' Schoenberg-Whitney\n',...
       '   conditions of the given order wrto the given knots.\n', ...
       'If you have trouble coming up with a good knot sequence for the', ...
       ' data\n', ...
       '   sites, try using spap2(l,k,x,y), with l the desired number of\n',...
       '   polynomial pieces in the spline fit.']))
  else
     error(lasterr)
  end
end
