COIN-OR::LEMON - Graph Library

source: lemon-0.x/src/work/marci/augmenting_flow.h @ 987:87f7c54892df

Last change on this file since 987:87f7c54892df was 987:87f7c54892df, checked in by Alpar Juttner, 15 years ago

Naming changes:

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