1.1 --- a/lemon/Makefile.am Thu Sep 07 14:04:31 2006 +0000
1.2 +++ b/lemon/Makefile.am Thu Sep 07 14:16:47 2006 +0000
1.3 @@ -57,6 +57,7 @@
1.4 lemon/graph_utils.h \
1.5 lemon/graph_writer.h \
1.6 lemon/grid_ugraph.h \
1.7 + lemon/hao_orlin.h \
1.8 lemon/hypercube_graph.h \
1.9 lemon/iterable_maps.h \
1.10 lemon/johnson.h \
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/lemon/hao_orlin.h Thu Sep 07 14:16:47 2006 +0000
2.3 @@ -0,0 +1,521 @@
2.4 +/* -*- C++ -*-
2.5 + * lemon/hao_orlin.h - Part of LEMON, a generic C++ optimization library
2.6 + *
2.7 + * Copyright (C) 2005 Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
2.8 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
2.9 + *
2.10 + * Permission to use, modify and distribute this software is granted
2.11 + * provided that this copyright notice appears in all copies. For
2.12 + * precise terms see the accompanying LICENSE file.
2.13 + *
2.14 + * This software is provided "AS IS" with no warranty of any kind,
2.15 + * express or implied, and with no claim as to its suitability for any
2.16 + * purpose.
2.17 + *
2.18 + */
2.19 +
2.20 +#ifndef LEMON_HAO_ORLIN_H
2.21 +#define LEMON_HAO_ORLIN_H
2.22 +
2.23 +#include <vector>
2.24 +#include <queue>
2.25 +#include <limits>
2.26 +
2.27 +#include <lemon/maps.h>
2.28 +#include <lemon/graph_utils.h>
2.29 +#include <lemon/graph_adaptor.h>
2.30 +#include <lemon/iterable_maps.h>
2.31 +
2.32 +
2.33 +/// \file
2.34 +/// \ingroup flowalgs
2.35 +/// Implementation of the Hao-Orlin algorithms class for testing network
2.36 +/// reliability.
2.37 +
2.38 +namespace lemon {
2.39 +
2.40 + /// \addtogroup flowalgs
2.41 + /// @{
2.42 +
2.43 + /// %Hao-Orlin algorithm for calculate minimum cut in directed graphs.
2.44 + ///
2.45 + /// Hao-Orlin calculates the minimum cut in directed graphs. It
2.46 + /// separates the nodes of the graph into two disjoint sets and the
2.47 + /// summary of the edge capacities go from the first set to the
2.48 + /// second set is the minimum. The algorithm is a modified
2.49 + /// push-relabel preflow algorithm and it calculates the minimum cat
2.50 + /// in \f$ O(n^3) \f$ time. The purpose of such algorithm is testing
2.51 + /// network reliability. For sparse undirected graph you can use the
2.52 + /// algorithm of Nagamochi and Ibraki which solves the undirected
2.53 + /// problem in \f$ O(n^3) \f$ time.
2.54 + ///
2.55 + /// \author Attila Bernath and Balazs Dezso
2.56 + template <typename _Graph,
2.57 + typename _CapacityMap = typename _Graph::template EdgeMap<int>,
2.58 + typename _Tolerance = Tolerance<typename _CapacityMap::Value> >
2.59 + class HaoOrlin {
2.60 + protected:
2.61 +
2.62 + typedef _Graph Graph;
2.63 + typedef _CapacityMap CapacityMap;
2.64 + typedef _Tolerance Tolerance;
2.65 +
2.66 + typedef typename CapacityMap::Value Value;
2.67 +
2.68 +
2.69 + typedef typename Graph::Node Node;
2.70 + typedef typename Graph::NodeIt NodeIt;
2.71 + typedef typename Graph::EdgeIt EdgeIt;
2.72 + typedef typename Graph::OutEdgeIt OutEdgeIt;
2.73 + typedef typename Graph::InEdgeIt InEdgeIt;
2.74 +
2.75 + const Graph* _graph;
2.76 + const CapacityMap* _capacity;
2.77 +
2.78 + typedef typename Graph::template EdgeMap<Value> FlowMap;
2.79 +
2.80 + FlowMap* _preflow;
2.81 +
2.82 + Node _source, _target;
2.83 + int _node_num;
2.84 +
2.85 + typedef ResGraphAdaptor<const Graph, Value, CapacityMap,
2.86 + FlowMap, Tolerance> ResGraph;
2.87 + typedef typename ResGraph::Node ResNode;
2.88 + typedef typename ResGraph::NodeIt ResNodeIt;
2.89 + typedef typename ResGraph::EdgeIt ResEdgeIt;
2.90 + typedef typename ResGraph::OutEdgeIt ResOutEdgeIt;
2.91 + typedef typename ResGraph::Edge ResEdge;
2.92 + typedef typename ResGraph::InEdgeIt ResInEdgeIt;
2.93 +
2.94 + ResGraph* _res_graph;
2.95 +
2.96 + typedef typename Graph::template NodeMap<ResEdge> CurrentArcMap;
2.97 + CurrentArcMap* _current_arc;
2.98 +
2.99 +
2.100 + typedef IterableBoolMap<Graph, Node> WakeMap;
2.101 + WakeMap* _wake;
2.102 +
2.103 + typedef typename Graph::template NodeMap<int> DistMap;
2.104 + DistMap* _dist;
2.105 +
2.106 + typedef typename Graph::template NodeMap<Value> ExcessMap;
2.107 + ExcessMap* _excess;
2.108 +
2.109 + typedef typename Graph::template NodeMap<bool> SourceSetMap;
2.110 + SourceSetMap* _source_set;
2.111 +
2.112 + std::vector<int> _level_size;
2.113 +
2.114 + int _highest_active;
2.115 + std::vector<std::list<Node> > _active_nodes;
2.116 +
2.117 + int _dormant_max;
2.118 + std::vector<std::list<Node> > _dormant;
2.119 +
2.120 +
2.121 + Value _min_cut;
2.122 +
2.123 + typedef typename Graph::template NodeMap<bool> MinCutMap;
2.124 + MinCutMap* _min_cut_map;
2.125 +
2.126 + Tolerance _tolerance;
2.127 +
2.128 + public:
2.129 +
2.130 + HaoOrlin(const Graph& graph, const CapacityMap& capacity,
2.131 + const Tolerance& tolerance = Tolerance()) :
2.132 + _graph(&graph), _capacity(&capacity),
2.133 + _preflow(0), _source(), _target(), _res_graph(0), _current_arc(0),
2.134 + _wake(0),_dist(0), _excess(0), _source_set(0),
2.135 + _highest_active(), _active_nodes(), _dormant_max(), _dormant(),
2.136 + _min_cut(), _min_cut_map(0), _tolerance(tolerance) {}
2.137 +
2.138 + ~HaoOrlin() {
2.139 + if (_res_graph) {
2.140 + delete _res_graph;
2.141 + }
2.142 + if (_min_cut_map) {
2.143 + delete _min_cut_map;
2.144 + }
2.145 + if (_current_arc) {
2.146 + delete _current_arc;
2.147 + }
2.148 + if (_source_set) {
2.149 + delete _source_set;
2.150 + }
2.151 + if (_excess) {
2.152 + delete _excess;
2.153 + }
2.154 + if (_dist) {
2.155 + delete _dist;
2.156 + }
2.157 + if (_wake) {
2.158 + delete _wake;
2.159 + }
2.160 + if (_preflow) {
2.161 + delete _preflow;
2.162 + }
2.163 + }
2.164 +
2.165 + private:
2.166 +
2.167 + void relabel(Node i) {
2.168 + int k = (*_dist)[i];
2.169 + if (_level_size[k] == 1) {
2.170 + ++_dormant_max;
2.171 + for (NodeIt it(*_graph); it != INVALID; ++it) {
2.172 + if ((*_wake)[it] && (*_dist)[it] >= k) {
2.173 + (*_wake)[it] = false;
2.174 + _dormant[_dormant_max].push_front(it);
2.175 + --_level_size[(*_dist)[it]];
2.176 + }
2.177 + }
2.178 + --_highest_active;
2.179 + } else {
2.180 + ResOutEdgeIt e(*_res_graph, i);
2.181 + while (e != INVALID && !(*_wake)[_res_graph->target(e)]) {
2.182 + ++e;
2.183 + }
2.184 +
2.185 + if (e == INVALID){
2.186 + ++_dormant_max;
2.187 + (*_wake)[i] = false;
2.188 + _dormant[_dormant_max].push_front(i);
2.189 + --_level_size[(*_dist)[i]];
2.190 + } else{
2.191 + Node j = _res_graph->target(e);
2.192 + int new_dist = (*_dist)[j];
2.193 + ++e;
2.194 + while (e != INVALID){
2.195 + Node j = _res_graph->target(e);
2.196 + if ((*_wake)[j] && new_dist > (*_dist)[j]) {
2.197 + new_dist = (*_dist)[j];
2.198 + }
2.199 + ++e;
2.200 + }
2.201 + --_level_size[(*_dist)[i]];
2.202 + (*_dist)[i] = new_dist + 1;
2.203 + _highest_active = (*_dist)[i];
2.204 + _active_nodes[_highest_active].push_front(i);
2.205 + ++_level_size[(*_dist)[i]];
2.206 + _res_graph->firstOut((*_current_arc)[i], i);
2.207 + }
2.208 + }
2.209 + }
2.210 +
2.211 + bool selectNewSink(){
2.212 + Node old_target = _target;
2.213 + (*_wake)[_target] = false;
2.214 + --_level_size[(*_dist)[_target]];
2.215 + _dormant[0].push_front(_target);
2.216 + (*_source_set)[_target] = true;
2.217 + if ((int)_dormant[0].size() == _node_num){
2.218 + _dormant[0].clear();
2.219 + return false;
2.220 + }
2.221 +
2.222 + bool wake_was_empty = false;
2.223 +
2.224 + if(_wake->trueNum() == 0) {
2.225 + while (!_dormant[_dormant_max].empty()){
2.226 + (*_wake)[_dormant[_dormant_max].front()] = true;
2.227 + ++_level_size[(*_dist)[_dormant[_dormant_max].front()]];
2.228 + _dormant[_dormant_max].pop_front();
2.229 + }
2.230 + --_dormant_max;
2.231 + wake_was_empty = true;
2.232 + }
2.233 +
2.234 + int min_dist = std::numeric_limits<int>::max();
2.235 + for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
2.236 + if (min_dist > (*_dist)[it]){
2.237 + _target = it;
2.238 + min_dist = (*_dist)[it];
2.239 + }
2.240 + }
2.241 +
2.242 + if (wake_was_empty){
2.243 + for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
2.244 + if (_tolerance.positive((*_excess)[it])) {
2.245 + if ((*_wake)[it] && it != _target) {
2.246 + _active_nodes[(*_dist)[it]].push_front(it);
2.247 + }
2.248 + if (_highest_active < (*_dist)[it]) {
2.249 + _highest_active = (*_dist)[it];
2.250 + }
2.251 + }
2.252 + }
2.253 + }
2.254 +
2.255 + for (ResOutEdgeIt e(*_res_graph, old_target); e!=INVALID; ++e){
2.256 + if (!(*_source_set)[_res_graph->target(e)]){
2.257 + push(e, _res_graph->rescap(e));
2.258 + }
2.259 + }
2.260 +
2.261 + return true;
2.262 + }
2.263 +
2.264 + Node findActiveNode() {
2.265 + while (_highest_active > 0 && _active_nodes[_highest_active].empty()){
2.266 + --_highest_active;
2.267 + }
2.268 + if( _highest_active > 0) {
2.269 + Node n = _active_nodes[_highest_active].front();
2.270 + _active_nodes[_highest_active].pop_front();
2.271 + return n;
2.272 + } else {
2.273 + return INVALID;
2.274 + }
2.275 + }
2.276 +
2.277 + ResEdge findAdmissibleEdge(const Node& n){
2.278 + ResEdge e = (*_current_arc)[n];
2.279 + while (e != INVALID &&
2.280 + ((*_dist)[n] <= (*_dist)[_res_graph->target(e)] ||
2.281 + !(*_wake)[_res_graph->target(e)])) {
2.282 + _res_graph->nextOut(e);
2.283 + }
2.284 + if (e != INVALID) {
2.285 + (*_current_arc)[n] = e;
2.286 + return e;
2.287 + } else {
2.288 + return INVALID;
2.289 + }
2.290 + }
2.291 +
2.292 + void push(ResEdge& e,const Value& delta){
2.293 + _res_graph->augment(e, delta);
2.294 + if (!_tolerance.positive(delta)) return;
2.295 +
2.296 + (*_excess)[_res_graph->source(e)] -= delta;
2.297 + Node a = _res_graph->target(e);
2.298 + if (!_tolerance.positive((*_excess)[a]) && (*_wake)[a] && a != _target) {
2.299 + _active_nodes[(*_dist)[a]].push_front(a);
2.300 + if (_highest_active < (*_dist)[a]) {
2.301 + _highest_active = (*_dist)[a];
2.302 + }
2.303 + }
2.304 + (*_excess)[a] += delta;
2.305 + }
2.306 +
2.307 + Value cutValue() {
2.308 + Value value = 0;
2.309 + for (typename WakeMap::TrueIt it(*_wake); it != INVALID; ++it) {
2.310 + for (InEdgeIt e(*_graph, it); e != INVALID; ++e) {
2.311 + if (!(*_wake)[_graph->source(e)]){
2.312 + value += (*_capacity)[e];
2.313 + }
2.314 + }
2.315 + }
2.316 + return value;
2.317 + }
2.318 +
2.319 + public:
2.320 +
2.321 + /// \brief Initializes the internal data structures.
2.322 + ///
2.323 + /// Initializes the internal data structures. It creates
2.324 + /// the maps, residual graph adaptor and some bucket structures
2.325 + /// for the algorithm.
2.326 + void init() {
2.327 + init(NodeIt(*_graph));
2.328 + }
2.329 +
2.330 + /// \brief Initializes the internal data structures.
2.331 + ///
2.332 + /// Initializes the internal data structures. It creates
2.333 + /// the maps, residual graph adaptor and some bucket structures
2.334 + /// for the algorithm. The \c source node is used as the push-relabel
2.335 + /// algorithm's source.
2.336 + void init(const Node& source) {
2.337 + _source = source;
2.338 + _node_num = countNodes(*_graph);
2.339 +
2.340 + _dormant.resize(_node_num);
2.341 + _level_size.resize(_node_num, 0);
2.342 + _active_nodes.resize(_node_num);
2.343 +
2.344 + if (!_preflow) {
2.345 + _preflow = new FlowMap(*_graph);
2.346 + }
2.347 + if (!_wake) {
2.348 + _wake = new WakeMap(*_graph);
2.349 + }
2.350 + if (!_dist) {
2.351 + _dist = new DistMap(*_graph);
2.352 + }
2.353 + if (!_excess) {
2.354 + _excess = new ExcessMap(*_graph);
2.355 + }
2.356 + if (!_source_set) {
2.357 + _source_set = new SourceSetMap(*_graph);
2.358 + }
2.359 + if (!_current_arc) {
2.360 + _current_arc = new CurrentArcMap(*_graph);
2.361 + }
2.362 + if (!_min_cut_map) {
2.363 + _min_cut_map = new MinCutMap(*_graph);
2.364 + }
2.365 + if (!_res_graph) {
2.366 + _res_graph = new ResGraph(*_graph, *_capacity, *_preflow);
2.367 + }
2.368 +
2.369 + _min_cut = std::numeric_limits<Value>::max();
2.370 + }
2.371 +
2.372 +
2.373 + /// \brief Calculates the minimum cut with the \c source node
2.374 + /// in the first partition.
2.375 + ///
2.376 + /// Calculates the minimum cut with the \c source node
2.377 + /// in the first partition.
2.378 + void calculateOut() {
2.379 + for (NodeIt it(*_graph); it != INVALID; ++it) {
2.380 + if (it != _source) {
2.381 + calculateOut(it);
2.382 + return;
2.383 + }
2.384 + }
2.385 + }
2.386 +
2.387 + /// \brief Calculates the minimum cut with the \c source node
2.388 + /// in the first partition.
2.389 + ///
2.390 + /// Calculates the minimum cut with the \c source node
2.391 + /// in the first partition. The \c target is the initial target
2.392 + /// for the push-relabel algorithm.
2.393 + void calculateOut(const Node& target) {
2.394 + for (NodeIt it(*_graph); it != INVALID; ++it) {
2.395 + (*_wake)[it] = true;
2.396 + (*_dist)[it] = 1;
2.397 + (*_excess)[it] = 0;
2.398 + (*_source_set)[it] = false;
2.399 +
2.400 + _res_graph->firstOut((*_current_arc)[it], it);
2.401 + }
2.402 +
2.403 + _target = target;
2.404 + (*_dist)[target] = 0;
2.405 +
2.406 + for (ResOutEdgeIt it(*_res_graph, _source); it != INVALID; ++it) {
2.407 + push(it, _res_graph->rescap(it));
2.408 + }
2.409 +
2.410 + _dormant[0].push_front(_source);
2.411 + (*_source_set)[_source] = true;
2.412 + _dormant_max = 0;
2.413 + (*_wake)[_source]=false;
2.414 +
2.415 + _level_size[0] = 1;
2.416 + _level_size[1] = _node_num - 1;
2.417 +
2.418 + do {
2.419 + Node n;
2.420 + while ((n = findActiveNode()) != INVALID) {
2.421 + ResEdge e;
2.422 + while (_tolerance.positive((*_excess)[n]) &&
2.423 + (e = findAdmissibleEdge(n)) != INVALID){
2.424 + Value delta;
2.425 + if ((*_excess)[n] < _res_graph->rescap(e)) {
2.426 + delta = (*_excess)[n];
2.427 + } else {
2.428 + delta = _res_graph->rescap(e);
2.429 + _res_graph->nextOut((*_current_arc)[n]);
2.430 + }
2.431 + if (!_tolerance.positive(delta)) continue;
2.432 + _res_graph->augment(e, delta);
2.433 + (*_excess)[_res_graph->source(e)] -= delta;
2.434 + Node a = _res_graph->target(e);
2.435 + if (!_tolerance.positive((*_excess)[a]) && a != _target) {
2.436 + _active_nodes[(*_dist)[a]].push_front(a);
2.437 + }
2.438 + (*_excess)[a] += delta;
2.439 + }
2.440 + if (_tolerance.positive((*_excess)[n])) {
2.441 + relabel(n);
2.442 + }
2.443 + }
2.444 +
2.445 + Value current_value = cutValue();
2.446 + if (_min_cut > current_value){
2.447 + for (NodeIt it(*_graph); it != INVALID; ++it) {
2.448 + _min_cut_map->set(it, !(*_wake)[it]);
2.449 + }
2.450 +
2.451 + _min_cut = current_value;
2.452 + }
2.453 +
2.454 + } while (selectNewSink());
2.455 + }
2.456 +
2.457 + void calculateIn() {
2.458 + for (NodeIt it(*_graph); it != INVALID; ++it) {
2.459 + if (it != _source) {
2.460 + calculateIn(it);
2.461 + return;
2.462 + }
2.463 + }
2.464 + }
2.465 +
2.466 + void run() {
2.467 + init();
2.468 + for (NodeIt it(*_graph); it != INVALID; ++it) {
2.469 + if (it != _source) {
2.470 + startOut(it);
2.471 + // startIn(it);
2.472 + return;
2.473 + }
2.474 + }
2.475 + }
2.476 +
2.477 + void run(const Node& s) {
2.478 + init(s);
2.479 + for (NodeIt it(*_graph); it != INVALID; ++it) {
2.480 + if (it != _source) {
2.481 + startOut(it);
2.482 + // startIn(it);
2.483 + return;
2.484 + }
2.485 + }
2.486 + }
2.487 +
2.488 + void run(const Node& s, const Node& t) {
2.489 + init(s);
2.490 + startOut(t);
2.491 + startIn(t);
2.492 + }
2.493 +
2.494 + /// \brief Returns the value of the minimum value cut with node \c
2.495 + /// source on the source side.
2.496 + ///
2.497 + /// Returns the value of the minimum value cut with node \c source
2.498 + /// on the source side. This function can be called after running
2.499 + /// \ref findMinCut.
2.500 + Value minCut() const {
2.501 + return _min_cut;
2.502 + }
2.503 +
2.504 +
2.505 + /// \brief Returns a minimum value cut.
2.506 + ///
2.507 + /// Sets \c nodeMap to the characteristic vector of a minimum
2.508 + /// value cut with node \c source on the source side. This
2.509 + /// function can be called after running \ref findMinCut.
2.510 + /// \pre nodeMap should be a bool-valued node-map.
2.511 + template <typename NodeMap>
2.512 + Value minCut(NodeMap& nodeMap) const {
2.513 + for (NodeIt it(*_graph); it != INVALID; ++it) {
2.514 + nodeMap.set(it, (*_min_cut_map)[it]);
2.515 + }
2.516 + return minCut();
2.517 + }
2.518 +
2.519 + }; //class HaoOrlin
2.520 +
2.521 +
2.522 +} //namespace lemon
2.523 +
2.524 +#endif //LEMON_HAO_ORLIN_H