Changeset 875:07ec2b52e53d in lemon-main

Timestamp:

03/17/10 10:29:57 (15 years ago)

Author:

Alpar Juttner <alpar@…>

Branch:

default

Parents:

874:d8ea85825e02 (diff), 867:5b926cc36a4b (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Phase:

public

Message:

Merge #356

Location:

lemon

Files:

• matching.h (modified) (64 diffs)
• matching.h (modified) (12 diffs)
• unionfind.h (modified) (1 diff)
• unionfind.h (modified) (1 diff)

lemon/matching.h

@@ -17 +17 @@
  */
 
-#ifndef LEMON_MAX_MATCHING_H
-#define LEMON_MAX_MATCHING_H
+#ifndef LEMON_MATCHING_H
+#define LEMON_MATCHING_H
 
 #include <vector>
@@ -29 +29 @@
 #include <lemon/bin_heap.h>
 #include <lemon/maps.h>
+#include <lemon/fractional_matching.h>
 
 ///\ingroup matching
@@ -42 +43 @@
   /// This class implements Edmonds' alternating forest matching algorithm
   /// for finding a maximum cardinality matching in a general undirected graph.
-  /// It can be started from an arbitrary initial matching 
+  /// It can be started from an arbitrary initial matching
   /// (the default is the empty one).
   ///
@@ -70 +71 @@
     ///\brief Status constants for Gallai-Edmonds decomposition.
     ///
-    ///These constants are used for indicating the Gallai-Edmonds 
+    ///These constants are used for indicating the Gallai-Edmonds
     ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
     ///induce a subgraph with factor-critical components, the nodes with
     ///status \c ODD (or \c A) form the canonical barrier, and the nodes
-    ///with status \c MATCHED (or \c C) induce a subgraph having a 
+    ///with status \c MATCHED (or \c C) induce a subgraph having a
     ///perfect matching.
     enum Status {
@@ -513 +514 @@
     }
 
-    /// \brief Start Edmonds' algorithm with a heuristic improvement 
+    /// \brief Start Edmonds' algorithm with a heuristic improvement
     /// for dense graphs
     ///
@@ -535 +536 @@
     /// \brief Run Edmonds' algorithm
     ///
-    /// This function runs Edmonds' algorithm. An additional heuristic of 
-    /// postponing shrinks is used for relatively dense graphs 
+    /// This function runs Edmonds' algorithm. An additional heuristic of
+    /// postponing shrinks is used for relatively dense graphs
     /// (for which <tt>m>=2*n</tt> holds).
     void run() {
@@ -557 +558 @@
     /// \brief Return the size (cardinality) of the matching.
     ///
-    /// This function returns the size (cardinality) of the current matching. 
+    /// This function returns the size (cardinality) of the current matching.
     /// After run() it returns the size of the maximum matching in the graph.
     int matchingSize() const {
@@ -571 +572 @@
     /// \brief Return \c true if the given edge is in the matching.
     ///
-    /// This function returns \c true if the given edge is in the current 
+    /// This function returns \c true if the given edge is in the current
     /// matching.
     bool matching(const Edge& edge) const {
@@ -580 +581 @@
     ///
     /// This function returns the matching arc (or edge) incident to the
-    /// given node in the current matching or \c INVALID if the node is 
+    /// given node in the current matching or \c INVALID if the node is
     /// not covered by the matching.
     Arc matching(const Node& n) const {
@@ -596 +597 @@
     /// \brief Return the mate of the given node.
     ///
-    /// This function returns the mate of the given node in the current 
+    /// This function returns the mate of the given node in the current
     /// matching or \c INVALID if the node is not covered by the matching.
     Node mate(const Node& n) const {
@@ -606 +607 @@
 
     /// \name Dual Solution
-    /// Functions to get the dual solution, i.e. the Gallai-Edmonds 
+    /// Functions to get the dual solution, i.e. the Gallai-Edmonds
     /// decomposition.
 
@@ -649 +650 @@
   /// \f$O(nm\log n)\f$ time complexity.
   ///
-  /// The maximum weighted matching problem is to find a subset of the 
-  /// edges in an undirected graph with maximum overall weight for which 
+  /// The maximum weighted matching problem is to find a subset of the
+  /// edges in an undirected graph with maximum overall weight for which
   /// each node has at most one incident edge.
   /// It can be formulated with the following linear program.
@@ -674 +675 @@
       \frac{\vert B \vert - 1}{2}z_B\f] */
   ///
-  /// The algorithm can be executed with the run() function. 
+  /// The algorithm can be executed with the run() function.
   /// After it the matching (the primal solution) and the dual solution
-  /// can be obtained using the query functions and the 
-  /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, 
-  /// which is able to iterate on the nodes of a blossom. 
+  /// can be obtained using the query functions and the
+  /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
+  /// which is able to iterate on the nodes of a blossom.
   /// If the value type is integer, then the dual solution is multiplied
   /// by \ref MaxWeightedMatching::dualScale "4".
   ///
   /// \tparam GR The undirected graph type the algorithm runs on.
-  /// \tparam WM The type edge weight map. The default type is 
+  /// \tparam WM The type edge weight map. The default type is
   /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
 #ifdef DOXYGEN
@@ -746 +747 @@
 
     enum Status {
-      EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
+      EVEN = -1, MATCHED = 0, ODD = 1
     };
 
@@ -798 +799 @@
 
     Value _delta_sum;
+    int _unmatched;
+
+    typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
+    FractionalMatching *_fractional;
 
     void createStructures() {
@@ -864 +869 @@
 
     void destroyStructures() {
-      _node_num = countNodes(_graph);
-      _blossom_num = _node_num * 3 / 2;
-
       if (_matching) {
         delete _matching;
@@ -941 +943 @@
             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
               _delta3->push(e, rw / 2);
-            }
-          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
-            if (_delta3->state(e) != _delta3->IN_HEAP) {
-              _delta3->push(e, rw);
             }
           } else {
@@ -969 +967 @@
                                (*_blossom_data)[vb].offset);
                 } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
-                           (*_blossom_data)[vb].offset){
+                           (*_blossom_data)[vb].offset) {
                   _delta2->decrease(vb, _blossom_set->classPrio(vb) -
                                    (*_blossom_data)[vb].offset);
@@ -988 +986 @@
       if (!_blossom_set->trivial(blossom)) {
         _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
-                     (*_blossom_data)[blossom].offset);
+                      (*_blossom_data)[blossom].offset);
       }
     }
@@ -1055 +1053 @@
                 }
               }
