COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/min_cost_flow.h @ 1836:1fee7c6b5129

Last change on this file since 1836:1fee7c6b5129 was 1763:49045f2d28d4, checked in by Balazs Dezso, 18 years ago

pred => predEdge rename

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