Ticket #48: 91d63b8b1a4c.patch

lemon/max_matching.h

@@ -31 +31 @@
 
 ///\ingroup matching
 ///\file
-///\brief Maximum matching algorithms in graph.
+///\brief Maximum matching algorithms in general graphs.
 
 namespace lemon {
 
-  ///\ingroup matching
+  /// \ingroup matching
   ///
-  ///\brief Edmonds' alternating forest maximum matching algorithm.
+  /// \brief Edmonds' alternating forest maximum matching algorithm.
   ///
-  ///This class provides Edmonds' alternating forest matching
-  ///algorithm. The starting matching (if any) can be passed to the
-  ///algorithm using some of init functions.
+  /// This class provides Edmonds' alternating forest matching
+  /// algorithm. The starting matching (if any) can be passed to the
+  /// algorithm using some of init functions.
   ///
-  ///The dual side of a matching is a map of the nodes to
-  ///MaxMatching::DecompType, having values \c D, \c A and \c C
-  ///showing the Gallai-Edmonds decomposition of the digraph. The nodes
-  ///in \c D induce a digraph with factor-critical components, the nodes
-  ///in \c A form the barrier, and the nodes in \c C induce a digraph
-  ///having a perfect matching. This decomposition can be attained by
-  ///calling \c decomposition() after running the algorithm.
+  /// The dual side of a matching is a map of the nodes to
+  /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c
+  /// MATCHED/C showing the Gallai-Edmonds decomposition of the
+  /// graph. The nodes in \c EVEN/D induce a graph with
+  /// factor-critical components, the nodes in \c ODD/A form the
+  /// barrier, and the nodes in \c MATCHED/C induce a graph having a
+  /// perfect matching. The number of the fractor critical components
+  /// minus the number of barrier nodes is a lower bound on the
+  /// unmatched nodes, and if the matching is optimal this bound is
+  /// tight. This decomposition can be attained by calling \c
+  /// decomposition() after running the algorithm.
   ///
-  ///\param Digraph The graph type the algorithm runs on.
-  template <typename Graph>
+  /// \param _Graph The graph type the algorithm runs on.
+  template <typename _Graph>
   class MaxMatching {
+  public:
 
-  protected:
+    typedef _Graph Graph;
+    typedef typename Graph::template NodeMap<typename Graph::Arc>
+    MatchingMap;
+
+    ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
+    ///
+    ///Indicates the Gallai-Edmonds decomposition of the graph, which
+    ///shows an upper bound on the size of a maximum matching. The
+    ///nodes with Status \c EVEN/D induce a graph with factor-critical
+    ///components, the nodes in \c ODD/A form the canonical barrier,
+    ///and the nodes in \c MATCHED/C induce a graph having a perfect
+    ///matching.
+    enum Status {
+      EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2
+    };
+
+    typedef typename Graph::template NodeMap<Status> StatusMap;
+
+  private:
 
     TEMPLATE_GRAPH_TYPEDEFS(Graph);
 
-    typedef typename Graph::template NodeMap<int> UFECrossRef;
-    typedef UnionFindEnum<UFECrossRef> UFE;
-    typedef std::vector<Node> NV;
+    typedef UnionFindEnum<IntNodeMap> BlossomSet;
+    typedef ExtendFindEnum<IntNodeMap> TreeSet;
+    typedef RangeMap<Node> NodeIntMap;
+    typedef MatchingMap EarMap;
+    typedef std::vector<Node> NodeQueue;
 
-    typedef typename Graph::template NodeMap<int> EFECrossRef;
-    typedef ExtendFindEnum<EFECrossRef> EFE;
+    const Graph& _graph;
+    MatchingMap* _matching;
+    StatusMap* _status;
+
+    EarMap* _ear;
+
+    IntNodeMap* _blossom_set_index;
+    BlossomSet* _blossom_set;
+    NodeIntMap* _blossom_rep;
+
+    IntNodeMap* _tree_set_index;
+    TreeSet* _tree_set;
+
+    NodeQueue _node_queue;
+    int _process, _postpone, _last;
+
+    int _node_num;
+
+  private:
+
+    void createStructures() {
+      _node_num = countNodes(_graph);
+      if (!_matching) {
+        _matching = new MatchingMap(_graph);
+      }
+      if (!_status) {
+        _status = new StatusMap(_graph);
+      }
+      if (!_ear) {
+        _ear = new EarMap(_graph);
+      }
+      if (!_blossom_set) {
+        _blossom_set_index = new IntNodeMap(_graph);
+        _blossom_set = new BlossomSet(*_blossom_set_index);
+      }
+      if (!_blossom_rep) {
+        _blossom_rep = new NodeIntMap(_node_num);
+      }
+      if (!_tree_set) {
+        _tree_set_index = new IntNodeMap(_graph);
+        _tree_set = new TreeSet(*_tree_set_index);
+      }
+      _node_queue.resize(_node_num);
+    }
+
+    void destroyStructures() {
+      if (_matching) {
+        delete _matching;
+      }
+      if (_status) {
+        delete _status;
+      }
+      if (_ear) {
+        delete _ear;
+      }
+      if (_blossom_set) {
+        delete _blossom_set;
+        delete _blossom_set_index;
+      }
+      if (_blossom_rep) {
+        delete _blossom_rep;
+      }
+      if (_tree_set) {
+        delete _tree_set_index;
+        delete _tree_set;
+      }
+    }
+
+    void processDense(const Node& n) {
+      _process = _postpone = _last = 0;
+      _node_queue[_last++] = n;
+
+      while (_process != _last) {
+        Node u = _node_queue[_process++];
+        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
+          Node v = _graph.target(a);
+          if ((*_status)[v] == MATCHED) {
+            extendOnArc(a);
+          } else if ((*_status)[v] == UNMATCHED) {
+            augmentOnArc(a);
+            return;
+          }
+        }
+      }
+
+      while (_postpone != _last) {
+        Node u = _node_queue[_postpone++];
+
+        for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
+          Node v = _graph.target(a);
+
+          if ((*_status)[v] == EVEN) {
+            if (_blossom_set->find(u) != _blossom_set->find(v)) {
+              shrinkOnEdge(a);
+            }
+          }
+
+          while (_process != _last) {
+            Node w = _node_queue[_process++];
+            for (OutArcIt b(_graph, w); b != INVALID; ++b) {
+              Node x = _graph.target(b);
+              if ((*_status)[x] == MATCHED) {
+                extendOnArc(b);
+              } else if ((*_status)[x] == UNMATCHED) {
+                augmentOnArc(b);
+                return;
+              }
+            }
+          }
+        }
+      }
+    }
+
+    void processSparse(const Node& n) {
+      _process = _last = 0;
+      _node_queue[_last++] = n;
+      while (_process != _last) {
+        Node u = _node_queue[_process++];
+        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
+          Node v = _graph.target(a);
+
+          if ((*_status)[v] == EVEN) {
+            if (_blossom_set->find(u) != _blossom_set->find(v)) {
+              shrinkOnEdge(a);
+            }
+          } else if ((*_status)[v] == MATCHED) {
+            extendOnArc(a);
+          } else if ((*_status)[v] == UNMATCHED) {
+            augmentOnArc(a);
+            return;
+          }
+        }
+      }
+    }
+
+    void shrinkOnEdge(const Edge& e) {
+      Node nca = INVALID;
+
+      {
+        std::set<Node> left_set, right_set;
+
+        Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
+        left_set.insert(left);
+
+        Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
+        right_set.insert(right);
+
+        while (true) {
+          if ((*_matching)[left] == INVALID) break;
+          left = _graph.target((*_matching)[left]);
+          left = (*_blossom_rep)[_blossom_set->
+                                 find(_graph.target((*_ear)[left]))];
+          if (right_set.find(left) != right_set.end()) {
+            nca = left;
+            break;
+          }
+          left_set.insert(left);
+
+          if ((*_matching)[right] == INVALID) break;
+          right = _graph.target((*_matching)[right]);
+          right = (*_blossom_rep)[_blossom_set->
+                                  find(_graph.target((*_ear)[right]))];
+          if (left_set.find(right) != left_set.end()) {
+            nca = right;
+            break;
+          }
+          right_set.insert(right);
+        }
+
+        if (nca == INVALID) {
+          if ((*_matching)[left] == INVALID) {
+            nca = right;
+            while (left_set.find(nca) == left_set.end()) {
+              nca = _graph.target((*_matching)[nca]);
+              nca =(*_blossom_rep)[_blossom_set->
+                                   find(_graph.target((*_ear)[nca]))];
+            }
+          } else {
+            nca = left;
+            while (right_set.find(nca) == right_set.end()) {
+              nca = _graph.target((*_matching)[nca]);
+              nca = (*_blossom_rep)[_blossom_set->
+                                   find(_graph.target((*_ear)[nca]))];
+            }
+          }
+        }
+      }
+
+      {
+
+        Node node = _graph.u(e);
+        Arc arc = _graph.direct(e, true);
+        Node base = (*_blossom_rep)[_blossom_set->find(node)];
+
+        while (base != nca) {
+          _ear->set(node, arc);
+
+          Node n = node;
+          while (n != base) {
+            n = _graph.target((*_matching)[n]);
+            Arc a = (*_ear)[n];
+            n = _graph.target(a);
+            _ear->set(n, _graph.oppositeArc(a));
+          }
+          node = _graph.target((*_matching)[base]);
+          _tree_set->erase(base);
+          _tree_set->erase(node);
+          _blossom_set->insert(node, _blossom_set->find(base));
+          _status->set(node, EVEN);
+          _node_queue[_last++] = node;
+          arc = _graph.oppositeArc((*_ear)[node]);
+          node = _graph.target((*_ear)[node]);
+          base = (*_blossom_rep)[_blossom_set->find(node)];
+          _blossom_set->join(_graph.target(arc), base);
+        }
+      }
+
+      _blossom_rep->set(_blossom_set->find(nca), nca);
+
+      {
+
+        Node node = _graph.v(e);
+        Arc arc = _graph.direct(e, false);
+        Node base = (*_blossom_rep)[_blossom_set->find(node)];
+
+        while (base != nca) {
+          _ear->set(node, arc);
+
+          Node n = node;
+          while (n != base) {
+            n = _graph.target((*_matching)[n]);
+            Arc a = (*_ear)[n];
+            n = _graph.target(a);
+            _ear->set(n, _graph.oppositeArc(a));
+          }
+          node = _graph.target((*_matching)[base]);
+          _tree_set->erase(base);
+          _tree_set->erase(node);
+          _blossom_set->insert(node, _blossom_set->find(base));
+          _status->set(node, EVEN);
+          _node_queue[_last++] = node;
+          arc = _graph.oppositeArc((*_ear)[node]);
+          node = _graph.target((*_ear)[node]);
+          base = (*_blossom_rep)[_blossom_set->find(node)];
+          _blossom_set->join(_graph.target(arc), base);
+        }
+      }
+
+      _blossom_rep->set(_blossom_set->find(nca), nca);
+    }
+
+
+
+    void extendOnArc(const Arc& a) {
+      Node base = _graph.source(a);
+      Node odd = _graph.target(a);
+
+      _ear->set(odd, _graph.oppositeArc(a));
+      Node even = _graph.target((*_matching)[odd]);
+      _blossom_rep->set(_blossom_set->insert(even), even);
+      _status->set(odd, ODD);
+      _status->set(even, EVEN);
+      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
+      _tree_set->insert(odd, tree);
+      _tree_set->insert(even, tree);
+      _node_queue[_last++] = even;
+
+    }
+
+    void augmentOnArc(const Arc& a) {
+      Node even = _graph.source(a);
+      Node odd = _graph.target(a);
+
+      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
+
+      _matching->set(odd, _graph.oppositeArc(a));
+      _status->set(odd, MATCHED);
+
+      Arc arc = (*_matching)[even];
+      _matching->set(even, a);
+
+      while (arc != INVALID) {
+        odd = _graph.target(arc);
+        arc = (*_ear)[odd];
+        even = _graph.target(arc);
+        _matching->set(odd, arc);
+        arc = (*_matching)[even];
+        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
+      }
+
+      for (typename TreeSet::ItemIt it(*_tree_set, tree);
+           it != INVALID; ++it) {
+        if ((*_status)[it] == ODD) {
+          _status->set(it, MATCHED);
+        } else {
+          int blossom = _blossom_set->find(it);
+          for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
+               jt != INVALID; ++jt) {
+            _status->set(jt, MATCHED);
+          }
+          _blossom_set->eraseClass(blossom);
+        }
+      }
+      _tree_set->eraseClass(tree);
+
+    }
 
   public:
 
