COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/athos/mincostflow.h @ 761:58243a389464

Last change on this file since 761:58243a389464 was 672:6c7bd0edd1d7, checked in by athos, 20 years ago

Seems to work. More tests required.

File size: 14.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 <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   
103
104    Cost total_cost;
105
106
107  public :
108
109
110   MinCostFlow(Graph& _graph, CostMap& _cost, SupplyDemandMap& _supply_demand):
111     graph(_graph),
112     cost(_cost),
113     supply_demand(_supply_demand),
114     flow(_graph),
115     potential(_graph){ }
116
117   
118    ///Runs the algorithm.
119
120    ///Runs the algorithm.
121
122    ///\todo May be it does make sense to be able to start with a nonzero
123    /// feasible primal-dual solution pair as well.
124    void run() {
125
126      //To store excess-deficit values
127      SupplyDemandMap excess_deficit(graph);
128
129      //Resetting variables from previous runs
130      //total_cost = 0;
131
132
133      typedef typename Graph::template NodeMap<int> HeapMap;
134      typedef BinHeap< Node, SupplyDemand, typename Graph::template NodeMap<int>,
135        std::greater<SupplyDemand> >    HeapType;
136
137      //A heap for the excess nodes
138      HeapMap excess_nodes_map(graph,-1);
139      HeapType excess_nodes(excess_nodes_map);
140
141      //A heap for the deficit nodes
142      HeapMap deficit_nodes_map(graph,-1);
143      HeapType deficit_nodes(deficit_nodes_map);
144
145      //A container to store nonabundant arcs
146      std::list<Edge> nonabundant_arcs;
147
148       
149      FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
150        flow.set(e,0);
151        nonabundant_arcs.push_back(e);
152      }
153
154      //Initial value for delta
155      SupplyDemand delta = 0;
156
157      typedef UnionFindEnum<Node, Graph::template NodeMap> UFE;
158
159      //A union-find structure to store the abundant components
160      typename UFE::MapType abund_comp_map(graph);
161      UFE abundant_components(abund_comp_map);
162
163
164
165      FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
166        excess_deficit.set(n,supply_demand[n]);
167        //A supply node
168        if (excess_deficit[n] > 0){
169          excess_nodes.push(n,excess_deficit[n]);
170        }
171        //A demand node
172        if (excess_deficit[n] < 0){
173          deficit_nodes.push(n, - excess_deficit[n]);
174        }
175        //Finding out starting value of delta
176        if (delta < abs(excess_deficit[n])){
177          delta = abs(excess_deficit[n]);
178        }
179        //Initialize the copy of the Dijkstra potential to zero
180        potential.set(n,0);
181        //Every single point is an abundant component initially
182        abundant_components.insert(n);
183      }
184
185      //It'll be allright as an initial value, though this value
186      //can be the maximum deficit here
187      SupplyDemand max_excess = delta;
188     
189      ///\bug This is a serious cheat here, before we have an uncapacitated ResGraph
190      ConstEdgeMap const_inf_map(MAXINT);
191     
192      //We need a residual graph which is uncapacitated
193      ResGraph res_graph(graph, const_inf_map, flow);
194     
195      //An EdgeMap to tell which arcs are abundant
196      typename Graph::template EdgeMap<bool> abundant_arcs(graph);
197
198      //Let's construct the sugraph consisting only of the abundant edges
199      typedef ConstMap< typename Graph::Node, bool > ConstNodeMap;
200
201      ConstNodeMap const_true_map(true);
202      typedef SubGraphWrapper< const Graph, ConstNodeMap,
203         typename Graph::template EdgeMap<bool> >
204        AbundantGraph;
205      AbundantGraph abundant_graph(graph, const_true_map, abundant_arcs );
206     
207      //Let's construct the residual graph for the abundant graph
208      typedef ResGraphWrapper<const AbundantGraph,int,ConstEdgeMap,FlowMap>
209        ResAbGraph;
210      //Again uncapacitated
211      ResAbGraph res_ab_graph(abundant_graph, const_inf_map, flow);
212     
213      //We need things for the bfs
214      typename ResAbGraph::template NodeMap<bool> bfs_reached(res_ab_graph);
215      typename ResAbGraph::template NodeMap<typename ResAbGraph::Edge>
216        bfs_pred(res_ab_graph);
217      NullMap<typename ResAbGraph::Node, int> bfs_dist_dummy;
218      //Teszt celbol:
219      //BfsIterator<ResAbGraph, typename ResAbGraph::template NodeMap<bool> >
220      //izebize(res_ab_graph, bfs_reached);
221      //We want to run bfs-es (more) on this graph 'res_ab_graph'
222      Bfs < const ResAbGraph ,
223        typename ResAbGraph::template NodeMap<bool>,
224        typename ResAbGraph::template NodeMap<typename ResAbGraph::Edge>,
225        NullMap<typename ResAbGraph::Node, int> >
226        bfs(res_ab_graph, bfs_reached, bfs_pred, bfs_dist_dummy);
227      /*This is what Marci wants for a bfs
228        template <typename Graph,
229            typename ReachedMap=typename Graph::template NodeMap<bool>,
230            typename PredMap
231            =typename Graph::template NodeMap<typename Graph::Edge>,
232            typename DistMap=typename Graph::template NodeMap<int> >
233            class Bfs : public BfsIterator<Graph, ReachedMap> {
234
235       */
236     
237      ModCostMap mod_cost(res_graph, cost, potential);
238
239      Dijkstra<ResGraph, ModCostMap> dijkstra(res_graph, mod_cost);
240
241      //We will use the number of the nodes of the graph often
242      int number_of_nodes = graph.nodeNum();
243
244      while (max_excess > 0){
245
246        //Reset delta if still too big
247        if (8*number_of_nodes*max_excess <= delta){
248          delta = max_excess;
249         
250        }
251
252        /*
253         * Beginning of the delta scaling phase
254        */
255        //Merge and stuff
256        {
257          SupplyDemand buf=8*number_of_nodes*delta;
258          typename std::list<Edge>::iterator i = nonabundant_arcs.begin();
259          while ( i != nonabundant_arcs.end() ){
260            if (flow[*i]>=buf){
261              Node a = abundant_components.find(res_graph.head(*i));
262              Node b = abundant_components.find(res_graph.tail(*i));
263              //Merge
264              if (a != b){
265                abundant_components.join(a,b);
266                //We want to push the smaller
267                //Which has greater absolut value excess/deficit
268                Node root=(abs(excess_deficit[a])>abs(excess_deficit[b]))?a:b;
269                //Which is the other
270                Node non_root = ( a == root ) ? b : a ;
271                abundant_components.makeRep(root);
272                SupplyDemand qty_to_augment = abs(excess_deficit[non_root]);
273                //Push the positive value
274                if (excess_deficit[non_root] < 0)
275                  swap(root, non_root);
276                //If the non_root node has excess/deficit at all
277                if (qty_to_augment>0){
278                  //Find path and augment
279                  bfs.run(typename AbundantGraph::Node(non_root));
280                  //root should be reached
281                 
282                  //Augmenting on the found path
283                  Node n=root;
284                  ResGraphEdge e;
285                  while (n!=non_root){
286                    e = bfs_pred[n];
287                    n = res_graph.tail(e);
288                    res_graph.augment(e,qty_to_augment);
289                  }
290         
291                  //We know that non_root had positive excess
292                  excess_nodes.set(non_root,
293                                   excess_nodes[non_root] - qty_to_augment);
294                  //But what about root node
295                  //It might have been positive and so became larger
296                  if (excess_deficit[root]>0){
297                    excess_nodes.set(root,
298                                     excess_nodes[root] + qty_to_augment);
299                  }
300                  else{
301                    //Or negative but not turned into positive
302                    deficit_nodes.set(root,
303                                      deficit_nodes[root] - qty_to_augment);
304                  }
305
306                  //Update the excess_deficit map
307                  excess_deficit[non_root] -= qty_to_augment;
308                  excess_deficit[root] += qty_to_augment;
309
310                 
311                }
312              }
313              //What happens to i?
314              //Marci and Zsolt says I shouldn't do such things
315              nonabundant_arcs.erase(i++);
316              abundant_arcs[*i] = true;
317            }
318            else
319              ++i;
320          }
321        }
322
323
324        Node s = excess_nodes.top();
325        max_excess = excess_nodes[s];
326        Node t = deficit_nodes.top();
327        if (max_excess < deficit_nodes[t]){
328          max_excess = deficit_nodes[t];
329        }
330
331
332        while(max_excess > (number_of_nodes-1)*delta/number_of_nodes){
333         
334         
335          //s es t valasztasa
336         
337          //Dijkstra part       
338          dijkstra.run(s);
339         
340          /*We know from theory that t can be reached
341          if (!dijkstra.reached(t)){
342            //There are no k paths from s to t
343            break;
344          };
345          */
346         
347          //We have to change the potential
348          FOR_EACH_LOC(typename ResGraph::NodeIt, n, res_graph){
349            potential[n] += dijkstra.distMap()[n];
350          }
351
352
353          //Augmenting on the sortest path
354          Node n=t;
355          ResGraphEdge e;
356          while (n!=s){
357            e = dijkstra.pred(n);
358            n = dijkstra.predNode(n);
359            res_graph.augment(e,delta);
360            /*
361            //Let's update the total cost
362            if (res_graph.forward(e))
363              total_cost += cost[e];
364            else
365              total_cost -= cost[e];       
366            */
367          }
368         
369          //Update the excess_deficit map
370          excess_deficit[s] -= delta;
371          excess_deficit[t] += delta;
372         
373
374          //Update the excess_nodes heap
375          if (delta > excess_nodes[s]){
376            if (delta > excess_nodes[s])
377              deficit_nodes.push(s,delta - excess_nodes[s]);
378            excess_nodes.pop();
379           
380          }
381          else{
382            excess_nodes.set(s, excess_nodes[s] - delta);
383          }
384          //Update the deficit_nodes heap
385          if (delta > deficit_nodes[t]){
386            if (delta > deficit_nodes[t])
387              excess_nodes.push(t,delta - deficit_nodes[t]);
388            deficit_nodes.pop();
389           
390          }
391          else{
392            deficit_nodes.set(t, deficit_nodes[t] - delta);
393          }
394          //Dijkstra part ends here
395         
396          //Choose s and t again
397          s = excess_nodes.top();
398          max_excess = excess_nodes[s];
399          t = deficit_nodes.top();
400          if (max_excess < deficit_nodes[t]){
401            max_excess = deficit_nodes[t];
402          }
403
404        }
405
406        /*
407         * End of the delta scaling phase
408        */
409
410        //Whatever this means
411        delta = delta / 2;
412
413        /*This is not necessary here
414        //Update the max_excess
415        max_excess = 0;
416        FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
417          if (max_excess < excess_deficit[n]){
418            max_excess = excess_deficit[n];
419          }
420        }
421        */
422
423         
424      }//while(max_excess > 0)
425     
426
427      //return i;
428    }
429
430
431
432
433    ///This function gives back the total cost of the found paths.
434    ///Assumes that \c run() has been run and nothing changed since then.
435    Cost totalCost(){
436      return total_cost;
437    }
438
439    ///Returns a const reference to the EdgeMap \c flow. \pre \ref run() must
440    ///be called before using this function.
441    const FlowMap &getFlow() const { return flow;}
442
443  ///Returns a const reference to the NodeMap \c potential (the dual solution).
444    /// \pre \ref run() must be called before using this function.
445    const PotentialMap &getPotential() const { return potential;}
446
447    ///This function checks, whether the given solution is optimal
448    ///Running after a \c run() should return with true
449    ///In this "state of the art" this only checks optimality, doesn't bother with feasibility
450    ///
451    ///\todo Is this OK here?
452    bool checkComplementarySlackness(){
453      Cost mod_pot;
454      Cost fl_e;
455      FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
456        //C^{\Pi}_{i,j}
457        mod_pot = cost[e]-potential[graph.head(e)]+potential[graph.tail(e)];
458        fl_e = flow[e];
459        //      std::cout << fl_e << std::endl;
460        if (mod_pot > 0 && fl_e != 0)
461          return false;
462
463      }
464      return true;
465    }
466
467    /*
468    //For testing purposes only
469    //Lists the node_properties
470    void write_property_vector(const SupplyDemandMap& a,
471                               char* prop_name="property"){
472      FOR_EACH_LOC(typename Graph::NodeIt, i, graph){
473        cout<<"Node id.: "<<graph.id(i)<<", "<<prop_name<<" value: "<<a[i]<<endl;
474      }
475      cout<<endl;
476    }
477    */
478    bool checkFeasibility(){
479      SupplyDemandMap supdem(graph);
480      FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
481
482        if ( flow[e] < 0){
483
484          return false;
485        }
486        supdem[graph.tail(e)] += flow[e];
487        supdem[graph.head(e)] -= flow[e];
488      }
489      //write_property_vector(supdem, "supdem");
490      //write_property_vector(supply_demand, "supply_demand");
491
492      FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
493
494        if ( supdem[n] != supply_demand[n]){
495          //cout<<"Node id.: "<<graph.id(n)<<" : "<<supdem[n]<<", should be: "<<supply_demand[n]<<endl;
496          return false;
497        }
498      }
499
500      return true;
501    }
502
503    bool checkOptimality(){
504      return checkFeasibility() && checkComplementarySlackness();
505    }
506
507  }; //class MinCostFlow
508
509  ///@}
510
511} //namespace hugo
512
513#endif //HUGO_MINCOSTFLOW_H
Note: See TracBrowser for help on using the repository browser.