lemon/circulation.h
author Alpar Juttner <alpar@cs.elte.hu>
Fri, 21 Nov 2008 14:42:47 +0000
changeset 414 5a7dbeaed70e
child 416 26fd85a3087e
permissions -rw-r--r--
Port Circulation from svn -r3516 (#175)
Namely,
- port the files
- apply the migrate script
- apply the unify script
- fix the compilation
- strip the demo input file
- break long lines
alpar@414
     1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
alpar@414
     2
 *
alpar@414
     3
 * This file is a part of LEMON, a generic C++ optimization library.
alpar@414
     4
 *
alpar@414
     5
 * Copyright (C) 2003-2008
alpar@414
     6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
alpar@414
     7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
alpar@414
     8
 *
alpar@414
     9
 * Permission to use, modify and distribute this software is granted
alpar@414
    10
 * provided that this copyright notice appears in all copies. For
alpar@414
    11
 * precise terms see the accompanying LICENSE file.
alpar@414
    12
 *
alpar@414
    13
 * This software is provided "AS IS" with no warranty of any kind,
alpar@414
    14
 * express or implied, and with no claim as to its suitability for any
alpar@414
    15
 * purpose.
alpar@414
    16
 *
alpar@414
    17
 */
alpar@414
    18
alpar@414
    19
#ifndef LEMON_CIRCULATION_H
alpar@414
    20
#define LEMON_CIRCULATION_H
alpar@414
    21
alpar@414
    22
#include <iostream>
alpar@414
    23
#include <queue>
alpar@414
    24
#include <lemon/tolerance.h>
alpar@414
    25
#include <lemon/elevator.h>
alpar@414
    26
alpar@414
    27
///\ingroup max_flow
alpar@414
    28
///\file
alpar@414
    29
///\brief Push-prelabel algorithm for finding a feasible circulation.
alpar@414
    30
///
alpar@414
    31
namespace lemon {
alpar@414
    32
alpar@414
    33
  /// \brief Default traits class of Circulation class.
alpar@414
    34
  ///
alpar@414
    35
  /// Default traits class of Circulation class.
alpar@414
    36
  /// \param _Graph Digraph type.
alpar@414
    37
  /// \param _CapacityMap Type of capacity map.
alpar@414
    38
  template <typename _Graph, typename _LCapMap,
alpar@414
    39
            typename _UCapMap, typename _DeltaMap>
alpar@414
    40
  struct CirculationDefaultTraits {
alpar@414
    41
alpar@414
    42
    /// \brief The digraph type the algorithm runs on.
alpar@414
    43
    typedef _Graph Digraph;
alpar@414
    44
alpar@414
    45
    /// \brief The type of the map that stores the circulation lower
alpar@414
    46
    /// bound.
alpar@414
    47
    ///
alpar@414
    48
    /// The type of the map that stores the circulation lower bound.
alpar@414
    49
    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
alpar@414
    50
    typedef _LCapMap LCapMap;
alpar@414
    51
alpar@414
    52
    /// \brief The type of the map that stores the circulation upper
alpar@414
    53
    /// bound.
alpar@414
    54
    ///
alpar@414
    55
    /// The type of the map that stores the circulation upper bound.
alpar@414
    56
    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
alpar@414
    57
    typedef _UCapMap UCapMap;
alpar@414
    58
alpar@414
    59
    /// \brief The type of the map that stores the upper bound of
alpar@414
    60
    /// node excess.
alpar@414
    61
    ///
alpar@414
    62
    /// The type of the map that stores the lower bound of node
alpar@414
    63
    /// excess. It must meet the \ref concepts::ReadMap "ReadMap"
alpar@414
    64
    /// concept.
alpar@414
    65
    typedef _DeltaMap DeltaMap;
alpar@414
    66
alpar@414
    67
    /// \brief The type of the length of the arcs.
alpar@414
    68
    typedef typename DeltaMap::Value Value;
alpar@414
    69
alpar@414
    70
    /// \brief The map type that stores the flow values.
alpar@414
    71
    ///
alpar@414
    72
    /// The map type that stores the flow values.
alpar@414
    73
    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
alpar@414
    74
    typedef typename Digraph::template ArcMap<Value> FlowMap;
alpar@414
    75
alpar@414
    76
    /// \brief Instantiates a FlowMap.
alpar@414
    77
    ///
alpar@414
    78
    /// This function instantiates a \ref FlowMap.
alpar@414
    79
    /// \param digraph The digraph, to which we would like to define
alpar@414
    80
    /// the flow map.
alpar@414
    81
    static FlowMap* createFlowMap(const Digraph& digraph) {
alpar@414
    82
      return new FlowMap(digraph);
alpar@414
    83
    }
alpar@414
    84
alpar@414
    85
    /// \brief The eleavator type used by Circulation algorithm.
alpar@414
    86
    ///
alpar@414
    87
    /// The elevator type used by Circulation algorithm.
alpar@414
    88
    ///
alpar@414
    89
    /// \sa Elevator
alpar@414
    90
    /// \sa LinkedElevator
alpar@414
    91
    typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
alpar@414
    92
alpar@414
    93
    /// \brief Instantiates an Elevator.
alpar@414
    94
    ///
alpar@414
    95
    /// This function instantiates a \ref Elevator.
alpar@414
    96
    /// \param digraph The digraph, to which we would like to define
alpar@414
    97
    /// the elevator.
alpar@414
    98
    /// \param max_level The maximum level of the elevator.
alpar@414
    99
    static Elevator* createElevator(const Digraph& digraph, int max_level) {
alpar@414
   100
      return new Elevator(digraph, max_level);
alpar@414
   101
    }
alpar@414
   102
alpar@414
   103
    /// \brief The tolerance used by the algorithm
alpar@414
   104
    ///
alpar@414
   105
    /// The tolerance used by the algorithm to handle inexact computation.
alpar@414
   106
    typedef lemon::Tolerance<Value> Tolerance;
alpar@414
   107
alpar@414
   108
  };
alpar@414
   109
alpar@414
   110
  ///Push-relabel algorithm for the Network Circulation Problem.
alpar@414
   111
alpar@414
   112
  /**
alpar@414
   113
     \ingroup max_flow
alpar@414
   114
     This class implements a push-relabel algorithm
alpar@414
   115
     or the Network Circulation Problem.
alpar@414
   116
     The exact formulation of this problem is the following.
alpar@414
   117
     \f[\sum_{e\in\rho(v)}x(e)-\sum_{e\in\delta(v)}x(e)\leq
alpar@414
   118
     -delta(v)\quad \forall v\in V \f]
alpar@414
   119
     \f[ lo(e)\leq x(e) \leq up(e) \quad \forall e\in E \f]
alpar@414
   120
  */
alpar@414
   121
  template<class _Graph,
alpar@414
   122
           class _LCapMap=typename _Graph::template ArcMap<int>,
alpar@414
   123
           class _UCapMap=_LCapMap,
alpar@414
   124
           class _DeltaMap=typename _Graph::template NodeMap<
alpar@414
   125
             typename _UCapMap::Value>,
alpar@414
   126
           class _Traits=CirculationDefaultTraits<_Graph, _LCapMap,
alpar@414
   127
                                                  _UCapMap, _DeltaMap> >
alpar@414
   128
  class Circulation {
alpar@414
   129
alpar@414
   130
    typedef _Traits Traits;
alpar@414
   131
    typedef typename Traits::Digraph Digraph;
alpar@414
   132
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
alpar@414
   133
alpar@414
   134
    typedef typename Traits::Value Value;
alpar@414
   135
alpar@414
   136
    typedef typename Traits::LCapMap LCapMap;
alpar@414
   137
    typedef typename Traits::UCapMap UCapMap;
alpar@414
   138
    typedef typename Traits::DeltaMap DeltaMap;
alpar@414
   139
    typedef typename Traits::FlowMap FlowMap;
alpar@414
   140
    typedef typename Traits::Elevator Elevator;
alpar@414
   141
    typedef typename Traits::Tolerance Tolerance;
alpar@414
   142
alpar@414
   143
    typedef typename Digraph::template NodeMap<Value> ExcessMap;
alpar@414
   144
alpar@414
   145
    const Digraph &_g;
alpar@414
   146
    int _node_num;
alpar@414
   147
alpar@414
   148
    const LCapMap *_lo;
alpar@414
   149
    const UCapMap *_up;
alpar@414
   150
    const DeltaMap *_delta;
alpar@414
   151
alpar@414
   152
    FlowMap *_flow;
alpar@414
   153
    bool _local_flow;
alpar@414
   154
alpar@414
   155
    Elevator* _level;
alpar@414
   156
    bool _local_level;
alpar@414
   157
alpar@414
   158
    ExcessMap* _excess;
alpar@414
   159
alpar@414
   160
    Tolerance _tol;
alpar@414
   161
    int _el;
alpar@414
   162
alpar@414
   163
  public:
alpar@414
   164
alpar@414
   165
    typedef Circulation Create;
alpar@414
   166
alpar@414
   167
    ///\name Named template parameters
alpar@414
   168
alpar@414
   169
    ///@{
alpar@414
   170
alpar@414
   171
    template <typename _FlowMap>
alpar@414
   172
    struct DefFlowMapTraits : public Traits {
alpar@414
   173
      typedef _FlowMap FlowMap;
alpar@414
   174
      static FlowMap *createFlowMap(const Digraph&) {
alpar@414
   175
        LEMON_ASSERT(false, "FlowMap is not initialized");
alpar@414
   176
        return 0; // ignore warnings
alpar@414
   177
      }
alpar@414
   178
    };
alpar@414
   179
alpar@414
   180
    /// \brief \ref named-templ-param "Named parameter" for setting
alpar@414
   181
    /// FlowMap type
alpar@414
   182
    ///
alpar@414
   183
    /// \ref named-templ-param "Named parameter" for setting FlowMap
alpar@414
   184
    /// type
alpar@414
   185
    template <typename _FlowMap>
alpar@414
   186
    struct DefFlowMap
alpar@414
   187
      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
alpar@414
   188
                           DefFlowMapTraits<_FlowMap> > {
alpar@414
   189
      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
alpar@414
   190
                          DefFlowMapTraits<_FlowMap> > Create;
alpar@414
   191
    };
alpar@414
   192
alpar@414
   193
    template <typename _Elevator>
alpar@414
   194
    struct DefElevatorTraits : public Traits {
alpar@414
   195
      typedef _Elevator Elevator;
alpar@414
   196
      static Elevator *createElevator(const Digraph&, int) {
alpar@414
   197
        LEMON_ASSERT(false, "Elevator is not initialized");
alpar@414
   198
        return 0; // ignore warnings
alpar@414
   199
      }
alpar@414
   200
    };
alpar@414
   201
alpar@414
   202
    /// \brief \ref named-templ-param "Named parameter" for setting
alpar@414
   203
    /// Elevator type
alpar@414
   204
    ///
alpar@414
   205
    /// \ref named-templ-param "Named parameter" for setting Elevator
alpar@414
   206
    /// type
alpar@414
   207
    template <typename _Elevator>
alpar@414
   208
    struct DefElevator
alpar@414
   209
      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
alpar@414
   210
                           DefElevatorTraits<_Elevator> > {
alpar@414
   211
      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
alpar@414
   212
                          DefElevatorTraits<_Elevator> > Create;
alpar@414
   213
    };
alpar@414
   214
alpar@414
   215
    template <typename _Elevator>
alpar@414
   216
    struct DefStandardElevatorTraits : public Traits {
alpar@414
   217
      typedef _Elevator Elevator;
alpar@414
   218
      static Elevator *createElevator(const Digraph& digraph, int max_level) {
alpar@414
   219
        return new Elevator(digraph, max_level);
alpar@414
   220
      }
alpar@414
   221
    };
alpar@414
   222
alpar@414
   223
    /// \brief \ref named-templ-param "Named parameter" for setting
alpar@414
   224
    /// Elevator type
alpar@414
   225
    ///
alpar@414
   226
    /// \ref named-templ-param "Named parameter" for setting Elevator
alpar@414
   227
    /// type. The Elevator should be standard constructor interface, ie.
alpar@414
   228
    /// the digraph and the maximum level should be passed to it.
alpar@414
   229
    template <typename _Elevator>
alpar@414
   230
    struct DefStandardElevator
alpar@414
   231
      : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
alpar@414
   232
                       DefStandardElevatorTraits<_Elevator> > {
alpar@414
   233
      typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
alpar@414
   234
                      DefStandardElevatorTraits<_Elevator> > Create;
alpar@414
   235
    };
alpar@414
   236
alpar@414
   237
    /// @}
alpar@414
   238
alpar@414
   239
  protected:
alpar@414
   240
alpar@414
   241
    Circulation() {}
alpar@414
   242
alpar@414
   243
  public:
alpar@414
   244
alpar@414
   245
    /// The constructor of the class.
alpar@414
   246
alpar@414
   247
    /// The constructor of the class.
alpar@414
   248
    /// \param g The digraph the algorithm runs on.
alpar@414
   249
    /// \param lo The lower bound capacity of the arcs.
alpar@414
   250
    /// \param up The upper bound capacity of the arcs.
alpar@414
   251
    /// \param delta The lower bound on node excess.
alpar@414
   252
    Circulation(const Digraph &g,const LCapMap &lo,
alpar@414
   253
                const UCapMap &up,const DeltaMap &delta)
alpar@414
   254
      : _g(g), _node_num(),
alpar@414
   255
        _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
alpar@414
   256
        _level(0), _local_level(false), _excess(0), _el() {}
alpar@414
   257
alpar@414
   258
    /// Destrcutor.
alpar@414
   259
    ~Circulation() {
alpar@414
   260
      destroyStructures();
alpar@414
   261
    }
alpar@414
   262
alpar@414
   263
  private:
alpar@414
   264
alpar@414
   265
    void createStructures() {
alpar@414
   266
      _node_num = _el = countNodes(_g);
alpar@414
   267
alpar@414
   268
      if (!_flow) {
alpar@414
   269
        _flow = Traits::createFlowMap(_g);
alpar@414
   270
        _local_flow = true;
alpar@414
   271
      }
alpar@414
   272
      if (!_level) {
alpar@414
   273
        _level = Traits::createElevator(_g, _node_num);
alpar@414
   274
        _local_level = true;
alpar@414
   275
      }
alpar@414
   276
      if (!_excess) {
alpar@414
   277
        _excess = new ExcessMap(_g);
alpar@414
   278
      }
alpar@414
   279
    }
alpar@414
   280
alpar@414
   281
    void destroyStructures() {
alpar@414
   282
      if (_local_flow) {
alpar@414
   283
        delete _flow;
alpar@414
   284
      }
alpar@414
   285
      if (_local_level) {
alpar@414
   286
        delete _level;
alpar@414
   287
      }
alpar@414
   288
      if (_excess) {
alpar@414
   289
        delete _excess;
alpar@414
   290
      }
alpar@414
   291
    }
alpar@414
   292
alpar@414
   293
  public:
alpar@414
   294
alpar@414
   295
    /// Sets the lower bound capacity map.
alpar@414
   296
alpar@414
   297
    /// Sets the lower bound capacity map.
alpar@414
   298
    /// \return \c (*this)
alpar@414
   299
    Circulation& lowerCapMap(const LCapMap& map) {
alpar@414
   300
      _lo = &map;
alpar@414
   301
      return *this;
alpar@414
   302
    }
alpar@414
   303
alpar@414
   304
    /// Sets the upper bound capacity map.
alpar@414
   305
alpar@414
   306
    /// Sets the upper bound capacity map.
alpar@414
   307
    /// \return \c (*this)
alpar@414
   308
    Circulation& upperCapMap(const LCapMap& map) {
alpar@414
   309
      _up = &map;
alpar@414
   310
      return *this;
alpar@414
   311
    }
alpar@414
   312
alpar@414
   313
    /// Sets the lower bound map on excess.
alpar@414
   314
alpar@414
   315
    /// Sets the lower bound map on excess.
alpar@414
   316
    /// \return \c (*this)
alpar@414
   317
    Circulation& deltaMap(const DeltaMap& map) {
alpar@414
   318
      _delta = &map;
alpar@414
   319
      return *this;
alpar@414
   320
    }
alpar@414
   321
alpar@414
   322
    /// Sets the flow map.
alpar@414
   323
alpar@414
   324
    /// Sets the flow map.
alpar@414
   325
    /// \return \c (*this)
alpar@414
   326
    Circulation& flowMap(FlowMap& map) {
alpar@414
   327
      if (_local_flow) {
alpar@414
   328
        delete _flow;
alpar@414
   329
        _local_flow = false;
alpar@414
   330
      }
alpar@414
   331
      _flow = &map;
alpar@414
   332
      return *this;
alpar@414
   333
    }
alpar@414
   334
alpar@414
   335
    /// Returns the flow map.
alpar@414
   336
alpar@414
   337
    /// \return The flow map.
alpar@414
   338
    ///
alpar@414
   339
    const FlowMap& flowMap() {
alpar@414
   340
      return *_flow;
alpar@414
   341
    }
alpar@414
   342
alpar@414
   343
    /// Sets the elevator.
alpar@414
   344
alpar@414
   345
    /// Sets the elevator.
alpar@414
   346
    /// \return \c (*this)
alpar@414
   347
    Circulation& elevator(Elevator& elevator) {
alpar@414
   348
      if (_local_level) {
alpar@414
   349
        delete _level;
alpar@414
   350
        _local_level = false;
alpar@414
   351
      }
alpar@414
   352
      _level = &elevator;
alpar@414
   353
      return *this;
alpar@414
   354
    }
alpar@414
   355
alpar@414
   356
    /// Returns the elevator.
alpar@414
   357
alpar@414
   358
    /// \return The elevator.
alpar@414
   359
    ///
alpar@414
   360
    const Elevator& elevator() {
alpar@414
   361
      return *_level;
alpar@414
   362
    }
alpar@414
   363
alpar@414
   364
    /// Sets the tolerance used by algorithm.
alpar@414
   365
alpar@414
   366
    /// Sets the tolerance used by algorithm.
alpar@414
   367
    ///
alpar@414
   368
    Circulation& tolerance(const Tolerance& tolerance) const {
alpar@414
   369
      _tol = tolerance;
alpar@414
   370
      return *this;
alpar@414
   371
    }
alpar@414
   372
alpar@414
   373
    /// Returns the tolerance used by algorithm.
alpar@414
   374
alpar@414
   375
    /// Returns the tolerance used by algorithm.
alpar@414
   376
    ///
alpar@414
   377
    const Tolerance& tolerance() const {
alpar@414
   378
      return tolerance;
alpar@414
   379
    }
alpar@414
   380
alpar@414
   381
    /// \name Execution control
alpar@414
   382
    /// The simplest way to execute the algorithm is to use one of the
alpar@414
   383
    /// member functions called \c run().
alpar@414
   384
    /// \n
alpar@414
   385
    /// If you need more control on initial solution or execution then
alpar@414
   386
    /// you have to call one \ref init() function and then the start()
alpar@414
   387
    /// function.
alpar@414
   388
alpar@414
   389
    ///@{
alpar@414
   390
alpar@414
   391
    /// Initializes the internal data structures.
alpar@414
   392
alpar@414
   393
    /// Initializes the internal data structures. This function sets
alpar@414
   394
    /// all flow values to the lower bound.
alpar@414
   395
    /// \return This function returns false if the initialization
alpar@414
   396
    /// process found a barrier.
alpar@414
   397
    void init()
alpar@414
   398
    {
alpar@414
   399
      createStructures();
alpar@414
   400
alpar@414
   401
      for(NodeIt n(_g);n!=INVALID;++n) {
alpar@414
   402
        _excess->set(n, (*_delta)[n]);
alpar@414
   403
      }
alpar@414
   404
alpar@414
   405
      for (ArcIt e(_g);e!=INVALID;++e) {
alpar@414
   406
        _flow->set(e, (*_lo)[e]);
alpar@414
   407
        _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
alpar@414
   408
        _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
alpar@414
   409
      }
alpar@414
   410
alpar@414
   411
      // global relabeling tested, but in general case it provides
alpar@414
   412
      // worse performance for random digraphs
alpar@414
   413
      _level->initStart();
alpar@414
   414
      for(NodeIt n(_g);n!=INVALID;++n)
alpar@414
   415
        _level->initAddItem(n);
alpar@414
   416
      _level->initFinish();
alpar@414
   417
      for(NodeIt n(_g);n!=INVALID;++n)
alpar@414
   418
        if(_tol.positive((*_excess)[n]))
alpar@414
   419
          _level->activate(n);
alpar@414
   420
    }
alpar@414
   421
alpar@414
   422
    /// Initializes the internal data structures.
alpar@414
   423
alpar@414
   424
    /// Initializes the internal data structures. This functions uses
alpar@414
   425
    /// greedy approach to construct the initial solution.
alpar@414
   426
    void greedyInit()
alpar@414
   427
    {
alpar@414
   428
      createStructures();
alpar@414
   429
alpar@414
   430
      for(NodeIt n(_g);n!=INVALID;++n) {
alpar@414
   431
        _excess->set(n, (*_delta)[n]);
alpar@414
   432
      }
alpar@414
   433
alpar@414
   434
      for (ArcIt e(_g);e!=INVALID;++e) {
alpar@414
   435
        if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
alpar@414
   436
          _flow->set(e, (*_up)[e]);
alpar@414
   437
          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
alpar@414
   438
          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
alpar@414
   439
        } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
alpar@414
   440
          _flow->set(e, (*_lo)[e]);
alpar@414
   441
          _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_lo)[e]);
alpar@414
   442
          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_lo)[e]);
