|
1 /* -*- C++ -*- |
|
2 * |
|
3 * This file is a part of LEMON, a generic C++ optimization library |
|
4 * |
|
5 * Copyright (C) 2003-2008 |
|
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport |
|
7 * (Egervary Research Group on Combinatorial Optimization, EGRES). |
|
8 * |
|
9 * Permission to use, modify and distribute this software is granted |
|
10 * provided that this copyright notice appears in all copies. For |
|
11 * precise terms see the accompanying LICENSE file. |
|
12 * |
|
13 * This software is provided "AS IS" with no warranty of any kind, |
|
14 * express or implied, and with no claim as to its suitability for any |
|
15 * purpose. |
|
16 * |
|
17 */ |
|
18 |
|
19 #ifndef LEMON_SUURBALLE_H |
|
20 #define LEMON_SUURBALLE_H |
|
21 |
|
22 ///\ingroup shortest_path |
|
23 ///\file |
|
24 ///\brief An algorithm for finding arc-disjoint paths between two |
|
25 /// nodes having minimum total length. |
|
26 |
|
27 #include <vector> |
|
28 #include <lemon/bin_heap.h> |
|
29 #include <lemon/path.h> |
|
30 |
|
31 namespace lemon { |
|
32 |
|
33 /// \addtogroup shortest_path |
|
34 /// @{ |
|
35 |
|
36 /// \brief Implementation of an algorithm for finding arc-disjoint |
|
37 /// paths between two nodes having minimum total length. |
|
38 /// |
|
39 /// \ref lemon::Suurballe "Suurballe" implements an algorithm for |
|
40 /// finding arc-disjoint paths having minimum total length (cost) |
|
41 /// from a given source node to a given target node in a directed |
|
42 /// digraph. |
|
43 /// |
|
44 /// In fact, this implementation is the specialization of the |
|
45 /// \ref CapacityScaling "successive shortest path" algorithm. |
|
46 /// |
|
47 /// \tparam Digraph The directed digraph type the algorithm runs on. |
|
48 /// \tparam LengthMap The type of the length (cost) map. |
|
49 /// |
|
50 /// \warning Length values should be \e non-negative \e integers. |
|
51 /// |
|
52 /// \note For finding node-disjoint paths this algorithm can be used |
|
53 /// with \ref SplitDigraphAdaptor. |
|
54 /// |
|
55 /// \author Attila Bernath and Peter Kovacs |
|
56 |
|
57 template < typename Digraph, |
|
58 typename LengthMap = typename Digraph::template ArcMap<int> > |
|
59 class Suurballe |
|
60 { |
|
61 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); |
|
62 |
|
63 typedef typename LengthMap::Value Length; |
|
64 typedef ConstMap<Arc, int> ConstArcMap; |
|
65 typedef typename Digraph::template NodeMap<Arc> PredMap; |
|
66 |
|
67 public: |
|
68 |
|
69 /// The type of the flow map. |
|
70 typedef typename Digraph::template ArcMap<int> FlowMap; |
|
71 /// The type of the potential map. |
|
72 typedef typename Digraph::template NodeMap<Length> PotentialMap; |
|
73 /// The type of the path structures. |
|
74 typedef SimplePath<Digraph> Path; |
|
75 |
|
76 private: |
|
77 |
|
78 /// \brief Special implementation of the \ref Dijkstra algorithm |
|
79 /// for finding shortest paths in the residual network. |
|
80 /// |
|
81 /// \ref ResidualDijkstra is a special implementation of the |
|
82 /// \ref Dijkstra algorithm for finding shortest paths in the |
|
83 /// residual network of the digraph with respect to the reduced arc |
|
84 /// lengths and modifying the node potentials according to the |
|
85 /// distance of the nodes. |
|
86 class ResidualDijkstra |
|
87 { |
|
88 typedef typename Digraph::template NodeMap<int> HeapCrossRef; |
|
89 typedef BinHeap<Length, HeapCrossRef> Heap; |
|
90 |
|
91 private: |
|
92 |
|
93 // The directed digraph the algorithm runs on |
|
94 const Digraph &_graph; |
|
95 |
|
96 // The main maps |
|
97 const FlowMap &_flow; |
|
98 const LengthMap &_length; |
|
99 PotentialMap &_potential; |
|
100 |
|
101 // The distance map |
|
102 PotentialMap _dist; |
|
103 // The pred arc map |
|
104 PredMap &_pred; |
|
105 // The processed (i.e. permanently labeled) nodes |
|
106 std::vector<Node> _proc_nodes; |
|
107 |
|
108 Node _s; |
|
109 Node _t; |
|
110 |
|
111 public: |
|
112 |
|
113 /// Constructor. |
|
114 ResidualDijkstra( const Digraph &digraph, |
|
115 const FlowMap &flow, |
|
116 const LengthMap &length, |
|
117 PotentialMap &potential, |
|
118 PredMap &pred, |
|
119 Node s, Node t ) : |
|
120 _graph(digraph), _flow(flow), _length(length), _potential(potential), |
|
121 _dist(digraph), _pred(pred), _s(s), _t(t) {} |
|
122 |
|
123 /// \brief Runs the algorithm. Returns \c true if a path is found |
|
124 /// from the source node to the target node. |
|
125 bool run() { |
|
126 HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP); |
|
127 Heap heap(heap_cross_ref); |
|
128 heap.push(_s, 0); |
|
129 _pred[_s] = INVALID; |
|
130 _proc_nodes.clear(); |
|
131 |
|
132 // Processing nodes |
|
133 while (!heap.empty() && heap.top() != _t) { |
|
134 Node u = heap.top(), v; |
|
135 Length d = heap.prio() + _potential[u], nd; |
|
136 _dist[u] = heap.prio(); |
|
137 heap.pop(); |
|
138 _proc_nodes.push_back(u); |
|
139 |
|
140 // Traversing outgoing arcs |
|
141 for (OutArcIt e(_graph, u); e != INVALID; ++e) { |
|
142 if (_flow[e] == 0) { |
|
143 v = _graph.target(e); |
|
144 switch(heap.state(v)) { |
|
145 case Heap::PRE_HEAP: |
|
146 heap.push(v, d + _length[e] - _potential[v]); |
|
147 _pred[v] = e; |
|
148 break; |
|
149 case Heap::IN_HEAP: |
|
150 nd = d + _length[e] - _potential[v]; |
|
151 if (nd < heap[v]) { |
|
152 heap.decrease(v, nd); |
|
153 _pred[v] = e; |
|
154 } |
|
155 break; |
|
156 case Heap::POST_HEAP: |
|
157 break; |
|
158 } |
|
159 } |
|
160 } |
|
161 |
|
162 // Traversing incoming arcs |
|
163 for (InArcIt e(_graph, u); e != INVALID; ++e) { |
|
164 if (_flow[e] == 1) { |
|
165 v = _graph.source(e); |
|
166 switch(heap.state(v)) { |
|
167 case Heap::PRE_HEAP: |
|
168 heap.push(v, d - _length[e] - _potential[v]); |
|
169 _pred[v] = e; |
|
170 break; |
|
171 case Heap::IN_HEAP: |
|
172 nd = d - _length[e] - _potential[v]; |
|
173 if (nd < heap[v]) { |
|
174 heap.decrease(v, nd); |
|
175 _pred[v] = e; |
|
176 } |
|
177 break; |
|
178 case Heap::POST_HEAP: |
|
179 break; |
|
180 } |
|
181 } |
|
182 } |
|
183 } |
|
184 if (heap.empty()) return false; |
|
185 |
|
186 // Updating potentials of processed nodes |
|
187 Length t_dist = heap.prio(); |
|
188 for (int i = 0; i < int(_proc_nodes.size()); ++i) |
|
189 _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist; |
|
190 return true; |
|
191 } |
|
192 |
|
193 }; //class ResidualDijkstra |
|
194 |
|
195 private: |
|
196 |
|
197 // The directed digraph the algorithm runs on |
|
198 const Digraph &_graph; |
|
199 // The length map |
|
200 const LengthMap &_length; |
|
201 |
|
202 // Arc map of the current flow |
|
203 FlowMap *_flow; |
|
204 bool _local_flow; |
|
205 // Node map of the current potentials |
|
206 PotentialMap *_potential; |
|
207 bool _local_potential; |
|
208 |
|
209 // The source node |
|
210 Node _source; |
|
211 // The target node |
|
212 Node _target; |
|
213 |
|
214 // Container to store the found paths |
|
215 std::vector< SimplePath<Digraph> > paths; |
|
216 int _path_num; |
|
217 |
|
218 // The pred arc map |
|
219 PredMap _pred; |
|
220 // Implementation of the Dijkstra algorithm for finding augmenting |
|
221 // shortest paths in the residual network |
|
222 ResidualDijkstra *_dijkstra; |
|
223 |
|
224 public: |
|
225 |
|
226 /// \brief Constructor. |
|
227 /// |
|
228 /// Constructor. |
|
229 /// |
|
230 /// \param digraph The directed digraph the algorithm runs on. |
|
231 /// \param length The length (cost) values of the arcs. |
|
232 /// \param s The source node. |
|
233 /// \param t The target node. |
|
234 Suurballe( const Digraph &digraph, |
|
235 const LengthMap &length, |
|
236 Node s, Node t ) : |
|
237 _graph(digraph), _length(length), _flow(0), _local_flow(false), |
|
238 _potential(0), _local_potential(false), _source(s), _target(t), |
|
239 _pred(digraph) {} |
|
240 |
|
241 /// Destructor. |
|
242 ~Suurballe() { |
|
243 if (_local_flow) delete _flow; |
|
244 if (_local_potential) delete _potential; |
|
245 delete _dijkstra; |
|
246 } |
|
247 |
|
248 /// \brief Sets the flow map. |
|
249 /// |
|
250 /// Sets the flow map. |
|
251 /// |
|
252 /// The found flow contains only 0 and 1 values. It is the union of |
|
253 /// the found arc-disjoint paths. |
|
254 /// |
|
255 /// \return \c (*this) |
|
256 Suurballe& flowMap(FlowMap &map) { |
|
257 if (_local_flow) { |
|
258 delete _flow; |
|
259 _local_flow = false; |
|
260 } |
|
261 _flow = ↦ |
|
262 return *this; |
|
263 } |
|
264 |
|
265 /// \brief Sets the potential map. |
|
266 /// |
|
267 /// Sets the potential map. |
|
268 /// |
|
269 /// The potentials provide the dual solution of the underlying |
|
270 /// minimum cost flow problem. |
|
271 /// |
|
272 /// \return \c (*this) |
|
273 Suurballe& potentialMap(PotentialMap &map) { |
|
274 if (_local_potential) { |
|
275 delete _potential; |
|
276 _local_potential = false; |
|
277 } |
|
278 _potential = ↦ |
|
279 return *this; |
|
280 } |
|
281 |
|
282 /// \name Execution control |
|
283 /// The simplest way to execute the algorithm is to call the run() |
|
284 /// function. |
|
285 /// \n |
|
286 /// If you only need the flow that is the union of the found |
|
287 /// arc-disjoint paths, you may call init() and findFlow(). |
|
288 |
|
289 /// @{ |
|
290 |
|
291 /// \brief Runs the algorithm. |
|
292 /// |
|
293 /// Runs the algorithm. |
|
294 /// |
|
295 /// \param k The number of paths to be found. |
|
296 /// |
|
297 /// \return \c k if there are at least \c k arc-disjoint paths |
|
298 /// from \c s to \c t. Otherwise it returns the number of |
|
299 /// arc-disjoint paths found. |
|
300 /// |
|
301 /// \note Apart from the return value, <tt>s.run(k)</tt> is just a |
|
302 /// shortcut of the following code. |
|
303 /// \code |
|
304 /// s.init(); |
|
305 /// s.findFlow(k); |
|
306 /// s.findPaths(); |
|
307 /// \endcode |
|
308 int run(int k = 2) { |
|
309 init(); |
|
310 findFlow(k); |
|
311 findPaths(); |
|
312 return _path_num; |
|
313 } |
|
314 |
|
315 /// \brief Initializes the algorithm. |
|
316 /// |
|
317 /// Initializes the algorithm. |
|
318 void init() { |
|
319 // Initializing maps |
|
320 if (!_flow) { |
|
321 _flow = new FlowMap(_graph); |
|
322 _local_flow = true; |
|
323 } |
|
324 if (!_potential) { |
|
325 _potential = new PotentialMap(_graph); |
|
326 _local_potential = true; |
|
327 } |
|
328 for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0; |
|
329 for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0; |
|
330 |
|
331 _dijkstra = new ResidualDijkstra( _graph, *_flow, _length, |
|
332 *_potential, _pred, |
|
333 _source, _target ); |
|
334 } |
|
335 |
|
336 /// \brief Executes the successive shortest path algorithm to find |
|
337 /// an optimal flow. |
|
338 /// |
|
339 /// Executes the successive shortest path algorithm to find a |
|
340 /// minimum cost flow, which is the union of \c k or less |
|
341 /// arc-disjoint paths. |
|
342 /// |
|
343 /// \return \c k if there are at least \c k arc-disjoint paths |
|
344 /// from \c s to \c t. Otherwise it returns the number of |
|
345 /// arc-disjoint paths found. |
|
346 /// |
|
347 /// \pre \ref init() must be called before using this function. |
|
348 int findFlow(int k = 2) { |
|
349 // Finding shortest paths |
|
350 _path_num = 0; |
|
351 while (_path_num < k) { |
|
352 // Running Dijkstra |
|
353 if (!_dijkstra->run()) break; |
|
354 ++_path_num; |
|
355 |
|
356 // Setting the flow along the found shortest path |
|
357 Node u = _target; |
|
358 Arc e; |
|
359 while ((e = _pred[u]) != INVALID) { |
|
360 if (u == _graph.target(e)) { |
|
361 (*_flow)[e] = 1; |
|
362 u = _graph.source(e); |
|
363 } else { |
|
364 (*_flow)[e] = 0; |
|
365 u = _graph.target(e); |
|
366 } |
|
367 } |
|
368 } |
|
369 return _path_num; |
|
370 } |
|
371 |
|
372 /// \brief Computes the paths from the flow. |
|
373 /// |
|
374 /// Computes the paths from the flow. |
|
375 /// |
|
376 /// \pre \ref init() and \ref findFlow() must be called before using |
|
377 /// this function. |
|
378 void findPaths() { |
|
379 // Creating the residual flow map (the union of the paths not |
|
380 // found so far) |
|
381 FlowMap res_flow(_graph); |
|
382 for(ArcIt a(_graph);a!=INVALID;++a) res_flow[a]=(*_flow)[a]; |
|
383 |
|
384 paths.clear(); |
|
385 paths.resize(_path_num); |
|
386 for (int i = 0; i < _path_num; ++i) { |
|
387 Node n = _source; |
|
388 while (n != _target) { |
|
389 OutArcIt e(_graph, n); |
|
390 for ( ; res_flow[e] == 0; ++e) ; |
|
391 n = _graph.target(e); |
|
392 paths[i].addBack(e); |
|
393 res_flow[e] = 0; |
|
394 } |
|
395 } |
|
396 } |
|
397 |
|
398 /// @} |
|
399 |
|
400 /// \name Query Functions |
|
401 /// The result of the algorithm can be obtained using these |
|
402 /// functions. |
|
403 /// \n The algorithm should be executed before using them. |
|
404 |
|
405 /// @{ |
|
406 |
|
407 /// \brief Returns a const reference to the arc map storing the |
|
408 /// found flow. |
|
409 /// |
|
410 /// Returns a const reference to the arc map storing the flow that |
|
411 /// is the union of the found arc-disjoint paths. |
|
412 /// |
|
413 /// \pre \ref run() or findFlow() must be called before using this |
|
414 /// function. |
|
415 const FlowMap& flowMap() const { |
|
416 return *_flow; |
|
417 } |
|
418 |
|
419 /// \brief Returns a const reference to the node map storing the |
|
420 /// found potentials (the dual solution). |
|
421 /// |
|
422 /// Returns a const reference to the node map storing the found |
|
423 /// potentials that provide the dual solution of the underlying |
|
424 /// minimum cost flow problem. |
|
425 /// |
|
426 /// \pre \ref run() or findFlow() must be called before using this |
|
427 /// function. |
|
428 const PotentialMap& potentialMap() const { |
|
429 return *_potential; |
|
430 } |
|
431 |
|
432 /// \brief Returns the flow on the given arc. |
|
433 /// |
|
434 /// Returns the flow on the given arc. |
|
435 /// It is \c 1 if the arc is involved in one of the found paths, |
|
436 /// otherwise it is \c 0. |
|
437 /// |
|
438 /// \pre \ref run() or findFlow() must be called before using this |
|
439 /// function. |
|
440 int flow(const Arc& arc) const { |
|
441 return (*_flow)[arc]; |
|
442 } |
|
443 |
|
444 /// \brief Returns the potential of the given node. |
|
445 /// |
|
446 /// Returns the potential of the given node. |
|
447 /// |
|
448 /// \pre \ref run() or findFlow() must be called before using this |
|
449 /// function. |
|
450 Length potential(const Node& node) const { |
|
451 return (*_potential)[node]; |
|
452 } |
|
453 |
|
454 /// \brief Returns the total length (cost) of the found paths (flow). |
|
455 /// |
|
456 /// Returns the total length (cost) of the found paths (flow). |
|
457 /// The complexity of the function is \f$ O(e) \f$. |
|
458 /// |
|
459 /// \pre \ref run() or findFlow() must be called before using this |
|
460 /// function. |
|
461 Length totalLength() const { |
|
462 Length c = 0; |
|
463 for (ArcIt e(_graph); e != INVALID; ++e) |
|
464 c += (*_flow)[e] * _length[e]; |
|
465 return c; |
|
466 } |
|
467 |
|
468 /// \brief Returns the number of the found paths. |
|
469 /// |
|
470 /// Returns the number of the found paths. |
|
471 /// |
|
472 /// \pre \ref run() or findFlow() must be called before using this |
|
473 /// function. |
|
474 int pathNum() const { |
|
475 return _path_num; |
|
476 } |
|
477 |
|
478 /// \brief Returns a const reference to the specified path. |
|
479 /// |
|
480 /// Returns a const reference to the specified path. |
|
481 /// |
|
482 /// \param i The function returns the \c i-th path. |
|
483 /// \c i must be between \c 0 and <tt>%pathNum()-1</tt>. |
|
484 /// |
|
485 /// \pre \ref run() or findPaths() must be called before using this |
|
486 /// function. |
|
487 Path path(int i) const { |
|
488 return paths[i]; |
|
489 } |
|
490 |
|
491 /// @} |
|
492 |
|
493 }; //class Suurballe |
|
494 |
|
495 ///@} |
|
496 |
|
497 } //namespace lemon |
|
498 |
|
499 #endif //LEMON_SUURBALLE_H |