#----------------------------------------------
#  AMPL Model File for the 
#  Link-Based Shared Protection Opaque Design
#  *** GLOBAL MODEL ***
#----------------------------------------------
#  option solver cplex;
#  option presolve 0;
param loop default 3;
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
param NoP;
set E within {N,N} default {};  
                     # set of arcs
set Q {L} within {1..NoP} default {} ordered;    
                     # protection paths for an (o,d)
set PQ {1..NoP} within {N,N};
                     # denotes the set of links in protection path p for (o,d) 
set AQ {N,N} within {1..NoP} default {};
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

#-------------------------------------------------------------
# Load Data File - DataEUSP.txt, DataNASP.txt, or DataUSSP.txt
data DataNASP.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;
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 {};
param objP default 1.0e20;
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;

#--------------------------------------------------------------------------------------------------------
# Determine candidate paths for routing

set P default {};    # set of paths {1..no of paths};

set W {D} within P default {} ordered;    
                     # 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}

#--------------------------------------------------------------------------------------------------------
# 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]}
		for {(o,d) in D: p in W[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];
			}	
		}
}
for {(i,j) in L}
	for {q in Q[i,j]}
		for {(ii,jj) in PQ[q]}
			if TPPlink[min(ii,jj),max(ii,jj)] < TPDarc[i,j] + TPDarc[j,i] then
				let TPPlink[min(ii,jj),max(ii,jj)] := TPDarc[i,j] + TPDarc[j,i];

display TPDarc, TPPlink;
#--------------------------------------------------------------------------------------------------------
# Routing and Provisioning Model
param TimeForPaths;
let TimeForPaths := _ampl_user_time + _total_solve_user_time - TimeForBestBudget;

var TF {(i,j) in L, ALB[i,j]} >= 0;    
                     # TF[i,j] denotes the total flow on arc (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 yp {(i,j) in L, ALB[i,j]} >=0;  # denotes protection traffic on link (i,j)
                     
var TErE {N} >=0;    # denotes the number of TEs needed at node n
var PO{(o,d) in D,W[o,d]} binary;  # PO[o,d,p] is 1 if the path p for demand (o,d) is used and 0 otherwise
var PP {1..NoP} binary;
                     # PP[p] is 1 if the protection path p is used and 0 otherwise

#var xp{(i,j) in L, Q[i,j]} >= 0;
                     # xp[i,j,p] denotes the protection traffic on path p for link (i,j)
                     
var xp{1..NoP} >= 0;
                     # xp[p] denotes the protection traffic on path p 
                     
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

subject to EnableFlowN {(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 BestLinkBudget {(i,j) in L}:
  sum {b in ALB[i,j]} LBp[i,j,b] = 1;
  
subject to OnePathOperational {(o,d) in D}:
  sum {p in W[o,d]} PO[o,d,p] = 1;
  
subject to OneProtectionOperational {(i,j) in L}:
  sum {p in Q[i,j]} PP[p] = 1;

subject to TFinProtectionPath {p in {1..NoP}}:
   xp [p] <= PP[p] * sum {(i,j) in L:p in Q[i,j]} (TPDarc[i,j] + TPDarc[j,i]);
   
subject to TFinArc {(i,j) in L}:
  sum {b in ALB[i,j]} TF[i,j,b] = sum {p in PathsForARC[i,j], (o,d) in D : p in W[o,d]} r[o,d] * PO[o,d,p] + 
  sum {p in PathsForARC[j,i], (o,d) in D : p in W[o,d]} r[o,d] * PO[o,d,p];

subject to ProtectionCapacity {(i,j) in L}:
    sum {p in Q[i,j]} xp[p] = sum {b in ALB[i,j]} (TF[i,j,b]);

subject to LCapacityNDirection {(i,j) in L, (ii,jj) in L}:
    sum {p in AQ[i,j]:p in Q[ii,jj]} xp[p] <= sum {b in ALB[i,j]} yp [i,j,b];
    
subject to LCapacityRDirection {(i,j) in L, (ii,jj) in L}:
    sum {p in AQ[j,i]:p in Q[ii,jj]} xp[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
    PO, PP, TF, xp, yp, A01_20, A21_40, A41_80, TErE, LBp,
#  constraints
    EnableFlowN, EnablePFlow, BestLinkBudget, OnePathOperational, 
    OneProtectionOperational, TFinProtectionPath, TFinArc, ProtectionCapacity, LCapacityNDirection, LCapacityRDirection, AsonL, SmallAs, TEsonN;

param TimeForRouting;
param TradA;
param TradR;
param TradTE;
param TradM;
param AVEnHOPS default 0;
param PathsInUse default 0;
param looprun default 1;

problem Routing;
option cplex_options 'timing=1 mipdisplay=2 mipgap=0.050 display=2';
option show_stats 1;
solve Routing;
#expand;
#display P, Path, PO, TF, yp, PP, LBp;
#display xp;
#display A01_20, A21_40, A41_80, TErE;
display PO;
display PP;

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;
printf"Opaque Cost       = %4d\n", OpaqueCost;
printf"Opaque AMPL Time  = %4.2lf\n", _ampl_time;
printf"Opaque CPLEX Time = %4.2lf\n\n", _solve_time;
#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;
let TimeInOpaque := _ampl_time + _total_solve_time;
	

# Number of Hops for Opaque 
for {(o,d) in D} { 
  for {p in W[o,d] : PO[o,d,p] = 1} { 
	    let PathsInUse := PathsInUse + 1; 
	    let AVEnHOPS := AVEnHOPS + card(Path[p]) - 1; 
	} 
}
let AVEnHOPS := AVEnHOPS / PathsInUse;
printf"Average Number of Hops = %4.2lf\n\n", AVEnHOPS;
#--------------------------------------------------------------------------------------------------------
display TotalDemand;
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;
printf"Time in Opaque           = %4.2lf\n", TimeInOpaque;