#  AMPL Model File for the 
#  Practical Integrated Design Strategies 
#  for Opaque and All-Optical DWDM Networks: 
#  Optimization Models and Solution Procedures

#  option solver cplex;
#  option presolve 0;

set   LB24;          # Set of Max Spans
let   LB24 := {1..24};     
param LinkBudget {1..24};  
                     # Link Budget Table
let LinkBudget[1] := 162;
for {i in {02..09}} let LinkBudget[i] := LinkBudget[i-1] - 4;
let LinkBudget[10] := 128;
for {i in {11..24}} let LinkBudget[i] := LinkBudget[i-1] - 2;
param NoPaths;       # max no of paths per (o,d) pair
set N ordered;       # set of cities
set L within {N,N};  # set of links
set D within {N,N} default {};  
                     # set of demand pairs
param r {D};         # r[o,d] denotes the working traffic for 
                     # demand (o,d)               
param DPMD {L};      # Polarization-Mode Dispersion parameter
param TotalSpaces;   # the length of the Spacing Array
param costTE;        # the cost of TE
param costR;         # the cost of R
param costA;         # the cost of A
param costA01To20; param costA21To40; param costA41To80;
param costM01To20; param costM21To40; param costM41To80;
param ModA;          # modA denotes the number of wavelengths per fiber
param Hut {N};       # hut numbers for the cities/nodes
param NoDemands;     # number of demands
param LoDemand;      # lower limit on demand
param HiDemand;      # upper limit on demand
param RandomSeed;    # random seed for the problem generator
# display NoDemands; # number of demand pairs
# display RandomSeed;# random seed used to generate demand pairs
set T := {1..TotalSpaces};
                     # set of hut numbers
param Spaces {T};    # distances between huts/nodes
param Begin  {L};    # beginning hut number of a link
param End    {L};    # ending hut number of a link
param OpaqueCost;    # Opaque design cost
data dataEU.txt;
#display Spaces, Begin, End;
param length {L} default 0;   
                     # length[j,k] denotes the length of link (j,k)
param TimeInOpaque;
for {(j,k) in L} {
# printf"(j,k) = (%s,%s)  ",j,k;
  for {s in {Begin[j,k]..End[j,k]}} {
#  printf"Spaces[%d] = %d\n",s,Spaces[s];
   let length[j,k] := length[j,k] + Spaces[s];
  }
}
set Cities default {};
                     # nodes where the traffic originates and terminates
for {n in N} {
  let Cities := Cities union {Hut[n]};
}

set   TH;            # set of hut & city numbers
let   TH := T union Cities;

param TotalDemand default 0;
                     # Total demand generated

# Initialize Random Number Generator
param CardN;
let   CardN := card(N);
param n1;
param n2;
for {i in {1..RandomSeed} } {
  let n1 := Uniform(1,CardN);
}

# Generate The Demands
for {i in {1..NoDemands}} {
  for {j in {1..1000}} {
    let n1 := floor(Uniform(1,CardN));
    let n2 := floor(Uniform(1,CardN));
    if n1 != n2 and ((member(n1,N),member(n2,N)) not in D and (member(n2,N),member(n1,N)) not in D)then break;
  }
  if n1 == n2 then {
    printf "Bug in Generator\n\n";
    quit;
  }
let D := D union {(member(n1,N),member(n2,N))};
let r[member(n1,N),member(n2,N)] := floor(Uniform(LoDemand,HiDemand));
let TotalDemand := TotalDemand + r[member(n1,N),member(n2,N)];
}
display D,r;
#display TotalDemand;
param TimeToGenerate;
let TimeToGenerate := _ampl_user_time;

set   HL within {TH} ordered; # Huts in a given link
let   HL := {1};
param distL {TH};           # distances between huts in a given link
param ndist {TH} default 0; # ndist[th] denotes the dist to the next h 
                            # with equipment
param LinkB {TH} default 0;  
param cardHL;
param cardH;
param Budget;

var y {TH} >= 0 <= Budget;   
                  # y[h] denotes the flow from hut h to the sink
var Y {TH} binary; # Y[h] = 1 if y[h] > 0 and = 0 otherwise


#--------------------------------------------------------------------------------------------------------
# Determine HL and distL {HL}and solve the hut selection model for each link budget and for each link
param counter2;
param distLBig;
param BestCost;
param BestR {L};     # number of huts with amplifiers for link l using the best link budget
param BestA {L};     # number of huts with regeneration for link l using the best link budget
param BestM {L};     # number of huts/nodes with a MUX/DMUX for link l using the best link budget
param R  {TH};
param A  {TH};
param TE {TH};
param MaxK default 0;
param rdist {1..card(L)};
param c;
param dpmdLastR;
param KLastR;
param key;  
param loop6;
param SetR;          # SetR = 1 a regen is needed and
                     # SetR = 0 if not
