% An outlier detection problem

clear all;
clc

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

fid_1 = fopen('outlier.txt','w');

repeats = 20;

n_list = 10000:2000:30000;
tol = 1e-8;
maxiter = 3000;
sigma = 1e-4;

%% Begin loop for n
for i = 1: length(n_list)
    n = n_list(i);
    m = n/5;
    s = n/20;
    r = m/20;
    
    progress_1 = []; progress_2= [];
    
    fopen(fid_1);
    
    for repeat = 1:repeats
        I = randperm(n);
        J = I(1:s);
        w = zeros(n,1);
        w(J) = randn(s,1);
        A = randn(m,n);
        
        for j = 1:n
            A(:,j) = A(:,j)/norm(A(:,j));
        end
        
        b = A*w;
        b(m-r+1:end) = b(m-r+1:end) + 10*sign(randn(r,1));
        
        if m > 2000
            clear opts
            opts.issym = 1;
            tstart = tic;
            L = eigs(A*A',1,'LM',opts);
            time_eig = toc(tstart);
        else
            tstart = tic;
            L = norm(A*A');
            time_eig = toc(tstart);
        end
        fprintf('\n Lipschitz constant L = %g, time_eig = %g\n', L, time_eig)
        
        fprintf('******** NPG_BB for SFP ********\n')
        
        clear opts
        opts.x0 = zeros(n,1);
        opts.maxiter = inf;
        opts.tol = tol;
        opts.freq = 5000;
        opts.disp = 0;
        opts.M = 5;
        tstart = tic;
        [x_1, fnew_1, iter_1, ~] = NPG_BB(A, b, r, s, L, opts);
        time_1 = toc(tstart);
        progress_1 = [progress_1; sqrt(fnew_1*2) iter_1 time_1];
        
        fprintf('******** NPG_constant for SFP ********\n')
        
        clear opts
        opts.x0 = zeros(n,1);
        opts.maxiter = maxiter;
        opts.tol = tol;
        opts.freq = 5000;
        opts.disp = 0;
        tstart = tic;
        [x_2, fnew_2, iter_2, ~] = NPG_constant(A, b, r, s, L, opts);
        time_2 = toc(tstart);
        progress_2 = [progress_2; sqrt(fnew_2*2) iter_2 time_2];
        
    end
    fprintf(fid_1, '%5d & %5d & %5d & %5d & %1.0e & %5.0f & %5.1f & %1.0e & %5.0f & %5.1f\n', n, m, s, r, mean(progress_1), mean(progress_2));
end
fclose(fid_1);