% Example computations related to one-dimensional deconvolution.
% Here we choose the regularization parameter using the L-curve method.
% Inversion is done by Tikhonov regularization.
%
% Samuli Siltanen February 2014

% Parameters for controlling the appearance of plots
lwidth = .5;
thickline = 2;
fsizeL = 30;
fsize = 20;
msize = 32;

% Load data computed using a very fine grid, avoiding inverse crime
load DC2_discretedata A x xx n m mn mIC sigma

% Choose data
data = mn;

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

% Vector of regularization parameters
alphavec = 10.^[-8:.02 : 2];

% Loop over regularization parameters and record corresponding norms
Lx_norms   = zeros(size(alphavec));
Ax_m_norms = zeros(size(alphavec));
for iii=1:length(alphavec)
    alpha = alphavec(iii);
    
    % Build stacked form of the problem
    syst  = [A; sqrt(alpha)*L];
    m_aug = [data(:);zeros(n,1)];
    
    % Reconstruct from noisy data using Tikhonov regularization and alpha as
    % regularization parameter
    rec = syst\m_aug;
    
    % Compute norms of the two competing terms in the regularization problem
    Lx_norms(iii)   = norm(L*rec);
    Ax_m_norms(iii) = norm(A*rec-data(:));
    
    % Take a look at the reconstruction
    figure(2)
    clf
    plot(xx,DC_target(xx),'k','linewidth',lwidth)
    hold on
    plot(x,rec,'r','linewidth',thickline)
    axis([0 1 -.2 2.2])
    set(gca,'xtick',[0 1/4 1/2 3/4 7/8 1],'fontsize',fsize)
    set(gca,'xticklabel',{'0' '1/4' '1/2' '3/4' '7/8' '1'},'fontsize',fsize)
    set(gca,'ytick',[0  1 2],'fontsize',fsize)
    axis square
    box off
    drawnow
    
    % Monitor the run
    disp(['Done ', num2str(iii), ' out of ', num2str(length(alphavec))])
    %pause
end

% Plot the L-curve. We choose the index of alphavec giving the
% alpha value at the corner of the L-curve. Here we do this just
% visually, not automatically.
best_index = 230; % This choice should be automated
figure(3)
clf
endind = round(.8*length(Ax_m_norms));
plot(log(Ax_m_norms(1:endind)),log(Lx_norms(1:endind)),'k','linewidth',lwidth)
hold on
plot(log(Ax_m_norms(best_index)),log(Lx_norms(best_index)),'r.','markersize',msize)
axis([min(log(Ax_m_norms(1:endind))) max(log(Ax_m_norms(1:endind))) min(log(Lx_norms(1:endind))) max(log(Lx_norms(1:endind)))])
xt1 = min(log(Ax_m_norms(1:endind)));
xt2 = log(Ax_m_norms(best_index));
xt3 = max(log(Ax_m_norms(1:endind)));
yt1 = min(log(Lx_norms(1:endind)));
yt2 = log(Lx_norms(best_index));
yt3 = max(log(Lx_norms(1:endind)));
axis square
box off

% Compute "best" reconstruction
alpha = alphavec(best_index);
syst  = [A; sqrt(alpha)*L];
m_aug = [data(:);zeros(n,1)];
rec   = syst\m_aug;

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

% Plot the result
figure(5)
clf
plot(xx,DC_target(xx),'k','linewidth',lwidth)
hold on
plot(x,rec,'r','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,'PlotBoxAspectRatio' , [2 1 1])
box off

