function [x,iter,fnew,time] = GIST(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
% epsilon - real number > 0
% opts    - contains the following structures variables, default values are in brackets:
%           x0      - starting point [0]
%           tol     - tolerance [1e-4]
%           maxiter - maximum number of iterations [inf]
%           disp    - display if equals 1 [0]
%           M       - memory of nonmonotone line search [2]
% 
% Output
%
% x       - approximate stationary point
% iter    - number of iterations
% fnew    - final function value
% time    - cputime

tstart = tic;

%Parameters
n = size(A,2);

if nargin <= 4
    x0 = zeros(n,1);
    tol = 1e-4;
    maxiter = inf;
    M = 5;
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,'maxiter'), maxiter = opts.maxiter; else maxiter = inf;  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


step_init = 1; % initial stepsize 

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

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 = A*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 = A*xnew - b;
        fnew = 0.5*norm(Atmp)^2 + mu*sum(log(1+abs(xnew)/epsilon));
        count = count + 1;
    end
    gnew = A'*Atmp;
    y = gnew - g1;
    s = xnew - x;
    
    % check for termination
    if (norm(y) + norm(s)/stepsize)/max(1,norm(xnew)) < tol
        x = xnew;
        break
    end
    
    % Barzilai-Borwein stepsize        
    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;
    
end

time = toc(tstart);