clear all;

rand('seed',2001);
randn('seed', 2022);

addpath('spgl1-2.1');
optsL1 = spgSetParms('verbosity',0);

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

nf = 0.01;
n = 1024;
m =  64;  % number of observations
sarray = [8 12];
Farray = [5 15];
Darray = [2  3];

repeat = 20;

a = clock;
fname = ['Results\LS-Lou-table' '-' date '-' int2str(a(4)) '-' int2str(a(5)) '.txt'];
fid = fopen(fname , 'w');
fprintf(fid, '%3s & %3s &%3s & %6s & %6s (%6s)  & %6s & %6s & %6s & %6s \n',...
    's', 'F', 'D',  't_spgl1', 't_mba', 't_sum',  'err_spgl1', 'err_mba', 'res_spgl1', 'res_mba');

for ss = 1:length(sarray)
    s = sarray(ss);
    for ff = 1:length(Farray)
        F = Farray(ff);
        for dd = 1:length(Darray)
            D = Darray(dd);
            
            process = [];
            
            for rr = 1 : repeat                
                % generate the original signal
                I = randperm(n);
                J = I(1:s);
                xorig = zeros(n,1);
                xorig(J) = sign(randn(s, 1)) .* 10.^(D*rand(s, 1));
                
                A = zeros(m,n);
                w = rand(m,1);
                for j = 1:n
                    A(:,j) = cos( 2*pi*j*w/F ) / sqrt(m);
                end
                
                error = nf*randn(m,1);
                b = A*xorig + error;
                sigma = 1.2*norm(error);
                if sigma >= norm(b)
                    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);            
                time_xf = toc(tstart1);
                fprintf('time_qr = %6.4f,  time_xf = %6.4f\n', time_qr, time_xf);
                
                fprintf('Start of SPGL1 \n');
                tstart2 = tic;
                x_spgl1 = spg_bpdn(A,b,sigma,optsL1);
                time_spgl1 = toc(tstart2);
                
                rec_err_spgl1 =  norm(x_spgl1 - xorig)/max(1, norm(xorig));
                residual_spgl1 = norm(A*x_spgl1-b)^2 - sigma^2;
                fprintf(' SPGL1 terminated for l1 : time = %4.1f, nnz = %d,  fval = %7.4e, rec_err = %g, residual = %g \n',...
                    time_spgl1, nnz(abs(x_spgl1) > 1e-6), norm(x_spgl1,1), rec_err_spgl1, residual_spgl1);
                
                normAxb = norm(A*x_spgl1 - b);
                % backtrack to feasibility
                if normAxb <= sigma % using the solution of SPGL1 (with slight backtracking) as the initial point of MBA_{l1/l2}
                    xfeas2 = x_spgl1;
                else
                    xfeas2 = xfeas1 + sigma*(x_spgl1 - xfeas1)/normAxb;
                end
                
                fprintf('Start of MBA_l1vl2_LS \n');
                tstart3 = tic;
                [xsoln, iter_mba]  = MBA_l1vl2_LS(A, b, sigma, xfeas2, alpha, maxiter, tol);
                time_mb = toc(tstart3);
                
                rec_err_mba = norm(xsoln - xorig)/max(1, norm(xorig));
                residual_mba = norm(A*xsoln-b)^2 - sigma^2;
                time_sum = sum( [time_qr, time_xf, time_spgl1, time_mb] );
                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_mba, residual_mba);
                
                process = [process;  time_spgl1, time_mb, time_sum, rec_err_spgl1,  rec_err_mba,...
                   residual_spgl1, residual_mba];

            end
            
            fprintf(fid,  '%3d & %3d & %3d & %6.2f & %6.2f (%6.2f)  & %6.1e  & %6.1e & %6.0e & %6.0e  \n', s ,F, D, mean(process));
            
        end
    end
end
fclose(fid);




