lemon/circulation.h
author Peter Kovacs <kpeter@inf.elte.hu>
Sun, 30 Nov 2008 14:51:05 +0100
changeset 417 235be9d4b6ab
parent 416 26fd85a3087e
child 437 6a2a33ad261b
permissions -rw-r--r--
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.
alpar@414
     1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
alpar@414
     2
 *
alpar@414
     3
 * This file is a part of LEMON, a generic C++ optimization library.
alpar@414
     4
 *
alpar@414
     5
 * Copyright (C) 2003-2008
alpar@414
     6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
alpar@414
     7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
alpar@414
     8
 *
alpar@414
     9
 * Permission to use, modify and distribute this software is granted
alpar@414
    10
 * provided that this copyright notice appears in all copies. For
alpar@414
    11
 * precise terms see the accompanying LICENSE file.
alpar@414
    12
 *
alpar@414
    13
 * This software is provided "AS IS" with no warranty of any kind,
alpar@414
    14
 * express or implied, and with no claim as to its suitability for any
alpar@414
    15
 * purpose.
alpar@414
    16
 *
alpar@414
    17
 */
alpar@414
    18
alpar@414
    19
#ifndef LEMON_CIRCULATION_H
alpar@414
    20
#define LEMON_CIRCULATION_H
alpar@414
    21
alpar@414
    22
#include <lemon/tolerance.h>
alpar@414
    23
#include <lemon/elevator.h>
alpar@414
    24
alpar@414
    25
///\ingroup max_flow
alpar@414
    26
///\file
kpeter@417
    27
///\brief Push-relabel algorithm for finding a feasible circulation.
alpar@414
    28
///
alpar@414
    29