param MaxDist;       # maximum distance before an o/e/o conversion
param TotalR;
param TotalA;
param TotalCost;
param Bestb {L};     # the best link budget for link l
param YY {TH};
param yy {TH};
param preH;
param curL;

for {(j,k) in L} {
  # display j,k;
  let HL := {Hut[j]};
  let distL[Hut[j]] := 0;
  let HL := HL union {Begin[j,k]..End[j,k]-1};
  for {h in HL diff {Hut[j]}}
    let distL[h] := Spaces[h];
  let HL := HL union {Hut[k]};
  let distL[Hut[k]] := Spaces[End[j,k]];
  # print huts and distances between the huts
  printf"HL\n";
  let counter2 := 0;
  for {h in HL} {
    printf"%6d ",h;
    let counter2 := counter2 + 1;
    if counter2 >= 10 then {
      printf"\n";
      let counter2 := 0;
    }
  }
  printf"\n\n";
  printf"distL\n";
  let counter2 := 0;
  let distLBig := 0.0;
  # determine the largest distance  
   for {h in HL} {
    printf"%6d ",distL[h];
    if distL[h] > distLBig then let distLBig := distL[h];
      let counter2 := counter2 + 1;
      if counter2 >= 10 then {
        printf"\n";
        let counter2 := 0;
      }
  }
  printf"\n\n";
  # solve the hut selection model for each link budget and for each link
  let BestCost     := 1.0e20;
  let BestR[j,k]   := 1.0e20;
  let BestA[j,k]   := 1.0e20;
  let BestM[j,k]   := 1.0e20;
  for {b in LB24} {
    let Budget := LinkBudget[b];
    if Budget < distLBig then {
      printf" Budget smaller than largest distance\n";
      break;
    }
    printf"Link Budget # %d   %d\n",b, LinkBudget[b];
    let preH:= 0;
    let curL:= 0;
    let {h in HL} y[h] := 0;
    let {h in HL} Y[h] := 0;
    for {h in HL} {
      if  curL + distL[h] <= Budget then let curL := curL + distL[h];
      else {
        let y[preH] := curL;
        let Y[preH] := 1;
        let curL:= distL[h];
      }
      if h = Hut[k] then {
        let y[h] := curL;
        let Y[h] := 1;
      }
      let preH:= h;
    }
    
    printf"Y\n";
    let counter2 := 0;
    for {h in HL} {
      printf"%6d ",Y[h];
      let counter2 := counter2 + 1;
      if counter2 >= 10 then {
        printf"\n";
        let counter2 := 0;
      }
    }
    printf"\n\n";
    let {h in HL} R[h]  := 0;
    let {h in HL} A[h]  := 0;
    let {h in HL} TE[h] := 0;

    let MaxK := 1;
    # Calculate ndist 
    let cardHL := card(HL);
    for {loop4 in {1..cardHL-1}} {
      for {loop5 in {loop4+1..cardHL}}
        if Y[member(loop5,HL)] == 1 then {
          let ndist[loop4] := y[member(loop5,HL)];
          break;
        }
      if Y[member(loop4,HL)] == 0 then let ndist[loop4] := 0;
    }  # for {loop4 in  }
    let ndist[cardHL] := 0;
    for {kk in {1..MaxK}} let rdist[kk] := 0;
    let c := 0;
    let dpmdLastR := DPMD[j,k];  
    let KLastR := 1;
    printf"loop   h dist   y Y ndist rdist K pmdLastR DPMD    R    A   TE\n";
    let loop6 := 0;
    for {h in HL} {
      let loop6 := loop6 + 1;
      let key   := 1;
      if h == first(HL) or h == last(HL) then {
        let TE[h] := 1;
        printf"%4d %3d %4d %3d %1d %5d %5d %1d",
          loop6,h,distL[h],y[h],Y[h],ndist[loop6],rdist[key],1;
        printf" %8.1f %4.1f %4d %4d %4d\n",
          dpmdLastR,dpmdLastR,R[h],A[h],TE[h];
        continue;
      }
      let rdist[key] := rdist[key] + distL[h];
      let SetR := 0;
      if Y[h] == 0 then {
        printf"%4d %3d %4d %3d %1d %5d %5d %1d",
        loop6,h,distL[h],y[h],Y[h],ndist[loop6],rdist[key],1;
        printf" %8.1f %4.1f %4d %4d %4d\n",
        dpmdLastR,dpmdLastR,R[h],A[h],TE[h];
        continue;
      }  # if Y[h]
      #Check Case 1:  Link Budget
      let c := c + 1;
      if c == b then {
        let SetR := 1;
        printf"------R Needed To Satisfy Link Budget----------\n";
      }
      #Check Case 2:  PMD Within A Link
      if SetR == 0 then {
        let MaxDist := 900/(DPMD[j,k]*DPMD[j,k]);
        if rdist[key]+ndist[loop6] > MaxDist then {
          let SetR := 1;
          printf"MaxDist: %f\n",MaxDist;
          printf"------R Needed To Satisfy PMD Within A Link----\n";
        }
      }

      if SetR == 1 then {
        let R[h]       := 1;
        let c          := 0;
        let rdist[key] := 0;
      }
      else 
        let A[h] := 1;
     
      printf"%4d %3d %4d %3d %1d %5d %5d %1d",
        loop6,h,distL[h],y[h],Y[h],ndist[loop6],rdist[key],1;
      printf" %8.1f %4.1f %4d %4d %4d\n",
        dpmdLastR,dpmdLastR,R[h],A[h],TE[h];
    }  # for {h in HL}

    let TotalR  := 0;
    let TotalA  := 0;
    for {th in HL} {
      let TotalR := TotalR + R[th];
      let TotalA := TotalA + A[th];
    }
    if BestR[j,k] = 0 and BestA[j,k] < (TotalA + 2 + 2 * TotalR) then break; 
    let TotalCost := 34*TotalR + TotalA; # 34x more expensive on average (1-80Lambdas)
    printf"\n\nLink Budget # %d, Value = %d\n\n",b, LinkBudget[b];
    if TotalCost <= BestCost then {
      let BestCost := TotalCost;
      let BestA[j,k]   := TotalA + 2 + 2 * TotalR;
      # 2 for the origin and destination & 2 for each R
      let BestR[j,k]   := TotalR;
      let BestM[j,k]   := 2 + 2 * TotalR;
      # 2 for the origin and destination & 2 for each R
      let Bestb[j,k]   := b;
      for {h in HL} {
        let YY[h] := Y[h];
        let yy[h] := y[h];
      }
    }
  }    # for {b in }
  # printf"\n(j,k,Bestb[j,k]): (%s,%s,%2d)\n",j,k,Bestb[j,k];
  # printf"\n(Best Link Budget): %d\n",LinkBudget[Bestb[j,k]];
}      # for {(j,k) in L}
param TimeForBestBudget;
let TimeForBestBudget := _ampl_user_time + _total_solve_user_time - TimeToGenerate;
#--------------------------------------------------------------------------------------------------------
# Determine candidate paths

