gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Bug fix in Preflow and Circulation (#307)
0 2 0
default
2 files changed with 4 insertions and 4 deletions:
↑ Collapse diff ↑
Ignore white space 768 line context
... ...
@@ -72,723 +72,723 @@
72 72
    /// The type of the map that stores the flow values.
73 73
    /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap"
74 74
    /// concept.
75 75
    typedef typename Digraph::template ArcMap<Value> FlowMap;
76 76

	
77 77
    /// \brief Instantiates a FlowMap.
78 78
    ///
79 79
    /// This function instantiates a \ref FlowMap.
80 80
    /// \param digraph The digraph for which we would like to define
81 81
    /// the flow map.
82 82
    static FlowMap* createFlowMap(const Digraph& digraph) {
83 83
      return new FlowMap(digraph);
84 84
    }
85 85

	
86 86
    /// \brief The elevator type used by the algorithm.
87 87
    ///
88 88
    /// The elevator type used by the algorithm.
89 89
    ///
90 90
    /// \sa Elevator
91 91
    /// \sa LinkedElevator
92 92
    typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
93 93

	
94 94
    /// \brief Instantiates an Elevator.
95 95
    ///
96 96
    /// This function instantiates an \ref Elevator.
97 97
    /// \param digraph The digraph for which we would like to define
98 98
    /// the elevator.
99 99
    /// \param max_level The maximum level of the elevator.
100 100
    static Elevator* createElevator(const Digraph& digraph, int max_level) {
101 101
      return new Elevator(digraph, max_level);
102 102
    }
103 103

	
104 104
    /// \brief The tolerance used by the algorithm
105 105
    ///
106 106
    /// The tolerance used by the algorithm to handle inexact computation.
107 107
    typedef lemon::Tolerance<Value> Tolerance;
108 108

	
109 109
  };
110 110

	
111 111
  /**
112 112
     \brief Push-relabel algorithm for the network circulation problem.
113 113

	
114 114
     \ingroup max_flow
115 115
     This class implements a push-relabel algorithm for the \e network
116 116
     \e circulation problem.
117 117
     It is to find a feasible circulation when lower and upper bounds
118 118
     are given for the flow values on the arcs and lower bounds are
119 119
     given for the difference between the outgoing and incoming flow
120 120
     at the nodes.
121 121

	
122 122
     The exact formulation of this problem is the following.
123 123
     Let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$
124 124
     \f$upper: A\rightarrow\mathbf{R}\cup\{\infty\}\f$ denote the lower and
125 125
     upper bounds on the arcs, for which \f$lower(uv) \leq upper(uv)\f$
126 126
     holds for all \f$uv\in A\f$, and \f$sup: V\rightarrow\mathbf{R}\f$
127 127
     denotes the signed supply values of the nodes.
128 128
     If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
129 129
     supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
130 130
     \f$-sup(u)\f$ demand.
131 131
     A feasible circulation is an \f$f: A\rightarrow\mathbf{R}\f$
132 132
     solution of the following problem.
133 133

	
134 134
     \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu)
135 135
     \geq sup(u) \quad \forall u\in V, \f]
136 136
     \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A. \f]
137 137
     
138 138
     The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
139 139
     zero or negative in order to have a feasible solution (since the sum
140 140
     of the expressions on the left-hand side of the inequalities is zero).
141 141
     It means that the total demand must be greater or equal to the total
142 142
     supply and all the supplies have to be carried out from the supply nodes,
143 143
     but there could be demands that are not satisfied.
144 144
     If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
145 145
     constraints have to be satisfied with equality, i.e. all demands
146 146
     have to be satisfied and all supplies have to be used.
147 147
     
148 148
     If you need the opposite inequalities in the supply/demand constraints
149 149
     (i.e. the total demand is less than the total supply and all the demands
150 150
     have to be satisfied while there could be supplies that are not used),
151 151
     then you could easily transform the problem to the above form by reversing
152 152
     the direction of the arcs and taking the negative of the supply values
153 153
     (e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
154 154

	
155 155
     This algorithm either calculates a feasible circulation, or provides
156 156
     a \ref barrier() "barrier", which prooves that a feasible soultion
157 157
     cannot exist.
158 158

	
159 159
     Note that this algorithm also provides a feasible solution for the
160 160
     \ref min_cost_flow "minimum cost flow problem".
161 161

	
162 162
     \tparam GR The type of the digraph the algorithm runs on.
163 163
     \tparam LM The type of the lower bound map. The default
164 164
     map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
165 165
     \tparam UM The type of the upper bound (capacity) map.
166 166
     The default map type is \c LM.
167 167
     \tparam SM The type of the supply map. The default map type is
168 168
     \ref concepts::Digraph::NodeMap "GR::NodeMap<UM::Value>".
169 169
  */
170 170
#ifdef DOXYGEN
171 171
template< typename GR,
172 172
          typename LM,
173 173
          typename UM,
174 174
          typename SM,
175 175
          typename TR >
176 176
#else
177 177
template< typename GR,
178 178
          typename LM = typename GR::template ArcMap<int>,
179 179
          typename UM = LM,
180 180
          typename SM = typename GR::template NodeMap<typename UM::Value>,
181 181
          typename TR = CirculationDefaultTraits<GR, LM, UM, SM> >
