src/lemon/min_cost_flow.h
author athos
Wed, 23 Mar 2005 11:51:40 +0000
changeset 1247 60708e1475ae
parent 1027 4ec35d1cd897
child 1270 806451fd084b
permissions -rw-r--r--
Completions.
     1 /* -*- C++ -*-
     2  * src/lemon/min_cost_flow.h - Part of LEMON, a generic C++ optimization library
     3  *
     4  * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     5  * (Egervary Combinatorial Optimization Research Group, EGRES).
     6  *
     7  * Permission to use, modify and distribute this software is granted
     8  * provided that this copyright notice appears in all copies. For
     9  * precise terms see the accompanying LICENSE file.
    10  *
    11  * This software is provided "AS IS" with no warranty of any kind,
    12  * express or implied, and with no claim as to its suitability for any
    13  * purpose.
    14  *
    15  */
    16 
    17 #ifndef LEMON_MIN_COST_FLOW_H
    18 #define LEMON_MIN_COST_FLOW_H
    19 
    20 ///\ingroup flowalgs
    21 ///\file
    22 ///\brief An algorithm for finding a flow of value \c k (for small values of \c k) having minimal total cost 
    23 
    24 
    25 #include <lemon/dijkstra.h>
    26 #include <lemon/graph_wrapper.h>
    27 #include <lemon/maps.h>
    28 #include <vector>
    29 
    30 namespace lemon {
    31 
    32 /// \addtogroup flowalgs
    33 /// @{
    34 
    35   ///\brief Implementation of an algorithm for finding a flow of value \c k 
    36   ///(for small values of \c k) having minimal total cost between 2 nodes 
    37   /// 
    38   ///
    39   /// The class \ref lemon::MinCostFlow "MinCostFlow" implements
    40   /// an algorithm for finding a flow of value \c k 
    41   /// having minimal total cost 
    42   /// from a given source node to a given target node in an
    43   /// edge-weighted directed graph. To this end, 
    44   /// the edge-capacities and edge-weitghs have to be nonnegative. 
    45   /// The edge-capacities should be integers, but the edge-weights can be 
    46   /// integers, reals or of other comparable numeric type.
    47   /// This algorithm is intended to use only for small values of \c k, 
    48   /// since it is only polynomial in k, 
    49   /// not in the length of k (which is log k). 
    50   /// In order to find the minimum cost flow of value \c k it 
    51   /// finds the minimum cost flow of value \c i for every 
    52   /// \c i between 0 and \c k. 
    53   ///
    54   ///\param Graph The directed graph type the algorithm runs on.
    55   ///\param LengthMap The type of the length map.
    56   ///\param CapacityMap The capacity map type.
    57   ///
    58   ///\author Attila Bernath
    59   template <typename Graph, typename LengthMap, typename CapacityMap>
    60   class MinCostFlow {
    61 
    62     typedef typename LengthMap::Value Length;
    63 
    64     //Warning: this should be integer type
    65     typedef typename CapacityMap::Value Capacity;
    66     
    67     typedef typename Graph::Node Node;
    68     typedef typename Graph::NodeIt NodeIt;
    69     typedef typename Graph::Edge Edge;
    70     typedef typename Graph::OutEdgeIt OutEdgeIt;
    71     typedef typename Graph::template EdgeMap<int> EdgeIntMap;
    72 
    73     typedef ResGraphWrapper<const Graph,int,CapacityMap,EdgeIntMap> ResGW;
    74     typedef typename ResGW::Edge ResGraphEdge;
    75 
    76   protected:
    77 
    78     const Graph& g;
    79     const LengthMap& length;
    80     const CapacityMap& capacity;
    81 
    82     EdgeIntMap flow; 
    83     typedef typename Graph::template NodeMap<Length> PotentialMap;
    84     PotentialMap potential;
    85 
    86     Node s;
    87     Node t;
    88     
    89     Length total_length;
    90 
    91     class ModLengthMap {   
    92       typedef typename Graph::template NodeMap<Length> NodeMap;
    93       const ResGW& g;
    94       const LengthMap &length;
    95       const NodeMap &pot;
    96     public :
    97       typedef typename LengthMap::Key Key;
    98       typedef typename LengthMap::Value Value;
    99 
   100       ModLengthMap(const ResGW& _g, 
   101 		   const LengthMap &_length, const NodeMap &_pot) : 
   102 	g(_g), /*rev(_rev),*/ length(_length), pot(_pot) { }
   103 	
   104       Value operator[](typename ResGW::Edge e) const {     
   105 	if (g.forward(e))
   106 	  return  length[e]-(pot[g.target(e)]-pot[g.source(e)]);   
   107 	else
   108 	  return -length[e]-(pot[g.target(e)]-pot[g.source(e)]);   
   109       }     
   110 	
   111     }; //ModLengthMap
   112 
   113     ResGW res_graph;
   114     ModLengthMap mod_length;
   115     Dijkstra<ResGW, ModLengthMap> dijkstra;
   116 
   117   public :
   118 
   119     /*! \brief The constructor of the class.
   120     
   121     \param _g The directed graph the algorithm runs on. 
   122     \param _length The length (weight or cost) of the edges. 
   123     \param _cap The capacity of the edges. 
   124     \param _s Source node.
   125     \param _t Target node.
   126     */
   127     MinCostFlow(Graph& _g, LengthMap& _length, CapacityMap& _cap, 
   128 		Node _s, Node _t) : 
   129       g(_g), length(_length), capacity(_cap), flow(_g), potential(_g), 
   130       s(_s), t(_t), 
   131       res_graph(g, capacity, flow), 
   132       mod_length(res_graph, length, potential),
   133       dijkstra(res_graph, mod_length) { 
   134       reset();
   135       }
   136 
   137     /*! Tries to augment the flow between s and t by 1.
   138       The return value shows if the augmentation is successful.
   139      */
   140     bool augment() {
   141       dijkstra.run(s);
   142       if (!dijkstra.reached(t)) {
   143 
   144 	//Unsuccessful augmentation.
   145 	return false;
   146       } else {
   147 
   148 	//We have to change the potential
   149 	for(typename ResGW::NodeIt n(res_graph); n!=INVALID; ++n)
   150 	  potential.set(n, potential[n]+dijkstra.distMap()[n]);
   151 	
   152 	//Augmenting on the sortest path
   153 	Node n=t;
   154 	ResGraphEdge e;
   155 	while (n!=s){
   156 	  e = dijkstra.pred(n);
   157 	  n = dijkstra.predNode(n);
   158 	  res_graph.augment(e,1);
   159 	  //Let's update the total length
   160 	  if (res_graph.forward(e))
   161 	    total_length += length[e];
   162 	  else 
   163 	    total_length -= length[e];	    
   164 	}
   165 
   166 	return true;
   167       }
   168     }
   169     
   170     /*! \brief Runs the algorithm.
   171     
   172     Runs the algorithm.
   173     Returns k if there is a flow of value at least k from s to t.
   174     Otherwise it returns the maximum value of a flow from s to t.
   175     
   176     \param k The value of the flow we are looking for.
   177     
   178     \todo May be it does make sense to be able to start with a nonzero 
   179     feasible primal-dual solution pair as well.
   180     
   181     \todo If the actual flow value is bigger than k, then everything is 
   182     cleared and the algorithm starts from zero flow. Is it a good approach?
   183     */
   184     int run(int k) {
   185       if (flowValue()>k) reset();
   186       while (flowValue()<k && augment()) { }
   187       return flowValue();
   188     }
   189 
   190     /*! \brief The class is reset to zero flow and potential.
   191       The class is reset to zero flow and potential.
   192      */
   193     void reset() {
   194       total_length=0;
   195       for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
   196       for (typename Graph::NodeIt n(g); n!=INVALID; ++n) potential.set(n, 0);  
   197     }
   198 
   199     /*! Returns the value of the actual flow. 
   200      */
   201     int flowValue() const {
   202       int i=0;
   203       for (typename Graph::OutEdgeIt e(g, s); e!=INVALID; ++e) i+=flow[e];
   204       for (typename Graph::InEdgeIt e(g, s); e!=INVALID; ++e) i-=flow[e];
   205       return i;
   206     }
   207 
   208     /// Total weight of the found flow.
   209 
   210     /// This function gives back the total weight of the found flow.
   211     Length totalLength(){
   212       return total_length;
   213     }
   214 
   215     ///Returns a const reference to the EdgeMap \c flow. 
   216 
   217     ///Returns a const reference to the EdgeMap \c flow. 
   218     const EdgeIntMap &getFlow() const { return flow;}
   219 
   220     /*! \brief Returns a const reference to the NodeMap \c potential (the dual solution).
   221 
   222     Returns a const reference to the NodeMap \c potential (the dual solution).
   223     */
   224     const PotentialMap &getPotential() const { return potential;}
   225 
   226     /*! \brief Checking the complementary slackness optimality criteria.
   227 
   228     This function checks, whether the given flow and potential 
   229     satisfiy the complementary slackness cnditions (i.e. these are optimal).
   230     This function only checks optimality, doesn't bother with feasibility.
   231     For testing purpose.
   232     */
   233     bool checkComplementarySlackness(){
   234       Length mod_pot;
   235       Length fl_e;
   236         for(typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
   237 	//C^{\Pi}_{i,j}
   238 	mod_pot = length[e]-potential[g.target(e)]+potential[g.source(e)];
   239 	fl_e = flow[e];
   240 	if (0<fl_e && fl_e<capacity[e]) {
   241 	  /// \todo better comparison is needed for real types, moreover, 
   242 	  /// this comparison here is superfluous.
   243 	  if (mod_pot != 0)
   244 	    return false;
   245 	} 
   246 	else {
   247 	  if (mod_pot > 0 && fl_e != 0)
   248 	    return false;
   249 	  if (mod_pot < 0 && fl_e != capacity[e])
   250 	    return false;
   251 	}
   252       }
   253       return true;
   254     }
   255     
   256   }; //class MinCostFlow
   257 
   258   ///@}
   259 
   260 } //namespace lemon
   261 
   262 #endif //LEMON_MIN_COST_FLOW_H