COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/hao_orlin.h @ 2476:059dcdda37c5

Last change on this file since 2476:059dcdda37c5 was 2391:14a343be7a5a, checked in by Alpar Juttner, 13 years ago

Happy New Year to all source files!

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