182 182
#endif
183 183
  class Circulation {
184 184
  public:
185 185

	
186 186
    ///The \ref CirculationDefaultTraits "traits class" of the algorithm.
187 187
    typedef TR Traits;
188 188
    ///The type of the digraph the algorithm runs on.
189 189
    typedef typename Traits::Digraph Digraph;
190 190
    ///The type of the flow and supply values.
191 191
    typedef typename Traits::Value Value;
192 192

	
193 193
    ///The type of the lower bound map.
194 194
    typedef typename Traits::LowerMap LowerMap;
195 195
    ///The type of the upper bound (capacity) map.
196 196
    typedef typename Traits::UpperMap UpperMap;
197 197
    ///The type of the supply map.
198 198
    typedef typename Traits::SupplyMap SupplyMap;
199 199
    ///The type of the flow map.
200 200
    typedef typename Traits::FlowMap FlowMap;
201 201

	
202 202
    ///The type of the elevator.
203 203
    typedef typename Traits::Elevator Elevator;
204 204
    ///The type of the tolerance.
205 205
    typedef typename Traits::Tolerance Tolerance;
206 206

	
207 207
  private:
208 208

	
209 209
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
210 210

	
211 211
    const Digraph &_g;
212 212
    int _node_num;
213 213

	
214 214
    const LowerMap *_lo;
215 215
    const UpperMap *_up;
216 216
    const SupplyMap *_supply;
217 217

	
218 218
    FlowMap *_flow;
219 219
    bool _local_flow;
220 220

	
221 221
    Elevator* _level;
222 222
    bool _local_level;
223 223

	
224 224
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
225 225
    ExcessMap* _excess;
226 226

	
227 227
    Tolerance _tol;
228 228
    int _el;
229 229

	
230 230
  public:
231 231

	
232 232
    typedef Circulation Create;
233 233

	
234 234
    ///\name Named Template Parameters
235 235

	
236 236
    ///@{
237 237

	
238 238
    template <typename T>
239 239
    struct SetFlowMapTraits : public Traits {
240 240
      typedef T FlowMap;
241 241
      static FlowMap *createFlowMap(const Digraph&) {
242 242
        LEMON_ASSERT(false, "FlowMap is not initialized");
243 243
        return 0; // ignore warnings
244 244
      }
245 245
    };
246 246

	
247 247
    /// \brief \ref named-templ-param "Named parameter" for setting
248 248
    /// FlowMap type
249 249
    ///
250 250
    /// \ref named-templ-param "Named parameter" for setting FlowMap
251 251
    /// type.
252 252
    template <typename T>
253 253
    struct SetFlowMap
254 254
      : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
255 255
                           SetFlowMapTraits<T> > {
256 256
      typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
257 257
                          SetFlowMapTraits<T> > Create;
258 258
    };
259 259

	
260 260
    template <typename T>
261 261
    struct SetElevatorTraits : public Traits {
262 262
      typedef T Elevator;
263 263
      static Elevator *createElevator(const Digraph&, int) {
264 264
        LEMON_ASSERT(false, "Elevator is not initialized");
265 265
        return 0; // ignore warnings
266 266
      }
267 267
    };
268 268

	
269 269
    /// \brief \ref named-templ-param "Named parameter" for setting
270 270
    /// Elevator type
271 271
    ///
272 272
    /// \ref named-templ-param "Named parameter" for setting Elevator
273 273
    /// type. If this named parameter is used, then an external
274 274
    /// elevator object must be passed to the algorithm using the
275 275
    /// \ref elevator(Elevator&) "elevator()" function before calling
276 276
    /// \ref run() or \ref init().
277 277
    /// \sa SetStandardElevator
278 278
    template <typename T>
279 279
    struct SetElevator
280 280
      : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
281 281
                           SetElevatorTraits<T> > {
282 282
      typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
283 283
                          SetElevatorTraits<T> > Create;
284 284
    };
285 285

	
286 286
    template <typename T>
287 287
    struct SetStandardElevatorTraits : public Traits {
288 288
      typedef T Elevator;
289 289
      static Elevator *createElevator(const Digraph& digraph, int max_level) {
290 290
        return new Elevator(digraph, max_level);
291 291
      }
292 292
    };
293 293

	
294 294
    /// \brief \ref named-templ-param "Named parameter" for setting
295 295
    /// Elevator type with automatic allocation
296 296
    ///
297 297
    /// \ref named-templ-param "Named parameter" for setting Elevator
298 298
    /// type with automatic allocation.
299 299
    /// The Elevator should have standard constructor interface to be
300 300
    /// able to automatically created by the algorithm (i.e. the
301 301
    /// digraph and the maximum level should be passed to it).
302 302
    /// However an external elevator object could also be passed to the
303 303
    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
304 304
    /// before calling \ref run() or \ref init().
305 305
    /// \sa SetElevator
306 306
    template <typename T>
307 307
    struct SetStandardElevator
308 308
      : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
309 309
                       SetStandardElevatorTraits<T> > {
310 310
      typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
311 311
                      SetStandardElevatorTraits<T> > Create;
312 312
    };
313 313

	
314 314
    /// @}
315 315

	
316 316
  protected:
317 317

	
318 318
    Circulation() {}
319 319

	
320 320
  public:
321 321

	
322 322
    /// Constructor.
323 323

	
324 324
    /// The constructor of the class.
325 325
    ///
326 326
    /// \param graph The digraph the algorithm runs on.
327 327
    /// \param lower The lower bounds for the flow values on the arcs.
328 328
    /// \param upper The upper bounds (capacities) for the flow values 
329 329
    /// on the arcs.
330 330
    /// \param supply The signed supply values of the nodes.
331 331
    Circulation(const Digraph &graph, const LowerMap &lower,
332 332
                const UpperMap &upper, const SupplyMap &supply)
333 333
      : _g(graph), _lo(&lower), _up(&upper), _supply(&supply),
334 334
        _flow(NULL), _local_flow(false), _level(NULL), _local_level(false),
335 335
        _excess(NULL) {}
336 336

	
337 337
    /// Destructor.
338 338
    ~Circulation() {
339 339
      destroyStructures();
340 340
    }
341 341

	
342 342

	
343 343
  private:
344 344

	
345 345
    bool checkBoundMaps() {
346 346
      for (ArcIt e(_g);e!=INVALID;++e) {
347 347
        if (_tol.less((*_up)[e], (*_lo)[e])) return false;
348 348
      }
349 349
      return true;
350 350
    }
351 351

	
352 352
    void createStructures() {
353 353
      _node_num = _el = countNodes(_g);
354 354

	
355 355
      if (!_flow) {
356 356
        _flow = Traits::createFlowMap(_g);
357 357
        _local_flow = true;
358 358
      }
359 359
      if (!_level) {
360 360
        _level = Traits::createElevator(_g, _node_num);
361 361
        _local_level = true;
362 362
      }
363 363
      if (!_excess) {
364 364
        _excess = new ExcessMap(_g);
365 365
      }
366 366
    }
367 367

	
368 368
    void destroyStructures() {
369 369
      if (_local_flow) {
370 370
        delete _flow;
371 371
      }
372 372
      if (_local_level) {
373 373
        delete _level;
374 374
      }
375 375
      if (_excess) {
376 376
        delete _excess;
377 377
      }
378 378
    }
379 379

	
380 380
  public:
381 381

	
382 382
    /// Sets the lower bound map.
383 383

	
384 384
    /// Sets the lower bound map.
385 385
    /// \return <tt>(*this)</tt>
386 386
    Circulation& lowerMap(const LowerMap& map) {
387 387
      _lo = &map;
388 388
      return *this;
389 389
    }
390 390

	
391 391
    /// Sets the upper bound (capacity) map.
392 392

	
393 393
    /// Sets the upper bound (capacity) map.
394 394
    /// \return <tt>(*this)</tt>
395 395
    Circulation& upperMap(const UpperMap& map) {
396 396
      _up = &map;
397 397
      return *this;
398 398
    }
399 399

	
400 400
    /// Sets the supply map.
401 401

	
402 402
    /// Sets the supply map.
403 403
    /// \return <tt>(*this)</tt>
404 404
    Circulation& supplyMap(const SupplyMap& map) {
405 405
      _supply = &map;
406 406
      return *this;
407 407
    }
408 408

	
409 409
    /// \brief Sets the flow map.
410 410
    ///
411 411
    /// Sets the flow map.
412 412
    /// If you don't use this function before calling \ref run() or
413 413
    /// \ref init(), an instance will be allocated automatically.
414 414
    /// The destructor deallocates this automatically allocated map,
415 415
    /// of course.
416 416
    /// \return <tt>(*this)</tt>
417 417
    Circulation& flowMap(FlowMap& map) {
418 418
      if (_local_flow) {
419 419
        delete _flow;
420 420
        _local_flow = false;
421 421
      }
422 422
      _flow = &map;
423 423
      return *this;
424 424
    }
425 425

	
426 426
    /// \brief Sets the elevator used by algorithm.
427 427
    ///
428 428
    /// Sets the elevator used by algorithm.
429 429
    /// If you don't use this function before calling \ref run() or
430 430
    /// \ref init(), an instance will be allocated automatically.
431 431
    /// The destructor deallocates this automatically allocated elevator,
432 432
    /// of course.
433 433
    /// \return <tt>(*this)</tt>
434 434
    Circulation& elevator(Elevator& elevator) {
435 435
      if (_local_level) {
436 436
        delete _level;
437 437
        _local_level = false;
438 438
      }
439 439
      _level = &elevator;
440 440
      return *this;
441 441
    }
442 442

	
443 443
    /// \brief Returns a const reference to the elevator.
444 444
    ///
445 445
    /// Returns a const reference to the elevator.
446 446
    ///
447 447
    /// \pre Either \ref run() or \ref init() must be called before
448 448
    /// using this function.
449 449
    const Elevator& elevator() const {
450 450
      return *_level;
451 451
    }
452 452

	
453 453
    /// \brief Sets the tolerance used by algorithm.
454 454
    ///
455 455
    /// Sets the tolerance used by algorithm.
456
    Circulation& tolerance(const Tolerance& tolerance) const {
456
    Circulation& tolerance(const Tolerance& tolerance) {
457 457
      _tol = tolerance;
458 458
      return *this;
459 459
    }
460 460

	
461 461
    /// \brief Returns a const reference to the tolerance.
462 462
    ///
463 463
    /// Returns a const reference to the tolerance.
464 464
    const Tolerance& tolerance() const {
465
      return tolerance;
465
      return _tol;
466 466
    }
467 467

	
468 468
    /// \name Execution Control
469 469
    /// The simplest way to execute the algorithm is to call \ref run().\n
470 470
    /// If you need more control on the initial solution or the execution,
471 471
    /// first you have to call one of the \ref init() functions, then
472 472
    /// the \ref start() function.
473 473

	
474 474
    ///@{
475 475

	
476 476
    /// Initializes the internal data structures.
477 477

	
478 478
    /// Initializes the internal data structures and sets all flow values
479 479
    /// to the lower bound.
480 480
    void init()
481 481
    {
482 482
      LEMON_DEBUG(checkBoundMaps(),
483 483
        "Upper bounds must be greater or equal to the lower bounds");
484 484

	
485 485
      createStructures();
486 486

	
487 487
      for(NodeIt n(_g);n!=INVALID;++n) {
488 488
        (*_excess)[n] = (*_supply)[n];
489 489
      }
490 490

	
491 491
      for (ArcIt e(_g);e!=INVALID;++e) {
492 492
        _flow->set(e, (*_lo)[e]);
493 493
        (*_excess)[_g.target(e)] += (*_flow)[e];
494 494
        (*_excess)[_g.source(e)] -= (*_flow)[e];
495 495
      }
496 496

	
497 497
      // global relabeling tested, but in general case it provides
498 498
      // worse performance for random digraphs
499 499
      _level->initStart();
500 500
      for(NodeIt n(_g);n!=INVALID;++n)
501 501
        _level->initAddItem(n);
502 502
      _level->initFinish();
503 503
      for(NodeIt n(_g);n!=INVALID;++n)
504 504
        if(_tol.positive((*_excess)[n]))
505 505
          _level->activate(n);
506 506
    }
507 507

	
508 508
    /// Initializes the internal data structures using a greedy approach.
509 509

	
510 510
    /// Initializes the internal data structures using a greedy approach
511 511
    /// to construct the initial solution.
512 512
    void greedyInit()
513 513
    {
514 514
      LEMON_DEBUG(checkBoundMaps(),
515 515
        "Upper bounds must be greater or equal to the lower bounds");
516 516

	
517 517
      createStructures();
518 518

	
519 519
      for(NodeIt n(_g);n!=INVALID;++n) {
520 520
        (*_excess)[n] = (*_supply)[n];
521 521
      }
522 522

	
523 523
      for (ArcIt e(_g);e!=INVALID;++e) {
524 524
        if (!_tol.less(-(*_excess)[_g.target(e)], (*_up)[e])) {
525 525
          _flow->set(e, (*_up)[e]);
526 526
          (*_excess)[_g.target(e)] += (*_up)[e];
527 527
          (*_excess)[_g.source(e)] -= (*_up)[e];
528 528
        } else if (_tol.less(-(*_excess)[_g.target(e)], (*_lo)[e])) {
529 529
          _flow->set(e, (*_lo)[e]);
530 530
          (*_excess)[_g.target(e)] += (*_lo)[e];
531 531
          (*_excess)[_g.source(e)] -= (*_lo)[e];
532 532
        } else {
533 533
          Value fc = -(*_excess)[_g.target(e)];
534 534
          _flow->set(e, fc);
535 535
          (*_excess)[_g.target(e)] = 0;
536 536
          (*_excess)[_g.source(e)] -= fc;
537 537
        }
538 538
      }
539 539

	
540 540
      _level->initStart();
541 541
      for(NodeIt n(_g);n!=INVALID;++n)
542 542
        _level->initAddItem(n);
543 543
      _level->initFinish();
544 544
      for(NodeIt n(_g);n!=INVALID;++n)
545 545
        if(_tol.positive((*_excess)[n]))
546 546
          _level->activate(n);
547 547
    }
548 548

	
549 549
    ///Executes the algorithm
550 550

	
551 551
    ///This function executes the algorithm.
552 552
    ///
553 553
    ///\return \c true if a feasible circulation is found.
554 554
    ///
555 555
    ///\sa barrier()
556 556
    ///\sa barrierMap()
557 557
    bool start()
558 558
    {
559 559

	
560 560
      Node act;
561 561
      Node bact=INVALID;
562 562
      Node last_activated=INVALID;
563 563
      while((act=_level->highestActive())!=INVALID) {
564 564
        int actlevel=(*_level)[act];
565 565
        int mlevel=_node_num;
566 566
        Value exc=(*_excess)[act];
567 567

	
568 568
        for(OutArcIt e(_g,act);e!=INVALID; ++e) {
569 569
          Node v = _g.target(e);
570 570
          Value fc=(*_up)[e]-(*_flow)[e];
571 571
          if(!_tol.positive(fc)) continue;
572 572
          if((*_level)[v]<actlevel) {
573 573
            if(!_tol.less(fc, exc)) {
574 574
              _flow->set(e, (*_flow)[e] + exc);
575 575
              (*_excess)[v] += exc;
576 576
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
577 577
                _level->activate(v);
578 578
              (*_excess)[act] = 0;
579 579
              _level->deactivate(act);
580 580
              goto next_l;
581 581
            }
582 582
            else {
583 583
              _flow->set(e, (*_up)[e]);
584 584
              (*_excess)[v] += fc;
585 585
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
586 586
                _level->activate(v);
587 587
              exc-=fc;
588 588
            }
589 589
          }
590 590
          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
591 591
        }
592 592
        for(InArcIt e(_g,act);e!=INVALID; ++e) {
593 593
          Node v = _g.source(e);
594 594
          Value fc=(*_flow)[e]-(*_lo)[e];
595 595
          if(!_tol.positive(fc)) continue;
596 596
          if((*_level)[v]<actlevel) {
597 597
            if(!_tol.less(fc, exc)) {
598 598
              _flow->set(e, (*_flow)[e] - exc);
599 599
              (*_excess)[v] += exc;
600 600
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
601 601
                _level->activate(v);
602 602
              (*_excess)[act] = 0;
603 603
              _level->deactivate(act);
604 604
              goto next_l;
605 605
            }
606 606
            else {
607 607
              _flow->set(e, (*_lo)[e]);
608 608
              (*_excess)[v] += fc;
609 609
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
610 610
                _level->activate(v);
611 611
              exc-=fc;
612 612
            }
613 613
          }
614 614
          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
615 615
        }
616 616

	
617 617
        (*_excess)[act] = exc;
618 618
        if(!_tol.positive(exc)) _level->deactivate(act);
619 619
        else if(mlevel==_node_num) {
620 620
          _level->liftHighestActiveToTop();
621 621
          _el = _node_num;
622 622
          return false;
623 623
        }
624 624
        else {
625 625
          _level->liftHighestActive(mlevel+1);
626 626
          if(_level->onLevel(actlevel)==0) {
627 627
            _el = actlevel;
628 628
            return false;
629 629
          }
630 630
        }
631 631
      next_l:
632 632
        ;
633 633
      }
634 634
      return true;
635 635
    }
636 636

	
637 637
    /// Runs the algorithm.
638 638

	
639 639
    /// This function runs the algorithm.
640 640
    ///
641 641
    /// \return \c true if a feasible circulation is found.
642 642
    ///
643 643
    /// \note Apart from the return value, c.run() is just a shortcut of
644 644
    /// the following code.
645 645
    /// \code
646 646
    ///   c.greedyInit();
647 647
    ///   c.start();
648 648
    /// \endcode
649 649
    bool run() {
650 650
      greedyInit();
651 651
      return start();
652 652
    }
653 653

	
654 654
    /// @}
655 655

	
656 656
    /// \name Query Functions
657 657
    /// The results of the circulation algorithm can be obtained using
658 658
    /// these functions.\n
659 659
    /// Either \ref run() or \ref start() should be called before
660 660
    /// using them.
661 661

	
662 662
    ///@{
663 663

	
664 664
    /// \brief Returns the flow value on the given arc.
665 665
    ///
666 666
    /// Returns the flow value on the given arc.
667 667
    ///
668 668
    /// \pre Either \ref run() or \ref init() must be called before
669 669
    /// using this function.
670 670
    Value flow(const Arc& arc) const {
671 671
      return (*_flow)[arc];
672 672
    }
673 673

	
674 674
    /// \brief Returns a const reference to the flow map.
675 675
    ///
676 676
    /// Returns a const reference to the arc map storing the found flow.
677 677
    ///
678 678
    /// \pre Either \ref run() or \ref init() must be called before
679 679
    /// using this function.
680 680
    const FlowMap& flowMap() const {
681 681
      return *_flow;
682 682
    }
683 683

	
684 684
    /**
685 685
       \brief Returns \c true if the given node is in a barrier.
686 686

	
687 687
       Barrier is a set \e B of nodes for which
688 688

	
689 689
       \f[ \sum_{uv\in A: u\in B} upper(uv) -
690 690
           \sum_{uv\in A: v\in B} lower(uv) < \sum_{v\in B} sup(v) \f]
691 691

	
692 692
       holds. The existence of a set with this property prooves that a
693 693
       feasible circualtion cannot exist.
694 694

	
695 695
       This function returns \c true if the given node is in the found
696 696
       barrier. If a feasible circulation is found, the function
697 697
       gives back \c false for every node.
698 698

	
699 699
       \pre Either \ref run() or \ref init() must be called before
700 700
       using this function.
701 701

	
702 702
       \sa barrierMap()
703 703
       \sa checkBarrier()
704 704
    */
705 705
    bool barrier(const Node& node) const
706 706
    {
707 707
      return (*_level)[node] >= _el;
708 708
    }
709 709

	
710 710
    /// \brief Gives back a barrier.
711 711
    ///
712 712
    /// This function sets \c bar to the characteristic vector of the
713 713
    /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
714 714
    /// node map with \c bool (or convertible) value type.
715 715
    ///
716 716
    /// If a feasible circulation is found, the function gives back an
717 717
    /// empty set, so \c bar[v] will be \c false for all nodes \c v.
718 718
    ///
719 719
    /// \note This function calls \ref barrier() for each node,
720 720
    /// so it runs in O(n) time.
721 721
    ///
722 722
    /// \pre Either \ref run() or \ref init() must be called before
723 723
    /// using this function.
724 724
    ///
725 725
    /// \sa barrier()
726 726
    /// \sa checkBarrier()
727 727
    template<class BarrierMap>
728 728
    void barrierMap(BarrierMap &bar) const
729 729
    {
730 730
      for(NodeIt n(_g);n!=INVALID;++n)
731 731
        bar.set(n, (*_level)[n] >= _el);
732 732
    }
733 733

	
734 734
    /// @}
735 735

	
736 736
    /// \name Checker Functions
737 737
    /// The feasibility of the results can be checked using
738 738
    /// these functions.\n
739 739
    /// Either \ref run() or \ref start() should be called before
740 740
    /// using them.
741 741

	
742 742
    ///@{
743 743

	
744 744
    ///Check if the found flow is a feasible circulation
745 745

	
746 746
    ///Check if the found flow is a feasible circulation,
747 747
    ///
748 748
    bool checkFlow() const {
749 749
      for(ArcIt e(_g);e!=INVALID;++e)
750 750
        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
751 751
      for(NodeIt n(_g);n!=INVALID;++n)
752 752
        {
753 753
          Value dif=-(*_supply)[n];
754 754
          for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
755 755
          for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
756 756
          if(_tol.negative(dif)) return false;
757 757
        }
758 758
      return true;
759 759
    }
760 760

	
761 761
    ///Check whether or not the last execution provides a barrier
762 762

	
763 763
    ///Check whether or not the last execution provides a barrier.
764 764
    ///\sa barrier()
765 765
    ///\sa barrierMap()
766 766
    bool checkBarrier() const
767 767
    {
768 768
      Value delta=0;
769 769
      Value inf_cap = std::numeric_limits<Value>::has_infinity ?
770 770
        std::numeric_limits<Value>::infinity() :
771 771
        std::numeric_limits<Value>::max();
772 772
      for(NodeIt n(_g);n!=INVALID;++n)
773 773
        if(barrier(n))
774 774
          delta-=(*_supply)[n];
775 775
      for(ArcIt e(_g);e!=INVALID;++e)
776 776
        {
777 777
          Node s=_g.source(e);
778 778
          Node t=_g.target(e);
779 779
          if(barrier(s)&&!barrier(t)) {
780 780
            if (_tol.less(inf_cap - (*_up)[e], delta)) return false;
781 781
            delta+=(*_up)[e];
782 782
          }
783 783
          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
784 784
        }
785 785
      return _tol.negative(delta);
786 786
    }
787 787

	
788 788
    /// @}
789 789

	
790 790
  };
791 791

	
792 792
}
793 793

	
794 794
#endif
Ignore white space 768 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
#ifndef LEMON_PREFLOW_H
20 20
#define LEMON_PREFLOW_H
21 21

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

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

	
29 29
namespace lemon {
30 30

	
31 31
  /// \brief Default traits class of Preflow class.
32 32
  ///
33 33
  /// Default traits class of Preflow class.
34 34
  /// \tparam GR Digraph type.
35 35
  /// \tparam CAP Capacity map type.
36 36
  template <typename GR, typename CAP>
37 37
  struct PreflowDefaultTraits {
38 38

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

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

	
48 48
    /// \brief The type of the flow values.
49 49
    typedef typename CapacityMap::Value Value;
50 50

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

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

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

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

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

	
89 89
  };
90 90

	
91 91

	
92 92
  /// \ingroup max_flow
93 93
  ///
94 94
  /// \brief %Preflow algorithm class.
95 95
  ///
96 96
  /// This class provides an implementation of Goldberg-Tarjan's \e preflow
97 97
  /// \e push-relabel algorithm producing a \ref max_flow
98 98
  /// "flow of maximum value" in a digraph.
99 99
  /// The preflow algorithms are the fastest known maximum
100 100
  /// flow algorithms. The current implementation use a mixture of the
101 101
  /// \e "highest label" and the \e "bound decrease" heuristics.
102 102
  /// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$.
103 103
  ///
104 104
  /// The algorithm consists of two phases. After the first phase
105 105
  /// the maximum flow value and the minimum cut is obtained. The
106 106
  /// second phase constructs a feasible maximum flow on each arc.
107 107
  ///
108 108
  /// \tparam GR The type of the digraph the algorithm runs on.
109 109
  /// \tparam CAP The type of the capacity map. The default map
110 110
  /// type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
111 111
#ifdef DOXYGEN
112 112
  template <typename GR, typename CAP, typename TR>
113 113
#else
114 114
  template <typename GR,
115 115
            typename CAP = typename GR::template ArcMap<int>,
116 116
            typename TR = PreflowDefaultTraits<GR, CAP> >
117 117
#endif
118 118
  class Preflow {
119 119
  public:
120 120

	
121 121
    ///The \ref PreflowDefaultTraits "traits class" of the algorithm.
122 122
    typedef TR Traits;
123 123
    ///The type of the digraph the algorithm runs on.
124 124
    typedef typename Traits::Digraph Digraph;
125 125
    ///The type of the capacity map.
126 126
    typedef typename Traits::CapacityMap CapacityMap;
127 127
    ///The type of the flow values.
128 128
    typedef typename Traits::Value Value;
129 129

	
130 130
    ///The type of the flow map.
131 131
    typedef typename Traits::FlowMap FlowMap;
132 132
    ///The type of the elevator.
133 133
    typedef typename Traits::Elevator Elevator;
134 134
    ///The type of the tolerance.
135 135
    typedef typename Traits::Tolerance Tolerance;
136 136

	
137 137
  private:
138 138

	
139 139
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
140 140

	
141 141
    const Digraph& _graph;
142 142
    const CapacityMap* _capacity;
143 143

	
144 144
    int _node_num;
145 145

	
146 146
    Node _source, _target;
147 147

	
148 148
    FlowMap* _flow;
149 149
    bool _local_flow;
150 150

	
151 151
    Elevator* _level;
152 152
    bool _local_level;
153 153

	
154 154
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
155 155
    ExcessMap* _excess;
156 156

	
157 157
    Tolerance _tolerance;
158 158

	
159 159
    bool _phase;
160 160

	
161 161

	
162 162
    void createStructures() {
163 163
      _node_num = countNodes(_graph);
164 164

	
165 165
      if (!_flow) {
166 166
        _flow = Traits::createFlowMap(_graph);
167 167
        _local_flow = true;
168 168
      }
169 169
      if (!_level) {
170 170
        _level = Traits::createElevator(_graph, _node_num);
171 171
        _local_level = true;
172 172
      }
173 173
      if (!_excess) {
174 174
        _excess = new ExcessMap(_graph);
175 175
      }
176 176
    }
177 177

	
178 178
    void destroyStructures() {
179 179
      if (_local_flow) {
180 180
        delete _flow;
181 181
      }
182 182
      if (_local_level) {
183 183
        delete _level;
184 184
      }
185 185
      if (_excess) {
186 186
        delete _excess;
187 187
      }
188 188
    }
189 189

	
190 190
  public:
191 191

	
192 192
    typedef Preflow Create;
193 193

	
194 194
    ///\name Named Template Parameters
195 195

	
196 196
    ///@{
197 197

	
198 198
    template <typename T>
199 199
    struct SetFlowMapTraits : public Traits {
200 200
      typedef T FlowMap;
201 201
      static FlowMap *createFlowMap(const Digraph&) {
202 202
        LEMON_ASSERT(false, "FlowMap is not initialized");
203 203
        return 0; // ignore warnings
204 204
      }
205 205
    };
206 206

	
207 207
    /// \brief \ref named-templ-param "Named parameter" for setting
208 208
    /// FlowMap type
209 209
    ///
210 210
    /// \ref named-templ-param "Named parameter" for setting FlowMap
211 211
    /// type.
212 212
    template <typename T>
213 213
    struct SetFlowMap
214 214
      : public Preflow<Digraph, CapacityMap, SetFlowMapTraits<T> > {
215 215
      typedef Preflow<Digraph, CapacityMap,
216 216
                      SetFlowMapTraits<T> > Create;
217 217
    };
218 218

	
219 219
    template <typename T>
220 220
    struct SetElevatorTraits : public Traits {
221 221
      typedef T Elevator;
222 222
      static Elevator *createElevator(const Digraph&, int) {
223 223
        LEMON_ASSERT(false, "Elevator is not initialized");
224 224
        return 0; // ignore warnings
225 225
      }
226 226
    };
227 227

	
228 228
    /// \brief \ref named-templ-param "Named parameter" for setting
229 229
    /// Elevator type
230 230
    ///
231 231
    /// \ref named-templ-param "Named parameter" for setting Elevator
232 232
    /// type. If this named parameter is used, then an external
233 233
    /// elevator object must be passed to the algorithm using the
234 234
    /// \ref elevator(Elevator&) "elevator()" function before calling
235 235
    /// \ref run() or \ref init().
236 236
    /// \sa SetStandardElevator
237 237
    template <typename T>
238 238
    struct SetElevator
239 239
      : public Preflow<Digraph, CapacityMap, SetElevatorTraits<T> > {
240 240
      typedef Preflow<Digraph, CapacityMap,
241 241
                      SetElevatorTraits<T> > Create;
242 242
    };
243 243

	
244 244
    template <typename T>
245 245
    struct SetStandardElevatorTraits : public Traits {
246 246
      typedef T Elevator;
247 247
      static Elevator *createElevator(const Digraph& digraph, int max_level) {
248 248
        return new Elevator(digraph, max_level);
249 249
      }
250 250
    };
251 251

	
252 252
    /// \brief \ref named-templ-param "Named parameter" for setting
253 253
    /// Elevator type with automatic allocation
254 254
    ///
255 255
    /// \ref named-templ-param "Named parameter" for setting Elevator
256 256
    /// type with automatic allocation.
257 257
    /// The Elevator should have standard constructor interface to be
258 258
    /// able to automatically created by the algorithm (i.e. the
259 259
    /// digraph and the maximum level should be passed to it).
260 260
    /// However an external elevator object could also be passed to the
261 261
    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
262 262
    /// before calling \ref run() or \ref init().
263 263
    /// \sa SetElevator
264 264
    template <typename T>
265 265
    struct SetStandardElevator
266 266
      : public Preflow<Digraph, CapacityMap,
267 267
                       SetStandardElevatorTraits<T> > {
268 268
      typedef Preflow<Digraph, CapacityMap,
269 269
                      SetStandardElevatorTraits<T> > Create;
270 270
    };
271 271

	
272 272
    /// @}
273 273

	
274 274
  protected:
275 275

	
276 276
    Preflow() {}
277 277

	
278 278
  public:
279 279

	
280 280

	
281 281
    /// \brief The constructor of the class.
282 282
    ///
283 283
    /// The constructor of the class.
284 284
    /// \param digraph The digraph the algorithm runs on.
285 285
    /// \param capacity The capacity of the arcs.
286 286
    /// \param source The source node.
287 287
    /// \param target The target node.
288 288
    Preflow(const Digraph& digraph, const CapacityMap& capacity,
289 289
            Node source, Node target)
290 290
      : _graph(digraph), _capacity(&capacity),
291 291
        _node_num(0), _source(source), _target(target),
292 292
        _flow(0), _local_flow(false),
293 293
        _level(0), _local_level(false),
294 294
        _excess(0), _tolerance(), _phase() {}
295 295

	
296 296
    /// \brief Destructor.
297 297
    ///
298 298
    /// Destructor.
299 299
    ~Preflow() {
300 300
      destroyStructures();
301 301
    }
302 302

	
303 303
    /// \brief Sets the capacity map.
304 304
    ///
305 305
    /// Sets the capacity map.
306 306
    /// \return <tt>(*this)</tt>
307 307
    Preflow& capacityMap(const CapacityMap& map) {
308 308
      _capacity = &map;
309 309
      return *this;
310 310
    }
311 311

	
312 312
    /// \brief Sets the flow map.
313 313
    ///
314 314
    /// Sets the flow map.
315 315
    /// If you don't use this function before calling \ref run() or
316 316
    /// \ref init(), an instance will be allocated automatically.
317 317
    /// The destructor deallocates this automatically allocated map,
318 318
    /// of course.
319 319
    /// \return <tt>(*this)</tt>
320 320
    Preflow& flowMap(FlowMap& map) {
321 321
      if (_local_flow) {
322 322
        delete _flow;
323 323
        _local_flow = false;
324 324
      }
325 325
      _flow = &map;
326 326
      return *this;
327 327
    }
328 328

	
329 329
    /// \brief Sets the source node.
330 330
    ///
331 331
    /// Sets the source node.
332 332
    /// \return <tt>(*this)</tt>
333 333
    Preflow& source(const Node& node) {
334 334
      _source = node;
335 335
      return *this;
336 336
    }
337 337

	
338 338
    /// \brief Sets the target node.
339 339
    ///
340 340
    /// Sets the target node.
341 341
    /// \return <tt>(*this)</tt>
342 342
    Preflow& target(const Node& node) {
343 343
      _target = node;
344 344
      return *this;
345 345
    }
346 346

	
347 347
    /// \brief Sets the elevator used by algorithm.
348 348
    ///
349 349
    /// Sets the elevator used by algorithm.
350 350
    /// If you don't use this function before calling \ref run() or
351 351
    /// \ref init(), an instance will be allocated automatically.
352 352
    /// The destructor deallocates this automatically allocated elevator,
353 353
    /// of course.
354 354
    /// \return <tt>(*this)</tt>
355 355
    Preflow& elevator(Elevator& elevator) {
356 356
      if (_local_level) {
357 357
        delete _level;
358 358
        _local_level = false;
359 359
      }
360 360
      _level = &elevator;
361 361
      return *this;
362 362
    }
363 363

	
364 364
    /// \brief Returns a const reference to the elevator.
365 365
    ///
366 366
    /// Returns a const reference to the elevator.
367 367
    ///
368 368
    /// \pre Either \ref run() or \ref init() must be called before
369 369
    /// using this function.
370 370
    const Elevator& elevator() const {
371 371
      return *_level;
372 372
    }
373 373

	
374 374
    /// \brief Sets the tolerance used by algorithm.
375 375
    ///
376 376
    /// Sets the tolerance used by algorithm.
377
    Preflow& tolerance(const Tolerance& tolerance) const {
377
    Preflow& tolerance(const Tolerance& tolerance) {
378 378
      _tolerance = tolerance;
379 379
      return *this;
380 380
    }
381 381

	
382 382
    /// \brief Returns a const reference to the tolerance.
383 383
    ///
384 384
    /// Returns a const reference to the tolerance.
385 385
    const Tolerance& tolerance() const {
386
      return tolerance;
386
      return _tolerance;
387 387
    }
388 388

	
389 389
    /// \name Execution Control
390 390
    /// The simplest way to execute the preflow algorithm is to use
391 391
    /// \ref run() or \ref runMinCut().\n
392 392
    /// If you need more control on the initial solution or the execution,
393 393
    /// first you have to call one of the \ref init() functions, then
394 394
    /// \ref startFirstPhase() and if you need it \ref startSecondPhase().
395 395

	
396 396
    ///@{
397 397

	
398 398
    /// \brief Initializes the internal data structures.
399 399
    ///
400 400
    /// Initializes the internal data structures and sets the initial
401 401
    /// flow to zero on each arc.
402 402
    void init() {
403 403
      createStructures();
404 404

	
405 405
      _phase = true;
406 406
      for (NodeIt n(_graph); n != INVALID; ++n) {
407 407
        (*_excess)[n] = 0;
408 408
      }
409 409

	
410 410
      for (ArcIt e(_graph); e != INVALID; ++e) {
411 411
        _flow->set(e, 0);
412 412
      }
413 413

	
414 414
      typename Digraph::template NodeMap<bool> reached(_graph, false);
415 415

	
416 416
      _level->initStart();
417 417
      _level->initAddItem(_target);
418 418

	
419 419
      std::vector<Node> queue;
420 420
      reached[_source] = true;
421 421

	
422 422
      queue.push_back(_target);
423 423
      reached[_target] = true;
424 424
      while (!queue.empty()) {
425 425
        _level->initNewLevel();
426 426
        std::vector<Node> nqueue;
427 427
        for (int i = 0; i < int(queue.size()); ++i) {
428 428
          Node n = queue[i];
429 429
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
430 430
            Node u = _graph.source(e);
431 431
            if (!reached[u] && _tolerance.positive((*_capacity)[e])) {
432 432
              reached[u] = true;
433 433
              _level->initAddItem(u);
434 434
              nqueue.push_back(u);
435 435
            }
436 436
          }
437 437
        }
438 438
        queue.swap(nqueue);
439 439
      }
440 440
      _level->initFinish();
441 441

	
442 442
      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
443 443
        if (_tolerance.positive((*_capacity)[e])) {
444 444
          Node u = _graph.target(e);
445 445
          if ((*_level)[u] == _level->maxLevel()) continue;
446 446
          _flow->set(e, (*_capacity)[e]);
447 447
          (*_excess)[u] += (*_capacity)[e];
448 448
          if (u != _target && !_level->active(u)) {
449 449
            _level->activate(u);
450 450
          }
451 451
        }
452 452
      }
453 453
    }
454 454

	
455 455
    /// \brief Initializes the internal data structures using the
456 456
    /// given flow map.
457 457
    ///
458 458
    /// Initializes the internal data structures and sets the initial
459 459
    /// flow to the given \c flowMap. The \c flowMap should contain a
460 460
    /// flow or at least a preflow, i.e. at each node excluding the
461 461
    /// source node the incoming flow should greater or equal to the
462 462
    /// outgoing flow.
463 463
    /// \return \c false if the given \c flowMap is not a preflow.
464 464
    template <typename FlowMap>
465 465
    bool init(const FlowMap& flowMap) {
466 466
      createStructures();
467 467

	
468 468
      for (ArcIt e(_graph); e != INVALID; ++e) {
469 469
        _flow->set(e, flowMap[e]);
470 470
      }
471 471

	
472 472
      for (NodeIt n(_graph); n != INVALID; ++n) {
473 473
        Value excess = 0;
474 474
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
475 475
          excess += (*_flow)[e];
476 476
        }
477 477
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
478 478
          excess -= (*_flow)[e];
479 479
        }
480 480
        if (excess < 0 && n != _source) return false;
481 481
        (*_excess)[n] = excess;
482 482
      }
