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