clear all
global dict Omega
MC=10;
j=[6:14];
maxIters=50;

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

% Parameters for DR.
options.gamma = 0.1;    % Stepsize parameter for Douglas-Rachford iteration.
options.tau   = 1;    % No relaxation.
options.niter = maxIters;
options.verb  = 0;
ProxF = @(x,ga)SoftThresh(x,ga);

OptTol  = 1E-6;

for ij=1:length(j)
  n=2^j(ij);  		% Number of measurements.
  p=n*4;      		% Coefficient vector length.
  k=ceil(0.05*n);  	% Sparsity level.

for rep=1:MC

  progressbar((ij-1)*MC+rep,MC*length(j));
  
  x = SparseVector(p, k, 'GAUSSIAN', true);

  % Random measurement (sensing) operator: Hadamard, Fourier, Real Fourier, Real sinusoid (RST), USE, etc.
  dict = 'Hadamard';
  tightFrame = 1;  % e.g. for Hadamard, Fourier and RST (tight frames), (Phi Phi') = I_n. Otherwise, if unknown or frame, set to 0.
  q = randperm(p);
  Omega = q(1:n)';

  % Observed data.
  y = FastMeasure(x, dict, Omega);

  K  = @(x)FastMeasure(x, dict, Omega);
  KS = @(x)FastMeasureAdjoint(x, p, dict, Omega);
  ProxG = @(x,ga)x + KS(y - K(x))/tightFrame;

  tic;xhatBPDR = dr(zeros(p,1),ProxF,ProxG,options);timeBPDR(rep,ij)=toc;
  tic;xhatlars = real(SolveLasso('FastCSOp', y, p, 'lars', maxIters, 0, 0, 0, OptTol));timelars(rep,ij)=toc;
  tic;xhatLP = real(SolveBP('FastCSOp', y, p, maxIters, 0, OptTol));timeLP(rep,ij)=toc;
  tic;xhatStOMP = real(SolveStOMP('FastCSOp', y, p, 'FAR', 0.01, maxIters, 0, OptTol));timeStOMP(rep,ij)=toc;
  errBPDR(rep,ij) = norm(xhatBPDR-x);
  errlars(rep,ij) = norm(xhatlars-x);
  errLP(rep,ij)   = norm(xhatLP-x);
  errStOMP(rep,ij)= norm(xhatStOMP-x);
  
end
end

[mean(errBPDR,1)' mean(errlars,1)' mean(errLP,1)' mean(errStOMP,1)']

figure
semilogy(j,mean(timeBPDR,1),'-b',j,mean(timelars,1),'--b',j,mean(timeLP,1),'-.b',j,mean(timeStOMP,1),'.-b');axis tight
set(gca,'FontSize',12);
legend('BP-DR','LARS','LP-Interior Point','StOMP',0);
xlabel('log_2(m)');ylabel('Computation time (s)');
title(sprintf('\\Phi=Dirac, H=%s, m/n=%d%%, sparsity=%.f%%',dict,100*n/p,100*k/n));

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