/* -*- 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_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 two nodes

   /// 

   ///

   /// The \ref lemon::SspMinCostFlow "Successive Shortest Path Minimum

   /// Cost Flow" 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 SspMinCostFlow {

     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.

     SspMinCostFlow(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();

     }

     /// \brief 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);  

     }

     /// \brief 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;

     }

     /// \brief Total cost of the found flow.

     ///

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

     Length totalLength(){

       return total_length;

     }

     /// \brief 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 SspMinCostFlow

   ///@}

 } //namespace lemon

 #endif //LEMON_MIN_COST_FLOW_H