tic
N=1500;	% length of the signal
k=45;	% number of nonzero terms in the signal
m=250;	% number of measurements
K=k+10;	% the K parameter of the adapted DDFG algorithm 

eps0=1;
Phi = randn(m,N); % draw a Gaussian measurement matrix
Phi = Phi./(ones(m,1)*sqrt(sum(Phi.^2))); % normalize each column

I = randperm(N); 
T = I(1:k); % generate k random locations for the support T of x
xr = randn(k,1); % random vector to be placed on T
x = zeros(N,1);
x(T) = xr;

y = Phi*x; % measured data

N_it = 49; % number of iterations to run the algorithm
X = zeros(N,N_it+31); % would like to store all the iterates
Y = zeros(N,N_it+31);
X(:,1) = 1; % this will imply that we start with uniform weights

Eps = eps0*ones(1,N_it+31); % initialize the epsilon sequence 

for l=1:N_it+30;
    X(:,l+1) = it_min_l2_w_fast(Phi,y,X(:,l),Eps(l));
    Y(:,l+1) = it_min_l2_w_fast(Phi,y,X(:,l),(1e-7));
    figure(1);
    plot(1:N,x,'rx',1:N,X(:,l+1),'bo');
    pause(0.001);
    sX = sort(abs(X(:,l+1)),'descend'); 
    Eps(l+1) = max(min(Eps(l),sX(K)/N),(1e-7)/N); % 1e-7 is arbitrary
end

figure(2);
subplot(121);
E = sum(abs(X(:,[1:N_it+1])-x*ones(1,N_it+1))); % convergence in l_1 to x
semilogy(E,'+-b');
hold on
%xf = Y(:,N_it+31); % exit state of the vector after N_it iterations
%E2 = sum(abs(Y(:,[1:N_it+1])-xf*ones(1,N_it+1))); % convergence in l_1 to the limit xf
%semilogy(E2,'.-b');
semilogy(Eps(1:N_it+1),'*-r');
%semilogy([1 N_it],[1 1]*N*Eps(end),'k'); 
%semilogy([1 N_it],[1 1]*min(abs(xr)),'g');
legend('||x^n-x||_1','\epsilon_n','Location','NorthEast');
title('decay of error and \epsilon_n');
hold off

subplot(122);
alpha = exp(diff(log(E(1:N_it))));
%alpha2 = exp(diff(log(E2(1:N_it))));
alpha3 = exp(diff(log(Eps(1:N_it))));
plot(1:N_it-1,alpha,'+-b',1:N_it-1,alpha3,'*-r');
legend('istantaneous rate for ||x^n-x||_1','istantaneous rate for \epsilon_n','Location','SouthEast');
title('instantaneous rate of convergence');
toc