COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/hugo/min_cost_flow.h @ 901:69a8e672acb1

Last change on this file since 901:69a8e672acb1 was 901:69a8e672acb1, checked in by marci, 20 years ago

correction of HUGO_... preproc defines.

File size: 6.9 KB
Line 
1// -*- c++ -*-
2#ifndef HUGO_MIN_COST_FLOW_H
3#define HUGO_MIN_COST_FLOW_H
4
5///\ingroup flowalgs
6///\file
7///\brief An algorithm for finding a flow of value \c k (for small values of \c k) having minimal total cost
8
9
10#include <hugo/dijkstra.h>
11#include <hugo/graph_wrapper.h>
12#include <hugo/maps.h>
13#include <vector>
14
15namespace hugo {
16
17/// \addtogroup flowalgs
18/// @{
19
20  ///\brief Implementation of an algorithm for finding a flow of value \c k
21  ///(for small values of \c k) having minimal total cost between 2 nodes
22  ///
23  ///
24  /// The class \ref hugo::MinCostFlow "MinCostFlow" implements
25  /// an algorithm for finding a flow of value \c k
26  /// having minimal total cost
27  /// from a given source node to a given target node in an
28  /// edge-weighted directed graph. To this end,
29  /// the edge-capacities and edge-weitghs have to be nonnegative.
30  /// The edge-capacities should be integers, but the edge-weights can be
31  /// integers, reals or of other comparable numeric type.
32  /// This algorithm is intended to use only for small values of \c k,
33  /// since it is only polynomial in k,
34  /// not in the length of k (which is log k).
35  /// In order to find the minimum cost flow of value \c k it
36  /// finds the minimum cost flow of value \c i for every
37  /// \c i between 0 and \c k.
38  ///
39  ///\param Graph The directed graph type the algorithm runs on.
40  ///\param LengthMap The type of the length map.
41  ///\param CapacityMap The capacity map type.
42  ///
43  ///\author Attila Bernath
44  template <typename Graph, typename LengthMap, typename CapacityMap>
45  class MinCostFlow {
46
47    typedef typename LengthMap::ValueType Length;
48
49    //Warning: this should be integer type
50    typedef typename CapacityMap::ValueType Capacity;
51   
52    typedef typename Graph::Node Node;
53    typedef typename Graph::NodeIt NodeIt;
54    typedef typename Graph::Edge Edge;
55    typedef typename Graph::OutEdgeIt OutEdgeIt;
56    typedef typename Graph::template EdgeMap<int> EdgeIntMap;
57
58
59    typedef ResGraphWrapper<const Graph,int,CapacityMap,EdgeIntMap> ResGraphType;
60    typedef typename ResGraphType::Edge ResGraphEdge;
61
62    class ModLengthMap {   
63      typedef typename Graph::template NodeMap<Length> NodeMap;
64      const ResGraphType& G;
65      const LengthMap &ol;
66      const NodeMap &pot;
67    public :
68      typedef typename LengthMap::KeyType KeyType;
69      typedef typename LengthMap::ValueType ValueType;
70       
71      ValueType operator[](typename ResGraphType::Edge e) const {     
72        if (G.forward(e))
73          return  ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);   
74        else
75          return -ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);   
76      }     
77       
78      ModLengthMap(const ResGraphType& _G,
79                   const LengthMap &o,  const NodeMap &p) :
80        G(_G), /*rev(_rev),*/ ol(o), pot(p){};
81    };//ModLengthMap
82
83
84  protected:
85   
86    //Input
87    const Graph& G;
88    const LengthMap& length;
89    const CapacityMap& capacity;
90
91
92    //auxiliary variables
93
94    //To store the flow
95    EdgeIntMap flow;
96    //To store the potential (dual variables)
97    typedef typename Graph::template NodeMap<Length> PotentialMap;
98    PotentialMap potential;
99   
100
101    Length total_length;
102
103
104  public :
105
106    /// The constructor of the class.
107   
108    ///\param _G The directed graph the algorithm runs on.
109    ///\param _length The length (weight or cost) of the edges.
110    ///\param _cap The capacity of the edges.
111    MinCostFlow(Graph& _G, LengthMap& _length, CapacityMap& _cap) : G(_G),
112      length(_length), capacity(_cap), flow(_G), potential(_G){ }
113
114   
115    ///Runs the algorithm.
116   
117    ///Runs the algorithm.
118    ///Returns k if there is a flow of value at least k edge-disjoint
119    ///from s to t.
120    ///Otherwise it returns the maximum value of a flow from s to t.
121    ///
122    ///\param s The source node.
123    ///\param t The target node.
124    ///\param k The value of the flow we are looking for.
125    ///
126    ///\todo May be it does make sense to be able to start with a nonzero
127    /// feasible primal-dual solution pair as well.
128    int run(Node s, Node t, int k) {
129
130      //Resetting variables from previous runs
131      total_length = 0;
132     
133      for (typename Graph::EdgeIt e(G); e!=INVALID; ++e) flow.set(e, 0);
134
135      //Initialize the potential to zero
136      for (typename Graph::NodeIt n(G); n!=INVALID; ++n) potential.set(n, 0);
137     
138     
139      //We need a residual graph
140      ResGraphType res_graph(G, capacity, flow);
141
142
143      ModLengthMap mod_length(res_graph, length, potential);
144
145      Dijkstra<ResGraphType, ModLengthMap> dijkstra(res_graph, mod_length);
146
147      int i;
148      for (i=0; i<k; ++i){
149        dijkstra.run(s);
150        if (!dijkstra.reached(t)){
151          //There are no flow of value k from s to t
152          break;
153        };
154       
155        //We have to change the potential
156        for(typename ResGraphType::NodeIt n(res_graph); n!=INVALID; ++n)
157          potential[n] += dijkstra.distMap()[n];
158
159
160        //Augmenting on the sortest path
161        Node n=t;
162        ResGraphEdge e;
163        while (n!=s){
164          e = dijkstra.pred(n);
165          n = dijkstra.predNode(n);
166          res_graph.augment(e,1);
167          //Let's update the total length
168          if (res_graph.forward(e))
169            total_length += length[e];
170          else
171            total_length -= length[e];     
172        }
173
174         
175      }
176     
177
178      return i;
179    }
180
181
182
183    /// Gives back the total weight of the found flow.
184
185    ///This function gives back the total weight of the found flow.
186    ///Assumes that \c run() has been run and nothing changed since then.
187    Length totalLength(){
188      return total_length;
189    }
190
191    ///Returns a const reference to the EdgeMap \c flow.
192
193    ///Returns a const reference to the EdgeMap \c flow.
194    ///\pre \ref run() must
195    ///be called before using this function.
196    const EdgeIntMap &getFlow() const { return flow;}
197
198    ///Returns a const reference to the NodeMap \c potential (the dual solution).
199
200    ///Returns a const reference to the NodeMap \c potential (the dual solution).
201    /// \pre \ref run() must be called before using this function.
202    const PotentialMap &getPotential() const { return potential;}
203
204    /// Checking the complementary slackness optimality criteria
205
206    ///This function checks, whether the given solution is optimal
207    ///If executed after the call of \c run() then it should return with true.
208    ///This function only checks optimality, doesn't bother with feasibility.
209    ///It is meant for testing purposes.
210    ///
211    bool checkComplementarySlackness(){
212      Length mod_pot;
213      Length fl_e;
214        for(typename Graph::EdgeIt e(G); e!=INVALID; ++e) {
215        //C^{\Pi}_{i,j}
216        mod_pot = length[e]-potential[G.head(e)]+potential[G.tail(e)];
217        fl_e = flow[e];
218        if (0<fl_e && fl_e<capacity[e]) {
219          /// \todo better comparison is needed for real types, moreover,
220          /// this comparison here is superfluous.
221          if (mod_pot != 0)
222            return false;
223        }
224        else {
225          if (mod_pot > 0 && fl_e != 0)
226            return false;
227          if (mod_pot < 0 && fl_e != capacity[e])
228            return false;
229        }
230      }
231      return true;
232    }
233   
234
235  }; //class MinCostFlow
236
237  ///@}
238
239} //namespace hugo
240
241#endif //HUGO_MIN_COST_FLOW_H
Note: See TracBrowser for help on using the repository browser.