%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Solution to Problem 4
%  Blind Deconvolution: A Matter of Norm
%  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 and K and E are Toeplitz.
%
%  We regularize the problem by controlling the norm
%  of the solution vector as well as the norm of the
%  residual r.
%
%  problem4.m   11/2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Modified 12-2008 to add axis labels to Figure 1.

close all
load spectdata
count = g;
[nrows,ncols]=size(K);

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

figure(1)
plot(bincenters,count,'LineWidth',1.5)
title('Observed counts')
xlabel('energy')
ylabel('particle count')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Total Least Squares with a Toeplitz constraint.
% The parameter lambda specifies the relative
% importance of keeping the norm of the solution vector
% small.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(' ')
disp('Total least squares solutions')
tol = 1.e-3;
figure(2)

nlambda = 0;
for lambda = [.02,.06,.16,.4]
      [f,e,r,itn] = stls(K,g,lambda,tol);
      nlambda = nlambda + 1;
      residnorm(nlambda) = norm(r);
      solnnorm(nlambda) = norm(f);
      disp(sprintf('Solution for lambda = %f',lambda))
      disp(sprintf(' residual norm = %f',residnorm(nlambda)))
      disp(sprintf(' solution norm = %f',solnnorm(nlambda)))
      disp(sprintf(' Number of Newton iterations = %d',itn))
      fpad = [0;0;f;0;0;0]; 
      if (nlambda == 2)
        ffinal = fpad;
      end
      subplot(2,2,nlambda)
      plot(bincenters,count,'r.', bincenters,fpad,'b','LineWidth',1.5)
      title(sprintf('2-norm, lambda = %5.2f',lambda))
      xlabel('energy')
      ylabel('particle count')
end

print -depsc stls.eps

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

disp(' ')
disp('The Best(??) Computed Solution: TLS, lambda=0.06')
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))
disp('Normalized counts for bins 7, 14, 16, and 21:')
disp(ffinal([7 14 16 21])/ffinal(7))