function [x,iter,time,count] = FB_env_ls_l1_l2_lbfgs(A,b,mu,lambda,opts)

% This function computes a stationary point for the following regularized l1 model
% 
% min 0.5*||Ax - b||^2 + mu*( ||x||_1 - lambda*||x||_2)
% by solving the following 2n - dimensional problem
% min 0.5*||Az - b||^2 - mu_2*<y,z> + mu_1*||z||_1 +...
%     delta_{||.|| <= 1}(y) 
% It uses the algorithmic framework as in L-BFGS and calls the subroutine
% soft_thresh and proj for solving the subproblem. It minimizes the FB envelope.
%
% Input
%
% A       - m by n matrix (m << n)
% b       - m by 1 vector measurement
% mu      - real number > 0
% lambda  - real number in  (0,1]
% opts    - contains the following structure variables, default values are in brackets:
%           x0      - starting point [0]
%           tol     - tolerance [1e-4]
%           maxiter - maximum number of iterations [inf]
%           L       - user input (upper) estimate of norm of A'*A [norm(A'*A)]
%           M       - number of former vectors used for L-BFGS
%           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);
step_init = 1;
eta = 2;   % parameter for Armijo line search 
c1 = 1e-4; % parameter for Armijo line search 
c2 = 1e5;


if nargin <= 4
  x0 = zeros(2*n,1);
  tol = 1e-6;
  maxiter = inf;
  L = 1; % does not require computation of norm(A'*A)
  M = 10; % memory of L-BFGS
  disp = 1;
  freq = 50;

else
  if isfield(opts,'x0'),      x0      = opts.x0;      else x0 = zeros(2*n,1);end
  if isfield(opts,'tol'),     tol     = opts.tol;     else tol = 1e-6;     end
  if isfield(opts,'maxiter'), maxiter = opts.maxiter; else maxiter = inf;  end
  if isfield(opts,'L'),       L       = opts.L;       else L = 1;          end
  if isfield(opts,'M'),       M       = opts.M;       else M = 10;         end
  if isfield(opts,'disp'),    disp    = opts.disp;    else disp = 1;       end
  if isfield(opts,'freq'),    freq    = opts.freq;    else freq = 50;       end
end

% Initialization
x = x0;
gamma = 0.95/L; % different from 1/eig value of A'*A
mu2 = mu*lambda;
% set matrix handles

Ahandle = @(z) A*z;
AThandle = @(z) A'*z;


% Compute function value and gradient
[fval,G_gamma] = fvalue(x,b,n,Ahandle,AThandle,gamma,mu,mu2);
grad = gvalue(n,G_gamma, Ahandle, AThandle,mu2,gamma);

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


% Define vectors for L-BFGS
S = zeros( 2*n,M );
Y = zeros( 2*n,M );

% Initialization
z = x;
g = grad;
seq = 1;
iter = 1;
count = 0;

while iter <= maxiter
    
    % L-BFGS direction    
    t = min( iter-1, M);
    if iter == 1
        p = -grad;
    else
        p = LBFGS_two_loop( t, grad, S, Y ,seq);
    end
    dderiv = grad'*p;
    % Check the angle condition
    ng = norm(grad);
    np = norm(p);
    if (dderiv > -1e-5*ng*np)||(ng/np > c2)||(ng/np < 1/c2)
        p = -grad;
        dderiv = grad'*p; % added on May 19, 2016; bug fixed. Only change after submission of paper in May 2016
        fprintf(' angle condition failed, steepest descent used, iter = %d\n', iter)
        count = count + 1;
    end
    
    stepsize = step_init;
    
    % Armijo line search
    xnew = x + stepsize*p;
    [fnew,G_gamma] = fvalue(xnew,b,n,Ahandle,AThandle,gamma,mu,mu2);
    while fnew > fval + c1*stepsize*dderiv
        stepsize = stepsize/eta;
        xnew = x + stepsize*p;
        [fnew,G_gamma] = fvalue(xnew,b,n,Ahandle,AThandle,gamma,mu,mu2);
    end
    gnew = gvalue(n,G_gamma, Ahandle, AThandle,mu2,gamma);
    if  mod(iter,freq) == 0 && disp == 1 
        fprintf(' iter = %5d nnz = %5d stepsize = %2.1e fval = %12.8f normg = %12.8f derec = %12.8f  \n',...
            iter, nnz(x(1:n)), stepsize, 0.5*norm(Ahandle(x(1:n)) - b)^2 + mu*norm(x(1:n),1) - mu2*norm(x(1:n),2),norm(grad), grad'*p)
    end
    
    if stepsize < 1e-10
        fprintf(' stepsize too small \n')
    end
    
    % check for termination
    if norm(gnew)/max(1,fnew) < tol
        if disp == 1
            fprintf(' TERMINATION \n')
            fprintf(' iter = %5d   fval = %12.6f  normg = %12.6f \n', iter, fval, norm(grad))
        end
        x = xnew;
        break 
    end
    
    % Update sequences for L-BFGS
    [z,g,S,Y,seq] = sequence_update(t,M,z,g,S,Y,xnew,gnew,seq);

    % Update iterations 
    x = xnew;
    fval = fnew;
    grad = gnew;
    iter = iter + 1;   
end
time = toc(tstart);

    
    
    
    
    
    
    
    
    
    
    
    
  