/* -*- C++ -*-

  * lemon/min_cost_flow.h - Part of LEMON, a generic C++ optimization library

  *

  * Copyright (C) 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_MIN_COST_FLOW_H

 #define LEMON_MIN_COST_FLOW_H

 ///\ingroup flowalgs

 ///\file

 ///\brief An algorithm for finding a flow of value \c k (for small values of \c k) having minimal total cost 

 #include <lemon/dijkstra.h>

 #include <lemon/graph_adaptor.h>

 #include <lemon/maps.h>

 #include <vector>

 namespace lemon {

 /// \addtogroup flowalgs

 /// @{

   ///\brief Implementation of an algorithm for finding a flow of value \c k 

   ///(for small values of \c k) having minimal total cost between 2 nodes 

   /// 

   ///

   /// The class \ref lemon::MinCostFlow "MinCostFlow" implements an

   /// algorithm for finding a flow of value \c k having minimal total

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

   /// directed graph with a cost function on the edges. To

   /// this end, the edge-capacities and edge-costs have to be

   /// nonnegative.  The edge-capacities should be integers, but the

   /// edge-costs can be integers, reals or of other comparable

   /// numeric type.  This algorithm is intended to be used only for

   /// small values of \c k, since it is only polynomial in k, not in

   /// the length of k (which is log k): in order to find the minimum

   /// cost flow of value \c k it finds the minimum cost flow of value

   /// \c i for every \c i between 0 and \c k.

   ///

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

   ///\param LengthMap The type of the length map.

   ///\param CapacityMap The capacity map type.

   ///

   ///\author Attila Bernath

   template <typename Graph, typename LengthMap, typename CapacityMap>

   class MinCostFlow {

     typedef typename LengthMap::Value Length;

     //Warning: this should be integer type

     typedef typename CapacityMap::Value Capacity;

     typedef typename Graph::Node Node;

     typedef typename Graph::NodeIt NodeIt;

     typedef typename Graph::Edge Edge;

     typedef typename Graph::OutEdgeIt OutEdgeIt;

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

     typedef ResGraphAdaptor<const Graph,int,CapacityMap,EdgeIntMap> ResGW;

     typedef typename ResGW::Edge ResGraphEdge;

   protected:

     const Graph& g;

     const LengthMap& length;

     const CapacityMap& capacity;

     EdgeIntMap flow; 

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

     PotentialMap potential;

     Node s;

     Node t;

     Length total_length;

     class ModLengthMap {   

       typedef typename Graph::template NodeMap<Length> NodeMap;

       const ResGW& g;

       const LengthMap &length;

       const NodeMap &pot;

     public :

       typedef typename LengthMap::Key Key;

       typedef typename LengthMap::Value Value;

       ModLengthMap(const ResGW& _g, 

 		   const LengthMap &_length, const NodeMap &_pot) : 

 	g(_g), /*rev(_rev),*/ length(_length), pot(_pot) { }

       Value operator[](typename ResGW::Edge e) const {     

 	if (g.forward(e))

 	  return  length[e]-(pot[g.target(e)]-pot[g.source(e)]);   

 	else

 	  return -length[e]-(pot[g.target(e)]-pot[g.source(e)]);   

       }     

     }; //ModLengthMap

     ResGW res_graph;

     ModLengthMap mod_length;

     Dijkstra<ResGW, ModLengthMap> dijkstra;

   public :

     /*! \brief The constructor of the class.

     \param _g The directed graph the algorithm runs on. 

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

     \param _cap The capacity of the edges. 

     \param _s Source node.

     \param _t Target node.

     */

     MinCostFlow(Graph& _g, LengthMap& _length, CapacityMap& _cap, 

 		Node _s, Node _t) : 

       g(_g), length(_length), capacity(_cap), flow(_g), potential(_g), 

       s(_s), t(_t), 

       res_graph(g, capacity, flow), 

       mod_length(res_graph, length, potential),

       dijkstra(res_graph, mod_length) { 

       reset();

       }

     /*! Tries to augment the flow between s and t by 1.

       The return value shows if the augmentation is successful.

      */

     bool augment() {

       dijkstra.run(s);

       if (!dijkstra.reached(t)) {

 	//Unsuccessful augmentation.

 	return false;

       } else {

 	//We have to change the potential

 	for(typename ResGW::NodeIt n(res_graph); n!=INVALID; ++n)

 	  potential.set(n, potential[n]+dijkstra.distMap()[n]);

 	//Augmenting on the shortest path

 	Node n=t;

 	ResGraphEdge e;

 	while (n!=s){

 	  e = dijkstra.predEdge(n);

 	  n = dijkstra.predNode(n);

 	  res_graph.augment(e,1);

 	  //Let's update the total length

 	  if (res_graph.forward(e))

 	    total_length += length[e];

 	  else 

 	    total_length -= length[e];	    

 	}

 	return true;

       }

     }

     /*! \brief Runs the algorithm.

     Runs the algorithm.

     Returns k if there is a flow of value at least k from s to t.

     Otherwise it returns the maximum value of a flow from s to t.

     \param k The value of the flow we are looking for.

     \todo May be it does make sense to be able to start with a nonzero 

     feasible primal-dual solution pair as well.

     \todo If the actual flow value is bigger than k, then everything is 

     cleared and the algorithm starts from zero flow. Is it a good approach?

     */

     int run(int k) {

       if (flowValue()>k) reset();

       while (flowValue()<k && augment()) { }

       return flowValue();

     }

     /*! \brief The class is reset to zero flow and potential.

       The class is reset to zero flow and potential.

      */

     void reset() {

       total_length=0;

       for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);

       for (typename Graph::NodeIt n(g); n!=INVALID; ++n) potential.set(n, 0);  

     }

     /*! Returns the value of the actual flow. 

      */

     int flowValue() const {

       int i=0;

       for (typename Graph::OutEdgeIt e(g, s); e!=INVALID; ++e) i+=flow[e];

       for (typename Graph::InEdgeIt e(g, s); e!=INVALID; ++e) i-=flow[e];

       return i;

     }

     /// Total cost of the found flow.

     /// This function gives back the total cost of the found flow.

     Length totalLength(){

       return total_length;

     }

     ///Returns a const reference to the EdgeMap \c flow. 

     ///Returns a const reference to the EdgeMap \c flow. 

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

     /*! \brief Returns a const reference to the NodeMap \c potential (the dual solution).

     Returns a const reference to the NodeMap \c potential (the dual solution).

     */

     const PotentialMap &getPotential() const { return potential;}

     /*! \brief Checking the complementary slackness optimality criteria.

     This function checks, whether the given flow and potential 

     satisfy the complementary slackness conditions (i.e. these are optimal).

     This function only checks optimality, doesn't bother with feasibility.

     For testing purpose.

     */

     bool checkComplementarySlackness(){

       Length mod_pot;

       Length fl_e;

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

 	//C^{\Pi}_{i,j}

 	mod_pot = length[e]-potential[g.target(e)]+potential[g.source(e)];

 	fl_e = flow[e];

 	if (0<fl_e && fl_e<capacity[e]) {

 	  /// \todo better comparison is needed for real types, moreover, 

 	  /// this comparison here is superfluous.

 	  if (mod_pot != 0)

 	    return false;

 	} 

 	else {

 	  if (mod_pot > 0 && fl_e != 0)

 	    return false;

 	  if (mod_pot < 0 && fl_e != capacity[e])

 	    return false;

 	}

       }

       return true;

     }

   }; //class MinCostFlow

   ///@}

 } //namespace lemon

 #endif //LEMON_MIN_COST_FLOW_H