Optimization: principles and algorithms, by Michel Bierlaire
simpleCycle.m
Go to the documentation of this file.
1 %> \file
2 %> Algorithm 21.1: Generation of a simple cycle flow from a circulation. Implementation of algorithm 21.1 of \cite Bier15-book
3 %>
4 %> @note Tested with \ref runSimpleCycle.m
5 %> @note Calls \ref prepareNetwork
6 %> @note Called by \ref circulationDecomposition
7 %>
8 %> @author Michel Bierlaire
9 %> @date Fri Mar 27 17:09:10 2015
10 %> @ingroup Algorithms
11 %> @ingroup chap21
12 
13 %> Extract a simple cycle flow from a circulation
14 %> @param adj the adjacency matrix of the network. It is a \f$m \times m\f$ matrix, such that the element at row i and column j corresponds to the id of the arc (i,j). The numbering should be from 1 to n.
15 %> @param circ the circulation flow. It is a vector of \f$\mathbb{R}^n\f$.
16 %> @param printlevel If different from 0, the cycle is printed. (Default: 0)
17 %> @return cycleflow a simple cycle flow. It is a vector of \f$\mathbb{R}^n\f$.
18 function cycleflow = simpleCycle(adj,circ,printlevel=0)
19  cycleflow = zeros(size(circ)) ;
20  nnodes = rows(adj) ;
21  if (columns(adj) != nnodes)
22  error("Adjacency matrix must be square") ;
23  endif
24  arcs = prepareNetwork(adj) ;
25  nonzeros = find(circ) ;
26  if(length(nonzeros)==0)
27  printf("The circ flow vector is zero")
28  return
29  endif
30  # Find an arc carrying positive flow
31  theArc = nonzeros(1) ;
32  theLoopArcs = [theArc] ;
33 
34  layerNumber = 1 ;
35  if (circ(theArc) > 0)
36  layer{layerNumber} = arcs(theArc,2) ;
37  otherNode = arcs(theArc,1) ;
38  theLoopOrientation = [1] ;
39  else
40  layer{layerNumber} = arcs(theArc,1) ;
41  otherNode = arcs(theArc,2) ;
42  theLoopOrientation = [-1] ;
43  endif
44  inserted = 1 ;
45  treatedNodes = zeros(nnodes,1) ;
46  treatedNodes(layer{layerNumber}) = 1 ;
47  while (inserted != 0)
48  inserted = 0 ;
49  currentlayer = zeros(0) ;
50  for i = layer{layerNumber}
51  for arcij = adj(i,:)
52  if (arcij != 0)
53  if (circ(arcij) > 0) && (treatedNodes(arcs(arcij,2)) == 0)
54 % printf("Forward arc: %d -> %d\n",arcs(arcij,1),arcs(arcij,2));
55  currentlayer(end + 1) = arcs(arcij,2) ;
56  inserted = 1 ;
57  treatedNodes(arcs(arcij,2)) = 1 ;
58  endif
59  endif
60  endfor
61  for arcij = adj(:,i)'
62  if (arcij != 0)
63  if (circ(arcij) < 0) && (treatedNodes(arcs(arcij,1)) == 0)
64 % printf("Backward arc: %d -> %d\n",arcs(arcij,1),arcs(arcij,2));
65  currentlayer(end + 1) = arcs(arcij,1) ;
66  inserted = 1 ;
67  treatedNodes(arcs(arcij,1)) = 1 ;
68  endif
69  endif
70  endfor
71  endfor
72  layerNumber += 1 ;
73  layer(layerNumber) = currentlayer;
74  endwhile
75 
76 % Here, layer contains the nodes organized into layers
77 % Search for the node to close the loop
78 
79  for key = 1:length(layer)
80  theLayer = layer{key} ;
81 % printf("Look for %d in \n",otherNode)
82  idx = find(theLayer == otherNode) ;
83  if (length(idx) > 0)
84  theLayerOfOtherNode = key ;
85  theIndexOfOtherNode = idx ;
86  endif
87  endfor
88 
89  theLoop = [otherNode] ;
90  flowMin = realmax ;
91  for k = theLayerOfOtherNode-1:-1:1
92 % Find an arc carrying positive flow. It must exist by construction
93  for ell = layer{k}
94  theForwardArc = adj(ell,theLoop(end)) ;
95  if (theForwardArc != 0)
96  forwardflow = circ(theForwardArc) ;
97  if (forwardflow > 0)
98  theLoop(end+1) = ell ;
99  theLoopArcs(end+1) = theForwardArc ;
100  theLoopOrientation(end+1) = 1 ;
101  if (forwardflow < flowMin)
102  flowMin = forwardflow ;
103  endif
104  break ;
105  endif
106  endif
107  theBackwardArc = adj(theLoop(end),ell) ;
108  if (theBackwardArc != 0)
109  backwardflow = circ(theBackwardArc) ;
110  if (backwardflow < 0)
111  theLoop(end+1) = ell ;
112  theLoopArcs(end+1) = theBackwardArc ;
113  theLoopOrientation(end+1) = -1 ;
114  if (-backwardflow < flowMin)
115  flowMin = -backwardflow;
116  endif
117  break ;
118  endif
119  endif
120  endfor
121  endfor
122  for i = length(theLoopArcs):-1:1
123  if (printlevel != 0)
124  printf("%d",theLoop(i))
125  if (theLoopOrientation(i) == 1)
126  printf("->") ;
127  else
128  printf("<-") ;
129  endif
130  endif
131  cycleflow(theLoopArcs(i)) = flowMin * theLoopOrientation(i) ;
132  endfor
133  if (printlevel != 0)
134  printf("%d [%f]\n",theLoop(end),flowMin)
135  endif
136 endfunction
137 
function circulationDecomposition(in adj, in circ)
Decompose a circulation into simple cycle flows.
function simpleCycle(in adj, in circ, in printlevel)
Extract a simple cycle flow from a circulation.
function prepareNetwork(in adj)
Identifies the upstream and downstream nodes of each arc.
Copyright 2015-2016 Michel Bierlaire