RhodeCode

      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();

      }

    }

  /// \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.

            _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);

/* -*- C++ -*-

 *

 * This file is a part of LEMON, a generic C++ optimization library

 *

 * Copyright (C) 2003-2008

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

 * (Egervary Research Group on Combinatorial Optimization, EGRES).

 *

 * Permission to use, modify and distribute this software is granted

 * provided that this copyright notice appears in all copies. For

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

#ifndef LEMON_SUURBALLE_H

#define LEMON_SUURBALLE_H

///\ingroup shortest_path

///\file

///\brief An algorithm for finding arc-disjoint paths between two

/// nodes having minimum total length.

#include <vector>

#include <lemon/bin_heap.h>

#include <lemon/path.h>

namespace lemon {

  /// \addtogroup shortest_path

  /// @{

  /// \brief Algorithm for finding arc-disjoint paths between two nodes

  /// having minimum total length.

  ///

  /// \ref lemon::Suurballe "Suurballe" implements an algorithm for

  /// finding arc-disjoint paths having minimum total length (cost)

  /// from a given source node to a given target node in a digraph.

  ///

  /// In fact, this implementation is the specialization of the

  /// \ref CapacityScaling "successive shortest path" algorithm.

  ///

  /// \tparam Digraph The digraph type the algorithm runs on.

  /// The default value is \c ListDigraph.

  /// \tparam LengthMap The type of the length (cost) map.

  /// The default value is <tt>Digraph::ArcMap<int></tt>.

  ///

  /// \warning Length values should be \e non-negative \e integers.

  ///

  /// \note For finding node-disjoint paths this algorithm can be used

#ifdef DOXYGEN

  template <typename Digraph, typename LengthMap>

#else

  template < typename Digraph = ListDigraph,

             typename LengthMap = typename Digraph::template ArcMap<int> >

#endif

  class Suurballe

  {

    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);

    typedef typename LengthMap::Value Length;

    typedef ConstMap<Arc, int> ConstArcMap;

    typedef typename Digraph::template NodeMap<Arc> PredMap;

  public:

    /// The type of the flow map.

    typedef typename Digraph::template ArcMap<int> FlowMap;

    /// The type of the potential map.

    typedef typename Digraph::template NodeMap<Length> PotentialMap;

    /// The type of the path structures.

    typedef SimplePath<Digraph> Path;

  private:

    /// \brief Special implementation of the Dijkstra algorithm

    /// for finding shortest paths in the residual network.

    ///

    /// \ref ResidualDijkstra is a special implementation of the

    /// \ref Dijkstra algorithm for finding shortest paths in the

    /// residual network of the digraph with respect to the reduced arc

    /// lengths and modifying the node potentials according to the

    /// distance of the nodes.

    class ResidualDijkstra

    {

      typedef typename Digraph::template NodeMap<int> HeapCrossRef;

      typedef BinHeap<Length, HeapCrossRef> Heap;

    private:

      // The digraph the algorithm runs on

      const Digraph &_graph;

      // The main maps

      const FlowMap &_flow;

      const LengthMap &_length;

      PotentialMap &_potential;

      // The distance map

      PotentialMap _dist;

      // The pred arc map

      PredMap &_pred;

      // The processed (i.e. permanently labeled) nodes

      std::vector<Node> _proc_nodes;

      Node _s;

      Node _t;

    public:

      /// Constructor.

      ResidualDijkstra( const Digraph &digraph,

                        const FlowMap &flow,

                        const LengthMap &length,

                        PotentialMap &potential,

                        PredMap &pred,

                        Node s, Node t ) :

        _graph(digraph), _flow(flow), _length(length), _potential(potential),

        _dist(digraph), _pred(pred), _s(s), _t(t) {}

      /// \brief Run the algorithm. It returns \c true if a path is found

      /// from the source node to the target node.

      bool run() {

        HeapCrossRef heap_cross_ref(_graph, Heap::PRE_HEAP);

        Heap heap(heap_cross_ref);

        heap.push(_s, 0);

        _pred[_s] = INVALID;

        _proc_nodes.clear();

        // Process nodes

        while (!heap.empty() && heap.top() != _t) {

          Node u = heap.top(), v;

          Length d = heap.prio() + _potential[u], nd;

          _dist[u] = heap.prio();

          heap.pop();

          _proc_nodes.push_back(u);

          // Traverse outgoing arcs

          for (OutArcIt e(_graph, u); e != INVALID; ++e) {

            if (_flow[e] == 0) {

              v = _graph.target(e);

              switch(heap.state(v)) {

              case Heap::PRE_HEAP:

                heap.push(v, d + _length[e] - _potential[v]);

                _pred[v] = e;

                break;