COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/hugo/mincostflows.h @ 860:3577b3db6089

Last change on this file since 860:3577b3db6089 was 860:3577b3db6089, checked in by athos, 16 years ago

Completed documentation for mincostflows and minlengthpaths.

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