# ------------------------------------------------------
# Candidate Path Generator For Links in the Networks
# For Link-Based Shared Protection
# Used to create DataEUSP, DataUSSP, DataNASP data files
# from DataEU, DataUS, DataNA data files
# Output is CandidatePathsforLBS.txt file
# ------------------------------------------------------
# option solver cplex;
option presolve 0;
param NoPaths;       # max no of paths per (o,d) pair
set N ordered;       # set of cities
set L within {N,N};  # set of links
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 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
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

#load Data file
# EU, US or NA
data DataEU.txt;

param length {L} default 0;   
                     # length[j,k] denotes the length of link (j,k)
for {(j,k) in L} {
  for {s in {Begin[j,k]..End[j,k]}}
   let length[j,k] := length[j,k] + Spaces[s];
}

# Determine candidate paths

set P default {};    # set of paths {1..no of paths};
set E within {N,N} default {};  
                     # set of arcs
set V within {N,N} default {};
set Q {N,N} within P default {};    
                     # working paths for an (o,d)
set PQ {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 AQ {E} within P default {};
param w default 0;

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

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 V} 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 V} cx[i,j] - sum{j in N: (j,i) in V} cx[j,i] = rhs[i];

subject to diff_path {(ii,jj) in L, p in Q[ii,jj]}:
  sum {(i,j) in PQ[p]} cx[i,j] <= card(PQ[p]) -1;

#---Shortest Path Problem-----
  problem WorkingPath {(ii,jj) in L}:
#  objective
    cost,
#  variables
    cx,
#  constraints
    flow_balance, {p in Q[ii,jj]} diff_path[ii,jj,p];
for {(ii,jj) in L} {
	let V := E diff {(ii,jj)};
	#let P:= P union {w};
  #let PQ[w]:= {(ii,jj)};
  #let W[ii,jj] := W[ii,jj] union {w};
  #let AQ[ii,jj] := AQ[ii,jj] union {w};
  for {spath in {1..NoPaths}} {
    printf"loop     #%d  ",spath;
    printf"Link (%d,%d)\n",ii,jj;
    display V;
    for {n in N} let rhs[n] := 0;
    let rhs[ii] := +1;
    let rhs[jj] := -1;
    problem WorkingPath[ii,jj];
    solve   WorkingPath[ii,jj];
    # 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 Q[ii,jj] := Q[ii,jj] union {w};
    let PQ[w] := {};
    for {(i,j) in V} {
      if cx[i,j] > 0 then{
        let PQ[w] := PQ[w] union {(i,j)};
        let AQ[i,j] := AQ[i,j] union {w};
      }
    }
    let Path[w] := {ii};
    for {(j,k) in PQ[w]:j == jj}{
      let Path[w] := Path[w] union {k};
      break;
    }
    repeat while card(Path[w]) < card(PQ[w]) + 1 {    
      for {k in N: (last(Path[w]),k) in PQ[w]}{
      let Path[w] := Path[w] union {k};
    }
    } # repeat   
  } 
} # for {(o,d) in D}
printf "\nparam NoP:= %d;\n",w>ebi.txt;
display Q, PQ, AQ>>ebi.txt;
quit;

let P:= P union {w+1..2*w};
for {(i,j) in L}{
  let Q[j,i] := {};
  let AQ[j,i] := {};
  for {p in Q[i,j]}
  	let Q[j,i] := Q[j,i] union {p+w};
  for {p in AQ[i,j]}
  	let AQ[j,i] := AQ[j,i] union {p+w};
}

for {p in {1..w}}{
  let PQ[p+w] := PQ[p];
}
printf "\nparam NoP:= %d;\n",2*w>CandidatePathsforLBS.txt;
display Q, PQ, AQ>>CandidatePathsforLBS.txt;