-    ///\brief Indicates the Gallai-Edmonds decomposition of the digraph.
+    /// \brief Constructor
     ///
-    ///Indicates the Gallai-Edmonds decomposition of the digraph, which
-    ///shows an upper bound on the size of a maximum matching. The
-    ///nodes with DecompType \c D induce a digraph with factor-critical
-    ///components, the nodes in \c A form the canonical barrier, and the
-    ///nodes in \c C induce a digraph having a perfect matching.
-    enum DecompType {
-      D=0,
-      A=1,
-      C=2
-    };
+    /// Constructor.
+    MaxMatching(const Graph& graph)
+      : _graph(graph), _matching(0), _status(0), _ear(0),
+        _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
+        _tree_set_index(0), _tree_set(0) {}
 
-  protected:
+    ~MaxMatching() {
+      destroyStructures();
+    }
 
-    static const int HEUR_density=2;
-    const Graph& g;
-    typename Graph::template NodeMap<Node> _mate;
-    typename Graph::template NodeMap<DecompType> position;
+    /// \name Execution control
+    /// The simplest way to execute the algorithm is to use the member
+    /// \c run() member function.
+    /// \n
 
-  public:
+    /// If you need more control on the execution, first you must call
+    /// \ref init(), \ref greedyInit() or \ref matchingInit()
+    /// functions, then you can start the algorithm with the \ref
+    /// startParse() or startDense() functions.
 
-    MaxMatching(const Graph& _g)
-      : g(_g), _mate(_g), position(_g) {}
+    ///@{
 
-    ///\brief Sets the actual matching to the empty matching.
+    /// \brief Sets the actual matching to the empty matching.
     ///
