gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Test file + doc improvements (#307)
0 4 0
default
4 files changed with 21 insertions and 7 deletions:
↑ Collapse diff ↑
Ignore white space 6 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_CIRCULATION_H
20 20
#define LEMON_CIRCULATION_H
21 21

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

	
26 26
///\ingroup max_flow
27 27
///\file
28 28
///\brief Push-relabel algorithm for finding a feasible circulation.
29 29
///
30 30
namespace lemon {
31 31

	
32 32
  /// \brief Default traits class of Circulation class.
33 33
  ///
34 34
  /// Default traits class of Circulation class.
35 35
  ///
36 36
  /// \tparam GR Type of the digraph the algorithm runs on.
37 37
  /// \tparam LM The type of the lower bound map.
38 38
  /// \tparam UM The type of the upper bound (capacity) map.
39 39
  /// \tparam SM The type of the supply map.
40 40
  template <typename GR, typename LM,
41 41
            typename UM, typename SM>
42 42
  struct CirculationDefaultTraits {
43 43

	
44 44
    /// \brief The type of the digraph the algorithm runs on.
45 45
    typedef GR Digraph;
46 46

	
47 47
    /// \brief The type of the lower bound map.
48 48
    ///
49 49
    /// The type of the map that stores the lower bounds on the arcs.
50 50
    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
51 51
    typedef LM LowerMap;
52 52

	
53 53
    /// \brief The type of the upper bound (capacity) map.
54 54
    ///
55 55
    /// The type of the map that stores the upper bounds (capacities)
56 56
    /// on the arcs.
57 57
    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
58 58
    typedef UM UpperMap;
59 59

	
60 60
    /// \brief The type of supply map.
61 61
    ///
62 62
    /// The type of the map that stores the signed supply values of the 
63 63
    /// nodes. 
64 64
    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
65 65
    typedef SM SupplyMap;
66 66

	
67 67
    /// \brief The type of the flow and supply values.
68 68
    typedef typename SupplyMap::Value Value;
69 69

	
70 70
    /// \brief The type of the map that stores the flow values.
71 71
    ///
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
    /// \brief Sets the tolerance used by algorithm.
453
    /// \brief Sets the tolerance used by the algorithm.
454 454
    ///
455
    /// Sets the tolerance used by algorithm.
455
    /// Sets the tolerance object used by the algorithm.
456
    /// \return <tt>(*this)</tt>
456 457
    Circulation& tolerance(const Tolerance& tolerance) {
457 458
      _tol = tolerance;
458 459
      return *this;
459 460
    }
460 461

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

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

	
474 476
    ///@{
475 477

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

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

	
485 487
      createStructures();
486 488

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

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

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

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

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

	
517 519
      createStructures();
518 520

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

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

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

	
549 551
    ///Executes the algorithm
550 552

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

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

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

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

	
637 639
    /// Runs the algorithm.
638 640

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

	
654 656
    /// @}
655 657

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

	
662 664
    ///@{
663 665

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

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

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

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

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

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

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

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

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

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

	
734 736
    /// @}
735 737

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

	
742 744
    ///@{
743 745

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

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

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

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

	
788 790
    /// @}
789 791

	
790 792
  };
791 793

	
792 794
}
793 795

	
794 796
#endif
Ignore white space 3072 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
  /// flow algorithms. The current implementation use a mixture of the
100
  /// flow algorithms. The current implementation uses 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
    /// \brief Sets the tolerance used by algorithm.
374
    /// \brief Sets the tolerance used by the algorithm.
375 375
    ///
376
    /// Sets the tolerance used by algorithm.
376
    /// Sets the tolerance object used by the algorithm.
377
    /// \return <tt>(*this)</tt>
377 378
    Preflow& tolerance(const Tolerance& tolerance) {
378 379
      _tolerance = tolerance;
379 380
      return *this;
380 381
    }
381 382

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

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

	
396 398
    ///@{
397 399

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

	
616 618
        no_more_push_1:
617 619

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

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

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

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

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

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

	
689 691
        no_more_push_2:
690 692

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

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

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

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

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

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

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

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

	
771 773
      for (NodeIt n(_graph); n != INVALID; ++n) {
772 774
        if (!reached[n]) {
773 775
          _level->dirtyTopButOne(n);
774 776
        } else if ((*_excess)[n] > 0 && _target != n) {
775 777
          _level->activate(n);
776 778
        }
777 779
      }
778 780

	
779 781
      Node n;
780 782
      while ((n = _level->highestActive()) != INVALID) {
781 783
        Value excess = (*_excess)[n];
782 784
        int level = _level->highestActiveLevel();
783 785
        int new_level = _level->maxLevel();
784 786

	
785 787
        for (OutArcIt e(_graph, n); e != INVALID; ++e) {
786 788
          Value rem = (*_capacity)[e] - (*_flow)[e];
787 789
          if (!_tolerance.positive(rem)) continue;
788 790
          Node v = _graph.target(e);
789 791
          if ((*_level)[v] < level) {
790 792
            if (!_level->active(v) && v != _source) {
791 793
              _level->activate(v);
792 794
            }
793 795
            if (!_tolerance.less(rem, excess)) {
794 796
              _flow->set(e, (*_flow)[e] + excess);
795 797
              (*_excess)[v] += excess;
796 798
              excess = 0;
797 799
              goto no_more_push;
798 800
            } else {
799 801
              excess -= rem;
800 802
              (*_excess)[v] += rem;
801 803
              _flow->set(e, (*_capacity)[e]);
802 804
            }
803 805
          } else if (new_level > (*_level)[v]) {
804 806
            new_level = (*_level)[v];
805 807
          }
806 808
        }
807 809

	
808 810
        for (InArcIt e(_graph, n); e != INVALID; ++e) {
809 811
          Value rem = (*_flow)[e];
810 812
          if (!_tolerance.positive(rem)) continue;
811 813
          Node v = _graph.source(e);
812 814
          if ((*_level)[v] < level) {
813 815
            if (!_level->active(v) && v != _source) {
814 816
              _level->activate(v);
815 817
            }
816 818
            if (!_tolerance.less(rem, excess)) {
817 819
              _flow->set(e, (*_flow)[e] - excess);
818 820
              (*_excess)[v] += excess;
819 821
              excess = 0;
820 822
              goto no_more_push;
821 823
            } else {
822 824
              excess -= rem;
823 825
              (*_excess)[v] += rem;
824 826
              _flow->set(e, 0);
825 827
            }
826 828
          } else if (new_level > (*_level)[v]) {
827 829
            new_level = (*_level)[v];
828 830
          }
829 831
        }
830 832

	
831 833
      no_more_push:
832 834

	
833 835
        (*_excess)[n] = excess;
834 836

	
835 837
        if (excess != 0) {
836 838
          if (new_level + 1 < _level->maxLevel()) {
837 839
            _level->liftHighestActive(new_level + 1);
838 840
          } else {
839 841
            // Calculation error
840 842
            _level->liftHighestActiveToTop();
841 843
          }
842 844
          if (_level->emptyLevel(level)) {
843 845
            // Calculation error
844 846
            _level->liftToTop(level);
845 847
          }
846 848
        } else {
847 849
          _level->deactivate(n);
848 850
        }
849 851

	
850 852
      }
851 853
    }
