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
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>
[661]14#include <hugo/for_each_macros.h>
15#include <hugo/unionfind.h>
[610]16
17namespace hugo {
18
19/// \addtogroup galgs
20/// @{
21
[661]22  ///\brief Implementation of an algorithm for solving the minimum cost general
23  /// flow problem in an uncapacitated network
[610]24  ///
25  ///
[633]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  ///
[661]33  /// The range of the cost (weight) function is nonnegative reals but
[610]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
[661]40  template <typename Graph, typename CostMap, typename SupplyDemandMap>
[633]41  class MinCostFlow {
[610]42
[661]43    typedef typename CostMap::ValueType Cost;
[610]44
[633]45
[635]46    typedef typename SupplyDemandMap::ValueType SupplyDemand;
[610]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;
[661]52    typedef typename Graph::template EdgeMap<SupplyDemand> FlowMap;
53    typedef ConstMap<Edge,SupplyDemand> ConstEdgeMap;
[610]54
55    //    typedef ConstMap<Edge,int> ConstMap;
56
[661]57    typedef ResGraphWrapper<const Graph,int,ConstEdgeMap,FlowMap> ResGraph;
58    typedef typename ResGraph::Edge ResGraphEdge;
[610]59
[661]60    class ModCostMap {   
61      //typedef typename ResGraph::template NodeMap<Cost> NodeMap;
62      typedef typename Graph::template NodeMap<Cost> NodeMap;
63      const ResGraph& res_graph;
[610]64      //      const EdgeIntMap& rev;
[661]65      const CostMap &ol;
[610]66      const NodeMap &pot;
67    public :
[661]68      typedef typename CostMap::KeyType KeyType;
69      typedef typename CostMap::ValueType ValueType;
[610]70       
[661]71      ValueType operator[](typename ResGraph::Edge e) const {     
[659]72        if (res_graph.forward(e))
73          return  ol[e]-(pot[res_graph.head(e)]-pot[res_graph.tail(e)]);   
[610]74        else
[659]75          return -ol[e]-(pot[res_graph.head(e)]-pot[res_graph.tail(e)]);   
[610]76      }     
77       
[661]78      ModCostMap(const ResGraph& _res_graph,
79                   const CostMap &o,  const NodeMap &p) :
[659]80        res_graph(_res_graph), /*rev(_rev),*/ ol(o), pot(p){};
[661]81    };//ModCostMap
[610]82
83
84  protected:
85   
86    //Input
[659]87    const Graph& graph;
[661]88    const CostMap& cost;
[635]89    const SupplyDemandMap& supply_demand;//supply or demand of nodes
[610]90
91
92    //auxiliary variables
93
94    //To store the flow
[661]95    FlowMap flow;
[610]96    //To store the potentila (dual variables)
[661]97    typename Graph::template NodeMap<Cost> potential;
[633]98    //To store excess-deficit values
[635]99    SupplyDemandMap excess_deficit;
[610]100   
101
[661]102    Cost total_cost;
[610]103
104
105  public :
106
107
[661]108    MinCostFlow(Graph& _graph, CostMap& _cost, SupplyDemandMap& _supply_demand) : graph(_graph),
109      cost(_cost), supply_demand(_supply_demand), flow(_graph), potential(_graph){ }
[610]110
111   
112    ///Runs the algorithm.
113
114    ///Runs the algorithm.
[635]115
[610]116    ///\todo May be it does make sense to be able to start with a nonzero
117    /// feasible primal-dual solution pair as well.
[659]118    void run() {
[610]119
120      //Resetting variables from previous runs
[661]121      //total_cost = 0;
[635]122
123      typedef typename Graph::template NodeMap<int> HeapMap;
[657]124      typedef Heap< Node, SupplyDemand, typename Graph::template NodeMap<int>,
[635]125        std::greater<SupplyDemand> >    HeapType;
126
127      //A heap for the excess nodes
[659]128      HeapMap excess_nodes_map(graph,-1);
[635]129      HeapType excess_nodes(excess_nodes_map);
130
131      //A heap for the deficit nodes
[659]132      HeapMap deficit_nodes_map(graph,-1);
[635]133      HeapType deficit_nodes(deficit_nodes_map);
134
[657]135      //A container to store nonabundant arcs
136      list<Edge> nonabundant_arcs;
[659]137
138       
139      FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
[610]140        flow.set(e,0);
[657]141        nonabundant_arcs.push_back(e);
[610]142      }
[633]143
144      //Initial value for delta
[635]145      SupplyDemand delta = 0;
146
[657]147      typedef UnionFindEnum<Node, Graph::template NodeMap> UFE;
148
149      //A union-find structure to store the abundant components
[659]150      UFE::MapType abund_comp_map(graph);
[657]151      UFE abundant_components(abund_comp_map);
152
153
154
[659]155      FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
[635]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]);
[633]160        }
[635]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        }
[633]169        //Initialize the copy of the Dijkstra potential to zero
[610]170        potential.set(n,0);
[657]171        //Every single point is an abundant component initially
172        abundant_components.insert(n);
[610]173      }
174
[635]175      //It'll be allright as an initial value, though this value
176      //can be the maximum deficit here
177      SupplyDemand max_excess = delta;
[610]178     
[661]179      ///\bug This is a serious cheat here, before we have an uncapacitated ResGraph
180      ConstEdgeMap const_inf_map(MAX_INT);
181     
[633]182      //We need a residual graph which is uncapacitated
[661]183      ResGraph res_graph(graph, const_inf_map, flow);
[659]184     
185      //An EdgeMap to tell which arcs are abundant
186      template typename Graph::EdgeMap<bool> abundant_arcs(graph);
[610]187
[659]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
[661]200      ResAbGraph res_ab_graph(abundant_graph, const_inf_map, flow);
[659]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);
[610]213     
[661]214      ModCostMap mod_cost(res_graph, cost, potential);
[610]215
[661]216      Dijkstra<ResGraph, ModCostMap> dijkstra(res_graph, mod_cost);
[610]217
[633]218
[635]219      while (max_excess > 0){
220
[657]221        //Reset delta if still too big
222        if (8*number_of_nodes*max_excess <= delta){
223          delta = max_excess;
224         
225        }
226
[645]227        /*
228         * Beginning of the delta scaling phase
229        */
[635]230        //Merge and stuff
[657]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);
[659]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                }
[657]284              }
285              //What happens to i?
[659]286              //Marci and Zsolt says I shouldn't do such things
287              nonabundant_arcs.erase(i++);
288              abundant_arcs[i] = true;
[657]289            }
290            else
291              ++i;
292          }
293        }
294
[635]295
296        Node s = excess_nodes.top();
297        SupplyDemand max_excess = excess_nodes[s];
298        Node t = deficit_nodes.top();
[659]299        if (max_excess < deficit_nodes[t]){
300          max_excess = deficit_nodes[t];
[635]301        }
302
303
[659]304        while(max_excess > (n-1)*delta/n){
305         
[635]306         
307          //s es t valasztasa
[659]308         
[635]309          //Dijkstra part       
310          dijkstra.run(s);
[659]311         
[635]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
[661]320          FOR_EACH_LOC(typename ResGraph::NodeIt, n, res_graph){
[635]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            /*
[661]333            //Let's update the total cost
[635]334            if (res_graph.forward(e))
[661]335              total_cost += cost[e];
[635]336            else
[661]337              total_cost -= cost[e];       
[635]338            */
339          }
[659]340         
341          //Update the excess_deficit map
342          excess_deficit[s] -= delta;
343          excess_deficit[t] += delta;
344         
[635]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
[659]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
[633]376        }
377
378        /*
[635]379         * End of the delta scaling phase
380        */
[633]381
[635]382        //Whatever this means
383        delta = delta / 2;
384
385        /*This is not necessary here
386        //Update the max_excess
387        max_excess = 0;
[659]388        FOR_EACH_LOC(typename Graph::NodeIt, n, graph){
[635]389          if (max_excess < excess_deficit[n]){
390            max_excess = excess_deficit[n];
[610]391          }
392        }
[633]393        */
[657]394
[610]395         
[635]396      }//while(max_excess > 0)
[610]397     
398
399      return i;
400    }
401
402
403
404
[661]405    ///This function gives back the total cost of the found paths.
[610]406    ///Assumes that \c run() has been run and nothing changed since then.
[661]407    Cost totalCost(){
408      return total_cost;
[610]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(){
[661]425      Cost mod_pot;
426      Cost fl_e;
[659]427      FOR_EACH_LOC(typename Graph::EdgeIt, e, graph){
[610]428        //C^{\Pi}_{i,j}
[661]429        mod_pot = cost[e]-potential[graph.head(e)]+potential[graph.tail(e)];
[610]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
[633]447  }; //class MinCostFlow
[610]448
449  ///@}
450
451} //namespace hugo
452
453#endif //HUGO_MINCOSTFLOW_H
Note: See TracBrowser for help on using the repository browser.