alpar@414
   443
        } else {
alpar@414
   444
          Value fc = -(*_excess)[_g.target(e)];
alpar@414
   445
          _flow->set(e, fc);
alpar@414
   446
          _excess->set(_g.target(e), 0);
alpar@414
   447
          _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
alpar@414
   448
        }
alpar@414
   449
      }
alpar@414
   450
alpar@414
   451
      _level->initStart();
alpar@414
   452
      for(NodeIt n(_g);n!=INVALID;++n)
alpar@414
   453
        _level->initAddItem(n);
alpar@414
   454
      _level->initFinish();
alpar@414
   455
      for(NodeIt n(_g);n!=INVALID;++n)
alpar@414
   456
        if(_tol.positive((*_excess)[n]))
alpar@414
   457
          _level->activate(n);
alpar@414
   458
    }
alpar@414
   459
alpar@414
   460
    ///Starts the algorithm
alpar@414
   461
alpar@414
   462
    ///This function starts the algorithm.
alpar@414
   463
    ///\return This function returns true if it found a feasible circulation.
alpar@414
   464
    ///
alpar@414
   465
    ///\sa barrier()
alpar@414
   466
    bool start()
alpar@414
   467
    {
alpar@414
   468
alpar@414
   469
      Node act;
alpar@414
   470
      Node bact=INVALID;
alpar@414
   471
      Node last_activated=INVALID;
alpar@414
   472
      while((act=_level->highestActive())!=INVALID) {
alpar@414
   473
        int actlevel=(*_level)[act];
alpar@414
   474
        int mlevel=_node_num;
alpar@414
   475
        Value exc=(*_excess)[act];
alpar@414
   476
alpar@414
   477
        for(OutArcIt e(_g,act);e!=INVALID; ++e) {
alpar@414
   478
          Node v = _g.target(e);
alpar@414
   479
          Value fc=(*_up)[e]-(*_flow)[e];
alpar@414
   480
          if(!_tol.positive(fc)) continue;
alpar@414
   481
          if((*_level)[v]<actlevel) {
alpar@414
   482
            if(!_tol.less(fc, exc)) {
alpar@414
   483
              _flow->set(e, (*_flow)[e] + exc);
alpar@414
   484
              _excess->set(v, (*_excess)[v] + exc);
alpar@414
   485
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
alpar@414
   486
                _level->activate(v);
alpar@414
   487
              _excess->set(act,0);
alpar@414
   488
              _level->deactivate(act);
alpar@414
   489
              goto next_l;
alpar@414
   490
            }
alpar@414
   491
            else {
alpar@414
   492
              _flow->set(e, (*_up)[e]);
alpar@414
   493
              _excess->set(v, (*_excess)[v] + fc);
alpar@414
   494
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
alpar@414
   495
                _level->activate(v);
alpar@414
   496
              exc-=fc;
alpar@414
   497
            }
alpar@414
   498
          }
alpar@414
   499
          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
alpar@414
   500
        }
alpar@414
   501
        for(InArcIt e(_g,act);e!=INVALID; ++e) {
alpar@414
   502
          Node v = _g.source(e);
alpar@414
   503
          Value fc=(*_flow)[e]-(*_lo)[e];
alpar@414
   504
          if(!_tol.positive(fc)) continue;
alpar@414
   505
          if((*_level)[v]<actlevel) {
alpar@414
   506
            if(!_tol.less(fc, exc)) {
alpar@414
   507
              _flow->set(e, (*_flow)[e] - exc);
alpar@414
   508
              _excess->set(v, (*_excess)[v] + exc);
alpar@414
   509
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
alpar@414
   510
                _level->activate(v);
alpar@414
   511
              _excess->set(act,0);
alpar@414
   512
              _level->deactivate(act);
alpar@414
   513
              goto next_l;
alpar@414
   514
            }
alpar@414
   515
            else {
alpar@414
   516
              _flow->set(e, (*_lo)[e]);
alpar@414
   517
              _excess->set(v, (*_excess)[v] + fc);
alpar@414
   518
              if(!_level->active(v) && _tol.positive((*_excess)[v]))
alpar@414
   519
                _level->activate(v);
alpar@414
   520
              exc-=fc;
alpar@414
   521
            }
alpar@414
   522
          }
alpar@414
   523
          else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
alpar@414
   524
        }
alpar@414
   525
alpar@414
   526
        _excess->set(act, exc);
alpar@414
   527
        if(!_tol.positive(exc)) _level->deactivate(act);
alpar@414
   528
        else if(mlevel==_node_num) {
alpar@414
   529
          _level->liftHighestActiveToTop();
alpar@414
   530
          _el = _node_num;
alpar@414
   531
          return false;
alpar@414
   532
        }
alpar@414
   533
        else {
alpar@414
   534
          _level->liftHighestActive(mlevel+1);
alpar@414
   535
          if(_level->onLevel(actlevel)==0) {
alpar@414
   536
            _el = actlevel;
alpar@414
   537
            return false;
alpar@414
   538
          }
alpar@414
   539
        }
alpar@414
   540
      next_l:
alpar@414
   541
        ;
alpar@414
   542
      }
alpar@414
   543
      return true;
alpar@414
   544
    }
