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
RevLine 
[610]1// -*- c++ -*-
[633]2#ifndef HUGO_MINCOSTFLOW_H
3#define HUGO_MINCOSTFLOW_H
[610]4
5///\ingroup galgs
6///\file
[645]7///\brief An algorithm for finding the minimum cost flow of given value in an uncapacitated network
[611]8
[610]9#include <hugo/dijkstra.h>
10#include <hugo/graph_wrapper.h>
11#include <hugo/maps.h>
12#include <vector>
[657]13#include <list>
[662]14#include <values.h>
[661]15#include <hugo/for_each_macros.h>
16#include <hugo/unionfind.h>
[662]17#include <hugo/bin_heap.h>
18#include <bfs_dfs.h>
[610]19
20namespace hugo {
21
22/// \addtogroup galgs
23/// @{
24
[661]25  ///\brief Implementation of an algorithm for solving the minimum cost general
26  /// flow problem in an uncapacitated network
[610]27  ///
28  ///
[633]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  ///
[661]36  /// The range of the cost (weight) function is nonnegative reals but
[610]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
[661]43  template <typename Graph, typename CostMap, typename SupplyDemandMap>
[633]44  class MinCostFlow {
[610]45
[661]46    typedef typename CostMap::ValueType Cost;
[610]47
[633]48
[635]49    typedef typename SupplyDemandMap::ValueType SupplyDemand;
[610]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;
[661]55    typedef typename Graph::template EdgeMap<SupplyDemand> FlowMap;
56    typedef ConstMap<Edge,SupplyDemand> ConstEdgeMap;
[610]57
58    //    typedef ConstMap<Edge,int> ConstMap;
59
[661]60    typedef ResGraphWrapper<const Graph,int,ConstEdgeMap,FlowMap> ResGraph;
61    typedef typename ResGraph::Edge ResGraphEdge;
[610]62
[661]63    class ModCostMap {   
64      //typedef typename ResGraph::template NodeMap<Cost> NodeMap;
65      typedef typename Graph::template NodeMap<Cost> NodeMap;
66      const ResGraph& res_graph;
[610]67      //      const EdgeIntMap& rev;
[661]68      const CostMap &ol;
[610]69      const NodeMap &pot;
70    public :
[661]71      typedef typename CostMap::KeyType KeyType;
72      typedef typename CostMap::ValueType ValueType;
[610]73       
[661]74      ValueType operator[](typename ResGraph::Edge e) const {     
[659]75        if (res_graph.forward(e))
76          return  ol[e]-(pot[res_graph.head(e)]-pot[res_graph.tail(e)]);   
[610]77        else
[659]78          return -ol[e]-(pot[res_graph.head(e)]-pot[res_graph.tail(e)]);   
[610]79      }     
80       
[661]81      ModCostMap(const ResGraph& _res_graph,
82                   const CostMap &o,  const NodeMap &p) :
[659]83        res_graph(_res_graph), /*rev(_rev),*/ ol(o), pot(p){};
[661]84    };//ModCostMap
[610]85
86
87  protected:
88   
89    //Input
[659]90    const Graph& graph;
[661]91    const CostMap& cost;
[635]92    const SupplyDemandMap& supply_demand;//supply or demand of nodes
[610]93
94
95    //auxiliary variables
96
97    //To store the flow
[661]98    FlowMap flow;
[662]99    //To store the potential (dual variables)
100    typedef typename Graph::template NodeMap<Cost> PotentialMap;
101    PotentialMap potential;
[633]102    //To store excess-deficit values
[635]103    SupplyDemandMap excess_deficit;
[610]104   
105
[661]106    Cost total_cost;
[610]107
108
109  public :
110
111
[662]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){ }
[610]119
120   
121    ///Runs the algorithm.
122
123    ///Runs the algorithm.
[635]124
[610]125    ///\todo May be it does make sense to be able to start with a nonzero
126    /// feasible primal-dual solution pair as well.
[659]127    void run() {
[610]128
129      //Resetting variables from previous runs
[661]130      //total_cost = 0;
[635]131
132      typedef typename Graph::template NodeMap<int> HeapMap;
[662]133      typedef BinHeap< Node, SupplyDemand, typename Graph::template NodeMap<int>,
[635]134        std::greater<SupplyDemand> >    HeapType;
135
136      //A heap for the excess nodes
[659]137      HeapMap excess_nodes_map(graph,-1);
[635]138      HeapType excess_nodes(excess_nodes_map);
139
140      //A heap for the deficit nodes
[659]141      HeapMap deficit_nodes_map(graph,-1);
[635]142      HeapType deficit_nodes(deficit_nodes_map);
143
[657]144      //A container to store nonabundant arcs
[662]145      std::list<Edge> nonabundant_arcs;
[659]146
147       
148      FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
[610]149        flow.set(e,0);
[657]150        nonabundant_arcs.push_back(e);
[610]151      }
[633]152
153      //Initial value for delta
[635]154      SupplyDemand delta = 0;
155
[657]156      typedef UnionFindEnum<Node, Graph::template NodeMap> UFE;
157
158      //A union-find structure to store the abundant components
[662]159      typename UFE::MapType abund_comp_map(graph);
[657]160      UFE abundant_components(abund_comp_map);
161
162
163
[659]164      FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
[635]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]);
[633]169        }
[635]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        }
[633]178        //Initialize the copy of the Dijkstra potential to zero
[610]179        potential.set(n,0);
[657]180        //Every single point is an abundant component initially
181        abundant_components.insert(n);
[610]182      }
183
[635]184      //It'll be allright as an initial value, though this value
185      //can be the maximum deficit here
186      SupplyDemand max_excess = delta;
[610]187     
[661]188      ///\bug This is a serious cheat here, before we have an uncapacitated ResGraph
[662]189      ConstEdgeMap const_inf_map(MAXINT);
[661]190     
[633]191      //We need a residual graph which is uncapacitated
[661]192      ResGraph res_graph(graph, const_inf_map, flow);
[659]193     
194      //An EdgeMap to tell which arcs are abundant
[662]195      typename Graph::template EdgeMap<bool> abundant_arcs(graph);
[610]196
[659]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);
[662]200      typedef SubGraphWrapper< const Graph, ConstNodeMap,
201         typename Graph::template EdgeMap<bool> >
[659]202        AbundantGraph;
203      AbundantGraph abundant_graph(graph, const_true_map, abundant_arcs );
204     
205      //Let's construct the residual graph for the abundant graph
[662]206      typedef ResGraphWrapper<const AbundantGraph,int,ConstEdgeMap,FlowMap>
[659]207        ResAbGraph;
208      //Again uncapacitated
[661]209      ResAbGraph res_ab_graph(abundant_graph, const_inf_map, flow);
[659]210     
211      //We need things for the bfs
[662]212      typename ResAbGraph::template NodeMap<bool> bfs_reached(res_ab_graph);
213      typename ResAbGraph::template NodeMap<typename ResAbGraph::Edge>
[659]214        bfs_pred(res_ab_graph);
[662]215      NullMap<typename ResAbGraph::Node, int> bfs_dist_dummy;
[659]216      //We want to run bfs-es (more) on this graph 'res_ab_graph'
217      Bfs < ResAbGraph ,
[662]218        typename ResAbGraph::template NodeMap<bool>,
219        typename ResAbGraph::template NodeMap<typename ResAbGraph::Edge>,
[659]220        NullMap<typename ResAbGraph::Node, int> >
221        bfs(res_ab_graph, bfs_reached, bfs_pred, bfs_dist_dummy);
[662]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       */
[610]231     
[661]232      ModCostMap mod_cost(res_graph, cost, potential);
[610]233
[661]234      Dijkstra<ResGraph, ModCostMap> dijkstra(res_graph, mod_cost);
[610]235
[633]236
[635]237      while (max_excess > 0){
238
[657]239        //Reset delta if still too big
240        if (8*number_of_nodes*max_excess <= delta){
241          delta = max_excess;
242         
243        }
244
[645]245        /*
246         * Beginning of the delta scaling phase
247        */
[635]248        //Merge and stuff
[657]249        {
250          SupplyDemand buf=8*number_of_nodes*delta;
[662]251          typename std::list<Edge>::iterator i = nonabundant_arcs.begin();
[657]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);
[659]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                }
[657]302              }
303              //What happens to i?
[659]304              //Marci and Zsolt says I shouldn't do such things
305              nonabundant_arcs.erase(i++);
306              abundant_arcs[i] = true;
[657]307            }
308            else
309              ++i;
310          }
311        }
312
[635]313
314        Node s = excess_nodes.top();
315        SupplyDemand max_excess = excess_nodes[s];
316        Node t = deficit_nodes.top();
[659]317        if (max_excess < deficit_nodes[t]){
318          max_excess = deficit_nodes[t];
[635]319        }
320
321
[662]322        while(max_excess > (number_of_nodes-1)*delta/number_of_nodes){
[659]323         
[635]324         
325          //s es t valasztasa
[659]326         
[635]327          //Dijkstra part       
328          dijkstra.run(s);
[659]329         
[635]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
[661]338          FOR_EACH_LOC(typename ResGraph::NodeIt, n, res_graph){
[635]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            /*
[661]351            //Let's update the total cost
[635]352            if (res_graph.forward(e))
[661]353              total_cost += cost[e];
[635]354            else
[661]355              total_cost -= cost[e];       
[635]356            */
357          }
[659]358         
359          //Update the excess_deficit map
360          excess_deficit[s] -= delta;
361          excess_deficit[t] += delta;
362         
[635]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
[659]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
[633]394        }
395
396        /*
[635]397         * End of the delta scaling phase
398        */
[633]399
[635]400        //Whatever this means
401        delta = delta / 2;
402
403        /*This is not necessary here
404        //Update the max_excess
405        max_excess = 0;
[659]406        FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
[635]407          if (max_excess < excess_deficit[n]){
408            max_excess = excess_deficit[n];
[610]409          }
410        }
[633]411        */
[657]412
[610]413         
[635]414      }//while(max_excess > 0)
[610]415     
416
417      return i;
418    }
419
420
421
422
[661]423    ///This function gives back the total cost of the found paths.
[610]424    ///Assumes that \c run() has been run and nothing changed since then.
[661]425    Cost totalCost(){
426      return total_cost;
[610]427    }
428
429    ///Returns a const reference to the EdgeMap \c flow. \pre \ref run() must
430    ///be called before using this function.
[662]431    const FlowMap &getFlow() const { return flow;}
[610]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.
[662]435    const PotentialMap &getPotential() const { return potential;}
[610]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(){
[661]443      Cost mod_pot;
444      Cost fl_e;
[659]445      FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
[610]446        //C^{\Pi}_{i,j}
[661]447        mod_pot = cost[e]-potential[graph.head(e)]+potential[graph.tail(e)];
[610]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
[633]465  }; //class MinCostFlow
[610]466
467  ///@}
468
469} //namespace hugo
470
471#endif //HUGO_MINCOSTFLOW_H
Note: See TracBrowser for help on using the repository browser.