852 854

	
853 855
    /// \brief Runs the preflow algorithm.
854 856
    ///
855 857
    /// Runs the preflow algorithm.
856 858
    /// \note pf.run() is just a shortcut of the following code.
857 859
    /// \code
858 860
    ///   pf.init();
859 861
    ///   pf.startFirstPhase();
860 862
    ///   pf.startSecondPhase();
861 863
    /// \endcode
862 864
    void run() {
863 865
      init();
864 866
      startFirstPhase();
865 867
      startSecondPhase();
866 868
    }
867 869

	
868 870
    /// \brief Runs the preflow algorithm to compute the minimum cut.
869 871
    ///
870 872
    /// Runs the preflow algorithm to compute the minimum cut.
871 873
    /// \note pf.runMinCut() is just a shortcut of the following code.
872 874
    /// \code
873 875
    ///   pf.init();
874 876
    ///   pf.startFirstPhase();
875 877
    /// \endcode
876 878
    void runMinCut() {
877 879
      init();
878 880
      startFirstPhase();
879 881
    }
880 882

	
881 883
    /// @}
882 884

	
883 885
    /// \name Query Functions
884 886
    /// The results of the preflow algorithm can be obtained using these
885 887
    /// functions.\n
886 888
    /// Either one of the \ref run() "run*()" functions or one of the
