gravatar
alpar (Alpar Juttner)
alpar@cs.elte.hu
Get rid of exceptions in Preflow
0 1 0
default
1 file changed with 4 insertions and 14 deletions:
↑ Collapse diff ↑
Ignore white space 4503599627370496 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_PREFLOW_H
20 20
#define LEMON_PREFLOW_H
21 21

	
22
#include <lemon/error.h>
23 22
#include <lemon/tolerance.h>
24 23
#include <lemon/elevator.h>
25 24

	
26 25
/// \file
27 26
/// \ingroup max_flow
28 27
/// \brief Implementation of the preflow algorithm.
29 28

	
30 29
namespace lemon {
31 30

	
32 31
  /// \brief Default traits class of Preflow class.
33 32
  ///
34 33
  /// Default traits class of Preflow class.
35 34
  /// \param _Graph Digraph type.
36 35
  /// \param _CapacityMap Type of capacity map.
37 36
  template <typename _Graph, typename _CapacityMap>
38 37
  struct PreflowDefaultTraits {
39 38

	
40 39
    /// \brief The digraph type the algorithm runs on.
41 40
    typedef _Graph Digraph;
42 41

	
43 42
    /// \brief The type of the map that stores the arc capacities.
44 43
    ///
45 44
    /// The type of the map that stores the arc capacities.
46 45
    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
47 46
    typedef _CapacityMap CapacityMap;
48 47

	
49 48
    /// \brief The type of the length of the arcs.
50 49
    typedef typename CapacityMap::Value Value;
51 50

	
52 51
    /// \brief The map type that stores the flow values.
53 52
    ///
54 53
    /// The map type that stores the flow values.
55 54
    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
56 55
    typedef typename Digraph::template ArcMap<Value> FlowMap;
57 56

	
58 57
    /// \brief Instantiates a FlowMap.
59 58
    ///
60 59
    /// This function instantiates a \ref FlowMap.
61 60
    /// \param digraph The digraph, to which we would like to define
62 61
    /// the flow map.
63 62
    static FlowMap* createFlowMap(const Digraph& digraph) {
64 63
      return new FlowMap(digraph);
65 64
    }
66 65

	
67 66
    /// \brief The eleavator type used by Preflow algorithm.
68 67
    ///
69 68
    /// The elevator type used by Preflow algorithm.
70 69
    ///
71 70
    /// \sa Elevator
72 71
    /// \sa LinkedElevator
73 72
    typedef LinkedElevator<Digraph, typename Digraph::Node> Elevator;
74 73

	
75 74
    /// \brief Instantiates an Elevator.
76 75
    ///
77 76
    /// This function instantiates a \ref Elevator.
78 77
    /// \param digraph The digraph, to which we would like to define
79 78
    /// the elevator.
80 79
    /// \param max_level The maximum level of the elevator.
81 80
    static Elevator* createElevator(const Digraph& digraph, int max_level) {
82 81
      return new Elevator(digraph, max_level);
83 82
    }
84 83

	
85 84
    /// \brief The tolerance used by the algorithm
86 85
    ///
87 86
    /// The tolerance used by the algorithm to handle inexact computation.
88 87
    typedef lemon::Tolerance<Value> Tolerance;
89 88

	
90 89
  };
91 90

	
92 91

	
93 92
  /// \ingroup max_flow
94 93
  ///
95 94
  /// \brief %Preflow algorithms class.
96 95
  ///
97 96
  /// This class provides an implementation of the Goldberg's \e
98 97
  /// preflow \e algorithm producing a flow of maximum value in a
99 98
  /// digraph. The preflow algorithms are the fastest known max
100 99
  /// flow algorithms. The current implementation use a mixture of the
101 100
  /// \e "highest label" and the \e "bound decrease" heuristics.
102 101
  /// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$.
103 102
  ///
104 103
  /// The algorithm consists from two phases. After the first phase
105 104
  /// the maximal flow value and the minimum cut can be obtained. The
106 105
  /// second phase constructs the feasible maximum flow on each arc.
107 106
  ///
108 107
  /// \param _Graph The digraph type the algorithm runs on.
109 108
  /// \param _CapacityMap The flow map type.
110 109
  /// \param _Traits Traits class to set various data types used by
111 110
  /// the algorithm.  The default traits class is \ref
112 111
  /// PreflowDefaultTraits.  See \ref PreflowDefaultTraits for the
113 112
  /// documentation of a %Preflow traits class.
114 113
  ///
115 114
  ///\author Jacint Szabo and Balazs Dezso
116 115
#ifdef DOXYGEN
117 116
  template <typename _Graph, typename _CapacityMap, typename _Traits>
118 117
#else
119 118
  template <typename _Graph,
120 119
            typename _CapacityMap = typename _Graph::template ArcMap<int>,
121 120
            typename _Traits = PreflowDefaultTraits<_Graph, _CapacityMap> >
122 121
#endif
123 122
  class Preflow {
124 123
  public:
125 124

	
126 125
    typedef _Traits Traits;
127 126
    typedef typename Traits::Digraph Digraph;
128 127
    typedef typename Traits::CapacityMap CapacityMap;
129 128
    typedef typename Traits::Value Value;
130 129

	
131 130
    typedef typename Traits::FlowMap FlowMap;
132 131
    typedef typename Traits::Elevator Elevator;
133 132
    typedef typename Traits::Tolerance Tolerance;
134 133

	
135
    /// \brief \ref Exception for uninitialized parameters.
136
    ///
137
    /// This error represents problems in the initialization
138
    /// of the parameters of the algorithms.
139
    class UninitializedParameter : public lemon::Exception {
140
    public:
141
      virtual const char* what() const throw() {
142
        return "lemon::Preflow::UninitializedParameter";
143
      }
144
    };
145

	
146 134
  private:
147 135

	
148 136
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
149 137

	
150 138
    const Digraph& _graph;
151 139
    const CapacityMap* _capacity;
152 140

	
153 141
    int _node_num;
154 142

	
155 143
    Node _source, _target;
156 144

	
157 145
    FlowMap* _flow;
158 146
    bool _local_flow;
159 147

	
160 148
    Elevator* _level;
161 149
    bool _local_level;
162 150

	
163 151
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
164 152
    ExcessMap* _excess;
165 153

	
166 154
    Tolerance _tolerance;
167 155

	
168 156
    bool _phase;
169 157

	
170 158

	
171 159
    void createStructures() {
172 160
      _node_num = countNodes(_graph);
173 161

	
174 162
      if (!_flow) {
175 163
        _flow = Traits::createFlowMap(_graph);
176 164
        _local_flow = true;
177 165
      }
178 166
      if (!_level) {
179 167
        _level = Traits::createElevator(_graph, _node_num);
180 168
        _local_level = true;
181 169
      }
182 170
      if (!_excess) {
183 171
        _excess = new ExcessMap(_graph);
184 172
      }
185 173
    }
186 174

	
187 175
    void destroyStructures() {
188 176
      if (_local_flow) {
189 177
        delete _flow;
190 178
      }
191 179
      if (_local_level) {
192 180
        delete _level;
193 181
      }
194 182
      if (_excess) {
195 183
        delete _excess;
196 184
      }
197 185
    }
198 186

	
199 187
  public:
200 188

	
201 189
    typedef Preflow Create;
202 190

	
203 191
    ///\name Named template parameters
204 192

	
205 193
    ///@{
206 194

	
207 195
    template <typename _FlowMap>
208 196
    struct DefFlowMapTraits : public Traits {
209 197
      typedef _FlowMap FlowMap;
210 198
      static FlowMap *createFlowMap(const Digraph&) {
211
        throw UninitializedParameter();
199
        LEMON_ASSERT(false, "FlowMap is not initialized");
200
        return 0; // ignore warnings
212 201
      }
213 202
    };
214 203

	
215 204
    /// \brief \ref named-templ-param "Named parameter" for setting
216 205
    /// FlowMap type
217 206
    ///
218 207
    /// \ref named-templ-param "Named parameter" for setting FlowMap
219 208
    /// type
220 209
    template <typename _FlowMap>
221 210
    struct DefFlowMap
222 211
      : public Preflow<Digraph, CapacityMap, DefFlowMapTraits<_FlowMap> > {
223 212
      typedef Preflow<Digraph, CapacityMap,
224 213
                      DefFlowMapTraits<_FlowMap> > Create;
225 214
    };
226 215

	
227 216
    template <typename _Elevator>
228 217
    struct DefElevatorTraits : public Traits {
229 218
      typedef _Elevator Elevator;
230 219
      static Elevator *createElevator(const Digraph&, int) {
231
        throw UninitializedParameter();
220
        LEMON_ASSERT(false, "Elevator is not initialized");
221
        return 0; // ignore warnings
232 222
      }
233 223
    };
234 224

	
235 225
    /// \brief \ref named-templ-param "Named parameter" for setting
236 226
    /// Elevator type
237 227
    ///
238 228
    /// \ref named-templ-param "Named parameter" for setting Elevator
239 229
    /// type
240 230
    template <typename _Elevator>
241 231
    struct DefElevator
242 232
      : public Preflow<Digraph, CapacityMap, DefElevatorTraits<_Elevator> > {
243 233
      typedef Preflow<Digraph, CapacityMap,
244 234
                      DefElevatorTraits<_Elevator> > Create;
245 235
    };
246 236

	
247 237
    template <typename _Elevator>
248 238
    struct DefStandardElevatorTraits : public Traits {
249 239
      typedef _Elevator Elevator;
250 240
      static Elevator *createElevator(const Digraph& digraph, int max_level) {
251 241
        return new Elevator(digraph, max_level);
252 242
      }
253 243
    };
254 244

	
255 245
    /// \brief \ref named-templ-param "Named parameter" for setting
256 246
    /// Elevator type
257 247
    ///
258 248
    /// \ref named-templ-param "Named parameter" for setting Elevator
259 249
    /// type. The Elevator should be standard constructor interface, ie.
260 250
    /// the digraph and the maximum level should be passed to it.
261 251
    template <typename _Elevator>
262 252
    struct DefStandardElevator
263 253
      : public Preflow<Digraph, CapacityMap,
264 254
                       DefStandardElevatorTraits<_Elevator> > {
265 255
      typedef Preflow<Digraph, CapacityMap,
266 256
                      DefStandardElevatorTraits<_Elevator> > Create;
267 257
    };
268 258

	
269 259
    /// @}
270 260

	
271 261
  protected:
272 262

	
273 263
    Preflow() {}
274 264

	
275 265
  public:
276 266

	
277 267

	
278 268
    /// \brief The constructor of the class.
279 269
    ///
280 270
    /// The constructor of the class.
281 271
    /// \param digraph The digraph the algorithm runs on.
282 272
    /// \param capacity The capacity of the arcs.
283 273
    /// \param source The source node.
284 274
    /// \param target The target node.
285 275
    Preflow(const Digraph& digraph, const CapacityMap& capacity,
286 276
               Node source, Node target)
287 277
      : _graph(digraph), _capacity(&capacity),
288 278
        _node_num(0), _source(source), _target(target),
289 279
        _flow(0), _local_flow(false),
290 280
        _level(0), _local_level(false),
291 281
        _excess(0), _tolerance(), _phase() {}
292 282

	
293 283
    /// \brief Destrcutor.
294 284
    ///
295 285
    /// Destructor.
296 286
    ~Preflow() {
297 287
      destroyStructures();
298 288
    }
299 289

	
300 290
    /// \brief Sets the capacity map.
301 291
    ///
302 292
    /// Sets the capacity map.
303 293
    /// \return \c (*this)
304 294
    Preflow& capacityMap(const CapacityMap& map) {
305 295
      _capacity = &map;
306 296
      return *this;
307 297
    }
308 298

	
309 299
    /// \brief Sets the flow map.
310 300
    ///
311 301
    /// Sets the flow map.
312 302
    /// \return \c (*this)
313 303
    Preflow& flowMap(FlowMap& map) {
314 304
      if (_local_flow) {
315 305
        delete _flow;
316 306
        _local_flow = false;
317 307
      }
318 308
      _flow = &map;
319 309
      return *this;
320 310
    }
321 311

	
322 312
    /// \brief Returns the flow map.
323 313
    ///
324 314
    /// \return The flow map.
325 315
    const FlowMap& flowMap() {
326 316
      return *_flow;
327 317
    }
328 318

	
329 319
    /// \brief Sets the elevator.
330 320
    ///
331 321
    /// Sets the elevator.
332 322
    /// \return \c (*this)
333 323
    Preflow& elevator(Elevator& elevator) {
334 324
      if (_local_level) {
335 325
        delete _level;
336 326
        _local_level = false;
337 327
      }
338 328
      _level = &elevator;
339 329
      return *this;
340 330
    }
341 331

	
342 332
    /// \brief Returns the elevator.
343 333
    ///
344 334
    /// \return The elevator.
345 335
    const Elevator& elevator() {
346 336
      return *_level;
347 337
    }
348 338

	
349 339
    /// \brief Sets the source node.
350 340
    ///
351 341
    /// Sets the source node.
352 342
    /// \return \c (*this)
353 343
    Preflow& source(const Node& node) {
354 344
      _source = node;
355 345
      return *this;
356 346
    }
357 347

	
358 348
    /// \brief Sets the target node.
359 349
    ///
360 350
    /// Sets the target node.
361 351
    /// \return \c (*this)
362 352
    Preflow& target(const Node& node) {
363 353
      _target = node;
364 354
      return *this;
365 355
    }
366 356

	
367 357
    /// \brief Sets the tolerance used by algorithm.
368 358
    ///
369 359
    /// Sets the tolerance used by algorithm.
370 360
    Preflow& tolerance(const Tolerance& tolerance) const {
371 361
      _tolerance = tolerance;
372 362
      return *this;
373 363
    }
374 364

	
375 365
    /// \brief Returns the tolerance used by algorithm.
376 366
    ///
377 367
    /// Returns the tolerance used by algorithm.
378 368
    const Tolerance& tolerance() const {
379 369
      return tolerance;
380 370
    }
381 371

	
382 372
    /// \name Execution control The simplest way to execute the
383 373
    /// algorithm is to use one of the member functions called \c
384 374
    /// run().
385 375
    /// \n
386 376
    /// If you need more control on initial solution or
387 377
    /// execution then you have to call one \ref init() function and then
388 378
    /// the startFirstPhase() and if you need the startSecondPhase().
389 379

	
390 380
    ///@{
391 381

	
392 382
    /// \brief Initializes the internal data structures.
393 383
    ///
394 384
    /// Initializes the internal data structures.
395 385
    ///
396 386
    void init() {
397 387
      createStructures();
398 388

	
399 389
      _phase = true;
400 390
      for (NodeIt n(_graph); n != INVALID; ++n) {
401 391
        _excess->set(n, 0);
402 392
      }
403 393

	
404 394
      for (ArcIt e(_graph); e != INVALID; ++e) {
405 395
        _flow->set(e, 0);
406 396
      }
407 397

	
408 398
      typename Digraph::template NodeMap<bool> reached(_graph, false);
409 399

	
410 400
      _level->initStart();
411 401
      _level->initAddItem(_target);
412 402

	
413 403
      std::vector<Node> queue;
414 404
      reached.set(_source, true);
415 405

	
416 406
      queue.push_back(_target);
417 407
      reached.set(_target, true);
418 408
      while (!queue.empty()) {
419 409
        _level->initNewLevel();
420 410
        std::vector<Node> nqueue;
421 411
        for (int i = 0; i < int(queue.size()); ++i) {
422 412
          Node n = queue[i];
423 413
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
424 414
            Node u = _graph.source(e);
425 415
            if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
426 416
              reached.set(u, true);
427 417
              _level->initAddItem(u);
428 418
              nqueue.push_back(u);
429 419
            }
430 420
          }
431 421
        }
432 422
        queue.swap(nqueue);
433 423
      }
434 424
      _level->initFinish();
435 425

	
436 426
      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
437 427
        if (_tolerance.positive((*_capacity)[e])) {
438 428
          Node u = _graph.target(e);
439 429
          if ((*_level)[u] == _level->maxLevel()) continue;
440 430
          _flow->set(e, (*_capacity)[e]);
441 431
          _excess->set(u, (*_excess)[u] + (*_capacity)[e]);
442 432
          if (u != _target && !_level->active(u)) {
443 433
            _level->activate(u);
444 434
          }
445 435
        }
446 436
      }
