COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/hao_orlin.h @ 2225:bb3d5e6f9fcb

Last change on this file since 2225:bb3d5e6f9fcb was 2225:bb3d5e6f9fcb, checked in by Balazs Dezso, 13 years ago

Doc fix

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