483 483

	
484 484
      typename Digraph::template NodeMap<bool> reached(_graph, false);
485 485

	
486 486
      _level->initStart();
487 487
      _level->initAddItem(_target);
488 488

	
489 489
      std::vector<Node> queue;
490 490
      reached[_source] = true;
491 491

	
492 492
      queue.push_back(_target);
493 493
      reached[_target] = true;
494 494
      while (!queue.empty()) {
495 495
        _level->initNewLevel();
496 496
        std::vector<Node> nqueue;
497 497
        for (int i = 0; i < int(queue.size()); ++i) {
498 498
          Node n = queue[i];
499 499
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
500 500
            Node u = _graph.source(e);
501 501
            if (!reached[u] &&
502 502
                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
503 503
              reached[u] = true;
504 504
              _level->initAddItem(u);
505 505
              nqueue.push_back(u);
506 506
            }
507 507
          }
508 508
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
509 509
            Node v = _graph.target(e);
510 510
            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
511 511
              reached[v] = true;
512 512
              _level->initAddItem(v);
513 513
              nqueue.push_back(v);
514 514
            }
515 515
          }
516 516
        }
517 517
        queue.swap(nqueue);
518 518
      }
519 519
      _level->initFinish();
520 520

	
521 521
      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {
522 522
        Value rem = (*_capacity)[e] - (*_flow)[e];
523 523
        if (_tolerance.positive(rem)) {
524 524
          Node u = _graph.target(e);
525 525
          if ((*_level)[u] == _level->maxLevel()) continue;
526 526
          _flow->set(e, (*_capacity)[e]);
527 527
          (*_excess)[u] += rem;
528 528
          if (u != _target && !_level->active(u)) {
529 529
            _level->activate(u);
530 530
          }
531 531
        }
532 532
      }
