lemon/fractional_matching.h
author Balazs Dezso <deba@inf.elte.hu>
Thu, 04 Mar 2010 10:17:02 +0100
changeset 950 86613aa28a0c
parent 948 636dadefe1e6
child 951 41d7ac528c3a
permissions -rw-r--r--
Fix documentation issues (#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 MaxFractionalMatching::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 uses priority queues and
   636   /// provides \f$O(nm\log n)\f$ time complexity.
   637   ///
   638   /// The maximum weighted fractional matching is a relaxation of the
   639   /// maximum weighted matching problem where the odd set constraints
   640   /// are omitted.
   641   /// It can be formulated with the following linear program.
   642   /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
   643   /// \f[x_e \ge 0\quad \forall e\in E\f]
   644   /// \f[\max \sum_{e\in E}x_ew_e\f]
   645   /// where \f$\delta(X)\f$ is the set of edges incident to a node in
   646   /// \f$X\f$. The result must be the union of a matching with one
   647   /// value edges and a set of odd length cycles with half value edges.
   648   ///
   649   /// The algorithm calculates an optimal fractional matching and a
   650   /// proof of the optimality. The solution of the dual problem can be
   651   /// used to check the result of the algorithm. The dual linear
   652   /// problem is the following.
   653   /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
   654   /// \f[y_u \ge 0 \quad \forall u \in V\f]
   655   /// \f[\min \sum_{u \in V}y_u \f]
   656   ///
   657   /// The algorithm can be executed with the run() function.
   658   /// After it the matching (the primal solution) and the dual solution
   659   /// can be obtained using the query functions.
   660   ///
   661   /// If the value type is integer, then the primal and the dual
   662   /// solutions are multiplied by
   663   /// \ref MaxWeightedFractionalMatching::primalScale "2" and
   664   /// \ref MaxWeightedFractionalMatching::dualScale "4" respectively.
   665   ///
   666   /// \tparam GR The undirected graph type the algorithm runs on.
   667   /// \tparam WM The type edge weight map. The default type is
   668   /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
   669 #ifdef DOXYGEN
   670   template <typename GR, typename WM>
   671 #else
   672   template <typename GR,
   673             typename WM = typename GR::template EdgeMap<int> >
   674 #endif
   675   class MaxWeightedFractionalMatching {
   676   public:
   677 
   678     /// The graph type of the algorithm
   679     typedef GR Graph;
   680     /// The type of the edge weight map
   681     typedef WM WeightMap;
   682     /// The value type of the edge weights
   683     typedef typename WeightMap::Value Value;
   684 
   685     /// The type of the matching map
   686     typedef typename Graph::template NodeMap<typename Graph::Arc>
   687     MatchingMap;
   688 
   689     /// \brief Scaling factor for primal solution
   690     ///
   691     /// Scaling factor for primal solution. It is equal to 2 or 1
   692     /// according to the value type.
   693     static const int primalScale =
   694       std::numeric_limits<Value>::is_integer ? 2 : 1;
   695 
   696     /// \brief Scaling factor for dual solution
   697     ///
   698     /// Scaling factor for dual solution. It is equal to 4 or 1
   699     /// according to the value type.
   700     static const int dualScale =
   701       std::numeric_limits<Value>::is_integer ? 4 : 1;
   702 
   703   private:
   704 
   705     TEMPLATE_GRAPH_TYPEDEFS(Graph);
   706 
   707     typedef typename Graph::template NodeMap<Value> NodePotential;
   708 
   709     const Graph& _graph;
   710     const WeightMap& _weight;
   711 
   712     MatchingMap* _matching;
   713     NodePotential* _node_potential;
   714 
   715     int _node_num;
   716     bool _allow_loops;
   717 
   718     enum Status {
   719       EVEN = -1, MATCHED = 0, ODD = 1
   720     };
   721 
   722     typedef typename Graph::template NodeMap<Status> StatusMap;
   723     StatusMap* _status;
   724 
   725     typedef typename Graph::template NodeMap<Arc> PredMap;
   726     PredMap* _pred;
   727 
   728     typedef ExtendFindEnum<IntNodeMap> TreeSet;
   729 
   730     IntNodeMap *_tree_set_index;
   731     TreeSet *_tree_set;
   732 
   733     IntNodeMap *_delta1_index;
   734     BinHeap<Value, IntNodeMap> *_delta1;
   735 
   736     IntNodeMap *_delta2_index;
   737     BinHeap<Value, IntNodeMap> *_delta2;
   738 
   739     IntEdgeMap *_delta3_index;
   740     BinHeap<Value, IntEdgeMap> *_delta3;
   741 
   742     Value _delta_sum;
   743 
   744     void createStructures() {
   745       _node_num = countNodes(_graph);
   746 
   747       if (!_matching) {
   748         _matching = new MatchingMap(_graph);
   749       }
   750       if (!_node_potential) {
   751         _node_potential = new NodePotential(_graph);
   752       }
   753       if (!_status) {
   754         _status = new StatusMap(_graph);
   755       }
   756       if (!_pred) {
   757         _pred = new PredMap(_graph);
   758       }
   759       if (!_tree_set) {
   760         _tree_set_index = new IntNodeMap(_graph);
   761         _tree_set = new TreeSet(*_tree_set_index);
   762       }
   763       if (!_delta1) {
   764         _delta1_index = new IntNodeMap(_graph);
   765         _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
   766       }
   767       if (!_delta2) {
   768         _delta2_index = new IntNodeMap(_graph);
   769         _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
   770       }
   771       if (!_delta3) {
   772         _delta3_index = new IntEdgeMap(_graph);
   773         _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
   774       }
   775     }
   776 
   777     void destroyStructures() {
   778       if (_matching) {
   779         delete _matching;
   780       }
   781       if (_node_potential) {
   782         delete _node_potential;
   783       }
   784       if (_status) {
   785         delete _status;
   786       }
   787       if (_pred) {
   788         delete _pred;
   789       }
   790       if (_tree_set) {
   791         delete _tree_set_index;
   792         delete _tree_set;
   793       }
   794       if (_delta1) {
   795         delete _delta1_index;
   796         delete _delta1;
   797       }
   798       if (_delta2) {
   799         delete _delta2_index;
   800         delete _delta2;
   801       }
   802       if (_delta3) {
   803         delete _delta3_index;
   804         delete _delta3;
   805       }
   806     }
   807 
   808     void matchedToEven(Node node, int tree) {
   809       _tree_set->insert(node, tree);
   810       _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
   811       _delta1->push(node, (*_node_potential)[node]);
   812 
   813       if (_delta2->state(node) == _delta2->IN_HEAP) {
   814         _delta2->erase(node);
   815       }
   816 
   817       for (InArcIt a(_graph, node); a != INVALID; ++a) {
   818         Node v = _graph.source(a);
   819         Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
   820           dualScale * _weight[a];
   821         if (node == v) {
   822           if (_allow_loops && _graph.direction(a)) {
   823             _delta3->push(a, rw / 2);
   824           }
   825         } else if ((*_status)[v] == EVEN) {
   826           _delta3->push(a, rw / 2);
   827         } else if ((*_status)[v] == MATCHED) {
   828           if (_delta2->state(v) != _delta2->IN_HEAP) {
   829             _pred->set(v, a);
   830             _delta2->push(v, rw);
   831           } else if ((*_delta2)[v] > rw) {
   832             _pred->set(v, a);
   833             _delta2->decrease(v, rw);
   834           }
   835         }
   836       }
   837     }
   838 
   839     void matchedToOdd(Node node, int tree) {
   840       _tree_set->insert(node, tree);
   841       _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
   842 
   843       if (_delta2->state(node) == _delta2->IN_HEAP) {
   844         _delta2->erase(node);
   845       }
   846     }
   847 
   848     void evenToMatched(Node node, int tree) {
   849       _delta1->erase(node);
   850       _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
   851       Arc min = INVALID;
   852       Value minrw = std::numeric_limits<Value>::max();
   853       for (InArcIt a(_graph, node); a != INVALID; ++a) {
   854         Node v = _graph.source(a);
   855         Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
   856           dualScale * _weight[a];
   857 
   858         if (node == v) {
   859           if (_allow_loops && _graph.direction(a)) {
   860             _delta3->erase(a);
   861           }
   862         } else if ((*_status)[v] == EVEN) {
   863           _delta3->erase(a);
   864           if (minrw > rw) {
   865             min = _graph.oppositeArc(a);
   866             minrw = rw;
   867           }
   868         } else if ((*_status)[v]  == MATCHED) {
   869           if ((*_pred)[v] == a) {
   870             Arc mina = INVALID;
   871             Value minrwa = std::numeric_limits<Value>::max();
   872             for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
   873               Node va = _graph.target(aa);
   874               if ((*_status)[va] != EVEN ||
   875                   _tree_set->find(va) == tree) continue;
   876               Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
   877                 dualScale * _weight[aa];
   878               if (minrwa > rwa) {
   879                 minrwa = rwa;
   880                 mina = aa;
   881               }
   882             }
   883             if (mina != INVALID) {
   884               _pred->set(v, mina);
   885               _delta2->increase(v, minrwa);
   886             } else {
   887               _pred->set(v, INVALID);
   888               _delta2->erase(v);
   889             }
   890           }
   891         }
   892       }
   893       if (min != INVALID) {
   894         _pred->set(node, min);
   895         _delta2->push(node, minrw);
   896       } else {
   897         _pred->set(node, INVALID);
   898       }
   899     }
   900 
   901     void oddToMatched(Node node) {
   902       _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
   903       Arc min = INVALID;
   904       Value minrw = std::numeric_limits<Value>::max();
   905       for (InArcIt a(_graph, node); a != INVALID; ++a) {
   906         Node v = _graph.source(a);
   907         if ((*_status)[v] != EVEN) continue;
   908         Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
   909           dualScale * _weight[a];
   910 
   911         if (minrw > rw) {
   912           min = _graph.oppositeArc(a);
   913           minrw = rw;
   914         }
   915       }
   916       if (min != INVALID) {
   917         _pred->set(node, min);
   918         _delta2->push(node, minrw);
   919       } else {
   920         _pred->set(node, INVALID);
   921       }
   922     }
   923 
   924     void alternatePath(Node even, int tree) {
   925       Node odd;
   926 
   927       _status->set(even, MATCHED);
   928       evenToMatched(even, tree);
   929 
   930       Arc prev = (*_matching)[even];
   931       while (prev != INVALID) {
   932         odd = _graph.target(prev);
   933         even = _graph.target((*_pred)[odd]);
   934         _matching->set(odd, (*_pred)[odd]);
   935         _status->set(odd, MATCHED);
   936         oddToMatched(odd);
   937 
   938         prev = (*_matching)[even];
   939         _status->set(even, MATCHED);
   940         _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
   941         evenToMatched(even, tree);
   942       }
   943     }
   944 
   945     void destroyTree(int tree) {
   946       for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
   947         if ((*_status)[n] == EVEN) {
   948           _status->set(n, MATCHED);
   949           evenToMatched(n, tree);
   950         } else if ((*_status)[n] == ODD) {
   951           _status->set(n, MATCHED);
   952           oddToMatched(n);
   953         }
   954       }
   955       _tree_set->eraseClass(tree);
   956     }
   957 
   958 
   959     void unmatchNode(const Node& node) {
   960       int tree = _tree_set->find(node);
   961 
   962       alternatePath(node, tree);
   963       destroyTree(tree);
   964 
   965       _matching->set(node, INVALID);
   966     }
   967 
   968 
   969     void augmentOnEdge(const Edge& edge) {
   970       Node left = _graph.u(edge);
   971       int left_tree = _tree_set->find(left);
   972 
   973       alternatePath(left, left_tree);
   974       destroyTree(left_tree);
   975       _matching->set(left, _graph.direct(edge, true));
   976 
   977       Node right = _graph.v(edge);
   978       int right_tree = _tree_set->find(right);
   979 
   980       alternatePath(right, right_tree);
   981       destroyTree(right_tree);
   982       _matching->set(right, _graph.direct(edge, false));
   983     }
   984 
   985     void augmentOnArc(const Arc& arc) {
   986       Node left = _graph.source(arc);
   987       _status->set(left, MATCHED);
   988       _matching->set(left, arc);
   989       _pred->set(left, arc);
   990 
   991       Node right = _graph.target(arc);
   992       int right_tree = _tree_set->find(right);
   993 
   994       alternatePath(right, right_tree);
   995       destroyTree(right_tree);
   996       _matching->set(right, _graph.oppositeArc(arc));
   997     }
   998 
   999     void extendOnArc(const Arc& arc) {
  1000       Node base = _graph.target(arc);
  1001       int tree = _tree_set->find(base);
  1002 
  1003       Node odd = _graph.source(arc);
  1004       _tree_set->insert(odd, tree);
  1005       _status->set(odd, ODD);
  1006       matchedToOdd(odd, tree);
  1007       _pred->set(odd, arc);
  1008 
  1009       Node even = _graph.target((*_matching)[odd]);
  1010       _tree_set->insert(even, tree);
  1011       _status->set(even, EVEN);
  1012       matchedToEven(even, tree);
  1013     }
  1014 
  1015     void cycleOnEdge(const Edge& edge, int tree) {
  1016       Node nca = INVALID;
  1017       std::vector<Node> left_path, right_path;
  1018 
  1019       {
  1020         std::set<Node> left_set, right_set;
  1021         Node left = _graph.u(edge);
  1022         left_path.push_back(left);
  1023         left_set.insert(left);
  1024 
  1025         Node right = _graph.v(edge);
  1026         right_path.push_back(right);
  1027         right_set.insert(right);
  1028 
  1029         while (true) {
  1030 
  1031           if (left_set.find(right) != left_set.end()) {
  1032             nca = right;
  1033             break;
  1034           }
  1035 
  1036           if ((*_matching)[left] == INVALID) break;
  1037 
  1038           left = _graph.target((*_matching)[left]);
  1039           left_path.push_back(left);
  1040           left = _graph.target((*_pred)[left]);
  1041           left_path.push_back(left);
  1042 
  1043           left_set.insert(left);
  1044 
  1045           if (right_set.find(left) != right_set.end()) {
  1046             nca = left;
  1047             break;
  1048           }
  1049 
  1050           if ((*_matching)[right] == INVALID) break;
  1051 
  1052           right = _graph.target((*_matching)[right]);
  1053           right_path.push_back(right);
  1054           right = _graph.target((*_pred)[right]);
  1055           right_path.push_back(right);
  1056 
  1057           right_set.insert(right);
  1058 
  1059         }
  1060 
  1061         if (nca == INVALID) {
  1062           if ((*_matching)[left] == INVALID) {
  1063             nca = right;
  1064             while (left_set.find(nca) == left_set.end()) {
  1065               nca = _graph.target((*_matching)[nca]);
  1066               right_path.push_back(nca);
  1067               nca = _graph.target((*_pred)[nca]);
  1068               right_path.push_back(nca);
  1069             }
  1070           } else {
  1071             nca = left;
  1072             while (right_set.find(nca) == right_set.end()) {
  1073               nca = _graph.target((*_matching)[nca]);
  1074               left_path.push_back(nca);
  1075               nca = _graph.target((*_pred)[nca]);
  1076               left_path.push_back(nca);
  1077             }
  1078           }
  1079         }
  1080       }
  1081 
  1082       alternatePath(nca, tree);
  1083       Arc prev;
  1084 
  1085       prev = _graph.direct(edge, true);
  1086       for (int i = 0; left_path[i] != nca; i += 2) {
  1087         _matching->set(left_path[i], prev);
  1088         _status->set(left_path[i], MATCHED);
  1089         evenToMatched(left_path[i], tree);
  1090 
  1091         prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
  1092         _status->set(left_path[i + 1], MATCHED);
  1093         oddToMatched(left_path[i + 1]);
  1094       }
  1095       _matching->set(nca, prev);
  1096 
  1097       for (int i = 0; right_path[i] != nca; i += 2) {
  1098         _status->set(right_path[i], MATCHED);
  1099         evenToMatched(right_path[i], tree);
  1100 
  1101         _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
  1102         _status->set(right_path[i + 1], MATCHED);
  1103         oddToMatched(right_path[i + 1]);
  1104       }
  1105 
  1106       destroyTree(tree);
  1107     }
  1108 
  1109     void extractCycle(const Arc &arc) {
  1110       Node left = _graph.source(arc);
  1111       Node odd = _graph.target((*_matching)[left]);
  1112       Arc prev;
  1113       while (odd != left) {
  1114         Node even = _graph.target((*_matching)[odd]);
  1115         prev = (*_matching)[odd];
  1116         odd = _graph.target((*_matching)[even]);
  1117         _matching->set(even, _graph.oppositeArc(prev));
  1118       }
  1119       _matching->set(left, arc);
  1120 
  1121       Node right = _graph.target(arc);
  1122       int right_tree = _tree_set->find(right);
  1123       alternatePath(right, right_tree);
  1124       destroyTree(right_tree);
  1125       _matching->set(right, _graph.oppositeArc(arc));
  1126     }
  1127 
  1128   public:
  1129 
  1130     /// \brief Constructor
  1131     ///
  1132     /// Constructor.
  1133     MaxWeightedFractionalMatching(const Graph& graph, const WeightMap& weight,
  1134                                   bool allow_loops = true)
  1135       : _graph(graph), _weight(weight), _matching(0),
  1136       _node_potential(0), _node_num(0), _allow_loops(allow_loops),
  1137       _status(0),  _pred(0),
  1138       _tree_set_index(0), _tree_set(0),
  1139 
  1140       _delta1_index(0), _delta1(0),
  1141       _delta2_index(0), _delta2(0),
  1142       _delta3_index(0), _delta3(0),
  1143 
  1144       _delta_sum() {}
  1145 
  1146     ~MaxWeightedFractionalMatching() {
  1147       destroyStructures();
  1148     }
  1149 
  1150     /// \name Execution Control
  1151     /// The simplest way to execute the algorithm is to use the
  1152     /// \ref run() member function.
  1153 
  1154     ///@{
  1155 
  1156     /// \brief Initialize the algorithm
  1157     ///
  1158     /// This function initializes the algorithm.
  1159     void init() {
  1160       createStructures();
  1161 
  1162       for (NodeIt n(_graph); n != INVALID; ++n) {
  1163         (*_delta1_index)[n] = _delta1->PRE_HEAP;
  1164         (*_delta2_index)[n] = _delta2->PRE_HEAP;
  1165       }
  1166       for (EdgeIt e(_graph); e != INVALID; ++e) {
  1167         (*_delta3_index)[e] = _delta3->PRE_HEAP;
  1168       }
  1169 
  1170       for (NodeIt n(_graph); n != INVALID; ++n) {
  1171         Value max = 0;
  1172         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  1173           if (_graph.target(e) == n && !_allow_loops) continue;
  1174           if ((dualScale * _weight[e]) / 2 > max) {
  1175             max = (dualScale * _weight[e]) / 2;
  1176           }
  1177         }
  1178         _node_potential->set(n, max);
  1179         _delta1->push(n, max);
  1180 
  1181         _tree_set->insert(n);
  1182 
  1183         _matching->set(n, INVALID);
  1184         _status->set(n, EVEN);
  1185       }
  1186 
  1187       for (EdgeIt e(_graph); e != INVALID; ++e) {
  1188         Node left = _graph.u(e);
  1189         Node right = _graph.v(e);
  1190         if (left == right && !_allow_loops) continue;
  1191         _delta3->push(e, ((*_node_potential)[left] +
  1192                           (*_node_potential)[right] -
  1193                           dualScale * _weight[e]) / 2);
  1194       }
  1195     }
  1196 
  1197     /// \brief Start the algorithm
  1198     ///
  1199     /// This function starts the algorithm.
  1200     ///
  1201     /// \pre \ref init() must be called before using this function.
  1202     void start() {
  1203       enum OpType {
  1204         D1, D2, D3
  1205       };
  1206 
  1207       int unmatched = _node_num;
  1208       while (unmatched > 0) {
  1209         Value d1 = !_delta1->empty() ?
  1210           _delta1->prio() : std::numeric_limits<Value>::max();
  1211 
  1212         Value d2 = !_delta2->empty() ?
  1213           _delta2->prio() : std::numeric_limits<Value>::max();
  1214 
  1215         Value d3 = !_delta3->empty() ?
  1216           _delta3->prio() : std::numeric_limits<Value>::max();
  1217 
  1218         _delta_sum = d3; OpType ot = D3;
  1219         if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
  1220         if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
  1221 
  1222         switch (ot) {
  1223         case D1:
  1224           {
  1225             Node n = _delta1->top();
  1226             unmatchNode(n);
  1227             --unmatched;
  1228           }
  1229           break;
  1230         case D2:
  1231           {
  1232             Node n = _delta2->top();
  1233             Arc a = (*_pred)[n];
  1234             if ((*_matching)[n] == INVALID) {
  1235               augmentOnArc(a);
  1236               --unmatched;
  1237             } else {
  1238               Node v = _graph.target((*_matching)[n]);
  1239               if ((*_matching)[n] !=
  1240                   _graph.oppositeArc((*_matching)[v])) {
  1241                 extractCycle(a);
  1242                 --unmatched;
  1243               } else {
  1244                 extendOnArc(a);
  1245               }
  1246             }
  1247           } break;
  1248         case D3:
  1249           {
  1250             Edge e = _delta3->top();
  1251 
  1252             Node left = _graph.u(e);
  1253             Node right = _graph.v(e);
  1254 
  1255             int left_tree = _tree_set->find(left);
  1256             int right_tree = _tree_set->find(right);
  1257 
  1258             if (left_tree == right_tree) {
  1259               cycleOnEdge(e, left_tree);
  1260               --unmatched;
  1261             } else {
  1262               augmentOnEdge(e);
  1263               unmatched -= 2;
  1264             }
  1265           } break;
  1266         }
  1267       }
  1268     }
  1269 
  1270     /// \brief Run the algorithm.
  1271     ///
  1272     /// This method runs the \c %MaxWeightedFractionalMatching algorithm.
  1273     ///
  1274     /// \note mwfm.run() is just a shortcut of the following code.
  1275     /// \code
  1276     ///   mwfm.init();
  1277     ///   mwfm.start();
  1278     /// \endcode
  1279     void run() {
  1280       init();
  1281       start();
  1282     }
  1283 
  1284     /// @}
  1285 
  1286     /// \name Primal Solution
  1287     /// Functions to get the primal solution, i.e. the maximum weighted
  1288     /// matching.\n
  1289     /// Either \ref run() or \ref start() function should be called before
  1290     /// using them.
  1291 
  1292     /// @{
  1293 
  1294     /// \brief Return the weight of the matching.
  1295     ///
  1296     /// This function returns the weight of the found matching. This
  1297     /// value is scaled by \ref primalScale "primal scale".
  1298     ///
  1299     /// \pre Either run() or start() must be called before using this function.
  1300     Value matchingWeight() const {
  1301       Value sum = 0;
  1302       for (NodeIt n(_graph); n != INVALID; ++n) {
  1303         if ((*_matching)[n] != INVALID) {
  1304           sum += _weight[(*_matching)[n]];
  1305         }
  1306       }
  1307       return sum * primalScale / 2;
  1308     }
  1309 
  1310     /// \brief Return the number of covered nodes in the matching.
  1311     ///
  1312     /// This function returns the number of covered nodes in the matching.
  1313     ///
  1314     /// \pre Either run() or start() must be called before using this function.
  1315     int matchingSize() const {
  1316       int num = 0;
  1317       for (NodeIt n(_graph); n != INVALID; ++n) {
  1318         if ((*_matching)[n] != INVALID) {
  1319           ++num;
  1320         }
  1321       }
  1322       return num;
  1323     }
  1324 
  1325     /// \brief Return \c true if the given edge is in the matching.
  1326     ///
  1327     /// This function returns \c true if the given edge is in the
  1328     /// found matching. The result is scaled by \ref primalScale
  1329     /// "primal scale".
  1330     ///
  1331     /// \pre Either run() or start() must be called before using this function.
  1332     Value matching(const Edge& edge) const {
  1333       return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
  1334         * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0)
  1335         * primalScale / 2;
  1336     }
  1337 
  1338     /// \brief Return the fractional matching arc (or edge) incident
  1339     /// to the given node.
  1340     ///
  1341     /// This function returns one of the fractional matching arc (or
  1342     /// edge) incident to the given node in the found matching or \c
  1343     /// INVALID if the node is not covered by the matching or if the
  1344     /// node is on an odd length cycle then it is the successor edge
  1345     /// on the cycle.
  1346     ///
  1347     /// \pre Either run() or start() must be called before using this function.
  1348     Arc matching(const Node& node) const {
  1349       return (*_matching)[node];
  1350     }
  1351 
  1352     /// \brief Return a const reference to the matching map.
  1353     ///
  1354     /// This function returns a const reference to a node map that stores
  1355     /// the matching arc (or edge) incident to each node.
  1356     const MatchingMap& matchingMap() const {
  1357       return *_matching;
  1358     }
  1359 
  1360     /// @}
  1361 
  1362     /// \name Dual Solution
  1363     /// Functions to get the dual solution.\n
  1364     /// Either \ref run() or \ref start() function should be called before
  1365     /// using them.
  1366 
  1367     /// @{
  1368 
  1369     /// \brief Return the value of the dual solution.
  1370     ///
  1371     /// This function returns the value of the dual solution.
  1372     /// It should be equal to the primal value scaled by \ref dualScale
  1373     /// "dual scale".
  1374     ///
  1375     /// \pre Either run() or start() must be called before using this function.
  1376     Value dualValue() const {
  1377       Value sum = 0;
  1378       for (NodeIt n(_graph); n != INVALID; ++n) {
  1379         sum += nodeValue(n);
  1380       }
  1381       return sum;
  1382     }
  1383 
  1384     /// \brief Return the dual value (potential) of the given node.
  1385     ///
  1386     /// This function returns the dual value (potential) of the given node.
  1387     ///
  1388     /// \pre Either run() or start() must be called before using this function.
  1389     Value nodeValue(const Node& n) const {
  1390       return (*_node_potential)[n];
  1391     }
  1392 
  1393     /// @}
  1394 
  1395   };
  1396 
  1397   /// \ingroup matching
  1398   ///
  1399   /// \brief Weighted fractional perfect matching in general graphs
  1400   ///
  1401   /// This class provides an efficient implementation of fractional
  1402   /// matching algorithm. The implementation uses priority queues and
  1403   /// provides \f$O(nm\log n)\f$ time complexity.
  1404   ///
  1405   /// The maximum weighted fractional perfect matching is a relaxation
  1406   /// of the maximum weighted perfect matching problem where the odd
  1407   /// set constraints are omitted.
  1408   /// It can be formulated with the following linear program.
  1409   /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
  1410   /// \f[x_e \ge 0\quad \forall e\in E\f]
  1411   /// \f[\max \sum_{e\in E}x_ew_e\f]
  1412   /// where \f$\delta(X)\f$ is the set of edges incident to a node in
  1413   /// \f$X\f$. The result must be the union of a matching with one
  1414   /// value edges and a set of odd length cycles with half value edges.
  1415   ///
  1416   /// The algorithm calculates an optimal fractional matching and a
  1417   /// proof of the optimality. The solution of the dual problem can be
  1418   /// used to check the result of the algorithm. The dual linear
  1419   /// problem is the following.
  1420   /// \f[ y_u + y_v \ge w_{uv} \quad \forall uv\in E\f]
  1421   /// \f[\min \sum_{u \in V}y_u \f]
  1422   ///
  1423   /// The algorithm can be executed with the run() function.
  1424   /// After it the matching (the primal solution) and the dual solution
  1425   /// can be obtained using the query functions.
  1426 
  1427   /// If the value type is integer, then the primal and the dual
  1428   /// solutions are multiplied by
  1429   /// \ref MaxWeightedPerfectFractionalMatching::primalScale "2" and
  1430   /// \ref MaxWeightedPerfectFractionalMatching::dualScale "4" respectively.
  1431   ///
  1432   /// \tparam GR The undirected graph type the algorithm runs on.
  1433   /// \tparam WM The type edge weight map. The default type is
  1434   /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
  1435 #ifdef DOXYGEN
  1436   template <typename GR, typename WM>
  1437 #else
  1438   template <typename GR,
  1439             typename WM = typename GR::template EdgeMap<int> >
  1440 #endif
  1441   class MaxWeightedPerfectFractionalMatching {
  1442   public:
  1443 
  1444     /// The graph type of the algorithm
  1445     typedef GR Graph;
  1446     /// The type of the edge weight map
  1447     typedef WM WeightMap;
  1448     /// The value type of the edge weights
  1449     typedef typename WeightMap::Value Value;
  1450 
  1451     /// The type of the matching map
  1452     typedef typename Graph::template NodeMap<typename Graph::Arc>
  1453     MatchingMap;
  1454 
  1455     /// \brief Scaling factor for primal solution
  1456     ///
  1457     /// Scaling factor for primal solution. It is equal to 2 or 1
  1458     /// according to the value type.
  1459     static const int primalScale =
  1460       std::numeric_limits<Value>::is_integer ? 2 : 1;
  1461 
  1462     /// \brief Scaling factor for dual solution
  1463     ///
  1464     /// Scaling factor for dual solution. It is equal to 4 or 1
  1465     /// according to the value type.
  1466     static const int dualScale =
  1467       std::numeric_limits<Value>::is_integer ? 4 : 1;
  1468 
  1469   private:
  1470 
  1471     TEMPLATE_GRAPH_TYPEDEFS(Graph);
  1472 
  1473     typedef typename Graph::template NodeMap<Value> NodePotential;
  1474 
  1475     const Graph& _graph;
  1476     const WeightMap& _weight;
  1477 
  1478     MatchingMap* _matching;
  1479     NodePotential* _node_potential;
  1480 
  1481     int _node_num;
  1482     bool _allow_loops;
  1483 
  1484     enum Status {
  1485       EVEN = -1, MATCHED = 0, ODD = 1
  1486     };
  1487 
  1488     typedef typename Graph::template NodeMap<Status> StatusMap;
  1489     StatusMap* _status;
  1490 
  1491     typedef typename Graph::template NodeMap<Arc> PredMap;
  1492     PredMap* _pred;
  1493 
  1494     typedef ExtendFindEnum<IntNodeMap> TreeSet;
  1495 
  1496     IntNodeMap *_tree_set_index;
  1497     TreeSet *_tree_set;
  1498 
  1499     IntNodeMap *_delta2_index;
  1500     BinHeap<Value, IntNodeMap> *_delta2;
  1501 
  1502     IntEdgeMap *_delta3_index;
  1503     BinHeap<Value, IntEdgeMap> *_delta3;
  1504 
  1505     Value _delta_sum;
  1506 
  1507     void createStructures() {
  1508       _node_num = countNodes(_graph);
  1509 
  1510       if (!_matching) {
  1511         _matching = new MatchingMap(_graph);
  1512       }
  1513       if (!_node_potential) {
  1514         _node_potential = new NodePotential(_graph);
  1515       }
  1516       if (!_status) {
  1517         _status = new StatusMap(_graph);
  1518       }
  1519       if (!_pred) {
  1520         _pred = new PredMap(_graph);
  1521       }
  1522       if (!_tree_set) {
  1523         _tree_set_index = new IntNodeMap(_graph);
  1524         _tree_set = new TreeSet(*_tree_set_index);
  1525       }
  1526       if (!_delta2) {
  1527         _delta2_index = new IntNodeMap(_graph);
  1528         _delta2 = new BinHeap<Value, IntNodeMap>(*_delta2_index);
  1529       }
  1530       if (!_delta3) {
  1531         _delta3_index = new IntEdgeMap(_graph);
  1532         _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
  1533       }
  1534     }
  1535 
  1536     void destroyStructures() {
  1537       if (_matching) {
  1538         delete _matching;
  1539       }
  1540       if (_node_potential) {
  1541         delete _node_potential;
  1542       }
  1543       if (_status) {
  1544         delete _status;
  1545       }
  1546       if (_pred) {
  1547         delete _pred;
  1548       }
  1549       if (_tree_set) {
  1550         delete _tree_set_index;
  1551         delete _tree_set;
  1552       }
  1553       if (_delta2) {
  1554         delete _delta2_index;
  1555         delete _delta2;
  1556       }
  1557       if (_delta3) {
  1558         delete _delta3_index;
  1559         delete _delta3;
  1560       }
  1561     }
  1562 
  1563     void matchedToEven(Node node, int tree) {
  1564       _tree_set->insert(node, tree);
  1565       _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
  1566 
  1567       if (_delta2->state(node) == _delta2->IN_HEAP) {
  1568         _delta2->erase(node);
  1569       }
  1570 
  1571       for (InArcIt a(_graph, node); a != INVALID; ++a) {
  1572         Node v = _graph.source(a);
  1573         Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
  1574           dualScale * _weight[a];
  1575         if (node == v) {
  1576           if (_allow_loops && _graph.direction(a)) {
  1577             _delta3->push(a, rw / 2);
  1578           }
  1579         } else if ((*_status)[v] == EVEN) {
  1580           _delta3->push(a, rw / 2);
  1581         } else if ((*_status)[v] == MATCHED) {
  1582           if (_delta2->state(v) != _delta2->IN_HEAP) {
  1583             _pred->set(v, a);
  1584             _delta2->push(v, rw);
  1585           } else if ((*_delta2)[v] > rw) {
  1586             _pred->set(v, a);
  1587             _delta2->decrease(v, rw);
  1588           }
  1589         }
  1590       }
  1591     }
  1592 
  1593     void matchedToOdd(Node node, int tree) {
  1594       _tree_set->insert(node, tree);
  1595       _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
  1596 
  1597       if (_delta2->state(node) == _delta2->IN_HEAP) {
  1598         _delta2->erase(node);
  1599       }
  1600     }
  1601 
  1602     void evenToMatched(Node node, int tree) {
  1603       _node_potential->set(node, (*_node_potential)[node] - _delta_sum);
  1604       Arc min = INVALID;
  1605       Value minrw = std::numeric_limits<Value>::max();
  1606       for (InArcIt a(_graph, node); a != INVALID; ++a) {
  1607         Node v = _graph.source(a);
  1608         Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
  1609           dualScale * _weight[a];
  1610 
  1611         if (node == v) {
  1612           if (_allow_loops && _graph.direction(a)) {
  1613             _delta3->erase(a);
  1614           }
  1615         } else if ((*_status)[v] == EVEN) {
  1616           _delta3->erase(a);
  1617           if (minrw > rw) {
  1618             min = _graph.oppositeArc(a);
  1619             minrw = rw;
  1620           }
  1621         } else if ((*_status)[v]  == MATCHED) {
  1622           if ((*_pred)[v] == a) {
  1623             Arc mina = INVALID;
  1624             Value minrwa = std::numeric_limits<Value>::max();
  1625             for (OutArcIt aa(_graph, v); aa != INVALID; ++aa) {
  1626               Node va = _graph.target(aa);
  1627               if ((*_status)[va] != EVEN ||
  1628                   _tree_set->find(va) == tree) continue;
  1629               Value rwa = (*_node_potential)[v] + (*_node_potential)[va] -
  1630                 dualScale * _weight[aa];
  1631               if (minrwa > rwa) {
  1632                 minrwa = rwa;
  1633                 mina = aa;
  1634               }
  1635             }
  1636             if (mina != INVALID) {
  1637               _pred->set(v, mina);
  1638               _delta2->increase(v, minrwa);
  1639             } else {
  1640               _pred->set(v, INVALID);
  1641               _delta2->erase(v);
  1642             }
  1643           }
  1644         }
  1645       }
  1646       if (min != INVALID) {
  1647         _pred->set(node, min);
  1648         _delta2->push(node, minrw);
  1649       } else {
  1650         _pred->set(node, INVALID);
  1651       }
  1652     }
  1653 
  1654     void oddToMatched(Node node) {
  1655       _node_potential->set(node, (*_node_potential)[node] + _delta_sum);
  1656       Arc min = INVALID;
  1657       Value minrw = std::numeric_limits<Value>::max();
  1658       for (InArcIt a(_graph, node); a != INVALID; ++a) {
  1659         Node v = _graph.source(a);
  1660         if ((*_status)[v] != EVEN) continue;
  1661         Value rw = (*_node_potential)[node] + (*_node_potential)[v] -
  1662           dualScale * _weight[a];
  1663 
  1664         if (minrw > rw) {
  1665           min = _graph.oppositeArc(a);
  1666           minrw = rw;
  1667         }
  1668       }
  1669       if (min != INVALID) {
  1670         _pred->set(node, min);
  1671         _delta2->push(node, minrw);
  1672       } else {
  1673         _pred->set(node, INVALID);
  1674       }
  1675     }
  1676 
  1677     void alternatePath(Node even, int tree) {
  1678       Node odd;
  1679 
  1680       _status->set(even, MATCHED);
  1681       evenToMatched(even, tree);
  1682 
  1683       Arc prev = (*_matching)[even];
  1684       while (prev != INVALID) {
  1685         odd = _graph.target(prev);
  1686         even = _graph.target((*_pred)[odd]);
  1687         _matching->set(odd, (*_pred)[odd]);
  1688         _status->set(odd, MATCHED);
  1689         oddToMatched(odd);
  1690 
  1691         prev = (*_matching)[even];
  1692         _status->set(even, MATCHED);
  1693         _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
  1694         evenToMatched(even, tree);
  1695       }
  1696     }
  1697 
  1698     void destroyTree(int tree) {
  1699       for (typename TreeSet::ItemIt n(*_tree_set, tree); n != INVALID; ++n) {
  1700         if ((*_status)[n] == EVEN) {
  1701           _status->set(n, MATCHED);
  1702           evenToMatched(n, tree);
  1703         } else if ((*_status)[n] == ODD) {
  1704           _status->set(n, MATCHED);
  1705           oddToMatched(n);
  1706         }
  1707       }
  1708       _tree_set->eraseClass(tree);
  1709     }
  1710 
  1711     void augmentOnEdge(const Edge& edge) {
  1712       Node left = _graph.u(edge);
  1713       int left_tree = _tree_set->find(left);
  1714 
  1715       alternatePath(left, left_tree);
  1716       destroyTree(left_tree);
  1717       _matching->set(left, _graph.direct(edge, true));
  1718 
  1719       Node right = _graph.v(edge);
  1720       int right_tree = _tree_set->find(right);
  1721 
  1722       alternatePath(right, right_tree);
  1723       destroyTree(right_tree);
  1724       _matching->set(right, _graph.direct(edge, false));
  1725     }
  1726 
  1727     void augmentOnArc(const Arc& arc) {
  1728       Node left = _graph.source(arc);
  1729       _status->set(left, MATCHED);
  1730       _matching->set(left, arc);
  1731       _pred->set(left, arc);
  1732 
  1733       Node right = _graph.target(arc);
  1734       int right_tree = _tree_set->find(right);
  1735 
  1736       alternatePath(right, right_tree);
  1737       destroyTree(right_tree);
  1738       _matching->set(right, _graph.oppositeArc(arc));
  1739     }
  1740 
  1741     void extendOnArc(const Arc& arc) {
  1742       Node base = _graph.target(arc);
  1743       int tree = _tree_set->find(base);
  1744 
  1745       Node odd = _graph.source(arc);
  1746       _tree_set->insert(odd, tree);
  1747       _status->set(odd, ODD);
  1748       matchedToOdd(odd, tree);
  1749       _pred->set(odd, arc);
  1750 
  1751       Node even = _graph.target((*_matching)[odd]);
  1752       _tree_set->insert(even, tree);
  1753       _status->set(even, EVEN);
  1754       matchedToEven(even, tree);
  1755     }
  1756 
  1757     void cycleOnEdge(const Edge& edge, int tree) {
  1758       Node nca = INVALID;
  1759       std::vector<Node> left_path, right_path;
  1760 
  1761       {
  1762         std::set<Node> left_set, right_set;
  1763         Node left = _graph.u(edge);
  1764         left_path.push_back(left);
  1765         left_set.insert(left);
  1766 
  1767         Node right = _graph.v(edge);
  1768         right_path.push_back(right);
  1769         right_set.insert(right);
  1770 
  1771         while (true) {
  1772 
  1773           if (left_set.find(right) != left_set.end()) {
  1774             nca = right;
  1775             break;
  1776           }
  1777 
  1778           if ((*_matching)[left] == INVALID) break;
  1779 
  1780           left = _graph.target((*_matching)[left]);
  1781           left_path.push_back(left);
  1782           left = _graph.target((*_pred)[left]);
  1783           left_path.push_back(left);
  1784 
  1785           left_set.insert(left);
  1786 
  1787           if (right_set.find(left) != right_set.end()) {
  1788             nca = left;
  1789             break;
  1790           }
  1791 
  1792           if ((*_matching)[right] == INVALID) break;
  1793 
  1794           right = _graph.target((*_matching)[right]);
  1795           right_path.push_back(right);
  1796           right = _graph.target((*_pred)[right]);
  1797           right_path.push_back(right);
  1798 
  1799           right_set.insert(right);
  1800 
  1801         }
  1802 
  1803         if (nca == INVALID) {
  1804           if ((*_matching)[left] == INVALID) {
  1805             nca = right;
  1806             while (left_set.find(nca) == left_set.end()) {
  1807               nca = _graph.target((*_matching)[nca]);
  1808               right_path.push_back(nca);
  1809               nca = _graph.target((*_pred)[nca]);
  1810               right_path.push_back(nca);
  1811             }
  1812           } else {
  1813             nca = left;
  1814             while (right_set.find(nca) == right_set.end()) {
  1815               nca = _graph.target((*_matching)[nca]);
  1816               left_path.push_back(nca);
  1817               nca = _graph.target((*_pred)[nca]);
  1818               left_path.push_back(nca);
  1819             }
  1820           }
  1821         }
  1822       }
  1823 
  1824       alternatePath(nca, tree);
  1825       Arc prev;
  1826 
  1827       prev = _graph.direct(edge, true);
  1828       for (int i = 0; left_path[i] != nca; i += 2) {
  1829         _matching->set(left_path[i], prev);
  1830         _status->set(left_path[i], MATCHED);
  1831         evenToMatched(left_path[i], tree);
  1832 
  1833         prev = _graph.oppositeArc((*_pred)[left_path[i + 1]]);
  1834         _status->set(left_path[i + 1], MATCHED);
  1835         oddToMatched(left_path[i + 1]);
  1836       }
  1837       _matching->set(nca, prev);
  1838 
  1839       for (int i = 0; right_path[i] != nca; i += 2) {
  1840         _status->set(right_path[i], MATCHED);
  1841         evenToMatched(right_path[i], tree);
  1842 
  1843         _matching->set(right_path[i + 1], (*_pred)[right_path[i + 1]]);
  1844         _status->set(right_path[i + 1], MATCHED);
  1845         oddToMatched(right_path[i + 1]);
  1846       }
  1847 
  1848       destroyTree(tree);
  1849     }
  1850 
  1851     void extractCycle(const Arc &arc) {
  1852       Node left = _graph.source(arc);
  1853       Node odd = _graph.target((*_matching)[left]);
  1854       Arc prev;
  1855       while (odd != left) {
  1856         Node even = _graph.target((*_matching)[odd]);
  1857         prev = (*_matching)[odd];
  1858         odd = _graph.target((*_matching)[even]);
  1859         _matching->set(even, _graph.oppositeArc(prev));
  1860       }
  1861       _matching->set(left, arc);
  1862 
  1863       Node right = _graph.target(arc);
  1864       int right_tree = _tree_set->find(right);
  1865       alternatePath(right, right_tree);
  1866       destroyTree(right_tree);
  1867       _matching->set(right, _graph.oppositeArc(arc));
  1868     }
  1869 
  1870   public:
  1871 
  1872     /// \brief Constructor
  1873     ///
  1874     /// Constructor.
  1875     MaxWeightedPerfectFractionalMatching(const Graph& graph,
  1876                                          const WeightMap& weight,
  1877                                          bool allow_loops = true)
  1878       : _graph(graph), _weight(weight), _matching(0),
  1879       _node_potential(0), _node_num(0), _allow_loops(allow_loops),
  1880       _status(0),  _pred(0),
  1881       _tree_set_index(0), _tree_set(0),
  1882 
  1883       _delta2_index(0), _delta2(0),
  1884       _delta3_index(0), _delta3(0),
  1885 
  1886       _delta_sum() {}
  1887 
  1888     ~MaxWeightedPerfectFractionalMatching() {
  1889       destroyStructures();
  1890     }
  1891 
  1892     /// \name Execution Control
  1893     /// The simplest way to execute the algorithm is to use the
  1894     /// \ref run() member function.
  1895 
  1896     ///@{
  1897 
  1898     /// \brief Initialize the algorithm
  1899     ///
  1900     /// This function initializes the algorithm.
  1901     void init() {
  1902       createStructures();
  1903 
  1904       for (NodeIt n(_graph); n != INVALID; ++n) {
  1905         (*_delta2_index)[n] = _delta2->PRE_HEAP;
  1906       }
  1907       for (EdgeIt e(_graph); e != INVALID; ++e) {
  1908         (*_delta3_index)[e] = _delta3->PRE_HEAP;
  1909       }
  1910 
  1911       for (NodeIt n(_graph); n != INVALID; ++n) {
  1912         Value max = - std::numeric_limits<Value>::max();
  1913         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  1914           if (_graph.target(e) == n && !_allow_loops) continue;
  1915           if ((dualScale * _weight[e]) / 2 > max) {
  1916             max = (dualScale * _weight[e]) / 2;
  1917           }
  1918         }
  1919         _node_potential->set(n, max);
  1920 
  1921         _tree_set->insert(n);
  1922 
  1923         _matching->set(n, INVALID);
  1924         _status->set(n, EVEN);
  1925       }
  1926 
  1927       for (EdgeIt e(_graph); e != INVALID; ++e) {
  1928         Node left = _graph.u(e);
  1929         Node right = _graph.v(e);
  1930         if (left == right && !_allow_loops) continue;
  1931         _delta3->push(e, ((*_node_potential)[left] +
  1932                           (*_node_potential)[right] -
  1933                           dualScale * _weight[e]) / 2);
  1934       }
  1935     }
  1936 
  1937     /// \brief Start the algorithm
  1938     ///
  1939     /// This function starts the algorithm.
  1940     ///
  1941     /// \pre \ref init() must be called before using this function.
  1942     bool start() {
  1943       enum OpType {
  1944         D2, D3
  1945       };
  1946 
  1947       int unmatched = _node_num;
  1948       while (unmatched > 0) {
  1949         Value d2 = !_delta2->empty() ?
  1950           _delta2->prio() : std::numeric_limits<Value>::max();
  1951 
  1952         Value d3 = !_delta3->empty() ?
  1953           _delta3->prio() : std::numeric_limits<Value>::max();
  1954 
  1955         _delta_sum = d3; OpType ot = D3;
  1956         if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
  1957 
  1958         if (_delta_sum == std::numeric_limits<Value>::max()) {
  1959           return false;
  1960         }
  1961 
  1962         switch (ot) {
  1963         case D2:
  1964           {
  1965             Node n = _delta2->top();
  1966             Arc a = (*_pred)[n];
  1967             if ((*_matching)[n] == INVALID) {
  1968               augmentOnArc(a);
  1969               --unmatched;
  1970             } else {
  1971               Node v = _graph.target((*_matching)[n]);
  1972               if ((*_matching)[n] !=
  1973                   _graph.oppositeArc((*_matching)[v])) {
  1974                 extractCycle(a);
  1975                 --unmatched;
  1976               } else {
  1977                 extendOnArc(a);
  1978               }
  1979             }
  1980           } break;
  1981         case D3:
  1982           {
  1983             Edge e = _delta3->top();
  1984 
  1985             Node left = _graph.u(e);
  1986             Node right = _graph.v(e);
  1987 
  1988             int left_tree = _tree_set->find(left);
  1989             int right_tree = _tree_set->find(right);
  1990 
  1991             if (left_tree == right_tree) {
  1992               cycleOnEdge(e, left_tree);
  1993               --unmatched;
  1994             } else {
  1995               augmentOnEdge(e);
  1996               unmatched -= 2;
  1997             }
  1998           } break;
  1999         }
  2000       }
  2001       return true;
  2002     }
  2003 
  2004     /// \brief Run the algorithm.
  2005     ///
  2006     /// This method runs the \c %MaxWeightedPerfectFractionalMatching 
  2007     /// algorithm.
  2008     ///
  2009     /// \note mwfm.run() is just a shortcut of the following code.
  2010     /// \code
  2011     ///   mwpfm.init();
  2012     ///   mwpfm.start();
  2013     /// \endcode
  2014     bool run() {
  2015       init();
  2016       return start();
  2017     }
  2018 
  2019     /// @}
  2020 
  2021     /// \name Primal Solution
  2022     /// Functions to get the primal solution, i.e. the maximum weighted
  2023     /// matching.\n
  2024     /// Either \ref run() or \ref start() function should be called before
  2025     /// using them.
  2026 
  2027     /// @{
  2028 
  2029     /// \brief Return the weight of the matching.
  2030     ///
  2031     /// This function returns the weight of the found matching. This
  2032     /// value is scaled by \ref primalScale "primal scale".
  2033     ///
  2034     /// \pre Either run() or start() must be called before using this function.
  2035     Value matchingWeight() const {
  2036       Value sum = 0;
  2037       for (NodeIt n(_graph); n != INVALID; ++n) {
  2038         if ((*_matching)[n] != INVALID) {
  2039           sum += _weight[(*_matching)[n]];
  2040         }
  2041       }
  2042       return sum * primalScale / 2;
  2043     }
  2044 
  2045     /// \brief Return the number of covered nodes in the matching.
  2046     ///
  2047     /// This function returns the number of covered nodes in the matching.
  2048     ///
  2049     /// \pre Either run() or start() must be called before using this function.
  2050     int matchingSize() const {
  2051       int num = 0;
  2052       for (NodeIt n(_graph); n != INVALID; ++n) {
  2053         if ((*_matching)[n] != INVALID) {
  2054           ++num;
  2055         }
  2056       }
  2057       return num;
  2058     }
  2059 
  2060     /// \brief Return \c true if the given edge is in the matching.
  2061     ///
  2062     /// This function returns \c true if the given edge is in the
  2063     /// found matching. The result is scaled by \ref primalScale
  2064     /// "primal scale".
  2065     ///
  2066     /// \pre Either run() or start() must be called before using this function.
  2067     Value matching(const Edge& edge) const {
  2068       return Value(edge == (*_matching)[_graph.u(edge)] ? 1 : 0)
  2069         * primalScale / 2 + Value(edge == (*_matching)[_graph.v(edge)] ? 1 : 0)
  2070         * primalScale / 2;
  2071     }
  2072 
  2073     /// \brief Return the fractional matching arc (or edge) incident
  2074     /// to the given node.
  2075     ///
  2076     /// This function returns one of the fractional matching arc (or
  2077     /// edge) incident to the given node in the found matching or \c
  2078     /// INVALID if the node is not covered by the matching or if the
  2079     /// node is on an odd length cycle then it is the successor edge
  2080     /// on the cycle.
  2081     ///
  2082     /// \pre Either run() or start() must be called before using this function.
  2083     Arc matching(const Node& node) const {
  2084       return (*_matching)[node];
  2085     }
  2086 
  2087     /// \brief Return a const reference to the matching map.
  2088     ///
  2089     /// This function returns a const reference to a node map that stores
  2090     /// the matching arc (or edge) incident to each node.
  2091     const MatchingMap& matchingMap() const {
  2092       return *_matching;
  2093     }
  2094 
  2095     /// @}
  2096 
  2097     /// \name Dual Solution
  2098     /// Functions to get the dual solution.\n
  2099     /// Either \ref run() or \ref start() function should be called before
  2100     /// using them.
  2101 
  2102     /// @{
  2103 
  2104     /// \brief Return the value of the dual solution.
  2105     ///
  2106     /// This function returns the value of the dual solution.
  2107     /// It should be equal to the primal value scaled by \ref dualScale
  2108     /// "dual scale".
  2109     ///
  2110     /// \pre Either run() or start() must be called before using this function.
  2111     Value dualValue() const {
  2112       Value sum = 0;
  2113       for (NodeIt n(_graph); n != INVALID; ++n) {
  2114         sum += nodeValue(n);
  2115       }
  2116       return sum;
  2117     }
  2118 
  2119     /// \brief Return the dual value (potential) of the given node.
  2120     ///
  2121     /// This function returns the dual value (potential) of the given node.
  2122     ///
  2123     /// \pre Either run() or start() must be called before using this function.
  2124     Value nodeValue(const Node& n) const {
  2125       return (*_node_potential)[n];
  2126     }
  2127 
  2128     /// @}
  2129 
  2130   };
  2131 
  2132 } //END OF NAMESPACE LEMON
  2133 
  2134 #endif //LEMON_FRACTIONAL_MATCHING_H