lemon/max_matching.h
author Alpar Juttner <alpar@cs.elte.hu>
Mon, 01 Dec 2008 14:11:31 +0000
changeset 400 fa341dd6ab23
parent 327 91d63b8b1a4c
child 425 cace3206223b
permissions -rw-r--r--
Transform circulation demo to test
     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-2008
     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 _Graph The graph type the algorithm runs on.
    59   template <typename _Graph>
    60   class MaxMatching {
    61   public:
    62 
    63     typedef _Graph 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     /// startParse() 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 %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 _Graph,
   651             typename _WeightMap = typename _Graph::template EdgeMap<int> >
   652   class MaxWeightedMatching {
   653   public:
   654 
   655     typedef _Graph Graph;
   656     typedef _WeightMap WeightMap;
   657     typedef typename WeightMap::Value Value;
   658 
   659     /// \brief Scaling factor for dual solution
   660     ///
   661     /// Scaling factor for dual solution, it is equal to 4 or 1
   662     /// according to the value type.
   663     static const int dualScale =
   664       std::numeric_limits<Value>::is_integer ? 4 : 1;
   665 
   666     typedef typename Graph::template NodeMap<typename Graph::Arc>
   667     MatchingMap;
   668 
   669   private:
   670 
   671     TEMPLATE_GRAPH_TYPEDEFS(Graph);
   672 
   673     typedef typename Graph::template NodeMap<Value> NodePotential;
   674     typedef std::vector<Node> BlossomNodeList;
   675 
   676     struct BlossomVariable {
   677       int begin, end;
   678       Value value;
   679 
   680       BlossomVariable(int _begin, int _end, Value _value)
   681         : begin(_begin), end(_end), value(_value) {}
   682 
   683     };
   684 
   685     typedef std::vector<BlossomVariable> BlossomPotential;
   686 
   687     const Graph& _graph;
   688     const WeightMap& _weight;
   689 
   690     MatchingMap* _matching;
   691 
   692     NodePotential* _node_potential;
   693 
   694     BlossomPotential _blossom_potential;
   695     BlossomNodeList _blossom_node_list;
   696 
   697     int _node_num;
   698     int _blossom_num;
   699 
   700     typedef RangeMap<int> IntIntMap;
   701 
   702     enum Status {
   703       EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
   704     };
   705 
   706     typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
   707     struct BlossomData {
   708       int tree;
   709       Status status;
   710       Arc pred, next;
   711       Value pot, offset;
   712       Node base;
   713     };
   714 
   715     IntNodeMap *_blossom_index;
   716     BlossomSet *_blossom_set;
   717     RangeMap<BlossomData>* _blossom_data;
   718 
   719     IntNodeMap *_node_index;
   720     IntArcMap *_node_heap_index;
   721 
   722     struct NodeData {
   723 
   724       NodeData(IntArcMap& node_heap_index)
   725         : heap(node_heap_index) {}
   726 
   727       int blossom;
   728       Value pot;
   729       BinHeap<Value, IntArcMap> heap;
   730       std::map<int, Arc> heap_index;
   731 
   732       int tree;
   733     };
   734 
   735     RangeMap<NodeData>* _node_data;
   736 
   737     typedef ExtendFindEnum<IntIntMap> TreeSet;
   738 
   739     IntIntMap *_tree_set_index;
   740     TreeSet *_tree_set;
   741 
   742     IntNodeMap *_delta1_index;
   743     BinHeap<Value, IntNodeMap> *_delta1;
   744 
   745     IntIntMap *_delta2_index;
   746     BinHeap<Value, IntIntMap> *_delta2;
   747 
   748     IntEdgeMap *_delta3_index;
   749     BinHeap<Value, IntEdgeMap> *_delta3;
   750 
   751     IntIntMap *_delta4_index;
   752     BinHeap<Value, IntIntMap> *_delta4;
   753 
   754     Value _delta_sum;
   755 
   756     void createStructures() {
   757       _node_num = countNodes(_graph);
   758       _blossom_num = _node_num * 3 / 2;
   759 
   760       if (!_matching) {
   761         _matching = new MatchingMap(_graph);
   762       }
   763       if (!_node_potential) {
   764         _node_potential = new NodePotential(_graph);
   765       }
   766       if (!_blossom_set) {
   767         _blossom_index = new IntNodeMap(_graph);
   768         _blossom_set = new BlossomSet(*_blossom_index);
   769         _blossom_data = new RangeMap<BlossomData>(_blossom_num);
   770       }
   771 
   772       if (!_node_index) {
   773         _node_index = new IntNodeMap(_graph);
   774         _node_heap_index = new IntArcMap(_graph);
   775         _node_data = new RangeMap<NodeData>(_node_num,
   776                                               NodeData(*_node_heap_index));
   777       }
   778 
   779       if (!_tree_set) {
   780         _tree_set_index = new IntIntMap(_blossom_num);
   781         _tree_set = new TreeSet(*_tree_set_index);
   782       }
   783       if (!_delta1) {
   784         _delta1_index = new IntNodeMap(_graph);
   785         _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
   786       }
   787       if (!_delta2) {
   788         _delta2_index = new IntIntMap(_blossom_num);
   789         _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
   790       }
   791       if (!_delta3) {
   792         _delta3_index = new IntEdgeMap(_graph);
   793         _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
   794       }
   795       if (!_delta4) {
   796         _delta4_index = new IntIntMap(_blossom_num);
   797         _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
   798       }
   799     }
   800 
   801     void destroyStructures() {
   802       _node_num = countNodes(_graph);
   803       _blossom_num = _node_num * 3 / 2;
   804 
   805       if (_matching) {
   806         delete _matching;
   807       }
   808       if (_node_potential) {
   809         delete _node_potential;
   810       }
   811       if (_blossom_set) {
   812         delete _blossom_index;
   813         delete _blossom_set;
   814         delete _blossom_data;
   815       }
   816 
   817       if (_node_index) {
   818         delete _node_index;
   819         delete _node_heap_index;
   820         delete _node_data;
   821       }
   822 
   823       if (_tree_set) {
   824         delete _tree_set_index;
   825         delete _tree_set;
   826       }
   827       if (_delta1) {
   828         delete _delta1_index;
   829         delete _delta1;
   830       }
   831       if (_delta2) {
   832         delete _delta2_index;
   833         delete _delta2;
   834       }
   835       if (_delta3) {
   836         delete _delta3_index;
   837         delete _delta3;
   838       }
   839       if (_delta4) {
   840         delete _delta4_index;
   841         delete _delta4;
   842       }
   843     }
   844 
   845     void matchedToEven(int blossom, int tree) {
   846       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
   847         _delta2->erase(blossom);
   848       }
   849 
   850       if (!_blossom_set->trivial(blossom)) {
   851         (*_blossom_data)[blossom].pot -=
   852           2 * (_delta_sum - (*_blossom_data)[blossom].offset);
   853       }
   854 
   855       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
   856            n != INVALID; ++n) {
   857 
   858         _blossom_set->increase(n, std::numeric_limits<Value>::max());
   859         int ni = (*_node_index)[n];
   860 
   861         (*_node_data)[ni].heap.clear();
   862         (*_node_data)[ni].heap_index.clear();
   863 
   864         (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
   865 
   866         _delta1->push(n, (*_node_data)[ni].pot);
   867 
   868         for (InArcIt e(_graph, n); e != INVALID; ++e) {
   869           Node v = _graph.source(e);
   870           int vb = _blossom_set->find(v);
   871           int vi = (*_node_index)[v];
   872 
   873           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   874             dualScale * _weight[e];
   875 
   876           if ((*_blossom_data)[vb].status == EVEN) {
   877             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
   878               _delta3->push(e, rw / 2);
   879             }
   880           } else if ((*_blossom_data)[vb].status == UNMATCHED) {
   881             if (_delta3->state(e) != _delta3->IN_HEAP) {
   882               _delta3->push(e, rw);
   883             }
   884           } else {
   885             typename std::map<int, Arc>::iterator it =
   886               (*_node_data)[vi].heap_index.find(tree);
   887 
   888             if (it != (*_node_data)[vi].heap_index.end()) {
   889               if ((*_node_data)[vi].heap[it->second] > rw) {
   890                 (*_node_data)[vi].heap.replace(it->second, e);
   891                 (*_node_data)[vi].heap.decrease(e, rw);
   892                 it->second = e;
   893               }
   894             } else {
   895               (*_node_data)[vi].heap.push(e, rw);
   896               (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
   897             }
   898 
   899             if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
   900               _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
   901 
   902               if ((*_blossom_data)[vb].status == MATCHED) {
   903                 if (_delta2->state(vb) != _delta2->IN_HEAP) {
   904                   _delta2->push(vb, _blossom_set->classPrio(vb) -
   905                                (*_blossom_data)[vb].offset);
   906                 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
   907                            (*_blossom_data)[vb].offset){
   908                   _delta2->decrease(vb, _blossom_set->classPrio(vb) -
   909                                    (*_blossom_data)[vb].offset);
   910                 }
   911               }
   912             }
   913           }
   914         }
   915       }
   916       (*_blossom_data)[blossom].offset = 0;
   917     }
   918 
   919     void matchedToOdd(int blossom) {
   920       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
   921         _delta2->erase(blossom);
   922       }
   923       (*_blossom_data)[blossom].offset += _delta_sum;
   924       if (!_blossom_set->trivial(blossom)) {
   925         _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
   926                      (*_blossom_data)[blossom].offset);
   927       }
   928     }
   929 
   930     void evenToMatched(int blossom, int tree) {
   931       if (!_blossom_set->trivial(blossom)) {
   932         (*_blossom_data)[blossom].pot += 2 * _delta_sum;
   933       }
   934 
   935       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
   936            n != INVALID; ++n) {
   937         int ni = (*_node_index)[n];
   938         (*_node_data)[ni].pot -= _delta_sum;
   939 
   940         _delta1->erase(n);
   941 
   942         for (InArcIt e(_graph, n); e != INVALID; ++e) {
   943           Node v = _graph.source(e);
   944           int vb = _blossom_set->find(v);
   945           int vi = (*_node_index)[v];
   946 
   947           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   948             dualScale * _weight[e];
   949 
   950           if (vb == blossom) {
   951             if (_delta3->state(e) == _delta3->IN_HEAP) {
   952               _delta3->erase(e);
   953             }
   954           } else if ((*_blossom_data)[vb].status == EVEN) {
   955 
   956             if (_delta3->state(e) == _delta3->IN_HEAP) {
   957               _delta3->erase(e);
   958             }
   959 
   960             int vt = _tree_set->find(vb);
   961 
   962             if (vt != tree) {
   963 
   964               Arc r = _graph.oppositeArc(e);
   965 
   966               typename std::map<int, Arc>::iterator it =
   967                 (*_node_data)[ni].heap_index.find(vt);
   968 
   969               if (it != (*_node_data)[ni].heap_index.end()) {
   970                 if ((*_node_data)[ni].heap[it->second] > rw) {
   971                   (*_node_data)[ni].heap.replace(it->second, r);
   972                   (*_node_data)[ni].heap.decrease(r, rw);
   973                   it->second = r;
   974                 }
   975               } else {
   976                 (*_node_data)[ni].heap.push(r, rw);
   977                 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
   978               }
   979 
   980               if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
   981                 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
   982 
   983                 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
   984                   _delta2->push(blossom, _blossom_set->classPrio(blossom) -
   985                                (*_blossom_data)[blossom].offset);
   986                 } else if ((*_delta2)[blossom] >
   987                            _blossom_set->classPrio(blossom) -
   988                            (*_blossom_data)[blossom].offset){
   989                   _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
   990                                    (*_blossom_data)[blossom].offset);
   991                 }
   992               }
   993             }
   994 
   995           } else if ((*_blossom_data)[vb].status == UNMATCHED) {
   996             if (_delta3->state(e) == _delta3->IN_HEAP) {
   997               _delta3->erase(e);
   998             }
   999           } else {
  1000 
  1001             typename std::map<int, Arc>::iterator it =
  1002               (*_node_data)[vi].heap_index.find(tree);
  1003 
  1004             if (it != (*_node_data)[vi].heap_index.end()) {
  1005               (*_node_data)[vi].heap.erase(it->second);
  1006               (*_node_data)[vi].heap_index.erase(it);
  1007               if ((*_node_data)[vi].heap.empty()) {
  1008                 _blossom_set->increase(v, std::numeric_limits<Value>::max());
  1009               } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
  1010                 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
  1011               }
  1012 
  1013               if ((*_blossom_data)[vb].status == MATCHED) {
  1014                 if (_blossom_set->classPrio(vb) ==
  1015                     std::numeric_limits<Value>::max()) {
  1016                   _delta2->erase(vb);
  1017                 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
  1018                            (*_blossom_data)[vb].offset) {
  1019                   _delta2->increase(vb, _blossom_set->classPrio(vb) -
  1020                                    (*_blossom_data)[vb].offset);
  1021                 }
  1022               }
  1023             }
  1024           }
  1025         }
  1026       }
  1027     }
  1028 
  1029     void oddToMatched(int blossom) {
  1030       (*_blossom_data)[blossom].offset -= _delta_sum;
  1031 
  1032       if (_blossom_set->classPrio(blossom) !=
  1033           std::numeric_limits<Value>::max()) {
  1034         _delta2->push(blossom, _blossom_set->classPrio(blossom) -
  1035                        (*_blossom_data)[blossom].offset);
  1036       }
  1037 
  1038       if (!_blossom_set->trivial(blossom)) {
  1039         _delta4->erase(blossom);
  1040       }
  1041     }
  1042 
  1043     void oddToEven(int blossom, int tree) {
  1044       if (!_blossom_set->trivial(blossom)) {
  1045         _delta4->erase(blossom);
  1046         (*_blossom_data)[blossom].pot -=
  1047           2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
  1048       }
  1049 
  1050       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  1051            n != INVALID; ++n) {
  1052         int ni = (*_node_index)[n];
  1053 
  1054         _blossom_set->increase(n, std::numeric_limits<Value>::max());
  1055 
  1056         (*_node_data)[ni].heap.clear();
  1057         (*_node_data)[ni].heap_index.clear();
  1058         (*_node_data)[ni].pot +=
  1059           2 * _delta_sum - (*_blossom_data)[blossom].offset;
  1060 
  1061         _delta1->push(n, (*_node_data)[ni].pot);
  1062 
  1063         for (InArcIt e(_graph, n); e != INVALID; ++e) {
  1064           Node v = _graph.source(e);
  1065           int vb = _blossom_set->find(v);
  1066           int vi = (*_node_index)[v];
  1067 
  1068           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  1069             dualScale * _weight[e];
  1070 
  1071           if ((*_blossom_data)[vb].status == EVEN) {
  1072             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
  1073               _delta3->push(e, rw / 2);
  1074             }
  1075           } else if ((*_blossom_data)[vb].status == UNMATCHED) {
  1076             if (_delta3->state(e) != _delta3->IN_HEAP) {
  1077               _delta3->push(e, rw);
  1078             }
  1079           } else {
  1080 
  1081             typename std::map<int, Arc>::iterator it =
  1082               (*_node_data)[vi].heap_index.find(tree);
  1083 
  1084             if (it != (*_node_data)[vi].heap_index.end()) {
  1085               if ((*_node_data)[vi].heap[it->second] > rw) {
  1086                 (*_node_data)[vi].heap.replace(it->second, e);
  1087                 (*_node_data)[vi].heap.decrease(e, rw);
  1088                 it->second = e;
  1089               }
  1090             } else {
  1091               (*_node_data)[vi].heap.push(e, rw);
  1092               (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
  1093             }
  1094 
  1095             if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
  1096               _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
  1097 
  1098               if ((*_blossom_data)[vb].status == MATCHED) {
  1099                 if (_delta2->state(vb) != _delta2->IN_HEAP) {
  1100                   _delta2->push(vb, _blossom_set->classPrio(vb) -
  1101                                (*_blossom_data)[vb].offset);
  1102                 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
  1103                            (*_blossom_data)[vb].offset) {
  1104                   _delta2->decrease(vb, _blossom_set->classPrio(vb) -
  1105                                    (*_blossom_data)[vb].offset);
  1106                 }
  1107               }
  1108             }
  1109           }
  1110         }
  1111       }
  1112       (*_blossom_data)[blossom].offset = 0;
  1113     }
  1114 
  1115 
  1116     void matchedToUnmatched(int blossom) {
  1117       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
  1118         _delta2->erase(blossom);
  1119       }
  1120 
  1121       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  1122            n != INVALID; ++n) {
  1123         int ni = (*_node_index)[n];
  1124 
  1125         _blossom_set->increase(n, std::numeric_limits<Value>::max());
  1126 
  1127         (*_node_data)[ni].heap.clear();
  1128         (*_node_data)[ni].heap_index.clear();
  1129 
  1130         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  1131           Node v = _graph.target(e);
  1132           int vb = _blossom_set->find(v);
  1133           int vi = (*_node_index)[v];
  1134 
  1135           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  1136             dualScale * _weight[e];
  1137 
  1138           if ((*_blossom_data)[vb].status == EVEN) {
  1139             if (_delta3->state(e) != _delta3->IN_HEAP) {
  1140               _delta3->push(e, rw);
  1141             }
  1142           }
  1143         }
  1144       }
  1145     }
  1146 
  1147     void unmatchedToMatched(int blossom) {
  1148       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  1149            n != INVALID; ++n) {
  1150         int ni = (*_node_index)[n];
  1151 
  1152         for (InArcIt e(_graph, n); e != INVALID; ++e) {
  1153           Node v = _graph.source(e);
  1154           int vb = _blossom_set->find(v);
  1155           int vi = (*_node_index)[v];
  1156 
  1157           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  1158             dualScale * _weight[e];
  1159 
  1160           if (vb == blossom) {
  1161             if (_delta3->state(e) == _delta3->IN_HEAP) {
  1162               _delta3->erase(e);
  1163             }
  1164           } else if ((*_blossom_data)[vb].status == EVEN) {
  1165 
  1166             if (_delta3->state(e) == _delta3->IN_HEAP) {
  1167               _delta3->erase(e);
  1168             }
  1169 
  1170             int vt = _tree_set->find(vb);
  1171 
  1172             Arc r = _graph.oppositeArc(e);
  1173 
  1174             typename std::map<int, Arc>::iterator it =
  1175               (*_node_data)[ni].heap_index.find(vt);
  1176 
  1177             if (it != (*_node_data)[ni].heap_index.end()) {
  1178               if ((*_node_data)[ni].heap[it->second] > rw) {
  1179                 (*_node_data)[ni].heap.replace(it->second, r);
  1180                 (*_node_data)[ni].heap.decrease(r, rw);
  1181                 it->second = r;
  1182               }
  1183             } else {
  1184               (*_node_data)[ni].heap.push(r, rw);
  1185               (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
  1186             }
  1187 
  1188             if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
  1189               _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
  1190 
  1191               if (_delta2->state(blossom) != _delta2->IN_HEAP) {
  1192                 _delta2->push(blossom, _blossom_set->classPrio(blossom) -
  1193                              (*_blossom_data)[blossom].offset);
  1194               } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
  1195                          (*_blossom_data)[blossom].offset){
  1196                 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
  1197                                  (*_blossom_data)[blossom].offset);
  1198               }
  1199             }
  1200 
  1201           } else if ((*_blossom_data)[vb].status == UNMATCHED) {
  1202             if (_delta3->state(e) == _delta3->IN_HEAP) {
  1203               _delta3->erase(e);
  1204             }
  1205           }
  1206         }
  1207       }
  1208     }
  1209 
  1210     void alternatePath(int even, int tree) {
  1211       int odd;
  1212 
  1213       evenToMatched(even, tree);
  1214       (*_blossom_data)[even].status = MATCHED;
  1215 
  1216       while ((*_blossom_data)[even].pred != INVALID) {
  1217         odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
  1218         (*_blossom_data)[odd].status = MATCHED;
  1219         oddToMatched(odd);
  1220         (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
  1221 
  1222         even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
  1223         (*_blossom_data)[even].status = MATCHED;
  1224         evenToMatched(even, tree);
  1225         (*_blossom_data)[even].next =
  1226           _graph.oppositeArc((*_blossom_data)[odd].pred);
  1227       }
  1228 
  1229     }
  1230 
  1231     void destroyTree(int tree) {
  1232       for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
  1233         if ((*_blossom_data)[b].status == EVEN) {
  1234           (*_blossom_data)[b].status = MATCHED;
  1235           evenToMatched(b, tree);
  1236         } else if ((*_blossom_data)[b].status == ODD) {
  1237           (*_blossom_data)[b].status = MATCHED;
  1238           oddToMatched(b);
  1239         }
  1240       }
  1241       _tree_set->eraseClass(tree);
  1242     }
  1243 
  1244 
  1245     void unmatchNode(const Node& node) {
  1246       int blossom = _blossom_set->find(node);
  1247       int tree = _tree_set->find(blossom);
  1248 
  1249       alternatePath(blossom, tree);
  1250       destroyTree(tree);
  1251 
  1252       (*_blossom_data)[blossom].status = UNMATCHED;
  1253       (*_blossom_data)[blossom].base = node;
  1254       matchedToUnmatched(blossom);
  1255     }
  1256 
  1257 
  1258     void augmentOnEdge(const Edge& edge) {
  1259 
  1260       int left = _blossom_set->find(_graph.u(edge));
  1261       int right = _blossom_set->find(_graph.v(edge));
  1262 
  1263       if ((*_blossom_data)[left].status == EVEN) {
  1264         int left_tree = _tree_set->find(left);
  1265         alternatePath(left, left_tree);
  1266         destroyTree(left_tree);
  1267       } else {
  1268         (*_blossom_data)[left].status = MATCHED;
  1269         unmatchedToMatched(left);
  1270       }
  1271 
  1272       if ((*_blossom_data)[right].status == EVEN) {
  1273         int right_tree = _tree_set->find(right);
  1274         alternatePath(right, right_tree);
  1275         destroyTree(right_tree);
  1276       } else {
  1277         (*_blossom_data)[right].status = MATCHED;
  1278         unmatchedToMatched(right);
  1279       }
  1280 
  1281       (*_blossom_data)[left].next = _graph.direct(edge, true);
  1282       (*_blossom_data)[right].next = _graph.direct(edge, false);
  1283     }
  1284 
  1285     void extendOnArc(const Arc& arc) {
  1286       int base = _blossom_set->find(_graph.target(arc));
  1287       int tree = _tree_set->find(base);
  1288 
  1289       int odd = _blossom_set->find(_graph.source(arc));
  1290       _tree_set->insert(odd, tree);
  1291       (*_blossom_data)[odd].status = ODD;
  1292       matchedToOdd(odd);
  1293       (*_blossom_data)[odd].pred = arc;
  1294 
  1295       int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
  1296       (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
  1297       _tree_set->insert(even, tree);
  1298       (*_blossom_data)[even].status = EVEN;
  1299       matchedToEven(even, tree);
  1300     }
  1301 
  1302     void shrinkOnEdge(const Edge& edge, int tree) {
  1303       int nca = -1;
  1304       std::vector<int> left_path, right_path;
  1305 
  1306       {
  1307         std::set<int> left_set, right_set;
  1308         int left = _blossom_set->find(_graph.u(edge));
  1309         left_path.push_back(left);
  1310         left_set.insert(left);
  1311 
  1312         int right = _blossom_set->find(_graph.v(edge));
  1313         right_path.push_back(right);
  1314         right_set.insert(right);
  1315 
  1316         while (true) {
  1317 
  1318           if ((*_blossom_data)[left].pred == INVALID) break;
  1319 
  1320           left =
  1321             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
  1322           left_path.push_back(left);
  1323           left =
  1324             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
  1325           left_path.push_back(left);
  1326 
  1327           left_set.insert(left);
  1328 
  1329           if (right_set.find(left) != right_set.end()) {
  1330             nca = left;
  1331             break;
  1332           }
  1333 
  1334           if ((*_blossom_data)[right].pred == INVALID) break;
  1335 
  1336           right =
  1337             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
  1338           right_path.push_back(right);
  1339           right =
  1340             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
  1341           right_path.push_back(right);
  1342 
  1343           right_set.insert(right);
  1344 
  1345           if (left_set.find(right) != left_set.end()) {
  1346             nca = right;
  1347             break;
  1348           }
  1349 
  1350         }
  1351 
  1352         if (nca == -1) {
  1353           if ((*_blossom_data)[left].pred == INVALID) {
  1354             nca = right;
  1355             while (left_set.find(nca) == left_set.end()) {
  1356               nca =
  1357                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  1358               right_path.push_back(nca);
  1359               nca =
  1360                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  1361               right_path.push_back(nca);
  1362             }
  1363           } else {
  1364             nca = left;
  1365             while (right_set.find(nca) == right_set.end()) {
  1366               nca =
  1367                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  1368               left_path.push_back(nca);
  1369               nca =
  1370                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  1371               left_path.push_back(nca);
  1372             }
  1373           }
  1374         }
  1375       }
  1376 
  1377       std::vector<int> subblossoms;
  1378       Arc prev;
  1379 
  1380       prev = _graph.direct(edge, true);
  1381       for (int i = 0; left_path[i] != nca; i += 2) {
  1382         subblossoms.push_back(left_path[i]);
  1383         (*_blossom_data)[left_path[i]].next = prev;
  1384         _tree_set->erase(left_path[i]);
  1385 
  1386         subblossoms.push_back(left_path[i + 1]);
  1387         (*_blossom_data)[left_path[i + 1]].status = EVEN;
  1388         oddToEven(left_path[i + 1], tree);
  1389         _tree_set->erase(left_path[i + 1]);
  1390         prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
  1391       }
  1392 
  1393       int k = 0;
  1394       while (right_path[k] != nca) ++k;
  1395 
  1396       subblossoms.push_back(nca);
  1397       (*_blossom_data)[nca].next = prev;
  1398 
  1399       for (int i = k - 2; i >= 0; i -= 2) {
  1400         subblossoms.push_back(right_path[i + 1]);
  1401         (*_blossom_data)[right_path[i + 1]].status = EVEN;
  1402         oddToEven(right_path[i + 1], tree);
  1403         _tree_set->erase(right_path[i + 1]);
  1404 
  1405         (*_blossom_data)[right_path[i + 1]].next =
  1406           (*_blossom_data)[right_path[i + 1]].pred;
  1407 
  1408         subblossoms.push_back(right_path[i]);
  1409         _tree_set->erase(right_path[i]);
  1410       }
  1411 
  1412       int surface =
  1413         _blossom_set->join(subblossoms.begin(), subblossoms.end());
  1414 
  1415       for (int i = 0; i < int(subblossoms.size()); ++i) {
  1416         if (!_blossom_set->trivial(subblossoms[i])) {
  1417           (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
  1418         }
  1419         (*_blossom_data)[subblossoms[i]].status = MATCHED;
  1420       }
  1421 
  1422       (*_blossom_data)[surface].pot = -2 * _delta_sum;
  1423       (*_blossom_data)[surface].offset = 0;
  1424       (*_blossom_data)[surface].status = EVEN;
  1425       (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
  1426       (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
  1427 
  1428       _tree_set->insert(surface, tree);
  1429       _tree_set->erase(nca);
  1430     }
  1431 
  1432     void splitBlossom(int blossom) {
  1433       Arc next = (*_blossom_data)[blossom].next;
  1434       Arc pred = (*_blossom_data)[blossom].pred;
  1435 
  1436       int tree = _tree_set->find(blossom);
  1437 
  1438       (*_blossom_data)[blossom].status = MATCHED;
  1439       oddToMatched(blossom);
  1440       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
  1441         _delta2->erase(blossom);
  1442       }
  1443 
  1444       std::vector<int> subblossoms;
  1445       _blossom_set->split(blossom, std::back_inserter(subblossoms));
  1446 
  1447       Value offset = (*_blossom_data)[blossom].offset;
  1448       int b = _blossom_set->find(_graph.source(pred));
  1449       int d = _blossom_set->find(_graph.source(next));
  1450 
  1451       int ib = -1, id = -1;
  1452       for (int i = 0; i < int(subblossoms.size()); ++i) {
  1453         if (subblossoms[i] == b) ib = i;
  1454         if (subblossoms[i] == d) id = i;
  1455 
  1456         (*_blossom_data)[subblossoms[i]].offset = offset;
  1457         if (!_blossom_set->trivial(subblossoms[i])) {
  1458           (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
  1459         }
  1460         if (_blossom_set->classPrio(subblossoms[i]) !=
  1461             std::numeric_limits<Value>::max()) {
  1462           _delta2->push(subblossoms[i],
  1463                         _blossom_set->classPrio(subblossoms[i]) -
  1464                         (*_blossom_data)[subblossoms[i]].offset);
  1465         }
  1466       }
  1467 
  1468       if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
  1469         for (int i = (id + 1) % subblossoms.size();
  1470              i != ib; i = (i + 2) % subblossoms.size()) {
  1471           int sb = subblossoms[i];
  1472           int tb = subblossoms[(i + 1) % subblossoms.size()];
  1473           (*_blossom_data)[sb].next =
  1474             _graph.oppositeArc((*_blossom_data)[tb].next);
  1475         }
  1476 
  1477         for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
  1478           int sb = subblossoms[i];
  1479           int tb = subblossoms[(i + 1) % subblossoms.size()];
  1480           int ub = subblossoms[(i + 2) % subblossoms.size()];
  1481 
  1482           (*_blossom_data)[sb].status = ODD;
  1483           matchedToOdd(sb);
  1484           _tree_set->insert(sb, tree);
  1485           (*_blossom_data)[sb].pred = pred;
  1486           (*_blossom_data)[sb].next =
  1487                            _graph.oppositeArc((*_blossom_data)[tb].next);
  1488 
  1489           pred = (*_blossom_data)[ub].next;
  1490 
  1491           (*_blossom_data)[tb].status = EVEN;
  1492           matchedToEven(tb, tree);
  1493           _tree_set->insert(tb, tree);
  1494           (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
  1495         }
  1496 
  1497         (*_blossom_data)[subblossoms[id]].status = ODD;
  1498         matchedToOdd(subblossoms[id]);
  1499         _tree_set->insert(subblossoms[id], tree);
  1500         (*_blossom_data)[subblossoms[id]].next = next;
  1501         (*_blossom_data)[subblossoms[id]].pred = pred;
  1502 
  1503       } else {
  1504 
  1505         for (int i = (ib + 1) % subblossoms.size();
  1506              i != id; i = (i + 2) % subblossoms.size()) {
  1507           int sb = subblossoms[i];
  1508           int tb = subblossoms[(i + 1) % subblossoms.size()];
  1509           (*_blossom_data)[sb].next =
  1510             _graph.oppositeArc((*_blossom_data)[tb].next);
  1511         }
  1512 
  1513         for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
  1514           int sb = subblossoms[i];
  1515           int tb = subblossoms[(i + 1) % subblossoms.size()];
  1516           int ub = subblossoms[(i + 2) % subblossoms.size()];
  1517 
  1518           (*_blossom_data)[sb].status = ODD;
  1519           matchedToOdd(sb);
  1520           _tree_set->insert(sb, tree);
  1521           (*_blossom_data)[sb].next = next;
  1522           (*_blossom_data)[sb].pred =
  1523             _graph.oppositeArc((*_blossom_data)[tb].next);
  1524 
  1525           (*_blossom_data)[tb].status = EVEN;
  1526           matchedToEven(tb, tree);
  1527           _tree_set->insert(tb, tree);
  1528           (*_blossom_data)[tb].pred =
  1529             (*_blossom_data)[tb].next =
  1530             _graph.oppositeArc((*_blossom_data)[ub].next);
  1531           next = (*_blossom_data)[ub].next;
  1532         }
  1533 
  1534         (*_blossom_data)[subblossoms[ib]].status = ODD;
  1535         matchedToOdd(subblossoms[ib]);
  1536         _tree_set->insert(subblossoms[ib], tree);
  1537         (*_blossom_data)[subblossoms[ib]].next = next;
  1538         (*_blossom_data)[subblossoms[ib]].pred = pred;
  1539       }
  1540       _tree_set->erase(blossom);
  1541     }
  1542 
  1543     void extractBlossom(int blossom, const Node& base, const Arc& matching) {
  1544       if (_blossom_set->trivial(blossom)) {
  1545         int bi = (*_node_index)[base];
  1546         Value pot = (*_node_data)[bi].pot;
  1547 
  1548         _matching->set(base, matching);
  1549         _blossom_node_list.push_back(base);
  1550         _node_potential->set(base, pot);
  1551       } else {
  1552 
  1553         Value pot = (*_blossom_data)[blossom].pot;
  1554         int bn = _blossom_node_list.size();
  1555 
  1556         std::vector<int> subblossoms;
  1557         _blossom_set->split(blossom, std::back_inserter(subblossoms));
  1558         int b = _blossom_set->find(base);
  1559         int ib = -1;
  1560         for (int i = 0; i < int(subblossoms.size()); ++i) {
  1561           if (subblossoms[i] == b) { ib = i; break; }
  1562         }
  1563 
  1564         for (int i = 1; i < int(subblossoms.size()); i += 2) {
  1565           int sb = subblossoms[(ib + i) % subblossoms.size()];
  1566           int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
  1567 
  1568           Arc m = (*_blossom_data)[tb].next;
  1569           extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
  1570           extractBlossom(tb, _graph.source(m), m);
  1571         }
  1572         extractBlossom(subblossoms[ib], base, matching);
  1573 
  1574         int en = _blossom_node_list.size();
  1575 
  1576         _blossom_potential.push_back(BlossomVariable(bn, en, pot));
  1577       }
  1578     }
  1579 
  1580     void extractMatching() {
  1581       std::vector<int> blossoms;
  1582       for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
  1583         blossoms.push_back(c);
  1584       }
  1585 
  1586       for (int i = 0; i < int(blossoms.size()); ++i) {
  1587         if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
  1588 
  1589           Value offset = (*_blossom_data)[blossoms[i]].offset;
  1590           (*_blossom_data)[blossoms[i]].pot += 2 * offset;
  1591           for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
  1592                n != INVALID; ++n) {
  1593             (*_node_data)[(*_node_index)[n]].pot -= offset;
  1594           }
  1595 
  1596           Arc matching = (*_blossom_data)[blossoms[i]].next;
  1597           Node base = _graph.source(matching);
  1598           extractBlossom(blossoms[i], base, matching);
  1599         } else {
  1600           Node base = (*_blossom_data)[blossoms[i]].base;
  1601           extractBlossom(blossoms[i], base, INVALID);
  1602         }
  1603       }
  1604     }
  1605 
  1606   public:
  1607 
  1608     /// \brief Constructor
  1609     ///
  1610     /// Constructor.
  1611     MaxWeightedMatching(const Graph& graph, const WeightMap& weight)
  1612       : _graph(graph), _weight(weight), _matching(0),
  1613         _node_potential(0), _blossom_potential(), _blossom_node_list(),
  1614         _node_num(0), _blossom_num(0),
  1615 
  1616         _blossom_index(0), _blossom_set(0), _blossom_data(0),
  1617         _node_index(0), _node_heap_index(0), _node_data(0),
  1618         _tree_set_index(0), _tree_set(0),
  1619 
  1620         _delta1_index(0), _delta1(0),
  1621         _delta2_index(0), _delta2(0),
  1622         _delta3_index(0), _delta3(0),
  1623         _delta4_index(0), _delta4(0),
  1624 
  1625         _delta_sum() {}
  1626 
  1627     ~MaxWeightedMatching() {
  1628       destroyStructures();
  1629     }
  1630 
  1631     /// \name Execution control
  1632     /// The simplest way to execute the algorithm is to use the
  1633     /// \c run() member function.
  1634 
  1635     ///@{
  1636 
  1637     /// \brief Initialize the algorithm
  1638     ///
  1639     /// Initialize the algorithm
  1640     void init() {
  1641       createStructures();
  1642 
  1643       for (ArcIt e(_graph); e != INVALID; ++e) {
  1644         _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
  1645       }
  1646       for (NodeIt n(_graph); n != INVALID; ++n) {
  1647         _delta1_index->set(n, _delta1->PRE_HEAP);
  1648       }
  1649       for (EdgeIt e(_graph); e != INVALID; ++e) {
  1650         _delta3_index->set(e, _delta3->PRE_HEAP);
  1651       }
  1652       for (int i = 0; i < _blossom_num; ++i) {
  1653         _delta2_index->set(i, _delta2->PRE_HEAP);
  1654         _delta4_index->set(i, _delta4->PRE_HEAP);
  1655       }
  1656 
  1657       int index = 0;
  1658       for (NodeIt n(_graph); n != INVALID; ++n) {
  1659         Value max = 0;
  1660         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  1661           if (_graph.target(e) == n) continue;
  1662           if ((dualScale * _weight[e]) / 2 > max) {
  1663             max = (dualScale * _weight[e]) / 2;
  1664           }
  1665         }
  1666         _node_index->set(n, index);
  1667         (*_node_data)[index].pot = max;
  1668         _delta1->push(n, max);
  1669         int blossom =
  1670           _blossom_set->insert(n, std::numeric_limits<Value>::max());
  1671 
  1672         _tree_set->insert(blossom);
  1673 
  1674         (*_blossom_data)[blossom].status = EVEN;
  1675         (*_blossom_data)[blossom].pred = INVALID;
  1676         (*_blossom_data)[blossom].next = INVALID;
  1677         (*_blossom_data)[blossom].pot = 0;
  1678         (*_blossom_data)[blossom].offset = 0;
  1679         ++index;
  1680       }
  1681       for (EdgeIt e(_graph); e != INVALID; ++e) {
  1682         int si = (*_node_index)[_graph.u(e)];
  1683         int ti = (*_node_index)[_graph.v(e)];
  1684         if (_graph.u(e) != _graph.v(e)) {
  1685           _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
  1686                             dualScale * _weight[e]) / 2);
  1687         }
  1688       }
  1689     }
  1690 
  1691     /// \brief Starts the algorithm
  1692     ///
  1693     /// Starts the algorithm
  1694     void start() {
  1695       enum OpType {
  1696         D1, D2, D3, D4
  1697       };
  1698 
  1699       int unmatched = _node_num;
  1700       while (unmatched > 0) {
  1701         Value d1 = !_delta1->empty() ?
  1702           _delta1->prio() : std::numeric_limits<Value>::max();
  1703 
  1704         Value d2 = !_delta2->empty() ?
  1705           _delta2->prio() : std::numeric_limits<Value>::max();
  1706 
  1707         Value d3 = !_delta3->empty() ?
  1708           _delta3->prio() : std::numeric_limits<Value>::max();
  1709 
  1710         Value d4 = !_delta4->empty() ?
  1711           _delta4->prio() : std::numeric_limits<Value>::max();
  1712 
  1713         _delta_sum = d1; OpType ot = D1;
  1714         if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
  1715         if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
  1716         if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
  1717 
  1718 
  1719         switch (ot) {
  1720         case D1:
  1721           {
  1722             Node n = _delta1->top();
  1723             unmatchNode(n);
  1724             --unmatched;
  1725           }
  1726           break;
  1727         case D2:
  1728           {
  1729             int blossom = _delta2->top();
  1730             Node n = _blossom_set->classTop(blossom);
  1731             Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
  1732             extendOnArc(e);
  1733           }
  1734           break;
  1735         case D3:
  1736           {
  1737             Edge e = _delta3->top();
  1738 
  1739             int left_blossom = _blossom_set->find(_graph.u(e));
  1740             int right_blossom = _blossom_set->find(_graph.v(e));
  1741 
  1742             if (left_blossom == right_blossom) {
  1743               _delta3->pop();
  1744             } else {
  1745               int left_tree;
  1746               if ((*_blossom_data)[left_blossom].status == EVEN) {
  1747                 left_tree = _tree_set->find(left_blossom);
  1748               } else {
  1749                 left_tree = -1;
  1750                 ++unmatched;
  1751               }
  1752               int right_tree;
  1753               if ((*_blossom_data)[right_blossom].status == EVEN) {
  1754                 right_tree = _tree_set->find(right_blossom);
  1755               } else {
  1756                 right_tree = -1;
  1757                 ++unmatched;
  1758               }
  1759 
  1760               if (left_tree == right_tree) {
  1761                 shrinkOnEdge(e, left_tree);
  1762               } else {
  1763                 augmentOnEdge(e);
  1764                 unmatched -= 2;
  1765               }
  1766             }
  1767           } break;
  1768         case D4:
  1769           splitBlossom(_delta4->top());
  1770           break;
  1771         }
  1772       }
  1773       extractMatching();
  1774     }
  1775 
  1776     /// \brief Runs %MaxWeightedMatching algorithm.
  1777     ///
  1778     /// This method runs the %MaxWeightedMatching algorithm.
  1779     ///
  1780     /// \note mwm.run() is just a shortcut of the following code.
  1781     /// \code
  1782     ///   mwm.init();
  1783     ///   mwm.start();
  1784     /// \endcode
  1785     void run() {
  1786       init();
  1787       start();
  1788     }
  1789 
  1790     /// @}
  1791 
  1792     /// \name Primal solution
  1793     /// Functions to get the primal solution, ie. the matching.
  1794 
  1795     /// @{
  1796 
  1797     /// \brief Returns the weight of the matching.
  1798     ///
  1799     /// Returns the weight of the matching.
  1800     Value matchingValue() const {
  1801       Value sum = 0;
  1802       for (NodeIt n(_graph); n != INVALID; ++n) {
  1803         if ((*_matching)[n] != INVALID) {
  1804           sum += _weight[(*_matching)[n]];
  1805         }
  1806       }
  1807       return sum /= 2;
  1808     }
  1809 
  1810     /// \brief Returns the cardinality of the matching.
  1811     ///
  1812     /// Returns the cardinality of the matching.
  1813     int matchingSize() const {
  1814       int num = 0;
  1815       for (NodeIt n(_graph); n != INVALID; ++n) {
  1816         if ((*_matching)[n] != INVALID) {
  1817           ++num;
  1818         }
  1819       }
  1820       return num /= 2;
  1821     }
  1822 
  1823     /// \brief Returns true when the edge is in the matching.
  1824     ///
  1825     /// Returns true when the edge is in the matching.
  1826     bool matching(const Edge& edge) const {
  1827       return edge == (*_matching)[_graph.u(edge)];
  1828     }
  1829 
  1830     /// \brief Returns the incident matching arc.
  1831     ///
  1832     /// Returns the incident matching arc from given node. If the
  1833     /// node is not matched then it gives back \c INVALID.
  1834     Arc matching(const Node& node) const {
  1835       return (*_matching)[node];
  1836     }
  1837 
  1838     /// \brief Returns the mate of the node.
  1839     ///
  1840     /// Returns the adjancent node in a mathcing arc. If the node is
  1841     /// not matched then it gives back \c INVALID.
  1842     Node mate(const Node& node) const {
  1843       return (*_matching)[node] != INVALID ?
  1844         _graph.target((*_matching)[node]) : INVALID;
  1845     }
  1846 
  1847     /// @}
  1848 
  1849     /// \name Dual solution
  1850     /// Functions to get the dual solution.
  1851 
  1852     /// @{
  1853 
  1854     /// \brief Returns the value of the dual solution.
  1855     ///
  1856     /// Returns the value of the dual solution. It should be equal to
  1857     /// the primal value scaled by \ref dualScale "dual scale".
  1858     Value dualValue() const {
  1859       Value sum = 0;
  1860       for (NodeIt n(_graph); n != INVALID; ++n) {
  1861         sum += nodeValue(n);
  1862       }
  1863       for (int i = 0; i < blossomNum(); ++i) {
  1864         sum += blossomValue(i) * (blossomSize(i) / 2);
  1865       }
  1866       return sum;
  1867     }
  1868 
  1869     /// \brief Returns the value of the node.
  1870     ///
  1871     /// Returns the the value of the node.
  1872     Value nodeValue(const Node& n) const {
  1873       return (*_node_potential)[n];
  1874     }
  1875 
  1876     /// \brief Returns the number of the blossoms in the basis.
  1877     ///
  1878     /// Returns the number of the blossoms in the basis.
  1879     /// \see BlossomIt
  1880     int blossomNum() const {
  1881       return _blossom_potential.size();
  1882     }
  1883 
  1884 
  1885     /// \brief Returns the number of the nodes in the blossom.
  1886     ///
  1887     /// Returns the number of the nodes in the blossom.
  1888     int blossomSize(int k) const {
  1889       return _blossom_potential[k].end - _blossom_potential[k].begin;
  1890     }
  1891 
  1892     /// \brief Returns the value of the blossom.
  1893     ///
  1894     /// Returns the the value of the blossom.
  1895     /// \see BlossomIt
  1896     Value blossomValue(int k) const {
  1897       return _blossom_potential[k].value;
  1898     }
  1899 
  1900     /// \brief Iterator for obtaining the nodes of the blossom.
  1901     ///
  1902     /// Iterator for obtaining the nodes of the blossom. This class
  1903     /// provides a common lemon style iterator for listing a
  1904     /// subset of the nodes.
  1905     class BlossomIt {
  1906     public:
  1907 
  1908       /// \brief Constructor.
  1909       ///
  1910       /// Constructor to get the nodes of the variable.
  1911       BlossomIt(const MaxWeightedMatching& algorithm, int variable)
  1912         : _algorithm(&algorithm)
  1913       {
  1914         _index = _algorithm->_blossom_potential[variable].begin;
  1915         _last = _algorithm->_blossom_potential[variable].end;
  1916       }
  1917 
  1918       /// \brief Conversion to node.
  1919       ///
  1920       /// Conversion to node.
  1921       operator Node() const {
  1922         return _algorithm->_blossom_node_list[_index];
  1923       }
  1924 
  1925       /// \brief Increment operator.
  1926       ///
  1927       /// Increment operator.
  1928       BlossomIt& operator++() {
  1929         ++_index;
  1930         return *this;
  1931       }
  1932 
  1933       /// \brief Validity checking
  1934       ///
  1935       /// Checks whether the iterator is invalid.
  1936       bool operator==(Invalid) const { return _index == _last; }
  1937 
  1938       /// \brief Validity checking
  1939       ///
  1940       /// Checks whether the iterator is valid.
  1941       bool operator!=(Invalid) const { return _index != _last; }
  1942 
  1943     private:
  1944       const MaxWeightedMatching* _algorithm;
  1945       int _last;
  1946       int _index;
  1947     };
  1948 
  1949     /// @}
  1950 
  1951   };
  1952 
  1953   /// \ingroup matching
  1954   ///
  1955   /// \brief Weighted perfect matching in general graphs
  1956   ///
  1957   /// This class provides an efficient implementation of Edmond's
  1958   /// maximum weighted perfect matching algorithm. The implementation
  1959   /// is based on extensive use of priority queues and provides
  1960   /// \f$O(nm\log(n))\f$ time complexity.
  1961   ///
  1962   /// The maximum weighted matching problem is to find undirected
  1963   /// edges in the graph with maximum overall weight and no two of
  1964   /// them shares their ends and covers all nodes. The problem can be
  1965   /// formulated with the following linear program.
  1966   /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
  1967   /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
  1968       \quad \forall B\in\mathcal{O}\f] */
  1969   /// \f[x_e \ge 0\quad \forall e\in E\f]
  1970   /// \f[\max \sum_{e\in E}x_ew_e\f]
  1971   /// where \f$\delta(X)\f$ is the set of edges incident to a node in
  1972   /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
  1973   /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
  1974   /// subsets of the nodes.
  1975   ///
  1976   /// The algorithm calculates an optimal matching and a proof of the
  1977   /// optimality. The solution of the dual problem can be used to check
  1978   /// the result of the algorithm. The dual linear problem is the
  1979   /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
  1980       w_{uv} \quad \forall uv\in E\f] */
  1981   /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
  1982   /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
  1983       \frac{\vert B \vert - 1}{2}z_B\f] */
  1984   ///
  1985   /// The algorithm can be executed with \c run() or the \c init() and
  1986   /// then the \c start() member functions. After it the matching can
  1987   /// be asked with \c matching() or mate() functions. The dual
  1988   /// solution can be get with \c nodeValue(), \c blossomNum() and \c
  1989   /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt
  1990   /// "BlossomIt" nested class which is able to iterate on the nodes
  1991   /// of a blossom. If the value type is integral then the dual
  1992   /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".
  1993   template <typename _Graph,
  1994             typename _WeightMap = typename _Graph::template EdgeMap<int> >
  1995   class MaxWeightedPerfectMatching {
  1996   public:
  1997 
  1998     typedef _Graph Graph;
  1999     typedef _WeightMap WeightMap;
  2000     typedef typename WeightMap::Value Value;
  2001 
  2002     /// \brief Scaling factor for dual solution
  2003     ///
  2004     /// Scaling factor for dual solution, it is equal to 4 or 1
  2005     /// according to the value type.
  2006     static const int dualScale =
  2007       std::numeric_limits<Value>::is_integer ? 4 : 1;
  2008 
  2009     typedef typename Graph::template NodeMap<typename Graph::Arc>
  2010     MatchingMap;
  2011 
  2012   private:
  2013 
  2014     TEMPLATE_GRAPH_TYPEDEFS(Graph);
  2015 
  2016     typedef typename Graph::template NodeMap<Value> NodePotential;
  2017     typedef std::vector<Node> BlossomNodeList;
  2018 
  2019     struct BlossomVariable {
  2020       int begin, end;
  2021       Value value;
  2022 
  2023       BlossomVariable(int _begin, int _end, Value _value)
  2024         : begin(_begin), end(_end), value(_value) {}
  2025 
  2026     };
  2027 
  2028     typedef std::vector<BlossomVariable> BlossomPotential;
  2029 
  2030     const Graph& _graph;
  2031     const WeightMap& _weight;
  2032 
  2033     MatchingMap* _matching;
  2034 
  2035     NodePotential* _node_potential;
  2036 
  2037     BlossomPotential _blossom_potential;
  2038     BlossomNodeList _blossom_node_list;
  2039 
  2040     int _node_num;
  2041     int _blossom_num;
  2042 
  2043     typedef RangeMap<int> IntIntMap;
  2044 
  2045     enum Status {
  2046       EVEN = -1, MATCHED = 0, ODD = 1
  2047     };
  2048 
  2049     typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
  2050     struct BlossomData {
  2051       int tree;
  2052       Status status;
  2053       Arc pred, next;
  2054       Value pot, offset;
  2055     };
  2056 
  2057     IntNodeMap *_blossom_index;
  2058     BlossomSet *_blossom_set;
  2059     RangeMap<BlossomData>* _blossom_data;
  2060 
  2061     IntNodeMap *_node_index;
  2062     IntArcMap *_node_heap_index;
  2063 
  2064     struct NodeData {
  2065 
  2066       NodeData(IntArcMap& node_heap_index)
  2067         : heap(node_heap_index) {}
  2068 
  2069       int blossom;
  2070       Value pot;
  2071       BinHeap<Value, IntArcMap> heap;
  2072       std::map<int, Arc> heap_index;
  2073 
  2074       int tree;
  2075     };
  2076 
  2077     RangeMap<NodeData>* _node_data;
  2078 
  2079     typedef ExtendFindEnum<IntIntMap> TreeSet;
  2080 
  2081     IntIntMap *_tree_set_index;
  2082     TreeSet *_tree_set;
  2083 
  2084     IntIntMap *_delta2_index;
  2085     BinHeap<Value, IntIntMap> *_delta2;
  2086 
  2087     IntEdgeMap *_delta3_index;
  2088     BinHeap<Value, IntEdgeMap> *_delta3;
  2089 
  2090     IntIntMap *_delta4_index;
  2091     BinHeap<Value, IntIntMap> *_delta4;
  2092 
  2093     Value _delta_sum;
  2094 
  2095     void createStructures() {
  2096       _node_num = countNodes(_graph);
  2097       _blossom_num = _node_num * 3 / 2;
  2098 
  2099       if (!_matching) {
  2100         _matching = new MatchingMap(_graph);
  2101       }
  2102       if (!_node_potential) {
  2103         _node_potential = new NodePotential(_graph);
  2104       }
  2105       if (!_blossom_set) {
  2106         _blossom_index = new IntNodeMap(_graph);
  2107         _blossom_set = new BlossomSet(*_blossom_index);
  2108         _blossom_data = new RangeMap<BlossomData>(_blossom_num);
  2109       }
  2110 
  2111       if (!_node_index) {
  2112         _node_index = new IntNodeMap(_graph);
  2113         _node_heap_index = new IntArcMap(_graph);
  2114         _node_data = new RangeMap<NodeData>(_node_num,
  2115                                             NodeData(*_node_heap_index));
  2116       }
  2117 
  2118       if (!_tree_set) {
  2119         _tree_set_index = new IntIntMap(_blossom_num);
  2120         _tree_set = new TreeSet(*_tree_set_index);
  2121       }
  2122       if (!_delta2) {
  2123         _delta2_index = new IntIntMap(_blossom_num);
  2124         _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
  2125       }
  2126       if (!_delta3) {
  2127         _delta3_index = new IntEdgeMap(_graph);
  2128         _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
  2129       }
  2130       if (!_delta4) {
  2131         _delta4_index = new IntIntMap(_blossom_num);
  2132         _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
  2133       }
  2134     }
  2135 
  2136     void destroyStructures() {
  2137       _node_num = countNodes(_graph);
  2138       _blossom_num = _node_num * 3 / 2;
  2139 
  2140       if (_matching) {
  2141         delete _matching;
  2142       }
  2143       if (_node_potential) {
  2144         delete _node_potential;
  2145       }
  2146       if (_blossom_set) {
  2147         delete _blossom_index;
  2148         delete _blossom_set;
  2149         delete _blossom_data;
  2150       }
  2151 
  2152       if (_node_index) {
  2153         delete _node_index;
  2154         delete _node_heap_index;
  2155         delete _node_data;
  2156       }
  2157 
  2158       if (_tree_set) {
  2159         delete _tree_set_index;
  2160         delete _tree_set;
  2161       }
  2162       if (_delta2) {
  2163         delete _delta2_index;
  2164         delete _delta2;
  2165       }
  2166       if (_delta3) {
  2167         delete _delta3_index;
  2168         delete _delta3;
  2169       }
  2170       if (_delta4) {
  2171         delete _delta4_index;
  2172         delete _delta4;
  2173       }
  2174     }
  2175 
  2176     void matchedToEven(int blossom, int tree) {
  2177       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
  2178         _delta2->erase(blossom);
  2179       }
  2180 
  2181       if (!_blossom_set->trivial(blossom)) {
  2182         (*_blossom_data)[blossom].pot -=
  2183           2 * (_delta_sum - (*_blossom_data)[blossom].offset);
  2184       }
  2185 
  2186       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  2187            n != INVALID; ++n) {
  2188 
  2189         _blossom_set->increase(n, std::numeric_limits<Value>::max());
  2190         int ni = (*_node_index)[n];
  2191 
  2192         (*_node_data)[ni].heap.clear();
  2193         (*_node_data)[ni].heap_index.clear();
  2194 
  2195         (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;
  2196 
  2197         for (InArcIt e(_graph, n); e != INVALID; ++e) {
  2198           Node v = _graph.source(e);
  2199           int vb = _blossom_set->find(v);
  2200           int vi = (*_node_index)[v];
  2201 
  2202           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  2203             dualScale * _weight[e];
  2204 
  2205           if ((*_blossom_data)[vb].status == EVEN) {
  2206             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
  2207               _delta3->push(e, rw / 2);
  2208             }
  2209           } else {
  2210             typename std::map<int, Arc>::iterator it =
  2211               (*_node_data)[vi].heap_index.find(tree);
  2212 
  2213             if (it != (*_node_data)[vi].heap_index.end()) {
  2214               if ((*_node_data)[vi].heap[it->second] > rw) {
  2215                 (*_node_data)[vi].heap.replace(it->second, e);
  2216                 (*_node_data)[vi].heap.decrease(e, rw);
  2217                 it->second = e;
  2218               }
  2219             } else {
  2220               (*_node_data)[vi].heap.push(e, rw);
  2221               (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
  2222             }
  2223 
  2224             if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
  2225               _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
  2226 
  2227               if ((*_blossom_data)[vb].status == MATCHED) {
  2228                 if (_delta2->state(vb) != _delta2->IN_HEAP) {
  2229                   _delta2->push(vb, _blossom_set->classPrio(vb) -
  2230                                (*_blossom_data)[vb].offset);
  2231                 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
  2232                            (*_blossom_data)[vb].offset){
  2233                   _delta2->decrease(vb, _blossom_set->classPrio(vb) -
  2234                                    (*_blossom_data)[vb].offset);
  2235                 }
  2236               }
  2237             }
  2238           }
  2239         }
  2240       }
  2241       (*_blossom_data)[blossom].offset = 0;
  2242     }
  2243 
  2244     void matchedToOdd(int blossom) {
  2245       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
  2246         _delta2->erase(blossom);
  2247       }
  2248       (*_blossom_data)[blossom].offset += _delta_sum;
  2249       if (!_blossom_set->trivial(blossom)) {
  2250         _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
  2251                      (*_blossom_data)[blossom].offset);
  2252       }
  2253     }
  2254 
  2255     void evenToMatched(int blossom, int tree) {
  2256       if (!_blossom_set->trivial(blossom)) {
  2257         (*_blossom_data)[blossom].pot += 2 * _delta_sum;
  2258       }
  2259 
  2260       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  2261            n != INVALID; ++n) {
  2262         int ni = (*_node_index)[n];
  2263         (*_node_data)[ni].pot -= _delta_sum;
  2264 
  2265         for (InArcIt e(_graph, n); e != INVALID; ++e) {
  2266           Node v = _graph.source(e);
  2267           int vb = _blossom_set->find(v);
  2268           int vi = (*_node_index)[v];
  2269 
  2270           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  2271             dualScale * _weight[e];
  2272 
  2273           if (vb == blossom) {
  2274             if (_delta3->state(e) == _delta3->IN_HEAP) {
  2275               _delta3->erase(e);
  2276             }
  2277           } else if ((*_blossom_data)[vb].status == EVEN) {
  2278 
  2279             if (_delta3->state(e) == _delta3->IN_HEAP) {
  2280               _delta3->erase(e);
  2281             }
  2282 
  2283             int vt = _tree_set->find(vb);
  2284 
  2285             if (vt != tree) {
  2286 
  2287               Arc r = _graph.oppositeArc(e);
  2288 
  2289               typename std::map<int, Arc>::iterator it =
  2290                 (*_node_data)[ni].heap_index.find(vt);
  2291 
  2292               if (it != (*_node_data)[ni].heap_index.end()) {
  2293                 if ((*_node_data)[ni].heap[it->second] > rw) {
  2294                   (*_node_data)[ni].heap.replace(it->second, r);
  2295                   (*_node_data)[ni].heap.decrease(r, rw);
  2296                   it->second = r;
  2297                 }
  2298               } else {
  2299                 (*_node_data)[ni].heap.push(r, rw);
  2300                 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
  2301               }
  2302 
  2303               if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
  2304                 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
  2305 
  2306                 if (_delta2->state(blossom) != _delta2->IN_HEAP) {
  2307                   _delta2->push(blossom, _blossom_set->classPrio(blossom) -
  2308                                (*_blossom_data)[blossom].offset);
  2309                 } else if ((*_delta2)[blossom] >
  2310                            _blossom_set->classPrio(blossom) -
  2311                            (*_blossom_data)[blossom].offset){
  2312                   _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
  2313                                    (*_blossom_data)[blossom].offset);
  2314                 }
  2315               }
  2316             }
  2317           } else {
  2318 
  2319             typename std::map<int, Arc>::iterator it =
  2320               (*_node_data)[vi].heap_index.find(tree);
  2321 
  2322             if (it != (*_node_data)[vi].heap_index.end()) {
  2323               (*_node_data)[vi].heap.erase(it->second);
  2324               (*_node_data)[vi].heap_index.erase(it);
  2325               if ((*_node_data)[vi].heap.empty()) {
  2326                 _blossom_set->increase(v, std::numeric_limits<Value>::max());
  2327               } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
  2328                 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
  2329               }
  2330 
  2331               if ((*_blossom_data)[vb].status == MATCHED) {
  2332                 if (_blossom_set->classPrio(vb) ==
  2333                     std::numeric_limits<Value>::max()) {
  2334                   _delta2->erase(vb);
  2335                 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
  2336                            (*_blossom_data)[vb].offset) {
  2337                   _delta2->increase(vb, _blossom_set->classPrio(vb) -
  2338                                    (*_blossom_data)[vb].offset);
  2339                 }
  2340               }
  2341             }
  2342           }
  2343         }
  2344       }
  2345     }
  2346 
  2347     void oddToMatched(int blossom) {
  2348       (*_blossom_data)[blossom].offset -= _delta_sum;
  2349 
  2350       if (_blossom_set->classPrio(blossom) !=
  2351           std::numeric_limits<Value>::max()) {
  2352         _delta2->push(blossom, _blossom_set->classPrio(blossom) -
  2353                        (*_blossom_data)[blossom].offset);
  2354       }
  2355 
  2356       if (!_blossom_set->trivial(blossom)) {
  2357         _delta4->erase(blossom);
  2358       }
  2359     }
  2360 
  2361     void oddToEven(int blossom, int tree) {
  2362       if (!_blossom_set->trivial(blossom)) {
  2363         _delta4->erase(blossom);
  2364         (*_blossom_data)[blossom].pot -=
  2365           2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
  2366       }
  2367 
  2368       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
  2369            n != INVALID; ++n) {
  2370         int ni = (*_node_index)[n];
  2371 
  2372         _blossom_set->increase(n, std::numeric_limits<Value>::max());
  2373 
  2374         (*_node_data)[ni].heap.clear();
  2375         (*_node_data)[ni].heap_index.clear();
  2376         (*_node_data)[ni].pot +=
  2377           2 * _delta_sum - (*_blossom_data)[blossom].offset;
  2378 
  2379         for (InArcIt e(_graph, n); e != INVALID; ++e) {
  2380           Node v = _graph.source(e);
  2381           int vb = _blossom_set->find(v);
  2382           int vi = (*_node_index)[v];
  2383 
  2384           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  2385             dualScale * _weight[e];
  2386 
  2387           if ((*_blossom_data)[vb].status == EVEN) {
  2388             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
  2389               _delta3->push(e, rw / 2);
  2390             }
  2391           } else {
  2392 
  2393             typename std::map<int, Arc>::iterator it =
  2394               (*_node_data)[vi].heap_index.find(tree);
  2395 
  2396             if (it != (*_node_data)[vi].heap_index.end()) {
  2397               if ((*_node_data)[vi].heap[it->second] > rw) {
  2398                 (*_node_data)[vi].heap.replace(it->second, e);
  2399                 (*_node_data)[vi].heap.decrease(e, rw);
  2400                 it->second = e;
  2401               }
  2402             } else {
  2403               (*_node_data)[vi].heap.push(e, rw);
  2404               (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
  2405             }
  2406 
  2407             if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
  2408               _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
  2409 
  2410               if ((*_blossom_data)[vb].status == MATCHED) {
  2411                 if (_delta2->state(vb) != _delta2->IN_HEAP) {
  2412                   _delta2->push(vb, _blossom_set->classPrio(vb) -
  2413                                (*_blossom_data)[vb].offset);
  2414                 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
  2415                            (*_blossom_data)[vb].offset) {
  2416                   _delta2->decrease(vb, _blossom_set->classPrio(vb) -
  2417                                    (*_blossom_data)[vb].offset);
  2418                 }
  2419               }
  2420             }
  2421           }
  2422         }
  2423       }
  2424       (*_blossom_data)[blossom].offset = 0;
  2425     }
  2426 
  2427     void alternatePath(int even, int tree) {
  2428       int odd;
  2429 
  2430       evenToMatched(even, tree);
  2431       (*_blossom_data)[even].status = MATCHED;
  2432 
  2433       while ((*_blossom_data)[even].pred != INVALID) {
  2434         odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));
  2435         (*_blossom_data)[odd].status = MATCHED;
  2436         oddToMatched(odd);
  2437         (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;
  2438 
  2439         even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));
  2440         (*_blossom_data)[even].status = MATCHED;
  2441         evenToMatched(even, tree);
  2442         (*_blossom_data)[even].next =
  2443           _graph.oppositeArc((*_blossom_data)[odd].pred);
  2444       }
  2445 
  2446     }
  2447 
  2448     void destroyTree(int tree) {
  2449       for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {
  2450         if ((*_blossom_data)[b].status == EVEN) {
  2451           (*_blossom_data)[b].status = MATCHED;
  2452           evenToMatched(b, tree);
  2453         } else if ((*_blossom_data)[b].status == ODD) {
  2454           (*_blossom_data)[b].status = MATCHED;
  2455           oddToMatched(b);
  2456         }
  2457       }
  2458       _tree_set->eraseClass(tree);
  2459     }
  2460 
  2461     void augmentOnEdge(const Edge& edge) {
  2462 
  2463       int left = _blossom_set->find(_graph.u(edge));
  2464       int right = _blossom_set->find(_graph.v(edge));
  2465 
  2466       int left_tree = _tree_set->find(left);
  2467       alternatePath(left, left_tree);
  2468       destroyTree(left_tree);
  2469 
  2470       int right_tree = _tree_set->find(right);
  2471       alternatePath(right, right_tree);
  2472       destroyTree(right_tree);
  2473 
  2474       (*_blossom_data)[left].next = _graph.direct(edge, true);
  2475       (*_blossom_data)[right].next = _graph.direct(edge, false);
  2476     }
  2477 
  2478     void extendOnArc(const Arc& arc) {
  2479       int base = _blossom_set->find(_graph.target(arc));
  2480       int tree = _tree_set->find(base);
  2481 
  2482       int odd = _blossom_set->find(_graph.source(arc));
  2483       _tree_set->insert(odd, tree);
  2484       (*_blossom_data)[odd].status = ODD;
  2485       matchedToOdd(odd);
  2486       (*_blossom_data)[odd].pred = arc;
  2487 
  2488       int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));
  2489       (*_blossom_data)[even].pred = (*_blossom_data)[even].next;
  2490       _tree_set->insert(even, tree);
  2491       (*_blossom_data)[even].status = EVEN;
  2492       matchedToEven(even, tree);
  2493     }
  2494 
  2495     void shrinkOnEdge(const Edge& edge, int tree) {
  2496       int nca = -1;
  2497       std::vector<int> left_path, right_path;
  2498 
  2499       {
  2500         std::set<int> left_set, right_set;
  2501         int left = _blossom_set->find(_graph.u(edge));
  2502         left_path.push_back(left);
  2503         left_set.insert(left);
  2504 
  2505         int right = _blossom_set->find(_graph.v(edge));
  2506         right_path.push_back(right);
  2507         right_set.insert(right);
  2508 
  2509         while (true) {
  2510 
  2511           if ((*_blossom_data)[left].pred == INVALID) break;
  2512 
  2513           left =
  2514             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
  2515           left_path.push_back(left);
  2516           left =
  2517             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));
  2518           left_path.push_back(left);
  2519 
  2520           left_set.insert(left);
  2521 
  2522           if (right_set.find(left) != right_set.end()) {
  2523             nca = left;
  2524             break;
  2525           }
  2526 
  2527           if ((*_blossom_data)[right].pred == INVALID) break;
  2528 
  2529           right =
  2530             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
  2531           right_path.push_back(right);
  2532           right =
  2533             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));
  2534           right_path.push_back(right);
  2535 
  2536           right_set.insert(right);
  2537 
  2538           if (left_set.find(right) != left_set.end()) {
  2539             nca = right;
  2540             break;
  2541           }
  2542 
  2543         }
  2544 
  2545         if (nca == -1) {
  2546           if ((*_blossom_data)[left].pred == INVALID) {
  2547             nca = right;
  2548             while (left_set.find(nca) == left_set.end()) {
  2549               nca =
  2550                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  2551               right_path.push_back(nca);
  2552               nca =
  2553                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  2554               right_path.push_back(nca);
  2555             }
  2556           } else {
  2557             nca = left;
  2558             while (right_set.find(nca) == right_set.end()) {
  2559               nca =
  2560                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  2561               left_path.push_back(nca);
  2562               nca =
  2563                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));
  2564               left_path.push_back(nca);
  2565             }
  2566           }
  2567         }
  2568       }
  2569 
  2570       std::vector<int> subblossoms;
  2571       Arc prev;
  2572 
  2573       prev = _graph.direct(edge, true);
  2574       for (int i = 0; left_path[i] != nca; i += 2) {
  2575         subblossoms.push_back(left_path[i]);
  2576         (*_blossom_data)[left_path[i]].next = prev;
  2577         _tree_set->erase(left_path[i]);
  2578 
  2579         subblossoms.push_back(left_path[i + 1]);
  2580         (*_blossom_data)[left_path[i + 1]].status = EVEN;
  2581         oddToEven(left_path[i + 1], tree);
  2582         _tree_set->erase(left_path[i + 1]);
  2583         prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);
  2584       }
  2585 
  2586       int k = 0;
  2587       while (right_path[k] != nca) ++k;
  2588 
  2589       subblossoms.push_back(nca);
  2590       (*_blossom_data)[nca].next = prev;
  2591 
  2592       for (int i = k - 2; i >= 0; i -= 2) {
  2593         subblossoms.push_back(right_path[i + 1]);
  2594         (*_blossom_data)[right_path[i + 1]].status = EVEN;
  2595         oddToEven(right_path[i + 1], tree);
  2596         _tree_set->erase(right_path[i + 1]);
  2597 
  2598         (*_blossom_data)[right_path[i + 1]].next =
  2599           (*_blossom_data)[right_path[i + 1]].pred;
  2600 
  2601         subblossoms.push_back(right_path[i]);
  2602         _tree_set->erase(right_path[i]);
  2603       }
  2604 
  2605       int surface =
  2606         _blossom_set->join(subblossoms.begin(), subblossoms.end());
  2607 
  2608       for (int i = 0; i < int(subblossoms.size()); ++i) {
  2609         if (!_blossom_set->trivial(subblossoms[i])) {
  2610           (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;
  2611         }
  2612         (*_blossom_data)[subblossoms[i]].status = MATCHED;
  2613       }
  2614 
  2615       (*_blossom_data)[surface].pot = -2 * _delta_sum;
  2616       (*_blossom_data)[surface].offset = 0;
  2617       (*_blossom_data)[surface].status = EVEN;
  2618       (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;
  2619       (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;
  2620 
  2621       _tree_set->insert(surface, tree);
  2622       _tree_set->erase(nca);
  2623     }
  2624 
  2625     void splitBlossom(int blossom) {
  2626       Arc next = (*_blossom_data)[blossom].next;
  2627       Arc pred = (*_blossom_data)[blossom].pred;
  2628 
  2629       int tree = _tree_set->find(blossom);
  2630 
  2631       (*_blossom_data)[blossom].status = MATCHED;
  2632       oddToMatched(blossom);
  2633       if (_delta2->state(blossom) == _delta2->IN_HEAP) {
  2634         _delta2->erase(blossom);
  2635       }
  2636 
  2637       std::vector<int> subblossoms;
  2638       _blossom_set->split(blossom, std::back_inserter(subblossoms));
  2639 
  2640       Value offset = (*_blossom_data)[blossom].offset;
  2641       int b = _blossom_set->find(_graph.source(pred));
  2642       int d = _blossom_set->find(_graph.source(next));
  2643 
  2644       int ib = -1, id = -1;
  2645       for (int i = 0; i < int(subblossoms.size()); ++i) {
  2646         if (subblossoms[i] == b) ib = i;
  2647         if (subblossoms[i] == d) id = i;
  2648 
  2649         (*_blossom_data)[subblossoms[i]].offset = offset;
  2650         if (!_blossom_set->trivial(subblossoms[i])) {
  2651           (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;
  2652         }
  2653         if (_blossom_set->classPrio(subblossoms[i]) !=
  2654             std::numeric_limits<Value>::max()) {
  2655           _delta2->push(subblossoms[i],
  2656                         _blossom_set->classPrio(subblossoms[i]) -
  2657                         (*_blossom_data)[subblossoms[i]].offset);
  2658         }
  2659       }
  2660 
  2661       if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {
  2662         for (int i = (id + 1) % subblossoms.size();
  2663              i != ib; i = (i + 2) % subblossoms.size()) {
  2664           int sb = subblossoms[i];
  2665           int tb = subblossoms[(i + 1) % subblossoms.size()];
  2666           (*_blossom_data)[sb].next =
  2667             _graph.oppositeArc((*_blossom_data)[tb].next);
  2668         }
  2669 
  2670         for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {
  2671           int sb = subblossoms[i];
  2672           int tb = subblossoms[(i + 1) % subblossoms.size()];
  2673           int ub = subblossoms[(i + 2) % subblossoms.size()];
  2674 
  2675           (*_blossom_data)[sb].status = ODD;
  2676           matchedToOdd(sb);
  2677           _tree_set->insert(sb, tree);
  2678           (*_blossom_data)[sb].pred = pred;
  2679           (*_blossom_data)[sb].next =
  2680                            _graph.oppositeArc((*_blossom_data)[tb].next);
  2681 
  2682           pred = (*_blossom_data)[ub].next;
  2683 
  2684           (*_blossom_data)[tb].status = EVEN;
  2685           matchedToEven(tb, tree);
  2686           _tree_set->insert(tb, tree);
  2687           (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;
  2688         }
  2689 
  2690         (*_blossom_data)[subblossoms[id]].status = ODD;
  2691         matchedToOdd(subblossoms[id]);
  2692         _tree_set->insert(subblossoms[id], tree);
  2693         (*_blossom_data)[subblossoms[id]].next = next;
  2694         (*_blossom_data)[subblossoms[id]].pred = pred;
  2695 
  2696       } else {
  2697 
  2698         for (int i = (ib + 1) % subblossoms.size();
  2699              i != id; i = (i + 2) % subblossoms.size()) {
  2700           int sb = subblossoms[i];
  2701           int tb = subblossoms[(i + 1) % subblossoms.size()];
  2702           (*_blossom_data)[sb].next =
  2703             _graph.oppositeArc((*_blossom_data)[tb].next);
  2704         }
  2705 
  2706         for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {
  2707           int sb = subblossoms[i];
  2708           int tb = subblossoms[(i + 1) % subblossoms.size()];
  2709           int ub = subblossoms[(i + 2) % subblossoms.size()];
  2710 
  2711           (*_blossom_data)[sb].status = ODD;
  2712           matchedToOdd(sb);
  2713           _tree_set->insert(sb, tree);
  2714           (*_blossom_data)[sb].next = next;
  2715           (*_blossom_data)[sb].pred =
  2716             _graph.oppositeArc((*_blossom_data)[tb].next);
  2717 
  2718           (*_blossom_data)[tb].status = EVEN;
  2719           matchedToEven(tb, tree);
  2720           _tree_set->insert(tb, tree);
  2721           (*_blossom_data)[tb].pred =
  2722             (*_blossom_data)[tb].next =
  2723             _graph.oppositeArc((*_blossom_data)[ub].next);
  2724           next = (*_blossom_data)[ub].next;
  2725         }
  2726 
  2727         (*_blossom_data)[subblossoms[ib]].status = ODD;
  2728         matchedToOdd(subblossoms[ib]);
  2729         _tree_set->insert(subblossoms[ib], tree);
  2730         (*_blossom_data)[subblossoms[ib]].next = next;
  2731         (*_blossom_data)[subblossoms[ib]].pred = pred;
  2732       }
  2733       _tree_set->erase(blossom);
  2734     }
  2735 
  2736     void extractBlossom(int blossom, const Node& base, const Arc& matching) {
  2737       if (_blossom_set->trivial(blossom)) {
  2738         int bi = (*_node_index)[base];
  2739         Value pot = (*_node_data)[bi].pot;
  2740 
  2741         _matching->set(base, matching);
  2742         _blossom_node_list.push_back(base);
  2743         _node_potential->set(base, pot);
  2744       } else {
  2745 
  2746         Value pot = (*_blossom_data)[blossom].pot;
  2747         int bn = _blossom_node_list.size();
  2748 
  2749         std::vector<int> subblossoms;
  2750         _blossom_set->split(blossom, std::back_inserter(subblossoms));
  2751         int b = _blossom_set->find(base);
  2752         int ib = -1;
  2753         for (int i = 0; i < int(subblossoms.size()); ++i) {
  2754           if (subblossoms[i] == b) { ib = i; break; }
  2755         }
  2756 
  2757         for (int i = 1; i < int(subblossoms.size()); i += 2) {
  2758           int sb = subblossoms[(ib + i) % subblossoms.size()];
  2759           int tb = subblossoms[(ib + i + 1) % subblossoms.size()];
  2760 
  2761           Arc m = (*_blossom_data)[tb].next;
  2762           extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));
  2763           extractBlossom(tb, _graph.source(m), m);
  2764         }
  2765         extractBlossom(subblossoms[ib], base, matching);
  2766 
  2767         int en = _blossom_node_list.size();
  2768 
  2769         _blossom_potential.push_back(BlossomVariable(bn, en, pot));
  2770       }
  2771     }
  2772 
  2773     void extractMatching() {
  2774       std::vector<int> blossoms;
  2775       for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {
  2776         blossoms.push_back(c);
  2777       }
  2778 
  2779       for (int i = 0; i < int(blossoms.size()); ++i) {
  2780 
  2781         Value offset = (*_blossom_data)[blossoms[i]].offset;
  2782         (*_blossom_data)[blossoms[i]].pot += 2 * offset;
  2783         for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);
  2784              n != INVALID; ++n) {
  2785           (*_node_data)[(*_node_index)[n]].pot -= offset;
  2786         }
  2787 
  2788         Arc matching = (*_blossom_data)[blossoms[i]].next;
  2789         Node base = _graph.source(matching);
  2790         extractBlossom(blossoms[i], base, matching);
  2791       }
  2792     }
  2793 
  2794   public:
  2795 
  2796     /// \brief Constructor
  2797     ///
  2798     /// Constructor.
  2799     MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)
  2800       : _graph(graph), _weight(weight), _matching(0),
  2801         _node_potential(0), _blossom_potential(), _blossom_node_list(),
  2802         _node_num(0), _blossom_num(0),
  2803 
  2804         _blossom_index(0), _blossom_set(0), _blossom_data(0),
  2805         _node_index(0), _node_heap_index(0), _node_data(0),
  2806         _tree_set_index(0), _tree_set(0),
  2807 
  2808         _delta2_index(0), _delta2(0),
  2809         _delta3_index(0), _delta3(0),
  2810         _delta4_index(0), _delta4(0),
  2811 
  2812         _delta_sum() {}
  2813 
  2814     ~MaxWeightedPerfectMatching() {
  2815       destroyStructures();
  2816     }
  2817 
  2818     /// \name Execution control
  2819     /// The simplest way to execute the algorithm is to use the
  2820     /// \c run() member function.
  2821 
  2822     ///@{
  2823 
  2824     /// \brief Initialize the algorithm
  2825     ///
  2826     /// Initialize the algorithm
  2827     void init() {
  2828       createStructures();
  2829 
  2830       for (ArcIt e(_graph); e != INVALID; ++e) {
  2831         _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
  2832       }
  2833       for (EdgeIt e(_graph); e != INVALID; ++e) {
  2834         _delta3_index->set(e, _delta3->PRE_HEAP);
  2835       }
  2836       for (int i = 0; i < _blossom_num; ++i) {
  2837         _delta2_index->set(i, _delta2->PRE_HEAP);
  2838         _delta4_index->set(i, _delta4->PRE_HEAP);
  2839       }
  2840 
  2841       int index = 0;
  2842       for (NodeIt n(_graph); n != INVALID; ++n) {
  2843         Value max = - std::numeric_limits<Value>::max();
  2844         for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  2845           if (_graph.target(e) == n) continue;
  2846           if ((dualScale * _weight[e]) / 2 > max) {
  2847             max = (dualScale * _weight[e]) / 2;
  2848           }
  2849         }
  2850         _node_index->set(n, index);
  2851         (*_node_data)[index].pot = max;
  2852         int blossom =
  2853           _blossom_set->insert(n, std::numeric_limits<Value>::max());
  2854 
  2855         _tree_set->insert(blossom);
  2856 
  2857         (*_blossom_data)[blossom].status = EVEN;
  2858         (*_blossom_data)[blossom].pred = INVALID;
  2859         (*_blossom_data)[blossom].next = INVALID;
  2860         (*_blossom_data)[blossom].pot = 0;
  2861         (*_blossom_data)[blossom].offset = 0;
  2862         ++index;
  2863       }
  2864       for (EdgeIt e(_graph); e != INVALID; ++e) {
  2865         int si = (*_node_index)[_graph.u(e)];
  2866         int ti = (*_node_index)[_graph.v(e)];
  2867         if (_graph.u(e) != _graph.v(e)) {
  2868           _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
  2869                             dualScale * _weight[e]) / 2);
  2870         }
  2871       }
  2872     }
  2873 
  2874     /// \brief Starts the algorithm
  2875     ///
  2876     /// Starts the algorithm
  2877     bool start() {
  2878       enum OpType {
  2879         D2, D3, D4
  2880       };
  2881 
  2882       int unmatched = _node_num;
  2883       while (unmatched > 0) {
  2884         Value d2 = !_delta2->empty() ?
  2885           _delta2->prio() : std::numeric_limits<Value>::max();
  2886 
  2887         Value d3 = !_delta3->empty() ?
  2888           _delta3->prio() : std::numeric_limits<Value>::max();
  2889 
  2890         Value d4 = !_delta4->empty() ?
  2891           _delta4->prio() : std::numeric_limits<Value>::max();
  2892 
  2893         _delta_sum = d2; OpType ot = D2;
  2894         if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
  2895         if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
  2896 
  2897         if (_delta_sum == std::numeric_limits<Value>::max()) {
  2898           return false;
  2899         }
  2900 
  2901         switch (ot) {
  2902         case D2:
  2903           {
  2904             int blossom = _delta2->top();
  2905             Node n = _blossom_set->classTop(blossom);
  2906             Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
  2907             extendOnArc(e);
  2908           }
  2909           break;
  2910         case D3:
  2911           {
  2912             Edge e = _delta3->top();
  2913 
  2914             int left_blossom = _blossom_set->find(_graph.u(e));
  2915             int right_blossom = _blossom_set->find(_graph.v(e));
  2916 
  2917             if (left_blossom == right_blossom) {
  2918               _delta3->pop();
  2919             } else {
  2920               int left_tree = _tree_set->find(left_blossom);
  2921               int right_tree = _tree_set->find(right_blossom);
  2922 
  2923               if (left_tree == right_tree) {
  2924                 shrinkOnEdge(e, left_tree);
  2925               } else {
  2926                 augmentOnEdge(e);
  2927                 unmatched -= 2;
  2928               }
  2929             }
  2930           } break;
  2931         case D4:
  2932           splitBlossom(_delta4->top());
  2933           break;
  2934         }
  2935       }
  2936       extractMatching();
  2937       return true;
  2938     }
  2939 
  2940     /// \brief Runs %MaxWeightedPerfectMatching algorithm.
  2941     ///
  2942     /// This method runs the %MaxWeightedPerfectMatching algorithm.
  2943     ///
  2944     /// \note mwm.run() is just a shortcut of the following code.
  2945     /// \code
  2946     ///   mwm.init();
  2947     ///   mwm.start();
  2948     /// \endcode
  2949     bool run() {
  2950       init();
  2951       return start();
  2952     }
  2953 
  2954     /// @}
  2955 
  2956     /// \name Primal solution
  2957     /// Functions to get the primal solution, ie. the matching.
  2958 
  2959     /// @{
  2960 
  2961     /// \brief Returns the matching value.
  2962     ///
  2963     /// Returns the matching value.
  2964     Value matchingValue() const {
  2965       Value sum = 0;
  2966       for (NodeIt n(_graph); n != INVALID; ++n) {
  2967         if ((*_matching)[n] != INVALID) {
  2968           sum += _weight[(*_matching)[n]];
  2969         }
  2970       }
  2971       return sum /= 2;
  2972     }
  2973 
  2974     /// \brief Returns true when the edge is in the matching.
  2975     ///
  2976     /// Returns true when the edge is in the matching.
  2977     bool matching(const Edge& edge) const {
  2978       return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
  2979     }
  2980 
  2981     /// \brief Returns the incident matching edge.
  2982     ///
  2983     /// Returns the incident matching arc from given edge.
  2984     Arc matching(const Node& node) const {
  2985       return (*_matching)[node];
  2986     }
  2987 
  2988     /// \brief Returns the mate of the node.
  2989     ///
  2990     /// Returns the adjancent node in a mathcing arc.
  2991     Node mate(const Node& node) const {
  2992       return _graph.target((*_matching)[node]);
  2993     }
  2994 
  2995     /// @}
  2996 
  2997     /// \name Dual solution
  2998     /// Functions to get the dual solution.
  2999 
  3000     /// @{
  3001 
  3002     /// \brief Returns the value of the dual solution.
  3003     ///
  3004     /// Returns the value of the dual solution. It should be equal to
  3005     /// the primal value scaled by \ref dualScale "dual scale".
  3006     Value dualValue() const {
  3007       Value sum = 0;
  3008       for (NodeIt n(_graph); n != INVALID; ++n) {
  3009         sum += nodeValue(n);
  3010       }
  3011       for (int i = 0; i < blossomNum(); ++i) {
  3012         sum += blossomValue(i) * (blossomSize(i) / 2);
  3013       }
  3014       return sum;
  3015     }
  3016 
  3017     /// \brief Returns the value of the node.
  3018     ///
  3019     /// Returns the the value of the node.
  3020     Value nodeValue(const Node& n) const {
  3021       return (*_node_potential)[n];
  3022     }
  3023 
  3024     /// \brief Returns the number of the blossoms in the basis.
  3025     ///
  3026     /// Returns the number of the blossoms in the basis.
  3027     /// \see BlossomIt
  3028     int blossomNum() const {
  3029       return _blossom_potential.size();
  3030     }
  3031 
  3032 
  3033     /// \brief Returns the number of the nodes in the blossom.
  3034     ///
  3035     /// Returns the number of the nodes in the blossom.
  3036     int blossomSize(int k) const {
  3037       return _blossom_potential[k].end - _blossom_potential[k].begin;
  3038     }
  3039 
  3040     /// \brief Returns the value of the blossom.
  3041     ///
  3042     /// Returns the the value of the blossom.
  3043     /// \see BlossomIt
  3044     Value blossomValue(int k) const {
  3045       return _blossom_potential[k].value;
  3046     }
  3047 
  3048     /// \brief Iterator for obtaining the nodes of the blossom.
  3049     ///
  3050     /// Iterator for obtaining the nodes of the blossom. This class
  3051     /// provides a common lemon style iterator for listing a
  3052     /// subset of the nodes.
  3053     class BlossomIt {
  3054     public:
  3055 
  3056       /// \brief Constructor.
  3057       ///
  3058       /// Constructor to get the nodes of the variable.
  3059       BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
  3060         : _algorithm(&algorithm)
  3061       {
  3062         _index = _algorithm->_blossom_potential[variable].begin;
  3063         _last = _algorithm->_blossom_potential[variable].end;
  3064       }
  3065 
  3066       /// \brief Conversion to node.
  3067       ///
  3068       /// Conversion to node.
  3069       operator Node() const {
  3070         return _algorithm->_blossom_node_list[_index];
  3071       }
  3072 
  3073       /// \brief Increment operator.
  3074       ///
  3075       /// Increment operator.
  3076       BlossomIt& operator++() {
  3077         ++_index;
  3078         return *this;
  3079       }
  3080 
  3081       /// \brief Validity checking
  3082       ///
  3083       /// Checks whether the iterator is invalid.
  3084       bool operator==(Invalid) const { return _index == _last; }
  3085 
  3086       /// \brief Validity checking
  3087       ///
  3088       /// Checks whether the iterator is valid.
  3089       bool operator!=(Invalid) const { return _index != _last; }
  3090 
  3091     private:
  3092       const MaxWeightedPerfectMatching* _algorithm;
  3093       int _last;
  3094       int _index;
  3095     };
  3096 
  3097     /// @}
  3098 
  3099   };
  3100 
  3101 
  3102 } //END OF NAMESPACE LEMON
  3103 
  3104 #endif //LEMON_MAX_MATCHING_H