# source:lemon-0.x/src/work/athos/mincostflow.h@659:c5984e925384

Last change on this file since 659:c5984e925384 was 659:c5984e925384, checked in by athos, 20 years ago

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