% Do polynomials fit to data points ( y(i), b(i) ), i=1:n, 
% using Least Squares, via truncated SVD
% (treat singular values less than tol*norm(A) as zero)
% 
% Set up data
y = (-5:.5:6)';
y = (-5:.25:6)';
b = sin(y/5*pi)+y/5;
% degrees contains polynomials degrees to be tried
degreeserr = [1:20];
degreeserr = [1:40];
degreesplot = [1,3,6,19,-1];

% Initialize internal and plotting variables
hold off, clg
n=length(y);
clear xsav, clear yeval, clear beval, clear err, clear condA
ymin=2*min(y);
ymax=2*max(y);
bmin=2*min(b);
bmax=2*max(b);
hold off, clg
switch = [-1,+1,-1,1];
labeltype = 0;
labels = ['-w -.w--w:w '];
plotptr=1;
subplot(121)

% try for polynmials of degree 1 to d
for di=degreeserr;
%  Build A matrix
   A=ones(n,1);
   for i=1:di, A=[A,A(:,i).*y]; end
   scl=sum(abs(A));
   A = A * diag(ones(size(scl))./scl);
%  Solve least squares problem to get polynomial coefficients
   [U,S,V] = svd(A,0);
   iS = zeros(di+1,di+1);
   for i=1:di+1,
      if ( S(i,i) >= tol*S(1,1) ) 
         iS(i,i) = 1/S(i,i);
      else
         iS(i,i) = 0;
      end
   end
   x = V*(iS*(U'*b));
%  x=A\b;
%  Save condition number
   condA = [condA; cond(A)];
%  save solution
   xsav = [xsav,[x;zeros(max(degreeserr)+1-di,1)]];
%  Compute error of fitted polynomial at data points
   err = [err;[di,norm(A*x-b,2)]];
%  Plot original data
   plot(y,b,'wo');
   axis([ymin,ymax,bmin,bmax])
   title ('Original data (o), and polynomials fits')
   xlabel(['Trunc. SVD tol = ',num2str(tol)])
   hold on
%  Evaluate and plot polynomial at k points
   if (degreesplot(plotptr) == di) 
      k = 100;
      yeval=(ymin:(ymax-ymin)/k:ymax);
      xud=flipud(diag(scl)\x);
      beval=polyval(xud,yeval);
      label = labels(labeltype*3+1:labeltype*3+3);
      hold on
      plot(yeval,beval,label)
%     Label graph
      if ( (beval(k+1) <= bmax ) & (beval(k+1) >= bmin) )
         text(yeval(k+1),beval(k+1),['deg = ',int2str(di)])
      elseif ( switch(plotptr) == 1 )
         kk = max(find( ( beval <= bmax ) & ( beval >= bmin ) ) );
         text(yeval(kk),beval(kk),['deg = ', int2str(di)])
      else 
         kk = min(find( ( beval <= bmax ) & ( beval >= bmin ) ) );
         text(yeval(kk),beval(kk),['deg = ', int2str(di)])
      end
%     disp('hit return to continue'), pause
      plotptr = plotptr+1;
      labeltype = rem(labeltype+1,4);
%     if (labeltype == 0) hold off, clg, end
   end
%
end
subplot(122)
semilogy(err(:,1),err(:,2),'w+'), grid
hold on
semilogy(err(:,1),eps*condA,'wo'), 
axis([0,max(degreeserr)+1,1e-16,10])
title('Least squares error')
xlabel('Degree of fitted polynomial')
% print -deps PolynomialFit.eps
% hold on
% title('Least squares error = +, reciprocal condition number = o')
% semilogy(err(:,1),ones(size(condA))./condA,'wo')