set P default {};    # set of paths {1..no of paths};
set E within {N,N} default {};  
                     # set of arcs
set W {D} within P default {};      
                     # working paths for a (o,d)
set ArcsInW {P} within {E};
                     # denotes the set of links in working path p for (o,d) 
set Path {P} within {N} ordered;
                     # Path[p] denotes the nodes in the nth path
set PathsForARC {E} within P default {};
param w default 0;

for {(j,k) in L}
  let E := E union {(j,k),(k,j)};

param rhs {N} default 0;
                     # supply/demand for node i
var cx {E} binary;   # flow on arc (i,j)

minimize cost: sum{(i,j) in E} length[min(i,j),max(i,j)] * cx[i,j];

# Flow Out(i) - Flow In(i) = b(i)
subject to flow_balance {i in N}:
sum{j in N: (i,j) in E} cx[i,j] - sum{j in N: (j,i) in E} cx[j,i] = rhs[i];

subject to diff_path {(o,d) in D, p in W[o,d]}:
    sum {(i,j) in ArcsInW[p]} cx[i,j] <= card(ArcsInW[p]) -1;

#---Shortest Path Problem-----
    problem WorkingPath {(o,d) in D}:
#  objective
    cost,
#  variables
    cx,
#  constraints
    flow_balance, {p in W[o,d]} diff_path[o,d,p];


for {(o,d) in D} {
  for {spath in {1..NoPaths}} {
    # printf"loop     #%d  ",spath;
    # printf"Demand pair (%d,%d)\n",o,d;
    for {n in N} let rhs[n] := 0;
    let rhs[o] := +1;
    let rhs[d] := -1;
    problem WorkingPath[o,d];
    solve   WorkingPath[o,d];
    # expand;
    # printf"solve_result_num = %d\n",solve_result_num;
    if 200 <= solve_result_num && solve_result_num <= 299 then break;
    let w := w + 1;
    let P:= P union {w};
    # display o,d;
    let W[o,d] := W[o,d] union {w};
    let ArcsInW[w] := {};
    for {(i,j) in E} {
      if cx[i,j] > 0 then{
        let ArcsInW[w] := ArcsInW[w] union {(i,j)};
        let PathsForARC[i,j] := PathsForARC[i,j] union {w};
      }
    }
    let Path[w] := {o};
    for {(j,k) in ArcsInW[w]:j == o}{
      let Path[w] := Path[w] union {k};
      break;
    }
    repeat while card(Path[w]) < card(ArcsInW[w]) + 1 {    
      for {k in N: (last(Path[w]),k) in ArcsInW[w]}{
        let Path[w] := Path[w] union {k};
      }
    } # repeat     
  }   
} # for {(o,d) in D}
#display ArcsInW;
#display Path;
#display PathsForARC;
#--------------------------------------------------------------------------------------------------------
# Routing and Provisioning Model
param TimeForPaths;
let TimeForPaths := _ampl_user_time + _total_solve_user_time - TimeForBestBudget;

