COIN-OR::LEMON - Graph Library

source: lemon/lemon/hao_orlin.h @ 451:09e416d35896

Last change on this file since 451:09e416d35896 was 428:7030149efed2, checked in by Alpar Juttner <alpar@…>, 15 years ago

Minor doc improvements in HaoOrlin? (#58)

File size: 29.5 KB
Line 
1/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library.
4 *
5 * Copyright (C) 2003-2008
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 <list>
24#include <limits>
25
26#include <lemon/maps.h>
27#include <lemon/core.h>
28#include <lemon/tolerance.h>
29
30/// \file
31/// \ingroup min_cut
32/// \brief Implementation of the Hao-Orlin algorithm.
33///
34/// Implementation of the Hao-Orlin algorithm class for testing network
35/// reliability.
36
37namespace lemon {
38
39  /// \ingroup min_cut
40  ///
41  /// \brief %Hao-Orlin algorithm to find a minimum cut in directed graphs.
42  ///
43  /// Hao-Orlin calculates a minimum cut in a directed graph
44  /// \f$D=(V,A)\f$. It takes a fixed node \f$ source \in V \f$ and
45  /// consists of two phases: in the first phase it determines a
46  /// minimum cut with \f$ source \f$ on the source-side (i.e. a set
47  /// \f$ X\subsetneq V \f$ with \f$ source \in X \f$ and minimal
48  /// out-degree) and in the second phase it determines a minimum cut
49  /// with \f$ source \f$ on the sink-side (i.e. a set
50  /// \f$ X\subsetneq V \f$ with \f$ source \notin X \f$ and minimal
51  /// out-degree). Obviously, the smaller of these two cuts will be a
52  /// minimum cut of \f$ D \f$. The algorithm is a modified
53  /// push-relabel preflow algorithm and our implementation calculates
54  /// the minimum cut in \f$ O(n^2\sqrt{m}) \f$ time (we use the
55  /// highest-label rule), or in \f$O(nm)\f$ for unit capacities. The
56  /// purpose of such algorithm is testing network reliability. For an
57  /// undirected graph you can run just the first phase of the
58  /// algorithm or you can use the algorithm of Nagamochi and Ibaraki
59  /// which solves the undirected problem in
60  /// \f$ O(nm + n^2 \log(n)) \f$ time: it is implemented in the
61  /// NagamochiIbaraki algorithm class.
62  ///
63  /// \param _Digraph 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 _Digraph::ArcMap<int>.
66  /// \param _Tolerance is the handler of the inexact computation. The
67  /// default type for this is Tolerance<CapacityMap::Value>.
68#ifdef DOXYGEN
69  template <typename _Digraph, typename _CapacityMap, typename _Tolerance>
70#else
71  template <typename _Digraph,
72            typename _CapacityMap = typename _Digraph::template ArcMap<int>,
73            typename _Tolerance = Tolerance<typename _CapacityMap::Value> >
74#endif
75  class HaoOrlin {
76  private:
77
78    typedef _Digraph Digraph;
79    typedef _CapacityMap CapacityMap;
80    typedef _Tolerance Tolerance;
81
82    typedef typename CapacityMap::Value Value;
83
84    TEMPLATE_GRAPH_TYPEDEFS(Digraph);
85
86    const Digraph& _graph;
87    const CapacityMap* _capacity;
88
89    typedef typename Digraph::template ArcMap<Value> FlowMap;
90    FlowMap* _flow;
91
92    Node _source;
93
94    int _node_num;
95
96    // Bucketing structure
97    std::vector<Node> _first, _last;
98    typename Digraph::template NodeMap<Node>* _next;
99    typename Digraph::template NodeMap<Node>* _prev;
100    typename Digraph::template NodeMap<bool>* _active;
101    typename Digraph::template NodeMap<int>* _bucket;
102
103    std::vector<bool> _dormant;
104
105    std::list<std::list<int> > _sets;
106    std::list<int>::iterator _highest;
107
108    typedef typename Digraph::template NodeMap<Value> ExcessMap;
109    ExcessMap* _excess;
110
111    typedef typename Digraph::template NodeMap<bool> SourceSetMap;
112    SourceSetMap* _source_set;
113
114    Value _min_cut;
115
116    typedef typename Digraph::template NodeMap<bool> MinCutMap;
117    MinCutMap* _min_cut_map;
118
119    Tolerance _tolerance;
120
121  public:
122
123    /// \brief Constructor
124    ///
125    /// Constructor of the algorithm class.
126    HaoOrlin(const Digraph& graph, const CapacityMap& capacity,
127             const Tolerance& tolerance = Tolerance()) :
128      _graph(graph), _capacity(&capacity), _flow(0), _source(),
129      _node_num(), _first(), _last(), _next(0), _prev(0),
130      _active(0), _bucket(0), _dormant(), _sets(), _highest(),
131      _excess(0), _source_set(0), _min_cut(), _min_cut_map(0),
132      _tolerance(tolerance) {}
133
134    ~HaoOrlin() {
135      if (_min_cut_map) {
136        delete _min_cut_map;
137      }
138      if (_source_set) {
139        delete _source_set;
140      }
141      if (_excess) {
142        delete _excess;
143      }
144      if (_next) {
145        delete _next;
146      }
147      if (_prev) {
148        delete _prev;
149      }
150      if (_active) {
151        delete _active;
152      }
153      if (_bucket) {
154        delete _bucket;
155      }
156      if (_flow) {
157        delete _flow;
158      }
159    }
160
161  private:
162
163    void activate(const Node& i) {
164      _active->set(i, true);
165
166      int bucket = (*_bucket)[i];
167
168      if ((*_prev)[i] == INVALID || (*_active)[(*_prev)[i]]) return;
169      //unlace
170      _next->set((*_prev)[i], (*_next)[i]);
171      if ((*_next)[i] != INVALID) {
172        _prev->set((*_next)[i], (*_prev)[i]);
173      } else {
174        _last[bucket] = (*_prev)[i];
175      }
176      //lace
177      _next->set(i, _first[bucket]);
178      _prev->set(_first[bucket], i);
179      _prev->set(i, INVALID);
180      _first[bucket] = i;
181    }
182
183    void deactivate(const Node& i) {
184      _active->set(i, false);
185      int bucket = (*_bucket)[i];
186
187      if ((*_next)[i] == INVALID || !(*_active)[(*_next)[i]]) return;
188
189      //unlace
190      _prev->set((*_next)[i], (*_prev)[i]);
191      if ((*_prev)[i] != INVALID) {
192        _next->set((*_prev)[i], (*_next)[i]);
193      } else {
194        _first[bucket] = (*_next)[i];
195      }
196      //lace
197      _prev->set(i, _last[bucket]);
198      _next->set(_last[bucket], i);
199      _next->set(i, INVALID);
200      _last[bucket] = i;
201    }
202
203    void addItem(const Node& i, int bucket) {
204      (*_bucket)[i] = bucket;
205      if (_last[bucket] != INVALID) {
206        _prev->set(i, _last[bucket]);
207        _next->set(_last[bucket], i);
208        _next->set(i, INVALID);
209        _last[bucket] = i;
210      } else {
211        _prev->set(i, INVALID);
212        _first[bucket] = i;
213        _next->set(i, INVALID);
214        _last[bucket] = i;
215      }
216    }
217
218    void findMinCutOut() {
219
220      for (NodeIt n(_graph); n != INVALID; ++n) {
221        _excess->set(n, 0);
222      }
223
224      for (ArcIt a(_graph); a != INVALID; ++a) {
225        _flow->set(a, 0);
226      }
227
228      int bucket_num = 0;
229      std::vector<Node> queue(_node_num);
230      int qfirst = 0, qlast = 0, qsep = 0;
231
232      {
233        typename Digraph::template NodeMap<bool> reached(_graph, false);
234
235        reached.set(_source, true);
236        bool first_set = true;
237
238        for (NodeIt t(_graph); t != INVALID; ++t) {
239          if (reached[t]) continue;
240          _sets.push_front(std::list<int>());
241         
242          queue[qlast++] = t;
243          reached.set(t, true);
244
245          while (qfirst != qlast) {
246            if (qsep == qfirst) {
247              ++bucket_num;
248              _sets.front().push_front(bucket_num);
249              _dormant[bucket_num] = !first_set;
250              _first[bucket_num] = _last[bucket_num] = INVALID;
251              qsep = qlast;
252            }
253
254            Node n = queue[qfirst++];
255            addItem(n, bucket_num);
256
257            for (InArcIt a(_graph, n); a != INVALID; ++a) {
258              Node u = _graph.source(a);
259              if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
260                reached.set(u, true);
261                queue[qlast++] = u;
262              }
263            }
264          }
265          first_set = false;
266        }
267
268        ++bucket_num;
269        _bucket->set(_source, 0);
270        _dormant[0] = true;
271      }
272      _source_set->set(_source, true);
273
274      Node target = _last[_sets.back().back()];
275      {
276        for (OutArcIt a(_graph, _source); a != INVALID; ++a) {
277          if (_tolerance.positive((*_capacity)[a])) {
278            Node u = _graph.target(a);
279            _flow->set(a, (*_capacity)[a]);
280            _excess->set(u, (*_excess)[u] + (*_capacity)[a]);
281            if (!(*_active)[u] && u != _source) {
282              activate(u);
283            }
284          }
285        }
286
287        if ((*_active)[target]) {
288          deactivate(target);
289        }
290
291        _highest = _sets.back().begin();
292        while (_highest != _sets.back().end() &&
293               !(*_active)[_first[*_highest]]) {
294          ++_highest;
295        }
296      }
297
298      while (true) {
299        while (_highest != _sets.back().end()) {
300          Node n = _first[*_highest];
301          Value excess = (*_excess)[n];
302          int next_bucket = _node_num;
303
304          int under_bucket;
305          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
306            under_bucket = -1;
307          } else {
308            under_bucket = *(++std::list<int>::iterator(_highest));
309          }
310
311          for (OutArcIt a(_graph, n); a != INVALID; ++a) {
312            Node v = _graph.target(a);
313            if (_dormant[(*_bucket)[v]]) continue;
314            Value rem = (*_capacity)[a] - (*_flow)[a];
315            if (!_tolerance.positive(rem)) continue;
316            if ((*_bucket)[v] == under_bucket) {
317              if (!(*_active)[v] && v != target) {
318                activate(v);
319              }
320              if (!_tolerance.less(rem, excess)) {
321                _flow->set(a, (*_flow)[a] + excess);
322                _excess->set(v, (*_excess)[v] + excess);
323                excess = 0;
324                goto no_more_push;
325              } else {
326                excess -= rem;
327                _excess->set(v, (*_excess)[v] + rem);
328                _flow->set(a, (*_capacity)[a]);
329              }
330            } else if (next_bucket > (*_bucket)[v]) {
331              next_bucket = (*_bucket)[v];
332            }
333          }
334
335          for (InArcIt a(_graph, n); a != INVALID; ++a) {
336            Node v = _graph.source(a);
337            if (_dormant[(*_bucket)[v]]) continue;
338            Value rem = (*_flow)[a];
339            if (!_tolerance.positive(rem)) continue;
340            if ((*_bucket)[v] == under_bucket) {
341              if (!(*_active)[v] && v != target) {
342                activate(v);
343              }
344              if (!_tolerance.less(rem, excess)) {
345                _flow->set(a, (*_flow)[a] - excess);
346                _excess->set(v, (*_excess)[v] + excess);
347                excess = 0;
348                goto no_more_push;
349              } else {
350                excess -= rem;
351                _excess->set(v, (*_excess)[v] + rem);
352                _flow->set(a, 0);
353              }
354            } else if (next_bucket > (*_bucket)[v]) {
355              next_bucket = (*_bucket)[v];
356            }
357          }
358
359        no_more_push:
360
361          _excess->set(n, excess);
362
363          if (excess != 0) {
364            if ((*_next)[n] == INVALID) {
365              typename std::list<std::list<int> >::iterator new_set =
366                _sets.insert(--_sets.end(), std::list<int>());
367              new_set->splice(new_set->end(), _sets.back(),
368                              _sets.back().begin(), ++_highest);
369              for (std::list<int>::iterator it = new_set->begin();
370                   it != new_set->end(); ++it) {
371                _dormant[*it] = true;
372              }
373              while (_highest != _sets.back().end() &&
374                     !(*_active)[_first[*_highest]]) {
375                ++_highest;
376              }
377            } else if (next_bucket == _node_num) {
378              _first[(*_bucket)[n]] = (*_next)[n];
379              _prev->set((*_next)[n], INVALID);
380
381              std::list<std::list<int> >::iterator new_set =
382                _sets.insert(--_sets.end(), std::list<int>());
383
384              new_set->push_front(bucket_num);
385              _bucket->set(n, bucket_num);
386              _first[bucket_num] = _last[bucket_num] = n;
387              _next->set(n, INVALID);
388              _prev->set(n, INVALID);
389              _dormant[bucket_num] = true;
390              ++bucket_num;
391
392              while (_highest != _sets.back().end() &&
393                     !(*_active)[_first[*_highest]]) {
394                ++_highest;
395              }
396            } else {
397              _first[*_highest] = (*_next)[n];
398              _prev->set((*_next)[n], INVALID);
399
400              while (next_bucket != *_highest) {
401                --_highest;
402              }
403
404              if (_highest == _sets.back().begin()) {
405                _sets.back().push_front(bucket_num);
406                _dormant[bucket_num] = false;
407                _first[bucket_num] = _last[bucket_num] = INVALID;
408                ++bucket_num;
409              }
410              --_highest;
411
412              _bucket->set(n, *_highest);
413              _next->set(n, _first[*_highest]);
414              if (_first[*_highest] != INVALID) {
415                _prev->set(_first[*_highest], n);
416              } else {
417                _last[*_highest] = n;
418              }
419              _first[*_highest] = n;
420            }
421          } else {
422
423            deactivate(n);
424            if (!(*_active)[_first[*_highest]]) {
425              ++_highest;
426              if (_highest != _sets.back().end() &&
427                  !(*_active)[_first[*_highest]]) {
428                _highest = _sets.back().end();
429              }
430            }
431          }
432        }
433
434        if ((*_excess)[target] < _min_cut) {
435          _min_cut = (*_excess)[target];
436          for (NodeIt i(_graph); i != INVALID; ++i) {
437            _min_cut_map->set(i, true);
438          }
439          for (std::list<int>::iterator it = _sets.back().begin();
440               it != _sets.back().end(); ++it) {
441            Node n = _first[*it];
442            while (n != INVALID) {
443              _min_cut_map->set(n, false);
444              n = (*_next)[n];
445            }
446          }
447        }
448
449        {
450          Node new_target;
451          if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
452            if ((*_next)[target] == INVALID) {
453              _last[(*_bucket)[target]] = (*_prev)[target];
454              new_target = (*_prev)[target];
455            } else {
456              _prev->set((*_next)[target], (*_prev)[target]);
457              new_target = (*_next)[target];
458            }
459            if ((*_prev)[target] == INVALID) {
460              _first[(*_bucket)[target]] = (*_next)[target];
461            } else {
462              _next->set((*_prev)[target], (*_next)[target]);
463            }
464          } else {
465            _sets.back().pop_back();
466            if (_sets.back().empty()) {
467              _sets.pop_back();
468              if (_sets.empty())
469                break;
470              for (std::list<int>::iterator it = _sets.back().begin();
471                   it != _sets.back().end(); ++it) {
472                _dormant[*it] = false;
473              }
474            }
475            new_target = _last[_sets.back().back()];
476          }
477
478          _bucket->set(target, 0);
479
480          _source_set->set(target, true);
481          for (OutArcIt a(_graph, target); a != INVALID; ++a) {
482            Value rem = (*_capacity)[a] - (*_flow)[a];
483            if (!_tolerance.positive(rem)) continue;
484            Node v = _graph.target(a);
485            if (!(*_active)[v] && !(*_source_set)[v]) {
486              activate(v);
487            }
488            _excess->set(v, (*_excess)[v] + rem);
489            _flow->set(a, (*_capacity)[a]);
490          }
491
492          for (InArcIt a(_graph, target); a != INVALID; ++a) {
493            Value rem = (*_flow)[a];
494            if (!_tolerance.positive(rem)) continue;
495            Node v = _graph.source(a);
496            if (!(*_active)[v] && !(*_source_set)[v]) {
497              activate(v);
498            }
499            _excess->set(v, (*_excess)[v] + rem);
500            _flow->set(a, 0);
501          }
502
503          target = new_target;
504          if ((*_active)[target]) {
505            deactivate(target);
506          }
507
508          _highest = _sets.back().begin();
509          while (_highest != _sets.back().end() &&
510                 !(*_active)[_first[*_highest]]) {
511            ++_highest;
512          }
513        }
514      }
515    }
516
517    void findMinCutIn() {
518
519      for (NodeIt n(_graph); n != INVALID; ++n) {
520        _excess->set(n, 0);
521      }
522
523      for (ArcIt a(_graph); a != INVALID; ++a) {
524        _flow->set(a, 0);
525      }
526
527      int bucket_num = 0;
528      std::vector<Node> queue(_node_num);
529      int qfirst = 0, qlast = 0, qsep = 0;
530
531      {
532        typename Digraph::template NodeMap<bool> reached(_graph, false);
533
534        reached.set(_source, true);
535
536        bool first_set = true;
537
538        for (NodeIt t(_graph); t != INVALID; ++t) {
539          if (reached[t]) continue;
540          _sets.push_front(std::list<int>());
541         
542          queue[qlast++] = t;
543          reached.set(t, true);
544
545          while (qfirst != qlast) {
546            if (qsep == qfirst) {
547              ++bucket_num;
548              _sets.front().push_front(bucket_num);
549              _dormant[bucket_num] = !first_set;
550              _first[bucket_num] = _last[bucket_num] = INVALID;
551              qsep = qlast;
552            }
553
554            Node n = queue[qfirst++];
555            addItem(n, bucket_num);
556
557            for (OutArcIt a(_graph, n); a != INVALID; ++a) {
558              Node u = _graph.target(a);
559              if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
560                reached.set(u, true);
561                queue[qlast++] = u;
562              }
563            }
564          }
565          first_set = false;
566        }
567
568        ++bucket_num;
569        _bucket->set(_source, 0);
570        _dormant[0] = true;
571      }
572      _source_set->set(_source, true);
573
574      Node target = _last[_sets.back().back()];
575      {
576        for (InArcIt a(_graph, _source); a != INVALID; ++a) {
577          if (_tolerance.positive((*_capacity)[a])) {
578            Node u = _graph.source(a);
579            _flow->set(a, (*_capacity)[a]);
580            _excess->set(u, (*_excess)[u] + (*_capacity)[a]);
581            if (!(*_active)[u] && u != _source) {
582              activate(u);
583            }
584          }
585        }
586        if ((*_active)[target]) {
587          deactivate(target);
588        }
589
590        _highest = _sets.back().begin();
591        while (_highest != _sets.back().end() &&
592               !(*_active)[_first[*_highest]]) {
593          ++_highest;
594        }
595      }
596
597
598      while (true) {
599        while (_highest != _sets.back().end()) {
600          Node n = _first[*_highest];
601          Value excess = (*_excess)[n];
602          int next_bucket = _node_num;
603
604          int under_bucket;
605          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
606            under_bucket = -1;
607          } else {
608            under_bucket = *(++std::list<int>::iterator(_highest));
609          }
610
611          for (InArcIt a(_graph, n); a != INVALID; ++a) {
612            Node v = _graph.source(a);
613            if (_dormant[(*_bucket)[v]]) continue;
614            Value rem = (*_capacity)[a] - (*_flow)[a];
615            if (!_tolerance.positive(rem)) continue;
616            if ((*_bucket)[v] == under_bucket) {
617              if (!(*_active)[v] && v != target) {
618                activate(v);
619              }
620              if (!_tolerance.less(rem, excess)) {
621                _flow->set(a, (*_flow)[a] + excess);
622                _excess->set(v, (*_excess)[v] + excess);
623                excess = 0;
624                goto no_more_push;
625              } else {
626                excess -= rem;
627                _excess->set(v, (*_excess)[v] + rem);
628                _flow->set(a, (*_capacity)[a]);
629              }
630            } else if (next_bucket > (*_bucket)[v]) {
631              next_bucket = (*_bucket)[v];
632            }
633          }
634
635          for (OutArcIt a(_graph, n); a != INVALID; ++a) {
636            Node v = _graph.target(a);
637            if (_dormant[(*_bucket)[v]]) continue;
638            Value rem = (*_flow)[a];
639            if (!_tolerance.positive(rem)) continue;
640            if ((*_bucket)[v] == under_bucket) {
641              if (!(*_active)[v] && v != target) {
642                activate(v);
643              }
644              if (!_tolerance.less(rem, excess)) {
645                _flow->set(a, (*_flow)[a] - excess);
646                _excess->set(v, (*_excess)[v] + excess);
647                excess = 0;
648                goto no_more_push;
649              } else {
650                excess -= rem;
651                _excess->set(v, (*_excess)[v] + rem);
652                _flow->set(a, 0);
653              }
654            } else if (next_bucket > (*_bucket)[v]) {
655              next_bucket = (*_bucket)[v];
656            }
657          }
658
659        no_more_push:
660
661          _excess->set(n, excess);
662
663          if (excess != 0) {
664            if ((*_next)[n] == INVALID) {
665              typename std::list<std::list<int> >::iterator new_set =
666                _sets.insert(--_sets.end(), std::list<int>());
667              new_set->splice(new_set->end(), _sets.back(),
668                              _sets.back().begin(), ++_highest);
669              for (std::list<int>::iterator it = new_set->begin();
670                   it != new_set->end(); ++it) {
671                _dormant[*it] = true;
672              }
673              while (_highest != _sets.back().end() &&
674                     !(*_active)[_first[*_highest]]) {
675                ++_highest;
676              }
677            } else if (next_bucket == _node_num) {
678              _first[(*_bucket)[n]] = (*_next)[n];
679              _prev->set((*_next)[n], INVALID);
680
681              std::list<std::list<int> >::iterator new_set =
682                _sets.insert(--_sets.end(), std::list<int>());
683
684              new_set->push_front(bucket_num);
685              _bucket->set(n, bucket_num);
686              _first[bucket_num] = _last[bucket_num] = n;
687              _next->set(n, INVALID);
688              _prev->set(n, INVALID);
689              _dormant[bucket_num] = true;
690              ++bucket_num;
691
692              while (_highest != _sets.back().end() &&
693                     !(*_active)[_first[*_highest]]) {
694                ++_highest;
695              }
696            } else {
697              _first[*_highest] = (*_next)[n];
698              _prev->set((*_next)[n], INVALID);
699
700              while (next_bucket != *_highest) {
701                --_highest;
702              }
703              if (_highest == _sets.back().begin()) {
704                _sets.back().push_front(bucket_num);
705                _dormant[bucket_num] = false;
706                _first[bucket_num] = _last[bucket_num] = INVALID;
707                ++bucket_num;
708              }
709              --_highest;
710
711              _bucket->set(n, *_highest);
712              _next->set(n, _first[*_highest]);
713              if (_first[*_highest] != INVALID) {
714                _prev->set(_first[*_highest], n);
715              } else {
716                _last[*_highest] = n;
717              }
718              _first[*_highest] = n;
719            }
720          } else {
721
722            deactivate(n);
723            if (!(*_active)[_first[*_highest]]) {
724              ++_highest;
725              if (_highest != _sets.back().end() &&
726                  !(*_active)[_first[*_highest]]) {
727                _highest = _sets.back().end();
728              }
729            }
730          }
731        }
732
733        if ((*_excess)[target] < _min_cut) {
734          _min_cut = (*_excess)[target];
735          for (NodeIt i(_graph); i != INVALID; ++i) {
736            _min_cut_map->set(i, false);
737          }
738          for (std::list<int>::iterator it = _sets.back().begin();
739               it != _sets.back().end(); ++it) {
740            Node n = _first[*it];
741            while (n != INVALID) {
742              _min_cut_map->set(n, true);
743              n = (*_next)[n];
744            }
745          }
746        }
747
748        {
749          Node new_target;
750          if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
751            if ((*_next)[target] == INVALID) {
752              _last[(*_bucket)[target]] = (*_prev)[target];
753              new_target = (*_prev)[target];
754            } else {
755              _prev->set((*_next)[target], (*_prev)[target]);
756              new_target = (*_next)[target];
757            }
758            if ((*_prev)[target] == INVALID) {
759              _first[(*_bucket)[target]] = (*_next)[target];
760            } else {
761              _next->set((*_prev)[target], (*_next)[target]);
762            }
763          } else {
764            _sets.back().pop_back();
765            if (_sets.back().empty()) {
766              _sets.pop_back();
767              if (_sets.empty())
768                break;
769              for (std::list<int>::iterator it = _sets.back().begin();
770                   it != _sets.back().end(); ++it) {
771                _dormant[*it] = false;
772              }
773            }
774            new_target = _last[_sets.back().back()];
775          }
776
777          _bucket->set(target, 0);
778
779          _source_set->set(target, true);
780          for (InArcIt a(_graph, target); a != INVALID; ++a) {
781            Value rem = (*_capacity)[a] - (*_flow)[a];
782            if (!_tolerance.positive(rem)) continue;
783            Node v = _graph.source(a);
784            if (!(*_active)[v] && !(*_source_set)[v]) {
785              activate(v);
786            }
787            _excess->set(v, (*_excess)[v] + rem);
788            _flow->set(a, (*_capacity)[a]);
789          }
790
791          for (OutArcIt a(_graph, target); a != INVALID; ++a) {
792            Value rem = (*_flow)[a];
793            if (!_tolerance.positive(rem)) continue;
794            Node v = _graph.target(a);
795            if (!(*_active)[v] && !(*_source_set)[v]) {
796              activate(v);
797            }
798            _excess->set(v, (*_excess)[v] + rem);
799            _flow->set(a, 0);
800          }
801
802          target = new_target;
803          if ((*_active)[target]) {
804            deactivate(target);
805          }
806
807          _highest = _sets.back().begin();
808          while (_highest != _sets.back().end() &&
809                 !(*_active)[_first[*_highest]]) {
810            ++_highest;
811          }
812        }
813      }
814    }
815
816  public:
817
818    /// \name Execution control
819    /// The simplest way to execute the algorithm is to use
820    /// one of the member functions called \c run(...).
821    /// \n
822    /// If you need more control on the execution,
823    /// first you must call \ref init(), then the \ref calculateIn() or
824    /// \ref calculateOut() functions.
825
826    /// @{
827
828    /// \brief Initializes the internal data structures.
829    ///
830    /// Initializes the internal data structures. It creates
831    /// the maps, residual graph adaptors and some bucket structures
832    /// for the algorithm.
833    void init() {
834      init(NodeIt(_graph));
835    }
836
837    /// \brief Initializes the internal data structures.
838    ///
839    /// Initializes the internal data structures. It creates
840    /// the maps, residual graph adaptor and some bucket structures
841    /// for the algorithm. Node \c source  is used as the push-relabel
842    /// algorithm's source.
843    void init(const Node& source) {
844      _source = source;
845
846      _node_num = countNodes(_graph);
847
848      _first.resize(_node_num);
849      _last.resize(_node_num);
850
851      _dormant.resize(_node_num);
852
853      if (!_flow) {
854        _flow = new FlowMap(_graph);
855      }
856      if (!_next) {
857        _next = new typename Digraph::template NodeMap<Node>(_graph);
858      }
859      if (!_prev) {
860        _prev = new typename Digraph::template NodeMap<Node>(_graph);
861      }
862      if (!_active) {
863        _active = new typename Digraph::template NodeMap<bool>(_graph);
864      }
865      if (!_bucket) {
866        _bucket = new typename Digraph::template NodeMap<int>(_graph);
867      }
868      if (!_excess) {
869        _excess = new ExcessMap(_graph);
870      }
871      if (!_source_set) {
872        _source_set = new SourceSetMap(_graph);
873      }
874      if (!_min_cut_map) {
875        _min_cut_map = new MinCutMap(_graph);
876      }
877
878      _min_cut = std::numeric_limits<Value>::max();
879    }
880
881
882    /// \brief Calculates a minimum cut with \f$ source \f$ on the
883    /// source-side.
884    ///
885    /// Calculates a minimum cut with \f$ source \f$ on the
886    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with
887    /// \f$ source \in X \f$ and minimal out-degree).
888    void calculateOut() {
889      findMinCutOut();
890    }
891
892    /// \brief Calculates a minimum cut with \f$ source \f$ on the
893    /// target-side.
894    ///
895    /// Calculates a minimum cut with \f$ source \f$ on the
896    /// target-side (i.e. a set \f$ X\subsetneq V \f$ with
897    /// \f$ source \in X \f$ and minimal out-degree).
898    void calculateIn() {
899      findMinCutIn();
900    }
901
902
903    /// \brief Runs the algorithm.
904    ///
905    /// Runs the algorithm. It finds nodes \c source and \c target
906    /// arbitrarily and then calls \ref init(), \ref calculateOut()
907    /// and \ref calculateIn().
908    void run() {
909      init();
910      calculateOut();
911      calculateIn();
912    }
913
914    /// \brief Runs the algorithm.
915    ///
916    /// Runs the algorithm. It uses the given \c source node, finds a
917    /// proper \c target and then calls the \ref init(), \ref
918    /// calculateOut() and \ref calculateIn().
919    void run(const Node& s) {
920      init(s);
921      calculateOut();
922      calculateIn();
923    }
924
925    /// @}
926
927    /// \name Query Functions
928    /// The result of the %HaoOrlin algorithm
929    /// can be obtained using these functions.
930    /// \n
931    /// Before using these functions, either \ref run(), \ref
932    /// calculateOut() or \ref calculateIn() must be called.
933
934    /// @{
935
936    /// \brief Returns the value of the minimum value cut.
937    ///
938    /// Returns the value of the minimum value cut.
939    Value minCutValue() const {
940      return _min_cut;
941    }
942
943
944    /// \brief Returns a minimum cut.
945    ///
946    /// Sets \c nodeMap to the characteristic vector of a minimum
947    /// value cut: it will give a nonempty set \f$ X\subsetneq V \f$
948    /// with minimal out-degree (i.e. \c nodeMap will be true exactly
949    /// for the nodes of \f$ X \f$).  \pre nodeMap should be a
950    /// bool-valued node-map.
951    template <typename NodeMap>
952    Value minCutMap(NodeMap& nodeMap) const {
953      for (NodeIt it(_graph); it != INVALID; ++it) {
954        nodeMap.set(it, (*_min_cut_map)[it]);
955      }
956      return _min_cut;
957    }
958
959    /// @}
960
961  }; //class HaoOrlin
962
963
964} //namespace lemon
965
966#endif //LEMON_HAO_ORLIN_H
Note: See TracBrowser for help on using the repository browser.