-    ///Sets the actual matching to the empty matching.
+    /// Sets the actual matching to the empty matching.
     ///
     void init() {
-      for(NodeIt v(g); v!=INVALID; ++v) {
-        _mate.set(v,INVALID);
-        position.set(v,C);
+      createStructures();
+      for(NodeIt n(_graph); n != INVALID; ++n) {
+        _matching->set(n, INVALID);
+        _status->set(n, UNMATCHED);
       }
     }
 
@@ -108 +437 @@
     ///
     ///For initial matchig it finds a maximal greedy matching.
     void greedyInit() {
-      for(NodeIt v(g); v!=INVALID; ++v) {
-        _mate.set(v,INVALID);
-        position.set(v,C);
+      createStructures();
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        _matching->set(n, INVALID);
+        _status->set(n, UNMATCHED);
       }
-      for(NodeIt v(g); v!=INVALID; ++v)
-        if ( _mate[v]==INVALID ) {
-          for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) {
-            Node y=g.runningNode(e);
-            if ( _mate[y]==INVALID && y!=v ) {
-              _mate.set(v,y);
-              _mate.set(y,v);
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        if ((*_matching)[n] == INVALID) {
+          for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
+            Node v = _graph.target(a);
+            if ((*_matching)[v] == INVALID && v != n) {
+              _matching->set(n, a);
+              _status->set(n, MATCHED);
+              _matching->set(v, _graph.oppositeArc(a));
+              _status->set(v, MATCHED);
               break;
             }
           }
         }
-    }
-
-    ///\brief Initialize the matching from each nodes' mate.
-    ///
-    ///Initialize the matching from a \c Node valued \c Node map. This
-    ///map must be \e symmetric, i.e. if \c map[u]==v then \c
-    ///map[v]==u must hold, and \c uv will be an arc of the initial
-    ///matching.
-    template <typename MateMap>
-    void mateMapInit(MateMap& map) {
-      for(NodeIt v(g); v!=INVALID; ++v) {
-        _mate.set(v,map[v]);
-        position.set(v,C);
       }
     }
 
-    ///\brief Initialize the matching from a node map with the
-    ///incident matching arcs.
+
+    /// \brief Initialize the matching from the map containing a matching.
     ///
-    ///Initialize the matching from an \c Edge valued \c Node map. \c
-    ///map[v] must be an \c Edge incident to \c v. This map must have
-    ///the property that if \c g.oppositeNode(u,map[u])==v then \c \c
-    ///g.oppositeNode(v,map[v])==u holds, and now some arc joining \c
-    ///u to \c v will be an arc of the matching.
-    template<typename MatchingMap>
-    void matchingMapInit(MatchingMap& map) {
-      for(NodeIt v(g); v!=INVALID; ++v) {
-        position.set(v,C);
-        Edge e=map[v];
-        if ( e!=INVALID )
-          _mate.set(v,g.oppositeNode(v,e));
-        else
-          _mate.set(v,INVALID);
+    /// Initialize the matching from a \c bool valued \c Edge map. This
+    /// map must have the property that there are no two incident edges
+    /// with true value, ie. it contains a matching.
+    /// \return %True if the map contains a matching.
+    template <typename MatchingMap>
+    bool matchingInit(const MatchingMap& matching) {
+      createStructures();
+
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        _matching->set(n, INVALID);
+        _status->set(n, UNMATCHED);
       }
+      for(EdgeIt e(_graph); e!=INVALID; ++e) {
+        if (matching[e]) {
+
+          Node u = _graph.u(e);
+          if ((*_matching)[u] != INVALID) return false;
+          _matching->set(u, _graph.direct(e, true));
+          _status->set(u, MATCHED);
+
+          Node v = _graph.v(e);
+          if ((*_matching)[v] != INVALID) return false;
+          _matching->set(v, _graph.direct(e, false));
+          _status->set(v, MATCHED);
+        }
+      }
+      return true;
     }
 
-    ///\brief Initialize the matching from the map containing the
-    ///undirected matching arcs.
+    /// \brief Starts Edmonds' algorithm
     ///
-    ///Initialize the matching from a \c bool valued \c Edge map. This
-    ///map must have the property that there are no two incident arcs
-    ///\c e, \c f with \c map[e]==map[f]==true. The arcs \c e with \c
-    ///map[e]==true form the matching.
-    template <typename MatchingMap>
-    void matchingInit(MatchingMap& map) {
-      for(NodeIt v(g); v!=INVALID; ++v) {
-        _mate.set(v,INVALID);
-        position.set(v,C);
-      }
-      for(EdgeIt e(g); e!=INVALID; ++e) {
-        if ( map[e] ) {
-          Node u=g.u(e);
-          Node v=g.v(e);
-          _mate.set(u,v);
-          _mate.set(v,u);
+    /// If runs the original Edmonds' algorithm.
+    void startSparse() {
+      for(NodeIt n(_graph); n != INVALID; ++n) {
+        if ((*_status)[n] == UNMATCHED) {
+          (*_blossom_rep)[_blossom_set->insert(n)] = n;
+          _tree_set->insert(n);
+          _status->set(n, EVEN);
+          processSparse(n);
         }
       }
     }
 
+    /// \brief Starts Edmonds' algorithm.
+    ///
+    /// It runs Edmonds' algorithm with a heuristic of postponing
+    /// shrinks, giving a faster algorithm for dense graphs.
+    void startDense() {
+      for(NodeIt n(_graph); n != INVALID; ++n) {
+        if ((*_status)[n] == UNMATCHED) {
+          (*_blossom_rep)[_blossom_set->insert(n)] = n;
+          _tree_set->insert(n);
+          _status->set(n, EVEN);
+          processDense(n);
+        }
+      }
+    }
 
-    ///\brief Runs Edmonds' algorithm.
+
+    /// \brief Runs Edmonds' algorithm
     ///
-    ///Runs Edmonds' algorithm for sparse digraphs (number of arcs <
-    ///2*number of nodes), and a heuristical Edmonds' algorithm with a
-    ///heuristic of postponing shrinks for dense digraphs.
+    /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>)
+    /// or Edmonds' algorithm with a heuristic of
+    /// postponing shrinks for dense graphs.
     void run() {
-      if (countEdges(g) < HEUR_density * countNodes(g)) {
+      if (countEdges(_graph) < 2 * countNodes(_graph)) {
         greedyInit();
         startSparse();
       } else {
@@ -198 +535 @@
       }
     }
 
+    /// @}
 
-    ///\brief Starts Edmonds' algorithm.
-    ///
-    ///If runs the original Edmonds' algorithm.
-    void startSparse() {
+    /// \name Primal solution
+    /// Functions for get the primal solution, ie. the matching.
 
-      typename Graph::template NodeMap<Node> ear(g,INVALID);
-      //undefined for the base nodes of the blossoms (i.e. for the
-      //representative elements of UFE blossom) and for the nodes in C
-
-      UFECrossRef blossom_base(g);
-      UFE blossom(blossom_base);
-      NV rep(countNodes(g));
-
-      EFECrossRef tree_base(g);
-      EFE tree(tree_base);
-
-      //If these UFE's would be members of the class then also
-      //blossom_base and tree_base should be a member.
-
-      //We build only one tree and the other vertices uncovered by the
-      //matching belong to C. (They can be considered as singleton
-      //trees.) If this tree can be augmented or no more
-      //grow/augmentation/shrink is possible then we return to this
-      //"for" cycle.
-      for(NodeIt v(g); v!=INVALID; ++v) {
-        if (position[v]==C && _mate[v]==INVALID) {
-          rep[blossom.insert(v)] = v;
-          tree.insert(v);
-          position.set(v,D);
-          normShrink(v, ear, blossom, rep, tree);
-        }
-      }
-    }
-
-    ///\brief Starts Edmonds' algorithm.
-    ///
-    ///It runs Edmonds' algorithm with a heuristic of postponing
-    ///shrinks, giving a faster algorithm for dense digraphs.
-    void startDense() {
-
-      typename Graph::template NodeMap<Node> ear(g,INVALID);
-      //undefined for the base nodes of the blossoms (i.e. for the
-      //representative elements of UFE blossom) and for the nodes in C
-
-      UFECrossRef blossom_base(g);
-      UFE blossom(blossom_base);
-      NV rep(countNodes(g));
-
-      EFECrossRef tree_base(g);
-      EFE tree(tree_base);
-
-      //If these UFE's would be members of the class then also
-      //blossom_base and tree_base should be a member.
-
-      //We build only one tree and the other vertices uncovered by the
-      //matching belong to C. (They can be considered as singleton
-      //trees.) If this tree can be augmented or no more
-      //grow/augmentation/shrink is possible then we return to this
-      //"for" cycle.
-      for(NodeIt v(g); v!=INVALID; ++v) {
-        if ( position[v]==C && _mate[v]==INVALID ) {
-          rep[blossom.insert(v)] = v;
-          tree.insert(v);
-          position.set(v,D);
-          lateShrink(v, ear, blossom, rep, tree);
-        }
-      }
-    }
-
-
+    /// @{
 
     ///\brief Returns the size of the actual matching stored.
     ///
     ///Returns the size of the actual matching stored. After \ref
-    ///run() it returns the size of a maximum matching in the digraph.
-    int size() const {
-      int s=0;
-      for(NodeIt v(g); v!=INVALID; ++v) {
-        if ( _mate[v]!=INVALID ) {
-          ++s;
+    ///run() it returns the size of the maximum matching in the graph.
+    int matchingSize() const {
+      int size = 0;
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        if ((*_matching)[n] != INVALID) {
+          ++size;
         }
       }
-      return s/2;
+      return size / 2;
     }
 
+    /// \brief Returns true when the edge is in the matching.
+    ///
+    /// Returns true when the edge is in the matching.
+    bool matching(const Edge& edge) const {
+      return edge == (*_matching)[_graph.u(edge)];
+    }
+
+    /// \brief Returns the matching edge incident to the given node.
+    ///
+    /// Returns the matching edge of a \c node in the actual matching or
+    /// INVALID if the \c node is not covered by the actual matching.
+    Arc matching(const Node& n) const {
+      return (*_matching)[n];
+    }
 
     ///\brief Returns the mate of a node in the actual matching.
     ///
-    ///Returns the mate of a \c node in the actual matching.
-    ///Returns INVALID if the \c node is not covered by the actual matching.
-    Node mate(const Node& node) const {
-      return _mate[node];
+    ///Returns the mate of a \c node in the actual matching or
+    ///INVALID if the \c node is not covered by the actual matching.
+    Node mate(const Node& n) const {
+      return (*_matching)[n] != INVALID ?
+        _graph.target((*_matching)[n]) : INVALID;
     }
 
-    ///\brief Returns the matching arc incident to the given node.
-    ///
-    ///Returns the matching arc of a \c node in the actual matching.
-    ///Returns INVALID if the \c node is not covered by the actual matching.
-    Edge matchingArc(const Node& node) const {
-      if (_mate[node] == INVALID) return INVALID;
-      Node n = node < _mate[node] ? node : _mate[node];
-      for (IncEdgeIt e(g, n); e != INVALID; ++e) {
-        if (g.oppositeNode(n, e) == _mate[n]) {
-          return e;
-        }
-      }
-      return INVALID;
-    }
+    /// @}
+
+    /// \name Dual solution
+    /// Functions for get the dual solution, ie. the decomposition.
+
+    /// @{
 
     /// \brief Returns the class of the node in the Edmonds-Gallai
     /// decomposition.
     ///
     /// Returns the class of the node in the Edmonds-Gallai
     /// decomposition.
-    DecompType decomposition(const Node& n) {
-      return position[n] == A;
+    Status decomposition(const Node& n) const {
+      return (*_status)[n];
     }
 
     /// \brief Returns true when the node is in the barrier.
     ///
     /// Returns true when the node is in the barrier.
-    bool barrier(const Node& n) {
-      return position[n] == A;
+    bool barrier(const Node& n) const {
+      return (*_status)[n] == ODD;
     }
 
-    ///\brief Gives back the matching in a \c Node of mates.
-    ///
-    ///Writes the stored matching to a \c Node valued \c Node map. The
-    ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c
-    ///map[v]==u will hold, and now \c uv is an arc of the matching.
-    template <typename MateMap>
-    void mateMap(MateMap& map) const {
-      for(NodeIt v(g); v!=INVALID; ++v) {
-        map.set(v,_mate[v]);
-      }
-    }
-
-    ///\brief Gives back the matching in an \c Edge valued \c Node
-    ///map.
-    ///
-    ///Writes the stored matching to an \c Edge valued \c Node
-    ///map. \c map[v] will be an \c Edge incident to \c v. This
-    ///map will have the property that if \c g.oppositeNode(u,map[u])
-    ///== v then \c map[u]==map[v] holds, and now this arc is an arc
-    ///of the matching.
-    template <typename MatchingMap>
-    void matchingMap(MatchingMap& map)  const {
-      typename Graph::template NodeMap<bool> todo(g,true);
-      for(NodeIt v(g); v!=INVALID; ++v) {
-        if (_mate[v]!=INVALID && v < _mate[v]) {
-          Node u=_mate[v];
-          for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
-            if ( g.runningNode(e) == u ) {
-              map.set(u,e);
-              map.set(v,e);
-              todo.set(u,false);
-              todo.set(v,false);
-              break;
-            }
-          }
-        }
-      }
-    }
-
-
-    ///\brief Gives back the matching in a \c bool valued \c Edge
-    ///map.
-    ///
-    ///Writes the matching stored to a \c bool valued \c Arc
-    ///map. This map will have the property that there are no two
-    ///incident arcs \c e, \c f with \c map[e]==map[f]==true. The
-    ///arcs \c e with \c map[e]==true form the matching.
-    template<typename MatchingMap>
-    void matching(MatchingMap& map) const {
-      for(EdgeIt e(g); e!=INVALID; ++e) map.set(e,false);
-
-      typename Graph::template NodeMap<bool> todo(g,true);
-      for(NodeIt v(g); v!=INVALID; ++v) {
-        if ( todo[v] && _mate[v]!=INVALID ) {
-          Node u=_mate[v];
-          for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
-            if ( g.runningNode(e) == u ) {
-              map.set(e,true);
-              todo.set(u,false);
-              todo.set(v,false);
-              break;
-            }
-          }
-        }
-      }
-    }
-
-
-    ///\brief Returns the canonical decomposition of the digraph after running
-    ///the algorithm.
-    ///
-    ///After calling any run methods of the class, it writes the
-    ///Gallai-Edmonds canonical decomposition of the digraph. \c map
-    ///must be a node map of \ref DecompType 's.
-    template <typename DecompositionMap>
-    void decomposition(DecompositionMap& map) const {
-      for(NodeIt v(g); v!=INVALID; ++v) map.set(v,position[v]);
-    }
-
-    ///\brief Returns a barrier on the nodes.
-    ///
-    ///After calling any run methods of the class, it writes a
-    ///canonical barrier on the nodes. The odd component number of the
-    ///remaining digraph minus the barrier size is a lower bound for the
-    ///uncovered nodes in the digraph. The \c map must be a node map of
-    ///bools.
-    template <typename BarrierMap>
-    void barrier(BarrierMap& barrier) {
-      for(NodeIt v(g); v!=INVALID; ++v) barrier.set(v,position[v] == A);
-    }
-
-  private:
-
-
-    void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear,
-                    UFE& blossom, NV& rep, EFE& tree) {
-      //We have one tree which we grow, and also shrink but only if it
-      //cannot be postponed. If we augment then we return to the "for"
-      //cycle of runEdmonds().
-
-      std::queue<Node> Q;   //queue of the totally unscanned nodes
-      Q.push(v);
-      std::queue<Node> R;
-      //queue of the nodes which must be scanned for a possible shrink
-
-      while ( !Q.empty() ) {
-        Node x=Q.front();
-        Q.pop();
-        for( IncEdgeIt e(g,x); e!= INVALID; ++e ) {
-          Node y=g.runningNode(e);
-          //growOrAugment grows if y is covered by the matching and
-          //augments if not. In this latter case it returns 1.
-          if (position[y]==C &&
-              growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
-        }
-        R.push(x);
-      }
-
-      while ( !R.empty() ) {
-        Node x=R.front();
-        R.pop();
-
-        for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) {
-          Node y=g.runningNode(e);
-
-          if ( position[y] == D && blossom.find(x) != blossom.find(y) )
-            //Recall that we have only one tree.
-            shrink( x, y, ear, blossom, rep, tree, Q);
-
-          while ( !Q.empty() ) {
-            Node z=Q.front();
-            Q.pop();
-            for( IncEdgeIt f(g,z); f!= INVALID; ++f ) {
-              Node w=g.runningNode(f);
-              //growOrAugment grows if y is covered by the matching and
-              //augments if not. In this latter case it returns 1.
-              if (position[w]==C &&
-                  growOrAugment(w, z, ear, blossom, rep, tree, Q)) return;
-            }
-            R.push(z);
-          }
-        } //for e
-      } // while ( !R.empty() )
-    }
-
-    void normShrink(Node v, typename Graph::template NodeMap<Node>& ear,
-                    UFE& blossom, NV& rep, EFE& tree) {
-      //We have one tree, which we grow and shrink. If we augment then we
-      //return to the "for" cycle of runEdmonds().
-
-      std::queue<Node> Q;   //queue of the unscanned nodes
-      Q.push(v);
-      while ( !Q.empty() ) {
-
-        Node x=Q.front();
-        Q.pop();
-
-        for( IncEdgeIt e(g,x); e!=INVALID; ++e ) {
-          Node y=g.runningNode(e);
-
-          switch ( position[y] ) {
-          case D:          //x and y must be in the same tree
-            if ( blossom.find(x) != blossom.find(y))
-              //x and y are in the same tree
-              shrink(x, y, ear, blossom, rep, tree, Q);
-            break;
-          case C:
-            //growOrAugment grows if y is covered by the matching and
-            //augments if not. In this latter case it returns 1.
-            if (growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
-            break;
-          default: break;
-          }
-        }
-      }
-    }
-
-    void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear,
-                UFE& blossom, NV& rep, EFE& tree,std::queue<Node>& Q) {
-      //x and y are the two adjacent vertices in two blossoms.
-
-      typename Graph::template NodeMap<bool> path(g,false);
-
-      Node b=rep[blossom.find(x)];
-      path.set(b,true);
-      b=_mate[b];
-      while ( b!=INVALID ) {
-        b=rep[blossom.find(ear[b])];
-        path.set(b,true);
-        b=_mate[b];
-      } //we go until the root through bases of blossoms and odd vertices
-
-      Node top=y;
-      Node middle=rep[blossom.find(top)];
-      Node bottom=x;
-      while ( !path[middle] )
-        shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
-      //Until we arrive to a node on the path, we update blossom, tree
-      //and the positions of the odd nodes.
-
-      Node base=middle;
-      top=x;
-      middle=rep[blossom.find(top)];
-      bottom=y;
-      Node blossom_base=rep[blossom.find(base)];
-      while ( middle!=blossom_base )
-        shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
-      //Until we arrive to a node on the path, we update blossom, tree
-      //and the positions of the odd nodes.
-
-      rep[blossom.find(base)] = base;
-    }
-
-    void shrinkStep(Node& top, Node& middle, Node& bottom,
-                    typename Graph::template NodeMap<Node>& ear,
-                    UFE& blossom, NV& rep, EFE& tree, std::queue<Node>& Q) {
-      //We traverse a blossom and update everything.
-
-      ear.set(top,bottom);
-      Node t=top;
-      while ( t!=middle ) {
-        Node u=_mate[t];
-        t=ear[u];
-        ear.set(t,u);
-      }
-      bottom=_mate[middle];
-      position.set(bottom,D);
-      Q.push(bottom);
-      top=ear[bottom];
-      Node oldmiddle=middle;
-      middle=rep[blossom.find(top)];
-      tree.erase(bottom);
-      tree.erase(oldmiddle);
-      blossom.insert(bottom);
-      blossom.join(bottom, oldmiddle);
-      blossom.join(top, oldmiddle);
-    }
-
-
-
-    bool growOrAugment(Node& y, Node& x, typename Graph::template
-                       NodeMap<Node>& ear, UFE& blossom, NV& rep, EFE& tree,
-                       std::queue<Node>& Q) {
-      //x is in a blossom in the tree, y is outside. If y is covered by
-      //the matching we grow, otherwise we augment. In this case we
-      //return 1.
-
-      if ( _mate[y]!=INVALID ) {       //grow
-        ear.set(y,x);
-        Node w=_mate[y];
-        rep[blossom.insert(w)] = w;
-        position.set(y,A);
-        position.set(w,D);
-        int t = tree.find(rep[blossom.find(x)]);
-        tree.insert(y,t);
-        tree.insert(w,t);
-        Q.push(w);
-      } else {                      //augment
-        augment(x, ear, blossom, rep, tree);
-        _mate.set(x,y);
-        _mate.set(y,x);
-        return true;
-      }
-      return false;
-    }
-
-    void augment(Node x, typename Graph::template NodeMap<Node>& ear,
-                 UFE& blossom, NV& rep, EFE& tree) {
-      Node v=_mate[x];
-      while ( v!=INVALID ) {
-
-        Node u=ear[v];
-        _mate.set(v,u);
-        Node tmp=v;
-        v=_mate[u];
-        _mate.set(u,tmp);
-      }
-      int y = tree.find(rep[blossom.find(x)]);
-      for (typename EFE::ItemIt tit(tree, y); tit != INVALID; ++tit) {
-        if ( position[tit] == D ) {
-          int b = blossom.find(tit);
-          for (typename UFE::ItemIt bit(blossom, b); bit != INVALID; ++bit) {
-            position.set(bit, C);
-          }
-          blossom.eraseClass(b);
-        } else position.set(tit, C);
-      }
-      tree.eraseClass(y);
-
-    }
+    /// @}
 
   };
 
@@ -627 +617 @@
   /// \f$O(nm\log(n))\f$ time complexity.
   ///
   /// The maximum weighted matching problem is to find undirected
-  /// arcs in the digraph with maximum overall weight and no two of
-  /// them shares their endpoints. The problem can be formulated with
-  /// the next linear program:
+  /// edges in the graph with maximum overall weight and no two of
+  /// them shares their ends. The problem can be formulated with the
+  /// following linear program.
   /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
-  ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
+  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
+      \quad \forall B\in\mathcal{O}\f] */
   /// \f[x_e \ge 0\quad \forall e\in E\f]
   /// \f[\max \sum_{e\in E}x_ew_e\f]
-  /// where \f$\delta(X)\f$ is the set of arcs incident to a node in
-  /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in
-  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
-  /// the nodes.
+  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
+  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
+  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
+  /// subsets of the nodes.
   ///
   /// The algorithm calculates an optimal matching and a proof of the
   /// optimality. The solution of the dual problem can be used to check
-  /// the result of the algorithm. The dual linear problem is the next:
-  /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f]
+  /// the result of the algorithm. The dual linear problem is the
+  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
+      z_B \ge w_{uv} \quad \forall uv\in E\f] */
   /// \f[y_u \ge 0 \quad \forall u \in V\f]
   /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
-  /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
+  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
+      \frac{\vert B \vert - 1}{2}z_B\f] */
   ///
   /// The algorithm can be executed with \c run() or the \c init() and
   /// then the \c start() member functions. After it the matching can
@@ -705 +698 @@
     int _node_num;
     int _blossom_num;
 
-    typedef typename Graph::template NodeMap<int> NodeIntMap;
-    typedef typename Graph::template ArcMap<int> ArcIntMap;
-    typedef typename Graph::template EdgeMap<int> EdgeIntMap;
     typedef RangeMap<int> IntIntMap;
 
     enum Status {
       EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
     };
 
-    typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
+    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
     struct BlossomData {
       int tree;
       Status status;
@@ -723 +713 @@
       Node base;
     };
 
-    NodeIntMap *_blossom_index;
+    IntNodeMap *_blossom_index;
     BlossomSet *_blossom_set;
     RangeMap<BlossomData>* _blossom_data;
 
-    NodeIntMap *_node_index;
-    ArcIntMap *_node_heap_index;
+    IntNodeMap *_node_index;
+    IntArcMap *_node_heap_index;
 
     struct NodeData {
 
-      NodeData(ArcIntMap& node_heap_index)
+      NodeData(IntArcMap& node_heap_index)
         : heap(node_heap_index) {}
 
       int blossom;
       Value pot;
-      BinHeap<Value, ArcIntMap> heap;
+      BinHeap<Value, IntArcMap> heap;
       std::map<int, Arc> heap_index;
 
       int tree;
@@ -750 +740 @@
     IntIntMap *_tree_set_index;
     TreeSet *_tree_set;
 
-    NodeIntMap *_delta1_index;
-    BinHeap<Value, NodeIntMap> *_delta1;
+    IntNodeMap *_delta1_index;
+    BinHeap<Value, IntNodeMap> *_delta1;
 
     IntIntMap *_delta2_index;
     BinHeap<Value, IntIntMap> *_delta2;
 
-    EdgeIntMap *_delta3_index;
-    BinHeap<Value, EdgeIntMap> *_delta3;
+    IntEdgeMap *_delta3_index;
+    BinHeap<Value, IntEdgeMap> *_delta3;
 
     IntIntMap *_delta4_index;
     BinHeap<Value, IntIntMap> *_delta4;
@@ -775 +765 @@
         _node_potential = new NodePotential(_graph);
       }
       if (!_blossom_set) {
-        _blossom_index = new NodeIntMap(_graph);
+        _blossom_index = new IntNodeMap(_graph);
         _blossom_set = new BlossomSet(*_blossom_index);
         _blossom_data = new RangeMap<BlossomData>(_blossom_num);
       }
 
       if (!_node_index) {
-        _node_index = new NodeIntMap(_graph);
-        _node_heap_index = new ArcIntMap(_graph);
+        _node_index = new IntNodeMap(_graph);
+        _node_heap_index = new IntArcMap(_graph);
         _node_data = new RangeMap<NodeData>(_node_num,
                                               NodeData(*_node_heap_index));
       }
@@ -792 +782 @@
         _tree_set = new TreeSet(*_tree_set_index);
       }
       if (!_delta1) {
-        _delta1_index = new NodeIntMap(_graph);
-        _delta1 = new BinHeap<Value, NodeIntMap>(*_delta1_index);
+        _delta1_index = new IntNodeMap(_graph);
+        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
       }
       if (!_delta2) {
         _delta2_index = new IntIntMap(_blossom_num);
         _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
       }
       if (!_delta3) {
-        _delta3_index = new EdgeIntMap(_graph);
-        _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index);
+        _delta3_index = new IntEdgeMap(_graph);
+        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
       }
       if (!_delta4) {
         _delta4_index = new IntIntMap(_blossom_num);
@@ -1266 +1256 @@
     }
 
 
-    void augmentOnArc(const Edge& arc) {
+    void augmentOnEdge(const Edge& edge) {
 
-      int left = _blossom_set->find(_graph.u(arc));
-      int right = _blossom_set->find(_graph.v(arc));
+      int left = _blossom_set->find(_graph.u(edge));
+      int right = _blossom_set->find(_graph.v(edge));
 
       if ((*_blossom_data)[left].status == EVEN) {
         int left_tree = _tree_set->find(left);
@@ -1289 +1279 @@
         unmatchedToMatched(right);
       }
 
-      (*_blossom_data)[left].next = _graph.direct(arc, true);
-      (*_blossom_data)[right].next = _graph.direct(arc, false);
+      (*_blossom_data)[left].next = _graph.direct(edge, true);
+      (*_blossom_data)[right].next = _graph.direct(edge, false);
     }
 
     void extendOnArc(const Arc& arc) {
@@ -1310 +1300 @@
       matchedToEven(even, tree);
     }
 
-    void shrinkOnArc(const Edge& edge, int tree) {
+    void shrinkOnEdge(const Edge& edge, int tree) {
       int nca = -1;
       std::vector<int> left_path, right_path;
 
@@ -1652 +1642 @@
       createStructures();
 
       for (ArcIt e(_graph); e != INVALID; ++e) {
-        _node_heap_index->set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP);
+        _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
       }
       for (NodeIt n(_graph); n != INVALID; ++n) {
         _delta1_index->set(n, _delta1->PRE_HEAP);
@@ -1769 +1759 @@
               }
 
               if (left_tree == right_tree) {
-                shrinkOnArc(e, left_tree);
+                shrinkOnEdge(e, left_tree);
               } else {
-                augmentOnArc(e);
+                augmentOnEdge(e);
                 unmatched -= 2;
               }
             }
@@ -1818 +1808 @@
       return sum /= 2;
     }
 
-    /// \brief Returns true when the arc is in the matching.
+    /// \brief Returns the cardinality of the matching.
     ///
-    /// Returns true when the arc is in the matching.
-    bool matching(const Edge& arc) const {
-      return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true);
+    /// Returns the cardinality of the matching.
+    int matchingSize() const {
+      int num = 0;
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        if ((*_matching)[n] != INVALID) {
+          ++num;
+        }
+      }
+      return num /= 2;
+    }
+
+    /// \brief Returns true when the edge is in the matching.
+    ///
+    /// Returns true when the edge is in the matching.
+    bool matching(const Edge& edge) const {
+      return edge == (*_matching)[_graph.u(edge)];
     }
 
     /// \brief Returns the incident matching arc.