887 889
    /// \ref startFirstPhase() "start*()" functions should be called
888 890
    /// before using them.
889 891

	
890 892
    ///@{
891 893

	
892 894
    /// \brief Returns the value of the maximum flow.
893 895
    ///
894 896
    /// Returns the value of the maximum flow by returning the excess
895 897
    /// of the target node. This value equals to the value of
896 898
    /// the maximum flow already after the first phase of the algorithm.
897 899
    ///
898 900
    /// \pre Either \ref run() or \ref init() must be called before
899 901
    /// using this function.
900 902
    Value flowValue() const {
901 903
      return (*_excess)[_target];
902 904
    }
903 905

	
904 906
    /// \brief Returns the flow value on the given arc.
905 907
    ///
906 908
    /// Returns the flow value on the given arc. This method can
907 909
    /// be called after the second phase of the algorithm.
908 910
    ///
909 911
    /// \pre Either \ref run() or \ref init() must be called before
910 912
    /// using this function.
911 913
    Value flow(const Arc& arc) const {
912 914
      return (*_flow)[arc];
913 915
    }
914 916

	
915 917
    /// \brief Returns a const reference to the flow map.
916 918
    ///
917 919
    /// Returns a const reference to the arc map storing the found flow.
918 920
    /// This method can be called after the second phase of the algorithm.
919 921
    ///
920 922
    /// \pre Either \ref run() or \ref init() must be called before
921 923
    /// using this function.
922 924
    const FlowMap& flowMap() const {
923 925
      return *_flow;
924 926
    }
925 927

	
926 928
    /// \brief Returns \c true when the node is on the source side of the
927 929
    /// minimum cut.
928 930
    ///
929 931
    /// Returns true when the node is on the source side of the found
930 932
    /// minimum cut. This method can be called both after running \ref
931 933
    /// startFirstPhase() and \ref startSecondPhase().
932 934
    ///
933 935
    /// \pre Either \ref run() or \ref init() must be called before
934 936
    /// using this function.
935 937
    bool minCut(const Node& node) const {
936 938
      return ((*_level)[node] == _level->maxLevel()) == _phase;
937 939
    }
938 940

	
939 941
    /// \brief Gives back a minimum value cut.
940 942
    ///
941 943
    /// Sets \c cutMap to the characteristic vector of a minimum value
942 944
    /// cut. \c cutMap should be a \ref concepts::WriteMap "writable"
943 945
    /// node map with \c bool (or convertible) value type.
944 946
    ///
945 947
    /// This method can be called both after running \ref startFirstPhase()
946 948
    /// and \ref startSecondPhase(). The result after the second phase
947 949
    /// could be slightly different if inexact computation is used.
948 950
    ///
949 951
    /// \note This function calls \ref minCut() for each node, so it runs in
950 952
    /// O(n) time.
951 953
    ///
952 954
    /// \pre Either \ref run() or \ref init() must be called before
953 955
    /// using this function.
954 956
    template <typename CutMap>
955 957
    void minCutMap(CutMap& cutMap) const {
956 958
      for (NodeIt n(_graph); n != INVALID; ++n) {
957 959
        cutMap.set(n, minCut(n));
958 960
      }
959 961
    }
960 962

	
961 963
    /// @}
962 964
  };
