clear all

n=512;	% Signal length.

% Signal.
y0 = InputSignal('Star');
y0 = y0(1:n);

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

% Representation dictionary: overcomplete DCT/DST+DIRAC.
% The dictrionary is not a tight frame because of the DST.
fineness = 2;
dict  = MakeList('DCT','DST','DIRAC');
pars1 = MakeList(fineness,fineness,0);
pars2 = MakeList(0,0,0);
pars3 = MakeList(0,0,0);

% Coefficient vector length.
p = SizeOfDict1(n,dict,pars1,pars2,pars3);
be= n/p; % 1/||\Phi||^2. 

% Observation vector.
PSNR = 50;
sigma = max(abs(y0))*10^(-PSNR/20);
y = y0 + sigma*randn(size(y0));

% 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 = 500;
options.verb  = 0;
lambda = 0.253;
K  = @(x)FastS1(x,n,dict,pars1,pars2,pars3);
KS = @(x)FastA1(x,dict,pars1,pars2,pars3);
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 = 500;
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 = 500;
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 = 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);
K  = @(x)FastA1(FastS1(x,n,dict,pars1,pars2,pars3),dict,pars1,pars2,pars3);
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

% Compute recovered signals and respective morphological components.
NumberOfDicts = LengthList(dict);
ind = 0;
for i = 1:NumberOfDicts,
	NAME = NthList(dict, i);
	PAR1 = NthList(pars1, i);
	PAR2 = NthList(pars2, i);
	PAR3 = NthList(pars3, i);
	DimDomain = SizeOfDict1(n, NAME, PAR1, PAR2, PAR3);
	C = xhatlassofb((ind+1):(ind+DimDomain));
	yhatlassofb(:,i) = eval(['Fast' NAME, 'Synthesis(C, n, PAR1, PAR2, PAR3)']);
	C = xhatlassofista((ind+1):(ind+DimDomain));
	yhatlassofista(:,i) = eval(['Fast' NAME, 'Synthesis(C, n, PAR1, PAR2, PAR3)']);
	C = xhatBPPD((ind+1):(ind+DimDomain));
	yhatBPPD(:,i) = eval(['Fast' NAME, 'Synthesis(C, n, PAR1, PAR2, PAR3)']);
	C = xhatBPDNPD((ind+1):(ind+DimDomain));
	yhatBPDNPD(:,i) = eval(['Fast' NAME, 'Synthesis(C, n, PAR1, PAR2, PAR3)']);
	C = xhatDSPD((ind+1):(ind+DimDomain));
	yhatDSPD(:,i) = eval(['Fast' NAME, 'Synthesis(C, n, PAR1, PAR2, PAR3)']);
	ind = ind + DimDomain;
end

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',' ','LASSO-ISTA','LASSO-FISTA','BP-PD','BPDN-PD','DS-PD');
fprintf('%-10s%-15g\t%-15g\t%-15g\t%-15g\t%-15g\n','||x||_0:',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\n','||x||_1:',norm(xhatlassofb,1),	...
							      norm(xhatlassofista,1),	...
							      norm(xhatBPPD,1), 	...
							      norm(xhatBPDNPD,1),	...
							      norm(xhatDSPD,1));
						     
fprintf('%-10s%-15g\t%-15g\t%-15g\t%-15g\t%-15g\n','CPU (s)', timelassofb, timelassofista, timeBPPD, timeBPDNPD, timeDSPD);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;clf					       
% Plot recovered signals.
subplot(6,1,1); 
plot(y);axis tight
title(sprintf('Noisy signal'));

subplot(6,3,4); 
plot([sum(yhatlassofb,2) y0]);axis tight;
title(sprintf('LASSO-ISTA PSNR=%g dB',psnr(y0,sum(yhatlassofb,2))));

subplot(6,3,7); 
plot([sum(yhatlassofista,2) y0]);axis tight;
title(sprintf('LASSO-FISTA PSNR=%g dB',psnr(y0,sum(yhatlassofista,2))));

subplot(6,3,10); 
plot([sum(yhatBPPD,2) y0]);axis tight;
title(sprintf('BP-PD PSNR=%g dB',psnr(y0,sum(yhatBPPD,2))));

subplot(6,3,13); 
plot([sum(yhatBPDNPD,2) y0]);axis tight;
title(sprintf('BPDN-PD PSNR=%g dB',psnr(y0,sum(yhatBPDNPD,2))));

subplot(6,3,16); 
plot([sum(yhatBPDNPD,2) y0]);axis tight;
title(sprintf('DS-PD PSNR=%g dB',psnr(y0,sum(yhatDSPD,2))));

% Plot DCT-DST recovered component.
subplot(6,3,5); 
plot(sum(yhatlassofb(:,1:2),2));axis tight;
title(sprintf('DCT-DST component LASSO-ISTA'));

