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