var F {P} >= 0 integer;     
                     # F[p] denotes the flow on path p
var TF {E} >= 0 integer;    
                     # TF[i,j] denotes the total flow on arc (i,j)
var A01_20 {L} >=0, <=1 integer; 
                     # denotes the number of A1-20s needed at a utilized hut on link l
var A21_40 {L} >=0, <=1 integer; 
                     # denotes the number of A21-40s needed at a utilized hut on link l
var A41_80 {L} >=0 integer; 
                     # denotes the number of A41-80s needed at a utilized hut on link l
var TErE {N} >=0;    # denotes the number of TEs needed at node n
var RE {L} >=0;      # denotes the number of Rs needed on link l

subject to Demand {(o,d) in D}:
    sum {p in W[o,d]} F[p] = r[o,d];

subject to TFinArc {(i,j) in E}:
    TF[i,j] = sum {p in PathsForARC[i,j]} F[p];

subject to AsonL {(i,j) in L}:
    TF[i,j] + TF[j,i] <= 20 * A01_20[i,j] + 40 * A21_40[i,j] + 80 * A41_80[i,j];
    
subject to SmallAs {(i,j) in L}:
    A01_20[i,j] + A21_40[i,j] <= 1;

subject to TEsonN {k in N}:
    sum {(k,j) in E} (TF[k,j] + TF[j,k]) = TErE[k];
    
subject to RsonL {(i,j) in L}:
    (TF[i,j] + TF[j,i]) * BestR[i,j] = RE[i,j];

minimize DesignCost: sum {k in N} (costTE * TErE[k]) 
+ sum {(i,j) in L} ((BestA[i,j] * (costA01To20 * A01_20[i,j] + costA21To40 * A21_40[i,j] + costA41To80 * A41_80[i,j])) 
+ costR * RE[i,j])
+ sum {(i,j) in L} ((costM01To20 * A01_20[i,j] + costM21To40 * A21_40[i,j] + costM41To80 * A41_80[i,j]) * BestM[i,j]);

#---Routing Problem-----
    problem Routing:
#  objective
    DesignCost,
#  variables
    F, TF, A01_20, A21_40, A41_80, TErE, RE,
#  constraints
    Demand, TFinArc, AsonL, SmallAs, TEsonN, RsonL;

problem Routing;
option cplex_options 'timing=1 mipdisplay=2 mipgap=0.010 display=2';
option show_stats 1;
solve Routing;
param TimeForRouting;
param TradA;
param TradR;
param TradTE;
param TradM;
let TimeForRouting := _solve_user_time;
let TradA := sum {(i,j) in L} (BestA[i,j] * (A01_20[i,j] + A21_40[i,j] + A41_80[i,j]));
let TradR := sum {(i,j) in L} RE[i,j];
let TradTE := sum {k in N} TErE[k];
let TradM := sum {(i,j) in L} (BestM[i,j] * (A01_20[i,j] + A21_40[i,j] + A41_80[i,j]));
let OpaqueCost := DesignCost;
printf"Opaque Cost       = %4d\n", OpaqueCost;
printf"Opaque AMPL Time  = %4.2lf\n", _ampl_user_time;
printf"Opaque CPLEX Time = %4.2lf\n\n", _total_solve_user_time;
let TimeInOpaque := _ampl_user_time + _total_solve_user_time;

param AVEnHOPS default 0;
param PathsInUse default 0;
# for {(o,d) in D} {
    # for {i in W[o,d] : F[i] > 0} {
        # let PathsInUse := PathsInUse + 1;
        # let AVEnHOPS := AVEnHOPS + card(Path[i]) - 1;
    # }
# }
# let AVEnHOPS := AVEnHOPS / PathsInUse;
# printf"Average Number of Hops = %4.2lf\n\n", AVEnHOPS;
#--------------------------------------------------------------------------------------------------------
# All-Optical DWDM NETWORK DESIGN
printf "\n*** All-Optical Design Starts ***\n\n";
param AllOpticalCost default 0;
param TimeInAllOptical;
param ysum;
set H {(o,d) in D,W[o,d]} within TH ordered;      
                     # H[o,d] denotes the huts in the path n from o to d
param K    {(o,d) in D, n in W[o,d], H[o,d,n]};  
                     # K[o,d,n,h] denotes the link number in path n from o to d
param dist {(o,d) in D, n in W[o,d], H[o,d,n]};  
                     # dist[o,d,n,h] denotes the dist from the previous hut to hut h
                     # for demand (o,d) path n
