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