-            }
-
-          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
-            if (_delta3->state(e) == _delta3->IN_HEAP) {
-              _delta3->erase(e);
             }
           } else {
@@ -1097 +1090 @@
           std::numeric_limits<Value>::max()) {
         _delta2->push(blossom, _blossom_set->classPrio(blossom) -
-                       (*_blossom_data)[blossom].offset);
+                      (*_blossom_data)[blossom].offset);
       }
 
@@ -1136 +1129 @@
             if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
               _delta3->push(e, rw / 2);
-            }
-          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
-            if (_delta3->state(e) != _delta3->IN_HEAP) {
-              _delta3->push(e, rw);
             }
           } else {
@@ -1177 +1166 @@
     }
 
-
-    void matchedToUnmatched(int blossom) {
-      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
-        _delta2->erase(blossom);
-      }
-
-      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
-           n != INVALID; ++n) {
-        int ni = (*_node_index)[n];
-
-        _blossom_set->increase(n, std::numeric_limits<Value>::max());
-
-        (*_node_data)[ni].heap.clear();
-        (*_node_data)[ni].heap_index.clear();
-
-        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
-          Node v = _graph.target(e);
-          int vb = _blossom_set->find(v);
-          int vi = (*_node_index)[v];
-
-          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
-            dualScale * _weight[e];
-
-          if ((*_blossom_data)[vb].status == EVEN) {
-            if (_delta3->state(e) != _delta3->IN_HEAP) {
-              _delta3->push(e, rw);
-            }
-          }
-        }
-      }
-    }
-
-    void unmatchedToMatched(int blossom) {
-      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
-           n != INVALID; ++n) {
-        int ni = (*_node_index)[n];
-
-        for (InArcIt e(_graph, n); e != INVALID; ++e) {
-          Node v = _graph.source(e);
-          int vb = _blossom_set->find(v);
-          int vi = (*_node_index)[v];
-
-          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
-            dualScale * _weight[e];
-
-          if (vb == blossom) {
-            if (_delta3->state(e) == _delta3->IN_HEAP) {
-              _delta3->erase(e);
-            }
-          } else if ((*_blossom_data)[vb].status == EVEN) {
-
-            if (_delta3->state(e) == _delta3->IN_HEAP) {
-              _delta3->erase(e);
-            }
-
-            int vt = _tree_set->find(vb);
-
-            Arc r = _graph.oppositeArc(e);
-
-            typename std::map<int, Arc>::iterator it =
-              (*_node_data)[ni].heap_index.find(vt);
-
-            if (it != (*_node_data)[ni].heap_index.end()) {
-              if ((*_node_data)[ni].heap[it->second] > rw) {
-                (*_node_data)[ni].heap.replace(it->second, r);
-                (*_node_data)[ni].heap.decrease(r, rw);
-                it->second = r;
-              }
-            } else {
-              (*_node_data)[ni].heap.push(r, rw);
-              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
-            }
-
-            if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
-              _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
-
-              if (_delta2->state(blossom) != _delta2->IN_HEAP) {
-                _delta2->push(blossom, _blossom_set->classPrio(blossom) -
-                             (*_blossom_data)[blossom].offset);
-              } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
-                         (*_blossom_data)[blossom].offset){
-                _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
-                                 (*_blossom_data)[blossom].offset);
-              }
-            }
-
-          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
-            if (_delta3->state(e) == _delta3->IN_HEAP) {
-              _delta3->erase(e);
-            }
-          }
-        }
-      }
-    }
-
     void alternatePath(int even, int tree) {
       int odd;
@@ -1314 +1208 @@
       destroyTree(tree);
 
-      (*_blossom_data)[blossom].status = UNMATCHED;
       (*_blossom_data)[blossom].base = node;
-      matchedToUnmatched(blossom);
-    }
-
+      (*_blossom_data)[blossom].next = INVALID;
+    }
 
     void augmentOnEdge(const Edge& edge) {
@@ -1325 +1217 @@
       int right = _blossom_set->find(_graph.v(edge));
 
-      if ((*_blossom_data)[left].status == EVEN) {
-        int left_tree = _tree_set->find(left);
-        alternatePath(left, left_tree);
-        destroyTree(left_tree);
-      } else {
-        (*_blossom_data)[left].status = MATCHED;
-        unmatchedToMatched(left);
-      }
-
-      if ((*_blossom_data)[right].status == EVEN) {
-        int right_tree = _tree_set->find(right);
-        alternatePath(right, right_tree);
-        destroyTree(right_tree);
-      } else {
-        (*_blossom_data)[right].status = MATCHED;
-        unmatchedToMatched(right);
-      }
+      int left_tree = _tree_set->find(left);
+      alternatePath(left, left_tree);
+      destroyTree(left_tree);
+
+      int right_tree = _tree_set->find(right);
+      alternatePath(right, right_tree);
+      destroyTree(right_tree);
 
       (*_blossom_data)[left].next = _graph.direct(edge, true);
       (*_blossom_data)[right].next = _graph.direct(edge, false);
+    }
+
+    void augmentOnArc(const Arc& arc) {
+
+      int left = _blossom_set->find(_graph.source(arc));
+      int right = _blossom_set->find(_graph.target(arc));
+
+      (*_blossom_data)[left].status = MATCHED;
+
+      int right_tree = _tree_set->find(right);
+      alternatePath(right, right_tree);
+      destroyTree(right_tree);
+
+      (*_blossom_data)[left].next = arc;
+      (*_blossom_data)[right].next = _graph.oppositeArc(arc);
     }
 
