%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Solution to Problem 6
%  Blind Deconvolution: A Matter of Norm
%  Dianne P. O'Leary 
%  Computing in Science and Engineering
%
%  In this problem, we use total
%  least squares (TLN) 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.
%
%  The norm used for minimization is either the infinity-norm
%  or the one-norm.
%
%  problem6.m   11/2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Modified 12-2008 to label axes in plots.

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

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

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

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Using the Infinity-norm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% lambdas don't regularize this one well

for lambda = [0 .06]
      [f,e,r,itn] = stlninfty(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]; 
      subplot(2,2,nlambda)
      plot(bincenters,count,'r.', bincenters,fpad,'b','LineWidth',1.5)
      title(sprintf('Infinity norm, lambda = %5.2f',lambda))
      xlabel('energy')
      ylabel('particle count')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Using the One-norm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for lambda = [.02 .06 ]
      [f,e,r,itn] = stlnone(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];
      subplot(2,2,nlambda)
      plot(bincenters,count,'r.', bincenters,fpad,'b','LineWidth',1.5)
      title(sprintf('1-norm, lambda = %5.2f',lambda))
      xlabel('energy')
      ylabel('particle count')
end

print -depsc stln.eps

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

fpad = [0;0;f;0;0;0];
disp(' ')
disp('The Best(??) Computed Solution: TLN, lambda=0.06')
disp('Bin No.  Bin Center  Est. Count')
disp([[1:27]',bincenters',fpad])
disp('Normalized counts for bins 7, 14, 17, and 20:')
disp(fpad([7 14 17 20])/fpad(7))
disp('Normalized counts for bins 7, 14, 16, and 21:')
disp(fpad([7 14 16 21])/fpad(7))