533 533
      for (InArcIt e(_graph, _source); e != INVALID; ++e) {
534 534
        Value rem = (*_flow)[e];
535 535
        if (_tolerance.positive(rem)) {
536 536
          Node v = _graph.source(e);
537 537
          if ((*_level)[v] == _level->maxLevel()) continue;
538 538
          _flow->set(e, 0);
539 539
          (*_excess)[v] += rem;
540 540
          if (v != _target && !_level->active(v)) {
541 541
            _level->activate(v);
542 542
          }
543 543
        }
544 544
      }
545 545
      return true;
546 546
    }
547 547

	
548 548
    /// \brief Starts the first phase of the preflow algorithm.
549 549
    ///
550 550
    /// The preflow algorithm consists of two phases, this method runs
551 551
    /// the first phase. After the first phase the maximum flow value
552 552
    /// and a minimum value cut can already be computed, although a
553 553
    /// maximum flow is not yet obtained. So after calling this method
554 554
    /// \ref flowValue() returns the value of a maximum flow and \ref
555 555
    /// minCut() returns a minimum cut.
556 556
    /// \pre One of the \ref init() functions must be called before
557 557
    /// using this function.
558 558
    void startFirstPhase() {
559 559
      _phase = true;
560 560

	
561 561
      Node n = _level->highestActive();
562 562
      int level = _level->highestActiveLevel();
563 563
      while (n != INVALID) {
564 564
        int num = _node_num;
565 565

	
566 566
        while (num > 0 && n != INVALID) {
567 567
          Value excess = (*_excess)[n];
568 568
          int new_level = _level->maxLevel();
569 569

	
570 570
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
571 571
            Value rem = (*_capacity)[e] - (*_flow)[e];
572 572
            if (!_tolerance.positive(rem)) continue;
573 573
            Node v = _graph.target(e);
574 574
            if ((*_level)[v] < level) {
575 575
              if (!_level->active(v) && v != _target) {
576 576
                _level->activate(v);
577 577
              }
578 578
              if (!_tolerance.less(rem, excess)) {
579 579
                _flow->set(e, (*_flow)[e] + excess);
580 580
                (*_excess)[v] += excess;
581 581
                excess = 0;
582 582
                goto no_more_push_1;
583 583
              } else {
584 584
                excess -= rem;
585 585
                (*_excess)[v] += rem;
586 586
                _flow->set(e, (*_capacity)[e]);
587 587
              }
588 588
            } else if (new_level > (*_level)[v]) {
589 589
              new_level = (*_level)[v];
590 590
            }
591 591
          }
592 592

	
593 593
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
594 594
            Value rem = (*_flow)[e];
595 595
            if (!_tolerance.positive(rem)) continue;
596 596
            Node v = _graph.source(e);
597 597
            if ((*_level)[v] < level) {
598 598
              if (!_level->active(v) && v != _target) {
599 599
                _level->activate(v);
600 600
              }
601 601
              if (!_tolerance.less(rem, excess)) {
602 602
                _flow->set(e, (*_flow)[e] - excess);
603 603
                (*_excess)[v] += excess;
604 604
                excess = 0;
605 605
                goto no_more_push_1;
606 606
              } else {
607 607
                excess -= rem;
608 608
                (*_excess)[v] += rem;
609 609
                _flow->set(e, 0);
610 610
              }
611 611
            } else if (new_level > (*_level)[v]) {
612 612
              new_level = (*_level)[v];
613 613
            }
614 614
          }
615 615

	
616 616
        no_more_push_1:
617 617

	
618 618
          (*_excess)[n] = excess;
619 619

	
620 620
          if (excess != 0) {
621 621
            if (new_level + 1 < _level->maxLevel()) {
622 622
              _level->liftHighestActive(new_level + 1);
623 623
            } else {
624 624
              _level->liftHighestActiveToTop();
625 625
            }
626 626
            if (_level->emptyLevel(level)) {
627 627
              _level->liftToTop(level);
628 628
            }
629 629
          } else {
630 630
            _level->deactivate(n);
631 631
          }
632 632

	
633 633
          n = _level->highestActive();
634 634
          level = _level->highestActiveLevel();
635 635
          --num;
636 636
        }
