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