namespace lemon {
alpar@414
    30
alpar@414
    31
  /// \brief Default traits class of Circulation class.
alpar@414
    32
  ///
alpar@414
    33
  /// Default traits class of Circulation class.
kpeter@417
    34
  /// \tparam _Diraph Digraph type.
kpeter@417
    35
  /// \tparam _LCapMap Lower bound capacity map type.
kpeter@417
    36
  /// \tparam _UCapMap Upper bound capacity map type.
kpeter@417
    37
  /// \tparam _DeltaMap Delta map type.
kpeter@417
    38
  template <typename _Diraph, typename _LCapMap,
alpar@414
    39
            typename _UCapMap, typename _DeltaMap>
alpar@414
    40
  struct CirculationDefaultTraits {
alpar@414
    41
kpeter@417
    42
    /// \brief The type of the digraph the algorithm runs on.
kpeter@417
    43
    typedef _Diraph Digraph;
alpar@414
    44
alpar@414
    45
    /// \brief The type of the map that stores the circulation lower
alpar@414
    46
    /// bound.
alpar@414
    47
    ///
alpar@414
    48
    /// The type of the map that stores the circulation lower bound.
alpar@414
    49
    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
alpar@414
    50
    typedef _LCapMap LCapMap;
alpar@414
    51
alpar@414
    52
    /// \brief The type of the map that stores the circulation upper
alpar@414
    53
    /// bound.
alpar@414
    54
    ///
alpar@414
    55
    /// The type of the map that stores the circulation upper bound.
alpar@414
    56
    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
alpar@414
    57
    typedef _UCapMap UCapMap;
alpar@414
    58
kpeter@417
    59
    /// \brief The type of the map that stores the lower bound for
kpeter@417
    60
    /// the supply of the nodes.
alpar@414
    61
    ///
kpeter@417
    62
    /// The type of the map that stores the lower bound for the supply
kpeter@417
    63
    /// of the nodes. It must meet the \ref concepts::ReadMap "ReadMap"
alpar@414
    64
    /// concept.
alpar@414
    65
    typedef _DeltaMap DeltaMap;
alpar@414
    66
kpeter@417
    67
    /// \brief The type of the flow values.
alpar@414
    68
    typedef typename DeltaMap::Value Value;
alpar@414
    69
kpeter@417
    70
    /// \brief The type of the map that stores the flow values.
alpar@414
    71
    ///
kpeter@417
    72
    /// The type of the map that stores the flow values.
alpar@414
    73
    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
alpar@414
    74
    typedef typename Digraph::template ArcMap<Value> FlowMap;
alpar@414
    75
alpar@414
    76
    /// \brief Instantiates a FlowMap.
alpar@414
    77
    ///
alpar@414
    78
    /// This function instantiates a \ref FlowMap.
alpar@414
    79
    /// \param digraph The digraph, to which we would like to define
alpar@414
    80
    /// the flow map.
alpar@414
    81
    static FlowMap* createFlowMap(const Digraph& digraph) {
alpar@414
    82
      return new FlowMap(digraph);
alpar@414
    83
    }
alpar@414
    84
kpeter@417
    85
    /// \brief The elevator type used by the algorithm.
alpar@414
    86
    ///
kpeter@417
    87
    /// The elevator type used by the algorithm.
alpar@414
    88
    ///
alpar@414
    89
    /// \sa Elevator
alpar@414
    90
    /// \sa LinkedElevator
alpar@414
    91
    typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
alpar@414
    92
alpar@414
    93
    /// \brief Instantiates an Elevator.
alpar@414
    94
    ///
kpeter@417
    95
    /// This function instantiates an \ref Elevator.
alpar@414
    96
    /// \param digraph The digraph, to which we would like to define
alpar@414
    97
    /// the elevator.
alpar@414
    98
    /// \param max_level The maximum level of the elevator.
alpar@414
    99
    static Elevator* createElevator(const Digraph& digraph, int max_level) {
alpar@414
   100
      return new Elevator(digraph, max_level);
alpar@414
   101
    }
alpar@414
   102
alpar@414
   103
    /// \brief The tolerance used by the algorithm
alpar@414
   104
    ///
alpar@414
   105
    /// The tolerance used by the algorithm to handle inexact computation.
alpar@414
   106
    typedef lemon::Tolerance<Value> Tolerance;
alpar@414
   107
alpar@414
   108
  };
alpar@414
   109
kpeter@417
   110
  /**
kpeter@417
   111
     \brief Push-relabel algorithm for the network circulation problem.
alpar@414
   112
alpar@414
   113
     \ingroup max_flow
kpeter@417
   114
     This class implements a push-relabel algorithm for the network
kpeter@417
   115
     circulation problem.
kpeter@417
   116
     It is to find a feasible circulation when lower and upper bounds
kpeter@417
   117
     are given for the flow values on the arcs and lower bounds
kpeter@417
   118
     are given for the supply values of the nodes.
kpeter@417
   119
alpar@414
   120
     The exact formulation of this problem is the following.
kpeter@417
   121
     Let \f$G=(V,A)\f$ be a digraph,
kpeter@417
   122
     \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$,
kpeter@417
   123
     \f$delta: V\rightarrow\mathbf{R}\f$. Find a feasible circulation
kpeter@417
   124
     \f$f: A\rightarrow\mathbf{R}^+_0\f$ so that
kpeter@417
   125
     \f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a)
kpeter@417
   126
     \geq delta(v) \quad \forall v\in V, \f]
kpeter@417
   127
     \f[ lower(a)\leq f(a) \leq upper(a) \quad \forall a\in A. \f]
kpeter@417
   128
     \note \f$delta(v)\f$ specifies a lower bound for the supply of node
kpeter@417
   129
     \f$v\f$. It can be either positive or negative, however note that
kpeter@417
   130
     \f$\sum_{v\in V}delta(v)\f$ should be zero or negative in order to
kpeter@417
   131
     have a feasible solution.
kpeter@417
   132
kpeter@417
   133
     \note A special case of this problem is when
kpeter@417
   134
     \f$\sum_{v\in V}delta(v) = 0\f$. Then the supply of each node \f$v\f$
kpeter@417
   135
     will be \e equal \e to \f$delta(v)\f$, if a circulation can be found.
kpeter@417
   136
     Thus a feasible solution for the
kpeter@417
   137
     \ref min_cost_flow "minimum cost flow" problem can be calculated
kpeter@417
   138
     in this way.
kpeter@417
   139
kpeter@417
   140
     \tparam _Digraph The type of the digraph the algorithm runs on.
kpeter@417
   141
     \tparam _LCapMap The type of the lower bound capacity map. The default
kpeter@417
   142
     map type is \ref concepts::Digraph::ArcMap "_Digraph::ArcMap<int>".
kpeter@417
   143
     \tparam _UCapMap The type of the upper bound capacity map. The default
kpeter@417
   144
     map type is \c _LCapMap.
kpeter@417
   145
     \tparam _DeltaMap The type of the map that stores the lower bound
kpeter@417
   146
     for the supply of the nodes. The default map type is
kpeter@417
   147
     \c _Digraph::ArcMap<_UCapMap::Value>.
alpar@414
   148
  */
kpeter@417
   149
#ifdef DOXYGEN
kpeter@417
   150
template< typename _Digraph,
kpeter@417
   151
          typename _LCapMap,
kpeter@417
   152
          typename _UCapMap,
kpeter@417
   153
          typename _DeltaMap,
kpeter@417
   154
          typename _Traits >
kpeter@417
   155
#else
kpeter@417
   156
template< typename _Digraph,
kpeter@417
   157
          typename _LCapMap = typename _Digraph::template ArcMap<int>,
kpeter@417
   158
          typename _UCapMap = _LCapMap,
kpeter@417
   159
          typename _DeltaMap = typename _Digraph::
kpeter@417
   160
                               template NodeMap<typename _UCapMap::Value>,
kpeter@417
   161
          typename _Traits=CirculationDefaultTraits<_Digraph, _LCapMap,
kpeter@417
   162
                                                    _UCapMap, _DeltaMap> >
kpeter@417
   163
#endif
alpar@414
   164
  class Circulation {
kpeter@417
   165
  public:
alpar@414
   166
kpeter@417
   167
    ///The \ref CirculationDefaultTraits "traits class" of the algorithm.
alpar@414
   168
    typedef _Traits Traits;
kpeter@417
   169
    ///The type of the digraph the algorithm runs on.
alpar@414
   170
    typedef typename Traits::Digraph Digraph;
kpeter@417
   171
    ///The type of the flow values.
alpar@414
   172
    typedef typename Traits::Value Value;
alpar@414
   173
kpeter@417
   174
    /// The type of the lower bound capacity map.
alpar@414
   175
    typedef typename Traits::LCapMap LCapMap;
kpeter@417
   176
    /// The type of the upper bound capacity map.
alpar@414
   177
    typedef typename Traits::UCapMap UCapMap;
kpeter@417
   178
    /// \brief The type of the map that stores the lower bound for
kpeter@417
   179
    /// the supply of the nodes.
alpar@414
   180
    typedef typename Traits::DeltaMap DeltaMap;
kpeter@417
   181
    ///The type of the flow map.
alpar@414
   182
    typedef typename Traits::FlowMap FlowMap;
kpeter@417
   183
kpeter@417
   184
    ///The type of the elevator.
alpar@414
   185
    typedef typename Traits::Elevator Elevator;
kpeter@417
   186
    ///The type of the tolerance.
alpar@414
   187
    typedef typename Traits::Tolerance Tolerance;
alpar@414
   188
kpeter@417
   189
  private:
kpeter@417
   190
kpeter@417
   191
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
alpar@414
   192
alpar@414
   193
    const Digraph &_g;
alpar@414
   194
    int _node_num;
alpar@414
   195
alpar@414
   196
    const LCapMap *_lo;
alpar@414
   197
    const UCapMap *_up;
alpar@414
   198
    const DeltaMap *_delta;
alpar@414
   199
alpar@414
   200
    FlowMap *_flow;
alpar@414
   201
    bool _local_flow;
alpar@414
   202
alpar@414
   203
    Elevator* _level;
alpar@414
   204
    bool _local_level;
alpar@414
   205
kpeter@417
   206
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
alpar@414
   207
    ExcessMap* _excess;
alpar@414
   208
alpar@414
   209
    Tolerance _tol;
alpar@414
   210
    int _el;
alpar@414
   211
alpar@414
   212
  public:
alpar@414
   213
alpar@414
   214
    typedef Circulation Create;
alpar@414
   215
kpeter@417
   216
    ///\name Named Template Parameters
alpar@414
   217
alpar@414
   218
    ///@{
alpar@414
   219
alpar@414
   220
    template <typename _FlowMap>
alpar@416
   221
    struct SetFlowMapTraits : public Traits {
alpar@414
   222
      typedef _FlowMap FlowMap;
alpar@414
   223
      static FlowMap *createFlowMap(const Digraph&) {
alpar@414
   224
        LEMON_ASSERT(false, "FlowMap is not initialized");
alpar@414
   225
        return 0; // ignore warnings
alpar@414
   226
      }
alpar@414
   227
    };
alpar@414
   228
alpar@414
   229
    /// \brief \ref named-templ-param "Named parameter" for setting
alpar@414
   230
    /// FlowMap type
alpar@414
   231
    ///
alpar@414
   232
    /// \ref named-templ-param "Named parameter" for setting FlowMap
kpeter@417
   233
    /// type.
alpar@414
   234
    template <typename _FlowMap>
alpar@416
   235
    struct SetFlowMap
alpar@414
   236
      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
alpar@416
   237
                           SetFlowMapTraits<_FlowMap> > {
alpar@414
   238
      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
alpar@416
   239
                          SetFlowMapTraits<_FlowMap> > Create;
alpar@414
   240
    };
alpar@414
   241
alpar@414
   242
    template <typename _Elevator>
alpar@416
   243
    struct SetElevatorTraits : public Traits {
alpar@414
   244
      typedef _Elevator Elevator;
alpar@414
   245
      static Elevator *createElevator(const Digraph&, int) {
alpar@414
   246
        LEMON_ASSERT(false, "Elevator is not initialized");
alpar@414
   247
        return 0; // ignore warnings
alpar@414
   248
      }
alpar@414
   249
    };
alpar@414
   250
alpar@414
   251
    /// \brief \ref named-templ-param "Named parameter" for setting
alpar@414
   252
    /// Elevator type
alpar@414
   253
    ///
alpar@414
   254
    /// \ref named-templ-param "Named parameter" for setting Elevator
kpeter@417
   255
    /// type. If this named parameter is used, then an external
kpeter@417
   256
    /// elevator object must be passed to the algorithm using the
kpeter@417
   257
    /// \ref elevator(Elevator&) "elevator()" function before calling
kpeter@417
   258
    /// \ref run() or \ref init().
kpeter@417
   259
    /// \sa SetStandardElevator
alpar@414
   260
    template <typename _Elevator>
alpar@416
   261
    struct SetElevator
alpar@414
   262
      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
alpar@416
   263
                           SetElevatorTraits<_Elevator> > {
alpar@414
   264
      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
alpar@416
   265
                          SetElevatorTraits<_Elevator> > Create;
alpar@414
   266
    };
alpar@414
   267
alpar@414
   268
    template <typename _Elevator>
alpar@416
   269
    struct SetStandardElevatorTraits : public Traits {
alpar@414
   270
      typedef _Elevator Elevator;
alpar@414
   271
      static Elevator *createElevator(const Digraph& digraph, int max_level) {
alpar@414
   272
        return new Elevator(digraph, max_level);
alpar@414
   273
      }
alpar@414
   274
    };
alpar@414
   275
alpar@414
   276
    /// \brief \ref named-templ-param "Named parameter" for setting
kpeter@417
   277
    /// Elevator type with automatic allocation
alpar@414
   278
    ///
alpar@414
   279
    /// \ref named-templ-param "Named parameter" for setting Elevator
kpeter@417
   280
    /// type with automatic allocation.
kpeter@417
   281
    /// The Elevator should have standard constructor interface to be
kpeter@417
   282
    /// able to automatically created by the algorithm (i.e. the
kpeter@417
   283
    /// digraph and the maximum level should be passed to it).
kpeter@417
   284
    /// However an external elevator object could also be passed to the
kpeter@417
   285
    /// algorithm with the \ref elevator(Elevator&) "elevator()" function
kpeter@417
   286
    /// before calling \ref run() or \ref init().
kpeter@417
   287
    /// \sa SetElevator
alpar@414
   288
    template <typename _Elevator>
alpar@416
   289
    struct SetStandardElevator
alpar@414
   290
      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
alpar@416
   291
                       SetStandardElevatorTraits<_Elevator> > {
alpar@414
   292
      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
alpar@416
   293
                      SetStandardElevatorTraits<_Elevator> > Create;
alpar@414
   294
    };
alpar@414
   295
alpar@414
   296
    /// @}
alpar@414
   297
alpar@414
   298
  protected:
alpar@414
   299
alpar@414
   300
    Circulation() {}
alpar@414
   301
alpar@414
   302
  public:
alpar@414
   303
alpar@414
   304
    /// The constructor of the class.
alpar@414
   305
alpar@414
   306
    /// The constructor of the class.
alpar@414
   307
    /// \param g The digraph the algorithm runs on.
alpar@414
   308
    /// \param lo The lower bound capacity of the arcs.
alpar@414
   309
    /// \param up The upper bound capacity of the arcs.
kpeter@417
   310
    /// \param delta The lower bound for the supply of the nodes.
alpar@414
   311
    Circulation(const Digraph &g,const LCapMap &lo,
alpar@414
   312
                const UCapMap &up,const DeltaMap &delta)
alpar@414
   313
      : _g(g), _node_num(),
alpar@414
   314
        _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
alpar@414
   315
        _level(0), _local_level(false), _excess(0), _el() {}
alpar@414
   316
kpeter@417
   317
    /// Destructor.
alpar@414
   318
    ~Circulation() {
alpar@414
   319
      destroyStructures();
alpar@414
   320
    }
alpar@414
   321
kpeter@417
   322
alpar@414
   323
  private:
alpar@414
   324
alpar@414
   325
    void createStructures() {
alpar@414
   326
      _node_num = _el = countNodes(_g);
alpar@414
   327
alpar@414
   328
      if (!_flow) {
alpar@414
   329
        _flow = Traits::createFlowMap(_g);
alpar@414
   330
        _local_flow = true;
alpar@414
   331
      }
alpar@414
   332
      if (!_level) {
alpar@414
   333
        _level = Traits::createElevator(_g, _node_num);
alpar@414
   334
        _local_level = true;
alpar@414
   335
      }
alpar@414
   336
      if (!_excess) {
alpar@414
   337
        _excess = new ExcessMap(_g);
alpar@414
   338
      }
alpar@414
   339
    }
alpar@414
   340
alpar@414
   341
    void destroyStructures() {
alpar@414
   342
      if (_local_flow) {
alpar@414
   343
        delete _flow;
alpar@414
   344
      }
alpar@414
   345
      if (_local_level) {
alpar@414
   346
        delete _level;
alpar@414
   347
      }
alpar@414
   348
      if (_excess) {
alpar@414
   349
        delete _excess;
alpar@414
   350
      }
alpar@414
   351
    }
alpar@414
   352
alpar@414
   353
  public:
alpar@414
   354
alpar@414
   355
    /// Sets the lower bound capacity map.
alpar@414
   356
alpar@414
   357
    /// Sets the lower bound capacity map.
kpeter@417
   358
    /// \return <tt>(*this)</tt>
alpar@414
   359
    Circulation& lowerCapMap(const LCapMap& map) {
alpar@414
   360
      _lo = &map;
alpar@414
   361
      return *this;
alpar@414
   362
    }
alpar@414
   363
alpar@414
   364
    /// Sets the upper bound capacity map.
alpar@414
   365
alpar@414
   366
    /// Sets the upper bound capacity map.
kpeter@417
   367
    /// \return <tt>(*this)</tt>
alpar@414
   368
    Circulation& upperCapMap(const LCapMap& map) {
alpar@414
   369
      _up = &map;
alpar@414
   370
      return *this;
alpar@414
   371
    }
alpar@414
   372
kpeter@417
   373
    /// Sets the lower bound map for the supply of the nodes.
alpar@414
   374
kpeter@417
   375
    /// Sets the lower bound map for the supply of the nodes.
kpeter@417
   376
    /// \return <tt>(*this)</tt>
alpar@414
   377
    Circulation& deltaMap(const DeltaMap& map) {
alpar@414
   378
      _delta = &map;
alpar@414
   379
      return *this;
alpar@414
   380
    }
alpar@414
   381
kpeter@417
   382
    /// \brief Sets the flow map.
kpeter@417
   383
    ///
alpar@414
   384
    /// Sets the flow map.
kpeter@417
   385
    /// If you don't use this function before calling \ref run() or
kpeter@417
   386
    /// \ref init(), an instance will be allocated automatically.
kpeter@417
   387
    /// The destructor deallocates this automatically allocated map,
kpeter@417
   388
    /// of course.
kpeter@417
   389
    /// \return <tt>(*this)</tt>
alpar@414
   390
    Circulation& flowMap(FlowMap& map) {
alpar@414
   391
      if (_local_flow) {
alpar@414
   392
        delete _flow;
alpar@414
   393
        _local_flow = false;
alpar@414
   394
      }
alpar@414
   395
      _flow = &map;
alpar@414
   396
      return *this;
alpar@414
   397
    }
alpar@414
   398
kpeter@417
   399
    /// \brief Sets the elevator used by algorithm.
alpar@414
   400
    ///
kpeter@417
   401
    /// Sets the elevator used by algorithm.
kpeter@417
   402
    /// If you don't use this function before calling \ref run() or
kpeter@417
   403
    /// \ref init(), an instance will be allocated automatically.
kpeter@417
   404
    /// The destructor deallocates this automatically allocated elevator,
kpeter@417
   405
    /// of course.
kpeter@417
   406
    /// \return <tt>(*this)</tt>
alpar@414
   407
    Circulation& elevator(Elevator& elevator) {
alpar@414
   408
      if (_local_level) {
alpar@414
   409
        delete _level;
alpar@414
   410
        _local_level = false;
alpar@414
   411
      }
alpar@414
   412
      _level = &elevator;
alpar@414
   413
      return *this;
alpar@414
   414
    }
alpar@414
   415
kpeter@417
   416
    /// \brief Returns a const reference to the elevator.
alpar@414
   417
    ///
kpeter@417
   418
    /// Returns a const reference to the elevator.
kpeter@417
   419
    ///
kpeter@417
   420
    /// \pre Either \ref run() or \ref init() must be called before
kpeter@417
   421
    /// using this function.
alpar@414
   422
    const Elevator& elevator() {
alpar@414
   423
      return *_level;
alpar@414
   424
    }
alpar@414
   425
kpeter@417
   426
    /// \brief Sets the tolerance used by algorithm.
kpeter@417
   427
    ///
alpar@414
   428
    /// Sets the tolerance used by algorithm.
alpar@414
   429
    Circulation& tolerance(const Tolerance& tolerance) const {
alpar@414
   430
      _tol = tolerance;
alpar@414
   431
      return *this;
alpar@414
   432
    }
alpar@414
   433
kpeter@417
   434
    /// \brief Returns a const reference to the tolerance.
alpar@414
   435
    ///
kpeter@417
   436
    /// Returns a const reference to the tolerance.
alpar@414
   437
    const Tolerance& tolerance() const {
alpar@414
   438
      return tolerance;
alpar@414
   439
    }
alpar@414
   440
kpeter@417
   441
    /// \name Execution Control
kpeter@417
   442
    /// The simplest way to execute the algorithm is to call \ref run().\n
kpeter@417
   443
    /// If you need more control on the initial solution or the execution,
kpeter@417
   444
    /// first you have to call one of the \ref init() functions, then
kpeter@417
   445
    /// the \ref start() function.
alpar@414
   446
alpar@414
   447
    ///@{
alpar@414
   448
alpar@414
   449
    /// Initializes the internal data structures.
alpar@414
   450
kpeter@417
   451
    /// Initializes the internal data structures and sets all flow values
kpeter@417
   452
    /// to the lower bound.
alpar@414
   453
    void init()
alpar@414
   454
    {
alpar@414
   455
      createStructures();
alpar@414
   456
alpar@414
   457
      for(NodeIt n(_g);n!=INVALID;++n) {
alpar@414
   458
        _excess->set(n, (*_delta)[n]);
alpar@414
   459
      }
alpar@414
   460
alpar@414
   461
      for (ArcIt e(_g);e!=INVALID;++e) {
alpar@414
   462
        _flow->set(e, (*_lo)[e]);
alpar@414
   463
        _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
alpar@414
   464
        _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
alpar@414
   465
      }
alpar@414
   466
alpar@414
   467
      // global relabeling tested, but in general case it provides
alpar@414
   468
      // worse performance for random digraphs
alpar@414
   469
      _level->initStart();
alpar@414
   470
      for(NodeIt n(_g);n!=INVALID;++n)
alpar@414
   471
        _level->initAddItem(n);
alpar@414
   472
      _level->initFinish();
alpar@414
   473
      for(NodeIt n(_g);n!=INVALID;++n)
alpar@414
   474
        if(_tol.positive((*_excess)[n]))
alpar@414
   475
          _level->activate(n);
alpar@414
   476
    }
alpar@414
   477
kpeter@417
   478
    /// Initializes the internal data structures using a greedy approach.
alpar@414
   479
kpeter@417
   480
    /// Initializes the internal data structures using a greedy approach
kpeter@417
   481
    /// to construct the initial solution.
alpar@414
   482
    void greedyInit()
alpar@414
   483
    {
alpar@414
   484
      createStructures();
alpar@414
   485
alpar@414
   486
      for(NodeIt n(_g);n!=INVALID;++n) {
alpar@414
   487
        _excess->set(n, (*_delta)[n]);
alpar@414
   488
      }
alpar@414
   489
alpar@414
   490
      for (ArcIt e(_g);e!=INVALID;++e) {
alpar@414
   491
        if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
alpar@414
   492
          _flow->set(e, (*_up)[e]);
alpar@414
   493
          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
alpar@414
   494
          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
alpar@414
   495
        } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
alpar@414
   496
          _flow->set(e, (*_lo)[e]);
alpar@414
   497
          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_lo)[e]);
alpar@414
   498
          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_lo)[e]);
