///\brief An algorithm for finding k paths of minimal total length.

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

 /// nodes having minimum total length.

 #include <lemon/maps.h>

 /// \addtogroup shortest_path

   /// \addtogroup shortest_path

 /// @{

   /// @{

   ///\brief Implementation of an algorithm for finding k edge-disjoint

   /// \brief Implementation of an algorithm for finding edge-disjoint

   /// paths between 2 nodes of minimal total length

   /// paths between two nodes having minimum total length.

   ///

   ///

   /// The class \ref lemon::Suurballe implements

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

   /// an algorithm for finding k edge-disjoint paths

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

   /// from a given source node to a given target node in an

   /// from a given source node to a given target node in a directed

   /// edge-weighted directed graph having minimal total weight (length).

   /// graph.

   ///

   ///

   ///\warning Length values should be nonnegative!

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

   /// 

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

   ///\param Graph The directed graph type the algorithm runs on.

   ///

   ///\param LengthMap The type of the length map (values should be nonnegative).

   /// \tparam Graph The directed graph type the algorithm runs on.

   ///

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

   ///\note It it questionable whether it is correct to call this method after

   ///

   ///%Suurballe for it is just a special case of Edmonds' and Karp's algorithm

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

   ///for finding minimum cost flows. In fact, this implementation just

   ///

   ///wraps the SspMinCostFlow algorithms. The paper of both %Suurballe and

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

   ///Edmonds-Karp published in 1972, therefore it is possibly right to

   /// with \ref SplitGraphAdaptor.

   ///state that they are

   ///

   ///independent results. Most frequently this special case is referred as

   /// \author Attila Bernath and Peter Kovacs

   ///%Suurballe method in the literature, especially in communication

   ///network context.

   template < typename Graph, 

   ///\author Attila Bernath

              typename LengthMap = typename Graph::template EdgeMap<int> >

   template <typename Graph, typename LengthMap>

   class Suurballe

   class Suurballe{

   {

     GRAPH_TYPEDEFS(typename Graph);

     typedef typename Graph::Node Node;

     // Edge map of the current flow

     typedef typename Graph::NodeIt NodeIt;

     FlowMap *_flow;

     typedef typename Graph::Edge Edge;

     bool _local_flow;

     typedef typename Graph::OutEdgeIt OutEdgeIt;

     // Node map of the current potentials

     typedef typename Graph::template EdgeMap<int> EdgeIntMap;

     PotentialMap *_potential;

     bool _local_potential;

     typedef ConstMap<Edge,int> ConstMap;

     // The source node

     const Graph& G;

     Node _source;

     // The target node

     Node s;

     Node _target;

     Node t;

     // Container to store the found paths

     //Auxiliary variables

     std::vector< SimplePath<Graph> > paths;

     //This is the capacity map for the mincostflow problem

     int _path_num;

     ConstMap const1map;

     //This MinCostFlow instance will actually solve the problem

     // The pred edge map

     SspMinCostFlow<Graph, LengthMap, ConstMap> min_cost_flow;

     PredMap _pred;

     // Implementation of the Dijkstra algorithm for finding augmenting

     //Container to store found paths

     // shortest paths in the residual network

     std::vector<SimplePath<Graph> > paths;

     ResidualDijkstra *_dijkstra;

   public :

   public:

     /// \brief Constructor.

     /// \brief The constructor of the class.

     ///

     ///

     /// Constructor.

     /// \param _G The directed graph the algorithm runs on. 

     ///

     /// \param _length The length (weight or cost) of the edges. 

     /// \param graph The directed graph the algorithm runs on.

     /// \param _s Source node.

     /// \param length The length (cost) values of the edges.

     /// \param _t Target node.

     /// \param s The source node.

     Suurballe(Graph& _G, LengthMap& _length, Node _s, Node _t) : 

     /// \param t The target node.

       G(_G), s(_s), t(_t), const1map(1), 

     Suurballe( const Graph &graph,

       min_cost_flow(_G, _length, const1map, _s, _t) { }

                const LengthMap &length,

                Node s, Node t ) :

       _graph(graph), _length(length), _flow(0), _local_flow(false),

       _potential(0), _local_potential(false), _source(s), _target(t),

       _pred(graph) {}

     /// Destructor.

     ~Suurballe() {

       if (_local_flow) delete _flow;

       if (_local_potential) delete _potential;

       delete _dijkstra;

     }

     /// \brief Sets the flow map.

     ///

     /// Sets the flow map.

     ///

     /// The found flow contains only 0 and 1 values. It is the union of

     /// the found edge-disjoint paths.

     ///

     /// \return \c (*this)

     Suurballe& flowMap(FlowMap &map) {

       if (_local_flow) {

         delete _flow;

         _local_flow = false;

       }

       _flow = ↦

       return *this;

     }

     /// \brief Sets the potential map.

     ///

     /// Sets the potential map.

     ///

     /// The potentials provide the dual solution of the underlying 

     /// minimum cost flow problem.

     ///

     /// \return \c (*this)

     Suurballe& potentialMap(PotentialMap &map) {

       if (_local_potential) {

         delete _potential;

         _local_potential = false;

       }

       _potential = ↦

       return *this;

     }

     /// \name Execution control

     /// The simplest way to execute the algorithm is to call the run()

     /// function.

     /// \n

     /// If you only need the flow that is the union of the found

     /// edge-disjoint paths, you may call init() and findFlow().

     /// @{

     /// Returns k if there are at least k edge-disjoint paths from s to t.

     ///

     /// Otherwise it returns the number of edge-disjoint paths found 

     /// \param k The number of paths to be found.

     /// from s to t.

     ///

     ///

     /// \return \c k if there are at least \c k edge-disjoint paths

     /// \param k How many paths are we looking for?

     /// from \c s to \c t. Otherwise it returns the number of

     ///

     /// edge-disjoint paths found.

     int run(int k) {

     ///

       int i = min_cost_flow.run(k);

     /// \note Apart from the return value, <tt>s.run(k)</tt> is just a

     /// shortcut of the following code.

       //Let's find the paths

     /// \code

       //We put the paths into stl vectors (as an inner representation). 

     ///   s.init();

       //In the meantime we lose the information stored in 'reversed'.

     ///   s.findFlow(k);

       //We suppose the lengths to be positive now.

     ///   s.findPaths();

     /// \endcode

       //We don't want to change the flow of min_cost_flow, so we make a copy

     int run(int k = 2) {

       //The name here suggests that the flow has only 0/1 values.

       init();

       EdgeIntMap reversed(G); 

       findFlow(k);

       findPaths();

       for(typename Graph::EdgeIt e(G); e!=INVALID; ++e) 

       return _path_num;

 	reversed[e] = min_cost_flow.getFlow()[e];

     }

     /// \brief Initializes the algorithm.

     ///

     /// Initializes the algorithm.

     void init() {

       // Initializing maps

       if (!_flow) {

         _flow = new FlowMap(_graph);

         _local_flow = true;

       }

       if (!_potential) {

         _potential = new PotentialMap(_graph);

         _local_potential = true;

       }

       for (EdgeIt e(_graph); e != INVALID; ++e) (*_flow)[e] = 0;

       for (NodeIt n(_graph); n != INVALID; ++n) (*_potential)[n] = 0;

       _dijkstra = new ResidualDijkstra( _graph, *_flow, _length, 

                                         *_potential, _pred,

                                         _source, _target );

     }

     /// \brief Executes the successive shortest path algorithm to find

     /// an optimal flow.

     ///

     /// Executes the successive shortest path algorithm to find a

     /// minimum cost flow, which is the union of \c k or less

     /// edge-disjoint paths.

     ///

     /// \return \c k if there are at least \c k edge-disjoint paths

     /// from \c s to \c t. Otherwise it returns the number of

     /// edge-disjoint paths found.

     ///

     /// \pre \ref init() must be called before using this function.

     int findFlow(int k = 2) {

       // Finding shortest paths

       _path_num = 0;

       while (_path_num < k) {

         // Running Dijkstra

         if (!_dijkstra->run()) break;

         ++_path_num;

         // Setting the flow along the found shortest path

         Node u = _target;

         Edge e;

         while ((e = _pred[u]) != INVALID) {

           if (u == _graph.target(e)) {

             (*_flow)[e] = 1;

             u = _graph.source(e);

           } else {

             (*_flow)[e] = 0;

             u = _graph.target(e);

           }

         }

       }

       return _path_num;

     }

     /// \brief Computes the paths from the flow.

     ///

     /// Computes the paths from the flow.

     ///

     /// \pre \ref init() and \ref findFlow() must be called before using

     /// this function.

     void findPaths() {

       // Creating the residual flow map (the union of the paths not

       // found so far)

       FlowMap res_flow(*_flow);

       paths.resize(k);

       paths.resize(_path_num);

       for (int j=0; j<i; ++j){

       for (int i = 0; i < _path_num; ++i) {

 	Node n=s;

         Node n = _source;

         while (n != _target) {

 	while (n!=t){

           OutEdgeIt e(_graph, n);

           for ( ; res_flow[e] == 0; ++e) ;

 	  OutEdgeIt e(G, n);

           n = _graph.target(e);

           paths[i].addBack(e);

 	  while (!reversed[e]){

           res_flow[e] = 0;

 	    ++e;

         }

 	  }

       }

 	  n = G.target(e);

     }

 	  paths[j].addBack(e);

 	  reversed[e] = 1-reversed[e];

     /// @}

 	}

     /// \name Query Functions

       }

     /// The result of the algorithm can be obtained using these

       return i;

     /// functions.

     }

     /// \n The algorithm should be executed before using them.

     /// @{

     /// \brief Returns the total length of the paths.

     ///

     /// \brief Returns a const reference to the edge map storing the

     /// This function gives back the total length of the found paths.

     /// found flow.

     Length totalLength(){

     ///

       return min_cost_flow.totalLength();

     /// Returns a const reference to the edge map storing the flow that

     }

     /// is the union of the found edge-disjoint paths.

     ///

     /// \brief Returns the found flow.

     /// \pre \ref run() or findFlow() must be called before using this

     ///

     /// function.

     /// This function returns a const reference to the EdgeMap \c flow.

     const FlowMap& flowMap() const {

     const EdgeIntMap &getFlow() const { return min_cost_flow.flow;}

       return *_flow;

     }

     /// \brief Returns the optimal dual solution

     ///

     /// \brief Returns a const reference to the node map storing the

     /// This function returns a const reference to the NodeMap \c

     /// found potentials (the dual solution).

     /// potential (the dual solution).

     ///

     const EdgeIntMap &getPotential() const { return min_cost_flow.potential;}

     /// Returns a const reference to the node map storing the found

     /// potentials that provide the dual solution of the underlying 

     /// \brief Checks whether the complementary slackness holds.

     /// minimum cost flow problem.

     ///

     ///

     /// This function checks, whether the given solution is optimal.

     /// \pre \ref run() or findFlow() must be called before using this

     /// Currently this function only checks optimality, doesn't bother

     /// function.

     /// with feasibility.  It is meant for testing purposes.

     const PotentialMap& potentialMap() const {

     bool checkComplementarySlackness(){

       return *_potential;

       return min_cost_flow.checkComplementarySlackness();

     }

     }

     /// \brief Returns the flow on the given edge.

     typedef SimplePath<Graph> Path; 

     ///

     /// Returns the flow on the given edge.

     /// \brief Read the found paths.

     /// It is \c 1 if the edge is involved in one of the found paths,

     ///

     /// otherwise it is \c 0.

     /// This function gives back the \c j-th path in argument p.

     ///

     /// Assumes that \c run() has been run and nothing has changed

     /// \pre \ref run() or findFlow() must be called before using this

     /// since then.

     /// function.

     ///

     int flow(const Edge& edge) const {

     /// \warning It is assumed that \c p is constructed to be a path

       return (*_flow)[edge];

     /// of graph \c G.  If \c j is not less than the result of

     }

     /// previous \c run, then the result here will be an empty path

     /// (\c j can be 0 as well).

     /// \brief Returns the potential of the given node.

     ///

     ///

     /// \param j Which path you want to get from the found paths (in a

     /// Returns the potential of the given node.

     /// real application you would get the found paths iteratively).

     ///

     Path path(int j) const {

     /// \pre \ref run() or findFlow() must be called before using this

       return paths[j];

     /// function.

     }

     Length potential(const Node& node) const {

       return (*_potential)[node];

     /// \brief Gives back the number of the paths.

     }

     ///

     /// Gives back the number of the constructed paths.

     /// \brief Returns the total length (cost) of the found paths (flow).

     ///

     /// Returns the total length (cost) of the found paths (flow).

     /// The complexity of the function is \f$ O(e) \f$.

     ///

     /// \pre \ref run() or findFlow() must be called before using this

     /// function.

     Length totalLength() const {

       Length c = 0;

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

         c += (*_flow)[e] * _length[e];

       return c;

     }

     /// \brief Returns the number of the found paths.

     ///

     /// Returns the number of the found paths.

     ///

     /// \pre \ref run() or findFlow() must be called before using this

     /// function.

       return paths.size();

       return _path_num;

     }

     }

     /// \brief Returns a const reference to the specified path.

     ///

     /// Returns a const reference to the specified path.

     ///

     /// \param i The function returns the \c i-th path.

     /// \c i must be between \c 0 and <tt>%pathNum()-1</tt>.

     ///

     /// \pre \ref run() or findPaths() must be called before using this

     /// function.

     Path path(int i) const {

       return paths[i];

     }

     /// @}