% generate data
clear all
clc

tol = 1e-8;
maxiter = inf;
freq = inf;
repeat = 50;

% parameter for PR
eta = 0.93;
% compare DR and PR
ratio = 150;
beta = 2.2;
% compare different DR
ratio_1 = 50;
ratio_2 = 100;


randn('seed',200); % fix seed
rand('seed',200); % fix seed

fid = fopen('0.93_2.2.txt','w');
 fid1 = fopen('50_100.txt','w');

for m = 100:100:500
  for n = 4000:1000:6000

    k = ceil(m/5);  % hard problem
    
    timeDR = 0;    timePR = 0;     timeDR_1 = 0;    timeDR_2 = 0;
    iterDR = 0;    iterPR = 0;     iterDR_1 = 0;    iterDR_2 = 0;
    fDRmax = 0;    fPRmax = 0;     fDRmax_1 = 0;    fDRmax_2 = 0;
    fDRmin = inf;  fPRmin = inf;   fDRmin_1 = inf;  fDRmin_2 = inf;
    successDR = 0; successPR = 0;  successDR_1 = 0; successDR_2 = 0;
    failDR = 0;    failPR = 0;     failDR_1 = 0;    failDR_2 = 0;
    iter_succDR = 0; iter_succPR = 0;   iter_succDR1 = 0; iter_succDR2 = 0;
    iter_failDR = 0; iter_failPR = 0;   iter_failDR1 = 0; iter_failDR2 = 0;

    for ii = 1:repeat

      A = randn(m,n);
      x0 = max(min(randn(k,1),1e6),-1e6);
      x1 = zeros(n,1);
      I = randperm(n); x1(I(1:k)) = x0;
      b = A*x1;

      fprintf(' Start DR 150\n')
      tDR = tic;
      [fval,z,x,y,iter] = DR_proj(A,b,k,tol,maxiter,ratio);
      timeDR = timeDR + toc(tDR);
      iterDR = iterDR + iter;
      fDRmax = max(fDRmax,fval);
      fDRmin = min(fDRmin,fval);
      if fval < 1e-12
        successDR = successDR + 1;
        iter_succDR = iter_succDR + iter;
      end
      if fval > 1e-6
        failDR = failDR + 1;
        iter_failDR = iter_failDR + iter;
      end


      fprintf(' Start PR 2.2\n')
      tPR = tic;
      [fval2,z2,x2,y2,iter2] = PR_proj_hp(A,b,k,tol,maxiter,beta,eta);
      timePR = timePR + toc(tPR);
      iterPR = iterPR + iter2;
      fPRmax = max(fPRmax,fval2);
      fPRmin = min(fPRmin,fval2);
      if fval2 < 1e-12
        successPR = successPR + 1;
        iter_succPR = iter_succPR + iter2;
      end
      if fval2 > 1e-6
        failPR = failPR + 1;
        iter_failPR = iter_failPR + iter2;
      end

      
      fprintf(' Start DR 50\n')
      tDR1 = tic;
      [fval3,z3,x3,y3,iter3] = DR_proj(A,b,k,tol,maxiter,ratio_1);
      timeDR_1 = timeDR_1 + toc(tDR1);
      iterDR_1 = iterDR_1 + iter3;
      fDRmax_1 = max(fDRmax_1,fval3);
      fDRmin_1 = min(fDRmin_1,fval3);
      if fval3 < 1e-12
        successDR_1 = successDR_1 + 1;
        iter_succDR1 = iter_succDR1 + iter3;
      end
      if fval3 > 1e-6
        failDR_1 = failDR_1 + 1;
        iter_failDR1 = iter_failDR1 + iter3;
      end
      
      
      fprintf(' Start DR 100\n')
      tDR2 = tic;
      [fval4,z4,x4,y4,iter4] = DR_proj(A,b,k,tol,maxiter,ratio_2);
      timeDR_2 = timeDR_2 + toc(tDR2);
      iterDR_2 = iterDR_2 + iter4;
      fDRmax_2 = max(fDRmax_2,fval4);
      fDRmin_2 = min(fDRmin_2,fval4);
      if fval4 < 1e-12
        successDR_2 = successDR_2 + 1;
        iter_succDR2 = iter_succDR2 + iter4;
      end
      if fval4 > 1e-6
        failDR_2 = failDR_2 + 1;
        iter_failDR2 = iter_failDR2 + iter4;
      end
      
      
 
    end

    fprintf(fid,' %5d & %5d &  %6.0f & %1.0e & %1.0e & %2.0f & %2.0f & %6.0f & %6.0f &  %6.0f & %1.0e & %1.0e & %2.0f & %2.0f &  %6.0f & %6.0f\n',...
      m, n, iterDR/repeat, fDRmax, fDRmin, successDR, failDR, iter_succDR/successDR, iter_failDR/failDR,...
            iterPR/repeat, fPRmax, fPRmin, successPR, failPR,  iter_succPR/successPR, iter_failPR/failPR);
        
    fprintf(fid1,' %5d & %5d &  %6.0f & %1.0e & %1.0e & %2.0f & %2.0f & %6.0f & %6.0f &  %6.0f & %1.0e & %1.0e & %2.0f & %2.0f &  %6.0f & %6.0f\n',...
      m, n, iterDR_1/repeat, fDRmax_1, fDRmin_1, successDR_1, failDR_1, iter_succDR1/successDR_1, iter_failDR1/failDR_1,...
            iterDR_2/repeat, fDRmax_2, fDRmin_2, successDR_2, failDR_2, iter_succDR2/successDR_2, iter_failDR2/failDR_2);
  end

end

fclose(fid);