gravatar
deba@inf.elte.hu
deba@inf.elte.hu
Fix in HaoOrlin (#264)
0 2 0
default
2 files changed with 14 insertions and 24 deletions:
↑ Collapse diff ↑
Ignore white space 6 line context
... ...
@@ -37,683 +37,685 @@
37 37
namespace lemon {
38 38

	
39 39
  /// \ingroup min_cut
40 40
  ///
41 41
  /// \brief Hao-Orlin algorithm for finding a minimum cut in a digraph.
42 42
  ///
43 43
  /// This class implements the Hao-Orlin algorithm for finding a minimum
44 44
  /// value cut in a directed graph \f$D=(V,A)\f$. 
45 45
  /// It takes a fixed node \f$ source \in V \f$ and
46 46
  /// consists of two phases: in the first phase it determines a
47 47
  /// minimum cut with \f$ source \f$ on the source-side (i.e. a set
48 48
  /// \f$ X\subsetneq V \f$ with \f$ source \in X \f$ and minimal outgoing
49 49
  /// capacity) and in the second phase it determines a minimum cut
50 50
  /// with \f$ source \f$ on the sink-side (i.e. a set
51 51
  /// \f$ X\subsetneq V \f$ with \f$ source \notin X \f$ and minimal outgoing
52 52
  /// capacity). Obviously, the smaller of these two cuts will be a
53 53
  /// minimum cut of \f$ D \f$. The algorithm is a modified
54 54
  /// preflow push-relabel algorithm. Our implementation calculates
55 55
  /// the minimum cut in \f$ O(n^2\sqrt{m}) \f$ time (we use the
56 56
  /// highest-label rule), or in \f$O(nm)\f$ for unit capacities. The
57 57
  /// purpose of such algorithm is e.g. testing network reliability.
58 58
  ///
59 59
  /// For an undirected graph you can run just the first phase of the
60 60
  /// algorithm or you can use the algorithm of Nagamochi and Ibaraki,
61 61
  /// which solves the undirected problem in \f$ O(nm + n^2 \log n) \f$ 
62 62
  /// time. It is implemented in the NagamochiIbaraki algorithm class.
63 63
  ///
64 64
  /// \tparam GR The type of the digraph the algorithm runs on.
65 65
  /// \tparam CAP The type of the arc map containing the capacities,
66 66
  /// which can be any numreric type. The default map type is
67 67
  /// \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
68 68
  /// \tparam TOL Tolerance class for handling inexact computations. The
69 69
  /// default tolerance type is \ref Tolerance "Tolerance<CAP::Value>".
70 70
#ifdef DOXYGEN
71 71
  template <typename GR, typename CAP, typename TOL>
72 72
#else
73 73
  template <typename GR,
74 74
            typename CAP = typename GR::template ArcMap<int>,
75 75
            typename TOL = Tolerance<typename CAP::Value> >
76 76
#endif
77 77
  class HaoOrlin {
78 78
  public:
79 79
   
80 80
    /// The digraph type of the algorithm
81 81
    typedef GR Digraph;
82 82
    /// The capacity map type of the algorithm
83 83
    typedef CAP CapacityMap;
84 84
    /// The tolerance type of the algorithm
85 85
    typedef TOL Tolerance;
86 86

	
87 87
  private:
88 88

	
89 89
    typedef typename CapacityMap::Value Value;
90 90

	
91 91
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
92 92

	
93 93
    const Digraph& _graph;
94 94
    const CapacityMap* _capacity;
95 95

	
96 96
    typedef typename Digraph::template ArcMap<Value> FlowMap;
97 97
    FlowMap* _flow;
98 98

	
99 99
    Node _source;
100 100

	
101 101
    int _node_num;
102 102

	
103 103
    // Bucketing structure
104 104
    std::vector<Node> _first, _last;
105 105
    typename Digraph::template NodeMap<Node>* _next;
106 106
    typename Digraph::template NodeMap<Node>* _prev;
107 107
    typename Digraph::template NodeMap<bool>* _active;
108 108
    typename Digraph::template NodeMap<int>* _bucket;
109 109

	
110 110
    std::vector<bool> _dormant;
111 111

	
112 112
    std::list<std::list<int> > _sets;
113 113
    std::list<int>::iterator _highest;
114 114

	
115 115
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
116 116
    ExcessMap* _excess;
117 117

	
118 118
    typedef typename Digraph::template NodeMap<bool> SourceSetMap;
119 119
    SourceSetMap* _source_set;
120 120

	
121 121
    Value _min_cut;
122 122

	
123 123
    typedef typename Digraph::template NodeMap<bool> MinCutMap;
124 124
    MinCutMap* _min_cut_map;
125 125

	
126 126
    Tolerance _tolerance;
127 127

	
128 128
  public:
129 129

	
130 130
    /// \brief Constructor
131 131
    ///
132 132
    /// Constructor of the algorithm class.
133 133
    HaoOrlin(const Digraph& graph, const CapacityMap& capacity,
134 134
             const Tolerance& tolerance = Tolerance()) :
135 135
      _graph(graph), _capacity(&capacity), _flow(0), _source(),
136 136
      _node_num(), _first(), _last(), _next(0), _prev(0),
137 137
      _active(0), _bucket(0), _dormant(), _sets(), _highest(),
138 138
      _excess(0), _source_set(0), _min_cut(), _min_cut_map(0),
139 139
      _tolerance(tolerance) {}
140 140

	
141 141
    ~HaoOrlin() {
142 142
      if (_min_cut_map) {
143 143
        delete _min_cut_map;
144 144
      }
145 145
      if (_source_set) {
146 146
        delete _source_set;
147 147
      }
148 148
      if (_excess) {
149 149
        delete _excess;
150 150
      }
151 151
      if (_next) {
152 152
        delete _next;
153 153
      }
154 154
      if (_prev) {
155 155
        delete _prev;
156 156
      }
157 157
      if (_active) {
158 158
        delete _active;
159 159
      }
160 160
      if (_bucket) {
161 161
        delete _bucket;
162 162
      }
163 163
      if (_flow) {
164 164
        delete _flow;
165 165
      }
166 166
    }
167 167

	
168 168
  private:
169 169

	
170 170
    void activate(const Node& i) {
171 171
      (*_active)[i] = true;
172 172

	
173 173
      int bucket = (*_bucket)[i];
174 174

	
175 175
      if ((*_prev)[i] == INVALID || (*_active)[(*_prev)[i]]) return;
176 176
      //unlace
177 177
      (*_next)[(*_prev)[i]] = (*_next)[i];
178 178
      if ((*_next)[i] != INVALID) {
179 179
        (*_prev)[(*_next)[i]] = (*_prev)[i];
180 180
      } else {
181 181
        _last[bucket] = (*_prev)[i];
182 182
      }
183 183
      //lace
184 184
      (*_next)[i] = _first[bucket];
185 185
      (*_prev)[_first[bucket]] = i;
186 186
      (*_prev)[i] = INVALID;
187 187
      _first[bucket] = i;
188 188
    }
189 189

	
190 190
    void deactivate(const Node& i) {
191 191
      (*_active)[i] = false;
192 192
      int bucket = (*_bucket)[i];
193 193

	
194 194
      if ((*_next)[i] == INVALID || !(*_active)[(*_next)[i]]) return;
195 195

	
196 196
      //unlace
197 197
      (*_prev)[(*_next)[i]] = (*_prev)[i];
198 198
      if ((*_prev)[i] != INVALID) {
199 199
        (*_next)[(*_prev)[i]] = (*_next)[i];
200 200
      } else {
201 201
        _first[bucket] = (*_next)[i];
202 202
      }
203 203
      //lace
204 204
      (*_prev)[i] = _last[bucket];
205 205
      (*_next)[_last[bucket]] = i;
206 206
      (*_next)[i] = INVALID;
207 207
      _last[bucket] = i;
208 208
    }
209 209

	
210 210
    void addItem(const Node& i, int bucket) {
211 211
      (*_bucket)[i] = bucket;
212 212
      if (_last[bucket] != INVALID) {
213 213
        (*_prev)[i] = _last[bucket];
214 214
        (*_next)[_last[bucket]] = i;
215 215
        (*_next)[i] = INVALID;
216 216
        _last[bucket] = i;
217 217
      } else {
218 218
        (*_prev)[i] = INVALID;
219 219
        _first[bucket] = i;
220 220
        (*_next)[i] = INVALID;
221 221
        _last[bucket] = i;
222 222
      }
223 223
    }
224 224

	
225 225
    void findMinCutOut() {
226 226

	
227 227
      for (NodeIt n(_graph); n != INVALID; ++n) {
228 228
        (*_excess)[n] = 0;
229
        (*_source_set)[n] = false;
229 230
      }
230 231

	
231 232
      for (ArcIt a(_graph); a != INVALID; ++a) {
232 233
        (*_flow)[a] = 0;
233 234
      }
234 235

	
235 236
      int bucket_num = 0;
236 237
      std::vector<Node> queue(_node_num);
237 238
      int qfirst = 0, qlast = 0, qsep = 0;
238 239

	
239 240
      {
240 241
        typename Digraph::template NodeMap<bool> reached(_graph, false);
241 242

	
242 243
        reached[_source] = true;
243 244
        bool first_set = true;
244 245

	
245 246
        for (NodeIt t(_graph); t != INVALID; ++t) {
246 247
          if (reached[t]) continue;
247 248
          _sets.push_front(std::list<int>());
248 249

	
249 250
          queue[qlast++] = t;
250 251
          reached[t] = true;
251 252

	
252 253
          while (qfirst != qlast) {
253 254
            if (qsep == qfirst) {
254 255
              ++bucket_num;
255 256
              _sets.front().push_front(bucket_num);
256 257
              _dormant[bucket_num] = !first_set;
257 258
              _first[bucket_num] = _last[bucket_num] = INVALID;
258 259
              qsep = qlast;
259 260
            }
260 261

	
261 262
            Node n = queue[qfirst++];
262 263
            addItem(n, bucket_num);
263 264

	
264 265
            for (InArcIt a(_graph, n); a != INVALID; ++a) {
265 266
              Node u = _graph.source(a);
266 267
              if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
267 268
                reached[u] = true;
268 269
                queue[qlast++] = u;
269 270
              }
270 271
            }
271 272
          }
272 273
          first_set = false;
273 274
        }
274 275

	
275 276
        ++bucket_num;
276 277
        (*_bucket)[_source] = 0;
277 278
        _dormant[0] = true;
278 279
      }
279 280
      (*_source_set)[_source] = true;
280 281

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

	
294 295
        if ((*_active)[target]) {
295 296
          deactivate(target);
296 297
        }
297 298

	
298 299
        _highest = _sets.back().begin();
299 300
        while (_highest != _sets.back().end() &&
300 301
               !(*_active)[_first[*_highest]]) {
301 302
          ++_highest;
302 303
        }
303 304
      }
304 305

	
305 306
      while (true) {
306 307
        while (_highest != _sets.back().end()) {
307 308
          Node n = _first[*_highest];
308 309
          Value excess = (*_excess)[n];
309 310
          int next_bucket = _node_num;
310 311

	
311 312
          int under_bucket;
312 313
          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
313 314
            under_bucket = -1;
314 315
          } else {
315 316
            under_bucket = *(++std::list<int>::iterator(_highest));
316 317
          }
317 318

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

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

	
366 367
        no_more_push:
367 368

	
368 369
          (*_excess)[n] = excess;
369 370

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

	
388 389
              std::list<std::list<int> >::iterator new_set =
389 390
                _sets.insert(--_sets.end(), std::list<int>());
390 391

	
391 392
              new_set->push_front(bucket_num);
392 393
              (*_bucket)[n] = bucket_num;
393 394
              _first[bucket_num] = _last[bucket_num] = n;
394 395
              (*_next)[n] = INVALID;
395 396
              (*_prev)[n] = INVALID;
396 397
              _dormant[bucket_num] = true;
397 398
              ++bucket_num;
398 399

	
399 400
              while (_highest != _sets.back().end() &&
400 401
                     !(*_active)[_first[*_highest]]) {
401 402
                ++_highest;
402 403
              }
403 404
            } else {
404 405
              _first[*_highest] = (*_next)[n];
405 406
              (*_prev)[(*_next)[n]] = INVALID;
406 407

	
407 408
              while (next_bucket != *_highest) {
408 409
                --_highest;
409 410
              }
410 411

	
411 412
              if (_highest == _sets.back().begin()) {
412 413
                _sets.back().push_front(bucket_num);
413 414
                _dormant[bucket_num] = false;
414 415
                _first[bucket_num] = _last[bucket_num] = INVALID;
415 416
                ++bucket_num;
416 417
              }
417 418
              --_highest;
418 419

	
419 420
              (*_bucket)[n] = *_highest;
420 421
              (*_next)[n] = _first[*_highest];
421 422
              if (_first[*_highest] != INVALID) {
422 423
                (*_prev)[_first[*_highest]] = n;
423 424
              } else {
424 425
                _last[*_highest] = n;
425 426
              }
426 427
              _first[*_highest] = n;
427 428
            }
428 429
          } else {
429 430

	
430 431
            deactivate(n);
431 432
            if (!(*_active)[_first[*_highest]]) {
432 433
              ++_highest;
433 434
              if (_highest != _sets.back().end() &&
434 435
                  !(*_active)[_first[*_highest]]) {
435 436
                _highest = _sets.back().end();
436 437
              }
437 438
            }
438 439
          }
439 440
        }
