function [x,iter,Fval,talg] = ReNPG(A,b,epsilon,lambda,tol,maxiter,M)
% This function computes a stationary point for the following model
% 
% min 0.5*||Ax - b||^2 + lambda*sum((abs(x)/epsilon + 1))
% 
% It uses Reweighted l1 with nonmonotone line search
%
% Input
%
% A       - m by n matrix (m << n)
% b       - m by 1 vector measurement
% lambda  - real number > 0
% epsilon - real number > 0
% tol     - tolerance
% maxiter - max number of iteration
% M       - nonmonotone linesearch memory
% 
% Output
%
% x       - approximate stationary point
% iter    - number of iterations
% Fval    - final function value
% talg    - cputime


tstart_alg = tic;
[~,n] = size(A);
x = zeros(n,1);
iter = 1;
L = 1;
c = 1e-4;
ell = lambda/epsilon^2;

tmp0 = A*x - b;
fval = 0.5*norm(tmp0)^2 + lambda*sum(log(abs(x)/epsilon + 1));
G = A'*tmp0;

f1 = -inf(M,1);
f1(1) = fval;
fmax = max(f1);

while iter < maxiter
    % update s
    s = lambda./(abs(x)+epsilon);
    
    % update x
    tmp = x - G/L;
    xnew = sign(tmp).*max(abs(tmp) - s/L,0);  
    tmp0 = A*xnew - b;
    fnew = 0.5*norm(tmp0)^2 + lambda*sum(log(abs(xnew)/epsilon + 1));
    
    while fmax - fnew < c/2*norm(xnew-x)^2;
        L = 2*L;
        tmp = x - G/L;
        xnew = sign(tmp).*max(abs(tmp) - s/L,0);
        tmp0 = A*xnew - b;
        fnew = 0.5*norm(tmp0)^2 + lambda*sum(log(abs(xnew)/epsilon + 1));
    end
    Gnew = A'*tmp0;
    S = xnew - x;
    Y = Gnew - G;
    
    %check for termination
    if (norm(Y) + (L + ell)*norm(S))/max(1,norm(xnew)) < tol
        x = xnew;
        break
    end
    
    %update L
    L = max(1e-8,min(1e8,(Y'*S)/norm(S)^2));
    
    % Update iterates
    x = xnew;
    G = Gnew;
    f1(mod(iter,M)+1) = fnew;
    fmax = max(f1);
    

    iter = iter + 1;
end
Fval = fnew;
talg = toc(tstart_alg);

