clear all

n=1024;	% Signal length.

% Signal.
k = 20;
x = SparseVector(n, k, 'GAUSSIAN', true);
I = find(x);

% Titles: 'book' (subfigures), or 'toolbox' (meaningful titles).
titles = 'book';

% Save figures: 1/0 (default 0).
savfig = 0;

% PSF
options.param = 0.9;
h = psfgenerate(3,n,options);

% Normalize to unit mass and FFT.
h = h/sum(h(:));
H = fft(h);
HS= conj(H);
  
% Observed noisy data.
SNR = 25;
y0 = real(ifft(fft(x).*H));
sigma = std(y0(:))*10^(-SNR/20);
y = y0 + sigma*randn(size(y0));

be = 1; % 1/max(abs(H)^2).

% Solve Lasso in l1-penalized form with FB and FISTA. 
% lambda chosen such that ||r|| <= epsilon.
options.mu    = be;	      	% Descent step-size for FB.
options.tau   = 1;		% No relaxation -> 1.
options.method= 'fb';
options.niter = 200;
options.verb  = 0;
lambda = 0.01;
K  = @(x)real(ifft(fft(x).*H));
KS = @(x)real(ifft(fft(x).*HS));
ProxF = @(x,mu)SoftThresh(x,mu*lambda);
GradG = @(x)KS(K(x) - y);
options.report = @(x)norm(y-K(x))^2/2+lambda*norm(x,1);
tic;[xhatlassofb,Efb] = fb(zeros(n,1), ProxF, GradG, be, options);timelassofb=toc
options.tau   = (4-options.mu/be)/2;  % Relaxed FB.
tic;[xhatlassofbrlx,Efbrlx] = fb(zeros(n,1), ProxF, GradG, be, options);timelassofb=toc
options.method = 'fista';
tic;[xhatlassofista,Efista] = fb(zeros(n,1), ProxF, GradG, be, options);timelassofista=toc

% Solve BP with PD.
options.mu    = .99*sqrt(be);	% Stepsize on the dual.
options.nu    = .99*sqrt(be);	% Stepsize on the primal.
options.tau   = 1;		% No relaxation.
options.niter = 200;
options.verb  = 0;
ProxG = @(x,mu)SoftThresh(x,mu);
ProxF = @(x,nu)y0;
ProxFS= prox_conjugate(ProxF);
GradH = @(x)0;
options.report = @(x,u)norm(x,1)+sum(y0.*u);
tic;[xhatBPPD,EBP] = primaldual(zeros(n,1), K,  KS, ProxFS, ProxG, GradH, options);timeBPPD=toc

% Solve BPDN with PD.
options.mu    = .99*sqrt(be);	% Stepsize on the dual.
options.nu    = .99*sqrt(be);	% Stepsize on the primal.
options.tau   = 1;		% No relaxation.
options.niter = 200;
options.verb  = 0;
epsilon = sqrt(n)*sigma*sqrt(1+2*sqrt(2/n));	      
ProxG = @(x,mu)SoftThresh(x,mu);
ProxF = @(x,nu)y+compute_l2ball_projection(x-y,epsilon);
ProxFS= prox_conjugate(ProxF);
GradH = @(x)0;
options.report = @(x,u)norm(x,1)+sum(y.*u)+epsilon*norm(u);
tic;[xhatBPDNPD,EBPDN] = primaldual(zeros(n,1), K,  KS, ProxFS, ProxG, GradH, options);timeBPDNPD=toc

% Solve Dantzig-Selector with PD.
options.mu    = .99*be;      % Stepsize on the dual.
options.nu    = .99*be;      % Stepsize on the primal.
options.tau   = 1;	     % No relaxation.
options.niter = 200;
options.verb  = 0;
delta = lambda;       % Same parameter as Lasso.      
ProxG = @(x,mu)SoftThresh(x,mu);
z = KS(y);
ProxF = @(x,nu)z+compute_linfball_projection(x-z,delta);
ProxFS= prox_conjugate(ProxF);
K  = @(x)real(ifft(fft(x).*(abs(H).^2)));
KS = K;
GradH = @(x)0;
options.report = @(x,u)norm(x,1)+sum(z.*u)+delta*norm(u,1);
tic;[xhatDSPD,EDS] = primaldual(zeros(n,1), K,  KS, ProxFS, ProxG, GradH, options);timeDSPD=toc

