lemon/max_matching.h
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 17 Apr 2009 18:04:36 +0200
changeset 609 e6927fe719e6
parent 425 cace3206223b
child 559 c5fd2d996909
permissions -rw-r--r--
Support >= and <= constraints in NetworkSimplex (#219, #234)

By default the same inequality constraints are supported as by
Circulation (the GEQ form), but the LEQ form can also be selected
using the problemType() function.

The documentation of the min. cost flow module is reworked and
extended with important notes and explanations about the different
variants of the problem and about the dual solution and optimality
conditions.
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     5  * Copyright (C) 2003-2009
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 
    19 #ifndef LEMON_MAX_MATCHING_H
    20 #define LEMON_MAX_MATCHING_H
    21 
    22 #include <vector>
    23 #include <queue>
    24 #include <set>
    25 #include <limits>
    26 
    27 #include <lemon/core.h>
    28 #include <lemon/unionfind.h>
    29 #include <lemon/bin_heap.h>
    30 #include <lemon/maps.h>
    31 
    32 ///\ingroup matching
    33 ///\file
    34 ///\brief Maximum matching algorithms in general graphs.
    35 
    36 namespace lemon {
    37 
    38   /// \ingroup matching
    39   ///
    40   /// \brief Edmonds' alternating forest maximum matching algorithm.
    41   ///
    42   /// This class implements Edmonds' alternating forest matching
    43   /// algorithm. The algorithm can be started from an arbitrary initial
    44   /// matching (the default is the empty one)
    45   ///
    46   /// The dual solution of the problem is a map of the nodes to
    47   /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c
    48   /// MATCHED/C showing the Gallai-Edmonds decomposition of the
    49   /// graph. The nodes in \c EVEN/D induce a graph with
    50   /// factor-critical components, the nodes in \c ODD/A form the
    51   /// barrier, and the nodes in \c MATCHED/C induce a graph having a
    52   /// perfect matching. The number of the factor-critical components
    53   /// minus the number of barrier nodes is a lower bound on the
    54   /// unmatched nodes, and the matching is optimal if and only if this bound is
    55   /// tight. This decomposition can be attained by calling \c
    56   /// decomposition() after running the algorithm.
    57   ///
    58   /// \param _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     /// startSparse() or startDense() functions.
   420 
   421     ///@{
   422 
   423     /// \brief Sets the actual matching to the empty matching.
   424     ///
   425     /// Sets the actual matching to the empty matching.
   426     ///
   427     void init() {
   428       createStructures();
   429       for(NodeIt n(_graph); n != INVALID; ++n) {
   430         _matching->set(n, INVALID);
   431         _status->set(n, UNMATCHED);
   432       }
   433     }
   434 
   435     ///\brief Finds an initial matching in a greedy way
   436     ///
   437     ///It finds an initial matching in a greedy way.
   438     void greedyInit() {
   439       createStructures();
   440       for (NodeIt n(_graph); n != INVALID; ++n) {
   441         _matching->set(n, INVALID);
   442         _status->set(n, UNMATCHED);
   443       }
   444       for (NodeIt n(_graph); n != INVALID; ++n) {
   445         if ((*_matching)[n] == INVALID) {
   446           for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
   447             Node v = _graph.target(a);
   448             if ((*_matching)[v] == INVALID && v != n) {
   449               _matching->set(n, a);
   450               _status->set(n, MATCHED);
   451               _matching->set(v, _graph.oppositeArc(a));
   452               _status->set(v, MATCHED);
   453               break;
   454             }
   455           }
   456         }
   457       }
   458     }
   459 
   460 
   461     /// \brief Initialize the matching from a map containing.
   462     ///
   463     /// Initialize the matching from a \c bool valued \c Edge map. This
   464     /// map must have the property that there are no two incident edges
   465     /// with true value, ie. it contains a matching.
   466     /// \return %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