lemon-main: lemon/max_matching.h@eda12d8ac953

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

     2  *

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

     4  *

     5  * Copyright (C) 2003-2009

     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

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

     8  *

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

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

    11  * precise terms see the accompanying LICENSE file.

    12  *

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

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

    15  * purpose.

    16  *

    17  */

    19 #ifndef LEMON_MAX_MATCHING_H

    20 #define LEMON_MAX_MATCHING_H

    22 #include <vector>

    23 #include <queue>

    24 #include <set>

    25 #include <limits>

    27 #include <lemon/core.h>

    28 #include <lemon/unionfind.h>

    29 #include <lemon/bin_heap.h>

    30 #include <lemon/maps.h>

    32 ///\ingroup matching

    33 ///\file

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

    36 namespace lemon {

    38   /// \ingroup matching

    39   ///

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

    41   ///

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

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

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

    45   ///

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

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

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

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

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

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

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

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

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

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

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

    57   ///

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

    59   template <typename GR>

    60   class MaxMatching {

    61   public:

    63     typedef GR Graph;

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

    65     MatchingMap;

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

    68     ///

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

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

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

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

    73     ///matching.

    74     enum Status {

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

    76     };

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

    80   private:

    82     TEMPLATE_GRAPH_TYPEDEFS(Graph);

    84     typedef UnionFindEnum<IntNodeMap> BlossomSet;

    85     typedef ExtendFindEnum<IntNodeMap> TreeSet;

    86     typedef RangeMap<Node> NodeIntMap;

    87     typedef MatchingMap EarMap;

    88     typedef std::vector<Node> NodeQueue;

    90     const Graph& _graph;

    91     MatchingMap* _matching;

    92     StatusMap* _status;

    94     EarMap* _ear;

    96     IntNodeMap* _blossom_set_index;

    97     BlossomSet* _blossom_set;

    98     NodeIntMap* _blossom_rep;

   100     IntNodeMap* _tree_set_index;

   101     TreeSet* _tree_set;

   103     NodeQueue _node_queue;

   104     int _process, _postpone, _last;

   106     int _node_num;

   108   private:

   110     void createStructures() {

   111       _node_num = countNodes(_graph);

   112       if (!_matching) {

   113         _matching = new MatchingMap(_graph);

   114       }

   115       if (!_status) {

   116         _status = new StatusMap(_graph);

   117       }

   118       if (!_ear) {

   119         _ear = new EarMap(_graph);

   120       }

   121       if (!_blossom_set) {

   122         _blossom_set_index = new IntNodeMap(_graph);

   123         _blossom_set = new BlossomSet(*_blossom_set_index);

   124       }

   125       if (!_blossom_rep) {

   126         _blossom_rep = new NodeIntMap(_node_num);

   127       }

   128       if (!_tree_set) {

   129         _tree_set_index = new IntNodeMap(_graph);

   130         _tree_set = new TreeSet(*_tree_set_index);

   131       }

   132       _node_queue.resize(_node_num);

   133     }

   135     void destroyStructures() {

   136       if (_matching) {

   137         delete _matching;

   138       }

   139       if (_status) {

   140         delete _status;

   141       }

   142       if (_ear) {

   143         delete _ear;

   144       }

   145       if (_blossom_set) {

   146         delete _blossom_set;

   147         delete _blossom_set_index;

   148       }

   149       if (_blossom_rep) {

   150         delete _blossom_rep;

   151       }

   152       if (_tree_set) {

   153         delete _tree_set_index;

   154         delete _tree_set;

   155       }

   156     }

   158     void processDense(const Node& n) {

   159       _process = _postpone = _last = 0;

   160       _node_queue[_last++] = n;

   162       while (_process != _last) {

   163         Node u = _node_queue[_process++];

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

   165           Node v = _graph.target(a);

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

   167             extendOnArc(a);

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

   169             augmentOnArc(a);

   170             return;

   171           }

   172         }

   173       }

   175       while (_postpone != _last) {

   176         Node u = _node_queue[_postpone++];

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

   179           Node v = _graph.target(a);

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

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

   183               shrinkOnEdge(a);

   184             }

   185           }

   187           while (_process != _last) {

   188             Node w = _node_queue[_process++];

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

   190               Node x = _graph.target(b);

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

   192                 extendOnArc(b);

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

   194                 augmentOnArc(b);

   195                 return;

   196               }

   197             }

   198           }

   199         }

   200       }

   201     }

   203     void processSparse(const Node& n) {

   204       _process = _last = 0;

   205       _node_queue[_last++] = n;

   206       while (_process != _last) {

   207         Node u = _node_queue[_process++];

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

   209           Node v = _graph.target(a);

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

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

   213               shrinkOnEdge(a);

   214             }

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

   216             extendOnArc(a);

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

   218             augmentOnArc(a);

   219             return;

   220           }

   221         }

   222       }

   223     }

   225     void shrinkOnEdge(const Edge& e) {

   226       Node nca = INVALID;

   228       {

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

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

   232         left_set.insert(left);

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

   235         right_set.insert(right);

   237         while (true) {

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

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

   240           left = (*_blossom_rep)[_blossom_set->

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

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

   243             nca = left;

   244             break;

   245           }

   246           left_set.insert(left);

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

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

   250           right = (*_blossom_rep)[_blossom_set->

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

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

   253             nca = right;

   254             break;

   255           }

   256           right_set.insert(right);

   257         }

   259         if (nca == INVALID) {

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

   261             nca = right;

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

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

   264               nca =(*_blossom_rep)[_blossom_set->

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

   266             }

   267           } else {

   268             nca = left;

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

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

   271               nca = (*_blossom_rep)[_blossom_set->

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

   273             }

   274           }

   275         }

   276       }

   278       {

   280         Node node = _graph.u(e);

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

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

   284         while (base != nca) {

   285           _ear->set(node, arc);

   287           Node n = node;

   288           while (n != base) {

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

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

   291             n = _graph.target(a);

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

   293           }

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

   295           _tree_set->erase(base);

   296           _tree_set->erase(node);

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

   298           _status->set(node, EVEN);

   299           _node_queue[_last++] = node;

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

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

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

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

   304         }

   305       }

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

   309       {

   311         Node node = _graph.v(e);

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

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

   315         while (base != nca) {

   316           _ear->set(node, arc);

   318           Node n = node;

   319           while (n != base) {

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

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

   322             n = _graph.target(a);

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

   324           }

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

   326           _tree_set->erase(base);

   327           _tree_set->erase(node);

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

   329           _status->set(node, EVEN);

   330           _node_queue[_last++] = node;

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

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

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

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

   335         }

   336       }

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

   339     }

   343     void extendOnArc(const Arc& a) {

   344       Node base = _graph.source(a);

   345       Node odd = _graph.target(a);

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

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

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

   350       _status->set(odd, ODD);

   351       _status->set(even, EVEN);

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

   353       _tree_set->insert(odd, tree);

   354       _tree_set->insert(even, tree);

   355       _node_queue[_last++] = even;

   357     }

   359     void augmentOnArc(const Arc& a) {

   360       Node even = _graph.source(a);

   361       Node odd = _graph.target(a);

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

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

   366       _status->set(odd, MATCHED);

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

   369       _matching->set(even, a);

   371       while (arc != INVALID) {

   372         odd = _graph.target(arc);

   373         arc = (*_ear)[odd];

   374         even = _graph.target(arc);

   375         _matching->set(odd, arc);

   376         arc = (*_matching)[even];

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

   378       }

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

   381            it != INVALID; ++it) {

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

   383           _status->set(it, MATCHED);

   384         } else {

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

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

   387                jt != INVALID; ++jt) {

   388             _status->set(jt, MATCHED);

   389           }

   390           _blossom_set->eraseClass(blossom);

   391         }

   392       }

   393       _tree_set->eraseClass(tree);

   395     }

   397   public:

   399     /// \brief Constructor

   400     ///

   401     /// Constructor.

   402     MaxMatching(const Graph& graph)

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

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

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

   407     ~MaxMatching() {

   408       destroyStructures();

   409     }

   411     /// \name Execution control

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

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

   414     /// \n

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

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

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

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

   421     ///@{

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

   424     ///

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

   426     ///

   427     void init() {

   428       createStructures();

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

   430         _matching->set(n, INVALID);

   431         _status->set(n, UNMATCHED);

   432       }

   433     }

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

   436     ///

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

   438     void greedyInit() {

   439       createStructures();

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

   441         _matching->set(n, INVALID);

   442         _status->set(n, UNMATCHED);

   443       }

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

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

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

   447             Node v = _graph.target(a);

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

   449               _matching->set(n, a);

   450               _status->set(n, MATCHED);

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

   452               _status->set(v, MATCHED);

   453               break;

   454             }

   455           }

   456         }

   457       }

   458     }

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

   462     ///

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

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

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

   466     /// \return \c true if the map contains a matching.

   467     template <typename MatchingMap>

   468     bool matchingInit(const MatchingMap& matching) {

   469       createStructures();

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

   472         _matching->set(n, INVALID);

   473         _status->set(n, UNMATCHED);

   474       }

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

   476         if (matching[e]) {

   478           Node u = _graph.u(e);

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

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

   481           _status->set(u, MATCHED);

   483           Node v = _graph.v(e);

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

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

   486           _status->set(v, MATCHED);

   487         }

   488       }

   489       return true;

   490     }

   492     /// \brief Starts Edmonds' algorithm

   493     ///

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

   495     void startSparse() {

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

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

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

   499           _tree_set->insert(n);

   500           _status->set(n, EVEN);

   501           processSparse(n);

   502         }

   503       }

   504     }

   506     /// \brief Starts Edmonds' algorithm.

   507     ///

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

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

   510     void startDense() {

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

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

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

   514           _tree_set->insert(n);

   515           _status->set(n, EVEN);

   516           processDense(n);

   517         }

   518       }

   519     }

   522     /// \brief Runs Edmonds' algorithm

   523     ///

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

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

   526     /// postponing shrinks for dense graphs.

   527     void run() {

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

   529         greedyInit();

   530         startSparse();

   531       } else {

   532         init();

   533         startDense();

   534       }

   535     }

   537     /// @}

   539     /// \name Primal solution

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

   542     /// @{

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

   545     ///

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

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

   548     int matchingSize() const {

   549       int size = 0;

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

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

   552           ++size;

   553         }

   554       }

   555       return size / 2;

   556     }

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

   559     ///

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

   561     bool matching(const Edge& edge) const {

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

   563     }

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

   566     ///

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

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

   569     Arc matching(const Node& n) const {

   570       return (*_matching)[n];

   571     }

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

   574     ///

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

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

   577     Node mate(const Node& n) const {

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

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

   580     }

   582     /// @}

   584     /// \name Dual solution

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

   587     /// @{

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

   590     /// decomposition.

   591     ///

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

   593     /// decomposition.

   594     Status decomposition(const Node& n) const {

   595       return (*_status)[n];

   596     }

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

   599     ///

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

   601     bool barrier(const Node& n) const {

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

   603     }

   605     /// @}

   607   };

   609   /// \ingroup matching

   610   ///

   611   /// \brief Weighted matching in general graphs

   612   ///

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

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

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

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

   617   ///

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

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

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

   621   /// following linear program.

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

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

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

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

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

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

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

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

   630   /// subsets of the nodes.

   631   ///

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

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

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

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

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

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

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

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

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

   641   ///

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

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

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

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

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

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

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

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

   650   template <typename GR,

   651             typename WM = typename GR::template EdgeMap<int> >

   652   class MaxWeightedMatching {

   653   public:

   655     ///\e

   656     typedef GR Graph;

   657     ///\e

   658     typedef WM WeightMap;

   659     ///\e

   660     typedef typename WeightMap::Value Value;

   662     /// \brief Scaling factor for dual solution

   663     ///

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

   665     /// according to the value type.

   666     static const int dualScale =

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

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

   670     MatchingMap;

   672   private:

   674     TEMPLATE_GRAPH_TYPEDEFS(Graph);

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

   677     typedef std::vector<Node> BlossomNodeList;

   679     struct BlossomVariable {

   680       int begin, end;

   681       Value value;

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

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

   686     };

   688     typedef std::vector<BlossomVariable> BlossomPotential;

   690     const Graph& _graph;

   691     const WeightMap& _weight;

   693     MatchingMap* _matching;

   695     NodePotential* _node_potential;

   697     BlossomPotential _blossom_potential;

   698     BlossomNodeList _blossom_node_list;

   700     int _node_num;

   701     int _blossom_num;

   703     typedef RangeMap<int> IntIntMap;

   705     enum Status {

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

   707     };

   709     typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;

   710     struct BlossomData {

   711       int tree;

   712       Status status;

   713       Arc pred, next;

   714       Value pot, offset;

   715       Node base;

   716     };

   718     IntNodeMap *_blossom_index;

   719     BlossomSet *_blossom_set;

   720     RangeMap<BlossomData>* _blossom_data;

   722     IntNodeMap *_node_index;

   723     IntArcMap *_node_heap_index;

   725     struct NodeData {

   727       NodeData(IntArcMap& node_heap_index)

   728         : heap(node_heap_index) {}

   730       int blossom;

   731       Value pot;

   732       BinHeap<Value, IntArcMap> heap;

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

   735       int tree;

   736     };

   738     RangeMap<NodeData>* _node_data;

   740     typedef ExtendFindEnum<IntIntMap> TreeSet;

   742     IntIntMap *_tree_set_index;

   743     TreeSet *_tree_set;

   745     IntNodeMap *_delta1_index;

   746     BinHeap<Value, IntNodeMap> *_delta1;

   748     IntIntMap *_delta2_index;

   749     BinHeap<Value, IntIntMap> *_delta2;

   751     IntEdgeMap *_delta3_index;

   752     BinHeap<Value, IntEdgeMap> *_delta3;

   754     IntIntMap *_delta4_index;

   755     BinHeap<Value, IntIntMap> *_delta4;

   757     Value _delta_sum;

   759     void createStructures() {

   760       _node_num = countNodes(_graph);

   761       _blossom_num = _node_num * 3 / 2;

   763       if (!_matching) {

   764         _matching = new MatchingMap(_graph);

   765       }

   766       if (!_node_potential) {

   767         _node_potential = new NodePotential(_graph);

   768       }

   769       if (!_blossom_set) {

   770         _blossom_index = new IntNodeMap(_graph);

   771         _blossom_set = new BlossomSet(*_blossom_index);

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

   773       }

   775       if (!_node_index) {

   776         _node_index = new IntNodeMap(_graph);

   777         _node_heap_index = new IntArcMap(_graph);

   778         _node_data = new RangeMap<NodeData>(_node_num,

   779                                               NodeData(*_node_heap_index));

   780       }

   782       if (!_tree_set) {

   783         _tree_set_index = new IntIntMap(_blossom_num);

   784         _tree_set = new TreeSet(*_tree_set_index);

   785       }

   786       if (!_delta1) {

   787         _delta1_index = new IntNodeMap(_graph);

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

   789       }

   790       if (!_delta2) {

   791         _delta2_index = new IntIntMap(_blossom_num);

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

   793       }

   794       if (!_delta3) {

   795         _delta3_index = new IntEdgeMap(_graph);

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

   797       }

   798       if (!_delta4) {

   799         _delta4_index = new IntIntMap(_blossom_num);

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

   801       }

   802     }

   804     void destroyStructures() {

   805       _node_num = countNodes(_graph);

   806       _blossom_num = _node_num * 3 / 2;

   808       if (_matching) {

   809         delete _matching;

   810       }

   811       if (_node_potential) {

   812         delete _node_potential;

   813       }

   814       if (_blossom_set) {

   815         delete _blossom_index;

   816         delete _blossom_set;

   817         delete _blossom_data;

   818       }

   820       if (_node_index) {

   821         delete _node_index;

   822         delete _node_heap_index;

   823         delete _node_data;

   824       }

   826       if (_tree_set) {

   827         delete _tree_set_index;

   828         delete _tree_set;

   829       }

   830       if (_delta1) {

   831         delete _delta1_index;

   832         delete _delta1;

   833       }

   834       if (_delta2) {

   835         delete _delta2_index;

   836         delete _delta2;

   837       }

   838       if (_delta3) {

   839         delete _delta3_index;

   840         delete _delta3;

   841       }

   842       if (_delta4) {

   843         delete _delta4_index;

   844         delete _delta4;

   845       }

   846     }

   848     void matchedToEven(int blossom, int tree) {

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

   850         _delta2->erase(blossom);

   851       }

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

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

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

   856       }

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

   859            n != INVALID; ++n) {

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

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

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

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

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

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

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

   872           Node v = _graph.source(e);

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

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

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

   877             dualScale * _weight[e];

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

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

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

   882             }

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

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

   885               _delta3->push(e, rw);

   886             }

   887           } else {

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

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

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

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

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

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

   895                 it->second = e;

   896               }

   897             } else {

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

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

   900             }

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

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

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

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

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

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

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

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

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

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

   913                 }

   914               }

   915             }

   916           }

   917         }

   918       }

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

   920     }

   922     void matchedToOdd(int blossom) {

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

   924         _delta2->erase(blossom);

   925       }

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

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

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

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

   930       }

   931     }

   933     void evenToMatched(int blossom, int tree) {

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

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

   936       }

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

   939            n != INVALID; ++n) {

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

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

   943         _delta1->erase(n);

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

   946           Node v = _graph.source(e);

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

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

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

   951             dualScale * _weight[e];

   953           if (vb == blossom) {

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

   955               _delta3->erase(e);

   956             }

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

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

   960               _delta3->erase(e);

   961             }

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

   965             if (vt != tree) {

   967               Arc r = _graph.oppositeArc(e);

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

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

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

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

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

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

   976                   it->second = r;

   977                 }

   978               } else {

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

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

   981               }

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

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

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

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

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

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

   990                            _blossom_set->classPrio(blossom) -

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

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

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

   994                 }

   995               }

   996             }

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

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

  1000               _delta3->erase(e);

  1001             }

  1002           } else {

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

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

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

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

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

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

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

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

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

  1014               }

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

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

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

  1019                   _delta2->erase(vb);

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

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

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

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

  1024                 }

  1025               }

  1026             }

  1027           }

  1028         }

  1029       }

  1030     }

  1032     void oddToMatched(int blossom) {

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

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

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

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

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

  1039       }

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

  1042         _delta4->erase(blossom);

  1043       }

  1044     }

  1046     void oddToEven(int blossom, int tree) {

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

  1048         _delta4->erase(blossom);

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

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

  1051       }

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

  1054            n != INVALID; ++n) {

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

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

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

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

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

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

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

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

  1067           Node v = _graph.source(e);

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

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

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

  1072             dualScale * _weight[e];

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

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

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

  1077             }

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

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

  1080               _delta3->push(e, rw);

  1081             }

  1082           } else {

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

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

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

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

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

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

  1091                 it->second = e;

  1092               }

  1093             } else {

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

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

  1096             }

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

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

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

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

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

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

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

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

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

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

  1109                 }

  1110               }

  1111             }

  1112           }

  1113         }

  1114       }

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

  1116     }

  1119     void matchedToUnmatched(int blossom) {

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

  1121         _delta2->erase(blossom);

  1122       }

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

  1125            n != INVALID; ++n) {

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

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

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

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

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

  1134           Node v = _graph.target(e);

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

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

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

  1139             dualScale * _weight[e];

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

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

  1143               _delta3->push(e, rw);

  1144             }

  1145           }

  1146         }

  1147       }

  1148     }

  1150     void unmatchedToMatched(int blossom) {

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

  1152            n != INVALID; ++n) {

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

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

  1156           Node v = _graph.source(e);

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

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

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

  1161             dualScale * _weight[e];

  1163           if (vb == blossom) {

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

  1165               _delta3->erase(e);

  1166             }

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

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

  1170               _delta3->erase(e);

  1171             }

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

  1175             Arc r = _graph.oppositeArc(e);

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

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

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

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

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

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

  1184                 it->second = r;

  1185               }

  1186             } else {

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

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

  1189             }

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

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

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

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

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

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

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

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

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

  1201               }

  1202             }

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

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

  1206               _delta3->erase(e);

  1207             }

  1208           }

  1209         }

  1210       }

  1211     }

  1213     void alternatePath(int even, int tree) {

  1214       int odd;

  1216       evenToMatched(even, tree);

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

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

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

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

  1222         oddToMatched(odd);

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

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

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

  1227         evenToMatched(even, tree);

  1228         (*_blossom_data)[even].next =

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

  1230       }

  1232     }

  1234     void destroyTree(int tree) {

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

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

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

  1238           evenToMatched(b, tree);

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

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

  1241           oddToMatched(b);

  1242         }

  1243       }

  1244       _tree_set->eraseClass(tree);

  1245     }

  1248     void unmatchNode(const Node& node) {

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

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

  1252       alternatePath(blossom, tree);

  1253       destroyTree(tree);

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

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

  1257       matchedToUnmatched(blossom);

  1258     }

  1261     void augmentOnEdge(const Edge& edge) {

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

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

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

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

  1268         alternatePath(left, left_tree);

  1269         destroyTree(left_tree);

  1270       } else {

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

  1272         unmatchedToMatched(left);

  1273       }

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

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

  1277         alternatePath(right, right_tree);

  1278         destroyTree(right_tree);

  1279       } else {

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

  1281         unmatchedToMatched(right);

  1282       }

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

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

  1286     }

  1288     void extendOnArc(const Arc& arc) {

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

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

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

  1293       _tree_set->insert(odd, tree);

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

  1295       matchedToOdd(odd);

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

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

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

  1300       _tree_set->insert(even, tree);

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

  1302       matchedToEven(even, tree);

  1303     }

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

  1306       int nca = -1;

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

  1309       {

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

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

  1312         left_path.push_back(left);

  1313         left_set.insert(left);

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

  1316         right_path.push_back(right);

  1317         right_set.insert(right);

  1319         while (true) {

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

  1323           left =

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

  1325           left_path.push_back(left);

  1326           left =

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

  1328           left_path.push_back(left);

  1330           left_set.insert(left);

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

  1333             nca = left;

  1334             break;

  1335           }

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

  1339           right =

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

  1341           right_path.push_back(right);

  1342           right =

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

  1344           right_path.push_back(right);

  1346           right_set.insert(right);

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

  1349             nca = right;

  1350             break;

  1351           }

  1353         }

  1355         if (nca == -1) {

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

  1357             nca = right;

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

  1359               nca =

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

  1361               right_path.push_back(nca);

  1362               nca =

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

  1364               right_path.push_back(nca);

  1365             }

  1366           } else {

  1367             nca = left;

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

  1369               nca =

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

  1371               left_path.push_back(nca);

  1372               nca =

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

  1374               left_path.push_back(nca);

  1375             }

  1376           }

  1377         }

  1378       }

  1380       std::vector<int> subblossoms;

  1381       Arc prev;

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

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

  1385         subblossoms.push_back(left_path[i]);

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

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

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

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

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

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

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

  1394       }

  1396       int k = 0;

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

  1399       subblossoms.push_back(nca);

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

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

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

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

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

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

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

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

  1411         subblossoms.push_back(right_path[i]);

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

  1413       }

  1415       int surface =

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

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

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

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

  1421         }

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

  1423       }

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

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

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

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

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

  1431       _tree_set->insert(surface, tree);

  1432       _tree_set->erase(nca);

  1433     }

  1435     void splitBlossom(int blossom) {

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

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

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

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

  1442       oddToMatched(blossom);

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

  1444         _delta2->erase(blossom);

  1445       }

  1447       std::vector<int> subblossoms;

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

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

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

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

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

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

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

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

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

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

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

  1462         }

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

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

  1465           _delta2->push(subblossoms[i],

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

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

  1468         }

  1469       }

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

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

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

  1474           int sb = subblossoms[i];

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

  1476           (*_blossom_data)[sb].next =

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

  1478         }

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

  1481           int sb = subblossoms[i];

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

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

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

  1486           matchedToOdd(sb);

  1487           _tree_set->insert(sb, tree);

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

  1489           (*_blossom_data)[sb].next =

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

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

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

  1495           matchedToEven(tb, tree);

  1496           _tree_set->insert(tb, tree);

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

  1498         }

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

  1501         matchedToOdd(subblossoms[id]);

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

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

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

  1506       } else {

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

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

  1510           int sb = subblossoms[i];

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

  1512           (*_blossom_data)[sb].next =

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

  1514         }

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

  1517           int sb = subblossoms[i];

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

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

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

  1522           matchedToOdd(sb);

  1523           _tree_set->insert(sb, tree);

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

  1525           (*_blossom_data)[sb].pred =

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

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

  1529           matchedToEven(tb, tree);

  1530           _tree_set->insert(tb, tree);

  1531           (*_blossom_data)[tb].pred =

  1532             (*_blossom_data)[tb].next =

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

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

  1535         }

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

  1538         matchedToOdd(subblossoms[ib]);

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

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

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

  1542       }

  1543       _tree_set->erase(blossom);

  1544     }

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

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

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

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

  1551         _matching->set(base, matching);

  1552         _blossom_node_list.push_back(base);

  1553         _node_potential->set(base, pot);

  1554       } else {

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

  1557         int bn = _blossom_node_list.size();

  1559         std::vector<int> subblossoms;

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

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

  1562         int ib = -1;

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

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

  1565         }

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

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

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

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

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

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

  1574         }

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

  1577         int en = _blossom_node_list.size();

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

  1580       }

  1581     }

  1583     void extractMatching() {

  1584       std::vector<int> blossoms;

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

  1586         blossoms.push_back(c);

  1587       }

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

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

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

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

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

  1595                n != INVALID; ++n) {

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

  1597           }

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

  1600           Node base = _graph.source(matching);

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

  1602         } else {

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

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

  1605         }

  1606       }

  1607     }

  1609   public:

  1611     /// \brief Constructor

  1612     ///

  1613     /// Constructor.

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

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

  1616         _node_potential(0), _blossom_potential(), _blossom_node_list(),

  1617         _node_num(0), _blossom_num(0),

  1619         _blossom_index(0), _blossom_set(0), _blossom_data(0),

  1620         _node_index(0), _node_heap_index(0), _node_data(0),

  1621         _tree_set_index(0), _tree_set(0),

  1623         _delta1_index(0), _delta1(0),

  1624         _delta2_index(0), _delta2(0),

  1625         _delta3_index(0), _delta3(0),

  1626         _delta4_index(0), _delta4(0),

  1628         _delta_sum() {}

  1630     ~MaxWeightedMatching() {

  1631       destroyStructures();

  1632     }

  1634     /// \name Execution control

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

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

  1638     ///@{

  1640     /// \brief Initialize the algorithm

  1641     ///

  1642     /// Initialize the algorithm

  1643     void init() {

  1644       createStructures();

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

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

  1648       }

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

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

  1651       }

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

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

  1654       }

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

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

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

  1658       }

  1660       int index = 0;

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

  1662         Value max = 0;

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

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

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

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

  1667           }

  1668         }

  1669         _node_index->set(n, index);

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

  1671         _delta1->push(n, max);

  1672         int blossom =

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

  1675         _tree_set->insert(blossom);

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

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

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

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

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

  1682         ++index;

  1683       }

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

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

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

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

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

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

  1690         }

  1691       }

  1692     }

  1694     /// \brief Starts the algorithm

  1695     ///

  1696     /// Starts the algorithm

  1697     void start() {

  1698       enum OpType {

  1699         D1, D2, D3, D4

  1700       };

  1702       int unmatched = _node_num;

  1703       while (unmatched > 0) {

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

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

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

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

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

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

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

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

  1716         _delta_sum = d1; OpType ot = D1;

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

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

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

  1722         switch (ot) {

  1723         case D1:

  1724           {

  1725             Node n = _delta1->top();

  1726             unmatchNode(n);

  1727             --unmatched;

  1728           }

  1729           break;

  1730         case D2:

  1731           {

  1732             int blossom = _delta2->top();

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

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

  1735             extendOnArc(e);

  1736           }

  1737           break;

  1738         case D3:

  1739           {

  1740             Edge e = _delta3->top();

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

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

  1745             if (left_blossom == right_blossom) {

  1746               _delta3->pop();

  1747             } else {

  1748               int left_tree;

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

  1750                 left_tree = _tree_set->find(left_blossom);

  1751               } else {

  1752                 left_tree = -1;

  1753                 ++unmatched;

  1754               }

  1755               int right_tree;

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

  1757                 right_tree = _tree_set->find(right_blossom);

  1758               } else {

  1759                 right_tree = -1;

  1760                 ++unmatched;

  1761               }

  1763               if (left_tree == right_tree) {

  1764                 shrinkOnEdge(e, left_tree);

  1765               } else {

  1766                 augmentOnEdge(e);

  1767                 unmatched -= 2;

  1768               }

  1769             }

  1770           } break;

  1771         case D4:

  1772           splitBlossom(_delta4->top());

  1773           break;

  1774         }

  1775       }

  1776       extractMatching();

  1777     }

  1779     /// \brief Runs %MaxWeightedMatching algorithm.

  1780     ///

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

  1782     ///

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

  1784     /// \code

  1785     ///   mwm.init();

  1786     ///   mwm.start();

  1787     /// \endcode

  1788     void run() {

  1789       init();

  1790       start();

  1791     }

  1793     /// @}

  1795     /// \name Primal solution

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

  1798     /// @{

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

  1801     ///

  1802     /// Returns the weight of the matching.

  1803     Value matchingValue() const {

  1804       Value sum = 0;

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

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

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

  1808         }

  1809       }

  1810       return sum /= 2;

  1811     }

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

  1814     ///

  1815     /// Returns the cardinality of the matching.

  1816     int matchingSize() const {

  1817       int num = 0;

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

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

  1820           ++num;

  1821         }

  1822       }

  1823       return num /= 2;

  1824     }

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

  1827     ///

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

  1829     bool matching(const Edge& edge) const {

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

  1831     }

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

  1834     ///

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

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

  1837     Arc matching(const Node& node) const {

  1838       return (*_matching)[node];

  1839     }

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

  1842     ///

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

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

  1845     Node mate(const Node& node) const {

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

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

  1848     }

  1850     /// @}

  1852     /// \name Dual solution

  1853     /// Functions to get the dual solution.

  1855     /// @{

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

  1858     ///

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

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

  1861     Value dualValue() const {

  1862       Value sum = 0;

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

  1864         sum += nodeValue(n);

  1865       }

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

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

  1868       }

  1869       return sum;

  1870     }

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

  1873     ///

  1874     /// Returns the the value of the node.

  1875     Value nodeValue(const Node& n) const {

  1876       return (*_node_potential)[n];

  1877     }

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

  1880     ///

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

  1882     /// \see BlossomIt

  1883     int blossomNum() const {

  1884       return _blossom_potential.size();

  1885     }

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

  1889     ///

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

  1891     int blossomSize(int k) const {

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

  1893     }

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

  1896     ///

  1897     /// Returns the the value of the blossom.

  1898     /// \see BlossomIt

  1899     Value blossomValue(int k) const {

  1900       return _blossom_potential[k].value;

  1901     }

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

  1904     ///

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

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

  1907     /// subset of the nodes.

  1908     class BlossomIt {

  1909     public:

  1911       /// \brief Constructor.

  1912       ///

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

  1914       BlossomIt(const MaxWeightedMatching& algorithm, int variable)

  1915         : _algorithm(&algorithm)

  1916       {

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

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

  1919       }

  1921       /// \brief Conversion to node.

  1922       ///

  1923       /// Conversion to node.

  1924       operator Node() const {

  1925         return _algorithm->_blossom_node_list[_index];

  1926       }

  1928       /// \brief Increment operator.

  1929       ///

  1930       /// Increment operator.

  1931       BlossomIt& operator++() {

  1932         ++_index;

  1933         return *this;

  1934       }

  1936       /// \brief Validity checking

  1937       ///

  1938       /// Checks whether the iterator is invalid.

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

  1941       /// \brief Validity checking

  1942       ///

  1943       /// Checks whether the iterator is valid.

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

  1946     private:

  1947       const MaxWeightedMatching* _algorithm;

  1948       int _last;

  1949       int _index;

  1950     };

  1952     /// @}

  1954   };

  1956   /// \ingroup matching

  1957   ///

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

  1959   ///

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

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

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

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

  1964   ///

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

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

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

  1968   /// formulated with the following linear program.

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

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

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

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

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

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

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

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

  1977   /// subsets of the nodes.

  1978   ///

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

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

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

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

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

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

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

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

  1987   ///

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

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

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

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

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

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

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

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

  1996   template <typename GR,

  1997             typename WM = typename GR::template EdgeMap<int> >

  1998   class MaxWeightedPerfectMatching {

  1999   public:

  2001     typedef GR Graph;

  2002     typedef WM WeightMap;

  2003     typedef typename WeightMap::Value Value;

  2005     /// \brief Scaling factor for dual solution

  2006     ///

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

  2008     /// according to the value type.

  2009     static const int dualScale =

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

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

  2013     MatchingMap;

  2015   private:

  2017     TEMPLATE_GRAPH_TYPEDEFS(Graph);

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

  2020     typedef std::vector<Node> BlossomNodeList;

  2022     struct BlossomVariable {

  2023       int begin, end;

  2024       Value value;

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

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

  2029     };

  2031     typedef std::vector<BlossomVariable> BlossomPotential;

  2033     const Graph& _graph;

  2034     const WeightMap& _weight;

  2036     MatchingMap* _matching;

  2038     NodePotential* _node_potential;

  2040     BlossomPotential _blossom_potential;

  2041     BlossomNodeList _blossom_node_list;

  2043     int _node_num;

  2044     int _blossom_num;

  2046     typedef RangeMap<int> IntIntMap;

  2048     enum Status {

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

  2050     };

  2052     typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;

  2053     struct BlossomData {

  2054       int tree;

  2055       Status status;

  2056       Arc pred, next;

  2057       Value pot, offset;

  2058     };

  2060     IntNodeMap *_blossom_index;

  2061     BlossomSet *_blossom_set;

  2062     RangeMap<BlossomData>* _blossom_data;

  2064     IntNodeMap *_node_index;

  2065     IntArcMap *_node_heap_index;

  2067     struct NodeData {

  2069       NodeData(IntArcMap& node_heap_index)

  2070         : heap(node_heap_index) {}

  2072       int blossom;

  2073       Value pot;

  2074       BinHeap<Value, IntArcMap> heap;

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

  2077       int tree;

  2078     };

  2080     RangeMap<NodeData>* _node_data;

  2082     typedef ExtendFindEnum<IntIntMap> TreeSet;

  2084     IntIntMap *_tree_set_index;

  2085     TreeSet *_tree_set;

  2087     IntIntMap *_delta2_index;

  2088     BinHeap<Value, IntIntMap> *_delta2;

  2090     IntEdgeMap *_delta3_index;

  2091     BinHeap<Value, IntEdgeMap> *_delta3;

  2093     IntIntMap *_delta4_index;

  2094     BinHeap<Value, IntIntMap> *_delta4;

  2096     Value _delta_sum;

  2098     void createStructures() {

  2099       _node_num = countNodes(_graph);

  2100       _blossom_num = _node_num * 3 / 2;

  2102       if (!_matching) {

  2103         _matching = new MatchingMap(_graph);

  2104       }

  2105       if (!_node_potential) {

  2106         _node_potential = new NodePotential(_graph);

  2107       }

  2108       if (!_blossom_set) {

  2109         _blossom_index = new IntNodeMap(_graph);

  2110         _blossom_set = new BlossomSet(*_blossom_index);

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

  2112       }

  2114       if (!_node_index) {

  2115         _node_index = new IntNodeMap(_graph);

  2116         _node_heap_index = new IntArcMap(_graph);

  2117         _node_data = new RangeMap<NodeData>(_node_num,

  2118                                             NodeData(*_node_heap_index));

  2119       }

  2121       if (!_tree_set) {

  2122         _tree_set_index = new IntIntMap(_blossom_num);

  2123         _tree_set = new TreeSet(*_tree_set_index);

  2124       }

  2125       if (!_delta2) {

  2126         _delta2_index = new IntIntMap(_blossom_num);

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

  2128       }

  2129       if (!_delta3) {

  2130         _delta3_index = new IntEdgeMap(_graph);

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

  2132       }

  2133       if (!_delta4) {

  2134         _delta4_index = new IntIntMap(_blossom_num);

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

  2136       }

  2137     }

  2139     void destroyStructures() {

  2140       _node_num = countNodes(_graph);

  2141       _blossom_num = _node_num * 3 / 2;

  2143       if (_matching) {

  2144         delete _matching;

  2145       }

  2146       if (_node_potential) {

  2147         delete _node_potential;

  2148       }

  2149       if (_blossom_set) {

  2150         delete _blossom_index;

  2151         delete _blossom_set;

  2152         delete _blossom_data;

  2153       }

  2155       if (_node_index) {

  2156         delete _node_index;

  2157         delete _node_heap_index;

  2158         delete _node_data;

  2159       }

  2161       if (_tree_set) {

  2162         delete _tree_set_index;

  2163         delete _tree_set;

  2164       }

  2165       if (_delta2) {

  2166         delete _delta2_index;

  2167         delete _delta2;

  2168       }

  2169       if (_delta3) {

  2170         delete _delta3_index;

  2171         delete _delta3;

  2172       }

  2173       if (_delta4) {

  2174         delete _delta4_index;

  2175         delete _delta4;

  2176       }

  2177     }

  2179     void matchedToEven(int blossom, int tree) {

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

  2181         _delta2->erase(blossom);

  2182       }

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

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

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

  2187       }

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

  2190            n != INVALID; ++n) {

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

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

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

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

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

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

  2201           Node v = _graph.source(e);

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

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

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

  2206             dualScale * _weight[e];

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

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

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

  2211             }

  2212           } else {

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

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

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

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

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

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

  2220                 it->second = e;

  2221               }

  2222             } else {

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

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

  2225             }

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

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

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

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

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

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

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

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

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

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

  2238                 }

  2239               }

  2240             }

  2241           }

  2242         }

  2243       }

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

  2245     }

  2247     void matchedToOdd(int blossom) {

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

  2249         _delta2->erase(blossom);

  2250       }

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

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

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

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

  2255       }

  2256     }

  2258     void evenToMatched(int blossom, int tree) {

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

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

  2261       }

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

  2264            n != INVALID; ++n) {

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

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

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

  2269           Node v = _graph.source(e);

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

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

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

  2274             dualScale * _weight[e];

  2276           if (vb == blossom) {

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

  2278               _delta3->erase(e);

  2279             }

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

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

  2283               _delta3->erase(e);

  2284             }

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

  2288             if (vt != tree) {

  2290               Arc r = _graph.oppositeArc(e);

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

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

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

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

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

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

  2299                   it->second = r;

  2300                 }

  2301               } else {

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

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

  2304               }

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

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

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

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

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

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

  2313                            _blossom_set->classPrio(blossom) -

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

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

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

  2317                 }

  2318               }

  2319             }

  2320           } else {

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

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

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

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

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

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

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

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

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

  2332               }

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

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

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

  2337                   _delta2->erase(vb);

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

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

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

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

  2342                 }

  2343               }

  2344             }

  2345           }

  2346         }

  2347       }

  2348     }

  2350     void oddToMatched(int blossom) {

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

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

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

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

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

  2357       }

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

  2360         _delta4->erase(blossom);

  2361       }

  2362     }

  2364     void oddToEven(int blossom, int tree) {

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

  2366         _delta4->erase(blossom);

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

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

  2369       }

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

  2372            n != INVALID; ++n) {

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

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

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

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

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

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

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

  2383           Node v = _graph.source(e);

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

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

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

  2388             dualScale * _weight[e];

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

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

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

  2393             }

  2394           } else {

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

  2397               (*_node_data)[vi].heap_index.find(tree);

  2399             if (it != (*_node_data)[vi].heap_index.end()) {

  2400               if ((*_node_data)[vi].heap[it->second] > rw) {

  2401                 (*_node_data)[vi].heap.replace(it->second, e);

  2402                 (*_node_data)[vi].heap.decrease(e, rw);

  2403                 it->second = e;

  2404               }

  2405             } else {

  2406               (*_node_data)[vi].heap.push(e, rw);

  2407               (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));

  2408             }

  2410             if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {

  2411               _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());

  2413               if ((*_blossom_data)[vb].status == MATCHED) {

  2414                 if (_delta2->state(vb) != _delta2->IN_HEAP) {

  2415                   _delta2->push(vb, _blossom_set->classPrio(vb) -

  2416                                (*_blossom_data)[vb].offset);

  2417                 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -

  2418                            (*_blossom_data)[vb].offset) {

  2419                   _delta2->decrease(vb, _blossom_set->classPrio(vb) -

  2420                                    (*_blossom_data)[vb].offset);

  2421                 }

  2422               }

  2423             }

  2424           }

  2425         }

  2426       }

  2427       (*_blossom_data)[blossom].offset = 0;

  2428     }

  2430     void alternatePath(int even, int tree) {

  2431       int odd;

  2433       evenToMatched(even, tree);

  2434       (*_blossom_data)[even].status = MATCHED;

  2436       while ((*_blossom_data)[even].pred != INVALID) {

  2437         odd = _blossom_set->find(_graph.target((*_blossom_data)[even].pred));

  2438         (*_blossom_data)[odd].status = MATCHED;

  2439         oddToMatched(odd);

  2440         (*_blossom_data)[odd].next = (*_blossom_data)[odd].pred;

  2442         even = _blossom_set->find(_graph.target((*_blossom_data)[odd].pred));

  2443         (*_blossom_data)[even].status = MATCHED;

  2444         evenToMatched(even, tree);

  2445         (*_blossom_data)[even].next =

  2446           _graph.oppositeArc((*_blossom_data)[odd].pred);

  2447       }

  2449     }

  2451     void destroyTree(int tree) {

  2452       for (TreeSet::ItemIt b(*_tree_set, tree); b != INVALID; ++b) {

  2453         if ((*_blossom_data)[b].status == EVEN) {

  2454           (*_blossom_data)[b].status = MATCHED;

  2455           evenToMatched(b, tree);

  2456         } else if ((*_blossom_data)[b].status == ODD) {

  2457           (*_blossom_data)[b].status = MATCHED;

  2458           oddToMatched(b);

  2459         }

  2460       }

  2461       _tree_set->eraseClass(tree);

  2462     }

  2464     void augmentOnEdge(const Edge& edge) {

  2466       int left = _blossom_set->find(_graph.u(edge));

  2467       int right = _blossom_set->find(_graph.v(edge));

  2469       int left_tree = _tree_set->find(left);

  2470       alternatePath(left, left_tree);

  2471       destroyTree(left_tree);

  2473       int right_tree = _tree_set->find(right);

  2474       alternatePath(right, right_tree);

  2475       destroyTree(right_tree);

  2477       (*_blossom_data)[left].next = _graph.direct(edge, true);

  2478       (*_blossom_data)[right].next = _graph.direct(edge, false);

  2479     }

  2481     void extendOnArc(const Arc& arc) {

  2482       int base = _blossom_set->find(_graph.target(arc));

  2483       int tree = _tree_set->find(base);

  2485       int odd = _blossom_set->find(_graph.source(arc));

  2486       _tree_set->insert(odd, tree);

  2487       (*_blossom_data)[odd].status = ODD;

  2488       matchedToOdd(odd);

  2489       (*_blossom_data)[odd].pred = arc;

  2491       int even = _blossom_set->find(_graph.target((*_blossom_data)[odd].next));

  2492       (*_blossom_data)[even].pred = (*_blossom_data)[even].next;

  2493       _tree_set->insert(even, tree);

  2494       (*_blossom_data)[even].status = EVEN;

  2495       matchedToEven(even, tree);

  2496     }

  2498     void shrinkOnEdge(const Edge& edge, int tree) {

  2499       int nca = -1;

  2500       std::vector<int> left_path, right_path;

  2502       {

  2503         std::set<int> left_set, right_set;

  2504         int left = _blossom_set->find(_graph.u(edge));

  2505         left_path.push_back(left);

  2506         left_set.insert(left);

  2508         int right = _blossom_set->find(_graph.v(edge));

  2509         right_path.push_back(right);

  2510         right_set.insert(right);

  2512         while (true) {

  2514           if ((*_blossom_data)[left].pred == INVALID) break;

  2516           left =

  2517             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));

  2518           left_path.push_back(left);

  2519           left =

  2520             _blossom_set->find(_graph.target((*_blossom_data)[left].pred));

  2521           left_path.push_back(left);

  2523           left_set.insert(left);

  2525           if (right_set.find(left) != right_set.end()) {

  2526             nca = left;

  2527             break;

  2528           }

  2530           if ((*_blossom_data)[right].pred == INVALID) break;

  2532           right =

  2533             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));

  2534           right_path.push_back(right);

  2535           right =

  2536             _blossom_set->find(_graph.target((*_blossom_data)[right].pred));

  2537           right_path.push_back(right);

  2539           right_set.insert(right);

  2541           if (left_set.find(right) != left_set.end()) {

  2542             nca = right;

  2543             break;

  2544           }

  2546         }

  2548         if (nca == -1) {

  2549           if ((*_blossom_data)[left].pred == INVALID) {

  2550             nca = right;

  2551             while (left_set.find(nca) == left_set.end()) {

  2552               nca =

  2553                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));

  2554               right_path.push_back(nca);

  2555               nca =

  2556                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));

  2557               right_path.push_back(nca);

  2558             }

  2559           } else {

  2560             nca = left;

  2561             while (right_set.find(nca) == right_set.end()) {

  2562               nca =

  2563                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));

  2564               left_path.push_back(nca);

  2565               nca =

  2566                 _blossom_set->find(_graph.target((*_blossom_data)[nca].pred));

  2567               left_path.push_back(nca);

  2568             }

  2569           }

  2570         }

  2571       }

  2573       std::vector<int> subblossoms;

  2574       Arc prev;

  2576       prev = _graph.direct(edge, true);

  2577       for (int i = 0; left_path[i] != nca; i += 2) {

  2578         subblossoms.push_back(left_path[i]);

  2579         (*_blossom_data)[left_path[i]].next = prev;

  2580         _tree_set->erase(left_path[i]);

  2582         subblossoms.push_back(left_path[i + 1]);

  2583         (*_blossom_data)[left_path[i + 1]].status = EVEN;

  2584         oddToEven(left_path[i + 1], tree);

  2585         _tree_set->erase(left_path[i + 1]);

  2586         prev = _graph.oppositeArc((*_blossom_data)[left_path[i + 1]].pred);

  2587       }

  2589       int k = 0;

  2590       while (right_path[k] != nca) ++k;

  2592       subblossoms.push_back(nca);

  2593       (*_blossom_data)[nca].next = prev;

  2595       for (int i = k - 2; i >= 0; i -= 2) {

  2596         subblossoms.push_back(right_path[i + 1]);

  2597         (*_blossom_data)[right_path[i + 1]].status = EVEN;

  2598         oddToEven(right_path[i + 1], tree);

  2599         _tree_set->erase(right_path[i + 1]);

  2601         (*_blossom_data)[right_path[i + 1]].next =

  2602           (*_blossom_data)[right_path[i + 1]].pred;

  2604         subblossoms.push_back(right_path[i]);

  2605         _tree_set->erase(right_path[i]);

  2606       }

  2608       int surface =

  2609         _blossom_set->join(subblossoms.begin(), subblossoms.end());

  2611       for (int i = 0; i < int(subblossoms.size()); ++i) {

  2612         if (!_blossom_set->trivial(subblossoms[i])) {

  2613           (*_blossom_data)[subblossoms[i]].pot += 2 * _delta_sum;

  2614         }

  2615         (*_blossom_data)[subblossoms[i]].status = MATCHED;

  2616       }

  2618       (*_blossom_data)[surface].pot = -2 * _delta_sum;

  2619       (*_blossom_data)[surface].offset = 0;

  2620       (*_blossom_data)[surface].status = EVEN;

  2621       (*_blossom_data)[surface].pred = (*_blossom_data)[nca].pred;

  2622       (*_blossom_data)[surface].next = (*_blossom_data)[nca].pred;

  2624       _tree_set->insert(surface, tree);

  2625       _tree_set->erase(nca);

  2626     }

  2628     void splitBlossom(int blossom) {

  2629       Arc next = (*_blossom_data)[blossom].next;

  2630       Arc pred = (*_blossom_data)[blossom].pred;

  2632       int tree = _tree_set->find(blossom);

  2634       (*_blossom_data)[blossom].status = MATCHED;

  2635       oddToMatched(blossom);

  2636       if (_delta2->state(blossom) == _delta2->IN_HEAP) {

  2637         _delta2->erase(blossom);

  2638       }

  2640       std::vector<int> subblossoms;

  2641       _blossom_set->split(blossom, std::back_inserter(subblossoms));

  2643       Value offset = (*_blossom_data)[blossom].offset;

  2644       int b = _blossom_set->find(_graph.source(pred));

  2645       int d = _blossom_set->find(_graph.source(next));

  2647       int ib = -1, id = -1;

  2648       for (int i = 0; i < int(subblossoms.size()); ++i) {

  2649         if (subblossoms[i] == b) ib = i;

  2650         if (subblossoms[i] == d) id = i;

  2652         (*_blossom_data)[subblossoms[i]].offset = offset;

  2653         if (!_blossom_set->trivial(subblossoms[i])) {

  2654           (*_blossom_data)[subblossoms[i]].pot -= 2 * offset;

  2655         }

  2656         if (_blossom_set->classPrio(subblossoms[i]) !=

  2657             std::numeric_limits<Value>::max()) {

  2658           _delta2->push(subblossoms[i],

  2659                         _blossom_set->classPrio(subblossoms[i]) -

  2660                         (*_blossom_data)[subblossoms[i]].offset);

  2661         }

  2662       }

  2664       if (id > ib ? ((id - ib) % 2 == 0) : ((ib - id) % 2 == 1)) {

  2665         for (int i = (id + 1) % subblossoms.size();

  2666              i != ib; i = (i + 2) % subblossoms.size()) {

  2667           int sb = subblossoms[i];

  2668           int tb = subblossoms[(i + 1) % subblossoms.size()];

  2669           (*_blossom_data)[sb].next =

  2670             _graph.oppositeArc((*_blossom_data)[tb].next);

  2671         }

  2673         for (int i = ib; i != id; i = (i + 2) % subblossoms.size()) {

  2674           int sb = subblossoms[i];

  2675           int tb = subblossoms[(i + 1) % subblossoms.size()];

  2676           int ub = subblossoms[(i + 2) % subblossoms.size()];

  2678           (*_blossom_data)[sb].status = ODD;

  2679           matchedToOdd(sb);

  2680           _tree_set->insert(sb, tree);

  2681           (*_blossom_data)[sb].pred = pred;

  2682           (*_blossom_data)[sb].next =

  2683                            _graph.oppositeArc((*_blossom_data)[tb].next);

  2685           pred = (*_blossom_data)[ub].next;

  2687           (*_blossom_data)[tb].status = EVEN;

  2688           matchedToEven(tb, tree);

  2689           _tree_set->insert(tb, tree);

  2690           (*_blossom_data)[tb].pred = (*_blossom_data)[tb].next;

  2691         }

  2693         (*_blossom_data)[subblossoms[id]].status = ODD;

  2694         matchedToOdd(subblossoms[id]);

  2695         _tree_set->insert(subblossoms[id], tree);

  2696         (*_blossom_data)[subblossoms[id]].next = next;

  2697         (*_blossom_data)[subblossoms[id]].pred = pred;

  2699       } else {

  2701         for (int i = (ib + 1) % subblossoms.size();

  2702              i != id; i = (i + 2) % subblossoms.size()) {

  2703           int sb = subblossoms[i];

  2704           int tb = subblossoms[(i + 1) % subblossoms.size()];

  2705           (*_blossom_data)[sb].next =

  2706             _graph.oppositeArc((*_blossom_data)[tb].next);

  2707         }

  2709         for (int i = id; i != ib; i = (i + 2) % subblossoms.size()) {

  2710           int sb = subblossoms[i];

  2711           int tb = subblossoms[(i + 1) % subblossoms.size()];

  2712           int ub = subblossoms[(i + 2) % subblossoms.size()];

  2714           (*_blossom_data)[sb].status = ODD;

  2715           matchedToOdd(sb);

  2716           _tree_set->insert(sb, tree);

  2717           (*_blossom_data)[sb].next = next;

  2718           (*_blossom_data)[sb].pred =

  2719             _graph.oppositeArc((*_blossom_data)[tb].next);

  2721           (*_blossom_data)[tb].status = EVEN;

  2722           matchedToEven(tb, tree);

  2723           _tree_set->insert(tb, tree);

  2724           (*_blossom_data)[tb].pred =

  2725             (*_blossom_data)[tb].next =

  2726             _graph.oppositeArc((*_blossom_data)[ub].next);

  2727           next = (*_blossom_data)[ub].next;

  2728         }

  2730         (*_blossom_data)[subblossoms[ib]].status = ODD;

  2731         matchedToOdd(subblossoms[ib]);

  2732         _tree_set->insert(subblossoms[ib], tree);

  2733         (*_blossom_data)[subblossoms[ib]].next = next;

  2734         (*_blossom_data)[subblossoms[ib]].pred = pred;

  2735       }

  2736       _tree_set->erase(blossom);

  2737     }

  2739     void extractBlossom(int blossom, const Node& base, const Arc& matching) {

  2740       if (_blossom_set->trivial(blossom)) {

  2741         int bi = (*_node_index)[base];

  2742         Value pot = (*_node_data)[bi].pot;

  2744         _matching->set(base, matching);

  2745         _blossom_node_list.push_back(base);

  2746         _node_potential->set(base, pot);

  2747       } else {

  2749         Value pot = (*_blossom_data)[blossom].pot;

  2750         int bn = _blossom_node_list.size();

  2752         std::vector<int> subblossoms;

  2753         _blossom_set->split(blossom, std::back_inserter(subblossoms));

  2754         int b = _blossom_set->find(base);

  2755         int ib = -1;

  2756         for (int i = 0; i < int(subblossoms.size()); ++i) {

  2757           if (subblossoms[i] == b) { ib = i; break; }

  2758         }

  2760         for (int i = 1; i < int(subblossoms.size()); i += 2) {

  2761           int sb = subblossoms[(ib + i) % subblossoms.size()];

  2762           int tb = subblossoms[(ib + i + 1) % subblossoms.size()];

  2764           Arc m = (*_blossom_data)[tb].next;

  2765           extractBlossom(sb, _graph.target(m), _graph.oppositeArc(m));

  2766           extractBlossom(tb, _graph.source(m), m);

  2767         }

  2768         extractBlossom(subblossoms[ib], base, matching);

  2770         int en = _blossom_node_list.size();

  2772         _blossom_potential.push_back(BlossomVariable(bn, en, pot));

  2773       }

  2774     }

  2776     void extractMatching() {

  2777       std::vector<int> blossoms;

  2778       for (typename BlossomSet::ClassIt c(*_blossom_set); c != INVALID; ++c) {

  2779         blossoms.push_back(c);

  2780       }

  2782       for (int i = 0; i < int(blossoms.size()); ++i) {

  2784         Value offset = (*_blossom_data)[blossoms[i]].offset;

  2785         (*_blossom_data)[blossoms[i]].pot += 2 * offset;

  2786         for (typename BlossomSet::ItemIt n(*_blossom_set, blossoms[i]);

  2787              n != INVALID; ++n) {

  2788           (*_node_data)[(*_node_index)[n]].pot -= offset;

  2789         }

  2791         Arc matching = (*_blossom_data)[blossoms[i]].next;

  2792         Node base = _graph.source(matching);

  2793         extractBlossom(blossoms[i], base, matching);

  2794       }

  2795     }

  2797   public:

  2799     /// \brief Constructor

  2800     ///

  2801     /// Constructor.

  2802     MaxWeightedPerfectMatching(const Graph& graph, const WeightMap& weight)

  2803       : _graph(graph), _weight(weight), _matching(0),

  2804         _node_potential(0), _blossom_potential(), _blossom_node_list(),

  2805         _node_num(0), _blossom_num(0),

  2807         _blossom_index(0), _blossom_set(0), _blossom_data(0),

  2808         _node_index(0), _node_heap_index(0), _node_data(0),

  2809         _tree_set_index(0), _tree_set(0),

  2811         _delta2_index(0), _delta2(0),

  2812         _delta3_index(0), _delta3(0),

  2813         _delta4_index(0), _delta4(0),

  2815         _delta_sum() {}

  2817     ~MaxWeightedPerfectMatching() {

  2818       destroyStructures();

  2819     }

  2821     /// \name Execution control

  2822     /// The simplest way to execute the algorithm is to use the

  2823     /// \c run() member function.

  2825     ///@{

  2827     /// \brief Initialize the algorithm

  2828     ///

  2829     /// Initialize the algorithm

  2830     void init() {

  2831       createStructures();

  2833       for (ArcIt e(_graph); e != INVALID; ++e) {

  2834         _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);

  2835       }

  2836       for (EdgeIt e(_graph); e != INVALID; ++e) {

  2837         _delta3_index->set(e, _delta3->PRE_HEAP);

  2838       }

  2839       for (int i = 0; i < _blossom_num; ++i) {

  2840         _delta2_index->set(i, _delta2->PRE_HEAP);

  2841         _delta4_index->set(i, _delta4->PRE_HEAP);

  2842       }

  2844       int index = 0;

  2845       for (NodeIt n(_graph); n != INVALID; ++n) {

  2846         Value max = - std::numeric_limits<Value>::max();

  2847         for (OutArcIt e(_graph, n); e != INVALID; ++e) {

  2848           if (_graph.target(e) == n) continue;

  2849           if ((dualScale * _weight[e]) / 2 > max) {

  2850             max = (dualScale * _weight[e]) / 2;

  2851           }

  2852         }

  2853         _node_index->set(n, index);

  2854         (*_node_data)[index].pot = max;

  2855         int blossom =

  2856           _blossom_set->insert(n, std::numeric_limits<Value>::max());

  2858         _tree_set->insert(blossom);

  2860         (*_blossom_data)[blossom].status = EVEN;

  2861         (*_blossom_data)[blossom].pred = INVALID;

  2862         (*_blossom_data)[blossom].next = INVALID;

  2863         (*_blossom_data)[blossom].pot = 0;

  2864         (*_blossom_data)[blossom].offset = 0;

  2865         ++index;

  2866       }

  2867       for (EdgeIt e(_graph); e != INVALID; ++e) {

  2868         int si = (*_node_index)[_graph.u(e)];

  2869         int ti = (*_node_index)[_graph.v(e)];

  2870         if (_graph.u(e) != _graph.v(e)) {

  2871           _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -

  2872                             dualScale * _weight[e]) / 2);

  2873         }

  2874       }

  2875     }

  2877     /// \brief Starts the algorithm

  2878     ///

  2879     /// Starts the algorithm

  2880     bool start() {

  2881       enum OpType {

  2882         D2, D3, D4

  2883       };

  2885       int unmatched = _node_num;

  2886       while (unmatched > 0) {

  2887         Value d2 = !_delta2->empty() ?

  2888           _delta2->prio() : std::numeric_limits<Value>::max();

  2890         Value d3 = !_delta3->empty() ?

  2891           _delta3->prio() : std::numeric_limits<Value>::max();

  2893         Value d4 = !_delta4->empty() ?

  2894           _delta4->prio() : std::numeric_limits<Value>::max();

  2896         _delta_sum = d2; OpType ot = D2;

  2897         if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }

  2898         if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }

  2900         if (_delta_sum == std::numeric_limits<Value>::max()) {

  2901           return false;

  2902         }

  2904         switch (ot) {

  2905         case D2:

  2906           {

  2907             int blossom = _delta2->top();

  2908             Node n = _blossom_set->classTop(blossom);

  2909             Arc e = (*_node_data)[(*_node_index)[n]].heap.top();

  2910             extendOnArc(e);

  2911           }

  2912           break;

  2913         case D3:

  2914           {

  2915             Edge e = _delta3->top();

  2917             int left_blossom = _blossom_set->find(_graph.u(e));

  2918             int right_blossom = _blossom_set->find(_graph.v(e));

  2920             if (left_blossom == right_blossom) {

  2921               _delta3->pop();

  2922             } else {

  2923               int left_tree = _tree_set->find(left_blossom);

  2924               int right_tree = _tree_set->find(right_blossom);

  2926               if (left_tree == right_tree) {

  2927                 shrinkOnEdge(e, left_tree);

  2928               } else {

  2929                 augmentOnEdge(e);

  2930                 unmatched -= 2;

  2931               }

  2932             }

  2933           } break;

  2934         case D4:

  2935           splitBlossom(_delta4->top());

  2936           break;

  2937         }

  2938       }

  2939       extractMatching();

  2940       return true;

  2941     }

  2943     /// \brief Runs %MaxWeightedPerfectMatching algorithm.

  2944     ///

  2945     /// This method runs the %MaxWeightedPerfectMatching algorithm.

  2946     ///

  2947     /// \note mwm.run() is just a shortcut of the following code.

  2948     /// \code

  2949     ///   mwm.init();

  2950     ///   mwm.start();

  2951     /// \endcode

  2952     bool run() {

  2953       init();

  2954       return start();

  2955     }

  2957     /// @}

  2959     /// \name Primal solution

  2960     /// Functions to get the primal solution, ie. the matching.

  2962     /// @{

  2964     /// \brief Returns the matching value.

  2965     ///

  2966     /// Returns the matching value.

  2967     Value matchingValue() const {

  2968       Value sum = 0;

  2969       for (NodeIt n(_graph); n != INVALID; ++n) {

  2970         if ((*_matching)[n] != INVALID) {

  2971           sum += _weight[(*_matching)[n]];

  2972         }

  2973       }

  2974       return sum /= 2;

  2975     }

  2977     /// \brief Returns true when the edge is in the matching.

  2978     ///

  2979     /// Returns true when the edge is in the matching.

  2980     bool matching(const Edge& edge) const {

  2981       return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;

  2982     }

  2984     /// \brief Returns the incident matching edge.

  2985     ///

  2986     /// Returns the incident matching arc from given edge.

  2987     Arc matching(const Node& node) const {

  2988       return (*_matching)[node];

  2989     }

  2991     /// \brief Returns the mate of the node.

  2992     ///

  2993     /// Returns the adjancent node in a mathcing arc.

  2994     Node mate(const Node& node) const {

  2995       return _graph.target((*_matching)[node]);

  2996     }

  2998     /// @}

  3000     /// \name Dual solution

  3001     /// Functions to get the dual solution.

  3003     /// @{

  3005     /// \brief Returns the value of the dual solution.

  3006     ///

  3007     /// Returns the value of the dual solution. It should be equal to

  3008     /// the primal value scaled by \ref dualScale "dual scale".

  3009     Value dualValue() const {

  3010       Value sum = 0;

  3011       for (NodeIt n(_graph); n != INVALID; ++n) {

  3012         sum += nodeValue(n);

  3013       }

  3014       for (int i = 0; i < blossomNum(); ++i) {

  3015         sum += blossomValue(i) * (blossomSize(i) / 2);

  3016       }

  3017       return sum;

  3018     }

  3020     /// \brief Returns the value of the node.

  3021     ///

  3022     /// Returns the the value of the node.

  3023     Value nodeValue(const Node& n) const {

  3024       return (*_node_potential)[n];

  3025     }

  3027     /// \brief Returns the number of the blossoms in the basis.

  3028     ///

  3029     /// Returns the number of the blossoms in the basis.

  3030     /// \see BlossomIt

  3031     int blossomNum() const {

  3032       return _blossom_potential.size();

  3033     }

  3036     /// \brief Returns the number of the nodes in the blossom.

  3037     ///

  3038     /// Returns the number of the nodes in the blossom.

  3039     int blossomSize(int k) const {

  3040       return _blossom_potential[k].end - _blossom_potential[k].begin;

  3041     }

  3043     /// \brief Returns the value of the blossom.

  3044     ///

  3045     /// Returns the the value of the blossom.

  3046     /// \see BlossomIt

  3047     Value blossomValue(int k) const {

  3048       return _blossom_potential[k].value;

  3049     }

  3051     /// \brief Iterator for obtaining the nodes of the blossom.

  3052     ///

  3053     /// Iterator for obtaining the nodes of the blossom. This class

  3054     /// provides a common lemon style iterator for listing a

  3055     /// subset of the nodes.

  3056     class BlossomIt {

  3057     public:

  3059       /// \brief Constructor.

  3060       ///

  3061       /// Constructor to get the nodes of the variable.

  3062       BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)

  3063         : _algorithm(&algorithm)

  3064       {

  3065         _index = _algorithm->_blossom_potential[variable].begin;

  3066         _last = _algorithm->_blossom_potential[variable].end;

  3067       }

  3069       /// \brief Conversion to node.

  3070       ///

  3071       /// Conversion to node.

  3072       operator Node() const {

  3073         return _algorithm->_blossom_node_list[_index];

  3074       }

  3076       /// \brief Increment operator.

  3077       ///

  3078       /// Increment operator.

  3079       BlossomIt& operator++() {

  3080         ++_index;

  3081         return *this;

  3082       }

  3084       /// \brief Validity checking

  3085       ///

  3086       /// Checks whether the iterator is invalid.

  3087       bool operator==(Invalid) const { return _index == _last; }

  3089       /// \brief Validity checking

  3090       ///

  3091       /// Checks whether the iterator is valid.

  3092       bool operator!=(Invalid) const { return _index != _last; }

  3094     private:

  3095       const MaxWeightedPerfectMatching* _algorithm;

  3096       int _last;

  3097       int _index;

  3098     };

  3100     /// @}

  3102   };

  3105 } //END OF NAMESPACE LEMON

  3107 #endif //LEMON_MAX_MATCHING_H

author	Akos Ladanyi <ladanyi@tmit.bme.hu>
	Wed, 01 Apr 2009 14:18:35 +0100
changeset 564	eda12d8ac953
parent 440	88ed40ad0d4f
child 581	aa1804409f29
permissions	-rw-r--r--