COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/athos/mincostflow.h @ 662:0155001b6f65

Last change on this file since 662:0155001b6f65 was 662:0155001b6f65, checked in by athos, 17 years ago

Almost compiles.

File size: 12.8 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      //We want to run bfs-es (more) on this graph 'res_ab_graph'
217      Bfs < ResAbGraph ,
218        typename ResAbGraph::template NodeMap<bool>,
219        typename ResAbGraph::template NodeMap<typename ResAbGraph::Edge>,
220        NullMap<typename ResAbGraph::Node, int> >
221        bfs(res_ab_graph, bfs_reached, bfs_pred, bfs_dist_dummy);
222      /*This is what Marci wants for a bfs
223        template <typename Graph,
224            typename ReachedMap=typename Graph::template NodeMap<bool>,
225            typename PredMap
226            =typename Graph::template NodeMap<typename Graph::Edge>,
227            typename DistMap=typename Graph::template NodeMap<int> >
228            class Bfs : public BfsIterator<Graph, ReachedMap> {
229
230       */
231     
232      ModCostMap mod_cost(res_graph, cost, potential);
233
234      Dijkstra<ResGraph, ModCostMap> dijkstra(res_graph, mod_cost);
235
236
237      while (max_excess > 0){
238
239        //Reset delta if still too big
240        if (8*number_of_nodes*max_excess <= delta){
241          delta = max_excess;
242         
243        }
244
245        /*
246         * Beginning of the delta scaling phase
247        */
248        //Merge and stuff
249        {
250          SupplyDemand buf=8*number_of_nodes*delta;
251          typename std::list<Edge>::iterator i = nonabundant_arcs.begin();
252          while ( i != nonabundant_arcs.end() ){
253            if (flow[i]>=buf){
254              Node a = abundant_components.find(res_graph.head(i));
255              Node b = abundant_components.find(res_graph.tail(i));
256              //Merge
257              if (a != b){
258                abundant_components.join(a,b);
259                //We want to push the smaller
260                //Which has greater absolut value excess/deficit
261                Node root=(abs(excess_deficit[a])>abs(excess_deficit[b]))?a:b;
262                //Which is the other
263                Node non_root = ( a == root ) ? b : a ;
264                abundant_components.makeRep(root);
265                SupplyDemand qty_to_augment = abs(excess_deficit[non_root]);
266                //Push the positive value
267                if (excess_deficit[non_root] < 0)
268                  swap(root, non_root);
269                //If the non_root node has excess/deficit at all
270                if (qty_to_augment>0){
271                  //Find path and augment
272                  bfs.run(non_root);
273                  //root should be reached
274                 
275                  //Augmenting on the found path
276                  Node n=root;
277                  ResGraphEdge e;
278                  while (n!=non_root){
279                    e = bfs_pred(n);
280                    n = res_graph.tail(e);
281                    res_graph.augment(e,qty_to_augment);
282                  }
283         
284                  //We know that non_root had positive excess
285                  excess_nodes[non_root] -= qty_to_augment;
286                  //But what about root node
287                  //It might have been positive and so became larger
288                  if (excess_deficit[root]>0){
289                    excess_nodes[root] += qty_to_augment;
290                  }
291                  else{
292                    //Or negative but not turned into positive
293                    deficit_nodes[root] -= qty_to_augment;
294                  }
295
296                  //Update the excess_deficit map
297                  excess_deficit[non_root] -= qty_to_augment;
298                  excess_deficit[root] += qty_to_augment;
299
300                 
301                }
302              }
303              //What happens to i?
304              //Marci and Zsolt says I shouldn't do such things
305              nonabundant_arcs.erase(i++);
306              abundant_arcs[i] = true;
307            }
308            else
309              ++i;
310          }
311        }
312
313
314        Node s = excess_nodes.top();
315        SupplyDemand max_excess = excess_nodes[s];
316        Node t = deficit_nodes.top();
317        if (max_excess < deficit_nodes[t]){
318          max_excess = deficit_nodes[t];
319        }
320
321
322        while(max_excess > (number_of_nodes-1)*delta/number_of_nodes){
323         
324         
325          //s es t valasztasa
326         
327          //Dijkstra part       
328          dijkstra.run(s);
329         
330          /*We know from theory that t can be reached
331          if (!dijkstra.reached(t)){
332            //There are no k paths from s to t
333            break;
334          };
335          */
336         
337          //We have to change the potential
338          FOR_EACH_LOC(typename ResGraph::NodeIt, n, res_graph){
339            potential[n] += dijkstra.distMap()[n];
340          }
341
342
343          //Augmenting on the sortest path
344          Node n=t;
345          ResGraphEdge e;
346          while (n!=s){
347            e = dijkstra.pred(n);
348            n = dijkstra.predNode(n);
349            res_graph.augment(e,delta);
350            /*
351            //Let's update the total cost
352            if (res_graph.forward(e))
353              total_cost += cost[e];
354            else
355              total_cost -= cost[e];       
356            */
357          }
358         
359          //Update the excess_deficit map
360          excess_deficit[s] -= delta;
361          excess_deficit[t] += delta;
362         
363
364          //Update the excess_nodes heap
365          if (delta >= excess_nodes[s]){
366            if (delta > excess_nodes[s])
367              deficit_nodes.push(s,delta - excess_nodes[s]);
368            excess_nodes.pop();
369           
370          }
371          else{
372            excess_nodes[s] -= delta;
373          }
374          //Update the deficit_nodes heap
375          if (delta >= deficit_nodes[t]){
376            if (delta > deficit_nodes[t])
377              excess_nodes.push(t,delta - deficit_nodes[t]);
378            deficit_nodes.pop();
379           
380          }
381          else{
382            deficit_nodes[t] -= delta;
383          }
384          //Dijkstra part ends here
385         
386          //Choose s and t again
387          s = excess_nodes.top();
388          max_excess = excess_nodes[s];
389          t = deficit_nodes.top();
390          if (max_excess < deficit_nodes[t]){
391            max_excess = deficit_nodes[t];
392          }
393
394        }
395
396        /*
397         * End of the delta scaling phase
398        */
399
400        //Whatever this means
401        delta = delta / 2;
402
403        /*This is not necessary here
404        //Update the max_excess
405        max_excess = 0;
406        FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
407          if (max_excess < excess_deficit[n]){
408            max_excess = excess_deficit[n];
409          }
410        }
411        */
412
413         
414      }//while(max_excess > 0)
415     
416
417      return i;
418    }
419
420
421
422
423    ///This function gives back the total cost of the found paths.
424    ///Assumes that \c run() has been run and nothing changed since then.
425    Cost totalCost(){
426      return total_cost;
427    }
428
429    ///Returns a const reference to the EdgeMap \c flow. \pre \ref run() must
430    ///be called before using this function.
431    const FlowMap &getFlow() const { return flow;}
432
433  ///Returns a const reference to the NodeMap \c potential (the dual solution).
434    /// \pre \ref run() must be called before using this function.
435    const PotentialMap &getPotential() const { return potential;}
436
437    ///This function checks, whether the given solution is optimal
438    ///Running after a \c run() should return with true
439    ///In this "state of the art" this only check optimality, doesn't bother with feasibility
440    ///
441    ///\todo Is this OK here?
442    bool checkComplementarySlackness(){
443      Cost mod_pot;
444      Cost fl_e;
445      FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
446        //C^{\Pi}_{i,j}
447        mod_pot = cost[e]-potential[graph.head(e)]+potential[graph.tail(e)];
448        fl_e = flow[e];
449        //      std::cout << fl_e << std::endl;
450        if (0<fl_e && fl_e<capacity[e]){
451          if (mod_pot != 0)
452            return false;
453        }
454        else{
455          if (mod_pot > 0 && fl_e != 0)
456            return false;
457          if (mod_pot < 0 && fl_e != capacity[e])
458            return false;
459        }
460      }
461      return true;
462    }
463   
464
465  }; //class MinCostFlow
466
467  ///@}
468
469} //namespace hugo
470
471#endif //HUGO_MINCOSTFLOW_H
Note: See TracBrowser for help on using the repository browser.