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