Turn off 32 bit only tests, cont'd.
3 * This file is a part of LEMON, a generic C++ optimization library
5 * Copyright (C) 2003-2006
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
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.
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
19 #ifndef LEMON_HAO_ORLIN_H
20 #define LEMON_HAO_ORLIN_H
26 #include <lemon/maps.h>
27 #include <lemon/graph_utils.h>
28 #include <lemon/graph_adaptor.h>
29 #include <lemon/iterable_maps.h>
34 /// \brief Implementation of the Hao-Orlin algorithm.
36 /// Implementation of the HaoOrlin algorithms class for testing network
43 /// \brief %Hao-Orlin algorithm to find a minimum cut in directed graphs.
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 consists
47 /// of two phases: in the first phase it determines a minimum cut
48 /// with \f$ source \f$ on the source-side (i.e. a set \f$ X\subsetneq V \f$
49 /// with \f$ source \in X \f$ and minimal out-degree) and in the
50 /// second phase it determines a minimum cut with \f$ source \f$ on the
51 /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \notin X \f$
52 /// and minimal out-degree). Obviously, the smaller of these two
53 /// cuts will be a minimum cut of \f$ D \f$. The algorithm is a
54 /// modified push-relabel preflow algorithm and our implementation
55 /// calculates the minimum cut in \f$ O(n^3) \f$ time (we use the
56 /// highest-label rule). The purpose of such an algorithm is testing
57 /// network reliability. For an undirected graph with \f$ n \f$
58 /// nodes and \f$ e \f$ edges you can use the algorithm of Nagamochi
59 /// and Ibaraki which solves the undirected problem in
60 /// \f$ O(ne + n^2 \log(n)) \f$ time: it is implemented in the MinCut
64 /// \param _Graph is the graph type of the algorithm.
65 /// \param _CapacityMap is an edge map of capacities which should
66 /// be any numreric type. The default type is _Graph::EdgeMap<int>.
67 /// \param _Tolerance is the handler of the inexact computation. The
68 /// default type for this is Tolerance<typename CapacityMap::Value>.
70 /// \author Attila Bernath and Balazs Dezso
72 template <typename _Graph, typename _CapacityMap, typename _Tolerance>
74 template <typename _Graph,
75 typename _CapacityMap = typename _Graph::template EdgeMap<int>,
76 typename _Tolerance = Tolerance<typename _CapacityMap::Value> >
82 typedef _CapacityMap CapacityMap;
83 typedef _Tolerance Tolerance;
85 typedef typename CapacityMap::Value Value;
88 typedef typename Graph::Node Node;
89 typedef typename Graph::NodeIt NodeIt;
90 typedef typename Graph::EdgeIt EdgeIt;
91 typedef typename Graph::OutEdgeIt OutEdgeIt;
92 typedef typename Graph::InEdgeIt InEdgeIt;
96 const CapacityMap* _capacity;
98 typedef typename Graph::template EdgeMap<Value> FlowMap;
102 Node _source, _target;
105 typedef ResGraphAdaptor<const Graph, Value, CapacityMap,
106 FlowMap, Tolerance> OutResGraph;
107 typedef typename OutResGraph::Edge OutResEdge;
109 OutResGraph* _out_res_graph;
111 typedef typename Graph::template NodeMap<OutResEdge> OutCurrentEdgeMap;
112 OutCurrentEdgeMap* _out_current_edge;
114 typedef RevGraphAdaptor<const Graph> RevGraph;
115 RevGraph* _rev_graph;
117 typedef ResGraphAdaptor<const RevGraph, Value, CapacityMap,
118 FlowMap, Tolerance> InResGraph;
119 typedef typename InResGraph::Edge InResEdge;
121 InResGraph* _in_res_graph;
123 typedef typename Graph::template NodeMap<InResEdge> InCurrentEdgeMap;
124 InCurrentEdgeMap* _in_current_edge;
127 typedef IterableBoolMap<Graph, Node> WakeMap;
130 typedef typename Graph::template NodeMap<int> DistMap;
133 typedef typename Graph::template NodeMap<Value> ExcessMap;
136 typedef typename Graph::template NodeMap<bool> SourceSetMap;
137 SourceSetMap* _source_set;
139 std::vector<int> _level_size;
142 std::vector<std::list<Node> > _active_nodes;
145 std::vector<std::list<Node> > _dormant;
150 typedef typename Graph::template NodeMap<bool> MinCutMap;
151 MinCutMap* _min_cut_map;
153 Tolerance _tolerance;
157 /// \brief Constructor
159 /// Constructor of the algorithm class.
160 HaoOrlin(const Graph& graph, const CapacityMap& capacity,
161 const Tolerance& tolerance = Tolerance()) :
162 _graph(&graph), _capacity(&capacity),
163 _preflow(0), _source(), _target(),
164 _out_res_graph(0), _out_current_edge(0),
165 _rev_graph(0), _in_res_graph(0), _in_current_edge(0),
166 _wake(0),_dist(0), _excess(0), _source_set(0),
167 _highest_active(), _active_nodes(), _dormant_max(), _dormant(),
168 _min_cut(), _min_cut_map(0), _tolerance(tolerance) {}
174 if (_in_current_edge) {
175 delete _in_current_edge;
178 delete _in_res_graph;
183 if (_out_current_edge) {
184 delete _out_current_edge;
186 if (_out_res_graph) {
187 delete _out_res_graph;
208 template <typename ResGraph, typename EdgeMap>
209 void findMinCut(const Node& target, bool out,
210 ResGraph& res_graph, EdgeMap& current_edge) {
211 typedef typename ResGraph::Edge ResEdge;
212 typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
214 for (typename Graph::EdgeIt it(*_graph); it != INVALID; ++it) {
217 for (NodeIt it(*_graph); it != INVALID; ++it) {
221 (*_source_set)[it] = false;
223 res_graph.firstOut(current_edge[it], it);
227 (*_dist)[target] = 0;
229 for (ResOutEdgeIt it(res_graph, _source); it != INVALID; ++it) {
230 Value delta = res_graph.rescap(it);
231 if (!_tolerance.positive(delta)) continue;
233 (*_excess)[res_graph.source(it)] -= delta;
234 res_graph.augment(it, delta);
235 Node a = res_graph.target(it);
236 if (!_tolerance.positive((*_excess)[a]) &&
237 (*_wake)[a] && a != _target) {
238 _active_nodes[(*_dist)[a]].push_front(a);
239 if (_highest_active < (*_dist)[a]) {
240 _highest_active = (*_dist)[a];
243 (*_excess)[a] += delta;
246 _dormant[0].push_front(_source);
247 (*_source_set)[_source] = true;
249 (*_wake)[_source] = false;
252 _level_size[1] = _node_num - 1;
256 while ((n = findActiveNode()) != INVALID) {
258 while (_tolerance.positive((*_excess)[n]) &&
259 (e = findAdmissibleEdge(n, res_graph, current_edge))
262 if ((*_excess)[n] < res_graph.rescap(e)) {
263 delta = (*_excess)[n];
265 delta = res_graph.rescap(e);
266 res_graph.nextOut(current_edge[n]);
268 if (!_tolerance.positive(delta)) continue;
269 res_graph.augment(e, delta);
270 (*_excess)[res_graph.source(e)] -= delta;
271 Node a = res_graph.target(e);
272 if (!_tolerance.positive((*_excess)[a]) && a != _target) {
273 _active_nodes[(*_dist)[a]].push_front(a);
275 (*_excess)[a] += delta;
277 if (_tolerance.positive((*_excess)[n])) {
278 relabel(n, res_graph, current_edge);
282 Value current_value = cutValue(out);
283 if (_min_cut > current_value){
285 for (NodeIt it(*_graph); it != INVALID; ++it) {
286 _min_cut_map->set(it, !(*_wake)[it]);
289 for (NodeIt it(*_graph); it != INVALID; ++it) {
290 _min_cut_map->set(it, (*_wake)[it]);
294 _min_cut = current_value;
297 } while (selectNewSink(res_graph));
300 template <typename ResGraph, typename EdgeMap>
301 void relabel(const Node& n, ResGraph& res_graph, EdgeMap& current_edge) {
302 typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
305 if (_level_size[k] == 1) {
307 for (NodeIt it(*_graph); it != INVALID; ++it) {
308 if ((*_wake)[it] && (*_dist)[it] >= k) {
309 (*_wake)[it] = false;
310 _dormant[_dormant_max].push_front(it);
311 --_level_size[(*_dist)[it]];
316 int new_dist = _node_num;
317 for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
318 Node t = res_graph.target(e);
319 if ((*_wake)[t] && new_dist > (*_dist)[t]) {
320 new_dist = (*_dist)[t];
323 if (new_dist == _node_num) {
326 _dormant[_dormant_max].push_front(n);
327 --_level_size[(*_dist)[n]];
329 --_level_size[(*_dist)[n]];
330 (*_dist)[n] = new_dist + 1;
331 _highest_active = (*_dist)[n];
332 _active_nodes[_highest_active].push_front(n);
333 ++_level_size[(*_dist)[n]];
334 res_graph.firstOut(current_edge[n], n);
339 template <typename ResGraph>
340 bool selectNewSink(ResGraph& res_graph) {
341 typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
343 Node old_target = _target;
344 (*_wake)[_target] = false;
345 --_level_size[(*_dist)[_target]];
346 _dormant[0].push_front(_target);
347 (*_source_set)[_target] = true;
348 if ((int)_dormant[0].size() == _node_num){
353 bool wake_was_empty = false;
355 if(_wake->trueNum() == 0) {
356 while (!_dormant[_dormant_max].empty()){
357 (*_wake)[_dormant[_dormant_max].front()] = true;
358 ++_level_size[(*_dist)[_dormant[_dormant_max].front()]];
359 _dormant[_dormant_max].pop_front();
362 wake_was_empty = true;
365 int min_dist = std::numeric_limits<int>::max();
366 for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
367 if (min_dist > (*_dist)[it]){
369 min_dist = (*_dist)[it];
374 for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
375 if (_tolerance.positive((*_excess)[it])) {
376 if ((*_wake)[it] && it != _target) {
377 _active_nodes[(*_dist)[it]].push_front(it);
379 if (_highest_active < (*_dist)[it]) {
380 _highest_active = (*_dist)[it];
386 for (ResOutEdgeIt e(res_graph, old_target); e!=INVALID; ++e){
387 if (!(*_source_set)[res_graph.target(e)]) {
388 Value delta = res_graph.rescap(e);
389 if (!_tolerance.positive(delta)) continue;
390 res_graph.augment(e, delta);
391 (*_excess)[res_graph.source(e)] -= delta;
392 Node a = res_graph.target(e);
393 if (!_tolerance.positive((*_excess)[a]) &&
394 (*_wake)[a] && a != _target) {
395 _active_nodes[(*_dist)[a]].push_front(a);
396 if (_highest_active < (*_dist)[a]) {
397 _highest_active = (*_dist)[a];
400 (*_excess)[a] += delta;
407 Node findActiveNode() {
408 while (_highest_active > 0 && _active_nodes[_highest_active].empty()){
411 if( _highest_active > 0) {
412 Node n = _active_nodes[_highest_active].front();
413 _active_nodes[_highest_active].pop_front();
420 template <typename ResGraph, typename EdgeMap>
421 typename ResGraph::Edge findAdmissibleEdge(const Node& n,
423 EdgeMap& current_edge) {
424 typedef typename ResGraph::Edge ResEdge;
425 ResEdge e = current_edge[n];
426 while (e != INVALID &&
427 ((*_dist)[n] <= (*_dist)[res_graph.target(e)] ||
428 !(*_wake)[res_graph.target(e)])) {
429 res_graph.nextOut(e);
439 Value cutValue(bool out) {
442 for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
443 for (InEdgeIt e(*_graph, it); e != INVALID; ++e) {
444 if (!(*_wake)[_graph->source(e)]){
445 value += (*_capacity)[e];
450 for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
451 for (OutEdgeIt e(*_graph, it); e != INVALID; ++e) {
452 if (!(*_wake)[_graph->target(e)]){
453 value += (*_capacity)[e];
464 /// \name Execution control
465 /// The simplest way to execute the algorithm is to use
466 /// one of the member functions called \c run(...).
468 /// If you need more control on the execution,
469 /// first you must call \ref init(), then the \ref calculateIn() or
470 /// \ref calculateIn() functions.
474 /// \brief Initializes the internal data structures.
476 /// Initializes the internal data structures. It creates
477 /// the maps, residual graph adaptors and some bucket structures
478 /// for the algorithm.
480 init(NodeIt(*_graph));
483 /// \brief Initializes the internal data structures.
485 /// Initializes the internal data structures. It creates
486 /// the maps, residual graph adaptor and some bucket structures
487 /// for the algorithm. Node \c source is used as the push-relabel
488 /// algorithm's source.
489 void init(const Node& source) {
491 _node_num = countNodes(*_graph);
493 _dormant.resize(_node_num);
494 _level_size.resize(_node_num, 0);
495 _active_nodes.resize(_node_num);
498 _preflow = new FlowMap(*_graph);
501 _wake = new WakeMap(*_graph);
504 _dist = new DistMap(*_graph);
507 _excess = new ExcessMap(*_graph);
510 _source_set = new SourceSetMap(*_graph);
512 if (!_out_res_graph) {
513 _out_res_graph = new OutResGraph(*_graph, *_capacity, *_preflow);
515 if (!_out_current_edge) {
516 _out_current_edge = new OutCurrentEdgeMap(*_graph);
519 _rev_graph = new RevGraph(*_graph);
521 if (!_in_res_graph) {
522 _in_res_graph = new InResGraph(*_rev_graph, *_capacity, *_preflow);
524 if (!_in_current_edge) {
525 _in_current_edge = new InCurrentEdgeMap(*_graph);
528 _min_cut_map = new MinCutMap(*_graph);
531 _min_cut = std::numeric_limits<Value>::max();
535 /// \brief Calculates a minimum cut with \f$ source \f$ on the
538 /// \brief Calculates a minimum cut with \f$ source \f$ on the
539 /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
540 /// and minimal out-degree).
541 void calculateOut() {
542 for (NodeIt it(*_graph); it != INVALID; ++it) {
550 /// \brief Calculates a minimum cut with \f$ source \f$ on the
553 /// \brief Calculates a minimum cut with \f$ source \f$ on the
554 /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
555 /// and minimal out-degree). The \c target is the initial target
556 /// for the push-relabel algorithm.
557 void calculateOut(const Node& target) {
558 findMinCut(target, true, *_out_res_graph, *_out_current_edge);
561 /// \brief Calculates a minimum cut with \f$ source \f$ on the
564 /// \brief Calculates a minimum cut with \f$ source \f$ on the
565 /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with
566 /// \f$ source \notin X \f$
567 /// and minimal out-degree).
569 for (NodeIt it(*_graph); it != INVALID; ++it) {
577 /// \brief Calculates a minimum cut with \f$ source \f$ on the
580 /// \brief Calculates a minimum cut with \f$ source \f$ on the
581 /// sink-side (i.e. a set \f$ X\subsetneq V
582 /// \f$ with \f$ source \notin X \f$ and minimal out-degree).
583 /// The \c target is the initial
584 /// target for the push-relabel algorithm.
585 void calculateIn(const Node& target) {
586 findMinCut(target, false, *_in_res_graph, *_in_current_edge);
589 /// \brief Runs the algorithm.
591 /// Runs the algorithm. It finds nodes \c source and \c target
592 /// arbitrarily and then calls \ref init(), \ref calculateOut()
593 /// and \ref calculateIn().
596 for (NodeIt it(*_graph); it != INVALID; ++it) {
605 /// \brief Runs the algorithm.
607 /// Runs the algorithm. It uses the given \c source node, finds a
608 /// proper \c target and then calls the \ref init(), \ref
609 /// calculateOut() and \ref calculateIn().
610 void run(const Node& s) {
612 for (NodeIt it(*_graph); it != INVALID; ++it) {
621 /// \brief Runs the algorithm.
623 /// Runs the algorithm. It just calls the \ref init() and then
624 /// \ref calculateOut() and \ref calculateIn().
625 void run(const Node& s, const Node& t) {
633 /// \name Query Functions The result of the %HaoOrlin algorithm
634 /// can be obtained using these functions.
636 /// Before the use of these functions, either \ref run(), \ref
637 /// calculateOut() or \ref calculateIn() must be called.
641 /// \brief Returns the value of the minimum value cut.
643 /// Returns the value of the minimum value cut.
644 Value minCut() const {
649 /// \brief Returns a minimum cut.
651 /// Sets \c nodeMap to the characteristic vector of a minimum
652 /// value cut: it will give a nonempty set \f$ X\subsetneq V \f$
653 /// with minimal out-degree (i.e. \c nodeMap will be true exactly
654 /// for the nodes of \f$ X \f$. \pre nodeMap should be a
655 /// bool-valued node-map.
656 template <typename NodeMap>
657 Value minCut(NodeMap& nodeMap) const {
658 for (NodeIt it(*_graph); it != INVALID; ++it) {
659 nodeMap.set(it, (*_min_cut_map)[it]);
671 #endif //LEMON_HAO_ORLIN_H