COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/min_cost_flow.h @ 2051:08652c1763f6

Last change on this file since 2051:08652c1763f6 was 1956:a055123339d5, checked in by Alpar Juttner, 18 years ago

Unified copyright notices

File size: 7.8 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2006
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
12 *
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_MIN_COST_FLOW_H
20#define LEMON_MIN_COST_FLOW_H
21
22///\ingroup flowalgs
23///\file
24///\brief An algorithm for finding a flow of value \c k (for small values of \c k) having minimal total cost
25
26
27#include <lemon/dijkstra.h>
28#include <lemon/graph_adaptor.h>
29#include <lemon/maps.h>
30#include <vector>
31
32namespace lemon {
33
34/// \addtogroup flowalgs
35/// @{
36
37  ///\brief Implementation of an algorithm for finding a flow of value \c k
38  ///(for small values of \c k) having minimal total cost between 2 nodes
39  ///
40  ///
41  /// The class \ref lemon::MinCostFlow "MinCostFlow" implements an
42  /// algorithm for finding a flow of value \c k having minimal total
43  /// cost from a given source node to a given target node in a
44  /// directed graph with a cost function on the edges. To
45  /// this end, the edge-capacities and edge-costs have to be
46  /// nonnegative.  The edge-capacities should be integers, but the
47  /// edge-costs can be integers, reals or of other comparable
48  /// numeric type.  This algorithm is intended to be used only for
49  /// small values of \c k, since it is only polynomial in k, not in
50  /// the length of k (which is log k): in order to find the minimum
51  /// cost flow of value \c k it finds the minimum cost flow of value
52  /// \c i for every \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 ResGraphAdaptor<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 (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 shortest path
153        Node n=t;
154        ResGraphEdge e;
155        while (n!=s){
156          e = dijkstra.predEdge(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 cost of the found flow.
209
210    /// This function gives back the total cost 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    satisfy the complementary slackness conditions (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
Note: See TracBrowser for help on using the repository browser.