Changeset 2530:f86f7e4eb2ba in lemon0.x for lemon/hao_orlin.h
 Timestamp:
 12/04/07 11:55:27 (16 years ago)
 Branch:
 default
 Phase:
 public
 Convert:
 svn:c9d7d8f590d60310b91f818b3a526b0e/lemon/trunk@3407
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/hao_orlin.h
r2391 r2530 20 20 #define LEMON_HAO_ORLIN_H 21 21 22 #include <cassert>23 24 25 26 22 #include <vector> 27 #include <queue>28 23 #include <list> 24 #include <ext/hash_set> 29 25 #include <limits> 30 26 … … 38 34 /// \brief Implementation of the HaoOrlin algorithm. 39 35 /// 40 /// Implementation of the Hao Orlin algorithmsclass for testing network36 /// Implementation of the HaoOrlin algorithm class for testing network 41 37 /// reliability. 42 38 … … 47 43 /// \brief %HaoOrlin algorithm to find a minimum cut in directed graphs. 48 44 /// 49 /// HaoOrlin 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 sourceside (i.e. a set \f$ X\subsetneq V \f$ 53 /// with \f$ source \in X \f$ and minimal outdegree) and in the 54 /// second phase it determines a minimum cut with \f$ source \f$ on the 55 /// sinkside (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \notin X \f$ 56 /// and minimal outdegree). Obviously, the smaller of these two 57 /// cuts will be a minimum cut of \f$ D \f$. The algorithm is a 58 /// modified pushrelabel preflow algorithm and our implementation 59 /// calculates the minimum cut in \f$ O(n^3) \f$ time (we use the 60 /// highestlabel 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 65 /// algorithm 66 /// class. 45 /// HaoOrlin 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 47 /// consists of two phases: in the first phase it determines a 48 /// minimum cut with \f$ source \f$ on the sourceside (i.e. a set 49 /// \f$ X\subsetneq V \f$ with \f$ source \in X \f$ and minimal 50 /// outdegree) and in the second phase it determines a minimum cut 51 /// with \f$ source \f$ on the sinkside (i.e. a set 52 /// \f$ X\subsetneq V \f$ with \f$ source \notin X \f$ and minimal 53 /// outdegree). Obviously, the smaller of these two cuts will be a 54 /// minimum cut of \f$ D \f$. The algorithm is a modified 55 /// pushrelabel preflow algorithm and our implementation calculates 56 /// the minimum cut in \f$ O(n^2\sqrt{m}) \f$ time (we use the 57 /// highestlabel rule), or in \f$O(nm)\f$ for unit capacities. The 58 /// purpose of such algorithm is testing network reliability. For an 59 /// undirected graph you can run just the first phase of the 60 /// algorithm or you can use the algorithm of Nagamochi and Ibaraki 61 /// which solves the undirected problem in 62 /// \f$ O(nm + n^2 \log(n)) \f$ time: it is implemented in the 63 /// NagamochiIbaraki algorithm class. 67 64 /// 68 65 /// \param _Graph is the graph type of the algorithm. … … 81 78 #endif 82 79 class HaoOrlin { 83 pr otected:80 private: 84 81 85 82 typedef _Graph Graph; … … 89 86 typedef typename CapacityMap::Value Value; 90 87 88 GRAPH_TYPEDEFS(typename Graph); 91 89 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; 97 98 const Graph* _graph; 99 90 const Graph& _graph; 100 91 const CapacityMap* _capacity; 101 92 102 93 typedef typename Graph::template EdgeMap<Value> FlowMap; 103 104 FlowMap* _preflow; 105 106 Node _source, _target; 94 FlowMap* _flow; 95 96 Node _source; 97 107 98 int _node_num; 108 99 109 typedef ResGraphAdaptor<const Graph, Value, CapacityMap, 110 FlowMap, Tolerance> OutResGraph; 111 typedef typename OutResGraph::Edge OutResEdge; 100 // Bucketing structure 101 std::vector<Node> _first, _last; 102 typename Graph::template NodeMap<Node>* _next; 103 typename Graph::template NodeMap<Node>* _prev; 104 typename Graph::template NodeMap<bool>* _active; 105 typename Graph::template NodeMap<int>* _bucket; 112 106 113 OutResGraph* _out_res_graph; 114 115 typedef RevGraphAdaptor<const Graph> RevGraph; 116 RevGraph* _rev_graph; 117 118 typedef ResGraphAdaptor<const RevGraph, Value, CapacityMap, 119 FlowMap, Tolerance> InResGraph; 120 typedef typename InResGraph::Edge InResEdge; 121 122 InResGraph* _in_res_graph; 123 124 typedef IterableBoolMap<Graph, Node> WakeMap; 125 WakeMap* _wake; 126 127 typedef typename Graph::template NodeMap<int> DistMap; 128 DistMap* _dist; 107 std::vector<bool> _dormant; 108 109 std::list<std::list<int> > _sets; 110 std::list<int>::iterator _highest; 129 111 130 112 typedef typename Graph::template NodeMap<Value> ExcessMap; … … 133 115 typedef typename Graph::template NodeMap<bool> SourceSetMap; 134 116 SourceSetMap* _source_set; 135 136 std::vector<int> _level_size;137 138 int _highest_active;139 std::vector<std::list<Node> > _active_nodes;140 141 int _dormant_max;142 std::vector<std::list<Node> > _dormant;143 144 117 145 118 Value _min_cut; … … 157 130 HaoOrlin(const Graph& graph, const CapacityMap& capacity, 158 131 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) {} 132 _graph(graph), _capacity(&capacity), _flow(0), _source(), 133 _node_num(), _first(), _last(), _next(0), _prev(0), 134 _active(0), _bucket(0), _dormant(), _sets(), _highest(), 135 _excess(0), _source_set(0), _min_cut(), _min_cut_map(0), 136 _tolerance(tolerance) {} 165 137 166 138 ~HaoOrlin() { … … 168 140 delete _min_cut_map; 169 141 } 170 if (_in_res_graph) {171 delete _in_res_graph;172 }173 if (_rev_graph) {174 delete _rev_graph;175 }176 if (_out_res_graph) {177 delete _out_res_graph;178 }179 142 if (_source_set) { 180 143 delete _source_set; … … 183 146 delete _excess; 184 147 } 185 if (_dist) { 186 delete _dist; 187 } 188 if (_wake) { 189 delete _wake; 190 } 191 if (_preflow) { 192 delete _preflow; 148 if (_next) { 149 delete _next; 150 } 151 if (_prev) { 152 delete _prev; 153 } 154 if (_active) { 155 delete _active; 156 } 157 if (_bucket) { 158 delete _bucket; 159 } 160 if (_flow) { 161 delete _flow; 193 162 } 194 163 } 195 164 196 165 private: 166 167 void activate(const Node& i) { 168 _active>set(i, true); 169 170 int bucket = (*_bucket)[i]; 171 172 if ((*_prev)[i] == INVALID  (*_active)[(*_prev)[i]]) return; 173 //unlace 174 _next>set((*_prev)[i], (*_next)[i]); 175 if ((*_next)[i] != INVALID) { 176 _prev>set((*_next)[i], (*_prev)[i]); 177 } else { 178 _last[bucket] = (*_prev)[i]; 179 } 180 //lace 181 _next>set(i, _first[bucket]); 182 _prev>set(_first[bucket], i); 183 _prev>set(i, INVALID); 184 _first[bucket] = i; 185 } 186 187 void deactivate(const Node& i) { 188 _active>set(i, false); 189 int bucket = (*_bucket)[i]; 190 191 if ((*_next)[i] == INVALID  !(*_active)[(*_next)[i]]) return; 192 193 //unlace 194 _prev>set((*_next)[i], (*_prev)[i]); 195 if ((*_prev)[i] != INVALID) { 196 _next>set((*_prev)[i], (*_next)[i]); 197 } else { 198 _first[bucket] = (*_next)[i]; 199 } 200 //lace 201 _prev>set(i, _last[bucket]); 202 _next>set(_last[bucket], i); 203 _next>set(i, INVALID); 204 _last[bucket] = i; 205 } 206 207 void addItem(const Node& i, int bucket) { 208 (*_bucket)[i] = bucket; 209 if (_last[bucket] != INVALID) { 210 _prev>set(i, _last[bucket]); 211 _next>set(_last[bucket], i); 212 _next>set(i, INVALID); 213 _last[bucket] = i; 214 } else { 215 _prev>set(i, INVALID); 216 _first[bucket] = i; 217 _next>set(i, INVALID); 218 _last[bucket] = i; 219 } 220 } 197 221 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; 202 203 for (typename Graph::EdgeIt it(*_graph); it != INVALID; ++it) { 204 (*_preflow)[it] = 0; 205 } 206 for (NodeIt it(*_graph); it != INVALID; ++it) { 207 (*_wake)[it] = true; 208 (*_dist)[it] = 1; 209 (*_excess)[it] = 0; 210 (*_source_set)[it] = false; 211 } 212 213 _dormant[0].push_front(_source); 214 (*_source_set)[_source] = true; 215 _dormant_max = 0; 216 (*_wake)[_source] = false; 217 218 _level_size[0] = 1; 219 _level_size[1] = _node_num  1; 220 221 _target = target; 222 (*_dist)[target] = 0; 223 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]; 233 } 234 } 235 (*_excess)[a] += delta; 236 } 237 238 239 do { 240 Node n; 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; 222 void findMinCutOut() { 223 224 for (NodeIt n(_graph); n != INVALID; ++n) { 225 _excess>set(n, 0); 226 } 227 228 for (EdgeIt e(_graph); e != INVALID; ++e) { 229 _flow>set(e, 0); 230 } 231 232 int bucket_num = 1; 233 234 { 235 typename Graph::template NodeMap<bool> reached(_graph, false); 236 237 reached.set(_source, true); 238 239 bool first_set = true; 240 241 for (NodeIt t(_graph); t != INVALID; ++t) { 242 if (reached[t]) continue; 243 _sets.push_front(std::list<int>()); 244 _sets.front().push_front(bucket_num); 245 _dormant[bucket_num] = !first_set; 246 247 _bucket>set(t, bucket_num); 248 _first[bucket_num] = _last[bucket_num] = t; 249 _next>set(t, INVALID); 250 _prev>set(t, INVALID); 251 252 ++bucket_num; 253 254 std::vector<Node> queue; 255 queue.push_back(t); 256 reached.set(t, true); 257 258 while (!queue.empty()) { 259 _sets.front().push_front(bucket_num); 260 _dormant[bucket_num] = !first_set; 261 _first[bucket_num] = _last[bucket_num] = INVALID; 262 263 std::vector<Node> nqueue; 264 for (int i = 0; i < int(queue.size()); ++i) { 265 Node n = queue[i]; 266 for (InEdgeIt e(_graph, n); e != INVALID; ++e) { 267 Node u = _graph.source(e); 268 if (!reached[u] && _tolerance.positive((*_capacity)[e])) { 269 reached.set(u, true); 270 addItem(u, bucket_num); 271 nqueue.push_back(u); 272 } 273 } 274 } 275 queue.swap(nqueue); 276 ++bucket_num; 277 } 278 _sets.front().pop_front(); 279 bucket_num; 280 first_set = false; 281 } 282 283 _bucket>set(_source, 0); 284 _dormant[0] = true; 285 } 286 _source_set>set(_source, true); 287 288 Node target = _last[_sets.back().back()]; 289 { 290 for (OutEdgeIt e(_graph, _source); e != INVALID; ++e) { 291 if (_tolerance.positive((*_capacity)[e])) { 292 Node u = _graph.target(e); 293 _flow>set(e, (*_capacity)[e]); 294 _excess>set(u, (*_excess)[u] + (*_capacity)[e]); 295 if (!(*_active)[u] && u != _source) { 296 activate(u); 297 } 298 } 299 } 300 if ((*_active)[target]) { 301 deactivate(target); 302 } 303 304 _highest = _sets.back().begin(); 305 while (_highest != _sets.back().end() && 306 !(*_active)[_first[*_highest]]) { 307 ++_highest; 308 } 309 } 310 311 312 while (true) { 313 while (_highest != _sets.back().end()) { 314 Node n = _first[*_highest]; 315 Value excess = (*_excess)[n]; 316 int next_bucket = _node_num; 317 318 int under_bucket; 319 if (++std::list<int>::iterator(_highest) == _sets.back().end()) { 320 under_bucket = 1; 321 } else { 322 under_bucket = *(++std::list<int>::iterator(_highest)); 323 } 324 325 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) { 326 Node v = _graph.target(e); 327 if (_dormant[(*_bucket)[v]]) continue; 328 Value rem = (*_capacity)[e]  (*_flow)[e]; 329 if (!_tolerance.positive(rem)) continue; 330 if ((*_bucket)[v] == under_bucket) { 331 if (!(*_active)[v] && v != target) { 332 activate(v); 333 } 334 if (!_tolerance.less(rem, excess)) { 335 _flow>set(e, (*_flow)[e] + excess); 336 _excess>set(v, (*_excess)[v] + excess); 337 excess = 0; 338 goto no_more_push; 339 } else { 340 excess = rem; 341 _excess>set(v, (*_excess)[v] + rem); 342 _flow>set(e, (*_capacity)[e]); 343 } 344 } else if (next_bucket > (*_bucket)[v]) { 345 next_bucket = (*_bucket)[v]; 346 } 347 } 348 349 for (InEdgeIt e(_graph, n); e != INVALID; ++e) { 350 Node v = _graph.source(e); 351 if (_dormant[(*_bucket)[v]]) continue; 352 Value rem = (*_flow)[e]; 353 if (!_tolerance.positive(rem)) continue; 354 if ((*_bucket)[v] == under_bucket) { 355 if (!(*_active)[v] && v != target) { 356 activate(v); 357 } 358 if (!_tolerance.less(rem, excess)) { 359 _flow>set(e, (*_flow)[e]  excess); 360 _excess>set(v, (*_excess)[v] + excess); 361 excess = 0; 362 goto no_more_push; 363 } else { 364 excess = rem; 365 _excess>set(v, (*_excess)[v] + rem); 366 _flow>set(e, 0); 367 } 368 } else if (next_bucket > (*_bucket)[v]) { 369 next_bucket = (*_bucket)[v]; 370 } 371 } 372 373 no_more_push: 374 375 _excess>set(n, excess); 376 377 if (excess != 0) { 378 if ((*_next)[n] == INVALID) { 379 typename std::list<std::list<int> >::iterator new_set = 380 _sets.insert(_sets.end(), std::list<int>()); 381 new_set>splice(new_set>end(), _sets.back(), 382 _sets.back().begin(), ++_highest); 383 for (std::list<int>::iterator it = new_set>begin(); 384 it != new_set>end(); ++it) { 385 _dormant[*it] = true; 386 } 387 while (_highest != _sets.back().end() && 388 !(*_active)[_first[*_highest]]) { 389 ++_highest; 390 } 391 } else if (next_bucket == _node_num) { 392 _first[(*_bucket)[n]] = (*_next)[n]; 393 _prev>set((*_next)[n], INVALID); 394 395 std::list<std::list<int> >::iterator new_set = 396 _sets.insert(_sets.end(), std::list<int>()); 397 398 new_set>push_front(bucket_num); 399 _bucket>set(n, bucket_num); 400 _first[bucket_num] = _last[bucket_num] = n; 401 _next>set(n, INVALID); 402 _prev>set(n, INVALID); 403 _dormant[bucket_num] = true; 404 ++bucket_num; 405 406 while (_highest != _sets.back().end() && 407 !(*_active)[_first[*_highest]]) { 408 ++_highest; 409 } 248 410 } else { 249 delta = (*_excess)[n]; 250 (*_excess)[n] = 0; 251 } 252 res_graph.augment(e, delta); 253 if ((*_excess)[a] == 0 && a != _target) { 254 _active_nodes[(*_dist)[a]].push_front(a); 255 } 256 (*_excess)[a] += delta; 257 if ((*_excess)[n] == 0) break; 258 } 259 if ((*_excess)[n] != 0) { 260 relabel(n, res_graph); 261 } 262 } 263 264 Value current_value = cutValue(out); 265 if (_min_cut > current_value){ 266 if (out) { 267 for (NodeIt it(*_graph); it != INVALID; ++it) { 268 _min_cut_map>set(it, !(*_wake)[it]); 269 } 270 } else { 271 for (NodeIt it(*_graph); it != INVALID; ++it) { 272 _min_cut_map>set(it, (*_wake)[it]); 273 } 274 } 275 276 _min_cut = current_value; 277 } 278 279 } while (selectNewSink(res_graph)); 280 } 281 282 template <typename ResGraph> 283 void relabel(const Node& n, ResGraph& res_graph) { 284 typedef typename ResGraph::OutEdgeIt ResOutEdgeIt; 285 286 int k = (*_dist)[n]; 287 if (_level_size[k] == 1) { 288 ++_dormant_max; 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]]; 294 } 295 } 296 _highest_active; 297 } else { 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]; 303 } 304 } 305 if (new_dist == _node_num) { 306 ++_dormant_max; 307 (*_wake)[n] = false; 308 _dormant[_dormant_max].push_front(n); 309 _level_size[(*_dist)[n]]; 310 } else { 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]]; 316 } 317 } 318 } 319 320 template <typename ResGraph> 321 bool selectNewSink(ResGraph& res_graph) { 322 typedef typename ResGraph::OutEdgeIt ResOutEdgeIt; 323 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){ 330 _dormant[0].clear(); 331 return false; 332 } 333 334 bool wake_was_empty = false; 335 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(); 341 } 342 _dormant_max; 343 wake_was_empty = true; 344 } 345 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]){ 349 _target = it; 350 min_dist = (*_dist)[it]; 351 } 352 } 353 354 if (wake_was_empty){ 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]; 360 } 361 } 362 } 363 } 364 365 Node n = old_target; 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]; 376 } 377 } 378 (*_excess)[a] += delta; 379 } 411 _first[*_highest] = (*_next)[n]; 412 _prev>set((*_next)[n], INVALID); 413 414 while (next_bucket != *_highest) { 415 _highest; 416 } 417 418 if (_highest == _sets.back().begin()) { 419 _sets.back().push_front(bucket_num); 420 _dormant[bucket_num] = false; 421 _first[bucket_num] = _last[bucket_num] = INVALID; 422 ++bucket_num; 423 } 424 _highest; 425 426 _bucket>set(n, *_highest); 427 _next>set(n, _first[*_highest]); 428 if (_first[*_highest] != INVALID) { 429 _prev>set(_first[*_highest], n); 430 } else { 431 _last[*_highest] = n; 432 } 433 _first[*_highest] = n; 434 } 435 } else { 436 437 deactivate(n); 438 if (!(*_active)[_first[*_highest]]) { 439 ++_highest; 440 if (_highest != _sets.back().end() && 441 !(*_active)[_first[*_highest]]) { 442 _highest = _sets.back().end(); 443 } 444 } 445 } 446 } 447 448 if ((*_excess)[target] < _min_cut) { 449 _min_cut = (*_excess)[target]; 450 for (NodeIt i(_graph); i != INVALID; ++i) { 451 _min_cut_map>set(i, true); 452 } 453 for (std::list<int>::iterator it = _sets.back().begin(); 454 it != _sets.back().end(); ++it) { 455 Node n = _first[*it]; 456 while (n != INVALID) { 457 _min_cut_map>set(n, false); 458 n = (*_next)[n]; 459 } 460 } 461 } 462 463 { 464 Node new_target; 465 if ((*_prev)[target] != INVALID) { 466 _last[(*_bucket)[target]] = (*_prev)[target]; 467 _next>set((*_prev)[target], INVALID); 468 new_target = (*_prev)[target]; 469 } else { 470 _sets.back().pop_back(); 471 if (_sets.back().empty()) { 472 _sets.pop_back(); 473 if (_sets.empty()) 474 break; 475 for (std::list<int>::iterator it = _sets.back().begin(); 476 it != _sets.back().end(); ++it) { 477 _dormant[*it] = false; 478 } 479 } 480 new_target = _last[_sets.back().back()]; 481 } 482 483 _bucket>set(target, 0); 484 485 _source_set>set(target, true); 486 for (OutEdgeIt e(_graph, target); e != INVALID; ++e) { 487 Value rem = (*_capacity)[e]  (*_flow)[e]; 488 if (!_tolerance.positive(rem)) continue; 489 Node v = _graph.target(e); 490 if (!(*_active)[v] && !(*_source_set)[v]) { 491 activate(v); 492 } 493 _excess>set(v, (*_excess)[v] + rem); 494 _flow>set(e, (*_capacity)[e]); 495 } 496 497 for (InEdgeIt e(_graph, target); e != INVALID; ++e) { 498 Value rem = (*_flow)[e]; 499 if (!_tolerance.positive(rem)) continue; 500 Node v = _graph.source(e); 501 if (!(*_active)[v] && !(*_source_set)[v]) { 502 activate(v); 503 } 504 _excess>set(v, (*_excess)[v] + rem); 505 _flow>set(e, 0); 506 } 507 508 target = new_target; 509 if ((*_active)[target]) { 510 deactivate(target); 511 } 512 513 _highest = _sets.back().begin(); 514 while (_highest != _sets.back().end() && 515 !(*_active)[_first[*_highest]]) { 516 ++_highest; 517 } 518 } 519 } 520 } 521 522 void findMinCutIn() { 523 524 for (NodeIt n(_graph); n != INVALID; ++n) { 525 _excess>set(n, 0); 526 } 527 528 for (EdgeIt e(_graph); e != INVALID; ++e) { 529 _flow>set(e, 0); 530 } 531 532 int bucket_num = 1; 380 533 381 return true; 382 } 383 384 Node findActiveNode() { 385 while (_highest_active > 0 && _active_nodes[_highest_active].empty()){ 386 _highest_active; 387 } 388 if( _highest_active > 0) { 389 Node n = _active_nodes[_highest_active].front(); 390 _active_nodes[_highest_active].pop_front(); 391 return n; 392 } else { 393 return INVALID; 394 } 395 } 396 397 Value cutValue(bool out) { 398 Value value = 0; 399 if (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]; 404 } 405 } 406 } 407 } else { 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]; 412 } 413 } 414 } 415 } 416 return value; 417 } 418 534 { 535 typename Graph::template NodeMap<bool> reached(_graph, false); 536 537 reached.set(_source, true); 538 539 bool first_set = true; 540 541 for (NodeIt t(_graph); t != INVALID; ++t) { 542 if (reached[t]) continue; 543 _sets.push_front(std::list<int>()); 544 _sets.front().push_front(bucket_num); 545 _dormant[bucket_num] = !first_set; 546 547 _bucket>set(t, bucket_num); 548 _first[bucket_num] = _last[bucket_num] = t; 549 _next>set(t, INVALID); 550 _prev>set(t, INVALID); 551 552 ++bucket_num; 553 554 std::vector<Node> queue; 555 queue.push_back(t); 556 reached.set(t, true); 557 558 while (!queue.empty()) { 559 _sets.front().push_front(bucket_num); 560 _dormant[bucket_num] = !first_set; 561 _first[bucket_num] = _last[bucket_num] = INVALID; 562 563 std::vector<Node> nqueue; 564 for (int i = 0; i < int(queue.size()); ++i) { 565 Node n = queue[i]; 566 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) { 567 Node u = _graph.target(e); 568 if (!reached[u] && _tolerance.positive((*_capacity)[e])) { 569 reached.set(u, true); 570 addItem(u, bucket_num); 571 nqueue.push_back(u); 572 } 573 } 574 } 575 queue.swap(nqueue); 576 ++bucket_num; 577 } 578 _sets.front().pop_front(); 579 bucket_num; 580 first_set = false; 581 } 582 583 _bucket>set(_source, 0); 584 _dormant[0] = true; 585 } 586 _source_set>set(_source, true); 587 588 Node target = _last[_sets.back().back()]; 589 { 590 for (InEdgeIt e(_graph, _source); e != INVALID; ++e) { 591 if (_tolerance.positive((*_capacity)[e])) { 592 Node u = _graph.source(e); 593 _flow>set(e, (*_capacity)[e]); 594 _excess>set(u, (*_excess)[u] + (*_capacity)[e]); 595 if (!(*_active)[u] && u != _source) { 596 activate(u); 597 } 598 } 599 } 600 if ((*_active)[target]) { 601 deactivate(target); 602 } 603 604 _highest = _sets.back().begin(); 605 while (_highest != _sets.back().end() && 606 !(*_active)[_first[*_highest]]) { 607 ++_highest; 608 } 609 } 610 611 612 while (true) { 613 while (_highest != _sets.back().end()) { 614 Node n = _first[*_highest]; 615 Value excess = (*_excess)[n]; 616 int next_bucket = _node_num; 617 618 int under_bucket; 619 if (++std::list<int>::iterator(_highest) == _sets.back().end()) { 620 under_bucket = 1; 621 } else { 622 under_bucket = *(++std::list<int>::iterator(_highest)); 623 } 624 625 for (InEdgeIt e(_graph, n); e != INVALID; ++e) { 626 Node v = _graph.source(e); 627 if (_dormant[(*_bucket)[v]]) continue; 628 Value rem = (*_capacity)[e]  (*_flow)[e]; 629 if (!_tolerance.positive(rem)) continue; 630 if ((*_bucket)[v] == under_bucket) { 631 if (!(*_active)[v] && v != target) { 632 activate(v); 633 } 634 if (!_tolerance.less(rem, excess)) { 635 _flow>set(e, (*_flow)[e] + excess); 636 _excess>set(v, (*_excess)[v] + excess); 637 excess = 0; 638 goto no_more_push; 639 } else { 640 excess = rem; 641 _excess>set(v, (*_excess)[v] + rem); 642 _flow>set(e, (*_capacity)[e]); 643 } 644 } else if (next_bucket > (*_bucket)[v]) { 645 next_bucket = (*_bucket)[v]; 646 } 647 } 648 649 for (OutEdgeIt e(_graph, n); e != INVALID; ++e) { 650 Node v = _graph.target(e); 651 if (_dormant[(*_bucket)[v]]) continue; 652 Value rem = (*_flow)[e]; 653 if (!_tolerance.positive(rem)) continue; 654 if ((*_bucket)[v] == under_bucket) { 655 if (!(*_active)[v] && v != target) { 656 activate(v); 657 } 658 if (!_tolerance.less(rem, excess)) { 659 _flow>set(e, (*_flow)[e]  excess); 660 _excess>set(v, (*_excess)[v] + excess); 661 excess = 0; 662 goto no_more_push; 663 } else { 664 excess = rem; 665 _excess>set(v, (*_excess)[v] + rem); 666 _flow>set(e, 0); 667 } 668 } else if (next_bucket > (*_bucket)[v]) { 669 next_bucket = (*_bucket)[v]; 670 } 671 } 672 673 no_more_push: 674 675 _excess>set(n, excess); 676 677 if (excess != 0) { 678 if ((*_next)[n] == INVALID) { 679 typename std::list<std::list<int> >::iterator new_set = 680 _sets.insert(_sets.end(), std::list<int>()); 681 new_set>splice(new_set>end(), _sets.back(), 682 _sets.back().begin(), ++_highest); 683 for (std::list<int>::iterator it = new_set>begin(); 684 it != new_set>end(); ++it) { 685 _dormant[*it] = true; 686 } 687 while (_highest != _sets.back().end() && 688 !(*_active)[_first[*_highest]]) { 689 ++_highest; 690 } 691 } else if (next_bucket == _node_num) { 692 _first[(*_bucket)[n]] = (*_next)[n]; 693 _prev>set((*_next)[n], INVALID); 694 695 std::list<std::list<int> >::iterator new_set = 696 _sets.insert(_sets.end(), std::list<int>()); 697 698 new_set>push_front(bucket_num); 699 _bucket>set(n, bucket_num); 700 _first[bucket_num] = _last[bucket_num] = n; 701 _next>set(n, INVALID); 702 _prev>set(n, INVALID); 703 _dormant[bucket_num] = true; 704 ++bucket_num; 705 706 while (_highest != _sets.back().end() && 707 !(*_active)[_first[*_highest]]) { 708 ++_highest; 709 } 710 } else { 711 _first[*_highest] = (*_next)[n]; 712 _prev>set((*_next)[n], INVALID); 713 714 while (next_bucket != *_highest) { 715 _highest; 716 } 717 if (_highest == _sets.back().begin()) { 718 _sets.back().push_front(bucket_num); 719 _dormant[bucket_num] = false; 720 _first[bucket_num] = _last[bucket_num] = INVALID; 721 ++bucket_num; 722 } 723 _highest; 724 725 _bucket>set(n, *_highest); 726 _next>set(n, _first[*_highest]); 727 if (_first[*_highest] != INVALID) { 728 _prev>set(_first[*_highest], n); 729 } else { 730 _last[*_highest] = n; 731 } 732 _first[*_highest] = n; 733 } 734 } else { 735 736 deactivate(n); 737 if (!(*_active)[_first[*_highest]]) { 738 ++_highest; 739 if (_highest != _sets.back().end() && 740 !(*_active)[_first[*_highest]]) { 741 _highest = _sets.back().end(); 742 } 743 } 744 } 745 } 746 747 if ((*_excess)[target] < _min_cut) { 748 _min_cut = (*_excess)[target]; 749 for (NodeIt i(_graph); i != INVALID; ++i) { 750 _min_cut_map>set(i, false); 751 } 752 for (std::list<int>::iterator it = _sets.back().begin(); 753 it != _sets.back().end(); ++it) { 754 Node n = _first[*it]; 755 while (n != INVALID) { 756 _min_cut_map>set(n, true); 757 n = (*_next)[n]; 758 } 759 } 760 } 761 762 { 763 Node new_target; 764 if ((*_prev)[target] != INVALID) { 765 _last[(*_bucket)[target]] = (*_prev)[target]; 766 _next>set((*_prev)[target], INVALID); 767 new_target = (*_prev)[target]; 768 } else { 769 _sets.back().pop_back(); 770 if (_sets.back().empty()) { 771 _sets.pop_back(); 772 if (_sets.empty()) 773 break; 774 for (std::list<int>::iterator it = _sets.back().begin(); 775 it != _sets.back().end(); ++it) { 776 _dormant[*it] = false; 777 } 778 } 779 new_target = _last[_sets.back().back()]; 780 } 781 782 _bucket>set(target, 0); 783 784 _source_set>set(target, true); 785 for (InEdgeIt e(_graph, target); e != INVALID; ++e) { 786 Value rem = (*_capacity)[e]  (*_flow)[e]; 787 if (!_tolerance.positive(rem)) continue; 788 Node v = _graph.source(e); 789 if (!(*_active)[v] && !(*_source_set)[v]) { 790 activate(v); 791 } 792 _excess>set(v, (*_excess)[v] + rem); 793 _flow>set(e, (*_capacity)[e]); 794 } 795 796 for (OutEdgeIt e(_graph, target); e != INVALID; ++e) { 797 Value rem = (*_flow)[e]; 798 if (!_tolerance.positive(rem)) continue; 799 Node v = _graph.target(e); 800 if (!(*_active)[v] && !(*_source_set)[v]) { 801 activate(v); 802 } 803 _excess>set(v, (*_excess)[v] + rem); 804 _flow>set(e, 0); 805 } 806 807 target = new_target; 808 if ((*_active)[target]) { 809 deactivate(target); 810 } 811 812 _highest = _sets.back().begin(); 813 while (_highest != _sets.back().end() && 814 !(*_active)[_first[*_highest]]) { 815 ++_highest; 816 } 817 } 818 } 819 } 419 820 420 821 public: … … 436 837 /// for the algorithm. 437 838 void init() { 438 init(NodeIt( *_graph));839 init(NodeIt(_graph)); 439 840 } 440 841 … … 447 848 void init(const Node& source) { 448 849 _source = source; 449 _node_num = countNodes(*_graph); 850 851 _node_num = countNodes(_graph); 852 853 _first.resize(_node_num); 854 _last.resize(_node_num); 450 855 451 856 _dormant.resize(_node_num); 452 _level_size.resize(_node_num, 0); 453 _active_nodes.resize(_node_num); 454 455 if (!_preflow) { 456 _preflow = new FlowMap(*_graph); 457 } 458 if (!_wake) { 459 _wake = new WakeMap(*_graph); 460 } 461 if (!_dist) { 462 _dist = new DistMap(*_graph); 857 858 if (!_flow) { 859 _flow = new FlowMap(_graph); 860 } 861 if (!_next) { 862 _next = new typename Graph::template NodeMap<Node>(_graph); 863 } 864 if (!_prev) { 865 _prev = new typename Graph::template NodeMap<Node>(_graph); 866 } 867 if (!_active) { 868 _active = new typename Graph::template NodeMap<bool>(_graph); 869 } 870 if (!_bucket) { 871 _bucket = new typename Graph::template NodeMap<int>(_graph); 463 872 } 464 873 if (!_excess) { 465 _excess = new ExcessMap(*_graph);874 _excess = new ExcessMap(_graph); 466 875 } 467 876 if (!_source_set) { 468 _source_set = new SourceSetMap(*_graph); 469 } 470 if (!_out_res_graph) { 471 _out_res_graph = new OutResGraph(*_graph, *_capacity, *_preflow); 472 } 473 if (!_rev_graph) { 474 _rev_graph = new RevGraph(*_graph); 475 } 476 if (!_in_res_graph) { 477 _in_res_graph = new InResGraph(*_rev_graph, *_capacity, *_preflow); 877 _source_set = new SourceSetMap(_graph); 478 878 } 479 879 if (!_min_cut_map) { 480 _min_cut_map = new MinCutMap(*_graph);880 _min_cut_map = new MinCutMap(_graph); 481 881 } 482 882 … … 488 888 /// sourceside. 489 889 /// 890 /// Calculates a minimum cut with \f$ source \f$ on the 891 /// sourceside (i.e. a set \f$ X\subsetneq V \f$ with \f$ source 892 /// \in X \f$ and minimal outdegree). 893 void calculateOut() { 894 findMinCutOut(); 895 } 896 490 897 /// \brief Calculates a minimum cut with \f$ source \f$ on the 491 /// sourceside (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$ 492 /// and minimal outdegree). 493 void calculateOut() { 494 for (NodeIt it(*_graph); it != INVALID; ++it) { 495 if (it != _source) { 496 calculateOut(it); 497 return; 498 } 499 } 500 } 501 502 /// \brief Calculates a minimum cut with \f$ source \f$ on the 503 /// sourceside. 898 /// targetside. 504 899 /// 505 /// \brief Calculates a minimum cut with \f$ source \f$ on the 506 /// sourceside (i.e. a set \f$ X\subsetneq V \f$ with \f$ source \in X \f$ 507 /// and minimal outdegree). The \c target is the initial target 508 /// for the pushrelabel algorithm. 509 void calculateOut(const Node& target) { 510 findMinCut(target, true, *_out_res_graph); 511 } 512 513 /// \brief Calculates a minimum cut with \f$ source \f$ on the 514 /// sinkside. 515 /// 516 /// \brief Calculates a minimum cut with \f$ source \f$ on the 517 /// sinkside (i.e. a set \f$ X\subsetneq V \f$ with 518 /// \f$ source \notin X \f$ 519 /// and minimal outdegree). 900 /// Calculates a minimum cut with \f$ source \f$ on the 901 /// targetside (i.e. a set \f$ X\subsetneq V \f$ with \f$ source 902 /// \in X \f$ and minimal outdegree). 520 903 void calculateIn() { 521 for (NodeIt it(*_graph); it != INVALID; ++it) { 522 if (it != _source) { 523 calculateIn(it); 524 return; 525 } 526 } 527 } 528 529 /// \brief Calculates a minimum cut with \f$ source \f$ on the 530 /// sinkside. 531 /// 532 /// \brief Calculates a minimum cut with \f$ source \f$ on the 533 /// sinkside (i.e. a set \f$ X\subsetneq V 534 /// \f$ with \f$ source \notin X \f$ and minimal outdegree). 535 /// The \c target is the initial 536 /// target for the pushrelabel algorithm. 537 void calculateIn(const Node& target) { 538 findMinCut(target, false, *_in_res_graph); 539 } 904 findMinCutIn(); 905 } 906 540 907 541 908 /// \brief Runs the algorithm. … … 546 913 void run() { 547 914 init(); 548 for (NodeIt it(*_graph); it != INVALID; ++it) { 549 if (it != _source) { 550 calculateOut(it); 551 calculateIn(it); 552 return; 553 } 554 } 915 calculateOut(); 916 calculateIn(); 555 917 } 556 918 … … 562 924 void run(const Node& s) { 563 925 init(s); 564 for (NodeIt it(*_graph); it != INVALID; ++it) { 565 if (it != _source) { 566 calculateOut(it); 567 calculateIn(it); 568 return; 569 } 570 } 571 } 572 573 /// \brief Runs the algorithm. 574 /// 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) { 578 init(s); 579 calculateOut(t); 580 calculateIn(t); 926 calculateOut(); 927 calculateIn(); 581 928 } 582 929 … … 609 956 template <typename NodeMap> 610 957 Value minCut(NodeMap& nodeMap) const { 611 for (NodeIt it( *_graph); it != INVALID; ++it) {958 for (NodeIt it(_graph); it != INVALID; ++it) { 612 959 nodeMap.set(it, (*_min_cut_map)[it]); 613 960 }
Note: See TracChangeset
for help on using the changeset viewer.