COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/hao_orlin.h @ 2376:0ed45a6c74b1

Last change on this file since 2376:0ed45a6c74b1 was 2376:0ed45a6c74b1, checked in by Balazs Dezso, 13 years ago

Reorganization of the modules and groups

File size: 17.9 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 <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 Graph::Node Node;
201      typedef typename ResGraph::Edge ResEdge;
202      typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
203
204      for (typename Graph::EdgeIt it(*_graph); it != INVALID; ++it) {
205        (*_preflow)[it] = 0;     
206      }
207      for (NodeIt it(*_graph); it != INVALID; ++it) {
208        (*_wake)[it] = true;
209        (*_dist)[it] = 1;
210        (*_excess)[it] = 0;
211        (*_source_set)[it] = false;
212      }
213
214      _dormant[0].push_front(_source);
215      (*_source_set)[_source] = true;
216      _dormant_max = 0;
217      (*_wake)[_source] = false;
218
219      _level_size[0] = 1;
220      _level_size[1] = _node_num - 1;
221
222      _target = target;
223      (*_dist)[target] = 0;
224
225      for (ResOutEdgeIt it(res_graph, _source); it != INVALID; ++it) {
226        Value delta = res_graph.rescap(it);
227        (*_excess)[_source] -= delta;
228        res_graph.augment(it, delta);
229        Node a = res_graph.target(it);
230        if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) {
231          _active_nodes[(*_dist)[a]].push_front(a);
232          if (_highest_active < (*_dist)[a]) {
233            _highest_active = (*_dist)[a];
234          }
235        }
236        (*_excess)[a] += delta;
237      }
238
239
240      do {
241        Node n;
242        while ((n = findActiveNode()) != INVALID) {
243          for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
244            Node a = res_graph.target(e);
245            if ((*_dist)[a] >= (*_dist)[n] || !(*_wake)[a]) continue;
246            Value delta = res_graph.rescap(e);
247            if (_tolerance.positive((*_excess)[n] - delta)) {
248              (*_excess)[n] -= delta;
249            } else {
250              delta = (*_excess)[n];
251              (*_excess)[n] = 0;
252            }
253            res_graph.augment(e, delta);
254            if ((*_excess)[a] == 0 && a != _target) {
255              _active_nodes[(*_dist)[a]].push_front(a);
256            }
257            (*_excess)[a] += delta;
258            if ((*_excess)[n] == 0) break;
259          }
260          if ((*_excess)[n] != 0) {
261            relabel(n, res_graph);
262          }
263        }
264
265        Value current_value = cutValue(out);
266        if (_min_cut > current_value){
267          if (out) {
268            for (NodeIt it(*_graph); it != INVALID; ++it) {
269              _min_cut_map->set(it, !(*_wake)[it]);
270            }
271          } else {
272            for (NodeIt it(*_graph); it != INVALID; ++it) {
273              _min_cut_map->set(it, (*_wake)[it]);
274            }
275          }
276
277          _min_cut = current_value;
278        }
279
280      } while (selectNewSink(res_graph));
281    }
282
283    template <typename ResGraph>
284    void relabel(const Node& n, ResGraph& res_graph) {
285      typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
286
287      int k = (*_dist)[n];
288      if (_level_size[k] == 1) {
289        ++_dormant_max;
290        for (NodeIt it(*_graph); it != INVALID; ++it) {
291          if ((*_wake)[it] && (*_dist)[it] >= k) {
292            (*_wake)[it] = false;
293            _dormant[_dormant_max].push_front(it);
294            --_level_size[(*_dist)[it]];
295          }
296        }
297        --_highest_active;
298      } else { 
299        int new_dist = _node_num;
300        for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
301          Node t = res_graph.target(e);
302          if ((*_wake)[t] && new_dist > (*_dist)[t]) {
303            new_dist = (*_dist)[t];
304          }
305        }
306        if (new_dist == _node_num) {
307          ++_dormant_max;
308          (*_wake)[n] = false;
309          _dormant[_dormant_max].push_front(n);
310          --_level_size[(*_dist)[n]];
311        } else {           
312          --_level_size[(*_dist)[n]];
313          (*_dist)[n] = new_dist + 1;
314          _highest_active = (*_dist)[n];
315          _active_nodes[_highest_active].push_front(n);
316          ++_level_size[(*_dist)[n]];
317        }
318      }
319    }
320
321    template <typename ResGraph>
322    bool selectNewSink(ResGraph& res_graph) {
323      typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
324
325      Node old_target = _target;
326      (*_wake)[_target] = false;
327      --_level_size[(*_dist)[_target]];
328      _dormant[0].push_front(_target);
329      (*_source_set)[_target] = true;
330      if ((int)_dormant[0].size() == _node_num){
331        _dormant[0].clear();
332        return false;
333      }
334
335      bool wake_was_empty = false;
336
337      if(_wake->trueNum() == 0) {
338        while (!_dormant[_dormant_max].empty()){
339          (*_wake)[_dormant[_dormant_max].front()] = true;
340          ++_level_size[(*_dist)[_dormant[_dormant_max].front()]];
341          _dormant[_dormant_max].pop_front();
342        }
343        --_dormant_max;
344        wake_was_empty = true;
345      }
346
347      int min_dist = std::numeric_limits<int>::max();
348      for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
349        if (min_dist > (*_dist)[it]){
350          _target = it;
351          min_dist = (*_dist)[it];
352        }
353      }
354
355      if (wake_was_empty){
356        for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
357          if ((*_excess)[it] != 0 && it != _target) {
358            _active_nodes[(*_dist)[it]].push_front(it);
359            if (_highest_active < (*_dist)[it]) {
360              _highest_active = (*_dist)[it];               
361            }
362          }
363        }
364      }
365
366      Node n = old_target;
367      for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e){
368        Node a = res_graph.target(e);
369        if (!(*_wake)[a]) continue;
370        Value delta = res_graph.rescap(e);
371        res_graph.augment(e, delta);
372        (*_excess)[n] -= delta;
373        if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) {
374          _active_nodes[(*_dist)[a]].push_front(a);
375          if (_highest_active < (*_dist)[a]) {
376            _highest_active = (*_dist)[a];
377          }
378        }
379        (*_excess)[a] += delta;
380      }
381     
382      return true;
383    }
384
385    Node findActiveNode() {
386      while (_highest_active > 0 && _active_nodes[_highest_active].empty()){
387        --_highest_active;
388      }
389      if( _highest_active > 0) {
390        Node n = _active_nodes[_highest_active].front();
391        _active_nodes[_highest_active].pop_front();
392        return n;
393      } else {
394        return INVALID;
395      }
396    }
397
398    Value cutValue(bool out) {
399      Value value = 0;
400      if (out) {
401        for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
402          for (InEdgeIt e(*_graph, it); e != INVALID; ++e) {
403            if (!(*_wake)[_graph->source(e)]){
404              value += (*_capacity)[e];
405            }   
406          }
407        }
408      } else {
409        for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
410          for (OutEdgeIt e(*_graph, it); e != INVALID; ++e) {
411            if (!(*_wake)[_graph->target(e)]){
412              value += (*_capacity)[e];
413            }   
414          }
415        }
416      }
417      return value;
418    }
419
420
421  public:
422
423    /// \name Execution control
424    /// The simplest way to execute the algorithm is to use
425    /// one of the member functions called \c run(...).
426    /// \n
427    /// If you need more control on the execution,
428    /// first you must call \ref init(), then the \ref calculateIn() or
429    /// \ref calculateIn() functions.
430
431    /// @{
432
433    /// \brief Initializes the internal data structures.
434    ///
435    /// Initializes the internal data structures. It creates
436    /// the maps, residual graph adaptors and some bucket structures
437    /// for the algorithm.
438    void init() {
439      init(NodeIt(*_graph));
440    }
441
442    /// \brief Initializes the internal data structures.
443    ///
444    /// Initializes the internal data structures. It creates
445    /// the maps, residual graph adaptor and some bucket structures
446    /// for the algorithm. Node \c source  is used as the push-relabel
447    /// algorithm's source.
448    void init(const Node& source) {
449      _source = source;
450      _node_num = countNodes(*_graph);
451
452      _dormant.resize(_node_num);
453      _level_size.resize(_node_num, 0);
454      _active_nodes.resize(_node_num);
455
456      if (!_preflow) {
457        _preflow = new FlowMap(*_graph);
458      }
459      if (!_wake) {
460        _wake = new WakeMap(*_graph);
461      }
462      if (!_dist) {
463        _dist = new DistMap(*_graph);
464      }
465      if (!_excess) {
466        _excess = new ExcessMap(*_graph);
467      }
468      if (!_source_set) {
469        _source_set = new SourceSetMap(*_graph);
470      }
471      if (!_out_res_graph) {
472        _out_res_graph = new OutResGraph(*_graph, *_capacity, *_preflow);
473      }
474      if (!_rev_graph) {
475        _rev_graph = new RevGraph(*_graph);
476      }
477      if (!_in_res_graph) {
478        _in_res_graph = new InResGraph(*_rev_graph, *_capacity, *_preflow);
479      }
480      if (!_min_cut_map) {
481        _min_cut_map = new MinCutMap(*_graph);
482      }
483
484      _min_cut = std::numeric_limits<Value>::max();
485    }
486
487
488    /// \brief Calculates a minimum cut with \f$ source \f$ on the
489    /// source-side.
490    ///
491    /// \brief Calculates a minimum cut with \f$ source \f$ on the
492    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
493    ///  and minimal out-degree).
494    void calculateOut() {
495      for (NodeIt it(*_graph); it != INVALID; ++it) {
496        if (it != _source) {
497          calculateOut(it);
498          return;
499        }
500      }
501    }
502
503    /// \brief Calculates a minimum cut with \f$ source \f$ on the
504    /// source-side.
505    ///
506    /// \brief Calculates a minimum cut with \f$ source \f$ on the
507    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
508    ///  and minimal out-degree). The \c target is the initial target
509    /// for the push-relabel algorithm.
510    void calculateOut(const Node& target) {
511      findMinCut(target, true, *_out_res_graph);
512    }
513
514    /// \brief Calculates a minimum cut with \f$ source \f$ on the
515    /// sink-side.
516    ///
517    /// \brief Calculates a minimum cut with \f$ source \f$ on the
518    /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with
519    /// \f$ source \notin X \f$
520    /// and minimal out-degree).
521    void calculateIn() {
522      for (NodeIt it(*_graph); it != INVALID; ++it) {
523        if (it != _source) {
524          calculateIn(it);
525          return;
526        }
527      }
528    }
529
530    /// \brief Calculates a minimum cut with \f$ source \f$ on the
531    /// sink-side.
532    ///
533    /// \brief Calculates a minimum cut with \f$ source \f$ on the
534    /// sink-side (i.e. a set \f$ X\subsetneq V
535    /// \f$ with \f$ source \notin X \f$ and minimal out-degree). 
536    /// The \c target is the initial
537    /// target for the push-relabel algorithm.
538    void calculateIn(const Node& target) {
539      findMinCut(target, false, *_in_res_graph);
540    }
541
542    /// \brief Runs the algorithm.
543    ///
544    /// Runs the algorithm. It finds nodes \c source and \c target
545    /// arbitrarily and then calls \ref init(), \ref calculateOut()
546    /// and \ref calculateIn().
547    void run() {
548      init();
549      for (NodeIt it(*_graph); it != INVALID; ++it) {
550        if (it != _source) {
551          calculateOut(it);
552          calculateIn(it);
553          return;
554        }
555      }
556    }
557
558    /// \brief Runs the algorithm.
559    ///
560    /// Runs the algorithm. It uses the given \c source node, finds a
561    /// proper \c target and then calls the \ref init(), \ref
562    /// calculateOut() and \ref calculateIn().
563    void run(const Node& s) {
564      init(s);
565      for (NodeIt it(*_graph); it != INVALID; ++it) {
566        if (it != _source) {
567          calculateOut(it);
568          calculateIn(it);
569          return;
570        }
571      }
572    }
573
574    /// \brief Runs the algorithm.
575    ///
576    /// Runs the algorithm. It just calls the \ref init() and then
577    /// \ref calculateOut() and \ref calculateIn().
578    void run(const Node& s, const Node& t) {
579      init(s);
580      calculateOut(t);
581      calculateIn(t);
582    }
583
584    /// @}
585   
586    /// \name Query Functions
587    /// The result of the %HaoOrlin algorithm
588    /// can be obtained using these functions.
589    /// \n
590    /// Before using these functions, either \ref run(), \ref
591    /// calculateOut() or \ref calculateIn() must be called.
592   
593    /// @{
594
595    /// \brief Returns the value of the minimum value cut.
596    ///
597    /// Returns the value of the minimum value cut.
598    Value minCut() const {
599      return _min_cut;
600    }
601
602
603    /// \brief Returns a minimum cut.
604    ///
605    /// Sets \c nodeMap to the characteristic vector of a minimum
606    /// value cut: it will give a nonempty set \f$ X\subsetneq V \f$
607    /// with minimal out-degree (i.e. \c nodeMap will be true exactly
608    /// for the nodes of \f$ X \f$).  \pre nodeMap should be a
609    /// bool-valued node-map.
610    template <typename NodeMap>
611    Value minCut(NodeMap& nodeMap) const {
612      for (NodeIt it(*_graph); it != INVALID; ++it) {
613        nodeMap.set(it, (*_min_cut_map)[it]);
614      }
615      return minCut();
616    }
617
618    /// @}
619   
620  }; //class HaoOrlin
621
622
623} //namespace lemon
624
625#endif //LEMON_HAO_ORLIN_H
Note: See TracBrowser for help on using the repository browser.