alpar@414
   545
alpar@414
   546
    /// Runs the circulation algorithm.
alpar@414
   547
alpar@414
   548
    /// Runs the circulation algorithm.
alpar@414
   549
    /// \note fc.run() is just a shortcut of the following code.
alpar@414
   550
    /// \code
alpar@414
   551
    ///   fc.greedyInit();
alpar@414
   552
    ///   return fc.start();
alpar@414
   553
    /// \endcode
alpar@414
   554
    bool run() {
alpar@414
   555
      greedyInit();
alpar@414
   556
      return start();
alpar@414
   557
    }
alpar@414
   558
alpar@414
   559
    /// @}
alpar@414
   560
alpar@414
   561
    /// \name Query Functions
alpar@414
   562
    /// The result of the %Circulation algorithm can be obtained using
alpar@414
   563
    /// these functions.
alpar@414
   564
    /// \n
alpar@414
   565
    /// Before the use of these functions,
alpar@414
   566
    /// either run() or start() must be called.
alpar@414
   567
alpar@414
   568
    ///@{
alpar@414
   569
alpar@414
   570
    /**
alpar@414
   571
       \brief Returns a barrier
alpar@414
   572
       
alpar@414
   573
       Barrier is a set \e B of nodes for which
alpar@414
   574
       \f[ \sum_{v\in B}-delta(v)<
alpar@414
   575
       \sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
alpar@414
   576
       holds. The existence of a set with this property prooves that a feasible
alpar@414
   577
       flow cannot exists.
alpar@414
   578
       \sa checkBarrier()
alpar@414
   579
       \sa run()
alpar@414
   580
    */