@@ -1913 +1916 @@
         _last = _algorithm->_blossom_potential[variable].end;
       }
 
-      /// \brief Invalid constructor.
-      ///
-      /// Invalid constructor.
-      BlossomIt(Invalid) : _index(-1) {}
-
       /// \brief Conversion to node.
       ///
       /// Conversion to node.
       operator Node() const {
-        return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
+        return _algorithm->_blossom_node_list[_index];
       }
 
       /// \brief Increment operator.
@@ -1930 +1928 @@
       /// Increment operator.
       BlossomIt& operator++() {
         ++_index;
-        if (_index == _last) {
-          _index = -1;
-        }
         return *this;
       }
 
-      bool operator==(const BlossomIt& it) const {
-        return _index == it._index;
-      }
-      bool operator!=(const BlossomIt& it) const {
-        return _index != it._index;
-      }
+      /// \brief Validity checking
+      ///
+      /// Checks whether the iterator is invalid.
+      bool operator==(Invalid) const { return _index == _last; }
+
+      /// \brief Validity checking
+      ///
+      /// Checks whether the iterator is valid.
+      bool operator!=(Invalid) const { return _index != _last; }
 
     private:
       const MaxWeightedMatching* _algorithm;
@@ -1958 +1956 @@
   /// \brief Weighted perfect matching in general graphs
   ///
   /// This class provides an efficient implementation of Edmond's
