Several improvements in maximum matching algorithms
authorBalazs Dezso <deba@inf.elte.hu>
Mon, 13 Oct 2008 14:00:11 +0200
changeset 33991d63b8b1a4c
parent 338 64ad48007fb2
child 342 5ba887b7def4
Several improvements in maximum matching algorithms
- The interface of MaxMatching is changed to be similar to the
weighted algorithms
- The internal data structure (the queue implementation and the
matching map) is changed in the MaxMatching algorithm, which
provides better runtime properties
- The Blossom iterators are changed slightly in the weighted matching
algorithms
- Several documentation improvments
- The test files are merged
lemon/max_matching.h
test/CMakeLists.txt
test/Makefile.am
test/max_matching_test.cc
test/max_weighted_matching_test.cc
     1.1 --- a/lemon/max_matching.h	Mon Oct 13 13:56:00 2008 +0200
     1.2 +++ b/lemon/max_matching.h	Mon Oct 13 14:00:11 2008 +0200
     1.3 @@ -31,76 +31,405 @@
     1.4  
     1.5  ///\ingroup matching
     1.6  ///\file
     1.7 -///\brief Maximum matching algorithms in graph.
     1.8 +///\brief Maximum matching algorithms in general graphs.
     1.9  
    1.10  namespace lemon {
    1.11  
    1.12 -  ///\ingroup matching
    1.13 +  /// \ingroup matching
    1.14    ///
    1.15 -  ///\brief Edmonds' alternating forest maximum matching algorithm.
    1.16 +  /// \brief Edmonds' alternating forest maximum matching algorithm.
    1.17    ///
    1.18 -  ///This class provides Edmonds' alternating forest matching
    1.19 -  ///algorithm. The starting matching (if any) can be passed to the
    1.20 -  ///algorithm using some of init functions.
    1.21 +  /// This class provides Edmonds' alternating forest matching
    1.22 +  /// algorithm. The starting matching (if any) can be passed to the
    1.23 +  /// algorithm using some of init functions.
    1.24    ///
    1.25 -  ///The dual side of a matching is a map of the nodes to
    1.26 -  ///MaxMatching::DecompType, having values \c D, \c A and \c C
    1.27 -  ///showing the Gallai-Edmonds decomposition of the digraph. The nodes
    1.28 -  ///in \c D induce a digraph with factor-critical components, the nodes
    1.29 -  ///in \c A form the barrier, and the nodes in \c C induce a digraph
    1.30 -  ///having a perfect matching. This decomposition can be attained by
    1.31 -  ///calling \c decomposition() after running the algorithm.
    1.32 +  /// The dual side of a matching is a map of the nodes to
    1.33 +  /// MaxMatching::Status, having values \c EVEN/D, \c ODD/A and \c
    1.34 +  /// MATCHED/C showing the Gallai-Edmonds decomposition of the
    1.35 +  /// graph. The nodes in \c EVEN/D induce a graph with
    1.36 +  /// factor-critical components, the nodes in \c ODD/A form the
    1.37 +  /// barrier, and the nodes in \c MATCHED/C induce a graph having a
    1.38 +  /// perfect matching. The number of the fractor critical components
    1.39 +  /// minus the number of barrier nodes is a lower bound on the
    1.40 +  /// unmatched nodes, and if the matching is optimal this bound is
    1.41 +  /// tight. This decomposition can be attained by calling \c
    1.42 +  /// decomposition() after running the algorithm.
    1.43    ///
    1.44 -  ///\param Digraph The graph type the algorithm runs on.
    1.45 -  template <typename Graph>
    1.46 +  /// \param _Graph The graph type the algorithm runs on.
    1.47 +  template <typename _Graph>
    1.48    class MaxMatching {
    1.49 -
    1.50 -  protected:
    1.51 +  public:
    1.52 +
    1.53 +    typedef _Graph Graph;
    1.54 +    typedef typename Graph::template NodeMap<typename Graph::Arc>
    1.55 +    MatchingMap;
    1.56 +
    1.57 +    ///\brief Indicates the Gallai-Edmonds decomposition of the graph.
    1.58 +    ///
    1.59 +    ///Indicates the Gallai-Edmonds decomposition of the graph, which
    1.60 +    ///shows an upper bound on the size of a maximum matching. The
    1.61 +    ///nodes with Status \c EVEN/D induce a graph with factor-critical
    1.62 +    ///components, the nodes in \c ODD/A form the canonical barrier,
    1.63 +    ///and the nodes in \c MATCHED/C induce a graph having a perfect
    1.64 +    ///matching.
    1.65 +    enum Status {
    1.66 +      EVEN = 1, D = 1, MATCHED = 0, C = 0, ODD = -1, A = -1, UNMATCHED = -2
    1.67 +    };
    1.68 +
    1.69 +    typedef typename Graph::template NodeMap<Status> StatusMap;
    1.70 +
    1.71 +  private:
    1.72  
    1.73      TEMPLATE_GRAPH_TYPEDEFS(Graph);
    1.74  
    1.75 -    typedef typename Graph::template NodeMap<int> UFECrossRef;
    1.76 -    typedef UnionFindEnum<UFECrossRef> UFE;
    1.77 -    typedef std::vector<Node> NV;
    1.78 -
    1.79 -    typedef typename Graph::template NodeMap<int> EFECrossRef;
    1.80 -    typedef ExtendFindEnum<EFECrossRef> EFE;
    1.81 +    typedef UnionFindEnum<IntNodeMap> BlossomSet;
    1.82 +    typedef ExtendFindEnum<IntNodeMap> TreeSet;
    1.83 +    typedef RangeMap<Node> NodeIntMap;
    1.84 +    typedef MatchingMap EarMap;
    1.85 +    typedef std::vector<Node> NodeQueue;
    1.86 +
    1.87 +    const Graph& _graph;
    1.88 +    MatchingMap* _matching;
    1.89 +    StatusMap* _status;
    1.90 +
    1.91 +    EarMap* _ear;
    1.92 +
    1.93 +    IntNodeMap* _blossom_set_index;
    1.94 +    BlossomSet* _blossom_set;
    1.95 +    NodeIntMap* _blossom_rep;
    1.96 +
    1.97 +    IntNodeMap* _tree_set_index;
    1.98 +    TreeSet* _tree_set;
    1.99 +
   1.100 +    NodeQueue _node_queue;
   1.101 +    int _process, _postpone, _last;
   1.102 +
   1.103 +    int _node_num;
   1.104 +
   1.105 +  private:
   1.106 +
   1.107 +    void createStructures() {
   1.108 +      _node_num = countNodes(_graph);
   1.109 +      if (!_matching) {
   1.110 +        _matching = new MatchingMap(_graph);
   1.111 +      }
   1.112 +      if (!_status) {
   1.113 +        _status = new StatusMap(_graph);
   1.114 +      }
   1.115 +      if (!_ear) {
   1.116 +        _ear = new EarMap(_graph);
   1.117 +      }
   1.118 +      if (!_blossom_set) {
   1.119 +        _blossom_set_index = new IntNodeMap(_graph);
   1.120 +        _blossom_set = new BlossomSet(*_blossom_set_index);
   1.121 +      }
   1.122 +      if (!_blossom_rep) {
   1.123 +        _blossom_rep = new NodeIntMap(_node_num);
   1.124 +      }
   1.125 +      if (!_tree_set) {
   1.126 +        _tree_set_index = new IntNodeMap(_graph);
   1.127 +        _tree_set = new TreeSet(*_tree_set_index);
   1.128 +      }
   1.129 +      _node_queue.resize(_node_num);
   1.130 +    }
   1.131 +
   1.132 +    void destroyStructures() {
   1.133 +      if (_matching) {
   1.134 +        delete _matching;
   1.135 +      }
   1.136 +      if (_status) {
   1.137 +        delete _status;
   1.138 +      }
   1.139 +      if (_ear) {
   1.140 +        delete _ear;
   1.141 +      }
   1.142 +      if (_blossom_set) {
   1.143 +        delete _blossom_set;
   1.144 +        delete _blossom_set_index;
   1.145 +      }
   1.146 +      if (_blossom_rep) {
   1.147 +        delete _blossom_rep;
   1.148 +      }
   1.149 +      if (_tree_set) {
   1.150 +        delete _tree_set_index;
   1.151 +        delete _tree_set;
   1.152 +      }
   1.153 +    }
   1.154 +
   1.155 +    void processDense(const Node& n) {
   1.156 +      _process = _postpone = _last = 0;
   1.157 +      _node_queue[_last++] = n;
   1.158 +
   1.159 +      while (_process != _last) {
   1.160 +        Node u = _node_queue[_process++];
   1.161 +        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
   1.162 +          Node v = _graph.target(a);
   1.163 +          if ((*_status)[v] == MATCHED) {
   1.164 +            extendOnArc(a);
   1.165 +          } else if ((*_status)[v] == UNMATCHED) {
   1.166 +            augmentOnArc(a);
   1.167 +            return;
   1.168 +          }
   1.169 +        }
   1.170 +      }
   1.171 +
   1.172 +      while (_postpone != _last) {
   1.173 +        Node u = _node_queue[_postpone++];
   1.174 +
   1.175 +        for (OutArcIt a(_graph, u); a != INVALID ; ++a) {
   1.176 +          Node v = _graph.target(a);
   1.177 +
   1.178 +          if ((*_status)[v] == EVEN) {
   1.179 +            if (_blossom_set->find(u) != _blossom_set->find(v)) {
   1.180 +              shrinkOnEdge(a);
   1.181 +            }
   1.182 +          }
   1.183 +
   1.184 +          while (_process != _last) {
   1.185 +            Node w = _node_queue[_process++];
   1.186 +            for (OutArcIt b(_graph, w); b != INVALID; ++b) {
   1.187 +              Node x = _graph.target(b);
   1.188 +              if ((*_status)[x] == MATCHED) {
   1.189 +                extendOnArc(b);
   1.190 +              } else if ((*_status)[x] == UNMATCHED) {
   1.191 +                augmentOnArc(b);
   1.192 +                return;
   1.193 +              }
   1.194 +            }
   1.195 +          }
   1.196 +        }
   1.197 +      }
   1.198 +    }
   1.199 +
   1.200 +    void processSparse(const Node& n) {
   1.201 +      _process = _last = 0;
   1.202 +      _node_queue[_last++] = n;
   1.203 +      while (_process != _last) {
   1.204 +        Node u = _node_queue[_process++];
   1.205 +        for (OutArcIt a(_graph, u); a != INVALID; ++a) {
   1.206 +          Node v = _graph.target(a);
   1.207 +
   1.208 +          if ((*_status)[v] == EVEN) {
   1.209 +            if (_blossom_set->find(u) != _blossom_set->find(v)) {
   1.210 +              shrinkOnEdge(a);
   1.211 +            }
   1.212 +          } else if ((*_status)[v] == MATCHED) {
   1.213 +            extendOnArc(a);
   1.214 +          } else if ((*_status)[v] == UNMATCHED) {
   1.215 +            augmentOnArc(a);
   1.216 +            return;
   1.217 +          }
   1.218 +        }
   1.219 +      }
   1.220 +    }
   1.221 +
   1.222 +    void shrinkOnEdge(const Edge& e) {
   1.223 +      Node nca = INVALID;
   1.224 +
   1.225 +      {
   1.226 +        std::set<Node> left_set, right_set;
   1.227 +
   1.228 +        Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];
   1.229 +        left_set.insert(left);
   1.230 +
   1.231 +        Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];
   1.232 +        right_set.insert(right);
   1.233 +
   1.234 +        while (true) {
   1.235 +          if ((*_matching)[left] == INVALID) break;
   1.236 +          left = _graph.target((*_matching)[left]);
   1.237 +          left = (*_blossom_rep)[_blossom_set->
   1.238 +                                 find(_graph.target((*_ear)[left]))];
   1.239 +          if (right_set.find(left) != right_set.end()) {
   1.240 +            nca = left;
   1.241 +            break;
   1.242 +          }
   1.243 +          left_set.insert(left);
   1.244 +
   1.245 +          if ((*_matching)[right] == INVALID) break;
   1.246 +          right = _graph.target((*_matching)[right]);
   1.247 +          right = (*_blossom_rep)[_blossom_set->
   1.248 +                                  find(_graph.target((*_ear)[right]))];
   1.249 +          if (left_set.find(right) != left_set.end()) {
   1.250 +            nca = right;
   1.251 +            break;
   1.252 +          }
   1.253 +          right_set.insert(right);
   1.254 +        }
   1.255 +
   1.256 +        if (nca == INVALID) {
   1.257 +          if ((*_matching)[left] == INVALID) {
   1.258 +            nca = right;
   1.259 +            while (left_set.find(nca) == left_set.end()) {
   1.260 +              nca = _graph.target((*_matching)[nca]);
   1.261 +              nca =(*_blossom_rep)[_blossom_set->
   1.262 +                                   find(_graph.target((*_ear)[nca]))];
   1.263 +            }
   1.264 +          } else {
   1.265 +            nca = left;
   1.266 +            while (right_set.find(nca) == right_set.end()) {
   1.267 +              nca = _graph.target((*_matching)[nca]);
   1.268 +              nca = (*_blossom_rep)[_blossom_set->
   1.269 +                                   find(_graph.target((*_ear)[nca]))];
   1.270 +            }
   1.271 +          }
   1.272 +        }
   1.273 +      }
   1.274 +
   1.275 +      {
   1.276 +
   1.277 +        Node node = _graph.u(e);
   1.278 +        Arc arc = _graph.direct(e, true);
   1.279 +        Node base = (*_blossom_rep)[_blossom_set->find(node)];
   1.280 +
   1.281 +        while (base != nca) {
   1.282 +          _ear->set(node, arc);
   1.283 +
   1.284 +          Node n = node;
   1.285 +          while (n != base) {
   1.286 +            n = _graph.target((*_matching)[n]);
   1.287 +            Arc a = (*_ear)[n];
   1.288 +            n = _graph.target(a);
   1.289 +            _ear->set(n, _graph.oppositeArc(a));
   1.290 +          }
   1.291 +          node = _graph.target((*_matching)[base]);
   1.292 +          _tree_set->erase(base);
   1.293 +          _tree_set->erase(node);
   1.294 +          _blossom_set->insert(node, _blossom_set->find(base));
   1.295 +          _status->set(node, EVEN);
   1.296 +          _node_queue[_last++] = node;
   1.297 +          arc = _graph.oppositeArc((*_ear)[node]);
   1.298 +          node = _graph.target((*_ear)[node]);
   1.299 +          base = (*_blossom_rep)[_blossom_set->find(node)];
   1.300 +          _blossom_set->join(_graph.target(arc), base);
   1.301 +        }
   1.302 +      }
   1.303 +
   1.304 +      _blossom_rep->set(_blossom_set->find(nca), nca);
   1.305 +
   1.306 +      {
   1.307 +
   1.308 +        Node node = _graph.v(e);
   1.309 +        Arc arc = _graph.direct(e, false);
   1.310 +        Node base = (*_blossom_rep)[_blossom_set->find(node)];
   1.311 +
   1.312 +        while (base != nca) {
   1.313 +          _ear->set(node, arc);
   1.314 +
   1.315 +          Node n = node;
   1.316 +          while (n != base) {
   1.317 +            n = _graph.target((*_matching)[n]);
   1.318 +            Arc a = (*_ear)[n];
   1.319 +            n = _graph.target(a);
   1.320 +            _ear->set(n, _graph.oppositeArc(a));
   1.321 +          }
   1.322 +          node = _graph.target((*_matching)[base]);
   1.323 +          _tree_set->erase(base);
   1.324 +          _tree_set->erase(node);
   1.325 +          _blossom_set->insert(node, _blossom_set->find(base));
   1.326 +          _status->set(node, EVEN);
   1.327 +          _node_queue[_last++] = node;
   1.328 +          arc = _graph.oppositeArc((*_ear)[node]);
   1.329 +          node = _graph.target((*_ear)[node]);
   1.330 +          base = (*_blossom_rep)[_blossom_set->find(node)];
   1.331 +          _blossom_set->join(_graph.target(arc), base);
   1.332 +        }
   1.333 +      }
   1.334 +
   1.335 +      _blossom_rep->set(_blossom_set->find(nca), nca);
   1.336 +    }
   1.337 +
   1.338 +
   1.339 +
   1.340 +    void extendOnArc(const Arc& a) {
   1.341 +      Node base = _graph.source(a);
   1.342 +      Node odd = _graph.target(a);
   1.343 +
   1.344 +      _ear->set(odd, _graph.oppositeArc(a));
   1.345 +      Node even = _graph.target((*_matching)[odd]);
   1.346 +      _blossom_rep->set(_blossom_set->insert(even), even);
   1.347 +      _status->set(odd, ODD);
   1.348 +      _status->set(even, EVEN);
   1.349 +      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);
   1.350 +      _tree_set->insert(odd, tree);
   1.351 +      _tree_set->insert(even, tree);
   1.352 +      _node_queue[_last++] = even;
   1.353 +
   1.354 +    }
   1.355 +
   1.356 +    void augmentOnArc(const Arc& a) {
   1.357 +      Node even = _graph.source(a);
   1.358 +      Node odd = _graph.target(a);
   1.359 +
   1.360 +      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);
   1.361 +
   1.362 +      _matching->set(odd, _graph.oppositeArc(a));
   1.363 +      _status->set(odd, MATCHED);
   1.364 +
   1.365 +      Arc arc = (*_matching)[even];
   1.366 +      _matching->set(even, a);
   1.367 +
   1.368 +      while (arc != INVALID) {
   1.369 +        odd = _graph.target(arc);
   1.370 +        arc = (*_ear)[odd];
   1.371 +        even = _graph.target(arc);
   1.372 +        _matching->set(odd, arc);
   1.373 +        arc = (*_matching)[even];
   1.374 +        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));
   1.375 +      }
   1.376 +
   1.377 +      for (typename TreeSet::ItemIt it(*_tree_set, tree);
   1.378 +           it != INVALID; ++it) {
   1.379 +        if ((*_status)[it] == ODD) {
   1.380 +          _status->set(it, MATCHED);
   1.381 +        } else {
   1.382 +          int blossom = _blossom_set->find(it);
   1.383 +          for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);
   1.384 +               jt != INVALID; ++jt) {
   1.385 +            _status->set(jt, MATCHED);
   1.386 +          }
   1.387 +          _blossom_set->eraseClass(blossom);
   1.388 +        }
   1.389 +      }
   1.390 +      _tree_set->eraseClass(tree);
   1.391 +
   1.392 +    }
   1.393  
   1.394    public:
   1.395  
   1.396 -    ///\brief Indicates the Gallai-Edmonds decomposition of the digraph.
   1.397 +    /// \brief Constructor
   1.398      ///
   1.399 -    ///Indicates the Gallai-Edmonds decomposition of the digraph, which
   1.400 -    ///shows an upper bound on the size of a maximum matching. The
   1.401 -    ///nodes with DecompType \c D induce a digraph with factor-critical
   1.402 -    ///components, the nodes in \c A form the canonical barrier, and the
   1.403 -    ///nodes in \c C induce a digraph having a perfect matching.
   1.404 -    enum DecompType {
   1.405 -      D=0,
   1.406 -      A=1,
   1.407 -      C=2
   1.408 -    };
   1.409 -
   1.410 -  protected:
   1.411 -
   1.412 -    static const int HEUR_density=2;
   1.413 -    const Graph& g;
   1.414 -    typename Graph::template NodeMap<Node> _mate;
   1.415 -    typename Graph::template NodeMap<DecompType> position;
   1.416 -
   1.417 -  public:
   1.418 -
   1.419 -    MaxMatching(const Graph& _g)
   1.420 -      : g(_g), _mate(_g), position(_g) {}
   1.421 -
   1.422 -    ///\brief Sets the actual matching to the empty matching.
   1.423 +    /// Constructor.
   1.424 +    MaxMatching(const Graph& graph)
   1.425 +      : _graph(graph), _matching(0), _status(0), _ear(0),
   1.426 +        _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),
   1.427 +        _tree_set_index(0), _tree_set(0) {}
   1.428 +
   1.429 +    ~MaxMatching() {
   1.430 +      destroyStructures();
   1.431 +    }
   1.432 +
   1.433 +    /// \name Execution control
   1.434 +    /// The simplest way to execute the algorithm is to use the member
   1.435 +    /// \c run() member function.
   1.436 +    /// \n
   1.437 +
   1.438 +    /// If you need more control on the execution, first you must call
   1.439 +    /// \ref init(), \ref greedyInit() or \ref matchingInit()
   1.440 +    /// functions, then you can start the algorithm with the \ref
   1.441 +    /// startParse() or startDense() functions.
   1.442 +
   1.443 +    ///@{
   1.444 +
   1.445 +    /// \brief Sets the actual matching to the empty matching.
   1.446      ///
   1.447 -    ///Sets the actual matching to the empty matching.
   1.448 +    /// Sets the actual matching to the empty matching.
   1.449      ///
   1.450      void init() {
   1.451 -      for(NodeIt v(g); v!=INVALID; ++v) {
   1.452 -        _mate.set(v,INVALID);
   1.453 -        position.set(v,C);
   1.454 +      createStructures();
   1.455 +      for(NodeIt n(_graph); n != INVALID; ++n) {
   1.456 +        _matching->set(n, INVALID);
   1.457 +        _status->set(n, UNMATCHED);
   1.458        }
   1.459      }
   1.460  
   1.461 @@ -108,88 +437,96 @@
   1.462      ///
   1.463      ///For initial matchig it finds a maximal greedy matching.
   1.464      void greedyInit() {
   1.465 -      for(NodeIt v(g); v!=INVALID; ++v) {
   1.466 -        _mate.set(v,INVALID);
   1.467 -        position.set(v,C);
   1.468 +      createStructures();
   1.469 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.470 +        _matching->set(n, INVALID);
   1.471 +        _status->set(n, UNMATCHED);
   1.472        }
   1.473 -      for(NodeIt v(g); v!=INVALID; ++v)
   1.474 -        if ( _mate[v]==INVALID ) {
   1.475 -          for( IncEdgeIt e(g,v); e!=INVALID ; ++e ) {
   1.476 -            Node y=g.runningNode(e);
   1.477 -            if ( _mate[y]==INVALID && y!=v ) {
   1.478 -              _mate.set(v,y);
   1.479 -              _mate.set(y,v);
   1.480 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.481 +        if ((*_matching)[n] == INVALID) {
   1.482 +          for (OutArcIt a(_graph, n); a != INVALID ; ++a) {
   1.483 +            Node v = _graph.target(a);
   1.484 +            if ((*_matching)[v] == INVALID && v != n) {
   1.485 +              _matching->set(n, a);
   1.486 +              _status->set(n, MATCHED);
   1.487 +              _matching->set(v, _graph.oppositeArc(a));
   1.488 +              _status->set(v, MATCHED);
   1.489                break;
   1.490              }
   1.491            }
   1.492          }
   1.493 -    }
   1.494 -
   1.495 -    ///\brief Initialize the matching from each nodes' mate.
   1.496 -    ///
   1.497 -    ///Initialize the matching from a \c Node valued \c Node map. This
   1.498 -    ///map must be \e symmetric, i.e. if \c map[u]==v then \c
   1.499 -    ///map[v]==u must hold, and \c uv will be an arc of the initial
   1.500 -    ///matching.
   1.501 -    template <typename MateMap>
   1.502 -    void mateMapInit(MateMap& map) {
   1.503 -      for(NodeIt v(g); v!=INVALID; ++v) {
   1.504 -        _mate.set(v,map[v]);
   1.505 -        position.set(v,C);
   1.506        }
   1.507      }
   1.508  
   1.509 -    ///\brief Initialize the matching from a node map with the
   1.510 -    ///incident matching arcs.
   1.511 +
   1.512 +    /// \brief Initialize the matching from the map containing a matching.
   1.513      ///
   1.514 -    ///Initialize the matching from an \c Edge valued \c Node map. \c
   1.515 -    ///map[v] must be an \c Edge incident to \c v. This map must have
   1.516 -    ///the property that if \c g.oppositeNode(u,map[u])==v then \c \c
   1.517 -    ///g.oppositeNode(v,map[v])==u holds, and now some arc joining \c
   1.518 -    ///u to \c v will be an arc of the matching.
   1.519 -    template<typename MatchingMap>
   1.520 -    void matchingMapInit(MatchingMap& map) {
   1.521 -      for(NodeIt v(g); v!=INVALID; ++v) {
   1.522 -        position.set(v,C);
   1.523 -        Edge e=map[v];
   1.524 -        if ( e!=INVALID )
   1.525 -          _mate.set(v,g.oppositeNode(v,e));
   1.526 -        else
   1.527 -          _mate.set(v,INVALID);
   1.528 +    /// Initialize the matching from a \c bool valued \c Edge map. This
   1.529 +    /// map must have the property that there are no two incident edges
   1.530 +    /// with true value, ie. it contains a matching.
   1.531 +    /// \return %True if the map contains a matching.
   1.532 +    template <typename MatchingMap>
   1.533 +    bool matchingInit(const MatchingMap& matching) {
   1.534 +      createStructures();
   1.535 +
   1.536 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.537 +        _matching->set(n, INVALID);
   1.538 +        _status->set(n, UNMATCHED);
   1.539        }
   1.540 +      for(EdgeIt e(_graph); e!=INVALID; ++e) {
   1.541 +        if (matching[e]) {
   1.542 +
   1.543 +          Node u = _graph.u(e);
   1.544 +          if ((*_matching)[u] != INVALID) return false;
   1.545 +          _matching->set(u, _graph.direct(e, true));
   1.546 +          _status->set(u, MATCHED);
   1.547 +
   1.548 +          Node v = _graph.v(e);
   1.549 +          if ((*_matching)[v] != INVALID) return false;
   1.550 +          _matching->set(v, _graph.direct(e, false));
   1.551 +          _status->set(v, MATCHED);
   1.552 +        }
   1.553 +      }
   1.554 +      return true;
   1.555      }
   1.556  
   1.557 -    ///\brief Initialize the matching from the map containing the
   1.558 -    ///undirected matching arcs.
   1.559 +    /// \brief Starts Edmonds' algorithm
   1.560      ///
   1.561 -    ///Initialize the matching from a \c bool valued \c Edge map. This
   1.562 -    ///map must have the property that there are no two incident arcs
   1.563 -    ///\c e, \c f with \c map[e]==map[f]==true. The arcs \c e with \c
   1.564 -    ///map[e]==true form the matching.
   1.565 -    template <typename MatchingMap>
   1.566 -    void matchingInit(MatchingMap& map) {
   1.567 -      for(NodeIt v(g); v!=INVALID; ++v) {
   1.568 -        _mate.set(v,INVALID);
   1.569 -        position.set(v,C);
   1.570 -      }
   1.571 -      for(EdgeIt e(g); e!=INVALID; ++e) {
   1.572 -        if ( map[e] ) {
   1.573 -          Node u=g.u(e);
   1.574 -          Node v=g.v(e);
   1.575 -          _mate.set(u,v);
   1.576 -          _mate.set(v,u);
   1.577 +    /// If runs the original Edmonds' algorithm.
   1.578 +    void startSparse() {
   1.579 +      for(NodeIt n(_graph); n != INVALID; ++n) {
   1.580 +        if ((*_status)[n] == UNMATCHED) {
   1.581 +          (*_blossom_rep)[_blossom_set->insert(n)] = n;
   1.582 +          _tree_set->insert(n);
   1.583 +          _status->set(n, EVEN);
   1.584 +          processSparse(n);
   1.585          }
   1.586        }
   1.587      }
   1.588  
   1.589 -
   1.590 -    ///\brief Runs Edmonds' algorithm.
   1.591 +    /// \brief Starts Edmonds' algorithm.
   1.592      ///
   1.593 -    ///Runs Edmonds' algorithm for sparse digraphs (number of arcs <
   1.594 -    ///2*number of nodes), and a heuristical Edmonds' algorithm with a
   1.595 -    ///heuristic of postponing shrinks for dense digraphs.
   1.596 +    /// It runs Edmonds' algorithm with a heuristic of postponing
   1.597 +    /// shrinks, giving a faster algorithm for dense graphs.
   1.598 +    void startDense() {
   1.599 +      for(NodeIt n(_graph); n != INVALID; ++n) {
   1.600 +        if ((*_status)[n] == UNMATCHED) {
   1.601 +          (*_blossom_rep)[_blossom_set->insert(n)] = n;
   1.602 +          _tree_set->insert(n);
   1.603 +          _status->set(n, EVEN);
   1.604 +          processDense(n);
   1.605 +        }
   1.606 +      }
   1.607 +    }
   1.608 +
   1.609 +
   1.610 +    /// \brief Runs Edmonds' algorithm
   1.611 +    ///
   1.612 +    /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>)
   1.613 +    /// or Edmonds' algorithm with a heuristic of
   1.614 +    /// postponing shrinks for dense graphs.
   1.615      void run() {
   1.616 -      if (countEdges(g) < HEUR_density * countNodes(g)) {
   1.617 +      if (countEdges(_graph) < 2 * countNodes(_graph)) {
   1.618          greedyInit();
   1.619          startSparse();
   1.620        } else {
   1.621 @@ -198,422 +535,75 @@
   1.622        }
   1.623      }
   1.624  
   1.625 -
   1.626 -    ///\brief Starts Edmonds' algorithm.
   1.627 -    ///
   1.628 -    ///If runs the original Edmonds' algorithm.
   1.629 -    void startSparse() {
   1.630 -
   1.631 -      typename Graph::template NodeMap<Node> ear(g,INVALID);
   1.632 -      //undefined for the base nodes of the blossoms (i.e. for the
   1.633 -      //representative elements of UFE blossom) and for the nodes in C
   1.634 -
   1.635 -      UFECrossRef blossom_base(g);
   1.636 -      UFE blossom(blossom_base);
   1.637 -      NV rep(countNodes(g));
   1.638 -
   1.639 -      EFECrossRef tree_base(g);
   1.640 -      EFE tree(tree_base);
   1.641 -
   1.642 -      //If these UFE's would be members of the class then also
   1.643 -      //blossom_base and tree_base should be a member.
   1.644 -
   1.645 -      //We build only one tree and the other vertices uncovered by the
   1.646 -      //matching belong to C. (They can be considered as singleton
   1.647 -      //trees.) If this tree can be augmented or no more
   1.648 -      //grow/augmentation/shrink is possible then we return to this
   1.649 -      //"for" cycle.
   1.650 -      for(NodeIt v(g); v!=INVALID; ++v) {
   1.651 -        if (position[v]==C && _mate[v]==INVALID) {
   1.652 -          rep[blossom.insert(v)] = v;
   1.653 -          tree.insert(v);
   1.654 -          position.set(v,D);
   1.655 -          normShrink(v, ear, blossom, rep, tree);
   1.656 -        }
   1.657 -      }
   1.658 -    }
   1.659 -
   1.660 -    ///\brief Starts Edmonds' algorithm.
   1.661 -    ///
   1.662 -    ///It runs Edmonds' algorithm with a heuristic of postponing
   1.663 -    ///shrinks, giving a faster algorithm for dense digraphs.
   1.664 -    void startDense() {
   1.665 -
   1.666 -      typename Graph::template NodeMap<Node> ear(g,INVALID);
   1.667 -      //undefined for the base nodes of the blossoms (i.e. for the
   1.668 -      //representative elements of UFE blossom) and for the nodes in C
   1.669 -
   1.670 -      UFECrossRef blossom_base(g);
   1.671 -      UFE blossom(blossom_base);
   1.672 -      NV rep(countNodes(g));
   1.673 -
   1.674 -      EFECrossRef tree_base(g);
   1.675 -      EFE tree(tree_base);
   1.676 -
   1.677 -      //If these UFE's would be members of the class then also
   1.678 -      //blossom_base and tree_base should be a member.
   1.679 -
   1.680 -      //We build only one tree and the other vertices uncovered by the
   1.681 -      //matching belong to C. (They can be considered as singleton
   1.682 -      //trees.) If this tree can be augmented or no more
   1.683 -      //grow/augmentation/shrink is possible then we return to this
   1.684 -      //"for" cycle.
   1.685 -      for(NodeIt v(g); v!=INVALID; ++v) {
   1.686 -        if ( position[v]==C && _mate[v]==INVALID ) {
   1.687 -          rep[blossom.insert(v)] = v;
   1.688 -          tree.insert(v);
   1.689 -          position.set(v,D);
   1.690 -          lateShrink(v, ear, blossom, rep, tree);
   1.691 -        }
   1.692 -      }
   1.693 -    }
   1.694 -
   1.695 -
   1.696 +    /// @}
   1.697 +
   1.698 +    /// \name Primal solution
   1.699 +    /// Functions for get the primal solution, ie. the matching.
   1.700 +
   1.701 +    /// @{
   1.702  
   1.703      ///\brief Returns the size of the actual matching stored.
   1.704      ///
   1.705      ///Returns the size of the actual matching stored. After \ref
   1.706 -    ///run() it returns the size of a maximum matching in the digraph.
   1.707 -    int size() const {
   1.708 -      int s=0;
   1.709 -      for(NodeIt v(g); v!=INVALID; ++v) {
   1.710 -        if ( _mate[v]!=INVALID ) {
   1.711 -          ++s;
   1.712 +    ///run() it returns the size of the maximum matching in the graph.
   1.713 +    int matchingSize() const {
   1.714 +      int size = 0;
   1.715 +      for (NodeIt n(_graph); n != INVALID; ++n) {
   1.716 +        if ((*_matching)[n] != INVALID) {
   1.717 +          ++size;
   1.718          }
   1.719        }
   1.720 -      return s/2;
   1.721 +      return size / 2;
   1.722      }
   1.723  
   1.724 +    /// \brief Returns true when the edge is in the matching.
   1.725 +    ///
   1.726 +    /// Returns true when the edge is in the matching.
   1.727 +    bool matching(const Edge& edge) const {
   1.728 +      return edge == (*_matching)[_graph.u(edge)];
   1.729 +    }
   1.730 +
   1.731 +    /// \brief Returns the matching edge incident to the given node.
   1.732 +    ///
   1.733 +    /// Returns the matching edge of a \c node in the actual matching or
   1.734 +    /// INVALID if the \c node is not covered by the actual matching.
   1.735 +    Arc matching(const Node& n) const {
   1.736 +      return (*_matching)[n];
   1.737 +    }
   1.738  
   1.739      ///\brief Returns the mate of a node in the actual matching.
   1.740      ///
   1.741 -    ///Returns the mate of a \c node in the actual matching.
   1.742 -    ///Returns INVALID if the \c node is not covered by the actual matching.
   1.743 -    Node mate(const Node& node) const {
   1.744 -      return _mate[node];
   1.745 +    ///Returns the mate of a \c node in the actual matching or
   1.746 +    ///INVALID if the \c node is not covered by the actual matching.
   1.747 +    Node mate(const Node& n) const {
   1.748 +      return (*_matching)[n] != INVALID ?
   1.749 +        _graph.target((*_matching)[n]) : INVALID;
   1.750      }
   1.751  
   1.752 -    ///\brief Returns the matching arc incident to the given node.
   1.753 -    ///
   1.754 -    ///Returns the matching arc of a \c node in the actual matching.
   1.755 -    ///Returns INVALID if the \c node is not covered by the actual matching.
   1.756 -    Edge matchingArc(const Node& node) const {
   1.757 -      if (_mate[node] == INVALID) return INVALID;
   1.758 -      Node n = node < _mate[node] ? node : _mate[node];
   1.759 -      for (IncEdgeIt e(g, n); e != INVALID; ++e) {
   1.760 -        if (g.oppositeNode(n, e) == _mate[n]) {
   1.761 -          return e;
   1.762 -        }
   1.763 -      }
   1.764 -      return INVALID;
   1.765 -    }
   1.766 +    /// @}
   1.767 +
   1.768 +    /// \name Dual solution
   1.769 +    /// Functions for get the dual solution, ie. the decomposition.
   1.770 +
   1.771 +    /// @{
   1.772  
   1.773      /// \brief Returns the class of the node in the Edmonds-Gallai
   1.774      /// decomposition.
   1.775      ///
   1.776      /// Returns the class of the node in the Edmonds-Gallai
   1.777      /// decomposition.
   1.778 -    DecompType decomposition(const Node& n) {
   1.779 -      return position[n] == A;
   1.780 +    Status decomposition(const Node& n) const {
   1.781 +      return (*_status)[n];
   1.782      }
   1.783  
   1.784      /// \brief Returns true when the node is in the barrier.
   1.785      ///
   1.786      /// Returns true when the node is in the barrier.
   1.787 -    bool barrier(const Node& n) {
   1.788 -      return position[n] == A;
   1.789 +    bool barrier(const Node& n) const {
   1.790 +      return (*_status)[n] == ODD;
   1.791      }
   1.792  
   1.793 -    ///\brief Gives back the matching in a \c Node of mates.
   1.794 -    ///
   1.795 -    ///Writes the stored matching to a \c Node valued \c Node map. The
   1.796 -    ///resulting map will be \e symmetric, i.e. if \c map[u]==v then \c
   1.797 -    ///map[v]==u will hold, and now \c uv is an arc of the matching.
   1.798 -    template <typename MateMap>
   1.799 -    void mateMap(MateMap& map) const {
   1.800 -      for(NodeIt v(g); v!=INVALID; ++v) {
   1.801 -        map.set(v,_mate[v]);
   1.802 -      }
   1.803 -    }
   1.804 -
   1.805 -    ///\brief Gives back the matching in an \c Edge valued \c Node
   1.806 -    ///map.
   1.807 -    ///
   1.808 -    ///Writes the stored matching to an \c Edge valued \c Node
   1.809 -    ///map. \c map[v] will be an \c Edge incident to \c v. This
   1.810 -    ///map will have the property that if \c g.oppositeNode(u,map[u])
   1.811 -    ///== v then \c map[u]==map[v] holds, and now this arc is an arc
   1.812 -    ///of the matching.
   1.813 -    template <typename MatchingMap>
   1.814 -    void matchingMap(MatchingMap& map)  const {
   1.815 -      typename Graph::template NodeMap<bool> todo(g,true);
   1.816 -      for(NodeIt v(g); v!=INVALID; ++v) {
   1.817 -        if (_mate[v]!=INVALID && v < _mate[v]) {
   1.818 -          Node u=_mate[v];
   1.819 -          for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
   1.820 -            if ( g.runningNode(e) == u ) {
   1.821 -              map.set(u,e);
   1.822 -              map.set(v,e);
   1.823 -              todo.set(u,false);
   1.824 -              todo.set(v,false);
   1.825 -              break;
   1.826 -            }
   1.827 -          }
   1.828 -        }
   1.829 -      }
   1.830 -    }
   1.831 -
   1.832 -
   1.833 -    ///\brief Gives back the matching in a \c bool valued \c Edge
   1.834 -    ///map.
   1.835 -    ///
   1.836 -    ///Writes the matching stored to a \c bool valued \c Arc
   1.837 -    ///map. This map will have the property that there are no two
   1.838 -    ///incident arcs \c e, \c f with \c map[e]==map[f]==true. The
   1.839 -    ///arcs \c e with \c map[e]==true form the matching.
   1.840 -    template<typename MatchingMap>
   1.841 -    void matching(MatchingMap& map) const {
   1.842 -      for(EdgeIt e(g); e!=INVALID; ++e) map.set(e,false);
   1.843 -
   1.844 -      typename Graph::template NodeMap<bool> todo(g,true);
   1.845 -      for(NodeIt v(g); v!=INVALID; ++v) {
   1.846 -        if ( todo[v] && _mate[v]!=INVALID ) {
   1.847 -          Node u=_mate[v];
   1.848 -          for(IncEdgeIt e(g,v); e!=INVALID; ++e) {
   1.849 -            if ( g.runningNode(e) == u ) {
   1.850 -              map.set(e,true);
   1.851 -              todo.set(u,false);
   1.852 -              todo.set(v,false);
   1.853 -              break;
   1.854 -            }
   1.855 -          }
   1.856 -        }
   1.857 -      }
   1.858 -    }
   1.859 -
   1.860 -
   1.861 -    ///\brief Returns the canonical decomposition of the digraph after running
   1.862 -    ///the algorithm.
   1.863 -    ///
   1.864 -    ///After calling any run methods of the class, it writes the
   1.865 -    ///Gallai-Edmonds canonical decomposition of the digraph. \c map
   1.866 -    ///must be a node map of \ref DecompType 's.
   1.867 -    template <typename DecompositionMap>
   1.868 -    void decomposition(DecompositionMap& map) const {
   1.869 -      for(NodeIt v(g); v!=INVALID; ++v) map.set(v,position[v]);
   1.870 -    }
   1.871 -
   1.872 -    ///\brief Returns a barrier on the nodes.
   1.873 -    ///
   1.874 -    ///After calling any run methods of the class, it writes a
   1.875 -    ///canonical barrier on the nodes. The odd component number of the
   1.876 -    ///remaining digraph minus the barrier size is a lower bound for the
   1.877 -    ///uncovered nodes in the digraph. The \c map must be a node map of
   1.878 -    ///bools.
   1.879 -    template <typename BarrierMap>
   1.880 -    void barrier(BarrierMap& barrier) {
   1.881 -      for(NodeIt v(g); v!=INVALID; ++v) barrier.set(v,position[v] == A);
   1.882 -    }
   1.883 -
   1.884 -  private:
   1.885 -
   1.886 -
   1.887 -    void lateShrink(Node v, typename Graph::template NodeMap<Node>& ear,
   1.888 -                    UFE& blossom, NV& rep, EFE& tree) {
   1.889 -      //We have one tree which we grow, and also shrink but only if it
   1.890 -      //cannot be postponed. If we augment then we return to the "for"
   1.891 -      //cycle of runEdmonds().
   1.892 -
   1.893 -      std::queue<Node> Q;   //queue of the totally unscanned nodes
   1.894 -      Q.push(v);
   1.895 -      std::queue<Node> R;
   1.896 -      //queue of the nodes which must be scanned for a possible shrink
   1.897 -
   1.898 -      while ( !Q.empty() ) {
   1.899 -        Node x=Q.front();
   1.900 -        Q.pop();
   1.901 -        for( IncEdgeIt e(g,x); e!= INVALID; ++e ) {
   1.902 -          Node y=g.runningNode(e);
   1.903 -          //growOrAugment grows if y is covered by the matching and
   1.904 -          //augments if not. In this latter case it returns 1.
   1.905 -          if (position[y]==C &&
   1.906 -              growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
   1.907 -        }
   1.908 -        R.push(x);
   1.909 -      }
   1.910 -
   1.911 -      while ( !R.empty() ) {
   1.912 -        Node x=R.front();
   1.913 -        R.pop();
   1.914 -
   1.915 -        for( IncEdgeIt e(g,x); e!=INVALID ; ++e ) {
   1.916 -          Node y=g.runningNode(e);
   1.917 -
   1.918 -          if ( position[y] == D && blossom.find(x) != blossom.find(y) )
   1.919 -            //Recall that we have only one tree.
   1.920 -            shrink( x, y, ear, blossom, rep, tree, Q);
   1.921 -
   1.922 -          while ( !Q.empty() ) {
   1.923 -            Node z=Q.front();
   1.924 -            Q.pop();
   1.925 -            for( IncEdgeIt f(g,z); f!= INVALID; ++f ) {
   1.926 -              Node w=g.runningNode(f);
   1.927 -              //growOrAugment grows if y is covered by the matching and
   1.928 -              //augments if not. In this latter case it returns 1.
   1.929 -              if (position[w]==C &&
   1.930 -                  growOrAugment(w, z, ear, blossom, rep, tree, Q)) return;
   1.931 -            }
   1.932 -            R.push(z);
   1.933 -          }
   1.934 -        } //for e
   1.935 -      } // while ( !R.empty() )
   1.936 -    }
   1.937 -
   1.938 -    void normShrink(Node v, typename Graph::template NodeMap<Node>& ear,
   1.939 -                    UFE& blossom, NV& rep, EFE& tree) {
   1.940 -      //We have one tree, which we grow and shrink. If we augment then we
   1.941 -      //return to the "for" cycle of runEdmonds().
   1.942 -
   1.943 -      std::queue<Node> Q;   //queue of the unscanned nodes
   1.944 -      Q.push(v);
   1.945 -      while ( !Q.empty() ) {
   1.946 -
   1.947 -        Node x=Q.front();
   1.948 -        Q.pop();
   1.949 -
   1.950 -        for( IncEdgeIt e(g,x); e!=INVALID; ++e ) {
   1.951 -          Node y=g.runningNode(e);
   1.952 -
   1.953 -          switch ( position[y] ) {
   1.954 -          case D:          //x and y must be in the same tree
   1.955 -            if ( blossom.find(x) != blossom.find(y))
   1.956 -              //x and y are in the same tree
   1.957 -              shrink(x, y, ear, blossom, rep, tree, Q);
   1.958 -            break;
   1.959 -          case C:
   1.960 -            //growOrAugment grows if y is covered by the matching and
   1.961 -            //augments if not. In this latter case it returns 1.
   1.962 -            if (growOrAugment(y, x, ear, blossom, rep, tree, Q)) return;
   1.963 -            break;
   1.964 -          default: break;
   1.965 -          }
   1.966 -        }
   1.967 -      }
   1.968 -    }
   1.969 -
   1.970 -    void shrink(Node x,Node y, typename Graph::template NodeMap<Node>& ear,
   1.971 -                UFE& blossom, NV& rep, EFE& tree,std::queue<Node>& Q) {
   1.972 -      //x and y are the two adjacent vertices in two blossoms.
   1.973 -
   1.974 -      typename Graph::template NodeMap<bool> path(g,false);
   1.975 -
   1.976 -      Node b=rep[blossom.find(x)];
   1.977 -      path.set(b,true);
   1.978 -      b=_mate[b];
   1.979 -      while ( b!=INVALID ) {
   1.980 -        b=rep[blossom.find(ear[b])];
   1.981 -        path.set(b,true);
   1.982 -        b=_mate[b];
   1.983 -      } //we go until the root through bases of blossoms and odd vertices
   1.984 -
   1.985 -      Node top=y;
   1.986 -      Node middle=rep[blossom.find(top)];
   1.987 -      Node bottom=x;
   1.988 -      while ( !path[middle] )
   1.989 -        shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
   1.990 -      //Until we arrive to a node on the path, we update blossom, tree
   1.991 -      //and the positions of the odd nodes.
   1.992 -
   1.993 -      Node base=middle;
   1.994 -      top=x;
   1.995 -      middle=rep[blossom.find(top)];
   1.996 -      bottom=y;
   1.997 -      Node blossom_base=rep[blossom.find(base)];
   1.998 -      while ( middle!=blossom_base )
   1.999 -        shrinkStep(top, middle, bottom, ear, blossom, rep, tree, Q);
  1.1000 -      //Until we arrive to a node on the path, we update blossom, tree
  1.1001 -      //and the positions of the odd nodes.
  1.1002 -
  1.1003 -      rep[blossom.find(base)] = base;
  1.1004 -    }
  1.1005 -
  1.1006 -    void shrinkStep(Node& top, Node& middle, Node& bottom,
  1.1007 -                    typename Graph::template NodeMap<Node>& ear,
  1.1008 -                    UFE& blossom, NV& rep, EFE& tree, std::queue<Node>& Q) {
  1.1009 -      //We traverse a blossom and update everything.
  1.1010 -
  1.1011 -      ear.set(top,bottom);
  1.1012 -      Node t=top;
  1.1013 -      while ( t!=middle ) {
  1.1014 -        Node u=_mate[t];
  1.1015 -        t=ear[u];
  1.1016 -        ear.set(t,u);
  1.1017 -      }
  1.1018 -      bottom=_mate[middle];
  1.1019 -      position.set(bottom,D);
  1.1020 -      Q.push(bottom);
  1.1021 -      top=ear[bottom];
  1.1022 -      Node oldmiddle=middle;
  1.1023 -      middle=rep[blossom.find(top)];
  1.1024 -      tree.erase(bottom);
  1.1025 -      tree.erase(oldmiddle);
  1.1026 -      blossom.insert(bottom);
  1.1027 -      blossom.join(bottom, oldmiddle);
  1.1028 -      blossom.join(top, oldmiddle);
  1.1029 -    }
  1.1030 -
  1.1031 -
  1.1032 -
  1.1033 -    bool growOrAugment(Node& y, Node& x, typename Graph::template
  1.1034 -                       NodeMap<Node>& ear, UFE& blossom, NV& rep, EFE& tree,
  1.1035 -                       std::queue<Node>& Q) {
  1.1036 -      //x is in a blossom in the tree, y is outside. If y is covered by
  1.1037 -      //the matching we grow, otherwise we augment. In this case we
  1.1038 -      //return 1.
  1.1039 -
  1.1040 -      if ( _mate[y]!=INVALID ) {       //grow
  1.1041 -        ear.set(y,x);
  1.1042 -        Node w=_mate[y];
  1.1043 -        rep[blossom.insert(w)] = w;
  1.1044 -        position.set(y,A);
  1.1045 -        position.set(w,D);
  1.1046 -        int t = tree.find(rep[blossom.find(x)]);
  1.1047 -        tree.insert(y,t);
  1.1048 -        tree.insert(w,t);
  1.1049 -        Q.push(w);
  1.1050 -      } else {                      //augment
  1.1051 -        augment(x, ear, blossom, rep, tree);
  1.1052 -        _mate.set(x,y);
  1.1053 -        _mate.set(y,x);
  1.1054 -        return true;
  1.1055 -      }
  1.1056 -      return false;
  1.1057 -    }
  1.1058 -
  1.1059 -    void augment(Node x, typename Graph::template NodeMap<Node>& ear,
  1.1060 -                 UFE& blossom, NV& rep, EFE& tree) {
  1.1061 -      Node v=_mate[x];
  1.1062 -      while ( v!=INVALID ) {
  1.1063 -
  1.1064 -        Node u=ear[v];
  1.1065 -        _mate.set(v,u);
  1.1066 -        Node tmp=v;
  1.1067 -        v=_mate[u];
  1.1068 -        _mate.set(u,tmp);
  1.1069 -      }
  1.1070 -      int y = tree.find(rep[blossom.find(x)]);
  1.1071 -      for (typename EFE::ItemIt tit(tree, y); tit != INVALID; ++tit) {
  1.1072 -        if ( position[tit] == D ) {
  1.1073 -          int b = blossom.find(tit);
  1.1074 -          for (typename UFE::ItemIt bit(blossom, b); bit != INVALID; ++bit) {
  1.1075 -            position.set(bit, C);
  1.1076 -          }
  1.1077 -          blossom.eraseClass(b);
  1.1078 -        } else position.set(tit, C);
  1.1079 -      }
  1.1080 -      tree.eraseClass(y);
  1.1081 -
  1.1082 -    }
  1.1083 +    /// @}
  1.1084  
  1.1085    };
  1.1086  
  1.1087 @@ -627,25 +617,28 @@
  1.1088    /// \f$O(nm\log(n))\f$ time complexity.
  1.1089    ///
  1.1090    /// The maximum weighted matching problem is to find undirected
  1.1091 -  /// arcs in the digraph with maximum overall weight and no two of
  1.1092 -  /// them shares their endpoints. The problem can be formulated with
  1.1093 -  /// the next linear program:
  1.1094 +  /// edges in the graph with maximum overall weight and no two of
  1.1095 +  /// them shares their ends. The problem can be formulated with the
  1.1096 +  /// following linear program.
  1.1097    /// \f[ \sum_{e \in \delta(u)}x_e \le 1 \quad \forall u\in V\f]
  1.1098 -  ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
  1.1099 +  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
  1.1100 +      \quad \forall B\in\mathcal{O}\f] */
  1.1101    /// \f[x_e \ge 0\quad \forall e\in E\f]
  1.1102    /// \f[\max \sum_{e\in E}x_ew_e\f]
  1.1103 -  /// where \f$\delta(X)\f$ is the set of arcs incident to a node in
  1.1104 -  /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in
  1.1105 -  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
  1.1106 -  /// the nodes.
  1.1107 +  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
  1.1108 +  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
  1.1109 +  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
  1.1110 +  /// subsets of the nodes.
  1.1111    ///
  1.1112    /// The algorithm calculates an optimal matching and a proof of the
  1.1113    /// optimality. The solution of the dual problem can be used to check
  1.1114 -  /// the result of the algorithm. The dual linear problem is the next:
  1.1115 -  /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f]
  1.1116 +  /// the result of the algorithm. The dual linear problem is the
  1.1117 +  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}
  1.1118 +      z_B \ge w_{uv} \quad \forall uv\in E\f] */
  1.1119    /// \f[y_u \ge 0 \quad \forall u \in V\f]
  1.1120    /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
  1.1121 -  /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
  1.1122 +  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
  1.1123 +      \frac{\vert B \vert - 1}{2}z_B\f] */
  1.1124    ///
  1.1125    /// The algorithm can be executed with \c run() or the \c init() and
  1.1126    /// then the \c start() member functions. After it the matching can
  1.1127 @@ -705,16 +698,13 @@
  1.1128      int _node_num;
  1.1129      int _blossom_num;
  1.1130  
  1.1131 -    typedef typename Graph::template NodeMap<int> NodeIntMap;
  1.1132 -    typedef typename Graph::template ArcMap<int> ArcIntMap;
  1.1133 -    typedef typename Graph::template EdgeMap<int> EdgeIntMap;
  1.1134      typedef RangeMap<int> IntIntMap;
  1.1135  
  1.1136      enum Status {
  1.1137        EVEN = -1, MATCHED = 0, ODD = 1, UNMATCHED = -2
  1.1138      };
  1.1139  
  1.1140 -    typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
  1.1141 +    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
  1.1142      struct BlossomData {
  1.1143        int tree;
  1.1144        Status status;
  1.1145 @@ -723,21 +713,21 @@
  1.1146        Node base;
  1.1147      };
  1.1148  
  1.1149 -    NodeIntMap *_blossom_index;
  1.1150 +    IntNodeMap *_blossom_index;
  1.1151      BlossomSet *_blossom_set;
  1.1152      RangeMap<BlossomData>* _blossom_data;
  1.1153  
  1.1154 -    NodeIntMap *_node_index;
  1.1155 -    ArcIntMap *_node_heap_index;
  1.1156 +    IntNodeMap *_node_index;
  1.1157 +    IntArcMap *_node_heap_index;
  1.1158  
  1.1159      struct NodeData {
  1.1160  
  1.1161 -      NodeData(ArcIntMap& node_heap_index)
  1.1162 +      NodeData(IntArcMap& node_heap_index)
  1.1163          : heap(node_heap_index) {}
  1.1164  
  1.1165        int blossom;
  1.1166        Value pot;
  1.1167 -      BinHeap<Value, ArcIntMap> heap;
  1.1168 +      BinHeap<Value, IntArcMap> heap;
  1.1169        std::map<int, Arc> heap_index;
  1.1170  
  1.1171        int tree;
  1.1172 @@ -750,14 +740,14 @@
  1.1173      IntIntMap *_tree_set_index;
  1.1174      TreeSet *_tree_set;
  1.1175  
  1.1176 -    NodeIntMap *_delta1_index;
  1.1177 -    BinHeap<Value, NodeIntMap> *_delta1;
  1.1178 +    IntNodeMap *_delta1_index;
  1.1179 +    BinHeap<Value, IntNodeMap> *_delta1;
  1.1180  
  1.1181      IntIntMap *_delta2_index;
  1.1182      BinHeap<Value, IntIntMap> *_delta2;
  1.1183  
  1.1184 -    EdgeIntMap *_delta3_index;
  1.1185 -    BinHeap<Value, EdgeIntMap> *_delta3;
  1.1186 +    IntEdgeMap *_delta3_index;
  1.1187 +    BinHeap<Value, IntEdgeMap> *_delta3;
  1.1188  
  1.1189      IntIntMap *_delta4_index;
  1.1190      BinHeap<Value, IntIntMap> *_delta4;
  1.1191 @@ -775,14 +765,14 @@
  1.1192          _node_potential = new NodePotential(_graph);
  1.1193        }
  1.1194        if (!_blossom_set) {
  1.1195 -        _blossom_index = new NodeIntMap(_graph);
  1.1196 +        _blossom_index = new IntNodeMap(_graph);
  1.1197          _blossom_set = new BlossomSet(*_blossom_index);
  1.1198          _blossom_data = new RangeMap<BlossomData>(_blossom_num);
  1.1199        }
  1.1200  
  1.1201        if (!_node_index) {
  1.1202 -        _node_index = new NodeIntMap(_graph);
  1.1203 -        _node_heap_index = new ArcIntMap(_graph);
  1.1204 +        _node_index = new IntNodeMap(_graph);
  1.1205 +        _node_heap_index = new IntArcMap(_graph);
  1.1206          _node_data = new RangeMap<NodeData>(_node_num,
  1.1207                                                NodeData(*_node_heap_index));
  1.1208        }
  1.1209 @@ -792,16 +782,16 @@
  1.1210          _tree_set = new TreeSet(*_tree_set_index);
  1.1211        }
  1.1212        if (!_delta1) {
  1.1213 -        _delta1_index = new NodeIntMap(_graph);
  1.1214 -        _delta1 = new BinHeap<Value, NodeIntMap>(*_delta1_index);
  1.1215 +        _delta1_index = new IntNodeMap(_graph);
  1.1216 +        _delta1 = new BinHeap<Value, IntNodeMap>(*_delta1_index);
  1.1217        }
  1.1218        if (!_delta2) {
  1.1219          _delta2_index = new IntIntMap(_blossom_num);
  1.1220          _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
  1.1221        }
  1.1222        if (!_delta3) {
  1.1223 -        _delta3_index = new EdgeIntMap(_graph);
  1.1224 -        _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index);
  1.1225 +        _delta3_index = new IntEdgeMap(_graph);
  1.1226 +        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
  1.1227        }
  1.1228        if (!_delta4) {
  1.1229          _delta4_index = new IntIntMap(_blossom_num);
  1.1230 @@ -1266,10 +1256,10 @@
  1.1231      }
  1.1232  
  1.1233  
  1.1234 -    void augmentOnArc(const Edge& arc) {
  1.1235 -
  1.1236 -      int left = _blossom_set->find(_graph.u(arc));
  1.1237 -      int right = _blossom_set->find(_graph.v(arc));
  1.1238 +    void augmentOnEdge(const Edge& edge) {
  1.1239 +
  1.1240 +      int left = _blossom_set->find(_graph.u(edge));
  1.1241 +      int right = _blossom_set->find(_graph.v(edge));
  1.1242  
  1.1243        if ((*_blossom_data)[left].status == EVEN) {
  1.1244          int left_tree = _tree_set->find(left);
  1.1245 @@ -1289,8 +1279,8 @@
  1.1246          unmatchedToMatched(right);
  1.1247        }
  1.1248  
  1.1249 -      (*_blossom_data)[left].next = _graph.direct(arc, true);
  1.1250 -      (*_blossom_data)[right].next = _graph.direct(arc, false);
  1.1251 +      (*_blossom_data)[left].next = _graph.direct(edge, true);
  1.1252 +      (*_blossom_data)[right].next = _graph.direct(edge, false);
  1.1253      }
  1.1254  
  1.1255      void extendOnArc(const Arc& arc) {
  1.1256 @@ -1310,7 +1300,7 @@
  1.1257        matchedToEven(even, tree);
  1.1258      }
  1.1259  
  1.1260 -    void shrinkOnArc(const Edge& edge, int tree) {
  1.1261 +    void shrinkOnEdge(const Edge& edge, int tree) {
  1.1262        int nca = -1;
  1.1263        std::vector<int> left_path, right_path;
  1.1264  
  1.1265 @@ -1652,7 +1642,7 @@
  1.1266        createStructures();
  1.1267  
  1.1268        for (ArcIt e(_graph); e != INVALID; ++e) {
  1.1269 -        _node_heap_index->set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP);
  1.1270 +        _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
  1.1271        }
  1.1272        for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1273          _delta1_index->set(n, _delta1->PRE_HEAP);
  1.1274 @@ -1769,9 +1759,9 @@
  1.1275                }
  1.1276  
  1.1277                if (left_tree == right_tree) {
  1.1278 -                shrinkOnArc(e, left_tree);
  1.1279 +                shrinkOnEdge(e, left_tree);
  1.1280                } else {
  1.1281 -                augmentOnArc(e);
  1.1282 +                augmentOnEdge(e);
  1.1283                  unmatched -= 2;
  1.1284                }
  1.1285              }
  1.1286 @@ -1818,11 +1808,24 @@
  1.1287        return sum /= 2;
  1.1288      }
  1.1289  
  1.1290 -    /// \brief Returns true when the arc is in the matching.
  1.1291 +    /// \brief Returns the cardinality of the matching.
  1.1292      ///
  1.1293 -    /// Returns true when the arc is in the matching.
  1.1294 -    bool matching(const Edge& arc) const {
  1.1295 -      return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true);
  1.1296 +    /// Returns the cardinality of the matching.
  1.1297 +    int matchingSize() const {
  1.1298 +      int num = 0;
  1.1299 +      for (NodeIt n(_graph); n != INVALID; ++n) {
  1.1300 +        if ((*_matching)[n] != INVALID) {
  1.1301 +          ++num;
  1.1302 +        }
  1.1303 +      }
  1.1304 +      return num /= 2;
  1.1305 +    }
  1.1306 +
  1.1307 +    /// \brief Returns true when the edge is in the matching.
  1.1308 +    ///
  1.1309 +    /// Returns true when the edge is in the matching.
  1.1310 +    bool matching(const Edge& edge) const {
  1.1311 +      return edge == (*_matching)[_graph.u(edge)];
  1.1312      }
  1.1313  
  1.1314      /// \brief Returns the incident matching arc.
  1.1315 @@ -1913,16 +1916,11 @@
  1.1316          _last = _algorithm->_blossom_potential[variable].end;
  1.1317        }
  1.1318  
  1.1319 -      /// \brief Invalid constructor.
  1.1320 -      ///
  1.1321 -      /// Invalid constructor.
  1.1322 -      BlossomIt(Invalid) : _index(-1) {}
  1.1323 -
  1.1324        /// \brief Conversion to node.
  1.1325        ///
  1.1326        /// Conversion to node.
  1.1327        operator Node() const {
  1.1328 -        return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
  1.1329 +        return _algorithm->_blossom_node_list[_index];
  1.1330        }
  1.1331  
  1.1332        /// \brief Increment operator.
  1.1333 @@ -1930,18 +1928,18 @@
  1.1334        /// Increment operator.
  1.1335        BlossomIt& operator++() {
  1.1336          ++_index;
  1.1337 -        if (_index == _last) {
  1.1338 -          _index = -1;
  1.1339 -        }
  1.1340          return *this;
  1.1341        }
  1.1342  
  1.1343 -      bool operator==(const BlossomIt& it) const {
  1.1344 -        return _index == it._index;
  1.1345 -      }
  1.1346 -      bool operator!=(const BlossomIt& it) const {
  1.1347 -        return _index != it._index;
  1.1348 -      }
  1.1349 +      /// \brief Validity checking
  1.1350 +      ///
  1.1351 +      /// Checks whether the iterator is invalid.
  1.1352 +      bool operator==(Invalid) const { return _index == _last; }
  1.1353 +
  1.1354 +      /// \brief Validity checking
  1.1355 +      ///
  1.1356 +      /// Checks whether the iterator is valid.
  1.1357 +      bool operator!=(Invalid) const { return _index != _last; }
  1.1358  
  1.1359      private:
  1.1360        const MaxWeightedMatching* _algorithm;
  1.1361 @@ -1958,29 +1956,32 @@
  1.1362    /// \brief Weighted perfect matching in general graphs
  1.1363    ///
  1.1364    /// This class provides an efficient implementation of Edmond's
  1.1365 -  /// maximum weighted perfecr matching algorithm. The implementation
  1.1366 +  /// maximum weighted perfect matching algorithm. The implementation
  1.1367    /// is based on extensive use of priority queues and provides
  1.1368    /// \f$O(nm\log(n))\f$ time complexity.
  1.1369    ///
  1.1370    /// The maximum weighted matching problem is to find undirected
  1.1371 -  /// arcs in the digraph with maximum overall weight and no two of
  1.1372 -  /// them shares their endpoints and covers all nodes. The problem
  1.1373 -  /// can be formulated with the next linear program:
  1.1374 +  /// edges in the graph with maximum overall weight and no two of
  1.1375 +  /// them shares their ends and covers all nodes. The problem can be
  1.1376 +  /// formulated with the following linear program.
  1.1377    /// \f[ \sum_{e \in \delta(u)}x_e = 1 \quad \forall u\in V\f]
  1.1378 -  ///\f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2} \quad \forall B\in\mathcal{O}\f]
  1.1379 +  /** \f[ \sum_{e \in \gamma(B)}x_e \le \frac{\vert B \vert - 1}{2}
  1.1380 +      \quad \forall B\in\mathcal{O}\f] */
  1.1381    /// \f[x_e \ge 0\quad \forall e\in E\f]
  1.1382    /// \f[\max \sum_{e\in E}x_ew_e\f]
  1.1383 -  /// where \f$\delta(X)\f$ is the set of arcs incident to a node in
  1.1384 -  /// \f$X\f$, \f$\gamma(X)\f$ is the set of arcs with both endpoints in
  1.1385 -  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality subsets of
  1.1386 -  /// the nodes.
  1.1387 +  /// where \f$\delta(X)\f$ is the set of edges incident to a node in
  1.1388 +  /// \f$X\f$, \f$\gamma(X)\f$ is the set of edges with both ends in
  1.1389 +  /// \f$X\f$ and \f$\mathcal{O}\f$ is the set of odd cardinality
  1.1390 +  /// subsets of the nodes.
  1.1391    ///
  1.1392    /// The algorithm calculates an optimal matching and a proof of the
  1.1393    /// optimality. The solution of the dual problem can be used to check
  1.1394 -  /// the result of the algorithm. The dual linear problem is the next:
  1.1395 -  /// \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge w_{uv} \quad \forall uv\in E\f]
  1.1396 +  /// the result of the algorithm. The dual linear problem is the
  1.1397 +  /** \f[ y_u + y_v + \sum_{B \in \mathcal{O}, uv \in \gamma(B)}z_B \ge
  1.1398 +      w_{uv} \quad \forall uv\in E\f] */
  1.1399    /// \f[z_B \ge 0 \quad \forall B \in \mathcal{O}\f]
  1.1400 -  /// \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}\frac{\vert B \vert - 1}{2}z_B\f]
  1.1401 +  /** \f[\min \sum_{u \in V}y_u + \sum_{B \in \mathcal{O}}
  1.1402 +      \frac{\vert B \vert - 1}{2}z_B\f] */
  1.1403    ///
  1.1404    /// The algorithm can be executed with \c run() or the \c init() and
  1.1405    /// then the \c start() member functions. After it the matching can
  1.1406 @@ -2040,16 +2041,13 @@
  1.1407      int _node_num;
  1.1408      int _blossom_num;
  1.1409  
  1.1410 -    typedef typename Graph::template NodeMap<int> NodeIntMap;
  1.1411 -    typedef typename Graph::template ArcMap<int> ArcIntMap;
  1.1412 -    typedef typename Graph::template EdgeMap<int> EdgeIntMap;
  1.1413      typedef RangeMap<int> IntIntMap;
  1.1414  
  1.1415      enum Status {
  1.1416        EVEN = -1, MATCHED = 0, ODD = 1
  1.1417      };
  1.1418  
  1.1419 -    typedef HeapUnionFind<Value, NodeIntMap> BlossomSet;
  1.1420 +    typedef HeapUnionFind<Value, IntNodeMap> BlossomSet;
  1.1421      struct BlossomData {
  1.1422        int tree;
  1.1423        Status status;
  1.1424 @@ -2057,21 +2055,21 @@
  1.1425        Value pot, offset;
  1.1426      };
  1.1427  
  1.1428 -    NodeIntMap *_blossom_index;
  1.1429 +    IntNodeMap *_blossom_index;
  1.1430      BlossomSet *_blossom_set;
  1.1431      RangeMap<BlossomData>* _blossom_data;
  1.1432  
  1.1433 -    NodeIntMap *_node_index;
  1.1434 -    ArcIntMap *_node_heap_index;
  1.1435 +    IntNodeMap *_node_index;
  1.1436 +    IntArcMap *_node_heap_index;
  1.1437  
  1.1438      struct NodeData {
  1.1439  
  1.1440 -      NodeData(ArcIntMap& node_heap_index)
  1.1441 +      NodeData(IntArcMap& node_heap_index)
  1.1442          : heap(node_heap_index) {}
  1.1443  
  1.1444        int blossom;
  1.1445        Value pot;
  1.1446 -      BinHeap<Value, ArcIntMap> heap;
  1.1447 +      BinHeap<Value, IntArcMap> heap;
  1.1448        std::map<int, Arc> heap_index;
  1.1449  
  1.1450        int tree;
  1.1451 @@ -2087,8 +2085,8 @@
  1.1452      IntIntMap *_delta2_index;
  1.1453      BinHeap<Value, IntIntMap> *_delta2;
  1.1454  
  1.1455 -    EdgeIntMap *_delta3_index;
  1.1456 -    BinHeap<Value, EdgeIntMap> *_delta3;
  1.1457 +    IntEdgeMap *_delta3_index;
  1.1458 +    BinHeap<Value, IntEdgeMap> *_delta3;
  1.1459  
  1.1460      IntIntMap *_delta4_index;
  1.1461      BinHeap<Value, IntIntMap> *_delta4;
  1.1462 @@ -2106,16 +2104,16 @@
  1.1463          _node_potential = new NodePotential(_graph);
  1.1464        }
  1.1465        if (!_blossom_set) {
  1.1466 -        _blossom_index = new NodeIntMap(_graph);
  1.1467 +        _blossom_index = new IntNodeMap(_graph);
  1.1468          _blossom_set = new BlossomSet(*_blossom_index);
  1.1469          _blossom_data = new RangeMap<BlossomData>(_blossom_num);
  1.1470        }
  1.1471  
  1.1472        if (!_node_index) {
  1.1473 -        _node_index = new NodeIntMap(_graph);
  1.1474 -        _node_heap_index = new ArcIntMap(_graph);
  1.1475 +        _node_index = new IntNodeMap(_graph);
  1.1476 +        _node_heap_index = new IntArcMap(_graph);
  1.1477          _node_data = new RangeMap<NodeData>(_node_num,
  1.1478 -                                              NodeData(*_node_heap_index));
  1.1479 +                                            NodeData(*_node_heap_index));
  1.1480        }
  1.1481  
  1.1482        if (!_tree_set) {
  1.1483 @@ -2127,8 +2125,8 @@
  1.1484          _delta2 = new BinHeap<Value, IntIntMap>(*_delta2_index);
  1.1485        }
  1.1486        if (!_delta3) {
  1.1487 -        _delta3_index = new EdgeIntMap(_graph);
  1.1488 -        _delta3 = new BinHeap<Value, EdgeIntMap>(*_delta3_index);
  1.1489 +        _delta3_index = new IntEdgeMap(_graph);
  1.1490 +        _delta3 = new BinHeap<Value, IntEdgeMap>(*_delta3_index);
  1.1491        }
  1.1492        if (!_delta4) {
  1.1493          _delta4_index = new IntIntMap(_blossom_num);
  1.1494 @@ -2461,10 +2459,10 @@
  1.1495        _tree_set->eraseClass(tree);
  1.1496      }
  1.1497  
  1.1498 -    void augmentOnArc(const Edge& arc) {
  1.1499 -
  1.1500 -      int left = _blossom_set->find(_graph.u(arc));
  1.1501 -      int right = _blossom_set->find(_graph.v(arc));
  1.1502 +    void augmentOnEdge(const Edge& edge) {
  1.1503 +
  1.1504 +      int left = _blossom_set->find(_graph.u(edge));
  1.1505 +      int right = _blossom_set->find(_graph.v(edge));
  1.1506  
  1.1507        int left_tree = _tree_set->find(left);
  1.1508        alternatePath(left, left_tree);
  1.1509 @@ -2474,8 +2472,8 @@
  1.1510        alternatePath(right, right_tree);
  1.1511        destroyTree(right_tree);
  1.1512  
  1.1513 -      (*_blossom_data)[left].next = _graph.direct(arc, true);
  1.1514 -      (*_blossom_data)[right].next = _graph.direct(arc, false);
  1.1515 +      (*_blossom_data)[left].next = _graph.direct(edge, true);
  1.1516 +      (*_blossom_data)[right].next = _graph.direct(edge, false);
  1.1517      }
  1.1518  
  1.1519      void extendOnArc(const Arc& arc) {
  1.1520 @@ -2495,7 +2493,7 @@
  1.1521        matchedToEven(even, tree);
  1.1522      }
  1.1523  
  1.1524 -    void shrinkOnArc(const Edge& edge, int tree) {
  1.1525 +    void shrinkOnEdge(const Edge& edge, int tree) {
  1.1526        int nca = -1;
  1.1527        std::vector<int> left_path, right_path;
  1.1528  
  1.1529 @@ -2831,7 +2829,7 @@
  1.1530        createStructures();
  1.1531  
  1.1532        for (ArcIt e(_graph); e != INVALID; ++e) {
  1.1533 -        _node_heap_index->set(e, BinHeap<Value, ArcIntMap>::PRE_HEAP);
  1.1534 +        _node_heap_index->set(e, BinHeap<Value, IntArcMap>::PRE_HEAP);
  1.1535        }
  1.1536        for (EdgeIt e(_graph); e != INVALID; ++e) {
  1.1537          _delta3_index->set(e, _delta3->PRE_HEAP);
  1.1538 @@ -2924,9 +2922,9 @@
  1.1539                int right_tree = _tree_set->find(right_blossom);
  1.1540  
  1.1541                if (left_tree == right_tree) {
  1.1542 -                shrinkOnArc(e, left_tree);
  1.1543 +                shrinkOnEdge(e, left_tree);
  1.1544                } else {
  1.1545 -                augmentOnArc(e);
  1.1546 +                augmentOnEdge(e);
  1.1547                  unmatched -= 2;
  1.1548                }
  1.1549              }
  1.1550 @@ -2974,16 +2972,16 @@
  1.1551        return sum /= 2;
  1.1552      }
  1.1553  
  1.1554 -    /// \brief Returns true when the arc is in the matching.
  1.1555 +    /// \brief Returns true when the edge is in the matching.
  1.1556      ///
  1.1557 -    /// Returns true when the arc is in the matching.
  1.1558 -    bool matching(const Edge& arc) const {
  1.1559 -      return (*_matching)[_graph.u(arc)] == _graph.direct(arc, true);
  1.1560 +    /// Returns true when the edge is in the matching.
  1.1561 +    bool matching(const Edge& edge) const {
  1.1562 +      return static_cast<const Edge&>((*_matching)[_graph.u(edge)]) == edge;
  1.1563      }
  1.1564  
  1.1565 -    /// \brief Returns the incident matching arc.
  1.1566 +    /// \brief Returns the incident matching edge.
  1.1567      ///
  1.1568 -    /// Returns the incident matching arc from given node.
  1.1569 +    /// Returns the incident matching arc from given edge.
  1.1570      Arc matching(const Node& node) const {
  1.1571        return (*_matching)[node];
  1.1572      }
  1.1573 @@ -3066,16 +3064,11 @@
  1.1574          _last = _algorithm->_blossom_potential[variable].end;
  1.1575        }
  1.1576  
  1.1577 -      /// \brief Invalid constructor.
  1.1578 -      ///
  1.1579 -      /// Invalid constructor.
  1.1580 -      BlossomIt(Invalid) : _index(-1) {}
  1.1581 -
  1.1582        /// \brief Conversion to node.
  1.1583        ///
  1.1584        /// Conversion to node.
  1.1585        operator Node() const {
  1.1586 -        return _algorithm ? _algorithm->_blossom_node_list[_index] : INVALID;
  1.1587 +        return _algorithm->_blossom_node_list[_index];
  1.1588        }
  1.1589  
  1.1590        /// \brief Increment operator.
  1.1591 @@ -3083,18 +3076,18 @@
  1.1592        /// Increment operator.
  1.1593        BlossomIt& operator++() {
  1.1594          ++_index;
  1.1595 -        if (_index == _last) {
  1.1596 -          _index = -1;
  1.1597 -        }
  1.1598          return *this;
  1.1599        }
  1.1600  
  1.1601 -      bool operator==(const BlossomIt& it) const {
  1.1602 -        return _index == it._index;
  1.1603 -      }
  1.1604 -      bool operator!=(const BlossomIt& it) const {
  1.1605 -        return _index != it._index;
  1.1606 -      }
  1.1607 +      /// \brief Validity checking
  1.1608 +      ///
  1.1609 +      /// Checks whether the iterator is invalid.
  1.1610 +      bool operator==(Invalid) const { return _index == _last; }
  1.1611 +
  1.1612 +      /// \brief Validity checking
  1.1613 +      ///
  1.1614 +      /// Checks whether the iterator is valid.
  1.1615 +      bool operator!=(Invalid) const { return _index != _last; }
  1.1616  
  1.1617      private:
  1.1618        const MaxWeightedPerfectMatching* _algorithm;
     2.1 --- a/test/CMakeLists.txt	Mon Oct 13 13:56:00 2008 +0200
     2.2 +++ b/test/CMakeLists.txt	Mon Oct 13 14:00:11 2008 +0200
     2.3 @@ -17,7 +17,6 @@
     2.4    kruskal_test
     2.5    maps_test
     2.6    max_matching_test
     2.7 -  max_weighted_matching_test
     2.8    random_test
     2.9    path_test
    2.10    time_measure_test
     3.1 --- a/test/Makefile.am	Mon Oct 13 13:56:00 2008 +0200
     3.2 +++ b/test/Makefile.am	Mon Oct 13 14:00:11 2008 +0200
     3.3 @@ -20,7 +20,6 @@
     3.4  	test/kruskal_test \
     3.5          test/maps_test \
     3.6  	test/max_matching_test \
     3.7 -	test/max_weighted_matching_test \
     3.8          test/random_test \
     3.9          test/path_test \
    3.10          test/test_tools_fail \
    3.11 @@ -45,7 +44,6 @@
    3.12  test_kruskal_test_SOURCES = test/kruskal_test.cc
    3.13  test_maps_test_SOURCES = test/maps_test.cc
    3.14  test_max_matching_test_SOURCES = test/max_matching_test.cc
    3.15 -test_max_weighted_matching_test_SOURCES = test/max_weighted_matching_test.cc
    3.16  test_path_test_SOURCES = test/path_test.cc
    3.17  test_random_test_SOURCES = test/random_test.cc
    3.18  test_test_tools_fail_SOURCES = test/test_tools_fail.cc
     4.1 --- a/test/max_matching_test.cc	Mon Oct 13 13:56:00 2008 +0200
     4.2 +++ b/test/max_matching_test.cc	Mon Oct 13 14:00:11 2008 +0200
     4.3 @@ -17,165 +17,294 @@
     4.4   */
     4.5  
     4.6  #include <iostream>
     4.7 +#include <sstream>
     4.8  #include <vector>
     4.9  #include <queue>
    4.10  #include <lemon/math.h>
    4.11  #include <cstdlib>
    4.12  
    4.13 +#include <lemon/max_matching.h>
    4.14 +#include <lemon/smart_graph.h>
    4.15 +#include <lemon/lgf_reader.h>
    4.16 +
    4.17  #include "test_tools.h"
    4.18 -#include <lemon/list_graph.h>
    4.19 -#include <lemon/max_matching.h>
    4.20  
    4.21  using namespace std;
    4.22  using namespace lemon;
    4.23  
    4.24 +GRAPH_TYPEDEFS(SmartGraph);
    4.25 +
    4.26 +
    4.27 +const int lgfn = 3;
    4.28 +const std::string lgf[lgfn] = {
    4.29 +  "@nodes\n"
    4.30 +  "label\n"
    4.31 +  "0\n"
    4.32 +  "1\n"
    4.33 +  "2\n"
    4.34 +  "3\n"
    4.35 +  "4\n"
    4.36 +  "5\n"
    4.37 +  "6\n"
    4.38 +  "7\n"
    4.39 +  "@edges\n"
    4.40 +  "     label  weight\n"
    4.41 +  "7 4  0      984\n"
    4.42 +  "0 7  1      73\n"
    4.43 +  "7 1  2      204\n"
    4.44 +  "2 3  3      583\n"
    4.45 +  "2 7  4      565\n"
    4.46 +  "2 1  5      582\n"
    4.47 +  "0 4  6      551\n"
    4.48 +  "2 5  7      385\n"
    4.49 +  "1 5  8      561\n"
    4.50 +  "5 3  9      484\n"
    4.51 +  "7 5  10     904\n"
    4.52 +  "3 6  11     47\n"
    4.53 +  "7 6  12     888\n"
    4.54 +  "3 0  13     747\n"
    4.55 +  "6 1  14     310\n",
    4.56 +
    4.57 +  "@nodes\n"
    4.58 +  "label\n"
    4.59 +  "0\n"
    4.60 +  "1\n"
    4.61 +  "2\n"
    4.62 +  "3\n"
    4.63 +  "4\n"
    4.64 +  "5\n"
    4.65 +  "6\n"
    4.66 +  "7\n"
    4.67 +  "@edges\n"
    4.68 +  "     label  weight\n"
    4.69 +  "2 5  0      710\n"
    4.70 +  "0 5  1      241\n"
    4.71 +  "2 4  2      856\n"
    4.72 +  "2 6  3      762\n"
    4.73 +  "4 1  4      747\n"
    4.74 +  "6 1  5      962\n"
    4.75 +  "4 7  6      723\n"
    4.76 +  "1 7  7      661\n"
    4.77 +  "2 3  8      376\n"
    4.78 +  "1 0  9      416\n"
    4.79 +  "6 7  10     391\n",
    4.80 +
    4.81 +  "@nodes\n"
    4.82 +  "label\n"
    4.83 +  "0\n"
    4.84 +  "1\n"
    4.85 +  "2\n"
    4.86 +  "3\n"
    4.87 +  "4\n"
    4.88 +  "5\n"
    4.89 +  "6\n"
    4.90 +  "7\n"
    4.91 +  "@edges\n"
    4.92 +  "     label  weight\n"
    4.93 +  "6 2  0      553\n"
    4.94 +  "0 7  1      653\n"
    4.95 +  "6 3  2      22\n"
    4.96 +  "4 7  3      846\n"
    4.97 +  "7 2  4      981\n"
    4.98 +  "7 6  5      250\n"
    4.99 +  "5 2  6      539\n",
   4.100 +};
   4.101 +
   4.102 +void checkMatching(const SmartGraph& graph,
   4.103 +                   const MaxMatching<SmartGraph>& mm) {
   4.104 +  int num = 0;
   4.105 +
   4.106 +  IntNodeMap comp_index(graph);
   4.107 +  UnionFind<IntNodeMap> comp(comp_index);
   4.108 +
   4.109 +  int barrier_num = 0;
   4.110 +
   4.111 +  for (NodeIt n(graph); n != INVALID; ++n) {
   4.112 +    check(mm.decomposition(n) == MaxMatching<SmartGraph>::EVEN ||
   4.113 +          mm.matching(n) != INVALID, "Wrong Gallai-Edmonds decomposition");
   4.114 +    if (mm.decomposition(n) == MaxMatching<SmartGraph>::ODD) {
   4.115 +      ++barrier_num;
   4.116 +    } else {
   4.117 +      comp.insert(n);
   4.118 +    }
   4.119 +  }
   4.120 +
   4.121 +  for (EdgeIt e(graph); e != INVALID; ++e) {
   4.122 +    if (mm.matching(e)) {
   4.123 +      check(e == mm.matching(graph.u(e)), "Wrong matching");
   4.124 +      check(e == mm.matching(graph.v(e)), "Wrong matching");
   4.125 +      ++num;
   4.126 +    }
   4.127 +    check(mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::EVEN ||
   4.128 +          mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::MATCHED,
   4.129 +          "Wrong Gallai-Edmonds decomposition");
   4.130 +
   4.131 +    check(mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::EVEN ||
   4.132 +          mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::MATCHED,
   4.133 +          "Wrong Gallai-Edmonds decomposition");
   4.134 +
   4.135 +    if (mm.decomposition(graph.u(e)) != MaxMatching<SmartGraph>::ODD &&
   4.136 +        mm.decomposition(graph.v(e)) != MaxMatching<SmartGraph>::ODD) {
   4.137 +      comp.join(graph.u(e), graph.v(e));
   4.138 +    }
   4.139 +  }
   4.140 +
   4.141 +  std::set<int> comp_root;
   4.142 +  int odd_comp_num = 0;
   4.143 +  for (NodeIt n(graph); n != INVALID; ++n) {
   4.144 +    if (mm.decomposition(n) != MaxMatching<SmartGraph>::ODD) {
   4.145 +      int root = comp.find(n);
   4.146 +      if (comp_root.find(root) == comp_root.end()) {
   4.147 +        comp_root.insert(root);
   4.148 +        if (comp.size(n) % 2 == 1) {
   4.149 +          ++odd_comp_num;
   4.150 +        }
   4.151 +      }
   4.152 +    }
   4.153 +  }
   4.154 +
   4.155 +  check(mm.matchingSize() == num, "Wrong matching");
   4.156 +  check(2 * num == countNodes(graph) - (odd_comp_num - barrier_num),
   4.157 +         "Wrong matching");
   4.158 +  return;
   4.159 +}
   4.160 +
   4.161 +void checkWeightedMatching(const SmartGraph& graph,
   4.162 +                   const SmartGraph::EdgeMap<int>& weight,
   4.163 +                   const MaxWeightedMatching<SmartGraph>& mwm) {
   4.164 +  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
   4.165 +    if (graph.u(e) == graph.v(e)) continue;
   4.166 +    int rw = mwm.nodeValue(graph.u(e)) + mwm.nodeValue(graph.v(e));
   4.167 +
   4.168 +    for (int i = 0; i < mwm.blossomNum(); ++i) {
   4.169 +      bool s = false, t = false;
   4.170 +      for (MaxWeightedMatching<SmartGraph>::BlossomIt n(mwm, i);
   4.171 +           n != INVALID; ++n) {
   4.172 +        if (graph.u(e) == n) s = true;
   4.173 +        if (graph.v(e) == n) t = true;
   4.174 +      }
   4.175 +      if (s == true && t == true) {
   4.176 +        rw += mwm.blossomValue(i);
   4.177 +      }
   4.178 +    }
   4.179 +    rw -= weight[e] * mwm.dualScale;
   4.180 +
   4.181 +    check(rw >= 0, "Negative reduced weight");
   4.182 +    check(rw == 0 || !mwm.matching(e),
   4.183 +          "Non-zero reduced weight on matching edge");
   4.184 +  }
   4.185 +
   4.186 +  int pv = 0;
   4.187 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   4.188 +    if (mwm.matching(n) != INVALID) {
   4.189 +      check(mwm.nodeValue(n) >= 0, "Invalid node value");
   4.190 +      pv += weight[mwm.matching(n)];
   4.191 +      SmartGraph::Node o = graph.target(mwm.matching(n));
   4.192 +      check(mwm.mate(n) == o, "Invalid matching");
   4.193 +      check(mwm.matching(n) == graph.oppositeArc(mwm.matching(o)),
   4.194 +            "Invalid matching");
   4.195 +    } else {
   4.196 +      check(mwm.mate(n) == INVALID, "Invalid matching");
   4.197 +      check(mwm.nodeValue(n) == 0, "Invalid matching");
   4.198 +    }
   4.199 +  }
   4.200 +
   4.201 +  int dv = 0;
   4.202 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   4.203 +    dv += mwm.nodeValue(n);
   4.204 +  }
   4.205 +
   4.206 +  for (int i = 0; i < mwm.blossomNum(); ++i) {
   4.207 +    check(mwm.blossomValue(i) >= 0, "Invalid blossom value");
   4.208 +    check(mwm.blossomSize(i) % 2 == 1, "Even blossom size");
   4.209 +    dv += mwm.blossomValue(i) * ((mwm.blossomSize(i) - 1) / 2);
   4.210 +  }
   4.211 +
   4.212 +  check(pv * mwm.dualScale == dv * 2, "Wrong duality");
   4.213 +
   4.214 +  return;
   4.215 +}
   4.216 +
   4.217 +void checkWeightedPerfectMatching(const SmartGraph& graph,
   4.218 +                          const SmartGraph::EdgeMap<int>& weight,
   4.219 +                          const MaxWeightedPerfectMatching<SmartGraph>& mwpm) {
   4.220 +  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
   4.221 +    if (graph.u(e) == graph.v(e)) continue;
   4.222 +    int rw = mwpm.nodeValue(graph.u(e)) + mwpm.nodeValue(graph.v(e));
   4.223 +
   4.224 +    for (int i = 0; i < mwpm.blossomNum(); ++i) {
   4.225 +      bool s = false, t = false;
   4.226 +      for (MaxWeightedPerfectMatching<SmartGraph>::BlossomIt n(mwpm, i);
   4.227 +           n != INVALID; ++n) {
   4.228 +        if (graph.u(e) == n) s = true;
   4.229 +        if (graph.v(e) == n) t = true;
   4.230 +      }
   4.231 +      if (s == true && t == true) {
   4.232 +        rw += mwpm.blossomValue(i);
   4.233 +      }
   4.234 +    }
   4.235 +    rw -= weight[e] * mwpm.dualScale;
   4.236 +
   4.237 +    check(rw >= 0, "Negative reduced weight");
   4.238 +    check(rw == 0 || !mwpm.matching(e),
   4.239 +          "Non-zero reduced weight on matching edge");
   4.240 +  }
   4.241 +
   4.242 +  int pv = 0;
   4.243 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   4.244 +    check(mwpm.matching(n) != INVALID, "Non perfect");
   4.245 +    pv += weight[mwpm.matching(n)];
   4.246 +    SmartGraph::Node o = graph.target(mwpm.matching(n));
   4.247 +    check(mwpm.mate(n) == o, "Invalid matching");
   4.248 +    check(mwpm.matching(n) == graph.oppositeArc(mwpm.matching(o)),
   4.249 +          "Invalid matching");
   4.250 +  }
   4.251 +
   4.252 +  int dv = 0;
   4.253 +  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   4.254 +    dv += mwpm.nodeValue(n);
   4.255 +  }
   4.256 +
   4.257 +  for (int i = 0; i < mwpm.blossomNum(); ++i) {
   4.258 +    check(mwpm.blossomValue(i) >= 0, "Invalid blossom value");
   4.259 +    check(mwpm.blossomSize(i) % 2 == 1, "Even blossom size");
   4.260 +    dv += mwpm.blossomValue(i) * ((mwpm.blossomSize(i) - 1) / 2);
   4.261 +  }
   4.262 +
   4.263 +  check(pv * mwpm.dualScale == dv * 2, "Wrong duality");
   4.264 +
   4.265 +  return;
   4.266 +}
   4.267 +
   4.268 +
   4.269  int main() {
   4.270  
   4.271 -  typedef ListGraph Graph;
   4.272 +  for (int i = 0; i < lgfn; ++i) {
   4.273 +    SmartGraph graph;
   4.274 +    SmartGraph::EdgeMap<int> weight(graph);
   4.275  
   4.276 -  GRAPH_TYPEDEFS(Graph);
   4.277 +    istringstream lgfs(lgf[i]);
   4.278 +    graphReader(graph, lgfs).
   4.279 +      edgeMap("weight", weight).run();
   4.280  
   4.281 -  Graph g;
   4.282 -  g.clear();
   4.283 +    MaxMatching<SmartGraph> mm(graph);
   4.284 +    mm.run();
   4.285 +    checkMatching(graph, mm);
   4.286  
   4.287 -  std::vector<Graph::Node> nodes;
   4.288 -  for (int i=0; i<13; ++i)
   4.289 -      nodes.push_back(g.addNode());
   4.290 +    MaxWeightedMatching<SmartGraph> mwm(graph, weight);
   4.291 +    mwm.run();
   4.292 +    checkWeightedMatching(graph, weight, mwm);
   4.293  
   4.294 -  g.addEdge(nodes[0], nodes[0]);
   4.295 -  g.addEdge(nodes[6], nodes[10]);
   4.296 -  g.addEdge(nodes[5], nodes[10]);
   4.297 -  g.addEdge(nodes[4], nodes[10]);
   4.298 -  g.addEdge(nodes[3], nodes[11]);
   4.299 -  g.addEdge(nodes[1], nodes[6]);
   4.300 -  g.addEdge(nodes[4], nodes[7]);
   4.301 -  g.addEdge(nodes[1], nodes[8]);
   4.302 -  g.addEdge(nodes[0], nodes[8]);
   4.303 -  g.addEdge(nodes[3], nodes[12]);
   4.304 -  g.addEdge(nodes[6], nodes[9]);
   4.305 -  g.addEdge(nodes[9], nodes[11]);
   4.306 -  g.addEdge(nodes[2], nodes[10]);
   4.307 -  g.addEdge(nodes[10], nodes[8]);
   4.308 -  g.addEdge(nodes[5], nodes[8]);
   4.309 -  g.addEdge(nodes[6], nodes[3]);
   4.310 -  g.addEdge(nodes[0], nodes[5]);
   4.311 -  g.addEdge(nodes[6], nodes[12]);
   4.312 +    MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
   4.313 +    bool perfect = mwpm.run();
   4.314  
   4.315 -  MaxMatching<Graph> max_matching(g);
   4.316 -  max_matching.init();
   4.317 -  max_matching.startDense();
   4.318 +    check(perfect == (mm.matchingSize() * 2 == countNodes(graph)),
   4.319 +          "Perfect matching found");
   4.320  
   4.321 -  int s=0;
   4.322 -  Graph::NodeMap<Node> mate(g,INVALID);
   4.323 -  max_matching.mateMap(mate);
   4.324 -  for(NodeIt v(g); v!=INVALID; ++v) {
   4.325 -    if ( mate[v]!=INVALID ) ++s;
   4.326 -  }
   4.327 -  int size=int(s/2);  //size will be used as the size of a maxmatching
   4.328 -
   4.329 -  for(NodeIt v(g); v!=INVALID; ++v) {
   4.330 -    max_matching.mate(v);
   4.331 -  }
   4.332 -
   4.333 -  check ( size == max_matching.size(), "mate() returns a different size matching than max_matching.size()" );
   4.334 -
   4.335 -  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos0(g);
   4.336 -  max_matching.decomposition(pos0);
   4.337 -
   4.338 -  max_matching.init();
   4.339 -  max_matching.startSparse();
   4.340 -  s=0;
   4.341 -  max_matching.mateMap(mate);
   4.342 -  for(NodeIt v(g); v!=INVALID; ++v) {
   4.343 -    if ( mate[v]!=INVALID ) ++s;
   4.344 -  }
   4.345 -  check ( int(s/2) == size, "The size does not equal!" );
   4.346 -
   4.347 -  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos1(g);
   4.348 -  max_matching.decomposition(pos1);
   4.349 -
   4.350 -  max_matching.run();
   4.351 -  s=0;
   4.352 -  max_matching.mateMap(mate);
   4.353 -  for(NodeIt v(g); v!=INVALID; ++v) {
   4.354 -    if ( mate[v]!=INVALID ) ++s;
   4.355 -  }
   4.356 -  check ( int(s/2) == size, "The size does not equal!" );
   4.357 -
   4.358 -  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos2(g);
   4.359 -  max_matching.decomposition(pos2);
   4.360 -
   4.361 -  max_matching.run();
   4.362 -  s=0;
   4.363 -  max_matching.mateMap(mate);
   4.364 -  for(NodeIt v(g); v!=INVALID; ++v) {
   4.365 -    if ( mate[v]!=INVALID ) ++s;
   4.366 -  }
   4.367 -  check ( int(s/2) == size, "The size does not equal!" );
   4.368 -
   4.369 -  Graph::NodeMap<MaxMatching<Graph>::DecompType> pos(g);
   4.370 -  max_matching.decomposition(pos);
   4.371 -
   4.372 -  bool ismatching=true;
   4.373 -  for(NodeIt v(g); v!=INVALID; ++v) {
   4.374 -    if ( mate[v]!=INVALID ) {
   4.375 -      Node u=mate[v];
   4.376 -      if (mate[u]!=v) ismatching=false;
   4.377 +    if (perfect) {
   4.378 +      checkWeightedPerfectMatching(graph, weight, mwpm);
   4.379      }
   4.380    }
   4.381 -  check ( ismatching, "It is not a matching!" );
   4.382 -
   4.383 -  bool coincide=true;
   4.384 -  for(NodeIt v(g); v!=INVALID; ++v) {
   4.385 -   if ( pos0[v] != pos1[v] || pos1[v]!=pos2[v] || pos2[v]!=pos[v] ) {
   4.386 -     coincide=false;
   4.387 -    }
   4.388 -  }
   4.389 -  check ( coincide, "The decompositions do not coincide! " );
   4.390 -
   4.391 -  bool noarc=true;
   4.392 -  for(EdgeIt e(g); e!=INVALID; ++e) {
   4.393 -   if ( (pos[g.v(e)]==max_matching.C &&
   4.394 -         pos[g.u(e)]==max_matching.D) ||
   4.395 -         (pos[g.v(e)]==max_matching.D &&
   4.396 -          pos[g.u(e)]==max_matching.C) )
   4.397 -      noarc=false;
   4.398 -  }
   4.399 -  check ( noarc, "There are arcs between D and C!" );
   4.400 -
   4.401 -  bool oddcomp=true;
   4.402 -  Graph::NodeMap<bool> todo(g,true);
   4.403 -  int num_comp=0;
   4.404 -  for(NodeIt v(g); v!=INVALID; ++v) {
   4.405 -   if ( pos[v]==max_matching.D && todo[v] ) {
   4.406 -      int comp_size=1;
   4.407 -      ++num_comp;
   4.408 -      std::queue<Node> Q;
   4.409 -      Q.push(v);
   4.410 -      todo.set(v,false);
   4.411 -      while (!Q.empty()) {
   4.412 -        Node w=Q.front();
   4.413 -        Q.pop();
   4.414 -        for(IncEdgeIt e(g,w); e!=INVALID; ++e) {
   4.415 -          Node u=g.runningNode(e);
   4.416 -          if ( pos[u]==max_matching.D && todo[u] ) {
   4.417 -            ++comp_size;
   4.418 -            Q.push(u);
   4.419 -            todo.set(u,false);
   4.420 -          }
   4.421 -        }
   4.422 -      }
   4.423 -      if ( !(comp_size % 2) ) oddcomp=false;
   4.424 -    }
   4.425 -  }
   4.426 -  check ( oddcomp, "A component of g[D] is not odd." );
   4.427 -
   4.428 -  int barrier=0;
   4.429 -  for(NodeIt v(g); v!=INVALID; ++v) {
   4.430 -    if ( pos[v]==max_matching.A ) ++barrier;
   4.431 -  }
   4.432 -  int expected_size=int( countNodes(g)-num_comp+barrier)/2;
   4.433 -  check ( size==expected_size, "The size of the matching is wrong." );
   4.434  
   4.435    return 0;
   4.436  }
     5.1 --- a/test/max_weighted_matching_test.cc	Mon Oct 13 13:56:00 2008 +0200
     5.2 +++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
     5.3 @@ -1,250 +0,0 @@
     5.4 -/* -*- mode: C++; indent-tabs-mode: nil; -*-
     5.5 - *
     5.6 - * This file is a part of LEMON, a generic C++ optimization library.
     5.7 - *
     5.8 - * Copyright (C) 2003-2008
     5.9 - * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    5.10 - * (Egervary Research Group on Combinatorial Optimization, EGRES).
    5.11 - *
    5.12 - * Permission to use, modify and distribute this software is granted
    5.13 - * provided that this copyright notice appears in all copies. For
    5.14 - * precise terms see the accompanying LICENSE file.
    5.15 - *
    5.16 - * This software is provided "AS IS" with no warranty of any kind,
    5.17 - * express or implied, and with no claim as to its suitability for any
    5.18 - * purpose.
    5.19 - *
    5.20 - */
    5.21 -
    5.22 -#include <iostream>
    5.23 -#include <sstream>
    5.24 -#include <vector>
    5.25 -#include <queue>
    5.26 -#include <lemon/math.h>
    5.27 -#include <cstdlib>
    5.28 -
    5.29 -#include "test_tools.h"
    5.30 -#include <lemon/smart_graph.h>
    5.31 -#include <lemon/max_matching.h>
    5.32 -#include <lemon/lgf_reader.h>
    5.33 -
    5.34 -using namespace std;
    5.35 -using namespace lemon;
    5.36 -
    5.37 -GRAPH_TYPEDEFS(SmartGraph);
    5.38 -
    5.39 -const int lgfn = 3;
    5.40 -const std::string lgf[lgfn] = {
    5.41 -  "@nodes\n"
    5.42 -  "label\n"
    5.43 -  "0\n"
    5.44 -  "1\n"
    5.45 -  "2\n"
    5.46 -  "3\n"
    5.47 -  "4\n"
    5.48 -  "5\n"
    5.49 -  "6\n"
    5.50 -  "7\n"
    5.51 -  "@edges\n"
    5.52 -  "label        weight\n"
    5.53 -  "7    4       0       984\n"
    5.54 -  "0    7       1       73\n"
    5.55 -  "7    1       2       204\n"
    5.56 -  "2    3       3       583\n"
    5.57 -  "2    7       4       565\n"
    5.58 -  "2    1       5       582\n"
    5.59 -  "0    4       6       551\n"
    5.60 -  "2    5       7       385\n"
    5.61 -  "1    5       8       561\n"
    5.62 -  "5    3       9       484\n"
    5.63 -  "7    5       10      904\n"
    5.64 -  "3    6       11      47\n"
    5.65 -  "7    6       12      888\n"
    5.66 -  "3    0       13      747\n"
    5.67 -  "6    1       14      310\n",
    5.68 -
    5.69 -  "@nodes\n"
    5.70 -  "label\n"
    5.71 -  "0\n"
    5.72 -  "1\n"
    5.73 -  "2\n"
    5.74 -  "3\n"
    5.75 -  "4\n"
    5.76 -  "5\n"
    5.77 -  "6\n"
    5.78 -  "7\n"
    5.79 -  "@edges\n"
    5.80 -  "label        weight\n"
    5.81 -  "2    5       0       710\n"
    5.82 -  "0    5       1       241\n"
    5.83 -  "2    4       2       856\n"
    5.84 -  "2    6       3       762\n"
    5.85 -  "4    1       4       747\n"
    5.86 -  "6    1       5       962\n"
    5.87 -  "4    7       6       723\n"
    5.88 -  "1    7       7       661\n"
    5.89 -  "2    3       8       376\n"
    5.90 -  "1    0       9       416\n"
    5.91 -  "6    7       10      391\n",
    5.92 -
    5.93 -  "@nodes\n"
    5.94 -  "label\n"
    5.95 -  "0\n"
    5.96 -  "1\n"
    5.97 -  "2\n"
    5.98 -  "3\n"
    5.99 -  "4\n"
   5.100 -  "5\n"
   5.101 -  "6\n"
   5.102 -  "7\n"
   5.103 -  "@edges\n"
   5.104 -  "label        weight\n"
   5.105 -  "6    2       0       553\n"
   5.106 -  "0    7       1       653\n"
   5.107 -  "6    3       2       22\n"
   5.108 -  "4    7       3       846\n"
   5.109 -  "7    2       4       981\n"
   5.110 -  "7    6       5       250\n"
   5.111 -  "5    2       6       539\n"
   5.112 -};
   5.113 -
   5.114 -void checkMatching(const SmartGraph& graph,
   5.115 -                   const SmartGraph::EdgeMap<int>& weight,
   5.116 -                   const MaxWeightedMatching<SmartGraph>& mwm) {
   5.117 -  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
   5.118 -    if (graph.u(e) == graph.v(e)) continue;
   5.119 -    int rw =
   5.120 -      mwm.nodeValue(graph.u(e)) + mwm.nodeValue(graph.v(e));
   5.121 -
   5.122 -    for (int i = 0; i < mwm.blossomNum(); ++i) {
   5.123 -      bool s = false, t = false;
   5.124 -      for (MaxWeightedMatching<SmartGraph>::BlossomIt n(mwm, i);
   5.125 -           n != INVALID; ++n) {
   5.126 -        if (graph.u(e) == n) s = true;
   5.127 -        if (graph.v(e) == n) t = true;
   5.128 -      }
   5.129 -      if (s == true && t == true) {
   5.130 -        rw += mwm.blossomValue(i);
   5.131 -      }
   5.132 -    }
   5.133 -    rw -= weight[e] * mwm.dualScale;
   5.134 -
   5.135 -    check(rw >= 0, "Negative reduced weight");
   5.136 -    check(rw == 0 || !mwm.matching(e),
   5.137 -          "Non-zero reduced weight on matching arc");
   5.138 -  }
   5.139 -
   5.140 -  int pv = 0;
   5.141 -  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   5.142 -    if (mwm.matching(n) != INVALID) {
   5.143 -      check(mwm.nodeValue(n) >= 0, "Invalid node value");
   5.144 -      pv += weight[mwm.matching(n)];
   5.145 -      SmartGraph::Node o = graph.target(mwm.matching(n));
   5.146 -      check(mwm.mate(n) == o, "Invalid matching");
   5.147 -      check(mwm.matching(n) == graph.oppositeArc(mwm.matching(o)),
   5.148 -            "Invalid matching");
   5.149 -    } else {
   5.150 -      check(mwm.mate(n) == INVALID, "Invalid matching");
   5.151 -      check(mwm.nodeValue(n) == 0, "Invalid matching");
   5.152 -    }
   5.153 -  }
   5.154 -
   5.155 -  int dv = 0;
   5.156 -  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   5.157 -    dv += mwm.nodeValue(n);
   5.158 -  }
   5.159 -
   5.160 -  for (int i = 0; i < mwm.blossomNum(); ++i) {
   5.161 -    check(mwm.blossomValue(i) >= 0, "Invalid blossom value");
   5.162 -    check(mwm.blossomSize(i) % 2 == 1, "Even blossom size");
   5.163 -    dv += mwm.blossomValue(i) * ((mwm.blossomSize(i) - 1) / 2);
   5.164 -  }
   5.165 -
   5.166 -  check(pv * mwm.dualScale == dv * 2, "Wrong duality");
   5.167 -
   5.168 -  return;
   5.169 -}
   5.170 -
   5.171 -void checkPerfectMatching(const SmartGraph& graph,
   5.172 -                          const SmartGraph::EdgeMap<int>& weight,
   5.173 -                          const MaxWeightedPerfectMatching<SmartGraph>& mwpm) {
   5.174 -  for (SmartGraph::EdgeIt e(graph); e != INVALID; ++e) {
   5.175 -    if (graph.u(e) == graph.v(e)) continue;
   5.176 -    int rw =
   5.177 -      mwpm.nodeValue(graph.u(e)) + mwpm.nodeValue(graph.v(e));
   5.178 -
   5.179 -    for (int i = 0; i < mwpm.blossomNum(); ++i) {
   5.180 -      bool s = false, t = false;
   5.181 -      for (MaxWeightedPerfectMatching<SmartGraph>::BlossomIt n(mwpm, i);
   5.182 -           n != INVALID; ++n) {
   5.183 -        if (graph.u(e) == n) s = true;
   5.184 -        if (graph.v(e) == n) t = true;
   5.185 -      }
   5.186 -      if (s == true && t == true) {
   5.187 -        rw += mwpm.blossomValue(i);
   5.188 -      }
   5.189 -    }
   5.190 -    rw -= weight[e] * mwpm.dualScale;
   5.191 -
   5.192 -    check(rw >= 0, "Negative reduced weight");
   5.193 -    check(rw == 0 || !mwpm.matching(e),
   5.194 -          "Non-zero reduced weight on matching arc");
   5.195 -  }
   5.196 -
   5.197 -  int pv = 0;
   5.198 -  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   5.199 -    check(mwpm.matching(n) != INVALID, "Non perfect");
   5.200 -    pv += weight[mwpm.matching(n)];
   5.201 -    SmartGraph::Node o = graph.target(mwpm.matching(n));
   5.202 -    check(mwpm.mate(n) == o, "Invalid matching");
   5.203 -    check(mwpm.matching(n) == graph.oppositeArc(mwpm.matching(o)),
   5.204 -          "Invalid matching");
   5.205 -  }
   5.206 -
   5.207 -  int dv = 0;
   5.208 -  for (SmartGraph::NodeIt n(graph); n != INVALID; ++n) {
   5.209 -    dv += mwpm.nodeValue(n);
   5.210 -  }
   5.211 -
   5.212 -  for (int i = 0; i < mwpm.blossomNum(); ++i) {
   5.213 -    check(mwpm.blossomValue(i) >= 0, "Invalid blossom value");
   5.214 -    check(mwpm.blossomSize(i) % 2 == 1, "Even blossom size");
   5.215 -    dv += mwpm.blossomValue(i) * ((mwpm.blossomSize(i) - 1) / 2);
   5.216 -  }
   5.217 -
   5.218 -  check(pv * mwpm.dualScale == dv * 2, "Wrong duality");
   5.219 -
   5.220 -  return;
   5.221 -}
   5.222 -
   5.223 -
   5.224 -int main() {
   5.225 -
   5.226 -  for (int i = 0; i < lgfn; ++i) {
   5.227 -    SmartGraph graph;
   5.228 -    SmartGraph::EdgeMap<int> weight(graph);
   5.229 -
   5.230 -    istringstream lgfs(lgf[i]);
   5.231 -    graphReader(graph, lgfs).
   5.232 -      edgeMap("weight", weight).run();
   5.233 -
   5.234 -    MaxWeightedMatching<SmartGraph> mwm(graph, weight);
   5.235 -    mwm.run();
   5.236 -    checkMatching(graph, weight, mwm);
   5.237 -
   5.238 -    MaxMatching<SmartGraph> mm(graph);
   5.239 -    mm.run();
   5.240 -
   5.241 -    MaxWeightedPerfectMatching<SmartGraph> mwpm(graph, weight);
   5.242 -    bool perfect = mwpm.run();
   5.243 -
   5.244 -    check(perfect == (mm.size() * 2 == countNodes(graph)),
   5.245 -          "Perfect matching found");
   5.246 -
   5.247 -    if (perfect) {
   5.248 -      checkPerfectMatching(graph, weight, mwpm);
   5.249 -    }
   5.250 -  }
   5.251 -
   5.252 -  return 0;
   5.253 -}