/* -*- mode: C++; indent-tabs-mode: nil; -*-

  *

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

  *

  * Copyright (C) 2003-2010

  * 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_EDMONDS_KARP_H

 #define LEMON_EDMONDS_KARP_H

 /// \file

 /// \ingroup max_flow

 /// \brief Implementation of the Edmonds-Karp algorithm.

 #include <lemon/tolerance.h>

 #include <vector>

 namespace lemon {

   /// \brief Default traits class of EdmondsKarp class.

   ///

   /// Default traits class of EdmondsKarp class.

   /// \param GR Digraph type.

   /// \param CAP Type of capacity map.

   template <typename GR, typename CAP>

   struct EdmondsKarpDefaultTraits {

     /// \brief The digraph type the algorithm runs on. 

     typedef GR Digraph;

     /// \brief The type of the map that stores the arc capacities.

     ///

     /// The type of the map that stores the arc capacities.

     /// It must meet the \ref concepts::ReadMap "ReadMap" concept.

     typedef CAP CapacityMap;

     /// \brief The type of the flow values.

     typedef typename CapacityMap::Value Value;

     /// \brief The type of the map that stores the flow values.

     ///

     /// The type of the map that stores the flow values.

     /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.

 #ifdef DOXYGEN

     typedef GR::ArcMap<Value> FlowMap;

 #else

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

 #endif

     /// \brief Instantiates a FlowMap.

     ///

     /// This function instantiates a \ref FlowMap. 

     /// \param digraph The digraph for which we would like to define

     /// the flow map.

     static FlowMap* createFlowMap(const Digraph& digraph) {

       return new FlowMap(digraph);

     }

     /// \brief The tolerance used by the algorithm

     ///

     /// The tolerance used by the algorithm to handle inexact computation.

     typedef lemon::Tolerance<Value> Tolerance;

   };

   /// \ingroup max_flow

   ///

   /// \brief Edmonds-Karp algorithms class.

   ///

   /// This class provides an implementation of the \e Edmonds-Karp \e

   /// algorithm producing a \ref max_flow "flow of maximum value" in a

   /// digraph \ref clrs01algorithms, \ref amo93networkflows,

   /// \ref edmondskarp72theoretical.

   /// The Edmonds-Karp algorithm is slower than the Preflow

   /// algorithm, but it has an advantage of the step-by-step execution

   /// control with feasible flow solutions. The \e source node, the \e

   /// target node, the \e capacity of the arcs and the \e starting \e

   /// flow value of the arcs should be passed to the algorithm

   /// through the constructor.

   ///

   /// The time complexity of the algorithm is \f$ O(nm^2) \f$ in

   /// worst case. Always try the Preflow algorithm instead of this if

   /// you just want to compute the optimal flow.

   ///

   /// \tparam GR The type of the digraph the algorithm runs on.

   /// \tparam CAP The type of the capacity map. The default map

   /// type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>". 

   /// \tparam TR The traits class that defines various types used by the

   /// algorithm. By default, it is \ref EdmondsKarpDefaultTraits

   /// "EdmondsKarpDefaultTraits<GR, CAP>".

   /// In most cases, this parameter should not be set directly,

   /// consider to use the named template parameters instead.

 #ifdef DOXYGEN

   template <typename GR, typename CAP, typename TR>

 #else 

   template <typename GR,

 	    typename CAP = typename GR::template ArcMap<int>,

             typename TR = EdmondsKarpDefaultTraits<GR, CAP> >

 #endif

   class EdmondsKarp {

   public:

     /// The \ref EdmondsKarpDefaultTraits "traits class" of the algorithm.

     typedef TR Traits;

     /// The type of the digraph the algorithm runs on.

     typedef typename Traits::Digraph Digraph;

     /// The type of the capacity map.

     typedef typename Traits::CapacityMap CapacityMap;

     /// The type of the flow values.

     typedef typename Traits::Value Value; 

     /// The type of the flow map.

     typedef typename Traits::FlowMap FlowMap;

     /// The type of the tolerance.

     typedef typename Traits::Tolerance Tolerance;

   private:

     TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);

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

     const Digraph& _graph;

     const CapacityMap* _capacity;

     Node _source, _target;

     FlowMap* _flow;

     bool _local_flow;

     PredMap* _pred;

     std::vector<Node> _queue;

     Tolerance _tolerance;

     Value _flow_value;

     void createStructures() {

       if (!_flow) {

 	_flow = Traits::createFlowMap(_graph);

 	_local_flow = true;

       }

       if (!_pred) {

 	_pred = new PredMap(_graph);

       }

       _queue.resize(countNodes(_graph));

     }

     void destroyStructures() {

       if (_local_flow) {

 	delete _flow;

       }

       if (_pred) {

 	delete _pred;

       }

     }

   public:

     ///\name Named template parameters

     ///@{

     template <typename T>

     struct DefFlowMapTraits : public Traits {

       typedef T FlowMap;

       static FlowMap *createFlowMap(const Digraph&) {

 	LEMON_ASSERT(false, "FlowMap is not initialized");

         return 0;

       }

     };

     /// \brief \ref named-templ-param "Named parameter" for setting

     /// FlowMap type

     ///

     /// \ref named-templ-param "Named parameter" for setting FlowMap

     /// type

     template <typename T>

     struct DefFlowMap 

       : public EdmondsKarp<Digraph, CapacityMap, DefFlowMapTraits<T> > {

       typedef EdmondsKarp<Digraph, CapacityMap, DefFlowMapTraits<T> >

       Create;

     };

     /// @}

   protected:

     EdmondsKarp() {}

   public:

     /// \brief The constructor of the class.

     ///

     /// The constructor of the class. 

     /// \param digraph The digraph the algorithm runs on. 

     /// \param capacity The capacity of the arcs. 

     /// \param source The source node.

     /// \param target The target node.

     EdmondsKarp(const Digraph& digraph, const CapacityMap& capacity,

 		Node source, Node target)

       : _graph(digraph), _capacity(&capacity), _source(source), _target(target),

 	_flow(0), _local_flow(false), _pred(0), _tolerance(), _flow_value()

     {

       LEMON_ASSERT(_source != _target,

                    "Flow source and target are the same nodes.");

     }

     /// \brief Destructor.

     ///

     /// Destructor.

     ~EdmondsKarp() {

       destroyStructures();

     }

     /// \brief Sets the capacity map.

     ///

     /// Sets the capacity map.

     /// \return <tt>(*this)</tt>

     EdmondsKarp& capacityMap(const CapacityMap& map) {

       _capacity = ↦

       return *this;

     }

     /// \brief Sets the flow map.

     ///

     /// Sets the flow map.

     /// If you don't use this function before calling \ref run() or

     /// \ref init(), an instance will be allocated automatically.

     /// The destructor deallocates this automatically allocated map,

     /// of course.

     /// \return <tt>(*this)</tt>

     EdmondsKarp& flowMap(FlowMap& map) {

       if (_local_flow) {

 	delete _flow;

 	_local_flow = false;

       }

       _flow = ↦

       return *this;

     }

     /// \brief Sets the source node.

     ///

     /// Sets the source node.

     /// \return <tt>(*this)</tt>

     EdmondsKarp& source(const Node& node) {

       _source = node;

       return *this;

     }

     /// \brief Sets the target node.

     ///

     /// Sets the target node.

     /// \return <tt>(*this)</tt>

     EdmondsKarp& target(const Node& node) {

       _target = node;

       return *this;

     }

     /// \brief Sets the tolerance used by algorithm.

     ///

     /// Sets the tolerance used by algorithm.

     /// \return <tt>(*this)</tt>

     EdmondsKarp& tolerance(const Tolerance& tolerance) {

       _tolerance = tolerance;

       return *this;

     } 

     /// \brief Returns a const reference to the tolerance.

     ///

     /// Returns a const reference to the tolerance object used by

     /// the algorithm.

     const Tolerance& tolerance() const {

       return _tolerance;

     } 

     /// \name Execution control

     /// The simplest way to execute the algorithm is to use \ref run().\n

     /// If you need better control on the initial solution or the execution,

     /// you have to call one of the \ref init() functions first, then

     /// \ref start() or multiple times the \ref augment() function.

     ///@{

     /// \brief Initializes the algorithm.

     ///

     /// Initializes the internal data structures and sets the initial

     /// flow to zero on each arc.

     void init() {

       createStructures();

       for (ArcIt it(_graph); it != INVALID; ++it) {

         _flow->set(it, 0);

       }

       _flow_value = 0;

     }

     /// \brief Initializes the algorithm using the given flow map.

     ///

     /// Initializes the internal data structures and sets the initial

     /// flow to the given \c flowMap. The \c flowMap should

     /// contain a feasible flow, i.e. at each node excluding the source

     /// and the target, the incoming flow should be equal to the

     /// outgoing flow.

     template <typename FlowMap>

     void flowInit(const FlowMap& flowMap) {

       createStructures();

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

 	_flow->set(e, flowMap[e]);

       }

       _flow_value = 0;

       for (OutArcIt jt(_graph, _source); jt != INVALID; ++jt) {

         _flow_value += (*_flow)[jt];

       }

       for (InArcIt jt(_graph, _source); jt != INVALID; ++jt) {

         _flow_value -= (*_flow)[jt];

       }

     }

     /// \brief Initializes the algorithm using the given flow map.

     ///

     /// Initializes the internal data structures and sets the initial

     /// flow to the given \c flowMap. The \c flowMap should

     /// contain a feasible flow, i.e. at each node excluding the source

     /// and the target, the incoming flow should be equal to the

     /// outgoing flow. 

     /// \return \c false when the given \c flowMap does not contain a

     /// feasible flow.

     template <typename FlowMap>

     bool checkedFlowInit(const FlowMap& flowMap) {

       createStructures();

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

 	_flow->set(e, flowMap[e]);

       }

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

         if (it == _source || it == _target) continue;

         Value outFlow = 0;

         for (OutArcIt jt(_graph, it); jt != INVALID; ++jt) {

           outFlow += (*_flow)[jt];

         }

         Value inFlow = 0;

         for (InArcIt jt(_graph, it); jt != INVALID; ++jt) {

           inFlow += (*_flow)[jt];

         }

         if (_tolerance.different(outFlow, inFlow)) {

           return false;

         }

       }

       for (ArcIt it(_graph); it != INVALID; ++it) {

         if (_tolerance.less((*_flow)[it], 0)) return false;

         if (_tolerance.less((*_capacity)[it], (*_flow)[it])) return false;

       }

       _flow_value = 0;

       for (OutArcIt jt(_graph, _source); jt != INVALID; ++jt) {

         _flow_value += (*_flow)[jt];

       }

       for (InArcIt jt(_graph, _source); jt != INVALID; ++jt) {

         _flow_value -= (*_flow)[jt];

       }

       return true;

     }

     /// \brief Augments the solution along a shortest path.

     /// 

     /// Augments the solution along a shortest path. This function searches a

     /// shortest path between the source and the target

     /// in the residual digraph by the Bfs algoritm.

     /// Then it increases the flow on this path with the minimal residual

     /// capacity on the path. If there is no such path, it gives back

     /// false.

     /// \return \c false when the augmenting did not success, i.e. the

     /// current flow is a feasible and optimal solution.

     bool augment() {

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

 	_pred->set(n, INVALID);

       }

       int first = 0, last = 1;

       _queue[0] = _source;

       _pred->set(_source, OutArcIt(_graph, _source));

       while (first != last && (*_pred)[_target] == INVALID) {

 	Node n = _queue[first++];

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

 	  Value rem = (*_capacity)[e] - (*_flow)[e];

 	  Node t = _graph.target(e);

 	  if (_tolerance.positive(rem) && (*_pred)[t] == INVALID) {

 	    _pred->set(t, e);

 	    _queue[last++] = t;

 	  }

 	}

 	for (InArcIt e(_graph, n); e != INVALID; ++e) {

 	  Value rem = (*_flow)[e];

 	  Node t = _graph.source(e);

 	  if (_tolerance.positive(rem) && (*_pred)[t] == INVALID) {

 	    _pred->set(t, e);

 	    _queue[last++] = t;

 	  }

 	}

       }

       if ((*_pred)[_target] != INVALID) {

 	Node n = _target;

 	Arc e = (*_pred)[n];

 	Value prem = (*_capacity)[e] - (*_flow)[e];

 	n = _graph.source(e);

 	while (n != _source) {

 	  e = (*_pred)[n];

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

 	    Value rem = (*_capacity)[e] - (*_flow)[e];

 	    if (rem < prem) prem = rem;

 	    n = _graph.source(e);

 	  } else {

 	    Value rem = (*_flow)[e];

 	    if (rem < prem) prem = rem;

 	    n = _graph.target(e);   

 	  } 

 	}

 	n = _target;

 	e = (*_pred)[n];

 	_flow->set(e, (*_flow)[e] + prem);

 	n = _graph.source(e);

 	while (n != _source) {

 	  e = (*_pred)[n];

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

 	    _flow->set(e, (*_flow)[e] + prem);

 	    n = _graph.source(e);

 	  } else {

 	    _flow->set(e, (*_flow)[e] - prem);

 	    n = _graph.target(e);   

 	  } 

 	}

 	_flow_value += prem;	

 	return true;

       } else {

 	return false;

       }

     }

     /// \brief Executes the algorithm

     ///

     /// Executes the algorithm by performing augmenting phases until the

     /// optimal solution is reached. 

     /// \pre One of the \ref init() functions must be called before

     /// using this function.

     void start() {

       while (augment()) {}

     }

     /// \brief Runs the algorithm.

     /// 

     /// Runs the Edmonds-Karp algorithm.

     /// \note ek.run() is just a shortcut of the following code.

     ///\code 

     /// ek.init();

     /// ek.start();

     ///\endcode

     void run() {

       init();

       start();

     }

     /// @}

     /// \name Query Functions

     /// The result of the Edmonds-Karp algorithm can be obtained using these

     /// functions.\n

     /// Either \ref run() or \ref start() should be called before using them.

     ///@{

     /// \brief Returns the value of the maximum flow.

     ///

     /// Returns the value of the maximum flow found by the algorithm.

     ///

     /// \pre Either \ref run() or \ref init() must be called before

     /// using this function.

     Value flowValue() const {

       return _flow_value;

     }

     /// \brief Returns the flow value on the given arc.

     ///

     /// Returns the flow value on the given arc.

     ///

     /// \pre Either \ref run() or \ref init() must be called before

     /// using this function.

     Value flow(const Arc& arc) const {

       return (*_flow)[arc];

     }

     /// \brief Returns a const reference to the flow map.

     ///

     /// Returns a const reference to the arc map storing the found flow.

     ///

     /// \pre Either \ref run() or \ref init() must be called before

     /// using this function.

     const FlowMap& flowMap() const {

       return *_flow;

     }

     /// \brief Returns \c true when the node is on the source side of the

     /// minimum cut.

     ///

     /// Returns true when the node is on the source side of the found

     /// minimum cut.

     ///

     /// \pre Either \ref run() or \ref init() must be called before

     /// using this function.

     bool minCut(const Node& node) const {

       return ((*_pred)[node] != INVALID) or node == _source;

     }

     /// \brief Gives back a minimum value cut.

     ///

     /// Sets \c cutMap to the characteristic vector of a minimum value

     /// cut. \c cutMap should be a \ref concepts::WriteMap "writable"

     /// node map with \c bool (or convertible) value type.

     ///

     /// \note This function calls \ref minCut() for each node, so it runs in

     /// O(n) time.

     ///

     /// \pre Either \ref run() or \ref init() must be called before

     /// using this function.

     template <typename CutMap>

     void minCutMap(CutMap& cutMap) const {

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

 	cutMap.set(n, (*_pred)[n] != INVALID);

       }

       cutMap.set(_source, true);

     }    

     /// @}

   };

 }

 #endif