Rework and improve Suurballe (#323)
authorPeter Kovacs <kpeter@inf.elte.hu>
Fri, 16 Oct 2009 01:06:16 +0200
changeset 853ec0b1b423b8b
parent 852 30c77d1c0cba
child 854 9a7e4e606f83
Rework and improve Suurballe (#323)

- Improve the implementation: use a specific, faster variant of
residual Dijkstra for the first search.
- Some reorganizatiopn to make the code simpler.
- Small doc improvements.
lemon/suurballe.h
     1.1 --- a/lemon/suurballe.h	Thu Oct 15 21:04:50 2009 +0200
     1.2 +++ b/lemon/suurballe.h	Fri Oct 16 01:06:16 2009 +0200
     1.3 @@ -46,7 +46,7 @@
     1.4    /// Note that this problem is a special case of the \ref min_cost_flow
     1.5    /// "minimum cost flow problem". This implementation is actually an
     1.6    /// efficient specialized version of the \ref CapacityScaling
     1.7 -  /// "Successive Shortest Path" algorithm directly for this problem.
     1.8 +  /// "successive shortest path" algorithm directly for this problem.
     1.9    /// Therefore this class provides query functions for flow values and
    1.10    /// node potentials (the dual solution) just like the minimum cost flow
    1.11    /// algorithms.
    1.12 @@ -57,7 +57,7 @@
    1.13    ///
    1.14    /// \warning Length values should be \e non-negative.
    1.15    ///
    1.16 -  /// \note For finding node-disjoint paths this algorithm can be used
    1.17 +  /// \note For finding \e node-disjoint paths, this algorithm can be used
    1.18    /// along with the \ref SplitNodes adaptor.
    1.19  #ifdef DOXYGEN
    1.20    template <typename GR, typename LEN>
    1.21 @@ -109,39 +109,36 @@
    1.22  
    1.23      private:
    1.24  
    1.25 -      // The digraph the algorithm runs on
    1.26        const Digraph &_graph;
    1.27 -
    1.28 -      // The main maps
    1.29 +      const LengthMap &_length;
    1.30        const FlowMap &_flow;
    1.31 -      const LengthMap &_length;
    1.32 -      PotentialMap &_potential;
    1.33 -
    1.34 -      // The distance map
    1.35 -      PotentialMap _dist;
    1.36 -      // The pred arc map
    1.37 +      PotentialMap &_pi;
    1.38        PredMap &_pred;
    1.39 -      // The processed (i.e. permanently labeled) nodes
    1.40 -      std::vector<Node> _proc_nodes;
    1.41 -
    1.42        Node _s;
    1.43        Node _t;
    1.44 +      
    1.45 +      PotentialMap _dist;
    1.46 +      std::vector<Node> _proc_nodes;
    1.47  
    1.48      public:
    1.49  
    1.50 -      /// Constructor.
    1.51 -      ResidualDijkstra( const Digraph &graph,
    1.52 -                        const FlowMap &flow,
    1.53 -                        const LengthMap &length,
    1.54 -                        PotentialMap &potential,
    1.55 -                        PredMap &pred,
    1.56 -                        Node s, Node t ) :
    1.57 -        _graph(graph), _flow(flow), _length(length), _potential(potential),
    1.58 -        _dist(graph), _pred(pred), _s(s), _t(t) {}
    1.59 +      // Constructor
    1.60 +      ResidualDijkstra(Suurballe &srb) :
    1.61 +        _graph(srb._graph), _length(srb._length),
    1.62 +        _flow(*srb._flow), _pi(*srb._potential), _pred(srb._pred), 
    1.63 +        _s(srb._s), _t(srb._t), _dist(_graph) {}
    1.64 +        
    1.65 +      // Run the algorithm and return true if a path is found
    1.66 +      // from the source node to the target node.
    1.67 +      bool run(int cnt) {
    1.68 +        return cnt == 0 ? startFirst() : start();
    1.69 +      }
    1.70  
    1.71 -      /// \brief Run the algorithm. It returns \c true if a path is found
    1.72 -      /// from the source node to the target node.
    1.73 -      bool run() {
    1.74 +    private:
    1.75 +    
    1.76 +      // Execute the algorithm for the first time (the flow and potential
    1.77 +      // functions have to be identically zero).
    1.78 +      bool startFirst() {
    1.79          HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
    1.80          Heap heap(heap_cross_ref);
    1.81          heap.push(_s, 0);
    1.82 @@ -151,29 +148,74 @@
    1.83          // Process nodes
    1.84          while (!heap.empty() && heap.top() != _t) {
    1.85            Node u = heap.top(), v;
    1.86 -          Length d = heap.prio() + _potential[u], nd;
    1.87 +          Length d = heap.prio(), dn;
    1.88            _dist[u] = heap.prio();
    1.89 +          _proc_nodes.push_back(u);
    1.90            heap.pop();
    1.91 +
    1.92 +          // Traverse outgoing arcs
    1.93 +          for (OutArcIt e(_graph, u); e != INVALID; ++e) {
    1.94 +            v = _graph.target(e);
    1.95 +            switch(heap.state(v)) {
    1.96 +              case Heap::PRE_HEAP:
    1.97 +                heap.push(v, d + _length[e]);
    1.98 +                _pred[v] = e;
    1.99 +                break;
   1.100 +              case Heap::IN_HEAP:
   1.101 +                dn = d + _length[e];
   1.102 +                if (dn < heap[v]) {
   1.103 +                  heap.decrease(v, dn);
   1.104 +                  _pred[v] = e;
   1.105 +                }
   1.106 +                break;
   1.107 +              case Heap::POST_HEAP:
   1.108 +                break;
   1.109 +            }
   1.110 +          }
   1.111 +        }
   1.112 +        if (heap.empty()) return false;
   1.113 +
   1.114 +        // Update potentials of processed nodes
   1.115 +        Length t_dist = heap.prio();
   1.116 +        for (int i = 0; i < int(_proc_nodes.size()); ++i)
   1.117 +          _pi[_proc_nodes[i]] = _dist[_proc_nodes[i]] - t_dist;
   1.118 +        return true;
   1.119 +      }
   1.120 +
   1.121 +      // Execute the algorithm.
   1.122 +      bool start() {
   1.123 +        HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
   1.124 +        Heap heap(heap_cross_ref);
   1.125 +        heap.push(_s, 0);
   1.126 +        _pred[_s] = INVALID;
   1.127 +        _proc_nodes.clear();
   1.128 +
   1.129 +        // Process nodes
   1.130 +        while (!heap.empty() && heap.top() != _t) {
   1.131 +          Node u = heap.top(), v;
   1.132 +          Length d = heap.prio() + _pi[u], dn;
   1.133 +          _dist[u] = heap.prio();
   1.134            _proc_nodes.push_back(u);
   1.135 +          heap.pop();
   1.136  
   1.137            // Traverse outgoing arcs
   1.138            for (OutArcIt e(_graph, u); e != INVALID; ++e) {
   1.139              if (_flow[e] == 0) {
   1.140                v = _graph.target(e);
   1.141                switch(heap.state(v)) {
   1.142 -              case Heap::PRE_HEAP:
   1.143 -                heap.push(v, d + _length[e] - _potential[v]);
   1.144 -                _pred[v] = e;
   1.145 -                break;
   1.146 -              case Heap::IN_HEAP:
   1.147 -                nd = d + _length[e] - _potential[v];
   1.148 -                if (nd < heap[v]) {
   1.149 -                  heap.decrease(v, nd);
   1.150 +                case Heap::PRE_HEAP:
   1.151 +                  heap.push(v, d + _length[e] - _pi[v]);
   1.152                    _pred[v] = e;
   1.153 -                }
   1.154 -                break;
   1.155 -              case Heap::POST_HEAP:
   1.156 -                break;
   1.157 +                  break;
   1.158 +                case Heap::IN_HEAP:
   1.159 +                  dn = d + _length[e] - _pi[v];
   1.160 +                  if (dn < heap[v]) {
   1.161 +                    heap.decrease(v, dn);
   1.162 +                    _pred[v] = e;
   1.163 +                  }
   1.164 +                  break;
   1.165 +                case Heap::POST_HEAP:
   1.166 +                  break;
   1.167                }
   1.168              }
   1.169            }
   1.170 @@ -183,19 +225,19 @@
   1.171              if (_flow[e] == 1) {
   1.172                v = _graph.source(e);
   1.173                switch(heap.state(v)) {
   1.174 -              case Heap::PRE_HEAP:
   1.175 -                heap.push(v, d - _length[e] - _potential[v]);
   1.176 -                _pred[v] = e;
   1.177 -                break;
   1.178 -              case Heap::IN_HEAP:
   1.179 -                nd = d - _length[e] - _potential[v];
   1.180 -                if (nd < heap[v]) {
   1.181 -                  heap.decrease(v, nd);
   1.182 +                case Heap::PRE_HEAP:
   1.183 +                  heap.push(v, d - _length[e] - _pi[v]);
   1.184                    _pred[v] = e;
   1.185 -                }
   1.186 -                break;
   1.187 -              case Heap::POST_HEAP:
   1.188 -                break;
   1.189 +                  break;
   1.190 +                case Heap::IN_HEAP:
   1.191 +                  dn = d - _length[e] - _pi[v];
   1.192 +                  if (dn < heap[v]) {
   1.193 +                    heap.decrease(v, dn);
   1.194 +                    _pred[v] = e;
   1.195 +                  }
   1.196 +                  break;
   1.197 +                case Heap::POST_HEAP:
   1.198 +                  break;
   1.199                }
   1.200              }
   1.201            }
   1.202 @@ -205,7 +247,7 @@
   1.203          // Update potentials of processed nodes
   1.204          Length t_dist = heap.prio();
   1.205          for (int i = 0; i < int(_proc_nodes.size()); ++i)
   1.206 -          _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
   1.207 +          _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
   1.208          return true;
   1.209        }
   1.210  
   1.211 @@ -226,19 +268,16 @@
   1.212      bool _local_potential;
   1.213  
   1.214      // The source node
   1.215 -    Node _source;
   1.216 +    Node _s;
   1.217      // The target node
   1.218 -    Node _target;
   1.219 +    Node _t;
   1.220  
   1.221      // Container to store the found paths
   1.222 -    std::vector< SimplePath<Digraph> > paths;
   1.223 +    std::vector<Path> _paths;
   1.224      int _path_num;
   1.225  
   1.226      // The pred arc map
   1.227      PredMap _pred;
   1.228 -    // Implementation of the Dijkstra algorithm for finding augmenting
   1.229 -    // shortest paths in the residual network
   1.230 -    ResidualDijkstra *_dijkstra;
   1.231  
   1.232    public:
   1.233  
   1.234 @@ -258,7 +297,6 @@
   1.235      ~Suurballe() {
   1.236        if (_local_flow) delete _flow;
   1.237        if (_local_potential) delete _potential;
   1.238 -      delete _dijkstra;
   1.239      }
   1.240  
   1.241      /// \brief Set the flow map.
   1.242 @@ -342,7 +380,7 @@
   1.243      ///
   1.244      /// \param s The source node.
   1.245      void init(const Node& s) {
   1.246 -      _source = s;
   1.247 +      _s = s;
   1.248  
   1.249        // Initialize maps
   1.250        if (!_flow) {
   1.251 @@ -372,20 +410,18 @@
   1.252      ///
   1.253      /// \pre \ref init() must be called before using this function.
   1.254      int findFlow(const Node& t, int k = 2) {
   1.255 -      _target = t;
   1.256 -      _dijkstra =
   1.257 -        new ResidualDijkstra( _graph, *_flow, _length, *_potential, _pred,
   1.258 -                              _source, _target );
   1.259 +      _t = t;
   1.260 +      ResidualDijkstra dijkstra(*this);
   1.261  
   1.262        // Find shortest paths
   1.263        _path_num = 0;
   1.264        while (_path_num < k) {
   1.265          // Run Dijkstra
   1.266 -        if (!_dijkstra->run()) break;
   1.267 +        if (!dijkstra.run(_path_num)) break;
   1.268          ++_path_num;
   1.269  
   1.270          // Set the flow along the found shortest path
   1.271 -        Node u = _target;
   1.272 +        Node u = _t;
   1.273          Arc e;
   1.274          while ((e = _pred[u]) != INVALID) {
   1.275            if (u == _graph.target(e)) {
   1.276 @@ -402,8 +438,8 @@
   1.277  
   1.278      /// \brief Compute the paths from the flow.
   1.279      ///
   1.280 -    /// This function computes the paths from the found minimum cost flow,
   1.281 -    /// which is the union of some arc-disjoint paths.
   1.282 +    /// This function computes arc-disjoint paths from the found minimum
   1.283 +    /// cost flow, which is the union of them.
   1.284      ///
   1.285      /// \pre \ref init() and \ref findFlow() must be called before using
   1.286      /// this function.
   1.287 @@ -411,15 +447,15 @@
   1.288        FlowMap res_flow(_graph);
   1.289        for(ArcIt a(_graph); a != INVALID; ++a) res_flow[a] = (*_flow)[a];
   1.290  
   1.291 -      paths.clear();
   1.292 -      paths.resize(_path_num);
   1.293 +      _paths.clear();
   1.294 +      _paths.resize(_path_num);
   1.295        for (int i = 0; i < _path_num; ++i) {
   1.296 -        Node n = _source;
   1.297 -        while (n != _target) {
   1.298 +        Node n = _s;
   1.299 +        while (n != _t) {
   1.300            OutArcIt e(_graph, n);
   1.301            for ( ; res_flow[e] == 0; ++e) ;
   1.302            n = _graph.target(e);
   1.303 -          paths[i].addBack(e);
   1.304 +          _paths[i].addBack(e);
   1.305            res_flow[e] = 0;
   1.306          }
   1.307        }
   1.308 @@ -518,7 +554,7 @@
   1.309      /// \pre \ref run() or \ref findPaths() must be called before using
   1.310      /// this function.
   1.311      const Path& path(int i) const {
   1.312 -      return paths[i];
   1.313 +      return _paths[i];
   1.314      }
   1.315  
   1.316      /// @}