COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/athos/mincostflow.h @ 671:708df4dc6ab6

Last change on this file since 671:708df4dc6ab6 was 671:708df4dc6ab6, checked in by athos, 17 years ago

Compiles now

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