clear all

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

n = 100;
s = [1:1:n-1];	% Sparsity.
m = [1:1:n-1];	% Number of measurements.
MC = 50;
PT = zeros(length(m),length(s));

for is=1:length(s)
  for im=1:length(m)
    for k=1:MC
      x0 = SparseVector(n,s(is),'Signs',1);
      A = randn(m(im),n);
      y = A*x0;
      x = SolveBP(A, y, n, 300, 0, 1E-6);
      %plot([x0 x]);title(sprintf('s=%d m=%d',s(is),m(im)));pause
      PT(im,is) = PT(im,is) + (norm(x-x0) <= 1E-5);
    end
   progressbar((is-1)*length(m)+im,length(s)*length(m));
   end
end
PT = PT/MC;


figure
imagesc(s,m,PT);axis xy;
colormap('gray')
hold on


% Theoretical phase transitions through the statistical dimension.
rho = @(x,t,a)x*(1+t^2)+((1-x)*((1+t^2)*erfc(t/sqrt(2))-sqrt(2/pi)*t*exp(-t^2/2)))/a;
drho= @(x,t,a)sqrt(2/pi)*exp(-t^2/2)/t - erfc(t/sqrt(2)) - a*x/(1-x);

r = linspace(1E-3,1-1E-3,100); % Relative sparsity r.
for ir=1:length(r)
 lB = erfcinv(min(sqrt(2)/(1-r(ir)),1-eps(1)));
 uB = sqrt(2/pi)*(1/(r(ir))-1); 
 ml1(ir)  = rho(r(ir),fzero(@(t)drho(r(ir),t,1),[lB uB]),1);
end

plot(r*n,n*ml1,'-',r*n,n*(2*r.*log(1./r)+5/4*r),'--b');axis([0.5 n-0.5 0.5 n-0.5])
xlabel('Sparsity');
ylabel('Number of measurements');
h = legend('Statistical dimension','$2s\log(n/s)+5/4s$');
set(h,'Interpreter','latex','Location','NorthWest');
legend boxoff;
hold off

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