function [x, iter] = MBA_l1vl2_LL(A, b, gamma, sigma, xfeas, alpha, maxiter, tol)

% Using MBA_{l1/l2} to solve the compressed sensing problem with Cauchy noise
% min \|x\|_1 / \|x\|_2  s.t.  \|Ax - b\|_{LL2,gamma} - delta<= 0, where {LL2, gamma} represents the Lorentzian norm

% Inputs
%   A                - m times n matrix; the sensing matrix
%   b                - m dimensional vector; the measurement
%   gamma            - the parameter in the Lorentzian norm; gamma > 0
%   sigma            - the  parameter for control the residual; 0 < sigma < \|b\|_{LL2, gamma}
%   xfeas            - the initial point (feasible)
%   alpha            - the parameter in the approximation of objective; alpha > 0
%   maxiter          - the maximal iteration number
%   tol              - the tolerance for termination 

% Outputs
%   x:               - the approximal solution given by MBA_{l1/l2}
%   iter:            - the total iterations 

x = xfeas;
L = 1;          % the initial L chosen in [ lmin, lmax ]
iter = 0;

normx = norm(x);
norm1x = norm(x,1);
omega = norm1x/normx;
y = (1 + omega/alpha/normx) * x;
Ax = A * x;
u = Ax - b;
qval = sum(log(1+u.^2/gamma^2)) - sigma;
grad = A' * ( 2*u./(gamma^2+u.^2) );
normgrad = norm(grad);
s = x - grad/L;
r = normgrad^2/(L^2) - 2*qval/L;

% start of the outer loop
while iter < maxiter
    % start of the inner loop
    while 1 == 1
        xtilde= SubP_alpha( y,s,r,alpha );
        Axtilde = A * xtilde;
        utilde = Axtilde - b;
        qvaltilde = sum(log(1+utilde.^2/gamma^2)) - sigma;
        
        if  qvaltilde <= 0    
            break;
        else
            L = 2 * L;
            s = x - grad/L;
            r = normgrad^2/L^2 - 2*qval/L;
        end
    end
    
    normxnew = norm(xtilde);
    norm1xnew = norm(xtilde, 1);
    gnew = A' * ( 2*utilde./(gamma^2+utilde.^2) );
    dx = xtilde - x;
    normdx = norm(dx);
    dg = gnew - grad;
    
    % termination criterion
    if  normdx <= tol * max(normxnew, 1)
        x = xtilde;
        fprintf(' Termination of MBA_l1vl2 for LL2 : iter=%d, fval=%g, residual=%6.4e,  sparsity = %g\n', iter, norm1xnew/normxnew, qvaltilde, nnz(abs(x) > 1e-6));
        break
    end
    
    % BB stepsize
    dgdx = dg'*dx;
    if dgdx > 1e-12
        L = max( 1e-8 , min(1e8 , (dgdx)/(normdx^2) ) );
    else
        L = max( 1e-8 , min(1e8 , L/2 ) );
    end
    
    % update   
    x = xtilde;
    normx = normxnew;
    norm1x = norm1xnew;
    omega = norm1x/normx;
    y = (1 + omega/alpha/normx) * x;
    qval = qvaltilde;
    grad = gnew;
    normgrad = norm(grad);
    s = x - grad/L;
    r = normgrad^2/(L^2) - 2*qval/L;
    
    iter = iter +1;
end


