% Example computations related to one-dimensional deconvolution.
% Demonstration of Tikhonov regularization and Morozov discrepancy
% principle with noise level 1%.
%
% Samuli Siltanen February 2012

% Plot parameters
lwidth = .5;
thickline = 2;
fsizeM = 30;
fsize = 20;
msize = 16;

% Load data computed using a very fine grid, avoiding inverse crime
load data/DC_cont_sampled_data128 cdata1 sigma1 
n  = length(cdata1);
mn = cdata1;

% Load precomputed measurement matrix A
load data/DC_discr_data x2 A2
A = A2;
recx = x2;

% Construct the regularization matrix L.
% In this case of classical Tikhonov regularization it is just identity.
L = eye(n);

% Construct original signal for comparison
xx = linspace(0,1,1024);
signal = DC_target(xx);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demonstration of inversion by Tikhonov regularization
% We use Morozov discrepancy principle to find the regularization parameter.

% How many singular values are zero?
[U,D,V] = svd(A);
svals   = diag(D);
r       = max(find(svals>1e-10));

% Here we implement finding the zero in a naive but simple way.
mpilkku  = U.'*mn;
% See if the noise level satisfies the condition necessary for the use of
% the Morozov method
minimum = norm(mpilkku(r+1:end));
maximum = norm(mpilkku);
delta = sqrt(n)*sigma1;
disp([minimum,delta,maximum]);
if (norm(mpilkku(r+1:end))>delta) | (delta>norm(mpilkku))
    error('Morozov condition is not satisfied; cannot continue computation')
end
alphavec = [0.0000001:.00001:.01];
objfun   = zeros(size(alphavec));
dvec     = zeros(size(mpilkku));
for iii = 1:length(alphavec)
    dvec(1:r)     = alphavec(iii)./(alphavec(iii) + svals(1:r).^2).*mpilkku(1:r);
    dvec(r+1:end) = mpilkku(r+1:end);
    objfun(iii)   = dvec.'*dvec - delta^2;
end
alpha = alphavec(min(find(objfun>0)));

% Reconstruct from noisy data using Tikhonov regularization and alpha as
% regularization parameter
syst  = [A; sqrt(alpha)*L];
m_aug = [mn;zeros(n,1)];
recn = syst\m_aug;

% Take a look at the reconstruction
figure(2)
clf
plot(alphavec,objfun,'k','linewidth',lwidth)
hold on
plot([alphavec(1),alpha],[0,0],'k--','linewidth',lwidth)
plot([alpha,alpha],[0,-10],'k--','linewidth',lwidth)
ylim([1.2*min(objfun) 1.2*max(objfun)])
set(gca,'ytick',[-0.01 0 0.01],'fontsize',fsizeM)
set(gca,'xtick',[alpha max(alphavec)],'fontsize',fsizeM)
axis square
box off
print -depsc ../../images/DC_Morozov1.eps

% Calculate relative error
truth  = DC_target(recx);
relerr = round(norm(recn(:)-truth(:))/norm(truth(:))*100);

% Plot the result
figure(5)
clf
plot(xx,signal,'k','linewidth',lwidth)
hold on
plot(recx,recn,'k','linewidth',thickline)
text(0.9,1.5,[num2str(relerr),'%'],'fontsize',fsize)
axis([0 1 -.5 2])
set(gca,'xtick',[0 1/2 1],'fontsize',fsize)
set(gca,'ytick',[0  1 2],'fontsize',fsize)
%set(gca,'xticklabel',{})
%set(gca,'yticklabel',{})
set(gca,'PlotBoxAspectRatio' , [2 1 1])
box off
print -depsc ../../images/DC_Morozov2.eps

