0
5
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-2010 |
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_CAPACITY_SCALING_H |
20 | 20 |
#define LEMON_CAPACITY_SCALING_H |
21 | 21 |
|
22 | 22 |
/// \ingroup min_cost_flow_algs |
23 | 23 |
/// |
24 | 24 |
/// \file |
25 | 25 |
/// \brief Capacity Scaling algorithm for finding a minimum cost flow. |
26 | 26 |
|
27 | 27 |
#include <vector> |
28 | 28 |
#include <limits> |
29 | 29 |
#include <lemon/core.h> |
30 | 30 |
#include <lemon/bin_heap.h> |
31 | 31 |
|
32 | 32 |
namespace lemon { |
33 | 33 |
|
34 | 34 |
/// \brief Default traits class of CapacityScaling algorithm. |
35 | 35 |
/// |
36 | 36 |
/// Default traits class of CapacityScaling algorithm. |
37 | 37 |
/// \tparam GR Digraph type. |
38 | 38 |
/// \tparam V The number type used for flow amounts, capacity bounds |
39 | 39 |
/// and supply values. By default it is \c int. |
40 | 40 |
/// \tparam C The number type used for costs and potentials. |
41 | 41 |
/// By default it is the same as \c V. |
42 | 42 |
template <typename GR, typename V = int, typename C = V> |
43 | 43 |
struct CapacityScalingDefaultTraits |
44 | 44 |
{ |
45 | 45 |
/// The type of the digraph |
46 | 46 |
typedef GR Digraph; |
47 | 47 |
/// The type of the flow amounts, capacity bounds and supply values |
48 | 48 |
typedef V Value; |
49 | 49 |
/// The type of the arc costs |
50 | 50 |
typedef C Cost; |
51 | 51 |
|
52 | 52 |
/// \brief The type of the heap used for internal Dijkstra computations. |
53 | 53 |
/// |
54 | 54 |
/// The type of the heap used for internal Dijkstra computations. |
55 | 55 |
/// It must conform to the \ref lemon::concepts::Heap "Heap" concept, |
56 | 56 |
/// its priority type must be \c Cost and its cross reference type |
57 | 57 |
/// must be \ref RangeMap "RangeMap<int>". |
58 | 58 |
typedef BinHeap<Cost, RangeMap<int> > Heap; |
59 | 59 |
}; |
60 | 60 |
|
61 | 61 |
/// \addtogroup min_cost_flow_algs |
62 | 62 |
/// @{ |
63 | 63 |
|
64 | 64 |
/// \brief Implementation of the Capacity Scaling algorithm for |
65 | 65 |
/// finding a \ref min_cost_flow "minimum cost flow". |
66 | 66 |
/// |
67 | 67 |
/// \ref CapacityScaling implements the capacity scaling version |
68 | 68 |
/// of the successive shortest path algorithm for finding a |
69 | 69 |
/// \ref min_cost_flow "minimum cost flow" \ref amo93networkflows, |
70 | 70 |
/// \ref edmondskarp72theoretical. It is an efficient dual |
71 | 71 |
/// solution method. |
72 | 72 |
/// |
73 | 73 |
/// Most of the parameters of the problem (except for the digraph) |
74 | 74 |
/// can be given using separate functions, and the algorithm can be |
75 | 75 |
/// executed using the \ref run() function. If some parameters are not |
76 | 76 |
/// specified, then default values will be used. |
77 | 77 |
/// |
78 | 78 |
/// \tparam GR The digraph type the algorithm runs on. |
79 | 79 |
/// \tparam V The number type used for flow amounts, capacity bounds |
80 | 80 |
/// and supply values in the algorithm. By default, it is \c int. |
81 | 81 |
/// \tparam C The number type used for costs and potentials in the |
82 | 82 |
/// algorithm. By default, it is the same as \c V. |
83 | 83 |
/// \tparam TR The traits class that defines various types used by the |
84 | 84 |
/// algorithm. By default, it is \ref CapacityScalingDefaultTraits |
85 | 85 |
/// "CapacityScalingDefaultTraits<GR, V, C>". |
86 | 86 |
/// In most cases, this parameter should not be set directly, |
87 | 87 |
/// consider to use the named template parameters instead. |
88 | 88 |
/// |
89 |
/// \warning Both |
|
89 |
/// \warning Both \c V and \c C must be signed number types. |
|
90 |
/// \warning All input data (capacities, supply values, and costs) must |
|
90 | 91 |
/// be integer. |
91 | 92 |
/// \warning This algorithm does not support negative costs for such |
92 | 93 |
/// arcs that have infinite upper bound. |
93 | 94 |
#ifdef DOXYGEN |
94 | 95 |
template <typename GR, typename V, typename C, typename TR> |
95 | 96 |
#else |
96 | 97 |
template < typename GR, typename V = int, typename C = V, |
97 | 98 |
typename TR = CapacityScalingDefaultTraits<GR, V, C> > |
98 | 99 |
#endif |
99 | 100 |
class CapacityScaling |
100 | 101 |
{ |
101 | 102 |
public: |
102 | 103 |
|
103 | 104 |
/// The type of the digraph |
104 | 105 |
typedef typename TR::Digraph Digraph; |
105 | 106 |
/// The type of the flow amounts, capacity bounds and supply values |
106 | 107 |
typedef typename TR::Value Value; |
107 | 108 |
/// The type of the arc costs |
108 | 109 |
typedef typename TR::Cost Cost; |
109 | 110 |
|
110 | 111 |
/// The type of the heap used for internal Dijkstra computations |
111 | 112 |
typedef typename TR::Heap Heap; |
112 | 113 |
|
113 | 114 |
/// The \ref CapacityScalingDefaultTraits "traits class" of the algorithm |
114 | 115 |
typedef TR Traits; |
115 | 116 |
|
116 | 117 |
public: |
117 | 118 |
|
118 | 119 |
/// \brief Problem type constants for the \c run() function. |
119 | 120 |
/// |
120 | 121 |
/// Enum type containing the problem type constants that can be |
121 | 122 |
/// returned by the \ref run() function of the algorithm. |
122 | 123 |
enum ProblemType { |
123 | 124 |
/// The problem has no feasible solution (flow). |
124 | 125 |
INFEASIBLE, |
125 | 126 |
/// The problem has optimal solution (i.e. it is feasible and |
126 | 127 |
/// bounded), and the algorithm has found optimal flow and node |
127 | 128 |
/// potentials (primal and dual solutions). |
128 | 129 |
OPTIMAL, |
129 | 130 |
/// The digraph contains an arc of negative cost and infinite |
130 | 131 |
/// upper bound. It means that the objective function is unbounded |
131 | 132 |
/// on that arc, however, note that it could actually be bounded |
132 | 133 |
/// over the feasible flows, but this algroithm cannot handle |
133 | 134 |
/// these cases. |
134 | 135 |
UNBOUNDED |
135 | 136 |
}; |
136 | 137 |
|
137 | 138 |
private: |
138 | 139 |
|
139 | 140 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
140 | 141 |
|
141 | 142 |
typedef std::vector<int> IntVector; |
142 | 143 |
typedef std::vector<Value> ValueVector; |
143 | 144 |
typedef std::vector<Cost> CostVector; |
144 | 145 |
typedef std::vector<char> BoolVector; |
145 | 146 |
// Note: vector<char> is used instead of vector<bool> for efficiency reasons |
146 | 147 |
|
147 | 148 |
private: |
148 | 149 |
|
149 | 150 |
// Data related to the underlying digraph |
150 | 151 |
const GR &_graph; |
151 | 152 |
int _node_num; |
152 | 153 |
int _arc_num; |
153 | 154 |
int _res_arc_num; |
154 | 155 |
int _root; |
155 | 156 |
|
156 | 157 |
// Parameters of the problem |
157 | 158 |
bool _have_lower; |
158 | 159 |
Value _sum_supply; |
159 | 160 |
|
160 | 161 |
// Data structures for storing the digraph |
161 | 162 |
IntNodeMap _node_id; |
162 | 163 |
IntArcMap _arc_idf; |
163 | 164 |
IntArcMap _arc_idb; |
164 | 165 |
IntVector _first_out; |
165 | 166 |
BoolVector _forward; |
166 | 167 |
IntVector _source; |
167 | 168 |
IntVector _target; |
168 | 169 |
IntVector _reverse; |
169 | 170 |
|
170 | 171 |
// Node and arc data |
171 | 172 |
ValueVector _lower; |
172 | 173 |
ValueVector _upper; |
173 | 174 |
CostVector _cost; |
174 | 175 |
ValueVector _supply; |
175 | 176 |
|
176 | 177 |
ValueVector _res_cap; |
177 | 178 |
CostVector _pi; |
178 | 179 |
ValueVector _excess; |
179 | 180 |
IntVector _excess_nodes; |
180 | 181 |
IntVector _deficit_nodes; |
181 | 182 |
|
182 | 183 |
Value _delta; |
183 | 184 |
int _factor; |
184 | 185 |
IntVector _pred; |
185 | 186 |
|
186 | 187 |
public: |
187 | 188 |
|
188 | 189 |
/// \brief Constant for infinite upper bounds (capacities). |
189 | 190 |
/// |
190 | 191 |
/// Constant for infinite upper bounds (capacities). |
191 | 192 |
/// It is \c std::numeric_limits<Value>::infinity() if available, |
192 | 193 |
/// \c std::numeric_limits<Value>::max() otherwise. |
193 | 194 |
const Value INF; |
194 | 195 |
|
195 | 196 |
private: |
196 | 197 |
|
197 | 198 |
// Special implementation of the Dijkstra algorithm for finding |
198 | 199 |
// shortest paths in the residual network of the digraph with |
199 | 200 |
// respect to the reduced arc costs and modifying the node |
200 | 201 |
// potentials according to the found distance labels. |
201 | 202 |
class ResidualDijkstra |
202 | 203 |
{ |
203 | 204 |
private: |
204 | 205 |
|
205 | 206 |
int _node_num; |
206 | 207 |
bool _geq; |
207 | 208 |
const IntVector &_first_out; |
208 | 209 |
const IntVector &_target; |
209 | 210 |
const CostVector &_cost; |
210 | 211 |
const ValueVector &_res_cap; |
211 | 212 |
const ValueVector &_excess; |
212 | 213 |
CostVector &_pi; |
213 | 214 |
IntVector &_pred; |
214 | 215 |
|
215 | 216 |
IntVector _proc_nodes; |
216 | 217 |
CostVector _dist; |
217 | 218 |
|
218 | 219 |
public: |
219 | 220 |
|
220 | 221 |
ResidualDijkstra(CapacityScaling& cs) : |
221 | 222 |
_node_num(cs._node_num), _geq(cs._sum_supply < 0), |
222 | 223 |
_first_out(cs._first_out), _target(cs._target), _cost(cs._cost), |
223 | 224 |
_res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi), |
224 | 225 |
_pred(cs._pred), _dist(cs._node_num) |
225 | 226 |
{} |
226 | 227 |
|
227 | 228 |
int run(int s, Value delta = 1) { |
228 | 229 |
RangeMap<int> heap_cross_ref(_node_num, Heap::PRE_HEAP); |
229 | 230 |
Heap heap(heap_cross_ref); |
230 | 231 |
heap.push(s, 0); |
231 | 232 |
_pred[s] = -1; |
232 | 233 |
_proc_nodes.clear(); |
233 | 234 |
|
234 | 235 |
// Process nodes |
235 | 236 |
while (!heap.empty() && _excess[heap.top()] > -delta) { |
236 | 237 |
int u = heap.top(), v; |
237 | 238 |
Cost d = heap.prio() + _pi[u], dn; |
238 | 239 |
_dist[u] = heap.prio(); |
239 | 240 |
_proc_nodes.push_back(u); |
240 | 241 |
heap.pop(); |
241 | 242 |
|
242 | 243 |
// Traverse outgoing residual arcs |
243 | 244 |
int last_out = _geq ? _first_out[u+1] : _first_out[u+1] - 1; |
244 | 245 |
for (int a = _first_out[u]; a != last_out; ++a) { |
245 | 246 |
if (_res_cap[a] < delta) continue; |
246 | 247 |
v = _target[a]; |
247 | 248 |
switch (heap.state(v)) { |
248 | 249 |
case Heap::PRE_HEAP: |
249 | 250 |
heap.push(v, d + _cost[a] - _pi[v]); |
250 | 251 |
_pred[v] = a; |
251 | 252 |
break; |
252 | 253 |
case Heap::IN_HEAP: |
253 | 254 |
dn = d + _cost[a] - _pi[v]; |
254 | 255 |
if (dn < heap[v]) { |
255 | 256 |
heap.decrease(v, dn); |
256 | 257 |
_pred[v] = a; |
257 | 258 |
} |
258 | 259 |
break; |
259 | 260 |
case Heap::POST_HEAP: |
260 | 261 |
break; |
261 | 262 |
} |
262 | 263 |
} |
263 | 264 |
} |
264 | 265 |
if (heap.empty()) return -1; |
265 | 266 |
|
266 | 267 |
// Update potentials of processed nodes |
267 | 268 |
int t = heap.top(); |
268 | 269 |
Cost dt = heap.prio(); |
269 | 270 |
for (int i = 0; i < int(_proc_nodes.size()); ++i) { |
270 | 271 |
_pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt; |
271 | 272 |
} |
272 | 273 |
|
273 | 274 |
return t; |
274 | 275 |
} |
275 | 276 |
|
276 | 277 |
}; //class ResidualDijkstra |
277 | 278 |
|
278 | 279 |
public: |
279 | 280 |
|
280 | 281 |
/// \name Named Template Parameters |
281 | 282 |
/// @{ |
282 | 283 |
|
283 | 284 |
template <typename T> |
284 | 285 |
struct SetHeapTraits : public Traits { |
285 | 286 |
typedef T Heap; |
286 | 287 |
}; |
287 | 288 |
|
288 | 289 |
/// \brief \ref named-templ-param "Named parameter" for setting |
289 | 290 |
/// \c Heap type. |
290 | 291 |
/// |
291 | 292 |
/// \ref named-templ-param "Named parameter" for setting \c Heap |
292 | 293 |
/// type, which is used for internal Dijkstra computations. |
293 | 294 |
/// It must conform to the \ref lemon::concepts::Heap "Heap" concept, |
294 | 295 |
/// its priority type must be \c Cost and its cross reference type |
295 | 296 |
/// must be \ref RangeMap "RangeMap<int>". |
296 | 297 |
template <typename T> |
297 | 298 |
struct SetHeap |
298 | 299 |
: public CapacityScaling<GR, V, C, SetHeapTraits<T> > { |
299 | 300 |
typedef CapacityScaling<GR, V, C, SetHeapTraits<T> > Create; |
300 | 301 |
}; |
301 | 302 |
|
302 | 303 |
/// @} |
303 | 304 |
|
304 | 305 |
protected: |
305 | 306 |
|
306 | 307 |
CapacityScaling() {} |
307 | 308 |
|
308 | 309 |
public: |
309 | 310 |
|
310 | 311 |
/// \brief Constructor. |
311 | 312 |
/// |
312 | 313 |
/// The constructor of the class. |
313 | 314 |
/// |
314 | 315 |
/// \param graph The digraph the algorithm runs on. |
315 | 316 |
CapacityScaling(const GR& graph) : |
316 | 317 |
_graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), |
317 | 318 |
INF(std::numeric_limits<Value>::has_infinity ? |
318 | 319 |
std::numeric_limits<Value>::infinity() : |
319 | 320 |
std::numeric_limits<Value>::max()) |
320 | 321 |
{ |
321 | 322 |
// Check the number types |
322 | 323 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
323 | 324 |
"The flow type of CapacityScaling must be signed"); |
324 | 325 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
325 | 326 |
"The cost type of CapacityScaling must be signed"); |
326 | 327 |
|
327 | 328 |
// Reset data structures |
328 | 329 |
reset(); |
329 | 330 |
} |
330 | 331 |
|
331 | 332 |
/// \name Parameters |
332 | 333 |
/// The parameters of the algorithm can be specified using these |
333 | 334 |
/// functions. |
334 | 335 |
|
335 | 336 |
/// @{ |
336 | 337 |
|
337 | 338 |
/// \brief Set the lower bounds on the arcs. |
338 | 339 |
/// |
339 | 340 |
/// This function sets the lower bounds on the arcs. |
340 | 341 |
/// If it is not used before calling \ref run(), the lower bounds |
341 | 342 |
/// will be set to zero on all arcs. |
342 | 343 |
/// |
343 | 344 |
/// \param map An arc map storing the lower bounds. |
344 | 345 |
/// Its \c Value type must be convertible to the \c Value type |
345 | 346 |
/// of the algorithm. |
346 | 347 |
/// |
347 | 348 |
/// \return <tt>(*this)</tt> |
348 | 349 |
template <typename LowerMap> |
349 | 350 |
CapacityScaling& lowerMap(const LowerMap& map) { |
350 | 351 |
_have_lower = true; |
351 | 352 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
352 | 353 |
_lower[_arc_idf[a]] = map[a]; |
353 | 354 |
_lower[_arc_idb[a]] = map[a]; |
354 | 355 |
} |
355 | 356 |
return *this; |
356 | 357 |
} |
357 | 358 |
|
358 | 359 |
/// \brief Set the upper bounds (capacities) on the arcs. |
359 | 360 |
/// |
360 | 361 |
/// This function sets the upper bounds (capacities) on the arcs. |
361 | 362 |
/// If it is not used before calling \ref run(), the upper bounds |
362 | 363 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be |
363 | 364 |
/// unbounded from above). |
364 | 365 |
/// |
365 | 366 |
/// \param map An arc map storing the upper bounds. |
366 | 367 |
/// Its \c Value type must be convertible to the \c Value type |
367 | 368 |
/// of the algorithm. |
368 | 369 |
/// |
369 | 370 |
/// \return <tt>(*this)</tt> |
370 | 371 |
template<typename UpperMap> |
371 | 372 |
CapacityScaling& upperMap(const UpperMap& map) { |
372 | 373 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
373 | 374 |
_upper[_arc_idf[a]] = map[a]; |
374 | 375 |
} |
375 | 376 |
return *this; |
376 | 377 |
} |
377 | 378 |
|
378 | 379 |
/// \brief Set the costs of the arcs. |
379 | 380 |
/// |
380 | 381 |
/// This function sets the costs of the arcs. |
381 | 382 |
/// If it is not used before calling \ref run(), the costs |
382 | 383 |
/// will be set to \c 1 on all arcs. |
383 | 384 |
/// |
384 | 385 |
/// \param map An arc map storing the costs. |
385 | 386 |
/// Its \c Value type must be convertible to the \c Cost type |
386 | 387 |
/// of the algorithm. |
387 | 388 |
/// |
388 | 389 |
/// \return <tt>(*this)</tt> |
389 | 390 |
template<typename CostMap> |
390 | 391 |
CapacityScaling& costMap(const CostMap& map) { |
391 | 392 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
392 | 393 |
_cost[_arc_idf[a]] = map[a]; |
393 | 394 |
_cost[_arc_idb[a]] = -map[a]; |
394 | 395 |
} |
395 | 396 |
return *this; |
396 | 397 |
} |
397 | 398 |
|
398 | 399 |
/// \brief Set the supply values of the nodes. |
399 | 400 |
/// |
400 | 401 |
/// This function sets the supply values of the nodes. |
401 | 402 |
/// If neither this function nor \ref stSupply() is used before |
402 | 403 |
/// calling \ref run(), the supply of each node will be set to zero. |
403 | 404 |
/// |
404 | 405 |
/// \param map A node map storing the supply values. |
405 | 406 |
/// Its \c Value type must be convertible to the \c Value type |
406 | 407 |
/// of the algorithm. |
407 | 408 |
/// |
408 | 409 |
/// \return <tt>(*this)</tt> |
409 | 410 |
template<typename SupplyMap> |
410 | 411 |
CapacityScaling& supplyMap(const SupplyMap& map) { |
411 | 412 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
412 | 413 |
_supply[_node_id[n]] = map[n]; |
413 | 414 |
} |
414 | 415 |
return *this; |
415 | 416 |
} |
416 | 417 |
|
417 | 418 |
/// \brief Set single source and target nodes and a supply value. |
418 | 419 |
/// |
419 | 420 |
/// This function sets a single source node and a single target node |
420 | 421 |
/// and the required flow value. |
421 | 422 |
/// If neither this function nor \ref supplyMap() is used before |
422 | 423 |
/// calling \ref run(), the supply of each node will be set to zero. |
423 | 424 |
/// |
424 | 425 |
/// Using this function has the same effect as using \ref supplyMap() |
425 | 426 |
/// with such a map in which \c k is assigned to \c s, \c -k is |
426 | 427 |
/// assigned to \c t and all other nodes have zero supply value. |
427 | 428 |
/// |
428 | 429 |
/// \param s The source node. |
429 | 430 |
/// \param t The target node. |
430 | 431 |
/// \param k The required amount of flow from node \c s to node \c t |
431 | 432 |
/// (i.e. the supply of \c s and the demand of \c t). |
432 | 433 |
/// |
433 | 434 |
/// \return <tt>(*this)</tt> |
434 | 435 |
CapacityScaling& stSupply(const Node& s, const Node& t, Value k) { |
435 | 436 |
for (int i = 0; i != _node_num; ++i) { |
436 | 437 |
_supply[i] = 0; |
437 | 438 |
} |
438 | 439 |
_supply[_node_id[s]] = k; |
439 | 440 |
_supply[_node_id[t]] = -k; |
440 | 441 |
return *this; |
441 | 442 |
} |
442 | 443 |
|
443 | 444 |
/// @} |
444 | 445 |
|
445 | 446 |
/// \name Execution control |
446 | 447 |
/// The algorithm can be executed using \ref run(). |
447 | 448 |
|
448 | 449 |
/// @{ |
449 | 450 |
|
450 | 451 |
/// \brief Run the algorithm. |
451 | 452 |
/// |
452 | 453 |
/// This function runs the algorithm. |
453 | 454 |
/// The paramters can be specified using functions \ref lowerMap(), |
454 | 455 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). |
455 | 456 |
/// For example, |
456 | 457 |
/// \code |
457 | 458 |
/// CapacityScaling<ListDigraph> cs(graph); |
458 | 459 |
/// cs.lowerMap(lower).upperMap(upper).costMap(cost) |
459 | 460 |
/// .supplyMap(sup).run(); |
460 | 461 |
/// \endcode |
461 | 462 |
/// |
462 | 463 |
/// This function can be called more than once. All the given parameters |
463 | 464 |
/// are kept for the next call, unless \ref resetParams() or \ref reset() |
464 | 465 |
/// is used, thus only the modified parameters have to be set again. |
465 | 466 |
/// If the underlying digraph was also modified after the construction |
466 | 467 |
/// of the class (or the last \ref reset() call), then the \ref reset() |
467 | 468 |
/// function must be called. |
468 | 469 |
/// |
469 | 470 |
/// \param factor The capacity scaling factor. It must be larger than |
470 | 471 |
/// one to use scaling. If it is less or equal to one, then scaling |
471 | 472 |
/// will be disabled. |
472 | 473 |
/// |
473 | 474 |
/// \return \c INFEASIBLE if no feasible flow exists, |
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-2010 |
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_COST_SCALING_H |
20 | 20 |
#define LEMON_COST_SCALING_H |
21 | 21 |
|
22 | 22 |
/// \ingroup min_cost_flow_algs |
23 | 23 |
/// \file |
24 | 24 |
/// \brief Cost scaling algorithm for finding a minimum cost flow. |
25 | 25 |
|
26 | 26 |
#include <vector> |
27 | 27 |
#include <deque> |
28 | 28 |
#include <limits> |
29 | 29 |
|
30 | 30 |
#include <lemon/core.h> |
31 | 31 |
#include <lemon/maps.h> |
32 | 32 |
#include <lemon/math.h> |
33 | 33 |
#include <lemon/static_graph.h> |
34 | 34 |
#include <lemon/circulation.h> |
35 | 35 |
#include <lemon/bellman_ford.h> |
36 | 36 |
|
37 | 37 |
namespace lemon { |
38 | 38 |
|
39 | 39 |
/// \brief Default traits class of CostScaling algorithm. |
40 | 40 |
/// |
41 | 41 |
/// Default traits class of CostScaling algorithm. |
42 | 42 |
/// \tparam GR Digraph type. |
43 | 43 |
/// \tparam V The number type used for flow amounts, capacity bounds |
44 | 44 |
/// and supply values. By default it is \c int. |
45 | 45 |
/// \tparam C The number type used for costs and potentials. |
46 | 46 |
/// By default it is the same as \c V. |
47 | 47 |
#ifdef DOXYGEN |
48 | 48 |
template <typename GR, typename V = int, typename C = V> |
49 | 49 |
#else |
50 | 50 |
template < typename GR, typename V = int, typename C = V, |
51 | 51 |
bool integer = std::numeric_limits<C>::is_integer > |
52 | 52 |
#endif |
53 | 53 |
struct CostScalingDefaultTraits |
54 | 54 |
{ |
55 | 55 |
/// The type of the digraph |
56 | 56 |
typedef GR Digraph; |
57 | 57 |
/// The type of the flow amounts, capacity bounds and supply values |
58 | 58 |
typedef V Value; |
59 | 59 |
/// The type of the arc costs |
60 | 60 |
typedef C Cost; |
61 | 61 |
|
62 | 62 |
/// \brief The large cost type used for internal computations |
63 | 63 |
/// |
64 | 64 |
/// The large cost type used for internal computations. |
65 | 65 |
/// It is \c long \c long if the \c Cost type is integer, |
66 | 66 |
/// otherwise it is \c double. |
67 | 67 |
/// \c Cost must be convertible to \c LargeCost. |
68 | 68 |
typedef double LargeCost; |
69 | 69 |
}; |
70 | 70 |
|
71 | 71 |
// Default traits class for integer cost types |
72 | 72 |
template <typename GR, typename V, typename C> |
73 | 73 |
struct CostScalingDefaultTraits<GR, V, C, true> |
74 | 74 |
{ |
75 | 75 |
typedef GR Digraph; |
76 | 76 |
typedef V Value; |
77 | 77 |
typedef C Cost; |
78 | 78 |
#ifdef LEMON_HAVE_LONG_LONG |
79 | 79 |
typedef long long LargeCost; |
80 | 80 |
#else |
81 | 81 |
typedef long LargeCost; |
82 | 82 |
#endif |
83 | 83 |
}; |
84 | 84 |
|
85 | 85 |
|
86 | 86 |
/// \addtogroup min_cost_flow_algs |
87 | 87 |
/// @{ |
88 | 88 |
|
89 | 89 |
/// \brief Implementation of the Cost Scaling algorithm for |
90 | 90 |
/// finding a \ref min_cost_flow "minimum cost flow". |
91 | 91 |
/// |
92 | 92 |
/// \ref CostScaling implements a cost scaling algorithm that performs |
93 | 93 |
/// push/augment and relabel operations for finding a \ref min_cost_flow |
94 | 94 |
/// "minimum cost flow" \ref amo93networkflows, \ref goldberg90approximation, |
95 | 95 |
/// \ref goldberg97efficient, \ref bunnagel98efficient. |
96 | 96 |
/// It is a highly efficient primal-dual solution method, which |
97 | 97 |
/// can be viewed as the generalization of the \ref Preflow |
98 | 98 |
/// "preflow push-relabel" algorithm for the maximum flow problem. |
99 | 99 |
/// |
100 | 100 |
/// Most of the parameters of the problem (except for the digraph) |
101 | 101 |
/// can be given using separate functions, and the algorithm can be |
102 | 102 |
/// executed using the \ref run() function. If some parameters are not |
103 | 103 |
/// specified, then default values will be used. |
104 | 104 |
/// |
105 | 105 |
/// \tparam GR The digraph type the algorithm runs on. |
106 | 106 |
/// \tparam V The number type used for flow amounts, capacity bounds |
107 | 107 |
/// and supply values in the algorithm. By default, it is \c int. |
108 | 108 |
/// \tparam C The number type used for costs and potentials in the |
109 | 109 |
/// algorithm. By default, it is the same as \c V. |
110 | 110 |
/// \tparam TR The traits class that defines various types used by the |
111 | 111 |
/// algorithm. By default, it is \ref CostScalingDefaultTraits |
112 | 112 |
/// "CostScalingDefaultTraits<GR, V, C>". |
113 | 113 |
/// In most cases, this parameter should not be set directly, |
114 | 114 |
/// consider to use the named template parameters instead. |
115 | 115 |
/// |
116 |
/// \warning Both |
|
116 |
/// \warning Both \c V and \c C must be signed number types. |
|
117 |
/// \warning All input data (capacities, supply values, and costs) must |
|
117 | 118 |
/// be integer. |
118 | 119 |
/// \warning This algorithm does not support negative costs for such |
119 | 120 |
/// arcs that have infinite upper bound. |
120 | 121 |
/// |
121 | 122 |
/// \note %CostScaling provides three different internal methods, |
122 | 123 |
/// from which the most efficient one is used by default. |
123 | 124 |
/// For more information, see \ref Method. |
124 | 125 |
#ifdef DOXYGEN |
125 | 126 |
template <typename GR, typename V, typename C, typename TR> |
126 | 127 |
#else |
127 | 128 |
template < typename GR, typename V = int, typename C = V, |
128 | 129 |
typename TR = CostScalingDefaultTraits<GR, V, C> > |
129 | 130 |
#endif |
130 | 131 |
class CostScaling |
131 | 132 |
{ |
132 | 133 |
public: |
133 | 134 |
|
134 | 135 |
/// The type of the digraph |
135 | 136 |
typedef typename TR::Digraph Digraph; |
136 | 137 |
/// The type of the flow amounts, capacity bounds and supply values |
137 | 138 |
typedef typename TR::Value Value; |
138 | 139 |
/// The type of the arc costs |
139 | 140 |
typedef typename TR::Cost Cost; |
140 | 141 |
|
141 | 142 |
/// \brief The large cost type |
142 | 143 |
/// |
143 | 144 |
/// The large cost type used for internal computations. |
144 | 145 |
/// By default, it is \c long \c long if the \c Cost type is integer, |
145 | 146 |
/// otherwise it is \c double. |
146 | 147 |
typedef typename TR::LargeCost LargeCost; |
147 | 148 |
|
148 | 149 |
/// The \ref CostScalingDefaultTraits "traits class" of the algorithm |
149 | 150 |
typedef TR Traits; |
150 | 151 |
|
151 | 152 |
public: |
152 | 153 |
|
153 | 154 |
/// \brief Problem type constants for the \c run() function. |
154 | 155 |
/// |
155 | 156 |
/// Enum type containing the problem type constants that can be |
156 | 157 |
/// returned by the \ref run() function of the algorithm. |
157 | 158 |
enum ProblemType { |
158 | 159 |
/// The problem has no feasible solution (flow). |
159 | 160 |
INFEASIBLE, |
160 | 161 |
/// The problem has optimal solution (i.e. it is feasible and |
161 | 162 |
/// bounded), and the algorithm has found optimal flow and node |
162 | 163 |
/// potentials (primal and dual solutions). |
163 | 164 |
OPTIMAL, |
164 | 165 |
/// The digraph contains an arc of negative cost and infinite |
165 | 166 |
/// upper bound. It means that the objective function is unbounded |
166 | 167 |
/// on that arc, however, note that it could actually be bounded |
167 | 168 |
/// over the feasible flows, but this algroithm cannot handle |
168 | 169 |
/// these cases. |
169 | 170 |
UNBOUNDED |
170 | 171 |
}; |
171 | 172 |
|
172 | 173 |
/// \brief Constants for selecting the internal method. |
173 | 174 |
/// |
174 | 175 |
/// Enum type containing constants for selecting the internal method |
175 | 176 |
/// for the \ref run() function. |
176 | 177 |
/// |
177 | 178 |
/// \ref CostScaling provides three internal methods that differ mainly |
178 | 179 |
/// in their base operations, which are used in conjunction with the |
179 | 180 |
/// relabel operation. |
180 | 181 |
/// By default, the so called \ref PARTIAL_AUGMENT |
181 | 182 |
/// "Partial Augment-Relabel" method is used, which proved to be |
182 | 183 |
/// the most efficient and the most robust on various test inputs. |
183 | 184 |
/// However, the other methods can be selected using the \ref run() |
184 | 185 |
/// function with the proper parameter. |
185 | 186 |
enum Method { |
186 | 187 |
/// Local push operations are used, i.e. flow is moved only on one |
187 | 188 |
/// admissible arc at once. |
188 | 189 |
PUSH, |
189 | 190 |
/// Augment operations are used, i.e. flow is moved on admissible |
190 | 191 |
/// paths from a node with excess to a node with deficit. |
191 | 192 |
AUGMENT, |
192 | 193 |
/// Partial augment operations are used, i.e. flow is moved on |
193 | 194 |
/// admissible paths started from a node with excess, but the |
194 | 195 |
/// lengths of these paths are limited. This method can be viewed |
195 | 196 |
/// as a combined version of the previous two operations. |
196 | 197 |
PARTIAL_AUGMENT |
197 | 198 |
}; |
198 | 199 |
|
199 | 200 |
private: |
200 | 201 |
|
201 | 202 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
202 | 203 |
|
203 | 204 |
typedef std::vector<int> IntVector; |
204 | 205 |
typedef std::vector<Value> ValueVector; |
205 | 206 |
typedef std::vector<Cost> CostVector; |
206 | 207 |
typedef std::vector<LargeCost> LargeCostVector; |
207 | 208 |
typedef std::vector<char> BoolVector; |
208 | 209 |
// Note: vector<char> is used instead of vector<bool> for efficiency reasons |
209 | 210 |
|
210 | 211 |
private: |
211 | 212 |
|
212 | 213 |
template <typename KT, typename VT> |
213 | 214 |
class StaticVectorMap { |
214 | 215 |
public: |
215 | 216 |
typedef KT Key; |
216 | 217 |
typedef VT Value; |
217 | 218 |
|
218 | 219 |
StaticVectorMap(std::vector<Value>& v) : _v(v) {} |
219 | 220 |
|
220 | 221 |
const Value& operator[](const Key& key) const { |
221 | 222 |
return _v[StaticDigraph::id(key)]; |
222 | 223 |
} |
223 | 224 |
|
224 | 225 |
Value& operator[](const Key& key) { |
225 | 226 |
return _v[StaticDigraph::id(key)]; |
226 | 227 |
} |
227 | 228 |
|
228 | 229 |
void set(const Key& key, const Value& val) { |
229 | 230 |
_v[StaticDigraph::id(key)] = val; |
230 | 231 |
} |
231 | 232 |
|
232 | 233 |
private: |
233 | 234 |
std::vector<Value>& _v; |
234 | 235 |
}; |
235 | 236 |
|
236 | 237 |
typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap; |
237 | 238 |
typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap; |
238 | 239 |
|
239 | 240 |
private: |
240 | 241 |
|
241 | 242 |
// Data related to the underlying digraph |
242 | 243 |
const GR &_graph; |
243 | 244 |
int _node_num; |
244 | 245 |
int _arc_num; |
245 | 246 |
int _res_node_num; |
246 | 247 |
int _res_arc_num; |
247 | 248 |
int _root; |
248 | 249 |
|
249 | 250 |
// Parameters of the problem |
250 | 251 |
bool _have_lower; |
251 | 252 |
Value _sum_supply; |
252 | 253 |
int _sup_node_num; |
253 | 254 |
|
254 | 255 |
// Data structures for storing the digraph |
255 | 256 |
IntNodeMap _node_id; |
256 | 257 |
IntArcMap _arc_idf; |
257 | 258 |
IntArcMap _arc_idb; |
258 | 259 |
IntVector _first_out; |
259 | 260 |
BoolVector _forward; |
260 | 261 |
IntVector _source; |
261 | 262 |
IntVector _target; |
262 | 263 |
IntVector _reverse; |
263 | 264 |
|
264 | 265 |
// Node and arc data |
265 | 266 |
ValueVector _lower; |
266 | 267 |
ValueVector _upper; |
267 | 268 |
CostVector _scost; |
268 | 269 |
ValueVector _supply; |
269 | 270 |
|
270 | 271 |
ValueVector _res_cap; |
271 | 272 |
LargeCostVector _cost; |
272 | 273 |
LargeCostVector _pi; |
273 | 274 |
ValueVector _excess; |
274 | 275 |
IntVector _next_out; |
275 | 276 |
std::deque<int> _active_nodes; |
276 | 277 |
|
277 | 278 |
// Data for scaling |
278 | 279 |
LargeCost _epsilon; |
279 | 280 |
int _alpha; |
280 | 281 |
|
281 | 282 |
IntVector _buckets; |
282 | 283 |
IntVector _bucket_next; |
283 | 284 |
IntVector _bucket_prev; |
284 | 285 |
IntVector _rank; |
285 | 286 |
int _max_rank; |
286 | 287 |
|
287 | 288 |
// Data for a StaticDigraph structure |
288 | 289 |
typedef std::pair<int, int> IntPair; |
289 | 290 |
StaticDigraph _sgr; |
290 | 291 |
std::vector<IntPair> _arc_vec; |
291 | 292 |
std::vector<LargeCost> _cost_vec; |
292 | 293 |
LargeCostArcMap _cost_map; |
293 | 294 |
LargeCostNodeMap _pi_map; |
294 | 295 |
|
295 | 296 |
public: |
296 | 297 |
|
297 | 298 |
/// \brief Constant for infinite upper bounds (capacities). |
298 | 299 |
/// |
299 | 300 |
/// Constant for infinite upper bounds (capacities). |
300 | 301 |
/// It is \c std::numeric_limits<Value>::infinity() if available, |
301 | 302 |
/// \c std::numeric_limits<Value>::max() otherwise. |
302 | 303 |
const Value INF; |
303 | 304 |
|
304 | 305 |
public: |
305 | 306 |
|
306 | 307 |
/// \name Named Template Parameters |
307 | 308 |
/// @{ |
308 | 309 |
|
309 | 310 |
template <typename T> |
310 | 311 |
struct SetLargeCostTraits : public Traits { |
311 | 312 |
typedef T LargeCost; |
312 | 313 |
}; |
313 | 314 |
|
314 | 315 |
/// \brief \ref named-templ-param "Named parameter" for setting |
315 | 316 |
/// \c LargeCost type. |
316 | 317 |
/// |
317 | 318 |
/// \ref named-templ-param "Named parameter" for setting \c LargeCost |
318 | 319 |
/// type, which is used for internal computations in the algorithm. |
319 | 320 |
/// \c Cost must be convertible to \c LargeCost. |
320 | 321 |
template <typename T> |
321 | 322 |
struct SetLargeCost |
322 | 323 |
: public CostScaling<GR, V, C, SetLargeCostTraits<T> > { |
323 | 324 |
typedef CostScaling<GR, V, C, SetLargeCostTraits<T> > Create; |
324 | 325 |
}; |
325 | 326 |
|
326 | 327 |
/// @} |
327 | 328 |
|
328 | 329 |
protected: |
329 | 330 |
|
330 | 331 |
CostScaling() {} |
331 | 332 |
|
332 | 333 |
public: |
333 | 334 |
|
334 | 335 |
/// \brief Constructor. |
335 | 336 |
/// |
336 | 337 |
/// The constructor of the class. |
337 | 338 |
/// |
338 | 339 |
/// \param graph The digraph the algorithm runs on. |
339 | 340 |
CostScaling(const GR& graph) : |
340 | 341 |
_graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), |
341 | 342 |
_cost_map(_cost_vec), _pi_map(_pi), |
342 | 343 |
INF(std::numeric_limits<Value>::has_infinity ? |
343 | 344 |
std::numeric_limits<Value>::infinity() : |
344 | 345 |
std::numeric_limits<Value>::max()) |
345 | 346 |
{ |
346 | 347 |
// Check the number types |
347 | 348 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
348 | 349 |
"The flow type of CostScaling must be signed"); |
349 | 350 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
350 | 351 |
"The cost type of CostScaling must be signed"); |
351 | 352 |
|
352 | 353 |
// Reset data structures |
353 | 354 |
reset(); |
354 | 355 |
} |
355 | 356 |
|
356 | 357 |
/// \name Parameters |
357 | 358 |
/// The parameters of the algorithm can be specified using these |
358 | 359 |
/// functions. |
359 | 360 |
|
360 | 361 |
/// @{ |
361 | 362 |
|
362 | 363 |
/// \brief Set the lower bounds on the arcs. |
363 | 364 |
/// |
364 | 365 |
/// This function sets the lower bounds on the arcs. |
365 | 366 |
/// If it is not used before calling \ref run(), the lower bounds |
366 | 367 |
/// will be set to zero on all arcs. |
367 | 368 |
/// |
368 | 369 |
/// \param map An arc map storing the lower bounds. |
369 | 370 |
/// Its \c Value type must be convertible to the \c Value type |
370 | 371 |
/// of the algorithm. |
371 | 372 |
/// |
372 | 373 |
/// \return <tt>(*this)</tt> |
373 | 374 |
template <typename LowerMap> |
374 | 375 |
CostScaling& lowerMap(const LowerMap& map) { |
375 | 376 |
_have_lower = true; |
376 | 377 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
377 | 378 |
_lower[_arc_idf[a]] = map[a]; |
378 | 379 |
_lower[_arc_idb[a]] = map[a]; |
379 | 380 |
} |
380 | 381 |
return *this; |
381 | 382 |
} |
382 | 383 |
|
383 | 384 |
/// \brief Set the upper bounds (capacities) on the arcs. |
384 | 385 |
/// |
385 | 386 |
/// This function sets the upper bounds (capacities) on the arcs. |
386 | 387 |
/// If it is not used before calling \ref run(), the upper bounds |
387 | 388 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be |
388 | 389 |
/// unbounded from above). |
389 | 390 |
/// |
390 | 391 |
/// \param map An arc map storing the upper bounds. |
391 | 392 |
/// Its \c Value type must be convertible to the \c Value type |
392 | 393 |
/// of the algorithm. |
393 | 394 |
/// |
394 | 395 |
/// \return <tt>(*this)</tt> |
395 | 396 |
template<typename UpperMap> |
396 | 397 |
CostScaling& upperMap(const UpperMap& map) { |
397 | 398 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
398 | 399 |
_upper[_arc_idf[a]] = map[a]; |
399 | 400 |
} |
400 | 401 |
return *this; |
401 | 402 |
} |
402 | 403 |
|
403 | 404 |
/// \brief Set the costs of the arcs. |
404 | 405 |
/// |
405 | 406 |
/// This function sets the costs of the arcs. |
406 | 407 |
/// If it is not used before calling \ref run(), the costs |
407 | 408 |
/// will be set to \c 1 on all arcs. |
408 | 409 |
/// |
409 | 410 |
/// \param map An arc map storing the costs. |
410 | 411 |
/// Its \c Value type must be convertible to the \c Cost type |
411 | 412 |
/// of the algorithm. |
412 | 413 |
/// |
413 | 414 |
/// \return <tt>(*this)</tt> |
414 | 415 |
template<typename CostMap> |
415 | 416 |
CostScaling& costMap(const CostMap& map) { |
416 | 417 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
417 | 418 |
_scost[_arc_idf[a]] = map[a]; |
418 | 419 |
_scost[_arc_idb[a]] = -map[a]; |
419 | 420 |
} |
420 | 421 |
return *this; |
421 | 422 |
} |
422 | 423 |
|
423 | 424 |
/// \brief Set the supply values of the nodes. |
424 | 425 |
/// |
425 | 426 |
/// This function sets the supply values of the nodes. |
426 | 427 |
/// If neither this function nor \ref stSupply() is used before |
427 | 428 |
/// calling \ref run(), the supply of each node will be set to zero. |
428 | 429 |
/// |
429 | 430 |
/// \param map A node map storing the supply values. |
430 | 431 |
/// Its \c Value type must be convertible to the \c Value type |
431 | 432 |
/// of the algorithm. |
432 | 433 |
/// |
433 | 434 |
/// \return <tt>(*this)</tt> |
434 | 435 |
template<typename SupplyMap> |
435 | 436 |
CostScaling& supplyMap(const SupplyMap& map) { |
436 | 437 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
437 | 438 |
_supply[_node_id[n]] = map[n]; |
438 | 439 |
} |
439 | 440 |
return *this; |
440 | 441 |
} |
441 | 442 |
|
442 | 443 |
/// \brief Set single source and target nodes and a supply value. |
443 | 444 |
/// |
444 | 445 |
/// This function sets a single source node and a single target node |
445 | 446 |
/// and the required flow value. |
446 | 447 |
/// If neither this function nor \ref supplyMap() is used before |
447 | 448 |
/// calling \ref run(), the supply of each node will be set to zero. |
448 | 449 |
/// |
449 | 450 |
/// Using this function has the same effect as using \ref supplyMap() |
450 | 451 |
/// with such a map in which \c k is assigned to \c s, \c -k is |
451 | 452 |
/// assigned to \c t and all other nodes have zero supply value. |
452 | 453 |
/// |
453 | 454 |
/// \param s The source node. |
454 | 455 |
/// \param t The target node. |
455 | 456 |
/// \param k The required amount of flow from node \c s to node \c t |
456 | 457 |
/// (i.e. the supply of \c s and the demand of \c t). |
457 | 458 |
/// |
458 | 459 |
/// \return <tt>(*this)</tt> |
459 | 460 |
CostScaling& stSupply(const Node& s, const Node& t, Value k) { |
460 | 461 |
for (int i = 0; i != _res_node_num; ++i) { |
461 | 462 |
_supply[i] = 0; |
462 | 463 |
} |
463 | 464 |
_supply[_node_id[s]] = k; |
464 | 465 |
_supply[_node_id[t]] = -k; |
465 | 466 |
return *this; |
466 | 467 |
} |
467 | 468 |
|
468 | 469 |
/// @} |
469 | 470 |
|
470 | 471 |
/// \name Execution control |
471 | 472 |
/// The algorithm can be executed using \ref run(). |
472 | 473 |
|
473 | 474 |
/// @{ |
474 | 475 |
|
475 | 476 |
/// \brief Run the algorithm. |
476 | 477 |
/// |
477 | 478 |
/// This function runs the algorithm. |
478 | 479 |
/// The paramters can be specified using functions \ref lowerMap(), |
479 | 480 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). |
480 | 481 |
/// For example, |
481 | 482 |
/// \code |
482 | 483 |
/// CostScaling<ListDigraph> cs(graph); |
483 | 484 |
/// cs.lowerMap(lower).upperMap(upper).costMap(cost) |
484 | 485 |
/// .supplyMap(sup).run(); |
485 | 486 |
/// \endcode |
486 | 487 |
/// |
487 | 488 |
/// This function can be called more than once. All the given parameters |
488 | 489 |
/// are kept for the next call, unless \ref resetParams() or \ref reset() |
489 | 490 |
/// is used, thus only the modified parameters have to be set again. |
490 | 491 |
/// If the underlying digraph was also modified after the construction |
491 | 492 |
/// of the class (or the last \ref reset() call), then the \ref reset() |
492 | 493 |
/// function must be called. |
493 | 494 |
/// |
494 | 495 |
/// \param method The internal method that will be used in the |
495 | 496 |
/// algorithm. For more information, see \ref Method. |
496 | 497 |
/// \param factor The cost scaling factor. It must be larger than one. |
497 | 498 |
/// |
498 | 499 |
/// \return \c INFEASIBLE if no feasible flow exists, |
499 | 500 |
/// \n \c OPTIMAL if the problem has optimal solution |
500 | 501 |
/// (i.e. it is feasible and bounded), and the algorithm has found |
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-2010 |
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_CYCLE_CANCELING_H |
20 | 20 |
#define LEMON_CYCLE_CANCELING_H |
21 | 21 |
|
22 | 22 |
/// \ingroup min_cost_flow_algs |
23 | 23 |
/// \file |
24 | 24 |
/// \brief Cycle-canceling algorithms for finding a minimum cost flow. |
25 | 25 |
|
26 | 26 |
#include <vector> |
27 | 27 |
#include <limits> |
28 | 28 |
|
29 | 29 |
#include <lemon/core.h> |
30 | 30 |
#include <lemon/maps.h> |
31 | 31 |
#include <lemon/path.h> |
32 | 32 |
#include <lemon/math.h> |
33 | 33 |
#include <lemon/static_graph.h> |
34 | 34 |
#include <lemon/adaptors.h> |
35 | 35 |
#include <lemon/circulation.h> |
36 | 36 |
#include <lemon/bellman_ford.h> |
37 | 37 |
#include <lemon/howard_mmc.h> |
38 | 38 |
|
39 | 39 |
namespace lemon { |
40 | 40 |
|
41 | 41 |
/// \addtogroup min_cost_flow_algs |
42 | 42 |
/// @{ |
43 | 43 |
|
44 | 44 |
/// \brief Implementation of cycle-canceling algorithms for |
45 | 45 |
/// finding a \ref min_cost_flow "minimum cost flow". |
46 | 46 |
/// |
47 | 47 |
/// \ref CycleCanceling implements three different cycle-canceling |
48 | 48 |
/// algorithms for finding a \ref min_cost_flow "minimum cost flow" |
49 | 49 |
/// \ref amo93networkflows, \ref klein67primal, |
50 | 50 |
/// \ref goldberg89cyclecanceling. |
51 | 51 |
/// The most efficent one (both theoretically and practically) |
52 | 52 |
/// is the \ref CANCEL_AND_TIGHTEN "Cancel and Tighten" algorithm, |
53 | 53 |
/// thus it is the default method. |
54 | 54 |
/// It is strongly polynomial, but in practice, it is typically much |
55 | 55 |
/// slower than the scaling algorithms and NetworkSimplex. |
56 | 56 |
/// |
57 | 57 |
/// Most of the parameters of the problem (except for the digraph) |
58 | 58 |
/// can be given using separate functions, and the algorithm can be |
59 | 59 |
/// executed using the \ref run() function. If some parameters are not |
60 | 60 |
/// specified, then default values will be used. |
61 | 61 |
/// |
62 | 62 |
/// \tparam GR The digraph type the algorithm runs on. |
63 | 63 |
/// \tparam V The number type used for flow amounts, capacity bounds |
64 | 64 |
/// and supply values in the algorithm. By default, it is \c int. |
65 | 65 |
/// \tparam C The number type used for costs and potentials in the |
66 | 66 |
/// algorithm. By default, it is the same as \c V. |
67 | 67 |
/// |
68 |
/// \warning Both |
|
68 |
/// \warning Both \c V and \c C must be signed number types. |
|
69 |
/// \warning All input data (capacities, supply values, and costs) must |
|
69 | 70 |
/// be integer. |
70 | 71 |
/// \warning This algorithm does not support negative costs for such |
71 | 72 |
/// arcs that have infinite upper bound. |
72 | 73 |
/// |
73 | 74 |
/// \note For more information about the three available methods, |
74 | 75 |
/// see \ref Method. |
75 | 76 |
#ifdef DOXYGEN |
76 | 77 |
template <typename GR, typename V, typename C> |
77 | 78 |
#else |
78 | 79 |
template <typename GR, typename V = int, typename C = V> |
79 | 80 |
#endif |
80 | 81 |
class CycleCanceling |
81 | 82 |
{ |
82 | 83 |
public: |
83 | 84 |
|
84 | 85 |
/// The type of the digraph |
85 | 86 |
typedef GR Digraph; |
86 | 87 |
/// The type of the flow amounts, capacity bounds and supply values |
87 | 88 |
typedef V Value; |
88 | 89 |
/// The type of the arc costs |
89 | 90 |
typedef C Cost; |
90 | 91 |
|
91 | 92 |
public: |
92 | 93 |
|
93 | 94 |
/// \brief Problem type constants for the \c run() function. |
94 | 95 |
/// |
95 | 96 |
/// Enum type containing the problem type constants that can be |
96 | 97 |
/// returned by the \ref run() function of the algorithm. |
97 | 98 |
enum ProblemType { |
98 | 99 |
/// The problem has no feasible solution (flow). |
99 | 100 |
INFEASIBLE, |
100 | 101 |
/// The problem has optimal solution (i.e. it is feasible and |
101 | 102 |
/// bounded), and the algorithm has found optimal flow and node |
102 | 103 |
/// potentials (primal and dual solutions). |
103 | 104 |
OPTIMAL, |
104 | 105 |
/// The digraph contains an arc of negative cost and infinite |
105 | 106 |
/// upper bound. It means that the objective function is unbounded |
106 | 107 |
/// on that arc, however, note that it could actually be bounded |
107 | 108 |
/// over the feasible flows, but this algroithm cannot handle |
108 | 109 |
/// these cases. |
109 | 110 |
UNBOUNDED |
110 | 111 |
}; |
111 | 112 |
|
112 | 113 |
/// \brief Constants for selecting the used method. |
113 | 114 |
/// |
114 | 115 |
/// Enum type containing constants for selecting the used method |
115 | 116 |
/// for the \ref run() function. |
116 | 117 |
/// |
117 | 118 |
/// \ref CycleCanceling provides three different cycle-canceling |
118 | 119 |
/// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel and Tighten" |
119 | 120 |
/// is used, which proved to be the most efficient and the most robust |
120 | 121 |
/// on various test inputs. |
121 | 122 |
/// However, the other methods can be selected using the \ref run() |
122 | 123 |
/// function with the proper parameter. |
123 | 124 |
enum Method { |
124 | 125 |
/// A simple cycle-canceling method, which uses the |
125 | 126 |
/// \ref BellmanFord "Bellman-Ford" algorithm with limited iteration |
126 | 127 |
/// number for detecting negative cycles in the residual network. |
127 | 128 |
SIMPLE_CYCLE_CANCELING, |
128 | 129 |
/// The "Minimum Mean Cycle-Canceling" algorithm, which is a |
129 | 130 |
/// well-known strongly polynomial method |
130 | 131 |
/// \ref goldberg89cyclecanceling. It improves along a |
131 | 132 |
/// \ref min_mean_cycle "minimum mean cycle" in each iteration. |
132 | 133 |
/// Its running time complexity is O(n<sup>2</sup>m<sup>3</sup>log(n)). |
133 | 134 |
MINIMUM_MEAN_CYCLE_CANCELING, |
134 | 135 |
/// The "Cancel And Tighten" algorithm, which can be viewed as an |
135 | 136 |
/// improved version of the previous method |
136 | 137 |
/// \ref goldberg89cyclecanceling. |
137 | 138 |
/// It is faster both in theory and in practice, its running time |
138 | 139 |
/// complexity is O(n<sup>2</sup>m<sup>2</sup>log(n)). |
139 | 140 |
CANCEL_AND_TIGHTEN |
140 | 141 |
}; |
141 | 142 |
|
142 | 143 |
private: |
143 | 144 |
|
144 | 145 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
145 | 146 |
|
146 | 147 |
typedef std::vector<int> IntVector; |
147 | 148 |
typedef std::vector<double> DoubleVector; |
148 | 149 |
typedef std::vector<Value> ValueVector; |
149 | 150 |
typedef std::vector<Cost> CostVector; |
150 | 151 |
typedef std::vector<char> BoolVector; |
151 | 152 |
// Note: vector<char> is used instead of vector<bool> for efficiency reasons |
152 | 153 |
|
153 | 154 |
private: |
154 | 155 |
|
155 | 156 |
template <typename KT, typename VT> |
156 | 157 |
class StaticVectorMap { |
157 | 158 |
public: |
158 | 159 |
typedef KT Key; |
159 | 160 |
typedef VT Value; |
160 | 161 |
|
161 | 162 |
StaticVectorMap(std::vector<Value>& v) : _v(v) {} |
162 | 163 |
|
163 | 164 |
const Value& operator[](const Key& key) const { |
164 | 165 |
return _v[StaticDigraph::id(key)]; |
165 | 166 |
} |
166 | 167 |
|
167 | 168 |
Value& operator[](const Key& key) { |
168 | 169 |
return _v[StaticDigraph::id(key)]; |
169 | 170 |
} |
170 | 171 |
|
171 | 172 |
void set(const Key& key, const Value& val) { |
172 | 173 |
_v[StaticDigraph::id(key)] = val; |
173 | 174 |
} |
174 | 175 |
|
175 | 176 |
private: |
176 | 177 |
std::vector<Value>& _v; |
177 | 178 |
}; |
178 | 179 |
|
179 | 180 |
typedef StaticVectorMap<StaticDigraph::Node, Cost> CostNodeMap; |
180 | 181 |
typedef StaticVectorMap<StaticDigraph::Arc, Cost> CostArcMap; |
181 | 182 |
|
182 | 183 |
private: |
183 | 184 |
|
184 | 185 |
|
185 | 186 |
// Data related to the underlying digraph |
186 | 187 |
const GR &_graph; |
187 | 188 |
int _node_num; |
188 | 189 |
int _arc_num; |
189 | 190 |
int _res_node_num; |
190 | 191 |
int _res_arc_num; |
191 | 192 |
int _root; |
192 | 193 |
|
193 | 194 |
// Parameters of the problem |
194 | 195 |
bool _have_lower; |
195 | 196 |
Value _sum_supply; |
196 | 197 |
|
197 | 198 |
// Data structures for storing the digraph |
198 | 199 |
IntNodeMap _node_id; |
199 | 200 |
IntArcMap _arc_idf; |
200 | 201 |
IntArcMap _arc_idb; |
201 | 202 |
IntVector _first_out; |
202 | 203 |
BoolVector _forward; |
203 | 204 |
IntVector _source; |
204 | 205 |
IntVector _target; |
205 | 206 |
IntVector _reverse; |
206 | 207 |
|
207 | 208 |
// Node and arc data |
208 | 209 |
ValueVector _lower; |
209 | 210 |
ValueVector _upper; |
210 | 211 |
CostVector _cost; |
211 | 212 |
ValueVector _supply; |
212 | 213 |
|
213 | 214 |
ValueVector _res_cap; |
214 | 215 |
CostVector _pi; |
215 | 216 |
|
216 | 217 |
// Data for a StaticDigraph structure |
217 | 218 |
typedef std::pair<int, int> IntPair; |
218 | 219 |
StaticDigraph _sgr; |
219 | 220 |
std::vector<IntPair> _arc_vec; |
220 | 221 |
std::vector<Cost> _cost_vec; |
221 | 222 |
IntVector _id_vec; |
222 | 223 |
CostArcMap _cost_map; |
223 | 224 |
CostNodeMap _pi_map; |
224 | 225 |
|
225 | 226 |
public: |
226 | 227 |
|
227 | 228 |
/// \brief Constant for infinite upper bounds (capacities). |
228 | 229 |
/// |
229 | 230 |
/// Constant for infinite upper bounds (capacities). |
230 | 231 |
/// It is \c std::numeric_limits<Value>::infinity() if available, |
231 | 232 |
/// \c std::numeric_limits<Value>::max() otherwise. |
232 | 233 |
const Value INF; |
233 | 234 |
|
234 | 235 |
public: |
235 | 236 |
|
236 | 237 |
/// \brief Constructor. |
237 | 238 |
/// |
238 | 239 |
/// The constructor of the class. |
239 | 240 |
/// |
240 | 241 |
/// \param graph The digraph the algorithm runs on. |
241 | 242 |
CycleCanceling(const GR& graph) : |
242 | 243 |
_graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), |
243 | 244 |
_cost_map(_cost_vec), _pi_map(_pi), |
244 | 245 |
INF(std::numeric_limits<Value>::has_infinity ? |
245 | 246 |
std::numeric_limits<Value>::infinity() : |
246 | 247 |
std::numeric_limits<Value>::max()) |
247 | 248 |
{ |
248 | 249 |
// Check the number types |
249 | 250 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
250 | 251 |
"The flow type of CycleCanceling must be signed"); |
251 | 252 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
252 | 253 |
"The cost type of CycleCanceling must be signed"); |
253 | 254 |
|
254 | 255 |
// Reset data structures |
255 | 256 |
reset(); |
256 | 257 |
} |
257 | 258 |
|
258 | 259 |
/// \name Parameters |
259 | 260 |
/// The parameters of the algorithm can be specified using these |
260 | 261 |
/// functions. |
261 | 262 |
|
262 | 263 |
/// @{ |
263 | 264 |
|
264 | 265 |
/// \brief Set the lower bounds on the arcs. |
265 | 266 |
/// |
266 | 267 |
/// This function sets the lower bounds on the arcs. |
267 | 268 |
/// If it is not used before calling \ref run(), the lower bounds |
268 | 269 |
/// will be set to zero on all arcs. |
269 | 270 |
/// |
270 | 271 |
/// \param map An arc map storing the lower bounds. |
271 | 272 |
/// Its \c Value type must be convertible to the \c Value type |
272 | 273 |
/// of the algorithm. |
273 | 274 |
/// |
274 | 275 |
/// \return <tt>(*this)</tt> |
275 | 276 |
template <typename LowerMap> |
276 | 277 |
CycleCanceling& lowerMap(const LowerMap& map) { |
277 | 278 |
_have_lower = true; |
278 | 279 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
279 | 280 |
_lower[_arc_idf[a]] = map[a]; |
280 | 281 |
_lower[_arc_idb[a]] = map[a]; |
281 | 282 |
} |
282 | 283 |
return *this; |
283 | 284 |
} |
284 | 285 |
|
285 | 286 |
/// \brief Set the upper bounds (capacities) on the arcs. |
286 | 287 |
/// |
287 | 288 |
/// This function sets the upper bounds (capacities) on the arcs. |
288 | 289 |
/// If it is not used before calling \ref run(), the upper bounds |
289 | 290 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be |
290 | 291 |
/// unbounded from above). |
291 | 292 |
/// |
292 | 293 |
/// \param map An arc map storing the upper bounds. |
293 | 294 |
/// Its \c Value type must be convertible to the \c Value type |
294 | 295 |
/// of the algorithm. |
295 | 296 |
/// |
296 | 297 |
/// \return <tt>(*this)</tt> |
297 | 298 |
template<typename UpperMap> |
298 | 299 |
CycleCanceling& upperMap(const UpperMap& map) { |
299 | 300 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
300 | 301 |
_upper[_arc_idf[a]] = map[a]; |
301 | 302 |
} |
302 | 303 |
return *this; |
303 | 304 |
} |
304 | 305 |
|
305 | 306 |
/// \brief Set the costs of the arcs. |
306 | 307 |
/// |
307 | 308 |
/// This function sets the costs of the arcs. |
308 | 309 |
/// If it is not used before calling \ref run(), the costs |
309 | 310 |
/// will be set to \c 1 on all arcs. |
310 | 311 |
/// |
311 | 312 |
/// \param map An arc map storing the costs. |
312 | 313 |
/// Its \c Value type must be convertible to the \c Cost type |
313 | 314 |
/// of the algorithm. |
314 | 315 |
/// |
315 | 316 |
/// \return <tt>(*this)</tt> |
316 | 317 |
template<typename CostMap> |
317 | 318 |
CycleCanceling& costMap(const CostMap& map) { |
318 | 319 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
319 | 320 |
_cost[_arc_idf[a]] = map[a]; |
320 | 321 |
_cost[_arc_idb[a]] = -map[a]; |
321 | 322 |
} |
322 | 323 |
return *this; |
323 | 324 |
} |
324 | 325 |
|
325 | 326 |
/// \brief Set the supply values of the nodes. |
326 | 327 |
/// |
327 | 328 |
/// This function sets the supply values of the nodes. |
328 | 329 |
/// If neither this function nor \ref stSupply() is used before |
329 | 330 |
/// calling \ref run(), the supply of each node will be set to zero. |
330 | 331 |
/// |
331 | 332 |
/// \param map A node map storing the supply values. |
332 | 333 |
/// Its \c Value type must be convertible to the \c Value type |
333 | 334 |
/// of the algorithm. |
334 | 335 |
/// |
335 | 336 |
/// \return <tt>(*this)</tt> |
336 | 337 |
template<typename SupplyMap> |
337 | 338 |
CycleCanceling& supplyMap(const SupplyMap& map) { |
338 | 339 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
339 | 340 |
_supply[_node_id[n]] = map[n]; |
340 | 341 |
} |
341 | 342 |
return *this; |
342 | 343 |
} |
343 | 344 |
|
344 | 345 |
/// \brief Set single source and target nodes and a supply value. |
345 | 346 |
/// |
346 | 347 |
/// This function sets a single source node and a single target node |
347 | 348 |
/// and the required flow value. |
348 | 349 |
/// If neither this function nor \ref supplyMap() is used before |
349 | 350 |
/// calling \ref run(), the supply of each node will be set to zero. |
350 | 351 |
/// |
351 | 352 |
/// Using this function has the same effect as using \ref supplyMap() |
352 | 353 |
/// with such a map in which \c k is assigned to \c s, \c -k is |
353 | 354 |
/// assigned to \c t and all other nodes have zero supply value. |
354 | 355 |
/// |
355 | 356 |
/// \param s The source node. |
356 | 357 |
/// \param t The target node. |
357 | 358 |
/// \param k The required amount of flow from node \c s to node \c t |
358 | 359 |
/// (i.e. the supply of \c s and the demand of \c t). |
359 | 360 |
/// |
360 | 361 |
/// \return <tt>(*this)</tt> |
361 | 362 |
CycleCanceling& stSupply(const Node& s, const Node& t, Value k) { |
362 | 363 |
for (int i = 0; i != _res_node_num; ++i) { |
363 | 364 |
_supply[i] = 0; |
364 | 365 |
} |
365 | 366 |
_supply[_node_id[s]] = k; |
366 | 367 |
_supply[_node_id[t]] = -k; |
367 | 368 |
return *this; |
368 | 369 |
} |
369 | 370 |
|
370 | 371 |
/// @} |
371 | 372 |
|
372 | 373 |
/// \name Execution control |
373 | 374 |
/// The algorithm can be executed using \ref run(). |
374 | 375 |
|
375 | 376 |
/// @{ |
376 | 377 |
|
377 | 378 |
/// \brief Run the algorithm. |
378 | 379 |
/// |
379 | 380 |
/// This function runs the algorithm. |
380 | 381 |
/// The paramters can be specified using functions \ref lowerMap(), |
381 | 382 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). |
382 | 383 |
/// For example, |
383 | 384 |
/// \code |
384 | 385 |
/// CycleCanceling<ListDigraph> cc(graph); |
385 | 386 |
/// cc.lowerMap(lower).upperMap(upper).costMap(cost) |
386 | 387 |
/// .supplyMap(sup).run(); |
387 | 388 |
/// \endcode |
388 | 389 |
/// |
389 | 390 |
/// This function can be called more than once. All the given parameters |
390 | 391 |
/// are kept for the next call, unless \ref resetParams() or \ref reset() |
391 | 392 |
/// is used, thus only the modified parameters have to be set again. |
392 | 393 |
/// If the underlying digraph was also modified after the construction |
393 | 394 |
/// of the class (or the last \ref reset() call), then the \ref reset() |
394 | 395 |
/// function must be called. |
395 | 396 |
/// |
396 | 397 |
/// \param method The cycle-canceling method that will be used. |
397 | 398 |
/// For more information, see \ref Method. |
398 | 399 |
/// |
399 | 400 |
/// \return \c INFEASIBLE if no feasible flow exists, |
400 | 401 |
/// \n \c OPTIMAL if the problem has optimal solution |
401 | 402 |
/// (i.e. it is feasible and bounded), and the algorithm has found |
402 | 403 |
/// optimal flow and node potentials (primal and dual solutions), |
403 | 404 |
/// \n \c UNBOUNDED if the digraph contains an arc of negative cost |
404 | 405 |
/// and infinite upper bound. It means that the objective function |
405 | 406 |
/// is unbounded on that arc, however, note that it could actually be |
406 | 407 |
/// bounded over the feasible flows, but this algroithm cannot handle |
407 | 408 |
/// these cases. |
408 | 409 |
/// |
409 | 410 |
/// \see ProblemType, Method |
410 | 411 |
/// \see resetParams(), reset() |
411 | 412 |
ProblemType run(Method method = CANCEL_AND_TIGHTEN) { |
412 | 413 |
ProblemType pt = init(); |
413 | 414 |
if (pt != OPTIMAL) return pt; |
414 | 415 |
start(method); |
415 | 416 |
return OPTIMAL; |
416 | 417 |
} |
417 | 418 |
|
418 | 419 |
/// \brief Reset all the parameters that have been given before. |
419 | 420 |
/// |
420 | 421 |
/// This function resets all the paramaters that have been given |
421 | 422 |
/// before using functions \ref lowerMap(), \ref upperMap(), |
422 | 423 |
/// \ref costMap(), \ref supplyMap(), \ref stSupply(). |
423 | 424 |
/// |
424 | 425 |
/// It is useful for multiple \ref run() calls. Basically, all the given |
425 | 426 |
/// parameters are kept for the next \ref run() call, unless |
426 | 427 |
/// \ref resetParams() or \ref reset() is used. |
427 | 428 |
/// If the underlying digraph was also modified after the construction |
428 | 429 |
/// of the class or the last \ref reset() call, then the \ref reset() |
429 | 430 |
/// function must be used, otherwise \ref resetParams() is sufficient. |
430 | 431 |
/// |
431 | 432 |
/// For example, |
432 | 433 |
/// \code |
433 | 434 |
/// CycleCanceling<ListDigraph> cs(graph); |
434 | 435 |
/// |
435 | 436 |
/// // First run |
436 | 437 |
/// cc.lowerMap(lower).upperMap(upper).costMap(cost) |
437 | 438 |
/// .supplyMap(sup).run(); |
438 | 439 |
/// |
439 | 440 |
/// // Run again with modified cost map (resetParams() is not called, |
440 | 441 |
/// // so only the cost map have to be set again) |
441 | 442 |
/// cost[e] += 100; |
442 | 443 |
/// cc.costMap(cost).run(); |
443 | 444 |
/// |
444 | 445 |
/// // Run again from scratch using resetParams() |
445 | 446 |
/// // (the lower bounds will be set to zero on all arcs) |
446 | 447 |
/// cc.resetParams(); |
447 | 448 |
/// cc.upperMap(capacity).costMap(cost) |
448 | 449 |
/// .supplyMap(sup).run(); |
449 | 450 |
/// \endcode |
450 | 451 |
/// |
451 | 452 |
/// \return <tt>(*this)</tt> |
452 | 453 |
/// |
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_KRUSKAL_H |
20 | 20 |
#define LEMON_KRUSKAL_H |
21 | 21 |
|
22 | 22 |
#include <algorithm> |
23 | 23 |
#include <vector> |
24 | 24 |
#include <lemon/unionfind.h> |
25 | 25 |
#include <lemon/maps.h> |
26 | 26 |
|
27 | 27 |
#include <lemon/core.h> |
28 | 28 |
#include <lemon/bits/traits.h> |
29 | 29 |
|
30 | 30 |
///\ingroup spantree |
31 | 31 |
///\file |
32 | 32 |
///\brief Kruskal's algorithm to compute a minimum cost spanning tree |
33 |
/// |
|
34 |
///Kruskal's algorithm to compute a minimum cost spanning tree. |
|
35 |
/// |
|
36 | 33 |
|
37 | 34 |
namespace lemon { |
38 | 35 |
|
39 | 36 |
namespace _kruskal_bits { |
40 | 37 |
|
41 | 38 |
// Kruskal for directed graphs. |
42 | 39 |
|
43 | 40 |
template <typename Digraph, typename In, typename Out> |
44 | 41 |
typename disable_if<lemon::UndirectedTagIndicator<Digraph>, |
45 | 42 |
typename In::value_type::second_type >::type |
46 | 43 |
kruskal(const Digraph& digraph, const In& in, Out& out,dummy<0> = 0) { |
47 | 44 |
typedef typename In::value_type::second_type Value; |
48 | 45 |
typedef typename Digraph::template NodeMap<int> IndexMap; |
49 | 46 |
typedef typename Digraph::Node Node; |
50 | 47 |
|
51 | 48 |
IndexMap index(digraph); |
52 | 49 |
UnionFind<IndexMap> uf(index); |
53 | 50 |
for (typename Digraph::NodeIt it(digraph); it != INVALID; ++it) { |
54 | 51 |
uf.insert(it); |
55 | 52 |
} |
56 | 53 |
|
57 | 54 |
Value tree_value = 0; |
58 | 55 |
for (typename In::const_iterator it = in.begin(); it != in.end(); ++it) { |
59 | 56 |
if (uf.join(digraph.target(it->first),digraph.source(it->first))) { |
60 | 57 |
out.set(it->first, true); |
61 | 58 |
tree_value += it->second; |
62 | 59 |
} |
63 | 60 |
else { |
64 | 61 |
out.set(it->first, false); |
65 | 62 |
} |
66 | 63 |
} |
67 | 64 |
return tree_value; |
68 | 65 |
} |
69 | 66 |
|
70 | 67 |
// Kruskal for undirected graphs. |
71 | 68 |
|
72 | 69 |
template <typename Graph, typename In, typename Out> |
73 | 70 |
typename enable_if<lemon::UndirectedTagIndicator<Graph>, |
74 | 71 |
typename In::value_type::second_type >::type |
75 | 72 |
kruskal(const Graph& graph, const In& in, Out& out,dummy<1> = 1) { |
76 | 73 |
typedef typename In::value_type::second_type Value; |
77 | 74 |
typedef typename Graph::template NodeMap<int> IndexMap; |
78 | 75 |
typedef typename Graph::Node Node; |
79 | 76 |
|
80 | 77 |
IndexMap index(graph); |
81 | 78 |
UnionFind<IndexMap> uf(index); |
82 | 79 |
for (typename Graph::NodeIt it(graph); it != INVALID; ++it) { |
83 | 80 |
uf.insert(it); |
84 | 81 |
} |
85 | 82 |
|
86 | 83 |
Value tree_value = 0; |
87 | 84 |
for (typename In::const_iterator it = in.begin(); it != in.end(); ++it) { |
88 | 85 |
if (uf.join(graph.u(it->first),graph.v(it->first))) { |
89 | 86 |
out.set(it->first, true); |
90 | 87 |
tree_value += it->second; |
91 | 88 |
} |
92 | 89 |
else { |
93 | 90 |
out.set(it->first, false); |
94 | 91 |
} |
95 | 92 |
} |
96 | 93 |
return tree_value; |
97 | 94 |
} |
98 | 95 |
|
99 | 96 |
|
100 | 97 |
template <typename Sequence> |
101 | 98 |
struct PairComp { |
102 | 99 |
typedef typename Sequence::value_type Value; |
103 | 100 |
bool operator()(const Value& left, const Value& right) { |
104 | 101 |
return left.second < right.second; |
105 | 102 |
} |
106 | 103 |
}; |
107 | 104 |
|
108 | 105 |
template <typename In, typename Enable = void> |
109 | 106 |
struct SequenceInputIndicator { |
110 | 107 |
static const bool value = false; |
111 | 108 |
}; |
112 | 109 |
|
113 | 110 |
template <typename In> |
114 | 111 |
struct SequenceInputIndicator<In, |
115 | 112 |
typename exists<typename In::value_type::first_type>::type> { |
116 | 113 |
static const bool value = true; |
117 | 114 |
}; |
118 | 115 |
|
119 | 116 |
template <typename In, typename Enable = void> |
120 | 117 |
struct MapInputIndicator { |
121 | 118 |
static const bool value = false; |
122 | 119 |
}; |
123 | 120 |
|
124 | 121 |
template <typename In> |
125 | 122 |
struct MapInputIndicator<In, |
126 | 123 |
typename exists<typename In::Value>::type> { |
127 | 124 |
static const bool value = true; |
128 | 125 |
}; |
129 | 126 |
|
130 | 127 |
template <typename In, typename Enable = void> |
131 | 128 |
struct SequenceOutputIndicator { |
132 | 129 |
static const bool value = false; |
133 | 130 |
}; |
134 | 131 |
|
135 | 132 |
template <typename Out> |
136 | 133 |
struct SequenceOutputIndicator<Out, |
137 | 134 |
typename exists<typename Out::value_type>::type> { |
138 | 135 |
static const bool value = true; |
139 | 136 |
}; |
140 | 137 |
|
141 | 138 |
template <typename Out, typename Enable = void> |
142 | 139 |
struct MapOutputIndicator { |
143 | 140 |
static const bool value = false; |
144 | 141 |
}; |
145 | 142 |
|
146 | 143 |
template <typename Out> |
147 | 144 |
struct MapOutputIndicator<Out, |
148 | 145 |
typename exists<typename Out::Value>::type> { |
149 | 146 |
static const bool value = true; |
150 | 147 |
}; |
151 | 148 |
|
152 | 149 |
template <typename In, typename InEnable = void> |
153 | 150 |
struct KruskalValueSelector {}; |
154 | 151 |
|
155 | 152 |
template <typename In> |
156 | 153 |
struct KruskalValueSelector<In, |
157 | 154 |
typename enable_if<SequenceInputIndicator<In>, void>::type> |
158 | 155 |
{ |
159 | 156 |
typedef typename In::value_type::second_type Value; |
160 | 157 |
}; |
161 | 158 |
|
162 | 159 |
template <typename In> |
163 | 160 |
struct KruskalValueSelector<In, |
164 | 161 |
typename enable_if<MapInputIndicator<In>, void>::type> |
165 | 162 |
{ |
166 | 163 |
typedef typename In::Value Value; |
167 | 164 |
}; |
168 | 165 |
|
169 | 166 |
template <typename Graph, typename In, typename Out, |
170 | 167 |
typename InEnable = void> |
171 | 168 |
struct KruskalInputSelector {}; |
172 | 169 |
|
173 | 170 |
template <typename Graph, typename In, typename Out, |
174 | 171 |
typename InEnable = void> |
175 | 172 |
struct KruskalOutputSelector {}; |
176 | 173 |
|
177 | 174 |
template <typename Graph, typename In, typename Out> |
178 | 175 |
struct KruskalInputSelector<Graph, In, Out, |
179 | 176 |
typename enable_if<SequenceInputIndicator<In>, void>::type > |
180 | 177 |
{ |
181 | 178 |
typedef typename In::value_type::second_type Value; |
182 | 179 |
|
183 | 180 |
static Value kruskal(const Graph& graph, const In& in, Out& out) { |
184 | 181 |
return KruskalOutputSelector<Graph, In, Out>:: |
185 | 182 |
kruskal(graph, in, out); |
186 | 183 |
} |
187 | 184 |
|
188 | 185 |
}; |
189 | 186 |
|
190 | 187 |
template <typename Graph, typename In, typename Out> |
191 | 188 |
struct KruskalInputSelector<Graph, In, Out, |
192 | 189 |
typename enable_if<MapInputIndicator<In>, void>::type > |
193 | 190 |
{ |
194 | 191 |
typedef typename In::Value Value; |
195 | 192 |
static Value kruskal(const Graph& graph, const In& in, Out& out) { |
196 | 193 |
typedef typename In::Key MapArc; |
197 | 194 |
typedef typename In::Value Value; |
198 | 195 |
typedef typename ItemSetTraits<Graph, MapArc>::ItemIt MapArcIt; |
199 | 196 |
typedef std::vector<std::pair<MapArc, Value> > Sequence; |
200 | 197 |
Sequence seq; |
201 | 198 |
|
202 | 199 |
for (MapArcIt it(graph); it != INVALID; ++it) { |
203 | 200 |
seq.push_back(std::make_pair(it, in[it])); |
204 | 201 |
} |
205 | 202 |
|
206 | 203 |
std::sort(seq.begin(), seq.end(), PairComp<Sequence>()); |
207 | 204 |
return KruskalOutputSelector<Graph, Sequence, Out>:: |
208 | 205 |
kruskal(graph, seq, out); |
209 | 206 |
} |
210 | 207 |
}; |
211 | 208 |
|
212 | 209 |
template <typename T> |
213 | 210 |
struct RemoveConst { |
214 | 211 |
typedef T type; |
215 | 212 |
}; |
216 | 213 |
|
217 | 214 |
template <typename T> |
218 | 215 |
struct RemoveConst<const T> { |
219 | 216 |
typedef T type; |
220 | 217 |
}; |
221 | 218 |
|
222 | 219 |
template <typename Graph, typename In, typename Out> |
223 | 220 |
struct KruskalOutputSelector<Graph, In, Out, |
224 | 221 |
typename enable_if<SequenceOutputIndicator<Out>, void>::type > |
225 | 222 |
{ |
226 | 223 |
typedef typename In::value_type::second_type Value; |
227 | 224 |
|
228 | 225 |
static Value kruskal(const Graph& graph, const In& in, Out& out) { |
229 | 226 |
typedef LoggerBoolMap<typename RemoveConst<Out>::type> Map; |
230 | 227 |
Map map(out); |
231 | 228 |
return _kruskal_bits::kruskal(graph, in, map); |
232 | 229 |
} |
233 | 230 |
|
234 | 231 |
}; |
235 | 232 |
|
236 | 233 |
template <typename Graph, typename In, typename Out> |
237 | 234 |
struct KruskalOutputSelector<Graph, In, Out, |
238 | 235 |
typename enable_if<MapOutputIndicator<Out>, void>::type > |
239 | 236 |
{ |
240 | 237 |
typedef typename In::value_type::second_type Value; |
241 | 238 |
|
242 | 239 |
static Value kruskal(const Graph& graph, const In& in, Out& out) { |
243 | 240 |
return _kruskal_bits::kruskal(graph, in, out); |
244 | 241 |
} |
245 | 242 |
}; |
246 | 243 |
|
247 | 244 |
} |
248 | 245 |
|
249 | 246 |
/// \ingroup spantree |
250 | 247 |
/// |
251 | 248 |
/// \brief Kruskal's algorithm for finding a minimum cost spanning tree of |
252 | 249 |
/// a graph. |
253 | 250 |
/// |
254 | 251 |
/// This function runs Kruskal's algorithm to find a minimum cost |
255 | 252 |
/// spanning tree of a graph. |
256 | 253 |
/// Due to some C++ hacking, it accepts various input and output types. |
257 | 254 |
/// |
258 | 255 |
/// \param g The graph the algorithm runs on. |
259 | 256 |
/// It can be either \ref concepts::Digraph "directed" or |
260 | 257 |
/// \ref concepts::Graph "undirected". |
261 | 258 |
/// If the graph is directed, the algorithm consider it to be |
262 | 259 |
/// undirected by disregarding the direction of the arcs. |
263 | 260 |
/// |
264 | 261 |
/// \param in This object is used to describe the arc/edge costs. |
265 | 262 |
/// It can be one of the following choices. |
266 | 263 |
/// - An STL compatible 'Forward Container' with |
267 | 264 |
/// <tt>std::pair<GR::Arc,C></tt> or |
268 | 265 |
/// <tt>std::pair<GR::Edge,C></tt> as its <tt>value_type</tt>, where |
269 | 266 |
/// \c C is the type of the costs. The pairs indicates the arcs/edges |
270 | 267 |
/// along with the assigned cost. <em>They must be in a |
271 | 268 |
/// cost-ascending order.</em> |
272 | 269 |
/// - Any readable arc/edge map. The values of the map indicate the |
273 | 270 |
/// arc/edge costs. |
274 | 271 |
/// |
275 | 272 |
/// \retval out Here we also have a choice. |
276 | 273 |
/// - It can be a writable arc/edge map with \c bool value type. After |
277 | 274 |
/// running the algorithm it will contain the found minimum cost spanning |
278 | 275 |
/// tree: the value of an arc/edge will be set to \c true if it belongs |
279 | 276 |
/// to the tree, otherwise it will be set to \c false. The value of |
280 | 277 |
/// each arc/edge will be set exactly once. |
281 | 278 |
/// - It can also be an iteraror of an STL Container with |
282 | 279 |
/// <tt>GR::Arc</tt> or <tt>GR::Edge</tt> as its |
283 | 280 |
/// <tt>value_type</tt>. The algorithm copies the elements of the |
284 | 281 |
/// found tree into this sequence. For example, if we know that the |
285 | 282 |
/// spanning tree of the graph \c g has say 53 arcs, then we can |
286 | 283 |
/// put its arcs into an STL vector \c tree with a code like this. |
287 | 284 |
///\code |
288 | 285 |
/// std::vector<Arc> tree(53); |
289 | 286 |
/// kruskal(g,cost,tree.begin()); |
290 | 287 |
///\endcode |
291 | 288 |
/// Or if we don't know in advance the size of the tree, we can |
292 | 289 |
/// write this. |
293 | 290 |
///\code |
294 | 291 |
/// std::vector<Arc> tree; |
295 | 292 |
/// kruskal(g,cost,std::back_inserter(tree)); |
296 | 293 |
///\endcode |
297 | 294 |
/// |
298 | 295 |
/// \return The total cost of the found spanning tree. |
299 | 296 |
/// |
300 | 297 |
/// \note If the input graph is not (weakly) connected, a spanning |
301 | 298 |
/// forest is calculated instead of a spanning tree. |
302 | 299 |
|
303 | 300 |
#ifdef DOXYGEN |
304 | 301 |
template <typename Graph, typename In, typename Out> |
305 | 302 |
Value kruskal(const Graph& g, const In& in, Out& out) |
306 | 303 |
#else |
307 | 304 |
template <class Graph, class In, class Out> |
308 | 305 |
inline typename _kruskal_bits::KruskalValueSelector<In>::Value |
309 | 306 |
kruskal(const Graph& graph, const In& in, Out& out) |
310 | 307 |
#endif |
311 | 308 |
{ |
312 | 309 |
return _kruskal_bits::KruskalInputSelector<Graph, In, Out>:: |
313 | 310 |
kruskal(graph, in, out); |
314 | 311 |
} |
315 | 312 |
|
316 | 313 |
|
317 | 314 |
template <class Graph, class In, class Out> |
318 | 315 |
inline typename _kruskal_bits::KruskalValueSelector<In>::Value |
319 | 316 |
kruskal(const Graph& graph, const In& in, const Out& out) |
320 | 317 |
{ |
321 | 318 |
return _kruskal_bits::KruskalInputSelector<Graph, In, const Out>:: |
322 | 319 |
kruskal(graph, in, out); |
323 | 320 |
} |
324 | 321 |
|
325 | 322 |
} //namespace lemon |
326 | 323 |
|
327 | 324 |
#endif //LEMON_KRUSKAL_H |
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-2010 |
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_NETWORK_SIMPLEX_H |
20 | 20 |
#define LEMON_NETWORK_SIMPLEX_H |
21 | 21 |
|
22 | 22 |
/// \ingroup min_cost_flow_algs |
23 | 23 |
/// |
24 | 24 |
/// \file |
25 | 25 |
/// \brief Network Simplex algorithm for finding a minimum cost flow. |
26 | 26 |
|
27 | 27 |
#include <vector> |
28 | 28 |
#include <limits> |
29 | 29 |
#include <algorithm> |
30 | 30 |
|
31 | 31 |
#include <lemon/core.h> |
32 | 32 |
#include <lemon/math.h> |
33 | 33 |
|
34 | 34 |
namespace lemon { |
35 | 35 |
|
36 | 36 |
/// \addtogroup min_cost_flow_algs |
37 | 37 |
/// @{ |
38 | 38 |
|
39 | 39 |
/// \brief Implementation of the primal Network Simplex algorithm |
40 | 40 |
/// for finding a \ref min_cost_flow "minimum cost flow". |
41 | 41 |
/// |
42 | 42 |
/// \ref NetworkSimplex implements the primal Network Simplex algorithm |
43 | 43 |
/// for finding a \ref min_cost_flow "minimum cost flow" |
44 | 44 |
/// \ref amo93networkflows, \ref dantzig63linearprog, |
45 | 45 |
/// \ref kellyoneill91netsimplex. |
46 | 46 |
/// This algorithm is a highly efficient specialized version of the |
47 | 47 |
/// linear programming simplex method directly for the minimum cost |
48 | 48 |
/// flow problem. |
49 | 49 |
/// |
50 | 50 |
/// In general, %NetworkSimplex is the fastest implementation available |
51 | 51 |
/// in LEMON for this problem. |
52 | 52 |
/// Moreover, it supports both directions of the supply/demand inequality |
53 | 53 |
/// constraints. For more information, see \ref SupplyType. |
54 | 54 |
/// |
55 | 55 |
/// Most of the parameters of the problem (except for the digraph) |
56 | 56 |
/// can be given using separate functions, and the algorithm can be |
57 | 57 |
/// executed using the \ref run() function. If some parameters are not |
58 | 58 |
/// specified, then default values will be used. |
59 | 59 |
/// |
60 | 60 |
/// \tparam GR The digraph type the algorithm runs on. |
61 | 61 |
/// \tparam V The number type used for flow amounts, capacity bounds |
62 | 62 |
/// and supply values in the algorithm. By default, it is \c int. |
63 | 63 |
/// \tparam C The number type used for costs and potentials in the |
64 | 64 |
/// algorithm. By default, it is the same as \c V. |
65 | 65 |
/// |
66 |
/// \warning Both |
|
66 |
/// \warning Both \c V and \c C must be signed number types. |
|
67 |
/// \warning All input data (capacities, supply values, and costs) must |
|
67 | 68 |
/// be integer. |
68 | 69 |
/// |
69 | 70 |
/// \note %NetworkSimplex provides five different pivot rule |
70 | 71 |
/// implementations, from which the most efficient one is used |
71 | 72 |
/// by default. For more information, see \ref PivotRule. |
72 | 73 |
template <typename GR, typename V = int, typename C = V> |
73 | 74 |
class NetworkSimplex |
74 | 75 |
{ |
75 | 76 |
public: |
76 | 77 |
|
77 | 78 |
/// The type of the flow amounts, capacity bounds and supply values |
78 | 79 |
typedef V Value; |
79 | 80 |
/// The type of the arc costs |
80 | 81 |
typedef C Cost; |
81 | 82 |
|
82 | 83 |
public: |
83 | 84 |
|
84 | 85 |
/// \brief Problem type constants for the \c run() function. |
85 | 86 |
/// |
86 | 87 |
/// Enum type containing the problem type constants that can be |
87 | 88 |
/// returned by the \ref run() function of the algorithm. |
88 | 89 |
enum ProblemType { |
89 | 90 |
/// The problem has no feasible solution (flow). |
90 | 91 |
INFEASIBLE, |
91 | 92 |
/// The problem has optimal solution (i.e. it is feasible and |
92 | 93 |
/// bounded), and the algorithm has found optimal flow and node |
93 | 94 |
/// potentials (primal and dual solutions). |
94 | 95 |
OPTIMAL, |
95 | 96 |
/// The objective function of the problem is unbounded, i.e. |
96 | 97 |
/// there is a directed cycle having negative total cost and |
97 | 98 |
/// infinite upper bound. |
98 | 99 |
UNBOUNDED |
99 | 100 |
}; |
100 | 101 |
|
101 | 102 |
/// \brief Constants for selecting the type of the supply constraints. |
102 | 103 |
/// |
103 | 104 |
/// Enum type containing constants for selecting the supply type, |
104 | 105 |
/// i.e. the direction of the inequalities in the supply/demand |
105 | 106 |
/// constraints of the \ref min_cost_flow "minimum cost flow problem". |
106 | 107 |
/// |
107 | 108 |
/// The default supply type is \c GEQ, the \c LEQ type can be |
108 | 109 |
/// selected using \ref supplyType(). |
109 | 110 |
/// The equality form is a special case of both supply types. |
110 | 111 |
enum SupplyType { |
111 | 112 |
/// This option means that there are <em>"greater or equal"</em> |
112 | 113 |
/// supply/demand constraints in the definition of the problem. |
113 | 114 |
GEQ, |
114 | 115 |
/// This option means that there are <em>"less or equal"</em> |
115 | 116 |
/// supply/demand constraints in the definition of the problem. |
116 | 117 |
LEQ |
117 | 118 |
}; |
118 | 119 |
|
119 | 120 |
/// \brief Constants for selecting the pivot rule. |
120 | 121 |
/// |
121 | 122 |
/// Enum type containing constants for selecting the pivot rule for |
122 | 123 |
/// the \ref run() function. |
123 | 124 |
/// |
124 | 125 |
/// \ref NetworkSimplex provides five different pivot rule |
125 | 126 |
/// implementations that significantly affect the running time |
126 | 127 |
/// of the algorithm. |
127 | 128 |
/// By default, \ref BLOCK_SEARCH "Block Search" is used, which |
128 | 129 |
/// proved to be the most efficient and the most robust on various |
129 | 130 |
/// test inputs. |
130 | 131 |
/// However, another pivot rule can be selected using the \ref run() |
131 | 132 |
/// function with the proper parameter. |
132 | 133 |
enum PivotRule { |
133 | 134 |
|
134 | 135 |
/// The \e First \e Eligible pivot rule. |
135 | 136 |
/// The next eligible arc is selected in a wraparound fashion |
136 | 137 |
/// in every iteration. |
137 | 138 |
FIRST_ELIGIBLE, |
138 | 139 |
|
139 | 140 |
/// The \e Best \e Eligible pivot rule. |
140 | 141 |
/// The best eligible arc is selected in every iteration. |
141 | 142 |
BEST_ELIGIBLE, |
142 | 143 |
|
143 | 144 |
/// The \e Block \e Search pivot rule. |
144 | 145 |
/// A specified number of arcs are examined in every iteration |
145 | 146 |
/// in a wraparound fashion and the best eligible arc is selected |
146 | 147 |
/// from this block. |
147 | 148 |
BLOCK_SEARCH, |
148 | 149 |
|
149 | 150 |
/// The \e Candidate \e List pivot rule. |
150 | 151 |
/// In a major iteration a candidate list is built from eligible arcs |
151 | 152 |
/// in a wraparound fashion and in the following minor iterations |
152 | 153 |
/// the best eligible arc is selected from this list. |
153 | 154 |
CANDIDATE_LIST, |
154 | 155 |
|
155 | 156 |
/// The \e Altering \e Candidate \e List pivot rule. |
156 | 157 |
/// It is a modified version of the Candidate List method. |
157 | 158 |
/// It keeps only the several best eligible arcs from the former |
158 | 159 |
/// candidate list and extends this list in every iteration. |
159 | 160 |
ALTERING_LIST |
160 | 161 |
}; |
161 | 162 |
|
162 | 163 |
private: |
163 | 164 |
|
164 | 165 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
165 | 166 |
|
166 | 167 |
typedef std::vector<int> IntVector; |
167 | 168 |
typedef std::vector<Value> ValueVector; |
168 | 169 |
typedef std::vector<Cost> CostVector; |
169 | 170 |
typedef std::vector<signed char> CharVector; |
170 | 171 |
// Note: vector<signed char> is used instead of vector<ArcState> and |
171 | 172 |
// vector<ArcDirection> for efficiency reasons |
172 | 173 |
|
173 | 174 |
// State constants for arcs |
174 | 175 |
enum ArcState { |
175 | 176 |
STATE_UPPER = -1, |
176 | 177 |
STATE_TREE = 0, |
177 | 178 |
STATE_LOWER = 1 |
178 | 179 |
}; |
179 | 180 |
|
180 | 181 |
// Direction constants for tree arcs |
181 | 182 |
enum ArcDirection { |
182 | 183 |
DIR_DOWN = -1, |
183 | 184 |
DIR_UP = 1 |
184 | 185 |
}; |
185 | 186 |
|
186 | 187 |
private: |
187 | 188 |
|
188 | 189 |
// Data related to the underlying digraph |
189 | 190 |
const GR &_graph; |
190 | 191 |
int _node_num; |
191 | 192 |
int _arc_num; |
192 | 193 |
int _all_arc_num; |
193 | 194 |
int _search_arc_num; |
194 | 195 |
|
195 | 196 |
// Parameters of the problem |
196 | 197 |
bool _have_lower; |
197 | 198 |
SupplyType _stype; |
198 | 199 |
Value _sum_supply; |
199 | 200 |
|
200 | 201 |
// Data structures for storing the digraph |
201 | 202 |
IntNodeMap _node_id; |
202 | 203 |
IntArcMap _arc_id; |
203 | 204 |
IntVector _source; |
204 | 205 |
IntVector _target; |
205 | 206 |
bool _arc_mixing; |
206 | 207 |
|
207 | 208 |
// Node and arc data |
208 | 209 |
ValueVector _lower; |
209 | 210 |
ValueVector _upper; |
210 | 211 |
ValueVector _cap; |
211 | 212 |
CostVector _cost; |
212 | 213 |
ValueVector _supply; |
213 | 214 |
ValueVector _flow; |
214 | 215 |
CostVector _pi; |
215 | 216 |
|
216 | 217 |
// Data for storing the spanning tree structure |
217 | 218 |
IntVector _parent; |
218 | 219 |
IntVector _pred; |
219 | 220 |
IntVector _thread; |
220 | 221 |
IntVector _rev_thread; |
221 | 222 |
IntVector _succ_num; |
222 | 223 |
IntVector _last_succ; |
223 | 224 |
CharVector _pred_dir; |
224 | 225 |
CharVector _state; |
225 | 226 |
IntVector _dirty_revs; |
226 | 227 |
int _root; |
227 | 228 |
|
228 | 229 |
// Temporary data used in the current pivot iteration |
229 | 230 |
int in_arc, join, u_in, v_in, u_out, v_out; |
230 | 231 |
Value delta; |
231 | 232 |
|
232 | 233 |
const Value MAX; |
233 | 234 |
|
234 | 235 |
public: |
235 | 236 |
|
236 | 237 |
/// \brief Constant for infinite upper bounds (capacities). |
237 | 238 |
/// |
238 | 239 |
/// Constant for infinite upper bounds (capacities). |
239 | 240 |
/// It is \c std::numeric_limits<Value>::infinity() if available, |
240 | 241 |
/// \c std::numeric_limits<Value>::max() otherwise. |
241 | 242 |
const Value INF; |
242 | 243 |
|
243 | 244 |
private: |
244 | 245 |
|
245 | 246 |
// Implementation of the First Eligible pivot rule |
246 | 247 |
class FirstEligiblePivotRule |
247 | 248 |
{ |
248 | 249 |
private: |
249 | 250 |
|
250 | 251 |
// References to the NetworkSimplex class |
251 | 252 |
const IntVector &_source; |
252 | 253 |
const IntVector &_target; |
253 | 254 |
const CostVector &_cost; |
254 | 255 |
const CharVector &_state; |
255 | 256 |
const CostVector &_pi; |
256 | 257 |
int &_in_arc; |
257 | 258 |
int _search_arc_num; |
258 | 259 |
|
259 | 260 |
// Pivot rule data |
260 | 261 |
int _next_arc; |
261 | 262 |
|
262 | 263 |
public: |
263 | 264 |
|
264 | 265 |
// Constructor |
265 | 266 |
FirstEligiblePivotRule(NetworkSimplex &ns) : |
266 | 267 |
_source(ns._source), _target(ns._target), |
267 | 268 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
268 | 269 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
269 | 270 |
_next_arc(0) |
270 | 271 |
{} |
271 | 272 |
|
272 | 273 |
// Find next entering arc |
273 | 274 |
bool findEnteringArc() { |
274 | 275 |
Cost c; |
275 | 276 |
for (int e = _next_arc; e != _search_arc_num; ++e) { |
276 | 277 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
277 | 278 |
if (c < 0) { |
278 | 279 |
_in_arc = e; |
279 | 280 |
_next_arc = e + 1; |
280 | 281 |
return true; |
281 | 282 |
} |
282 | 283 |
} |
283 | 284 |
for (int e = 0; e != _next_arc; ++e) { |
284 | 285 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
285 | 286 |
if (c < 0) { |
286 | 287 |
_in_arc = e; |
287 | 288 |
_next_arc = e + 1; |
288 | 289 |
return true; |
289 | 290 |
} |
290 | 291 |
} |
291 | 292 |
return false; |
292 | 293 |
} |
293 | 294 |
|
294 | 295 |
}; //class FirstEligiblePivotRule |
295 | 296 |
|
296 | 297 |
|
297 | 298 |
// Implementation of the Best Eligible pivot rule |
298 | 299 |
class BestEligiblePivotRule |
299 | 300 |
{ |
300 | 301 |
private: |
301 | 302 |
|
302 | 303 |
// References to the NetworkSimplex class |
303 | 304 |
const IntVector &_source; |
304 | 305 |
const IntVector &_target; |
305 | 306 |
const CostVector &_cost; |
306 | 307 |
const CharVector &_state; |
307 | 308 |
const CostVector &_pi; |
308 | 309 |
int &_in_arc; |
309 | 310 |
int _search_arc_num; |
310 | 311 |
|
311 | 312 |
public: |
312 | 313 |
|
313 | 314 |
// Constructor |
314 | 315 |
BestEligiblePivotRule(NetworkSimplex &ns) : |
315 | 316 |
_source(ns._source), _target(ns._target), |
316 | 317 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
317 | 318 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num) |
318 | 319 |
{} |
319 | 320 |
|
320 | 321 |
// Find next entering arc |
321 | 322 |
bool findEnteringArc() { |
322 | 323 |
Cost c, min = 0; |
323 | 324 |
for (int e = 0; e != _search_arc_num; ++e) { |
324 | 325 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
325 | 326 |
if (c < min) { |
326 | 327 |
min = c; |
327 | 328 |
_in_arc = e; |
328 | 329 |
} |
329 | 330 |
} |
330 | 331 |
return min < 0; |
331 | 332 |
} |
332 | 333 |
|
333 | 334 |
}; //class BestEligiblePivotRule |
334 | 335 |
|
335 | 336 |
|
336 | 337 |
// Implementation of the Block Search pivot rule |
337 | 338 |
class BlockSearchPivotRule |
338 | 339 |
{ |
339 | 340 |
private: |
340 | 341 |
|
341 | 342 |
// References to the NetworkSimplex class |
342 | 343 |
const IntVector &_source; |
343 | 344 |
const IntVector &_target; |
344 | 345 |
const CostVector &_cost; |
345 | 346 |
const CharVector &_state; |
346 | 347 |
const CostVector &_pi; |
347 | 348 |
int &_in_arc; |
348 | 349 |
int _search_arc_num; |
349 | 350 |
|
350 | 351 |
// Pivot rule data |
351 | 352 |
int _block_size; |
352 | 353 |
int _next_arc; |
353 | 354 |
|
354 | 355 |
public: |
355 | 356 |
|
356 | 357 |
// Constructor |
357 | 358 |
BlockSearchPivotRule(NetworkSimplex &ns) : |
358 | 359 |
_source(ns._source), _target(ns._target), |
359 | 360 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
360 | 361 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
361 | 362 |
_next_arc(0) |
362 | 363 |
{ |
363 | 364 |
// The main parameters of the pivot rule |
364 | 365 |
const double BLOCK_SIZE_FACTOR = 1.0; |
365 | 366 |
const int MIN_BLOCK_SIZE = 10; |
366 | 367 |
|
367 | 368 |
_block_size = std::max( int(BLOCK_SIZE_FACTOR * |
368 | 369 |
std::sqrt(double(_search_arc_num))), |
369 | 370 |
MIN_BLOCK_SIZE ); |
370 | 371 |
} |
371 | 372 |
|
372 | 373 |
// Find next entering arc |
373 | 374 |
bool findEnteringArc() { |
374 | 375 |
Cost c, min = 0; |
375 | 376 |
int cnt = _block_size; |
376 | 377 |
int e; |
377 | 378 |
for (e = _next_arc; e != _search_arc_num; ++e) { |
378 | 379 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
379 | 380 |
if (c < min) { |
380 | 381 |
min = c; |
381 | 382 |
_in_arc = e; |
382 | 383 |
} |
383 | 384 |
if (--cnt == 0) { |
384 | 385 |
if (min < 0) goto search_end; |
385 | 386 |
cnt = _block_size; |
386 | 387 |
} |
387 | 388 |
} |
388 | 389 |
for (e = 0; e != _next_arc; ++e) { |
389 | 390 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
390 | 391 |
if (c < min) { |
391 | 392 |
min = c; |
392 | 393 |
_in_arc = e; |
393 | 394 |
} |
394 | 395 |
if (--cnt == 0) { |
395 | 396 |
if (min < 0) goto search_end; |
396 | 397 |
cnt = _block_size; |
397 | 398 |
} |
398 | 399 |
} |
399 | 400 |
if (min >= 0) return false; |
400 | 401 |
|
401 | 402 |
search_end: |
402 | 403 |
_next_arc = e; |
403 | 404 |
return true; |
404 | 405 |
} |
405 | 406 |
|
406 | 407 |
}; //class BlockSearchPivotRule |
407 | 408 |
|
408 | 409 |
|
409 | 410 |
// Implementation of the Candidate List pivot rule |
410 | 411 |
class CandidateListPivotRule |
411 | 412 |
{ |
412 | 413 |
private: |
413 | 414 |
|
414 | 415 |
// References to the NetworkSimplex class |
415 | 416 |
const IntVector &_source; |
416 | 417 |
const IntVector &_target; |
417 | 418 |
const CostVector &_cost; |
418 | 419 |
const CharVector &_state; |
419 | 420 |
const CostVector &_pi; |
420 | 421 |
int &_in_arc; |
421 | 422 |
int _search_arc_num; |
422 | 423 |
|
423 | 424 |
// Pivot rule data |
424 | 425 |
IntVector _candidates; |
425 | 426 |
int _list_length, _minor_limit; |
426 | 427 |
int _curr_length, _minor_count; |
427 | 428 |
int _next_arc; |
428 | 429 |
|
429 | 430 |
public: |
430 | 431 |
|
431 | 432 |
/// Constructor |
432 | 433 |
CandidateListPivotRule(NetworkSimplex &ns) : |
433 | 434 |
_source(ns._source), _target(ns._target), |
434 | 435 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
435 | 436 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
436 | 437 |
_next_arc(0) |
437 | 438 |
{ |
438 | 439 |
// The main parameters of the pivot rule |
439 | 440 |
const double LIST_LENGTH_FACTOR = 0.25; |
440 | 441 |
const int MIN_LIST_LENGTH = 10; |
441 | 442 |
const double MINOR_LIMIT_FACTOR = 0.1; |
442 | 443 |
const int MIN_MINOR_LIMIT = 3; |
443 | 444 |
|
444 | 445 |
_list_length = std::max( int(LIST_LENGTH_FACTOR * |
445 | 446 |
std::sqrt(double(_search_arc_num))), |
446 | 447 |
MIN_LIST_LENGTH ); |
447 | 448 |
_minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length), |
448 | 449 |
MIN_MINOR_LIMIT ); |
449 | 450 |
_curr_length = _minor_count = 0; |
450 | 451 |
_candidates.resize(_list_length); |
0 comments (0 inline)