@@ -1549 +1446 @@
           (*_blossom_data)[sb].pred = pred;
           (*_blossom_data)[sb].next =
-                           _graph.oppositeArc((*_blossom_data)[tb].next);
+            _graph.oppositeArc((*_blossom_data)[tb].next);
 
           pred = (*_blossom_data)[ub].next;
@@ -1649 +1546 @@
 
       for (int i = 0; i < int(blossoms.size()); ++i) {
-        if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
+        if ((*_blossom_data)[blossoms[i]].next != INVALID) {
 
           Value offset = (*_blossom_data)[blossoms[i]].offset;
@@ -1687 +1584 @@
         _delta4_index(0), _delta4(0),
 
-        _delta_sum() {}
+        _delta_sum(), _unmatched(0),
+
+        _fractional(0)
+    {}
 
     ~MaxWeightedMatching() {
       destroyStructures();
+      if (_fractional) {
+        delete _fractional;
+      }
     }
 
@@ -1722 +1625 @@
       }
       
+      _unmatched = _node_num;
+
       _delta1->clear();
       _delta2->clear();
@@ -1765 +1670 @@
     }
 
+    /// \brief Initialize the algorithm with fractional matching
+    ///
+    /// This function initializes the algorithm with a fractional
+    /// matching. This initialization is also called jumpstart heuristic.
+    void fractionalInit() {
+      createStructures();
+      
+      if (_fractional == 0) {
+        _fractional = new FractionalMatching(_graph, _weight, false);
+      }
+      _fractional->run();
+
+      for (ArcIt e(_graph); e != INVALID; ++e) {
+        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
+      }
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        (*_delta1_index)[n] = _delta1->PRE_HEAP;
+      }
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+      }
+      for (int i = 0; i < _blossom_num; ++i) {
+        (*_delta2_index)[i] = _delta2->PRE_HEAP;
+        (*_delta4_index)[i] = _delta4->PRE_HEAP;
+      }
+
+      _unmatched = 0;
+
+      int index = 0;
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        Value pot = _fractional->nodeValue(n);
+        (*_node_index)[n] = index;
+        (*_node_data)[index].pot = pot;
+        int blossom =
+          _blossom_set->insert(n, std::numeric_limits<Value>::max());
+
+        (*_blossom_data)[blossom].status = MATCHED;
+        (*_blossom_data)[blossom].pred = INVALID;
+        (*_blossom_data)[blossom].next = _fractional->matching(n);
+        if (_fractional->matching(n) == INVALID) {
+          (*_blossom_data)[blossom].base = n;
+        }
+        (*_blossom_data)[blossom].pot = 0;
+        (*_blossom_data)[blossom].offset = 0;
+        ++index;
+      }
+
+      typename Graph::template NodeMap<bool> processed(_graph, false);
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        if (processed[n]) continue;
+        processed[n] = true;
+        if (_fractional->matching(n) == INVALID) continue;
+        int num = 1;
+        Node v = _graph.target(_fractional->matching(n));
+        while (n != v) {
+          processed[v] = true;
+          v = _graph.target(_fractional->matching(v));
+          ++num;
+        }
+
+        if (num % 2 == 1) {
+          std::vector<int> subblossoms(num);
+
+          subblossoms[--num] = _blossom_set->find(n);
+          _delta1->push(n, _fractional->nodeValue(n));
+          v = _graph.target(_fractional->matching(n));
+          while (n != v) {
+            subblossoms[--num] = _blossom_set->find(v);
+            _delta1->push(v, _fractional->nodeValue(v));
+            v = _graph.target(_fractional->matching(v));            
+          }
+          
+          int surface = 
+            _blossom_set->join(subblossoms.begin(), subblossoms.end());
+          (*_blossom_data)[surface].status = EVEN;
+          (*_blossom_data)[surface].pred = INVALID;
+          (*_blossom_data)[surface].next = INVALID;
+          (*_blossom_data)[surface].pot = 0;
+          (*_blossom_data)[surface].offset = 0;
+          
+          _tree_set->insert(surface);
+          ++_unmatched;
+        }
+      }
+
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+        int si = (*_node_index)[_graph.u(e)];
+        int sb = _blossom_set->find(_graph.u(e));
+        int ti = (*_node_index)[_graph.v(e)];
+        int tb = _blossom_set->find(_graph.v(e));
+        if ((*_blossom_data)[sb].status == EVEN &&
+            (*_blossom_data)[tb].status == EVEN && sb != tb) {
+          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
+                            dualScale * _weight[e]) / 2);
+        }
+      }
+
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        int nb = _blossom_set->find(n);
+        if ((*_blossom_data)[nb].status != MATCHED) continue;
+        int ni = (*_node_index)[n];
+
+        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
+          Node v = _graph.target(e);
+          int vb = _blossom_set->find(v);
+          int vi = (*_node_index)[v];
+
+          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
+            dualScale * _weight[e];
+
+          if ((*_blossom_data)[vb].status == EVEN) {
+
+            int vt = _tree_set->find(vb);
+
+            typename std::map<int, Arc>::iterator it =
+              (*_node_data)[ni].heap_index.find(vt);
+
+            if (it != (*_node_data)[ni].heap_index.end()) {
+              if ((*_node_data)[ni].heap[it->second] > rw) {
+                (*_node_data)[ni].heap.replace(it->second, e);
+                (*_node_data)[ni].heap.decrease(e, rw);
+                it->second = e;
+              }
+            } else {
+              (*_node_data)[ni].heap.push(e, rw);
+              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
+            }
+          }
+        }
+            
+        if (!(*_node_data)[ni].heap.empty()) {
+          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
+          _delta2->push(nb, _blossom_set->classPrio(nb));
+        }
+      }
+    }
+
     /// \brief Start the algorithm
     ///
     /// This function starts the algorithm.
     ///