param dpmd {(o,d) in D, n in W[o,d], {1..card(L)}};     
                     # dpmd[o,d,n,k] denotes the DPMD value for link number k in path n from o to d
set   Link {(o,d) in D, n in W[o,d], {1..card(L)}} within L;     
                     # Link[o,d,s] denotes the link in segment s from o to d path n
                                 
set AA within N;
set BB within N;
set CC within {N,N};
set DD within {N,N};
set Temp1 ordered;
set Temp2 ordered;
param Max {L} default 0;
param lasthut;
param lastdist;
param index;
# param counter2;
param bbb;
param LinkV {1..card(L)};
param MinLB;
param dpmdv {1..card(L)};
param Value1; param Value2; param Value3;
param Value4; param Value5; param Value6;
param TotalA01To20 default 0; param TotalA21To40 default 0;
param TotalA41To80 default 0; param TotalM01To20 default 0;
param TotalM21To40 default 0; param TotalM41To80 default 0;
param TotalFiber default 0; 
                     # Fiber length
param TotalTE default 0;
param FF; param FFF;
param RR; param RRR;
param TT; param TTT;
param TTTT; param TTTTT;
param CostReduction;        
let {h in Cities} YY[h] := 1;
let {th in TH} R[th]  := 0;
let {th in TH} A[th]  := 0;
let {th in TH} TE[th] := 0;
                               