alpar@414
   499
        } else {
alpar@414
   500
          Value fc = -(*_excess)[_g.target(e)];
alpar@414
   501
          _flow->set(e, fc);
alpar@414
   502
          _excess->set(_g.target(e), 0);
alpar@414
   503
          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
alpar@414
   504
        }
alpar@414
   505
      }
alpar@414
   506
alpar@414
   507
      _level->initStart();
alpar@414
   508
      for(NodeIt n(_g);n!=INVALID;++n)
alpar@414
   509
        _level->initAddItem(n);
alpar@414
   510
      _level->initFinish();
alpar@414
   511
      for(NodeIt n(_g);n!=INVALID;++n)
alpar@414
   512
        if(_tol.positive((*_excess)[n]))
alpar@414
   513
          _level->activate(n);
alpar@414
   514
    }
alpar@414
   515
kpeter@417
   516
    ///Executes the algorithm
alpar@414
   517
kpeter@417
   518
    ///This function executes the algorithm.
kpeter@417
   519
    ///
kpeter@417
   520
    ///\return \c true if a feasible circulation is found.
alpar@414
   521
    ///
alpar@414
   522
    ///\sa barrier()
kpeter@417
   523
    ///\sa barrierMap()
alpar@414
   524
    bool start()
alpar@414
   525
    {
alpar@414
   526
alpar@414
   527
      Node act;
alpar@414
   528
      Node bact=INVALID;
alpar@414
   529
      Node last_activated=INVALID;
alpar@414
   530
      while((act=_level->highestActive())!=INVALID) {
alpar@414
   531
        int actlevel=(*_level)[act];
alpar@414
   532
        int mlevel=_node_num;
alpar@414
   533
        Value exc=(*_excess)[act];
alpar@414
   534
alpar@414
   535
        for(OutArcIt e(_g,act);e!=INVALID; ++e) {
alpar@414
   536
          Node v = _g.target(e);
alpar@414
   537
          Value fc=(*_up)[e]-(*_flow)[e];
alpar@414
   538
          if(!_tol.positive(fc)) continue;
alpar@414
   539
          if((*_level)[v]<actlevel) {
alpar@414
   540
            if(!_tol.less(fc, exc)) {
alpar@414
   541
              _flow->set(e, (*_flow)[e] + exc);
alpar@414
   542
              _excess->set(v, (*_excess)[v] + exc);
alpar@414
   543
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
alpar@414
   544
                _level->activate(v);
alpar@414
   545
              _excess->set(act,0);
alpar@414
   546
              _level->deactivate(act);
alpar@414
   547
              goto next_l;
alpar@414
   548
            }
alpar@414
   549
            else {
alpar@414
   550
              _flow->set(e, (*_up)[e]);
alpar@414
   551
              _excess->set(v, (*_excess)[v] + fc);
alpar@414
   552
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
alpar@414
   553
                _level->activate(v);
alpar@414
   554
              exc-=fc;
alpar@414
   555
            }
alpar@414
   556
          }
alpar@414
   557
          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
alpar@414
   558
        }
alpar@414
   559
        for(InArcIt e(_g,act);e!=INVALID; ++e) {
alpar@414
   560
          Node v = _g.source(e);
alpar@414
   561
          Value fc=(*_flow)[e]-(*_lo)[e];
alpar@414
   562
          if(!_tol.positive(fc)) continue;
alpar@414
   563
          if((*_level)[v]<actlevel) {
alpar@414
   564
            if(!_tol.less(fc, exc)) {
alpar@414
   565
              _flow->set(e, (*_flow)[e] - exc);
alpar@414
   566
              _excess->set(v, (*_excess)[v] + exc);
alpar@414
   567
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
alpar@414
   568
                _level->activate(v);
alpar@414
   569
              _excess->set(act,0);
alpar@414
   570
              _level->deactivate(act);
alpar@414
   571
              goto next_l;
alpar@414
   572
            }
alpar@414
   573
            else {
alpar@414
   574
              _flow->set(e, (*_lo)[e]);
alpar@414
   575
              _excess->set(v, (*_excess)[v] + fc);
alpar@414
   576
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
alpar@414
   577
                _level->activate(v);
alpar@414
   578
              exc-=fc;
alpar@414
   579
            }
alpar@414
   580
          }
alpar@414
   581
          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
alpar@414
   582
        }
alpar@414
   583
alpar@414
   584
        _excess->set(act, exc);
alpar@414
   585
        if(!_tol.positive(exc)) _level->deactivate(act);
alpar@414
   586
        else if(mlevel==_node_num) {
alpar@414
   587
          _level->liftHighestActiveToTop();
alpar@414
   588
          _el = _node_num;
alpar@414
   589
          return false;
alpar@414
   590
        }
alpar@414
   591
        else {
alpar@414
   592
          _level->liftHighestActive(mlevel+1);
alpar@414
   593
          if(_level->onLevel(actlevel)==0) {
alpar@414
   594
            _el = actlevel;
alpar@414
   595
            return false;
alpar@414
   596
          }
alpar@414
   597
        }
alpar@414
   598
      next_l:
alpar@414
   599
        ;
alpar@414
   600
      }
alpar@414
   601
      return true;
alpar@414
   602
    }