-    /// \pre \ref init() must be called before using this function.
+    /// \pre \ref init() or \ref fractionalInit() must be called
+    /// before using this function.
     void start() {
       enum OpType {
@@ -1775 +1818 @@
       };
 
-      int unmatched = _node_num;
-      while (unmatched > 0) {
+      while (_unmatched > 0) {
         Value d1 = !_delta1->empty() ?
           _delta1->prio() : std::numeric_limits<Value>::max();
@@ -1789 +1831 @@
           _delta4->prio() : std::numeric_limits<Value>::max();
 
-        _delta_sum = d1; OpType ot = D1;
+        _delta_sum = d3; OpType ot = D3;
+        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
         if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
-        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
         if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
-
 
         switch (ot) {
@@ -1800 +1841 @@
             Node n = _delta1->top();
             unmatchNode(n);
-            --unmatched;
+            --_unmatched;
           }
           break;
@@ -1807 +1848 @@
             int blossom = _delta2->top();
             Node n = _blossom_set->classTop(blossom);
-            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
-            extendOnArc(e);
+            Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
+            if ((*_blossom_data)[blossom].next == INVALID) {
+              augmentOnArc(a);
+              --_unmatched;
+            } else {
+              extendOnArc(a);
+            }
           }
           break;
@@ -1821 +1867 @@
               _delta3->pop();
             } else {
-              int left_tree;
-              if ((*_blossom_data)[left_blossom].status == EVEN) {
-                left_tree = _tree_set->find(left_blossom);
-              } else {
-                left_tree = -1;
-                ++unmatched;
-              }
-              int right_tree;
-              if ((*_blossom_data)[right_blossom].status == EVEN) {
-                right_tree = _tree_set->find(right_blossom);
-              } else {
-                right_tree = -1;
-                ++unmatched;
-              }
+              int left_tree = _tree_set->find(left_blossom);
+              int right_tree = _tree_set->find(right_blossom);
 
               if (left_tree == right_tree) {
@@ -1840 +1874 @@
               } else {
                 augmentOnEdge(e);
-                unmatched -= 2;
+                _unmatched -= 2;
               }
             }
@@ -1858 +1892 @@
     /// \note mwm.run() is just a shortcut of the following code.
     /// \code
-    ///   mwm.init();
+    ///   mwm.fractionalInit();
     ///   mwm.start();
     /// \endcode
     void run() {
-      init();
+      fractionalInit();
       start();
     }
@@ -1869 +1903 @@
 
     /// \name Primal Solution
-    /// Functions to get the primal solution, i.e. the maximum weighted 
+    /// Functions to get the primal solution, i.e. the maximum weighted
     /// matching.\n
     /// Either \ref run() or \ref start() function should be called before
@@ -1888 +1922 @@
         }
       }
