% SPAPIDEM	Demonstrate spline interpolation.
% Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
                                          %latest change: November 28, 1990
% latest change: May 16, 1991 change x to tau, ll  to  n , and add material 
% on optimal spline interpolation.

echo on;clc;clg;home
%	Here are some sample data, much used for testing spline approximation
%  with variable knots, the so-called Titanium Heat Data, which record some
%  property of titanium measured as a function of temperature. We'll use it
%  merely to illustrate the use of spline interpolation.

[xx,yy]=titanium;

%	The plot of the data shows a rather sharp peak.

pause				% Touch any key to continue

frame=[-1 1 -.1 .1]+[min(xx),max(xx),min(yy),max(yy)];
plot(xx,yy,'x');
axis(frame);pause
hold on

%     We'll pick a few data points from these somewhat rough data, since we
%  want to interpolate it. 


pick=[1 5 11 21  27 29 31 33 35 40 45 49];
tau=xx(pick);y=yy(pick);

% Here is a picture of the data, with the selected data points marked.
pause				% Touch any key to continue
plot(tau,y,'o');pause
 
%	Since a spline of order  k  with  n+k  knots has  n  degrees of
%  freedom, and we have 12 data points, a fit with a cubic spline, i.e., a
%  fourth order spline, requires  12+4  knots.

%	Moreover, this knot sequence  t  must satisfy the Schoenberg-Whitney
%  conditions, i.e., must be such that the i-th data
%  abscissa lies in the support of the i-th B-spline, i.e.,

%                t(i) < tau(i) < t(i+k) , all  i ,

%  (with equality allowed in case of a knot of multiplicity  k ).
%	We can achieve this condition for our example by using the data
%  abscissa as knots, but adding two knots at either end.

n=length(tau);dl=tau(2)-tau(1);dr=tau(n)-tau(n-1);
t=[tau(1)-dl*[2 1] tau tau(n)+dr*[1 2]];  % construct the knot sequence
sp=spapi(t,tau,y);     % This constructs the spline.

pause				% Touch any key to continue

%	Now, for the plot:

pp=sp2pp(sp);		% converts to pp form

%    We don't care about the part of the curve outside the data interval,
%  so we cut the  pp-rep. down to that interval....

pc=ppcut(pp,[tau(1) tau(n)]);

%           ....  and only then plot it,

values=fnplt(pc,'-.');pause  % saving the plotted points for later use.

% (The statement  values=fnplt(pp,'-.',[tau(1) tau(n)])  would have had the
% same effect, without the necessity of introducing  pc .)
hold off

%	A closer look at the left part of the spline fit shows some
%  undulations.

pause				% Touch any key to continue

xxx=[0:40]/40*(tau(5)-tau(1))+tau(1);
clf;
plot(xxx,fnval(pc,xxx),'-.',tau,y,'o');
axis([tau(1) tau(5) .6 1.2]);
pause;

%     The unreasonable bump in the first interval stems from the fact that
%  our spline goes smoothly to zero at its first knot. Here is a picture
%  of the entire spline.

values2=fnplt(pp,'-.');pause
hold on

%     Here are again the data points as well.

pause 				% Touch any key to continue

plot(tau,y,'o');pause
hold off

%     There are many ways to enforce a more reasonable boundary behavior,
%  e.g., by prescribing the slope (or higher derivatives) at the boundary,
%  but they all come down to making sure that the interval of real interest
%  to us does not extend all the way to the endpoints of the data interval.

%	Here is a simple idea: We add two more data points outside the given
%  data interval and choose as our data there the values of the straight
%  line through the first two data points.

pause				% Touch any key to continue

tt=[tau(1)-[4 3 2 1]*dl tau tau(n)+[1 2 3 4]*dr];
xx=[tau(1)-[2 1]*dl tau tau(n)+[1 2]*dr];
yy=[y(1)-[2 1]*(y(2)-y(1)) y y(n)+[1 2]*(y(n)-y(n-1))];
sp2=spapi(tt,xx,yy);
pp2=sp2pp(sp2);
plot(values(1,:),fnval(pp2,values(1,:)),'-',tau,y,'o');
axis(frame);pause

%     Here is also the original spline fit, dash-dotted, to show the reduction 
%  of the undulation in the first and last interval.

pause				% Touch any key to continue

