Consider the CPXMIP_OPTIMAL_TOL status as OPTIMAL too.
3 * This file is a part of LEMON, a generic C++ optimization library
5 * Copyright (C) 2003-2007
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
31 #include <lemon/maps.h>
32 #include <lemon/graph_utils.h>
33 #include <lemon/graph_adaptor.h>
34 #include <lemon/iterable_maps.h>
38 /// \brief Implementation of the Hao-Orlin algorithm.
40 /// Implementation of the HaoOrlin algorithms class for testing network
47 /// \brief %Hao-Orlin algorithm to find a minimum cut in directed graphs.
49 /// Hao-Orlin calculates a minimum cut in a directed graph
50 /// \f$ D=(V,A) \f$. It takes a fixed node \f$ source \in V \f$ and consists
51 /// of two phases: in the first phase it determines a minimum cut
52 /// with \f$ source \f$ on the source-side (i.e. a set \f$ X\subsetneq V \f$
53 /// with \f$ source \in X \f$ and minimal out-degree) and in the
54 /// second phase it determines a minimum cut with \f$ source \f$ on the
55 /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \notin X \f$
56 /// and minimal out-degree). Obviously, the smaller of these two
57 /// cuts will be a minimum cut of \f$ D \f$. The algorithm is a
58 /// modified push-relabel preflow algorithm and our implementation
59 /// calculates the minimum cut in \f$ O(n^3) \f$ time (we use the
60 /// highest-label rule). The purpose of such an algorithm is testing
61 /// network reliability. For an undirected graph with \f$ n \f$
62 /// nodes and \f$ e \f$ edges you can use the algorithm of Nagamochi
63 /// and Ibaraki which solves the undirected problem in
64 /// \f$ O(ne + n^2 \log(n)) \f$ time: it is implemented in the MinCut
68 /// \param _Graph is the graph type of the algorithm.
69 /// \param _CapacityMap is an edge map of capacities which should
70 /// be any numreric type. The default type is _Graph::EdgeMap<int>.
71 /// \param _Tolerance is the handler of the inexact computation. The
72 /// default type for this is Tolerance<typename CapacityMap::Value>.
74 /// \author Attila Bernath and Balazs Dezso
76 template <typename _Graph, typename _CapacityMap, typename _Tolerance>
78 template <typename _Graph,
79 typename _CapacityMap = typename _Graph::template EdgeMap<int>,
80 typename _Tolerance = Tolerance<typename _CapacityMap::Value> >
86 typedef _CapacityMap CapacityMap;
87 typedef _Tolerance Tolerance;
89 typedef typename CapacityMap::Value Value;
92 typedef typename Graph::Node Node;
93 typedef typename Graph::NodeIt NodeIt;
94 typedef typename Graph::EdgeIt EdgeIt;
95 typedef typename Graph::OutEdgeIt OutEdgeIt;
96 typedef typename Graph::InEdgeIt InEdgeIt;
100 const CapacityMap* _capacity;
102 typedef typename Graph::template EdgeMap<Value> FlowMap;
106 Node _source, _target;
109 typedef ResGraphAdaptor<const Graph, Value, CapacityMap,
110 FlowMap, Tolerance> OutResGraph;
111 typedef typename OutResGraph::Edge OutResEdge;
113 OutResGraph* _out_res_graph;
115 typedef RevGraphAdaptor<const Graph> RevGraph;
116 RevGraph* _rev_graph;
118 typedef ResGraphAdaptor<const RevGraph, Value, CapacityMap,
119 FlowMap, Tolerance> InResGraph;
120 typedef typename InResGraph::Edge InResEdge;
122 InResGraph* _in_res_graph;
124 typedef IterableBoolMap<Graph, Node> WakeMap;
127 typedef typename Graph::template NodeMap<int> DistMap;
130 typedef typename Graph::template NodeMap<Value> ExcessMap;
133 typedef typename Graph::template NodeMap<bool> SourceSetMap;
134 SourceSetMap* _source_set;
136 std::vector<int> _level_size;
139 std::vector<std::list<Node> > _active_nodes;
142 std::vector<std::list<Node> > _dormant;
147 typedef typename Graph::template NodeMap<bool> MinCutMap;
148 MinCutMap* _min_cut_map;
150 Tolerance _tolerance;
154 /// \brief Constructor
156 /// Constructor of the algorithm class.
157 HaoOrlin(const Graph& graph, const CapacityMap& capacity,
158 const Tolerance& tolerance = Tolerance()) :
159 _graph(&graph), _capacity(&capacity),
160 _preflow(0), _source(), _target(),
161 _out_res_graph(0), _rev_graph(0), _in_res_graph(0),
162 _wake(0),_dist(0), _excess(0), _source_set(0),
163 _highest_active(), _active_nodes(), _dormant_max(), _dormant(),
164 _min_cut(), _min_cut_map(0), _tolerance(tolerance) {}
171 delete _in_res_graph;
176 if (_out_res_graph) {
177 delete _out_res_graph;
198 template <typename ResGraph>
199 void findMinCut(const Node& target, bool out, ResGraph& res_graph) {
200 typedef typename ResGraph::Edge ResEdge;
201 typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
203 for (typename Graph::EdgeIt it(*_graph); it != INVALID; ++it) {
206 for (NodeIt it(*_graph); it != INVALID; ++it) {
210 (*_source_set)[it] = false;
213 _dormant[0].push_front(_source);
214 (*_source_set)[_source] = true;
216 (*_wake)[_source] = false;
219 _level_size[1] = _node_num - 1;
222 (*_dist)[target] = 0;
224 for (ResOutEdgeIt it(res_graph, _source); it != INVALID; ++it) {
225 Value delta = res_graph.rescap(it);
226 (*_excess)[_source] -= delta;
227 res_graph.augment(it, delta);
228 Node a = res_graph.target(it);
229 if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) {
230 _active_nodes[(*_dist)[a]].push_front(a);
231 if (_highest_active < (*_dist)[a]) {
232 _highest_active = (*_dist)[a];
235 (*_excess)[a] += delta;
241 while ((n = findActiveNode()) != INVALID) {
242 for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
243 Node a = res_graph.target(e);
244 if ((*_dist)[a] >= (*_dist)[n] || !(*_wake)[a]) continue;
245 Value delta = res_graph.rescap(e);
246 if (_tolerance.positive((*_excess)[n] - delta)) {
247 (*_excess)[n] -= delta;
249 delta = (*_excess)[n];
252 res_graph.augment(e, delta);
253 if ((*_excess)[a] == 0 && a != _target) {
254 _active_nodes[(*_dist)[a]].push_front(a);
256 (*_excess)[a] += delta;
257 if ((*_excess)[n] == 0) break;
259 if ((*_excess)[n] != 0) {
260 relabel(n, res_graph);
264 Value current_value = cutValue(out);
265 if (_min_cut > current_value){
267 for (NodeIt it(*_graph); it != INVALID; ++it) {
268 _min_cut_map->set(it, !(*_wake)[it]);
271 for (NodeIt it(*_graph); it != INVALID; ++it) {
272 _min_cut_map->set(it, (*_wake)[it]);
276 _min_cut = current_value;
279 } while (selectNewSink(res_graph));
282 template <typename ResGraph>
283 void relabel(const Node& n, ResGraph& res_graph) {
284 typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
287 if (_level_size[k] == 1) {
289 for (NodeIt it(*_graph); it != INVALID; ++it) {
290 if ((*_wake)[it] && (*_dist)[it] >= k) {
291 (*_wake)[it] = false;
292 _dormant[_dormant_max].push_front(it);
293 --_level_size[(*_dist)[it]];
298 int new_dist = _node_num;
299 for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) {
300 Node t = res_graph.target(e);
301 if ((*_wake)[t] && new_dist > (*_dist)[t]) {
302 new_dist = (*_dist)[t];
305 if (new_dist == _node_num) {
308 _dormant[_dormant_max].push_front(n);
309 --_level_size[(*_dist)[n]];
311 --_level_size[(*_dist)[n]];
312 (*_dist)[n] = new_dist + 1;
313 _highest_active = (*_dist)[n];
314 _active_nodes[_highest_active].push_front(n);
315 ++_level_size[(*_dist)[n]];
320 template <typename ResGraph>
321 bool selectNewSink(ResGraph& res_graph) {
322 typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
324 Node old_target = _target;
325 (*_wake)[_target] = false;
326 --_level_size[(*_dist)[_target]];
327 _dormant[0].push_front(_target);
328 (*_source_set)[_target] = true;
329 if (int(_dormant[0].size()) == _node_num){
334 bool wake_was_empty = false;
336 if(_wake->trueNum() == 0) {
337 while (!_dormant[_dormant_max].empty()){
338 (*_wake)[_dormant[_dormant_max].front()] = true;
339 ++_level_size[(*_dist)[_dormant[_dormant_max].front()]];
340 _dormant[_dormant_max].pop_front();
343 wake_was_empty = true;
346 int min_dist = std::numeric_limits<int>::max();
347 for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
348 if (min_dist > (*_dist)[it]){
350 min_dist = (*_dist)[it];
355 for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
356 if ((*_excess)[it] != 0 && it != _target) {
357 _active_nodes[(*_dist)[it]].push_front(it);
358 if (_highest_active < (*_dist)[it]) {
359 _highest_active = (*_dist)[it];
366 for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e){
367 Node a = res_graph.target(e);
368 if (!(*_wake)[a]) continue;
369 Value delta = res_graph.rescap(e);
370 res_graph.augment(e, delta);
371 (*_excess)[n] -= delta;
372 if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) {
373 _active_nodes[(*_dist)[a]].push_front(a);
374 if (_highest_active < (*_dist)[a]) {
375 _highest_active = (*_dist)[a];
378 (*_excess)[a] += delta;
384 Node findActiveNode() {
385 while (_highest_active > 0 && _active_nodes[_highest_active].empty()){
388 if( _highest_active > 0) {
389 Node n = _active_nodes[_highest_active].front();
390 _active_nodes[_highest_active].pop_front();
397 Value cutValue(bool out) {
400 for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
401 for (InEdgeIt e(*_graph, it); e != INVALID; ++e) {
402 if (!(*_wake)[_graph->source(e)]){
403 value += (*_capacity)[e];
408 for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
409 for (OutEdgeIt e(*_graph, it); e != INVALID; ++e) {
410 if (!(*_wake)[_graph->target(e)]){
411 value += (*_capacity)[e];
422 /// \name Execution control
423 /// The simplest way to execute the algorithm is to use
424 /// one of the member functions called \c run(...).
426 /// If you need more control on the execution,
427 /// first you must call \ref init(), then the \ref calculateIn() or
428 /// \ref calculateIn() functions.
432 /// \brief Initializes the internal data structures.
434 /// Initializes the internal data structures. It creates
435 /// the maps, residual graph adaptors and some bucket structures
436 /// for the algorithm.
438 init(NodeIt(*_graph));
441 /// \brief Initializes the internal data structures.
443 /// Initializes the internal data structures. It creates
444 /// the maps, residual graph adaptor and some bucket structures
445 /// for the algorithm. Node \c source is used as the push-relabel
446 /// algorithm's source.
447 void init(const Node& source) {
449 _node_num = countNodes(*_graph);
451 _dormant.resize(_node_num);
452 _level_size.resize(_node_num, 0);
453 _active_nodes.resize(_node_num);
456 _preflow = new FlowMap(*_graph);
459 _wake = new WakeMap(*_graph);
462 _dist = new DistMap(*_graph);
465 _excess = new ExcessMap(*_graph);
468 _source_set = new SourceSetMap(*_graph);
470 if (!_out_res_graph) {
471 _out_res_graph = new OutResGraph(*_graph, *_capacity, *_preflow);
474 _rev_graph = new RevGraph(*_graph);
476 if (!_in_res_graph) {
477 _in_res_graph = new InResGraph(*_rev_graph, *_capacity, *_preflow);
480 _min_cut_map = new MinCutMap(*_graph);
483 _min_cut = std::numeric_limits<Value>::max();
487 /// \brief Calculates a minimum cut with \f$ source \f$ on the
490 /// \brief Calculates a minimum cut with \f$ source \f$ on the
491 /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
492 /// and minimal out-degree).
493 void calculateOut() {
494 for (NodeIt it(*_graph); it != INVALID; ++it) {
502 /// \brief Calculates a minimum cut with \f$ source \f$ on the
505 /// \brief Calculates a minimum cut with \f$ source \f$ on the
506 /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$
507 /// and minimal out-degree). The \c target is the initial target
508 /// for the push-relabel algorithm.
509 void calculateOut(const Node& target) {
510 findMinCut(target, true, *_out_res_graph);
513 /// \brief Calculates a minimum cut with \f$ source \f$ on the
516 /// \brief Calculates a minimum cut with \f$ source \f$ on the
517 /// sink-side (i.e. a set \f$ X\subsetneq V \f$ with
518 /// \f$ source \notin X \f$
519 /// and minimal out-degree).
521 for (NodeIt it(*_graph); it != INVALID; ++it) {
529 /// \brief Calculates a minimum cut with \f$ source \f$ on the
532 /// \brief Calculates a minimum cut with \f$ source \f$ on the
533 /// sink-side (i.e. a set \f$ X\subsetneq V
534 /// \f$ with \f$ source \notin X \f$ and minimal out-degree).
535 /// The \c target is the initial
536 /// target for the push-relabel algorithm.
537 void calculateIn(const Node& target) {
538 findMinCut(target, false, *_in_res_graph);
541 /// \brief Runs the algorithm.
543 /// Runs the algorithm. It finds nodes \c source and \c target
544 /// arbitrarily and then calls \ref init(), \ref calculateOut()
545 /// and \ref calculateIn().
548 for (NodeIt it(*_graph); it != INVALID; ++it) {
557 /// \brief Runs the algorithm.
559 /// Runs the algorithm. It uses the given \c source node, finds a
560 /// proper \c target and then calls the \ref init(), \ref
561 /// calculateOut() and \ref calculateIn().
562 void run(const Node& s) {
564 for (NodeIt it(*_graph); it != INVALID; ++it) {
573 /// \brief Runs the algorithm.
575 /// Runs the algorithm. It just calls the \ref init() and then
576 /// \ref calculateOut() and \ref calculateIn().
577 void run(const Node& s, const Node& t) {
585 /// \name Query Functions
586 /// The result of the %HaoOrlin algorithm
587 /// can be obtained using these functions.
589 /// Before using these functions, either \ref run(), \ref
590 /// calculateOut() or \ref calculateIn() must be called.
594 /// \brief Returns the value of the minimum value cut.
596 /// Returns the value of the minimum value cut.
597 Value minCut() const {
602 /// \brief Returns a minimum cut.
604 /// Sets \c nodeMap to the characteristic vector of a minimum
605 /// value cut: it will give a nonempty set \f$ X\subsetneq V \f$
606 /// with minimal out-degree (i.e. \c nodeMap will be true exactly
607 /// for the nodes of \f$ X \f$). \pre nodeMap should be a
608 /// bool-valued node-map.
609 template <typename NodeMap>
610 Value minCut(NodeMap& nodeMap) const {
611 for (NodeIt it(*_graph); it != INVALID; ++it) {
612 nodeMap.set(it, (*_min_cut_map)[it]);
624 #endif //LEMON_HAO_ORLIN_H