/* -*- C++ -*-

 *

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

 *

 * Copyright (C) 2003-2007

 * 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/graph_adaptor.h>

#include <lemon/tolerance.h>

#include <lemon/bfs.h>

namespace lemon {

  /// \ingroup max_flow

  /// \brief Edmonds-Karp algorithms class.

  ///

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

  /// algorithm producing a flow of maximum value in a directed

  /// graph. 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 edges and the \e starting \e flow value of the

  /// edges should be passed to the algorithm through the

  /// constructor.

  ///

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

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

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

  ///

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

  /// \param _Number The number type of the capacities and the flow values.

  /// \param _CapacityMap The capacity map type.

  /// \param _FlowMap The flow map type.

  /// \param _Tolerance The tolerance class to handle computation problems.

  ///

  /// \author Balazs Dezso 

#ifdef DOXYGEN

  template <typename _Graph, typename _Number,

	    typename _CapacityMap, typename _FlowMap, typename _Tolerance>

#else

  template <typename _Graph, typename _Number,

	    typename _CapacityMap = typename _Graph::template EdgeMap<_Number>,

            typename _FlowMap = typename _Graph::template EdgeMap<_Number>,

	    typename _Tolerance = Tolerance<_Number> >

#endif

  class EdmondsKarp {

  public:

    /// \brief \ref Exception for the case when the source equals the target.

    ///

    /// \ref Exception for the case when the source equals the target.

    ///

    class InvalidArgument : public lemon::LogicError {

    public:

      virtual const char* what() const throw() {

	return "lemon::EdmondsKarp::InvalidArgument";

      }

    };

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

    typedef _Graph Graph;

    /// \brief The value type of the algorithms.

    typedef _Number Number;

    /// \brief The capacity map on the edges.

    typedef _CapacityMap CapacityMap;

    /// \brief The flow map on the edges.

    typedef _FlowMap FlowMap;

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

    typedef _Tolerance Tolerance;

    typedef ResGraphAdaptor<Graph, Number, CapacityMap, 

                            FlowMap, Tolerance> ResGraph;

  private:

    typedef typename Graph::Node Node;

    typedef typename Graph::Edge Edge;

    typedef typename Graph::NodeIt NodeIt;

    typedef typename Graph::EdgeIt EdgeIt;

    typedef typename Graph::InEdgeIt InEdgeIt;

    typedef typename Graph::OutEdgeIt OutEdgeIt;

  public:

    /// \brief The constructor of the class.

    ///

    /// The constructor of the class. 

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

    /// \param source The source node.

    /// \param target The target node.

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

    /// \param flow The flow of the edges. 

    /// \param tolerance Tolerance class.

    EdmondsKarp(const Graph& graph, Node source, Node target,

                const CapacityMap& capacity, FlowMap& flow, 

                const Tolerance& tolerance = Tolerance())

      : _graph(graph), _capacity(capacity), _flow(flow), 

        _tolerance(tolerance), _resgraph(graph, capacity, flow, tolerance), 

        _source(source), _target(target) 

    {

      if (_source == _target) {

        throw InvalidArgument();

      }

    }

    /// \brief Initializes the algorithm

    /// 

    /// It sets the flow to empty flow.

    void init() {

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

        _flow.set(it, 0);

      }

      _value = 0;

    }

    /// \brief Initializes the algorithm

    /// 

    /// If the flow map initially flow this let the flow map

    /// unchanged but the flow value will be set by the flow

    /// on the outedges from the source.

    void flowInit() {

      _value = 0;

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

        _value += _flow[jt];

      }

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

        _value -= _flow[jt];

      }

    }

    /// \brief Initializes the algorithm

    /// 

    /// If the flow map initially flow this let the flow map

    /// unchanged but the flow value will be set by the flow

    /// on the outedges from the source. It also checks that

    /// the flow map really contains a flow.

    /// \return %True when the flow map really a flow.

    bool checkedFlowInit() {

      _value = 0;

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

        _value += _flow[jt];

      }

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

        _value -= _flow[jt];

      }

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

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

        Number outFlow = 0;

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

          outFlow += _flow[jt];

        }

        Number inFlow = 0;

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

          inFlow += _flow[jt];

        }

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

          return false;

        }

      }

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

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

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

      }

      return true;

    }

    /// \brief Augment the solution on an edge shortest path.

    /// 

    /// Augment the solution on an edge shortest path. It search an

    /// edge shortest path between the source and the target

    /// in the residual graph with the bfs algoritm.

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

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

    /// false.

    /// \return %False when the augmenting is not success so the

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

    bool augment() {

      typename Bfs<ResGraph>

      ::template DefDistMap<NullMap<Node, int> >

      ::Create bfs(_resgraph);

      NullMap<Node, int> distMap;

      bfs.distMap(distMap);

      bfs.init();

      bfs.addSource(_source);

      bfs.start(_target);

      if (!bfs.reached(_target)) {

        return false;

      }

      Number min = _resgraph.rescap(bfs.predEdge(_target));

      for (Node it = bfs.predNode(_target); it != _source; 

           it = bfs.predNode(it)) {

        if (min > _resgraph.rescap(bfs.predEdge(it))) {

          min = _resgraph.rescap(bfs.predEdge(it));

        }

      } 

      for (Node it = _target; it != _source; it = bfs.predNode(it)) {

        _resgraph.augment(bfs.predEdge(it), min);

      }

      _value += min;

      return true;

    }

    /// \brief Executes the algorithm

    ///

    /// It runs augmenting phases until the optimal solution is reached. 

    void start() {

      while (augment()) {}

    }

    /// \brief Gives back the current flow value.

    ///

    /// Gives back the current flow _value.

    Number flowValue() const {

      return _value;

    }

    /// \brief runs the algorithm.

    /// 

    /// It is just a shorthand for:

    ///

    ///\code 

    /// ek.init();

    /// ek.start();

    ///\endcode

    void run() {

      init();

      start();

    }

    /// \brief Returns a minimum value cut.

    ///

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

    /// It simply calls the minMinCut member.

    /// \retval cut Write node bool map. 

    template <typename CutMap>

    void minCut(CutMap& cut) const {

      minMinCut(cut);

    }

    /// \brief Returns the inclusionwise minimum of the minimum value cuts.

    ///

    /// Sets \c cut to the characteristic vector of the minimum value cut

    /// which is inclusionwise minimum. It is computed by processing a

    /// bfs from the source node \c source in the residual graph.  

    /// \retval cut Write node bool map. 

    template <typename CutMap>

    void minMinCut(CutMap& cut) const {

      typename Bfs<ResGraph>

      ::template DefDistMap<NullMap<Node, int> >

      ::template DefProcessedMap<CutMap>

      ::Create bfs(_resgraph);

      NullMap<Node, int> distMap;

      bfs.distMap(distMap);

      bfs.processedMap(cut);

      bfs.run(_source);

    }

    /// \brief Returns the inclusionwise minimum of the minimum value cuts.

    ///

    /// Sets \c cut to the characteristic vector of the minimum value cut

    /// which is inclusionwise minimum. It is computed by processing a

    /// bfs from the source node \c source in the residual graph.  

    /// \retval cut Write node bool map. 

    template <typename CutMap>

    void maxMinCut(CutMap& cut) const {

      typedef RevGraphAdaptor<const ResGraph> RevGraph;

      RevGraph revgraph(_resgraph);

      typename Bfs<RevGraph>

      ::template DefDistMap<NullMap<Node, int> >

      ::template DefPredMap<NullMap<Node, Edge> >

      ::template DefProcessedMap<NotWriteMap<CutMap> >

      ::Create bfs(revgraph);

      NullMap<Node, int> distMap;

      bfs.distMap(distMap);

      NullMap<Node, Edge> predMap;

      bfs.predMap(predMap);

      NotWriteMap<CutMap> notcut(cut);

      bfs.processedMap(notcut);

      bfs.run(_target);

    }

    /// \brief Returns the source node.

    ///

    /// Returns the source node.

    /// 

    Node source() const { 

      return _source;

    }

    /// \brief Returns the target node.

    ///

    /// Returns the target node.

    /// 

    Node target() const { 

      return _target;

    }

    /// \brief Returns a reference to capacity map.

    ///

    /// Returns a reference to capacity map.

    /// 

    const CapacityMap &capacityMap() const { 

      return *_capacity;

    }

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

    ///

    /// Returns a reference to flow map.

    /// 

    const FlowMap &flowMap() const { 

      return *_flow;

    }

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

    ///

    /// Returns the tolerance used by algorithm.

    const Tolerance& tolerance() const {

      return tolerance;

    } 

  private:

    const Graph& _graph;

    const CapacityMap& _capacity;

    FlowMap& _flow;

    Tolerance _tolerance;

    ResGraph _resgraph;

    Node _source, _target;

    Number _value;

  };

}

#endif