clear all;

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

% algorithm parameter
sigma = 1e-2;
M2 = 5;
mu = 0.001;
re_freq = 200;
dis_freq = inf;
maxiter = 5000;
tol = 1e-5;

fid = fopen('datarandmu1e3.txt','w');

%% problem setting
for index = 1:10
    progress = [];   
    for repeats = 1:30
        
        
       %% Problem generation
        % signal length
        n = 2560*index;
        % number of spikes to put down
        k = 80*index; % originally 80*index
        % number of observations to make
        m = 720*index;
        lambda = 1;
        
        I = randperm(n);
        J = I(1:k);
        x = zeros(n,1);
        x(J) = randn(k,1);
        
        A = randn(m,n);
        for j = 1:n
            A(:,j) = A(:,j)/norm(A(:,j));
        end
        b = A*x + sigma*randn(m,1);
        
        if mu >= 0.5*norm(A'*b,inf) % check the differentiability assumption
            fprintf(fid,' Conditions violated\n');
            pause
        end
        
        if m > 2000
            clear opts
            opts.issym = 1;
            tstart = tic;
            nmA = eigs(A*A',1,'LM',opts);
            time_lambda = toc(tstart);
        else
            tstart = tic;
            nmA = norm(A*A');
            time_lambda = toc(tstart);
        end
        fprintf('\n Lipschitz constant L = %g\n', nmA)
        

      %% NPG
        clear opts
        opts.x0 = zeros(n,1);
        opts.freq = 5000;
        opts.maxiter = inf;
        opts.disp = 0;
        opts.M = M2;
        opts.tol = tol;
        fprintf(' ********** L1 - L2: Majorization with line search**********\n')
        
        [x1,iter1,time1] = NPG_ls_l1_l2(A,b,mu,lambda,opts);
        fprintf(' iter = %5d  time = %5.1f  sparsity = %3d\n', iter1, time1, sum(abs(x1) > 1e-5));
        
        fval1 = 0.5*norm(A*x1 - b)^2 + mu*(norm(x1,1) - norm(x1));
        fprintf(' fval = %12.8f\n', fval1)   
        
        %% Nesterov with restart
        clear opts
        tic
        [fval2,x_new2,iter2] = pro_gra_extra_l1l2( A, b, nmA, mu, tol, inf, re_freq, dis_freq);
        time2 = toc;
        
        fprintf(' iter = %5d  time = %5.1f  sparsity = %3d\n', iter2, time2, sum(abs(x_new2) > 1e-5));
        fprintf(' fval = %12.8f\n', fval2)
        
       %% Proximal Gradient Algorithm
        clear opts
        tic
        [fval3,x_new3,iter3] = pro_gra_l1l2( A, b, nmA, mu, tol, maxiter);
        time3 = toc;
        
        fprintf(' iter = %5d  time = %5.1f  sparsity = %3d\n', iter3, time3, sum(abs(x_new3) > 1e-5));
        fprintf(' fval = %12.8f\n', fval3)
        
        progress = [progress; time_lambda iter1 iter2 iter3 time1 time2 time3 fval1 fval2 fval3];
        
    end
    fprintf(fid,' %5.0f & %5.0f & %5.0f & %5.1f & %5.0f & %5.0f & %5.0f & %5.1f & %5.1f & %5.1f & %5.4e & %5.4e & %5.4e \n', n, m, k, mean(progress));
    
end
fclose(fid);