963 965
}
964 966

	
965 967
#endif
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#include <iostream>
20 20

	
21 21
#include "test_tools.h"
22 22
#include <lemon/list_graph.h>
23 23
#include <lemon/circulation.h>
24 24
#include <lemon/lgf_reader.h>
25 25
#include <lemon/concepts/digraph.h>
26 26
#include <lemon/concepts/maps.h>
27 27

	
28 28
using namespace lemon;
29 29

	
30 30
char test_lgf[] =
31 31
  "@nodes\n"
32 32
  "label\n"
33 33
  "0\n"
34 34
  "1\n"
35 35
  "2\n"
36 36
  "3\n"
37 37
  "4\n"
38 38
  "5\n"
39 39
  "@arcs\n"
40 40
  "     lcap  ucap\n"
41 41
  "0 1  2  10\n"
42 42
  "0 2  2  6\n"
43 43
  "1 3  4  7\n"
44 44
  "1 4  0  5\n"
45 45
  "2 4  1  3\n"
46 46
  "3 5  3  8\n"
47 47
  "4 5  3  7\n"
48 48
  "@attributes\n"
49 49
  "source 0\n"
50 50
  "sink   5\n";
51 51

	
52 52
void checkCirculationCompile()
53 53
{
54 54
  typedef int VType;
55 55
  typedef concepts::Digraph Digraph;
56 56

	
57 57
  typedef Digraph::Node Node;
58 58
  typedef Digraph::Arc Arc;
59 59
  typedef concepts::ReadMap<Arc,VType> CapMap;
60 60
  typedef concepts::ReadMap<Node,VType> SupplyMap;
61 61
  typedef concepts::ReadWriteMap<Arc,VType> FlowMap;
62 62
  typedef concepts::WriteMap<Node,bool> BarrierMap;
63 63

	
64 64
  typedef Elevator<Digraph, Digraph::Node> Elev;
65 65
  typedef LinkedElevator<Digraph, Digraph::Node> LinkedElev;
66 66

	
67 67
  Digraph g;
68 68
  Node n;
69 69
  Arc a;
70 70
  CapMap lcap, ucap;
71 71
  SupplyMap supply;
72 72
  FlowMap flow;
73 73
  BarrierMap bar;
74 74
  VType v;
75 75
  bool b;
76 76

	
77 77
  typedef Circulation<Digraph, CapMap, CapMap, SupplyMap>
78 78
            ::SetFlowMap<FlowMap>
79 79
            ::SetElevator<Elev>
80 80
            ::SetStandardElevator<LinkedElev>
81 81
            ::Create CirculationType;
82 82
  CirculationType circ_test(g, lcap, ucap, supply);
83 83
  const CirculationType& const_circ_test = circ_test;
84 84
   
85 85
  circ_test
86 86
    .lowerMap(lcap)
87 87
    .upperMap(ucap)
88 88
    .supplyMap(supply)
89 89
    .flowMap(flow);
90
  
91
  const CirculationType::Elevator& elev = const_circ_test.elevator();
92
  circ_test.elevator(const_cast<CirculationType::Elevator&>(elev));
93
  CirculationType::Tolerance tol = const_circ_test.tolerance();
94
  circ_test.tolerance(tol);
90 95

	
91 96
  circ_test.init();
92 97
  circ_test.greedyInit();
93 98
  circ_test.start();
94 99
  circ_test.run();
95 100

	
96 101
  v = const_circ_test.flow(a);
97 102
  const FlowMap& fm = const_circ_test.flowMap();
98 103
  b = const_circ_test.barrier(n);
99 104
  const_circ_test.barrierMap(bar);
100 105
  
101 106
  ignore_unused_variable_warning(fm);
102 107
}
103 108

	
104 109
template <class G, class LM, class UM, class DM>
105 110
void checkCirculation(const G& g, const LM& lm, const UM& um,
106 111
                      const DM& dm, bool find)
