COIN-OR::LEMON - Graph Library

Ticket #60: 3774c5a08cb0.patch

File 3774c5a08cb0.patch, 29.5 KB (added by Balazs Dezso, 12 years ago)

Port min cost arborescence

  • lemon/Makefile.am

    # HG changeset patch
    # User Balazs Dezso <deba@inf.elte.hu>
    # Date 1228257227 -3600
    # Node ID 3774c5a08cb0cb58eecdd7b007a3f2ede9ab7599
    # Parent  ad483acf1654bd4d0c093b807b858ed59529f78e
    Port MinCostArborescence algorithm from SVN #3509
    
    diff -r ad483acf1654 -r 3774c5a08cb0 lemon/Makefile.am
    a b  
    4444        lemon/maps.h \
    4545        lemon/math.h \
    4646        lemon/max_matching.h \
     47        lemon/min_cost_arborescence.h \
    4748        lemon/nauty_reader.h \
    4849        lemon/path.h \
    4950        lemon/preflow.h \
  • new file lemon/min_cost_arborescence.h

    diff -r ad483acf1654 -r 3774c5a08cb0 lemon/min_cost_arborescence.h
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#ifndef LEMON_MIN_COST_ARBORESCENCE_H
     20#define LEMON_MIN_COST_ARBORESCENCE_H
     21
     22///\ingroup spantree
     23///\file
     24///\brief Minimum Cost Arborescence algorithm.
     25
     26#include <vector>
     27
     28#include <lemon/list_graph.h>
     29#include <lemon/bin_heap.h>
     30#include <lemon/assert.h>
     31
     32namespace lemon {
     33
     34
     35  /// \brief Default traits class for MinCostArborescence class.
     36  ///
     37  /// Default traits class for MinCostArborescence class.
     38  /// \param _Digraph Digraph type.
     39  /// \param _CostMap Type of cost map.
     40  template <class _Digraph, class _CostMap>
     41  struct MinCostArborescenceDefaultTraits{
     42
     43    /// \brief The digraph type the algorithm runs on.
     44    typedef _Digraph Digraph;
     45
     46    /// \brief The type of the map that stores the arc costs.
     47    ///
     48    /// The type of the map that stores the arc costs.
     49    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
     50    typedef _CostMap CostMap;
     51
     52    /// \brief The value type of the costs.
     53    ///
     54    /// The value type of the costs.
     55    typedef typename CostMap::Value Value;
     56
     57    /// \brief The type of the map that stores which arcs are in the
     58    /// arborescence.
     59    ///
     60    /// The type of the map that stores which arcs are in the
     61    /// arborescence.  It must meet the \ref concepts::WriteMap
     62    /// "WriteMap" concept.  Initially it will be set to false on each
     63    /// arc. After it will set all arborescence arcs once.
     64    typedef typename Digraph::template ArcMap<bool> ArborescenceMap;
     65
     66    /// \brief Instantiates a ArborescenceMap.
     67    ///
     68    /// This function instantiates a \ref ArborescenceMap.
     69    /// \param digraph is the graph, to which we would like to
     70    /// calculate the ArborescenceMap.
     71    static ArborescenceMap *createArborescenceMap(const Digraph &digraph){
     72      return new ArborescenceMap(digraph);
     73    }
     74
     75    /// \brief The type of the PredMap
     76    ///
     77    /// The type of the PredMap. It is a node map with an arc value type.
     78    typedef typename Digraph::template NodeMap<typename Digraph::Arc> PredMap;
     79
     80    /// \brief Instantiates a PredMap.
     81    ///
     82    /// This function instantiates a \ref PredMap.
     83    /// \param _digraph is the digraph, to which we would like to define the
     84    /// PredMap.
     85    static PredMap *createPredMap(const Digraph &digraph){
     86      return new PredMap(digraph);
     87    }
     88
     89  };
     90
     91  /// \ingroup spantree
     92  ///
     93  /// \brief %MinCostArborescence algorithm class.
     94  ///
     95  /// This class provides an efficient implementation of
     96  /// %MinCostArborescence algorithm. The arborescence is a tree
     97  /// which is directed from a given source node of the digraph. One or
     98  /// more sources should be given for the algorithm and it will calculate
     99  /// the minimum cost subgraph which are union of arborescences with the
     100  /// given sources and spans all the nodes which are reachable from the
     101  /// sources. The time complexity of the algorithm is \f$ O(n^2+e) \f$.
     102  ///
     103  /// The algorithm provides also an optimal dual solution, therefore
     104  /// the optimality of the solution can be checked.
     105  ///
     106  /// \param _Digraph The digraph type the algorithm runs on. The default value
     107  /// is \ref ListDigraph.
     108  /// \param _CostMap This read-only ArcMap determines the costs of the
     109  /// arcs. It is read once for each arc, so the map may involve in
     110  /// relatively time consuming process to compute the arc cost if
     111  /// it is necessary. The default map type is \ref
     112  /// concepts::Digraph::ArcMap "Digraph::ArcMap<int>".
     113  /// \param _Traits Traits class to set various data types used
     114  /// by the algorithm. The default traits class is
     115  /// \ref MinCostArborescenceDefaultTraits
     116  /// "MinCostArborescenceDefaultTraits<_Digraph, _CostMap>".  See \ref
     117  /// MinCostArborescenceDefaultTraits for the documentation of a
     118  /// MinCostArborescence traits class.
     119  ///
     120  /// \author Balazs Dezso
     121#ifndef DOXYGEN
     122  template <typename _Digraph = ListDigraph,
     123            typename _CostMap = typename _Digraph::template ArcMap<int>,
     124            typename _Traits =
     125            MinCostArborescenceDefaultTraits<_Digraph, _CostMap> >
     126#else
     127  template <typename _Digraph, typename _CostMap, typedef _Traits>
     128#endif
     129  class MinCostArborescence {
     130  public:
     131
     132    /// The traits.
     133    typedef _Traits Traits;
     134    /// The type of the underlying digraph.
     135    typedef typename Traits::Digraph Digraph;
     136    /// The type of the map that stores the arc costs.
     137    typedef typename Traits::CostMap CostMap;
     138    ///The type of the costs of the arcs.
     139    typedef typename Traits::Value Value;
     140    ///The type of the predecessor map.
     141    typedef typename Traits::PredMap PredMap;
     142    ///The type of the map that stores which arcs are in the arborescence.
     143    typedef typename Traits::ArborescenceMap ArborescenceMap;
     144
     145    typedef MinCostArborescence Create;
     146
     147  private:
     148
     149    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
     150
     151    struct CostArc {
     152
     153      Arc arc;
     154      Value value;
     155
     156      CostArc() {}
     157      CostArc(Arc _arc, Value _value) : arc(_arc), value(_value) {}
     158
     159    };
     160
     161    const Digraph *_digraph;
     162    const CostMap *_cost;
     163
     164    PredMap *_pred;
     165    bool local_pred;
     166
     167    ArborescenceMap *_arborescence;
     168    bool local_arborescence;
     169
     170    typedef typename Digraph::template ArcMap<int> ArcOrder;
     171    ArcOrder *_arc_order;
     172
     173    typedef typename Digraph::template NodeMap<int> NodeOrder;
     174    NodeOrder *_node_order;
     175
     176    typedef typename Digraph::template NodeMap<CostArc> CostArcMap;
     177    CostArcMap *_cost_arcs;
     178
     179    struct StackLevel {
     180
     181      std::vector<CostArc> arcs;
     182      int node_level;
     183
     184    };
     185
     186    std::vector<StackLevel> level_stack;
     187    std::vector<Node> queue;
     188
     189    typedef std::vector<typename Digraph::Node> DualNodeList;
     190
     191    DualNodeList _dual_node_list;
     192
     193    struct DualVariable {
     194      int begin, end;
     195      Value value;
     196
     197      DualVariable(int _begin, int _end, Value _value)
     198        : begin(_begin), end(_end), value(_value) {}
     199
     200    };
     201
     202    typedef std::vector<DualVariable> DualVariables;
     203
     204    DualVariables _dual_variables;
     205
     206    typedef typename Digraph::template NodeMap<int> HeapCrossRef;
     207
     208    HeapCrossRef *_heap_cross_ref;
     209
     210    typedef BinHeap<int, HeapCrossRef> Heap;
     211
     212    Heap *_heap;
     213
     214  protected:
     215
     216    MinCostArborescence() {}
     217
     218  private:
     219
     220    void createStructures() {
     221      if (!_pred) {
     222        local_pred = true;
     223        _pred = Traits::createPredMap(*_digraph);
     224      }
     225      if (!_arborescence) {
     226        local_arborescence = true;
     227        _arborescence = Traits::createArborescenceMap(*_digraph);
     228      }
     229      if (!_arc_order) {
     230        _arc_order = new ArcOrder(*_digraph);
     231      }
     232      if (!_node_order) {
     233        _node_order = new NodeOrder(*_digraph);
     234      }
     235      if (!_cost_arcs) {
     236        _cost_arcs = new CostArcMap(*_digraph);
     237      }
     238      if (!_heap_cross_ref) {
     239        _heap_cross_ref = new HeapCrossRef(*_digraph, -1);
     240      }
     241      if (!_heap) {
     242        _heap = new Heap(*_heap_cross_ref);
     243      }
     244    }
     245
     246    void destroyStructures() {
     247      if (local_arborescence) {
     248        delete _arborescence;
     249      }
     250      if (local_pred) {
     251        delete _pred;
     252      }
     253      if (_arc_order) {
     254        delete _arc_order;
     255      }
     256      if (_node_order) {
     257        delete _node_order;
     258      }
     259      if (_cost_arcs) {
     260        delete _cost_arcs;
     261      }
     262      if (_heap) {
     263        delete _heap;
     264      }
     265      if (_heap_cross_ref) {
     266        delete _heap_cross_ref;
     267      }
     268    }
     269
     270    Arc prepare(Node node) {
     271      std::vector<Node> nodes;
     272      (*_node_order)[node] = _dual_node_list.size();
     273      StackLevel level;
     274      level.node_level = _dual_node_list.size();
     275      _dual_node_list.push_back(node);
     276      for (InArcIt it(*_digraph, node); it != INVALID; ++it) {
     277        Arc arc = it;
     278        Node source = _digraph->source(arc);
     279        Value value = (*_cost)[it];
     280        if (source == node || (*_node_order)[source] == -3) continue;
     281        if ((*_cost_arcs)[source].arc == INVALID) {
     282          (*_cost_arcs)[source].arc = arc;
     283          (*_cost_arcs)[source].value = value;
     284          nodes.push_back(source);
     285        } else {
     286          if ((*_cost_arcs)[source].value > value) {
     287            (*_cost_arcs)[source].arc = arc;
     288            (*_cost_arcs)[source].value = value;
     289          }
     290        }
     291      }
     292      CostArc minimum = (*_cost_arcs)[nodes[0]];
     293      for (int i = 1; i < int(nodes.size()); ++i) {
     294        if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
     295          minimum = (*_cost_arcs)[nodes[i]];
     296        }
     297      }
     298      _arc_order->set(minimum.arc, _dual_variables.size());
     299      DualVariable var(_dual_node_list.size() - 1,
     300                       _dual_node_list.size(), minimum.value);
     301      _dual_variables.push_back(var);
     302      for (int i = 0; i < int(nodes.size()); ++i) {
     303        (*_cost_arcs)[nodes[i]].value -= minimum.value;
     304        level.arcs.push_back((*_cost_arcs)[nodes[i]]);
     305        (*_cost_arcs)[nodes[i]].arc = INVALID;
     306      }
     307      level_stack.push_back(level);
     308      return minimum.arc;
     309    }
     310
     311    Arc contract(Node node) {
     312      int node_bottom = bottom(node);
     313      std::vector<Node> nodes;
     314      while (!level_stack.empty() &&
     315             level_stack.back().node_level >= node_bottom) {
     316        for (int i = 0; i < int(level_stack.back().arcs.size()); ++i) {
     317          Arc arc = level_stack.back().arcs[i].arc;
     318          Node source = _digraph->source(arc);
     319          Value value = level_stack.back().arcs[i].value;
     320          if ((*_node_order)[source] >= node_bottom) continue;
     321          if ((*_cost_arcs)[source].arc == INVALID) {
     322            (*_cost_arcs)[source].arc = arc;
     323            (*_cost_arcs)[source].value = value;
     324            nodes.push_back(source);
     325          } else {
     326            if ((*_cost_arcs)[source].value > value) {
     327              (*_cost_arcs)[source].arc = arc;
     328              (*_cost_arcs)[source].value = value;
     329            }
     330          }
     331        }
     332        level_stack.pop_back();
     333      }
     334      CostArc minimum = (*_cost_arcs)[nodes[0]];
     335      for (int i = 1; i < int(nodes.size()); ++i) {
     336        if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
     337          minimum = (*_cost_arcs)[nodes[i]];
     338        }
     339      }
     340      _arc_order->set(minimum.arc, _dual_variables.size());
     341      DualVariable var(node_bottom, _dual_node_list.size(), minimum.value);
     342      _dual_variables.push_back(var);
     343      StackLevel level;
     344      level.node_level = node_bottom;
     345      for (int i = 0; i < int(nodes.size()); ++i) {
     346        (*_cost_arcs)[nodes[i]].value -= minimum.value;
     347        level.arcs.push_back((*_cost_arcs)[nodes[i]]);
     348        (*_cost_arcs)[nodes[i]].arc = INVALID;
     349      }
     350      level_stack.push_back(level);
     351      return minimum.arc;
     352    }
     353
     354    int bottom(Node node) {
     355      int k = level_stack.size() - 1;
     356      while (level_stack[k].node_level > (*_node_order)[node]) {
     357        --k;
     358      }
     359      return level_stack[k].node_level;
     360    }
     361
     362    void finalize(Arc arc) {
     363      Node node = _digraph->target(arc);
     364      _heap->push(node, (*_arc_order)[arc]);
     365      _pred->set(node, arc);
     366      while (!_heap->empty()) {
     367        Node source = _heap->top();
     368        _heap->pop();
     369        _node_order->set(source, -1);
     370        for (OutArcIt it(*_digraph, source); it != INVALID; ++it) {
     371          if ((*_arc_order)[it] < 0) continue;
     372          Node target = _digraph->target(it);
     373          switch(_heap->state(target)) {
     374          case Heap::PRE_HEAP:
     375            _heap->push(target, (*_arc_order)[it]);
     376            _pred->set(target, it);
     377            break;
     378          case Heap::IN_HEAP:
     379            if ((*_arc_order)[it] < (*_heap)[target]) {
     380              _heap->decrease(target, (*_arc_order)[it]);
     381              _pred->set(target, it);
     382            }
     383            break;
     384          case Heap::POST_HEAP:
     385            break;
     386          }
     387        }
     388        _arborescence->set((*_pred)[source], true);
     389      }
     390    }
     391
     392
     393  public:
     394
     395    /// \name Named template parameters
     396
     397    /// @{
     398
     399    template <class T>
     400    struct DefArborescenceMapTraits : public Traits {
     401      typedef T ArborescenceMap;
     402      static ArborescenceMap *createArborescenceMap(const Digraph &)
     403      {
     404        LEMON_ASSERT(false, "ArborescenceMap is not initialized");
     405        return 0; // ignore warnings
     406      }
     407    };
     408
     409    /// \brief \ref named-templ-param "Named parameter" for
     410    /// setting ArborescenceMap type
     411    ///
     412    /// \ref named-templ-param "Named parameter" for setting
     413    /// ArborescenceMap type
     414    template <class T>
     415    struct DefArborescenceMap
     416      : public MinCostArborescence<Digraph, CostMap,
     417                                   DefArborescenceMapTraits<T> > {
     418    };
     419
     420    template <class T>
     421    struct DefPredMapTraits : public Traits {
     422      typedef T PredMap;
     423      static PredMap *createPredMap(const Digraph &)
     424      {
     425        LEMON_ASSERT(false, "PredMap is not initialized");
     426      }
     427    };
     428
     429    /// \brief \ref named-templ-param "Named parameter" for
     430    /// setting PredMap type
     431    ///
     432    /// \ref named-templ-param "Named parameter" for setting
     433    /// PredMap type
     434    template <class T>
     435    struct DefPredMap
     436      : public MinCostArborescence<Digraph, CostMap, DefPredMapTraits<T> > {
     437    };
     438
     439    /// @}
     440
     441    /// \brief Constructor.
     442    ///
     443    /// \param _digraph The digraph the algorithm will run on.
     444    /// \param _cost The cost map used by the algorithm.
     445    MinCostArborescence(const Digraph& digraph, const CostMap& cost)
     446      : _digraph(&digraph), _cost(&cost), _pred(0), local_pred(false),
     447        _arborescence(0), local_arborescence(false),
     448        _arc_order(0), _node_order(0), _cost_arcs(0),
     449        _heap_cross_ref(0), _heap(0) {}
     450
     451    /// \brief Destructor.
     452    ~MinCostArborescence() {
     453      destroyStructures();
     454    }
     455
     456    /// \brief Sets the arborescence map.
     457    ///
     458    /// Sets the arborescence map.
     459    /// \return \c (*this)
     460    MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
     461      if (local_arborescence) {
     462        delete _arborescence;
     463      }
     464      local_arborescence = false;
     465      _arborescence = &m;
     466      return *this;
     467    }
     468
     469    /// \brief Sets the arborescence map.
     470    ///
     471    /// Sets the arborescence map.
     472    /// \return \c (*this)
     473    MinCostArborescence& predMap(PredMap& m) {
     474      if (local_pred) {
     475        delete _pred;
     476      }
     477      local_pred = false;
     478      _pred = &m;
     479      return *this;
     480    }
     481
     482    /// \name Query Functions
     483    /// The result of the %MinCostArborescence algorithm can be obtained
     484    /// using these functions.\n
     485    /// Before the use of these functions,
     486    /// either run() or start() must be called.
     487
     488    /// @{
     489
     490    /// \brief Returns a reference to the arborescence map.
     491    ///
     492    /// Returns a reference to the arborescence map.
     493    const ArborescenceMap& arborescenceMap() const {
     494      return *_arborescence;
     495    }
     496
     497    /// \brief Returns true if the arc is in the arborescence.
     498    ///
     499    /// Returns true if the arc is in the arborescence.
     500    /// \param arc The arc of the digraph.
     501    /// \pre \ref run() must be called before using this function.
     502    bool arborescence(Arc arc) const {
     503      return (*_pred)[_digraph->target(arc)] == arc;
     504    }
     505
     506    /// \brief Returns a reference to the pred map.
     507    ///
     508    /// Returns a reference to the pred map.
     509    const PredMap& predMap() const {
     510      return *_pred;
     511    }
     512
     513    /// \brief Returns the predecessor arc of the given node.
     514    ///
     515    /// Returns the predecessor arc of the given node.
     516    Arc pred(Node node) const {
     517      return (*_pred)[node];
     518    }
     519
     520    /// \brief Returns the cost of the arborescence.
     521    ///
     522    /// Returns the cost of the arborescence.
     523    Value arborescenceValue() const {
     524      Value sum = 0;
     525      for (ArcIt it(*_digraph); it != INVALID; ++it) {
     526        if (arborescence(it)) {
     527          sum += (*_cost)[it];
     528        }
     529      }
     530      return sum;
     531    }
     532
     533    /// \brief Indicates that a node is reachable from the sources.
     534    ///
     535    /// Indicates that a node is reachable from the sources.
     536    bool reached(Node node) const {
     537      return (*_node_order)[node] != -3;
     538    }
     539
     540    /// \brief Indicates that a node is processed.
     541    ///
     542    /// Indicates that a node is processed. The arborescence path exists
     543    /// from the source to the given node.
     544    bool processed(Node node) const {
     545      return (*_node_order)[node] == -1;
     546    }
     547
     548    /// \brief Returns the number of the dual variables in basis.
     549    ///
     550    /// Returns the number of the dual variables in basis.
     551    int dualNum() const {
     552      return _dual_variables.size();
     553    }
     554
     555    /// \brief Returns the value of the dual solution.
     556    ///
     557    /// Returns the value of the dual solution. It should be
     558    /// equal to the arborescence value.
     559    Value dualValue() const {
     560      Value sum = 0;
     561      for (int i = 0; i < int(_dual_variables.size()); ++i) {
     562        sum += _dual_variables[i].value;
     563      }
     564      return sum;
     565    }
     566
     567    /// \brief Returns the number of the nodes in the dual variable.
     568    ///
     569    /// Returns the number of the nodes in the dual variable.
     570    int dualSize(int k) const {
     571      return _dual_variables[k].end - _dual_variables[k].begin;
     572    }
     573
     574    /// \brief Returns the value of the dual variable.
     575    ///
     576    /// Returns the the value of the dual variable.
     577    const Value& dualValue(int k) const {
     578      return _dual_variables[k].value;
     579    }
     580
     581    /// \brief Lemon iterator for get a dual variable.
     582    ///
     583    /// Lemon iterator for get a dual variable. This class provides
     584    /// a common style lemon iterator which gives back a subset of
     585    /// the nodes.
     586    class DualIt {
     587    public:
     588
     589      /// \brief Constructor.
     590      ///
     591      /// Constructor for get the nodeset of the variable.
     592      DualIt(const MinCostArborescence& algorithm, int variable)
     593        : _algorithm(&algorithm)
     594      {
     595        _index = _algorithm->_dual_variables[variable].begin;
     596        _last = _algorithm->_dual_variables[variable].end;
     597      }
     598
     599      /// \brief Conversion to node.
     600      ///
     601      /// Conversion to node.
     602      operator Node() const {
     603        return _algorithm->_dual_node_list[_index];
     604      }
     605
     606      /// \brief Increment operator.
     607      ///
     608      /// Increment operator.
     609      DualIt& operator++() {
     610        ++_index;
     611        return *this;
     612      }
     613
     614      /// \brief Validity checking
     615      ///
     616      /// Checks whether the iterator is invalid.
     617      bool operator==(Invalid) const {
     618        return _index == _last;
     619      }
     620
     621      /// \brief Validity checking
     622      ///
     623      /// Checks whether the iterator is valid.
     624      bool operator!=(Invalid) const {
     625        return _index != _last;
     626      }
     627
     628    private:
     629      const MinCostArborescence* _algorithm;
     630      int _index, _last;
     631    };
     632
     633    /// @}
     634
     635    /// \name Execution control
     636    /// The simplest way to execute the algorithm is to use
     637    /// one of the member functions called \c run(...). \n
     638    /// If you need more control on the execution,
     639    /// first you must call \ref init(), then you can add several
     640    /// source nodes with \ref addSource().
     641    /// Finally \ref start() will perform the arborescence
     642    /// computation.
     643
     644    ///@{
     645
     646    /// \brief Initializes the internal data structures.
     647    ///
     648    /// Initializes the internal data structures.
     649    ///
     650    void init() {
     651      createStructures();
     652      _heap->clear();
     653      for (NodeIt it(*_digraph); it != INVALID; ++it) {
     654        (*_cost_arcs)[it].arc = INVALID;
     655        _node_order->set(it, -3);
     656        _heap_cross_ref->set(it, Heap::PRE_HEAP);
     657        _pred->set(it, INVALID);
     658      }
     659      for (ArcIt it(*_digraph); it != INVALID; ++it) {
     660        _arborescence->set(it, false);
     661        _arc_order->set(it, -1);
     662      }
     663      _dual_node_list.clear();
     664      _dual_variables.clear();
     665    }
     666
     667    /// \brief Adds a new source node.
     668    ///
     669    /// Adds a new source node to the algorithm.
     670    void addSource(Node source) {
     671      std::vector<Node> nodes;
     672      nodes.push_back(source);
     673      while (!nodes.empty()) {
     674        Node node = nodes.back();
     675        nodes.pop_back();
     676        for (OutArcIt it(*_digraph, node); it != INVALID; ++it) {
     677          Node target = _digraph->target(it);
     678          if ((*_node_order)[target] == -3) {
     679            (*_node_order)[target] = -2;
     680            nodes.push_back(target);
     681            queue.push_back(target);
     682          }
     683        }
     684      }
     685      (*_node_order)[source] = -1;
     686    }
     687
     688    /// \brief Processes the next node in the priority queue.
     689    ///
     690    /// Processes the next node in the priority queue.
     691    ///
     692    /// \return The processed node.
     693    ///
     694    /// \warning The queue must not be empty!
     695    Node processNextNode() {
     696      Node node = queue.back();
     697      queue.pop_back();
     698      if ((*_node_order)[node] == -2) {
     699        Arc arc = prepare(node);
     700        Node source = _digraph->source(arc);
     701        while ((*_node_order)[source] != -1) {
     702          if ((*_node_order)[source] >= 0) {
     703            arc = contract(source);
     704          } else {
     705            arc = prepare(source);
     706          }
     707          source = _digraph->source(arc);
     708        }
     709        finalize(arc);
     710        level_stack.clear();
     711      }
     712      return node;
     713    }
     714
     715    /// \brief Returns the number of the nodes to be processed.
     716    ///
     717    /// Returns the number of the nodes to be processed.
     718    int queueSize() const {
     719      return queue.size();
     720    }
     721
     722    /// \brief Returns \c false if there are nodes to be processed.
     723    ///
     724    /// Returns \c false if there are nodes to be processed.
     725    bool emptyQueue() const {
     726      return queue.empty();
     727    }
     728
     729    /// \brief Executes the algorithm.
     730    ///
     731    /// Executes the algorithm.
     732    ///
     733    /// \pre init() must be called and at least one node should be added
     734    /// with addSource() before using this function.
     735    ///
     736    ///\note mca.start() is just a shortcut of the following code.
     737    ///\code
     738    ///while (!mca.emptyQueue()) {
     739    ///  mca.processNextNode();
     740    ///}
     741    ///\endcode
     742    void start() {
     743      while (!emptyQueue()) {
     744        processNextNode();
     745      }
     746    }
     747
     748    /// \brief Runs %MinCostArborescence algorithm from node \c s.
     749    ///
     750    /// This method runs the %MinCostArborescence algorithm from
     751    /// a root node \c s.
     752    ///
     753    /// \note mca.run(s) is just a shortcut of the following code.
     754    /// \code
     755    /// mca.init();
     756    /// mca.addSource(s);
     757    /// mca.start();
     758    /// \endcode
     759    void run(Node node) {
     760      init();
     761      addSource(node);
     762      start();
     763    }
     764
     765    ///@}
     766
     767  };
     768
     769  /// \ingroup spantree
     770  ///
     771  /// \brief Function type interface for MinCostArborescence algorithm.
     772  ///
     773  /// Function type interface for MinCostArborescence algorithm.
     774  /// \param digraph The Digraph that the algorithm runs on.
     775  /// \param cost The CostMap of the arcs.
     776  /// \param source The source of the arborescence.
     777  /// \retval arborescence The bool ArcMap which stores the arborescence.
     778  /// \return The cost of the arborescence.
     779  ///
     780  /// \sa MinCostArborescence
     781  template <typename Digraph, typename CostMap, typename ArborescenceMap>
     782  typename CostMap::Value minCostArborescence(const Digraph& digraph,
     783                                              const CostMap& cost,
     784                                              typename Digraph::Node source,
     785                                              ArborescenceMap& arborescence) {
     786    typename MinCostArborescence<Digraph, CostMap>
     787      ::template DefArborescenceMap<ArborescenceMap>
     788      ::Create mca(digraph, cost);
     789    mca.arborescenceMap(arborescence);
     790    mca.run(source);
     791    return mca.arborescenceValue();
     792  }
     793
     794}
     795
     796#endif
  • test/Makefile.am

    diff -r ad483acf1654 -r 3774c5a08cb0 test/Makefile.am
    a b  
    2525        test/hao_orlin_test \
    2626        test/maps_test \
    2727        test/max_matching_test \
     28        test/min_cost_arborescence_test \
    2829        test/random_test \
    2930        test/path_test \
    3031        test/preflow_test \
     
    5455test_hao_orlin_test_SOURCES = test/hao_orlin_test.cc
    5556test_maps_test_SOURCES = test/maps_test.cc
    5657test_max_matching_test_SOURCES = test/max_matching_test.cc
     58test_min_cost_arborescence_test_SOURCES = test/min_cost_arborescence_test.cc
    5759test_path_test_SOURCES = test/path_test.cc
    5860test_preflow_test_SOURCES = test/preflow_test.cc
    5961test_suurballe_test_SOURCES = test/suurballe_test.cc
  • new file test/min_cost_arborescence_test.cc

    diff -r ad483acf1654 -r 3774c5a08cb0 test/min_cost_arborescence_test.cc
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#include <iostream>
     20#include <set>
     21#include <vector>
     22#include <iterator>
     23
     24#include <lemon/smart_graph.h>
     25#include <lemon/min_cost_arborescence.h>
     26#include <lemon/lgf_reader.h>
     27
     28#include "test_tools.h"
     29
     30using namespace lemon;
     31using namespace std;
     32
     33const char test_lgf[] =
     34  "@nodes\n"
     35  "label\n"
     36  "0\n"
     37  "1\n"
     38  "2\n"
     39  "3\n"
     40  "4\n"
     41  "5\n"
     42  "6\n"
     43  "7\n"
     44  "8\n"
     45  "9\n"
     46  "@arcs\n"
     47  "     label  cost\n"
     48  "1 8  0      107\n"
     49  "0 3  1      70\n"
     50  "2 1  2      46\n"
     51  "4 1  3      28\n"
     52  "4 4  4      91\n"
     53  "3 9  5      76\n"
     54  "9 8  6      61\n"
     55  "8 1  7      39\n"
     56  "9 8  8      74\n"
     57  "8 0  9      39\n"
     58  "4 3  10     45\n"
     59  "2 2  11     34\n"
     60  "0 1  12     100\n"
     61  "6 3  13     95\n"
     62  "4 1  14     22\n"
     63  "1 1  15     31\n"
     64  "7 2  16     51\n"
     65  "2 6  17     29\n"
     66  "8 3  18     115\n"
     67  "6 9  19     32\n"
     68  "1 1  20     60\n"
     69  "0 3  21     40\n"
     70  "@attributes\n"
     71  "source 0\n";
     72
     73int main() {
     74  typedef SmartDigraph Digraph;
     75  DIGRAPH_TYPEDEFS(Digraph);
     76
     77  typedef Digraph::ArcMap<double> CostMap;
     78
     79  Digraph digraph;
     80  CostMap cost(digraph);
     81  Node source;
     82
     83  std::istringstream is(test_lgf);
     84  digraphReader(digraph, is).
     85    arcMap("cost", cost).
     86    node("source", source).run();
     87
     88  MinCostArborescence<Digraph, CostMap> mca(digraph, cost);
     89  mca.run(source);
     90
     91  vector<pair<double, set<Node> > > dualSolution(mca.dualNum());
     92
     93  for (int i = 0; i < mca.dualNum(); ++i) {
     94    dualSolution[i].first = mca.dualValue(i);
     95    for (MinCostArborescence<Digraph, CostMap>::DualIt it(mca, i);
     96         it != INVALID; ++it) {
     97      dualSolution[i].second.insert(it);
     98    }
     99  }
     100
     101  for (ArcIt it(digraph); it != INVALID; ++it) {
     102    if (mca.reached(digraph.source(it))) {
     103      double sum = 0.0;
     104      for (int i = 0; i < int(dualSolution.size()); ++i) {
     105        if (dualSolution[i].second.find(digraph.target(it))
     106            != dualSolution[i].second.end() &&
     107            dualSolution[i].second.find(digraph.source(it))
     108            == dualSolution[i].second.end()) {
     109          sum += dualSolution[i].first;
     110        }
     111      }
     112      if (mca.arborescence(it)) {
     113        check(sum == cost[it], "INVALID DUAL");
     114      }
     115      check(sum <= cost[it], "INVALID DUAL");
     116    }
     117  }
     118
     119
     120  check(mca.dualValue() == mca.arborescenceValue(), "INVALID DUAL");
     121
     122  check(mca.reached(source), "INVALID ARBORESCENCE");
     123  for (ArcIt a(digraph); a != INVALID; ++a) {
     124    check(!mca.reached(digraph.source(a)) ||
     125          mca.reached(digraph.target(a)), "INVALID ARBORESCENCE");
     126  }
     127
     128  for (NodeIt n(digraph); n != INVALID; ++n) {
     129    if (!mca.reached(n)) continue;
     130    int cnt = 0;
     131    for (InArcIt a(digraph, n); a != INVALID; ++a) {
     132      if (mca.arborescence(a)) {
     133        check(mca.pred(n) == a, "INVALID ARBORESCENCE");
     134        ++cnt;
     135      }
     136    }
     137    check((n == source ? cnt == 0 : cnt == 1), "INVALID ARBORESCENCE");
     138  }
     139
     140  Digraph::ArcMap<bool> arborescence(digraph);
     141  check(mca.arborescenceValue() ==
     142        minCostArborescence(digraph, cost, source, arborescence),
     143        "WRONG FUNCTION");
     144
     145  return 0;
     146}