clear all

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

load simu_sky.mat
n = length(x);

% PSF
load simu_psf.mat

% Normalize to unit mass and FFT.
h = h/sum(h(:));
H = fft2(h);
HS= conj(H);

% Dictionary (UDWT).
qmf=MakeONFilter('Symmlet',6);
dict='UDWT2';pars1=2;pars2=qmf;pars3=0;

% Coefficient vector length.
p = SizeOfDict2(n,dict,pars1,pars2,pars3);
be = (n*n)/p; % Inverse Lipschitz constant of the gradient, because ||H||=1.

% Observed noisy data.
BSNR   = 30;
y0 = real(ifft2(fft2(x).*H));
sigma = std(y0(:))*10^(-BSNR/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.1*sigma;
K  = @(x)real(ifft2(fft2(FastS2(x,n,dict,pars1,pars2,pars3)).*H));
KS = @(x)FastA2(real(ifft2(fft2(x).*HS)),dict,pars1,pars2,pars3);
ProxF = @(x,mu)SoftThresh(x,mu*lambda);
GradG = @(x)KS(K(x) - y);
%options.report = @(x)sum(sum((y-K(x)).^2))/2+lambda*norm(x,1);
tic;[xhatlassofista,Efista] = fb(zeros(p,1), ProxF, GradG, be, options);timelassofista=toc

% Solve BP with DR.
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 = 100;
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(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 = 100;
options.verb  = 0;
epsilon = 0.1*n*sigma;	      
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 = 500;
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);
KK  = @(x)KS(K(x));
KKS = KK;
GradH = @(x)0;
options.report = @(x,u)norm(x,1)+sum(z(:).*u(:))+delta*norm(u(:),1);
tic;[xhatDSPD,EDS] = primaldual(zeros(p,1), KK,  KKS, ProxFS, ProxG, GradH, options);timeDSPD=toc

% Compute recovered images.
yhatlassofista = FastS2(xhatlassofista,n,dict,pars1,pars2,pars3);
yhatBPPD = FastS2(xhatBPPD,n,dict,pars1,pars2,pars3);
yhatBPDNPD = FastS2(xhatBPDNPD,n,dict,pars1,pars2,pars3);
yhatDSPD = FastS2(xhatDSPD,n,dict,pars1,pars2,pars3);

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

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

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
% Display recovered images.
subplot(321); 
imagesc(x);axis image;axis off
title('Original image');

subplot(322); 
imagesc(y);axis image;axis off
title(sprintf('Degraded image: blurred and noisy BSNR=%g dB',BSNR));

subplot(323); 
imagesc(yhatlassofista);axis image;axis off
title(sprintf('Deconvolved LASSO-FISTA PSNR=%g dB',psnr(x(:),yhatlassofista(:))));

subplot(324); 
imagesc(yhatBPPD);axis image;axis off
title(sprintf('Deconvolved BP-PD PSNR=%g dB',psnr(x(:),yhatBPPD(:))));

subplot(325); 
imagesc(yhatBPDNPD);axis image;axis off
title(sprintf('Deconvolved BPDN-PD PSNR=%g dB',psnr(x(:),yhatBPDNPD(:))));

subplot(326); 
imagesc(yhatDSPD);axis image;axis off
title(sprintf('Deconvolved DS-PD PSNR=%g dB',psnr(x(:),yhatDSPD(:))));

colormap('hsv');

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