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