clear all

%Signal list.
signals = {'Transients','Greasy','Carbon','WernerSorrows','FM','TwoChirp'};

% Read, input or make the signal.
readtype = {'r','r','i','i','i','m'};

% Signal length (valid for input and make).
lenlist = [512 8192 1024 1024 1024 1024];

% Representation dictionary.
qmf = MakeONFilter('Symmlet', 8);
dictlist  = {MakeList('WP'),MakeList('CP'),MakeList('WP'),MakeList('CP'),MakeList('WP'),MakeList('WP')};
pars1list = {MakeList(log2(lenlist(1))),MakeList(log2(lenlist(2))-4),MakeList(log2(lenlist(3))),MakeList(log2(lenlist(4))-4), ... 	
	     MakeList(log2(lenlist(5))-4),MakeList(log2(lenlist(6))-4)};
pars2list = {MakeList(qmf),MakeList(0),MakeList(qmf),MakeList(0),MakeList(qmf),MakeList(qmf)};
pars3list = {MakeList(0),MakeList(0),MakeList(0),MakeList(0),MakeList(0),MakeList(0)};

% Radius of the l1-ball for the Lasso. It is chosen equal to the l1 norm of the l2-constrained version.
%Rlist     = [76 5.4E5 8 360 150 330];

% Regularization parameter for the penalized version of the Lasso. 
% It is chosen such that the residual norm is <= l2-constraint radius (BPDN/Lasso l2-constrained form).
lambdalist     = [0.33 130 0.12 0.7 0.161 0.161];

% Titles: 'book' (subfigures), or 'toolbox' (meaningful titles).
titles = 'book';

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

for is=1:length(signals)
  n = lenlist(is);
  switch readtype{is}
  	case 'r'
  	 y0 = ReadSignal(signals{is});
	case 'i'
	 y0 = InputSignal(signals{is},n);
	case 'm'
	 y0 = MakeSignal(signals{is},n);
  end
  y0 = y0(:);

  dict  = dictlist{is};
  pars1 = pars1list{is};
  pars2 = pars2list{is};
  pars3 = pars3list{is};

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

  % Observation vector.
  PSNR = 25;
  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 	= 1/tightFrame; 	% Descent step-size for FB.
  options.tau	= 1;			% No relaxation -> 1.
  options.method= 'fb';
  options.niter	= 200;
  options.verb 	= 0;
  %R = Rlist(is);
  lambda = lambdalist(is);
  K  = @(x)FastS1(x,n,dict,pars1,pars2,pars3);
  KS = @(x)FastA1(x,dict,pars1,pars2,pars3);
  %ProxF = @(x,mu)compute_l1ball_projection(x,R);
  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, 1/tightFrame, options);timelassofb=toc
  options.tau	= (4-options.mu*tightFrame)/2;	% Relaxed FB.
  tic;[xhatlassofbrlx,Efbrlx] = fb(zeros(p,1), ProxF, GradG, 1/tightFrame, options);timelassofb=toc
  options.method = 'fista';
  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	= 200;
  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	= 200;
  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 = 200;
  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

  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-DR','BPDN-DR','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(xhatBPDR))),	     	...
							 	length(find(abs(xhatBPDNDR))),	     	...
							 	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(xhatBPDR,1), 	  ...
								norm(xhatBPDNDR,1), 	  ...
							 	norm(xhatDSPD,1));
						       
  fprintf('%-10s%-15g\t%-15g\t%-15g\t%-15g\t%-15g\n','CPU (s)', timelassofb, timelassofista, timeBPDR, timeBPDNDR, timeDSPD);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  figure;clf
  set(gcf,'Name',signals{is},'NumberTitle','off');
  subplot(3,1,1);
  plot(y0);axis tight
  xlabel('Time');
  if strcmp(titles,'book'),	h=title(['(a)']);set(h,'fontsize',18);
  else				title(['Signal:' signals{is}])
  end
  

  subplot(3,3,4);
  if strcmp(dict,'WP'), PhasePlane(xhatlassofb, 'WP', n, qmf, 256);
  else			PhasePlane(xhatlassofb, 'CP', n, 256);
  end
  if strcmp(titles,'book'),	h=title(['(b)']);set(h,'fontsize',18);
  else				title(['Lasso-ISTA'])
  end

  subplot(3,3,5);
  if strcmp(dict,'WP'), PhasePlane(xhatlassofista, 'WP', n, qmf, 256);
  else			PhasePlane(xhatlassofista, 'CP', n, 256);
  end
  if strcmp(titles,'book'),	h=title(['(c)']);set(h,'fontsize',18);
  else				title(['Lasso-FISTA'])
  end
  
  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
  if strcmp(titles,'book'),	h=title(['(d)']); set(h,'fontsize',18); end

  subplot(3,3,7);
  if strcmp(dict,'WP'), PhasePlane(xhatBPDR, 'WP', n, qmf, 256);
  else			PhasePlane(xhatBPDR, 'CP', n, 256);
  end			
  if strcmp(titles,'book'),	h=title(['(e)']);set(h,'fontsize',18);
  else				title(['BP-DR'])
  end

  subplot(3,3,8);
  if strcmp(dict,'WP'), PhasePlane(xhatBPDNDR, 'WP', n, qmf, 256);
  else			PhasePlane(xhatBPDNDR, 'CP', n, 256);
  end			
  if strcmp(titles,'book'),	h=title(['(f)']);set(h,'fontsize',18);
  else				title(['BPDN-DR'])
  end
  
  subplot(3,3,9);
  if strcmp(dict,'WP'), PhasePlane(xhatDSPD, 'WP', n, qmf, 256);
  else			PhasePlane(xhatDSPD, 'CP', n, 256);
  end			
  if strcmp(titles,'book'),	h=title(['(g)']);set(h,'fontsize',18);
  else				title(['DS-PD'])
  end

  colormap(1-gray)
  drawnow

  if savfig,	
  	saveas(gcf,['1D/Datasets/testsSparseApproxSynthesis1D' signals{is} '.fig'],'fig');
	saveas(gcf,['1D/Datasets/testsSparseApproxSynthesis1D' signals{is} '.eps'],'epsc');
  end
end

tilefigs