function [x,iter,time] = NPG_ls_log(A,b,mu,epsilon,opts)

% This function computes a stationary point for the following regularized l1-l2 model
% 
% min 0.5*||Ax - b||^2 + mu*sum(log(1+abs(x)/epsilon))
% 
% It uses NPG
%
% Input
%
% A       - m by n matrix (m << n)
% b       - m by 1 vector measurement
% mu      - real number > 0
% gamma   - real number between 0 and 1
% opts    - contains the following structures variables, default values are in brackets:
%           x0      - starting point [0]
%           tol     - tolerance [1e-4]
%           freq    - frequency for displaying output [5000]
%           maxiter - maximum number of iterations [inf]
%           L       - user input (upper) estimate of norm of A'*A [norm(A'*A)]
%           k       - as explained in the problem [1]
%           disp    - display [1]: 1 - display per freq iterations
%                                  0 - no display
% 
% Output
%
% x       - approximate stationary point
% iter    - number of iterations
% time    - cputime

tstart = tic;

% Parameters
n = size(A,2);

if nargin <= 4
    x0 = zeros(n,1);
    tol = 1e-4;
    freq = 5000;
    maxiter = inf;
    disp = 1;
else
    if isfield(opts,'x0'),      x0      = opts.x0;      else x0 = zeros(n,1);end
    if isfield(opts,'tol'),     tol     = opts.tol;     else tol = 1e-4;     end
    if isfield(opts,'freq'),    freq    = opts.freq;    else freq = 5000;    end
    if isfield(opts,'maxiter'), maxiter = opts.maxiter; else maxiter = inf;  end
    if isfield(opts,'disp'),    disp    = opts.disp;    else disp = 0;       end
    if isfield(opts,'M'),       M       = opts.M;       else M = 2; end
end

% initialization
x = x0;
eta = 2; % line search parameter
sigma = 1e-4; % line search parameter

% mg = mu*gamma;
step_init = 1; % initial stepsize 

% set matrix handles
Ahandle = @(x) A*x;
AThandle = @(x) A'*x;

% compute function value
fnew = 0.5*norm(Ahandle(x) - b)^2 + mu*sum(log(1+abs(x)/epsilon));
% compute gradient
g1 = AThandle(Ahandle(x) - b);

if disp == 1
  fprintf(' ********** start **********\n')
end

iter = 1;
f1 = -inf(M,1); % set_up for nonmonotone line search
f1(1) = fnew;
fmax = max(f1);

while iter < maxiter
    
    % line search
    stepsize = step_init;
    
    xnew = ls_log_sub(x - stepsize*g1,stepsize,epsilon,mu,n);
    Atmp = Ahandle(xnew) - b;
    fnew = 0.5*norm(Atmp)^2 + mu*sum(log(1+abs(xnew)/epsilon));
    count = 0;
    while fnew > fmax - sigma/2*norm(xnew - x)^2
        stepsize = stepsize/eta;
        xnew = ls_log_sub(x - stepsize*g1,stepsize,epsilon,mu,n);
        Atmp = Ahandle(xnew) - b;
        fnew = 0.5*norm(Atmp)^2 + mu*sum(log(1+abs(xnew)/epsilon));
        count = count + 1;
    end
    
    % check for termination
    if norm(xnew - x)/max(1,norm(xnew)) < tol
        if disp == 1
            fprintf(' TERMINATION \n')
            fprintf(' iter = %5d   fval = %12.6f\n', iter, fnew)
        end
        x = xnew;
        break
    end
    
    % Barzilai-Borwein stepsize
    s = xnew - x;
    gnew = AThandle(Atmp);
    y = gnew - g1;
    step_init = max(min(norm(s)^2/(s'*y),1e8),1e-8);
    
    % update iterate, gradient and function value
    x = xnew;
    g1 = gnew;
    f1(mod(iter,M)+1,1) = fnew;
    fmax = max(f1);
    iter = iter + 1;
    
    if mod(iter,freq) == 0 && disp == 1
        fprintf(' iter = %5d   fval = %12.6f  ls_used = %3d\n', iter, fnew, count)
    end
    
end

time = toc(tstart);