-      return sum /= 2;
+      return sum / 2;
     }
 
@@ -1908 +1942 @@
     /// \brief Return \c true if the given edge is in the matching.
     ///
-    /// This function returns \c true if the given edge is in the found 
+    /// This function returns \c true if the given edge is in the found
     /// matching.
     ///
@@ -1919 +1953 @@
     ///
     /// This function returns the matching arc (or edge) incident to the
-    /// given node in the found matching or \c INVALID if the node is 
+    /// given node in the found matching or \c INVALID if the node is
     /// not covered by the matching.
     ///
@@ -1937 +1971 @@
     /// \brief Return the mate of the given node.
     ///
-    /// This function returns the mate of the given node in the found 
+    /// This function returns the mate of the given node in the found
     /// matching or \c INVALID if the node is not covered by the matching.
     ///
@@ -1957 +1991 @@
     /// \brief Return the value of the dual solution.
     ///
-    /// This function returns the value of the dual solution. 
-    /// It should be equal to the primal value scaled by \ref dualScale 
+    /// This function returns the value of the dual solution.
+    /// It should be equal to the primal value scaled by \ref dualScale
     /// "dual scale".
     ///
@@ -2013 +2047 @@
     /// \brief Iterator for obtaining the nodes of a blossom.
     ///
-    /// This class provides an iterator for obtaining the nodes of the 
+    /// This class provides an iterator for obtaining the nodes of the
     /// given blossom. It lists a subset of the nodes.
-    /// Before using this iterator, you must allocate a 
+    /// Before using this iterator, you must allocate a
     /// MaxWeightedMatching class and execute it.
     class BlossomIt {
@@ -2024 +2058 @@
       /// Constructor to get the nodes of the given variable.
       ///
-      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or 
-      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be 
+      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
+      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
       /// called before initializing this iterator.
       BlossomIt(const MaxWeightedMatching& algorithm, int variable)
@@ -2078 +2112 @@
   /// \f$O(nm\log n)\f$ time complexity.
   ///
-  /// The maximum weighted perfect matching problem is to find a subset of 
-  /// the edges in an undirected graph with maximum overall weight for which 
+  /// The maximum weighted perfect matching problem is to find a subset of
+  /// the edges in an undirected graph with maximum overall weight for which
   /// each node has exactly one incident edge.
   /// It can be formulated with the following linear program.
@@ -2102 +2136 @@
       \frac{\vert B \vert - 1}{2}z_B\f] */
   ///
-  /// The algorithm can be executed with the run() function. 
+  /// The algorithm can be executed with the run() function.
   /// After it the matching (the primal solution) and the dual solution
-  /// can be obtained using the query functions and the 
-  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, 
-  /// which is able to iterate on the nodes of a blossom. 
+  /// can be obtained using the query functions and the
+  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
+  /// which is able to iterate on the nodes of a blossom.
   /// If the value type is integer, then the dual solution is multiplied
   /// by \ref MaxWeightedMatching::dualScale "4".
   ///
   /// \tparam GR The undirected graph type the algorithm runs on.
-  /// \tparam WM The type edge weight map. The default type is 
+  /// \tparam WM The type edge weight map. The default type is
   /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
 #ifdef DOXYGEN
@@ -2222 +2256 @@
 
     Value _delta_sum;
