function [ x,iter,talg ] = IRL1e1( A,b,epsilon,lambda, L, tol,maxiter, re_freq)
% This function computes a stationary point for the following model
% 
% min 0.5*||Ax - b||^2 + lambda*sum((abs(x)/epsilon + 1))
% 
% It uses IRL1e1
%
% Input
%
% A       - m by n matrix (m << n)
% b       - m by 1 vector measurement
% lambda  - real number > 0
% epsilon - real number > 0
% L       - user supply value of max(eig(A*A'))
% tol     - tolerance
% maxiter - max number of iteration
% re_freq - restart frequency for fixed restart
% 
% Output
%
% x       - approximate stationary point
% iter    - number of iterations
% talg    - cputime

tstart_alg = tic;
[~,n] = size(A);
xold = zeros(n,1);
x = xold;
iter = 1;
theta = 1;
thetaold = 1;
num_restart = 0;
ell = lambda/epsilon^2;


%updates in IRL1e1
while iter < maxiter
    %update s
    s = lambda./(abs(x)+epsilon);
    
    %update beta and y
    beta = theta*(1/thetaold - 1);
    y = x + beta*(x - xold);
    xold = x;
    
    % update x
    G = A'*(A*y-b);
    tmp = y - G/L;
    x = sign(tmp).*max(abs(tmp) - s/L,0);
     
    if (2*L*norm(x - y) + ell*norm(x - xold))/max(1,norm(x)) < tol
        break
    end
    
    %update theta
    thetaold = theta;
    theta = 2/(1 + sqrt(1 + 4/theta^2));
    
    %restart theta
    if re_freq < inf
        if iter > 0 && (mod(iter, re_freq) == 0 || (y-x)'*(x-xold) > 0)
            num_restart = num_restart + 1;
            theta = 1;
            thetaold = 1;
        end
    end
    
    iter = iter + 1;
end


talg = toc(tstart_alg);