for {(o,d) in D} {
  for {i in W[o,d] : F[i] > 0} {
    let PathsInUse := PathsInUse + 1;
    let AVEnHOPS := AVEnHOPS + card(Path[i]) - 1;
    # determine H[o,d] and dist[o,d,h]
    let H[o,d,i] := {};
    for {loop in {2..card(Path[i])}} {
      let AA  := {member(loop-1,Path[i])};
      let BB  := {member(loop  ,Path[i])};
      let CC  := {AA,BB};
      let DD  := {BB,AA};
      let key := 0;
      if CC within L then let key := 1;
      if DD within L then let key := 2;
      # display AA,BB,CC,DD,key;
      if key == 1 then {
        for {(j,k) in CC} {
          let Temp1    := {Begin[j,k]..End[j,k]-1};
          let lasthut  := Hut[k];
          let lastdist := Spaces[End[j,k]];
          let dpmd[o,d,i,loop-1] := DPMD[j,k];
          let Link[o,d,i,loop-1] := CC;
          let Max[j,k] := Max[j,k] + F[i];
        }
      }
      if key == 2 then {
        for {(j,k) in DD} {
          let Temp1    := {Begin[j,k]..End[j,k]-1};
          let lasthut  := Hut[j];
          let lastdist := Spaces[Begin[j,k]];
          let dpmd[o,d,i,loop-1] := DPMD[j,k];
          let Link[o,d,i,loop-1] := DD;
          let Max[j,k] := Max[j,k] + F[i];
        }
      }
      # display Temp1, lasthut, lastdist;
      let Temp2 := {1..card(Temp1)};
      # display Temp2;
      if loop == 2 then {
        # display H[o,d,i];
        let H[o,d,i] := {Hut[o]};
        let dist[o,d,i,Hut[o]] := 0;
        # display dist[o,d,i,Hut[o]];
        let K [o,d,i,Hut[o]] := loop - 1;
      }
      for {loop2 in Temp2} {
        if key == 1 then let index := loop2;
        if key == 2 then let index := card(Temp1)-loop2+1;
        let H[o,d,i] := H[o,d,i] union {member(index,Temp1)};
        if key == 1 then 
          let dist[o,d,i,member(index,Temp1)] := Spaces[member(index,Temp1)];
        if key == 2 then 
          let dist[o,d,i,member(index,Temp1)] := Spaces[member(index,Temp1)+1];
        let K[o,d,i,member(index,Temp1)] := loop - 1;
      }
      let H[o,d,i] := H[o,d,i] union {lasthut};
      let dist[o,d,i,lasthut] := lastdist;
      let K   [o,d,i,lasthut] := loop - 1;
    }  # for {loop in   }
    printf"H[o,d,i]\n";
    let counter2 := 0;
    for {h in H[o,d,i]} {
      printf"%6d ",h;
      let counter2 := counter2 + 1;
      if counter2 >= 10 then {
        printf"\n";
        let counter2 := 0;
      }
    }
    printf"\n\n";
    printf"dist[o,d,i,*]\n";
    let counter2 := 0;
    for {h in H[o,d,i]} {
      printf"%6d ",dist[o,d,i,h];
      let counter2 := counter2 + 1;
      if counter2 >= 10 then {
        printf"\n";
        let counter2 := 0;
      }
    }
    printf"\n\n";
    printf"   K[o,d,i,*]\n";
    let counter2 := 0;
    for {h in H[o,d,i]} {
      printf"%6d ",   K[o,d,i,h];
      let counter2 := counter2 + 1;
      if counter2 >= 10 then {
        printf"\n";
        let counter2 := 0;
      }
    }
    printf"\n\n";
    printf"\n\n"; 
  
    let MaxK := 1;
    let ysum := 0;
    for {h in H[o,d,i]} {
      if K[o,d,i,h] > 1 then let MaxK := K[o,d,i,h];
      let ysum   := ysum + dist[o,d,i,h];
      let yy[h]  := ysum;
      if YY[h] == 1 then let ysum := 0;
    } 
    for {kk in {1..MaxK}} {
      for {(jjj,kkk) in Link[o,d,i,kk]} {
        let bbb := Bestb[jjj,kkk];
        let LinkV[kk] := bbb;
        printf"(Link, Bestb, Link Budget): (%s %s %3d %3d) \n",
          jjj, kkk, bbb, LinkBudget[bbb];
      }
    }
    # Calculate ndist 
    let cardH := card(H[o,d,i]);
    for {loop4 in {1..cardH-1}} {
      for {loop5 in {loop4+1..cardH}}
        if YY[member(loop5,H[o,d,i])] == 1 then {
          let ndist[loop4] := yy[member(loop5,H[o,d,i])];
          break;
        }
        if YY[member(loop4,H[o,d,i])] == 0 then let ndist[loop4] := 0;
    }  # for {loop4 in  }
    let ndist[cardH] := 0;
    let loop6 := 0;
    for {h in H[o,d,i]} {
      let loop6        := loop6 + 1;
      let key          := K[o,d,i,h];
      let LinkB[loop6] := LinkV[key];
    }
    for {kk in {1..MaxK}} let rdist[kk] := 0; 
    let c := 0;
    let dpmdLastR := dpmd[o,d,i,1];
    let KLastR := 1;
    printf"loop   h dist   y Y ndist rdist K pmdLastR DPMD    R    A   TE Link Budget\n";
    let loop6 := 0;
    let MinLB := LinkB[1];
    for {h in H[o,d,i]} {
      let loop6 := loop6 + 1;
      let key   := K[o,d,i,h];
      if h == first(H[o,d,i]) or h == last(H[o,d,i]) then {
        let TE[h] := TE[h] + F[i];
        printf"%4d %3d %4d %3d %1d %5d %5d %1d",
          loop6,h,dist[o,d,i,h],yy[h],YY[h],ndist[loop6],rdist[key],K[o,d,i,h];
        printf" %8.1f %4.1f %4d %4d %4d %11d\n",
          dpmdLastR,dpmd[o,d,i,key],R[h],A[h],TE[h],LinkB[loop6];
        continue;
      }        
	  let rdist[key] := rdist[key] + dist[o,d,i,h];
	  let dpmdv[key] := dpmd[o,d,i,key];
	  let MinLB := min(MinLB, LinkB[loop6]);
	  let SetR := 0;
	  if YY[h] == 0 then {
        printf"%4d %3d %4d %3d %1d %5d %5d %1d",
          loop6,h,dist[o,d,i,h],yy[h],YY[h],ndist[loop6],rdist[key],K[o,d,i,h];
        printf" %8.1f %4.1f %4d %4d %4d %11d\n",
          dpmdLastR,dpmd[o,d,i,key],R[h],A[h],TE[h],LinkB[loop6];
        continue;
      }
      # Check Case 1:  Link Budget
      let c := c + 1;
      if c >= MinLB or c >= LinkB[loop6+1] then {
        let SetR := 1;
        printf"------R Needed To Satisfy Link Budget----------\n";
      }
	  # Check Case 2:  PMD Within A Link
	  if SetR == 0 && KLastR == key then {
	    let MaxDist := 900/(dpmd[o,d,i,key]*dpmd[o,d,i,key]);
	    if rdist[key]+ndist[loop6] > MaxDist then {
	      let SetR := 1;
	      printf"MaxDist: %f\n",MaxDist;
	      printf"------R Needed To Satisfy PMD Within A Link----\n";
	    }
	  }
	  # Check Case 3: PMD Across A City      
	  if SetR == 0 && KLastR+1 == key then {
	    let Value1 := dpmdLastR*dpmdLastR*rdist[KLastR];
	    let Value2 := dpmd[o,d,i,key]*dpmd[o,d,i,key]*(rdist[key]+ndist[loop6]);
	    let Value3 := Value1 + Value2;
	    if  Value3 >= 900 then {
	      let SetR := 1;
	      printf"Value1, Value2, Value3: %f %f %f\n",Value1, Value2, Value3;
	      printf"------R Needed To Satisfy PMD Across A City----\n";
	    }
	  }
	  # Check Case 4: PMD Across 2 Cities   
	  if SetR == 0 && KLastR+2 == key then {
	    let Value1 := dpmdv[KLastR]*dpmdv[KLastR]*rdist[KLastR];
	    let Value2 := dpmdv[KLastR+1]*dpmdv[KLastR+1]*rdist[KLastR+1];
	    let Value3 := dpmd[o,d,i,key]*dpmd[o,d,i,key]*(rdist[key]+ndist[loop6]);
	    let Value4 := Value1 + Value2 + Value3;
	    if  Value4 >= 900 then {
	      let SetR := 1;
	      printf"Value1, Value2, Value3, Value4: %f %f %f %f\n",Value1, Value2, Value3, Value4;
	      printf"------R Needed To Satisfy PMD Across 2 Cities----\n";
	    }
	  }
	  # Check Case 5: PMD Across 3 Cities   
	  if SetR == 0 && KLastR+3 == key then {
	    let Value1 := dpmdv[KLastR]*dpmdv[KLastR]*rdist[KLastR];
	    let Value2 := dpmdv[KLastR+1]*dpmdv[KLastR+1]*rdist[KLastR+1];
	    let Value3 := dpmdv[KLastR+2]*dpmdv[KLastR+2]*rdist[KLastR+2];
	    let Value4 := dpmd[o,d,i,key]*dpmd[o,d,i,key]*(rdist[key]+ndist[loop6]);
	    let Value5 := Value1 + Value2 + Value3 + Value4;
	    if  Value5 >= 900 then {
	      let SetR := 1;
	      printf"Value1, Value2, Value3, Value4, Value5: %f %f %f %f %f\n",Value1, Value2, Value3, Value4, Value5;
	      printf"------R Needed To Satisfy PMD Across 3 Cities----\n";
	    }
	  }
	  # Check Case 6: PMD Across 4 Cities   
	  if SetR == 0 && KLastR+4 == key then {
	    let Value1 := dpmdv[KLastR]*dpmdv[KLastR]*rdist[KLastR];
	    let Value2 := dpmdv[KLastR+1]*dpmdv[KLastR+1]*rdist[KLastR+1];
	    let Value3 := dpmdv[KLastR+2]*dpmdv[KLastR+2]*rdist[KLastR+2];
	    let Value4 := dpmdv[KLastR+3]*dpmdv[KLastR+3]*rdist[KLastR+3];
	    let Value5 := dpmd[o,d,i,key]*dpmd[o,d,i,key]*(rdist[key]+ndist[loop6]);
	    let Value6 := Value1 + Value2 + Value3 + Value4 + Value5;
	    if  Value6 >= 900 then {
	      let SetR := 1;
	      printf"Value1, Value2, Value3, Value4, Value5, Value6: %f %f %f %f %f %f\n",Value1, Value2, Value3, Value4, Value5, Value6;
	      printf"------R Needed To Satisfy PMD Across 4 Cities----\n";
	    }
	  }
	  if KLastR+4 < key then {
	    printf"$$$$$$$$$$$$ Bug #1 $$$$$$$$$$$$$$\n\n";
	    quit;
	  } 
	  #display h, YY[h], SetR;
	  if SetR == 1 then {
	    let R[h]       := R[h] + F[i];
	    let c          := 0;
	    let dpmdLastR  := dpmd[o,d,i,key];
	    let KLastR     := key;
	    let rdist[key] := 0;
	    let MinLB      := LinkB[loop6+1];
	  }
	  if SetR == 0  && h not in Cities then {
	    let A[h] := A[h] + F[i];
	  }
      printf"%4d %3d %4d %3d %1d %5d %5d %1d",
        loop6,h,dist[o,d,i,h],yy[h],YY[h],ndist[loop6],rdist[key],K[o,d,i,h];
      printf" %8.1f %4.1f %4d %4d %4d %11d\n",
        dpmdLastR,dpmd[o,d,i,key],R[h],A[h],TE[h],LinkB[loop6];
    } # for {h in H  }
  }
}
param TimeInHdist;
let TimeInHdist :=  _ampl_user_time + _total_solve_user_time - TimeInOpaque;
#set TL within {N,N}; # set of links to analyze for equipment placement     

