function [ z,iter,talg ] = IRL1e2( 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 IRL1e2
%
% 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 - 2*refreq is the cycle length for repeating theta
% 
% Output
%
% z       - approximate stationary point
% iter    - number of iterations
% talg    - cputime

tstart_alg = tic;
[~,n] = size(A);
x = zeros(n,1);
z = x;
iter = 1;
% generate theta for use
theta = zeros(2*re_freq,1); % generate a theta vector
theta(1) = 1;
count = 1;
while count < 2*re_freq
    if  count <= re_freq - 1
        theta(count+1) = 2/(1 + sqrt(1 + 4/theta(count)^2));
    elseif count >= re_freq + 1
        theta(count+1) = theta(count)/sqrt(1 - theta(count)); % inverse of the formula above
    else
        theta(count+1) = theta(count);
    end
    count = count + 1;
end

ell = lambda/epsilon^2;

%updates in IRL1e2
while iter < maxiter
    %update s
    s = lambda./(abs(x)+epsilon);
    
    %update y
    k = mod(iter-1,2*re_freq)+1;%counter for theta
    y = (1-theta(k))*x + theta(k)*z; 
    
    % update x
    G = A'*(A*y-b);
    tmp = z - G/theta(k)/L;
    z = sign(tmp).*max(abs(tmp) - s/theta(k)/L,0);
    
    %update x
    xold = x;
    x = (1-theta(k))*x + theta(k)*z;
 
    
    %check for termination
    if (L*norm(z - y) + L*norm(x - y) + ell*norm(z - xold))/max(1,norm(z)) < tol 
        break
    end
    
    iter = iter + 1;    
end

talg = toc(tstart_alg);