alpar@414
   581
    template<class GT>
alpar@414
   582
    void barrierMap(GT &bar)
alpar@414
   583
    {
alpar@414
   584
      for(NodeIt n(_g);n!=INVALID;++n)
alpar@414
   585
        bar.set(n, (*_level)[n] >= _el);
alpar@414
   586
    }
alpar@414
   587
alpar@414
   588
    ///Returns true if the node is in the barrier
alpar@414
   589
alpar@414
   590
    ///Returns true if the node is in the barrier
alpar@414
   591
    ///\sa barrierMap()
alpar@414
   592
    bool barrier(const Node& node)
alpar@414
   593
    {
alpar@414
   594
      return (*_level)[node] >= _el;
alpar@414
   595
    }
alpar@414
   596
alpar@414
   597
    /// \brief Returns the flow on the arc.
alpar@414
   598
    ///
alpar@414
   599
    /// Sets the \c flowMap to the flow on the arcs. This method can
alpar@414
   600
    /// be called after the second phase of algorithm.
alpar@414
   601
    Value flow(const Arc& arc) const {
alpar@414
   602
      return (*_flow)[arc];
alpar@414
   603
    }
alpar@414
   604
alpar@414
   605
    /// @}
alpar@414
   606
alpar@414
   607
    /// \name Checker Functions
