clear all;

rand('seed',2000);
randn('seed',2021);

alpha = 1;     % regularization parameter for MBA_{l1/l2}
maxiter = inf;
tol = 1e-6;

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

repeat = 20;

a = clock;
fname = ['Results\Lorentzian_table' '-'  date '-' int2str(a(4)) '-' int2str(a(5)) '.txt'];
fid = fopen(fname , 'w');
fprintf(fid, '%6s & %6s  & %6s  & %6s (%6s) & %6s & %6s & %6s & %6s \n', ...
    'index', 't_qr',  't_scp', 't_mba', 't_sum', 'err_scp' ,  'err_mba', 'res_scp', 'res_mba');

for ii = 1: length(indexarray)
    index = indexarray(ii);  % problem size
    n = 2560 * index;
    m = 720 * index;
    k = 80 * index;     % sparsity of xorig
    
    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,n);
        for j = 1:n
            A(:,j) = A(:,j)/norm(A(:,j));
        end
        error = nf*tan(pi*(rand(m,1) - 1/2));
        b = A*xorig + error;
        gamma = 0.02;
        tmpvalue =  sum(log(1 + error.^2/gamma^2));
        sigma = 1.2 * tmpvalue;         
        if sigma >= sum(log(1 + b.^2/gamma^2))
            error('0 is included in the feasible set. \n');
        end
        
        tstart0  = tic;
        [Q,R] = qr(A',0);
        time_qr = toc(tstart0);
        
        tstart1 = tic;
        xfeas1 = Q*(R'\b);       % taking A^{dagger}b as the initial point for SCP_{ls}_{LL2}
        time_xf1= toc(tstart1);
        fprintf('time_qr = %6.4f ; time_xf1 = %6.4f \n', time_qr, time_xf1);
        
        fprintf('Start of SCP_ls_LL2 \n');
        tstart3 = tic;
        [ x_scp, iter_scp, talg] = SCP_ls_LL2( A,b,sigma,1,xfeas1,0,gamma,tol,maxiter );
        time_scp = toc(tstart3);
        rec_err_scp = norm(x_scp - xorig)/max(1, norm(xorig));
        residual_scp = sum(log(1 + (A*x_scp - b).^2/gamma^2)) - sigma;
        fprintf(' SCP_ls_LL2 terminated for l1 : time = %4.1f, iter = %d,  nnz = %d,  fval = %g, rec_err = %g, residual = %7.4e  \n',...
            time_scp, iter_scp, nnz(abs(x_scp) > 1e-6), norm(x_scp,1), rec_err_scp, residual_scp);
        
        xfeas2 = x_scp;     % using the solution of SCP_{ls}_{LL2} as the initial point of MBA_{l1/l2}
        
        fprintf('Start of MBA_l1vl2 \n');
        tstart = tic;
        [xsoln, iter_mba] = MBA_l1vl2_LL(A, b, gamma, sigma, xfeas2, alpha, maxiter, tol);
        time_mba = toc(tstart);
        rec_err_mba = norm(xsoln - xorig)/max(1, norm(xorig));
        residual_mba =  sum(log(1 + (A*xsoln - b).^2/gamma^2)) - sigma;
        time_sum = sum( [time_qr  time_xf1  time_scp  time_mba] );
        fprintf('  MBA terminated:  time_mb = %4.1f, iter = %d, nnz %d,  fval = %g, rec_err = %g, residual = %7.4e \n',...
            time_mba,iter_mba, nnz(abs(xsoln) > 1e-6), norm(xsoln, 1)/norm(xsoln), rec_err_mba, residual_mba);
        
        process = [process; time_qr,time_scp, time_mba, time_sum, rec_err_scp, rec_err_mba, residual_scp, residual_mba];
    end
    
    fprintf(fid, '%6d & %6.1f   & %6.1f &  %6.1f (%6.1f) & %6.1e & %6.1e & %6.0e& %6.0e \n' , index, mean(process) );
end
fclose(fid);