447 437
    }
448 438

	
449 439
    /// \brief Initializes the internal data structures.
450 440
    ///
451 441
    /// Initializes the internal data structures and sets the initial
452 442
    /// flow to the given \c flowMap. The \c flowMap should contain a
453 443
    /// flow or at least a preflow, ie. in each node excluding the
454 444
    /// target the incoming flow should greater or equal to the
455 445
    /// outgoing flow.
456 446
    /// \return %False when the given \c flowMap is not a preflow.
457 447
    template <typename FlowMap>
458 448
    bool flowInit(const FlowMap& flowMap) {
459 449
      createStructures();
460 450

	
461 451
      for (ArcIt e(_graph); e != INVALID; ++e) {
462 452
        _flow->set(e, flowMap[e]);
463 453
      }
464 454

	
465 455
      for (NodeIt n(_graph); n != INVALID; ++n) {
466 456
        Value excess = 0;
467 457
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
468 458
          excess += (*_flow)[e];
469 459
        }
470 460
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
471 461
          excess -= (*_flow)[e];
472 462
        }
473 463
        if (excess < 0 && n != _source) return false;
474 464
        _excess->set(n, excess);
475 465
      }
476 466

	
477 467
      typename Digraph::template NodeMap<bool> reached(_graph, false);
478 468

	
479 469
      _level->initStart();