alpar@414
   603
kpeter@417
   604
    /// Runs the algorithm.
alpar@414
   605
kpeter@417
   606
    /// This function runs the algorithm.
kpeter@417
   607
    ///
kpeter@417
   608
    /// \return \c true if a feasible circulation is found.
kpeter@417
   609
    ///
kpeter@417
   610
    /// \note Apart from the return value, c.run() is just a shortcut of
kpeter@417
   611
    /// the following code.
alpar@414
   612
    /// \code
kpeter@417
   613
    ///   c.greedyInit();
kpeter@417
   614
    ///   c.start();
alpar@414
   615
    /// \endcode
alpar@414
   616
    bool run() {
alpar@414
   617
      greedyInit();
alpar@414
   618
      return start();
alpar@414
   619
    }
alpar@414
   620
alpar@414
   621
    /// @}
alpar@414
   622
alpar@414
   623
    /// \name Query Functions
kpeter@417
   624
    /// The results of the circulation algorithm can be obtained using
kpeter@417
   625
    /// these functions.\n
kpeter@417
   626
    /// Either \ref run() or \ref start() should be called before
kpeter@417
   627
    /// using them.
alpar@414
   628
alpar@414
   629
    ///@{
alpar@414
   630
kpeter@417
   631
    /// \brief Returns the flow on the given arc.