fprintf('%s\n','*'*ones(1,90));
fprintf('%40sSummary\n',' ');
fprintf('%s\n','*'*ones(1,90));
fprintf('%-10s%-15s\t%-15s\t%-15s\t%-15s\t%-15s\t%-15s\n',' ','Original','LASSO-ISTA','LASSO-FISTA','BP-DR','BPDN-DR','DS-PD');
fprintf('%-10s%-15g\t%-15g\t%-15g\t%-15g\t%-15g\t%-15g\n','||x||_0:',length(find(abs(x))),		 ...
							      	     length(find(abs(xhatlassofb))), 	 ...
							      	     length(find(abs(xhatlassofista))),  ...
							      	     length(find(abs(xhatBPPD))),	 ...
							      	     length(find(abs(xhatBPDNPD))),	 ...
							      	     length(find(abs(xhatDSPD))));   

fprintf('%-10s%-15g\t%-15g\t%-15g\t%-15g\t%-15g\t%-15g\n','||x||_1:',norm(x,1),			...
							      	     norm(xhatlassofb,1), 	...
							      	     norm(xhatlassofista,1), 	...
							      	     norm(xhatBPPD,1), 		...
							      	     norm(xhatBPDNPD,1),	...
							      	     norm(xhatDSPD,1));

fprintf('%-25s\t%-15g\t%-15g\t%-15g\t%-15g\t%-15g\n','CPU (s):', timelassofb, timelassofista, timeBPPD, timeBPDNPD, timeDSPD);

figure
t = [1:n];
subplot(3,1,1); 
plot(t,y);axis tight

if strcmp(titles,'book'),     
	h=title(['(a)']);set(h,'fontsize',18);
else			      
	title(sprintf('Noisy blurred spikes PSNR_{in}=%g dB',PSNR));
end

subplot(3,3,4);
stem(t(xhatlassofb~=0),xhatlassofb(xhatlassofb~=0),'.b');hold on;stem(t(I),x(I),'or');axis([1 n min(x)-0.01 max(x)+0.01]);hold off
if strcmp(titles,'book'),     h=title(['(b)']);set(h,'fontsize',18);
else			      title(sprintf('LASSO-ISTA PSNR=%g dB',psnr(x(:),xhatlassofb(:))));
end

subplot(3,3,5);
stem(t(xhatlassofista~=0),xhatlassofista(xhatlassofista~=0),'.b');hold on;stem(t(I),x(I),'or');axis([1 n min(x)-0.01 max(x)+0.01]);hold off
if strcmp(titles,'book'),     h=title(['(c)']);set(h,'fontsize',18);
else			      title(sprintf('LASSO-FISTA PSNR=%g dB',psnr(x(:),xhatlassofista(:))));
end

subplot(3,3,6);
plot(1:options.niter,Efb,'-k',1:options.niter,Efbrlx,'-.k',1:options.niter,Efista,'--k');axis([1 options.niter min(Efista) max(Efista)])
xlabel('Iteration');ylabel('Objective function');
legend('ISTA (no relaxation)','ISTA (over-relaxation)','FISTA');legend boxoff
if strcmp(titles,'book'),     h=title(['(d)']);set(h,'fontsize',18); end

subplot(3,3,7); 
stem(t(xhatBPPD~=0),xhatBPPD(xhatBPPD~=0),'.b');hold on;stem(t(I),x(I),'or');axis([1 n min(x)-0.01 max(x)+0.01]);hold off
if strcmp(titles,'book'),     h=title(['(e)']);set(h,'fontsize',18);
else			      title(sprintf('BP-DR PSNR=%g dB',psnr(x(:),xhatBPPD(:))));;
end

subplot(3,3,8);
stem(t(xhatBPDNPD~=0),xhatBPDNPD(xhatBPDNPD~=0),'.b');hold on;stem(t(I),x(I),'or');axis([1 n min(x)-0.01 max(x)+0.01]);hold off
if strcmp(titles,'book'),     h=title(['(f)']);set(h,'fontsize',18);
else			      title(sprintf('BPDN-DR PSNR=%g dB',psnr(x(:),xhatBPDNPD(:))));
end

subplot(3,3,9);
stem(t(xhatDSPD~=0),xhatBPDNPD(xhatDSPD~=0),'.b');hold on;stem(t(I),x(I),'or');axis([1 n min(x)-0.01 max(x)+0.01]);hold off
if strcmp(titles,'book'),     h=title(['(g)']);set(h,'fontsize',18);
else			      title(sprintf('DS-PD PSNR=%g dB',psnr(x(:),xhatDSPD(:))));
end

if savfig,
  saveas(gcf,'1D/Datasets/testsSparseSpikeDeconv1D.fig','fig');
  saveas(gcf,'1D/Datasets/testsSparseSpikeDeconv1D.eps','epsc');
end