echo off
%TSPDEM Demo of bivariate tensor product spline work.

%   latest change: May 24, 1991
%   cb :  9 may '95 (use .' instead of ')
%   cb : 31 sep '95 (use updated spbrk, replace use of shg )
%                   (check results against tspap2.m and again tspapi.m)
%                   (take account of the fact that MESH treats arrays differently
%                    from the way bivariate Approximation Theory does)
%   cb : 6 apr 96 (change view settings, to agree with tsp example in doc)
% cb : 15nov97; 22dec97 (correct mistakes pointed out by a user)
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc. 
%   $Revision: 1.10 $

home, close, echo on
% 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, as is now possible in version 2 of this toolbox without any 
% knowledge of the details.
%
% For example, all the calculations to follow below can now be carried out by
% the statements:

%   fnplt( spap2({knotsx,knotsy}, [kx ky], {xx, yy}, z) ), view(150,50)

% See the demonstration M-file EXAMPLES for some simple examples.

pause                                            % Touch any key to continue

% However, it might be of interest to see a detailed description of how 
% approximation to gridded data by tensor products is actually done for 
% bivariate data.  This will also come in handy when some tensor product 
% construct not handled directly by the M-files in this toolbox is needed.
%
% 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]);
[xx,yy] = ndgrid(x,y); z = franke(xx,yy);

%  NOTE the use of NDGRID in preference to MESHGRID.
%  Here is a plot of the data ..

pause                                            % Touch any key to continue

fig1 = figure; mesh(x,y,z.'), view(150,50)
title('data from Franke function'),pause

% NOTE the need to use here  z.'  because of MATLAB's matrix-oriented view
% when plotting an array. This can be a serious hassle in bivariate
% approximation since there it is customary to think of  z(i,j)  as the
% value at  (x(i),y(j))  rather than at  (x(j),y(i)) , as MATLAB's MESH does.
%
% I chose point density somewhat higher near the boundary, to pin down the
% approximation there.

pause                                            % Touch any key to continue

% Next, I choose some 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(x,yy,vals.'), view(150,50)
title('simultaneous approximation to all curves in y-direction'),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                                            % Touch any key to continue

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

coefsy = fnbrk(sp,'c');

% 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(i,r)  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 value  coefsy(i,:)  at  x(i) , all  i , 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' :

coefs = fnbrk(sp2,'c').';

% 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];

% one could 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

fig2 = figure; mesh(xv,yv,values.'), view(150,50)
title('the spline approximant'), 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).';

% In fact, starting with version 2 of this toolbox, fnval and spmak can deal
% directly with multivariate splines, hence this statement can be replaced by

value3 = fnval( spmak({knotsx,knotsy},coefs), {xv,yv} );

% Better yet, the construction of the approximation can be done by  o n e 
% call to  spap2, therefore we can these values directly from the given data
% by the statement

value4 = fnval( spap2({knotsx,knotsy},[kx ky],{x,y},z), {xv,yv} );

pause                                            % Touch any key to continue

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

disp(max(max(abs(values-value2)+abs(values-value3)+abs(values-value4)))/...
     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(x,y,errors.'), view(150,50)
title('error at the given data points'),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(xv,y,valsb.'), view(150,50)
title('simultaneous approximation of all curves in the x-direction'),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:


figure(fig1)
pause                                            % Touch any key to continue
close(fig1)
%    Now comes the second step, to get the actual surface. First, extract
% the coefficients:

coefsx = fnbrk(spb,'c');

% 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':

coefsb = fnbrk(spb2,'c');

%    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

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

% The above steps can now be carried out by just two commands, one for the
% construction of the least-squares approximant, the other for its evaluation
% at a rectangular mesh.

tsp = spap2({knotsx,knotsy},[kx,ky],{x,y},z);
valuet = fnval(tsp,{xv,yv});
% Here, for a check, is the relative difference between the values
% computed earlier and those computed now:
disp( max(max(abs(values-valuet)))/max(max(abs(values))) )

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( [0 1 (y(2:(ly-2))+y(3:(ly-1)))/2 ], ky);

spi = spapi(knotsy,y,z);
coefsy = fnbrk(spi,'c');
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.');

icoefs = fnbrk(spi2,'c').';

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(xv,yv,ivalues.'), view(150,50)
title('the spline interpolant'),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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The above steps can now be carried out by just two commands, one for the
% construction of the interpolant, the other for its evaluation
% at a rectangular mesh.

tsp = spapi({knotsx,knotsy},{x,y},z);
valuet = fnval(tsp,{xv,yv});
% Here, for a check, is the relative difference between the values
% computed earlier and those computed now:
disp( max(max(abs(ivalues-valuet)))/max(max(abs(ivalues))) )

pause                                            % Touch any key to continue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
close(fig2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
