echo off
%PCKKNTDM Knot choices demo.

%   cb: 29nov97 (part of former SPAPIDEM elaborated) 20may98
%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   $Revision: 1.10 $

home, close, echo on
%   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
%   to illustrate some methods for knot selection.

[xx,yy]=titanium;

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

frame=[-1 1 -.1 .1]+[min(xx),max(xx),min(yy),max(yy)];axis(frame);
plot(xx,yy,'x'), grid off
title('the titanium heat data show a pronounced peak')
hold on
pause            % Touch any key to continue


%     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.

plot(tau,y,'ro')
title('titanium heat data with selected points marked')

pause            % Touch any key to continue

%   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
%  site 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 ).

%  One way to choose a knot sequence satisfying all these conditions is
%  as the optimal knots, of Gaffney/Powell and Micchelli/Rivlin/Winograd.

pause                 % Touch any key to finish

%     In  o p t i m a l  spline interpolation, to values 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 .

% Here is the beginning of the help from OPTKNT:

pause            % Touch any key to continue

%OPTKNT Optimal knot distribution.
%
%   OPTKNT(TAU,K)  returns the optimal knot sequence (according to the 
%   optimal recovery theory of Micchelli/Rivlin/Winograd and Gaffney/Powell)
%   for interpolation at data points TAU(1), ..., TAU(n) by splines
%   of order K.
%   OPTKNT(TAU,K,MAXITER) permits specification of the number MAXITER of
%   iterations to be tried, the default being 10.
%
%   The interior knots are the  n-K  sign-changes in the absolutely constant
%   function  h  that satisfies
%
%          integral{ f(x)h(x) : TAU(1) < x < TAU(n) } = 0
%
%   for all splines  f  of order K with knot sequence TAU.
%
%   For example, Micchelli/Rivlin/Winograd and Gaffney/Powell would 
%   approve of the following way of constructing a spline interpolant
%   of order k to given data  x, y :
%
%      sp = spapi(optknt(x,k),x,y);

pause            % Touch any key to continue

% We try this for interpolation on our example, in which we want to
% interpolate by cubic splines to data  (tau(i),y(i)),  i=1,...,n :

k = 4; 
osp = spapi( optknt(tau,k), tau,y);

% Here is the plot (after you ...

pause            % ... touch any key to continue)

fnplt(osp,'r'), hold on, plot(tau,y,'ro'), grid off
title('spline interpolant to the selected data using optimal knots')
xlabel('...  a bit disconcerting')

%  ... a bit disconcerting!

pause            % Touch any key to continue

%  Here are also the (interior) optimal knots:
xi = fnbrk(osp,'knots'); xi([1:k end+1-(1:k)]) = [];
plot(xi,zeros(size(xi)),'x')
xlabel('...  and the optimal knots (x)')

pause            %Touch any key to continue

%     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.

% Actually, for these data, the ordinary cubic spline interpolant provided by
% CSAPI does quite well:

cs = csapi(tau,y); fnplt(cs,'b',1.2)
xlabel('...  and the optimal knots (x), and the interpolant provided by CSAPI.')
hold off

pause            %Touch any key to continue
clf

%     Knots must be selected when doing least-squares approximation by splines.
% One approach is to use equally spaced knots to begin with, then use NEWKNT
% with the approximation obtained for a better knot distribution.

% This is now illustrated with the full titanium heat data set.

pause            %Touch any key to continue

sp = spap2(augknt(linspace(xx(1), xx(end), 2+fix(length(xx)/4)),k),k, xx,yy);

fnplt(sp,'r'), hold on, plot(xx,yy,'x'), axis(frame)

title('Least-squares fit by cubic spline (red) with uniform knot sequence,')

pause            %Touch any key to continue

% This is not at all satisfactory. So we use newknt for a spline approximation
% of the same order and with the same number of polynomial pieces, but better
% distributed:

spgood = spap2( newknt(sp), k, xx,yy);
fnplt(spgood,'g',1.5)
xlabel('... and (green) with a knot sequence chosen from it by NEWKNT.')
hold off, pause

% This is quite good. Incidentally, one less interior knot would not have
% sufficed in this case.

pause            %Touch any key to finish
echo off, close

