function [output,p] = csaps(x,y,p,xx,w)
%CSAPS Cubic smoothing spline.
%
%   CSAPS(X,Y,P)  returns the ppform of the cubic smoothing spline for the
%   given data (X,Y)  and depending on the smoothing parameter P from 
%   [0 .. 1].  This smoothing spline  f  minimizes
%
%   P * sum_i W(i)(Y(i) - f(X(i)))^2  +  (1-P) * integral (D^2 f)^2
%
%   with W = ones(length(X),1) the default value for  W .
%   For  P=0, the smoothing spline is the least-squares straight line fit to
%   the data, while, at the other extreme, i.e., for  P=1, it is the `natural' 
%   or variational cubic spline interpolant.  The transition region between 
%   these two extremes is usually only a  rather small range of values for P 
%   and its location strongly depends on the data sites.
%
%    P  will be chosen in an ad hoc fashion somewhere in that small range, in
%   case no P is specified by the input or a negative P is input. In any case,
%   the value of P actually used is optionally returned.
%   If you have difficulty choosing  P but have some feeling for the size
%   of the noise in  Y , consider using instead  spaps(X,Y,tol)  which, in
%   effect, chooses  P  in such a way that
%           sum_i W(i)(Y(i) - f(X(i)))^2  =  tol  
%   thereby obtaining the smoothest function (in the sense that the roughness
%   measure  integral (D^2 f)^2  is smallest) among all functions  f  with
%   sum_i W(i)(Y(i) - f(X(i)))^2  <=  tol  .
%
%   CSAPS(X,Y,P,XX) returns instead the value(s) at XX of the cubic smoothing
%   spline, unless XX is empty, in which case the ppform of the cubic smoothing
%   spline is returned. This latter option is important when the user wants
%   to supply the weights W, as in the following:
%
%   CSAPS(X,Y,P,XX,W) returns, depending on whether or not XX is empty, the
%   ppform, or the values at XX,  of the cubic smoothing spline for the 
%   specified weights W.
%   
%   For example, 
%
%      x = linspace(0,2*pi,21); y = sin(x)+(rand(1,21)-.5)*.1;
%      pp = csaps(x,y,.4,[],[ones(1,11), repmat(5,1,10)]);
%
%   returns a smoothed version of the data that is much closer to the data
%   in the right half, because of the much larger weight there.
%
%   It is also possible to vary the smoothness requirement, by having P be a 
%   sequence rather than a scalar. In that case, the roughness measure is
%   taken to be
%                         integral lambda (D^2 f)^2 ,
%   with  lambda  the piecewise constant function with interior breaks X
%   whose value on  (X(i-1) .. X(i))  is  P(i), all i, while  P(1) continues
%   to be taken as the smoothing parameter, P.
%
%   It is also possible to smooth data on a rectangular grid and
%   obtain smoothed values on a rectangular grid or at scattered
%   points, by the calls
%   VALUES = CSAPS({X1, ...,Xm},Y,P,XX)
%   VALUES = CSAPS({X1, ...,Xm},Y,P,XX,W)
%   or the ppform of the corresponding cubic smoothing spline by the calls
%   PP = CSAPS({X1, ...,Xm},Y,P)
%   PP = CSAPS({X1, ...,Xm},Y,P,[],W)
%   in which Y is expected to have size [d,length(X1),...,.length(Xm)]
%   (or [length(X1),...,.length(Xm)] if the function is to be scalar-valued),
%   and P is either a scalar or an m-vector or a cell of length m with P{i}
%   either a scalar or else a length(Xi)-vector, and XX is either a list
%   of m-vectors XX(:,j) or else a cell-array {XX1, ..., XXm} specifying
%   the m-dimensional grid at which to evaluate the interpolant, and,
%   correspondingly, W, if given, is a cell array of weight sequences for the
%   m dimensions (with an empty W{i} indicating the default choice).
%
%   See also SPAPS, CSAPSDEM.

%   Carl de Boor 2 sep 89
%   cb :  9 may '95 (use .' instead of ')
%   cb : 23 oct '95 (use sparse matrices, handle vector-valued ordinates)
%   cb : 09 mar 96 (correct mistake in sparse matrix formula for R)
%   cb : 03 mar 97 (optionally provide a weight vector)
%   cb : 06oct97 (improve the help)
%   cb : 26oct97 (also handle gridded data)
%   cb : 19may98 (standardize the help)
%   cb : 12sep98 (improve the help)
%   cb : 20sep98 (piecewise constant weight function in roughness measure)
%   cb : 28nov98 (handle sparse Y)
%   cb : 24apr99 (protect against `near-zero' weights)
%   cb : 02aug99 (optionally return the P used; improve the help)
%   cb : 06sep99 (improve the help)
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   $Revision: 1.13 $

if nargin<3, p = -1; end
if nargin<4, xx = []; end
if nargin<5, w = []; end
   
