gravatar
deba@inf.elte.hu
deba@inf.elte.hu
New queue implementation for HaoOrlin class (#58)
0 1 0
default
1 file changed with 47 insertions and 66 deletions:
↑ Collapse diff ↑
Ignore white space 1536 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2008
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_HAO_ORLIN_H
20 20
#define LEMON_HAO_ORLIN_H
21 21

	
22 22
#include <vector>
23 23
#include <list>
24 24
#include <limits>
25 25

	
26 26
#include <lemon/maps.h>
27 27
#include <lemon/core.h>
28 28
#include <lemon/tolerance.h>
29 29

	
30 30
/// \file
31 31
/// \ingroup min_cut
32 32
/// \brief Implementation of the Hao-Orlin algorithm.
33 33
///
34 34
/// Implementation of the Hao-Orlin algorithm class for testing network
35 35
/// reliability.
36 36

	
37 37
namespace lemon {
38 38

	
39 39
  /// \ingroup min_cut
40 40
  ///
41 41
  /// \brief %Hao-Orlin algorithm to find a minimum cut in directed graphs.
42 42
  ///
43 43
  /// Hao-Orlin calculates a minimum cut in a directed graph
44 44
  /// \f$D=(V,A)\f$. It takes a fixed node \f$ source \in V \f$ and
45 45
  /// consists of two phases: in the first phase it determines a
46 46
  /// minimum cut with \f$ source \f$ on the source-side (i.e. a set
47 47
  /// \f$ X\subsetneq V \f$ with \f$ source \in X \f$ and minimal
48 48
  /// out-degree) and in the second phase it determines a minimum cut
49 49
  /// with \f$ source \f$ on the sink-side (i.e. a set
50 50
  /// \f$ X\subsetneq V \f$ with \f$ source \notin X \f$ and minimal
51 51
  /// out-degree). Obviously, the smaller of these two cuts will be a
52 52
  /// minimum cut of \f$ D \f$. The algorithm is a modified
53 53
  /// push-relabel preflow algorithm and our implementation calculates
54 54
  /// the minimum cut in \f$ O(n^2\sqrt{m}) \f$ time (we use the
55 55
  /// highest-label rule), or in \f$O(nm)\f$ for unit capacities. The
56 56
  /// purpose of such algorithm is testing network reliability. For an
57 57
  /// undirected graph you can run just the first phase of the
58 58
  /// algorithm or you can use the algorithm of Nagamochi and Ibaraki
59 59
  /// which solves the undirected problem in
60 60
  /// \f$ O(nm + n^2 \log(n)) \f$ time: it is implemented in the
61 61
  /// NagamochiIbaraki algorithm class.
62 62
  ///
63 63
  /// \param _Digraph is the graph type of the algorithm.
64 64
  /// \param _CapacityMap is an edge map of capacities which should
65 65
  /// be any numreric type. The default type is _Digraph::ArcMap<int>.
66 66
  /// \param _Tolerance is the handler of the inexact computation. The
67 67
  /// default type for this is Tolerance<CapacityMap::Value>.
68 68
#ifdef DOXYGEN
69 69
  template <typename _Digraph, typename _CapacityMap, typename _Tolerance>
70 70
#else
71 71
  template <typename _Digraph,
72 72
            typename _CapacityMap = typename _Digraph::template ArcMap<int>,
73 73
            typename _Tolerance = Tolerance<typename _CapacityMap::Value> >
74 74
#endif
75 75
  class HaoOrlin {
76 76
  private:
77 77

	
78 78
    typedef _Digraph Digraph;
79 79
    typedef _CapacityMap CapacityMap;
80 80
    typedef _Tolerance Tolerance;
81 81

	
82 82
    typedef typename CapacityMap::Value Value;
83 83

	
84 84
    TEMPLATE_GRAPH_TYPEDEFS(Digraph);
85 85

	
86 86
    const Digraph& _graph;
87 87
    const CapacityMap* _capacity;
88 88

	
89 89
    typedef typename Digraph::template ArcMap<Value> FlowMap;
90 90
    FlowMap* _flow;
91 91

	
92 92
    Node _source;
93 93

	
94 94
    int _node_num;
95 95

	
96 96
    // Bucketing structure
97 97
    std::vector<Node> _first, _last;
98 98
    typename Digraph::template NodeMap<Node>* _next;
99 99
    typename Digraph::template NodeMap<Node>* _prev;
100 100
    typename Digraph::template NodeMap<bool>* _active;
101 101
    typename Digraph::template NodeMap<int>* _bucket;
102 102

	
103 103
    std::vector<bool> _dormant;
104 104

	
105 105
    std::list<std::list<int> > _sets;
106 106
    std::list<int>::iterator _highest;
107 107

	
108 108
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
109 109
    ExcessMap* _excess;
110 110

	
111 111
    typedef typename Digraph::template NodeMap<bool> SourceSetMap;
112 112
    SourceSetMap* _source_set;
113 113

	
114 114
    Value _min_cut;
115 115

	
116 116
    typedef typename Digraph::template NodeMap<bool> MinCutMap;
117 117
    MinCutMap* _min_cut_map;
118 118

	
119 119
    Tolerance _tolerance;
120 120

	
121 121
  public:
122 122

	
123 123
    /// \brief Constructor
124 124
    ///
125 125
    /// Constructor of the algorithm class.
126 126
    HaoOrlin(const Digraph& graph, const CapacityMap& capacity,
127 127
             const Tolerance& tolerance = Tolerance()) :
128 128
      _graph(graph), _capacity(&capacity), _flow(0), _source(),
129 129
      _node_num(), _first(), _last(), _next(0), _prev(0),
130 130
      _active(0), _bucket(0), _dormant(), _sets(), _highest(),
131 131
      _excess(0), _source_set(0), _min_cut(), _min_cut_map(0),
132 132
      _tolerance(tolerance) {}
133 133

	
134 134
    ~HaoOrlin() {
135 135
      if (_min_cut_map) {
136 136
        delete _min_cut_map;
137 137
      }
138 138
      if (_source_set) {
139 139
        delete _source_set;
140 140
      }
141 141
      if (_excess) {
142 142
        delete _excess;
143 143
      }
144 144
      if (_next) {
145 145
        delete _next;
146 146
      }
147 147
      if (_prev) {
148 148
        delete _prev;
149 149
      }
150 150
      if (_active) {
151 151
        delete _active;
152 152
      }
153 153
      if (_bucket) {
154 154
        delete _bucket;
155 155
      }
156 156
      if (_flow) {
157 157
        delete _flow;
158 158
      }
159 159
    }
160 160

	
161 161
  private:
162 162

	
163 163
    void activate(const Node& i) {
164 164
      _active->set(i, true);
165 165

	
166 166
      int bucket = (*_bucket)[i];
167 167

	
168 168
      if ((*_prev)[i] == INVALID || (*_active)[(*_prev)[i]]) return;
169 169
      //unlace
170 170
      _next->set((*_prev)[i], (*_next)[i]);
171 171
      if ((*_next)[i] != INVALID) {
172 172
        _prev->set((*_next)[i], (*_prev)[i]);
173 173
      } else {
174 174
        _last[bucket] = (*_prev)[i];
175 175
      }
176 176
      //lace
177 177
      _next->set(i, _first[bucket]);
178 178
      _prev->set(_first[bucket], i);
179 179
      _prev->set(i, INVALID);
180 180
      _first[bucket] = i;
181 181
    }
182 182

	
183 183
    void deactivate(const Node& i) {
184 184
      _active->set(i, false);
185 185
      int bucket = (*_bucket)[i];
186 186

	
187 187
      if ((*_next)[i] == INVALID || !(*_active)[(*_next)[i]]) return;
188 188

	
189 189
      //unlace
190 190
      _prev->set((*_next)[i], (*_prev)[i]);
191 191
      if ((*_prev)[i] != INVALID) {
192 192
        _next->set((*_prev)[i], (*_next)[i]);
193 193
      } else {
194 194
        _first[bucket] = (*_next)[i];
195 195
      }
196 196
      //lace
197 197
      _prev->set(i, _last[bucket]);
198 198
      _next->set(_last[bucket], i);
199 199
      _next->set(i, INVALID);
200 200
      _last[bucket] = i;
201 201
    }
202 202

	
203 203
    void addItem(const Node& i, int bucket) {
204 204
      (*_bucket)[i] = bucket;
205 205
      if (_last[bucket] != INVALID) {
206 206
        _prev->set(i, _last[bucket]);
207 207
        _next->set(_last[bucket], i);
208 208
        _next->set(i, INVALID);
209 209
        _last[bucket] = i;
210 210
      } else {
211 211
        _prev->set(i, INVALID);
212 212
        _first[bucket] = i;
213 213
        _next->set(i, INVALID);
214 214
        _last[bucket] = i;
215 215
      }
216 216
    }
217 217

	
218 218
    void findMinCutOut() {
219 219

	
220 220
      for (NodeIt n(_graph); n != INVALID; ++n) {
221 221
        _excess->set(n, 0);
222 222
      }
223 223

	
224 224
      for (ArcIt a(_graph); a != INVALID; ++a) {
225 225
        _flow->set(a, 0);
226 226
      }
227 227

	
228
      int bucket_num = 1;
228
      int bucket_num = 0;
229
      std::vector<Node> queue(_node_num);
230
      int qfirst = 0, qlast = 0, qsep = 0;
229 231

	
230 232
      {
231 233
        typename Digraph::template NodeMap<bool> reached(_graph, false);
232 234

	
233 235
        reached.set(_source, true);
234

	
235 236
        bool first_set = true;
236 237

	
237 238
        for (NodeIt t(_graph); t != INVALID; ++t) {
238 239
          if (reached[t]) continue;
239 240
          _sets.push_front(std::list<int>());
240
          _sets.front().push_front(bucket_num);
241
          _dormant[bucket_num] = !first_set;
242

	
243
          _bucket->set(t, bucket_num);
244
          _first[bucket_num] = _last[bucket_num] = t;
245
          _next->set(t, INVALID);
246
          _prev->set(t, INVALID);
247

	
248
          ++bucket_num;
249

	
250
          std::vector<Node> queue;
251
          queue.push_back(t);
241
          
242
          queue[qlast++] = t;
252 243
          reached.set(t, true);
253 244

	
254
          while (!queue.empty()) {
255
            _sets.front().push_front(bucket_num);
256
            _dormant[bucket_num] = !first_set;
257
            _first[bucket_num] = _last[bucket_num] = INVALID;
245
          while (qfirst != qlast) {
246
            if (qsep == qfirst) {
247
              ++bucket_num;
248
              _sets.front().push_front(bucket_num);
249
              _dormant[bucket_num] = !first_set;
250
              _first[bucket_num] = _last[bucket_num] = INVALID;
251
              qsep = qlast;
252
            }
258 253

	
259
            std::vector<Node> nqueue;
260
            for (int i = 0; i < int(queue.size()); ++i) {
261
              Node n = queue[i];
262
              for (InArcIt a(_graph, n); a != INVALID; ++a) {
263
                Node u = _graph.source(a);
264
                if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
265
                  reached.set(u, true);
266
                  addItem(u, bucket_num);
267
                  nqueue.push_back(u);
268
                }
254
            Node n = queue[qfirst++];
255
            addItem(n, bucket_num);
256

	
257
            for (InArcIt a(_graph, n); a != INVALID; ++a) {
258
              Node u = _graph.source(a);
259
              if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
260
                reached.set(u, true);
261
                queue[qlast++] = u;
269 262
              }
270 263
            }
271
            queue.swap(nqueue);
272
            ++bucket_num;
273 264
          }
274
          _sets.front().pop_front();
275
          --bucket_num;
276 265
          first_set = false;
277 266
        }
278 267

	
268
        ++bucket_num;
279 269
        _bucket->set(_source, 0);
280 270
        _dormant[0] = true;
281 271
      }
282 272
      _source_set->set(_source, true);
283 273

	
284 274
      Node target = _last[_sets.back().back()];
285 275
      {
286 276
        for (OutArcIt a(_graph, _source); a != INVALID; ++a) {
287 277
          if (_tolerance.positive((*_capacity)[a])) {
288 278
            Node u = _graph.target(a);
289 279
            _flow->set(a, (*_capacity)[a]);
290 280
            _excess->set(u, (*_excess)[u] + (*_capacity)[a]);
291 281
            if (!(*_active)[u] && u != _source) {
292 282
              activate(u);
293 283
            }
294 284
          }
295 285
        }
296 286

	
297 287
        if ((*_active)[target]) {
298 288
          deactivate(target);
299 289
        }
300 290

	
301 291
        _highest = _sets.back().begin();
302 292
        while (_highest != _sets.back().end() &&
303 293
               !(*_active)[_first[*_highest]]) {
304 294
          ++_highest;
305 295
        }
306 296
      }
307 297

	
308 298
      while (true) {
309 299
        while (_highest != _sets.back().end()) {
310 300
          Node n = _first[*_highest];
311 301
          Value excess = (*_excess)[n];
312 302
          int next_bucket = _node_num;
313 303

	
314 304
          int under_bucket;
315 305
          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
316 306
            under_bucket = -1;
317 307
          } else {
318 308
            under_bucket = *(++std::list<int>::iterator(_highest));
319 309
          }
320 310

	
321 311
          for (OutArcIt a(_graph, n); a != INVALID; ++a) {
322 312
            Node v = _graph.target(a);
323 313
            if (_dormant[(*_bucket)[v]]) continue;
324 314
            Value rem = (*_capacity)[a] - (*_flow)[a];
325 315
            if (!_tolerance.positive(rem)) continue;
326 316
            if ((*_bucket)[v] == under_bucket) {
327 317
              if (!(*_active)[v] && v != target) {
328 318
                activate(v);
329 319
              }
330 320
              if (!_tolerance.less(rem, excess)) {
331 321
                _flow->set(a, (*_flow)[a] + excess);
332 322
                _excess->set(v, (*_excess)[v] + excess);
333 323
                excess = 0;
334 324
                goto no_more_push;
335 325
              } else {
336 326
                excess -= rem;
337 327
                _excess->set(v, (*_excess)[v] + rem);
338 328
                _flow->set(a, (*_capacity)[a]);
339 329
              }
340 330
            } else if (next_bucket > (*_bucket)[v]) {
341 331
              next_bucket = (*_bucket)[v];
342 332
            }
343 333
          }
344 334

	
345 335
          for (InArcIt a(_graph, n); a != INVALID; ++a) {
346 336
            Node v = _graph.source(a);
347 337
            if (_dormant[(*_bucket)[v]]) continue;
348 338
            Value rem = (*_flow)[a];
349 339
            if (!_tolerance.positive(rem)) continue;
350 340
            if ((*_bucket)[v] == under_bucket) {
351 341
              if (!(*_active)[v] && v != target) {
352 342
                activate(v);
353 343
              }
354 344
              if (!_tolerance.less(rem, excess)) {
355 345
                _flow->set(a, (*_flow)[a] - excess);
356 346
                _excess->set(v, (*_excess)[v] + excess);
357 347
                excess = 0;
358 348
                goto no_more_push;
359 349
              } else {
360 350
                excess -= rem;
361 351
                _excess->set(v, (*_excess)[v] + rem);
362 352
                _flow->set(a, 0);
363 353
              }
364 354
            } else if (next_bucket > (*_bucket)[v]) {
365 355
              next_bucket = (*_bucket)[v];
366 356
            }
367 357
          }
368 358

	
369 359
        no_more_push:
370 360

	
371 361
          _excess->set(n, excess);
372 362

	
373 363
          if (excess != 0) {
374 364
            if ((*_next)[n] == INVALID) {
375 365
              typename std::list<std::list<int> >::iterator new_set =
376 366
                _sets.insert(--_sets.end(), std::list<int>());
377 367
              new_set->splice(new_set->end(), _sets.back(),
378 368
                              _sets.back().begin(), ++_highest);
379 369
              for (std::list<int>::iterator it = new_set->begin();
380 370
                   it != new_set->end(); ++it) {
381 371
                _dormant[*it] = true;
382 372
              }
383 373
              while (_highest != _sets.back().end() &&
384 374
                     !(*_active)[_first[*_highest]]) {
385 375
                ++_highest;
386 376
              }
387 377
            } else if (next_bucket == _node_num) {
388 378
              _first[(*_bucket)[n]] = (*_next)[n];
389 379
              _prev->set((*_next)[n], INVALID);
390 380

	
391 381
              std::list<std::list<int> >::iterator new_set =
392 382
                _sets.insert(--_sets.end(), std::list<int>());
393 383

	
394 384
              new_set->push_front(bucket_num);
395 385
              _bucket->set(n, bucket_num);
396 386
              _first[bucket_num] = _last[bucket_num] = n;
397 387
              _next->set(n, INVALID);
398 388
              _prev->set(n, INVALID);
399 389
              _dormant[bucket_num] = true;
400 390
              ++bucket_num;
401 391

	
402 392
              while (_highest != _sets.back().end() &&
403 393
                     !(*_active)[_first[*_highest]]) {
404 394
                ++_highest;
405 395
              }
406 396
            } else {
407 397
              _first[*_highest] = (*_next)[n];
408 398
              _prev->set((*_next)[n], INVALID);
409 399

	
410 400
              while (next_bucket != *_highest) {
411 401
                --_highest;
412 402
              }
413 403

	
414 404
              if (_highest == _sets.back().begin()) {
415 405
                _sets.back().push_front(bucket_num);
416 406
                _dormant[bucket_num] = false;
417 407
                _first[bucket_num] = _last[bucket_num] = INVALID;
418 408
                ++bucket_num;
419 409
              }
420 410
              --_highest;
421 411

	
422 412
              _bucket->set(n, *_highest);
423 413
              _next->set(n, _first[*_highest]);
424 414
              if (_first[*_highest] != INVALID) {
425 415
                _prev->set(_first[*_highest], n);
426 416
              } else {
427 417
                _last[*_highest] = n;
428 418
              }
429 419
              _first[*_highest] = n;
430 420
            }
431 421
          } else {
432 422

	
433 423
            deactivate(n);
434 424
            if (!(*_active)[_first[*_highest]]) {
435 425
              ++_highest;
436 426
              if (_highest != _sets.back().end() &&
437 427
                  !(*_active)[_first[*_highest]]) {
438 428
                _highest = _sets.back().end();
439 429
              }
440 430
            }
441 431
          }
442 432
        }
443 433

	
444 434
        if ((*_excess)[target] < _min_cut) {
445 435
          _min_cut = (*_excess)[target];
446 436
          for (NodeIt i(_graph); i != INVALID; ++i) {
447 437
            _min_cut_map->set(i, true);
448 438
          }
449 439
          for (std::list<int>::iterator it = _sets.back().begin();
450 440
               it != _sets.back().end(); ++it) {
451 441
            Node n = _first[*it];
452 442
            while (n != INVALID) {
453 443
              _min_cut_map->set(n, false);
454 444
              n = (*_next)[n];
455 445
            }
456 446
          }
457 447
        }
458 448

	
459 449
        {
460 450
          Node new_target;
461 451
          if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
462 452
            if ((*_next)[target] == INVALID) {
463 453
              _last[(*_bucket)[target]] = (*_prev)[target];
464 454
              new_target = (*_prev)[target];
465 455
            } else {
466 456
              _prev->set((*_next)[target], (*_prev)[target]);
467 457
              new_target = (*_next)[target];
468 458
            }
469 459
            if ((*_prev)[target] == INVALID) {
470 460
              _first[(*_bucket)[target]] = (*_next)[target];
471 461
            } else {
472 462
              _next->set((*_prev)[target], (*_next)[target]);
473 463
            }
474 464
          } else {
475 465
            _sets.back().pop_back();
476 466
            if (_sets.back().empty()) {
477 467
              _sets.pop_back();
478 468
              if (_sets.empty())
479 469
                break;
480 470
              for (std::list<int>::iterator it = _sets.back().begin();
481 471
                   it != _sets.back().end(); ++it) {
482 472
                _dormant[*it] = false;
483 473
              }
484 474
            }
485 475
            new_target = _last[_sets.back().back()];
486 476
          }
487 477

	
488 478
          _bucket->set(target, 0);
489 479

	
490 480
          _source_set->set(target, true);
491 481
          for (OutArcIt a(_graph, target); a != INVALID; ++a) {
492 482
            Value rem = (*_capacity)[a] - (*_flow)[a];
493 483
            if (!_tolerance.positive(rem)) continue;
494 484
            Node v = _graph.target(a);
495 485
            if (!(*_active)[v] && !(*_source_set)[v]) {
496 486
              activate(v);
497 487
            }
498 488
            _excess->set(v, (*_excess)[v] + rem);
499 489
            _flow->set(a, (*_capacity)[a]);
500 490
          }
501 491

	
502 492
          for (InArcIt a(_graph, target); a != INVALID; ++a) {
503 493
            Value rem = (*_flow)[a];
504 494
            if (!_tolerance.positive(rem)) continue;
505 495
            Node v = _graph.source(a);
506 496
            if (!(*_active)[v] && !(*_source_set)[v]) {
507 497
              activate(v);
508 498
            }
509 499
            _excess->set(v, (*_excess)[v] + rem);
510 500
            _flow->set(a, 0);
511 501
          }
512 502

	
513 503
          target = new_target;
514 504
          if ((*_active)[target]) {
515 505
            deactivate(target);
516 506
          }
517 507

	
518 508
          _highest = _sets.back().begin();
519 509
          while (_highest != _sets.back().end() &&
520 510
                 !(*_active)[_first[*_highest]]) {
521 511
            ++_highest;
522 512
          }
523 513
        }
524 514
      }
525 515
    }
526 516

	
527 517
    void findMinCutIn() {
528 518

	
529 519
      for (NodeIt n(_graph); n != INVALID; ++n) {
530 520
        _excess->set(n, 0);
531 521
      }
532 522

	
533 523
      for (ArcIt a(_graph); a != INVALID; ++a) {
534 524
        _flow->set(a, 0);
535 525
      }
536 526

	
537
      int bucket_num = 1;
527
      int bucket_num = 0;
528
      std::vector<Node> queue(_node_num);
529
      int qfirst = 0, qlast = 0, qsep = 0;
538 530

	
539 531
      {
540 532
        typename Digraph::template NodeMap<bool> reached(_graph, false);
541 533

	
542 534
        reached.set(_source, true);
543 535

	
544 536
        bool first_set = true;
545 537

	
546 538
        for (NodeIt t(_graph); t != INVALID; ++t) {
547 539
          if (reached[t]) continue;
548 540
          _sets.push_front(std::list<int>());
549
          _sets.front().push_front(bucket_num);
550
          _dormant[bucket_num] = !first_set;
551

	
552
          _bucket->set(t, bucket_num);
553
          _first[bucket_num] = _last[bucket_num] = t;
554
          _next->set(t, INVALID);
555
          _prev->set(t, INVALID);
556

	
557
          ++bucket_num;
558

	
559
          std::vector<Node> queue;
560
          queue.push_back(t);
541
          
542
          queue[qlast++] = t;
561 543
          reached.set(t, true);
562 544

	
563
          while (!queue.empty()) {
564
            _sets.front().push_front(bucket_num);
565
            _dormant[bucket_num] = !first_set;
566
            _first[bucket_num] = _last[bucket_num] = INVALID;
545
          while (qfirst != qlast) {
546
            if (qsep == qfirst) {
547
              ++bucket_num;
548
              _sets.front().push_front(bucket_num);
549
              _dormant[bucket_num] = !first_set;
550
              _first[bucket_num] = _last[bucket_num] = INVALID;
551
              qsep = qlast;
552
            }
567 553

	
568
            std::vector<Node> nqueue;
569
            for (int i = 0; i < int(queue.size()); ++i) {
570
              Node n = queue[i];
571
              for (OutArcIt a(_graph, n); a != INVALID; ++a) {
572
                Node u = _graph.target(a);
573
                if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
574
                  reached.set(u, true);
575
                  addItem(u, bucket_num);
576
                  nqueue.push_back(u);
577
                }
554
            Node n = queue[qfirst++];
555
            addItem(n, bucket_num);
556

	
557
            for (OutArcIt a(_graph, n); a != INVALID; ++a) {
558
              Node u = _graph.target(a);
559
              if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
560
                reached.set(u, true);
561
                queue[qlast++] = u;
578 562
              }
579 563
            }
580
            queue.swap(nqueue);
581
            ++bucket_num;
582 564
          }
583
          _sets.front().pop_front();
584
          --bucket_num;
585 565
          first_set = false;
586 566
        }
587 567

	
568
        ++bucket_num;
588 569
        _bucket->set(_source, 0);
589 570
        _dormant[0] = true;
590 571
      }
591 572
      _source_set->set(_source, true);
592 573

	
593 574
      Node target = _last[_sets.back().back()];
594 575
      {
595 576
        for (InArcIt a(_graph, _source); a != INVALID; ++a) {
596 577
          if (_tolerance.positive((*_capacity)[a])) {
597 578
            Node u = _graph.source(a);
598 579
            _flow->set(a, (*_capacity)[a]);
599 580
            _excess->set(u, (*_excess)[u] + (*_capacity)[a]);
600 581
            if (!(*_active)[u] && u != _source) {
601 582
              activate(u);
602 583
            }
603 584
          }
604 585
        }
605 586
        if ((*_active)[target]) {
606 587
          deactivate(target);
607 588
        }
608 589

	
609 590
        _highest = _sets.back().begin();
610 591
        while (_highest != _sets.back().end() &&
611 592
               !(*_active)[_first[*_highest]]) {
612 593
          ++_highest;
613 594
        }
614 595
      }
615 596

	
616 597

	
617 598
      while (true) {
618 599
        while (_highest != _sets.back().end()) {
619 600
          Node n = _first[*_highest];
620 601
          Value excess = (*_excess)[n];
621 602
          int next_bucket = _node_num;
622 603

	
623 604
          int under_bucket;
624 605
          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
625 606
            under_bucket = -1;
626 607
          } else {
627 608
            under_bucket = *(++std::list<int>::iterator(_highest));
628 609
          }
629 610

	
630 611
          for (InArcIt a(_graph, n); a != INVALID; ++a) {
631 612
            Node v = _graph.source(a);
632 613
            if (_dormant[(*_bucket)[v]]) continue;
633 614
            Value rem = (*_capacity)[a] - (*_flow)[a];
634 615
            if (!_tolerance.positive(rem)) continue;
635 616
            if ((*_bucket)[v] == under_bucket) {
636 617
              if (!(*_active)[v] && v != target) {
637 618
                activate(v);
638 619
              }
639 620
              if (!_tolerance.less(rem, excess)) {
640 621
                _flow->set(a, (*_flow)[a] + excess);
641 622
                _excess->set(v, (*_excess)[v] + excess);
642 623
                excess = 0;
643 624
                goto no_more_push;
644 625
              } else {
645 626
                excess -= rem;
646 627
                _excess->set(v, (*_excess)[v] + rem);
647 628
                _flow->set(a, (*_capacity)[a]);
648 629
              }
649 630
            } else if (next_bucket > (*_bucket)[v]) {
650 631
              next_bucket = (*_bucket)[v];
651 632
            }
652 633
          }
653 634

	
654 635
          for (OutArcIt a(_graph, n); a != INVALID; ++a) {
655 636
            Node v = _graph.target(a);
656 637
            if (_dormant[(*_bucket)[v]]) continue;
657 638
            Value rem = (*_flow)[a];
658 639
            if (!_tolerance.positive(rem)) continue;
659 640
            if ((*_bucket)[v] == under_bucket) {
660 641
              if (!(*_active)[v] && v != target) {
661 642
                activate(v);
662 643
              }
663 644
              if (!_tolerance.less(rem, excess)) {
664 645
                _flow->set(a, (*_flow)[a] - excess);
665 646
                _excess->set(v, (*_excess)[v] + excess);
666 647
                excess = 0;
667 648
                goto no_more_push;
668 649
              } else {
669 650
                excess -= rem;
670 651
                _excess->set(v, (*_excess)[v] + rem);
671 652
                _flow->set(a, 0);
672 653
              }
673 654
            } else if (next_bucket > (*_bucket)[v]) {
674 655
              next_bucket = (*_bucket)[v];
675 656
            }
676 657
          }
677 658

	
678 659
        no_more_push:
679 660

	
680 661
          _excess->set(n, excess);
681 662

	
682 663
          if (excess != 0) {
683 664
            if ((*_next)[n] == INVALID) {
684 665
              typename std::list<std::list<int> >::iterator new_set =
685 666
                _sets.insert(--_sets.end(), std::list<int>());
686 667
              new_set->splice(new_set->end(), _sets.back(),
687 668
                              _sets.back().begin(), ++_highest);
688 669
              for (std::list<int>::iterator it = new_set->begin();
689 670
                   it != new_set->end(); ++it) {
690 671
                _dormant[*it] = true;
691 672
              }
692 673
              while (_highest != _sets.back().end() &&
693 674
                     !(*_active)[_first[*_highest]]) {
694 675
                ++_highest;
695 676
              }
696 677
            } else if (next_bucket == _node_num) {
697 678
              _first[(*_bucket)[n]] = (*_next)[n];
698 679
              _prev->set((*_next)[n], INVALID);
699 680

	
700 681
              std::list<std::list<int> >::iterator new_set =
701 682
                _sets.insert(--_sets.end(), std::list<int>());
702 683

	
703 684
              new_set->push_front(bucket_num);
704 685
              _bucket->set(n, bucket_num);
705 686
              _first[bucket_num] = _last[bucket_num] = n;
706 687
              _next->set(n, INVALID);
707 688
              _prev->set(n, INVALID);
708 689
              _dormant[bucket_num] = true;
709 690
              ++bucket_num;
710 691

	
711 692
              while (_highest != _sets.back().end() &&
712 693
                     !(*_active)[_first[*_highest]]) {
713 694
                ++_highest;
714 695
              }
715 696
            } else {
716 697
              _first[*_highest] = (*_next)[n];
717 698
              _prev->set((*_next)[n], INVALID);
718 699

	
719 700
              while (next_bucket != *_highest) {
720 701
                --_highest;
721 702
              }
722 703
              if (_highest == _sets.back().begin()) {
723 704
                _sets.back().push_front(bucket_num);
724 705
                _dormant[bucket_num] = false;
725 706
                _first[bucket_num] = _last[bucket_num] = INVALID;
726 707
                ++bucket_num;
727 708
              }
728 709
              --_highest;
729 710

	
730 711
              _bucket->set(n, *_highest);
731 712
              _next->set(n, _first[*_highest]);
732 713
              if (_first[*_highest] != INVALID) {
733 714
                _prev->set(_first[*_highest], n);
734 715
              } else {
735 716
                _last[*_highest] = n;
736 717
              }
737 718
              _first[*_highest] = n;
738 719
            }
739 720
          } else {
740 721

	
741 722
            deactivate(n);
742 723
            if (!(*_active)[_first[*_highest]]) {
743 724
              ++_highest;
744 725
              if (_highest != _sets.back().end() &&
745 726
                  !(*_active)[_first[*_highest]]) {
746 727
                _highest = _sets.back().end();
747 728
              }
748 729
            }
749 730
          }
750 731
        }
751 732

	
752 733
        if ((*_excess)[target] < _min_cut) {
753 734
          _min_cut = (*_excess)[target];
754 735
          for (NodeIt i(_graph); i != INVALID; ++i) {
755 736
            _min_cut_map->set(i, false);
756 737
          }
757 738
          for (std::list<int>::iterator it = _sets.back().begin();
758 739
               it != _sets.back().end(); ++it) {
759 740
            Node n = _first[*it];
760 741
            while (n != INVALID) {
761 742
              _min_cut_map->set(n, true);
762 743
              n = (*_next)[n];
763 744
            }
764 745
          }
765 746
        }
766 747

	
767 748
        {
768 749
          Node new_target;
769 750
          if ((*_prev)[target] != INVALID || (*_next)[target] != INVALID) {
770 751
            if ((*_next)[target] == INVALID) {
771 752
              _last[(*_bucket)[target]] = (*_prev)[target];
772 753
              new_target = (*_prev)[target];
773 754
            } else {
774 755
              _prev->set((*_next)[target], (*_prev)[target]);
775 756
              new_target = (*_next)[target];
776 757
            }
777 758
            if ((*_prev)[target] == INVALID) {
778 759
              _first[(*_bucket)[target]] = (*_next)[target];
779 760
            } else {
780 761
              _next->set((*_prev)[target], (*_next)[target]);
781 762
            }
782 763
          } else {
783 764
            _sets.back().pop_back();
784 765
            if (_sets.back().empty()) {
785 766
              _sets.pop_back();
786 767
              if (_sets.empty())
787 768
                break;
788 769
              for (std::list<int>::iterator it = _sets.back().begin();
789 770
                   it != _sets.back().end(); ++it) {
790 771
                _dormant[*it] = false;
791 772
              }
792 773
            }
793 774
            new_target = _last[_sets.back().back()];
794 775
          }
795 776

	
796 777
          _bucket->set(target, 0);
797 778

	
798 779
          _source_set->set(target, true);
799 780
          for (InArcIt a(_graph, target); a != INVALID; ++a) {
800 781
            Value rem = (*_capacity)[a] - (*_flow)[a];
801 782
            if (!_tolerance.positive(rem)) continue;
802 783
            Node v = _graph.source(a);
803 784
            if (!(*_active)[v] && !(*_source_set)[v]) {
804 785
              activate(v);
805 786
            }
806 787
            _excess->set(v, (*_excess)[v] + rem);
807 788
            _flow->set(a, (*_capacity)[a]);
808 789
          }
809 790

	
810 791
          for (OutArcIt a(_graph, target); a != INVALID; ++a) {
811 792
            Value rem = (*_flow)[a];
812 793
            if (!_tolerance.positive(rem)) continue;
813 794
            Node v = _graph.target(a);
814 795
            if (!(*_active)[v] && !(*_source_set)[v]) {
815 796
              activate(v);
816 797
            }
817 798
            _excess->set(v, (*_excess)[v] + rem);
818 799
            _flow->set(a, 0);
819 800
          }
820 801

	
821 802
          target = new_target;
822 803
          if ((*_active)[target]) {
823 804
            deactivate(target);
824 805
          }
825 806

	
826 807
          _highest = _sets.back().begin();
827 808
          while (_highest != _sets.back().end() &&
828 809
                 !(*_active)[_first[*_highest]]) {
829 810
            ++_highest;
830 811
          }
831 812
        }
832 813
      }
833 814
    }
834 815

	
835 816
  public:
836 817

	
837 818
    /// \name Execution control
838 819
    /// The simplest way to execute the algorithm is to use
839 820
    /// one of the member functions called \c run(...).
840 821
    /// \n
841 822
    /// If you need more control on the execution,
842 823
    /// first you must call \ref init(), then the \ref calculateIn() or
843 824
    /// \ref calculateIn() functions.
844 825

	
845 826
    /// @{
846 827

	
847 828
    /// \brief Initializes the internal data structures.
848 829
    ///
849 830
    /// Initializes the internal data structures. It creates
850 831
    /// the maps, residual graph adaptors and some bucket structures
851 832
    /// for the algorithm.
852 833
    void init() {
853 834
      init(NodeIt(_graph));
854 835
    }
855 836

	
856 837
    /// \brief Initializes the internal data structures.
857 838
    ///
858 839
    /// Initializes the internal data structures. It creates
859 840
    /// the maps, residual graph adaptor and some bucket structures
860 841
    /// for the algorithm. Node \c source  is used as the push-relabel
861 842
    /// algorithm's source.
862 843
    void init(const Node& source) {
863 844
      _source = source;
864 845

	
865 846
      _node_num = countNodes(_graph);
866 847

	
867
      _first.resize(_node_num + 1);
868
      _last.resize(_node_num + 1);
848
      _first.resize(_node_num);
849
      _last.resize(_node_num);
869 850

	
870
      _dormant.resize(_node_num + 1);
851
      _dormant.resize(_node_num);
871 852

	
872 853
      if (!_flow) {
873 854
        _flow = new FlowMap(_graph);
874 855
      }
875 856
      if (!_next) {
876 857
        _next = new typename Digraph::template NodeMap<Node>(_graph);
877 858
      }
878 859
      if (!_prev) {
879 860
        _prev = new typename Digraph::template NodeMap<Node>(_graph);
880 861
      }
881 862
      if (!_active) {
882 863
        _active = new typename Digraph::template NodeMap<bool>(_graph);
883 864
      }
884 865
      if (!_bucket) {
885 866
        _bucket = new typename Digraph::template NodeMap<int>(_graph);
886 867
      }
887 868
      if (!_excess) {
888 869
        _excess = new ExcessMap(_graph);
889 870
      }
890 871
      if (!_source_set) {
891 872
        _source_set = new SourceSetMap(_graph);
892 873
      }
893 874
      if (!_min_cut_map) {
894 875
        _min_cut_map = new MinCutMap(_graph);
895 876
      }
896 877

	
897 878
      _min_cut = std::numeric_limits<Value>::max();
898 879
    }
899 880

	
900 881

	
901 882
    /// \brief Calculates a minimum cut with \f$ source \f$ on the
902 883
    /// source-side.
903 884
    ///
904 885
    /// Calculates a minimum cut with \f$ source \f$ on the
905 886
    /// source-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source
906 887
    /// \in X \f$ and minimal out-degree).
907 888
    void calculateOut() {
908 889
      findMinCutOut();
909 890
    }
910 891

	
911 892
    /// \brief Calculates a minimum cut with \f$ source \f$ on the
912 893
    /// target-side.
913 894
    ///
914 895
    /// Calculates a minimum cut with \f$ source \f$ on the
915 896
    /// target-side (i.e. a set \f$ X\subsetneq V \f$ with \f$ source
916 897
    /// \in X \f$ and minimal out-degree).
917 898
    void calculateIn() {
918 899
      findMinCutIn();
919 900
    }
920 901

	
921 902

	
922 903
    /// \brief Runs the algorithm.
923 904
    ///
924 905
    /// Runs the algorithm. It finds nodes \c source and \c target
925 906
    /// arbitrarily and then calls \ref init(), \ref calculateOut()
926 907
    /// and \ref calculateIn().
927 908
    void run() {
928 909
      init();
929 910
      calculateOut();
930 911
      calculateIn();
931 912
    }
932 913

	
933 914
    /// \brief Runs the algorithm.
934 915
    ///
935 916
    /// Runs the algorithm. It uses the given \c source node, finds a
936 917
    /// proper \c target and then calls the \ref init(), \ref
937 918
    /// calculateOut() and \ref calculateIn().
938 919
    void run(const Node& s) {
939 920
      init(s);
940 921
      calculateOut();
941 922
      calculateIn();
942 923
    }
943 924

	
944 925
    /// @}
945 926

	
946 927
    /// \name Query Functions
947 928
    /// The result of the %HaoOrlin algorithm
948 929
    /// can be obtained using these functions.
949 930
    /// \n
950 931
    /// Before using these functions, either \ref run(), \ref
951 932
    /// calculateOut() or \ref calculateIn() must be called.
952 933

	
953 934
    /// @{
954 935

	
955 936
    /// \brief Returns the value of the minimum value cut.
956 937
    ///
957 938
    /// Returns the value of the minimum value cut.
958 939
    Value minCutValue() const {
959 940
      return _min_cut;
960 941
    }
961 942

	
962 943

	
963 944
    /// \brief Returns a minimum cut.
964 945
    ///
965 946
    /// Sets \c nodeMap to the characteristic vector of a minimum
966 947
    /// value cut: it will give a nonempty set \f$ X\subsetneq V \f$
967 948
    /// with minimal out-degree (i.e. \c nodeMap will be true exactly
968 949
    /// for the nodes of \f$ X \f$).  \pre nodeMap should be a
969 950
    /// bool-valued node-map.
970 951
    template <typename NodeMap>
971 952
    Value minCutMap(NodeMap& nodeMap) const {
972 953
      for (NodeIt it(_graph); it != INVALID; ++it) {
973 954
        nodeMap.set(it, (*_min_cut_map)[it]);
974 955
      }
975 956
      return _min_cut;
976 957
    }
977 958

	
978 959
    /// @}
979 960

	
980 961
  }; //class HaoOrlin
981 962

	
982 963

	
983 964
} //namespace lemon
984 965

	
985 966
#endif //LEMON_HAO_ORLIN_H
0 comments (0 inline)