COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/marci/augmenting_flow.h @ 862:732f2acb7239

Last change on this file since 862:732f2acb7239 was 862:732f2acb7239, checked in by marci, 17 years ago

bug correction

File size: 19.1 KB
Line 
1// -*- C++ -*-
2#ifndef HUGO_AUGMENTING_FLOW_H
3#define HUGO_AUGMENTING_FLOW_H
4
5#include <vector>
6//#include <queue>
7//#include <stack>
8#include <iostream>
9
10#include <hugo/graph_wrapper.h>
11#include <bfs_dfs.h>
12#include <hugo/invalid.h>
13#include <hugo/maps.h>
14
15/// \file
16/// \brief Maximum flow algorithms.
17/// \ingroup galgs
18
19namespace hugo {
20
21  /// \brief A map for filtering the edge-set to those edges
22  /// which are tight w.r.t. some node_potential map and
23  /// edge_distance map.
24  ///
25  /// A node-map node_potential is said to be a potential w.r.t.
26  /// an edge-map edge_distance
27  /// if and only if for each edge e, node_potential[g.head(e)]
28  /// <= edge_distance[e]+node_potential[g.tail(e)]
29  /// (or the reverse inequality holds for each edge).
30  /// An edge is said to be tight if this inequality holds with equality,
31  /// and the map returns true exactly for those edges.
32  /// To avoid rounding errors, it is recommended to use this class with exact
33  /// types, e.g. with int.
34  template<typename Graph,
35           typename NodePotentialMap, typename EdgeDistanceMap>
36  class TightEdgeFilterMap {
37  protected:
38    const Graph* g;
39    NodePotentialMap* node_potential;
40    EdgeDistanceMap* edge_distance;
41  public:
42    TightEdgeFilterMap(Graph& _g, NodePotentialMap& _node_potential,
43                       EdgeDistanceMap& _edge_distance) :
44      g(&_g), node_potential(&_node_potential),
45      edge_distance(&_edge_distance) { }
46//     void set(const typename Graph::Node& n, int a) {
47//       pot->set(n, a);
48//     }
49//     int operator[](const typename Graph::Node& n) const {
50//       return (*node_potential)[n];
51//     }
52    bool operator[](const typename Graph::Edge& e) const {
53      return ((*node_potential)[g->head(e)] ==
54              (*edge_distance)[e]+(*node_potential)[g->tail(e)]);
55    }
56  };
57
58  /// \addtogroup galgs
59  /// @{                                                                                                                                       
60  /// Class for augmenting path flow algorithms.
61
62  /// This class provides various algorithms for finding a flow of
63  /// maximum value in a directed graph. The \e source node, the \e
64  /// target node, the \e capacity of the edges and the \e starting \e
65  /// flow value of the edges should be passed to the algorithm through the
66  /// constructor.
67//   /// It is possible to change these quantities using the
68//   /// functions \ref resetSource, \ref resetTarget, \ref resetCap and
69//   /// \ref resetFlow. Before any subsequent runs of any algorithm of
70//   /// the class \ref resetFlow should be called.
71
72  /// After running an algorithm of the class, the actual flow value
73  /// can be obtained by calling \ref flowValue(). The minimum
74  /// value cut can be written into a \c node map of \c bools by
75  /// calling \ref minCut. (\ref minMinCut and \ref maxMinCut writes
76  /// the inclusionwise minimum and maximum of the minimum value
77  /// cuts, resp.)                                                                                                                               
78  ///\param Graph The directed graph type the algorithm runs on.
79  ///\param Num The number type of the capacities and the flow values.
80  ///\param CapMap The capacity map type.
81  ///\param FlowMap The flow map type.                                                                                                           
82  ///\author Marton Makai
83  template <typename Graph, typename Num,
84            typename CapMap=typename Graph::template EdgeMap<Num>,
85            typename FlowMap=typename Graph::template EdgeMap<Num> >
86  class AugmentingFlow {
87  protected:
88    typedef typename Graph::Node Node;
89    typedef typename Graph::NodeIt NodeIt;
90    typedef typename Graph::EdgeIt EdgeIt;
91    typedef typename Graph::OutEdgeIt OutEdgeIt;
92    typedef typename Graph::InEdgeIt InEdgeIt;
93
94    const Graph* g;
95    Node s;
96    Node t;
97    const CapMap* capacity;
98    FlowMap* flow;
99//    int n;      //the number of nodes of G
100    typedef ResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;   
101    //typedef ExpResGraphWrapper<const Graph, Num, CapMap, FlowMap> ResGW;
102    typedef typename ResGW::OutEdgeIt ResGWOutEdgeIt;
103    typedef typename ResGW::Edge ResGWEdge;
104    //typedef typename ResGW::template NodeMap<bool> ReachedMap;
105    typedef typename Graph::template NodeMap<int> ReachedMap;
106
107    //level works as a bool map in augmenting path algorithms and is
108    //used by bfs for storing reached information.  In preflow, it
109    //shows the levels of nodes.     
110    ReachedMap level;
111
112  public:
113    ///Indicates the property of the starting flow.
114
115    ///Indicates the property of the starting flow. The meanings are as follows:
116    ///- \c ZERO_FLOW: constant zero flow
117    ///- \c GEN_FLOW: any flow, i.e. the sum of the in-flows equals to
118    ///the sum of the out-flows in every node except the \e source and
119    ///the \e target.
120    ///- \c PRE_FLOW: any preflow, i.e. the sum of the in-flows is at
121    ///least the sum of the out-flows in every node except the \e source.
122    ///- \c NO_FLOW: indicates an unspecified edge map. \ref flow will be
123    ///set to the constant zero flow in the beginning of the algorithm in this case.
124    enum FlowEnum{
125      ZERO_FLOW,
126      GEN_FLOW,
127      PRE_FLOW,
128      NO_FLOW
129    };
130
131    enum StatusEnum {
132      AFTER_NOTHING,
133      AFTER_AUGMENTING,
134      AFTER_FAST_AUGMENTING,
135      AFTER_PRE_FLOW_PHASE_1,     
136      AFTER_PRE_FLOW_PHASE_2
137    };
138
139    /// Don not needle this flag only if necessary.
140    StatusEnum status;
141    int number_of_augmentations;
142
143
144    template<typename IntMap>
145    class TrickyReachedMap {
146    protected:
147      IntMap* map;
148      int* number_of_augmentations;
149    public:
150      TrickyReachedMap(IntMap& _map, int& _number_of_augmentations) :
151        map(&_map), number_of_augmentations(&_number_of_augmentations) { }
152      void set(const Node& n, bool b) {
153        if (b)
154          map->set(n, *number_of_augmentations);
155        else
156          map->set(n, *number_of_augmentations-1);
157      }
158      bool operator[](const Node& n) const {
159        return (*map)[n]==*number_of_augmentations;
160      }
161    };
162   
163    AugmentingFlow(const Graph& _G, Node _s, Node _t, const CapMap& _capacity,
164                   FlowMap& _flow) :
165      g(&_G), s(_s), t(_t), capacity(&_capacity),
166      flow(&_flow), //n(_G.nodeNum()),
167      level(_G), //excess(_G,0),
168      status(AFTER_NOTHING), number_of_augmentations(0) { }
169
170    /// Starting from a flow, this method searches for an augmenting path
171    /// according to the Edmonds-Karp algorithm
172    /// and augments the flow on if any.
173    /// The return value shows if the augmentation was succesful.
174    bool augmentOnShortestPath();
175    bool augmentOnShortestPath2();
176
177    /// Starting from a flow, this method searches for an augmenting blocking
178    /// flow according to Dinits' algorithm and augments the flow on if any.
179    /// The blocking flow is computed in a physically constructed
180    /// residual graph of type \c Mutablegraph.
181    /// The return value show sif the augmentation was succesful.
182    template<typename MutableGraph> bool augmentOnBlockingFlow();
183
184    /// The same as \c augmentOnBlockingFlow<MutableGraph> but the
185    /// residual graph is not constructed physically.
186    /// The return value shows if the augmentation was succesful.
187    bool augmentOnBlockingFlow2();
188
189    template<typename _CutMap>
190    void actMinCut(_CutMap& M) const {
191      NodeIt v;
192      switch (status) {
193        case AFTER_PRE_FLOW_PHASE_1:
194//      std::cout << "AFTER_PRE_FLOW_PHASE_1" << std::endl;
195//      for(g->first(v); g->valid(v); g->next(v)) {
196//        if (level[v] < n) {
197//          M.set(v, false);
198//        } else {
199//          M.set(v, true);
200//        }
201//      }
202        break;
203      case AFTER_PRE_FLOW_PHASE_2:
204//      std::cout << "AFTER_PRE_FLOW_PHASE_2" << std::endl;
205        break;
206      case AFTER_NOTHING:
207//      std::cout << "AFTER_NOTHING" << std::endl;
208        minMinCut(M);
209        break;
210      case AFTER_AUGMENTING:
211//      std::cout << "AFTER_AUGMENTING" << std::endl;
212        for(g->first(v); v!=INVALID; ++v) {
213          if (level[v]) {
214            M.set(v, true);
215          } else {
216            M.set(v, false);
217          }
218        }
219        break;
220      case AFTER_FAST_AUGMENTING:
221//      std::cout << "AFTER_FAST_AUGMENTING" << std::endl;
222        for(g->first(v); v!=INVALID; ++v) {
223          if (level[v]==number_of_augmentations) {
224            M.set(v, true);
225          } else {
226            M.set(v, false);
227          }
228        }
229        break;
230      }
231    }
232
233    template<typename _CutMap>
234    void minMinCut(_CutMap& M) const {
235      std::queue<Node> queue;
236
237      M.set(s,true);
238      queue.push(s);
239
240      while (!queue.empty()) {
241        Node w=queue.front();
242        queue.pop();
243
244        OutEdgeIt e;
245        for(g->first(e,w) ; e!=INVALID; ++e) {
246          Node v=g->head(e);
247          if (!M[v] && (*flow)[e] < (*capacity)[e] ) {
248            queue.push(v);
249            M.set(v, true);
250          }
251        }
252
253        InEdgeIt f;
254        for(g->first(f,w) ; f!=INVALID; ++f) {
255          Node v=g->tail(f);
256          if (!M[v] && (*flow)[f] > 0 ) {
257            queue.push(v);
258            M.set(v, true);
259          }
260        }
261      }
262    }
263
264    template<typename _CutMap>
265    void minMinCut2(_CutMap& M) const {
266      ResGW res_graph(*g, *capacity, *flow);
267      BfsIterator<ResGW, _CutMap> bfs(res_graph, M);
268      bfs.pushAndSetReached(s);
269      while (!bfs.finished()) ++bfs;
270    }
271
272    Num flowValue() const {
273      Num a=0;
274      for (InEdgeIt e(*g, t); e!=INVALID; ++e) a+=(*flow)[e];
275      for (OutEdgeIt e(*g, t); e!=INVALID; ++e) a-=(*flow)[e];
276      return a;
277      //marci figyu: excess[t] epp ezt adja preflow 1. fazisa utan   
278    }
279
280  };
281
282
283
284  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
285  bool AugmentingFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath()
286  {
287    ResGW res_graph(*g, *capacity, *flow);
288    bool _augment=false;
289
290    //ReachedMap level(res_graph);
291    for (typename Graph::NodeIt n(*g); n!=INVALID; ++n) level.set(n, 0);
292    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
293    bfs.pushAndSetReached(s);
294
295    typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
296    pred.set(s, INVALID);
297
298    typename ResGW::template NodeMap<Num> free(res_graph);
299
300    //searching for augmenting path
301    while ( !bfs.finished() ) {
302      ResGWEdge e=bfs;
303      if (e!=INVALID && bfs.isBNodeNewlyReached()) {
304        Node v=res_graph.tail(e);
305        Node w=res_graph.head(e);
306        pred.set(w, e);
307        if (pred[v]!=INVALID) {
308          free.set(w, std::min(free[v], res_graph.resCap(e)));
309        } else {
310          free.set(w, res_graph.resCap(e));
311        }
312        if (res_graph.head(e)==t) { _augment=true; break; }
313      }
314
315      ++bfs;
316    } //end of searching augmenting path
317
318    if (_augment) {
319      Node n=t;
320      Num augment_value=free[t];
321      while (pred[n]!=INVALID) {
322        ResGWEdge e=pred[n];
323        res_graph.augment(e, augment_value);
324        n=res_graph.tail(e);
325      }
326    }
327
328    status=AFTER_AUGMENTING;
329    return _augment;
330  }
331
332  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
333  bool AugmentingFlow<Graph, Num, CapMap, FlowMap>::augmentOnShortestPath2()
334  {
335    ResGW res_graph(*g, *capacity, *flow);
336    bool _augment=false;
337
338    if (status!=AFTER_FAST_AUGMENTING) {
339      for (typename Graph::NodeIt n(*g); n!=INVALID; ++n) level.set(n, 0);
340      number_of_augmentations=1;
341    } else {
342      ++number_of_augmentations;
343    }
344    TrickyReachedMap<ReachedMap>
345      tricky_reached_map(level, number_of_augmentations);
346    //ReachedMap level(res_graph);
347//    FOR_EACH_LOC(typename Graph::NodeIt, e, *g) level.set(e, 0);
348    BfsIterator<ResGW, TrickyReachedMap<ReachedMap> >
349      bfs(res_graph, tricky_reached_map);
350    bfs.pushAndSetReached(s);
351
352    typename ResGW::template NodeMap<ResGWEdge> pred(res_graph);
353    pred.set(s, INVALID);
354
355    typename ResGW::template NodeMap<Num> free(res_graph);
356
357    //searching for augmenting path
358    while ( !bfs.finished() ) {
359      ResGWEdge e=bfs;
360      if (e!=INVALID && bfs.isBNodeNewlyReached()) {
361        Node v=res_graph.tail(e);
362        Node w=res_graph.head(e);
363        pred.set(w, e);
364        if (pred[v]!=INVALID) {
365          free.set(w, std::min(free[v], res_graph.resCap(e)));
366        } else {
367          free.set(w, res_graph.resCap(e));
368        }
369        if (res_graph.head(e)==t) { _augment=true; break; }
370      }
371
372      ++bfs;
373    } //end of searching augmenting path
374
375    if (_augment) {
376      Node n=t;
377      Num augment_value=free[t];
378      while (pred[n]!=INVALID) {
379        ResGWEdge e=pred[n];
380        res_graph.augment(e, augment_value);
381        n=res_graph.tail(e);
382      }
383    }
384
385    status=AFTER_FAST_AUGMENTING;
386    return _augment;
387  }
388
389
390  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
391  template<typename MutableGraph>
392  bool AugmentingFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow()
393  {
394    typedef MutableGraph MG;
395    bool _augment=false;
396
397    ResGW res_graph(*g, *capacity, *flow);
398
399    //bfs for distances on the residual graph
400    //ReachedMap level(res_graph);
401    for (typename Graph::NodeIt n(*g); n!=INVALID; ++n) level.set(n, 0);
402    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
403    bfs.pushAndSetReached(s);
404    typename ResGW::template NodeMap<int>
405      dist(res_graph); //filled up with 0's
406
407    //F will contain the physical copy of the residual graph
408    //with the set of edges which are on shortest paths
409    MG F;
410    typename ResGW::template NodeMap<typename MG::Node>
411      res_graph_to_F(res_graph);
412    {
413      typename ResGW::NodeIt n;
414      for(res_graph.first(n); n!=INVALID; ++n)
415        res_graph_to_F.set(n, F.addNode());
416    }
417
418    typename MG::Node sF=res_graph_to_F[s];
419    typename MG::Node tF=res_graph_to_F[t];
420    typename MG::template EdgeMap<ResGWEdge> original_edge(F);
421    typename MG::template EdgeMap<Num> residual_capacity(F);
422
423    while ( !bfs.finished() ) {
424      ResGWEdge e=bfs;
425      if (e!=INVALID) {
426        if (bfs.isBNodeNewlyReached()) {
427          dist.set(res_graph.head(e), dist[res_graph.tail(e)]+1);
428          typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)],
429                                        res_graph_to_F[res_graph.head(e)]);
430          //original_edge.update();
431          original_edge.set(f, e);
432          //residual_capacity.update();
433          residual_capacity.set(f, res_graph.resCap(e));
434        } else {
435          if (dist[res_graph.head(e)]==(dist[res_graph.tail(e)]+1)) {
436            typename MG::Edge f=F.addEdge(res_graph_to_F[res_graph.tail(e)],
437                                          res_graph_to_F[res_graph.head(e)]);
438            //original_edge.update();
439            original_edge.set(f, e);
440            //residual_capacity.update();
441            residual_capacity.set(f, res_graph.resCap(e));
442          }
443        }
444      }
445      ++bfs;
446    } //computing distances from s in the residual graph
447
448    bool __augment=true;
449
450    while (__augment) {
451      __augment=false;
452      //computing blocking flow with dfs
453      DfsIterator< MG, typename MG::template NodeMap<bool> > dfs(F);
454      typename MG::template NodeMap<typename MG::Edge> pred(F);
455      pred.set(sF, INVALID);
456      //invalid iterators for sources
457
458      typename MG::template NodeMap<Num> free(F);
459
460      dfs.pushAndSetReached(sF);
461      while (!dfs.finished()) {
462        ++dfs;
463        if (typename MG::Edge(dfs)!=INVALID) {
464          if (dfs.isBNodeNewlyReached()) {
465            typename MG::Node v=F.tail(dfs);
466            typename MG::Node w=F.head(dfs);
467            pred.set(w, dfs);
468            if (pred[v]!=INVALID) {
469              free.set(w, std::min(free[v], residual_capacity[dfs]));
470            } else {
471              free.set(w, residual_capacity[dfs]);
472            }
473            if (w==tF) {
474              __augment=true;
475              _augment=true;
476              break;
477            }
478
479          } else {
480            F.erase(typename MG::Edge(dfs));
481          }
482        }
483      }
484
485      if (__augment) {
486        typename MG::Node n=tF;
487        Num augment_value=free[tF];
488        while (pred[n]!=INVALID) {
489          typename MG::Edge e=pred[n];
490          res_graph.augment(original_edge[e], augment_value);
491          n=F.tail(e);
492          if (residual_capacity[e]==augment_value)
493            F.erase(e);
494          else
495            residual_capacity.set(e, residual_capacity[e]-augment_value);
496        }
497      }
498
499    }
500
501    status=AFTER_AUGMENTING;
502    return _augment;
503  }
504
505  /// Blocking flow augmentation without constructing the layered
506  /// graph physically in which the blocking flow is computed.
507  template <typename Graph, typename Num, typename CapMap, typename FlowMap>
508  bool AugmentingFlow<Graph, Num, CapMap, FlowMap>::augmentOnBlockingFlow2()
509  {
510    bool _augment=false;
511
512    ResGW res_graph(*g, *capacity, *flow);
513
514    //Potential map, for distances from s
515    typename ResGW::template NodeMap<int> potential(res_graph, 0);
516    typedef ConstMap<typename ResGW::Edge, int> Const1Map;
517    Const1Map const_1_map(1);
518    TightEdgeFilterMap<ResGW, typename ResGW::template NodeMap<int>,
519      Const1Map> tight_edge_filter(res_graph, potential, const_1_map);
520
521    for (typename Graph::NodeIt n(*g); n!=INVALID; ++n) level.set(n, 0);
522    BfsIterator<ResGW, ReachedMap> bfs(res_graph, level);
523    bfs.pushAndSetReached(s);
524
525    //computing distances from s in the residual graph
526    while ( !bfs.finished() ) {
527      ResGWEdge e=bfs;
528      if (e!=INVALID && bfs.isBNodeNewlyReached())
529        potential.set(res_graph.head(e), potential[res_graph.tail(e)]+1);
530      ++bfs;
531    }
532
533    //Subgraph containing the edges on some shortest paths
534    //(i.e. tight edges)
535    ConstMap<typename ResGW::Node, bool> true_map(true);
536    typedef SubGraphWrapper<ResGW, ConstMap<typename ResGW::Node, bool>,
537      TightEdgeFilterMap<ResGW, typename ResGW::template NodeMap<int>,
538      Const1Map> > FilterResGW;
539    FilterResGW filter_res_graph(res_graph, true_map, tight_edge_filter);
540
541    //Subgraph, which is able to delete edges which are already
542    //met by the dfs
543    typename FilterResGW::template NodeMap<typename FilterResGW::Edge>
544      first_out_edges(filter_res_graph);
545    for (typename FilterResGW::NodeIt v(filter_res_graph); v!=INVALID; ++v)
546      first_out_edges.set
547        (v, typename FilterResGW::OutEdgeIt(filter_res_graph, v));
548
549    typedef ErasingFirstGraphWrapper<FilterResGW, typename FilterResGW::
550      template NodeMap<typename FilterResGW::Edge> > ErasingResGW;
551    ErasingResGW erasing_res_graph(filter_res_graph, first_out_edges);
552
553    bool __augment=true;
554
555    while (__augment) {
556
557      __augment=false;
558      //computing blocking flow with dfs
559      DfsIterator< ErasingResGW,
560        typename ErasingResGW::template NodeMap<bool> >
561        dfs(erasing_res_graph);
562      typename ErasingResGW::
563        template NodeMap<typename ErasingResGW::Edge> pred(erasing_res_graph);
564      pred.set(s, INVALID);
565      //invalid iterators for sources
566
567      typename ErasingResGW::template NodeMap<Num>
568        free1(erasing_res_graph);
569
570      dfs.pushAndSetReached
571        /// \bug hugo 0.2
572        (typename ErasingResGW::Node
573         (typename FilterResGW::Node
574          (typename ResGW::Node(s)
575           )
576          )
577         );
578       
579      while (!dfs.finished()) {
580        ++dfs;
581        if (typename ErasingResGW::Edge(dfs)!=INVALID) {
582          if (dfs.isBNodeNewlyReached()) {
583           
584            typename ErasingResGW::Node v=erasing_res_graph.tail(dfs);
585            typename ErasingResGW::Node w=erasing_res_graph.head(dfs);
586
587            pred.set(w, typename ErasingResGW::Edge(dfs));
588            if (pred[v]!=INVALID) {
589              free1.set
590                (w, std::min(free1[v], res_graph.resCap
591                             (typename ErasingResGW::Edge(dfs))));
592            } else {
593              free1.set
594                (w, res_graph.resCap
595                 (typename ErasingResGW::Edge(dfs)));
596            }
597
598            if (w==t) {
599              __augment=true;
600              _augment=true;
601              break;
602            }
603          } else {
604            erasing_res_graph.erase(dfs);
605          }
606        }
607      }
608
609      if (__augment) {
610        typename ErasingResGW::Node
611          n=typename FilterResGW::Node(typename ResGW::Node(t));
612        Num augment_value=free1[n];
613        while (pred[n]!=INVALID) {
614          typename ErasingResGW::Edge e=pred[n];
615          res_graph.augment(e, augment_value);
616          n=erasing_res_graph.tail(e);
617          if (res_graph.resCap(e)==0)
618            erasing_res_graph.erase(e);
619        }
620      }
621
622    } //while (__augment)
623
624    status=AFTER_AUGMENTING;
625    return _augment;
626  }
627
628
629} //namespace hugo
630
631#endif //HUGO_AUGMENTING_FLOW_H
632
633
Note: See TracBrowser for help on using the repository browser.