if iscell(x)     % we are to handle gridded data

   m = length(x);
   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 X is a cell-array of length m, then Y must be an ', ...
               'm- or (m+1)-dimensional array.'])
   end
   
   if ~iscell(p)  % because of the possibility of weighted roughness measures
                  % must have P be a cell array in the multivariate case.
      if length(p)~=m, p = repmat(p(1),1,m); end
      p = num2cell(p);
   end
   if isempty(w), w = cell(1,m); end

   v = y; sizev = sizey;
   for i=m:-1:1   % carry out coordinatewise smoothing
      [cs,p{i}] = csaps1(x{i}, reshape(v,prod(sizev(1:m)),sizev(m+1)), ...
                  p{i}, [], w{i});
      [b,v,l,k] = ppbrk(cs);
      breaks{i} = b;
      sizev(m+1) = l*k; 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 pp coefficients;
   % It remains to make up the formal description:
   if isempty(xx)
      output = ppmak(breaks, v);
   else
      output = fnval(ppmak(breaks,v),xx);
   end

else             % we have univariate data

   [output,p] = csaps1(x,y,p,xx,w);

end

function [output,p] = csaps1(x,y,p,xx,w)
%CSAPS1 univariate cubic smoothing spline

n=length(x);[xi,ind]=sort(x);xi=xi(:);
output=[];
if n<2, error('There should be at least two data points.'), end
if all(diff(xi))==0, error('The data sites should be distinct.'), end

[yd,yn] = size(y); % if y happens to be a one-column matrix, change it to
                   % a one-row matrix.
if yn==1, yn=yd; y=reshape(y,1,yn); yd=1; end

if n~=yn
   error('Number of data sites and data values should match.')
end

yi=y(:,ind).'; dd = ones(1,yd);

% remove all points corresponding to relatively small weights (since a
% (near-)zero weight in effect asks for the corresponding datum to be dis-
% regarded while, at the same time, leading to bad condition and even
% division by zero).
if isempty(w)
   w = ones(n,1);
else
   w = w(ind); maxw = max(abs(w));
   ignorep = find( w <= (1e-13)*maxw );
   if ~isempty(ignorep)
      w(ignorep) = []; xi(ignorep) = []; yi(ignorep,:) = [];
      n = length(xi);
      if n<2
        error('There should be at least two data points with positive weights.')
      end
   end
end

dx=diff(xi); divdif=diff(yi)./dx(:,dd);
if n==2 % the smoothing spline is the straight line interpolant
   pp=ppmak(xi.',[divdif.' yi(1,:).'],yd);
else % set up the linear system for solving for the 2nd derivatives at  xi .
     % this is taken from (XIV.6)ff of the `Practical Guide to Splines'
     % with the diagonal matrix  D^2 there equal to diag(1/w) here.
     % Make use of sparsity of the system.

   dxol = dx;
   if length(p)>1    % check for weights in roughness measure
      if length(p)~=n
         error('P should either be a scalar, or a length(X)-vector.')
      end
      lam = p(2:end).'; p = p(1);
      dxol = dx./lam;
   end
      
   R = spdiags([dxol(2:n-1), 2*(dxol(2:n-1)+dxol(1:n-2)), dxol(1:n-2)],...
                                         -1:1, n-2,n-2);
   odx=ones(n-1,1)./dx;
   Qt = spdiags([odx(1:n-2), -(odx(2:n-1)+odx(1:n-2)), odx(2:n-1)], ...
                                                0:2, n-2,n);
   % solve for the 2nd derivatives
   if isempty(w), w = ones(n,1); end
   W = spdiags(ones(n,1)./w(:),0,n,n);
   if p<0 % we are to determine an appropriate P
      QtWQ = Qt*W*Qt.'; p = 1/(1+trace(R)/(6*trace(QtWQ)));
      u=(6*(1-p)*QtWQ+p*R)\diff(divdif);
   else
      u=(6*(1-p)*Qt*W*Qt.'+p*R)\diff(divdif);
   end
   % ... and convert to pp form
   % Qt.'*u=diff([0;diff([0;u;0])./dx;0])
   yi = yi - ...
    (6*(1-p))*W*diff([zeros(1,yd)
                 diff([zeros(1,yd);u;zeros(1,yd)])./dx(:,dd)
                 zeros(1,yd)]);
   c3 = [zeros(1,yd);p*u;zeros(1,yd)];
   c2=diff(yi)./dx(:,dd)-dxol(:,dd).*(2*c3(1:n-1,:)+c3(2:n,:));
   if exist('lam')
      dxtl = dx.*lam;
      pp=ppmak(xi.',...
      reshape([(diff(c3)./dxtl(:,dd)).',3*(c3(1:n-1,:)./lam(:,dd)).', ...
                                    c2.',yi(1:n-1,:).'], (n-1)*yd,4),yd);
   else
      pp=ppmak(xi.',...
      reshape([(diff(c3)./dx(:,dd)).',3*c3(1:n-1,:).',c2.',yi(1:n-1,:).'],...
                                                            (n-1)*yd,4),yd);
   end
end

if isempty(xx)
   output=pp;
else
   output=ppual(pp,xx);
end