480 470
      _level->initAddItem(_target);
481 471

	
482 472
      std::vector<Node> queue;
483 473
      reached.set(_source, true);
484 474

	
485 475
      queue.push_back(_target);
486 476
      reached.set(_target, true);
487 477
      while (!queue.empty()) {
488 478
        _level->initNewLevel();
489 479
        std::vector<Node> nqueue;
490 480
        for (int i = 0; i < int(queue.size()); ++i) {
491 481
          Node n = queue[i];
492 482
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
493 483
            Node u = _graph.source(e);
494 484
            if (!reached[u] &&
495 485
                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
496 486
              reached.set(u, true);
497 487
              _level->initAddItem(u);
498 488
              nqueue.push_back(u);
499 489
            }
500 490
          }
501 491
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
502 492
            Node v = _graph.target(e);
503 493
            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
504 494
              reached.set(v, true);
505 495
              _level->initAddItem(v);
506 496
              nqueue.push_back(v);
507 497
            }
508 498
          }
509 499
        }
510 500
        queue.swap(nqueue);
511 501
      }
512 502
      _level->initFinish();
513 503

	
514 504
      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
515 505
        Value rem = (*_capacity)[e] - (*_flow)[e];
516 506
        if (_tolerance.positive(rem)) {
517 507
          Node u = _graph.target(e);
518 508
          if ((*_level)[u] == _level->maxLevel()) continue;
519 509
          _flow->set(e, (*_capacity)[e]);
520 510
          _excess->set(u, (*_excess)[u] + rem);
521 511
          if (u != _target && !_level->active(u)) {
522 512
            _level->activate(u);
523 513
          }
524 514
        }
