gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Many doc improvements for Circulation (#175) - More precise doc for members. - Several doc fixes. - Add doc for public types. - Better formulations. - Add useful notes to the problem description. - Use supply instead of excess in the doc. - Hide the doc of the traits class parameter. - Use \tparam for template parameters.
0 1 0
default
1 file changed with 230 insertions and 131 deletions:
↑ Collapse diff ↑
Ignore white space 48 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2008
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_CIRCULATION_H
20 20
#define LEMON_CIRCULATION_H
21 21

	
22
#include <iostream>
23
#include <queue>
24 22
#include <lemon/tolerance.h>
25 23
#include <lemon/elevator.h>
26 24

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

	
33 31
  /// \brief Default traits class of Circulation class.
34 32
  ///
35 33
  /// Default traits class of Circulation class.
36
  /// \param _Graph Digraph type.
37
  /// \param _CapacityMap Type of capacity map.
38
  template <typename _Graph, typename _LCapMap,
34
  /// \tparam _Diraph Digraph type.
35
  /// \tparam _LCapMap Lower bound capacity map type.
36
  /// \tparam _UCapMap Upper bound capacity map type.
37
  /// \tparam _DeltaMap Delta map type.
38
  template <typename _Diraph, typename _LCapMap,
39 39
            typename _UCapMap, typename _DeltaMap>
40 40
  struct CirculationDefaultTraits {
41 41

	
42
    /// \brief The digraph type the algorithm runs on.
43
    typedef _Graph Digraph;
42
    /// \brief The type of the digraph the algorithm runs on.
43
    typedef _Diraph Digraph;
44 44

	
45 45
    /// \brief The type of the map that stores the circulation lower
46 46
    /// bound.
47 47
    ///
48 48
    /// The type of the map that stores the circulation lower bound.
49 49
    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
50 50
    typedef _LCapMap LCapMap;
51 51

	
52 52
    /// \brief The type of the map that stores the circulation upper
53 53
    /// bound.
54 54
    ///
55 55
    /// The type of the map that stores the circulation upper bound.
56 56
    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
57 57
    typedef _UCapMap UCapMap;
58 58

	
59
    /// \brief The type of the map that stores the upper bound of
60
    /// node excess.
59
    /// \brief The type of the map that stores the lower bound for
60
    /// the supply of the nodes.
61 61
    ///
62
    /// The type of the map that stores the lower bound of node
63
    /// excess. It must meet the \ref concepts::ReadMap "ReadMap"
62
    /// The type of the map that stores the lower bound for the supply
63
    /// of the nodes. It must meet the \ref concepts::ReadMap "ReadMap"
64 64
    /// concept.
65 65
    typedef _DeltaMap DeltaMap;
66 66

	
67
    /// \brief The type of the length of the arcs.
67
    /// \brief The type of the flow values.
68 68
    typedef typename DeltaMap::Value Value;
69 69

	
70
    /// \brief The map type that stores the flow values.
70
    /// \brief The type of the map that stores the flow values.
71 71
    ///
72
    /// The map type that stores the flow values.
72
    /// The type of the map that stores the flow values.
73 73
    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
74 74
    typedef typename Digraph::template ArcMap<Value> FlowMap;
75 75

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

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

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

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

	
108 108
  };
109 109

	
110
  ///Push-relabel algorithm for the Network Circulation Problem.
110
  /**
111
     \brief Push-relabel algorithm for the network circulation problem.
111 112

	
112
  /**
113 113
     \ingroup max_flow
114
     This class implements a push-relabel algorithm
115
     or the Network Circulation Problem.
114
     This class implements a push-relabel algorithm for the network
115
     circulation problem.
116
     It is to find a feasible circulation when lower and upper bounds
117
     are given for the flow values on the arcs and lower bounds
118
     are given for the supply values of the nodes.
119

	
116 120
     The exact formulation of this problem is the following.
117
     \f[\sum_{e\in\rho(v)}x(e)-\sum_{e\in\delta(v)}x(e)\leq
118
     -delta(v)\quad \forall v\in V \f]
119
     \f[ lo(e)\leq x(e) \leq up(e) \quad \forall e\in E \f]
121
     Let \f$G=(V,A)\f$ be a digraph,
122
     \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$,
123
     \f$delta: V\rightarrow\mathbf{R}\f$. Find a feasible circulation
124
     \f$f: A\rightarrow\mathbf{R}^+_0\f$ so that
125
     \f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a)
126
     \geq delta(v) \quad \forall v\in V, \f]
127
     \f[ lower(a)\leq f(a) \leq upper(a) \quad \forall a\in A. \f]
128
     \note \f$delta(v)\f$ specifies a lower bound for the supply of node
129
     \f$v\f$. It can be either positive or negative, however note that
130
     \f$\sum_{v\in V}delta(v)\f$ should be zero or negative in order to
131
     have a feasible solution.
132

	
133
     \note A special case of this problem is when
134
     \f$\sum_{v\in V}delta(v) = 0\f$. Then the supply of each node \f$v\f$
135
     will be \e equal \e to \f$delta(v)\f$, if a circulation can be found.
136
     Thus a feasible solution for the
137
     \ref min_cost_flow "minimum cost flow" problem can be calculated
138
     in this way.
139

	
140
     \tparam _Digraph The type of the digraph the algorithm runs on.
141
     \tparam _LCapMap The type of the lower bound capacity map. The default
142
     map type is \ref concepts::Digraph::ArcMap "_Digraph::ArcMap<int>".
143
     \tparam _UCapMap The type of the upper bound capacity map. The default
144
     map type is \c _LCapMap.
145
     \tparam _DeltaMap The type of the map that stores the lower bound
146
     for the supply of the nodes. The default map type is
147
     \c _Digraph::ArcMap<_UCapMap::Value>.
120 148
  */
121
  template<class _Graph,
122
           class _LCapMap=typename _Graph::template ArcMap<int>,
123
           class _UCapMap=_LCapMap,
124
           class _DeltaMap=typename _Graph::template NodeMap<
125
             typename _UCapMap::Value>,
126
           class _Traits=CirculationDefaultTraits<_Graph, _LCapMap,
127
                                                  _UCapMap, _DeltaMap> >
149
#ifdef DOXYGEN
150
template< typename _Digraph,
151
          typename _LCapMap,
152
          typename _UCapMap,
153
          typename _DeltaMap,
154
          typename _Traits >
155
#else
156
template< typename _Digraph,
157
          typename _LCapMap = typename _Digraph::template ArcMap<int>,
158
          typename _UCapMap = _LCapMap,
159
          typename _DeltaMap = typename _Digraph::
160
                               template NodeMap<typename _UCapMap::Value>,
161
          typename _Traits=CirculationDefaultTraits<_Digraph, _LCapMap,
162
                                                    _UCapMap, _DeltaMap> >
163
#endif
128 164
  class Circulation {
165
  public:
129 166

	
167
    ///The \ref CirculationDefaultTraits "traits class" of the algorithm.
130 168
    typedef _Traits Traits;
169
    ///The type of the digraph the algorithm runs on.
131 170
    typedef typename Traits::Digraph Digraph;
132
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
133

	
171
    ///The type of the flow values.
134 172
    typedef typename Traits::Value Value;
135 173

	
174
    /// The type of the lower bound capacity map.
136 175
    typedef typename Traits::LCapMap LCapMap;
176
    /// The type of the upper bound capacity map.
137 177
    typedef typename Traits::UCapMap UCapMap;
178
    /// \brief The type of the map that stores the lower bound for
179
    /// the supply of the nodes.
138 180
    typedef typename Traits::DeltaMap DeltaMap;
181
    ///The type of the flow map.
139 182
    typedef typename Traits::FlowMap FlowMap;
183

	
184
    ///The type of the elevator.
140 185
    typedef typename Traits::Elevator Elevator;
186
    ///The type of the tolerance.
141 187
    typedef typename Traits::Tolerance Tolerance;
142 188

	
143
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
189
  private:
190

	
191
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
144 192

	
145 193
    const Digraph &_g;
146 194
    int _node_num;
147 195

	
148 196
    const LCapMap *_lo;
149 197
    const UCapMap *_up;
150 198
    const DeltaMap *_delta;
151 199

	
152 200
    FlowMap *_flow;
153 201
    bool _local_flow;
154 202

	
155 203
    Elevator* _level;
156 204
    bool _local_level;
157 205

	
206
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
158 207
    ExcessMap* _excess;
159 208

	
160 209
    Tolerance _tol;
161 210
    int _el;
162 211

	
163 212
  public:
164 213

	
165 214
    typedef Circulation Create;
166 215

	
167
    ///\name Named template parameters
216
    ///\name Named Template Parameters
168 217

	
169 218
    ///@{
170 219

	
171 220
    template <typename _FlowMap>
172 221
    struct SetFlowMapTraits : public Traits {
173 222
      typedef _FlowMap FlowMap;
174 223
      static FlowMap *createFlowMap(const Digraph&) {
175 224
        LEMON_ASSERT(false, "FlowMap is not initialized");
176 225
        return 0; // ignore warnings
177 226
      }
178 227
    };
179 228

	
180 229
    /// \brief \ref named-templ-param "Named parameter" for setting
181 230
    /// FlowMap type
182 231
    ///
183 232
    /// \ref named-templ-param "Named parameter" for setting FlowMap
184
    /// type
233
    /// type.
185 234
    template <typename _FlowMap>
186 235
    struct SetFlowMap
187 236
      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
188 237
                           SetFlowMapTraits<_FlowMap> > {
189 238
      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
190 239
                          SetFlowMapTraits<_FlowMap> > Create;
191 240
    };
192 241

	
193 242
    template <typename _Elevator>
194 243
    struct SetElevatorTraits : public Traits {
195 244
      typedef _Elevator Elevator;
196 245
      static Elevator *createElevator(const Digraph&, int) {
197 246
        LEMON_ASSERT(false, "Elevator is not initialized");
198 247
        return 0; // ignore warnings
199 248
      }
200 249
    };
201 250

	
202 251
    /// \brief \ref named-templ-param "Named parameter" for setting
203 252
    /// Elevator type
204 253
    ///
205 254
    /// \ref named-templ-param "Named parameter" for setting Elevator
206
    /// type
255
    /// type. If this named parameter is used, then an external
256
    /// elevator object must be passed to the algorithm using the
257
    /// \ref elevator(Elevator&) "elevator()" function before calling
258
    /// \ref run() or \ref init().
259
    /// \sa SetStandardElevator
207 260
    template <typename _Elevator>
208 261
    struct SetElevator
209 262
      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
210 263
                           SetElevatorTraits<_Elevator> > {
211 264
      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
212 265
                          SetElevatorTraits<_Elevator> > Create;
213 266
    };
214 267

	
215 268
    template <typename _Elevator>
216 269
    struct SetStandardElevatorTraits : public Traits {
217 270
      typedef _Elevator Elevator;
218 271
      static Elevator *createElevator(const Digraph& digraph, int max_level) {
219 272
        return new Elevator(digraph, max_level);
220 273
      }
221 274
    };
222 275

	
223 276
    /// \brief \ref named-templ-param "Named parameter" for setting
224
    /// Elevator type
277
    /// Elevator type with automatic allocation
225 278
    ///
226 279
    /// \ref named-templ-param "Named parameter" for setting Elevator
227
    /// type. The Elevator should be standard constructor interface, ie.
228
    /// the digraph and the maximum level should be passed to it.
280
    /// type with automatic allocation.
281
    /// The Elevator should have standard constructor interface to be
282
    /// able to automatically created by the algorithm (i.e. the
283
    /// digraph and the maximum level should be passed to it).
284
    /// However an external elevator object could also be passed to the
285
    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
286
    /// before calling \ref run() or \ref init().
287
    /// \sa SetElevator
229 288
    template <typename _Elevator>
230 289
    struct SetStandardElevator
231 290
      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
232 291
                       SetStandardElevatorTraits<_Elevator> > {
233 292
      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
234 293
                      SetStandardElevatorTraits<_Elevator> > Create;
235 294
    };
236 295

	
237 296
    /// @}
238 297

	
239 298
  protected:
240 299

	
241 300
    Circulation() {}
242 301

	
243 302
  public:
244 303

	
245 304
    /// The constructor of the class.
246 305

	
247 306
    /// The constructor of the class.
248 307
    /// \param g The digraph the algorithm runs on.
249 308
    /// \param lo The lower bound capacity of the arcs.
250 309
    /// \param up The upper bound capacity of the arcs.
251
    /// \param delta The lower bound on node excess.
310
    /// \param delta The lower bound for the supply of the nodes.
252 311
    Circulation(const Digraph &g,const LCapMap &lo,
253 312
                const UCapMap &up,const DeltaMap &delta)
254 313
      : _g(g), _node_num(),
255 314
        _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
256 315
        _level(0), _local_level(false), _excess(0), _el() {}
257 316

	
258
    /// Destrcutor.
317
    /// Destructor.
259 318
    ~Circulation() {
260 319
      destroyStructures();
261 320
    }
262 321

	
322

	
263 323
  private:
264 324

	
265 325
    void createStructures() {
266 326
      _node_num = _el = countNodes(_g);
267 327

	
268 328
      if (!_flow) {
269 329
        _flow = Traits::createFlowMap(_g);
270 330
        _local_flow = true;
271 331
      }
272 332
      if (!_level) {
273 333
        _level = Traits::createElevator(_g, _node_num);
274 334
        _local_level = true;
275 335
      }
276 336
      if (!_excess) {
277 337
        _excess = new ExcessMap(_g);
278 338
      }
279 339
    }
280 340

	
281 341
    void destroyStructures() {
282 342
      if (_local_flow) {
283 343
        delete _flow;
284 344
      }
285 345
      if (_local_level) {
286 346
        delete _level;
287 347
      }
288 348
      if (_excess) {
289 349
        delete _excess;
290 350
      }
291 351
    }
292 352

	
293 353
  public:
294 354

	
295 355
    /// Sets the lower bound capacity map.
296 356

	
297 357
    /// Sets the lower bound capacity map.
298
    /// \return \c (*this)
358
    /// \return <tt>(*this)</tt>
299 359
    Circulation& lowerCapMap(const LCapMap& map) {
300 360
      _lo = &map;
301 361
      return *this;
302 362
    }
303 363

	
304 364
    /// Sets the upper bound capacity map.
305 365

	
306 366
    /// Sets the upper bound capacity map.
307
    /// \return \c (*this)
367
    /// \return <tt>(*this)</tt>
308 368
    Circulation& upperCapMap(const LCapMap& map) {
309 369
      _up = &map;
310 370
      return *this;
311 371
    }
312 372

	
313
    /// Sets the lower bound map on excess.
373
    /// Sets the lower bound map for the supply of the nodes.
314 374

	
315
    /// Sets the lower bound map on excess.
316
    /// \return \c (*this)
375
    /// Sets the lower bound map for the supply of the nodes.
376
    /// \return <tt>(*this)</tt>
317 377
    Circulation& deltaMap(const DeltaMap& map) {
318 378
      _delta = &map;
319 379
      return *this;
320 380
    }
321 381

	
382
    /// \brief Sets the flow map.
383
    ///
322 384
    /// Sets the flow map.
323

	
324
    /// Sets the flow map.
325
    /// \return \c (*this)
385
    /// If you don't use this function before calling \ref run() or
386
    /// \ref init(), an instance will be allocated automatically.
387
    /// The destructor deallocates this automatically allocated map,
388
    /// of course.
389
    /// \return <tt>(*this)</tt>
326 390
    Circulation& flowMap(FlowMap& map) {
327 391
      if (_local_flow) {
328 392
        delete _flow;
329 393
        _local_flow = false;
330 394
      }
331 395
      _flow = &map;
332 396
      return *this;
333 397
    }
334 398

	
335
    /// Returns the flow map.
336

	
337
    /// \return The flow map.
399
    /// \brief Sets the elevator used by algorithm.
338 400
    ///
339
    const FlowMap& flowMap() {
340
      return *_flow;
341
    }
342

	
343
    /// Sets the elevator.
344

	
345
    /// Sets the elevator.
346
    /// \return \c (*this)
401
    /// Sets the elevator used by algorithm.
402
    /// If you don't use this function before calling \ref run() or
403
    /// \ref init(), an instance will be allocated automatically.
404
    /// The destructor deallocates this automatically allocated elevator,
405
    /// of course.
406
    /// \return <tt>(*this)</tt>
347 407
    Circulation& elevator(Elevator& elevator) {
348 408
      if (_local_level) {
349 409
        delete _level;
350 410
        _local_level = false;
351 411
      }
352 412
      _level = &elevator;
353 413
      return *this;
354 414
    }
355 415

	
356
    /// Returns the elevator.
357

	
358
    /// \return The elevator.
416
    /// \brief Returns a const reference to the elevator.
359 417
    ///
418
    /// Returns a const reference to the elevator.
419
    ///
420
    /// \pre Either \ref run() or \ref init() must be called before
421
    /// using this function.
360 422
    const Elevator& elevator() {
361 423
      return *_level;
362 424
    }
363 425

	
426
    /// \brief Sets the tolerance used by algorithm.
427
    ///
364 428
    /// Sets the tolerance used by algorithm.
365

	
366
    /// Sets the tolerance used by algorithm.
367
    ///
368 429
    Circulation& tolerance(const Tolerance& tolerance) const {
369 430
      _tol = tolerance;
370 431
      return *this;
371 432
    }
372 433

	
373
    /// Returns the tolerance used by algorithm.
374

	
375
    /// Returns the tolerance used by algorithm.
434
    /// \brief Returns a const reference to the tolerance.
376 435
    ///
436
    /// Returns a const reference to the tolerance.
377 437
    const Tolerance& tolerance() const {
378 438
      return tolerance;
379 439
    }
380 440

	
381
    /// \name Execution control
382
    /// The simplest way to execute the algorithm is to use one of the
383
    /// member functions called \c run().
384
    /// \n
385
    /// If you need more control on initial solution or execution then
386
    /// you have to call one \ref init() function and then the start()
387
    /// function.
441
    /// \name Execution Control
442
    /// The simplest way to execute the algorithm is to call \ref run().\n
443
    /// If you need more control on the initial solution or the execution,
444
    /// first you have to call one of the \ref init() functions, then
445
    /// the \ref start() function.
388 446

	
389 447
    ///@{
390 448

	
391 449
    /// Initializes the internal data structures.
392 450

	
393
    /// Initializes the internal data structures. This function sets
394
    /// all flow values to the lower bound.
395
    /// \return This function returns false if the initialization
396
    /// process found a barrier.
451
    /// Initializes the internal data structures and sets all flow values
452
    /// to the lower bound.
397 453
    void init()
398 454
    {
399 455
      createStructures();
400 456

	
401 457
      for(NodeIt n(_g);n!=INVALID;++n) {
402 458
        _excess->set(n, (*_delta)[n]);
403 459
      }
404 460

	
405 461
      for (ArcIt e(_g);e!=INVALID;++e) {
406 462
        _flow->set(e, (*_lo)[e]);
407 463
        _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
408 464
        _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
409 465
      }
410 466

	
411 467
      // global relabeling tested, but in general case it provides
412 468
      // worse performance for random digraphs
413 469
      _level->initStart();
414 470
      for(NodeIt n(_g);n!=INVALID;++n)
415 471
        _level->initAddItem(n);
416 472
      _level->initFinish();
417 473
      for(NodeIt n(_g);n!=INVALID;++n)
418 474
        if(_tol.positive((*_excess)[n]))
419 475
          _level->activate(n);
420 476
    }
421 477

	
422
    /// Initializes the internal data structures.
478
    /// Initializes the internal data structures using a greedy approach.
423 479

	
424
    /// Initializes the internal data structures. This functions uses
425
    /// greedy approach to construct the initial solution.
480
    /// Initializes the internal data structures using a greedy approach
481
    /// to construct the initial solution.
426 482
    void greedyInit()
427 483
    {
428 484
      createStructures();
429 485

	
430 486
      for(NodeIt n(_g);n!=INVALID;++n) {
431 487
        _excess->set(n, (*_delta)[n]);
432 488
      }
433 489

	
434 490
      for (ArcIt e(_g);e!=INVALID;++e) {
435 491
        if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
436 492
          _flow->set(e, (*_up)[e]);
437 493
          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
438 494
          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
439 495
        } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
440 496
          _flow->set(e, (*_lo)[e]);
441 497
          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_lo)[e]);
442 498
          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_lo)[e]);
