COIN-OR::LEMON - Graph Library

Ticket #48: a2b149cc05d2.patch

File a2b149cc05d2.patch, 109.2 KB (added by Balazs Dezso, 13 years ago)
  • lemon/Makefile.am

    # HG changeset patch
    # User Balazs Dezso <deba@inf.elte.hu>
    # Date 1223890161 -7200
    # Node ID a2b149cc05d2edb2beba5259af637f71902005c8
    # Parent  6dbd5184c6a99e936a7de899e64c114e30c0f94f
    Port maximum matching algorithms from svn 3498 (ticket #48)
    
    The MaxMatching algorithm is changed in interface and in
    internal data structures
    
    diff -r 6dbd5184c6a9 -r a2b149cc05d2 lemon/Makefile.am
    a b  
    3535        lemon/list_graph.h \
    3636        lemon/maps.h \
    3737        lemon/math.h \
     38        lemon/max_matching.h \
    3839        lemon/path.h \
    3940        lemon/random.h \
    4041        lemon/smart_graph.h \
  • new file lemon/max_matching.h

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

    diff -r 6dbd5184c6a9 -r a2b149cc05d2 test/CMakeLists.txt
    a b  
    1616  heap_test
    1717  kruskal_test
    1818  maps_test
     19  max_matching_test
    1920  random_test
    2021  path_test
    2122  time_measure_test
  • test/Makefile.am

    diff -r 6dbd5184c6a9 -r a2b149cc05d2 test/Makefile.am
    a b  
    1919        test/heap_test \
    2020        test/kruskal_test \
    2121        test/maps_test \
     22        test/max_matching_test \
    2223        test/random_test \
    2324        test/path_test \
    2425        test/test_tools_fail \
     
    4243test_heap_test_SOURCES = test/heap_test.cc
    4344test_kruskal_test_SOURCES = test/kruskal_test.cc
    4445test_maps_test_SOURCES = test/maps_test.cc
     46test_max_matching_test_SOURCES = test/max_matching_test.cc
    4547test_path_test_SOURCES = test/path_test.cc
    4648test_random_test_SOURCES = test/random_test.cc
    4749test_test_tools_fail_SOURCES = test/test_tools_fail.cc
  • new file test/max_matching_test.cc

    diff -r 6dbd5184c6a9 -r a2b149cc05d2 test/max_matching_test.cc
    - +  
     1/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2 *
     3 * This file is a part of LEMON, a generic C++ optimization library.
     4 *
     5 * Copyright (C) 2003-2008
     6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8 *
     9 * Permission to use, modify and distribute this software is granted
     10 * provided that this copyright notice appears in all copies. For
     11 * precise terms see the accompanying LICENSE file.
     12 *
     13 * This software is provided "AS IS" with no warranty of any kind,
     14 * express or implied, and with no claim as to its suitability for any
     15 * purpose.
     16 *
     17 */
     18
     19#include <iostream>
     20#include <sstream>
     21#include <vector>
     22#include <queue>
     23#include <lemon/math.h>
     24#include <cstdlib>
     25
     26#include <lemon/max_matching.h>
     27#include <lemon/smart_graph.h>
     28#include <lemon/lgf_reader.h>
     29
     30#include "test_tools.h"
     31
     32using namespace std;
     33using namespace lemon;
     34
     35GRAPH_TYPEDEFS(SmartGraph);
     36
     37
     38const int lgfn = 3;
     39const std::string lgf[lgfn] = {
     40  "@nodes\n"
     41  "label\n"
     42  "0\n"
     43  "1\n"
     44  "2\n"
     45  "3\n"
     46  "4\n"
     47  "5\n"
     48  "6\n"
     49  "7\n"
     50  "@edges\n"
     51  "     label  weight\n"
     52  "7 4  0      984\n"
     53  "0 7  1      73\n"
     54  "7 1  2      204\n"
     55  "2 3  3      583\n"
     56  "2 7  4      565\n"
     57  "2 1  5      582\n"
     58  "0 4  6      551\n"
     59  "2 5  7      385\n"
     60  "1 5  8      561\n"
     61  "5 3  9      484\n"
     62  "7 5  10     904\n"
     63  "3 6  11     47\n"
     64  "7 6  12     888\n"
     65  "3 0  13     747\n"
     66  "6 1  14     310\n",
     67
     68  "@nodes\n"
     69  "label\n"
     70  "0\n"
     71  "1\n"
     72  "2\n"
     73  "3\n"
     74  "4\n"
     75  "5\n"
     76  "6\n"
     77  "7\n"
     78  "@edges\n"
     79  "     label  weight\n"
     80  "2 5  0      710\n"
     81  "0 5  1      241\n"
     82  "2 4  2      856\n"
     83  "2 6  3      762\n"
     84  "4 1  4      747\n"
     85  "6 1  5      962\n"
     86  "4 7  6      723\n"
     87  "1 7  7      661\n"
     88  "2 3  8      376\n"
     89  "1 0  9      416\n"
     90  "6 7  10     391\n",
     91
     92  "@nodes\n"
     93  "label\n"
     94  "0\n"
     95  "1\n"
     96  "2\n"
     97  "3\n"
     98  "4\n"
     99  "5\n"
     100  "6\n"
     101  "7\n"
     102  "@edges\n"
     103  "     label  weight\n"
     104  "6 2  0      553\n"
     105  "0 7  1      653\n"
     106  "6 3  2      22\n"
     107  "4 7  3      846\n"
     108  "7 2  4      981\n"
     109  "7 6  5      250\n"
     110  "5 2  6      539\n",
     111};
     112
     113void checkMatching(const SmartGraph& graph,
     114                   const MaxMatching<SmartGraph>& mm) {
     115  int num = 0;
     116
     117  IntNodeMap comp_index(graph);
     118  UnionFind<IntNodeMap> comp(comp_index);
     119
     120  int barrier_num = 0;
     121
     122  for (NodeIt n(graph); n != INVALID; ++n) {
     123    check(mm.decomposition(n) == MaxMatching<SmartGraph>::EVEN ||
     124          mm.matching(n) != INVALID, "Wrong Gallai-Edmonds decomposition");
     125    if (mm.decomposition(n) == MaxMatching<SmartGraph>::ODD) {
     126      ++barrier_num;
     127    } else {
     128      comp.insert(n);
     129    }
     130  }
     131
     132  for (EdgeIt e(graph); e != INVALID; ++e) {
     133    if (mm.matching(e)) {
     134      check(e == mm.matching(graph.u(e)), "Wrong matching");
     135      check(e == mm.matching(graph.v(e)), "Wrong matching");
     136      ++num;
     137    }
     138    check(mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::EVEN ||
     139          mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::MATCHED,
     140          "Wrong Gallai-Edmonds decomposition");
     141
     142    check(mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::EVEN ||
     143          mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::MATCHED,
     144          "Wrong Gallai-Edmonds decomposition");
     145
     146    if (mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::ODD &&
     147        mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::ODD) {
     148      comp.join(graph.u(e), graph.v(e));
     149    }
     150  }
     151
     152  std::set<int> comp_root;
     153  int odd_comp_num = 0;
     154  for (NodeIt n(graph); n != INVALID; ++n) {
     155    if (mm.decomposition(n) != MaxMatching<SmartGraph>::ODD) {
     156      int root = comp.find(n);
     157      if (comp_root.find(root) == comp_root.end()) {
     158        comp_root.insert(root);
     159        if (comp.size(n) % 2 == 1) {
     160          ++odd_comp_num;
     161        }
     162      }
     163    }
     164  }
     165
     166  check(mm.matchingSize() == num, "Wrong matching");
     167  check(2 * num == countNodes(graph) - (odd_comp_num - barrier_num),
     168         "Wrong matching");
     169  return;
     170}
     171
     172void checkWeightedMatching(const SmartGraph& graph,
     173                   const SmartGraph::EdgeMap<int>& weight,
     174                   const MaxWeightedMatching<SmartGraph>& mwm) {
     175  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
     176    if (graph.u(e) == graph.v(e)) continue;
     177    int rw = mwm.nodeValue(graph.u(e)) + mwm.nodeValue(graph.v(e));
     178
     179    for (int i = 0; i < mwm.blossomNum(); ++i) {
     180      bool s = false, t = false;
     181      for (MaxWeightedMatching<SmartGraph>::BlossomIt n(mwm, i);
     182           n != INVALID; ++n) {
     183        if (graph.u(e) == n) s = true;
     184        if (graph.v(e) == n) t = true;
     185      }
     186      if (s == true && t == true) {
     187        rw += mwm.blossomValue(i);
     188      }
     189    }
     190    rw -= weight[e] * mwm.dualScale;
     191
     192    check(rw >= 0, "Negative reduced weight");
     193    check(rw == 0 || !mwm.matching(e),
     194          "Non-zero reduced weight on matching edge");
     195  }
     196
     197  int pv = 0;
     198  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     199    if (mwm.matching(n) != INVALID) {
     200      check(mwm.nodeValue(n) >= 0, "Invalid node value");
     201      pv += weight[mwm.matching(n)];
     202      SmartGraph::Node o = graph.target(mwm.matching(n));
     203      check(mwm.mate(n) == o, "Invalid matching");
     204      check(mwm.matching(n) == graph.oppositeArc(mwm.matching(o)),
     205            "Invalid matching");
     206    } else {
     207      check(mwm.mate(n) == INVALID, "Invalid matching");
     208      check(mwm.nodeValue(n) == 0, "Invalid matching");
     209    }
     210  }
     211
     212  int dv = 0;
     213  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     214    dv += mwm.nodeValue(n);
     215  }
     216
     217  for (int i = 0; i < mwm.blossomNum(); ++i) {
     218    check(mwm.blossomValue(i) >= 0, "Invalid blossom value");
     219    check(mwm.blossomSize(i) % 2 == 1, "Even blossom size");
     220    dv += mwm.blossomValue(i) * ((mwm.blossomSize(i) - 1) / 2);
     221  }
     222
     223  check(pv * mwm.dualScale == dv * 2, "Wrong duality");
     224
     225  return;
     226}
     227
     228void checkWeightedPerfectMatching(const SmartGraph& graph,
     229                          const SmartGraph::EdgeMap<int>& weight,
     230                          const MaxWeightedPerfectMatching<SmartGraph>& mwpm) {
     231  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
     232    if (graph.u(e) == graph.v(e)) continue;
     233    int rw = mwpm.nodeValue(graph.u(e)) + mwpm.nodeValue(graph.v(e));
     234
     235    for (int i = 0; i < mwpm.blossomNum(); ++i) {
     236      bool s = false, t = false;
     237      for (MaxWeightedPerfectMatching<SmartGraph>::BlossomIt n(mwpm, i);
     238           n != INVALID; ++n) {
     239        if (graph.u(e) == n) s = true;
     240        if (graph.v(e) == n) t = true;
     241      }
     242      if (s == true && t == true) {
     243        rw += mwpm.blossomValue(i);
     244      }
     245    }
     246    rw -= weight[e] * mwpm.dualScale;
     247
     248    check(rw >= 0, "Negative reduced weight");
     249    check(rw == 0 || !mwpm.matching(e),
     250          "Non-zero reduced weight on matching edge");
     251  }
     252
     253  int pv = 0;
     254  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     255    check(mwpm.matching(n) != INVALID, "Non perfect");
     256    pv += weight[mwpm.matching(n)];
     257    SmartGraph::Node o = graph.target(mwpm.matching(n));
     258    check(mwpm.mate(n) == o, "Invalid matching");
     259    check(mwpm.matching(n) == graph.oppositeArc(mwpm.matching(o)),
     260          "Invalid matching");
     261  }
     262
     263  int dv = 0;
     264  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
     265    dv += mwpm.nodeValue(n);
     266  }
     267
     268  for (int i = 0; i < mwpm.blossomNum(); ++i) {
     269    check(mwpm.blossomValue(i) >= 0, "Invalid blossom value");
     270    check(mwpm.blossomSize(i) % 2 == 1, "Even blossom size");
     271    dv += mwpm.blossomValue(i) * ((mwpm.blossomSize(i) - 1) / 2);
     272  }
     273
     274  check(pv * mwpm.dualScale == dv * 2, "Wrong duality");
     275
     276  return;
     277}
     278
     279
     280int main() {
     281
     282  for (int i = 0; i < lgfn; ++i) {
     283    SmartGraph graph;
     284    SmartGraph::EdgeMap<int> weight(graph);
     285
     286    istringstream lgfs(lgf[i]);
     287    graphReader(graph, lgfs).
     288      edgeMap("weight", weight).run();
     289
     290    MaxMatching<SmartGraph> mm(graph);
     291    mm.run();
     292    checkMatching(graph, mm);
     293
     294    MaxWeightedMatching<SmartGraph> mwm(graph, weight);
     295    mwm.run();
     296    checkWeightedMatching(graph, weight, mwm);
     297
     298    MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
     299    bool perfect = mwpm.run();
     300
     301    check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
     302          "Perfect matching found");
     303
     304    if (perfect) {
     305      checkWeightedPerfectMatching(graph, weight, mwpm);
     306    }
     307  }
     308
     309  return 0;
     310}