lemon: lemon/max_matching.h@91d63b8b1a4c

     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-

     2  *

     3  * This file is a part of LEMON, a generic C++ optimization library.

     4  *

     5  * Copyright (C) 2003-2008

     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).

     8  *

     9  * Permission to use, modify and distribute this software is granted

    10  * provided that this copyright notice appears in all copies. For

    11  * precise terms see the accompanying LICENSE file.

    12  *

    13  * This software is provided "AS IS" with no warranty of any kind,

    14  * express or implied, and with no claim as to its suitability for any

    15  * purpose.

    16  *

    17  */

    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 provides Edmonds' alternating forest matching

    43   /// algorithm. The starting matching (if any) can be passed to the

    44   /// algorithm using some of init functions.

    45   ///

    46   /// The dual side of a matching 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 fractor critical components

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

    54   /// unmatched nodes, and if the matching is optimal 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, which

    70     ///shows an upper bound on the size of a maximum matching. The

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

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

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

    74     ///matching.

    75     enum Status {

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

    77     };

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

    81   private:

    83     TEMPLATE_GRAPH_TYPEDEFS(Graph);

    85     typedef UnionFindEnum<IntNodeMap> BlossomSet;

    86     typedef ExtendFindEnum<IntNodeMap> TreeSet;

    87     typedef RangeMap<Node> NodeIntMap;

    88     typedef MatchingMap EarMap;

    89     typedef std::vector<Node> NodeQueue;

    91     const Graph& _graph;

    92     MatchingMap* _matching;

    93     StatusMap* _status;

    95     EarMap* _ear;

    97     IntNodeMap* _blossom_set_index;

    98     BlossomSet* _blossom_set;

    99     NodeIntMap* _blossom_rep;

   101     IntNodeMap* _tree_set_index;

   102     TreeSet* _tree_set;

   104     NodeQueue _node_queue;

   105     int _process, _postpone, _last;

   107     int _node_num;

   109   private:

   111     void createStructures() {

   112       _node_num = countNodes(_graph);

   113       if (!_matching) {

   114         _matching = new MatchingMap(_graph);

   115       }

   116       if (!_status) {

   117         _status = new StatusMap(_graph);

   118       }

   119       if (!_ear) {

   120         _ear = new EarMap(_graph);

   121       }

   122       if (!_blossom_set) {

   123         _blossom_set_index = new IntNodeMap(_graph);

   124         _blossom_set = new BlossomSet(*_blossom_set_index);

   125       }

   126       if (!_blossom_rep) {

   127         _blossom_rep = new NodeIntMap(_node_num);

   128       }

   129       if (!_tree_set) {

   130         _tree_set_index = new IntNodeMap(_graph);

   131         _tree_set = new TreeSet(*_tree_set_index);

   132       }

   133       _node_queue.resize(_node_num);

   134     }

   136     void destroyStructures() {

   137       if (_matching) {

   138         delete _matching;

   139       }

   140       if (_status) {

   141         delete _status;

   142       }

   143       if (_ear) {

   144         delete _ear;

   145       }

   146       if (_blossom_set) {

   147         delete _blossom_set;

   148         delete _blossom_set_index;

   149       }

   150       if (_blossom_rep) {

   151         delete _blossom_rep;

   152       }

   153       if (_tree_set) {

   154         delete _tree_set_index;

   155         delete _tree_set;

   156       }

   157     }

   159     void processDense(const Node& n) {

   160       _process = _postpone = _last = 0;

   161       _node_queue[_last++] = n;

   163       while (_process != _last) {

   164         Node u = _node_queue[_process++];

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

   166           Node v = _graph.target(a);

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

   168             extendOnArc(a);

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

   170             augmentOnArc(a);

   171             return;

   172           }

   173         }

   174       }

   176       while (_postpone != _last) {

   177         Node u = _node_queue[_postpone++];

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

   180           Node v = _graph.target(a);

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

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

   184               shrinkOnEdge(a);

   185             }

   186           }

   188           while (_process != _last) {

   189             Node w = _node_queue[_process++];

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

   191               Node x = _graph.target(b);

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

   193                 extendOnArc(b);

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

   195                 augmentOnArc(b);

   196                 return;

   197               }

   198             }

   199           }

   200         }

   201       }

   202     }

   204     void processSparse(const Node& n) {

   205       _process = _last = 0;

   206       _node_queue[_last++] = n;

   207       while (_process != _last) {

   208         Node u = _node_queue[_process++];

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

   210           Node v = _graph.target(a);

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

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

   214               shrinkOnEdge(a);

   215             }

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

   217             extendOnArc(a);

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

   219             augmentOnArc(a);

   220             return;

   221           }

   222         }

   223       }

   224     }

   226     void shrinkOnEdge(const Edge& e) {

   227       Node nca = INVALID;

   229       {

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

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

   233         left_set.insert(left);

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

   236         right_set.insert(right);

   238         while (true) {

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

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

   241           left = (*_blossom_rep)[_blossom_set->

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

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

   244             nca = left;

   245             break;

   246           }

   247           left_set.insert(left);

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

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

   251           right = (*_blossom_rep)[_blossom_set->

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

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

   254             nca = right;

   255             break;

   256           }

   257           right_set.insert(right);

   258         }

   260         if (nca == INVALID) {

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

   262             nca = right;

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

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

   265               nca =(*_blossom_rep)[_blossom_set->

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

   267             }

   268           } else {

   269             nca = left;

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

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

   272               nca = (*_blossom_rep)[_blossom_set->

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

   274             }

   275           }

   276         }

   277       }

   279       {

   281         Node node = _graph.u(e);

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

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

   285         while (base != nca) {

   286           _ear->set(node, arc);

   288           Node n = node;

   289           while (n != base) {

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

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

   292             n = _graph.target(a);

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

   294           }

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

   296           _tree_set->erase(base);

   297           _tree_set->erase(node);

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

   299           _status->set(node, EVEN);

   300           _node_queue[_last++] = node;

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

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

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

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

   305         }

   306       }

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

   310       {

   312         Node node = _graph.v(e);

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

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

   316         while (base != nca) {

   317           _ear->set(node, arc);

   319           Node n = node;

   320           while (n != base) {

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

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

   323             n = _graph.target(a);

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

   325           }

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

   327           _tree_set->erase(base);

   328           _tree_set->erase(node);

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

   330           _status->set(node, EVEN);

   331           _node_queue[_last++] = node;

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

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

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

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

   336         }

   337       }

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

   340     }

   344     void extendOnArc(const Arc& a) {

   345       Node base = _graph.source(a);

   346       Node odd = _graph.target(a);

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

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

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

   351       _status->set(odd, ODD);

   352       _status->set(even, EVEN);

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

   354       _tree_set->insert(odd, tree);

   355       _tree_set->insert(even, tree);

   356       _node_queue[_last++] = even;

   358     }

   360     void augmentOnArc(const Arc& a) {

   361       Node even = _graph.source(a);

   362       Node odd = _graph.target(a);

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

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

   367       _status->set(odd, MATCHED);

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

   370       _matching->set(even, a);

   372       while (arc != INVALID) {

   373         odd = _graph.target(arc);

   374         arc = (*_ear)[odd];

   375         even = _graph.target(arc);

   376         _matching->set(odd, arc);

   377         arc = (*_matching)[even];

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

   379       }

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

   382            it != INVALID; ++it) {

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

   384           _status->set(it, MATCHED);

   385         } else {

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

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

   388                jt != INVALID; ++jt) {

   389             _status->set(jt, MATCHED);

   390           }

   391           _blossom_set->eraseClass(blossom);

   392         }

   393       }

   394       _tree_set->eraseClass(tree);

   396     }

   398   public:

   400     /// \brief Constructor

   401     ///

   402     /// Constructor.

   403     MaxMatching(const Graph& graph)

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

   405         _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),

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

   408     ~MaxMatching() {

   409       destroyStructures();

   410     }

   412     /// \name Execution control

   413     /// The simplest way to execute the algorithm is to use the member

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

   415     /// \n

   417     /// If you need more control on the execution, first you must call

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

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

   420     /// startParse() or startDense() functions.

   422     ///@{

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

   425     ///

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

   427     ///

   428     void init() {

   429       createStructures();

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

   431         _matching->set(n, INVALID);

   432         _status->set(n, UNMATCHED);

   433       }

   434     }

   436     ///\brief Finds a greedy matching for initial matching.

   437     ///

   438     ///For initial matchig it finds a maximal greedy matching.

   439     void greedyInit() {

   440       createStructures();

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

   442         _matching->set(n, INVALID);

   443         _status->set(n, UNMATCHED);

   444       }

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

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

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

   448             Node v = _graph.target(a);

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

   450               _matching->set(n, a);

   451               _status->set(n, MATCHED);

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

   453               _status->set(v, MATCHED);

   454               break;

   455             }

   456           }

   457         }

   458       }

   459     }

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

   463     ///

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

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

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

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

   468     template <typename MatchingMap>

   469     bool matchingInit(const MatchingMap& matching) {

   470       createStructures();

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

   473         _matching->set(n, INVALID);

   474         _status->set(n, UNMATCHED);

   475       }

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

   477         if (matching[e]) {

   479           Node u = _graph.u(e);

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

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

   482           _status->set(u, MATCHED);

   484           Node v = _graph.v(e);

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

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

   487           _status->set(v, MATCHED);

   488         }

   489       }

   490       return true;

   491     }

   493     /// \brief Starts Edmonds' algorithm

   494     ///

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

   496     void startSparse() {

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

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

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

   500           _tree_set->insert(n);

   501           _status->set(n, EVEN);

   502           processSparse(n);

   503         }

   504       }

   505     }

   507     /// \brief Starts Edmonds' algorithm.

   508     ///

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

   510     /// shrinks, giving a faster algorithm for dense graphs.

   511     void startDense() {

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

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

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

   515           _tree_set->insert(n);

   516           _status->set(n, EVEN);

   517           processDense(n);

   518         }

   519       }

   520     }

   523     /// \brief Runs Edmonds' algorithm

   524     ///

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

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

   527     /// postponing shrinks for dense graphs.

   528     void run() {

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

   530         greedyInit();

   531         startSparse();

   532       } else {

   533         init();

   534         startDense();

   535       }

   536     }

   538     /// @}

   540     /// \name Primal solution

   541     /// Functions for get the primal solution, ie. the matching.

   543     /// @{

   545     ///\brief Returns the size of the actual matching stored.

   546     ///

   547     ///Returns the size of the actual matching stored. After \ref

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

   549     int matchingSize() const {

   550       int size = 0;

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

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

   553           ++size;

   554         }

   555       }

   556       return size / 2;

   557     }

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

   560     ///

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

   562     bool matching(const Edge& edge) const {

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

   564     }

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

   567     ///

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

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

   570     Arc matching(const Node& n) const {

   571       return (*_matching)[n];

   572     }

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

   575     ///

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

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

   578     Node mate(const Node& n) const {

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

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

   581     }

   583     /// @}

   585     /// \name Dual solution

   586     /// Functions for get the dual solution, ie. the decomposition.

   588     /// @{

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

   591     /// decomposition.

   592     ///

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

   594     /// decomposition.

   595     Status decomposition(const Node& n) const {

   596       return (*_status)[n];

   597     }

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

   600     ///

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

   602     bool barrier(const Node& n) const {

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

   604     }

   606     /// @}

   608   };

   610   /// \ingroup matching

   611   ///

   612   /// \brief Weighted matching in general graphs

   613   ///

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

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

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

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

   618   ///

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

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

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

   622   /// following linear program.

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

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

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

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

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

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

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

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

   631   /// subsets of the nodes.

   632   ///

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

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

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

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

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

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

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

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

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

   642   ///

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

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

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

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

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

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

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

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

   651   template <typename _Graph,

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

   653   class MaxWeightedMatching {

   654   public:

   656     typedef _Graph Graph;

   657     typedef _WeightMap WeightMap;

   658     typedef typename WeightMap::Value Value;

   660     /// \brief Scaling factor for dual solution

   661     ///

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

   663     /// according to the value type.

   664     static const int dualScale =

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

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

   668     MatchingMap;

   670   private:

   672     TEMPLATE_GRAPH_TYPEDEFS(Graph);

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

   675     typedef std::vector<Node> BlossomNodeList;

   677     struct BlossomVariable {

   678       int begin, end;

   679       Value value;

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

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

   684     };

   686     typedef std::vector<BlossomVariable> BlossomPotential;

   688     const Graph& _graph;

   689     const WeightMap& _weight;

   691     MatchingMap* _matching;

   693     NodePotential* _node_potential;

   695     BlossomPotential _blossom_potential;

   696     BlossomNodeList _blossom_node_list;

   698     int _node_num;

   699     int _blossom_num;

   701     typedef RangeMap<int> IntIntMap;

   703     enum Status {

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

   705     };

   707     typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;

   708     struct BlossomData {

   709       int tree;

   710       Status status;

   711       Arc pred, next;

   712       Value pot, offset;

   713       Node base;

   714     };

   716     IntNodeMap *_blossom_index;

   717     BlossomSet *_blossom_set;

   718     RangeMap<BlossomData>* _blossom_data;

   720     IntNodeMap *_node_index;

   721     IntArcMap *_node_heap_index;

   723     struct NodeData {

   725       NodeData(IntArcMap& node_heap_index)

   726         : heap(node_heap_index) {}

   728       int blossom;

   729       Value pot;

   730       BinHeap<Value, IntArcMap> heap;

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

   733       int tree;

   734     };

   736     RangeMap<NodeData>* _node_data;

   738     typedef ExtendFindEnum<IntIntMap> TreeSet;

   740     IntIntMap *_tree_set_index;

   741     TreeSet *_tree_set;

   743     IntNodeMap *_delta1_index;

   744     BinHeap<Value, IntNodeMap> *_delta1;

   746     IntIntMap *_delta2_index;

   747     BinHeap<Value, IntIntMap> *_delta2;

   749     IntEdgeMap *_delta3_index;

   750     BinHeap<Value, IntEdgeMap> *_delta3;

   752     IntIntMap *_delta4_index;

   753     BinHeap<Value, IntIntMap> *_delta4;

   755     Value _delta_sum;

   757     void createStructures() {

   758       _node_num = countNodes(_graph);

   759       _blossom_num = _node_num * 3 / 2;

   761       if (!_matching) {

   762         _matching = new MatchingMap(_graph);

   763       }

   764       if (!_node_potential) {

   765         _node_potential = new NodePotential(_graph);

   766       }

   767       if (!_blossom_set) {

   768         _blossom_index = new IntNodeMap(_graph);

   769         _blossom_set = new BlossomSet(*_blossom_index);

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

   771       }

   773       if (!_node_index) {

   774         _node_index = new IntNodeMap(_graph);

   775         _node_heap_index = new IntArcMap(_graph);

   776         _node_data = new RangeMap<NodeData>(_node_num,

   777                                               NodeData(*_node_heap_index));

   778       }

   780       if (!_tree_set) {

   781         _tree_set_index = new IntIntMap(_blossom_num);

   782         _tree_set = new TreeSet(*_tree_set_index);

   783       }

   784       if (!_delta1) {

   785         _delta1_index = new IntNodeMap(_graph);

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

   787       }

   788       if (!_delta2) {

   789         _delta2_index = new IntIntMap(_blossom_num);

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

   791       }

   792       if (!_delta3) {

   793         _delta3_index = new IntEdgeMap(_graph);

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

   795       }

   796       if (!_delta4) {

   797         _delta4_index = new IntIntMap(_blossom_num);

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

   799       }

   800     }

   802     void destroyStructures() {

   803       _node_num = countNodes(_graph);

   804       _blossom_num = _node_num * 3 / 2;

   806       if (_matching) {

   807         delete _matching;

   808       }

   809       if (_node_potential) {

   810         delete _node_potential;

   811       }

   812       if (_blossom_set) {

   813         delete _blossom_index;

   814         delete _blossom_set;

   815         delete _blossom_data;

   816       }

   818       if (_node_index) {

   819         delete _node_index;

   820         delete _node_heap_index;

   821         delete _node_data;

   822       }

   824       if (_tree_set) {

   825         delete _tree_set_index;

   826         delete _tree_set;

   827       }

   828       if (_delta1) {

   829         delete _delta1_index;

   830         delete _delta1;

   831       }

   832       if (_delta2) {

   833         delete _delta2_index;

   834         delete _delta2;

   835       }

   836       if (_delta3) {

   837         delete _delta3_index;

   838         delete _delta3;

   839       }

   840       if (_delta4) {

   841         delete _delta4_index;

   842         delete _delta4;

   843       }

   844     }

   846     void matchedToEven(int blossom, int tree) {

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

   848         _delta2->erase(blossom);

   849       }

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

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

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

   854       }

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

   857            n != INVALID; ++n) {

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

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

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

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

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

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

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

   870           Node v = _graph.source(e);

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

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

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

   875             dualScale * _weight[e];

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

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

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

   880             }

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

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

   883               _delta3->push(e, rw);

   884             }

   885           } else {

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

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

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

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

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

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

   893                 it->second = e;

   894               }

   895             } else {

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

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

   898             }

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

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

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

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

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

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

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

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

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

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

   911                 }

   912               }

   913             }

   914           }

   915         }

   916       }

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

   918     }

   920     void matchedToOdd(int blossom) {

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

   922         _delta2->erase(blossom);

   923       }

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

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

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

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

   928       }

   929     }

   931     void evenToMatched(int blossom, int tree) {

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

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

   934       }

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

   937            n != INVALID; ++n) {

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

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

   941         _delta1->erase(n);

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

   944           Node v = _graph.source(e);

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

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

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

   949             dualScale * _weight[e];

   951           if (vb == blossom) {

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

   953               _delta3->erase(e);

   954             }

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

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

   958               _delta3->erase(e);

   959             }

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

   963             if (vt != tree) {

   965               Arc r = _graph.oppositeArc(e);

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

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

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

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

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

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

   974                   it->second = r;

   975                 }

   976               } else {

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

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

   979               }

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

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

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

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

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

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

   988                            _blossom_set->classPrio(blossom) -

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

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

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

   992                 }

   993               }

   994             }

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

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

   998               _delta3->erase(e);

   999             }

  1000           } else {

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

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

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

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

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

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

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

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

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

  1012               }

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

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

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

  1017                   _delta2->erase(vb);

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

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

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

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

  1022                 }

  1023               }

  1024             }

  1025           }

  1026         }

  1027       }

  1028     }

  1030     void oddToMatched(int blossom) {

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

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

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

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

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

  1037       }

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

  1040         _delta4->erase(blossom);

  1041       }

  1042     }

  1044     void oddToEven(int blossom, int tree) {

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

  1046         _delta4->erase(blossom);

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

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

  1049       }

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

  1052            n != INVALID; ++n) {

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

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

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

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

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

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

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

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

  1065           Node v = _graph.source(e);

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

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

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

  1070             dualScale * _weight[e];

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

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

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

  1075             }

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

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

  1078               _delta3->push(e, rw);

  1079             }

  1080           } else {

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

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

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

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

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

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

  1089                 it->second = e;

  1090               }

  1091             } else {

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

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

  1094             }

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

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

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

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

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

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

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

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

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

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

  1107                 }

  1108               }

  1109             }

  1110           }

  1111         }

  1112       }

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

  1114     }

  1117     void matchedToUnmatched(int blossom) {

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

  1119         _delta2->erase(blossom);

  1120       }

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

  1123            n != INVALID; ++n) {

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

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

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

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

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

  1132           Node v = _graph.target(e);

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

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

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

  1137             dualScale * _weight[e];

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

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

  1141               _delta3->push(e, rw);

  1142             }

  1143           }

  1144         }

  1145       }

  1146     }

  1148     void unmatchedToMatched(int blossom) {

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

  1150            n != INVALID; ++n) {

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

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

  1154           Node v = _graph.source(e);

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

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

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

  1159             dualScale * _weight[e];

  1161           if (vb == blossom) {

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

  1163               _delta3->erase(e);

  1164             }

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

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

  1168               _delta3->erase(e);

  1169             }

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

  1173             Arc r = _graph.oppositeArc(e);

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

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

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

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

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

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

  1182                 it->second = r;

  1183               }

  1184             } else {

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

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

  1187             }

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

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

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

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

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

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

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

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

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

  1199               }

  1200             }

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

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

  1204               _delta3->erase(e);

  1205             }

  1206           }

  1207         }

  1208       }

  1209     }

  1211     void alternatePath(int even, int tree) {

  1212       int odd;

  1214       evenToMatched(even, tree);

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

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

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

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

  1220         oddToMatched(odd);

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

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

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

  1225         evenToMatched(even, tree);

  1226         (*_blossom_data)[even].next =

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

  1228       }

  1230     }

  1232     void destroyTree(int tree) {

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

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

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

  1236           evenToMatched(b, tree);

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

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

  1239           oddToMatched(b);

  1240         }

  1241       }

  1242       _tree_set->eraseClass(tree);

  1243     }

  1246     void unmatchNode(const Node& node) {

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

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

  1250       alternatePath(blossom, tree);

  1251       destroyTree(tree);

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

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

  1255       matchedToUnmatched(blossom);

  1256     }

  1259     void augmentOnEdge(const Edge& edge) {

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

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

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

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

  1266         alternatePath(left, left_tree);

  1267         destroyTree(left_tree);

  1268       } else {

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

  1270         unmatchedToMatched(left);

  1271       }

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

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

  1275         alternatePath(right, right_tree);

  1276         destroyTree(right_tree);

  1277       } else {

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

  1279         unmatchedToMatched(right);

  1280       }

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

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

  1284     }

  1286     void extendOnArc(const Arc& arc) {

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

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

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

  1291       _tree_set->insert(odd, tree);

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

  1293       matchedToOdd(odd);

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

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

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

  1298       _tree_set->insert(even, tree);

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

  1300       matchedToEven(even, tree);

  1301     }

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

  1304       int nca = -1;

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

  1307       {

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

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

  1310         left_path.push_back(left);

  1311         left_set.insert(left);

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

  1314         right_path.push_back(right);

  1315         right_set.insert(right);

  1317         while (true) {

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

  1321           left =

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

  1323           left_path.push_back(left);

  1324           left =

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

  1326           left_path.push_back(left);

  1328           left_set.insert(left);

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

  1331             nca = left;

  1332             break;

  1333           }

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

  1337           right =

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

  1339           right_path.push_back(right);

  1340           right =

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

  1342           right_path.push_back(right);

  1344           right_set.insert(right);

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

  1347             nca = right;

  1348             break;

  1349           }

  1351         }

  1353         if (nca == -1) {

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

  1355             nca = right;

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

  1357               nca =

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

  1359               right_path.push_back(nca);

  1360               nca =

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

  1362               right_path.push_back(nca);

  1363             }

  1364           } else {

  1365             nca = left;

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

  1367               nca =

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

  1369               left_path.push_back(nca);

  1370               nca =

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

  1372               left_path.push_back(nca);

  1373             }

  1374           }

  1375         }

  1376       }

  1378       std::vector<int> subblossoms;

  1379       Arc prev;

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

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

  1383         subblossoms.push_back(left_path[i]);

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

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

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

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

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

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

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

  1392       }

  1394       int k = 0;

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

  1397       subblossoms.push_back(nca);

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

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

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

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

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

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

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

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

  1409         subblossoms.push_back(right_path[i]);

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

  1411       }

  1413       int surface =

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

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

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

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

  1419         }

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

  1421       }

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

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

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

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

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

  1429       _tree_set->insert(surface, tree);

  1430       _tree_set->erase(nca);

  1431     }

  1433     void splitBlossom(int blossom) {

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

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

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

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

  1440       oddToMatched(blossom);

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

  1442         _delta2->erase(blossom);

  1443       }

  1445       std::vector<int> subblossoms;

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

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

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

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

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

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

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

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

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

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

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

  1460         }

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

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

  1463           _delta2->push(subblossoms[i],

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

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

  1466         }

  1467       }

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

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

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

  1472           int sb = subblossoms[i];

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

  1474           (*_blossom_data)[sb].next =

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

  1476         }

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

  1479           int sb = subblossoms[i];

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

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

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

  1484           matchedToOdd(sb);

  1485           _tree_set->insert(sb, tree);

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

  1487           (*_blossom_data)[sb].next =

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

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

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

  1493           matchedToEven(tb, tree);

  1494           _tree_set->insert(tb, tree);

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

  1496         }

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

  1499         matchedToOdd(subblossoms[id]);

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

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

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

  1504       } else {

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

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

  1508           int sb = subblossoms[i];

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

  1510           (*_blossom_data)[sb].next =

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

  1512         }

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

  1515           int sb = subblossoms[i];

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

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

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

  1520           matchedToOdd(sb);

  1521           _tree_set->insert(sb, tree);

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

  1523           (*_blossom_data)[sb].pred =

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

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

  1527           matchedToEven(tb, tree);

  1528           _tree_set->insert(tb, tree);

  1529           (*_blossom_data)[tb].pred =

  1530             (*_blossom_data)[tb].next =

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

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

  1533         }

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

  1536         matchedToOdd(subblossoms[ib]);

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

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

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

  1540       }

  1541       _tree_set->erase(blossom);

  1542     }

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

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

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

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

  1549         _matching->set(base, matching);

  1550         _blossom_node_list.push_back(base);

  1551         _node_potential->set(base, pot);

  1552       } else {

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

  1555         int bn = _blossom_node_list.size();

  1557         std::vector<int> subblossoms;

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

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

  1560         int ib = -1;

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

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

  1563         }

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

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

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

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

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

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

  1572         }

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

  1575         int en = _blossom_node_list.size();

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

  1578       }

  1579     }

  1581     void extractMatching() {

  1582       std::vector<int> blossoms;

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

  1584         blossoms.push_back(c);

  1585       }

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

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

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

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

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

  1593                n != INVALID; ++n) {

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

  1595           }

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

  1598           Node base = _graph.source(matching);

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

  1600         } else {

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

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

  1603         }

  1604       }

  1605     }

  1607   public:

  1609     /// \brief Constructor

  1610     ///

  1611     /// Constructor.

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

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

  1614         _node_potential(0), _blossom_potential(), _blossom_node_list(),

  1615         _node_num(0), _blossom_num(0),

  1617         _blossom_index(0), _blossom_set(0), _blossom_data(0),

  1618         _node_index(0), _node_heap_index(0), _node_data(0),

  1619         _tree_set_index(0), _tree_set(0),

  1621         _delta1_index(0), _delta1(0),

  1622         _delta2_index(0), _delta2(0),

  1623         _delta3_index(0), _delta3(0),

  1624         _delta4_index(0), _delta4(0),

  1626         _delta_sum() {}

  1628     ~MaxWeightedMatching() {

  1629       destroyStructures();

  1630     }

  1632     /// \name Execution control

  1633     /// The simplest way to execute the algorithm is to use the member

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

  1636     ///@{

  1638     /// \brief Initialize the algorithm

  1639     ///

  1640     /// Initialize the algorithm

  1641     void init() {

  1642       createStructures();

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

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

  1646       }

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

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

  1649       }

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

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

  1652       }

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

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

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

  1656       }

  1658       int index = 0;

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

  1660         Value max = 0;

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

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

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

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

  1665           }

  1666         }

  1667         _node_index->set(n, index);

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

  1669         _delta1->push(n, max);

  1670         int blossom =

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

  1673         _tree_set->insert(blossom);

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

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

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

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

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

  1680         ++index;

  1681       }

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

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

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

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

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

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

  1688         }

  1689       }

  1690     }

  1692     /// \brief Starts the algorithm

  1693     ///

  1694     /// Starts the algorithm

  1695     void start() {

  1696       enum OpType {

  1697         D1, D2, D3, D4

  1698       };

  1700       int unmatched = _node_num;

  1701       while (unmatched > 0) {

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

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

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

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

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

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

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

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

  1714         _delta_sum = d1; OpType ot = D1;

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

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

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

  1720         switch (ot) {

  1721         case D1:

  1722           {

  1723             Node n = _delta1->top();

  1724             unmatchNode(n);

  1725             --unmatched;

  1726           }

  1727           break;

  1728         case D2:

  1729           {

  1730             int blossom = _delta2->top();

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

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

  1733             extendOnArc(e);

  1734           }

  1735           break;

  1736         case D3:

  1737           {

  1738             Edge e = _delta3->top();

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

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

  1743             if (left_blossom == right_blossom) {

  1744               _delta3->pop();

  1745             } else {

  1746               int left_tree;

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

  1748                 left_tree = _tree_set->find(left_blossom);

  1749               } else {

  1750                 left_tree = -1;

  1751                 ++unmatched;

  1752               }

  1753               int right_tree;

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

  1755                 right_tree = _tree_set->find(right_blossom);

  1756               } else {

  1757                 right_tree = -1;

  1758                 ++unmatched;

  1759               }

  1761               if (left_tree == right_tree) {

  1762                 shrinkOnEdge(e, left_tree);

  1763               } else {

  1764                 augmentOnEdge(e);

  1765                 unmatched -= 2;

  1766               }

  1767             }

  1768           } break;

  1769         case D4:

  1770           splitBlossom(_delta4->top());

  1771           break;

  1772         }

  1773       }

  1774       extractMatching();

  1775     }

  1777     /// \brief Runs %MaxWeightedMatching algorithm.

  1778     ///

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

  1780     ///

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

  1782     /// \code

  1783     ///   mwm.init();

  1784     ///   mwm.start();

  1785     /// \endcode

  1786     void run() {

  1787       init();

  1788       start();

  1789     }

  1791     /// @}

  1793     /// \name Primal solution

  1794     /// Functions for get the primal solution, ie. the matching.

  1796     /// @{

  1798     /// \brief Returns the matching value.

  1799     ///

  1800     /// Returns the matching value.

  1801     Value matchingValue() const {

  1802       Value sum = 0;

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

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

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

  1806         }

  1807       }

  1808       return sum /= 2;

  1809     }

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

  1812     ///

  1813     /// Returns the cardinality of the matching.

  1814     int matchingSize() const {

  1815       int num = 0;

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

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

  1818           ++num;

  1819         }

  1820       }

  1821       return num /= 2;

  1822     }

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

  1825     ///

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

  1827     bool matching(const Edge& edge) const {

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

  1829     }

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

  1832     ///

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

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

  1835     Arc matching(const Node& node) const {

  1836       return (*_matching)[node];

  1837     }

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

  1840     ///

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

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

  1843     Node mate(const Node& node) const {

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

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

  1846     }

  1848     /// @}

  1850     /// \name Dual solution

  1851     /// Functions for get the dual solution.

  1853     /// @{

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

  1856     ///

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

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

  1859     Value dualValue() const {

  1860       Value sum = 0;

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

  1862         sum += nodeValue(n);

  1863       }

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

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

  1866       }

  1867       return sum;

  1868     }

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

  1871     ///

  1872     /// Returns the the value of the node.

  1873     Value nodeValue(const Node& n) const {

  1874       return (*_node_potential)[n];

  1875     }

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

  1878     ///

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

  1880     /// \see BlossomIt

  1881     int blossomNum() const {

  1882       return _blossom_potential.size();

  1883     }

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

  1887     ///

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

  1889     int blossomSize(int k) const {

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

  1891     }

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

  1894     ///

  1895     /// Returns the the value of the blossom.

  1896     /// \see BlossomIt

  1897     Value blossomValue(int k) const {

  1898       return _blossom_potential[k].value;

  1899     }

  1901     /// \brief Lemon iterator for get the items of the blossom.

  1902     ///

  1903     /// Lemon iterator for get the nodes of the blossom. This class

  1904     /// provides a common style lemon iterator which gives back a

  1905     /// subset of the nodes.

  1906     class BlossomIt {

  1907     public:

  1909       /// \brief Constructor.

  1910       ///

  1911       /// Constructor for get the nodes of the variable.

  1912       BlossomIt(const MaxWeightedMatching& algorithm, int variable)

  1913         : _algorithm(&algorithm)

  1914       {

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

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

  1917       }

  1919       /// \brief Conversion to node.

  1920       ///

  1921       /// Conversion to node.

  1922       operator Node() const {

  1923         return _algorithm->_blossom_node_list[_index];

  1924       }

  1926       /// \brief Increment operator.

  1927       ///

  1928       /// Increment operator.

  1929       BlossomIt& operator++() {

  1930         ++_index;

  1931         return *this;

  1932       }

  1934       /// \brief Validity checking

  1935       ///

  1936       /// Checks whether the iterator is invalid.

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

  1939       /// \brief Validity checking

  1940       ///

  1941       /// Checks whether the iterator is valid.

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

  1944     private:

  1945       const MaxWeightedMatching* _algorithm;

  1946       int _last;

  1947       int _index;

  1948     };

  1950     /// @}

  1952   };

  1954   /// \ingroup matching

  1955   ///

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

  1957   ///

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

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

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

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

  1962   ///

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

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

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

  1966   /// formulated with the following linear program.

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

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

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

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

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

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

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

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

  1975   /// subsets of the nodes.

  1976   ///

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

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

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

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

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

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

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

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

  1985   ///

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

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

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

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

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

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

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

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

  1994   template <typename _Graph,

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

  1996   class MaxWeightedPerfectMatching {

  1997   public:

  1999     typedef _Graph Graph;

  2000     typedef _WeightMap WeightMap;

  2001     typedef typename WeightMap::Value Value;

  2003     /// \brief Scaling factor for dual solution

  2004     ///

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

  2006     /// according to the value type.

  2007     static const int dualScale =

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

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

  2011     MatchingMap;

  2013   private:

  2015     TEMPLATE_GRAPH_TYPEDEFS(Graph);

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

  2018     typedef std::vector<Node> BlossomNodeList;

  2020     struct BlossomVariable {

  2021       int begin, end;

  2022       Value value;

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

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

  2027     };

  2029     typedef std::vector<BlossomVariable> BlossomPotential;

  2031     const Graph& _graph;

  2032     const WeightMap& _weight;

  2034     MatchingMap* _matching;

  2036     NodePotential* _node_potential;

  2038     BlossomPotential _blossom_potential;

  2039     BlossomNodeList _blossom_node_list;

  2041     int _node_num;

  2042     int _blossom_num;

  2044     typedef RangeMap<int> IntIntMap;

  2046     enum Status {

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

  2048     };

  2050     typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;

  2051     struct BlossomData {

  2052       int tree;

  2053       Status status;

  2054       Arc pred, next;

  2055       Value pot, offset;

  2056     };

  2058     IntNodeMap *_blossom_index;

  2059     BlossomSet *_blossom_set;

  2060     RangeMap<BlossomData>* _blossom_data;

  2062     IntNodeMap *_node_index;

  2063     IntArcMap *_node_heap_index;

  2065     struct NodeData {

  2067       NodeData(IntArcMap& node_heap_index)

  2068         : heap(node_heap_index) {}

  2070       int blossom;

  2071       Value pot;

  2072       BinHeap<Value, IntArcMap> heap;

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

  2075       int tree;

  2076     };

  2078     RangeMap<NodeData>* _node_data;

  2080     typedef ExtendFindEnum<IntIntMap> TreeSet;

  2082     IntIntMap *_tree_set_index;

  2083     TreeSet *_tree_set;

  2085     IntIntMap *_delta2_index;

  2086     BinHeap<Value, IntIntMap> *_delta2;

  2088     IntEdgeMap *_delta3_index;

  2089     BinHeap<Value, IntEdgeMap> *_delta3;

  2091     IntIntMap *_delta4_index;

  2092     BinHeap<Value, IntIntMap> *_delta4;

  2094     Value _delta_sum;

  2096     void createStructures() {

  2097       _node_num = countNodes(_graph);

  2098       _blossom_num = _node_num * 3 / 2;

  2100       if (!_matching) {

  2101         _matching = new MatchingMap(_graph);

  2102       }

  2103       if (!_node_potential) {

  2104         _node_potential = new NodePotential(_graph);

  2105       }

  2106       if (!_blossom_set) {

  2107         _blossom_index = new IntNodeMap(_graph);

  2108         _blossom_set = new BlossomSet(*_blossom_index);

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

  2110       }

  2112       if (!_node_index) {

  2113         _node_index = new IntNodeMap(_graph);

  2114         _node_heap_index = new IntArcMap(_graph);

  2115         _node_data = new RangeMap<NodeData>(_node_num,

  2116                                             NodeData(*_node_heap_index));

  2117       }

  2119       if (!_tree_set) {

  2120         _tree_set_index = new IntIntMap(_blossom_num);

  2121         _tree_set = new TreeSet(*_tree_set_index);

  2122       }

  2123       if (!_delta2) {

  2124         _delta2_index = new IntIntMap(_blossom_num);

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

  2126       }

  2127       if (!_delta3) {

  2128         _delta3_index = new IntEdgeMap(_graph);

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

  2130       }

  2131       if (!_delta4) {

  2132         _delta4_index = new IntIntMap(_blossom_num);

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

  2134       }

  2135     }

  2137     void destroyStructures() {

  2138       _node_num = countNodes(_graph);

  2139       _blossom_num = _node_num * 3 / 2;

  2141       if (_matching) {

  2142         delete _matching;

  2143       }

  2144       if (_node_potential) {

  2145         delete _node_potential;

  2146       }

  2147       if (_blossom_set) {

  2148         delete _blossom_index;

  2149         delete _blossom_set;

  2150         delete _blossom_data;

  2151       }

  2153       if (_node_index) {

  2154         delete _node_index;

  2155         delete _node_heap_index;

  2156         delete _node_data;

  2157       }

  2159       if (_tree_set) {

  2160         delete _tree_set_index;

  2161         delete _tree_set;

  2162       }

  2163       if (_delta2) {

  2164         delete _delta2_index;

  2165         delete _delta2;

  2166       }

  2167       if (_delta3) {

  2168         delete _delta3_index;

  2169         delete _delta3;

  2170       }

  2171       if (_delta4) {

  2172         delete _delta4_index;

  2173         delete _delta4;

  2174       }

  2175     }

  2177     void matchedToEven(int blossom, int tree) {

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

  2179         _delta2->erase(blossom);

  2180       }

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

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

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

  2185       }

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

  2188            n != INVALID; ++n) {

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

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

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

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

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

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

  2199           Node v = _graph.source(e);

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

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

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

  2204             dualScale * _weight[e];

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

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

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

  2209             }

  2210           } else {

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

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

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

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

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

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

  2218                 it->second = e;

  2219               }

  2220             } else {

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

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

  2223             }

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

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

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

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

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

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

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

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

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

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

  2236                 }

  2237               }

  2238             }

  2239           }

  2240         }

  2241       }

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

  2243     }

  2245     void matchedToOdd(int blossom) {

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

  2247         _delta2->erase(blossom);

  2248       }

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

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

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

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

  2253       }

  2254     }

  2256     void evenToMatched(int blossom, int tree) {

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

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

  2259       }

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

  2262            n != INVALID; ++n) {

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

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

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

  2267           Node v = _graph.source(e);

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

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

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

  2272             dualScale * _weight[e];

  2274           if (vb == blossom) {

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

  2276               _delta3->erase(e);

  2277             }

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

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

  2281               _delta3->erase(e);

  2282             }

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

  2286             if (vt != tree) {

  2288               Arc r = _graph.oppositeArc(e);

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

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

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

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

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

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

  2297                   it->second = r;

  2298                 }

  2299               } else {

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

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

  2302               }

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

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

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

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

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

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

  2311                            _blossom_set->classPrio(blossom) -

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

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

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

  2315                 }

  2316               }

  2317             }

  2318           } else {

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

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

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

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

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

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

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

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

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

  2330               }

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

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

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

  2335                   _delta2->erase(vb);

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

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

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

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

  2340                 }

  2341               }

  2342             }

  2343           }

  2344         }

  2345       }

  2346     }

  2348     void oddToMatched(int blossom) {

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

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

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

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

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

  2355       }

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

  2358         _delta4->erase(blossom);

  2359       }

  2360     }

  2362     void oddToEven(int blossom, int tree) {

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

  2364         _delta4->erase(blossom);

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

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

  2367       }

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

  2370            n != INVALID; ++n) {

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

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

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

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

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

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

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

  2381           Node v = _graph.source(e);

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

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

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

  2386             dualScale * _weight[e];

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

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

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

  2391             }

  2392           } else {

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

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

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

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

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

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

  2401                 it->second = e;

  2402               }

  2403             } else {

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

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

  2406             }

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

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

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

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

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

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

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

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

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

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

  2419                 }

  2420               }

  2421             }

  2422           }

  2423         }

  2424       }

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

  2426     }

  2428     void alternatePath(int even, int tree) {

  2429       int odd;

  2431       evenToMatched(even, tree);

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

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

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

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

  2437         oddToMatched(odd);

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

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

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

  2442         evenToMatched(even, tree);

  2443         (*_blossom_data)[even].next =

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

  2445       }

  2447     }

  2449     void destroyTree(int tree) {

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

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

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

  2453           evenToMatched(b, tree);

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

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

  2456           oddToMatched(b);

  2457         }

  2458       }

  2459       _tree_set->eraseClass(tree);

  2460     }

  2462     void augmentOnEdge(const Edge& edge) {

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

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

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

  2468       alternatePath(left, left_tree);

  2469       destroyTree(left_tree);

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

  2472       alternatePath(right, right_tree);

  2473       destroyTree(right_tree);

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

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

  2477     }

  2479     void extendOnArc(const Arc& arc) {

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

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

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

  2484       _tree_set->insert(odd, tree);

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

  2486       matchedToOdd(odd);

  2487       (*_blossom_data)[odd].pred = arc;

  2489       int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));

  2490       (*_blossom_data)[even].pred = (*_blossom_data)[even].next;

  2491       _tree_set->insert(even, tree);

  2492       (*_blossom_data)[even].status = EVEN;

  2493       matchedToEven(even, tree);

  2494     }

  2496     void shrinkOnEdge(const Edge& edge, int tree) {

  2497       int nca = -1;

  2498       std::vector<int> left_path, right_path;

  2500       {

  2501         std::set<int> left_set, right_set;

  2502         int left = _blossom_set->find(_graph.u(edge));

  2503         left_path.push_back(left);

  2504         left_set.insert(left);

  2506         int right = _blossom_set->find(_graph.v(edge));

  2507         right_path.push_back(right);

  2508         right_set.insert(right);

  2510         while (true) {

  2512           if ((*_blossom_data)[left].pred == INVALID) break;

  2514           left =

  2515             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));

  2516           left_path.push_back(left);

  2517           left =

  2518             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));

  2519           left_path.push_back(left);

  2521           left_set.insert(left);

  2523           if (right_set.find(left) != right_set.end()) {

  2524             nca = left;

  2525             break;

  2526           }

  2528           if ((*_blossom_data)[right].pred == INVALID) break;

  2530           right =

  2531             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));

  2532           right_path.push_back(right);

  2533           right =

  2534             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));

  2535           right_path.push_back(right);

  2537           right_set.insert(right);

  2539           if (left_set.find(right) != left_set.end()) {

  2540             nca = right;

  2541             break;

  2542           }

  2544         }

  2546         if (nca == -1) {

  2547           if ((*_blossom_data)[left].pred == INVALID) {

  2548             nca = right;

  2549             while (left_set.find(nca) == left_set.end()) {

  2550               nca =

  2551                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));

  2552               right_path.push_back(nca);

  2553               nca =

  2554                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));

  2555               right_path.push_back(nca);

  2556             }

  2557           } else {

  2558             nca = left;

  2559             while (right_set.find(nca) == right_set.end()) {

  2560               nca =

  2561                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));

  2562               left_path.push_back(nca);

  2563               nca =

  2564                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));

  2565               left_path.push_back(nca);

  2566             }

  2567           }

  2568         }

  2569       }

  2571       std::vector<int> subblossoms;

  2572       Arc prev;

  2574       prev = _graph.direct(edge, true);

  2575       for (int i = 0; left_path[i] != nca; i += 2) {

  2576         subblossoms.push_back(left_path[i]);

  2577         (*_blossom_data)[left_path[i]].next = prev;

  2578         _tree_set->erase(left_path[i]);

  2580         subblossoms.push_back(left_path[i + 1]);

  2581         (*_blossom_data)[left_path[i + 1]].status = EVEN;

  2582         oddToEven(left_path[i + 1], tree);

  2583         _tree_set->erase(left_path[i + 1]);

  2584         prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);

  2585       }

  2587       int k = 0;

  2588       while (right_path[k] != nca) ++k;

  2590       subblossoms.push_back(nca);

  2591       (*_blossom_data)[nca].next = prev;

  2593       for (int i = k - 2; i >= 0; i -= 2) {

  2594         subblossoms.push_back(right_path[i + 1]);

  2595         (*_blossom_data)[right_path[i + 1]].status = EVEN;

  2596         oddToEven(right_path[i + 1], tree);

  2597         _tree_set->erase(right_path[i + 1]);

  2599         (*_blossom_data)[right_path[i + 1]].next =

  2600           (*_blossom_data)[right_path[i + 1]].pred;

  2602         subblossoms.push_back(right_path[i]);

  2603         _tree_set->erase(right_path[i]);

  2604       }

  2606       int surface =

  2607         _blossom_set->join(subblossoms.begin(), subblossoms.end());

  2609       for (int i = 0; i < int(subblossoms.size()); ++i) {

  2610         if (!_blossom_set->trivial(subblossoms[i])) {

  2611           (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;

  2612         }

  2613         (*_blossom_data)[subblossoms[i]].status = MATCHED;

  2614       }

  2616       (*_blossom_data)[surface].pot = -2 * _delta_sum;

  2617       (*_blossom_data)[surface].offset = 0;

  2618       (*_blossom_data)[surface].status = EVEN;

  2619       (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;

  2620       (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;

  2622       _tree_set->insert(surface, tree);

  2623       _tree_set->erase(nca);

  2624     }

  2626     void splitBlossom(int blossom) {

  2627       Arc next = (*_blossom_data)[blossom].next;

  2628       Arc pred = (*_blossom_data)[blossom].pred;

  2630       int tree = _tree_set->find(blossom);

  2632       (*_blossom_data)[blossom].status = MATCHED;

  2633       oddToMatched(blossom);

  2634       if (_delta2->state(blossom) == _delta2->IN_HEAP) {

  2635         _delta2->erase(blossom);

  2636       }

  2638       std::vector<int> subblossoms;

  2639       _blossom_set->split(blossom, std::back_inserter(subblossoms));

  2641       Value offset = (*_blossom_data)[blossom].offset;

  2642       int b = _blossom_set->find(_graph.source(pred));

  2643       int d = _blossom_set->find(_graph.source(next));

  2645       int ib = -1, id = -1;

  2646       for (int i = 0; i < int(subblossoms.size()); ++i) {

  2647         if (subblossoms[i] == b) ib = i;

  2648         if (subblossoms[i] == d) id = i;

  2650         (*_blossom_data)[subblossoms[i]].offset = offset;

  2651         if (!_blossom_set->trivial(subblossoms[i])) {

  2652           (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;

  2653         }

  2654         if (_blossom_set->classPrio(subblossoms[i]) !=

  2655             std::numeric_limits<Value>::max()) {

  2656           _delta2->push(subblossoms[i],

  2657                         _blossom_set->classPrio(subblossoms[i]) -

  2658                         (*_blossom_data)[subblossoms[i]].offset);

  2659         }

  2660       }

  2662       if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {

  2663         for (int i = (id + 1) % subblossoms.size();

  2664              i != ib; i = (i + 2) % subblossoms.size()) {

  2665           int sb = subblossoms[i];

  2666           int tb = subblossoms[(i + 1) % subblossoms.size()];

  2667           (*_blossom_data)[sb].next =

  2668             _graph.oppositeArc((*_blossom_data)[tb].next);

  2669         }

  2671         for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {

  2672           int sb = subblossoms[i];

  2673           int tb = subblossoms[(i + 1) % subblossoms.size()];

  2674           int ub = subblossoms[(i + 2) % subblossoms.size()];

  2676           (*_blossom_data)[sb].status = ODD;

  2677           matchedToOdd(sb);

  2678           _tree_set->insert(sb, tree);

  2679           (*_blossom_data)[sb].pred = pred;

  2680           (*_blossom_data)[sb].next =

  2681                            _graph.oppositeArc((*_blossom_data)[tb].next);

  2683           pred = (*_blossom_data)[ub].next;

  2685           (*_blossom_data)[tb].status = EVEN;

  2686           matchedToEven(tb, tree);

  2687           _tree_set->insert(tb, tree);

  2688           (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;

  2689         }

  2691         (*_blossom_data)[subblossoms[id]].status = ODD;

  2692         matchedToOdd(subblossoms[id]);

  2693         _tree_set->insert(subblossoms[id], tree);

  2694         (*_blossom_data)[subblossoms[id]].next = next;

  2695         (*_blossom_data)[subblossoms[id]].pred = pred;

  2697       } else {

  2699         for (int i = (ib + 1) % subblossoms.size();

  2700              i != id; i = (i + 2) % subblossoms.size()) {

  2701           int sb = subblossoms[i];

  2702           int tb = subblossoms[(i + 1) % subblossoms.size()];

  2703           (*_blossom_data)[sb].next =

  2704             _graph.oppositeArc((*_blossom_data)[tb].next);

  2705         }

  2707         for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {

  2708           int sb = subblossoms[i];

  2709           int tb = subblossoms[(i + 1) % subblossoms.size()];

  2710           int ub = subblossoms[(i + 2) % subblossoms.size()];

  2712           (*_blossom_data)[sb].status = ODD;

  2713           matchedToOdd(sb);

  2714           _tree_set->insert(sb, tree);

  2715           (*_blossom_data)[sb].next = next;

  2716           (*_blossom_data)[sb].pred =

  2717             _graph.oppositeArc((*_blossom_data)[tb].next);

  2719           (*_blossom_data)[tb].status = EVEN;

  2720           matchedToEven(tb, tree);

  2721           _tree_set->insert(tb, tree);

  2722           (*_blossom_data)[tb].pred =

  2723             (*_blossom_data)[tb].next =

  2724             _graph.oppositeArc((*_blossom_data)[ub].next);

  2725           next = (*_blossom_data)[ub].next;

  2726         }

  2728         (*_blossom_data)[subblossoms[ib]].status = ODD;

  2729         matchedToOdd(subblossoms[ib]);

  2730         _tree_set->insert(subblossoms[ib], tree);

  2731         (*_blossom_data)[subblossoms[ib]].next = next;

  2732         (*_blossom_data)[subblossoms[ib]].pred = pred;

  2733       }

  2734       _tree_set->erase(blossom);

  2735     }

  2737     void extractBlossom(int blossom, const Node& base, const Arc& matching) {

  2738       if (_blossom_set->trivial(blossom)) {

  2739         int bi = (*_node_index)[base];

  2740         Value pot = (*_node_data)[bi].pot;

  2742         _matching->set(base, matching);

  2743         _blossom_node_list.push_back(base);

  2744         _node_potential->set(base, pot);

  2745       } else {

  2747         Value pot = (*_blossom_data)[blossom].pot;

  2748         int bn = _blossom_node_list.size();

  2750         std::vector<int> subblossoms;

  2751         _blossom_set->split(blossom, std::back_inserter(subblossoms));

  2752         int b = _blossom_set->find(base);

  2753         int ib = -1;

  2754         for (int i = 0; i < int(subblossoms.size()); ++i) {

  2755           if (subblossoms[i] == b) { ib = i; break; }

  2756         }

  2758         for (int i = 1; i < int(subblossoms.size()); i += 2) {

  2759           int sb = subblossoms[(ib + i) % subblossoms.size()];

  2760           int tb = subblossoms[(ib + i + 1) % subblossoms.size()];

  2762           Arc m = (*_blossom_data)[tb].next;

  2763           extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));

  2764           extractBlossom(tb, _graph.source(m), m);

  2765         }

  2766         extractBlossom(subblossoms[ib], base, matching);

  2768         int en = _blossom_node_list.size();

  2770         _blossom_potential.push_back(BlossomVariable(bn, en, pot));

  2771       }

  2772     }

  2774     void extractMatching() {

  2775       std::vector<int> blossoms;

  2776       for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {

  2777         blossoms.push_back(c);

  2778       }

  2780       for (int i = 0; i < int(blossoms.size()); ++i) {

  2782         Value offset = (*_blossom_data)[blossoms[i]].offset;

  2783         (*_blossom_data)[blossoms[i]].pot += 2 * offset;

  2784         for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);

  2785              n != INVALID; ++n) {

  2786           (*_node_data)[(*_node_index)[n]].pot -= offset;

  2787         }

  2789         Arc matching = (*_blossom_data)[blossoms[i]].next;

  2790         Node base = _graph.source(matching);

  2791         extractBlossom(blossoms[i], base, matching);

  2792       }

  2793     }

  2795   public:

  2797     /// \brief Constructor

  2798     ///

  2799     /// Constructor.

  2800     MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)

  2801       : _graph(graph), _weight(weight), _matching(0),

  2802         _node_potential(0), _blossom_potential(), _blossom_node_list(),

  2803         _node_num(0), _blossom_num(0),

  2805         _blossom_index(0), _blossom_set(0), _blossom_data(0),

  2806         _node_index(0), _node_heap_index(0), _node_data(0),

  2807         _tree_set_index(0), _tree_set(0),

  2809         _delta2_index(0), _delta2(0),

  2810         _delta3_index(0), _delta3(0),

  2811         _delta4_index(0), _delta4(0),

  2813         _delta_sum() {}

  2815     ~MaxWeightedPerfectMatching() {

  2816       destroyStructures();

  2817     }

  2819     /// \name Execution control

  2820     /// The simplest way to execute the algorithm is to use the member

  2821     /// \c run() member function.

  2823     ///@{

  2825     /// \brief Initialize the algorithm

  2826     ///

  2827     /// Initialize the algorithm

  2828     void init() {

  2829       createStructures();

  2831       for (ArcIt e(_graph); e != INVALID; ++e) {

  2832         _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);

  2833       }

  2834       for (EdgeIt e(_graph); e != INVALID; ++e) {

  2835         _delta3_index->set(e, _delta3->PRE_HEAP);

  2836       }

  2837       for (int i = 0; i < _blossom_num; ++i) {

  2838         _delta2_index->set(i, _delta2->PRE_HEAP);

  2839         _delta4_index->set(i, _delta4->PRE_HEAP);

  2840       }

  2842       int index = 0;

  2843       for (NodeIt n(_graph); n != INVALID; ++n) {

  2844         Value max = - std::numeric_limits<Value>::max();

  2845         for (OutArcIt e(_graph, n); e != INVALID; ++e) {

  2846           if (_graph.target(e) == n) continue;

  2847           if ((dualScale * _weight[e]) / 2 > max) {

  2848             max = (dualScale * _weight[e]) / 2;

  2849           }

  2850         }

  2851         _node_index->set(n, index);

  2852         (*_node_data)[index].pot = max;

  2853         int blossom =

  2854           _blossom_set->insert(n, std::numeric_limits<Value>::max());

  2856         _tree_set->insert(blossom);

  2858         (*_blossom_data)[blossom].status = EVEN;

  2859         (*_blossom_data)[blossom].pred = INVALID;

  2860         (*_blossom_data)[blossom].next = INVALID;

  2861         (*_blossom_data)[blossom].pot = 0;

  2862         (*_blossom_data)[blossom].offset = 0;

  2863         ++index;

  2864       }

  2865       for (EdgeIt e(_graph); e != INVALID; ++e) {

  2866         int si = (*_node_index)[_graph.u(e)];

  2867         int ti = (*_node_index)[_graph.v(e)];

  2868         if (_graph.u(e) != _graph.v(e)) {

  2869           _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -

  2870                             dualScale * _weight[e]) / 2);

  2871         }

  2872       }

  2873     }

  2875     /// \brief Starts the algorithm

  2876     ///

  2877     /// Starts the algorithm

  2878     bool start() {

  2879       enum OpType {

  2880         D2, D3, D4

  2881       };

  2883       int unmatched = _node_num;

  2884       while (unmatched > 0) {

  2885         Value d2 = !_delta2->empty() ?

  2886           _delta2->prio() : std::numeric_limits<Value>::max();

  2888         Value d3 = !_delta3->empty() ?

  2889           _delta3->prio() : std::numeric_limits<Value>::max();

  2891         Value d4 = !_delta4->empty() ?

  2892           _delta4->prio() : std::numeric_limits<Value>::max();

  2894         _delta_sum = d2; OpType ot = D2;

  2895         if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }

  2896         if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }

  2898         if (_delta_sum == std::numeric_limits<Value>::max()) {

  2899           return false;

  2900         }

  2902         switch (ot) {

  2903         case D2:

  2904           {

  2905             int blossom = _delta2->top();

  2906             Node n = _blossom_set->classTop(blossom);

  2907             Arc e = (*_node_data)[(*_node_index)[n]].heap.top();

  2908             extendOnArc(e);

  2909           }

  2910           break;

  2911         case D3:

  2912           {

  2913             Edge e = _delta3->top();

  2915             int left_blossom = _blossom_set->find(_graph.u(e));

  2916             int right_blossom = _blossom_set->find(_graph.v(e));

  2918             if (left_blossom == right_blossom) {

  2919               _delta3->pop();

  2920             } else {

  2921               int left_tree = _tree_set->find(left_blossom);

  2922               int right_tree = _tree_set->find(right_blossom);

  2924               if (left_tree == right_tree) {

  2925                 shrinkOnEdge(e, left_tree);

  2926               } else {

  2927                 augmentOnEdge(e);

  2928                 unmatched -= 2;

  2929               }

  2930             }

  2931           } break;

  2932         case D4:

  2933           splitBlossom(_delta4->top());

  2934           break;

  2935         }

  2936       }

  2937       extractMatching();

  2938       return true;

  2939     }

  2941     /// \brief Runs %MaxWeightedPerfectMatching algorithm.

  2942     ///

  2943     /// This method runs the %MaxWeightedPerfectMatching algorithm.

  2944     ///

  2945     /// \note mwm.run() is just a shortcut of the following code.

  2946     /// \code

  2947     ///   mwm.init();

  2948     ///   mwm.start();

  2949     /// \endcode

  2950     bool run() {

  2951       init();

  2952       return start();

  2953     }

  2955     /// @}

  2957     /// \name Primal solution

  2958     /// Functions for get the primal solution, ie. the matching.

  2960     /// @{

  2962     /// \brief Returns the matching value.

  2963     ///

  2964     /// Returns the matching value.

  2965     Value matchingValue() const {

  2966       Value sum = 0;

  2967       for (NodeIt n(_graph); n != INVALID; ++n) {

  2968         if ((*_matching)[n] != INVALID) {

  2969           sum += _weight[(*_matching)[n]];

  2970         }

  2971       }

  2972       return sum /= 2;

  2973     }

  2975     /// \brief Returns true when the edge is in the matching.

  2976     ///

  2977     /// Returns true when the edge is in the matching.

  2978     bool matching(const Edge& edge) const {

  2979       return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;

  2980     }

  2982     /// \brief Returns the incident matching edge.

  2983     ///

  2984     /// Returns the incident matching arc from given edge.

  2985     Arc matching(const Node& node) const {

  2986       return (*_matching)[node];

  2987     }

  2989     /// \brief Returns the mate of the node.

  2990     ///

  2991     /// Returns the adjancent node in a mathcing arc.

  2992     Node mate(const Node& node) const {

  2993       return _graph.target((*_matching)[node]);

  2994     }

  2996     /// @}

  2998     /// \name Dual solution

  2999     /// Functions for get the dual solution.

  3001     /// @{

  3003     /// \brief Returns the value of the dual solution.

  3004     ///

  3005     /// Returns the value of the dual solution. It should be equal to

  3006     /// the primal value scaled by \ref dualScale "dual scale".

  3007     Value dualValue() const {

  3008       Value sum = 0;

  3009       for (NodeIt n(_graph); n != INVALID; ++n) {

  3010         sum += nodeValue(n);

  3011       }

  3012       for (int i = 0; i < blossomNum(); ++i) {

  3013         sum += blossomValue(i) * (blossomSize(i) / 2);

  3014       }

  3015       return sum;

  3016     }

  3018     /// \brief Returns the value of the node.

  3019     ///

  3020     /// Returns the the value of the node.

  3021     Value nodeValue(const Node& n) const {

  3022       return (*_node_potential)[n];

  3023     }

  3025     /// \brief Returns the number of the blossoms in the basis.

  3026     ///

  3027     /// Returns the number of the blossoms in the basis.

  3028     /// \see BlossomIt

  3029     int blossomNum() const {

  3030       return _blossom_potential.size();

  3031     }

  3034     /// \brief Returns the number of the nodes in the blossom.

  3035     ///

  3036     /// Returns the number of the nodes in the blossom.

  3037     int blossomSize(int k) const {

  3038       return _blossom_potential[k].end - _blossom_potential[k].begin;

  3039     }

  3041     /// \brief Returns the value of the blossom.

  3042     ///

  3043     /// Returns the the value of the blossom.

  3044     /// \see BlossomIt

  3045     Value blossomValue(int k) const {

  3046       return _blossom_potential[k].value;

  3047     }

  3049     /// \brief Lemon iterator for get the items of the blossom.

  3050     ///

  3051     /// Lemon iterator for get the nodes of the blossom. This class

  3052     /// provides a common style lemon iterator which gives back a

  3053     /// subset of the nodes.

  3054     class BlossomIt {

  3055     public:

  3057       /// \brief Constructor.

  3058       ///

  3059       /// Constructor for get the nodes of the variable.

  3060       BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)

  3061         : _algorithm(&algorithm)

  3062       {

  3063         _index = _algorithm->_blossom_potential[variable].begin;

  3064         _last = _algorithm->_blossom_potential[variable].end;

  3065       }

  3067       /// \brief Conversion to node.

  3068       ///

  3069       /// Conversion to node.

  3070       operator Node() const {

  3071         return _algorithm->_blossom_node_list[_index];

  3072       }

  3074       /// \brief Increment operator.

  3075       ///

  3076       /// Increment operator.

  3077       BlossomIt& operator++() {

  3078         ++_index;

  3079         return *this;

  3080       }

  3082       /// \brief Validity checking

  3083       ///

  3084       /// Checks whether the iterator is invalid.

  3085       bool operator==(Invalid) const { return _index == _last; }

  3087       /// \brief Validity checking

  3088       ///

  3089       /// Checks whether the iterator is valid.

  3090       bool operator!=(Invalid) const { return _index != _last; }

  3092     private:

  3093       const MaxWeightedPerfectMatching* _algorithm;

  3094       int _last;

  3095       int _index;

  3096     };

  3098     /// @}

  3100   };

  3103 } //END OF NAMESPACE LEMON

  3105 #endif //LEMON_MAX_MATCHING_H

author	Balazs Dezso <deba@inf.elte.hu>
	Mon, 13 Oct 2008 14:00:11 +0200
changeset 339	91d63b8b1a4c
parent 338	64ad48007fb2
child 342	5ba887b7def4
permissions	-rw-r--r--