% Example computations related to one-dimensional deconvolution.
% Here we choose the regularization parameter using the L-curve method.
% Inversion is done by generalized Tikhonov regularization
% (derivative penalty) using stacked form.
%
% Samuli Siltanen March 2010

% Plot parameters
lwidth = .2;
redlwidth = .7;
fsize = 10;
msize = 7;
sigma = .01; % Noise level

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create data on finer grid to avoid inverse crime

% Construct signal
Nf  = 1499;
sf  = linspace(0,1,Nf);
Dsf = sf(2)-sf(1);
xf  = zeros(size(sf));
xf((sf>(1/4))&(sf<(1/2))) = 1;
xf((sf>(3/4))&(sf<(7/8))) = 2;

psf2 = conv2(ones(1,20),ones(1,20));
psf2 = conv2(psf2,psf2);
psf2 = psf2/sum(psf2);
xf = conv2(xf,psf2,'same')

%keyboard


% Construct normalized point spread function
a    = .05;
nPSFf = ceil(a/Dsf);
sPSFf = [-nPSFf:nPSFf]*Dsf;
PSFf = zeros(size(sPSFf));
PSFf((sPSFf>(-a))&(sPSFf<a)) = 1/(2*a);
PSFf = PSFf/(Dsf*sum(PSFf));

% Construct convolved signal
mf = Dsf*conv2(xf,PSFf,'same');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct objects needed for reconstruction
N  = 128;
s  = linspace(0,1,N);
Ds = s(2)-s(1);

% Construct normalized point spread function
nPSF = ceil(a/Ds);
sPSF = [-nPSF:nPSF]*Ds;
PSF = zeros(size(sPSF));
PSF((sPSF>(-a))&(sPSF<a)) = 1/(2*a);
PSF = PSF/(Ds*sum(PSF));

% Construct measurement matrix A
A = convmtx(PSF,N);
A = Ds*A(:,nPSF+1:end-nPSF);

% Construct measurements with and without noise
m = interp1(sf,mf,s,'linear');
m = m(:);
mn = m + sigma*randn(size(m));

% Construct continuum quantities for visual comparison
% Construct signal
Nc = 5000;
sc = linspace(-2,2,2*Nc+1);
Xc = zeros(size(sc));
Xc((sc>(1/4))&(sc<(1/2))) = 1;
Xc((sc>(3/4))&(sc<(7/8))) = 2;
PSFc = zeros(size(sc));
PSFc((sc>(-a))&(sc<a)) = 1/(2*a);
Mc = (sc(2)-sc(1))*conv2(Xc,PSFc,'full');
Mc = Mc(Nc+[1:length(sc)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demonstration of inversion by Tikhonov regularization

% Build discrete derivative matrix of size (N+1)x(N).
% This implements difference between consecutive values.
L = eye(N);
L = -[L;zeros(1,N)]+[zeros(1,N);L];
L = (1/Ds)*L;
L = L(2:(end-1),:);

%L = eye(N); %ALTERNATIVE: basic Tikhonov regularization

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

% Loop over regularization parameters and record corresponding norms
Lx_norms   = zeros(size(deltavec));
Ax_m_norms = zeros(size(deltavec));
for iii=1:length(deltavec)
    delta = deltavec(iii);
    % Build stacked form of the problem
    syst  = [A; sqrt(delta)*L];
    m_aug = [mn;zeros(N-1,1)];
    
    %m_aug = [mn;zeros(N,1)]; %ALTERNATIVE: basic Tikhonov regularization
    
    % Reconstruct from noisy data using Tikhonov regularization and delta as
    % regularization parameter
    recn = syst\m_aug;
    
    % Compute norms of the two competing terms in the regularization problem
    Lx_norms(iii)   = norm(L*recn);
    Ax_m_norms(iii) = norm(A*recn-mn);
    
    % Take a look at the reconstruction
    figure(2)
    clf
    plot(sc,Xc,'r','linewidth',redlwidth)
    hold on
    plot(s,recn,'k','linewidth',lwidth)
    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
    %pause
end

% Plot the L-curve
figure(3)
clf
plot(log(Ax_m_norms),log(Lx_norms))
title('L-curve analysis','fontsize',16)
xlabel('log ||Ax-m||','fontsize',16)
ylabel('log ||Lx||','fontsize',16)
hold on
% Choose the index of deltavec giving the delta value at the corner of the
% L-curve. Here we do this just visually, not automatically.
best_index = 90;
plot(log(Ax_m_norms(best_index)),log(Lx_norms(best_index)),'r*')

pause

% Compute "best" reconstruction
delta1 = deltavec(best_index);
% Build stacked form of the problem
syst  = [A;  sqrt(delta)*L];

% Reconstruct from noisy data using Tikhonov regularization and delta1 as
% regularization parameter. We use stacked form of the problem
syst  = [A; sqrt(delta1)*L];
m_aug = [mn(:);zeros(N-1,1)];
%m_aug = [mn;zeros(N,1)]; %ALTERNATIVE: basic Tikhonov regularization
rec1  = syst\m_aug;
x  = zeros(size(s));
x((s>(1/4))&(s<(1/2))) = 1;
x((s>(3/4))&(s<(7/8))) = 2;
relerr1 = round(norm(rec1(:)-x(:))/norm(x(:))*100);

% Take a look at the reconstruction
figure(7)
clf
plot(sc,Xc,'r','linewidth',redlwidth)
hold on
plot(s,rec1,'k','linewidth',lwidth)
title(['\delta=', num2str(delta1), ', error=', num2str(relerr1), '%'],'fontsize',fsize)
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

