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