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