COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/athos/mincostflow.h @ 657:531fc5f575ef

Last change on this file since 657:531fc5f575ef was 657:531fc5f575ef, checked in by athos, 20 years ago

Not ready yet.

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