clear all

n=1024;	% Number of measurements.
p=n*4;	% Coefficient vector length.
k=20;	% Sparsity level.

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

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

% Stricly or compressible signal.
Compress = 0;
if ~Compress,
	x = SparseVector(p, k, 'GAUSSIAN', true);
	I = find(x);
else
	s = 0.7;
	x = SparseVector(p, k, 'Signs', true).*(1:p)'.^(-1/s);
	x = x(randperm(p));
	x = x/max(abs(x));
	I = find(x);
end

% Random measurement (sensing) operator: Hadamard, Fourier, RealFourier, RealSinusoid (RST).
dict = 'RST';
q = randperm(p);
Omega = q(1:n)';
tightFrame = 1;  % e.g. for Hadamard, Fourier and RST (tight frames), (Phi Phi') = I_n. Otherwise, if unknown or frame, set to 0.
be = 1/tightFrame;

% Observed noisy data.
SNR = 40;
y0 = FastMeasure(x, dict, Omega);
sigma = std(y0)*10^(-SNR/20);
y = y0 + sigma*randn(n,1);

% 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;
%R = Rlist(is);
lambda = 0.005;
K  = @(x)FastMeasure(x, dict, Omega);
KS = @(x)FastMeasureAdjoint(x, p, dict, Omega);
%ProxF = @(x,mu)compute_l1ball_projection(x,R);
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(p,1), ProxF, GradG, be, options);timelassofb=toc
options.tau   = (4-options.mu/be)/2;  % Relaxed FB.
tic;[xhatlassofbrlx,Efbrlx] = fb(zeros(p,1), ProxF, GradG, be, options);timelassofb=toc
options.method = 'fista';
tic;[xhatlassofista,Efista] = fb(zeros(p,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(y.*u);
tic;[xhatBPPD,EBP] = primaldual(zeros(p,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(p,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)FastMeasureAdjoint(FastMeasure(x, dict, Omega), p, dict, Omega);
KS = K;
GradH = @(x)0;
options.report = @(x,u)norm(x,1)+sum(z.*u)+delta*norm(u,1);
tic;[xhatDSPD,EDS] = primaldual(zeros(p,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-PD','BPDN-PD','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:p];
subplot(3,1,1); 
plot(y);axis tight
title(sprintf('Noisy measurements (%s) PSNR=%g dB',dict,SNR));axis tight

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-PD 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-PD 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/testsCS1implicit.fig','fig');
  saveas(gcf,'1D/Datasets/testsCS1implicit.eps','epsc');
end