COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/athos/mincostflow.h @ 645:d93d8b9906d1

Last change on this file since 645:d93d8b9906d1 was 645:d93d8b9906d1, checked in by athos, 17 years ago

I don't really feel like working on this at the moment.

File size: 8.2 KB
Line 
1// -*- c++ -*-
2#ifndef HUGO_MINCOSTFLOW_H
3#define HUGO_MINCOSTFLOW_H
4
5///\ingroup galgs
6///\file
7///\brief An algorithm for finding the minimum cost flow of given value in an uncapacitated network
8
9#include <hugo/dijkstra.h>
10#include <hugo/graph_wrapper.h>
11#include <hugo/maps.h>
12#include <vector>
13#include <for_each_macros.h>
14
15namespace hugo {
16
17/// \addtogroup galgs
18/// @{
19
20  ///\brief Implementation of an algorithm for finding the minimum cost flow
21  /// of given value in an uncapacitated network
22  ///
23  ///
24  /// The class \ref hugo::MinCostFlow "MinCostFlow" implements
25  /// an algorithm for solving the following general minimum cost flow problem>
26  ///
27  ///
28  ///
29  /// \warning It is assumed here that the problem has a feasible solution
30  ///
31  /// The range of the length (weight) function is nonnegative reals but
32  /// the range of capacity function is the set of nonnegative integers.
33  /// It is not a polinomial time algorithm for counting the minimum cost
34  /// maximal flow, since it counts the minimum cost flow for every value 0..M
35  /// where \c M is the value of the maximal flow.
36  ///
37  ///\author Attila Bernath
38  template <typename Graph, typename LengthMap, typename SupplyDemandMap>
39  class MinCostFlow {
40
41    typedef typename LengthMap::ValueType Length;
42
43
44    typedef typename SupplyDemandMap::ValueType SupplyDemand;
45   
46    typedef typename Graph::Node Node;
47    typedef typename Graph::NodeIt NodeIt;
48    typedef typename Graph::Edge Edge;
49    typedef typename Graph::OutEdgeIt OutEdgeIt;
50    typedef typename Graph::template EdgeMap<int> EdgeIntMap;
51
52    //    typedef ConstMap<Edge,int> ConstMap;
53
54    typedef ResGraphWrapper<const Graph,int,CapacityMap,EdgeIntMap> ResGraphType;
55    typedef typename ResGraphType::Edge ResGraphEdge;
56
57    class ModLengthMap {   
58      //typedef typename ResGraphType::template NodeMap<Length> NodeMap;
59      typedef typename Graph::template NodeMap<Length> NodeMap;
60      const ResGraphType& G;
61      //      const EdgeIntMap& rev;
62      const LengthMap &ol;
63      const NodeMap &pot;
64    public :
65      typedef typename LengthMap::KeyType KeyType;
66      typedef typename LengthMap::ValueType ValueType;
67       
68      ValueType operator[](typename ResGraphType::Edge e) const {     
69        if (G.forward(e))
70          return  ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);   
71        else
72          return -ol[e]-(pot[G.head(e)]-pot[G.tail(e)]);   
73      }     
74       
75      ModLengthMap(const ResGraphType& _G,
76                   const LengthMap &o,  const NodeMap &p) :
77        G(_G), /*rev(_rev),*/ ol(o), pot(p){};
78    };//ModLengthMap
79
80
81  protected:
82   
83    //Input
84    const Graph& G;
85    const LengthMap& length;
86    const SupplyDemandMap& supply_demand;//supply or demand of nodes
87
88
89    //auxiliary variables
90
91    //To store the flow
92    EdgeIntMap flow;
93    //To store the potentila (dual variables)
94    typename Graph::template NodeMap<Length> potential;
95    //To store excess-deficit values
96    SupplyDemandMap excess_deficit;
97   
98
99    Length total_length;
100
101
102  public :
103
104
105    MinCostFlow(Graph& _G, LengthMap& _length, SupplyDemandMap& _supply_demand) : G(_G),
106      length(_length), supply_demand(_supply_demand), flow(_G), potential(_G){ }
107
108   
109    ///Runs the algorithm.
110
111    ///Runs the algorithm.
112
113    ///\todo May be it does make sense to be able to start with a nonzero
114    /// feasible primal-dual solution pair as well.
115    int run() {
116
117      //Resetting variables from previous runs
118      //total_length = 0;
119
120      typedef typename Graph::template NodeMap<int> HeapMap;
121      typedef Heap<Node, SupplyDemand, typename Graph::template NodeMap<int>,
122        std::greater<SupplyDemand> >    HeapType;
123
124      //A heap for the excess nodes
125      HeapMap excess_nodes_map(G,-1);
126      HeapType excess_nodes(excess_nodes_map);
127
128      //A heap for the deficit nodes
129      HeapMap deficit_nodes_map(G,-1);
130      HeapType deficit_nodes(deficit_nodes_map);
131
132     
133      FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
134        flow.set(e,0);
135      }
136
137      //Initial value for delta
138      SupplyDemand delta = 0;
139
140      FOR_EACH_LOC(typename Graph::NodeIt, n, G){
141        excess_deficit.set(n,supply_demand[n]);
142        //A supply node
143        if (excess_deficit[n] > 0){
144          excess_nodes.push(n,excess_deficit[n]);
145        }
146        //A demand node
147        if (excess_deficit[n] < 0){
148          deficit_nodes.push(n, - excess_deficit[n]);
149        }
150        //Finding out starting value of delta
151        if (delta < abs(excess_deficit[n])){
152          delta = abs(excess_deficit[n]);
153        }
154        //Initialize the copy of the Dijkstra potential to zero
155        potential.set(n,0);
156      }
157
158      //It'll be allright as an initial value, though this value
159      //can be the maximum deficit here
160      SupplyDemand max_excess = delta;
161     
162      //We need a residual graph which is uncapacitated
163      ResGraphType res_graph(G, flow);
164
165
166     
167      ModLengthMap mod_length(res_graph, length, potential);
168
169      Dijkstra<ResGraphType, ModLengthMap> dijkstra(res_graph, mod_length);
170
171
172      while (max_excess > 0){
173
174        /*
175         * Beginning of the delta scaling phase
176        */
177       
178        //Merge and stuff
179
180        Node s = excess_nodes.top();
181        SupplyDemand max_excess = excess_nodes[s];
182        Node t = deficit_nodes.top();
183        if (max_excess < dificit_nodes[t]){
184          max_excess = dificit_nodes[t];
185        }
186
187
188        while(max_excess > 0){
189
190         
191          //s es t valasztasa
192
193          //Dijkstra part       
194          dijkstra.run(s);
195
196          /*We know from theory that t can be reached
197          if (!dijkstra.reached(t)){
198            //There are no k paths from s to t
199            break;
200          };
201          */
202         
203          //We have to change the potential
204          FOR_EACH_LOC(typename ResGraphType::NodeIt, n, res_graph){
205            potential[n] += dijkstra.distMap()[n];
206          }
207
208
209          //Augmenting on the sortest path
210          Node n=t;
211          ResGraphEdge e;
212          while (n!=s){
213            e = dijkstra.pred(n);
214            n = dijkstra.predNode(n);
215            res_graph.augment(e,delta);
216            /*
217            //Let's update the total length
218            if (res_graph.forward(e))
219              total_length += length[e];
220            else
221              total_length -= length[e];           
222            */
223          }
224
225          //Update the excess_nodes heap
226          if (delta >= excess_nodes[s]){
227            if (delta > excess_nodes[s])
228              deficit_nodes.push(s,delta - excess_nodes[s]);
229            excess_nodes.pop();
230           
231          }
232          else{
233            excess_nodes[s] -= delta;
234          }
235          //Update the deficit_nodes heap
236          if (delta >= deficit_nodes[t]){
237            if (delta > deficit_nodes[t])
238              excess_nodes.push(t,delta - deficit_nodes[t]);
239            deficit_nodes.pop();
240           
241          }
242          else{
243            deficit_nodes[t] -= delta;
244          }
245          //Dijkstra part ends here
246        }
247
248        /*
249         * End of the delta scaling phase
250        */
251
252        //Whatever this means
253        delta = delta / 2;
254
255        /*This is not necessary here
256        //Update the max_excess
257        max_excess = 0;
258        FOR_EACH_LOC(typename Graph::NodeIt, n, G){
259          if (max_excess < excess_deficit[n]){
260            max_excess = excess_deficit[n];
261          }
262        }
263        */
264        //Reset delta if still too big
265        if (8*number_of_nodes*max_excess <= delta){
266          delta = max_excess;
267         
268        }
269         
270      }//while(max_excess > 0)
271     
272
273      return i;
274    }
275
276
277
278
279    ///This function gives back the total length of the found paths.
280    ///Assumes that \c run() has been run and nothing changed since then.
281    Length totalLength(){
282      return total_length;
283    }
284
285    ///Returns a const reference to the EdgeMap \c flow. \pre \ref run() must
286    ///be called before using this function.
287    const EdgeIntMap &getFlow() const { return flow;}
288
289  ///Returns a const reference to the NodeMap \c potential (the dual solution).
290    /// \pre \ref run() must be called before using this function.
291    const EdgeIntMap &getPotential() const { return potential;}
292
293    ///This function checks, whether the given solution is optimal
294    ///Running after a \c run() should return with true
295    ///In this "state of the art" this only check optimality, doesn't bother with feasibility
296    ///
297    ///\todo Is this OK here?
298    bool checkComplementarySlackness(){
299      Length mod_pot;
300      Length fl_e;
301      FOR_EACH_LOC(typename Graph::EdgeIt, e, G){
302        //C^{\Pi}_{i,j}
303        mod_pot = length[e]-potential[G.head(e)]+potential[G.tail(e)];
304        fl_e = flow[e];
305        //      std::cout << fl_e << std::endl;
306        if (0<fl_e && fl_e<capacity[e]){
307          if (mod_pot != 0)
308            return false;
309        }
310        else{
311          if (mod_pot > 0 && fl_e != 0)
312            return false;
313          if (mod_pot < 0 && fl_e != capacity[e])
314            return false;
315        }
316      }
317      return true;
318    }
319   
320
321  }; //class MinCostFlow
322
323  ///@}
324
325} //namespace hugo
326
327#endif //HUGO_MINCOSTFLOW_H
Note: See TracBrowser for help on using the repository browser.