COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/hao_orlin.h @ 2228:f71b0f9a7c3a

Last change on this file since 2228:f71b0f9a7c3a was 2228:f71b0f9a7c3a, checked in by athos, 17 years ago

Improved documentation.

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