subplot(6,3,8); 
plot(sum(yhatlassofista(:,1:2),2));axis tight;
title(sprintf('DCT-DST component LASSO-FISTA'));

subplot(6,3,11); 
plot(sum(yhatBPPD(:,1:2),2));axis tight;
title(sprintf('DCT-DST component BP-PD'));

subplot(6,3,14); 
plot(sum(yhatBPDNPD(:,1:2),2));axis tight;
title(sprintf('DCT-DST component BPDN-PD'));

subplot(6,3,17); 
plot(sum(yhatDSPD(:,1:2),2));axis tight;
title(sprintf('DCT-DST component DS-PD'));

% Plot energies.
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

subplot(3,3,9); 
loglog(1:options.niter,EBP/EBP(2),'-k',1:options.niter,EBPDN/EBPDN(2),'--b',1:options.niter,EDS/EDS(2),'-.r',1:options.niter,1./(1:options.niter),':k');axis tight
xlabel('Iteration');ylabel('Ergodic partial duality gap');
legend('BP-PD','BPDN-PD','DS-PD','1/iteration');legend boxoff

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;clf
% Plot recovered coeffs.
subplot(5,3,1); 
plot(xhatlassofb);axis([1 p min(xhatlassofb) max(xhatlassofb)]);
title(sprintf('Coeffs of LASSO-ISTA'));

subplot(5,3,4); 
plot(xhatlassofista);axis([1 p min(xhatlassofista) max(xhatlassofista)]);
title(sprintf('Coeffs of LASSO-FISTA'));

subplot(5,3,7); 
plot(xhatBPPD);axis([1 p min(xhatBPPD) max(xhatBPPD)]);
title(sprintf('Coeffs of BP-PD'));

subplot(5,3,10); 
plot(xhatBPDNPD);axis([1 p min(xhatBPDNPD) max(xhatBPDNPD)]);
title(sprintf('Coeffs of BPDN-PD'));

subplot(5,3,13); 
plot(xhatDSPD);axis([1 p min(xhatDSPD) max(xhatDSPD)]);
title(sprintf('Coeffs of DS-PD'));

% Plot DCT-DST recovered coeffs.
subplot(5,3,2); 
plot(xhatlassofb(1:2*n*fineness));
axis tight;
axis1 = axis;
axis([1 2*n*fineness axis1(3) axis1(4)]);
title(sprintf('DCT-DST coeffs LASSO-ISTA'));

subplot(5,3,5); 
plot(xhatlassofista(1:2*n*fineness));
axis tight;
axis1 = axis;
axis([1 2*n*fineness axis1(3) axis1(4)]);
title(sprintf('DCT-DST coeffs LASSO-FISTA'));

subplot(5,3,8); 
plot(xhatBPPD(1:2*n*fineness));
axis tight;
axis1 = axis;
axis([1 2*n*fineness axis1(3) axis1(4)]);
title(sprintf('DCT coeffs BP-PD'));

subplot(5,3,11); 
plot(xhatBPDNPD(1:2*n*fineness));
axis tight;
axis1 = axis;
axis([1 2*n*fineness axis1(3) axis1(4)]);
title(sprintf('DCT coeffs BPDN-PD'));

subplot(5,3,14); 
plot(xhatDSPD(1:2*n*fineness));
axis tight;
axis1 = axis;
axis([1 2*n*fineness axis1(3) axis1(4)]);
title(sprintf('DCT coeffs DS-PD'));

% Plot Dirac recovered coeffs.
subplot(5,3,3); 
plot(xhatlassofb(end-n+1:end));
axis tight;
axis1 = axis;
axis([1 n axis1(3) axis1(4)]);
title(sprintf('Dirac coeffs LASSO-ISTA'));

subplot(5,3,6); 
plot(xhatlassofista(end-n+1:end));
axis tight;
axis1 = axis;
axis([1 n axis1(3) axis1(4)]);
title(sprintf('Dirac coeffs lLASSO-FISTA'));

subplot(5,3,9); 
plot(xhatBPPD(end-n+1:end));
axis tight;
axis1 = axis;
axis([1 n axis1(3) axis1(4)]);
title(sprintf('Dirac coeffs BP-PD'));

subplot(5,3,12); 
plot(xhatBPDNPD(end-n+1:end));
axis tight;
axis1 = axis;
axis([1 n axis1(3) axis1(4)]);
title(sprintf('Dirac coeffs BPDN-PD'));

subplot(5,3,15); 
plot(xhatDSPD(end-n+1:end));
axis tight;
axis1 = axis;
axis([1 n axis1(3) axis1(4)]);
title(sprintf('Dirac coeffs DS-PD'));

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

tilefigs