RhodeCode

  }

  /// \ingroup connectivity

  ///

  /// \brief Find the connected components of an undirected graph

  ///

  /// Find the connected components of an undirected graph.

  ///

  /// \param graph The graph. It must be undirected.

  /// \retval compMap A writable node map. The values will be set from 0 to

  /// the number of the connected components minus one. Each values of the map

  /// will be set exactly once, the values of a certain component will be

  /// set continuously.

  /// \return The number of components

  ///

  template <class Graph, class NodeMap>

  int connectedComponents(const Graph &graph, NodeMap &compMap) {

    checkConcept<concepts::Graph, Graph>();

    typedef typename Graph::Node Node;

    typedef typename Graph::Arc Arc;

    checkConcept<concepts::WriteMap<Node, int>, NodeMap>();

    typedef NullMap<Node, Arc> PredMap;

    typedef NullMap<Node, int> DistMap;

    int compNum = 0;

    typename Bfs<Graph>::

      template SetPredMap<PredMap>::

      template SetDistMap<DistMap>::

      Create bfs(graph);

    PredMap predMap;

    bfs.predMap(predMap);

    DistMap distMap;

    bfs.distMap(distMap);

    bfs.init();

    for(typename Graph::NodeIt n(graph); n != INVALID; ++n) {

      if(!bfs.reached(n)) {

        bfs.addSource(n);

        while (!bfs.emptyQueue()) {

          compMap.set(bfs.nextNode(), compNum);

          bfs.processNextNode();

        }

        ++compNum;

      }

    }

    return compNum;

  }

  namespace _connectivity_bits {

    template <typename Digraph, typename Iterator >

    struct LeaveOrderVisitor : public DfsVisitor<Digraph> {

    public:

      typedef typename Digraph::Node Node;

      LeaveOrderVisitor(Iterator it) : _it(it) {}

      void leave(const Node& node) {

        *(_it++) = node;

      }

    private:

      Iterator _it;

    };

    template <typename Digraph, typename Map>

    struct FillMapVisitor : public DfsVisitor<Digraph> {

    public:

      typedef typename Digraph::Node Node;

      typedef typename Map::Value Value;

      FillMapVisitor(Map& map, Value& value)

        : _map(map), _value(value) {}

      void reach(const Node& node) {

        _map.set(node, _value);

      }

    private:

      Map& _map;

      Value& _value;

    };

    template <typename Digraph, typename ArcMap>

    struct StronglyConnectedCutArcsVisitor : public DfsVisitor<Digraph> {

    public:

      typedef typename Digraph::Node Node;

      typedef typename Digraph::Arc Arc;

      StronglyConnectedCutArcsVisitor(const Digraph& digraph,

                                      ArcMap& cutMap,

                                      int& cutNum)

        : _digraph(digraph), _cutMap(cutMap), _cutNum(cutNum),

          _compMap(digraph, -1), _num(-1) {

      }

      void start(const Node&) {

        ++_num;

      }

      void reach(const Node& node) {

        _compMap.set(node, _num);

      }

      void examine(const Arc& arc) {

         if (_compMap[_digraph.source(arc)] !=

             _compMap[_digraph.target(arc)]) {

           _cutMap.set(arc, true);

           ++_cutNum;

         }

      }

    private:

      const Digraph& _digraph;

      ArcMap& _cutMap;

      int& _cutNum;

      typename Digraph::template NodeMap<int> _compMap;

      int _num;

    };

  }

  /// \ingroup connectivity

  ///

  /// \brief Check whether the given directed graph is strongly connected.

  ///

  /// Check whether the given directed graph is strongly connected. The

  /// graph is strongly connected when any two nodes of the graph are

  /// connected with directed paths in both direction.

  /// \return %False when the graph is not strongly connected.

  /// \see connected

  ///

  /// \note By definition, the empty graph is strongly connected.

  template <typename Digraph>

  bool stronglyConnected(const Digraph& digraph) {

    checkConcept<concepts::Digraph, Digraph>();

    typedef typename Digraph::Node Node;

    typedef typename Digraph::NodeIt NodeIt;

    typename Digraph::Node source = NodeIt(digraph);

    if (source == INVALID) return true;

    using namespace _connectivity_bits;

    typedef DfsVisitor<Digraph> Visitor;

    Visitor visitor;

    DfsVisit<Digraph, Visitor> dfs(digraph, visitor);

    dfs.init();

    dfs.addSource(source);

    dfs.start();

    for (NodeIt it(digraph); it != INVALID; ++it) {

      if (!dfs.reached(it)) {

        return false;

      }

    }

    typedef ReverseDigraph<const Digraph> RDigraph;

    typedef typename RDigraph::NodeIt RNodeIt;

    RDigraph rdigraph(digraph);

    typedef DfsVisitor<Digraph> RVisitor;

    RVisitor rvisitor;

    DfsVisit<RDigraph, RVisitor> rdfs(rdigraph, rvisitor);

    rdfs.init();

    rdfs.addSource(source);

    rdfs.start();

    for (RNodeIt it(rdigraph); it != INVALID; ++it) {

      if (!rdfs.reached(it)) {

        return false;

      }

    }

    return true;

  }

  /// \ingroup connectivity

  ///

  /// \brief Count the strongly connected components of a directed graph

  ///

  /// Count the strongly connected components of a directed graph.

  /// The strongly connected components are the classes of an

  /// equivalence relation on the nodes of the graph. Two nodes are in

  /// the same class if they are connected with directed paths in both

  /// direction.

  ///

  /// \return The number of components

  /// \note By definition, the empty graph has zero

  /// strongly connected components.

  template <typename Digraph>

  int countStronglyConnectedComponents(const Digraph& digraph) {

    checkConcept<concepts::Digraph, Digraph>();

    using namespace _connectivity_bits;

    typedef typename Digraph::Node Node;

    typedef typename Digraph::Arc Arc;

    typedef typename Digraph::NodeIt NodeIt;

    typedef typename Digraph::ArcIt ArcIt;

    typedef std::vector<Node> Container;

    typedef typename Container::iterator Iterator;

    Container nodes(countNodes(digraph));

    typedef LeaveOrderVisitor<Digraph, Iterator> Visitor;

    Visitor visitor(nodes.begin());

    DfsVisit<Digraph, Visitor> dfs(digraph, visitor);

    dfs.init();

    for (NodeIt it(digraph); it != INVALID; ++it) {

      if (!dfs.reached(it)) {

        dfs.addSource(it);

        dfs.start();

      }

    }

    typedef typename Container::reverse_iterator RIterator;

    typedef ReverseDigraph<const Digraph> RDigraph;

    RDigraph rdigraph(digraph);

    typedef DfsVisitor<Digraph> RVisitor;

    RVisitor rvisitor;

    DfsVisit<RDigraph, RVisitor> rdfs(rdigraph, rvisitor);

    int compNum = 0;

    rdfs.init();

    for (RIterator it = nodes.rbegin(); it != nodes.rend(); ++it) {

      if (!rdfs.reached(*it)) {

        rdfs.addSource(*it);

        rdfs.start();

        ++compNum;

      }

    }

    return compNum;

  }

  /// \ingroup connectivity

  ///

  /// \brief Find the strongly connected components of a directed graph

  ///

  /// Find the strongly connected components of a directed graph.  The

  /// strongly connected components are the classes of an equivalence

  /// relation on the nodes of the graph. Two nodes are in

  /// relationship when there are directed paths between them in both

  /// direction. In addition, the numbering of components will satisfy

  /// that there is no arc going from a higher numbered component to

  /// a lower.

  ///

  /// \param digraph The digraph.

  /// \retval compMap A writable node map. The values will be set from 0 to

  /// the number of the strongly connected components minus one. Each value

  /// of the map will be set exactly once, the values of a certain component

  /// will be set continuously.

  /// \return The number of components

  ///

  template <typename Digraph, typename NodeMap>

  int stronglyConnectedComponents(const Digraph& digraph, NodeMap& compMap) {

    checkConcept<concepts::Digraph, Digraph>();

    typedef typename Digraph::Node Node;

    typedef typename Digraph::NodeIt NodeIt;

    checkConcept<concepts::WriteMap<Node, int>, NodeMap>();

    using namespace _connectivity_bits;

    typedef std::vector<Node> Container;

    typedef typename Container::iterator Iterator;

    Container nodes(countNodes(digraph));

    typedef LeaveOrderVisitor<Digraph, Iterator> Visitor;

    Visitor visitor(nodes.begin());

    DfsVisit<Digraph, Visitor> dfs(digraph, visitor);

    dfs.init();

    for (NodeIt it(digraph); it != INVALID; ++it) {

      if (!dfs.reached(it)) {

        dfs.addSource(it);

        dfs.start();

      }

    }

    typedef typename Container::reverse_iterator RIterator;

    typedef ReverseDigraph<const Digraph> RDigraph;

    RDigraph rdigraph(digraph);

    int compNum = 0;

    typedef FillMapVisitor<RDigraph, NodeMap> RVisitor;

    RVisitor rvisitor(compMap, compNum);

    DfsVisit<RDigraph, RVisitor> rdfs(rdigraph, rvisitor);

    rdfs.init();

    for (RIterator it = nodes.rbegin(); it != nodes.rend(); ++it) {

      if (!rdfs.reached(*it)) {

        rdfs.addSource(*it);

        rdfs.start();

        ++compNum;

      }

    }

    return compNum;

  }

  /// \ingroup connectivity

  ///

  /// \brief Find the cut arcs of the strongly connected components.

  ///

  /// Find the cut arcs of the strongly connected components.

  /// The strongly connected components are the classes of an equivalence

  /// relation on the nodes of the graph. Two nodes are in relationship

  /// when there are directed paths between them in both direction.

  /// The strongly connected components are separated by the cut arcs.

  ///

  /// \param graph The graph.

  /// \retval cutMap A writable node map. The values will be set true when the

  /// arc is a cut arc.

  ///

  /// \return The number of cut arcs

  template <typename Digraph, typename ArcMap>

  int stronglyConnectedCutArcs(const Digraph& graph, ArcMap& cutMap) {

    checkConcept<concepts::Digraph, Digraph>();

    typedef typename Digraph::Node Node;

    typedef typename Digraph::Arc Arc;

    typedef typename Digraph::NodeIt NodeIt;

    checkConcept<concepts::WriteMap<Arc, bool>, ArcMap>();

    using namespace _connectivity_bits;

    typedef std::vector<Node> Container;

    typedef typename Container::iterator Iterator;

    Container nodes(countNodes(graph));

    typedef LeaveOrderVisitor<Digraph, Iterator> Visitor;

    Visitor visitor(nodes.begin());

    DfsVisit<Digraph, Visitor> dfs(graph, visitor);

    dfs.init();

    for (NodeIt it(graph); it != INVALID; ++it) {

      if (!dfs.reached(it)) {

        dfs.addSource(it);

        dfs.start();

      }

    }

    typedef typename Container::reverse_iterator RIterator;

    typedef ReverseDigraph<const Digraph> RDigraph;

    RDigraph rgraph(graph);

    int cutNum = 0;

    typedef StronglyConnectedCutArcsVisitor<RDigraph, ArcMap> RVisitor;

    RVisitor rvisitor(rgraph, cutMap, cutNum);

    DfsVisit<RDigraph, RVisitor> rdfs(rgraph, rvisitor);

    rdfs.init();

    for (RIterator it = nodes.rbegin(); it != nodes.rend(); ++it) {

      if (!rdfs.reached(*it)) {

        rdfs.addSource(*it);

        rdfs.start();

      }

    }

    return cutNum;

  }

  namespace _connectivity_bits {

    template <typename Digraph>

    class CountBiNodeConnectedComponentsVisitor : public DfsVisitor<Digraph> {

    public:

      typedef typename Digraph::Node Node;

      typedef typename Digraph::Arc Arc;

      typedef typename Digraph::Edge Edge;

  template <typename Graph>

  bool biEdgeConnected(const Graph& graph) {

    return countBiEdgeConnectedComponents(graph) <= 1;

  }

  /// \ingroup connectivity

  ///

  /// \brief Count the bi-edge-connected components.

  ///

  /// This function count the bi-edge-connected components in an undirected

  /// graph. The bi-edge-connected components are the classes of an equivalence

  /// relation on the nodes. Two nodes are in relationship when they are

  /// connected with at least two edge-disjoint paths.

  ///

  /// \param graph The undirected graph.

  /// \return The number of components.

  template <typename Graph>

  int countBiEdgeConnectedComponents(const Graph& graph) {

    checkConcept<concepts::Graph, Graph>();

    typedef typename Graph::NodeIt NodeIt;

    using namespace _connectivity_bits;

    typedef CountBiEdgeConnectedComponentsVisitor<Graph> Visitor;

    int compNum = 0;

    Visitor visitor(graph, compNum);

    DfsVisit<Graph, Visitor> dfs(graph, visitor);

    dfs.init();

    for (NodeIt it(graph); it != INVALID; ++it) {

      if (!dfs.reached(it)) {

        dfs.addSource(it);

        dfs.start();

      }

    }

    return compNum;

  }

  /// \ingroup connectivity

  ///

  /// \brief Find the bi-edge-connected components.

  ///

  /// This function finds the bi-edge-connected components in an undirected

  /// graph. The bi-edge-connected components are the classes of an equivalence

  /// relation on the nodes. Two nodes are in relationship when they are

  /// connected at least two edge-disjoint paths.

  ///

  /// \param graph The graph.

  /// \retval compMap A writable node map. The values will be set from 0 to

  /// the number of the biconnected components minus one. Each values

  /// of the map will be set exactly once, the values of a certain component

  /// will be set continuously.

  /// \return The number of components.

  ///

  template <typename Graph, typename NodeMap>

  int biEdgeConnectedComponents(const Graph& graph, NodeMap& compMap) {

    checkConcept<concepts::Graph, Graph>();

    typedef typename Graph::NodeIt NodeIt;

    typedef typename Graph::Node Node;

    checkConcept<concepts::WriteMap<Node, int>, NodeMap>();

    using namespace _connectivity_bits;

    typedef BiEdgeConnectedComponentsVisitor<Graph, NodeMap> Visitor;

    int compNum = 0;

    Visitor visitor(graph, compMap, compNum);

    DfsVisit<Graph, Visitor> dfs(graph, visitor);

    dfs.init();

    for (NodeIt it(graph); it != INVALID; ++it) {

      if (!dfs.reached(it)) {

        dfs.addSource(it);

        dfs.start();

      }

    }

    return compNum;

  }

  /// \ingroup connectivity

  ///

  /// \brief Find the bi-edge-connected cut edges.

  ///

  /// This function finds the bi-edge-connected components in an undirected

  /// graph. The bi-edge-connected components are the classes of an equivalence

  /// relation on the nodes. Two nodes are in relationship when they are

  /// connected with at least two edge-disjoint paths. The bi-edge-connected

  /// components are separted by edges which are the cut edges of the

  /// components.

  ///

  /// \param graph The graph.

  /// \retval cutMap A writable node map. The values will be set true when the

  /// edge is a cut edge.

  /// \return The number of cut edges.

  template <typename Graph, typename EdgeMap>

  int biEdgeConnectedCutEdges(const Graph& graph, EdgeMap& cutMap) {

    checkConcept<concepts::Graph, Graph>();

    typedef typename Graph::NodeIt NodeIt;

    typedef typename Graph::Edge Edge;

    checkConcept<concepts::WriteMap<Edge, bool>, EdgeMap>();

    using namespace _connectivity_bits;

    typedef BiEdgeConnectedCutEdgesVisitor<Graph, EdgeMap> Visitor;

    int cutNum = 0;

    Visitor visitor(graph, cutMap, cutNum);

    DfsVisit<Graph, Visitor> dfs(graph, visitor);

    dfs.init();

    for (NodeIt it(graph); it != INVALID; ++it) {

      if (!dfs.reached(it)) {

        dfs.addSource(it);

        dfs.start();

      }

    }

    return cutNum;

  }

  namespace _connectivity_bits {

    template <typename Digraph, typename IntNodeMap>

    class TopologicalSortVisitor : public DfsVisitor<Digraph> {

    public:

      typedef typename Digraph::Node Node;

      typedef typename Digraph::Arc edge;

      TopologicalSortVisitor(IntNodeMap& order, int num)

        : _order(order), _num(num) {}

      void leave(const Node& node) {

        _order.set(node, --_num);

      }

    private:

      IntNodeMap& _order;

      int _num;

    };

  }

  /// \ingroup connectivity

  ///

  /// \brief Sort the nodes of a DAG into topolgical order.

  ///

  /// Sort the nodes of a DAG into topolgical order.

  ///

  /// \param graph The graph. It must be directed and acyclic.

  /// \retval order A writable node map. The values will be set from 0 to

  /// the number of the nodes in the graph minus one. Each values of the map

  /// will be set exactly once, the values  will be set descending order.

  ///

  /// \see checkedTopologicalSort

  /// \see dag

  template <typename Digraph, typename NodeMap>

  void topologicalSort(const Digraph& graph, NodeMap& order) {

    using namespace _connectivity_bits;

    checkConcept<concepts::Digraph, Digraph>();

    checkConcept<concepts::WriteMap<typename Digraph::Node, int>, NodeMap>();

    typedef typename Digraph::Node Node;

    typedef typename Digraph::NodeIt NodeIt;

    typedef typename Digraph::Arc Arc;

    TopologicalSortVisitor<Digraph, NodeMap>

      visitor(order, countNodes(graph));

    DfsVisit<Digraph, TopologicalSortVisitor<Digraph, NodeMap> >

      dfs(graph, visitor);

    dfs.init();

    for (NodeIt it(graph); it != INVALID; ++it) {

      if (!dfs.reached(it)) {

        dfs.addSource(it);

        dfs.start();

      }

    }

  }

  /// \ingroup connectivity

  ///

  /// \brief Sort the nodes of a DAG into topolgical order.

  ///

  /// Sort the nodes of a DAG into topolgical order. It also checks

  /// that the given graph is DAG.

  ///

  /// \retval order A readable - writable node map. The values will be set

  /// from 0 to the number of the nodes in the graph minus one. Each values

  /// of the map will be set exactly once, the values will be set descending

  /// order.

  /// \return %False when the graph is not DAG.

  ///

  /// \see topologicalSort

  /// \see dag

  template <typename Digraph, typename NodeMap>

  bool checkedTopologicalSort(const Digraph& digraph, NodeMap& order) {

    using namespace _connectivity_bits;

    checkConcept<concepts::Digraph, Digraph>();

    checkConcept<concepts::ReadWriteMap<typename Digraph::Node, int>,

      NodeMap>();

    typedef typename Digraph::Node Node;

    typedef typename Digraph::NodeIt NodeIt;

    typedef typename Digraph::Arc Arc;

    for (NodeIt it(digraph); it != INVALID; ++it) {

      order.set(it, -1);

    }

    TopologicalSortVisitor<Digraph, NodeMap>

      visitor(order, countNodes(digraph));

    DfsVisit<Digraph, TopologicalSortVisitor<Digraph, NodeMap> >

      dfs(digraph, visitor);

    dfs.init();

    for (NodeIt it(digraph); it != INVALID; ++it) {

      if (!dfs.reached(it)) {

        dfs.addSource(it);

        while (!dfs.emptyQueue()) {

           Arc arc = dfs.nextArc();

           Node target = digraph.target(arc);

           if (dfs.reached(target) && order[target] == -1) {

             return false;

           }

           dfs.processNextArc();

         }

      }

    }

    return true;

  }

  /// \ingroup connectivity

  ///

  /// \brief Check that the given directed graph is a DAG.

  ///

  /// Check that the given directed graph is a DAG. The DAG is

  /// an Directed Acyclic Digraph.

  /// \return %False when the graph is not DAG.

  /// \see acyclic

  template <typename Digraph>

  bool dag(const Digraph& digraph) {

    checkConcept<concepts::Digraph, Digraph>();

    typedef typename Digraph::Node Node;

    typedef typename Digraph::NodeIt NodeIt;

    typedef typename Digraph::Arc Arc;

    typedef typename Digraph::template NodeMap<bool> ProcessedMap;

    typename Dfs<Digraph>::template SetProcessedMap<ProcessedMap>::

      Create dfs(digraph);

    ProcessedMap processed(digraph);

    dfs.processedMap(processed);

    dfs.init();

    for (NodeIt it(digraph); it != INVALID; ++it) {

      if (!dfs.reached(it)) {

        dfs.addSource(it);

        while (!dfs.emptyQueue()) {

          Arc edge = dfs.nextArc();

          Node target = digraph.target(edge);

          if (dfs.reached(target) && !processed[target]) {

            return false;

          }

          dfs.processNextArc();

        }

      }

    }

    return true;

  }

  /// \ingroup connectivity

  ///

  /// \brief Check that the given undirected graph is acyclic.

  ///

  /// Check that the given undirected graph acyclic.

  /// \param graph The undirected graph.

  /// \return %True when there is no circle in the graph.

  /// \see dag

  template <typename Graph>

  bool acyclic(const Graph& graph) {

    checkConcept<concepts::Graph, Graph>();

    typedef typename Graph::Node Node;

    typedef typename Graph::NodeIt NodeIt;

    typedef typename Graph::Arc Arc;

    Dfs<Graph> dfs(graph);

    dfs.init();

    for (NodeIt it(graph); it != INVALID; ++it) {

      if (!dfs.reached(it)) {

        dfs.addSource(it);

        while (!dfs.emptyQueue()) {

          Arc edge = dfs.nextArc();

          Node source = graph.source(edge);

          Node target = graph.target(edge);

          if (dfs.reached(target) &&

              dfs.predArc(source) != graph.oppositeArc(edge)) {

            return false;

          }

          dfs.processNextArc();

        }

      }

    }

    return true;

  }

  /// \ingroup connectivity

  ///

  /// \brief Check that the given undirected graph is tree.

  ///

  /// Check that the given undirected graph is tree.

  /// \param graph The undirected graph.

  /// \return %True when the graph is acyclic and connected.

  template <typename Graph>

  bool tree(const Graph& graph) {

    checkConcept<concepts::Graph, Graph>();

    typedef typename Graph::Node Node;

    typedef typename Graph::NodeIt NodeIt;

    typedef typename Graph::Arc Arc;

    Dfs<Graph> dfs(graph);

    dfs.init();

    dfs.addSource(NodeIt(graph));

    while (!dfs.emptyQueue()) {

      Arc edge = dfs.nextArc();

      Node source = graph.source(edge);

      Node target = graph.target(edge);

      if (dfs.reached(target) &&

          dfs.predArc(source) != graph.oppositeArc(edge)) {

        return false;

      }

      dfs.processNextArc();

    }

    for (NodeIt it(graph); it != INVALID; ++it) {

      if (!dfs.reached(it)) {

        return false;

      }

    }

    return true;

  }

  namespace _connectivity_bits {

    template <typename Digraph>

    class BipartiteVisitor : public BfsVisitor<Digraph> {

    public:

      typedef typename Digraph::Arc Arc;

      typedef typename Digraph::Node Node;

      BipartiteVisitor(const Digraph& graph, bool& bipartite)

        : _graph(graph), _part(graph), _bipartite(bipartite) {}

      void start(const Node& node) {

        _part[node] = true;

      }

      void discover(const Arc& edge) {

        _part.set(_graph.target(edge), !_part[_graph.source(edge)]);

      }

      void examine(const Arc& edge) {

        _bipartite = _bipartite &&

          _part[_graph.target(edge)] != _part[_graph.source(edge)];

      }

    private:

      const Digraph& _graph;

      typename Digraph::template NodeMap<bool> _part;

      bool& _bipartite;

    };

    template <typename Digraph, typename PartMap>

    class BipartitePartitionsVisitor : public BfsVisitor<Digraph> {

    public:

      typedef typename Digraph::Arc Arc;

      typedef typename Digraph::Node Node;

      {

        std::set<Node> left_set, right_set;

        Node left = (*_blossom_rep)[_blossom_set->find(_graph.u(e))];

        left_set.insert(left);

        Node right = (*_blossom_rep)[_blossom_set->find(_graph.v(e))];

        right_set.insert(right);

        while (true) {

          if ((*_matching)[left] == INVALID) break;

          left = _graph.target((*_matching)[left]);

          left = (*_blossom_rep)[_blossom_set->

                                 find(_graph.target((*_ear)[left]))];

          if (right_set.find(left) != right_set.end()) {

            nca = left;

            break;

          }

          left_set.insert(left);

          if ((*_matching)[right] == INVALID) break;

          right = _graph.target((*_matching)[right]);

          right = (*_blossom_rep)[_blossom_set->

                                  find(_graph.target((*_ear)[right]))];

          if (left_set.find(right) != left_set.end()) {

            nca = right;

            break;

          }

          right_set.insert(right);

        }

        if (nca == INVALID) {

          if ((*_matching)[left] == INVALID) {

            nca = right;

            while (left_set.find(nca) == left_set.end()) {

              nca = _graph.target((*_matching)[nca]);

              nca =(*_blossom_rep)[_blossom_set->

                                   find(_graph.target((*_ear)[nca]))];

            }

          } else {

            nca = left;

            while (right_set.find(nca) == right_set.end()) {

              nca = _graph.target((*_matching)[nca]);

              nca = (*_blossom_rep)[_blossom_set->

                                   find(_graph.target((*_ear)[nca]))];

            }

          }

        }

      }

      {

        Node node = _graph.u(e);

        Arc arc = _graph.direct(e, true);

        Node base = (*_blossom_rep)[_blossom_set->find(node)];

        while (base != nca) {

          _ear->set(node, arc);

          Node n = node;

          while (n != base) {

            n = _graph.target((*_matching)[n]);

            Arc a = (*_ear)[n];

            n = _graph.target(a);

            _ear->set(n, _graph.oppositeArc(a));

          }

          node = _graph.target((*_matching)[base]);

          _tree_set->erase(base);

          _tree_set->erase(node);

          _blossom_set->insert(node, _blossom_set->find(base));

          _status->set(node, EVEN);

          _node_queue[_last++] = node;

          arc = _graph.oppositeArc((*_ear)[node]);

          node = _graph.target((*_ear)[node]);

          base = (*_blossom_rep)[_blossom_set->find(node)];

          _blossom_set->join(_graph.target(arc), base);

        }

      }

      _blossom_rep->set(_blossom_set->find(nca), nca);

      {

        Node node = _graph.v(e);

        Arc arc = _graph.direct(e, false);

        Node base = (*_blossom_rep)[_blossom_set->find(node)];

        while (base != nca) {

          _ear->set(node, arc);

          Node n = node;

          while (n != base) {

            n = _graph.target((*_matching)[n]);

            Arc a = (*_ear)[n];

            n = _graph.target(a);

            _ear->set(n, _graph.oppositeArc(a));

          }

          node = _graph.target((*_matching)[base]);

          _tree_set->erase(base);

          _tree_set->erase(node);

          _blossom_set->insert(node, _blossom_set->find(base));

          _status->set(node, EVEN);

          _node_queue[_last++] = node;

          arc = _graph.oppositeArc((*_ear)[node]);

          node = _graph.target((*_ear)[node]);

          base = (*_blossom_rep)[_blossom_set->find(node)];

          _blossom_set->join(_graph.target(arc), base);

        }

      }

      _blossom_rep->set(_blossom_set->find(nca), nca);

    }

    void extendOnArc(const Arc& a) {

      Node base = _graph.source(a);

      Node odd = _graph.target(a);

      _ear->set(odd, _graph.oppositeArc(a));

      Node even = _graph.target((*_matching)[odd]);

      _blossom_rep->set(_blossom_set->insert(even), even);

      _status->set(odd, ODD);

      _status->set(even, EVEN);

      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(base)]);

      _tree_set->insert(odd, tree);

      _tree_set->insert(even, tree);

      _node_queue[_last++] = even;

    }

    void augmentOnArc(const Arc& a) {

      Node even = _graph.source(a);

      Node odd = _graph.target(a);

      int tree = _tree_set->find((*_blossom_rep)[_blossom_set->find(even)]);

      _matching->set(odd, _graph.oppositeArc(a));

      _status->set(odd, MATCHED);

      Arc arc = (*_matching)[even];

      _matching->set(even, a);

      while (arc != INVALID) {

        odd = _graph.target(arc);

        arc = (*_ear)[odd];

        even = _graph.target(arc);

        _matching->set(odd, arc);

        arc = (*_matching)[even];

        _matching->set(even, _graph.oppositeArc((*_matching)[odd]));

      }

      for (typename TreeSet::ItemIt it(*_tree_set, tree);

           it != INVALID; ++it) {

        if ((*_status)[it] == ODD) {

          _status->set(it, MATCHED);

        } else {

          int blossom = _blossom_set->find(it);

          for (typename BlossomSet::ItemIt jt(*_blossom_set, blossom);

               jt != INVALID; ++jt) {

            _status->set(jt, MATCHED);

          }

          _blossom_set->eraseClass(blossom);

        }

      }

      _tree_set->eraseClass(tree);

    }

  public:

    /// \brief Constructor

    ///

    /// Constructor.

    MaxMatching(const Graph& graph)

      : _graph(graph), _matching(0), _status(0), _ear(0),

        _blossom_set_index(0), _blossom_set(0), _blossom_rep(0),

        _tree_set_index(0), _tree_set(0) {}

    ~MaxMatching() {

      destroyStructures();

    }

    /// \name Execution control

    /// The simplest way to execute the algorithm is to use the

    /// \c run() member function.

    /// \n

    /// If you need better control on the execution, you must call

    /// \ref init(), \ref greedyInit() or \ref matchingInit()

    /// functions first, then you can start the algorithm with the \ref

    ///@{

    /// \brief Sets the actual matching to the empty matching.

    ///

    /// Sets the actual matching to the empty matching.

    ///

    void init() {

      createStructures();

      for(NodeIt n(_graph); n != INVALID; ++n) {

        _matching->set(n, INVALID);

        _status->set(n, UNMATCHED);

      }

    }

    ///\brief Finds an initial matching in a greedy way

    ///

    ///It finds an initial matching in a greedy way.

    void greedyInit() {

      createStructures();

      for (NodeIt n(_graph); n != INVALID; ++n) {

        _matching->set(n, INVALID);

        _status->set(n, UNMATCHED);

      }

      for (NodeIt n(_graph); n != INVALID; ++n) {

        if ((*_matching)[n] == INVALID) {

          for (OutArcIt a(_graph, n); a != INVALID ; ++a) {

            Node v = _graph.target(a);

            if ((*_matching)[v] == INVALID && v != n) {

              _matching->set(n, a);

              _status->set(n, MATCHED);

              _matching->set(v, _graph.oppositeArc(a));

              _status->set(v, MATCHED);

              break;

            }

          }

        }

      }

    }

    /// \brief Initialize the matching from a map containing.

    ///

    /// Initialize the matching from a \c bool valued \c Edge map. This

    /// map must have the property that there are no two incident edges

    /// with true value, ie. it contains a matching.

    /// \return %True if the map contains a matching.

    template <typename MatchingMap>

    bool matchingInit(const MatchingMap& matching) {

      createStructures();

      for (NodeIt n(_graph); n != INVALID; ++n) {

        _matching->set(n, INVALID);

        _status->set(n, UNMATCHED);

      }

      for(EdgeIt e(_graph); e!=INVALID; ++e) {

        if (matching[e]) {

          Node u = _graph.u(e);

          if ((*_matching)[u] != INVALID) return false;

          _matching->set(u, _graph.direct(e, true));

          _status->set(u, MATCHED);

          Node v = _graph.v(e);

          if ((*_matching)[v] != INVALID) return false;

          _matching->set(v, _graph.direct(e, false));

          _status->set(v, MATCHED);

        }

      }

      return true;

    }

    /// \brief Starts Edmonds' algorithm

    ///

    /// If runs the original Edmonds' algorithm.

    void startSparse() {

      for(NodeIt n(_graph); n != INVALID; ++n) {

        if ((*_status)[n] == UNMATCHED) {

          (*_blossom_rep)[_blossom_set->insert(n)] = n;

          _tree_set->insert(n);

          _status->set(n, EVEN);

          processSparse(n);

        }

      }

    }

    /// \brief Starts Edmonds' algorithm.

    ///

    /// It runs Edmonds' algorithm with a heuristic of postponing

    /// shrinks, therefore resulting in a faster algorithm for dense graphs.

    void startDense() {

      for(NodeIt n(_graph); n != INVALID; ++n) {

        if ((*_status)[n] == UNMATCHED) {

          (*_blossom_rep)[_blossom_set->insert(n)] = n;

          _tree_set->insert(n);

          _status->set(n, EVEN);

          processDense(n);

        }

      }

    }

    /// \brief Runs Edmonds' algorithm

    ///

    /// Runs Edmonds' algorithm for sparse graphs (<tt>m<2*n</tt>)

    /// or Edmonds' algorithm with a heuristic of

    /// postponing shrinks for dense graphs.

    void run() {

      if (countEdges(_graph) < 2 * countNodes(_graph)) {

        greedyInit();

        startSparse();

      } else {

        init();

        startDense();

      }

    }

    /// @}

    /// \name Primal solution

    /// Functions to get the primal solution, ie. the matching.

    /// @{

    ///\brief Returns the size of the current matching.

    ///

    ///Returns the size of the current matching. After \ref

    ///run() it returns the size of the maximum matching in the graph.

    int matchingSize() const {

      int size = 0;

      for (NodeIt n(_graph); n != INVALID; ++n) {

        if ((*_matching)[n] != INVALID) {

          ++size;

        }

      }

      return size / 2;

    }

    /// \brief Returns true when the edge is in the matching.

    ///

    /// Returns true when the edge is in the matching.

    bool matching(const Edge& edge) const {

      return edge == (*_matching)[_graph.u(edge)];

    }