hold on;plot(values(1,:),values(2,:),'-.');pause;hold off;
%    Finally, here is a closer look at the first four data intervals which
%  shows more clearly the reduction of the undulation near the left end.

pause				% Touch any key to continue

plot(xxx,fnval(pp2,xxx),'-',tau,y,'o',xxx,fnval(pp,xxx),'-.');
axis([tau(1) tau(5) .6 2.2]); pause;
hold off

pause				% Touch any key to finish

%     In  o p t i m a l  spline interpolation,  at points  
%                         tau(1), ..., tau(n)
% say, one chooses the knots so as to minimize the constant in a standard 
% error formula. Specifically, one chooses the first and the last data point 
% as a k-fold knot. The remaining  n-k  knots are supplied by  optknt

help optknt

pause				% Touch any key to continue

% We try this alternative way for choosing knots for interpolation on our
% example, in which we were interpolating by cubic splines to data 
%  (tau(i),y(i)),  i=1,...,n :

 k = 4;   xi = optknt(tau,k);
osp = spapi(augknt([tau(1) xi tau(n)],k),tau,y);

% Here is the plot

pause				% Touch any key to continue

fnplt(osp); pause

%  ... a bit disconcerting!

% And here is also the earlier interpolant, for comparison:

pause				% Touch any key to continue

hold on;
plot(values(1,:), values(2,:),'-.'); pause

%  One can clearly see the given data, as the points were these two
% interpolants coincide, but here are the data as well, to be sure:

pause				% Touch any key to continue

plot(tau,y,'o');pause

% Finally, here are also the (interior) optimal knots:

pause				% Touch any key to continue

plot(xi,0*xi,'x');pause
hold off

%     The knot choice for optimal interpolation is designed to make the 
% maximum over  a l l  functions  f  of the ratio   norm{f - If}/norm{D^k f}  
% of the norm of the interpolation error  f - If  to the norm of the  k-th 
% derivative  D^k f  of the interpoland as small as  possible.  Since our data 
% imply that  D^k f  is rather large, the interpolation error near the flat 
% part of the data is of acceptable size for such an `optimal' scheme.

%  But, to be sure that the optimal knots were calculated correctly, here is
% a test. The (interior) optimal knots  xi  are determined in such a way 
% that the spline function of order  k+1  with interior knots  xi  which
% vanishes at all the data points  tau  is  p e r f e c t , i.e., has an 
% absolutely constant  k-th derivative. 

%   Here is the verification: 
pause				% Touch any key to continue


data = [tau, tau(n)+1]; % add that data point to the right, to pin down a
                        % unique spline of order  k+1  and vanishing at all
                        % the tau(i),
datay=[0*tau,1];   % ... and make the spline nonzero at that point.

perfect=spapi(augknt(sort([xi,data([1 n+1])]),k+1),data,datay);

%  Now plot the supposedly perfect spline, along with its zeros  tau :
% Areas where it is relatively large are places where, all else being equal,
% the error in optimal spline interpolation can be expected to be relatively 
% large. Ideally, one would want to choose the data points  tau  so as to make
% this perfect spline equi-oscillate.

pause				% Touch any key to continue

fnplt(perfect,'-',tau([1 n]));
   title('the perfect spline for the given tau')
hold on
plot(tau,zeros(1,n),'o');pause
hold off
% now plot also the derivatives, including the piecewise constant one
% that is supposed to be absolutely constant, always showing  tau  as well:
dp = perfect;
for j=1:k
   dp=fnder(dp);
   fnplt(dp,'-',tau([1 n]));
   title([int2str(j),'. derivative'])
   hold on
   plot(tau,zeros(1,n),'o');pause
   hold off
end

%   This looks indeed absolutely constant. Note that its breakpoints  xi 
% satisfy the Schoenberg-Whitney conditions with respect to the  tau  and for
%  k = 4 , i.e.
%
%                 tau(i) < xi(i) < tau(i+4),   i = 1,...,n-4 
%
% as is made explicit in the next overlay, in which the intervals  
%  [tau(i)..tau(i+4)]  are plotted, with  xi(i)  marked by an 'x'.
% As you will see, each  xi(i)  lies even very close to the center of the
% corresponding interval [tau(i)..tau(i+4)] :

pause				% Touch any key to continue

hold on
for i=1:(n-k)
   plot([tau(i),tau(i+k)],i*.5e-5*[1 1],'-')
   plot(xi(i),i*.5e-5,'x')
end
hold off