107 112
{
108 113
  Circulation<G, LM, UM, DM> circ(g, lm, um, dm);
109 114
  bool ret = circ.run();
110 115
  if (find) {
111 116
    check(ret, "A feasible solution should have been found.");
112 117
    check(circ.checkFlow(), "The found flow is corrupt.");
113 118
    check(!circ.checkBarrier(), "A barrier should not have been found.");
114 119
  } else {
115 120
    check(!ret, "A feasible solution should not have been found.");
116 121
    check(circ.checkBarrier(), "The found barrier is corrupt.");
117 122
  }
118 123
}
119 124

	
120 125
int main (int, char*[])
121 126
{
122 127
  typedef ListDigraph Digraph;
123 128
  DIGRAPH_TYPEDEFS(Digraph);
124 129

	
125 130
  Digraph g;
126 131
  IntArcMap lo(g), up(g);
127 132
  IntNodeMap delta(g, 0);
128 133
  Node s, t;
129 134

	
130 135
  std::istringstream input(test_lgf);
131 136
  DigraphReader<Digraph>(g,input).
132 137
    arcMap("lcap", lo).
133 138
    arcMap("ucap", up).
134 139
    node("source",s).
135 140
    node("sink",t).
136 141
    run();
137 142

	
138 143
  delta[s] = 7; delta[t] = -7;
139 144
  checkCirculation(g, lo, up, delta, true);
140 145

	
141 146
  delta[s] = 13; delta[t] = -13;
142 147
  checkCirculation(g, lo, up, delta, true);
143 148

	
144 149
  delta[s] = 6; delta[t] = -6;
145 150
  checkCirculation(g, lo, up, delta, false);
146 151

	
147 152
  delta[s] = 14; delta[t] = -14;
148 153
  checkCirculation(g, lo, up, delta, false);
149 154

	
150 155
  delta[s] = 7; delta[t] = -13;
151 156
  checkCirculation(g, lo, up, delta, true);
152 157

	
153 158
  delta[s] = 5; delta[t] = -15;
154 159
  checkCirculation(g, lo, up, delta, true);
155 160

	
156 161
  delta[s] = 10; delta[t] = -11;
157 162
  checkCirculation(g, lo, up, delta, true);
158 163

	
159 164
  delta[s] = 11; delta[t] = -10;
160 165
  checkCirculation(g, lo, up, delta, false);
161 166

	
162 167
  return 0;
163 168
}
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#include <iostream>
20 20

	
21 21
#include "test_tools.h"
22 22
#include <lemon/smart_graph.h>
23 23
#include <lemon/preflow.h>
24 24
#include <lemon/concepts/digraph.h>
25 25
#include <lemon/concepts/maps.h>
26 26
#include <lemon/lgf_reader.h>
27 27
#include <lemon/elevator.h>
28 28

	
29 29
using namespace lemon;
30 30

	
31 31
char test_lgf[] =
32 32
  "@nodes\n"
33 33
  "label\n"
34 34
  "0\n"
35 35
  "1\n"
36 36
  "2\n"
37 37
  "3\n"
38 38
  "4\n"
39 39
  "5\n"
40 40
  "6\n"
41 41
  "7\n"
42 42
  "8\n"
43 43
  "9\n"
44 44
  "@arcs\n"
45 45
  "    label capacity\n"
46 46
  "0 1 0     20\n"
47 47
  "0 2 1     0\n"
48 48
  "1 1 2     3\n"
49 49
  "1 2 3     8\n"
50 50
  "1 3 4     8\n"
51 51
  "2 5 5     5\n"
52 52
  "3 2 6     5\n"
53 53
  "3 5 7     5\n"
54 54
  "3 6 8     5\n"
55 55
  "4 3 9     3\n"
56 56
  "5 7 10    3\n"
57 57
  "5 6 11    10\n"
58 58
  "5 8 12    10\n"
59 59
  "6 8 13    8\n"
60 60
  "8 9 14    20\n"
61 61
  "8 1 15    5\n"
62 62
  "9 5 16    5\n"
63 63
  "@attributes\n"
64 64
  "source 1\n"
65 65
  "target 8\n";
