COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/marci/augmenting_flow.h @ 888:cc3590763f7f

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