COIN-OR::LEMON - Graph Library

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