clear all;

rand('seed',1999);
randn('seed',2020);

alpha = 1;     % regularization parameter for MBA
maxiter = inf;
tol = 1e-6;

indexarray = [2, 4, 6, 8, 10];      % problem size
nf = 0.01;

repeat = 20;

a = clock;
fname = ['Results\RCS_table' '-' date '-' int2str(a(4)) '-' int2str(a(5)) '.txt'];
fid = fopen(fname, 'w');
fprintf(fid, '%6s & %6s & %6s (%6s) & %6s & %6s \n' , 'index', 't_qr', 't_mba', 't_sum', 'rec_err', 'residual');

for ii = 1:length(indexarray)
    index = indexarray(ii);
    n = 2560*index;
    m = 720*index;
    k = 80*index;
    r = 10*index;
    
    process = [];
    
    for rr = 1 : repeat        
        % generate the original signal
        I = randperm(n);
        J = I(1:k);
        xorig = zeros(n,1);
        xorig(J) = randn(k,1);
        
        A = randn(m+r,n);
        for j = 1:n
            A(:,j) = A(:,j)/norm(A(:,j));
        end
        z = [zeros(m,1); 2*sign(randn(r,1)) ];      
        error = nf*randn(m+r,1);
        
        b = A*xorig - z + error;
        sigma = 1.2*norm(error);       
        sp_outlier = 2*r;  
        [~, bI] = sort(abs(b), 'descend');
        bI = bI(sp_outlier+1:end);
        distbS = norm(b(bI));
        if distbS <= sigma
            error('0 is included in the feasible set. \n');
        end
        
        tstart = tic;
        [Q,R] = qr(A',0);
        time_qr = toc(tstart);
        
        tstart1 = tic;
        xfeas = Q*(R'\b);        % taking A^{dagger}b as the initial point
        time_xf = toc(tstart1);
        fprintf('time_qr  =  %6.4f, time_xf  = %6.4f \n',  time_qr, time_xf);
        
        tstart3 = tic;
        [xsoln, iter_mb] = MBA_l1vl2_RCS(A, b, sp_outlier, sigma, xfeas, alpha, maxiter, tol);
        time_mb = toc(tstart3);
        
        rec_err = norm(xsoln - xorig)/max(1, norm(xorig));
        u = A*xsoln - b;
        [v, I] = sort(abs(u), 'descend');
        I = I(sp_outlier+1:end);
        residual = norm(u(I))^2 - sigma^2;
        time_sum = sum( [time_qr, time_xf,  time_mb] );
        process = [process; time_qr,  time_mb, time_sum, rec_err, residual];
        fprintf(' MBA terminated: time_mb = %4.2f, nnz = %d,  fval = %7.4e, rec_err = %g, residual = %7.4e \n',...
            time_mb, nnz(abs(xsoln) > 1e-6), norm(xsoln, 1)/norm(xsoln), rec_err, residual);
    end
        fprintf(fid, '%6d & %6.1f & %6.1f (%6.1f)  & %6.1e & %6.0e \n', index, mean(process));
    
end
fclose(fid);

