0
3
0
1 | 1 |
/* -*- C++ -*- |
2 | 2 |
* |
3 | 3 |
* This file is a part of LEMON, a generic C++ optimization library |
4 | 4 |
* |
5 | 5 |
* Copyright (C) 2003-2008 |
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 |
/// \tparam V The |
|
38 |
/// \tparam V The number type used for flow amounts, capacity bounds |
|
39 | 39 |
/// and supply values. By default it is \c int. |
40 |
/// \tparam C The |
|
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". It is an efficient dual |
70 | 70 |
/// solution method. |
71 | 71 |
/// |
72 | 72 |
/// Most of the parameters of the problem (except for the digraph) |
73 | 73 |
/// can be given using separate functions, and the algorithm can be |
74 | 74 |
/// executed using the \ref run() function. If some parameters are not |
75 | 75 |
/// specified, then default values will be used. |
76 | 76 |
/// |
77 | 77 |
/// \tparam GR The digraph type the algorithm runs on. |
78 |
/// \tparam V The |
|
78 |
/// \tparam V The number type used for flow amounts, capacity bounds |
|
79 | 79 |
/// and supply values in the algorithm. By default it is \c int. |
80 |
/// \tparam C The |
|
80 |
/// \tparam C The number type used for costs and potentials in the |
|
81 | 81 |
/// algorithm. By default it is the same as \c V. |
82 | 82 |
/// |
83 |
/// \warning Both |
|
83 |
/// \warning Both number types must be signed and all input data must |
|
84 | 84 |
/// be integer. |
85 | 85 |
/// \warning This algorithm does not support negative costs for such |
86 | 86 |
/// arcs that have infinite upper bound. |
87 | 87 |
#ifdef DOXYGEN |
88 | 88 |
template <typename GR, typename V, typename C, typename TR> |
89 | 89 |
#else |
90 | 90 |
template < typename GR, typename V = int, typename C = V, |
91 | 91 |
typename TR = CapacityScalingDefaultTraits<GR, V, C> > |
92 | 92 |
#endif |
93 | 93 |
class CapacityScaling |
94 | 94 |
{ |
95 | 95 |
public: |
96 | 96 |
|
97 | 97 |
/// The type of the digraph |
98 | 98 |
typedef typename TR::Digraph Digraph; |
99 | 99 |
/// The type of the flow amounts, capacity bounds and supply values |
100 | 100 |
typedef typename TR::Value Value; |
101 | 101 |
/// The type of the arc costs |
102 | 102 |
typedef typename TR::Cost Cost; |
103 | 103 |
|
104 | 104 |
/// The type of the heap used for internal Dijkstra computations |
105 | 105 |
typedef typename TR::Heap Heap; |
106 | 106 |
|
107 | 107 |
/// The \ref CapacityScalingDefaultTraits "traits class" of the algorithm |
108 | 108 |
typedef TR Traits; |
109 | 109 |
|
110 | 110 |
public: |
111 | 111 |
|
112 | 112 |
/// \brief Problem type constants for the \c run() function. |
113 | 113 |
/// |
114 | 114 |
/// Enum type containing the problem type constants that can be |
115 | 115 |
/// returned by the \ref run() function of the algorithm. |
116 | 116 |
enum ProblemType { |
117 | 117 |
/// The problem has no feasible solution (flow). |
118 | 118 |
INFEASIBLE, |
119 | 119 |
/// The problem has optimal solution (i.e. it is feasible and |
120 | 120 |
/// bounded), and the algorithm has found optimal flow and node |
121 | 121 |
/// potentials (primal and dual solutions). |
122 | 122 |
OPTIMAL, |
123 | 123 |
/// The digraph contains an arc of negative cost and infinite |
124 | 124 |
/// upper bound. It means that the objective function is unbounded |
125 |
/// on that arc, however note that it could actually be bounded |
|
125 |
/// on that arc, however, note that it could actually be bounded |
|
126 | 126 |
/// over the feasible flows, but this algroithm cannot handle |
127 | 127 |
/// these cases. |
128 | 128 |
UNBOUNDED |
129 | 129 |
}; |
130 | 130 |
|
131 | 131 |
private: |
132 | 132 |
|
133 | 133 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
134 | 134 |
|
135 | 135 |
typedef std::vector<int> IntVector; |
136 | 136 |
typedef std::vector<char> BoolVector; |
137 | 137 |
typedef std::vector<Value> ValueVector; |
138 | 138 |
typedef std::vector<Cost> CostVector; |
139 | 139 |
|
140 | 140 |
private: |
141 | 141 |
|
142 | 142 |
// Data related to the underlying digraph |
143 | 143 |
const GR &_graph; |
144 | 144 |
int _node_num; |
145 | 145 |
int _arc_num; |
146 | 146 |
int _res_arc_num; |
147 | 147 |
int _root; |
148 | 148 |
|
149 | 149 |
// Parameters of the problem |
150 | 150 |
bool _have_lower; |
151 | 151 |
Value _sum_supply; |
152 | 152 |
|
153 | 153 |
// Data structures for storing the digraph |
154 | 154 |
IntNodeMap _node_id; |
155 | 155 |
IntArcMap _arc_idf; |
156 | 156 |
IntArcMap _arc_idb; |
157 | 157 |
IntVector _first_out; |
158 | 158 |
BoolVector _forward; |
159 | 159 |
IntVector _source; |
160 | 160 |
IntVector _target; |
161 | 161 |
IntVector _reverse; |
162 | 162 |
|
163 | 163 |
// Node and arc data |
164 | 164 |
ValueVector _lower; |
165 | 165 |
ValueVector _upper; |
166 | 166 |
CostVector _cost; |
167 | 167 |
ValueVector _supply; |
168 | 168 |
|
169 | 169 |
ValueVector _res_cap; |
170 | 170 |
CostVector _pi; |
171 | 171 |
ValueVector _excess; |
172 | 172 |
IntVector _excess_nodes; |
173 | 173 |
IntVector _deficit_nodes; |
174 | 174 |
|
175 | 175 |
Value _delta; |
176 | 176 |
int _factor; |
177 | 177 |
IntVector _pred; |
178 | 178 |
|
179 | 179 |
public: |
180 | 180 |
|
181 | 181 |
/// \brief Constant for infinite upper bounds (capacities). |
182 | 182 |
/// |
183 | 183 |
/// Constant for infinite upper bounds (capacities). |
184 | 184 |
/// It is \c std::numeric_limits<Value>::infinity() if available, |
185 | 185 |
/// \c std::numeric_limits<Value>::max() otherwise. |
186 | 186 |
const Value INF; |
187 | 187 |
|
188 | 188 |
private: |
189 | 189 |
|
190 | 190 |
// Special implementation of the Dijkstra algorithm for finding |
191 | 191 |
// shortest paths in the residual network of the digraph with |
192 | 192 |
// respect to the reduced arc costs and modifying the node |
193 | 193 |
// potentials according to the found distance labels. |
194 | 194 |
class ResidualDijkstra |
195 | 195 |
{ |
196 | 196 |
private: |
197 | 197 |
|
198 | 198 |
int _node_num; |
199 | 199 |
bool _geq; |
200 | 200 |
const IntVector &_first_out; |
201 | 201 |
const IntVector &_target; |
202 | 202 |
const CostVector &_cost; |
203 | 203 |
const ValueVector &_res_cap; |
204 | 204 |
const ValueVector &_excess; |
205 | 205 |
CostVector &_pi; |
206 | 206 |
IntVector &_pred; |
207 | 207 |
|
208 | 208 |
IntVector _proc_nodes; |
209 | 209 |
CostVector _dist; |
210 | 210 |
|
211 | 211 |
public: |
212 | 212 |
|
213 | 213 |
ResidualDijkstra(CapacityScaling& cs) : |
214 | 214 |
_node_num(cs._node_num), _geq(cs._sum_supply < 0), |
215 | 215 |
_first_out(cs._first_out), _target(cs._target), _cost(cs._cost), |
216 | 216 |
_res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi), |
217 | 217 |
_pred(cs._pred), _dist(cs._node_num) |
218 | 218 |
{} |
219 | 219 |
|
220 | 220 |
int run(int s, Value delta = 1) { |
221 | 221 |
RangeMap<int> heap_cross_ref(_node_num, Heap::PRE_HEAP); |
222 | 222 |
Heap heap(heap_cross_ref); |
223 | 223 |
heap.push(s, 0); |
224 | 224 |
_pred[s] = -1; |
225 | 225 |
_proc_nodes.clear(); |
226 | 226 |
|
227 | 227 |
// Process nodes |
228 | 228 |
while (!heap.empty() && _excess[heap.top()] > -delta) { |
229 | 229 |
int u = heap.top(), v; |
230 | 230 |
Cost d = heap.prio() + _pi[u], dn; |
231 | 231 |
_dist[u] = heap.prio(); |
232 | 232 |
_proc_nodes.push_back(u); |
233 | 233 |
heap.pop(); |
234 | 234 |
|
235 | 235 |
// Traverse outgoing residual arcs |
236 | 236 |
int last_out = _geq ? _first_out[u+1] : _first_out[u+1] - 1; |
237 | 237 |
for (int a = _first_out[u]; a != last_out; ++a) { |
238 | 238 |
if (_res_cap[a] < delta) continue; |
239 | 239 |
v = _target[a]; |
240 | 240 |
switch (heap.state(v)) { |
241 | 241 |
case Heap::PRE_HEAP: |
242 | 242 |
heap.push(v, d + _cost[a] - _pi[v]); |
243 | 243 |
_pred[v] = a; |
244 | 244 |
break; |
245 | 245 |
case Heap::IN_HEAP: |
246 | 246 |
dn = d + _cost[a] - _pi[v]; |
247 | 247 |
if (dn < heap[v]) { |
248 | 248 |
heap.decrease(v, dn); |
249 | 249 |
_pred[v] = a; |
250 | 250 |
} |
251 | 251 |
break; |
252 | 252 |
case Heap::POST_HEAP: |
253 | 253 |
break; |
254 | 254 |
} |
255 | 255 |
} |
256 | 256 |
} |
257 | 257 |
if (heap.empty()) return -1; |
258 | 258 |
|
259 | 259 |
// Update potentials of processed nodes |
260 | 260 |
int t = heap.top(); |
261 | 261 |
Cost dt = heap.prio(); |
262 | 262 |
for (int i = 0; i < int(_proc_nodes.size()); ++i) { |
263 | 263 |
_pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt; |
264 | 264 |
} |
265 | 265 |
|
266 | 266 |
return t; |
267 | 267 |
} |
268 | 268 |
|
269 | 269 |
}; //class ResidualDijkstra |
270 | 270 |
|
271 | 271 |
public: |
272 | 272 |
|
273 | 273 |
/// \name Named Template Parameters |
274 | 274 |
/// @{ |
275 | 275 |
|
276 | 276 |
template <typename T> |
277 | 277 |
struct SetHeapTraits : public Traits { |
278 | 278 |
typedef T Heap; |
279 | 279 |
}; |
280 | 280 |
|
281 | 281 |
/// \brief \ref named-templ-param "Named parameter" for setting |
282 | 282 |
/// \c Heap type. |
283 | 283 |
/// |
284 | 284 |
/// \ref named-templ-param "Named parameter" for setting \c Heap |
285 | 285 |
/// type, which is used for internal Dijkstra computations. |
286 | 286 |
/// It must conform to the \ref lemon::concepts::Heap "Heap" concept, |
287 | 287 |
/// its priority type must be \c Cost and its cross reference type |
288 | 288 |
/// must be \ref RangeMap "RangeMap<int>". |
289 | 289 |
template <typename T> |
290 | 290 |
struct SetHeap |
291 | 291 |
: public CapacityScaling<GR, V, C, SetHeapTraits<T> > { |
292 | 292 |
typedef CapacityScaling<GR, V, C, SetHeapTraits<T> > Create; |
293 | 293 |
}; |
294 | 294 |
|
295 | 295 |
/// @} |
296 | 296 |
|
297 | 297 |
public: |
298 | 298 |
|
299 | 299 |
/// \brief Constructor. |
300 | 300 |
/// |
301 | 301 |
/// The constructor of the class. |
302 | 302 |
/// |
303 | 303 |
/// \param graph The digraph the algorithm runs on. |
304 | 304 |
CapacityScaling(const GR& graph) : |
305 | 305 |
_graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), |
306 | 306 |
INF(std::numeric_limits<Value>::has_infinity ? |
307 | 307 |
std::numeric_limits<Value>::infinity() : |
308 | 308 |
std::numeric_limits<Value>::max()) |
309 | 309 |
{ |
310 |
// Check the |
|
310 |
// Check the number types |
|
311 | 311 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
312 | 312 |
"The flow type of CapacityScaling must be signed"); |
313 | 313 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
314 | 314 |
"The cost type of CapacityScaling must be signed"); |
315 | 315 |
|
316 | 316 |
// Resize vectors |
317 | 317 |
_node_num = countNodes(_graph); |
318 | 318 |
_arc_num = countArcs(_graph); |
319 | 319 |
_res_arc_num = 2 * (_arc_num + _node_num); |
320 | 320 |
_root = _node_num; |
321 | 321 |
++_node_num; |
322 | 322 |
|
323 | 323 |
_first_out.resize(_node_num + 1); |
324 | 324 |
_forward.resize(_res_arc_num); |
325 | 325 |
_source.resize(_res_arc_num); |
326 | 326 |
_target.resize(_res_arc_num); |
327 | 327 |
_reverse.resize(_res_arc_num); |
328 | 328 |
|
329 | 329 |
_lower.resize(_res_arc_num); |
330 | 330 |
_upper.resize(_res_arc_num); |
331 | 331 |
_cost.resize(_res_arc_num); |
332 | 332 |
_supply.resize(_node_num); |
333 | 333 |
|
334 | 334 |
_res_cap.resize(_res_arc_num); |
335 | 335 |
_pi.resize(_node_num); |
336 | 336 |
_excess.resize(_node_num); |
337 | 337 |
_pred.resize(_node_num); |
338 | 338 |
|
339 | 339 |
// Copy the graph |
340 | 340 |
int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1; |
341 | 341 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
342 | 342 |
_node_id[n] = i; |
343 | 343 |
} |
344 | 344 |
i = 0; |
345 | 345 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
346 | 346 |
_first_out[i] = j; |
347 | 347 |
for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
348 | 348 |
_arc_idf[a] = j; |
349 | 349 |
_forward[j] = true; |
350 | 350 |
_source[j] = i; |
351 | 351 |
_target[j] = _node_id[_graph.runningNode(a)]; |
352 | 352 |
} |
353 | 353 |
for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
354 | 354 |
_arc_idb[a] = j; |
355 | 355 |
_forward[j] = false; |
356 | 356 |
_source[j] = i; |
357 | 357 |
_target[j] = _node_id[_graph.runningNode(a)]; |
358 | 358 |
} |
359 | 359 |
_forward[j] = false; |
360 | 360 |
_source[j] = i; |
361 | 361 |
_target[j] = _root; |
362 | 362 |
_reverse[j] = k; |
363 | 363 |
_forward[k] = true; |
364 | 364 |
_source[k] = _root; |
365 | 365 |
_target[k] = i; |
366 | 366 |
_reverse[k] = j; |
367 | 367 |
++j; ++k; |
368 | 368 |
} |
369 | 369 |
_first_out[i] = j; |
370 | 370 |
_first_out[_node_num] = k; |
371 | 371 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
372 | 372 |
int fi = _arc_idf[a]; |
373 | 373 |
int bi = _arc_idb[a]; |
374 | 374 |
_reverse[fi] = bi; |
375 | 375 |
_reverse[bi] = fi; |
376 | 376 |
} |
377 | 377 |
|
378 | 378 |
// Reset parameters |
379 | 379 |
reset(); |
380 | 380 |
} |
381 | 381 |
|
382 | 382 |
/// \name Parameters |
383 | 383 |
/// The parameters of the algorithm can be specified using these |
384 | 384 |
/// functions. |
385 | 385 |
|
386 | 386 |
/// @{ |
387 | 387 |
|
388 | 388 |
/// \brief Set the lower bounds on the arcs. |
389 | 389 |
/// |
390 | 390 |
/// This function sets the lower bounds on the arcs. |
391 | 391 |
/// If it is not used before calling \ref run(), the lower bounds |
392 | 392 |
/// will be set to zero on all arcs. |
393 | 393 |
/// |
394 | 394 |
/// \param map An arc map storing the lower bounds. |
395 | 395 |
/// Its \c Value type must be convertible to the \c Value type |
396 | 396 |
/// of the algorithm. |
397 | 397 |
/// |
398 | 398 |
/// \return <tt>(*this)</tt> |
399 | 399 |
template <typename LowerMap> |
400 | 400 |
CapacityScaling& lowerMap(const LowerMap& map) { |
401 | 401 |
_have_lower = true; |
402 | 402 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
403 | 403 |
_lower[_arc_idf[a]] = map[a]; |
404 | 404 |
_lower[_arc_idb[a]] = map[a]; |
405 | 405 |
} |
406 | 406 |
return *this; |
407 | 407 |
} |
408 | 408 |
|
409 | 409 |
/// \brief Set the upper bounds (capacities) on the arcs. |
410 | 410 |
/// |
411 | 411 |
/// This function sets the upper bounds (capacities) on the arcs. |
412 | 412 |
/// If it is not used before calling \ref run(), the upper bounds |
413 | 413 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be |
414 |
/// unbounded from above |
|
414 |
/// unbounded from above). |
|
415 | 415 |
/// |
416 | 416 |
/// \param map An arc map storing the upper bounds. |
417 | 417 |
/// Its \c Value type must be convertible to the \c Value type |
418 | 418 |
/// of the algorithm. |
419 | 419 |
/// |
420 | 420 |
/// \return <tt>(*this)</tt> |
421 | 421 |
template<typename UpperMap> |
422 | 422 |
CapacityScaling& upperMap(const UpperMap& map) { |
423 | 423 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
424 | 424 |
_upper[_arc_idf[a]] = map[a]; |
425 | 425 |
} |
426 | 426 |
return *this; |
427 | 427 |
} |
428 | 428 |
|
429 | 429 |
/// \brief Set the costs of the arcs. |
430 | 430 |
/// |
431 | 431 |
/// This function sets the costs of the arcs. |
432 | 432 |
/// If it is not used before calling \ref run(), the costs |
433 | 433 |
/// will be set to \c 1 on all arcs. |
434 | 434 |
/// |
435 | 435 |
/// \param map An arc map storing the costs. |
436 | 436 |
/// Its \c Value type must be convertible to the \c Cost type |
437 | 437 |
/// of the algorithm. |
438 | 438 |
/// |
439 | 439 |
/// \return <tt>(*this)</tt> |
440 | 440 |
template<typename CostMap> |
441 | 441 |
CapacityScaling& costMap(const CostMap& map) { |
442 | 442 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
443 | 443 |
_cost[_arc_idf[a]] = map[a]; |
444 | 444 |
_cost[_arc_idb[a]] = -map[a]; |
445 | 445 |
} |
446 | 446 |
return *this; |
447 | 447 |
} |
448 | 448 |
|
449 | 449 |
/// \brief Set the supply values of the nodes. |
450 | 450 |
/// |
451 | 451 |
/// This function sets the supply values of the nodes. |
452 | 452 |
/// If neither this function nor \ref stSupply() is used before |
453 | 453 |
/// calling \ref run(), the supply of each node will be set to zero. |
454 | 454 |
/// |
455 | 455 |
/// \param map A node map storing the supply values. |
456 | 456 |
/// Its \c Value type must be convertible to the \c Value type |
457 | 457 |
/// of the algorithm. |
458 | 458 |
/// |
459 | 459 |
/// \return <tt>(*this)</tt> |
460 | 460 |
template<typename SupplyMap> |
461 | 461 |
CapacityScaling& supplyMap(const SupplyMap& map) { |
462 | 462 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
463 | 463 |
_supply[_node_id[n]] = map[n]; |
464 | 464 |
} |
465 | 465 |
return *this; |
466 | 466 |
} |
467 | 467 |
|
468 | 468 |
/// \brief Set single source and target nodes and a supply value. |
469 | 469 |
/// |
470 | 470 |
/// This function sets a single source node and a single target node |
471 | 471 |
/// and the required flow value. |
472 | 472 |
/// If neither this function nor \ref supplyMap() is used before |
473 | 473 |
/// calling \ref run(), the supply of each node will be set to zero. |
474 | 474 |
/// |
475 | 475 |
/// Using this function has the same effect as using \ref supplyMap() |
476 | 476 |
/// with such a map in which \c k is assigned to \c s, \c -k is |
477 | 477 |
/// assigned to \c t and all other nodes have zero supply value. |
478 | 478 |
/// |
479 | 479 |
/// \param s The source node. |
480 | 480 |
/// \param t The target node. |
481 | 481 |
/// \param k The required amount of flow from node \c s to node \c t |
482 | 482 |
/// (i.e. the supply of \c s and the demand of \c t). |
483 | 483 |
/// |
484 | 484 |
/// \return <tt>(*this)</tt> |
485 | 485 |
CapacityScaling& stSupply(const Node& s, const Node& t, Value k) { |
486 | 486 |
for (int i = 0; i != _node_num; ++i) { |
487 | 487 |
_supply[i] = 0; |
488 | 488 |
} |
489 | 489 |
_supply[_node_id[s]] = k; |
490 | 490 |
_supply[_node_id[t]] = -k; |
491 | 491 |
return *this; |
492 | 492 |
} |
493 | 493 |
|
494 | 494 |
/// @} |
495 | 495 |
|
496 | 496 |
/// \name Execution control |
497 | 497 |
/// The algorithm can be executed using \ref run(). |
498 | 498 |
|
499 | 499 |
/// @{ |
500 | 500 |
|
501 | 501 |
/// \brief Run the algorithm. |
502 | 502 |
/// |
503 | 503 |
/// This function runs the algorithm. |
504 | 504 |
/// The paramters can be specified using functions \ref lowerMap(), |
505 | 505 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). |
506 | 506 |
/// For example, |
507 | 507 |
/// \code |
508 | 508 |
/// CapacityScaling<ListDigraph> cs(graph); |
509 | 509 |
/// cs.lowerMap(lower).upperMap(upper).costMap(cost) |
510 | 510 |
/// .supplyMap(sup).run(); |
511 | 511 |
/// \endcode |
512 | 512 |
/// |
513 | 513 |
/// This function can be called more than once. All the parameters |
514 | 514 |
/// that have been given are kept for the next call, unless |
515 | 515 |
/// \ref reset() is called, thus only the modified parameters |
516 | 516 |
/// have to be set again. See \ref reset() for examples. |
517 |
/// However the underlying digraph must not be modified after this |
|
517 |
/// However, the underlying digraph must not be modified after this |
|
518 | 518 |
/// class have been constructed, since it copies and extends the graph. |
519 | 519 |
/// |
520 | 520 |
/// \param factor The capacity scaling factor. It must be larger than |
521 | 521 |
/// one to use scaling. If it is less or equal to one, then scaling |
522 | 522 |
/// will be disabled. |
523 | 523 |
/// |
524 | 524 |
/// \return \c INFEASIBLE if no feasible flow exists, |
525 | 525 |
/// \n \c OPTIMAL if the problem has optimal solution |
526 | 526 |
/// (i.e. it is feasible and bounded), and the algorithm has found |
527 | 527 |
/// optimal flow and node potentials (primal and dual solutions), |
528 | 528 |
/// \n \c UNBOUNDED if the digraph contains an arc of negative cost |
529 | 529 |
/// and infinite upper bound. It means that the objective function |
530 |
/// is unbounded on that arc, however note that it could actually be |
|
530 |
/// is unbounded on that arc, however, note that it could actually be |
|
531 | 531 |
/// bounded over the feasible flows, but this algroithm cannot handle |
532 | 532 |
/// these cases. |
533 | 533 |
/// |
534 | 534 |
/// \see ProblemType |
535 | 535 |
ProblemType run(int factor = 4) { |
536 | 536 |
_factor = factor; |
537 | 537 |
ProblemType pt = init(); |
538 | 538 |
if (pt != OPTIMAL) return pt; |
539 | 539 |
return start(); |
540 | 540 |
} |
541 | 541 |
|
542 | 542 |
/// \brief Reset all the parameters that have been given before. |
543 | 543 |
/// |
544 | 544 |
/// This function resets all the paramaters that have been given |
545 | 545 |
/// before using functions \ref lowerMap(), \ref upperMap(), |
546 | 546 |
/// \ref costMap(), \ref supplyMap(), \ref stSupply(). |
547 | 547 |
/// |
548 | 548 |
/// It is useful for multiple run() calls. If this function is not |
549 | 549 |
/// used, all the parameters given before are kept for the next |
550 | 550 |
/// \ref run() call. |
551 | 551 |
/// However, the underlying digraph must not be modified after this |
552 | 552 |
/// class have been constructed, since it copies and extends the graph. |
553 | 553 |
/// |
554 | 554 |
/// For example, |
555 | 555 |
/// \code |
556 | 556 |
/// CapacityScaling<ListDigraph> cs(graph); |
557 | 557 |
/// |
558 | 558 |
/// // First run |
559 | 559 |
/// cs.lowerMap(lower).upperMap(upper).costMap(cost) |
560 | 560 |
/// .supplyMap(sup).run(); |
561 | 561 |
/// |
562 | 562 |
/// // Run again with modified cost map (reset() is not called, |
563 | 563 |
/// // so only the cost map have to be set again) |
564 | 564 |
/// cost[e] += 100; |
565 | 565 |
/// cs.costMap(cost).run(); |
566 | 566 |
/// |
567 | 567 |
/// // Run again from scratch using reset() |
568 | 568 |
/// // (the lower bounds will be set to zero on all arcs) |
569 | 569 |
/// cs.reset(); |
570 | 570 |
/// cs.upperMap(capacity).costMap(cost) |
571 | 571 |
/// .supplyMap(sup).run(); |
572 | 572 |
/// \endcode |
573 | 573 |
/// |
574 | 574 |
/// \return <tt>(*this)</tt> |
575 | 575 |
CapacityScaling& reset() { |
576 | 576 |
for (int i = 0; i != _node_num; ++i) { |
577 | 577 |
_supply[i] = 0; |
578 | 578 |
} |
579 | 579 |
for (int j = 0; j != _res_arc_num; ++j) { |
580 | 580 |
_lower[j] = 0; |
581 | 581 |
_upper[j] = INF; |
582 | 582 |
_cost[j] = _forward[j] ? 1 : -1; |
583 | 583 |
} |
584 | 584 |
_have_lower = false; |
585 | 585 |
return *this; |
586 | 586 |
} |
587 | 587 |
|
588 | 588 |
/// @} |
589 | 589 |
|
590 | 590 |
/// \name Query Functions |
591 | 591 |
/// The results of the algorithm can be obtained using these |
592 | 592 |
/// functions.\n |
593 | 593 |
/// The \ref run() function must be called before using them. |
594 | 594 |
|
595 | 595 |
/// @{ |
596 | 596 |
|
597 | 597 |
/// \brief Return the total cost of the found flow. |
598 | 598 |
/// |
599 | 599 |
/// This function returns the total cost of the found flow. |
600 | 600 |
/// Its complexity is O(e). |
601 | 601 |
/// |
602 | 602 |
/// \note The return type of the function can be specified as a |
603 | 603 |
/// template parameter. For example, |
604 | 604 |
/// \code |
605 | 605 |
/// cs.totalCost<double>(); |
606 | 606 |
/// \endcode |
607 | 607 |
/// It is useful if the total cost cannot be stored in the \c Cost |
608 | 608 |
/// type of the algorithm, which is the default return type of the |
609 | 609 |
/// function. |
610 | 610 |
/// |
611 | 611 |
/// \pre \ref run() must be called before using this function. |
612 | 612 |
template <typename Number> |
613 | 613 |
Number totalCost() const { |
614 | 614 |
Number c = 0; |
615 | 615 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
616 | 616 |
int i = _arc_idb[a]; |
617 | 617 |
c += static_cast<Number>(_res_cap[i]) * |
618 | 618 |
(-static_cast<Number>(_cost[i])); |
619 | 619 |
} |
620 | 620 |
return c; |
621 | 621 |
} |
622 | 622 |
|
623 | 623 |
#ifndef DOXYGEN |
624 | 624 |
Cost totalCost() const { |
625 | 625 |
return totalCost<Cost>(); |
626 | 626 |
} |
627 | 627 |
#endif |
628 | 628 |
|
629 | 629 |
/// \brief Return the flow on the given arc. |
630 | 630 |
/// |
631 | 631 |
/// This function returns the flow on the given arc. |
632 | 632 |
/// |
633 | 633 |
/// \pre \ref run() must be called before using this function. |
634 | 634 |
Value flow(const Arc& a) const { |
635 | 635 |
return _res_cap[_arc_idb[a]]; |
636 | 636 |
} |
637 | 637 |
|
638 | 638 |
/// \brief Return the flow map (the primal solution). |
639 | 639 |
/// |
640 | 640 |
/// This function copies the flow value on each arc into the given |
641 | 641 |
/// map. The \c Value type of the algorithm must be convertible to |
642 | 642 |
/// the \c Value type of the map. |
643 | 643 |
/// |
644 | 644 |
/// \pre \ref run() must be called before using this function. |
645 | 645 |
template <typename FlowMap> |
646 | 646 |
void flowMap(FlowMap &map) const { |
647 | 647 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
648 | 648 |
map.set(a, _res_cap[_arc_idb[a]]); |
649 | 649 |
} |
650 | 650 |
} |
651 | 651 |
|
652 | 652 |
/// \brief Return the potential (dual value) of the given node. |
653 | 653 |
/// |
654 | 654 |
/// This function returns the potential (dual value) of the |
655 | 655 |
/// given node. |
656 | 656 |
/// |
657 | 657 |
/// \pre \ref run() must be called before using this function. |
658 | 658 |
Cost potential(const Node& n) const { |
659 | 659 |
return _pi[_node_id[n]]; |
660 | 660 |
} |
661 | 661 |
|
662 | 662 |
/// \brief Return the potential map (the dual solution). |
663 | 663 |
/// |
664 | 664 |
/// This function copies the potential (dual value) of each node |
665 | 665 |
/// into the given map. |
666 | 666 |
/// The \c Cost type of the algorithm must be convertible to the |
667 | 667 |
/// \c Value type of the map. |
668 | 668 |
/// |
669 | 669 |
/// \pre \ref run() must be called before using this function. |
670 | 670 |
template <typename PotentialMap> |
671 | 671 |
void potentialMap(PotentialMap &map) const { |
672 | 672 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
673 | 673 |
map.set(n, _pi[_node_id[n]]); |
674 | 674 |
} |
675 | 675 |
} |
676 | 676 |
|
677 | 677 |
/// @} |
678 | 678 |
|
679 | 679 |
private: |
680 | 680 |
|
681 | 681 |
// Initialize the algorithm |
682 | 682 |
ProblemType init() { |
683 | 683 |
if (_node_num == 0) return INFEASIBLE; |
684 | 684 |
|
685 | 685 |
// Check the sum of supply values |
686 | 686 |
_sum_supply = 0; |
687 | 687 |
for (int i = 0; i != _root; ++i) { |
688 | 688 |
_sum_supply += _supply[i]; |
689 | 689 |
} |
690 | 690 |
if (_sum_supply > 0) return INFEASIBLE; |
691 | 691 |
|
692 | 692 |
// Initialize vectors |
693 | 693 |
for (int i = 0; i != _root; ++i) { |
694 | 694 |
_pi[i] = 0; |
695 | 695 |
_excess[i] = _supply[i]; |
696 | 696 |
} |
697 | 697 |
|
698 | 698 |
// Remove non-zero lower bounds |
699 | 699 |
const Value MAX = std::numeric_limits<Value>::max(); |
700 | 700 |
int last_out; |
701 | 701 |
if (_have_lower) { |
702 | 702 |
for (int i = 0; i != _root; ++i) { |
703 | 703 |
last_out = _first_out[i+1]; |
704 | 704 |
for (int j = _first_out[i]; j != last_out; ++j) { |
705 | 705 |
if (_forward[j]) { |
706 | 706 |
Value c = _lower[j]; |
707 | 707 |
if (c >= 0) { |
708 | 708 |
_res_cap[j] = _upper[j] < MAX ? _upper[j] - c : INF; |
709 | 709 |
} else { |
710 | 710 |
_res_cap[j] = _upper[j] < MAX + c ? _upper[j] - c : INF; |
711 | 711 |
} |
712 | 712 |
_excess[i] -= c; |
713 | 713 |
_excess[_target[j]] += c; |
714 | 714 |
} else { |
715 | 715 |
_res_cap[j] = 0; |
716 | 716 |
} |
717 | 717 |
} |
718 | 718 |
} |
719 | 719 |
} else { |
720 | 720 |
for (int j = 0; j != _res_arc_num; ++j) { |
721 | 721 |
_res_cap[j] = _forward[j] ? _upper[j] : 0; |
722 | 722 |
} |
1 | 1 |
/* -*- C++ -*- |
2 | 2 |
* |
3 | 3 |
* This file is a part of LEMON, a generic C++ optimization library |
4 | 4 |
* |
5 | 5 |
* Copyright (C) 2003-2008 |
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 |
/// \tparam V The |
|
43 |
/// \tparam V The number type used for flow amounts, capacity bounds |
|
44 | 44 |
/// and supply values. By default it is \c int. |
45 |
/// \tparam C The |
|
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 minimum cost |
94 | 94 |
/// flow. It is an efficient primal-dual solution method, which |
95 | 95 |
/// can be viewed as the generalization of the \ref Preflow |
96 | 96 |
/// "preflow push-relabel" algorithm for the maximum flow problem. |
97 | 97 |
/// |
98 | 98 |
/// Most of the parameters of the problem (except for the digraph) |
99 | 99 |
/// can be given using separate functions, and the algorithm can be |
100 | 100 |
/// executed using the \ref run() function. If some parameters are not |
101 | 101 |
/// specified, then default values will be used. |
102 | 102 |
/// |
103 | 103 |
/// \tparam GR The digraph type the algorithm runs on. |
104 |
/// \tparam V The |
|
104 |
/// \tparam V The number type used for flow amounts, capacity bounds |
|
105 | 105 |
/// and supply values in the algorithm. By default it is \c int. |
106 |
/// \tparam C The |
|
106 |
/// \tparam C The number type used for costs and potentials in the |
|
107 | 107 |
/// algorithm. By default it is the same as \c V. |
108 | 108 |
/// |
109 |
/// \warning Both |
|
109 |
/// \warning Both number types must be signed and all input data must |
|
110 | 110 |
/// be integer. |
111 | 111 |
/// \warning This algorithm does not support negative costs for such |
112 | 112 |
/// arcs that have infinite upper bound. |
113 | 113 |
/// |
114 | 114 |
/// \note %CostScaling provides three different internal methods, |
115 | 115 |
/// from which the most efficient one is used by default. |
116 | 116 |
/// For more information, see \ref Method. |
117 | 117 |
#ifdef DOXYGEN |
118 | 118 |
template <typename GR, typename V, typename C, typename TR> |
119 | 119 |
#else |
120 | 120 |
template < typename GR, typename V = int, typename C = V, |
121 | 121 |
typename TR = CostScalingDefaultTraits<GR, V, C> > |
122 | 122 |
#endif |
123 | 123 |
class CostScaling |
124 | 124 |
{ |
125 | 125 |
public: |
126 | 126 |
|
127 | 127 |
/// The type of the digraph |
128 | 128 |
typedef typename TR::Digraph Digraph; |
129 | 129 |
/// The type of the flow amounts, capacity bounds and supply values |
130 | 130 |
typedef typename TR::Value Value; |
131 | 131 |
/// The type of the arc costs |
132 | 132 |
typedef typename TR::Cost Cost; |
133 | 133 |
|
134 | 134 |
/// \brief The large cost type |
135 | 135 |
/// |
136 | 136 |
/// The large cost type used for internal computations. |
137 | 137 |
/// Using the \ref CostScalingDefaultTraits "default traits class", |
138 | 138 |
/// it is \c long \c long if the \c Cost type is integer, |
139 | 139 |
/// otherwise it is \c double. |
140 | 140 |
typedef typename TR::LargeCost LargeCost; |
141 | 141 |
|
142 | 142 |
/// The \ref CostScalingDefaultTraits "traits class" of the algorithm |
143 | 143 |
typedef TR Traits; |
144 | 144 |
|
145 | 145 |
public: |
146 | 146 |
|
147 | 147 |
/// \brief Problem type constants for the \c run() function. |
148 | 148 |
/// |
149 | 149 |
/// Enum type containing the problem type constants that can be |
150 | 150 |
/// returned by the \ref run() function of the algorithm. |
151 | 151 |
enum ProblemType { |
152 | 152 |
/// The problem has no feasible solution (flow). |
153 | 153 |
INFEASIBLE, |
154 | 154 |
/// The problem has optimal solution (i.e. it is feasible and |
155 | 155 |
/// bounded), and the algorithm has found optimal flow and node |
156 | 156 |
/// potentials (primal and dual solutions). |
157 | 157 |
OPTIMAL, |
158 | 158 |
/// The digraph contains an arc of negative cost and infinite |
159 | 159 |
/// upper bound. It means that the objective function is unbounded |
160 |
/// on that arc, however note that it could actually be bounded |
|
160 |
/// on that arc, however, note that it could actually be bounded |
|
161 | 161 |
/// over the feasible flows, but this algroithm cannot handle |
162 | 162 |
/// these cases. |
163 | 163 |
UNBOUNDED |
164 | 164 |
}; |
165 | 165 |
|
166 | 166 |
/// \brief Constants for selecting the internal method. |
167 | 167 |
/// |
168 | 168 |
/// Enum type containing constants for selecting the internal method |
169 | 169 |
/// for the \ref run() function. |
170 | 170 |
/// |
171 | 171 |
/// \ref CostScaling provides three internal methods that differ mainly |
172 | 172 |
/// in their base operations, which are used in conjunction with the |
173 | 173 |
/// relabel operation. |
174 | 174 |
/// By default, the so called \ref PARTIAL_AUGMENT |
175 | 175 |
/// "Partial Augment-Relabel" method is used, which proved to be |
176 | 176 |
/// the most efficient and the most robust on various test inputs. |
177 | 177 |
/// However, the other methods can be selected using the \ref run() |
178 | 178 |
/// function with the proper parameter. |
179 | 179 |
enum Method { |
180 | 180 |
/// Local push operations are used, i.e. flow is moved only on one |
181 | 181 |
/// admissible arc at once. |
182 | 182 |
PUSH, |
183 | 183 |
/// Augment operations are used, i.e. flow is moved on admissible |
184 | 184 |
/// paths from a node with excess to a node with deficit. |
185 | 185 |
AUGMENT, |
186 | 186 |
/// Partial augment operations are used, i.e. flow is moved on |
187 | 187 |
/// admissible paths started from a node with excess, but the |
188 | 188 |
/// lengths of these paths are limited. This method can be viewed |
189 | 189 |
/// as a combined version of the previous two operations. |
190 | 190 |
PARTIAL_AUGMENT |
191 | 191 |
}; |
192 | 192 |
|
193 | 193 |
private: |
194 | 194 |
|
195 | 195 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
196 | 196 |
|
197 | 197 |
typedef std::vector<int> IntVector; |
198 | 198 |
typedef std::vector<char> BoolVector; |
199 | 199 |
typedef std::vector<Value> ValueVector; |
200 | 200 |
typedef std::vector<Cost> CostVector; |
201 | 201 |
typedef std::vector<LargeCost> LargeCostVector; |
202 | 202 |
|
203 | 203 |
private: |
204 | 204 |
|
205 | 205 |
template <typename KT, typename VT> |
206 | 206 |
class VectorMap { |
207 | 207 |
public: |
208 | 208 |
typedef KT Key; |
209 | 209 |
typedef VT Value; |
210 | 210 |
|
211 | 211 |
VectorMap(std::vector<Value>& v) : _v(v) {} |
212 | 212 |
|
213 | 213 |
const Value& operator[](const Key& key) const { |
214 | 214 |
return _v[StaticDigraph::id(key)]; |
215 | 215 |
} |
216 | 216 |
|
217 | 217 |
Value& operator[](const Key& key) { |
218 | 218 |
return _v[StaticDigraph::id(key)]; |
219 | 219 |
} |
220 | 220 |
|
221 | 221 |
void set(const Key& key, const Value& val) { |
222 | 222 |
_v[StaticDigraph::id(key)] = val; |
223 | 223 |
} |
224 | 224 |
|
225 | 225 |
private: |
226 | 226 |
std::vector<Value>& _v; |
227 | 227 |
}; |
228 | 228 |
|
229 | 229 |
typedef VectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap; |
230 | 230 |
typedef VectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap; |
231 | 231 |
|
232 | 232 |
private: |
233 | 233 |
|
234 | 234 |
// Data related to the underlying digraph |
235 | 235 |
const GR &_graph; |
236 | 236 |
int _node_num; |
237 | 237 |
int _arc_num; |
238 | 238 |
int _res_node_num; |
239 | 239 |
int _res_arc_num; |
240 | 240 |
int _root; |
241 | 241 |
|
242 | 242 |
// Parameters of the problem |
243 | 243 |
bool _have_lower; |
244 | 244 |
Value _sum_supply; |
245 | 245 |
|
246 | 246 |
// Data structures for storing the digraph |
247 | 247 |
IntNodeMap _node_id; |
248 | 248 |
IntArcMap _arc_idf; |
249 | 249 |
IntArcMap _arc_idb; |
250 | 250 |
IntVector _first_out; |
251 | 251 |
BoolVector _forward; |
252 | 252 |
IntVector _source; |
253 | 253 |
IntVector _target; |
254 | 254 |
IntVector _reverse; |
255 | 255 |
|
256 | 256 |
// Node and arc data |
257 | 257 |
ValueVector _lower; |
258 | 258 |
ValueVector _upper; |
259 | 259 |
CostVector _scost; |
260 | 260 |
ValueVector _supply; |
261 | 261 |
|
262 | 262 |
ValueVector _res_cap; |
263 | 263 |
LargeCostVector _cost; |
264 | 264 |
LargeCostVector _pi; |
265 | 265 |
ValueVector _excess; |
266 | 266 |
IntVector _next_out; |
267 | 267 |
std::deque<int> _active_nodes; |
268 | 268 |
|
269 | 269 |
// Data for scaling |
270 | 270 |
LargeCost _epsilon; |
271 | 271 |
int _alpha; |
272 | 272 |
|
273 | 273 |
// Data for a StaticDigraph structure |
274 | 274 |
typedef std::pair<int, int> IntPair; |
275 | 275 |
StaticDigraph _sgr; |
276 | 276 |
std::vector<IntPair> _arc_vec; |
277 | 277 |
std::vector<LargeCost> _cost_vec; |
278 | 278 |
LargeCostArcMap _cost_map; |
279 | 279 |
LargeCostNodeMap _pi_map; |
280 | 280 |
|
281 | 281 |
public: |
282 | 282 |
|
283 | 283 |
/// \brief Constant for infinite upper bounds (capacities). |
284 | 284 |
/// |
285 | 285 |
/// Constant for infinite upper bounds (capacities). |
286 | 286 |
/// It is \c std::numeric_limits<Value>::infinity() if available, |
287 | 287 |
/// \c std::numeric_limits<Value>::max() otherwise. |
288 | 288 |
const Value INF; |
289 | 289 |
|
290 | 290 |
public: |
291 | 291 |
|
292 | 292 |
/// \name Named Template Parameters |
293 | 293 |
/// @{ |
294 | 294 |
|
295 | 295 |
template <typename T> |
296 | 296 |
struct SetLargeCostTraits : public Traits { |
297 | 297 |
typedef T LargeCost; |
298 | 298 |
}; |
299 | 299 |
|
300 | 300 |
/// \brief \ref named-templ-param "Named parameter" for setting |
301 | 301 |
/// \c LargeCost type. |
302 | 302 |
/// |
303 | 303 |
/// \ref named-templ-param "Named parameter" for setting \c LargeCost |
304 | 304 |
/// type, which is used for internal computations in the algorithm. |
305 | 305 |
/// \c Cost must be convertible to \c LargeCost. |
306 | 306 |
template <typename T> |
307 | 307 |
struct SetLargeCost |
308 | 308 |
: public CostScaling<GR, V, C, SetLargeCostTraits<T> > { |
309 | 309 |
typedef CostScaling<GR, V, C, SetLargeCostTraits<T> > Create; |
310 | 310 |
}; |
311 | 311 |
|
312 | 312 |
/// @} |
313 | 313 |
|
314 | 314 |
public: |
315 | 315 |
|
316 | 316 |
/// \brief Constructor. |
317 | 317 |
/// |
318 | 318 |
/// The constructor of the class. |
319 | 319 |
/// |
320 | 320 |
/// \param graph The digraph the algorithm runs on. |
321 | 321 |
CostScaling(const GR& graph) : |
322 | 322 |
_graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), |
323 | 323 |
_cost_map(_cost_vec), _pi_map(_pi), |
324 | 324 |
INF(std::numeric_limits<Value>::has_infinity ? |
325 | 325 |
std::numeric_limits<Value>::infinity() : |
326 | 326 |
std::numeric_limits<Value>::max()) |
327 | 327 |
{ |
328 |
// Check the |
|
328 |
// Check the number types |
|
329 | 329 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
330 | 330 |
"The flow type of CostScaling must be signed"); |
331 | 331 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
332 | 332 |
"The cost type of CostScaling must be signed"); |
333 | 333 |
|
334 | 334 |
// Resize vectors |
335 | 335 |
_node_num = countNodes(_graph); |
336 | 336 |
_arc_num = countArcs(_graph); |
337 | 337 |
_res_node_num = _node_num + 1; |
338 | 338 |
_res_arc_num = 2 * (_arc_num + _node_num); |
339 | 339 |
_root = _node_num; |
340 | 340 |
|
341 | 341 |
_first_out.resize(_res_node_num + 1); |
342 | 342 |
_forward.resize(_res_arc_num); |
343 | 343 |
_source.resize(_res_arc_num); |
344 | 344 |
_target.resize(_res_arc_num); |
345 | 345 |
_reverse.resize(_res_arc_num); |
346 | 346 |
|
347 | 347 |
_lower.resize(_res_arc_num); |
348 | 348 |
_upper.resize(_res_arc_num); |
349 | 349 |
_scost.resize(_res_arc_num); |
350 | 350 |
_supply.resize(_res_node_num); |
351 | 351 |
|
352 | 352 |
_res_cap.resize(_res_arc_num); |
353 | 353 |
_cost.resize(_res_arc_num); |
354 | 354 |
_pi.resize(_res_node_num); |
355 | 355 |
_excess.resize(_res_node_num); |
356 | 356 |
_next_out.resize(_res_node_num); |
357 | 357 |
|
358 | 358 |
_arc_vec.reserve(_res_arc_num); |
359 | 359 |
_cost_vec.reserve(_res_arc_num); |
360 | 360 |
|
361 | 361 |
// Copy the graph |
362 | 362 |
int i = 0, j = 0, k = 2 * _arc_num + _node_num; |
363 | 363 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
364 | 364 |
_node_id[n] = i; |
365 | 365 |
} |
366 | 366 |
i = 0; |
367 | 367 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
368 | 368 |
_first_out[i] = j; |
369 | 369 |
for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
370 | 370 |
_arc_idf[a] = j; |
371 | 371 |
_forward[j] = true; |
372 | 372 |
_source[j] = i; |
373 | 373 |
_target[j] = _node_id[_graph.runningNode(a)]; |
374 | 374 |
} |
375 | 375 |
for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
376 | 376 |
_arc_idb[a] = j; |
377 | 377 |
_forward[j] = false; |
378 | 378 |
_source[j] = i; |
379 | 379 |
_target[j] = _node_id[_graph.runningNode(a)]; |
380 | 380 |
} |
381 | 381 |
_forward[j] = false; |
382 | 382 |
_source[j] = i; |
383 | 383 |
_target[j] = _root; |
384 | 384 |
_reverse[j] = k; |
385 | 385 |
_forward[k] = true; |
386 | 386 |
_source[k] = _root; |
387 | 387 |
_target[k] = i; |
388 | 388 |
_reverse[k] = j; |
389 | 389 |
++j; ++k; |
390 | 390 |
} |
391 | 391 |
_first_out[i] = j; |
392 | 392 |
_first_out[_res_node_num] = k; |
393 | 393 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
394 | 394 |
int fi = _arc_idf[a]; |
395 | 395 |
int bi = _arc_idb[a]; |
396 | 396 |
_reverse[fi] = bi; |
397 | 397 |
_reverse[bi] = fi; |
398 | 398 |
} |
399 | 399 |
|
400 | 400 |
// Reset parameters |
401 | 401 |
reset(); |
402 | 402 |
} |
403 | 403 |
|
404 | 404 |
/// \name Parameters |
405 | 405 |
/// The parameters of the algorithm can be specified using these |
406 | 406 |
/// functions. |
407 | 407 |
|
408 | 408 |
/// @{ |
409 | 409 |
|
410 | 410 |
/// \brief Set the lower bounds on the arcs. |
411 | 411 |
/// |
412 | 412 |
/// This function sets the lower bounds on the arcs. |
413 | 413 |
/// If it is not used before calling \ref run(), the lower bounds |
414 | 414 |
/// will be set to zero on all arcs. |
415 | 415 |
/// |
416 | 416 |
/// \param map An arc map storing the lower bounds. |
417 | 417 |
/// Its \c Value type must be convertible to the \c Value type |
418 | 418 |
/// of the algorithm. |
419 | 419 |
/// |
420 | 420 |
/// \return <tt>(*this)</tt> |
421 | 421 |
template <typename LowerMap> |
422 | 422 |
CostScaling& lowerMap(const LowerMap& map) { |
423 | 423 |
_have_lower = true; |
424 | 424 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
425 | 425 |
_lower[_arc_idf[a]] = map[a]; |
426 | 426 |
_lower[_arc_idb[a]] = map[a]; |
427 | 427 |
} |
428 | 428 |
return *this; |
429 | 429 |
} |
430 | 430 |
|
431 | 431 |
/// \brief Set the upper bounds (capacities) on the arcs. |
432 | 432 |
/// |
433 | 433 |
/// This function sets the upper bounds (capacities) on the arcs. |
434 | 434 |
/// If it is not used before calling \ref run(), the upper bounds |
435 | 435 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be |
436 |
/// unbounded from above |
|
436 |
/// unbounded from above). |
|
437 | 437 |
/// |
438 | 438 |
/// \param map An arc map storing the upper bounds. |
439 | 439 |
/// Its \c Value type must be convertible to the \c Value type |
440 | 440 |
/// of the algorithm. |
441 | 441 |
/// |
442 | 442 |
/// \return <tt>(*this)</tt> |
443 | 443 |
template<typename UpperMap> |
444 | 444 |
CostScaling& upperMap(const UpperMap& map) { |
445 | 445 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
446 | 446 |
_upper[_arc_idf[a]] = map[a]; |
447 | 447 |
} |
448 | 448 |
return *this; |
449 | 449 |
} |
450 | 450 |
|
451 | 451 |
/// \brief Set the costs of the arcs. |
452 | 452 |
/// |
453 | 453 |
/// This function sets the costs of the arcs. |
454 | 454 |
/// If it is not used before calling \ref run(), the costs |
455 | 455 |
/// will be set to \c 1 on all arcs. |
456 | 456 |
/// |
457 | 457 |
/// \param map An arc map storing the costs. |
458 | 458 |
/// Its \c Value type must be convertible to the \c Cost type |
459 | 459 |
/// of the algorithm. |
460 | 460 |
/// |
461 | 461 |
/// \return <tt>(*this)</tt> |
462 | 462 |
template<typename CostMap> |
463 | 463 |
CostScaling& costMap(const CostMap& map) { |
464 | 464 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
465 | 465 |
_scost[_arc_idf[a]] = map[a]; |
466 | 466 |
_scost[_arc_idb[a]] = -map[a]; |
467 | 467 |
} |
468 | 468 |
return *this; |
469 | 469 |
} |
470 | 470 |
|
471 | 471 |
/// \brief Set the supply values of the nodes. |
472 | 472 |
/// |
473 | 473 |
/// This function sets the supply values of the nodes. |
474 | 474 |
/// If neither this function nor \ref stSupply() is used before |
475 | 475 |
/// calling \ref run(), the supply of each node will be set to zero. |
476 | 476 |
/// |
477 | 477 |
/// \param map A node map storing the supply values. |
478 | 478 |
/// Its \c Value type must be convertible to the \c Value type |
479 | 479 |
/// of the algorithm. |
480 | 480 |
/// |
481 | 481 |
/// \return <tt>(*this)</tt> |
482 | 482 |
template<typename SupplyMap> |
483 | 483 |
CostScaling& supplyMap(const SupplyMap& map) { |
484 | 484 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
485 | 485 |
_supply[_node_id[n]] = map[n]; |
486 | 486 |
} |
487 | 487 |
return *this; |
488 | 488 |
} |
489 | 489 |
|
490 | 490 |
/// \brief Set single source and target nodes and a supply value. |
491 | 491 |
/// |
492 | 492 |
/// This function sets a single source node and a single target node |
493 | 493 |
/// and the required flow value. |
494 | 494 |
/// If neither this function nor \ref supplyMap() is used before |
495 | 495 |
/// calling \ref run(), the supply of each node will be set to zero. |
496 | 496 |
/// |
497 | 497 |
/// Using this function has the same effect as using \ref supplyMap() |
498 | 498 |
/// with such a map in which \c k is assigned to \c s, \c -k is |
499 | 499 |
/// assigned to \c t and all other nodes have zero supply value. |
500 | 500 |
/// |
501 | 501 |
/// \param s The source node. |
502 | 502 |
/// \param t The target node. |
503 | 503 |
/// \param k The required amount of flow from node \c s to node \c t |
504 | 504 |
/// (i.e. the supply of \c s and the demand of \c t). |
505 | 505 |
/// |
506 | 506 |
/// \return <tt>(*this)</tt> |
507 | 507 |
CostScaling& stSupply(const Node& s, const Node& t, Value k) { |
508 | 508 |
for (int i = 0; i != _res_node_num; ++i) { |
509 | 509 |
_supply[i] = 0; |
510 | 510 |
} |
511 | 511 |
_supply[_node_id[s]] = k; |
512 | 512 |
_supply[_node_id[t]] = -k; |
513 | 513 |
return *this; |
514 | 514 |
} |
515 | 515 |
|
516 | 516 |
/// @} |
517 | 517 |
|
518 | 518 |
/// \name Execution control |
519 | 519 |
/// The algorithm can be executed using \ref run(). |
520 | 520 |
|
521 | 521 |
/// @{ |
522 | 522 |
|
523 | 523 |
/// \brief Run the algorithm. |
524 | 524 |
/// |
525 | 525 |
/// This function runs the algorithm. |
526 | 526 |
/// The paramters can be specified using functions \ref lowerMap(), |
527 | 527 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). |
528 | 528 |
/// For example, |
529 | 529 |
/// \code |
530 | 530 |
/// CostScaling<ListDigraph> cs(graph); |
531 | 531 |
/// cs.lowerMap(lower).upperMap(upper).costMap(cost) |
532 | 532 |
/// .supplyMap(sup).run(); |
533 | 533 |
/// \endcode |
534 | 534 |
/// |
535 | 535 |
/// This function can be called more than once. All the parameters |
536 | 536 |
/// that have been given are kept for the next call, unless |
537 | 537 |
/// \ref reset() is called, thus only the modified parameters |
538 | 538 |
/// have to be set again. See \ref reset() for examples. |
539 | 539 |
/// However, the underlying digraph must not be modified after this |
540 | 540 |
/// class have been constructed, since it copies and extends the graph. |
541 | 541 |
/// |
542 | 542 |
/// \param method The internal method that will be used in the |
543 | 543 |
/// algorithm. For more information, see \ref Method. |
544 | 544 |
/// \param factor The cost scaling factor. It must be larger than one. |
545 | 545 |
/// |
546 | 546 |
/// \return \c INFEASIBLE if no feasible flow exists, |
547 | 547 |
/// \n \c OPTIMAL if the problem has optimal solution |
548 | 548 |
/// (i.e. it is feasible and bounded), and the algorithm has found |
549 | 549 |
/// optimal flow and node potentials (primal and dual solutions), |
550 | 550 |
/// \n \c UNBOUNDED if the digraph contains an arc of negative cost |
551 | 551 |
/// and infinite upper bound. It means that the objective function |
552 |
/// is unbounded on that arc, however note that it could actually be |
|
552 |
/// is unbounded on that arc, however, note that it could actually be |
|
553 | 553 |
/// bounded over the feasible flows, but this algroithm cannot handle |
554 | 554 |
/// these cases. |
555 | 555 |
/// |
556 | 556 |
/// \see ProblemType, Method |
557 | 557 |
ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) { |
558 | 558 |
_alpha = factor; |
559 | 559 |
ProblemType pt = init(); |
560 | 560 |
if (pt != OPTIMAL) return pt; |
561 | 561 |
start(method); |
562 | 562 |
return OPTIMAL; |
563 | 563 |
} |
564 | 564 |
|
565 | 565 |
/// \brief Reset all the parameters that have been given before. |
566 | 566 |
/// |
567 | 567 |
/// This function resets all the paramaters that have been given |
568 | 568 |
/// before using functions \ref lowerMap(), \ref upperMap(), |
569 | 569 |
/// \ref costMap(), \ref supplyMap(), \ref stSupply(). |
570 | 570 |
/// |
571 | 571 |
/// It is useful for multiple run() calls. If this function is not |
572 | 572 |
/// used, all the parameters given before are kept for the next |
573 | 573 |
/// \ref run() call. |
574 |
/// However the underlying digraph must not be modified after this |
|
574 |
/// However, the underlying digraph must not be modified after this |
|
575 | 575 |
/// class have been constructed, since it copies and extends the graph. |
576 | 576 |
/// |
577 | 577 |
/// For example, |
578 | 578 |
/// \code |
579 | 579 |
/// CostScaling<ListDigraph> cs(graph); |
580 | 580 |
/// |
581 | 581 |
/// // First run |
582 | 582 |
/// cs.lowerMap(lower).upperMap(upper).costMap(cost) |
583 | 583 |
/// .supplyMap(sup).run(); |
584 | 584 |
/// |
585 | 585 |
/// // Run again with modified cost map (reset() is not called, |
586 | 586 |
/// // so only the cost map have to be set again) |
587 | 587 |
/// cost[e] += 100; |
588 | 588 |
/// cs.costMap(cost).run(); |
589 | 589 |
/// |
590 | 590 |
/// // Run again from scratch using reset() |
591 | 591 |
/// // (the lower bounds will be set to zero on all arcs) |
592 | 592 |
/// cs.reset(); |
593 | 593 |
/// cs.upperMap(capacity).costMap(cost) |
594 | 594 |
/// .supplyMap(sup).run(); |
595 | 595 |
/// \endcode |
596 | 596 |
/// |
597 | 597 |
/// \return <tt>(*this)</tt> |
598 | 598 |
CostScaling& reset() { |
599 | 599 |
for (int i = 0; i != _res_node_num; ++i) { |
600 | 600 |
_supply[i] = 0; |
601 | 601 |
} |
602 | 602 |
int limit = _first_out[_root]; |
603 | 603 |
for (int j = 0; j != limit; ++j) { |
604 | 604 |
_lower[j] = 0; |
605 | 605 |
_upper[j] = INF; |
606 | 606 |
_scost[j] = _forward[j] ? 1 : -1; |
607 | 607 |
} |
608 | 608 |
for (int j = limit; j != _res_arc_num; ++j) { |
609 | 609 |
_lower[j] = 0; |
610 | 610 |
_upper[j] = INF; |
611 | 611 |
_scost[j] = 0; |
612 | 612 |
_scost[_reverse[j]] = 0; |
613 | 613 |
} |
614 | 614 |
_have_lower = false; |
615 | 615 |
return *this; |
616 | 616 |
} |
617 | 617 |
|
618 | 618 |
/// @} |
619 | 619 |
|
620 | 620 |
/// \name Query Functions |
621 | 621 |
/// The results of the algorithm can be obtained using these |
622 | 622 |
/// functions.\n |
623 | 623 |
/// The \ref run() function must be called before using them. |
624 | 624 |
|
625 | 625 |
/// @{ |
626 | 626 |
|
627 | 627 |
/// \brief Return the total cost of the found flow. |
628 | 628 |
/// |
629 | 629 |
/// This function returns the total cost of the found flow. |
630 | 630 |
/// Its complexity is O(e). |
631 | 631 |
/// |
632 | 632 |
/// \note The return type of the function can be specified as a |
633 | 633 |
/// template parameter. For example, |
634 | 634 |
/// \code |
635 | 635 |
/// cs.totalCost<double>(); |
636 | 636 |
/// \endcode |
637 | 637 |
/// It is useful if the total cost cannot be stored in the \c Cost |
638 | 638 |
/// type of the algorithm, which is the default return type of the |
639 | 639 |
/// function. |
640 | 640 |
/// |
641 | 641 |
/// \pre \ref run() must be called before using this function. |
642 | 642 |
template <typename Number> |
643 | 643 |
Number totalCost() const { |
644 | 644 |
Number c = 0; |
645 | 645 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
646 | 646 |
int i = _arc_idb[a]; |
647 | 647 |
c += static_cast<Number>(_res_cap[i]) * |
648 | 648 |
(-static_cast<Number>(_scost[i])); |
649 | 649 |
} |
650 | 650 |
return c; |
651 | 651 |
} |
652 | 652 |
|
653 | 653 |
#ifndef DOXYGEN |
654 | 654 |
Cost totalCost() const { |
655 | 655 |
return totalCost<Cost>(); |
656 | 656 |
} |
657 | 657 |
#endif |
658 | 658 |
|
659 | 659 |
/// \brief Return the flow on the given arc. |
660 | 660 |
/// |
661 | 661 |
/// This function returns the flow on the given arc. |
662 | 662 |
/// |
663 | 663 |
/// \pre \ref run() must be called before using this function. |
664 | 664 |
Value flow(const Arc& a) const { |
665 | 665 |
return _res_cap[_arc_idb[a]]; |
666 | 666 |
} |
667 | 667 |
|
668 | 668 |
/// \brief Return the flow map (the primal solution). |
669 | 669 |
/// |
670 | 670 |
/// This function copies the flow value on each arc into the given |
671 | 671 |
/// map. The \c Value type of the algorithm must be convertible to |
672 | 672 |
/// the \c Value type of the map. |
673 | 673 |
/// |
674 | 674 |
/// \pre \ref run() must be called before using this function. |
675 | 675 |
template <typename FlowMap> |
676 | 676 |
void flowMap(FlowMap &map) const { |
677 | 677 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
678 | 678 |
map.set(a, _res_cap[_arc_idb[a]]); |
679 | 679 |
} |
680 | 680 |
} |
681 | 681 |
|
682 | 682 |
/// \brief Return the potential (dual value) of the given node. |
683 | 683 |
/// |
684 | 684 |
/// This function returns the potential (dual value) of the |
685 | 685 |
/// given node. |
686 | 686 |
/// |
687 | 687 |
/// \pre \ref run() must be called before using this function. |
688 | 688 |
Cost potential(const Node& n) const { |
689 | 689 |
return static_cast<Cost>(_pi[_node_id[n]]); |
690 | 690 |
} |
691 | 691 |
|
692 | 692 |
/// \brief Return the potential map (the dual solution). |
693 | 693 |
/// |
694 | 694 |
/// This function copies the potential (dual value) of each node |
695 | 695 |
/// into the given map. |
696 | 696 |
/// The \c Cost type of the algorithm must be convertible to the |
697 | 697 |
/// \c Value type of the map. |
698 | 698 |
/// |
699 | 699 |
/// \pre \ref run() must be called before using this function. |
700 | 700 |
template <typename PotentialMap> |
701 | 701 |
void potentialMap(PotentialMap &map) const { |
702 | 702 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
703 | 703 |
map.set(n, static_cast<Cost>(_pi[_node_id[n]])); |
704 | 704 |
} |
705 | 705 |
} |
706 | 706 |
|
707 | 707 |
/// @} |
708 | 708 |
|
709 | 709 |
private: |
710 | 710 |
|
711 | 711 |
// Initialize the algorithm |
712 | 712 |
ProblemType init() { |
713 | 713 |
if (_res_node_num == 0) return INFEASIBLE; |
714 | 714 |
|
715 | 715 |
// Check the sum of supply values |
716 | 716 |
_sum_supply = 0; |
717 | 717 |
for (int i = 0; i != _root; ++i) { |
718 | 718 |
_sum_supply += _supply[i]; |
719 | 719 |
} |
720 | 720 |
if (_sum_supply > 0) return INFEASIBLE; |
721 | 721 |
|
722 | 722 |
|
723 | 723 |
// Initialize vectors |
724 | 724 |
for (int i = 0; i != _res_node_num; ++i) { |
725 | 725 |
_pi[i] = 0; |
726 | 726 |
_excess[i] = _supply[i]; |
727 | 727 |
} |
728 | 728 |
|
729 | 729 |
// Remove infinite upper bounds and check negative arcs |
730 | 730 |
const Value MAX = std::numeric_limits<Value>::max(); |
731 | 731 |
int last_out; |
732 | 732 |
if (_have_lower) { |
733 | 733 |
for (int i = 0; i != _root; ++i) { |
734 | 734 |
last_out = _first_out[i+1]; |
735 | 735 |
for (int j = _first_out[i]; j != last_out; ++j) { |
736 | 736 |
if (_forward[j]) { |
737 | 737 |
Value c = _scost[j] < 0 ? _upper[j] : _lower[j]; |
738 | 738 |
if (c >= MAX) return UNBOUNDED; |
739 | 739 |
_excess[i] -= c; |
740 | 740 |
_excess[_target[j]] += c; |
741 | 741 |
} |
742 | 742 |
} |
743 | 743 |
} |
744 | 744 |
} else { |
745 | 745 |
for (int i = 0; i != _root; ++i) { |
746 | 746 |
last_out = _first_out[i+1]; |
747 | 747 |
for (int j = _first_out[i]; j != last_out; ++j) { |
748 | 748 |
if (_forward[j] && _scost[j] < 0) { |
749 | 749 |
Value c = _upper[j]; |
750 | 750 |
if (c >= MAX) return UNBOUNDED; |
751 | 751 |
_excess[i] -= c; |
752 | 752 |
_excess[_target[j]] += c; |
753 | 753 |
} |
754 | 754 |
} |
755 | 755 |
} |
756 | 756 |
} |
757 | 757 |
Value ex, max_cap = 0; |
758 | 758 |
for (int i = 0; i != _res_node_num; ++i) { |
759 | 759 |
ex = _excess[i]; |
760 | 760 |
_excess[i] = 0; |
761 | 761 |
if (ex < 0) max_cap -= ex; |
762 | 762 |
} |
763 | 763 |
for (int j = 0; j != _res_arc_num; ++j) { |
764 | 764 |
if (_upper[j] >= MAX) _upper[j] = max_cap; |
765 | 765 |
} |
766 | 766 |
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_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 |
/// This algorithm is a specialized version of the linear programming |
|
47 |
/// simplex method directly for the minimum cost flow problem. |
|
48 |
/// |
|
46 |
/// This algorithm is a highly efficient specialized version of the |
|
47 |
/// linear programming simplex method directly for the minimum cost |
|
48 |
/// flow problem. |
|
49 | 49 |
/// |
50 |
/// In general this class is the fastest implementation available |
|
51 |
/// in LEMON for the minimum cost flow problem. |
|
52 |
/// |
|
50 |
/// In general, %NetworkSimplex is the fastest implementation available |
|
51 |
/// in LEMON for this problem. |
|
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 |
/// \tparam V The |
|
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 |
/// \tparam C The |
|
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 number types must be signed and all input data must |
|
67 | 67 |
/// be integer. |
68 | 68 |
/// |
69 | 69 |
/// \note %NetworkSimplex provides five different pivot rule |
70 | 70 |
/// implementations, from which the most efficient one is used |
71 | 71 |
/// by default. For more information, see \ref PivotRule. |
72 | 72 |
template <typename GR, typename V = int, typename C = V> |
73 | 73 |
class NetworkSimplex |
74 | 74 |
{ |
75 | 75 |
public: |
76 | 76 |
|
77 | 77 |
/// The type of the flow amounts, capacity bounds and supply values |
78 | 78 |
typedef V Value; |
79 | 79 |
/// The type of the arc costs |
80 | 80 |
typedef C Cost; |
81 | 81 |
|
82 | 82 |
public: |
83 | 83 |
|
84 | 84 |
/// \brief Problem type constants for the \c run() function. |
85 | 85 |
/// |
86 | 86 |
/// Enum type containing the problem type constants that can be |
87 | 87 |
/// returned by the \ref run() function of the algorithm. |
88 | 88 |
enum ProblemType { |
89 | 89 |
/// The problem has no feasible solution (flow). |
90 | 90 |
INFEASIBLE, |
91 | 91 |
/// The problem has optimal solution (i.e. it is feasible and |
92 | 92 |
/// bounded), and the algorithm has found optimal flow and node |
93 | 93 |
/// potentials (primal and dual solutions). |
94 | 94 |
OPTIMAL, |
95 | 95 |
/// The objective function of the problem is unbounded, i.e. |
96 | 96 |
/// there is a directed cycle having negative total cost and |
97 | 97 |
/// infinite upper bound. |
98 | 98 |
UNBOUNDED |
99 | 99 |
}; |
100 | 100 |
|
101 | 101 |
/// \brief Constants for selecting the type of the supply constraints. |
102 | 102 |
/// |
103 | 103 |
/// Enum type containing constants for selecting the supply type, |
104 | 104 |
/// i.e. the direction of the inequalities in the supply/demand |
105 | 105 |
/// constraints of the \ref min_cost_flow "minimum cost flow problem". |
106 | 106 |
/// |
107 | 107 |
/// The default supply type is \c GEQ, the \c LEQ type can be |
108 | 108 |
/// selected using \ref supplyType(). |
109 | 109 |
/// The equality form is a special case of both supply types. |
110 | 110 |
enum SupplyType { |
111 | 111 |
/// This option means that there are <em>"greater or equal"</em> |
112 | 112 |
/// supply/demand constraints in the definition of the problem. |
113 | 113 |
GEQ, |
114 | 114 |
/// This option means that there are <em>"less or equal"</em> |
115 | 115 |
/// supply/demand constraints in the definition of the problem. |
116 | 116 |
LEQ |
117 | 117 |
}; |
118 | 118 |
|
119 | 119 |
/// \brief Constants for selecting the pivot rule. |
120 | 120 |
/// |
121 | 121 |
/// Enum type containing constants for selecting the pivot rule for |
122 | 122 |
/// the \ref run() function. |
123 | 123 |
/// |
124 | 124 |
/// \ref NetworkSimplex provides five different pivot rule |
125 | 125 |
/// implementations that significantly affect the running time |
126 | 126 |
/// of the algorithm. |
127 | 127 |
/// By default, \ref BLOCK_SEARCH "Block Search" is used, which |
128 | 128 |
/// proved to be the most efficient and the most robust on various |
129 |
/// test inputs |
|
129 |
/// test inputs. |
|
130 | 130 |
/// However, another pivot rule can be selected using the \ref run() |
131 | 131 |
/// function with the proper parameter. |
132 | 132 |
enum PivotRule { |
133 | 133 |
|
134 | 134 |
/// The \e First \e Eligible pivot rule. |
135 | 135 |
/// The next eligible arc is selected in a wraparound fashion |
136 | 136 |
/// in every iteration. |
137 | 137 |
FIRST_ELIGIBLE, |
138 | 138 |
|
139 | 139 |
/// The \e Best \e Eligible pivot rule. |
140 | 140 |
/// The best eligible arc is selected in every iteration. |
141 | 141 |
BEST_ELIGIBLE, |
142 | 142 |
|
143 | 143 |
/// The \e Block \e Search pivot rule. |
144 | 144 |
/// A specified number of arcs are examined in every iteration |
145 | 145 |
/// in a wraparound fashion and the best eligible arc is selected |
146 | 146 |
/// from this block. |
147 | 147 |
BLOCK_SEARCH, |
148 | 148 |
|
149 | 149 |
/// The \e Candidate \e List pivot rule. |
150 | 150 |
/// In a major iteration a candidate list is built from eligible arcs |
151 | 151 |
/// in a wraparound fashion and in the following minor iterations |
152 | 152 |
/// the best eligible arc is selected from this list. |
153 | 153 |
CANDIDATE_LIST, |
154 | 154 |
|
155 | 155 |
/// The \e Altering \e Candidate \e List pivot rule. |
156 | 156 |
/// It is a modified version of the Candidate List method. |
157 | 157 |
/// It keeps only the several best eligible arcs from the former |
158 | 158 |
/// candidate list and extends this list in every iteration. |
159 | 159 |
ALTERING_LIST |
160 | 160 |
}; |
161 | 161 |
|
162 | 162 |
private: |
163 | 163 |
|
164 | 164 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
165 | 165 |
|
166 | 166 |
typedef std::vector<int> IntVector; |
167 | 167 |
typedef std::vector<char> CharVector; |
168 | 168 |
typedef std::vector<Value> ValueVector; |
169 | 169 |
typedef std::vector<Cost> CostVector; |
170 | 170 |
|
171 | 171 |
// State constants for arcs |
172 | 172 |
enum ArcStateEnum { |
173 | 173 |
STATE_UPPER = -1, |
174 | 174 |
STATE_TREE = 0, |
175 | 175 |
STATE_LOWER = 1 |
176 | 176 |
}; |
177 | 177 |
|
178 | 178 |
private: |
179 | 179 |
|
180 | 180 |
// Data related to the underlying digraph |
181 | 181 |
const GR &_graph; |
182 | 182 |
int _node_num; |
183 | 183 |
int _arc_num; |
184 | 184 |
int _all_arc_num; |
185 | 185 |
int _search_arc_num; |
186 | 186 |
|
187 | 187 |
// Parameters of the problem |
188 | 188 |
bool _have_lower; |
189 | 189 |
SupplyType _stype; |
190 | 190 |
Value _sum_supply; |
191 | 191 |
|
192 | 192 |
// Data structures for storing the digraph |
193 | 193 |
IntNodeMap _node_id; |
194 | 194 |
IntArcMap _arc_id; |
195 | 195 |
IntVector _source; |
196 | 196 |
IntVector _target; |
197 | 197 |
|
198 | 198 |
// Node and arc data |
199 | 199 |
ValueVector _lower; |
200 | 200 |
ValueVector _upper; |
201 | 201 |
ValueVector _cap; |
202 | 202 |
CostVector _cost; |
203 | 203 |
ValueVector _supply; |
204 | 204 |
ValueVector _flow; |
205 | 205 |
CostVector _pi; |
206 | 206 |
|
207 | 207 |
// Data for storing the spanning tree structure |
208 | 208 |
IntVector _parent; |
209 | 209 |
IntVector _pred; |
210 | 210 |
IntVector _thread; |
211 | 211 |
IntVector _rev_thread; |
212 | 212 |
IntVector _succ_num; |
213 | 213 |
IntVector _last_succ; |
214 | 214 |
IntVector _dirty_revs; |
215 | 215 |
CharVector _forward; |
216 | 216 |
CharVector _state; |
217 | 217 |
int _root; |
218 | 218 |
|
219 | 219 |
// Temporary data used in the current pivot iteration |
220 | 220 |
int in_arc, join, u_in, v_in, u_out, v_out; |
221 | 221 |
int first, second, right, last; |
222 | 222 |
int stem, par_stem, new_stem; |
223 | 223 |
Value delta; |
224 | 224 |
|
225 | 225 |
const Value MAX; |
226 | 226 |
|
227 | 227 |
public: |
228 | 228 |
|
229 | 229 |
/// \brief Constant for infinite upper bounds (capacities). |
230 | 230 |
/// |
231 | 231 |
/// Constant for infinite upper bounds (capacities). |
232 | 232 |
/// It is \c std::numeric_limits<Value>::infinity() if available, |
233 | 233 |
/// \c std::numeric_limits<Value>::max() otherwise. |
234 | 234 |
const Value INF; |
235 | 235 |
|
236 | 236 |
private: |
237 | 237 |
|
238 | 238 |
// Implementation of the First Eligible pivot rule |
239 | 239 |
class FirstEligiblePivotRule |
240 | 240 |
{ |
241 | 241 |
private: |
242 | 242 |
|
243 | 243 |
// References to the NetworkSimplex class |
244 | 244 |
const IntVector &_source; |
245 | 245 |
const IntVector &_target; |
246 | 246 |
const CostVector &_cost; |
247 | 247 |
const CharVector &_state; |
248 | 248 |
const CostVector &_pi; |
249 | 249 |
int &_in_arc; |
250 | 250 |
int _search_arc_num; |
251 | 251 |
|
252 | 252 |
// Pivot rule data |
253 | 253 |
int _next_arc; |
254 | 254 |
|
255 | 255 |
public: |
256 | 256 |
|
257 | 257 |
// Constructor |
258 | 258 |
FirstEligiblePivotRule(NetworkSimplex &ns) : |
259 | 259 |
_source(ns._source), _target(ns._target), |
260 | 260 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
261 | 261 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
262 | 262 |
_next_arc(0) |
263 | 263 |
{} |
264 | 264 |
|
265 | 265 |
// Find next entering arc |
266 | 266 |
bool findEnteringArc() { |
267 | 267 |
Cost c; |
268 | 268 |
for (int e = _next_arc; e < _search_arc_num; ++e) { |
269 | 269 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
270 | 270 |
if (c < 0) { |
271 | 271 |
_in_arc = e; |
272 | 272 |
_next_arc = e + 1; |
273 | 273 |
return true; |
274 | 274 |
} |
275 | 275 |
} |
276 | 276 |
for (int e = 0; e < _next_arc; ++e) { |
277 | 277 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
278 | 278 |
if (c < 0) { |
279 | 279 |
_in_arc = e; |
280 | 280 |
_next_arc = e + 1; |
281 | 281 |
return true; |
282 | 282 |
} |
283 | 283 |
} |
284 | 284 |
return false; |
285 | 285 |
} |
286 | 286 |
|
287 | 287 |
}; //class FirstEligiblePivotRule |
288 | 288 |
|
289 | 289 |
|
290 | 290 |
// Implementation of the Best Eligible pivot rule |
291 | 291 |
class BestEligiblePivotRule |
292 | 292 |
{ |
293 | 293 |
private: |
294 | 294 |
|
295 | 295 |
// References to the NetworkSimplex class |
296 | 296 |
const IntVector &_source; |
297 | 297 |
const IntVector &_target; |
298 | 298 |
const CostVector &_cost; |
299 | 299 |
const CharVector &_state; |
300 | 300 |
const CostVector &_pi; |
301 | 301 |
int &_in_arc; |
302 | 302 |
int _search_arc_num; |
303 | 303 |
|
304 | 304 |
public: |
305 | 305 |
|
306 | 306 |
// Constructor |
307 | 307 |
BestEligiblePivotRule(NetworkSimplex &ns) : |
308 | 308 |
_source(ns._source), _target(ns._target), |
309 | 309 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
310 | 310 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num) |
311 | 311 |
{} |
312 | 312 |
|
313 | 313 |
// Find next entering arc |
314 | 314 |
bool findEnteringArc() { |
315 | 315 |
Cost c, min = 0; |
316 | 316 |
for (int e = 0; e < _search_arc_num; ++e) { |
317 | 317 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
318 | 318 |
if (c < min) { |
319 | 319 |
min = c; |
320 | 320 |
_in_arc = e; |
321 | 321 |
} |
... | ... |
@@ -448,477 +448,477 @@ |
448 | 448 |
Cost min, c; |
449 | 449 |
int e; |
450 | 450 |
if (_curr_length > 0 && _minor_count < _minor_limit) { |
451 | 451 |
// Minor iteration: select the best eligible arc from the |
452 | 452 |
// current candidate list |
453 | 453 |
++_minor_count; |
454 | 454 |
min = 0; |
455 | 455 |
for (int i = 0; i < _curr_length; ++i) { |
456 | 456 |
e = _candidates[i]; |
457 | 457 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
458 | 458 |
if (c < min) { |
459 | 459 |
min = c; |
460 | 460 |
_in_arc = e; |
461 | 461 |
} |
462 | 462 |
else if (c >= 0) { |
463 | 463 |
_candidates[i--] = _candidates[--_curr_length]; |
464 | 464 |
} |
465 | 465 |
} |
466 | 466 |
if (min < 0) return true; |
467 | 467 |
} |
468 | 468 |
|
469 | 469 |
// Major iteration: build a new candidate list |
470 | 470 |
min = 0; |
471 | 471 |
_curr_length = 0; |
472 | 472 |
for (e = _next_arc; e < _search_arc_num; ++e) { |
473 | 473 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
474 | 474 |
if (c < 0) { |
475 | 475 |
_candidates[_curr_length++] = e; |
476 | 476 |
if (c < min) { |
477 | 477 |
min = c; |
478 | 478 |
_in_arc = e; |
479 | 479 |
} |
480 | 480 |
if (_curr_length == _list_length) goto search_end; |
481 | 481 |
} |
482 | 482 |
} |
483 | 483 |
for (e = 0; e < _next_arc; ++e) { |
484 | 484 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
485 | 485 |
if (c < 0) { |
486 | 486 |
_candidates[_curr_length++] = e; |
487 | 487 |
if (c < min) { |
488 | 488 |
min = c; |
489 | 489 |
_in_arc = e; |
490 | 490 |
} |
491 | 491 |
if (_curr_length == _list_length) goto search_end; |
492 | 492 |
} |
493 | 493 |
} |
494 | 494 |
if (_curr_length == 0) return false; |
495 | 495 |
|
496 | 496 |
search_end: |
497 | 497 |
_minor_count = 1; |
498 | 498 |
_next_arc = e; |
499 | 499 |
return true; |
500 | 500 |
} |
501 | 501 |
|
502 | 502 |
}; //class CandidateListPivotRule |
503 | 503 |
|
504 | 504 |
|
505 | 505 |
// Implementation of the Altering Candidate List pivot rule |
506 | 506 |
class AlteringListPivotRule |
507 | 507 |
{ |
508 | 508 |
private: |
509 | 509 |
|
510 | 510 |
// References to the NetworkSimplex class |
511 | 511 |
const IntVector &_source; |
512 | 512 |
const IntVector &_target; |
513 | 513 |
const CostVector &_cost; |
514 | 514 |
const CharVector &_state; |
515 | 515 |
const CostVector &_pi; |
516 | 516 |
int &_in_arc; |
517 | 517 |
int _search_arc_num; |
518 | 518 |
|
519 | 519 |
// Pivot rule data |
520 | 520 |
int _block_size, _head_length, _curr_length; |
521 | 521 |
int _next_arc; |
522 | 522 |
IntVector _candidates; |
523 | 523 |
CostVector _cand_cost; |
524 | 524 |
|
525 | 525 |
// Functor class to compare arcs during sort of the candidate list |
526 | 526 |
class SortFunc |
527 | 527 |
{ |
528 | 528 |
private: |
529 | 529 |
const CostVector &_map; |
530 | 530 |
public: |
531 | 531 |
SortFunc(const CostVector &map) : _map(map) {} |
532 | 532 |
bool operator()(int left, int right) { |
533 | 533 |
return _map[left] > _map[right]; |
534 | 534 |
} |
535 | 535 |
}; |
536 | 536 |
|
537 | 537 |
SortFunc _sort_func; |
538 | 538 |
|
539 | 539 |
public: |
540 | 540 |
|
541 | 541 |
// Constructor |
542 | 542 |
AlteringListPivotRule(NetworkSimplex &ns) : |
543 | 543 |
_source(ns._source), _target(ns._target), |
544 | 544 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
545 | 545 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
546 | 546 |
_next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost) |
547 | 547 |
{ |
548 | 548 |
// The main parameters of the pivot rule |
549 | 549 |
const double BLOCK_SIZE_FACTOR = 1.0; |
550 | 550 |
const int MIN_BLOCK_SIZE = 10; |
551 | 551 |
const double HEAD_LENGTH_FACTOR = 0.1; |
552 | 552 |
const int MIN_HEAD_LENGTH = 3; |
553 | 553 |
|
554 | 554 |
_block_size = std::max( int(BLOCK_SIZE_FACTOR * |
555 | 555 |
std::sqrt(double(_search_arc_num))), |
556 | 556 |
MIN_BLOCK_SIZE ); |
557 | 557 |
_head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size), |
558 | 558 |
MIN_HEAD_LENGTH ); |
559 | 559 |
_candidates.resize(_head_length + _block_size); |
560 | 560 |
_curr_length = 0; |
561 | 561 |
} |
562 | 562 |
|
563 | 563 |
// Find next entering arc |
564 | 564 |
bool findEnteringArc() { |
565 | 565 |
// Check the current candidate list |
566 | 566 |
int e; |
567 | 567 |
for (int i = 0; i < _curr_length; ++i) { |
568 | 568 |
e = _candidates[i]; |
569 | 569 |
_cand_cost[e] = _state[e] * |
570 | 570 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
571 | 571 |
if (_cand_cost[e] >= 0) { |
572 | 572 |
_candidates[i--] = _candidates[--_curr_length]; |
573 | 573 |
} |
574 | 574 |
} |
575 | 575 |
|
576 | 576 |
// Extend the list |
577 | 577 |
int cnt = _block_size; |
578 | 578 |
int limit = _head_length; |
579 | 579 |
|
580 | 580 |
for (e = _next_arc; e < _search_arc_num; ++e) { |
581 | 581 |
_cand_cost[e] = _state[e] * |
582 | 582 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
583 | 583 |
if (_cand_cost[e] < 0) { |
584 | 584 |
_candidates[_curr_length++] = e; |
585 | 585 |
} |
586 | 586 |
if (--cnt == 0) { |
587 | 587 |
if (_curr_length > limit) goto search_end; |
588 | 588 |
limit = 0; |
589 | 589 |
cnt = _block_size; |
590 | 590 |
} |
591 | 591 |
} |
592 | 592 |
for (e = 0; e < _next_arc; ++e) { |
593 | 593 |
_cand_cost[e] = _state[e] * |
594 | 594 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
595 | 595 |
if (_cand_cost[e] < 0) { |
596 | 596 |
_candidates[_curr_length++] = e; |
597 | 597 |
} |
598 | 598 |
if (--cnt == 0) { |
599 | 599 |
if (_curr_length > limit) goto search_end; |
600 | 600 |
limit = 0; |
601 | 601 |
cnt = _block_size; |
602 | 602 |
} |
603 | 603 |
} |
604 | 604 |
if (_curr_length == 0) return false; |
605 | 605 |
|
606 | 606 |
search_end: |
607 | 607 |
|
608 | 608 |
// Make heap of the candidate list (approximating a partial sort) |
609 | 609 |
make_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
610 | 610 |
_sort_func ); |
611 | 611 |
|
612 | 612 |
// Pop the first element of the heap |
613 | 613 |
_in_arc = _candidates[0]; |
614 | 614 |
_next_arc = e; |
615 | 615 |
pop_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
616 | 616 |
_sort_func ); |
617 | 617 |
_curr_length = std::min(_head_length, _curr_length - 1); |
618 | 618 |
return true; |
619 | 619 |
} |
620 | 620 |
|
621 | 621 |
}; //class AlteringListPivotRule |
622 | 622 |
|
623 | 623 |
public: |
624 | 624 |
|
625 | 625 |
/// \brief Constructor. |
626 | 626 |
/// |
627 | 627 |
/// The constructor of the class. |
628 | 628 |
/// |
629 | 629 |
/// \param graph The digraph the algorithm runs on. |
630 | 630 |
/// \param arc_mixing Indicate if the arcs have to be stored in a |
631 | 631 |
/// mixed order in the internal data structure. |
632 | 632 |
/// In special cases, it could lead to better overall performance, |
633 | 633 |
/// but it is usually slower. Therefore it is disabled by default. |
634 | 634 |
NetworkSimplex(const GR& graph, bool arc_mixing = false) : |
635 | 635 |
_graph(graph), _node_id(graph), _arc_id(graph), |
636 | 636 |
MAX(std::numeric_limits<Value>::max()), |
637 | 637 |
INF(std::numeric_limits<Value>::has_infinity ? |
638 | 638 |
std::numeric_limits<Value>::infinity() : MAX) |
639 | 639 |
{ |
640 |
// Check the |
|
640 |
// Check the number types |
|
641 | 641 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
642 | 642 |
"The flow type of NetworkSimplex must be signed"); |
643 | 643 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
644 | 644 |
"The cost type of NetworkSimplex must be signed"); |
645 | 645 |
|
646 | 646 |
// Resize vectors |
647 | 647 |
_node_num = countNodes(_graph); |
648 | 648 |
_arc_num = countArcs(_graph); |
649 | 649 |
int all_node_num = _node_num + 1; |
650 | 650 |
int max_arc_num = _arc_num + 2 * _node_num; |
651 | 651 |
|
652 | 652 |
_source.resize(max_arc_num); |
653 | 653 |
_target.resize(max_arc_num); |
654 | 654 |
|
655 | 655 |
_lower.resize(_arc_num); |
656 | 656 |
_upper.resize(_arc_num); |
657 | 657 |
_cap.resize(max_arc_num); |
658 | 658 |
_cost.resize(max_arc_num); |
659 | 659 |
_supply.resize(all_node_num); |
660 | 660 |
_flow.resize(max_arc_num); |
661 | 661 |
_pi.resize(all_node_num); |
662 | 662 |
|
663 | 663 |
_parent.resize(all_node_num); |
664 | 664 |
_pred.resize(all_node_num); |
665 | 665 |
_forward.resize(all_node_num); |
666 | 666 |
_thread.resize(all_node_num); |
667 | 667 |
_rev_thread.resize(all_node_num); |
668 | 668 |
_succ_num.resize(all_node_num); |
669 | 669 |
_last_succ.resize(all_node_num); |
670 | 670 |
_state.resize(max_arc_num); |
671 | 671 |
|
672 | 672 |
// Copy the graph |
673 | 673 |
int i = 0; |
674 | 674 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
675 | 675 |
_node_id[n] = i; |
676 | 676 |
} |
677 | 677 |
if (arc_mixing) { |
678 | 678 |
// Store the arcs in a mixed order |
679 | 679 |
int k = std::max(int(std::sqrt(double(_arc_num))), 10); |
680 | 680 |
int i = 0, j = 0; |
681 | 681 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
682 | 682 |
_arc_id[a] = i; |
683 | 683 |
_source[i] = _node_id[_graph.source(a)]; |
684 | 684 |
_target[i] = _node_id[_graph.target(a)]; |
685 | 685 |
if ((i += k) >= _arc_num) i = ++j; |
686 | 686 |
} |
687 | 687 |
} else { |
688 | 688 |
// Store the arcs in the original order |
689 | 689 |
int i = 0; |
690 | 690 |
for (ArcIt a(_graph); a != INVALID; ++a, ++i) { |
691 | 691 |
_arc_id[a] = i; |
692 | 692 |
_source[i] = _node_id[_graph.source(a)]; |
693 | 693 |
_target[i] = _node_id[_graph.target(a)]; |
694 | 694 |
} |
695 | 695 |
} |
696 | 696 |
|
697 | 697 |
// Reset parameters |
698 | 698 |
reset(); |
699 | 699 |
} |
700 | 700 |
|
701 | 701 |
/// \name Parameters |
702 | 702 |
/// The parameters of the algorithm can be specified using these |
703 | 703 |
/// functions. |
704 | 704 |
|
705 | 705 |
/// @{ |
706 | 706 |
|
707 | 707 |
/// \brief Set the lower bounds on the arcs. |
708 | 708 |
/// |
709 | 709 |
/// This function sets the lower bounds on the arcs. |
710 | 710 |
/// If it is not used before calling \ref run(), the lower bounds |
711 | 711 |
/// will be set to zero on all arcs. |
712 | 712 |
/// |
713 | 713 |
/// \param map An arc map storing the lower bounds. |
714 | 714 |
/// Its \c Value type must be convertible to the \c Value type |
715 | 715 |
/// of the algorithm. |
716 | 716 |
/// |
717 | 717 |
/// \return <tt>(*this)</tt> |
718 | 718 |
template <typename LowerMap> |
719 | 719 |
NetworkSimplex& lowerMap(const LowerMap& map) { |
720 | 720 |
_have_lower = true; |
721 | 721 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
722 | 722 |
_lower[_arc_id[a]] = map[a]; |
723 | 723 |
} |
724 | 724 |
return *this; |
725 | 725 |
} |
726 | 726 |
|
727 | 727 |
/// \brief Set the upper bounds (capacities) on the arcs. |
728 | 728 |
/// |
729 | 729 |
/// This function sets the upper bounds (capacities) on the arcs. |
730 | 730 |
/// If it is not used before calling \ref run(), the upper bounds |
731 | 731 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be |
732 |
/// unbounded from above |
|
732 |
/// unbounded from above). |
|
733 | 733 |
/// |
734 | 734 |
/// \param map An arc map storing the upper bounds. |
735 | 735 |
/// Its \c Value type must be convertible to the \c Value type |
736 | 736 |
/// of the algorithm. |
737 | 737 |
/// |
738 | 738 |
/// \return <tt>(*this)</tt> |
739 | 739 |
template<typename UpperMap> |
740 | 740 |
NetworkSimplex& upperMap(const UpperMap& map) { |
741 | 741 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
742 | 742 |
_upper[_arc_id[a]] = map[a]; |
743 | 743 |
} |
744 | 744 |
return *this; |
745 | 745 |
} |
746 | 746 |
|
747 | 747 |
/// \brief Set the costs of the arcs. |
748 | 748 |
/// |
749 | 749 |
/// This function sets the costs of the arcs. |
750 | 750 |
/// If it is not used before calling \ref run(), the costs |
751 | 751 |
/// will be set to \c 1 on all arcs. |
752 | 752 |
/// |
753 | 753 |
/// \param map An arc map storing the costs. |
754 | 754 |
/// Its \c Value type must be convertible to the \c Cost type |
755 | 755 |
/// of the algorithm. |
756 | 756 |
/// |
757 | 757 |
/// \return <tt>(*this)</tt> |
758 | 758 |
template<typename CostMap> |
759 | 759 |
NetworkSimplex& costMap(const CostMap& map) { |
760 | 760 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
761 | 761 |
_cost[_arc_id[a]] = map[a]; |
762 | 762 |
} |
763 | 763 |
return *this; |
764 | 764 |
} |
765 | 765 |
|
766 | 766 |
/// \brief Set the supply values of the nodes. |
767 | 767 |
/// |
768 | 768 |
/// This function sets the supply values of the nodes. |
769 | 769 |
/// If neither this function nor \ref stSupply() is used before |
770 | 770 |
/// calling \ref run(), the supply of each node will be set to zero. |
771 | 771 |
/// |
772 | 772 |
/// \param map A node map storing the supply values. |
773 | 773 |
/// Its \c Value type must be convertible to the \c Value type |
774 | 774 |
/// of the algorithm. |
775 | 775 |
/// |
776 | 776 |
/// \return <tt>(*this)</tt> |
777 | 777 |
template<typename SupplyMap> |
778 | 778 |
NetworkSimplex& supplyMap(const SupplyMap& map) { |
779 | 779 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
780 | 780 |
_supply[_node_id[n]] = map[n]; |
781 | 781 |
} |
782 | 782 |
return *this; |
783 | 783 |
} |
784 | 784 |
|
785 | 785 |
/// \brief Set single source and target nodes and a supply value. |
786 | 786 |
/// |
787 | 787 |
/// This function sets a single source node and a single target node |
788 | 788 |
/// and the required flow value. |
789 | 789 |
/// If neither this function nor \ref supplyMap() is used before |
790 | 790 |
/// calling \ref run(), the supply of each node will be set to zero. |
791 | 791 |
/// |
792 | 792 |
/// Using this function has the same effect as using \ref supplyMap() |
793 | 793 |
/// with such a map in which \c k is assigned to \c s, \c -k is |
794 | 794 |
/// assigned to \c t and all other nodes have zero supply value. |
795 | 795 |
/// |
796 | 796 |
/// \param s The source node. |
797 | 797 |
/// \param t The target node. |
798 | 798 |
/// \param k The required amount of flow from node \c s to node \c t |
799 | 799 |
/// (i.e. the supply of \c s and the demand of \c t). |
800 | 800 |
/// |
801 | 801 |
/// \return <tt>(*this)</tt> |
802 | 802 |
NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) { |
803 | 803 |
for (int i = 0; i != _node_num; ++i) { |
804 | 804 |
_supply[i] = 0; |
805 | 805 |
} |
806 | 806 |
_supply[_node_id[s]] = k; |
807 | 807 |
_supply[_node_id[t]] = -k; |
808 | 808 |
return *this; |
809 | 809 |
} |
810 | 810 |
|
811 | 811 |
/// \brief Set the type of the supply constraints. |
812 | 812 |
/// |
813 | 813 |
/// This function sets the type of the supply/demand constraints. |
814 | 814 |
/// If it is not used before calling \ref run(), the \ref GEQ supply |
815 | 815 |
/// type will be used. |
816 | 816 |
/// |
817 | 817 |
/// For more information, see \ref SupplyType. |
818 | 818 |
/// |
819 | 819 |
/// \return <tt>(*this)</tt> |
820 | 820 |
NetworkSimplex& supplyType(SupplyType supply_type) { |
821 | 821 |
_stype = supply_type; |
822 | 822 |
return *this; |
823 | 823 |
} |
824 | 824 |
|
825 | 825 |
/// @} |
826 | 826 |
|
827 | 827 |
/// \name Execution Control |
828 | 828 |
/// The algorithm can be executed using \ref run(). |
829 | 829 |
|
830 | 830 |
/// @{ |
831 | 831 |
|
832 | 832 |
/// \brief Run the algorithm. |
833 | 833 |
/// |
834 | 834 |
/// This function runs the algorithm. |
835 | 835 |
/// The paramters can be specified using functions \ref lowerMap(), |
836 | 836 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), |
837 | 837 |
/// \ref supplyType(). |
838 | 838 |
/// For example, |
839 | 839 |
/// \code |
840 | 840 |
/// NetworkSimplex<ListDigraph> ns(graph); |
841 | 841 |
/// ns.lowerMap(lower).upperMap(upper).costMap(cost) |
842 | 842 |
/// .supplyMap(sup).run(); |
843 | 843 |
/// \endcode |
844 | 844 |
/// |
845 | 845 |
/// This function can be called more than once. All the parameters |
846 | 846 |
/// that have been given are kept for the next call, unless |
847 | 847 |
/// \ref reset() is called, thus only the modified parameters |
848 | 848 |
/// have to be set again. See \ref reset() for examples. |
849 | 849 |
/// However, the underlying digraph must not be modified after this |
850 | 850 |
/// class have been constructed, since it copies and extends the graph. |
851 | 851 |
/// |
852 | 852 |
/// \param pivot_rule The pivot rule that will be used during the |
853 | 853 |
/// algorithm. For more information, see \ref PivotRule. |
854 | 854 |
/// |
855 | 855 |
/// \return \c INFEASIBLE if no feasible flow exists, |
856 | 856 |
/// \n \c OPTIMAL if the problem has optimal solution |
857 | 857 |
/// (i.e. it is feasible and bounded), and the algorithm has found |
858 | 858 |
/// optimal flow and node potentials (primal and dual solutions), |
859 | 859 |
/// \n \c UNBOUNDED if the objective function of the problem is |
860 | 860 |
/// unbounded, i.e. there is a directed cycle having negative total |
861 | 861 |
/// cost and infinite upper bound. |
862 | 862 |
/// |
863 | 863 |
/// \see ProblemType, PivotRule |
864 | 864 |
ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) { |
865 | 865 |
if (!init()) return INFEASIBLE; |
866 | 866 |
return start(pivot_rule); |
867 | 867 |
} |
868 | 868 |
|
869 | 869 |
/// \brief Reset all the parameters that have been given before. |
870 | 870 |
/// |
871 | 871 |
/// This function resets all the paramaters that have been given |
872 | 872 |
/// before using functions \ref lowerMap(), \ref upperMap(), |
873 | 873 |
/// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(). |
874 | 874 |
/// |
875 | 875 |
/// It is useful for multiple run() calls. If this function is not |
876 | 876 |
/// used, all the parameters given before are kept for the next |
877 | 877 |
/// \ref run() call. |
878 | 878 |
/// However, the underlying digraph must not be modified after this |
879 | 879 |
/// class have been constructed, since it copies and extends the graph. |
880 | 880 |
/// |
881 | 881 |
/// For example, |
882 | 882 |
/// \code |
883 | 883 |
/// NetworkSimplex<ListDigraph> ns(graph); |
884 | 884 |
/// |
885 | 885 |
/// // First run |
886 | 886 |
/// ns.lowerMap(lower).upperMap(upper).costMap(cost) |
887 | 887 |
/// .supplyMap(sup).run(); |
888 | 888 |
/// |
889 | 889 |
/// // Run again with modified cost map (reset() is not called, |
890 | 890 |
/// // so only the cost map have to be set again) |
891 | 891 |
/// cost[e] += 100; |
892 | 892 |
/// ns.costMap(cost).run(); |
893 | 893 |
/// |
894 | 894 |
/// // Run again from scratch using reset() |
895 | 895 |
/// // (the lower bounds will be set to zero on all arcs) |
896 | 896 |
/// ns.reset(); |
897 | 897 |
/// ns.upperMap(capacity).costMap(cost) |
898 | 898 |
/// .supplyMap(sup).run(); |
899 | 899 |
/// \endcode |
900 | 900 |
/// |
901 | 901 |
/// \return <tt>(*this)</tt> |
902 | 902 |
NetworkSimplex& reset() { |
903 | 903 |
for (int i = 0; i != _node_num; ++i) { |
904 | 904 |
_supply[i] = 0; |
905 | 905 |
} |
906 | 906 |
for (int i = 0; i != _arc_num; ++i) { |
907 | 907 |
_lower[i] = 0; |
908 | 908 |
_upper[i] = INF; |
909 | 909 |
_cost[i] = 1; |
910 | 910 |
} |
911 | 911 |
_have_lower = false; |
912 | 912 |
_stype = GEQ; |
913 | 913 |
return *this; |
914 | 914 |
} |
915 | 915 |
|
916 | 916 |
/// @} |
917 | 917 |
|
918 | 918 |
/// \name Query Functions |
919 | 919 |
/// The results of the algorithm can be obtained using these |
920 | 920 |
/// functions.\n |
921 | 921 |
/// The \ref run() function must be called before using them. |
922 | 922 |
|
923 | 923 |
/// @{ |
924 | 924 |
|
0 comments (0 inline)