COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/hao_orlin.h @ 2530:f86f7e4eb2ba

Last change on this file since 2530:f86f7e4eb2ba was 2530:f86f7e4eb2ba, checked in by Balazs Dezso, 11 years ago

Reimplementation of Hao-Orlin algorithm
Little modifictaion in NagamochiIbaraki?
More docs for minimum cut algorithms

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