/* -*- C++ -*- * * This file is a part of LEMON, a generic C++ optimization library * * Copyright (C) 2003-2007 * 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 #include #include #include /// \file /// \ingroup min_cut /// \brief Implementation of the Hao-Orlin algorithm. /// /// Implementation of the HaoOrlin algorithms 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^3) \f$ time (we use the /// highest-label rule). The purpose of such an algorithm is testing /// network reliability. For an undirected graph with \f$n \f$ /// nodes and \f$e \f$ edges you can use the algorithm of Nagamochi /// and Ibaraki which solves the undirected problem in /// \f$O(ne + n^2 \log(n)) \f$ time: it is implemented in the MinCut /// algorithm /// class. /// /// \param _Graph 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 _Graph::EdgeMap. /// \param _Tolerance is the handler of the inexact computation. The /// default type for this is Tolerance. /// /// \author Attila Bernath and Balazs Dezso #ifdef DOXYGEN template #else template , typename _Tolerance = Tolerance > #endif class HaoOrlin { protected: typedef _Graph Graph; typedef _CapacityMap CapacityMap; typedef _Tolerance Tolerance; typedef typename CapacityMap::Value Value; typedef typename Graph::Node Node; typedef typename Graph::NodeIt NodeIt; typedef typename Graph::EdgeIt EdgeIt; typedef typename Graph::OutEdgeIt OutEdgeIt; typedef typename Graph::InEdgeIt InEdgeIt; const Graph* _graph; const CapacityMap* _capacity; typedef typename Graph::template EdgeMap FlowMap; FlowMap* _preflow; Node _source, _target; int _node_num; typedef ResGraphAdaptor OutResGraph; typedef typename OutResGraph::Edge OutResEdge; OutResGraph* _out_res_graph; typedef RevGraphAdaptor RevGraph; RevGraph* _rev_graph; typedef ResGraphAdaptor InResGraph; typedef typename InResGraph::Edge InResEdge; InResGraph* _in_res_graph; typedef IterableBoolMap WakeMap; WakeMap* _wake; typedef typename Graph::template NodeMap DistMap; DistMap* _dist; typedef typename Graph::template NodeMap ExcessMap; ExcessMap* _excess; typedef typename Graph::template NodeMap SourceSetMap; SourceSetMap* _source_set; std::vector _level_size; int _highest_active; std::vector > _active_nodes; int _dormant_max; std::vector > _dormant; Value _min_cut; typedef typename Graph::template NodeMap MinCutMap; MinCutMap* _min_cut_map; Tolerance _tolerance; public: /// \brief Constructor /// /// Constructor of the algorithm class. HaoOrlin(const Graph& graph, const CapacityMap& capacity, const Tolerance& tolerance = Tolerance()) : _graph(&graph), _capacity(&capacity), _preflow(0), _source(), _target(), _out_res_graph(0), _rev_graph(0), _in_res_graph(0), _wake(0),_dist(0), _excess(0), _source_set(0), _highest_active(), _active_nodes(), _dormant_max(), _dormant(), _min_cut(), _min_cut_map(0), _tolerance(tolerance) {} ~HaoOrlin() { if (_min_cut_map) { delete _min_cut_map; } if (_in_res_graph) { delete _in_res_graph; } if (_rev_graph) { delete _rev_graph; } if (_out_res_graph) { delete _out_res_graph; } if (_source_set) { delete _source_set; } if (_excess) { delete _excess; } if (_dist) { delete _dist; } if (_wake) { delete _wake; } if (_preflow) { delete _preflow; } } private: template void findMinCut(const Node& target, bool out, ResGraph& res_graph) { typedef typename ResGraph::Edge ResEdge; typedef typename ResGraph::OutEdgeIt ResOutEdgeIt; for (typename Graph::EdgeIt it(*_graph); it != INVALID; ++it) { (*_preflow)[it] = 0; } for (NodeIt it(*_graph); it != INVALID; ++it) { (*_wake)[it] = true; (*_dist)[it] = 1; (*_excess)[it] = 0; (*_source_set)[it] = false; } _dormant[0].push_front(_source); (*_source_set)[_source] = true; _dormant_max = 0; (*_wake)[_source] = false; _level_size[0] = 1; _level_size[1] = _node_num - 1; _target = target; (*_dist)[target] = 0; for (ResOutEdgeIt it(res_graph, _source); it != INVALID; ++it) { Value delta = res_graph.rescap(it); (*_excess)[_source] -= delta; res_graph.augment(it, delta); Node a = res_graph.target(it); if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) { _active_nodes[(*_dist)[a]].push_front(a); if (_highest_active < (*_dist)[a]) { _highest_active = (*_dist)[a]; } } (*_excess)[a] += delta; } do { Node n; while ((n = findActiveNode()) != INVALID) { for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) { Node a = res_graph.target(e); if ((*_dist)[a] >= (*_dist)[n] || !(*_wake)[a]) continue; Value delta = res_graph.rescap(e); if (_tolerance.positive((*_excess)[n] - delta)) { (*_excess)[n] -= delta; } else { delta = (*_excess)[n]; (*_excess)[n] = 0; } res_graph.augment(e, delta); if ((*_excess)[a] == 0 && a != _target) { _active_nodes[(*_dist)[a]].push_front(a); } (*_excess)[a] += delta; if ((*_excess)[n] == 0) break; } if ((*_excess)[n] != 0) { relabel(n, res_graph); } } Value current_value = cutValue(out); if (_min_cut > current_value){ if (out) { for (NodeIt it(*_graph); it != INVALID; ++it) { _min_cut_map->set(it, !(*_wake)[it]); } } else { for (NodeIt it(*_graph); it != INVALID; ++it) { _min_cut_map->set(it, (*_wake)[it]); } } _min_cut = current_value; } } while (selectNewSink(res_graph)); } template void relabel(const Node& n, ResGraph& res_graph) { typedef typename ResGraph::OutEdgeIt ResOutEdgeIt; int k = (*_dist)[n]; if (_level_size[k] == 1) { ++_dormant_max; for (NodeIt it(*_graph); it != INVALID; ++it) { if ((*_wake)[it] && (*_dist)[it] >= k) { (*_wake)[it] = false; _dormant[_dormant_max].push_front(it); --_level_size[(*_dist)[it]]; } } --_highest_active; } else { int new_dist = _node_num; for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e) { Node t = res_graph.target(e); if ((*_wake)[t] && new_dist > (*_dist)[t]) { new_dist = (*_dist)[t]; } } if (new_dist == _node_num) { ++_dormant_max; (*_wake)[n] = false; _dormant[_dormant_max].push_front(n); --_level_size[(*_dist)[n]]; } else { --_level_size[(*_dist)[n]]; (*_dist)[n] = new_dist + 1; _highest_active = (*_dist)[n]; _active_nodes[_highest_active].push_front(n); ++_level_size[(*_dist)[n]]; } } } template bool selectNewSink(ResGraph& res_graph) { typedef typename ResGraph::OutEdgeIt ResOutEdgeIt; Node old_target = _target; (*_wake)[_target] = false; --_level_size[(*_dist)[_target]]; _dormant[0].push_front(_target); (*_source_set)[_target] = true; if (int(_dormant[0].size()) == _node_num){ _dormant[0].clear(); return false; } bool wake_was_empty = false; if(_wake->trueNum() == 0) { while (!_dormant[_dormant_max].empty()){ (*_wake)[_dormant[_dormant_max].front()] = true; ++_level_size[(*_dist)[_dormant[_dormant_max].front()]]; _dormant[_dormant_max].pop_front(); } --_dormant_max; wake_was_empty = true; } int min_dist = std::numeric_limits::max(); for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) { if (min_dist > (*_dist)[it]){ _target = it; min_dist = (*_dist)[it]; } } if (wake_was_empty){ for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) { if ((*_excess)[it] != 0 && it != _target) { _active_nodes[(*_dist)[it]].push_front(it); if (_highest_active < (*_dist)[it]) { _highest_active = (*_dist)[it]; } } } } Node n = old_target; for (ResOutEdgeIt e(res_graph, n); e != INVALID; ++e){ Node a = res_graph.target(e); if (!(*_wake)[a]) continue; Value delta = res_graph.rescap(e); res_graph.augment(e, delta); (*_excess)[n] -= delta; if ((*_excess)[a] == 0 && (*_wake)[a] && a != _target) { _active_nodes[(*_dist)[a]].push_front(a); if (_highest_active < (*_dist)[a]) { _highest_active = (*_dist)[a]; } } (*_excess)[a] += delta; } return true; } Node findActiveNode() { while (_highest_active > 0 && _active_nodes[_highest_active].empty()){ --_highest_active; } if( _highest_active > 0) { Node n = _active_nodes[_highest_active].front(); _active_nodes[_highest_active].pop_front(); return n; } else { return INVALID; } } Value cutValue(bool out) { Value value = 0; if (out) { for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) { for (InEdgeIt e(*_graph, it); e != INVALID; ++e) { if (!(*_wake)[_graph->source(e)]){ value += (*_capacity)[e]; } } } } else { for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) { for (OutEdgeIt e(*_graph, it); e != INVALID; ++e) { if (!(*_wake)[_graph->target(e)]){ value += (*_capacity)[e]; } } } } return value; } 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); _dormant.resize(_node_num); _level_size.resize(_node_num, 0); _active_nodes.resize(_node_num); if (!_preflow) { _preflow = new FlowMap(*_graph); } if (!_wake) { _wake = new WakeMap(*_graph); } if (!_dist) { _dist = new DistMap(*_graph); } if (!_excess) { _excess = new ExcessMap(*_graph); } if (!_source_set) { _source_set = new SourceSetMap(*_graph); } if (!_out_res_graph) { _out_res_graph = new OutResGraph(*_graph, *_capacity, *_preflow); } if (!_rev_graph) { _rev_graph = new RevGraph(*_graph); } if (!_in_res_graph) { _in_res_graph = new InResGraph(*_rev_graph, *_capacity, *_preflow); } 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. /// /// \brief 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() { for (NodeIt it(*_graph); it != INVALID; ++it) { if (it != _source) { calculateOut(it); return; } } } /// \brief Calculates a minimum cut with \f$source \f$ on the /// source-side. /// /// \brief 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). The \c target is the initial target /// for the push-relabel algorithm. void calculateOut(const Node& target) { findMinCut(target, true, *_out_res_graph); } /// \brief Calculates a minimum cut with \f$source \f$ on the /// sink-side. /// /// \brief Calculates 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). void calculateIn() { for (NodeIt it(*_graph); it != INVALID; ++it) { if (it != _source) { calculateIn(it); return; } } } /// \brief Calculates a minimum cut with \f$source \f$ on the /// sink-side. /// /// \brief Calculates 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). /// The \c target is the initial /// target for the push-relabel algorithm. void calculateIn(const Node& target) { findMinCut(target, false, *_in_res_graph); } /// \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(); for (NodeIt it(*_graph); it != INVALID; ++it) { if (it != _source) { calculateOut(it); calculateIn(it); return; } } } /// \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); for (NodeIt it(*_graph); it != INVALID; ++it) { if (it != _source) { calculateOut(it); calculateIn(it); return; } } } /// \brief Runs the algorithm. /// /// Runs the algorithm. It just calls the \ref init() and then /// \ref calculateOut() and \ref calculateIn(). void run(const Node& s, const Node& t) { init(s); calculateOut(t); calculateIn(t); } /// @} /// \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 minCut() 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 minCut(NodeMap& nodeMap) const { for (NodeIt it(*_graph); it != INVALID; ++it) { nodeMap.set(it, (*_min_cut_map)[it]); } return minCut(); } /// @} }; //class HaoOrlin } //namespace lemon #endif //LEMON_HAO_ORLIN_H