reset;
param ver symbolic = "ac36";

param T integer, >=1;
param N integer, >=2;

param layers{1..T} default 0;
param pump{t in 1..T, 1..layers[t]};
set BOTTLENECKS{t in 1..T,1..layers[t]} dimen 2;
set PATH{t in 1..T,1..layers[t]} dimen 2;
set INITLINKS{1..T} dimen 2;
param load{t in 1..T, l in 1..layers[t], (i,j) in PATH[t,l]} >=0;

data dump/path.txt;

commands aron-config.txt;

param bottlnum{t in 1..T, l in 1..layers[t]} = card(BOTTLENECKS[t,l]);

param incrabs{t in 1..T, l in 1..layers[t]} = (if l>=2 then incrabs[t,l-1] else 1)*pump[t,l];

set TC =
# {0} union
# setof{i in 5..3 by -1} round(10^-i,10) union
  setof{i in 0.033..0.08 by 0.003} round(i,10)
    within interval [0,1);

set P =
  TC union
    setof{t in 1..T, l in 1..layers[t]: incrabs[t,l]>1}
      1/incrabs[t,l]
  union
    setof{t in 1..T, l in 1..layers[t],(i,j) in PATH[t,l] diff BOTTLENECKS[t,l]: incrabs[t,l]>1}
      load[t,l,i,j]/incrabs[t,l];

#===FEC===
param minN{p in P} = max(M,if p=0 then 0 else floor(log10(DER/M)/log10(p)));

param maxN_probe{p in P, pr in 1..4} = minN[p]*2^pr;

param Binomial_probe{p in P, pr in 1..4} =
  if M>=2 then
    prod{i in 1..M-1}
    (
      (maxN_probe[p,pr]-M+1+i)/i *
        p^((maxN_probe[p,pr]-M+1)/(M-1))*(1-p)
    )
  else
    p^maxN_probe[p,pr];

param maxN{p in P} = maxN_probe[p, min{pr in 1..4: Binomial_probe[p,pr]<=DER/M} pr];


set NN{p in P} = {minN[p]..maxN[p]};


param Binomial1{p in P, Ntr in NN[p]}=
  if M>=2 then
    prod{i in 1..M-1}((Ntr-M+1+i)/i*p^((Ntr-M+1)/(M-1))*(1-p))
  else
    p^Ntr;

param pbyq{p in P} = p/(1-p);

param Binomial{p in P, Ntr in NN[p], n in Ntr-M+1..Ntr} =
  if n=Ntr-M+1
  then
    Binomial1[p,Ntr]
  else
    Binomial[p,Ntr,n-1]*(Ntr-n+1)/n*pbyq[p];

param Failure{p in P, Ntr in NN[p]} = sum{i in Ntr-M+1..Ntr}Binomial[p,Ntr,i];

param goodN{p in P} = min{Ntr in NN[p]: Failure[p,Ntr]<=DER} Ntr;


param FEC{p in P} default
 if goodN[p]-1 in NN[p] and Failure[p,goodN[p]-1]>DER
   then
     goodN[p]-1 +
       (log(Failure[p,goodN[p]-1]) - log(DER)) /
       (log(Failure[p,goodN[p]-1]) - log(Failure[p,goodN[p]]))
   else goodN[p];
#===

#param FEC{p in P} = M + if p=0 then 0 else log(DER/M)/log(p);

param ARN{tc in TC, p in P} =
  if p<=tc
   then
     0
   else
     (FEC[p]-FEC[tc])/FEC[tc];

param benchmark_start;
param benchmark_stop;
param benchmark_info symbolic;
param logrec symbolic =
  sprintf("%s, %s, %s, done in %f seconds",
    ver,sub(ctime(),"^ *[^ ]* *",""),benchmark_info,benchmark_stop-benchmark_start);
param logfile symbolic = "aron.log";


param aron_b{tc in TC, t in 1..T, l in 1..layers[t]} =
  (if l >=2 then aron_b[tc,t,l-1]) +
  if incrabs[t,l]>1
    then
      bottlnum[t,l] * ARN[tc,1/incrabs[t,l]];

param aron_p{tc in TC, t in 1..T, l in 1..layers[t]} =
  sum{(i,j) in PATH[t,l] diff BOTTLENECKS[t,l]: incrabs[t,l]>1}
    ARN[tc,load[t,l,i,j]/incrabs[t,l]];

param aron{tc in TC, t in 1..T, l in 1..layers[t]} =
  aron_b[tc,t,l] + aron_p[tc,t,l];

param shift{t in 1..T} = if pump[t,1]<=1 then 1 else 0;
param chunk = (T + chunks-1) div chunks;
set FRAMES{i in 1..chunks} = {t in (i-1)*chunk+1..min(i*chunk,T): 1+shift[t] in 1..layers[t]};
param frames{i in 1..chunks} = card(FRAMES[i]);
set CHUNKS = {i in 1..chunks: frames[i]>=1};

param maxlayers{i in CHUNKS} = max{t in FRAMES[i]} (layers[t] - shift[t]);

#average ARON per frame
param aron_f{tc in TC, i in CHUNKS, l in 1..maxlayers[i]} =
  1/frames[i] * sum{t in FRAMES[i]} aron[tc,t,min(l+shift[t],layers[t])];

param links{t in 1..T, l in 1..layers[t]} = if pump[t,l] > 1 then
  card(PATH[t,l]) + sum{ll in 1..l-1: pump[t,ll]>1} card(BOTTLENECKS[t,ll]);

param aron_l{tc in TC, i in CHUNKS, l in 1..maxlayers[i]} =
    ( sum{t in FRAMES[i]} aron[tc,t,min(l+shift[t],layers[t])] )
          /
      ( sum{t in FRAMES[i]} links[t,min(l+shift[t],layers[t])] );

param dir symbolic;
param file symbolic;

let dir := sprintf("%s-aron",ver);
shell ("rmdir /s /q " & dir);
shell ("mkdir " & dir);
shell ("del /q " & dir & "\*");
shell ("copy " & ver & "*.txt " & dir);
shell ("copy config.txt " & dir);
shell ("copy aron-config.txt " & dir);

for{benchmark in {"ARON-Binomial", "ARON-Binomial-perlink", "ARON-Shannon", "ARON-Shannon-perlink"}}
{
  let file := dir & "/aron-" & benchmark & ".csv";
  let benchmark_info := benchmark;
  let benchmark_start := _ampl_time;
  print "Benchmarking",benchmark_info,"...";
  
  if benchmark == "ARON-Shannon" then
  {
    print "Reseting FEC to Shannon limit";
    reset data FEC;
    for{p in P}
      let FEC[p] := M/(1-p);
  }
  
  for{tc in TC}
  {
    printf ",t=%s",
      round(tc,if tc>0 then max( round(-log(tc)/log(10))+1 ,3))
        > (file);
  }
  printf "\n" > (file);
  for{i in CHUNKS}
  {
    for{l in 1..min(10,maxlayers[i])}
    {
      printf "layer%d",l > (file);
      for{tc in TC}
      { 
        printf ",%f",
          if match(benchmark,"perlink")
            then aron_l[tc,i,l]
            else aron_f[tc,i,l] > (file);
      }
      printf "\n" > (file);
    }
    printf "\n" > (file);
  }
  close (file);
  
  let benchmark_stop := _ampl_time;
  reset data logrec;
  print logrec;
  print logrec >> (logfile);
  close (logfile);
}