Changeset 2586:37fb2c384c78 in lemon-0.x for lemon
- Timestamp:
- 02/29/08 16:55:13 (16 years ago)
- Branch:
- default
- Phase:
- public
- Convert:
- svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3468
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
lemon/suurballe.h
r2553 r2586 22 22 ///\ingroup shortest_path 23 23 ///\file 24 ///\brief An algorithm for finding k paths of minimal total length. 25 26 27 #include <lemon/maps.h> 24 ///\brief An algorithm for finding edge-disjoint paths between two 25 /// nodes having minimum total length. 26 28 27 #include <vector> 28 #include <lemon/bin_heap.h> 29 29 #include <lemon/path.h> 30 #include <lemon/ssp_min_cost_flow.h>31 30 32 31 namespace lemon { 33 32 34 /// \addtogroup shortest_path35 /// @{36 37 /// \brief Implementation of an algorithm for finding kedge-disjoint38 /// paths between 2 nodes of minimal total length39 /// 40 /// The class \ref lemon::Suurballe implements41 /// an algorithm for finding k edge-disjoint paths42 /// from a given source node to a given target node in a n43 /// edge-weighted directed graph having minimal total weight (length).44 /// 45 /// \warning Length values should be nonnegative!46 /// 47 /// \param Graph The directed graph type the algorithm runs on.48 /// \param LengthMap The type of the length map (values should be nonnegative).49 /// 50 /// \note It it questionable whether it is correct to call this method after51 /// %Suurballe for it is just a special case of Edmonds' and Karp's algorithm52 /// for finding minimum cost flows. In fact, this implementation just53 /// wraps the SspMinCostFlow algorithms. The paper of both %Suurballe and54 /// Edmonds-Karp published in 1972, therefore it is possibly right to55 /// state that they are56 /// independent results. Most frequently this special case is referred as57 ///%Suurballe method in the literature, especially in communication58 ///network context.59 ///\author Attila Bernath60 template <typename Graph, typename LengthMap>61 class Suurballe{62 33 /// \addtogroup shortest_path 34 /// @{ 35 36 /// \brief Implementation of an algorithm for finding edge-disjoint 37 /// paths between two nodes having minimum total length. 38 /// 39 /// \ref lemon::Suurballe "Suurballe" implements an algorithm for 40 /// finding edge-disjoint paths having minimum total length (cost) 41 /// from a given source node to a given target node in a directed 42 /// graph. 43 /// 44 /// In fact, this implementation is the specialization of the 45 /// \ref CapacityScaling "successive shortest path" algorithm. 46 /// 47 /// \tparam Graph The directed graph 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 SplitGraphAdaptor. 54 /// 55 /// \author Attila Bernath and Peter Kovacs 56 57 template < typename Graph, 58 typename LengthMap = typename Graph::template EdgeMap<int> > 59 class Suurballe 60 { 61 GRAPH_TYPEDEFS(typename Graph); 63 62 64 63 typedef typename LengthMap::Value Length; 64 typedef ConstMap<Edge, int> ConstEdgeMap; 65 typedef typename Graph::template NodeMap<Edge> PredMap; 66 67 public: 68 69 /// The type of the flow map. 70 typedef typename Graph::template EdgeMap<int> FlowMap; 71 /// The type of the potential map. 72 typedef typename Graph::template NodeMap<Length> PotentialMap; 73 /// The type of the path structures. 74 typedef SimplePath<Graph> 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 graph with respect to the reduced edge 84 /// lengths and modifying the node potentials according to the 85 /// distance of the nodes. 86 class ResidualDijkstra 87 { 88 typedef typename Graph::template NodeMap<int> HeapCrossRef; 89 typedef BinHeap<Length, HeapCrossRef> Heap; 90 91 private: 92 93 // The directed graph the algorithm runs on 94 const Graph &_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 edge 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 Graph &graph, 115 const FlowMap &flow, 116 const LengthMap &length, 117 PotentialMap &potential, 118 PredMap &pred, 119 Node s, Node t ) : 120 _graph(graph), _flow(flow), _length(length), _potential(potential), 121 _dist(graph), _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 edges 141 for (OutEdgeIt 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 edges 163 for (InEdgeIt 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 graph the algorithm runs on 198 const Graph &_graph; 199 // The length map 200 const LengthMap &_length; 65 201 66 typedef typename Graph::Node Node; 67 typedef typename Graph::NodeIt NodeIt; 68 typedef typename Graph::Edge Edge; 69 typedef typename Graph::OutEdgeIt OutEdgeIt; 70 typedef typename Graph::template EdgeMap<int> EdgeIntMap; 71 72 typedef ConstMap<Edge,int> ConstMap; 73 74 const Graph& G; 75 76 Node s; 77 Node t; 78 79 //Auxiliary variables 80 //This is the capacity map for the mincostflow problem 81 ConstMap const1map; 82 //This MinCostFlow instance will actually solve the problem 83 SspMinCostFlow<Graph, LengthMap, ConstMap> min_cost_flow; 84 85 //Container to store found paths 86 std::vector<SimplePath<Graph> > paths; 87 88 public : 89 90 91 /// \brief The constructor of the class. 92 /// 93 /// \param _G The directed graph the algorithm runs on. 94 /// \param _length The length (weight or cost) of the edges. 95 /// \param _s Source node. 96 /// \param _t Target node. 97 Suurballe(Graph& _G, LengthMap& _length, Node _s, Node _t) : 98 G(_G), s(_s), t(_t), const1map(1), 99 min_cost_flow(_G, _length, const1map, _s, _t) { } 202 // Edge 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<Graph> > paths; 216 int _path_num; 217 218 // The pred edge 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 graph The directed graph the algorithm runs on. 231 /// \param length The length (cost) values of the edges. 232 /// \param s The source node. 233 /// \param t The target node. 234 Suurballe( const Graph &graph, 235 const LengthMap &length, 236 Node s, Node t ) : 237 _graph(graph), _length(length), _flow(0), _local_flow(false), 238 _potential(0), _local_potential(false), _source(s), _target(t), 239 _pred(graph) {} 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 edge-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 /// edge-disjoint paths, you may call init() and findFlow(). 288 289 /// @{ 100 290 101 291 /// \brief Runs the algorithm. 102 292 /// 103 293 /// Runs the algorithm. 104 /// Returns k if there are at least k edge-disjoint paths from s to t. 105 /// Otherwise it returns the number of edge-disjoint paths found 106 /// from s to t. 107 /// 108 /// \param k How many paths are we looking for? 109 /// 110 int run(int k) { 111 int i = min_cost_flow.run(k); 112 113 //Let's find the paths 114 //We put the paths into stl vectors (as an inner representation). 115 //In the meantime we lose the information stored in 'reversed'. 116 //We suppose the lengths to be positive now. 117 118 //We don't want to change the flow of min_cost_flow, so we make a copy 119 //The name here suggests that the flow has only 0/1 values. 120 EdgeIntMap reversed(G); 121 122 for(typename Graph::EdgeIt e(G); e!=INVALID; ++e) 123 reversed[e] = min_cost_flow.getFlow()[e]; 124 294 /// 295 /// \param k The number of paths to be found. 296 /// 297 /// \return \c k if there are at least \c k edge-disjoint paths 298 /// from \c s to \c t. Otherwise it returns the number of 299 /// edge-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 (EdgeIt 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 /// edge-disjoint paths. 342 /// 343 /// \return \c k if there are at least \c k edge-disjoint paths 344 /// from \c s to \c t. Otherwise it returns the number of 345 /// edge-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 Edge 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(*_flow); 382 125 383 paths.clear(); 126 paths.resize(k); 127 for (int j=0; j<i; ++j){ 128 Node n=s; 129 130 while (n!=t){ 131 132 OutEdgeIt e(G, n); 133 134 while (!reversed[e]){ 135 ++e; 136 } 137 n = G.target(e); 138 paths[j].addBack(e); 139 reversed[e] = 1-reversed[e]; 140 } 141 142 } 143 return i; 144 } 145 146 147 /// \brief Returns the total length of the paths. 148 /// 149 /// This function gives back the total length of the found paths. 150 Length totalLength(){ 151 return min_cost_flow.totalLength(); 152 } 153 154 /// \brief Returns the found flow. 155 /// 156 /// This function returns a const reference to the EdgeMap \c flow. 157 const EdgeIntMap &getFlow() const { return min_cost_flow.flow;} 158 159 /// \brief Returns the optimal dual solution 160 /// 161 /// This function returns a const reference to the NodeMap \c 162 /// potential (the dual solution). 163 const EdgeIntMap &getPotential() const { return min_cost_flow.potential;} 164 165 /// \brief Checks whether the complementary slackness holds. 166 /// 167 /// This function checks, whether the given solution is optimal. 168 /// Currently this function only checks optimality, doesn't bother 169 /// with feasibility. It is meant for testing purposes. 170 bool checkComplementarySlackness(){ 171 return min_cost_flow.checkComplementarySlackness(); 172 } 173 174 typedef SimplePath<Graph> Path; 175 176 /// \brief Read the found paths. 177 /// 178 /// This function gives back the \c j-th path in argument p. 179 /// Assumes that \c run() has been run and nothing has changed 180 /// since then. 181 /// 182 /// \warning It is assumed that \c p is constructed to be a path 183 /// of graph \c G. If \c j is not less than the result of 184 /// previous \c run, then the result here will be an empty path 185 /// (\c j can be 0 as well). 186 /// 187 /// \param j Which path you want to get from the found paths (in a 188 /// real application you would get the found paths iteratively). 189 Path path(int j) const { 190 return paths[j]; 191 } 192 193 /// \brief Gives back the number of the paths. 194 /// 195 /// Gives back the number of the constructed paths. 384 paths.resize(_path_num); 385 for (int i = 0; i < _path_num; ++i) { 386 Node n = _source; 387 while (n != _target) { 388 OutEdgeIt e(_graph, n); 389 for ( ; res_flow[e] == 0; ++e) ; 390 n = _graph.target(e); 391 paths[i].addBack(e); 392 res_flow[e] = 0; 393 } 394 } 395 } 396 397 /// @} 398 399 /// \name Query Functions 400 /// The result of the algorithm can be obtained using these 401 /// functions. 402 /// \n The algorithm should be executed before using them. 403 404 /// @{ 405 406 /// \brief Returns a const reference to the edge map storing the 407 /// found flow. 408 /// 409 /// Returns a const reference to the edge map storing the flow that 410 /// is the union of the found edge-disjoint paths. 411 /// 412 /// \pre \ref run() or findFlow() must be called before using this 413 /// function. 414 const FlowMap& flowMap() const { 415 return *_flow; 416 } 417 418 /// \brief Returns a const reference to the node map storing the 419 /// found potentials (the dual solution). 420 /// 421 /// Returns a const reference to the node map storing the found 422 /// potentials that provide the dual solution of the underlying 423 /// minimum cost flow problem. 424 /// 425 /// \pre \ref run() or findFlow() must be called before using this 426 /// function. 427 const PotentialMap& potentialMap() const { 428 return *_potential; 429 } 430 431 /// \brief Returns the flow on the given edge. 432 /// 433 /// Returns the flow on the given edge. 434 /// It is \c 1 if the edge is involved in one of the found paths, 435 /// otherwise it is \c 0. 436 /// 437 /// \pre \ref run() or findFlow() must be called before using this 438 /// function. 439 int flow(const Edge& edge) const { 440 return (*_flow)[edge]; 441 } 442 443 /// \brief Returns the potential of the given node. 444 /// 445 /// Returns the potential of the given node. 446 /// 447 /// \pre \ref run() or findFlow() must be called before using this 448 /// function. 449 Length potential(const Node& node) const { 450 return (*_potential)[node]; 451 } 452 453 /// \brief Returns the total length (cost) of the found paths (flow). 454 /// 455 /// Returns the total length (cost) of the found paths (flow). 456 /// The complexity of the function is \f$ O(e) \f$. 457 /// 458 /// \pre \ref run() or findFlow() must be called before using this 459 /// function. 460 Length totalLength() const { 461 Length c = 0; 462 for (EdgeIt e(_graph); e != INVALID; ++e) 463 c += (*_flow)[e] * _length[e]; 464 return c; 465 } 466 467 /// \brief Returns the number of the found paths. 468 /// 469 /// Returns the number of the found paths. 470 /// 471 /// \pre \ref run() or findFlow() must be called before using this 472 /// function. 196 473 int pathNum() const { 197 return paths.size(); 198 } 474 return _path_num; 475 } 476 477 /// \brief Returns a const reference to the specified path. 478 /// 479 /// Returns a const reference to the specified path. 480 /// 481 /// \param i The function returns the \c i-th path. 482 /// \c i must be between \c 0 and <tt>%pathNum()-1</tt>. 483 /// 484 /// \pre \ref run() or findPaths() must be called before using this 485 /// function. 486 Path path(int i) const { 487 return paths[i]; 488 } 489 490 /// @} 199 491 200 492 }; //class Suurballe
Note: See TracChangeset
for help on using the changeset viewer.