/* -*- C++ -*-
 *
 * This file is a part of LEMON, a generic C++ optimization library
 *
 * Copyright (C) 2003-2006
 * 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 flowalgs
/// \brief Implementation of the Edmonds-Karp algorithm.

#include <lemon/graph_adaptor.h>
#include <lemon/tolerance.h>
#include <lemon/bfs.h>

namespace lemon {

  /// \ingroup flowalgs
  /// \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