#----------------------------------------------
#  AMPL Model File for the 
#  Path-Based Shared Protection Opaque Design
#  *** MULTIPLE-PATH MODEL ***
#----------------------------------------------
# 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;

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
param NoPaths;       # max number of cand. paths per demand
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;

#-------------------------------------------------------------
# Load Data File - DataEU.txt, DataNA.txt, or DataUS.txt
data DataEU.txt;
#-------------------------------------------------------------

display NoDemands;   # number of demand pairs
display RandomSeed;  # random seed used to generate demand pairs
param length {L} default 0;   
                     # length[j,k] denotes the length of link (j,k)
param TimeInOpaque 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)];
}
#let D:= D union {(3,4)};
#let r[3,4]:= 37;
#let TotalDemand := TotalDemand + r[3,4];
#let D:= D union {(4,1)};
#let r[4,1]:= 63;
#let TotalDemand := TotalDemand + r[4,1];
#let NoDemands:= NoDemands + 2;
display D,r;
#display TotalDemand;
param TimeToGenerate;
let TimeToGenerate := _ampl_user_time;
#--------------------------------------------------------------------------------------------------------
# Determine HL and distL {HL}
# and solve the hut selection model for each link budget and for each link

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};  # the best link budget for link l
param YY {TH};
param yy {TH};
param preH;
param curL;
param tNewPaths default 0;
param AMPLoc{LB24,L} default 1.0e20;
param MUXLoc{LB24,L} default 1.0e20;
param RLoc{LB24,L} default 1.0e20;
set ALB{L} within LB24 default {};
param flag default 0;
set RLB default {};

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;
  let flag := 0;
  let RLB:= {};
  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];
    }
    let AMPLoc[b,j,k] := TotalA + 2 + 2 * TotalR;
    let MUXLoc[b,j,k] := 2 + 2 * TotalR;
    let RLoc[b,j,k]   := TotalR;
    let flag:= 0;
    let RLB := {};
    if b = 1 then let ALB[j,k] := ALB[j,k] union {b};
    else{ 
		   for {lb in ALB[j,k]}{ 
		     if AMPLoc[b,j,k] >= AMPLoc[lb,j,k] and MUXLoc[b,j,k] >= MUXLoc[lb,j,k] and RLoc[b,j,k] >= RLoc[lb,j,k] then{
		       let flag:= 1;
		       break;    
		     }
		     else if AMPLoc[b,j,k] <= AMPLoc[lb,j,k] and MUXLoc[b,j,k] <= MUXLoc[lb,j,k] and RLoc[b,j,k] <= RLoc[lb,j,k] then
		       let RLB:= RLB union {lb};
		   }
		   if card(RLB) > 0 then let ALB[j,k]:= {ALB[j,k] union {b}} diff RLB;
		   else if flag = 0 then let ALB[j,k]:= ALB[j,k] union {b};
    }
  }    # 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}
#for {(i,j) in L}
#	for {b in LB24}{
#	  if AMPLoc[b,i,j] = 1.0e20 then continue;
#	  printf"\nLink : (%s,%s), Budget : %d, AMP : %d, Mux : %d, R : %d\n", i, j, b, AMPLoc[b,i,j] , MUXLoc[b,i,j] , RLoc[b,i,j];
#  }
#display AMPLoc, MUXLoc, RLoc;
#display ALB;
param TimeForBestBudget;
let TimeForBestBudget := _ampl_user_time + _total_solve_user_time - TimeToGenerate;
#--------------------------------------------------------------------------------------------------------
# Candidate paths model
set P default {};   # set of paths {1..no of paths};
param pn default 2; # denotes the no of individual paths for path p
set E within {N,N} default {};  
            # set of arcs
set CD {D} within P default {};   
            # working & protection paths for an (o,d)
            
set ArcsInPath {P,1..pn} within {E};
            # new denotes the set of links in 
            # path p,n for (o,d)
            
set ArcsWP {P} within {E};

set MPath {P,1..pn} within {N} ordered;
                     # WPath[p] denotes the nodes in the nth working path

