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