440 441

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

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

	
485 486
          (*_bucket)[target] = 0;
486 487

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

	
499 500
          for (InArcIt a(_graph, target); a != INVALID; ++a) {
500 501
            Value rem = (*_flow)[a];
501 502
            if (!_tolerance.positive(rem)) continue;
502 503
            Node v = _graph.source(a);
503 504
            if (!(*_active)[v] && !(*_source_set)[v]) {
504 505
              activate(v);
505 506
            }
506 507
            (*_excess)[v] += rem;
507 508
            (*_flow)[a] = 0;
508 509
          }
509 510

	
510 511
          target = new_target;
511 512
          if ((*_active)[target]) {
512 513
            deactivate(target);
513 514
          }
514 515

	
515 516
          _highest = _sets.back().begin();
516 517
          while (_highest != _sets.back().end() &&
517 518
                 !(*_active)[_first[*_highest]]) {
518 519
            ++_highest;
519 520
          }
520 521
        }
521 522
      }
522 523
    }
523 524

	
524 525
    void findMinCutIn() {
525 526

	
526 527
      for (NodeIt n(_graph); n != INVALID; ++n) {
527 528
        (*_excess)[n] = 0;
529
        (*_source_set)[n] = false;
528 530
      }
529 531

	
530 532
      for (ArcIt a(_graph); a != INVALID; ++a) {
531 533
        (*_flow)[a] = 0;
532 534
      }
533 535

	
534 536
      int bucket_num = 0;
535 537
      std::vector<Node> queue(_node_num);
536 538
      int qfirst = 0, qlast = 0, qsep = 0;
537 539

	
538 540
      {
539 541
        typename Digraph::template NodeMap<bool> reached(_graph, false);
540 542

	
541 543
        reached[_source] = true;
542 544

	
543 545
        bool first_set = true;
544 546

	
545 547
        for (NodeIt t(_graph); t != INVALID; ++t) {
546 548
          if (reached[t]) continue;
547 549
          _sets.push_front(std::list<int>());
548 550

	
549 551
          queue[qlast++] = t;
550 552
          reached[t] = true;
551 553

	
552 554
          while (qfirst != qlast) {
553 555
            if (qsep == qfirst) {
554 556
              ++bucket_num;
555 557
              _sets.front().push_front(bucket_num);
556 558
              _dormant[bucket_num] = !first_set;
557 559
              _first[bucket_num] = _last[bucket_num] = INVALID;
558 560
              qsep = qlast;
559 561
            }
560 562

	
561 563
            Node n = queue[qfirst++];
562 564
            addItem(n, bucket_num);
563 565

	
564 566
            for (OutArcIt a(_graph, n); a != INVALID; ++a) {
565 567
              Node u = _graph.target(a);
566 568
              if (!reached[u] && _tolerance.positive((*_capacity)[a])) {
567 569
                reached[u] = true;
568 570
                queue[qlast++] = u;
569 571
              }
570 572
            }
571 573
          }
572 574
          first_set = false;
573 575
        }
574 576

	
575 577
        ++bucket_num;
576 578
        (*_bucket)[_source] = 0;
577 579
        _dormant[0] = true;
578 580
      }
579 581
      (*_source_set)[_source] = true;
580 582

	
581 583
      Node target = _last[_sets.back().back()];
582 584
      {
583 585
        for (InArcIt a(_graph, _source); a != INVALID; ++a) {
584 586
          if (_tolerance.positive((*_capacity)[a])) {
585 587
            Node u = _graph.source(a);
586 588
            (*_flow)[a] = (*_capacity)[a];
587 589
            (*_excess)[u] += (*_capacity)[a];
588 590
            if (!(*_active)[u] && u != _source) {
589 591
              activate(u);
590 592
            }
591 593
          }
592 594
        }
593 595
        if ((*_active)[target]) {
594 596
          deactivate(target);
595 597
        }
596 598

	
597 599
        _highest = _sets.back().begin();
598 600
        while (_highest != _sets.back().end() &&
599 601
               !(*_active)[_first[*_highest]]) {
600 602
          ++_highest;
601 603
        }
602 604
      }
603 605

	
604 606

	
605 607
      while (true) {
606 608
        while (_highest != _sets.back().end()) {
607 609
          Node n = _first[*_highest];
608 610
          Value excess = (*_excess)[n];
609 611
          int next_bucket = _node_num;
610 612

	
611 613
          int under_bucket;
612 614
          if (++std::list<int>::iterator(_highest) == _sets.back().end()) {
613 615
            under_bucket = -1;
614 616
          } else {
615 617
            under_bucket = *(++std::list<int>::iterator(_highest));
616 618
          }
617 619

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

	
642 644
          for (OutArcIt a(_graph, n); a != INVALID; ++a) {
643 645
            Node v = _graph.target(a);
644 646
            if (_dormant[(*_bucket)[v]]) continue;
645 647
            Value rem = (*_flow)[a];
646 648
            if (!_tolerance.positive(rem)) continue;
647 649
            if ((*_bucket)[v] == under_bucket) {
648 650
              if (!(*_active)[v] && v != target) {
649 651
                activate(v);
650 652
              }
651 653
              if (!_tolerance.less(rem, excess)) {
652 654
                (*_flow)[a] -= excess;
653 655
                (*_excess)[v] += excess;
654 656
                excess = 0;
655 657
                goto no_more_push;
656 658
              } else {
657 659
                excess -= rem;
658 660
                (*_excess)[v] += rem;
659 661
                (*_flow)[a] = 0;
660 662
              }
661 663
            } else if (next_bucket > (*_bucket)[v]) {
662 664
              next_bucket = (*_bucket)[v];
663 665
            }
664 666
          }
665 667

	
666 668
        no_more_push:
667 669

	
668 670
          (*_excess)[n] = excess;
669 671

	
670 672
          if (excess != 0) {
671 673
            if ((*_next)[n] == INVALID) {
672 674
              typename std::list<std::list<int> >::iterator new_set =
673 675
                _sets.insert(--_sets.end(), std::list<int>());
674 676
              new_set->splice(new_set->end(), _sets.back(),
675 677
                              _sets.back().begin(), ++_highest);
676 678
              for (std::list<int>::iterator it = new_set->begin();
677 679
                   it != new_set->end(); ++it) {
678 680
                _dormant[*it] = true;
679 681
              }
680 682
              while (_highest != _sets.back().end() &&
681 683
                     !(*_active)[_first[*_highest]]) {
682 684
                ++_highest;
683 685
              }
684 686
            } else if (next_bucket == _node_num) {
685 687
              _first[(*_bucket)[n]] = (*_next)[n];
686 688
              (*_prev)[(*_next)[n]] = INVALID;
687 689

	
688 690
              std::list<std::list<int> >::iterator new_set =
689 691
                _sets.insert(--_sets.end(), std::list<int>());
690 692

	
691 693
              new_set->push_front(bucket_num);
692 694
              (*_bucket)[n] = bucket_num;
693 695
              _first[bucket_num] = _last[bucket_num] = n;
694 696
              (*_next)[n] = INVALID;
695 697
              (*_prev)[n] = INVALID;
696 698
              _dormant[bucket_num] = true;
697 699
              ++bucket_num;
698 700

	
699 701
              while (_highest != _sets.back().end() &&
700 702
                     !(*_active)[_first[*_highest]]) {
701 703
                ++_highest;
702 704
              }
703 705
            } else {
704 706
              _first[*_highest] = (*_next)[n];
705 707
              (*_prev)[(*_next)[n]] = INVALID;
706 708

	
707 709
              while (next_bucket != *_highest) {
708 710
                --_highest;
709 711
              }
710 712
              if (_highest == _sets.back().begin()) {
711 713
                _sets.back().push_front(bucket_num);
712 714
                _dormant[bucket_num] = false;
713 715
                _first[bucket_num] = _last[bucket_num] = INVALID;
714 716
                ++bucket_num;
715 717
              }
716 718
              --_highest;
717 719

	
718 720
              (*_bucket)[n] = *_highest;
719 721
              (*_next)[n] = _first[*_highest];
Ignore white space 384 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-2009
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
#include <sstream>
20 20

	
21 21
#include <lemon/smart_graph.h>
22 22
#include <lemon/adaptors.h>
23 23
#include <lemon/concepts/digraph.h>
24 24
#include <lemon/concepts/maps.h>
25 25
#include <lemon/lgf_reader.h>
26 26
#include <lemon/hao_orlin.h>
27 27

	
28 28
#include "test_tools.h"
29 29

	
30 30
using namespace lemon;
31 31
using namespace std;
32 32

	
33 33
const std::string lgf =
34 34
  "@nodes\n"
35 35
  "label\n"
36 36
  "0\n"
37 37
  "1\n"
38 38
  "2\n"
39 39
  "3\n"
40 40
  "4\n"
41 41
  "5\n"
42 42
  "@edges\n"
43 43
  "     cap1 cap2 cap3\n"
44 44
  "0 1  1    1    1   \n"
45 45
  "0 2  2    2    4   \n"
46 46
  "1 2  4    4    4   \n"
47 47
  "3 4  1    1    1   \n"
48 48
  "3 5  2    2    4   \n"
49 49
  "4 5  4    4    4   \n"
50 50
  "5 4  4    4    4   \n"
51 51
  "2 3  1    6    6   \n"
52 52
  "4 0  1    6    6   \n";
53 53

	
54 54
void checkHaoOrlinCompile()
55 55
{
56 56
  typedef int Value;
57 57
  typedef concepts::Digraph Digraph;
58 58

	
59 59
  typedef Digraph::Node Node;
60 60
  typedef Digraph::Arc Arc;
61 61
  typedef concepts::ReadMap<Arc, Value> CapMap;
62 62
  typedef concepts::WriteMap<Node, bool> CutMap;
63 63

	
64 64
  Digraph g;
65 65
  Node n;
66 66
  CapMap cap;
67 67
  CutMap cut;
68 68
  Value v;
69 69

	
70 70
  HaoOrlin<Digraph, CapMap> ho_test(g, cap);
71 71
  const HaoOrlin<Digraph, CapMap>&
72 72
    const_ho_test = ho_test;
73 73

	
74 74
  ho_test.init();
75 75
  ho_test.init(n);
76 76
  ho_test.calculateOut();
77 77
  ho_test.calculateIn();
78 78
  ho_test.run();
79 79
  ho_test.run(n);
80 80

	
81 81
  v = const_ho_test.minCutValue();
82 82
  v = const_ho_test.minCutMap(cut);
83 83
}
84 84

	
85 85
template <typename Graph, typename CapMap, typename CutMap>
86 86
typename CapMap::Value 
87 87
  cutValue(const Graph& graph, const CapMap& cap, const CutMap& cut)
88 88
{
89 89
  typename CapMap::Value sum = 0;
90 90
  for (typename Graph::ArcIt a(graph); a != INVALID; ++a) {
91 91
    if (cut[graph.source(a)] && !cut[graph.target(a)])
92 92
      sum += cap[a];
93 93
  }
94 94
  return sum;
95 95
}
96 96

	
97 97
int main() {
98 98
  SmartDigraph graph;
99 99
  SmartDigraph::ArcMap<int> cap1(graph), cap2(graph), cap3(graph);
100 100
  SmartDigraph::NodeMap<bool> cut(graph);
101 101

	
102 102
  istringstream input(lgf);
103 103
  digraphReader(graph, input)
104 104
    .arcMap("cap1", cap1)
105 105
    .arcMap("cap2", cap2)
106 106
    .arcMap("cap3", cap3)
107 107
    .run();
108 108

	
109 109
  {
110 110
    HaoOrlin<SmartDigraph> ho(graph, cap1);
111 111
    ho.run();
112 112
    ho.minCutMap(cut);
113 113
    
114
    // BUG: The cut value should be positive
115
    check(ho.minCutValue() == 0, "Wrong cut value");
116
    // BUG: It should work
117
    //check(ho.minCutValue() == cutValue(graph, cap1, cut), "Wrong cut value");
114
    check(ho.minCutValue() == 1, "Wrong cut value");
115
    check(ho.minCutValue() == cutValue(graph, cap1, cut), "Wrong cut value");
118 116
  }
119 117
  {
120 118
    HaoOrlin<SmartDigraph> ho(graph, cap2);
121 119
    ho.run();
122 120
    ho.minCutMap(cut);
123
    
124
    // BUG: The cut value should be positive
125
    check(ho.minCutValue() == 0, "Wrong cut value");
126
    // BUG: It should work
127
    //check(ho.minCutValue() == cutValue(graph, cap2, cut), "Wrong cut value");
121

	
122
    check(ho.minCutValue() == 1, "Wrong cut value");
123
    check(ho.minCutValue() == cutValue(graph, cap2, cut), "Wrong cut value");
128 124
  }
129 125
  {
130 126
    HaoOrlin<SmartDigraph> ho(graph, cap3);
131 127
    ho.run();
132 128
    ho.minCutMap(cut);
133 129
    
134
    // BUG: The cut value should be positive
135
    check(ho.minCutValue() == 0, "Wrong cut value");
136
    // BUG: It should work
137
    //check(ho.minCutValue() == cutValue(graph, cap3, cut), "Wrong cut value");
130
    check(ho.minCutValue() == 1, "Wrong cut value");
131
    check(ho.minCutValue() == cutValue(graph, cap3, cut), "Wrong cut value");
138 132
  }
139 133
  
140 134
  typedef Undirector<SmartDigraph> UGraph;
141 135
  UGraph ugraph(graph);
142 136
  
143 137
  {
144 138
    HaoOrlin<UGraph, SmartDigraph::ArcMap<int> > ho(ugraph, cap1);
145 139
    ho.run();
146 140
    ho.minCutMap(cut);
147 141
    
148
    // BUG: The cut value should be 2
149
    check(ho.minCutValue() == 1, "Wrong cut value");
150
    // BUG: It should work
151
    //check(ho.minCutValue() == cutValue(ugraph, cap1, cut), "Wrong cut value");
142
    check(ho.minCutValue() == 2, "Wrong cut value");
143
    check(ho.minCutValue() == cutValue(ugraph, cap1, cut), "Wrong cut value");
152 144
  }
153 145
  {
154 146
    HaoOrlin<UGraph, SmartDigraph::ArcMap<int> > ho(ugraph, cap2);
155 147
    ho.run();
156 148
    ho.minCutMap(cut);
157 149
    
158
    // TODO: Check this cut value
159
    check(ho.minCutValue() == 4, "Wrong cut value");
160
    // BUG: It should work
161
    //check(ho.minCutValue() == cutValue(ugraph, cap2, cut), "Wrong cut value");
150
    check(ho.minCutValue() == 5, "Wrong cut value");
151
    check(ho.minCutValue() == cutValue(ugraph, cap2, cut), "Wrong cut value");
162 152
  }
163 153
  {
164 154
    HaoOrlin<UGraph, SmartDigraph::ArcMap<int> > ho(ugraph, cap3);
165 155
    ho.run();
166 156
    ho.minCutMap(cut);
167 157
    
168
    // TODO: Check this cut value
169 158
    check(ho.minCutValue() == 5, "Wrong cut value");
170
    // BUG: It should work
171
    //check(ho.minCutValue() == cutValue(ugraph, cap3, cut), "Wrong cut value");
159
    check(ho.minCutValue() == cutValue(ugraph, cap3, cut), "Wrong cut value");
172 160
  }
173 161

	
174 162
  return 0;
175 163
}
0 comments (0 inline)