alpar@440: /* -*- mode: C++; indent-tabs-mode: nil; -*- alpar@345: * alpar@440: * This file is a part of LEMON, a generic C++ optimization library. alpar@345: * alpar@440: * Copyright (C) 2003-2009 alpar@345: * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport alpar@345: * (Egervary Research Group on Combinatorial Optimization, EGRES). alpar@345: * alpar@345: * Permission to use, modify and distribute this software is granted alpar@345: * provided that this copyright notice appears in all copies. For alpar@345: * precise terms see the accompanying LICENSE file. alpar@345: * alpar@345: * This software is provided "AS IS" with no warranty of any kind, alpar@345: * express or implied, and with no claim as to its suitability for any alpar@345: * purpose. alpar@345: * alpar@345: */ alpar@345: alpar@345: #ifndef LEMON_SUURBALLE_H alpar@345: #define LEMON_SUURBALLE_H alpar@345: alpar@345: ///\ingroup shortest_path alpar@345: ///\file alpar@345: ///\brief An algorithm for finding arc-disjoint paths between two alpar@345: /// nodes having minimum total length. alpar@345: alpar@345: #include kpeter@615: #include alpar@345: #include alpar@345: #include deba@503: #include deba@503: #include alpar@345: alpar@345: namespace lemon { alpar@345: alpar@345: /// \addtogroup shortest_path alpar@345: /// @{ alpar@345: kpeter@346: /// \brief Algorithm for finding arc-disjoint paths between two nodes kpeter@346: /// having minimum total length. alpar@345: /// alpar@345: /// \ref lemon::Suurballe "Suurballe" implements an algorithm for alpar@345: /// finding arc-disjoint paths having minimum total length (cost) kpeter@346: /// from a given source node to a given target node in a digraph. alpar@345: /// kpeter@615: /// Note that this problem is a special case of the \ref min_cost_flow kpeter@615: /// "minimum cost flow problem". This implementation is actually an kpeter@615: /// efficient specialized version of the \ref CapacityScaling kpeter@615: /// "Successive Shortest Path" algorithm directly for this problem. kpeter@615: /// Therefore this class provides query functions for flow values and kpeter@615: /// node potentials (the dual solution) just like the minimum cost flow kpeter@615: /// algorithms. alpar@345: /// kpeter@550: /// \tparam GR The digraph type the algorithm runs on. kpeter@615: /// \tparam LEN The type of the length map. kpeter@615: /// The default value is GR::ArcMap. alpar@345: /// kpeter@718: /// \warning Length values should be \e non-negative. alpar@345: /// alpar@345: /// \note For finding node-disjoint paths this algorithm can be used kpeter@615: /// along with the \ref SplitNodes adaptor. kpeter@346: #ifdef DOXYGEN kpeter@550: template kpeter@346: #else kpeter@615: template < typename GR, kpeter@550: typename LEN = typename GR::template ArcMap > kpeter@346: #endif alpar@345: class Suurballe alpar@345: { kpeter@550: TEMPLATE_DIGRAPH_TYPEDEFS(GR); alpar@345: alpar@345: typedef ConstMap ConstArcMap; kpeter@550: typedef typename GR::template NodeMap PredMap; alpar@345: alpar@345: public: alpar@345: kpeter@550: /// The type of the digraph the algorithm runs on. kpeter@550: typedef GR Digraph; kpeter@550: /// The type of the length map. kpeter@550: typedef LEN LengthMap; kpeter@550: /// The type of the lengths. kpeter@550: typedef typename LengthMap::Value Length; kpeter@615: #ifdef DOXYGEN kpeter@615: /// The type of the flow map. kpeter@615: typedef GR::ArcMap FlowMap; kpeter@615: /// The type of the potential map. kpeter@615: typedef GR::NodeMap PotentialMap; kpeter@615: #else alpar@345: /// The type of the flow map. alpar@345: typedef typename Digraph::template ArcMap FlowMap; alpar@345: /// The type of the potential map. alpar@345: typedef typename Digraph::template NodeMap PotentialMap; kpeter@615: #endif kpeter@615: alpar@345: /// The type of the path structures. kpeter@615: typedef SimplePath Path; alpar@345: alpar@345: private: alpar@440: kpeter@615: // ResidualDijkstra is a special implementation of the kpeter@615: // Dijkstra algorithm for finding shortest paths in the kpeter@615: // residual network with respect to the reduced arc lengths kpeter@615: // and modifying the node potentials according to the kpeter@615: // distance of the nodes. alpar@345: class ResidualDijkstra alpar@345: { alpar@345: typedef typename Digraph::template NodeMap HeapCrossRef; alpar@345: typedef BinHeap Heap; alpar@345: alpar@345: private: alpar@345: kpeter@346: // The digraph the algorithm runs on alpar@345: const Digraph &_graph; alpar@345: alpar@345: // The main maps alpar@345: const FlowMap &_flow; alpar@345: const LengthMap &_length; alpar@345: PotentialMap &_potential; alpar@345: alpar@345: // The distance map alpar@345: PotentialMap _dist; alpar@345: // The pred arc map alpar@345: PredMap &_pred; alpar@345: // The processed (i.e. permanently labeled) nodes alpar@345: std::vector _proc_nodes; alpar@440: alpar@345: Node _s; alpar@345: Node _t; alpar@345: alpar@345: public: alpar@345: alpar@345: /// Constructor. kpeter@615: ResidualDijkstra( const Digraph &graph, alpar@345: const FlowMap &flow, alpar@345: const LengthMap &length, alpar@345: PotentialMap &potential, alpar@345: PredMap &pred, alpar@345: Node s, Node t ) : kpeter@615: _graph(graph), _flow(flow), _length(length), _potential(potential), kpeter@615: _dist(graph), _pred(pred), _s(s), _t(t) {} alpar@345: kpeter@346: /// \brief Run the algorithm. It returns \c true if a path is found alpar@345: /// from the source node to the target node. alpar@345: bool run() { alpar@345: HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP); alpar@345: Heap heap(heap_cross_ref); alpar@345: heap.push(_s, 0); alpar@345: _pred[_s] = INVALID; alpar@345: _proc_nodes.clear(); alpar@345: kpeter@346: // Process nodes alpar@345: while (!heap.empty() && heap.top() != _t) { alpar@345: Node u = heap.top(), v; alpar@345: Length d = heap.prio() + _potential[u], nd; alpar@345: _dist[u] = heap.prio(); alpar@345: heap.pop(); alpar@345: _proc_nodes.push_back(u); alpar@345: kpeter@346: // Traverse outgoing arcs alpar@345: for (OutArcIt e(_graph, u); e != INVALID; ++e) { alpar@345: if (_flow[e] == 0) { alpar@345: v = _graph.target(e); alpar@345: switch(heap.state(v)) { alpar@345: case Heap::PRE_HEAP: alpar@345: heap.push(v, d + _length[e] - _potential[v]); alpar@345: _pred[v] = e; alpar@345: break; alpar@345: case Heap::IN_HEAP: alpar@345: nd = d + _length[e] - _potential[v]; alpar@345: if (nd < heap[v]) { alpar@345: heap.decrease(v, nd); alpar@345: _pred[v] = e; alpar@345: } alpar@345: break; alpar@345: case Heap::POST_HEAP: alpar@345: break; alpar@345: } alpar@345: } alpar@345: } alpar@345: kpeter@346: // Traverse incoming arcs alpar@345: for (InArcIt e(_graph, u); e != INVALID; ++e) { alpar@345: if (_flow[e] == 1) { alpar@345: v = _graph.source(e); alpar@345: switch(heap.state(v)) { alpar@345: case Heap::PRE_HEAP: alpar@345: heap.push(v, d - _length[e] - _potential[v]); alpar@345: _pred[v] = e; alpar@345: break; alpar@345: case Heap::IN_HEAP: alpar@345: nd = d - _length[e] - _potential[v]; alpar@345: if (nd < heap[v]) { alpar@345: heap.decrease(v, nd); alpar@345: _pred[v] = e; alpar@345: } alpar@345: break; alpar@345: case Heap::POST_HEAP: alpar@345: break; alpar@345: } alpar@345: } alpar@345: } alpar@345: } alpar@345: if (heap.empty()) return false; alpar@345: kpeter@346: // Update potentials of processed nodes alpar@345: Length t_dist = heap.prio(); alpar@345: for (int i = 0; i < int(_proc_nodes.size()); ++i) alpar@345: _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist; alpar@345: return true; alpar@345: } alpar@345: alpar@345: }; //class ResidualDijkstra alpar@345: alpar@345: private: alpar@345: kpeter@346: // The digraph the algorithm runs on alpar@345: const Digraph &_graph; alpar@345: // The length map alpar@345: const LengthMap &_length; alpar@440: alpar@345: // Arc map of the current flow alpar@345: FlowMap *_flow; alpar@345: bool _local_flow; alpar@345: // Node map of the current potentials alpar@345: PotentialMap *_potential; alpar@345: bool _local_potential; alpar@345: alpar@345: // The source node alpar@345: Node _source; alpar@345: // The target node alpar@345: Node _target; alpar@345: alpar@345: // Container to store the found paths alpar@345: std::vector< SimplePath > paths; alpar@345: int _path_num; alpar@345: alpar@345: // The pred arc map alpar@345: PredMap _pred; alpar@345: // Implementation of the Dijkstra algorithm for finding augmenting alpar@345: // shortest paths in the residual network alpar@345: ResidualDijkstra *_dijkstra; alpar@345: alpar@345: public: alpar@345: alpar@345: /// \brief Constructor. alpar@345: /// alpar@345: /// Constructor. alpar@345: /// kpeter@615: /// \param graph The digraph the algorithm runs on. alpar@345: /// \param length The length (cost) values of the arcs. kpeter@615: Suurballe( const Digraph &graph, kpeter@615: const LengthMap &length ) : kpeter@615: _graph(graph), _length(length), _flow(0), _local_flow(false), kpeter@615: _potential(0), _local_potential(false), _pred(graph) kpeter@718: {} alpar@345: alpar@345: /// Destructor. alpar@345: ~Suurballe() { alpar@345: if (_local_flow) delete _flow; alpar@345: if (_local_potential) delete _potential; alpar@345: delete _dijkstra; alpar@345: } alpar@345: kpeter@346: /// \brief Set the flow map. alpar@345: /// kpeter@346: /// This function sets the flow map. kpeter@615: /// If it is not used before calling \ref run() or \ref init(), kpeter@615: /// an instance will be allocated automatically. The destructor kpeter@615: /// deallocates this automatically allocated map, of course. alpar@345: /// kpeter@615: /// The found flow contains only 0 and 1 values, since it is the kpeter@615: /// union of the found arc-disjoint paths. alpar@345: /// kpeter@550: /// \return (*this) alpar@345: Suurballe& flowMap(FlowMap &map) { alpar@345: if (_local_flow) { alpar@345: delete _flow; alpar@345: _local_flow = false; alpar@345: } alpar@345: _flow = ↦ alpar@345: return *this; alpar@345: } alpar@345: kpeter@346: /// \brief Set the potential map. alpar@345: /// kpeter@346: /// This function sets the potential map. kpeter@615: /// If it is not used before calling \ref run() or \ref init(), kpeter@615: /// an instance will be allocated automatically. The destructor kpeter@615: /// deallocates this automatically allocated map, of course. alpar@345: /// kpeter@615: /// The node potentials provide the dual solution of the underlying kpeter@615: /// \ref min_cost_flow "minimum cost flow problem". alpar@345: /// kpeter@550: /// \return (*this) alpar@345: Suurballe& potentialMap(PotentialMap &map) { alpar@345: if (_local_potential) { alpar@345: delete _potential; alpar@345: _local_potential = false; alpar@345: } alpar@345: _potential = ↦ alpar@345: return *this; alpar@345: } alpar@345: kpeter@576: /// \name Execution Control alpar@345: /// The simplest way to execute the algorithm is to call the run() alpar@345: /// function. alpar@345: /// \n alpar@345: /// If you only need the flow that is the union of the found alpar@345: /// arc-disjoint paths, you may call init() and findFlow(). alpar@345: alpar@345: /// @{ alpar@345: kpeter@346: /// \brief Run the algorithm. alpar@345: /// kpeter@346: /// This function runs the algorithm. alpar@345: /// kpeter@615: /// \param s The source node. kpeter@615: /// \param t The target node. alpar@345: /// \param k The number of paths to be found. alpar@345: /// kpeter@346: /// \return \c k if there are at least \c k arc-disjoint paths from kpeter@346: /// \c s to \c t in the digraph. Otherwise it returns the number of alpar@345: /// arc-disjoint paths found. alpar@345: /// kpeter@615: /// \note Apart from the return value, s.run(s, t, k) is kpeter@615: /// just a shortcut of the following code. alpar@345: /// \code kpeter@615: /// s.init(s); kpeter@615: /// s.findFlow(t, k); alpar@345: /// s.findPaths(); alpar@345: /// \endcode kpeter@615: int run(const Node& s, const Node& t, int k = 2) { kpeter@615: init(s); kpeter@615: findFlow(t, k); alpar@345: findPaths(); alpar@345: return _path_num; alpar@345: } alpar@345: kpeter@346: /// \brief Initialize the algorithm. alpar@345: /// kpeter@346: /// This function initializes the algorithm. kpeter@615: /// kpeter@615: /// \param s The source node. kpeter@615: void init(const Node& s) { kpeter@615: _source = s; kpeter@615: kpeter@346: // Initialize maps alpar@345: if (!_flow) { alpar@345: _flow = new FlowMap(_graph); alpar@345: _local_flow = true; alpar@345: } alpar@345: if (!_potential) { alpar@345: _potential = new PotentialMap(_graph); alpar@345: _local_potential = true; alpar@345: } alpar@345: for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0; alpar@345: for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0; alpar@345: } alpar@345: kpeter@615: /// \brief Execute the algorithm to find an optimal flow. alpar@345: /// kpeter@346: /// This function executes the successive shortest path algorithm to kpeter@615: /// find a minimum cost flow, which is the union of \c k (or less) alpar@345: /// arc-disjoint paths. alpar@345: /// kpeter@615: /// \param t The target node. kpeter@615: /// \param k The number of paths to be found. kpeter@615: /// kpeter@346: /// \return \c k if there are at least \c k arc-disjoint paths from kpeter@615: /// the source node to the given node \c t in the digraph. kpeter@615: /// Otherwise it returns the number of arc-disjoint paths found. alpar@345: /// alpar@345: /// \pre \ref init() must be called before using this function. kpeter@615: int findFlow(const Node& t, int k = 2) { kpeter@615: _target = t; kpeter@615: _dijkstra = kpeter@615: new ResidualDijkstra( _graph, *_flow, _length, *_potential, _pred, kpeter@615: _source, _target ); kpeter@615: kpeter@346: // Find shortest paths alpar@345: _path_num = 0; alpar@345: while (_path_num < k) { kpeter@346: // Run Dijkstra alpar@345: if (!_dijkstra->run()) break; alpar@345: ++_path_num; alpar@345: kpeter@346: // Set the flow along the found shortest path alpar@345: Node u = _target; alpar@345: Arc e; alpar@345: while ((e = _pred[u]) != INVALID) { alpar@345: if (u == _graph.target(e)) { alpar@345: (*_flow)[e] = 1; alpar@345: u = _graph.source(e); alpar@345: } else { alpar@345: (*_flow)[e] = 0; alpar@345: u = _graph.target(e); alpar@345: } alpar@345: } alpar@345: } alpar@345: return _path_num; alpar@345: } alpar@440: kpeter@346: /// \brief Compute the paths from the flow. alpar@345: /// kpeter@615: /// This function computes the paths from the found minimum cost flow, kpeter@615: /// which is the union of some arc-disjoint paths. alpar@345: /// alpar@345: /// \pre \ref init() and \ref findFlow() must be called before using alpar@345: /// this function. alpar@345: void findPaths() { alpar@345: FlowMap res_flow(_graph); kpeter@346: for(ArcIt a(_graph); a != INVALID; ++a) res_flow[a] = (*_flow)[a]; alpar@345: alpar@345: paths.clear(); alpar@345: paths.resize(_path_num); alpar@345: for (int i = 0; i < _path_num; ++i) { alpar@345: Node n = _source; alpar@345: while (n != _target) { alpar@345: OutArcIt e(_graph, n); alpar@345: for ( ; res_flow[e] == 0; ++e) ; alpar@345: n = _graph.target(e); alpar@345: paths[i].addBack(e); alpar@345: res_flow[e] = 0; alpar@345: } alpar@345: } alpar@345: } alpar@345: alpar@345: /// @} alpar@345: alpar@345: /// \name Query Functions kpeter@346: /// The results of the algorithm can be obtained using these alpar@345: /// functions. alpar@345: /// \n The algorithm should be executed before using them. alpar@345: alpar@345: /// @{ alpar@345: kpeter@615: /// \brief Return the total length of the found paths. kpeter@615: /// kpeter@615: /// This function returns the total length of the found paths, i.e. kpeter@615: /// the total cost of the found flow. kpeter@615: /// The complexity of the function is O(e). kpeter@615: /// kpeter@615: /// \pre \ref run() or \ref findFlow() must be called before using kpeter@615: /// this function. kpeter@615: Length totalLength() const { kpeter@615: Length c = 0; kpeter@615: for (ArcIt e(_graph); e != INVALID; ++e) kpeter@615: c += (*_flow)[e] * _length[e]; kpeter@615: return c; kpeter@615: } kpeter@615: kpeter@615: /// \brief Return the flow value on the given arc. kpeter@615: /// kpeter@615: /// This function returns the flow value on the given arc. kpeter@615: /// It is \c 1 if the arc is involved in one of the found arc-disjoint kpeter@615: /// paths, otherwise it is \c 0. kpeter@615: /// kpeter@615: /// \pre \ref run() or \ref findFlow() must be called before using kpeter@615: /// this function. kpeter@615: int flow(const Arc& arc) const { kpeter@615: return (*_flow)[arc]; kpeter@615: } kpeter@615: kpeter@615: /// \brief Return a const reference to an arc map storing the alpar@345: /// found flow. alpar@345: /// kpeter@615: /// This function returns a const reference to an arc map storing kpeter@346: /// the flow that is the union of the found arc-disjoint paths. alpar@345: /// kpeter@346: /// \pre \ref run() or \ref findFlow() must be called before using kpeter@346: /// this function. alpar@345: const FlowMap& flowMap() const { alpar@345: return *_flow; alpar@345: } alpar@345: kpeter@346: /// \brief Return the potential of the given node. alpar@345: /// kpeter@346: /// This function returns the potential of the given node. kpeter@615: /// The node potentials provide the dual solution of the kpeter@615: /// underlying \ref min_cost_flow "minimum cost flow problem". alpar@345: /// kpeter@346: /// \pre \ref run() or \ref findFlow() must be called before using kpeter@346: /// this function. alpar@345: Length potential(const Node& node) const { alpar@345: return (*_potential)[node]; alpar@345: } alpar@345: kpeter@615: /// \brief Return a const reference to a node map storing the kpeter@615: /// found potentials (the dual solution). alpar@345: /// kpeter@615: /// This function returns a const reference to a node map storing kpeter@615: /// the found potentials that provide the dual solution of the kpeter@615: /// underlying \ref min_cost_flow "minimum cost flow problem". alpar@345: /// kpeter@346: /// \pre \ref run() or \ref findFlow() must be called before using kpeter@346: /// this function. kpeter@615: const PotentialMap& potentialMap() const { kpeter@615: return *_potential; alpar@345: } alpar@345: kpeter@346: /// \brief Return the number of the found paths. alpar@345: /// kpeter@346: /// This function returns the number of the found paths. alpar@345: /// kpeter@346: /// \pre \ref run() or \ref findFlow() must be called before using kpeter@346: /// this function. alpar@345: int pathNum() const { alpar@345: return _path_num; alpar@345: } alpar@345: kpeter@346: /// \brief Return a const reference to the specified path. alpar@345: /// kpeter@346: /// This function returns a const reference to the specified path. alpar@345: /// kpeter@615: /// \param i The function returns the i-th path. alpar@345: /// \c i must be between \c 0 and %pathNum()-1. alpar@345: /// kpeter@346: /// \pre \ref run() or \ref findPaths() must be called before using kpeter@346: /// this function. kpeter@717: const Path& path(int i) const { alpar@345: return paths[i]; alpar@345: } alpar@345: alpar@345: /// @} alpar@345: alpar@345: }; //class Suurballe alpar@345: alpar@345: ///@} alpar@345: alpar@345: } //namespace lemon alpar@345: alpar@345: #endif //LEMON_SUURBALLE_H