1.1 --- a/lemon/circulation.h Mon Jan 12 23:11:39 2009 +0100
1.2 +++ b/lemon/circulation.h Thu Nov 05 15:48:01 2009 +0100
1.3 @@ -21,6 +21,7 @@
1.4
1.5 #include <lemon/tolerance.h>
1.6 #include <lemon/elevator.h>
1.7 +#include <limits>
1.8
1.9 ///\ingroup max_flow
1.10 ///\file
1.11 @@ -31,52 +32,56 @@
1.12 /// \brief Default traits class of Circulation class.
1.13 ///
1.14 /// Default traits class of Circulation class.
1.15 - /// \tparam _Diraph Digraph type.
1.16 - /// \tparam _LCapMap Lower bound capacity map type.
1.17 - /// \tparam _UCapMap Upper bound capacity map type.
1.18 - /// \tparam _DeltaMap Delta map type.
1.19 - template <typename _Diraph, typename _LCapMap,
1.20 - typename _UCapMap, typename _DeltaMap>
1.21 + ///
1.22 + /// \tparam GR Type of the digraph the algorithm runs on.
1.23 + /// \tparam LM The type of the lower bound map.
1.24 + /// \tparam UM The type of the upper bound (capacity) map.
1.25 + /// \tparam SM The type of the supply map.
1.26 + template <typename GR, typename LM,
1.27 + typename UM, typename SM>
1.28 struct CirculationDefaultTraits {
1.29
1.30 /// \brief The type of the digraph the algorithm runs on.
1.31 - typedef _Diraph Digraph;
1.32 + typedef GR Digraph;
1.33
1.34 - /// \brief The type of the map that stores the circulation lower
1.35 - /// bound.
1.36 + /// \brief The type of the lower bound map.
1.37 ///
1.38 - /// The type of the map that stores the circulation lower bound.
1.39 - /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
1.40 - typedef _LCapMap LCapMap;
1.41 + /// The type of the map that stores the lower bounds on the arcs.
1.42 + /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
1.43 + typedef LM LowerMap;
1.44
1.45 - /// \brief The type of the map that stores the circulation upper
1.46 - /// bound.
1.47 + /// \brief The type of the upper bound (capacity) map.
1.48 ///
1.49 - /// The type of the map that stores the circulation upper bound.
1.50 - /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
1.51 - typedef _UCapMap UCapMap;
1.52 + /// The type of the map that stores the upper bounds (capacities)
1.53 + /// on the arcs.
1.54 + /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
1.55 + typedef UM UpperMap;
1.56
1.57 - /// \brief The type of the map that stores the lower bound for
1.58 - /// the supply of the nodes.
1.59 + /// \brief The type of supply map.
1.60 ///
1.61 - /// The type of the map that stores the lower bound for the supply
1.62 - /// of the nodes. It must meet the \ref concepts::ReadMap "ReadMap"
1.63 - /// concept.
1.64 - typedef _DeltaMap DeltaMap;
1.65 + /// The type of the map that stores the signed supply values of the
1.66 + /// nodes.
1.67 + /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
1.68 + typedef SM SupplyMap;
1.69
1.70 - /// \brief The type of the flow values.
1.71 - typedef typename DeltaMap::Value Value;
1.72 + /// \brief The type of the flow and supply values.
1.73 + typedef typename SupplyMap::Value Value;
1.74
1.75 /// \brief The type of the map that stores the flow values.
1.76 ///
1.77 /// The type of the map that stores the flow values.
1.78 - /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
1.79 + /// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap"
1.80 + /// concept.
1.81 +#ifdef DOXYGEN
1.82 + typedef GR::ArcMap<Value> FlowMap;
1.83 +#else
1.84 typedef typename Digraph::template ArcMap<Value> FlowMap;
1.85 +#endif
1.86
1.87 /// \brief Instantiates a FlowMap.
1.88 ///
1.89 /// This function instantiates a \ref FlowMap.
1.90 - /// \param digraph The digraph, to which we would like to define
1.91 + /// \param digraph The digraph for which we would like to define
1.92 /// the flow map.
1.93 static FlowMap* createFlowMap(const Digraph& digraph) {
1.94 return new FlowMap(digraph);
1.95 @@ -86,14 +91,17 @@
1.96 ///
1.97 /// The elevator type used by the algorithm.
1.98 ///
1.99 - /// \sa Elevator
1.100 - /// \sa LinkedElevator
1.101 + /// \sa Elevator, LinkedElevator
1.102 +#ifdef DOXYGEN
1.103 + typedef lemon::Elevator<GR, GR::Node> Elevator;
1.104 +#else
1.105 typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
1.106 +#endif
1.107
1.108 /// \brief Instantiates an Elevator.
1.109 ///
1.110 /// This function instantiates an \ref Elevator.
1.111 - /// \param digraph The digraph, to which we would like to define
1.112 + /// \param digraph The digraph for which we would like to define
1.113 /// the elevator.
1.114 /// \param max_level The maximum level of the elevator.
1.115 static Elevator* createElevator(const Digraph& digraph, int max_level) {
1.116 @@ -111,73 +119,90 @@
1.117 \brief Push-relabel algorithm for the network circulation problem.
1.118
1.119 \ingroup max_flow
1.120 - This class implements a push-relabel algorithm for the network
1.121 - circulation problem.
1.122 + This class implements a push-relabel algorithm for the \e network
1.123 + \e circulation problem.
1.124 It is to find a feasible circulation when lower and upper bounds
1.125 - are given for the flow values on the arcs and lower bounds
1.126 - are given for the supply values of the nodes.
1.127 + are given for the flow values on the arcs and lower bounds are
1.128 + given for the difference between the outgoing and incoming flow
1.129 + at the nodes.
1.130
1.131 The exact formulation of this problem is the following.
1.132 - Let \f$G=(V,A)\f$ be a digraph,
1.133 - \f$lower, upper: A\rightarrow\mathbf{R}^+_0\f$,
1.134 - \f$delta: V\rightarrow\mathbf{R}\f$. Find a feasible circulation
1.135 - \f$f: A\rightarrow\mathbf{R}^+_0\f$ so that
1.136 - \f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a)
1.137 - \geq delta(v) \quad \forall v\in V, \f]
1.138 - \f[ lower(a)\leq f(a) \leq upper(a) \quad \forall a\in A. \f]
1.139 - \note \f$delta(v)\f$ specifies a lower bound for the supply of node
1.140 - \f$v\f$. It can be either positive or negative, however note that
1.141 - \f$\sum_{v\in V}delta(v)\f$ should be zero or negative in order to
1.142 - have a feasible solution.
1.143 + Let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$
1.144 + \f$upper: A\rightarrow\mathbf{R}\cup\{\infty\}\f$ denote the lower and
1.145 + upper bounds on the arcs, for which \f$lower(uv) \leq upper(uv)\f$
1.146 + holds for all \f$uv\in A\f$, and \f$sup: V\rightarrow\mathbf{R}\f$
1.147 + denotes the signed supply values of the nodes.
1.148 + If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
1.149 + supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
1.150 + \f$-sup(u)\f$ demand.
1.151 + A feasible circulation is an \f$f: A\rightarrow\mathbf{R}\f$
1.152 + solution of the following problem.
1.153
1.154 - \note A special case of this problem is when
1.155 - \f$\sum_{v\in V}delta(v) = 0\f$. Then the supply of each node \f$v\f$
1.156 - will be \e equal \e to \f$delta(v)\f$, if a circulation can be found.
1.157 - Thus a feasible solution for the
1.158 - \ref min_cost_flow "minimum cost flow" problem can be calculated
1.159 - in this way.
1.160 + \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu)
1.161 + \geq sup(u) \quad \forall u\in V, \f]
1.162 + \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A. \f]
1.163 +
1.164 + The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
1.165 + zero or negative in order to have a feasible solution (since the sum
1.166 + of the expressions on the left-hand side of the inequalities is zero).
1.167 + It means that the total demand must be greater or equal to the total
1.168 + supply and all the supplies have to be carried out from the supply nodes,
1.169 + but there could be demands that are not satisfied.
1.170 + If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
1.171 + constraints have to be satisfied with equality, i.e. all demands
1.172 + have to be satisfied and all supplies have to be used.
1.173 +
1.174 + If you need the opposite inequalities in the supply/demand constraints
1.175 + (i.e. the total demand is less than the total supply and all the demands
1.176 + have to be satisfied while there could be supplies that are not used),
1.177 + then you could easily transform the problem to the above form by reversing
1.178 + the direction of the arcs and taking the negative of the supply values
1.179 + (e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
1.180
1.181 - \tparam _Digraph The type of the digraph the algorithm runs on.
1.182 - \tparam _LCapMap The type of the lower bound capacity map. The default
1.183 - map type is \ref concepts::Digraph::ArcMap "_Digraph::ArcMap<int>".
1.184 - \tparam _UCapMap The type of the upper bound capacity map. The default
1.185 - map type is \c _LCapMap.
1.186 - \tparam _DeltaMap The type of the map that stores the lower bound
1.187 - for the supply of the nodes. The default map type is
1.188 - \c _Digraph::ArcMap<_UCapMap::Value>.
1.189 + This algorithm either calculates a feasible circulation, or provides
1.190 + a \ref barrier() "barrier", which prooves that a feasible soultion
1.191 + cannot exist.
1.192 +
1.193 + Note that this algorithm also provides a feasible solution for the
1.194 + \ref min_cost_flow "minimum cost flow problem".
1.195 +
1.196 + \tparam GR The type of the digraph the algorithm runs on.
1.197 + \tparam LM The type of the lower bound map. The default
1.198 + map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".
1.199 + \tparam UM The type of the upper bound (capacity) map.
1.200 + The default map type is \c LM.
1.201 + \tparam SM The type of the supply map. The default map type is
1.202 + \ref concepts::Digraph::NodeMap "GR::NodeMap<UM::Value>".
1.203 */
1.204 #ifdef DOXYGEN
1.205 -template< typename _Digraph,
1.206 - typename _LCapMap,
1.207 - typename _UCapMap,
1.208 - typename _DeltaMap,
1.209 - typename _Traits >
1.210 +template< typename GR,
1.211 + typename LM,
1.212 + typename UM,
1.213 + typename SM,
1.214 + typename TR >
1.215 #else
1.216 -template< typename _Digraph,
1.217 - typename _LCapMap = typename _Digraph::template ArcMap<int>,
1.218 - typename _UCapMap = _LCapMap,
1.219 - typename _DeltaMap = typename _Digraph::
1.220 - template NodeMap<typename _UCapMap::Value>,
1.221 - typename _Traits=CirculationDefaultTraits<_Digraph, _LCapMap,
1.222 - _UCapMap, _DeltaMap> >
1.223 +template< typename GR,
1.224 + typename LM = typename GR::template ArcMap<int>,
1.225 + typename UM = LM,
1.226 + typename SM = typename GR::template NodeMap<typename UM::Value>,
1.227 + typename TR = CirculationDefaultTraits<GR, LM, UM, SM> >
1.228 #endif
1.229 class Circulation {
1.230 public:
1.231
1.232 ///The \ref CirculationDefaultTraits "traits class" of the algorithm.
1.233 - typedef _Traits Traits;
1.234 + typedef TR Traits;
1.235 ///The type of the digraph the algorithm runs on.
1.236 typedef typename Traits::Digraph Digraph;
1.237 - ///The type of the flow values.
1.238 + ///The type of the flow and supply values.
1.239 typedef typename Traits::Value Value;
1.240
1.241 - /// The type of the lower bound capacity map.
1.242 - typedef typename Traits::LCapMap LCapMap;
1.243 - /// The type of the upper bound capacity map.
1.244 - typedef typename Traits::UCapMap UCapMap;
1.245 - /// \brief The type of the map that stores the lower bound for
1.246 - /// the supply of the nodes.
1.247 - typedef typename Traits::DeltaMap DeltaMap;
1.248 + ///The type of the lower bound map.
1.249 + typedef typename Traits::LowerMap LowerMap;
1.250 + ///The type of the upper bound (capacity) map.
1.251 + typedef typename Traits::UpperMap UpperMap;
1.252 + ///The type of the supply map.
1.253 + typedef typename Traits::SupplyMap SupplyMap;
1.254 ///The type of the flow map.
1.255 typedef typename Traits::FlowMap FlowMap;
1.256
1.257 @@ -193,9 +218,9 @@
1.258 const Digraph &_g;
1.259 int _node_num;
1.260
1.261 - const LCapMap *_lo;
1.262 - const UCapMap *_up;
1.263 - const DeltaMap *_delta;
1.264 + const LowerMap *_lo;
1.265 + const UpperMap *_up;
1.266 + const SupplyMap *_supply;
1.267
1.268 FlowMap *_flow;
1.269 bool _local_flow;
1.270 @@ -217,9 +242,9 @@
1.271
1.272 ///@{
1.273
1.274 - template <typename _FlowMap>
1.275 + template <typename T>
1.276 struct SetFlowMapTraits : public Traits {
1.277 - typedef _FlowMap FlowMap;
1.278 + typedef T FlowMap;
1.279 static FlowMap *createFlowMap(const Digraph&) {
1.280 LEMON_ASSERT(false, "FlowMap is not initialized");
1.281 return 0; // ignore warnings
1.282 @@ -231,17 +256,17 @@
1.283 ///
1.284 /// \ref named-templ-param "Named parameter" for setting FlowMap
1.285 /// type.
1.286 - template <typename _FlowMap>
1.287 + template <typename T>
1.288 struct SetFlowMap
1.289 - : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
1.290 - SetFlowMapTraits<_FlowMap> > {
1.291 - typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
1.292 - SetFlowMapTraits<_FlowMap> > Create;
1.293 + : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
1.294 + SetFlowMapTraits<T> > {
1.295 + typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
1.296 + SetFlowMapTraits<T> > Create;
1.297 };
1.298
1.299 - template <typename _Elevator>
1.300 + template <typename T>
1.301 struct SetElevatorTraits : public Traits {
1.302 - typedef _Elevator Elevator;
1.303 + typedef T Elevator;
1.304 static Elevator *createElevator(const Digraph&, int) {
1.305 LEMON_ASSERT(false, "Elevator is not initialized");
1.306 return 0; // ignore warnings
1.307 @@ -257,17 +282,17 @@
1.308 /// \ref elevator(Elevator&) "elevator()" function before calling
1.309 /// \ref run() or \ref init().
1.310 /// \sa SetStandardElevator
1.311 - template <typename _Elevator>
1.312 + template <typename T>
1.313 struct SetElevator
1.314 - : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
1.315 - SetElevatorTraits<_Elevator> > {
1.316 - typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
1.317 - SetElevatorTraits<_Elevator> > Create;
1.318 + : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
1.319 + SetElevatorTraits<T> > {
1.320 + typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
1.321 + SetElevatorTraits<T> > Create;
1.322 };
1.323
1.324 - template <typename _Elevator>
1.325 + template <typename T>
1.326 struct SetStandardElevatorTraits : public Traits {
1.327 - typedef _Elevator Elevator;
1.328 + typedef T Elevator;
1.329 static Elevator *createElevator(const Digraph& digraph, int max_level) {
1.330 return new Elevator(digraph, max_level);
1.331 }
1.332 @@ -285,12 +310,12 @@
1.333 /// algorithm with the \ref elevator(Elevator&) "elevator()" function
1.334 /// before calling \ref run() or \ref init().
1.335 /// \sa SetElevator
1.336 - template <typename _Elevator>
1.337 + template <typename T>
1.338 struct SetStandardElevator
1.339 - : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
1.340 - SetStandardElevatorTraits<_Elevator> > {
1.341 - typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
1.342 - SetStandardElevatorTraits<_Elevator> > Create;
1.343 + : public Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
1.344 + SetStandardElevatorTraits<T> > {
1.345 + typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap,
1.346 + SetStandardElevatorTraits<T> > Create;
1.347 };
1.348
1.349 /// @}
1.350 @@ -301,18 +326,20 @@
1.351
1.352 public:
1.353
1.354 - /// The constructor of the class.
1.355 + /// Constructor.
1.356
1.357 /// The constructor of the class.
1.358 - /// \param g The digraph the algorithm runs on.
1.359 - /// \param lo The lower bound capacity of the arcs.
1.360 - /// \param up The upper bound capacity of the arcs.
1.361 - /// \param delta The lower bound for the supply of the nodes.
1.362 - Circulation(const Digraph &g,const LCapMap &lo,
1.363 - const UCapMap &up,const DeltaMap &delta)
1.364 - : _g(g), _node_num(),
1.365 - _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
1.366 - _level(0), _local_level(false), _excess(0), _el() {}
1.367 + ///
1.368 + /// \param graph The digraph the algorithm runs on.
1.369 + /// \param lower The lower bounds for the flow values on the arcs.
1.370 + /// \param upper The upper bounds (capacities) for the flow values
1.371 + /// on the arcs.
1.372 + /// \param supply The signed supply values of the nodes.
1.373 + Circulation(const Digraph &graph, const LowerMap &lower,
1.374 + const UpperMap &upper, const SupplyMap &supply)
1.375 + : _g(graph), _lo(&lower), _up(&upper), _supply(&supply),
1.376 + _flow(NULL), _local_flow(false), _level(NULL), _local_level(false),
1.377 + _excess(NULL) {}
1.378
1.379 /// Destructor.
1.380 ~Circulation() {
1.381 @@ -322,6 +349,13 @@
1.382
1.383 private:
1.384
1.385 + bool checkBoundMaps() {
1.386 + for (ArcIt e(_g);e!=INVALID;++e) {
1.387 + if (_tol.less((*_up)[e], (*_lo)[e])) return false;
1.388 + }
1.389 + return true;
1.390 + }
1.391 +
1.392 void createStructures() {
1.393 _node_num = _el = countNodes(_g);
1.394
1.395 @@ -352,30 +386,30 @@
1.396
1.397 public:
1.398
1.399 - /// Sets the lower bound capacity map.
1.400 + /// Sets the lower bound map.
1.401
1.402 - /// Sets the lower bound capacity map.
1.403 + /// Sets the lower bound map.
1.404 /// \return <tt>(*this)</tt>
1.405 - Circulation& lowerCapMap(const LCapMap& map) {
1.406 + Circulation& lowerMap(const LowerMap& map) {
1.407 _lo = ↦
1.408 return *this;
1.409 }
1.410
1.411 - /// Sets the upper bound capacity map.
1.412 + /// Sets the upper bound (capacity) map.
1.413
1.414 - /// Sets the upper bound capacity map.
1.415 + /// Sets the upper bound (capacity) map.
1.416 /// \return <tt>(*this)</tt>
1.417 - Circulation& upperCapMap(const LCapMap& map) {
1.418 + Circulation& upperMap(const UpperMap& map) {
1.419 _up = ↦
1.420 return *this;
1.421 }
1.422
1.423 - /// Sets the lower bound map for the supply of the nodes.
1.424 + /// Sets the supply map.
1.425
1.426 - /// Sets the lower bound map for the supply of the nodes.
1.427 + /// Sets the supply map.
1.428 /// \return <tt>(*this)</tt>
1.429 - Circulation& deltaMap(const DeltaMap& map) {
1.430 - _delta = ↦
1.431 + Circulation& supplyMap(const SupplyMap& map) {
1.432 + _supply = ↦
1.433 return *this;
1.434 }
1.435
1.436 @@ -423,25 +457,27 @@
1.437 return *_level;
1.438 }
1.439
1.440 - /// \brief Sets the tolerance used by algorithm.
1.441 + /// \brief Sets the tolerance used by the algorithm.
1.442 ///
1.443 - /// Sets the tolerance used by algorithm.
1.444 - Circulation& tolerance(const Tolerance& tolerance) const {
1.445 + /// Sets the tolerance object used by the algorithm.
1.446 + /// \return <tt>(*this)</tt>
1.447 + Circulation& tolerance(const Tolerance& tolerance) {
1.448 _tol = tolerance;
1.449 return *this;
1.450 }
1.451
1.452 /// \brief Returns a const reference to the tolerance.
1.453 ///
1.454 - /// Returns a const reference to the tolerance.
1.455 + /// Returns a const reference to the tolerance object used by
1.456 + /// the algorithm.
1.457 const Tolerance& tolerance() const {
1.458 - return tolerance;
1.459 + return _tol;
1.460 }
1.461
1.462 /// \name Execution Control
1.463 /// The simplest way to execute the algorithm is to call \ref run().\n
1.464 - /// If you need more control on the initial solution or the execution,
1.465 - /// first you have to call one of the \ref init() functions, then
1.466 + /// If you need better control on the initial solution or the execution,
1.467 + /// you have to call one of the \ref init() functions first, then
1.468 /// the \ref start() function.
1.469
1.470 ///@{
1.471 @@ -452,16 +488,19 @@
1.472 /// to the lower bound.
1.473 void init()
1.474 {
1.475 + LEMON_DEBUG(checkBoundMaps(),
1.476 + "Upper bounds must be greater or equal to the lower bounds");
1.477 +
1.478 createStructures();
1.479
1.480 for(NodeIt n(_g);n!=INVALID;++n) {
1.481 - _excess->set(n, (*_delta)[n]);
1.482 + (*_excess)[n] = (*_supply)[n];
1.483 }
1.484
1.485 for (ArcIt e(_g);e!=INVALID;++e) {
1.486 _flow->set(e, (*_lo)[e]);
1.487 - _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
1.488 - _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
1.489 + (*_excess)[_g.target(e)] += (*_flow)[e];
1.490 + (*_excess)[_g.source(e)] -= (*_flow)[e];
1.491 }
1.492
1.493 // global relabeling tested, but in general case it provides
1.494 @@ -481,26 +520,29 @@
1.495 /// to construct the initial solution.
1.496 void greedyInit()
1.497 {
1.498 + LEMON_DEBUG(checkBoundMaps(),
1.499 + "Upper bounds must be greater or equal to the lower bounds");
1.500 +
1.501 createStructures();
1.502
1.503 for(NodeIt n(_g);n!=INVALID;++n) {
1.504 - _excess->set(n, (*_delta)[n]);
1.505 + (*_excess)[n] = (*_supply)[n];
1.506 }
1.507
1.508 for (ArcIt e(_g);e!=INVALID;++e) {
1.509 - if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
1.510 + if (!_tol.less(-(*_excess)[_g.target(e)], (*_up)[e])) {
1.511 _flow->set(e, (*_up)[e]);
1.512 - _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
1.513 - _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
1.514 - } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
1.515 + (*_excess)[_g.target(e)] += (*_up)[e];
1.516 + (*_excess)[_g.source(e)] -= (*_up)[e];
1.517 + } else if (_tol.less(-(*_excess)[_g.target(e)], (*_lo)[e])) {
1.518 _flow->set(e, (*_lo)[e]);
1.519 - _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_lo)[e]);
1.520 - _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_lo)[e]);
1.521 + (*_excess)[_g.target(e)] += (*_lo)[e];
1.522 + (*_excess)[_g.source(e)] -= (*_lo)[e];
1.523 } else {
1.524 Value fc = -(*_excess)[_g.target(e)];
1.525 _flow->set(e, fc);
1.526 - _excess->set(_g.target(e), 0);
1.527 - _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
1.528 + (*_excess)[_g.target(e)] = 0;
1.529 + (*_excess)[_g.source(e)] -= fc;
1.530 }
1.531 }
1.532
1.533 @@ -539,16 +581,16 @@
1.534 if((*_level)[v]<actlevel) {
1.535 if(!_tol.less(fc, exc)) {
1.536 _flow->set(e, (*_flow)[e] + exc);
1.537 - _excess->set(v, (*_excess)[v] + exc);
1.538 + (*_excess)[v] += exc;
1.539 if(!_level->active(v) && _tol.positive((*_excess)[v]))
1.540 _level->activate(v);
1.541 - _excess->set(act,0);
1.542 + (*_excess)[act] = 0;
1.543 _level->deactivate(act);
1.544 goto next_l;
1.545 }
1.546 else {
1.547 _flow->set(e, (*_up)[e]);
1.548 - _excess->set(v, (*_excess)[v] + fc);
1.549 + (*_excess)[v] += fc;
1.550 if(!_level->active(v) && _tol.positive((*_excess)[v]))
1.551 _level->activate(v);
1.552 exc-=fc;
1.553 @@ -563,16 +605,16 @@
1.554 if((*_level)[v]<actlevel) {
1.555 if(!_tol.less(fc, exc)) {
1.556 _flow->set(e, (*_flow)[e] - exc);
1.557 - _excess->set(v, (*_excess)[v] + exc);
1.558 + (*_excess)[v] += exc;
1.559 if(!_level->active(v) && _tol.positive((*_excess)[v]))
1.560 _level->activate(v);
1.561 - _excess->set(act,0);
1.562 + (*_excess)[act] = 0;
1.563 _level->deactivate(act);
1.564 goto next_l;
1.565 }
1.566 else {
1.567 _flow->set(e, (*_lo)[e]);
1.568 - _excess->set(v, (*_excess)[v] + fc);
1.569 + (*_excess)[v] += fc;
1.570 if(!_level->active(v) && _tol.positive((*_excess)[v]))
1.571 _level->activate(v);
1.572 exc-=fc;
1.573 @@ -581,7 +623,7 @@
1.574 else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
1.575 }
1.576
1.577 - _excess->set(act, exc);
1.578 + (*_excess)[act] = exc;
1.579 if(!_tol.positive(exc)) _level->deactivate(act);
1.580 else if(mlevel==_node_num) {
1.581 _level->liftHighestActiveToTop();
1.582 @@ -628,9 +670,9 @@
1.583
1.584 ///@{
1.585
1.586 - /// \brief Returns the flow on the given arc.
1.587 + /// \brief Returns the flow value on the given arc.
1.588 ///
1.589 - /// Returns the flow on the given arc.
1.590 + /// Returns the flow value on the given arc.
1.591 ///
1.592 /// \pre Either \ref run() or \ref init() must be called before
1.593 /// using this function.
1.594 @@ -653,8 +695,8 @@
1.595
1.596 Barrier is a set \e B of nodes for which
1.597
1.598 - \f[ \sum_{a\in\delta_{out}(B)} upper(a) -
1.599 - \sum_{a\in\delta_{in}(B)} lower(a) < \sum_{v\in B}delta(v) \f]
1.600 + \f[ \sum_{uv\in A: u\in B} upper(uv) -
1.601 + \sum_{uv\in A: v\in B} lower(uv) < \sum_{v\in B} sup(v) \f]
1.602
1.603 holds. The existence of a set with this property prooves that a
1.604 feasible circualtion cannot exist.
1.605 @@ -684,7 +726,7 @@
1.606 /// empty set, so \c bar[v] will be \c false for all nodes \c v.
1.607 ///
1.608 /// \note This function calls \ref barrier() for each node,
1.609 - /// so it runs in \f$O(n)\f$ time.
1.610 + /// so it runs in O(n) time.
1.611 ///
1.612 /// \pre Either \ref run() or \ref init() must be called before
1.613 /// using this function.
1.614 @@ -717,7 +759,7 @@
1.615 if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
1.616 for(NodeIt n(_g);n!=INVALID;++n)
1.617 {
1.618 - Value dif=-(*_delta)[n];
1.619 + Value dif=-(*_supply)[n];
1.620 for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
1.621 for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
1.622 if(_tol.negative(dif)) return false;
1.623 @@ -733,14 +775,20 @@
1.624 bool checkBarrier() const
1.625 {
1.626 Value delta=0;
1.627 + Value inf_cap = std::numeric_limits<Value>::has_infinity ?
1.628 + std::numeric_limits<Value>::infinity() :
1.629 + std::numeric_limits<Value>::max();
1.630 for(NodeIt n(_g);n!=INVALID;++n)
1.631 if(barrier(n))
1.632 - delta-=(*_delta)[n];
1.633 + delta-=(*_supply)[n];
1.634 for(ArcIt e(_g);e!=INVALID;++e)
1.635 {
1.636 Node s=_g.source(e);
1.637 Node t=_g.target(e);
1.638 - if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
1.639 + if(barrier(s)&&!barrier(t)) {
1.640 + if (_tol.less(inf_cap - (*_up)[e], delta)) return false;
1.641 + delta+=(*_up)[e];
1.642 + }
1.643 else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
1.644 }
1.645 return _tol.negative(delta);