lemon/matching.h
changeset 874 d8ea85825e02
parent 868 0513ccfea967
child 875 07ec2b52e53d
     1.1 --- a/lemon/matching.h	Tue Mar 16 21:18:39 2010 +0100
     1.2 +++ b/lemon/matching.h	Tue Mar 16 21:27:35 2010 +0100
     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 @@ -28,6 +28,7 @@
    1.15  #include <lemon/unionfind.h>
    1.16  #include <lemon/bin_heap.h>
    1.17  #include <lemon/maps.h>
    1.18 +#include <lemon/fractional_matching.h>
    1.19  
    1.20  ///\ingroup matching
    1.21  ///\file
    1.22 @@ -41,7 +42,7 @@
    1.23    ///
    1.24    /// This class implements Edmonds' alternating forest matching algorithm
    1.25    /// for finding a maximum cardinality matching in a general undirected graph.
    1.26 -  /// It can be started from an arbitrary initial matching 
    1.27 +  /// It can be started from an arbitrary initial matching
    1.28    /// (the default is the empty one).
    1.29    ///
    1.30    /// The dual solution of the problem is a map of the nodes to
    1.31 @@ -69,11 +70,11 @@
    1.32  
    1.33      ///\brief Status constants for Gallai-Edmonds decomposition.
    1.34      ///
    1.35 -    ///These constants are used for indicating the Gallai-Edmonds 
    1.36 +    ///These constants are used for indicating the Gallai-Edmonds
    1.37      ///decomposition of a graph. The nodes with status \c EVEN (or \c D)
    1.38      ///induce a subgraph with factor-critical components, the nodes with
    1.39      ///status \c ODD (or \c A) form the canonical barrier, and the nodes
    1.40 -    ///with status \c MATCHED (or \c C) induce a subgraph having a 
    1.41 +    ///with status \c MATCHED (or \c C) induce a subgraph having a
    1.42      ///perfect matching.
    1.43      enum Status {
    1.44        EVEN = 1,       ///< = 1. (\c D is an alias for \c EVEN.)
    1.45 @@ -512,7 +513,7 @@
    1.46        }
    1.47      }
    1.48  
    1.49 -    /// \brief Start Edmonds' algorithm with a heuristic improvement 
    1.50 +    /// \brief Start Edmonds' algorithm with a heuristic improvement
    1.51      /// for dense graphs
    1.52      ///
    1.53      /// This function runs Edmonds' algorithm with a heuristic of postponing
    1.54 @@ -534,8 +535,8 @@
    1.55  
    1.56      /// \brief Run Edmonds' algorithm
    1.57      ///
    1.58 -    /// This function runs Edmonds' algorithm. An additional heuristic of 
    1.59 -    /// postponing shrinks is used for relatively dense graphs 
    1.60 +    /// This function runs Edmonds' algorithm. An additional heuristic of
    1.61 +    /// postponing shrinks is used for relatively dense graphs
    1.62      /// (for which <tt>m>=2*n</tt> holds).
    1.63      void run() {
    1.64        if (countEdges(_graph) < 2 * countNodes(_graph)) {
    1.65 @@ -556,7 +557,7 @@
    1.66  
    1.67      /// \brief Return the size (cardinality) of the matching.
    1.68      ///
    1.69 -    /// This function returns the size (cardinality) of the current matching. 
    1.70 +    /// This function returns the size (cardinality) of the current matching.
    1.71      /// After run() it returns the size of the maximum matching in the graph.
    1.72      int matchingSize() const {
    1.73        int size = 0;
    1.74 @@ -570,7 +571,7 @@
    1.75  
    1.76      /// \brief Return \c true if the given edge is in the matching.
    1.77      ///
    1.78 -    /// This function returns \c true if the given edge is in the current 
    1.79 +    /// This function returns \c true if the given edge is in the current
    1.80      /// matching.
    1.81      bool matching(const Edge& edge) const {
    1.82        return edge == (*_matching)[_graph.u(edge)];
    1.83 @@ -579,7 +580,7 @@
    1.84      /// \brief Return the matching arc (or edge) incident to the given node.
    1.85      ///
    1.86      /// This function returns the matching arc (or edge) incident to the
    1.87 -    /// given node in the current matching or \c INVALID if the node is 
    1.88 +    /// given node in the current matching or \c INVALID if the node is
    1.89      /// not covered by the matching.
    1.90      Arc matching(const Node& n) const {
    1.91        return (*_matching)[n];
    1.92 @@ -595,7 +596,7 @@
    1.93  
    1.94      /// \brief Return the mate of the given node.
    1.95      ///
    1.96 -    /// This function returns the mate of the given node in the current 
    1.97 +    /// This function returns the mate of the given node in the current
    1.98      /// matching or \c INVALID if the node is not covered by the matching.
    1.99      Node mate(const Node& n) const {
   1.100        return (*_matching)[n] != INVALID ?
   1.101 @@ -605,7 +606,7 @@
   1.102      /// @}
   1.103  
   1.104      /// \name Dual Solution
   1.105 -    /// Functions to get the dual solution, i.e. the Gallai-Edmonds 
   1.106 +    /// Functions to get the dual solution, i.e. the Gallai-Edmonds
   1.107      /// decomposition.
   1.108  
   1.109      /// @{
   1.110 @@ -648,8 +649,8 @@
   1.111    /// on extensive use of priority queues and provides
   1.112    /// \f$O(nm\log n)\f$ time complexity.
   1.113    ///
   1.114 -  /// The maximum weighted matching problem is to find a subset of the 
   1.115 -  /// edges in an undirected graph with maximum overall weight for which 
   1.116 +  /// The maximum weighted matching problem is to find a subset of the
   1.117 +  /// edges in an undirected graph with maximum overall weight for which
   1.118    /// each node has at most one incident edge.
   1.119    /// It can be formulated with the following linear program.
   1.120    /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
   1.121 @@ -673,16 +674,16 @@
   1.122    /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
   1.123        \frac{\vert B \vert - 1}{2}z_B\f] */
   1.124    ///
   1.125 -  /// The algorithm can be executed with the run() function. 
   1.126 +  /// The algorithm can be executed with the run() function.
   1.127    /// After it the matching (the primal solution) and the dual solution
   1.128 -  /// can be obtained using the query functions and the 
   1.129 -  /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class, 
   1.130 -  /// which is able to iterate on the nodes of a blossom. 
   1.131 +  /// can be obtained using the query functions and the
   1.132 +  /// \ref MaxWeightedMatching::BlossomIt "BlossomIt" nested class,
   1.133 +  /// which is able to iterate on the nodes of a blossom.
   1.134    /// If the value type is integer, then the dual solution is multiplied
   1.135    /// by \ref MaxWeightedMatching::dualScale "4".
   1.136    ///
   1.137    /// \tparam GR The undirected graph type the algorithm runs on.
   1.138 -  /// \tparam WM The type edge weight map. The default type is 
   1.139 +  /// \tparam WM The type edge weight map. The default type is
   1.140    /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
   1.141  #ifdef DOXYGEN
   1.142    template <typename GR, typename WM>
   1.143 @@ -745,7 +746,7 @@
   1.144      typedef RangeMap<int> IntIntMap;
   1.145  
   1.146      enum Status {
   1.147 -      EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
   1.148 +      EVEN = -1, MATCHED = 0, ODD = 1
   1.149      };
   1.150  
   1.151      typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
   1.152 @@ -797,6 +798,10 @@
   1.153      BinHeap<Value, IntIntMap> *_delta4;
   1.154  
   1.155      Value _delta_sum;
   1.156 +    int _unmatched;
   1.157 +
   1.158 +    typedef MaxWeightedFractionalMatching<Graph, WeightMap> FractionalMatching;
   1.159 +    FractionalMatching *_fractional;
   1.160  
   1.161      void createStructures() {
   1.162        _node_num = countNodes(_graph);
   1.163 @@ -844,9 +849,6 @@
   1.164      }
   1.165  
   1.166      void destroyStructures() {
   1.167 -      _node_num = countNodes(_graph);
   1.168 -      _blossom_num = _node_num * 3 / 2;
   1.169 -
   1.170        if (_matching) {
   1.171          delete _matching;
   1.172        }
   1.173 @@ -922,10 +924,6 @@
   1.174              if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
   1.175                _delta3->push(e, rw / 2);
   1.176              }
   1.177 -          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
   1.178 -            if (_delta3->state(e) != _delta3->IN_HEAP) {
   1.179 -              _delta3->push(e, rw);
   1.180 -            }
   1.181            } else {
   1.182              typename std::map<int, Arc>::iterator it =
   1.183                (*_node_data)[vi].heap_index.find(tree);
   1.184 @@ -949,202 +947,6 @@
   1.185                    _delta2->push(vb, _blossom_set->classPrio(vb) -
   1.186                                 (*_blossom_data)[vb].offset);
   1.187                  } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
   1.188 -                           (*_blossom_data)[vb].offset){
   1.189 -                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
   1.190 -                                   (*_blossom_data)[vb].offset);
   1.191 -                }
   1.192 -              }
   1.193 -            }
   1.194 -          }
   1.195 -        }
   1.196 -      }
   1.197 -      (*_blossom_data)[blossom].offset = 0;
   1.198 -    }
   1.199 -
   1.200 -    void matchedToOdd(int blossom) {
   1.201 -      if (_delta2->state(blossom) == _delta2->IN_HEAP) {
   1.202 -        _delta2->erase(blossom);
   1.203 -      }
   1.204 -      (*_blossom_data)[blossom].offset += _delta_sum;
   1.205 -      if (!_blossom_set->trivial(blossom)) {
   1.206 -        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
   1.207 -                     (*_blossom_data)[blossom].offset);
   1.208 -      }
   1.209 -    }
   1.210 -
   1.211 -    void evenToMatched(int blossom, int tree) {
   1.212 -      if (!_blossom_set->trivial(blossom)) {
   1.213 -        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
   1.214 -      }
   1.215 -
   1.216 -      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
   1.217 -           n != INVALID; ++n) {
   1.218 -        int ni = (*_node_index)[n];
   1.219 -        (*_node_data)[ni].pot -= _delta_sum;
   1.220 -
   1.221 -        _delta1->erase(n);
   1.222 -
   1.223 -        for (InArcIt e(_graph, n); e != INVALID; ++e) {
   1.224 -          Node v = _graph.source(e);
   1.225 -          int vb = _blossom_set->find(v);
   1.226 -          int vi = (*_node_index)[v];
   1.227 -
   1.228 -          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   1.229 -            dualScale * _weight[e];
   1.230 -
   1.231 -          if (vb == blossom) {
   1.232 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   1.233 -              _delta3->erase(e);
   1.234 -            }
   1.235 -          } else if ((*_blossom_data)[vb].status == EVEN) {
   1.236 -
   1.237 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   1.238 -              _delta3->erase(e);
   1.239 -            }
   1.240 -
   1.241 -            int vt = _tree_set->find(vb);
   1.242 -
   1.243 -            if (vt != tree) {
   1.244 -
   1.245 -              Arc r = _graph.oppositeArc(e);
   1.246 -
   1.247 -              typename std::map<int, Arc>::iterator it =
   1.248 -                (*_node_data)[ni].heap_index.find(vt);
   1.249 -
   1.250 -              if (it != (*_node_data)[ni].heap_index.end()) {
   1.251 -                if ((*_node_data)[ni].heap[it->second] > rw) {
   1.252 -                  (*_node_data)[ni].heap.replace(it->second, r);
   1.253 -                  (*_node_data)[ni].heap.decrease(r, rw);
   1.254 -                  it->second = r;
   1.255 -                }
   1.256 -              } else {
   1.257 -                (*_node_data)[ni].heap.push(r, rw);
   1.258 -                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
   1.259 -              }
   1.260 -
   1.261 -              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
   1.262 -                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
   1.263 -
   1.264 -                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
   1.265 -                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
   1.266 -                               (*_blossom_data)[blossom].offset);
   1.267 -                } else if ((*_delta2)[blossom] >
   1.268 -                           _blossom_set->classPrio(blossom) -
   1.269 -                           (*_blossom_data)[blossom].offset){
   1.270 -                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
   1.271 -                                   (*_blossom_data)[blossom].offset);
   1.272 -                }
   1.273 -              }
   1.274 -            }
   1.275 -
   1.276 -          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
   1.277 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   1.278 -              _delta3->erase(e);
   1.279 -            }
   1.280 -          } else {
   1.281 -
   1.282 -            typename std::map<int, Arc>::iterator it =
   1.283 -              (*_node_data)[vi].heap_index.find(tree);
   1.284 -
   1.285 -            if (it != (*_node_data)[vi].heap_index.end()) {
   1.286 -              (*_node_data)[vi].heap.erase(it->second);
   1.287 -              (*_node_data)[vi].heap_index.erase(it);
   1.288 -              if ((*_node_data)[vi].heap.empty()) {
   1.289 -                _blossom_set->increase(v, std::numeric_limits<Value>::max());
   1.290 -              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
   1.291 -                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
   1.292 -              }
   1.293 -
   1.294 -              if ((*_blossom_data)[vb].status == MATCHED) {
   1.295 -                if (_blossom_set->classPrio(vb) ==
   1.296 -                    std::numeric_limits<Value>::max()) {
   1.297 -                  _delta2->erase(vb);
   1.298 -                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
   1.299 -                           (*_blossom_data)[vb].offset) {
   1.300 -                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
   1.301 -                                   (*_blossom_data)[vb].offset);
   1.302 -                }
   1.303 -              }
   1.304 -            }
   1.305 -          }
   1.306 -        }
   1.307 -      }
   1.308 -    }
   1.309 -
   1.310 -    void oddToMatched(int blossom) {
   1.311 -      (*_blossom_data)[blossom].offset -= _delta_sum;
   1.312 -
   1.313 -      if (_blossom_set->classPrio(blossom) !=
   1.314 -          std::numeric_limits<Value>::max()) {
   1.315 -        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
   1.316 -                       (*_blossom_data)[blossom].offset);
   1.317 -      }
   1.318 -
   1.319 -      if (!_blossom_set->trivial(blossom)) {
   1.320 -        _delta4->erase(blossom);
   1.321 -      }
   1.322 -    }
   1.323 -
   1.324 -    void oddToEven(int blossom, int tree) {
   1.325 -      if (!_blossom_set->trivial(blossom)) {
   1.326 -        _delta4->erase(blossom);
   1.327 -        (*_blossom_data)[blossom].pot -=
   1.328 -          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
   1.329 -      }
   1.330 -
   1.331 -      for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
   1.332 -           n != INVALID; ++n) {
   1.333 -        int ni = (*_node_index)[n];
   1.334 -
   1.335 -        _blossom_set->increase(n, std::numeric_limits<Value>::max());
   1.336 -
   1.337 -        (*_node_data)[ni].heap.clear();
   1.338 -        (*_node_data)[ni].heap_index.clear();
   1.339 -        (*_node_data)[ni].pot +=
   1.340 -          2 * _delta_sum - (*_blossom_data)[blossom].offset;
   1.341 -
   1.342 -        _delta1->push(n, (*_node_data)[ni].pot);
   1.343 -
   1.344 -        for (InArcIt e(_graph, n); e != INVALID; ++e) {
   1.345 -          Node v = _graph.source(e);
   1.346 -          int vb = _blossom_set->find(v);
   1.347 -          int vi = (*_node_index)[v];
   1.348 -
   1.349 -          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   1.350 -            dualScale * _weight[e];
   1.351 -
   1.352 -          if ((*_blossom_data)[vb].status == EVEN) {
   1.353 -            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
   1.354 -              _delta3->push(e, rw / 2);
   1.355 -            }
   1.356 -          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
   1.357 -            if (_delta3->state(e) != _delta3->IN_HEAP) {
   1.358 -              _delta3->push(e, rw);
   1.359 -            }
   1.360 -          } else {
   1.361 -
   1.362 -            typename std::map<int, Arc>::iterator it =
   1.363 -              (*_node_data)[vi].heap_index.find(tree);
   1.364 -
   1.365 -            if (it != (*_node_data)[vi].heap_index.end()) {
   1.366 -              if ((*_node_data)[vi].heap[it->second] > rw) {
   1.367 -                (*_node_data)[vi].heap.replace(it->second, e);
   1.368 -                (*_node_data)[vi].heap.decrease(e, rw);
   1.369 -                it->second = e;
   1.370 -              }
   1.371 -            } else {
   1.372 -              (*_node_data)[vi].heap.push(e, rw);
   1.373 -              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
   1.374 -            }
   1.375 -
   1.376 -            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
   1.377 -              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
   1.378 -
   1.379 -              if ((*_blossom_data)[vb].status == MATCHED) {
   1.380 -                if (_delta2->state(vb) != _delta2->IN_HEAP) {
   1.381 -                  _delta2->push(vb, _blossom_set->classPrio(vb) -
   1.382 -                               (*_blossom_data)[vb].offset);
   1.383 -                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
   1.384                             (*_blossom_data)[vb].offset) {
   1.385                    _delta2->decrease(vb, _blossom_set->classPrio(vb) -
   1.386                                     (*_blossom_data)[vb].offset);
   1.387 @@ -1157,43 +959,145 @@
   1.388        (*_blossom_data)[blossom].offset = 0;
   1.389      }
   1.390  
   1.391 -
   1.392 -    void matchedToUnmatched(int blossom) {
   1.393 +    void matchedToOdd(int blossom) {
   1.394        if (_delta2->state(blossom) == _delta2->IN_HEAP) {
   1.395          _delta2->erase(blossom);
   1.396        }
   1.397 +      (*_blossom_data)[blossom].offset += _delta_sum;
   1.398 +      if (!_blossom_set->trivial(blossom)) {
   1.399 +        _delta4->push(blossom, (*_blossom_data)[blossom].pot / 2 +
   1.400 +                      (*_blossom_data)[blossom].offset);
   1.401 +      }
   1.402 +    }
   1.403 +
   1.404 +    void evenToMatched(int blossom, int tree) {
   1.405 +      if (!_blossom_set->trivial(blossom)) {
   1.406 +        (*_blossom_data)[blossom].pot += 2 * _delta_sum;
   1.407 +      }
   1.408  
   1.409        for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
   1.410             n != INVALID; ++n) {
   1.411          int ni = (*_node_index)[n];
   1.412 -
   1.413 -        _blossom_set->increase(n, std::numeric_limits<Value>::max());
   1.414 -
   1.415 -        (*_node_data)[ni].heap.clear();
   1.416 -        (*_node_data)[ni].heap_index.clear();
   1.417 -
   1.418 -        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
   1.419 -          Node v = _graph.target(e);
   1.420 +        (*_node_data)[ni].pot -= _delta_sum;
   1.421 +
   1.422 +        _delta1->erase(n);
   1.423 +
   1.424 +        for (InArcIt e(_graph, n); e != INVALID; ++e) {
   1.425 +          Node v = _graph.source(e);
   1.426            int vb = _blossom_set->find(v);
   1.427            int vi = (*_node_index)[v];
   1.428  
   1.429            Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   1.430              dualScale * _weight[e];
   1.431  
   1.432 -          if ((*_blossom_data)[vb].status == EVEN) {
   1.433 -            if (_delta3->state(e) != _delta3->IN_HEAP) {
   1.434 -              _delta3->push(e, rw);
   1.435 +          if (vb == blossom) {
   1.436 +            if (_delta3->state(e) == _delta3->IN_HEAP) {
   1.437 +              _delta3->erase(e);
   1.438 +            }
   1.439 +          } else if ((*_blossom_data)[vb].status == EVEN) {
   1.440 +
   1.441 +            if (_delta3->state(e) == _delta3->IN_HEAP) {
   1.442 +              _delta3->erase(e);
   1.443 +            }
   1.444 +
   1.445 +            int vt = _tree_set->find(vb);
   1.446 +
   1.447 +            if (vt != tree) {
   1.448 +
   1.449 +              Arc r = _graph.oppositeArc(e);
   1.450 +
   1.451 +              typename std::map<int, Arc>::iterator it =
   1.452 +                (*_node_data)[ni].heap_index.find(vt);
   1.453 +
   1.454 +              if (it != (*_node_data)[ni].heap_index.end()) {
   1.455 +                if ((*_node_data)[ni].heap[it->second] > rw) {
   1.456 +                  (*_node_data)[ni].heap.replace(it->second, r);
   1.457 +                  (*_node_data)[ni].heap.decrease(r, rw);
   1.458 +                  it->second = r;
   1.459 +                }
   1.460 +              } else {
   1.461 +                (*_node_data)[ni].heap.push(r, rw);
   1.462 +                (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
   1.463 +              }
   1.464 +
   1.465 +              if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
   1.466 +                _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
   1.467 +
   1.468 +                if (_delta2->state(blossom) != _delta2->IN_HEAP) {
   1.469 +                  _delta2->push(blossom, _blossom_set->classPrio(blossom) -
   1.470 +                               (*_blossom_data)[blossom].offset);
   1.471 +                } else if ((*_delta2)[blossom] >
   1.472 +                           _blossom_set->classPrio(blossom) -
   1.473 +                           (*_blossom_data)[blossom].offset){
   1.474 +                  _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
   1.475 +                                   (*_blossom_data)[blossom].offset);
   1.476 +                }
   1.477 +              }
   1.478 +            }
   1.479 +          } else {
   1.480 +
   1.481 +            typename std::map<int, Arc>::iterator it =
   1.482 +              (*_node_data)[vi].heap_index.find(tree);
   1.483 +
   1.484 +            if (it != (*_node_data)[vi].heap_index.end()) {
   1.485 +              (*_node_data)[vi].heap.erase(it->second);
   1.486 +              (*_node_data)[vi].heap_index.erase(it);
   1.487 +              if ((*_node_data)[vi].heap.empty()) {
   1.488 +                _blossom_set->increase(v, std::numeric_limits<Value>::max());
   1.489 +              } else if ((*_blossom_set)[v] < (*_node_data)[vi].heap.prio()) {
   1.490 +                _blossom_set->increase(v, (*_node_data)[vi].heap.prio());
   1.491 +              }
   1.492 +
   1.493 +              if ((*_blossom_data)[vb].status == MATCHED) {
   1.494 +                if (_blossom_set->classPrio(vb) ==
   1.495 +                    std::numeric_limits<Value>::max()) {
   1.496 +                  _delta2->erase(vb);
   1.497 +                } else if ((*_delta2)[vb] < _blossom_set->classPrio(vb) -
   1.498 +                           (*_blossom_data)[vb].offset) {
   1.499 +                  _delta2->increase(vb, _blossom_set->classPrio(vb) -
   1.500 +                                   (*_blossom_data)[vb].offset);
   1.501 +                }
   1.502 +              }
   1.503              }
   1.504            }
   1.505          }
   1.506        }
   1.507      }
   1.508  
   1.509 -    void unmatchedToMatched(int blossom) {
   1.510 +    void oddToMatched(int blossom) {
   1.511 +      (*_blossom_data)[blossom].offset -= _delta_sum;
   1.512 +
   1.513 +      if (_blossom_set->classPrio(blossom) !=
   1.514 +          std::numeric_limits<Value>::max()) {
   1.515 +        _delta2->push(blossom, _blossom_set->classPrio(blossom) -
   1.516 +                      (*_blossom_data)[blossom].offset);
   1.517 +      }
   1.518 +
   1.519 +      if (!_blossom_set->trivial(blossom)) {
   1.520 +        _delta4->erase(blossom);
   1.521 +      }
   1.522 +    }
   1.523 +
   1.524 +    void oddToEven(int blossom, int tree) {
   1.525 +      if (!_blossom_set->trivial(blossom)) {
   1.526 +        _delta4->erase(blossom);
   1.527 +        (*_blossom_data)[blossom].pot -=
   1.528 +          2 * (2 * _delta_sum - (*_blossom_data)[blossom].offset);
   1.529 +      }
   1.530 +
   1.531        for (typename BlossomSet::ItemIt n(*_blossom_set, blossom);
   1.532             n != INVALID; ++n) {
   1.533          int ni = (*_node_index)[n];
   1.534  
   1.535 +        _blossom_set->increase(n, std::numeric_limits<Value>::max());
   1.536 +
   1.537 +        (*_node_data)[ni].heap.clear();
   1.538 +        (*_node_data)[ni].heap_index.clear();
   1.539 +        (*_node_data)[ni].pot +=
   1.540 +          2 * _delta_sum - (*_blossom_data)[blossom].offset;
   1.541 +
   1.542 +        _delta1->push(n, (*_node_data)[ni].pot);
   1.543 +
   1.544          for (InArcIt e(_graph, n); e != INVALID; ++e) {
   1.545            Node v = _graph.source(e);
   1.546            int vb = _blossom_set->find(v);
   1.547 @@ -1202,54 +1106,44 @@
   1.548            Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   1.549              dualScale * _weight[e];
   1.550  
   1.551 -          if (vb == blossom) {
   1.552 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   1.553 -              _delta3->erase(e);
   1.554 +          if ((*_blossom_data)[vb].status == EVEN) {
   1.555 +            if (_delta3->state(e) != _delta3->IN_HEAP && blossom != vb) {
   1.556 +              _delta3->push(e, rw / 2);
   1.557              }
   1.558 -          } else if ((*_blossom_data)[vb].status == EVEN) {
   1.559 -
   1.560 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   1.561 -              _delta3->erase(e);
   1.562 -            }
   1.563 -
   1.564 -            int vt = _tree_set->find(vb);
   1.565 -
   1.566 -            Arc r = _graph.oppositeArc(e);
   1.567 +          } else {
   1.568  
   1.569              typename std::map<int, Arc>::iterator it =
   1.570 -              (*_node_data)[ni].heap_index.find(vt);
   1.571 -
   1.572 -            if (it != (*_node_data)[ni].heap_index.end()) {
   1.573 -              if ((*_node_data)[ni].heap[it->second] > rw) {
   1.574 -                (*_node_data)[ni].heap.replace(it->second, r);
   1.575 -                (*_node_data)[ni].heap.decrease(r, rw);
   1.576 -                it->second = r;
   1.577 +              (*_node_data)[vi].heap_index.find(tree);
   1.578 +
   1.579 +            if (it != (*_node_data)[vi].heap_index.end()) {
   1.580 +              if ((*_node_data)[vi].heap[it->second] > rw) {
   1.581 +                (*_node_data)[vi].heap.replace(it->second, e);
   1.582 +                (*_node_data)[vi].heap.decrease(e, rw);
   1.583 +                it->second = e;
   1.584                }
   1.585              } else {
   1.586 -              (*_node_data)[ni].heap.push(r, rw);
   1.587 -              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, r));
   1.588 +              (*_node_data)[vi].heap.push(e, rw);
   1.589 +              (*_node_data)[vi].heap_index.insert(std::make_pair(tree, e));
   1.590              }
   1.591  
   1.592 -            if ((*_blossom_set)[n] > (*_node_data)[ni].heap.prio()) {
   1.593 -              _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
   1.594 -
   1.595 -              if (_delta2->state(blossom) != _delta2->IN_HEAP) {
   1.596 -                _delta2->push(blossom, _blossom_set->classPrio(blossom) -
   1.597 -                             (*_blossom_data)[blossom].offset);
   1.598 -              } else if ((*_delta2)[blossom] > _blossom_set->classPrio(blossom)-
   1.599 -                         (*_blossom_data)[blossom].offset){
   1.600 -                _delta2->decrease(blossom, _blossom_set->classPrio(blossom) -
   1.601 -                                 (*_blossom_data)[blossom].offset);
   1.602 +            if ((*_blossom_set)[v] > (*_node_data)[vi].heap.prio()) {
   1.603 +              _blossom_set->decrease(v, (*_node_data)[vi].heap.prio());
   1.604 +
   1.605 +              if ((*_blossom_data)[vb].status == MATCHED) {
   1.606 +                if (_delta2->state(vb) != _delta2->IN_HEAP) {
   1.607 +                  _delta2->push(vb, _blossom_set->classPrio(vb) -
   1.608 +                               (*_blossom_data)[vb].offset);
   1.609 +                } else if ((*_delta2)[vb] > _blossom_set->classPrio(vb) -
   1.610 +                           (*_blossom_data)[vb].offset) {
   1.611 +                  _delta2->decrease(vb, _blossom_set->classPrio(vb) -
   1.612 +                                   (*_blossom_data)[vb].offset);
   1.613 +                }
   1.614                }
   1.615              }
   1.616 -
   1.617 -          } else if ((*_blossom_data)[vb].status == UNMATCHED) {
   1.618 -            if (_delta3->state(e) == _delta3->IN_HEAP) {
   1.619 -              _delta3->erase(e);
   1.620 -            }
   1.621            }
   1.622          }
   1.623        }
   1.624 +      (*_blossom_data)[blossom].offset = 0;
   1.625      }
   1.626  
   1.627      void alternatePath(int even, int tree) {
   1.628 @@ -1294,39 +1188,42 @@
   1.629        alternatePath(blossom, tree);
   1.630        destroyTree(tree);
   1.631  
   1.632 -      (*_blossom_data)[blossom].status = UNMATCHED;
   1.633        (*_blossom_data)[blossom].base = node;
   1.634 -      matchedToUnmatched(blossom);
   1.635 +      (*_blossom_data)[blossom].next = INVALID;
   1.636      }
   1.637  
   1.638 -
   1.639      void augmentOnEdge(const Edge& edge) {
   1.640  
   1.641        int left = _blossom_set->find(_graph.u(edge));
   1.642        int right = _blossom_set->find(_graph.v(edge));
   1.643  
   1.644 -      if ((*_blossom_data)[left].status == EVEN) {
   1.645 -        int left_tree = _tree_set->find(left);
   1.646 -        alternatePath(left, left_tree);
   1.647 -        destroyTree(left_tree);
   1.648 -      } else {
   1.649 -        (*_blossom_data)[left].status = MATCHED;
   1.650 -        unmatchedToMatched(left);
   1.651 -      }
   1.652 -
   1.653 -      if ((*_blossom_data)[right].status == EVEN) {
   1.654 -        int right_tree = _tree_set->find(right);
   1.655 -        alternatePath(right, right_tree);
   1.656 -        destroyTree(right_tree);
   1.657 -      } else {
   1.658 -        (*_blossom_data)[right].status = MATCHED;
   1.659 -        unmatchedToMatched(right);
   1.660 -      }
   1.661 +      int left_tree = _tree_set->find(left);
   1.662 +      alternatePath(left, left_tree);
   1.663 +      destroyTree(left_tree);
   1.664 +
   1.665 +      int right_tree = _tree_set->find(right);
   1.666 +      alternatePath(right, right_tree);
   1.667 +      destroyTree(right_tree);
   1.668  
   1.669        (*_blossom_data)[left].next = _graph.direct(edge, true);
   1.670        (*_blossom_data)[right].next = _graph.direct(edge, false);
   1.671      }
   1.672  
   1.673 +    void augmentOnArc(const Arc& arc) {
   1.674 +
   1.675 +      int left = _blossom_set->find(_graph.source(arc));
   1.676 +      int right = _blossom_set->find(_graph.target(arc));
   1.677 +
   1.678 +      (*_blossom_data)[left].status = MATCHED;
   1.679 +
   1.680 +      int right_tree = _tree_set->find(right);
   1.681 +      alternatePath(right, right_tree);
   1.682 +      destroyTree(right_tree);
   1.683 +
   1.684 +      (*_blossom_data)[left].next = arc;
   1.685 +      (*_blossom_data)[right].next = _graph.oppositeArc(arc);
   1.686 +    }
   1.687 +
   1.688      void extendOnArc(const Arc& arc) {
   1.689        int base = _blossom_set->find(_graph.target(arc));
   1.690        int tree = _tree_set->find(base);
   1.691 @@ -1529,7 +1426,7 @@
   1.692            _tree_set->insert(sb, tree);
   1.693            (*_blossom_data)[sb].pred = pred;
   1.694            (*_blossom_data)[sb].next =
   1.695 -                           _graph.oppositeArc((*_blossom_data)[tb].next);
   1.696 +            _graph.oppositeArc((*_blossom_data)[tb].next);
   1.697  
   1.698            pred = (*_blossom_data)[ub].next;
   1.699  
   1.700 @@ -1629,7 +1526,7 @@
   1.701        }
   1.702  
   1.703        for (int i = 0; i < int(blossoms.size()); ++i) {
   1.704 -        if ((*_blossom_data)[blossoms[i]].status == MATCHED) {
   1.705 +        if ((*_blossom_data)[blossoms[i]].next != INVALID) {
   1.706  
   1.707            Value offset = (*_blossom_data)[blossoms[i]].offset;
   1.708            (*_blossom_data)[blossoms[i]].pot += 2 * offset;
   1.709 @@ -1667,10 +1564,16 @@
   1.710          _delta3_index(0), _delta3(0),
   1.711          _delta4_index(0), _delta4(0),
   1.712  
   1.713 -        _delta_sum() {}
   1.714 +        _delta_sum(), _unmatched(0),
   1.715 +
   1.716 +        _fractional(0)
   1.717 +    {}
   1.718  
   1.719      ~MaxWeightedMatching() {
   1.720        destroyStructures();
   1.721 +      if (_fractional) {
   1.722 +        delete _fractional;
   1.723 +      }
   1.724      }
   1.725  
   1.726      /// \name Execution Control
   1.727 @@ -1699,6 +1602,8 @@
   1.728          (*_delta4_index)[i] = _delta4->PRE_HEAP;
   1.729        }
   1.730  
   1.731 +      _unmatched = _node_num;
   1.732 +
   1.733        int index = 0;
   1.734        for (NodeIt n(_graph); n != INVALID; ++n) {
   1.735          Value max = 0;
   1.736 @@ -1733,18 +1638,155 @@
   1.737        }
   1.738      }
   1.739  
   1.740 +    /// \brief Initialize the algorithm with fractional matching
   1.741 +    ///
   1.742 +    /// This function initializes the algorithm with a fractional
   1.743 +    /// matching. This initialization is also called jumpstart heuristic.
   1.744 +    void fractionalInit() {
   1.745 +      createStructures();
   1.746 +      
   1.747 +      if (_fractional == 0) {
   1.748 +        _fractional = new FractionalMatching(_graph, _weight, false);
   1.749 +      }
   1.750 +      _fractional->run();
   1.751 +
   1.752 +      for (ArcIt e(_graph); e != INVALID; ++e) {
   1.753 +        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
   1.754 +      }
   1.755 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.756 +        (*_delta1_index)[n] = _delta1->PRE_HEAP;
   1.757 +      }
   1.758 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
   1.759 +        (*_delta3_index)[e] = _delta3->PRE_HEAP;
   1.760 +      }
   1.761 +      for (int i = 0; i < _blossom_num; ++i) {
   1.762 +        (*_delta2_index)[i] = _delta2->PRE_HEAP;
   1.763 +        (*_delta4_index)[i] = _delta4->PRE_HEAP;
   1.764 +      }
   1.765 +
   1.766 +      _unmatched = 0;
   1.767 +
   1.768 +      int index = 0;
   1.769 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.770 +        Value pot = _fractional->nodeValue(n);
   1.771 +        (*_node_index)[n] = index;
   1.772 +        (*_node_data)[index].pot = pot;
   1.773 +        int blossom =
   1.774 +          _blossom_set->insert(n, std::numeric_limits<Value>::max());
   1.775 +
   1.776 +        (*_blossom_data)[blossom].status = MATCHED;
   1.777 +        (*_blossom_data)[blossom].pred = INVALID;
   1.778 +        (*_blossom_data)[blossom].next = _fractional->matching(n);
   1.779 +        if (_fractional->matching(n) == INVALID) {
   1.780 +          (*_blossom_data)[blossom].base = n;
   1.781 +        }
   1.782 +        (*_blossom_data)[blossom].pot = 0;
   1.783 +        (*_blossom_data)[blossom].offset = 0;
   1.784 +        ++index;
   1.785 +      }
   1.786 +
   1.787 +      typename Graph::template NodeMap<bool> processed(_graph, false);
   1.788 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.789 +        if (processed[n]) continue;
   1.790 +        processed[n] = true;
   1.791 +        if (_fractional->matching(n) == INVALID) continue;
   1.792 +        int num = 1;
   1.793 +        Node v = _graph.target(_fractional->matching(n));
   1.794 +        while (n != v) {
   1.795 +          processed[v] = true;
   1.796 +          v = _graph.target(_fractional->matching(v));
   1.797 +          ++num;
   1.798 +        }
   1.799 +
   1.800 +        if (num % 2 == 1) {
   1.801 +          std::vector<int> subblossoms(num);
   1.802 +
   1.803 +          subblossoms[--num] = _blossom_set->find(n);
   1.804 +          _delta1->push(n, _fractional->nodeValue(n));
   1.805 +          v = _graph.target(_fractional->matching(n));
   1.806 +          while (n != v) {
   1.807 +            subblossoms[--num] = _blossom_set->find(v);
   1.808 +            _delta1->push(v, _fractional->nodeValue(v));
   1.809 +            v = _graph.target(_fractional->matching(v));            
   1.810 +          }
   1.811 +          
   1.812 +          int surface = 
   1.813 +            _blossom_set->join(subblossoms.begin(), subblossoms.end());
   1.814 +          (*_blossom_data)[surface].status = EVEN;
   1.815 +          (*_blossom_data)[surface].pred = INVALID;
   1.816 +          (*_blossom_data)[surface].next = INVALID;
   1.817 +          (*_blossom_data)[surface].pot = 0;
   1.818 +          (*_blossom_data)[surface].offset = 0;
   1.819 +          
   1.820 +          _tree_set->insert(surface);
   1.821 +          ++_unmatched;
   1.822 +        }
   1.823 +      }
   1.824 +
   1.825 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
   1.826 +        int si = (*_node_index)[_graph.u(e)];
   1.827 +        int sb = _blossom_set->find(_graph.u(e));
   1.828 +        int ti = (*_node_index)[_graph.v(e)];
   1.829 +        int tb = _blossom_set->find(_graph.v(e));
   1.830 +        if ((*_blossom_data)[sb].status == EVEN &&
   1.831 +            (*_blossom_data)[tb].status == EVEN && sb != tb) {
   1.832 +          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
   1.833 +                            dualScale * _weight[e]) / 2);
   1.834 +        }
   1.835 +      }
   1.836 +
   1.837 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.838 +        int nb = _blossom_set->find(n);
   1.839 +        if ((*_blossom_data)[nb].status != MATCHED) continue;
   1.840 +        int ni = (*_node_index)[n];
   1.841 +
   1.842 +        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
   1.843 +          Node v = _graph.target(e);
   1.844 +          int vb = _blossom_set->find(v);
   1.845 +          int vi = (*_node_index)[v];
   1.846 +
   1.847 +          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
   1.848 +            dualScale * _weight[e];
   1.849 +
   1.850 +          if ((*_blossom_data)[vb].status == EVEN) {
   1.851 +
   1.852 +            int vt = _tree_set->find(vb);
   1.853 +
   1.854 +            typename std::map<int, Arc>::iterator it =
   1.855 +              (*_node_data)[ni].heap_index.find(vt);
   1.856 +
   1.857 +            if (it != (*_node_data)[ni].heap_index.end()) {
   1.858 +              if ((*_node_data)[ni].heap[it->second] > rw) {
   1.859 +                (*_node_data)[ni].heap.replace(it->second, e);
   1.860 +                (*_node_data)[ni].heap.decrease(e, rw);
   1.861 +                it->second = e;
   1.862 +              }
   1.863 +            } else {
   1.864 +              (*_node_data)[ni].heap.push(e, rw);
   1.865 +              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
   1.866 +            }
   1.867 +          }
   1.868 +        }
   1.869 +            
   1.870 +        if (!(*_node_data)[ni].heap.empty()) {
   1.871 +          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
   1.872 +          _delta2->push(nb, _blossom_set->classPrio(nb));
   1.873 +        }
   1.874 +      }
   1.875 +    }
   1.876 +
   1.877      /// \brief Start the algorithm
   1.878      ///
   1.879      /// This function starts the algorithm.
   1.880      ///
   1.881 -    /// \pre \ref init() must be called before using this function.
   1.882 +    /// \pre \ref init() or \ref fractionalInit() must be called
   1.883 +    /// before using this function.
   1.884      void start() {
   1.885        enum OpType {
   1.886          D1, D2, D3, D4
   1.887        };
   1.888  
   1.889 -      int unmatched = _node_num;
   1.890 -      while (unmatched > 0) {
   1.891 +      while (_unmatched > 0) {
   1.892          Value d1 = !_delta1->empty() ?
   1.893            _delta1->prio() : std::numeric_limits<Value>::max();
   1.894  
   1.895 @@ -1757,26 +1799,30 @@
   1.896          Value d4 = !_delta4->empty() ?
   1.897            _delta4->prio() : std::numeric_limits<Value>::max();
   1.898  
   1.899 -        _delta_sum = d1; OpType ot = D1;
   1.900 +        _delta_sum = d3; OpType ot = D3;
   1.901 +        if (d1 < _delta_sum) { _delta_sum = d1; ot = D1; }
   1.902          if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
   1.903 -        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
   1.904          if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
   1.905  
   1.906 -
   1.907          switch (ot) {
   1.908          case D1:
   1.909            {
   1.910              Node n = _delta1->top();
   1.911              unmatchNode(n);
   1.912 -            --unmatched;
   1.913 +            --_unmatched;
   1.914            }
   1.915            break;
   1.916          case D2:
   1.917            {
   1.918              int blossom = _delta2->top();
   1.919              Node n = _blossom_set->classTop(blossom);
   1.920 -            Arc e = (*_node_data)[(*_node_index)[n]].heap.top();
   1.921 -            extendOnArc(e);
   1.922 +            Arc a = (*_node_data)[(*_node_index)[n]].heap.top();
   1.923 +            if ((*_blossom_data)[blossom].next == INVALID) {
   1.924 +              augmentOnArc(a);
   1.925 +              --_unmatched;
   1.926 +            } else {
   1.927 +              extendOnArc(a);
   1.928 +            }
   1.929            }
   1.930            break;
   1.931          case D3:
   1.932 @@ -1789,26 +1835,14 @@
   1.933              if (left_blossom == right_blossom) {
   1.934                _delta3->pop();
   1.935              } else {
   1.936 -              int left_tree;
   1.937 -              if ((*_blossom_data)[left_blossom].status == EVEN) {
   1.938 -                left_tree = _tree_set->find(left_blossom);
   1.939 -              } else {
   1.940 -                left_tree = -1;
   1.941 -                ++unmatched;
   1.942 -              }
   1.943 -              int right_tree;
   1.944 -              if ((*_blossom_data)[right_blossom].status == EVEN) {
   1.945 -                right_tree = _tree_set->find(right_blossom);
   1.946 -              } else {
   1.947 -                right_tree = -1;
   1.948 -                ++unmatched;
   1.949 -              }
   1.950 +              int left_tree = _tree_set->find(left_blossom);
   1.951 +              int right_tree = _tree_set->find(right_blossom);
   1.952  
   1.953                if (left_tree == right_tree) {
   1.954                  shrinkOnEdge(e, left_tree);
   1.955                } else {
   1.956                  augmentOnEdge(e);
   1.957 -                unmatched -= 2;
   1.958 +                _unmatched -= 2;
   1.959                }
   1.960              }
   1.961            } break;
   1.962 @@ -1826,18 +1860,18 @@
   1.963      ///
   1.964      /// \note mwm.run() is just a shortcut of the following code.
   1.965      /// \code
   1.966 -    ///   mwm.init();
   1.967 +    ///   mwm.fractionalInit();
   1.968      ///   mwm.start();
   1.969      /// \endcode
   1.970      void run() {
   1.971 -      init();
   1.972 +      fractionalInit();
   1.973        start();
   1.974      }
   1.975  
   1.976      /// @}
   1.977  
   1.978      /// \name Primal Solution
   1.979 -    /// Functions to get the primal solution, i.e. the maximum weighted 
   1.980 +    /// Functions to get the primal solution, i.e. the maximum weighted
   1.981      /// matching.\n
   1.982      /// Either \ref run() or \ref start() function should be called before
   1.983      /// using them.
   1.984 @@ -1856,7 +1890,7 @@
   1.985            sum += _weight[(*_matching)[n]];
   1.986          }
   1.987        }
   1.988 -      return sum /= 2;
   1.989 +      return sum / 2;
   1.990      }
   1.991  
   1.992      /// \brief Return the size (cardinality) of the matching.
   1.993 @@ -1876,7 +1910,7 @@
   1.994  
   1.995      /// \brief Return \c true if the given edge is in the matching.
   1.996      ///
   1.997 -    /// This function returns \c true if the given edge is in the found 
   1.998 +    /// This function returns \c true if the given edge is in the found
   1.999      /// matching.
  1.1000      ///
  1.1001      /// \pre Either run() or start() must be called before using this function.
  1.1002 @@ -1887,7 +1921,7 @@
  1.1003      /// \brief Return the matching arc (or edge) incident to the given node.
  1.1004      ///
  1.1005      /// This function returns the matching arc (or edge) incident to the
  1.1006 -    /// given node in the found matching or \c INVALID if the node is 
  1.1007 +    /// given node in the found matching or \c INVALID if the node is
  1.1008      /// not covered by the matching.
  1.1009      ///
  1.1010      /// \pre Either run() or start() must be called before using this function.
  1.1011 @@ -1905,7 +1939,7 @@
  1.1012  
  1.1013      /// \brief Return the mate of the given node.
  1.1014      ///
  1.1015 -    /// This function returns the mate of the given node in the found 
  1.1016 +    /// This function returns the mate of the given node in the found
  1.1017      /// matching or \c INVALID if the node is not covered by the matching.
  1.1018      ///
  1.1019      /// \pre Either run() or start() must be called before using this function.
  1.1020 @@ -1925,8 +1959,8 @@
  1.1021  
  1.1022      /// \brief Return the value of the dual solution.
  1.1023      ///
  1.1024 -    /// This function returns the value of the dual solution. 
  1.1025 -    /// It should be equal to the primal value scaled by \ref dualScale 
  1.1026 +    /// This function returns the value of the dual solution.
  1.1027 +    /// It should be equal to the primal value scaled by \ref dualScale
  1.1028      /// "dual scale".
  1.1029      ///
  1.1030      /// \pre Either run() or start() must be called before using this function.
  1.1031 @@ -1981,9 +2015,9 @@
  1.1032  
  1.1033      /// \brief Iterator for obtaining the nodes of a blossom.
  1.1034      ///
  1.1035 -    /// This class provides an iterator for obtaining the nodes of the 
  1.1036 +    /// This class provides an iterator for obtaining the nodes of the
  1.1037      /// given blossom. It lists a subset of the nodes.
  1.1038 -    /// Before using this iterator, you must allocate a 
  1.1039 +    /// Before using this iterator, you must allocate a
  1.1040      /// MaxWeightedMatching class and execute it.
  1.1041      class BlossomIt {
  1.1042      public:
  1.1043 @@ -1992,8 +2026,8 @@
  1.1044        ///
  1.1045        /// Constructor to get the nodes of the given variable.
  1.1046        ///
  1.1047 -      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or 
  1.1048 -      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be 
  1.1049 +      /// \pre Either \ref MaxWeightedMatching::run() "algorithm.run()" or
  1.1050 +      /// \ref MaxWeightedMatching::start() "algorithm.start()" must be
  1.1051        /// called before initializing this iterator.
  1.1052        BlossomIt(const MaxWeightedMatching& algorithm, int variable)
  1.1053          : _algorithm(&algorithm)
  1.1054 @@ -2046,8 +2080,8 @@
  1.1055    /// is based on extensive use of priority queues and provides
  1.1056    /// \f$O(nm\log n)\f$ time complexity.
  1.1057    ///
  1.1058 -  /// The maximum weighted perfect matching problem is to find a subset of 
  1.1059 -  /// the edges in an undirected graph with maximum overall weight for which 
  1.1060 +  /// The maximum weighted perfect matching problem is to find a subset of
  1.1061 +  /// the edges in an undirected graph with maximum overall weight for which
  1.1062    /// each node has exactly one incident edge.
  1.1063    /// It can be formulated with the following linear program.
  1.1064    /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
  1.1065 @@ -2070,16 +2104,16 @@
  1.1066    /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
  1.1067        \frac{\vert B \vert - 1}{2}z_B\f] */
  1.1068    ///
  1.1069 -  /// The algorithm can be executed with the run() function. 
  1.1070 +  /// The algorithm can be executed with the run() function.
  1.1071    /// After it the matching (the primal solution) and the dual solution
  1.1072 -  /// can be obtained using the query functions and the 
  1.1073 -  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class, 
  1.1074 -  /// which is able to iterate on the nodes of a blossom. 
  1.1075 +  /// can be obtained using the query functions and the
  1.1076 +  /// \ref MaxWeightedPerfectMatching::BlossomIt "BlossomIt" nested class,
  1.1077 +  /// which is able to iterate on the nodes of a blossom.
  1.1078    /// If the value type is integer, then the dual solution is multiplied
  1.1079    /// by \ref MaxWeightedMatching::dualScale "4".
  1.1080    ///
  1.1081    /// \tparam GR The undirected graph type the algorithm runs on.
  1.1082 -  /// \tparam WM The type edge weight map. The default type is 
  1.1083 +  /// \tparam WM The type edge weight map. The default type is
  1.1084    /// \ref concepts::Graph::EdgeMap "GR::EdgeMap<int>".
  1.1085  #ifdef DOXYGEN
  1.1086    template <typename GR, typename WM>
  1.1087 @@ -2190,6 +2224,11 @@
  1.1088      BinHeap<Value, IntIntMap> *_delta4;
  1.1089  
  1.1090      Value _delta_sum;
  1.1091 +    int _unmatched;
  1.1092 +
  1.1093 +    typedef MaxWeightedPerfectFractionalMatching<Graph, WeightMap> 
  1.1094 +    FractionalMatching;
  1.1095 +    FractionalMatching *_fractional;
  1.1096  
  1.1097      void createStructures() {
  1.1098        _node_num = countNodes(_graph);
  1.1099 @@ -2233,9 +2272,6 @@
  1.1100      }
  1.1101  
  1.1102      void destroyStructures() {
  1.1103 -      _node_num = countNodes(_graph);
  1.1104 -      _blossom_num = _node_num * 3 / 2;
  1.1105 -
  1.1106        if (_matching) {
  1.1107          delete _matching;
  1.1108        }
  1.1109 @@ -2908,10 +2944,16 @@
  1.1110          _delta3_index(0), _delta3(0),
  1.1111          _delta4_index(0), _delta4(0),
  1.1112  
  1.1113 -        _delta_sum() {}
  1.1114 +        _delta_sum(), _unmatched(0),
  1.1115 +
  1.1116 +        _fractional(0)
  1.1117 +    {}
  1.1118  
  1.1119      ~MaxWeightedPerfectMatching() {
  1.1120        destroyStructures();
  1.1121 +      if (_fractional) {
  1.1122 +        delete _fractional;
  1.1123 +      }
  1.1124      }
  1.1125  
  1.1126      /// \name Execution Control
  1.1127 @@ -2937,6 +2979,8 @@
  1.1128          (*_delta4_index)[i] = _delta4->PRE_HEAP;
  1.1129        }
  1.1130  
  1.1131 +      _unmatched = _node_num;
  1.1132 +
  1.1133        int index = 0;
  1.1134        for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1135          Value max = - std::numeric_limits<Value>::max();
  1.1136 @@ -2970,18 +3014,152 @@
  1.1137        }
  1.1138      }
  1.1139  
  1.1140 +    /// \brief Initialize the algorithm with fractional matching
  1.1141 +    ///
  1.1142 +    /// This function initializes the algorithm with a fractional
  1.1143 +    /// matching. This initialization is also called jumpstart heuristic.
  1.1144 +    void fractionalInit() {
  1.1145 +      createStructures();
  1.1146 +      
  1.1147 +      if (_fractional == 0) {
  1.1148 +        _fractional = new FractionalMatching(_graph, _weight, false);
  1.1149 +      }
  1.1150 +      if (!_fractional->run()) {
  1.1151 +        _unmatched = -1;
  1.1152 +        return;
  1.1153 +      }
  1.1154 +
  1.1155 +      for (ArcIt e(_graph); e != INVALID; ++e) {
  1.1156 +        (*_node_heap_index)[e] = BinHeap<Value, IntArcMap>::PRE_HEAP;
  1.1157 +      }
  1.1158 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
  1.1159 +        (*_delta3_index)[e] = _delta3->PRE_HEAP;
  1.1160 +      }
  1.1161 +      for (int i = 0; i < _blossom_num; ++i) {
  1.1162 +        (*_delta2_index)[i] = _delta2->PRE_HEAP;
  1.1163 +        (*_delta4_index)[i] = _delta4->PRE_HEAP;
  1.1164 +      }
  1.1165 +
  1.1166 +      _unmatched = 0;
  1.1167 +
  1.1168 +      int index = 0;
  1.1169 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1170 +        Value pot = _fractional->nodeValue(n);
  1.1171 +        (*_node_index)[n] = index;
  1.1172 +        (*_node_data)[index].pot = pot;
  1.1173 +        int blossom =
  1.1174 +          _blossom_set->insert(n, std::numeric_limits<Value>::max());
  1.1175 +
  1.1176 +        (*_blossom_data)[blossom].status = MATCHED;
  1.1177 +        (*_blossom_data)[blossom].pred = INVALID;
  1.1178 +        (*_blossom_data)[blossom].next = _fractional->matching(n);
  1.1179 +        (*_blossom_data)[blossom].pot = 0;
  1.1180 +        (*_blossom_data)[blossom].offset = 0;
  1.1181 +        ++index;
  1.1182 +      }
  1.1183 +
  1.1184 +      typename Graph::template NodeMap<bool> processed(_graph, false);
  1.1185 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1186 +        if (processed[n]) continue;
  1.1187 +        processed[n] = true;
  1.1188 +        if (_fractional->matching(n) == INVALID) continue;
  1.1189 +        int num = 1;
  1.1190 +        Node v = _graph.target(_fractional->matching(n));
  1.1191 +        while (n != v) {
  1.1192 +          processed[v] = true;
  1.1193 +          v = _graph.target(_fractional->matching(v));
  1.1194 +          ++num;
  1.1195 +        }
  1.1196 +
  1.1197 +        if (num % 2 == 1) {
  1.1198 +          std::vector<int> subblossoms(num);
  1.1199 +
  1.1200 +          subblossoms[--num] = _blossom_set->find(n);
  1.1201 +          v = _graph.target(_fractional->matching(n));
  1.1202 +          while (n != v) {
  1.1203 +            subblossoms[--num] = _blossom_set->find(v);
  1.1204 +            v = _graph.target(_fractional->matching(v));            
  1.1205 +          }
  1.1206 +          
  1.1207 +          int surface = 
  1.1208 +            _blossom_set->join(subblossoms.begin(), subblossoms.end());
  1.1209 +          (*_blossom_data)[surface].status = EVEN;
  1.1210 +          (*_blossom_data)[surface].pred = INVALID;
  1.1211 +          (*_blossom_data)[surface].next = INVALID;
  1.1212 +          (*_blossom_data)[surface].pot = 0;
  1.1213 +          (*_blossom_data)[surface].offset = 0;
  1.1214 +          
  1.1215 +          _tree_set->insert(surface);
  1.1216 +          ++_unmatched;
  1.1217 +        }
  1.1218 +      }
  1.1219 +
  1.1220 +      for (EdgeIt e(_graph); e != INVALID; ++e) {
  1.1221 +        int si = (*_node_index)[_graph.u(e)];
  1.1222 +        int sb = _blossom_set->find(_graph.u(e));
  1.1223 +        int ti = (*_node_index)[_graph.v(e)];
  1.1224 +        int tb = _blossom_set->find(_graph.v(e));
  1.1225 +        if ((*_blossom_data)[sb].status == EVEN &&
  1.1226 +            (*_blossom_data)[tb].status == EVEN && sb != tb) {
  1.1227 +          _delta3->push(e, ((*_node_data)[si].pot + (*_node_data)[ti].pot -
  1.1228 +                            dualScale * _weight[e]) / 2);
  1.1229 +        }
  1.1230 +      }
  1.1231 +
  1.1232 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1233 +        int nb = _blossom_set->find(n);
  1.1234 +        if ((*_blossom_data)[nb].status != MATCHED) continue;
  1.1235 +        int ni = (*_node_index)[n];
  1.1236 +
  1.1237 +        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
  1.1238 +          Node v = _graph.target(e);
  1.1239 +          int vb = _blossom_set->find(v);
  1.1240 +          int vi = (*_node_index)[v];
  1.1241 +
  1.1242 +          Value rw = (*_node_data)[ni].pot + (*_node_data)[vi].pot -
  1.1243 +            dualScale * _weight[e];
  1.1244 +
  1.1245 +          if ((*_blossom_data)[vb].status == EVEN) {
  1.1246 +
  1.1247 +            int vt = _tree_set->find(vb);
  1.1248 +
  1.1249 +            typename std::map<int, Arc>::iterator it =
  1.1250 +              (*_node_data)[ni].heap_index.find(vt);
  1.1251 +
  1.1252 +            if (it != (*_node_data)[ni].heap_index.end()) {
  1.1253 +              if ((*_node_data)[ni].heap[it->second] > rw) {
  1.1254 +                (*_node_data)[ni].heap.replace(it->second, e);
  1.1255 +                (*_node_data)[ni].heap.decrease(e, rw);
  1.1256 +                it->second = e;
  1.1257 +              }
  1.1258 +            } else {
  1.1259 +              (*_node_data)[ni].heap.push(e, rw);
  1.1260 +              (*_node_data)[ni].heap_index.insert(std::make_pair(vt, e));
  1.1261 +            }
  1.1262 +          }
  1.1263 +        }
  1.1264 +            
  1.1265 +        if (!(*_node_data)[ni].heap.empty()) {
  1.1266 +          _blossom_set->decrease(n, (*_node_data)[ni].heap.prio());
  1.1267 +          _delta2->push(nb, _blossom_set->classPrio(nb));
  1.1268 +        }
  1.1269 +      }
  1.1270 +    }
  1.1271 +
  1.1272      /// \brief Start the algorithm
  1.1273      ///
  1.1274      /// This function starts the algorithm.
  1.1275      ///
  1.1276 -    /// \pre \ref init() must be called before using this function.
  1.1277 +    /// \pre \ref init() or \ref fractionalInit() must be called before
  1.1278 +    /// using this function.
  1.1279      bool start() {
  1.1280        enum OpType {
  1.1281          D2, D3, D4
  1.1282        };
  1.1283  
  1.1284 -      int unmatched = _node_num;
  1.1285 -      while (unmatched > 0) {
  1.1286 +      if (_unmatched == -1) return false;
  1.1287 +
  1.1288 +      while (_unmatched > 0) {
  1.1289          Value d2 = !_delta2->empty() ?
  1.1290            _delta2->prio() : std::numeric_limits<Value>::max();
  1.1291  
  1.1292 @@ -2991,8 +3169,8 @@
  1.1293          Value d4 = !_delta4->empty() ?
  1.1294            _delta4->prio() : std::numeric_limits<Value>::max();
  1.1295  
  1.1296 -        _delta_sum = d2; OpType ot = D2;
  1.1297 -        if (d3 < _delta_sum) { _delta_sum = d3; ot = D3; }
  1.1298 +        _delta_sum = d3; OpType ot = D3;
  1.1299 +        if (d2 < _delta_sum) { _delta_sum = d2; ot = D2; }
  1.1300          if (d4 < _delta_sum) { _delta_sum = d4; ot = D4; }
  1.1301  
  1.1302          if (_delta_sum == std::numeric_limits<Value>::max()) {
  1.1303 @@ -3025,7 +3203,7 @@
  1.1304                  shrinkOnEdge(e, left_tree);
  1.1305                } else {
  1.1306                  augmentOnEdge(e);
  1.1307 -                unmatched -= 2;
  1.1308 +                _unmatched -= 2;
  1.1309                }
  1.1310              }
  1.1311            } break;
  1.1312 @@ -3044,18 +3222,18 @@
  1.1313      ///
  1.1314      /// \note mwpm.run() is just a shortcut of the following code.
  1.1315      /// \code
  1.1316 -    ///   mwpm.init();
  1.1317 +    ///   mwpm.fractionalInit();
  1.1318      ///   mwpm.start();
  1.1319      /// \endcode
  1.1320      bool run() {
  1.1321 -      init();
  1.1322 +      fractionalInit();
  1.1323        return start();
  1.1324      }
  1.1325  
  1.1326      /// @}
  1.1327  
  1.1328      /// \name Primal Solution
  1.1329 -    /// Functions to get the primal solution, i.e. the maximum weighted 
  1.1330 +    /// Functions to get the primal solution, i.e. the maximum weighted
  1.1331      /// perfect matching.\n
  1.1332      /// Either \ref run() or \ref start() function should be called before
  1.1333      /// using them.
  1.1334 @@ -3074,12 +3252,12 @@
  1.1335            sum += _weight[(*_matching)[n]];
  1.1336          }
  1.1337        }
  1.1338 -      return sum /= 2;
  1.1339 +      return sum / 2;
  1.1340      }
  1.1341  
  1.1342      /// \brief Return \c true if the given edge is in the matching.
  1.1343      ///
  1.1344 -    /// This function returns \c true if the given edge is in the found 
  1.1345 +    /// This function returns \c true if the given edge is in the found
  1.1346      /// matching.
  1.1347      ///
  1.1348      /// \pre Either run() or start() must be called before using this function.
  1.1349 @@ -3090,7 +3268,7 @@
  1.1350      /// \brief Return the matching arc (or edge) incident to the given node.
  1.1351      ///
  1.1352      /// This function returns the matching arc (or edge) incident to the
  1.1353 -    /// given node in the found matching or \c INVALID if the node is 
  1.1354 +    /// given node in the found matching or \c INVALID if the node is
  1.1355      /// not covered by the matching.
  1.1356      ///
  1.1357      /// \pre Either run() or start() must be called before using this function.
  1.1358 @@ -3108,7 +3286,7 @@
  1.1359  
  1.1360      /// \brief Return the mate of the given node.
  1.1361      ///
  1.1362 -    /// This function returns the mate of the given node in the found 
  1.1363 +    /// This function returns the mate of the given node in the found
  1.1364      /// matching or \c INVALID if the node is not covered by the matching.
  1.1365      ///
  1.1366      /// \pre Either run() or start() must be called before using this function.
  1.1367 @@ -3127,8 +3305,8 @@
  1.1368  
  1.1369      /// \brief Return the value of the dual solution.
  1.1370      ///
  1.1371 -    /// This function returns the value of the dual solution. 
  1.1372 -    /// It should be equal to the primal value scaled by \ref dualScale 
  1.1373 +    /// This function returns the value of the dual solution.
  1.1374 +    /// It should be equal to the primal value scaled by \ref dualScale
  1.1375      /// "dual scale".
  1.1376      ///
  1.1377      /// \pre Either run() or start() must be called before using this function.
  1.1378 @@ -3183,9 +3361,9 @@
  1.1379  
  1.1380      /// \brief Iterator for obtaining the nodes of a blossom.
  1.1381      ///
  1.1382 -    /// This class provides an iterator for obtaining the nodes of the 
  1.1383 +    /// This class provides an iterator for obtaining the nodes of the
  1.1384      /// given blossom. It lists a subset of the nodes.
  1.1385 -    /// Before using this iterator, you must allocate a 
  1.1386 +    /// Before using this iterator, you must allocate a
  1.1387      /// MaxWeightedPerfectMatching class and execute it.
  1.1388      class BlossomIt {
  1.1389      public:
  1.1390 @@ -3194,8 +3372,8 @@
  1.1391        ///
  1.1392        /// Constructor to get the nodes of the given variable.
  1.1393        ///
  1.1394 -      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()" 
  1.1395 -      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()" 
  1.1396 +      /// \pre Either \ref MaxWeightedPerfectMatching::run() "algorithm.run()"
  1.1397 +      /// or \ref MaxWeightedPerfectMatching::start() "algorithm.start()"
  1.1398        /// must be called before initializing this iterator.
  1.1399        BlossomIt(const MaxWeightedPerfectMatching& algorithm, int variable)
  1.1400          : _algorithm(&algorithm)
  1.1401 @@ -3241,4 +3419,4 @@
  1.1402  
  1.1403  } //END OF NAMESPACE LEMON
  1.1404  
  1.1405 -#endif //LEMON_MAX_MATCHING_H
  1.1406 +#endif //LEMON_MATCHING_H