525 515
      }
526 516
      for (InArcIt e(_graph, _source); e != INVALID; ++e) {
527 517
        Value rem = (*_flow)[e];
528 518
        if (_tolerance.positive(rem)) {
529 519
          Node v = _graph.source(e);
530 520
          if ((*_level)[v] == _level->maxLevel()) continue;
531 521
          _flow->set(e, 0);
532 522
          _excess->set(v, (*_excess)[v] + rem);
533 523
          if (v != _target && !_level->active(v)) {
534 524
            _level->activate(v);
535 525
          }
536 526
        }
537 527
      }
538 528
      return true;
539 529
    }
540 530

	
541 531
    /// \brief Starts the first phase of the preflow algorithm.
542 532
    ///
543 533
    /// The preflow algorithm consists of two phases, this method runs
544 534
    /// the first phase. After the first phase the maximum flow value
545 535
    /// and a minimum value cut can already be computed, although a
546 536
    /// maximum flow is not yet obtained. So after calling this method
547 537
    /// \ref flowValue() returns the value of a maximum flow and \ref
548 538
    /// minCut() returns a minimum cut.
549 539
    /// \pre One of the \ref init() functions should be called.
550 540
    void startFirstPhase() {
551 541
      _phase = true;
552 542

	
553 543
      Node n = _level->highestActive();
554 544
      int level = _level->highestActiveLevel();
555 545
      while (n != INVALID) {
556 546
        int num = _node_num;
557 547

	
558 548
        while (num > 0 && n != INVALID) {
559 549
          Value excess = (*_excess)[n];
560 550
          int new_level = _level->maxLevel();
561 551

	
562 552
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
563 553
            Value rem = (*_capacity)[e] - (*_flow)[e];
564 554
            if (!_tolerance.positive(rem)) continue;
565 555
            Node v = _graph.target(e);
566 556
            if ((*_level)[v] < level) {
567 557
              if (!_level->active(v) && v != _target) {
568 558
                _level->activate(v);
569 559
              }
570 560
              if (!_tolerance.less(rem, excess)) {
571 561
                _flow->set(e, (*_flow)[e] + excess);
572 562
                _excess->set(v, (*_excess)[v] + excess);
573 563
                excess = 0;
574 564
                goto no_more_push_1;
575 565
              } else {
576 566
                excess -= rem;
577 567
                _excess->set(v, (*_excess)[v] + rem);
578 568
                _flow->set(e, (*_capacity)[e]);
579 569
              }
580 570
            } else if (new_level > (*_level)[v]) {
581 571
              new_level = (*_level)[v];
582 572
            }
583 573
          }
584 574

	
585 575
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
586 576
            Value rem = (*_flow)[e];
587 577
            if (!_tolerance.positive(rem)) continue;
588 578
            Node v = _graph.source(e);
589 579
            if ((*_level)[v] < level) {
590 580
              if (!_level->active(v) && v != _target) {
591 581
                _level->activate(v);
592 582
              }
593 583
              if (!_tolerance.less(rem, excess)) {
594 584
                _flow->set(e, (*_flow)[e] - excess);
595 585
                _excess->set(v, (*_excess)[v] + excess);
596 586
                excess = 0;
597 587
                goto no_more_push_1;
598 588
              } else {
599 589
                excess -= rem;
600 590
                _excess->set(v, (*_excess)[v] + rem);
601 591
                _flow->set(e, 0);
602 592
              }
603 593
            } else if (new_level > (*_level)[v]) {
604 594
              new_level = (*_level)[v];
605 595
            }
606 596
          }
607 597

	
608 598
        no_more_push_1:
609 599

	
610 600
          _excess->set(n, excess);
611 601

	
612 602
          if (excess != 0) {
613 603
            if (new_level + 1 < _level->maxLevel()) {
614 604
              _level->liftHighestActive(new_level + 1);
615 605
            } else {
616 606
              _level->liftHighestActiveToTop();
617 607
            }
618 608
            if (_level->emptyLevel(level)) {
619 609
              _level->liftToTop(level);
620 610
            }
621 611
          } else {
622 612
            _level->deactivate(n);
623 613
          }
624 614

	
625 615
          n = _level->highestActive();
626 616
          level = _level->highestActiveLevel();
627 617
          --num;
628 618
        }
