COIN-OR::LEMON - Graph Library

Ticket #323: 323-rework-ec0b1b423b8b.patch

File 323-rework-ec0b1b423b8b.patch, 10.4 KB (added by Peter Kovacs, 15 years ago)
  • lemon/suurballe.h

    # HG changeset patch
    # User Peter Kovacs <kpeter@inf.elte.hu>
    # Date 1255647976 -7200
    # Node ID ec0b1b423b8b100883ab1ef766204eb978363c80
    # Parent  30c77d1c0cba549cb3f5309642d576cfa9c16534
    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.
    
    diff --git a/lemon/suurballe.h b/lemon/suurballe.h
    a b  
    4646  /// Note that this problem is a special case of the \ref min_cost_flow
    4747  /// "minimum cost flow problem". This implementation is actually an
    4848  /// efficient specialized version of the \ref CapacityScaling
    49   /// "Successive Shortest Path" algorithm directly for this problem.
     49  /// "successive shortest path" algorithm directly for this problem.
    5050  /// Therefore this class provides query functions for flow values and
    5151  /// node potentials (the dual solution) just like the minimum cost flow
    5252  /// algorithms.
     
    5757  ///
    5858  /// \warning Length values should be \e non-negative.
    5959  ///
    60   /// \note For finding node-disjoint paths this algorithm can be used
     60  /// \note For finding \e node-disjoint paths, this algorithm can be used
    6161  /// along with the \ref SplitNodes adaptor.
    6262#ifdef DOXYGEN
    6363  template <typename GR, typename LEN>
     
    109109
    110110    private:
    111111
    112       // The digraph the algorithm runs on
    113112      const Digraph &_graph;
    114 
    115       // The main maps
     113      const LengthMap &_length;
    116114      const FlowMap &_flow;
    117       const LengthMap &_length;
    118       PotentialMap &_potential;
    119 
    120       // The distance map
    121       PotentialMap _dist;
    122       // The pred arc map
     115      PotentialMap &_pi;
    123116      PredMap &_pred;
    124       // The processed (i.e. permanently labeled) nodes
    125       std::vector<Node> _proc_nodes;
    126 
    127117      Node _s;
    128118      Node _t;
     119     
     120      PotentialMap _dist;
     121      std::vector<Node> _proc_nodes;
    129122
    130123    public:
    131124
    132       /// Constructor.
    133       ResidualDijkstra( const Digraph &graph,
    134                         const FlowMap &flow,
    135                         const LengthMap &length,
    136                         PotentialMap &potential,
    137                         PredMap &pred,
    138                         Node s, Node t ) :
    139         _graph(graph), _flow(flow), _length(length), _potential(potential),
    140         _dist(graph), _pred(pred), _s(s), _t(t) {}
     125      // Constructor
     126      ResidualDijkstra(Suurballe &srb) :
     127        _graph(srb._graph), _length(srb._length),
     128        _flow(*srb._flow), _pi(*srb._potential), _pred(srb._pred),
     129        _s(srb._s), _t(srb._t), _dist(_graph) {}
     130       
     131      // Run the algorithm and return true if a path is found
     132      // from the source node to the target node.
     133      bool run(int cnt) {
     134        return cnt == 0 ? startFirst() : start();
     135      }
    141136
    142       /// \brief Run the algorithm. It returns \c true if a path is found
    143       /// from the source node to the target node.
    144       bool run() {
     137    private:
     138   
     139      // Execute the algorithm for the first time (the flow and potential
     140      // functions have to be identically zero).
     141      bool startFirst() {
    145142        HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
    146143        Heap heap(heap_cross_ref);
    147144        heap.push(_s, 0);
     
    151148        // Process nodes
    152149        while (!heap.empty() && heap.top() != _t) {
    153150          Node u = heap.top(), v;
    154           Length d = heap.prio() + _potential[u], nd;
     151          Length d = heap.prio(), dn;
    155152          _dist[u] = heap.prio();
     153          _proc_nodes.push_back(u);
    156154          heap.pop();
     155
     156          // Traverse outgoing arcs
     157          for (OutArcIt e(_graph, u); e != INVALID; ++e) {
     158            v = _graph.target(e);
     159            switch(heap.state(v)) {
     160              case Heap::PRE_HEAP:
     161                heap.push(v, d + _length[e]);
     162                _pred[v] = e;
     163                break;
     164              case Heap::IN_HEAP:
     165                dn = d + _length[e];
     166                if (dn < heap[v]) {
     167                  heap.decrease(v, dn);
     168                  _pred[v] = e;
     169                }
     170                break;
     171              case Heap::POST_HEAP:
     172                break;
     173            }
     174          }
     175        }
     176        if (heap.empty()) return false;
     177
     178        // Update potentials of processed nodes
     179        Length t_dist = heap.prio();
     180        for (int i = 0; i < int(_proc_nodes.size()); ++i)
     181          _pi[_proc_nodes[i]] = _dist[_proc_nodes[i]] - t_dist;
     182        return true;
     183      }
     184
     185      // Execute the algorithm.
     186      bool start() {
     187        HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
     188        Heap heap(heap_cross_ref);
     189        heap.push(_s, 0);
     190        _pred[_s] = INVALID;
     191        _proc_nodes.clear();
     192
     193        // Process nodes
     194        while (!heap.empty() && heap.top() != _t) {
     195          Node u = heap.top(), v;
     196          Length d = heap.prio() + _pi[u], dn;
     197          _dist[u] = heap.prio();
    157198          _proc_nodes.push_back(u);
     199          heap.pop();
    158200
    159201          // Traverse outgoing arcs
    160202          for (OutArcIt e(_graph, u); e != INVALID; ++e) {
    161203            if (_flow[e] == 0) {
    162204              v = _graph.target(e);
    163205              switch(heap.state(v)) {
    164               case Heap::PRE_HEAP:
    165                 heap.push(v, d + _length[e] - _potential[v]);
    166                 _pred[v] = e;
    167                 break;
    168               case Heap::IN_HEAP:
    169                 nd = d + _length[e] - _potential[v];
    170                 if (nd < heap[v]) {
    171                   heap.decrease(v, nd);
     206                case Heap::PRE_HEAP:
     207                  heap.push(v, d + _length[e] - _pi[v]);
    172208                  _pred[v] = e;
    173                 }
    174                 break;
    175               case Heap::POST_HEAP:
    176                 break;
     209                  break;
     210                case Heap::IN_HEAP:
     211                  dn = d + _length[e] - _pi[v];
     212                  if (dn < heap[v]) {
     213                    heap.decrease(v, dn);
     214                    _pred[v] = e;
     215                  }
     216                  break;
     217                case Heap::POST_HEAP:
     218                  break;
    177219              }
    178220            }
    179221          }
     
    183225            if (_flow[e] == 1) {
    184226              v = _graph.source(e);
    185227              switch(heap.state(v)) {
    186               case Heap::PRE_HEAP:
    187                 heap.push(v, d - _length[e] - _potential[v]);
    188                 _pred[v] = e;
    189                 break;
    190               case Heap::IN_HEAP:
    191                 nd = d - _length[e] - _potential[v];
    192                 if (nd < heap[v]) {
    193                   heap.decrease(v, nd);
     228                case Heap::PRE_HEAP:
     229                  heap.push(v, d - _length[e] - _pi[v]);
    194230                  _pred[v] = e;
    195                 }
    196                 break;
    197               case Heap::POST_HEAP:
    198                 break;
     231                  break;
     232                case Heap::IN_HEAP:
     233                  dn = d - _length[e] - _pi[v];
     234                  if (dn < heap[v]) {
     235                    heap.decrease(v, dn);
     236                    _pred[v] = e;
     237                  }
     238                  break;
     239                case Heap::POST_HEAP:
     240                  break;
    199241              }
    200242            }
    201243          }
     
    205247        // Update potentials of processed nodes
    206248        Length t_dist = heap.prio();
    207249        for (int i = 0; i < int(_proc_nodes.size()); ++i)
    208           _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
     250          _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
    209251        return true;
    210252      }
    211253
     
    226268    bool _local_potential;
    227269
    228270    // The source node
    229     Node _source;
     271    Node _s;
    230272    // The target node
    231     Node _target;
     273    Node _t;
    232274
    233275    // Container to store the found paths
    234     std::vector< SimplePath<Digraph> > paths;
     276    std::vector<Path> _paths;
    235277    int _path_num;
    236278
    237279    // The pred arc map
    238280    PredMap _pred;
    239     // Implementation of the Dijkstra algorithm for finding augmenting
    240     // shortest paths in the residual network
    241     ResidualDijkstra *_dijkstra;
    242281
    243282  public:
    244283
     
    258297    ~Suurballe() {
    259298      if (_local_flow) delete _flow;
    260299      if (_local_potential) delete _potential;
    261       delete _dijkstra;
    262300    }
    263301
    264302    /// \brief Set the flow map.
     
    342380    ///
    343381    /// \param s The source node.
    344382    void init(const Node& s) {
    345       _source = s;
     383      _s = s;
    346384
    347385      // Initialize maps
    348386      if (!_flow) {
     
    372410    ///
    373411    /// \pre \ref init() must be called before using this function.
    374412    int findFlow(const Node& t, int k = 2) {
    375       _target = t;
    376       _dijkstra =
    377         new ResidualDijkstra( _graph, *_flow, _length, *_potential, _pred,
    378                               _source, _target );
     413      _t = t;
     414      ResidualDijkstra dijkstra(*this);
    379415
    380416      // Find shortest paths
    381417      _path_num = 0;
    382418      while (_path_num < k) {
    383419        // Run Dijkstra
    384         if (!_dijkstra->run()) break;
     420        if (!dijkstra.run(_path_num)) break;
    385421        ++_path_num;
    386422
    387423        // Set the flow along the found shortest path
    388         Node u = _target;
     424        Node u = _t;
    389425        Arc e;
    390426        while ((e = _pred[u]) != INVALID) {
    391427          if (u == _graph.target(e)) {
     
    402438
    403439    /// \brief Compute the paths from the flow.
    404440    ///
    405     /// This function computes the paths from the found minimum cost flow,
    406     /// which is the union of some arc-disjoint paths.
     441    /// This function computes arc-disjoint paths from the found minimum
     442    /// cost flow, which is the union of them.
    407443    ///
    408444    /// \pre \ref init() and \ref findFlow() must be called before using
    409445    /// this function.
     
    411447      FlowMap res_flow(_graph);
    412448      for(ArcIt a(_graph); a != INVALID; ++a) res_flow[a] = (*_flow)[a];
    413449
    414       paths.clear();
    415       paths.resize(_path_num);
     450      _paths.clear();
     451      _paths.resize(_path_num);
    416452      for (int i = 0; i < _path_num; ++i) {
    417         Node n = _source;
    418         while (n != _target) {
     453        Node n = _s;
     454        while (n != _t) {
    419455          OutArcIt e(_graph, n);
    420456          for ( ; res_flow[e] == 0; ++e) ;
    421457          n = _graph.target(e);
    422           paths[i].addBack(e);
     458          _paths[i].addBack(e);
    423459          res_flow[e] = 0;
    424460        }
    425461      }
     
    518554    /// \pre \ref run() or \ref findPaths() must be called before using
    519555    /// this function.
    520556    const Path& path(int i) const {
    521       return paths[i];
     557      return _paths[i];
    522558    }
    523559
    524560    /// @}