lemon/matching.h
author Alpar Juttner <alpar@cs.elte.hu>
Wed, 28 Nov 2012 11:41:40 +0100
changeset 1007 7e368d9b67f7
parent 651 3adf5e2d1e62
child 875 07ec2b52e53d
permissions -rw-r--r--
Avoid GCC 4.7 compiler warnings (#453)
     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_MAX_MATCHING_H
    20 #define LEMON_MAX_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 
    32 ///\ingroup matching
    33 ///\file
    34 ///\brief Maximum matching algorithms in general graphs.
    35 
    36 namespace lemon {
    37 
    38   /// \ingroup matching
    39   ///
    40   /// \brief Maximum cardinality matching in general graphs
    41   ///
    42   /// This class implements Edmonds' alternating forest matching algorithm
    43   /// for finding a maximum cardinality matching in a general undirected graph.
    44   /// It can be started from an arbitrary initial matching 
    45   /// (the default is the empty one).
    46   ///
    47   /// The dual solution of the problem is a map of the nodes to
    48   /// \ref MaxMatching::Status "Status", having values \c EVEN (or \c D),
    49   /// \c ODD (or \c A) and \c MATCHED (or \c C) defining the Gallai-Edmonds
    50   /// decomposition of the graph. The nodes in \c EVEN/D induce a subgraph
    51   /// with factor-critical components, the nodes in \c ODD/A form the
    52   /// canonical barrier, and the nodes in \c MATCHED/C induce a graph having
    53   /// a perfect matching. The number of the factor-critical components
    54   /// minus the number of barrier nodes is a lower bound on the
    55   /// unmatched nodes, and the matching is optimal if and only if this bound is
    56   /// tight. This decomposition can be obtained using \ref status() or
    57   /// \ref statusMap() after running the algorithm.
    58   ///
    59   /// \tparam GR The undirected graph type the algorithm runs on.
    60   template <typename GR>
    61   class MaxMatching {
    62   public:
    63 
    64     /// The graph type of the algorithm
    65     typedef GR Graph;
    66     /// The type of the matching map
    67     typedef typename Graph::template NodeMap<typename Graph::Arc>
    68     MatchingMap;
    69 
    70     ///\brief Status constants for Gallai-Edmonds decomposition.
    71     ///
    72     ///These constants are used for indicating the Gallai-Edmonds 
    73     ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
    74     ///induce a subgraph with factor-critical components, the nodes with
    75     ///status \c ODD (or \c A) form the canonical barrier, and the nodes
    76     ///with status \c MATCHED (or \c C) induce a subgraph having a 
    77     ///perfect matching.
    78     enum Status {
    79       EVEN = 1,       ///< = 1. (\c D is an alias for \c EVEN.)
    80       D = 1,
    81       MATCHED = 0,    ///< = 0. (\c C is an alias for \c MATCHED.)
    82       C = 0,
    83       ODD = -1,       ///< = -1. (\c A is an alias for \c ODD.)
    84       A = -1,
    85       UNMATCHED = -2  ///< = -2.
    86     };
    87 
    88     /// The type of the status map
    89     typedef typename Graph::template NodeMap<Status> StatusMap;
    90 
    91   private:
    92 
    93     TEMPLATE_GRAPH_TYPEDEFS(Graph);
    94 
    95     typedef UnionFindEnum<IntNodeMap> BlossomSet;
    96     typedef ExtendFindEnum<IntNodeMap> TreeSet;
    97     typedef RangeMap<Node> NodeIntMap;
    98     typedef MatchingMap EarMap;
    99     typedef std::vector<Node> NodeQueue;
   100 
   101     const Graph& _graph;
   102     MatchingMap* _matching;
   103     StatusMap* _status;
   104 
   105     EarMap* _ear;
   106 
   107     IntNodeMap* _blossom_set_index;
   108     BlossomSet* _blossom_set;
   109     NodeIntMap* _blossom_rep;
   110 
   111     IntNodeMap* _tree_set_index;
   112     TreeSet* _tree_set;
   113 
   114     NodeQueue _node_queue;
   115     int _process, _postpone, _last;
   116 
   117     int _node_num;
   118 
   119   private:
   120 
   121     void createStructures() {
   122       _node_num = countNodes(_graph);
   123       if (!_matching) {
   124         _matching = new MatchingMap(_graph);
   125       }
   126       if (!_status) {
   127         _status = new StatusMap(_graph);
   128       }
   129       if (!_ear) {
   130         _ear = new EarMap(_graph);
   131       }
   132       if (!_blossom_set) {
   133         _blossom_set_index = new IntNodeMap(_graph);
   134         _blossom_set = new BlossomSet(*_blossom_set_index);
   135       }
   136       if (!_blossom_rep) {
   137         _blossom_rep = new NodeIntMap(_node_num);
   138       }
   139       if (!_tree_set) {
   140         _tree_set_index = new IntNodeMap(_graph);
   141         _tree_set = new TreeSet(*_tree_set_index);
   142       }
   143       _node_queue.resize(_node_num);
   144     }
   145 
   146     void destroyStructures() {
   147       if (_matching) {
   148         delete _matching;
   149       }
   150       if (_status) {
   151         delete _status;
   152       }
   153       if (_ear) {
   154         delete _ear;
   155       }
   156       if (_blossom_set) {
   157         delete _blossom_set;
   158         delete _blossom_set_index;
   159       }
   160       if (_blossom_rep) {
   161         delete _blossom_rep;
   162       }
   163       if (_tree_set) {
   164         delete _tree_set_index;
   165         delete _tree_set;
   166       }
   167     }
   168 
   169     void processDense(const Node& n) {
   170       _process = _postpone = _last = 0;
   171       _node_queue[_last++] = n;
   172 
   173       while (_process != _last) {
   174         Node u = _node_queue[_process++];
   175         for (OutArcIt a(_graph, u); a != INVALID; ++a) {
   176           Node v = _graph.target(a);
   177           if ((*_status)[v] == MATCHED) {
   178             extendOnArc(a);
   179           } else if ((*_status)[v] == UNMATCHED) {
   180             augmentOnArc(a);
   181             return;
   182           }
   183         }
   184       }
   185 
   186       while (_postpone != _last) {
   187         Node u = _node_queue[_postpone++];
   188 
   189         for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
   190           Node v = _graph.target(a);
   191 
   192           if ((*_status)[v] == EVEN) {
   193             if (_blossom_set->find(u) != _blossom_set->find(v)) {
   194               shrinkOnEdge(a);
   195             }
   196           }
   197 
   198           while (_process != _last) {
   199             Node w = _node_queue[_process++];
   200             for (OutArcIt b(_graph, w); b != INVALID; ++b) {
   201               Node x = _graph.target(b);
   202               if ((*_status)[x] == MATCHED) {
   203                 extendOnArc(b);
   204               } else if ((*_status)[x] == UNMATCHED) {
   205                 augmentOnArc(b);
   206                 return;
   207               }
   208             }
   209           }
   210         }
   211       }
   212     }
   213 
   214     void processSparse(const Node& n) {
   215       _process = _last = 0;
   216       _node_queue[_last++] = n;
   217       while (_process != _last) {
   218         Node u = _node_queue[_process++];
   219         for (OutArcIt a(_graph, u); a != INVALID; ++a) {
   220           Node v = _graph.target(a);
   221 
   222           if ((*_status)[v] == EVEN) {
   223             if (_blossom_set->find(u) != _blossom_set->find(v)) {
   224               shrinkOnEdge(a);
   225             }
   226           } else if ((*_status)[v] == MATCHED) {
   227             extendOnArc(a);
   228           } else if ((*_status)[v] == UNMATCHED) {
   229             augmentOnArc(a);
   230             return;
   231           }
   232         }
   233       }
   234     }
   235 
   236     void shrinkOnEdge(const Edge& e) {
   237       Node nca = INVALID;
   238 
   239       {
   240         std::set<Node> left_set, right_set;
   241 
   242         Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
   243         left_set.insert(left);
   244 
   245         Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
   246         right_set.insert(right);
   247 
   248         while (true) {
   249           if ((*_matching)[left] == INVALID) break;
   250           left = _graph.target((*_matching)[left]);
   251           left = (*_blossom_rep)[_blossom_set->
   252                                  find(_graph.target((*_ear)[left]))];
   253           if (right_set.find(left) != right_set.end()) {
   254             nca = left;
   255             break;
   256           }
   257           left_set.insert(left);
   258 
   259           if ((*_matching)[right] == INVALID) break;
   260           right = _graph.target((*_matching)[right]);
   261           right = (*_blossom_rep)[_blossom_set->
   262                                   find(_graph.target((*_ear)[right]))];
   263           if (left_set.find(right) != left_set.end()) {
   264             nca = right;
   265             break;
   266           }
   267           right_set.insert(right);
   268         }
   269 
   270         if (nca == INVALID) {
   271           if ((*_matching)[left] == INVALID) {
   272             nca = right;
   273             while (left_set.find(nca) == left_set.end()) {
   274               nca = _graph.target((*_matching)[nca]);
   275               nca =(*_blossom_rep)[_blossom_set->
   276                                    find(_graph.target((*_ear)[nca]))];
   277             }
   278           } else {
   279             nca = left;
   280             while (right_set.find(nca) == right_set.end()) {
   281               nca = _graph.target((*_matching)[nca]);
   282               nca = (*_blossom_rep)[_blossom_set->
   283                                    find(_graph.target((*_ear)[nca]))];
   284             }
   285           }
   286         }
   287       }
   288 
   289       {
   290 
   291         Node node = _graph.u(e);
   292         Arc arc = _graph.direct(e, true);
   293         Node base = (*_blossom_rep)[_blossom_set->find(node)];
   294 
   295         while (base != nca) {
   296           (*_ear)[node] = arc;
   297 
   298           Node n = node;
   299           while (n != base) {
   300             n = _graph.target((*_matching)[n]);
   301             Arc a = (*_ear)[n];
   302             n = _graph.target(a);
   303             (*_ear)[n] = _graph.oppositeArc(a);
   304           }
   305           node = _graph.target((*_matching)[base]);
   306           _tree_set->erase(base);
   307           _tree_set->erase(node);
   308           _blossom_set->insert(node, _blossom_set->find(base));
   309           (*_status)[node] = EVEN;
   310           _node_queue[_last++] = node;
   311           arc = _graph.oppositeArc((*_ear)[node]);
   312           node = _graph.target((*_ear)[node]);
   313           base = (*_blossom_rep)[_blossom_set->find(node)];
   314           _blossom_set->join(_graph.target(arc), base);
   315         }
   316       }
   317 
   318       (*_blossom_rep)[_blossom_set->find(nca)] = nca;
   319 
   320       {
   321 
   322         Node node = _graph.v(e);
   323         Arc arc = _graph.direct(e, false);
   324         Node base = (*_blossom_rep)[_blossom_set->find(node)];
   325 
   326         while (base != nca) {
   327           (*_ear)[node] = arc;
   328 
   329           Node n = node;
   330           while (n != base) {
   331             n = _graph.target((*_matching)[n]);
   332             Arc a = (*_ear)[n];
   333             n = _graph.target(a);
   334             (*_ear)[n] = _graph.oppositeArc(a);
   335           }
   336           node = _graph.target((*_matching)[base]);
   337           _tree_set->erase(base);
   338           _tree_set->erase(node);
   339           _blossom_set->insert(node, _blossom_set->find(base));
   340           (*_status)[node] = EVEN;
   341           _node_queue[_last++] = node;
   342           arc = _graph.oppositeArc((*_ear)[node]);
   343           node = _graph.target((*_ear)[node]);
   344           base = (*_blossom_rep)[_blossom_set->find(node)];
   345           _blossom_set->join(_graph.target(arc), base);
   346         }
   347       }
   348 
   349       (*_blossom_rep)[_blossom_set->find(nca)] = nca;
   350     }
   351 
   352     void extendOnArc(const Arc& a) {
   353       Node base = _graph.source(a);
   354       Node odd = _graph.target(a);
   355 
   356       (*_ear)[odd] = _graph.oppositeArc(a);
   357       Node even = _graph.target((*_matching)[odd]);
   358       (*_blossom_rep)[_blossom_set->insert(even)] = even;
   359       (*_status)[odd] = ODD;
   360       (*_status)[even] = EVEN;
   361       int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
   362       _tree_set->insert(odd, tree);
   363       _tree_set->insert(even, tree);
   364       _node_queue[_last++] = even;
   365 
   366     }
   367 
   368     void augmentOnArc(const Arc& a) {
   369       Node even = _graph.source(a);
   370       Node odd = _graph.target(a);
   371 
   372       int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
   373 
   374       (*_matching)[odd] = _graph.oppositeArc(a);
   375       (*_status)[odd] = MATCHED;
   376 
   377       Arc arc = (*_matching)[even];
   378       (*_matching)[even] = a;
   379 
   380       while (arc != INVALID) {
   381         odd = _graph.target(arc);
   382         arc = (*_ear)[odd];
   383         even = _graph.target(arc);
   384         (*_matching)[odd] = arc;
   385         arc = (*_matching)[even];
   386         (*_matching)[even] = _graph.oppositeArc((*_matching)[odd]);
   387       }
   388 
   389       for (typename TreeSet::ItemIt it(*_tree_set, tree);
   390            it != INVALID; ++it) {
   391         if ((*_status)[it] == ODD) {
   392           (*_status)[it] = MATCHED;
   393         } else {
   394           int blossom = _blossom_set->find(it);
   395           for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
   396                jt != INVALID; ++jt) {
   397             (*_status)[jt] = MATCHED;
   398           }
   399           _blossom_set->eraseClass(blossom);
   400         }
   401       }
   402       _tree_set->eraseClass(tree);
   403 
   404     }
   405 
   406   public:
   407 
   408     /// \brief Constructor
   409     ///
   410     /// Constructor.
   411     MaxMatching(const Graph& graph)
   412       : _graph(graph), _matching(0), _status(0), _ear(0),
   413         _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
   414         _tree_set_index(0), _tree_set(0) {}
   415 
   416     ~MaxMatching() {
   417       destroyStructures();
   418     }
   419 
   420     /// \name Execution Control
   421     /// The simplest way to execute the algorithm is to use the
   422     /// \c run() member function.\n
   423     /// If you need better control on the execution, you have to call
   424     /// one of the functions \ref init(), \ref greedyInit() or
   425     /// \ref matchingInit() first, then you can start the algorithm with
   426     /// \ref startSparse() or \ref startDense().
   427 
   428     ///@{
   429 
   430     /// \brief Set the initial matching to the empty matching.
   431     ///
   432     /// This function sets the initial matching to the empty matching.
   433     void init() {
   434       createStructures();
   435       for(NodeIt n(_graph); n != INVALID; ++n) {
   436         (*_matching)[n] = INVALID;
   437         (*_status)[n] = UNMATCHED;
   438       }
   439     }
   440 
   441     /// \brief Find an initial matching in a greedy way.
   442     ///
   443     /// This function finds an initial matching in a greedy way.
   444     void greedyInit() {
   445       createStructures();
   446       for (NodeIt n(_graph); n != INVALID; ++n) {
   447         (*_matching)[n] = INVALID;
   448         (*_status)[n] = UNMATCHED;
   449       }
   450       for (NodeIt n(_graph); n != INVALID; ++n) {
   451         if ((*_matching)[n] == INVALID) {
   452           for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
   453             Node v = _graph.target(a);
   454             if ((*_matching)[v] == INVALID && v != n) {
   455               (*_matching)[n] = a;
   456               (*_status)[n] = MATCHED;
   457               (*_matching)[v] = _graph.oppositeArc(a);
   458               (*_status)[v] = MATCHED;
   459               break;
   460             }
   461           }
   462         }
   463       }
   464     }
   465 
   466 
   467     /// \brief Initialize the matching from a map.
   468     ///
   469     /// This function initializes the matching from a \c bool valued edge
   470     /// map. This map should have the property that there are no two incident
   471     /// edges with \c true value, i.e. it really contains a matching.
   472     /// \return \c true if the map contains a matching.
   473     template <typename MatchingMap>
   474     bool matchingInit(const MatchingMap& matching) {
   475       createStructures();
   476 
   477       for (NodeIt n(_graph); n != INVALID; ++n) {
   478         (*_matching)[n] = INVALID;
   479         (*_status)[n] = UNMATCHED;
   480       }
   481       for(EdgeIt e(_graph); e!=INVALID; ++e) {
   482         if (matching[e]) {
   483 
   484           Node u = _graph.u(e);
   485           if ((*_matching)[u] != INVALID) return false;
   486           (*_matching)[u] = _graph.direct(e, true);
   487           (*_status)[u] = MATCHED;
   488 
   489           Node v = _graph.v(e);
   490           if ((*_matching)[v] != INVALID) return false;
   491           (*_matching)[v] = _graph.direct(e, false);
   492           (*_status)[v] = MATCHED;
   493         }
   494       }
   495       return true;
   496     }
   497 
   498     /// \brief Start Edmonds' algorithm
   499     ///
   500     /// This function runs the original Edmonds' algorithm.
   501     ///
   502     /// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be
   503     /// called before using this function.
   504     void startSparse() {
   505       for(NodeIt n(_graph); n != INVALID; ++n) {
   506         if ((*_status)[n] == UNMATCHED) {
   507           (*_blossom_rep)[_blossom_set->insert(n)] = n;
   508           _tree_set->insert(n);
   509           (*_status)[n] = EVEN;
   510           processSparse(n);
   511         }
   512       }
   513     }
   514 
   515     /// \brief Start Edmonds' algorithm with a heuristic improvement 
   516     /// for dense graphs
   517     ///
   518     /// This function runs Edmonds' algorithm with a heuristic of postponing
   519     /// shrinks, therefore resulting in a faster algorithm for dense graphs.
   520     ///
   521     /// \pre \ref init(), \ref greedyInit() or \ref matchingInit() must be
   522     /// called before using this function.
   523     void startDense() {
   524       for(NodeIt n(_graph); n != INVALID; ++n) {
   525         if ((*_status)[n] == UNMATCHED) {
   526           (*_blossom_rep)[_blossom_set->insert(n)] = n;
   527           _tree_set->insert(n);
   528           (*_status)[n] = EVEN;
   529           processDense(n);
   530         }
   531       }
   532     }
   533 
   534 
   535     /// \brief Run Edmonds' algorithm
   536     ///
   537     /// This function runs Edmonds' algorithm. An additional heuristic of 
   538     /// postponing shrinks is used for relatively dense graphs 
   539     /// (for which <tt>m>=2*n</tt> holds).
   540     void run() {
   541       if (countEdges(_graph) < 2 * countNodes(_graph)) {
   542         greedyInit();
   543         startSparse();
   544       } else {
   545         init();
   546         startDense();
   547       }
   548     }
   549 
   550     /// @}
   551 
   552     /// \name Primal Solution
   553     /// Functions to get the primal solution, i.e. the maximum matching.
   554 
   555     /// @{
   556 
   557     /// \brief Return the size (cardinality) of the matching.
   558     ///
   559     /// This function returns the size (cardinality) of the current matching. 
   560     /// After run() it returns the size of the maximum matching in the graph.
   561     int matchingSize() const {
   562       int size = 0;
   563       for (NodeIt n(_graph); n != INVALID; ++n) {
   564         if ((*_matching)[n] != INVALID) {
   565           ++size;
   566         }
   567       }
   568       return size / 2;
   569     }
   570 
   571     /// \brief Return \c true if the given edge is in the matching.
   572     ///
   573     /// This function returns \c true if the given edge is in the current 
   574     /// matching.
   575     bool matching(const Edge& edge) const {
   576       return edge == (*_matching)[_graph.u(edge)];
   577     }
   578 
   579     /// \brief Return the matching arc (or edge) incident to the given node.
   580     ///
   581     /// This function returns the matching arc (or edge) incident to the
   582     /// given node in the current matching or \c INVALID if the node is 
   583     /// not covered by the matching.
   584     Arc matching(const Node& n) const {
   585       return (*_matching)[n];
   586     }
   587 
   588     /// \brief Return a const reference to the matching map.
   589     ///
   590     /// This function returns a const reference to a node map that stores
   591     /// the matching arc (or edge) incident to each node.
   592     const MatchingMap& matchingMap() const {
   593       return *_matching;
   594     }
   595 
   596     /// \brief Return the mate of the given node.
   597     ///
   598     /// This function returns the mate of the given node in the current 
   599     /// matching or \c INVALID if the node is not covered by the matching.
   600     Node mate(const Node& n) const {
   601       return (*_matching)[n] != INVALID ?
   602         _graph.target((*_matching)[n]) : INVALID;
   603     }
   604 
   605     /// @}
   606 
   607     /// \name Dual Solution
   608     /// Functions to get the dual solution, i.e. the Gallai-Edmonds 
   609     /// decomposition.
   610 
   611     /// @{
   612 
   613     /// \brief Return the status of the given node in the Edmonds-Gallai
   614     /// decomposition.
   615     ///
   616     /// This function returns the \ref Status "status" of the given node
   617     /// in the Edmonds-Gallai decomposition.
   618     Status status(const Node& n) const {
   619       return (*_status)[n];
   620     }
   621 
   622     /// \brief Return a const reference to the status map, which stores
   623     /// the Edmonds-Gallai decomposition.
   624     ///
   625     /// This function returns a const reference to a node map that stores the
   626     /// \ref Status "status" of each node in the Edmonds-Gallai decomposition.
   627     const StatusMap& statusMap() const {
   628       return *_status;
   629     }
   630 
   631     /// \brief Return \c true if the given node is in the barrier.
   632     ///
   633     /// This function returns \c true if the given node is in the barrier.
   634     bool barrier(const Node& n) const {
   635       return (*_status)[n] == ODD;
   636     }
   637 
   638     /// @}
   639 
   640   };
   641 
   642   /// \ingroup matching
   643   ///
   644   /// \brief Weighted matching in general graphs
   645   ///
   646   /// This class provides an efficient implementation of Edmond's
   647   /// maximum weighted matching algorithm. The implementation is based
   648   /// on extensive use of priority queues and provides
   649   /// \f$O(nm\log n)\f$ time complexity.
   650   ///
   651   /// The maximum weighted matching problem is to find a subset of the 
   652   /// edges in an undirected graph with maximum overall weight for which 
   653   /// each node has at most one incident edge.
   654   /// It can be formulated with the following linear program.
   655   /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
   656   /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
   657       \quad \forall B\in\mathcal{O}\f] */
   658   /// \f[x_e \ge 0\quad \forall e\in E\f]
   659   /// \f[\max \sum_{e\in E}x_ew_e\f]
   660   /// where \f$\delta(X)\f$ is the set of edges incident to a node in
   661   /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
   662   /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
   663   /// subsets of the nodes.
   664   ///
   665   /// The algorithm calculates an optimal matching and a proof of the
   666   /// optimality. The solution of the dual problem can be used to check
   667   /// the result of the algorithm. The dual linear problem is the
   668   /// following.
   669   /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
   670       z_B \ge w_{uv} \quad \forall uv\in E\f] */
   671   /// \f[y_u \ge 0 \quad \forall u \in V\f]
   672   /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
   673   /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
   674       \frac{\vert B \vert - 1}{2}z_B\f] */
   675   ///
   676   /// The algorithm can be executed with the run() function. 
   677   /// After it the matching (the primal solution) and the dual solution
   678   /// can be obtained using the query functions and the 
   679   /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, 
   680   /// which is able to iterate on the nodes of a blossom. 
   681   /// If the value type is integer, then the dual solution is multiplied
   682   /// by \ref MaxWeightedMatching::dualScale "4".
   683   ///
   684   /// \tparam GR The undirected graph type the algorithm runs on.
   685   /// \tparam WM The type edge weight map. The default type is 
   686   /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
   687 #ifdef DOXYGEN
   688   template <typename GR, typename WM>
   689 #else
   690   template <typename GR,
   691             typename WM = typename GR::template EdgeMap<int> >
   692 #endif
   693   class MaxWeightedMatching {
   694   public:
   695 
   696     /// The graph type of the algorithm
   697     typedef GR Graph;
   698     /// The type of the edge weight map
   699     typedef WM WeightMap;
   700     /// The value type of the edge weights
   701     typedef typename WeightMap::Value Value;
   702 
   703     /// The type of the matching map
   704     typedef typename Graph::template NodeMap<typename Graph::Arc>
   705     MatchingMap;
   706 
   707     /// \brief Scaling factor for dual solution
   708     ///
   709     /// Scaling factor for dual solution. It is equal to 4 or 1
   710     /// according to the value type.
   711     static const int dualScale =
   712       std::numeric_limits<Value>::is_integer ? 4 : 1;
   713 
   714   private:
   715 
   716     TEMPLATE_GRAPH_TYPEDEFS(Graph);
   717 
   718     typedef typename Graph::template NodeMap<Value> NodePotential;
   719     typedef std::vector<Node> BlossomNodeList;
   720 
   721     struct BlossomVariable {
   722       int begin, end;
   723       Value value;
   724 
   725       BlossomVariable(int _begin, int _end, Value _value)
   726         : begin(_begin), end(_end), value(_value) {}
   727 
   728     };
   729 
   730     typedef std::vector<BlossomVariable> BlossomPotential;
   731 
   732     const Graph& _graph;
   733     const WeightMap& _weight;
   734 
   735     MatchingMap* _matching;
   736 
   737     NodePotential* _node_potential;
   738 
   739     BlossomPotential _blossom_potential;
   740     BlossomNodeList _blossom_node_list;
   741 
   742     int _node_num;
   743     int _blossom_num;
   744 
   745     typedef RangeMap<int> IntIntMap;
   746 
   747     enum Status {
   748       EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
   749     };
   750 
   751     typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
   752     struct BlossomData {
   753       int tree;
   754       Status status;
   755       Arc pred, next;
   756       Value pot, offset;
   757       Node base;
   758     };
   759 
   760     IntNodeMap *_blossom_index;
   761     BlossomSet *_blossom_set;
   762     RangeMap<BlossomData>* _blossom_data;
   763 
   764     IntNodeMap *_node_index;
   765     IntArcMap *_node_heap_index;
   766 
   767     struct NodeData {
   768 
   769       NodeData(IntArcMap& node_heap_index)
   770         : heap(node_heap_index) {}
   771 
   772       int blossom;
   773       Value pot;
   774       BinHeap<Value, IntArcMap> heap;
   775       std::map<int, Arc> heap_index;
   776 
   777       int tree;
   778     };
   779 
   780     RangeMap<NodeData>* _node_data;
   781 
   782     typedef ExtendFindEnum<IntIntMap> TreeSet;
   783 
   784     IntIntMap *_tree_set_index;
   785     TreeSet *_tree_set;
   786 
   787     IntNodeMap *_delta1_index;
   788     BinHeap<Value, IntNodeMap> *_delta1;
   789 
   790     IntIntMap *_delta2_index;
   791     BinHeap<Value, IntIntMap> *_delta2;
   792 
   793     IntEdgeMap *_delta3_index;
   794     BinHeap<Value, IntEdgeMap> *_delta3;
   795 
   796     IntIntMap *_delta4_index;
   797     BinHeap<Value, IntIntMap> *_delta4;
   798 
   799     Value _delta_sum;
   800 
   801     void createStructures() {
   802       _node_num = countNodes(_graph);
   803       _blossom_num = _node_num * 3 / 2;
   804 
   805       if (!_matching) {
   806         _matching = new MatchingMap(_graph);
   807       }
   808 
   809       if (!_node_potential) {
   810         _node_potential = new NodePotential(_graph);
   811       }
   812 
   813       if (!_blossom_set) {
   814         _blossom_index = new IntNodeMap(_graph);
   815         _blossom_set = new BlossomSet(*_blossom_index);
   816         _blossom_data = new RangeMap<BlossomData>(_blossom_num);
   817       } else if (_blossom_data->size() != _blossom_num) {
   818         delete _blossom_data;
   819         _blossom_data = new RangeMap<BlossomData>(_blossom_num);
   820       }
   821 
   822       if (!_node_index) {
   823         _node_index = new IntNodeMap(_graph);
   824         _node_heap_index = new IntArcMap(_graph);
   825         _node_data = new RangeMap<NodeData>(_node_num,
   826                                             NodeData(*_node_heap_index));
   827       } else {
   828         delete _node_data;
   829         _node_data = new RangeMap<NodeData>(_node_num,
   830                                             NodeData(*_node_heap_index));
   831       }
   832 
   833       if (!_tree_set) {
   834         _tree_set_index = new IntIntMap(_blossom_num);
   835         _tree_set = new TreeSet(*_tree_set_index);
   836       } else {
   837         _tree_set_index->resize(_blossom_num);
   838       }
   839 
   840       if (!_delta1) {
   841         _delta1_index = new IntNodeMap(_graph);
   842         _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
   843       }
   844 
   845       if (!_delta2) {
   846         _delta2_index = new IntIntMap(_blossom_num);
   847         _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
   848       } else {
   849         _delta2_index->resize(_blossom_num);
   850       }
   851 
   852       if (!_delta3) {
   853         _delta3_index = new IntEdgeMap(_graph);
   854         _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
   855       }
   856 
   857       if (!_delta4) {
   858         _delta4_index = new IntIntMap(_blossom_num);
   859         _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
   860       } else {
   861         _delta4_index->resize(_blossom_num);
   862       }
   863     }
   864 
   865     void destroyStructures() {
   866       _node_num = countNodes(_graph);
   867       _blossom_num = _node_num * 3 / 2;
   868 
   869       if (_matching) {
   870         delete _matching;
   871       }
   872       if (_node_potential) {
   873         delete _node_potential;
   874       }
   875       if (_blossom_set) {
   876         delete _blossom_index;
   877         delete _blossom_set;
   878         delete _blossom_data;
   879       }
   880 
   881       if (_node_index) {
   882         delete _node_index;
   883         delete _node_heap_index;
   884         delete _node_data;
   885       }
   886 
   887       if (_tree_set) {
   888         delete _tree_set_index;
   889         delete _tree_set;
   890       }
   891       if (_delta1) {
   892         delete _delta1_index;
   893         delete _delta1;
   894       }
   895       if (_delta2) {
   896         delete _delta2_index;
   897         delete _delta2;
   898       }
   899       if (_delta3) {
   900         delete _delta3_index;
   901         delete _delta3;
   902       }
   903       if (_delta4) {
   904         delete _delta4_index;
   905         delete _delta4;
   906       }
   907     }
   908 
   909     void matchedToEven(int blossom, int tree) {
   910       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
   911         _delta2->erase(blossom);
   912       }
   913 
   914       if (!_blossom_set->trivial(blossom)) {
   915         (*_blossom_data)[blossom].pot -=
   916           2 * (_delta_sum - (*_blossom_data)[blossom].offset);
   917       }
   918 
   919       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
   920            n != INVALID; ++n) {
   921 
   922         _blossom_set->increase(n, std::numeric_limits<Value>::max());
   923         int ni = (*_node_index)[n];
   924 
   925         (*_node_data)[ni].heap.clear();
   926         (*_node_data)[ni].heap_index.clear();
   927 
   928         (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
   929 
   930         _delta1->push(n, (*_node_data)[ni].pot);
   931 
   932         for (InArcIt e(_graph, n); e != INVALID; ++e) {
   933           Node v = _graph.source(e);
   934           int vb = _blossom_set->find(v);
   935           int vi = (*_node_index)[v];
   936 
   937           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   938             dualScale * _weight[e];
   939 
   940           if ((*_blossom_data)[vb].status == EVEN) {
   941             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
   942               _delta3->push(e, rw / 2);
   943             }
   944           } else if ((*_blossom_data)[vb].status == UNMATCHED) {
   945             if (_delta3->state(e) != _delta3->IN_HEAP) {
   946               _delta3->push(e, rw);
   947             }
   948           } else {
   949             typename std::map<int, Arc>::iterator it =
   950               (*_node_data)[vi].heap_index.find(tree);
   951 
   952             if (it != (*_node_data)[vi].heap_index.end()) {
   953               if ((*_node_data)[vi].heap[it->second] > rw) {
   954                 (*_node_data)[vi].heap.replace(it->second, e);
   955                 (*_node_data)[vi].heap.decrease(e, rw);
   956                 it->second = e;
   957               }
   958             } else {
   959               (*_node_data)[vi].heap.push(e, rw);
   960               (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
   961             }
   962 
   963             if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
   964               _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
   965 
   966               if ((*_blossom_data)[vb].status == MATCHED) {
   967                 if (_delta2->state(vb) != _delta2->IN_HEAP) {
   968                   _delta2->push(vb, _blossom_set->classPrio(vb) -
   969                                (*_blossom_data)[vb].offset);
   970                 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
   971                            (*_blossom_data)[vb].offset){
   972                   _delta2->decrease(vb, _blossom_set->classPrio(vb) -
   973                                    (*_blossom_data)[vb].offset);
   974                 }
   975               }
   976             }
   977           }
   978         }
   979       }
   980       (*_blossom_data)[blossom].offset = 0;
   981     }
   982 
   983     void matchedToOdd(int blossom) {
   984       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
   985         _delta2->erase(blossom);
   986       }
   987       (*_blossom_data)[blossom].offset += _delta_sum;
   988       if (!_blossom_set->trivial(blossom)) {
   989         _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
   990                      (*_blossom_data)[blossom].offset);
   991       }
   992     }
   993 
   994     void evenToMatched(int blossom, int tree) {
   995       if (!_blossom_set->trivial(blossom)) {
   996         (*_blossom_data)[blossom].pot += 2 * _delta_sum;
   997       }
   998 
   999       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  1000            n != INVALID; ++n) {
  1001         int ni = (*_node_index)[n];
  1002         (*_node_data)[ni].pot -= _delta_sum;
  1003 
  1004         _delta1->erase(n);
  1005 
  1006         for (InArcIt e(_graph, n); e != INVALID; ++e) {
  1007           Node v = _graph.source(e);
  1008           int vb = _blossom_set->find(v);
  1009           int vi = (*_node_index)[v];
  1010 
  1011           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  1012             dualScale * _weight[e];
  1013 
  1014           if (vb == blossom) {
  1015             if (_delta3->state(e) == _delta3->IN_HEAP) {
  1016               _delta3->erase(e);
  1017             }
  1018           } else if ((*_blossom_data)[vb].status == EVEN) {
  1019 
  1020             if (_delta3->state(e) == _delta3->IN_HEAP) {
  1021               _delta3->erase(e);
  1022             }
  1023 
  1024             int vt = _tree_set->find(vb);
  1025 
  1026             if (vt != tree) {
  1027 
  1028               Arc r = _graph.oppositeArc(e);
  1029 
  1030               typename std::map<int, Arc>::iterator it =
  1031                 (*_node_data)[ni].heap_index.find(vt);
  1032 
  1033               if (it != (*_node_data)[ni].heap_index.end()) {
  1034                 if ((*_node_data)[ni].heap[it->second] > rw) {
  1035                   (*_node_data)[ni].heap.replace(it->second, r);
  1036                   (*_node_data)[ni].heap.decrease(r, rw);
  1037                   it->second = r;
  1038                 }
  1039               } else {
  1040                 (*_node_data)[ni].heap.push(r, rw);
  1041                 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
  1042               }
  1043 
  1044               if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
  1045                 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
  1046 
  1047                 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
  1048                   _delta2->push(blossom, _blossom_set->classPrio(blossom) -
  1049                                (*_blossom_data)[blossom].offset);
  1050                 } else if ((*_delta2)[blossom] >
  1051                            _blossom_set->classPrio(blossom) -
  1052                            (*_blossom_data)[blossom].offset){
  1053                   _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
  1054                                    (*_blossom_data)[blossom].offset);
  1055                 }
  1056               }
  1057             }
  1058 
  1059           } else if ((*_blossom_data)[vb].status == UNMATCHED) {
  1060             if (_delta3->state(e) == _delta3->IN_HEAP) {
  1061               _delta3->erase(e);
  1062             }
  1063           } else {
  1064 
  1065             typename std::map<int, Arc>::iterator it =
  1066               (*_node_data)[vi].heap_index.find(tree);
  1067 
  1068             if (it != (*_node_data)[vi].heap_index.end()) {
  1069               (*_node_data)[vi].heap.erase(it->second);
  1070               (*_node_data)[vi].heap_index.erase(it);
  1071               if ((*_node_data)[vi].heap.empty()) {
  1072                 _blossom_set->increase(v, std::numeric_limits<Value>::max());
  1073               } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
  1074                 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
  1075               }
  1076 
  1077               if ((*_blossom_data)[vb].status == MATCHED) {
  1078                 if (_blossom_set->classPrio(vb) ==
  1079                     std::numeric_limits<Value>::max()) {
  1080                   _delta2->erase(vb);
  1081                 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
  1082                            (*_blossom_data)[vb].offset) {
  1083                   _delta2->increase(vb, _blossom_set->classPrio(vb) -
  1084                                    (*_blossom_data)[vb].offset);
  1085                 }
  1086               }
  1087             }
  1088           }
  1089         }
  1090       }
  1091     }
  1092 
  1093     void oddToMatched(int blossom) {
  1094       (*_blossom_data)[blossom].offset -= _delta_sum;
  1095 
  1096       if (_blossom_set->classPrio(blossom) !=
  1097           std::numeric_limits<Value>::max()) {
  1098         _delta2->push(blossom, _blossom_set->classPrio(blossom) -
  1099                        (*_blossom_data)[blossom].offset);
  1100       }
  1101 
  1102       if (!_blossom_set->trivial(blossom)) {
  1103         _delta4->erase(blossom);
  1104       }
  1105     }
  1106 
  1107     void oddToEven(int blossom, int tree) {
  1108       if (!_blossom_set->trivial(blossom)) {
  1109         _delta4->erase(blossom);
  1110         (*_blossom_data)[blossom].pot -=
  1111           2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
  1112       }
  1113 
  1114       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  1115            n != INVALID; ++n) {
  1116         int ni = (*_node_index)[n];
  1117 
  1118         _blossom_set->increase(n, std::numeric_limits<Value>::max());
  1119 
  1120         (*_node_data)[ni].heap.clear();
  1121         (*_node_data)[ni].heap_index.clear();
  1122         (*_node_data)[ni].pot +=
  1123           2 * _delta_sum - (*_blossom_data)[blossom].offset;
  1124 
  1125         _delta1->push(n, (*_node_data)[ni].pot);
  1126 
  1127         for (InArcIt e(_graph, n); e != INVALID; ++e) {
  1128           Node v = _graph.source(e);
  1129           int vb = _blossom_set->find(v);
  1130           int vi = (*_node_index)[v];
  1131 
  1132           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  1133             dualScale * _weight[e];
  1134 
  1135           if ((*_blossom_data)[vb].status == EVEN) {
  1136             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
  1137               _delta3->push(e, rw / 2);
  1138             }
  1139           } else if ((*_blossom_data)[vb].status == UNMATCHED) {
  1140             if (_delta3->state(e) != _delta3->IN_HEAP) {
  1141               _delta3->push(e, rw);
  1142             }
  1143           } else {
  1144 
  1145             typename std::map<int, Arc>::iterator it =
  1146               (*_node_data)[vi].heap_index.find(tree);
  1147 
  1148             if (it != (*_node_data)[vi].heap_index.end()) {
  1149               if ((*_node_data)[vi].heap[it->second] > rw) {
  1150                 (*_node_data)[vi].heap.replace(it->second, e);
  1151                 (*_node_data)[vi].heap.decrease(e, rw);
  1152                 it->second = e;
  1153               }
  1154             } else {
  1155               (*_node_data)[vi].heap.push(e, rw);
  1156               (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
  1157             }
  1158 
  1159             if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
  1160               _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
  1161 
  1162               if ((*_blossom_data)[vb].status == MATCHED) {
  1163                 if (_delta2->state(vb) != _delta2->IN_HEAP) {
  1164                   _delta2->push(vb, _blossom_set->classPrio(vb) -
  1165                                (*_blossom_data)[vb].offset);
  1166                 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
  1167                            (*_blossom_data)[vb].offset) {
  1168                   _delta2->decrease(vb, _blossom_set->classPrio(vb) -
  1169                                    (*_blossom_data)[vb].offset);
  1170                 }
  1171               }
  1172             }
  1173           }
  1174         }
  1175       }
  1176       (*_blossom_data)[blossom].offset = 0;
  1177     }
  1178 
  1179 
  1180     void matchedToUnmatched(int blossom) {
  1181       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
  1182         _delta2->erase(blossom);
  1183       }
  1184 
  1185       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  1186            n != INVALID; ++n) {
  1187         int ni = (*_node_index)[n];
  1188 
  1189         _blossom_set->increase(n, std::numeric_limits<Value>::max());
  1190 
  1191         (*_node_data)[ni].heap.clear();
  1192         (*_node_data)[ni].heap_index.clear();
  1193 
  1194         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  1195           Node v = _graph.target(e);
  1196           int vb = _blossom_set->find(v);
  1197           int vi = (*_node_index)[v];
  1198 
  1199           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  1200             dualScale * _weight[e];
  1201 
  1202           if ((*_blossom_data)[vb].status == EVEN) {
  1203             if (_delta3->state(e) != _delta3->IN_HEAP) {
  1204               _delta3->push(e, rw);
  1205             }
  1206           }
  1207         }
  1208       }
  1209     }
  1210 
  1211     void unmatchedToMatched(int blossom) {
  1212       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  1213            n != INVALID; ++n) {
  1214         int ni = (*_node_index)[n];
  1215 
  1216         for (InArcIt e(_graph, n); e != INVALID; ++e) {
  1217           Node v = _graph.source(e);
  1218           int vb = _blossom_set->find(v);
  1219           int vi = (*_node_index)[v];
  1220 
  1221           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  1222             dualScale * _weight[e];
  1223 
  1224           if (vb == blossom) {
  1225             if (_delta3->state(e) == _delta3->IN_HEAP) {
  1226               _delta3->erase(e);
  1227             }
  1228           } else if ((*_blossom_data)[vb].status == EVEN) {
  1229 
  1230             if (_delta3->state(e) == _delta3->IN_HEAP) {
  1231               _delta3->erase(e);
  1232             }
  1233 
  1234             int vt = _tree_set->find(vb);
  1235 
  1236             Arc r = _graph.oppositeArc(e);
  1237 
  1238             typename std::map<int, Arc>::iterator it =
  1239               (*_node_data)[ni].heap_index.find(vt);
  1240 
  1241             if (it != (*_node_data)[ni].heap_index.end()) {
  1242               if ((*_node_data)[ni].heap[it->second] > rw) {
  1243                 (*_node_data)[ni].heap.replace(it->second, r);
  1244                 (*_node_data)[ni].heap.decrease(r, rw);
  1245                 it->second = r;
  1246               }
  1247             } else {
  1248               (*_node_data)[ni].heap.push(r, rw);
  1249               (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
  1250             }
  1251 
  1252             if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
  1253               _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
  1254 
  1255               if (_delta2->state(blossom) != _delta2->IN_HEAP) {
  1256                 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
  1257                              (*_blossom_data)[blossom].offset);
  1258               } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
  1259                          (*_blossom_data)[blossom].offset){
  1260                 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
  1261                                  (*_blossom_data)[blossom].offset);
  1262               }
  1263             }
  1264 
  1265           } else if ((*_blossom_data)[vb].status == UNMATCHED) {
  1266             if (_delta3->state(e) == _delta3->IN_HEAP) {
  1267               _delta3->erase(e);
  1268             }
  1269           }
  1270         }
  1271       }
  1272     }
  1273 
  1274     void alternatePath(int even, int tree) {
  1275       int odd;
  1276 
  1277       evenToMatched(even, tree);
  1278       (*_blossom_data)[even].status = MATCHED;
  1279 
  1280       while ((*_blossom_data)[even].pred != INVALID) {
  1281         odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
  1282         (*_blossom_data)[odd].status = MATCHED;
  1283         oddToMatched(odd);
  1284         (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
  1285 
  1286         even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
  1287         (*_blossom_data)[even].status = MATCHED;
  1288         evenToMatched(even, tree);
  1289         (*_blossom_data)[even].next =
  1290           _graph.oppositeArc((*_blossom_data)[odd].pred);
  1291       }
  1292 
  1293     }
  1294 
  1295     void destroyTree(int tree) {
  1296       for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
  1297         if ((*_blossom_data)[b].status == EVEN) {
  1298           (*_blossom_data)[b].status = MATCHED;
  1299           evenToMatched(b, tree);
  1300         } else if ((*_blossom_data)[b].status == ODD) {
  1301           (*_blossom_data)[b].status = MATCHED;
  1302           oddToMatched(b);
  1303         }
  1304       }
  1305       _tree_set->eraseClass(tree);
  1306     }
  1307 
  1308 
  1309     void unmatchNode(const Node& node) {
  1310       int blossom = _blossom_set->find(node);
  1311       int tree = _tree_set->find(blossom);
  1312 
  1313       alternatePath(blossom, tree);
  1314       destroyTree(tree);
  1315 
  1316       (*_blossom_data)[blossom].status = UNMATCHED;
  1317       (*_blossom_data)[blossom].base = node;
  1318       matchedToUnmatched(blossom);
  1319     }
  1320 
  1321 
  1322     void augmentOnEdge(const Edge& edge) {
  1323 
  1324       int left = _blossom_set->find(_graph.u(edge));
  1325       int right = _blossom_set->find(_graph.v(edge));
  1326 
  1327       if ((*_blossom_data)[left].status == EVEN) {
  1328         int left_tree = _tree_set->find(left);
  1329         alternatePath(left, left_tree);
  1330         destroyTree(left_tree);
  1331       } else {
  1332         (*_blossom_data)[left].status = MATCHED;
  1333         unmatchedToMatched(left);
  1334       }
  1335 
  1336       if ((*_blossom_data)[right].status == EVEN) {
  1337         int right_tree = _tree_set->find(right);
  1338         alternatePath(right, right_tree);
  1339         destroyTree(right_tree);
  1340       } else {
  1341         (*_blossom_data)[right].status = MATCHED;
  1342         unmatchedToMatched(right);
  1343       }
  1344 
  1345       (*_blossom_data)[left].next = _graph.direct(edge, true);
  1346       (*_blossom_data)[right].next = _graph.direct(edge, false);
  1347     }
  1348 
  1349     void extendOnArc(const Arc& arc) {
  1350       int base = _blossom_set->find(_graph.target(arc));
  1351       int tree = _tree_set->find(base);
  1352 
  1353       int odd = _blossom_set->find(_graph.source(arc));
  1354       _tree_set->insert(odd, tree);
  1355       (*_blossom_data)[odd].status = ODD;
  1356       matchedToOdd(odd);
  1357       (*_blossom_data)[odd].pred = arc;
  1358 
  1359       int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
  1360       (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
  1361       _tree_set->insert(even, tree);
  1362       (*_blossom_data)[even].status = EVEN;
  1363       matchedToEven(even, tree);
  1364     }
  1365 
  1366     void shrinkOnEdge(const Edge& edge, int tree) {
  1367       int nca = -1;
  1368       std::vector<int> left_path, right_path;
  1369 
  1370       {
  1371         std::set<int> left_set, right_set;
  1372         int left = _blossom_set->find(_graph.u(edge));
  1373         left_path.push_back(left);
  1374         left_set.insert(left);
  1375 
  1376         int right = _blossom_set->find(_graph.v(edge));
  1377         right_path.push_back(right);
  1378         right_set.insert(right);
  1379 
  1380         while (true) {
  1381 
  1382           if ((*_blossom_data)[left].pred == INVALID) break;
  1383 
  1384           left =
  1385             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
  1386           left_path.push_back(left);
  1387           left =
  1388             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
  1389           left_path.push_back(left);
  1390 
  1391           left_set.insert(left);
  1392 
  1393           if (right_set.find(left) != right_set.end()) {
  1394             nca = left;
  1395             break;
  1396           }
  1397 
  1398           if ((*_blossom_data)[right].pred == INVALID) break;
  1399 
  1400           right =
  1401             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
  1402           right_path.push_back(right);
  1403           right =
  1404             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
  1405           right_path.push_back(right);
  1406 
  1407           right_set.insert(right);
  1408 
  1409           if (left_set.find(right) != left_set.end()) {
  1410             nca = right;
  1411             break;
  1412           }
  1413 
  1414         }
  1415 
  1416         if (nca == -1) {
  1417           if ((*_blossom_data)[left].pred == INVALID) {
  1418             nca = right;
  1419             while (left_set.find(nca) == left_set.end()) {
  1420               nca =
  1421                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  1422               right_path.push_back(nca);
  1423               nca =
  1424                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  1425               right_path.push_back(nca);
  1426             }
  1427           } else {
  1428             nca = left;
  1429             while (right_set.find(nca) == right_set.end()) {
  1430               nca =
  1431                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  1432               left_path.push_back(nca);
  1433               nca =
  1434                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  1435               left_path.push_back(nca);
  1436             }
  1437           }
  1438         }
  1439       }
  1440 
  1441       std::vector<int> subblossoms;
  1442       Arc prev;
  1443 
  1444       prev = _graph.direct(edge, true);
  1445       for (int i = 0; left_path[i] != nca; i += 2) {
  1446         subblossoms.push_back(left_path[i]);
  1447         (*_blossom_data)[left_path[i]].next = prev;
  1448         _tree_set->erase(left_path[i]);
  1449 
  1450         subblossoms.push_back(left_path[i + 1]);
  1451         (*_blossom_data)[left_path[i + 1]].status = EVEN;
  1452         oddToEven(left_path[i + 1], tree);
  1453         _tree_set->erase(left_path[i + 1]);
  1454         prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
  1455       }
  1456 
  1457       int k = 0;
  1458       while (right_path[k] != nca) ++k;
  1459 
  1460       subblossoms.push_back(nca);
  1461       (*_blossom_data)[nca].next = prev;
  1462 
  1463       for (int i = k - 2; i >= 0; i -= 2) {
  1464         subblossoms.push_back(right_path[i + 1]);
  1465         (*_blossom_data)[right_path[i + 1]].status = EVEN;
  1466         oddToEven(right_path[i + 1], tree);
  1467         _tree_set->erase(right_path[i + 1]);
  1468 
  1469         (*_blossom_data)[right_path[i + 1]].next =
  1470           (*_blossom_data)[right_path[i + 1]].pred;
  1471 
  1472         subblossoms.push_back(right_path[i]);
  1473         _tree_set->erase(right_path[i]);
  1474       }
  1475 
  1476       int surface =
  1477         _blossom_set->join(subblossoms.begin(), subblossoms.end());
  1478 
  1479       for (int i = 0; i < int(subblossoms.size()); ++i) {
  1480         if (!_blossom_set->trivial(subblossoms[i])) {
  1481           (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
  1482         }
  1483         (*_blossom_data)[subblossoms[i]].status = MATCHED;
  1484       }
  1485 
  1486       (*_blossom_data)[surface].pot = -2 * _delta_sum;
  1487       (*_blossom_data)[surface].offset = 0;
  1488       (*_blossom_data)[surface].status = EVEN;
  1489       (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
  1490       (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
  1491 
  1492       _tree_set->insert(surface, tree);
  1493       _tree_set->erase(nca);
  1494     }
  1495 
  1496     void splitBlossom(int blossom) {
  1497       Arc next = (*_blossom_data)[blossom].next;
  1498       Arc pred = (*_blossom_data)[blossom].pred;
  1499 
  1500       int tree = _tree_set->find(blossom);
  1501 
  1502       (*_blossom_data)[blossom].status = MATCHED;
  1503       oddToMatched(blossom);
  1504       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
  1505         _delta2->erase(blossom);
  1506       }
  1507 
  1508       std::vector<int> subblossoms;
  1509       _blossom_set->split(blossom, std::back_inserter(subblossoms));
  1510 
  1511       Value offset = (*_blossom_data)[blossom].offset;
  1512       int b = _blossom_set->find(_graph.source(pred));
  1513       int d = _blossom_set->find(_graph.source(next));
  1514 
  1515       int ib = -1, id = -1;
  1516       for (int i = 0; i < int(subblossoms.size()); ++i) {
  1517         if (subblossoms[i] == b) ib = i;
  1518         if (subblossoms[i] == d) id = i;
  1519 
  1520         (*_blossom_data)[subblossoms[i]].offset = offset;
  1521         if (!_blossom_set->trivial(subblossoms[i])) {
  1522           (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
  1523         }
  1524         if (_blossom_set->classPrio(subblossoms[i]) !=
  1525             std::numeric_limits<Value>::max()) {
  1526           _delta2->push(subblossoms[i],
  1527                         _blossom_set->classPrio(subblossoms[i]) -
  1528                         (*_blossom_data)[subblossoms[i]].offset);
  1529         }
  1530       }
  1531 
  1532       if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
  1533         for (int i = (id + 1) % subblossoms.size();
  1534              i != ib; i = (i + 2) % subblossoms.size()) {
  1535           int sb = subblossoms[i];
  1536           int tb = subblossoms[(i + 1) % subblossoms.size()];
  1537           (*_blossom_data)[sb].next =
  1538             _graph.oppositeArc((*_blossom_data)[tb].next);
  1539         }
  1540 
  1541         for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
  1542           int sb = subblossoms[i];
  1543           int tb = subblossoms[(i + 1) % subblossoms.size()];
  1544           int ub = subblossoms[(i + 2) % subblossoms.size()];
  1545 
  1546           (*_blossom_data)[sb].status = ODD;
  1547           matchedToOdd(sb);
  1548           _tree_set->insert(sb, tree);
  1549           (*_blossom_data)[sb].pred = pred;
  1550           (*_blossom_data)[sb].next =
  1551                            _graph.oppositeArc((*_blossom_data)[tb].next);
  1552 
  1553           pred = (*_blossom_data)[ub].next;
  1554 
  1555           (*_blossom_data)[tb].status = EVEN;
  1556           matchedToEven(tb, tree);
  1557           _tree_set->insert(tb, tree);
  1558           (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
  1559         }
  1560 
  1561         (*_blossom_data)[subblossoms[id]].status = ODD;
  1562         matchedToOdd(subblossoms[id]);
  1563         _tree_set->insert(subblossoms[id], tree);
  1564         (*_blossom_data)[subblossoms[id]].next = next;
  1565         (*_blossom_data)[subblossoms[id]].pred = pred;
  1566 
  1567       } else {
  1568 
  1569         for (int i = (ib + 1) % subblossoms.size();
  1570              i != id; i = (i + 2) % subblossoms.size()) {
  1571           int sb = subblossoms[i];
  1572           int tb = subblossoms[(i + 1) % subblossoms.size()];
  1573           (*_blossom_data)[sb].next =
  1574             _graph.oppositeArc((*_blossom_data)[tb].next);
  1575         }
  1576 
  1577         for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
  1578           int sb = subblossoms[i];
  1579           int tb = subblossoms[(i + 1) % subblossoms.size()];
  1580           int ub = subblossoms[(i + 2) % subblossoms.size()];
  1581 
  1582           (*_blossom_data)[sb].status = ODD;
  1583           matchedToOdd(sb);
  1584           _tree_set->insert(sb, tree);
  1585           (*_blossom_data)[sb].next = next;
  1586           (*_blossom_data)[sb].pred =
  1587             _graph.oppositeArc((*_blossom_data)[tb].next);
  1588 
  1589           (*_blossom_data)[tb].status = EVEN;
  1590           matchedToEven(tb, tree);
  1591           _tree_set->insert(tb, tree);
  1592           (*_blossom_data)[tb].pred =
  1593             (*_blossom_data)[tb].next =
  1594             _graph.oppositeArc((*_blossom_data)[ub].next);
  1595           next = (*_blossom_data)[ub].next;
  1596         }
  1597 
  1598         (*_blossom_data)[subblossoms[ib]].status = ODD;
  1599         matchedToOdd(subblossoms[ib]);
  1600         _tree_set->insert(subblossoms[ib], tree);
  1601         (*_blossom_data)[subblossoms[ib]].next = next;
  1602         (*_blossom_data)[subblossoms[ib]].pred = pred;
  1603       }
  1604       _tree_set->erase(blossom);
  1605     }
  1606 
  1607     void extractBlossom(int blossom, const Node& base, const Arc& matching) {
  1608       if (_blossom_set->trivial(blossom)) {
  1609         int bi = (*_node_index)[base];
  1610         Value pot = (*_node_data)[bi].pot;
  1611 
  1612         (*_matching)[base] = matching;
  1613         _blossom_node_list.push_back(base);
  1614         (*_node_potential)[base] = pot;
  1615       } else {
  1616 
  1617         Value pot = (*_blossom_data)[blossom].pot;
  1618         int bn = _blossom_node_list.size();
  1619 
  1620         std::vector<int> subblossoms;
  1621         _blossom_set->split(blossom, std::back_inserter(subblossoms));
  1622         int b = _blossom_set->find(base);
  1623         int ib = -1;
  1624         for (int i = 0; i < int(subblossoms.size()); ++i) {
  1625           if (subblossoms[i] == b) { ib = i; break; }
  1626         }
  1627 
  1628         for (int i = 1; i < int(subblossoms.size()); i += 2) {
  1629           int sb = subblossoms[(ib + i) % subblossoms.size()];
  1630           int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
  1631 
  1632           Arc m = (*_blossom_data)[tb].next;
  1633           extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
  1634           extractBlossom(tb, _graph.source(m), m);
  1635         }
  1636         extractBlossom(subblossoms[ib], base, matching);
  1637 
  1638         int en = _blossom_node_list.size();
  1639 
  1640         _blossom_potential.push_back(BlossomVariable(bn, en, pot));
  1641       }
  1642     }
  1643 
  1644     void extractMatching() {
  1645       std::vector<int> blossoms;
  1646       for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
  1647         blossoms.push_back(c);
  1648       }
  1649 
  1650       for (int i = 0; i < int(blossoms.size()); ++i) {
  1651         if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
  1652 
  1653           Value offset = (*_blossom_data)[blossoms[i]].offset;
  1654           (*_blossom_data)[blossoms[i]].pot += 2 * offset;
  1655           for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
  1656                n != INVALID; ++n) {
  1657             (*_node_data)[(*_node_index)[n]].pot -= offset;
  1658           }
  1659 
  1660           Arc matching = (*_blossom_data)[blossoms[i]].next;
  1661           Node base = _graph.source(matching);
  1662           extractBlossom(blossoms[i], base, matching);
  1663         } else {
  1664           Node base = (*_blossom_data)[blossoms[i]].base;
  1665           extractBlossom(blossoms[i], base, INVALID);
  1666         }
  1667       }
  1668     }
  1669 
  1670   public:
  1671 
  1672     /// \brief Constructor
  1673     ///
  1674     /// Constructor.
  1675     MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
  1676       : _graph(graph), _weight(weight), _matching(0),
  1677         _node_potential(0), _blossom_potential(), _blossom_node_list(),
  1678         _node_num(0), _blossom_num(0),
  1679 
  1680         _blossom_index(0), _blossom_set(0), _blossom_data(0),
  1681         _node_index(0), _node_heap_index(0), _node_data(0),
  1682         _tree_set_index(0), _tree_set(0),
  1683 
  1684         _delta1_index(0), _delta1(0),
  1685         _delta2_index(0), _delta2(0),
  1686         _delta3_index(0), _delta3(0),
  1687         _delta4_index(0), _delta4(0),
  1688 
  1689         _delta_sum() {}
  1690 
  1691     ~MaxWeightedMatching() {
  1692       destroyStructures();
  1693     }
  1694 
  1695     /// \name Execution Control
  1696     /// The simplest way to execute the algorithm is to use the
  1697     /// \ref run() member function.
  1698 
  1699     ///@{
  1700 
  1701     /// \brief Initialize the algorithm
  1702     ///
  1703     /// This function initializes the algorithm.
  1704     void init() {
  1705       createStructures();
  1706 
  1707       _blossom_node_list.clear();
  1708       _blossom_potential.clear();
  1709 
  1710       for (ArcIt e(_graph); e != INVALID; ++e) {
  1711         (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
  1712       }
  1713       for (NodeIt n(_graph); n != INVALID; ++n) {
  1714         (*_delta1_index)[n] = _delta1->PRE_HEAP;
  1715       }
  1716       for (EdgeIt e(_graph); e != INVALID; ++e) {
  1717         (*_delta3_index)[e] = _delta3->PRE_HEAP;
  1718       }
  1719       for (int i = 0; i < _blossom_num; ++i) {
  1720         (*_delta2_index)[i] = _delta2->PRE_HEAP;
  1721         (*_delta4_index)[i] = _delta4->PRE_HEAP;
  1722       }
  1723       
  1724       _delta1->clear();
  1725       _delta2->clear();
  1726       _delta3->clear();
  1727       _delta4->clear();
  1728       _blossom_set->clear();
  1729       _tree_set->clear();
  1730 
  1731       int index = 0;
  1732       for (NodeIt n(_graph); n != INVALID; ++n) {
  1733         Value max = 0;
  1734         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  1735           if (_graph.target(e) == n) continue;
  1736           if ((dualScale * _weight[e]) / 2 > max) {
  1737             max = (dualScale * _weight[e]) / 2;
  1738           }
  1739         }
  1740         (*_node_index)[n] = index;
  1741         (*_node_data)[index].heap_index.clear();
  1742         (*_node_data)[index].heap.clear();
  1743         (*_node_data)[index].pot = max;
  1744         _delta1->push(n, max);
  1745         int blossom =
  1746           _blossom_set->insert(n, std::numeric_limits<Value>::max());
  1747 
  1748         _tree_set->insert(blossom);
  1749 
  1750         (*_blossom_data)[blossom].status = EVEN;
  1751         (*_blossom_data)[blossom].pred = INVALID;
  1752         (*_blossom_data)[blossom].next = INVALID;
  1753         (*_blossom_data)[blossom].pot = 0;
  1754         (*_blossom_data)[blossom].offset = 0;
  1755         ++index;
  1756       }
  1757       for (EdgeIt e(_graph); e != INVALID; ++e) {
  1758         int si = (*_node_index)[_graph.u(e)];
  1759         int ti = (*_node_index)[_graph.v(e)];
  1760         if (_graph.u(e) != _graph.v(e)) {
  1761           _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
  1762                             dualScale * _weight[e]) / 2);
  1763         }
  1764       }
  1765     }
  1766 
  1767     /// \brief Start the algorithm
  1768     ///
  1769     /// This function starts the algorithm.
  1770     ///
  1771     /// \pre \ref init() must be called before using this function.
  1772     void start() {
  1773       enum OpType {
  1774         D1, D2, D3, D4
  1775       };
  1776 
  1777       int unmatched = _node_num;
  1778       while (unmatched > 0) {
  1779         Value d1 = !_delta1->empty() ?
  1780           _delta1->prio() : std::numeric_limits<Value>::max();
  1781 
  1782         Value d2 = !_delta2->empty() ?
  1783           _delta2->prio() : std::numeric_limits<Value>::max();
  1784 
  1785         Value d3 = !_delta3->empty() ?
  1786           _delta3->prio() : std::numeric_limits<Value>::max();
  1787 
  1788         Value d4 = !_delta4->empty() ?
  1789           _delta4->prio() : std::numeric_limits<Value>::max();
  1790 
  1791         _delta_sum = d1; OpType ot = D1;
  1792         if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
  1793         if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
  1794         if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
  1795 
  1796 
  1797         switch (ot) {
  1798         case D1:
  1799           {
  1800             Node n = _delta1->top();
  1801             unmatchNode(n);
  1802             --unmatched;
  1803           }
  1804           break;
  1805         case D2:
  1806           {
  1807             int blossom = _delta2->top();
  1808             Node n = _blossom_set->classTop(blossom);
  1809             Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
  1810             extendOnArc(e);
  1811           }
  1812           break;
  1813         case D3:
  1814           {
  1815             Edge e = _delta3->top();
  1816 
  1817             int left_blossom = _blossom_set->find(_graph.u(e));
  1818             int right_blossom = _blossom_set->find(_graph.v(e));
  1819 
  1820             if (left_blossom == right_blossom) {
  1821               _delta3->pop();
  1822             } else {
  1823               int left_tree;
  1824               if ((*_blossom_data)[left_blossom].status == EVEN) {
  1825                 left_tree = _tree_set->find(left_blossom);
  1826               } else {
  1827                 left_tree = -1;
  1828                 ++unmatched;
  1829               }
  1830               int right_tree;
  1831               if ((*_blossom_data)[right_blossom].status == EVEN) {
  1832                 right_tree = _tree_set->find(right_blossom);
  1833               } else {
  1834                 right_tree = -1;
  1835                 ++unmatched;
  1836               }
  1837 
  1838               if (left_tree == right_tree) {
  1839                 shrinkOnEdge(e, left_tree);
  1840               } else {
  1841                 augmentOnEdge(e);
  1842                 unmatched -= 2;
  1843               }
  1844             }
  1845           } break;
  1846         case D4:
  1847           splitBlossom(_delta4->top());
  1848           break;
  1849         }
  1850       }
  1851       extractMatching();
  1852     }
  1853 
  1854     /// \brief Run the algorithm.
  1855     ///
  1856     /// This method runs the \c %MaxWeightedMatching algorithm.
  1857     ///
  1858     /// \note mwm.run() is just a shortcut of the following code.
  1859     /// \code
  1860     ///   mwm.init();
  1861     ///   mwm.start();
  1862     /// \endcode
  1863     void run() {
  1864       init();
  1865       start();
  1866     }
  1867 
  1868     /// @}
  1869 
  1870     /// \name Primal Solution
  1871     /// Functions to get the primal solution, i.e. the maximum weighted 
  1872     /// matching.\n
  1873     /// Either \ref run() or \ref start() function should be called before
  1874     /// using them.
  1875 
  1876     /// @{
  1877 
  1878     /// \brief Return the weight of the matching.
  1879     ///
  1880     /// This function returns the weight of the found matching.
  1881     ///
  1882     /// \pre Either run() or start() must be called before using this function.
  1883     Value matchingWeight() const {
  1884       Value sum = 0;
  1885       for (NodeIt n(_graph); n != INVALID; ++n) {
  1886         if ((*_matching)[n] != INVALID) {
  1887           sum += _weight[(*_matching)[n]];
  1888         }
  1889       }
  1890       return sum /= 2;
  1891     }
  1892 
  1893     /// \brief Return the size (cardinality) of the matching.
  1894     ///
  1895     /// This function returns the size (cardinality) of the found matching.
  1896     ///
  1897     /// \pre Either run() or start() must be called before using this function.
  1898     int matchingSize() const {
  1899       int num = 0;
  1900       for (NodeIt n(_graph); n != INVALID; ++n) {
  1901         if ((*_matching)[n] != INVALID) {
  1902           ++num;
  1903         }
  1904       }
  1905       return num /= 2;
  1906     }
  1907 
  1908     /// \brief Return \c true if the given edge is in the matching.
  1909     ///
  1910     /// This function returns \c true if the given edge is in the found 
  1911     /// matching.
  1912     ///
  1913     /// \pre Either run() or start() must be called before using this function.
  1914     bool matching(const Edge& edge) const {
  1915       return edge == (*_matching)[_graph.u(edge)];
  1916     }
  1917 
  1918     /// \brief Return the matching arc (or edge) incident to the given node.
  1919     ///
  1920     /// This function returns the matching arc (or edge) incident to the
  1921     /// given node in the found matching or \c INVALID if the node is 
  1922     /// not covered by the matching.
  1923     ///
  1924     /// \pre Either run() or start() must be called before using this function.
  1925     Arc matching(const Node& node) const {
  1926       return (*_matching)[node];
  1927     }
  1928 
  1929     /// \brief Return a const reference to the matching map.
  1930     ///
  1931     /// This function returns a const reference to a node map that stores
  1932     /// the matching arc (or edge) incident to each node.
  1933     const MatchingMap& matchingMap() const {
  1934       return *_matching;
  1935     }
  1936 
  1937     /// \brief Return the mate of the given node.
  1938     ///
  1939     /// This function returns the mate of the given node in the found 
  1940     /// matching or \c INVALID if the node is not covered by the matching.
  1941     ///
  1942     /// \pre Either run() or start() must be called before using this function.
  1943     Node mate(const Node& node) const {
  1944       return (*_matching)[node] != INVALID ?
  1945         _graph.target((*_matching)[node]) : INVALID;
  1946     }
  1947 
  1948     /// @}
  1949 
  1950     /// \name Dual Solution
  1951     /// Functions to get the dual solution.\n
  1952     /// Either \ref run() or \ref start() function should be called before
  1953     /// using them.
  1954 
  1955     /// @{
  1956 
  1957     /// \brief Return the value of the dual solution.
  1958     ///
  1959     /// This function returns the value of the dual solution. 
  1960     /// It should be equal to the primal value scaled by \ref dualScale 
  1961     /// "dual scale".
  1962     ///
  1963     /// \pre Either run() or start() must be called before using this function.
  1964     Value dualValue() const {
  1965       Value sum = 0;
  1966       for (NodeIt n(_graph); n != INVALID; ++n) {
  1967         sum += nodeValue(n);
  1968       }
  1969       for (int i = 0; i < blossomNum(); ++i) {
  1970         sum += blossomValue(i) * (blossomSize(i) / 2);
  1971       }
  1972       return sum;
  1973     }
  1974 
  1975     /// \brief Return the dual value (potential) of the given node.
  1976     ///
  1977     /// This function returns the dual value (potential) of the given node.
  1978     ///
  1979     /// \pre Either run() or start() must be called before using this function.
  1980     Value nodeValue(const Node& n) const {
  1981       return (*_node_potential)[n];
  1982     }
  1983 
  1984     /// \brief Return the number of the blossoms in the basis.
  1985     ///
  1986     /// This function returns the number of the blossoms in the basis.
  1987     ///
  1988     /// \pre Either run() or start() must be called before using this function.
  1989     /// \see BlossomIt
  1990     int blossomNum() const {
  1991       return _blossom_potential.size();
  1992     }
  1993 
  1994     /// \brief Return the number of the nodes in the given blossom.
  1995     ///
  1996     /// This function returns the number of the nodes in the given blossom.
  1997     ///
  1998     /// \pre Either run() or start() must be called before using this function.
  1999     /// \see BlossomIt
  2000     int blossomSize(int k) const {
  2001       return _blossom_potential[k].end - _blossom_potential[k].begin;
  2002     }
  2003 
  2004     /// \brief Return the dual value (ptential) of the given blossom.
  2005     ///
  2006     /// This function returns the dual value (ptential) of the given blossom.
  2007     ///
  2008     /// \pre Either run() or start() must be called before using this function.
  2009     Value blossomValue(int k) const {
  2010       return _blossom_potential[k].value;
  2011     }
  2012 
  2013     /// \brief Iterator for obtaining the nodes of a blossom.
  2014     ///
  2015     /// This class provides an iterator for obtaining the nodes of the 
  2016     /// given blossom. It lists a subset of the nodes.
  2017     /// Before using this iterator, you must allocate a 
  2018     /// MaxWeightedMatching class and execute it.
  2019     class BlossomIt {
  2020     public:
  2021 
  2022       /// \brief Constructor.
  2023       ///
  2024       /// Constructor to get the nodes of the given variable.
  2025       ///
  2026       /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or 
  2027       /// \ref MaxWeightedMatching::start() "algorithm.start()" must be 
  2028       /// called before initializing this iterator.
  2029       BlossomIt(const MaxWeightedMatching& algorithm, int variable)
  2030         : _algorithm(&algorithm)
  2031       {
  2032         _index = _algorithm->_blossom_potential[variable].begin;
  2033         _last = _algorithm->_blossom_potential[variable].end;
  2034       }
  2035 
  2036       /// \brief Conversion to \c Node.
  2037       ///
  2038       /// Conversion to \c Node.
  2039       operator Node() const {
  2040         return _algorithm->_blossom_node_list[_index];
  2041       }
  2042 
  2043       /// \brief Increment operator.
  2044       ///
  2045       /// Increment operator.
  2046       BlossomIt& operator++() {
  2047         ++_index;
  2048         return *this;
  2049       }
  2050 
  2051       /// \brief Validity checking
  2052       ///
  2053       /// Checks whether the iterator is invalid.
  2054       bool operator==(Invalid) const { return _index == _last; }
  2055 
  2056       /// \brief Validity checking
  2057       ///
  2058       /// Checks whether the iterator is valid.
  2059       bool operator!=(Invalid) const { return _index != _last; }
  2060 
  2061     private:
  2062       const MaxWeightedMatching* _algorithm;
  2063       int _last;
  2064       int _index;
  2065     };
  2066 
  2067     /// @}
  2068 
  2069   };
  2070 
  2071   /// \ingroup matching
  2072   ///
  2073   /// \brief Weighted perfect matching in general graphs
  2074   ///
  2075   /// This class provides an efficient implementation of Edmond's
  2076   /// maximum weighted perfect matching algorithm. The implementation
  2077   /// is based on extensive use of priority queues and provides
  2078   /// \f$O(nm\log n)\f$ time complexity.
  2079   ///
  2080   /// The maximum weighted perfect matching problem is to find a subset of 
  2081   /// the edges in an undirected graph with maximum overall weight for which 
  2082   /// each node has exactly one incident edge.
  2083   /// It can be formulated with the following linear program.
  2084   /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
  2085   /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
  2086       \quad \forall B\in\mathcal{O}\f] */
  2087   /// \f[x_e \ge 0\quad \forall e\in E\f]
  2088   /// \f[\max \sum_{e\in E}x_ew_e\f]
  2089   /// where \f$\delta(X)\f$ is the set of edges incident to a node in
  2090   /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
  2091   /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
  2092   /// subsets of the nodes.
  2093   ///
  2094   /// The algorithm calculates an optimal matching and a proof of the
  2095   /// optimality. The solution of the dual problem can be used to check
  2096   /// the result of the algorithm. The dual linear problem is the
  2097   /// following.
  2098   /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
  2099       w_{uv} \quad \forall uv\in E\f] */
  2100   /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
  2101   /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
  2102       \frac{\vert B \vert - 1}{2}z_B\f] */
  2103   ///
  2104   /// The algorithm can be executed with the run() function. 
  2105   /// After it the matching (the primal solution) and the dual solution
  2106   /// can be obtained using the query functions and the 
  2107   /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, 
  2108   /// which is able to iterate on the nodes of a blossom. 
  2109   /// If the value type is integer, then the dual solution is multiplied
  2110   /// by \ref MaxWeightedMatching::dualScale "4".
  2111   ///
  2112   /// \tparam GR The undirected graph type the algorithm runs on.
  2113   /// \tparam WM The type edge weight map. The default type is 
  2114   /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
  2115 #ifdef DOXYGEN
  2116   template <typename GR, typename WM>
  2117 #else
  2118   template <typename GR,
  2119             typename WM = typename GR::template EdgeMap<int> >
  2120 #endif
  2121   class MaxWeightedPerfectMatching {
  2122   public:
  2123 
  2124     /// The graph type of the algorithm
  2125     typedef GR Graph;
  2126     /// The type of the edge weight map
  2127     typedef WM WeightMap;
  2128     /// The value type of the edge weights
  2129     typedef typename WeightMap::Value Value;
  2130 
  2131     /// \brief Scaling factor for dual solution
  2132     ///
  2133     /// Scaling factor for dual solution, it is equal to 4 or 1
  2134     /// according to the value type.
  2135     static const int dualScale =
  2136       std::numeric_limits<Value>::is_integer ? 4 : 1;
  2137 
  2138     /// The type of the matching map
  2139     typedef typename Graph::template NodeMap<typename Graph::Arc>
  2140     MatchingMap;
  2141 
  2142   private:
  2143 
  2144     TEMPLATE_GRAPH_TYPEDEFS(Graph);
  2145 
  2146     typedef typename Graph::template NodeMap<Value> NodePotential;
  2147     typedef std::vector<Node> BlossomNodeList;
  2148 
  2149     struct BlossomVariable {
  2150       int begin, end;
  2151       Value value;
  2152 
  2153       BlossomVariable(int _begin, int _end, Value _value)
  2154         : begin(_begin), end(_end), value(_value) {}
  2155 
  2156     };
  2157 
  2158     typedef std::vector<BlossomVariable> BlossomPotential;
  2159 
  2160     const Graph& _graph;
  2161     const WeightMap& _weight;
  2162 
  2163     MatchingMap* _matching;
  2164 
  2165     NodePotential* _node_potential;
  2166 
  2167     BlossomPotential _blossom_potential;
  2168     BlossomNodeList _blossom_node_list;
  2169 
  2170     int _node_num;
  2171     int _blossom_num;
  2172 
  2173     typedef RangeMap<int> IntIntMap;
  2174 
  2175     enum Status {
  2176       EVEN = -1, MATCHED = 0, ODD = 1
  2177     };
  2178 
  2179     typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
  2180     struct BlossomData {
  2181       int tree;
  2182       Status status;
  2183       Arc pred, next;
  2184       Value pot, offset;
  2185     };
  2186 
  2187     IntNodeMap *_blossom_index;
  2188     BlossomSet *_blossom_set;
  2189     RangeMap<BlossomData>* _blossom_data;
  2190 
  2191     IntNodeMap *_node_index;
  2192     IntArcMap *_node_heap_index;
  2193 
  2194     struct NodeData {
  2195 
  2196       NodeData(IntArcMap& node_heap_index)
  2197         : heap(node_heap_index) {}
  2198 
  2199       int blossom;
  2200       Value pot;
  2201       BinHeap<Value, IntArcMap> heap;
  2202       std::map<int, Arc> heap_index;
  2203 
  2204       int tree;
  2205     };
  2206 
  2207     RangeMap<NodeData>* _node_data;
  2208 
  2209     typedef ExtendFindEnum<IntIntMap> TreeSet;
  2210 
  2211     IntIntMap *_tree_set_index;
  2212     TreeSet *_tree_set;
  2213 
  2214     IntIntMap *_delta2_index;
  2215     BinHeap<Value, IntIntMap> *_delta2;
  2216 
  2217     IntEdgeMap *_delta3_index;
  2218     BinHeap<Value, IntEdgeMap> *_delta3;
  2219 
  2220     IntIntMap *_delta4_index;
  2221     BinHeap<Value, IntIntMap> *_delta4;
  2222 
  2223     Value _delta_sum;
  2224 
  2225     void createStructures() {
  2226       _node_num = countNodes(_graph);
  2227       _blossom_num = _node_num * 3 / 2;
  2228 
  2229       if (!_matching) {
  2230         _matching = new MatchingMap(_graph);
  2231       }
  2232 
  2233       if (!_node_potential) {
  2234         _node_potential = new NodePotential(_graph);
  2235       }
  2236 
  2237       if (!_blossom_set) {
  2238         _blossom_index = new IntNodeMap(_graph);
  2239         _blossom_set = new BlossomSet(*_blossom_index);
  2240         _blossom_data = new RangeMap<BlossomData>(_blossom_num);
  2241       } else if (_blossom_data->size() != _blossom_num) {
  2242         delete _blossom_data;
  2243         _blossom_data = new RangeMap<BlossomData>(_blossom_num);
  2244       }
  2245 
  2246       if (!_node_index) {
  2247         _node_index = new IntNodeMap(_graph);
  2248         _node_heap_index = new IntArcMap(_graph);
  2249         _node_data = new RangeMap<NodeData>(_node_num,
  2250                                             NodeData(*_node_heap_index));
  2251       } else if (_node_data->size() != _node_num) {
  2252         delete _node_data;
  2253         _node_data = new RangeMap<NodeData>(_node_num,
  2254                                             NodeData(*_node_heap_index));
  2255       }
  2256 
  2257       if (!_tree_set) {
  2258         _tree_set_index = new IntIntMap(_blossom_num);
  2259         _tree_set = new TreeSet(*_tree_set_index);
  2260       } else {
  2261         _tree_set_index->resize(_blossom_num);
  2262       }
  2263 
  2264       if (!_delta2) {
  2265         _delta2_index = new IntIntMap(_blossom_num);
  2266         _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
  2267       } else {
  2268         _delta2_index->resize(_blossom_num);
  2269       }
  2270 
  2271       if (!_delta3) {
  2272         _delta3_index = new IntEdgeMap(_graph);
  2273         _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
  2274       }
  2275 
  2276       if (!_delta4) {
  2277         _delta4_index = new IntIntMap(_blossom_num);
  2278         _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
  2279       } else {
  2280         _delta4_index->resize(_blossom_num);
  2281       }
  2282     }
  2283 
  2284     void destroyStructures() {
  2285       _node_num = countNodes(_graph);
  2286       _blossom_num = _node_num * 3 / 2;
  2287 
  2288       if (_matching) {
  2289         delete _matching;
  2290       }
  2291       if (_node_potential) {
  2292         delete _node_potential;
  2293       }
  2294       if (_blossom_set) {
  2295         delete _blossom_index;
  2296         delete _blossom_set;
  2297         delete _blossom_data;
  2298       }
  2299 
  2300       if (_node_index) {
  2301         delete _node_index;
  2302         delete _node_heap_index;
  2303         delete _node_data;
  2304       }
  2305 
  2306       if (_tree_set) {
  2307         delete _tree_set_index;
  2308         delete _tree_set;
  2309       }
  2310       if (_delta2) {
  2311         delete _delta2_index;
  2312         delete _delta2;
  2313       }
  2314       if (_delta3) {
  2315         delete _delta3_index;
  2316         delete _delta3;
  2317       }
  2318       if (_delta4) {
  2319         delete _delta4_index;
  2320         delete _delta4;
  2321       }
  2322     }
  2323 
  2324     void matchedToEven(int blossom, int tree) {
  2325       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
  2326         _delta2->erase(blossom);
  2327       }
  2328 
  2329       if (!_blossom_set->trivial(blossom)) {
  2330         (*_blossom_data)[blossom].pot -=
  2331           2 * (_delta_sum - (*_blossom_data)[blossom].offset);
  2332       }
  2333 
  2334       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  2335            n != INVALID; ++n) {
  2336 
  2337         _blossom_set->increase(n, std::numeric_limits<Value>::max());
  2338         int ni = (*_node_index)[n];
  2339 
  2340         (*_node_data)[ni].heap.clear();
  2341         (*_node_data)[ni].heap_index.clear();
  2342 
  2343         (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
  2344 
  2345         for (InArcIt e(_graph, n); e != INVALID; ++e) {
  2346           Node v = _graph.source(e);
  2347           int vb = _blossom_set->find(v);
  2348           int vi = (*_node_index)[v];
  2349 
  2350           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  2351             dualScale * _weight[e];
  2352 
  2353           if ((*_blossom_data)[vb].status == EVEN) {
  2354             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
  2355               _delta3->push(e, rw / 2);
  2356             }
  2357           } else {
  2358             typename std::map<int, Arc>::iterator it =
  2359               (*_node_data)[vi].heap_index.find(tree);
  2360 
  2361             if (it != (*_node_data)[vi].heap_index.end()) {
  2362               if ((*_node_data)[vi].heap[it->second] > rw) {
  2363                 (*_node_data)[vi].heap.replace(it->second, e);
  2364                 (*_node_data)[vi].heap.decrease(e, rw);
  2365                 it->second = e;
  2366               }
  2367             } else {
  2368               (*_node_data)[vi].heap.push(e, rw);
  2369               (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
  2370             }
  2371 
  2372             if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
  2373               _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
  2374 
  2375               if ((*_blossom_data)[vb].status == MATCHED) {
  2376                 if (_delta2->state(vb) != _delta2->IN_HEAP) {
  2377                   _delta2->push(vb, _blossom_set->classPrio(vb) -
  2378                                (*_blossom_data)[vb].offset);
  2379                 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
  2380                            (*_blossom_data)[vb].offset){
  2381                   _delta2->decrease(vb, _blossom_set->classPrio(vb) -
  2382                                    (*_blossom_data)[vb].offset);
  2383                 }
  2384               }
  2385             }
  2386           }
  2387         }
  2388       }
  2389       (*_blossom_data)[blossom].offset = 0;
  2390     }
  2391 
  2392     void matchedToOdd(int blossom) {
  2393       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
  2394         _delta2->erase(blossom);
  2395       }
  2396       (*_blossom_data)[blossom].offset += _delta_sum;
  2397       if (!_blossom_set->trivial(blossom)) {
  2398         _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
  2399                      (*_blossom_data)[blossom].offset);
  2400       }
  2401     }
  2402 
  2403     void evenToMatched(int blossom, int tree) {
  2404       if (!_blossom_set->trivial(blossom)) {
  2405         (*_blossom_data)[blossom].pot += 2 * _delta_sum;
  2406       }
  2407 
  2408       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  2409            n != INVALID; ++n) {
  2410         int ni = (*_node_index)[n];
  2411         (*_node_data)[ni].pot -= _delta_sum;
  2412 
  2413         for (InArcIt e(_graph, n); e != INVALID; ++e) {
  2414           Node v = _graph.source(e);
  2415           int vb = _blossom_set->find(v);
  2416           int vi = (*_node_index)[v];
  2417 
  2418           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  2419             dualScale * _weight[e];
  2420 
  2421           if (vb == blossom) {
  2422             if (_delta3->state(e) == _delta3->IN_HEAP) {
  2423               _delta3->erase(e);
  2424             }
  2425           } else if ((*_blossom_data)[vb].status == EVEN) {
  2426 
  2427             if (_delta3->state(e) == _delta3->IN_HEAP) {
  2428               _delta3->erase(e);
  2429             }
  2430 
  2431             int vt = _tree_set->find(vb);
  2432 
  2433             if (vt != tree) {
  2434 
  2435               Arc r = _graph.oppositeArc(e);
  2436 
  2437               typename std::map<int, Arc>::iterator it =
  2438                 (*_node_data)[ni].heap_index.find(vt);
  2439 
  2440               if (it != (*_node_data)[ni].heap_index.end()) {
  2441                 if ((*_node_data)[ni].heap[it->second] > rw) {
  2442                   (*_node_data)[ni].heap.replace(it->second, r);
  2443                   (*_node_data)[ni].heap.decrease(r, rw);
  2444                   it->second = r;
  2445                 }
  2446               } else {
  2447                 (*_node_data)[ni].heap.push(r, rw);
  2448                 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
  2449               }
  2450 
  2451               if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
  2452                 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
  2453 
  2454                 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
  2455                   _delta2->push(blossom, _blossom_set->classPrio(blossom) -
  2456                                (*_blossom_data)[blossom].offset);
  2457                 } else if ((*_delta2)[blossom] >
  2458                            _blossom_set->classPrio(blossom) -
  2459                            (*_blossom_data)[blossom].offset){
  2460                   _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
  2461                                    (*_blossom_data)[blossom].offset);
  2462                 }
  2463               }
  2464             }
  2465           } else {
  2466 
  2467             typename std::map<int, Arc>::iterator it =
  2468               (*_node_data)[vi].heap_index.find(tree);
  2469 
  2470             if (it != (*_node_data)[vi].heap_index.end()) {
  2471               (*_node_data)[vi].heap.erase(it->second);
  2472               (*_node_data)[vi].heap_index.erase(it);
  2473               if ((*_node_data)[vi].heap.empty()) {
  2474                 _blossom_set->increase(v, std::numeric_limits<Value>::max());
  2475               } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
  2476                 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
  2477               }
  2478 
  2479               if ((*_blossom_data)[vb].status == MATCHED) {
  2480                 if (_blossom_set->classPrio(vb) ==
  2481                     std::numeric_limits<Value>::max()) {
  2482                   _delta2->erase(vb);
  2483                 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
  2484                            (*_blossom_data)[vb].offset) {
  2485                   _delta2->increase(vb, _blossom_set->classPrio(vb) -
  2486                                    (*_blossom_data)[vb].offset);
  2487                 }
  2488               }
  2489             }
  2490           }
  2491         }
  2492       }
  2493     }
  2494 
  2495     void oddToMatched(int blossom) {
  2496       (*_blossom_data)[blossom].offset -= _delta_sum;
  2497 
  2498       if (_blossom_set->classPrio(blossom) !=
  2499           std::numeric_limits<Value>::max()) {
  2500         _delta2->push(blossom, _blossom_set->classPrio(blossom) -
  2501                        (*_blossom_data)[blossom].offset);
  2502       }
  2503 
  2504       if (!_blossom_set->trivial(blossom)) {
  2505         _delta4->erase(blossom);
  2506       }
  2507     }
  2508 
  2509     void oddToEven(int blossom, int tree) {
  2510       if (!_blossom_set->trivial(blossom)) {
  2511         _delta4->erase(blossom);
  2512         (*_blossom_data)[blossom].pot -=
  2513           2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
  2514       }
  2515 
  2516       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  2517            n != INVALID; ++n) {
  2518         int ni = (*_node_index)[n];
  2519 
  2520         _blossom_set->increase(n, std::numeric_limits<Value>::max());
  2521 
  2522         (*_node_data)[ni].heap.clear();
  2523         (*_node_data)[ni].heap_index.clear();
  2524         (*_node_data)[ni].pot +=
  2525           2 * _delta_sum - (*_blossom_data)[blossom].offset;
  2526 
  2527         for (InArcIt e(_graph, n); e != INVALID; ++e) {
  2528           Node v = _graph.source(e);
  2529           int vb = _blossom_set->find(v);
  2530           int vi = (*_node_index)[v];
  2531 
  2532           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  2533             dualScale * _weight[e];
  2534 
  2535           if ((*_blossom_data)[vb].status == EVEN) {
  2536             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
  2537               _delta3->push(e, rw / 2);
  2538             }
  2539           } else {
  2540 
  2541             typename std::map<int, Arc>::iterator it =
  2542               (*_node_data)[vi].heap_index.find(tree);
  2543 
  2544             if (it != (*_node_data)[vi].heap_index.end()) {
  2545               if ((*_node_data)[vi].heap[it->second] > rw) {
  2546                 (*_node_data)[vi].heap.replace(it->second, e);
  2547                 (*_node_data)[vi].heap.decrease(e, rw);
  2548                 it->second = e;
  2549               }
  2550             } else {
  2551               (*_node_data)[vi].heap.push(e, rw);
  2552               (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
  2553             }
  2554 
  2555             if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
  2556               _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
  2557 
  2558               if ((*_blossom_data)[vb].status == MATCHED) {
  2559                 if (_delta2->state(vb) != _delta2->IN_HEAP) {
  2560                   _delta2->push(vb, _blossom_set->classPrio(vb) -
  2561                                (*_blossom_data)[vb].offset);
  2562                 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
  2563                            (*_blossom_data)[vb].offset) {
  2564                   _delta2->decrease(vb, _blossom_set->classPrio(vb) -
  2565                                    (*_blossom_data)[vb].offset);
  2566                 }
  2567               }
  2568             }
  2569           }
  2570         }
  2571       }
  2572       (*_blossom_data)[blossom].offset = 0;
  2573     }
  2574 
  2575     void alternatePath(int even, int tree) {
  2576       int odd;
  2577 
  2578       evenToMatched(even, tree);
  2579       (*_blossom_data)[even].status = MATCHED;
  2580 
  2581       while ((*_blossom_data)[even].pred != INVALID) {
  2582         odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
  2583         (*_blossom_data)[odd].status = MATCHED;
  2584         oddToMatched(odd);
  2585         (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
  2586 
  2587         even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
  2588         (*_blossom_data)[even].status = MATCHED;
  2589         evenToMatched(even, tree);
  2590         (*_blossom_data)[even].next =
  2591           _graph.oppositeArc((*_blossom_data)[odd].pred);
  2592       }
  2593 
  2594     }
  2595 
  2596     void destroyTree(int tree) {
  2597       for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
  2598         if ((*_blossom_data)[b].status == EVEN) {
  2599           (*_blossom_data)[b].status = MATCHED;
  2600           evenToMatched(b, tree);
  2601         } else if ((*_blossom_data)[b].status == ODD) {
  2602           (*_blossom_data)[b].status = MATCHED;
  2603           oddToMatched(b);
  2604         }
  2605       }
  2606       _tree_set->eraseClass(tree);
  2607     }
  2608 
  2609     void augmentOnEdge(const Edge& edge) {
  2610 
  2611       int left = _blossom_set->find(_graph.u(edge));
  2612       int right = _blossom_set->find(_graph.v(edge));
  2613 
  2614       int left_tree = _tree_set->find(left);
  2615       alternatePath(left, left_tree);
  2616       destroyTree(left_tree);
  2617 
  2618       int right_tree = _tree_set->find(right);
  2619       alternatePath(right, right_tree);
  2620       destroyTree(right_tree);
  2621 
  2622       (*_blossom_data)[left].next = _graph.direct(edge, true);
  2623       (*_blossom_data)[right].next = _graph.direct(edge, false);
  2624     }
  2625 
  2626     void extendOnArc(const Arc& arc) {
  2627       int base = _blossom_set->find(_graph.target(arc));
  2628       int tree = _tree_set->find(base);
  2629 
  2630       int odd = _blossom_set->find(_graph.source(arc));
  2631       _tree_set->insert(odd, tree);
  2632       (*_blossom_data)[odd].status = ODD;
  2633       matchedToOdd(odd);
  2634       (*_blossom_data)[odd].pred = arc;
  2635 
  2636       int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
  2637       (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
  2638       _tree_set->insert(even, tree);
  2639       (*_blossom_data)[even].status = EVEN;
  2640       matchedToEven(even, tree);
  2641     }
  2642 
  2643     void shrinkOnEdge(const Edge& edge, int tree) {
  2644       int nca = -1;
  2645       std::vector<int> left_path, right_path;
  2646 
  2647       {
  2648         std::set<int> left_set, right_set;
  2649         int left = _blossom_set->find(_graph.u(edge));
  2650         left_path.push_back(left);
  2651         left_set.insert(left);
  2652 
  2653         int right = _blossom_set->find(_graph.v(edge));
  2654         right_path.push_back(right);
  2655         right_set.insert(right);
  2656 
  2657         while (true) {
  2658 
  2659           if ((*_blossom_data)[left].pred == INVALID) break;
  2660 
  2661           left =
  2662             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
  2663           left_path.push_back(left);
  2664           left =
  2665             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
  2666           left_path.push_back(left);
  2667 
  2668           left_set.insert(left);
  2669 
  2670           if (right_set.find(left) != right_set.end()) {
  2671             nca = left;
  2672             break;
  2673           }
  2674 
  2675           if ((*_blossom_data)[right].pred == INVALID) break;
  2676 
  2677           right =
  2678             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
  2679           right_path.push_back(right);
  2680           right =
  2681             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
  2682           right_path.push_back(right);
  2683 
  2684           right_set.insert(right);
  2685 
  2686           if (left_set.find(right) != left_set.end()) {
  2687             nca = right;
  2688             break;
  2689           }
  2690 
  2691         }
  2692 
  2693         if (nca == -1) {
  2694           if ((*_blossom_data)[left].pred == INVALID) {
  2695             nca = right;
  2696             while (left_set.find(nca) == left_set.end()) {
  2697               nca =
  2698                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  2699               right_path.push_back(nca);
  2700               nca =
  2701                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  2702               right_path.push_back(nca);
  2703             }
  2704           } else {
  2705             nca = left;
  2706             while (right_set.find(nca) == right_set.end()) {
  2707               nca =
  2708                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  2709               left_path.push_back(nca);
  2710               nca =
  2711                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  2712               left_path.push_back(nca);
  2713             }
  2714           }
  2715         }
  2716       }
  2717 
  2718       std::vector<int> subblossoms;
  2719       Arc prev;
  2720 
  2721       prev = _graph.direct(edge, true);
  2722       for (int i = 0; left_path[i] != nca; i += 2) {
  2723         subblossoms.push_back(left_path[i]);
  2724         (*_blossom_data)[left_path[i]].next = prev;
  2725         _tree_set->erase(left_path[i]);
  2726 
  2727         subblossoms.push_back(left_path[i + 1]);
  2728         (*_blossom_data)[left_path[i + 1]].status = EVEN;
  2729         oddToEven(left_path[i + 1], tree);
  2730         _tree_set->erase(left_path[i + 1]);
  2731         prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
  2732       }
  2733 
  2734       int k = 0;
  2735       while (right_path[k] != nca) ++k;
  2736 
  2737       subblossoms.push_back(nca);
  2738       (*_blossom_data)[nca].next = prev;
  2739 
  2740       for (int i = k - 2; i >= 0; i -= 2) {
  2741         subblossoms.push_back(right_path[i + 1]);
  2742         (*_blossom_data)[right_path[i + 1]].status = EVEN;
  2743         oddToEven(right_path[i + 1], tree);
  2744         _tree_set->erase(right_path[i + 1]);
  2745 
  2746         (*_blossom_data)[right_path[i + 1]].next =
  2747           (*_blossom_data)[right_path[i + 1]].pred;
  2748 
  2749         subblossoms.push_back(right_path[i]);
  2750         _tree_set->erase(right_path[i]);
  2751       }
  2752 
  2753       int surface =
  2754         _blossom_set->join(subblossoms.begin(), subblossoms.end());
  2755 
  2756       for (int i = 0; i < int(subblossoms.size()); ++i) {
  2757         if (!_blossom_set->trivial(subblossoms[i])) {
  2758           (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
  2759         }
  2760         (*_blossom_data)[subblossoms[i]].status = MATCHED;
  2761       }
  2762 
  2763       (*_blossom_data)[surface].pot = -2 * _delta_sum;
  2764       (*_blossom_data)[surface].offset = 0;
  2765       (*_blossom_data)[surface].status = EVEN;
  2766       (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
  2767       (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
  2768 
  2769       _tree_set->insert(surface, tree);
  2770       _tree_set->erase(nca);
  2771     }
  2772 
  2773     void splitBlossom(int blossom) {
  2774       Arc next = (*_blossom_data)[blossom].next;
  2775       Arc pred = (*_blossom_data)[blossom].pred;
  2776 
  2777       int tree = _tree_set->find(blossom);
  2778 
  2779       (*_blossom_data)[blossom].status = MATCHED;
  2780       oddToMatched(blossom);
  2781       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
  2782         _delta2->erase(blossom);
  2783       }
  2784 
  2785       std::vector<int> subblossoms;
  2786       _blossom_set->split(blossom, std::back_inserter(subblossoms));
  2787 
  2788       Value offset = (*_blossom_data)[blossom].offset;
  2789       int b = _blossom_set->find(_graph.source(pred));
  2790       int d = _blossom_set->find(_graph.source(next));
  2791 
  2792       int ib = -1, id = -1;
  2793       for (int i = 0; i < int(subblossoms.size()); ++i) {
  2794         if (subblossoms[i] == b) ib = i;
  2795         if (subblossoms[i] == d) id = i;
  2796 
  2797         (*_blossom_data)[subblossoms[i]].offset = offset;
  2798         if (!_blossom_set->trivial(subblossoms[i])) {
  2799           (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
  2800         }
  2801         if (_blossom_set->classPrio(subblossoms[i]) !=
  2802             std::numeric_limits<Value>::max()) {
  2803           _delta2->push(subblossoms[i],
  2804                         _blossom_set->classPrio(subblossoms[i]) -
  2805                         (*_blossom_data)[subblossoms[i]].offset);
  2806         }
  2807       }
  2808 
  2809       if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
  2810         for (int i = (id + 1) % subblossoms.size();
  2811              i != ib; i = (i + 2) % subblossoms.size()) {
  2812           int sb = subblossoms[i];
  2813           int tb = subblossoms[(i + 1) % subblossoms.size()];
  2814           (*_blossom_data)[sb].next =
  2815             _graph.oppositeArc((*_blossom_data)[tb].next);
  2816         }
  2817 
  2818         for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
  2819           int sb = subblossoms[i];
  2820           int tb = subblossoms[(i + 1) % subblossoms.size()];
  2821           int ub = subblossoms[(i + 2) % subblossoms.size()];
  2822 
  2823           (*_blossom_data)[sb].status = ODD;
  2824           matchedToOdd(sb);
  2825           _tree_set->insert(sb, tree);
  2826           (*_blossom_data)[sb].pred = pred;
  2827           (*_blossom_data)[sb].next =
  2828                            _graph.oppositeArc((*_blossom_data)[tb].next);
  2829 
  2830           pred = (*_blossom_data)[ub].next;
  2831 
  2832           (*_blossom_data)[tb].status = EVEN;
  2833           matchedToEven(tb, tree);
  2834           _tree_set->insert(tb, tree);
  2835           (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
  2836         }
  2837 
  2838         (*_blossom_data)[subblossoms[id]].status = ODD;
  2839         matchedToOdd(subblossoms[id]);
  2840         _tree_set->insert(subblossoms[id], tree);
  2841         (*_blossom_data)[subblossoms[id]].next = next;
  2842         (*_blossom_data)[subblossoms[id]].pred = pred;
  2843 
  2844       } else {
  2845 
  2846         for (int i = (ib + 1) % subblossoms.size();
  2847              i != id; i = (i + 2) % subblossoms.size()) {
  2848           int sb = subblossoms[i];
  2849           int tb = subblossoms[(i + 1) % subblossoms.size()];
  2850           (*_blossom_data)[sb].next =
  2851             _graph.oppositeArc((*_blossom_data)[tb].next);
  2852         }
  2853 
  2854         for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
  2855           int sb = subblossoms[i];
  2856           int tb = subblossoms[(i + 1) % subblossoms.size()];
  2857           int ub = subblossoms[(i + 2) % subblossoms.size()];
  2858 
  2859           (*_blossom_data)[sb].status = ODD;
  2860           matchedToOdd(sb);
  2861           _tree_set->insert(sb, tree);
  2862           (*_blossom_data)[sb].next = next;
  2863           (*_blossom_data)[sb].pred =
  2864             _graph.oppositeArc((*_blossom_data)[tb].next);
  2865 
  2866           (*_blossom_data)[tb].status = EVEN;
  2867           matchedToEven(tb, tree);
  2868           _tree_set->insert(tb, tree);
  2869           (*_blossom_data)[tb].pred =
  2870             (*_blossom_data)[tb].next =
  2871             _graph.oppositeArc((*_blossom_data)[ub].next);
  2872           next = (*_blossom_data)[ub].next;
  2873         }
  2874 
  2875         (*_blossom_data)[subblossoms[ib]].status = ODD;
  2876         matchedToOdd(subblossoms[ib]);
  2877         _tree_set->insert(subblossoms[ib], tree);
  2878         (*_blossom_data)[subblossoms[ib]].next = next;
  2879         (*_blossom_data)[subblossoms[ib]].pred = pred;
  2880       }
  2881       _tree_set->erase(blossom);
  2882     }
  2883 
  2884     void extractBlossom(int blossom, const Node& base, const Arc& matching) {
  2885       if (_blossom_set->trivial(blossom)) {
  2886         int bi = (*_node_index)[base];
  2887         Value pot = (*_node_data)[bi].pot;
  2888 
  2889         (*_matching)[base] = matching;
  2890         _blossom_node_list.push_back(base);
  2891         (*_node_potential)[base] = pot;
  2892       } else {
  2893 
  2894         Value pot = (*_blossom_data)[blossom].pot;
  2895         int bn = _blossom_node_list.size();
  2896 
  2897         std::vector<int> subblossoms;
  2898         _blossom_set->split(blossom, std::back_inserter(subblossoms));
  2899         int b = _blossom_set->find(base);
  2900         int ib = -1;
  2901         for (int i = 0; i < int(subblossoms.size()); ++i) {
  2902           if (subblossoms[i] == b) { ib = i; break; }
  2903         }
  2904 
  2905         for (int i = 1; i < int(subblossoms.size()); i += 2) {
  2906           int sb = subblossoms[(ib + i) % subblossoms.size()];
  2907           int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
  2908 
  2909           Arc m = (*_blossom_data)[tb].next;
  2910           extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
  2911           extractBlossom(tb, _graph.source(m), m);
  2912         }
  2913         extractBlossom(subblossoms[ib], base, matching);
  2914 
  2915         int en = _blossom_node_list.size();
  2916 
  2917         _blossom_potential.push_back(BlossomVariable(bn, en, pot));
  2918       }
  2919     }
  2920 
  2921     void extractMatching() {
  2922       std::vector<int> blossoms;
  2923       for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
  2924         blossoms.push_back(c);
  2925       }
  2926 
  2927       for (int i = 0; i < int(blossoms.size()); ++i) {
  2928 
  2929         Value offset = (*_blossom_data)[blossoms[i]].offset;
  2930         (*_blossom_data)[blossoms[i]].pot += 2 * offset;
  2931         for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
  2932              n != INVALID; ++n) {
  2933           (*_node_data)[(*_node_index)[n]].pot -= offset;
  2934         }
  2935 
  2936         Arc matching = (*_blossom_data)[blossoms[i]].next;
  2937         Node base = _graph.source(matching);
  2938         extractBlossom(blossoms[i], base, matching);
  2939       }
  2940     }
  2941 
  2942   public:
  2943 
  2944     /// \brief Constructor
  2945     ///
  2946     /// Constructor.
  2947     MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
  2948       : _graph(graph), _weight(weight), _matching(0),
  2949         _node_potential(0), _blossom_potential(), _blossom_node_list(),
  2950         _node_num(0), _blossom_num(0),
  2951 
  2952         _blossom_index(0), _blossom_set(0), _blossom_data(0),
  2953         _node_index(0), _node_heap_index(0), _node_data(0),
  2954         _tree_set_index(0), _tree_set(0),
  2955 
  2956         _delta2_index(0), _delta2(0),
  2957         _delta3_index(0), _delta3(0),
  2958         _delta4_index(0), _delta4(0),
  2959 
  2960         _delta_sum() {}
  2961 
  2962     ~MaxWeightedPerfectMatching() {
  2963       destroyStructures();
  2964     }
  2965 
  2966     /// \name Execution Control
  2967     /// The simplest way to execute the algorithm is to use the
  2968     /// \ref run() member function.
  2969 
  2970     ///@{
  2971 
  2972     /// \brief Initialize the algorithm
  2973     ///
  2974     /// This function initializes the algorithm.
  2975     void init() {
  2976       createStructures();
  2977 
  2978       _blossom_node_list.clear();
  2979       _blossom_potential.clear();
  2980 
  2981       for (ArcIt e(_graph); e != INVALID; ++e) {
  2982         (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
  2983       }
  2984       for (EdgeIt e(_graph); e != INVALID; ++e) {
  2985         (*_delta3_index)[e] = _delta3->PRE_HEAP;
  2986       }
  2987       for (int i = 0; i < _blossom_num; ++i) {
  2988         (*_delta2_index)[i] = _delta2->PRE_HEAP;
  2989         (*_delta4_index)[i] = _delta4->PRE_HEAP;
  2990       }
  2991 
  2992       _delta2->clear();
  2993       _delta3->clear();
  2994       _delta4->clear();
  2995       _blossom_set->clear();
  2996       _tree_set->clear();
  2997 
  2998       int index = 0;
  2999       for (NodeIt n(_graph); n != INVALID; ++n) {
  3000         Value max = - std::numeric_limits<Value>::max();
  3001         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  3002           if (_graph.target(e) == n) continue;
  3003           if ((dualScale * _weight[e]) / 2 > max) {
  3004             max = (dualScale * _weight[e]) / 2;
  3005           }
  3006         }
  3007         (*_node_index)[n] = index;
  3008         (*_node_data)[index].heap_index.clear();
  3009         (*_node_data)[index].heap.clear();
  3010         (*_node_data)[index].pot = max;
  3011         int blossom =
  3012           _blossom_set->insert(n, std::numeric_limits<Value>::max());
  3013 
  3014         _tree_set->insert(blossom);
  3015 
  3016         (*_blossom_data)[blossom].status = EVEN;
  3017         (*_blossom_data)[blossom].pred = INVALID;
  3018         (*_blossom_data)[blossom].next = INVALID;
  3019         (*_blossom_data)[blossom].pot = 0;
  3020         (*_blossom_data)[blossom].offset = 0;
  3021         ++index;
  3022       }
  3023       for (EdgeIt e(_graph); e != INVALID; ++e) {
  3024         int si = (*_node_index)[_graph.u(e)];
  3025         int ti = (*_node_index)[_graph.v(e)];
  3026         if (_graph.u(e) != _graph.v(e)) {
  3027           _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
  3028                             dualScale * _weight[e]) / 2);
  3029         }
  3030       }
  3031     }
  3032 
  3033     /// \brief Start the algorithm
  3034     ///
  3035     /// This function starts the algorithm.
  3036     ///
  3037     /// \pre \ref init() must be called before using this function.
  3038     bool start() {
  3039       enum OpType {
  3040         D2, D3, D4
  3041       };
  3042 
  3043       int unmatched = _node_num;
  3044       while (unmatched > 0) {
  3045         Value d2 = !_delta2->empty() ?
  3046           _delta2->prio() : std::numeric_limits<Value>::max();
  3047 
  3048         Value d3 = !_delta3->empty() ?
  3049           _delta3->prio() : std::numeric_limits<Value>::max();
  3050 
  3051         Value d4 = !_delta4->empty() ?
  3052           _delta4->prio() : std::numeric_limits<Value>::max();
  3053 
  3054         _delta_sum = d2; OpType ot = D2;
  3055         if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
  3056         if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
  3057 
  3058         if (_delta_sum == std::numeric_limits<Value>::max()) {
  3059           return false;
  3060         }
  3061 
  3062         switch (ot) {
  3063         case D2:
  3064           {
  3065             int blossom = _delta2->top();
  3066             Node n = _blossom_set->classTop(blossom);
  3067             Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
  3068             extendOnArc(e);
  3069           }
  3070           break;
  3071         case D3:
  3072           {
  3073             Edge e = _delta3->top();
  3074 
  3075             int left_blossom = _blossom_set->find(_graph.u(e));
  3076             int right_blossom = _blossom_set->find(_graph.v(e));
  3077 
  3078             if (left_blossom == right_blossom) {
  3079               _delta3->pop();
  3080             } else {
  3081               int left_tree = _tree_set->find(left_blossom);
  3082               int right_tree = _tree_set->find(right_blossom);
  3083 
  3084               if (left_tree == right_tree) {
  3085                 shrinkOnEdge(e, left_tree);
  3086               } else {
  3087                 augmentOnEdge(e);
  3088                 unmatched -= 2;
  3089               }
  3090             }
  3091           } break;
  3092         case D4:
  3093           splitBlossom(_delta4->top());
  3094           break;
  3095         }
  3096       }
  3097       extractMatching();
  3098       return true;
  3099     }
  3100 
  3101     /// \brief Run the algorithm.
  3102     ///
  3103     /// This method runs the \c %MaxWeightedPerfectMatching algorithm.
  3104     ///
  3105     /// \note mwpm.run() is just a shortcut of the following code.
  3106     /// \code
  3107     ///   mwpm.init();
  3108     ///   mwpm.start();
  3109     /// \endcode
  3110     bool run() {
  3111       init();
  3112       return start();
  3113     }
  3114 
  3115     /// @}
  3116 
  3117     /// \name Primal Solution
  3118     /// Functions to get the primal solution, i.e. the maximum weighted 
  3119     /// perfect matching.\n
  3120     /// Either \ref run() or \ref start() function should be called before
  3121     /// using them.
  3122 
  3123     /// @{
  3124 
  3125     /// \brief Return the weight of the matching.
  3126     ///
  3127     /// This function returns the weight of the found matching.
  3128     ///
  3129     /// \pre Either run() or start() must be called before using this function.
  3130     Value matchingWeight() const {
  3131       Value sum = 0;
  3132       for (NodeIt n(_graph); n != INVALID; ++n) {
  3133         if ((*_matching)[n] != INVALID) {
  3134           sum += _weight[(*_matching)[n]];
  3135         }
  3136       }
  3137       return sum /= 2;
  3138     }
  3139 
  3140     /// \brief Return \c true if the given edge is in the matching.
  3141     ///
  3142     /// This function returns \c true if the given edge is in the found 
  3143     /// matching.
  3144     ///
  3145     /// \pre Either run() or start() must be called before using this function.
  3146     bool matching(const Edge& edge) const {
  3147       return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
  3148     }
  3149 
  3150     /// \brief Return the matching arc (or edge) incident to the given node.
  3151     ///
  3152     /// This function returns the matching arc (or edge) incident to the
  3153     /// given node in the found matching or \c INVALID if the node is 
  3154     /// not covered by the matching.
  3155     ///
  3156     /// \pre Either run() or start() must be called before using this function.
  3157     Arc matching(const Node& node) const {
  3158       return (*_matching)[node];
  3159     }
  3160 
  3161     /// \brief Return a const reference to the matching map.
  3162     ///
  3163     /// This function returns a const reference to a node map that stores
  3164     /// the matching arc (or edge) incident to each node.
  3165     const MatchingMap& matchingMap() const {
  3166       return *_matching;
  3167     }
  3168 
  3169     /// \brief Return the mate of the given node.
  3170     ///
  3171     /// This function returns the mate of the given node in the found 
  3172     /// matching or \c INVALID if the node is not covered by the matching.
  3173     ///
  3174     /// \pre Either run() or start() must be called before using this function.
  3175     Node mate(const Node& node) const {
  3176       return _graph.target((*_matching)[node]);
  3177     }
  3178 
  3179     /// @}
  3180 
  3181     /// \name Dual Solution
  3182     /// Functions to get the dual solution.\n
  3183     /// Either \ref run() or \ref start() function should be called before
  3184     /// using them.
  3185 
  3186     /// @{
  3187 
  3188     /// \brief Return the value of the dual solution.
  3189     ///
  3190     /// This function returns the value of the dual solution. 
  3191     /// It should be equal to the primal value scaled by \ref dualScale 
  3192     /// "dual scale".
  3193     ///
  3194     /// \pre Either run() or start() must be called before using this function.
  3195     Value dualValue() const {
  3196       Value sum = 0;
  3197       for (NodeIt n(_graph); n != INVALID; ++n) {
  3198         sum += nodeValue(n);
  3199       }
  3200       for (int i = 0; i < blossomNum(); ++i) {
  3201         sum += blossomValue(i) * (blossomSize(i) / 2);
  3202       }
  3203       return sum;
  3204     }
  3205 
  3206     /// \brief Return the dual value (potential) of the given node.
  3207     ///
  3208     /// This function returns the dual value (potential) of the given node.
  3209     ///
  3210     /// \pre Either run() or start() must be called before using this function.
  3211     Value nodeValue(const Node& n) const {
  3212       return (*_node_potential)[n];
  3213     }
  3214 
  3215     /// \brief Return the number of the blossoms in the basis.
  3216     ///
  3217     /// This function returns the number of the blossoms in the basis.
  3218     ///
  3219     /// \pre Either run() or start() must be called before using this function.
  3220     /// \see BlossomIt
  3221     int blossomNum() const {
  3222       return _blossom_potential.size();
  3223     }
  3224 
  3225     /// \brief Return the number of the nodes in the given blossom.
  3226     ///
  3227     /// This function returns the number of the nodes in the given blossom.
  3228     ///
  3229     /// \pre Either run() or start() must be called before using this function.
  3230     /// \see BlossomIt
  3231     int blossomSize(int k) const {
  3232       return _blossom_potential[k].end - _blossom_potential[k].begin;
  3233     }
  3234 
  3235     /// \brief Return the dual value (ptential) of the given blossom.
  3236     ///
  3237     /// This function returns the dual value (ptential) of the given blossom.
  3238     ///
  3239     /// \pre Either run() or start() must be called before using this function.
  3240     Value blossomValue(int k) const {
  3241       return _blossom_potential[k].value;
  3242     }
  3243 
  3244     /// \brief Iterator for obtaining the nodes of a blossom.
  3245     ///
  3246     /// This class provides an iterator for obtaining the nodes of the 
  3247     /// given blossom. It lists a subset of the nodes.
  3248     /// Before using this iterator, you must allocate a 
  3249     /// MaxWeightedPerfectMatching class and execute it.
  3250     class BlossomIt {
  3251     public:
  3252 
  3253       /// \brief Constructor.
  3254       ///
  3255       /// Constructor to get the nodes of the given variable.
  3256       ///
  3257       /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" 
  3258       /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" 
  3259       /// must be called before initializing this iterator.
  3260       BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
  3261         : _algorithm(&algorithm)
  3262       {
  3263         _index = _algorithm->_blossom_potential[variable].begin;
  3264         _last = _algorithm->_blossom_potential[variable].end;
  3265       }
  3266 
  3267       /// \brief Conversion to \c Node.
  3268       ///
  3269       /// Conversion to \c Node.
  3270       operator Node() const {
  3271         return _algorithm->_blossom_node_list[_index];
  3272       }
  3273 
  3274       /// \brief Increment operator.
  3275       ///
  3276       /// Increment operator.
  3277       BlossomIt& operator++() {
  3278         ++_index;
  3279         return *this;
  3280       }
  3281 
  3282       /// \brief Validity checking
  3283       ///
  3284       /// This function checks whether the iterator is invalid.
  3285       bool operator==(Invalid) const { return _index == _last; }
  3286 
  3287       /// \brief Validity checking
  3288       ///
  3289       /// This function checks whether the iterator is valid.
  3290       bool operator!=(Invalid) const { return _index != _last; }
  3291 
  3292     private:
  3293       const MaxWeightedPerfectMatching* _algorithm;
  3294       int _last;
  3295       int _index;
  3296     };
  3297 
  3298     /// @}
  3299 
  3300   };
  3301 
  3302 } //END OF NAMESPACE LEMON
  3303 
  3304 #endif //LEMON_MAX_MATCHING_H