COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/goldberg_tarjan.h @ 2553:bfced05fa852

Last change on this file since 2553:bfced05fa852 was 2553:bfced05fa852, checked in by Alpar Juttner, 11 years ago

Happy New Year to LEMON (+ better update-copyright-header script)

File size: 27.2 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2008
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 *
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
12 *
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
15 * purpose.
16 *
17 */
18
19#ifndef LEMON_GOLDBERG_TARJAN_H
20#define LEMON_GOLDBERG_TARJAN_H
21
22#include <vector>
23#include <queue>
24
25#include <lemon/error.h>
26#include <lemon/bits/invalid.h>
27#include <lemon/tolerance.h>
28#include <lemon/maps.h>
29#include <lemon/graph_utils.h>
30#include <lemon/dynamic_tree.h>
31#include <limits>
32
33/// \file
34/// \ingroup max_flow
35/// \brief Implementation of the preflow algorithm.
36
37namespace lemon {
38
39  /// \brief Default traits class of GoldbergTarjan class.
40  ///
41  /// Default traits class of GoldbergTarjan class.
42  /// \param _Graph Graph type.
43  /// \param _CapacityMap Type of capacity map.
44  template <typename _Graph, typename _CapacityMap>
45  struct GoldbergTarjanDefaultTraits {
46
47    /// \brief The graph type the algorithm runs on.
48    typedef _Graph Graph;
49
50    /// \brief The type of the map that stores the edge capacities.
51    ///
52    /// The type of the map that stores the edge capacities.
53    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
54    typedef _CapacityMap CapacityMap;
55
56    /// \brief The type of the length of the edges.
57    typedef typename CapacityMap::Value Value;
58
59    /// \brief The map type that stores the flow values.
60    ///
61    /// The map type that stores the flow values.
62    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
63    typedef typename Graph::template EdgeMap<Value> FlowMap;
64
65    /// \brief Instantiates a FlowMap.
66    ///
67    /// This function instantiates a \ref FlowMap.
68    /// \param graph The graph, to which we would like to define the flow map.
69    static FlowMap* createFlowMap(const Graph& graph) {
70      return new FlowMap(graph);
71    }
72
73    /// \brief The eleavator type used by GoldbergTarjan algorithm.
74    ///
75    /// The elevator type used by GoldbergTarjan algorithm.
76    ///
77    /// \sa Elevator
78    /// \sa LinkedElevator
79    typedef LinkedElevator<Graph, typename Graph::Node> Elevator;
80   
81    /// \brief Instantiates an Elevator.
82    ///
83    /// This function instantiates a \ref Elevator.
84    /// \param graph The graph, to which we would like to define the elevator.
85    /// \param max_level The maximum level of the elevator.
86    static Elevator* createElevator(const Graph& graph, int max_level) {
87      return new Elevator(graph, max_level);
88    }
89
90    /// \brief The tolerance used by the algorithm
91    ///
92    /// The tolerance used by the algorithm to handle inexact computation.
93    typedef Tolerance<Value> Tolerance;
94
95  };
96
97  /// \ingroup max_flow
98  /// \brief Goldberg-Tarjan algorithms class.
99  ///
100  /// This class provides an implementation of the \e GoldbergTarjan
101  /// \e algorithm producing a flow of maximum value in a directed
102  /// graph. The GoldbergTarjan algorithm is a theoretical improvement
103  /// of the Goldberg's \ref Preflow "preflow" algorithm by using the \ref
104  /// DynamicTree "dynamic tree" data structure of Sleator and Tarjan.
105  ///
106  /// The original preflow algorithm with \e highest \e label
107  /// heuristic has \f$O(n^2\sqrt{e})\f$ or with \e FIFO heuristic has
108  /// \f$O(n^3)\f$ time complexity. The current algorithm improved
109  /// this complexity to \f$O(nm\log(\frac{n^2}{m}))\f$. The algorithm
110  /// builds limited size dynamic trees on the residual graph, which
111  /// can be used to make multilevel push operations and gives a
112  /// better bound on the number of non-saturating pushes.
113  ///
114  /// \param Graph The directed graph type the algorithm runs on.
115  /// \param CapacityMap The capacity map type.
116  /// \param _Traits Traits class to set various data types used by
117  /// the algorithm.  The default traits class is \ref
118  /// GoldbergTarjanDefaultTraits.  See \ref
119  /// GoldbergTarjanDefaultTraits for the documentation of a
120  /// Goldberg-Tarjan traits class.
121  ///
122  /// \author Tamas Hamori and Balazs Dezso
123#ifdef DOXYGEN
124  template <typename _Graph, typename _CapacityMap, typename _Traits>
125#else
126  template <typename _Graph,
127            typename _CapacityMap = typename _Graph::template EdgeMap<int>,
128            typename _Traits =
129            GoldbergTarjanDefaultTraits<_Graph, _CapacityMap> >
130#endif
131  class GoldbergTarjan {
132  public:
133
134    typedef _Traits Traits;
135    typedef typename Traits::Graph Graph;
136    typedef typename Traits::CapacityMap CapacityMap;
137    typedef typename Traits::Value Value;
138
139    typedef typename Traits::FlowMap FlowMap;
140    typedef typename Traits::Elevator Elevator;
141    typedef typename Traits::Tolerance Tolerance;
142   
143  protected:
144
145    GRAPH_TYPEDEFS(typename Graph);
146
147    typedef typename Graph::template NodeMap<Node> NodeNodeMap;
148    typedef typename Graph::template NodeMap<int> IntNodeMap;
149
150    typedef typename Graph::template NodeMap<Edge> EdgeNodeMap;
151    typedef typename Graph::template EdgeMap<Edge> EdgeEdgeMap;
152
153    typedef typename std::vector<Node> VecNode;
154 
155    typedef DynamicTree<Value,IntNodeMap,Tolerance> DynTree;
156
157    const Graph& _graph;
158    const CapacityMap* _capacity;
159    int _node_num; //the number of nodes of G
160
161    Node _source;
162    Node _target;
163
164    FlowMap* _flow;
165    bool _local_flow;
166
167    Elevator* _level;
168    bool _local_level;
169
170    typedef typename Graph::template NodeMap<Value> ExcessMap;
171    ExcessMap* _excess;
172   
173    Tolerance _tolerance;
174
175    bool _phase;
176
177    // constant for treesize
178    static const double _tree_bound = 2;
179    int _max_tree_size;
180   
181    //tags for the dynamic tree
182    DynTree *_dt;
183    //datastructure of dyanamic tree
184    IntNodeMap *_dt_index;
185    //datastrucure for solution of the communication between the two class
186    EdgeNodeMap *_dt_edges;
187    //nodeMap for storing the outgoing edge from the node in the tree 
188
189    //max of the type Value
190    const Value _max_value;
191           
192  public:
193
194    typedef GoldbergTarjan Create;
195
196    ///\name Named template parameters
197
198    ///@{
199
200    template <typename _FlowMap>
201    struct DefFlowMapTraits : public Traits {
202      typedef _FlowMap FlowMap;
203      static FlowMap *createFlowMap(const Graph&) {
204        throw UninitializedParameter();
205      }
206    };
207
208    /// \brief \ref named-templ-param "Named parameter" for setting
209    /// FlowMap type
210    ///
211    /// \ref named-templ-param "Named parameter" for setting FlowMap
212    /// type
213    template <typename _FlowMap>
214    struct DefFlowMap
215      : public GoldbergTarjan<Graph, CapacityMap,
216                              DefFlowMapTraits<_FlowMap> > {
217      typedef GoldbergTarjan<Graph, CapacityMap,
218                             DefFlowMapTraits<_FlowMap> > Create;
219    };
220
221    template <typename _Elevator>
222    struct DefElevatorTraits : public Traits {
223      typedef _Elevator Elevator;
224      static Elevator *createElevator(const Graph&, int) {
225        throw UninitializedParameter();
226      }
227    };
228
229    /// \brief \ref named-templ-param "Named parameter" for setting
230    /// Elevator type
231    ///
232    /// \ref named-templ-param "Named parameter" for setting Elevator
233    /// type
234    template <typename _Elevator>
235    struct DefElevator
236      : public GoldbergTarjan<Graph, CapacityMap,
237                              DefElevatorTraits<_Elevator> > {
238      typedef GoldbergTarjan<Graph, CapacityMap,
239                             DefElevatorTraits<_Elevator> > Create;
240    };
241
242    template <typename _Elevator>
243    struct DefStandardElevatorTraits : public Traits {
244      typedef _Elevator Elevator;
245      static Elevator *createElevator(const Graph& graph, int max_level) {
246        return new Elevator(graph, max_level);
247      }
248    };
249
250    /// \brief \ref named-templ-param "Named parameter" for setting
251    /// Elevator type
252    ///
253    /// \ref named-templ-param "Named parameter" for setting Elevator
254    /// type. The Elevator should be standard constructor interface, ie.
255    /// the graph and the maximum level should be passed to it.
256    template <typename _Elevator>
257    struct DefStandardElevator
258      : public GoldbergTarjan<Graph, CapacityMap,
259                              DefStandardElevatorTraits<_Elevator> > {
260      typedef GoldbergTarjan<Graph, CapacityMap,
261                             DefStandardElevatorTraits<_Elevator> > Create;
262    };   
263
264
265    ///\ref Exception for the case when s=t.
266
267    ///\ref Exception for the case when the source equals the target.
268    class InvalidArgument : public lemon::LogicError {
269    public:
270      virtual const char* what() const throw() {
271        return "lemon::GoldbergTarjan::InvalidArgument";
272      }
273    };
274
275  protected:
276   
277    GoldbergTarjan() {}
278
279  public:
280
281    /// \brief The constructor of the class.
282    ///
283    /// The constructor of the class.
284    /// \param graph The directed graph the algorithm runs on.
285    /// \param capacity The capacity of the edges.
286    /// \param source The source node.
287    /// \param target The target node.
288    /// Except the graph, all of these parameters can be reset by
289    /// calling \ref source, \ref target, \ref capacityMap and \ref
290    /// flowMap, resp.
291    GoldbergTarjan(const Graph& graph, const CapacityMap& capacity,
292                   Node source, Node target)
293      : _graph(graph), _capacity(&capacity),
294        _node_num(), _source(source), _target(target),
295        _flow(0), _local_flow(false),
296        _level(0), _local_level(false),
297        _excess(0), _tolerance(),
298        _phase(), _max_tree_size(),
299        _dt(0), _dt_index(0), _dt_edges(0),
300        _max_value(std::numeric_limits<Value>::max()) {
301      if (_source == _target) throw InvalidArgument();
302    }
303
304    /// \brief Destrcutor.
305    ///
306    /// Destructor.
307    ~GoldbergTarjan() {
308      destroyStructures();
309    }
310   
311    /// \brief Sets the capacity map.
312    ///
313    /// Sets the capacity map.
314    /// \return \c (*this)
315    GoldbergTarjan& capacityMap(const CapacityMap& map) {
316      _capacity = &map;
317      return *this;
318    }
319
320    /// \brief Sets the flow map.
321    ///
322    /// Sets the flow map.
323    /// \return \c (*this)
324    GoldbergTarjan& flowMap(FlowMap& map) {
325      if (_local_flow) {
326        delete _flow;
327        _local_flow = false;
328      }
329      _flow = &map;
330      return *this;
331    }
332
333    /// \brief Returns the flow map.
334    ///
335    /// \return The flow map.
336    const FlowMap& flowMap() {
337      return *_flow;
338    }
339
340    /// \brief Sets the elevator.
341    ///
342    /// Sets the elevator.
343    /// \return \c (*this)
344    GoldbergTarjan& elevator(Elevator& elevator) {
345      if (_local_level) {
346        delete _level;
347        _local_level = false;
348      }
349      _level = &elevator;
350      return *this;
351    }
352
353    /// \brief Returns the elevator.
354    ///
355    /// \return The elevator.
356    const Elevator& elevator() {
357      return *_level;
358    }
359
360    /// \brief Sets the source node.
361    ///
362    /// Sets the source node.
363    /// \return \c (*this)
364    GoldbergTarjan& source(const Node& node) {
365      _source = node;
366      return *this;
367    }
368
369    /// \brief Sets the target node.
370    ///
371    /// Sets the target node.
372    /// \return \c (*this)
373    GoldbergTarjan& target(const Node& node) {
374      _target = node;
375      return *this;
376    }
377 
378    /// \brief Sets the tolerance used by algorithm.
379    ///
380    /// Sets the tolerance used by algorithm.
381    GoldbergTarjan& tolerance(const Tolerance& tolerance) const {
382      _tolerance = tolerance;
383      if (_dt) {
384        _dt->tolerance(_tolerance);
385      }
386      return *this;
387    }
388
389    /// \brief Returns the tolerance used by algorithm.
390    ///
391    /// Returns the tolerance used by algorithm.
392    const Tolerance& tolerance() const {
393      return tolerance;
394    }
395
396         
397  private:
398   
399    void createStructures() {
400      _node_num = countNodes(_graph);
401
402      _max_tree_size = int((double(_node_num) * double(_node_num)) /
403                           double(countEdges(_graph)));
404
405      if (!_flow) {
406        _flow = Traits::createFlowMap(_graph);
407        _local_flow = true;
408      }
409      if (!_level) {
410        _level = Traits::createElevator(_graph, _node_num);
411        _local_level = true;
412      }
413      if (!_excess) {
414        _excess = new ExcessMap(_graph);
415      }
416      if (!_dt_index && !_dt) {
417        _dt_index = new IntNodeMap(_graph);
418        _dt = new DynTree(*_dt_index, _tolerance);
419      }
420      if (!_dt_edges) {
421        _dt_edges = new EdgeNodeMap(_graph);
422      }
423    }
424
425    void destroyStructures() {
426      if (_local_flow) {
427        delete _flow;
428      }
429      if (_local_level) {
430        delete _level;
431      }
432      if (_excess) {
433        delete _excess;
434      }
435      if (_dt) {
436        delete _dt;
437      }
438      if (_dt_index) {
439        delete _dt_index;
440      }
441      if (_dt_edges) {
442        delete _dt_edges;
443      }
444    }
445
446    bool sendOut(Node n, Value& excess) {
447
448      Node w = _dt->findRoot(n);
449     
450      while (w != n) {
451     
452        Value rem;     
453        Node u = _dt->findCost(n, rem);
454
455        if (_tolerance.positive(rem) && !_level->active(w) && w != _target) {
456          _level->activate(w);
457        }
458
459        if (!_tolerance.less(rem, excess)) {
460          _dt->addCost(n, - excess);
461          _excess->set(w, (*_excess)[w] + excess);
462          excess = 0;
463          return true;
464        }
465       
466       
467        _dt->addCost(n, - rem);
468
469        excess -= rem;
470        _excess->set(w, (*_excess)[w] + rem);
471       
472        _dt->cut(u);
473        _dt->addCost(u, _max_value);
474         
475        Edge e = (*_dt_edges)[u];
476        _dt_edges->set(u, INVALID);
477       
478        if (u == _graph.source(e)) {
479          _flow->set(e, (*_capacity)[e]);
480        } else {
481          _flow->set(e, 0);
482        }           
483
484        w = _dt->findRoot(n);
485      }     
486      return false;
487    }
488
489    bool sendIn(Node n, Value& excess) {
490
491      Node w = _dt->findRoot(n);
492     
493      while (w != n) {
494     
495        Value rem;     
496        Node u = _dt->findCost(n, rem);
497
498        if (_tolerance.positive(rem) && !_level->active(w) && w != _source) {
499          _level->activate(w);
500        }
501
502        if (!_tolerance.less(rem, excess)) {
503          _dt->addCost(n, - excess);
504          _excess->set(w, (*_excess)[w] + excess);
505          excess = 0;
506          return true;
507        }
508       
509       
510        _dt->addCost(n, - rem);
511
512        excess -= rem;
513        _excess->set(w, (*_excess)[w] + rem);
514       
515        _dt->cut(u);
516        _dt->addCost(u, _max_value);
517         
518        Edge e = (*_dt_edges)[u];
519        _dt_edges->set(u, INVALID);
520       
521        if (u == _graph.source(e)) {
522          _flow->set(e, (*_capacity)[e]);
523        } else {
524          _flow->set(e, 0);
525        }           
526
527        w = _dt->findRoot(n);
528      }     
529      return false;
530    }
531
532    void cutChildren(Node n) {
533   
534      for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
535       
536        Node v = _graph.target(e);
537       
538        if ((*_dt_edges)[v] != INVALID && (*_dt_edges)[v] == e) {
539          _dt->cut(v);
540          _dt_edges->set(v, INVALID);
541          Value rem;
542          _dt->findCost(v, rem);
543          _dt->addCost(v, - rem);
544          _dt->addCost(v, _max_value);
545          _flow->set(e, rem);
546
547        }     
548      }
549
550      for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
551       
552        Node v = _graph.source(e);
553       
554        if ((*_dt_edges)[v] != INVALID && (*_dt_edges)[v] == e) {
555          _dt->cut(v);
556          _dt_edges->set(v, INVALID);
557          Value rem;
558          _dt->findCost(v, rem);
559          _dt->addCost(v, - rem);
560          _dt->addCost(v, _max_value);
561          _flow->set(e, (*_capacity)[e] - rem);
562
563        }     
564      }
565    }
566
567    void extractTrees() {
568      for (NodeIt n(_graph); n != INVALID; ++n) {
569       
570        Node w = _dt->findRoot(n);
571     
572        while (w != n) {
573     
574          Value rem;     
575          Node u = _dt->findCost(n, rem);
576
577          _dt->cut(u);
578          _dt->addCost(u, - rem);
579          _dt->addCost(u, _max_value);
580         
581          Edge e = (*_dt_edges)[u];
582          _dt_edges->set(u, INVALID);
583         
584          if (u == _graph.source(e)) {
585            _flow->set(e, (*_capacity)[e] - rem);
586          } else {
587            _flow->set(e, rem);
588          }
589         
590          w = _dt->findRoot(n);
591        }     
592      }
593    }
594
595  public:   
596
597    /// \name Execution control The simplest way to execute the
598    /// algorithm is to use one of the member functions called \c
599    /// run(). 
600    /// \n
601    /// If you need more control on initial solution or
602    /// execution then you have to call one \ref init() function and then
603    /// the startFirstPhase() and if you need the startSecondPhase(). 
604   
605    ///@{
606
607    /// \brief Initializes the internal data structures.
608    ///
609    /// Initializes the internal data structures.
610    ///
611    void init() {
612      createStructures();
613
614      for (NodeIt n(_graph); n != INVALID; ++n) {
615        _excess->set(n, 0);
616      }
617
618      for (EdgeIt e(_graph); e != INVALID; ++e) {
619        _flow->set(e, 0);
620      }
621
622      _dt->clear();
623      for (NodeIt v(_graph); v != INVALID; ++v) {
624        (*_dt_edges)[v] = INVALID;
625        _dt->makeTree(v);
626        _dt->addCost(v, _max_value);
627      }
628
629      typename Graph::template NodeMap<bool> reached(_graph, false);
630
631      _level->initStart();
632      _level->initAddItem(_target);
633
634      std::vector<Node> queue;
635      reached.set(_source, true);
636
637      queue.push_back(_target);
638      reached.set(_target, true);
639      while (!queue.empty()) {
640        _level->initNewLevel();
641        std::vector<Node> nqueue;
642        for (int i = 0; i < int(queue.size()); ++i) {
643          Node n = queue[i];
644          for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
645            Node u = _graph.source(e);
646            if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
647              reached.set(u, true);
648              _level->initAddItem(u);
649              nqueue.push_back(u);
650            }
651          }
652        }
653        queue.swap(nqueue);
654      }
655      _level->initFinish();
656
657      for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) {
658        if (_tolerance.positive((*_capacity)[e])) {
659          Node u = _graph.target(e);
660          if ((*_level)[u] == _level->maxLevel()) continue;
661          _flow->set(e, (*_capacity)[e]);
662          _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
663          if (u != _target && !_level->active(u)) {
664            _level->activate(u);
665          }
666        }
667      }
668    }
669
670    /// \brief Starts the first phase of the preflow algorithm.
671    ///
672    /// The preflow algorithm consists of two phases, this method runs
673    /// the first phase. After the first phase the maximum flow value
674    /// and a minimum value cut can already be computed, although a
675    /// maximum flow is not yet obtained. So after calling this method
676    /// \ref flowValue() returns the value of a maximum flow and \ref
677    /// minCut() returns a minimum cut.     
678    /// \pre One of the \ref init() functions should be called.
679    void startFirstPhase() {
680      _phase = true;
681      Node n;
682
683      while ((n = _level->highestActive()) != INVALID) {
684        Value excess = (*_excess)[n];
685        int level = _level->highestActiveLevel();
686        int new_level = _level->maxLevel();
687
688        if (_dt->findRoot(n) != n) {
689          if (sendOut(n, excess)) goto no_more_push;
690        }
691       
692        for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
693          Value rem = (*_capacity)[e] - (*_flow)[e];
694          Node v = _graph.target(e);
695         
696          if (!_tolerance.positive(rem) && (*_dt_edges)[v] != e) continue;
697
698          if ((*_level)[v] < level) {
699           
700            if (_dt->findSize(n) + _dt->findSize(v) <
701                _tree_bound * _max_tree_size) {
702              _dt->addCost(n, -_max_value);
703              _dt->addCost(n, rem);
704              _dt->link(n, v);
705              _dt_edges->set(n, e);
706              if (sendOut(n, excess)) goto no_more_push;
707            } else {
708              if (!_level->active(v) && v != _target) {
709                _level->activate(v);
710              }
711              if (!_tolerance.less(rem, excess)) {
712                _flow->set(e, (*_flow)[e] + excess);
713                _excess->set(v, (*_excess)[v] + excess);
714                excess = 0;               
715                goto no_more_push;
716              } else {
717                excess -= rem;
718                _excess->set(v, (*_excess)[v] + rem);
719                _flow->set(e, (*_capacity)[e]);
720              }         
721            }
722          } else if (new_level > (*_level)[v]) {
723            new_level = (*_level)[v];
724          }
725        }
726
727        for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
728          Value rem = (*_flow)[e];
729          Node v = _graph.source(e);
730          if (!_tolerance.positive(rem) && (*_dt_edges)[v] != e) continue;
731
732          if ((*_level)[v] < level) {
733           
734            if (_dt->findSize(n) + _dt->findSize(v) <
735                _tree_bound * _max_tree_size) {
736              _dt->addCost(n, - _max_value);
737              _dt->addCost(n, rem);
738              _dt->link(n, v);
739              _dt_edges->set(n, e);
740              if (sendOut(n, excess)) goto no_more_push;
741            } else {
742              if (!_level->active(v) && v != _target) {
743                _level->activate(v);
744              }
745              if (!_tolerance.less(rem, excess)) {
746                _flow->set(e, (*_flow)[e] - excess);
747                _excess->set(v, (*_excess)[v] + excess);
748                excess = 0;               
749                goto no_more_push;
750              } else {
751                excess -= rem;
752                _excess->set(v, (*_excess)[v] + rem);
753                _flow->set(e, 0);
754              }         
755            }
756          } else if (new_level > (*_level)[v]) {
757            new_level = (*_level)[v];
758          }
759        }
760               
761      no_more_push:
762
763        _excess->set(n, excess);
764       
765        if (excess != 0) {
766          cutChildren(n);
767          if (new_level + 1 < _level->maxLevel()) {
768            _level->liftHighestActive(new_level + 1);
769          } else {
770            _level->liftHighestActiveToTop();
771          }
772          if (_level->emptyLevel(level)) {
773            _level->liftToTop(level);
774          }
775        } else {
776          _level->deactivate(n);
777        }       
778      }
779      extractTrees();
780    }
781
782    /// \brief Starts the second phase of the preflow algorithm.
783    ///
784    /// The preflow algorithm consists of two phases, this method runs
785    /// the second phase. After calling \ref init() and \ref
786    /// startFirstPhase() and then \ref startSecondPhase(), \ref
787    /// flowMap() return a maximum flow, \ref flowValue() returns the
788    /// value of a maximum flow, \ref minCut() returns a minimum cut
789    /// \pre The \ref init() and startFirstPhase() functions should be
790    /// called before.
791    void startSecondPhase() {
792      _phase = false;
793     
794      typename Graph::template NodeMap<bool> reached(_graph, false);
795      for (NodeIt n(_graph); n != INVALID; ++n) {
796        reached.set(n, (*_level)[n] < _level->maxLevel());
797      }
798
799      _level->initStart();
800      _level->initAddItem(_source);
801 
802      std::vector<Node> queue;
803      queue.push_back(_source);
804      reached.set(_source, true);
805
806      while (!queue.empty()) {
807        _level->initNewLevel();
808        std::vector<Node> nqueue;
809        for (int i = 0; i < int(queue.size()); ++i) {
810          Node n = queue[i];
811          for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
812            Node v = _graph.target(e);
813            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
814              reached.set(v, true);
815              _level->initAddItem(v);
816              nqueue.push_back(v);
817            }
818          }
819          for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
820            Node u = _graph.source(e);
821            if (!reached[u] &&
822                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
823              reached.set(u, true);
824              _level->initAddItem(u);
825              nqueue.push_back(u);
826            }
827          }
828        }
829        queue.swap(nqueue);
830      }
831      _level->initFinish();
832
833      for (NodeIt n(_graph); n != INVALID; ++n) {
834        if (!reached[n]) {
835          _level->markToBottom(n);
836        } else if ((*_excess)[n] > 0 && _target != n) {
837          _level->activate(n);
838        }
839      }
840
841      Node n;
842
843      while ((n = _level->highestActive()) != INVALID) {
844        Value excess = (*_excess)[n];
845        int level = _level->highestActiveLevel();
846        int new_level = _level->maxLevel();
847
848        if (_dt->findRoot(n) != n) {
849          if (sendIn(n, excess)) goto no_more_push;
850        }
851       
852        for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
853          Value rem = (*_capacity)[e] - (*_flow)[e];
854          Node v = _graph.target(e);
855         
856          if (!_tolerance.positive(rem) && (*_dt_edges)[v] != e) continue;
857
858          if ((*_level)[v] < level) {
859           
860            if (_dt->findSize(n) + _dt->findSize(v) <
861                _tree_bound * _max_tree_size) {
862              _dt->addCost(n, -_max_value);
863              _dt->addCost(n, rem);
864              _dt->link(n, v);
865              _dt_edges->set(n, e);
866              if (sendIn(n, excess)) goto no_more_push;
867            } else {
868              if (!_level->active(v) && v != _source) {
869                _level->activate(v);
870              }
871              if (!_tolerance.less(rem, excess)) {
872                _flow->set(e, (*_flow)[e] + excess);
873                _excess->set(v, (*_excess)[v] + excess);
874                excess = 0;               
875                goto no_more_push;
876              } else {
877                excess -= rem;
878                _excess->set(v, (*_excess)[v] + rem);
879                _flow->set(e, (*_capacity)[e]);
880              }         
881            }
882          } else if (new_level > (*_level)[v]) {
883            new_level = (*_level)[v];
884          }
885        }
886
887        for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
888          Value rem = (*_flow)[e];
889          Node v = _graph.source(e);
890          if (!_tolerance.positive(rem) && (*_dt_edges)[v] != e) continue;
891
892          if ((*_level)[v] < level) {
893           
894            if (_dt->findSize(n) + _dt->findSize(v) <
895                _tree_bound * _max_tree_size) {
896              _dt->addCost(n, - _max_value);
897              _dt->addCost(n, rem);
898              _dt->link(n, v);
899              _dt_edges->set(n, e);
900              if (sendIn(n, excess)) goto no_more_push;
901            } else {
902              if (!_level->active(v) && v != _source) {
903                _level->activate(v);
904              }
905              if (!_tolerance.less(rem, excess)) {
906                _flow->set(e, (*_flow)[e] - excess);
907                _excess->set(v, (*_excess)[v] + excess);
908                excess = 0;               
909                goto no_more_push;
910              } else {
911                excess -= rem;
912                _excess->set(v, (*_excess)[v] + rem);
913                _flow->set(e, 0);
914              }         
915            }
916          } else if (new_level > (*_level)[v]) {
917            new_level = (*_level)[v];
918          }
919        }
920               
921      no_more_push:
922
923        _excess->set(n, excess);
924       
925        if (excess != 0) {
926          cutChildren(n);
927          if (new_level + 1 < _level->maxLevel()) {
928            _level->liftHighestActive(new_level + 1);
929          } else {
930            _level->liftHighestActiveToTop();
931          }
932          if (_level->emptyLevel(level)) {
933            _level->liftToTop(level);
934          }
935        } else {
936          _level->deactivate(n);
937        }       
938      }
939      extractTrees();
940    }
941
942    /// \brief Runs the Goldberg-Tarjan algorithm. 
943    ///
944    /// Runs the Goldberg-Tarjan algorithm.
945    /// \note pf.run() is just a shortcut of the following code.
946    /// \code
947    ///   pf.init();
948    ///   pf.startFirstPhase();
949    ///   pf.startSecondPhase();
950    /// \endcode
951    void run() {
952      init();
953      startFirstPhase();
954      startSecondPhase();
955    }
956
957    /// \brief Runs the Goldberg-Tarjan algorithm to compute the minimum cut. 
958    ///
959    /// Runs the Goldberg-Tarjan algorithm to compute the minimum cut.
960    /// \note pf.runMinCut() is just a shortcut of the following code.
961    /// \code
962    ///   pf.init();
963    ///   pf.startFirstPhase();
964    /// \endcode
965    void runMinCut() {
966      init();
967      startFirstPhase();
968    }
969
970    /// @}
971
972    /// \name Query Functions
973    /// The result of the Goldberg-Tarjan algorithm can be obtained
974    /// using these functions.
975    /// \n
976    /// Before the use of these functions, either run() or start() must
977    /// be called.
978   
979    ///@{
980
981    /// \brief Returns the value of the maximum flow.
982    ///
983    /// Returns the value of the maximum flow by returning the excess
984    /// of the target node \c t. This value equals to the value of
985    /// the maximum flow already after the first phase.
986    Value flowValue() const {
987      return (*_excess)[_target];
988    }
989
990    /// \brief Returns true when the node is on the source side of minimum cut.
991    ///
992    /// Returns true when the node is on the source side of minimum
993    /// cut. This method can be called both after running \ref
994    /// startFirstPhase() and \ref startSecondPhase().
995    bool minCut(const Node& node) const {
996      return ((*_level)[node] == _level->maxLevel()) == _phase;
997    }
998 
999    /// \brief Returns a minimum value cut.
1000    ///
1001    /// Sets the \c cutMap to the characteristic vector of a minimum value
1002    /// cut. This method can be called both after running \ref
1003    /// startFirstPhase() and \ref startSecondPhase(). The result after second
1004    /// phase could be changed slightly if inexact computation is used.
1005    /// \pre The \c cutMap should be a bool-valued node-map.
1006    template <typename CutMap>
1007    void minCutMap(CutMap& cutMap) const {
1008      for (NodeIt n(_graph); n != INVALID; ++n) {
1009        cutMap.set(n, minCut(n));
1010      }
1011    }
1012
1013    /// \brief Returns the flow on the edge.
1014    ///
1015    /// Sets the \c flowMap to the flow on the edges. This method can
1016    /// be called after the second phase of algorithm.
1017    Value flow(const Edge& edge) const {
1018      return (*_flow)[edge];
1019    }
1020
1021    /// @}
1022
1023  };
1024 
1025} //namespace lemon
1026
1027#endif
Note: See TracBrowser for help on using the repository browser.