# Determine All-Optical COST

for {(j,k) in L} {
  # display j, k, Max[j,k];
  if Max[j,k] == 0 then continue;
  let FF  := floor(Max[j,k]/ModA);
  let FFF := Max[j,k] - FF *ModA;
  let TotalA41To80 := TotalA41To80 + 2*FF;
  let TotalM41To80 := TotalM41To80 + 2*FF;
  if 01 <= FFF && FFF <= 20 then {
    let TotalA01To20 := TotalA01To20 + 2;
    let TotalM01To20 := TotalM01To20 + 2;
  }
  if 21 <= FFF && FFF <= 40 then {
    let TotalA21To40 := TotalA21To40 + 2;
    let TotalM21To40 := TotalM21To40 + 2;
  }
  if 41 <= FFF && FFF <= 80 then {
    let TotalA41To80 := TotalA41To80 + 2;
    let TotalM41To80 := TotalM41To80 + 2;
  }
} # for {(j,k) in L }  

# Huts
let TotalTE := 0;
let TotalR := 0;
for {th in TH} {
  if TE[th] == 0 and A[th] == 0 and R[th] == 0 then continue;
  # display th, TE[th], A[th], R[th];
  let TotalTE := TotalTE + TE[th];
  let TotalR := TotalR + R[th];
  if th in Cities then continue;
  let RR   := floor(R[th]/ModA);
  let RRR  := R[th] - RR*ModA;
  let TTTT := A[th];
  let TotalA41To80 := TotalA41To80 + 2*RR;
  let TotalM41To80 := TotalM41To80 + 2*RR;
  if RRR > 0 then {
    let TTT  := min(TTTT,80-RRR);
    let RRR  := RRR + TTT;
    let TTTT := TTTT - TTT;
  }
  # Case 1: 01 <= RRR <= 20 => 
  #      update 01-20, TTTT := max(0,TTTT-(20-RRR))
  if 01 <= RRR && RRR <= 20 then {
    let TotalA01To20 := TotalA01To20 + 2;
    let TotalM01To20 := TotalM01To20 + 2;
    let TTTT := max(0,TTTT-(20-RRR));
  }  # Case 1
  # Case 2: 21 <= RRR <= 40 => 
  #      update 21-40, TTTT := max(0,TTTT-(40-RRR))
  if 21 <= RRR && RRR <= 40 then {
    let TotalA21To40 := TotalA21To40 + 2;
    let TotalM21To40 := TotalM21To40 + 2;
    let TTTT := max(0,TTTT-(40-RRR));
  }  # Case 2
  # Case 3: 41 <= RRR <= 80 => 
  #      update 41-80, TTTT := max(0,TTTT-(80-RRR))
  if 41 <= RRR && RRR <= 80 then {
    let TotalA41To80 := TotalA41To80 + 2;
    let TotalM41To80 := TotalM41To80 + 2;
    let TTTT := max(0,TTTT-(80-RRR));
  }  # Case 3
  let TT  := floor(TTTT/ModA);
  let TTT := TTTT - TT*ModA;
  let TotalA41To80 := TotalA41To80 + TT;
  if 01 <= TTT && TTT <= 20 then
    let TotalA01To20 := TotalA01To20 + 1;
  if 21 <= TTT && TTT <= 40 then
    let TotalA21To40 := TotalA21To40 + 1;
  if 41 <= TTT && TTT <= 80 then
    let TotalA41To80 := TotalA41To80 + 1;
}  # for {th in TH}  #135

