# HG changeset patch # User Balazs Dezso # Date 2008-12-01 23:12:16 # Node ID b8ce15103485867c017d78e50739fac04e1c9b4a # Parent 59d3aa4f921f4ec4b4c4e24e0edaae4ca7038c5b Port Hao-Orlin algorithm from SVN -r3509 (#58) diff --git a/lemon/Makefile.am b/lemon/Makefile.am --- a/lemon/Makefile.am +++ b/lemon/Makefile.am @@ -36,6 +36,7 @@ lemon/grid_graph.h \ lemon/hypercube_graph.h \ lemon/kruskal.h \ + lemon/hao_orlin.h \ lemon/lgf_reader.h \ lemon/lgf_writer.h \ lemon/list_graph.h \ diff --git a/lemon/hao_orlin.h b/lemon/hao_orlin.h new file mode 100644 --- /dev/null +++ b/lemon/hao_orlin.h @@ -0,0 +1,985 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2008 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#ifndef LEMON_HAO_ORLIN_H +#define LEMON_HAO_ORLIN_H + +#include +#include +#include + +#include +#include +#include + +/// \file +/// \ingroup min_cut +/// \brief Implementation of the Hao-Orlin algorithm. +/// +/// Implementation of the Hao-Orlin algorithm class for testing network +/// reliability. + +namespace lemon { + + /// \ingroup min_cut + /// + /// \brief %Hao-Orlin algorithm to find a minimum cut in directed graphs. + /// + /// Hao-Orlin calculates a minimum cut in a directed graph + /// \f$D=(V,A)\f$. It takes a fixed node \f$ source \in V \f$ and + /// consists of two phases: in the first phase it determines a + /// minimum cut with \f$ source \f$ on the source-side (i.e. a set + /// \f$ X\subsetneq V \f$ with \f$ source \in X \f$ and minimal + /// out-degree) and in the second phase it determines a minimum cut + /// with \f$ source \f$ on the sink-side (i.e. a set + /// \f$ X\subsetneq V \f$ with \f$ source \notin X \f$ and minimal + /// out-degree). Obviously, the smaller of these two cuts will be a + /// minimum cut of \f$ D \f$. The algorithm is a modified + /// push-relabel preflow algorithm and our implementation calculates + /// the minimum cut in \f$ O(n^2\sqrt{m}) \f$ time (we use the + /// highest-label rule), or in \f$O(nm)\f$ for unit capacities. The + /// purpose of such algorithm is testing network reliability. For an + /// undirected graph you can run just the first phase of the + /// algorithm or you can use the algorithm of Nagamochi and Ibaraki + /// which solves the undirected problem in + /// \f$ O(nm + n^2 \log(n)) \f$ time: it is implemented in the + /// NagamochiIbaraki algorithm class. + /// + /// \param _Digraph is the graph type of the algorithm. + /// \param _CapacityMap is an edge map of capacities which should + /// be any numreric type. The default type is _Digraph::ArcMap. + /// \param _Tolerance is the handler of the inexact computation. The + /// default type for this is Tolerance. +#ifdef DOXYGEN + template +#else + template , + typename _Tolerance = Tolerance > +#endif + class HaoOrlin { + private: + + typedef _Digraph Digraph; + typedef _CapacityMap CapacityMap; + typedef _Tolerance Tolerance; + + typedef typename CapacityMap::Value Value; + + TEMPLATE_GRAPH_TYPEDEFS(Digraph); + + const Digraph& _graph; + const CapacityMap* _capacity; + + typedef typename Digraph::template ArcMap FlowMap; + FlowMap* _flow; + + Node _source; + + int _node_num; + + // Bucketing structure + std::vector _first, _last; + typename Digraph::template NodeMap* _next; + typename Digraph::template NodeMap* _prev; + typename Digraph::template NodeMap* _active; + typename Digraph::template NodeMap* _bucket; + + std::vector _dormant; + + std::list > _sets; + std::list::iterator _highest; + + typedef typename Digraph::template NodeMap ExcessMap; + ExcessMap* _excess; + + typedef typename Digraph::template NodeMap SourceSetMap; + SourceSetMap* _source_set; + + Value _min_cut; + + typedef typename Digraph::template NodeMap MinCutMap; + MinCutMap* _min_cut_map; + + Tolerance _tolerance; + + public: + + /// \brief Constructor + /// + /// Constructor of the algorithm class. + HaoOrlin(const Digraph& graph, const CapacityMap& capacity, + const Tolerance& tolerance = Tolerance()) : + _graph(graph), _capacity(&capacity), _flow(0), _source(), + _node_num(), _first(), _last(), _next(0), _prev(0), + _active(0), _bucket(0), _dormant(), _sets(), _highest(), + _excess(0), _source_set(0), _min_cut(), _min_cut_map(0), + _tolerance(tolerance) {} + + ~HaoOrlin() { + if (_min_cut_map) { + delete _min_cut_map; + } + if (_source_set) { + delete _source_set; + } + if (_excess) { + delete _excess; + } + if (_next) { + delete _next; + } + if (_prev) { + delete _prev; + } + if (_active) { + delete _active; + } + if (_bucket) { + delete _bucket; + } + if (_flow) { + delete _flow; + } + } + + private: + + void activate(const Node& i) { + _active->set(i, true); + + int bucket = (*_bucket)[i]; + + if ((*_prev)[i] == INVALID || (*_active)[(*_prev)[i]]) return; + //unlace + _next->set((*_prev)[i], (*_next)[i]); + if ((*_next)[i] != INVALID) { + _prev->set((*_next)[i], (*_prev)[i]); + } else { + _last[bucket] = (*_prev)[i]; + } + //lace + _next->set(i, _first[bucket]); + _prev->set(_first[bucket], i); + _prev->set(i, INVALID); + _first[bucket] = i; + } + + void deactivate(const Node& i) { + _active->set(i, false); + int bucket = (*_bucket)[i]; + + if ((*_next)[i] == INVALID || !(*_active)[(*_next)[i]]) return; + + //unlace + _prev->set((*_next)[i], (*_prev)[i]); + if ((*_prev)[i] != INVALID) { + _next->set((*_prev)[i], (*_next)[i]); + } else { + _first[bucket] = (*_next)[i]; + } + //lace + _prev->set(i, _last[bucket]); + _next->set(_last[bucket], i); + _next->set(i, INVALID); + _last[bucket] = i; + } + + void addItem(const Node& i, int bucket) { + (*_bucket)[i] = bucket; + if (_last[bucket] != INVALID) { + _prev->set(i, _last[bucket]); + _next->set(_last[bucket], i); + _next->set(i, INVALID); + _last[bucket] = i; + } else { + _prev->set(i, INVALID); + _first[bucket] = i; + _next->set(i, INVALID); + _last[bucket] = i; + } + } + + void findMinCutOut() { + + for (NodeIt n(_graph); n != INVALID; ++n) { + _excess->set(n, 0); + } + + for (ArcIt a(_graph); a != INVALID; ++a) { + _flow->set(a, 0); + } + + int bucket_num = 1; + + { + typename Digraph::template NodeMap reached(_graph, false); + + reached.set(_source, true); + + bool first_set = true; + + for (NodeIt t(_graph); t != INVALID; ++t) { + if (reached[t]) continue; + _sets.push_front(std::list()); + _sets.front().push_front(bucket_num); + _dormant[bucket_num] = !first_set; + + _bucket->set(t, bucket_num); + _first[bucket_num] = _last[bucket_num] = t; + _next->set(t, INVALID); + _prev->set(t, INVALID); + + ++bucket_num; + + std::vector queue; + queue.push_back(t); + reached.set(t, true); + + while (!queue.empty()) { + _sets.front().push_front(bucket_num); + _dormant[bucket_num] = !first_set; + _first[bucket_num] = _last[bucket_num] = INVALID; + + std::vector nqueue; + for (int i = 0; i < int(queue.size()); ++i) { + Node n = queue[i]; + for (InArcIt a(_graph, n); a != INVALID; ++a) { + Node u = _graph.source(a); + if (!reached[u] && _tolerance.positive((*_capacity)[a])) { + reached.set(u, true); + addItem(u, bucket_num); + nqueue.push_back(u); + } + } + } + queue.swap(nqueue); + ++bucket_num; + } + _sets.front().pop_front(); + --bucket_num; + first_set = false; + } + + _bucket->set(_source, 0); + _dormant[0] = true; + } + _source_set->set(_source, true); + + Node target = _last[_sets.back().back()]; + { + for (OutArcIt a(_graph, _source); a != INVALID; ++a) { + if (_tolerance.positive((*_capacity)[a])) { + Node u = _graph.target(a); + _flow->set(a, (*_capacity)[a]); + _excess->set(u, (*_excess)[u] + (*_capacity)[a]); + if (!(*_active)[u] && u != _source) { + activate(u); + } + } + } + + if ((*_active)[target]) { + deactivate(target); + } + + _highest = _sets.back().begin(); + while (_highest != _sets.back().end() && + !(*_active)[_first[*_highest]]) { + ++_highest; + } + } + + while (true) { + while (_highest != _sets.back().end()) { + Node n = _first[*_highest]; + Value excess = (*_excess)[n]; + int next_bucket = _node_num; + + int under_bucket; + if (++std::list::iterator(_highest) == _sets.back().end()) { + under_bucket = -1; + } else { + under_bucket = *(++std::list::iterator(_highest)); + } + + for (OutArcIt a(_graph, n); a != INVALID; ++a) { + Node v = _graph.target(a); + if (_dormant[(*_bucket)[v]]) continue; + Value rem = (*_capacity)[a] - (*_flow)[a]; + if (!_tolerance.positive(rem)) continue; + if ((*_bucket)[v] == under_bucket) { + if (!(*_active)[v] && v != target) { + activate(v); + } + if (!_tolerance.less(rem, excess)) { + _flow->set(a, (*_flow)[a] + excess); + _excess->set(v, (*_excess)[v] + excess); + excess = 0; + goto no_more_push; + } else { + excess -= rem; + _excess->set(v, (*_excess)[v] + rem); + _flow->set(a, (*_capacity)[a]); + } + } else if (next_bucket > (*_bucket)[v]) { + next_bucket = (*_bucket)[v]; + } + } + + for (InArcIt a(_graph, n); a != INVALID; ++a) { + Node v = _graph.source(a); + if (_dormant[(*_bucket)[v]]) continue; + Value rem = (*_flow)[a]; + if (!_tolerance.positive(rem)) continue; + if ((*_bucket)[v] == under_bucket) { + if (!(*_active)[v] && v != target) { + activate(v); + } + if (!_tolerance.less(rem, excess)) { + _flow->set(a, (*_flow)[a] - excess); + _excess->set(v, (*_excess)[v] + excess); + excess = 0; + goto no_more_push; + } else { + excess -= rem; + _excess->set(v, (*_excess)[v] + rem); + _flow->set(a, 0); + } + } else if (next_bucket > (*_bucket)[v]) { + next_bucket = (*_bucket)[v]; + } + } + + no_more_push: + + _excess->set(n, excess); + + if (excess != 0) { + if ((*_next)[n] == INVALID) { + typename std::list >::iterator new_set = + _sets.insert(--_sets.end(), std::list()); + new_set->splice(new_set->end(), _sets.back(), + _sets.back().begin(), ++_highest); + for (std::list::iterator it = new_set->begin(); + it != new_set->end(); ++it) { + _dormant[*it] = true; + } + while (_highest != _sets.back().end() && + !(*_active)[_first[*_highest]]) { + ++_highest; + } + } else if (next_bucket == _node_num) { + _first[(*_bucket)[n]] = (*_next)[n]; + _prev->set((*_next)[n], INVALID); + + std::list >::iterator new_set = + _sets.insert(--_sets.end(), std::list()); + + new_set->push_front(bucket_num); + _bucket->set(n, bucket_num); + _first[bucket_num] = _last[bucket_num] = n; + _next->set(n, INVALID); + _prev->set(n, INVALID); + _dormant[bucket_num] = true; + ++bucket_num; + + while (_highest != _sets.back().end() && + !(*_active)[_first[*_highest]]) { + ++_highest; + } + } else { + _first[*_highest] = (*_next)[n]; + _prev->set((*_next)[n], INVALID); + + while (next_bucket != *_highest) { + --_highest; + } + + if (_highest == _sets.back().begin()) { + _sets.back().push_front(bucket_num); + _dormant[bucket_num] = false; + _first[bucket_num] = _last[bucket_num] = INVALID; + ++bucket_num; + } + --_highest; + + _bucket->set(n, *_highest); + _next->set(n, _first[*_highest]); + if (_first[*_highest] != INVALID) { + _prev->set(_first[*_highest], n); + } else { + _last[*_highest] = n; + } + _first[*_highest] = n; + } + } else { + + deactivate(n); + if (!(*_active)[_first[*_highest]]) { + ++_highest; + if (_highest != _sets.back().end() && + !(*_active)[_first[*_highest]]) { + _highest = _sets.back().end(); + } + } + } + } + + if ((*_excess)[target] < _min_cut) { + _min_cut = (*_excess)[target]; + for (NodeIt i(_graph); i != INVALID; ++i) { + _min_cut_map->set(i, true); + } + for (std::list::iterator it = _sets.back().begin(); + it != _sets.back().end(); ++it) { + Node n = _first[*it]; + while (n != INVALID) { + _min_cut_map->set(n, false); + n = (*_next)[n]; + } + } + } + + { + Node new_target; + if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) { + if ((*_next)[target] == INVALID) { + _last[(*_bucket)[target]] = (*_prev)[target]; + new_target = (*_prev)[target]; + } else { + _prev->set((*_next)[target], (*_prev)[target]); + new_target = (*_next)[target]; + } + if ((*_prev)[target] == INVALID) { + _first[(*_bucket)[target]] = (*_next)[target]; + } else { + _next->set((*_prev)[target], (*_next)[target]); + } + } else { + _sets.back().pop_back(); + if (_sets.back().empty()) { + _sets.pop_back(); + if (_sets.empty()) + break; + for (std::list::iterator it = _sets.back().begin(); + it != _sets.back().end(); ++it) { + _dormant[*it] = false; + } + } + new_target = _last[_sets.back().back()]; + } + + _bucket->set(target, 0); + + _source_set->set(target, true); + for (OutArcIt a(_graph, target); a != INVALID; ++a) { + Value rem = (*_capacity)[a] - (*_flow)[a]; + if (!_tolerance.positive(rem)) continue; + Node v = _graph.target(a); + if (!(*_active)[v] && !(*_source_set)[v]) { + activate(v); + } + _excess->set(v, (*_excess)[v] + rem); + _flow->set(a, (*_capacity)[a]); + } + + for (InArcIt a(_graph, target); a != INVALID; ++a) { + Value rem = (*_flow)[a]; + if (!_tolerance.positive(rem)) continue; + Node v = _graph.source(a); + if (!(*_active)[v] && !(*_source_set)[v]) { + activate(v); + } + _excess->set(v, (*_excess)[v] + rem); + _flow->set(a, 0); + } + + target = new_target; + if ((*_active)[target]) { + deactivate(target); + } + + _highest = _sets.back().begin(); + while (_highest != _sets.back().end() && + !(*_active)[_first[*_highest]]) { + ++_highest; + } + } + } + } + + void findMinCutIn() { + + for (NodeIt n(_graph); n != INVALID; ++n) { + _excess->set(n, 0); + } + + for (ArcIt a(_graph); a != INVALID; ++a) { + _flow->set(a, 0); + } + + int bucket_num = 1; + + { + typename Digraph::template NodeMap reached(_graph, false); + + reached.set(_source, true); + + bool first_set = true; + + for (NodeIt t(_graph); t != INVALID; ++t) { + if (reached[t]) continue; + _sets.push_front(std::list()); + _sets.front().push_front(bucket_num); + _dormant[bucket_num] = !first_set; + + _bucket->set(t, bucket_num); + _first[bucket_num] = _last[bucket_num] = t; + _next->set(t, INVALID); + _prev->set(t, INVALID); + + ++bucket_num; + + std::vector queue; + queue.push_back(t); + reached.set(t, true); + + while (!queue.empty()) { + _sets.front().push_front(bucket_num); + _dormant[bucket_num] = !first_set; + _first[bucket_num] = _last[bucket_num] = INVALID; + + std::vector nqueue; + for (int i = 0; i < int(queue.size()); ++i) { + Node n = queue[i]; + for (OutArcIt a(_graph, n); a != INVALID; ++a) { + Node u = _graph.target(a); + if (!reached[u] && _tolerance.positive((*_capacity)[a])) { + reached.set(u, true); + addItem(u, bucket_num); + nqueue.push_back(u); + } + } + } + queue.swap(nqueue); + ++bucket_num; + } + _sets.front().pop_front(); + --bucket_num; + first_set = false; + } + + _bucket->set(_source, 0); + _dormant[0] = true; + } + _source_set->set(_source, true); + + Node target = _last[_sets.back().back()]; + { + for (InArcIt a(_graph, _source); a != INVALID; ++a) { + if (_tolerance.positive((*_capacity)[a])) { + Node u = _graph.source(a); + _flow->set(a, (*_capacity)[a]); + _excess->set(u, (*_excess)[u] + (*_capacity)[a]); + if (!(*_active)[u] && u != _source) { + activate(u); + } + } + } + if ((*_active)[target]) { + deactivate(target); + } + + _highest = _sets.back().begin(); + while (_highest != _sets.back().end() && + !(*_active)[_first[*_highest]]) { + ++_highest; + } + } + + + while (true) { + while (_highest != _sets.back().end()) { + Node n = _first[*_highest]; + Value excess = (*_excess)[n]; + int next_bucket = _node_num; + + int under_bucket; + if (++std::list::iterator(_highest) == _sets.back().end()) { + under_bucket = -1; + } else { + under_bucket = *(++std::list::iterator(_highest)); + } + + for (InArcIt a(_graph, n); a != INVALID; ++a) { + Node v = _graph.source(a); + if (_dormant[(*_bucket)[v]]) continue; + Value rem = (*_capacity)[a] - (*_flow)[a]; + if (!_tolerance.positive(rem)) continue; + if ((*_bucket)[v] == under_bucket) { + if (!(*_active)[v] && v != target) { + activate(v); + } + if (!_tolerance.less(rem, excess)) { + _flow->set(a, (*_flow)[a] + excess); + _excess->set(v, (*_excess)[v] + excess); + excess = 0; + goto no_more_push; + } else { + excess -= rem; + _excess->set(v, (*_excess)[v] + rem); + _flow->set(a, (*_capacity)[a]); + } + } else if (next_bucket > (*_bucket)[v]) { + next_bucket = (*_bucket)[v]; + } + } + + for (OutArcIt a(_graph, n); a != INVALID; ++a) { + Node v = _graph.target(a); + if (_dormant[(*_bucket)[v]]) continue; + Value rem = (*_flow)[a]; + if (!_tolerance.positive(rem)) continue; + if ((*_bucket)[v] == under_bucket) { + if (!(*_active)[v] && v != target) { + activate(v); + } + if (!_tolerance.less(rem, excess)) { + _flow->set(a, (*_flow)[a] - excess); + _excess->set(v, (*_excess)[v] + excess); + excess = 0; + goto no_more_push; + } else { + excess -= rem; + _excess->set(v, (*_excess)[v] + rem); + _flow->set(a, 0); + } + } else if (next_bucket > (*_bucket)[v]) { + next_bucket = (*_bucket)[v]; + } + } + + no_more_push: + + _excess->set(n, excess); + + if (excess != 0) { + if ((*_next)[n] == INVALID) { + typename std::list >::iterator new_set = + _sets.insert(--_sets.end(), std::list()); + new_set->splice(new_set->end(), _sets.back(), + _sets.back().begin(), ++_highest); + for (std::list::iterator it = new_set->begin(); + it != new_set->end(); ++it) { + _dormant[*it] = true; + } + while (_highest != _sets.back().end() && + !(*_active)[_first[*_highest]]) { + ++_highest; + } + } else if (next_bucket == _node_num) { + _first[(*_bucket)[n]] = (*_next)[n]; + _prev->set((*_next)[n], INVALID); + + std::list >::iterator new_set = + _sets.insert(--_sets.end(), std::list()); + + new_set->push_front(bucket_num); + _bucket->set(n, bucket_num); + _first[bucket_num] = _last[bucket_num] = n; + _next->set(n, INVALID); + _prev->set(n, INVALID); + _dormant[bucket_num] = true; + ++bucket_num; + + while (_highest != _sets.back().end() && + !(*_active)[_first[*_highest]]) { + ++_highest; + } + } else { + _first[*_highest] = (*_next)[n]; + _prev->set((*_next)[n], INVALID); + + while (next_bucket != *_highest) { + --_highest; + } + if (_highest == _sets.back().begin()) { + _sets.back().push_front(bucket_num); + _dormant[bucket_num] = false; + _first[bucket_num] = _last[bucket_num] = INVALID; + ++bucket_num; + } + --_highest; + + _bucket->set(n, *_highest); + _next->set(n, _first[*_highest]); + if (_first[*_highest] != INVALID) { + _prev->set(_first[*_highest], n); + } else { + _last[*_highest] = n; + } + _first[*_highest] = n; + } + } else { + + deactivate(n); + if (!(*_active)[_first[*_highest]]) { + ++_highest; + if (_highest != _sets.back().end() && + !(*_active)[_first[*_highest]]) { + _highest = _sets.back().end(); + } + } + } + } + + if ((*_excess)[target] < _min_cut) { + _min_cut = (*_excess)[target]; + for (NodeIt i(_graph); i != INVALID; ++i) { + _min_cut_map->set(i, false); + } + for (std::list::iterator it = _sets.back().begin(); + it != _sets.back().end(); ++it) { + Node n = _first[*it]; + while (n != INVALID) { + _min_cut_map->set(n, true); + n = (*_next)[n]; + } + } + } + + { + Node new_target; + if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) { + if ((*_next)[target] == INVALID) { + _last[(*_bucket)[target]] = (*_prev)[target]; + new_target = (*_prev)[target]; + } else { + _prev->set((*_next)[target], (*_prev)[target]); + new_target = (*_next)[target]; + } + if ((*_prev)[target] == INVALID) { + _first[(*_bucket)[target]] = (*_next)[target]; + } else { + _next->set((*_prev)[target], (*_next)[target]); + } + } else { + _sets.back().pop_back(); + if (_sets.back().empty()) { + _sets.pop_back(); + if (_sets.empty()) + break; + for (std::list::iterator it = _sets.back().begin(); + it != _sets.back().end(); ++it) { + _dormant[*it] = false; + } + } + new_target = _last[_sets.back().back()]; + } + + _bucket->set(target, 0); + + _source_set->set(target, true); + for (InArcIt a(_graph, target); a != INVALID; ++a) { + Value rem = (*_capacity)[a] - (*_flow)[a]; + if (!_tolerance.positive(rem)) continue; + Node v = _graph.source(a); + if (!(*_active)[v] && !(*_source_set)[v]) { + activate(v); + } + _excess->set(v, (*_excess)[v] + rem); + _flow->set(a, (*_capacity)[a]); + } + + for (OutArcIt a(_graph, target); a != INVALID; ++a) { + Value rem = (*_flow)[a]; + if (!_tolerance.positive(rem)) continue; + Node v = _graph.target(a); + if (!(*_active)[v] && !(*_source_set)[v]) { + activate(v); + } + _excess->set(v, (*_excess)[v] + rem); + _flow->set(a, 0); + } + + target = new_target; + if ((*_active)[target]) { + deactivate(target); + } + + _highest = _sets.back().begin(); + while (_highest != _sets.back().end() && + !(*_active)[_first[*_highest]]) { + ++_highest; + } + } + } + } + + public: + + /// \name Execution control + /// The simplest way to execute the algorithm is to use + /// one of the member functions called \c run(...). + /// \n + /// If you need more control on the execution, + /// first you must call \ref init(), then the \ref calculateIn() or + /// \ref calculateIn() functions. + + /// @{ + + /// \brief Initializes the internal data structures. + /// + /// Initializes the internal data structures. It creates + /// the maps, residual graph adaptors and some bucket structures + /// for the algorithm. + void init() { + init(NodeIt(_graph)); + } + + /// \brief Initializes the internal data structures. + /// + /// Initializes the internal data structures. It creates + /// the maps, residual graph adaptor and some bucket structures + /// for the algorithm. Node \c source is used as the push-relabel + /// algorithm's source. + void init(const Node& source) { + _source = source; + + _node_num = countNodes(_graph); + + _first.resize(_node_num + 1); + _last.resize(_node_num + 1); + + _dormant.resize(_node_num + 1); + + if (!_flow) { + _flow = new FlowMap(_graph); + } + if (!_next) { + _next = new typename Digraph::template NodeMap(_graph); + } + if (!_prev) { + _prev = new typename Digraph::template NodeMap(_graph); + } + if (!_active) { + _active = new typename Digraph::template NodeMap(_graph); + } + if (!_bucket) { + _bucket = new typename Digraph::template NodeMap(_graph); + } + if (!_excess) { + _excess = new ExcessMap(_graph); + } + if (!_source_set) { + _source_set = new SourceSetMap(_graph); + } + if (!_min_cut_map) { + _min_cut_map = new MinCutMap(_graph); + } + + _min_cut = std::numeric_limits::max(); + } + + + /// \brief Calculates a minimum cut with \f$ source \f$ on the + /// source-side. + /// + /// Calculates a minimum cut with \f$ source \f$ on the + /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source + /// \in X \f$ and minimal out-degree). + void calculateOut() { + findMinCutOut(); + } + + /// \brief Calculates a minimum cut with \f$ source \f$ on the + /// target-side. + /// + /// Calculates a minimum cut with \f$ source \f$ on the + /// target-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source + /// \in X \f$ and minimal out-degree). + void calculateIn() { + findMinCutIn(); + } + + + /// \brief Runs the algorithm. + /// + /// Runs the algorithm. It finds nodes \c source and \c target + /// arbitrarily and then calls \ref init(), \ref calculateOut() + /// and \ref calculateIn(). + void run() { + init(); + calculateOut(); + calculateIn(); + } + + /// \brief Runs the algorithm. + /// + /// Runs the algorithm. It uses the given \c source node, finds a + /// proper \c target and then calls the \ref init(), \ref + /// calculateOut() and \ref calculateIn(). + void run(const Node& s) { + init(s); + calculateOut(); + calculateIn(); + } + + /// @} + + /// \name Query Functions + /// The result of the %HaoOrlin algorithm + /// can be obtained using these functions. + /// \n + /// Before using these functions, either \ref run(), \ref + /// calculateOut() or \ref calculateIn() must be called. + + /// @{ + + /// \brief Returns the value of the minimum value cut. + /// + /// Returns the value of the minimum value cut. + Value minCutValue() const { + return _min_cut; + } + + + /// \brief Returns a minimum cut. + /// + /// Sets \c nodeMap to the characteristic vector of a minimum + /// value cut: it will give a nonempty set \f$ X\subsetneq V \f$ + /// with minimal out-degree (i.e. \c nodeMap will be true exactly + /// for the nodes of \f$ X \f$). \pre nodeMap should be a + /// bool-valued node-map. + template + Value minCutMap(NodeMap& nodeMap) const { + for (NodeIt it(_graph); it != INVALID; ++it) { + nodeMap.set(it, (*_min_cut_map)[it]); + } + return _min_cut; + } + + /// @} + + }; //class HaoOrlin + + +} //namespace lemon + +#endif //LEMON_HAO_ORLIN_H