Optimization: principles and algorithms, by Michel Bierlaire
unsaturatedPath.m
Go to the documentation of this file.
1 %> \file
2 %> Algorithm 24.1: unsaturated path. Implementation of algorithm 24.1 of \cite Bier15-book
3 %>
4 %> @note Tested with \ref runUnsaturatedPath.m
5 %> @note Called by \ref fordFulkerson
6 %> @note Calls \ref prepareNetwork
7 %>
8 %> @author Michel Bierlaire
9 %> @date Thu Apr 9 18:08:09 2015
10 %> @ingroup Algorithms
11 %> @ingroup chap24
12 
13 
14 
15 %> Find an unsaturated path from o to d
16 %>
17 %> @note Calls \ref buildPathBackward
18 %> @note Calls \ref printLayers
19 %> @param adj adjacency matrix of the network. Each nonzero entry corresponds to an arc. The value of the entry is the ID of the arc.
20 %> @param orig the origin node
21 %> @param dest the destination node
22 %> @param lb the lower bound on the flow
23 %> @param ub the upper bound on the flow
24 %> @param flow the current flow vector
25 %> @param printlevel If different from 0, the layers are printed. (Default: 0)
26 %> @return path.arcs: sequence of arcs on the path
27 %> @return path.dir: +1 if forward, -1 if backward
28 %> @return cut: array of nodes defining the cut
29 function [path, cut, pathFound] = unsaturatedPath(adj,orig,dest,lb,ub,flow,printlevel=0)
30  addpath("../chap21") ;
31  pathFound = 0 ;
32  cut = zeros(0) ;
33  arcs = prepareNetwork(adj) ;
34  nnodes = rows(adj) ;
35  treatedNodes = zeros(nnodes,1) ;
36 
37  currentlayer.arcs = -1 ;
38  currentlayer.nodes = orig ;
39  currentlayer.dir = 0 ;
40 
41  layer{1} = currentlayer ;
42  treatedNodes(orig) = 1 ;
43  layerNumber = 1 ;
44  inserted = 1 ;
45  while (inserted != 0)
46  inserted = 0 ;
47 # printf("Layer %d\n=========\n",layerNumber);
48  currentlayer.arcs = zeros(0) ;
49  currentlayer.nodes = zeros(0) ;
50  currentlayer.dir = zeros(0) ;
51  for i = layer{layerNumber}.nodes
52 # printf("Treat node %d\n",i) ;
53  for arcij = adj(:,i)'
54  if (arcij != 0)
55  if (flow(arcij) > lb(arcij) && (treatedNodes(arcs(arcij,1)) == 0))
56 # printf("Backward arc: %d -> %d\n",arcs(arcij,1),arcs(arcij,2));
57  currentlayer.arcs(end + 1) = arcij ;
58  currentlayer.nodes(end + 1) = arcs(arcij,1) ;
59  currentlayer.dir(end + 1) = -1 ;
60  if (currentlayer.nodes(end) == dest)
61  layer(end+1) = currentlayer ;
62  path = buildPathBackward(layer,dest,arcs) ;
63  pathFound = 1 ;
64  if (printlevel != 0)
65  printLayers(layer,arcs) ;
66  endif
67  return ;
68  endif
69  inserted = 1 ;
70  treatedNodes(arcs(arcij,1)) = 1 ;
71  endif
72  endif
73  endfor
74  for arcij = adj(i,:)
75  if (arcij != 0)
76  if (flow(arcij) < ub(arcij) && (treatedNodes(arcs(arcij,2)) == 0))
77 # printf("Forward arc: %d -> %d\n",arcs(arcij,1),arcs(arcij,2));
78  currentlayer.arcs(end + 1) = arcij ;
79  currentlayer.nodes(end + 1) = arcs(arcij,2) ;
80  currentlayer.dir(end + 1) = 1 ;
81  if (currentlayer.nodes(end) == dest)
82  layer(end+1) = currentlayer ;
83  path = buildPathBackward(layer,dest,arcs) ;
84  pathFound = 1 ;
85  if (printlevel != 0)
86  printLayers(layer,arcs) ;
87  endif
88  return ;
89  endif
90  inserted = 1 ;
91  treatedNodes(arcs(arcij,2)) = 1 ;
92  endif
93  endif
94  endfor
95  endfor
96  layerNumber += 1 ;
97  layer(layerNumber) = currentlayer ;
98  endwhile
99 % Saturated cut identified
100  path.arcs = zeros(0) ;
101  path.dir = zeros(0) ;
102  if (printlevel != 0)
103  printLayers(layer,arcs) ;
104  endif
105  cut = find(treatedNodes) ;
106 endfunction
107 
108 %> Build the path by traversing the layers backward.
109 %>
110 %> @note Called by \ref unsaturatedPath
111 %> @param layer list of layers. Each layer contains an array of nodes, arcs and directions.
112 %> @param dest destination node to reach
113 %> @param arcs A matrix containing the upstream and downstream nodes of each arc. Each row corresponds to an arc. Column 1 contains the upstream node and column 2 the downstream. Output of \ref prepareNetwork.
114 %> @return path.arcs: sequence of arcs on the path
115 %> @return path.dir: +1 if forward, -1 if backward
116 function path = buildPathBackward(layer,dest,arcs)
117  path.arcs = zeros(0) ;
118  path.dir = zeros(0) ;
119  currentNode = dest ;
120  for i = length(layer):-1:2
121  theLayer = layer{i} ;
122  idx = find(theLayer.nodes == currentNode) ;
123  path.dir = [theLayer.dir(idx) path.dir] ;
124  path.arcs = [theLayer.arcs(idx) path.arcs] ;
125  if (theLayer.dir(idx) == 1)
126  currentNode = arcs(theLayer.arcs(idx),1) ;
127  else
128  currentNode = arcs(theLayer.arcs(idx),2) ;
129  endif
130  endfor
131 endfunction
132 
133 %> Print the layers
134 %>
135 %> @note Called by \ref unsaturatedPath
136 %> @param layer list of layers. Each layer contains an array of nodes, arcs and directions.
137 %> @param arcs A matrix containing the upstream and downstream nodes of each arc. Each row corresponds to an arc. Column 1 contains the upstream node and column 2 the downstream. Output of \ref prepareNetwork.
138 function printLayers(layer,arcs)
139  for i = 1:length(layer)
140  printf("S_%d = {",i) ;
141  theLayer = layer{i} ;
142  for j = 1:length(theLayer.dir)
143  if (theLayer.arcs(j) == -1)
144  printf("%d ",theLayer.nodes(j)) ;
145  else
146  if (theLayer.dir(j) == 1)
147  printf("%d[%d ->], ",arcs(theLayer.arcs(j),2),arcs(theLayer.arcs(j),1)) ;
148  else
149  printf("%d[<- %d], ",arcs(theLayer.arcs(j),1),arcs(theLayer.arcs(j),2)) ;
150  endif
151  endif ;
152  endfor
153  printf("}\n") ;
154  endfor
155 endfunction
function buildPathBackward(in layer, in dest, in arcs)
Build the path by traversing the layers backward.
function unsaturatedPath(in adj, in orig, in dest, in lb, in ub, in flow, in printlevel)
Find an unsaturated path from o to d.
function fordFulkerson(in adj, in orig, in dest, in lb, in ub)
Find the maximum flow between two nodes using the Ford Fulkerson algorithm.
function prepareNetwork(in adj)
Identifies the upstream and downstream nodes of each arc.
function printLayers(in layer, in arcs)
Print the layers.
Copyright 2015-2016 Michel Bierlaire