lemon/fractional_matching.h
author Balazs Dezso <deba@inf.elte.hu>
Sat, 26 Sep 2009 10:17:31 +0200
changeset 870 61120524af27
child 871 86613aa28a0c
permissions -rw-r--r--
Fractional matching initialization of weighted matchings (#314)
     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-2009
     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_FRACTIONAL_MATCHING_H
    20 #define LEMON_FRACTIONAL_MATCHING_H
    21 
    22 #include <vector>
    23 #include <queue>
    24 #include <set>
    25 #include <limits>
    26 
    27 #include <lemon/core.h>
    28 #include <lemon/unionfind.h>
    29 #include <lemon/bin_heap.h>
    30 #include <lemon/maps.h>
    31 #include <lemon/assert.h>
    32 #include <lemon/elevator.h>
    33 
    34 ///\ingroup matching
    35 ///\file
    36 ///\brief Fractional matching algorithms in general graphs.
    37 
    38 namespace lemon {
    39 
    40   /// \brief Default traits class of MaxFractionalMatching class.
    41   ///
    42   /// Default traits class of MaxFractionalMatching class.
    43   /// \tparam GR Graph type.
    44   template <typename GR>
    45   struct MaxFractionalMatchingDefaultTraits {
    46 
    47     /// \brief The type of the graph the algorithm runs on.
    48     typedef GR Graph;
    49 
    50     /// \brief The type of the map that stores the matching.
    51     ///
    52     /// The type of the map that stores the matching arcs.
    53     /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
    54     typedef typename Graph::template NodeMap<typename GR::Arc> MatchingMap;
    55 
    56     /// \brief Instantiates a MatchingMap.
    57     ///
    58     /// This function instantiates a \ref MatchingMap.
    59     /// \param graph The graph for which we would like to define
    60     /// the matching map.
    61     static MatchingMap* createMatchingMap(const Graph& graph) {
    62       return new MatchingMap(graph);
    63     }
    64 
    65     /// \brief The elevator type used by MaxFractionalMatching algorithm.
    66     ///
    67     /// The elevator type used by MaxFractionalMatching algorithm.
    68     ///
    69     /// \sa Elevator
    70     /// \sa LinkedElevator
    71     typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
    72 
    73     /// \brief Instantiates an Elevator.
    74     ///
    75     /// This function instantiates an \ref Elevator.
    76     /// \param graph The graph for which we would like to define
    77     /// the elevator.
    78     /// \param max_level The maximum level of the elevator.
    79     static Elevator* createElevator(const Graph& graph, int max_level) {
    80       return new Elevator(graph, max_level);
    81     }
    82   };
    83 
    84   /// \ingroup matching
    85   ///
    86   /// \brief Max cardinality fractional matching
    87   ///
    88   /// This class provides an implementation of fractional matching
    89   /// algorithm based on push-relabel principle.
    90   ///
    91   /// The maximum cardinality fractional matching is a relaxation of the
    92   /// maximum cardinality matching problem where the odd set constraints
    93   /// are omitted.
    94   /// It can be formulated with the following linear program.
    95   /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
    96   /// \f[x_e \ge 0\quad \forall e\in E\f]
    97   /// \f[\max \sum_{e\in E}x_e\f]
    98   /// where \f$\delta(X)\f$ is the set of edges incident to a node in
    99   /// \f$X\f$. The result can be represented as the union of a
   100   /// matching with one value edges and a set of odd length cycles
   101   /// with half value edges.
   102   ///
   103   /// The algorithm calculates an optimal fractional matching and a
   104   /// barrier. The number of adjacents of any node set minus the size
   105   /// of node set is a lower bound on the uncovered nodes in the
   106   /// graph. For maximum matching a barrier is computed which
   107   /// maximizes this difference.
   108   ///
   109   /// The algorithm can be executed with the run() function.  After it
   110   /// the matching (the primal solution) and the barrier (the dual
   111   /// solution) can be obtained using the query functions.
   112   ///
   113   /// The primal solution is multiplied by
   114   /// \ref MaxWeightedMatching::primalScale "2".
   115   ///
   116   /// \tparam GR The undirected graph type the algorithm runs on.
   117 #ifdef DOXYGEN
   118   template <typename GR, typename TR>
   119 #else
   120   template <typename GR,
   121             typename TR = MaxFractionalMatchingDefaultTraits<GR> >
   122 #endif
   123   class MaxFractionalMatching {
   124   public:
   125 
   126     /// \brief The \ref MaxFractionalMatchingDefaultTraits "traits
   127     /// class" of the algorithm.
   128     typedef TR Traits;
   129     /// The type of the graph the algorithm runs on.
   130     typedef typename TR::Graph Graph;
   131     /// The type of the matching map.
   132     typedef typename TR::MatchingMap MatchingMap;
   133     /// The type of the elevator.
   134     typedef typename TR::Elevator Elevator;
   135 
   136     /// \brief Scaling factor for primal solution
   137     ///
   138     /// Scaling factor for primal solution.
   139     static const int primalScale = 2;
   140 
   141   private:
   142 
   143     const Graph &_graph;
   144     int _node_num;
   145     bool _allow_loops;
   146     int _empty_level;
   147 
   148     TEMPLATE_GRAPH_TYPEDEFS(Graph);
   149 
   150     bool _local_matching;
   151     MatchingMap *_matching;
   152 
   153     bool _local_level;
   154     Elevator *_level;
   155 
   156     typedef typename Graph::template NodeMap<int> InDegMap;
   157     InDegMap *_indeg;
   158 
   159     void createStructures() {
   160       _node_num = countNodes(_graph);
   161 
   162       if (!_matching) {
   163         _local_matching = true;
   164         _matching = Traits::createMatchingMap(_graph);
   165       }
   166       if (!_level) {
   167         _local_level = true;
   168         _level = Traits::createElevator(_graph, _node_num);
   169       }
   170       if (!_indeg) {
   171         _indeg = new InDegMap(_graph);
   172       }
   173     }
   174 
   175     void destroyStructures() {
   176       if (_local_matching) {
   177         delete _matching;
   178       }
   179       if (_local_level) {
   180         delete _level;
   181       }
   182       if (_indeg) {
   183         delete _indeg;
   184       }
   185     }
   186 
   187     void postprocessing() {
   188       for (NodeIt n(_graph); n != INVALID; ++n) {
   189         if ((*_indeg)[n] != 0) continue;
   190         _indeg->set(n, -1);
   191         Node u = n;
   192         while ((*_matching)[u] != INVALID) {
   193           Node v = _graph.target((*_matching)[u]);
   194           _indeg->set(v, -1);
   195           Arc a = _graph.oppositeArc((*_matching)[u]);
   196           u = _graph.target((*_matching)[v]);
   197           _indeg->set(u, -1);
   198           _matching->set(v, a);
   199         }
   200       }
   201 
   202       for (NodeIt n(_graph); n != INVALID; ++n) {
   203         if ((*_indeg)[n] != 1) continue;
   204         _indeg->set(n, -1);
   205 
   206         int num = 1;
   207         Node u = _graph.target((*_matching)[n]);
   208         while (u != n) {
   209           _indeg->set(u, -1);
   210           u = _graph.target((*_matching)[u]);
   211           ++num;
   212         }
   213         if (num % 2 == 0 && num > 2) {
   214           Arc prev = _graph.oppositeArc((*_matching)[n]);
   215           Node v = _graph.target((*_matching)[n]);
   216           u = _graph.target((*_matching)[v]);
   217           _matching->set(v, prev);
   218           while (u != n) {
   219             prev = _graph.oppositeArc((*_matching)[u]);
   220             v = _graph.target((*_matching)[u]);
   221             u = _graph.target((*_matching)[v]);
   222             _matching->set(v, prev);
   223           }
   224         }
   225       }
   226     }
   227 
   228   public:
   229 
   230     typedef MaxFractionalMatching Create;
   231 
   232     ///\name Named Template Parameters
   233 
   234     ///@{
   235 
   236     template <typename T>
   237     struct SetMatchingMapTraits : public Traits {
   238       typedef T MatchingMap;
   239       static MatchingMap *createMatchingMap(const Graph&) {
   240         LEMON_ASSERT(false, "MatchingMap is not initialized");
   241         return 0; // ignore warnings
   242       }
   243     };
   244 
   245     /// \brief \ref named-templ-param "Named parameter" for setting
   246     /// MatchingMap type
   247     ///
   248     /// \ref named-templ-param "Named parameter" for setting MatchingMap
   249     /// type.
   250     template <typename T>
   251     struct SetMatchingMap
   252       : public MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > {
   253       typedef MaxFractionalMatching<Graph, SetMatchingMapTraits<T> > Create;
   254     };
   255 
   256     template <typename T>
   257     struct SetElevatorTraits : public Traits {
   258       typedef T Elevator;
   259       static Elevator *createElevator(const Graph&, int) {
   260         LEMON_ASSERT(false, "Elevator is not initialized");
   261         return 0; // ignore warnings
   262       }
   263     };
   264 
   265     /// \brief \ref named-templ-param "Named parameter" for setting
   266     /// Elevator type
   267     ///
   268     /// \ref named-templ-param "Named parameter" for setting Elevator
   269     /// type. If this named parameter is used, then an external
   270     /// elevator object must be passed to the algorithm using the
   271     /// \ref elevator(Elevator&) "elevator()" function before calling
   272     /// \ref run() or \ref init().
   273     /// \sa SetStandardElevator
   274     template <typename T>
   275     struct SetElevator
   276       : public MaxFractionalMatching<Graph, SetElevatorTraits<T> > {
   277       typedef MaxFractionalMatching<Graph, SetElevatorTraits<T> > Create;
   278     };
   279 
   280     template <typename T>
   281     struct SetStandardElevatorTraits : public Traits {
   282       typedef T Elevator;
   283       static Elevator *createElevator(const Graph& graph, int max_level) {
   284         return new Elevator(graph, max_level);
   285       }
   286     };
   287 
   288     /// \brief \ref named-templ-param "Named parameter" for setting
   289     /// Elevator type with automatic allocation
   290     ///
   291     /// \ref named-templ-param "Named parameter" for setting Elevator
   292     /// type with automatic allocation.
   293     /// The Elevator should have standard constructor interface to be
   294     /// able to automatically created by the algorithm (i.e. the
   295     /// graph and the maximum level should be passed to it).
   296     /// However an external elevator object could also be passed to the
   297     /// algorithm with the \ref elevator(Elevator&) "elevator()" function
   298     /// before calling \ref run() or \ref init().
   299     /// \sa SetElevator
   300     template <typename T>
   301     struct SetStandardElevator
   302       : public MaxFractionalMatching<Graph, SetStandardElevatorTraits<T> > {
   303       typedef MaxFractionalMatching<Graph,
   304                                     SetStandardElevatorTraits<T> > Create;
   305     };
   306 
   307     /// @}
   308 
   309   protected:
   310 
   311     MaxFractionalMatching() {}
   312 
   313   public:
   314 
   315     /// \brief Constructor
   316     ///
   317     /// Constructor.
   318     ///
   319     MaxFractionalMatching(const Graph &graph, bool allow_loops = true)
   320       : _graph(graph), _allow_loops(allow_loops),
   321         _local_matching(false), _matching(0),
   322         _local_level(false), _level(0),  _indeg(0)
   323     {}
   324 
   325     ~MaxFractionalMatching() {
   326       destroyStructures();
   327     }
   328 
   329     /// \brief Sets the matching map.
   330     ///
   331     /// Sets the matching map.
   332     /// If you don't use this function before calling \ref run() or
   333     /// \ref init(), an instance will be allocated automatically.
   334     /// The destructor deallocates this automatically allocated map,
   335     /// of course.
   336     /// \return <tt>(*this)</tt>
   337     MaxFractionalMatching& matchingMap(MatchingMap& map) {
   338       if (_local_matching) {
   339         delete _matching;
   340         _local_matching = false;
   341       }
   342       _matching = &map;
   343       return *this;
   344     }
   345 
   346     /// \brief Sets the elevator used by algorithm.
   347     ///
   348     /// Sets the elevator used by algorithm.
   349     /// If you don't use this function before calling \ref run() or
   350     /// \ref init(), an instance will be allocated automatically.
   351     /// The destructor deallocates this automatically allocated elevator,
   352     /// of course.
   353     /// \return <tt>(*this)</tt>
   354     MaxFractionalMatching& elevator(Elevator& elevator) {
   355       if (_local_level) {
   356         delete _level;
   357         _local_level = false;
   358       }
   359       _level = &elevator;
   360       return *this;
   361     }
   362 
   363     /// \brief Returns a const reference to the elevator.
   364     ///
   365     /// Returns a const reference to the elevator.
   366     ///
   367     /// \pre Either \ref run() or \ref init() must be called before
   368     /// using this function.
   369     const Elevator& elevator() const {
   370       return *_level;
   371     }
   372 
   373     /// \name Execution control
   374     /// The simplest way to execute the algorithm is to use one of the
   375     /// member functions called \c run(). \n
   376     /// If you need more control on the execution, first
   377     /// you must call \ref init() and then one variant of the start()
   378     /// member.
   379 
   380     /// @{
   381 
   382     /// \brief Initializes the internal data structures.
   383     ///
   384     /// Initializes the internal data structures and sets the initial
   385     /// matching.
   386     void init() {
   387       createStructures();
   388 
   389       _level->initStart();
   390       for (NodeIt n(_graph); n != INVALID; ++n) {
   391         _indeg->set(n, 0);
   392         _matching->set(n, INVALID);
   393         _level->initAddItem(n);
   394       }
   395       _level->initFinish();
   396 
   397       _empty_level = _node_num;
   398       for (NodeIt n(_graph); n != INVALID; ++n) {
   399         for (OutArcIt a(_graph, n); a != INVALID; ++a) {
   400           if (_graph.target(a) == n && !_allow_loops) continue;
   401           _matching->set(n, a);
   402           Node v = _graph.target((*_matching)[n]);
   403           _indeg->set(v, (*_indeg)[v] + 1);
   404           break;
   405         }
   406       }
   407 
   408       for (NodeIt n(_graph); n != INVALID; ++n) {
   409         if ((*_indeg)[n] == 0) {
   410           _level->activate(n);
   411         }
   412       }
   413     }
   414 
   415     /// \brief Starts the algorithm and computes a fractional matching
   416     ///
   417     /// The algorithm computes a maximum fractional matching.
   418     ///
   419     /// \param postprocess The algorithm computes first a matching
   420     /// which is a union of a matching with one value edges, cycles
   421     /// with half value edges and even length paths with half value
   422     /// edges. If the parameter is true, then after the push-relabel
   423     /// algorithm it postprocesses the matching to contain only
   424     /// matching edges and half value odd cycles.
   425     void start(bool postprocess = true) {
   426       Node n;
   427       while ((n = _level->highestActive()) != INVALID) {
   428         int level = _level->highestActiveLevel();
   429         int new_level = _level->maxLevel();
   430         for (InArcIt a(_graph, n); a != INVALID; ++a) {
   431           Node u = _graph.source(a);
   432           if (n == u && !_allow_loops) continue;
   433           Node v = _graph.target((*_matching)[u]);
   434           if ((*_level)[v] < level) {
   435             _indeg->set(v, (*_indeg)[v] - 1);
   436             if ((*_indeg)[v] == 0) {
   437               _level->activate(v);
   438             }
   439             _matching->set(u, a);
   440             _indeg->set(n, (*_indeg)[n] + 1);
   441             _level->deactivate(n);
   442             goto no_more_push;
   443           } else if (new_level > (*_level)[v]) {
   444             new_level = (*_level)[v];
   445           }
   446         }
   447 
   448         if (new_level + 1 < _level->maxLevel()) {
   449           _level->liftHighestActive(new_level + 1);
   450         } else {
   451           _level->liftHighestActiveToTop();
   452         }
   453         if (_level->emptyLevel(level)) {
   454           _level->liftToTop(level);
   455         }
   456       no_more_push:
   457         ;
   458       }
   459       for (NodeIt n(_graph); n != INVALID; ++n) {
   460         if ((*_matching)[n] == INVALID) continue;
   461         Node u = _graph.target((*_matching)[n]);
   462         if ((*_indeg)[u] > 1) {
   463           _indeg->set(u, (*_indeg)[u] - 1);
   464           _matching->set(n, INVALID);
   465         }
   466       }
   467       if (postprocess) {
   468         postprocessing();
   469       }
   470     }
   471 
   472     /// \brief Starts the algorithm and computes a perfect fractional
   473     /// matching
   474     ///
   475     /// The algorithm computes a perfect fractional matching. If it
   476     /// does not exists, then the algorithm returns false and the
   477     /// matching is undefined and the barrier.
   478     ///
   479     /// \param postprocess The algorithm computes first a matching
   480     /// which is a union of a matching with one value edges, cycles
   481     /// with half value edges and even length paths with half value
   482     /// edges. If the parameter is true, then after the push-relabel
   483     /// algorithm it postprocesses the matching to contain only
   484     /// matching edges and half value odd cycles.
   485     bool startPerfect(bool postprocess = true) {
   486       Node n;
   487       while ((n = _level->highestActive()) != INVALID) {
   488         int level = _level->highestActiveLevel();
   489         int new_level = _level->maxLevel();
   490         for (InArcIt a(_graph, n); a != INVALID; ++a) {
   491           Node u = _graph.source(a);
   492           if (n == u && !_allow_loops) continue;
   493           Node v = _graph.target((*_matching)[u]);
   494           if ((*_level)[v] < level) {
   495             _indeg->set(v, (*_indeg)[v] - 1);
   496             if ((*_indeg)[v] == 0) {
   497               _level->activate(v);
   498             }
   499             _matching->set(u, a);
   500             _indeg->set(n, (*_indeg)[n] + 1);
   501             _level->deactivate(n);
   502             goto no_more_push;
   503           } else if (new_level > (*_level)[v]) {
   504             new_level = (*_level)[v];
   505           }
   506         }
   507 
   508         if (new_level + 1 < _level->maxLevel()) {
   509           _level->liftHighestActive(new_level + 1);
   510         } else {
   511           _level->liftHighestActiveToTop();
   512           _empty_level = _level->maxLevel() - 1;
   513           return false;
   514         }
   515         if (_level->emptyLevel(level)) {
   516           _level->liftToTop(level);
   517           _empty_level = level;
   518           return false;
   519         }
   520       no_more_push:
   521         ;
   522       }
   523       if (postprocess) {
   524         postprocessing();
   525       }
   526       return true;
   527     }
   528 
   529     /// \brief Runs the algorithm
   530     ///
   531     /// Just a shortcut for the next code:
   532     ///\code
   533     /// init();
   534     /// start();
   535     ///\endcode
   536     void run(bool postprocess = true) {
   537       init();
   538       start(postprocess);
   539     }
   540 
   541     /// \brief Runs the algorithm to find a perfect fractional matching
   542     ///
   543     /// Just a shortcut for the next code:
   544     ///\code
   545     /// init();
   546     /// startPerfect();
   547     ///\endcode
   548     bool runPerfect(bool postprocess = true) {
   549       init();
   550       return startPerfect(postprocess);
   551     }
   552 
   553     ///@}
   554 
   555     /// \name Query Functions
   556     /// The result of the %Matching algorithm can be obtained using these
   557     /// functions.\n
   558     /// Before the use of these functions,
   559     /// either run() or start() must be called.
   560     ///@{
   561 
   562 
   563     /// \brief Return the number of covered nodes in the matching.
   564     ///
   565     /// This function returns the number of covered nodes in the matching.
   566     ///
   567     /// \pre Either run() or start() must be called before using this function.
   568     int matchingSize() const {
   569       int num = 0;
   570       for (NodeIt n(_graph); n != INVALID; ++n) {
   571         if ((*_matching)[n] != INVALID) {
   572           ++num;
   573         }
   574       }
   575       return num;
   576     }
   577 
   578     /// \brief Returns a const reference to the matching map.
   579     ///
   580     /// Returns a const reference to the node map storing the found
   581     /// fractional matching. This method can be called after
   582     /// running the algorithm.
   583     ///
   584     /// \pre Either \ref run() or \ref init() must be called before
   585     /// using this function.
   586     const MatchingMap& matchingMap() const {
   587       return *_matching;
   588     }
   589 
   590     /// \brief Return \c true if the given edge is in the matching.
   591     ///
   592     /// This function returns \c true if the given edge is in the
   593     /// found matching. The result is scaled by \ref primalScale
   594     /// "primal scale".
   595     ///
   596     /// \pre Either run() or start() must be called before using this function.
   597     int matching(const Edge& edge) const {
   598       return (edge == (*_matching)[_graph.u(edge)] ? 1 : 0) +
   599         (edge == (*_matching)[_graph.v(edge)] ? 1 : 0);
   600     }
   601 
   602     /// \brief Return the fractional matching arc (or edge) incident
   603     /// to the given node.
   604     ///
   605     /// This function returns one of the fractional matching arc (or
   606     /// edge) incident to the given node in the found matching or \c
   607     /// INVALID if the node is not covered by the matching or if the
   608     /// node is on an odd length cycle then it is the successor edge
   609     /// on the cycle.
   610     ///
   611     /// \pre Either run() or start() must be called before using this function.
   612     Arc matching(const Node& node) const {
   613       return (*_matching)[node];
   614     }
   615 
   616     /// \brief Returns true if the node is in the barrier
   617     ///
   618     /// The barrier is a subset of the nodes. If the nodes in the
   619     /// barrier have less adjacent nodes than the size of the barrier,
   620     /// then at least as much nodes cannot be covered as the
   621     /// difference of the two subsets.
   622     bool barrier(const Node& node) const {
   623       return (*_level)[node] >= _empty_level;
   624     }
   625 
   626     /// @}
   627 
   628   };
   629 
   630   /// \ingroup matching
   631   ///
   632   /// \brief Weighted fractional matching in general graphs
   633   ///
   634   /// This class provides an efficient implementation of fractional
   635   /// matching algorithm. The implementation is based on extensive use
   636   /// of priority queues and provides \f$O(nm\log n)\f$ time
   637   /// complexity.
   638   ///
   639   /// The maximum weighted fractional matching is a relaxation of the
   640   /// maximum weighted matching problem where the odd set constraints
   641   /// are omitted.
   642   /// It can be formulated with the following linear program.
   643   /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
   644   /// \f[x_e \ge 0\quad \forall e\in E\f]
   645   /// \f[\max \sum_{e\in E}x_ew_e\f]
   646   /// where \f$\delta(X)\f$ is the set of edges incident to a node in
   647   /// \f$X\f$. The result must be the union of a matching with one
   648   /// value edges and a set of odd length cycles with half value edges.
   649   ///
   650   /// The algorithm calculates an optimal fractional matching and a
   651   /// proof of the optimality. The solution of the dual problem can be
   652   /// used to check the result of the algorithm. The dual linear
   653   /// problem is the following.
   654   /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
   655   /// \f[y_u \ge 0 \quad \forall u \in V\f]
   656   /// \f[\min \sum_{u \in V}y_u \f] ///
   657   ///
   658   /// The algorithm can be executed with the run() function.
   659   /// After it the matching (the primal solution) and the dual solution
   660   /// can be obtained using the query functions.
   661   ///
   662   /// If the value type is integer, then the primal and the dual
   663   /// solutions are multiplied by
   664   /// \ref MaxWeightedMatching::primalScale "2" and
   665   /// \ref MaxWeightedMatching::dualScale "4" respectively.
   666   ///
   667   /// \tparam GR The undirected graph type the algorithm runs on.
   668   /// \tparam WM The type edge weight map. The default type is
   669   /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
   670 #ifdef DOXYGEN
   671   template <typename GR, typename WM>
   672 #else
   673   template <typename GR,
   674             typename WM = typename GR::template EdgeMap<int> >
   675 #endif
   676   class MaxWeightedFractionalMatching {
   677   public:
   678 
   679     /// The graph type of the algorithm
   680     typedef GR Graph;
   681     /// The type of the edge weight map
   682     typedef WM WeightMap;
   683     /// The value type of the edge weights
   684     typedef typename WeightMap::Value Value;
   685 
   686     /// The type of the matching map
   687     typedef typename Graph::template NodeMap<typename Graph::Arc>
   688     MatchingMap;
   689 
   690     /// \brief Scaling factor for primal solution
   691     ///
   692     /// Scaling factor for primal solution. It is equal to 2 or 1
   693     /// according to the value type.
   694     static const int primalScale =
   695       std::numeric_limits<Value>::is_integer ? 2 : 1;
   696 
   697     /// \brief Scaling factor for dual solution
   698     ///
   699     /// Scaling factor for dual solution. It is equal to 4 or 1
   700     /// according to the value type.
   701     static const int dualScale =
   702       std::numeric_limits<Value>::is_integer ? 4 : 1;
   703 
   704   private:
   705 
   706     TEMPLATE_GRAPH_TYPEDEFS(Graph);
   707 
   708     typedef typename Graph::template NodeMap<Value> NodePotential;
   709 
   710     const Graph& _graph;
   711     const WeightMap& _weight;
   712 
   713     MatchingMap* _matching;
   714     NodePotential* _node_potential;
   715 
   716     int _node_num;
   717     bool _allow_loops;
   718 
   719     enum Status {
   720       EVEN = -1, MATCHED = 0, ODD = 1
   721     };
   722 
   723     typedef typename Graph::template NodeMap<Status> StatusMap;
   724     StatusMap* _status;
   725 
   726     typedef typename Graph::template NodeMap<Arc> PredMap;
   727     PredMap* _pred;
   728 
   729     typedef ExtendFindEnum<IntNodeMap> TreeSet;
   730 
   731     IntNodeMap *_tree_set_index;
   732     TreeSet *_tree_set;
   733 
   734     IntNodeMap *_delta1_index;
   735     BinHeap<Value, IntNodeMap> *_delta1;
   736 
   737     IntNodeMap *_delta2_index;
   738     BinHeap<Value, IntNodeMap> *_delta2;
   739 
   740     IntEdgeMap *_delta3_index;
   741     BinHeap<Value, IntEdgeMap> *_delta3;
   742 
   743     Value _delta_sum;
   744 
   745     void createStructures() {
   746       _node_num = countNodes(_graph);
   747 
   748       if (!_matching) {
   749         _matching = new MatchingMap(_graph);
   750       }
   751       if (!_node_potential) {
   752         _node_potential = new NodePotential(_graph);
   753       }
   754       if (!_status) {
   755         _status = new StatusMap(_graph);
   756       }
   757       if (!_pred) {
   758         _pred = new PredMap(_graph);
   759       }
   760       if (!_tree_set) {
   761         _tree_set_index = new IntNodeMap(_graph);
   762         _tree_set = new TreeSet(*_tree_set_index);
   763       }
   764       if (!_delta1) {
   765         _delta1_index = new IntNodeMap(_graph);
   766         _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
   767       }
   768       if (!_delta2) {
   769         _delta2_index = new IntNodeMap(_graph);
   770         _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
   771       }
   772       if (!_delta3) {
   773         _delta3_index = new IntEdgeMap(_graph);
   774         _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
   775       }
   776     }
   777 
   778     void destroyStructures() {
   779       if (_matching) {
   780         delete _matching;
   781       }
   782       if (_node_potential) {
   783         delete _node_potential;
   784       }
   785       if (_status) {
   786         delete _status;
   787       }
   788       if (_pred) {
   789         delete _pred;
   790       }
   791       if (_tree_set) {
   792         delete _tree_set_index;
   793         delete _tree_set;
   794       }
   795       if (_delta1) {
   796         delete _delta1_index;
   797         delete _delta1;
   798       }
   799       if (_delta2) {
   800         delete _delta2_index;
   801         delete _delta2;
   802       }
   803       if (_delta3) {
   804         delete _delta3_index;
   805         delete _delta3;
   806       }
   807     }
   808 
   809     void matchedToEven(Node node, int tree) {
   810       _tree_set->insert(node, tree);
   811       _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
   812       _delta1->push(node, (*_node_potential)[node]);
   813 
   814       if (_delta2->state(node) == _delta2->IN_HEAP) {
   815         _delta2->erase(node);
   816       }
   817 
   818       for (InArcIt a(_graph, node); a != INVALID; ++a) {
   819         Node v = _graph.source(a);
   820         Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
   821           dualScale * _weight[a];
   822         if (node == v) {
   823           if (_allow_loops && _graph.direction(a)) {
   824             _delta3->push(a, rw / 2);
   825           }
   826         } else if ((*_status)[v] == EVEN) {
   827           _delta3->push(a, rw / 2);
   828         } else if ((*_status)[v] == MATCHED) {
   829           if (_delta2->state(v) != _delta2->IN_HEAP) {
   830             _pred->set(v, a);
   831             _delta2->push(v, rw);
   832           } else if ((*_delta2)[v] > rw) {
   833             _pred->set(v, a);
   834             _delta2->decrease(v, rw);
   835           }
   836         }
   837       }
   838     }
   839 
   840     void matchedToOdd(Node node, int tree) {
   841       _tree_set->insert(node, tree);
   842       _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
   843 
   844       if (_delta2->state(node) == _delta2->IN_HEAP) {
   845         _delta2->erase(node);
   846       }
   847     }
   848 
   849     void evenToMatched(Node node, int tree) {
   850       _delta1->erase(node);
   851       _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
   852       Arc min = INVALID;
   853       Value minrw = std::numeric_limits<Value>::max();
   854       for (InArcIt a(_graph, node); a != INVALID; ++a) {
   855         Node v = _graph.source(a);
   856         Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
   857           dualScale * _weight[a];
   858 
   859         if (node == v) {
   860           if (_allow_loops && _graph.direction(a)) {
   861             _delta3->erase(a);
   862           }
   863         } else if ((*_status)[v] == EVEN) {
   864           _delta3->erase(a);
   865           if (minrw > rw) {
   866             min = _graph.oppositeArc(a);
   867             minrw = rw;
   868           }
   869         } else if ((*_status)[v]  == MATCHED) {
   870           if ((*_pred)[v] == a) {
   871             Arc mina = INVALID;
   872             Value minrwa = std::numeric_limits<Value>::max();
   873             for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
   874               Node va = _graph.target(aa);
   875               if ((*_status)[va] != EVEN ||
   876                   _tree_set->find(va) == tree) continue;
   877               Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
   878                 dualScale * _weight[aa];
   879               if (minrwa > rwa) {
   880                 minrwa = rwa;
   881                 mina = aa;
   882               }
   883             }
   884             if (mina != INVALID) {
   885               _pred->set(v, mina);
   886               _delta2->increase(v, minrwa);
   887             } else {
   888               _pred->set(v, INVALID);
   889               _delta2->erase(v);
   890             }
   891           }
   892         }
   893       }
   894       if (min != INVALID) {
   895         _pred->set(node, min);
   896         _delta2->push(node, minrw);
   897       } else {
   898         _pred->set(node, INVALID);
   899       }
   900     }
   901 
   902     void oddToMatched(Node node) {
   903       _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
   904       Arc min = INVALID;
   905       Value minrw = std::numeric_limits<Value>::max();
   906       for (InArcIt a(_graph, node); a != INVALID; ++a) {
   907         Node v = _graph.source(a);
   908         if ((*_status)[v] != EVEN) continue;
   909         Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
   910           dualScale * _weight[a];
   911 
   912         if (minrw > rw) {
   913           min = _graph.oppositeArc(a);
   914           minrw = rw;
   915         }
   916       }
   917       if (min != INVALID) {
   918         _pred->set(node, min);
   919         _delta2->push(node, minrw);
   920       } else {
   921         _pred->set(node, INVALID);
   922       }
   923     }
   924 
   925     void alternatePath(Node even, int tree) {
   926       Node odd;
   927 
   928       _status->set(even, MATCHED);
   929       evenToMatched(even, tree);
   930 
   931       Arc prev = (*_matching)[even];
   932       while (prev != INVALID) {
   933         odd = _graph.target(prev);
   934         even = _graph.target((*_pred)[odd]);
   935         _matching->set(odd, (*_pred)[odd]);
   936         _status->set(odd, MATCHED);
   937         oddToMatched(odd);
   938 
   939         prev = (*_matching)[even];
   940         _status->set(even, MATCHED);
   941         _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
   942         evenToMatched(even, tree);
   943       }
   944     }
   945 
   946     void destroyTree(int tree) {
   947       for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
   948         if ((*_status)[n] == EVEN) {
   949           _status->set(n, MATCHED);
   950           evenToMatched(n, tree);
   951         } else if ((*_status)[n] == ODD) {
   952           _status->set(n, MATCHED);
   953           oddToMatched(n);
   954         }
   955       }
   956       _tree_set->eraseClass(tree);
   957     }
   958 
   959 
   960     void unmatchNode(const Node& node) {
   961       int tree = _tree_set->find(node);
   962 
   963       alternatePath(node, tree);
   964       destroyTree(tree);
   965 
   966       _matching->set(node, INVALID);
   967     }
   968 
   969 
   970     void augmentOnEdge(const Edge& edge) {
   971       Node left = _graph.u(edge);
   972       int left_tree = _tree_set->find(left);
   973 
   974       alternatePath(left, left_tree);
   975       destroyTree(left_tree);
   976       _matching->set(left, _graph.direct(edge, true));
   977 
   978       Node right = _graph.v(edge);
   979       int right_tree = _tree_set->find(right);
   980 
   981       alternatePath(right, right_tree);
   982       destroyTree(right_tree);
   983       _matching->set(right, _graph.direct(edge, false));
   984     }
   985 
   986     void augmentOnArc(const Arc& arc) {
   987       Node left = _graph.source(arc);
   988       _status->set(left, MATCHED);
   989       _matching->set(left, arc);
   990       _pred->set(left, arc);
   991 
   992       Node right = _graph.target(arc);
   993       int right_tree = _tree_set->find(right);
   994 
   995       alternatePath(right, right_tree);
   996       destroyTree(right_tree);
   997       _matching->set(right, _graph.oppositeArc(arc));
   998     }
   999 
  1000     void extendOnArc(const Arc& arc) {
  1001       Node base = _graph.target(arc);
  1002       int tree = _tree_set->find(base);
  1003 
  1004       Node odd = _graph.source(arc);
  1005       _tree_set->insert(odd, tree);
  1006       _status->set(odd, ODD);
  1007       matchedToOdd(odd, tree);
  1008       _pred->set(odd, arc);
  1009 
  1010       Node even = _graph.target((*_matching)[odd]);
  1011       _tree_set->insert(even, tree);
  1012       _status->set(even, EVEN);
  1013       matchedToEven(even, tree);
  1014     }
  1015 
  1016     void cycleOnEdge(const Edge& edge, int tree) {
  1017       Node nca = INVALID;
  1018       std::vector<Node> left_path, right_path;
  1019 
  1020       {
  1021         std::set<Node> left_set, right_set;
  1022         Node left = _graph.u(edge);
  1023         left_path.push_back(left);
  1024         left_set.insert(left);
  1025 
  1026         Node right = _graph.v(edge);
  1027         right_path.push_back(right);
  1028         right_set.insert(right);
  1029 
  1030         while (true) {
  1031 
  1032           if (left_set.find(right) != left_set.end()) {
  1033             nca = right;
  1034             break;
  1035           }
  1036 
  1037           if ((*_matching)[left] == INVALID) break;
  1038 
  1039           left = _graph.target((*_matching)[left]);
  1040           left_path.push_back(left);
  1041           left = _graph.target((*_pred)[left]);
  1042           left_path.push_back(left);
  1043 
  1044           left_set.insert(left);
  1045 
  1046           if (right_set.find(left) != right_set.end()) {
  1047             nca = left;
  1048             break;
  1049           }
  1050 
  1051           if ((*_matching)[right] == INVALID) break;
  1052 
  1053           right = _graph.target((*_matching)[right]);
  1054           right_path.push_back(right);
  1055           right = _graph.target((*_pred)[right]);
  1056           right_path.push_back(right);
  1057 
  1058           right_set.insert(right);
  1059 
  1060         }
  1061 
  1062         if (nca == INVALID) {
  1063           if ((*_matching)[left] == INVALID) {
  1064             nca = right;
  1065             while (left_set.find(nca) == left_set.end()) {
  1066               nca = _graph.target((*_matching)[nca]);
  1067               right_path.push_back(nca);
  1068               nca = _graph.target((*_pred)[nca]);
  1069               right_path.push_back(nca);
  1070             }
  1071           } else {
  1072             nca = left;
  1073             while (right_set.find(nca) == right_set.end()) {
  1074               nca = _graph.target((*_matching)[nca]);
  1075               left_path.push_back(nca);
  1076               nca = _graph.target((*_pred)[nca]);
  1077               left_path.push_back(nca);
  1078             }
  1079           }
  1080         }
  1081       }
  1082 
  1083       alternatePath(nca, tree);
  1084       Arc prev;
  1085 
  1086       prev = _graph.direct(edge, true);
  1087       for (int i = 0; left_path[i] != nca; i += 2) {
  1088         _matching->set(left_path[i], prev);
  1089         _status->set(left_path[i], MATCHED);
  1090         evenToMatched(left_path[i], tree);
  1091 
  1092         prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
  1093         _status->set(left_path[i + 1], MATCHED);
  1094         oddToMatched(left_path[i + 1]);
  1095       }
  1096       _matching->set(nca, prev);
  1097 
  1098       for (int i = 0; right_path[i] != nca; i += 2) {
  1099         _status->set(right_path[i], MATCHED);
  1100         evenToMatched(right_path[i], tree);
  1101 
  1102         _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
  1103         _status->set(right_path[i + 1], MATCHED);
  1104         oddToMatched(right_path[i + 1]);
  1105       }
  1106 
  1107       destroyTree(tree);
  1108     }
  1109 
  1110     void extractCycle(const Arc &arc) {
  1111       Node left = _graph.source(arc);
  1112       Node odd = _graph.target((*_matching)[left]);
  1113       Arc prev;
  1114       while (odd != left) {
  1115         Node even = _graph.target((*_matching)[odd]);
  1116         prev = (*_matching)[odd];
  1117         odd = _graph.target((*_matching)[even]);
  1118         _matching->set(even, _graph.oppositeArc(prev));
  1119       }
  1120       _matching->set(left, arc);
  1121 
  1122       Node right = _graph.target(arc);
  1123       int right_tree = _tree_set->find(right);
  1124       alternatePath(right, right_tree);
  1125       destroyTree(right_tree);
  1126       _matching->set(right, _graph.oppositeArc(arc));
  1127     }
  1128 
  1129   public:
  1130 
  1131     /// \brief Constructor
  1132     ///
  1133     /// Constructor.
  1134     MaxWeightedFractionalMatching(const Graph& graph, const WeightMap& weight,
  1135                                   bool allow_loops = true)
  1136       : _graph(graph), _weight(weight), _matching(0),
  1137       _node_potential(0), _node_num(0), _allow_loops(allow_loops),
  1138       _status(0),  _pred(0),
  1139       _tree_set_index(0), _tree_set(0),
  1140 
  1141       _delta1_index(0), _delta1(0),
  1142       _delta2_index(0), _delta2(0),
  1143       _delta3_index(0), _delta3(0),
  1144 
  1145       _delta_sum() {}
  1146 
  1147     ~MaxWeightedFractionalMatching() {
  1148       destroyStructures();
  1149     }
  1150 
  1151     /// \name Execution Control
  1152     /// The simplest way to execute the algorithm is to use the
  1153     /// \ref run() member function.
  1154 
  1155     ///@{
  1156 
  1157     /// \brief Initialize the algorithm
  1158     ///
  1159     /// This function initializes the algorithm.
  1160     void init() {
  1161       createStructures();
  1162 
  1163       for (NodeIt n(_graph); n != INVALID; ++n) {
  1164         (*_delta1_index)[n] = _delta1->PRE_HEAP;
  1165         (*_delta2_index)[n] = _delta2->PRE_HEAP;
  1166       }
  1167       for (EdgeIt e(_graph); e != INVALID; ++e) {
  1168         (*_delta3_index)[e] = _delta3->PRE_HEAP;
  1169       }
  1170 
  1171       for (NodeIt n(_graph); n != INVALID; ++n) {
  1172         Value max = 0;
  1173         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  1174           if (_graph.target(e) == n && !_allow_loops) continue;
  1175           if ((dualScale * _weight[e]) / 2 > max) {
  1176             max = (dualScale * _weight[e]) / 2;
  1177           }
  1178         }
  1179         _node_potential->set(n, max);
  1180         _delta1->push(n, max);
  1181 
  1182         _tree_set->insert(n);
  1183 
  1184         _matching->set(n, INVALID);
  1185         _status->set(n, EVEN);
  1186       }
  1187 
  1188       for (EdgeIt e(_graph); e != INVALID; ++e) {
  1189         Node left = _graph.u(e);
  1190         Node right = _graph.v(e);
  1191         if (left == right && !_allow_loops) continue;
  1192         _delta3->push(e, ((*_node_potential)[left] +
  1193                           (*_node_potential)[right] -
  1194                           dualScale * _weight[e]) / 2);
  1195       }
  1196     }
  1197 
  1198     /// \brief Start the algorithm
  1199     ///
  1200     /// This function starts the algorithm.
  1201     ///
  1202     /// \pre \ref init() must be called before using this function.
  1203     void start() {
  1204       enum OpType {
  1205         D1, D2, D3
  1206       };
  1207 
  1208       int unmatched = _node_num;
  1209       while (unmatched > 0) {
  1210         Value d1 = !_delta1->empty() ?
  1211           _delta1->prio() : std::numeric_limits<Value>::max();
  1212 
  1213         Value d2 = !_delta2->empty() ?
  1214           _delta2->prio() : std::numeric_limits<Value>::max();
  1215 
  1216         Value d3 = !_delta3->empty() ?
  1217           _delta3->prio() : std::numeric_limits<Value>::max();
  1218 
  1219         _delta_sum = d3; OpType ot = D3;
  1220         if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
  1221         if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
  1222 
  1223         switch (ot) {
  1224         case D1:
  1225           {
  1226             Node n = _delta1->top();
  1227             unmatchNode(n);
  1228             --unmatched;
  1229           }
  1230           break;
  1231         case D2:
  1232           {
  1233             Node n = _delta2->top();
  1234             Arc a = (*_pred)[n];
  1235             if ((*_matching)[n] == INVALID) {
  1236               augmentOnArc(a);
  1237               --unmatched;
  1238             } else {
  1239               Node v = _graph.target((*_matching)[n]);
  1240               if ((*_matching)[n] !=
  1241                   _graph.oppositeArc((*_matching)[v])) {
  1242                 extractCycle(a);
  1243                 --unmatched;
  1244               } else {
  1245                 extendOnArc(a);
  1246               }
  1247             }
  1248           } break;
  1249         case D3:
  1250           {
  1251             Edge e = _delta3->top();
  1252 
  1253             Node left = _graph.u(e);
  1254             Node right = _graph.v(e);
  1255 
  1256             int left_tree = _tree_set->find(left);
  1257             int right_tree = _tree_set->find(right);
  1258 
  1259             if (left_tree == right_tree) {
  1260               cycleOnEdge(e, left_tree);
  1261               --unmatched;
  1262             } else {
  1263               augmentOnEdge(e);
  1264               unmatched -= 2;
  1265             }
  1266           } break;
  1267         }
  1268       }
  1269     }
  1270 
  1271     /// \brief Run the algorithm.
  1272     ///
  1273     /// This method runs the \c %MaxWeightedMatching algorithm.
  1274     ///
  1275     /// \note mwfm.run() is just a shortcut of the following code.
  1276     /// \code
  1277     ///   mwfm.init();
  1278     ///   mwfm.start();
  1279     /// \endcode
  1280     void run() {
  1281       init();
  1282       start();
  1283     }
  1284 
  1285     /// @}
  1286 
  1287     /// \name Primal Solution
  1288     /// Functions to get the primal solution, i.e. the maximum weighted
  1289     /// matching.\n
  1290     /// Either \ref run() or \ref start() function should be called before
  1291     /// using them.
  1292 
  1293     /// @{
  1294 
  1295     /// \brief Return the weight of the matching.
  1296     ///
  1297     /// This function returns the weight of the found matching. This
  1298     /// value is scaled by \ref primalScale "primal scale".
  1299     ///
  1300     /// \pre Either run() or start() must be called before using this function.
  1301     Value matchingWeight() const {
  1302       Value sum = 0;
  1303       for (NodeIt n(_graph); n != INVALID; ++n) {
  1304         if ((*_matching)[n] != INVALID) {
  1305           sum += _weight[(*_matching)[n]];
  1306         }
  1307       }
  1308       return sum * primalScale / 2;
  1309     }
  1310 
  1311     /// \brief Return the number of covered nodes in the matching.
  1312     ///
  1313     /// This function returns the number of covered nodes in the matching.
  1314     ///
  1315     /// \pre Either run() or start() must be called before using this function.
  1316     int matchingSize() const {
  1317       int num = 0;
  1318       for (NodeIt n(_graph); n != INVALID; ++n) {
  1319         if ((*_matching)[n] != INVALID) {
  1320           ++num;
  1321         }
  1322       }
  1323       return num;
  1324     }
  1325 
  1326     /// \brief Return \c true if the given edge is in the matching.
  1327     ///
  1328     /// This function returns \c true if the given edge is in the
  1329     /// found matching. The result is scaled by \ref primalScale
  1330     /// "primal scale".
  1331     ///
  1332     /// \pre Either run() or start() must be called before using this function.
  1333     Value matching(const Edge& edge) const {
  1334       return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
  1335         * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0)
  1336         * primalScale / 2;
  1337     }
  1338 
  1339     /// \brief Return the fractional matching arc (or edge) incident
  1340     /// to the given node.
  1341     ///
  1342     /// This function returns one of the fractional matching arc (or
  1343     /// edge) incident to the given node in the found matching or \c
  1344     /// INVALID if the node is not covered by the matching or if the
  1345     /// node is on an odd length cycle then it is the successor edge
  1346     /// on the cycle.
  1347     ///
  1348     /// \pre Either run() or start() must be called before using this function.
  1349     Arc matching(const Node& node) const {
  1350       return (*_matching)[node];
  1351     }
  1352 
  1353     /// \brief Return a const reference to the matching map.
  1354     ///
  1355     /// This function returns a const reference to a node map that stores
  1356     /// the matching arc (or edge) incident to each node.
  1357     const MatchingMap& matchingMap() const {
  1358       return *_matching;
  1359     }
  1360 
  1361     /// @}
  1362 
  1363     /// \name Dual Solution
  1364     /// Functions to get the dual solution.\n
  1365     /// Either \ref run() or \ref start() function should be called before
  1366     /// using them.
  1367 
  1368     /// @{
  1369 
  1370     /// \brief Return the value of the dual solution.
  1371     ///
  1372     /// This function returns the value of the dual solution.
  1373     /// It should be equal to the primal value scaled by \ref dualScale
  1374     /// "dual scale".
  1375     ///
  1376     /// \pre Either run() or start() must be called before using this function.
  1377     Value dualValue() const {
  1378       Value sum = 0;
  1379       for (NodeIt n(_graph); n != INVALID; ++n) {
  1380         sum += nodeValue(n);
  1381       }
  1382       return sum;
  1383     }
  1384 
  1385     /// \brief Return the dual value (potential) of the given node.
  1386     ///
  1387     /// This function returns the dual value (potential) of the given node.
  1388     ///
  1389     /// \pre Either run() or start() must be called before using this function.
  1390     Value nodeValue(const Node& n) const {
  1391       return (*_node_potential)[n];
  1392     }
  1393 
  1394     /// @}
  1395 
  1396   };
  1397 
  1398   /// \ingroup matching
  1399   ///
  1400   /// \brief Weighted fractional perfect matching in general graphs
  1401   ///
  1402   /// This class provides an efficient implementation of fractional
  1403   /// matching algorithm. The implementation is based on extensive use
  1404   /// of priority queues and provides \f$O(nm\log n)\f$ time
  1405   /// complexity.
  1406   ///
  1407   /// The maximum weighted fractional perfect matching is a relaxation
  1408   /// of the maximum weighted perfect matching problem where the odd
  1409   /// set constraints are omitted.
  1410   /// It can be formulated with the following linear program.
  1411   /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
  1412   /// \f[x_e \ge 0\quad \forall e\in E\f]
  1413   /// \f[\max \sum_{e\in E}x_ew_e\f]
  1414   /// where \f$\delta(X)\f$ is the set of edges incident to a node in
  1415   /// \f$X\f$. The result must be the union of a matching with one
  1416   /// value edges and a set of odd length cycles with half value edges.
  1417   ///
  1418   /// The algorithm calculates an optimal fractional matching and a
  1419   /// proof of the optimality. The solution of the dual problem can be
  1420   /// used to check the result of the algorithm. The dual linear
  1421   /// problem is the following.
  1422   /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
  1423   /// \f[\min \sum_{u \in V}y_u \f] ///
  1424   ///
  1425   /// The algorithm can be executed with the run() function.
  1426   /// After it the matching (the primal solution) and the dual solution
  1427   /// can be obtained using the query functions.
  1428 
  1429   /// If the value type is integer, then the primal and the dual
  1430   /// solutions are multiplied by
  1431   /// \ref MaxWeightedMatching::primalScale "2" and
  1432   /// \ref MaxWeightedMatching::dualScale "4" respectively.
  1433   ///
  1434   /// \tparam GR The undirected graph type the algorithm runs on.
  1435   /// \tparam WM The type edge weight map. The default type is
  1436   /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
  1437 #ifdef DOXYGEN
  1438   template <typename GR, typename WM>
  1439 #else
  1440   template <typename GR,
  1441             typename WM = typename GR::template EdgeMap<int> >
  1442 #endif
  1443   class MaxWeightedPerfectFractionalMatching {
  1444   public:
  1445 
  1446     /// The graph type of the algorithm
  1447     typedef GR Graph;
  1448     /// The type of the edge weight map
  1449     typedef WM WeightMap;
  1450     /// The value type of the edge weights
  1451     typedef typename WeightMap::Value Value;
  1452 
  1453     /// The type of the matching map
  1454     typedef typename Graph::template NodeMap<typename Graph::Arc>
  1455     MatchingMap;
  1456 
  1457     /// \brief Scaling factor for primal solution
  1458     ///
  1459     /// Scaling factor for primal solution. It is equal to 2 or 1
  1460     /// according to the value type.
  1461     static const int primalScale =
  1462       std::numeric_limits<Value>::is_integer ? 2 : 1;
  1463 
  1464     /// \brief Scaling factor for dual solution
  1465     ///
  1466     /// Scaling factor for dual solution. It is equal to 4 or 1
  1467     /// according to the value type.
  1468     static const int dualScale =
  1469       std::numeric_limits<Value>::is_integer ? 4 : 1;
  1470 
  1471   private:
  1472 
  1473     TEMPLATE_GRAPH_TYPEDEFS(Graph);
  1474 
  1475     typedef typename Graph::template NodeMap<Value> NodePotential;
  1476 
  1477     const Graph& _graph;
  1478     const WeightMap& _weight;
  1479 
  1480     MatchingMap* _matching;
  1481     NodePotential* _node_potential;
  1482 
  1483     int _node_num;
  1484     bool _allow_loops;
  1485 
  1486     enum Status {
  1487       EVEN = -1, MATCHED = 0, ODD = 1
  1488     };
  1489 
  1490     typedef typename Graph::template NodeMap<Status> StatusMap;
  1491     StatusMap* _status;
  1492 
  1493     typedef typename Graph::template NodeMap<Arc> PredMap;
  1494     PredMap* _pred;
  1495 
  1496     typedef ExtendFindEnum<IntNodeMap> TreeSet;
  1497 
  1498     IntNodeMap *_tree_set_index;
  1499     TreeSet *_tree_set;
  1500 
  1501     IntNodeMap *_delta2_index;
  1502     BinHeap<Value, IntNodeMap> *_delta2;
  1503 
  1504     IntEdgeMap *_delta3_index;
  1505     BinHeap<Value, IntEdgeMap> *_delta3;
  1506 
  1507     Value _delta_sum;
  1508 
  1509     void createStructures() {
  1510       _node_num = countNodes(_graph);
  1511 
  1512       if (!_matching) {
  1513         _matching = new MatchingMap(_graph);
  1514       }
  1515       if (!_node_potential) {
  1516         _node_potential = new NodePotential(_graph);
  1517       }
  1518       if (!_status) {
  1519         _status = new StatusMap(_graph);
  1520       }
  1521       if (!_pred) {
  1522         _pred = new PredMap(_graph);
  1523       }
  1524       if (!_tree_set) {
  1525         _tree_set_index = new IntNodeMap(_graph);
  1526         _tree_set = new TreeSet(*_tree_set_index);
  1527       }
  1528       if (!_delta2) {
  1529         _delta2_index = new IntNodeMap(_graph);
  1530         _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
  1531       }
  1532       if (!_delta3) {
  1533         _delta3_index = new IntEdgeMap(_graph);
  1534         _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
  1535       }
  1536     }
  1537 
  1538     void destroyStructures() {
  1539       if (_matching) {
  1540         delete _matching;
  1541       }
  1542       if (_node_potential) {
  1543         delete _node_potential;
  1544       }
  1545       if (_status) {
  1546         delete _status;
  1547       }
  1548       if (_pred) {
  1549         delete _pred;
  1550       }
  1551       if (_tree_set) {
  1552         delete _tree_set_index;
  1553         delete _tree_set;
  1554       }
  1555       if (_delta2) {
  1556         delete _delta2_index;
  1557         delete _delta2;
  1558       }
  1559       if (_delta3) {
  1560         delete _delta3_index;
  1561         delete _delta3;
  1562       }
  1563     }
  1564 
  1565     void matchedToEven(Node node, int tree) {
  1566       _tree_set->insert(node, tree);
  1567       _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
  1568 
  1569       if (_delta2->state(node) == _delta2->IN_HEAP) {
  1570         _delta2->erase(node);
  1571       }
  1572 
  1573       for (InArcIt a(_graph, node); a != INVALID; ++a) {
  1574         Node v = _graph.source(a);
  1575         Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
  1576           dualScale * _weight[a];
  1577         if (node == v) {
  1578           if (_allow_loops && _graph.direction(a)) {
  1579             _delta3->push(a, rw / 2);
  1580           }
  1581         } else if ((*_status)[v] == EVEN) {
  1582           _delta3->push(a, rw / 2);
  1583         } else if ((*_status)[v] == MATCHED) {
  1584           if (_delta2->state(v) != _delta2->IN_HEAP) {
  1585             _pred->set(v, a);
  1586             _delta2->push(v, rw);
  1587           } else if ((*_delta2)[v] > rw) {
  1588             _pred->set(v, a);
  1589             _delta2->decrease(v, rw);
  1590           }
  1591         }
  1592       }
  1593     }
  1594 
  1595     void matchedToOdd(Node node, int tree) {
  1596       _tree_set->insert(node, tree);
  1597       _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
  1598 
  1599       if (_delta2->state(node) == _delta2->IN_HEAP) {
  1600         _delta2->erase(node);
  1601       }
  1602     }
  1603 
  1604     void evenToMatched(Node node, int tree) {
  1605       _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
  1606       Arc min = INVALID;
  1607       Value minrw = std::numeric_limits<Value>::max();
  1608       for (InArcIt a(_graph, node); a != INVALID; ++a) {
  1609         Node v = _graph.source(a);
  1610         Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
  1611           dualScale * _weight[a];
  1612 
  1613         if (node == v) {
  1614           if (_allow_loops && _graph.direction(a)) {
  1615             _delta3->erase(a);
  1616           }
  1617         } else if ((*_status)[v] == EVEN) {
  1618           _delta3->erase(a);
  1619           if (minrw > rw) {
  1620             min = _graph.oppositeArc(a);
  1621             minrw = rw;
  1622           }
  1623         } else if ((*_status)[v]  == MATCHED) {
  1624           if ((*_pred)[v] == a) {
  1625             Arc mina = INVALID;
  1626             Value minrwa = std::numeric_limits<Value>::max();
  1627             for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
  1628               Node va = _graph.target(aa);
  1629               if ((*_status)[va] != EVEN ||
  1630                   _tree_set->find(va) == tree) continue;
  1631               Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
  1632                 dualScale * _weight[aa];
  1633               if (minrwa > rwa) {
  1634                 minrwa = rwa;
  1635                 mina = aa;
  1636               }
  1637             }
  1638             if (mina != INVALID) {
  1639               _pred->set(v, mina);
  1640               _delta2->increase(v, minrwa);
  1641             } else {
  1642               _pred->set(v, INVALID);
  1643               _delta2->erase(v);
  1644             }
  1645           }
  1646         }
  1647       }
  1648       if (min != INVALID) {
  1649         _pred->set(node, min);
  1650         _delta2->push(node, minrw);
  1651       } else {
  1652         _pred->set(node, INVALID);
  1653       }
  1654     }
  1655 
  1656     void oddToMatched(Node node) {
  1657       _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
  1658       Arc min = INVALID;
  1659       Value minrw = std::numeric_limits<Value>::max();
  1660       for (InArcIt a(_graph, node); a != INVALID; ++a) {
  1661         Node v = _graph.source(a);
  1662         if ((*_status)[v] != EVEN) continue;
  1663         Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
  1664           dualScale * _weight[a];
  1665 
  1666         if (minrw > rw) {
  1667           min = _graph.oppositeArc(a);
  1668           minrw = rw;
  1669         }
  1670       }
  1671       if (min != INVALID) {
  1672         _pred->set(node, min);
  1673         _delta2->push(node, minrw);
  1674       } else {
  1675         _pred->set(node, INVALID);
  1676       }
  1677     }
  1678 
  1679     void alternatePath(Node even, int tree) {
  1680       Node odd;
  1681 
  1682       _status->set(even, MATCHED);
  1683       evenToMatched(even, tree);
  1684 
  1685       Arc prev = (*_matching)[even];
  1686       while (prev != INVALID) {
  1687         odd = _graph.target(prev);
  1688         even = _graph.target((*_pred)[odd]);
  1689         _matching->set(odd, (*_pred)[odd]);
  1690         _status->set(odd, MATCHED);
  1691         oddToMatched(odd);
  1692 
  1693         prev = (*_matching)[even];
  1694         _status->set(even, MATCHED);
  1695         _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
  1696         evenToMatched(even, tree);
  1697       }
  1698     }
  1699 
  1700     void destroyTree(int tree) {
  1701       for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
  1702         if ((*_status)[n] == EVEN) {
  1703           _status->set(n, MATCHED);
  1704           evenToMatched(n, tree);
  1705         } else if ((*_status)[n] == ODD) {
  1706           _status->set(n, MATCHED);
  1707           oddToMatched(n);
  1708         }
  1709       }
  1710       _tree_set->eraseClass(tree);
  1711     }
  1712 
  1713     void augmentOnEdge(const Edge& edge) {
  1714       Node left = _graph.u(edge);
  1715       int left_tree = _tree_set->find(left);
  1716 
  1717       alternatePath(left, left_tree);
  1718       destroyTree(left_tree);
  1719       _matching->set(left, _graph.direct(edge, true));
  1720 
  1721       Node right = _graph.v(edge);
  1722       int right_tree = _tree_set->find(right);
  1723 
  1724       alternatePath(right, right_tree);
  1725       destroyTree(right_tree);
  1726       _matching->set(right, _graph.direct(edge, false));
  1727     }
  1728 
  1729     void augmentOnArc(const Arc& arc) {
  1730       Node left = _graph.source(arc);
  1731       _status->set(left, MATCHED);
  1732       _matching->set(left, arc);
  1733       _pred->set(left, arc);
  1734 
  1735       Node right = _graph.target(arc);
  1736       int right_tree = _tree_set->find(right);
  1737 
  1738       alternatePath(right, right_tree);
  1739       destroyTree(right_tree);
  1740       _matching->set(right, _graph.oppositeArc(arc));
  1741     }
  1742 
  1743     void extendOnArc(const Arc& arc) {
  1744       Node base = _graph.target(arc);
  1745       int tree = _tree_set->find(base);
  1746 
  1747       Node odd = _graph.source(arc);
  1748       _tree_set->insert(odd, tree);
  1749       _status->set(odd, ODD);
  1750       matchedToOdd(odd, tree);
  1751       _pred->set(odd, arc);
  1752 
  1753       Node even = _graph.target((*_matching)[odd]);
  1754       _tree_set->insert(even, tree);
  1755       _status->set(even, EVEN);
  1756       matchedToEven(even, tree);
  1757     }
  1758 
  1759     void cycleOnEdge(const Edge& edge, int tree) {
  1760       Node nca = INVALID;
  1761       std::vector<Node> left_path, right_path;
  1762 
  1763       {
  1764         std::set<Node> left_set, right_set;
  1765         Node left = _graph.u(edge);
  1766         left_path.push_back(left);
  1767         left_set.insert(left);
  1768 
  1769         Node right = _graph.v(edge);
  1770         right_path.push_back(right);
  1771         right_set.insert(right);
  1772 
  1773         while (true) {
  1774 
  1775           if (left_set.find(right) != left_set.end()) {
  1776             nca = right;
  1777             break;
  1778           }
  1779 
  1780           if ((*_matching)[left] == INVALID) break;
  1781 
  1782           left = _graph.target((*_matching)[left]);
  1783           left_path.push_back(left);
  1784           left = _graph.target((*_pred)[left]);
  1785           left_path.push_back(left);
  1786 
  1787           left_set.insert(left);
  1788 
  1789           if (right_set.find(left) != right_set.end()) {
  1790             nca = left;
  1791             break;
  1792           }
  1793 
  1794           if ((*_matching)[right] == INVALID) break;
  1795 
  1796           right = _graph.target((*_matching)[right]);
  1797           right_path.push_back(right);
  1798           right = _graph.target((*_pred)[right]);
  1799           right_path.push_back(right);
  1800 
  1801           right_set.insert(right);
  1802 
  1803         }
  1804 
  1805         if (nca == INVALID) {
  1806           if ((*_matching)[left] == INVALID) {
  1807             nca = right;
  1808             while (left_set.find(nca) == left_set.end()) {
  1809               nca = _graph.target((*_matching)[nca]);
  1810               right_path.push_back(nca);
  1811               nca = _graph.target((*_pred)[nca]);
  1812               right_path.push_back(nca);
  1813             }
  1814           } else {
  1815             nca = left;
  1816             while (right_set.find(nca) == right_set.end()) {
  1817               nca = _graph.target((*_matching)[nca]);
  1818               left_path.push_back(nca);
  1819               nca = _graph.target((*_pred)[nca]);
  1820               left_path.push_back(nca);
  1821             }
  1822           }
  1823         }
  1824       }
  1825 
  1826       alternatePath(nca, tree);
  1827       Arc prev;
  1828 
  1829       prev = _graph.direct(edge, true);
  1830       for (int i = 0; left_path[i] != nca; i += 2) {
  1831         _matching->set(left_path[i], prev);
  1832         _status->set(left_path[i], MATCHED);
  1833         evenToMatched(left_path[i], tree);
  1834 
  1835         prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
  1836         _status->set(left_path[i + 1], MATCHED);
  1837         oddToMatched(left_path[i + 1]);
  1838       }
  1839       _matching->set(nca, prev);
  1840 
  1841       for (int i = 0; right_path[i] != nca; i += 2) {
  1842         _status->set(right_path[i], MATCHED);
  1843         evenToMatched(right_path[i], tree);
  1844 
  1845         _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
  1846         _status->set(right_path[i + 1], MATCHED);
  1847         oddToMatched(right_path[i + 1]);
  1848       }
  1849 
  1850       destroyTree(tree);
  1851     }
  1852 
  1853     void extractCycle(const Arc &arc) {
  1854       Node left = _graph.source(arc);
  1855       Node odd = _graph.target((*_matching)[left]);
  1856       Arc prev;
  1857       while (odd != left) {
  1858         Node even = _graph.target((*_matching)[odd]);
  1859         prev = (*_matching)[odd];
  1860         odd = _graph.target((*_matching)[even]);
  1861         _matching->set(even, _graph.oppositeArc(prev));
  1862       }
  1863       _matching->set(left, arc);
  1864 
  1865       Node right = _graph.target(arc);
  1866       int right_tree = _tree_set->find(right);
  1867       alternatePath(right, right_tree);
  1868       destroyTree(right_tree);
  1869       _matching->set(right, _graph.oppositeArc(arc));
  1870     }
  1871 
  1872   public:
  1873 
  1874     /// \brief Constructor
  1875     ///
  1876     /// Constructor.
  1877     MaxWeightedPerfectFractionalMatching(const Graph& graph,
  1878                                          const WeightMap& weight,
  1879                                          bool allow_loops = true)
  1880       : _graph(graph), _weight(weight), _matching(0),
  1881       _node_potential(0), _node_num(0), _allow_loops(allow_loops),
  1882       _status(0),  _pred(0),
  1883       _tree_set_index(0), _tree_set(0),
  1884 
  1885       _delta2_index(0), _delta2(0),
  1886       _delta3_index(0), _delta3(0),
  1887 
  1888       _delta_sum() {}
  1889 
  1890     ~MaxWeightedPerfectFractionalMatching() {
  1891       destroyStructures();
  1892     }
  1893 
  1894     /// \name Execution Control
  1895     /// The simplest way to execute the algorithm is to use the
  1896     /// \ref run() member function.
  1897 
  1898     ///@{
  1899 
  1900     /// \brief Initialize the algorithm
  1901     ///
  1902     /// This function initializes the algorithm.
  1903     void init() {
  1904       createStructures();
  1905 
  1906       for (NodeIt n(_graph); n != INVALID; ++n) {
  1907         (*_delta2_index)[n] = _delta2->PRE_HEAP;
  1908       }
  1909       for (EdgeIt e(_graph); e != INVALID; ++e) {
  1910         (*_delta3_index)[e] = _delta3->PRE_HEAP;
  1911       }
  1912 
  1913       for (NodeIt n(_graph); n != INVALID; ++n) {
  1914         Value max = - std::numeric_limits<Value>::max();
  1915         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  1916           if (_graph.target(e) == n && !_allow_loops) continue;
  1917           if ((dualScale * _weight[e]) / 2 > max) {
  1918             max = (dualScale * _weight[e]) / 2;
  1919           }
  1920         }
  1921         _node_potential->set(n, max);
  1922 
  1923         _tree_set->insert(n);
  1924 
  1925         _matching->set(n, INVALID);
  1926         _status->set(n, EVEN);
  1927       }
  1928 
  1929       for (EdgeIt e(_graph); e != INVALID; ++e) {
  1930         Node left = _graph.u(e);
  1931         Node right = _graph.v(e);
  1932         if (left == right && !_allow_loops) continue;
  1933         _delta3->push(e, ((*_node_potential)[left] +
  1934                           (*_node_potential)[right] -
  1935                           dualScale * _weight[e]) / 2);
  1936       }
  1937     }
  1938 
  1939     /// \brief Start the algorithm
  1940     ///
  1941     /// This function starts the algorithm.
  1942     ///
  1943     /// \pre \ref init() must be called before using this function.
  1944     bool start() {
  1945       enum OpType {
  1946         D2, D3
  1947       };
  1948 
  1949       int unmatched = _node_num;
  1950       while (unmatched > 0) {
  1951         Value d2 = !_delta2->empty() ?
  1952           _delta2->prio() : std::numeric_limits<Value>::max();
  1953 
  1954         Value d3 = !_delta3->empty() ?
  1955           _delta3->prio() : std::numeric_limits<Value>::max();
  1956 
  1957         _delta_sum = d3; OpType ot = D3;
  1958         if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
  1959 
  1960         if (_delta_sum == std::numeric_limits<Value>::max()) {
  1961           return false;
  1962         }
  1963 
  1964         switch (ot) {
  1965         case D2:
  1966           {
  1967             Node n = _delta2->top();
  1968             Arc a = (*_pred)[n];
  1969             if ((*_matching)[n] == INVALID) {
  1970               augmentOnArc(a);
  1971               --unmatched;
  1972             } else {
  1973               Node v = _graph.target((*_matching)[n]);
  1974               if ((*_matching)[n] !=
  1975                   _graph.oppositeArc((*_matching)[v])) {
  1976                 extractCycle(a);
  1977                 --unmatched;
  1978               } else {
  1979                 extendOnArc(a);
  1980               }
  1981             }
  1982           } break;
  1983         case D3:
  1984           {
  1985             Edge e = _delta3->top();
  1986 
  1987             Node left = _graph.u(e);
  1988             Node right = _graph.v(e);
  1989 
  1990             int left_tree = _tree_set->find(left);
  1991             int right_tree = _tree_set->find(right);
  1992 
  1993             if (left_tree == right_tree) {
  1994               cycleOnEdge(e, left_tree);
  1995               --unmatched;
  1996             } else {
  1997               augmentOnEdge(e);
  1998               unmatched -= 2;
  1999             }
  2000           } break;
  2001         }
  2002       }
  2003       return true;
  2004     }
  2005 
  2006     /// \brief Run the algorithm.
  2007     ///
  2008     /// This method runs the \c %MaxWeightedMatching algorithm.
  2009     ///
  2010     /// \note mwfm.run() is just a shortcut of the following code.
  2011     /// \code
  2012     ///   mwpfm.init();
  2013     ///   mwpfm.start();
  2014     /// \endcode
  2015     bool run() {
  2016       init();
  2017       return start();
  2018     }
  2019 
  2020     /// @}
  2021 
  2022     /// \name Primal Solution
  2023     /// Functions to get the primal solution, i.e. the maximum weighted
  2024     /// matching.\n
  2025     /// Either \ref run() or \ref start() function should be called before
  2026     /// using them.
  2027 
  2028     /// @{
  2029 
  2030     /// \brief Return the weight of the matching.
  2031     ///
  2032     /// This function returns the weight of the found matching. This
  2033     /// value is scaled by \ref primalScale "primal scale".
  2034     ///
  2035     /// \pre Either run() or start() must be called before using this function.
  2036     Value matchingWeight() const {
  2037       Value sum = 0;
  2038       for (NodeIt n(_graph); n != INVALID; ++n) {
  2039         if ((*_matching)[n] != INVALID) {
  2040           sum += _weight[(*_matching)[n]];
  2041         }
  2042       }
  2043       return sum * primalScale / 2;
  2044     }
  2045 
  2046     /// \brief Return the number of covered nodes in the matching.
  2047     ///
  2048     /// This function returns the number of covered nodes in the matching.
  2049     ///
  2050     /// \pre Either run() or start() must be called before using this function.
  2051     int matchingSize() const {
  2052       int num = 0;
  2053       for (NodeIt n(_graph); n != INVALID; ++n) {
  2054         if ((*_matching)[n] != INVALID) {
  2055           ++num;
  2056         }
  2057       }
  2058       return num;
  2059     }
  2060 
  2061     /// \brief Return \c true if the given edge is in the matching.
  2062     ///
  2063     /// This function returns \c true if the given edge is in the
  2064     /// found matching. The result is scaled by \ref primalScale
  2065     /// "primal scale".
  2066     ///
  2067     /// \pre Either run() or start() must be called before using this function.
  2068     Value matching(const Edge& edge) const {
  2069       return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
  2070         * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0)
  2071         * primalScale / 2;
  2072     }
  2073 
  2074     /// \brief Return the fractional matching arc (or edge) incident
  2075     /// to the given node.
  2076     ///
  2077     /// This function returns one of the fractional matching arc (or
  2078     /// edge) incident to the given node in the found matching or \c
  2079     /// INVALID if the node is not covered by the matching or if the
  2080     /// node is on an odd length cycle then it is the successor edge
  2081     /// on the cycle.
  2082     ///
  2083     /// \pre Either run() or start() must be called before using this function.
  2084     Arc matching(const Node& node) const {
  2085       return (*_matching)[node];
  2086     }
  2087 
  2088     /// \brief Return a const reference to the matching map.
  2089     ///
  2090     /// This function returns a const reference to a node map that stores
  2091     /// the matching arc (or edge) incident to each node.
  2092     const MatchingMap& matchingMap() const {
  2093       return *_matching;
  2094     }
  2095 
  2096     /// @}
  2097 
  2098     /// \name Dual Solution
  2099     /// Functions to get the dual solution.\n
  2100     /// Either \ref run() or \ref start() function should be called before
  2101     /// using them.
  2102 
  2103     /// @{
  2104 
  2105     /// \brief Return the value of the dual solution.
  2106     ///
  2107     /// This function returns the value of the dual solution.
  2108     /// It should be equal to the primal value scaled by \ref dualScale
  2109     /// "dual scale".
  2110     ///
  2111     /// \pre Either run() or start() must be called before using this function.
  2112     Value dualValue() const {
  2113       Value sum = 0;
  2114       for (NodeIt n(_graph); n != INVALID; ++n) {
  2115         sum += nodeValue(n);
  2116       }
  2117       return sum;
  2118     }
  2119 
  2120     /// \brief Return the dual value (potential) of the given node.
  2121     ///
  2122     /// This function returns the dual value (potential) of the given node.
  2123     ///
  2124     /// \pre Either run() or start() must be called before using this function.
  2125     Value nodeValue(const Node& n) const {
  2126       return (*_node_potential)[n];
  2127     }
  2128 
  2129     /// @}
  2130 
  2131   };
  2132 
  2133 } //END OF NAMESPACE LEMON
  2134 
  2135 #endif //LEMON_FRACTIONAL_MATCHING_H