kpeter@417
   632
    ///
kpeter@417
   633
    /// Returns the flow on the given arc.
kpeter@417
   634
    ///
kpeter@417
   635
    /// \pre Either \ref run() or \ref init() must be called before
kpeter@417
   636
    /// using this function.
kpeter@417
   637
    Value flow(const Arc& arc) const {
kpeter@417
   638
      return (*_flow)[arc];
kpeter@417
   639
    }
kpeter@417
   640
kpeter@417
   641
    /// \brief Returns a const reference to the flow map.
kpeter@417
   642
    ///
kpeter@417
   643
    /// Returns a const reference to the arc map storing the found flow.
kpeter@417
   644
    ///
kpeter@417
   645
    /// \pre Either \ref run() or \ref init() must be called before
kpeter@417
   646
    /// using this function.
kpeter@417
   647
    const FlowMap& flowMap() {
kpeter@417
   648
      return *_flow;
kpeter@417
   649
    }
kpeter@417
   650
alpar@414
   651
    /**
kpeter@417
   652
       \brief Returns \c true if the given node is in a barrier.
kpeter@417
   653
alpar@414
   654
       Barrier is a set \e B of nodes for which
kpeter@417
   655
kpeter@417
   656
       \f[ \sum_{a\in\delta_{out}(B)} upper(a) -
kpeter@417
   657
           \sum_{a\in\delta_{in}(B)} lower(a) < \sum_{v\in B}delta(v) \f]
kpeter@417
   658
kpeter@417
   659
       holds. The existence of a set with this property prooves that a
kpeter@417
   660
       feasible circualtion cannot exist.
kpeter@417
   661
kpeter@417
   662
       This function returns \c true if the given node is in the found
kpeter@417
   663
       barrier. If a feasible circulation is found, the function
kpeter@417
   664
       gives back \c false for every node.
kpeter@417
   665
kpeter@417
   666
       \pre Either \ref run() or \ref init() must be called before
kpeter@417
   667
       using this function.
kpeter@417
   668
kpeter@417
   669
       \sa barrierMap()
alpar@414
   670
       \sa checkBarrier()
alpar@414
   671
    */