+    int _unmatched;
+
+    typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap> 
+    FractionalMatching;
+    FractionalMatching *_fractional;
 
     void createStructures() {
@@ -2283 +2322 @@
 
     void destroyStructures() {
-      _node_num = countNodes(_graph);
-      _blossom_num = _node_num * 3 / 2;
-
       if (_matching) {
         delete _matching;
@@ -2958 +2994 @@
         _delta4_index(0), _delta4(0),
 
-        _delta_sum() {}
+        _delta_sum(), _unmatched(0),
+
+        _fractional(0)
+    {}
 
     ~MaxWeightedPerfectMatching() {
       destroyStructures();
+      if (_fractional) {
+        delete _fractional;
+      }
     }
 
@@ -2989 +3031 @@
         (*_delta4_index)[i] = _delta4->PRE_HEAP;
       }
+
+      _unmatched = _node_num;
 
       _delta2->clear();
@@ -3031 +3075 @@
     }
 
+    /// \brief Initialize the algorithm with fractional matching
+    ///
+    /// This function initializes the algorithm with a fractional
+    /// matching. This initialization is also called jumpstart heuristic.
+    void fractionalInit() {
+      createStructures();
+      
+      if (_fractional == 0) {
+        _fractional = new FractionalMatching(_graph, _weight, false);
+      }
+      if (!_fractional->run()) {
+        _unmatched = -1;
+        return;
+      }
+
+      for (ArcIt e(_graph); e != INVALID; ++e) {
+        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
+      }
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+        (*_delta3_index)[e] = _delta3->PRE_HEAP;
+      }
+      for (int i = 0; i < _blossom_num; ++i) {
+        (*_delta2_index)[i] = _delta2->PRE_HEAP;
+        (*_delta4_index)[i] = _delta4->PRE_HEAP;
+      }
+
+      _unmatched = 0;
+
+      int index = 0;
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        Value pot = _fractional->nodeValue(n);
+        (*_node_index)[n] = index;
+        (*_node_data)[index].pot = pot;
+        int blossom =
+          _blossom_set->insert(n, std::numeric_limits<Value>::max());
+
+        (*_blossom_data)[blossom].status = MATCHED;
+        (*_blossom_data)[blossom].pred = INVALID;
+        (*_blossom_data)[blossom].next = _fractional->matching(n);
+        (*_blossom_data)[blossom].pot = 0;
+        (*_blossom_data)[blossom].offset = 0;
+        ++index;
+      }
+
+      typename Graph::template NodeMap<bool> processed(_graph, false);
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        if (processed[n]) continue;
+        processed[n] = true;
+        if (_fractional->matching(n) == INVALID) continue;
+        int num = 1;
+        Node v = _graph.target(_fractional->matching(n));
+        while (n != v) {
+          processed[v] = true;
+          v = _graph.target(_fractional->matching(v));
+          ++num;
+        }
+
+        if (num % 2 == 1) {
+          std::vector<int> subblossoms(num);
+
+          subblossoms[--num] = _blossom_set->find(n);
+          v = _graph.target(_fractional->matching(n));
+          while (n != v) {
+            subblossoms[--num] = _blossom_set->find(v);
+            v = _graph.target(_fractional->matching(v));            
+          }
+          
+          int surface = 
+            _blossom_set->join(subblossoms.begin(), subblossoms.end());
+          (*_blossom_data)[surface].status = EVEN;
+          (*_blossom_data)[surface].pred = INVALID;
+          (*_blossom_data)[surface].next = INVALID;
+          (*_blossom_data)[surface].pot = 0;
+          (*_blossom_data)[surface].offset = 0;
+          
+          _tree_set->insert(surface);
+          ++_unmatched;
+        }
+      }
+
+      for (EdgeIt e(_graph); e != INVALID; ++e) {
+        int si = (*_node_index)[_graph.u(e)];
+        int sb = _blossom_set->find(_graph.u(e));
+        int ti = (*_node_index)[_graph.v(e)];
+        int tb = _blossom_set->find(_graph.v(e));
+        if ((*_blossom_data)[sb].status == EVEN &&
+            (*_blossom_data)[tb].status == EVEN && sb != tb) {
+          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
+                            dualScale * _weight[e]) / 2);
+        }
+      }
+
+      for (NodeIt n(_graph); n != INVALID; ++n) {
+        int nb = _blossom_set->find(n);
+        if ((*_blossom_data)[nb].status != MATCHED) continue;
+        int ni = (*_node_index)[n];
+
+        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
+          Node v = _graph.target(e);
+          int vb = _blossom_set->find(v);
+          int vi = (*_node_index)[v];
+
+          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
+            dualScale * _weight[e];
+
+          if ((*_blossom_data)[vb].status == EVEN) {
+
+            int vt = _tree_set->find(vb);
+
+            typename std::map<int, Arc>::iterator it =
+              (*_node_data)[ni].heap_index.find(vt);
+
+            if (it != (*_node_data)[ni].heap_index.end()) {
+              if ((*_node_data)[ni].heap[it->second] > rw) {
+                (*_node_data)[ni].heap.replace(it->second, e);
+                (*_node_data)[ni].heap.decrease(e, rw);
+                it->second = e;
+              }
+            } else {
+              (*_node_data)[ni].heap.push(e, rw);
+              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
+            }
+          }
+        }
+            
+        if (!(*_node_data)[ni].heap.empty()) {
+          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
+          _delta2->push(nb, _blossom_set->classPrio(nb));
+        }
+      }
+    }
+
     /// \brief Start the algorithm
     ///
     /// This function starts the algorithm.
     ///