443 499
        } else {
444 500
          Value fc = -(*_excess)[_g.target(e)];
445 501
          _flow->set(e, fc);
446 502
          _excess->set(_g.target(e), 0);
447 503
          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
448 504
        }
449 505
      }
450 506

	
451 507
      _level->initStart();
452 508
      for(NodeIt n(_g);n!=INVALID;++n)
453 509
        _level->initAddItem(n);
454 510
      _level->initFinish();
455 511
      for(NodeIt n(_g);n!=INVALID;++n)
456 512
        if(_tol.positive((*_excess)[n]))
457 513
          _level->activate(n);
458 514
    }
459 515

	
460
    ///Starts the algorithm
516
    ///Executes the algorithm
461 517

	
462
    ///This function starts the algorithm.
463
    ///\return This function returns true if it found a feasible circulation.
518
    ///This function executes the algorithm.
519
    ///
520
    ///\return \c true if a feasible circulation is found.
464 521
    ///
465 522
    ///\sa barrier()
523
    ///\sa barrierMap()
466 524
    bool start()
467 525
    {
468 526

	
469 527
      Node act;
470 528
      Node bact=INVALID;
471 529
      Node last_activated=INVALID;
472 530
      while((act=_level->highestActive())!=INVALID) {
473 531
        int actlevel=(*_level)[act];
474 532
        int mlevel=_node_num;
475 533
        Value exc=(*_excess)[act];
476 534

	
477 535
        for(OutArcIt e(_g,act);e!=INVALID; ++e) {
478 536
          Node v = _g.target(e);
479 537
          Value fc=(*_up)[e]-(*_flow)[e];
480 538
          if(!_tol.positive(fc)) continue;
481 539
          if((*_level)[v]<actlevel) {
482 540
            if(!_tol.less(fc, exc)) {
483 541
              _flow->set(e, (*_flow)[e] + exc);
484 542
              _excess->set(v, (*_excess)[v] + exc);
485 543
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
486 544
                _level->activate(v);
487 545
              _excess->set(act,0);
488 546
              _level->deactivate(act);
489 547
              goto next_l;
... ...
@@ -522,135 +580,176 @@
522 580
          }
523 581
          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
524 582
        }
525 583

	
526 584
        _excess->set(act, exc);
527 585
        if(!_tol.positive(exc)) _level->deactivate(act);
528 586
        else if(mlevel==_node_num) {
529 587
          _level->liftHighestActiveToTop();
530 588
          _el = _node_num;
531 589
          return false;
532 590
        }
533 591
        else {
534 592
          _level->liftHighestActive(mlevel+1);
535 593
          if(_level->onLevel(actlevel)==0) {
536 594
            _el = actlevel;
537 595
            return false;
538 596
          }
539 597
        }
540 598
      next_l:
541 599
        ;
542 600
      }
543 601
      return true;
544 602
    }
545 603

	
546
    /// Runs the circulation algorithm.
604
    /// Runs the algorithm.
547 605

	
548
    /// Runs the circulation algorithm.
549
    /// \note fc.run() is just a shortcut of the following code.
606
    /// This function runs the algorithm.
607
    ///
608
    /// \return \c true if a feasible circulation is found.
609
    ///
610
    /// \note Apart from the return value, c.run() is just a shortcut of
611
    /// the following code.
550 612
    /// \code
551
    ///   fc.greedyInit();
552
    ///   return fc.start();
613
    ///   c.greedyInit();
614
    ///   c.start();
553 615
    /// \endcode
554 616
    bool run() {
555 617
      greedyInit();
556 618
      return start();
557 619
    }
558 620

	
559 621
    /// @}
560 622

	
561 623
    /// \name Query Functions
562
    /// The result of the %Circulation algorithm can be obtained using
563
    /// these functions.
564
    /// \n
565
    /// Before the use of these functions,
566
    /// either run() or start() must be called.
624
    /// The results of the circulation algorithm can be obtained using
625
    /// these functions.\n
626
    /// Either \ref run() or \ref start() should be called before
627
    /// using them.
567 628

	
568 629
    ///@{
569 630

	
631
    /// \brief Returns the flow on the given arc.
632
    ///
633
    /// Returns the flow on the given arc.
634
    ///
635
    /// \pre Either \ref run() or \ref init() must be called before
636
    /// using this function.
637
    Value flow(const Arc& arc) const {
638
      return (*_flow)[arc];
639
    }
640

	
641
    /// \brief Returns a const reference to the flow map.
642
    ///
643
    /// Returns a const reference to the arc map storing the found flow.
644
    ///
645
    /// \pre Either \ref run() or \ref init() must be called before
646
    /// using this function.
647
    const FlowMap& flowMap() {
648
      return *_flow;
649
    }
650

	
570 651
    /**
571
       \brief Returns a barrier
572
       
652
       \brief Returns \c true if the given node is in a barrier.
653

	
573 654
       Barrier is a set \e B of nodes for which
574
       \f[ \sum_{v\in B}-delta(v)<
575
       \sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
576
       holds. The existence of a set with this property prooves that a feasible
577
       flow cannot exists.
655

	
656
       \f[ \sum_{a\in\delta_{out}(B)} upper(a) -
657
           \sum_{a\in\delta_{in}(B)} lower(a) < \sum_{v\in B}delta(v) \f]
658

	
659
       holds. The existence of a set with this property prooves that a
660
       feasible circualtion cannot exist.
661

	
662
       This function returns \c true if the given node is in the found
663
       barrier. If a feasible circulation is found, the function
664
       gives back \c false for every node.
665

	
666
       \pre Either \ref run() or \ref init() must be called before
667
       using this function.
668

	
669
       \sa barrierMap()
578 670
       \sa checkBarrier()
579
       \sa run()
580 671
    */
581
    template<class GT>
582
    void barrierMap(GT &bar)
672
    bool barrier(const Node& node)
673
    {
674
      return (*_level)[node] >= _el;
675
    }
676

	
677
    /// \brief Gives back a barrier.
678
    ///
679
    /// This function sets \c bar to the characteristic vector of the
680
    /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
681
    /// node map with \c bool (or convertible) value type.
682
    ///
683
    /// If a feasible circulation is found, the function gives back an
684
    /// empty set, so \c bar[v] will be \c false for all nodes \c v.
685
    ///
686
    /// \note This function calls \ref barrier() for each node,
687
    /// so it runs in \f$O(n)\f$ time.
688
    ///
689
    /// \pre Either \ref run() or \ref init() must be called before
690
    /// using this function.
691
    ///
692
    /// \sa barrier()
693
    /// \sa checkBarrier()
694
    template<class BarrierMap>
695
    void barrierMap(BarrierMap &bar)
583 696
    {
584 697
      for(NodeIt n(_g);n!=INVALID;++n)
585 698
        bar.set(n, (*_level)[n] >= _el);
586 699
    }
587 700

	
588
    ///Returns true if the node is in the barrier
589

	
590
    ///Returns true if the node is in the barrier
591
    ///\sa barrierMap()
592
    bool barrier(const Node& node)
593
    {
594
      return (*_level)[node] >= _el;
595
    }
596

	
597
    /// \brief Returns the flow on the arc.
598
    ///
599
    /// Sets the \c flowMap to the flow on the arcs. This method can
600
    /// be called after the second phase of algorithm.
601
    Value flow(const Arc& arc) const {
602
      return (*_flow)[arc];
603
    }
604

	
605 701
    /// @}
606 702

	
607 703
    /// \name Checker Functions
608
    /// The feasibility  of the results can be checked using
609
    /// these functions.
610
    /// \n
611
    /// Before the use of these functions,
612
    /// either run() or start() must be called.
704
    /// The feasibility of the results can be checked using
705
    /// these functions.\n
706
    /// Either \ref run() or \ref start() should be called before
707
    /// using them.
613 708

	
614 709
    ///@{
615 710

	
616
    ///Check if the  \c flow is a feasible circulation
711
    ///Check if the found flow is a feasible circulation
712

	
713
    ///Check if the found flow is a feasible circulation,
714
    ///
617 715
    bool checkFlow() {
618 716
      for(ArcIt e(_g);e!=INVALID;++e)
619 717
        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
620 718
      for(NodeIt n(_g);n!=INVALID;++n)
621 719
        {
622 720
          Value dif=-(*_delta)[n];
623 721
          for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
624 722
          for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
625 723
          if(_tol.negative(dif)) return false;
626 724
        }
627 725
      return true;
628 726
    }
629 727

	
630 728
    ///Check whether or not the last execution provides a barrier
631 729

	
632
    ///Check whether or not the last execution provides a barrier
730
    ///Check whether or not the last execution provides a barrier.
633 731
    ///\sa barrier()
732
    ///\sa barrierMap()
634 733
    bool checkBarrier()
635 734
    {
636 735
      Value delta=0;
637 736
      for(NodeIt n(_g);n!=INVALID;++n)
638 737
        if(barrier(n))
639 738
          delta-=(*_delta)[n];
640 739
      for(ArcIt e(_g);e!=INVALID;++e)
641 740
        {
642 741
          Node s=_g.source(e);
643 742
          Node t=_g.target(e);
644 743
          if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
645 744
          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
646 745
        }
647 746
      return _tol.negative(delta);
648 747
    }
649 748

	
650 749
    /// @}
651 750

	
652 751
  };
653 752

	
654 753
}
655 754

	
656 755
#endif
0 comments (0 inline)