kpeter@417
   672
    bool barrier(const Node& node)
kpeter@417
   673
    {
kpeter@417
   674
      return (*_level)[node] >= _el;
kpeter@417
   675
    }
kpeter@417
   676
kpeter@417
   677
    /// \brief Gives back a barrier.
kpeter@417
   678
    ///
kpeter@417
   679
    /// This function sets \c bar to the characteristic vector of the
kpeter@417
   680
    /// found barrier. \c bar should be a \ref concepts::WriteMap "writable"
kpeter@417
   681
    /// node map with \c bool (or convertible) value type.
kpeter@417
   682
    ///
kpeter@417
   683
    /// If a feasible circulation is found, the function gives back an
kpeter@417
   684
    /// empty set, so \c bar[v] will be \c false for all nodes \c v.
kpeter@417
   685
    ///
kpeter@417
   686
    /// \note This function calls \ref barrier() for each node,
kpeter@417
   687
    /// so it runs in \f$O(n)\f$ time.
kpeter@417
   688
    ///
kpeter@417
   689
    /// \pre Either \ref run() or \ref init() must be called before
kpeter@417
   690
    /// using this function.
kpeter@417
   691
    ///
kpeter@417
   692
    /// \sa barrier()
kpeter@417
   693
    /// \sa checkBarrier()
