COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/lemon/min_cost_flow.h @ 941:186aa53d2802

Last change on this file since 941:186aa53d2802 was 941:186aa53d2802, checked in by marci, 19 years ago

Suurballe and MinCostFlow? classes are now able to increase the flow 1 by 1 with
this->augment()

File size: 7.8 KB
Line 
1/* -*- C++ -*-
2 * src/lemon/min_cost_flow.h - Part of LEMON, a generic C++ optimization library
3 *
4 * Copyright (C) 2004 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
30namespace 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::ValueType Length;
63
64    //Warning: this should be integer type
65    typedef typename CapacityMap::ValueType 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::KeyType KeyType;
98      typedef typename LengthMap::ValueType ValueType;
99
100      ModLengthMap(const ResGW& _g,
101                   const LengthMap &_length, const NodeMap &_pot) :
102        g(_g), /*rev(_rev),*/ length(_length), pot(_pot) { }
103       
104      ValueType operator[](typename ResGW::Edge e) const {     
105        if (g.forward(e))
106          return  length[e]-(pot[g.head(e)]-pot[g.tail(e)]);   
107        else
108          return -length[e]-(pot[g.head(e)]-pot[g.tail(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[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.head(e)]+potential[g.tail(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
Note: See TracBrowser for help on using the repository browser.