% Example computations related to one-dimensional deconvolution.
% Demonstration of truncated SVD solution applied to simulated noisy data
% involving no inverse crime.
%
% The following routines must be computed before this one:
% DC2_discretedata_comp.m
%
% Jennifer Mueller and Samuli Siltanen, October 2015

% Plot parameters
lwidth = .5;
thickline = 2;
fsize = 12;
msize = 4;

% Load precomputed results
load DC2_discretedata A x xx n m mn mIC sigma

% Choose data by commenting two unwanted lines among the three lines below.
% There is no inverse crime in any of these choices.
%data = m;   % No added noise
data = mn(:); % Small added noise

% Compute singular value decomposition (SVD)
[U,D,V] = svd(A);
svals = diag(D);

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demonstration of inversion by truncated SVD with several singular vectors
iiivec = ceil([.05 .1 .2 .4 .6 .8 1]*n);
for iii = iiivec
    
    % Plot current singular vector
    figure(2)
    clf
    subplot(2,1,1)
    plot(1:n,V(:,iii),'k','linewidth',lwidth)
    axis([1 n -.25 .3])
    set(gca,'ytick',[-.2 0 .2],'fontsize',fsize)
    set(gca,'xtick',[1 round(n/2) n],'fontsize',fsize)
    set(gca,'xticklabel',{})
    set(gca,'yticklabel',{})
    set(gca,'PlotBoxAspectRatio' , [2 1 1])
    box off
    title(['Singular vector number ', num2str(iii)])
    
    % Reconstruct from noisy data using only large singular values.
    ns               = iii;
    [row,col]        = size(D.');
    Dplus            = sparse(row,col);
    svals            = diag(D);
    Dplus(1:ns,1:ns) = diag(1./svals(1:ns));
    rec              = V*Dplus*U.'*data(:);
    % Calculate relative error
    truth  = DC_target(x);
    relerr = round(norm(rec(:)-truth(:))/norm(truth(:))*100);
    
    % Take a look at the reconstruction
    subplot(2,1,2)
    plot(xx,signal,'k','linewidth',lwidth)
    hold on
    plot(x,rec,'b','linewidth',lwidth)
    text(0.9,1.5,[num2str(relerr),'%'],'fontsize',fsize)
    
    %title(['SVD truncated at ', num2str(ns)])
    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
    title(['Reconstruction using ', num2str(iii), ' singular vectors'])
    
    % Plot the used and unused singular values
    figure(1)
    clf
    semilogy([1:n],svals,'k','linewidth',lwidth)
    hold on
    semilogy([1:iii],svals(1:iii),'r','linewidth',thickline)
    xlim([1 n])
    axis square
    title('Used (red) and unused (black) singular values')
    
    drawnow
    pause
end