-  /// maximum weighted perfecr matching algorithm. The implementation
+  /// maximum weighted perfect matching algorithm. The implementation
   /// is based on extensive use of priority queues and provides
   /// \f$O(nm\log(n))\f$ time complexity.
   ///
   /// The maximum weighted matching problem is to find undirected
-  /// arcs in the digraph with maximum overall weight and no two of
-  /// them shares their endpoints and covers all nodes. The problem
-  /// can be formulated with the next linear program:
+  /// edges in the graph with maximum overall weight and no two of
+  /// them shares their ends and covers all nodes. The problem can be
+  /// formulated with the following linear program.
   /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
-  ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
+  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
+      \quad \forall B\in\mathcal{O}\f] */
   /// \f[x_e \ge 0\quad \forall e\in E\f]
   /// \f[\max \sum_{e\in E}x_ew_e\f]
-  /// where \f$\delta(X)\f$ is the set of arcs incident to a node in
-  /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in
-  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
-  /// the nodes.
+  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
+  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
+  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
+  /// subsets of the nodes.
   ///
   /// The algorithm calculates an optimal matching and a proof of the
   /// optimality. The solution of the dual problem can be used to check
-  /// the result of the algorithm. The dual linear problem is the next:
-  /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f]
+  /// the result of the algorithm. The dual linear problem is the
+  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
+      w_{uv} \quad \forall uv\in E\f] */
   /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
-  /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
+  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
+      \frac{\vert B \vert - 1}{2}z_B\f] */
   ///
   /// The algorithm can be executed with \c run() or the \c init() and
   /// then the \c start() member functions. After it the matching can
@@ -2040 +2041 @@
     int _node_num;
     int _blossom_num;
 
-    typedef typename Graph::template NodeMap<int> NodeIntMap;
-    typedef typename Graph::template ArcMap<int> ArcIntMap;
-    typedef typename Graph::template EdgeMap<int> EdgeIntMap;
     typedef RangeMap<int> IntIntMap;
 
     enum Status {
       EVEN = -1, MATCHED = 0, ODD = 1
     };
 
-    typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
+    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
     struct BlossomData {
       int tree;
       Status status;
@@ -2057 +2055 @@
       Value pot, offset;
     };
 
-    NodeIntMap *_blossom_index;
+    IntNodeMap *_blossom_index;
     BlossomSet *_blossom_set;
     RangeMap<BlossomData>* _blossom_data;
 
