% Test problem

clear all
clc

randn('seed',1);
rand('seed',1);

fid_1 = fopen('Table_1.txt','w');
fid_2 = fopen('Table_2.txt','w');

repeats = 20;

lambda_list = [5e-3 1e-3 5e-4];
n_list = [3000 6000 9000];
mu = 0.99;
sigma = 0.01;
tol = 1e-4;


%% Begin loop for lambda and n
for h = 1: length(lambda_list)
    lambda = lambda_list(h);
    for i = 1: length(n_list)
        n = n_list(i);
        m = n/5;
        k = m/4;
        r = m/20;
        
        Table_1 = []; Table_2 = []; 
        
        for repeat = 1: repeats
            
            A = randn(m+r,n);
            for ii = 1:n
                A(:,ii) = A(:,ii)/norm(A(:,ii));
            end
            
            z = [zeros(m,1);8*ones(r,1)];
            err = sigma*randn(m+r,1);
            
            I = randperm(n);
            J = I(1:k);
            x = zeros(n,1);
            x(J) = randn(k,1);
            
            b = A*x - z + err;
            
            time_eig = 0;
            if m > 2000
                clear opts
                opts.issym = 1;
                tstart = tic;
                L = eigs(A*A',1,'LM',opts);
                time_eig = time_eig + toc(tstart);
            else
                tstart = tic;
                L = norm(A*A');
                time_eig = time_eig + toc(tstart);
            end
            fprintf('\n Lipschitz constant L = %g, time_eig = %g\n', L, time_eig)
            
            
            %% For Table 1
            ktest = 0.8*k;
            rtest = r;
            fprintf('******** pDCAe for Least squares and Truncated regularization for Table 1 ********\n')
            clear opts
            opts.x0 = zeros(n,1);
            opts.maxiter = inf;
            opts.tol = tol;
            opts.freq = 20000;
            opts.disp = 1;
            [x_1, z_1, iter_1, fval_1, time_1] = p_DCAe(A, b, ktest, rtest, lambda, mu, L, opts);
            
            fprintf('******** NPG for Least squares and Truncated regularization for Table 1********\n')
            clear opts
            opts.x0 = zeros(n,1);
            opts.maxiter = inf;
            opts.tol = tol;
            opts.freq = 20000;
            opts.disp = 1;
            opts.M = 5;
            [x_npg_1, z_npg_1, iter_npg_1, fval_npg_1, time_npg_1] = NPG_DC(A, b, ktest, rtest, lambda, mu, L, opts);
               
            
            %% For Table 2
            ktest = 0.8*k;
            rtest = 1.1*r;
            fprintf('******** pDCAe for Least squares and Truncated regularization for Table 2 ********\n')
            clear opts
            opts.x0 = zeros(n,1);
            opts.maxiter = inf;
            opts.tol = tol;
            opts.freq = 20000;
            opts.disp = 1;
            [x_2, z_2, iter_2, fval_2, time_2] = p_DCAe(A, b, ktest, rtest, lambda, mu, L, opts);
            
            fprintf('******** NPG for Least squares and Truncated regularization for Table 2********\n')
            clear opts
            opts.x0 = zeros(n,1);
            opts.maxiter = inf;
            opts.tol = tol;
            opts.freq = 20000;
            opts.disp = 1;
            opts.M = 5;
            [x_npg_2, z_npg_2, iter_npg_2, fval_npg_2, time_npg_2] = NPG_DC(A, b, ktest, rtest, lambda, mu, L, opts);
            
            Table_1 = [Table_1; norm(x_1 - x)/sqrt(n) iter_1 iter_npg_1 fval_1 fval_npg_1 time_1 time_npg_1];
            Table_2 = [Table_2; norm(x_2 - x)/sqrt(n) iter_2 iter_npg_2 fval_2 fval_npg_2 time_2 time_npg_2];
        end
        fprintf(fid_1,'%4.0e & %5.0d & %5.0d & %5.0d & %5.0d & %3.1e & %5.0f & %5.0f & %7.4e & %7.4e & %5.1f & %5.1f \n', lambda, m, n, k, r, mean(Table_1));
        fprintf(fid_2,'%4.0e & %5.0d & %5.0d & %5.0d & %5.0d & %3.1e & %5.0f & %5.0f & %7.4e & %7.4e & %5.1f & %5.1f \n', lambda, m, n, k, r, mean(Table_2));
    end
end

fclose(fid_1);
fclose(fid_2);



