function runcode_DIRL_bian

% Random instance
% Give the values of parameters needed

addpath('spgl1-2.1')
randn('seed',2020) % These fix the instance we test
rand('seed',2002)

repeats_index = 5;                        
table1_ADMM = []; table1_SPGL1 = []; table1_SCP = []; table1_L_time = [];

for j = 1:repeats_index
    
 repeats = 30;                             
 table2_ADMM = []; table2_SPGL1 = []; table2_SCP = []; table2_L_time = [];


for repeat = 1:repeats
   
    
    index = 2*j;                           
    n = 2560*index;
    m = 540*index; % originally 720*index
    s = 80*index; % nonzero entries
    sigma0 = 0.01;
    delta = 5*sigma0; % parameter for Lorentzian norm
    eps = 0.1; % parameter for logpen
    tol_min = 1e-4;
    maxiter = inf;
    
    % Generate x_orig
    I = randperm(n);
    x0 = zeros(n,1);
    x0(I(1:s)) = randn(s,1);
    
    % Generate A and b
    A = randn(m,n);
    
    noise = sigma0*(randn(m,1)./randn(m,1));
    b = A*x0 + noise;
    
    % Generate sigma by Cauchy loss function
    [confunc,~] = Lorentzian((noise).^2,delta);
    sigma = confunc*1.2;
    
    % Generate Slater point
    tstart1 = tic;
    [Q,R] = qr(A',0);
    time_QR = toc(tstart1);
    
    
    tstart2 = tic;
    xstar = Q*(R'\b);
    time_xstar = toc(tstart2);
    
    
    % Generate L
    if m > 2000
        clear opts
        opts.issym = 1;
        tstart3 = tic;
        L = eigs(A*A',1,'LM',opts);
        time_lambda = toc(tstart3);
    else
        tstart3 = tic;
        L = norm(A*A');
        time_lambda = toc(tstart3);
    end
     fprintf('\n Lipschitz constant L = %g, time_lambda = %g, time_QR = %g, time_xstar = %g \n', L, time_lambda, time_QR, time_xstar);
     
    
    % Save original data
    
    L0 = L;
    A0 = A;
    b0 = b;
    table2_L_time = [table2_L_time; L0, time_lambda, time_QR, time_xstar];
    
    % ADMM%    
    time_start = tic;   
    x = xstar;
    w = logpen(x,eps);
    
    bAx2 = zeros(m,1);
    [confunc,congrad] = Lorentzian(bAx2,delta);
    
    xADMM = zeros(n,1); % initial point
    zADMM = zeros(m,1);
    yADMM = zeros(m,1);
    
    tol_tau1 = 1;
    tol_mu1 = 1;
    
    inner_iter_totalADMM = 0;
    iter_totalADMM = 0;
    
    while iter_totalADMM < maxiter
        xold1 = x;
        
        v = sqrt(congrad);
        b = v.*b0;
        const = confunc - sum(congrad.*bAx2);
        barsigma = sigma - const;
        mu = sqrt(barsigma);
        
        L = max(congrad)*L0;
        beta = sqrt(1/L);
        
        tol_tau1 = tol_tau1/5;
        tol_tau = max(tol_tau1, 1e-8);
        sigmabd = min([mu, barsigma]);
        tol_mu1 = tol_mu1/1.2;
        tol_mu = max(tol_mu1, 1e-8);
        
        fprintf(' ***** ADMM beta = %g, L = %g, tol = %g *****\n', beta, L, sigmabd);
        [x,xsparse,iter2,~,xADMM,zADMM,yADMM] = ADMM_CS_LY_S(A0,v,b,mu,w,L,beta,sigmabd,inf,5000,xADMM,zADMM,yADMM, xstar,tol_mu,xold1,tol_tau);
        inner_iter_totalADMM = inner_iter_totalADMM + iter2;
        iter_totalADMM = iter_totalADMM +1;
        w = logpen(x,eps);
        bAx2 = (b0 - A0*x).^2;
        [confunc,congrad] = Lorentzian(bAx2,delta);
        
        if norm(x - xold1)/max(norm(xold1),1) <= tol_min
            fprintf(' D-IRL terminates\n')
            break
        end
        
        
    end
    
    time_ADMM = toc(time_start);
    fval_ADMM = sum(log(1 + abs(xsparse)/eps));
    rec_err_ADMM = norm(xsparse - x0)/max(norm(x0),1);
    residual_ADMM = (sum(log(1 + (A0*xsparse - b0).^2/delta^2)) - sigma)/sigma;
   
    if rec_err_ADMM <= 0.01
   
      iter_s_ADMM = inner_iter_totalADMM;      iter_f_ADMM = 0;
      time_s_ADMM = time_ADMM;                 time_f_ADMM = 0;
      rec_err_s_ADMM = rec_err_ADMM;           rec_err_f_ADMM = 0;
      s_ADMM = 1;
      
    else
      iter_s_ADMM = 0;                         iter_f_ADMM = inner_iter_totalADMM;      
      time_s_ADMM = 0;                         time_f_ADMM = time_ADMM;                 
      rec_err_s_ADMM = 0;                      rec_err_f_ADMM = rec_err_ADMM;           
      s_ADMM = 0;
    
    end
    
   
    fprintf(' D-IRL method terminated: time = %4.1f, total_inneriter = %d, total_outeriter = %d, nnz = %d,  fval = %8.7e, rec_err = %3.2e, residual = %2.1e  \n',...
        time_ADMM, inner_iter_totalADMM, iter_totalADMM, nnz(xsparse), fval_ADMM, rec_err_ADMM, residual_ADMM);
    
    table2_ADMM = [table2_ADMM; iter_s_ADMM, iter_f_ADMM, time_s_ADMM, time_f_ADMM, rec_err_s_ADMM, rec_err_f_ADMM, residual_ADMM, s_ADMM];   
    
 % SPGL1%
    fprintf(' ********* Start of SPGL1 ********** \n');
    time_start2 = tic;
    xspg = xstar;
    sw = logpen(xspg,eps);
    
    sbAx2 = zeros(m,1);
    [sconfunc,scongrad] = Lorentzian(sbAx2,delta);
    inner_iter_totalSpg = 0;
    iter_totalspg = 0;

    while iter_totalspg < maxiter
        sxold = xspg;

        sv = sqrt(scongrad);
        sb = sv.*b0;
        sconst = sconfunc - sum(scongrad.*sbAx2);
        
        smu = sqrt(sigma - sconst); 
        
        opts = spgSetParms('weights', sw,'primal_norm',@NormL1_primal,'dual_norm', @NormL1_dual, 'project',@NormL1_project, 'verbosity',0,'history','true');
        
        [spgu,~,~,info] = spg_bpdn( @spgA, sb, smu, opts);
        inner_iter_totalSpg = inner_iter_totalSpg + info.iter;
        iter_totalspg = iter_totalspg +1;
        stmp = norm(sb - spgA(spgu,1));
        if stmp > smu
            xspg = (1 - smu/stmp)*xstar + smu/stmp*spgu;
        else
            xspg = spgu;
        end
        
        sw = logpen(xspg,eps);
        sbAx2 = (b0 - A0*xspg).^2;
        [sconfunc,scongrad] = Lorentzian(sbAx2,delta);
        
        if norm(xspg - sxold)/max(norm(sxold),1) <= tol_min
            fprintf(' D-IRL terminates\n')
            break
        end
          
    end
       
   
        
    time_SPGL1 = toc(time_start2);
    fval_SPGL1 = sum(log(1 + abs(spgu)/eps));
    rec_err_SPGL1 = norm(spgu - x0)/max(norm(x0),1);
    residual_SPGL1 = (sum(log(1 + (A0*spgu - b0).^2/delta^2)) - sigma)/sigma;
    
     if rec_err_SPGL1 <= 0.01
   
      iter_s_SPGL1 = inner_iter_totalSpg;      iter_f_SPGL1 = 0;
      time_s_SPGL1 = time_SPGL1;                 time_f_SPGL1 = 0;
      rec_err_s_SPGL1 = rec_err_SPGL1;           rec_err_f_SPGL1 = 0;
      s_SPGL1 = 1;
      
    else
      iter_s_SPGL1 = 0;                         iter_f_SPGL1 = inner_iter_totalSpg;      
      time_s_SPGL1 = 0;                         time_f_SPGL1 = time_SPGL1;                 
      rec_err_s_SPGL1 = 0;                      rec_err_f_SPGL1 = rec_err_SPGL1;             
      s_SPGL1 = 0;
    
    end
    
    
   
    fprintf(' D-IRL method terminated: time = %4.1f, total_inneriter = %d, total_outeriter = %d,  nnz = %d,  fval = %8.7e, rec_err = %3.2e, residual = %2.1e  \n',...
        time_SPGL1, inner_iter_totalSpg, iter_totalspg, nnz(spgu), fval_SPGL1, rec_err_SPGL1, residual_SPGL1);
    
table2_SPGL1 = [table2_SPGL1; iter_s_SPGL1, iter_f_SPGL1, time_s_SPGL1, time_f_SPGL1, rec_err_s_SPGL1, rec_err_f_SPGL1, residual_SPGL1, s_SPGL1];   
    
   

 % SCP_ls%
    fprintf(' ********* Start of SCP_log_LL2 ********** \n');
    tstart3 = tic;
    [x_scp, iter_scp, talg] = SCP_log_LL2( A0,b0,sigma,eps,xstar,eps,delta,1e-5,inf);
    
    time_SCP = toc(tstart3);
    fval_SCP = sum(log(1 + abs(x_scp)/eps));
    rec_err_SCP = norm(x_scp - x0)/max(1, norm(x0));
    residual_SCP = (sum(log(1 + (A0*x_scp - b0).^2/delta^2)) - sigma)/sigma;
   
    if rec_err_SCP <= 0.01
   
      iter_s_SCP = iter_scp;                 iter_f_SCP = 0;
      time_s_SCP = time_SCP;                 time_f_SCP = 0;
      rec_err_s_SCP = rec_err_SCP;           rec_err_f_SCP = 0;
      s_SCP = 1;
      
    else
      iter_s_SCP = 0;                         iter_f_SCP = iter_scp;      
      time_s_SCP = 0;                         time_f_SCP = time_SCP;                 
      rec_err_s_SCP = 0;                      rec_err_f_SCP = rec_err_SCP;              
      s_SCP = 0;
    
    end
    
    fprintf(' SCP_log_LL2 terminated: time = %4.1f, iter = %d,  nnz = %d,  fval = %8.7e, rec_err = %3.2e, residual = %2.1e  \n',...
        time_SCP, iter_scp, nnz(x_scp), fval_SCP, rec_err_SCP, residual_SCP);
    
table2_SCP = [table2_SCP; iter_s_SCP, iter_f_SCP, time_s_SCP, time_f_SCP, rec_err_s_SCP, rec_err_f_SCP, residual_SCP, s_SCP];   
    
    
    
    
    
 end
sumnum_s_ADMM = sum(table2_ADMM(:,8));
sumnum_f_ADMM = repeats - sumnum_s_ADMM; 
table1_ADMM = [table1_ADMM; sumnum_s_ADMM/repeats*100, sum(table2_ADMM(:,1))/sumnum_s_ADMM,  sum(table2_ADMM(:,2))/sumnum_f_ADMM, sum(table2_ADMM(:,3))/sumnum_s_ADMM,  sum(table2_ADMM(:,4))/sumnum_f_ADMM, ...
                            sum(table2_ADMM(:,5))/sumnum_s_ADMM,  sum(table2_ADMM(:,6))/sumnum_f_ADMM, min(table2_ADMM(:,7)),  max(table2_ADMM(:,7))];
                        
sumnum_s_SPGL1 = sum(table2_SPGL1(:,8));
sumnum_f_SPGL1 = repeats - sumnum_s_SPGL1;                        
table1_SPGL1 = [table1_SPGL1; sumnum_s_SPGL1/repeats*100, sum(table2_SPGL1(:,1))/sumnum_s_SPGL1,  sum(table2_SPGL1(:,2))/sumnum_f_SPGL1, sum(table2_SPGL1(:,3))/sumnum_s_SPGL1,  sum(table2_SPGL1(:,4))/sumnum_f_SPGL1, ...
                              sum(table2_SPGL1(:,5))/sumnum_s_SPGL1,  sum(table2_SPGL1(:,6))/sumnum_f_SPGL1, min(table2_SPGL1(:,7)),  max(table2_SPGL1(:,7))];                        
                        
sumnum_s_SCP = sum(table2_SCP(:,8));
sumnum_f_SCP = repeats - sumnum_s_SCP;                        
table1_SCP = [table1_SCP; sumnum_s_SCP/repeats*100, sum(table2_SCP(:,1))/sumnum_s_SCP,  sum(table2_SCP(:,2))/sumnum_f_SCP, sum(table2_SCP(:,3))/sumnum_s_SCP,  sum(table2_SCP(:,4))/sumnum_f_SCP, ...
                          sum(table2_SCP(:,5))/sumnum_s_SCP,  sum(table2_SCP(:,6))/sumnum_f_SCP, min(table2_SCP(:,7)),  max(table2_SCP(:,7))];     
                      
table1_L_time = [table1_L_time; mean(table2_L_time)];                      
   


end

% Results = 'D:\Matlab1';
% mkdir(Results);

a = clock;
% fname = [Results, 'Results_SSQ_ADMM' '_' date '_' int2str(a(4)) '-' int2str(a(5)) '.txt'];
fname = ['Results_SSQ_ADMM' '_' date '_' int2str(a(4)) '-' int2str(a(5)) '.txt'];
fid1 = fopen(fname, 'w');

for ii = 1:repeats_index
    fprintf(fid1, '& %3.0f  & %13.0f  & %13.0f  & %6.1f  & %6.1f  & %2.1e  & %2.1e  & %2.1e  & %2.1e \n', table1_ADMM(ii,:));
end

fclose(fid1);  




b = clock;
% fname = [Results, 'Results_SSQ_SPGL1' '_' date '_' int2str(b(4)) '-' int2str(b(5)) '.txt'];
fname = ['Results_SSQ_SPGL1' '_' date '_' int2str(b(4)) '-' int2str(b(5)) '.txt'];
fid2 = fopen(fname, 'w');

for ii = 1:repeats_index
    fprintf(fid2, '& %3.0f  & %13.0f  & %13.0f  & %6.1f  & %6.1f  & %2.1e  & %2.1e  & %2.1e  & %2.1e \n', table1_SPGL1(ii,:));
end

fclose(fid2);  



c = clock;
% fname = [Results, 'Results_SSQ_SCP' '_' date '_' int2str(c(4)) '-' int2str(c(5)) '.txt'];
fname =  ['Results_SSQ_SCP' '_' date '_' int2str(c(4)) '-' int2str(c(5)) '.txt'];
fid3 = fopen(fname, 'w');

for ii = 1:repeats_index
    fprintf(fid3, '& %3.0f  & %13.0f  & %13.0f  & %6.1f  & %6.1f  & %2.1e  & %2.1e  & %2.1e  & %2.1e \n', table1_SCP(ii,:));
end

fclose(fid3);  



d = clock;
% fname = [Results, 'Results_SSQ_L_time' '_' date '_' int2str(d(4)) '-' int2str(d(5)) '.txt'];
fname =  ['Results_SSQ_L_time' '_' date '_' int2str(d(4)) '-' int2str(d(5)) '.txt'];
fid4 = fopen(fname, 'w');

for ii = 1:repeats_index
    fprintf(fid4, '& %3.2e  & %6.1f  & %6.1f  & %6.1f  \n', table1_L_time(ii,:));
end

fclose(fid4);  

save results table1_ADMM table1_SPGL1 table1_SCP table1_L_time



% nested function

function out = spgA(xin,cases)

if cases == 1
    out = sv.*(A0*xin);
else 
    out = A0'*(sv.*xin);
end
end

end