General improvements in weighted matching algorithms (#314)
authorBalazs Dezso <deba@inf.elte.hu>
Sun, 20 Sep 2009 21:38:24 +0200
changeset 8680513ccfea967
parent 712 6d5f547e5bfb
child 869 636dadefe1e6
General improvements in weighted matching algorithms (#314)

- Fix include guard
- Uniform handling of MATCHED and UNMATCHED blossoms
- Prefer operations which decrease the number of trees
- Fix improper use of '/='

The solved problems did not cause wrong solution.
lemon/matching.h
     1.1 --- a/lemon/matching.h	Mon Aug 31 20:27:38 2009 +0200
     1.2 +++ b/lemon/matching.h	Sun Sep 20 21:38:24 2009 +0200
     1.3 @@ -16,8 +16,8 @@
     1.4   *
     1.5   */
     1.6  
     1.7 -#ifndef LEMON_MAX_MATCHING_H
     1.8 -#define LEMON_MAX_MATCHING_H
     1.9 +#ifndef LEMON_MATCHING_H
    1.10 +#define LEMON_MATCHING_H
    1.11  
    1.12  #include <vector>
    1.13  #include <queue>
    1.14 @@ -41,7 +41,7 @@
    1.15    ///
    1.16    /// This class implements Edmonds' alternating forest matching algorithm
    1.17    /// for finding a maximum cardinality matching in a general undirected graph.
    1.18 -  /// It can be started from an arbitrary initial matching 
    1.19 +  /// It can be started from an arbitrary initial matching
    1.20    /// (the default is the empty one).
    1.21    ///
    1.22    /// The dual solution of the problem is a map of the nodes to
    1.23 @@ -69,11 +69,11 @@
    1.24  
    1.25      ///\brief Status constants for Gallai-Edmonds decomposition.
    1.26      ///
    1.27 -    ///These constants are used for indicating the Gallai-Edmonds 
    1.28 +    ///These constants are used for indicating the Gallai-Edmonds
    1.29      ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
    1.30      ///induce a subgraph with factor-critical components, the nodes with
    1.31      ///status \c ODD (or \c A) form the canonical barrier, and the nodes
    1.32 -    ///with status \c MATCHED (or \c C) induce a subgraph having a 
    1.33 +    ///with status \c MATCHED (or \c C) induce a subgraph having a
    1.34      ///perfect matching.
    1.35      enum Status {
    1.36        EVEN = 1,       ///< = 1. (\c D is an alias for \c EVEN.)
    1.37 @@ -512,7 +512,7 @@
    1.38        }
    1.39      }
    1.40  
    1.41 -    /// \brief Start Edmonds' algorithm with a heuristic improvement 
    1.42 +    /// \brief Start Edmonds' algorithm with a heuristic improvement
    1.43      /// for dense graphs
    1.44      ///
    1.45      /// This function runs Edmonds' algorithm with a heuristic of postponing
    1.46 @@ -534,8 +534,8 @@
    1.47  
    1.48      /// \brief Run Edmonds' algorithm
    1.49      ///
    1.50 -    /// This function runs Edmonds' algorithm. An additional heuristic of 
    1.51 -    /// postponing shrinks is used for relatively dense graphs 
    1.52 +    /// This function runs Edmonds' algorithm. An additional heuristic of
    1.53 +    /// postponing shrinks is used for relatively dense graphs
    1.54      /// (for which <tt>m>=2*n</tt> holds).
    1.55      void run() {
    1.56        if (countEdges(_graph) < 2 * countNodes(_graph)) {
    1.57 @@ -556,7 +556,7 @@
    1.58  
    1.59      /// \brief Return the size (cardinality) of the matching.
    1.60      ///
    1.61 -    /// This function returns the size (cardinality) of the current matching. 
    1.62 +    /// This function returns the size (cardinality) of the current matching.
    1.63      /// After run() it returns the size of the maximum matching in the graph.
    1.64      int matchingSize() const {
    1.65        int size = 0;
    1.66 @@ -570,7 +570,7 @@
    1.67  
    1.68      /// \brief Return \c true if the given edge is in the matching.
    1.69      ///
    1.70 -    /// This function returns \c true if the given edge is in the current 
    1.71 +    /// This function returns \c true if the given edge is in the current
    1.72      /// matching.
    1.73      bool matching(const Edge& edge) const {
    1.74        return edge == (*_matching)[_graph.u(edge)];
    1.75 @@ -579,7 +579,7 @@
    1.76      /// \brief Return the matching arc (or edge) incident to the given node.
    1.77      ///
    1.78      /// This function returns the matching arc (or edge) incident to the
    1.79 -    /// given node in the current matching or \c INVALID if the node is 
    1.80 +    /// given node in the current matching or \c INVALID if the node is
    1.81      /// not covered by the matching.
    1.82      Arc matching(const Node& n) const {
    1.83        return (*_matching)[n];
    1.84 @@ -595,7 +595,7 @@
    1.85  
    1.86      /// \brief Return the mate of the given node.
    1.87      ///
    1.88 -    /// This function returns the mate of the given node in the current 
    1.89 +    /// This function returns the mate of the given node in the current
    1.90      /// matching or \c INVALID if the node is not covered by the matching.
    1.91      Node mate(const Node& n) const {
    1.92        return (*_matching)[n] != INVALID ?
    1.93 @@ -605,7 +605,7 @@
    1.94      /// @}
    1.95  
    1.96      /// \name Dual Solution
    1.97 -    /// Functions to get the dual solution, i.e. the Gallai-Edmonds 
    1.98 +    /// Functions to get the dual solution, i.e. the Gallai-Edmonds
    1.99      /// decomposition.
   1.100  
   1.101      /// @{
   1.102 @@ -648,8 +648,8 @@
   1.103    /// on extensive use of priority queues and provides
   1.104    /// \f$O(nm\log n)\f$ time complexity.
   1.105    ///
   1.106 -  /// The maximum weighted matching problem is to find a subset of the 
   1.107 -  /// edges in an undirected graph with maximum overall weight for which 
   1.108 +  /// The maximum weighted matching problem is to find a subset of the
   1.109 +  /// edges in an undirected graph with maximum overall weight for which
   1.110    /// each node has at most one incident edge.
   1.111    /// It can be formulated with the following linear program.
   1.112    /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
   1.113 @@ -673,16 +673,16 @@
   1.114    /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
   1.115        \frac{\vert B \vert - 1}{2}z_B\f] */
   1.116    ///
   1.117 -  /// The algorithm can be executed with the run() function. 
   1.118 +  /// The algorithm can be executed with the run() function.
   1.119    /// After it the matching (the primal solution) and the dual solution
   1.120 -  /// can be obtained using the query functions and the 
   1.121 -  /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, 
   1.122 -  /// which is able to iterate on the nodes of a blossom. 
   1.123 +  /// can be obtained using the query functions and the
   1.124 +  /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
   1.125 +  /// which is able to iterate on the nodes of a blossom.
   1.126    /// If the value type is integer, then the dual solution is multiplied
   1.127    /// by \ref MaxWeightedMatching::dualScale "4".
   1.128    ///
   1.129    /// \tparam GR The undirected graph type the algorithm runs on.
   1.130 -  /// \tparam WM The type edge weight map. The default type is 
   1.131 +  /// \tparam WM The type edge weight map. The default type is
   1.132    /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
   1.133  #ifdef DOXYGEN
   1.134    template <typename GR, typename WM>
   1.135 @@ -745,7 +745,7 @@
   1.136      typedef RangeMap<int> IntIntMap;
   1.137  
   1.138      enum Status {
   1.139 -      EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
   1.140 +      EVEN = -1, MATCHED = 0, ODD = 1
   1.141      };
   1.142  
   1.143      typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
   1.144 @@ -844,9 +844,6 @@
   1.145      }
   1.146  
   1.147      void destroyStructures() {
   1.148 -      _node_num = countNodes(_graph);
   1.149 -      _blossom_num = _node_num * 3 / 2;
   1.150 -
   1.151        if (_matching) {
   1.152          delete _matching;
   1.153        }
   1.154 @@ -922,10 +919,6 @@
   1.155              if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
   1.156                _delta3->push(e, rw / 2);
   1.157              }
   1.158 -          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
   1.159 -            if (_delta3->state(e) != _delta3->IN_HEAP) {
   1.160 -              _delta3->push(e, rw);
   1.161 -            }
   1.162            } else {
   1.163              typename std::map<int, Arc>::iterator it =
   1.164                (*_node_data)[vi].heap_index.find(tree);
   1.165 @@ -949,202 +942,6 @@
   1.166                    _delta2->push(vb, _blossom_set->classPrio(vb) -
   1.167                                 (*_blossom_data)[vb].offset);
   1.168                  } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
   1.169 -                           (*_blossom_data)[vb].offset){
   1.170 -                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
   1.171 -                                   (*_blossom_data)[vb].offset);
   1.172 -                }
   1.173 -              }
   1.174 -            }
   1.175 -          }
   1.176 -        }
   1.177 -      }
   1.178 -      (*_blossom_data)[blossom].offset = 0;
   1.179 -    }
   1.180 -
   1.181 -    void matchedToOdd(int blossom) {
   1.182 -      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
   1.183 -        _delta2->erase(blossom);
   1.184 -      }
   1.185 -      (*_blossom_data)[blossom].offset += _delta_sum;
   1.186 -      if (!_blossom_set->trivial(blossom)) {
   1.187 -        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
   1.188 -                     (*_blossom_data)[blossom].offset);
   1.189 -      }
   1.190 -    }
   1.191 -
   1.192 -    void evenToMatched(int blossom, int tree) {
   1.193 -      if (!_blossom_set->trivial(blossom)) {
   1.194 -        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
   1.195 -      }
   1.196 -
   1.197 -      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
   1.198 -           n != INVALID; ++n) {
   1.199 -        int ni = (*_node_index)[n];
   1.200 -        (*_node_data)[ni].pot -= _delta_sum;
   1.201 -
   1.202 -        _delta1->erase(n);
   1.203 -
   1.204 -        for (InArcIt e(_graph, n); e != INVALID; ++e) {
   1.205 -          Node v = _graph.source(e);
   1.206 -          int vb = _blossom_set->find(v);
   1.207 -          int vi = (*_node_index)[v];
   1.208 -
   1.209 -          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   1.210 -            dualScale * _weight[e];
   1.211 -
   1.212 -          if (vb == blossom) {
   1.213 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   1.214 -              _delta3->erase(e);
   1.215 -            }
   1.216 -          } else if ((*_blossom_data)[vb].status == EVEN) {
   1.217 -
   1.218 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   1.219 -              _delta3->erase(e);
   1.220 -            }
   1.221 -
   1.222 -            int vt = _tree_set->find(vb);
   1.223 -
   1.224 -            if (vt != tree) {
   1.225 -
   1.226 -              Arc r = _graph.oppositeArc(e);
   1.227 -
   1.228 -              typename std::map<int, Arc>::iterator it =
   1.229 -                (*_node_data)[ni].heap_index.find(vt);
   1.230 -
   1.231 -              if (it != (*_node_data)[ni].heap_index.end()) {
   1.232 -                if ((*_node_data)[ni].heap[it->second] > rw) {
   1.233 -                  (*_node_data)[ni].heap.replace(it->second, r);
   1.234 -                  (*_node_data)[ni].heap.decrease(r, rw);
   1.235 -                  it->second = r;
   1.236 -                }
   1.237 -              } else {
   1.238 -                (*_node_data)[ni].heap.push(r, rw);
   1.239 -                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
   1.240 -              }
   1.241 -
   1.242 -              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
   1.243 -                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
   1.244 -
   1.245 -                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
   1.246 -                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
   1.247 -                               (*_blossom_data)[blossom].offset);
   1.248 -                } else if ((*_delta2)[blossom] >
   1.249 -                           _blossom_set->classPrio(blossom) -
   1.250 -                           (*_blossom_data)[blossom].offset){
   1.251 -                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
   1.252 -                                   (*_blossom_data)[blossom].offset);
   1.253 -                }
   1.254 -              }
   1.255 -            }
   1.256 -
   1.257 -          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
   1.258 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   1.259 -              _delta3->erase(e);
   1.260 -            }
   1.261 -          } else {
   1.262 -
   1.263 -            typename std::map<int, Arc>::iterator it =
   1.264 -              (*_node_data)[vi].heap_index.find(tree);
   1.265 -
   1.266 -            if (it != (*_node_data)[vi].heap_index.end()) {
   1.267 -              (*_node_data)[vi].heap.erase(it->second);
   1.268 -              (*_node_data)[vi].heap_index.erase(it);
   1.269 -              if ((*_node_data)[vi].heap.empty()) {
   1.270 -                _blossom_set->increase(v, std::numeric_limits<Value>::max());
   1.271 -              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
   1.272 -                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
   1.273 -              }
   1.274 -
   1.275 -              if ((*_blossom_data)[vb].status == MATCHED) {
   1.276 -                if (_blossom_set->classPrio(vb) ==
   1.277 -                    std::numeric_limits<Value>::max()) {
   1.278 -                  _delta2->erase(vb);
   1.279 -                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
   1.280 -                           (*_blossom_data)[vb].offset) {
   1.281 -                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
   1.282 -                                   (*_blossom_data)[vb].offset);
   1.283 -                }
   1.284 -              }
   1.285 -            }
   1.286 -          }
   1.287 -        }
   1.288 -      }
   1.289 -    }
   1.290 -
   1.291 -    void oddToMatched(int blossom) {
   1.292 -      (*_blossom_data)[blossom].offset -= _delta_sum;
   1.293 -
   1.294 -      if (_blossom_set->classPrio(blossom) !=
   1.295 -          std::numeric_limits<Value>::max()) {
   1.296 -        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
   1.297 -                       (*_blossom_data)[blossom].offset);
   1.298 -      }
   1.299 -
   1.300 -      if (!_blossom_set->trivial(blossom)) {
   1.301 -        _delta4->erase(blossom);
   1.302 -      }
   1.303 -    }
   1.304 -
   1.305 -    void oddToEven(int blossom, int tree) {
   1.306 -      if (!_blossom_set->trivial(blossom)) {
   1.307 -        _delta4->erase(blossom);
   1.308 -        (*_blossom_data)[blossom].pot -=
   1.309 -          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
   1.310 -      }
   1.311 -
   1.312 -      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
   1.313 -           n != INVALID; ++n) {
   1.314 -        int ni = (*_node_index)[n];
   1.315 -
   1.316 -        _blossom_set->increase(n, std::numeric_limits<Value>::max());
   1.317 -
   1.318 -        (*_node_data)[ni].heap.clear();
   1.319 -        (*_node_data)[ni].heap_index.clear();
   1.320 -        (*_node_data)[ni].pot +=
   1.321 -          2 * _delta_sum - (*_blossom_data)[blossom].offset;
   1.322 -
   1.323 -        _delta1->push(n, (*_node_data)[ni].pot);
   1.324 -
   1.325 -        for (InArcIt e(_graph, n); e != INVALID; ++e) {
   1.326 -          Node v = _graph.source(e);
   1.327 -          int vb = _blossom_set->find(v);
   1.328 -          int vi = (*_node_index)[v];
   1.329 -
   1.330 -          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   1.331 -            dualScale * _weight[e];
   1.332 -
   1.333 -          if ((*_blossom_data)[vb].status == EVEN) {
   1.334 -            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
   1.335 -              _delta3->push(e, rw / 2);
   1.336 -            }
   1.337 -          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
   1.338 -            if (_delta3->state(e) != _delta3->IN_HEAP) {
   1.339 -              _delta3->push(e, rw);
   1.340 -            }
   1.341 -          } else {
   1.342 -
   1.343 -            typename std::map<int, Arc>::iterator it =
   1.344 -              (*_node_data)[vi].heap_index.find(tree);
   1.345 -
   1.346 -            if (it != (*_node_data)[vi].heap_index.end()) {
   1.347 -              if ((*_node_data)[vi].heap[it->second] > rw) {
   1.348 -                (*_node_data)[vi].heap.replace(it->second, e);
   1.349 -                (*_node_data)[vi].heap.decrease(e, rw);
   1.350 -                it->second = e;
   1.351 -              }
   1.352 -            } else {
   1.353 -              (*_node_data)[vi].heap.push(e, rw);
   1.354 -              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
   1.355 -            }
   1.356 -
   1.357 -            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
   1.358 -              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
   1.359 -
   1.360 -              if ((*_blossom_data)[vb].status == MATCHED) {
   1.361 -                if (_delta2->state(vb) != _delta2->IN_HEAP) {
   1.362 -                  _delta2->push(vb, _blossom_set->classPrio(vb) -
   1.363 -                               (*_blossom_data)[vb].offset);
   1.364 -                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
   1.365                             (*_blossom_data)[vb].offset) {
   1.366                    _delta2->decrease(vb, _blossom_set->classPrio(vb) -
   1.367                                     (*_blossom_data)[vb].offset);
   1.368 @@ -1157,43 +954,145 @@
   1.369        (*_blossom_data)[blossom].offset = 0;
   1.370      }
   1.371  
   1.372 -
   1.373 -    void matchedToUnmatched(int blossom) {
   1.374 +    void matchedToOdd(int blossom) {
   1.375        if (_delta2->state(blossom) == _delta2->IN_HEAP) {
   1.376          _delta2->erase(blossom);
   1.377        }
   1.378 +      (*_blossom_data)[blossom].offset += _delta_sum;
   1.379 +      if (!_blossom_set->trivial(blossom)) {
   1.380 +        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
   1.381 +                      (*_blossom_data)[blossom].offset);
   1.382 +      }
   1.383 +    }
   1.384 +
   1.385 +    void evenToMatched(int blossom, int tree) {
   1.386 +      if (!_blossom_set->trivial(blossom)) {
   1.387 +        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
   1.388 +      }
   1.389  
   1.390        for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
   1.391             n != INVALID; ++n) {
   1.392          int ni = (*_node_index)[n];
   1.393 -
   1.394 -        _blossom_set->increase(n, std::numeric_limits<Value>::max());
   1.395 -
   1.396 -        (*_node_data)[ni].heap.clear();
   1.397 -        (*_node_data)[ni].heap_index.clear();
   1.398 -
   1.399 -        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
   1.400 -          Node v = _graph.target(e);
   1.401 +        (*_node_data)[ni].pot -= _delta_sum;
   1.402 +
   1.403 +        _delta1->erase(n);
   1.404 +
   1.405 +        for (InArcIt e(_graph, n); e != INVALID; ++e) {
   1.406 +          Node v = _graph.source(e);
   1.407            int vb = _blossom_set->find(v);
   1.408            int vi = (*_node_index)[v];
   1.409  
   1.410            Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   1.411              dualScale * _weight[e];
   1.412  
   1.413 -          if ((*_blossom_data)[vb].status == EVEN) {
   1.414 -            if (_delta3->state(e) != _delta3->IN_HEAP) {
   1.415 -              _delta3->push(e, rw);
   1.416 +          if (vb == blossom) {
   1.417 +            if (_delta3->state(e) == _delta3->IN_HEAP) {
   1.418 +              _delta3->erase(e);
   1.419 +            }
   1.420 +          } else if ((*_blossom_data)[vb].status == EVEN) {
   1.421 +
   1.422 +            if (_delta3->state(e) == _delta3->IN_HEAP) {
   1.423 +              _delta3->erase(e);
   1.424 +            }
   1.425 +
   1.426 +            int vt = _tree_set->find(vb);
   1.427 +
   1.428 +            if (vt != tree) {
   1.429 +
   1.430 +              Arc r = _graph.oppositeArc(e);
   1.431 +
   1.432 +              typename std::map<int, Arc>::iterator it =
   1.433 +                (*_node_data)[ni].heap_index.find(vt);
   1.434 +
   1.435 +              if (it != (*_node_data)[ni].heap_index.end()) {
   1.436 +                if ((*_node_data)[ni].heap[it->second] > rw) {
   1.437 +                  (*_node_data)[ni].heap.replace(it->second, r);
   1.438 +                  (*_node_data)[ni].heap.decrease(r, rw);
   1.439 +                  it->second = r;
   1.440 +                }
   1.441 +              } else {
   1.442 +                (*_node_data)[ni].heap.push(r, rw);
   1.443 +                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
   1.444 +              }
   1.445 +
   1.446 +              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
   1.447 +                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
   1.448 +
   1.449 +                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
   1.450 +                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
   1.451 +                               (*_blossom_data)[blossom].offset);
   1.452 +                } else if ((*_delta2)[blossom] >
   1.453 +                           _blossom_set->classPrio(blossom) -
   1.454 +                           (*_blossom_data)[blossom].offset){
   1.455 +                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
   1.456 +                                   (*_blossom_data)[blossom].offset);
   1.457 +                }
   1.458 +              }
   1.459 +            }
   1.460 +          } else {
   1.461 +
   1.462 +            typename std::map<int, Arc>::iterator it =
   1.463 +              (*_node_data)[vi].heap_index.find(tree);
   1.464 +
   1.465 +            if (it != (*_node_data)[vi].heap_index.end()) {
   1.466 +              (*_node_data)[vi].heap.erase(it->second);
   1.467 +              (*_node_data)[vi].heap_index.erase(it);
   1.468 +              if ((*_node_data)[vi].heap.empty()) {
   1.469 +                _blossom_set->increase(v, std::numeric_limits<Value>::max());
   1.470 +              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
   1.471 +                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
   1.472 +              }
   1.473 +
   1.474 +              if ((*_blossom_data)[vb].status == MATCHED) {
   1.475 +                if (_blossom_set->classPrio(vb) ==
   1.476 +                    std::numeric_limits<Value>::max()) {
   1.477 +                  _delta2->erase(vb);
   1.478 +                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
   1.479 +                           (*_blossom_data)[vb].offset) {
   1.480 +                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
   1.481 +                                   (*_blossom_data)[vb].offset);
   1.482 +                }
   1.483 +              }
   1.484              }
   1.485            }
   1.486          }
   1.487        }
   1.488      }
   1.489  
   1.490 -    void unmatchedToMatched(int blossom) {
   1.491 +    void oddToMatched(int blossom) {
   1.492 +      (*_blossom_data)[blossom].offset -= _delta_sum;
   1.493 +
   1.494 +      if (_blossom_set->classPrio(blossom) !=
   1.495 +          std::numeric_limits<Value>::max()) {
   1.496 +        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
   1.497 +                      (*_blossom_data)[blossom].offset);
   1.498 +      }
   1.499 +
   1.500 +      if (!_blossom_set->trivial(blossom)) {
   1.501 +        _delta4->erase(blossom);
   1.502 +      }
   1.503 +    }
   1.504 +
   1.505 +    void oddToEven(int blossom, int tree) {
   1.506 +      if (!_blossom_set->trivial(blossom)) {
   1.507 +        _delta4->erase(blossom);
   1.508 +        (*_blossom_data)[blossom].pot -=
   1.509 +          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
   1.510 +      }
   1.511 +
   1.512        for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
   1.513             n != INVALID; ++n) {
   1.514          int ni = (*_node_index)[n];
   1.515  
   1.516 +        _blossom_set->increase(n, std::numeric_limits<Value>::max());
   1.517 +
   1.518 +        (*_node_data)[ni].heap.clear();
   1.519 +        (*_node_data)[ni].heap_index.clear();
   1.520 +        (*_node_data)[ni].pot +=
   1.521 +          2 * _delta_sum - (*_blossom_data)[blossom].offset;
   1.522 +
   1.523 +        _delta1->push(n, (*_node_data)[ni].pot);
   1.524 +
   1.525          for (InArcIt e(_graph, n); e != INVALID; ++e) {
   1.526            Node v = _graph.source(e);
   1.527            int vb = _blossom_set->find(v);
   1.528 @@ -1202,54 +1101,44 @@
   1.529            Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   1.530              dualScale * _weight[e];
   1.531  
   1.532 -          if (vb == blossom) {
   1.533 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   1.534 -              _delta3->erase(e);
   1.535 +          if ((*_blossom_data)[vb].status == EVEN) {
   1.536 +            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
   1.537 +              _delta3->push(e, rw / 2);
   1.538              }
   1.539 -          } else if ((*_blossom_data)[vb].status == EVEN) {
   1.540 -
   1.541 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   1.542 -              _delta3->erase(e);
   1.543 -            }
   1.544 -
   1.545 -            int vt = _tree_set->find(vb);
   1.546 -
   1.547 -            Arc r = _graph.oppositeArc(e);
   1.548 +          } else {
   1.549  
   1.550              typename std::map<int, Arc>::iterator it =
   1.551 -              (*_node_data)[ni].heap_index.find(vt);
   1.552 -
   1.553 -            if (it != (*_node_data)[ni].heap_index.end()) {
   1.554 -              if ((*_node_data)[ni].heap[it->second] > rw) {
   1.555 -                (*_node_data)[ni].heap.replace(it->second, r);
   1.556 -                (*_node_data)[ni].heap.decrease(r, rw);
   1.557 -                it->second = r;
   1.558 +              (*_node_data)[vi].heap_index.find(tree);
   1.559 +
   1.560 +            if (it != (*_node_data)[vi].heap_index.end()) {
   1.561 +              if ((*_node_data)[vi].heap[it->second] > rw) {
   1.562 +                (*_node_data)[vi].heap.replace(it->second, e);
   1.563 +                (*_node_data)[vi].heap.decrease(e, rw);
   1.564 +                it->second = e;
   1.565                }
   1.566              } else {
   1.567 -              (*_node_data)[ni].heap.push(r, rw);
   1.568 -              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
   1.569 +              (*_node_data)[vi].heap.push(e, rw);
   1.570 +              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
   1.571              }
   1.572  
   1.573 -            if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
   1.574 -              _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
   1.575 -
   1.576 -              if (_delta2->state(blossom) != _delta2->IN_HEAP) {
   1.577 -                _delta2->push(blossom, _blossom_set->classPrio(blossom) -
   1.578 -                             (*_blossom_data)[blossom].offset);
   1.579 -              } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
   1.580 -                         (*_blossom_data)[blossom].offset){
   1.581 -                _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
   1.582 -                                 (*_blossom_data)[blossom].offset);
   1.583 +            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
   1.584 +              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
   1.585 +
   1.586 +              if ((*_blossom_data)[vb].status == MATCHED) {
   1.587 +                if (_delta2->state(vb) != _delta2->IN_HEAP) {
   1.588 +                  _delta2->push(vb, _blossom_set->classPrio(vb) -
   1.589 +                               (*_blossom_data)[vb].offset);
   1.590 +                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
   1.591 +                           (*_blossom_data)[vb].offset) {
   1.592 +                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
   1.593 +                                   (*_blossom_data)[vb].offset);
   1.594 +                }
   1.595                }
   1.596              }
   1.597 -
   1.598 -          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
   1.599 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   1.600 -              _delta3->erase(e);
   1.601 -            }
   1.602            }
   1.603          }
   1.604        }
   1.605 +      (*_blossom_data)[blossom].offset = 0;
   1.606      }
   1.607  
   1.608      void alternatePath(int even, int tree) {
   1.609 @@ -1294,39 +1183,42 @@
   1.610        alternatePath(blossom, tree);
   1.611        destroyTree(tree);
   1.612  
   1.613 -      (*_blossom_data)[blossom].status = UNMATCHED;
   1.614        (*_blossom_data)[blossom].base = node;
   1.615 -      matchedToUnmatched(blossom);
   1.616 +      (*_blossom_data)[blossom].next = INVALID;
   1.617      }
   1.618  
   1.619 -
   1.620      void augmentOnEdge(const Edge& edge) {
   1.621  
   1.622        int left = _blossom_set->find(_graph.u(edge));
   1.623        int right = _blossom_set->find(_graph.v(edge));
   1.624  
   1.625 -      if ((*_blossom_data)[left].status == EVEN) {
   1.626 -        int left_tree = _tree_set->find(left);
   1.627 -        alternatePath(left, left_tree);
   1.628 -        destroyTree(left_tree);
   1.629 -      } else {
   1.630 -        (*_blossom_data)[left].status = MATCHED;
   1.631 -        unmatchedToMatched(left);
   1.632 -      }
   1.633 -
   1.634 -      if ((*_blossom_data)[right].status == EVEN) {
   1.635 -        int right_tree = _tree_set->find(right);
   1.636 -        alternatePath(right, right_tree);
   1.637 -        destroyTree(right_tree);
   1.638 -      } else {
   1.639 -        (*_blossom_data)[right].status = MATCHED;
   1.640 -        unmatchedToMatched(right);
   1.641 -      }
   1.642 +      int left_tree = _tree_set->find(left);
   1.643 +      alternatePath(left, left_tree);
   1.644 +      destroyTree(left_tree);
   1.645 +
   1.646 +      int right_tree = _tree_set->find(right);
   1.647 +      alternatePath(right, right_tree);
   1.648 +      destroyTree(right_tree);
   1.649  
   1.650        (*_blossom_data)[left].next = _graph.direct(edge, true);
   1.651        (*_blossom_data)[right].next = _graph.direct(edge, false);
   1.652      }
   1.653  
   1.654 +    void augmentOnArc(const Arc& arc) {
   1.655 +
   1.656 +      int left = _blossom_set->find(_graph.source(arc));
   1.657 +      int right = _blossom_set->find(_graph.target(arc));
   1.658 +
   1.659 +      (*_blossom_data)[left].status = MATCHED;
   1.660 +
   1.661 +      int right_tree = _tree_set->find(right);
   1.662 +      alternatePath(right, right_tree);
   1.663 +      destroyTree(right_tree);
   1.664 +
   1.665 +      (*_blossom_data)[left].next = arc;
   1.666 +      (*_blossom_data)[right].next = _graph.oppositeArc(arc);
   1.667 +    }
   1.668 +
   1.669      void extendOnArc(const Arc& arc) {
   1.670        int base = _blossom_set->find(_graph.target(arc));
   1.671        int tree = _tree_set->find(base);
   1.672 @@ -1529,7 +1421,7 @@
   1.673            _tree_set->insert(sb, tree);
   1.674            (*_blossom_data)[sb].pred = pred;
   1.675            (*_blossom_data)[sb].next =
   1.676 -                           _graph.oppositeArc((*_blossom_data)[tb].next);
   1.677 +            _graph.oppositeArc((*_blossom_data)[tb].next);
   1.678  
   1.679            pred = (*_blossom_data)[ub].next;
   1.680  
   1.681 @@ -1629,7 +1521,7 @@
   1.682        }
   1.683  
   1.684        for (int i = 0; i < int(blossoms.size()); ++i) {
   1.685 -        if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
   1.686 +        if ((*_blossom_data)[blossoms[i]].next != INVALID) {
   1.687  
   1.688            Value offset = (*_blossom_data)[blossoms[i]].offset;
   1.689            (*_blossom_data)[blossoms[i]].pot += 2 * offset;
   1.690 @@ -1757,12 +1649,11 @@
   1.691          Value d4 = !_delta4->empty() ?
   1.692            _delta4->prio() : std::numeric_limits<Value>::max();
   1.693  
   1.694 -        _delta_sum = d1; OpType ot = D1;
   1.695 +        _delta_sum = d3; OpType ot = D3;
   1.696 +        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
   1.697          if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
   1.698 -        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
   1.699          if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
   1.700  
   1.701 -
   1.702          switch (ot) {
   1.703          case D1:
   1.704            {
   1.705 @@ -1775,8 +1666,13 @@
   1.706            {
   1.707              int blossom = _delta2->top();
   1.708              Node n = _blossom_set->classTop(blossom);
   1.709 -            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
   1.710 -            extendOnArc(e);
   1.711 +            Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
   1.712 +            if ((*_blossom_data)[blossom].next == INVALID) {
   1.713 +              augmentOnArc(a);
   1.714 +              --unmatched;
   1.715 +            } else {
   1.716 +              extendOnArc(a);
   1.717 +            }
   1.718            }
   1.719            break;
   1.720          case D3:
   1.721 @@ -1789,20 +1685,8 @@
   1.722              if (left_blossom == right_blossom) {
   1.723                _delta3->pop();
   1.724              } else {
   1.725 -              int left_tree;
   1.726 -              if ((*_blossom_data)[left_blossom].status == EVEN) {
   1.727 -                left_tree = _tree_set->find(left_blossom);
   1.728 -              } else {
   1.729 -                left_tree = -1;
   1.730 -                ++unmatched;
   1.731 -              }
   1.732 -              int right_tree;
   1.733 -              if ((*_blossom_data)[right_blossom].status == EVEN) {
   1.734 -                right_tree = _tree_set->find(right_blossom);
   1.735 -              } else {
   1.736 -                right_tree = -1;
   1.737 -                ++unmatched;
   1.738 -              }
   1.739 +              int left_tree = _tree_set->find(left_blossom);
   1.740 +              int right_tree = _tree_set->find(right_blossom);
   1.741  
   1.742                if (left_tree == right_tree) {
   1.743                  shrinkOnEdge(e, left_tree);
   1.744 @@ -1837,7 +1721,7 @@
   1.745      /// @}
   1.746  
   1.747      /// \name Primal Solution
   1.748 -    /// Functions to get the primal solution, i.e. the maximum weighted 
   1.749 +    /// Functions to get the primal solution, i.e. the maximum weighted
   1.750      /// matching.\n
   1.751      /// Either \ref run() or \ref start() function should be called before
   1.752      /// using them.
   1.753 @@ -1856,7 +1740,7 @@
   1.754            sum += _weight[(*_matching)[n]];
   1.755          }
   1.756        }
   1.757 -      return sum /= 2;
   1.758 +      return sum / 2;
   1.759      }
   1.760  
   1.761      /// \brief Return the size (cardinality) of the matching.
   1.762 @@ -1876,7 +1760,7 @@
   1.763  
   1.764      /// \brief Return \c true if the given edge is in the matching.
   1.765      ///
   1.766 -    /// This function returns \c true if the given edge is in the found 
   1.767 +    /// This function returns \c true if the given edge is in the found
   1.768      /// matching.
   1.769      ///
   1.770      /// \pre Either run() or start() must be called before using this function.
   1.771 @@ -1887,7 +1771,7 @@
   1.772      /// \brief Return the matching arc (or edge) incident to the given node.
   1.773      ///
   1.774      /// This function returns the matching arc (or edge) incident to the
   1.775 -    /// given node in the found matching or \c INVALID if the node is 
   1.776 +    /// given node in the found matching or \c INVALID if the node is
   1.777      /// not covered by the matching.
   1.778      ///
   1.779      /// \pre Either run() or start() must be called before using this function.
   1.780 @@ -1905,7 +1789,7 @@
   1.781  
   1.782      /// \brief Return the mate of the given node.
   1.783      ///
   1.784 -    /// This function returns the mate of the given node in the found 
   1.785 +    /// This function returns the mate of the given node in the found
   1.786      /// matching or \c INVALID if the node is not covered by the matching.
   1.787      ///
   1.788      /// \pre Either run() or start() must be called before using this function.
   1.789 @@ -1925,8 +1809,8 @@
   1.790  
   1.791      /// \brief Return the value of the dual solution.
   1.792      ///
   1.793 -    /// This function returns the value of the dual solution. 
   1.794 -    /// It should be equal to the primal value scaled by \ref dualScale 
   1.795 +    /// This function returns the value of the dual solution.
   1.796 +    /// It should be equal to the primal value scaled by \ref dualScale
   1.797      /// "dual scale".
   1.798      ///
   1.799      /// \pre Either run() or start() must be called before using this function.
   1.800 @@ -1981,9 +1865,9 @@
   1.801  
   1.802      /// \brief Iterator for obtaining the nodes of a blossom.
   1.803      ///
   1.804 -    /// This class provides an iterator for obtaining the nodes of the 
   1.805 +    /// This class provides an iterator for obtaining the nodes of the
   1.806      /// given blossom. It lists a subset of the nodes.
   1.807 -    /// Before using this iterator, you must allocate a 
   1.808 +    /// Before using this iterator, you must allocate a
   1.809      /// MaxWeightedMatching class and execute it.
   1.810      class BlossomIt {
   1.811      public:
   1.812 @@ -1992,8 +1876,8 @@
   1.813        ///
   1.814        /// Constructor to get the nodes of the given variable.
   1.815        ///
   1.816 -      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or 
   1.817 -      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be 
   1.818 +      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
   1.819 +      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
   1.820        /// called before initializing this iterator.
   1.821        BlossomIt(const MaxWeightedMatching& algorithm, int variable)
   1.822          : _algorithm(&algorithm)
   1.823 @@ -2046,8 +1930,8 @@
   1.824    /// is based on extensive use of priority queues and provides
   1.825    /// \f$O(nm\log n)\f$ time complexity.
   1.826    ///
   1.827 -  /// The maximum weighted perfect matching problem is to find a subset of 
   1.828 -  /// the edges in an undirected graph with maximum overall weight for which 
   1.829 +  /// The maximum weighted perfect matching problem is to find a subset of
   1.830 +  /// the edges in an undirected graph with maximum overall weight for which
   1.831    /// each node has exactly one incident edge.
   1.832    /// It can be formulated with the following linear program.
   1.833    /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
   1.834 @@ -2070,16 +1954,16 @@
   1.835    /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
   1.836        \frac{\vert B \vert - 1}{2}z_B\f] */
   1.837    ///
   1.838 -  /// The algorithm can be executed with the run() function. 
   1.839 +  /// The algorithm can be executed with the run() function.
   1.840    /// After it the matching (the primal solution) and the dual solution
   1.841 -  /// can be obtained using the query functions and the 
   1.842 -  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, 
   1.843 -  /// which is able to iterate on the nodes of a blossom. 
   1.844 +  /// can be obtained using the query functions and the
   1.845 +  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
   1.846 +  /// which is able to iterate on the nodes of a blossom.
   1.847    /// If the value type is integer, then the dual solution is multiplied
   1.848    /// by \ref MaxWeightedMatching::dualScale "4".
   1.849    ///
   1.850    /// \tparam GR The undirected graph type the algorithm runs on.
   1.851 -  /// \tparam WM The type edge weight map. The default type is 
   1.852 +  /// \tparam WM The type edge weight map. The default type is
   1.853    /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
   1.854  #ifdef DOXYGEN
   1.855    template <typename GR, typename WM>
   1.856 @@ -2233,9 +2117,6 @@
   1.857      }
   1.858  
   1.859      void destroyStructures() {
   1.860 -      _node_num = countNodes(_graph);
   1.861 -      _blossom_num = _node_num * 3 / 2;
   1.862 -
   1.863        if (_matching) {
   1.864          delete _matching;
   1.865        }
   1.866 @@ -2991,8 +2872,8 @@
   1.867          Value d4 = !_delta4->empty() ?
   1.868            _delta4->prio() : std::numeric_limits<Value>::max();
   1.869  
   1.870 -        _delta_sum = d2; OpType ot = D2;
   1.871 -        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
   1.872 +        _delta_sum = d3; OpType ot = D3;
   1.873 +        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
   1.874          if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
   1.875  
   1.876          if (_delta_sum == std::numeric_limits<Value>::max()) {
   1.877 @@ -3055,7 +2936,7 @@
   1.878      /// @}
   1.879  
   1.880      /// \name Primal Solution
   1.881 -    /// Functions to get the primal solution, i.e. the maximum weighted 
   1.882 +    /// Functions to get the primal solution, i.e. the maximum weighted
   1.883      /// perfect matching.\n
   1.884      /// Either \ref run() or \ref start() function should be called before
   1.885      /// using them.
   1.886 @@ -3074,12 +2955,12 @@
   1.887            sum += _weight[(*_matching)[n]];
   1.888          }
   1.889        }
   1.890 -      return sum /= 2;
   1.891 +      return sum / 2;
   1.892      }
   1.893  
   1.894      /// \brief Return \c true if the given edge is in the matching.
   1.895      ///
   1.896 -    /// This function returns \c true if the given edge is in the found 
   1.897 +    /// This function returns \c true if the given edge is in the found
   1.898      /// matching.
   1.899      ///
   1.900      /// \pre Either run() or start() must be called before using this function.
   1.901 @@ -3090,7 +2971,7 @@
   1.902      /// \brief Return the matching arc (or edge) incident to the given node.
   1.903      ///
   1.904      /// This function returns the matching arc (or edge) incident to the
   1.905 -    /// given node in the found matching or \c INVALID if the node is 
   1.906 +    /// given node in the found matching or \c INVALID if the node is
   1.907      /// not covered by the matching.
   1.908      ///
   1.909      /// \pre Either run() or start() must be called before using this function.
   1.910 @@ -3108,7 +2989,7 @@
   1.911  
   1.912      /// \brief Return the mate of the given node.
   1.913      ///
   1.914 -    /// This function returns the mate of the given node in the found 
   1.915 +    /// This function returns the mate of the given node in the found
   1.916      /// matching or \c INVALID if the node is not covered by the matching.
   1.917      ///
   1.918      /// \pre Either run() or start() must be called before using this function.
   1.919 @@ -3127,8 +3008,8 @@
   1.920  
   1.921      /// \brief Return the value of the dual solution.
   1.922      ///
   1.923 -    /// This function returns the value of the dual solution. 
   1.924 -    /// It should be equal to the primal value scaled by \ref dualScale 
   1.925 +    /// This function returns the value of the dual solution.
   1.926 +    /// It should be equal to the primal value scaled by \ref dualScale
   1.927      /// "dual scale".
   1.928      ///
   1.929      /// \pre Either run() or start() must be called before using this function.
   1.930 @@ -3183,9 +3064,9 @@
   1.931  
   1.932      /// \brief Iterator for obtaining the nodes of a blossom.
   1.933      ///
   1.934 -    /// This class provides an iterator for obtaining the nodes of the 
   1.935 +    /// This class provides an iterator for obtaining the nodes of the
   1.936      /// given blossom. It lists a subset of the nodes.
   1.937 -    /// Before using this iterator, you must allocate a 
   1.938 +    /// Before using this iterator, you must allocate a
   1.939      /// MaxWeightedPerfectMatching class and execute it.
   1.940      class BlossomIt {
   1.941      public:
   1.942 @@ -3194,8 +3075,8 @@
   1.943        ///
   1.944        /// Constructor to get the nodes of the given variable.
   1.945        ///
   1.946 -      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" 
   1.947 -      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" 
   1.948 +      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
   1.949 +      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
   1.950        /// must be called before initializing this iterator.
   1.951        BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
   1.952          : _algorithm(&algorithm)
   1.953 @@ -3241,4 +3122,4 @@
   1.954  
   1.955  } //END OF NAMESPACE LEMON
   1.956  
   1.957 -#endif //LEMON_MAX_MATCHING_H
   1.958 +#endif //LEMON_MATCHING_H