-    NodeIntMap *_node_index;
-    ArcIntMap *_node_heap_index;
+    IntNodeMap *_node_index;
+    IntArcMap *_node_heap_index;
 
     struct NodeData {
 
-      NodeData(ArcIntMap& node_heap_index)
+      NodeData(IntArcMap& node_heap_index)
         : heap(node_heap_index) {}
 
       int blossom;
       Value pot;
-      BinHeap<Value, ArcIntMap> heap;
+      BinHeap<Value, IntArcMap> heap;
       std::map<int, Arc> heap_index;
 
       int tree;
@@ -2087 +2085 @@
     IntIntMap *_delta2_index;
     BinHeap<Value, IntIntMap> *_delta2;
 
-    EdgeIntMap *_delta3_index;
-    BinHeap<Value, EdgeIntMap> *_delta3;
+    IntEdgeMap *_delta3_index;
+    BinHeap<Value, IntEdgeMap> *_delta3;
 
     IntIntMap *_delta4_index;
     BinHeap<Value, IntIntMap> *_delta4;
@@ -2106 +2104 @@
         _node_potential = new NodePotential(_graph);
       }
       if (!_blossom_set) {
-        _blossom_index = new NodeIntMap(_graph);
+        _blossom_index = new IntNodeMap(_graph);
         _blossom_set = new BlossomSet(*_blossom_index);
         _blossom_data = new RangeMap<BlossomData>(_blossom_num);
       }
 
       if (!_node_index) {
-        _node_index = new NodeIntMap(_graph);
-        _node_heap_index = new ArcIntMap(_graph);
+        _node_index = new IntNodeMap(_graph);
+        _node_heap_index = new IntArcMap(_graph);
         _node_data = new RangeMap<NodeData>(_node_num,
-                                              NodeData(*_node_heap_index));
+                                            NodeData(*_node_heap_index));
       }
 
       if (!_tree_set) {
@@ -2127 +2125 @@
         _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
       }
       if (!_delta3) {
-        _delta3_index = new EdgeIntMap(_graph);
-        _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index);
+        _delta3_index = new IntEdgeMap(_graph);
+        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
       }
       if (!_delta4) {
         _delta4_index = new IntIntMap(_blossom_num);
@@ -2461 +2459 @@
       _tree_set->eraseClass(tree);
     }
 
-    void augmentOnArc(const Edge& arc) {
+    void augmentOnEdge(const Edge& edge) {
 
-      int left = _blossom_set->find(_graph.u(arc));
-      int right = _blossom_set->find(_graph.v(arc));
+      int left = _blossom_set->find(_graph.u(edge));
+      int right = _blossom_set->find(_graph.v(edge));
 
       int left_tree = _tree_set->find(left);
       alternatePath(left, left_tree);
@@ -2474 +2472 @@
       alternatePath(right, right_tree);
       destroyTree(right_tree);
 
-      (*_blossom_data)[left].next = _graph.direct(arc, true);
-      (*_blossom_data)[right].next = _graph.direct(arc, false);
+      (*_blossom_data)[left].next = _graph.direct(edge, true);
+      (*_blossom_data)[right].next = _graph.direct(edge, false);
     }
 
     void extendOnArc(const Arc& arc) {
@@ -2495 +2493 @@
       matchedToEven(even, tree);
     }
 
-    void shrinkOnArc(const Edge& edge, int tree) {
+    void shrinkOnEdge(const Edge& edge, int tree) {
       int nca = -1;
       std::vector<int> left_path, right_path;
 
@@ -2831 +2829 @@
       createStructures();
 
       for (ArcIt e(_graph); e != INVALID; ++e) {
-        _node_heap_index->set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP);
+        _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
       }
       for (EdgeIt e(_graph); e != INVALID; ++e) {
         _delta3_index->set(e, _delta3->PRE_HEAP);
@@ -2924 +2922 @@
               int right_tree = _tree_set->find(right_blossom);
 
               if (left_tree == right_tree) {
-                shrinkOnArc(e, left_tree);
+                shrinkOnEdge(e, left_tree);
               } else {
-                augmentOnArc(e);
+                augmentOnEdge(e);
                 unmatched -= 2;
               }
             }
@@ -2974 +2972 @@
       return sum /= 2;
     }
 
-    /// \brief Returns true when the arc is in the matching.
+    /// \brief Returns true when the edge is in the matching.
     ///
-    /// Returns true when the arc is in the matching.
-    bool matching(const Edge& arc) const {
-      return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true);
+    /// Returns true when the edge is in the matching.
+    bool matching(const Edge& edge) const {
+      return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
     }
 
-    /// \brief Returns the incident matching arc.
+    /// \brief Returns the incident matching edge.
     ///
-    /// Returns the incident matching arc from given node.
+    /// Returns the incident matching arc from given edge.
     Arc matching(const Node& node) const {
       return (*_matching)[node];
     }
@@ -3066 +3064 @@
         _last = _algorithm->_blossom_potential[variable].end;
       }
 
-      /// \brief Invalid constructor.
-      ///
-      /// Invalid constructor.
-      BlossomIt(Invalid) : _index(-1) {}
-
       /// \brief Conversion to node.
       ///
       /// Conversion to node.
       operator Node() const {
-        return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
+        return _algorithm->_blossom_node_list[_index];
       }
 
       /// \brief Increment operator.
@@ -3083 +3076 @@
       /// Increment operator.
       BlossomIt& operator++() {
         ++_index;
-        if (_index == _last) {
-          _index = -1;
-        }
         return *this;
       }
 
-      bool operator==(const BlossomIt& it) const {
-        return _index == it._index;
-      }
-      bool operator!=(const BlossomIt& it) const {
-        return _index != it._index;
-      }
+      /// \brief Validity checking
+      ///
+      /// Checks whether the iterator is invalid.
+      bool operator==(Invalid) const { return _index == _last; }
+
+      /// \brief Validity checking
+      ///
+      /// Checks whether the iterator is valid.
+      bool operator!=(Invalid) const { return _index != _last; }
 
     private:
       const MaxWeightedPerfectMatching* _algorithm;

test/CMakeLists.txt

@@ -17 +17 @@
   kruskal_test
   maps_test
   max_matching_test
-  max_weighted_matching_test
   random_test
   path_test
   time_measure_test

test/Makefile.am

@@ -20 +20 @@
         test/kruskal_test \
         test/maps_test \
         test/max_matching_test \
-        test/max_weighted_matching_test \
         test/random_test \
         test/path_test \
         test/test_tools_fail \
@@ -45 +44 @@
 test_kruskal_test_SOURCES = test/kruskal_test.cc
 test_maps_test_SOURCES = test/maps_test.cc
 test_max_matching_test_SOURCES = test/max_matching_test.cc
-test_max_weighted_matching_test_SOURCES = test/max_weighted_matching_test.cc
 test_path_test_SOURCES = test/path_test.cc
 test_random_test_SOURCES = test/random_test.cc
 test_test_tools_fail_SOURCES = test/test_tools_fail.cc

test/max_matching_test.cc

@@ -17 +17 @@
  */
 
 #include <iostream>
+#include <sstream>
 #include <vector>
 #include <queue>
 #include <lemon/math.h>
 #include <cstdlib>
 
+#include <lemon/max_matching.h>
+#include <lemon/smart_graph.h>
+#include <lemon/lgf_reader.h>
+
 #include "test_tools.h"
-#include <lemon/list_graph.h>
-#include <lemon/max_matching.h>
 
 using namespace std;
 using namespace lemon;
 