629 619

	
630 620
        num = _node_num * 20;
631 621
        while (num > 0 && n != INVALID) {
632 622
          Value excess = (*_excess)[n];
633 623
          int new_level = _level->maxLevel();
634 624

	
635 625
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
636 626
            Value rem = (*_capacity)[e] - (*_flow)[e];
637 627
            if (!_tolerance.positive(rem)) continue;
638 628
            Node v = _graph.target(e);
639 629
            if ((*_level)[v] < level) {
640 630
              if (!_level->active(v) && v != _target) {
641 631
                _level->activate(v);
642 632
              }
643 633
              if (!_tolerance.less(rem, excess)) {
644 634
                _flow->set(e, (*_flow)[e] + excess);
645 635
                _excess->set(v, (*_excess)[v] + excess);
646 636
                excess = 0;
647 637
                goto no_more_push_2;
648 638
              } else {
649 639
                excess -= rem;
650 640
                _excess->set(v, (*_excess)[v] + rem);
651 641
                _flow->set(e, (*_capacity)[e]);
652 642
              }
653 643
            } else if (new_level > (*_level)[v]) {
654 644
              new_level = (*_level)[v];
655 645
            }
656 646
          }
657 647

	
658 648
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
659 649
            Value rem = (*_flow)[e];
660 650
            if (!_tolerance.positive(rem)) continue;
