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@623: #include alpar@345: #include alpar@345: #include deba@519: #include deba@519: #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@623: /// Note that this problem is a special case of the \ref min_cost_flow kpeter@623: /// "minimum cost flow problem". This implementation is actually an kpeter@623: /// efficient specialized version of the \ref CapacityScaling kpeter@853: /// "successive shortest path" algorithm directly for this problem. kpeter@623: /// Therefore this class provides query functions for flow values and kpeter@623: /// node potentials (the dual solution) just like the minimum cost flow kpeter@623: /// algorithms. alpar@345: /// kpeter@559: /// \tparam GR The digraph type the algorithm runs on. kpeter@623: /// \tparam LEN The type of the length map. kpeter@623: /// The default value is GR::ArcMap. alpar@345: /// kpeter@852: /// \warning Length values should be \e non-negative. alpar@345: /// kpeter@853: /// \note For finding \e node-disjoint paths, this algorithm can be used kpeter@623: /// along with the \ref SplitNodes adaptor. kpeter@346: #ifdef DOXYGEN kpeter@559: template kpeter@346: #else kpeter@623: template < typename GR, kpeter@559: typename LEN = typename GR::template ArcMap > kpeter@346: #endif alpar@345: class Suurballe alpar@345: { kpeter@559: TEMPLATE_DIGRAPH_TYPEDEFS(GR); alpar@345: alpar@345: typedef ConstMap ConstArcMap; kpeter@559: typedef typename GR::template NodeMap PredMap; alpar@345: alpar@345: public: alpar@345: kpeter@559: /// The type of the digraph the algorithm runs on. kpeter@559: typedef GR Digraph; kpeter@559: /// The type of the length map. kpeter@559: typedef LEN LengthMap; kpeter@559: /// The type of the lengths. kpeter@559: typedef typename LengthMap::Value Length; kpeter@623: #ifdef DOXYGEN kpeter@623: /// The type of the flow map. kpeter@623: typedef GR::ArcMap FlowMap; kpeter@623: /// The type of the potential map. kpeter@623: typedef GR::NodeMap PotentialMap; kpeter@623: #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@623: #endif kpeter@623: alpar@345: /// The type of the path structures. kpeter@623: typedef SimplePath Path; alpar@345: alpar@345: private: alpar@440: kpeter@623: // ResidualDijkstra is a special implementation of the kpeter@623: // Dijkstra algorithm for finding shortest paths in the kpeter@623: // residual network with respect to the reduced arc lengths kpeter@623: // and modifying the node potentials according to the kpeter@623: // 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: alpar@345: const Digraph &_graph; kpeter@853: const LengthMap &_length; alpar@345: const FlowMap &_flow; kpeter@853: PotentialMap &_pi; alpar@345: PredMap &_pred; alpar@345: Node _s; alpar@345: Node _t; kpeter@853: kpeter@853: PotentialMap _dist; kpeter@853: std::vector _proc_nodes; alpar@345: alpar@345: public: alpar@345: kpeter@853: // Constructor kpeter@853: ResidualDijkstra(Suurballe &srb) : kpeter@853: _graph(srb._graph), _length(srb._length), kpeter@853: _flow(*srb._flow), _pi(*srb._potential), _pred(srb._pred), kpeter@853: _s(srb._s), _t(srb._t), _dist(_graph) {} kpeter@853: kpeter@853: // Run the algorithm and return true if a path is found kpeter@853: // from the source node to the target node. kpeter@853: bool run(int cnt) { kpeter@853: return cnt == 0 ? startFirst() : start(); kpeter@853: } alpar@345: kpeter@853: private: kpeter@853: kpeter@853: // Execute the algorithm for the first time (the flow and potential kpeter@853: // functions have to be identically zero). kpeter@853: bool startFirst() { 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; kpeter@853: Length d = heap.prio(), dn; alpar@345: _dist[u] = heap.prio(); kpeter@853: _proc_nodes.push_back(u); alpar@345: heap.pop(); kpeter@853: kpeter@853: // Traverse outgoing arcs kpeter@853: for (OutArcIt e(_graph, u); e != INVALID; ++e) { kpeter@853: v = _graph.target(e); kpeter@853: switch(heap.state(v)) { kpeter@853: case Heap::PRE_HEAP: kpeter@853: heap.push(v, d + _length[e]); kpeter@853: _pred[v] = e; kpeter@853: break; kpeter@853: case Heap::IN_HEAP: kpeter@853: dn = d + _length[e]; kpeter@853: if (dn < heap[v]) { kpeter@853: heap.decrease(v, dn); kpeter@853: _pred[v] = e; kpeter@853: } kpeter@853: break; kpeter@853: case Heap::POST_HEAP: kpeter@853: break; kpeter@853: } kpeter@853: } kpeter@853: } kpeter@853: if (heap.empty()) return false; kpeter@853: kpeter@853: // Update potentials of processed nodes kpeter@853: Length t_dist = heap.prio(); kpeter@853: for (int i = 0; i < int(_proc_nodes.size()); ++i) kpeter@853: _pi[_proc_nodes[i]] = _dist[_proc_nodes[i]] - t_dist; kpeter@853: return true; kpeter@853: } kpeter@853: kpeter@853: // Execute the algorithm. kpeter@853: bool start() { kpeter@853: HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP); kpeter@853: Heap heap(heap_cross_ref); kpeter@853: heap.push(_s, 0); kpeter@853: _pred[_s] = INVALID; kpeter@853: _proc_nodes.clear(); kpeter@853: kpeter@853: // Process nodes kpeter@853: while (!heap.empty() && heap.top() != _t) { kpeter@853: Node u = heap.top(), v; kpeter@853: Length d = heap.prio() + _pi[u], dn; kpeter@853: _dist[u] = heap.prio(); alpar@345: _proc_nodes.push_back(u); kpeter@853: heap.pop(); 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)) { kpeter@853: case Heap::PRE_HEAP: kpeter@853: heap.push(v, d + _length[e] - _pi[v]); alpar@345: _pred[v] = e; kpeter@853: break; kpeter@853: case Heap::IN_HEAP: kpeter@853: dn = d + _length[e] - _pi[v]; kpeter@853: if (dn < heap[v]) { kpeter@853: heap.decrease(v, dn); kpeter@853: _pred[v] = e; kpeter@853: } kpeter@853: break; kpeter@853: case Heap::POST_HEAP: kpeter@853: 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)) { kpeter@853: case Heap::PRE_HEAP: kpeter@853: heap.push(v, d - _length[e] - _pi[v]); alpar@345: _pred[v] = e; kpeter@853: break; kpeter@853: case Heap::IN_HEAP: kpeter@853: dn = d - _length[e] - _pi[v]; kpeter@853: if (dn < heap[v]) { kpeter@853: heap.decrease(v, dn); kpeter@853: _pred[v] = e; kpeter@853: } kpeter@853: break; kpeter@853: case Heap::POST_HEAP: kpeter@853: 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) kpeter@853: _pi[_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 kpeter@853: Node _s; alpar@345: // The target node kpeter@853: Node _t; alpar@345: alpar@345: // Container to store the found paths kpeter@853: std::vector _paths; alpar@345: int _path_num; alpar@345: alpar@345: // The pred arc map alpar@345: PredMap _pred; alpar@345: alpar@345: public: alpar@345: alpar@345: /// \brief Constructor. alpar@345: /// alpar@345: /// Constructor. alpar@345: /// kpeter@623: /// \param graph The digraph the algorithm runs on. alpar@345: /// \param length The length (cost) values of the arcs. kpeter@623: Suurballe( const Digraph &graph, kpeter@623: const LengthMap &length ) : kpeter@623: _graph(graph), _length(length), _flow(0), _local_flow(false), kpeter@623: _potential(0), _local_potential(false), _pred(graph) kpeter@852: {} alpar@345: alpar@345: /// Destructor. alpar@345: ~Suurballe() { alpar@345: if (_local_flow) delete _flow; alpar@345: if (_local_potential) delete _potential; alpar@345: } alpar@345: kpeter@346: /// \brief Set the flow map. alpar@345: /// kpeter@346: /// This function sets the flow map. kpeter@623: /// If it is not used before calling \ref run() or \ref init(), kpeter@623: /// an instance will be allocated automatically. The destructor kpeter@623: /// deallocates this automatically allocated map, of course. alpar@345: /// kpeter@623: /// The found flow contains only 0 and 1 values, since it is the kpeter@623: /// union of the found arc-disjoint paths. alpar@345: /// kpeter@559: /// \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@623: /// If it is not used before calling \ref run() or \ref init(), kpeter@623: /// an instance will be allocated automatically. The destructor kpeter@623: /// deallocates this automatically allocated map, of course. alpar@345: /// kpeter@623: /// The node potentials provide the dual solution of the underlying kpeter@623: /// \ref min_cost_flow "minimum cost flow problem". alpar@345: /// kpeter@559: /// \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@584: /// \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@623: /// \param s The source node. kpeter@623: /// \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@623: /// \note Apart from the return value, s.run(s, t, k) is kpeter@623: /// just a shortcut of the following code. alpar@345: /// \code kpeter@623: /// s.init(s); kpeter@623: /// s.findFlow(t, k); alpar@345: /// s.findPaths(); alpar@345: /// \endcode kpeter@623: int run(const Node& s, const Node& t, int k = 2) { kpeter@623: init(s); kpeter@623: 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@623: /// kpeter@623: /// \param s The source node. kpeter@623: void init(const Node& s) { kpeter@853: _s = s; kpeter@623: 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@623: /// \brief Execute the algorithm to find an optimal flow. alpar@345: /// kpeter@346: /// This function executes the successive shortest path algorithm to kpeter@623: /// find a minimum cost flow, which is the union of \c k (or less) alpar@345: /// arc-disjoint paths. alpar@345: /// kpeter@623: /// \param t The target node. kpeter@623: /// \param k The number of paths to be found. kpeter@623: /// kpeter@346: /// \return \c k if there are at least \c k arc-disjoint paths from kpeter@623: /// the source node to the given node \c t in the digraph. kpeter@623: /// 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@623: int findFlow(const Node& t, int k = 2) { kpeter@853: _t = t; kpeter@853: ResidualDijkstra dijkstra(*this); kpeter@623: kpeter@346: // Find shortest paths alpar@345: _path_num = 0; alpar@345: while (_path_num < k) { kpeter@346: // Run Dijkstra kpeter@853: if (!dijkstra.run(_path_num)) break; alpar@345: ++_path_num; alpar@345: kpeter@346: // Set the flow along the found shortest path kpeter@853: Node u = _t; 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@853: /// This function computes arc-disjoint paths from the found minimum kpeter@853: /// cost flow, which is the union of them. 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: kpeter@853: _paths.clear(); kpeter@853: _paths.resize(_path_num); alpar@345: for (int i = 0; i < _path_num; ++i) { kpeter@853: Node n = _s; kpeter@853: while (n != _t) { alpar@345: OutArcIt e(_graph, n); alpar@345: for ( ; res_flow[e] == 0; ++e) ; alpar@345: n = _graph.target(e); kpeter@853: _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@623: /// \brief Return the total length of the found paths. kpeter@623: /// kpeter@623: /// This function returns the total length of the found paths, i.e. kpeter@623: /// the total cost of the found flow. kpeter@623: /// The complexity of the function is O(e). kpeter@623: /// kpeter@623: /// \pre \ref run() or \ref findFlow() must be called before using kpeter@623: /// this function. kpeter@623: Length totalLength() const { kpeter@623: Length c = 0; kpeter@623: for (ArcIt e(_graph); e != INVALID; ++e) kpeter@623: c += (*_flow)[e] * _length[e]; kpeter@623: return c; kpeter@623: } kpeter@623: kpeter@623: /// \brief Return the flow value on the given arc. kpeter@623: /// kpeter@623: /// This function returns the flow value on the given arc. kpeter@623: /// It is \c 1 if the arc is involved in one of the found arc-disjoint kpeter@623: /// paths, otherwise it is \c 0. kpeter@623: /// kpeter@623: /// \pre \ref run() or \ref findFlow() must be called before using kpeter@623: /// this function. kpeter@623: int flow(const Arc& arc) const { kpeter@623: return (*_flow)[arc]; kpeter@623: } kpeter@623: kpeter@623: /// \brief Return a const reference to an arc map storing the alpar@345: /// found flow. alpar@345: /// kpeter@623: /// 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@623: /// The node potentials provide the dual solution of the kpeter@623: /// 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@623: /// \brief Return a const reference to a node map storing the kpeter@623: /// found potentials (the dual solution). alpar@345: /// kpeter@623: /// This function returns a const reference to a node map storing kpeter@623: /// the found potentials that provide the dual solution of the kpeter@623: /// 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@623: const PotentialMap& potentialMap() const { kpeter@623: 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@623: /// \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@851: const Path& path(int i) const { kpeter@853: 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