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

% Using MBA_{l1/l2} to solve the robust compressed sensing problem
% min \|x\|_1 / \|x\|_2  s.t. dist^2(Ax-b, S) <= \sigma^2 ,  where S = {u: \|u\|_0 <= sp_outlier}

% Inputs
%   A                - m times n matrix; the sensing matrix
%   b                - m dimensional vector; the measurement
%   sp_outlier       - the no. of outliers; positive interger 
%   sigma            - the  parameter for control the residual; satisfying dist(b, S) > sigma >0
%   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;
[~, I] = sort(abs(u), 'descend');
I = I(sp_outlier+1:end);
uI = u(I);
qval = norm(uI)^2 - sigma^2;          
v = zeros( length(b) ,1);
v(I) = uI;        % v = Ax - b - Proj_S(Ax-b)
subgq =  A' * (2 * v);                  % subgradient of q at x
normsubgq = norm(subgq);
s = x - subgq/L;
r = normsubgq^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;
        [~, Itilde] = sort(abs(utilde), 'descend');
        Itilde = Itilde(sp_outlier+1:end);
        uItilde = utilde(Itilde);
        qvaltilde = norm(uItilde)^2 - sigma^2;
        
        if qvaltilde <= 0     
            break;
        else
            L = 2 * L;
            s = x - subgq/L;
            r = normsubgq^2/L^2 - 2*qval/L;
        end
    end
    
    normxnew = norm(xtilde);
    norm1xnew = norm(xtilde, 1);
    v = zeros( length(b),1);
    v(Itilde) = uItilde;
    subgqtilde =  A' *(2* v);
    dx = xtilde - x;
    normdx = norm(dx);
    dg = subgqtilde - subgq;

    % termination criterion
    if  normdx <= tol * max(normxnew, 1)
        x = xtilde;
        fprintf(' Termination of MBA_l1vl2 for RCS: 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;
    subgq = subgqtilde;
    normsubgq = norm(subgq);
    s = x - subgq/L;
    r = normsubgq^2/(L^2) - 2*qval/L;
    
    iter = iter +1;
end