COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/athos/mincostflow.h @ 661:d306e777117e

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

Corrected some obvious errors.

File size: 12.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 <list>
14#include <hugo/for_each_macros.h>
15#include <hugo/unionfind.h>
16
17namespace hugo {
18
19/// \addtogroup galgs
20/// @{
21
22  ///\brief Implementation of an algorithm for solving the minimum cost general
23  /// flow problem 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 cost (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 CostMap, typename SupplyDemandMap>
41  class MinCostFlow {
42
43    typedef typename CostMap::ValueType Cost;
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<SupplyDemand> FlowMap;
53    typedef ConstMap<Edge,SupplyDemand> ConstEdgeMap;
54
55    //    typedef ConstMap<Edge,int> ConstMap;
56
57    typedef ResGraphWrapper<const Graph,int,ConstEdgeMap,FlowMap> ResGraph;
58    typedef typename ResGraph::Edge ResGraphEdge;
59
60    class ModCostMap {   
61      //typedef typename ResGraph::template NodeMap<Cost> NodeMap;
62      typedef typename Graph::template NodeMap<Cost> NodeMap;
63      const ResGraph& res_graph;
64      //      const EdgeIntMap& rev;
65      const CostMap &ol;
66      const NodeMap &pot;
67    public :
68      typedef typename CostMap::KeyType KeyType;
69      typedef typename CostMap::ValueType ValueType;
70       
71      ValueType operator[](typename ResGraph::Edge e) const {     
72        if (res_graph.forward(e))
73          return  ol[e]-(pot[res_graph.head(e)]-pot[res_graph.tail(e)]);   
74        else
75          return -ol[e]-(pot[res_graph.head(e)]-pot[res_graph.tail(e)]);   
76      }     
77       
78      ModCostMap(const ResGraph& _res_graph,
79                   const CostMap &o,  const NodeMap &p) :
80        res_graph(_res_graph), /*rev(_rev),*/ ol(o), pot(p){};
81    };//ModCostMap
82
83
84  protected:
85   
86    //Input
87    const Graph& graph;
88    const CostMap& cost;
89    const SupplyDemandMap& supply_demand;//supply or demand of nodes
90
91
92    //auxiliary variables
93
94    //To store the flow
95    FlowMap flow;
96    //To store the potentila (dual variables)
97    typename Graph::template NodeMap<Cost> potential;
98    //To store excess-deficit values
99    SupplyDemandMap excess_deficit;
100   
101
102    Cost total_cost;
103
104
105  public :
106
107
108    MinCostFlow(Graph& _graph, CostMap& _cost, SupplyDemandMap& _supply_demand) : graph(_graph),
109      cost(_cost), supply_demand(_supply_demand), flow(_graph), potential(_graph){ }
110
111   
112    ///Runs the algorithm.
113
114    ///Runs the algorithm.
115
116    ///\todo May be it does make sense to be able to start with a nonzero
117    /// feasible primal-dual solution pair as well.
118    void run() {
119
120      //Resetting variables from previous runs
121      //total_cost = 0;
122
123      typedef typename Graph::template NodeMap<int> HeapMap;
124      typedef Heap< Node, SupplyDemand, typename Graph::template NodeMap<int>,
125        std::greater<SupplyDemand> >    HeapType;
126
127      //A heap for the excess nodes
128      HeapMap excess_nodes_map(graph,-1);
129      HeapType excess_nodes(excess_nodes_map);
130
131      //A heap for the deficit nodes
132      HeapMap deficit_nodes_map(graph,-1);
133      HeapType deficit_nodes(deficit_nodes_map);
134
135      //A container to store nonabundant arcs
136      list<Edge> nonabundant_arcs;
137
138       
139      FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
140        flow.set(e,0);
141        nonabundant_arcs.push_back(e);
142      }
143
144      //Initial value for delta
145      SupplyDemand delta = 0;
146
147      typedef UnionFindEnum<Node, Graph::template NodeMap> UFE;
148
149      //A union-find structure to store the abundant components
150      UFE::MapType abund_comp_map(graph);
151      UFE abundant_components(abund_comp_map);
152
153
154
155      FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
156        excess_deficit.set(n,supply_demand[n]);
157        //A supply node
158        if (excess_deficit[n] > 0){
159          excess_nodes.push(n,excess_deficit[n]);
160        }
161        //A demand node
162        if (excess_deficit[n] < 0){
163          deficit_nodes.push(n, - excess_deficit[n]);
164        }
165        //Finding out starting value of delta
166        if (delta < abs(excess_deficit[n])){
167          delta = abs(excess_deficit[n]);
168        }
169        //Initialize the copy of the Dijkstra potential to zero
170        potential.set(n,0);
171        //Every single point is an abundant component initially
172        abundant_components.insert(n);
173      }
174
175      //It'll be allright as an initial value, though this value
176      //can be the maximum deficit here
177      SupplyDemand max_excess = delta;
178     
179      ///\bug This is a serious cheat here, before we have an uncapacitated ResGraph
180      ConstEdgeMap const_inf_map(MAX_INT);
181     
182      //We need a residual graph which is uncapacitated
183      ResGraph res_graph(graph, const_inf_map, flow);
184     
185      //An EdgeMap to tell which arcs are abundant
186      template typename Graph::EdgeMap<bool> abundant_arcs(graph);
187
188      //Let's construct the sugraph consisting only of the abundant edges
189      typedef ConstMap< typename Graph::Node, bool > ConstNodeMap;
190      ConstNodeMap const_true_map(true);
191      typedef SubGraphWrapper< Graph, ConstNodeMap,
192         template typename Graph::EdgeMap<bool> >
193        AbundantGraph;
194      AbundantGraph abundant_graph(graph, const_true_map, abundant_arcs );
195     
196      //Let's construct the residual graph for the abundant graph
197      typedef ResGraphWrapper<const AbundantGraph,int,CapacityMap,EdgeIntMap>
198        ResAbGraph;
199      //Again uncapacitated
200      ResAbGraph res_ab_graph(abundant_graph, const_inf_map, flow);
201     
202      //We need things for the bfs
203      typename ResAbGraph::NodeMap<bool> bfs_reached(res_ab_graph);
204      typename ResAbGraph::NodeMap<typename ResAbGraph::Edge>
205        bfs_pred(res_ab_graph);
206      NullMap<typename ResAbGraph::Node, int> bfs_dist_dummy(res_ab_graph);
207      //We want to run bfs-es (more) on this graph 'res_ab_graph'
208      Bfs < ResAbGraph ,
209        typename ResAbGraph::NodeMap<bool>,
210        typename ResAbGraph::NodeMap<typename ResAbGraph::Edge>,
211        NullMap<typename ResAbGraph::Node, int> >
212        bfs(res_ab_graph, bfs_reached, bfs_pred, bfs_dist_dummy);
213     
214      ModCostMap mod_cost(res_graph, cost, potential);
215
216      Dijkstra<ResGraph, ModCostMap> dijkstra(res_graph, mod_cost);
217
218
219      while (max_excess > 0){
220
221        //Reset delta if still too big
222        if (8*number_of_nodes*max_excess <= delta){
223          delta = max_excess;
224         
225        }
226
227        /*
228         * Beginning of the delta scaling phase
229        */
230        //Merge and stuff
231        {
232          SupplyDemand buf=8*number_of_nodes*delta;
233          list<Edge>::iterator i = nonabundant_arcs.begin();
234          while ( i != nonabundant_arcs.end() ){
235            if (flow[i]>=buf){
236              Node a = abundant_components.find(res_graph.head(i));
237              Node b = abundant_components.find(res_graph.tail(i));
238              //Merge
239              if (a != b){
240                abundant_components.join(a,b);
241                //We want to push the smaller
242                //Which has greater absolut value excess/deficit
243                Node root=(abs(excess_deficit[a])>abs(excess_deficit[b]))?a:b;
244                //Which is the other
245                Node non_root = ( a == root ) ? b : a ;
246                abundant_components.makeRep(root);
247                SupplyDemand qty_to_augment = abs(excess_deficit[non_root]);
248                //Push the positive value
249                if (excess_deficit[non_root] < 0)
250                  swap(root, non_root);
251                //If the non_root node has excess/deficit at all
252                if (qty_to_augment>0){
253                  //Find path and augment
254                  bfs.run(non_root);
255                  //root should be reached
256                 
257                  //Augmenting on the found path
258                  Node n=root;
259                  ResGraphEdge e;
260                  while (n!=non_root){
261                    e = bfs_pred(n);
262                    n = res_graph.tail(e);
263                    res_graph.augment(e,qty_to_augment);
264                  }
265         
266                  //We know that non_root had positive excess
267                  excess_nodes[non_root] -= qty_to_augment;
268                  //But what about root node
269                  //It might have been positive and so became larger
270                  if (excess_deficit[root]>0){
271                    excess_nodes[root] += qty_to_augment;
272                  }
273                  else{
274                    //Or negative but not turned into positive
275                    deficit_nodes[root] -= qty_to_augment;
276                  }
277
278                  //Update the excess_deficit map
279                  excess_deficit[non_root] -= qty_to_augment;
280                  excess_deficit[root] += qty_to_augment;
281
282                 
283                }
284              }
285              //What happens to i?
286              //Marci and Zsolt says I shouldn't do such things
287              nonabundant_arcs.erase(i++);
288              abundant_arcs[i] = true;
289            }
290            else
291              ++i;
292          }
293        }
294
295
296        Node s = excess_nodes.top();
297        SupplyDemand max_excess = excess_nodes[s];
298        Node t = deficit_nodes.top();
299        if (max_excess < deficit_nodes[t]){
300          max_excess = deficit_nodes[t];
301        }
302
303
304        while(max_excess > (n-1)*delta/n){
305         
306         
307          //s es t valasztasa
308         
309          //Dijkstra part       
310          dijkstra.run(s);
311         
312          /*We know from theory that t can be reached
313          if (!dijkstra.reached(t)){
314            //There are no k paths from s to t
315            break;
316          };
317          */
318         
319          //We have to change the potential
320          FOR_EACH_LOC(typename ResGraph::NodeIt, n, res_graph){
321            potential[n] += dijkstra.distMap()[n];
322          }
323
324
325          //Augmenting on the sortest path
326          Node n=t;
327          ResGraphEdge e;
328          while (n!=s){
329            e = dijkstra.pred(n);
330            n = dijkstra.predNode(n);
331            res_graph.augment(e,delta);
332            /*
333            //Let's update the total cost
334            if (res_graph.forward(e))
335              total_cost += cost[e];
336            else
337              total_cost -= cost[e];       
338            */
339          }
340         
341          //Update the excess_deficit map
342          excess_deficit[s] -= delta;
343          excess_deficit[t] += delta;
344         
345
346          //Update the excess_nodes heap
347          if (delta >= excess_nodes[s]){
348            if (delta > excess_nodes[s])
349              deficit_nodes.push(s,delta - excess_nodes[s]);
350            excess_nodes.pop();
351           
352          }
353          else{
354            excess_nodes[s] -= delta;
355          }
356          //Update the deficit_nodes heap
357          if (delta >= deficit_nodes[t]){
358            if (delta > deficit_nodes[t])
359              excess_nodes.push(t,delta - deficit_nodes[t]);
360            deficit_nodes.pop();
361           
362          }
363          else{
364            deficit_nodes[t] -= delta;
365          }
366          //Dijkstra part ends here
367         
368          //Choose s and t again
369          s = excess_nodes.top();
370          max_excess = excess_nodes[s];
371          t = deficit_nodes.top();
372          if (max_excess < deficit_nodes[t]){
373            max_excess = deficit_nodes[t];
374          }
375
376        }
377
378        /*
379         * End of the delta scaling phase
380        */
381
382        //Whatever this means
383        delta = delta / 2;
384
385        /*This is not necessary here
386        //Update the max_excess
387        max_excess = 0;
388        FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
389          if (max_excess < excess_deficit[n]){
390            max_excess = excess_deficit[n];
391          }
392        }
393        */
394
395         
396      }//while(max_excess > 0)
397     
398
399      return i;
400    }
401
402
403
404
405    ///This function gives back the total cost of the found paths.
406    ///Assumes that \c run() has been run and nothing changed since then.
407    Cost totalCost(){
408      return total_cost;
409    }
410
411    ///Returns a const reference to the EdgeMap \c flow. \pre \ref run() must
412    ///be called before using this function.
413    const EdgeIntMap &getFlow() const { return flow;}
414
415  ///Returns a const reference to the NodeMap \c potential (the dual solution).
416    /// \pre \ref run() must be called before using this function.
417    const EdgeIntMap &getPotential() const { return potential;}
418
419    ///This function checks, whether the given solution is optimal
420    ///Running after a \c run() should return with true
421    ///In this "state of the art" this only check optimality, doesn't bother with feasibility
422    ///
423    ///\todo Is this OK here?
424    bool checkComplementarySlackness(){
425      Cost mod_pot;
426      Cost fl_e;
427      FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
428        //C^{\Pi}_{i,j}
429        mod_pot = cost[e]-potential[graph.head(e)]+potential[graph.tail(e)];
430        fl_e = flow[e];
431        //      std::cout << fl_e << std::endl;
432        if (0<fl_e && fl_e<capacity[e]){
433          if (mod_pot != 0)
434            return false;
435        }
436        else{
437          if (mod_pot > 0 && fl_e != 0)
438            return false;
439          if (mod_pot < 0 && fl_e != capacity[e])
440            return false;
441        }
442      }
443      return true;
444    }
445   
446
447  }; //class MinCostFlow
448
449  ///@}
450
451} //namespace hugo
452
453#endif //HUGO_MINCOSTFLOW_H
Note: See TracBrowser for help on using the repository browser.