COIN-OR::LEMON - Graph Library

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/suurballe.h

    r927 r670  
    3030#include <lemon/path.h>
    3131#include <lemon/list_graph.h>
    32 #include <lemon/dijkstra.h>
    3332#include <lemon/maps.h>
    3433
     
    4847  /// "minimum cost flow problem". This implementation is actually an
    4948  /// efficient specialized version of the \ref CapacityScaling
    50   /// "successive shortest path" algorithm directly for this problem.
     49  /// "Successive Shortest Path" algorithm directly for this problem.
    5150  /// Therefore this class provides query functions for flow values and
    5251  /// node potentials (the dual solution) just like the minimum cost flow
     
    5756  /// The default value is <tt>GR::ArcMap<int></tt>.
    5857  ///
    59   /// \warning Length values should be \e non-negative.
     58  /// \warning Length values should be \e non-negative \e integers.
    6059  ///
    61   /// \note For finding \e node-disjoint paths, this algorithm can be used
     60  /// \note For finding node-disjoint paths this algorithm can be used
    6261  /// along with the \ref SplitNodes adaptor.
    6362#ifdef DOXYGEN
     
    9998  private:
    10099
    101     typedef typename Digraph::template NodeMap<int> HeapCrossRef;
    102     typedef BinHeap<Length, HeapCrossRef> Heap;
    103 
    104100    // ResidualDijkstra is a special implementation of the
    105101    // Dijkstra algorithm for finding shortest paths in the
     
    109105    class ResidualDijkstra
    110106    {
     107      typedef typename Digraph::template NodeMap<int> HeapCrossRef;
     108      typedef BinHeap<Length, HeapCrossRef> Heap;
     109
    111110    private:
    112111
     112      // The digraph the algorithm runs on
    113113      const Digraph &_graph;
     114
     115      // The main maps
     116      const FlowMap &_flow;
    114117      const LengthMap &_length;
    115       const FlowMap &_flow;
    116       PotentialMap &_pi;
     118      PotentialMap &_potential;
     119
     120      // The distance map
     121      PotentialMap _dist;
     122      // The pred arc map
    117123      PredMap &_pred;
     124      // The processed (i.e. permanently labeled) nodes
     125      std::vector<Node> _proc_nodes;
     126
    118127      Node _s;
    119128      Node _t;
    120      
    121       PotentialMap _dist;
    122       std::vector<Node> _proc_nodes;
    123129
    124130    public:
    125131
    126       // Constructor
    127       ResidualDijkstra(Suurballe &srb) :
    128         _graph(srb._graph), _length(srb._length),
    129         _flow(*srb._flow), _pi(*srb._potential), _pred(srb._pred),
    130         _s(srb._s), _t(srb._t), _dist(_graph) {}
    131        
    132       // Run the algorithm and return true if a path is found
    133       // from the source node to the target node.
    134       bool run(int cnt) {
    135         return cnt == 0 ? startFirst() : start();
    136       }
    137 
    138     private:
    139    
    140       // Execute the algorithm for the first time (the flow and potential
    141       // functions have to be identically zero).
    142       bool startFirst() {
     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) {}
     141
     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() {
    143145        HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
    144146        Heap heap(heap_cross_ref);
     
    150152        while (!heap.empty() && heap.top() != _t) {
    151153          Node u = heap.top(), v;
    152           Length d = heap.prio(), dn;
     154          Length d = heap.prio() + _potential[u], nd;
    153155          _dist[u] = heap.prio();
     156          heap.pop();
    154157          _proc_nodes.push_back(u);
    155           heap.pop();
    156158
    157159          // Traverse outgoing arcs
    158160          for (OutArcIt e(_graph, u); e != INVALID; ++e) {
    159             v = _graph.target(e);
    160             switch(heap.state(v)) {
     161            if (_flow[e] == 0) {
     162              v = _graph.target(e);
     163              switch(heap.state(v)) {
    161164              case Heap::PRE_HEAP:
    162                 heap.push(v, d + _length[e]);
     165                heap.push(v, d + _length[e] - _potential[v]);
    163166                _pred[v] = e;
    164167                break;
    165168              case Heap::IN_HEAP:
    166                 dn = d + _length[e];
    167                 if (dn < heap[v]) {
    168                   heap.decrease(v, dn);
     169                nd = d + _length[e] - _potential[v];
     170                if (nd < heap[v]) {
     171                  heap.decrease(v, nd);
    169172                  _pred[v] = e;
    170173                }
     
    172175              case Heap::POST_HEAP:
    173176                break;
    174             }
    175           }
    176         }
    177         if (heap.empty()) return false;
    178 
    179         // Update potentials of processed nodes
    180         Length t_dist = heap.prio();
    181         for (int i = 0; i < int(_proc_nodes.size()); ++i)
    182           _pi[_proc_nodes[i]] = _dist[_proc_nodes[i]] - t_dist;
    183         return true;
    184       }
    185 
    186       // Execute the algorithm.
    187       bool start() {
    188         HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
    189         Heap heap(heap_cross_ref);
    190         heap.push(_s, 0);
    191         _pred[_s] = INVALID;
    192         _proc_nodes.clear();
    193 
    194         // Process nodes
    195         while (!heap.empty() && heap.top() != _t) {
    196           Node u = heap.top(), v;
    197           Length d = heap.prio() + _pi[u], dn;
    198           _dist[u] = heap.prio();
    199           _proc_nodes.push_back(u);
    200           heap.pop();
    201 
    202           // Traverse outgoing arcs
    203           for (OutArcIt e(_graph, u); e != INVALID; ++e) {
    204             if (_flow[e] == 0) {
    205               v = _graph.target(e);
    206               switch(heap.state(v)) {
    207                 case Heap::PRE_HEAP:
    208                   heap.push(v, d + _length[e] - _pi[v]);
    209                   _pred[v] = e;
    210                   break;
    211                 case Heap::IN_HEAP:
    212                   dn = d + _length[e] - _pi[v];
    213                   if (dn < heap[v]) {
    214                     heap.decrease(v, dn);
    215                     _pred[v] = e;
    216                   }
    217                   break;
    218                 case Heap::POST_HEAP:
    219                   break;
    220177              }
    221178            }
     
    227184              v = _graph.source(e);
    228185              switch(heap.state(v)) {
    229                 case Heap::PRE_HEAP:
    230                   heap.push(v, d - _length[e] - _pi[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);
    231194                  _pred[v] = e;
    232                   break;
    233                 case Heap::IN_HEAP:
    234                   dn = d - _length[e] - _pi[v];
    235                   if (dn < heap[v]) {
    236                     heap.decrease(v, dn);
    237                     _pred[v] = e;
    238                   }
    239                   break;
    240                 case Heap::POST_HEAP:
    241                   break;
     195                }
     196                break;
     197              case Heap::POST_HEAP:
     198                break;
    242199              }
    243200            }
     
    249206        Length t_dist = heap.prio();
    250207        for (int i = 0; i < int(_proc_nodes.size()); ++i)
    251           _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
     208          _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
    252209        return true;
    253210      }
     
    270227
    271228    // The source node
    272     Node _s;
     229    Node _source;
    273230    // The target node
    274     Node _t;
     231    Node _target;
    275232
    276233    // Container to store the found paths
    277     std::vector<Path> _paths;
     234    std::vector< SimplePath<Digraph> > paths;
    278235    int _path_num;
    279236
    280237    // The pred arc map
    281238    PredMap _pred;
    282    
    283     // Data for full init
    284     PotentialMap *_init_dist;
    285     PredMap *_init_pred;
    286     bool _full_init;
     239    // Implementation of the Dijkstra algorithm for finding augmenting
     240    // shortest paths in the residual network
     241    ResidualDijkstra *_dijkstra;
    287242
    288243  public:
     
    297252               const LengthMap &length ) :
    298253      _graph(graph), _length(length), _flow(0), _local_flow(false),
    299       _potential(0), _local_potential(false), _pred(graph),
    300       _init_dist(0), _init_pred(0)
    301     {}
     254      _potential(0), _local_potential(false), _pred(graph)
     255    {
     256      LEMON_ASSERT(std::numeric_limits<Length>::is_integer,
     257        "The length type of Suurballe must be integer");
     258    }
    302259
    303260    /// Destructor.
     
    305262      if (_local_flow) delete _flow;
    306263      if (_local_potential) delete _potential;
    307       delete _init_dist;
    308       delete _init_pred;
     264      delete _dijkstra;
    309265    }
    310266
     
    351307    /// \name Execution Control
    352308    /// The simplest way to execute the algorithm is to call the run()
    353     /// function.\n
    354     /// If you need to execute the algorithm many times using the same
    355     /// source node, then you may call fullInit() once and start()
    356     /// for each target node.\n
     309    /// function.
     310    /// \n
    357311    /// If you only need the flow that is the union of the found
    358     /// arc-disjoint paths, then you may call findFlow() instead of
    359     /// start().
     312    /// arc-disjoint paths, you may call init() and findFlow().
    360313
    361314    /// @{
     
    377330    /// \code
    378331    ///   s.init(s);
    379     ///   s.start(t, k);
     332    ///   s.findFlow(t, k);
     333    ///   s.findPaths();
    380334    /// \endcode
    381335    int run(const Node& s, const Node& t, int k = 2) {
    382336      init(s);
    383       start(t, k);
     337      findFlow(t, k);
     338      findPaths();
    384339      return _path_num;
    385340    }
     
    387342    /// \brief Initialize the algorithm.
    388343    ///
    389     /// This function initializes the algorithm with the given source node.
     344    /// This function initializes the algorithm.
    390345    ///
    391346    /// \param s The source node.
    392347    void init(const Node& s) {
    393       _s = s;
     348      _source = s;
    394349
    395350      // Initialize maps
     
    402357        _local_potential = true;
    403358      }
    404       _full_init = false;
    405     }
    406 
    407     /// \brief Initialize the algorithm and perform Dijkstra.
    408     ///
    409     /// This function initializes the algorithm and performs a full
    410     /// Dijkstra search from the given source node. It makes consecutive
    411     /// executions of \ref start() "start(t, k)" faster, since they
    412     /// have to perform %Dijkstra only k-1 times.
    413     ///
    414     /// This initialization is usually worth using instead of \ref init()
    415     /// if the algorithm is executed many times using the same source node.
    416     ///
    417     /// \param s The source node.
    418     void fullInit(const Node& s) {
    419       // Initialize maps
    420       init(s);
    421       if (!_init_dist) {
    422         _init_dist = new PotentialMap(_graph);
    423       }
    424       if (!_init_pred) {
    425         _init_pred = new PredMap(_graph);
    426       }
    427 
    428       // Run a full Dijkstra
    429       typename Dijkstra<Digraph, LengthMap>
    430         ::template SetStandardHeap<Heap>
    431         ::template SetDistMap<PotentialMap>
    432         ::template SetPredMap<PredMap>
    433         ::Create dijk(_graph, _length);
    434       dijk.distMap(*_init_dist).predMap(*_init_pred);
    435       dijk.run(s);
    436      
    437       _full_init = true;
    438     }
    439 
    440     /// \brief Execute the algorithm.
    441     ///
    442     /// This function executes the algorithm.
    443     ///
    444     /// \param t The target node.
    445     /// \param k The number of paths to be found.
    446     ///
    447     /// \return \c k if there are at least \c k arc-disjoint paths from
    448     /// \c s to \c t in the digraph. Otherwise it returns the number of
    449     /// arc-disjoint paths found.
    450     ///
    451     /// \note Apart from the return value, <tt>s.start(t, k)</tt> is
    452     /// just a shortcut of the following code.
    453     /// \code
    454     ///   s.findFlow(t, k);
    455     ///   s.findPaths();
    456     /// \endcode
    457     int start(const Node& t, int k = 2) {
    458       findFlow(t, k);
    459       findPaths();
    460       return _path_num;
     359      for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
     360      for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
    461361    }
    462362
     
    476376    /// \pre \ref init() must be called before using this function.
    477377    int findFlow(const Node& t, int k = 2) {
    478       _t = t;
    479       ResidualDijkstra dijkstra(*this);
    480      
    481       // Initialization
    482       for (ArcIt e(_graph); e != INVALID; ++e) {
    483         (*_flow)[e] = 0;
    484       }
    485       if (_full_init) {
    486         for (NodeIt n(_graph); n != INVALID; ++n) {
    487           (*_potential)[n] = (*_init_dist)[n];
    488         }
    489         Node u = _t;
    490         Arc e;
    491         while ((e = (*_init_pred)[u]) != INVALID) {
    492           (*_flow)[e] = 1;
    493           u = _graph.source(e);
    494         }       
    495         _path_num = 1;
    496       } else {
    497         for (NodeIt n(_graph); n != INVALID; ++n) {
    498           (*_potential)[n] = 0;
    499         }
    500         _path_num = 0;
    501       }
     378      _target = t;
     379      _dijkstra =
     380        new ResidualDijkstra( _graph, *_flow, _length, *_potential, _pred,
     381                              _source, _target );
    502382
    503383      // Find shortest paths
     384      _path_num = 0;
    504385      while (_path_num < k) {
    505386        // Run Dijkstra
    506         if (!dijkstra.run(_path_num)) break;
     387        if (!_dijkstra->run()) break;
    507388        ++_path_num;
    508389
    509390        // Set the flow along the found shortest path
    510         Node u = _t;
     391        Node u = _target;
    511392        Arc e;
    512393        while ((e = _pred[u]) != INVALID) {
     
    525406    /// \brief Compute the paths from the flow.
    526407    ///
    527     /// This function computes arc-disjoint paths from the found minimum
    528     /// cost flow, which is the union of them.
     408    /// This function computes the paths from the found minimum cost flow,
     409    /// which is the union of some arc-disjoint paths.
    529410    ///
    530411    /// \pre \ref init() and \ref findFlow() must be called before using
     
    534415      for(ArcIt a(_graph); a != INVALID; ++a) res_flow[a] = (*_flow)[a];
    535416
    536       _paths.clear();
    537       _paths.resize(_path_num);
     417      paths.clear();
     418      paths.resize(_path_num);
    538419      for (int i = 0; i < _path_num; ++i) {
    539         Node n = _s;
    540         while (n != _t) {
     420        Node n = _source;
     421        while (n != _target) {
    541422          OutArcIt e(_graph, n);
    542423          for ( ; res_flow[e] == 0; ++e) ;
    543424          n = _graph.target(e);
    544           _paths[i].addBack(e);
     425          paths[i].addBack(e);
    545426          res_flow[e] = 0;
    546427        }
     
    640521    /// \pre \ref run() or \ref findPaths() must be called before using
    641522    /// this function.
    642     const Path& path(int i) const {
    643       return _paths[i];
     523    Path path(int i) const {
     524      return paths[i];
    644525    }
    645526
Note: See TracChangeset for help on using the changeset viewer.