66 66

	
67 67
void checkPreflowCompile()
68 68
{
69 69
  typedef int VType;
70 70
  typedef concepts::Digraph Digraph;
71 71

	
72 72
  typedef Digraph::Node Node;
73 73
  typedef Digraph::Arc Arc;
74 74
  typedef concepts::ReadMap<Arc,VType> CapMap;
75 75
  typedef concepts::ReadWriteMap<Arc,VType> FlowMap;
76 76
  typedef concepts::WriteMap<Node,bool> CutMap;
77 77

	
78 78
  typedef Elevator<Digraph, Digraph::Node> Elev;
79 79
  typedef LinkedElevator<Digraph, Digraph::Node> LinkedElev;
80 80

	
81 81
  Digraph g;
82 82
  Node n;
83 83
  Arc e;
84 84
  CapMap cap;
85 85
  FlowMap flow;
86 86
  CutMap cut;
87 87
  VType v;
88 88
  bool b;
89 89

	
90 90
  typedef Preflow<Digraph, CapMap>
91 91
            ::SetFlowMap<FlowMap>
92 92
            ::SetElevator<Elev>
93 93
            ::SetStandardElevator<LinkedElev>
94 94
            ::Create PreflowType;
95 95
  PreflowType preflow_test(g, cap, n, n);
96 96
  const PreflowType& const_preflow_test = preflow_test;
97
  
98
  const PreflowType::Elevator& elev = const_preflow_test.elevator();
99
  preflow_test.elevator(const_cast<PreflowType::Elevator&>(elev));
100
  PreflowType::Tolerance tol = const_preflow_test.tolerance();
101
  preflow_test.tolerance(tol);
97 102

	
98 103
  preflow_test
99 104
    .capacityMap(cap)
100 105
    .flowMap(flow)
101 106
    .source(n)
102 107
    .target(n);
103 108

	
104 109
  preflow_test.init();
105 110
  preflow_test.init(cap);
106 111
  preflow_test.startFirstPhase();
107 112
  preflow_test.startSecondPhase();
108 113
  preflow_test.run();
109 114
  preflow_test.runMinCut();
110 115

	
111 116
  v = const_preflow_test.flowValue();
112 117
  v = const_preflow_test.flow(e);
113 118
  const FlowMap& fm = const_preflow_test.flowMap();
114 119
  b = const_preflow_test.minCut(n);
115 120
  const_preflow_test.minCutMap(cut);
116 121
  
117 122
  ignore_unused_variable_warning(fm);
118 123
}
119 124

	
120 125
int cutValue (const SmartDigraph& g,
121 126
              const SmartDigraph::NodeMap<bool>& cut,
122 127
              const SmartDigraph::ArcMap<int>& cap) {
123 128

	
124 129
  int c=0;
125 130
  for(SmartDigraph::ArcIt e(g); e!=INVALID; ++e) {
126 131
    if (cut[g.source(e)] && !cut[g.target(e)]) c+=cap[e];
127 132
  }
128 133
  return c;
129 134
}
130 135

	
131 136
bool checkFlow(const SmartDigraph& g,
132 137
               const SmartDigraph::ArcMap<int>& flow,
133 138
               const SmartDigraph::ArcMap<int>& cap,
134 139
               SmartDigraph::Node s, SmartDigraph::Node t) {
135 140

	
136 141
  for (SmartDigraph::ArcIt e(g); e != INVALID; ++e) {
137 142
    if (flow[e] < 0 || flow[e] > cap[e]) return false;
138 143
  }
139 144

	
140 145
  for (SmartDigraph::NodeIt n(g); n != INVALID; ++n) {
141 146
    if (n == s || n == t) continue;
142 147
    int sum = 0;
143 148
    for (SmartDigraph::OutArcIt e(g, n); e != INVALID; ++e) {
144 149
      sum += flow[e];
145 150
    }
146 151
    for (SmartDigraph::InArcIt e(g, n); e != INVALID; ++e) {
147 152
      sum -= flow[e];
148 153
    }
149 154
    if (sum != 0) return false;
150 155
  }
151 156
  return true;
152 157
}
153 158

	
154 159
int main() {
155 160

	
156 161
  typedef SmartDigraph Digraph;
157 162

	
158 163
  typedef Digraph::Node Node;
159 164
  typedef Digraph::NodeIt NodeIt;
160 165
  typedef Digraph::ArcIt ArcIt;
161 166
  typedef Digraph::ArcMap<int> CapMap;
162 167
  typedef Digraph::ArcMap<int> FlowMap;
163 168
  typedef Digraph::NodeMap<bool> CutMap;
164 169

	
165 170
  typedef Preflow<Digraph, CapMap> PType;
166 171

	
167 172
  Digraph g;
168 173
  Node s, t;
169 174
  CapMap cap(g);
170 175
  std::istringstream input(test_lgf);
171 176
  DigraphReader<Digraph>(g,input).
172 177
    arcMap("capacity", cap).
173 178
    node("source",s).
174 179
    node("target",t).
175 180
    run();
176 181

	
177 182
  PType preflow_test(g, cap, s, t);
178 183
  preflow_test.run();
179 184

	
180 185
  check(checkFlow(g, preflow_test.flowMap(), cap, s, t),
181 186
        "The flow is not feasible.");
182 187

	
183 188
  CutMap min_cut(g);
184 189
  preflow_test.minCutMap(min_cut);
185 190
  int min_cut_value=cutValue(g,min_cut,cap);
186 191

	
187 192
  check(preflow_test.flowValue() == min_cut_value,
188 193
        "The max flow value is not equal to the three min cut values.");
189 194

	
190 195
  FlowMap flow(g);
191 196
  for(ArcIt e(g); e!=INVALID; ++e) flow[e] = preflow_test.flowMap()[e];
192 197

	
193 198
  int flow_value=preflow_test.flowValue();
194 199

	
195 200
  for(ArcIt e(g); e!=INVALID; ++e) cap[e]=2*cap[e];
196 201
  preflow_test.init(flow);
197 202
  preflow_test.startFirstPhase();
198 203

	
199 204
  CutMap min_cut1(g);
200 205
  preflow_test.minCutMap(min_cut1);
201 206
  min_cut_value=cutValue(g,min_cut1,cap);
202 207

	
203 208
  check(preflow_test.flowValue() == min_cut_value &&
204 209
        min_cut_value == 2*flow_value,
205 210
        "The max flow value or the min cut value is wrong.");
206 211

	
207 212
  preflow_test.startSecondPhase();
208 213

	
209 214
  check(checkFlow(g, preflow_test.flowMap(), cap, s, t),
210 215
        "The flow is not feasible.");
211 216

	
212 217
  CutMap min_cut2(g);
213 218
  preflow_test.minCutMap(min_cut2);
214 219
  min_cut_value=cutValue(g,min_cut2,cap);
215 220

	
216 221
  check(preflow_test.flowValue() == min_cut_value &&
217 222
        min_cut_value == 2*flow_value,
218 223
        "The max flow value or the three min cut values were not doubled");
219 224

	
220 225

	
221 226
  preflow_test.flowMap(flow);
222 227

	
223 228
  NodeIt tmp1(g,s);
224 229
  ++tmp1;
225 230
  if ( tmp1 != INVALID ) s=tmp1;
226 231

	
227 232
  NodeIt tmp2(g,t);
228 233
  ++tmp2;
229 234
  if ( tmp2 != INVALID ) t=tmp2;
230 235

	
231 236
  preflow_test.source(s);
232 237
  preflow_test.target(t);
233 238

	
234 239
  preflow_test.run();
235 240

	
236 241
  CutMap min_cut3(g);
237 242
  preflow_test.minCutMap(min_cut3);
238 243
  min_cut_value=cutValue(g,min_cut3,cap);
239 244

	
240 245

	
241 246
  check(preflow_test.flowValue() == min_cut_value,
242 247
        "The max flow value or the three min cut values are incorrect.");
243 248

	
244 249
  return 0;
245 250
}
0 comments (0 inline)