alpar@906: /* -*- C++ -*-
ladanyi@1435:  * lemon/min_cost_flow.h - Part of LEMON, a generic C++ optimization library
alpar@906:  *
alpar@1164:  * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
alpar@1359:  * (Egervary Research Group on Combinatorial Optimization, EGRES).
alpar@906:  *
alpar@906:  * Permission to use, modify and distribute this software is granted
alpar@906:  * provided that this copyright notice appears in all copies. For
alpar@906:  * precise terms see the accompanying LICENSE file.
alpar@906:  *
alpar@906:  * This software is provided "AS IS" with no warranty of any kind,
alpar@906:  * express or implied, and with no claim as to its suitability for any
alpar@906:  * purpose.
alpar@906:  *
alpar@906:  */
alpar@906: 
alpar@921: #ifndef LEMON_MIN_COST_FLOW_H
alpar@921: #define LEMON_MIN_COST_FLOW_H
alpar@899: 
alpar@899: ///\ingroup flowalgs
alpar@899: ///\file
alpar@899: ///\brief An algorithm for finding a flow of value \c k (for small values of \c k) having minimal total cost 
alpar@899: 
alpar@899: 
alpar@921: #include <lemon/dijkstra.h>
alpar@1401: #include <lemon/graph_adaptor.h>
alpar@921: #include <lemon/maps.h>
alpar@899: #include <vector>
alpar@899: 
alpar@921: namespace lemon {
alpar@899: 
alpar@899: /// \addtogroup flowalgs
alpar@899: /// @{
alpar@899: 
alpar@899:   ///\brief Implementation of an algorithm for finding a flow of value \c k 
alpar@899:   ///(for small values of \c k) having minimal total cost between 2 nodes 
alpar@899:   /// 
alpar@899:   ///
jacint@1270:   /// The class \ref lemon::MinCostFlow "MinCostFlow" implements an
jacint@1270:   /// algorithm for finding a flow of value \c k having minimal total
athos@1527:   /// cost from a given source node to a given target node in a
athos@1527:   /// directed graph with a cost function on the edges. To
athos@1527:   /// this end, the edge-capacities and edge-costs have to be
athos@1527:   /// nonnegative.  The edge-capacities should be integers, but the
athos@1527:   /// edge-costs can be integers, reals or of other comparable
athos@1527:   /// numeric type.  This algorithm is intended to be used only for
athos@1527:   /// small values of \c k, since it is only polynomial in k, not in
athos@1527:   /// the length of k (which is log k): in order to find the minimum
athos@1527:   /// cost flow of value \c k it finds the minimum cost flow of value
athos@1527:   /// \c i for every \c i between 0 and \c k.
alpar@899:   ///
alpar@899:   ///\param Graph The directed graph type the algorithm runs on.
alpar@899:   ///\param LengthMap The type of the length map.
alpar@899:   ///\param CapacityMap The capacity map type.
alpar@899:   ///
alpar@899:   ///\author Attila Bernath
alpar@899:   template <typename Graph, typename LengthMap, typename CapacityMap>
alpar@899:   class MinCostFlow {
alpar@899: 
alpar@987:     typedef typename LengthMap::Value Length;
alpar@899: 
alpar@899:     //Warning: this should be integer type
alpar@987:     typedef typename CapacityMap::Value Capacity;
alpar@899:     
alpar@899:     typedef typename Graph::Node Node;
alpar@899:     typedef typename Graph::NodeIt NodeIt;
alpar@899:     typedef typename Graph::Edge Edge;
alpar@899:     typedef typename Graph::OutEdgeIt OutEdgeIt;
alpar@899:     typedef typename Graph::template EdgeMap<int> EdgeIntMap;
alpar@899: 
alpar@1401:     typedef ResGraphAdaptor<const Graph,int,CapacityMap,EdgeIntMap> ResGW;
marci@910:     typedef typename ResGW::Edge ResGraphEdge;
alpar@899: 
marci@941:   protected:
marci@941: 
marci@941:     const Graph& g;
marci@941:     const LengthMap& length;
marci@941:     const CapacityMap& capacity;
marci@941: 
marci@941:     EdgeIntMap flow; 
marci@941:     typedef typename Graph::template NodeMap<Length> PotentialMap;
marci@941:     PotentialMap potential;
marci@941: 
marci@941:     Node s;
marci@941:     Node t;
marci@941:     
marci@941:     Length total_length;
marci@941: 
alpar@899:     class ModLengthMap {   
alpar@899:       typedef typename Graph::template NodeMap<Length> NodeMap;
marci@941:       const ResGW& g;
marci@941:       const LengthMap &length;
alpar@899:       const NodeMap &pot;
alpar@899:     public :
alpar@987:       typedef typename LengthMap::Key Key;
alpar@987:       typedef typename LengthMap::Value Value;
marci@941: 
marci@941:       ModLengthMap(const ResGW& _g, 
marci@941: 		   const LengthMap &_length, const NodeMap &_pot) : 
marci@941: 	g(_g), /*rev(_rev),*/ length(_length), pot(_pot) { }
alpar@899: 	
alpar@987:       Value operator[](typename ResGW::Edge e) const {     
marci@941: 	if (g.forward(e))
alpar@986: 	  return  length[e]-(pot[g.target(e)]-pot[g.source(e)]);   
alpar@899: 	else
alpar@986: 	  return -length[e]-(pot[g.target(e)]-pot[g.source(e)]);   
alpar@899:       }     
alpar@899: 	
marci@941:     }; //ModLengthMap
alpar@899: 
marci@941:     ResGW res_graph;
marci@941:     ModLengthMap mod_length;
marci@941:     Dijkstra<ResGW, ModLengthMap> dijkstra;
alpar@899: 
alpar@899:   public :
alpar@899: 
marci@941:     /*! \brief The constructor of the class.
alpar@899:     
marci@941:     \param _g The directed graph the algorithm runs on. 
athos@1527:     \param _length The length (cost) of the edges. 
marci@941:     \param _cap The capacity of the edges. 
marci@941:     \param _s Source node.
marci@941:     \param _t Target node.
marci@941:     */
marci@941:     MinCostFlow(Graph& _g, LengthMap& _length, CapacityMap& _cap, 
marci@941: 		Node _s, Node _t) : 
marci@941:       g(_g), length(_length), capacity(_cap), flow(_g), potential(_g), 
marci@941:       s(_s), t(_t), 
marci@941:       res_graph(g, capacity, flow), 
marci@941:       mod_length(res_graph, length, potential),
marci@941:       dijkstra(res_graph, mod_length) { 
marci@941:       reset();
marci@941:       }
alpar@899: 
marci@941:     /*! Tries to augment the flow between s and t by 1.
marci@941:       The return value shows if the augmentation is successful.
marci@941:      */
marci@941:     bool augment() {
marci@941:       dijkstra.run(s);
marci@941:       if (!dijkstra.reached(t)) {
alpar@899: 
marci@941: 	//Unsuccessful augmentation.
marci@941: 	return false;
marci@941:       } else {
alpar@899: 
marci@941: 	//We have to change the potential
marci@941: 	for(typename ResGW::NodeIt n(res_graph); n!=INVALID; ++n)
marci@1027: 	  potential.set(n, potential[n]+dijkstra.distMap()[n]);
alpar@899: 	
jacint@1270: 	//Augmenting on the shortest path
alpar@899: 	Node n=t;
alpar@899: 	ResGraphEdge e;
alpar@899: 	while (n!=s){
alpar@899: 	  e = dijkstra.pred(n);
alpar@899: 	  n = dijkstra.predNode(n);
alpar@899: 	  res_graph.augment(e,1);
alpar@899: 	  //Let's update the total length
alpar@899: 	  if (res_graph.forward(e))
alpar@899: 	    total_length += length[e];
alpar@899: 	  else 
alpar@899: 	    total_length -= length[e];	    
alpar@899: 	}
alpar@899: 
marci@941: 	return true;
alpar@899:       }
marci@941:     }
marci@941:     
marci@941:     /*! \brief Runs the algorithm.
marci@941:     
marci@941:     Runs the algorithm.
marci@941:     Returns k if there is a flow of value at least k from s to t.
marci@941:     Otherwise it returns the maximum value of a flow from s to t.
marci@941:     
marci@941:     \param k The value of the flow we are looking for.
marci@941:     
marci@941:     \todo May be it does make sense to be able to start with a nonzero 
marci@941:     feasible primal-dual solution pair as well.
marci@941:     
marci@941:     \todo If the actual flow value is bigger than k, then everything is 
marci@941:     cleared and the algorithm starts from zero flow. Is it a good approach?
marci@941:     */
marci@941:     int run(int k) {
marci@941:       if (flowValue()>k) reset();
marci@941:       while (flowValue()<k && augment()) { }
marci@941:       return flowValue();
marci@941:     }
alpar@899: 
marci@941:     /*! \brief The class is reset to zero flow and potential.
marci@941:       The class is reset to zero flow and potential.
marci@941:      */
marci@941:     void reset() {
marci@941:       total_length=0;
marci@941:       for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
marci@941:       for (typename Graph::NodeIt n(g); n!=INVALID; ++n) potential.set(n, 0);  
marci@941:     }
marci@941: 
marci@941:     /*! Returns the value of the actual flow. 
marci@941:      */
marci@941:     int flowValue() const {
marci@941:       int i=0;
marci@941:       for (typename Graph::OutEdgeIt e(g, s); e!=INVALID; ++e) i+=flow[e];
marci@941:       for (typename Graph::InEdgeIt e(g, s); e!=INVALID; ++e) i-=flow[e];
alpar@899:       return i;
alpar@899:     }
alpar@899: 
athos@1527:     /// Total cost of the found flow.
alpar@899: 
athos@1527:     /// This function gives back the total cost of the found flow.
alpar@899:     Length totalLength(){
alpar@899:       return total_length;
alpar@899:     }
alpar@899: 
alpar@899:     ///Returns a const reference to the EdgeMap \c flow. 
alpar@899: 
alpar@899:     ///Returns a const reference to the EdgeMap \c flow. 
alpar@899:     const EdgeIntMap &getFlow() const { return flow;}
alpar@899: 
marci@941:     /*! \brief Returns a const reference to the NodeMap \c potential (the dual solution).
alpar@899: 
marci@941:     Returns a const reference to the NodeMap \c potential (the dual solution).
marci@941:     */
alpar@899:     const PotentialMap &getPotential() const { return potential;}
alpar@899: 
marci@941:     /*! \brief Checking the complementary slackness optimality criteria.
alpar@899: 
marci@941:     This function checks, whether the given flow and potential 
jacint@1270:     satisfy the complementary slackness conditions (i.e. these are optimal).
marci@941:     This function only checks optimality, doesn't bother with feasibility.
marci@941:     For testing purpose.
marci@941:     */
alpar@899:     bool checkComplementarySlackness(){
alpar@899:       Length mod_pot;
alpar@899:       Length fl_e;
marci@941:         for(typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
alpar@899: 	//C^{\Pi}_{i,j}
alpar@986: 	mod_pot = length[e]-potential[g.target(e)]+potential[g.source(e)];
alpar@899: 	fl_e = flow[e];
alpar@899: 	if (0<fl_e && fl_e<capacity[e]) {
alpar@899: 	  /// \todo better comparison is needed for real types, moreover, 
alpar@899: 	  /// this comparison here is superfluous.
alpar@899: 	  if (mod_pot != 0)
alpar@899: 	    return false;
alpar@899: 	} 
alpar@899: 	else {
alpar@899: 	  if (mod_pot > 0 && fl_e != 0)
alpar@899: 	    return false;
alpar@899: 	  if (mod_pot < 0 && fl_e != capacity[e])
alpar@899: 	    return false;
alpar@899: 	}
alpar@899:       }
alpar@899:       return true;
alpar@899:     }
alpar@899:     
alpar@899:   }; //class MinCostFlow
alpar@899: 
alpar@899:   ///@}
alpar@899: 
alpar@921: } //namespace lemon
alpar@899: 
alpar@921: #endif //LEMON_MIN_COST_FLOW_H