kpeter@417
   694
    template<class BarrierMap>
kpeter@417
   695
    void barrierMap(BarrierMap &bar)
alpar@414
   696
    {
alpar@414
   697
      for(NodeIt n(_g);n!=INVALID;++n)
alpar@414
   698
        bar.set(n, (*_level)[n] >= _el);
alpar@414
   699
    }
alpar@414
   700
alpar@414
   701
    /// @}
alpar@414
   702
alpar@414
   703
    /// \name Checker Functions
kpeter@417
   704
    /// The feasibility of the results can be checked using
kpeter@417
   705
    /// these functions.\n
kpeter@417
   706
    /// Either \ref run() or \ref start() should be called before
kpeter@417
   707
    /// using them.
alpar@414
   708
alpar@414
   709
    ///@{
alpar@414
   710
kpeter@417
   711
    ///Check if the found flow is a feasible circulation
kpeter@417
   712
kpeter@417
   713
    ///Check if the found flow is a feasible circulation,
kpeter@417
   714
    ///
alpar@414
   715
    bool checkFlow() {
alpar@414
   716
      for(ArcIt e(_g);e!=INVALID;++e)
alpar@414
   717
        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
alpar@414
   718
      for(NodeIt n(_g);n!=INVALID;++n)
alpar@414
   719
        {
alpar@414
   720
          Value dif=-(*_delta)[n];
alpar@414
   721
          for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
alpar@414
   722
          for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
alpar@414
   723
          if(_tol.negative(dif)) return false;
alpar@414
   724
        }
alpar@414
   725
      return true;
alpar@414
   726
    }
alpar@414
   727
alpar@414
   728
    ///Check whether or not the last execution provides a barrier
alpar@414
   729
kpeter@417
   730
    ///Check whether or not the last execution provides a barrier.
alpar@414
   731
    ///\sa barrier()
kpeter@417
   732
    ///\sa barrierMap()
alpar@414
   733
    bool checkBarrier()
alpar@414
   734
    {
alpar@414
   735
      Value delta=0;
alpar@414
   736
      for(NodeIt n(_g);n!=INVALID;++n)
alpar@414
   737
        if(barrier(n))
alpar@414
   738
          delta-=(*_delta)[n];
alpar@414
   739
      for(ArcIt e(_g);e!=INVALID;++e)
alpar@414
   740
        {
alpar@414
   741
          Node s=_g.source(e);
alpar@414
   742
          Node t=_g.target(e);
alpar@414
   743
          if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
alpar@414
   744
          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
alpar@414
   745
        }
alpar@414
   746
      return _tol.negative(delta);
alpar@414
   747
    }
alpar@414
   748
alpar@414
   749
    /// @}
alpar@414
   750
alpar@414
   751
  };
alpar@414
   752
alpar@414
   753
}
alpar@414
   754
alpar@414
   755
#endif