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