function main_leastsquares

%DRIVER_line
%DRIVER_plane
%DRIVER_plane_svd_wrong
%DRIVER_plane_svd
DRIVER_leastsq
DRIVER_interp

return

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

function DRIVER_line

rng(0)

xx = linspace(-1,1,100)';
c1 = 1;
c2 = 2;
yy = c1 + c2*xx + 0.5*randn(size(xx));

plot(xx,yy,'r.',...
     'MarkerSize',20,'LineWidth',3)
legend('given function samples')
keyboard

A      = zeros(length(xx),2);
A(:,1) = ones(size(xx));
A(:,2) = xx;
cc     = A\yy;   % Matlab notation for solving the least squares problem
fprintf(1,'Exact     values:   c1 = %8.5f   c2 = %8.5f\n',c1,c2)
fprintf(1,'Estimates values:   c1 = %8.5f   c2 = %8.5f\n',cc(1),cc(2))

plot(xx,yy,'r.',...
     xx,cc(1) + cc(2)*xx,'k-',...
     'MarkerSize',20,'LineWidth',3)
legend('given function samples','best fit line')
keyboard

return


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

function DRIVER_plane

nsamp = 500;
%[TT1,TT2] = meshgrid(linspace(-1,1,21));
%tt1 = reshape(TT1,numel(TT1),1);
%tt2 = reshape(TT2,numel(TT2),1);
tt1 = 1 - 2*rand(nsamp,1);
tt2 = 1 - 2*rand(nsamp,1);
c0  = 2;
c1  = 2;
c2  = -1;
yy  = c0 + c1*tt1 + c2*tt2 + 0.5*randn(size(tt1));

hold off
plot3(tt1,tt2,yy,'r.','MarkerSize',20)
legend('given function samples')
keyboard

A      = zeros(length(tt1),3);
A(:,1) = ones(size(tt1));
A(:,2) = tt1;
A(:,3) = tt2;
cc     = A\yy;
fprintf(1,'Exact     values:   c0 = %8.5f   c1 = %8.5f   c2 = %8.5f\n',c0,c1,c2)
fprintf(1,'Estimates values:   c0 = %8.5f   c1 = %8.5f   c2 = %8.5f\n',cc(1),cc(2),cc(3))

[XX1,XX2] = meshgrid(linspace(-1,1,11));

hold off
plot3(tt1,tt2,yy,'r.','MarkerSize',20)
hold on
mesh(XX1,XX2,cc(1) + cc(2)*XX1 + cc(3)*XX2,...
     'LineWidth',2)
hold off
legend('given function samples','best fit plane')
keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function uses the SVD to find a "best fit" plane to a cloud of
% points in three dimensions. It contains an intentional bug in that it
% does not shift the cloud of points to "center" it around the origin.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_plane_svd_wrong

nsamp = 500;
tt1   = 1 - 2*rand(nsamp,1);
tt2   = 1 - 2*rand(nsamp,1);
yy    = 2 + 2*tt1 - tt2 + 0.5*randn(size(tt1));

hold off
plot3(tt1,tt2,yy,'r.','MarkerSize',20)
legend('given function samples')
keyboard

A       = [tt1,tt2,yy];
[U,D,V] = svd(A);
Aproj   = A*V(:,1:2)*V(:,1:2)';

plot3(tt1,tt2,yy,'r.',...
      Aproj(:,1),Aproj(:,2),Aproj(:,3),'b.',...
      'MarkerSize',20)
legend('given function samples','best fit plane')
keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function uses the SVD to find a "best fit" plane to a cloud of
% points in three dimensions.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_plane_svd

nsamp = 500;
tt1   = 1 - 2*rand(nsamp,1);
tt2   = 1 - 2*rand(nsamp,1);
c1    = 2;
c2    = 2;
c3    = -1;
yy    = c1 + c2*tt1 + c3*tt2 + 0.2*randn(size(tt1));

hold off
plot3(tt1,tt2,yy,'r.','MarkerSize',20)
legend('given function samples')
keyboard

A       = [tt1,tt2,yy];
xxmean  = (1/nsamp)*sum(A,1);
Ashift  = A - ones(nsamp,1)*xxmean;
[U,D,V] = svd(Ashift);
Aproj   = ones(nsamp,1)*xxmean + Ashift*V(:,1:2)*V(:,1:2)';
fprintf(1,'Exact     values:   c1 = %8.5f   c2 = %8.5f   c3 = %8.5f\n',c1,c2,c3)
fprintf(1,'Estimates values:   c1 = %8.5f   c2 = %8.5f   c3 = %8.5f\n',...
        xxmean(3) - V(3,1:2)*(V(1:2,1:2)\xxmean(1:2)'),-V(1,3)/V(3,3),-V(2,3)/V(3,3))

plot3(tt1,tt2,yy,'r.',...
      Aproj(:,1),Aproj(:,2),Aproj(:,3),'b.',...
      'MarkerSize',20)
legend('given function samples','best fit plane')
keyboard

return

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

function DRIVER_leastsq

nsamp = 200;
tt    = linspace(-1,1,nsamp)';
%fun = @(t) exp(-.5*t).*cos(5*t);
fun = @(t) 1./(1 + 10*t.*t);

plot(tt,fun(tt),'r',...
     'MarkerSize',20,'LineWidth',2)
legend('given function')
keyboard

n   = 5;
xx  = linspace(-1,1,n)';
%xx  = cos(linspace(0,pi,n)');
A   = zeros(length(tt),length(xx));
for i = 0:(n-1)
  A(:,i+1) = tt.^(n-1-i);
end
cc  = A\fun(tt);

plot(tt,fun(tt),'r',...
     tt,polyval(cc,tt),'b-',...
     'MarkerSize',20,'LineWidth',2)
legend('given function',sprintf('degree %d best fit poly',n))
keyboard

return

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

function DRIVER_interp

tt  = linspace(-1,1)';
%fun = @(t) exp(-.5*t).*cos(5*t);
%fun = @(t) 1./(1 + t.*t);
fun = @(t) 1./(1 + 10*t.*t);

plot(tt,fun(tt),'r',...
     'MarkerSize',20,'LineWidth',2)
legend('given function')
keyboard

n   = 10;
xx  = linspace(-1,1,n)';
%xx  = cos(linspace(0,pi,n)');
A   = vander(xx);
cc  = A\fun(xx);

plot(tt,fun(tt),'r',...
     tt,polyval(cc,tt),'b-',...
     xx,fun(xx),'b.',...
     'MarkerSize',20,'LineWidth',2)
legend('given function',sprintf('degree %d interpolant',n))
keyboard

return

