COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/hao_orlin.h @ 2282:9d7b12f83daa

Last change on this file since 2282:9d7b12f83daa was 2275:ff46747676ed, checked in by athos, 13 years ago

Small bugs in the documentation 2.

File size: 19.6 KB
Line 
1/* -*- C++ -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library
4 *
5 * Copyright (C) 2003-2006
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_HAO_ORLIN_H
20#define LEMON_HAO_ORLIN_H
21
22#include <vector>
23#include <queue>
24#include <limits>
25
26#include <lemon/maps.h>
27#include <lemon/graph_utils.h>
28#include <lemon/graph_adaptor.h>
29#include <lemon/iterable_maps.h>
30 
31
32/// \file
33/// \ingroup flowalgs
34/// \brief Implementation of the Hao-Orlin algorithm.
35///
36/// Implementation of the HaoOrlin algorithms class for testing network
37/// reliability.
38
39namespace lemon {
40
41  /// \ingroup flowalgs
42  ///
43  /// \brief %Hao-Orlin algorithm to find a minimum cut in directed graphs.
44  ///
45  /// Hao-Orlin calculates a minimum cut in a directed graph
46  /// \f$ D=(V,A) \f$. It takes a fixed node \f$ source \in V \f$ and consists
47  /// of two phases: in the first phase it determines a minimum cut
48  /// with \f$ source \f$ on the source-side (i.e. a set \f$ X\subsetneq V \f$
49  /// with \f$ source \in X \f$ and minimal out-degree) and in the
50  /// second phase it determines a minimum cut with \f$ source \f$ on the
51  /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \notin X \f$
52  /// and minimal out-degree). Obviously, the smaller of these two
53  /// cuts will be a minimum cut of \f$ D \f$. The algorithm is a
54  /// modified push-relabel preflow algorithm and our implementation
55  /// calculates the minimum cut in \f$ O(n^3) \f$ time (we use the
56  /// highest-label rule). The purpose of such an algorithm is testing
57  /// network reliability. For an undirected graph with \f$ n \f$
58  /// nodes and \f$ e \f$ edges you can use the algorithm of Nagamochi
59  /// and Ibaraki which solves the undirected problem in
60  /// \f$ O(ne + n^2 \log(n)) \f$ time: it is implemented in the MinCut
61  /// algorithm
62  /// class.
63  ///
64  /// \param _Graph is the graph type of the algorithm.
65  /// \param _CapacityMap is an edge map of capacities which should
66  /// be any numreric type. The default type is _Graph::EdgeMap<int>.
67  /// \param _Tolerance is the handler of the inexact computation. The
68  /// default type for this is Tolerance<typename CapacityMap::Value>.
69  ///
70  /// \author Attila Bernath and Balazs Dezso
71#ifdef DOXYGEN
72  template <typename _Graph, typename _CapacityMap, typename _Tolerance>
73#else
74  template <typename _Graph,
75            typename _CapacityMap = typename _Graph::template EdgeMap<int>,
76            typename _Tolerance = Tolerance<typename _CapacityMap::Value> >
77#endif
78  class HaoOrlin {
79  protected:
80
81    typedef _Graph Graph;
82    typedef _CapacityMap CapacityMap;
83    typedef _Tolerance Tolerance;
84
85    typedef typename CapacityMap::Value Value;
86
87   
88    typedef typename Graph::Node Node;
89    typedef typename Graph::NodeIt NodeIt;
90    typedef typename Graph::EdgeIt EdgeIt;
91    typedef typename Graph::OutEdgeIt OutEdgeIt;
92    typedef typename Graph::InEdgeIt InEdgeIt;
93
94    const Graph* _graph;
95
96    const CapacityMap* _capacity;
97
98    typedef typename Graph::template EdgeMap<Value> FlowMap;
99
100    FlowMap* _preflow;
101
102    Node _source, _target;
103    int _node_num;
104
105    typedef ResGraphAdaptor<const Graph, Value, CapacityMap,
106                            FlowMap, Tolerance> OutResGraph;
107    typedef typename OutResGraph::Edge OutResEdge;
108   
109    OutResGraph* _out_res_graph;
110
111    typedef typename Graph::template NodeMap<OutResEdge> OutCurrentEdgeMap;
112    OutCurrentEdgeMap* _out_current_edge; 
113
114    typedef RevGraphAdaptor<const Graph> RevGraph;
115    RevGraph* _rev_graph;
116
117    typedef ResGraphAdaptor<const RevGraph, Value, CapacityMap,
118                            FlowMap, Tolerance> InResGraph;
119    typedef typename InResGraph::Edge InResEdge;
120   
121    InResGraph* _in_res_graph;
122
123    typedef typename Graph::template NodeMap<InResEdge> InCurrentEdgeMap;
124    InCurrentEdgeMap* _in_current_edge; 
125
126
127    typedef IterableBoolMap<Graph, Node> WakeMap;
128    WakeMap* _wake;
129
130    typedef typename Graph::template NodeMap<int> DistMap;
131    DistMap* _dist; 
132   
133    typedef typename Graph::template NodeMap<Value> ExcessMap;
134    ExcessMap* _excess;
135
136    typedef typename Graph::template NodeMap<bool> SourceSetMap;
137    SourceSetMap* _source_set;
138
139    std::vector<int> _level_size;
140
141    int _highest_active;
142    std::vector<std::list<Node> > _active_nodes;
143
144    int _dormant_max;
145    std::vector<std::list<Node> > _dormant;
146
147
148    Value _min_cut;
149
150    typedef typename Graph::template NodeMap<bool> MinCutMap;
151    MinCutMap* _min_cut_map;
152
153    Tolerance _tolerance;
154
155  public:
156
157    /// \brief Constructor
158    ///
159    /// Constructor of the algorithm class.
160    HaoOrlin(const Graph& graph, const CapacityMap& capacity,
161             const Tolerance& tolerance = Tolerance()) :
162      _graph(&graph), _capacity(&capacity),
163      _preflow(0), _source(), _target(),
164      _out_res_graph(0), _out_current_edge(0),
165      _rev_graph(0), _in_res_graph(0), _in_current_edge(0),
166      _wake(0),_dist(0), _excess(0), _source_set(0),
167      _highest_active(), _active_nodes(), _dormant_max(), _dormant(),
168      _min_cut(), _min_cut_map(0), _tolerance(tolerance) {}
169
170    ~HaoOrlin() {
171      if (_min_cut_map) {
172        delete _min_cut_map;
173      }
174      if (_in_current_edge) {
175        delete _in_current_edge;
176      }
177      if (_in_res_graph) {
178        delete _in_res_graph;
179      }
180      if (_rev_graph) {
181        delete _rev_graph;
182      }
183      if (_out_current_edge) {
184        delete _out_current_edge;
185      }
186      if (_out_res_graph) {
187        delete _out_res_graph;
188      }
189      if (_source_set) {
190        delete _source_set;
191      }
192      if (_excess) {
193        delete _excess;
194      }
195      if (_dist) {
196        delete _dist;
197      }
198      if (_wake) {
199        delete _wake;
200      }
201      if (_preflow) {
202        delete _preflow;
203      }
204    }
205   
206  private:
207   
208    template <typename ResGraph, typename EdgeMap>
209    void findMinCut(const Node& target, bool out,
210                    ResGraph& res_graph, EdgeMap& current_edge) {
211      typedef typename ResGraph::Edge ResEdge;
212      typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
213
214      for (typename Graph::EdgeIt it(*_graph); it != INVALID; ++it) {
215        (*_preflow)[it] = 0;     
216      }
217      for (NodeIt it(*_graph); it != INVALID; ++it) {
218        (*_wake)[it] = true;
219        (*_dist)[it] = 1;
220        (*_excess)[it] = 0;
221        (*_source_set)[it] = false;
222
223        res_graph.firstOut(current_edge[it], it);
224      }
225
226      _target = target;
227      (*_dist)[target] = 0;
228
229      for (ResOutEdgeIt it(res_graph, _source); it != INVALID; ++it) {
230        Value delta = res_graph.rescap(it);
231        if (!_tolerance.positive(delta)) continue;
232       
233        (*_excess)[res_graph.source(it)] -= delta;
234        res_graph.augment(it, delta);
235        Node a = res_graph.target(it);
236        if (!_tolerance.positive((*_excess)[a]) &&
237            (*_wake)[a] && a != _target) {
238          _active_nodes[(*_dist)[a]].push_front(a);
239          if (_highest_active < (*_dist)[a]) {
240            _highest_active = (*_dist)[a];
241          }
242        }
243        (*_excess)[a] += delta;
244      }
245
246      _dormant[0].push_front(_source);
247      (*_source_set)[_source] = true;
248      _dormant_max = 0;
249      (*_wake)[_source] = false;
250
251      _level_size[0] = 1;
252      _level_size[1] = _node_num - 1;
253
254      do {
255        Node n;
256        while ((n = findActiveNode()) != INVALID) {
257          ResEdge e;
258          while (_tolerance.positive((*_excess)[n]) &&
259                 (e = findAdmissibleEdge(n, res_graph, current_edge))
260                 != INVALID){
261            Value delta;
262            if ((*_excess)[n] < res_graph.rescap(e)) {
263              delta = (*_excess)[n];
264            } else {
265              delta = res_graph.rescap(e);
266              res_graph.nextOut(current_edge[n]);
267            }
268            if (!_tolerance.positive(delta)) continue;
269            res_graph.augment(e, delta);
270            (*_excess)[res_graph.source(e)] -= delta;
271            Node a = res_graph.target(e);
272            if (!_tolerance.positive((*_excess)[a]) && a != _target) {
273              _active_nodes[(*_dist)[a]].push_front(a);
274            }
275            (*_excess)[a] += delta;
276          }
277          if (_tolerance.positive((*_excess)[n])) {
278            relabel(n, res_graph, current_edge);
279          }
280        }
281
282        Value current_value = cutValue(out);
283        if (_min_cut > current_value){
284          if (out) {
285            for (NodeIt it(*_graph); it != INVALID; ++it) {
286              _min_cut_map->set(it, !(*_wake)[it]);
287            }
288          } else {
289            for (NodeIt it(*_graph); it != INVALID; ++it) {
290              _min_cut_map->set(it, (*_wake)[it]);
291            }
292          }
293
294          _min_cut = current_value;
295        }
296
297      } while (selectNewSink(res_graph));
298    }
299
300    template <typename ResGraph, typename EdgeMap>
301    void relabel(const Node& n, ResGraph& res_graph, EdgeMap& current_edge) {
302      typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
303
304      int k = (*_dist)[n];
305      if (_level_size[k] == 1) {
306        ++_dormant_max;
307        for (NodeIt it(*_graph); it != INVALID; ++it) {
308          if ((*_wake)[it] && (*_dist)[it] >= k) {
309            (*_wake)[it] = false;
310            _dormant[_dormant_max].push_front(it);
311            --_level_size[(*_dist)[it]];
312          }
313        }
314        --_highest_active;
315      } else { 
316        int new_dist = _node_num;
317        for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
318          Node t = res_graph.target(e);
319          if ((*_wake)[t] && new_dist > (*_dist)[t]) {
320            new_dist = (*_dist)[t];
321          }
322        }
323        if (new_dist == _node_num) {
324          ++_dormant_max;
325          (*_wake)[n] = false;
326          _dormant[_dormant_max].push_front(n);
327          --_level_size[(*_dist)[n]];
328        } else {           
329          --_level_size[(*_dist)[n]];
330          (*_dist)[n] = new_dist + 1;
331          _highest_active = (*_dist)[n];
332          _active_nodes[_highest_active].push_front(n);
333          ++_level_size[(*_dist)[n]];
334          res_graph.firstOut(current_edge[n], n);
335        }
336      }
337    }
338
339    template <typename ResGraph>
340    bool selectNewSink(ResGraph& res_graph) {
341      typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
342
343      Node old_target = _target;
344      (*_wake)[_target] = false;
345      --_level_size[(*_dist)[_target]];
346      _dormant[0].push_front(_target);
347      (*_source_set)[_target] = true;
348      if ((int)_dormant[0].size() == _node_num){
349        _dormant[0].clear();
350        return false;
351      }
352
353      bool wake_was_empty = false;
354
355      if(_wake->trueNum() == 0) {
356        while (!_dormant[_dormant_max].empty()){
357          (*_wake)[_dormant[_dormant_max].front()] = true;
358          ++_level_size[(*_dist)[_dormant[_dormant_max].front()]];
359          _dormant[_dormant_max].pop_front();
360        }
361        --_dormant_max;
362        wake_was_empty = true;
363      }
364
365      int min_dist = std::numeric_limits<int>::max();
366      for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
367        if (min_dist > (*_dist)[it]){
368          _target = it;
369          min_dist = (*_dist)[it];
370        }
371      }
372
373      if (wake_was_empty){
374        for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
375          if (_tolerance.positive((*_excess)[it])) {
376            if ((*_wake)[it] && it != _target) {
377              _active_nodes[(*_dist)[it]].push_front(it);
378            }
379            if (_highest_active < (*_dist)[it]) {
380              _highest_active = (*_dist)[it];               
381            }
382          }
383        }
384      }
385
386      for (ResOutEdgeIt e(res_graph, old_target); e!=INVALID; ++e){
387        if (!(*_source_set)[res_graph.target(e)]) {
388          Value delta = res_graph.rescap(e);
389          if (!_tolerance.positive(delta)) continue;
390          res_graph.augment(e, delta);
391          (*_excess)[res_graph.source(e)] -= delta;
392          Node a = res_graph.target(e);
393          if (!_tolerance.positive((*_excess)[a]) &&
394              (*_wake)[a] && a != _target) {
395            _active_nodes[(*_dist)[a]].push_front(a);
396            if (_highest_active < (*_dist)[a]) {
397              _highest_active = (*_dist)[a];
398            }
399          }
400          (*_excess)[a] += delta;
401        }
402      }
403     
404      return true;
405    }
406   
407    Node findActiveNode() {
408      while (_highest_active > 0 && _active_nodes[_highest_active].empty()){
409        --_highest_active;
410      }
411      if( _highest_active > 0) {
412        Node n = _active_nodes[_highest_active].front();
413        _active_nodes[_highest_active].pop_front();
414        return n;
415      } else {
416        return INVALID;
417      }
418    }
419
420    template <typename ResGraph, typename EdgeMap>
421    typename ResGraph::Edge findAdmissibleEdge(const Node& n,
422                                               ResGraph& res_graph,
423                                               EdgeMap& current_edge) {
424      typedef typename ResGraph::Edge ResEdge;
425      ResEdge e = current_edge[n];
426      while (e != INVALID &&
427             ((*_dist)[n] <= (*_dist)[res_graph.target(e)] ||
428              !(*_wake)[res_graph.target(e)])) {
429        res_graph.nextOut(e);
430      }
431      if (e != INVALID) {
432        current_edge[n] = e;   
433        return e;
434      } else {
435        return INVALID;
436      }
437    }
438
439    Value cutValue(bool out) {
440      Value value = 0;
441      if (out) {
442        for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
443          for (InEdgeIt e(*_graph, it); e != INVALID; ++e) {
444            if (!(*_wake)[_graph->source(e)]){
445              value += (*_capacity)[e];
446            }   
447          }
448        }
449      } else {
450        for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
451          for (OutEdgeIt e(*_graph, it); e != INVALID; ++e) {
452            if (!(*_wake)[_graph->target(e)]){
453              value += (*_capacity)[e];
454            }   
455          }
456        }
457      }
458      return value;
459    }
460
461
462  public:
463
464    /// \name Execution control
465    /// The simplest way to execute the algorithm is to use
466    /// one of the member functions called \c run(...).
467    /// \n
468    /// If you need more control on the execution,
469    /// first you must call \ref init(), then the \ref calculateIn() or
470    /// \ref calculateIn() functions.
471
472    /// @{
473
474    /// \brief Initializes the internal data structures.
475    ///
476    /// Initializes the internal data structures. It creates
477    /// the maps, residual graph adaptors and some bucket structures
478    /// for the algorithm.
479    void init() {
480      init(NodeIt(*_graph));
481    }
482
483    /// \brief Initializes the internal data structures.
484    ///
485    /// Initializes the internal data structures. It creates
486    /// the maps, residual graph adaptor and some bucket structures
487    /// for the algorithm. Node \c source  is used as the push-relabel
488    /// algorithm's source.
489    void init(const Node& source) {
490      _source = source;
491      _node_num = countNodes(*_graph);
492
493      _dormant.resize(_node_num);
494      _level_size.resize(_node_num, 0);
495      _active_nodes.resize(_node_num);
496
497      if (!_preflow) {
498        _preflow = new FlowMap(*_graph);
499      }
500      if (!_wake) {
501        _wake = new WakeMap(*_graph);
502      }
503      if (!_dist) {
504        _dist = new DistMap(*_graph);
505      }
506      if (!_excess) {
507        _excess = new ExcessMap(*_graph);
508      }
509      if (!_source_set) {
510        _source_set = new SourceSetMap(*_graph);
511      }
512      if (!_out_res_graph) {
513        _out_res_graph = new OutResGraph(*_graph, *_capacity, *_preflow);
514      }
515      if (!_out_current_edge) {
516        _out_current_edge = new OutCurrentEdgeMap(*_graph);
517      }
518      if (!_rev_graph) {
519        _rev_graph = new RevGraph(*_graph);
520      }
521      if (!_in_res_graph) {
522        _in_res_graph = new InResGraph(*_rev_graph, *_capacity, *_preflow);
523      }
524      if (!_in_current_edge) {
525        _in_current_edge = new InCurrentEdgeMap(*_graph);
526      }
527      if (!_min_cut_map) {
528        _min_cut_map = new MinCutMap(*_graph);
529      }
530
531      _min_cut = std::numeric_limits<Value>::max();
532    }
533
534
535    /// \brief Calculates a minimum cut with \f$ source \f$ on the
536    /// source-side.
537    ///
538    /// \brief Calculates a minimum cut with \f$ source \f$ on the
539    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
540    ///  and minimal out-degree).
541    void calculateOut() {
542      for (NodeIt it(*_graph); it != INVALID; ++it) {
543        if (it != _source) {
544          calculateOut(it);
545          return;
546        }
547      }
548    }
549
550    /// \brief Calculates a minimum cut with \f$ source \f$ on the
551    /// source-side.
552    ///
553    /// \brief Calculates a minimum cut with \f$ source \f$ on the
554    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
555    ///  and minimal out-degree). The \c target is the initial target
556    /// for the push-relabel algorithm.
557    void calculateOut(const Node& target) {
558      findMinCut(target, true, *_out_res_graph, *_out_current_edge);
559    }
560
561    /// \brief Calculates a minimum cut with \f$ source \f$ on the
562    /// sink-side.
563    ///
564    /// \brief Calculates a minimum cut with \f$ source \f$ on the
565    /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with
566    /// \f$ source \notin X \f$
567    /// and minimal out-degree).
568    void calculateIn() {
569      for (NodeIt it(*_graph); it != INVALID; ++it) {
570        if (it != _source) {
571          calculateIn(it);
572          return;
573        }
574      }
575    }
576
577    /// \brief Calculates a minimum cut with \f$ source \f$ on the
578    /// sink-side.
579    ///
580    /// \brief Calculates a minimum cut with \f$ source \f$ on the
581    /// sink-side (i.e. a set \f$ X\subsetneq V
582    /// \f$ with \f$ source \notin X \f$ and minimal out-degree). 
583    /// The \c target is the initial
584    /// target for the push-relabel algorithm.
585    void calculateIn(const Node& target) {
586      findMinCut(target, false, *_in_res_graph, *_in_current_edge);
587    }
588
589    /// \brief Runs the algorithm.
590    ///
591    /// Runs the algorithm. It finds nodes \c source and \c target
592    /// arbitrarily and then calls \ref init(), \ref calculateOut()
593    /// and \ref calculateIn().
594    void run() {
595      init();
596      for (NodeIt it(*_graph); it != INVALID; ++it) {
597        if (it != _source) {
598          calculateOut(it);
599          calculateIn(it);
600          return;
601        }
602      }
603    }
604
605    /// \brief Runs the algorithm.
606    ///
607    /// Runs the algorithm. It uses the given \c source node, finds a
608    /// proper \c target and then calls the \ref init(), \ref
609    /// calculateOut() and \ref calculateIn().
610    void run(const Node& s) {
611      init(s);
612      for (NodeIt it(*_graph); it != INVALID; ++it) {
613        if (it != _source) {
614          calculateOut(it);
615          calculateIn(it);
616          return;
617        }
618      }
619    }
620
621    /// \brief Runs the algorithm.
622    ///
623    /// Runs the algorithm. It just calls the \ref init() and then
624    /// \ref calculateOut() and \ref calculateIn().
625    void run(const Node& s, const Node& t) {
626      init(s);
627      calculateOut(t);
628      calculateIn(t);
629    }
630
631    /// @}
632   
633    /// \name Query Functions
634    /// The result of the %HaoOrlin algorithm
635    /// can be obtained using these functions.
636    /// \n
637    /// Before using these functions, either \ref run(), \ref
638    /// calculateOut() or \ref calculateIn() must be called.
639   
640    /// @{
641
642    /// \brief Returns the value of the minimum value cut.
643    ///
644    /// Returns the value of the minimum value cut.
645    Value minCut() const {
646      return _min_cut;
647    }
648
649
650    /// \brief Returns a minimum cut.
651    ///
652    /// Sets \c nodeMap to the characteristic vector of a minimum
653    /// value cut: it will give a nonempty set \f$ X\subsetneq V \f$
654    /// with minimal out-degree (i.e. \c nodeMap will be true exactly
655    /// for the nodes of \f$ X \f$).  \pre nodeMap should be a
656    /// bool-valued node-map.
657    template <typename NodeMap>
658    Value minCut(NodeMap& nodeMap) const {
659      for (NodeIt it(*_graph); it != INVALID; ++it) {
660        nodeMap.set(it, (*_min_cut_map)[it]);
661      }
662      return minCut();
663    }
664
665    /// @}
666   
667  }; //class HaoOrlin
668
669
670} //namespace lemon
671
672#endif //LEMON_HAO_ORLIN_H
Note: See TracBrowser for help on using the repository browser.