#   All-Optical Design with Leasing Enabled
#   Table 9
#   option presolve 0;
#   option solver cplex;

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;

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} default 0.1;      
                     # 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
param NoPaths;       # max number of cand. paths per demand
# 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
param QFib;          # unavailability of a fiber cable/duct
param QReg;          # unavailability of a regen
param QTe;           # unavailability of a TE
param QMUX;          # unavailability of a MUX
param QOXC;          # unavailability of a OXC
param QAMP;          # unavailability of a Amp
param runNo default 0;
param cc {N,N} default 0;
set LG within {N,N} default {};  
                     # set of links
set EG within {N,N} default {};  
                     # set of links
data DataEL.txt;

for {i in N}{
  for {j in N:(min(i,j),max(i,j)) not in L and i<>j}{
    let LG:=LG union {(min(i,j),max(i,j))};
    let EG:= EG union {(i,j), (j,i)};
  }
}
let L:=L union LG;
display LG;
#display Spaces, Begin, End;
param length {L} default 0;   
                     # length[j,k] denotes the length of link (j,k)
param TimeInOpaque{1..runNo} default 0;
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];
  }
}
param totallength default 0;
let {(j,k) in L} totallength := totallength + length[j,k];
let totallength:= totallength / card(L);
#display card(L),totallength;

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 := round(Uniform(1,CardN));
    let n2 := round(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


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} default 100;  
                     # the best link budget for link l
param YY {TH};
param yy {TH};
param preH;
param curL;
for {(j,k) in LG}{
  let BestA[j,k]:= ceil(Spaces[Begin[j,k]]/LinkBudget[24]);
  let BestR[j,k]:= floor(BestA[j,k]/24);
  let BestM[j,k]:= 2 + 2 * BestR[j,k];
}
# display BestA, BestR, BestM;

for {(j,k) in L: (j,k) not in LG} {
  #display j,k;
  # determine HL and distL {HL}
  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;
  # also 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];
    # *** Hut Selection Begins ***
    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;
    }
    # *** Hut Selection Ends ***
    #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)
    #display TotalCost;
    #display TotalA, TotalR;
    #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;

#------------------------------------------------------------

set P default {};    # set of paths {1..no of paths};
param pm {D} default 2;
                     # denotes the no of individual paths for path p 
param pn {P} default 2; 
                     # denotes the no of individual paths for path p
set E within {N,N} default {};  
                     # set of arcs
set PD {D} within P default {};   
                     # working & protection paths for an (o,d)          
set ArcsInPath {pk in P,1..pn[pk]} within {E};
                     # new denotes the set of links in 
                     # path p,n for (o,d)           
set ArcsWP {P} within {E};

set MPath {pk in P,1..pn[pk]} within {N} ordered;
                     # WPath[p] denotes the nodes in the nth working 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)

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 TEleased {N} >=0;# denotes the number of TEs leased at node n
var RE {L} >=0;      # denotes the number of Rs needed on link l
var PO{D,P} binary;  # PO[o,d,p] is 1 if the path p for demand (o,d) is used and 0 otherwise

param TimeForPaths{1..runNo} default 0;
param TimeForRouting{1..runNo} default 0;
param TradA;
param TradR;
param TradTE;
param AVEnHOPS default 0;
param PathsInUse default 0;
param UN{D} default 0;
param UNw;
param zP;

