COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/hao_orlin.h @ 2211:c790d04e192a

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

Hao-Orlin algorithm

It is based on Attila's work
It is tested on all dimacs files in data directory

It may need more execution control

  • possible interruption after each findNewSink
File size: 13.8 KB
Line 
1/* -*- C++ -*-
2 * lemon/hao_orlin.h - Part of LEMON, a generic C++ optimization library
3 *
4 * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
5 * (Egervary Research Group on Combinatorial Optimization, EGRES).
6 *
7 * Permission to use, modify and distribute this software is granted
8 * provided that this copyright notice appears in all copies. For
9 * precise terms see the accompanying LICENSE file.
10 *
11 * This software is provided "AS IS" with no warranty of any kind,
12 * express or implied, and with no claim as to its suitability for any
13 * purpose.
14 *
15 */
16
17#ifndef LEMON_HAO_ORLIN_H
18#define LEMON_HAO_ORLIN_H
19
20#include <vector>
21#include <queue>
22#include <limits>
23
24#include <lemon/maps.h>
25#include <lemon/graph_utils.h>
26#include <lemon/graph_adaptor.h>
27#include <lemon/iterable_maps.h>
28 
29
30/// \file
31/// \ingroup flowalgs
32/// Implementation of the Hao-Orlin algorithms class for testing network
33/// reliability.
34
35namespace lemon {
36
37  /// \addtogroup flowalgs
38  /// @{                                                   
39
40  /// %Hao-Orlin algorithm for calculate minimum cut in directed graphs.
41  ///
42  /// Hao-Orlin calculates the minimum cut in directed graphs. It
43  /// separates the nodes of the graph into two disjoint sets and the
44  /// summary of the edge capacities go from the first set to the
45  /// second set is the minimum.  The algorithm is a modified
46  /// push-relabel preflow algorithm and it calculates the minimum cat
47  /// in \f$ O(n^3) \f$ time. The purpose of such algorithm is testing
48  /// network reliability. For sparse undirected graph you can use the
49  /// algorithm of Nagamochi and Ibraki which solves the undirected
50  /// problem in \f$ O(n^3) \f$ time.
51  ///
52  /// \author Attila Bernath and Balazs Dezso
53  template <typename _Graph,
54            typename _CapacityMap = typename _Graph::template EdgeMap<int>,
55            typename _Tolerance = Tolerance<typename _CapacityMap::Value> >
56  class HaoOrlin {
57  protected:
58
59    typedef _Graph Graph;
60    typedef _CapacityMap CapacityMap;
61    typedef _Tolerance Tolerance;
62
63    typedef typename CapacityMap::Value Value;
64
65   
66    typedef typename Graph::Node Node;
67    typedef typename Graph::NodeIt NodeIt;
68    typedef typename Graph::EdgeIt EdgeIt;
69    typedef typename Graph::OutEdgeIt OutEdgeIt;
70    typedef typename Graph::InEdgeIt InEdgeIt;
71
72    const Graph* _graph;
73    const CapacityMap* _capacity;
74
75    typedef typename Graph::template EdgeMap<Value> FlowMap;
76
77    FlowMap* _preflow;
78
79    Node _source, _target;
80    int _node_num;
81
82    typedef ResGraphAdaptor<const Graph, Value, CapacityMap,
83                            FlowMap, Tolerance> ResGraph;
84    typedef typename ResGraph::Node ResNode;
85    typedef typename ResGraph::NodeIt ResNodeIt;
86    typedef typename ResGraph::EdgeIt ResEdgeIt;
87    typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
88    typedef typename ResGraph::Edge ResEdge;
89    typedef typename ResGraph::InEdgeIt ResInEdgeIt;
90
91    ResGraph* _res_graph;
92
93    typedef typename Graph::template NodeMap<ResEdge> CurrentArcMap;
94    CurrentArcMap* _current_arc; 
95
96
97    typedef IterableBoolMap<Graph, Node> WakeMap;
98    WakeMap* _wake;
99
100    typedef typename Graph::template NodeMap<int> DistMap;
101    DistMap* _dist; 
102   
103    typedef typename Graph::template NodeMap<Value> ExcessMap;
104    ExcessMap* _excess;
105
106    typedef typename Graph::template NodeMap<bool> SourceSetMap;
107    SourceSetMap* _source_set;
108
109    std::vector<int> _level_size;
110
111    int _highest_active;
112    std::vector<std::list<Node> > _active_nodes;
113
114    int _dormant_max;
115    std::vector<std::list<Node> > _dormant;
116
117
118    Value _min_cut;
119
120    typedef typename Graph::template NodeMap<bool> MinCutMap;
121    MinCutMap* _min_cut_map;
122
123    Tolerance _tolerance;
124
125  public:
126
127    HaoOrlin(const Graph& graph, const CapacityMap& capacity,
128             const Tolerance& tolerance = Tolerance()) :
129      _graph(&graph), _capacity(&capacity),
130      _preflow(0), _source(), _target(), _res_graph(0), _current_arc(0),
131      _wake(0),_dist(0), _excess(0), _source_set(0),
132      _highest_active(), _active_nodes(), _dormant_max(), _dormant(),
133      _min_cut(), _min_cut_map(0), _tolerance(tolerance) {}
134
135    ~HaoOrlin() {
136      if (_res_graph) {
137        delete _res_graph;
138      }
139      if (_min_cut_map) {
140        delete _min_cut_map;
141      }
142      if (_current_arc) {
143        delete _current_arc;
144      }
145      if (_source_set) {
146        delete _source_set;
147      }
148      if (_excess) {
149        delete _excess;
150      }
151      if (_dist) {
152        delete _dist;
153      }
154      if (_wake) {
155        delete _wake;
156      }
157      if (_preflow) {
158        delete _preflow;
159      }
160    }
161   
162  private:
163   
164    void relabel(Node i) {
165      int k = (*_dist)[i];
166      if (_level_size[k] == 1) {
167        ++_dormant_max;
168        for (NodeIt it(*_graph); it != INVALID; ++it) {
169          if ((*_wake)[it] && (*_dist)[it] >= k) {
170            (*_wake)[it] = false;
171            _dormant[_dormant_max].push_front(it);
172            --_level_size[(*_dist)[it]];
173          }
174        }
175        --_highest_active;
176      } else {
177        ResOutEdgeIt e(*_res_graph, i);
178        while (e != INVALID && !(*_wake)[_res_graph->target(e)]) {
179          ++e;
180        }
181
182        if (e == INVALID){
183          ++_dormant_max;
184          (*_wake)[i] = false;
185          _dormant[_dormant_max].push_front(i);
186          --_level_size[(*_dist)[i]];
187        } else{     
188          Node j = _res_graph->target(e);
189          int new_dist = (*_dist)[j];
190          ++e;
191          while (e != INVALID){
192            Node j = _res_graph->target(e);
193            if ((*_wake)[j] && new_dist > (*_dist)[j]) {
194              new_dist = (*_dist)[j];
195            }
196            ++e;
197          }
198          --_level_size[(*_dist)[i]];
199          (*_dist)[i] = new_dist + 1;
200          _highest_active = (*_dist)[i];
201          _active_nodes[_highest_active].push_front(i);
202          ++_level_size[(*_dist)[i]];
203          _res_graph->firstOut((*_current_arc)[i], i);
204        }
205      }
206    }
207
208    bool selectNewSink(){
209      Node old_target = _target;
210      (*_wake)[_target] = false;
211      --_level_size[(*_dist)[_target]];
212      _dormant[0].push_front(_target);
213      (*_source_set)[_target] = true;
214      if ((int)_dormant[0].size() == _node_num){
215        _dormant[0].clear();
216        return false;
217      }
218
219      bool wake_was_empty = false;
220
221      if(_wake->trueNum() == 0) {
222        while (!_dormant[_dormant_max].empty()){
223          (*_wake)[_dormant[_dormant_max].front()] = true;
224          ++_level_size[(*_dist)[_dormant[_dormant_max].front()]];
225          _dormant[_dormant_max].pop_front();
226        }
227        --_dormant_max;
228        wake_was_empty = true;
229      }
230
231      int min_dist = std::numeric_limits<int>::max();
232      for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
233        if (min_dist > (*_dist)[it]){
234          _target = it;
235          min_dist = (*_dist)[it];
236        }
237      }
238
239      if (wake_was_empty){
240        for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
241          if (_tolerance.positive((*_excess)[it])) {
242            if ((*_wake)[it] && it != _target) {
243              _active_nodes[(*_dist)[it]].push_front(it);
244            }
245            if (_highest_active < (*_dist)[it]) {
246              _highest_active = (*_dist)[it];               
247            }
248          }
249        }
250      }
251
252      for (ResOutEdgeIt e(*_res_graph, old_target); e!=INVALID; ++e){
253        if (!(*_source_set)[_res_graph->target(e)]){
254          push(e, _res_graph->rescap(e));
255        }
256      }
257     
258      return true;
259    }
260   
261    Node findActiveNode() {
262      while (_highest_active > 0 && _active_nodes[_highest_active].empty()){
263        --_highest_active;
264      }
265      if( _highest_active > 0) {
266        Node n = _active_nodes[_highest_active].front();
267        _active_nodes[_highest_active].pop_front();
268        return n;
269      } else {
270        return INVALID;
271      }
272    }
273
274    ResEdge findAdmissibleEdge(const Node& n){
275      ResEdge e = (*_current_arc)[n];
276      while (e != INVALID &&
277             ((*_dist)[n] <= (*_dist)[_res_graph->target(e)] ||
278              !(*_wake)[_res_graph->target(e)])) {
279        _res_graph->nextOut(e);
280      }
281      if (e != INVALID) {
282        (*_current_arc)[n] = e;
283        return e;
284      } else {
285        return INVALID;
286      }
287    }
288
289    void push(ResEdge& e,const Value& delta){
290      _res_graph->augment(e, delta);
291      if (!_tolerance.positive(delta)) return;
292     
293      (*_excess)[_res_graph->source(e)] -= delta;
294      Node a = _res_graph->target(e);
295      if (!_tolerance.positive((*_excess)[a]) && (*_wake)[a] && a != _target) {
296        _active_nodes[(*_dist)[a]].push_front(a);
297        if (_highest_active < (*_dist)[a]) {
298          _highest_active = (*_dist)[a];
299        }
300      }
301      (*_excess)[a] += delta;
302    }
303   
304    Value cutValue() {
305      Value value = 0;
306      for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
307        for (InEdgeIt e(*_graph, it); e != INVALID; ++e) {
308          if (!(*_wake)[_graph->source(e)]){
309            value += (*_capacity)[e];
310          }     
311        }
312      }
313      return value;
314    }   
315
316  public:
317
318    /// \brief Initializes the internal data structures.
319    ///
320    /// Initializes the internal data structures. It creates
321    /// the maps, residual graph adaptor and some bucket structures
322    /// for the algorithm.
323    void init() {
324      init(NodeIt(*_graph));
325    }
326
327    /// \brief Initializes the internal data structures.
328    ///
329    /// Initializes the internal data structures. It creates
330    /// the maps, residual graph adaptor and some bucket structures
331    /// for the algorithm. The \c source node is used as the push-relabel
332    /// algorithm's source.
333    void init(const Node& source) {
334      _source = source;
335      _node_num = countNodes(*_graph);
336
337      _dormant.resize(_node_num);
338      _level_size.resize(_node_num, 0);
339      _active_nodes.resize(_node_num);
340
341      if (!_preflow) {
342        _preflow = new FlowMap(*_graph);
343      }
344      if (!_wake) {
345        _wake = new WakeMap(*_graph);
346      }
347      if (!_dist) {
348        _dist = new DistMap(*_graph);
349      }
350      if (!_excess) {
351        _excess = new ExcessMap(*_graph);
352      }
353      if (!_source_set) {
354        _source_set = new SourceSetMap(*_graph);
355      }
356      if (!_current_arc) {
357        _current_arc = new CurrentArcMap(*_graph);
358      }
359      if (!_min_cut_map) {
360        _min_cut_map = new MinCutMap(*_graph);
361      }
362      if (!_res_graph) {
363        _res_graph = new ResGraph(*_graph, *_capacity, *_preflow);
364      }
365
366      _min_cut = std::numeric_limits<Value>::max();
367    }
368
369
370    /// \brief Calculates the minimum cut with the \c source node
371    /// in the first partition.
372    ///
373    /// Calculates the minimum cut with the \c source node
374    /// in the first partition.
375    void calculateOut() {
376      for (NodeIt it(*_graph); it != INVALID; ++it) {
377        if (it != _source) {
378          calculateOut(it);
379          return;
380        }
381      }
382    }
383
384    /// \brief Calculates the minimum cut with the \c source node
385    /// in the first partition.
386    ///
387    /// Calculates the minimum cut with the \c source node
388    /// in the first partition. The \c target is the initial target
389    /// for the push-relabel algorithm.
390    void calculateOut(const Node& target) {
391      for (NodeIt it(*_graph); it != INVALID; ++it) {
392        (*_wake)[it] = true;
393        (*_dist)[it] = 1;
394        (*_excess)[it] = 0;
395        (*_source_set)[it] = false;
396
397        _res_graph->firstOut((*_current_arc)[it], it);
398      }
399
400      _target = target;
401      (*_dist)[target] = 0;
402
403      for (ResOutEdgeIt it(*_res_graph, _source); it != INVALID; ++it) {
404        push(it, _res_graph->rescap(it));
405      }
406
407      _dormant[0].push_front(_source);
408      (*_source_set)[_source] = true;
409      _dormant_max = 0;
410      (*_wake)[_source]=false;
411
412      _level_size[0] = 1;
413      _level_size[1] = _node_num - 1;
414
415      do {
416        Node n;
417        while ((n = findActiveNode()) != INVALID) {
418          ResEdge e;
419          while (_tolerance.positive((*_excess)[n]) &&
420                 (e = findAdmissibleEdge(n)) != INVALID){
421            Value delta;
422            if ((*_excess)[n] < _res_graph->rescap(e)) {
423              delta = (*_excess)[n];
424            } else {
425              delta = _res_graph->rescap(e);
426              _res_graph->nextOut((*_current_arc)[n]);
427            }
428            if (!_tolerance.positive(delta)) continue;
429            _res_graph->augment(e, delta);
430            (*_excess)[_res_graph->source(e)] -= delta;
431            Node a = _res_graph->target(e);
432            if (!_tolerance.positive((*_excess)[a]) && a != _target) {
433              _active_nodes[(*_dist)[a]].push_front(a);
434            }
435            (*_excess)[a] += delta;
436          }
437          if (_tolerance.positive((*_excess)[n])) {
438            relabel(n);
439          }
440        }
441
442        Value current_value = cutValue();
443        if (_min_cut > current_value){
444          for (NodeIt it(*_graph); it != INVALID; ++it) {
445            _min_cut_map->set(it, !(*_wake)[it]);
446          }
447
448          _min_cut = current_value;
449        }
450
451      } while (selectNewSink());
452    }
453
454    void calculateIn() {
455      for (NodeIt it(*_graph); it != INVALID; ++it) {
456        if (it != _source) {
457          calculateIn(it);
458          return;
459        }
460      }
461    }
462
463    void run() {
464      init();
465      for (NodeIt it(*_graph); it != INVALID; ++it) {
466        if (it != _source) {
467          startOut(it);
468          //          startIn(it);
469          return;
470        }
471      }
472    }
473
474    void run(const Node& s) {
475      init(s);
476      for (NodeIt it(*_graph); it != INVALID; ++it) {
477        if (it != _source) {
478          startOut(it);
479          //          startIn(it);
480          return;
481        }
482      }
483    }
484
485    void run(const Node& s, const Node& t) {
486      init(s);
487      startOut(t);
488      startIn(t);
489    }
490   
491    /// \brief Returns the value of the minimum value cut with node \c
492    /// source on the source side.
493    ///
494    /// Returns the value of the minimum value cut with node \c source
495    /// on the source side. This function can be called after running
496    /// \ref findMinCut.
497    Value minCut() const {
498      return _min_cut;
499    }
500
501
502    /// \brief Returns a minimum value cut.
503    ///
504    /// Sets \c nodeMap to the characteristic vector of a minimum
505    /// value cut with node \c source on the source side. This
506    /// function can be called after running \ref findMinCut. 
507    /// \pre nodeMap should be a bool-valued node-map.
508    template <typename NodeMap>
509    Value minCut(NodeMap& nodeMap) const {
510      for (NodeIt it(*_graph); it != INVALID; ++it) {
511        nodeMap.set(it, (*_min_cut_map)[it]);
512      }
513      return minCut();
514    }
515   
516  }; //class HaoOrlin
517
518
519} //namespace lemon
520
521#endif //LEMON_HAO_ORLIN_H
Note: See TracBrowser for help on using the repository browser.