637 637

	
638 638
        num = _node_num * 20;
639 639
        while (num > 0 && n != INVALID) {
640 640
          Value excess = (*_excess)[n];
641 641
          int new_level = _level->maxLevel();
642 642

	
643 643
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
644 644
            Value rem = (*_capacity)[e] - (*_flow)[e];
645 645
            if (!_tolerance.positive(rem)) continue;
646 646
            Node v = _graph.target(e);
647 647
            if ((*_level)[v] < level) {
648 648
              if (!_level->active(v) && v != _target) {
649 649
                _level->activate(v);
650 650
              }
651 651
              if (!_tolerance.less(rem, excess)) {
652 652
                _flow->set(e, (*_flow)[e] + excess);
653 653
                (*_excess)[v] += excess;
654 654
                excess = 0;
655 655
                goto no_more_push_2;
656 656
              } else {
657 657
                excess -= rem;
658 658
                (*_excess)[v] += rem;
659 659
                _flow->set(e, (*_capacity)[e]);
660 660
              }
661 661
            } else if (new_level > (*_level)[v]) {
662 662
              new_level = (*_level)[v];
663 663
            }
664 664
          }
665 665

	
666 666
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
667 667
            Value rem = (*_flow)[e];
668 668
            if (!_tolerance.positive(rem)) continue;