# for All-Optical DWDM NETWORK DESIGN
param AllOpticalCost default 0;
param TimeInAllOptical{1..runNo} default 0;
param ysum;
set H {(o,d) in D,{1..pm[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 {1..pm[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 {1..pm[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 {1..pm[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 {1..pm[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 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;        


#UN for All-Optical
set TL within {N,N}; # set of links to analyze for equipment placement
param svR  {TH};
param svA  {TH};
param svTE {TH};
param svMax {L} default 0;
set Path within {N} ordered;
param newpaths;
param trsh default 9.50643E-06; 
                     # denotes the threshold value for unavailability
set UNt within {N,N} default {}; 
                     # set of demand pairs with unavailabilities higher than trsh

subject to PathOperational {(o,d) in D}:
  sum {p in PD[o,d]} PO[o,d,p] = 1;

subject to TFinArc {(i,j) in E}:
  TF[i,j] = sum {p in PathsForARC[i,j], (o,d) in D : p in PD[o,d]} r[o,d] * PO[o,d,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 TEsLeased {k in N}:
  sum {(k,j) in E:(j,k) in LG or (k,j) in LG} (TF[k,j] + TF[j,k]) = TEleased[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] + 100 * TEleased[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]) + 100 * sum {(i,j) in LG} ((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])
  + 100 * sum {(i,j) in LG} ((costM01To20 * A01_20[i,j] + costM21To40 * A21_40[i,j] + costM41To80 * A41_80[i,j]) * BestM[i,j]);

minimize cost: sum{(i,j) in E} (length[min(i,j),max(i,j)] * cx[i,j]) + sum{(i,j) in EG} (100000 * 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 nodelimit_in {j in N: rhs[j] = 0}: sum{i in N:(i,j) in E} cx[i,j] <= 1;

subject to nodelimit_out {j in N: rhs[j] = 0}: sum{i in N:(j,i) in E} cx[j,i] <= 1;

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


#---Shortest Path Problem-----
  problem WorkingPath {(o,d) in D}:
#  objective
    cost,
#  variables
    cx,
#  constraints
    flow_balance, {p in PD[o,d]} diff_path[o,d,p], nodelimit_in, nodelimit_out;
    
#---Routing Problem-----
  problem Routing:
#  objective
    DesignCost,
#  variables
    TF, PO, A01_20, A21_40, A41_80, TErE, TEleased, RE,
#  constraints
   TFinArc, PathOperational, AsonL, SmallAs, TEsonN, TEsLeased, RsonL;
   
show problems;
display _PROBS;
param TotalTime default 0;
param UNmax;
param cardLeased default 0;
# for processing leased
param chkL;
param chkL2;
param LsdTE;
set LEASEDL ordered;
set LEASEDN;
set LEASEDLINK dimen 2;
param LeasedDistance;
param ALeased;
param BLeased;
param CLeased;
param HLeased;
param GLeased;
param FLeased;
param TotalOXC;
param LeasedLineCost;
param LeasedA20;
param LeasedA40;
param LeasedA80;
param LeasedM20;
param LeasedM40;
param LeasedM80;
param TotalLeasedCost;
param LCdistance default 0;
param LeasedR;
param LeasedRi;
param LeasedTE;
param KLeased;
param nodeL default 0;
# end of processing leased
param nonewp default 0;

repeat{ 
  let TotalLeasedCost := 0;
  let runNo:= runNo + 1;
  let UNmax:= 0;
  printf" \n -------------------------------- \n ";
  printf" \n ******* Run Number : %5d ******* \n ", runNo;
  printf" \n -------------------------------- \n ";
  for {(o,d) in D} {
    let newpaths:= w;
    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] := pm[o,d];
      let rhs[d] := -1*pm[o,d];
      #display P, PD, rhs, ArcsInPath, ArcsWP, MPath, PathsForARC, pn;
      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};
      let pn[w] := pm[o,d];
      # display o,d;
      let PD[o,d] := PD[o,d] union {w};
      let ArcsWP[w] := {};
      let {i in 1..pm[o,d]} ArcsInPath[w,i] := {};
      for {(i,j) in E}
        if cx[i,j] > 0 then let ArcsWP[w] := ArcsWP[w] union {(i,j)};
      let {i in 1..pm[o,d]} MPath[w,i] := {o};
      for {(j,k) in ArcsWP[w]:j = o}{
        for {i in 1..pm[o,d]}{
          if card(MPath[w,i]) = 1 then {
            let MPath[w,i] := MPath[w,i] union {k};
            let PathsForARC[j,k] := PathsForARC[j,k] union {w};
            let ArcsInPath[w,i] := ArcsInPath[w,i] union {(j,k)};
            break;
          }
          else continue;
        }
      }
      repeat while (sum {i in 1..pm[o,d]} card(MPath[w,i])) < card(ArcsWP[w]) + 2 {
        for {j in 1..pm[o,d]}{
          for {k in N: (last(MPath[w,j]),k) in ArcsWP[w]}{
              let ArcsInPath[w,j] := ArcsInPath[w,j] union {(last(MPath[w,j]),k)};
              let PathsForARC[last(MPath[w,j]),k] := PathsForARC[last(MPath[w,j]),k] union {w};
              let MPath[w,j] := MPath[w,j] union {k};
          }
        }   
      }
    } 
    if newpaths = w then quit;
  } # for {(o,d) in D}
  #  display ArcsWP, MPath;
  # display PathsForARC, ArcsInPath;
  if runNo = 1 then let TimeForPaths[1] := _ampl_user_time - TimeForBestBudget;
  else let TimeForPaths[runNo] := _ampl_user_time - TotalTime;
  # determine Routing for Opaque Design
  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;
  problem Routing;
  option cplex_options 'timing=1 mipdisplay=2 mipgap=0.050 display=2'; 
  option show_stats 1;
  problem Routing;
  solve Routing;
  # display PO, TF;
  let TimeForRouting[runNo] := _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 OpaqueCost := sum {k in N} (costTE * (TErE[k] + 0.25 * TEleased[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]) + 0.25 * sum {(i,j) in LG} ((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])
  + 0.25 * sum {(i,j) in LG} ((costM01To20 * A01_20[i,j] + costM21To40 * A21_40[i,j] + costM41To80 * A41_80[i,j]) * BestM[i,j]);
  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[runNo] := _ampl_user_time - TotalTime;
  display TimeInOpaque,TimeForRouting;
  # Unavailability for Opaque
  for {(o,d) in D}{
    for {z in PD[o,d] : PO[o,d,z] = 1} let zP:= z;
    let UNw := 0;
    for {k in 1..pn[zP]} {
      let UNw := 0;
      for {(i,j) in ArcsInPath[zP,k]}
        let UNw := UNw + ((ceil(r[o,d]/80)) * (BestA[min(i,j),max(i,j)] * QAMP + BestM[min(i,j),max(i,j)] * QMUX)) + r[o,d] * BestR[min(i,j),max(i,j)] * QReg + length[min(i,j),max(i,j)] * (QFib);
      let UNw:= UNw + (card(MPath[zP,k]) - 2) * (QOXC + 2 *  QTe * r[o,d]) + 2 * QTe * r[o,d];
      #display (card(MPath[zP,k]) - 2);
      if UN[o,d] = 0  then  let UN[o,d] := UNw;
      else let UN[o,d] := UN[o,d] * UNw;
    }
    let UN[o,d] := UN[o,d] + 2 * QOXC;  
  }
  # display UN;
  printf"Links to be leased:\n";
  printf"From   To   Distance   R\n";
  for {(j,k) in L: (min(j,k),max(j,k)) in LG} {
    if TF[j,k] > 0 or TF[k,j] > 0 then{
      printf"%3d   %3d    %4d     %3d\n",min(j,k),max(j,k),cc[min(j,k),max(j,k)],TF[j,k]+TF[k,j]; 
      let cardLeased:= cardLeased + 1;     
    }
  }
  for {(o,d) in D} if UNmax < UN[o,d] then let UNmax:= UN[o,d];
  let UNmax := UNmax * 525960;
  printf" \nMax Unavailability for All-Optical (min): %.8f\n\n ", UNmax;
  display cardLeased;
  let cardLeased:= 0;
  let UNmax:= 0;
  let {(o,d) in D} UN[o,d]:= 0;
  # Number of Hops for Opaque
  # for {(o,d) in D} {
  # for {p in PD[o,d] : PO[o,d,p] = 1} {
  #   for {i in 1..pm[o,d]} {
  #     let PathsInUse := PathsInUse + 1;
  #     let AVEnHOPS := AVEnHOPS + card(MPath[p,i]) - 1;
  #   }
  # }
  # }
  # let AVEnHOPS := AVEnHOPS / PathsInUse;
  # printf"Average Number of Hops = %4.2lf\n\n", AVEnHOPS;
  let AVEnHOPS := 0; let PathsInUse := 0;
  
  # AllOptical Starts
  printf "\n*** All-Optical Design Starts ***\n\n";
  let {th in TH} svR[th]  := 0;
  let {th in TH} svA[th]  := 0;
  let {th in TH} svTE[th] := 0;
  for {(o,d) in D} {
    let ALeased:= floor(r[o,d]/ModA);
    let BLeased:= r[o,d] mod ModA;
    let LeasedLineCost:=0;
    let LCdistance:= 0;
    for {i in {1..pm[o,d]}}{
      #printf"\n\n";
      #display o,d;
      let TL:={};
      let LeasedDistance:= 0;
      let LEASEDL:= {};
      let LEASEDN:= {};
      let LEASEDLINK:= {};
      let HLeased:= 0;
      let GLeased:= 0;
      let KLeased:= 0;
      let TotalOXC:= 0;
      let LeasedA20:=0;
      let LeasedA40:=0;
      let LeasedA80:=0;
      let LeasedM20:=0;
      let LeasedM40:=0;
      let LeasedM80:=0;
      let LeasedR:=0;
      let LeasedRi:=0;
      let LeasedTE:=0;
      for {z in PD[o,d] : PO[o,d,z] = 1} let Path := MPath[z,i];      
      let PathsInUse := PathsInUse + 1;
      let AVEnHOPS := AVEnHOPS + card(Path) - 1;
      # determine H[o,d] and dist[o,d,h]
      let H[o,d,i] := {};
      for {loop in {2..card(Path)}} {
        let AA  := {member(loop-1,Path)};
        let BB  := {member(loop  ,Path)};
        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] + r[o,d];
          }
        }
        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] + r[o,d];
          }
        }
        # 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";
      #display o,d;  
      for {j in {1..card(Path)-1}}{
        let TL:=TL union {(min(member(j,Path),member(j+1,Path)),max(member(j,Path),member(j+1,Path)))};
        if (min(member(j,Path),member(j+1,Path)),max(member(j,Path),member(j+1,Path))) in LG then{
          let LEASEDL:= LEASEDL union {member(j,Path),member(j+1,Path)};
          let LEASEDLINK:= LEASEDLINK union {(min(member(j,Path),member(j+1,Path)),max(member(j,Path),member(j+1,Path)))};
        }
      }
      if card(LEASEDL) > 0 then {
        for {(jjj,kkk) in LEASEDLINK}{
          let LeasedDistance:= LeasedDistance + cc[jjj,kkk];
          let LEASEDN:= LEASEDN union {Hut[jjj],Hut[kkk]};
        }
        let CLeased:= ceil(LeasedDistance/LinkBudget[card(LB24)]) - 1; # no of huts
        let FLeased:= card(LEASEDL);
        #display CLeased,LEASEDN,LEASEDL;
      }
      # Assuming a hut in every 100km on a leased line
      # display LeasedDistance;
      # display TL;
      # display LEASEDL;

      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;
      }
      let TL:={};
      for {kk in {1..MaxK}} {
        for {(jjj,kkk) in Link[o,d,i,kk]} {
          let TL:= TL union {(jjj,kkk)};
          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];
      #UN
      for {h in TH} {
        let svR[h]  := R[h];
        let svA[h]  := A[h];
        let svTE[h] := TE[h];
        let R[h]  := 0;
        let A[h]  := 0;
        let TE[h] := 0;
      }
      for {(ii,jj) in TL}{
        let svMax[ii,jj] := Max[ii,jj];
        let Max [ii,jj] := r[o,d];
      }
      let chkL:=0;
      let chkL2:=0;
      let LsdTE:=0;
      #UN
      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]) or h in LEASEDN then {
          if chkL = 1 and chkL2 = 1 then {
            let TE[nodeL] := TE[nodeL] - r[o,d];
            #printf"%4d %3d %4d %3d %1d %5d %5d %1d", loop6-1,nodeL,dist[o,d,i,nodeL],yy[nodeL],YY[nodeL],ndist[loop6-1],rdist[key],K[o,d,i,nodeL];
            #printf" %8.1f %4.1f %4d %4d %4d %11d *Updated\n", dpmdLastR,dpmd[o,d,i,key],R[nodeL],A[nodeL],TE[nodeL],LinkB[loop6-1];
            let GLeased:= GLeased + 1;
            let LCdistance:= LCdistance + dist[o,d,i,h];
            #display LCdistance;
            #let LsdTE:= LsdTE + 1; 
          }
          if h in LEASEDN and h == first(H[o,d,i]) then {
            let TE[h] := TE[h] - r[o,d];
            let HLeased:= HLeased + 1;
            #let LsdTE:= LsdTE + 1;
          }
          else if h in LEASEDN and h == last(H[o,d,i]) then {
            let TE[h] := TE[h] - r[o,d];
            let HLeased:= HLeased + 1;
            #let LsdTE:= LsdTE + 1;
          }
          if chkL = 1 and chkL2 <> 1 then {
            let chkL2:= 1;
            let LCdistance:= dist[o,d,i,h];
          }
          if h in LEASEDN then {
            let chkL:=1; let nodeL:= h;
          }
          let TE[h] := TE[h] + r[o,d];
          #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;
        }
        # if  LCdistance > 0 then display LCdistance;
        let LeasedRi:= floor(LCdistance/ (card(LB24)*LinkBudget[card(LB24)])) * r[o,d];
        #need an R for every 2400Km for Leased Lines
        if  LeasedRi > 0 then display LeasedRi;
        let LCdistance:= 0;
        if chkL = 1 then let KLeased:= KLeased + 1;
        let chkL:=0;
        let chkL2:=0;
        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;
        } #27
        #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";
        } #24
        #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] + r[o,d];
          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] + r[o,d];
        }
        #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  }
      #UN
      let AllOpticalCost   := 0;
      let TotalA01To20 := 0;
      let TotalM01To20 := 0;
      let TotalA21To40 := 0;
      let TotalM21To40 := 0;
      let TotalA41To80 := 0;
      let TotalM41To80 := 0;
      let TotalFiber   := 0;
      let TotalTE      := 0;
      let TotalR       := 0;
      let TotalA       := 0;
      
      for {(j,k) in TL} {
      # 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;
        }
      # display TotalA01To20, TotalA21To40, TotalA41To80;
      } # 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
      for {(ii,jj) in TL}
       let TotalFiber:= TotalFiber + length[ii,jj];
      # Leased Equipment
      if  LCdistance > 0 then {
        #display LCdistance;
        let LeasedR:= (LeasedR + LeasedRi + floor(LCdistance/2400)) * r[o,d];
        #display LeasedR;
      }
      if card(LEASEDL) > 0 then {
        if BLeased <= 20 then {
          let LeasedA20:= LeasedA20 + CLeased + (2*FLeased-HLeased-KLeased) + floor(LeasedR/r[o,d]);
          let LeasedM20:= LeasedM20 + (2*FLeased-HLeased-KLeased) + 2 * floor(LeasedR/r[o,d]);    
        }
        else if BLeased <= 40 then{
          let LeasedA40:= LeasedA40 + CLeased + (2*FLeased-HLeased-KLeased) + floor(LeasedR/r[o,d]);
          let LeasedM40 := LeasedM40 + (2*FLeased-HLeased-KLeased) + 2 * floor(LeasedR/r[o,d]);
        }
        else if BLeased < 80 then{
          let LeasedA80 := LeasedA80 + CLeased + (2*FLeased-HLeased-KLeased) + floor(LeasedR/r[o,d]);
          let LeasedM80 := LeasedM80 + (2*FLeased-HLeased-KLeased) + 2 * floor(LeasedR/r[o,d]);
        }
        let LeasedA80 := LeasedA80 + ALeased*(CLeased + (2*FLeased-HLeased-KLeased-KLeased) + floor(LeasedR/r[o,d]));
        let LeasedM80 := LeasedM80 + ALeased*((2*FLeased-HLeased-KLeased) + 2 * floor(LeasedR/r[o,d]));
        
        let TotalOXC := FLeased;
        let TotalR := TotalR + LeasedR;
        let LeasedTE := LeasedTE + (FLeased - GLeased) * r[o,d];
        #display LeasedA20, LeasedA40, LeasedA80;
        #display LeasedM20, LeasedM40, LeasedM80;
        #display LeasedTE;
        #display LeasedR;
        #display TotalOXC;
        let LeasedLineCost:= LeasedLineCost + LeasedA20*costA21To40 + LeasedA40*costA21To40 + LeasedA80*costA41To80+
        LeasedM20*costM01To20 + LeasedM40*costM21To40 + LeasedM80*costM41To80 + LeasedTE * costTE + LeasedR * costR;
        let TotalA := TotalA + LeasedA20 + LeasedA40 + LeasedA80;
        let TotalM01To20:= TotalM01To20 + LeasedM20;
        let TotalM21To40:= TotalM21To40 + LeasedM40;
        let TotalM41To80:= TotalM41To80 + LeasedM80;
        let TotalTE:= TotalTE + LeasedTE;
      }  
      if UN[o,d] == 0 then let UN[o,d]:= (QFib  * TotalFiber + QReg * TotalR + QTe * TotalTE + QMUX * (TotalM01To20 + TotalM21To40 + TotalM41To80) + QOXC * (card(Path)-2) + QAMP * (TotalA01To20 + TotalA21To40 + TotalA41To80));
      else let UN[o,d] := UN[o,d] * (QFib  * TotalFiber + QReg * TotalR + QTe * TotalTE + QMUX * (TotalM01To20 + TotalM21To40 + TotalM41To80) + QOXC * (card(Path)-2) + QAMP * (TotalA01To20 + TotalA21To40 + TotalA41To80));
      #display TotalFiber, TotalR, TotalTE, TotalM01To20 + TotalM21To40 + TotalM41To80, TotalA, (card(Path)-2);
      #display UN;
      for {h in TH} {
        let R[h]  := R[h] + svR[h];
        let A[h]  := A[h] + svA[h];
        let TE[h] := TE[h] + svTE[h];
      }
      for {(ii,jj) in TL}{
        let Max[ii,jj]:= svMax[ii,jj];
      }
    }
    let TotalLeasedCost:= TotalLeasedCost + 1.25 * LeasedLineCost; 
  }
  let {(o,d) in D} UN[o,d]:= UN[o,d] + 2 * QOXC;
  let TL:=L;
  let AllOpticalCost   := 0;
  let TotalA01To20 := 0;
  let TotalM01To20 := 0;
  let TotalA21To40 := 0;
  let TotalM21To40 := 0;
  let TotalA41To80 := 0;
  let TotalM41To80 := 0;
  let TotalFiber   := 0;
  let TotalTE      := 0;
  let TotalR       := 0;
  let TotalA       := 0;
  
  for {(j,k) in TL} {
  # 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;
    }
  # display TotalA01To20, TotalA21To40, TotalA41To80;
  } # 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
  for {(ii,jj) in TL}
   let TotalFiber:= TotalFiber + length[ii,jj];
  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;
  display TotalLeasedCost;
  let AllOpticalCost:= (AllOpticalCost + TotalLeasedCost);
  display AllOpticalCost;
  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;
  let TimeInAllOptical[runNo] := _ampl_user_time - TotalTime - TimeInOpaque[runNo];
  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[runNo];
  printf"Time for Routing(CPLEX)    = %4.2lf\n", TimeForRouting[runNo];
  printf"Time in Opaque       = %4.2lf\n", TimeInOpaque[runNo];
  printf"Time In AllOptical     = %4.2lf\n", TimeInAllOptical[runNo];
  printf"Total Time = %4.2lf\n\n", _ampl_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;
  #display UN;
  #display BestA;
  for {(o,d) in D} if UNmax < UN[o,d] then let UNmax:= UN[o,d];
  let UNmax := UNmax * 525960;
  printf" \nMax Unavailability for All-Optical (min): %.8f\n\n ", UNmax;
  let UNt:= {};
  for {(o,d) in D}
    if UN[o,d] > trsh then let UNt:= UNt union {(o,d)};
  #display UNt;
  display _ampl_user_time;
  display nonewp;
  if card(UNt) = 0 then quit;
  let w:= 0;
  reset data P, PD, rhs, ArcsInPath, ArcsWP, MPath, PathsForARC, pn;
  for {(o,d) in UNt}{
    let pm[o,d]:= pm[o,d] + 1;
    let nonewp:= nonewp + 1;
  }
  #display pn, pm;
  let {(o,d) in D} UN[o,d]:=0;
  reset data K,dist,H,dpmd,Link,Max;
  let TotalTime :=_ampl_user_time;
}