COIN-OR::LEMON - Graph Library

Ignore:
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • lemon/suurballe.h

    r670 r927  
    3030#include <lemon/path.h>
    3131#include <lemon/list_graph.h>
     32#include <lemon/dijkstra.h>
    3233#include <lemon/maps.h>
    3334
     
    4748  /// "minimum cost flow problem". This implementation is actually an
    4849  /// efficient specialized version of the \ref CapacityScaling
    49   /// "Successive Shortest Path" algorithm directly for this problem.
     50  /// "successive shortest path" algorithm directly for this problem.
    5051  /// Therefore this class provides query functions for flow values and
    5152  /// node potentials (the dual solution) just like the minimum cost flow
     
    5657  /// The default value is <tt>GR::ArcMap<int></tt>.
    5758  ///
    58   /// \warning Length values should be \e non-negative \e integers.
     59  /// \warning Length values should be \e non-negative.
    5960  ///
    60   /// \note For finding node-disjoint paths this algorithm can be used
     61  /// \note For finding \e node-disjoint paths, this algorithm can be used
    6162  /// along with the \ref SplitNodes adaptor.
    6263#ifdef DOXYGEN
     
    9899  private:
    99100
     101    typedef typename Digraph::template NodeMap<int> HeapCrossRef;
     102    typedef BinHeap<Length, HeapCrossRef> Heap;
     103
    100104    // ResidualDijkstra is a special implementation of the
    101105    // Dijkstra algorithm for finding shortest paths in the
     
    105109    class ResidualDijkstra
    106110    {
    107       typedef typename Digraph::template NodeMap<int> HeapCrossRef;
    108       typedef BinHeap<Length, HeapCrossRef> Heap;
    109 
    110111    private:
    111112
    112       // The digraph the algorithm runs on
    113113      const Digraph &_graph;
    114 
    115       // The main maps
     114      const LengthMap &_length;
    116115      const FlowMap &_flow;
    117       const LengthMap &_length;
    118       PotentialMap &_potential;
    119 
    120       // The distance map
    121       PotentialMap _dist;
    122       // The pred arc map
     116      PotentialMap &_pi;
    123117      PredMap &_pred;
    124       // The processed (i.e. permanently labeled) nodes
    125       std::vector<Node> _proc_nodes;
    126 
    127118      Node _s;
    128119      Node _t;
     120     
     121      PotentialMap _dist;
     122      std::vector<Node> _proc_nodes;
    129123
    130124    public:
    131125
    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() {
     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() {
    145143        HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);
    146144        Heap heap(heap_cross_ref);
     
    152150        while (!heap.empty() && heap.top() != _t) {
    153151          Node u = heap.top(), v;
    154           Length d = heap.prio() + _potential[u], nd;
     152          Length d = heap.prio(), dn;
    155153          _dist[u] = heap.prio();
     154          _proc_nodes.push_back(u);
    156155          heap.pop();
     156
     157          // Traverse outgoing arcs
     158          for (OutArcIt e(_graph, u); e != INVALID; ++e) {
     159            v = _graph.target(e);
     160            switch(heap.state(v)) {
     161              case Heap::PRE_HEAP:
     162                heap.push(v, d + _length[e]);
     163                _pred[v] = e;
     164                break;
     165              case Heap::IN_HEAP:
     166                dn = d + _length[e];
     167                if (dn < heap[v]) {
     168                  heap.decrease(v, dn);
     169                  _pred[v] = e;
     170                }
     171                break;
     172              case Heap::POST_HEAP:
     173                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();
    157199          _proc_nodes.push_back(u);
     200          heap.pop();
    158201
    159202          // Traverse outgoing arcs
     
    162205              v = _graph.target(e);
    163206              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);
     207                case Heap::PRE_HEAP:
     208                  heap.push(v, d + _length[e] - _pi[v]);
    172209                  _pred[v] = e;
    173                 }
    174                 break;
    175               case Heap::POST_HEAP:
    176                 break;
     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;
    177220              }
    178221            }
     
    184227              v = _graph.source(e);
    185228              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);
     229                case Heap::PRE_HEAP:
     230                  heap.push(v, d - _length[e] - _pi[v]);
    194231                  _pred[v] = e;
    195                 }
    196                 break;
    197               case Heap::POST_HEAP:
    198                 break;
     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;
    199242              }
    200243            }
     
    206249        Length t_dist = heap.prio();
    207250        for (int i = 0; i < int(_proc_nodes.size()); ++i)
    208           _potential[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
     251          _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - t_dist;
    209252        return true;
    210253      }
     
    227270
    228271    // The source node
    229     Node _source;
     272    Node _s;
    230273    // The target node
    231     Node _target;
     274    Node _t;
    232275
    233276    // Container to store the found paths
    234     std::vector< SimplePath<Digraph> > paths;
     277    std::vector<Path> _paths;
    235278    int _path_num;
    236279
    237280    // The pred arc map
    238281    PredMap _pred;
    239     // Implementation of the Dijkstra algorithm for finding augmenting
    240     // shortest paths in the residual network
    241     ResidualDijkstra *_dijkstra;
     282   
     283    // Data for full init
     284    PotentialMap *_init_dist;
     285    PredMap *_init_pred;
     286    bool _full_init;
    242287
    243288  public:
     
    252297               const LengthMap &length ) :
    253298      _graph(graph), _length(length), _flow(0), _local_flow(false),
    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     }
     299      _potential(0), _local_potential(false), _pred(graph),
     300      _init_dist(0), _init_pred(0)
     301    {}
    259302
    260303    /// Destructor.
     
    262305      if (_local_flow) delete _flow;
    263306      if (_local_potential) delete _potential;
    264       delete _dijkstra;
     307      delete _init_dist;
     308      delete _init_pred;
    265309    }
    266310
     
    307351    /// \name Execution Control
    308352    /// The simplest way to execute the algorithm is to call the run()
    309     /// function.
    310     /// \n
     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
    311357    /// If you only need the flow that is the union of the found
    312     /// arc-disjoint paths, you may call init() and findFlow().
     358    /// arc-disjoint paths, then you may call findFlow() instead of
     359    /// start().
    313360
    314361    /// @{
     
    330377    /// \code
    331378    ///   s.init(s);
    332     ///   s.findFlow(t, k);
    333     ///   s.findPaths();
     379    ///   s.start(t, k);
    334380    /// \endcode
    335381    int run(const Node& s, const Node& t, int k = 2) {
    336382      init(s);
    337       findFlow(t, k);
    338       findPaths();
     383      start(t, k);
    339384      return _path_num;
    340385    }
     
    342387    /// \brief Initialize the algorithm.
    343388    ///
    344     /// This function initializes the algorithm.
     389    /// This function initializes the algorithm with the given source node.
    345390    ///
    346391    /// \param s The source node.
    347392    void init(const Node& s) {
    348       _source = s;
     393      _s = s;
    349394
    350395      // Initialize maps
     
    357402        _local_potential = true;
    358403      }
    359       for (ArcIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;
    360       for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;
     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;
    361461    }
    362462
     
    376476    /// \pre \ref init() must be called before using this function.
    377477    int findFlow(const Node& t, int k = 2) {
    378       _target = t;
    379       _dijkstra =
    380         new ResidualDijkstra( _graph, *_flow, _length, *_potential, _pred,
    381                               _source, _target );
     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      }
    382502
    383503      // Find shortest paths
    384       _path_num = 0;
    385504      while (_path_num < k) {
    386505        // Run Dijkstra
    387         if (!_dijkstra->run()) break;
     506        if (!dijkstra.run(_path_num)) break;
    388507        ++_path_num;
    389508
    390509        // Set the flow along the found shortest path
    391         Node u = _target;
     510        Node u = _t;
    392511        Arc e;
    393512        while ((e = _pred[u]) != INVALID) {
     
    406525    /// \brief Compute the paths from the flow.
    407526    ///
    408     /// This function computes the paths from the found minimum cost flow,
    409     /// which is the union of some arc-disjoint paths.
     527    /// This function computes arc-disjoint paths from the found minimum
     528    /// cost flow, which is the union of them.
    410529    ///
    411530    /// \pre \ref init() and \ref findFlow() must be called before using
     
    415534      for(ArcIt a(_graph); a != INVALID; ++a) res_flow[a] = (*_flow)[a];
    416535
    417       paths.clear();
    418       paths.resize(_path_num);
     536      _paths.clear();
     537      _paths.resize(_path_num);
    419538      for (int i = 0; i < _path_num; ++i) {
    420         Node n = _source;
    421         while (n != _target) {
     539        Node n = _s;
     540        while (n != _t) {
    422541          OutArcIt e(_graph, n);
    423542          for ( ; res_flow[e] == 0; ++e) ;
    424543          n = _graph.target(e);
    425           paths[i].addBack(e);
     544          _paths[i].addBack(e);
    426545          res_flow[e] = 0;
    427546        }
     
    521640    /// \pre \ref run() or \ref findPaths() must be called before using
    522641    /// this function.
    523     Path path(int i) const {
    524       return paths[i];
     642    const Path& path(int i) const {
     643      return _paths[i];
    525644    }
    526645
  • test/suurballe_test.cc

    r670 r927  
    102102  k = suurb_test.run(n, n, k);
    103103  suurb_test.init(n);
     104  suurb_test.fullInit(n);
     105  suurb_test.start(n);
     106  suurb_test.start(n, k);
    104107  k = suurb_test.findFlow(n);
    105108  k = suurb_test.findFlow(n, k);
Note: See TracChangeset for help on using the changeset viewer.