clear all

n=1024;	% Signal length.

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

% Signal.
% Cosine signal
t = ((-n/2):(n/2-1))'/n;
%y0dct = cos(pi * (n/16) * t);
y0dct=zeros(n,1);
y0dct(1:n/4)=cos(2*pi*t(1:n/4)*32);
y0dct(n-n/4+1:n)=cos(2*pi*t(1:n/4)*32);
y0dct(n/2-n/8+1:n/2+n/8)=cos(2*pi*t(1:n/4)*16);
% Make the gaussian part.
y0wav=exp(-1E3.*t.^2)+exp(-2E3.*(t+0.25).^2)+exp(-1.5E3.*(t-0.25).^2);
y0 = y0dct + y0wav;

% Representation dictionary (here local DCT + DWT).
qmf=MakeONFilter('Symmlet',6);
%[qmf,dqmf]=MakeBSFilter('CDF',[4 4]);
dict1='RealFourier';pars11=n/8;pars12=0;pars13=0;
dict2='UDWT';pars21=2;pars22=qmf;pars23=0;
dict=MakeList(dict1,dict2);
pars1=MakeList(pars11,pars21);
pars2=MakeList(pars12,pars22);
pars3=MakeList(pars13,pars23);

% Coefficient vector length.
p = SizeOfDict1(n,dict,pars1,pars2,pars3);
tightFrame = p/n; % The dictionary is a tight frame with constant p/n.

% Observed noisy data.
SNR = 15;
sigma = std(y0)*10^(-SNR/20);
y = y0 + sigma*randn(size(y0));

% Solve Lasso in l1-penalized form with FISTA. 
% lambda chosen such that ||r|| <= epsilon.
options.niter = 300;
options.method = 'fista';
options.verb  = 0;
lambda = 0.45;
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;[xhatlassofista,Efista] = fb(zeros(p,1), ProxF, GradG, 1/tightFrame, options);timelassofista=toc

% Solve BP with DR.
options.gamma = 1;    % Stepsize parameter for Douglas-Rachford iteration.
options.tau   = 1;    % No relaxation.
options.niter = 300;
options.verb  = 0;
ProxF = @(x,ga)SoftThresh(x,ga);
ProxG = @(x,ga)x + KS(y0 - K(x))/tightFrame;
tic;[xhatBPDR,EBP] = dr(zeros(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 = 300;
options.verb  = 0;
epsilon = sqrt(n)*sigma*sqrt(1+2*sqrt(2/n));	      
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,1),ProxF,ProxG,options);timeBPDNDR=toc

% Solve Dantzig-Selector with PD.
options.mu    = .99*n/p;      % Stepsize on the dual.
options.nu    = .99*n/p;      % Stepsize on the primal.
options.tau   = 1;	      % No relaxation.
options.niter = 300;
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

% Solve with MCA.
tic;yhatMCA=SolveMCA(y,dict,pars1,MakeList(1/pars11,pars22),pars3,options.niter,0,3*sigma,epsilon,0);timeMCA=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 = xhatlassofista((ind+1):(ind+DimDomain));
	yhatlassofista(:,i) = eval(['Fast' NAME, 'Synthesis(C, n, PAR1, PAR2, PAR3)']);
	C = xhatBPDR((ind+1):(ind+DimDomain));
	yhatBPDR(:,i) = eval(['Fast' NAME, 'Synthesis(C, n, PAR1, PAR2, PAR3)']);
	C = xhatBPDNDR((ind+1):(ind+DimDomain));
	yhatBPDNDR(:,i) = eval(['Fast' NAME, 'Synthesis(C, n, PAR1, PAR2, PAR3)']);
	C = eval(['Fast' NAME, 'Analysis(squeeze(yhatMCA(i,:,:)), PAR1, PAR2, PAR3)']);
	xhatMCA((ind+1):(ind+DimDomain)) = C;
	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\t%-15s\t%-15s\t%-15s\t%-15s\t%-15s\n',' ','LASSO-FISTA','BP-DR','BPDN-DR','DS-PD','MCA');
fprintf('%-10s\t%-15g\t%-15g\t%-15g\t%-15g\t%-15g\n','||x||_0:',length(find(abs(xhatlassofista))), 	...
							      length(find(abs(xhatBPDR))), 		...
							      length(find(abs(xhatBPDNDR))), 		...
							      length(find(abs(xhatDSPD))), 		...
							      length(find(abs(xhatMCA))));

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot recovered signals.
figure;
subplot(6,1,1); 
plot([y0 y]);axis tight
legend('Noiseless',sprintf('Noisy SNR_{in}=%g dB',SNR));legend boxoff

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

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

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

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

subplot(6,3,16); 
plot([squeeze(sum(yhatMCA,1)) y0]);axis tight;
title(sprintf('MCA PSNR=%g dB',psnr(y0,squeeze(sum(yhatMCA,1)))));

% Plot DCT recovered component.
subplot(6,3,5); 
plot([yhatlassofista(:,1) y0dct]);axis tight
title(sprintf('LASSO-FISTA PSNR=%g dB',psnr(y0dct(:),yhatlassofista(:,1))));

subplot(6,3,8); 
plot([yhatBPDR(:,1) y0dct]);axis tight
title(sprintf('BP-DR PSNR=%g dB',psnr(y0dct(:),yhatBPDR(:,1))));

subplot(6,3,11); 
plot([yhatBPDNDR(:,1) y0dct]);axis tight
title(sprintf('BPDN-DR PSNR=%g dB',psnr(y0dct(:),yhatBPDNDR(:,1))));

subplot(6,3,14); 
plot([yhatDSPD(:,1) y0dct]);axis tight
title(sprintf('DS-PD PSNR=%g dB',psnr(y0dct(:),yhatDSPD(:,1))));

subplot(6,3,17); 
plot([squeeze(yhatMCA(1,:,:)) y0dct]);axis tight
title(sprintf('MCA PSNR=%g dB',psnr(y0dct(:),squeeze(yhatMCA(1,:,:)))));

% Plot Gaussian recovered component.
subplot(6,3,6); 
plot([yhatlassofista(:,2) y0wav]);axis tight
title(sprintf('LASSO-FISTA PSNR=%g dB',psnr(y0wav(:),yhatlassofista(:,2))));

subplot(6,3,9); 
plot([yhatBPDR(:,2) y0wav]);axis tight
title(sprintf('BP-DR PSNR=%g dB',psnr(y0wav(:),yhatBPDR(:,2))));

subplot(6,3,12); 
plot([yhatBPDNDR(:,2) y0wav]);axis tight
title(sprintf('BPDN-DR PSNR=%g dB',psnr(y0wav(:),yhatBPDNDR(:,2))));

subplot(6,3,15); 
plot([yhatDSPD(:,2) y0wav]);axis tight
title(sprintf('DS-PD PSNR=%g dB',psnr(y0wav(:),yhatDSPD(:,2))));

subplot(6,3,18); 
plot([squeeze(yhatMCA(2,:,:)) y0wav]);axis tight
title(sprintf('MCA PSNR=%g dB',psnr(y0wav(:),squeeze(yhatMCA(2,:,:)))));

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