661 651
            Node v = _graph.source(e);
662 652
            if ((*_level)[v] < level) {
663 653
              if (!_level->active(v) && v != _target) {
664 654
                _level->activate(v);
665 655
              }
666 656
              if (!_tolerance.less(rem, excess)) {
667 657
                _flow->set(e, (*_flow)[e] - excess);
668 658
                _excess->set(v, (*_excess)[v] + excess);
669 659
                excess = 0;
670 660
                goto no_more_push_2;
671 661
              } else {
672 662
                excess -= rem;
673 663
                _excess->set(v, (*_excess)[v] + rem);
674 664
                _flow->set(e, 0);
675 665
              }
676 666
            } else if (new_level > (*_level)[v]) {
677 667
              new_level = (*_level)[v];
678 668
            }
679 669
          }
680 670

	
681 671
        no_more_push_2:
682 672

	
683 673
          _excess->set(n, excess);
684 674

	
685 675
          if (excess != 0) {
686 676
            if (new_level + 1 < _level->maxLevel()) {
687 677
              _level->liftActiveOn(level, new_level + 1);
688 678
            } else {
689 679
              _level->liftActiveToTop(level);
690 680
            }
691 681
            if (_level->emptyLevel(level)) {
692 682
              _level->liftToTop(level);
693 683
            }
694 684
          } else {
695 685
            _level->deactivate(n);
696 686
          }
697 687

	
698 688
          while (level >= 0 && _level->activeFree(level)) {
699 689
            --level;
700 690
          }
