44 /// from a given source node to a given target node in a digraph. |
44 /// from a given source node to a given target node in a digraph. |
45 /// |
45 /// |
46 /// Note that this problem is a special case of the \ref min_cost_flow |
46 /// Note that this problem is a special case of the \ref min_cost_flow |
47 /// "minimum cost flow problem". This implementation is actually an |
47 /// "minimum cost flow problem". This implementation is actually an |
48 /// efficient specialized version of the \ref CapacityScaling |
48 /// efficient specialized version of the \ref CapacityScaling |
49 /// "Successive Shortest Path" algorithm directly for this problem. |
49 /// "successive shortest path" algorithm directly for this problem. |
50 /// Therefore this class provides query functions for flow values and |
50 /// Therefore this class provides query functions for flow values and |
51 /// node potentials (the dual solution) just like the minimum cost flow |
51 /// node potentials (the dual solution) just like the minimum cost flow |
52 /// algorithms. |
52 /// algorithms. |
53 /// |
53 /// |
54 /// \tparam GR The digraph type the algorithm runs on. |
54 /// \tparam GR The digraph type the algorithm runs on. |
55 /// \tparam LEN The type of the length map. |
55 /// \tparam LEN The type of the length map. |
56 /// The default value is <tt>GR::ArcMap<int></tt>. |
56 /// The default value is <tt>GR::ArcMap<int></tt>. |
57 /// |
57 /// |
58 /// \warning Length values should be \e non-negative. |
58 /// \warning Length values should be \e non-negative. |
59 /// |
59 /// |
60 /// \note For finding node-disjoint paths this algorithm can be used |
60 /// \note For finding \e node-disjoint paths, this algorithm can be used |
61 /// along with the \ref SplitNodes adaptor. |
61 /// along with the \ref SplitNodes adaptor. |
62 #ifdef DOXYGEN |
62 #ifdef DOXYGEN |
63 template <typename GR, typename LEN> |
63 template <typename GR, typename LEN> |
64 #else |
64 #else |
65 template < typename GR, |
65 template < typename GR, |
107 typedef typename Digraph::template NodeMap<int> HeapCrossRef; |
107 typedef typename Digraph::template NodeMap<int> HeapCrossRef; |
108 typedef BinHeap<Length, HeapCrossRef> Heap; |
108 typedef BinHeap<Length, HeapCrossRef> Heap; |
109 |
109 |
110 private: |
110 private: |
111 |
111 |
112 // The digraph the algorithm runs on |
|
113 const Digraph &_graph; |
112 const Digraph &_graph; |
114 |
113 const LengthMap &_length; |
115 // The main maps |
|
116 const FlowMap &_flow; |
114 const FlowMap &_flow; |
117 const LengthMap &_length; |
115 PotentialMap &_pi; |
118 PotentialMap &_potential; |
|
119 |
|
120 // The distance map |
|
121 PotentialMap _dist; |
|
122 // The pred arc map |
|
123 PredMap &_pred; |
116 PredMap &_pred; |
124 // The processed (i.e. permanently labeled) nodes |
|
125 std::vector<Node> _proc_nodes; |
|
126 |
|
127 Node _s; |
117 Node _s; |
128 Node _t; |
118 Node _t; |
|
119 |
|
120 PotentialMap _dist; |
|
121 std::vector<Node> _proc_nodes; |
129 |
122 |
130 public: |
123 public: |
131 |
124 |
132 /// Constructor. |
125 // Constructor |
133 ResidualDijkstra( const Digraph &graph, |
126 ResidualDijkstra(Suurballe &srb) : |
134 const FlowMap &flow, |
127 _graph(srb._graph), _length(srb._length), |
135 const LengthMap &length, |
128 _flow(*srb._flow), _pi(*srb._potential), _pred(srb._pred), |
136 PotentialMap &potential, |
129 _s(srb._s), _t(srb._t), _dist(_graph) {} |
137 PredMap &pred, |
130 |
138 Node s, Node t ) : |
131 // Run the algorithm and return true if a path is found |
139 _graph(graph), _flow(flow), _length(length), _potential(potential), |
132 // from the source node to the target node. |
140 _dist(graph), _pred(pred), _s(s), _t(t) {} |
133 bool run(int cnt) { |
141 |
134 return cnt == 0 ? startFirst() : start(); |
142 /// \brief Run the algorithm. It returns \c true if a path is found |
135 } |
143 /// from the source node to the target node. |
136 |
144 bool run() { |
137 private: |
|
138 |
|
139 // Execute the algorithm for the first time (the flow and potential |
|
140 // functions have to be identically zero). |
|
141 bool startFirst() { |
145 HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP); |
142 HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP); |
146 Heap heap(heap_cross_ref); |
143 Heap heap(heap_cross_ref); |
147 heap.push(_s, 0); |
144 heap.push(_s, 0); |
148 _pred[_s] = INVALID; |
145 _pred[_s] = INVALID; |
149 _proc_nodes.clear(); |
146 _proc_nodes.clear(); |
150 |
147 |
151 // Process nodes |
148 // Process nodes |
152 while (!heap.empty() && heap.top() != _t) { |
149 while (!heap.empty() && heap.top() != _t) { |
153 Node u = heap.top(), v; |
150 Node u = heap.top(), v; |
154 Length d = heap.prio() + _potential[u], nd; |
151 Length d = heap.prio(), dn; |
155 _dist[u] = heap.prio(); |
152 _dist[u] = heap.prio(); |
|
153 _proc_nodes.push_back(u); |
156 heap.pop(); |
154 heap.pop(); |
|
155 |
|
156 // Traverse outgoing arcs |
|
157 for (OutArcIt e(_graph, u); e != INVALID; ++e) { |
|
158 v = _graph.target(e); |
|
159 switch(heap.state(v)) { |
|
160 case Heap::PRE_HEAP: |
|
161 heap.push(v, d + _length[e]); |
|
162 _pred[v] = e; |
|
163 break; |
|
164 case Heap::IN_HEAP: |
|
165 dn = d + _length[e]; |
|
166 if (dn < heap[v]) { |
|
167 heap.decrease(v, dn); |
|
168 _pred[v] = e; |
|
169 } |
|
170 break; |
|
171 case Heap::POST_HEAP: |
|
172 break; |
|
173 } |
|
174 } |
|
175 } |
|
176 if (heap.empty()) return false; |
|
177 |
|
178 // Update potentials of processed nodes |
|
179 Length t_dist = heap.prio(); |
|
180 for (int i = 0; i < int(_proc_nodes.size()); ++i) |
|
181 _pi[_proc_nodes[i]] = _dist[_proc_nodes[i]] - t_dist; |
|
182 return true; |
|
183 } |
|
184 |
|
185 // Execute the algorithm. |
|
186 bool start() { |
|
187 HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP); |
|
188 Heap heap(heap_cross_ref); |
|
189 heap.push(_s, 0); |
|
190 _pred[_s] = INVALID; |
|
191 _proc_nodes.clear(); |
|
192 |
|
193 // Process nodes |
|
194 while (!heap.empty() && heap.top() != _t) { |
|
195 Node u = heap.top(), v; |
|
196 Length d = heap.prio() + _pi[u], dn; |
|
197 _dist[u] = heap.prio(); |
157 _proc_nodes.push_back(u); |
198 _proc_nodes.push_back(u); |
|
199 heap.pop(); |
158 |
200 |
159 // Traverse outgoing arcs |
201 // Traverse outgoing arcs |
160 for (OutArcIt e(_graph, u); e != INVALID; ++e) { |
202 for (OutArcIt e(_graph, u); e != INVALID; ++e) { |
161 if (_flow[e] == 0) { |
203 if (_flow[e] == 0) { |
162 v = _graph.target(e); |
204 v = _graph.target(e); |
163 switch(heap.state(v)) { |
205 switch(heap.state(v)) { |
164 case Heap::PRE_HEAP: |
206 case Heap::PRE_HEAP: |
165 heap.push(v, d + _length[e] - _potential[v]); |
207 heap.push(v, d + _length[e] - _pi[v]); |
166 _pred[v] = e; |
|
167 break; |
|
168 case Heap::IN_HEAP: |
|
169 nd = d + _length[e] - _potential[v]; |
|
170 if (nd < heap[v]) { |
|
171 heap.decrease(v, nd); |
|
172 _pred[v] = e; |
208 _pred[v] = e; |
173 } |
209 break; |
174 break; |
210 case Heap::IN_HEAP: |
175 case Heap::POST_HEAP: |
211 dn = d + _length[e] - _pi[v]; |
176 break; |
212 if (dn < heap[v]) { |
|
213 heap.decrease(v, dn); |
|
214 _pred[v] = e; |
|
215 } |
|
216 break; |
|
217 case Heap::POST_HEAP: |
|
218 break; |
177 } |
219 } |
178 } |
220 } |
179 } |
221 } |
180 |
222 |
181 // Traverse incoming arcs |
223 // Traverse incoming arcs |
182 for (InArcIt e(_graph, u); e != INVALID; ++e) { |
224 for (InArcIt e(_graph, u); e != INVALID; ++e) { |
183 if (_flow[e] == 1) { |
225 if (_flow[e] == 1) { |
184 v = _graph.source(e); |
226 v = _graph.source(e); |
185 switch(heap.state(v)) { |
227 switch(heap.state(v)) { |
186 case Heap::PRE_HEAP: |
228 case Heap::PRE_HEAP: |
187 heap.push(v, d - _length[e] - _potential[v]); |
229 heap.push(v, d - _length[e] - _pi[v]); |
188 _pred[v] = e; |
|
189 break; |
|
190 case Heap::IN_HEAP: |
|
191 nd = d - _length[e] - _potential[v]; |
|
192 if (nd < heap[v]) { |
|
193 heap.decrease(v, nd); |
|
194 _pred[v] = e; |
230 _pred[v] = e; |
195 } |
231 break; |
196 break; |
232 case Heap::IN_HEAP: |
197 case Heap::POST_HEAP: |
233 dn = d - _length[e] - _pi[v]; |
198 break; |
234 if (dn < heap[v]) { |
|
235 heap.decrease(v, dn); |
|
236 _pred[v] = e; |
|
237 } |
|
238 break; |
|
239 case Heap::POST_HEAP: |
|
240 break; |
199 } |
241 } |
200 } |
242 } |
201 } |
243 } |
202 } |
244 } |
203 if (heap.empty()) return false; |
245 if (heap.empty()) return false; |
204 |
246 |
205 // Update potentials of processed nodes |
247 // Update potentials of processed nodes |
206 Length t_dist = heap.prio(); |
248 Length t_dist = heap.prio(); |
207 for (int i = 0; i < int(_proc_nodes.size()); ++i) |
249 for (int i = 0; i < int(_proc_nodes.size()); ++i) |
208 _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist; |
250 _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist; |
209 return true; |
251 return true; |
210 } |
252 } |
211 |
253 |
212 }; //class ResidualDijkstra |
254 }; //class ResidualDijkstra |
213 |
255 |
370 /// the source node to the given node \c t in the digraph. |
408 /// the source node to the given node \c t in the digraph. |
371 /// Otherwise it returns the number of arc-disjoint paths found. |
409 /// Otherwise it returns the number of arc-disjoint paths found. |
372 /// |
410 /// |
373 /// \pre \ref init() must be called before using this function. |
411 /// \pre \ref init() must be called before using this function. |
374 int findFlow(const Node& t, int k = 2) { |
412 int findFlow(const Node& t, int k = 2) { |
375 _target = t; |
413 _t = t; |
376 _dijkstra = |
414 ResidualDijkstra dijkstra(*this); |
377 new ResidualDijkstra( _graph, *_flow, _length, *_potential, _pred, |
|
378 _source, _target ); |
|
379 |
415 |
380 // Find shortest paths |
416 // Find shortest paths |
381 _path_num = 0; |
417 _path_num = 0; |
382 while (_path_num < k) { |
418 while (_path_num < k) { |
383 // Run Dijkstra |
419 // Run Dijkstra |
384 if (!_dijkstra->run()) break; |
420 if (!dijkstra.run(_path_num)) break; |
385 ++_path_num; |
421 ++_path_num; |
386 |
422 |
387 // Set the flow along the found shortest path |
423 // Set the flow along the found shortest path |
388 Node u = _target; |
424 Node u = _t; |
389 Arc e; |
425 Arc e; |
390 while ((e = _pred[u]) != INVALID) { |
426 while ((e = _pred[u]) != INVALID) { |
391 if (u == _graph.target(e)) { |
427 if (u == _graph.target(e)) { |
392 (*_flow)[e] = 1; |
428 (*_flow)[e] = 1; |
393 u = _graph.source(e); |
429 u = _graph.source(e); |
400 return _path_num; |
436 return _path_num; |
401 } |
437 } |
402 |
438 |
403 /// \brief Compute the paths from the flow. |
439 /// \brief Compute the paths from the flow. |
404 /// |
440 /// |
405 /// This function computes the paths from the found minimum cost flow, |
441 /// This function computes arc-disjoint paths from the found minimum |
406 /// which is the union of some arc-disjoint paths. |
442 /// cost flow, which is the union of them. |
407 /// |
443 /// |
408 /// \pre \ref init() and \ref findFlow() must be called before using |
444 /// \pre \ref init() and \ref findFlow() must be called before using |
409 /// this function. |
445 /// this function. |
410 void findPaths() { |
446 void findPaths() { |
411 FlowMap res_flow(_graph); |
447 FlowMap res_flow(_graph); |
412 for(ArcIt a(_graph); a != INVALID; ++a) res_flow[a] = (*_flow)[a]; |
448 for(ArcIt a(_graph); a != INVALID; ++a) res_flow[a] = (*_flow)[a]; |
413 |
449 |
414 paths.clear(); |
450 _paths.clear(); |
415 paths.resize(_path_num); |
451 _paths.resize(_path_num); |
416 for (int i = 0; i < _path_num; ++i) { |
452 for (int i = 0; i < _path_num; ++i) { |
417 Node n = _source; |
453 Node n = _s; |
418 while (n != _target) { |
454 while (n != _t) { |
419 OutArcIt e(_graph, n); |
455 OutArcIt e(_graph, n); |
420 for ( ; res_flow[e] == 0; ++e) ; |
456 for ( ; res_flow[e] == 0; ++e) ; |
421 n = _graph.target(e); |
457 n = _graph.target(e); |
422 paths[i].addBack(e); |
458 _paths[i].addBack(e); |
423 res_flow[e] = 0; |
459 res_flow[e] = 0; |
424 } |
460 } |
425 } |
461 } |
426 } |
462 } |
427 |
463 |