669 669
            Node v = _graph.source(e);
670 670
            if ((*_level)[v] < level) {
671 671
              if (!_level->active(v) && v != _target) {
672 672
                _level->activate(v);
673 673
              }
674 674
              if (!_tolerance.less(rem, excess)) {
675 675
                _flow->set(e, (*_flow)[e] - excess);
676 676
                (*_excess)[v] += excess;
677 677
                excess = 0;
678 678
                goto no_more_push_2;
679 679
              } else {
680 680
                excess -= rem;
681 681
                (*_excess)[v] += rem;
682 682
                _flow->set(e, 0);
683 683
              }
684 684
            } else if (new_level > (*_level)[v]) {
685 685
              new_level = (*_level)[v];
686 686
            }
687 687
          }
688 688

	
689 689
        no_more_push_2:
690 690

	
691 691
          (*_excess)[n] = excess;
692 692

	
693 693
          if (excess != 0) {
694 694
            if (new_level + 1 < _level->maxLevel()) {
695 695
              _level->liftActiveOn(level, new_level + 1);
696 696
            } else {
697 697
              _level->liftActiveToTop(level);
698 698
            }
699 699
            if (_level->emptyLevel(level)) {
700 700
              _level->liftToTop(level);
701 701
            }
702 702
          } else {
703 703
            _level->deactivate(n);
704 704
          }