701 691
          if (level == -1) {
702 692
            n = _level->highestActive();
703 693
            level = _level->highestActiveLevel();
704 694
          } else {
705 695
            n = _level->activeOn(level);
706 696
          }
707 697
          --num;
708 698
        }
709 699
      }
710 700
    }
711 701

	
712 702
    /// \brief Starts the second phase of the preflow algorithm.
713 703
    ///
714 704
    /// The preflow algorithm consists of two phases, this method runs
715 705
    /// the second phase. After calling \ref init() and \ref
716 706
    /// startFirstPhase() and then \ref startSecondPhase(), \ref
717 707
    /// flowMap() return a maximum flow, \ref flowValue() returns the
718 708
    /// value of a maximum flow, \ref minCut() returns a minimum cut
719 709
    /// \pre The \ref init() and startFirstPhase() functions should be
720 710
    /// called before.
721 711
    void startSecondPhase() {
722 712
      _phase = false;
723 713

	
724 714
      typename Digraph::template NodeMap<bool> reached(_graph);
725 715
      for (NodeIt n(_graph); n != INVALID; ++n) {
726 716
        reached.set(n, (*_level)[n] < _level->maxLevel());
727 717
      }
728 718

	
729 719
      _level->initStart();
730 720
      _level->initAddItem(_source);
731 721

	
732 722
      std::vector<Node> queue;
733 723
      queue.push_back(_source);
734 724
      reached.set(_source, true);
735 725

	
736 726
      while (!queue.empty()) {
737 727
        _level->initNewLevel();
738 728
        std::vector<Node> nqueue;
739 729
        for (int i = 0; i < int(queue.size()); ++i) {
740 730
          Node n = queue[i];
741 731
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
742 732
            Node v = _graph.target(e);
743 733
            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
744 734
              reached.set(v, true);
745 735
              _level->initAddItem(v);
746 736
              nqueue.push_back(v);
747 737
            }
748 738
          }
749 739
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
750 740
            Node u = _graph.source(e);
751 741
            if (!reached[u] &&
752 742
                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
753 743
              reached.set(u, true);
754 744
              _level->initAddItem(u);
755 745
              nqueue.push_back(u);
756 746
            }
757 747
          }
758 748
        }
759 749
        queue.swap(nqueue);
760 750
      }
761 751
      _level->initFinish();
762 752

	
763 753
      for (NodeIt n(_graph); n != INVALID; ++n) {
764 754
        if (!reached[n]) {
765 755
          _level->dirtyTopButOne(n);
766 756
        } else if ((*_excess)[n] > 0 && _target != n) {
767 757
          _level->activate(n);
768 758
        }
769 759
      }
770 760

	
771 761
      Node n;
772 762
      while ((n = _level->highestActive()) != INVALID) {
773 763
        Value excess = (*_excess)[n];
774 764
        int level = _level->highestActiveLevel();
775 765
        int new_level = _level->maxLevel();
776 766

	
777 767
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
778 768
          Value rem = (*_capacity)[e] - (*_flow)[e];
779 769
          if (!_tolerance.positive(rem)) continue;
780 770
          Node v = _graph.target(e);
781 771
          if ((*_level)[v] < level) {
782 772
            if (!_level->active(v) && v != _source) {
783 773
              _level->activate(v);
784 774
            }
785 775
            if (!_tolerance.less(rem, excess)) {
786 776
              _flow->set(e, (*_flow)[e] + excess);
787 777
              _excess->set(v, (*_excess)[v] + excess);
788 778
              excess = 0;
789 779
              goto no_more_push;
790 780
            } else {
791 781
              excess -= rem;
792 782
              _excess->set(v, (*_excess)[v] + rem);
793 783
              _flow->set(e, (*_capacity)[e]);
794 784
            }
795 785
          } else if (new_level > (*_level)[v]) {
796 786
            new_level = (*_level)[v];
797 787
          }
798 788
        }
799 789

	
800 790
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
801 791
          Value rem = (*_flow)[e];
802 792
          if (!_tolerance.positive(rem)) continue;
803 793
          Node v = _graph.source(e);