-    /// \pre \ref init() must be called before using this function.
+    /// \pre \ref init() or \ref fractionalInit() must be called before
+    /// using this function.
     bool start() {
       enum OpType {
@@ -3041 +3218 @@
       };
 
-      int unmatched = _node_num;
-      while (unmatched > 0) {
+      if (_unmatched == -1) return false;
+
+      while (_unmatched > 0) {
         Value d2 = !_delta2->empty() ?
           _delta2->prio() : std::numeric_limits<Value>::max();
@@ -3052 +3230 @@
           _delta4->prio() : std::numeric_limits<Value>::max();
 
-        _delta_sum = d2; OpType ot = D2;
-        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
+        _delta_sum = d3; OpType ot = D3;
+        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
         if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
 
@@ -3086 +3264 @@
               } else {
                 augmentOnEdge(e);
-                unmatched -= 2;
+                _unmatched -= 2;
               }
             }
@@ -3105 +3283 @@
     /// \note mwpm.run() is just a shortcut of the following code.
     /// \code
-    ///   mwpm.init();
+    ///   mwpm.fractionalInit();
     ///   mwpm.start();
     /// \endcode
     bool run() {
-      init();
+      fractionalInit();
       return start();
     }
@@ -3116 +3294 @@
 
     /// \name Primal Solution
-    /// Functions to get the primal solution, i.e. the maximum weighted 
+    /// Functions to get the primal solution, i.e. the maximum weighted
     /// perfect matching.\n
     /// Either \ref run() or \ref start() function should be called before
@@ -3135 +3313 @@
         }
       }
-      return sum /= 2;
+      return sum / 2;
     }
 
     /// \brief Return \c true if the given edge is in the matching.
     ///
-    /// This function returns \c true if the given edge is in the found 
+    /// This function returns \c true if the given edge is in the found
     /// matching.
     ///
@@ -3151 +3329 @@
     ///
     /// This function returns the matching arc (or edge) incident to the
-    /// given node in the found matching or \c INVALID if the node is 
+    /// given node in the found matching or \c INVALID if the node is
     /// not covered by the matching.
     ///
@@ -3169 +3347 @@
     /// \brief Return the mate of the given node.
     ///
-    /// This function returns the mate of the given node in the found 
+    /// This function returns the mate of the given node in the found
     /// matching or \c INVALID if the node is not covered by the matching.
     ///
@@ -3188 +3366 @@
     /// \brief Return the value of the dual solution.
     ///
-    /// This function returns the value of the dual solution. 
-    /// It should be equal to the primal value scaled by \ref dualScale 
+    /// This function returns the value of the dual solution.
+    /// It should be equal to the primal value scaled by \ref dualScale
     /// "dual scale".
     ///
@@ -3244 +3422 @@
     /// \brief Iterator for obtaining the nodes of a blossom.
     ///
-    /// This class provides an iterator for obtaining the nodes of the 
+    /// This class provides an iterator for obtaining the nodes of the
     /// given blossom. It lists a subset of the nodes.
-    /// Before using this iterator, you must allocate a 
+    /// Before using this iterator, you must allocate a
     /// MaxWeightedPerfectMatching class and execute it.
     class BlossomIt {
@@ -3255 +3433 @@
       /// Constructor to get the nodes of the given variable.
       ///
-      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" 
-      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" 
+      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
+      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
       /// must be called before initializing this iterator.
       BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
@@ -3302 +3480 @@
 } //END OF NAMESPACE LEMON
 
-#endif //LEMON_MAX_MATCHING_H
+#endif //LEMON_MATCHING_H

lemon/matching.h

@@ -811 +811 @@
         _matching = new MatchingMap(_graph);
       }
+
       if (!_node_potential) {
         _node_potential = new NodePotential(_graph);
       }
+
       if (!_blossom_set) {
         _blossom_index = new IntNodeMap(_graph);
         _blossom_set = new BlossomSet(*_blossom_index);
         _blossom_data = new RangeMap<BlossomData>(_blossom_num);
+      } else if (_blossom_data->size() != _blossom_num) {
+        delete _blossom_data;
+        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
       }
 
@@ -824 +829 @@
         _node_heap_index = new IntArcMap(_graph);
         _node_data = new RangeMap<NodeData>(_node_num,
-                                              NodeData(*_node_heap_index));
+                                            NodeData(*_node_heap_index));
+      } else {
+        delete _node_data;
+        _node_data = new RangeMap<NodeData>(_node_num,
+                                            NodeData(*_node_heap_index));
       }
 
