clear all

n=64;   % Number of measurements n x n.
p=n*2;  % Coefficient image length.
k=50;   % Sparsity level.

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

% Stricly sparse image.
x = SparseVector(p*p, k, 'GAUSSIAN', true);

% Random measurement (sensing) operator: Hadamard, Fourier, Real Fourier (RST), USE, etc.
dict = 'RST';
tightFrame = 1;  % For Hadamard, RST and Fourier (Parseval tight frames), (Phi Phi') = I_n. Otherwise, if unknown or frame, set to 0.
q = randperm(p*p);
Omega = q(1:n*n)';
be = 1/tightFrame;

% Observed noisy data.
PSNR = 20;
y0 = FastMeasure2D(x, dict, Omega);
sigma = max(abs(y0))*10^(-PSNR/20);
y = y0 + sigma*randn(size(y0));

% Solve Lasso in l1-penalized form with FISTA. 
% lambda chosen such that ||r|| <= epsilon.
options.mu    = be;		% Descent step-size for FB.
options.method = 'fista';
options.niter = 100;
options.verb  = 0;
lambda = 0.03;
K  = @(x)FastMeasure2D(x, dict, Omega);
KS = @(x)FastMeasureAdjoint2D(x, p*p, dict, Omega);
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;[xhatlassofista,Efista] = fb(zeros(p*p,1), ProxF, GradG, be, options);timelassofista=toc

% Solve BP with DR.
options.gamma = 0.1;    % Stepsize parameter for Douglas-Rachford iteration.
options.tau   = 1;    % No relaxation.
options.niter = 100;
options.verb  = 0;
ProxF = @(x,ga)SoftThresh(x,ga);
ProxG = @(x,ga)x + KS(y0 - K(x))/tightFrame;
tic;[xhatBPDR,EBP] = dr(zeros(p*p,1),ProxF,ProxG,options);timeBPDR=toc

% Solve BPDN with DR.
options.gamma = 0.1;    % Stepsize parameter for Douglas-Rachford iteration.
options.tau   = 1;    % No relaxation.
options.niter = 100;
options.verb  = 0;
%epsilon = n*sigma*sqrt(1+2*sqrt(2/(n*n)));	      
epsilon = n*sigma;
ProxF = @(x,ga)SoftThresh(x,ga);
ProxG = prox_F_tightframe(@(x,ga)compute_l2ball_projection(x,epsilon),K,KS,tightFrame,y);
tic;[xhatBPDNDR,EBPDN] = dr(zeros(p*p,1),ProxF,ProxG,options);timeBPDNDR=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 = 100;
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)FastMeasureAdjoint2D(FastMeasure2D(x, dict, Omega), p*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*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\n',' ','Original','LASSO-FISTA','BP-DR','BPDN-DR','DS-PD');
fprintf('%-10s%-15g\t%-15g\t%-15g\t%-15g\t%-15g\n','||x||_0:',length(find(abs(x))),		...
							      length(find(abs(xhatlassofista))), ...
							      length(find(abs(xhatBPDR))),	...
							      length(find(abs(xhatBPDNDR))),	...
							      length(find(abs(xhatDSPD))));   

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

fprintf('%-25s\t%-15g\t%-15g\t%-15g\t%-15g\n','CPU (s):', timelassofista, timeBPDR, timeBPDNDR, timeDSPD);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
subplot(231)
imagesc(reshape(real(y(1:n*n)),n,n));axis image;axis off
title(sprintf('Projections %s matrix PSNR=%g dB',dict,PSNR));

subplot(232)
imagesc(reshape(x,p,p));axis image;axis off
title(sprintf('Original image'));

subplot(233)
imagesc(reshape(xhatlassofista,p,p));axis image;axis off
title(sprintf('LASSO-FISTA MSE=%e',sqrt(norm(xhatlassofista-x).^2/(p*p))));

subplot(234)
imagesc(reshape(xhatBPDR,p,p));axis image;axis off
title(sprintf('BP-DR (noiseless) MSE=%e',sqrt(norm(xhatBPDR-x).^2/(p*p))));

subplot(235)
imagesc(reshape(xhatBPDNDR,p,p));axis image;axis off
title(sprintf('BPDN-DR MSE=%e',sqrt(norm(xhatBPDNDR-x).^2/(p*p))));

subplot(236)
imagesc(reshape(xhatDSPD,p,p));axis image;axis off
title(sprintf('DS-PDF MSE=%e',sqrt(norm(xhatDSPD-x).^2/(p*p))));

colormap('gray')

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