% SPALLDEM	Demonstrate form and usage of splines.
% Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
                                          % latest change: 24 January 1992
clc;clg;home;echo on

%  A  s p l i n e  is specified by its (nondecreasing)  k n o t  sequence

%                      t

% and by its B-spline coefficient sequence

%                      a

% With  n := length(a), and  n+k := length(t), the spline is
% of  o r d e r   k . This means that its polynomial pieces have
% degree < k .

pause;                % Touch any key to continue

%     This means, explicitly, that

%             n
%  spline := sum  B_(i,k)*a(i) ,
%            i=1

% with  B_(i,k)  the i-th B-spline of order  k  for the given knot
% sequence  t  , i.e., the B-spline with knots  t(i), ..., t(i+k) .


%      The B-spline with knots  t(i) <= ... <= t(i+k)  is
% positive on the interval  (t(i) .. t(i+k))  and is zero outside the interval.
% It is pp of order  k  with breaks at the points  t(i), ..., t(i+k) . These
% knots may coincide, and the precise  m u l t i p l i c i t y  of a knot
% governs the smoothness with which the two polynomial pieces join there.

%     Here is the picture of a cubic B-spline, which is typical of B-splines
% with simple knots.

pause;                % Touch any key to continue

fnplt(spmak([0 1.3 1.9 3.1 4],1));pause

% See the demo  bspldem  for various examples, particularly for the interplay
% between knot multiplicity and smoothness.

%	To repeat: A spline is a linear combination of B-splines. Usually, a
%  spline is constructed from some information, like function values and/or
%  derivative values, or as the approximate solution of some ordinary
%  differential equation.  But it is also possible to make up a spline from
%  scratch, by providing its knot sequence and its coefficient sequence to
%  the M-file  spmak .

pause;				%  Touch any key to continue

%	For example, we might supply the uniform knot sequence [1:10] and the
%  coefficient sequence [3:8] . Since there are 10 knots and 6 coefficients,
%  the order must be 4, i.e., we get a cubic spline.

sp=spmak([1:10],[3:8])
echo off;
ttype=11;td=1;tn=6;tt=[1:10];tk=4;tcoefs=[3:8];
zcheck(sp-[ttype td tn tcoefs tk tt]);echo on
pause;				%  Touch any key to continue

%  The output may be a bit bewildering but, as with pp's, you don't really
%  have to know just how the spline is stored in the vector  sp . If you
%  do ever want to know the constituents, you can always unmake the spline 
%  using  spbrk  as follows.
pause					% Touch any key to continue

[knots,coefs,n,k]=spbrk(sp)

pause                                  % Touch any key to continue
echo off;
zcheck(knots-tt);zcheck(coefs-tcoefs);zcheck([n k]-[tn tk]);echo on;

%  But the point of the spline toolbox is that there shouldn't be any need
%  for you to look up these details. You simply use  sp  as an argument to
%  M-files that evaluate, differentiate, integrate, convert, or plot the 
%  spline whose description is contained in  sp .

%      For example, here is a plot of the spline we just made up. Notice
% how this spline vanishes at the endpoints of its basic interval, due to
% the fact that all its knots are simple.

pause                                  % Touch any key to continue
fnplt(sp);pause

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


%     A spline is usually constructed as an approximation to some function.
%  Information about that function is typically available as function values.
%  Here, for example, is a data set which has become a standard test case, the
%  t i t a n i u m  h e a t  data. It is challenging since it is relatively
%  flat except for one sharp peak.
pause                                  % Touch any key to continue
[x,y]=titanium;
plot(x,y,'o');
frame=[x(1),x(length(x)),min(y)-1,max(y)+1]; axis(frame); pause
hold on
%    As a first attempt at fitting this data, we try a C^1 piecewise parabolic
%  with six pieces and uniformly placed breaks:

k=3; l=6; breaks=x(1)+[0:l]*(x(length(x))-x(1))/l;

%    We use the M-file  augknt  to convert the given breakpoint sequence into
%  the knot sequence whose corresponding B-splines provide a basis for the pp
%  space chosen:

knots=augknt(breaks,k)
echo off;
zcheck(knots-[595,595,595,675,755,835,915,995,1075,1075,1075]);echo on

pause                                  % Touch any key to continue

%     The M-file  spap2  provides a least-squares spline fit to given data:

splin0=spap2(knots,k,x,y);

%  We plot this approximation on top of the data:
pause                                  % Touch any key to continue
fnplt(splin0);pause

%   The approximation is not very good, but it is at least correct in the 
% sense that the spline crosses the given function eight times, i.e., at
% least as many times as there are degrees of freedom (i.e., B-splines) in
% our approximating family. Recall that that number is
length(knots)-k
% and here is the plot for another look.
pause                                  % Touch any key to continue
hold off

%    The best way to judge the quality of an approximation is to look at the
% error.  Specifically, one looks for systematic behavior and, when using
% pp functions, tries to reduce the error by placing additional breakpoints
% in the area of largest error. Here is the error plot:
pause                                  % Touch any key to continue
error0=y-fnval(splin0,x);
plot(x,error0);hold on;plot(frame(1:2),[0 0]);pause

%    Now we add breaks at the point or points of greatest error:
break1=sort([breaks,860,900,950]);
splin1=spap2(augknt(break1,k),k,x,y);
erro1=y-fnval(splin1,x);
pause                                  % Touch any key to continue
plot(x,erro1,'-.');pause

%    There is only one systematic part still left. 
% To fit it, we make  900  a double knot, thus allowing the approximating
% spline a discontinuous first derivative (since it is quadratic):
break2=sort([break1,900]);
splin2=spap2(augknt(break2,k),k,x,y);
erro2=y-fnval(splin2,x);
pause                                  % Touch any key to continue

plot(x,erro2,':');pause

%    There is still some systematic oscillation in the error. We place an 
% additional breakpoint at its maximum:
break3=sort([break2,880]);
splin3=spap2(augknt(break3,k),k,x,y);
erro3=y-fnval(splin3,x);
pause                                  % Touch any key to continue
plot(x,erro3,'-');pause

%   This could be improved further, but you get the idea. 

%   Here are the data and the last spline fit constructed:
hold off
plot(x,fnval(splin3,x))
axis(frame);
hold on
plot(x,y,'o');pause

hold off
splst
