%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  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.
%
%  This program shows how the data was generated.
%  It writes file spectdata.mat,
%  containing  K, g, nbins, bincenters, bnds,
%  and truesoln.mat, containing true_count.
%
%  makedata.m   11/2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all
randn('state',0)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The true solution: 
% the 4 energies are in center, and
% the relative peak heights are in proportion.
% The array sd controls the peak widths.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

center =     [2.54, 3.53, 3.85, 3.25];
proportion = [1,    2,    1,    1.5];
sd = [.1 .01, .05, .05];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Particles have energies true_energy, centered around
% the 4 energies, but with a normal distribution.
% We use 5000 particles per proportion.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ninc = 5000;
npart = proportion * ninc;
breaks = [1,npart(1),sum(npart(1:2)), ...
            sum(npart(1:3)),sum(npart(1:4))];

smear = .1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct the true energies.
% true_energy(i) = energy of particle i
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

true_energy(1:breaks(2)) = sd(1)*randn(npart(1),1)  ...
            + center(1);
for i=2:length(center),
  true_energy(breaks(i)+1:breaks(i+1)) =  ...
         sd(1)*randn(npart(i),1) + center(i);
end

nparticles = length(true_energy);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate observed energies: blurring by normal(1,smear)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

observed_energy = true_energy+smear*randn(1,nparticles);

minenergy = floor(10*min(observed_energy))/10 - .1;
maxenergy =  ceil(10*max(observed_energy))/10 + .1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now put the particles in energy bins, incrementing
% once per 100 particles.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

eincrement = .1;
nbins = round(10*(maxenergy-minenergy));
bnds = minenergy:.1:maxenergy;
bincenters = bnds(1:nbins) + eincrement/2;

for i=1:nbins,
  true_count(i,1) = ...
    sum((true_energy> bnds(i))&(true_energy <= bnds(i+1)));
  count(i,1) = ...
    sum((observed_energy > bnds(i))&(observed_energy <= bnds(i+1)));
end

count = floor(count/100)/10;
true_count = true_count/1000;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct the matrix K.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

normalrow = normpdf([0:eincrement:(nbins-1)*eincrement],0,eincrement);
K = toeplitz(normalrow,normalrow);

for i=1:nbins,
  K(:,i) = K(:,i)/sum(K(:,i));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Cut off 1st 2 and last 3 columns, where counts
% are assumed to be zero.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

K = K(:,3:nbins-3);
g = count;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Save the data and the true solution, and plot
% the data, as in Figure 1 of the homework.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

save spectdata.mat K g nbins bincenters bnds
save truesoln.mat true_count

figure(1)
plot(bincenters,count)
xlabel('energy')
ylabel('thousands of particles')
print -depsc data2.eps