@@ -830 +839 @@
         _tree_set_index = new IntIntMap(_blossom_num);
         _tree_set = new TreeSet(*_tree_set_index);
-      }
+      } else {
+        _tree_set_index->resize(_blossom_num);
+      }
+
       if (!_delta1) {
         _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);
-      }
+      } else {
+        _delta2_index->resize(_blossom_num);
+      }
+
       if (!_delta3) {
         _delta3_index = new IntEdgeMap(_graph);
         _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
       }
+
       if (!_delta4) {
         _delta4_index = new IntIntMap(_blossom_num);
         _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
+      } else {
+        _delta4_index->resize(_blossom_num);
       }
     }
@@ -1589 +1608 @@
       createStructures();
 
+      _blossom_node_list.clear();
+      _blossom_potential.clear();
+
       for (ArcIt e(_graph); e != INVALID; ++e) {
         (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
@@ -1602 +1624 @@
         (*_delta4_index)[i] = _delta4->PRE_HEAP;
       }
-
+      
       _unmatched = _node_num;
+
+      _delta1->clear();
+      _delta2->clear();
+      _delta3->clear();
+      _delta4->clear();
+      _blossom_set->clear();
+      _tree_set->clear();
 
       int index = 0;
@@ -1615 +1644 @@
         }
         (*_node_index)[n] = index;
+        (*_node_data)[index].heap_index.clear();
+        (*_node_data)[index].heap.clear();
         (*_node_data)[index].pot = max;
         _delta1->push(n, max);
@@ -2238 +2269 @@
         _matching = new MatchingMap(_graph);
       }
+
       if (!_node_potential) {
         _node_potential = new NodePotential(_graph);
       }
+
       if (!_blossom_set) {
         _blossom_index = new IntNodeMap(_graph);
         _blossom_set = new BlossomSet(*_blossom_index);
+        _blossom_data = new RangeMap<BlossomData>(_blossom_num);
+      } else if (_blossom_data->size() != _blossom_num) {
+        delete _blossom_data;
         _blossom_data = new RangeMap<BlossomData>(_blossom_num);
       }
@@ -2252 +2288 @@
         _node_data = new RangeMap<NodeData>(_node_num,
                                             NodeData(*_node_heap_index));
+      } else if (_node_data->size() != _node_num) {
+        delete _node_data;
+        _node_data = new RangeMap<NodeData>(_node_num,
+                                            NodeData(*_node_heap_index));
       }
 
@@ -2257 +2297 @@
         _tree_set_index = new IntIntMap(_blossom_num);
         _tree_set = new TreeSet(*_tree_set_index);
-      }
+      } else {
+        _tree_set_index->resize(_blossom_num);
+      }
+
       if (!_delta2) {
         _delta2_index = new IntIntMap(_blossom_num);
         _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
-      }
+      } else {
+        _delta2_index->resize(_blossom_num);
+      }
+
       if (!_delta3) {
         _delta3_index = new IntEdgeMap(_graph);
         _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
       }
+
       if (!_delta4) {
         _delta4_index = new IntIntMap(_blossom_num);
         _delta4 = new BinHeap<Value, IntIntMap>(*_delta4_index);
+      } else {
+        _delta4_index->resize(_blossom_num);
       }
     }
@@ -2969 +3018 @@
       createStructures();
 
+      _blossom_node_list.clear();
+      _blossom_potential.clear();
+
       for (ArcIt e(_graph); e != INVALID; ++e) {
         (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
@@ -2981 +3033 @@
 
       _unmatched = _node_num;
+
+      _delta2->clear();
+      _delta3->clear();
+      _delta4->clear();
+      _blossom_set->clear();
+      _tree_set->clear();
 
       int index = 0;
@@ -2992 +3050 @@
         }
         (*_node_index)[n] = index;
+        (*_node_data)[index].heap_index.clear();
+        (*_node_data)[index].heap.clear();
         (*_node_data)[index].pot = max;
         int blossom =

lemon/unionfind.h

@@ -44 +44 @@
   /// This is a very simple but efficient implementation, providing
   /// only four methods: join (union), find, insert and size.
-  /// For more features see the \ref UnionFindEnum class.
+  /// For more features, see the \ref UnionFindEnum class.
   ///
   /// It is primarily used in Kruskal algorithm for finding minimal

lemon/unionfind.h

@@ -1289 +1289 @@
         first_free_class(-1), first_free_node(-1) {}
 
+    /// \brief Clears the union-find data structure
+    ///
+    /// Erase each item from the data structure.
+    void clear() {
+      nodes.clear();
+      classes.clear();
+      first_free_node = first_free_class = first_class = -1;
+    }
+
     /// \brief Insert a new node into a new component.
     ///