set PathsForARC {E,1..pn} within P default {};
set Path within {N} ordered;
param newpaths;
param DistForPath {P,1..pn} default 0;
set TempSet1 dimen 2;
set TempSet2;
param temp;

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 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 CD[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 CD[o,d]} diff_path[o,d,p], nodelimit_in, nodelimit_out;

#--------------------------------------------------------------------------------------------------------
# Determine candidate paths
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] := pn;
    let rhs[d] := -1*pn;
    #display P, CD, 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};
    # display o,d;
    let CD[o,d] := CD[o,d] union {w};
    let ArcsWP[w] := {};
    let {i in 1..pn} 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..pn} MPath[w,i] := {o};
    for {(j,k) in ArcsWP[w]:j = o}{
      for {i in 1..pn}{
        if card(MPath[w,i]) = 1 then {
            let MPath[w,i] := MPath[w,i] union {k};
            #let PathsForARC[j,k,i] := PathsForARC[j,k,i] union {w};
            let ArcsInPath[w,i] := ArcsInPath[w,i] union {(j,k)};
            let DistForPath[w,i] := DistForPath[w,i] + length[min(j,k),max(j,k)];
            break;
        }
        else continue;
      }
    }
    repeat while (sum {i in 1..pn} card(MPath[w,i])) < card(ArcsWP[w]) + 2 {
      for {j in 1..pn}{
        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,j] := PathsForARC[last(MPath[w,j]),k,j] union {w};
            let DistForPath[w,j] := DistForPath[w,j] + length[min(last(MPath[w,j]),k),max(last(MPath[w,j]),k)];
            let MPath[w,j] := MPath[w,j] union {k};
        }
      }   
    }
    # Only 2 paths per cycle for the Global Path-Based Shared Protection Model
    let flag := 0;
    if DistForPath[w,1] > DistForPath[w,2] then {
      let TempSet1 := ArcsInPath[w,1];
      let ArcsInPath[w,1] := ArcsInPath[w,2];
      let ArcsInPath[w,2] := TempSet1;
      let TempSet2 := MPath[w,1];
      let MPath[w,1] := MPath[w,2];
      let MPath[w,1] := TempSet2;
      let flag := 1;
    }
    if flag = 1 then {
      let temp:= DistForPath[w,1];
      let DistForPath[w,1] := DistForPath[w,2];
      let DistForPath[w,2] := temp;
    }
  } 
  if newpaths = w then quit;
} # for {(o,d) in D}
for {p in P, q in {1..2}}
  for {(i,j) in ArcsInPath[p,q]}
    let PathsForARC[i,j,q] := PathsForARC[i,j,q] union {p};

#display P, CD, rhs, ArcsInPath, ArcsWP, MPath, PathsForARC, pn, DistForPath;
param TimeForPaths;
let TimeForPaths := _ampl_user_time + _total_solve_user_time - TimeForBestBudget;
#--------------------------------------------------------------------------------------------------------
# Pre-processing

param TPDarc {E} default 0;
set Tset dimen 2 default {};
param TPPlink {L} default 0;

for {(i,j) in E}{
  #display i,j;
	let Tset := {};
	for {p in PathsForARC[i,j,1]}
		for {(o,d) in D: p in CD[o,d]} {
			if (o,d) in Tset then continue;
			else {
			  #display o,d,r[o,d];
				let Tset:= Tset union {(o,d)};
				let TPDarc[i,j] := TPDarc[i,j] + r[o,d];
			}	
		}
}
set Tmp;
set Tod dimen 2 default {};
set Pod{L} dimen 2 default {};
param TPl {L} default 0;

