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