let AllOpticalCost := costTE*TotalTE + costR*TotalR 
  + costA01To20*TotalA01To20 + costM01To20*TotalM01To20
  + costA21To40*TotalA21To40 + costM21To40*TotalM21To40
  + costA41To80*TotalA41To80 + costM41To80*TotalM41To80;
if AllOpticalCost == 0 then let AllOpticalCost := OpaqueCost;
display AllOpticalCost,OpaqueCost;
let TotalA := TotalA01To20 + TotalA21To40 + TotalA41To80;
display TotalTE, TotalR, TotalA;
display TotalA01To20, TotalA21To40, TotalA41To80;
display TotalM01To20, TotalM21To40, TotalM41To80;
let CostReduction := 100.0*(OpaqueCost-AllOpticalCost)/AllOpticalCost;
printf"Cost Reduction = %4.2lf\n\n",CostReduction;
display TotalDemand;
#display _ampl_user_time;
#display _total_solve_user_time;
let TimeInAllOptical := (_ampl_user_time + _total_solve_user_time) - TimeInOpaque;
printf"Time to Generate         =   %4.2lf\n", TimeToGenerate;
printf"Time for Best Budget     =   %4.2lf\n", TimeForBestBudget;
printf"Time for Paths           =   %4.2lf\n", TimeForPaths;
printf"Time for Routing(CPLEX)  =   %4.2lf\n", TimeForRouting;
printf"Time in Opaque           =   %4.2lf\n", TimeInOpaque;
printf"Time in Hdist(Part1)     =   %4.2lf\n", TimeInHdist;
printf"Time In AllOptical       =   %4.2lf\n", TimeInAllOptical;
printf"Total Time = %4.2lf\n\n", (_ampl_user_time + _total_solve_user_time);
let AVEnHOPS := AVEnHOPS / PathsInUse;
printf"Average Number of Hops   = %4.2lf\n\n", AVEnHOPS;
printf"Opaque Design As         = %4d\n", TradA;
printf"Opaque Design Rs         = %4d\n", TradR;
printf"Opaque Design TEs        = %4d\n", TradTE;
printf"Opaque Design MUX/DMUXs  = %4d\n", TradM;