author Balazs Dezso Sun, 20 Sep 2009 21:38:24 +0200 changeset 947 0513ccfea967 parent 759 6d5f547e5bfb child 948 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 file | annotate | diff | comparison | revisions
     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