%TSPDEM	Demonstrate tensor product splines.
% Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
                                    % latest change: May 24, 1991
clc;clg;echo on;home
% Since the toolbox can handle splines with  v e c t o r  coefficients, it is
% easy to implement interpolation or approximation to gridded data by tensor
% product splines.
% 
% Consider, for example, least-squares approximation to given data
%  z(i,j) = f(x(i),y(j)), i=1,...,I, j=1,...,J , e.g.,

x = sort([[0:10]/10,.03 .07, .93 .97]);
y = sort([[0:6]/6,.03 .07, .93 .97]);
z = franke(x'*ones(1,length(y)),ones(length(x),1)*y);

%  Here is a plot of the data ..

pause                                               %Touch any key to continue

mesh(z,[50,50]);pause

% It doesn't give an accurate image of the data because the MATLAB function
%  mesh  has no facilities for data from a  n o n uniform rectangular grid.
% I chose point density somewhat higher near the boundary, to pin down the
% approximation there.

pause                                               %Touch any key to continue

% Next, we choose some appropriate knot sequence and spline order for the 
% y-direction, e.g.,

ky = 3; knotsy = augknt([0,.25,.5,.75,1],ky);

% Then

sp=spap2(knotsy,ky,y,z);

% provides a spline approximation to all the curves  (y,z(i,:)) .
% In particular, the statement

yy=[-.1:.05:1.1];
vals = fnval(sp,yy);

% provides the array  vals , whose entry  vals(i,j)  can be taken as an
% approximation to the value  f(x(i),yy(j))  of the underlying function
%  f  at the mesh-point (x(i),yy(j)) .  This is evident when we plot  vals  
% using  mesh :

pause                                               %Touch any key to continue

mesh(vals,[50,50]);pause

% Note that, for each  x(i) , both the first two and the last two values
% are zero since both the first two and the last two points in  yy  are
% outside the basic interval for the spline in  sp .

% Note also the ridges. They confirm that we are plotting here smooth curves in
% one direction only.

pause

% To get an actual surface, you now have to go a step further.  Look at the
% coefficients  coefsy  of the spline in  sp :

[ignored,coefsy]=spbrk(sp);

% Abstractly, you can think of the spline in  sp  as the function

%       y |--> sum_r coefsy(r,:) B_{r,ky}(y) 

% with the i-th entry  coefsy(r,i)  of the vector coefficient  coefsy(r,:)
% corresponding to  x(i) , all  i .  This suggests approximating each 
% coefficient vector  coefsy(r,:)  by a spline of the same order  kx  and with 
% the same appropriate knot sequence  knotsx :

kx = 4; knotsx=augknt([0:.2:1],kx);
sp2=spap2(knotsx,kx,x,coefsy');

pause                                               %Touch any key to continue

% Note that  spap2(knots,k,x,fx)  expects  fx(:,j)  to be the datum at  x(j) ,
% i.e., expects each  c o l u m n  of  fx  to be a datum. Since we wanted to 
% fit the data  coefsy(r,:)  at  x(r) , all  r , we have to present  spap2  
% with the  t r a n s p o s e  of  coefsy .

% Now consider the transpose   coefs = cxy'  of the coefficients  cxy  of the 
% resulting spline `curve' :

[ignored,cxy]=spbrk(sp2);
coefs=cxy';

% It provides the  b i v a r i a t e  spline approximation
% 
%     (x,y) |--> sum_q sum_r coefs(q,r) B_{q,kx}(x) B_{r,ky}(y)
% 
% to the original data
% 
%    (x(i),y(j)) |-->  z(x(i),y(j)) .

pause                                               %Touch any key to continue

% To evaluate this spline surface at the grid (xv(i),yv(j)), e.g., the grid

xv=[0:.025:1];yv=[0:.025:1];

% you can do the following:

values = spcol(knotsx,kx,xv)*coefs*spcol(knotsy,ky,yv)';

% Here is the mesh-plot of these values: 

pause                                               %Touch any key to continue

mesh(values,[50,50]);pause

% This makes good sense since  spcol(knotsx,kx,xv)  is the matrix whose
% (i,q)-entry equals the value  B_{q,kx}(xv(i))  at  xv(i)  of the  q-th 
% B-spline of order  kx  for the knot sequence  knotsx .
 
pause                                               %Touch any key to continue

% Since the matrices  spcol(knotsx,kx,xv)  and  spcol(knotsy,ky,yv)  are 
% banded, it may be more efficient (though perhaps more memory-consuming)
% for `large'  xv  and  yv  to make use of   fnval , as follows:

value2=fnval(spmak(knotsx,fnval(spmak(knotsy,coefs),yv)'),xv)';
 

% Here is a check, viz. the  r e l a t i v e  difference between the values 
% computed in these two different ways:

disp( max(max(abs(values-value2)))/max(max(abs(values))) )

pause                                               %Touch any key to continue

% Here is also a plot of the error, i.e., the difference between the given data
% and the value of the approximation at those data points:

errors = z -  spcol(knotsx,kx,x)*coefs*spcol(knotsy,ky,y)';

pause                                               %Touch any key to continue

mesh(errors,[50,50]);pause

% The  r e l a t i v e  error is of size
disp( max(max(abs(errors)))/max(max(abs(z))) )

% This is perhaps not too impressive. On the other hand, we used only ...
disp(size(coefs))
% ... coefficients to fit a data array of size
disp(size(z))


pause                                               %Touch any key to continue

%    The approach followed here seems  b i a s e d :  We first think of the
% given data  z  as describing a vector-valued function of  y , and then we
% treat the matrix formed by the vector coefficients of the approximating 
% curve as describing a vector-valued function of  x .

%    What happens when we take things in the opposite order, i.e., think of
%  z  as describing a vector-valued function of  x , and then treat the matrix
% made up from the vector coefficients of the approximating curve as describing
% a vector-valued function of  y ?

%    Perhaps surprisingly, the final approximation is the same (up to
% roundoff). Here is the numerical experiment.

%    First, we fit a spline curve to the data, but this time with  x  as the
% independent variable, hence it is the   r o w s  of  z  which now become the
% data points. Correspondingly, we must supply  z' (rather than  z ) to spap2:

pause                                               %Touch any key to continue

spb=spap2(knotsx,kx,x,z');

% ... thus obtaining a spline approximation to all the curves  (x,z(:,j)) .
% In particular, the statement

valsb = (fnval(spb,xv))';

% provides the array  valsb , whose entry  valsb(i,j)  can be taken as an
% approximation to the value  f(xv(i),y(j))  of the underlying function
%  f  at the mesh-point (xv(i),y(j)) .  This is evident when we plot  valsb  
% using  mesh :

pause                                               %Touch any key to continue

mesh(valsb,[50,50]);pause

% Note the ridges. They confirm that we are plotting here smooth curves in
% one direction only.

pause                                               %Touch any key to continue

% Here, for a quick comparison is the earlier mesh of curves, with the smooth
% curves running in the other direction:

pause                                               %Touch any key to continue

mesh(vals,[50,50]);pause

%    Now comes the second step, to get the actual surface. First, extract 
% the coefficients:

[ignored,coefsx]=spbrk(spb);

% Abstractly, you can think of the spline in  spb  as the function

%       x |--> sum_r coefsx(r,:) B_{r,kx}(x) 

% with the j-th entry  coefsy(r,j)  of the vector coefficient  coefsx(r,:)
% corresponding to  y(j) , all  j .  Thus we now fit each
% coefficient vector  coefsx(r,:)  by a spline of the same order  ky  and with 
% the same (appropriate) knot sequence  knotsy :

spb2=spap2(knotsy,ky,y,coefsx');

pause                                               %Touch any key to continue

% Note that we need again to transpose the coefficient array from  spb , 
% since  spap2  takes the columns of that last input argument as the data 
% points.  

% Correspondingly, there is no need now to transpose the coefficient array
%  coefsb  of the resulting  `curve':

[ignored,coefsb]=spbrk(spb2);

%    The claim is that  coefsb  equals the earlier coefficient array  coefs
% (up to round-off), and here is the test:

disp( max(max(abs(coefs - coefsb))) )

%     Not bad, eh? - See the User's Guide for an explanation.

%  Thus, the  b i v a r i a t e  spline approximation
% 
%     (x,y) |--> sum_q sum_r coefsb(q,r) B_{q,kx}(x) B_{r,ky}(y)
% 
% to the original data
% 
%    (x(i),y(j)) |-->  z(x(i),y(j)) 
%
% obtained coincides with the earlier one (which generated  coefs  rather than
%  coefsb ).


pause                                               %Touch any key to continue

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%    Since the data come from a smooth function, we should be interpolating
% it, i.e., use  spapi  instead of  spap2 , or, equivalently, use  spap2
% with the appropriate knot sequences. For illustration, here is the same 
% process done with  spapi .

pause                                               %Touch any key to continue

% To recall, the data points were
% x = sort([[0:10]/10,.03 .07, .93 .97]);
% y = sort([[0:6]/6,.03 .07, .93 .97]);
% 

% We use again quadratic splines in  y , hence use knots midway between data
% points:

ly = length(y);
knotsy = augknt(sort([0 1 (y(2:(ly-2))+y(3:(ly-1)))/2]),ky);

spi=spapi(knotsy,y,z);
[ignored,coefsy] = spbrk(spi);
pause                                               %Touch any key to continue

% We use again cubics in  x , and use the not-a-knot condition, therefore use
% all but the second and the second-last data point as knots:

lx = length(x);
knotsx = augknt(x([1,[3:(lx-2)],lx]),kx);
spi2 = spapi(knotsx,x,coefsy');

[ignored,cxy] = spbrk(spi2);
icoefs = cxy';

pause                                               %Touch any key to continue

% We are ready to evaluate the interpolant at a fine mesh:

ivalues = spcol(knotsx,kx,xv)*icoefs*spcol(knotsy,ky,yv)';

% Here is the mesh-plot of these values: 

pause                                               %Touch any key to continue

mesh(ivalues,[50,50]);pause


% Its error, as an approximation to the Franke function, is computed next:

fvalues = franke(xv'*ones(1,length(yv)),ones(length(xv),1)*yv);
error =  fvalues - ivalues;

% The error is of  r e l a t i v e  size ... 
disp( max(max(abs(error)))/max(max(abs(fvalues))) )

%  ... which is better than the error at the data points of the
% earlier least-squares approximation. For comparison, here is also the 
% overall relative error of that least-squares approximation:

disp( max(max(abs(fvalues-values)))/max(max(abs(fvalues))) )
pause                                               %Touch any key to continue

