COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/dinitz_sleator_tarjan.h @ 2517:d9cfac072869

Last change on this file since 2517:d9cfac072869 was 2514:57143c09dc20, checked in by Balazs Dezso, 12 years ago

Redesign the maximum flow algorithms

Redesigned interface
Preflow changed to use elevator
Edmonds-Karp does not use the ResGraphAdaptor?
Goldberg-Tarjan algorithm (Preflow with Dynamic Trees)
Dinitz-Sleator-Tarjan (Blocking flow with Dynamic Tree)

File size: 19.7 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2007
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_DINITZ_SLEATOR_TARJAN_H
20#define LEMON_DINITZ_SLEATOR_TARJAN_H
21
22/// \file
23/// \ingroup max_flow
24/// \brief Implementation the dynamic tree data structure of Sleator
25/// and Tarjan.
26
27#include <lemon/time_measure.h>
28#include <queue>
29#include <lemon/graph_utils.h>
30#include <lemon/tolerance.h>
31#include <lemon/graph_adaptor.h>
32#include <lemon/bfs.h>
33#include <lemon/edge_set.h>
34#include <lemon/dynamic_tree.h>
35
36#include <vector>
37#include <limits>
38#include <fstream>
39
40
41namespace lemon {
42
43  /// \brief Default traits class of DinitzSleatorTarjan class.
44  ///
45  /// Default traits class of DinitzSleatorTarjan class.
46  /// \param _Graph Graph type.
47  /// \param _CapacityMap Type of capacity map.
48  template <typename _Graph, typename _CapacityMap>
49  struct DinitzSleatorTarjanDefaultTraits {
50
51    /// \brief The graph type the algorithm runs on.
52    typedef _Graph Graph;
53
54    /// \brief The type of the map that stores the edge capacities.
55    ///
56    /// The type of the map that stores the edge capacities.
57    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
58    typedef _CapacityMap CapacityMap;
59
60    /// \brief The type of the length of the edges.
61    typedef typename CapacityMap::Value Value;
62
63    /// \brief The map type that stores the flow values.
64    ///
65    /// The map type that stores the flow values.
66    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
67    typedef typename Graph::template EdgeMap<Value> FlowMap;
68
69    /// \brief Instantiates a FlowMap.
70    ///
71    /// This function instantiates a \ref FlowMap.
72    /// \param graph The graph, to which we would like to define the flow map.
73    static FlowMap* createFlowMap(const Graph& graph) {
74      return new FlowMap(graph);
75    }
76
77    /// \brief The tolerance used by the algorithm
78    ///
79    /// The tolerance used by the algorithm to handle inexact computation.
80    typedef Tolerance<Value> Tolerance;
81
82  };
83
84  /// \ingroup max_flow
85  ///
86  /// \brief Dinitz-Sleator-Tarjan algorithms class.
87  ///
88  /// This class provides an implementation of the \e
89  /// Dinitz-Sleator-Tarjan \e algorithm producing a flow of maximum
90  /// value in a directed graph. The DinitzSleatorTarjan algorithm is
91  /// the fastest known max flow algorithms wich using blocking
92  /// flow. It is an improvement of the Dinitz algorithm by using the
93  /// \ref DynamicTree "dynamic tree" data structure of Sleator and
94  /// Tarjan.
95  ///
96  /// This blocking flow algorithms builds a layered graph according
97  /// to \ref Bfs "breadth-first search" distance from the target node
98  /// in the reversed residual graph. The layered graph contains each
99  /// residual edge which steps one level down. After that the
100  /// algorithm constructs a blocking flow on the layered graph and
101  /// augments the overall flow with it. The number of the levels in
102  /// the layered graph is strictly increasing in each augmenting
103  /// phase therefore the number of the augmentings is at most
104  /// \f$n-1\f$.  The length of each phase is at most
105  /// \f$O(m\log(n))\f$, that the overall time complexity is
106  /// \f$O(nm\log(n))\f$.
107  ///
108  /// \param _Graph The directed graph type the algorithm runs on.
109  /// \param _CapacityMap The capacity map type.
110  /// \param _Traits Traits class to set various data types used by
111  /// the algorithm.  The default traits class is \ref
112  /// DinitzSleatorTarjanDefaultTraits.  See \ref
113  /// DinitzSleatorTarjanDefaultTraits for the documentation of a
114  /// Dinitz-Sleator-Tarjan traits class.
115  ///
116  /// \author Tamas Hamori and Balazs Dezso
117#ifdef DOXYGEN
118  template <typename _Graph, typename _CapacityMap, typename _Traits>
119#else
120  template <typename _Graph,
121            typename _CapacityMap = typename _Graph::template EdgeMap<int>,
122            typename _Traits =
123            DinitzSleatorTarjanDefaultTraits<_Graph, _CapacityMap> >
124#endif
125  class DinitzSleatorTarjan {
126  public:
127
128    typedef _Traits Traits;
129    typedef typename Traits::Graph Graph;
130    typedef typename Traits::CapacityMap CapacityMap;
131    typedef typename Traits::Value Value;
132
133    typedef typename Traits::FlowMap FlowMap;
134    typedef typename Traits::Tolerance Tolerance;
135
136
137  private:
138
139    GRAPH_TYPEDEFS(typename Graph);
140
141
142    typedef typename Graph::template NodeMap<int> LevelMap;
143    typedef typename Graph::template NodeMap<int> IntNodeMap;
144    typedef typename Graph::template NodeMap<Edge> EdgeNodeMap;
145    typedef DynamicTree<Value, IntNodeMap, Tolerance, false> DynTree;
146
147  private:
148   
149    const Graph& _graph;
150    const CapacityMap* _capacity;
151
152    Node _source, _target;
153
154    FlowMap* _flow;
155    bool _local_flow;
156
157    IntNodeMap* _level;
158    EdgeNodeMap* _dt_edges;
159   
160    IntNodeMap* _dt_index;
161    DynTree* _dt;
162
163    Tolerance _tolerance;
164   
165    Value _flow_value;
166    Value _max_value;
167
168
169  public:
170
171    typedef DinitzSleatorTarjan Create;
172
173    ///\name Named template parameters
174
175    ///@{
176
177    template <typename _FlowMap>
178    struct DefFlowMapTraits : public Traits {
179      typedef _FlowMap FlowMap;
180      static FlowMap *createFlowMap(const Graph&) {
181        throw UninitializedParameter();
182      }
183    };
184
185    /// \brief \ref named-templ-param "Named parameter" for setting
186    /// FlowMap type
187    ///
188    /// \ref named-templ-param "Named parameter" for setting FlowMap
189    /// type
190    template <typename _FlowMap>
191    struct DefFlowMap
192      : public DinitzSleatorTarjan<Graph, CapacityMap,
193                              DefFlowMapTraits<_FlowMap> > {
194      typedef DinitzSleatorTarjan<Graph, CapacityMap,
195                             DefFlowMapTraits<_FlowMap> > Create;
196    };
197
198    template <typename _Elevator>
199    struct DefElevatorTraits : public Traits {
200      typedef _Elevator Elevator;
201      static Elevator *createElevator(const Graph&, int) {
202        throw UninitializedParameter();
203      }
204    };
205
206    /// @}
207
208    /// \brief \ref Exception for the case when the source equals the target.
209    ///
210    /// \ref Exception for the case when the source equals the target.
211    ///
212    class InvalidArgument : public lemon::LogicError {
213    public:
214      virtual const char* what() const throw() {
215        return "lemon::DinitzSleatorTarjan::InvalidArgument";
216      }
217    };
218
219    /// \brief The constructor of the class.
220    ///
221    /// The constructor of the class.
222    /// \param graph The directed graph the algorithm runs on.
223    /// \param capacity The capacity of the edges.
224    /// \param source The source node.
225    /// \param target The target node.
226    DinitzSleatorTarjan(const Graph& graph, const CapacityMap& capacity,
227                        Node source, Node target)
228      : _graph(graph), _capacity(&capacity),
229        _source(source), _target(target),
230        _flow(0), _local_flow(false),
231        _level(0), _dt_edges(0),
232        _dt_index(0), _dt(0),
233        _tolerance(), _flow_value(), _max_value()
234    {
235      if (_source == _target) {
236        throw InvalidArgument();
237      }
238    }
239
240    /// \brief Destrcutor.
241    ///
242    /// Destructor.
243    ~DinitzSleatorTarjan() {
244      destroyStructures();
245    }
246
247    /// \brief Sets the capacity map.
248    ///
249    /// Sets the capacity map.
250    /// \return \c (*this)
251    DinitzSleatorTarjan& capacityMap(const CapacityMap& map) {
252      _capacity = &map;
253      return *this;
254    }
255
256    /// \brief Sets the flow map.
257    ///
258    /// Sets the flow map.
259    /// \return \c (*this)
260    DinitzSleatorTarjan& flowMap(FlowMap& map) {
261      if (_local_flow) {
262        delete _flow;
263        _local_flow = false;
264      }
265      _flow = &map;
266      return *this;
267    }
268
269    /// \brief Returns the flow map.
270    ///
271    /// \return The flow map.
272    const FlowMap& flowMap() {
273      return *_flow;
274    }
275
276    /// \brief Sets the source node.
277    ///
278    /// Sets the source node.
279    /// \return \c (*this)
280    DinitzSleatorTarjan& source(const Node& node) {
281      _source = node;
282      return *this;
283    }
284
285    /// \brief Sets the target node.
286    ///
287    /// Sets the target node.
288    /// \return \c (*this)
289    DinitzSleatorTarjan& target(const Node& node) {
290      _target = node;
291      return *this;
292    }
293
294    /// \brief Sets the tolerance used by algorithm.
295    ///
296    /// Sets the tolerance used by algorithm.
297    DinitzSleatorTarjan& tolerance(const Tolerance& tolerance) const {
298      _tolerance = tolerance;
299      if (_dt) {
300        _dt.tolerance(_tolerance);
301      }
302      return *this;
303    }
304
305    /// \brief Returns the tolerance used by algorithm.
306    ///
307    /// Returns the tolerance used by algorithm.
308    const Tolerance& tolerance() const {
309      return tolerance;
310    }
311
312  private:
313       
314    void createStructures() {
315      if (!_flow) {
316        _flow = Traits::createFlowMap(_graph);
317        _local_flow = true;
318      }
319      if (!_level) {
320        _level = new LevelMap(_graph);
321      }
322      if (!_dt_index && !_dt) {
323        _dt_index = new IntNodeMap(_graph);
324        _dt = new DynTree(*_dt_index, _tolerance);
325      }
326      if (!_dt_edges) {
327        _dt_edges = new EdgeNodeMap(_graph);
328      }
329      _max_value = _dt->maxValue();
330    }
331
332    void destroyStructures() {
333      if (_local_flow) {
334        delete _flow;
335      }
336      if (_level) {
337        delete _level;
338      }
339      if (_dt) {
340        delete _dt;
341      }
342      if (_dt_index) {
343        delete _dt_index;
344      }
345      if (_dt_edges) {
346        delete _dt_edges;
347      }
348    }
349
350    bool createLayeredGraph() {
351
352      for (NodeIt n(_graph); n != INVALID; ++n) {
353        _level->set(n, -2);
354      }
355     
356      int level = 0;
357
358      std::vector<Node> queue;
359      queue.push_back(_target);
360      _level->set(_target, level);
361     
362      while ((*_level)[_source] == -2 && !queue.empty()) {
363        std::vector<Node> nqueue;
364        ++level;
365       
366        for (int i = 0; i < int(queue.size()); ++i) {
367          Node n = queue[i];
368         
369          for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
370            Node v = _graph.target(e);
371            if ((*_level)[v] != -2) continue;
372            Value rem = (*_flow)[e];
373            if (!_tolerance.positive(rem)) continue;
374            _level->set(v, level);
375            nqueue.push_back(v);
376          }
377
378          for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
379            Node v = _graph.source(e);
380            if ((*_level)[v] != -2) continue;
381            Value rem = (*_capacity)[e] - (*_flow)[e];
382            if (!_tolerance.positive(rem)) continue;
383            _level->set(v, level);
384            nqueue.push_back(v);
385          }
386
387        }
388        queue.swap(nqueue);
389      }
390     
391      return (*_level)[_source] != -2;
392    }
393
394    void initEdges() {
395      for (NodeIt n(_graph); n != INVALID; ++n) {
396        _graph.firstOut((*_dt_edges)[n], n);
397      }
398    }
399       
400   
401    void augmentPath() {
402      Value rem;
403      Node n = _dt->findCost(_source, rem);
404      _flow_value += rem;
405      _dt->addCost(_source, - rem);
406
407      _dt->cut(n);
408      _dt->addCost(n, _max_value);
409
410      Edge e = (*_dt_edges)[n];
411      if (_graph.source(e) == n) {
412        _flow->set(e, (*_capacity)[e]);
413       
414        _graph.nextOut(e);
415        if (e == INVALID) {
416          _graph.firstIn(e, n);
417        }
418      } else {
419        _flow->set(e, 0);
420        _graph.nextIn(e);
421      }
422      _dt_edges->set(n, e);
423
424    }
425
426    bool advance(Node n) {
427      Edge e = (*_dt_edges)[n];
428      if (e == INVALID) return false;
429
430      Node u;
431      Value rem;     
432      if (_graph.source(e) == n) {
433        u = _graph.target(e);
434        while ((*_level)[n] != (*_level)[u] + 1 ||
435               !_tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
436          _graph.nextOut(e);
437          if (e == INVALID) break;
438          u = _graph.target(e);
439        }
440        if (e != INVALID) {
441          rem = (*_capacity)[e] - (*_flow)[e];
442        } else {
443          _graph.firstIn(e, n);
444          if (e == INVALID) {
445            _dt_edges->set(n, INVALID);
446            return false;
447          }
448          u = _graph.source(e);
449          while ((*_level)[n] != (*_level)[u] + 1 ||
450                 !_tolerance.positive((*_flow)[e])) {
451            _graph.nextIn(e);
452            if (e == INVALID) {
453              _dt_edges->set(n, INVALID);
454              return false;
455            }
456            u = _graph.source(e);
457          }
458          rem = (*_flow)[e];
459        }
460      } else {
461        u = _graph.source(e);
462        while ((*_level)[n] != (*_level)[u] + 1 ||
463               !_tolerance.positive((*_flow)[e])) {
464          _graph.nextIn(e);
465          if (e == INVALID) {
466            _dt_edges->set(n, INVALID);
467            return false;
468          }
469          u = _graph.source(e);
470        }
471        rem = (*_flow)[e];
472      }
473
474      _dt->addCost(n, - std::numeric_limits<Value>::max());
475      _dt->addCost(n, rem);
476      _dt->link(n, u);
477      _dt_edges->set(n, e);
478      return true;
479    }
480
481    void retreat(Node n) {
482      _level->set(n, -1);
483     
484      for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
485        Node u = _graph.target(e);
486        if ((*_dt_edges)[u] == e && _dt->findRoot(u) == n) {
487          Value rem;
488          _dt->findCost(u, rem);
489          _flow->set(e, rem);
490          _dt->cut(u);
491          _dt->addCost(u, - rem);
492          _dt->addCost(u, _max_value);
493        }
494      }
495      for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
496        Node u = _graph.source(e);
497        if ((*_dt_edges)[u] == e && _dt->findRoot(u) == n) {
498          Value rem;
499          _dt->findCost(u, rem);
500          _flow->set(e, (*_capacity)[e] - rem);
501          _dt->cut(u);
502          _dt->addCost(u, - rem);
503          _dt->addCost(u, _max_value);
504        }
505      }
506    }
507
508    void extractTrees() {
509      for (NodeIt n(_graph); n != INVALID; ++n) {
510       
511        Node w = _dt->findRoot(n);
512     
513        while (w != n) {
514     
515          Value rem;     
516          Node u = _dt->findCost(n, rem);
517
518          _dt->cut(u);
519          _dt->addCost(u, - rem);
520          _dt->addCost(u, _max_value);
521         
522          Edge e = (*_dt_edges)[u];
523          _dt_edges->set(u, INVALID);
524         
525          if (u == _graph.source(e)) {
526            _flow->set(e, (*_capacity)[e] - rem);
527          } else {
528            _flow->set(e, rem);
529          }
530         
531          w = _dt->findRoot(n);
532        }     
533      }
534    }
535
536
537  public:
538   
539    /// \name Execution control The simplest way to execute the
540    /// algorithm is to use the \c run() member functions.
541    /// \n
542    /// If you need more control on initial solution or
543    /// execution then you have to call one \ref init() function and then
544    /// the start() or multiple times the \c augment() member function. 
545   
546    ///@{
547
548    /// \brief Initializes the algorithm
549    ///
550    /// It sets the flow to empty flow.
551    void init() {
552      createStructures();
553
554      _dt->clear();
555      for (NodeIt n(_graph); n != INVALID; ++n) {
556        _dt->makeTree(n);
557        _dt->addCost(n, _max_value);
558      }
559
560      for (EdgeIt it(_graph); it != INVALID; ++it) {
561        _flow->set(it, 0);
562      }
563      _flow_value = 0;
564    }
565   
566    /// \brief Initializes the algorithm
567    ///
568    /// Initializes the flow to the \c flowMap. The \c flowMap should
569    /// contain a feasible flow, ie. in each node excluding the source
570    /// and the target the incoming flow should be equal to the
571    /// outgoing flow.
572    template <typename FlowMap>
573    void flowInit(const FlowMap& flowMap) {
574      createStructures();
575
576      _dt->clear();
577      for (NodeIt n(_graph); n != INVALID; ++n) {
578        _dt->makeTree(n);
579        _dt->addCost(n, _max_value);
580      }
581
582      for (EdgeIt e(_graph); e != INVALID; ++e) {
583        _flow->set(e, flowMap[e]);
584      }
585      _flow_value = 0;
586      for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
587        _flow_value += (*_flow)[jt];
588      }
589      for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
590        _flow_value -= (*_flow)[jt];
591      }
592    }
593
594    /// \brief Initializes the algorithm
595    ///
596    /// Initializes the flow to the \c flowMap. The \c flowMap should
597    /// contain a feasible flow, ie. in each node excluding the source
598    /// and the target the incoming flow should be equal to the
599    /// outgoing flow. 
600    /// \return %False when the given flowMap does not contain
601    /// feasible flow.
602    template <typename FlowMap>
603    bool checkedFlowInit(const FlowMap& flowMap) {
604      createStructures();
605
606      _dt->clear();
607      for (NodeIt n(_graph); n != INVALID; ++n) {
608        _dt->makeTree(n);
609        _dt->addCost(n, _max_value);
610      }
611
612      for (EdgeIt e(_graph); e != INVALID; ++e) {
613        _flow->set(e, flowMap[e]);
614      }
615      for (NodeIt it(_graph); it != INVALID; ++it) {
616        if (it == _source || it == _target) continue;
617        Value outFlow = 0;
618        for (OutEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
619          outFlow += (*_flow)[jt];
620        }
621        Value inFlow = 0;
622        for (InEdgeIt jt(_graph, it); jt != INVALID; ++jt) {
623          inFlow += (*_flow)[jt];
624        }
625        if (_tolerance.different(outFlow, inFlow)) {
626          return false;
627        }
628      }
629      for (EdgeIt it(_graph); it != INVALID; ++it) {
630        if (_tolerance.less((*_flow)[it], 0)) return false;
631        if (_tolerance.less((*_capacity)[it], (*_flow)[it])) return false;
632      }
633      _flow_value = 0;
634      for (OutEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
635        _flow_value += (*_flow)[jt];
636      }
637      for (InEdgeIt jt(_graph, _source); jt != INVALID; ++jt) {
638        _flow_value -= (*_flow)[jt];
639      }
640      return true;
641    }
642
643    /// \brief Executes the algorithm
644    ///
645    /// It runs augmenting phases by adding blocking flow until the
646    /// optimal solution is reached.
647    void start() {
648      while (augment());
649    }
650
651    /// \brief Augments the flow with a blocking flow on a layered
652    /// graph.
653    ///
654    /// This function builds a layered graph and then find a blocking
655    /// flow on this graph. The number of the levels in the layered
656    /// graph is strictly increasing in each augmenting phase
657    /// therefore the number of the augmentings is at most \f$ n-1
658    /// \f$.  The length of each phase is at most \f$ O(m \log(n))
659    /// \f$, that the overall time complexity is \f$ O(nm \log(n)) \f$.
660    /// \return %False when there is not residual path between the
661    /// source and the target so the current flow is a feasible and
662    /// optimal solution.
663    bool augment() {
664      Node n;
665
666      if (createLayeredGraph()) {
667       
668        Timer bf_timer;
669        initEdges();
670
671        n = _dt->findRoot(_source);
672        while (true) {
673          Edge e;
674          if (n == _target) {
675            augmentPath();
676          } else if (!advance(n)) {
677            if (n != _source) {
678              retreat(n);
679            } else {
680              break;
681            }
682          }
683          n = _dt->findRoot(_source);
684        }     
685        extractTrees();
686
687        return true;
688      } else {
689        return false;
690      }
691    }
692   
693    /// \brief runs the algorithm.
694    ///
695    /// It is just a shorthand for:
696    ///
697    ///\code
698    /// ek.init();
699    /// ek.start();
700    ///\endcode
701    void run() {
702      init();
703      start();
704    }
705
706    /// @}
707
708    /// \name Query Functions
709    /// The result of the %Dijkstra algorithm can be obtained using these
710    /// functions.\n
711    /// Before the use of these functions,
712    /// either run() or start() must be called.
713   
714    ///@{
715
716    /// \brief Returns the value of the maximum flow.
717    ///
718    /// Returns the value of the maximum flow by returning the excess
719    /// of the target node \c t. This value equals to the value of
720    /// the maximum flow already after the first phase.
721    Value flowValue() const {
722      return _flow_value;
723    }
724
725
726    /// \brief Returns the flow on the edge.
727    ///
728    /// Sets the \c flowMap to the flow on the edges. This method can
729    /// be called after the second phase of algorithm.
730    Value flow(const Edge& edge) const {
731      return (*_flow)[edge];
732    }
733
734    /// \brief Returns true when the node is on the source side of minimum cut.
735    ///
736
737    /// Returns true when the node is on the source side of minimum
738    /// cut. This method can be called both after running \ref
739    /// startFirstPhase() and \ref startSecondPhase().
740    bool minCut(const Node& node) const {
741      return (*_level)[node] == -2;
742    }
743
744    /// \brief Returns a minimum value cut.
745    ///
746    /// Sets \c cut to the characteristic vector of a minimum value cut
747    /// It simply calls the minMinCut member.
748    /// \retval cut Write node bool map.
749    template <typename CutMap>
750    void minCutMap(CutMap& cutMap) const {
751      for (NodeIt n(_graph); n != INVALID; ++n) {
752        cutMap.set(n, (*_level)[n] == -2);
753      }
754      cutMap.set(_source, true);
755    }   
756
757    /// @}
758
759  };
760}
761
762#endif
Note: See TracBrowser for help on using the repository browser.