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

arg_list = argv ();
if (size(arg_list)(1)!=7)
  printf("Usage: compare_system_order_with_fault [m-files-dir out-dir graph-file b-file v kstart kstop]\n");
  exit(0);
endif

dir = arg_list{2};
addpath(arg_list{1});
load(arg_list{3},"A");
load(arg_list{4},"b");
n = size(A)(1);
A1=A=make_col_stochastic(A);
printf("compute gap\n");
g=1 - flipud(sort(abs(eigs(A,2,"lm"))))(2);
printf("compute fault\n");
v=str2num(arg_list{5});
v_input_set = find(A1(v,:));
for j = v_input_set
  r = 1 - A1(v,j);
  if r == 0
    printf("error: killing node %d makes node %d a dead end.\n", v, j);
  else
    A1(:,j) = A1(:,j)/r;
  endif
endfor
A1(:,v) = zeros(n,1);
A1(v,:) = zeros(1,n);
printf("compute gap after fault\n");
g1 = 1 - flipud(sort(abs(eigs(A1,2,"lm"))))(2);

kstart=str2num(arg_list{6});
kstop=str2num(arg_list{7});
printf("A gap = %f\n",g);
printf("gap after the failure = %f\n",g1);

#save gaps
if(kstart==4)
  printf("kstart = 4\n");
  gaps=[g g1];
  save(strcat(dir,"gaps.mat"),"gaps");
endif

#execute the experiment for k=4:4:120 (instead of 4:2:30,34:4:120)

for k1=kstart:4:kstop
  kmax=120-k1;
  VTol =[];
  Vrank= [];
  H=[];
  for k2=0:4:kmax
    printf("k1=%d k2=%d\n",k1,k2);
    H = impulse_response_with_detected_failure(A,b,k1,k2,v);
    GTol=[];
    Grank=[];
    for i=1:n
      if i!=v
        GTol = [ GTol ; gap_estimate_nonzero_b(H(i,:)) ];
        Grank = [ Grank ; gap_estimate_nonzero_b_with_rank(H(i,:)) ];
      endif
    endfor
    VTol=[VTol GTol];
    Vrank=[Vrank Grank];
  endfor
  printf("saving data\n");
  save(strcat(dir,"/k1-",num2str(k1),"-TOL.mat"),"VTol");
  save(strcat(dir,"/k1-",num2str(k1),"-rank.mat"),"Vrank");
endfor

