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