diff -r 30c77d1c0cba -r ec0b1b423b8b lemon/suurballe.h --- a/lemon/suurballe.h Thu Oct 15 21:04:50 2009 +0200 +++ b/lemon/suurballe.h Fri Oct 16 01:06:16 2009 +0200 @@ -46,7 +46,7 @@ /// Note that this problem is a special case of the \ref min_cost_flow /// "minimum cost flow problem". This implementation is actually an /// efficient specialized version of the \ref CapacityScaling - /// "Successive Shortest Path" algorithm directly for this problem. + /// "successive shortest path" algorithm directly for this problem. /// Therefore this class provides query functions for flow values and /// node potentials (the dual solution) just like the minimum cost flow /// algorithms. @@ -57,7 +57,7 @@ /// /// \warning Length values should be \e non-negative. /// - /// \note For finding node-disjoint paths this algorithm can be used + /// \note For finding \e node-disjoint paths, this algorithm can be used /// along with the \ref SplitNodes adaptor. #ifdef DOXYGEN template @@ -109,39 +109,36 @@ private: - // The digraph the algorithm runs on const Digraph &_graph; - - // The main maps + const LengthMap &_length; const FlowMap &_flow; - const LengthMap &_length; - PotentialMap &_potential; - - // The distance map - PotentialMap _dist; - // The pred arc map + PotentialMap &_pi; PredMap &_pred; - // The processed (i.e. permanently labeled) nodes - std::vector _proc_nodes; - Node _s; Node _t; + + PotentialMap _dist; + std::vector _proc_nodes; public: - /// Constructor. - ResidualDijkstra( const Digraph &graph, - const FlowMap &flow, - const LengthMap &length, - PotentialMap &potential, - PredMap &pred, - Node s, Node t ) : - _graph(graph), _flow(flow), _length(length), _potential(potential), - _dist(graph), _pred(pred), _s(s), _t(t) {} + // Constructor + ResidualDijkstra(Suurballe &srb) : + _graph(srb._graph), _length(srb._length), + _flow(*srb._flow), _pi(*srb._potential), _pred(srb._pred), + _s(srb._s), _t(srb._t), _dist(_graph) {} + + // Run the algorithm and return true if a path is found + // from the source node to the target node. + bool run(int cnt) { + return cnt == 0 ? startFirst() : start(); + } - /// \brief Run the algorithm. It returns \c true if a path is found - /// from the source node to the target node. - bool run() { + private: + + // Execute the algorithm for the first time (the flow and potential + // functions have to be identically zero). + bool startFirst() { HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP); Heap heap(heap_cross_ref); heap.push(_s, 0); @@ -151,29 +148,74 @@ // Process nodes while (!heap.empty() && heap.top() != _t) { Node u = heap.top(), v; - Length d = heap.prio() + _potential[u], nd; + Length d = heap.prio(), dn; _dist[u] = heap.prio(); + _proc_nodes.push_back(u); heap.pop(); + + // Traverse outgoing arcs + for (OutArcIt e(_graph, u); e != INVALID; ++e) { + v = _graph.target(e); + switch(heap.state(v)) { + case Heap::PRE_HEAP: + heap.push(v, d + _length[e]); + _pred[v] = e; + break; + case Heap::IN_HEAP: + dn = d + _length[e]; + if (dn < heap[v]) { + heap.decrease(v, dn); + _pred[v] = e; + } + break; + case Heap::POST_HEAP: + break; + } + } + } + if (heap.empty()) return false; + + // Update potentials of processed nodes + Length t_dist = heap.prio(); + for (int i = 0; i < int(_proc_nodes.size()); ++i) + _pi[_proc_nodes[i]] = _dist[_proc_nodes[i]] - t_dist; + return true; + } + + // Execute the algorithm. + bool start() { + HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP); + Heap heap(heap_cross_ref); + heap.push(_s, 0); + _pred[_s] = INVALID; + _proc_nodes.clear(); + + // Process nodes + while (!heap.empty() && heap.top() != _t) { + Node u = heap.top(), v; + Length d = heap.prio() + _pi[u], dn; + _dist[u] = heap.prio(); _proc_nodes.push_back(u); + heap.pop(); // Traverse outgoing arcs for (OutArcIt e(_graph, u); e != INVALID; ++e) { if (_flow[e] == 0) { v = _graph.target(e); switch(heap.state(v)) { - case Heap::PRE_HEAP: - heap.push(v, d + _length[e] - _potential[v]); - _pred[v] = e; - break; - case Heap::IN_HEAP: - nd = d + _length[e] - _potential[v]; - if (nd < heap[v]) { - heap.decrease(v, nd); + case Heap::PRE_HEAP: + heap.push(v, d + _length[e] - _pi[v]); _pred[v] = e; - } - break; - case Heap::POST_HEAP: - break; + break; + case Heap::IN_HEAP: + dn = d + _length[e] - _pi[v]; + if (dn < heap[v]) { + heap.decrease(v, dn); + _pred[v] = e; + } + break; + case Heap::POST_HEAP: + break; } } } @@ -183,19 +225,19 @@ if (_flow[e] == 1) { v = _graph.source(e); switch(heap.state(v)) { - case Heap::PRE_HEAP: - heap.push(v, d - _length[e] - _potential[v]); - _pred[v] = e; - break; - case Heap::IN_HEAP: - nd = d - _length[e] - _potential[v]; - if (nd < heap[v]) { - heap.decrease(v, nd); + case Heap::PRE_HEAP: + heap.push(v, d - _length[e] - _pi[v]); _pred[v] = e; - } - break; - case Heap::POST_HEAP: - break; + break; + case Heap::IN_HEAP: + dn = d - _length[e] - _pi[v]; + if (dn < heap[v]) { + heap.decrease(v, dn); + _pred[v] = e; + } + break; + case Heap::POST_HEAP: + break; } } } @@ -205,7 +247,7 @@ // Update potentials of processed nodes Length t_dist = heap.prio(); for (int i = 0; i < int(_proc_nodes.size()); ++i) - _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist; + _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist; return true; } @@ -226,19 +268,16 @@ bool _local_potential; // The source node - Node _source; + Node _s; // The target node - Node _target; + Node _t; // Container to store the found paths - std::vector< SimplePath > paths; + std::vector _paths; int _path_num; // The pred arc map PredMap _pred; - // Implementation of the Dijkstra algorithm for finding augmenting - // shortest paths in the residual network - ResidualDijkstra *_dijkstra; public: @@ -258,7 +297,6 @@ ~Suurballe() { if (_local_flow) delete _flow; if (_local_potential) delete _potential; - delete _dijkstra; } /// \brief Set the flow map. @@ -342,7 +380,7 @@ /// /// \param s The source node. void init(const Node& s) { - _source = s; + _s = s; // Initialize maps if (!_flow) { @@ -372,20 +410,18 @@ /// /// \pre \ref init() must be called before using this function. int findFlow(const Node& t, int k = 2) { - _target = t; - _dijkstra = - new ResidualDijkstra( _graph, *_flow, _length, *_potential, _pred, - _source, _target ); + _t = t; + ResidualDijkstra dijkstra(*this); // Find shortest paths _path_num = 0; while (_path_num < k) { // Run Dijkstra - if (!_dijkstra->run()) break; + if (!dijkstra.run(_path_num)) break; ++_path_num; // Set the flow along the found shortest path - Node u = _target; + Node u = _t; Arc e; while ((e = _pred[u]) != INVALID) { if (u == _graph.target(e)) { @@ -402,8 +438,8 @@ /// \brief Compute the paths from the flow. /// - /// This function computes the paths from the found minimum cost flow, - /// which is the union of some arc-disjoint paths. + /// This function computes arc-disjoint paths from the found minimum + /// cost flow, which is the union of them. /// /// \pre \ref init() and \ref findFlow() must be called before using /// this function. @@ -411,15 +447,15 @@ FlowMap res_flow(_graph); for(ArcIt a(_graph); a != INVALID; ++a) res_flow[a] = (*_flow)[a]; - paths.clear(); - paths.resize(_path_num); + _paths.clear(); + _paths.resize(_path_num); for (int i = 0; i < _path_num; ++i) { - Node n = _source; - while (n != _target) { + Node n = _s; + while (n != _t) { OutArcIt e(_graph, n); for ( ; res_flow[e] == 0; ++e) ; n = _graph.target(e); - paths[i].addBack(e); + _paths[i].addBack(e); res_flow[e] = 0; } } @@ -518,7 +554,7 @@ /// \pre \ref run() or \ref findPaths() must be called before using /// this function. const Path& path(int i) const { - return paths[i]; + return _paths[i]; } /// @}