705 705

	
706 706
          while (level >= 0 && _level->activeFree(level)) {
707 707
            --level;
708 708
          }
709 709
          if (level == -1) {
710 710
            n = _level->highestActive();
711 711
            level = _level->highestActiveLevel();
712 712
          } else {
713 713
            n = _level->activeOn(level);
714 714
          }
715 715
          --num;
716 716
        }
717 717
      }
718 718
    }
719 719

	
720 720
    /// \brief Starts the second phase of the preflow algorithm.
721 721
    ///
722 722
    /// The preflow algorithm consists of two phases, this method runs
723 723
    /// the second phase. After calling one of the \ref init() functions
724 724
    /// and \ref startFirstPhase() and then \ref startSecondPhase(),
725 725
    /// \ref flowMap() returns a maximum flow, \ref flowValue() returns the
726 726
    /// value of a maximum flow, \ref minCut() returns a minimum cut
727 727
    /// \pre One of the \ref init() functions and \ref startFirstPhase()
728 728
    /// must be called before using this function.
729 729
    void startSecondPhase() {
730 730
      _phase = false;
731 731

	
732 732
      typename Digraph::template NodeMap<bool> reached(_graph);
733 733
      for (NodeIt n(_graph); n != INVALID; ++n) {
734 734
        reached[n] = (*_level)[n] < _level->maxLevel();
735 735
      }
736 736

	
737 737
      _level->initStart();
738 738
      _level->initAddItem(_source);
739 739

	
740 740
      std::vector<Node> queue;
741 741
      queue.push_back(_source);
742 742
      reached[_source] = true;
743 743

	
744 744
      while (!queue.empty()) {
745 745
        _level->initNewLevel();
746 746
        std::vector<Node> nqueue;
747 747
        for (int i = 0; i < int(queue.size()); ++i) {
748 748
          Node n = queue[i];
749 749
          for (OutArcIt e(_graph, n); e != INVALID; ++e) {
750 750
            Node v = _graph.target(e);
751 751
            if (!reached[v] && _tolerance.positive((*_flow)[e])) {
752 752
              reached[v] = true;
753 753
              _level->initAddItem(v);
754 754
              nqueue.push_back(v);
755 755
            }
756 756
          }
757 757
          for (InArcIt e(_graph, n); e != INVALID; ++e) {
758 758
            Node u = _graph.source(e);
759 759
            if (!reached[u] &&
760 760
                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {
761 761
              reached[u] = true;
762 762
              _level->initAddItem(u);
763 763
              nqueue.push_back(u);
764 764
            }
765 765
          }
766 766
        }
767 767
        queue.swap(nqueue);
768 768
      }
769 769
      _level->initFinish();
770 770

	
0 comments (0 inline)