for {(i,j) in L}{
  #printf"\n(i,j)  = %d,%d\n",i,j;
  let Tod := {};
  let {(ii,jj) in L} Pod[ii,jj] := {};
  let {(ii,jj) in L} TPl[ii,jj] := 0;
  let Tmp := {PathsForARC[i,j,1] union PathsForARC[j,i,1]};
  #display Tmp;
  for {t in Tmp}
  	for {(o,d) in D: t in CD[o,d]}
  	  let Tod:= Tod union {(o,d)};
  #display Tod; 
  if card(Tod) > 1 then {
    for {t in Tmp, (o,d) in D:t in CD[o,d]}
      for {(ii,jj) in L}
        if ((ii,jj) in ArcsInPath[t,2] or (jj,ii) in ArcsInPath[t,2]) and (o,d) in Tod then {
          if (o,d) not in Pod[ii,jj] then {
            let Pod[ii,jj]:= Pod[ii,jj] union {(o,d)};
            #display Pod[ii,jj];
            #printf"\n(ii,jj)  = %d,%d\n",ii,jj;
            let TPl[ii,jj] := TPl[ii,jj] + r[o,d];
          }
        }   
  }
  for {(ii,jj) in L}
    if TPl[ii,jj] > TPPlink[ii,jj] then let TPPlink[ii,jj] := TPl[ii,jj];
  #display TPPlink;
  for {p in {PathsForARC[i,j,2] union PathsForARC[j,i,2]}}
    for {(o,d) in D: p in CD[o,d]}
      if r[o,d] > TPPlink[i,j] then let TPPlink[i,j] := r[o,d];
}
for {(i,j) in L}
  if TPPlink[i,j] > TotalDemand then{
    printf"\nBug in Preprocessing!";
    quit;
  }
#display TPDarc, TPPlink;
#display TotalDemand;
param TimeForPreprocess;
let TimeForPreprocess := _ampl_user_time + _total_solve_user_time - TimeForPaths;
#--------------------------------------------------------------------------------------------------------
# Routing and Provisioning Model
var LBp {(i,j) in L, ALB[i,j]} binary;
                     # LBp[b,i,j] is 1 if the link budget b is used for link (i,j) and 0 otherwise
var TF {(i,j) in L,  ALB[i,j]} >= 0 integer;    
                     # TF[i,j] denotes the total flow on arc (i,j)
var yp {(i,j) in L, ALB[i,j]} >=0;  
                     # denotes protection traffic on link (i,j)
var A01_20 {(i,j) in L, ALB[i,j]} binary; 
                     # denotes the number of A1-20s needed at a utilized hut on link l
var A21_40 {(i,j) in L, ALB[i,j]} binary; 
                     # denotes the number of A21-40s needed at a utilized hut on link l
var A41_80 {(i,j) in L, ALB[i,j]} >=0 integer; 
                     # denotes the number of A41-80s needed at a utilized hut on link l
                                  
var PF {(i,j) in E,PathsForARC[i,j,2]} >= 0 integer;  # PF[i,j,p] denotes the total protection flow on arc (i,j) for path p

var TErE {N} >=0;    # denotes the number of TEs needed at node n
var CO{P} >=0, <=1;  # PO[o,d,p] is 1 if the path p for demand (o,d) is used and 0 otherwise

                     
subject to BestLinkBudget {(i,j) in L}:
  sum {b in ALB[i,j]} LBp[i,j,b] = 1;

subject to EnableWFlow {(i,j) in L, b in ALB[i,j]}:  TF[i,j,b]  <= (TPDarc[i,j] + TPDarc[j,i]) * LBp[i,j,b];

subject to EnablePFlow {(i,j) in L, b in ALB[i,j]}:  yp[i,j,b]  <= TPPlink[i,j] * LBp[i,j,b];

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

subject to TFinArc {(i,j) in L}:
  sum {b in ALB[i,j]} TF[i,j,b]  = sum {(o,d) in D, p in CD[o,d]:(i,j) in ArcsInPath[p,1]} r[o,d] * CO[p]
    + sum {(o,d) in D, p in CD[o,d]:(j,i) in ArcsInPath[p,1]} r[o,d] * CO[p];

subject to PFinArc {(i,j) in E,p in PathsForARC[i,j,2]}:
  PF[i,j,p] = sum {(o,d) in D:p in CD[o,d]} r[o,d] * CO[p];