alpar@414
   608
    /// The feasibility  of the results can be checked using
alpar@414
   609
    /// these functions.
alpar@414
   610
    /// \n
alpar@414
   611
    /// Before the use of these functions,
alpar@414
   612
    /// either run() or start() must be called.
alpar@414
   613
alpar@414
   614
    ///@{
alpar@414
   615
alpar@414
   616
    ///Check if the  \c flow is a feasible circulation
alpar@414
   617
    bool checkFlow() {
alpar@414
   618
      for(ArcIt e(_g);e!=INVALID;++e)
alpar@414
   619
        if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
alpar@414
   620
      for(NodeIt n(_g);n!=INVALID;++n)
alpar@414
   621
        {
alpar@414
   622
          Value dif=-(*_delta)[n];
alpar@414
   623
          for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
alpar@414
   624
          for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
alpar@414
   625
          if(_tol.negative(dif)) return false;
alpar@414
   626
        }
alpar@414
   627
      return true;
alpar@414
   628
    }
alpar@414
   629
alpar@414
   630
    ///Check whether or not the last execution provides a barrier
alpar@414
   631
alpar@414
   632
    ///Check whether or not the last execution provides a barrier
alpar@414
   633
    ///\sa barrier()
alpar@414
   634
    bool checkBarrier()
alpar@414
   635
    {
alpar@414
   636
      Value delta=0;
alpar@414
   637
      for(NodeIt n(_g);n!=INVALID;++n)
alpar@414
   638
        if(barrier(n))
alpar@414
   639
          delta-=(*_delta)[n];
alpar@414
   640
      for(ArcIt e(_g);e!=INVALID;++e)
alpar@414
   641
        {
alpar@414
   642
          Node s=_g.source(e);
alpar@414
   643
          Node t=_g.target(e);
alpar@414
   644
          if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
alpar@414
   645
          else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
alpar@414
   646
        }
alpar@414
   647
      return _tol.negative(delta);
alpar@414
   648
    }
alpar@414
   649
alpar@414
   650
    /// @}
alpar@414
   651
alpar@414
   652
  };
alpar@414
   653
alpar@414
   654
}
alpar@414
   655
alpar@414
   656
#endif