COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/hao_orlin.h

Last change on this file was 2624:dc4dd5fc0e25, checked in by Balazs Dezso, 15 years ago

Bug fixes is HaoOrlin? and MinCostArborescence?

MinCostArborescence?

  • proper deallocation

HaoOrlin?

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