% run code
clc
clear
maxiter = inf; % maximum number of iterations
m = 250;
repeats = 10;
fid = fopen('flassoppg.txt','w');
fclose(fid);
fid2 = fopen('flassoprox.txt','w');
fclose(fid2);

for n = [10000 20000 30000]
  for alpha = [1e-4 3e-4 5e-4]
    
    PPG = [];
    PROX = [];
    
    for repeat = 1:repeats
      lambda1 = alpha*m; 
      lambda2 = 100*lambda1; % differ by 100 times
      
      diff = 0;
      while diff == 0
        A = randn(m,n);
        for i = 1:n-1
          A(:,i) = A(:,i)/norm(A(:,i)); % normalize each column
        end
     
        x = sparse(n,1);
        x(1:20) = 20*randn(1);
        x(41) = 30*randn(1);
        x(71:85) = 10*randn(1);
        x(121:125) = 20*randn(1);
      
        A(:,end) = -sign(A*x + rand(1)); % -b_i
        for i = 1:n-1
          A(:,i) = A(:,i).*A(:,end);
        end
        diff = max(A(:,end)) - min(A(:,end));
      end
      
      [zbest,fbest,dval,dfeas,iter,time] = ppg_logistic_flasso(A,lambda1,lambda2,1e-4,maxiter,500);
      PPG = [PPG;iter time fbest dval dfeas];
      
      [zbest,fbest,dval,dfeas,iter,time] = mFBS_logistic_flasso_uv(A,lambda1,lambda2,1e-4,20000,500);      
      PROX = [PROX;iter time fbest dval dfeas];
    end
    
    PPG = mean(PPG,1);
    PROX = mean(PROX,1);
    
    fid = fopen('flassoppg.txt','a');
    fprintf(fid,' %3d & %6d & %2.1e & %5.0f & %5.1f & %0.3e/%0.3e/%0.1e\n', m, n, alpha, PPG);
    fclose(fid);
    fid2 = fopen('flassoprox.txt','a');
    fprintf(fid2,' %3d & %6d & %2.1e & %5.0f & %5.1f & %0.3e/%0.3e/%0.1e\n', m, n, alpha, PROX);
    fclose(fid2);
  end
end
















% fprintf(' nnz(zbest) %d, %d, nnz(M zbest), %d, %d\n', sum(abs(zbest(1:end-1)) > 1e-3), sum(abs(zbest(1:end-1)) > 1e-5),...
%   sum(abs(Mmap(zbest,1)) > 1e-3), sum(abs(Mmap(zbest,1)) > 1e-5))

% profile on
% [zbest,fbest,dval,dfeas,iter,time] = prox_logistic_flasso_yz(A,lambda1,lambda2,1e-4,inf,500);
% profile off
% 
% fprintf(' nnz(zbest) %d, %d, nnz(M zbest), %d, %d\n', sum(abs(zbest(1:end-1)) > 1e-4), sum(abs(zbest(1:end-1)) > 1e-5),...
%   sum(abs(Mmap(zbest,1)) > 1e-4), sum(abs(Mmap(zbest,1)) > 1e-5))