subject to Capacity {(i,j) in L,(ii,jj) in L diff {(i,j)}}:
  sum {p in {PathsForARC[ii,jj,1] union PathsForARC[jj,ii,1]}:p in PathsForARC[i,j,2]} PF[i,j,p] + sum {p in {PathsForARC[ii,jj,1] union PathsForARC[jj,ii,1]}:p in PathsForARC[j,i,2]} PF[j,i,p] <= sum {b in ALB[i,j]} yp[i,j,b];

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

subject to TEsonN {k in N}:
  sum {(i,j) in L, b in ALB[i,j]: i=k or j=k} (TF[i,j,b] + yp [i,j,b]) = TErE[k];
  
minimize DesignCost: sum {k in N} (costTE * TErE[k]) 
 + sum {(i,j) in L} ((sum {b in ALB[i,j]} (AMPLoc[b,i,j] * (costA01To20 * A01_20[i,j,b] + costA21To40 * A21_40[i,j,b] + costA41To80 * A41_80[i,j,b]))) 
 + costR * sum {b in ALB[i,j]} ((TF[i,j,b] + yp [i,j,b]) * RLoc [b,min(i,j),max(i,j)]))
 + sum {(i,j) in L} (sum{b in ALB[i,j]} ((costM01To20 * A01_20[i,j,b] + costM21To40 * A21_40[i,j,b] + costM41To80 * A41_80[i,j,b]) * MUXLoc[b,i,j]));
    
#---Routing Problem-----
  problem Routing:
#  objective
    DesignCost,
#  variables
    TF, PF, CO, LBp, yp, A01_20, A21_40, A41_80, TErE,
#  constraints
   BestLinkBudget, EnableWFlow, EnablePFlow, PathOperational, TFinArc, PFinArc, Capacity,  AsonL, SmallAs, TEsonN;

problem Routing;
option cplex_options 'timing=1 mipdisplay=2 mipgap=0.050 display=2';
#option show_stats 1;
solve Routing;
#expand;
display CO, TF, yp, PF;
display LBp;
param TimeForRouting;
param TradA;
param TradR;
param TradTE;
param TradM;
param AVEnHOPS default 0;
param PathsInUse default 0;
let TimeForRouting := _solve_user_time;
let TradA := sum {(i,j) in L, b in ALB[i,j]} (AMPLoc[b,i,j] * (A01_20[i,j,b] + A21_40[i,j,b] + A41_80[i,j,b]));
let TradR := sum {(i,j) in L} (sum {b in ALB[i,j]} ((TF[i,j,b] + yp [i,j,b]) * RLoc [b,min(i,j),max(i,j)]));
let TradTE := sum {k in N} TErE[k];
let TradM := sum {(i,j) in L, b in ALB[i,j]} (MUXLoc[b,i,j] * (A01_20[i,j,b] + A21_40[i,j,b] + A41_80[i,j,b]));
let OpaqueCost := DesignCost;
let TimeInOpaque := _ampl_time + _total_solve_time;
	
# Number of Hops for Opaque 
for {(o,d) in D} { 
  for {p in CD[o,d] : CO[p] = 1} { 
	    let PathsInUse := PathsInUse + 2; 
	    let AVEnHOPS := AVEnHOPS + card(MPath[p,1]) + card(MPath[p,2]) - 2; 
	} 
}
let AVEnHOPS := AVEnHOPS / PathsInUse;
printf"Average Number of Hops = %4.2lf\n\n", AVEnHOPS;
#--------------------------------------------------------------------------------------------------------
display TotalDemand;
printf"Opaque Cost       = %4d\n", OpaqueCost;
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;
printf"Opaque AMPL Time         =   %4.2lf\n", _ampl_time;
printf"Time in Opaque           =   %4.2lf\n", TimeInOpaque;
printf"Time to Generate         =   %4.2lf\n", TimeToGenerate;
printf"Time for Candidate Paths =   %4.2lf\n", TimeForPaths;
printf"Time for Optical Reach   =   %4.2lf\n", TimeForBestBudget;
printf"Time for Preprocessing   =   %4.2lf\n", TimeForPreprocess;
printf"Time for Routing         =   %4.2lf\n", TimeForRouting;