#! /usr/bin/octave-3.2.3 -qf

arg_list = argv ();
if (size(arg_list)(1)!=2)
  printf("Usage: compare_system_order [m-files-dir out-dir]\n");
  exit(0);
endif

dir = arg_list{2};
addpath(arg_list{1});

load(strcat(dir,"/b.mat"),"b");
load(strcat(dir,"/graph.mat"),"A");
A=make_col_stochastic(A);
gap=1 - flipud(sort(abs(eigs(A,2,"lm"))))(2);
save(strcat(dir,"/gap.mat"),"gap");

Vtol=[];
Vn=[];
Vrank=[];
for k=[4:4:120]
  printf("k=%d",k);
  H = impulse_response(A,b,k);
  #with TOL
  printf("...TOL");
  G = [];
  for i = 1:10000
    G = [ G ; gap_estimate_nonzero_b(H(i,:)) ];
  endfor
  Vtol=[Vtol G];
  printf(" done");
  #without TOL
  #printf("...noTOL");
  #G = [];
  #for i = 1:1000
  #  G = [ G ; gap_estimate_nonzero_b_without_TOL(H(i,:)) ];
  #endfor
  #Vn=[Vn G];
  #printf(" done");
  #with Rank
  printf("...RANK");
  G = [];
  for i = 1:10000
    G = [ G ; gap_estimate_nonzero_b_with_rank(H(i,:)) ];
  endfor
  Vrank=[Vrank G];
  printf(" done\n");
endfor

save(strcat(dir,"/data_TOL.mat"),"Vtol");
#save(strcat(dir,"/er/data_noTOL.mat"),"Vn");
save(strcat(dir,"/data_rank.mat"),"Vrank");


printf("process results\n");

#load(strcat(dir,"er/data_noTOL.mat"),"Vn");
load(strcat(dir,"/data_TOL.mat"),"Vtol");
load(strcat(dir,"/data_rank.mat"),"Vrank");
#load(strcat(dir,"er/graph.mat"),"A");
e=log(10^-6);

load(strcat(dir,"/gap.mat"),"gap");
g=gap;
l=1-g;
t=e/log(l);

for z=1:2
  if (z==1)
    G=Vtol;
    printf("TOL\n");
  elseif (z==2)
    G=Vrank;
    printf("rank\n");
  endif

  L=1-G;
  T=[];
  for i=1:size(L)(1)
    for j=1:size(L)(2)
      T(i,j)=e/log(L(i,j));
    endfor
    #T(:,i)=e./log(L(:,i));
  endfor

  #save("check.dat","T");

  i=0;
  RelG=[];
  RelL=[];
  RelT=[];
  AbsG=[];
  AbsL=[];
  AbsT=[];

  #RelGgossip=[];
  #RelLgossip=[];
  #RelTgossip=[];
  #AbsGgossip=[];
  #AbsLgossip=[];
  #AbsTgossip=[];

  for k1=[4:4:120]
    printf("k1=%d\n",k1);
    i=i+1;
    ##gap
    data=G(:,i);
    err=abs(g-data);
    q=quantile(err, [0.1 0.5 0.9]);
    q=[k1 q'];
    AbsG=[AbsG; q];
    err=(abs(data-g)/g)*100;
    q=quantile(err, [0.1 0.5 0.9]);
    q=[k1 q'];
    RelG=[RelG; q];
  
    
    ##lambda
    data=L(:,i);
    err=abs(l-data);
    q=quantile(err, [0.1 0.5 0.9]);
    q=[k1 q'];
    AbsL=[AbsL; q];
    err=(abs(data-l)/l)*100;
    q=quantile(err, [0.1 0.5 0.9]);
    q=[k1 q'];
    RelL=[RelL; q];
    ##tau
    data=T(:,i);
    err=abs(t-data);
    q=quantile(err, [0.1 0.5 0.9]);
    q=[k1 q'];
    AbsT=[AbsT; q];
    err=(abs(data-t)/t)*100;
    q=quantile(err, [0.1 0.5 0.9]);
    q=[k1 q'];
    RelT=[RelT; q];
  endfor

  if (z==1)
    #printf("save noTOL\n");
    #save(strcat(dir,"er/rel_gap_noTOL.mat"),"RelG");
    #save(strcat(dir,"er/rel_lambda_noTOL.mat"),"RelL");
    #save(strcat(dir,"er/rel_tau_noTOL.mat"),"RelT");
    #save(strcat(dir,"er/abs_gap_noTOL.mat"),"AbsG");
    #save(strcat(dir,"er/abs_lambda_noTOL.mat"),"AbsL");
    #save(strcat(dir,"er/abs_tau_noTOL.mat"),"AbsT");
  #elseif (z==2)
    printf("save TOL\n");
    save(strcat(dir,"/rel_gap_TOL.mat"),"RelG");
    save(strcat(dir,"/rel_lambda_TOL.mat"),"RelL");
    save(strcat(dir,"/rel_tau_TOL.mat"),"RelT");
    save(strcat(dir,"/abs_gap_TOL.mat"),"AbsG");
    save(strcat(dir,"/abs_lambda_TOL.mat"),"AbsL");
    save(strcat(dir,"/abs_tau_TOL.mat"),"AbsT");
  elseif (z==2)
    printf("save rank\n");
    save(strcat(dir,"/rel_gap_rank.mat"),"RelG");
    save(strcat(dir,"/rel_lambda_rank.mat"),"RelL");
    save(strcat(dir,"/rel_tau_rank.mat"),"RelT");
    save(strcat(dir,"/abs_gap_rank.mat"),"AbsG");
    save(strcat(dir,"/abs_lambda_rank.mat"),"AbsL");
    save(strcat(dir,"/abs_tau_rank.mat"),"AbsT");
  endif
endfor

