lemon: lemon/max_matching.h@1c5d6e47921f

     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  */

    19 #ifndef LEMON_MAX_MATCHING_H

    20 #define LEMON_MAX_MATCHING_H

    22 #include <vector>

    23 #include <queue>

    24 #include <set>

    25 #include <limits>

    27 #include <lemon/core.h>

    28 #include <lemon/unionfind.h>

    29 #include <lemon/bin_heap.h>

    30 #include <lemon/maps.h>

    32 ///\ingroup matching

    33 ///\file

    34 ///\brief Maximum matching algorithms in general graphs.

    36 namespace lemon {

    38   /// \ingroup matching

    39   ///

    40   /// \brief Edmonds' alternating forest maximum matching algorithm.

    41   ///

    42   /// This class implements Edmonds' alternating forest matching

    43   /// algorithm. The algorithm can be started from an arbitrary initial

    44   /// matching (the default is the empty one)

    45   ///

    46   /// The dual solution of the problem is a map of the nodes to

    47   /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c

    48   /// MATCHED/C showing the Gallai-Edmonds decomposition of the

    49   /// graph. The nodes in \c EVEN/D induce a graph with

    50   /// factor-critical components, the nodes in \c ODD/A form the

    51   /// barrier, and the nodes in \c MATCHED/C induce a graph having a

    52   /// perfect matching. The number of the factor-critical components

    53   /// minus the number of barrier nodes is a lower bound on the

    54   /// unmatched nodes, and the matching is optimal if and only if this bound is

    55   /// tight. This decomposition can be attained by calling \c

    56   /// decomposition() after running the algorithm.

    57   ///

    58   /// \param _Graph The graph type the algorithm runs on.

    59   template <typename _Graph>

    60   class MaxMatching {

    61   public:

    63     typedef _Graph Graph;

    64     typedef typename Graph::template NodeMap<typename Graph::Arc>

    65     MatchingMap;

    67     ///\brief Indicates the Gallai-Edmonds decomposition of the graph.

    68     ///

    69     ///Indicates the Gallai-Edmonds decomposition of the graph. The

    70     ///nodes with Status \c EVEN/D induce a graph with factor-critical

    71     ///components, the nodes in \c ODD/A form the canonical barrier,

    72     ///and the nodes in \c MATCHED/C induce a graph having a perfect

    73     ///matching.

    74     enum Status {

    75       EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2

    76     };

    78     typedef typename Graph::template NodeMap<Status> StatusMap;

    80   private:

    82     TEMPLATE_GRAPH_TYPEDEFS(Graph);

    84     typedef UnionFindEnum<IntNodeMap> BlossomSet;

    85     typedef ExtendFindEnum<IntNodeMap> TreeSet;

    86     typedef RangeMap<Node> NodeIntMap;

    87     typedef MatchingMap EarMap;

    88     typedef std::vector<Node> NodeQueue;

    90     const Graph& _graph;

    91     MatchingMap* _matching;

    92     StatusMap* _status;

    94     EarMap* _ear;

    96     IntNodeMap* _blossom_set_index;

    97     BlossomSet* _blossom_set;

    98     NodeIntMap* _blossom_rep;

   100     IntNodeMap* _tree_set_index;

   101     TreeSet* _tree_set;

   103     NodeQueue _node_queue;

   104     int _process, _postpone, _last;

   106     int _node_num;

   108   private:

   110     void createStructures() {

   111       _node_num = countNodes(_graph);

   112       if (!_matching) {

   113         _matching = new MatchingMap(_graph);

   114       }

   115       if (!_status) {

   116         _status = new StatusMap(_graph);

   117       }

   118       if (!_ear) {

   119         _ear = new EarMap(_graph);

   120       }

   121       if (!_blossom_set) {

   122         _blossom_set_index = new IntNodeMap(_graph);

   123         _blossom_set = new BlossomSet(*_blossom_set_index);

   124       }

   125       if (!_blossom_rep) {

   126         _blossom_rep = new NodeIntMap(_node_num);

   127       }

   128       if (!_tree_set) {

   129         _tree_set_index = new IntNodeMap(_graph);

   130         _tree_set = new TreeSet(*_tree_set_index);

   131       }

   132       _node_queue.resize(_node_num);

   133     }

   135     void destroyStructures() {

   136       if (_matching) {

   137         delete _matching;

   138       }

   139       if (_status) {

   140         delete _status;

   141       }

   142       if (_ear) {

   143         delete _ear;

   144       }

   145       if (_blossom_set) {

   146         delete _blossom_set;

   147         delete _blossom_set_index;

   148       }

   149       if (_blossom_rep) {

   150         delete _blossom_rep;

   151       }

   152       if (_tree_set) {

   153         delete _tree_set_index;

   154         delete _tree_set;

   155       }

   156     }

   158     void processDense(const Node& n) {

   159       _process = _postpone = _last = 0;

   160       _node_queue[_last++] = n;

   162       while (_process != _last) {

   163         Node u = _node_queue[_process++];

   164         for (OutArcIt a(_graph, u); a != INVALID; ++a) {

   165           Node v = _graph.target(a);

   166           if ((*_status)[v] == MATCHED) {

   167             extendOnArc(a);

   168           } else if ((*_status)[v] == UNMATCHED) {

   169             augmentOnArc(a);

   170             return;

   171           }

   172         }

   173       }

   175       while (_postpone != _last) {

   176         Node u = _node_queue[_postpone++];

   178         for (OutArcIt a(_graph, u); a != INVALID ; ++a) {

   179           Node v = _graph.target(a);

   181           if ((*_status)[v] == EVEN) {

   182             if (_blossom_set->find(u) != _blossom_set->find(v)) {

   183               shrinkOnEdge(a);

   184             }

   185           }

   187           while (_process != _last) {

   188             Node w = _node_queue[_process++];

   189             for (OutArcIt b(_graph, w); b != INVALID; ++b) {

   190               Node x = _graph.target(b);

   191               if ((*_status)[x] == MATCHED) {

   192                 extendOnArc(b);

   193               } else if ((*_status)[x] == UNMATCHED) {

   194                 augmentOnArc(b);

   195                 return;

   196               }

   197             }

   198           }

   199         }

   200       }

   201     }

   203     void processSparse(const Node& n) {

   204       _process = _last = 0;

   205       _node_queue[_last++] = n;

   206       while (_process != _last) {

   207         Node u = _node_queue[_process++];

   208         for (OutArcIt a(_graph, u); a != INVALID; ++a) {

   209           Node v = _graph.target(a);

   211           if ((*_status)[v] == EVEN) {

   212             if (_blossom_set->find(u) != _blossom_set->find(v)) {

   213               shrinkOnEdge(a);

   214             }

   215           } else if ((*_status)[v] == MATCHED) {

   216             extendOnArc(a);

   217           } else if ((*_status)[v] == UNMATCHED) {

   218             augmentOnArc(a);

   219             return;

   220           }

   221         }

   222       }

   223     }

   225     void shrinkOnEdge(const Edge& e) {

   226       Node nca = INVALID;

   228       {

   229         std::set<Node> left_set, right_set;

   231         Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];

   232         left_set.insert(left);

   234         Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];

   235         right_set.insert(right);

   237         while (true) {

   238           if ((*_matching)[left] == INVALID) break;

   239           left = _graph.target((*_matching)[left]);

   240           left = (*_blossom_rep)[_blossom_set->

   241                                  find(_graph.target((*_ear)[left]))];

   242           if (right_set.find(left) != right_set.end()) {

   243             nca = left;

   244             break;

   245           }

   246           left_set.insert(left);

   248           if ((*_matching)[right] == INVALID) break;

   249           right = _graph.target((*_matching)[right]);

   250           right = (*_blossom_rep)[_blossom_set->

   251                                   find(_graph.target((*_ear)[right]))];

   252           if (left_set.find(right) != left_set.end()) {

   253             nca = right;

   254             break;

   255           }

   256           right_set.insert(right);

   257         }

   259         if (nca == INVALID) {

   260           if ((*_matching)[left] == INVALID) {

   261             nca = right;

   262             while (left_set.find(nca) == left_set.end()) {

   263               nca = _graph.target((*_matching)[nca]);

   264               nca =(*_blossom_rep)[_blossom_set->

   265                                    find(_graph.target((*_ear)[nca]))];

   266             }

   267           } else {

   268             nca = left;

   269             while (right_set.find(nca) == right_set.end()) {

   270               nca = _graph.target((*_matching)[nca]);

   271               nca = (*_blossom_rep)[_blossom_set->

   272                                    find(_graph.target((*_ear)[nca]))];

   273             }

   274           }

   275         }

   276       }

   278       {

   280         Node node = _graph.u(e);

   281         Arc arc = _graph.direct(e, true);

   282         Node base = (*_blossom_rep)[_blossom_set->find(node)];

   284         while (base != nca) {

   285           _ear->set(node, arc);

   287           Node n = node;

   288           while (n != base) {

   289             n = _graph.target((*_matching)[n]);

   290             Arc a = (*_ear)[n];

   291             n = _graph.target(a);

   292             _ear->set(n, _graph.oppositeArc(a));

   293           }

   294           node = _graph.target((*_matching)[base]);

   295           _tree_set->erase(base);

   296           _tree_set->erase(node);

   297           _blossom_set->insert(node, _blossom_set->find(base));

   298           _status->set(node, EVEN);

   299           _node_queue[_last++] = node;

   300           arc = _graph.oppositeArc((*_ear)[node]);

   301           node = _graph.target((*_ear)[node]);

   302           base = (*_blossom_rep)[_blossom_set->find(node)];

   303           _blossom_set->join(_graph.target(arc), base);

   304         }

   305       }

   307       _blossom_rep->set(_blossom_set->find(nca), nca);

   309       {

   311         Node node = _graph.v(e);

   312         Arc arc = _graph.direct(e, false);

   313         Node base = (*_blossom_rep)[_blossom_set->find(node)];

   315         while (base != nca) {

   316           _ear->set(node, arc);

   318           Node n = node;

   319           while (n != base) {

   320             n = _graph.target((*_matching)[n]);

   321             Arc a = (*_ear)[n];

   322             n = _graph.target(a);

   323             _ear->set(n, _graph.oppositeArc(a));

   324           }

   325           node = _graph.target((*_matching)[base]);

   326           _tree_set->erase(base);

   327           _tree_set->erase(node);

   328           _blossom_set->insert(node, _blossom_set->find(base));

   329           _status->set(node, EVEN);

   330           _node_queue[_last++] = node;

   331           arc = _graph.oppositeArc((*_ear)[node]);

   332           node = _graph.target((*_ear)[node]);

   333           base = (*_blossom_rep)[_blossom_set->find(node)];

   334           _blossom_set->join(_graph.target(arc), base);

   335         }

   336       }

   338       _blossom_rep->set(_blossom_set->find(nca), nca);

   339     }

   343     void extendOnArc(const Arc& a) {

   344       Node base = _graph.source(a);

   345       Node odd = _graph.target(a);

   347       _ear->set(odd, _graph.oppositeArc(a));

   348       Node even = _graph.target((*_matching)[odd]);

   349       _blossom_rep->set(_blossom_set->insert(even), even);

   350       _status->set(odd, ODD);

   351       _status->set(even, EVEN);

   352       int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);

   353       _tree_set->insert(odd, tree);

   354       _tree_set->insert(even, tree);

   355       _node_queue[_last++] = even;

   357     }

   359     void augmentOnArc(const Arc& a) {

   360       Node even = _graph.source(a);

   361       Node odd = _graph.target(a);

   363       int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);

   365       _matching->set(odd, _graph.oppositeArc(a));

   366       _status->set(odd, MATCHED);

   368       Arc arc = (*_matching)[even];

   369       _matching->set(even, a);

   371       while (arc != INVALID) {

   372         odd = _graph.target(arc);

   373         arc = (*_ear)[odd];

   374         even = _graph.target(arc);

   375         _matching->set(odd, arc);

   376         arc = (*_matching)[even];

   377         _matching->set(even, _graph.oppositeArc((*_matching)[odd]));

   378       }

   380       for (typename TreeSet::ItemIt it(*_tree_set, tree);

   381            it != INVALID; ++it) {

   382         if ((*_status)[it] == ODD) {

   383           _status->set(it, MATCHED);

   384         } else {

   385           int blossom = _blossom_set->find(it);

   386           for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);

   387                jt != INVALID; ++jt) {

   388             _status->set(jt, MATCHED);

   389           }

   390           _blossom_set->eraseClass(blossom);

   391         }

   392       }

   393       _tree_set->eraseClass(tree);

   395     }

   397   public:

   399     /// \brief Constructor

   400     ///

   401     /// Constructor.

   402     MaxMatching(const Graph& graph)

   403       : _graph(graph), _matching(0), _status(0), _ear(0),

   404         _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),

   405         _tree_set_index(0), _tree_set(0) {}

   407     ~MaxMatching() {

   408       destroyStructures();

   409     }

   411     /// \name Execution control

   412     /// The simplest way to execute the algorithm is to use the

   413     /// \c run() member function.

   414     /// \n

   416     /// If you need better control on the execution, you must call

   417     /// \ref init(), \ref greedyInit() or \ref matchingInit()

   418     /// functions first, then you can start the algorithm with the \ref

   419     /// startSparse() or startDense() functions.

   421     ///@{

   423     /// \brief Sets the actual matching to the empty matching.

   424     ///

   425     /// Sets the actual matching to the empty matching.

   426     ///

   427     void init() {

   428       createStructures();

   429       for(NodeIt n(_graph); n != INVALID; ++n) {

   430         _matching->set(n, INVALID);

   431         _status->set(n, UNMATCHED);

   432       }

   433     }

   435     ///\brief Finds an initial matching in a greedy way

   436     ///

   437     ///It finds an initial matching in a greedy way.

   438     void greedyInit() {

   439       createStructures();

   440       for (NodeIt n(_graph); n != INVALID; ++n) {

   441         _matching->set(n, INVALID);

   442         _status->set(n, UNMATCHED);

   443       }

   444       for (NodeIt n(_graph); n != INVALID; ++n) {

   445         if ((*_matching)[n] == INVALID) {

   446           for (OutArcIt a(_graph, n); a != INVALID ; ++a) {

   447             Node v = _graph.target(a);

   448             if ((*_matching)[v] == INVALID && v != n) {

   449               _matching->set(n, a);

   450               _status->set(n, MATCHED);

   451               _matching->set(v, _graph.oppositeArc(a));

   452               _status->set(v, MATCHED);

   453               break;

   454             }

   455           }

   456         }

   457       }

   458     }

   461     /// \brief Initialize the matching from a map containing.

   462     ///

   463     /// Initialize the matching from a \c bool valued \c Edge map. This

   464     /// map must have the property that there are no two incident edges

   465     /// with true value, ie. it contains a matching.

   466     /// \return %True if the map contains a matching.

   467     template <typename MatchingMap>

   468     bool matchingInit(const MatchingMap& matching) {

   469       createStructures();

   471       for (NodeIt n(_graph); n != INVALID; ++n) {

   472         _matching->set(n, INVALID);

   473         _status->set(n, UNMATCHED);

   474       }

   475       for(EdgeIt e(_graph); e!=INVALID; ++e) {

   476         if (matching[e]) {

   478           Node u = _graph.u(e);

   479           if ((*_matching)[u] != INVALID) return false;

   480           _matching->set(u, _graph.direct(e, true));

   481           _status->set(u, MATCHED);

   483           Node v = _graph.v(e);

   484           if ((*_matching)[v] != INVALID) return false;

   485           _matching->set(v, _graph.direct(e, false));

   486           _status->set(v, MATCHED);

   487         }

   488       }

   489       return true;

   490     }

   492     /// \brief Starts Edmonds' algorithm

   493     ///

   494     /// If runs the original Edmonds' algorithm.

   495     void startSparse() {

   496       for(NodeIt n(_graph); n != INVALID; ++n) {

   497         if ((*_status)[n] == UNMATCHED) {

   498           (*_blossom_rep)[_blossom_set->insert(n)] = n;

   499           _tree_set->insert(n);

   500           _status->set(n, EVEN);

   501           processSparse(n);

   502         }

   503       }

   504     }

   506     /// \brief Starts Edmonds' algorithm.

   507     ///

   508     /// It runs Edmonds' algorithm with a heuristic of postponing

   509     /// shrinks, therefore resulting in a faster algorithm for dense graphs.

   510     void startDense() {

   511       for(NodeIt n(_graph); n != INVALID; ++n) {

   512         if ((*_status)[n] == UNMATCHED) {

   513           (*_blossom_rep)[_blossom_set->insert(n)] = n;

   514           _tree_set->insert(n);

   515           _status->set(n, EVEN);

   516           processDense(n);

   517         }

   518       }

   519     }

   522     /// \brief Runs Edmonds' algorithm

   523     ///

   524     /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>)

   525     /// or Edmonds' algorithm with a heuristic of

   526     /// postponing shrinks for dense graphs.

   527     void run() {

   528       if (countEdges(_graph) < 2 * countNodes(_graph)) {

   529         greedyInit();

   530         startSparse();

   531       } else {

   532         init();

   533         startDense();

   534       }

   535     }

   537     /// @}

   539     /// \name Primal solution

   540     /// Functions to get the primal solution, ie. the matching.

   542     /// @{

   544     ///\brief Returns the size of the current matching.

   545     ///

   546     ///Returns the size of the current matching. After \ref

   547     ///run() it returns the size of the maximum matching in the graph.

   548     int matchingSize() const {

   549       int size = 0;

   550       for (NodeIt n(_graph); n != INVALID; ++n) {

   551         if ((*_matching)[n] != INVALID) {

   552           ++size;

   553         }

   554       }

   555       return size / 2;

   556     }

   558     /// \brief Returns true when the edge is in the matching.

   559     ///

   560     /// Returns true when the edge is in the matching.

   561     bool matching(const Edge& edge) const {

   562       return edge == (*_matching)[_graph.u(edge)];

   563     }

   565     /// \brief Returns the matching edge incident to the given node.

   566     ///

   567     /// Returns the matching edge of a \c node in the actual matching or

   568     /// INVALID if the \c node is not covered by the actual matching.

   569     Arc matching(const Node& n) const {

   570       return (*_matching)[n];

   571     }

   573     ///\brief Returns the mate of a node in the actual matching.

   574     ///

   575     ///Returns the mate of a \c node in the actual matching or

   576     ///INVALID if the \c node is not covered by the actual matching.

   577     Node mate(const Node& n) const {

   578       return (*_matching)[n] != INVALID ?

   579         _graph.target((*_matching)[n]) : INVALID;

   580     }

   582     /// @}

   584     /// \name Dual solution

   585     /// Functions to get the dual solution, ie. the decomposition.

   587     /// @{

   589     /// \brief Returns the class of the node in the Edmonds-Gallai

   590     /// decomposition.

   591     ///

   592     /// Returns the class of the node in the Edmonds-Gallai

   593     /// decomposition.

   594     Status decomposition(const Node& n) const {

   595       return (*_status)[n];

   596     }

   598     /// \brief Returns true when the node is in the barrier.

   599     ///

   600     /// Returns true when the node is in the barrier.

   601     bool barrier(const Node& n) const {

   602       return (*_status)[n] == ODD;

   603     }

   605     /// @}

   607   };

   609   /// \ingroup matching

   610   ///

   611   /// \brief Weighted matching in general graphs

   612   ///

   613   /// This class provides an efficient implementation of Edmond's

   614   /// maximum weighted matching algorithm. The implementation is based

   615   /// on extensive use of priority queues and provides

   616   /// \f$O(nm\log(n))\f$ time complexity.

   617   ///

   618   /// The maximum weighted matching problem is to find undirected

   619   /// edges in the graph with maximum overall weight and no two of

   620   /// them shares their ends. The problem can be formulated with the

   621   /// following linear program.

   622   /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]

   623   /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}

   624       \quad \forall B\in\mathcal{O}\f] */

   625   /// \f[x_e \ge 0\quad \forall e\in E\f]

   626   /// \f[\max \sum_{e\in E}x_ew_e\f]

   627   /// where \f$\delta(X)\f$ is the set of edges incident to a node in

   628   /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in

   629   /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality

   630   /// subsets of the nodes.

   631   ///

   632   /// The algorithm calculates an optimal matching and a proof of the

   633   /// optimality. The solution of the dual problem can be used to check

   634   /// the result of the algorithm. The dual linear problem is the

   635   /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}

   636       z_B \ge w_{uv} \quad \forall uv\in E\f] */

   637   /// \f[y_u \ge 0 \quad \forall u \in V\f]

   638   /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]

   639   /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}

   640       \frac{\vert B \vert - 1}{2}z_B\f] */

   641   ///

   642   /// The algorithm can be executed with \c run() or the \c init() and

   643   /// then the \c start() member functions. After it the matching can

   644   /// be asked with \c matching() or mate() functions. The dual

   645   /// solution can be get with \c nodeValue(), \c blossomNum() and \c

   646   /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt

   647   /// "BlossomIt" nested class, which is able to iterate on the nodes

   648   /// of a blossom. If the value type is integral then the dual

   649   /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".

   650   template <typename _Graph,

   651             typename _WeightMap = typename _Graph::template EdgeMap<int> >

   652   class MaxWeightedMatching {

   653   public:

   655     typedef _Graph Graph;

   656     typedef _WeightMap WeightMap;

   657     typedef typename WeightMap::Value Value;

   659     /// \brief Scaling factor for dual solution

   660     ///

   661     /// Scaling factor for dual solution, it is equal to 4 or 1

   662     /// according to the value type.

   663     static const int dualScale =

   664       std::numeric_limits<Value>::is_integer ? 4 : 1;

   666     typedef typename Graph::template NodeMap<typename Graph::Arc>

   667     MatchingMap;

   669   private:

   671     TEMPLATE_GRAPH_TYPEDEFS(Graph);

   673     typedef typename Graph::template NodeMap<Value> NodePotential;

   674     typedef std::vector<Node> BlossomNodeList;

   676     struct BlossomVariable {

   677       int begin, end;

   678       Value value;

   680       BlossomVariable(int _begin, int _end, Value _value)

   681         : begin(_begin), end(_end), value(_value) {}

   683     };

   685     typedef std::vector<BlossomVariable> BlossomPotential;

   687     const Graph& _graph;

   688     const WeightMap& _weight;

   690     MatchingMap* _matching;

   692     NodePotential* _node_potential;

   694     BlossomPotential _blossom_potential;

   695     BlossomNodeList _blossom_node_list;

   697     int _node_num;

   698     int _blossom_num;

   700     typedef RangeMap<int> IntIntMap;

   702     enum Status {

   703       EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2

   704     };

   706     typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;

   707     struct BlossomData {

   708       int tree;

   709       Status status;

   710       Arc pred, next;

   711       Value pot, offset;

   712       Node base;

   713     };

   715     IntNodeMap *_blossom_index;

   716     BlossomSet *_blossom_set;

   717     RangeMap<BlossomData>* _blossom_data;

   719     IntNodeMap *_node_index;

   720     IntArcMap *_node_heap_index;

   722     struct NodeData {

   724       NodeData(IntArcMap& node_heap_index)

   725         : heap(node_heap_index) {}

   727       int blossom;

   728       Value pot;

   729       BinHeap<Value, IntArcMap> heap;

   730       std::map<int, Arc> heap_index;

   732       int tree;

   733     };

   735     RangeMap<NodeData>* _node_data;

   737     typedef ExtendFindEnum<IntIntMap> TreeSet;

   739     IntIntMap *_tree_set_index;

   740     TreeSet *_tree_set;

   742     IntNodeMap *_delta1_index;

   743     BinHeap<Value, IntNodeMap> *_delta1;

   745     IntIntMap *_delta2_index;

   746     BinHeap<Value, IntIntMap> *_delta2;

   748     IntEdgeMap *_delta3_index;

   749     BinHeap<Value, IntEdgeMap> *_delta3;

   751     IntIntMap *_delta4_index;

   752     BinHeap<Value, IntIntMap> *_delta4;

   754     Value _delta_sum;

   756     void createStructures() {

   757       _node_num = countNodes(_graph);

   758       _blossom_num = _node_num * 3 / 2;

   760       if (!_matching) {

   761         _matching = new MatchingMap(_graph);

   762       }

   763       if (!_node_potential) {

   764         _node_potential = new NodePotential(_graph);

   765       }

   766       if (!_blossom_set) {

   767         _blossom_index = new IntNodeMap(_graph);

   768         _blossom_set = new BlossomSet(*_blossom_index);

   769         _blossom_data = new RangeMap<BlossomData>(_blossom_num);

   770       }

   772       if (!_node_index) {

   773         _node_index = new IntNodeMap(_graph);

   774         _node_heap_index = new IntArcMap(_graph);

   775         _node_data = new RangeMap<NodeData>(_node_num,

   776                                               NodeData(*_node_heap_index));

   777       }

   779       if (!_tree_set) {

   780         _tree_set_index = new IntIntMap(_blossom_num);

   781         _tree_set = new TreeSet(*_tree_set_index);

   782       }

   783       if (!_delta1) {

   784         _delta1_index = new IntNodeMap(_graph);

   785         _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);

   786       }

   787       if (!_delta2) {

   788         _delta2_index = new IntIntMap(_blossom_num);

   789         _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);

   790       }

   791       if (!_delta3) {

   792         _delta3_index = new IntEdgeMap(_graph);

   793         _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);

   794       }

   795       if (!_delta4) {

   796         _delta4_index = new IntIntMap(_blossom_num);

   797         _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);

   798       }

   799     }

   801     void destroyStructures() {

   802       _node_num = countNodes(_graph);

   803       _blossom_num = _node_num * 3 / 2;

   805       if (_matching) {

   806         delete _matching;

   807       }

   808       if (_node_potential) {

   809         delete _node_potential;

   810       }

   811       if (_blossom_set) {

   812         delete _blossom_index;

   813         delete _blossom_set;

   814         delete _blossom_data;

   815       }

   817       if (_node_index) {

   818         delete _node_index;

   819         delete _node_heap_index;

   820         delete _node_data;

   821       }

   823       if (_tree_set) {

   824         delete _tree_set_index;

   825         delete _tree_set;

   826       }

   827       if (_delta1) {

   828         delete _delta1_index;

   829         delete _delta1;

   830       }

   831       if (_delta2) {

   832         delete _delta2_index;

   833         delete _delta2;

   834       }

   835       if (_delta3) {

   836         delete _delta3_index;

   837         delete _delta3;

   838       }

   839       if (_delta4) {

   840         delete _delta4_index;

   841         delete _delta4;

   842       }

   843     }

   845     void matchedToEven(int blossom, int tree) {

   846       if (_delta2->state(blossom) == _delta2->IN_HEAP) {

   847         _delta2->erase(blossom);

   848       }

   850       if (!_blossom_set->trivial(blossom)) {

   851         (*_blossom_data)[blossom].pot -=

   852           2 * (_delta_sum - (*_blossom_data)[blossom].offset);

   853       }

   855       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);

   856            n != INVALID; ++n) {

   858         _blossom_set->increase(n, std::numeric_limits<Value>::max());

   859         int ni = (*_node_index)[n];

   861         (*_node_data)[ni].heap.clear();

   862         (*_node_data)[ni].heap_index.clear();

   864         (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;

   866         _delta1->push(n, (*_node_data)[ni].pot);

   868         for (InArcIt e(_graph, n); e != INVALID; ++e) {

   869           Node v = _graph.source(e);

   870           int vb = _blossom_set->find(v);

   871           int vi = (*_node_index)[v];

   873           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -

   874             dualScale * _weight[e];

   876           if ((*_blossom_data)[vb].status == EVEN) {

   877             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {

   878               _delta3->push(e, rw / 2);

   879             }

   880           } else if ((*_blossom_data)[vb].status == UNMATCHED) {

   881             if (_delta3->state(e) != _delta3->IN_HEAP) {

   882               _delta3->push(e, rw);

   883             }

   884           } else {

   885             typename std::map<int, Arc>::iterator it =

   886               (*_node_data)[vi].heap_index.find(tree);

   888             if (it != (*_node_data)[vi].heap_index.end()) {

   889               if ((*_node_data)[vi].heap[it->second] > rw) {

   890                 (*_node_data)[vi].heap.replace(it->second, e);

   891                 (*_node_data)[vi].heap.decrease(e, rw);

   892                 it->second = e;

   893               }

   894             } else {

   895               (*_node_data)[vi].heap.push(e, rw);

   896               (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));

   897             }

   899             if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {

   900               _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());

   902               if ((*_blossom_data)[vb].status == MATCHED) {

   903                 if (_delta2->state(vb) != _delta2->IN_HEAP) {

   904                   _delta2->push(vb, _blossom_set->classPrio(vb) -

   905                                (*_blossom_data)[vb].offset);

   906                 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -

   907                            (*_blossom_data)[vb].offset){

   908                   _delta2->decrease(vb, _blossom_set->classPrio(vb) -

   909                                    (*_blossom_data)[vb].offset);

   910                 }

   911               }

   912             }

   913           }

   914         }

   915       }

   916       (*_blossom_data)[blossom].offset = 0;

   917     }

   919     void matchedToOdd(int blossom) {

   920       if (_delta2->state(blossom) == _delta2->IN_HEAP) {

   921         _delta2->erase(blossom);

   922       }

   923       (*_blossom_data)[blossom].offset += _delta_sum;

   924       if (!_blossom_set->trivial(blossom)) {

   925         _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +

   926                      (*_blossom_data)[blossom].offset);

   927       }

   928     }

   930     void evenToMatched(int blossom, int tree) {

   931       if (!_blossom_set->trivial(blossom)) {

   932         (*_blossom_data)[blossom].pot += 2 * _delta_sum;

   933       }

   935       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);

   936            n != INVALID; ++n) {

   937         int ni = (*_node_index)[n];

   938         (*_node_data)[ni].pot -= _delta_sum;

   940         _delta1->erase(n);

   942         for (InArcIt e(_graph, n); e != INVALID; ++e) {

   943           Node v = _graph.source(e);

   944           int vb = _blossom_set->find(v);

   945           int vi = (*_node_index)[v];

   947           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -

   948             dualScale * _weight[e];

   950           if (vb == blossom) {

   951             if (_delta3->state(e) == _delta3->IN_HEAP) {

   952               _delta3->erase(e);

   953             }

   954           } else if ((*_blossom_data)[vb].status == EVEN) {

   956             if (_delta3->state(e) == _delta3->IN_HEAP) {

   957               _delta3->erase(e);

   958             }

   960             int vt = _tree_set->find(vb);

   962             if (vt != tree) {

   964               Arc r = _graph.oppositeArc(e);

   966               typename std::map<int, Arc>::iterator it =

   967                 (*_node_data)[ni].heap_index.find(vt);

   969               if (it != (*_node_data)[ni].heap_index.end()) {

   970                 if ((*_node_data)[ni].heap[it->second] > rw) {

   971                   (*_node_data)[ni].heap.replace(it->second, r);

   972                   (*_node_data)[ni].heap.decrease(r, rw);

   973                   it->second = r;

   974                 }

   975               } else {

   976                 (*_node_data)[ni].heap.push(r, rw);

   977                 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));

   978               }

   980               if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {

   981                 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());

   983                 if (_delta2->state(blossom) != _delta2->IN_HEAP) {

   984                   _delta2->push(blossom, _blossom_set->classPrio(blossom) -

   985                                (*_blossom_data)[blossom].offset);

   986                 } else if ((*_delta2)[blossom] >

   987                            _blossom_set->classPrio(blossom) -

   988                            (*_blossom_data)[blossom].offset){

   989                   _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -

   990                                    (*_blossom_data)[blossom].offset);

   991                 }

   992               }

   993             }

   995           } else if ((*_blossom_data)[vb].status == UNMATCHED) {

   996             if (_delta3->state(e) == _delta3->IN_HEAP) {

   997               _delta3->erase(e);

   998             }

   999           } else {

  1001             typename std::map<int, Arc>::iterator it =

  1002               (*_node_data)[vi].heap_index.find(tree);

  1004             if (it != (*_node_data)[vi].heap_index.end()) {

  1005               (*_node_data)[vi].heap.erase(it->second);

  1006               (*_node_data)[vi].heap_index.erase(it);

  1007               if ((*_node_data)[vi].heap.empty()) {

  1008                 _blossom_set->increase(v, std::numeric_limits<Value>::max());

  1009               } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {

  1010                 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());

  1011               }

  1013               if ((*_blossom_data)[vb].status == MATCHED) {

  1014                 if (_blossom_set->classPrio(vb) ==

  1015                     std::numeric_limits<Value>::max()) {

  1016                   _delta2->erase(vb);

  1017                 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -

  1018                            (*_blossom_data)[vb].offset) {

  1019                   _delta2->increase(vb, _blossom_set->classPrio(vb) -

  1020                                    (*_blossom_data)[vb].offset);

  1021                 }

  1022               }

  1023             }

  1024           }

  1025         }

  1026       }

  1027     }

  1029     void oddToMatched(int blossom) {

  1030       (*_blossom_data)[blossom].offset -= _delta_sum;

  1032       if (_blossom_set->classPrio(blossom) !=

  1033           std::numeric_limits<Value>::max()) {

  1034         _delta2->push(blossom, _blossom_set->classPrio(blossom) -

  1035                        (*_blossom_data)[blossom].offset);

  1036       }

  1038       if (!_blossom_set->trivial(blossom)) {

  1039         _delta4->erase(blossom);

  1040       }

  1041     }

  1043     void oddToEven(int blossom, int tree) {

  1044       if (!_blossom_set->trivial(blossom)) {

  1045         _delta4->erase(blossom);

  1046         (*_blossom_data)[blossom].pot -=

  1047           2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);

  1048       }

  1050       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);

  1051            n != INVALID; ++n) {

  1052         int ni = (*_node_index)[n];

  1054         _blossom_set->increase(n, std::numeric_limits<Value>::max());

  1056         (*_node_data)[ni].heap.clear();

  1057         (*_node_data)[ni].heap_index.clear();

  1058         (*_node_data)[ni].pot +=

  1059           2 * _delta_sum - (*_blossom_data)[blossom].offset;

  1061         _delta1->push(n, (*_node_data)[ni].pot);

  1063         for (InArcIt e(_graph, n); e != INVALID; ++e) {

  1064           Node v = _graph.source(e);

  1065           int vb = _blossom_set->find(v);

  1066           int vi = (*_node_index)[v];

  1068           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -

  1069             dualScale * _weight[e];

  1071           if ((*_blossom_data)[vb].status == EVEN) {

  1072             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {

  1073               _delta3->push(e, rw / 2);

  1074             }

  1075           } else if ((*_blossom_data)[vb].status == UNMATCHED) {

  1076             if (_delta3->state(e) != _delta3->IN_HEAP) {

  1077               _delta3->push(e, rw);

  1078             }

  1079           } else {

  1081             typename std::map<int, Arc>::iterator it =

  1082               (*_node_data)[vi].heap_index.find(tree);

  1084             if (it != (*_node_data)[vi].heap_index.end()) {

  1085               if ((*_node_data)[vi].heap[it->second] > rw) {

  1086                 (*_node_data)[vi].heap.replace(it->second, e);

  1087                 (*_node_data)[vi].heap.decrease(e, rw);

  1088                 it->second = e;

  1089               }

  1090             } else {

  1091               (*_node_data)[vi].heap.push(e, rw);

  1092               (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));

  1093             }

  1095             if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {

  1096               _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());

  1098               if ((*_blossom_data)[vb].status == MATCHED) {

  1099                 if (_delta2->state(vb) != _delta2->IN_HEAP) {

  1100                   _delta2->push(vb, _blossom_set->classPrio(vb) -

  1101                                (*_blossom_data)[vb].offset);

  1102                 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -

  1103                            (*_blossom_data)[vb].offset) {

  1104                   _delta2->decrease(vb, _blossom_set->classPrio(vb) -

  1105                                    (*_blossom_data)[vb].offset);

  1106                 }

  1107               }

  1108             }

  1109           }

  1110         }

  1111       }

  1112       (*_blossom_data)[blossom].offset = 0;

  1113     }

  1116     void matchedToUnmatched(int blossom) {

  1117       if (_delta2->state(blossom) == _delta2->IN_HEAP) {

  1118         _delta2->erase(blossom);

  1119       }

  1121       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);

  1122            n != INVALID; ++n) {

  1123         int ni = (*_node_index)[n];

  1125         _blossom_set->increase(n, std::numeric_limits<Value>::max());

  1127         (*_node_data)[ni].heap.clear();

  1128         (*_node_data)[ni].heap_index.clear();

  1130         for (OutArcIt e(_graph, n); e != INVALID; ++e) {

  1131           Node v = _graph.target(e);

  1132           int vb = _blossom_set->find(v);

  1133           int vi = (*_node_index)[v];

  1135           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -

  1136             dualScale * _weight[e];

  1138           if ((*_blossom_data)[vb].status == EVEN) {

  1139             if (_delta3->state(e) != _delta3->IN_HEAP) {

  1140               _delta3->push(e, rw);

  1141             }

  1142           }

  1143         }

  1144       }

  1145     }

  1147     void unmatchedToMatched(int blossom) {

  1148       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);

  1149            n != INVALID; ++n) {

  1150         int ni = (*_node_index)[n];

  1152         for (InArcIt e(_graph, n); e != INVALID; ++e) {

  1153           Node v = _graph.source(e);

  1154           int vb = _blossom_set->find(v);

  1155           int vi = (*_node_index)[v];

  1157           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -

  1158             dualScale * _weight[e];

  1160           if (vb == blossom) {

  1161             if (_delta3->state(e) == _delta3->IN_HEAP) {

  1162               _delta3->erase(e);

  1163             }

  1164           } else if ((*_blossom_data)[vb].status == EVEN) {

  1166             if (_delta3->state(e) == _delta3->IN_HEAP) {

  1167               _delta3->erase(e);

  1168             }

  1170             int vt = _tree_set->find(vb);

  1172             Arc r = _graph.oppositeArc(e);

  1174             typename std::map<int, Arc>::iterator it =

  1175               (*_node_data)[ni].heap_index.find(vt);

  1177             if (it != (*_node_data)[ni].heap_index.end()) {

  1178               if ((*_node_data)[ni].heap[it->second] > rw) {

  1179                 (*_node_data)[ni].heap.replace(it->second, r);

  1180                 (*_node_data)[ni].heap.decrease(r, rw);

  1181                 it->second = r;

  1182               }

  1183             } else {

  1184               (*_node_data)[ni].heap.push(r, rw);

  1185               (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));

  1186             }

  1188             if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {

  1189               _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());

  1191               if (_delta2->state(blossom) != _delta2->IN_HEAP) {

  1192                 _delta2->push(blossom, _blossom_set->classPrio(blossom) -

  1193                              (*_blossom_data)[blossom].offset);

  1194               } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-

  1195                          (*_blossom_data)[blossom].offset){

  1196                 _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -

  1197                                  (*_blossom_data)[blossom].offset);

  1198               }

  1199             }

  1201           } else if ((*_blossom_data)[vb].status == UNMATCHED) {

  1202             if (_delta3->state(e) == _delta3->IN_HEAP) {

  1203               _delta3->erase(e);

  1204             }

  1205           }

  1206         }

  1207       }

  1208     }

  1210     void alternatePath(int even, int tree) {

  1211       int odd;

  1213       evenToMatched(even, tree);

  1214       (*_blossom_data)[even].status = MATCHED;

  1216       while ((*_blossom_data)[even].pred != INVALID) {

  1217         odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));

  1218         (*_blossom_data)[odd].status = MATCHED;

  1219         oddToMatched(odd);

  1220         (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;

  1222         even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));

  1223         (*_blossom_data)[even].status = MATCHED;

  1224         evenToMatched(even, tree);

  1225         (*_blossom_data)[even].next =

  1226           _graph.oppositeArc((*_blossom_data)[odd].pred);

  1227       }

  1229     }

  1231     void destroyTree(int tree) {

  1232       for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {

  1233         if ((*_blossom_data)[b].status == EVEN) {

  1234           (*_blossom_data)[b].status = MATCHED;

  1235           evenToMatched(b, tree);

  1236         } else if ((*_blossom_data)[b].status == ODD) {

  1237           (*_blossom_data)[b].status = MATCHED;

  1238           oddToMatched(b);

  1239         }

  1240       }

  1241       _tree_set->eraseClass(tree);

  1242     }

  1245     void unmatchNode(const Node& node) {

  1246       int blossom = _blossom_set->find(node);

  1247       int tree = _tree_set->find(blossom);

  1249       alternatePath(blossom, tree);

  1250       destroyTree(tree);

  1252       (*_blossom_data)[blossom].status = UNMATCHED;

  1253       (*_blossom_data)[blossom].base = node;

  1254       matchedToUnmatched(blossom);

  1255     }

  1258     void augmentOnEdge(const Edge& edge) {

  1260       int left = _blossom_set->find(_graph.u(edge));

  1261       int right = _blossom_set->find(_graph.v(edge));

  1263       if ((*_blossom_data)[left].status == EVEN) {

  1264         int left_tree = _tree_set->find(left);

  1265         alternatePath(left, left_tree);

  1266         destroyTree(left_tree);

  1267       } else {

  1268         (*_blossom_data)[left].status = MATCHED;

  1269         unmatchedToMatched(left);

  1270       }

  1272       if ((*_blossom_data)[right].status == EVEN) {

  1273         int right_tree = _tree_set->find(right);

  1274         alternatePath(right, right_tree);

  1275         destroyTree(right_tree);

  1276       } else {

  1277         (*_blossom_data)[right].status = MATCHED;

  1278         unmatchedToMatched(right);

  1279       }

  1281       (*_blossom_data)[left].next = _graph.direct(edge, true);

  1282       (*_blossom_data)[right].next = _graph.direct(edge, false);

  1283     }

  1285     void extendOnArc(const Arc& arc) {

  1286       int base = _blossom_set->find(_graph.target(arc));

  1287       int tree = _tree_set->find(base);

  1289       int odd = _blossom_set->find(_graph.source(arc));

  1290       _tree_set->insert(odd, tree);

  1291       (*_blossom_data)[odd].status = ODD;

  1292       matchedToOdd(odd);

  1293       (*_blossom_data)[odd].pred = arc;

  1295       int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));

  1296       (*_blossom_data)[even].pred = (*_blossom_data)[even].next;

  1297       _tree_set->insert(even, tree);

  1298       (*_blossom_data)[even].status = EVEN;

  1299       matchedToEven(even, tree);

  1300     }

  1302     void shrinkOnEdge(const Edge& edge, int tree) {

  1303       int nca = -1;

  1304       std::vector<int> left_path, right_path;

  1306       {

  1307         std::set<int> left_set, right_set;

  1308         int left = _blossom_set->find(_graph.u(edge));

  1309         left_path.push_back(left);

  1310         left_set.insert(left);

  1312         int right = _blossom_set->find(_graph.v(edge));

  1313         right_path.push_back(right);

  1314         right_set.insert(right);

  1316         while (true) {

  1318           if ((*_blossom_data)[left].pred == INVALID) break;

  1320           left =

  1321             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));

  1322           left_path.push_back(left);

  1323           left =

  1324             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));

  1325           left_path.push_back(left);

  1327           left_set.insert(left);

  1329           if (right_set.find(left) != right_set.end()) {

  1330             nca = left;

  1331             break;

  1332           }

  1334           if ((*_blossom_data)[right].pred == INVALID) break;

  1336           right =

  1337             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));

  1338           right_path.push_back(right);

  1339           right =

  1340             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));

  1341           right_path.push_back(right);

  1343           right_set.insert(right);

  1345           if (left_set.find(right) != left_set.end()) {

  1346             nca = right;

  1347             break;

  1348           }

  1350         }

  1352         if (nca == -1) {

  1353           if ((*_blossom_data)[left].pred == INVALID) {

  1354             nca = right;

  1355             while (left_set.find(nca) == left_set.end()) {

  1356               nca =

  1357                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));

  1358               right_path.push_back(nca);

  1359               nca =

  1360                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));

  1361               right_path.push_back(nca);

  1362             }

  1363           } else {

  1364             nca = left;

  1365             while (right_set.find(nca) == right_set.end()) {

  1366               nca =

  1367                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));

  1368               left_path.push_back(nca);

  1369               nca =

  1370                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));

  1371               left_path.push_back(nca);

  1372             }

  1373           }

  1374         }

  1375       }

  1377       std::vector<int> subblossoms;

  1378       Arc prev;

  1380       prev = _graph.direct(edge, true);

  1381       for (int i = 0; left_path[i] != nca; i += 2) {

  1382         subblossoms.push_back(left_path[i]);

  1383         (*_blossom_data)[left_path[i]].next = prev;

  1384         _tree_set->erase(left_path[i]);

  1386         subblossoms.push_back(left_path[i + 1]);

  1387         (*_blossom_data)[left_path[i + 1]].status = EVEN;

  1388         oddToEven(left_path[i + 1], tree);

  1389         _tree_set->erase(left_path[i + 1]);

  1390         prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);

  1391       }

  1393       int k = 0;

  1394       while (right_path[k] != nca) ++k;

  1396       subblossoms.push_back(nca);

  1397       (*_blossom_data)[nca].next = prev;

  1399       for (int i = k - 2; i >= 0; i -= 2) {

  1400         subblossoms.push_back(right_path[i + 1]);

  1401         (*_blossom_data)[right_path[i + 1]].status = EVEN;

  1402         oddToEven(right_path[i + 1], tree);

  1403         _tree_set->erase(right_path[i + 1]);

  1405         (*_blossom_data)[right_path[i + 1]].next =

  1406           (*_blossom_data)[right_path[i + 1]].pred;

  1408         subblossoms.push_back(right_path[i]);

  1409         _tree_set->erase(right_path[i]);

  1410       }

  1412       int surface =

  1413         _blossom_set->join(subblossoms.begin(), subblossoms.end());

  1415       for (int i = 0; i < int(subblossoms.size()); ++i) {

  1416         if (!_blossom_set->trivial(subblossoms[i])) {

  1417           (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;

  1418         }

  1419         (*_blossom_data)[subblossoms[i]].status = MATCHED;

  1420       }

  1422       (*_blossom_data)[surface].pot = -2 * _delta_sum;

  1423       (*_blossom_data)[surface].offset = 0;

  1424       (*_blossom_data)[surface].status = EVEN;

  1425       (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;

  1426       (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;

  1428       _tree_set->insert(surface, tree);

  1429       _tree_set->erase(nca);

  1430     }

  1432     void splitBlossom(int blossom) {

  1433       Arc next = (*_blossom_data)[blossom].next;

  1434       Arc pred = (*_blossom_data)[blossom].pred;

  1436       int tree = _tree_set->find(blossom);

  1438       (*_blossom_data)[blossom].status = MATCHED;

  1439       oddToMatched(blossom);

  1440       if (_delta2->state(blossom) == _delta2->IN_HEAP) {

  1441         _delta2->erase(blossom);

  1442       }

  1444       std::vector<int> subblossoms;

  1445       _blossom_set->split(blossom, std::back_inserter(subblossoms));

  1447       Value offset = (*_blossom_data)[blossom].offset;

  1448       int b = _blossom_set->find(_graph.source(pred));

  1449       int d = _blossom_set->find(_graph.source(next));

  1451       int ib = -1, id = -1;

  1452       for (int i = 0; i < int(subblossoms.size()); ++i) {

  1453         if (subblossoms[i] == b) ib = i;

  1454         if (subblossoms[i] == d) id = i;

  1456         (*_blossom_data)[subblossoms[i]].offset = offset;

  1457         if (!_blossom_set->trivial(subblossoms[i])) {

  1458           (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;

  1459         }

  1460         if (_blossom_set->classPrio(subblossoms[i]) !=

  1461             std::numeric_limits<Value>::max()) {

  1462           _delta2->push(subblossoms[i],

  1463                         _blossom_set->classPrio(subblossoms[i]) -

  1464                         (*_blossom_data)[subblossoms[i]].offset);

  1465         }

  1466       }

  1468       if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {

  1469         for (int i = (id + 1) % subblossoms.size();

  1470              i != ib; i = (i + 2) % subblossoms.size()) {

  1471           int sb = subblossoms[i];

  1472           int tb = subblossoms[(i + 1) % subblossoms.size()];

  1473           (*_blossom_data)[sb].next =

  1474             _graph.oppositeArc((*_blossom_data)[tb].next);

  1475         }

  1477         for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {

  1478           int sb = subblossoms[i];

  1479           int tb = subblossoms[(i + 1) % subblossoms.size()];

  1480           int ub = subblossoms[(i + 2) % subblossoms.size()];

  1482           (*_blossom_data)[sb].status = ODD;

  1483           matchedToOdd(sb);

  1484           _tree_set->insert(sb, tree);

  1485           (*_blossom_data)[sb].pred = pred;

  1486           (*_blossom_data)[sb].next =

  1487                            _graph.oppositeArc((*_blossom_data)[tb].next);

  1489           pred = (*_blossom_data)[ub].next;

  1491           (*_blossom_data)[tb].status = EVEN;

  1492           matchedToEven(tb, tree);

  1493           _tree_set->insert(tb, tree);

  1494           (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;

  1495         }

  1497         (*_blossom_data)[subblossoms[id]].status = ODD;

  1498         matchedToOdd(subblossoms[id]);

  1499         _tree_set->insert(subblossoms[id], tree);

  1500         (*_blossom_data)[subblossoms[id]].next = next;

  1501         (*_blossom_data)[subblossoms[id]].pred = pred;

  1503       } else {

  1505         for (int i = (ib + 1) % subblossoms.size();

  1506              i != id; i = (i + 2) % subblossoms.size()) {

  1507           int sb = subblossoms[i];

  1508           int tb = subblossoms[(i + 1) % subblossoms.size()];

  1509           (*_blossom_data)[sb].next =

  1510             _graph.oppositeArc((*_blossom_data)[tb].next);

  1511         }

  1513         for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {

  1514           int sb = subblossoms[i];

  1515           int tb = subblossoms[(i + 1) % subblossoms.size()];

  1516           int ub = subblossoms[(i + 2) % subblossoms.size()];

  1518           (*_blossom_data)[sb].status = ODD;

  1519           matchedToOdd(sb);

  1520           _tree_set->insert(sb, tree);

  1521           (*_blossom_data)[sb].next = next;

  1522           (*_blossom_data)[sb].pred =

  1523             _graph.oppositeArc((*_blossom_data)[tb].next);

  1525           (*_blossom_data)[tb].status = EVEN;

  1526           matchedToEven(tb, tree);

  1527           _tree_set->insert(tb, tree);

  1528           (*_blossom_data)[tb].pred =

  1529             (*_blossom_data)[tb].next =

  1530             _graph.oppositeArc((*_blossom_data)[ub].next);

  1531           next = (*_blossom_data)[ub].next;

  1532         }

  1534         (*_blossom_data)[subblossoms[ib]].status = ODD;

  1535         matchedToOdd(subblossoms[ib]);

  1536         _tree_set->insert(subblossoms[ib], tree);

  1537         (*_blossom_data)[subblossoms[ib]].next = next;

  1538         (*_blossom_data)[subblossoms[ib]].pred = pred;

  1539       }

  1540       _tree_set->erase(blossom);

  1541     }

  1543     void extractBlossom(int blossom, const Node& base, const Arc& matching) {

  1544       if (_blossom_set->trivial(blossom)) {

  1545         int bi = (*_node_index)[base];

  1546         Value pot = (*_node_data)[bi].pot;

  1548         _matching->set(base, matching);

  1549         _blossom_node_list.push_back(base);

  1550         _node_potential->set(base, pot);

  1551       } else {

  1553         Value pot = (*_blossom_data)[blossom].pot;

  1554         int bn = _blossom_node_list.size();

  1556         std::vector<int> subblossoms;

  1557         _blossom_set->split(blossom, std::back_inserter(subblossoms));

  1558         int b = _blossom_set->find(base);

  1559         int ib = -1;

  1560         for (int i = 0; i < int(subblossoms.size()); ++i) {

  1561           if (subblossoms[i] == b) { ib = i; break; }

  1562         }

  1564         for (int i = 1; i < int(subblossoms.size()); i += 2) {

  1565           int sb = subblossoms[(ib + i) % subblossoms.size()];

  1566           int tb = subblossoms[(ib + i + 1) % subblossoms.size()];

  1568           Arc m = (*_blossom_data)[tb].next;

  1569           extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));

  1570           extractBlossom(tb, _graph.source(m), m);

  1571         }

  1572         extractBlossom(subblossoms[ib], base, matching);

  1574         int en = _blossom_node_list.size();

  1576         _blossom_potential.push_back(BlossomVariable(bn, en, pot));

  1577       }

  1578     }

  1580     void extractMatching() {

  1581       std::vector<int> blossoms;

  1582       for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {

  1583         blossoms.push_back(c);

  1584       }

  1586       for (int i = 0; i < int(blossoms.size()); ++i) {

  1587         if ((*_blossom_data)[blossoms[i]].status == MATCHED) {

  1589           Value offset = (*_blossom_data)[blossoms[i]].offset;

  1590           (*_blossom_data)[blossoms[i]].pot += 2 * offset;

  1591           for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);

  1592                n != INVALID; ++n) {

  1593             (*_node_data)[(*_node_index)[n]].pot -= offset;

  1594           }

  1596           Arc matching = (*_blossom_data)[blossoms[i]].next;

  1597           Node base = _graph.source(matching);

  1598           extractBlossom(blossoms[i], base, matching);

  1599         } else {

  1600           Node base = (*_blossom_data)[blossoms[i]].base;

  1601           extractBlossom(blossoms[i], base, INVALID);

  1602         }

  1603       }

  1604     }

  1606   public:

  1608     /// \brief Constructor

  1609     ///

  1610     /// Constructor.

  1611     MaxWeightedMatching(const Graph& graph, const WeightMap& weight)

  1612       : _graph(graph), _weight(weight), _matching(0),

  1613         _node_potential(0), _blossom_potential(), _blossom_node_list(),

  1614         _node_num(0), _blossom_num(0),

  1616         _blossom_index(0), _blossom_set(0), _blossom_data(0),

  1617         _node_index(0), _node_heap_index(0), _node_data(0),

  1618         _tree_set_index(0), _tree_set(0),

  1620         _delta1_index(0), _delta1(0),

  1621         _delta2_index(0), _delta2(0),

  1622         _delta3_index(0), _delta3(0),

  1623         _delta4_index(0), _delta4(0),

  1625         _delta_sum() {}

  1627     ~MaxWeightedMatching() {

  1628       destroyStructures();

  1629     }

  1631     /// \name Execution control

  1632     /// The simplest way to execute the algorithm is to use the

  1633     /// \c run() member function.

  1635     ///@{

  1637     /// \brief Initialize the algorithm

  1638     ///

  1639     /// Initialize the algorithm

  1640     void init() {

  1641       createStructures();

  1643       for (ArcIt e(_graph); e != INVALID; ++e) {

  1644         _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);

  1645       }

  1646       for (NodeIt n(_graph); n != INVALID; ++n) {

  1647         _delta1_index->set(n, _delta1->PRE_HEAP);

  1648       }

  1649       for (EdgeIt e(_graph); e != INVALID; ++e) {

  1650         _delta3_index->set(e, _delta3->PRE_HEAP);

  1651       }

  1652       for (int i = 0; i < _blossom_num; ++i) {

  1653         _delta2_index->set(i, _delta2->PRE_HEAP);

  1654         _delta4_index->set(i, _delta4->PRE_HEAP);

  1655       }

  1657       int index = 0;

  1658       for (NodeIt n(_graph); n != INVALID; ++n) {

  1659         Value max = 0;

  1660         for (OutArcIt e(_graph, n); e != INVALID; ++e) {

  1661           if (_graph.target(e) == n) continue;

  1662           if ((dualScale * _weight[e]) / 2 > max) {

  1663             max = (dualScale * _weight[e]) / 2;

  1664           }

  1665         }

  1666         _node_index->set(n, index);

  1667         (*_node_data)[index].pot = max;

  1668         _delta1->push(n, max);

  1669         int blossom =

  1670           _blossom_set->insert(n, std::numeric_limits<Value>::max());

  1672         _tree_set->insert(blossom);

  1674         (*_blossom_data)[blossom].status = EVEN;

  1675         (*_blossom_data)[blossom].pred = INVALID;

  1676         (*_blossom_data)[blossom].next = INVALID;

  1677         (*_blossom_data)[blossom].pot = 0;

  1678         (*_blossom_data)[blossom].offset = 0;

  1679         ++index;

  1680       }

  1681       for (EdgeIt e(_graph); e != INVALID; ++e) {

  1682         int si = (*_node_index)[_graph.u(e)];

  1683         int ti = (*_node_index)[_graph.v(e)];

  1684         if (_graph.u(e) != _graph.v(e)) {

  1685           _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -

  1686                             dualScale * _weight[e]) / 2);

  1687         }

  1688       }

  1689     }

  1691     /// \brief Starts the algorithm

  1692     ///

  1693     /// Starts the algorithm

  1694     void start() {

  1695       enum OpType {

  1696         D1, D2, D3, D4

  1697       };

  1699       int unmatched = _node_num;

  1700       while (unmatched > 0) {

  1701         Value d1 = !_delta1->empty() ?

  1702           _delta1->prio() : std::numeric_limits<Value>::max();

  1704         Value d2 = !_delta2->empty() ?

  1705           _delta2->prio() : std::numeric_limits<Value>::max();

  1707         Value d3 = !_delta3->empty() ?

  1708           _delta3->prio() : std::numeric_limits<Value>::max();

  1710         Value d4 = !_delta4->empty() ?

  1711           _delta4->prio() : std::numeric_limits<Value>::max();

  1713         _delta_sum = d1; OpType ot = D1;

  1714         if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }

  1715         if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }

  1716         if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }

  1719         switch (ot) {

  1720         case D1:

  1721           {

  1722             Node n = _delta1->top();

  1723             unmatchNode(n);

  1724             --unmatched;

  1725           }

  1726           break;

  1727         case D2:

  1728           {

  1729             int blossom = _delta2->top();

  1730             Node n = _blossom_set->classTop(blossom);

  1731             Arc e = (*_node_data)[(*_node_index)[n]].heap.top();

  1732             extendOnArc(e);

  1733           }

  1734           break;

  1735         case D3:

  1736           {

  1737             Edge e = _delta3->top();

  1739             int left_blossom = _blossom_set->find(_graph.u(e));

  1740             int right_blossom = _blossom_set->find(_graph.v(e));

  1742             if (left_blossom == right_blossom) {

  1743               _delta3->pop();

  1744             } else {

  1745               int left_tree;

  1746               if ((*_blossom_data)[left_blossom].status == EVEN) {

  1747                 left_tree = _tree_set->find(left_blossom);

  1748               } else {

  1749                 left_tree = -1;

  1750                 ++unmatched;

  1751               }

  1752               int right_tree;

  1753               if ((*_blossom_data)[right_blossom].status == EVEN) {

  1754                 right_tree = _tree_set->find(right_blossom);

  1755               } else {

  1756                 right_tree = -1;

  1757                 ++unmatched;

  1758               }

  1760               if (left_tree == right_tree) {

  1761                 shrinkOnEdge(e, left_tree);

  1762               } else {

  1763                 augmentOnEdge(e);

  1764                 unmatched -= 2;

  1765               }

  1766             }

  1767           } break;

  1768         case D4:

  1769           splitBlossom(_delta4->top());

  1770           break;

  1771         }

  1772       }

  1773       extractMatching();

  1774     }

  1776     /// \brief Runs %MaxWeightedMatching algorithm.

  1777     ///

  1778     /// This method runs the %MaxWeightedMatching algorithm.

  1779     ///

  1780     /// \note mwm.run() is just a shortcut of the following code.

  1781     /// \code

  1782     ///   mwm.init();

  1783     ///   mwm.start();

  1784     /// \endcode

  1785     void run() {

  1786       init();

  1787       start();

  1788     }

  1790     /// @}

  1792     /// \name Primal solution

  1793     /// Functions to get the primal solution, ie. the matching.

  1795     /// @{

  1797     /// \brief Returns the weight of the matching.

  1798     ///

  1799     /// Returns the weight of the matching.

  1800     Value matchingValue() const {

  1801       Value sum = 0;

  1802       for (NodeIt n(_graph); n != INVALID; ++n) {

  1803         if ((*_matching)[n] != INVALID) {

  1804           sum += _weight[(*_matching)[n]];

  1805         }

  1806       }

  1807       return sum /= 2;

  1808     }

  1810     /// \brief Returns the cardinality of the matching.

  1811     ///

  1812     /// Returns the cardinality of the matching.

  1813     int matchingSize() const {

  1814       int num = 0;

  1815       for (NodeIt n(_graph); n != INVALID; ++n) {

  1816         if ((*_matching)[n] != INVALID) {

  1817           ++num;

  1818         }

  1819       }

  1820       return num /= 2;

  1821     }

  1823     /// \brief Returns true when the edge is in the matching.

  1824     ///

  1825     /// Returns true when the edge is in the matching.

  1826     bool matching(const Edge& edge) const {

  1827       return edge == (*_matching)[_graph.u(edge)];

  1828     }

  1830     /// \brief Returns the incident matching arc.

  1831     ///

  1832     /// Returns the incident matching arc from given node. If the

  1833     /// node is not matched then it gives back \c INVALID.

  1834     Arc matching(const Node& node) const {

  1835       return (*_matching)[node];

  1836     }

  1838     /// \brief Returns the mate of the node.

  1839     ///

  1840     /// Returns the adjancent node in a mathcing arc. If the node is

  1841     /// not matched then it gives back \c INVALID.

  1842     Node mate(const Node& node) const {

  1843       return (*_matching)[node] != INVALID ?

  1844         _graph.target((*_matching)[node]) : INVALID;

  1845     }

  1847     /// @}

  1849     /// \name Dual solution

  1850     /// Functions to get the dual solution.

  1852     /// @{

  1854     /// \brief Returns the value of the dual solution.

  1855     ///

  1856     /// Returns the value of the dual solution. It should be equal to

  1857     /// the primal value scaled by \ref dualScale "dual scale".

  1858     Value dualValue() const {

  1859       Value sum = 0;

  1860       for (NodeIt n(_graph); n != INVALID; ++n) {

  1861         sum += nodeValue(n);

  1862       }

  1863       for (int i = 0; i < blossomNum(); ++i) {

  1864         sum += blossomValue(i) * (blossomSize(i) / 2);

  1865       }

  1866       return sum;

  1867     }

  1869     /// \brief Returns the value of the node.

  1870     ///

  1871     /// Returns the the value of the node.

  1872     Value nodeValue(const Node& n) const {

  1873       return (*_node_potential)[n];

  1874     }

  1876     /// \brief Returns the number of the blossoms in the basis.

  1877     ///

  1878     /// Returns the number of the blossoms in the basis.

  1879     /// \see BlossomIt

  1880     int blossomNum() const {

  1881       return _blossom_potential.size();

  1882     }

  1885     /// \brief Returns the number of the nodes in the blossom.

  1886     ///

  1887     /// Returns the number of the nodes in the blossom.

  1888     int blossomSize(int k) const {

  1889       return _blossom_potential[k].end - _blossom_potential[k].begin;

  1890     }

  1892     /// \brief Returns the value of the blossom.

  1893     ///

  1894     /// Returns the the value of the blossom.

  1895     /// \see BlossomIt

  1896     Value blossomValue(int k) const {

  1897       return _blossom_potential[k].value;

  1898     }

  1900     /// \brief Iterator for obtaining the nodes of the blossom.

  1901     ///

  1902     /// Iterator for obtaining the nodes of the blossom. This class

  1903     /// provides a common lemon style iterator for listing a

  1904     /// subset of the nodes.

  1905     class BlossomIt {

  1906     public:

  1908       /// \brief Constructor.

  1909       ///

  1910       /// Constructor to get the nodes of the variable.

  1911       BlossomIt(const MaxWeightedMatching& algorithm, int variable)

  1912         : _algorithm(&algorithm)

  1913       {

  1914         _index = _algorithm->_blossom_potential[variable].begin;

  1915         _last = _algorithm->_blossom_potential[variable].end;

  1916       }

  1918       /// \brief Conversion to node.

  1919       ///

  1920       /// Conversion to node.

  1921       operator Node() const {

  1922         return _algorithm->_blossom_node_list[_index];

  1923       }

  1925       /// \brief Increment operator.

  1926       ///

  1927       /// Increment operator.

  1928       BlossomIt& operator++() {

  1929         ++_index;

  1930         return *this;

  1931       }

  1933       /// \brief Validity checking

  1934       ///

  1935       /// Checks whether the iterator is invalid.

  1936       bool operator==(Invalid) const { return _index == _last; }

  1938       /// \brief Validity checking

  1939       ///

  1940       /// Checks whether the iterator is valid.

  1941       bool operator!=(Invalid) const { return _index != _last; }

  1943     private:

  1944       const MaxWeightedMatching* _algorithm;

  1945       int _last;

  1946       int _index;

  1947     };

  1949     /// @}

  1951   };

  1953   /// \ingroup matching

  1954   ///

  1955   /// \brief Weighted perfect matching in general graphs

  1956   ///

  1957   /// This class provides an efficient implementation of Edmond's

  1958   /// maximum weighted perfect matching algorithm. The implementation

  1959   /// is based on extensive use of priority queues and provides

  1960   /// \f$O(nm\log(n))\f$ time complexity.

  1961   ///

  1962   /// The maximum weighted matching problem is to find undirected

  1963   /// edges in the graph with maximum overall weight and no two of

  1964   /// them shares their ends and covers all nodes. The problem can be

  1965   /// formulated with the following linear program.

  1966   /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]

  1967   /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}

  1968       \quad \forall B\in\mathcal{O}\f] */

  1969   /// \f[x_e \ge 0\quad \forall e\in E\f]

  1970   /// \f[\max \sum_{e\in E}x_ew_e\f]

  1971   /// where \f$\delta(X)\f$ is the set of edges incident to a node in

  1972   /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in

  1973   /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality

  1974   /// subsets of the nodes.

  1975   ///

  1976   /// The algorithm calculates an optimal matching and a proof of the

  1977   /// optimality. The solution of the dual problem can be used to check

  1978   /// the result of the algorithm. The dual linear problem is the

  1979   /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge

  1980       w_{uv} \quad \forall uv\in E\f] */

  1981   /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]

  1982   /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}

  1983       \frac{\vert B \vert - 1}{2}z_B\f] */

  1984   ///

  1985   /// The algorithm can be executed with \c run() or the \c init() and

  1986   /// then the \c start() member functions. After it the matching can

  1987   /// be asked with \c matching() or mate() functions. The dual

  1988   /// solution can be get with \c nodeValue(), \c blossomNum() and \c

  1989   /// blossomValue() members and \ref MaxWeightedMatching::BlossomIt

  1990   /// "BlossomIt" nested class which is able to iterate on the nodes

  1991   /// of a blossom. If the value type is integral then the dual

  1992   /// solution is multiplied by \ref MaxWeightedMatching::dualScale "4".

  1993   template <typename _Graph,

  1994             typename _WeightMap = typename _Graph::template EdgeMap<int> >

  1995   class MaxWeightedPerfectMatching {

  1996   public:

  1998     typedef _Graph Graph;

  1999     typedef _WeightMap WeightMap;

  2000     typedef typename WeightMap::Value Value;

  2002     /// \brief Scaling factor for dual solution

  2003     ///

  2004     /// Scaling factor for dual solution, it is equal to 4 or 1

  2005     /// according to the value type.

  2006     static const int dualScale =

  2007       std::numeric_limits<Value>::is_integer ? 4 : 1;

  2009     typedef typename Graph::template NodeMap<typename Graph::Arc>

  2010     MatchingMap;

  2012   private:

  2014     TEMPLATE_GRAPH_TYPEDEFS(Graph);

  2016     typedef typename Graph::template NodeMap<Value> NodePotential;

  2017     typedef std::vector<Node> BlossomNodeList;

  2019     struct BlossomVariable {

  2020       int begin, end;

  2021       Value value;

  2023       BlossomVariable(int _begin, int _end, Value _value)

  2024         : begin(_begin), end(_end), value(_value) {}

  2026     };

  2028     typedef std::vector<BlossomVariable> BlossomPotential;

  2030     const Graph& _graph;

  2031     const WeightMap& _weight;

  2033     MatchingMap* _matching;

  2035     NodePotential* _node_potential;

  2037     BlossomPotential _blossom_potential;

  2038     BlossomNodeList _blossom_node_list;

  2040     int _node_num;

  2041     int _blossom_num;

  2043     typedef RangeMap<int> IntIntMap;

  2045     enum Status {

  2046       EVEN = -1, MATCHED = 0, ODD = 1

  2047     };

  2049     typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;

  2050     struct BlossomData {

  2051       int tree;

  2052       Status status;

  2053       Arc pred, next;

  2054       Value pot, offset;

  2055     };

  2057     IntNodeMap *_blossom_index;

  2058     BlossomSet *_blossom_set;

  2059     RangeMap<BlossomData>* _blossom_data;

  2061     IntNodeMap *_node_index;

  2062     IntArcMap *_node_heap_index;

  2064     struct NodeData {

  2066       NodeData(IntArcMap& node_heap_index)

  2067         : heap(node_heap_index) {}

  2069       int blossom;

  2070       Value pot;

  2071       BinHeap<Value, IntArcMap> heap;

  2072       std::map<int, Arc> heap_index;

  2074       int tree;

  2075     };

  2077     RangeMap<NodeData>* _node_data;

  2079     typedef ExtendFindEnum<IntIntMap> TreeSet;

  2081     IntIntMap *_tree_set_index;

  2082     TreeSet *_tree_set;

  2084     IntIntMap *_delta2_index;

  2085     BinHeap<Value, IntIntMap> *_delta2;

  2087     IntEdgeMap *_delta3_index;

  2088     BinHeap<Value, IntEdgeMap> *_delta3;

  2090     IntIntMap *_delta4_index;

  2091     BinHeap<Value, IntIntMap> *_delta4;

  2093     Value _delta_sum;

  2095     void createStructures() {

  2096       _node_num = countNodes(_graph);

  2097       _blossom_num = _node_num * 3 / 2;

  2099       if (!_matching) {

  2100         _matching = new MatchingMap(_graph);

  2101       }

  2102       if (!_node_potential) {

  2103         _node_potential = new NodePotential(_graph);

  2104       }

  2105       if (!_blossom_set) {

  2106         _blossom_index = new IntNodeMap(_graph);

  2107         _blossom_set = new BlossomSet(*_blossom_index);

  2108         _blossom_data = new RangeMap<BlossomData>(_blossom_num);

  2109       }

  2111       if (!_node_index) {

  2112         _node_index = new IntNodeMap(_graph);

  2113         _node_heap_index = new IntArcMap(_graph);

  2114         _node_data = new RangeMap<NodeData>(_node_num,

  2115                                             NodeData(*_node_heap_index));

  2116       }

  2118       if (!_tree_set) {

  2119         _tree_set_index = new IntIntMap(_blossom_num);

  2120         _tree_set = new TreeSet(*_tree_set_index);

  2121       }

  2122       if (!_delta2) {

  2123         _delta2_index = new IntIntMap(_blossom_num);

  2124         _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);

  2125       }

  2126       if (!_delta3) {

  2127         _delta3_index = new IntEdgeMap(_graph);

  2128         _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);

  2129       }

  2130       if (!_delta4) {

  2131         _delta4_index = new IntIntMap(_blossom_num);

  2132         _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);

  2133       }

  2134     }

  2136     void destroyStructures() {

  2137       _node_num = countNodes(_graph);

  2138       _blossom_num = _node_num * 3 / 2;

  2140       if (_matching) {

  2141         delete _matching;

  2142       }

  2143       if (_node_potential) {

  2144         delete _node_potential;

  2145       }

  2146       if (_blossom_set) {

  2147         delete _blossom_index;

  2148         delete _blossom_set;

  2149         delete _blossom_data;

  2150       }

  2152       if (_node_index) {

  2153         delete _node_index;

  2154         delete _node_heap_index;

  2155         delete _node_data;

  2156       }

  2158       if (_tree_set) {

  2159         delete _tree_set_index;

  2160         delete _tree_set;

  2161       }

  2162       if (_delta2) {

  2163         delete _delta2_index;

  2164         delete _delta2;

  2165       }

  2166       if (_delta3) {

  2167         delete _delta3_index;

  2168         delete _delta3;

  2169       }

  2170       if (_delta4) {

  2171         delete _delta4_index;

  2172         delete _delta4;

  2173       }

  2174     }

  2176     void matchedToEven(int blossom, int tree) {

  2177       if (_delta2->state(blossom) == _delta2->IN_HEAP) {

  2178         _delta2->erase(blossom);

  2179       }

  2181       if (!_blossom_set->trivial(blossom)) {

  2182         (*_blossom_data)[blossom].pot -=

  2183           2 * (_delta_sum - (*_blossom_data)[blossom].offset);

  2184       }

  2186       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);

  2187            n != INVALID; ++n) {

  2189         _blossom_set->increase(n, std::numeric_limits<Value>::max());

  2190         int ni = (*_node_index)[n];

  2192         (*_node_data)[ni].heap.clear();

  2193         (*_node_data)[ni].heap_index.clear();

  2195         (*_node_data)[ni].pot += _delta_sum - (*_blossom_data)[blossom].offset;

  2197         for (InArcIt e(_graph, n); e != INVALID; ++e) {

  2198           Node v = _graph.source(e);

  2199           int vb = _blossom_set->find(v);

  2200           int vi = (*_node_index)[v];

  2202           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -

  2203             dualScale * _weight[e];

  2205           if ((*_blossom_data)[vb].status == EVEN) {

  2206             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {

  2207               _delta3->push(e, rw / 2);

  2208             }

  2209           } else {

  2210             typename std::map<int, Arc>::iterator it =

  2211               (*_node_data)[vi].heap_index.find(tree);

  2213             if (it != (*_node_data)[vi].heap_index.end()) {

  2214               if ((*_node_data)[vi].heap[it->second] > rw) {

  2215                 (*_node_data)[vi].heap.replace(it->second, e);

  2216                 (*_node_data)[vi].heap.decrease(e, rw);

  2217                 it->second = e;

  2218               }

  2219             } else {

  2220               (*_node_data)[vi].heap.push(e, rw);

  2221               (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));

  2222             }

  2224             if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {

  2225               _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());

  2227               if ((*_blossom_data)[vb].status == MATCHED) {

  2228                 if (_delta2->state(vb) != _delta2->IN_HEAP) {

  2229                   _delta2->push(vb, _blossom_set->classPrio(vb) -

  2230                                (*_blossom_data)[vb].offset);

  2231                 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -

  2232                            (*_blossom_data)[vb].offset){

  2233                   _delta2->decrease(vb, _blossom_set->classPrio(vb) -

  2234                                    (*_blossom_data)[vb].offset);

  2235                 }

  2236               }

  2237             }

  2238           }

  2239         }

  2240       }

  2241       (*_blossom_data)[blossom].offset = 0;

  2242     }

  2244     void matchedToOdd(int blossom) {

  2245       if (_delta2->state(blossom) == _delta2->IN_HEAP) {

  2246         _delta2->erase(blossom);

  2247       }

  2248       (*_blossom_data)[blossom].offset += _delta_sum;

  2249       if (!_blossom_set->trivial(blossom)) {

  2250         _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +

  2251                      (*_blossom_data)[blossom].offset);

  2252       }

  2253     }

  2255     void evenToMatched(int blossom, int tree) {

  2256       if (!_blossom_set->trivial(blossom)) {

  2257         (*_blossom_data)[blossom].pot += 2 * _delta_sum;

  2258       }

  2260       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);

  2261            n != INVALID; ++n) {

  2262         int ni = (*_node_index)[n];

  2263         (*_node_data)[ni].pot -= _delta_sum;

  2265         for (InArcIt e(_graph, n); e != INVALID; ++e) {

  2266           Node v = _graph.source(e);

  2267           int vb = _blossom_set->find(v);

  2268           int vi = (*_node_index)[v];

  2270           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -

  2271             dualScale * _weight[e];

  2273           if (vb == blossom) {

  2274             if (_delta3->state(e) == _delta3->IN_HEAP) {

  2275               _delta3->erase(e);

  2276             }

  2277           } else if ((*_blossom_data)[vb].status == EVEN) {

  2279             if (_delta3->state(e) == _delta3->IN_HEAP) {

  2280               _delta3->erase(e);

  2281             }

  2283             int vt = _tree_set->find(vb);

  2285             if (vt != tree) {

  2287               Arc r = _graph.oppositeArc(e);

  2289               typename std::map<int, Arc>::iterator it =

  2290                 (*_node_data)[ni].heap_index.find(vt);

  2292               if (it != (*_node_data)[ni].heap_index.end()) {

  2293                 if ((*_node_data)[ni].heap[it->second] > rw) {

  2294                   (*_node_data)[ni].heap.replace(it->second, r);

  2295                   (*_node_data)[ni].heap.decrease(r, rw);

  2296                   it->second = r;

  2297                 }

  2298               } else {

  2299                 (*_node_data)[ni].heap.push(r, rw);

  2300                 (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));

  2301               }

  2303               if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {

  2304                 _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());

  2306                 if (_delta2->state(blossom) != _delta2->IN_HEAP) {

  2307                   _delta2->push(blossom, _blossom_set->classPrio(blossom) -

  2308                                (*_blossom_data)[blossom].offset);

  2309                 } else if ((*_delta2)[blossom] >

  2310                            _blossom_set->classPrio(blossom) -

  2311                            (*_blossom_data)[blossom].offset){

  2312                   _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -

  2313                                    (*_blossom_data)[blossom].offset);

  2314                 }

  2315               }

  2316             }

  2317           } else {

  2319             typename std::map<int, Arc>::iterator it =

  2320               (*_node_data)[vi].heap_index.find(tree);

  2322             if (it != (*_node_data)[vi].heap_index.end()) {

  2323               (*_node_data)[vi].heap.erase(it->second);

  2324               (*_node_data)[vi].heap_index.erase(it);

  2325               if ((*_node_data)[vi].heap.empty()) {

  2326                 _blossom_set->increase(v, std::numeric_limits<Value>::max());

  2327               } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {

  2328                 _blossom_set->increase(v, (*_node_data)[vi].heap.prio());

  2329               }

  2331               if ((*_blossom_data)[vb].status == MATCHED) {

  2332                 if (_blossom_set->classPrio(vb) ==

  2333                     std::numeric_limits<Value>::max()) {

  2334                   _delta2->erase(vb);

  2335                 } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -

  2336                            (*_blossom_data)[vb].offset) {

  2337                   _delta2->increase(vb, _blossom_set->classPrio(vb) -

  2338                                    (*_blossom_data)[vb].offset);

  2339                 }

  2340               }

  2341             }

  2342           }

  2343         }

  2344       }

  2345     }

  2347     void oddToMatched(int blossom) {

  2348       (*_blossom_data)[blossom].offset -= _delta_sum;

  2350       if (_blossom_set->classPrio(blossom) !=

  2351           std::numeric_limits<Value>::max()) {

  2352         _delta2->push(blossom, _blossom_set->classPrio(blossom) -

  2353                        (*_blossom_data)[blossom].offset);

  2354       }

  2356       if (!_blossom_set->trivial(blossom)) {

  2357         _delta4->erase(blossom);

  2358       }

  2359     }

  2361     void oddToEven(int blossom, int tree) {

  2362       if (!_blossom_set->trivial(blossom)) {

  2363         _delta4->erase(blossom);

  2364         (*_blossom_data)[blossom].pot -=

  2365           2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);

  2366       }

  2368       for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);

  2369            n != INVALID; ++n) {

  2370         int ni = (*_node_index)[n];

  2372         _blossom_set->increase(n, std::numeric_limits<Value>::max());

  2374         (*_node_data)[ni].heap.clear();

  2375         (*_node_data)[ni].heap_index.clear();

  2376         (*_node_data)[ni].pot +=

  2377           2 * _delta_sum - (*_blossom_data)[blossom].offset;

  2379         for (InArcIt e(_graph, n); e != INVALID; ++e) {

  2380           Node v = _graph.source(e);

  2381           int vb = _blossom_set->find(v);

  2382           int vi = (*_node_index)[v];

  2384           Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -

  2385             dualScale * _weight[e];

  2387           if ((*_blossom_data)[vb].status == EVEN) {

  2388             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {

  2389               _delta3->push(e, rw / 2);

  2390             }

  2391           } else {

  2393             typename std::map<int, Arc>::iterator it =

  2394               (*_node_data)[vi].heap_index.find(tree);

  2396             if (it != (*_node_data)[vi].heap_index.end()) {

  2397               if ((*_node_data)[vi].heap[it->second] > rw) {

  2398                 (*_node_data)[vi].heap.replace(it->second, e);

  2399                 (*_node_data)[vi].heap.decrease(e, rw);

  2400                 it->second = e;

  2401               }

  2402             } else {

  2403               (*_node_data)[vi].heap.push(e, rw);

  2404               (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));

  2405             }

  2407             if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {

  2408               _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());

  2410               if ((*_blossom_data)[vb].status == MATCHED) {

  2411                 if (_delta2->state(vb) != _delta2->IN_HEAP) {

  2412                   _delta2->push(vb, _blossom_set->classPrio(vb) -

  2413                                (*_blossom_data)[vb].offset);

  2414                 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -

  2415                            (*_blossom_data)[vb].offset) {

  2416                   _delta2->decrease(vb, _blossom_set->classPrio(vb) -

  2417                                    (*_blossom_data)[vb].offset);

  2418                 }

  2419               }

  2420             }

  2421           }

  2422         }

  2423       }

  2424       (*_blossom_data)[blossom].offset = 0;

  2425     }

  2427     void alternatePath(int even, int tree) {

  2428       int odd;

  2430       evenToMatched(even, tree);

  2431       (*_blossom_data)[even].status = MATCHED;

  2433       while ((*_blossom_data)[even].pred != INVALID) {

  2434         odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));

  2435         (*_blossom_data)[odd].status = MATCHED;

  2436         oddToMatched(odd);

  2437         (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;

  2439         even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));

  2440         (*_blossom_data)[even].status = MATCHED;

  2441         evenToMatched(even, tree);

  2442         (*_blossom_data)[even].next =

  2443           _graph.oppositeArc((*_blossom_data)[odd].pred);

  2444       }

  2446     }

  2448     void destroyTree(int tree) {

  2449       for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {

  2450         if ((*_blossom_data)[b].status == EVEN) {

  2451           (*_blossom_data)[b].status = MATCHED;

  2452           evenToMatched(b, tree);

  2453         } else if ((*_blossom_data)[b].status == ODD) {

  2454           (*_blossom_data)[b].status = MATCHED;

  2455           oddToMatched(b);

  2456         }

  2457       }

  2458       _tree_set->eraseClass(tree);

  2459     }

  2461     void augmentOnEdge(const Edge& edge) {

  2463       int left = _blossom_set->find(_graph.u(edge));

  2464       int right = _blossom_set->find(_graph.v(edge));

  2466       int left_tree = _tree_set->find(left);

  2467       alternatePath(left, left_tree);

  2468       destroyTree(left_tree);

  2470       int right_tree = _tree_set->find(right);

  2471       alternatePath(right, right_tree);

  2472       destroyTree(right_tree);

  2474       (*_blossom_data)[left].next = _graph.direct(edge, true);

  2475       (*_blossom_data)[right].next = _graph.direct(edge, false);

  2476     }

  2478     void extendOnArc(const Arc& arc) {

  2479       int base = _blossom_set->find(_graph.target(arc));

  2480       int tree = _tree_set->find(base);

  2482       int odd = _blossom_set->find(_graph.source(arc));

  2483       _tree_set->insert(odd, tree);

  2484       (*_blossom_data)[odd].status = ODD;

  2485       matchedToOdd(odd);

  2486       (*_blossom_data)[odd].pred = arc;

  2488       int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));

  2489       (*_blossom_data)[even].pred = (*_blossom_data)[even].next;

  2490       _tree_set->insert(even, tree);

  2491       (*_blossom_data)[even].status = EVEN;

  2492       matchedToEven(even, tree);

  2493     }

  2495     void shrinkOnEdge(const Edge& edge, int tree) {

  2496       int nca = -1;

  2497       std::vector<int> left_path, right_path;

  2499       {

  2500         std::set<int> left_set, right_set;

  2501         int left = _blossom_set->find(_graph.u(edge));

  2502         left_path.push_back(left);

  2503         left_set.insert(left);

  2505         int right = _blossom_set->find(_graph.v(edge));

  2506         right_path.push_back(right);

  2507         right_set.insert(right);

  2509         while (true) {

  2511           if ((*_blossom_data)[left].pred == INVALID) break;

  2513           left =

  2514             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));

  2515           left_path.push_back(left);

  2516           left =

  2517             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));

  2518           left_path.push_back(left);

  2520           left_set.insert(left);

  2522           if (right_set.find(left) != right_set.end()) {

  2523             nca = left;

  2524             break;

  2525           }

  2527           if ((*_blossom_data)[right].pred == INVALID) break;

  2529           right =

  2530             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));

  2531           right_path.push_back(right);

  2532           right =

  2533             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));

  2534           right_path.push_back(right);

  2536           right_set.insert(right);

  2538           if (left_set.find(right) != left_set.end()) {

  2539             nca = right;

  2540             break;

  2541           }

  2543         }

  2545         if (nca == -1) {

  2546           if ((*_blossom_data)[left].pred == INVALID) {

  2547             nca = right;

  2548             while (left_set.find(nca) == left_set.end()) {

  2549               nca =

  2550                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));

  2551               right_path.push_back(nca);

  2552               nca =

  2553                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));

  2554               right_path.push_back(nca);

  2555             }

  2556           } else {

  2557             nca = left;

  2558             while (right_set.find(nca) == right_set.end()) {

  2559               nca =

  2560                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));

  2561               left_path.push_back(nca);

  2562               nca =

  2563                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));

  2564               left_path.push_back(nca);

  2565             }

  2566           }

  2567         }

  2568       }

  2570       std::vector<int> subblossoms;

  2571       Arc prev;

  2573       prev = _graph.direct(edge, true);

  2574       for (int i = 0; left_path[i] != nca; i += 2) {

  2575         subblossoms.push_back(left_path[i]);

  2576         (*_blossom_data)[left_path[i]].next = prev;

  2577         _tree_set->erase(left_path[i]);

  2579         subblossoms.push_back(left_path[i + 1]);

  2580         (*_blossom_data)[left_path[i + 1]].status = EVEN;

  2581         oddToEven(left_path[i + 1], tree);

  2582         _tree_set->erase(left_path[i + 1]);

  2583         prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);

  2584       }

  2586       int k = 0;

  2587       while (right_path[k] != nca) ++k;

  2589       subblossoms.push_back(nca);

  2590       (*_blossom_data)[nca].next = prev;

  2592       for (int i = k - 2; i >= 0; i -= 2) {

  2593         subblossoms.push_back(right_path[i + 1]);

  2594         (*_blossom_data)[right_path[i + 1]].status = EVEN;

  2595         oddToEven(right_path[i + 1], tree);

  2596         _tree_set->erase(right_path[i + 1]);

  2598         (*_blossom_data)[right_path[i + 1]].next =

  2599           (*_blossom_data)[right_path[i + 1]].pred;

  2601         subblossoms.push_back(right_path[i]);

  2602         _tree_set->erase(right_path[i]);

  2603       }

  2605       int surface =

  2606         _blossom_set->join(subblossoms.begin(), subblossoms.end());

  2608       for (int i = 0; i < int(subblossoms.size()); ++i) {

  2609         if (!_blossom_set->trivial(subblossoms[i])) {

  2610           (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;

  2611         }

  2612         (*_blossom_data)[subblossoms[i]].status = MATCHED;

  2613       }

  2615       (*_blossom_data)[surface].pot = -2 * _delta_sum;

  2616       (*_blossom_data)[surface].offset = 0;

  2617       (*_blossom_data)[surface].status = EVEN;

  2618       (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;

  2619       (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;

  2621       _tree_set->insert(surface, tree);

  2622       _tree_set->erase(nca);

  2623     }

  2625     void splitBlossom(int blossom) {

  2626       Arc next = (*_blossom_data)[blossom].next;

  2627       Arc pred = (*_blossom_data)[blossom].pred;

  2629       int tree = _tree_set->find(blossom);

  2631       (*_blossom_data)[blossom].status = MATCHED;

  2632       oddToMatched(blossom);

  2633       if (_delta2->state(blossom) == _delta2->IN_HEAP) {

  2634         _delta2->erase(blossom);

  2635       }

  2637       std::vector<int> subblossoms;

  2638       _blossom_set->split(blossom, std::back_inserter(subblossoms));

  2640       Value offset = (*_blossom_data)[blossom].offset;

  2641       int b = _blossom_set->find(_graph.source(pred));

  2642       int d = _blossom_set->find(_graph.source(next));

  2644       int ib = -1, id = -1;

  2645       for (int i = 0; i < int(subblossoms.size()); ++i) {

  2646         if (subblossoms[i] == b) ib = i;

  2647         if (subblossoms[i] == d) id = i;

  2649         (*_blossom_data)[subblossoms[i]].offset = offset;

  2650         if (!_blossom_set->trivial(subblossoms[i])) {

  2651           (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;

  2652         }

  2653         if (_blossom_set->classPrio(subblossoms[i]) !=

  2654             std::numeric_limits<Value>::max()) {

  2655           _delta2->push(subblossoms[i],

  2656                         _blossom_set->classPrio(subblossoms[i]) -

  2657                         (*_blossom_data)[subblossoms[i]].offset);

  2658         }

  2659       }

  2661       if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {

  2662         for (int i = (id + 1) % subblossoms.size();

  2663              i != ib; i = (i + 2) % subblossoms.size()) {

  2664           int sb = subblossoms[i];

  2665           int tb = subblossoms[(i + 1) % subblossoms.size()];

  2666           (*_blossom_data)[sb].next =

  2667             _graph.oppositeArc((*_blossom_data)[tb].next);

  2668         }

  2670         for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {

  2671           int sb = subblossoms[i];

  2672           int tb = subblossoms[(i + 1) % subblossoms.size()];

  2673           int ub = subblossoms[(i + 2) % subblossoms.size()];

  2675           (*_blossom_data)[sb].status = ODD;

  2676           matchedToOdd(sb);

  2677           _tree_set->insert(sb, tree);

  2678           (*_blossom_data)[sb].pred = pred;

  2679           (*_blossom_data)[sb].next =

  2680                            _graph.oppositeArc((*_blossom_data)[tb].next);

  2682           pred = (*_blossom_data)[ub].next;

  2684           (*_blossom_data)[tb].status = EVEN;

  2685           matchedToEven(tb, tree);

  2686           _tree_set->insert(tb, tree);

  2687           (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;

  2688         }

  2690         (*_blossom_data)[subblossoms[id]].status = ODD;

  2691         matchedToOdd(subblossoms[id]);

  2692         _tree_set->insert(subblossoms[id], tree);

  2693         (*_blossom_data)[subblossoms[id]].next = next;

  2694         (*_blossom_data)[subblossoms[id]].pred = pred;

  2696       } else {

  2698         for (int i = (ib + 1) % subblossoms.size();

  2699              i != id; i = (i + 2) % subblossoms.size()) {

  2700           int sb = subblossoms[i];

  2701           int tb = subblossoms[(i + 1) % subblossoms.size()];

  2702           (*_blossom_data)[sb].next =

  2703             _graph.oppositeArc((*_blossom_data)[tb].next);

  2704         }

  2706         for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {

  2707           int sb = subblossoms[i];

  2708           int tb = subblossoms[(i + 1) % subblossoms.size()];

  2709           int ub = subblossoms[(i + 2) % subblossoms.size()];

  2711           (*_blossom_data)[sb].status = ODD;

  2712           matchedToOdd(sb);

  2713           _tree_set->insert(sb, tree);

  2714           (*_blossom_data)[sb].next = next;

  2715           (*_blossom_data)[sb].pred =

  2716             _graph.oppositeArc((*_blossom_data)[tb].next);

  2718           (*_blossom_data)[tb].status = EVEN;

  2719           matchedToEven(tb, tree);

  2720           _tree_set->insert(tb, tree);

  2721           (*_blossom_data)[tb].pred =

  2722             (*_blossom_data)[tb].next =

  2723             _graph.oppositeArc((*_blossom_data)[ub].next);

  2724           next = (*_blossom_data)[ub].next;

  2725         }

  2727         (*_blossom_data)[subblossoms[ib]].status = ODD;

  2728         matchedToOdd(subblossoms[ib]);

  2729         _tree_set->insert(subblossoms[ib], tree);

  2730         (*_blossom_data)[subblossoms[ib]].next = next;

  2731         (*_blossom_data)[subblossoms[ib]].pred = pred;

  2732       }

  2733       _tree_set->erase(blossom);

  2734     }

  2736     void extractBlossom(int blossom, const Node& base, const Arc& matching) {

  2737       if (_blossom_set->trivial(blossom)) {

  2738         int bi = (*_node_index)[base];

  2739         Value pot = (*_node_data)[bi].pot;

  2741         _matching->set(base, matching);

  2742         _blossom_node_list.push_back(base);

  2743         _node_potential->set(base, pot);

  2744       } else {

  2746         Value pot = (*_blossom_data)[blossom].pot;

  2747         int bn = _blossom_node_list.size();

  2749         std::vector<int> subblossoms;

  2750         _blossom_set->split(blossom, std::back_inserter(subblossoms));

  2751         int b = _blossom_set->find(base);

  2752         int ib = -1;

  2753         for (int i = 0; i < int(subblossoms.size()); ++i) {

  2754           if (subblossoms[i] == b) { ib = i; break; }

  2755         }

  2757         for (int i = 1; i < int(subblossoms.size()); i += 2) {

  2758           int sb = subblossoms[(ib + i) % subblossoms.size()];

  2759           int tb = subblossoms[(ib + i + 1) % subblossoms.size()];

  2761           Arc m = (*_blossom_data)[tb].next;

  2762           extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));

  2763           extractBlossom(tb, _graph.source(m), m);

  2764         }

  2765         extractBlossom(subblossoms[ib], base, matching);

  2767         int en = _blossom_node_list.size();

  2769         _blossom_potential.push_back(BlossomVariable(bn, en, pot));

  2770       }

  2771     }

  2773     void extractMatching() {

  2774       std::vector<int> blossoms;

  2775       for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {

  2776         blossoms.push_back(c);

  2777       }

  2779       for (int i = 0; i < int(blossoms.size()); ++i) {

  2781         Value offset = (*_blossom_data)[blossoms[i]].offset;

  2782         (*_blossom_data)[blossoms[i]].pot += 2 * offset;

  2783         for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);

  2784              n != INVALID; ++n) {

  2785           (*_node_data)[(*_node_index)[n]].pot -= offset;

  2786         }

  2788         Arc matching = (*_blossom_data)[blossoms[i]].next;

  2789         Node base = _graph.source(matching);

  2790         extractBlossom(blossoms[i], base, matching);

  2791       }

  2792     }

  2794   public:

  2796     /// \brief Constructor

  2797     ///

  2798     /// Constructor.

  2799     MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)

  2800       : _graph(graph), _weight(weight), _matching(0),

  2801         _node_potential(0), _blossom_potential(), _blossom_node_list(),

  2802         _node_num(0), _blossom_num(0),

  2804         _blossom_index(0), _blossom_set(0), _blossom_data(0),

  2805         _node_index(0), _node_heap_index(0), _node_data(0),

  2806         _tree_set_index(0), _tree_set(0),

  2808         _delta2_index(0), _delta2(0),

  2809         _delta3_index(0), _delta3(0),

  2810         _delta4_index(0), _delta4(0),

  2812         _delta_sum() {}

  2814     ~MaxWeightedPerfectMatching() {

  2815       destroyStructures();

  2816     }

  2818     /// \name Execution control

  2819     /// The simplest way to execute the algorithm is to use the

  2820     /// \c run() member function.

  2822     ///@{

  2824     /// \brief Initialize the algorithm

  2825     ///

  2826     /// Initialize the algorithm

  2827     void init() {

  2828       createStructures();

  2830       for (ArcIt e(_graph); e != INVALID; ++e) {

  2831         _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);

  2832       }

  2833       for (EdgeIt e(_graph); e != INVALID; ++e) {

  2834         _delta3_index->set(e, _delta3->PRE_HEAP);

  2835       }

  2836       for (int i = 0; i < _blossom_num; ++i) {

  2837         _delta2_index->set(i, _delta2->PRE_HEAP);

  2838         _delta4_index->set(i, _delta4->PRE_HEAP);

  2839       }

  2841       int index = 0;

  2842       for (NodeIt n(_graph); n != INVALID; ++n) {

  2843         Value max = - std::numeric_limits<Value>::max();

  2844         for (OutArcIt e(_graph, n); e != INVALID; ++e) {

  2845           if (_graph.target(e) == n) continue;

  2846           if ((dualScale * _weight[e]) / 2 > max) {

  2847             max = (dualScale * _weight[e]) / 2;

  2848           }

  2849         }

  2850         _node_index->set(n, index);

  2851         (*_node_data)[index].pot = max;

  2852         int blossom =

  2853           _blossom_set->insert(n, std::numeric_limits<Value>::max());

  2855         _tree_set->insert(blossom);

  2857         (*_blossom_data)[blossom].status = EVEN;

  2858         (*_blossom_data)[blossom].pred = INVALID;

  2859         (*_blossom_data)[blossom].next = INVALID;

  2860         (*_blossom_data)[blossom].pot = 0;

  2861         (*_blossom_data)[blossom].offset = 0;

  2862         ++index;

  2863       }

  2864       for (EdgeIt e(_graph); e != INVALID; ++e) {

  2865         int si = (*_node_index)[_graph.u(e)];

  2866         int ti = (*_node_index)[_graph.v(e)];

  2867         if (_graph.u(e) != _graph.v(e)) {

  2868           _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -

  2869                             dualScale * _weight[e]) / 2);

  2870         }

  2871       }

  2872     }

  2874     /// \brief Starts the algorithm

  2875     ///

  2876     /// Starts the algorithm

  2877     bool start() {

  2878       enum OpType {

  2879         D2, D3, D4

  2880       };

  2882       int unmatched = _node_num;

  2883       while (unmatched > 0) {

  2884         Value d2 = !_delta2->empty() ?

  2885           _delta2->prio() : std::numeric_limits<Value>::max();

  2887         Value d3 = !_delta3->empty() ?

  2888           _delta3->prio() : std::numeric_limits<Value>::max();

  2890         Value d4 = !_delta4->empty() ?

  2891           _delta4->prio() : std::numeric_limits<Value>::max();

  2893         _delta_sum = d2; OpType ot = D2;

  2894         if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }

  2895         if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }

  2897         if (_delta_sum == std::numeric_limits<Value>::max()) {

  2898           return false;

  2899         }

  2901         switch (ot) {

  2902         case D2:

  2903           {

  2904             int blossom = _delta2->top();

  2905             Node n = _blossom_set->classTop(blossom);

  2906             Arc e = (*_node_data)[(*_node_index)[n]].heap.top();

  2907             extendOnArc(e);

  2908           }

  2909           break;

  2910         case D3:

  2911           {

  2912             Edge e = _delta3->top();

  2914             int left_blossom = _blossom_set->find(_graph.u(e));

  2915             int right_blossom = _blossom_set->find(_graph.v(e));

  2917             if (left_blossom == right_blossom) {

  2918               _delta3->pop();

  2919             } else {

  2920               int left_tree = _tree_set->find(left_blossom);

  2921               int right_tree = _tree_set->find(right_blossom);

  2923               if (left_tree == right_tree) {

  2924                 shrinkOnEdge(e, left_tree);

  2925               } else {

  2926                 augmentOnEdge(e);

  2927                 unmatched -= 2;

  2928               }

  2929             }

  2930           } break;

  2931         case D4:

  2932           splitBlossom(_delta4->top());

  2933           break;

  2934         }

  2935       }

  2936       extractMatching();

  2937       return true;

  2938     }

  2940     /// \brief Runs %MaxWeightedPerfectMatching algorithm.

  2941     ///

  2942     /// This method runs the %MaxWeightedPerfectMatching algorithm.

  2943     ///

  2944     /// \note mwm.run() is just a shortcut of the following code.

  2945     /// \code

  2946     ///   mwm.init();

  2947     ///   mwm.start();

  2948     /// \endcode

  2949     bool run() {

  2950       init();

  2951       return start();

  2952     }

  2954     /// @}

  2956     /// \name Primal solution

  2957     /// Functions to get the primal solution, ie. the matching.

  2959     /// @{

  2961     /// \brief Returns the matching value.

  2962     ///

  2963     /// Returns the matching value.

  2964     Value matchingValue() const {

  2965       Value sum = 0;

  2966       for (NodeIt n(_graph); n != INVALID; ++n) {

  2967         if ((*_matching)[n] != INVALID) {

  2968           sum += _weight[(*_matching)[n]];

  2969         }

  2970       }

  2971       return sum /= 2;

  2972     }

  2974     /// \brief Returns true when the edge is in the matching.

  2975     ///

  2976     /// Returns true when the edge is in the matching.

  2977     bool matching(const Edge& edge) const {

  2978       return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;

  2979     }

  2981     /// \brief Returns the incident matching edge.

  2982     ///

  2983     /// Returns the incident matching arc from given edge.

  2984     Arc matching(const Node& node) const {

  2985       return (*_matching)[node];

  2986     }

  2988     /// \brief Returns the mate of the node.

  2989     ///

  2990     /// Returns the adjancent node in a mathcing arc.

  2991     Node mate(const Node& node) const {

  2992       return _graph.target((*_matching)[node]);

  2993     }

  2995     /// @}

  2997     /// \name Dual solution

  2998     /// Functions to get the dual solution.

  3000     /// @{

  3002     /// \brief Returns the value of the dual solution.

  3003     ///

  3004     /// Returns the value of the dual solution. It should be equal to

  3005     /// the primal value scaled by \ref dualScale "dual scale".

  3006     Value dualValue() const {

  3007       Value sum = 0;

  3008       for (NodeIt n(_graph); n != INVALID; ++n) {

  3009         sum += nodeValue(n);

  3010       }

  3011       for (int i = 0; i < blossomNum(); ++i) {

  3012         sum += blossomValue(i) * (blossomSize(i) / 2);

  3013       }

  3014       return sum;

  3015     }

  3017     /// \brief Returns the value of the node.

  3018     ///

  3019     /// Returns the the value of the node.

  3020     Value nodeValue(const Node& n) const {

  3021       return (*_node_potential)[n];

  3022     }

  3024     /// \brief Returns the number of the blossoms in the basis.

  3025     ///

  3026     /// Returns the number of the blossoms in the basis.

  3027     /// \see BlossomIt

  3028     int blossomNum() const {

  3029       return _blossom_potential.size();

  3030     }

  3033     /// \brief Returns the number of the nodes in the blossom.

  3034     ///

  3035     /// Returns the number of the nodes in the blossom.

  3036     int blossomSize(int k) const {

  3037       return _blossom_potential[k].end - _blossom_potential[k].begin;

  3038     }

  3040     /// \brief Returns the value of the blossom.

  3041     ///

  3042     /// Returns the the value of the blossom.

  3043     /// \see BlossomIt

  3044     Value blossomValue(int k) const {

  3045       return _blossom_potential[k].value;

  3046     }

  3048     /// \brief Iterator for obtaining the nodes of the blossom.

  3049     ///

  3050     /// Iterator for obtaining the nodes of the blossom. This class

  3051     /// provides a common lemon style iterator for listing a

  3052     /// subset of the nodes.

  3053     class BlossomIt {

  3054     public:

  3056       /// \brief Constructor.

  3057       ///

  3058       /// Constructor to get the nodes of the variable.

  3059       BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)

  3060         : _algorithm(&algorithm)

  3061       {

  3062         _index = _algorithm->_blossom_potential[variable].begin;

  3063         _last = _algorithm->_blossom_potential[variable].end;

  3064       }

  3066       /// \brief Conversion to node.

  3067       ///

  3068       /// Conversion to node.

  3069       operator Node() const {

  3070         return _algorithm->_blossom_node_list[_index];

  3071       }

  3073       /// \brief Increment operator.

  3074       ///

  3075       /// Increment operator.

  3076       BlossomIt& operator++() {

  3077         ++_index;

  3078         return *this;

  3079       }

  3081       /// \brief Validity checking

  3082       ///

  3083       /// Checks whether the iterator is invalid.

  3084       bool operator==(Invalid) const { return _index == _last; }

  3086       /// \brief Validity checking

  3087       ///

  3088       /// Checks whether the iterator is valid.

  3089       bool operator!=(Invalid) const { return _index != _last; }

  3091     private:

  3092       const MaxWeightedPerfectMatching* _algorithm;

  3093       int _last;

  3094       int _index;

  3095     };

  3097     /// @}

  3099   };

  3102 } //END OF NAMESPACE LEMON

  3104 #endif //LEMON_MAX_MATCHING_H

author	Alpar Juttner <alpar@cs.elte.hu>
	Mon, 23 Feb 2009 12:33:40 +0000
changeset 577	1c5d6e47921f
parent 444	cace3206223b
child 606	c5fd2d996909
permissions	-rw-r--r--