+GRAPH_TYPEDEFS(SmartGraph);
+
+
+const int lgfn = 3;
+const std::string lgf[lgfn] = {
+  "@nodes\n"
+  "label\n"
+  "0\n"
+  "1\n"
+  "2\n"
+  "3\n"
+  "4\n"
+  "5\n"
+  "6\n"
+  "7\n"
+  "@edges\n"
+  "     label  weight\n"
+  "7 4  0      984\n"
+  "0 7  1      73\n"
+  "7 1  2      204\n"
+  "2 3  3      583\n"
+  "2 7  4      565\n"
+  "2 1  5      582\n"
+  "0 4  6      551\n"
+  "2 5  7      385\n"
+  "1 5  8      561\n"
+  "5 3  9      484\n"
+  "7 5  10     904\n"
+  "3 6  11     47\n"
+  "7 6  12     888\n"
+  "3 0  13     747\n"
+  "6 1  14     310\n",
+
+  "@nodes\n"
+  "label\n"
+  "0\n"
+  "1\n"
+  "2\n"
+  "3\n"
+  "4\n"
+  "5\n"
+  "6\n"
+  "7\n"
+  "@edges\n"
+  "     label  weight\n"
+  "2 5  0      710\n"
+  "0 5  1      241\n"
+  "2 4  2      856\n"
+  "2 6  3      762\n"
+  "4 1  4      747\n"
+  "6 1  5      962\n"
+  "4 7  6      723\n"
+  "1 7  7      661\n"
+  "2 3  8      376\n"
+  "1 0  9      416\n"
+  "6 7  10     391\n",
+
+  "@nodes\n"
+  "label\n"
+  "0\n"
+  "1\n"
+  "2\n"
+  "3\n"
+  "4\n"
+  "5\n"
+  "6\n"
+  "7\n"
+  "@edges\n"
+  "     label  weight\n"
+  "6 2  0      553\n"
+  "0 7  1      653\n"
+  "6 3  2      22\n"
+  "4 7  3      846\n"
+  "7 2  4      981\n"
+  "7 6  5      250\n"
+  "5 2  6      539\n",
+};
+
+void checkMatching(const SmartGraph& graph,
+                   const MaxMatching<SmartGraph>& mm) {
+  int num = 0;
+
+  IntNodeMap comp_index(graph);
+  UnionFind<IntNodeMap> comp(comp_index);
+
+  int barrier_num = 0;
+
+  for (NodeIt n(graph); n != INVALID; ++n) {
+    check(mm.decomposition(n) == MaxMatching<SmartGraph>::EVEN ||
+          mm.matching(n) != INVALID, "Wrong Gallai-Edmonds decomposition");
+    if (mm.decomposition(n) == MaxMatching<SmartGraph>::ODD) {
+      ++barrier_num;
+    } else {
+      comp.insert(n);
+    }
+  }
+
+  for (EdgeIt e(graph); e != INVALID; ++e) {
+    if (mm.matching(e)) {
+      check(e == mm.matching(graph.u(e)), "Wrong matching");
+      check(e == mm.matching(graph.v(e)), "Wrong matching");
+      ++num;
+    }
+    check(mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::EVEN ||
+          mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::MATCHED,
+          "Wrong Gallai-Edmonds decomposition");
+
+    check(mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::EVEN ||
+          mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::MATCHED,
+          "Wrong Gallai-Edmonds decomposition");
+
+    if (mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::ODD &&
+        mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::ODD) {
+      comp.join(graph.u(e), graph.v(e));
+    }
+  }
+
+  std::set<int> comp_root;
+  int odd_comp_num = 0;
+  for (NodeIt n(graph); n != INVALID; ++n) {
+    if (mm.decomposition(n) != MaxMatching<SmartGraph>::ODD) {
+      int root = comp.find(n);
+      if (comp_root.find(root) == comp_root.end()) {
+        comp_root.insert(root);
+        if (comp.size(n) % 2 == 1) {
+          ++odd_comp_num;
+        }
+      }
+    }
+  }
+
+  check(mm.matchingSize() == num, "Wrong matching");
+  check(2 * num == countNodes(graph) - (odd_comp_num - barrier_num),
+         "Wrong matching");
+  return;
+}
+
+void checkWeightedMatching(const SmartGraph& graph,
+                   const SmartGraph::EdgeMap<int>& weight,
+                   const MaxWeightedMatching<SmartGraph>& mwm) {
+  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
+    if (graph.u(e) == graph.v(e)) continue;
+    int rw = mwm.nodeValue(graph.u(e)) + mwm.nodeValue(graph.v(e));
+
+    for (int i = 0; i < mwm.blossomNum(); ++i) {
+      bool s = false, t = false;
+      for (MaxWeightedMatching<SmartGraph>::BlossomIt n(mwm, i);
+           n != INVALID; ++n) {
+        if (graph.u(e) == n) s = true;
+        if (graph.v(e) == n) t = true;
+      }
+      if (s == true && t == true) {
+        rw += mwm.blossomValue(i);
+      }
+    }
+    rw -= weight[e] * mwm.dualScale;
+
+    check(rw >= 0, "Negative reduced weight");
+    check(rw == 0 || !mwm.matching(e),
+          "Non-zero reduced weight on matching edge");
+  }
+
+  int pv = 0;
+  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
+    if (mwm.matching(n) != INVALID) {
+      check(mwm.nodeValue(n) >= 0, "Invalid node value");
+      pv += weight[mwm.matching(n)];
+      SmartGraph::Node o = graph.target(mwm.matching(n));
+      check(mwm.mate(n) == o, "Invalid matching");
+      check(mwm.matching(n) == graph.oppositeArc(mwm.matching(o)),
+            "Invalid matching");
+    } else {
+      check(mwm.mate(n) == INVALID, "Invalid matching");
+      check(mwm.nodeValue(n) == 0, "Invalid matching");
+    }
+  }
+
+  int dv = 0;
+  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
+    dv += mwm.nodeValue(n);
+  }
+
+  for (int i = 0; i < mwm.blossomNum(); ++i) {
+    check(mwm.blossomValue(i) >= 0, "Invalid blossom value");
+    check(mwm.blossomSize(i) % 2 == 1, "Even blossom size");
+    dv += mwm.blossomValue(i) * ((mwm.blossomSize(i) - 1) / 2);
+  }
+
+  check(pv * mwm.dualScale == dv * 2, "Wrong duality");
+
+  return;
+}
+
+void checkWeightedPerfectMatching(const SmartGraph& graph,
+                          const SmartGraph::EdgeMap<int>& weight,
+                          const MaxWeightedPerfectMatching<SmartGraph>& mwpm) {
+  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
+    if (graph.u(e) == graph.v(e)) continue;
+    int rw = mwpm.nodeValue(graph.u(e)) + mwpm.nodeValue(graph.v(e));
+
+    for (int i = 0; i < mwpm.blossomNum(); ++i) {
+      bool s = false, t = false;
+      for (MaxWeightedPerfectMatching<SmartGraph>::BlossomIt n(mwpm, i);
+           n != INVALID; ++n) {
+        if (graph.u(e) == n) s = true;
+        if (graph.v(e) == n) t = true;
+      }
+      if (s == true && t == true) {
+        rw += mwpm.blossomValue(i);
+      }
+    }
+    rw -= weight[e] * mwpm.dualScale;
+
+    check(rw >= 0, "Negative reduced weight");
+    check(rw == 0 || !mwpm.matching(e),
+          "Non-zero reduced weight on matching edge");
+  }
+
+  int pv = 0;
+  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
+    check(mwpm.matching(n) != INVALID, "Non perfect");
+    pv += weight[mwpm.matching(n)];
+    SmartGraph::Node o = graph.target(mwpm.matching(n));
+    check(mwpm.mate(n) == o, "Invalid matching");
+    check(mwpm.matching(n) == graph.oppositeArc(mwpm.matching(o)),
+          "Invalid matching");
+  }
+
+  int dv = 0;
+  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
+    dv += mwpm.nodeValue(n);
+  }
+
+  for (int i = 0; i < mwpm.blossomNum(); ++i) {
+    check(mwpm.blossomValue(i) >= 0, "Invalid blossom value");
+    check(mwpm.blossomSize(i) % 2 == 1, "Even blossom size");
+    dv += mwpm.blossomValue(i) * ((mwpm.blossomSize(i) - 1) / 2);
+  }
+
+  check(pv * mwpm.dualScale == dv * 2, "Wrong duality");
+
+  return;
+}
+
+
 int main() {
 
-  typedef ListGraph Graph;
+  for (int i = 0; i < lgfn; ++i) {
+    SmartGraph graph;
+    SmartGraph::EdgeMap<int> weight(graph);
 
-  GRAPH_TYPEDEFS(Graph);
+    istringstream lgfs(lgf[i]);
+    graphReader(graph, lgfs).
+      edgeMap("weight", weight).run();
 
-  Graph g;
-  g.clear();
+    MaxMatching<SmartGraph> mm(graph);
+    mm.run();
+    checkMatching(graph, mm);
 
-  std::vector<Graph::Node> nodes;
-  for (int i=0; i<13; ++i)
-      nodes.push_back(g.addNode());
+    MaxWeightedMatching<SmartGraph> mwm(graph, weight);
+    mwm.run();
+    checkWeightedMatching(graph, weight, mwm);
 
-  g.addEdge(nodes[0], nodes[0]);
-  g.addEdge(nodes[6], nodes[10]);
-  g.addEdge(nodes[5], nodes[10]);
-  g.addEdge(nodes[4], nodes[10]);
-  g.addEdge(nodes[3], nodes[11]);
-  g.addEdge(nodes[1], nodes[6]);
-  g.addEdge(nodes[4], nodes[7]);
-  g.addEdge(nodes[1], nodes[8]);
-  g.addEdge(nodes[0], nodes[8]);
-  g.addEdge(nodes[3], nodes[12]);
-  g.addEdge(nodes[6], nodes[9]);
-  g.addEdge(nodes[9], nodes[11]);
-  g.addEdge(nodes[2], nodes[10]);
-  g.addEdge(nodes[10], nodes[8]);
-  g.addEdge(nodes[5], nodes[8]);
-  g.addEdge(nodes[6], nodes[3]);
-  g.addEdge(nodes[0], nodes[5]);
-  g.addEdge(nodes[6], nodes[12]);
+    MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
+    bool perfect = mwpm.run();
 
-  MaxMatching<Graph> max_matching(g);
-  max_matching.init();
-  max_matching.startDense();
+    check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
+          "Perfect matching found");
 
-  int s=0;
-  Graph::NodeMap<Node> mate(g,INVALID);
-  max_matching.mateMap(mate);
-  for(NodeIt v(g); v!=INVALID; ++v) {
-    if ( mate[v]!=INVALID ) ++s;
-  }
-  int size=int(s/2);  //size will be used as the size of a maxmatching
-
-  for(NodeIt v(g); v!=INVALID; ++v) {
-    max_matching.mate(v);
-  }
-
-  check ( size == max_matching.size(), "mate() returns a different size matching than max_matching.size()" );
-
-  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos0(g);
-  max_matching.decomposition(pos0);
-
-  max_matching.init();
-  max_matching.startSparse();
-  s=0;
-  max_matching.mateMap(mate);
-  for(NodeIt v(g); v!=INVALID; ++v) {
-    if ( mate[v]!=INVALID ) ++s;
-  }
-  check ( int(s/2) == size, "The size does not equal!" );
-
-  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos1(g);
-  max_matching.decomposition(pos1);
-
-  max_matching.run();
-  s=0;
-  max_matching.mateMap(mate);
-  for(NodeIt v(g); v!=INVALID; ++v) {
-    if ( mate[v]!=INVALID ) ++s;
-  }
-  check ( int(s/2) == size, "The size does not equal!" );
-
-  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos2(g);
-  max_matching.decomposition(pos2);
-
-  max_matching.run();
-  s=0;
-  max_matching.mateMap(mate);
-  for(NodeIt v(g); v!=INVALID; ++v) {
-    if ( mate[v]!=INVALID ) ++s;
-  }
-  check ( int(s/2) == size, "The size does not equal!" );
-
-  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos(g);
-  max_matching.decomposition(pos);
-
-  bool ismatching=true;
-  for(NodeIt v(g); v!=INVALID; ++v) {
-    if ( mate[v]!=INVALID ) {
-      Node u=mate[v];
-      if (mate[u]!=v) ismatching=false;
+    if (perfect) {
+      checkWeightedPerfectMatching(graph, weight, mwpm);
     }
   }
-  check ( ismatching, "It is not a matching!" );
-
-  bool coincide=true;
-  for(NodeIt v(g); v!=INVALID; ++v) {
-   if ( pos0[v] != pos1[v] || pos1[v]!=pos2[v] || pos2[v]!=pos[v] ) {
-     coincide=false;
-    }
-  }
-  check ( coincide, "The decompositions do not coincide! " );
-
-  bool noarc=true;
-  for(EdgeIt e(g); e!=INVALID; ++e) {
-   if ( (pos[g.v(e)]==max_matching.C &&
-         pos[g.u(e)]==max_matching.D) ||
-         (pos[g.v(e)]==max_matching.D &&
-          pos[g.u(e)]==max_matching.C) )
-      noarc=false;
-  }
-  check ( noarc, "There are arcs between D and C!" );
-
-  bool oddcomp=true;
-  Graph::NodeMap<bool> todo(g,true);
-  int num_comp=0;
-  for(NodeIt v(g); v!=INVALID; ++v) {
-   if ( pos[v]==max_matching.D && todo[v] ) {
-      int comp_size=1;
-      ++num_comp;
-      std::queue<Node> Q;
-      Q.push(v);
-      todo.set(v,false);
-      while (!Q.empty()) {
-        Node w=Q.front();
-        Q.pop();
-        for(IncEdgeIt e(g,w); e!=INVALID; ++e) {
-          Node u=g.runningNode(e);
-          if ( pos[u]==max_matching.D && todo[u] ) {
-            ++comp_size;
-            Q.push(u);
-            todo.set(u,false);
-          }
-        }
-      }
-      if ( !(comp_size % 2) ) oddcomp=false;
-    }
-  }
-  check ( oddcomp, "A component of g[D] is not odd." );
-
-  int barrier=0;
-  for(NodeIt v(g); v!=INVALID; ++v) {
-    if ( pos[v]==max_matching.A ) ++barrier;
-  }
-  int expected_size=int( countNodes(g)-num_comp+barrier)/2;
-  check ( size==expected_size, "The size of the matching is wrong." );
 
   return 0;
 }

deleted file test/max_weighted_matching_test.cc

@@ -1 @@
-/* -*- mode: C++; indent-tabs-mode: nil; -*-
- *
- * This file is a part of LEMON, a generic C++ optimization library.
- *
- * Copyright (C) 2003-2008
- * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
- * (Egervary Research Group on Combinatorial Optimization, EGRES).
- *
- * Permission to use, modify and distribute this software is granted
- * provided that this copyright notice appears in all copies. For
- * precise terms see the accompanying LICENSE file.
- *
- * This software is provided "AS IS" with no warranty of any kind,
- * express or implied, and with no claim as to its suitability for any
- * purpose.
- *
- */
-
-#include <iostream>
-#include <sstream>
-#include <vector>
-#include <queue>
-#include <lemon/math.h>
-#include <cstdlib>
-
-#include "test_tools.h"
-#include <lemon/smart_graph.h>
-#include <lemon/max_matching.h>
-#include <lemon/lgf_reader.h>
-
-using namespace std;
-using namespace lemon;
-
-GRAPH_TYPEDEFS(SmartGraph);
-
-const int lgfn = 3;
-const std::string lgf[lgfn] = {
-  "@nodes\n"
-  "label\n"
-  "0\n"
-  "1\n"
-  "2\n"
-  "3\n"
-  "4\n"
-  "5\n"
-  "6\n"
-  "7\n"
-  "@edges\n"
-  "label        weight\n"
-  "7    4       0       984\n"
-  "0    7       1       73\n"
-  "7    1       2       204\n"
-  "2    3       3       583\n"
-  "2    7       4       565\n"
-  "2    1       5       582\n"
-  "0    4       6       551\n"
-  "2    5       7       385\n"
-  "1    5       8       561\n"
-  "5    3       9       484\n"
-  "7    5       10      904\n"
-  "3    6       11      47\n"
-  "7    6       12      888\n"
-  "3    0       13      747\n"
-  "6    1       14      310\n",
-
-  "@nodes\n"
-  "label\n"
-  "0\n"
-  "1\n"
-  "2\n"
-  "3\n"
-  "4\n"
-  "5\n"
-  "6\n"
-  "7\n"
-  "@edges\n"
-  "label        weight\n"
-  "2    5       0       710\n"
-  "0    5       1       241\n"
-  "2    4       2       856\n"
-  "2    6       3       762\n"
-  "4    1       4       747\n"
-  "6    1       5       962\n"
-  "4    7       6       723\n"
-  "1    7       7       661\n"
-  "2    3       8       376\n"
-  "1    0       9       416\n"
-  "6    7       10      391\n",
-
-  "@nodes\n"
-  "label\n"
-  "0\n"
-  "1\n"
-  "2\n"
-  "3\n"
-  "4\n"
-  "5\n"
-  "6\n"
-  "7\n"
-  "@edges\n"
-  "label        weight\n"
-  "6    2       0       553\n"
-  "0    7       1       653\n"
-  "6    3       2       22\n"
-  "4    7       3       846\n"
-  "7    2       4       981\n"
-  "7    6       5       250\n"
-  "5    2       6       539\n"
-};
-
-void checkMatching(const SmartGraph& graph,
-                   const SmartGraph::EdgeMap<int>& weight,
-                   const MaxWeightedMatching<SmartGraph>& mwm) {
-  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
-    if (graph.u(e) == graph.v(e)) continue;
-    int rw =
-      mwm.nodeValue(graph.u(e)) + mwm.nodeValue(graph.v(e));
-
-    for (int i = 0; i < mwm.blossomNum(); ++i) {
-      bool s = false, t = false;
-      for (MaxWeightedMatching<SmartGraph>::BlossomIt n(mwm, i);
-           n != INVALID; ++n) {
-        if (graph.u(e) == n) s = true;
-        if (graph.v(e) == n) t = true;
-      }
-      if (s == true && t == true) {
-        rw += mwm.blossomValue(i);
-      }
-    }
-    rw -= weight[e] * mwm.dualScale;
-
-    check(rw >= 0, "Negative reduced weight");
-    check(rw == 0 || !mwm.matching(e),
-          "Non-zero reduced weight on matching arc");
-  }
-
-  int pv = 0;
-  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
-    if (mwm.matching(n) != INVALID) {
-      check(mwm.nodeValue(n) >= 0, "Invalid node value");
-      pv += weight[mwm.matching(n)];
-      SmartGraph::Node o = graph.target(mwm.matching(n));
-      check(mwm.mate(n) == o, "Invalid matching");
-      check(mwm.matching(n) == graph.oppositeArc(mwm.matching(o)),
-            "Invalid matching");
-    } else {
-      check(mwm.mate(n) == INVALID, "Invalid matching");
-      check(mwm.nodeValue(n) == 0, "Invalid matching");
-    }
-  }
-
-  int dv = 0;
-  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
-    dv += mwm.nodeValue(n);
-  }
-
-  for (int i = 0; i < mwm.blossomNum(); ++i) {
-    check(mwm.blossomValue(i) >= 0, "Invalid blossom value");
-    check(mwm.blossomSize(i) % 2 == 1, "Even blossom size");
-    dv += mwm.blossomValue(i) * ((mwm.blossomSize(i) - 1) / 2);
-  }
-
-  check(pv * mwm.dualScale == dv * 2, "Wrong duality");
-
-  return;
-}
-
-void checkPerfectMatching(const SmartGraph& graph,
-                          const SmartGraph::EdgeMap<int>& weight,
-                          const MaxWeightedPerfectMatching<SmartGraph>& mwpm) {
-  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
-    if (graph.u(e) == graph.v(e)) continue;
-    int rw =
-      mwpm.nodeValue(graph.u(e)) + mwpm.nodeValue(graph.v(e));
-
-    for (int i = 0; i < mwpm.blossomNum(); ++i) {
-      bool s = false, t = false;
-      for (MaxWeightedPerfectMatching<SmartGraph>::BlossomIt n(mwpm, i);
-           n != INVALID; ++n) {
-        if (graph.u(e) == n) s = true;
-        if (graph.v(e) == n) t = true;
-      }
-      if (s == true && t == true) {
-        rw += mwpm.blossomValue(i);
-      }
-    }
-    rw -= weight[e] * mwpm.dualScale;
-
-    check(rw >= 0, "Negative reduced weight");
-    check(rw == 0 || !mwpm.matching(e),
-          "Non-zero reduced weight on matching arc");
-  }
-
-  int pv = 0;
-  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
-    check(mwpm.matching(n) != INVALID, "Non perfect");
-    pv += weight[mwpm.matching(n)];
-    SmartGraph::Node o = graph.target(mwpm.matching(n));
-    check(mwpm.mate(n) == o, "Invalid matching");
-    check(mwpm.matching(n) == graph.oppositeArc(mwpm.matching(o)),
-          "Invalid matching");
-  }
-
-  int dv = 0;
-  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
-    dv += mwpm.nodeValue(n);
-  }
-
-  for (int i = 0; i < mwpm.blossomNum(); ++i) {
-    check(mwpm.blossomValue(i) >= 0, "Invalid blossom value");
-    check(mwpm.blossomSize(i) % 2 == 1, "Even blossom size");
-    dv += mwpm.blossomValue(i) * ((mwpm.blossomSize(i) - 1) / 2);
-  }
-
-  check(pv * mwpm.dualScale == dv * 2, "Wrong duality");
-
-  return;
-}
-
-
-int main() {
-
-  for (int i = 0; i < lgfn; ++i) {
-    SmartGraph graph;
-    SmartGraph::EdgeMap<int> weight(graph);
-
-    istringstream lgfs(lgf[i]);
-    graphReader(graph, lgfs).
-      edgeMap("weight", weight).run();
-
-    MaxWeightedMatching<SmartGraph> mwm(graph, weight);
-    mwm.run();
-    checkMatching(graph, weight, mwm);
-
-    MaxMatching<SmartGraph> mm(graph);
-    mm.run();
-
-    MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
-    bool perfect = mwpm.run();
-
-    check(perfect == (mm.size() * 2 == countNodes(graph)),
-          "Perfect matching found");
-
-    if (perfect) {
-      checkPerfectMatching(graph, weight, mwpm);
-    }
-  }
-
-  return 0;
-}