804 794
          if ((*_level)[v] < level) {
805 795
            if (!_level->active(v) && v != _source) {
806 796
              _level->activate(v);
807 797
            }
808 798
            if (!_tolerance.less(rem, excess)) {
809 799
              _flow->set(e, (*_flow)[e] - excess);
810 800
              _excess->set(v, (*_excess)[v] + excess);
811 801
              excess = 0;
812 802
              goto no_more_push;
813 803
            } else {
814 804
              excess -= rem;
815 805
              _excess->set(v, (*_excess)[v] + rem);
816 806
              _flow->set(e, 0);
817 807
            }
818 808
          } else if (new_level > (*_level)[v]) {
819 809
            new_level = (*_level)[v];
820 810
          }
821 811
        }
822 812

	
823 813
      no_more_push:
824 814

	
825 815
        _excess->set(n, excess);
826 816

	
827 817
        if (excess != 0) {
828 818
          if (new_level + 1 < _level->maxLevel()) {
829 819
            _level->liftHighestActive(new_level + 1);
830 820
          } else {
831 821
            // Calculation error
832 822
            _level->liftHighestActiveToTop();
833 823
          }
834 824
          if (_level->emptyLevel(level)) {
835 825
            // Calculation error
836 826
            _level->liftToTop(level);
837 827
          }
838 828
        } else {
839 829
          _level->deactivate(n);
840 830
        }
841 831

	
842 832
      }
843 833
    }
844 834

	
845 835
    /// \brief Runs the preflow algorithm.
846 836
    ///
847 837
    /// Runs the preflow algorithm.
848 838
    /// \note pf.run() is just a shortcut of the following code.
849 839
    /// \code
850 840
    ///   pf.init();
851 841
    ///   pf.startFirstPhase();
852 842
    ///   pf.startSecondPhase();
853 843
    /// \endcode
854 844
    void run() {
855 845
      init();
856 846
      startFirstPhase();
857 847
      startSecondPhase();
858 848
    }
859 849

	
860 850
    /// \brief Runs the preflow algorithm to compute the minimum cut.
861 851
    ///
862 852
    /// Runs the preflow algorithm to compute the minimum cut.
863 853
    /// \note pf.runMinCut() is just a shortcut of the following code.
864 854
    /// \code
865 855
    ///   pf.init();
866 856
    ///   pf.startFirstPhase();
867 857
    /// \endcode
868 858
    void runMinCut() {
869 859
      init();
870 860
      startFirstPhase();
871 861
    }
872 862

	
873 863
    /// @}
874 864

	
875 865
    /// \name Query Functions
876 866
    /// The result of the %Preflow algorithm can be obtained using these
877 867
    /// functions.\n
878 868
    /// Before the use of these functions,
879 869
    /// either run() or start() must be called.
880 870

	
881 871
    ///@{
882 872

	
883 873
    /// \brief Returns the value of the maximum flow.
884 874
    ///
885 875
    /// Returns the value of the maximum flow by returning the excess
886 876
    /// of the target node \c t. This value equals to the value of
887 877
    /// the maximum flow already after the first phase.
888 878
    Value flowValue() const {
889 879
      return (*_excess)[_target];
890 880
    }
891 881

	
892 882
    /// \brief Returns true when the node is on the source side of minimum cut.
893 883
    ///
894 884
    /// Returns true when the node is on the source side of minimum
895 885
    /// cut. This method can be called both after running \ref
896 886
    /// startFirstPhase() and \ref startSecondPhase().
897 887
    bool minCut(const Node& node) const {
898 888
      return ((*_level)[node] == _level->maxLevel()) == _phase;
899 889
    }
900 890

	
901 891
    /// \brief Returns a minimum value cut.
902 892
    ///
903 893
    /// Sets the \c cutMap to the characteristic vector of a minimum value
904 894
    /// cut. This method can be called both after running \ref
905 895
    /// startFirstPhase() and \ref startSecondPhase(). The result after second
906 896
    /// phase could be changed slightly if inexact computation is used.
907 897
    /// \pre The \c cutMap should be a bool-valued node-map.
908 898
    template <typename CutMap>
909 899
    void minCutMap(CutMap& cutMap) const {
910 900
      for (NodeIt n(_graph); n != INVALID; ++n) {
911 901
        cutMap.set(n, minCut(n));
912 902
      }
913 903
    }
914 904

	
915 905
    /// \brief Returns the flow on the arc.
916 906
    ///
917 907
    /// Sets the \c flowMap to the flow on the arcs. This method can
918 908
    /// be called after the second phase of algorithm.
919 909
    Value flow(const Arc& arc) const {
920 910
      return (*_flow)[arc];
921 911
    }
922 912

	
923 913
    /// @}
924 914
  };
925 915
}
926 916

	
927 917
#endif
0 comments (0 inline)