0
4
0
1 | 1 |
/* -*- mode: C++; indent-tabs-mode: nil; -*- |
2 | 2 |
* |
3 | 3 |
* This file is a part of LEMON, a generic C++ optimization library. |
4 | 4 |
* |
5 | 5 |
* Copyright (C) 2003-2009 |
6 | 6 |
* Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport |
7 | 7 |
* (Egervary Research Group on Combinatorial Optimization, EGRES). |
8 | 8 |
* |
9 | 9 |
* Permission to use, modify and distribute this software is granted |
10 | 10 |
* provided that this copyright notice appears in all copies. For |
11 | 11 |
* precise terms see the accompanying LICENSE file. |
12 | 12 |
* |
13 | 13 |
* This software is provided "AS IS" with no warranty of any kind, |
14 | 14 |
* express or implied, and with no claim as to its suitability for any |
15 | 15 |
* purpose. |
16 | 16 |
* |
17 | 17 |
*/ |
18 | 18 |
|
19 | 19 |
#ifndef LEMON_CIRCULATION_H |
20 | 20 |
#define LEMON_CIRCULATION_H |
21 | 21 |
|
22 | 22 |
#include <lemon/tolerance.h> |
23 | 23 |
#include <lemon/elevator.h> |
24 | 24 |
#include <limits> |
25 | 25 |
|
26 | 26 |
///\ingroup max_flow |
27 | 27 |
///\file |
28 | 28 |
///\brief Push-relabel algorithm for finding a feasible circulation. |
29 | 29 |
/// |
30 | 30 |
namespace lemon { |
31 | 31 |
|
32 | 32 |
/// \brief Default traits class of Circulation class. |
33 | 33 |
/// |
34 | 34 |
/// Default traits class of Circulation class. |
35 | 35 |
/// |
36 | 36 |
/// \tparam GR Type of the digraph the algorithm runs on. |
37 | 37 |
/// \tparam LM The type of the lower bound map. |
38 | 38 |
/// \tparam UM The type of the upper bound (capacity) map. |
39 | 39 |
/// \tparam SM The type of the supply map. |
40 | 40 |
template <typename GR, typename LM, |
41 | 41 |
typename UM, typename SM> |
42 | 42 |
struct CirculationDefaultTraits { |
43 | 43 |
|
44 | 44 |
/// \brief The type of the digraph the algorithm runs on. |
45 | 45 |
typedef GR Digraph; |
46 | 46 |
|
47 | 47 |
/// \brief The type of the lower bound map. |
48 | 48 |
/// |
49 | 49 |
/// The type of the map that stores the lower bounds on the arcs. |
50 | 50 |
/// It must conform to the \ref concepts::ReadMap "ReadMap" concept. |
51 | 51 |
typedef LM LowerMap; |
52 | 52 |
|
53 | 53 |
/// \brief The type of the upper bound (capacity) map. |
54 | 54 |
/// |
55 | 55 |
/// The type of the map that stores the upper bounds (capacities) |
56 | 56 |
/// on the arcs. |
57 | 57 |
/// It must conform to the \ref concepts::ReadMap "ReadMap" concept. |
58 | 58 |
typedef UM UpperMap; |
59 | 59 |
|
60 | 60 |
/// \brief The type of supply map. |
61 | 61 |
/// |
62 | 62 |
/// The type of the map that stores the signed supply values of the |
63 | 63 |
/// nodes. |
64 | 64 |
/// It must conform to the \ref concepts::ReadMap "ReadMap" concept. |
65 | 65 |
typedef SM SupplyMap; |
66 | 66 |
|
67 | 67 |
/// \brief The type of the flow and supply values. |
68 | 68 |
typedef typename SupplyMap::Value Value; |
69 | 69 |
|
70 | 70 |
/// \brief The type of the map that stores the flow values. |
71 | 71 |
/// |
72 | 72 |
/// The type of the map that stores the flow values. |
73 | 73 |
/// It must conform to the \ref concepts::ReadWriteMap "ReadWriteMap" |
74 | 74 |
/// concept. |
75 | 75 |
typedef typename Digraph::template ArcMap<Value> FlowMap; |
76 | 76 |
|
77 | 77 |
/// \brief Instantiates a FlowMap. |
78 | 78 |
/// |
79 | 79 |
/// This function instantiates a \ref FlowMap. |
80 | 80 |
/// \param digraph The digraph for which we would like to define |
81 | 81 |
/// the flow map. |
82 | 82 |
static FlowMap* createFlowMap(const Digraph& digraph) { |
83 | 83 |
return new FlowMap(digraph); |
84 | 84 |
} |
85 | 85 |
|
86 | 86 |
/// \brief The elevator type used by the algorithm. |
87 | 87 |
/// |
88 | 88 |
/// The elevator type used by the algorithm. |
89 | 89 |
/// |
90 | 90 |
/// \sa Elevator |
91 | 91 |
/// \sa LinkedElevator |
92 | 92 |
typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator; |
93 | 93 |
|
94 | 94 |
/// \brief Instantiates an Elevator. |
95 | 95 |
/// |
96 | 96 |
/// This function instantiates an \ref Elevator. |
97 | 97 |
/// \param digraph The digraph for which we would like to define |
98 | 98 |
/// the elevator. |
99 | 99 |
/// \param max_level The maximum level of the elevator. |
100 | 100 |
static Elevator* createElevator(const Digraph& digraph, int max_level) { |
101 | 101 |
return new Elevator(digraph, max_level); |
102 | 102 |
} |
103 | 103 |
|
104 | 104 |
/// \brief The tolerance used by the algorithm |
105 | 105 |
/// |
106 | 106 |
/// The tolerance used by the algorithm to handle inexact computation. |
107 | 107 |
typedef lemon::Tolerance<Value> Tolerance; |
108 | 108 |
|
109 | 109 |
}; |
110 | 110 |
|
111 | 111 |
/** |
112 | 112 |
\brief Push-relabel algorithm for the network circulation problem. |
113 | 113 |
|
114 | 114 |
\ingroup max_flow |
115 | 115 |
This class implements a push-relabel algorithm for the \e network |
116 | 116 |
\e circulation problem. |
117 | 117 |
It is to find a feasible circulation when lower and upper bounds |
118 | 118 |
are given for the flow values on the arcs and lower bounds are |
119 | 119 |
given for the difference between the outgoing and incoming flow |
120 | 120 |
at the nodes. |
121 | 121 |
|
122 | 122 |
The exact formulation of this problem is the following. |
123 | 123 |
Let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{R}\f$ |
124 | 124 |
\f$upper: A\rightarrow\mathbf{R}\cup\{\infty\}\f$ denote the lower and |
125 | 125 |
upper bounds on the arcs, for which \f$lower(uv) \leq upper(uv)\f$ |
126 | 126 |
holds for all \f$uv\in A\f$, and \f$sup: V\rightarrow\mathbf{R}\f$ |
127 | 127 |
denotes the signed supply values of the nodes. |
128 | 128 |
If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$ |
129 | 129 |
supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with |
130 | 130 |
\f$-sup(u)\f$ demand. |
131 | 131 |
A feasible circulation is an \f$f: A\rightarrow\mathbf{R}\f$ |
132 | 132 |
solution of the following problem. |
133 | 133 |
|
134 | 134 |
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) |
135 | 135 |
\geq sup(u) \quad \forall u\in V, \f] |
136 | 136 |
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A. \f] |
137 | 137 |
|
138 | 138 |
The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be |
139 | 139 |
zero or negative in order to have a feasible solution (since the sum |
140 | 140 |
of the expressions on the left-hand side of the inequalities is zero). |
141 | 141 |
It means that the total demand must be greater or equal to the total |
142 | 142 |
supply and all the supplies have to be carried out from the supply nodes, |
143 | 143 |
but there could be demands that are not satisfied. |
144 | 144 |
If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand |
145 | 145 |
constraints have to be satisfied with equality, i.e. all demands |
146 | 146 |
have to be satisfied and all supplies have to be used. |
147 | 147 |
|
148 | 148 |
If you need the opposite inequalities in the supply/demand constraints |
149 | 149 |
(i.e. the total demand is less than the total supply and all the demands |
150 | 150 |
have to be satisfied while there could be supplies that are not used), |
151 | 151 |
then you could easily transform the problem to the above form by reversing |
152 | 152 |
the direction of the arcs and taking the negative of the supply values |
153 | 153 |
(e.g. using \ref ReverseDigraph and \ref NegMap adaptors). |
154 | 154 |
|
155 | 155 |
This algorithm either calculates a feasible circulation, or provides |
156 | 156 |
a \ref barrier() "barrier", which prooves that a feasible soultion |
157 | 157 |
cannot exist. |
158 | 158 |
|
159 | 159 |
Note that this algorithm also provides a feasible solution for the |
160 | 160 |
\ref min_cost_flow "minimum cost flow problem". |
161 | 161 |
|
162 | 162 |
\tparam GR The type of the digraph the algorithm runs on. |
163 | 163 |
\tparam LM The type of the lower bound map. The default |
164 | 164 |
map type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>". |
165 | 165 |
\tparam UM The type of the upper bound (capacity) map. |
166 | 166 |
The default map type is \c LM. |
167 | 167 |
\tparam SM The type of the supply map. The default map type is |
168 | 168 |
\ref concepts::Digraph::NodeMap "GR::NodeMap<UM::Value>". |
169 | 169 |
*/ |
170 | 170 |
#ifdef DOXYGEN |
171 | 171 |
template< typename GR, |
172 | 172 |
typename LM, |
173 | 173 |
typename UM, |
174 | 174 |
typename SM, |
175 | 175 |
typename TR > |
176 | 176 |
#else |
177 | 177 |
template< typename GR, |
178 | 178 |
typename LM = typename GR::template ArcMap<int>, |
179 | 179 |
typename UM = LM, |
180 | 180 |
typename SM = typename GR::template NodeMap<typename UM::Value>, |
181 | 181 |
typename TR = CirculationDefaultTraits<GR, LM, UM, SM> > |
182 | 182 |
#endif |
183 | 183 |
class Circulation { |
184 | 184 |
public: |
185 | 185 |
|
186 | 186 |
///The \ref CirculationDefaultTraits "traits class" of the algorithm. |
187 | 187 |
typedef TR Traits; |
188 | 188 |
///The type of the digraph the algorithm runs on. |
189 | 189 |
typedef typename Traits::Digraph Digraph; |
190 | 190 |
///The type of the flow and supply values. |
191 | 191 |
typedef typename Traits::Value Value; |
192 | 192 |
|
193 | 193 |
///The type of the lower bound map. |
194 | 194 |
typedef typename Traits::LowerMap LowerMap; |
195 | 195 |
///The type of the upper bound (capacity) map. |
196 | 196 |
typedef typename Traits::UpperMap UpperMap; |
197 | 197 |
///The type of the supply map. |
198 | 198 |
typedef typename Traits::SupplyMap SupplyMap; |
199 | 199 |
///The type of the flow map. |
200 | 200 |
typedef typename Traits::FlowMap FlowMap; |
201 | 201 |
|
202 | 202 |
///The type of the elevator. |
203 | 203 |
typedef typename Traits::Elevator Elevator; |
204 | 204 |
///The type of the tolerance. |
205 | 205 |
typedef typename Traits::Tolerance Tolerance; |
206 | 206 |
|
207 | 207 |
private: |
208 | 208 |
|
209 | 209 |
TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); |
210 | 210 |
|
211 | 211 |
const Digraph &_g; |
212 | 212 |
int _node_num; |
213 | 213 |
|
214 | 214 |
const LowerMap *_lo; |
215 | 215 |
const UpperMap *_up; |
216 | 216 |
const SupplyMap *_supply; |
217 | 217 |
|
218 | 218 |
FlowMap *_flow; |
219 | 219 |
bool _local_flow; |
220 | 220 |
|
221 | 221 |
Elevator* _level; |
222 | 222 |
bool _local_level; |
223 | 223 |
|
224 | 224 |
typedef typename Digraph::template NodeMap<Value> ExcessMap; |
225 | 225 |
ExcessMap* _excess; |
226 | 226 |
|
227 | 227 |
Tolerance _tol; |
228 | 228 |
int _el; |
229 | 229 |
|
230 | 230 |
public: |
231 | 231 |
|
232 | 232 |
typedef Circulation Create; |
233 | 233 |
|
234 | 234 |
///\name Named Template Parameters |
235 | 235 |
|
236 | 236 |
///@{ |
237 | 237 |
|
238 | 238 |
template <typename T> |
239 | 239 |
struct SetFlowMapTraits : public Traits { |
240 | 240 |
typedef T FlowMap; |
241 | 241 |
static FlowMap *createFlowMap(const Digraph&) { |
242 | 242 |
LEMON_ASSERT(false, "FlowMap is not initialized"); |
243 | 243 |
return 0; // ignore warnings |
244 | 244 |
} |
245 | 245 |
}; |
246 | 246 |
|
247 | 247 |
/// \brief \ref named-templ-param "Named parameter" for setting |
248 | 248 |
/// FlowMap type |
249 | 249 |
/// |
250 | 250 |
/// \ref named-templ-param "Named parameter" for setting FlowMap |
251 | 251 |
/// type. |
252 | 252 |
template <typename T> |
253 | 253 |
struct SetFlowMap |
254 | 254 |
: public Circulation<Digraph, LowerMap, UpperMap, SupplyMap, |
255 | 255 |
SetFlowMapTraits<T> > { |
256 | 256 |
typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap, |
257 | 257 |
SetFlowMapTraits<T> > Create; |
258 | 258 |
}; |
259 | 259 |
|
260 | 260 |
template <typename T> |
261 | 261 |
struct SetElevatorTraits : public Traits { |
262 | 262 |
typedef T Elevator; |
263 | 263 |
static Elevator *createElevator(const Digraph&, int) { |
264 | 264 |
LEMON_ASSERT(false, "Elevator is not initialized"); |
265 | 265 |
return 0; // ignore warnings |
266 | 266 |
} |
267 | 267 |
}; |
268 | 268 |
|
269 | 269 |
/// \brief \ref named-templ-param "Named parameter" for setting |
270 | 270 |
/// Elevator type |
271 | 271 |
/// |
272 | 272 |
/// \ref named-templ-param "Named parameter" for setting Elevator |
273 | 273 |
/// type. If this named parameter is used, then an external |
274 | 274 |
/// elevator object must be passed to the algorithm using the |
275 | 275 |
/// \ref elevator(Elevator&) "elevator()" function before calling |
276 | 276 |
/// \ref run() or \ref init(). |
277 | 277 |
/// \sa SetStandardElevator |
278 | 278 |
template <typename T> |
279 | 279 |
struct SetElevator |
280 | 280 |
: public Circulation<Digraph, LowerMap, UpperMap, SupplyMap, |
281 | 281 |
SetElevatorTraits<T> > { |
282 | 282 |
typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap, |
283 | 283 |
SetElevatorTraits<T> > Create; |
284 | 284 |
}; |
285 | 285 |
|
286 | 286 |
template <typename T> |
287 | 287 |
struct SetStandardElevatorTraits : public Traits { |
288 | 288 |
typedef T Elevator; |
289 | 289 |
static Elevator *createElevator(const Digraph& digraph, int max_level) { |
290 | 290 |
return new Elevator(digraph, max_level); |
291 | 291 |
} |
292 | 292 |
}; |
293 | 293 |
|
294 | 294 |
/// \brief \ref named-templ-param "Named parameter" for setting |
295 | 295 |
/// Elevator type with automatic allocation |
296 | 296 |
/// |
297 | 297 |
/// \ref named-templ-param "Named parameter" for setting Elevator |
298 | 298 |
/// type with automatic allocation. |
299 | 299 |
/// The Elevator should have standard constructor interface to be |
300 | 300 |
/// able to automatically created by the algorithm (i.e. the |
301 | 301 |
/// digraph and the maximum level should be passed to it). |
302 | 302 |
/// However an external elevator object could also be passed to the |
303 | 303 |
/// algorithm with the \ref elevator(Elevator&) "elevator()" function |
304 | 304 |
/// before calling \ref run() or \ref init(). |
305 | 305 |
/// \sa SetElevator |
306 | 306 |
template <typename T> |
307 | 307 |
struct SetStandardElevator |
308 | 308 |
: public Circulation<Digraph, LowerMap, UpperMap, SupplyMap, |
309 | 309 |
SetStandardElevatorTraits<T> > { |
310 | 310 |
typedef Circulation<Digraph, LowerMap, UpperMap, SupplyMap, |
311 | 311 |
SetStandardElevatorTraits<T> > Create; |
312 | 312 |
}; |
313 | 313 |
|
314 | 314 |
/// @} |
315 | 315 |
|
316 | 316 |
protected: |
317 | 317 |
|
318 | 318 |
Circulation() {} |
319 | 319 |
|
320 | 320 |
public: |
321 | 321 |
|
322 | 322 |
/// Constructor. |
323 | 323 |
|
324 | 324 |
/// The constructor of the class. |
325 | 325 |
/// |
326 | 326 |
/// \param graph The digraph the algorithm runs on. |
327 | 327 |
/// \param lower The lower bounds for the flow values on the arcs. |
328 | 328 |
/// \param upper The upper bounds (capacities) for the flow values |
329 | 329 |
/// on the arcs. |
330 | 330 |
/// \param supply The signed supply values of the nodes. |
331 | 331 |
Circulation(const Digraph &graph, const LowerMap &lower, |
332 | 332 |
const UpperMap &upper, const SupplyMap &supply) |
333 | 333 |
: _g(graph), _lo(&lower), _up(&upper), _supply(&supply), |
334 | 334 |
_flow(NULL), _local_flow(false), _level(NULL), _local_level(false), |
335 | 335 |
_excess(NULL) {} |
336 | 336 |
|
337 | 337 |
/// Destructor. |
338 | 338 |
~Circulation() { |
339 | 339 |
destroyStructures(); |
340 | 340 |
} |
341 | 341 |
|
342 | 342 |
|
343 | 343 |
private: |
344 | 344 |
|
345 | 345 |
bool checkBoundMaps() { |
346 | 346 |
for (ArcIt e(_g);e!=INVALID;++e) { |
347 | 347 |
if (_tol.less((*_up)[e], (*_lo)[e])) return false; |
348 | 348 |
} |
349 | 349 |
return true; |
350 | 350 |
} |
351 | 351 |
|
352 | 352 |
void createStructures() { |
353 | 353 |
_node_num = _el = countNodes(_g); |
354 | 354 |
|
355 | 355 |
if (!_flow) { |
356 | 356 |
_flow = Traits::createFlowMap(_g); |
357 | 357 |
_local_flow = true; |
358 | 358 |
} |
359 | 359 |
if (!_level) { |
360 | 360 |
_level = Traits::createElevator(_g, _node_num); |
361 | 361 |
_local_level = true; |
362 | 362 |
} |
363 | 363 |
if (!_excess) { |
364 | 364 |
_excess = new ExcessMap(_g); |
365 | 365 |
} |
366 | 366 |
} |
367 | 367 |
|
368 | 368 |
void destroyStructures() { |
369 | 369 |
if (_local_flow) { |
370 | 370 |
delete _flow; |
371 | 371 |
} |
372 | 372 |
if (_local_level) { |
373 | 373 |
delete _level; |
374 | 374 |
} |
375 | 375 |
if (_excess) { |
376 | 376 |
delete _excess; |
377 | 377 |
} |
378 | 378 |
} |
379 | 379 |
|
380 | 380 |
public: |
381 | 381 |
|
382 | 382 |
/// Sets the lower bound map. |
383 | 383 |
|
384 | 384 |
/// Sets the lower bound map. |
385 | 385 |
/// \return <tt>(*this)</tt> |
386 | 386 |
Circulation& lowerMap(const LowerMap& map) { |
387 | 387 |
_lo = ↦ |
388 | 388 |
return *this; |
389 | 389 |
} |
390 | 390 |
|
391 | 391 |
/// Sets the upper bound (capacity) map. |
392 | 392 |
|
393 | 393 |
/// Sets the upper bound (capacity) map. |
394 | 394 |
/// \return <tt>(*this)</tt> |
395 | 395 |
Circulation& upperMap(const UpperMap& map) { |
396 | 396 |
_up = ↦ |
397 | 397 |
return *this; |
398 | 398 |
} |
399 | 399 |
|
400 | 400 |
/// Sets the supply map. |
401 | 401 |
|
402 | 402 |
/// Sets the supply map. |
403 | 403 |
/// \return <tt>(*this)</tt> |
404 | 404 |
Circulation& supplyMap(const SupplyMap& map) { |
405 | 405 |
_supply = ↦ |
406 | 406 |
return *this; |
407 | 407 |
} |
408 | 408 |
|
409 | 409 |
/// \brief Sets the flow map. |
410 | 410 |
/// |
411 | 411 |
/// Sets the flow map. |
412 | 412 |
/// If you don't use this function before calling \ref run() or |
413 | 413 |
/// \ref init(), an instance will be allocated automatically. |
414 | 414 |
/// The destructor deallocates this automatically allocated map, |
415 | 415 |
/// of course. |
416 | 416 |
/// \return <tt>(*this)</tt> |
417 | 417 |
Circulation& flowMap(FlowMap& map) { |
418 | 418 |
if (_local_flow) { |
419 | 419 |
delete _flow; |
420 | 420 |
_local_flow = false; |
421 | 421 |
} |
422 | 422 |
_flow = ↦ |
423 | 423 |
return *this; |
424 | 424 |
} |
425 | 425 |
|
426 | 426 |
/// \brief Sets the elevator used by algorithm. |
427 | 427 |
/// |
428 | 428 |
/// Sets the elevator used by algorithm. |
429 | 429 |
/// If you don't use this function before calling \ref run() or |
430 | 430 |
/// \ref init(), an instance will be allocated automatically. |
431 | 431 |
/// The destructor deallocates this automatically allocated elevator, |
432 | 432 |
/// of course. |
433 | 433 |
/// \return <tt>(*this)</tt> |
434 | 434 |
Circulation& elevator(Elevator& elevator) { |
435 | 435 |
if (_local_level) { |
436 | 436 |
delete _level; |
437 | 437 |
_local_level = false; |
438 | 438 |
} |
439 | 439 |
_level = &elevator; |
440 | 440 |
return *this; |
441 | 441 |
} |
442 | 442 |
|
443 | 443 |
/// \brief Returns a const reference to the elevator. |
444 | 444 |
/// |
445 | 445 |
/// Returns a const reference to the elevator. |
446 | 446 |
/// |
447 | 447 |
/// \pre Either \ref run() or \ref init() must be called before |
448 | 448 |
/// using this function. |
449 | 449 |
const Elevator& elevator() const { |
450 | 450 |
return *_level; |
451 | 451 |
} |
452 | 452 |
|
453 |
/// \brief Sets the tolerance used by algorithm. |
|
453 |
/// \brief Sets the tolerance used by the algorithm. |
|
454 | 454 |
/// |
455 |
/// Sets the tolerance used by algorithm. |
|
455 |
/// Sets the tolerance object used by the algorithm. |
|
456 |
/// \return <tt>(*this)</tt> |
|
456 | 457 |
Circulation& tolerance(const Tolerance& tolerance) { |
457 | 458 |
_tol = tolerance; |
458 | 459 |
return *this; |
459 | 460 |
} |
460 | 461 |
|
461 | 462 |
/// \brief Returns a const reference to the tolerance. |
462 | 463 |
/// |
463 |
/// Returns a const reference to the tolerance |
|
464 |
/// Returns a const reference to the tolerance object used by |
|
465 |
/// the algorithm. |
|
464 | 466 |
const Tolerance& tolerance() const { |
465 | 467 |
return _tol; |
466 | 468 |
} |
467 | 469 |
|
468 | 470 |
/// \name Execution Control |
469 | 471 |
/// The simplest way to execute the algorithm is to call \ref run().\n |
470 | 472 |
/// If you need more control on the initial solution or the execution, |
471 | 473 |
/// first you have to call one of the \ref init() functions, then |
472 | 474 |
/// the \ref start() function. |
473 | 475 |
|
474 | 476 |
///@{ |
475 | 477 |
|
476 | 478 |
/// Initializes the internal data structures. |
477 | 479 |
|
478 | 480 |
/// Initializes the internal data structures and sets all flow values |
479 | 481 |
/// to the lower bound. |
480 | 482 |
void init() |
481 | 483 |
{ |
482 | 484 |
LEMON_DEBUG(checkBoundMaps(), |
483 | 485 |
"Upper bounds must be greater or equal to the lower bounds"); |
484 | 486 |
|
485 | 487 |
createStructures(); |
486 | 488 |
|
487 | 489 |
for(NodeIt n(_g);n!=INVALID;++n) { |
488 | 490 |
(*_excess)[n] = (*_supply)[n]; |
489 | 491 |
} |
490 | 492 |
|
491 | 493 |
for (ArcIt e(_g);e!=INVALID;++e) { |
492 | 494 |
_flow->set(e, (*_lo)[e]); |
493 | 495 |
(*_excess)[_g.target(e)] += (*_flow)[e]; |
494 | 496 |
(*_excess)[_g.source(e)] -= (*_flow)[e]; |
495 | 497 |
} |
496 | 498 |
|
497 | 499 |
// global relabeling tested, but in general case it provides |
498 | 500 |
// worse performance for random digraphs |
499 | 501 |
_level->initStart(); |
500 | 502 |
for(NodeIt n(_g);n!=INVALID;++n) |
501 | 503 |
_level->initAddItem(n); |
502 | 504 |
_level->initFinish(); |
503 | 505 |
for(NodeIt n(_g);n!=INVALID;++n) |
504 | 506 |
if(_tol.positive((*_excess)[n])) |
505 | 507 |
_level->activate(n); |
506 | 508 |
} |
507 | 509 |
|
508 | 510 |
/// Initializes the internal data structures using a greedy approach. |
509 | 511 |
|
510 | 512 |
/// Initializes the internal data structures using a greedy approach |
511 | 513 |
/// to construct the initial solution. |
512 | 514 |
void greedyInit() |
513 | 515 |
{ |
514 | 516 |
LEMON_DEBUG(checkBoundMaps(), |
515 | 517 |
"Upper bounds must be greater or equal to the lower bounds"); |
516 | 518 |
|
517 | 519 |
createStructures(); |
518 | 520 |
|
519 | 521 |
for(NodeIt n(_g);n!=INVALID;++n) { |
520 | 522 |
(*_excess)[n] = (*_supply)[n]; |
521 | 523 |
} |
522 | 524 |
|
523 | 525 |
for (ArcIt e(_g);e!=INVALID;++e) { |
524 | 526 |
if (!_tol.less(-(*_excess)[_g.target(e)], (*_up)[e])) { |
525 | 527 |
_flow->set(e, (*_up)[e]); |
526 | 528 |
(*_excess)[_g.target(e)] += (*_up)[e]; |
527 | 529 |
(*_excess)[_g.source(e)] -= (*_up)[e]; |
528 | 530 |
} else if (_tol.less(-(*_excess)[_g.target(e)], (*_lo)[e])) { |
529 | 531 |
_flow->set(e, (*_lo)[e]); |
530 | 532 |
(*_excess)[_g.target(e)] += (*_lo)[e]; |
531 | 533 |
(*_excess)[_g.source(e)] -= (*_lo)[e]; |
532 | 534 |
} else { |
533 | 535 |
Value fc = -(*_excess)[_g.target(e)]; |
534 | 536 |
_flow->set(e, fc); |
535 | 537 |
(*_excess)[_g.target(e)] = 0; |
536 | 538 |
(*_excess)[_g.source(e)] -= fc; |
537 | 539 |
} |
538 | 540 |
} |
539 | 541 |
|
540 | 542 |
_level->initStart(); |
541 | 543 |
for(NodeIt n(_g);n!=INVALID;++n) |
542 | 544 |
_level->initAddItem(n); |
543 | 545 |
_level->initFinish(); |
544 | 546 |
for(NodeIt n(_g);n!=INVALID;++n) |
545 | 547 |
if(_tol.positive((*_excess)[n])) |
546 | 548 |
_level->activate(n); |
547 | 549 |
} |
548 | 550 |
|
549 | 551 |
///Executes the algorithm |
550 | 552 |
|
551 | 553 |
///This function executes the algorithm. |
552 | 554 |
/// |
553 | 555 |
///\return \c true if a feasible circulation is found. |
554 | 556 |
/// |
555 | 557 |
///\sa barrier() |
556 | 558 |
///\sa barrierMap() |
557 | 559 |
bool start() |
558 | 560 |
{ |
559 | 561 |
|
560 | 562 |
Node act; |
561 | 563 |
Node bact=INVALID; |
562 | 564 |
Node last_activated=INVALID; |
563 | 565 |
while((act=_level->highestActive())!=INVALID) { |
564 | 566 |
int actlevel=(*_level)[act]; |
565 | 567 |
int mlevel=_node_num; |
566 | 568 |
Value exc=(*_excess)[act]; |
567 | 569 |
|
568 | 570 |
for(OutArcIt e(_g,act);e!=INVALID; ++e) { |
569 | 571 |
Node v = _g.target(e); |
570 | 572 |
Value fc=(*_up)[e]-(*_flow)[e]; |
571 | 573 |
if(!_tol.positive(fc)) continue; |
572 | 574 |
if((*_level)[v]<actlevel) { |
573 | 575 |
if(!_tol.less(fc, exc)) { |
574 | 576 |
_flow->set(e, (*_flow)[e] + exc); |
575 | 577 |
(*_excess)[v] += exc; |
576 | 578 |
if(!_level->active(v) && _tol.positive((*_excess)[v])) |
577 | 579 |
_level->activate(v); |
578 | 580 |
(*_excess)[act] = 0; |
579 | 581 |
_level->deactivate(act); |
580 | 582 |
goto next_l; |
581 | 583 |
} |
582 | 584 |
else { |
583 | 585 |
_flow->set(e, (*_up)[e]); |
584 | 586 |
(*_excess)[v] += fc; |
585 | 587 |
if(!_level->active(v) && _tol.positive((*_excess)[v])) |
586 | 588 |
_level->activate(v); |
587 | 589 |
exc-=fc; |
588 | 590 |
} |
589 | 591 |
} |
590 | 592 |
else if((*_level)[v]<mlevel) mlevel=(*_level)[v]; |
591 | 593 |
} |
592 | 594 |
for(InArcIt e(_g,act);e!=INVALID; ++e) { |
593 | 595 |
Node v = _g.source(e); |
594 | 596 |
Value fc=(*_flow)[e]-(*_lo)[e]; |
595 | 597 |
if(!_tol.positive(fc)) continue; |
596 | 598 |
if((*_level)[v]<actlevel) { |
597 | 599 |
if(!_tol.less(fc, exc)) { |
598 | 600 |
_flow->set(e, (*_flow)[e] - exc); |
599 | 601 |
(*_excess)[v] += exc; |
600 | 602 |
if(!_level->active(v) && _tol.positive((*_excess)[v])) |
601 | 603 |
_level->activate(v); |
602 | 604 |
(*_excess)[act] = 0; |
603 | 605 |
_level->deactivate(act); |
604 | 606 |
goto next_l; |
605 | 607 |
} |
606 | 608 |
else { |
607 | 609 |
_flow->set(e, (*_lo)[e]); |
608 | 610 |
(*_excess)[v] += fc; |
609 | 611 |
if(!_level->active(v) && _tol.positive((*_excess)[v])) |
610 | 612 |
_level->activate(v); |
611 | 613 |
exc-=fc; |
612 | 614 |
} |
613 | 615 |
} |
614 | 616 |
else if((*_level)[v]<mlevel) mlevel=(*_level)[v]; |
615 | 617 |
} |
616 | 618 |
|
617 | 619 |
(*_excess)[act] = exc; |
618 | 620 |
if(!_tol.positive(exc)) _level->deactivate(act); |
619 | 621 |
else if(mlevel==_node_num) { |
620 | 622 |
_level->liftHighestActiveToTop(); |
621 | 623 |
_el = _node_num; |
622 | 624 |
return false; |
623 | 625 |
} |
624 | 626 |
else { |
625 | 627 |
_level->liftHighestActive(mlevel+1); |
626 | 628 |
if(_level->onLevel(actlevel)==0) { |
627 | 629 |
_el = actlevel; |
628 | 630 |
return false; |
629 | 631 |
} |
630 | 632 |
} |
631 | 633 |
next_l: |
632 | 634 |
; |
633 | 635 |
} |
634 | 636 |
return true; |
635 | 637 |
} |
636 | 638 |
|
637 | 639 |
/// Runs the algorithm. |
638 | 640 |
|
639 | 641 |
/// This function runs the algorithm. |
640 | 642 |
/// |
641 | 643 |
/// \return \c true if a feasible circulation is found. |
642 | 644 |
/// |
643 | 645 |
/// \note Apart from the return value, c.run() is just a shortcut of |
644 | 646 |
/// the following code. |
645 | 647 |
/// \code |
646 | 648 |
/// c.greedyInit(); |
647 | 649 |
/// c.start(); |
648 | 650 |
/// \endcode |
649 | 651 |
bool run() { |
650 | 652 |
greedyInit(); |
651 | 653 |
return start(); |
652 | 654 |
} |
653 | 655 |
|
654 | 656 |
/// @} |
655 | 657 |
|
656 | 658 |
/// \name Query Functions |
657 | 659 |
/// The results of the circulation algorithm can be obtained using |
658 | 660 |
/// these functions.\n |
659 | 661 |
/// Either \ref run() or \ref start() should be called before |
660 | 662 |
/// using them. |
661 | 663 |
|
662 | 664 |
///@{ |
663 | 665 |
|
664 | 666 |
/// \brief Returns the flow value on the given arc. |
665 | 667 |
/// |
666 | 668 |
/// Returns the flow value on the given arc. |
667 | 669 |
/// |
668 | 670 |
/// \pre Either \ref run() or \ref init() must be called before |
669 | 671 |
/// using this function. |
670 | 672 |
Value flow(const Arc& arc) const { |
671 | 673 |
return (*_flow)[arc]; |
672 | 674 |
} |
673 | 675 |
|
674 | 676 |
/// \brief Returns a const reference to the flow map. |
675 | 677 |
/// |
676 | 678 |
/// Returns a const reference to the arc map storing the found flow. |
677 | 679 |
/// |
678 | 680 |
/// \pre Either \ref run() or \ref init() must be called before |
679 | 681 |
/// using this function. |
680 | 682 |
const FlowMap& flowMap() const { |
681 | 683 |
return *_flow; |
682 | 684 |
} |
683 | 685 |
|
684 | 686 |
/** |
685 | 687 |
\brief Returns \c true if the given node is in a barrier. |
686 | 688 |
|
687 | 689 |
Barrier is a set \e B of nodes for which |
688 | 690 |
|
689 | 691 |
\f[ \sum_{uv\in A: u\in B} upper(uv) - |
690 | 692 |
\sum_{uv\in A: v\in B} lower(uv) < \sum_{v\in B} sup(v) \f] |
691 | 693 |
|
692 | 694 |
holds. The existence of a set with this property prooves that a |
693 | 695 |
feasible circualtion cannot exist. |
694 | 696 |
|
695 | 697 |
This function returns \c true if the given node is in the found |
696 | 698 |
barrier. If a feasible circulation is found, the function |
697 | 699 |
gives back \c false for every node. |
698 | 700 |
|
699 | 701 |
\pre Either \ref run() or \ref init() must be called before |
700 | 702 |
using this function. |
701 | 703 |
|
702 | 704 |
\sa barrierMap() |
703 | 705 |
\sa checkBarrier() |
704 | 706 |
*/ |
705 | 707 |
bool barrier(const Node& node) const |
706 | 708 |
{ |
707 | 709 |
return (*_level)[node] >= _el; |
708 | 710 |
} |
709 | 711 |
|
710 | 712 |
/// \brief Gives back a barrier. |
711 | 713 |
/// |
712 | 714 |
/// This function sets \c bar to the characteristic vector of the |
713 | 715 |
/// found barrier. \c bar should be a \ref concepts::WriteMap "writable" |
714 | 716 |
/// node map with \c bool (or convertible) value type. |
715 | 717 |
/// |
716 | 718 |
/// If a feasible circulation is found, the function gives back an |
717 | 719 |
/// empty set, so \c bar[v] will be \c false for all nodes \c v. |
718 | 720 |
/// |
719 | 721 |
/// \note This function calls \ref barrier() for each node, |
720 | 722 |
/// so it runs in O(n) time. |
721 | 723 |
/// |
722 | 724 |
/// \pre Either \ref run() or \ref init() must be called before |
723 | 725 |
/// using this function. |
724 | 726 |
/// |
725 | 727 |
/// \sa barrier() |
726 | 728 |
/// \sa checkBarrier() |
727 | 729 |
template<class BarrierMap> |
728 | 730 |
void barrierMap(BarrierMap &bar) const |
729 | 731 |
{ |
730 | 732 |
for(NodeIt n(_g);n!=INVALID;++n) |
731 | 733 |
bar.set(n, (*_level)[n] >= _el); |
732 | 734 |
} |
733 | 735 |
|
734 | 736 |
/// @} |
735 | 737 |
|
736 | 738 |
/// \name Checker Functions |
737 | 739 |
/// The feasibility of the results can be checked using |
738 | 740 |
/// these functions.\n |
739 | 741 |
/// Either \ref run() or \ref start() should be called before |
740 | 742 |
/// using them. |
741 | 743 |
|
742 | 744 |
///@{ |
743 | 745 |
|
744 | 746 |
///Check if the found flow is a feasible circulation |
745 | 747 |
|
746 | 748 |
///Check if the found flow is a feasible circulation, |
747 | 749 |
/// |
748 | 750 |
bool checkFlow() const { |
749 | 751 |
for(ArcIt e(_g);e!=INVALID;++e) |
750 | 752 |
if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false; |
751 | 753 |
for(NodeIt n(_g);n!=INVALID;++n) |
752 | 754 |
{ |
753 | 755 |
Value dif=-(*_supply)[n]; |
754 | 756 |
for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e]; |
755 | 757 |
for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e]; |
756 | 758 |
if(_tol.negative(dif)) return false; |
757 | 759 |
} |
758 | 760 |
return true; |
759 | 761 |
} |
760 | 762 |
|
761 | 763 |
///Check whether or not the last execution provides a barrier |
762 | 764 |
|
763 | 765 |
///Check whether or not the last execution provides a barrier. |
764 | 766 |
///\sa barrier() |
765 | 767 |
///\sa barrierMap() |
766 | 768 |
bool checkBarrier() const |
767 | 769 |
{ |
768 | 770 |
Value delta=0; |
769 | 771 |
Value inf_cap = std::numeric_limits<Value>::has_infinity ? |
770 | 772 |
std::numeric_limits<Value>::infinity() : |
771 | 773 |
std::numeric_limits<Value>::max(); |
772 | 774 |
for(NodeIt n(_g);n!=INVALID;++n) |
773 | 775 |
if(barrier(n)) |
774 | 776 |
delta-=(*_supply)[n]; |
775 | 777 |
for(ArcIt e(_g);e!=INVALID;++e) |
776 | 778 |
{ |
777 | 779 |
Node s=_g.source(e); |
778 | 780 |
Node t=_g.target(e); |
779 | 781 |
if(barrier(s)&&!barrier(t)) { |
780 | 782 |
if (_tol.less(inf_cap - (*_up)[e], delta)) return false; |
781 | 783 |
delta+=(*_up)[e]; |
782 | 784 |
} |
783 | 785 |
else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e]; |
784 | 786 |
} |
785 | 787 |
return _tol.negative(delta); |
786 | 788 |
} |
787 | 789 |
|
788 | 790 |
/// @} |
789 | 791 |
|
790 | 792 |
}; |
791 | 793 |
|
792 | 794 |
} |
793 | 795 |
|
794 | 796 |
#endif |
1 | 1 |
/* -*- mode: C++; indent-tabs-mode: nil; -*- |
2 | 2 |
* |
3 | 3 |
* This file is a part of LEMON, a generic C++ optimization library. |
4 | 4 |
* |
5 | 5 |
* Copyright (C) 2003-2009 |
6 | 6 |
* Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport |
7 | 7 |
* (Egervary Research Group on Combinatorial Optimization, EGRES). |
8 | 8 |
* |
9 | 9 |
* Permission to use, modify and distribute this software is granted |
10 | 10 |
* provided that this copyright notice appears in all copies. For |
11 | 11 |
* precise terms see the accompanying LICENSE file. |
12 | 12 |
* |
13 | 13 |
* This software is provided "AS IS" with no warranty of any kind, |
14 | 14 |
* express or implied, and with no claim as to its suitability for any |
15 | 15 |
* purpose. |
16 | 16 |
* |
17 | 17 |
*/ |
18 | 18 |
|
19 | 19 |
#ifndef LEMON_PREFLOW_H |
20 | 20 |
#define LEMON_PREFLOW_H |
21 | 21 |
|
22 | 22 |
#include <lemon/tolerance.h> |
23 | 23 |
#include <lemon/elevator.h> |
24 | 24 |
|
25 | 25 |
/// \file |
26 | 26 |
/// \ingroup max_flow |
27 | 27 |
/// \brief Implementation of the preflow algorithm. |
28 | 28 |
|
29 | 29 |
namespace lemon { |
30 | 30 |
|
31 | 31 |
/// \brief Default traits class of Preflow class. |
32 | 32 |
/// |
33 | 33 |
/// Default traits class of Preflow class. |
34 | 34 |
/// \tparam GR Digraph type. |
35 | 35 |
/// \tparam CAP Capacity map type. |
36 | 36 |
template <typename GR, typename CAP> |
37 | 37 |
struct PreflowDefaultTraits { |
38 | 38 |
|
39 | 39 |
/// \brief The type of the digraph the algorithm runs on. |
40 | 40 |
typedef GR Digraph; |
41 | 41 |
|
42 | 42 |
/// \brief The type of the map that stores the arc capacities. |
43 | 43 |
/// |
44 | 44 |
/// The type of the map that stores the arc capacities. |
45 | 45 |
/// It must meet the \ref concepts::ReadMap "ReadMap" concept. |
46 | 46 |
typedef CAP CapacityMap; |
47 | 47 |
|
48 | 48 |
/// \brief The type of the flow values. |
49 | 49 |
typedef typename CapacityMap::Value Value; |
50 | 50 |
|
51 | 51 |
/// \brief The type of the map that stores the flow values. |
52 | 52 |
/// |
53 | 53 |
/// The type of the map that stores the flow values. |
54 | 54 |
/// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept. |
55 | 55 |
typedef typename Digraph::template ArcMap<Value> FlowMap; |
56 | 56 |
|
57 | 57 |
/// \brief Instantiates a FlowMap. |
58 | 58 |
/// |
59 | 59 |
/// This function instantiates a \ref FlowMap. |
60 | 60 |
/// \param digraph The digraph for which we would like to define |
61 | 61 |
/// the flow map. |
62 | 62 |
static FlowMap* createFlowMap(const Digraph& digraph) { |
63 | 63 |
return new FlowMap(digraph); |
64 | 64 |
} |
65 | 65 |
|
66 | 66 |
/// \brief The elevator type used by Preflow algorithm. |
67 | 67 |
/// |
68 | 68 |
/// The elevator type used by Preflow algorithm. |
69 | 69 |
/// |
70 | 70 |
/// \sa Elevator |
71 | 71 |
/// \sa LinkedElevator |
72 | 72 |
typedef LinkedElevator<Digraph, typename Digraph::Node> Elevator; |
73 | 73 |
|
74 | 74 |
/// \brief Instantiates an Elevator. |
75 | 75 |
/// |
76 | 76 |
/// This function instantiates an \ref Elevator. |
77 | 77 |
/// \param digraph The digraph for which we would like to define |
78 | 78 |
/// the elevator. |
79 | 79 |
/// \param max_level The maximum level of the elevator. |
80 | 80 |
static Elevator* createElevator(const Digraph& digraph, int max_level) { |
81 | 81 |
return new Elevator(digraph, max_level); |
82 | 82 |
} |
83 | 83 |
|
84 | 84 |
/// \brief The tolerance used by the algorithm |
85 | 85 |
/// |
86 | 86 |
/// The tolerance used by the algorithm to handle inexact computation. |
87 | 87 |
typedef lemon::Tolerance<Value> Tolerance; |
88 | 88 |
|
89 | 89 |
}; |
90 | 90 |
|
91 | 91 |
|
92 | 92 |
/// \ingroup max_flow |
93 | 93 |
/// |
94 | 94 |
/// \brief %Preflow algorithm class. |
95 | 95 |
/// |
96 | 96 |
/// This class provides an implementation of Goldberg-Tarjan's \e preflow |
97 | 97 |
/// \e push-relabel algorithm producing a \ref max_flow |
98 | 98 |
/// "flow of maximum value" in a digraph. |
99 | 99 |
/// The preflow algorithms are the fastest known maximum |
100 |
/// flow algorithms. The current implementation |
|
100 |
/// flow algorithms. The current implementation uses a mixture of the |
|
101 | 101 |
/// \e "highest label" and the \e "bound decrease" heuristics. |
102 | 102 |
/// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$. |
103 | 103 |
/// |
104 | 104 |
/// The algorithm consists of two phases. After the first phase |
105 | 105 |
/// the maximum flow value and the minimum cut is obtained. The |
106 | 106 |
/// second phase constructs a feasible maximum flow on each arc. |
107 | 107 |
/// |
108 | 108 |
/// \tparam GR The type of the digraph the algorithm runs on. |
109 | 109 |
/// \tparam CAP The type of the capacity map. The default map |
110 | 110 |
/// type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>". |
111 | 111 |
#ifdef DOXYGEN |
112 | 112 |
template <typename GR, typename CAP, typename TR> |
113 | 113 |
#else |
114 | 114 |
template <typename GR, |
115 | 115 |
typename CAP = typename GR::template ArcMap<int>, |
116 | 116 |
typename TR = PreflowDefaultTraits<GR, CAP> > |
117 | 117 |
#endif |
118 | 118 |
class Preflow { |
119 | 119 |
public: |
120 | 120 |
|
121 | 121 |
///The \ref PreflowDefaultTraits "traits class" of the algorithm. |
122 | 122 |
typedef TR Traits; |
123 | 123 |
///The type of the digraph the algorithm runs on. |
124 | 124 |
typedef typename Traits::Digraph Digraph; |
125 | 125 |
///The type of the capacity map. |
126 | 126 |
typedef typename Traits::CapacityMap CapacityMap; |
127 | 127 |
///The type of the flow values. |
128 | 128 |
typedef typename Traits::Value Value; |
129 | 129 |
|
130 | 130 |
///The type of the flow map. |
131 | 131 |
typedef typename Traits::FlowMap FlowMap; |
132 | 132 |
///The type of the elevator. |
133 | 133 |
typedef typename Traits::Elevator Elevator; |
134 | 134 |
///The type of the tolerance. |
135 | 135 |
typedef typename Traits::Tolerance Tolerance; |
136 | 136 |
|
137 | 137 |
private: |
138 | 138 |
|
139 | 139 |
TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); |
140 | 140 |
|
141 | 141 |
const Digraph& _graph; |
142 | 142 |
const CapacityMap* _capacity; |
143 | 143 |
|
144 | 144 |
int _node_num; |
145 | 145 |
|
146 | 146 |
Node _source, _target; |
147 | 147 |
|
148 | 148 |
FlowMap* _flow; |
149 | 149 |
bool _local_flow; |
150 | 150 |
|
151 | 151 |
Elevator* _level; |
152 | 152 |
bool _local_level; |
153 | 153 |
|
154 | 154 |
typedef typename Digraph::template NodeMap<Value> ExcessMap; |
155 | 155 |
ExcessMap* _excess; |
156 | 156 |
|
157 | 157 |
Tolerance _tolerance; |
158 | 158 |
|
159 | 159 |
bool _phase; |
160 | 160 |
|
161 | 161 |
|
162 | 162 |
void createStructures() { |
163 | 163 |
_node_num = countNodes(_graph); |
164 | 164 |
|
165 | 165 |
if (!_flow) { |
166 | 166 |
_flow = Traits::createFlowMap(_graph); |
167 | 167 |
_local_flow = true; |
168 | 168 |
} |
169 | 169 |
if (!_level) { |
170 | 170 |
_level = Traits::createElevator(_graph, _node_num); |
171 | 171 |
_local_level = true; |
172 | 172 |
} |
173 | 173 |
if (!_excess) { |
174 | 174 |
_excess = new ExcessMap(_graph); |
175 | 175 |
} |
176 | 176 |
} |
177 | 177 |
|
178 | 178 |
void destroyStructures() { |
179 | 179 |
if (_local_flow) { |
180 | 180 |
delete _flow; |
181 | 181 |
} |
182 | 182 |
if (_local_level) { |
183 | 183 |
delete _level; |
184 | 184 |
} |
185 | 185 |
if (_excess) { |
186 | 186 |
delete _excess; |
187 | 187 |
} |
188 | 188 |
} |
189 | 189 |
|
190 | 190 |
public: |
191 | 191 |
|
192 | 192 |
typedef Preflow Create; |
193 | 193 |
|
194 | 194 |
///\name Named Template Parameters |
195 | 195 |
|
196 | 196 |
///@{ |
197 | 197 |
|
198 | 198 |
template <typename T> |
199 | 199 |
struct SetFlowMapTraits : public Traits { |
200 | 200 |
typedef T FlowMap; |
201 | 201 |
static FlowMap *createFlowMap(const Digraph&) { |
202 | 202 |
LEMON_ASSERT(false, "FlowMap is not initialized"); |
203 | 203 |
return 0; // ignore warnings |
204 | 204 |
} |
205 | 205 |
}; |
206 | 206 |
|
207 | 207 |
/// \brief \ref named-templ-param "Named parameter" for setting |
208 | 208 |
/// FlowMap type |
209 | 209 |
/// |
210 | 210 |
/// \ref named-templ-param "Named parameter" for setting FlowMap |
211 | 211 |
/// type. |
212 | 212 |
template <typename T> |
213 | 213 |
struct SetFlowMap |
214 | 214 |
: public Preflow<Digraph, CapacityMap, SetFlowMapTraits<T> > { |
215 | 215 |
typedef Preflow<Digraph, CapacityMap, |
216 | 216 |
SetFlowMapTraits<T> > Create; |
217 | 217 |
}; |
218 | 218 |
|
219 | 219 |
template <typename T> |
220 | 220 |
struct SetElevatorTraits : public Traits { |
221 | 221 |
typedef T Elevator; |
222 | 222 |
static Elevator *createElevator(const Digraph&, int) { |
223 | 223 |
LEMON_ASSERT(false, "Elevator is not initialized"); |
224 | 224 |
return 0; // ignore warnings |
225 | 225 |
} |
226 | 226 |
}; |
227 | 227 |
|
228 | 228 |
/// \brief \ref named-templ-param "Named parameter" for setting |
229 | 229 |
/// Elevator type |
230 | 230 |
/// |
231 | 231 |
/// \ref named-templ-param "Named parameter" for setting Elevator |
232 | 232 |
/// type. If this named parameter is used, then an external |
233 | 233 |
/// elevator object must be passed to the algorithm using the |
234 | 234 |
/// \ref elevator(Elevator&) "elevator()" function before calling |
235 | 235 |
/// \ref run() or \ref init(). |
236 | 236 |
/// \sa SetStandardElevator |
237 | 237 |
template <typename T> |
238 | 238 |
struct SetElevator |
239 | 239 |
: public Preflow<Digraph, CapacityMap, SetElevatorTraits<T> > { |
240 | 240 |
typedef Preflow<Digraph, CapacityMap, |
241 | 241 |
SetElevatorTraits<T> > Create; |
242 | 242 |
}; |
243 | 243 |
|
244 | 244 |
template <typename T> |
245 | 245 |
struct SetStandardElevatorTraits : public Traits { |
246 | 246 |
typedef T Elevator; |
247 | 247 |
static Elevator *createElevator(const Digraph& digraph, int max_level) { |
248 | 248 |
return new Elevator(digraph, max_level); |
249 | 249 |
} |
250 | 250 |
}; |
251 | 251 |
|
252 | 252 |
/// \brief \ref named-templ-param "Named parameter" for setting |
253 | 253 |
/// Elevator type with automatic allocation |
254 | 254 |
/// |
255 | 255 |
/// \ref named-templ-param "Named parameter" for setting Elevator |
256 | 256 |
/// type with automatic allocation. |
257 | 257 |
/// The Elevator should have standard constructor interface to be |
258 | 258 |
/// able to automatically created by the algorithm (i.e. the |
259 | 259 |
/// digraph and the maximum level should be passed to it). |
260 | 260 |
/// However an external elevator object could also be passed to the |
261 | 261 |
/// algorithm with the \ref elevator(Elevator&) "elevator()" function |
262 | 262 |
/// before calling \ref run() or \ref init(). |
263 | 263 |
/// \sa SetElevator |
264 | 264 |
template <typename T> |
265 | 265 |
struct SetStandardElevator |
266 | 266 |
: public Preflow<Digraph, CapacityMap, |
267 | 267 |
SetStandardElevatorTraits<T> > { |
268 | 268 |
typedef Preflow<Digraph, CapacityMap, |
269 | 269 |
SetStandardElevatorTraits<T> > Create; |
270 | 270 |
}; |
271 | 271 |
|
272 | 272 |
/// @} |
273 | 273 |
|
274 | 274 |
protected: |
275 | 275 |
|
276 | 276 |
Preflow() {} |
277 | 277 |
|
278 | 278 |
public: |
279 | 279 |
|
280 | 280 |
|
281 | 281 |
/// \brief The constructor of the class. |
282 | 282 |
/// |
283 | 283 |
/// The constructor of the class. |
284 | 284 |
/// \param digraph The digraph the algorithm runs on. |
285 | 285 |
/// \param capacity The capacity of the arcs. |
286 | 286 |
/// \param source The source node. |
287 | 287 |
/// \param target The target node. |
288 | 288 |
Preflow(const Digraph& digraph, const CapacityMap& capacity, |
289 | 289 |
Node source, Node target) |
290 | 290 |
: _graph(digraph), _capacity(&capacity), |
291 | 291 |
_node_num(0), _source(source), _target(target), |
292 | 292 |
_flow(0), _local_flow(false), |
293 | 293 |
_level(0), _local_level(false), |
294 | 294 |
_excess(0), _tolerance(), _phase() {} |
295 | 295 |
|
296 | 296 |
/// \brief Destructor. |
297 | 297 |
/// |
298 | 298 |
/// Destructor. |
299 | 299 |
~Preflow() { |
300 | 300 |
destroyStructures(); |
301 | 301 |
} |
302 | 302 |
|
303 | 303 |
/// \brief Sets the capacity map. |
304 | 304 |
/// |
305 | 305 |
/// Sets the capacity map. |
306 | 306 |
/// \return <tt>(*this)</tt> |
307 | 307 |
Preflow& capacityMap(const CapacityMap& map) { |
308 | 308 |
_capacity = ↦ |
309 | 309 |
return *this; |
310 | 310 |
} |
311 | 311 |
|
312 | 312 |
/// \brief Sets the flow map. |
313 | 313 |
/// |
314 | 314 |
/// Sets the flow map. |
315 | 315 |
/// If you don't use this function before calling \ref run() or |
316 | 316 |
/// \ref init(), an instance will be allocated automatically. |
317 | 317 |
/// The destructor deallocates this automatically allocated map, |
318 | 318 |
/// of course. |
319 | 319 |
/// \return <tt>(*this)</tt> |
320 | 320 |
Preflow& flowMap(FlowMap& map) { |
321 | 321 |
if (_local_flow) { |
322 | 322 |
delete _flow; |
323 | 323 |
_local_flow = false; |
324 | 324 |
} |
325 | 325 |
_flow = ↦ |
326 | 326 |
return *this; |
327 | 327 |
} |
328 | 328 |
|
329 | 329 |
/// \brief Sets the source node. |
330 | 330 |
/// |
331 | 331 |
/// Sets the source node. |
332 | 332 |
/// \return <tt>(*this)</tt> |
333 | 333 |
Preflow& source(const Node& node) { |
334 | 334 |
_source = node; |
335 | 335 |
return *this; |
336 | 336 |
} |
337 | 337 |
|
338 | 338 |
/// \brief Sets the target node. |
339 | 339 |
/// |
340 | 340 |
/// Sets the target node. |
341 | 341 |
/// \return <tt>(*this)</tt> |
342 | 342 |
Preflow& target(const Node& node) { |
343 | 343 |
_target = node; |
344 | 344 |
return *this; |
345 | 345 |
} |
346 | 346 |
|
347 | 347 |
/// \brief Sets the elevator used by algorithm. |
348 | 348 |
/// |
349 | 349 |
/// Sets the elevator used by algorithm. |
350 | 350 |
/// If you don't use this function before calling \ref run() or |
351 | 351 |
/// \ref init(), an instance will be allocated automatically. |
352 | 352 |
/// The destructor deallocates this automatically allocated elevator, |
353 | 353 |
/// of course. |
354 | 354 |
/// \return <tt>(*this)</tt> |
355 | 355 |
Preflow& elevator(Elevator& elevator) { |
356 | 356 |
if (_local_level) { |
357 | 357 |
delete _level; |
358 | 358 |
_local_level = false; |
359 | 359 |
} |
360 | 360 |
_level = &elevator; |
361 | 361 |
return *this; |
362 | 362 |
} |
363 | 363 |
|
364 | 364 |
/// \brief Returns a const reference to the elevator. |
365 | 365 |
/// |
366 | 366 |
/// Returns a const reference to the elevator. |
367 | 367 |
/// |
368 | 368 |
/// \pre Either \ref run() or \ref init() must be called before |
369 | 369 |
/// using this function. |
370 | 370 |
const Elevator& elevator() const { |
371 | 371 |
return *_level; |
372 | 372 |
} |
373 | 373 |
|
374 |
/// \brief Sets the tolerance used by algorithm. |
|
374 |
/// \brief Sets the tolerance used by the algorithm. |
|
375 | 375 |
/// |
376 |
/// Sets the tolerance used by algorithm. |
|
376 |
/// Sets the tolerance object used by the algorithm. |
|
377 |
/// \return <tt>(*this)</tt> |
|
377 | 378 |
Preflow& tolerance(const Tolerance& tolerance) { |
378 | 379 |
_tolerance = tolerance; |
379 | 380 |
return *this; |
380 | 381 |
} |
381 | 382 |
|
382 | 383 |
/// \brief Returns a const reference to the tolerance. |
383 | 384 |
/// |
384 |
/// Returns a const reference to the tolerance |
|
385 |
/// Returns a const reference to the tolerance object used by |
|
386 |
/// the algorithm. |
|
385 | 387 |
const Tolerance& tolerance() const { |
386 | 388 |
return _tolerance; |
387 | 389 |
} |
388 | 390 |
|
389 | 391 |
/// \name Execution Control |
390 | 392 |
/// The simplest way to execute the preflow algorithm is to use |
391 | 393 |
/// \ref run() or \ref runMinCut().\n |
392 | 394 |
/// If you need more control on the initial solution or the execution, |
393 | 395 |
/// first you have to call one of the \ref init() functions, then |
394 | 396 |
/// \ref startFirstPhase() and if you need it \ref startSecondPhase(). |
395 | 397 |
|
396 | 398 |
///@{ |
397 | 399 |
|
398 | 400 |
/// \brief Initializes the internal data structures. |
399 | 401 |
/// |
400 | 402 |
/// Initializes the internal data structures and sets the initial |
401 | 403 |
/// flow to zero on each arc. |
402 | 404 |
void init() { |
403 | 405 |
createStructures(); |
404 | 406 |
|
405 | 407 |
_phase = true; |
406 | 408 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
407 | 409 |
(*_excess)[n] = 0; |
408 | 410 |
} |
409 | 411 |
|
410 | 412 |
for (ArcIt e(_graph); e != INVALID; ++e) { |
411 | 413 |
_flow->set(e, 0); |
412 | 414 |
} |
413 | 415 |
|
414 | 416 |
typename Digraph::template NodeMap<bool> reached(_graph, false); |
415 | 417 |
|
416 | 418 |
_level->initStart(); |
417 | 419 |
_level->initAddItem(_target); |
418 | 420 |
|
419 | 421 |
std::vector<Node> queue; |
420 | 422 |
reached[_source] = true; |
421 | 423 |
|
422 | 424 |
queue.push_back(_target); |
423 | 425 |
reached[_target] = true; |
424 | 426 |
while (!queue.empty()) { |
425 | 427 |
_level->initNewLevel(); |
426 | 428 |
std::vector<Node> nqueue; |
427 | 429 |
for (int i = 0; i < int(queue.size()); ++i) { |
428 | 430 |
Node n = queue[i]; |
429 | 431 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
430 | 432 |
Node u = _graph.source(e); |
431 | 433 |
if (!reached[u] && _tolerance.positive((*_capacity)[e])) { |
432 | 434 |
reached[u] = true; |
433 | 435 |
_level->initAddItem(u); |
434 | 436 |
nqueue.push_back(u); |
435 | 437 |
} |
436 | 438 |
} |
437 | 439 |
} |
438 | 440 |
queue.swap(nqueue); |
439 | 441 |
} |
440 | 442 |
_level->initFinish(); |
441 | 443 |
|
442 | 444 |
for (OutArcIt e(_graph, _source); e != INVALID; ++e) { |
443 | 445 |
if (_tolerance.positive((*_capacity)[e])) { |
444 | 446 |
Node u = _graph.target(e); |
445 | 447 |
if ((*_level)[u] == _level->maxLevel()) continue; |
446 | 448 |
_flow->set(e, (*_capacity)[e]); |
447 | 449 |
(*_excess)[u] += (*_capacity)[e]; |
448 | 450 |
if (u != _target && !_level->active(u)) { |
449 | 451 |
_level->activate(u); |
450 | 452 |
} |
451 | 453 |
} |
452 | 454 |
} |
453 | 455 |
} |
454 | 456 |
|
455 | 457 |
/// \brief Initializes the internal data structures using the |
456 | 458 |
/// given flow map. |
457 | 459 |
/// |
458 | 460 |
/// Initializes the internal data structures and sets the initial |
459 | 461 |
/// flow to the given \c flowMap. The \c flowMap should contain a |
460 | 462 |
/// flow or at least a preflow, i.e. at each node excluding the |
461 | 463 |
/// source node the incoming flow should greater or equal to the |
462 | 464 |
/// outgoing flow. |
463 | 465 |
/// \return \c false if the given \c flowMap is not a preflow. |
464 | 466 |
template <typename FlowMap> |
465 | 467 |
bool init(const FlowMap& flowMap) { |
466 | 468 |
createStructures(); |
467 | 469 |
|
468 | 470 |
for (ArcIt e(_graph); e != INVALID; ++e) { |
469 | 471 |
_flow->set(e, flowMap[e]); |
470 | 472 |
} |
471 | 473 |
|
472 | 474 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
473 | 475 |
Value excess = 0; |
474 | 476 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
475 | 477 |
excess += (*_flow)[e]; |
476 | 478 |
} |
477 | 479 |
for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
478 | 480 |
excess -= (*_flow)[e]; |
479 | 481 |
} |
480 | 482 |
if (excess < 0 && n != _source) return false; |
481 | 483 |
(*_excess)[n] = excess; |
482 | 484 |
} |
483 | 485 |
|
484 | 486 |
typename Digraph::template NodeMap<bool> reached(_graph, false); |
485 | 487 |
|
486 | 488 |
_level->initStart(); |
487 | 489 |
_level->initAddItem(_target); |
488 | 490 |
|
489 | 491 |
std::vector<Node> queue; |
490 | 492 |
reached[_source] = true; |
491 | 493 |
|
492 | 494 |
queue.push_back(_target); |
493 | 495 |
reached[_target] = true; |
494 | 496 |
while (!queue.empty()) { |
495 | 497 |
_level->initNewLevel(); |
496 | 498 |
std::vector<Node> nqueue; |
497 | 499 |
for (int i = 0; i < int(queue.size()); ++i) { |
498 | 500 |
Node n = queue[i]; |
499 | 501 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
500 | 502 |
Node u = _graph.source(e); |
501 | 503 |
if (!reached[u] && |
502 | 504 |
_tolerance.positive((*_capacity)[e] - (*_flow)[e])) { |
503 | 505 |
reached[u] = true; |
504 | 506 |
_level->initAddItem(u); |
505 | 507 |
nqueue.push_back(u); |
506 | 508 |
} |
507 | 509 |
} |
508 | 510 |
for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
509 | 511 |
Node v = _graph.target(e); |
510 | 512 |
if (!reached[v] && _tolerance.positive((*_flow)[e])) { |
511 | 513 |
reached[v] = true; |
512 | 514 |
_level->initAddItem(v); |
513 | 515 |
nqueue.push_back(v); |
514 | 516 |
} |
515 | 517 |
} |
516 | 518 |
} |
517 | 519 |
queue.swap(nqueue); |
518 | 520 |
} |
519 | 521 |
_level->initFinish(); |
520 | 522 |
|
521 | 523 |
for (OutArcIt e(_graph, _source); e != INVALID; ++e) { |
522 | 524 |
Value rem = (*_capacity)[e] - (*_flow)[e]; |
523 | 525 |
if (_tolerance.positive(rem)) { |
524 | 526 |
Node u = _graph.target(e); |
525 | 527 |
if ((*_level)[u] == _level->maxLevel()) continue; |
526 | 528 |
_flow->set(e, (*_capacity)[e]); |
527 | 529 |
(*_excess)[u] += rem; |
528 | 530 |
if (u != _target && !_level->active(u)) { |
529 | 531 |
_level->activate(u); |
530 | 532 |
} |
531 | 533 |
} |
532 | 534 |
} |
533 | 535 |
for (InArcIt e(_graph, _source); e != INVALID; ++e) { |
534 | 536 |
Value rem = (*_flow)[e]; |
535 | 537 |
if (_tolerance.positive(rem)) { |
536 | 538 |
Node v = _graph.source(e); |
537 | 539 |
if ((*_level)[v] == _level->maxLevel()) continue; |
538 | 540 |
_flow->set(e, 0); |
539 | 541 |
(*_excess)[v] += rem; |
540 | 542 |
if (v != _target && !_level->active(v)) { |
541 | 543 |
_level->activate(v); |
542 | 544 |
} |
543 | 545 |
} |
544 | 546 |
} |
545 | 547 |
return true; |
546 | 548 |
} |
547 | 549 |
|
548 | 550 |
/// \brief Starts the first phase of the preflow algorithm. |
549 | 551 |
/// |
550 | 552 |
/// The preflow algorithm consists of two phases, this method runs |
551 | 553 |
/// the first phase. After the first phase the maximum flow value |
552 | 554 |
/// and a minimum value cut can already be computed, although a |
553 | 555 |
/// maximum flow is not yet obtained. So after calling this method |
554 | 556 |
/// \ref flowValue() returns the value of a maximum flow and \ref |
555 | 557 |
/// minCut() returns a minimum cut. |
556 | 558 |
/// \pre One of the \ref init() functions must be called before |
557 | 559 |
/// using this function. |
558 | 560 |
void startFirstPhase() { |
559 | 561 |
_phase = true; |
560 | 562 |
|
561 | 563 |
Node n = _level->highestActive(); |
562 | 564 |
int level = _level->highestActiveLevel(); |
563 | 565 |
while (n != INVALID) { |
564 | 566 |
int num = _node_num; |
565 | 567 |
|
566 | 568 |
while (num > 0 && n != INVALID) { |
567 | 569 |
Value excess = (*_excess)[n]; |
568 | 570 |
int new_level = _level->maxLevel(); |
569 | 571 |
|
570 | 572 |
for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
571 | 573 |
Value rem = (*_capacity)[e] - (*_flow)[e]; |
572 | 574 |
if (!_tolerance.positive(rem)) continue; |
573 | 575 |
Node v = _graph.target(e); |
574 | 576 |
if ((*_level)[v] < level) { |
575 | 577 |
if (!_level->active(v) && v != _target) { |
576 | 578 |
_level->activate(v); |
577 | 579 |
} |
578 | 580 |
if (!_tolerance.less(rem, excess)) { |
579 | 581 |
_flow->set(e, (*_flow)[e] + excess); |
580 | 582 |
(*_excess)[v] += excess; |
581 | 583 |
excess = 0; |
582 | 584 |
goto no_more_push_1; |
583 | 585 |
} else { |
584 | 586 |
excess -= rem; |
585 | 587 |
(*_excess)[v] += rem; |
586 | 588 |
_flow->set(e, (*_capacity)[e]); |
587 | 589 |
} |
588 | 590 |
} else if (new_level > (*_level)[v]) { |
589 | 591 |
new_level = (*_level)[v]; |
590 | 592 |
} |
591 | 593 |
} |
592 | 594 |
|
593 | 595 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
594 | 596 |
Value rem = (*_flow)[e]; |
595 | 597 |
if (!_tolerance.positive(rem)) continue; |
596 | 598 |
Node v = _graph.source(e); |
597 | 599 |
if ((*_level)[v] < level) { |
598 | 600 |
if (!_level->active(v) && v != _target) { |
599 | 601 |
_level->activate(v); |
600 | 602 |
} |
601 | 603 |
if (!_tolerance.less(rem, excess)) { |
602 | 604 |
_flow->set(e, (*_flow)[e] - excess); |
603 | 605 |
(*_excess)[v] += excess; |
604 | 606 |
excess = 0; |
605 | 607 |
goto no_more_push_1; |
606 | 608 |
} else { |
607 | 609 |
excess -= rem; |
608 | 610 |
(*_excess)[v] += rem; |
609 | 611 |
_flow->set(e, 0); |
610 | 612 |
} |
611 | 613 |
} else if (new_level > (*_level)[v]) { |
612 | 614 |
new_level = (*_level)[v]; |
613 | 615 |
} |
614 | 616 |
} |
615 | 617 |
|
616 | 618 |
no_more_push_1: |
617 | 619 |
|
618 | 620 |
(*_excess)[n] = excess; |
619 | 621 |
|
620 | 622 |
if (excess != 0) { |
621 | 623 |
if (new_level + 1 < _level->maxLevel()) { |
622 | 624 |
_level->liftHighestActive(new_level + 1); |
623 | 625 |
} else { |
624 | 626 |
_level->liftHighestActiveToTop(); |
625 | 627 |
} |
626 | 628 |
if (_level->emptyLevel(level)) { |
627 | 629 |
_level->liftToTop(level); |
628 | 630 |
} |
629 | 631 |
} else { |
630 | 632 |
_level->deactivate(n); |
631 | 633 |
} |
632 | 634 |
|
633 | 635 |
n = _level->highestActive(); |
634 | 636 |
level = _level->highestActiveLevel(); |
635 | 637 |
--num; |
636 | 638 |
} |
637 | 639 |
|
638 | 640 |
num = _node_num * 20; |
639 | 641 |
while (num > 0 && n != INVALID) { |
640 | 642 |
Value excess = (*_excess)[n]; |
641 | 643 |
int new_level = _level->maxLevel(); |
642 | 644 |
|
643 | 645 |
for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
644 | 646 |
Value rem = (*_capacity)[e] - (*_flow)[e]; |
645 | 647 |
if (!_tolerance.positive(rem)) continue; |
646 | 648 |
Node v = _graph.target(e); |
647 | 649 |
if ((*_level)[v] < level) { |
648 | 650 |
if (!_level->active(v) && v != _target) { |
649 | 651 |
_level->activate(v); |
650 | 652 |
} |
651 | 653 |
if (!_tolerance.less(rem, excess)) { |
652 | 654 |
_flow->set(e, (*_flow)[e] + excess); |
653 | 655 |
(*_excess)[v] += excess; |
654 | 656 |
excess = 0; |
655 | 657 |
goto no_more_push_2; |
656 | 658 |
} else { |
657 | 659 |
excess -= rem; |
658 | 660 |
(*_excess)[v] += rem; |
659 | 661 |
_flow->set(e, (*_capacity)[e]); |
660 | 662 |
} |
661 | 663 |
} else if (new_level > (*_level)[v]) { |
662 | 664 |
new_level = (*_level)[v]; |
663 | 665 |
} |
664 | 666 |
} |
665 | 667 |
|
666 | 668 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
667 | 669 |
Value rem = (*_flow)[e]; |
668 | 670 |
if (!_tolerance.positive(rem)) continue; |
669 | 671 |
Node v = _graph.source(e); |
670 | 672 |
if ((*_level)[v] < level) { |
671 | 673 |
if (!_level->active(v) && v != _target) { |
672 | 674 |
_level->activate(v); |
673 | 675 |
} |
674 | 676 |
if (!_tolerance.less(rem, excess)) { |
675 | 677 |
_flow->set(e, (*_flow)[e] - excess); |
676 | 678 |
(*_excess)[v] += excess; |
677 | 679 |
excess = 0; |
678 | 680 |
goto no_more_push_2; |
679 | 681 |
} else { |
680 | 682 |
excess -= rem; |
681 | 683 |
(*_excess)[v] += rem; |
682 | 684 |
_flow->set(e, 0); |
683 | 685 |
} |
684 | 686 |
} else if (new_level > (*_level)[v]) { |
685 | 687 |
new_level = (*_level)[v]; |
686 | 688 |
} |
687 | 689 |
} |
688 | 690 |
|
689 | 691 |
no_more_push_2: |
690 | 692 |
|
691 | 693 |
(*_excess)[n] = excess; |
692 | 694 |
|
693 | 695 |
if (excess != 0) { |
694 | 696 |
if (new_level + 1 < _level->maxLevel()) { |
695 | 697 |
_level->liftActiveOn(level, new_level + 1); |
696 | 698 |
} else { |
697 | 699 |
_level->liftActiveToTop(level); |
698 | 700 |
} |
699 | 701 |
if (_level->emptyLevel(level)) { |
700 | 702 |
_level->liftToTop(level); |
701 | 703 |
} |
702 | 704 |
} else { |
703 | 705 |
_level->deactivate(n); |
704 | 706 |
} |
705 | 707 |
|
706 | 708 |
while (level >= 0 && _level->activeFree(level)) { |
707 | 709 |
--level; |
708 | 710 |
} |
709 | 711 |
if (level == -1) { |
710 | 712 |
n = _level->highestActive(); |
711 | 713 |
level = _level->highestActiveLevel(); |
712 | 714 |
} else { |
713 | 715 |
n = _level->activeOn(level); |
714 | 716 |
} |
715 | 717 |
--num; |
716 | 718 |
} |
717 | 719 |
} |
718 | 720 |
} |
719 | 721 |
|
720 | 722 |
/// \brief Starts the second phase of the preflow algorithm. |
721 | 723 |
/// |
722 | 724 |
/// The preflow algorithm consists of two phases, this method runs |
723 | 725 |
/// the second phase. After calling one of the \ref init() functions |
724 | 726 |
/// and \ref startFirstPhase() and then \ref startSecondPhase(), |
725 | 727 |
/// \ref flowMap() returns a maximum flow, \ref flowValue() returns the |
726 | 728 |
/// value of a maximum flow, \ref minCut() returns a minimum cut |
727 | 729 |
/// \pre One of the \ref init() functions and \ref startFirstPhase() |
728 | 730 |
/// must be called before using this function. |
729 | 731 |
void startSecondPhase() { |
730 | 732 |
_phase = false; |
731 | 733 |
|
732 | 734 |
typename Digraph::template NodeMap<bool> reached(_graph); |
733 | 735 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
734 | 736 |
reached[n] = (*_level)[n] < _level->maxLevel(); |
735 | 737 |
} |
736 | 738 |
|
737 | 739 |
_level->initStart(); |
738 | 740 |
_level->initAddItem(_source); |
739 | 741 |
|
740 | 742 |
std::vector<Node> queue; |
741 | 743 |
queue.push_back(_source); |
742 | 744 |
reached[_source] = true; |
743 | 745 |
|
744 | 746 |
while (!queue.empty()) { |
745 | 747 |
_level->initNewLevel(); |
746 | 748 |
std::vector<Node> nqueue; |
747 | 749 |
for (int i = 0; i < int(queue.size()); ++i) { |
748 | 750 |
Node n = queue[i]; |
749 | 751 |
for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
750 | 752 |
Node v = _graph.target(e); |
751 | 753 |
if (!reached[v] && _tolerance.positive((*_flow)[e])) { |
752 | 754 |
reached[v] = true; |
753 | 755 |
_level->initAddItem(v); |
754 | 756 |
nqueue.push_back(v); |
755 | 757 |
} |
756 | 758 |
} |
757 | 759 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
758 | 760 |
Node u = _graph.source(e); |
759 | 761 |
if (!reached[u] && |
760 | 762 |
_tolerance.positive((*_capacity)[e] - (*_flow)[e])) { |
761 | 763 |
reached[u] = true; |
762 | 764 |
_level->initAddItem(u); |
763 | 765 |
nqueue.push_back(u); |
764 | 766 |
} |
765 | 767 |
} |
766 | 768 |
} |
767 | 769 |
queue.swap(nqueue); |
768 | 770 |
} |
769 | 771 |
_level->initFinish(); |
770 | 772 |
|
771 | 773 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
772 | 774 |
if (!reached[n]) { |
773 | 775 |
_level->dirtyTopButOne(n); |
774 | 776 |
} else if ((*_excess)[n] > 0 && _target != n) { |
775 | 777 |
_level->activate(n); |
776 | 778 |
} |
777 | 779 |
} |
778 | 780 |
|
779 | 781 |
Node n; |
780 | 782 |
while ((n = _level->highestActive()) != INVALID) { |
781 | 783 |
Value excess = (*_excess)[n]; |
782 | 784 |
int level = _level->highestActiveLevel(); |
783 | 785 |
int new_level = _level->maxLevel(); |
784 | 786 |
|
785 | 787 |
for (OutArcIt e(_graph, n); e != INVALID; ++e) { |
786 | 788 |
Value rem = (*_capacity)[e] - (*_flow)[e]; |
787 | 789 |
if (!_tolerance.positive(rem)) continue; |
788 | 790 |
Node v = _graph.target(e); |
789 | 791 |
if ((*_level)[v] < level) { |
790 | 792 |
if (!_level->active(v) && v != _source) { |
791 | 793 |
_level->activate(v); |
792 | 794 |
} |
793 | 795 |
if (!_tolerance.less(rem, excess)) { |
794 | 796 |
_flow->set(e, (*_flow)[e] + excess); |
795 | 797 |
(*_excess)[v] += excess; |
796 | 798 |
excess = 0; |
797 | 799 |
goto no_more_push; |
798 | 800 |
} else { |
799 | 801 |
excess -= rem; |
800 | 802 |
(*_excess)[v] += rem; |
801 | 803 |
_flow->set(e, (*_capacity)[e]); |
802 | 804 |
} |
803 | 805 |
} else if (new_level > (*_level)[v]) { |
804 | 806 |
new_level = (*_level)[v]; |
805 | 807 |
} |
806 | 808 |
} |
807 | 809 |
|
808 | 810 |
for (InArcIt e(_graph, n); e != INVALID; ++e) { |
809 | 811 |
Value rem = (*_flow)[e]; |
810 | 812 |
if (!_tolerance.positive(rem)) continue; |
811 | 813 |
Node v = _graph.source(e); |
812 | 814 |
if ((*_level)[v] < level) { |
813 | 815 |
if (!_level->active(v) && v != _source) { |
814 | 816 |
_level->activate(v); |
815 | 817 |
} |
816 | 818 |
if (!_tolerance.less(rem, excess)) { |
817 | 819 |
_flow->set(e, (*_flow)[e] - excess); |
818 | 820 |
(*_excess)[v] += excess; |
819 | 821 |
excess = 0; |
820 | 822 |
goto no_more_push; |
821 | 823 |
} else { |
822 | 824 |
excess -= rem; |
823 | 825 |
(*_excess)[v] += rem; |
824 | 826 |
_flow->set(e, 0); |
825 | 827 |
} |
826 | 828 |
} else if (new_level > (*_level)[v]) { |
827 | 829 |
new_level = (*_level)[v]; |
828 | 830 |
} |
829 | 831 |
} |
830 | 832 |
|
831 | 833 |
no_more_push: |
832 | 834 |
|
833 | 835 |
(*_excess)[n] = excess; |
834 | 836 |
|
835 | 837 |
if (excess != 0) { |
836 | 838 |
if (new_level + 1 < _level->maxLevel()) { |
837 | 839 |
_level->liftHighestActive(new_level + 1); |
838 | 840 |
} else { |
839 | 841 |
// Calculation error |
840 | 842 |
_level->liftHighestActiveToTop(); |
841 | 843 |
} |
842 | 844 |
if (_level->emptyLevel(level)) { |
843 | 845 |
// Calculation error |
844 | 846 |
_level->liftToTop(level); |
845 | 847 |
} |
846 | 848 |
} else { |
847 | 849 |
_level->deactivate(n); |
848 | 850 |
} |
849 | 851 |
|
850 | 852 |
} |
851 | 853 |
} |
852 | 854 |
|
853 | 855 |
/// \brief Runs the preflow algorithm. |
854 | 856 |
/// |
855 | 857 |
/// Runs the preflow algorithm. |
856 | 858 |
/// \note pf.run() is just a shortcut of the following code. |
857 | 859 |
/// \code |
858 | 860 |
/// pf.init(); |
859 | 861 |
/// pf.startFirstPhase(); |
860 | 862 |
/// pf.startSecondPhase(); |
861 | 863 |
/// \endcode |
862 | 864 |
void run() { |
863 | 865 |
init(); |
864 | 866 |
startFirstPhase(); |
865 | 867 |
startSecondPhase(); |
866 | 868 |
} |
867 | 869 |
|
868 | 870 |
/// \brief Runs the preflow algorithm to compute the minimum cut. |
869 | 871 |
/// |
870 | 872 |
/// Runs the preflow algorithm to compute the minimum cut. |
871 | 873 |
/// \note pf.runMinCut() is just a shortcut of the following code. |
872 | 874 |
/// \code |
873 | 875 |
/// pf.init(); |
874 | 876 |
/// pf.startFirstPhase(); |
875 | 877 |
/// \endcode |
876 | 878 |
void runMinCut() { |
877 | 879 |
init(); |
878 | 880 |
startFirstPhase(); |
879 | 881 |
} |
880 | 882 |
|
881 | 883 |
/// @} |
882 | 884 |
|
883 | 885 |
/// \name Query Functions |
884 | 886 |
/// The results of the preflow algorithm can be obtained using these |
885 | 887 |
/// functions.\n |
886 | 888 |
/// Either one of the \ref run() "run*()" functions or one of the |
887 | 889 |
/// \ref startFirstPhase() "start*()" functions should be called |
888 | 890 |
/// before using them. |
889 | 891 |
|
890 | 892 |
///@{ |
891 | 893 |
|
892 | 894 |
/// \brief Returns the value of the maximum flow. |
893 | 895 |
/// |
894 | 896 |
/// Returns the value of the maximum flow by returning the excess |
895 | 897 |
/// of the target node. This value equals to the value of |
896 | 898 |
/// the maximum flow already after the first phase of the algorithm. |
897 | 899 |
/// |
898 | 900 |
/// \pre Either \ref run() or \ref init() must be called before |
899 | 901 |
/// using this function. |
900 | 902 |
Value flowValue() const { |
901 | 903 |
return (*_excess)[_target]; |
902 | 904 |
} |
903 | 905 |
|
904 | 906 |
/// \brief Returns the flow value on the given arc. |
905 | 907 |
/// |
906 | 908 |
/// Returns the flow value on the given arc. This method can |
907 | 909 |
/// be called after the second phase of the algorithm. |
908 | 910 |
/// |
909 | 911 |
/// \pre Either \ref run() or \ref init() must be called before |
910 | 912 |
/// using this function. |
911 | 913 |
Value flow(const Arc& arc) const { |
912 | 914 |
return (*_flow)[arc]; |
913 | 915 |
} |
914 | 916 |
|
915 | 917 |
/// \brief Returns a const reference to the flow map. |
916 | 918 |
/// |
917 | 919 |
/// Returns a const reference to the arc map storing the found flow. |
918 | 920 |
/// This method can be called after the second phase of the algorithm. |
919 | 921 |
/// |
920 | 922 |
/// \pre Either \ref run() or \ref init() must be called before |
921 | 923 |
/// using this function. |
922 | 924 |
const FlowMap& flowMap() const { |
923 | 925 |
return *_flow; |
924 | 926 |
} |
925 | 927 |
|
926 | 928 |
/// \brief Returns \c true when the node is on the source side of the |
927 | 929 |
/// minimum cut. |
928 | 930 |
/// |
929 | 931 |
/// Returns true when the node is on the source side of the found |
930 | 932 |
/// minimum cut. This method can be called both after running \ref |
931 | 933 |
/// startFirstPhase() and \ref startSecondPhase(). |
932 | 934 |
/// |
933 | 935 |
/// \pre Either \ref run() or \ref init() must be called before |
934 | 936 |
/// using this function. |
935 | 937 |
bool minCut(const Node& node) const { |
936 | 938 |
return ((*_level)[node] == _level->maxLevel()) == _phase; |
937 | 939 |
} |
938 | 940 |
|
939 | 941 |
/// \brief Gives back a minimum value cut. |
940 | 942 |
/// |
941 | 943 |
/// Sets \c cutMap to the characteristic vector of a minimum value |
942 | 944 |
/// cut. \c cutMap should be a \ref concepts::WriteMap "writable" |
943 | 945 |
/// node map with \c bool (or convertible) value type. |
944 | 946 |
/// |
945 | 947 |
/// This method can be called both after running \ref startFirstPhase() |
946 | 948 |
/// and \ref startSecondPhase(). The result after the second phase |
947 | 949 |
/// could be slightly different if inexact computation is used. |
948 | 950 |
/// |
949 | 951 |
/// \note This function calls \ref minCut() for each node, so it runs in |
950 | 952 |
/// O(n) time. |
951 | 953 |
/// |
952 | 954 |
/// \pre Either \ref run() or \ref init() must be called before |
953 | 955 |
/// using this function. |
954 | 956 |
template <typename CutMap> |
955 | 957 |
void minCutMap(CutMap& cutMap) const { |
956 | 958 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
957 | 959 |
cutMap.set(n, minCut(n)); |
958 | 960 |
} |
959 | 961 |
} |
960 | 962 |
|
961 | 963 |
/// @} |
962 | 964 |
}; |
963 | 965 |
} |
964 | 966 |
|
965 | 967 |
#endif |
1 | 1 |
/* -*- mode: C++; indent-tabs-mode: nil; -*- |
2 | 2 |
* |
3 | 3 |
* This file is a part of LEMON, a generic C++ optimization library. |
4 | 4 |
* |
5 | 5 |
* Copyright (C) 2003-2009 |
6 | 6 |
* Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport |
7 | 7 |
* (Egervary Research Group on Combinatorial Optimization, EGRES). |
8 | 8 |
* |
9 | 9 |
* Permission to use, modify and distribute this software is granted |
10 | 10 |
* provided that this copyright notice appears in all copies. For |
11 | 11 |
* precise terms see the accompanying LICENSE file. |
12 | 12 |
* |
13 | 13 |
* This software is provided "AS IS" with no warranty of any kind, |
14 | 14 |
* express or implied, and with no claim as to its suitability for any |
15 | 15 |
* purpose. |
16 | 16 |
* |
17 | 17 |
*/ |
18 | 18 |
|
19 | 19 |
#include <iostream> |
20 | 20 |
|
21 | 21 |
#include "test_tools.h" |
22 | 22 |
#include <lemon/list_graph.h> |
23 | 23 |
#include <lemon/circulation.h> |
24 | 24 |
#include <lemon/lgf_reader.h> |
25 | 25 |
#include <lemon/concepts/digraph.h> |
26 | 26 |
#include <lemon/concepts/maps.h> |
27 | 27 |
|
28 | 28 |
using namespace lemon; |
29 | 29 |
|
30 | 30 |
char test_lgf[] = |
31 | 31 |
"@nodes\n" |
32 | 32 |
"label\n" |
33 | 33 |
"0\n" |
34 | 34 |
"1\n" |
35 | 35 |
"2\n" |
36 | 36 |
"3\n" |
37 | 37 |
"4\n" |
38 | 38 |
"5\n" |
39 | 39 |
"@arcs\n" |
40 | 40 |
" lcap ucap\n" |
41 | 41 |
"0 1 2 10\n" |
42 | 42 |
"0 2 2 6\n" |
43 | 43 |
"1 3 4 7\n" |
44 | 44 |
"1 4 0 5\n" |
45 | 45 |
"2 4 1 3\n" |
46 | 46 |
"3 5 3 8\n" |
47 | 47 |
"4 5 3 7\n" |
48 | 48 |
"@attributes\n" |
49 | 49 |
"source 0\n" |
50 | 50 |
"sink 5\n"; |
51 | 51 |
|
52 | 52 |
void checkCirculationCompile() |
53 | 53 |
{ |
54 | 54 |
typedef int VType; |
55 | 55 |
typedef concepts::Digraph Digraph; |
56 | 56 |
|
57 | 57 |
typedef Digraph::Node Node; |
58 | 58 |
typedef Digraph::Arc Arc; |
59 | 59 |
typedef concepts::ReadMap<Arc,VType> CapMap; |
60 | 60 |
typedef concepts::ReadMap<Node,VType> SupplyMap; |
61 | 61 |
typedef concepts::ReadWriteMap<Arc,VType> FlowMap; |
62 | 62 |
typedef concepts::WriteMap<Node,bool> BarrierMap; |
63 | 63 |
|
64 | 64 |
typedef Elevator<Digraph, Digraph::Node> Elev; |
65 | 65 |
typedef LinkedElevator<Digraph, Digraph::Node> LinkedElev; |
66 | 66 |
|
67 | 67 |
Digraph g; |
68 | 68 |
Node n; |
69 | 69 |
Arc a; |
70 | 70 |
CapMap lcap, ucap; |
71 | 71 |
SupplyMap supply; |
72 | 72 |
FlowMap flow; |
73 | 73 |
BarrierMap bar; |
74 | 74 |
VType v; |
75 | 75 |
bool b; |
76 | 76 |
|
77 | 77 |
typedef Circulation<Digraph, CapMap, CapMap, SupplyMap> |
78 | 78 |
::SetFlowMap<FlowMap> |
79 | 79 |
::SetElevator<Elev> |
80 | 80 |
::SetStandardElevator<LinkedElev> |
81 | 81 |
::Create CirculationType; |
82 | 82 |
CirculationType circ_test(g, lcap, ucap, supply); |
83 | 83 |
const CirculationType& const_circ_test = circ_test; |
84 | 84 |
|
85 | 85 |
circ_test |
86 | 86 |
.lowerMap(lcap) |
87 | 87 |
.upperMap(ucap) |
88 | 88 |
.supplyMap(supply) |
89 | 89 |
.flowMap(flow); |
90 | 90 |
|
91 |
const CirculationType::Elevator& elev = const_circ_test.elevator(); |
|
92 |
circ_test.elevator(const_cast<CirculationType::Elevator&>(elev)); |
|
93 |
CirculationType::Tolerance tol = const_circ_test.tolerance(); |
|
94 |
circ_test.tolerance(tol); |
|
95 |
|
|
91 | 96 |
circ_test.init(); |
92 | 97 |
circ_test.greedyInit(); |
93 | 98 |
circ_test.start(); |
94 | 99 |
circ_test.run(); |
95 | 100 |
|
96 | 101 |
v = const_circ_test.flow(a); |
97 | 102 |
const FlowMap& fm = const_circ_test.flowMap(); |
98 | 103 |
b = const_circ_test.barrier(n); |
99 | 104 |
const_circ_test.barrierMap(bar); |
100 | 105 |
|
101 | 106 |
ignore_unused_variable_warning(fm); |
102 | 107 |
} |
103 | 108 |
|
104 | 109 |
template <class G, class LM, class UM, class DM> |
105 | 110 |
void checkCirculation(const G& g, const LM& lm, const UM& um, |
106 | 111 |
const DM& dm, bool find) |
107 | 112 |
{ |
108 | 113 |
Circulation<G, LM, UM, DM> circ(g, lm, um, dm); |
109 | 114 |
bool ret = circ.run(); |
110 | 115 |
if (find) { |
111 | 116 |
check(ret, "A feasible solution should have been found."); |
112 | 117 |
check(circ.checkFlow(), "The found flow is corrupt."); |
113 | 118 |
check(!circ.checkBarrier(), "A barrier should not have been found."); |
114 | 119 |
} else { |
115 | 120 |
check(!ret, "A feasible solution should not have been found."); |
116 | 121 |
check(circ.checkBarrier(), "The found barrier is corrupt."); |
117 | 122 |
} |
118 | 123 |
} |
119 | 124 |
|
120 | 125 |
int main (int, char*[]) |
121 | 126 |
{ |
122 | 127 |
typedef ListDigraph Digraph; |
123 | 128 |
DIGRAPH_TYPEDEFS(Digraph); |
124 | 129 |
|
125 | 130 |
Digraph g; |
126 | 131 |
IntArcMap lo(g), up(g); |
127 | 132 |
IntNodeMap delta(g, 0); |
128 | 133 |
Node s, t; |
129 | 134 |
|
130 | 135 |
std::istringstream input(test_lgf); |
131 | 136 |
DigraphReader<Digraph>(g,input). |
132 | 137 |
arcMap("lcap", lo). |
133 | 138 |
arcMap("ucap", up). |
134 | 139 |
node("source",s). |
135 | 140 |
node("sink",t). |
136 | 141 |
run(); |
137 | 142 |
|
138 | 143 |
delta[s] = 7; delta[t] = -7; |
139 | 144 |
checkCirculation(g, lo, up, delta, true); |
140 | 145 |
|
141 | 146 |
delta[s] = 13; delta[t] = -13; |
142 | 147 |
checkCirculation(g, lo, up, delta, true); |
143 | 148 |
|
144 | 149 |
delta[s] = 6; delta[t] = -6; |
145 | 150 |
checkCirculation(g, lo, up, delta, false); |
146 | 151 |
|
147 | 152 |
delta[s] = 14; delta[t] = -14; |
148 | 153 |
checkCirculation(g, lo, up, delta, false); |
149 | 154 |
|
150 | 155 |
delta[s] = 7; delta[t] = -13; |
151 | 156 |
checkCirculation(g, lo, up, delta, true); |
152 | 157 |
|
153 | 158 |
delta[s] = 5; delta[t] = -15; |
154 | 159 |
checkCirculation(g, lo, up, delta, true); |
155 | 160 |
|
156 | 161 |
delta[s] = 10; delta[t] = -11; |
157 | 162 |
checkCirculation(g, lo, up, delta, true); |
158 | 163 |
|
159 | 164 |
delta[s] = 11; delta[t] = -10; |
160 | 165 |
checkCirculation(g, lo, up, delta, false); |
161 | 166 |
|
162 | 167 |
return 0; |
163 | 168 |
} |
1 | 1 |
/* -*- mode: C++; indent-tabs-mode: nil; -*- |
2 | 2 |
* |
3 | 3 |
* This file is a part of LEMON, a generic C++ optimization library. |
4 | 4 |
* |
5 | 5 |
* Copyright (C) 2003-2009 |
6 | 6 |
* Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport |
7 | 7 |
* (Egervary Research Group on Combinatorial Optimization, EGRES). |
8 | 8 |
* |
9 | 9 |
* Permission to use, modify and distribute this software is granted |
10 | 10 |
* provided that this copyright notice appears in all copies. For |
11 | 11 |
* precise terms see the accompanying LICENSE file. |
12 | 12 |
* |
13 | 13 |
* This software is provided "AS IS" with no warranty of any kind, |
14 | 14 |
* express or implied, and with no claim as to its suitability for any |
15 | 15 |
* purpose. |
16 | 16 |
* |
17 | 17 |
*/ |
18 | 18 |
|
19 | 19 |
#include <iostream> |
20 | 20 |
|
21 | 21 |
#include "test_tools.h" |
22 | 22 |
#include <lemon/smart_graph.h> |
23 | 23 |
#include <lemon/preflow.h> |
24 | 24 |
#include <lemon/concepts/digraph.h> |
25 | 25 |
#include <lemon/concepts/maps.h> |
26 | 26 |
#include <lemon/lgf_reader.h> |
27 | 27 |
#include <lemon/elevator.h> |
28 | 28 |
|
29 | 29 |
using namespace lemon; |
30 | 30 |
|
31 | 31 |
char test_lgf[] = |
32 | 32 |
"@nodes\n" |
33 | 33 |
"label\n" |
34 | 34 |
"0\n" |
35 | 35 |
"1\n" |
36 | 36 |
"2\n" |
37 | 37 |
"3\n" |
38 | 38 |
"4\n" |
39 | 39 |
"5\n" |
40 | 40 |
"6\n" |
41 | 41 |
"7\n" |
42 | 42 |
"8\n" |
43 | 43 |
"9\n" |
44 | 44 |
"@arcs\n" |
45 | 45 |
" label capacity\n" |
46 | 46 |
"0 1 0 20\n" |
47 | 47 |
"0 2 1 0\n" |
48 | 48 |
"1 1 2 3\n" |
49 | 49 |
"1 2 3 8\n" |
50 | 50 |
"1 3 4 8\n" |
51 | 51 |
"2 5 5 5\n" |
52 | 52 |
"3 2 6 5\n" |
53 | 53 |
"3 5 7 5\n" |
54 | 54 |
"3 6 8 5\n" |
55 | 55 |
"4 3 9 3\n" |
56 | 56 |
"5 7 10 3\n" |
57 | 57 |
"5 6 11 10\n" |
58 | 58 |
"5 8 12 10\n" |
59 | 59 |
"6 8 13 8\n" |
60 | 60 |
"8 9 14 20\n" |
61 | 61 |
"8 1 15 5\n" |
62 | 62 |
"9 5 16 5\n" |
63 | 63 |
"@attributes\n" |
64 | 64 |
"source 1\n" |
65 | 65 |
"target 8\n"; |
66 | 66 |
|
67 | 67 |
void checkPreflowCompile() |
68 | 68 |
{ |
69 | 69 |
typedef int VType; |
70 | 70 |
typedef concepts::Digraph Digraph; |
71 | 71 |
|
72 | 72 |
typedef Digraph::Node Node; |
73 | 73 |
typedef Digraph::Arc Arc; |
74 | 74 |
typedef concepts::ReadMap<Arc,VType> CapMap; |
75 | 75 |
typedef concepts::ReadWriteMap<Arc,VType> FlowMap; |
76 | 76 |
typedef concepts::WriteMap<Node,bool> CutMap; |
77 | 77 |
|
78 | 78 |
typedef Elevator<Digraph, Digraph::Node> Elev; |
79 | 79 |
typedef LinkedElevator<Digraph, Digraph::Node> LinkedElev; |
80 | 80 |
|
81 | 81 |
Digraph g; |
82 | 82 |
Node n; |
83 | 83 |
Arc e; |
84 | 84 |
CapMap cap; |
85 | 85 |
FlowMap flow; |
86 | 86 |
CutMap cut; |
87 | 87 |
VType v; |
88 | 88 |
bool b; |
89 | 89 |
|
90 | 90 |
typedef Preflow<Digraph, CapMap> |
91 | 91 |
::SetFlowMap<FlowMap> |
92 | 92 |
::SetElevator<Elev> |
93 | 93 |
::SetStandardElevator<LinkedElev> |
94 | 94 |
::Create PreflowType; |
95 | 95 |
PreflowType preflow_test(g, cap, n, n); |
96 | 96 |
const PreflowType& const_preflow_test = preflow_test; |
97 | 97 |
|
98 |
const PreflowType::Elevator& elev = const_preflow_test.elevator(); |
|
99 |
preflow_test.elevator(const_cast<PreflowType::Elevator&>(elev)); |
|
100 |
PreflowType::Tolerance tol = const_preflow_test.tolerance(); |
|
101 |
preflow_test.tolerance(tol); |
|
102 |
|
|
98 | 103 |
preflow_test |
99 | 104 |
.capacityMap(cap) |
100 | 105 |
.flowMap(flow) |
101 | 106 |
.source(n) |
102 | 107 |
.target(n); |
103 | 108 |
|
104 | 109 |
preflow_test.init(); |
105 | 110 |
preflow_test.init(cap); |
106 | 111 |
preflow_test.startFirstPhase(); |
107 | 112 |
preflow_test.startSecondPhase(); |
108 | 113 |
preflow_test.run(); |
109 | 114 |
preflow_test.runMinCut(); |
110 | 115 |
|
111 | 116 |
v = const_preflow_test.flowValue(); |
112 | 117 |
v = const_preflow_test.flow(e); |
113 | 118 |
const FlowMap& fm = const_preflow_test.flowMap(); |
114 | 119 |
b = const_preflow_test.minCut(n); |
115 | 120 |
const_preflow_test.minCutMap(cut); |
116 | 121 |
|
117 | 122 |
ignore_unused_variable_warning(fm); |
118 | 123 |
} |
119 | 124 |
|
120 | 125 |
int cutValue (const SmartDigraph& g, |
121 | 126 |
const SmartDigraph::NodeMap<bool>& cut, |
122 | 127 |
const SmartDigraph::ArcMap<int>& cap) { |
123 | 128 |
|
124 | 129 |
int c=0; |
125 | 130 |
for(SmartDigraph::ArcIt e(g); e!=INVALID; ++e) { |
126 | 131 |
if (cut[g.source(e)] && !cut[g.target(e)]) c+=cap[e]; |
127 | 132 |
} |
128 | 133 |
return c; |
129 | 134 |
} |
130 | 135 |
|
131 | 136 |
bool checkFlow(const SmartDigraph& g, |
132 | 137 |
const SmartDigraph::ArcMap<int>& flow, |
133 | 138 |
const SmartDigraph::ArcMap<int>& cap, |
134 | 139 |
SmartDigraph::Node s, SmartDigraph::Node t) { |
135 | 140 |
|
136 | 141 |
for (SmartDigraph::ArcIt e(g); e != INVALID; ++e) { |
137 | 142 |
if (flow[e] < 0 || flow[e] > cap[e]) return false; |
138 | 143 |
} |
139 | 144 |
|
140 | 145 |
for (SmartDigraph::NodeIt n(g); n != INVALID; ++n) { |
141 | 146 |
if (n == s || n == t) continue; |
142 | 147 |
int sum = 0; |
143 | 148 |
for (SmartDigraph::OutArcIt e(g, n); e != INVALID; ++e) { |
144 | 149 |
sum += flow[e]; |
145 | 150 |
} |
146 | 151 |
for (SmartDigraph::InArcIt e(g, n); e != INVALID; ++e) { |
147 | 152 |
sum -= flow[e]; |
148 | 153 |
} |
149 | 154 |
if (sum != 0) return false; |
150 | 155 |
} |
151 | 156 |
return true; |
152 | 157 |
} |
153 | 158 |
|
154 | 159 |
int main() { |
155 | 160 |
|
156 | 161 |
typedef SmartDigraph Digraph; |
157 | 162 |
|
158 | 163 |
typedef Digraph::Node Node; |
159 | 164 |
typedef Digraph::NodeIt NodeIt; |
160 | 165 |
typedef Digraph::ArcIt ArcIt; |
161 | 166 |
typedef Digraph::ArcMap<int> CapMap; |
162 | 167 |
typedef Digraph::ArcMap<int> FlowMap; |
163 | 168 |
typedef Digraph::NodeMap<bool> CutMap; |
164 | 169 |
|
165 | 170 |
typedef Preflow<Digraph, CapMap> PType; |
166 | 171 |
|
167 | 172 |
Digraph g; |
168 | 173 |
Node s, t; |
169 | 174 |
CapMap cap(g); |
170 | 175 |
std::istringstream input(test_lgf); |
171 | 176 |
DigraphReader<Digraph>(g,input). |
172 | 177 |
arcMap("capacity", cap). |
173 | 178 |
node("source",s). |
174 | 179 |
node("target",t). |
175 | 180 |
run(); |
176 | 181 |
|
177 | 182 |
PType preflow_test(g, cap, s, t); |
178 | 183 |
preflow_test.run(); |
179 | 184 |
|
180 | 185 |
check(checkFlow(g, preflow_test.flowMap(), cap, s, t), |
181 | 186 |
"The flow is not feasible."); |
182 | 187 |
|
183 | 188 |
CutMap min_cut(g); |
184 | 189 |
preflow_test.minCutMap(min_cut); |
185 | 190 |
int min_cut_value=cutValue(g,min_cut,cap); |
186 | 191 |
|
187 | 192 |
check(preflow_test.flowValue() == min_cut_value, |
188 | 193 |
"The max flow value is not equal to the three min cut values."); |
189 | 194 |
|
190 | 195 |
FlowMap flow(g); |
191 | 196 |
for(ArcIt e(g); e!=INVALID; ++e) flow[e] = preflow_test.flowMap()[e]; |
192 | 197 |
|
193 | 198 |
int flow_value=preflow_test.flowValue(); |
194 | 199 |
|
195 | 200 |
for(ArcIt e(g); e!=INVALID; ++e) cap[e]=2*cap[e]; |
196 | 201 |
preflow_test.init(flow); |
197 | 202 |
preflow_test.startFirstPhase(); |
198 | 203 |
|
199 | 204 |
CutMap min_cut1(g); |
200 | 205 |
preflow_test.minCutMap(min_cut1); |
201 | 206 |
min_cut_value=cutValue(g,min_cut1,cap); |
202 | 207 |
|
203 | 208 |
check(preflow_test.flowValue() == min_cut_value && |
204 | 209 |
min_cut_value == 2*flow_value, |
205 | 210 |
"The max flow value or the min cut value is wrong."); |
206 | 211 |
|
207 | 212 |
preflow_test.startSecondPhase(); |
208 | 213 |
|
209 | 214 |
check(checkFlow(g, preflow_test.flowMap(), cap, s, t), |
210 | 215 |
"The flow is not feasible."); |
211 | 216 |
|
212 | 217 |
CutMap min_cut2(g); |
213 | 218 |
preflow_test.minCutMap(min_cut2); |
214 | 219 |
min_cut_value=cutValue(g,min_cut2,cap); |
215 | 220 |
|
216 | 221 |
check(preflow_test.flowValue() == min_cut_value && |
217 | 222 |
min_cut_value == 2*flow_value, |
218 | 223 |
"The max flow value or the three min cut values were not doubled"); |
219 | 224 |
|
220 | 225 |
|
221 | 226 |
preflow_test.flowMap(flow); |
222 | 227 |
|
223 | 228 |
NodeIt tmp1(g,s); |
224 | 229 |
++tmp1; |
225 | 230 |
if ( tmp1 != INVALID ) s=tmp1; |
226 | 231 |
|
227 | 232 |
NodeIt tmp2(g,t); |
228 | 233 |
++tmp2; |
229 | 234 |
if ( tmp2 != INVALID ) t=tmp2; |
230 | 235 |
|
231 | 236 |
preflow_test.source(s); |
232 | 237 |
preflow_test.target(t); |
233 | 238 |
|
234 | 239 |
preflow_test.run(); |
235 | 240 |
|
236 | 241 |
CutMap min_cut3(g); |
237 | 242 |
preflow_test.minCutMap(min_cut3); |
238 | 243 |
min_cut_value=cutValue(g,min_cut3,cap); |
239 | 244 |
|
240 | 245 |
|
241 | 246 |
check(preflow_test.flowValue() == min_cut_value, |
242 | 247 |
"The max flow value or the three min cut values are incorrect."); |
243 | 248 |
|
244 | 249 |
return 0; |
245 | 250 |
} |
0 comments (0 inline)