clear all

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

% Real-life data.
load ima_Jupiter1_256.mat
n = length(y);
sigma = 1.8; % This is the std for Jupiter Deep Impact image (estimated by the MAD).

% PSF
load psf_star21.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);
tightFrame = p/(n*n);
be = (n*n)/p; % Inverse Lipschitz constant of the gradient, because ||H||=1.

% Solve Lasso + positivity constraint on \Phi x with GFB. 
% lambda chosen such that ||r|| <= epsilon.
options.mu    = be;		% Descent step-size for FB.
options.tau   = 1;		% No relaxation.
options.niter = 100;
options.verb  = 0;
lambda = 0.1*sigma;
Phi  = @(x)FastS2(x,n,dict,pars1,pars2,pars3);
PhiS = @(x)FastA2(x,dict,pars1,pars2,pars3);
K  = @(x)real(ifft2(fft2(Phi(x)).*H));
KS = @(x)PhiS(real(ifft2(fft2(x).*HS)));
ProxFi{1} = @(x,mu)SoftThresh(x,mu*lambda);
ProxFi{2} = prox_F_tightframe(@(x,ga)max(x,0),Phi,PhiS,tightFrame,0);
GradG = @(x)KS(K(x) - y);
%options.report = @(x)sum(sum((y-K(x)).^2))/2+lambda*norm(x,1);
tic;[xhatlassogfb,Egfb] = gfb(zeros(p,2), GradG, ProxFi, options);timelassogfb=toc

% Solve BPDN + positivity constraint on \Phi x with PD.
options.mu    = .99*sqrt(be/2);	% Stepsize on the dual.
options.nu    = .99*sqrt(be/2);	% Stepsize on the primal.
options.tau   = 1;		% No relaxation.
options.niter = 100;
options.verb  = 0;
epsilon = 3*n*sigma;
K  = @(x)[real(ifft2(fft2(Phi(x)).*H)) Phi(x)];
KS = @(x)PhiS(real(ifft2(fft2(x(:,1:n)).*HS)))+PhiS(x(:,n+1:end));	      
ProxG = @(x,mu)SoftThresh(x,mu);
ProxF = @(x,nu)[y+compute_l2ball_projection(x(:,1:n)-y,epsilon) max(x(:,n+1:end),0)];
ProxFS= prox_conjugate(ProxF);
GradH = @(x)0;
%options.report = @(x,u)norm(x,1)+sum(y(:).*u(:))+epsilon*norm(u(:));
tic;[xhatBPDNPD,EBPDN,u] = primaldual(zeros(p,1), K,  KS, ProxFS, ProxG, GradH, options);timeBPDNPD=toc

% Solve Dantzig-Selector + positivity constraint on \Phi x with PD.
options.mu    = .99*be/sqrt(be+1);      % Stepsize on the dual.
options.nu    = .99*be/sqrt(be+1);      % 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);
K  = @(x)real(ifft2(fft2(Phi(x)).*H));
KS = @(x)PhiS(real(ifft2(fft2(x).*HS)));
z = KS(y);
KK  = @(x)[KS(K(x)) x];
KKS = @(x)KS(K(x(:,1)))+x(:,2);
ProxF = @(x,nu)[z+compute_linfball_projection(x(:,1)-z,delta) x(:,2)+PhiS(max(Phi(x(:,2)),0)-Phi(x(:,2)))/tightFrame];
ProxFS= prox_conjugate(ProxF);
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.
yhatlassogfb = Phi(xhatlassogfb);
yhatBPDNPD = Phi(xhatBPDNPD);
yhatDSPD = Phi(xhatDSPD);

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

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

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
% Display recovered images.
subplot(221); 
imagesc(y);axis image;axis off
title('Real image');

subplot(222); 
imagesc(yhatlassogfb);axis image;axis off
title('Deconvolved LASSO+Pos-GFB');

subplot(223); 
imagesc(yhatBPDNPD);axis image;axis off
title('Deconvolved BPDN+Pos-PD');

subplot(224); 
imagesc(yhatDSPD);axis image;axis off
title('Deconvolved DS+Pos-PD');

colormap('gray');

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