%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Solution to Problems 1 and 3 of
%  Blind Deconvolution: Errors, Errors Everywhere     
%  Dianne P. O'Leary 
%  Computing in Science and Engineering
%
%  In this problem, we use least squares (LS) and total
%  least squares (TLS) to deconvolve a spectrum.
%  The model is (K+E)f = g+r, where K and g are given
%  in the file spectdata.mat.
%  For LS we force E to be zero.
%
%  problem1_and_3.m   11/2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all
load spectdata
count = g;
[nrows,ncols]=size(K);
startcut = ncols - 10;
keep = [12 15 17 21]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the data in Figure 1.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)
plot(bincenters,count)
title('Observed counts')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem 1: Truncated Least Squares 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('Least squares')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Do the SVD once.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[u,s,v]=svd(K);
strial = diag(s);
ghat = u'*count;
estvariance = sum(ghat(ncols+1:nrows).^2)/(nrows-ncols)
expectedresidnorm = sqrt(estvariance*nrows)
condition_number = strial(1)/strial(ncols)

figure(3)
subp = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute LS solutions for various values of the cutoff
% parameter, the number of singular values retained.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for cutoff = startcut:ncols
   y = zeros(ncols,1);
   y(1:cutoff) = ghat(1:cutoff)./strial(1:cutoff);
   trialcounts = v*y;
   residnorm(cutoff) = norm(K*trialcounts-count);
   solnnorm(cutoff) = norm(trialcounts);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pad with zeros:
% assume 1st two and last 3 true counts are zero.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   trialcounts = [0;0;trialcounts;0;0;0];
   plot(bincenters,count,'g.',bincenters,trialcounts,'b')
   legend('Observed counts', 'Reconstructed counts')
   disp(sprintf( ...
    'cutoff=%d, solnnorm=%f, residnorm=%f,expectedresidnorm=%f' ...
    ,cutoff,solnnorm(cutoff), residnorm(cutoff), expectedresidnorm))
   disp(sprintf('Pausing for display.  Press any key to continue.'))
   pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot sample solutions for Least Squares in Figure 2.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   if (sum(cutoff==keep)==1)
     figure(2)
     subp = subp + 1;
     subplot(2,2,subp)
     plot(bincenters,trialcounts,'b')
     title(sprintf('Retaining %d singular values',cutoff))
     figure(3)
     trialcountssav = trialcounts;
   end
end

figure(2)
print -depsc ls.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the L-curve for Least Squares in Figure 3.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(3)
loglog(residnorm(startcut:ncols),solnnorm(startcut:ncols))
hold on
loglog(residnorm(startcut:ncols),solnnorm(startcut:ncols),'o')
axis([min(residnorm(startcut:ncols)) max(residnorm(startcut:ncols)) ...
      min(solnnorm(startcut:ncols)) max(solnnorm(startcut:ncols))])
xlabel('residual norm')
ylabel('solution norm')
title('L-curve for LS')
hold off
print -depsc lcurvels.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem 3: Truncated Total Least Squares
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(' ')
disp('Total least squares solutions')

subp = 0;
figure(5)
[ut,st,vt] = svd([K,g]);
stdiag = diag(st);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute TLS solutions for various values of the cutoff
% parameter, the number of singular values retained.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for cutoff = startcut:ncols, 
   vnorm = norm(vt(ncols+1,cutoff+1:ncols+1));
   f = -vt(1:ncols,cutoff+1:ncols+1) ...
            *vt(ncols+1,cutoff+1:ncols+1)'/vnorm^2;

   f = [0;0;f;0;0;0];
   if (cutoff == 15)
     ffinal = f;
   end
   hold off
   plot(bincenters,count,'g.', bincenters,f,'b')
   disp(sprintf( ...
    'cutoff=%d, solnnorm=%f, residnorm=%f' ...
    ,cutoff,solnnorm(cutoff), residnorm(cutoff)))
   disp(sprintf('Pausing for display.  Press any key to continue.'))
   pause
   residnorm(cutoff) = norm(stdiag(cutoff+1:ncols+1));
   solnnorm(cutoff) = norm(f);
   if (sum(cutoff==keep)==1)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot sample solutions for Total Least Squares in Figure 4.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     figure(4)
     subp = subp + 1;
     subplot(2,2,subp)
     plot( bincenters,f,'b')
     title(sprintf('Retaining %d singular values',cutoff))
     figure(5)
     fsav = f;
   end
end

figure(4)
print -depsc tls.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the L-curve for Total Least Squares in Figure 5.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(5)
loglog(residnorm(startcut:ncols),solnnorm(startcut:ncols))
hold on
loglog(residnorm(startcut:ncols),solnnorm(startcut:ncols),'o')
axis([min(residnorm(startcut:ncols)) max(residnorm(startcut:ncols)) ...
      min(solnnorm(startcut:ncols)) max(solnnorm(startcut:ncols))])
xlabel('residual norm')
ylabel('solution norm')
title('L-curve for TLS')
hold off
print -depsc lcurvetls.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Print the "best" solution and the corresponding 
% relative counts.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(' ')
disp('The Best(??) Computed Solution: TLS, cutoff=15')
disp('Bin No.  Bin Center  Est. Count')
disp([[1:27]',bincenters',ffinal])
disp('Normalized counts for bins 7, 14, 17, and 20:')
disp(ffinal([7 14 17 20])/ffinal(7))